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

   1       subroutine elecont(lprint,ncont,icont,ist,ien,ipermmin)
   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 iperm,ipermmin,ii,jj
  14       integer i,j,k,ist,ien,iteli,itelj,ind,i1,i2,it1,it2,ic1,ic2
  15       double precision rri,xi,yi,zi,dxi,dyi,dzi,xmedi,ymedi,zmedi,
  16      &  xj,yj,zj,dxj,dyj,dzj,aaa,bbb,ael6i,ael3i,rrmij,rmij,r3ij,r6ij,
  17      &  vrmij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,evdwij,el1,el2,
  18      &  eesij,ees,evdw,ene, rij,zj_temp,xj_temp,yj_temp,
  19      & sscale,sscagrad,dist_temp,xj_safe,yj_safe,zj_safe,dist_init
  20       double precision elpp6c(2,2),elpp3c(2,2),ael6c(2,2),ael3c(2,2),
  21      &  appc(2,2),bppc(2,2),epp_(2,2),rpp_(2,2)
  22       double precision elcutoff,elecutoff_14
  23       integer ncont,icont(2,maxcont),xshift,yshift,zshift,isubchap
  24       double precision econt(maxcont)
  25       double precision boxshift
  26 *
  27 * Load the constants of peptide bond - peptide bond interactions.
  28 * Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
  29 * proline) - determined by averaging ECEPP energy.
  30 *
  31 * as of 7/06/91.
  32 *
  33 c      data epp_   / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
  34       data rpp_   / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
  35       data elpp6c  /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
  36       data elpp3c  / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
  37       data elcutoff /-0.3d0/,elecutoff_14 /-0.5d0/
  38       ees=0.0d0
  39       evdw=0.0d0
  40       if (lprint) write (iout,'(a)')
  41      &  "Constants of electrostatic interaction energy expression."
  42       do i=1,2
  43         do j=1,2
  44         rri=rpp_(i,j)**6
  45         appc(i,j)=epp(i,j)*rri*rri
  46         bppc(i,j)=-2.0*epp(i,j)*rri
  47         ael6c(i,j)=elpp6c(i,j)*4.2**6
  48         ael3c(i,j)=elpp3c(i,j)*4.2**3
  49         if (lprint)
  50      &  write (iout,'(2i2,4e15.4)') i,j,appc(i,j),bppc(i,j),ael6c(i,j),
  51      &                               ael3c(i,j)
  52         enddo
  53       enddo
  54       ncont=0
  55       do 1 i=ist,ien-2
  56         ii=iperm(i,ipermmin)
  57         xi=c(1,ii)
  58         yi=c(2,ii)
  59         zi=c(3,ii)
  60         dxi=c(1,ii+1)-c(1,ii)
  61         dyi=c(2,ii+1)-c(2,ii)
  62         dzi=c(3,ii+1)-c(3,ii)
  63         xmedi=xi+0.5*dxi
  64         ymedi=yi+0.5*dyi
  65         zmedi=zi+0.5*dzi
  66         call to_box(xmedi,ymedi,zmedi)
  67 c        write (iout,*) "i",xmedi,ymedi,zmedi
  68         do 4 j=i+2,ien-1
  69           jj=iperm(j,ipermmin)
  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,jj+1)-c(1,jj)
  81           dyj=c(2,jj+1)-c(2,jj)
  82           dzj=c(3,jj+1)-c(3,jj)
  83           xj=c(1,jj)+0.5*dxj
  84           yj=c(2,jj)+0.5*dyj
  85           zj=c(3,jj)+0.5*dzj
  86 c          write (iout,*) "j",xj,yj,zj
  87           call to_box(xj,yj,zj)
  88           xj=boxshift(xj-xmedi,boxxsize)
  89           yj=boxshift(yj-ymedi,boxysize)
  90           zj=boxshift(zj-zmedi,boxzsize)
  91 c          write (iout,*) "j",xj,yj,zj
  92           rij=xj*xj+yj*yj+zj*zj
  93           rrmij=1.0/(xj*xj+yj*yj+zj*zj)
  94           rmij=sqrt(rrmij)
  95           r3ij=rrmij*rmij
  96           r6ij=r3ij*r3ij
  97           vrmij=vblinv*rmij
  98           cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
  99           cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
 100           cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
 101           fac=cosa-3.0*cosb*cosg
 102           ev1=aaa*r6ij*r6ij
 103           ev2=bbb*r6ij
 104           fac3=ael6i*r6ij
 105           fac4=ael3i*r3ij
 106           evdwij=ev1+ev2
 107           el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
 108           el2=fac4*fac
 109           eesij=el1+el2
 110           if (j.gt.i+2 .and. eesij.le.elcutoff .or.
 111      &        j.eq.i+2 .and. eesij.le.elecutoff_14) then
 112              ncont=ncont+1
 113              icont(1,ncont)=i
 114              icont(2,ncont)=j
 115              econt(ncont)=eesij
 116           endif
 117           ees=ees+eesij
 118           evdw=evdw+evdwij*sss
 119 c          write (iout,*) "i",i," j",j," rij",dsqrt(rij)," eesij",eesij
 120     4   continue
 121     1 continue
 122       if (lprint) then
 123         write (iout,*) 'Total average electrostatic energy: ',ees
 124         write (iout,*) 'VDW energy between peptide-group centers: ',evdw
 125         write (iout,*)
 126         write (iout,*) 'Electrostatic contacts before pruning: '
 127         do i=1,ncont
 128           i1=icont(1,i)
 129           i2=icont(2,i)
 130           it1=itype(i1)
 131           it2=itype(i2)
 132           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
 133      &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
 134         enddo
 135       endif
 136 c For given residues keep only the contacts with the greatest energy.
 137       i=0
 138       do while (i.lt.ncont)
 139         i=i+1
 140         ene=econt(i)
 141         ic1=icont(1,i)
 142         ic2=icont(2,i)
 143         j=i
 144         do while (j.lt.ncont)
 145           j=j+1
 146           if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or.
 147      &        ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
 148 c            write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
 149 c     &       " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
 150             if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
 151               if (ic1.eq.icont(1,j)) then
 152                 do k=1,ncont
 153                   if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)
 154      &               .and. iabs(icont(1,k)-ic1).le.2 .and.
 155      &               econt(k).lt.econt(j) ) goto 21
 156                 enddo
 157               else if (ic2.eq.icont(2,j) ) then
 158                 do k=1,ncont
 159                   if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)
 160      &               .and. iabs(icont(2,k)-ic2).le.2 .and.
 161      &               econt(k).lt.econt(j) ) goto 21
 162                 enddo
 163               endif
 164 c Remove ith contact
 165               do k=i+1,ncont
 166                 icont(1,k-1)=icont(1,k)
 167                 icont(2,k-1)=icont(2,k)
 168                 econt(k-1)=econt(k)
 169               enddo
 170               i=i-1
 171               ncont=ncont-1
 172 c              write (iout,*) "ncont",ncont
 173 c              do k=1,ncont
 174 c                write (iout,*) icont(1,k),icont(2,k)
 175 c              enddo
 176               goto 20
 177             else if (econt(j).gt.ene .and. ic2.ne.ic1+2)
 178      &      then
 179               if (ic1.eq.icont(1,j)) then
 180                 do k=1,ncont
 181                   if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2
 182      &               .and. iabs(icont(1,k)-icont(1,j)).le.2 .and.
 183      &               econt(k).lt.econt(i) ) goto 21
 184                 enddo
 185               else if (ic2.eq.icont(2,j) ) then
 186                 do k=1,ncont
 187                   if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1
 188      &               .and. iabs(icont(2,k)-icont(2,j)).le.2 .and.
 189      &               econt(k).lt.econt(i) ) goto 21
 190                 enddo
 191               endif
 192 c Remove jth contact
 193               do k=j+1,ncont
 194                 icont(1,k-1)=icont(1,k)
 195                 icont(2,k-1)=icont(2,k)
 196                 econt(k-1)=econt(k)
 197               enddo
 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               j=j-1
 204             endif
 205           endif
 206    21     continue
 207         enddo
 208    20   continue
 209       enddo
 210       if (lprint) then
 211         write (iout,*)
 212         write (iout,*) 'Electrostatic contacts after pruning: '
 213         do i=1,ncont
 214           i1=icont(1,i)
 215           i2=icont(2,i)
 216           it1=itype(i1)
 217           it2=itype(i2)
 218           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
 219      &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
 220         enddo
 221       endif
 222       return
 223       end