+
+!c! tail distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+ Rtail = dsqrt( &
+ (Rtail_distance(1)*Rtail_distance(1)) &
+ + (Rtail_distance(2)*Rtail_distance(2)) &
+ + (Rtail_distance(3)*Rtail_distance(3)))
+! tail location and distance calculations
+! dhead1
+ d1 = dheadcat(1, 1, itypi, itypj)
+! print *,"d1",d1
+! d1=0.0d0
+! d2 = dhead(2, 1, itypi, itypj)
+ 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 + d1 * dc_norm(k, i)
+ chead(k,2) = c(k, j)
+ ENDDO
+! distance
+! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ call to_box(chead(1,1),chead(2,1),chead(3,1))
+ call to_box(chead(1,2),chead(2,2),chead(3,2))
+
+! distance
+! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ do k=1,3
+ Rhead_distance(k) = boxshift(chead(k,2) - chead(k,1),locbox(k))
+ 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)))
+!-------------------------------------------------------------------
+! 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
+ 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 = 0.0d0 ! vbld_inv(j+nres)
+! 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)
+ sss_ele_cut=sscale_ele(1.0d0/(rij))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij))
+! print *,sss_ele_cut,sss_ele_grad,&
+! 1.0d0/(rij),r_cut_ele,rlamb_ele
+ if (sss_ele_cut.le.0.0) cycle
+ 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
+ om2=0.0d0
+ om12=0.0d0
+ 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 = Rtail - sig + sig0ij
+ IF (rij_shift.le.0.0D0) THEN
+ evdw = 1.0D20
+! if (evdw.gt.1.0d6) then
+! write (*,'(2(1x,a3,i3),6f6.2)') &
+! restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
+! 1.0d0/rij,Rtail,Rhead,rij_shift, sig, sig0ij
+!evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,&
+! Equad,evdwij+Fcav+eheadtail,evdw
+! endif
+ RETURN
+ END IF
+ sigder = -sig * sigsq
+ rij_shift = 1.0D0 / rij_shift
+ fac = rij_shift**expon
+ c1 = fac * fac * aa_aq_cat(itypi,itypj)
+! print *,"ADAM",aa_aq(itypi,itypj)
+
+! c1 = 0.0d0
+ c2 = fac * bb_aq_cat(itypi,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
+!#ifdef TSCSC
+! IF (bb_aq(itypi,itypj).gt.0) THEN
+! evdw_p = evdw_p + evdwij
+! ELSE
+! evdw_m = evdw_m + evdwij
+! END IF
+!#else
+ evdw = evdw &
+ + evdwij*sss_ele_cut
+!#endif
+ c1 = c1 * eps1 * eps2rt**2 * eps3rt**2
+ fac = -expon * (c1 + evdwij) * rij_shift
+ sigder = fac * sigder
+! Calculate distance derivative
+ gg(1) = fac*sss_ele_cut+evdwij*sss_ele_grad
+ gg(2) = fac*sss_ele_cut+evdwij*sss_ele_grad
+ gg(3) = fac*sss_ele_cut+evdwij*sss_ele_grad
+
+ fac = chis1 * sqom1 + chis2 * sqom2 &
+ - 2.0d0 * chis12 * om1 * om2 * om12
+
+ pom = 1.0d0 - chis1 * chis2 * sqom12
+! print *,"TUT2",fac,chis1,sqom1,pom
+ Lambf = (1.0d0 - (fac / pom))
+ Lambf = dsqrt(Lambf)
+ sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0)
+ Chif = Rtail * sparrow
+ ChiLambf = Chif * Lambf
+ eagle = dsqrt(ChiLambf)
+ bat = ChiLambf ** 11.0d0
+ top = b1cav * ( eagle + b2cav * ChiLambf - b3cav )
+ bot = 1.0d0 + b4cav * (ChiLambf ** 12.0d0)
+ botsq = bot * bot
+ Fcav = top / bot
+
+ dtop = b1cav * ((Lambf / (2.0d0 * eagle)) + (b2cav * Lambf))
+ dbot = 12.0d0 * b4cav * bat * Lambf
+ dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow*sss_ele_cut+&
+ Fcav*sss_ele_grad
+ Fcav=Fcav*sss_ele_cut
+ dtop = b1cav * ((Chif / (2.0d0 * eagle)) + (b2cav * Chif))
+ dbot = 12.0d0 * b4cav * 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)
+ dCAVdOM1 = dFdL * ( dFdOM1 )
+! dCAVdOM2 = dFdL * ( dFdOM2 )
+! dCAVdOM12 = dFdL * ( dFdOM12 )
+ dCAVdOM2=0.0d0
+ dCAVdOM12=0.0d0
+
+ DO k= 1, 3
+ ertail(k) = Rtail_distance(k)/Rtail
+ END DO
+ erdxi = scalar( ertail(1), dC_norm(1,i) )
+ erdxj = scalar( ertail(1), dC_norm(1,j) )
+ facd1 = dtailcat(1,itypi,itypj) * vbld_inv(i)
+ facd2 = dtailcat(2,itypi,itypj) * vbld_inv(j+nres)
+ DO k = 1, 3
+ pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i))
+! gradpepcatx(k,i) = gradpepcatx(k,i) &
+! - (( dFdR + gg(k) ) * pom)
+ pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres))
+! gvdwx(k,j) = gvdwx(k,j) &
+! + (( dFdR + gg(k) ) * pom)
+ gradpepcat(k,i) = gradpepcat(k,i) &
+ - (( dFdR + gg(k) ) * ertail(k))/2.0d0
+ gradpepcat(k,i+1) = gradpepcat(k,i+1) &
+ - (( dFdR + gg(k) ) * ertail(k))/2.0d0
+
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ + (( dFdR + gg(k) ) * ertail(k))
+ gg(k) = 0.0d0
+ ENDDO
+ if (itype(j,5).gt.0) then
+!c! Compute head-head and head-tail energies for each state
+ isel = 3
+!c! Dipole-charge interactions
+ CALL edq_cat_pep(ecl, elj, epol)
+ eheadtail = ECL + elj + epol
+! print *,"i,",i,eheadtail
+! eheadtail = 0.0d0
+ else
+!HERE WATER and other types of molecules solvents will be added
+! write(iout,*) "not yet implemented"
+ CALL edd_cat_pep(ecl)
+ eheadtail=ecl
+! CALL edd_cat_pep
+! eheadtail=0.0d0
+ endif
+ evdw = evdw + Fcav + eheadtail
+! if (evdw.gt.1.0d6) then
+! write (*,'(2(1x,a3,i3),3f6.2,10f16.7)') &
+! restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
+! 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,&
+! Equad,evdwij+Fcav+eheadtail,evdw
+! endif
+ IF (energy_dec) write (iout,'(2(1x,a3,i3),3f6.2,10f16.7)') &
+ restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
+ 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,&
+ Equad,evdwij+Fcav+eheadtail,evdw
+! evdw = evdw + Fcav + eheadtail
+
+! iF (nstate(itypi,itypj).eq.1) THEN
+ CALL sc_grad_cat_pep
+! END IF
+!c!-------------------------------------------------------------------
+!c! NAPISY KONCOWE
+ END DO ! j
+! END DO ! i
+!c write (iout,*) "Number of loop steps in EGB:",ind
+!c energy_dec=.false.
+! print *,"EVDW KURW",evdw,nres
+ 23 continue
+! print *,"before leave sc_grad_cat", i,j, gradpepcat(1,nres-1)
+
+ return
+ end subroutine ecats_prot_amber
+
+!---------------------------------------------------------------------------
+! old for Ca2+
+ subroutine ecat_prot(ecation_prot)
+! use calc_data
+! use comm_momo
+ integer i,j,k,subchap,itmp,inum
+ real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,&
+ r7,r4
+ real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
+ dist_init,dist_temp,ecation_prot,rcal,rocal, &
+ Evan1,Evan2,EC,cm1mag,DASGL,delta,r0p,Epepcat, &
+ catl,cml,calpl, Etotal_p, Etotal_m,rtab,wdip,wmodquad,wquad1, &
+ wquad2,wvan1,E1,E2,wconst,wvan2,rcpm,dcmag,sin2thet,sinthet, &
+ costhet,v1m,v2m,wh2o,wc,rsecp,Ir,Irsecp,Irthrp,Irfourp,Irfiftp,&
+ Irsistp,Irseven,Irtwelv,Irthir,dE1dr,dE2dr,dEdcos,wquad2p,opt, &
+ rs,rthrp,rfourp,rsixp,reight,Irsixp,Ireight,Irtw,Irfourt, &
+ opt1,opt2,opt3,opt4,opt5,opt6,opt7,opt8,opt9,opt10,opt11,opt12,&
+ opt13,opt14,opt15,opt16,opt17,opt18,opt19, &
+ Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip,&
+ ndiv,ndivi
+ real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,&
+ gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, &
+ dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, &
+ tab1,tab2,tab3,diff,cm1,sc,p,tcat,talp,cm,drcp,drcp_norm,vcat, &
+ v1,v2,v3,myd_norm,dx,vcm,valpha,drdpep,dcosdpep,dcosddci,dEdpep,&
+ dEcCat,dEdipCm,dEdipCalp,dEquad1Cat,dEquad1Cm,dEquad1Calp, &
+ dEquad2Cat,dEquad2Cm,dEquad2Calpd,Evan1Cat,dEvan1Calp,dEvan2Cat,&
+ dEvan2Calp,dEtotalCat,dscvec,dEcCm,dEcCalp,dEdipCat,dEquad2Calp,&
+ dEvan1Cat
+ real(kind=8),dimension(6) :: vcatprm
+ ecation_prot=0.0d0
+! first lets calculate interaction with peptide groups
+ if (nres_molec(5).eq.0) return
+ itmp=0
+ do i=1,4
+ itmp=itmp+nres_molec(i)
+ enddo
+! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization
+ do i=ibond_start,ibond_end
+! cycle
+
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle ! leave dummy atoms
+ xi=0.5d0*(c(1,i)+c(1,i+1))
+ yi=0.5d0*(c(2,i)+c(2,i+1))
+ zi=0.5d0*(c(3,i)+c(3,i+1))
+ call to_box(xi,yi,zi)
+
+ do j=itmp+1,itmp+nres_molec(5)
+! print *,"WTF",itmp,j,i
+! all parameters were for Ca2+ to approximate single charge divide by two
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+ wconst=78*ndiv
+ wdip =1.092777950857032D2
+ wdip=wdip/wconst
+ wmodquad=-2.174122713004870D4
+ wmodquad=wmodquad/wconst
+ wquad1 = 3.901232068562804D1
+ wquad1=wquad1/wconst
+ wquad2 = 3
+ wquad2=wquad2/wconst
+ wvan1 = 0.1
+ wvan2 = 6
+! itmp=0
+
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ call to_box(xj,yj,zj)
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+! enddo
+! enddo
+ rcpm = sqrt(xj**2+yj**2+zj**2)
+ drcp_norm(1)=xj/rcpm
+ drcp_norm(2)=yj/rcpm
+ drcp_norm(3)=zj/rcpm
+ dcmag=0.0
+ do k=1,3
+ dcmag=dcmag+dc(k,i)**2
+ enddo
+ dcmag=dsqrt(dcmag)
+ do k=1,3
+ myd_norm(k)=dc(k,i)/dcmag
+ enddo
+ costhet=drcp_norm(1)*myd_norm(1)+drcp_norm(2)*myd_norm(2)+&
+ drcp_norm(3)*myd_norm(3)
+ rsecp = rcpm**2
+ Ir = 1.0d0/rcpm
+ Irsecp = 1.0d0/rsecp
+ Irthrp = Irsecp/rcpm
+ Irfourp = Irthrp/rcpm
+ Irfiftp = Irfourp/rcpm
+ Irsistp=Irfiftp/rcpm
+ Irseven=Irsistp/rcpm
+ Irtwelv=Irsistp*Irsistp
+ Irthir=Irtwelv/rcpm
+ sin2thet = (1-costhet*costhet)
+ sinthet=sqrt(sin2thet)
+ E1 = wdip*Irsecp*costhet+(wmodquad*Irfourp+wquad1*Irthrp)&
+ *sin2thet
+ E2 = -wquad1*Irthrp*wquad2+wvan1*(wvan2**12*Irtwelv-&
+ 2*wvan2**6*Irsistp)
+ ecation_prot = ecation_prot+E1+E2
+! print *,"ecatprot",i,j,ecation_prot,rcpm
+ dE1dr = -2*costhet*wdip*Irthrp-&
+ (4*wmodquad*Irfiftp+3*wquad1*Irfourp)*sin2thet
+ dE2dr = 3*wquad1*wquad2*Irfourp- &
+ 12*wvan1*wvan2**6*(wvan2**6*Irthir-Irseven)
+ dEdcos = wdip*Irsecp-2*(wmodquad*Irfourp+wquad1*Irthrp)*costhet
+ do k=1,3
+ drdpep(k) = -drcp_norm(k)
+ dcosdpep(k) = Ir*(costhet*drcp_norm(k)-myd_norm(k))
+ dcosddci(k) = drcp_norm(k)/dcmag-costhet*myd_norm(k)/dcmag
+ dEdpep(k) = (dE1dr+dE2dr)*drdpep(k)+dEdcos*dcosdpep(k)
+ dEddci(k) = dEdcos*dcosddci(k)
+ enddo
+ do k=1,3
+ gradpepcat(k,i)=gradpepcat(k,i)+0.5D0*dEdpep(k)-dEddci(k)
+ gradpepcat(k,i+1)=gradpepcat(k,i+1)+0.5D0*dEdpep(k)+dEddci(k)
+ gradpepcat(k,j)=gradpepcat(k,j)-dEdpep(k)
+ enddo
+ enddo ! j
+ enddo ! i
+!------------------------------------------sidechains
+! do i=1,nres_molec(1)
+ do i=ibond_start,ibond_end
+ if ((itype(i,1).eq.ntyp1)) cycle ! leave dummy atoms
+! cycle
+! print *,i,ecation_prot
+ xi=(c(1,i+nres))
+ yi=(c(2,i+nres))
+ zi=(c(3,i+nres))
+ call to_box(xi,yi,zi)
+ do k=1,3
+ cm1(k)=dc(k,i+nres)
+ enddo
+ cm1mag=sqrt(cm1(1)**2+cm1(2)**2+cm1(3)**2)
+ do j=itmp+1,itmp+nres_molec(5)
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ call to_box(xj,yj,zj)
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+! enddo
+! enddo
+! 15- Glu 16-Asp
+ if((itype(i,1).eq.15.or.itype(i,1).eq.16).or.&
+ ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.&
+ (itype(i,1).eq.25))) then
+ if(itype(i,1).eq.16) then
+ inum=1
+ else
+ inum=2
+ endif
+ do k=1,6
+ vcatprm(k)=catprm(k,inum)
+ enddo
+ dASGL=catprm(7,inum)
+! do k=1,3
+! vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres)
+ vcm(1)=(cm1(1)/cm1mag)*dASGL+xi
+ vcm(2)=(cm1(2)/cm1mag)*dASGL+yi
+ vcm(3)=(cm1(3)/cm1mag)*dASGL+zi
+
+! valpha(k)=c(k,i)
+! vcat(k)=c(k,j)
+ if (subchap.eq.1) then
+ vcat(1)=xj_temp
+ vcat(2)=yj_temp
+ vcat(3)=zj_temp
+ else
+ vcat(1)=xj_safe
+ vcat(2)=yj_safe
+ vcat(3)=zj_safe
+ endif
+ valpha(1)=xi-c(1,i+nres)+c(1,i)
+ valpha(2)=yi-c(2,i+nres)+c(2,i)
+ valpha(3)=zi-c(3,i+nres)+c(3,i)
+
+! enddo
+ do k=1,3
+ dx(k) = vcat(k)-vcm(k)
+ enddo
+ do k=1,3
+ v1(k)=(vcm(k)-valpha(k))
+ v2(k)=(vcat(k)-valpha(k))
+ enddo
+ v1m = sqrt(v1(1)**2+v1(2)**2+v1(3)**2)
+ v2m = sqrt(v2(1)**2+v2(2)**2+v2(3)**2)
+ v1dpv2 = v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
+
+! The weights of the energy function calculated from
+!The quantum mechanical GAMESS simulations of calcium with ASP/GLU
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ ndivi=0.5
+ else
+ ndivi=1.0
+ endif
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
+ wh2o=78*ndivi*ndiv
+ wc = vcatprm(1)
+ wc=wc/wh2o
+ wdip =vcatprm(2)
+ wdip=wdip/wh2o
+ wquad1 =vcatprm(3)
+ wquad1=wquad1/wh2o
+ wquad2 = vcatprm(4)
+ wquad2=wquad2/wh2o
+ wquad2p = 1.0d0-wquad2
+ wvan1 = vcatprm(5)
+ wvan2 =vcatprm(6)
+ opt = dx(1)**2+dx(2)**2
+ rsecp = opt+dx(3)**2
+ rs = sqrt(rsecp)
+ rthrp = rsecp*rs
+ rfourp = rthrp*rs
+ rsixp = rfourp*rsecp
+ reight=rsixp*rsecp
+ Ir = 1.0d0/rs
+ Irsecp = 1.0d0/rsecp
+ Irthrp = Irsecp/rs
+ Irfourp = Irthrp/rs
+ Irsixp = 1.0d0/rsixp
+ Ireight=1.0d0/reight
+ Irtw=Irsixp*Irsixp
+ Irthir=Irtw/rs
+ Irfourt=Irthir/rs
+ opt1 = (4*rs*dx(3)*wdip)
+ opt2 = 6*rsecp*wquad1*opt
+ opt3 = wquad1*wquad2p*Irsixp
+ opt4 = (wvan1*wvan2**12)
+ opt5 = opt4*12*Irfourt
+ opt6 = 2*wvan1*wvan2**6
+ opt7 = 6*opt6*Ireight
+ opt8 = wdip/v1m
+ opt10 = wdip/v2m
+ opt11 = (rsecp*v2m)**2
+ opt12 = (rsecp*v1m)**2
+ opt14 = (v1m*v2m*rsecp)**2
+ opt15 = -wquad1/v2m**2
+ opt16 = (rthrp*(v1m*v2m)**2)**2
+ opt17 = (v1m**2*rthrp)**2
+ opt18 = -wquad1/rthrp
+ opt19 = (v1m**2*v2m**2)**2
+ Ec = wc*Ir
+ do k=1,3
+ dEcCat(k) = -(dx(k)*wc)*Irthrp
+ dEcCm(k)=(dx(k)*wc)*Irthrp
+ dEcCalp(k)=0.0d0
+ enddo
+ Edip=opt8*(v1dpv2)/(rsecp*v2m)
+ do k=1,3
+ dEdipCat(k)=opt8*(v1(k)*rsecp*v2m-((v2(k)/v2m &
+ *rsecp+2*dx(k)*v2m)*v1dpv2))/opt11
+ dEdipCm(k)=opt10*(v2(k)*rsecp*v1m-((v1(k)/v1m &
+ *rsecp-2*dx(k)*v1m)*v1dpv2))/opt12
+ dEdipCalp(k)=wdip*((-v1(k)-v2(k))*rsecp*v1m &
+ *v2m-(-v1(k)/v1m*v2m*rsecp-v2(k)/v2m*v1m*rsecp) &
+ *v1dpv2)/opt14
+ enddo
+ Equad1=-wquad1*v1dpv2**2/(rthrp*(v1m*v2m)**2)
+ do k=1,3
+ dEquad1Cat(k)=-wquad1*(2*v1(k)*v1dpv2*(rthrp* &
+ (v1m*v2m)**2)-(3*dx(k)*rs*(v1m*v2m)**2+2*v1m*2* &
+ v2(k)*1/2*1/v2m*v1m*v2m*rthrp)*v1dpv2**2)/opt16
+ dEquad1Cm(k)=-wquad1*(2*v2(k)*v1dpv2*(rthrp* &
+ (v1m*v2m)**2)-(-3*dx(k)*rs*(v1m*v2m)**2+2*v2m*2* &
+ v1(k)*1/2*1/v1m*v2m*v1m*rthrp)*v1dpv2**2)/opt16
+ dEquad1Calp(k)=opt18*(2*(-v1(k)-v2(k))*v1dpv2* &
+ v1m**2*v2m**2-(-2*v1(k)*v2m**2-2*v2(k)*v1m**2)* &
+ v1dpv2**2)/opt19
+ enddo
+ Equad2=wquad1*wquad2p*Irthrp
+ do k=1,3
+ dEquad2Cat(k)=-3*dx(k)*rs*opt3
+ dEquad2Cm(k)=3*dx(k)*rs*opt3
+ dEquad2Calp(k)=0.0d0
+ enddo
+ Evan1=opt4*Irtw
+ do k=1,3
+ dEvan1Cat(k)=-dx(k)*opt5
+ dEvan1Cm(k)=dx(k)*opt5
+ dEvan1Calp(k)=0.0d0
+ enddo
+ Evan2=-opt6*Irsixp
+ do k=1,3
+ dEvan2Cat(k)=dx(k)*opt7
+ dEvan2Cm(k)=-dx(k)*opt7
+ dEvan2Calp(k)=0.0d0
+ enddo
+ ecation_prot=ecation_prot+Ec+Edip+Equad1+Equad2+Evan1+Evan2
+! print *,ecation_prot,Ec+Edip+Equad1+Equad2+Evan1+Evan2
+
+ do k=1,3
+ dEtotalCat(k)=dEcCat(k)+dEdipCat(k)+dEquad1Cat(k)+ &
+ dEquad2Cat(k)+dEvan1Cat(k)+dEvan2Cat(k)
+!c write(*,*) 'dEtotalCat inside', (dEtotalCat(l),l=1,3)
+ dEtotalCm(k)=dEcCm(k)+dEdipCm(k)+dEquad1Cm(k)+ &
+ dEquad2Cm(k)+dEvan1Cm(k)+dEvan2Cm(k)
+ dEtotalCalp(k)=dEcCalp(k)+dEdipCalp(k)+dEquad1Calp(k) &
+ +dEquad2Calp(k)+dEvan1Calp(k)+dEvan2Calp(k)
+ enddo
+ dscmag = 0.0d0
+ do k=1,3
+ dscvec(k) = dc(k,i+nres)
+ dscmag = dscmag+dscvec(k)*dscvec(k)
+ enddo
+ dscmag3 = dscmag
+ dscmag = sqrt(dscmag)
+ dscmag3 = dscmag3*dscmag
+ constA = 1.0d0+dASGL/dscmag
+ constB = 0.0d0
+ do k=1,3
+ constB = constB+dscvec(k)*dEtotalCm(k)
+ enddo
+ constB = constB*dASGL/dscmag3
+ do k=1,3
+ gg(k) = dEtotalCm(k)+dEtotalCalp(k)
+ gradpepcatx(k,i)=gradpepcatx(k,i)+ &
+ constA*dEtotalCm(k)-constB*dscvec(k)
+! print *,j,constA,dEtotalCm(k),constB,dscvec(k)
+ gradpepcat(k,i)=gradpepcat(k,i)+gg(k)
+ gradpepcat(k,j)=gradpepcat(k,j)+dEtotalCat(k)
+ enddo
+ else if (itype(i,1).eq.13.or.itype(i,1).eq.14) then
+ if(itype(i,1).eq.14) then
+ inum=3
+ else
+ inum=4
+ endif
+ do k=1,6
+ vcatprm(k)=catprm(k,inum)
+ enddo
+ dASGL=catprm(7,inum)
+! do k=1,3
+! vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres)
+! valpha(k)=c(k,i)
+! vcat(k)=c(k,j)
+! enddo
+ vcm(1)=(cm1(1)/cm1mag)*dASGL+xi
+ vcm(2)=(cm1(2)/cm1mag)*dASGL+yi
+ vcm(3)=(cm1(3)/cm1mag)*dASGL+zi
+ if (subchap.eq.1) then
+ vcat(1)=xj_temp
+ vcat(2)=yj_temp
+ vcat(3)=zj_temp
+ else
+ vcat(1)=xj_safe
+ vcat(2)=yj_safe
+ vcat(3)=zj_safe
+ endif
+ valpha(1)=xi-c(1,i+nres)+c(1,i)
+ valpha(2)=yi-c(2,i+nres)+c(2,i)
+ valpha(3)=zi-c(3,i+nres)+c(3,i)
+
+
+ do k=1,3
+ dx(k) = vcat(k)-vcm(k)
+ enddo
+ do k=1,3
+ v1(k)=(vcm(k)-valpha(k))
+ v2(k)=(vcat(k)-valpha(k))
+ enddo
+ v1m = sqrt(v1(1)**2+v1(2)**2+v1(3)**2)
+ v2m = sqrt(v2(1)**2+v2(2)**2+v2(3)**2)
+ v1dpv2 = v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
+! The weights of the energy function calculated from
+!The quantum mechanical GAMESS simulations of ASN/GLN with calcium
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
+ wh2o=78*ndiv
+ wdip =vcatprm(2)
+ wdip=wdip/wh2o
+ wquad1 =vcatprm(3)
+ wquad1=wquad1/wh2o
+ wquad2 = vcatprm(4)
+ wquad2=wquad2/wh2o
+ wquad2p = 1-wquad2
+ wvan1 = vcatprm(5)
+ wvan2 =vcatprm(6)
+ opt = dx(1)**2+dx(2)**2
+ rsecp = opt+dx(3)**2
+ rs = sqrt(rsecp)
+ rthrp = rsecp*rs
+ rfourp = rthrp*rs
+ rsixp = rfourp*rsecp
+ reight=rsixp*rsecp
+ Ir = 1.0d0/rs
+ Irsecp = 1/rsecp
+ Irthrp = Irsecp/rs
+ Irfourp = Irthrp/rs
+ Irsixp = 1/rsixp
+ Ireight=1/reight
+ Irtw=Irsixp*Irsixp
+ Irthir=Irtw/rs
+ Irfourt=Irthir/rs
+ opt1 = (4*rs*dx(3)*wdip)
+ opt2 = 6*rsecp*wquad1*opt
+ opt3 = wquad1*wquad2p*Irsixp
+ opt4 = (wvan1*wvan2**12)
+ opt5 = opt4*12*Irfourt
+ opt6 = 2*wvan1*wvan2**6
+ opt7 = 6*opt6*Ireight
+ opt8 = wdip/v1m
+ opt10 = wdip/v2m
+ opt11 = (rsecp*v2m)**2
+ opt12 = (rsecp*v1m)**2
+ opt14 = (v1m*v2m*rsecp)**2
+ opt15 = -wquad1/v2m**2
+ opt16 = (rthrp*(v1m*v2m)**2)**2
+ opt17 = (v1m**2*rthrp)**2
+ opt18 = -wquad1/rthrp
+ opt19 = (v1m**2*v2m**2)**2
+ Edip=opt8*(v1dpv2)/(rsecp*v2m)
+ do k=1,3
+ dEdipCat(k)=opt8*(v1(k)*rsecp*v2m-((v2(k)/v2m&
+ *rsecp+2*dx(k)*v2m)*v1dpv2))/opt11
+ dEdipCm(k)=opt10*(v2(k)*rsecp*v1m-((v1(k)/v1m&
+ *rsecp-2*dx(k)*v1m)*v1dpv2))/opt12
+ dEdipCalp(k)=wdip*((-v1(k)-v2(k))*rsecp*v1m&
+ *v2m-(-v1(k)/v1m*v2m*rsecp-v2(k)/v2m*v1m*rsecp)&
+ *v1dpv2)/opt14
+ enddo
+ Equad1=-wquad1*v1dpv2**2/(rthrp*(v1m*v2m)**2)
+ do k=1,3
+ dEquad1Cat(k)=-wquad1*(2*v1(k)*v1dpv2*(rthrp*&
+ (v1m*v2m)**2)-(3*dx(k)*rs*(v1m*v2m)**2+2*v1m*2*&
+ v2(k)*1/2*1/v2m*v1m*v2m*rthrp)*v1dpv2**2)/opt16
+ dEquad1Cm(k)=-wquad1*(2*v2(k)*v1dpv2*(rthrp*&
+ (v1m*v2m)**2)-(-3*dx(k)*rs*(v1m*v2m)**2+2*v2m*2*&
+ v1(k)*1/2*1/v1m*v2m*v1m*rthrp)*v1dpv2**2)/opt16
+ dEquad1Calp(k)=opt18*(2*(-v1(k)-v2(k))*v1dpv2* &
+ v1m**2*v2m**2-(-2*v1(k)*v2m**2-2*v2(k)*v1m**2)*&
+ v1dpv2**2)/opt19
+ enddo
+ Equad2=wquad1*wquad2p*Irthrp
+ do k=1,3
+ dEquad2Cat(k)=-3*dx(k)*rs*opt3
+ dEquad2Cm(k)=3*dx(k)*rs*opt3
+ dEquad2Calp(k)=0.0d0
+ enddo
+ Evan1=opt4*Irtw
+ do k=1,3
+ dEvan1Cat(k)=-dx(k)*opt5
+ dEvan1Cm(k)=dx(k)*opt5
+ dEvan1Calp(k)=0.0d0
+ enddo
+ Evan2=-opt6*Irsixp
+ do k=1,3
+ dEvan2Cat(k)=dx(k)*opt7
+ dEvan2Cm(k)=-dx(k)*opt7
+ dEvan2Calp(k)=0.0d0
+ enddo
+ ecation_prot = ecation_prot+Edip+Equad1+Equad2+Evan1+Evan2
+ do k=1,3
+ dEtotalCat(k)=dEdipCat(k)+dEquad1Cat(k)+ &
+ dEquad2Cat(k)+dEvan1Cat(k)+dEvan2Cat(k)
+ dEtotalCm(k)=dEdipCm(k)+dEquad1Cm(k)+ &
+ dEquad2Cm(k)+dEvan1Cm(k)+dEvan2Cm(k)
+ dEtotalCalp(k)=dEdipCalp(k)+dEquad1Calp(k) &
+ +dEquad2Calp(k)+dEvan1Calp(k)+dEvan2Calp(k)
+ enddo
+ dscmag = 0.0d0
+ do k=1,3
+ dscvec(k) = c(k,i+nres)-c(k,i)
+! TU SPRAWDZ???
+! dscvec(1) = xj
+! dscvec(2) = yj
+! dscvec(3) = zj
+
+ dscmag = dscmag+dscvec(k)*dscvec(k)
+ enddo
+ dscmag3 = dscmag
+ dscmag = sqrt(dscmag)
+ dscmag3 = dscmag3*dscmag
+ constA = 1+dASGL/dscmag
+ constB = 0.0d0
+ do k=1,3
+ constB = constB+dscvec(k)*dEtotalCm(k)
+ enddo
+ constB = constB*dASGL/dscmag3
+ do k=1,3
+ gg(k) = dEtotalCm(k)+dEtotalCalp(k)
+ gradpepcatx(k,i)=gradpepcatx(k,i)+ &
+ constA*dEtotalCm(k)-constB*dscvec(k)
+ gradpepcat(k,i)=gradpepcat(k,i)+gg(k)
+ gradpepcat(k,j)=gradpepcat(k,j)+dEtotalCat(k)
+ enddo
+ else
+ rcal = 0.0d0
+ do k=1,3
+! r(k) = c(k,j)-c(k,i+nres)
+ r(1) = xj
+ r(2) = yj
+ r(3) = zj
+ rcal = rcal+r(k)*r(k)
+ enddo
+ ract=sqrt(rcal)
+ rocal=1.5
+ epscalc=0.2
+ r0p=0.5*(rocal+sig0(itype(i,1)))
+ r06 = r0p**6
+ r012 = r06*r06
+ Evan1=epscalc*(r012/rcal**6)
+ Evan2=epscalc*2*(r06/rcal**3)
+ r4 = rcal**4
+ r7 = rcal**7
+ do k=1,3
+ dEvan1Cm(k) = 12*r(k)*epscalc*r012/r7
+ dEvan2Cm(k) = 12*r(k)*epscalc*r06/r4
+ enddo
+ do k=1,3
+ dEtotalCm(k)=dEvan1Cm(k)+dEvan2Cm(k)
+ enddo
+ ecation_prot = ecation_prot+ Evan1+Evan2
+ do k=1,3
+ gradpepcatx(k,i)=gradpepcatx(k,i)+ &
+ dEtotalCm(k)
+ gradpepcat(k,i)=gradpepcat(k,i)+dEtotalCm(k)
+ gradpepcat(k,j)=gradpepcat(k,j)-dEtotalCm(k)
+ enddo
+ endif ! 13-16 residues
+ enddo !j
+ enddo !i
+ return
+ end subroutine ecat_prot
+
+!----------------------------------------------------------------------------
+!---------------------------------------------------------------------------
+ subroutine ecat_nucl(ecation_nucl)
+ integer i,j,k,subchap,itmp,inum,itypi,itypj
+ real(kind=8) :: xi,yi,zi,xj,yj,zj
+ real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
+ dist_init,dist_temp,ecation_nucl,Evan1,Evan2,Ecav,Egb,wdip1,wdip2, &
+ wvan1,wvan2,wgbsig,wgbeps,wgbchi,wgbchip,wcav1,wcav2,wcav3,wcav4, &
+ wcavsig,wcavchi,v1m,v1dpdx,wh2o,wc,Edip,rcs2,invrcs6,invrcs8,invrcs12, &
+ invrcs14,rcb,rcb2,invrcb,invrcb2,invrcb4,invrcb6,cosinus,cos2,dcosdcatconst, &
+ dcosdcalpconst,dcosdcmconst,rcav,rcav11,rcav12,constcav1,constcav2, &
+ constgb1,constgb2,constdvan1,constdvan2,sgb,sgb6,sgb7,sgb12,sgb13, &
+ cavnum,cavdenom,invcavdenom2,dcavnumdcos,dcavnumdr,dcavdenomdcos, &
+ dcavdenomdr,sslipi,ssgradlipi,sslipj,ssgradlipj,aa,bb
+ real(kind=8),dimension(3) ::gg,r,dEtotalCm,dEtotalCalp,dEvan1Cm,&
+ dEvan2Cm,cm1,cm,vcat,vsug,v1,v2,dx,vcm,dEdipCm,dEdipCalp, &
+ dEvan1Calp,dEvan2Cat,dEvan2Calp,dEtotalCat,dEdipCat,dEvan1Cat,dcosdcat, &
+ dcosdcalp,dcosdcm,dEgbdCat,dEgbdCalp,dEgbdCm,dEcavdCat,dEcavdCalp, &
+ dEcavdCm,boxik
+ real(kind=8),dimension(14) :: vcatnuclprm
+ ecation_nucl=0.0d0
+ boxik(1)=boxxsize
+ boxik(2)=boxysize
+ boxik(3)=boxzsize
+
+ if (nres_molec(5).eq.0) return
+ itmp=0
+ do i=1,4
+ itmp=itmp+nres_molec(i)
+ enddo
+! print *,nres_molec(2),"nres2"
+ do i=ibond_nucl_start,ibond_nucl_end
+! do i=iatsc_s_nucl,iatsc_e_nucl
+ if ((itype(i,2).eq.ntyp1_molec(2))) cycle ! leave dummy atoms
+ xi=(c(1,i+nres))
+ yi=(c(2,i+nres))
+ zi=(c(3,i+nres))
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+ do k=1,3
+ cm1(k)=dc(k,i+nres)
+ enddo
+ do j=itmp+1,itmp+nres_molec(5)
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ call to_box(xj,yj,zj)
+! print *,i,j,itmp
+! write(iout,*) "xi,yi,zi,xj,yj,zj", xi,yi,zi,xj,yj,zj
+! call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+! aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+! +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+! bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+! +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+! write(iout,*) 'after shift', xj,yj,zj
+ dist_init=xj**2+yj**2+zj**2
+
+ itypi=itype(i,2)
+ itypj=itype(j,5)
+ do k=1,13
+ vcatnuclprm(k)=catnuclprm(k,itypi,itypj)
+ enddo
+ do k=1,3
+ vcm(k)=c(k,i+nres)
+ vsug(k)=c(k,i)
+ vcat(k)=c(k,j)
+ enddo
+ call to_box(vcm(1),vcm(2),vcm(3))
+ call to_box(vsug(1),vsug(2),vsug(3))
+ call to_box(vcat(1),vcat(2),vcat(3))
+ do k=1,3
+! dx(k) = vcat(k)-vcm(k)
+! enddo
+ dx(k)=boxshift(vcat(k)-vcm(k),boxik(k))
+! do k=1,3
+ v1(k)=dc(k,i+nres)
+ v2(k)=boxshift(vcat(k)-vsug(k),boxik(k))
+ enddo
+ v1m = sqrt(v1(1)**2+v1(2)**2+v1(3)**2)
+ v1dpdx = v1(1)*dx(1)+v1(2)*dx(2)+v1(3)*dx(3)
+! The weights of the energy function calculated from
+!The quantum mechanical Gaussian simulations of potassium and sodium with deoxynucleosides
+ wh2o=78
+ wdip1 = vcatnuclprm(1)
+ wdip1 = wdip1/wh2o !w1
+ wdip2 = vcatnuclprm(2)
+ wdip2 = wdip2/wh2o !w2
+ wvan1 = vcatnuclprm(3)
+ wvan2 = vcatnuclprm(4) !pis1
+ wgbsig = vcatnuclprm(5) !sigma0
+ wgbeps = vcatnuclprm(6) !epsi0
+ wgbchi = vcatnuclprm(7) !chi1
+ wgbchip = vcatnuclprm(8) !chip1
+ wcavsig = vcatnuclprm(9) !sig
+ wcav1 = vcatnuclprm(10) !b1
+ wcav2 = vcatnuclprm(11) !b2
+ wcav3 = vcatnuclprm(12) !b3
+ wcav4 = vcatnuclprm(13) !b4
+ wcavchi = vcatnuclprm(14) !chis1
+ rcs2 = v2(1)**2+v2(2)**2+v2(3)**2
+ invrcs6 = 1/rcs2**3
+ invrcs8 = invrcs6/rcs2
+ invrcs12 = invrcs6**2
+ invrcs14 = invrcs12/rcs2
+ rcb2 = dx(1)**2+dx(2)**2+dx(3)**2
+ rcb = sqrt(rcb2)
+ invrcb = 1/rcb
+ invrcb2 = invrcb**2
+ invrcb4 = invrcb2**2
+ invrcb6 = invrcb4*invrcb2
+ cosinus = v1dpdx/(v1m*rcb)
+ cos2 = cosinus**2
+ dcosdcatconst = invrcb2/v1m
+ dcosdcalpconst = invrcb/v1m**2
+ dcosdcmconst = invrcb2/v1m**2
+ do k=1,3
+ dcosdcat(k) = (v1(k)*rcb-dx(k)*v1m*cosinus)*dcosdcatconst
+ dcosdcalp(k) = (v1(k)*rcb*cosinus-dx(k)*v1m)*dcosdcalpconst
+ dcosdcm(k) = ((dx(k)-v1(k))*v1m*rcb+ &
+ cosinus*(dx(k)*v1m**2-v1(k)*rcb2))*dcosdcmconst
+ enddo
+ rcav = rcb/wcavsig
+ rcav11 = rcav**11
+ rcav12 = rcav11*rcav
+ constcav1 = 1-wcavchi*cos2
+ constcav2 = sqrt(constcav1)
+ constgb1 = 1/sqrt(1-wgbchi*cos2)
+ constgb2 = wgbeps*(1-wgbchip*cos2)**2
+ constdvan1 = 12*wvan1*wvan2**12*invrcs14
+ constdvan2 = 6*wvan1*wvan2**6*invrcs8
+!----------------------------------------------------------------------------
+!Gay-Berne term
+!---------------------------------------------------------------------------
+ sgb = 1/(1-constgb1+(rcb/wgbsig))
+ sgb6 = sgb**6
+ sgb7 = sgb6*sgb
+ sgb12 = sgb6**2
+ sgb13 = sgb12*sgb
+ Egb = constgb2*(sgb12-sgb6)
+ do k=1,3
+ dEgbdCat(k) = -constgb2/wgbsig*(12*sgb13-6*sgb7)*invrcb*dx(k) &
+ +(constgb1**3*constgb2*wgbchi*cosinus*(12*sgb13-6*sgb7) &
+ -4*wgbeps*wgbchip*cosinus*(1-wgbchip*cos2)*(sgb12-sgb6))*dcosdcat(k)
+ dEgbdCm(k) = constgb2/wgbsig*(12*sgb13-6*sgb7)*invrcb*dx(k) &
+ +(constgb1**3*constgb2*wgbchi*cosinus*(12*sgb13-6*sgb7) &
+ -4*wgbeps*wgbchip*cosinus*(1-wgbchip*cos2)*(sgb12-sgb6))*dcosdcm(k)
+ dEgbdCalp(k) = (constgb1**3*constgb2*wgbchi*cosinus &
+ *(12*sgb13-6*sgb7) &
+ -4*wgbeps*wgbchip*cosinus*(1-wgbchip*cos2)*(sgb12-sgb6))*dcosdcalp(k)
+ enddo
+!----------------------------------------------------------------------------
+!cavity term
+!---------------------------------------------------------------------------
+ cavnum = sqrt(rcav*constcav2)+wcav2*rcav*constcav2-wcav3
+ cavdenom = 1+wcav4*rcav12*constcav1**6
+ Ecav = wcav1*cavnum/cavdenom
+ invcavdenom2 = 1/cavdenom**2
+ dcavnumdcos = -wcavchi*cosinus/constcav2 &
+ *(sqrt(rcav/constcav2)/2+wcav2*rcav)
+ dcavnumdr = (0.5*sqrt(constcav2/rcav)+wcav2*constcav2)/wcavsig
+ dcavdenomdcos = -12*wcav4*wcavchi*rcav12*constcav1**5*cosinus
+ dcavdenomdr = 12*wcav4/wcavsig*rcav11*constcav1**6
+ do k=1,3
+ dEcavdCat(k) = ((dcavnumdcos*cavdenom-dcavdenomdcos*cavnum) &
+ *dcosdcat(k)+(dcavnumdr*cavdenom-dcavdenomdr*cavnum)/rcb*dx(k))*wcav1*invcavdenom2
+ dEcavdCm(k) = ((dcavnumdcos*cavdenom-dcavdenomdcos*cavnum) &
+ *dcosdcm(k)-(dcavnumdr*cavdenom-dcavdenomdr*cavnum)/rcb*dx(k))*wcav1*invcavdenom2
+ dEcavdCalp(k) = (dcavnumdcos*cavdenom-dcavdenomdcos*cavnum) &
+ *dcosdcalp(k)*wcav1*invcavdenom2
+ enddo
+!----------------------------------------------------------------------------
+!van der Waals and dipole-charge interaction energy
+!---------------------------------------------------------------------------
+ Evan1 = wvan1*wvan2**12*invrcs12
+ do k=1,3
+ dEvan1Cat(k) = -v2(k)*constdvan1
+ dEvan1Cm(k) = 0.0d0
+ dEvan1Calp(k) = v2(k)*constdvan1
+ enddo
+ Evan2 = -wvan1*wvan2**6*invrcs6
+ do k=1,3
+ dEvan2Cat(k) = v2(k)*constdvan2
+ dEvan2Cm(k) = 0.0d0
+ dEvan2Calp(k) = -v2(k)*constdvan2
+ enddo
+ Edip = wdip1*cosinus*invrcb2-wdip2*(1-cos2)*invrcb4
+ do k=1,3
+ dEdipCat(k) = (-2*wdip1*cosinus*invrcb4 &
+ +4*wdip2*(1-cos2)*invrcb6)*dx(k) &
+ +dcosdcat(k)*(wdip1*invrcb2+2*wdip2*cosinus*invrcb4)
+ dEdipCm(k) = (2*wdip1*cosinus*invrcb4 &
+ -4*wdip2*(1-cos2)*invrcb6)*dx(k) &
+ +dcosdcm(k)*(wdip1*invrcb2+2*wdip2*cosinus*invrcb4)
+ dEdipCalp(k) = dcosdcalp(k)*(wdip1*invrcb2 &
+ +2*wdip2*cosinus*invrcb4)
+ enddo
+ if (energy_dec) write (iout,'(2i5,4(a6,f7.3))') i,j, &
+ ' E GB ',Egb,' ECav ',Ecav,' Evdw ',Evan1+Evan2,' Edip ',Edip
+ ecation_nucl=ecation_nucl+Ecav+Egb+Edip+Evan1+Evan2
+ do k=1,3
+ dEtotalCat(k) = dEcavdCat(k)+dEvan1Cat(k)+dEvan2Cat(k) &
+ +dEgbdCat(k)+dEdipCat(k)
+ dEtotalCm(k) = dEcavdCm(k)+dEvan1Cm(k)+dEvan2Cm(k) &
+ +dEgbdCm(k)+dEdipCm(k)
+ dEtotalCalp(k) = dEcavdCalp(k)+dEgbdCalp(k)+dEvan1Calp(k) &
+ +dEdipCalp(k)+dEvan2Calp(k)
+ enddo
+ do k=1,3
+ gg(k) = dEtotalCm(k)+dEtotalCalp(k)
+ gradnuclcatx(k,i)=gradnuclcatx(k,i)+dEtotalCm(k)
+ gradnuclcat(k,i)=gradnuclcat(k,i)+gg(k)
+ gradnuclcat(k,j)=gradnuclcat(k,j)+dEtotalCat(k)
+ enddo
+ enddo !j
+ enddo !i
+ return
+ end subroutine ecat_nucl
+
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+ subroutine eprot_sc_base(escbase)
+ use calc_data
+! implicit real(kind=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) :: escbase
+ 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
+ escbase=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)
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+ 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)
+ call to_box(xj,yj,zj)
+! call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+! aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+! +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+! bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+! +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+
+ dxj = dc_norm( 1, nres+j )
+ dyj = dc_norm( 2, nres+j )
+ dzj = dc_norm( 3, nres+j )
+! print *,i,j,itypi,itypj
+ d1i = dhead_scbasei(itypi,itypj) !this is shift of dipole/charge
+ d1j = dhead_scbasej(itypi,itypj) !this is shift of dipole/charge
+! d1i=0.0d0
+! d1j=0.0d0
+! BetaT = 1.0d0 / (298.0d0 * Rb)
+! Gay-berne var's
+ sig0ij = sigma_scbase( itypi,itypj )
+ if (sig0ij.lt.0.2) print *,"KURWA",sig0ij,itypi,itypj
+ chi1 = chi_scbase( itypi, itypj,1 )
+ chi2 = chi_scbase( itypi, itypj,2 )
+! chi1=0.0d0
+! chi2=0.0d0
+ chi12 = chi1 * chi2
+ chip1 = chipp_scbase( itypi, itypj,1 )
+ chip2 = chipp_scbase( itypi, itypj,2 )
+! chip1=0.0d0
+! chip2=0.0d0
+ chip12 = chip1 * chip2
+! not used by momo potential, but needed by sc_angular which is shared
+! by all energy_potential subroutines
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+ a12sq = rborn_scbasei(itypi,itypj) * rborn_scbasej(itypi,itypj)
+! a12sq = a12sq * a12sq
+! charge of amino acid itypi is...
+ chis1 = chis_scbase(itypi,itypj,1)
+ chis2 = chis_scbase(itypi,itypj,2)
+ chis12 = chis1 * chis2
+ sig1 = sigmap1_scbase(itypi,itypj)
+ sig2 = sigmap2_scbase(itypi,itypj)
+! write (*,*) "sig1 = ", sig1
+! write (*,*) "sig2 = ", sig2
+! alpha factors from Fcav/Gcav
+ b1 = alphasur_scbase(1,itypi,itypj)
+! b1=0.0d0
+ b2 = alphasur_scbase(2,itypi,itypj)
+ b3 = alphasur_scbase(3,itypi,itypj)
+ b4 = alphasur_scbase(4,itypi,itypj)
+! used to determine whether we want to do quadrupole calculations
+! used by Fgb
+ eps_in = epsintab_scbase(itypi,itypj)
+ 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
+ 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+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)
+ 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)))
+!-------------------------------------------------------------------
+! 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
+ 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)
+! 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_scbase(itypi,itypj)
+! c1 = 0.0d0
+ c2 = fac * bb_scbase(itypi,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
+! if (b2.gt.0.0) then
+ 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=dsqrt(sig1**2.0d0 + sig2**2.0d0)
+ if (b1.eq.0.0d0) sparrow=1.0d0
+ sparrow = 1.0d0 / sparrow
+! write (*,*) "sparrow = ", sparrow,sig1,sig2,b1
+ 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
+! eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
+! eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
+! eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 &
+! -2.0D0*alf12*eps3der+sigder*sigsq_om12
+! print *,"EOMY",eom1,eom2,eom12
+! erdxi = scalar( ertail(1), dC_norm(1,i+nres) )
+! erdxj = scalar( ertail(1), dC_norm(1,j+nres) )
+! here dtail=0.0
+! facd1 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+! facd2 = dtail(2,itypi,itypj) * vbld_inv(j+nres)
+ 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))
+ gvdwx_scbase(k,i) = gvdwx_scbase(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_scbase(k,j) = gvdwx_scbase(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_scbase(k,i) = gvdwc_scbase(k,i) &
+ - (( dFdR + gg(k) ) * ertail(k))
+!c! & - ( dFdR * ertail(k))
+
+ gvdwc_scbase(k,j) = gvdwc_scbase(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
+
+! else
+
+! endif
+!Now dipole-dipole
+ if (wdipdip_scbase(2,itypi,itypj).gt.0.0d0) then
+ w1 = wdipdip_scbase(1,itypi,itypj)
+ w2 = -wdipdip_scbase(3,itypi,itypj)/2.0
+ w3 = wdipdip_scbase(2,itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! ECL
+ 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
+!c! write (*,*) "w1 = ", w1
+!c! write (*,*) "w2 = ", w2
+!c! write (*,*) "om1 = ", om1
+!c! write (*,*) "om2 = ", om2
+!c! write (*,*) "om12 = ", om12
+!c! write (*,*) "fac = ", fac
+!c! write (*,*) "c1 = ", c1
+!c! write (*,*) "c2 = ", c2
+!c! write (*,*) "Ecl = ", Ecl
+!c! write (*,*) "c2_1 = ", (w2 / Rhead ** 6.0d0)
+!c! write (*,*) "c2_2 = ",
+!c! & (4.0d0 + fac * fac -3.0d0 * (sqom1 + sqom2))
+!c!-------------------------------------------------------------------
+!c! dervative of ECL is GCL...
+!c! dECL/dr
+ 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 = d1i * vbld_inv(i+nres)
+ facd2 = d1j * vbld_inv(j+nres)
+ DO k = 1, 3
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx_scbase(k,i) = gvdwx_scbase(k,i) &
+ - dGCLdR * pom
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx_scbase(k,j) = gvdwx_scbase(k,j) &
+ + dGCLdR * pom
+
+ gvdwc_scbase(k,i) = gvdwc_scbase(k,i) &
+ - dGCLdR * erhead(k)
+ gvdwc_scbase(k,j) = gvdwc_scbase(k,j) &
+ + dGCLdR * erhead(k)
+ END DO
+ endif
+!now charge with dipole eg. ARG-dG
+ if (wqdip_scbase(2,itypi,itypj).gt.0.0d0) then
+ alphapol1 = alphapol_scbase(itypi,itypj)
+ w1 = wqdip_scbase(1,itypi,itypj)
+ w2 = wqdip_scbase(2,itypi,itypj)
+! w1=0.0d0
+! w2=0.0d0
+! pis = sig0head_scbase(itypi,itypj)
+! eps_head = epshead_scbase(itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! R1 - distance between head of ith side chain and tail of jth sidechain
+ R1 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances tail is center of side-chain
+ R1=R1+(c(k,j+nres)-chead(k,1))**2
+ END DO
+!c! Pitagoras
+ R1 = dsqrt(R1)
+
+!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!-------------------------------------------------------------------
+!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)
+!c--------------------------------------------------------------------
+!c Polarization energy
+!c Epol
+ MomoFac1 = (1.0d0 - chi1 * sqom2)
+ RR1 = R1 * R1 / MomoFac1
+ ee1 = exp(-( RR1 / (4.0d0 * a12sq) ))
+ 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
+ dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,1) = ((c(k,j+nres)-chead(k,1))/R1)
+ END DO
+
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j+nres) )
+ bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) )
+! bat=0.0d0
+ federmaus = scalar(erhead_tail(1,1),dC_norm(1,j+nres))
+ facd1 = d1i * vbld_inv(i+nres)
+ facd2 = d1j * vbld_inv(j+nres)
+! 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
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx_scbase(k,i) = gvdwx_scbase(k,i) &
+ - dGCLdR * pom &
+ - dPOLdR1 * (erhead_tail(k,1))
+! & - dGLJdR * pom
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx_scbase(k,j) = gvdwx_scbase(k,j) &
+ + dGCLdR * pom &
+ + dPOLdR1 * (erhead_tail(k,1))
+! & + dGLJdR * pom
+
+
+ gvdwc_scbase(k,i) = gvdwc_scbase(k,i) &
+ - dGCLdR * erhead(k) &
+ - dPOLdR1 * erhead_tail(k,1)
+! & - dGLJdR * erhead(k)
+
+ gvdwc_scbase(k,j) = gvdwc_scbase(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1)
+! & + dGLJdR * erhead(k)
+
+ END DO
+ endif
+! print *,i,j,evdwij,epol,Fcav,ECL
+ escbase=escbase+evdwij+epol+Fcav+ECL
+ if (energy_dec) write (iout,'(a22,2i5,4f8.3,f16.3)'), &
+ "escbase:evdw,pol,cav,CL",i,j,evdwij,epol,Fcav,ECL,escbase
+ if (energy_dec) write (iout,*) "evdwij,", evdwij, 1.0/rij, sig, sig0ij
+ call sc_grad_scbase
+ enddo
+ enddo
+
+ return
+ end subroutine eprot_sc_base
+ SUBROUTINE sc_grad_scbase
+ 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
+
+! print *,eom1,eom2,eom12,i,j,"eom1,2,12",erij(1),erij(2),erij(3)
+! 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,nres+j) - om2 * erij(k))
+ gg(k) = gg(k) + eom1 * dcosom1(k) + eom2 * dcosom2(k)
+ gvdwx_scbase(k,i)= gvdwx_scbase(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
+ gvdwx_scbase(k,j)= gvdwx_scbase(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
+ gvdwc_scbase(k,i)=gvdwc_scbase(k,i)-gg(k)
+ gvdwc_scbase(k,j)=gvdwc_scbase(k,j)+gg(k)
+ END DO
+
+ RETURN
+ END SUBROUTINE sc_grad_scbase
+
+
+ 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
+ call to_box(xi,yi,zi)
+ 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)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ dist_init=xj**2+yj**2+zj**2
+ 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
+ if (energy_dec) write (iout,'(a22,2i5,4f8.3,f16.3)'), &
+ "epepbase:evdw,pol,cav,CL",i,j,evdwij,epol,Fcav,ECL,epepbase
+ 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(kind=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,aa,bb
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,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)
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+ 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
+ call to_box(xj,yj,zj)
+! call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+! aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+! +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+! bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+! +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ 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(kind=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
+ call to_box(xi,yi,zi)
+
+ 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
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+
+ dist_init=xj**2+yj**2+zj**2
+ 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
+ if (energy_dec) write (iout,'(a22,2i5,4f8.3,f16.3)'), &
+ "epeppho:evdw,pol,cav,CL",i,j,evdwij,epol,Fcav,ECL,epeppho
+
+ epeppho=epeppho+evdwij+Fcav+ECL
+! print *,i,j,evdwij,Fcav,ECL,rij_shift
+ enddo
+ enddo
+ end subroutine eprot_pep_phosphate
+!!!!!!!!!!!!!!!!-------------------------------------------------------------
+ subroutine emomo(evdw)
+ use calc_data
+ use comm_momo
+! implicit real(kind=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,itypi1,subchap,isel,countss
+ real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,e1,e2,sigm,epsi
+ real(kind=8) :: evdw,aa,bb
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,ssgradlipi,ssgradlipj, &
+ sslipi,sslipj,faclip,alpha_sco
+ integer :: ii,icont
+ real(kind=8) :: fracinbuf
+ real (kind=8) :: escpho
+ real (kind=8),dimension(4):: ener
+ real(kind=8) :: b1,b2,egb
+ real(kind=8) :: Fisocav,ECL,Elj,Equad,Epol,eheadtail,&
+ Lambf,&
+ Chif,ChiLambf,Fcav,dFdR,dFdOM1,&
+ dFdOM2,dFdL,dFdOM12,&
+ federmaus,&
+ d1i,d1j
+! real(kind=8),dimension(3,2)::erhead_tail
+! real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead,Rtail_distance
+ real(kind=8) :: facd4, adler, Fgb, facd3
+ integer troll,jj,istate
+ real (kind=8) :: dcosom1(3),dcosom2(3)
+ evdw=0.0d0
+ eps_out=80.0d0
+ sss_ele_cut=1.0d0
+ countss=0
+! print *,"EVDW KURW",evdw,nres
+! do i=iatsc_s,iatsc_e
+! print *,"I am in EVDW",i
+ do icont=g_listscsc_start,g_listscsc_end
+ i=newcontlisti(icont)
+ j=newcontlistj(icont)
+
+ itypi=iabs(itype(i,1))
+! if (i.ne.47) cycle
+ if (itypi.eq.ntyp1) cycle
+ itypi1=iabs(itype(i+1,1))
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+! endif
+! print *, sslipi,ssgradlipi
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+! dsci_inv=dsc_inv(itypi)
+ dsci_inv=vbld_inv(i+nres)
+! write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
+! write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
+!
+! Calculate SC interaction energy.
+!
+! do iint=1,nint_gr(i)
+! do j=istart(i,iint),iend(i,iint)
+! print *,"JA PIER",i,j,iint,istart(i,iint),iend(i,iint)
+ IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+ call dyn_ssbond_ene(i,j,evdwij,countss)
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
+ 'evdw',i,j,evdwij,' ss'
+! if (energy_dec) write (iout,*) &
+! 'evdw',i,j,evdwij,' ss'
+ do k=j+1,iend(i,iint)
+!C search over all next residues
+ if (dyn_ss_mask(k)) then
+!C check if they are cysteins
+!C write(iout,*) 'k=',k
+
+!c write(iout,*) "PRZED TRI", evdwij
+! evdwij_przed_tri=evdwij
+ call triple_ssbond_ene(i,j,k,evdwij)
+!c if(evdwij_przed_tri.ne.evdwij) then
+!c write (iout,*) "TRI:", evdwij, evdwij_przed_tri
+!c endif
+
+!c write(iout,*) "PO TRI", evdwij
+!C call the energy function that removes the artifical triple disulfide
+!C bond the soubroutine is located in ssMD.F
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
+ 'evdw',i,j,evdwij,'tss'
+ endif!dyn_ss_mask(k)
+ enddo! k
+ ELSE
+!el ind=ind+1
+ itypj=iabs(itype(j,1))
+ if (itypj.eq.ntyp1) cycle
+ CALL elgrad_init(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+
+! if (j.ne.78) cycle
+! dscj_inv=dsc_inv(itypj)
+ dscj_inv=vbld_inv(j+nres)
+ xj=c(1,j+nres)
+ yj=c(2,j+nres)
+ zj=c(3,j+nres)
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+! write(iout,*) "KRUWA", i,j
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+ +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ Rreal(1)=xj
+ Rreal(2)=yj
+ Rreal(3)=zj
+ dxj = dc_norm( 1, nres+j )
+ dyj = dc_norm( 2, nres+j )
+ dzj = dc_norm( 3, nres+j )
+! print *,i,j,itypi,itypj
+! d1i=0.0d0
+! d1j=0.0d0
+! BetaT = 1.0d0 / (298.0d0 * Rb)
+! Gay-berne var's
+!1! sig0ij = sigma_scsc( itypi,itypj )
+! chi1=0.0d0
+! chi2=0.0d0
+! chip1=0.0d0
+! chip2=0.0d0
+! not used by momo potential, but needed by sc_angular which is shared
+! by all energy_potential subroutines
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+ a12sq = rborn(itypi,itypj) * rborn(itypj,itypi)
+! a12sq = a12sq * a12sq
+! charge of amino acid itypi is...
+ chis1 = chis(itypi,itypj)
+ chis2 = chis(itypj,itypi)
+ chis12 = chis1 * chis2
+ sig1 = sigmap1(itypi,itypj)
+ sig2 = sigmap2(itypi,itypj)
+! write (*,*) "sig1 = ", sig1
+! chis1=0.0
+! chis2=0.0
+! chis12 = chis1 * chis2
+! sig1=0.0
+! sig2=0.0
+! write (*,*) "sig2 = ", sig2
+! alpha factors from Fcav/Gcav
+ b1cav = alphasur(1,itypi,itypj)
+! b1cav=0.0d0
+ b2cav = alphasur(2,itypi,itypj)
+ b3cav = alphasur(3,itypi,itypj)
+ b4cav = alphasur(4,itypi,itypj)
+! used to determine whether we want to do quadrupole calculations
+ eps_in = epsintab(itypi,itypj)
+ if (eps_in.eq.0.0) eps_in=1.0
+
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+ Rtail = 0.0d0
+! dtail(1,itypi,itypj)=0.0
+! dtail(2,itypi,itypj)=0.0
+
+ DO k = 1, 3
+ ctail(k,1)=c(k,i+nres)-dtail(1,itypi,itypj)*dc_norm(k,nres+i)
+ ctail(k,2)=c(k,j+nres)-dtail(2,itypi,itypj)*dc_norm(k,nres+j)
+ END DO
+ call to_box (ctail(1,1),ctail(2,1),ctail(3,1))
+ call to_box (ctail(1,2),ctail(2,2),ctail(3,2))
+
+!c! tail distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+ Rtail_distance(1)=boxshift(ctail( 1, 2 ) - ctail( 1,1 ),boxxsize)
+ Rtail_distance(2)=boxshift(ctail( 2, 2 ) - ctail( 2,1 ),boxysize)
+ Rtail_distance(3)=boxshift(ctail( 3, 2 ) - ctail( 3,1 ),boxzsize)
+ Rtail = dsqrt( &
+ (Rtail_distance(1)*Rtail_distance(1)) &
+ + (Rtail_distance(2)*Rtail_distance(2)) &
+ + (Rtail_distance(3)*Rtail_distance(3)))
+
+! write (*,*) "eps_inout_fac = ", eps_inout_fac
+!-------------------------------------------------------------------
+! tail location and distance calculations
+ d1 = dhead(1, 1, itypi, itypj)
+ d2 = dhead(2, 1, itypi, itypj)
+
+ 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) + d1 * dc_norm(k, i+nres)
+ chead(k,2) = c(k, j+nres) + d2 * dc_norm(k, j+nres)
+! distance
+ enddo
+ if (energy_dec) write(iout,*) "before",chead(1,1),chead(2,1),chead(3,1)
+ if (energy_dec) write(iout,*) "before",chead(1,2),chead(2,2),chead(3,2)
+ call to_box (chead(1,1),chead(2,1),chead(3,1))
+ call to_box (chead(1,2),chead(2,2),chead(3,2))
+
+!c! head distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+ if (energy_dec) write(iout,*) "after",chead(1,1),chead(2,1),chead(3,1)
+ if (energy_dec) write(iout,*) "after",chead(1,2),chead(2,2),chead(3,2)
+
+ Rhead_distance(1)=boxshift(chead( 1, 2 ) - chead( 1,1 ),boxxsize)
+ Rhead_distance(2)=boxshift(chead( 2, 2 ) - chead( 2,1 ),boxysize)
+ Rhead_distance(3)=boxshift(chead( 3, 2 ) - chead( 3,1 ),boxzsize)
+ if (energy_dec) write(iout,*) "after,rdi",(Rhead_distance(k),k=1,3)
+! 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)))
+!-------------------------------------------------------------------
+! 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
+ 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)
+! 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)
+ sss_ele_cut=sscale_ele(1.0d0/(rij))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij))
+! sss_ele_cut=1.0d0
+! sss_ele_grad=0.0d0
+! print *,sss_ele_cut,sss_ele_grad,&
+! 1.0d0/(rij),r_cut_ele,rlamb_ele
+ if (sss_ele_cut.le.0.0) cycle
+
+!----------------------------
+ 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 = Rtail - 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_aq(itypi,itypj)
+! print *,"ADAM",aa_aq(itypi,itypj)
+
+! c1 = 0.0d0
+ c2 = fac * bb_aq(itypi,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
+!#ifdef TSCSC
+! IF (bb_aq(itypi,itypj).gt.0) THEN
+! evdw_p = evdw_p + evdwij
+! ELSE
+! evdw_m = evdw_m + evdwij
+! END IF
+!#else
+ evdw = evdw &
+ + evdwij*sss_ele_cut
+!#endif
+
+ 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*sss_ele_cut
+ gg(2) = fac*sss_ele_cut
+ gg(3) = fac*sss_ele_cut
+! if (b2.gt.0.0) then
+ 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))
+! print *,"fac,pom",fac,pom,Lambf
+ Lambf = dsqrt(Lambf)
+ sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0)
+! print *,"sig1,sig2",sig1,sig2,itypi,itypj
+! write (*,*) "sparrow = ", sparrow
+ Chif = Rtail * sparrow
+! print *,"rij,sparrow",rij , sparrow
+ ChiLambf = Chif * Lambf
+ eagle = dsqrt(ChiLambf)
+ bat = ChiLambf ** 11.0d0
+ top = b1cav * ( eagle + b2cav * ChiLambf - b3cav )
+ bot = 1.0d0 + b4cav * (ChiLambf ** 12.0d0)
+ botsq = bot * bot
+! print *,top,bot,"bot,top",ChiLambf,Chif
+ Fcav = top / bot
+
+ dtop = b1cav * ((Lambf / (2.0d0 * eagle)) + (b2cav * Lambf))
+ dbot = 12.0d0 * b4cav * bat * Lambf
+ dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow*sss_ele_cut
+ dtop = b1cav * ((Chif / (2.0d0 * eagle)) + (b2cav * Chif))
+ dbot = 12.0d0 * b4cav * 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 )
+
+ DO k= 1, 3
+ ertail(k) = Rtail_distance(k)/Rtail
+ END DO
+ erdxi = scalar( ertail(1), dC_norm(1,i+nres) )
+ erdxj = scalar( ertail(1), dC_norm(1,j+nres) )
+ facd1 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+ facd2 = dtail(2,itypi,itypj) * vbld_inv(j+nres)
+ DO k = 1, 3
+!c! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
+!c! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
+ pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx(k,i) = gvdwx(k,i) &
+ - (( dFdR + gg(k) ) * pom)&
+ -sss_ele_grad*Rreal(k)*rij*(Fcav+evdwij)
+!c! & - ( dFdR * pom )
+ pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx(k,j) = gvdwx(k,j) &
+ + (( dFdR + gg(k) ) * pom) &
+ +sss_ele_grad*Rreal(k)*rij*(Fcav+evdwij)
+
+!c! & + ( dFdR * pom )
+
+ gvdwc(k,i) = gvdwc(k,i) &
+ - (( dFdR + gg(k) ) * ertail(k)) &
+ -sss_ele_grad*Rreal(k)*rij*(Fcav+evdwij)
+
+!c! & - ( dFdR * ertail(k))
+
+ gvdwc(k,j) = gvdwc(k,j) &
+ + (( dFdR + gg(k) ) * ertail(k)) &
+ +sss_ele_grad*Rreal(k)*rij*(Fcav+evdwij)
+
+!c! & + ( dFdR * ertail(k))
+
+ gg(k) = 0.0d0
+! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
+! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
+ END DO
+
+
+!c! Compute head-head and head-tail energies for each state
+
+ isel = iabs(Qi) + iabs(Qj)
+! double charge for Phophorylated! itype - 25,27,27
+! if ((itype(i).eq.27).or.(itype(i).eq.26).or.(itype(i).eq.25)) then
+! Qi=Qi*2
+! Qij=Qij*2
+! endif
+! if ((itype(j).eq.27).or.(itype(j).eq.26).or.(itype(j).eq.25)) then
+! Qj=Qj*2
+! Qij=Qij*2
+! endif
+
+! isel=0
+! if (isel.eq.2) isel=0
+! if (isel.eq.3) isel=0
+! if (iabs(Qj).eq.1) isel=0
+! nstate(itypi,itypj)=1
+ IF (isel.eq.0) THEN
+!c! No charges - do nothing
+ eheadtail = 0.0d0
+
+ ELSE IF (isel.eq.4) THEN
+!c! Calculate dipole-dipole interactions
+ CALL edd(ecl)
+ eheadtail = ECL
+! eheadtail = 0.0d0
+
+ ELSE IF (isel.eq.1 .and. iabs(Qi).eq.1) THEN
+!c! Charge-nonpolar interactions
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ Qi=Qi*2
+ Qij=Qij*2
+ endif
+ if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+ Qj=Qj*2
+ Qij=Qij*2
+ endif
+
+ CALL eqn(epol)
+ eheadtail = epol
+! eheadtail = 0.0d0
+
+ ELSE IF (isel.eq.1 .and. iabs(Qj).eq.1) THEN
+!c! Nonpolar-charge interactions
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ Qi=Qi*2
+ Qij=Qij*2
+ endif
+ if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+ Qj=Qj*2
+ Qij=Qij*2
+ endif
+
+ CALL enq(epol)
+ eheadtail = epol
+! eheadtail = 0.0d0
+
+ ELSE IF (isel.eq.3 .and. icharge(itypj).eq.2) THEN
+!c! Charge-dipole interactions
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ Qi=Qi*2
+ Qij=Qij*2
+ endif
+ if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+ Qj=Qj*2
+ Qij=Qij*2
+ endif
+
+ CALL eqd(ecl, elj, epol)
+ eheadtail = ECL + elj + epol
+! eheadtail = 0.0d0
+
+ ELSE IF (isel.eq.3 .and. icharge(itypi).eq.2) THEN
+!c! Dipole-charge interactions
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ Qi=Qi*2
+ Qij=Qij*2
+ endif
+ if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+ Qj=Qj*2
+ Qij=Qij*2
+ endif
+ CALL edq(ecl, elj, epol)
+ eheadtail = ECL + elj + epol
+! eheadtail = 0.0d0
+
+ ELSE IF ((isel.eq.2.and. &
+ iabs(Qi).eq.1).and. &
+ nstate(itypi,itypj).eq.1) THEN
+!c! Same charge-charge interaction ( +/+ or -/- )
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ Qi=Qi*2
+ Qij=Qij*2
+ endif
+ if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+ Qj=Qj*2
+ Qij=Qij*2
+ endif
+
+ CALL eqq(Ecl,Egb,Epol,Fisocav,Elj)
+ eheadtail = ECL + Egb + Epol + Fisocav + Elj
+! eheadtail = 0.0d0
+
+ ELSE IF ((isel.eq.2.and. &
+ iabs(Qi).eq.1).and. &
+ nstate(itypi,itypj).ne.1) THEN
+!c! Different charge-charge interaction ( +/- or -/+ )
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ Qi=Qi*2
+ Qij=Qij*2
+ endif
+ if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+ Qj=Qj*2
+ Qij=Qij*2
+ endif
+
+ CALL energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad)
+ END IF
+ END IF ! this endif ends the "catch the gly-gly" at the beggining of Fcav
+ evdw = evdw + Fcav*sss_ele_cut + eheadtail*sss_ele_cut
+
+ IF (energy_dec) write (iout,'(2(1x,a3,i3),3f6.2,10f16.7)') &
+ restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
+ 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,&
+ Equad,evdwij+Fcav+eheadtail,evdw
+! evdw = evdw + Fcav + eheadtail
+
+ iF (nstate(itypi,itypj).eq.1) THEN
+ CALL sc_grad
+ END IF
+!c!-------------------------------------------------------------------
+!c! NAPISY KONCOWE
+ ! END DO ! j
+ !END DO ! iint
+ END DO ! i
+!c write (iout,*) "Number of loop steps in EGB:",ind
+!c energy_dec=.false.
+! print *,"EVDW KURW",evdw,nres
+
+ RETURN
+ END SUBROUTINE emomo
+!C------------------------------------------------------------------------------------
+ SUBROUTINE eqq(Ecl,Egb,Epol,Fisocav,Elj)
+ use calc_data
+ use comm_momo
+ real (kind=8) :: facd3, facd4, federmaus, adler,&
+ Ecl,Egb,Epol,Fisocav,Elj,Fgb,debkap,sgrad
+! integer :: k
+!c! Epol and Gpol analytical parameters
+ alphapol1 = alphapol(itypi,itypj)
+ alphapol2 = alphapol(itypj,itypi)
+!c! Fisocav and Gisocav analytical parameters
+ al1 = alphiso(1,itypi,itypj)
+ al2 = alphiso(2,itypi,itypj)
+ al3 = alphiso(3,itypi,itypj)
+ al4 = alphiso(4,itypi,itypj)
+ csig = (1.0d0 &
+ / dsqrt(sigiso1(itypi, itypj)**2.0d0 &
+ + sigiso2(itypi,itypj)**2.0d0))
+!c!
+ pis = sig0head(itypi,itypj)
+ eps_head = epshead(itypi,itypj)
+ Rhead_sq = Rhead * Rhead
+!c! R1 - distance between head of ith side chain and tail of jth sidechain
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R1 = 0.0d0
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances needed by Epol
+ R1=R1+(ctail(k,2)-chead(k,1))**2
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R1 = dsqrt(R1)
+ R2 = dsqrt(R2)
+
+!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! Coulomb electrostatic interaction
+ Ecl = (332.0d0 * Qij) / Rhead
+!c! derivative of Ecl is Gcl...
+ dGCLdR = (-332.0d0 * Qij ) / Rhead_sq*sss_ele_cut
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq))
+ Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0)
+ debkap=debaykap(itypi,itypj)
+ Egb = -(332.0d0 * Qij *&
+ (1.0/eps_in-dexp(-debkap*Fgb)/eps_out)) / Fgb
+! print *,"EGB WTF",Qij,eps_inout_fac,Fgb,itypi,itypj,eps_in,eps_out
+!c! Derivative of Egb is Ggb...
+ dGGBdFGB = -(-332.0d0 * Qij * &
+ (1.0/eps_in-dexp(-debkap*Fgb)/eps_out))/(Fgb*Fgb)&
+ -(332.0d0 * Qij *&
+ (dexp(-debkap*Fgb)*debkap/eps_out))/ Fgb
+ dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )/ ( 2.0d0 * Fgb )
+ dGGBdR = dGGBdFGB * dFGBdR*sss_ele_cut
+!c!-------------------------------------------------------------------
+!c! Fisocav - isotropic cavity creation term
+!c! or "how much energy it costs to put charged head in water"
+ pom = Rhead * csig
+ top = al1 * (dsqrt(pom) + al2 * pom - al3)
+ bot = (1.0d0 + al4 * pom**12.0d0)
+ botsq = bot * bot
+ FisoCav = top / bot
+! write (*,*) "Rhead = ",Rhead
+! write (*,*) "csig = ",csig
+! write (*,*) "pom = ",pom
+! write (*,*) "al1 = ",al1
+! write (*,*) "al2 = ",al2
+! write (*,*) "al3 = ",al3
+! write (*,*) "al4 = ",al4
+! write (*,*) "top = ",top
+! write (*,*) "bot = ",bot
+!c! Derivative of Fisocav is GCV...
+ dtop = al1 * ((1.0d0 / (2.0d0 * dsqrt(pom))) + al2)
+ dbot = 12.0d0 * al4 * pom ** 11.0d0
+ dGCVdR = ((dtop * bot - top * dbot) / botsq) * csig*sss_ele_cut
+!c!-------------------------------------------------------------------
+!c! Epol
+!c! Polarization energy - charged heads polarize hydrophobic "neck"
+ MomoFac1 = (1.0d0 - chi1 * sqom2)
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR1 = ( R1 * R1 ) / MomoFac1
+ RR2 = ( R2 * R2 ) / MomoFac2
+ ee1 = exp(-( RR1 / (4.0d0 * a12sq) ))
+ ee2 = exp(-( RR2 / (4.0d0 * a12sq) ))
+ fgb1 = sqrt( RR1 + a12sq * ee1 )
+ fgb2 = sqrt( RR2 + a12sq * ee2 )
+ epol = 332.0d0 * eps_inout_fac * ( &
+ (( alphapol1 / fgb1 )**4.0d0)+((alphapol2/fgb2) ** 4.0d0 ))
+!c! epol = 0.0d0
+ dPOLdFGB1 = -(1328.0d0 * eps_inout_fac * alphapol1 ** 4.0d0)&
+ / (fgb1 ** 5.0d0)
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0)&
+ / (fgb2 ** 5.0d0)
+ dFGBdR1 = ( (R1 / MomoFac1)* ( 2.0d0 - (0.5d0 * ee1) ) )&
+ / ( 2.0d0 * fgb1 )
+ dFGBdR2 = ( (R2 / MomoFac2)* ( 2.0d0 - (0.5d0 * ee2) ) )&
+ / ( 2.0d0 * fgb2 )
+ dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1))&
+ * ( 2.0d0 - 0.5d0 * ee1) ) / ( 2.0d0 * fgb1 )
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2))&
+ * ( 2.0d0 - 0.5d0 * ee2) ) / ( 2.0d0 * fgb2 )
+ dPOLdR1 = dPOLdFGB1 * dFGBdR1*sss_ele_cut
+!c! dPOLdR1 = 0.0d0
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2*sss_ele_cut
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
+!c! dPOLdOM2 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Elj
+!c! Lennard-Jones 6-12 interaction between heads
+ pom = (pis / Rhead)**6.0d0
+ Elj = 4.0d0 * eps_head * pom * (pom-1.0d0)
+!c! derivative of Elj is Glj
+ dGLJdR = 4.0d0 * eps_head*(((-12.0d0*pis**12.0d0)/(Rhead**13.0d0))&
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))*sss_ele_cut
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! These things do the dRdX derivatives, that is
+!c! allow us to change what we see from function that changes with
+!c! distance to function that changes with LOCATION (of the interaction
+!c! site)
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1)
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ END DO
+
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j+nres) )
+ bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) )
+ federmaus = scalar(erhead_tail(1,1),dC_norm(1,j+nres))
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
+ facd1 = d1 * vbld_inv(i+nres)
+ facd2 = d2 * vbld_inv(j+nres)
+ facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+ facd4 = dtail(2,itypi,itypj) * vbld_inv(j+nres)
+
+!c! Now we add appropriate partial derivatives (one in each dimension)
+ DO k = 1, 3
+ hawk = (erhead_tail(k,1) + &
+ facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres)))
+ condor = (erhead_tail(k,2) + &
+ facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j+nres)))
+ sgrad=(Ecl+Egb+Epol+Fisocav+Elj)*sss_ele_grad*rreal(k)*rij
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx(k,i) = gvdwx(k,i) &
+ - dGCLdR * pom&
+ - dGGBdR * pom&
+ - dGCVdR * pom&
+ - dPOLdR1 * hawk&
+ - dPOLdR2 * (erhead_tail(k,2)&
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))&
+ - dGLJdR * pom-sgrad
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx(k,j) = gvdwx(k,j)+ dGCLdR * pom&
+ + dGGBdR * pom+ dGCVdR * pom&
+ + dPOLdR1 * (erhead_tail(k,1)&
+ -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres)))&
+ + dPOLdR2 * condor + dGLJdR * pom+sgrad
+
+ gvdwc(k,i) = gvdwc(k,i) &
+ - dGCLdR * erhead(k)&
+ - dGGBdR * erhead(k)&
+ - dGCVdR * erhead(k)&
+ - dPOLdR1 * erhead_tail(k,1)&
+ - dPOLdR2 * erhead_tail(k,2)&
+ - dGLJdR * erhead(k)-sgrad
+
+ gvdwc(k,j) = gvdwc(k,j) &
+ + dGCLdR * erhead(k) &
+ + dGGBdR * erhead(k) &
+ + dGCVdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1) &
+ + dPOLdR2 * erhead_tail(k,2)&
+ + dGLJdR * erhead(k)+sgrad
+
+ END DO
+ RETURN
+ END SUBROUTINE eqq
+
+ SUBROUTINE eqq_cat(Ecl,Egb,Epol,Fisocav,Elj)
+ use calc_data
+ use comm_momo
+ real (kind=8) :: facd3, facd4, federmaus, adler,&
+ Ecl,Egb,Epol,Fisocav,Elj,Fgb,debkap
+! integer :: k
+!c! Epol and Gpol analytical parameters
+ alphapol1 = alphapolcat(itypi,itypj)
+ alphapol2 = alphapolcat2(itypj,itypi)
+!c! Fisocav and Gisocav analytical parameters
+ al1 = alphisocat(1,itypi,itypj)
+ al2 = alphisocat(2,itypi,itypj)
+ al3 = alphisocat(3,itypi,itypj)
+ al4 = alphisocat(4,itypi,itypj)
+ csig = (1.0d0 &
+ / dsqrt(sigiso1cat(itypi, itypj)**2.0d0 &
+ + sigiso2cat(itypi,itypj)**2.0d0))
+!c!
+ pis = sig0headcat(itypi,itypj)
+ eps_head = epsheadcat(itypi,itypj)
+ Rhead_sq = Rhead * Rhead
+!c! R1 - distance between head of ith side chain and tail of jth sidechain
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R1 = 0.0d0
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances needed by Epol
+ R1=R1+(ctail(k,2)-chead(k,1))**2
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R1 = dsqrt(R1)
+ R2 = dsqrt(R2)
+
+!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! Coulomb electrostatic interaction
+ Ecl = (332.0d0 * Qij) / Rhead
+!c! derivative of Ecl is Gcl...
+ dGCLdR = (-332.0d0 * Qij ) / Rhead_sq*sss_ele_cut+ECL*sss_ele_grad
+ ECL=ECL*sss_ele_cut
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+
+ ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq))
+ Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0)
+ debkap=debaykapcat(itypi,itypj)
+ if (energy_dec) write(iout,*) "egb",Qij,debkap,Fgb,a12sq,ee0
+ Egb = -(332.0d0 * Qij *&
+ (1.0/eps_in-dexp(-debkap*Fgb)/eps_out)) / Fgb
+! print *,"EGB WTF",Qij,eps_inout_fac,Fgb,itypi,itypj,eps_in,eps_out
+!c! Derivative of Egb is Ggb...
+ dGGBdFGB = -(-332.0d0 * Qij * &
+ (1.0/eps_in-dexp(-debkap*Fgb)/eps_out))/(Fgb*Fgb)&
+ -(332.0d0 * Qij *&
+ (dexp(-debkap*Fgb)*debkap/eps_out))/ Fgb
+ dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )/ ( 2.0d0 * Fgb )
+ dGGBdR = dGGBdFGB * dFGBdR*sss_ele_cut+Egb*sss_ele_grad
+ Egb=Egb*sss_ele_grad
+!c!-------------------------------------------------------------------
+!c! Fisocav - isotropic cavity creation term
+!c! or "how much energy it costs to put charged head in water"
+ pom = Rhead * csig
+ top = al1 * (dsqrt(pom) + al2 * pom - al3)
+ bot = (1.0d0 + al4 * pom**12.0d0)
+ botsq = bot * bot
+ FisoCav = top / bot
+! write (*,*) "Rhead = ",Rhead
+! write (*,*) "csig = ",csig
+! write (*,*) "pom = ",pom
+! write (*,*) "al1 = ",al1
+! write (*,*) "al2 = ",al2
+! write (*,*) "al3 = ",al3
+! write (*,*) "al4 = ",al4
+! write (*,*) "top = ",top
+! write (*,*) "bot = ",bot
+!c! Derivative of Fisocav is GCV...
+ dtop = al1 * ((1.0d0 / (2.0d0 * dsqrt(pom))) + al2)
+ dbot = 12.0d0 * al4 * pom ** 11.0d0
+ dGCVdR = ((dtop * bot - top * dbot) / botsq) * csig*sss_ele_cut&
+ +FisoCav*sss_ele_grad
+ FisoCav=FisoCav*sss_ele_cut
+!c!-------------------------------------------------------------------
+!c! Epol
+!c! Polarization energy - charged heads polarize hydrophobic "neck"
+ MomoFac1 = (1.0d0 - chi1 * sqom2)
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR1 = ( R1 * R1 ) / MomoFac1
+ RR2 = ( R2 * R2 ) / MomoFac2
+ ee1 = exp(-( RR1 / (4.0d0 * a12sq) ))
+ ee2 = exp(-( RR2 / (4.0d0 * a12sq) ))
+ fgb1 = sqrt( RR1 + a12sq * ee1 )
+ fgb2 = sqrt( RR2 + a12sq * ee2 )
+ epol = 332.0d0 * eps_inout_fac * ( &
+ (( alphapol1 / fgb1 )**4.0d0)+((alphapol2/fgb2) ** 4.0d0 ))
+!c! epol = 0.0d0
+ dPOLdFGB1 = -(1328.0d0 * eps_inout_fac * alphapol1 ** 4.0d0)&
+ / (fgb1 ** 5.0d0)
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0)&
+ / (fgb2 ** 5.0d0)
+ dFGBdR1 = ( (R1 / MomoFac1)* ( 2.0d0 - (0.5d0 * ee1) ) )&
+ / ( 2.0d0 * fgb1 )
+ dFGBdR2 = ( (R2 / MomoFac2)* ( 2.0d0 - (0.5d0 * ee2) ) )&
+ / ( 2.0d0 * fgb2 )
+ dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1))&
+ * ( 2.0d0 - 0.5d0 * ee1) ) / ( 2.0d0 * fgb1 )
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2))&
+ * ( 2.0d0 - 0.5d0 * ee2) ) / ( 2.0d0 * fgb2 )
+ dPOLdR1 = dPOLdFGB1 * dFGBdR1!*sss_ele_cut+epol*sss_ele_grad
+!c! dPOLdR1 = 0.0d0
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2!*sss_ele_cut+epol*sss_ele_grad
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
+! epol=epol*sss_ele_cut
+!c! dPOLdOM2 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Elj
+!c! Lennard-Jones 6-12 interaction between heads
+ pom = (pis / Rhead)**6.0d0
+ Elj = 4.0d0 * eps_head * pom * (pom-1.0d0)
+!c! derivative of Elj is Glj
+ dGLJdR = 4.0d0 * eps_head*(((-12.0d0*pis**12.0d0)/(Rhead**13.0d0))&
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))*sss_ele_cut&
+ +(Elj+epol)*sss_ele_grad
+ Elj=Elj*sss_ele_cut
+ epol=epol*sss_ele_cut
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! These things do the dRdX derivatives, that is
+!c! allow us to change what we see from function that changes with
+!c! distance to function that changes with LOCATION (of the interaction
+!c! site)
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1)
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ 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) )
+ federmaus = scalar(erhead_tail(1,1),dC_norm(1,j))
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
+ facd1 = d1 * vbld_inv(i+nres)
+ facd2 = d2 * vbld_inv(j)
+ facd3 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres)
+ facd4 = dtailcat(2,itypi,itypj) * vbld_inv(j)
+
+!c! Now we add appropriate partial derivatives (one in each dimension)
+ DO k = 1, 3
+ hawk = (erhead_tail(k,1) + &
+ facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres)))
+ condor = (erhead_tail(k,2) + &
+ facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j)))
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gradpepcatx(k,i) = gradpepcatx(k,i) &
+ - dGCLdR * pom&
+ - dGGBdR * pom&
+ - dGCVdR * pom&
+ - dPOLdR1 * hawk&
+ - dPOLdR2 * (erhead_tail(k,2)&
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))&
+ - dGLJdR * pom
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
+! gradpepcatx(k,j) = gradpepcatx(k,j)+ dGCLdR * pom&
+! + dGGBdR * pom+ dGCVdR * pom&
+! + dPOLdR1 * (erhead_tail(k,1)&
+! -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j)))&
+! + dPOLdR2 * condor + dGLJdR * pom
+
+ gradpepcat(k,i) = gradpepcat(k,i) &
+ - dGCLdR * erhead(k)&
+ - dGGBdR * erhead(k)&
+ - dGCVdR * erhead(k)&
+ - dPOLdR1 * erhead_tail(k,1)&
+ - dPOLdR2 * erhead_tail(k,2)&
+ - dGLJdR * erhead(k)
+
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ + dGCLdR * erhead(k) &
+ + dGGBdR * erhead(k) &
+ + dGCVdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1) &
+ + dPOLdR2 * erhead_tail(k,2)&
+ + dGLJdR * erhead(k)
+
+ END DO
+ RETURN
+ END SUBROUTINE eqq_cat
+!c!-------------------------------------------------------------------
+ SUBROUTINE energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad)
+ use comm_momo
+ use calc_data
+
+ double precision eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad
+ double precision ener(4)
+ double precision dcosom1(3),dcosom2(3)
+!c! used in Epol derivatives
+ double precision facd3, facd4
+ double precision federmaus, adler,sgrad
+ integer istate,ii,jj
+ real (kind=8) :: Fgb
+! print *,"CALLING EQUAD"
+!c! Epol and Gpol analytical parameters
+ alphapol1 = alphapol(itypi,itypj)
+ alphapol2 = alphapol(itypj,itypi)
+!c! Fisocav and Gisocav analytical parameters
+ al1 = alphiso(1,itypi,itypj)
+ al2 = alphiso(2,itypi,itypj)
+ al3 = alphiso(3,itypi,itypj)
+ al4 = alphiso(4,itypi,itypj)
+ csig = (1.0d0 / dsqrt(sigiso1(itypi, itypj)**2.0d0&
+ + sigiso2(itypi,itypj)**2.0d0))
+!c!
+ w1 = wqdip(1,itypi,itypj)
+ w2 = wqdip(2,itypi,itypj)
+ pis = sig0head(itypi,itypj)
+ eps_head = epshead(itypi,itypj)
+!c! First things first:
+!c! We need to do sc_grad's job with GB and Fcav
+ eom1 = eps2der * eps2rt_om1 &
+ - 2.0D0 * alf1 * eps3der&
+ + sigder * sigsq_om1&
+ + dCAVdOM1
+ eom2 = eps2der * eps2rt_om2 &
+ + 2.0D0 * alf2 * eps3der&
+ + sigder * sigsq_om2&
+ + dCAVdOM2
+ eom12 = evdwij * eps1_om12 &
+ + eps2der * eps2rt_om12 &
+ - 2.0D0 * alf12 * eps3der&
+ + sigder *sigsq_om12&
+ + dCAVdOM12
+!c! now some magical transformations to project gradient into
+!c! three cartesian vectors
+ 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))
+ gg(k) = gg(k) + eom1 * dcosom1(k) + eom2 * dcosom2(k)
+!c! this acts on hydrophobic center of interaction
+ gvdwx(k,i)= gvdwx(k,i) - gg(k)*sss_ele_cut &
+ + (eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))&
+ + eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*sss_ele_cut
+ gvdwx(k,j)= gvdwx(k,j) + gg(k)*sss_ele_cut &
+ + (eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))&
+ + eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*sss_ele_cut
+!c! this acts on Calpha
+ gvdwc(k,i)=gvdwc(k,i)-gg(k)*sss_ele_cut
+ gvdwc(k,j)=gvdwc(k,j)+gg(k)*sss_ele_cut
+ END DO
+!c! sc_grad is done, now we will compute
+ eheadtail = 0.0d0
+ eom1 = 0.0d0
+ eom2 = 0.0d0
+ eom12 = 0.0d0
+ DO istate = 1, nstate(itypi,itypj)
+!c*************************************************************
+ IF (istate.ne.1) THEN
+ IF (istate.lt.3) THEN
+ ii = 1
+ ELSE
+ ii = 2
+ END IF
+ jj = istate/ii
+ d1 = dhead(1,ii,itypi,itypj)
+ d2 = dhead(2,jj,itypi,itypj)
+ do k=1,3
+ chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres)
+ chead(k,2) = c(k, j+nres) + d2 * dc_norm(k, j+nres)
+! distance
+ enddo
+ call to_box (chead(1,1),chead(2,1),chead(3,1))
+ call to_box (chead(1,2),chead(2,2),chead(3,2))
+
+!c! head distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+
+ Rhead_distance(1)=boxshift(chead( 1, 2 ) - chead( 1,1 ),boxxsize)
+ Rhead_distance(2)=boxshift(chead( 2, 2 ) - chead( 2,1 ),boxysize)
+ Rhead_distance(3)=boxshift(chead( 3, 2 ) - chead( 3,1 ),boxzsize)
+! 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)))
+
+! DO k = 1,3
+! chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres)
+! chead(k,2) = c(k, j+nres) + d2 * dc_norm(k, j+nres)
+! Rhead_distance(k) = chead(k,2) - chead(k,1)
+! END DO
+!c! 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)))
+ END IF
+ Rhead_sq = Rhead * Rhead
+
+!c! R1 - distance between head of ith side chain and tail of jth sidechain
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R1 = 0.0d0
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R1=R1+(ctail(k,2)-chead(k,1))**2
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R1 = dsqrt(R1)
+ R2 = dsqrt(R2)
+ Ecl = (332.0d0 * Qij) / (Rhead * eps_in)
+!c! Ecl = 0.0d0
+!c! write (*,*) "Ecl = ", Ecl
+!c! derivative of Ecl is Gcl...
+ dGCLdR = (-332.0d0 * Qij ) / (Rhead_sq * eps_in)
+!c! dGCLdR = 0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Generalised Born Solvent Polarization
+ ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq))
+ Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0)
+ Egb = -(332.0d0 * Qij * eps_inout_fac) / Fgb
+!c! Egb = 0.0d0
+!c! write (*,*) "a1*a2 = ", a12sq
+!c! write (*,*) "Rhead = ", Rhead
+!c! write (*,*) "Rhead_sq = ", Rhead_sq
+!c! write (*,*) "ee = ", ee
+!c! write (*,*) "Fgb = ", Fgb
+!c! write (*,*) "fac = ", eps_inout_fac
+!c! write (*,*) "Qij = ", Qij
+!c! write (*,*) "Egb = ", Egb
+!c! Derivative of Egb is Ggb...
+!c! dFGBdR is used by Quad's later...
+ dGGBdFGB = -(-332.0d0 * Qij * eps_inout_fac) / (Fgb * Fgb)
+ dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )&
+ / ( 2.0d0 * Fgb )
+ dGGBdR = dGGBdFGB * dFGBdR
+!c! dGGBdR = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Fisocav - isotropic cavity creation term
+ pom = Rhead * csig
+ top = al1 * (dsqrt(pom) + al2 * pom - al3)
+ bot = (1.0d0 + al4 * pom**12.0d0)
+ botsq = bot * bot
+ FisoCav = top / bot
+ dtop = al1 * ((1.0d0 / (2.0d0 * dsqrt(pom))) + al2)
+ dbot = 12.0d0 * al4 * pom ** 11.0d0
+ dGCVdR = ((dtop * bot - top * dbot) / botsq) * csig
+
+!c! dGCVdR = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Polarization energy
+!c! Epol
+ MomoFac1 = (1.0d0 - chi1 * sqom2)
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR1 = ( R1 * R1 ) / MomoFac1
+ RR2 = ( R2 * R2 ) / MomoFac2
+ ee1 = exp(-( RR1 / (4.0d0 * a12sq) ))
+ ee2 = exp(-( RR2 / (4.0d0 * a12sq) ))
+ fgb1 = sqrt( RR1 + a12sq * ee1 )
+ fgb2 = sqrt( RR2 + a12sq * ee2 )
+ epol = 332.0d0 * eps_inout_fac * (&
+ (( alphapol1 / fgb1 )**4.0d0)+((alphapol2/fgb2) ** 4.0d0 ))
+!c! epol = 0.0d0
+!c! derivative of Epol is Gpol...
+ dPOLdFGB1 = -(1328.0d0 * eps_inout_fac * alphapol1 ** 4.0d0)&
+ / (fgb1 ** 5.0d0)
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0)&
+ / (fgb2 ** 5.0d0)
+ dFGBdR1 = ( (R1 / MomoFac1) &
+ * ( 2.0d0 - (0.5d0 * ee1) ) )&
+ / ( 2.0d0 * fgb1 )
+ dFGBdR2 = ( (R2 / MomoFac2) &
+ * ( 2.0d0 - (0.5d0 * ee2) ) ) &
+ / ( 2.0d0 * fgb2 )
+ dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1)) &
+ * ( 2.0d0 - 0.5d0 * ee1) ) &
+ / ( 2.0d0 * fgb1 )
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) &
+ * ( 2.0d0 - 0.5d0 * ee2) ) &
+ / ( 2.0d0 * fgb2 )
+ dPOLdR1 = dPOLdFGB1 * dFGBdR1
+!c! dPOLdR1 = 0.0d0
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
+ pom = (pis / Rhead)**6.0d0
+ Elj = 4.0d0 * eps_head * pom * (pom-1.0d0)
+!c! Elj = 0.0d0
+!c! derivative of Elj is Glj
+ dGLJdR = 4.0d0 * eps_head &
+ * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) &
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))
+!c! dGLJdR = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Equad
+ IF (Wqd.ne.0.0d0) THEN
+ Beta1 = 5.0d0 + 3.0d0 * (sqom12 - 1.0d0) &
+ - 37.5d0 * ( sqom1 + sqom2 ) &
+ + 157.5d0 * ( sqom1 * sqom2 ) &
+ - 45.0d0 * om1*om2*om12
+ fac = -( Wqd / (2.0d0 * Fgb**5.0d0) )
+ Equad = fac * Beta1
+!c! Equad = 0.0d0
+!c! derivative of Equad...
+ dQUADdR = ((2.5d0 * Wqd * Beta1) / (Fgb**6.0d0)) * dFGBdR
+!c! dQUADdR = 0.0d0
+ dQUADdOM1 = fac* (-75.0d0*om1 + 315.0d0*om1*sqom2 - 45.0d0*om2*om12)
+!c! dQUADdOM1 = 0.0d0
+ dQUADdOM2 = fac* (-75.0d0*om2 + 315.0d0*sqom1*om2 - 45.0d0*om1*om12)
+!c! dQUADdOM2 = 0.0d0
+ dQUADdOM12 = fac * ( 6.0d0*om12 - 45.0d0*om1*om2 )
+ ELSE
+ Beta1 = 0.0d0
+ Equad = 0.0d0
+ END IF
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! Angular stuff
+ eom1 = dPOLdOM1 + dQUADdOM1
+ eom2 = dPOLdOM2 + dQUADdOM2
+ eom12 = dQUADdOM12
+!c! now some magical transformations to project gradient into
+!c! three cartesian vectors
+ 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))
+ tuna(k) = eom1 * dcosom1(k) + eom2 * dcosom2(k)
+ END DO
+!c! Radial stuff
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1)
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ END DO
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j+nres) )
+ bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) )
+ federmaus = scalar(erhead_tail(1,1),dC_norm(1,j+nres))
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
+ facd1 = d1 * vbld_inv(i+nres)
+ facd2 = d2 * vbld_inv(j+nres)
+ facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+ 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))
+ condor = erhead_tail(k,2) + &
+ facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j+nres))
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+!c! this acts on hydrophobic center of interaction
+! sgrad=sss_ele_grad*(Ecl+Egb+FisoCav+epol+Elj)*rij*rreal(k)
+ gheadtail(k,1,1) = gheadtail(k,1,1) &
+ - dGCLdR * pom &
+ - dGGBdR * pom &
+ - dGCVdR * pom &
+ - dPOLdR1 * hawk &
+ - dPOLdR2 * (erhead_tail(k,2) &
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))&
+ - dGLJdR * pom &
+ - dQUADdR * pom&
+ - tuna(k) &
+ + (eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))&
+ + eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*sss_ele_cut
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+!c! this acts on hydrophobic center of interaction
+ gheadtail(k,2,1) = gheadtail(k,2,1) &
+ + dGCLdR * pom &
+ + dGGBdR * pom &
+ + dGCVdR * pom &
+ + dPOLdR1 * (erhead_tail(k,1) &
+ -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) &
+ + dPOLdR2 * condor &
+ + dGLJdR * pom &
+ + dQUADdR * pom &
+ + tuna(k) &
+ + (eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
+ + eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*sss_ele_cut
+
+!c! this acts on Calpha
+ gheadtail(k,3,1) = gheadtail(k,3,1) &
+ - dGCLdR * erhead(k)&
+ - dGGBdR * erhead(k)&
+ - dGCVdR * erhead(k)&
+ - dPOLdR1 * erhead_tail(k,1)&
+ - dPOLdR2 * erhead_tail(k,2)&
+ - dGLJdR * erhead(k) &
+ - dQUADdR * erhead(k)&
+ - tuna(k)
+!c! this acts on Calpha
+ gheadtail(k,4,1) = gheadtail(k,4,1) &
+ + dGCLdR * erhead(k) &
+ + dGGBdR * erhead(k) &
+ + dGCVdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1) &
+ + dPOLdR2 * erhead_tail(k,2) &
+ + dGLJdR * erhead(k) &
+ + dQUADdR * erhead(k)&
+ + tuna(k)
+ END DO
+ ener(istate) = ECL + Egb + Epol + Fisocav + Elj + Equad
+ eheadtail = eheadtail &
+ + wstate(istate, itypi, itypj) &
+ * dexp(-betaT * ener(istate))
+!c! foreach cartesian dimension
+ DO k = 1, 3
+!c! foreach of two gvdwx and gvdwc
+ DO l = 1, 4
+ gheadtail(k,l,2) = gheadtail(k,l,2) &
+ + wstate( istate, itypi, itypj ) &
+ * dexp(-betaT * ener(istate)) &
+ * gheadtail(k,l,1)
+ gheadtail(k,l,1) = 0.0d0
+ END DO
+ END DO
+ END DO
+!c! Here ended the gigantic DO istate = 1, 4, which starts
+!c! at the beggining of the subroutine
+
+ DO k = 1, 3
+ DO l = 1, 4
+ gheadtail(k,l,2) = gheadtail(k,l,2) / eheadtail
+ END DO
+ gvdwx(k,i) = gvdwx(k,i) + gheadtail(k,1,2)*sss_ele_cut
+ gvdwx(k,j) = gvdwx(k,j) + gheadtail(k,2,2)*sss_ele_cut
+ gvdwc(k,i) = gvdwc(k,i) + gheadtail(k,3,2)*sss_ele_cut
+ gvdwc(k,j) = gvdwc(k,j) + gheadtail(k,4,2)*sss_ele_cut
+ DO l = 1, 4
+ gheadtail(k,l,1) = 0.0d0
+ gheadtail(k,l,2) = 0.0d0
+ END DO
+ END DO
+ eheadtail = (-dlog(eheadtail)) / betaT
+ do k=1,3
+ gvdwx(k,i) = gvdwx(k,i) - eheadtail*sss_ele_grad*rreal(k)*rij
+ gvdwx(k,j) = gvdwx(k,j) + eheadtail*sss_ele_grad*rreal(k)*rij
+ gvdwc(k,i) = gvdwc(k,i) - eheadtail*sss_ele_grad*rreal(k)*rij
+ gvdwc(k,j) = gvdwc(k,j) + eheadtail*sss_ele_grad*rreal(k)*rij
+ enddo
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ dQUADdOM1 = 0.0d0
+ dQUADdOM2 = 0.0d0
+ dQUADdOM12 = 0.0d0
+ RETURN
+ END SUBROUTINE energy_quad
+!!-----------------------------------------------------------
+ SUBROUTINE eqn(Epol)
+ use comm_momo
+ use calc_data
+
+ double precision facd4, federmaus,epol
+ alphapol1 = alphapol(itypi,itypj)
+!c! R1 - distance between head of ith side chain and tail of jth sidechain
+ R1 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R1=R1+(ctail(k,2)-chead(k,1))**2
+ END DO
+!c! Pitagoras
+ R1 = dsqrt(R1)
+
+!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 Polarization energy
+!c Epol
+ MomoFac1 = (1.0d0 - chi1 * sqom2)
+ RR1 = R1 * R1 / MomoFac1
+ ee1 = exp(-( RR1 / (4.0d0 * a12sq) ))
+ fgb1 = sqrt( RR1 + a12sq * ee1)
+ epol = 332.0d0 * eps_inout_fac * (( alphapol1 / fgb1 )**4.0d0)
+ 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*sss_ele_cut
+! epol=epol*sss_ele_cut
+!c! dPOLdR1 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
+ DO k = 1, 3
+ erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1)
+ END DO
+ bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) )
+ federmaus = scalar(erhead_tail(1,1),dC_norm(1,j+nres))
+ facd1 = d1 * vbld_inv(i+nres)
+ 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)))
+
+ gvdwx(k,i) = gvdwx(k,i) &
+ - dPOLdR1 * hawk-epol*sss_ele_grad*rreal(k)*rij
+ gvdwx(k,j) = gvdwx(k,j) &
+ + dPOLdR1 * (erhead_tail(k,1) &
+ -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres)))&
+ +epol*sss_ele_grad*rreal(k)*rij
+
+ gvdwc(k,i) = gvdwc(k,i) - dPOLdR1 * erhead_tail(k,1)&
+ -epol*sss_ele_grad*rreal(k)*rij
+ gvdwc(k,j) = gvdwc(k,j) + dPOLdR1 * erhead_tail(k,1)&
+ +epol*sss_ele_grad*rreal(k)*rij
+
+ END DO
+ RETURN
+ END SUBROUTINE eqn
+ SUBROUTINE enq(Epol)
+ use calc_data
+ use comm_momo
+ double precision facd3, adler,epol
+ alphapol2 = alphapol(itypj,itypi)
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R2 = dsqrt(R2)
+
+!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 Polarization energy
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR2 = R2 * R2 / MomoFac2
+ ee2 = exp(-(RR2 / (4.0d0 * a12sq)))
+ fgb2 = sqrt(RR2 + a12sq * ee2)
+ epol = 332.0d0 * eps_inout_fac * ((alphapol2/fgb2) ** 4.0d0 )
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0) &
+ / (fgb2 ** 5.0d0)
+ dFGBdR2 = ( (R2 / MomoFac2) &
+ * ( 2.0d0 - (0.5d0 * ee2) ) ) &
+ / (2.0d0 * fgb2)
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) &
+ * (2.0d0 - 0.5d0 * ee2) ) &
+ / (2.0d0 * fgb2)
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2*sss_ele_cut
+! epol=epol*sss_ele_cut
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! (See comments in Eqq)
+ DO k = 1, 3
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ END DO
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
+ facd2 = d2 * vbld_inv(j+nres)
+ facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+ DO k = 1, 3
+ condor = (erhead_tail(k,2) &
+ + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j+nres)))
+
+ gvdwx(k,i) = gvdwx(k,i) &
+ - dPOLdR2 * (erhead_tail(k,2) &
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))&
+ -epol*sss_ele_grad*rreal(k)*rij
+ gvdwx(k,j) = gvdwx(k,j) &
+ + dPOLdR2 * condor+epol*sss_ele_grad*rreal(k)*rij
+
+
+ gvdwc(k,i) = gvdwc(k,i) &
+ - dPOLdR2 * erhead_tail(k,2)-epol*sss_ele_grad*rreal(k)*rij
+
+ gvdwc(k,j) = gvdwc(k,j) &
+ + dPOLdR2 * erhead_tail(k,2)+epol*sss_ele_grad*rreal(k)*rij
+
+
+ END DO
+ RETURN
+ END SUBROUTINE enq
+
+ SUBROUTINE enq_cat(Epol)
+ use calc_data
+ use comm_momo
+ double precision facd3, adler,epol
+ alphapol2 = alphapolcat(itypi,itypj)
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R2 = dsqrt(R2)
+
+!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 Polarization energy
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR2 = R2 * R2 / MomoFac2
+ ee2 = exp(-(RR2 / (4.0d0 * a12sq)))
+ fgb2 = sqrt(RR2 + a12sq * ee2)
+ epol = 332.0d0 * eps_inout_fac * ((alphapol2/fgb2) ** 4.0d0 )
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0) &
+ / (fgb2 ** 5.0d0)
+ dFGBdR2 = ( (R2 / MomoFac2) &
+ * ( 2.0d0 - (0.5d0 * ee2) ) ) &
+ / (2.0d0 * fgb2)
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) &
+ * (2.0d0 - 0.5d0 * ee2) ) &
+ / (2.0d0 * fgb2)
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2*sss_ele_cut+epol*sss_ele_grad
+ epol=epol*sss_ele_cut
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! (See comments in Eqq)
+ DO k = 1, 3
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ END DO
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
+ facd2 = d2 * vbld_inv(j+nres)
+ facd3 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres)
+ DO k = 1, 3
+ condor = (erhead_tail(k,2) &
+ + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j)))
+
+ gradpepcatx(k,i) = gradpepcatx(k,i) &
+ - dPOLdR2 * (erhead_tail(k,2) &
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))
+! gradpepcatx(k,j) = gradpepcatx(k,j) &
+! + dPOLdR2 * condor
+
+ gradpepcat(k,i) = gradpepcat(k,i) &
+ - dPOLdR2 * erhead_tail(k,2)
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ + dPOLdR2 * erhead_tail(k,2)
+
+ END DO
+ RETURN
+ END SUBROUTINE enq_cat
+
+ SUBROUTINE eqd(Ecl,Elj,Epol)
+ use calc_data
+ use comm_momo
+ double precision facd4, federmaus,ecl,elj,epol,sgrad
+ alphapol1 = alphapol(itypi,itypj)
+ w1 = wqdip(1,itypi,itypj)
+ w2 = wqdip(2,itypi,itypj)
+ pis = sig0head(itypi,itypj)
+ eps_head = epshead(itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! R1 - distance between head of ith side chain and tail of jth sidechain
+ R1 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R1=R1+(ctail(k,2)-chead(k,1))**2
+ END DO
+!c! Pitagoras
+ R1 = dsqrt(R1)
+
+!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 * Qi * om1
+ hawk = w2 * Qi * Qi * (1.0d0 - sqom2)
+ Ecl = sparrow / Rhead**2.0d0 &
+ - hawk / Rhead**4.0d0
+ dGCLdR = (- 2.0d0 * sparrow / Rhead**3.0d0 &
+ + 4.0d0 * hawk / Rhead**5.0d0)*sss_ele_cut
+!c! dF/dom1
+ dGCLdOM1 = (w1 * Qi) / (Rhead**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = (2.0d0 * w2 * Qi * Qi * om2) / (Rhead ** 4.0d0)
+!c--------------------------------------------------------------------
+!c Polarization energy
+!c Epol
+ MomoFac1 = (1.0d0 - chi1 * sqom2)
+ RR1 = R1 * R1 / MomoFac1
+ ee1 = exp(-( RR1 / (4.0d0 * a12sq) ))
+ fgb1 = sqrt( RR1 + a12sq * ee1)
+ epol = 332.0d0 * eps_inout_fac * (( alphapol1 / fgb1 )**4.0d0)
+!c! epol = 0.0d0
+!c!------------------------------------------------------------------
+!c! 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*sss_ele_cut
+!c! dPOLdR1 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
+!c! dPOLdOM2 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Elj
+ pom = (pis / Rhead)**6.0d0
+ Elj = 4.0d0 * eps_head * pom * (pom-1.0d0)
+!c! derivative of Elj is Glj
+ dGLJdR = 4.0d0 * eps_head &
+ * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) &
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))*sss_ele_cut
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1)
+ END DO
+
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j+nres) )
+ bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) )
+ federmaus = scalar(erhead_tail(1,1),dC_norm(1,j+nres))
+ facd1 = d1 * vbld_inv(i+nres)
+ facd2 = d2 * vbld_inv(j+nres)
+ 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)))
+ sgrad=(epol+elj+ecl)*sss_ele_grad*rreal(k)*rij
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx(k,i) = gvdwx(k,i) &
+ - dGCLdR * pom&
+ - dPOLdR1 * hawk &
+ - dGLJdR * pom &
+ -sgrad
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx(k,j) = gvdwx(k,j) &
+ + dGCLdR * pom &
+ + dPOLdR1 * (erhead_tail(k,1) &
+ -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) &
+ + dGLJdR * pom+sgrad
+
+
+ gvdwc(k,i) = gvdwc(k,i) &
+ - dGCLdR * erhead(k) &
+ - dPOLdR1 * erhead_tail(k,1) &
+ - dGLJdR * erhead(k)-sgrad
+
+ gvdwc(k,j) = gvdwc(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1) &
+ + dGLJdR * erhead(k)+sgrad
+
+ END DO
+ RETURN
+ END SUBROUTINE eqd
+
+ SUBROUTINE eqd_cat(Ecl,Elj,Epol)
+ use calc_data
+ use comm_momo
+ double precision facd4, federmaus,ecl,elj,epol
+ alphapol1 = alphapolcat(itypi,itypj)
+ w1 = wqdipcat(1,itypi,itypj)
+ w2 = wqdipcat(2,itypi,itypj)
+ pis = sig0headcat(itypi,itypj)
+ eps_head = epsheadcat(itypi,itypj)
+! eps_head=0.0d0
+! w2=0.0d0
+! alphapol1=0.0d0
+!c!-------------------------------------------------------------------
+!c! R1 - distance between head of ith side chain and tail of jth sidechain
+ R1 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R1=R1+(ctail(k,2)-chead(k,1))**2
+ END DO
+!c! Pitagoras
+ R1 = dsqrt(R1)
+
+!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 * Qi * om1
+ hawk = w2 * Qi * Qi * (1.0d0 - sqom2)
+ Ecl = sparrow / Rhead**2.0d0 &
+ - hawk / Rhead**4.0d0
+ dGCLdR =sss_ele_cut*(-2.0d0 * sparrow / Rhead**3.0d0 &
+ + 4.0d0 * hawk / Rhead**5.0d0)+sss_ele_grad*ECL
+ ECL=ECL*sss_ele_cut
+!c! dF/dom1
+ dGCLdOM1 = (w1 * Qi) / (Rhead**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = 0.0d0 !
+
+!(2.0d0 * w2 * Qi * Qi * om2) / (Rhead ** 4.0d0)
+
+!c--------------------------------------------------------------------
+!c Polarization energy
+!c Epol
+ MomoFac1 = (1.0d0 - chi1 * sqom2)
+ RR1 = R1 * R1 / MomoFac1
+ ee1 = exp(-( RR1 / (4.0d0 * a12sq) ))
+ fgb1 = sqrt( RR1 + a12sq * ee1)
+ epol = 332.0d0 * eps_inout_fac * (( alphapol1 / fgb1 )**4.0d0)
+!c! epol = 0.0d0
+!c!------------------------------------------------------------------
+!c! 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 = 0.0d0 ! as om2 is 0
+! (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1)) &
+! * (2.0d0 - 0.5d0 * ee1) ) &
+! / (2.0d0 * fgb1)
+ dPOLdR1 = dPOLdFGB1 * dFGBdR1*sss_ele_cut+epol*sss_ele_grad
+!c! dPOLdR1 = 0.0d0
+ dPOLdOM1 = 0.0d0
+! dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
+ dPOLdOM2 = 0.0d0
+ epol=epol*sss_ele_cut
+!c!-------------------------------------------------------------------
+!c! Elj
+ pom = (pis / Rhead)**6.0d0
+ Elj = 4.0d0 * eps_head * pom * (pom-1.0d0)
+!c! derivative of Elj is Glj
+ dGLJdR = 4.0d0 * eps_head*sss_ele_cut &
+ * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) &
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))+Elj*sss_ele_grad
+ Elj=Elj*sss_ele_cut
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1)
+ END DO
+
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) )
+ facd1 = d1 * vbld_inv(i+nres)
+
+ DO k = 1, 3
+ hawk = (erhead_tail(k,1) + &
+ facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres)))
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gradpepcatx(k,i) = gradpepcatx(k,i) &
+ - dGCLdR * pom&
+ - dPOLdR1 * hawk &
+ - dGLJdR * pom
+
+! pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+! gradpepcatx(k,j) = gradpepcatx(k,j) &
+! + dGCLdR * pom &
+! + dPOLdR1 * (erhead_tail(k,1) &
+! -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) &
+! + dGLJdR * pom
+
+
+ gradpepcat(k,i) = gradpepcat(k,i) &
+ - dGCLdR * erhead(k) &
+ - dPOLdR1 * erhead_tail(k,1) &
+ - dGLJdR * erhead(k)
+
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1) &
+ + dGLJdR * erhead(k)
+
+ END DO
+ RETURN
+ END SUBROUTINE eqd_cat
+
+ SUBROUTINE edq(Ecl,Elj,Epol)
+! IMPLICIT NONE
+ use comm_momo
+ use calc_data
+
+ double precision facd3, adler,ecl,elj,epol,sgrad
+ alphapol2 = alphapol(itypj,itypi)
+ w1 = wqdip(1,itypi,itypj)
+ w2 = wqdip(2,itypi,itypj)
+ pis = sig0head(itypi,itypj)
+ eps_head = epshead(itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R2 = dsqrt(R2)
+
+!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 * Qj * om1
+ hawk = w2 * Qj * Qj * (1.0d0 - sqom2)
+ ECL = sparrow / Rhead**2.0d0 &
+ - hawk / Rhead**4.0d0
+!c!-------------------------------------------------------------------
+!c! derivative of ecl is Gcl
+!c! dF/dr part
+ dGCLdR =sss_ele_cut*(- 2.0d0 * sparrow / Rhead**3.0d0 &
+ + 4.0d0 * hawk / Rhead**5.0d0)
+!c! dF/dom1
+ dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0)
+!c--------------------------------------------------------------------
+!c Polarization energy
+!c Epol
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR2 = R2 * R2 / MomoFac2
+ ee2 = exp(-(RR2 / (4.0d0 * a12sq)))
+ fgb2 = sqrt(RR2 + a12sq * ee2)
+ epol = 332.0d0 * eps_inout_fac * ((alphapol2/fgb2) ** 4.0d0 )
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0) &
+ / (fgb2 ** 5.0d0)
+ dFGBdR2 = ( (R2 / MomoFac2) &
+ * ( 2.0d0 - (0.5d0 * ee2) ) ) &
+ / (2.0d0 * fgb2)
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) &
+ * (2.0d0 - 0.5d0 * ee2) ) &
+ / (2.0d0 * fgb2)
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2*sss_ele_cut
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Elj
+ pom = (pis / Rhead)**6.0d0
+ Elj = 4.0d0 * eps_head * pom * (pom-1.0d0)
+!c! derivative of Elj is Glj
+ dGLJdR = 4.0d0 * eps_head &
+ * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) &
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))*sss_ele_cut
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! (see comments in Eqq)
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ END DO
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j+nres) )
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
+ facd1 = d1 * vbld_inv(i+nres)
+ facd2 = d2 * vbld_inv(j+nres)
+ facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+ DO k = 1, 3
+ condor = (erhead_tail(k,2) &
+ + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j+nres)))
+ sgrad=(epol+elj+ecl)*sss_ele_grad*rreal(k)*rij
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx(k,i) = gvdwx(k,i) &
+ - dGCLdR * pom &
+ - dPOLdR2 * (erhead_tail(k,2) &
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) &
+ - dGLJdR * pom-sgrad
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx(k,j) = gvdwx(k,j) &
+ + dGCLdR * pom &
+ + dPOLdR2 * condor &
+ + dGLJdR * pom+sgrad
+
+
+ gvdwc(k,i) = gvdwc(k,i) &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k)-sgrad
+
+ gvdwc(k,j) = gvdwc(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR2 * erhead_tail(k,2) &
+ + dGLJdR * erhead(k)+sgrad
+
+ END DO
+ RETURN
+ END SUBROUTINE edq
+
+ SUBROUTINE edq_cat(Ecl,Elj,Epol)
+ use comm_momo
+ use calc_data
+
+ double precision facd3, adler,ecl,elj,epol
+ alphapol2 = alphapolcat(itypi,itypj)
+ w1 = wqdipcat(1,itypi,itypj)
+ w2 = wqdipcat(2,itypi,itypj)
+ pis = sig0headcat(itypi,itypj)
+ eps_head = epsheadcat(itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R2 = dsqrt(R2)
+
+!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
+! write(iout,*) "KURWA2",Rhead
+ sparrow = w1 * Qj * om1
+ hawk = w2 * Qj * Qj * (1.0d0 - sqom2)
+ ECL = sparrow / Rhead**2.0d0 &
+ - hawk / Rhead**4.0d0
+!c!-------------------------------------------------------------------
+!c! derivative of ecl is Gcl
+!c! dF/dr part
+ dGCLdR =( - 2.0d0 * sparrow / Rhead**3.0d0 &
+ + 4.0d0 * hawk / Rhead**5.0d0)*sss_ele_cut+ECL*sss_ele_grad
+!c! dF/dom1
+ dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0)
+ ECL=ECL*sss_ele_cut
+!c--------------------------------------------------------------------
+!c--------------------------------------------------------------------
+!c Polarization energy
+!c Epol
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR2 = R2 * R2 / MomoFac2
+ ee2 = exp(-(RR2 / (4.0d0 * a12sq)))
+ fgb2 = sqrt(RR2 + a12sq * ee2)
+ epol = 332.0d0 * eps_inout_fac * ((alphapol2/fgb2) ** 4.0d0 )
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0) &
+ / (fgb2 ** 5.0d0)
+ dFGBdR2 = ( (R2 / MomoFac2) &
+ * ( 2.0d0 - (0.5d0 * ee2) ) ) &
+ / (2.0d0 * fgb2)
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) &
+ * (2.0d0 - 0.5d0 * ee2) ) &
+ / (2.0d0 * fgb2)
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2*sss_ele_cut+epol*sss_ele_grad
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ epol=epol*sss_ele_cut
+!c!-------------------------------------------------------------------
+!c! Elj
+ pom = (pis / Rhead)**6.0d0
+ Elj = 4.0d0 * eps_head * pom * (pom-1.0d0)
+!c! derivative of Elj is Glj
+ dGLJdR = 4.0d0 * eps_head &
+ * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) &
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))*sss_ele_cut+&
+ Elj*sss_ele_grad
+ Elj=Elj*sss_ele_cut
+!c!-------------------------------------------------------------------
+
+!c! Return the results
+!c! (see comments in Eqq)
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ END DO
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j) )
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
+ facd1 = d1 * vbld_inv(i+nres)
+ facd2 = d2 * vbld_inv(j)
+ facd3 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres)
+ DO k = 1, 3
+ condor = (erhead_tail(k,2) &
+ + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j)))
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gradpepcatx(k,i) = gradpepcatx(k,i) &
+ - dGCLdR * pom &
+ - dPOLdR2 * (erhead_tail(k,2) &
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) &
+ - dGLJdR * pom
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
+! gradpepcatx(k,j) = gradpepcatx(k,j) &
+! + dGCLdR * pom &
+! + dPOLdR2 * condor &
+! + dGLJdR * pom
+
+
+ gradpepcat(k,i) = gradpepcat(k,i) &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k)
+
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR2 * erhead_tail(k,2) &
+ + dGLJdR * erhead(k)
+
+ END DO
+ RETURN
+ END SUBROUTINE edq_cat
+
+ SUBROUTINE edq_cat_pep(Ecl,Elj,Epol)
+ use comm_momo
+ use calc_data
+
+ double precision facd3, adler,ecl,elj,epol
+ alphapol2 = alphapolcat(itypi,itypj)
+ w1 = wqdipcat(1,itypi,itypj)
+ w2 = wqdipcat(2,itypi,itypj)
+ pis = sig0headcat(itypi,itypj)
+ eps_head = epsheadcat(itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R2 = dsqrt(R2)
+
+!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 * Qj * om1
+ hawk = w2 * Qj * Qj * (1.0d0 - sqom2)
+! print *,"CO2", itypi,itypj
+! print *,"CO?!.", w1,w2,Qj,om1
+ ECL = sparrow / Rhead**2.0d0 &
+ - hawk / Rhead**4.0d0
+!c!-------------------------------------------------------------------
+!c! derivative of ecl is Gcl
+!c! dF/dr part
+ dGCLdR = (- 2.0d0 * sparrow / Rhead**3.0d0 &
+ + 4.0d0 * hawk / Rhead**5.0d0)*sss_ele_cut+&
+ ECL*sss_ele_grad
+ ECL=ECL*sss_ele_cut
+!c! dF/dom1
+ dGCLdOM1 = (w1 * Qj) / (Rhead**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = (2.0d0 * w2 * Qj * Qj * om2) / (Rhead ** 4.0d0)
+!c--------------------------------------------------------------------
+!c--------------------------------------------------------------------
+!c Polarization energy
+!c Epol
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR2 = R2 * R2 / MomoFac2
+ ee2 = exp(-(RR2 / (4.0d0 * a12sq)))
+ fgb2 = sqrt(RR2 + a12sq * ee2)
+ epol = 332.0d0 * eps_inout_fac * ((alphapol2/fgb2) ** 4.0d0 )
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0) &
+ / (fgb2 ** 5.0d0)
+ dFGBdR2 = ( (R2 / MomoFac2) &
+ * ( 2.0d0 - (0.5d0 * ee2) ) ) &
+ / (2.0d0 * fgb2)
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) &
+ * (2.0d0 - 0.5d0 * ee2) ) &
+ / (2.0d0 * fgb2)
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2*sss_ele_cut+epol*sss_ele_grad
+ epol=epol*sss_ele_grad
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Elj
+ pom = (pis / Rhead)**6.0d0
+ Elj = 4.0d0 * eps_head * pom * (pom-1.0d0)
+!c! derivative of Elj is Glj
+ dGLJdR = 4.0d0 * eps_head*sss_ele_cut &
+ * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) &
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))+Elj*sss_ele_grad
+ Elj=Elj*sss_ele_cut
+!c!-------------------------------------------------------------------
+
+!c! Return the results
+!c! (see comments in Eqq)
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ END DO
+ erdxi = scalar( erhead(1), dC_norm(1,i) )
+ erdxj = scalar( erhead(1), dC_norm(1,j) )
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i) )
+ facd1 = d1 * vbld_inv(i+1)/2.0
+ facd2 = d2 * vbld_inv(j)
+ facd3 = dtailcat(1,itypi,itypj) * vbld_inv(i+1)/2.0
+ DO k = 1, 3
+ condor = (erhead_tail(k,2) &
+ + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j)))
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i))
+! gradpepcatx(k,i) = gradpepcatx(k,i) &
+! - dGCLdR * pom &
+! - dPOLdR2 * (erhead_tail(k,2) &
+! -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) &
+! - dGLJdR * pom
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
+! gradpepcatx(k,j) = gradpepcatx(k,j) &
+! + dGCLdR * pom &
+! + dPOLdR2 * condor &
+! + dGLJdR * pom
+
+
+ gradpepcat(k,i) = gradpepcat(k,i) +0.5d0*( &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k))
+ gradpepcat(k,i+1) = gradpepcat(k,i+1) +0.5d0*( &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k))
+
+
+ gradpepcat(k,j) = gradpepcat(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR2 * erhead_tail(k,2) &
+ + dGLJdR * erhead(k)
+
+ END DO
+ RETURN
+ END SUBROUTINE edq_cat_pep
+
+ SUBROUTINE edd(ECL)
+! IMPLICIT NONE
+ use comm_momo
+ use calc_data
+
+ double precision ecl
+!c! csig = sigiso(itypi,itypj)
+ w1 = wqdip(1,itypi,itypj)
+ w2 = wqdip(2,itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! ECL
+ 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))
+ ECL = c1 - c2
+!c! write (*,*) "w1 = ", w1
+!c! write (*,*) "w2 = ", w2
+!c! write (*,*) "om1 = ", om1
+!c! write (*,*) "om2 = ", om2
+!c! write (*,*) "om12 = ", om12
+!c! write (*,*) "fac = ", fac
+!c! write (*,*) "c1 = ", c1
+!c! write (*,*) "c2 = ", c2
+!c! write (*,*) "Ecl = ", Ecl
+!c! write (*,*) "c2_1 = ", (w2 / Rhead ** 6.0d0)
+!c! write (*,*) "c2_2 = ",
+!c! & (4.0d0 + fac * fac -3.0d0 * (sqom1 + sqom2))
+!c!-------------------------------------------------------------------
+!c! dervative of ECL is GCL...
+!c! dECL/dr
+ c1 = (-3.0d0 * w1 * fac) / (Rhead ** 4.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead ** 7.0d0) &
+ * (4.0d0 + fac * fac - 3.0d0 * (sqom1 + sqom2))
+ dGCLdR = (c1 - c2)*sss_ele_cut!+ECL*sss_ele_grad
+!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 )
+ dGCLdOM1 = c1 - c2
+!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 )
+ dGCLdOM2 = c1 - c2
+!c! dECL/dom12
+ c1 = w1 / (Rhead ** 3.0d0)
+ c2 = ( 2.0d0 * w2 * fac ) / Rhead ** 6.0d0
+ dGCLdOM12 = c1 - c2
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! (see comments in Eqq)
+ 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(k,i) = gvdwx(k,i)- dGCLdR * pom-(ecl*sss_ele_grad*Rreal(k)*rij)
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx(k,j) = gvdwx(k,j)+ dGCLdR * pom+(ecl*sss_ele_grad*Rreal(k)*rij)
+
+ gvdwc(k,i) = gvdwc(k,i)- dGCLdR * erhead(k)-(ecl*sss_ele_grad*Rreal(k)*rij)
+ gvdwc(k,j) = gvdwc(k,j)+ dGCLdR * erhead(k)+(ecl*sss_ele_grad*Rreal(k)*rij)
+ END DO
+ RETURN
+ END SUBROUTINE edd
+ SUBROUTINE edd_cat(ECL)
+! IMPLICIT NONE
+ use comm_momo
+ use calc_data
+
+ double precision ecl
+!c! csig = sigiso(itypi,itypj)
+ w1 = wqdipcat(1,itypi,itypj)
+ w2 = wqdipcat(2,itypi,itypj)
+! w2=0.0d0
+!c!-------------------------------------------------------------------
+!c! ECL
+! print *,"om1",om1,om2,om12
+ fac = - 3.0d0 * om1 !after integer and simplify
+ c1 = (w1 / (Rhead**3.0d0)) * fac
+ c2 = (w2 / Rhead ** 6.0d0) &
+ * (4.0d0 + 6.0d0*sqom1 ) !after integration and simplification
+ ECL = c1 - c2
+!c! dervative of ECL is GCL...
+!c! dECL/dr
+ c1 = (-3.0d0 * w1 * fac) / (Rhead ** 4.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead ** 7.0d0) &
+ * (4.0d0 + 6.0d0*sqom1)
+ dGCLdR = (c1 - c2)*sss_ele_cut+ECL*sss_ele_grad
+!c! dECL/dom1
+ c1 = (-3.0d0 * w1) / (Rhead**3.0d0)
+ c2 = (12.0d0 * w2*om1) / (Rhead**6.0d0)
+ dGCLdOM1 = c1 - c2
+!c! dECL/dom2
+! c1 = (-3.0d0 * w1 * om1 ) / (Rhead**3.0d0)
+ c1=0.0 ! this is because om2 is 0
+! c2 = (-6.0d0 * w2) / (Rhead**6.0d0) &
+! * ( om1 * om12 - 3.0d0 * sqom1 * om2 + om2 )
+ c2=0.0 !om is 0
+ dGCLdOM2 = c1 - c2
+!c! dECL/dom12
+! c1 = w1 / (Rhead ** 3.0d0)
+ c1=0.0d0 ! this is because om12 is 0
+! c2 = ( 2.0d0 * w2 * fac ) / Rhead ** 6.0d0
+ c2=0.0d0 !om12 is 0
+ dGCLdOM12 = c1 - c2
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! (see comments in Eqq)
+ 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))
+ gradpepcatx(k,i) = gradpepcatx(k,i) - dGCLdR * pom
+! pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+! gradpepcatx(k,j) = gradpepcatx(k,j) + dGCLdR * pom
+
+ gradpepcat(k,i) = gradpepcat(k,i) - dGCLdR * erhead(k)
+ gradpepcat(k,j) = gradpepcat(k,j) + dGCLdR * erhead(k)
+ END DO
+ RETURN
+ END SUBROUTINE edd_cat
+ SUBROUTINE edd_cat_pep(ECL)
+! IMPLICIT NONE
+ use comm_momo
+ use calc_data
+
+ double precision ecl
+!c! csig = sigiso(itypi,itypj)
+ w1 = wqdipcat(1,itypi,itypj)
+ w2 = wqdipcat(2,itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! ECL
+ 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))
+ ECL = c1 - c2
+!c! dECL/dr
+ c1 = (-3.0d0 * w1 * fac) / (Rhead ** 4.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead ** 7.0d0) &
+ * (4.0d0 + fac * fac - 3.0d0 * (sqom1 + sqom2))
+ dGCLdR = (c1 - c2)*sss_ele_cut+ECL*sss_ele_grad
+ ECL=ECL*sss_ele_cut
+!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 )
+ dGCLdOM1 = c1 - c2
+!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 )
+ dGCLdOM2 = c1 - c2
+ dGCLdOM2=0.0d0 ! this is because om2=0
+!c! dECL/dom12
+ c1 = w1 / (Rhead ** 3.0d0)
+ c2 = ( 2.0d0 * w2 * fac ) / Rhead ** 6.0d0
+ dGCLdOM12 = c1 - c2
+ dGCLdOM12=0.0d0 !this is because om12=0.0
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! (see comments in Eqq)
+ DO k= 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ END DO
+ erdxi = scalar( erhead(1), dC_norm(1,i) )
+ erdxj = scalar( erhead(1), dC_norm(1,j+nres) )
+ facd1 = d1 * vbld_inv(i)
+ facd2 = d2 * vbld_inv(j+nres)
+ DO k = 1, 3
+
+ pom = facd1*(erhead(k)-erdxi*dC_norm(k,i))
+ gradpepcat(k,i) = gradpepcat(k,i) + dGCLdR * pom
+ gradpepcat(k,i+1) = gradpepcat(k,i+1) - dGCLdR * pom
+! pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+! gradpepcatx(k,j) = gradpepcatx(k,j) + dGCLdR * pom
+
+ gradpepcat(k,i) = gradpepcat(k,i) - dGCLdR * erhead(k)*0.5d0
+ gradpepcat(k,i+1) = gradpepcat(k,i+1)- dGCLdR * erhead(k)*0.5d0
+ gradpepcat(k,j) = gradpepcat(k,j) + dGCLdR * erhead(k)
+ END DO
+ RETURN
+ END SUBROUTINE edd_cat_pep
+
+ SUBROUTINE elgrad_init(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+! IMPLICIT NONE
+ use comm_momo
+ use calc_data
+
+ real(kind=8) :: eheadtail,Egb,Ecl,Elj,Equad,Epol,Rb
+ eps_out=80.0d0
+ itypi = itype(i,1)
+ itypj = itype(j,1)
+!c! 1/(Gas Constant * Thermostate temperature) = BetaT
+!c! ENABLE THIS LINE WHEN USING CHECKGRAD!!!
+!c! t_bath = 300
+!c! BetaT = 1.0d0 / (t_bath * Rb)i
+ Rb=0.001986d0
+ BetaT = 1.0d0 / (298.0d0 * Rb)
+!c! Gay-berne var's
+ sig0ij = sigma( itypi,itypj )
+ chi1 = chi( itypi, itypj )
+ chi2 = chi( itypj, itypi )
+ chi12 = chi1 * chi2
+ chip1 = chipp( itypi, itypj )
+ chip2 = chipp( itypj, itypi )
+ chip12 = chip1 * chip2
+! chi1=0.0
+! chi2=0.0
+! chi12=0.0
+! chip1=0.0
+! chip2=0.0
+! chip12=0.0
+!c! not used by momo potential, but needed by sc_angular which is shared
+!c! by all energy_potential subroutines
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+!c! location, location, location
+! xj = c( 1, nres+j ) - xi
+! yj = c( 2, nres+j ) - yi
+! zj = c( 3, nres+j ) - zi
+ dxj = dc_norm( 1, nres+j )
+ dyj = dc_norm( 2, nres+j )
+ dzj = dc_norm( 3, nres+j )
+!c! distance from center of chain(?) to polar/charged head
+!c! write (*,*) "istate = ", 1
+!c! write (*,*) "ii = ", 1
+!c! write (*,*) "jj = ", 1
+ d1 = dhead(1, 1, itypi, itypj)
+ d2 = dhead(2, 1, itypi, itypj)
+!c! ai*aj from Fgb
+ a12sq = rborn(itypi,itypj) * rborn(itypj,itypi)
+!c! a12sq = a12sq * a12sq
+!c! charge of amino acid itypi is...
+ Qi = icharge(itypi)
+ Qj = icharge(itypj)
+ Qij = Qi * Qj
+!c! chis1,2,12
+ chis1 = chis(itypi,itypj)
+ chis2 = chis(itypj,itypi)
+ chis12 = chis1 * chis2
+ sig1 = sigmap1(itypi,itypj)
+ sig2 = sigmap2(itypi,itypj)
+!c! write (*,*) "sig1 = ", sig1
+!c! write (*,*) "sig2 = ", sig2
+!c! alpha factors from Fcav/Gcav
+ b1cav = alphasur(1,itypi,itypj)
+! b1cav=0.0
+ b2cav = alphasur(2,itypi,itypj)
+ b3cav = alphasur(3,itypi,itypj)
+ b4cav = alphasur(4,itypi,itypj)
+ wqd = wquad(itypi, itypj)
+!c! used by Fgb
+ eps_in = epsintab(itypi,itypj)
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+!c! write (*,*) "eps_inout_fac = ", eps_inout_fac
+!c!-------------------------------------------------------------------
+!c! tail location and distance calculations
+ Rtail = 0.0d0
+ DO k = 1, 3
+ ctail(k,1)=c(k,i+nres)-dtail(1,itypi,itypj)*dc_norm(k,nres+i)
+ ctail(k,2)=c(k,j+nres)-dtail(2,itypi,itypj)*dc_norm(k,nres+j)
+ END DO
+!c! tail distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+ Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 )
+ Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 )
+ Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 )
+ Rtail = dsqrt( &
+ (Rtail_distance(1)*Rtail_distance(1)) &
+ + (Rtail_distance(2)*Rtail_distance(2)) &
+ + (Rtail_distance(3)*Rtail_distance(3)))
+!c!-------------------------------------------------------------------
+!c! Calculate location and distance between polar heads
+!c! distance between heads
+!c! for each one of our three dimensional space...
+ d1 = dhead(1, 1, itypi, itypj)
+ d2 = dhead(2, 1, itypi, itypj)
+
+ DO k = 1,3
+!c! location of polar head is computed by taking hydrophobic centre
+!c! and moving by a d1 * dc_norm vector
+!c! see unres publications for very informative images
+ chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres)
+ chead(k,2) = c(k, j+nres) + d2 * dc_norm(k, j+nres)
+!c! distance
+!c! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+!c! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ Rhead_distance(k) = chead(k,2) - chead(k,1)
+ END DO
+!c! 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)))
+!c!-------------------------------------------------------------------
+!c! zero everything that should be zero'ed
+ Egb = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ eheadtail = 0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ RETURN
+ END SUBROUTINE elgrad_init
+
+
+ SUBROUTINE elgrad_init_cat(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+ use comm_momo
+ use calc_data
+ real(kind=8) :: eheadtail,Egb,Ecl,Elj,Equad,Epol,Rb
+ eps_out=80.0d0
+ itypi = itype(i,1)
+ itypj = itype(j,5)
+!c! 1/(Gas Constant * Thermostate temperature) = BetaT
+!c! ENABLE THIS LINE WHEN USING CHECKGRAD!!!
+!c! t_bath = 300
+!c! BetaT = 1.0d0 / (t_bath * Rb)i
+ Rb=0.001986d0
+ BetaT = 1.0d0 / (298.0d0 * Rb)
+!c! Gay-berne var's
+ sig0ij = sigmacat( itypi,itypj )
+ chi1 = chi1cat( itypi, itypj )
+ chi2 = 0.0d0
+ chi12 = 0.0d0
+ chip1 = chipp1cat( itypi, itypj )
+ chip2 = 0.0d0
+ chip12 = 0.0d0
+!c! not used by momo potential, but needed by sc_angular which is shared
+!c! by all energy_potential subroutines
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+ dxj = 0.0d0 !dc_norm( 1, nres+j )
+ dyj = 0.0d0 !dc_norm( 2, nres+j )
+ dzj = 0.0d0 !dc_norm( 3, nres+j )
+!c! distance from center of chain(?) to polar/charged head
+ d1 = dheadcat(1, 1, itypi, itypj)
+ d2 = dheadcat(2, 1, itypi, itypj)
+!c! ai*aj from Fgb
+ a12sq = rborn1cat(itypi,itypj) * rborn2cat(itypi,itypj)
+!c! a12sq = a12sq * a12sq
+!c! charge of amino acid itypi is...
+ Qi = icharge(itypi)
+ Qj = ichargecat(itypj)
+ Qij = Qi * Qj
+!c! chis1,2,12
+ chis1 = chis1cat(itypi,itypj)
+ chis2 = 0.0d0
+ chis12 = 0.0d0
+ sig1 = sigmap1cat(itypi,itypj)
+ sig2 = sigmap2cat(itypi,itypj)
+!c! alpha factors from Fcav/Gcav
+ b1cav = alphasurcat(1,itypi,itypj)
+ b2cav = alphasurcat(2,itypi,itypj)
+ b3cav = alphasurcat(3,itypi,itypj)
+ b4cav = alphasurcat(4,itypi,itypj)
+ wqd = wquadcat(itypi, itypj)
+!c! used by Fgb
+ eps_in = epsintabcat(itypi,itypj)
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+!c!-------------------------------------------------------------------
+!c! tail location and distance calculations
+ Rtail = 0.0d0
+ DO k = 1, 3
+ ctail(k,1)=c(k,i+nres)-dtailcat(1,itypi,itypj)*dc_norm(k,nres+i)
+ ctail(k,2)=c(k,j)!-dtailcat(2,itypi,itypj)*dc_norm(k,nres+j)
+ END DO
+!c! tail distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+ Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 )
+ Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 )
+ Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 )
+ Rtail = dsqrt( &
+ (Rtail_distance(1)*Rtail_distance(1)) &
+ + (Rtail_distance(2)*Rtail_distance(2)) &
+ + (Rtail_distance(3)*Rtail_distance(3)))
+!c!-------------------------------------------------------------------
+!c! Calculate location and distance between polar heads
+!c! distance between heads
+!c! for each one of our three dimensional space...
+ d1 = dheadcat(1, 1, itypi, itypj)
+ d2 = dheadcat(2, 1, itypi, itypj)
+
+ DO k = 1,3
+!c! location of polar head is computed by taking hydrophobic centre
+!c! and moving by a d1 * dc_norm vector
+!c! see unres publications for very informative images
+ chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres)
+ chead(k,2) = c(k, j)
+!c! distance
+!c! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+!c! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ Rhead_distance(k) = chead(k,2) - chead(k,1)
+ END DO
+!c! 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)))
+!c!-------------------------------------------------------------------
+!c! zero everything that should be zero'ed
+ Egb = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ eheadtail = 0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ RETURN
+ END SUBROUTINE elgrad_init_cat
+
+ SUBROUTINE elgrad_init_cat_pep(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+ use comm_momo
+ use calc_data
+ real(kind=8) :: eheadtail,Egb,Ecl,Elj,Equad,Epol,Rb
+ eps_out=80.0d0
+ itypi = 10
+ itypj = itype(j,5)
+!c! 1/(Gas Constant * Thermostate temperature) = BetaT
+!c! ENABLE THIS LINE WHEN USING CHECKGRAD!!!
+!c! t_bath = 300
+!c! BetaT = 1.0d0 / (t_bath * Rb)i
+ Rb=0.001986d0
+ BetaT = 1.0d0 / (298.0d0 * Rb)
+!c! Gay-berne var's
+ sig0ij = sigmacat( itypi,itypj )
+ chi1 = chi1cat( itypi, itypj )
+ chi2 = 0.0d0
+ chi12 = 0.0d0
+ chip1 = chipp1cat( itypi, itypj )
+ chip2 = 0.0d0
+ chip12 = 0.0d0
+!c! not used by momo potential, but needed by sc_angular which is shared
+!c! by all energy_potential subroutines
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+ dxj = 0.0d0 !dc_norm( 1, nres+j )
+ dyj = 0.0d0 !dc_norm( 2, nres+j )
+ dzj = 0.0d0 !dc_norm( 3, nres+j )
+!c! distance from center of chain(?) to polar/charged head
+ d1 = dheadcat(1, 1, itypi, itypj)
+ d2 = dheadcat(2, 1, itypi, itypj)
+!c! ai*aj from Fgb
+ a12sq = rborn1cat(itypi,itypj) * rborn2cat(itypi,itypj)
+!c! a12sq = a12sq * a12sq
+!c! charge of amino acid itypi is...
+ Qi = 0
+ Qj = ichargecat(itypj)
+! Qij = Qi * Qj
+!c! chis1,2,12
+ chis1 = chis1cat(itypi,itypj)
+ chis2 = 0.0d0
+ chis12 = 0.0d0
+ sig1 = sigmap1cat(itypi,itypj)
+ sig2 = sigmap2cat(itypi,itypj)
+!c! alpha factors from Fcav/Gcav
+ b1cav = alphasurcat(1,itypi,itypj)
+ b2cav = alphasurcat(2,itypi,itypj)
+ b3cav = alphasurcat(3,itypi,itypj)
+ b4cav = alphasurcat(4,itypi,itypj)
+ wqd = wquadcat(itypi, itypj)
+!c! used by Fgb
+ eps_in = epsintabcat(itypi,itypj)
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+!c!-------------------------------------------------------------------
+!c! tail location and distance calculations
+ Rtail = 0.0d0
+ DO k = 1, 3
+ ctail(k,1)=(c(k,i)+c(k,i+1))/2.0-dtailcat(1,itypi,itypj)*dc_norm(k,i)
+ ctail(k,2)=c(k,j)!-dtailcat(2,itypi,itypj)*dc_norm(k,nres+j)
+ END DO
+!c! tail distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+ Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 )
+ Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 )
+ Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 )
+ Rtail = dsqrt( &
+ (Rtail_distance(1)*Rtail_distance(1)) &
+ + (Rtail_distance(2)*Rtail_distance(2)) &
+ + (Rtail_distance(3)*Rtail_distance(3)))
+!c!-------------------------------------------------------------------
+!c! Calculate location and distance between polar heads
+!c! distance between heads
+!c! for each one of our three dimensional space...
+ d1 = dheadcat(1, 1, itypi, itypj)
+ d2 = dheadcat(2, 1, itypi, itypj)
+
+ DO k = 1,3
+!c! location of polar head is computed by taking hydrophobic centre
+!c! and moving by a d1 * dc_norm vector
+!c! see unres publications for very informative images
+ chead(k,1) = (c(k, i)+c(k,i+1))/2.0 + d1 * dc_norm(k, i)
+ chead(k,2) = c(k, j)
+!c! distance
+!c! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+!c! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ Rhead_distance(k) = chead(k,2) - chead(k,1)
+ END DO
+!c! 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)))
+!c!-------------------------------------------------------------------
+!c! zero everything that should be zero'ed
+ Egb = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ eheadtail = 0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ RETURN
+ END SUBROUTINE elgrad_init_cat_pep
+
+ double precision function tschebyshev(m,n,x,y)
+ implicit none
+ integer i,m,n
+ double precision x(n),y,yy(0:maxvar),aux
+!c Tschebyshev polynomial. Note that the first term is omitted
+!c m=0: the constant term is included
+!c m=1: the constant term is not included
+ yy(0)=1.0d0
+ yy(1)=y
+ do i=2,n
+ yy(i)=2*yy(1)*yy(i-1)-yy(i-2)
+ enddo
+ aux=0.0d0
+ do i=m,n
+ aux=aux+x(i)*yy(i)
+ enddo
+ tschebyshev=aux
+ return
+ end function tschebyshev
+!C--------------------------------------------------------------------------
+ double precision function gradtschebyshev(m,n,x,y)
+ implicit none
+ integer i,m,n
+ double precision x(n+1),y,yy(0:maxvar),aux
+!c Tschebyshev polynomial. Note that the first term is omitted
+!c m=0: the constant term is included
+!c m=1: the constant term is not included
+ yy(0)=1.0d0
+ yy(1)=2.0d0*y
+ do i=2,n
+ yy(i)=2*y*yy(i-1)-yy(i-2)
+ enddo
+ aux=0.0d0
+ do i=m,n
+ aux=aux+x(i+1)*yy(i)*(i+1)
+!C print *, x(i+1),yy(i),i
+ enddo
+ gradtschebyshev=aux
+ return
+ end function gradtschebyshev
+!!!!!!!!!--------------------------------------------------------------
+ subroutine lipid_bond(elipbond)
+ real(kind=8) :: elipbond,fac,dist_sub,sumdist
+ real(kind=8), dimension(3):: dist
+ integer(kind=8) :: i,j,k,ibra,ityp,jtyp,ityp1
+ elipbond=0.0d0
+! print *,"before",ilipbond_start,ilipbond_end
+ do i=ilipbond_start,ilipbond_end
+! print *,i,i+1,"i,i+1"
+ ityp=itype(i,4)
+ ityp1=itype(i+1,4)
+! print *,ityp,ityp1,"itype"
+ j=i+1
+ if (ityp.eq.12) ibra=i
+ if ((ityp.eq.ntyp1_molec(4)).or.(ityp1.ge.ntyp1_molec(4)-1)) cycle
+ if (ityp.eq.(ntyp1_molec(4)-1)) then
+ !cofniecie do ostatnie GL1
+! i=ibra
+ j=ibra