+ subroutine lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+ implicit none
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CHAIN'
+ double precision xi,yi,zi,sslipi,ssgradlipi
+ double precision fracinbuf
+! double precision sscalelip,sscagradlip
+#ifdef DEBUG
+ write (iout,*) "bordlipbot",bordlipbot," bordliptop",bordliptop
+ write (iout,*) "buflipbot",buflipbot," lipbufthick",lipbufthick
+ write (iout,*) "xi yi zi",xi,yi,zi
+#endif
+ if ((zi.gt.bordlipbot).and.(zi.lt.bordliptop)) then
+! the energy transfer exist
+ if (zi.lt.buflipbot) then
+! what fraction I am in
+ fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
+! lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
+#ifdef DEBUG
+ write (iout,*) "sslipi",sslipi," ssgradlipi",ssgradlipi
+#endif
+ return
+ end subroutine lipid_layer
+!-------------------------------------------------------------
+ subroutine ecat_prot_transition(ecation_prottran)
+ integer:: itypi,itypj,ityptrani,ityptranj,k,l,i,j
+ real(kind=8),dimension(3):: cjtemp,citemp,diff,dsctemp,vecsc,&
+ diffnorm,boxx,r,dEvan1Cm,dEvan2Cm,dEtotalCm
+ real(kind=8):: ecation_prottran,dista,sdist,De,ene,x0left,&
+ alphac,grad,sumvec,simplesum,pom,erdxi,facd1,&
+ sss_ele_cut,sss_ele_cut_grad,sss2min,sss2mingrad,&
+ ene1,ene2,grad1,grad2,evan1,evan2,rcal,r4,r7,r0p,&
+ r06,r012,epscalc,rocal,ract
+ ecation_prottran=0.0d0
+ boxx(1)=boxxsize
+ boxx(2)=boxysize
+ boxx(3)=boxzsize
+ write(iout,*) "start ecattran",g_listcatsctran_start,g_listcatsctran_end
+ do k=g_listcatsctran_start,g_listcatsctran_end
+ i=newcontlistcatsctrani(k)
+ j=newcontlistcatsctranj(k)
+! print *,i,j,"in new tran"
+ do l=1,3
+ citemp(l)=c(l,i+nres)
+ cjtemp(l)=c(l,j)
+ enddo
+
+ itypi=itype(i,1) !as the first is the protein part
+ itypj=itype(j,5) !as the second part is always cation
+! remapping to internal types
+! read (iiontran,*,err=123,end=123) (agamacattran(k,j,i),k=1,3),&
+! (athetacattran(k,j,i),k=1,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),&
+! demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),&
+! x0cattrans(j,i)
+
+ if (itypj.eq.6) then
+ ityptranj=1 !as now only Zn2+ is this needs to be modified for other ions
+ endif
+ if (itypi.eq.16) then
+ ityptrani=1
+ elseif (itypi.eq.1) then
+ ityptrani=2
+ elseif (itypi.eq.15) then
+ ityptrani=3
+ elseif (itypi.eq.17) then
+ ityptrani=4
+ elseif (itypi.eq.2) then
+ ityptrani=5
+ else
+ ityptrani=6
+ endif
+
+ if (ityptrani.gt.ntrantyp(ityptranj)) then
+! do l=1,3
+! write(iout,*),gradcattranc(l,j),gradcattranx(l,i)
+! enddo
+!volume excluded
+ call to_box(cjtemp(1),cjtemp(2),cjtemp(3))
+ call to_box(citemp(1),citemp(2),citemp(3))
+ rcal=0.0d0
+ do l=1,3
+ r(l)=boxshift(cjtemp(l)-citemp(l),boxx(l))
+ rcal=rcal+r(l)*r(l)
+ enddo
+ ract=sqrt(rcal)
+ if (ract.gt.r_cut_ele) cycle
+ sss_ele_cut=sscale_ele(ract)
+ sss_ele_cut_grad=sscagrad_ele(ract)
+ 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 l=1,3
+ dEvan1Cm(l) = 12*r(l)*epscalc*r012/r7
+ dEvan2Cm(l) = 12*r(l)*epscalc*r06/r4
+ enddo
+ do l=1,3
+ dEtotalCm(l)=(dEvan1Cm(l)+dEvan2Cm(l))*sss_ele_cut-&
+ (Evan1+Evan2)*sss_ele_cut_grad*r(l)/ract
+ enddo
+ ecation_prottran = ecation_prottran+&
+ (Evan1+Evan2)*sss_ele_cut
+ do l=1,3
+ gradcattranx(l,i)=gradcattranx(l,i)+dEtotalCm(l)
+ gradcattranc(l,i)=gradcattranc(l,i)+dEtotalCm(l)
+ gradcattranc(l,j)=gradcattranc(l,j)-dEtotalCm(l)
+ enddo
+
+ ene=0.0d0
+ else
+! cycle
+ sumvec=0.0d0
+ simplesum=0.0d0
+ do l=1,3
+ vecsc(l)=citemp(l)-c(l,i)
+ sumvec=sumvec+vecsc(l)**2
+ simplesum=simplesum+vecsc(l)
+ enddo
+ sumvec=dsqrt(sumvec)
+ call to_box(cjtemp(1),cjtemp(2),cjtemp(3))
+ call to_box(citemp(1),citemp(2),citemp(3))
+! sumvec=2.0d0
+ do l=1,3
+ dsctemp(l)=c(l,i+nres)&
+ +(acatshiftdsc(ityptrani,ityptranj)-1.0d0)*vecsc(l)&
+ +bcatshiftdsc(ityptrani,ityptranj)*vecsc(l)/sumvec
+ enddo
+ call to_box(dsctemp(1),dsctemp(2),dsctemp(3))
+ sdist=0.0d0
+ do l=1,3
+ diff(l)=boxshift(dsctemp(l)-cjtemp(l),boxx(l))
+ sdist=sdist+diff(l)*diff(l)
+ enddo
+ dista=sqrt(sdist)
+ if (dista.gt.r_cut_ele) cycle
+
+ sss_ele_cut=sscale_ele(dista)
+ sss_ele_cut_grad=sscagrad_ele(dista)
+ sss2min=sscale2(dista,x0cattrans(ityptrani,ityptranj)-0.1d0,0.2d0)
+ De=demorsecat(ityptrani,ityptranj)
+ alphac=alphamorsecat(ityptrani,ityptranj)
+ if (sss2min.eq.1.0d0) then
+! print *,"ityptrani",ityptrani,ityptranj
+ x0left=x0catleft(ityptrani,ityptranj) ! to mn
+ ene=sss_ele_cut*(-De+De*(1.0d0-dexp(-alphac*(dista-x0left)))**2)
+ grad=2.0d0*alphac*De*dexp(-alphac*(dista-x0left))*&
+ (1.0d0-dexp(-alphac*(dista-x0left)))*sss_ele_cut&
+ +ene/sss_ele_cut*sss_ele_cut_grad
+ else if (sss2min.eq.0.0d0) then
+ x0left=x0catright(ityptrani,ityptranj)
+ ene=sss_ele_cut*(-De+De*(1.0d0-dexp(-alphac*(dista-x0left)))**2)
+ grad=2.0d0*alphac*De*dexp(-alphac*(dista-x0left))*&
+ (1.0d0-dexp(-alphac*(dista-x0left)))*sss_ele_cut&
+ +ene/sss_ele_cut*sss_ele_cut_grad
+ else
+ sss2mingrad=sscagrad2(dista,x0cattrans(ityptrani,ityptranj)-0.1d0,0.2d0)
+ x0left=x0catleft(ityptrani,ityptranj)
+ ene1=sss_ele_cut*(-De+De*(1.0d0-dexp(-alphac*(dista-x0left)))**2)
+ grad1=2.0d0*alphac*De*dexp(-alphac*(dista-x0left))*&
+ (1.0d0-dexp(-alphac*(dista-x0left)))*sss_ele_cut&
+ +ene/sss_ele_cut*sss_ele_cut_grad
+ x0left=x0catright(ityptrani,ityptranj)
+ ene2=sss_ele_cut*(-De+De*(1.0d0-dexp(-alphac*(dista-x0left)))**2)
+ grad2=2.0d0*alphac*De*dexp(-alphac*(dista-x0left))*&
+ (1.0d0-dexp(-alphac*(dista-x0left)))*sss_ele_cut&
+ +ene/sss_ele_cut*sss_ele_cut_grad
+ ene=sss2min*ene1+(1.0d0-sss2min)*ene2
+ grad=sss2min*grad1+(1.0d0-sss2min)*grad2+sss2mingrad*(ene1-ene2)
+ endif
+ do l=1,3
+ diffnorm(l)= diff(l)/dista
+ enddo
+ erdxi=scalar(diffnorm(1),dc_norm(1,i+nres))
+ facd1=bcatshiftdsc(ityptrani,ityptranj)/sumvec
+
+ do l=1,3
+! 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)
+ pom=diffnorm(l)+facd1*(diffnorm(l)-erdxi*dc_norm(l,i+nres))
+! write(iout,*),gradcattranc(l,j),gradcattranx(l,i),grad*diff(l)/dista
+
+ gradcattranx(l,i)=gradcattranx(l,i)+grad*pom&
+ +grad*diffnorm(l)*(acatshiftdsc(ityptrani,ityptranj)-1.0d0)
+! *( bcatshiftdsc(ityptrani,ityptranj)*&
+! (1.0d0/sumvec-(vecsc(l)*simplesum)*(sumvec**(-3.0d0))))
+ gradcattranc(l,i)=gradcattranc(l,i)+grad*diff(l)/dista
+! +sss_ele_cut_grad*ene*diff(l)/dista/sss_ele_cut
+ gradcattranc(l,j)=gradcattranc(l,j)-grad*diff(l)/dista
+! -sss_ele_cut_grad*ene*diff(l)/dista/sss_ele_cut
+ enddo
+ ecation_prottran=ecation_prottran+ene
+ if (energy_dec) write(iout,*) "etrancat",i,j,ene,x0left,De,dista,&
+ alphac
+ endif
+ enddo
+! do k=g_listcatptran_start,g_listcatptran_end
+! ene=0.0d0 this will be used if peptide group interaction is needed
+! enddo
+
+
+ return
+ end subroutine
+ subroutine ecat_prot_ang(ecation_protang)
+ integer:: itypi,itypj,ityptrani,ityptranj,k,l,i,j,n,m,&
+ ityptrani1,ityptranj1,ityptrani2,ityptranj2,&
+ i1,i2,j1,j2,k1,k2,k3,i3,j3,ityptrani3,ityptranj3
+
+ real(kind=8),dimension(3):: cjtemp,citemp,diff,dsctemp,vecsc,&
+ diffnorm,boxx,dscvec,dscvecnorm,diffnorm2,&
+ dscvec2,dscvecnorm2,cjtemp2,citemp2,diff2,dsctemp2,&
+ vecsc2,diff1,diffnorm1,diff3,mindiffnorm2
+ real(kind=8),dimension(3):: dscvec1,dscvecnorm1,cjtemp1,citemp1,vecsc1,dsctemp1,&
+ dscvec3,dscvecnorm3,cjtemp3,citemp3,vecsc3,dsctemp3,&
+ diffnorm3,diff4,diffnorm4
+
+ real(kind=8):: ecation_protang,dista,sdist,De,ene,x0left,&
+ alphac,grad,sumvec,sumdscvec,pom,erdxi,facd1,&
+ sss_ele_cut,sss_ele_cut_grad,sss2min,sss2mingrad,&
+ simplesum,cosval,part1,part2a,part2,part2b,part3,&
+ part4a,part4b,part4,bottom,dista2,sdist2,sumvec2,&
+ sumdscvec2,simplesum2,dista1,sdist1,sumvec1,simplesum1,&
+ sumdscvec1,facd2,scal1a,scal1b,scal2a,scal2b,&
+ sss2mingrad1,sss2mingrad2,sss2min1,sss2min2,pom1,pom2,&
+ det1ij,det2ij,cosom1,cosom2,cosom12,cosphij,dista3,&
+ sumvec3
+ real(kind=8):: sinom1,sinom2,sinaux,dephiij,sumdscvec3,sumscvec3,&
+ cosphi,sdist3,simplesum3,det1t2ij,sss2mingrad3,sss2min3,&
+ scal1c,scal2c,scal3a,scal3b,scal3c,facd3,facd2b,scal3d,&
+ scal3e,dista4,sdist4,pom3,sssmintot
+
+ ecation_protang=0.0d0
+ boxx(1)=boxxsize
+ boxx(2)=boxysize
+ boxx(3)=boxzsize
+! print *,"KUR**3",g_listcatscang_start,g_listcatscang_end
+! go to 19
+! go to 21
+ do k=g_listcatscang_start,g_listcatscang_end
+ ene=0.0d0
+ i=newcontlistcatscangi(k)
+ j=newcontlistcatscangj(k)
+ itypi=itype(i,1) !as the first is the protein part
+ itypj=itype(j,5) !as the second part is always cation
+! print *,"KUR**4",i,j,itypi,itypj
+! remapping to internal types
+! read (iiontran,*,err=123,end=123) (agamacattran(k,j,i),k=1,3),&
+! (athetacattran(k,j,i),k=1,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),&
+! demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),&
+! x0cattrans(j,i)
+ if (itypj.eq.6) then
+ ityptranj=1 !as now only Zn2+ is this needs to be modified for other ions
+ endif
+ if (itypi.eq.16) then
+ ityptrani=1
+ elseif (itypi.eq.1) then
+ ityptrani=2
+ elseif (itypi.eq.15) then
+ ityptrani=3
+ elseif (itypi.eq.17) then
+ ityptrani=4
+ elseif (itypi.eq.2) then
+ ityptrani=5
+ else
+ ityptrani=6
+ endif
+ if (ityptrani.gt.ntrantyp(ityptranj)) cycle
+ do l=1,3
+ citemp(l)=c(l,i+nres)
+ cjtemp(l)=c(l,j)
+ enddo
+ sumvec=0.0d0
+ simplesum=0.0d0
+ do l=1,3
+ vecsc(l)=citemp(l)-c(l,i)
+ sumvec=sumvec+vecsc(l)**2
+ simplesum=simplesum+vecsc(l)
+ enddo
+ sumvec=dsqrt(sumvec)
+ sumdscvec=0.0d0
+ do l=1,3
+ dsctemp(l)=c(l,i)&
+! +1.0d0
+ +(acatshiftdsc(ityptrani,ityptranj))*vecsc(l)&
+ +bcatshiftdsc(ityptrani,ityptranj)*vecsc(l)/sumvec
+ dscvec(l)= &
+!1.0d0
+ (acatshiftdsc(ityptrani,ityptranj))*vecsc(l)&
+ +bcatshiftdsc(ityptrani,ityptranj)*vecsc(l)/sumvec
+ sumdscvec=sumdscvec+dscvec(l)**2
+ enddo
+ sumdscvec=dsqrt(sumdscvec)
+ do l=1,3
+ dscvecnorm(l)=dscvec(l)/sumdscvec
+ enddo
+ call to_box(dsctemp(1),dsctemp(2),dsctemp(3))
+ call to_box(cjtemp(1),cjtemp(2),cjtemp(3))
+ sdist=0.0d0
+ do l=1,3
+ diff(l)=boxshift(dsctemp(l)-cjtemp(l),boxx(l))
+ sdist=sdist+diff(l)*diff(l)
+ enddo
+ dista=sqrt(sdist)
+ do l=1,3
+ diffnorm(l)= diff(l)/dista
+ enddo
+ cosval=scalar(diffnorm(1),dc_norm(1,i+nres))
+ grad=0.0d0
+ sss2min=sscale2(dista,r_cut_ang,1.0d0)
+ sss2mingrad=sscagrad2(dista,r_cut_ang,1.0d0)
+ ene=ene&
+ +tschebyshev(1,6,athetacattran(1,ityptrani,ityptranj),cosval)
+ grad=gradtschebyshev(0,5,athetacattran(1,ityptrani,ityptranj),cosval)*sss2min
+
+ facd1=bcatshiftdsc(ityptrani,ityptranj)/sumvec
+ erdxi=scalar(diffnorm(1),dc_norm(1,i+nres))
+ part1=0.0d0
+ part2=0.0d0
+ part3=0.0d0
+ part4=0.0d0
+ do l=1,3
+ bottom=sumvec**2*sdist
+ part1=diff(l)*sumvec*dista
+ part2a=(acatshiftdsc(ityptrani,ityptranj))*vecsc(l)
+ part2b=0.0d0
+ !bcatshiftdsc(ityptrani,ityptranj)/sumvec*&
+ !(vecsc(l)-cosval*dista*dc_norm(l,i+nres))
+ part2=(part2a+part2b)*sumvec*dista
+ part3=cosval*sumvec*dista*dc_norm(l,i+nres)*dista
+ part4a=diff(l)*acatshiftdsc(ityptrani,ityptranj)
+ part4b=bcatshiftdsc(ityptrani,ityptranj)/sumvec*&
+ (diff(l)-cosval*dista*dc_norm(l,i+nres))
+ part4=cosval*sumvec*(part4a+part4b)*sumvec
+! gradlipang(m,l)=gradlipang(m,l)+(fac &
+! *(xa(m)-scalar*vnorm*xb(m)/wnorm)&
+! /(vnorm*wnorm))
+
+! 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)
+ pom=diffnorm(l)+facd1*(diffnorm(l)-erdxi*dc_norm(l,i+nres))
+
+ gradcatangc(l,j)=gradcatangc(l,j)-grad*&
+ (dscvec(l)-cosval*diffnorm(l)*sumdscvec)/(sumdscvec*dista)-&
+ ene*sss2mingrad*diffnorm(l)
+
+ gradcatangc(l,i)=gradcatangc(l,i)+grad*&
+ (dscvec(l)-cosval*diffnorm(l)*sumdscvec)/(sumdscvec*dista)+&
+ ene*sss2mingrad*diffnorm(l)
+
+ gradcatangx(l,i)=gradcatangx(l,i)+grad*&
+ (part1+part2-part3-part4)/bottom+&
+ ene*sss2mingrad*pom+&
+ ene*sss2mingrad*diffnorm(l)*(acatshiftdsc(ityptrani,ityptranj)-1.0d0)
+! +grad*(dscvec(l)-cosval*diffnorm(l)*sumdscvec)/(sumdscvec*dista)&
+! +grad*pom+grad*diffnorm(l)*(acatshiftdsc(ityptrani,ityptranj)-1.0d0)
+!&
+! (diff(l)-cosval*dscvecnorm(l)*dista)/(sumdscvec*dista)
+
+
+
+
+
+ enddo
+! print *,i,j,cosval,tschebyshev(1,3,aomicattr(1,ityptranj),cosval)&
+! ,aomicattr(0,ityptranj),ene
+ if (energy_dec) write(iout,*) i,j,ityptrani,ityptranj,ene,cosval
+ ecation_protang=ecation_protang+ene*sss2min
+ enddo
+ 19 continue
+! print *,"KUR**",g_listcatscangf_start,g_listcatscangf_end
+ do k=g_listcatscangf_start,g_listcatscangf_end
+ ene=0.0d0
+ i1=newcontlistcatscangfi(k)
+ j1=newcontlistcatscangfj(k)
+ itypi=itype(i1,1) !as the first is the protein part
+ itypj=itype(j1,5) !as the second part is always cation
+ if (itypj.eq.6) then
+ ityptranj1=1 !as now only Zn2+ is this needs to be modified for other ions
+ endif
+ if (itypi.eq.16) then
+ ityptrani1=1
+ elseif (itypi.eq.1) then
+ ityptrani1=2
+ elseif (itypi.eq.15) then
+ ityptrani1=3
+ elseif (itypi.eq.17) then
+ ityptrani1=4
+ elseif (itypi.eq.2) then
+ ityptrani1=5
+ else
+ ityptrani1=6
+ endif
+ do l=1,3
+ citemp1(l)=c(l,i1+nres)
+ cjtemp1(l)=c(l,j1)
+ enddo
+ sumvec1=0.0d0
+ simplesum1=0.0d0
+ do l=1,3
+ vecsc1(l)=citemp1(l)-c(l,i1)
+ sumvec1=sumvec1+vecsc1(l)**2
+ simplesum1=simplesum1+vecsc1(l)
+ enddo
+ sumvec1=dsqrt(sumvec1)
+ sumdscvec1=0.0d0
+ do l=1,3
+ dsctemp1(l)=c(l,i1)&
+! +1.0d0
+ +(acatshiftdsc(ityptrani1,ityptranj1))*vecsc1(l)&
+ +bcatshiftdsc(ityptrani1,ityptranj1)*vecsc1(l)/sumvec1
+ dscvec1(l)= &
+!1.0d0
+ (acatshiftdsc(ityptrani1,ityptranj1))*vecsc1(l)&
+ +bcatshiftdsc(ityptrani1,ityptranj1)*vecsc1(l)/sumvec1
+ sumdscvec1=sumdscvec1+dscvec1(l)**2
+ enddo
+ sumdscvec1=dsqrt(sumdscvec1)
+ do l=1,3
+ dscvecnorm1(l)=dscvec1(l)/sumdscvec1
+ enddo
+ call to_box(dsctemp1(1),dsctemp1(2),dsctemp1(3))
+ call to_box(cjtemp1(1),cjtemp1(2),cjtemp1(3))
+ sdist1=0.0d0
+ do l=1,3
+ diff1(l)=boxshift(dsctemp1(l)-cjtemp1(l),boxx(l))
+ sdist1=sdist1+diff1(l)*diff1(l)
+ enddo
+ dista1=sqrt(sdist1)
+ do l=1,3
+ diffnorm1(l)= diff1(l)/dista1
+ enddo
+ sss2min1=sscale2(dista1,r_cut_ang,1.0d0)
+ sss2mingrad1=sscagrad2(dista1,r_cut_ang,1.0d0)
+ if (ityptrani1.gt.ntrantyp(ityptranj1)) cycle
+
+!-----------------------------------------------------------------
+! do m=k+1,g_listcatscang_end
+ ene=0.0d0
+ i2=newcontlistcatscangfk(k)
+ j2=j1
+ if (j1.ne.j2) cycle
+ itypi=itype(i2,1) !as the first is the protein part
+ itypj=itype(j2,5) !as the second part is always cation
+ if (itypj.eq.6) then
+ ityptranj2=1 !as now only Zn2+ is this needs to be modified for other ions
+ endif
+ if (itypi.eq.16) then
+ ityptrani2=1
+ elseif (itypi.eq.1) then
+ ityptrani2=2
+ elseif (itypi.eq.15) then
+ ityptrani2=3
+ elseif (itypi.eq.17) then
+ ityptrani2=4
+ elseif (itypi.eq.2) then
+ ityptrani2=5
+ else
+ ityptrani2=6
+ endif
+ if (ityptrani2.gt.ntrantyp(ityptranj2)) cycle
+
+ do l=1,3
+ citemp2(l)=c(l,i2+nres)
+ cjtemp2(l)=c(l,j2)
+ enddo
+ sumvec2=0.0d0
+ simplesum2=0.0d0
+ do l=1,3
+ vecsc2(l)=citemp2(l)-c(l,i2)
+ sumvec2=sumvec2+vecsc2(l)**2
+ simplesum2=simplesum2+vecsc2(l)
+ enddo
+ sumvec2=dsqrt(sumvec2)
+ sumdscvec2=0.0d0
+ do l=1,3
+ dsctemp2(l)=c(l,i2)&
+! +1.0d0
+ +(acatshiftdsc(ityptrani2,ityptranj2))*vecsc2(l)&
+ +bcatshiftdsc(ityptrani2,ityptranj2)*vecsc2(l)/sumvec2
+ dscvec2(l)= &
+!1.0d0
+ (acatshiftdsc(ityptrani2,ityptranj2))*vecsc2(l)&
+ +bcatshiftdsc(ityptrani2,ityptranj2)*vecsc2(l)/sumvec2
+ sumdscvec2=sumdscvec2+dscvec2(l)**2
+ enddo
+ sumdscvec2=dsqrt(sumdscvec2)
+ do l=1,3
+ dscvecnorm2(l)=dscvec2(l)/sumdscvec2
+ enddo
+ call to_box(dsctemp2(1),dsctemp2(2),dsctemp2(3))
+ call to_box(cjtemp2(1),cjtemp2(2),cjtemp2(3))
+ sdist2=0.0d0
+ do l=1,3
+ diff2(l)=boxshift(dsctemp2(l)-cjtemp2(l),boxx(l))
+! diff2(l)=1.0d0
+ sdist2=sdist2+diff2(l)*diff2(l)
+ enddo
+ dista2=sqrt(sdist2)
+ do l=1,3
+ diffnorm2(l)= diff2(l)/dista2
+ enddo
+! print *,i1,i2,diffnorm2(1)
+ cosval=scalar(diffnorm1(1),diffnorm2(1))
+ grad=0.0d0
+ sss2min2=sscale2(dista2,r_cut_ang,1.0d0)
+ sss2mingrad2=sscagrad2(dista2,r_cut_ang,1.0d0)
+ ene=ene+tschebyshev(1,3,aomicattr(1,ityptranj1),cosval)
+ grad=gradtschebyshev(0,2,aomicattr(1,ityptranj1),cosval)*sss2min2*sss2min1
+ part1=0.0d0
+ part2=0.0d0
+ part3=0.0d0
+ part4=0.0d0
+ ecation_protang=ecation_protang+ene*sss2min2*sss2min1
+ facd1=bcatshiftdsc(ityptrani1,ityptranj1)/sumvec1
+ facd2=bcatshiftdsc(ityptrani2,ityptranj2)/sumvec2
+ scal1a=scalar(diffnorm1(1),dc_norm(1,i1+nres))
+ scal1b=scalar(diffnorm2(1),dc_norm(1,i1+nres))
+ scal2a=scalar(diffnorm1(1),dc_norm(1,i2+nres))
+ scal2b=scalar(diffnorm2(1),dc_norm(1,i2+nres))
+
+ if (energy_dec) write(iout,*) "omi", i,j,ityptrani,ityptranj,ene,cosval,aomicattr(1,ityptranj1),&
+ aomicattr(2,ityptranj1),aomicattr(3,ityptranj1),tschebyshev(1,3,aomicattr(1,ityptranj1),cosval)
+
+!*sss2min
+ do l=1,3
+ pom1=diffnorm1(l)+facd1*(diffnorm1(l)-scal1a*dc_norm(l,i1+nres))
+ pom2=diffnorm2(l)+facd2*(diffnorm2(l)-scal2b*dc_norm(l,i2+nres))
+
+
+ gradcatangc(l,i1)=gradcatangc(l,i1)+grad*(diff2(l)-&
+ cosval*diffnorm1(l)*dista2)/(dista2*dista1)+&
+ ene*sss2mingrad1*diffnorm1(l)*sss2min2
+
+
+ gradcatangx(l,i1)=gradcatangx(l,i1)+grad/(dista2*dista1)*&
+ (acatshiftdsc(ityptrani1,ityptranj1)*diff2(l)+&
+ facd1*(diff2(l)-scal1b*dc_norm(l,i1+nres)*dista2)-&
+ cosval*dista2/dista1*&
+ (acatshiftdsc(ityptrani1,ityptranj1)*diff1(l)+&
+ facd1*(diff1(l)-scal1a*dc_norm(l,i1+nres)*dista1)))+&
+ ene*sss2mingrad1*sss2min2*(pom1+&
+ diffnorm1(l)*(acatshiftdsc(ityptrani1,ityptranj1)-1.0d0))
+
+
+ gradcatangx(l,i2)=gradcatangx(l,i2)+grad/(dista2*dista1)*&
+ (acatshiftdsc(ityptrani2,ityptranj2)*diff1(l)+&
+ facd2*(diff1(l)-scal2a*dc_norm(l,i2+nres)*dista1)-&
+ cosval*dista1/dista2*&
+ (acatshiftdsc(ityptrani2,ityptranj2)*diff2(l)+&
+ facd2*(diff2(l)-scal2b*dc_norm(l,i2+nres)*dista2)))+&
+ ene*sss2mingrad2*sss2min1*(pom2+&
+ diffnorm2(l)*(acatshiftdsc(ityptrani2,ityptranj2)-1.0d0))
+
+
+ gradcatangx(l,i2)=gradcatangx(l,i2)
+ gradcatangc(l,i2)=gradcatangc(l,i2)+grad*(diff1(l)-&
+ cosval*diffnorm2(l)*dista1)/(dista2*dista1)+&
+ ene*sss2mingrad2*diffnorm2(l)*sss2min1
+
+ gradcatangc(l,j2)=gradcatangc(l,j2)-grad*(diff2(l)/dista2/dista1-&
+ cosval*diff1(l)/dista1/dista1+diff1(l)/dista2/dista1-&
+ cosval*diff2(l)/dista2/dista2)-&
+ ene*sss2mingrad1*diffnorm1(l)*sss2min2-&
+ ene*sss2mingrad2*diffnorm2(l)*sss2min1
+
+
+ enddo
+
+ enddo
+! enddo
+!#ifdef DUBUG
+ 21 continue
+! do k1=g_listcatscang_start,g_listcatscang_end
+! print *,"KURNA",g_listcatscangt_start,g_listcatscangt_end
+ do k1=g_listcatscangt_start,g_listcatscangt_end
+ i1=newcontlistcatscangti(k1)
+ j1=newcontlistcatscangtj(k1)
+ itypi=itype(i1,1) !as the first is the protein part
+ itypj=itype(j1,5) !as the second part is always cation
+ if (itypj.eq.6) then
+ ityptranj1=1 !as now only Zn2+ is this needs to be modified for other ions
+ endif
+ if (itypi.eq.16) then
+ ityptrani1=1
+ elseif (itypi.eq.1) then
+ ityptrani1=2
+ elseif (itypi.eq.15) then
+ ityptrani1=3
+ elseif (itypi.eq.17) then
+ ityptrani1=4
+ elseif (itypi.eq.2) then
+ ityptrani1=5
+ else
+ ityptrani1=6
+ endif
+ do l=1,3
+ citemp1(l)=c(l,i1+nres)
+ cjtemp1(l)=c(l,j1)
+ enddo
+ sumvec1=0.0d0
+ simplesum1=0.0d0
+ do l=1,3
+ vecsc1(l)=citemp1(l)-c(l,i1)
+ sumvec1=sumvec1+vecsc1(l)**2
+ simplesum1=simplesum1+vecsc1(l)
+ enddo
+ sumvec1=dsqrt(sumvec1)
+ sumdscvec1=0.0d0
+ do l=1,3
+ dsctemp1(l)=c(l,i1)&
+ +(acatshiftdsc(ityptrani1,ityptranj1))*vecsc1(l)&
+ +bcatshiftdsc(ityptrani1,ityptranj1)*vecsc1(l)/sumvec1
+ dscvec1(l)= &
+ (acatshiftdsc(ityptrani1,ityptranj1))*vecsc1(l)&
+ +bcatshiftdsc(ityptrani1,ityptranj1)*vecsc1(l)/sumvec1
+ sumdscvec1=sumdscvec1+dscvec1(l)**2
+ enddo
+ sumdscvec1=dsqrt(sumdscvec1)
+ do l=1,3
+ dscvecnorm1(l)=dscvec1(l)/sumdscvec1
+ enddo
+ call to_box(dsctemp1(1),dsctemp1(2),dsctemp1(3))
+ call to_box(cjtemp1(1),cjtemp1(2),cjtemp1(3))
+ sdist1=0.0d0
+ do l=1,3
+ diff1(l)=boxshift(dsctemp1(l)-cjtemp1(l),boxx(l))
+ sdist1=sdist1+diff1(l)*diff1(l)
+ enddo
+ dista1=sqrt(sdist1)
+ do l=1,3
+ diffnorm1(l)= diff1(l)/dista1
+ enddo
+ sss2min1=sscale2(dista1,r_cut_ang,1.0d0)
+ sss2mingrad1=sscagrad2(dista1,r_cut_ang,1.0d0)
+ if (ityptrani1.gt.ntrantyp(ityptranj1)) cycle
+!---------------before second loop
+! do k2=k1+1,g_listcatscang_end
+ i2=newcontlistcatscangtk(k1)
+ j2=j1
+! print *,"TUTU3",i1,i2,j1,j2
+ if (i2.eq.i1) cycle
+ if (j2.ne.j1) cycle
+ itypi=itype(i2,1) !as the first is the protein part
+ itypj=itype(j2,5) !as the second part is always cation
+ if (itypj.eq.6) then
+ ityptranj2=1 !as now only Zn2+ is this needs to be modified for other ions
+ endif
+ if (itypi.eq.16) then
+ ityptrani2=1
+ elseif (itypi.eq.1) then
+ ityptrani2=2
+ elseif (itypi.eq.15) then
+ ityptrani2=3
+ elseif (itypi.eq.17) then
+ ityptrani2=4
+ elseif (itypi.eq.2) then
+ ityptrani2=5
+ else
+ ityptrani2=6
+ endif
+ if (ityptrani2.gt.ntrantyp(ityptranj2)) cycle
+ do l=1,3
+ citemp2(l)=c(l,i2+nres)
+ cjtemp2(l)=c(l,j2)
+ enddo
+ sumvec2=0.0d0
+ simplesum2=0.0d0
+ do l=1,3
+ vecsc2(l)=citemp2(l)-c(l,i2)
+ sumvec2=sumvec2+vecsc2(l)**2
+ simplesum2=simplesum2+vecsc2(l)
+ enddo
+ sumvec2=dsqrt(sumvec2)
+ sumdscvec2=0.0d0
+ do l=1,3
+ dsctemp2(l)=c(l,i2)&
+ +(acatshiftdsc(ityptrani2,ityptranj2))*vecsc2(l)&
+ +bcatshiftdsc(ityptrani2,ityptranj2)*vecsc2(l)/sumvec2
+ dscvec2(l)= &
+ (acatshiftdsc(ityptrani2,ityptranj2))*vecsc2(l)&
+ +bcatshiftdsc(ityptrani2,ityptranj2)*vecsc2(l)/sumvec2
+ sumdscvec2=sumdscvec2+dscvec2(l)**2
+ enddo
+ sumdscvec2=dsqrt(sumdscvec2)
+ do l=1,3
+ dscvecnorm2(l)=dscvec2(l)/sumdscvec2
+ enddo
+ call to_box(dsctemp2(1),dsctemp2(2),dsctemp2(3))
+ call to_box(cjtemp2(1),cjtemp2(2),cjtemp2(3))
+ sdist2=0.0d0
+ do l=1,3
+ diff2(l)=boxshift(dsctemp2(l)-cjtemp2(l),boxx(l))
+! diff2(l)=1.0d0
+ sdist2=sdist2+diff2(l)*diff2(l)
+ enddo
+ dista2=sqrt(sdist2)
+ do l=1,3
+ diffnorm2(l)= diff2(l)/dista2
+ mindiffnorm2(l)=-diffnorm2(l)
+ enddo
+! print *,i1,i2,diffnorm2(1)
+ cosom1=scalar(diffnorm1(1),diffnorm2(1))
+ sss2min2=sscale2(dista2,r_cut_ang,1.0d0)
+ sss2mingrad2=sscagrad2(dista2,r_cut_ang,1.0d0)
+
+!---------------- before third loop
+! do k3=g_listcatscang_start,g_listcatscang_end
+ ene=0.0d0
+ i3=newcontlistcatscangtl(k1)
+ j3=j1
+! print *,"TUTU4",i1,i2,i3,j1,j2,j3
+
+ if (i3.eq.i2) cycle
+ if (i3.eq.i1) cycle
+ if (j3.ne.j1) cycle
+ itypi=itype(i3,1) !as the first is the protein part
+ itypj=itype(j3,5) !as the second part is always cation
+ if (itypj.eq.6) then
+ ityptranj3=1 !as now only Zn2+ is this needs to be modified for other ions
+ endif
+ if (itypi.eq.16) then
+ ityptrani3=1
+ elseif (itypi.eq.1) then
+ ityptrani3=2
+ elseif (itypi.eq.15) then
+ ityptrani3=3
+ elseif (itypi.eq.17) then
+ ityptrani3=4
+ elseif (itypi.eq.2) then
+ ityptrani3=5
+ else
+ ityptrani3=6
+ endif
+ if (ityptrani3.gt.ntrantyp(ityptranj3)) cycle
+ do l=1,3
+ citemp3(l)=c(l,i3+nres)
+ cjtemp3(l)=c(l,j3)
+ enddo
+ sumvec3=0.0d0
+ simplesum3=0.0d0
+ do l=1,3
+ vecsc3(l)=citemp3(l)-c(l,i3)
+ sumvec3=sumvec3+vecsc3(l)**2
+ simplesum3=simplesum3+vecsc3(l)
+ enddo
+ sumvec3=dsqrt(sumvec3)
+ sumdscvec3=0.0d0
+ do l=1,3
+ dsctemp3(l)=c(l,i3)&
+ +(acatshiftdsc(ityptrani3,ityptranj3))*vecsc3(l)&
+ +bcatshiftdsc(ityptrani3,ityptranj3)*vecsc3(l)/sumvec3
+ dscvec3(l)= &
+ (acatshiftdsc(ityptrani3,ityptranj3))*vecsc3(l)&
+ +bcatshiftdsc(ityptrani3,ityptranj3)*vecsc3(l)/sumvec3
+ sumdscvec3=sumdscvec3+dscvec3(l)**2
+ enddo
+ sumdscvec3=dsqrt(sumdscvec3)
+ do l=1,3
+ dscvecnorm3(l)=dscvec3(l)/sumdscvec3
+ enddo
+ call to_box(dsctemp3(1),dsctemp3(2),dsctemp3(3))
+ call to_box(cjtemp3(1),cjtemp3(2),cjtemp3(3))
+ sdist3=0.0d0
+ do l=1,3
+ diff3(l)=boxshift(dsctemp3(l)-dsctemp2(l),boxx(l))
+ sdist3=sdist3+diff3(l)*diff3(l)
+ enddo
+ dista3=sqrt(sdist3)
+ do l=1,3
+ diffnorm3(l)= diff3(l)/dista3
+ enddo
+ sdist4=0.0d0
+ do l=1,3
+ diff4(l)=boxshift(dsctemp3(l)-cjtemp2(l),boxx(l))
+! diff2(l)=1.0d0
+ sdist4=sdist4+diff4(l)*diff4(l)
+ enddo
+ dista4=sqrt(sdist4)
+ do l=1,3
+ diffnorm4(l)= diff4(l)/dista4
+ enddo
+
+ sss2min3=sscale2(dista4,r_cut_ang,1.0d0)
+ sss2mingrad3=sscagrad2(dista4,r_cut_ang,1.0d0)
+ sssmintot=sss2min3*sss2min2*sss2min1
+ if (ityptrani3.gt.ntrantyp(ityptranj3)) cycle
+ cosom12=scalar(diffnorm3(1),diffnorm1(1))
+ cosom2=scalar(diffnorm3(1),mindiffnorm2(1))
+ sinom1=dsqrt(1.0d0-cosom1*cosom1)
+ sinom2=dsqrt(1.0d0-cosom2*cosom2)
+ cosphi=cosom12-cosom1*cosom2
+ sinaux=sinom1*sinom2
+ ene=ene+mytschebyshev(1,3,agamacattran(1,ityptrani2,ityptranj2),cosphi,sinaux)
+ call mygradtschebyshev(1,3,agamacattran(1,ityptrani2,ityptranj2)&
+ ,cosphi,sinaux,dephiij,det1t2ij)
+
+ det1ij=-det1t2ij*sinom2*cosom1/sinom1-dephiij*cosom2
+ det2ij=-det1t2ij*sinom1*cosom2/sinom2-dephiij*cosom1
+ facd1=bcatshiftdsc(ityptrani1,ityptranj1)/sumvec1
+ facd2=bcatshiftdsc(ityptrani2,ityptranj2)/sumvec2
+! facd2b=bcatshiftdsc(ityptrani2,ityptranj2)/sumvec3
+ facd3=bcatshiftdsc(ityptrani3,ityptranj3)/sumvec3
+ scal1a=scalar(diffnorm1(1),dc_norm(1,i1+nres))
+ scal1b=scalar(diffnorm2(1),dc_norm(1,i1+nres))
+ scal1c=scalar(diffnorm3(1),dc_norm(1,i1+nres))
+ scal2a=scalar(diffnorm1(1),dc_norm(1,i2+nres))
+ scal2b=scalar(diffnorm2(1),dc_norm(1,i2+nres))
+ scal2c=scalar(diffnorm3(1),dc_norm(1,i2+nres))
+ scal3a=scalar(diffnorm1(1),dc_norm(1,i3+nres))
+ scal3b=scalar(mindiffnorm2(1),dc_norm(1,i3+nres))
+ scal3d=scalar(diffnorm2(1),dc_norm(1,i3+nres))
+ scal3c=scalar(diffnorm3(1),dc_norm(1,i3+nres))
+ scal3e=scalar(diffnorm4(1),dc_norm(1,i3+nres))
+
+
+ do l=1,3
+ pom1=diffnorm1(l)+facd1*(diffnorm1(l)-scal1a*dc_norm(l,i1+nres))
+ pom2=diffnorm2(l)+facd2*(diffnorm2(l)-scal2b*dc_norm(l,i2+nres))
+ pom3=diffnorm4(l)+facd3*(diffnorm4(l)-scal3e*dc_norm(l,i3+nres))
+
+ gradcatangc(l,i1)=gradcatangc(l,i1)&
+ +det1ij*sssmintot*(diff2(l)-diffnorm1(l)*cosom1*dista2)/(dista2*dista1)+&
+ dephiij*sssmintot*(diff3(l)-diffnorm1(l)*cosom12*dista3)/(dista3*dista1)&
+ +ene*sss2mingrad1*diffnorm1(l)*sss2min2*sss2min3
+
+
+ gradcatangc(l,i2)=gradcatangc(l,i2)+(&
+ det1ij*(diff1(l)-diffnorm2(l)*cosom1*dista1)/(dista2*dista1)+&
+ det2ij*(-diff3(l)+mindiffnorm2(l)*cosom2*dista3)/(dista3*dista2)&
+ -det2ij*(-diff2(l)-diffnorm3(l)*cosom2*dista2)/(dista3*dista2)&
+ -dephiij*(diff1(l)-diffnorm3(l)*cosom12*dista1)/(dista3*dista1))*sssmintot&
+ +ene*sss2mingrad2*diffnorm2(l)*sss2min1*sss2min3
+
+
+
+ gradcatangc(l,i3)=gradcatangc(l,i3)&
+ +det2ij*(-diff2(l)-diffnorm3(l)*cosom2*dista2)/(dista3*dista2)*sssmintot&
+ +dephiij*(diff1(l)-diffnorm3(l)*cosom12*dista1)/(dista3*dista1)*sssmintot&
+ +ene*sss2mingrad3*diffnorm4(l)*sss2min1*sss2min2
+
+
+ gradcatangc(l,j1)=gradcatangc(l,j1)-&
+ sssmintot*(det1ij*(diff2(l)-diffnorm1(l)*cosom1*dista2)/(dista2*dista1)+&
+ dephiij*(diff3(l)-diffnorm1(l)*cosom12*dista3)/(dista3*dista1))&
+ -(det1ij*(diff1(l)-diffnorm2(l)*cosom1*dista1)/(dista1*dista2)+&
+ det2ij*(-diff3(l)+mindiffnorm2(l)*cosom2*dista3)/(dista3*dista2))*sssmintot&
+ -ene*sss2mingrad1*diffnorm1(l)*sss2min2*sss2min3&
+ -ene*sss2mingrad2*diffnorm2(l)*sss2min1*sss2min3&
+ -ene*sss2mingrad3*diffnorm4(l)*sss2min1*sss2min2
+
+
+ gradcatangx(l,i1)=gradcatangx(l,i1)+(det1ij/(dista2*dista1)*&
+ (acatshiftdsc(ityptrani1,ityptranj1)*diff2(l)+&
+ facd1*(diff2(l)-scal1b*dc_norm(l,i1+nres)*dista2)-&
+ cosom1*dista2/dista1*&
+ (acatshiftdsc(ityptrani1,ityptranj1)*diff1(l)+&
+ facd1*(diff1(l)-scal1a*dc_norm(l,i1+nres)*dista1)))&
+ +dephiij/(dista3*dista1)*&
+ (acatshiftdsc(ityptrani1,ityptranj1)*diff3(l)+&
+ facd1*(diff3(l)-scal1c*dc_norm(l,i1+nres)*dista3)-&
+ cosom12*dista3/dista1*&
+ (acatshiftdsc(ityptrani1,ityptranj1)*diff1(l)+&
+ facd1*(diff1(l)-scal1a*dc_norm(l,i1+nres)*dista1))))*sssmintot&
+ +ene*sss2mingrad1*sss2min2*sss2min3*(pom1+&
+ diffnorm1(l)*(acatshiftdsc(ityptrani1,ityptranj1)-1.0d0))
+
+
+ gradcatangx(l,i3)=gradcatangx(l,i3)+(&
+ det2ij/(dista3*dista2)*&
+ (acatshiftdsc(ityptrani3,ityptranj3)*(-diff2(l))+&
+ facd3*(-diff2(l)-scal3b*dc_norm(l,i3+nres)*dista2)-&
+ cosom2*dista2/dista3*&
+ (acatshiftdsc(ityptrani3,ityptranj3)*diff3(l)+&
+ facd3*(diff3(l)-scal3c*dc_norm(l,i3+nres)*dista3)))&
+ +dephiij/(dista3*dista1)*&
+ (acatshiftdsc(ityptrani3,ityptranj3)*diff1(l)+&
+ facd3*(diff1(l)-scal3a*dc_norm(l,i3+nres)*dista1)-&
+ cosom12*dista1/dista3*&
+ (acatshiftdsc(ityptrani3,ityptranj3)*diff3(l)+&
+ facd3*(diff3(l)-scal3c*dc_norm(l,i3+nres)*dista3))))*sssmintot&
+ +ene*sss2mingrad3*sss2min2*sss2min1*(pom3+&
+ diffnorm4(l)*(acatshiftdsc(ityptrani3,ityptranj3)-1.0d0))
+
+
+ gradcatangx(l,i2)=gradcatangx(l,i2)+(&!
+ det1ij/(dista2*dista1)*&!
+ (acatshiftdsc(ityptrani2,ityptranj2)*diff1(l)&!
+ +facd2*(diff1(l)-scal2a*dc_norm(l,i2+nres)*dista1)&
+ -cosom1*dista1/dista2*&!
+ (acatshiftdsc(ityptrani2,ityptranj2)*diff2(l)+&!
+ facd2*(diff2(l)-scal2b*dc_norm(l,i2+nres)*dista2)))+&
+ det2ij/(dista3*dista2)*&!
+ (acatshiftdsc(ityptrani2,ityptranj2)*diff2(l)+&!
+ facd2*(diff2(l)-scal2b*dc_norm(l,i2+nres)*dista2)&
+ -(acatshiftdsc(ityptrani2,ityptranj2)*diff3(l)+&!
+ facd2*(diff3(l)-scal2c*dc_norm(l,i2+nres)*dista3))&
+ -cosom2*dista3/dista2*&!
+ (acatshiftdsc(ityptrani2,ityptranj2)*diff2(l)+&!
+ facd2*(diff2(l)-scal2b*dc_norm(l,i2+nres)*dista2))&
+ +cosom2*dista2/dista3*&!
+ (acatshiftdsc(ityptrani2,ityptranj2)*diff3(l)+&!
+ facd2*(diff3(l)-scal2c*dc_norm(l,i2+nres)*dista3)))&
+ +dephiij/(dista3*dista1)*&!
+ (-(acatshiftdsc(ityptrani2,ityptranj2)*diff1(l)+&!
+ facd2*(diff1(l)-scal2a*dc_norm(l,i2+nres)*dista1))+&
+ cosom12*dista1/dista3*&!
+ (acatshiftdsc(ityptrani2,ityptranj2)*diff3(l)+&!
+ facd2*(diff3(l)-scal2c*dc_norm(l,i2+nres)*dista3))))*sssmintot&
+ +ene*sss2mingrad2*sss2min3*sss2min1*(pom2+&
+ diffnorm2(l)*(acatshiftdsc(ityptrani2,ityptranj2)-1.0d0))
+
+
+ enddo
+! print *,i1,i2,i3,j1,j2,j3,"tors",ene,sinaux,cosphi
+! print *,"param",agamacattran(1,ityptrani2,ityptranj2),ityptranj2,ityptrani2
+ ecation_protang=ecation_protang+ene*sssmintot
+ enddo
+! enddo
+! enddo
+!#endif
+ return
+ end subroutine
+!--------------------------------------------------------------------------
+!c------------------------------------------------------------------------------
+ double precision function mytschebyshev(m,n,x,y,yt)
+ implicit none
+ integer i,m,n
+ double precision x(n),y,yt,yy(0:100),aux
+!c Tschebyshev polynomial in y multiplied by sin(t1)sin(t2) (yt).
+!c 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)*yt*yt
+ enddo
+ aux=0.0d0
+ do i=m,n
+ aux=aux+x(i)*yy(i)
+ enddo
+!c print *,(yy(i),i=1,n)
+ mytschebyshev=aux
+ return
+ end function
+!C--------------------------------------------------------------------------
+!C--------------------------------------------------------------------------
+ subroutine mygradtschebyshev(m,n,x,y,yt,fy,fyt)
+ implicit none
+ integer i,m,n
+ double precision x(n+1),y,yt,fy,fyt,yy(0:100),yb(0:100), &
+ ybt(0:100)
+!c Derivative of Tschebyshev polynomial in y multiplied by sin(t1)sin(t2) (yt).
+!c 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
+ yb(0)=0.0d0
+ yb(1)=1.0d0
+ ybt(0)=0.0d0
+ ybt(1)=0.0d0
+ do i=2,n
+ yy(i)=2*yy(1)*yy(i-1)-yy(i-2)*yt*yt
+ yb(i)=2*yy(i-1)+2*yy(1)*yb(i-1)-yb(i-2)*yt*yt
+ ybt(i)=2*yy(1)*ybt(i-1)-ybt(i-2)*yt*yt-2*yy(i-2)*yt
+ enddo
+ fy=0.0d0
+ fyt=0.0d0
+ do i=m,n
+ fy=fy+x(i)*yb(i)
+ fyt=fyt+x(i)*ybt(i)
+ enddo
+ return
+ end subroutine
+ subroutine fodstep(nsteps)
+ use geometry_data, only: c, nres, theta, alph
+ use geometry, only:alpha,beta,dist
+ integer, intent(in) :: nsteps
+ integer idxtomod, j, i
+ double precision RD0, RD1, fi
+! double precision alpha
+! double precision beta
+! double precision dist
+! double precision compute_RD
+ double precision TT
+ real :: r21(5)
+!c ! Założenia: dla łańcucha zapisanego w tablicy c zawierającego
+!c ! nres elementów CA i CB da się wyznaczyć kąty płaskie
+!c ! theta (procedura Alpha) i kÄ…ty torsyjne (procedura beta),
+!c ! zapisywane w tablicach theta i alph.
+!c ! Na podstawie danych z tych tablic da się odtworzyć
+!c ! strukturę 3D łańcucha procedurą chainbuild.
+!c !
+! print *,"fodstep: nres=",nres
+ RD0 = compute_RD()
+! print *, "RD0before step: ",RD0
+ do j=1,nsteps
+!c ! Wyznaczenie kątów theta na podstawie struktury
+!c ! zapisanej w tablicy c
+ do i=3,nres
+ TT=alpha(i-2,i-1,i)
+ theta(i)=TT
+!c print *,"TT=",TT
+ end do
+!c ! Wyznaczenie kątów phi na podstawie struktury
+!c ! zapisanej w tablicy c
+ do i=4,nres
+ phi(i)=beta(i-3,i-2,i-1,i)
+ end do
+!c ! Wyznaczenie odległości między atomami
+!c ! vbld(i)=dist(i-1,i)
+ do i=2,nres
+ vbld(i)=dist(i-1,i)
+ end do
+!c ! losujemy kilka liczb
+ call random_number(r21)
+!c ! r21(1): indeks pozycji do zmiany
+!c ! r21(2): kÄ…t (r21(2)/20.0-1/40.0)
+!c ! r21(3): wybór tablicy
+ RD0 = compute_RD()
+!c print *, "RD before step: ",RD0
+ fi = (r21(2)/20.0-1.0/40.0) ! o tyle radianów zmienimy losowy kąt
+ if (r21(3) .le. 0.5) then
+ idxtomod = 3+r21(1)*(nres - 2)
+ theta(idxtomod) = theta(idxtomod)+fi
+! print *,"Zmiana kÄ…ta theta(",&
+! idxtomod,") o fi = ",fi
+ else
+ idxtomod = 4+r21(1)*(nres - 3)
+ phi(idxtomod) = phi(idxtomod)+fi
+! print *,"Zmiana kÄ…ta phi(",&
+! idxtomod,") o fi = ",fi
+ end if
+!c ! odtwarzamy łańcuch
+ call chainbuild
+!c ! czy coś się polepszyło?
+ RD1 = compute_RD()
+ if (RD1 .gt. RD0) then ! nie, wycofujemy zmianÄ™
+! print *, "RD after step: ",RD1," rejected"
+ if (r21(3) .le. 0.5) then
+ theta(idxtomod) = theta(idxtomod)-fi
+ else
+ phi(idxtomod) = phi(idxtomod)-fi
+ end if
+ call chainbuild ! odtworzenie pierwotnej wersji (bez zmienionego kÄ…ta)
+ else
+! print *, "RD after step: ",RD1," accepted"
+ continue
+ end if
+ end do
+ end subroutine
+!c-----------------------------------------------------------------------------------------
+ subroutine orientation_matrix(res) ! obliczenie macierzy oraz przygotowanie ea z tymi przeksztalceniami
+ use geometry_data, only: c, nres
+ use energy_data, only: itype
+ double precision, intent(out) :: res(4,4)
+ double precision resM(4,4)
+ double precision M(4,4)
+ double precision M2(4,4)
+ integer i, j, maxi, maxj
+! double precision sq
+ double precision maxd, dd
+ double precision v1(3)
+ double precision v2(3)
+ double precision vecnea(3)
+ double precision mean_ea(3)
+ double precision fi
+!c ! liczymy atomy efektywne i zapisujemy w tablicy ea
+ do i=1,nres
+!c if (itype(i,1) .ne. 10) then
+ if (itype(i,1) .ne. 10) then
+ ea(1,i) = c(1,i+nres)
+ ea(2,i) = c(2,i+nres)
+ ea(3,i) = c(3,i+nres)
+ else
+ ea(1,i) = c(1,i)
+ ea(2,i) = c(2,i)
+ ea(3,i) = c(3,i)
+ end if
+ end do
+ call IdentityM(resM)
+ if (nres .le. 2) then
+ print *, "nres too small (should be at least 2), stopping"
+ stop
+ end if
+ do i=1,3
+ v1(i)=ea(i,1)
+ v2(i)=ea(i,2)
+ end do
+!c ! szukamy najwiekszej odleglosci miedzy atomami efektywnymi ea
+ call Dist3d(maxd,v1,v2)
+!c ! odleglosc miedzy pierwsza para atomow efektywnych
+ maxi = 1
+ maxj = 2
+ do i=1,nres-1
+ do j=i+1,nres
+ v1(1)=ea(1,i)
+ v1(2)=ea(2,i)
+ v1(3)=ea(3,i)
+ v2(1)=ea(1,j)
+ v2(2)=ea(2,j)
+ v2(3)=ea(3,j)
+ call Dist3d(dd,v1,v2)
+ if (dd .gt. maxd) then
+ maxd = dd
+ maxi = i
+ maxj = j
+ end if
+ end do
+ end do
+ vecnea(1)=ea(1,maxi)-ea(1,maxj)
+ vecnea(2)=ea(2,maxi)-ea(2,maxj)
+ vecnea(3)=ea(3,maxi)-ea(3,maxj)
+ if (vecnea(1) .lt. 0) then
+ vecnea(1) = -vecnea(1)
+ vecnea(2) = -vecnea(2)
+ vecnea(3) = -vecnea(3)
+ end if
+!c ! obliczenie kata obrotu wokol osi Z
+ fi = -atan2(vecnea(2),vecnea(1))
+ call RotateZ(M,fi)
+!c ! obliczenie kata obrotu wokol osi Y
+ fi = atan2(vecnea(3), sqrt(sq(vecnea(1))+sq(vecnea(2))))
+ call RotateY(M2,fi)
+ M = matmul(M2,M)
+!c ! Przeksztalcamy wszystkie atomy efektywne
+!c ! uzyskujac najwieksza odleglosc ulożona wzdluz OX
+!c ! ea = transform_eatoms(ea,M)
+ do i=1,nres
+ v1(1)=ea(1,i)
+ v1(2)=ea(2,i)
+ v1(3)=ea(3,i)
+ call tranform_point(v2,v1,M)
+ ea(1,i)=v2(1)
+ ea(2,i)=v2(2)
+ ea(3,i)=v2(3)
+ end do
+ resM = M
+!c ! Teraz szukamy najdluzszego rzutu na plaszczyzne YZ
+!c ! (czyli w liczeniu odleglosci bierzemy pod uwage tylko wsp. y, z)
+ maxd = sqrt( sq(ea(2,1)-ea(2,2)) + sq(ea(3,1)-ea(3,2))) ! aktualnie max odl
+ maxi = 1 ! indeksy atomow
+ maxj = 2 ! miedzy ktorymi jest max odl (chwilowe)
+ do i=1,nres-1
+ do j=i+1,nres
+ dd = sqrt( (ea(2,i)-ea(2,j))**2 + (ea(3,i)-ea(3,j))**2)
+ if (dd .gt. maxd) then
+ maxd = dd
+ maxi = i
+ maxj = j
+ end if
+ end do
+ end do
+!c ! Teraz obrocimy wszystko wokol OX tak, zeby znaleziony rzut
+!c ! byl rownolegly do OY
+ vecnea(1) = ea(1,maxi)-ea(1,maxj)
+ vecnea(2) = ea(2,maxi)-ea(2,maxj)
+ vecnea(3) = ea(3,maxi)-ea(3,maxj)
+!c ! jeśli współrzędna vecnea.y < 0, to robimy odwrotnie
+ if (vecnea(2) .lt. 0) then
+ vecnea(1) = -vecnea(1)
+ vecnea(2) = -vecnea(2)
+ vecnea(3) = -vecnea(3)
+ end if
+!c ! obliczenie kąta obrotu wokół osi X
+ fi = -atan2(vecnea(3),vecnea(2))
+ call RotateX(M,fi)
+!c ! Przeksztalcamy wszystkie atomy efektywne
+ do i=1,nres
+ v1(1)=ea(1,i)
+ v1(2)=ea(2,i)
+ v1(3)=ea(3,i)
+ call tranform_point(v2,v1,M)
+ ea(1,i)=v2(1)
+ ea(2,i)=v2(2)
+ ea(3,i)=v2(3)
+ end do
+ resM = matmul(M,resM) ! zbieramy wynik (sprawdzic kolejnosc M,resM)
+!c ! centrujemy
+ mean_ea(1) = 0
+ mean_ea(2) = 0
+ mean_ea(3) = 0
+ do i=1,nres
+ mean_ea(1) = mean_ea(1) + ea(1,i)
+ mean_ea(2) = mean_ea(2) + ea(2,i)
+ mean_ea(3) = mean_ea(3) + ea(3,i)
+ end do
+ v1(1) = -mean_ea(1)/nres
+ v1(2) = -mean_ea(2)/nres
+ v1(3) = -mean_ea(3)/nres
+ call TranslateV(M,v1)
+ resM = matmul(M,resM)
+!c ! przesuwamy
+ do i=1,nres
+ ea(1,i) = ea(1,i) + v1(1)
+ ea(2,i) = ea(2,i) + v1(2)
+ ea(3,i) = ea(3,i) + v1(3)
+ end do
+ res = resM
+!c ! wynikowa macierz przeksztalcenia lancucha
+!c ! (ale lancuch w ea juz mamy przeksztalcony)
+ return
+ end subroutine
+ double precision function compute_rd
+ use geometry_data, only: nres
+ use energy_data, only: itype
+ implicit none
+ double precision or_mat(4,4)
+! double precision hydrophobicity
+ integer neatoms
+ double precision cutoff
+ double precision ho(70000)
+ double precision ht(70000)
+ double precision hosum, htsum
+ double precision marg, sigmax, sigmay, sigmaz
+ integer i, j
+ double precision v1(3)
+ double precision v2(3)
+ double precision rijdivc, coll, tmpkwadrat, tmppotega, dist
+ double precision OdivT, OdivR, ot_one, or_one, RD_classic
+ call orientation_matrix(or_mat)
+!c ! tam juz liczy sie tablica ea
+ neatoms = nres
+ cutoff = 8.99d0
+!c ! granica oddzialywania w A (powyzej ignorujemy oddzialywanie)
+!c ! Najpierw liczymy "obserwowana hydrofobowosc"
+ hosum = 0.0d0 ! na sume pol ho, do celow pozniejszej normalizacji
+ do j=1,neatoms
+ ho(j)=0.0d0
+ do i=1,neatoms
+ if (j .eq. i) then ! nie uwzgledniamy oddzialywania atomu z samym soba
+ cycle
+ end if
+ v1(1)=ea(1,i)
+ v1(2)=ea(2,i)
+ v1(3)=ea(3,i)
+ v2(1)=ea(1,j)
+ v2(2)=ea(2,j)
+ v2(3)=ea(3,j)
+ call Dist3d(dist,v1,v2) ! odleglosc miedzy atomami
+ if (dist .gt. cutoff) then ! za daleko, nie uwzgledniamy
+ cycle
+ end if
+ rijdivc = dist / cutoff
+ coll = 0.0d0
+ tmppotega = rijdivc*rijdivc
+ tmpkwadrat = tmppotega
+ coll = coll + 7*tmpkwadrat
+ tmppotega = tmppotega * tmpkwadrat ! do potęgi 4
+ coll = coll - 9*tmppotega
+ tmppotega = tmppotega * tmpkwadrat ! do potęgi 6
+ coll = coll + 5*tmppotega
+ tmppotega = tmppotega * tmpkwadrat ! do potęgi 8
+ coll = coll - tmppotega
+!c ! Wersja: Bryliński 2007
+!c ! EAtoms[j].collectedhp += EAtoms[i].hyphob*(1 - 0.5 * coll);
+!c ! ea$ho[j] = ea$ho[j] + hydrophobicity(ea$resid[i])*(1-0.5*coll)
+!c ! Wersja: Banach Konieczny Roterman 2014
+!c ! EAtoms[j].collectedhp += (EAtoms[i].hyphob+EAtoms[j].hyphob)*(1 - 0.5 * coll);
+!c ponizej bylo itype(i,1) w miejscu itype(i) oraz itype(j,1) w miejscu itype(j)
+ ho(j) = ho(j) + (hydrophobicity(itype(i,1))+&
+ hydrophobicity(itype(j,1)))*(1.0d0-0.5_8*coll)
+ end do
+ hosum = hosum + ho(j)
+ end do
+!c ! Normalizujemy
+ do i=1,neatoms
+ ho(i) = ho(i) / hosum
+ end do
+!c ! Koniec liczenia hydrofobowosci obserwowanej (profil ho)
+!c ! Teraz liczymy "teoretyczna hydrofobowosc", wedlug kropli i rozkladu Gaussa
+ htsum = 0.0d0
+!c ! tu zbieramy sume ht, uzyjemy potem do normalizacji
+!c ! Ustalimy teraz parametry rozkladu Gaussa, czyli sigmy (srodek jest w (0,0,0)).
+!c ! To bedzie (max odl od srodka + margines) / 3, oddzielnie dla kazdej wspolrzednej.
+ marg = 9.0d0
+ htsum = 0.0d0
+!c ! jeszcze raz zerujemy
+!c ! szukamy ekstremalnej wartosci wspolrzednej x (max wart bezwzgl)
+ sigmax = ea(1,1)
+ do i=2,neatoms
+ if (abs(ea(1,i))>sigmax) then
+ sigmax = abs(ea(1,i))
+ end if
+ end do
+ sigmax = (marg + sigmax) / 3.0d0
+!c ! szukamy ekstremalnej wartosci wspolrzednej y (max wart bezwzgl)
+ sigmay = ea(2,1)
+ do i=2,neatoms
+ if (abs(ea(2,i))>sigmay) then
+ sigmay = abs(ea(2,i))
+ end if
+ end do
+ sigmay = (marg + sigmay) / 3.0d0
+!c ! szukamy ekstremalnej wartosci wspolrzednej z (max wart bezwzgl)
+ sigmaz = ea(3,1)
+ do i=2,neatoms
+ if (abs(ea(3,i))>sigmaz) then
+ sigmaz = abs(ea(3,i))
+ end if
+ end do
+ sigmaz = (marg + sigmaz) / 3.0d0
+!c !sigmax = (marg + max(abs(max(ea$acoor[,1])), abs(min(ea$acoor[,1]))))/3.0
+!c !sigmay = (marg + max(abs(max(ea$acoor[,2])), abs(min(ea$acoor[,2]))))/3.0
+!c !sigmaz = (marg + max(abs(max(ea$acoor[,3])), abs(min(ea$acoor[,3]))))/3.0
+!c ! print *,"sigmax =",sigmax," sigmay =",sigmay," sigmaz = ",sigmaz
+ do j=1,neatoms
+ ht(j)= exp(-(ea(1,j))**2/(2*sigmax**2))&
+ * exp(-(ea(2,j))**2/(2*sigmay**2)) &
+ * exp(-(ea(3,j))**2/(2*sigmaz**2))
+ htsum = htsum + ht(j)
+ end do
+!c ! Normalizujemy
+ do i=1, neatoms
+ ht(i) = ht(i) / htsum
+ end do
+!c ! Teraz liczymy RD
+ OdivT = 0.0d0
+ OdivR = 0.0d0
+ do j=1,neatoms
+ if (ho(j) .ne. 0) then
+ ot_one = ho(j) * log(ho(j)/ht(j)) / log(2.0d0)
+ OdivT = OdivT + ot_one
+ or_one = ho(j) * log(ho(j)/ (1.0d0/neatoms)) / log(2.0_8)
+ OdivR = OdivR + or_one
+ endif
+ end do
+ RD_classic = OdivT / (OdivT+OdivR)
+ compute_rd = RD_classic
+ return
+ end function
+ function hydrophobicity(id) ! do przepisania (bylo: identyfikowanie aa po nazwach)
+ integer id
+ double precision hydrophobicity
+ hydrophobicity = 0.0d0
+ if (id .eq. 1) then
+ hydrophobicity = 1.000d0 ! CYS
+ return
+ endif
+ if (id .eq. 2) then
+ hydrophobicity = 0.828d0 ! MET
+ return
+ endif
+ if (id .eq. 3) then
+ hydrophobicity = 0.906d0 ! PHE
+ return
+ endif
+ if (id .eq. 4) then
+ hydrophobicity = 0.883d0 ! ILE
+ return
+ endif
+ if (id .eq. 5) then
+ hydrophobicity = 0.783d0 ! LEU
+ return
+ endif
+ if (id .eq. 6) then
+ hydrophobicity = 0.811d0 ! VAL
+ return
+ endif
+ if (id .eq. 7) then
+ hydrophobicity = 0.856d0 ! TRP
+ return
+ endif
+ if (id .eq. 8) then
+ hydrophobicity = 0.700d0 ! TYR
+ return
+ endif
+ if (id .eq. 9) then
+ hydrophobicity = 0.572d0 ! ALA
+ return
+ endif
+ if (id .eq. 10) then
+ hydrophobicity = 0.550d0 ! GLY
+ return
+ endif
+ if (id .eq. 11) then
+ hydrophobicity = 0.478d0 ! THR
+ return
+ endif
+ if (id .eq. 12) then
+ hydrophobicity = 0.422d0 ! SER
+ return
+ endif
+ if (id .eq. 13) then
+ hydrophobicity = 0.250d0 ! GLN
+ return
+ endif
+ if (id .eq. 14) then
+ hydrophobicity = 0.278d0 ! ASN
+ return
+ endif
+ if (id .eq. 15) then
+ hydrophobicity = 0.083d0 ! GLU
+ return
+ endif
+ if (id .eq. 16) then
+ hydrophobicity = 0.167d0 ! ASP
+ return
+ endif
+ if (id .eq. 17) then
+ hydrophobicity = 0.628d0 ! HIS
+ return
+ endif
+ if (id .eq. 18) then
+ hydrophobicity = 0.272d0 ! ARG
+ return
+ endif
+ if (id .eq. 19) then
+ hydrophobicity = 0.000d0 ! LYS
+ return
+ endif
+ if (id .eq. 20) then
+ hydrophobicity = 0.300d0 ! PRO
+ return
+ endif
+ return
+ end function hydrophobicity
+ subroutine mycrossprod(res,b,c)
+ implicit none
+ double precision, intent(out) :: res(3)
+ double precision, intent(in) :: b(3)
+ double precision, intent(in) :: c(3)
+!c ! Tylko dla wektorów trójwymiarowych, ale nie sprawdzamy tego tutaj
+ res(1) = b(2)*c(3)-b(3)*c(2)
+ res(2) = b(3)*c(1)-b(1)*c(3)
+ res(3) = b(1)*c(2)-b(2)*c(1)
+ return
+ end subroutine
+ subroutine mydotprod(res,b,c)
+ implicit none
+ double precision, intent(out) :: res
+ double precision, intent(in) :: b(3)
+ double precision, intent(in) :: c(3)
+!c ! Tylko dla wektorów trójwymiarowych, ale nie sprawdzamy tego tutaj
+ res = b(1)*c(1)+b(2)*c(2)+b(3)*c(3)
+ return
+ end subroutine
+!c ! cosinus k¹ta miêdzy wektorami trójwymiarowymi
+ subroutine cosfi(res, x, y)
+ implicit none
+ double precision, intent(out) :: res
+ double precision, intent(in) :: x(3)
+ double precision, intent(in) :: y(3)
+ double precision LxLy
+ LxLy=sqrt(x(1)*x(1)+x(2)*x(2)+x(3)*x(3)) *&
+ sqrt(y(1)*y(1)+y(2)*y(2)+y(3)*y(3))
+ if (LxLy==0.0) then
+ res = 0.0d0
+ else
+ call mydotprod(res,x,y)
+ res = res / LxLy
+ end if
+ return
+ end subroutine
+
+
+ subroutine Dist3d(res,v1,v2)
+ implicit none
+ double precision, intent(out) :: res
+ double precision, intent(in) :: v1(3)
+ double precision, intent(in) :: v2(3)
+! double precision sq
+ res = sqrt( sq(v1(1)-v2(1)) + sq(v1(2)-v2(2)) + sq(v1(3)-v2(3)))
+ return
+ end subroutine
+!c ! Przeksztalca wsp. 3d uzywajac macierzy przeksztalcenia M (4x4)
+ subroutine tranform_point(res,v3d,M)
+ implicit none
+ double precision, intent(out) :: res(3)
+ double precision, intent(in) :: v3d(3)
+ double precision, intent(in) :: M(4,4)
+
+ res(1) = M(1,1)*v3d(1) + M(1,2)*v3d(2) + M(1,3)*v3d(3) + M(1,4)
+ res(2) = M(2,1)*v3d(1) + M(2,2)*v3d(2) + M(2,3)*v3d(3) + M(2,4)
+ res(3) = M(3,1)*v3d(1) + M(3,2)*v3d(2) + M(3,3)*v3d(3) + M(3,4)
+ return
+ end subroutine
+!c ! TranslateV: macierz translacji o wektor V
+ subroutine TranslateV(res,V)
+ implicit none
+ double precision, intent(out) :: res(4,4)
+ double precision, intent(in) :: v(3)
+ res(1,1) = 1.0d0
+ res(1,2) = 0
+ res(1,3) = 0
+ res(1,4) = v(1)
+ res(2,1) = 0
+ res(2,2) = 1.0d0
+ res(2,3) = 0
+ res(2,4) = v(2)
+ res(3,1) = 0
+ res(3,2) = 0
+ res(3,3) = 1.0d0
+ res(3,4) = v(3)
+ res(4,1) = 0
+ res(4,2) = 0
+ res(4,3) = 0
+ res(4,4) = 1.0d0
+ return
+ end subroutine
+!c ! RotateX: macierz obrotu wokol osi OX o kat fi
+ subroutine RotateX(res,fi)
+ implicit none
+ double precision, intent(out) :: res(4,4)
+ double precision, intent(in) :: fi
+ res(1,1) = 1.0d0
+ res(1,2) = 0
+ res(1,3) = 0
+ res(1,4) = 0
+ res(2,1) = 0
+ res(2,2) = cos(fi)
+ res(2,3) = -sin(fi)
+ res(2,4) = 0
+ res(3,1) = 0
+ res(3,2) = sin(fi)
+ res(3,3) = cos(fi)
+ res(3,4) = 0
+ res(4,1) = 0
+ res(4,2) = 0
+ res(4,3) = 0
+ res(4,4) = 1.0d0
+ return
+ end subroutine
+!c ! RotateY: macierz obrotu wokol osi OY o kat fi
+ subroutine RotateY(res,fi)
+ implicit none
+ double precision, intent(out) :: res(4,4)
+ double precision, intent(in) :: fi
+ res(1,1) = cos(fi)
+ res(1,2) = 0
+ res(1,3) = sin(fi)
+ res(1,4) = 0
+ res(2,1) = 0
+ res(2,2) = 1.0d0
+ res(2,3) = 0
+ res(2,4) = 0
+ res(3,1) = -sin(fi)
+ res(3,2) = 0
+ res(3,3) = cos(fi)
+ res(3,4) = 0
+ res(4,1) = 0
+ res(4,2) = 0
+ res(4,3) = 0
+ res(4,4) = 1.0d0
+ return
+ end subroutine
+!c ! RotateZ: macierz obrotu wokol osi OZ o kat fi
+ subroutine RotateZ(res,fi)
+ implicit none
+ double precision, intent(out) :: res(4,4)
+ double precision, intent(in) :: fi
+ res(1,1) = cos(fi)
+ res(1,2) = -sin(fi)
+ res(1,3) = 0
+ res(1,4) = 0
+ res(2,1) = sin(fi)
+ res(2,2) = cos(fi)
+ res(2,3) = 0
+ res(2,4) = 0
+ res(3,1) = 0
+ res(3,2) = 0
+ res(3,3) = 1.0d0
+ res(3,4) = 0
+ res(4,1) = 0
+ res(4,2) = 0
+ res(4,3) = 0
+ res(4,4) = 1.0d0
+ return
+ end subroutine
+!c ! IdentityM
+ subroutine IdentityM(res)
+ implicit none
+ double precision, intent(out) :: res(4,4)
+ res(1,1) = 1.0d0
+ res(1,2) = 0
+ res(1,3) = 0
+ res(1,4) = 0
+ res(2,1) = 0
+ res(2,2) = 1.0d0
+ res(2,3) = 0
+ res(2,4) = 0
+ res(3,1) = 0
+ res(3,2) = 0
+ res(3,3) = 1.0d0
+ res(3,4) = 0
+ res(4,1) = 0
+ res(4,2) = 0
+ res(4,3) = 0
+ res(4,4) = 1.0d0
+ return
+ end subroutine
+ double precision function sq(x)
+ double precision x
+ sq = x*x
+ return
+ end function sq
+
+#ifdef LBFGS
+ double precision function funcgrad(x,g)
+ use MD_data, only: totT,usampl
+ implicit none
+ double precision energia(0:n_ene)
+ double precision x(nvar),g(nvar)
+ integer i
+ call var_to_geom(nvar,x)
+ call zerograd
+ call chainbuild
+ call etotal(energia(0))
+ call sum_gradient
+ funcgrad=energia(0)
+ call cart2intgrad(nvar,g)
+ if (usampl) then
+ do i=1,nres-3
+ gloc(i,icg)=gloc(i,icg)+dugamma(i)
+ enddo
+ do i=1,nres-2
+ gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
+ enddo
+ endif
+ do i=1,nvar
+ g(i)=g(i)+gloc(i,icg)
+ enddo
+ return
+ end function funcgrad
+ subroutine cart2intgrad(n,g)
+ integer n
+ double precision g(n)
+ double precision drt(3,3,nres),rdt(3,3,nres),dp(3,3),&
+ temp(3,3),prordt(3,3,nres),prodrt(3,3,nres)
+ double precision xx(3),xx1(3),alphi,omegi,xj,dpjk,yp,xp,xxp,yyp
+ double precision cosalphi,sinalphi,cosomegi,sinomegi,theta2,&
+ cost2,sint2,rj,dxoiij,tempkl,dxoijk,dsci,zzp,dj,dpkl
+ double precision fromto(3,3),aux(6)
+ integer i,ii,j,jjj,k,l,m,indi,ind,ind1
+ logical sideonly
+ sideonly=.false.
+ g=0.0d0
+ if (sideonly) goto 10
+ do i=1,nres-2
+ rdt(1,1,i)=-rt(1,2,i)
+ rdt(1,2,i)= rt(1,1,i)
+ rdt(1,3,i)= 0.0d0
+ rdt(2,1,i)=-rt(2,2,i)
+ rdt(2,2,i)= rt(2,1,i)
+ rdt(2,3,i)= 0.0d0
+ rdt(3,1,i)=-rt(3,2,i)
+ rdt(3,2,i)= rt(3,1,i)
+ rdt(3,3,i)= 0.0d0
+ enddo
+ do i=2,nres-2
+ drt(1,1,i)= 0.0d0
+ drt(1,2,i)= 0.0d0
+ drt(1,3,i)= 0.0d0
+ drt(2,1,i)= rt(3,1,i)
+ drt(2,2,i)= rt(3,2,i)
+ drt(2,3,i)= rt(3,3,i)
+ drt(3,1,i)=-rt(2,1,i)
+ drt(3,2,i)=-rt(2,2,i)
+ drt(3,3,i)=-rt(2,3,i)
+ enddo
+ ind1=0
+ do i=1,nres-2
+ ind1=ind1+1
+ if (n.gt.nphi) then
+
+ do j=1,3
+ do k=1,2
+ dpjk=0.0D0
+ do l=1,3
+ dpjk=dpjk+prod(j,l,i)*rdt(l,k,i)
+ enddo
+ dp(j,k)=dpjk
+ prordt(j,k,i)=dp(j,k)
+ enddo
+ dp(j,3)=0.0D0
+ g(nphi+i)=g(nphi+i)+vbld(i+2)*dp(j,1)*gradc(j,i+1,icg)
+ enddo
+ xx1(1)=-0.5D0*xloc(2,i+1)
+ xx1(2)= 0.5D0*xloc(1,i+1)
+ do j=1,3
+ xj=0.0D0
+ do k=1,2
+ xj=xj+r(j,k,i)*xx1(k)
+ enddo
+ xx(j)=xj
+ enddo
+ do j=1,3
+ rj=0.0D0
+ do k=1,3
+ rj=rj+prod(j,k,i)*xx(k)
+ enddo
+ g(nphi+i)=g(nphi+i)+rj*gradx(j,i+1,icg)
+ enddo
+ if (i.lt.nres-2) then
+ do j=1,3
+ dxoiij=0.0D0
+ do k=1,3
+ dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
+ enddo
+ g(nphi+i)=g(nphi+i)+dxoiij*gradx(j,i+2,icg)
+ enddo
+ endif
+
+ endif
+
+
+ if (i.gt.1) then
+ do j=1,3
+ do k=1,3
+ dpjk=0.0
+ do l=2,3
+ dpjk=dpjk+prod(j,l,i)*drt(l,k,i)
+ enddo
+ dp(j,k)=dpjk
+ prodrt(j,k,i)=dp(j,k)
+ enddo
+ g(i-1)=g(i-1)+vbld(i+2)*dp(j,1)*gradc(j,i+1,icg)
+ enddo
+ endif
+ xx(1)= 0.0D0
+ xx(3)= xloc(2,i+1)*r(2,2,i)+xloc(3,i+1)*r(2,3,i)
+ xx(2)=-xloc(2,i+1)*r(3,2,i)-xloc(3,i+1)*r(3,3,i)
+ if (i.gt.1) then
+ do j=1,3
+ rj=0.0D0
+ do k=2,3
+ rj=rj+prod(j,k,i)*xx(k)
+ enddo
+ g(i-1)=g(i-1)-rj*gradx(j,i+1,icg)
+ enddo
+ endif
+ if (i.gt.1) then
+ do j=1,3
+ dxoiij=0.0D0
+ do k=1,3
+ dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
+ enddo
+ g(i-1)=g(i-1)+dxoiij*gradx(j,i+2,icg)
+ enddo
+ endif
+ do j=i+1,nres-2
+ ind1=ind1+1
+ call build_fromto(i+1,j+1,fromto)
+ do k=1,3
+ do l=1,3
+ tempkl=0.0D0
+ do m=1,2
+ tempkl=tempkl+prordt(k,m,i)*fromto(m,l)
+ enddo
+ temp(k,l)=tempkl
+ enddo
+ enddo
+ if (n.gt.nphi) then
+ do k=1,3
+ g(nphi+i)=g(nphi+i)+vbld(j+2)*temp(k,1)*gradc(k,j+1,icg)
+ enddo
+ do k=1,3
+ dxoijk=0.0D0
+ do l=1,3
+ dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
+ enddo
+ g(nphi+i)=g(nphi+i)+dxoijk*gradx(k,j+2,icg)
+ enddo
+ endif
+ do k=1,3
+ do l=1,3
+ tempkl=0.0D0
+ do m=1,3
+ tempkl=tempkl+prodrt(k,m,i)*fromto(m,l)
+ enddo
+ temp(k,l)=tempkl
+ enddo
+ enddo
+ if (i.gt.1) then
+ do k=1,3
+ g(i-1)=g(i-1)+vbld(j+2)*temp(k,1)*gradc(k,j+1,icg)
+ enddo
+ do k=1,3
+ dxoijk=0.0D0
+ do l=1,3
+ dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
+ enddo
+ g(i-1)=g(i-1)+dxoijk*gradx(k,j+2,icg)
+ enddo
+ endif
+ enddo
+ enddo
+
+ if (nvar.le.nphi+ntheta) return
+
+ 10 continue
+ do i=2,nres-1
+ if (iabs(itype(i,1)).eq.10 .or. itype(i,1).eq.ntyp1& !) cycle
+ .or. mask_side(i).eq.0 ) cycle
+ ii=ialph(i,1)
+ dsci=vbld(i+nres)
+#ifdef OSF
+ alphi=alph(i)
+ omegi=omeg(i)
+ if(alphi.ne.alphi) alphi=100.0
+ if(omegi.ne.omegi) omegi=-100.0
+#else
+ alphi=alph(i)
+ omegi=omeg(i)
+#endif
+ cosalphi=dcos(alphi)
+ sinalphi=dsin(alphi)
+ cosomegi=dcos(omegi)
+ sinomegi=dsin(omegi)
+ temp(1,1)=-dsci*sinalphi
+ temp(2,1)= dsci*cosalphi*cosomegi
+ temp(3,1)=-dsci*cosalphi*sinomegi
+ temp(1,2)=0.0D0
+ temp(2,2)=-dsci*sinalphi*sinomegi
+ temp(3,2)=-dsci*sinalphi*cosomegi
+ theta2=pi-0.5D0*theta(i+1)
+ cost2=dcos(theta2)
+ sint2=dsin(theta2)
+ jjj=0
+ do j=1,2
+ xp=temp(1,j)
+ yp=temp(2,j)
+ xxp= xp*cost2+yp*sint2
+ yyp=-xp*sint2+yp*cost2
+ zzp=temp(3,j)
+ xx(1)=xxp
+ xx(2)=yyp*r(2,2,i-1)+zzp*r(2,3,i-1)
+ xx(3)=yyp*r(3,2,i-1)+zzp*r(3,3,i-1)
+ do k=1,3
+ dj=0.0D0
+ do l=1,3
+ dj=dj+prod(k,l,i-1)*xx(l)
+ enddo
+ aux(jjj+k)=dj
+ enddo
+ jjj=jjj+3
+ enddo
+ do k=1,3
+ g(ii)=g(ii)+aux(k)*gradx(k,i,icg)
+ g(ii+nside)=g(ii+nside)+aux(k+3)*gradx(k,i,icg)
+ enddo
+ enddo
+ return
+ end subroutine cart2intgrad
+
+
+#endif
+
+!-----------LIPID-MARTINI-UNRES-PROTEIN
+
+! new for K+
+ subroutine elip_prot(evdw)
+! subroutine emart_prot2(emartion_prot)
+ use calc_data
+ use comm_momo
+
+ logical :: lprn
+!el local variables
+ integer :: iint,itypi1,subchap,isel,itmp
+ 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,ki
+ 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,&
+ emartions_prot_amber,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)
+ real(kind=8) ::locbox(3)
+ locbox(1)=boxxsize
+ locbox(2)=boxysize
+ locbox(3)=boxzsize
+
+ evdw=0.0D0
+ if (nres_molec(4).eq.0) return
+ eps_out=80.0d0
+! sss_ele_cut=1.0d0
+
+ itmp=0
+ do i=1,4
+ itmp=itmp+nres_molec(i)
+ enddo
+! go to 17
+! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization
+! do i=ibond_start,ibond_end
+ do ki=g_listmartsc_start,g_listmartsc_end
+ i=newcontlistmartsci(ki)
+ j=newcontlistmartscj(ki)
+
+! print *,"I am in EVDW",i
+ itypi=iabs(itype(i,1))
+
+! if (i.ne.47) cycle
+ if ((itypi.eq.ntyp1).or.(itypi.eq.10)) 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)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+ dsci_inv=vbld_inv(i+nres)
+! do j=itmp+1,itmp+nres_molec(5)
+
+! Calculate SC interaction energy.
+ itypj=iabs(itype(j,4))
+ if ((itypj.gt.ntyp_molec(4))) cycle
+ CALL elgrad_init_mart(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+! print *,i,j,"after elgrad"
+ dscj_inv=0.0
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+
+ call to_box(xj,yj,zj)
+! 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,*) "xj,yj,zj", xj,yj,zj,boxxsize
+ rreal(1)=xj
+ rreal(2)=yj
+ rreal(3)=zj
+ dxj=0.0
+ dyj=0.0
+ dzj=0.0
+! dxj = dc_norm( 1, nres+j )
+! dyj = dc_norm( 2, nres+j )
+! dzj = dc_norm( 3, nres+j )
+
+ itypi = itype(i,1)
+ itypj = itype(j,4)
+! Parameters from fitting the analitical expressions to the PMF obtained by umbrella
+! sampling performed with amber package
+! alf1 = 0.0d0
+! alf2 = 0.0d0
+! alf12 = 0.0d0
+! a12sq = rborn(itypi,itypj) * rborn(itypj,itypi)
+ chi1 = chi1mart(itypi,itypj)
+ chis1 = chis1mart(itypi,itypj)
+ chip1 = chipp1mart(itypi,itypj)
+! chi1=0.0d0
+! chis1=0.0d0
+! chip1=0.0d0
+ chi2=0.0
+ chip2=0.0
+ chis2=0.0
+! chis2 = chis(itypj,itypi)
+ chis12 = chis1 * chis2
+ sig1 = sigmap1mart(itypi,itypj)
+ sig2=0.0d0
+! sig2 = sigmap2(itypi,itypj)
+! alpha factors from Fcav/Gcav
+ b1cav = alphasurmart(1,itypi,itypj)
+ b2cav = alphasurmart(2,itypi,itypj)
+ b3cav = alphasurmart(3,itypi,itypj)
+ b4cav = alphasurmart(4,itypi,itypj)
+
+! b1cav=0.0d0
+! b2cav=0.0d0
+! b3cav=0.0d0
+! b4cav=0.0d0
+
+! used to determine whether we want to do quadrupole calculations
+ eps_in = epsintabmart(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
+
+ DO k = 1, 3
+ ctail(k,1)=c(k,i+nres)
+ ctail(k,2)=c(k,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)
+ do k=1,3
+ Rtail_distance(k) = boxshift(ctail(k,2) - ctail(k,1),locbox(k))
+ enddo
+ Rtail = dsqrt( &
+ (Rtail_distance(1)*Rtail_distance(1)) &
+ + (Rtail_distance(2)*Rtail_distance(2)) &
+ + (Rtail_distance(3)*Rtail_distance(3)))
+! tail lomartion and distance calculations
+! dhead1
+ d1 = dheadmart(1, 1, itypi, itypj)
+! d2 = dhead(2, 1, itypi, itypj)
+ DO k = 1,3
+! lomartion of polar head is computed by taking hydrophobic centre
+! and moving by a d1 * dc_norm vector
+! see unres publimartions for very informative images
+ chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres)
+ chead(k,2) = c(k, j)
+ 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))
+! write(iout,*) "TEST",chead(1,1),chead(2,1),chead(3,1),dc_norm(k, i+nres),d1
+! 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
+ Fisocav=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))
+! 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
+ if (evdw.gt.1.0d6) then
+ write (*,'(2(1x,a3,i3),7f7.2)') &
+ restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
+ 1.0d0/rij,Rtail,Rhead,rij_shift, sig, sig0ij,sigsq
+ write(*,*) facsig,faceps1_inv,om1,chiom1,chi1
+ write(*,*) "ANISO?!",chi1
+!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_mart(itypi,itypj)
+! print *,"ADAM",aa_aq(itypi,itypj)
+
+! c1 = 0.0d0
+ c2 = fac * bb_aq_mart(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
+! print *,"GG(1),distance grad",gg(1)
+ fac = chis1 * sqom1 + chis2 * sqom2 &
+ - 2.0d0 * chis12 * om1 * om2 * om12
+ pom = 1.0d0 - chis1 * chis2 * sqom12
+ 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
+ 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 )
+
+ 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) )
+ facd1 = dtailmart(1,itypi,itypj) * vbld_inv(i+nres)
+ facd2 = dtailmart(2,itypi,itypj) * vbld_inv(j)
+ DO k = 1, 3
+ pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres))
+ gradpepmartx(k,i) = gradpepmartx(k,i) &
+ - (( dFdR + gg(k) ) * pom)*sss_ele_cut&
+ -(evdwij+Fcav)*rij*sss_ele_grad*rreal(k)
+
+ pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j))
+! gvdwx(k,j) = gvdwx(k,j) &
+! + (( dFdR + gg(k) ) * pom)
+ gradpepmart(k,i) = gradpepmart(k,i) &
+ - (( dFdR + gg(k) ) * ertail(k))*sss_ele_cut&
+ -(evdwij+Fcav)*rij*sss_ele_grad*rreal(k)
+
+ gradpepmart(k,j) = gradpepmart(k,j) &
+ + (( dFdR + gg(k) ) * ertail(k))*sss_ele_cut&
+ +(evdwij+Fcav)*rij*sss_ele_grad*rreal(k)
+
+ gg(k) = 0.0d0
+ ENDDO
+!c! Compute head-head and head-tail energies for each state
+!! if (.false.) then ! turn off electrostatic
+ isel = iabs(Qi)+iabs(Qj)
+ if ((itype(j,4).gt.4).and.(itype(j,4).lt.14)) isel=isel+2
+! isel=0
+! if (isel.eq.2) isel=0
+ IF (isel.le.1) THEN
+ eheadtail = 0.0d0
+ ELSE IF (isel.eq.3) THEN
+ if (iabs(Qj).eq.1) then
+ CALL edq_mart(ecl, elj, epol)
+ eheadtail = ECL + elj + epol
+ else
+ 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
+ call eqd_mart(ecl,elj,epol)
+ eheadtail = ECL + elj + epol
+ endif
+ ELSE IF ((isel.eq.2)) THEN
+ if (iabs(Qi).ne.1) then
+ eheadtail=0.0d0
+ else
+ 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
+ CALL eqq_mart(Ecl,Egb,Epol,Fisocav,Elj)
+ eheadtail = ECL + Egb + Epol + Fisocav + Elj
+ endif
+ ELSE IF (isel.eq.4) then
+ call edd_mart(ecl)
+ eheadtail = ECL
+ ENDIF
+! write(iout,*) "not yet implemented",j,itype(j,5)
+!! endif ! turn off electrostatic
+ evdw = evdw + (Fcav + eheadtail)*sss_ele_cut
+! 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 (energy_dec) write(iout,*) "FCAV", &
+ sig1,sig2,b1cav,b2cav,b3cav,b4cav
+! print *,"before sc_grad_mart", i,j, gradpepmart(1,j)
+! iF (nstate(itypi,itypj).eq.1) THEN
+ CALL sc_grad_mart
+! print *,"after sc_grad_mart", i,j, gradpepmart(1,j)
+
+! 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
+!!! return
+ 17 continue
+! go to 23
+! do i=ibond_start,ibond_end
+
+ do ki=g_listmartp_start,g_listmartp_end
+ i=newcontlistmartpi(ki)
+ j=newcontlistmartpj(ki)
+
+! print *,"I am in EVDW",i
+ itypi=10 ! the peptide group parameters are for glicine
+
+! if (i.ne.47) cycle
+ if ((itype(i,1).eq.ntyp1).or.itype(i+1,1).eq.ntyp1) cycle
+ itypi1=iabs(itype(i+1,1))
+ 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)
+ dxi=dc_norm(1,i)
+ dyi=dc_norm(2,i)
+ dzi=dc_norm(3,i)
+ dsci_inv=vbld_inv(i+1)/2.0
+! do j=itmp+1,itmp+nres_molec(5)
+
+! Calculate SC interaction energy.
+ itypj=iabs(itype(j,4))
+ if ((itypj.gt.ntyp_molec(4))) cycle
+ CALL elgrad_init_mart_pep(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+
+ dscj_inv=0.0
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ rreal(1)=xj
+ rreal(2)=yj
+ rreal(3)=zj
+
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+
+ dxj = 0.0d0! dc_norm( 1, nres+j )
+ dyj = 0.0d0!dc_norm( 2, nres+j )
+ dzj = 0.0d0! dc_norm( 3, nres+j )
+
+ itypi = 10
+ itypj = itype(j,4)
+! Parameters from fitting the analitical expressions to the PMF obtained by umbrella
+! sampling performed with amber package
+! alf1 = 0.0d0
+! alf2 = 0.0d0
+! alf12 = 0.0d0
+! a12sq = rborn(itypi,itypj) * rborn(itypj,itypi)
+ chi1 = chi1mart(itypi,itypj)
+ chis1 = chis1mart(itypi,itypj)
+ chip1 = chipp1mart(itypi,itypj)
+! chi1=0.0d0
+! chis1=0.0d0
+! chip1=0.0d0
+ chi2=0.0
+ chip2=0.0
+ chis2=0.0
+! chis2 = chis(itypj,itypi)
+ chis12 = chis1 * chis2
+ sig1 = sigmap1mart(itypi,itypj)
+ sig2=0.0
+! sig2 = sigmap2(itypi,itypj)
+! alpha factors from Fcav/Gcav
+ b1cav = alphasurmart(1,itypi,itypj)
+ b2cav = alphasurmart(2,itypi,itypj)
+ b3cav = alphasurmart(3,itypi,itypj)
+ b4cav = alphasurmart(4,itypi,itypj)
+
+! used to determine whether we want to do quadrupole calculations
+ eps_in = epsintabmart(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
+
+ DO k = 1, 3
+ ctail(k,1)=(c(k,i)+c(k,i+1))/2.0
+ ctail(k,2)=c(k,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)
+ do k=1,3
+ Rtail_distance(k) = boxshift(ctail(k,2) - ctail(k,1),locbox(k))
+ enddo
+
+!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 lomartion and distance calculations
+! dhead1
+ d1 = dheadmart(1, 1, itypi, itypj)
+! print *,"d1",d1
+! d1=0.0d0
+! d2 = dhead(2, 1, itypi, itypj)
+ DO k = 1,3
+! lomartion of polar head is computed by taking hydrophobic centre
+! and moving by a d1 * dc_norm vector
+! see unres publimartions 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_mart(itypi,itypj)
+! print *,"ADAM",aa_aq(itypi,itypj)
+
+! c1 = 0.0d0
+ c2 = fac * bb_aq_mart(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
+ gg(2) = fac
+ gg(3) = fac
+
+ 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
+ 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 = dtailmart(1,itypi,itypj) * vbld_inv(i)
+ facd2 = dtailmart(2,itypi,itypj) * vbld_inv(j+nres)
+ DO k = 1, 3
+ pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i))
+! gradpepmartx(k,i) = gradpepmartx(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)
+ gradpepmart(k,i) = gradpepmart(k,i) &
+ - (( dFdR + gg(k) ) * ertail(k))/2.0d0*sss_ele_cut&
+ -(evdwij+Fcav)*rij*sss_ele_grad*rreal(k)*0.5d0
+ gradpepmart(k,i+1) = gradpepmart(k,i+1) &
+ - (( dFdR + gg(k) ) * ertail(k))/2.0d0*sss_ele_cut&
+ -(evdwij+Fcav)*rij*sss_ele_grad*rreal(k)*0.5d0
+
+ gradpepmart(k,j) = gradpepmart(k,j) &
+ + (( dFdR + gg(k) ) * ertail(k))*sss_ele_cut&
+ +(evdwij+Fcav)*rij*sss_ele_grad*rreal(k)
+
+ gg(k) = 0.0d0
+ ENDDO
+!c! Compute head-head and head-tail energies for each state
+!c! Dipole-charge interactions
+ isel = 2+iabs(Qj)
+ if ((itype(j,4).gt.4).and.(itype(j,4).lt.14)) isel=isel+2
+! if (isel.eq.4) isel=0
+ if (isel.le.2) then
+ eheadtail=0.0d0
+ ELSE if (isel.eq.3) then
+ CALL edq_mart_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_mart_pep(ecl)
+ eheadtail=ecl
+! CALL edd_mart_pep
+! eheadtail=0.0d0
+ endif
+ evdw = evdw +( Fcav + eheadtail)*sss_ele_cut
+! 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_mart_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_mart", i,j, gradpepmart(1,nres-1)
+
+ return
+ end subroutine elip_prot
+
+ SUBROUTINE eqq_mart(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 = alphapolmart(itypi,itypj)
+ alphapol2 = alphapolmart2(itypj,itypi)
+!c! Fisocav and Gisocav analytical parameters
+ al1 = alphisomart(1,itypi,itypj)
+ al2 = alphisomart(2,itypi,itypj)
+ al3 = alphisomart(3,itypi,itypj)
+ al4 = alphisomart(4,itypi,itypj)
+ csig = (1.0d0 &
+ / dsqrt(sigiso1mart(itypi, itypj)**2.0d0 &
+ + sigiso2mart(itypi,itypj)**2.0d0))
+!c!
+ pis = sig0headmart(itypi,itypj)
+ eps_head = epsheadmart(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
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+
+ ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq))
+ Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0)
+ debkap=debaykapmart(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
+!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
+!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)))
+!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 = dtailmart(1,itypi,itypj) * vbld_inv(i+nres)
+ facd4 = dtailmart(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))
+ gradpepmartx(k,i) = gradpepmartx(k,i) &
+ +sss_ele_cut*(- 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)-&
+ sss_ele_grad*rij*rreal(k)*(Ecl+Egb+Epol+Fisocav+Elj)
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
+! gradpepmartx(k,j) = gradpepmartx(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
+
+ gradpepmart(k,i) = gradpepmart(k,i) + &
+ sss_ele_cut*(- dGCLdR * erhead(k)&
+ - dGGBdR * erhead(k)&
+ - dGCVdR * erhead(k)&
+ - dPOLdR1 * erhead_tail(k,1)&
+ - dPOLdR2 * erhead_tail(k,2)&
+ - dGLJdR * erhead(k))&
+ - sss_ele_grad*rij*rreal(k)*(Ecl+Egb+Epol+Fisocav+Elj)
+
+
+ gradpepmart(k,j) = gradpepmart(k,j) + &
+ sss_ele_cut*( dGCLdR * erhead(k) &
+ + dGGBdR * erhead(k) &
+ + dGCVdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1) &
+ + dPOLdR2 * erhead_tail(k,2)&
+ + dGLJdR * erhead(k))&
+ +sss_ele_grad*rij*rreal(k)*(Ecl+Egb+Epol+Fisocav+Elj)
+ END DO
+ RETURN
+ END SUBROUTINE eqq_mart
+
+ SUBROUTINE eqd_mart(Ecl,Elj,Epol)
+ use calc_data
+ use comm_momo
+ double precision facd4, federmaus,ecl,elj,epol
+ alphapol1 = alphapolmart(itypi,itypj)
+ w1 = wqdipmart(1,itypi,itypj)
+ w2 = wqdipmart(2,itypi,itypj)
+ pis = sig0headmart(itypi,itypj)
+ eps_head = epsheadmart(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)
+!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
+!c! dPOLdR1 = 0.0d0
+ dPOLdOM1 = 0.0d0
+! dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
+ 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)))
+ 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))
+ gradpepmartx(k,i) = gradpepmartx(k,i) &
+ - dGCLdR * pom&
+ - dPOLdR1 * hawk &
+ - dGLJdR * pom&
+ -(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+
+
+! pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+! gradpepmartx(k,j) = gradpepmartx(k,j) &
+! + dGCLdR * pom &
+! + dPOLdR1 * (erhead_tail(k,1) &
+! -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) &
+! + dGLJdR * pom
+
+
+ gradpepmart(k,i) = gradpepmart(k,i) &
+ - dGCLdR * erhead(k) &
+ - dPOLdR1 * erhead_tail(k,1) &
+ - dGLJdR * erhead(k)&
+ -(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+
+
+ gradpepmart(k,j) = gradpepmart(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1) &
+ + dGLJdR * erhead(k)&
+ +(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+
+
+ END DO
+ RETURN
+ END SUBROUTINE eqd_mart
+
+ SUBROUTINE edq_mart(Ecl,Elj,Epol)
+ use comm_momo
+ use calc_data
+
+ double precision facd3, adler,ecl,elj,epol
+ alphapol2 = alphapolmart(itypi,itypj)
+ w1 = wqdipmart(1,itypi,itypj)
+ w2 = wqdipmart(2,itypi,itypj)
+ pis = sig0headmart(itypi,itypj)
+ eps_head = epsheadmart(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
+!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
+!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) )
+ 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 = dtailmart(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))
+ gradpepmartx(k,i) = gradpepmartx(k,i) &
+ - dGCLdR * pom &
+ - dPOLdR2 * (erhead_tail(k,2) &
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) &
+ - dGLJdR * pom&
+ -(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
+! gradpepmartx(k,j) = gradpepmartx(k,j) &
+! + dGCLdR * pom &
+! + dPOLdR2 * condor &
+! + dGLJdR * pom
+
+
+ gradpepmart(k,i) = gradpepmart(k,i) &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k)&
+ -(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+
+
+ gradpepmart(k,j) = gradpepmart(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR2 * erhead_tail(k,2) &
+ + dGLJdR * erhead(k)&
+ +(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+
+ END DO
+ RETURN
+ END SUBROUTINE edq_mart
+
+ SUBROUTINE edq_mart_pep(Ecl,Elj,Epol)
+ use comm_momo
+ use calc_data
+
+ double precision facd3, adler,ecl,elj,epol
+ alphapol2 = alphapolmart(itypi,itypj)
+ w1 = wqdipmart(1,itypi,itypj)
+ w2 = wqdipmart(2,itypi,itypj)
+ pis = sig0headmart(itypi,itypj)
+ eps_head = epsheadmart(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
+!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
+!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)))
+!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) )
+ facd1 = d1 * vbld_inv(i+1)
+ DO k = 1, 3
+ pom = facd1*(erhead(k)-erdxi*dC_norm(k,i))
+! gradpepmartx(k,i) = gradpepmartx(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))
+! gradpepmartx(k,j) = gradpepmartx(k,j) &
+! + dGCLdR * pom &
+! + dPOLdR2 * condor &
+! + dGLJdR * pom
+
+ gradpepmart(k,i) = gradpepmart(k,i)+pom*(dGCLdR+dGLJdR)
+ gradpepmart(k,i+1) = gradpepmart(k,i+1)-pom*(dGCLdR+dGLJdR)
+
+ gradpepmart(k,i) = gradpepmart(k,i) +0.5d0*( &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k))&
+ -(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+ gradpepmart(k,i+1) = gradpepmart(k,i+1) +0.5d0*( &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k))&
+ -(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+
+
+
+ gradpepmart(k,j) = gradpepmart(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR2 * erhead_tail(k,2) &
+ + dGLJdR * erhead(k)&
+ +(Ecl+Elj+Epol)*sss_ele_grad*rreal(k)*rij
+
+
+ END DO
+ RETURN
+ END SUBROUTINE edq_mart_pep
+!--------------------------------------------------------------------------
+
+ SUBROUTINE edd_mart(ECL)
+! IMPLICIT NONE
+ use comm_momo
+ use calc_data
+
+ double precision ecl
+!c! csig = sigiso(itypi,itypj)
+ w1 = wqdipmart(1,itypi,itypj)
+ w2 = wqdipmart(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 simplifimartion
+ 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
+!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))
+ gradpepmartx(k,i) = gradpepmartx(k,i) - dGCLdR * pom&
+ -ecl*sss_ele_grad*rij*rreal(k)
+! pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+! gradpepmartx(k,j) = gradpepmartx(k,j) + dGCLdR * pom
+
+ gradpepmart(k,i) = gradpepmart(k,i) - dGCLdR * erhead(k)&
+ -ecl*sss_ele_grad*rij*rreal(k)
+
+ gradpepmart(k,j) = gradpepmart(k,j) + dGCLdR * erhead(k)&
+ +ecl*sss_ele_grad*rij*rreal(k)
+
+ END DO
+ RETURN
+ END SUBROUTINE edd_mart
+ SUBROUTINE edd_mart_pep(ECL)
+! IMPLICIT NONE
+ use comm_momo
+ use calc_data
+
+ double precision ecl
+!c! csig = sigiso(itypi,itypj)
+ w1 = wqdipmart(1,itypi,itypj)
+ w2 = wqdipmart(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
+!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))
+ gradpepmart(k,i) = gradpepmart(k,i) + dGCLdR * pom
+ gradpepmart(k,i+1) = gradpepmart(k,i+1) - dGCLdR * pom
+! pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+! gradpepmartx(k,j) = gradpepmartx(k,j) + dGCLdR * pom
+
+ gradpepmart(k,i) = gradpepmart(k,i) - dGCLdR * erhead(k)*0.5d0&
+ -ECL*sss_ele_grad*rreal(k)*rij
+ gradpepmart(k,i+1) = gradpepmart(k,i+1)- dGCLdR * erhead(k)*0.5d0&
+ -ECL*sss_ele_grad*rreal(k)*rij
+
+ gradpepmart(k,j) = gradpepmart(k,j) + dGCLdR * erhead(k)&
+ +ECL*sss_ele_grad*rreal(k)*rij
+
+ END DO
+ RETURN
+ END SUBROUTINE edd_mart_pep
+
+ SUBROUTINE elgrad_init_mart(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,4)
+! print *,"in elegrad",i,j,itypi,itypj
+!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 = sigmamart( itypi,itypj )
+ chi1 = chi1mart( itypi, itypj )
+ chi2 = 0.0d0
+ chi12 = 0.0d0
+ chip1 = chipp1mart( 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 )
+! print *,"before dheadmart"
+!c! distance from center of chain(?) to polar/charged head
+ d1 = dheadmart(1, 1, itypi, itypj)
+ d2 = dheadmart(2, 1, itypi, itypj)
+!c! ai*aj from Fgb
+ a12sq = rborn1mart(itypi,itypj) * rborn2mart(itypi,itypj)
+!c! a12sq = a12sq * a12sq
+!c! charge of amino acid itypi is...
+! print *,"after dheadmart"
+ Qi = icharge(itypi)
+ Qj = ichargelipid(itypj)
+ Qij = Qi * Qj
+! print *,"after icharge"
+
+!c! chis1,2,12
+ chis1 = chis1mart(itypi,itypj)
+ chis2 = 0.0d0
+ chis12 = 0.0d0
+ sig1 = sigmap1mart(itypi,itypj)
+ sig2 = sigmap2mart(itypi,itypj)
+! print *,"before alphasurmart"
+!c! alpha factors from Fcav/Gcav
+ b1cav = alphasurmart(1,itypi,itypj)
+ b2cav = alphasurmart(2,itypi,itypj)
+ b3cav = alphasurmart(3,itypi,itypj)
+ b4cav = alphasurmart(4,itypi,itypj)
+ wqd = wquadmart(itypi, itypj)
+! print *,"after alphasurmar n wquad"
+!c! used by Fgb
+ eps_in = epsintabmart(itypi,itypj)
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+!c!-------------------------------------------------------------------
+!c! tail lomartion and distance calculations
+ Rtail = 0.0d0
+ DO k = 1, 3
+ ctail(k,1)=c(k,i+nres)-dtailmart(1,itypi,itypj)*dc_norm(k,nres+i)
+ ctail(k,2)=c(k,j)!-dtailmart(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 lomartion and distance between polar heads
+!c! distance between heads
+!c! for each one of our three dimensional space...
+ d1 = dheadmart(1, 1, itypi, itypj)
+ d2 = dheadmart(2, 1, itypi, itypj)
+
+ DO k = 1,3
+!c! lomartion of polar head is computed by taking hydrophobic centre
+!c! and moving by a d1 * dc_norm vector
+!c! see unres publimartions 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_mart
+
+ SUBROUTINE elgrad_init_mart_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,4)
+!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 = sigmamart( itypi,itypj )
+ chi1 = chi1mart( itypi, itypj )
+ chi2 = 0.0d0
+ chi12 = 0.0d0
+ chip1 = chipp1mart( 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 = dheadmart(1, 1, itypi, itypj)
+ d2 = dheadmart(2, 1, itypi, itypj)
+!c! ai*aj from Fgb
+ a12sq = rborn1mart(itypi,itypj) * rborn2mart(itypi,itypj)
+!c! a12sq = a12sq * a12sq
+!c! charge of amino acid itypi is...
+ Qi = 0
+ Qj = ichargelipid(itypj)
+! Qij = Qi * Qj
+!c! chis1,2,12
+ chis1 = chis1mart(itypi,itypj)
+ chis2 = 0.0d0
+ chis12 = 0.0d0
+ sig1 = sigmap1mart(itypi,itypj)
+ sig2 = sigmap2mart(itypi,itypj)
+!c! alpha factors from Fcav/Gcav
+ b1cav = alphasurmart(1,itypi,itypj)
+ b2cav = alphasurmart(2,itypi,itypj)
+ b3cav = alphasurmart(3,itypi,itypj)
+ b4cav = alphasurmart(4,itypi,itypj)
+ wqd = wquadmart(itypi, itypj)
+!c! used by Fgb
+ eps_in = epsintabmart(itypi,itypj)
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+!c!-------------------------------------------------------------------
+!c! tail lomartion and distance calculations
+ Rtail = 0.0d0
+ DO k = 1, 3
+ ctail(k,1)=(c(k,i)+c(k,i+1))/2.0-dtailmart(1,itypi,itypj)*dc_norm(k,i)
+ ctail(k,2)=c(k,j)!-dtailmart(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 lomartion and distance between polar heads
+!c! distance between heads
+!c! for each one of our three dimensional space...
+ d1 = dheadmart(1, 1, itypi, itypj)
+ d2 = dheadmart(2, 1, itypi, itypj)
+
+ DO k = 1,3
+!c! lomartion of polar head is computed by taking hydrophobic centre
+!c! and moving by a d1 * dc_norm vector
+!c! see unres publimartions 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_mart_pep
+
+ subroutine sc_grad_mart
+ use calc_data
+ real(kind=8), dimension(3) :: dcosom1,dcosom2
+ 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
+! diagnostics only
+! eom1=0.0d0
+! eom2=0.0d0
+! eom12=evdwij*eps1_om12
+! end diagnostics
+
+ 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))
+ enddo
+ do k=1,3
+ gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))
+! print *,'gg',k,gg(k)
+ enddo
+! print *,i,j,gg_lipi(3),gg_lipj(3),sss_ele_cut
+! write (iout,*) "gg",(gg(k),k=1,3)
+ do k=1,3
+ gradpepmartx(k,i)=gradpepmartx(k,i)-gg(k)*sss_ele_cut &
+ +(eom12*(dc_norm(k,j)-om12*dc_norm(k,nres+i)) &
+ +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*sss_ele_cut
+
+! gradpepcatx(k,j)=gradpepcatx(k,j)+gg(k) &
+! +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,j)) &
+! +eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv
+
+! write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
+! +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+! write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
+! +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ enddo
+!
+! Calculate the components of the gradient in DC and X
+!
+ do l=1,3
+ gradpepmart(l,i)=gradpepmart(l,i)-gg(l)*sss_ele_cut
+ gradpepmart(l,j)=gradpepmart(l,j)+gg(l)*sss_ele_cut
+ enddo
+ end subroutine sc_grad_mart
+
+ subroutine sc_grad_mart_pep
+ use calc_data
+ real(kind=8), dimension(3) :: dcosom1,dcosom2
+ 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
+! diagnostics only
+! eom1=0.0d0
+! eom2=0.0d0
+! eom12=evdwij*eps1_om12
+! end diagnostics
+! write (iout,*) "gg",(gg(k),k=1,3)
+
+ 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)
+ gradpepmart(k,i)= gradpepmart(k,i) +sss_ele_cut*(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)
+ gradpepmart(k,i+1)= gradpepmart(k,i+1) +sss_ele_cut*(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)
+ gradpepmart(k,j)=gradpepmart(k,j)+gg(k)*sss_ele_cut
+ enddo
+ end subroutine sc_grad_mart_pep