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