source/unres/src-HCD-5D/eelecij_scale.F

   1       subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4 #ifdef MPI
   5       include "mpif.h"
   6 #endif
   7       include 'COMMON.CONTROL'
   8       include 'COMMON.IOUNITS'
   9       include 'COMMON.GEO'
  10       include 'COMMON.VAR'
  11       include 'COMMON.LOCAL'
  12       include 'COMMON.CHAIN'
  13       include 'COMMON.DERIV'
  14       include 'COMMON.INTERACT'
  15       include 'COMMON.CONTACTS'
  16       include 'COMMON.TORSION'
  17       include 'COMMON.VECTORS'
  18       include 'COMMON.FFIELD'
  19       include 'COMMON.TIME1'
  20       include 'COMMON.SHIELD'
  21       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
  22      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
  23       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
  24      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
  25       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
  26      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
  27      &    num_conti,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 c          time00=MPI_Wtime()
  40 cd      write (iout,*) "eelecij",i,j
  41 C      print *,"WCHODZE2"
  42           ind=ind+1
  43           iteli=itel(i)
  44           itelj=itel(j)
  45           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
  46           aaa=app(iteli,itelj)
  47           bbb=bpp(iteli,itelj)
  48           ael6i=ael6(iteli,itelj)
  49           ael3i=ael3(iteli,itelj)
  50           dxj=dc(1,j)
  51           dyj=dc(2,j)
  52           dzj=dc(3,j)
  53           dx_normj=dc_norm(1,j)
  54           dy_normj=dc_norm(2,j)
  55           dz_normj=dc_norm(3,j)
  56           xj=c(1,j)+0.5D0*dxj
  57           yj=c(2,j)+0.5D0*dyj
  58           zj=c(3,j)+0.5D0*dzj
  59          xj=mod(xj,boxxsize)
  60           if (xj.lt.0) xj=xj+boxxsize
  61           yj=mod(yj,boxysize)
  62           if (yj.lt.0) yj=yj+boxysize
  63           zj=mod(zj,boxzsize)
  64           if (zj.lt.0) zj=zj+boxzsize
  65       dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
  66       xj_safe=xj
  67       yj_safe=yj
  68       zj_safe=zj
  69       isubchap=0
  70       do xshift=-1,1
  71       do yshift=-1,1
  72       do zshift=-1,1
  73           xj=xj_safe+xshift*boxxsize
  74           yj=yj_safe+yshift*boxysize
  75           zj=zj_safe+zshift*boxzsize
  76           dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
  77           if(dist_temp.lt.dist_init) then
  78             dist_init=dist_temp
  79             xj_temp=xj
  80             yj_temp=yj
  81             zj_temp=zj
  82             isubchap=1
  83           endif
  84        enddo
  85        enddo
  86        enddo
  87        if (isubchap.eq.1) then
  88           xj=xj_temp-xmedi
  89           yj=yj_temp-ymedi
  90           zj=zj_temp-zmedi
  91        else
  92           xj=xj_safe-xmedi
  93           yj=yj_safe-ymedi
  94           zj=zj_safe-zmedi
  95        endif
  96
  97           rij=xj*xj+yj*yj+zj*zj
  98           rrmij=1.0D0/rij
  99           rij=dsqrt(rij)
 100           rmij=1.0D0/rij
 101 c For extracting the short-range part of Evdwpp
 102           sss=sscale(rij/rpp(iteli,itelj))
 103           sssgrad=sscagrad(rij/rpp(iteli,itelj))
 104           r3ij=rrmij*rmij
 105           r6ij=r3ij*r3ij
 106           cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
 107           cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
 108           cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
 109           fac=cosa-3.0D0*cosb*cosg
 110           ev1=aaa*r6ij*r6ij
 111 c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
 112           if (j.eq.i+2) ev1=scal_el*ev1
 113           ev2=bbb*r6ij
 114           fac3=ael6i*r6ij
 115           fac4=ael3i*r3ij
 116           evdwij=ev1+ev2
 117           if (shield_mode.eq.0) then
 118           fac_shield(i)=1.0
 119           fac_shield(j)=1.0
 120           endif
 121           el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
 122           el2=fac4*fac
 123           el1=el1*fac_shield(i)**2*fac_shield(j)**2
 124           el2=el2*fac_shield(i)**2*fac_shield(j)**2
 125           eesij=el1+el2
 126 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
 127           ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
 128           ees=ees+eesij
 129           evdw1=evdw1+evdwij*(1.0d0-sss)
 130 cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
 131 cd     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
 132 cd     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
 133 cd     &      xmedi,ymedi,zmedi,xj,yj,zj
 134
 135           if (energy_dec) then
 136               write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
 137               write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
 138           endif
 139
 140 C
 141 C Calculate contributions to the Cartesian gradient.
 142 C
 143 #ifdef SPLITELE
 144           facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
 145           facel=-3*rrmij*(el1+eesij)
 146           fac1=fac
 147           erij(1)=xj*rmij
 148           erij(2)=yj*rmij
 149           erij(3)=zj*rmij
 150 *
 151 * Radial derivatives. First process both termini of the fragment (i,j)
 152 *
 153           ggg(1)=facel*xj
 154           ggg(2)=facel*yj
 155           ggg(3)=facel*zj
 156           if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
 157      &  (shield_mode.gt.0)) then
 158 C          print *,i,j
 159           do ilist=1,ishield_list(i)
 160            iresshield=shield_list(ilist,i)
 161            do k=1,3
 162            rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
 163      &      *2.0
 164            gshieldx(k,iresshield)=gshieldx(k,iresshield)+
 165      &              rlocshield
 166      & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0
 167             gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
 168 C           gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
 169 C     & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
 170 C             if (iresshield.gt.i) then
 171 C               do ishi=i+1,iresshield-1
 172 C                gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
 173 C     & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
 174 C
 175 C              enddo
 176 C             else
 177 C               do ishi=iresshield,i
 178 C                gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
 179 C     & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
 180 C
 181 C               enddo
 182 C              endif
 183            enddo
 184           enddo
 185           do ilist=1,ishield_list(j)
 186            iresshield=shield_list(ilist,j)
 187            do k=1,3
 188            rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
 189      &     *2.0
 190            gshieldx(k,iresshield)=gshieldx(k,iresshield)+
 191      &              rlocshield
 192      & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
 193            gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
 194            enddo
 195           enddo
 196
 197           do k=1,3
 198             gshieldc(k,i)=gshieldc(k,i)+
 199      &              grad_shield(k,i)*eesij/fac_shield(i)*2.0
 200             gshieldc(k,j)=gshieldc(k,j)+
 201      &              grad_shield(k,j)*eesij/fac_shield(j)*2.0
 202             gshieldc(k,i-1)=gshieldc(k,i-1)+
 203      &              grad_shield(k,i)*eesij/fac_shield(i)*2.0
 204             gshieldc(k,j-1)=gshieldc(k,j-1)+
 205      &              grad_shield(k,j)*eesij/fac_shield(j)*2.0
 206
 207            enddo
 208            endif
 209
 210 c          do k=1,3
 211 c            ghalf=0.5D0*ggg(k)
 212 c            gelc(k,i)=gelc(k,i)+ghalf
 213 c            gelc(k,j)=gelc(k,j)+ghalf
 214 c          enddo
 215 c 9/28/08 AL Gradient compotents will be summed only at the end
 216           do k=1,3
 217             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
 218             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
 219           enddo
 220 *
 221 * Loop over residues i+1 thru j-1.
 222 *
 223 cgrad          do k=i+1,j-1
 224 cgrad            do l=1,3
 225 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
 226 cgrad            enddo
 227 cgrad          enddo
 228           ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
 229           ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
 230           ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
 231 c          do k=1,3
 232 c            ghalf=0.5D0*ggg(k)
 233 c            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
 234 c            gvdwpp(k,j)=gvdwpp(k,j)+ghalf
 235 c          enddo
 236 c 9/28/08 AL Gradient compotents will be summed only at the end
 237           do k=1,3
 238             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
 239             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
 240           enddo
 241 *
 242 * Loop over residues i+1 thru j-1.
 243 *
 244 cgrad          do k=i+1,j-1
 245 cgrad            do l=1,3
 246 cgrad              gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
 247 cgrad            enddo
 248 cgrad          enddo
 249 #else
 250           facvdw=ev1+evdwij*(1.0d0-sss)
 251           facel=el1+eesij
 252           fac1=fac
 253           fac=-3*rrmij*(facvdw+facvdw+facel)
 254           erij(1)=xj*rmij
 255           erij(2)=yj*rmij
 256           erij(3)=zj*rmij
 257 *
 258 * Radial derivatives. First process both termini of the fragment (i,j)
 259 *
 260           ggg(1)=fac*xj
 261           ggg(2)=fac*yj
 262           ggg(3)=fac*zj
 263 c          do k=1,3
 264 c            ghalf=0.5D0*ggg(k)
 265 c            gelc(k,i)=gelc(k,i)+ghalf
 266 c            gelc(k,j)=gelc(k,j)+ghalf
 267 c          enddo
 268 c 9/28/08 AL Gradient compotents will be summed only at the end
 269           do k=1,3
 270             gelc_long(k,j)=gelc(k,j)+ggg(k)
 271             gelc_long(k,i)=gelc(k,i)-ggg(k)
 272           enddo
 273 *
 274 * Loop over residues i+1 thru j-1.
 275 *
 276 cgrad          do k=i+1,j-1
 277 cgrad            do l=1,3
 278 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
 279 cgrad            enddo
 280 cgrad          enddo
 281 c 9/28/08 AL Gradient compotents will be summed only at the end
 282 C          ggg(1)=facvdw*xj
 283 C          ggg(2)=facvdw*yj
 284 C          ggg(3)=facvdw*zj
 285           ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
 286           ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
 287           ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
 288           do k=1,3
 289             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
 290             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
 291           enddo
 292 #endif
 293 *
 294 * Angular part
 295 *
 296           ecosa=2.0D0*fac3*fac1+fac4
 297           fac4=-3.0D0*fac4
 298           fac3=-6.0D0*fac3
 299           ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
 300           ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
 301           do k=1,3
 302             dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
 303             dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
 304           enddo
 305 cd        print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
 306 cd   &          (dcosg(k),k=1,3)
 307           do k=1,3
 308             ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
 309      &      fac_shield(i)**2*fac_shield(j)**2
 310
 311           enddo
 312 c          do k=1,3
 313 c            ghalf=0.5D0*ggg(k)
 314 c            gelc(k,i)=gelc(k,i)+ghalf
 315 c     &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
 316 c     &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
 317 c            gelc(k,j)=gelc(k,j)+ghalf
 318 c     &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
 319 c     &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
 320 c          enddo
 321 cgrad          do k=i+1,j-1
 322 cgrad            do l=1,3
 323 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
 324 cgrad            enddo
 325 cgrad          enddo
 326           do k=1,3
 327             gelc(k,i)=gelc(k,i)
 328      &               +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
 329      &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
 330      &           *fac_shield(i)**2*fac_shield(j)**2
 331
 332             gelc(k,j)=gelc(k,j)
 333      &               +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
 334      &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
 335      &           *fac_shield(i)**2*fac_shield(j)**2
 336             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
 337             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
 338           enddo
 339           IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
 340      &        .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
 341      &        .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
 342 C
 343 C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
 344 C   energy of a peptide unit is assumed in the form of a second-order
 345 C   Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
 346 C   Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
 347 C   are computed for EVERY pair of non-contiguous peptide groups.
 348 C
 349           if (j.lt.nres-1) then
 350             j1=j+1
 351             j2=j-1
 352           else
 353             j1=j-1
 354             j2=j-2
 355           endif
 356           kkk=0
 357           do k=1,2
 358             do l=1,2
 359               kkk=kkk+1
 360               muij(kkk)=mu(k,i)*mu(l,j)
 361             enddo
 362           enddo
 363 cd         write (iout,*) 'EELEC: i',i,' j',j
 364 cd          write (iout,*) 'j',j,' j1',j1,' j2',j2
 365 cd          write(iout,*) 'muij',muij
 366           ury=scalar(uy(1,i),erij)
 367           urz=scalar(uz(1,i),erij)
 368           vry=scalar(uy(1,j),erij)
 369           vrz=scalar(uz(1,j),erij)
 370           a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
 371           a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
 372           a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
 373           a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
 374           fac=dsqrt(-ael6i)*r3ij
 375           a22=a22*fac
 376           a23=a23*fac
 377           a32=a32*fac
 378           a33=a33*fac
 379 cd          write (iout,'(4i5,4f10.5)')
 380 cd     &     i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
 381 cd          write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
 382 cd          write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
 383 cd     &      uy(:,j),uz(:,j)
 384 cd          write (iout,'(4f10.5)')
 385 cd     &      scalar(uy(1,i),uy(1,j)),scalar(uy(1,i),uz(1,j)),
 386 cd     &      scalar(uz(1,i),uy(1,j)),scalar(uz(1,i),uz(1,j))
 387 cd          write (iout,'(4f10.5)') ury,urz,vry,vrz
 388 cd           write (iout,'(9f10.5/)')
 389 cd     &      fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
 390 C Derivatives of the elements of A in virtual-bond vectors
 391           call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
 392           do k=1,3
 393             uryg(k,1)=scalar(erder(1,k),uy(1,i))
 394             uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
 395             uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
 396             urzg(k,1)=scalar(erder(1,k),uz(1,i))
 397             urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
 398             urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
 399             vryg(k,1)=scalar(erder(1,k),uy(1,j))
 400             vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
 401             vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
 402             vrzg(k,1)=scalar(erder(1,k),uz(1,j))
 403             vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
 404             vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
 405           enddo
 406 C Compute radial contributions to the gradient
 407           facr=-3.0d0*rrmij
 408           a22der=a22*facr
 409           a23der=a23*facr
 410           a32der=a32*facr
 411           a33der=a33*facr
 412           agg(1,1)=a22der*xj
 413           agg(2,1)=a22der*yj
 414           agg(3,1)=a22der*zj
 415           agg(1,2)=a23der*xj
 416           agg(2,2)=a23der*yj
 417           agg(3,2)=a23der*zj
 418           agg(1,3)=a32der*xj
 419           agg(2,3)=a32der*yj
 420           agg(3,3)=a32der*zj
 421           agg(1,4)=a33der*xj
 422           agg(2,4)=a33der*yj
 423           agg(3,4)=a33der*zj
 424 C Add the contributions coming from er
 425           fac3=-3.0d0*fac
 426           do k=1,3
 427             agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
 428             agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
 429             agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
 430             agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
 431           enddo
 432           do k=1,3
 433 C Derivatives in DC(i)
 434 cgrad            ghalf1=0.5d0*agg(k,1)
 435 cgrad            ghalf2=0.5d0*agg(k,2)
 436 cgrad            ghalf3=0.5d0*agg(k,3)
 437 cgrad            ghalf4=0.5d0*agg(k,4)
 438             aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
 439      &      -3.0d0*uryg(k,2)*vry)!+ghalf1
 440             aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
 441      &      -3.0d0*uryg(k,2)*vrz)!+ghalf2
 442             aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
 443      &      -3.0d0*urzg(k,2)*vry)!+ghalf3
 444             aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
 445      &      -3.0d0*urzg(k,2)*vrz)!+ghalf4
 446 C Derivatives in DC(i+1)
 447             aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
 448      &      -3.0d0*uryg(k,3)*vry)!+agg(k,1)
 449             aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
 450      &      -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
 451             aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
 452      &      -3.0d0*urzg(k,3)*vry)!+agg(k,3)
 453             aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
 454      &      -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
 455 C Derivatives in DC(j)
 456             aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
 457      &      -3.0d0*vryg(k,2)*ury)!+ghalf1
 458             aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
 459      &      -3.0d0*vrzg(k,2)*ury)!+ghalf2
 460             aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
 461      &      -3.0d0*vryg(k,2)*urz)!+ghalf3
 462             aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
 463      &      -3.0d0*vrzg(k,2)*urz)!+ghalf4
 464 C Derivatives in DC(j+1) or DC(nres-1)
 465             aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
 466      &      -3.0d0*vryg(k,3)*ury)
 467             aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
 468      &      -3.0d0*vrzg(k,3)*ury)
 469             aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
 470      &      -3.0d0*vryg(k,3)*urz)
 471             aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
 472      &      -3.0d0*vrzg(k,3)*urz)
 473 cgrad            if (j.eq.nres-1 .and. i.lt.j-2) then
 474 cgrad              do l=1,4
 475 cgrad                aggj1(k,l)=aggj1(k,l)+agg(k,l)
 476 cgrad              enddo
 477 cgrad            endif
 478           enddo
 479           acipa(1,1)=a22
 480           acipa(1,2)=a23
 481           acipa(2,1)=a32
 482           acipa(2,2)=a33
 483           a22=-a22
 484           a23=-a23
 485           do l=1,2
 486             do k=1,3
 487               agg(k,l)=-agg(k,l)
 488               aggi(k,l)=-aggi(k,l)
 489               aggi1(k,l)=-aggi1(k,l)
 490               aggj(k,l)=-aggj(k,l)
 491               aggj1(k,l)=-aggj1(k,l)
 492             enddo
 493           enddo
 494           if (j.lt.nres-1) then
 495             a22=-a22
 496             a32=-a32
 497             do l=1,3,2
 498               do k=1,3
 499                 agg(k,l)=-agg(k,l)
 500                 aggi(k,l)=-aggi(k,l)
 501                 aggi1(k,l)=-aggi1(k,l)
 502                 aggj(k,l)=-aggj(k,l)
 503                 aggj1(k,l)=-aggj1(k,l)
 504               enddo
 505             enddo
 506           else
 507             a22=-a22
 508             a23=-a23
 509             a32=-a32
 510             a33=-a33
 511             do l=1,4
 512               do k=1,3
 513                 agg(k,l)=-agg(k,l)
 514                 aggi(k,l)=-aggi(k,l)
 515                 aggi1(k,l)=-aggi1(k,l)
 516                 aggj(k,l)=-aggj(k,l)
 517                 aggj1(k,l)=-aggj1(k,l)
 518               enddo
 519             enddo
 520           endif
 521           ENDIF ! WCORR
 522           IF (wel_loc.gt.0.0d0) THEN
 523 C Contribution to the local-electrostatic energy coming from the i-j pair
 524           eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
 525      &     +a33*muij(4)
 526 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
 527
 528           if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
 529      &            'eelloc',i,j,eel_loc_ij
 530
 531
 532           if (shield_mode.eq.0) then
 533            fac_shield(i)=1.0
 534            fac_shield(j)=1.0
 535 C          else
 536 C           fac_shield(i)=0.4
 537 C           fac_shield(j)=0.6
 538           endif
 539           eel_loc_ij=eel_loc_ij
 540      &    *fac_shield(i)*fac_shield(j)
 541           eel_loc=eel_loc+eel_loc_ij
 542
 543           if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
 544      &  (shield_mode.gt.0)) then
 545 C          print *,i,j
 546
 547           do ilist=1,ishield_list(i)
 548            iresshield=shield_list(ilist,i)
 549            do k=1,3
 550            rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
 551      &                                          /fac_shield(i)
 552 C     &      *2.0
 553            gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
 554      &              rlocshield
 555      & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
 556             gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
 557      &      +rlocshield
 558            enddo
 559           enddo
 560           do ilist=1,ishield_list(j)
 561            iresshield=shield_list(ilist,j)
 562            do k=1,3
 563            rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
 564      &                                       /fac_shield(j)
 565 C     &     *2.0
 566            gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
 567      &              rlocshield
 568      & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
 569            gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
 570      &             +rlocshield
 571
 572            enddo
 573           enddo
 574
 575           do k=1,3
 576             gshieldc_ll(k,i)=gshieldc_ll(k,i)+
 577      &              grad_shield(k,i)*eel_loc_ij/fac_shield(i)
 578             gshieldc_ll(k,j)=gshieldc_ll(k,j)+
 579      &              grad_shield(k,j)*eel_loc_ij/fac_shield(j)
 580             gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
 581      &              grad_shield(k,i)*eel_loc_ij/fac_shield(i)
 582             gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
 583      &              grad_shield(k,j)*eel_loc_ij/fac_shield(j)
 584            enddo
 585            endif
 586
 587 C Partial derivatives in virtual-bond dihedral angles gamma
 588           if (i.gt.1)
 589      &    gel_loc_loc(i-1)=gel_loc_loc(i-1)+
 590      &            (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
 591      &           +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
 592      &    *fac_shield(i)*fac_shield(j)
 593
 594           gel_loc_loc(j-1)=gel_loc_loc(j-1)+
 595      &            (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
 596      &           +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
 597      &    *fac_shield(i)*fac_shield(j)
 598
 599 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
 600           do l=1,3
 601             ggg(l)=(agg(l,1)*muij(1)+
 602      &          agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
 603      &    *fac_shield(i)*fac_shield(j)
 604
 605             gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
 606             gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
 607 cgrad            ghalf=0.5d0*ggg(l)
 608 cgrad            gel_loc(l,i)=gel_loc(l,i)+ghalf
 609 cgrad            gel_loc(l,j)=gel_loc(l,j)+ghalf
 610           enddo
 611 cgrad          do k=i+1,j2
 612 cgrad            do l=1,3
 613 cgrad              gel_loc(l,k)=gel_loc(l,k)+ggg(l)
 614 cgrad            enddo
 615 cgrad          enddo
 616 C Remaining derivatives of eello
 617           do l=1,3
 618             gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
 619      &          aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
 620      &    *fac_shield(i)*fac_shield(j)
 621
 622             gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
 623      &         aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
 624      &    *fac_shield(i)*fac_shield(j)
 625
 626             gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
 627      &          aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
 628      &    *fac_shield(i)*fac_shield(j)
 629
 630             gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
 631      &         aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
 632      &    *fac_shield(i)*fac_shield(j)
 633
 634           enddo
 635           ENDIF
 636 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
 637 c          if (j.gt.i+1 .and. num_conti.le.maxconts) then
 638           if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
 639      &       .and. num_conti.le.maxconts) then
 640 c            write (iout,*) i,j," entered corr"
 641 C
 642 C Calculate the contact function. The ith column of the array JCONT will
 643 C contain the numbers of atoms that make contacts with the atom I (of numbers
 644 C greater than I). The arrays FACONT and GACONT will contain the values of
 645 C the contact function and its derivative.
 646 c           r0ij=1.02D0*rpp(iteli,itelj)
 647 c           r0ij=1.11D0*rpp(iteli,itelj)
 648             r0ij=2.20D0*rpp(iteli,itelj)
 649 c           r0ij=1.55D0*rpp(iteli,itelj)
 650             call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
 651             if (fcont.gt.0.0D0) then
 652               num_conti=num_conti+1
 653               if (num_conti.gt.maxconts) then
 654                 write (iout,*) 'WARNING - max. # of contacts exceeded;',
 655      &                         ' will skip next contacts for this conf.'
 656               else
 657                 jcont_hb(num_conti,i)=j
 658 cd                write (iout,*) "i",i," j",j," num_conti",num_conti,
 659 cd     &           " jcont_hb",jcont_hb(num_conti,i)
 660                 IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
 661      &          wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
 662 C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
 663 C  terms.
 664                 d_cont(num_conti,i)=rij
 665 cd                write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
 666 C     --- Electrostatic-interaction matrix ---
 667                 a_chuj(1,1,num_conti,i)=a22
 668                 a_chuj(1,2,num_conti,i)=a23
 669                 a_chuj(2,1,num_conti,i)=a32
 670                 a_chuj(2,2,num_conti,i)=a33
 671 C     --- Gradient of rij
 672                 do kkk=1,3
 673                   grij_hb_cont(kkk,num_conti,i)=erij(kkk)
 674                 enddo
 675                 kkll=0
 676                 do k=1,2
 677                   do l=1,2
 678                     kkll=kkll+1
 679                     do m=1,3
 680                       a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
 681                       a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
 682                       a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
 683                       a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
 684                       a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
 685                     enddo
 686                   enddo
 687                 enddo
 688                 ENDIF
 689                 IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
 690 C Calculate contact energies
 691                 cosa4=4.0D0*cosa
 692                 wij=cosa-3.0D0*cosb*cosg
 693                 cosbg1=cosb+cosg
 694                 cosbg2=cosb-cosg
 695 c               fac3=dsqrt(-ael6i)/r0ij**3
 696                 fac3=dsqrt(-ael6i)*r3ij
 697 c                 ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
 698                 ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
 699                 if (ees0tmp.gt.0) then
 700                   ees0pij=dsqrt(ees0tmp)
 701                 else
 702                   ees0pij=0
 703                 endif
 704 c                ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
 705                 ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
 706                 if (ees0tmp.gt.0) then
 707                   ees0mij=dsqrt(ees0tmp)
 708                 else
 709                   ees0mij=0
 710                 endif
 711 c               ees0mij=0.0D0
 712                 if (shield_mode.eq.0) then
 713                 fac_shield(i)=1.0d0
 714                 fac_shield(j)=1.0d0
 715                 else
 716                 ees0plist(num_conti,i)=j
 717 C                fac_shield(i)=0.4d0
 718 C                fac_shield(j)=0.6d0
 719                 endif
 720                 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
 721      &          *fac_shield(i)*fac_shield(j)
 722                 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
 723      &          *fac_shield(i)*fac_shield(j)
 724
 725 C Diagnostics. Comment out or remove after debugging!
 726 c               ees0p(num_conti,i)=0.5D0*fac3*ees0pij
 727 c               ees0m(num_conti,i)=0.5D0*fac3*ees0mij
 728 c               ees0m(num_conti,i)=0.0D0
 729 C End diagnostics.
 730 c               write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
 731 c    & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
 732 C Angular derivatives of the contact function
 733                 ees0pij1=fac3/ees0pij
 734                 ees0mij1=fac3/ees0mij
 735                 fac3p=-3.0D0*fac3*rrmij
 736                 ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
 737                 ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
 738 c               ees0mij1=0.0D0
 739                 ecosa1=       ees0pij1*( 1.0D0+0.5D0*wij)
 740                 ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
 741                 ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
 742                 ecosa2=       ees0mij1*(-1.0D0+0.5D0*wij)
 743                 ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
 744                 ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
 745                 ecosap=ecosa1+ecosa2
 746                 ecosbp=ecosb1+ecosb2
 747                 ecosgp=ecosg1+ecosg2
 748                 ecosam=ecosa1-ecosa2
 749                 ecosbm=ecosb1-ecosb2
 750                 ecosgm=ecosg1-ecosg2
 751 C Diagnostics
 752 c               ecosap=ecosa1
 753 c               ecosbp=ecosb1
 754 c               ecosgp=ecosg1
 755 c               ecosam=0.0D0
 756 c               ecosbm=0.0D0
 757 c               ecosgm=0.0D0
 758 C End diagnostics
 759                 facont_hb(num_conti,i)=fcont
 760                 fprimcont=fprimcont/rij
 761 cd              facont_hb(num_conti,i)=1.0D0
 762 C Following line is for diagnostics.
 763 cd              fprimcont=0.0D0
 764                 do k=1,3
 765                   dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
 766                   dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
 767                 enddo
 768                 do k=1,3
 769                   gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
 770                   gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
 771                 enddo
 772                 gggp(1)=gggp(1)+ees0pijp*xj
 773                 gggp(2)=gggp(2)+ees0pijp*yj
 774                 gggp(3)=gggp(3)+ees0pijp*zj
 775                 gggm(1)=gggm(1)+ees0mijp*xj
 776                 gggm(2)=gggm(2)+ees0mijp*yj
 777                 gggm(3)=gggm(3)+ees0mijp*zj
 778 C Derivatives due to the contact function
 779                 gacont_hbr(1,num_conti,i)=fprimcont*xj
 780                 gacont_hbr(2,num_conti,i)=fprimcont*yj
 781                 gacont_hbr(3,num_conti,i)=fprimcont*zj
 782                 do k=1,3
 783 c
 784 c 10/24/08 cgrad and ! comments indicate the parts of the code removed
 785 c          following the change of gradient-summation algorithm.
 786 c
 787 cgrad                  ghalfp=0.5D0*gggp(k)
 788 cgrad                  ghalfm=0.5D0*gggm(k)
 789                   gacontp_hb1(k,num_conti,i)=!ghalfp
 790      &              +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
 791      &              + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
 792      &          *fac_shield(i)*fac_shield(j)
 793
 794                   gacontp_hb2(k,num_conti,i)=!ghalfp
 795      &              +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
 796      &              + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
 797      &          *fac_shield(i)*fac_shield(j)
 798
 799                   gacontp_hb3(k,num_conti,i)=gggp(k)
 800      &          *fac_shield(i)*fac_shield(j)
 801
 802                   gacontm_hb1(k,num_conti,i)=!ghalfm
 803      &              +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
 804      &              + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
 805      &          *fac_shield(i)*fac_shield(j)
 806
 807                   gacontm_hb2(k,num_conti,i)=!ghalfm
 808      &              +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
 809      &              + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
 810      &          *fac_shield(i)*fac_shield(j)
 811
 812                   gacontm_hb3(k,num_conti,i)=gggm(k)
 813      &          *fac_shield(i)*fac_shield(j)
 814
 815                 enddo
 816               ENDIF ! wcorr
 817               endif  ! num_conti.le.maxconts
 818             endif  ! fcont.gt.0
 819           endif    ! j.gt.i+1
 820           if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
 821             do k=1,4
 822               do l=1,3
 823                 ghalf=0.5d0*agg(l,k)
 824                 aggi(l,k)=aggi(l,k)+ghalf
 825                 aggi1(l,k)=aggi1(l,k)+agg(l,k)
 826                 aggj(l,k)=aggj(l,k)+ghalf
 827               enddo
 828             enddo
 829             if (j.eq.nres-1 .and. i.lt.j-2) then
 830               do k=1,4
 831                 do l=1,3
 832                   aggj1(l,k)=aggj1(l,k)+agg(l,k)
 833                 enddo
 834               enddo
 835             endif
 836           endif
 837 c          t_eelecij=t_eelecij+MPI_Wtime()-time00
 838       return
 839       end