source/unres/src_MD-M/eelec.F

   1       subroutine eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
   2 C
   3 C This subroutine calculates the average interaction energy and its gradient
   4 C in the virtual-bond vectors between non-adjacent peptide groups, based on
   5 C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
   6 C The potential depends both on the distance of peptide-group centers and on
   7 C the orientation of the CA-CA virtual bonds.
   8 C
   9       implicit real*8 (a-h,o-z)
  10       include 'DIMENSIONS'
  11       include 'COMMON.CONTROL'
  12       include 'COMMON.IOUNITS'
  13       include 'COMMON.GEO'
  14       include 'COMMON.VAR'
  15       include 'COMMON.LOCAL'
  16       include 'COMMON.CHAIN'
  17       include 'COMMON.DERIV'
  18       include 'COMMON.INTERACT'
  19       include 'COMMON.CONTACTS'
  20       include 'COMMON.TORSION'
  21       include 'COMMON.VECTORS'
  22       include 'COMMON.FFIELD'
  23       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
  24      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
  25       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
  26      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
  27       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2
  28 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
  29 #ifdef MOMENT
  30       double precision scal_el /1.0d0/
  31 #else
  32       double precision scal_el /0.5d0/
  33 #endif
  34 C 12/13/98
  35 C 13-go grudnia roku pamietnego...
  36       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
  37      &                   0.0d0,1.0d0,0.0d0,
  38      &                   0.0d0,0.0d0,1.0d0/
  39 cd      write(iout,*) 'In EELEC'
  40 cd      do i=1,nloctyp
  41 cd        write(iout,*) 'Type',i
  42 cd        write(iout,*) 'B1',B1(:,i)
  43 cd        write(iout,*) 'B2',B2(:,i)
  44 cd        write(iout,*) 'CC',CC(:,:,i)
  45 cd        write(iout,*) 'DD',DD(:,:,i)
  46 cd        write(iout,*) 'EE',EE(:,:,i)
  47 cd      enddo
  48 cd      call check_vecgrad
  49 cd      stop
  50       if (icheckgrad.eq.1) then
  51         do i=1,nres-1
  52           fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
  53           do k=1,3
  54             dc_norm(k,i)=dc(k,i)*fac
  55           enddo
  56 c          write (iout,*) 'i',i,' fac',fac
  57         enddo
  58       endif
  59       if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
  60      &    .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
  61      &    wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
  62 c        call vec_and_deriv
  63         call set_matrices
  64       endif
  65 cd      do i=1,nres-1
  66 cd        write (iout,*) 'i=',i
  67 cd        do k=1,3
  68 cd        write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
  69 cd        enddo
  70 cd        do k=1,3
  71 cd          write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
  72 cd     &     uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
  73 cd        enddo
  74 cd      enddo
  75       num_conti_hb=0
  76       ees=0.0D0
  77       evdw1=0.0D0
  78       eel_loc=0.0d0
  79       eello_turn3=0.0d0
  80       eello_turn4=0.0d0
  81       ind=0
  82       do i=1,nres
  83         num_cont_hb(i)=0
  84       enddo
  85 cd      print '(a)','Enter EELEC'
  86 cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
  87       do i=1,nres
  88         gel_loc_loc(i)=0.0d0
  89         gcorr_loc(i)=0.0d0
  90       enddo
  91       do i=iatel_s,iatel_e
  92         dxi=dc(1,i)
  93         dyi=dc(2,i)
  94         dzi=dc(3,i)
  95         dx_normi=dc_norm(1,i)
  96         dy_normi=dc_norm(2,i)
  97         dz_normi=dc_norm(3,i)
  98         xmedi=c(1,i)+0.5d0*dxi
  99         ymedi=c(2,i)+0.5d0*dyi
 100         zmedi=c(3,i)+0.5d0*dzi
 101         num_conti=0
 102 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
 103         do j=ielstart(i),ielend(i)
 104           ind=ind+1
 105           iteli=itel(i)
 106           itelj=itel(j)
 107           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
 108           aaa=app(iteli,itelj)
 109           bbb=bpp(iteli,itelj)
 110           ael6i=ael6(iteli,itelj)
 111           ael3i=ael3(iteli,itelj)
 112 C Diagnostics only!!!
 113 c         aaa=0.0D0
 114 c         bbb=0.0D0
 115 c         ael6i=0.0D0
 116 c         ael3i=0.0D0
 117 C End diagnostics
 118           dxj=dc(1,j)
 119           dyj=dc(2,j)
 120           dzj=dc(3,j)
 121           dx_normj=dc_norm(1,j)
 122           dy_normj=dc_norm(2,j)
 123           dz_normj=dc_norm(3,j)
 124           xj=c(1,j)+0.5D0*dxj-xmedi
 125           yj=c(2,j)+0.5D0*dyj-ymedi
 126           zj=c(3,j)+0.5D0*dzj-zmedi
 127           rij=xj*xj+yj*yj+zj*zj
 128           rrmij=1.0D0/rij
 129           rij=dsqrt(rij)
 130           rmij=1.0D0/rij
 131           r3ij=rrmij*rmij
 132           r6ij=r3ij*r3ij
 133           cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
 134           cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
 135           cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
 136           fac=cosa-3.0D0*cosb*cosg
 137           ev1=aaa*r6ij*r6ij
 138 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
 139           if (j.eq.i+2) ev1=scal_el*ev1
 140           ev2=bbb*r6ij
 141           fac3=ael6i*r6ij
 142           fac4=ael3i*r3ij
 143           evdwij=ev1+ev2
 144           el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
 145           el2=fac4*fac
 146           eesij=el1+el2
 147 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
 148           ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
 149           ees=ees+eesij
 150           evdw1=evdw1+evdwij
 151 cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
 152 cd     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
 153 cd     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
 154 cd     &      xmedi,ymedi,zmedi,xj,yj,zj
 155
 156           if (energy_dec) then
 157               write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
 158               write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
 159           endif
 160
 161 C
 162 C Calculate contributions to the Cartesian gradient.
 163 C
 164 #ifdef SPLITELE
 165           facvdw=-6*rrmij*(ev1+evdwij)
 166           facel=-3*rrmij*(el1+eesij)
 167           fac1=fac
 168           erij(1)=xj*rmij
 169           erij(2)=yj*rmij
 170           erij(3)=zj*rmij
 171 *
 172 * Radial derivatives. First process both termini of the fragment (i,j)
 173 *
 174           ggg(1)=facel*xj
 175           ggg(2)=facel*yj
 176           ggg(3)=facel*zj
 177           do k=1,3
 178             ghalf=0.5D0*ggg(k)
 179             gelc(k,i)=gelc(k,i)+ghalf
 180             gelc(k,j)=gelc(k,j)+ghalf
 181           enddo
 182 *
 183 * Loop over residues i+1 thru j-1.
 184 *
 185 caug8          do k=i+1,j-1
 186 caug8            do l=1,3
 187 caug8              gelc(l,k)=gelc(l,k)+ggg(l)
 188 caug8            enddo
 189 caug8          enddo
 190           ggg(1)=facvdw*xj
 191           ggg(2)=facvdw*yj
 192           ggg(3)=facvdw*zj
 193           do k=1,3
 194             ghalf=0.5D0*ggg(k)
 195             gvdwpp(k,i)=gvdwpp(k,i)+ghalf
 196             gvdwpp(k,j)=gvdwpp(k,j)+ghalf
 197           enddo
 198 *
 199 * Loop over residues i+1 thru j-1.
 200 *
 201 caug8          do k=i+1,j-1
 202 caug8            do l=1,3
 203 caug8              gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
 204 caug8            enddo
 205 caug8          enddo
 206 #else
 207           facvdw=ev1+evdwij
 208           facel=el1+eesij
 209           fac1=fac
 210           fac=-3*rrmij*(facvdw+facvdw+facel)
 211           erij(1)=xj*rmij
 212           erij(2)=yj*rmij
 213           erij(3)=zj*rmij
 214 *
 215 * Radial derivatives. First process both termini of the fragment (i,j)
 216 *
 217           ggg(1)=fac*xj
 218           ggg(2)=fac*yj
 219           ggg(3)=fac*zj
 220           do k=1,3
 221             ghalf=0.5D0*ggg(k)
 222             gelc(k,i)=gelc(k,i)+ghalf
 223             gelc(k,j)=gelc(k,j)+ghalf
 224           enddo
 225 *
 226 * Loop over residues i+1 thru j-1.
 227 *
 228 caug8          do k=i+1,j-1
 229 caug8            do l=1,3
 230 caug8              gelc(l,k)=gelc(l,k)+ggg(l)
 231 caug8            enddo
 232 caug8          enddo
 233 #endif
 234 *
 235 * Angular part
 236 *
 237           ecosa=2.0D0*fac3*fac1+fac4
 238           fac4=-3.0D0*fac4
 239           fac3=-6.0D0*fac3
 240           ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
 241           ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
 242           do k=1,3
 243             dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
 244             dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
 245           enddo
 246 cd        print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
 247 cd   &          (dcosg(k),k=1,3)
 248           do k=1,3
 249             ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
 250           enddo
 251           do k=1,3
 252             ghalf=0.5D0*ggg(k)
 253             gelc(k,i)=gelc(k,i)+ghalf
 254      &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
 255      &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
 256             gelc(k,j)=gelc(k,j)+ghalf
 257      &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
 258      &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
 259           enddo
 260 caug8          do k=i+1,j-1
 261 caug8            do l=1,3
 262 caug8              gelc(l,k)=gelc(l,k)+ggg(l)
 263 caug8            enddo
 264 caug8          enddo
 265
 266         enddo ! j
 267         num_cont_hb(i)=num_conti
 268       enddo   ! i
 269 c      write (iout,*) "Number of loop steps in EELEC:",ind
 270 cd      do i=1,nres
 271 cd        write (iout,'(i3,3f10.5,5x,3f10.5)')
 272 cd     &     i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
 273 cd      enddo
 274 c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
 275 ccc      eel_loc=eel_loc+eello_turn3
 276       return
 277       end
 278