source/wham/src-M/elecont.f

   1       subroutine elecont(lprint,ncont,icont,ist,ien)
   2       implicit none
   3       include 'DIMENSIONS'
   4       include 'DIMENSIONS.ZSCOPT'
   5       include 'DIMENSIONS.COMPAR'
   6       include 'COMMON.IOUNITS'
   7       include 'COMMON.CHAIN'
   8       include 'COMMON.INTERACT'
   9       include 'COMMON.FFIELD'
  10       include 'COMMON.NAMES'
  11       include 'COMMON.LOCAL'
  12       logical lprint
  13       integer i,j,k,ist,ien,iteli,itelj,ind,i1,i2,it1,it2,ic1,ic2
  14       double precision rri,xi,yi,zi,dxi,dyi,dzi,xmedi,ymedi,zmedi,
  15      &  xj,yj,zj,dxj,dyj,dzj,aaa,bbb,ael6i,ael3i,rrmij,rmij,r3ij,r6ij,
  16      &  vrmij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,evdwij,el1,el2,
  17      &  eesij,ees,evdw,ene, rij,zj_temp,xj_temp,yj_temp,
  18      & sscale,sscagrad,dist_temp,xj_safe,yj_safe,zj_safe,dist_init
  19       double precision elpp6c(2,2),elpp3c(2,2),ael6c(2,2),ael3c(2,2),
  20      &  appc(2,2),bppc(2,2)
  21       double precision elcutoff,elecutoff_14
  22       integer ncont,icont(2,maxcont),xshift,yshift,zshift,isubchap
  23       double precision econt(maxcont)
  24 *
  25 * Load the constants of peptide bond - peptide bond interactions.
  26 * Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
  27 * proline) - determined by averaging ECEPP energy.
  28 *
  29 * as of 7/06/91.
  30 *
  31 c      data epp    / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
  32 c      data rpp    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
  33       data elpp6c  /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
  34       data elpp3c  / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
  35       data elcutoff /-0.3d0/,elecutoff_14 /-0.5d0/
  36       ees=0.0d0
  37       evdw=0.0d0
  38       if (lprint) write (iout,'(a)')
  39      &  "Constants of electrostatic interaction energy expression."
  40       do i=1,2
  41         do j=1,2
  42         rri=rpp(i,j)**6
  43         appc(i,j)=epp(i,j)*rri*rri
  44         bppc(i,j)=-2.0*epp(i,j)*rri
  45         ael6c(i,j)=elpp6c(i,j)*4.2**6
  46         ael3c(i,j)=elpp3c(i,j)*4.2**3
  47         if (lprint)
  48      &  write (iout,'(2i2,4e15.4)') i,j,appc(i,j),bppc(i,j),ael6c(i,j),
  49      &                               ael3c(i,j)
  50         enddo
  51       enddo
  52       ncont=0
  53       do 1 i=ist,ien-2
  54         xi=c(1,i)
  55         yi=c(2,i)
  56         zi=c(3,i)
  57         dxi=c(1,i+1)-c(1,i)
  58         dyi=c(2,i+1)-c(2,i)
  59         dzi=c(3,i+1)-c(3,i)
  60         xmedi=xi+0.5*dxi
  61         ymedi=yi+0.5*dyi
  62         zmedi=zi+0.5*dzi
  63           xmedi=mod(xmedi,boxxsize)
  64           if (xmedi.lt.0) xmedi=xmedi+boxxsize
  65           ymedi=mod(ymedi,boxysize)
  66           if (ymedi.lt.0) ymedi=ymedi+boxysize
  67           zmedi=mod(zmedi,boxzsize)
  68           if (zmedi.lt.0) zmedi=zmedi+boxzsize
  69         do 4 j=i+2,ien-1
  70           ind=ind+1
  71           iteli=itel(i)
  72           itelj=itel(j)
  73           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
  74           if (iteli.eq.2 .and. itelj.eq.2
  75      &      .or.iteli.eq.0 .or.itelj.eq.0) goto 4
  76           aaa=appc(iteli,itelj)
  77           bbb=bppc(iteli,itelj)
  78           ael6i=ael6c(iteli,itelj)
  79           ael3i=ael3c(iteli,itelj)
  80           dxj=c(1,j+1)-c(1,j)
  81           dyj=c(2,j+1)-c(2,j)
  82           dzj=c(3,j+1)-c(3,j)
  83           xj=c(1,j)+0.5*dxj
  84           yj=c(2,j)+0.5*dyj
  85           zj=c(3,j)+0.5*dzj
  86           xj=mod(xj,boxxsize)
  87           if (xj.lt.0) xj=xj+boxxsize
  88           yj=mod(yj,boxysize)
  89           if (yj.lt.0) yj=yj+boxysize
  90           zj=mod(zj,boxzsize)
  91           if (zj.lt.0) zj=zj+boxzsize
  92       dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
  93       xj_safe=xj
  94       yj_safe=yj
  95       zj_safe=zj
  96       isubchap=0
  97       do xshift=-1,1
  98       do yshift=-1,1
  99       do zshift=-1,1
 100           xj=xj_safe+xshift*boxxsize
 101           yj=yj_safe+yshift*boxysize
 102           zj=zj_safe+zshift*boxzsize
 103           dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
 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=ael6i*r6ij
 138           fac4=ael3i*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     4   continue
 153     1 continue
 154       if (lprint) then
 155         write (iout,*) 'Total average electrostatic energy: ',ees
 156         write (iout,*) 'VDW energy between peptide-group centers: ',evdw
 157         write (iout,*)
 158         write (iout,*) 'Electrostatic contacts before pruning: '
 159         do i=1,ncont
 160           i1=icont(1,i)
 161           i2=icont(2,i)
 162           it1=itype(i1)
 163           it2=itype(i2)
 164           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
 165      &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
 166         enddo
 167       endif
 168 c For given residues keep only the contacts with the greatest energy.
 169       i=0
 170       do while (i.lt.ncont)
 171         i=i+1
 172         ene=econt(i)
 173         ic1=icont(1,i)
 174         ic2=icont(2,i)
 175         j=i
 176         do while (j.lt.ncont)
 177           j=j+1
 178           if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or.
 179      &        ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
 180 c            write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
 181 c     &       " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
 182             if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
 183               if (ic1.eq.icont(1,j)) then
 184                 do k=1,ncont
 185                   if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)
 186      &               .and. iabs(icont(1,k)-ic1).le.2 .and.
 187      &               econt(k).lt.econt(j) ) goto 21
 188                 enddo
 189               else if (ic2.eq.icont(2,j) ) then
 190                 do k=1,ncont
 191                   if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)
 192      &               .and. iabs(icont(2,k)-ic2).le.2 .and.
 193      &               econt(k).lt.econt(j) ) goto 21
 194                 enddo
 195               endif
 196 c Remove ith contact
 197               do k=i+1,ncont
 198                 icont(1,k-1)=icont(1,k)
 199                 icont(2,k-1)=icont(2,k)
 200                 econt(k-1)=econt(k)
 201               enddo
 202               i=i-1
 203               ncont=ncont-1
 204 c              write (iout,*) "ncont",ncont
 205 c              do k=1,ncont
 206 c                write (iout,*) icont(1,k),icont(2,k)
 207 c              enddo
 208               goto 20
 209             else if (econt(j).gt.ene .and. ic2.ne.ic1+2)
 210      &      then
 211               if (ic1.eq.icont(1,j)) then
 212                 do k=1,ncont
 213                   if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2
 214      &               .and. iabs(icont(1,k)-icont(1,j)).le.2 .and.
 215      &               econt(k).lt.econt(i) ) goto 21
 216                 enddo
 217               else if (ic2.eq.icont(2,j) ) then
 218                 do k=1,ncont
 219                   if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1
 220      &               .and. iabs(icont(2,k)-icont(2,j)).le.2 .and.
 221      &               econt(k).lt.econt(i) ) goto 21
 222                 enddo
 223               endif
 224 c Remove jth contact
 225               do k=j+1,ncont
 226                 icont(1,k-1)=icont(1,k)
 227                 icont(2,k-1)=icont(2,k)
 228                 econt(k-1)=econt(k)
 229               enddo
 230               ncont=ncont-1
 231 c              write (iout,*) "ncont",ncont
 232 c              do k=1,ncont
 233 c                write (iout,*) icont(1,k),icont(2,k)
 234 c              enddo
 235               j=j-1
 236             endif
 237           endif
 238    21     continue
 239         enddo
 240    20   continue
 241       enddo
 242       if (lprint) then
 243         write (iout,*)
 244         write (iout,*) 'Electrostatic contacts after pruning: '
 245         do i=1,ncont
 246           i1=icont(1,i)
 247           i2=icont(2,i)
 248           it1=itype(i1)
 249           it2=itype(i2)
 250           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
 251      &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
 252         enddo
 253       endif
 254       return
 255       end