+
+ subroutine epep_sc_base(epepbase)
+ use calc_data
+ logical :: lprn
+!el local variables
+ integer :: iint,itypi,itypi1,itypj,subchap
+ real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
+ real(kind=8) :: evdw,sig0ij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,aa,bb,ssgradlipi,ssgradlipj, &
+ sslipi,sslipj,faclip
+ integer :: ii
+ real(kind=8) :: fracinbuf
+ real (kind=8) :: epepbase
+ real (kind=8),dimension(4):: ener
+ real(kind=8) :: b1,b2,b3,b4,egb,eps_in,eps_inout_fac,eps_out
+ real(kind=8) :: ECL,Elj,Equad,Epol,eheadtail,rhead,dGCLOM2,&
+ sqom1,sqom2,sqom12,c1,c2,c3,pom,Lambf,sparrow,&
+ Chif,ChiLambf,bat,eagle,top,bot,botsq,Fcav,dtop,dFdR,dFdOM1,&
+ dFdOM2,w1,w2,w3,dGCLdR,dFdL,dFdOM12,dbot ,&
+ r1,eps_head,alphapol1,pis,facd2,d2,facd1,d1,erdxj,erdxi,federmaus,&
+ dPOLdR1,dFGBdOM2,dFGBdR1,dPOLdFGB1,RR1,MomoFac1,hawk,d1i,d1j,&
+ sig1,sig2,chis12,chis2,ee1,fgb1,a12sq,chis1
+ real(kind=8),dimension(3,2)::chead,erhead_tail
+ real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead
+ integer troll
+ eps_out=80.0d0
+ epepbase=0.0d0
+! do i=1,nres_molec(1)-1
+ do i=ibond_start,ibond_end
+ if (itype(i,1).eq.ntyp1_molec(1).or.itype(i+1,1).eq.ntyp1_molec(1)) cycle
+!C itypi = itype(i,1)
+ dxi = dc_norm(1,i)
+ dyi = dc_norm(2,i)
+ dzi = dc_norm(3,i)
+! print *,dxi,(-c(1,i)+c(1,i+1))*vbld_inv(i+1)
+ dsci_inv = vbld_inv(i+1)/2.0
+ xi=(c(1,i)+c(1,i+1))/2.0
+ yi=(c(2,i)+c(2,i+1))/2.0
+ zi=(c(3,i)+c(3,i+1))/2.0
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+ do j=nres_molec(1)+1,nres_molec(2)+nres_molec(1)
+ itypj= itype(j,2)
+ if (itype(j,2).eq.ntyp1_molec(2))cycle
+ xj=c(1,j+nres)
+ yj=c(2,j+nres)
+ zj=c(3,j+nres)
+ 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 )
+! d1i = dhead_scbasei(itypi) !this is shift of dipole/charge
+! d1j = dhead_scbasej(itypi) !this is shift of dipole/charge
+
+! Gay-berne var's
+ sig0ij = sigma_pepbase(itypj )
+ chi1 = chi_pepbase(itypj,1 )
+ chi2 = chi_pepbase(itypj,2 )
+! chi1=0.0d0
+! chi2=0.0d0
+ chi12 = chi1 * chi2
+ chip1 = chipp_pepbase(itypj,1 )
+ chip2 = chipp_pepbase(itypj,2 )
+! chip1=0.0d0
+! chip2=0.0d0
+ chip12 = chip1 * chip2
+ chis1 = chis_pepbase(itypj,1)
+ chis2 = chis_pepbase(itypj,2)
+ chis12 = chis1 * chis2
+ sig1 = sigmap1_pepbase(itypj)
+ sig2 = sigmap2_pepbase(itypj)
+! write (*,*) "sig1 = ", sig1
+! write (*,*) "sig2 = ", sig2
+ DO k = 1,3
+! location of polar head is computed by taking hydrophobic centre
+! and moving by a d1 * dc_norm vector
+! see unres publications for very informative images
+ chead(k,1) = (c(k,i)+c(k,i+1))/2.0
+! + d1i * dc_norm(k, i+nres)
+ chead(k,2) = c(k, j+nres)
+! + d1j * dc_norm(k, j+nres)
+! distance
+! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ Rhead_distance(k) = chead(k,2) - chead(k,1)
+! print *,gvdwc_pepbase(k,i)
+
+ END DO
+ Rhead = dsqrt( &
+ (Rhead_distance(1)*Rhead_distance(1)) &
+ + (Rhead_distance(2)*Rhead_distance(2)) &
+ + (Rhead_distance(3)*Rhead_distance(3)))
+
+! alpha factors from Fcav/Gcav
+ b1 = alphasur_pepbase(1,itypj)
+! b1=0.0d0
+ b2 = alphasur_pepbase(2,itypj)
+ b3 = alphasur_pepbase(3,itypj)
+ b4 = alphasur_pepbase(4,itypj)
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+ rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj)
+! print *,i,j,rrij
+ rij = dsqrt(rrij)
+!----------------------------
+ evdwij = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ Fcav=0.0d0
+ eheadtail = 0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ Fcav = 0.0d0
+ dFdR = 0.0d0
+ dCAVdOM1 = 0.0d0
+ dCAVdOM2 = 0.0d0
+ dCAVdOM12 = 0.0d0
+ dscj_inv = vbld_inv(j+nres)
+ CALL sc_angular
+! this should be in elgrad_init but om's are calculated by sc_angular
+! which in turn is used by older potentials
+! om = omega, sqom = om^2
+ sqom1 = om1 * om1
+ sqom2 = om2 * om2
+ sqom12 = om12 * om12
+
+! now we calculate EGB - Gey-Berne
+! It will be summed up in evdwij and saved in evdw
+ sigsq = 1.0D0 / sigsq
+ sig = sig0ij * dsqrt(sigsq)
+ rij_shift = 1.0/rij - sig + sig0ij
+ IF (rij_shift.le.0.0D0) THEN
+ evdw = 1.0D20
+ RETURN
+ END IF
+ sigder = -sig * sigsq
+ rij_shift = 1.0D0 / rij_shift
+ fac = rij_shift**expon
+ c1 = fac * fac * aa_pepbase(itypj)
+! c1 = 0.0d0
+ c2 = fac * bb_pepbase(itypj)
+! c2 = 0.0d0
+ evdwij = eps1 * eps2rt * eps3rt * ( c1 + c2 )
+ eps2der = eps3rt * evdwij
+ eps3der = eps2rt * evdwij
+! evdwij = 4.0d0 * eps2rt * eps3rt * evdwij
+ evdwij = eps2rt * eps3rt * evdwij
+ c1 = c1 * eps1 * eps2rt**2 * eps3rt**2
+ fac = -expon * (c1 + evdwij) * rij_shift
+ sigder = fac * sigder
+! fac = rij * fac
+! Calculate distance derivative
+ gg(1) = fac
+ gg(2) = fac
+ gg(3) = fac
+ fac = chis1 * sqom1 + chis2 * sqom2 &
+ - 2.0d0 * chis12 * om1 * om2 * om12
+! we will use pom later in Gcav, so dont mess with it!
+ pom = 1.0d0 - chis1 * chis2 * sqom12
+ Lambf = (1.0d0 - (fac / pom))
+ Lambf = dsqrt(Lambf)
+ sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0)
+! write (*,*) "sparrow = ", sparrow
+ Chif = 1.0d0/rij * sparrow
+ ChiLambf = Chif * Lambf
+ eagle = dsqrt(ChiLambf)
+ bat = ChiLambf ** 11.0d0
+ top = b1 * ( eagle + b2 * ChiLambf - b3 )
+ bot = 1.0d0 + b4 * (ChiLambf ** 12.0d0)
+ botsq = bot * bot
+ Fcav = top / bot
+! print *,i,j,Fcav
+ dtop = b1 * ((Lambf / (2.0d0 * eagle)) + (b2 * Lambf))
+ dbot = 12.0d0 * b4 * bat * Lambf
+ dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow
+! dFdR = 0.0d0
+! write (*,*) "dFcav/dR = ", dFdR
+ dtop = b1 * ((Chif / (2.0d0 * eagle)) + (b2 * Chif))
+ dbot = 12.0d0 * b4 * bat * Chif
+ eagle = Lambf * pom
+ dFdOM1 = -(chis1 * om1 - chis12 * om2 * om12) / (eagle)
+ dFdOM2 = -(chis2 * om2 - chis12 * om1 * om12) / (eagle)
+ dFdOM12 = chis12 * (chis1 * om1 * om12 - om2) &
+ * (chis2 * om2 * om12 - om1) / (eagle * pom)
+
+ dFdL = ((dtop * bot - top * dbot) / botsq)
+! dFdL = 0.0d0
+ dCAVdOM1 = dFdL * ( dFdOM1 )
+ dCAVdOM2 = dFdL * ( dFdOM2 )
+ dCAVdOM12 = dFdL * ( dFdOM12 )
+
+ ertail(1) = xj*rij
+ ertail(2) = yj*rij
+ ertail(3) = zj*rij
+ DO k = 1, 3
+! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
+! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
+ pom = ertail(k)
+!-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres))
+ gvdwc_pepbase(k,i) = gvdwc_pepbase(k,i) &
+ - (( dFdR + gg(k) ) * pom)/2.0
+! print *,gvdwc_pepbase(k,i),i,(( dFdR + gg(k) ) * pom)/2.0
+! +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
+! +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+! & - ( dFdR * pom )
+ pom = ertail(k)
+!-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx_pepbase(k,j) = gvdwx_pepbase(k,j) &
+ + (( dFdR + gg(k) ) * pom)
+! +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
+! +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+!c! & + ( dFdR * pom )
+
+ gvdwc_pepbase(k,i+1) = gvdwc_pepbase(k,i+1) &
+ - (( dFdR + gg(k) ) * ertail(k))/2.0
+! print *,gvdwc_pepbase(k,i+1),i+1,(( dFdR + gg(k) ) * pom)/2.0
+
+!c! & - ( dFdR * ertail(k))
+
+ gvdwc_pepbase(k,j) = gvdwc_pepbase(k,j) &
+ + (( dFdR + gg(k) ) * ertail(k))
+!c! & + ( dFdR * ertail(k))
+
+ gg(k) = 0.0d0
+!c! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
+!c! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
+ END DO
+
+
+ w1 = wdipdip_pepbase(1,itypj)
+ w2 = -wdipdip_pepbase(3,itypj)/2.0
+ w3 = wdipdip_pepbase(2,itypj)
+! w1=0.0d0
+! w2=0.0d0
+!c!-------------------------------------------------------------------
+!c! ECL
+! w3=0.0d0
+ fac = (om12 - 3.0d0 * om1 * om2)
+ c1 = (w1 / (Rhead**3.0d0)) * fac
+ c2 = (w2 / Rhead ** 6.0d0) &
+ * (4.0d0 + fac * fac -3.0d0 * (sqom1 + sqom2))
+ c3= (w3/ Rhead ** 6.0d0) &
+ * (2.0d0 - 2.0d0*fac*fac +3.0d0*(sqom1 + sqom2))
+
+ ECL = c1 - c2 + c3
+
+ c1 = (-3.0d0 * w1 * fac) / (Rhead ** 4.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead ** 7.0d0) &
+ * (4.0d0 + fac * fac - 3.0d0 * (sqom1 + sqom2))
+ c3= (-6.0d0 * w3) / (Rhead ** 7.0d0) &
+ * (2.0d0 - 2.0d0*fac*fac +3.0d0*(sqom1 + sqom2))
+
+ dGCLdR = c1 - c2 + c3
+!c! dECL/dom1
+ c1 = (-3.0d0 * w1 * om2 ) / (Rhead**3.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead**6.0d0) &
+ * ( om2 * om12 - 3.0d0 * om1 * sqom2 + om1 )
+ c3 =(6.0d0*w3/ Rhead ** 6.0d0)*(om1-2.0d0*(fac)*(-om2))
+ dGCLdOM1 = c1 - c2 + c3
+!c! dECL/dom2
+ c1 = (-3.0d0 * w1 * om1 ) / (Rhead**3.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead**6.0d0) &
+ * ( om1 * om12 - 3.0d0 * sqom1 * om2 + om2 )
+ c3 =(6.0d0*w3/ Rhead ** 6.0d0)*(om2-2.0d0*(fac)*(-om1))
+
+ dGCLdOM2 = c1 - c2 + c3
+!c! dECL/dom12
+ c1 = w1 / (Rhead ** 3.0d0)
+ c2 = ( 2.0d0 * w2 * fac ) / Rhead ** 6.0d0
+ c3 = (w3/ Rhead ** 6.0d0)*(-4.0d0*fac)
+ dGCLdOM12 = c1 - c2 + c3
+ DO k= 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ END DO
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j+nres) )
+! facd1 = d1 * vbld_inv(i+nres)
+! facd2 = d2 * vbld_inv(j+nres)
+ DO k = 1, 3
+
+! pom = erhead(k)
+!+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+! gvdwx_pepbase(k,i) = gvdwx_scbase(k,i) &
+! - dGCLdR * pom
+ pom = erhead(k)
+!+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx_pepbase(k,j) = gvdwx_pepbase(k,j) &
+ + dGCLdR * pom
+
+ gvdwc_pepbase(k,i) = gvdwc_pepbase(k,i) &
+ - dGCLdR * erhead(k)/2.0d0
+! print *,gvdwc_pepbase(k,i+1),i+1,- dGCLdR * erhead(k)/2.0d0
+ gvdwc_pepbase(k,i+1) = gvdwc_pepbase(k,i+1) &
+ - dGCLdR * erhead(k)/2.0d0
+! print *,gvdwc_pepbase(k,i+1),i+1,- dGCLdR * erhead(k)/2.0d0
+ gvdwc_pepbase(k,j) = gvdwc_pepbase(k,j) &
+ + dGCLdR * erhead(k)
+ END DO
+! print *,i,j,evdwij,Fcav,ECL,"vdw,cav,ecl"
+ epepbase=epepbase+evdwij+Fcav+ECL
+ call sc_grad_pepbase
+ enddo
+ enddo
+ END SUBROUTINE epep_sc_base
+ SUBROUTINE sc_grad_pepbase
+ use calc_data
+
+ real (kind=8) :: dcosom1(3),dcosom2(3)
+ eom1 = &
+ eps2der * eps2rt_om1 &
+ - 2.0D0 * alf1 * eps3der &
+ + sigder * sigsq_om1 &
+ + dCAVdOM1 &
+ + dGCLdOM1 &
+ + dPOLdOM1
+
+ eom2 = &
+ eps2der * eps2rt_om2 &
+ + 2.0D0 * alf2 * eps3der &
+ + sigder * sigsq_om2 &
+ + dCAVdOM2 &
+ + dGCLdOM2 &
+ + dPOLdOM2
+
+ eom12 = &
+ evdwij * eps1_om12 &
+ + eps2der * eps2rt_om12 &
+ - 2.0D0 * alf12 * eps3der &
+ + sigder *sigsq_om12 &
+ + dCAVdOM12 &
+ + dGCLdOM12
+! om12=0.0
+! eom12=0.0
+! print *,eom1,eom2,eom12,om12,i,j,"eom1,2,12",erij(1),erij(2),erij(3)
+! if (i.eq.30) print *,gvdwc_pepbase(k,i),- gg(k),&
+! (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i)))&
+! *dsci_inv*2.0
+! print *,dsci_inv,dscj_inv,dc_norm(2,nres+j),dc_norm(2,nres+i),&
+! gg(1),gg(2),"rozne"
+ DO k = 1, 3
+ dcosom1(k) = rij * (dc_norm(k,i) - om1 * erij(k))
+ dcosom2(k) = rij * (dc_norm(k,nres+j) - om2 * erij(k))
+ gg(k) = gg(k) + eom1 * dcosom1(k) + eom2 * dcosom2(k)
+ gvdwc_pepbase(k,i)= gvdwc_pepbase(k,i) +0.5*(- gg(k)) &
+ + (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i)))&
+ *dsci_inv*2.0 &
+ - (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0
+ gvdwc_pepbase(k,i+1)= gvdwc_pepbase(k,i+1) +0.5*(- gg(k)) &
+ - (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i))) &
+ *dsci_inv*2.0 &
+ + (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0
+! print *,eom12,eom2,om12,om2
+!eom12*(-dc_norm(k,i)/2.0-om12*dc_norm(k,nres+j)),&
+! (eom2*(erij(k)-om2*dc_norm(k,nres+j)))
+ gvdwx_pepbase(k,j)= gvdwx_pepbase(k,j) + gg(k) &
+ + (eom12*(dc_norm(k,i)-om12*dc_norm(k,nres+j))&
+ + eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ gvdwc_pepbase(k,j)=gvdwc_pepbase(k,j)+gg(k)
+ END DO
+ RETURN
+ END SUBROUTINE sc_grad_pepbase
+ subroutine eprot_sc_phosphate(escpho)
+ use calc_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.LOCAL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.NAMES'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CALC'
+! include 'COMMON.CONTROL'
+! include 'COMMON.SBRIDGE'
+ logical :: lprn
+!el local variables
+ integer :: iint,itypi,itypi1,itypj,subchap
+ real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
+ real(kind=8) :: evdw,sig0ij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,aa,bb,ssgradlipi,ssgradlipj, &
+ sslipi,sslipj,faclip,alpha_sco
+ integer :: ii
+ real(kind=8) :: fracinbuf
+ real (kind=8) :: escpho
+ real (kind=8),dimension(4):: ener
+ real(kind=8) :: b1,b2,b3,b4,egb,eps_in,eps_inout_fac,eps_out
+ real(kind=8) :: ECL,Elj,Equad,Epol,eheadtail,rhead,dGCLOM2,&
+ sqom1,sqom2,sqom12,c1,c2,pom,Lambf,sparrow,&
+ Chif,ChiLambf,bat,eagle,top,bot,botsq,Fcav,dtop,dFdR,dFdOM1,&
+ dFdOM2,w1,w2,dGCLdR,dFdL,dFdOM12,dbot ,&
+ r1,eps_head,alphapol1,pis,facd2,d2,facd1,d1,erdxj,erdxi,federmaus,&
+ dPOLdR1,dFGBdOM2,dFGBdR1,dPOLdFGB1,RR1,MomoFac1,hawk,d1i,d1j,&
+ sig1,sig2,chis12,chis2,ee1,fgb1,a12sq,chis1,Rhead_sq,Qij,dFGBdOM1
+ real(kind=8),dimension(3,2)::chead,erhead_tail
+ real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead
+ integer troll
+ eps_out=80.0d0
+ escpho=0.0d0
+! do i=1,nres_molec(1)
+ do i=ibond_start,ibond_end
+ if (itype(i,1).eq.ntyp1_molec(1)) cycle
+ itypi = itype(i,1)
+ dxi = dc_norm(1,nres+i)
+ dyi = dc_norm(2,nres+i)
+ dzi = dc_norm(3,nres+i)
+ dsci_inv = vbld_inv(i+nres)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+ do j=nres_molec(1)+1,nres_molec(2)+nres_molec(1)-1
+ itypj= itype(j,2)
+ if ((itype(j,2).eq.ntyp1_molec(2)).or.&
+ (itype(j+1,2).eq.ntyp1_molec(2))) cycle
+ xj=(c(1,j)+c(1,j+1))/2.0
+ yj=(c(2,j)+c(2,j+1))/2.0
+ zj=(c(3,j)+c(3,j+1))/2.0
+ 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,j )
+ dyj = dc_norm( 2,j )
+ dzj = dc_norm( 3,j )
+ dscj_inv = vbld_inv(j+1)
+
+! Gay-berne var's
+ sig0ij = sigma_scpho(itypi )
+ chi1 = chi_scpho(itypi,1 )
+ chi2 = chi_scpho(itypi,2 )
+! chi1=0.0d0
+! chi2=0.0d0
+ chi12 = chi1 * chi2
+ chip1 = chipp_scpho(itypi,1 )
+ chip2 = chipp_scpho(itypi,2 )
+! chip1=0.0d0
+! chip2=0.0d0
+ chip12 = chip1 * chip2
+ chis1 = chis_scpho(itypi,1)
+ chis2 = chis_scpho(itypi,2)
+ chis12 = chis1 * chis2
+ sig1 = sigmap1_scpho(itypi)
+ sig2 = sigmap2_scpho(itypi)
+! write (*,*) "sig1 = ", sig1
+! write (*,*) "sig1 = ", sig1
+! write (*,*) "sig2 = ", sig2
+! alpha factors from Fcav/Gcav
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+ a12sq = rborn_scphoi(itypi) * rborn_scphoj(itypi)
+
+ b1 = alphasur_scpho(1,itypi)
+! b1=0.0d0
+ b2 = alphasur_scpho(2,itypi)
+ b3 = alphasur_scpho(3,itypi)
+ b4 = alphasur_scpho(4,itypi)
+! used to determine whether we want to do quadrupole calculations
+! used by Fgb
+ eps_in = epsintab_scpho(itypi)
+ if (eps_in.eq.0.0) eps_in=1.0
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+! write (*,*) "eps_inout_fac = ", eps_inout_fac
+!-------------------------------------------------------------------
+! tail location and distance calculations
+ d1i = dhead_scphoi(itypi) !this is shift of dipole/charge
+ d1j = 0.0
+ DO k = 1,3
+! location of polar head is computed by taking hydrophobic centre
+! and moving by a d1 * dc_norm vector
+! see unres publications for very informative images
+ chead(k,1) = c(k, i+nres) + d1i * dc_norm(k, i+nres)
+ chead(k,2) = (c(k, j) + c(k, j+1))/2.0
+! distance
+! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ Rhead_distance(k) = chead(k,2) - chead(k,1)
+ END DO
+! pitagoras (root of sum of squares)
+ Rhead = dsqrt( &
+ (Rhead_distance(1)*Rhead_distance(1)) &
+ + (Rhead_distance(2)*Rhead_distance(2)) &
+ + (Rhead_distance(3)*Rhead_distance(3)))
+ Rhead_sq=Rhead**2.0
+!-------------------------------------------------------------------
+! zero everything that should be zero'ed
+ evdwij = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ Fcav=0.0d0
+ eheadtail = 0.0d0
+ dGCLdR=0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ Fcav = 0.0d0
+ dFdR = 0.0d0
+ dCAVdOM1 = 0.0d0
+ dCAVdOM2 = 0.0d0
+ dCAVdOM12 = 0.0d0
+ dscj_inv = vbld_inv(j+1)/2.0
+!dhead_scbasej(itypi,itypj)
+! print *,i,j,dscj_inv,dsci_inv
+! rij holds 1/(distance of Calpha atoms)
+ rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj)
+ rij = dsqrt(rrij)
+!----------------------------
+ CALL sc_angular
+! this should be in elgrad_init but om's are calculated by sc_angular
+! which in turn is used by older potentials
+! om = omega, sqom = om^2
+ sqom1 = om1 * om1
+ sqom2 = om2 * om2
+ sqom12 = om12 * om12
+
+! now we calculate EGB - Gey-Berne
+! It will be summed up in evdwij and saved in evdw
+ sigsq = 1.0D0 / sigsq
+ sig = sig0ij * dsqrt(sigsq)
+! rij_shift = 1.0D0 / rij - sig + sig0ij
+ rij_shift = 1.0/rij - sig + sig0ij
+ IF (rij_shift.le.0.0D0) THEN
+ evdw = 1.0D20
+ RETURN
+ END IF
+ sigder = -sig * sigsq
+ rij_shift = 1.0D0 / rij_shift
+ fac = rij_shift**expon
+ c1 = fac * fac * aa_scpho(itypi)
+! c1 = 0.0d0
+ c2 = fac * bb_scpho(itypi)
+! c2 = 0.0d0
+ evdwij = eps1 * eps2rt * eps3rt * ( c1 + c2 )
+ eps2der = eps3rt * evdwij
+ eps3der = eps2rt * evdwij
+! evdwij = 4.0d0 * eps2rt * eps3rt * evdwij
+ evdwij = eps2rt * eps3rt * evdwij
+ c1 = c1 * eps1 * eps2rt**2 * eps3rt**2
+ fac = -expon * (c1 + evdwij) * rij_shift
+ sigder = fac * sigder
+! fac = rij * fac
+! Calculate distance derivative
+ gg(1) = fac
+ gg(2) = fac
+ gg(3) = fac
+ fac = chis1 * sqom1 + chis2 * sqom2 &
+ - 2.0d0 * chis12 * om1 * om2 * om12
+! we will use pom later in Gcav, so dont mess with it!
+ pom = 1.0d0 - chis1 * chis2 * sqom12
+ Lambf = (1.0d0 - (fac / pom))
+ Lambf = dsqrt(Lambf)
+ sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0)
+! write (*,*) "sparrow = ", sparrow
+ Chif = 1.0d0/rij * sparrow
+ ChiLambf = Chif * Lambf
+ eagle = dsqrt(ChiLambf)
+ bat = ChiLambf ** 11.0d0
+ top = b1 * ( eagle + b2 * ChiLambf - b3 )
+ bot = 1.0d0 + b4 * (ChiLambf ** 12.0d0)
+ botsq = bot * bot
+ Fcav = top / bot
+ dtop = b1 * ((Lambf / (2.0d0 * eagle)) + (b2 * Lambf))
+ dbot = 12.0d0 * b4 * bat * Lambf
+ dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow
+! dFdR = 0.0d0
+! write (*,*) "dFcav/dR = ", dFdR
+ dtop = b1 * ((Chif / (2.0d0 * eagle)) + (b2 * Chif))
+ dbot = 12.0d0 * b4 * bat * Chif
+ eagle = Lambf * pom
+ dFdOM1 = -(chis1 * om1 - chis12 * om2 * om12) / (eagle)
+ dFdOM2 = -(chis2 * om2 - chis12 * om1 * om12) / (eagle)
+ dFdOM12 = chis12 * (chis1 * om1 * om12 - om2) &
+ * (chis2 * om2 * om12 - om1) / (eagle * pom)
+
+ dFdL = ((dtop * bot - top * dbot) / botsq)
+! dFdL = 0.0d0
+ dCAVdOM1 = dFdL * ( dFdOM1 )
+ dCAVdOM2 = dFdL * ( dFdOM2 )
+ dCAVdOM12 = dFdL * ( dFdOM12 )
+
+ ertail(1) = xj*rij
+ ertail(2) = yj*rij
+ ertail(3) = zj*rij
+ DO k = 1, 3
+! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
+! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
+! if (i.eq.3) print *,'decl0',gvdwx_scpho(k,i),i
+
+ pom = ertail(k)
+! print *,pom,gg(k),dFdR
+!-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx_scpho(k,i) = gvdwx_scpho(k,i) &
+ - (( dFdR + gg(k) ) * pom)
+! +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
+! +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+! & - ( dFdR * pom )
+! pom = ertail(k)
+!-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres))
+! gvdwx_scpho(k,j) = gvdwx_scpho(k,j) &
+! + (( dFdR + gg(k) ) * pom)
+! +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
+! +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+!c! & + ( dFdR * pom )
+
+ gvdwc_scpho(k,i) = gvdwc_scpho(k,i) &
+ - (( dFdR + gg(k) ) * ertail(k))
+!c! & - ( dFdR * ertail(k))
+
+ gvdwc_scpho(k,j) = gvdwc_scpho(k,j) &
+ + (( dFdR + gg(k) ) * ertail(k))/2.0
+
+ gvdwc_scpho(k,j+1) = gvdwc_scpho(k,j+1) &
+ + (( dFdR + gg(k) ) * ertail(k))/2.0
+
+!c! & + ( dFdR * ertail(k))
+
+ gg(k) = 0.0d0
+ ENDDO
+!c! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
+!c! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
+! alphapol1 = alphapol_scpho(itypi)
+ if (wqq_scpho(itypi).ne.0.0) then
+ Qij=wqq_scpho(itypi)/eps_in
+ alpha_sco=1.d0/alphi_scpho(itypi)
+! Qij=0.0
+ Ecl = (332.0d0 * Qij*dexp(-Rhead*alpha_sco)) / Rhead
+!c! derivative of Ecl is Gcl...
+ dGCLdR = (-332.0d0 * Qij*dexp(-Rhead*alpha_sco)* &
+ (Rhead*alpha_sco+1) ) / Rhead_sq
+ if (energy_dec) write(iout,*) "ECL",ECL,Rhead,1.0/rij
+ else if (wqdip_scpho(2,itypi).gt.0.0d0) then
+ w1 = wqdip_scpho(1,itypi)
+ w2 = wqdip_scpho(2,itypi)
+! w1=0.0d0
+! w2=0.0d0
+! pis = sig0head_scbase(itypi,itypj)
+! eps_head = epshead_scbase(itypi,itypj)
+!c!-------------------------------------------------------------------
+
+!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj)
+!c! & +dhead(1,1,itypi,itypj))**2))
+!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj)
+!c! & +dhead(2,1,itypi,itypj))**2))
+
+!c!-------------------------------------------------------------------
+!c! ecl
+ sparrow = w1 * om1
+ hawk = w2 * (1.0d0 - sqom2)
+ Ecl = sparrow / Rhead**2.0d0 &
+ - hawk / Rhead**4.0d0
+!c!-------------------------------------------------------------------
+ if (energy_dec) write(iout,*) "ECLdipdip",ECL,Rhead,&
+ 1.0/rij,sparrow
+
+!c! derivative of ecl is Gcl
+!c! dF/dr part
+ dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 &
+ + 4.0d0 * hawk / Rhead**5.0d0
+!c! dF/dom1
+ dGCLdOM1 = (w1) / (Rhead**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = (2.0d0 * w2 * om2) / (Rhead ** 4.0d0)
+ endif
+
+!c--------------------------------------------------------------------
+!c Polarization energy
+!c Epol
+ R1 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances tail is center of side-chain
+ R1=R1+((c(k,j)+c(k,j+1))/2.0-chead(k,1))**2
+ END DO
+!c! Pitagoras
+ R1 = dsqrt(R1)
+
+ alphapol1 = alphapol_scpho(itypi)
+! alphapol1=0.0
+ MomoFac1 = (1.0d0 - chi2 * sqom1)
+ RR1 = R1 * R1 / MomoFac1
+ ee1 = exp(-( RR1 / (4.0d0 * a12sq) ))
+! print *,"ee1",ee1,a12sq,alphapol1,eps_inout_fac
+ fgb1 = sqrt( RR1 + a12sq * ee1)
+! eps_inout_fac=0.0d0
+ epol = 332.0d0 * eps_inout_fac * (( alphapol1 / fgb1 )**4.0d0)
+! derivative of Epol is Gpol...
+ dPOLdFGB1 = -(1328.0d0 * eps_inout_fac * alphapol1 ** 4.0d0) &
+ / (fgb1 ** 5.0d0)
+ dFGBdR1 = ( (R1 / MomoFac1) &
+ * ( 2.0d0 - (0.5d0 * ee1) ) ) &
+ / ( 2.0d0 * fgb1 )
+ dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1)) &
+ * (2.0d0 - 0.5d0 * ee1) ) &
+ / (2.0d0 * fgb1)
+ dPOLdR1 = dPOLdFGB1 * dFGBdR1
+! dPOLdR1 = 0.0d0
+! dPOLdOM1 = 0.0d0
+ dFGBdOM1 = (((R1 * R1 * chi2 * om1) / (MomoFac1 * MomoFac1)) &
+ * (2.0d0 - 0.5d0 * ee1) ) &
+ / (2.0d0 * fgb1)
+
+ dPOLdOM1 = dPOLdFGB1 * dFGBdOM1
+ dPOLdOM2 = 0.0
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,1) = (((c(k,j)+c(k,j+1))/2.0-chead(k,1))/R1)
+ END DO
+
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j) )
+ bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) )
+! bat=0.0d0
+ federmaus = scalar(erhead_tail(1,1),dC_norm(1,j))
+ facd1 = d1i * vbld_inv(i+nres)
+ facd2 = d1j * vbld_inv(j)
+! facd4 = dtail(2,itypi,itypj) * vbld_inv(j+nres)
+
+ DO k = 1, 3
+ hawk = (erhead_tail(k,1) + &
+ facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres)))
+! facd1=0.0d0
+! facd2=0.0d0
+! if (i.eq.3) print *,'decl1',dGCLdR,dPOLdR1,gvdwc_scpho(k,i),i,&
+! pom,(erhead_tail(k,1))
+
+! print *,'decl',dGCLdR,dPOLdR1,gvdwc_scpho(k,i)
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx_scpho(k,i) = gvdwx_scpho(k,i) &
+ - dGCLdR * pom &
+ - dPOLdR1 * (erhead_tail(k,1))
+! & - dGLJdR * pom
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
+! gvdwx_scpho(k,j) = gvdwx_scpho(k,j) &
+! + dGCLdR * pom &
+! + dPOLdR1 * (erhead_tail(k,1))
+! & + dGLJdR * pom
+
+
+ gvdwc_scpho(k,i) = gvdwc_scpho(k,i) &
+ - dGCLdR * erhead(k) &
+ - dPOLdR1 * erhead_tail(k,1)
+! & - dGLJdR * erhead(k)
+
+ gvdwc_scpho(k,j) = gvdwc_scpho(k,j) &
+ + (dGCLdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1))/2.0
+ gvdwc_scpho(k,j+1) = gvdwc_scpho(k,j+1) &
+ + (dGCLdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1))/2.0
+
+! & + dGLJdR * erhead(k)
+! if (i.eq.3) print *,'decl2',dGCLdR,dPOLdR1,gvdwc_scpho(k,i),i
+
+ END DO
+! if (i.eq.3) print *,i,j,evdwij,epol,Fcav,ECL
+ if (energy_dec) write (iout,'(a22,2i5,4f8.3,f16.3)'), &
+ "escpho:evdw,pol,cav,CL",i,j,evdwij,epol,Fcav,ECL,escpho
+ escpho=escpho+evdwij+epol+Fcav+ECL
+ call sc_grad_scpho
+ enddo
+
+ enddo
+
+ return
+ end subroutine eprot_sc_phosphate
+ SUBROUTINE sc_grad_scpho
+ use calc_data
+
+ real (kind=8) :: dcosom1(3),dcosom2(3)
+ eom1 = &
+ eps2der * eps2rt_om1 &
+ - 2.0D0 * alf1 * eps3der &
+ + sigder * sigsq_om1 &
+ + dCAVdOM1 &
+ + dGCLdOM1 &
+ + dPOLdOM1
+
+ eom2 = &
+ eps2der * eps2rt_om2 &
+ + 2.0D0 * alf2 * eps3der &
+ + sigder * sigsq_om2 &
+ + dCAVdOM2 &
+ + dGCLdOM2 &
+ + dPOLdOM2
+
+ eom12 = &
+ evdwij * eps1_om12 &
+ + eps2der * eps2rt_om12 &
+ - 2.0D0 * alf12 * eps3der &
+ + sigder *sigsq_om12 &
+ + dCAVdOM12 &
+ + dGCLdOM12
+! om12=0.0
+! eom12=0.0
+! print *,eom1,eom2,eom12,om12,i,j,"eom1,2,12",erij(1),erij(2),erij(3)
+! if (i.eq.30) print *,gvdwc_scpho(k,i),- gg(k),&
+! (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i)))&
+! *dsci_inv*2.0
+! print *,dsci_inv,dscj_inv,dc_norm(2,nres+j),dc_norm(2,nres+i),&
+! gg(1),gg(2),"rozne"
+ DO k = 1, 3
+ dcosom1(k) = rij * (dc_norm(k,nres+i) - om1 * erij(k))
+ dcosom2(k) = rij * (dc_norm(k,j) - om2 * erij(k))
+ gg(k) = gg(k) + eom1 * dcosom1(k) + eom2 * dcosom2(k)
+ gvdwc_scpho(k,j)= gvdwc_scpho(k,j) +0.5*( gg(k)) &
+ + (-eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,j)))&
+ *dscj_inv*2.0 &
+ - (eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv*2.0
+ gvdwc_scpho(k,j+1)= gvdwc_scpho(k,j+1) +0.5*( gg(k)) &
+ - (-eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,j))) &
+ *dscj_inv*2.0 &
+ + (eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv*2.0
+ gvdwx_scpho(k,i)= gvdwx_scpho(k,i) - gg(k) &
+ + (eom12*(dc_norm(k,j)-om12*dc_norm(k,nres+i)) &
+ + eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+
+! print *,eom12,eom2,om12,om2
+!eom12*(-dc_norm(k,i)/2.0-om12*dc_norm(k,nres+j)),&
+! (eom2*(erij(k)-om2*dc_norm(k,nres+j)))
+! gvdwx_scpho(k,j)= gvdwx_scpho(k,j) + gg(k) &
+! + (eom12*(dc_norm(k,i)-om12*dc_norm(k,nres+j))&
+! + eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ gvdwc_scpho(k,i)=gvdwc_scpho(k,i)-gg(k)
+ END DO
+ RETURN
+ END SUBROUTINE sc_grad_scpho
+ subroutine eprot_pep_phosphate(epeppho)
+ use calc_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.LOCAL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.NAMES'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CALC'
+! include 'COMMON.CONTROL'
+! include 'COMMON.SBRIDGE'
+ logical :: lprn
+!el local variables
+ integer :: iint,itypi,itypi1,itypj,subchap
+ real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
+ real(kind=8) :: evdw,sig0ij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,aa,bb,ssgradlipi,ssgradlipj, &
+ sslipi,sslipj,faclip
+ integer :: ii
+ real(kind=8) :: fracinbuf
+ real (kind=8) :: epeppho
+ real (kind=8),dimension(4):: ener
+ real(kind=8) :: b1,b2,b3,b4,egb,eps_in,eps_inout_fac,eps_out
+ real(kind=8) :: ECL,Elj,Equad,Epol,eheadtail,rhead,dGCLOM2,&
+ sqom1,sqom2,sqom12,c1,c2,pom,Lambf,sparrow,&
+ Chif,ChiLambf,bat,eagle,top,bot,botsq,Fcav,dtop,dFdR,dFdOM1,&
+ dFdOM2,w1,w2,dGCLdR,dFdL,dFdOM12,dbot ,&
+ r1,eps_head,alphapol1,pis,facd2,d2,facd1,d1,erdxj,erdxi,federmaus,&
+ dPOLdR1,dFGBdOM2,dFGBdR1,dPOLdFGB1,RR1,MomoFac1,hawk,d1i,d1j,&
+ sig1,sig2,chis12,chis2,ee1,fgb1,a12sq,chis1,Rhead_sq,Qij,dFGBdOM1
+ real(kind=8),dimension(3,2)::chead,erhead_tail
+ real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead
+ integer troll
+ real (kind=8) :: dcosom1(3),dcosom2(3)
+ epeppho=0.0d0
+! do i=1,nres_molec(1)
+ do i=ibond_start,ibond_end
+ if (itype(i,1).eq.ntyp1_molec(1)) cycle
+ itypi = itype(i,1)
+ dsci_inv = vbld_inv(i+1)/2.0
+ dxi = dc_norm(1,i)
+ dyi = dc_norm(2,i)
+ dzi = dc_norm(3,i)
+ xi=(c(1,i)+c(1,i+1))/2.0
+ yi=(c(2,i)+c(2,i+1))/2.0
+ zi=(c(3,i)+c(3,i+1))/2.0
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+ do j=nres_molec(1)+1,nres_molec(2)+nres_molec(1)-1
+ itypj= itype(j,2)
+ if ((itype(j,2).eq.ntyp1_molec(2)).or.&
+ (itype(j+1,2).eq.ntyp1_molec(2))) cycle
+ xj=(c(1,j)+c(1,j+1))/2.0
+ yj=(c(2,j)+c(2,j+1))/2.0
+ zj=(c(3,j)+c(3,j+1))/2.0
+ 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
+ rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj)
+ rij = dsqrt(rrij)
+ dxj = dc_norm( 1,j )
+ dyj = dc_norm( 2,j )
+ dzj = dc_norm( 3,j )
+ dscj_inv = vbld_inv(j+1)/2.0
+! Gay-berne var's
+ sig0ij = sigma_peppho
+ chi1=0.0d0
+ chi2=0.0d0
+ chi12 = chi1 * chi2
+ chip1=0.0d0
+ chip2=0.0d0
+ chip12 = chip1 * chip2
+ chis1 = 0.0d0
+ chis2 = 0.0d0
+ chis12 = chis1 * chis2
+ sig1 = sigmap1_peppho
+ sig2 = sigmap2_peppho
+! write (*,*) "sig1 = ", sig1
+! write (*,*) "sig1 = ", sig1
+! write (*,*) "sig2 = ", sig2
+! alpha factors from Fcav/Gcav
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+ b1 = alphasur_peppho(1)
+! b1=0.0d0
+ b2 = alphasur_peppho(2)
+ b3 = alphasur_peppho(3)
+ b4 = alphasur_peppho(4)
+ CALL sc_angular
+ sqom1=om1*om1
+ evdwij = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ Fcav=0.0d0
+ eheadtail = 0.0d0
+ dGCLdR=0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ Fcav = 0.0d0
+ dFdR = 0.0d0
+ dCAVdOM1 = 0.0d0
+ dCAVdOM2 = 0.0d0
+ dCAVdOM12 = 0.0d0
+ rij_shift = rij
+ fac = rij_shift**expon
+ c1 = fac * fac * aa_peppho
+! c1 = 0.0d0
+ c2 = fac * bb_peppho
+! c2 = 0.0d0
+ evdwij = c1 + c2
+! Now cavity....................
+ eagle = dsqrt(1.0/rij_shift)
+ top = b1 * ( eagle + b2 * 1.0/rij_shift - b3 )
+ bot = 1.0d0 + b4 * (1.0/rij_shift ** 12.0d0)
+ botsq = bot * bot
+ Fcav = top / bot
+ dtop = b1 * ((1.0/ (2.0d0 * eagle)) + (b2))
+ dbot = 12.0d0 * b4 * (1.0/rij_shift) ** 11.0d0
+ dFdR = ((dtop * bot - top * dbot) / botsq)
+ w1 = wqdip_peppho(1)
+ w2 = wqdip_peppho(2)
+! w1=0.0d0
+! w2=0.0d0
+! pis = sig0head_scbase(itypi,itypj)
+! eps_head = epshead_scbase(itypi,itypj)
+!c!-------------------------------------------------------------------
+
+!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj)
+!c! & +dhead(1,1,itypi,itypj))**2))
+!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj)
+!c! & +dhead(2,1,itypi,itypj))**2))
+
+!c!-------------------------------------------------------------------
+!c! ecl
+ sparrow = w1 * om1
+ hawk = w2 * (1.0d0 - sqom1)
+ Ecl = sparrow * rij_shift**2.0d0 &
+ - hawk * rij_shift**4.0d0
+!c!-------------------------------------------------------------------
+!c! derivative of ecl is Gcl
+!c! dF/dr part
+! rij_shift=5.0
+ dGCLdR = - 2.0d0 * sparrow * rij_shift**3.0d0 &
+ + 4.0d0 * hawk * rij_shift**5.0d0
+!c! dF/dom1
+ dGCLdOM1 = (w1) * (rij_shift**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = (2.0d0 * w2 * om1) * (rij_shift ** 4.0d0)
+ eom1 = dGCLdOM1+dGCLdOM2
+ eom2 = 0.0
+
+ fac = -expon * (c1 + evdwij) * rij_shift+dFdR+dGCLdR
+! fac=0.0
+ gg(1) = fac*xj*rij
+ gg(2) = fac*yj*rij
+ gg(3) = fac*zj*rij
+ do k=1,3
+ gvdwc_peppho(k,j) = gvdwc_peppho(k,j) +gg(k)/2.0
+ gvdwc_peppho(k,j+1) = gvdwc_peppho(k,j+1) +gg(k)/2.0
+ gvdwc_peppho(k,i) = gvdwc_peppho(k,i) -gg(k)/2.0
+ gvdwc_peppho(k,i+1) = gvdwc_peppho(k,i+1) -gg(k)/2.0
+ gg(k)=0.0
+ enddo
+
+ DO k = 1, 3
+ dcosom1(k) = rij* (dc_norm(k,i) - om1 * erij(k))
+ dcosom2(k) = rij* (dc_norm(k,j) - om2 * erij(k))
+ gg(k) = gg(k) + eom1 * dcosom1(k)! + eom2 * dcosom2(k)
+ gvdwc_peppho(k,j)= gvdwc_peppho(k,j) +0.5*( gg(k)) !&
+! - (eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv*2.0
+ gvdwc_peppho(k,j+1)= gvdwc_peppho(k,j+1) +0.5*( gg(k)) !&
+! + (eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv*2.0
+ gvdwc_peppho(k,i)= gvdwc_peppho(k,i) -0.5*( gg(k)) &
+ - (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0
+ gvdwc_peppho(k,i+1)= gvdwc_peppho(k,i+1) - 0.5*( gg(k)) &
+ + (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0
+ enddo
+ epeppho=epeppho+evdwij+Fcav+ECL
+! print *,i,j,evdwij,Fcav,ECL,rij_shift
+ enddo
+ enddo
+ end subroutine eprot_pep_phosphate