+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+ dsci_inv=vbld_inv(i+nres)
+!C
+!C Calculate SC interaction energy.
+!C
+ do iint=1,nint_gr_nucl(i)
+! print *,"tu?",i,istart_nucl(i,iint),iend_nucl(i,iint)
+ do j=istart_nucl(i,iint),iend_nucl(i,iint)
+ ind=ind+1
+! print *,"JESTEM"
+ itypj=itype(j,2)
+ if (itypj.eq.ntyp1_molec(2)) cycle
+ dscj_inv=vbld_inv(j+nres)
+ sig0ij=sigma_nucl(itypi,itypj)
+ chi1=chi_nucl(itypi,itypj)
+ chi2=chi_nucl(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip_nucl(itypi,itypj)
+ chip2=chip_nucl(itypj,itypi)
+ chip12=chip1*chip2
+! xj=c(1,nres+j)-xi
+! yj=c(2,nres+j)-yi
+! zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+!C Calculate angle-dependent terms of energy and contributions to their
+!C derivatives.
+ erij(1)=xj*rij
+ erij(2)=yj*rij
+ erij(3)=zj*rij
+ om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
+ om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
+ om12=dxi*dxj+dyi*dyj+dzi*dzj
+ call sc_angular_nucl
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+sig0ij
+! print *,rij_shift,"rij_shift"
+!c write (2,*) " rij",1.0D0/rij," sig",sig," sig0ij",sig0ij,
+!c & " rij_shift",rij_shift
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D20
+ return
+ endif
+ sigder=-sig*sigsq
+!c---------------------------------------------------------------
+ rij_shift=1.0D0/rij_shift
+ fac=rij_shift**expon
+ e1=fac*fac*aa_nucl(itypi,itypj)
+ e2=fac*bb_nucl(itypi,itypj)
+ evdwij=eps1*eps2rt*(e1+e2)
+!c write (2,*) "eps1",eps1," eps2rt",eps2rt,
+!c & " e1",e1," e2",e2," evdwij",evdwij
+ eps2der=evdwij
+ evdwij=evdwij*eps2rt
+ evdwsb=evdwsb+evdwij
+ if (lprn) then
+ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
+ restyp(itypi,2),i,restyp(itypj,2),j, &
+ epsi,sigm,chi1,chi2,chip1,chip2, &
+ eps1,eps2rt**2,sig,sig0ij, &
+ om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,&
+ evdwij
+ write (iout,*) "aa",aa_nucl(itypi,itypj)," bb",bb_nucl(itypi,itypj)
+ endif
+
+ if (energy_dec) write (iout,'(a6,2i5,e15.3,a4)') &
+ 'evdw',i,j,evdwij,"tu3"
+
+
+!C Calculate gradient components.
+ e1=e1*eps1*eps2rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ fac=rij*fac
+!c fac=0.0d0
+!C Calculate the radial part of the gradient
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+!C Calculate angular part of the gradient.
+ call sc_grad_nucl
+ call eelsbij(eelij,num_conti2)
+ if (energy_dec .and. &
+ (j.eq.i+1.or.j.eq.nres-i+1.or.j.eq.nres-i.or.j.eq.nres-i+2)) &
+ write (istat,'(e14.5)') evdwij
+ eelsb=eelsb+eelij
+ enddo ! j
+ enddo ! iint
+ num_cont_hb(i)=num_conti2
+ enddo ! i
+!c write (iout,*) "Number of loop steps in EGB:",ind
+!cccc energy_dec=.false.
+ return
+ end subroutine esb_gb
+!-------------------------------------------------------------------------------
+ subroutine eelsbij(eesij,num_conti2)
+ use comm_locel
+ use calc_data_nucl
+ real(kind=8),dimension(3) :: ggg,gggp,gggm,dcosb,dcosg
+ real(kind=8),dimension(3,3) :: erder,uryg,urzg,vryg,vrzg
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,rlocshield,fracinbuf
+ integer xshift,yshift,zshift,ilist,iresshield,num_conti2
+
+!c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
+ real(kind=8) scal_el /0.5d0/
+ integer :: iteli,itelj,kkk,kkll,m,isubchap
+ real(kind=8) :: ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp,facfac
+ real(kind=8) :: ees,evdw1,eel_loc,aaa,bbb,ael3i,ael63i,ael32i
+ real(kind=8) :: dx_normj,dy_normj,dz_normj,&
+ r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,fac5,fac6,&
+ el1,el2,el3,el4,eesij,ees0ij,facvdw,facel,fac1,ecosa,&
+ ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der,&
+ a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij,&
+ ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp,&
+ ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,&
+ ecosgp,ecosam,ecosbm,ecosgm,ghalf,itypi,itypj
+ ind=ind+1
+ itypi=itype(i,2)
+ itypj=itype(j,2)
+! print *,i,j,itypi,itypj,istype(i),istype(j),"????"
+ ael6i=ael6_nucl(itypi,itypj)
+ ael3i=ael3_nucl(itypi,itypj)
+ ael63i=ael63_nucl(itypi,itypj)
+ ael32i=ael32_nucl(itypi,itypj)
+!c write (iout,*) "eelecij",i,j,itype(i),itype(j),
+!c & ael6i,ael3i,ael63i,al32i,rij,rrij
+ dxj=dc(1,j+nres)
+ dyj=dc(2,j+nres)
+ dzj=dc(3,j+nres)
+ dx_normi=dc_norm(1,i+nres)
+ dy_normi=dc_norm(2,i+nres)
+ dz_normi=dc_norm(3,i+nres)
+ dx_normj=dc_norm(1,j+nres)
+ dy_normj=dc_norm(2,j+nres)
+ dz_normj=dc_norm(3,j+nres)
+!c xj=c(1,j)+0.5D0*dxj-xmedi
+!c yj=c(2,j)+0.5D0*dyj-ymedi
+!c zj=c(3,j)+0.5D0*dzj-zmedi
+ if (ipot_nucl.ne.2) then
+ cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
+ cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
+ cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
+ else
+ cosa=om12
+ cosb=om1
+ cosg=om2
+ endif
+ r3ij=rij*rrij
+ r6ij=r3ij*r3ij
+ fac=cosa-3.0D0*cosb*cosg
+ facfac=fac*fac
+ fac1=3.0d0*(cosb*cosb+cosg*cosg)
+ fac3=ael6i*r6ij
+ fac4=ael3i*r3ij
+ fac5=ael63i*r6ij
+ fac6=ael32i*r6ij
+!c write (iout,*) "r3ij",r3ij," r6ij",r6ij," fac",fac," fac1",fac1,
+!c & " fac2",fac2," fac3",fac3," fac4",fac4," fac5",fac5," fac6",fac6
+ el1=fac3*(4.0D0+facfac-fac1)
+ el2=fac4*fac
+ el3=fac5*(2.0d0-2.0d0*facfac+fac1)
+ el4=fac6*facfac
+ eesij=el1+el2+el3+el4
+!C 12/26/95 - for the evaluation of multi-body H-bonding interactions
+ ees0ij=4.0D0+facfac-fac1
+
+ if (energy_dec) then
+ if(j.eq.i+1.or.j.eq.nres-i+1.or.j.eq.nres-i.or.j.eq.nres-i+2) &
+ write (istat,'(2a1,i4,1x,2a1,i4,4f10.5,3e12.5,$)') &
+ sugartyp(istype(i)),restyp(itypi,2),i,sugartyp(istype(j)),&
+ restyp(itypj,2),j,1.0d0/rij,cosa,cosb,cosg,fac*r3ij, &
+ (4.0D0+facfac-fac1)*r6ij,(2.0d0-2.0d0*facfac+fac1)*r6ij
+ write (iout,'(a6,2i5,e15.3)') 'ees',i,j,eesij
+ endif
+
+!C
+!C Calculate contributions to the Cartesian gradient.
+!C
+ facel=-3.0d0*rrij*(eesij+el1+el3+el4)
+ fac1=fac
+!c erij(1)=xj*rmij
+!c erij(2)=yj*rmij
+!c erij(3)=zj*rmij
+!*
+!* Radial derivatives. First process both termini of the fragment (i,j)
+!*
+ ggg(1)=facel*xj
+ ggg(2)=facel*yj
+ ggg(3)=facel*zj
+ do k=1,3
+ gelsbc(k,j)=gelsbc(k,j)+ggg(k)
+ gelsbc(k,i)=gelsbc(k,i)-ggg(k)
+ gelsbx(k,j)=gelsbx(k,j)+ggg(k)
+ gelsbx(k,i)=gelsbx(k,i)-ggg(k)
+ enddo
+!*
+!* Angular part
+!*
+ ecosa=2.0D0*fac3*fac1+fac4+(-4.0d0*fac5+2.0d0*fac6)*fac1
+ fac4=-3.0D0*fac4
+ fac3=-6.0D0*fac3
+ fac5= 6.0d0*fac5
+ fac6=-6.0d0*fac6
+ ecosb=fac3*(fac1*cosg+cosb)+cosg*fac4+(cosb+2*fac1*cosg)*fac5+&
+ fac6*fac1*cosg
+ ecosg=fac3*(fac1*cosb+cosg)+cosb*fac4+(cosg+2*fac1*cosb)*fac5+&
+ fac6*fac1*cosb
+ do k=1,3
+ dcosb(k)=rij*(dc_norm(k,i+nres)-erij(k)*cosb)
+ dcosg(k)=rij*(dc_norm(k,j+nres)-erij(k)*cosg)
+ enddo
+ do k=1,3
+ ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
+ enddo
+ do k=1,3
+ gelsbx(k,i)=gelsbx(k,i)-ggg(k) &
+ +(ecosa*(dc_norm(k,j+nres)-cosa*dc_norm(k,i+nres))&
+ + ecosb*(erij(k)-cosb*dc_norm(k,i+nres)))*vbld_inv(i+nres)
+ gelsbx(k,j)=gelsbx(k,j)+ggg(k) &
+ +(ecosa*(dc_norm(k,i+nres)-cosa*dc_norm(k,j+nres))&
+ + ecosg*(erij(k)-cosg*dc_norm(k,j+nres)))*vbld_inv(j+nres)
+ gelsbc(k,j)=gelsbc(k,j)+ggg(k)
+ gelsbc(k,i)=gelsbc(k,i)-ggg(k)
+ enddo
+! IF ( (wcorr_nucl.gt.0.0d0.or.wcorr3_nucl.gt.0.0d0) .and.
+ IF ( j.gt.i+1 .and.&
+ num_conti.le.maxconts) THEN
+!C
+!C Calculate the contact function. The ith column of the array JCONT will
+!C contain the numbers of atoms that make contacts with the atom I (of numbers
+!C greater than I). The arrays FACONT and GACONT will contain the values of
+!C the contact function and its derivative.
+ r0ij=2.20D0*sigma(itypi,itypj)
+!c write (2,*) "ij",i,j," rij",1.0d0/rij," r0ij",r0ij
+ call gcont(rij,r0ij,1.0D0,0.2d0/r0ij,fcont,fprimcont)
+!c write (2,*) "fcont",fcont
+ if (fcont.gt.0.0D0) then
+ num_conti=num_conti+1
+ num_conti2=num_conti2+1
+
+ if (num_conti.gt.maxconts) then
+ write (iout,*) 'WARNING - max. # of contacts exceeded;',&
+ ' will skip next contacts for this conf.'
+ else
+ jcont_hb(num_conti,i)=j
+!c write (iout,*) "num_conti",num_conti,
+!c & " jcont_hb",jcont_hb(num_conti,i)
+!C Calculate contact energies
+ cosa4=4.0D0*cosa
+ wij=cosa-3.0D0*cosb*cosg
+ cosbg1=cosb+cosg
+ cosbg2=cosb-cosg
+ fac3=dsqrt(-ael6i)*r3ij
+!c write (2,*) "ael6i",ael6i," r3ij",r3ij," fac3",fac3
+ ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
+ if (ees0tmp.gt.0) then
+ ees0pij=dsqrt(ees0tmp)
+ else
+ ees0pij=0
+ endif
+ ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
+ if (ees0tmp.gt.0) then
+ ees0mij=dsqrt(ees0tmp)
+ else
+ ees0mij=0
+ endif
+ ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
+ ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
+!c write (iout,*) "i",i," j",j,
+!c & " ees0m",ees0m(num_conti,i)," ees0p",ees0p(num_conti,i)
+ ees0pij1=fac3/ees0pij
+ ees0mij1=fac3/ees0mij
+ fac3p=-3.0D0*fac3*rrij
+ ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
+ ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
+ ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
+ ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
+ ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
+ ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
+ ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
+ ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
+ ecosap=ecosa1+ecosa2
+ ecosbp=ecosb1+ecosb2
+ ecosgp=ecosg1+ecosg2
+ ecosam=ecosa1-ecosa2
+ ecosbm=ecosb1-ecosb2
+ ecosgm=ecosg1-ecosg2
+!C End diagnostics
+ facont_hb(num_conti,i)=fcont
+ fprimcont=fprimcont/rij
+ do k=1,3
+ gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
+ gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
+ enddo
+ gggp(1)=gggp(1)+ees0pijp*xj
+ gggp(2)=gggp(2)+ees0pijp*yj
+ gggp(3)=gggp(3)+ees0pijp*zj
+ gggm(1)=gggm(1)+ees0mijp*xj
+ gggm(2)=gggm(2)+ees0mijp*yj
+ gggm(3)=gggm(3)+ees0mijp*zj
+!C Derivatives due to the contact function
+ gacont_hbr(1,num_conti,i)=fprimcont*xj
+ gacont_hbr(2,num_conti,i)=fprimcont*yj
+ gacont_hbr(3,num_conti,i)=fprimcont*zj
+ do k=1,3
+!c
+!c Gradient of the correlation terms
+!c
+ gacontp_hb1(k,num_conti,i)= &
+ (ecosap*(dc_norm(k,j+nres)-cosa*dc_norm(k,i+nres)) &
+ + ecosbp*(erij(k)-cosb*dc_norm(k,i+nres)))*vbld_inv(i+nres)
+ gacontp_hb2(k,num_conti,i)= &
+ (ecosap*(dc_norm(k,i+nres)-cosa*dc_norm(k,j+nres)) &
+ + ecosgp*(erij(k)-cosg*dc_norm(k,j+nres)))*vbld_inv(j+nres)
+ gacontp_hb3(k,num_conti,i)=gggp(k)
+ gacontm_hb1(k,num_conti,i)= &
+ (ecosam*(dc_norm(k,j+nres)-cosa*dc_norm(k,i+nres)) &
+ + ecosbm*(erij(k)-cosb*dc_norm(k,i+nres)))*vbld_inv(i+nres)
+ gacontm_hb2(k,num_conti,i)= &
+ (ecosam*(dc_norm(k,i+nres)-cosa*dc_norm(k,j+nres))&
+ + ecosgm*(erij(k)-cosg*dc_norm(k,j+nres)))*vbld_inv(j+nres)
+ gacontm_hb3(k,num_conti,i)=gggm(k)
+ enddo
+ endif
+ endif
+ ENDIF
+ return
+ end subroutine eelsbij
+!------------------------------------------------------------------
+ subroutine sc_grad_nucl
+ use comm_locel
+ use calc_data_nucl
+ real(kind=8),dimension(3) :: dcosom1,dcosom2
+ eom1=eps2der*eps2rt_om1+sigder*sigsq_om1
+ eom2=eps2der*eps2rt_om2+sigder*sigsq_om2
+ eom12=evdwij*eps1_om12+eps2der*eps2rt_om12+sigder*sigsq_om12
+ do k=1,3
+ dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
+ dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
+ enddo
+ do k=1,3
+ gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
+ enddo
+ do k=1,3
+ gvdwsbx(k,i)=gvdwsbx(k,i)-gg(k) &
+ +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))&
+ +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+ gvdwsbx(k,j)=gvdwsbx(k,j)+gg(k) &
+ +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
+ +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ enddo
+!C
+!C Calculate the components of the gradient in DC and X
+!C
+ do l=1,3
+ gvdwsbc(l,i)=gvdwsbc(l,i)-gg(l)
+ gvdwsbc(l,j)=gvdwsbc(l,j)+gg(l)
+ enddo
+ return
+ end subroutine sc_grad_nucl
+!-----------------------------------------------------------------------
+ subroutine esb(esbloc)
+!C Calculate the local energy of a side chain and its derivatives in the
+!C corresponding virtual-bond valence angles THETA and the spherical angles
+!C ALPHA and OMEGA derived from AM1 all-atom calculations.
+!C added by Urszula Kozlowska. 07/11/2007
+!C
+ real(kind=8),dimension(3):: x_prime,y_prime,z_prime
+ real(kind=8),dimension(9):: x
+ real(kind=8) :: sumene,dsc_i,dp2_i,xx,yy,zz,sumene1, &
+ sumene2,sumene3,sumene4,s1,s1_6,s2,s2_6,&
+ de_dxx,de_dyy,de_dzz,de_dt,s1_t,s1_6_t,s2_t,s2_6_t
+ real(kind=8),dimension(3):: dXX_Ci1,dYY_Ci1,dZZ_Ci1,dXX_Ci,&
+ dYY_Ci,dZZ_Ci,dXX_XYZ,dYY_XYZ,dZZ_XYZ,dt_dCi,dt_dCi1
+ real(kind=8) :: esbloc,delta,cosfac2,cosfac,sinfac2,sinfac,de_dtt,&
+ cossc,cossc1,cosfac2xx,sinfac2yy,pom1,pom
+ integer::it,nlobit,i,j,k
+! common /sccalc/ time11,time12,time112,theti,it,nlobit
+ delta=0.02d0*pi
+ esbloc=0.0D0
+ do i=loc_start_nucl,loc_end_nucl
+ if (itype(i,2).eq.ntyp1_molec(2)) cycle
+ costtab(i+1) =dcos(theta(i+1))
+ sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
+ cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
+ sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
+ cosfac2=0.5d0/(1.0d0+costtab(i+1))
+ cosfac=dsqrt(cosfac2)
+ sinfac2=0.5d0/(1.0d0-costtab(i+1))
+ sinfac=dsqrt(sinfac2)
+ it=itype(i,2)
+ if (it.eq.10) goto 1
+
+!c
+!C Compute the axes of tghe local cartesian coordinates system; store in
+!c x_prime, y_prime and z_prime
+!c
+ do j=1,3
+ x_prime(j) = 0.00
+ y_prime(j) = 0.00
+ z_prime(j) = 0.00
+ enddo
+!C write(2,*) "dc_norm", dc_norm(1,i+nres),dc_norm(2,i+nres),
+!C & dc_norm(3,i+nres)
+ do j = 1,3
+ x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
+ y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
+ enddo
+ do j = 1,3
+ z_prime(j) = -uz(j,i-1)
+ enddo
+
+ xx=0.0d0
+ yy=0.0d0
+ zz=0.0d0
+ do j = 1,3
+ xx = xx + x_prime(j)*dc_norm(j,i+nres)
+ yy = yy + y_prime(j)*dc_norm(j,i+nres)
+ zz = zz + z_prime(j)*dc_norm(j,i+nres)
+ enddo
+
+ xxtab(i)=xx
+ yytab(i)=yy
+ zztab(i)=zz
+ it=itype(i,2)
+ do j = 1,9
+ x(j) = sc_parmin_nucl(j,it)
+ enddo
+#ifdef CHECK_COORD
+!Cc diagnostics - remove later
+ xx1 = dcos(alph(2))
+ yy1 = dsin(alph(2))*dcos(omeg(2))
+ zz1 = -dsin(alph(2))*dsin(omeg(2))
+ write(2,'(3f8.1,3f9.3,1x,3f9.3)') &
+ alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,&
+ xx1,yy1,zz1
+!C," --- ", xx_w,yy_w,zz_w
+!c end diagnostics
+#endif
+ sumene = enesc_nucl(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ esbloc = esbloc + sumene
+ if (energy_dec) write(iout,*) "i",i," esbloc",sumene,esbloc,xx,yy,zz
+ if (energy_dec) write(iout,*) "x",(x(k),k=1,9)
+#ifdef DEBUG
+ write (2,*) "x",(x(k),k=1,9)
+!C
+!C This section to check the numerical derivatives of the energy of ith side
+!C chain in xx, yy, zz, and theta. Use the -DDEBUG compiler option or insert
+!C #define DEBUG in the code to turn it on.
+!C
+ write (2,*) "sumene =",sumene
+ aincr=1.0d-7
+ xxsave=xx
+ xx=xx+aincr
+ write (2,*) xx,yy,zz
+ sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ de_dxx_num=(sumenep-sumene)/aincr
+ xx=xxsave
+ write (2,*) "xx+ sumene from enesc=",sumenep,sumene
+ yysave=yy
+ yy=yy+aincr
+ write (2,*) xx,yy,zz
+ sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ de_dyy_num=(sumenep-sumene)/aincr
+ yy=yysave
+ write (2,*) "yy+ sumene from enesc=",sumenep,sumene
+ zzsave=zz
+ zz=zz+aincr
+ write (2,*) xx,yy,zz
+ sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ de_dzz_num=(sumenep-sumene)/aincr
+ zz=zzsave
+ write (2,*) "zz+ sumene from enesc=",sumenep,sumene
+ costsave=cost2tab(i+1)
+ sintsave=sint2tab(i+1)
+ cost2tab(i+1)=dcos(0.5d0*(theta(i+1)+aincr))
+ sint2tab(i+1)=dsin(0.5d0*(theta(i+1)+aincr))
+ sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ de_dt_num=(sumenep-sumene)/aincr
+ write (2,*) " t+ sumene from enesc=",sumenep,sumene
+ cost2tab(i+1)=costsave
+ sint2tab(i+1)=sintsave
+!C End of diagnostics section.
+#endif
+!C
+!C Compute the gradient of esc
+!C
+ de_dxx=x(1)+2*x(4)*xx+x(7)*zz+x(8)*yy
+ de_dyy=x(2)+2*x(5)*yy+x(8)*xx+x(9)*zz
+ de_dzz=x(3)+2*x(6)*zz+x(7)*xx+x(9)*yy
+ de_dtt=0.0d0
+#ifdef DEBUG
+ write (2,*) "x",(x(k),k=1,9)
+ write (2,*) "xx",xx," yy",yy," zz",zz
+ write (2,*) "de_xx ",de_xx," de_yy ",de_yy,&
+ " de_zz ",de_zz," de_tt ",de_tt
+ write (2,*) "de_xx_num",de_dxx_num," de_yy_num",de_dyy_num,&
+ " de_zz_num",de_dzz_num," de_dt_num",de_dt_num
+#endif
+!C
+ cossc=scalar(dc_norm(1,i),dc_norm(1,i+nres))
+ cossc1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres))
+ cosfac2xx=cosfac2*xx
+ sinfac2yy=sinfac2*yy
+ do k = 1,3
+ dt_dCi(k) = -(dc_norm(k,i-1)+costtab(i+1)*dc_norm(k,i))*&
+ vbld_inv(i+1)
+ dt_dCi1(k)= -(dc_norm(k,i)+costtab(i+1)*dc_norm(k,i-1))*&
+ vbld_inv(i)
+ pom=(dC_norm(k,i+nres)-cossc*dC_norm(k,i))*vbld_inv(i+1)
+ pom1=(dC_norm(k,i+nres)-cossc1*dC_norm(k,i-1))*vbld_inv(i)
+!c write (iout,*) "i",i," k",k," pom",pom," pom1",pom1,
+!c & " dt_dCi",dt_dCi(k)," dt_dCi1",dt_dCi1(k)
+!c write (iout,*) "dC_norm",(dC_norm(j,i),j=1,3),
+!c & (dC_norm(j,i-1),j=1,3)," vbld_inv",vbld_inv(i+1),vbld_inv(i)
+ dXX_Ci(k)=pom*cosfac-dt_dCi(k)*cosfac2xx
+ dXX_Ci1(k)=-pom1*cosfac-dt_dCi1(k)*cosfac2xx
+ dYY_Ci(k)=pom*sinfac+dt_dCi(k)*sinfac2yy
+ dYY_Ci1(k)=pom1*sinfac+dt_dCi1(k)*sinfac2yy
+ dZZ_Ci1(k)=0.0d0
+ dZZ_Ci(k)=0.0d0
+ do j=1,3
+ dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres)
+ dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres)
+ enddo
+
+ dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
+ dYY_XYZ(k)=vbld_inv(i+nres)*(y_prime(k)-yy*dC_norm(k,i+nres))
+ dZZ_XYZ(k)=vbld_inv(i+nres)*(z_prime(k)-zz*dC_norm(k,i+nres))
+!c
+ dt_dCi(k) = -dt_dCi(k)/sinttab(i+1)
+ dt_dCi1(k)= -dt_dCi1(k)/sinttab(i+1)
+ enddo
+
+ do k=1,3
+ dXX_Ctab(k,i)=dXX_Ci(k)
+ dXX_C1tab(k,i)=dXX_Ci1(k)
+ dYY_Ctab(k,i)=dYY_Ci(k)
+ dYY_C1tab(k,i)=dYY_Ci1(k)
+ dZZ_Ctab(k,i)=dZZ_Ci(k)
+ dZZ_C1tab(k,i)=dZZ_Ci1(k)
+ dXX_XYZtab(k,i)=dXX_XYZ(k)
+ dYY_XYZtab(k,i)=dYY_XYZ(k)
+ dZZ_XYZtab(k,i)=dZZ_XYZ(k)
+ enddo
+ do k = 1,3
+!c write (iout,*) "k",k," dxx_ci1",dxx_ci1(k)," dyy_ci1",
+!c & dyy_ci1(k)," dzz_ci1",dzz_ci1(k)
+!c write (iout,*) "k",k," dxx_ci",dxx_ci(k)," dyy_ci",
+!c & dyy_ci(k)," dzz_ci",dzz_ci(k)
+!c write (iout,*) "k",k," dt_dci",dt_dci(k)," dt_dci",
+!c & dt_dci(k)
+!c write (iout,*) "k",k," dxx_XYZ",dxx_XYZ(k)," dyy_XYZ",
+!c & dyy_XYZ(k)," dzz_XYZ",dzz_XYZ(k)
+ gsbloc(k,i-1)=gsbloc(k,i-1)+de_dxx*dxx_ci1(k) &
+ +de_dyy*dyy_ci1(k)+de_dzz*dzz_ci1(k)+de_dt*dt_dCi1(k)
+ gsbloc(k,i)=gsbloc(k,i)+de_dxx*dxx_Ci(k) &
+ +de_dyy*dyy_Ci(k)+de_dzz*dzz_Ci(k)+de_dt*dt_dCi(k)
+ gsblocx(k,i)= de_dxx*dxx_XYZ(k)&
+ +de_dyy*dyy_XYZ(k)+de_dzz*dzz_XYZ(k)
+ enddo
+!c write(iout,*) "ENERGY GRAD = ", (gsbloc(k,i-1),k=1,3),
+!c & (gsbloc(k,i),k=1,3),(gsblocx(k,i),k=1,3)
+
+!C to check gradient call subroutine check_grad
+
+ 1 continue
+ enddo
+ return
+ end subroutine esb
+!=-------------------------------------------------------
+ real(kind=8) function enesc_nucl(x,xx,yy,zz,cost2,sint2)
+! implicit none
+ real(kind=8),dimension(9):: x(9)
+ real(kind=8) :: xx,yy,zz,cost2,sint2,sumene1,sumene2, &
+ sumene3,sumene4,sumene,dsc_i,dp2_i,dscp1,dscp2,s1,s1_6,s2,s2_6
+ integer i
+!c write (2,*) "enesc"
+!c write (2,*) "x",(x(i),i=1,9)
+!c write(2,*)"xx",xx," yy",yy," zz",zz," cost2",cost2," sint2",sint2
+ sumene=x(1)*xx+x(2)*yy+x(3)*zz+x(4)*xx**2 &
+ + x(5)*yy**2+x(6)*zz**2+x(7)*xx*zz+x(8)*xx*yy &
+ + x(9)*yy*zz
+ enesc_nucl=sumene
+ return
+ end function enesc_nucl
+!-----------------------------------------------------------------------------
+ subroutine multibody_hb_nucl(ecorr,ecorr3,n_corr,n_corr1)
+#ifdef MPI
+ include 'mpif.h'
+ integer,parameter :: max_cont=2000
+ integer,parameter:: max_dim=2*(8*3+6)
+ integer, parameter :: msglen1=max_cont*max_dim
+ integer,parameter :: msglen2=2*msglen1
+ integer source,CorrelType,CorrelID,Error
+ real(kind=8) :: buffer(max_cont,max_dim)
+ integer status(MPI_STATUS_SIZE)
+ integer :: ierror,nbytes
+#endif
+ real(kind=8),dimension(3):: gx(3),gx1(3)
+ real(kind=8) :: time00
+ logical lprn,ldone
+ integer i,j,i1,j1,jj,kk,num_conti,num_conti1,nn
+ real(kind=8) ecorr,ecorr3
+ integer :: n_corr,n_corr1,mm,msglen
+!C Set lprn=.true. for debugging
+ lprn=.false.
+ n_corr=0
+ n_corr1=0
+#ifdef MPI
+ if(.not.allocated(zapas2)) allocate(zapas2(3,maxconts,nres,8))
+
+ if (nfgtasks.le.1) goto 30
+ if (lprn) then
+ write (iout,'(a)') 'Contact function values:'
+ do i=nnt,nct-1
+ write (iout,'(2i3,50(1x,i2,f5.2))') &
+ i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), &
+ j=1,num_cont_hb(i))
+ enddo
+ endif
+!C Caution! Following code assumes that electrostatic interactions concerning
+!C a given atom are split among at most two processors!
+ CorrelType=477
+ CorrelID=fg_rank+1
+ ldone=.false.
+ do i=1,max_cont
+ do j=1,max_dim
+ buffer(i,j)=0.0D0
+ enddo
+ enddo
+ mm=mod(fg_rank,2)
+!c write (*,*) 'MyRank',MyRank,' mm',mm
+ if (mm) 20,20,10
+ 10 continue
+!c write (*,*) 'Sending: MyRank',MyRank,' mm',mm,' ldone',ldone
+ if (fg_rank.gt.0) then
+!C Send correlation contributions to the preceding processor
+ msglen=msglen1
+ nn=num_cont_hb(iatel_s_nucl)
+ call pack_buffer(max_cont,max_dim,iatel_s,0,buffer)
+!c write (*,*) 'The BUFFER array:'
+!c do i=1,nn
+!c write (*,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,30)
+!c enddo
+ if (ielstart_nucl(iatel_s_nucl).gt.iatel_s_nucl+ispp) then
+ msglen=msglen2
+ call pack_buffer(max_cont,max_dim,iatel_s+1,30,buffer)
+!C Clear the contacts of the atom passed to the neighboring processor
+ nn=num_cont_hb(iatel_s_nucl+1)
+!c do i=1,nn
+!c write (*,'(i2,9(3f8.3,2x))') i,(buffer(i,j+30),j=1,30)
+!c enddo
+ num_cont_hb(iatel_s_nucl)=0
+ endif
+!cd write (iout,*) 'Processor ',fg_rank,MyRank,
+!cd & ' is sending correlation contribution to processor',fg_rank-1,
+!cd & ' msglen=',msglen
+!c write (*,*) 'Processor ',fg_rank,MyRank,
+!c & ' is sending correlation contribution to processor',fg_rank-1,
+!c & ' msglen=',msglen,' CorrelType=',CorrelType
+ time00=MPI_Wtime()
+ call MPI_Send(buffer,msglen,MPI_DOUBLE_PRECISION,fg_rank-1, &
+ CorrelType,FG_COMM,IERROR)
+ time_sendrecv=time_sendrecv+MPI_Wtime()-time00
+!cd write (iout,*) 'Processor ',fg_rank,
+!cd & ' has sent correlation contribution to processor',fg_rank-1,
+!cd & ' msglen=',msglen,' CorrelID=',CorrelID
+!c write (*,*) 'Processor ',fg_rank,
+!c & ' has sent correlation contribution to processor',fg_rank-1,
+!c & ' msglen=',msglen,' CorrelID=',CorrelID
+!c msglen=msglen1
+ endif ! (fg_rank.gt.0)
+ if (ldone) goto 30
+ ldone=.true.
+ 20 continue
+!c write (*,*) 'Receiving: MyRank',MyRank,' mm',mm,' ldone',ldone
+ if (fg_rank.lt.nfgtasks-1) then
+!C Receive correlation contributions from the next processor
+ msglen=msglen1
+ if (ielend_nucl(iatel_e_nucl).lt.nct_molec(2)-1) msglen=msglen2
+!cd write (iout,*) 'Processor',fg_rank,
+!cd & ' is receiving correlation contribution from processor',fg_rank+1,
+!cd & ' msglen=',msglen,' CorrelType=',CorrelType
+!c write (*,*) 'Processor',fg_rank,
+!c &' is receiving correlation contribution from processor',fg_rank+1,
+!c & ' msglen=',msglen,' CorrelType=',CorrelType
+ time00=MPI_Wtime()
+ nbytes=-1
+ do while (nbytes.le.0)
+ call MPI_Probe(fg_rank+1,CorrelType,FG_COMM,status,IERROR)
+ call MPI_Get_count(status,MPI_DOUBLE_PRECISION,nbytes,IERROR)
+ enddo
+!c print *,'Processor',myrank,' msglen',msglen,' nbytes',nbytes
+ call MPI_Recv(buffer,nbytes,MPI_DOUBLE_PRECISION, &
+ fg_rank+1,CorrelType,FG_COMM,status,IERROR)
+ time_sendrecv=time_sendrecv+MPI_Wtime()-time00
+!c write (*,*) 'Processor',fg_rank,
+!c &' has received correlation contribution from processor',fg_rank+1,
+!c & ' msglen=',msglen,' nbytes=',nbytes
+!c write (*,*) 'The received BUFFER array:'
+!c do i=1,max_cont
+!c write (*,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,60)
+!c enddo
+ if (msglen.eq.msglen1) then
+ call unpack_buffer(max_cont,max_dim,iatel_e_nucl+1,0,buffer)
+ else if (msglen.eq.msglen2) then
+ call unpack_buffer(max_cont,max_dim,iatel_e_nucl,0,buffer)
+ call unpack_buffer(max_cont,max_dim,iatel_e_nucl+1,30,buffer)
+ else
+ write (iout,*) &
+ 'ERROR!!!! message length changed while processing correlations.'
+ write (*,*) &
+ 'ERROR!!!! message length changed while processing correlations.'
+ call MPI_Abort(MPI_COMM_WORLD,Error,IERROR)
+ endif ! msglen.eq.msglen1
+ endif ! fg_rank.lt.nfgtasks-1
+ if (ldone) goto 30
+ ldone=.true.
+ goto 10
+ 30 continue
+#endif
+ if (lprn) then
+ write (iout,'(a)') 'Contact function values:'
+ do i=nnt_molec(2),nct_molec(2)-1
+ write (iout,'(2i3,50(1x,i2,f5.2))') &
+ i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), &
+ j=1,num_cont_hb(i))
+ enddo
+ endif
+ ecorr=0.0D0
+ ecorr3=0.0d0
+!C Remove the loop below after debugging !!!
+ do i=nnt_molec(2),nct_molec(2)
+ do j=1,3
+ gradcorr_nucl(j,i)=0.0D0
+ gradxorr_nucl(j,i)=0.0D0
+ gradcorr3_nucl(j,i)=0.0D0
+ gradxorr3_nucl(j,i)=0.0D0
+ enddo
+ enddo
+ print *,"iatsc_s_nucl,iatsc_e_nucl",iatsc_s_nucl,iatsc_e_nucl
+!C Calculate the local-electrostatic correlation terms
+ do i=iatsc_s_nucl,iatsc_e_nucl
+ i1=i+1
+ num_conti=num_cont_hb(i)
+ num_conti1=num_cont_hb(i+1)
+ print *,i,num_conti,num_conti1
+ do jj=1,num_conti
+ j=jcont_hb(jj,i)
+ do kk=1,num_conti1
+ j1=jcont_hb(kk,i1)
+!c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
+!c & ' jj=',jj,' kk=',kk
+ if (j1.eq.j+1 .or. j1.eq.j-1) then
+!C
+!C Contacts I-J and (I+1)-(J+1) or (I+1)-(J-1) occur simultaneously.
+!C The system gains extra energy.
+!C Tentative expression & coefficients; assumed d(stacking)=4.5 A,
+!C parallel dipoles of stacknig bases and sin(mui)sin(muj)/eps/d^3=0.7
+!C Need to implement full formulas 34 and 35 from Liwo et al., 1998.
+!C
+ ecorr=ecorr+ehbcorr_nucl(i,j,i+1,j1,jj,kk,0.528D0,0.132D0)
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
+ 'ecorrh',i,j,ehbcorr_nucl(i,j,i+1,j1,jj,kk,0.528D0,0.132D0)
+ n_corr=n_corr+1
+ else if (j1.eq.j) then
+!C
+!C Contacts I-J and I-(J+1) occur simultaneously.
+!C The system loses extra energy.
+!C Tentative expression & c?oefficients; assumed d(stacking)=4.5 A,
+!C parallel dipoles of stacknig bases and sin(mui)sin(muj)/eps/d^3=0.7
+!C Need to implement full formulas 32 from Liwo et al., 1998.
+!C
+!c write (iout,*) 'ecorr3: i=',i,' j=',j,' i1=',i1,' j1=',j1,
+!c & ' jj=',jj,' kk=',kk
+ ecorr3=ecorr3+ehbcorr3_nucl(i,j,i+1,j,jj,kk,0.310D0,-0.155D0)
+ endif
+ enddo ! kk
+ do kk=1,num_conti
+ j1=jcont_hb(kk,i)
+!c write (iout,*) 'ecorr3: i=',i,' j=',j,' i1=',i1,' j1=',j1,
+!c & ' jj=',jj,' kk=',kk
+ if (j1.eq.j+1) then
+!C Contacts I-J and (I+1)-J occur simultaneously.
+!C The system loses extra energy.
+ ecorr3=ecorr3+ehbcorr3_nucl(i,j,i,j+1,jj,kk,0.310D0,-0.155D0)
+ endif ! j1==j+1
+ enddo ! kk
+ enddo ! jj
+ enddo ! i
+ return
+ end subroutine multibody_hb_nucl
+!-----------------------------------------------------------
+ real(kind=8) function ehbcorr_nucl(i,j,k,l,jj,kk,coeffp,coeffm)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.DERIV'
+! include 'COMMON.INTERACT'
+! include 'COMMON.CONTACTS'
+ real(kind=8),dimension(3) :: gx,gx1
+ logical :: lprn
+!el local variables
+ integer :: i,j,k,l,jj,kk,ll,ilist,m, iresshield
+ real(kind=8) :: coeffp,coeffm,eij,ekl,ees0pij,ees0pkl,ees0mij,&
+ ees0mkl,ees,coeffpees0pij,coeffmees0mij,&
+ coeffpees0pkl,coeffmees0mkl,gradlongij,gradlongkl, &
+ rlocshield
+
+ lprn=.false.
+ eij=facont_hb(jj,i)
+ ekl=facont_hb(kk,k)
+ ees0pij=ees0p(jj,i)
+ ees0pkl=ees0p(kk,k)
+ ees0mij=ees0m(jj,i)
+ ees0mkl=ees0m(kk,k)
+ ekont=eij*ekl
+ ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl)
+! print *,"ehbcorr_nucl",ekont,ees
+!cd ees=-(coeffp*ees0pkl+coeffm*ees0mkl)
+!C Following 4 lines for diagnostics.
+!cd ees0pkl=0.0D0
+!cd ees0pij=1.0D0
+!cd ees0mkl=0.0D0
+!cd ees0mij=1.0D0
+!cd write (iout,*)'Contacts have occurred for nucleic bases',
+!cd & i,j,' fcont:',eij,' eij',' eesij',ees0pij,ees0mij,' and ',k,l
+!cd & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' ees=',ees
+!C Calculate the multi-body contribution to energy.
+! ecorr_nucl=ecorr_nucl+ekont*ees
+!C Calculate multi-body contributions to the gradient.
+ coeffpees0pij=coeffp*ees0pij
+ coeffmees0mij=coeffm*ees0mij
+ coeffpees0pkl=coeffp*ees0pkl
+ coeffmees0mkl=coeffm*ees0mkl
+ do ll=1,3
+ gradxorr_nucl(ll,i)=gradxorr_nucl(ll,i) &
+ -ekont*(coeffpees0pkl*gacontp_hb1(ll,jj,i)+&
+ coeffmees0mkl*gacontm_hb1(ll,jj,i))
+ gradxorr_nucl(ll,j)=gradxorr_nucl(ll,j) &
+ -ekont*(coeffpees0pkl*gacontp_hb2(ll,jj,i)+&
+ coeffmees0mkl*gacontm_hb2(ll,jj,i))
+ gradxorr_nucl(ll,k)=gradxorr_nucl(ll,k) &
+ -ekont*(coeffpees0pij*gacontp_hb1(ll,kk,k)+&
+ coeffmees0mij*gacontm_hb1(ll,kk,k))
+ gradxorr_nucl(ll,l)=gradxorr_nucl(ll,l) &
+ -ekont*(coeffpees0pij*gacontp_hb2(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb2(ll,kk,k))
+ gradlongij=ees*ekl*gacont_hbr(ll,jj,i)- &
+ ekont*(coeffpees0pkl*gacontp_hb3(ll,jj,i)+ &
+ coeffmees0mkl*gacontm_hb3(ll,jj,i))
+ gradcorr_nucl(ll,j)=gradcorr_nucl(ll,j)+gradlongij
+ gradcorr_nucl(ll,i)=gradcorr_nucl(ll,i)-gradlongij
+ gradlongkl=ees*eij*gacont_hbr(ll,kk,k)- &
+ ekont*(coeffpees0pij*gacontp_hb3(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb3(ll,kk,k))
+ gradcorr_nucl(ll,l)=gradcorr_nucl(ll,l)+gradlongkl
+ gradcorr_nucl(ll,k)=gradcorr_nucl(ll,k)-gradlongkl
+ gradxorr_nucl(ll,i)=gradxorr_nucl(ll,i)-gradlongij
+ gradxorr_nucl(ll,j)=gradxorr_nucl(ll,j)+gradlongij
+ gradxorr_nucl(ll,k)=gradxorr_nucl(ll,k)-gradlongkl
+ gradxorr_nucl(ll,l)=gradxorr_nucl(ll,l)+gradlongkl
+ enddo
+ ehbcorr_nucl=ekont*ees
+ return
+ end function ehbcorr_nucl
+!-------------------------------------------------------------------------
+
+ real(kind=8) function ehbcorr3_nucl(i,j,k,l,jj,kk,coeffp,coeffm)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.DERIV'
+! include 'COMMON.INTERACT'
+! include 'COMMON.CONTACTS'
+ real(kind=8),dimension(3) :: gx,gx1
+ logical :: lprn
+!el local variables
+ integer :: i,j,k,l,jj,kk,ll,ilist,m, iresshield
+ real(kind=8) :: coeffp,coeffm,eij,ekl,ees0pij,ees0pkl,ees0mij,&
+ ees0mkl,ees,coeffpees0pij,coeffmees0mij,&
+ coeffpees0pkl,coeffmees0mkl,gradlongij,gradlongkl, &
+ rlocshield
+
+ lprn=.false.
+ eij=facont_hb(jj,i)
+ ekl=facont_hb(kk,k)
+ ees0pij=ees0p(jj,i)
+ ees0pkl=ees0p(kk,k)
+ ees0mij=ees0m(jj,i)
+ ees0mkl=ees0m(kk,k)
+ ekont=eij*ekl
+ ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl)
+!cd ees=-(coeffp*ees0pkl+coeffm*ees0mkl)
+!C Following 4 lines for diagnostics.
+!cd ees0pkl=0.0D0
+!cd ees0pij=1.0D0
+!cd ees0mkl=0.0D0
+!cd ees0mij=1.0D0
+!cd write (iout,*)'Contacts have occurred for nucleic bases',
+!cd & i,j,' fcont:',eij,' eij',' eesij',ees0pij,ees0mij,' and ',k,l
+!cd & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' ees=',ees
+!C Calculate the multi-body contribution to energy.
+! ecorr=ecorr+ekont*ees
+!C Calculate multi-body contributions to the gradient.
+ coeffpees0pij=coeffp*ees0pij
+ coeffmees0mij=coeffm*ees0mij
+ coeffpees0pkl=coeffp*ees0pkl
+ coeffmees0mkl=coeffm*ees0mkl
+ do ll=1,3
+ gradxorr3_nucl(ll,i)=gradxorr3_nucl(ll,i) &
+ -ekont*(coeffpees0pkl*gacontp_hb1(ll,jj,i)+&
+ coeffmees0mkl*gacontm_hb1(ll,jj,i))
+ gradxorr3_nucl(ll,j)=gradxorr3_nucl(ll,j) &
+ -ekont*(coeffpees0pkl*gacontp_hb2(ll,jj,i)+ &
+ coeffmees0mkl*gacontm_hb2(ll,jj,i))
+ gradxorr3_nucl(ll,k)=gradxorr3_nucl(ll,k) &
+ -ekont*(coeffpees0pij*gacontp_hb1(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb1(ll,kk,k))
+ gradxorr3_nucl(ll,l)=gradxorr3_nucl(ll,l) &
+ -ekont*(coeffpees0pij*gacontp_hb2(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb2(ll,kk,k))
+ gradlongij=ees*ekl*gacont_hbr(ll,jj,i)- &
+ ekont*(coeffpees0pkl*gacontp_hb3(ll,jj,i)+ &
+ coeffmees0mkl*gacontm_hb3(ll,jj,i))
+ gradcorr3_nucl(ll,j)=gradcorr3_nucl(ll,j)+gradlongij
+ gradcorr3_nucl(ll,i)=gradcorr3_nucl(ll,i)-gradlongij
+ gradlongkl=ees*eij*gacont_hbr(ll,kk,k)- &
+ ekont*(coeffpees0pij*gacontp_hb3(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb3(ll,kk,k))
+ gradcorr3_nucl(ll,l)=gradcorr3_nucl(ll,l)+gradlongkl
+ gradcorr3_nucl(ll,k)=gradcorr3_nucl(ll,k)-gradlongkl
+ gradxorr3_nucl(ll,i)=gradxorr3_nucl(ll,i)-gradlongij
+ gradxorr3_nucl(ll,j)=gradxorr3_nucl(ll,j)+gradlongij
+ gradxorr3_nucl(ll,k)=gradxorr3_nucl(ll,k)-gradlongkl
+ gradxorr3_nucl(ll,l)=gradxorr3_nucl(ll,l)+gradlongkl
+ enddo
+ ehbcorr3_nucl=ekont*ees
+ return
+ end function ehbcorr3_nucl
+#ifdef MPI
+ subroutine pack_buffer(dimen1,dimen2,atom,indx,buffer)
+ integer dimen1,dimen2,atom,indx,numcont,i,ii,k,j,num_kont,num_kont_old
+ real(kind=8):: buffer(dimen1,dimen2)
+ num_kont=num_cont_hb(atom)
+ do i=1,num_kont
+ do k=1,8
+ do j=1,3
+ buffer(i,indx+(k-1)*3+j)=zapas2(j,i,atom,k)
+ enddo ! j
+ enddo ! k
+ buffer(i,indx+25)=facont_hb(i,atom)
+ buffer(i,indx+26)=ees0p(i,atom)
+ buffer(i,indx+27)=ees0m(i,atom)
+ buffer(i,indx+28)=d_cont(i,atom)
+ buffer(i,indx+29)=dfloat(jcont_hb(i,atom))
+ enddo ! i
+ buffer(1,indx+30)=dfloat(num_kont)
+ return
+ end subroutine pack_buffer
+!c------------------------------------------------------------------------------
+ subroutine unpack_buffer(dimen1,dimen2,atom,indx,buffer)
+ integer dimen1,dimen2,atom,indx,numcont,i,ii,k,j,num_kont,num_kont_old
+ real(kind=8):: buffer(dimen1,dimen2)
+! double precision zapas
+! common /contacts_hb/ zapas(3,maxconts,maxres,8),
+! & facont_hb(maxconts,maxres),ees0p(maxconts,maxres),
+! & ees0m(maxconts,maxres),d_cont(maxconts,maxres),
+! & num_cont_hb(maxres),jcont_hb(maxconts,maxres)
+ num_kont=buffer(1,indx+30)
+ num_kont_old=num_cont_hb(atom)
+ num_cont_hb(atom)=num_kont+num_kont_old
+ do i=1,num_kont
+ ii=i+num_kont_old
+ do k=1,8
+ do j=1,3
+ zapas2(j,ii,atom,k)=buffer(i,indx+(k-1)*3+j)
+ enddo ! j
+ enddo ! k
+ facont_hb(ii,atom)=buffer(i,indx+25)
+ ees0p(ii,atom)=buffer(i,indx+26)
+ ees0m(ii,atom)=buffer(i,indx+27)
+ d_cont(i,atom)=buffer(i,indx+28)
+ jcont_hb(ii,atom)=buffer(i,indx+29)
+ enddo ! i
+ return
+ end subroutine unpack_buffer
+!c------------------------------------------------------------------------------
+#endif
+
+!----------------------------------------------------------------------------