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

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