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