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