source/unres/src-HCD-5D/elecont.f

   1       subroutine elecont(lprint,ncont,icont)
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'COMMON.IOUNITS'
   5       include 'COMMON.CHAIN'
   6       include 'COMMON.INTERACT'
   7       include 'COMMON.LOCAL'
   8       include 'COMMON.FFIELD'
   9       include 'COMMON.NAMES'
  10       logical lprint
  11       double precision elpp_6(2,2),elpp_3(2,2),ael6_(2,2),ael3_(2,2)
  12       double precision app_(2,2),bpp_(2,2),rpp_(2,2)
  13       integer ncont,icont(2,maxcont)
  14       double precision econt(maxcont)
  15 *
  16 * Load the constants of peptide bond - peptide bond interactions.
  17 * Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
  18 * proline) - determined by averaging ECEPP energy.
  19 *
  20 * as of 7/06/91.
  21 *
  22 c      data epp    / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
  23       data rpp_    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
  24       data elpp_6  /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
  25       data elpp_3  / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
  26       data elcutoff /-0.3d0/,elecutoff_14 /-0.5d0/
  27       if (lprint) write (iout,'(a)')
  28      &  "Constants of electrostatic interaction energy expression."
  29       do i=1,2
  30         do j=1,2
  31         rri=rpp_(i,j)**6
  32         app_(i,j)=epp(i,j)*rri*rri
  33         bpp_(i,j)=-2.0*epp(i,j)*rri
  34         ael6_(i,j)=elpp_6(i,j)*4.2**6
  35         ael3_(i,j)=elpp_3(i,j)*4.2**3
  36         if (lprint)
  37      &  write (iout,'(2i2,4e15.4)') i,j,app_(i,j),bpp_(i,j),ael6_(i,j),
  38      &                               ael3_(i,j)
  39         enddo
  40       enddo
  41       ncont=0
  42       ees=0.0
  43       evdw=0.0
  44       do 1 i=nnt,nct-2
  45         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) goto 1
  46         xi=c(1,i)
  47         yi=c(2,i)
  48         zi=c(3,i)
  49         dxi=c(1,i+1)-c(1,i)
  50         dyi=c(2,i+1)-c(2,i)
  51         dzi=c(3,i+1)-c(3,i)
  52         xmedi=xi+0.5*dxi
  53         ymedi=yi+0.5*dyi
  54         zmedi=zi+0.5*dzi
  55           xmedi=mod(xmedi,boxxsize)
  56           if (xmedi.lt.0) xmedi=xmedi+boxxsize
  57           ymedi=mod(ymedi,boxysize)
  58           if (ymedi.lt.0) ymedi=ymedi+boxysize
  59           zmedi=mod(zmedi,boxzsize)
  60           if (zmedi.lt.0) zmedi=zmedi+boxzsize
  61         do 4 j=i+2,nct-1
  62           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) goto 4
  63           ind=ind+1
  64           iteli=itel(i)
  65           itelj=itel(j)
  66           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
  67           if (iteli.eq.2 .and. itelj.eq.2) goto 4
  68           aaa=app_(iteli,itelj)
  69           bbb=bpp_(iteli,itelj)
  70           ael6_i=ael6_(iteli,itelj)
  71           ael3_i=ael3_(iteli,itelj)
  72           dxj=c(1,j+1)-c(1,j)
  73           dyj=c(2,j+1)-c(2,j)
  74           dzj=c(3,j+1)-c(3,j)
  75           xj=c(1,j)+0.5*dxj
  76           yj=c(2,j)+0.5*dyj
  77           zj=c(3,j)+0.5*dzj
  78           xj=mod(xj,boxxsize)
  79           if (xj.lt.0) xj=xj+boxxsize
  80           yj=mod(yj,boxysize)
  81           if (yj.lt.0) yj=yj+boxysize
  82           zj=mod(zj,boxzsize)
  83           if (zj.lt.0) zj=zj+boxzsize
  84       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
  85       write (iout,*) "i",i,xi,yi,zi," j",j,xj,yj,xj,"dist",
  86      &   dsqrt(dist_init)
  87       xj_safe=xj
  88       yj_safe=yj
  89       zj_safe=zj
  90       isubchap=0
  91       do xshift=-1,1
  92       do yshift=-1,1
  93       do zshift=-1,1
  94           xj=xj_safe+xshift*boxxsize
  95           yj=yj_safe+yshift*boxysize
  96           zj=zj_safe+zshift*boxzsize
  97           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
  98           if(dist_temp.lt.dist_init) then
  99             dist_init=dist_temp
 100             xj_temp=xj
 101             yj_temp=yj
 102             zj_temp=zj
 103             isubchap=1
 104           endif
 105        enddo
 106        enddo
 107        enddo
 108        if (isubchap.eq.1) then
 109           xj=xj_temp-xmedi
 110           yj=yj_temp-ymedi
 111           zj=zj_temp-zmedi
 112        else
 113           xj=xj_safe-xmedi
 114           yj=yj_safe-ymedi
 115           zj=zj_safe-zmedi
 116        endif
 117           rij=xj*xj+yj*yj+zj*zj
 118             sss=sscale(sqrt(rij))
 119             sssgrad=sscagrad(sqrt(rij))
 120           rrmij=1.0/(xj*xj+yj*yj+zj*zj)
 121           rmij=sqrt(rrmij)
 122           r3ij=rrmij*rmij
 123           r6ij=r3ij*r3ij
 124           vrmij=vblinv*rmij
 125           cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
 126           cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
 127           cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
 128           fac=cosa-3.0*cosb*cosg
 129           ev1=aaa*r6ij*r6ij
 130           ev2=bbb*r6ij
 131           fac3=ael6_i*r6ij
 132           fac4=ael3_i*r3ij
 133           evdwij=ev1+ev2
 134           el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
 135           el2=fac4*fac
 136           eesij=el1+el2
 137           if (j.gt.i+2 .and. eesij.le.elcutoff .or.
 138      &        j.eq.i+2 .and. eesij.le.elecutoff_14) then
 139              ncont=ncont+1
 140              icont(1,ncont)=i
 141              icont(2,ncont)=j
 142              econt(ncont)=eesij
 143           endif
 144           ees=ees+eesij
 145           evdw=evdw+evdwij*sss
 146           write (iout,*) "i"," j",j," rij",dsqrt(rij)," eesij",eesij
 147     4   continue
 148     1 continue
 149       if (lprint) then
 150         write (iout,*) 'Total average electrostatic energy: ',ees
 151         write (iout,*) 'VDW energy between peptide-group centers: ',evdw
 152         write (iout,*)
 153         write (iout,*) 'Electrostatic contacts before pruning: '
 154         do i=1,ncont
 155           i1=icont(1,i)
 156           i2=icont(2,i)
 157           it1=itype(i1)
 158           it2=itype(i2)
 159           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
 160      &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
 161         enddo
 162       endif
 163 c For given residues keep only the contacts with the greatest energy.
 164       i=0
 165       do while (i.lt.ncont)
 166         i=i+1
 167         ene=econt(i)
 168         ic1=icont(1,i)
 169         ic2=icont(2,i)
 170         j=i
 171         do while (j.lt.ncont)
 172           j=j+1
 173           if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or.
 174      &        ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
 175 c            write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
 176 c     &       " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
 177             if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
 178               if (ic1.eq.icont(1,j)) then
 179                 do k=1,ncont
 180                   if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)
 181      &               .and. iabs(icont(1,k)-ic1).le.2 .and.
 182      &               econt(k).lt.econt(j) ) goto 21
 183                 enddo
 184               else if (ic2.eq.icont(2,j) ) then
 185                 do k=1,ncont
 186                   if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)
 187      &               .and. iabs(icont(2,k)-ic2).le.2 .and.
 188      &               econt(k).lt.econt(j) ) goto 21
 189                 enddo
 190               endif
 191 c Remove ith contact
 192               do k=i+1,ncont
 193                 icont(1,k-1)=icont(1,k)
 194                 icont(2,k-1)=icont(2,k)
 195                 econt(k-1)=econt(k)
 196               enddo
 197               i=i-1
 198               ncont=ncont-1
 199 c              write (iout,*) "ncont",ncont
 200 c              do k=1,ncont
 201 c                write (iout,*) icont(1,k),icont(2,k)
 202 c              enddo
 203               goto 20
 204             else if (econt(j).gt.ene .and. ic2.ne.ic1+2)
 205      &      then
 206               if (ic1.eq.icont(1,j)) then
 207                 do k=1,ncont
 208                   if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2
 209      &               .and. iabs(icont(1,k)-icont(1,j)).le.2 .and.
 210      &               econt(k).lt.econt(i) ) goto 21
 211                 enddo
 212               else if (ic2.eq.icont(2,j) ) then
 213                 do k=1,ncont
 214                   if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1
 215      &               .and. iabs(icont(2,k)-icont(2,j)).le.2 .and.
 216      &               econt(k).lt.econt(i) ) goto 21
 217                 enddo
 218               endif
 219 c Remove jth contact
 220               do k=j+1,ncont
 221                 icont(1,k-1)=icont(1,k)
 222                 icont(2,k-1)=icont(2,k)
 223                 econt(k-1)=econt(k)
 224               enddo
 225               ncont=ncont-1
 226 c              write (iout,*) "ncont",ncont
 227 c              do k=1,ncont
 228 c                write (iout,*) icont(1,k),icont(2,k)
 229 c              enddo
 230               j=j-1
 231             endif
 232           endif
 233    21     continue
 234         enddo
 235    20   continue
 236       enddo
 237       if (lprint) then
 238         write (iout,*)
 239         write (iout,*) 'Electrostatic contacts after pruning: '
 240         do i=1,ncont
 241           i1=icont(1,i)
 242           i2=icont(2,i)
 243           it1=itype(i1)
 244           it2=itype(i2)
 245           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
 246      &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
 247         enddo
 248       endif
 249       return
 250       end
 251 c--------------------------------------------
 252       subroutine secondary2(lprint)
 253       implicit real*8 (a-h,o-z)
 254       include 'DIMENSIONS'
 255       include 'COMMON.CHAIN'
 256       include 'COMMON.IOUNITS'
 257       include 'COMMON.FRAG'
 258       include 'COMMON.VAR'
 259       include 'COMMON.GEO'
 260       include 'COMMON.CONTROL'
 261       integer ncont,icont(2,maxcont),isec(maxres,4),nsec(maxres)
 262       logical lprint,not_done,freeres
 263       double precision p1,p2
 264       external freeres
 265
 266 cc????      if(.not.dccart) call chainbuild
 267 cd      call write_pdb(99,'sec structure',0d0)
 268       ncont=0
 269       nbfrag=0
 270       nhfrag=0
 271       do i=1,nres
 272         isec(i,1)=0
 273         isec(i,2)=0
 274         nsec(i)=0
 275       enddo
 276
 277       call elecont(lprint,ncont,icont)
 278
 279 c finding parallel beta
 280 cd      write (iout,*) '------- looking for parallel beta -----------'
 281       nbeta=0
 282       nstrand=0
 283       do i=1,ncont
 284         i1=icont(1,i)
 285         j1=icont(2,i)
 286         if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
 287           ii1=i1
 288           jj1=j1
 289 cd          write (iout,*) i1,j1
 290           not_done=.true.
 291           do while (not_done)
 292            i1=i1+1
 293            j1=j1+1
 294             do j=1,ncont
 295               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and.
 296      &             freeres(i1,j1,nsec,isec)) goto 5
 297             enddo
 298             not_done=.false.
 299   5         continue
 300 cd            write (iout,*) i1,j1,not_done
 301           enddo
 302           j1=j1-1
 303           i1=i1-1
 304           if (i1-ii1.gt.1) then
 305             ii1=max0(ii1-1,1)
 306             jj1=max0(jj1-1,1)
 307             nbeta=nbeta+1
 308             if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',
 309      &               nbeta,ii1,i1,jj1,j1
 310
 311             nbfrag=nbfrag+1
 312             bfrag(1,nbfrag)=ii1+1
 313             bfrag(2,nbfrag)=i1+1
 314             bfrag(3,nbfrag)=jj1+1
 315             bfrag(4,nbfrag)=min0(j1+1,nres)
 316
 317             do ij=ii1,i1
 318              nsec(ij)=nsec(ij)+1
 319              isec(ij,nsec(ij))=nbeta
 320             enddo
 321             do ij=jj1,j1
 322              nsec(ij)=nsec(ij)+1
 323              isec(ij,nsec(ij))=nbeta
 324             enddo
 325
 326            if(lprint) then
 327             nstrand=nstrand+1
 328             if (nbeta.le.9) then
 329               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
 330      &          "DefPropRes 'strand",nstrand,
 331      &          "' 'num = ",ii1-1,"..",i1-1,"'"
 332             else
 333               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
 334      &          "DefPropRes 'strand",nstrand,
 335      &          "' 'num = ",ii1-1,"..",i1-1,"'"
 336             endif
 337             nstrand=nstrand+1
 338             if (nbeta.le.9) then
 339               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
 340      &          "DefPropRes 'strand",nstrand,
 341      &          "' 'num = ",jj1-1,"..",j1-1,"'"
 342             else
 343               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
 344      &          "DefPropRes 'strand",nstrand,
 345      &          "' 'num = ",jj1-1,"..",j1-1,"'"
 346             endif
 347               write(12,'(a8,4i4)')
 348      &          "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
 349            endif
 350           endif
 351         endif
 352       enddo
 353
 354 c finding alpha or 310 helix
 355
 356       nhelix=0
 357       do i=1,ncont
 358         i1=icont(1,i)
 359         j1=icont(2,i)
 360         p1=phi(i1+2)*rad2deg
 361         p2=0.0
 362         if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
 363
 364
 365         if (j1.eq.i1+3 .and.
 366      &       ((p1.ge.10.and.p1.le.80).or.i1.le.2).and.
 367      &       ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
 368 cd          if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
 369 co          if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
 370           ii1=i1
 371           jj1=j1
 372           if (nsec(ii1).eq.0) then
 373             not_done=.true.
 374           else
 375             not_done=.false.
 376           endif
 377           do while (not_done)
 378             i1=i1+1
 379             j1=j1+1
 380             do j=1,ncont
 381               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
 382             enddo
 383             not_done=.false.
 384   10        continue
 385             p1=phi(i1+2)*rad2deg
 386             p2=phi(j1+2)*rad2deg
 387             if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80)
 388      &                              not_done=.false.
 389 cd           write (iout,*) i1,j1,not_done,p1,p2
 390           enddo
 391           j1=j1+1
 392           if (j1-ii1.gt.5) then
 393             nhelix=nhelix+1
 394 cd            write (iout,*)'helix',nhelix,ii1,j1
 395
 396             nhfrag=nhfrag+1
 397             hfrag(1,nhfrag)=ii1
 398             hfrag(2,nhfrag)=j1
 399
 400             do ij=ii1,j1
 401              nsec(ij)=-1
 402             enddo
 403            if (lprint) then
 404             write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
 405             if (nhelix.le.9) then
 406               write(12,'(a17,i1,a9,i3,a2,i3,a1)')
 407      &          "DefPropRes 'helix",nhelix,
 408      &          "' 'num = ",ii1-1,"..",j1-2,"'"
 409             else
 410               write(12,'(a17,i2,a9,i3,a2,i3,a1)')
 411      &          "DefPropRes 'helix",nhelix,
 412      &          "' 'num = ",ii1-1,"..",j1-2,"'"
 413             endif
 414            endif
 415           endif
 416         endif
 417       enddo
 418
 419       if (nhelix.gt.0.and.lprint) then
 420         write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
 421         do i=2,nhelix
 422          if (nhelix.le.9) then
 423           write(12,'(a8,i1,$)') " | helix",i
 424          else
 425           write(12,'(a8,i2,$)') " | helix",i
 426          endif
 427         enddo
 428         write(12,'(a1)') "'"
 429       endif
 430
 431
 432 c finding antiparallel beta
 433 cd      write (iout,*) '--------- looking for antiparallel beta ---------'
 434
 435       do i=1,ncont
 436         i1=icont(1,i)
 437         j1=icont(2,i)
 438         if (freeres(i1,j1,nsec,isec)) then
 439           ii1=i1
 440           jj1=j1
 441 cd          write (iout,*) i1,j1
 442
 443           not_done=.true.
 444           do while (not_done)
 445            i1=i1+1
 446            j1=j1-1
 447             do j=1,ncont
 448               if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and.
 449      &             freeres(i1,j1,nsec,isec)) goto 6
 450             enddo
 451             not_done=.false.
 452   6         continue
 453 cd            write (iout,*) i1,j1,not_done
 454           enddo
 455           i1=i1-1
 456           j1=j1+1
 457           if (i1-ii1.gt.1) then
 458
 459             nbfrag=nbfrag+1
 460             bfrag(1,nbfrag)=ii1
 461             bfrag(2,nbfrag)=min0(i1+1,nres)
 462             bfrag(3,nbfrag)=min0(jj1+1,nres)
 463             bfrag(4,nbfrag)=j1
 464
 465             nbeta=nbeta+1
 466             iii1=max0(ii1-1,1)
 467             do ij=iii1,i1
 468              nsec(ij)=nsec(ij)+1
 469              if (nsec(ij).le.2) then
 470               isec(ij,nsec(ij))=nbeta
 471              endif
 472             enddo
 473             jjj1=max0(j1-1,1)
 474             do ij=jjj1,jj1
 475              nsec(ij)=nsec(ij)+1
 476              if (nsec(ij).le.2 .and. nsec(ij).gt.0) then
 477               isec(ij,nsec(ij))=nbeta
 478              endif
 479             enddo
 480
 481
 482            if (lprint) then
 483             write (iout,'(a,i3,4i4)')'antiparallel beta',
 484      &                   nbeta,ii1-1,i1,jj1,j1-1
 485             nstrand=nstrand+1
 486             if (nstrand.le.9) then
 487               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
 488      &          "DefPropRes 'strand",nstrand,
 489      &          "' 'num = ",ii1-2,"..",i1-1,"'"
 490             else
 491               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
 492      &          "DefPropRes 'strand",nstrand,
 493      &          "' 'num = ",ii1-2,"..",i1-1,"'"
 494             endif
 495             nstrand=nstrand+1
 496             if (nstrand.le.9) then
 497               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
 498      &          "DefPropRes 'strand",nstrand,
 499      &          "' 'num = ",j1-2,"..",jj1-1,"'"
 500             else
 501               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
 502      &          "DefPropRes 'strand",nstrand,
 503      &          "' 'num = ",j1-2,"..",jj1-1,"'"
 504             endif
 505               write(12,'(a8,4i4)')
 506      &          "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
 507            endif
 508           endif
 509         endif
 510       enddo
 511
 512       if (nstrand.gt.0.and.lprint) then
 513         write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
 514         do i=2,nstrand
 515          if (i.le.9) then
 516           write(12,'(a9,i1,$)') " | strand",i
 517          else
 518           write(12,'(a9,i2,$)') " | strand",i
 519          endif
 520         enddo
 521         write(12,'(a1)') "'"
 522       endif
 523
 524
 525
 526       if (lprint) then
 527        write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
 528        write(12,'(a20)') "XMacStand ribbon.mac"
 529
 530
 531        write(iout,*) 'UNRES seq:'
 532        do j=1,nbfrag
 533         write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
 534        enddo
 535
 536        do j=1,nhfrag
 537         write(iout,*) 'helix ',(hfrag(i,j),i=1,2)
 538        enddo
 539       endif
 540
 541       return
 542       end
 543 c-------------------------------------------------
 544       logical function freeres(i,j,nsec,isec)
 545       implicit real*8 (a-h,o-z)
 546       include 'DIMENSIONS'
 547       integer isec(maxres,4),nsec(maxres)
 548       freeres=.false.
 549
 550       if (nsec(i).lt.0.or.nsec(j).lt.0) return
 551       if (nsec(i).gt.1.or.nsec(j).gt.1) return
 552       do k=1,nsec(i)
 553         do l=1,nsec(j)
 554           if (isec(i,k).eq.isec(j,l)) return
 555         enddo
 556       enddo
 557       freeres=.true.
 558       return
 559       end
 560