+ subroutine ebend_nucl(etheta_nucl)
+ real(kind=8),dimension(nntheterm_nucl+1) :: coskt,sinkt !mmaxtheterm
+ real(kind=8),dimension(nsingle_nucl+1) :: cosph1,sinph1,cosph2,sinph2 !maxsingle
+ real(kind=8),dimension(ndouble_nucl+1,ndouble_nucl+1) :: cosph1ph2,sinph1ph2 !maxdouble,maxdouble
+ logical :: lprn=.false., lprn1=.false.
+!el local variables
+ integer :: i,k,iblock,ityp1,ityp2,ityp3,l,m
+ real(kind=8) :: dethetai,dephii,dephii1,theti2,phii,phii1,ethetai
+ real(kind=8) :: aux,etheta_nucl,ccl,ssl,scl,csl,ethetacnstr
+! local variables for constrains
+ real(kind=8) :: difi,thetiii
+ integer itheta
+ etheta_nucl=0.0D0
+! print *,"ithet_start",ithet_nucl_start," ithet_end",ithet_nucl_end,nres
+ do i=ithet_nucl_start,ithet_nucl_end
+ if ((itype(i-1,2).eq.ntyp1_molec(2)).or.&
+ (itype(i-2,2).eq.ntyp1_molec(2)).or. &
+ (itype(i,2).eq.ntyp1_molec(2))) cycle
+ dethetai=0.0d0
+ dephii=0.0d0
+ dephii1=0.0d0
+ theti2=0.5d0*theta(i)
+ ityp2=ithetyp_nucl(itype(i-1,2))
+ do k=1,nntheterm_nucl
+ coskt(k)=dcos(k*theti2)
+ sinkt(k)=dsin(k*theti2)
+ enddo
+ if (i.gt.3 .and. itype(i-2,2).ne.ntyp1_molec(2)) then
+#ifdef OSF
+ phii=phi(i)
+ if (phii.ne.phii) phii=150.0
+#else
+ phii=phi(i)
+#endif
+ ityp1=ithetyp_nucl(itype(i-2,2))
+ do k=1,nsingle_nucl
+ cosph1(k)=dcos(k*phii)
+ sinph1(k)=dsin(k*phii)
+ enddo
+ else
+ phii=0.0d0
+ ityp1=nthetyp_nucl+1
+ do k=1,nsingle_nucl
+ cosph1(k)=0.0d0
+ sinph1(k)=0.0d0
+ enddo
+ endif
+
+ if (i.lt.nres .and. itype(i,2).ne.ntyp1_molec(2)) then
+#ifdef OSF
+ phii1=phi(i+1)
+ if (phii1.ne.phii1) phii1=150.0
+ phii1=pinorm(phii1)
+#else
+ phii1=phi(i+1)
+#endif
+ ityp3=ithetyp_nucl(itype(i,2))
+ do k=1,nsingle_nucl
+ cosph2(k)=dcos(k*phii1)
+ sinph2(k)=dsin(k*phii1)
+ enddo
+ else
+ phii1=0.0d0
+ ityp3=nthetyp_nucl+1
+ do k=1,nsingle_nucl
+ cosph2(k)=0.0d0
+ sinph2(k)=0.0d0
+ enddo
+ endif
+ ethetai=aa0thet_nucl(ityp1,ityp2,ityp3)
+ do k=1,ndouble_nucl
+ do l=1,k-1
+ ccl=cosph1(l)*cosph2(k-l)
+ ssl=sinph1(l)*sinph2(k-l)
+ scl=sinph1(l)*cosph2(k-l)
+ csl=cosph1(l)*sinph2(k-l)
+ cosph1ph2(l,k)=ccl-ssl
+ cosph1ph2(k,l)=ccl+ssl
+ sinph1ph2(l,k)=scl+csl
+ sinph1ph2(k,l)=scl-csl
+ enddo
+ enddo
+ if (lprn) then
+ write (iout,*) "i",i," ityp1",ityp1," ityp2",ityp2,&
+ " ityp3",ityp3," theti2",theti2," phii",phii," phii1",phii1
+ write (iout,*) "coskt and sinkt",nntheterm_nucl
+ do k=1,nntheterm_nucl
+ write (iout,*) k,coskt(k),sinkt(k)
+ enddo
+ endif
+ do k=1,ntheterm_nucl
+ ethetai=ethetai+aathet_nucl(k,ityp1,ityp2,ityp3)*sinkt(k)
+ dethetai=dethetai+0.5d0*k*aathet_nucl(k,ityp1,ityp2,ityp3)&
+ *coskt(k)
+ if (lprn)&
+ write (iout,*) "k",k," aathet",aathet_nucl(k,ityp1,ityp2,ityp3),&
+ " ethetai",ethetai
+ enddo
+ if (lprn) then
+ write (iout,*) "cosph and sinph"
+ do k=1,nsingle_nucl
+ write (iout,*) k,cosph1(k),sinph1(k),cosph2(k),sinph2(k)
+ enddo
+ write (iout,*) "cosph1ph2 and sinph2ph2"
+ do k=2,ndouble_nucl
+ do l=1,k-1
+ write (iout,*) l,k,cosph1ph2(l,k),cosph1ph2(k,l),&
+ sinph1ph2(l,k),sinph1ph2(k,l)
+ enddo
+ enddo
+ write(iout,*) "ethetai",ethetai
+ endif
+ do m=1,ntheterm2_nucl
+ do k=1,nsingle_nucl
+ aux=bbthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph1(k)&
+ +ccthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph1(k)&
+ +ddthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph2(k)&
+ +eethet_nucl(k,m,ityp1,ityp2,ityp3)*sinph2(k)
+ ethetai=ethetai+sinkt(m)*aux
+ dethetai=dethetai+0.5d0*m*aux*coskt(m)
+ dephii=dephii+k*sinkt(m)*(&
+ ccthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph1(k)-&
+ bbthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph1(k))
+ dephii1=dephii1+k*sinkt(m)*(&
+ eethet_nucl(k,m,ityp1,ityp2,ityp3)*cosph2(k)-&
+ ddthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph2(k))
+ if (lprn) &
+ write (iout,*) "m",m," k",k," bbthet",&
+ bbthet_nucl(k,m,ityp1,ityp2,ityp3)," ccthet",&
+ ccthet_nucl(k,m,ityp1,ityp2,ityp3)," ddthet",&
+ ddthet_nucl(k,m,ityp1,ityp2,ityp3)," eethet",&
+ eethet_nucl(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+ enddo
+ enddo
+ if (lprn) &
+ write(iout,*) "ethetai",ethetai
+ do m=1,ntheterm3_nucl
+ do k=2,ndouble_nucl
+ do l=1,k-1
+ aux=ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+&
+ ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+&
+ ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+&
+ ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)
+ ethetai=ethetai+sinkt(m)*aux
+ dethetai=dethetai+0.5d0*m*coskt(m)*aux
+ dephii=dephii+l*sinkt(m)*(&
+ -ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)-&
+ ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+&
+ ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+&
+ ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+ dephii1=dephii1+(k-l)*sinkt(m)*( &
+ -ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+&
+ ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+&
+ ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)-&
+ ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+ if (lprn) then
+ write (iout,*) "m",m," k",k," l",l," ffthet", &
+ ffthet_nucl(l,k,m,ityp1,ityp2,ityp3), &
+ ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)," ggthet",&
+ ggthet_nucl(l,k,m,ityp1,ityp2,ityp3),&
+ ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+ write (iout,*) cosph1ph2(l,k)*sinkt(m), &
+ cosph1ph2(k,l)*sinkt(m),&
+ sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m)
+ endif
+ enddo
+ enddo
+ enddo
+10 continue
+ if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') &
+ i,theta(i)*rad2deg,phii*rad2deg, &
+ phii1*rad2deg,ethetai
+ etheta_nucl=etheta_nucl+ethetai
+! print *,i,"partial sum",etheta_nucl
+ if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang_nucl*dephii
+ if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang_nucl*dephii1
+ gloc(nphi+i-2,icg)=wang_nucl*dethetai
+ enddo
+ return
+ end subroutine ebend_nucl
+!----------------------------------------------------
+ subroutine etor_nucl(etors_nucl)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.VAR'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCAL'
+! include 'COMMON.TORSION'
+! include 'COMMON.INTERACT'
+! include 'COMMON.DERIV'
+! include 'COMMON.CHAIN'
+! include 'COMMON.NAMES'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.FFIELD'
+! include 'COMMON.TORCNSTR'
+! include 'COMMON.CONTROL'
+ real(kind=8) :: etors_nucl,edihcnstr
+ logical :: lprn
+!el local variables
+ integer :: i,j,iblock,itori,itori1
+ real(kind=8) :: phii,gloci,v1ij,v2ij,cosphi,sinphi,&
+ vl1ij,vl2ij,vl3ij,pom1,difi,etors_ii,pom
+! Set lprn=.true. for debugging
+ lprn=.false.
+! lprn=.true.
+ etors_nucl=0.0D0
+! print *,"iphi_nucl_start/end", iphi_nucl_start,iphi_nucl_end
+ do i=iphi_nucl_start,iphi_nucl_end
+ if (itype(i-2,2).eq.ntyp1_molec(2) .or. itype(i-1,2).eq.ntyp1_molec(2) &
+ .or. itype(i-3,2).eq.ntyp1_molec(2) &
+ .or. itype(i,2).eq.ntyp1_molec(2)) cycle
+ etors_ii=0.0D0
+ itori=itortyp_nucl(itype(i-2,2))
+ itori1=itortyp_nucl(itype(i-1,2))
+ phii=phi(i)
+! print *,i,itori,itori1
+ gloci=0.0D0
+!C Regular cosine and sine terms
+ do j=1,nterm_nucl(itori,itori1)
+ v1ij=v1_nucl(j,itori,itori1)
+ v2ij=v2_nucl(j,itori,itori1)
+ cosphi=dcos(j*phii)
+ sinphi=dsin(j*phii)
+ etors_nucl=etors_nucl+v1ij*cosphi+v2ij*sinphi
+ if (energy_dec) etors_ii=etors_ii+&
+ v1ij*cosphi+v2ij*sinphi
+ gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
+ enddo
+!C Lorentz terms
+!C v1
+!C E = SUM ----------------------------------- - v1
+!C [v2 cos(phi/2)+v3 sin(phi/2)]^2 + 1
+!C
+ cosphi=dcos(0.5d0*phii)
+ sinphi=dsin(0.5d0*phii)
+ do j=1,nlor_nucl(itori,itori1)
+ vl1ij=vlor1_nucl(j,itori,itori1)
+ vl2ij=vlor2_nucl(j,itori,itori1)
+ vl3ij=vlor3_nucl(j,itori,itori1)
+ pom=vl2ij*cosphi+vl3ij*sinphi
+ pom1=1.0d0/(pom*pom+1.0d0)
+ etors_nucl=etors_nucl+vl1ij*pom1
+ if (energy_dec) etors_ii=etors_ii+ &
+ vl1ij*pom1
+ pom=-pom*pom1*pom1
+ gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom
+ enddo
+!C Subtract the constant term
+ etors_nucl=etors_nucl-v0_nucl(itori,itori1)
+ if (energy_dec) write (iout,'(a6,i5,0pf7.3)') &
+ 'etor',i,etors_ii-v0_nucl(itori,itori1)
+ if (lprn) &
+ write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') &
+ restyp(itype(i-2,2),2),i-2,restyp(itype(i-1,2),2),i-1,itori,itori1, &
+ (v1_nucl(j,itori,itori1),j=1,6),(v2_nucl(j,itori,itori1),j=1,6)
+ gloc(i-3,icg)=gloc(i-3,icg)+wtor_nucl*gloci
+!c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
+ enddo
+ return
+ end subroutine etor_nucl
+!------------------------------------------------------------
+ subroutine epp_nucl_sub(evdw1,ees)
+!C
+!C This subroutine calculates the average interaction energy and its gradient
+!C in the virtual-bond vectors between non-adjacent peptide groups, based on
+!C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
+!C The potential depends both on the distance of peptide-group centers and on
+!C the orientation of the CA-CA virtual bonds.
+!C
+ integer :: i,j,k,iteli,itelj,num_conti,isubchap,ind
+ real(kind=8) :: dxi,dyi,dzi,dxj,dyj,dzj,aaa,bbb
+ real(kind=8) :: xj,yj,zj,rij,rrmij,sss,r3ij,r6ij,evdw1,&
+ dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,&
+ dx_normj,dy_normj,dz_normj,rmij,ev1,ev2,evdwij,facvdw
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,sss_grad,fac,evdw1ij
+ integer xshift,yshift,zshift
+ real(kind=8),dimension(3):: ggg,gggp,gggm,erij
+ real(kind=8) :: ees,eesij
+!c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
+ real(kind=8) scal_el /0.5d0/
+ t_eelecij=0.0d0
+ ees=0.0D0
+ evdw1=0.0D0
+ ind=0
+!c
+!c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
+!c
+! print *,"iatel_s_nucl,iatel_e_nucl",iatel_s_nucl,iatel_e_nucl
+ do i=iatel_s_nucl,iatel_e_nucl
+ if (itype(i,2).eq.ntyp1_molec(2) .or. itype(i+1,2).eq.ntyp1_molec(2)) cycle
+ dxi=dc(1,i)
+ dyi=dc(2,i)
+ dzi=dc(3,i)
+ dx_normi=dc_norm(1,i)
+ dy_normi=dc_norm(2,i)
+ dz_normi=dc_norm(3,i)
+ xmedi=c(1,i)+0.5d0*dxi
+ ymedi=c(2,i)+0.5d0*dyi
+ zmedi=c(3,i)+0.5d0*dzi
+ xmedi=dmod(xmedi,boxxsize)
+ if (xmedi.lt.0) xmedi=xmedi+boxxsize
+ ymedi=dmod(ymedi,boxysize)
+ if (ymedi.lt.0) ymedi=ymedi+boxysize
+ zmedi=dmod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
+
+ do j=ielstart_nucl(i),ielend_nucl(i)
+ if (itype(j,2).eq.ntyp1_molec(2) .or. itype(j+1,2).eq.ntyp1_molec(2)) cycle
+ ind=ind+1
+ dxj=dc(1,j)
+ dyj=dc(2,j)
+ dzj=dc(3,j)
+! xj=c(1,j)+0.5D0*dxj-xmedi
+! yj=c(2,j)+0.5D0*dyj-ymedi
+! zj=c(3,j)+0.5D0*dzj-zmedi
+ xj=c(1,j)+0.5D0*dxj
+ yj=c(2,j)+0.5D0*dyj
+ zj=c(3,j)+0.5D0*dzj
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ isubchap=0
+ dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ isubchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (isubchap.eq.1) then
+!C print *,i,j
+ xj=xj_temp-xmedi
+ yj=yj_temp-ymedi
+ zj=zj_temp-zmedi
+ else
+ xj=xj_safe-xmedi
+ yj=yj_safe-ymedi
+ zj=zj_safe-zmedi
+ endif
+
+ rij=xj*xj+yj*yj+zj*zj
+!c write (2,*)"ij",i,j," r0pp",r0pp," rij",rij," epspp",epspp
+ fac=(r0pp**2/rij)**3
+ ev1=epspp*fac*fac
+ ev2=epspp*fac
+ evdw1ij=ev1-2*ev2
+ fac=(-ev1-evdw1ij)/rij
+! write (2,*)"fac",fac," ev1",ev1," ev2",ev2," evdw1ij",evdw1ij
+ if (energy_dec) write(iout,'(2i5,a9,f10.4)') i,j,"evdw1ij",evdw1ij
+ evdw1=evdw1+evdw1ij
+!C
+!C Calculate contributions to the Cartesian gradient.
+!C
+ ggg(1)=fac*xj
+ ggg(2)=fac*yj
+ ggg(3)=fac*zj
+ do k=1,3
+ gvdwpp_nucl(k,i)=gvdwpp_nucl(k,i)-ggg(k)
+ gvdwpp_nucl(k,j)=gvdwpp_nucl(k,j)+ggg(k)
+ enddo
+!c phoshate-phosphate electrostatic interactions
+ rij=dsqrt(rij)
+ fac=1.0d0/rij
+ eesij=dexp(-BEES*rij)*fac
+! write (2,*)"fac",fac," eesijpp",eesij
+ if (energy_dec) write(iout,'(2i5,a9,f10.4)') i,j,"eesijpp",eesij
+ ees=ees+eesij
+!c fac=-eesij*fac
+ fac=-(fac+BEES)*eesij*fac
+ ggg(1)=fac*xj
+ ggg(2)=fac*yj
+ ggg(3)=fac*zj
+!c write(2,*) "ggg",i,j,ggg(1),ggg(2),ggg(3)
+!c write(2,*) "gelpp",i,(gelpp(k,i),k=1,3)
+!c write(2,*) "gelpp",j,(gelpp(k,j),k=1,3)
+ do k=1,3
+ gelpp(k,i)=gelpp(k,i)-ggg(k)
+ gelpp(k,j)=gelpp(k,j)+ggg(k)
+ enddo
+ enddo ! j
+ enddo ! i
+!c ees=332.0d0*ees
+ ees=AEES*ees
+ do i=nnt,nct
+!c write (2,*) "i",i," gelpp",(gelpp(k,i),k=1,3)
+ do k=1,3
+ gvdwpp_nucl(k,i)=6*gvdwpp_nucl(k,i)
+!c gelpp(k,i)=332.0d0*gelpp(k,i)
+ gelpp(k,i)=AEES*gelpp(k,i)
+ enddo
+!c write (2,*) "i",i," gelpp",(gelpp(k,i),k=1,3)
+ enddo
+!c write (2,*) "total EES",ees
+ return
+ end subroutine epp_nucl_sub
+!---------------------------------------------------------------------
+ subroutine epsb(evdwpsb,eelpsb)
+! use comm_locel
+!C
+!C This subroutine calculates the excluded-volume interaction energy between
+!C peptide-group centers and side chains and its gradient in virtual-bond and
+!C side-chain vectors.
+!C
+ real(kind=8),dimension(3):: ggg
+ integer :: i,iint,j,k,iteli,itypj,subchap
+ real(kind=8) :: evdw2,evdw2_14,xi,yi,zi,xj,yj,zj,rrij,fac,&
+ e1,e2,evdwij,rij,evdwpsb,eelpsb
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
+ integer xshift,yshift,zshift
+
+!cd print '(a)','Enter ESCP'
+!cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
+ eelpsb=0.0d0
+ evdwpsb=0.0d0
+! print *,"iatscp_s_nucl,iatscp_e_nucl",iatscp_s_nucl,iatscp_e_nucl
+ do i=iatscp_s_nucl,iatscp_e_nucl
+ if (itype(i,2).eq.ntyp1_molec(2) &
+ .or. itype(i+1,2).eq.ntyp1_molec(2)) cycle
+ xi=0.5D0*(c(1,i)+c(1,i+1))
+ yi=0.5D0*(c(2,i)+c(2,i+1))
+ zi=0.5D0*(c(3,i)+c(3,i+1))
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+
+ do iint=1,nscp_gr_nucl(i)
+
+ do j=iscpstart_nucl(i,iint),iscpend_nucl(i,iint)
+ itypj=itype(j,2)
+ if (itypj.eq.ntyp1_molec(2)) cycle
+!C Uncomment following three lines for SC-p interactions
+!c xj=c(1,nres+j)-xi
+!c yj=c(2,nres+j)-yi
+!c zj=c(3,nres+j)-zi
+!C Uncomment following three lines for Ca-p interactions
+! xj=c(1,j)-xi
+! yj=c(2,j)-yi
+! zj=c(3,j)-zi
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ fac=rrij**expon2
+ e1=fac*fac*aad_nucl(itypj)
+ e2=fac*bad_nucl(itypj)
+ if (iabs(j-i) .le. 2) then
+ e1=scal14*e1
+ e2=scal14*e2
+ endif
+ evdwij=e1+e2
+ evdwpsb=evdwpsb+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a4)') &
+ 'evdw2',i,j,evdwij,"tu4"
+!C
+!C Calculate contributions to the gradient in the virtual-bond and SC vectors.
+!C
+ fac=-(evdwij+e1)*rrij
+ ggg(1)=xj*fac
+ ggg(2)=yj*fac
+ ggg(3)=zj*fac
+ do k=1,3
+ gvdwpsb1(k,i)=gvdwpsb1(k,i)-ggg(k)
+ gvdwpsb(k,j)=gvdwpsb(k,j)+ggg(k)
+ enddo
+ enddo
+
+ enddo ! iint
+ enddo ! i
+ do i=1,nct
+ do j=1,3
+ gvdwpsb(j,i)=expon*gvdwpsb(j,i)
+ gvdwpsb1(j,i)=expon*gvdwpsb1(j,i)
+ enddo
+ enddo
+ return
+ end subroutine epsb
+
+!------------------------------------------------------
+ subroutine esb_gb(evdwsb,eelsb)
+ use comm_locel
+ use calc_data_nucl
+ integer :: iint,itypi,itypi1,itypj,subchap,num_conti2
+ real(kind=8) :: xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
+ real(kind=8) :: evdw,sig0iji,evdwsb,eelsb,ecorr,eelij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,aa,bb,faclip,sig0ij
+ integer :: ii
+ logical lprn
+ evdw=0.0D0
+ eelsb=0.0d0
+ ecorr=0.0d0
+ evdwsb=0.0D0
+ lprn=.false.
+ ind=0
+! print *,"iastsc_nucl",iatsc_s_nucl,iatsc_e_nucl
+ do i=iatsc_s_nucl,iatsc_e_nucl
+ num_conti=0
+ num_conti2=0
+ itypi=itype(i,2)
+! PRINT *,"I=",i,itypi
+ if (itypi.eq.ntyp1_molec(2)) cycle
+ itypi1=itype(i+1,2)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ xi=dmod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=dmod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=dmod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+ dsci_inv=vbld_inv(i+nres)
+!C
+!C Calculate SC interaction energy.
+!C
+ do iint=1,nint_gr_nucl(i)
+! print *,"tu?",i,istart_nucl(i,iint),iend_nucl(i,iint)
+ do j=istart_nucl(i,iint),iend_nucl(i,iint)
+ ind=ind+1
+! print *,"JESTEM"
+ itypj=itype(j,2)
+ if (itypj.eq.ntyp1_molec(2)) cycle
+ dscj_inv=vbld_inv(j+nres)
+ sig0ij=sigma_nucl(itypi,itypj)
+ chi1=chi_nucl(itypi,itypj)
+ chi2=chi_nucl(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip_nucl(itypi,itypj)
+ chip2=chip_nucl(itypj,itypi)
+ chip12=chip1*chip2
+! xj=c(1,nres+j)-xi
+! yj=c(2,nres+j)-yi
+! zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+!C Calculate angle-dependent terms of energy and contributions to their
+!C derivatives.
+ erij(1)=xj*rij
+ erij(2)=yj*rij
+ erij(3)=zj*rij
+ om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
+ om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
+ om12=dxi*dxj+dyi*dyj+dzi*dzj
+ call sc_angular_nucl
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+sig0ij
+! print *,rij_shift,"rij_shift"
+!c write (2,*) " rij",1.0D0/rij," sig",sig," sig0ij",sig0ij,
+!c & " rij_shift",rij_shift
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D20
+ return
+ endif
+ sigder=-sig*sigsq
+!c---------------------------------------------------------------
+ rij_shift=1.0D0/rij_shift
+ fac=rij_shift**expon
+ e1=fac*fac*aa_nucl(itypi,itypj)
+ e2=fac*bb_nucl(itypi,itypj)
+ evdwij=eps1*eps2rt*(e1+e2)
+!c write (2,*) "eps1",eps1," eps2rt",eps2rt,
+!c & " e1",e1," e2",e2," evdwij",evdwij
+ eps2der=evdwij
+ evdwij=evdwij*eps2rt
+ evdwsb=evdwsb+evdwij
+ if (lprn) then
+ sigm=dabs(aa_nucl(itypi,itypj)/bb_nucl(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb_nucl(itypi,itypj)**2/aa_nucl(itypi,itypj)
+ write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
+ restyp(itypi,2),i,restyp(itypj,2),j, &
+ epsi,sigm,chi1,chi2,chip1,chip2, &
+ eps1,eps2rt**2,sig,sig0ij, &
+ om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,&
+ evdwij
+ write (iout,*) "aa",aa_nucl(itypi,itypj)," bb",bb_nucl(itypi,itypj)
+ endif
+
+ if (energy_dec) write (iout,'(a6,2i5,e15.3,a4)') &
+ 'evdw',i,j,evdwij,"tu3"
+
+
+!C Calculate gradient components.
+ e1=e1*eps1*eps2rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ fac=rij*fac
+!c fac=0.0d0
+!C Calculate the radial part of the gradient
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+!C Calculate angular part of the gradient.
+ call sc_grad_nucl
+ call eelsbij(eelij,num_conti2)
+ if (energy_dec .and. &
+ (j.eq.i+1.or.j.eq.nres-i+1.or.j.eq.nres-i.or.j.eq.nres-i+2)) &
+ write (istat,'(e14.5)') evdwij
+ eelsb=eelsb+eelij
+ enddo ! j
+ enddo ! iint
+ num_cont_hb(i)=num_conti2
+ enddo ! i
+!c write (iout,*) "Number of loop steps in EGB:",ind
+!cccc energy_dec=.false.
+ return
+ end subroutine esb_gb
+!-------------------------------------------------------------------------------
+ subroutine eelsbij(eesij,num_conti2)
+ use comm_locel
+ use calc_data_nucl
+ real(kind=8),dimension(3) :: ggg,gggp,gggm,dcosb,dcosg
+ real(kind=8),dimension(3,3) :: erder,uryg,urzg,vryg,vrzg
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,rlocshield,fracinbuf
+ integer xshift,yshift,zshift,ilist,iresshield,num_conti2
+
+!c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
+ real(kind=8) scal_el /0.5d0/
+ integer :: iteli,itelj,kkk,kkll,m,isubchap
+ real(kind=8) :: ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp,facfac
+ real(kind=8) :: ees,evdw1,eel_loc,aaa,bbb,ael3i,ael63i,ael32i
+ real(kind=8) :: dx_normj,dy_normj,dz_normj,&
+ r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,fac5,fac6,&
+ el1,el2,el3,el4,eesij,ees0ij,facvdw,facel,fac1,ecosa,&
+ ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der,&
+ a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij,&
+ ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp,&
+ ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,&
+ ecosgp,ecosam,ecosbm,ecosgm,ghalf,itypi,itypj
+ ind=ind+1
+ itypi=itype(i,2)
+ itypj=itype(j,2)
+! print *,i,j,itypi,itypj,istype(i),istype(j),"????"
+ ael6i=ael6_nucl(itypi,itypj)
+ ael3i=ael3_nucl(itypi,itypj)
+ ael63i=ael63_nucl(itypi,itypj)
+ ael32i=ael32_nucl(itypi,itypj)
+!c write (iout,*) "eelecij",i,j,itype(i),itype(j),
+!c & ael6i,ael3i,ael63i,al32i,rij,rrij
+ dxj=dc(1,j+nres)
+ dyj=dc(2,j+nres)
+ dzj=dc(3,j+nres)
+ dx_normi=dc_norm(1,i+nres)
+ dy_normi=dc_norm(2,i+nres)
+ dz_normi=dc_norm(3,i+nres)
+ dx_normj=dc_norm(1,j+nres)
+ dy_normj=dc_norm(2,j+nres)
+ dz_normj=dc_norm(3,j+nres)
+!c xj=c(1,j)+0.5D0*dxj-xmedi
+!c yj=c(2,j)+0.5D0*dyj-ymedi
+!c zj=c(3,j)+0.5D0*dzj-zmedi
+ if (ipot_nucl.ne.2) then
+ cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
+ cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
+ cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
+ else
+ cosa=om12
+ cosb=om1
+ cosg=om2
+ endif
+ r3ij=rij*rrij
+ r6ij=r3ij*r3ij
+ fac=cosa-3.0D0*cosb*cosg
+ facfac=fac*fac
+ fac1=3.0d0*(cosb*cosb+cosg*cosg)
+ fac3=ael6i*r6ij
+ fac4=ael3i*r3ij
+ fac5=ael63i*r6ij
+ fac6=ael32i*r6ij
+!c write (iout,*) "r3ij",r3ij," r6ij",r6ij," fac",fac," fac1",fac1,
+!c & " fac2",fac2," fac3",fac3," fac4",fac4," fac5",fac5," fac6",fac6
+ el1=fac3*(4.0D0+facfac-fac1)
+ el2=fac4*fac
+ el3=fac5*(2.0d0-2.0d0*facfac+fac1)
+ el4=fac6*facfac
+ eesij=el1+el2+el3+el4
+!C 12/26/95 - for the evaluation of multi-body H-bonding interactions
+ ees0ij=4.0D0+facfac-fac1
+
+ if (energy_dec) then
+ if(j.eq.i+1.or.j.eq.nres-i+1.or.j.eq.nres-i.or.j.eq.nres-i+2) &
+ write (istat,'(2a1,i4,1x,2a1,i4,4f10.5,3e12.5,$)') &
+ sugartyp(istype(i)),restyp(itypi,2),i,sugartyp(istype(j)),&
+ restyp(itypj,2),j,1.0d0/rij,cosa,cosb,cosg,fac*r3ij, &
+ (4.0D0+facfac-fac1)*r6ij,(2.0d0-2.0d0*facfac+fac1)*r6ij
+ write (iout,'(a6,2i5,e15.3)') 'ees',i,j,eesij
+ endif
+
+!C
+!C Calculate contributions to the Cartesian gradient.
+!C
+ facel=-3.0d0*rrij*(eesij+el1+el3+el4)
+ fac1=fac
+!c erij(1)=xj*rmij
+!c erij(2)=yj*rmij
+!c erij(3)=zj*rmij
+!*
+!* Radial derivatives. First process both termini of the fragment (i,j)
+!*
+ ggg(1)=facel*xj
+ ggg(2)=facel*yj
+ ggg(3)=facel*zj
+ do k=1,3
+ gelsbc(k,j)=gelsbc(k,j)+ggg(k)
+ gelsbc(k,i)=gelsbc(k,i)-ggg(k)
+ gelsbx(k,j)=gelsbx(k,j)+ggg(k)
+ gelsbx(k,i)=gelsbx(k,i)-ggg(k)
+ enddo
+!*
+!* Angular part
+!*
+ ecosa=2.0D0*fac3*fac1+fac4+(-4.0d0*fac5+2.0d0*fac6)*fac1
+ fac4=-3.0D0*fac4
+ fac3=-6.0D0*fac3
+ fac5= 6.0d0*fac5
+ fac6=-6.0d0*fac6
+ ecosb=fac3*(fac1*cosg+cosb)+cosg*fac4+(cosb+2*fac1*cosg)*fac5+&
+ fac6*fac1*cosg
+ ecosg=fac3*(fac1*cosb+cosg)+cosb*fac4+(cosg+2*fac1*cosb)*fac5+&
+ fac6*fac1*cosb
+ do k=1,3
+ dcosb(k)=rij*(dc_norm(k,i+nres)-erij(k)*cosb)
+ dcosg(k)=rij*(dc_norm(k,j+nres)-erij(k)*cosg)
+ enddo
+ do k=1,3
+ ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
+ enddo
+ do k=1,3
+ gelsbx(k,i)=gelsbx(k,i)-ggg(k) &
+ +(ecosa*(dc_norm(k,j+nres)-cosa*dc_norm(k,i+nres))&
+ + ecosb*(erij(k)-cosb*dc_norm(k,i+nres)))*vbld_inv(i+nres)
+ gelsbx(k,j)=gelsbx(k,j)+ggg(k) &
+ +(ecosa*(dc_norm(k,i+nres)-cosa*dc_norm(k,j+nres))&
+ + ecosg*(erij(k)-cosg*dc_norm(k,j+nres)))*vbld_inv(j+nres)
+ gelsbc(k,j)=gelsbc(k,j)+ggg(k)
+ gelsbc(k,i)=gelsbc(k,i)-ggg(k)
+ enddo
+! IF ( (wcorr_nucl.gt.0.0d0.or.wcorr3_nucl.gt.0.0d0) .and.
+ IF ( j.gt.i+1 .and.&
+ num_conti.le.maxconts) THEN
+!C
+!C Calculate the contact function. The ith column of the array JCONT will
+!C contain the numbers of atoms that make contacts with the atom I (of numbers
+!C greater than I). The arrays FACONT and GACONT will contain the values of
+!C the contact function and its derivative.
+ r0ij=2.20D0*sigma(itypi,itypj)
+!c write (2,*) "ij",i,j," rij",1.0d0/rij," r0ij",r0ij
+ call gcont(rij,r0ij,1.0D0,0.2d0/r0ij,fcont,fprimcont)
+!c write (2,*) "fcont",fcont
+ if (fcont.gt.0.0D0) then
+ num_conti=num_conti+1
+ num_conti2=num_conti2+1
+
+ if (num_conti.gt.maxconts) then
+ write (iout,*) 'WARNING - max. # of contacts exceeded;',&
+ ' will skip next contacts for this conf.'
+ else
+ jcont_hb(num_conti,i)=j
+!c write (iout,*) "num_conti",num_conti,
+!c & " jcont_hb",jcont_hb(num_conti,i)
+!C Calculate contact energies
+ cosa4=4.0D0*cosa
+ wij=cosa-3.0D0*cosb*cosg
+ cosbg1=cosb+cosg
+ cosbg2=cosb-cosg
+ fac3=dsqrt(-ael6i)*r3ij
+!c write (2,*) "ael6i",ael6i," r3ij",r3ij," fac3",fac3
+ ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
+ if (ees0tmp.gt.0) then
+ ees0pij=dsqrt(ees0tmp)
+ else
+ ees0pij=0
+ endif
+ ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
+ if (ees0tmp.gt.0) then
+ ees0mij=dsqrt(ees0tmp)
+ else
+ ees0mij=0
+ endif
+ ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
+ ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
+!c write (iout,*) "i",i," j",j,
+!c & " ees0m",ees0m(num_conti,i)," ees0p",ees0p(num_conti,i)
+ ees0pij1=fac3/ees0pij
+ ees0mij1=fac3/ees0mij
+ fac3p=-3.0D0*fac3*rrij
+ ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
+ ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
+ ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
+ ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
+ ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
+ ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
+ ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
+ ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
+ ecosap=ecosa1+ecosa2
+ ecosbp=ecosb1+ecosb2
+ ecosgp=ecosg1+ecosg2
+ ecosam=ecosa1-ecosa2
+ ecosbm=ecosb1-ecosb2
+ ecosgm=ecosg1-ecosg2
+!C End diagnostics
+ facont_hb(num_conti,i)=fcont
+ fprimcont=fprimcont/rij
+ do k=1,3
+ gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
+ gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
+ enddo
+ gggp(1)=gggp(1)+ees0pijp*xj
+ gggp(2)=gggp(2)+ees0pijp*yj
+ gggp(3)=gggp(3)+ees0pijp*zj
+ gggm(1)=gggm(1)+ees0mijp*xj
+ gggm(2)=gggm(2)+ees0mijp*yj
+ gggm(3)=gggm(3)+ees0mijp*zj
+!C Derivatives due to the contact function
+ gacont_hbr(1,num_conti,i)=fprimcont*xj
+ gacont_hbr(2,num_conti,i)=fprimcont*yj
+ gacont_hbr(3,num_conti,i)=fprimcont*zj
+ do k=1,3
+!c
+!c Gradient of the correlation terms
+!c
+ gacontp_hb1(k,num_conti,i)= &
+ (ecosap*(dc_norm(k,j+nres)-cosa*dc_norm(k,i+nres)) &
+ + ecosbp*(erij(k)-cosb*dc_norm(k,i+nres)))*vbld_inv(i+nres)
+ gacontp_hb2(k,num_conti,i)= &
+ (ecosap*(dc_norm(k,i+nres)-cosa*dc_norm(k,j+nres)) &
+ + ecosgp*(erij(k)-cosg*dc_norm(k,j+nres)))*vbld_inv(j+nres)
+ gacontp_hb3(k,num_conti,i)=gggp(k)
+ gacontm_hb1(k,num_conti,i)= &
+ (ecosam*(dc_norm(k,j+nres)-cosa*dc_norm(k,i+nres)) &
+ + ecosbm*(erij(k)-cosb*dc_norm(k,i+nres)))*vbld_inv(i+nres)
+ gacontm_hb2(k,num_conti,i)= &
+ (ecosam*(dc_norm(k,i+nres)-cosa*dc_norm(k,j+nres))&
+ + ecosgm*(erij(k)-cosg*dc_norm(k,j+nres)))*vbld_inv(j+nres)
+ gacontm_hb3(k,num_conti,i)=gggm(k)
+ enddo
+ endif
+ endif
+ ENDIF
+ return
+ end subroutine eelsbij
+!------------------------------------------------------------------
+ subroutine sc_grad_nucl
+ use comm_locel
+ use calc_data_nucl
+ real(kind=8),dimension(3) :: dcosom1,dcosom2
+ eom1=eps2der*eps2rt_om1+sigder*sigsq_om1
+ eom2=eps2der*eps2rt_om2+sigder*sigsq_om2
+ eom12=evdwij*eps1_om12+eps2der*eps2rt_om12+sigder*sigsq_om12
+ do k=1,3
+ dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
+ dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
+ enddo
+ do k=1,3
+ gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
+ enddo
+ do k=1,3
+ gvdwsbx(k,i)=gvdwsbx(k,i)-gg(k) &
+ +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))&
+ +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+ gvdwsbx(k,j)=gvdwsbx(k,j)+gg(k) &
+ +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
+ +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ enddo
+!C
+!C Calculate the components of the gradient in DC and X
+!C
+ do l=1,3
+ gvdwsbc(l,i)=gvdwsbc(l,i)-gg(l)
+ gvdwsbc(l,j)=gvdwsbc(l,j)+gg(l)
+ enddo
+ return
+ end subroutine sc_grad_nucl
+!-----------------------------------------------------------------------
+ subroutine esb(esbloc)
+!C Calculate the local energy of a side chain and its derivatives in the
+!C corresponding virtual-bond valence angles THETA and the spherical angles
+!C ALPHA and OMEGA derived from AM1 all-atom calculations.
+!C added by Urszula Kozlowska. 07/11/2007
+!C
+ real(kind=8),dimension(3):: x_prime,y_prime,z_prime
+ real(kind=8),dimension(9):: x
+ real(kind=8) :: sumene,dsc_i,dp2_i,xx,yy,zz,sumene1, &
+ sumene2,sumene3,sumene4,s1,s1_6,s2,s2_6,&
+ de_dxx,de_dyy,de_dzz,de_dt,s1_t,s1_6_t,s2_t,s2_6_t
+ real(kind=8),dimension(3):: dXX_Ci1,dYY_Ci1,dZZ_Ci1,dXX_Ci,&
+ dYY_Ci,dZZ_Ci,dXX_XYZ,dYY_XYZ,dZZ_XYZ,dt_dCi,dt_dCi1
+ real(kind=8) :: esbloc,delta,cosfac2,cosfac,sinfac2,sinfac,de_dtt,&
+ cossc,cossc1,cosfac2xx,sinfac2yy,pom1,pom
+ integer::it,nlobit,i,j,k
+! common /sccalc/ time11,time12,time112,theti,it,nlobit
+ delta=0.02d0*pi
+ esbloc=0.0D0
+ do i=loc_start_nucl,loc_end_nucl
+ if (itype(i,2).eq.ntyp1_molec(2)) cycle
+ costtab(i+1) =dcos(theta(i+1))
+ sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
+ cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
+ sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
+ cosfac2=0.5d0/(1.0d0+costtab(i+1))
+ cosfac=dsqrt(cosfac2)
+ sinfac2=0.5d0/(1.0d0-costtab(i+1))
+ sinfac=dsqrt(sinfac2)
+ it=itype(i,2)
+ if (it.eq.10) goto 1
+
+!c
+!C Compute the axes of tghe local cartesian coordinates system; store in
+!c x_prime, y_prime and z_prime
+!c
+ do j=1,3
+ x_prime(j) = 0.00
+ y_prime(j) = 0.00
+ z_prime(j) = 0.00
+ enddo
+!C write(2,*) "dc_norm", dc_norm(1,i+nres),dc_norm(2,i+nres),
+!C & dc_norm(3,i+nres)
+ do j = 1,3
+ x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
+ y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
+ enddo
+ do j = 1,3
+ z_prime(j) = -uz(j,i-1)
+! z_prime(j)=0.0
+ enddo
+
+ xx=0.0d0
+ yy=0.0d0
+ zz=0.0d0
+ do j = 1,3
+ xx = xx + x_prime(j)*dc_norm(j,i+nres)
+ yy = yy + y_prime(j)*dc_norm(j,i+nres)
+ zz = zz + z_prime(j)*dc_norm(j,i+nres)
+ enddo
+
+ xxtab(i)=xx
+ yytab(i)=yy
+ zztab(i)=zz
+ it=itype(i,2)
+ do j = 1,9
+ x(j) = sc_parmin_nucl(j,it)
+ enddo
+#ifdef CHECK_COORD
+!Cc diagnostics - remove later
+ xx1 = dcos(alph(2))
+ yy1 = dsin(alph(2))*dcos(omeg(2))
+ zz1 = -dsin(alph(2))*dsin(omeg(2))
+ write(2,'(3f8.1,3f9.3,1x,3f9.3)') &
+ alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,&
+ xx1,yy1,zz1
+!C," --- ", xx_w,yy_w,zz_w
+!c end diagnostics
+#endif
+ sumene = enesc_nucl(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ esbloc = esbloc + sumene
+ sumene2= enesc_nucl(x,xx,yy,0.0d0,cost2tab(i+1),sint2tab(i+1))
+! print *,"enecomp",sumene,sumene2
+! if (energy_dec) write(iout,*) "i",i," esbloc",sumene,esbloc,xx,yy,zz
+! if (energy_dec) write(iout,*) "x",(x(k),k=1,9)
+#ifdef DEBUG
+ write (2,*) "x",(x(k),k=1,9)
+!C
+!C This section to check the numerical derivatives of the energy of ith side
+!C chain in xx, yy, zz, and theta. Use the -DDEBUG compiler option or insert
+!C #define DEBUG in the code to turn it on.
+!C
+ write (2,*) "sumene =",sumene
+ aincr=1.0d-7
+ xxsave=xx
+ xx=xx+aincr
+ write (2,*) xx,yy,zz
+ sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ de_dxx_num=(sumenep-sumene)/aincr
+ xx=xxsave
+ write (2,*) "xx+ sumene from enesc=",sumenep,sumene
+ yysave=yy
+ yy=yy+aincr
+ write (2,*) xx,yy,zz
+ sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ de_dyy_num=(sumenep-sumene)/aincr
+ yy=yysave
+ write (2,*) "yy+ sumene from enesc=",sumenep,sumene
+ zzsave=zz
+ zz=zz+aincr
+ write (2,*) xx,yy,zz
+ sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ de_dzz_num=(sumenep-sumene)/aincr
+ zz=zzsave
+ write (2,*) "zz+ sumene from enesc=",sumenep,sumene
+ costsave=cost2tab(i+1)
+ sintsave=sint2tab(i+1)
+ cost2tab(i+1)=dcos(0.5d0*(theta(i+1)+aincr))
+ sint2tab(i+1)=dsin(0.5d0*(theta(i+1)+aincr))
+ sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ de_dt_num=(sumenep-sumene)/aincr
+ write (2,*) " t+ sumene from enesc=",sumenep,sumene
+ cost2tab(i+1)=costsave
+ sint2tab(i+1)=sintsave
+!C End of diagnostics section.
+#endif
+!C
+!C Compute the gradient of esc
+!C
+ de_dxx=x(1)+2*x(4)*xx+x(7)*zz+x(8)*yy
+ de_dyy=x(2)+2*x(5)*yy+x(8)*xx+x(9)*zz
+ de_dzz=x(3)+2*x(6)*zz+x(7)*xx+x(9)*yy
+ de_dtt=0.0d0
+#ifdef DEBUG
+ write (2,*) "x",(x(k),k=1,9)
+ write (2,*) "xx",xx," yy",yy," zz",zz
+ write (2,*) "de_xx ",de_xx," de_yy ",de_yy,&
+ " de_zz ",de_zz," de_tt ",de_tt
+ write (2,*) "de_xx_num",de_dxx_num," de_yy_num",de_dyy_num,&
+ " de_zz_num",de_dzz_num," de_dt_num",de_dt_num
+#endif
+!C
+ cossc=scalar(dc_norm(1,i),dc_norm(1,i+nres))
+ cossc1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres))
+ cosfac2xx=cosfac2*xx
+ sinfac2yy=sinfac2*yy
+ do k = 1,3
+ dt_dCi(k) = -(dc_norm(k,i-1)+costtab(i+1)*dc_norm(k,i))*&
+ vbld_inv(i+1)
+ dt_dCi1(k)= -(dc_norm(k,i)+costtab(i+1)*dc_norm(k,i-1))*&
+ vbld_inv(i)
+ pom=(dC_norm(k,i+nres)-cossc*dC_norm(k,i))*vbld_inv(i+1)
+ pom1=(dC_norm(k,i+nres)-cossc1*dC_norm(k,i-1))*vbld_inv(i)
+!c write (iout,*) "i",i," k",k," pom",pom," pom1",pom1,
+!c & " dt_dCi",dt_dCi(k)," dt_dCi1",dt_dCi1(k)
+!c write (iout,*) "dC_norm",(dC_norm(j,i),j=1,3),
+!c & (dC_norm(j,i-1),j=1,3)," vbld_inv",vbld_inv(i+1),vbld_inv(i)
+ dXX_Ci(k)=pom*cosfac-dt_dCi(k)*cosfac2xx
+ dXX_Ci1(k)=-pom1*cosfac-dt_dCi1(k)*cosfac2xx
+ dYY_Ci(k)=pom*sinfac+dt_dCi(k)*sinfac2yy
+ dYY_Ci1(k)=pom1*sinfac+dt_dCi1(k)*sinfac2yy
+ dZZ_Ci1(k)=0.0d0
+ dZZ_Ci(k)=0.0d0
+ do j=1,3
+ dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres)
+ dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres)
+ enddo
+
+ dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
+ dYY_XYZ(k)=vbld_inv(i+nres)*(y_prime(k)-yy*dC_norm(k,i+nres))
+ dZZ_XYZ(k)=vbld_inv(i+nres)*(z_prime(k)-zz*dC_norm(k,i+nres))
+!c
+ dt_dCi(k) = -dt_dCi(k)/sinttab(i+1)
+ dt_dCi1(k)= -dt_dCi1(k)/sinttab(i+1)
+ enddo
+
+ do k=1,3
+ dXX_Ctab(k,i)=dXX_Ci(k)
+ dXX_C1tab(k,i)=dXX_Ci1(k)
+ dYY_Ctab(k,i)=dYY_Ci(k)
+ dYY_C1tab(k,i)=dYY_Ci1(k)
+ dZZ_Ctab(k,i)=dZZ_Ci(k)
+ dZZ_C1tab(k,i)=dZZ_Ci1(k)
+ dXX_XYZtab(k,i)=dXX_XYZ(k)
+ dYY_XYZtab(k,i)=dYY_XYZ(k)
+ dZZ_XYZtab(k,i)=dZZ_XYZ(k)
+ enddo
+ do k = 1,3
+!c write (iout,*) "k",k," dxx_ci1",dxx_ci1(k)," dyy_ci1",
+!c & dyy_ci1(k)," dzz_ci1",dzz_ci1(k)
+!c write (iout,*) "k",k," dxx_ci",dxx_ci(k)," dyy_ci",
+!c & dyy_ci(k)," dzz_ci",dzz_ci(k)
+!c write (iout,*) "k",k," dt_dci",dt_dci(k)," dt_dci",
+!c & dt_dci(k)
+!c write (iout,*) "k",k," dxx_XYZ",dxx_XYZ(k)," dyy_XYZ",
+!c & dyy_XYZ(k)," dzz_XYZ",dzz_XYZ(k)
+ gsbloc(k,i-1)=gsbloc(k,i-1)+(de_dxx*dxx_ci1(k) &
+ +de_dyy*dyy_ci1(k)+de_dzz*dzz_ci1(k)+de_dt*dt_dCi1(k))
+ gsbloc(k,i)=gsbloc(k,i)+(de_dxx*dxx_Ci(k) &
+ +de_dyy*dyy_Ci(k)+de_dzz*dzz_Ci(k)+de_dt*dt_dCi(k))
+ gsblocx(k,i)= de_dxx*dxx_XYZ(k)&
+ +de_dyy*dyy_XYZ(k)+de_dzz*dzz_XYZ(k)
+! print *,i,de_dxx*dxx_ci1(k)+de_dyy*dyy_ci1(k),de_dzz*dzz_ci1(k)*2
+ enddo
+!c write(iout,*) "ENERGY GRAD = ", (gsbloc(k,i-1),k=1,3),
+!c & (gsbloc(k,i),k=1,3),(gsblocx(k,i),k=1,3)
+
+!C to check gradient call subroutine check_grad
+
+ 1 continue
+ enddo
+ return
+ end subroutine esb
+!=-------------------------------------------------------
+ real(kind=8) function enesc_nucl(x,xx,yy,zz,cost2,sint2)
+! implicit none
+ real(kind=8),dimension(9):: x(9)
+ real(kind=8) :: xx,yy,zz,cost2,sint2,sumene1,sumene2, &
+ sumene3,sumene4,sumene,dsc_i,dp2_i,dscp1,dscp2,s1,s1_6,s2,s2_6
+ integer i
+!c write (2,*) "enesc"
+!c write (2,*) "x",(x(i),i=1,9)
+!c write(2,*)"xx",xx," yy",yy," zz",zz," cost2",cost2," sint2",sint2
+ sumene=x(1)*xx+x(2)*yy+x(3)*zz+x(4)*xx**2 &
+ + x(5)*yy**2+x(6)*zz**2+x(7)*xx*zz+x(8)*xx*yy &
+ + x(9)*yy*zz
+ enesc_nucl=sumene
+ return
+ end function enesc_nucl
+!-----------------------------------------------------------------------------
+ subroutine multibody_hb_nucl(ecorr,ecorr3,n_corr,n_corr1)
+#ifdef MPI
+ include 'mpif.h'
+ integer,parameter :: max_cont=2000
+ integer,parameter:: max_dim=2*(8*3+6)
+ integer, parameter :: msglen1=max_cont*max_dim
+ integer,parameter :: msglen2=2*msglen1
+ integer source,CorrelType,CorrelID,Error
+ real(kind=8) :: buffer(max_cont,max_dim)
+ integer status(MPI_STATUS_SIZE)
+ integer :: ierror,nbytes
+#endif
+ real(kind=8),dimension(3):: gx(3),gx1(3)
+ real(kind=8) :: time00
+ logical lprn,ldone
+ integer i,j,i1,j1,jj,kk,num_conti,num_conti1,nn
+ real(kind=8) ecorr,ecorr3
+ integer :: n_corr,n_corr1,mm,msglen
+!C Set lprn=.true. for debugging
+ lprn=.false.
+ n_corr=0
+ n_corr1=0
+#ifdef MPI
+ if(.not.allocated(zapas2)) allocate(zapas2(3,maxconts,nres,8))
+
+ if (nfgtasks.le.1) goto 30
+ if (lprn) then
+ write (iout,'(a)') 'Contact function values:'
+ do i=nnt,nct-1
+ write (iout,'(2i3,50(1x,i2,f5.2))') &
+ i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), &
+ j=1,num_cont_hb(i))
+ enddo
+ endif
+!C Caution! Following code assumes that electrostatic interactions concerning
+!C a given atom are split among at most two processors!
+ CorrelType=477
+ CorrelID=fg_rank+1
+ ldone=.false.
+ do i=1,max_cont
+ do j=1,max_dim
+ buffer(i,j)=0.0D0
+ enddo
+ enddo
+ mm=mod(fg_rank,2)
+!c write (*,*) 'MyRank',MyRank,' mm',mm
+ if (mm) 20,20,10
+ 10 continue
+!c write (*,*) 'Sending: MyRank',MyRank,' mm',mm,' ldone',ldone
+ if (fg_rank.gt.0) then
+!C Send correlation contributions to the preceding processor
+ msglen=msglen1
+ nn=num_cont_hb(iatel_s_nucl)
+ call pack_buffer(max_cont,max_dim,iatel_s,0,buffer)
+!c write (*,*) 'The BUFFER array:'
+!c do i=1,nn
+!c write (*,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,30)
+!c enddo
+ if (ielstart_nucl(iatel_s_nucl).gt.iatel_s_nucl+ispp) then
+ msglen=msglen2
+ call pack_buffer(max_cont,max_dim,iatel_s+1,30,buffer)
+!C Clear the contacts of the atom passed to the neighboring processor
+ nn=num_cont_hb(iatel_s_nucl+1)
+!c do i=1,nn
+!c write (*,'(i2,9(3f8.3,2x))') i,(buffer(i,j+30),j=1,30)
+!c enddo
+ num_cont_hb(iatel_s_nucl)=0
+ endif
+!cd write (iout,*) 'Processor ',fg_rank,MyRank,
+!cd & ' is sending correlation contribution to processor',fg_rank-1,
+!cd & ' msglen=',msglen
+!c write (*,*) 'Processor ',fg_rank,MyRank,
+!c & ' is sending correlation contribution to processor',fg_rank-1,
+!c & ' msglen=',msglen,' CorrelType=',CorrelType
+ time00=MPI_Wtime()
+ call MPI_Send(buffer,msglen,MPI_DOUBLE_PRECISION,fg_rank-1, &
+ CorrelType,FG_COMM,IERROR)
+ time_sendrecv=time_sendrecv+MPI_Wtime()-time00
+!cd write (iout,*) 'Processor ',fg_rank,
+!cd & ' has sent correlation contribution to processor',fg_rank-1,
+!cd & ' msglen=',msglen,' CorrelID=',CorrelID
+!c write (*,*) 'Processor ',fg_rank,
+!c & ' has sent correlation contribution to processor',fg_rank-1,
+!c & ' msglen=',msglen,' CorrelID=',CorrelID
+!c msglen=msglen1
+ endif ! (fg_rank.gt.0)
+ if (ldone) goto 30
+ ldone=.true.
+ 20 continue
+!c write (*,*) 'Receiving: MyRank',MyRank,' mm',mm,' ldone',ldone
+ if (fg_rank.lt.nfgtasks-1) then
+!C Receive correlation contributions from the next processor
+ msglen=msglen1
+ if (ielend_nucl(iatel_e_nucl).lt.nct_molec(2)-1) msglen=msglen2
+!cd write (iout,*) 'Processor',fg_rank,
+!cd & ' is receiving correlation contribution from processor',fg_rank+1,
+!cd & ' msglen=',msglen,' CorrelType=',CorrelType
+!c write (*,*) 'Processor',fg_rank,
+!c &' is receiving correlation contribution from processor',fg_rank+1,
+!c & ' msglen=',msglen,' CorrelType=',CorrelType
+ time00=MPI_Wtime()
+ nbytes=-1
+ do while (nbytes.le.0)
+ call MPI_Probe(fg_rank+1,CorrelType,FG_COMM,status,IERROR)
+ call MPI_Get_count(status,MPI_DOUBLE_PRECISION,nbytes,IERROR)
+ enddo
+!c print *,'Processor',myrank,' msglen',msglen,' nbytes',nbytes
+ call MPI_Recv(buffer,nbytes,MPI_DOUBLE_PRECISION, &
+ fg_rank+1,CorrelType,FG_COMM,status,IERROR)
+ time_sendrecv=time_sendrecv+MPI_Wtime()-time00
+!c write (*,*) 'Processor',fg_rank,
+!c &' has received correlation contribution from processor',fg_rank+1,
+!c & ' msglen=',msglen,' nbytes=',nbytes
+!c write (*,*) 'The received BUFFER array:'
+!c do i=1,max_cont
+!c write (*,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,60)
+!c enddo
+ if (msglen.eq.msglen1) then
+ call unpack_buffer(max_cont,max_dim,iatel_e_nucl+1,0,buffer)
+ else if (msglen.eq.msglen2) then
+ call unpack_buffer(max_cont,max_dim,iatel_e_nucl,0,buffer)
+ call unpack_buffer(max_cont,max_dim,iatel_e_nucl+1,30,buffer)
+ else
+ write (iout,*) &
+ 'ERROR!!!! message length changed while processing correlations.'
+ write (*,*) &
+ 'ERROR!!!! message length changed while processing correlations.'
+ call MPI_Abort(MPI_COMM_WORLD,Error,IERROR)
+ endif ! msglen.eq.msglen1
+ endif ! fg_rank.lt.nfgtasks-1
+ if (ldone) goto 30
+ ldone=.true.
+ goto 10
+ 30 continue
+#endif
+ if (lprn) then
+ write (iout,'(a)') 'Contact function values:'
+ do i=nnt_molec(2),nct_molec(2)-1
+ write (iout,'(2i3,50(1x,i2,f5.2))') &
+ i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), &
+ j=1,num_cont_hb(i))
+ enddo
+ endif
+ ecorr=0.0D0
+ ecorr3=0.0d0
+!C Remove the loop below after debugging !!!
+! do i=nnt_molec(2),nct_molec(2)
+! do j=1,3
+! gradcorr_nucl(j,i)=0.0D0
+! gradxorr_nucl(j,i)=0.0D0
+! gradcorr3_nucl(j,i)=0.0D0
+! gradxorr3_nucl(j,i)=0.0D0
+! enddo
+! enddo
+! print *,"iatsc_s_nucl,iatsc_e_nucl",iatsc_s_nucl,iatsc_e_nucl
+!C Calculate the local-electrostatic correlation terms
+ do i=iatsc_s_nucl,iatsc_e_nucl
+ i1=i+1
+ num_conti=num_cont_hb(i)
+ num_conti1=num_cont_hb(i+1)
+! print *,i,num_conti,num_conti1
+ do jj=1,num_conti
+ j=jcont_hb(jj,i)
+ do kk=1,num_conti1
+ j1=jcont_hb(kk,i1)
+!c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
+!c & ' jj=',jj,' kk=',kk
+ if (j1.eq.j+1 .or. j1.eq.j-1) then
+!C
+!C Contacts I-J and (I+1)-(J+1) or (I+1)-(J-1) occur simultaneously.
+!C The system gains extra energy.
+!C Tentative expression & coefficients; assumed d(stacking)=4.5 A,
+!C parallel dipoles of stacknig bases and sin(mui)sin(muj)/eps/d^3=0.7
+!C Need to implement full formulas 34 and 35 from Liwo et al., 1998.
+!C
+ ecorr=ecorr+ehbcorr_nucl(i,j,i+1,j1,jj,kk,0.528D0,0.132D0)
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
+ 'ecorrh',i,j,ehbcorr_nucl(i,j,i+1,j1,jj,kk,0.528D0,0.132D0)
+ n_corr=n_corr+1
+ else if (j1.eq.j) then
+!C
+!C Contacts I-J and I-(J+1) occur simultaneously.
+!C The system loses extra energy.
+!C Tentative expression & c?oefficients; assumed d(stacking)=4.5 A,
+!C parallel dipoles of stacknig bases and sin(mui)sin(muj)/eps/d^3=0.7
+!C Need to implement full formulas 32 from Liwo et al., 1998.
+!C
+!c write (iout,*) 'ecorr3: i=',i,' j=',j,' i1=',i1,' j1=',j1,
+!c & ' jj=',jj,' kk=',kk
+ ecorr3=ecorr3+ehbcorr3_nucl(i,j,i+1,j,jj,kk,0.310D0,-0.155D0)
+ endif
+ enddo ! kk
+ do kk=1,num_conti
+ j1=jcont_hb(kk,i)
+!c write (iout,*) 'ecorr3: i=',i,' j=',j,' i1=',i1,' j1=',j1,
+!c & ' jj=',jj,' kk=',kk
+ if (j1.eq.j+1) then
+!C Contacts I-J and (I+1)-J occur simultaneously.
+!C The system loses extra energy.
+ ecorr3=ecorr3+ehbcorr3_nucl(i,j,i,j+1,jj,kk,0.310D0,-0.155D0)
+ endif ! j1==j+1
+ enddo ! kk
+ enddo ! jj
+ enddo ! i
+ return
+ end subroutine multibody_hb_nucl
+!-----------------------------------------------------------
+ real(kind=8) function ehbcorr_nucl(i,j,k,l,jj,kk,coeffp,coeffm)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.DERIV'
+! include 'COMMON.INTERACT'
+! include 'COMMON.CONTACTS'
+ real(kind=8),dimension(3) :: gx,gx1
+ logical :: lprn
+!el local variables
+ integer :: i,j,k,l,jj,kk,ll,ilist,m, iresshield
+ real(kind=8) :: coeffp,coeffm,eij,ekl,ees0pij,ees0pkl,ees0mij,&
+ ees0mkl,ees,coeffpees0pij,coeffmees0mij,&
+ coeffpees0pkl,coeffmees0mkl,gradlongij,gradlongkl, &
+ rlocshield
+
+ lprn=.false.
+ eij=facont_hb(jj,i)
+ ekl=facont_hb(kk,k)
+ ees0pij=ees0p(jj,i)
+ ees0pkl=ees0p(kk,k)
+ ees0mij=ees0m(jj,i)
+ ees0mkl=ees0m(kk,k)
+ ekont=eij*ekl
+ ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl)
+! print *,"ehbcorr_nucl",ekont,ees
+!cd ees=-(coeffp*ees0pkl+coeffm*ees0mkl)
+!C Following 4 lines for diagnostics.
+!cd ees0pkl=0.0D0
+!cd ees0pij=1.0D0
+!cd ees0mkl=0.0D0
+!cd ees0mij=1.0D0
+!cd write (iout,*)'Contacts have occurred for nucleic bases',
+!cd & i,j,' fcont:',eij,' eij',' eesij',ees0pij,ees0mij,' and ',k,l
+!cd & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' ees=',ees
+!C Calculate the multi-body contribution to energy.
+! ecorr_nucl=ecorr_nucl+ekont*ees
+!C Calculate multi-body contributions to the gradient.
+ coeffpees0pij=coeffp*ees0pij
+ coeffmees0mij=coeffm*ees0mij
+ coeffpees0pkl=coeffp*ees0pkl
+ coeffmees0mkl=coeffm*ees0mkl
+ do ll=1,3
+ gradxorr_nucl(ll,i)=gradxorr_nucl(ll,i) &
+ -ekont*(coeffpees0pkl*gacontp_hb1(ll,jj,i)+&
+ coeffmees0mkl*gacontm_hb1(ll,jj,i))
+ gradxorr_nucl(ll,j)=gradxorr_nucl(ll,j) &
+ -ekont*(coeffpees0pkl*gacontp_hb2(ll,jj,i)+&
+ coeffmees0mkl*gacontm_hb2(ll,jj,i))
+ gradxorr_nucl(ll,k)=gradxorr_nucl(ll,k) &
+ -ekont*(coeffpees0pij*gacontp_hb1(ll,kk,k)+&
+ coeffmees0mij*gacontm_hb1(ll,kk,k))
+ gradxorr_nucl(ll,l)=gradxorr_nucl(ll,l) &
+ -ekont*(coeffpees0pij*gacontp_hb2(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb2(ll,kk,k))
+ gradlongij=ees*ekl*gacont_hbr(ll,jj,i)- &
+ ekont*(coeffpees0pkl*gacontp_hb3(ll,jj,i)+ &
+ coeffmees0mkl*gacontm_hb3(ll,jj,i))
+ gradcorr_nucl(ll,j)=gradcorr_nucl(ll,j)+gradlongij
+ gradcorr_nucl(ll,i)=gradcorr_nucl(ll,i)-gradlongij
+ gradlongkl=ees*eij*gacont_hbr(ll,kk,k)- &
+ ekont*(coeffpees0pij*gacontp_hb3(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb3(ll,kk,k))
+ gradcorr_nucl(ll,l)=gradcorr_nucl(ll,l)+gradlongkl
+ gradcorr_nucl(ll,k)=gradcorr_nucl(ll,k)-gradlongkl
+ gradxorr_nucl(ll,i)=gradxorr_nucl(ll,i)-gradlongij
+ gradxorr_nucl(ll,j)=gradxorr_nucl(ll,j)+gradlongij
+ gradxorr_nucl(ll,k)=gradxorr_nucl(ll,k)-gradlongkl
+ gradxorr_nucl(ll,l)=gradxorr_nucl(ll,l)+gradlongkl
+ enddo
+ ehbcorr_nucl=ekont*ees
+ return
+ end function ehbcorr_nucl
+!-------------------------------------------------------------------------
+
+ real(kind=8) function ehbcorr3_nucl(i,j,k,l,jj,kk,coeffp,coeffm)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.DERIV'
+! include 'COMMON.INTERACT'
+! include 'COMMON.CONTACTS'
+ real(kind=8),dimension(3) :: gx,gx1
+ logical :: lprn
+!el local variables
+ integer :: i,j,k,l,jj,kk,ll,ilist,m, iresshield
+ real(kind=8) :: coeffp,coeffm,eij,ekl,ees0pij,ees0pkl,ees0mij,&
+ ees0mkl,ees,coeffpees0pij,coeffmees0mij,&
+ coeffpees0pkl,coeffmees0mkl,gradlongij,gradlongkl, &
+ rlocshield
+
+ lprn=.false.
+ eij=facont_hb(jj,i)
+ ekl=facont_hb(kk,k)
+ ees0pij=ees0p(jj,i)
+ ees0pkl=ees0p(kk,k)
+ ees0mij=ees0m(jj,i)
+ ees0mkl=ees0m(kk,k)
+ ekont=eij*ekl
+ ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl)
+!cd ees=-(coeffp*ees0pkl+coeffm*ees0mkl)
+!C Following 4 lines for diagnostics.
+!cd ees0pkl=0.0D0
+!cd ees0pij=1.0D0
+!cd ees0mkl=0.0D0
+!cd ees0mij=1.0D0
+!cd write (iout,*)'Contacts have occurred for nucleic bases',
+!cd & i,j,' fcont:',eij,' eij',' eesij',ees0pij,ees0mij,' and ',k,l
+!cd & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' ees=',ees
+!C Calculate the multi-body contribution to energy.
+! ecorr=ecorr+ekont*ees
+!C Calculate multi-body contributions to the gradient.
+ coeffpees0pij=coeffp*ees0pij
+ coeffmees0mij=coeffm*ees0mij
+ coeffpees0pkl=coeffp*ees0pkl
+ coeffmees0mkl=coeffm*ees0mkl
+ do ll=1,3
+ gradxorr3_nucl(ll,i)=gradxorr3_nucl(ll,i) &
+ -ekont*(coeffpees0pkl*gacontp_hb1(ll,jj,i)+&
+ coeffmees0mkl*gacontm_hb1(ll,jj,i))
+ gradxorr3_nucl(ll,j)=gradxorr3_nucl(ll,j) &
+ -ekont*(coeffpees0pkl*gacontp_hb2(ll,jj,i)+ &
+ coeffmees0mkl*gacontm_hb2(ll,jj,i))
+ gradxorr3_nucl(ll,k)=gradxorr3_nucl(ll,k) &
+ -ekont*(coeffpees0pij*gacontp_hb1(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb1(ll,kk,k))
+ gradxorr3_nucl(ll,l)=gradxorr3_nucl(ll,l) &
+ -ekont*(coeffpees0pij*gacontp_hb2(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb2(ll,kk,k))
+ gradlongij=ees*ekl*gacont_hbr(ll,jj,i)- &
+ ekont*(coeffpees0pkl*gacontp_hb3(ll,jj,i)+ &
+ coeffmees0mkl*gacontm_hb3(ll,jj,i))
+ gradcorr3_nucl(ll,j)=gradcorr3_nucl(ll,j)+gradlongij
+ gradcorr3_nucl(ll,i)=gradcorr3_nucl(ll,i)-gradlongij
+ gradlongkl=ees*eij*gacont_hbr(ll,kk,k)- &
+ ekont*(coeffpees0pij*gacontp_hb3(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb3(ll,kk,k))
+ gradcorr3_nucl(ll,l)=gradcorr3_nucl(ll,l)+gradlongkl
+ gradcorr3_nucl(ll,k)=gradcorr3_nucl(ll,k)-gradlongkl
+ gradxorr3_nucl(ll,i)=gradxorr3_nucl(ll,i)-gradlongij
+ gradxorr3_nucl(ll,j)=gradxorr3_nucl(ll,j)+gradlongij
+ gradxorr3_nucl(ll,k)=gradxorr3_nucl(ll,k)-gradlongkl
+ gradxorr3_nucl(ll,l)=gradxorr3_nucl(ll,l)+gradlongkl
+ enddo
+ ehbcorr3_nucl=ekont*ees
+ return
+ end function ehbcorr3_nucl
+#ifdef MPI
+ subroutine pack_buffer(dimen1,dimen2,atom,indx,buffer)
+ integer dimen1,dimen2,atom,indx,numcont,i,ii,k,j,num_kont,num_kont_old
+ real(kind=8):: buffer(dimen1,dimen2)
+ num_kont=num_cont_hb(atom)
+ do i=1,num_kont
+ do k=1,8
+ do j=1,3
+ buffer(i,indx+(k-1)*3+j)=zapas2(j,i,atom,k)
+ enddo ! j
+ enddo ! k
+ buffer(i,indx+25)=facont_hb(i,atom)
+ buffer(i,indx+26)=ees0p(i,atom)
+ buffer(i,indx+27)=ees0m(i,atom)
+ buffer(i,indx+28)=d_cont(i,atom)
+ buffer(i,indx+29)=dfloat(jcont_hb(i,atom))
+ enddo ! i
+ buffer(1,indx+30)=dfloat(num_kont)
+ return
+ end subroutine pack_buffer
+!c------------------------------------------------------------------------------
+ subroutine unpack_buffer(dimen1,dimen2,atom,indx,buffer)
+ integer dimen1,dimen2,atom,indx,numcont,i,ii,k,j,num_kont,num_kont_old
+ real(kind=8):: buffer(dimen1,dimen2)
+! double precision zapas
+! common /contacts_hb/ zapas(3,maxconts,maxres,8),
+! & facont_hb(maxconts,maxres),ees0p(maxconts,maxres),
+! & ees0m(maxconts,maxres),d_cont(maxconts,maxres),
+! & num_cont_hb(maxres),jcont_hb(maxconts,maxres)
+ num_kont=buffer(1,indx+30)
+ num_kont_old=num_cont_hb(atom)
+ num_cont_hb(atom)=num_kont+num_kont_old
+ do i=1,num_kont
+ ii=i+num_kont_old
+ do k=1,8
+ do j=1,3
+ zapas2(j,ii,atom,k)=buffer(i,indx+(k-1)*3+j)
+ enddo ! j
+ enddo ! k
+ facont_hb(ii,atom)=buffer(i,indx+25)
+ ees0p(ii,atom)=buffer(i,indx+26)
+ ees0m(ii,atom)=buffer(i,indx+27)
+ d_cont(i,atom)=buffer(i,indx+28)
+ jcont_hb(ii,atom)=buffer(i,indx+29)
+ enddo ! i
+ return
+ end subroutine unpack_buffer
+!c------------------------------------------------------------------------------
+#endif
+ subroutine ecatcat(ecationcation)
+ integer :: i,j,itmp,xshift,yshift,zshift,subchap,k
+ real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,&
+ r7,r4,ecationcation,k0,rcal
+ real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
+ dist_init,dist_temp,Evan1cat,Evan2cat,Eeleccat
+ real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,&
+ gg,r
+
+ ecationcation=0.0d0
+ if (nres_molec(5).eq.0) return
+ rcat0=3.472
+ epscalc=0.05
+ r06 = rcat0**6
+ r012 = r06**2
+ k0 = 332.0*(2.0*2.0)/80.0
+ itmp=0
+
+ do i=1,4
+ itmp=itmp+nres_molec(i)
+ enddo
+! write(iout,*) "itmp",itmp
+ do i=itmp+1,itmp+nres_molec(5)-1
+
+ xi=c(1,i)
+ yi=c(2,i)
+ zi=c(3,i)
+
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+
+ do j=i+1,itmp+nres_molec(5)
+! print *,i,j,'catcat'
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+! write(iout,*) c(1,i),xi,xj,"xy",boxxsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+ rcal =xj**2+yj**2+zj**2
+ ract=sqrt(rcal)
+! rcat0=3.472
+! epscalc=0.05
+! r06 = rcat0**6
+! r012 = r06**2
+! k0 = 332*(2*2)/80
+ Evan1cat=epscalc*(r012/rcal**6)
+ Evan2cat=epscalc*2*(r06/rcal**3)
+ Eeleccat=k0/ract
+ r7 = rcal**7
+ r4 = rcal**4
+ r(1)=xj
+ r(2)=yj
+ r(3)=zj
+ do k=1,3
+ dEvan1Cmcat(k)=-12*r(k)*epscalc*r012/r7
+ dEvan2Cmcat(k)=-12*r(k)*epscalc*r06/r4
+ dEeleccat(k)=-k0*r(k)/ract**3
+ enddo
+ do k=1,3
+ gg(k) = dEvan1Cmcat(k)+dEvan2Cmcat(k)+dEeleccat(k)
+ gradcatcat(k,i)=gradcatcat(k,i)-gg(k)
+ gradcatcat(k,j)=gradcatcat(k,j)+gg(k)
+ enddo
+
+! write(iout,*) "ecatcat",i,j, ecationcation,xj,yj,zj
+ ecationcation=ecationcation+Evan1cat+Evan2cat+Eeleccat
+ enddo
+ enddo
+ return
+ end subroutine ecatcat
+!---------------------------------------------------------------------------
+ subroutine ecat_prot(ecation_prot)
+ integer i,j,k,subchap,itmp,inum
+ real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,&
+ r7,r4,ecationcation
+ real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
+ dist_init,dist_temp,ecation_prot,rcal,rocal, &
+ Evan1,Evan2,EC,cm1mag,DASGL,delta,r0p,Epepcat, &
+ catl,cml,calpl, Etotal_p, Etotal_m,rtab,wdip,wmodquad,wquad1, &
+ wquad2,wvan1,E1,E2,wconst,wvan2,rcpm,dcmag,sin2thet,sinthet, &
+ costhet,v1m,v2m,wh2o,wc,rsecp,Ir,Irsecp,Irthrp,Irfourp,Irfiftp,&
+ Irsistp,Irseven,Irtwelv,Irthir,dE1dr,dE2dr,dEdcos,wquad2p,opt, &
+ rs,rthrp,rfourp,rsixp,reight,Irsixp,Ireight,Irtw,Irfourt, &
+ opt1,opt2,opt3,opt4,opt5,opt6,opt7,opt8,opt9,opt10,opt11,opt12,&
+ opt13,opt14,opt15,opt16,opt17,opt18,opt19, &
+ Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip
+ real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,&
+ gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, &
+ dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, &
+ tab1,tab2,tab3,diff,cm1,sc,p,tcat,talp,cm,drcp,drcp_norm,vcat, &
+ v1,v2,v3,myd_norm,dx,vcm,valpha,drdpep,dcosdpep,dcosddci,dEdpep,&
+ dEcCat,dEdipCm,dEdipCalp,dEquad1Cat,dEquad1Cm,dEquad1Calp, &
+ dEquad2Cat,dEquad2Cm,dEquad2Calpd,Evan1Cat,dEvan1Calp,dEvan2Cat,&
+ dEvan2Calp,dEtotalCat,dscvec,dEcCm,dEcCalp,dEdipCat,dEquad2Calp,&
+ dEvan1Cat
+ real(kind=8),dimension(6) :: vcatprm
+ ecation_prot=0.0d0
+! first lets calculate interaction with peptide groups
+ if (nres_molec(5).eq.0) return
+ wconst=78
+ wdip =1.092777950857032D2
+ wdip=wdip/wconst
+ wmodquad=-2.174122713004870D4
+ wmodquad=wmodquad/wconst
+ wquad1 = 3.901232068562804D1
+ wquad1=wquad1/wconst
+ wquad2 = 3
+ wquad2=wquad2/wconst
+ wvan1 = 0.1
+ wvan2 = 6
+ itmp=0
+ do i=1,4
+ itmp=itmp+nres_molec(i)
+ enddo
+! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization
+ do i=ibond_start,ibond_end
+! cycle
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle ! leave dummy atoms
+ xi=0.5d0*(c(1,i)+c(1,i+1))
+ yi=0.5d0*(c(2,i)+c(2,i+1))
+ zi=0.5d0*(c(3,i)+c(3,i+1))
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+
+ do j=itmp+1,itmp+nres_molec(5)
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+! enddo
+! enddo
+ rcpm = sqrt(xj**2+yj**2+zj**2)
+ drcp_norm(1)=xj/rcpm
+ drcp_norm(2)=yj/rcpm
+ drcp_norm(3)=zj/rcpm
+ dcmag=0.0
+ do k=1,3
+ dcmag=dcmag+dc(k,i)**2
+ enddo
+ dcmag=dsqrt(dcmag)
+ do k=1,3
+ myd_norm(k)=dc(k,i)/dcmag
+ enddo
+ costhet=drcp_norm(1)*myd_norm(1)+drcp_norm(2)*myd_norm(2)+&
+ drcp_norm(3)*myd_norm(3)
+ rsecp = rcpm**2
+ Ir = 1.0d0/rcpm
+ Irsecp = 1.0d0/rsecp
+ Irthrp = Irsecp/rcpm
+ Irfourp = Irthrp/rcpm
+ Irfiftp = Irfourp/rcpm
+ Irsistp=Irfiftp/rcpm
+ Irseven=Irsistp/rcpm
+ Irtwelv=Irsistp*Irsistp
+ Irthir=Irtwelv/rcpm
+ sin2thet = (1-costhet*costhet)
+ sinthet=sqrt(sin2thet)
+ E1 = wdip*Irsecp*costhet+(wmodquad*Irfourp+wquad1*Irthrp)&
+ *sin2thet
+ E2 = -wquad1*Irthrp*wquad2+wvan1*(wvan2**12*Irtwelv-&
+ 2*wvan2**6*Irsistp)
+ ecation_prot = ecation_prot+E1+E2
+ dE1dr = -2*costhet*wdip*Irthrp-&
+ (4*wmodquad*Irfiftp+3*wquad1*Irfourp)*sin2thet
+ dE2dr = 3*wquad1*wquad2*Irfourp- &
+ 12*wvan1*wvan2**6*(wvan2**6*Irthir-Irseven)
+ dEdcos = wdip*Irsecp-2*(wmodquad*Irfourp+wquad1*Irthrp)*costhet
+ do k=1,3
+ drdpep(k) = -drcp_norm(k)
+ dcosdpep(k) = Ir*(costhet*drcp_norm(k)-myd_norm(k))
+ dcosddci(k) = drcp_norm(k)/dcmag-costhet*myd_norm(k)/dcmag
+ dEdpep(k) = (dE1dr+dE2dr)*drdpep(k)+dEdcos*dcosdpep(k)
+ dEddci(k) = dEdcos*dcosddci(k)
+ enddo
+ do k=1,3
+ gradpepcat(k,i)=gradpepcat(k,i)+0.5D0*dEdpep(k)-dEddci(k)
+ gradpepcat(k,i+1)=gradpepcat(k,i+1)+0.5D0*dEdpep(k)+dEddci(k)
+ gradpepcat(k,j)=gradpepcat(k,j)-dEdpep(k)
+ enddo
+ enddo ! j
+ enddo ! i
+!------------------------------------------sidechains
+! do i=1,nres_molec(1)
+ do i=ibond_start,ibond_end
+ if ((itype(i,1).eq.ntyp1)) cycle ! leave dummy atoms
+! cycle
+! print *,i,ecation_prot
+ xi=(c(1,i+nres))
+ yi=(c(2,i+nres))
+ zi=(c(3,i+nres))
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+ do k=1,3
+ cm1(k)=dc(k,i+nres)
+ enddo
+ cm1mag=sqrt(cm1(1)**2+cm1(2)**2+cm1(3)**2)
+ do j=itmp+1,itmp+nres_molec(5)
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+! enddo
+! enddo
+ if(itype(i,1).eq.15.or.itype(i,1).eq.16) then
+ if(itype(i,1).eq.16) then
+ inum=1
+ else
+ inum=2
+ endif
+ do k=1,6
+ vcatprm(k)=catprm(k,inum)
+ enddo
+ dASGL=catprm(7,inum)
+ do k=1,3
+ vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres)
+ valpha(k)=c(k,i)
+ vcat(k)=c(k,j)
+ enddo
+ do k=1,3
+ dx(k) = vcat(k)-vcm(k)
+ enddo
+ do k=1,3
+ v1(k)=(vcm(k)-valpha(k))
+ v2(k)=(vcat(k)-valpha(k))
+ enddo
+ v1m = sqrt(v1(1)**2+v1(2)**2+v1(3)**2)
+ v2m = sqrt(v2(1)**2+v2(2)**2+v2(3)**2)
+ v1dpv2 = v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
+
+! The weights of the energy function calculated from
+!The quantum mechanical GAMESS simulations of calcium with ASP/GLU
+ wh2o=78
+ wc = vcatprm(1)
+ wc=wc/wh2o
+ wdip =vcatprm(2)
+ wdip=wdip/wh2o
+ wquad1 =vcatprm(3)
+ wquad1=wquad1/wh2o
+ wquad2 = vcatprm(4)
+ wquad2=wquad2/wh2o
+ wquad2p = 1-wquad2
+ wvan1 = vcatprm(5)
+ wvan2 =vcatprm(6)
+ opt = dx(1)**2+dx(2)**2
+ rsecp = opt+dx(3)**2
+ rs = sqrt(rsecp)
+ rthrp = rsecp*rs
+ rfourp = rthrp*rs
+ rsixp = rfourp*rsecp
+ reight=rsixp*rsecp
+ Ir = 1.0d0/rs
+ Irsecp = 1/rsecp
+ Irthrp = Irsecp/rs
+ Irfourp = Irthrp/rs
+ Irsixp = 1/rsixp
+ Ireight=1/reight
+ Irtw=Irsixp*Irsixp
+ Irthir=Irtw/rs
+ Irfourt=Irthir/rs
+ opt1 = (4*rs*dx(3)*wdip)
+ opt2 = 6*rsecp*wquad1*opt
+ opt3 = wquad1*wquad2p*Irsixp
+ opt4 = (wvan1*wvan2**12)
+ opt5 = opt4*12*Irfourt
+ opt6 = 2*wvan1*wvan2**6
+ opt7 = 6*opt6*Ireight
+ opt8 = wdip/v1m
+ opt10 = wdip/v2m
+ opt11 = (rsecp*v2m)**2
+ opt12 = (rsecp*v1m)**2
+ opt14 = (v1m*v2m*rsecp)**2
+ opt15 = -wquad1/v2m**2
+ opt16 = (rthrp*(v1m*v2m)**2)**2
+ opt17 = (v1m**2*rthrp)**2
+ opt18 = -wquad1/rthrp
+ opt19 = (v1m**2*v2m**2)**2
+ Ec = wc*Ir
+ do k=1,3
+ dEcCat(k) = -(dx(k)*wc)*Irthrp
+ dEcCm(k)=(dx(k)*wc)*Irthrp
+ dEcCalp(k)=0.0d0
+ enddo
+ Edip=opt8*(v1dpv2)/(rsecp*v2m)
+ do k=1,3
+ dEdipCat(k)=opt8*(v1(k)*rsecp*v2m-((v2(k)/v2m &
+ *rsecp+2*dx(k)*v2m)*v1dpv2))/opt11
+ dEdipCm(k)=opt10*(v2(k)*rsecp*v1m-((v1(k)/v1m &
+ *rsecp-2*dx(k)*v1m)*v1dpv2))/opt12
+ dEdipCalp(k)=wdip*((-v1(k)-v2(k))*rsecp*v1m &
+ *v2m-(-v1(k)/v1m*v2m*rsecp-v2(k)/v2m*v1m*rsecp) &
+ *v1dpv2)/opt14
+ enddo
+ Equad1=-wquad1*v1dpv2**2/(rthrp*(v1m*v2m)**2)
+ do k=1,3
+ dEquad1Cat(k)=-wquad1*(2*v1(k)*v1dpv2*(rthrp* &
+ (v1m*v2m)**2)-(3*dx(k)*rs*(v1m*v2m)**2+2*v1m*2* &
+ v2(k)*1/2*1/v2m*v1m*v2m*rthrp)*v1dpv2**2)/opt16
+ dEquad1Cm(k)=-wquad1*(2*v2(k)*v1dpv2*(rthrp* &
+ (v1m*v2m)**2)-(-3*dx(k)*rs*(v1m*v2m)**2+2*v2m*2* &
+ v1(k)*1/2*1/v1m*v2m*v1m*rthrp)*v1dpv2**2)/opt16
+ dEquad1Calp(k)=opt18*(2*(-v1(k)-v2(k))*v1dpv2* &
+ v1m**2*v2m**2-(-2*v1(k)*v2m**2-2*v2(k)*v1m**2)* &
+ v1dpv2**2)/opt19
+ enddo
+ Equad2=wquad1*wquad2p*Irthrp
+ do k=1,3
+ dEquad2Cat(k)=-3*dx(k)*rs*opt3
+ dEquad2Cm(k)=3*dx(k)*rs*opt3
+ dEquad2Calp(k)=0.0d0
+ enddo
+ Evan1=opt4*Irtw
+ do k=1,3
+ dEvan1Cat(k)=-dx(k)*opt5
+ dEvan1Cm(k)=dx(k)*opt5
+ dEvan1Calp(k)=0.0d0
+ enddo
+ Evan2=-opt6*Irsixp
+ do k=1,3
+ dEvan2Cat(k)=dx(k)*opt7
+ dEvan2Cm(k)=-dx(k)*opt7
+ dEvan2Calp(k)=0.0d0
+ enddo
+ ecation_prot=ecation_prot+Ec+Edip+Equad1+Equad2+Evan1+Evan2
+! print *,ecation_prot,Ec+Edip+Equad1+Equad2+Evan1+Evan2
+
+ do k=1,3
+ dEtotalCat(k)=dEcCat(k)+dEdipCat(k)+dEquad1Cat(k)+ &
+ dEquad2Cat(k)+dEvan1Cat(k)+dEvan2Cat(k)
+!c write(*,*) 'dEtotalCat inside', (dEtotalCat(l),l=1,3)
+ dEtotalCm(k)=dEcCm(k)+dEdipCm(k)+dEquad1Cm(k)+ &
+ dEquad2Cm(k)+dEvan1Cm(k)+dEvan2Cm(k)
+ dEtotalCalp(k)=dEcCalp(k)+dEdipCalp(k)+dEquad1Calp(k) &
+ +dEquad2Calp(k)+dEvan1Calp(k)+dEvan2Calp(k)
+ enddo
+ dscmag = 0.0d0
+ do k=1,3
+ dscvec(k) = dc(k,i+nres)
+ dscmag = dscmag+dscvec(k)*dscvec(k)
+ enddo
+ dscmag3 = dscmag
+ dscmag = sqrt(dscmag)
+ dscmag3 = dscmag3*dscmag
+ constA = 1.0d0+dASGL/dscmag
+ constB = 0.0d0
+ do k=1,3
+ constB = constB+dscvec(k)*dEtotalCm(k)
+ enddo
+ constB = constB*dASGL/dscmag3
+ do k=1,3
+ gg(k) = dEtotalCm(k)+dEtotalCalp(k)
+ gradpepcatx(k,i)=gradpepcatx(k,i)+ &
+ constA*dEtotalCm(k)-constB*dscvec(k)
+! print *,j,constA,dEtotalCm(k),constB,dscvec(k)
+ gradpepcat(k,i)=gradpepcat(k,i)+gg(k)
+ gradpepcat(k,j)=gradpepcat(k,j)+dEtotalCat(k)
+ enddo
+ else if (itype(i,1).eq.13.or.itype(i,1).eq.14) then
+ if(itype(i,1).eq.14) then
+ inum=3
+ else
+ inum=4
+ endif
+ do k=1,6
+ vcatprm(k)=catprm(k,inum)
+ enddo
+ dASGL=catprm(7,inum)
+ do k=1,3
+ vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres)
+ valpha(k)=c(k,i)
+ vcat(k)=c(k,j)
+ enddo
+
+ do k=1,3
+ dx(k) = vcat(k)-vcm(k)
+ enddo
+ do k=1,3
+ v1(k)=(vcm(k)-valpha(k))
+ v2(k)=(vcat(k)-valpha(k))
+ enddo
+ v1m = sqrt(v1(1)**2+v1(2)**2+v1(3)**2)
+ v2m = sqrt(v2(1)**2+v2(2)**2+v2(3)**2)
+ v1dpv2 = v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
+! The weights of the energy function calculated from
+!The quantum mechanical GAMESS simulations of ASN/GLN with calcium
+ wh2o=78
+ wdip =vcatprm(2)
+ wdip=wdip/wh2o
+ wquad1 =vcatprm(3)
+ wquad1=wquad1/wh2o
+ wquad2 = vcatprm(4)
+ wquad2=wquad2/wh2o
+ wquad2p = 1-wquad2
+ wvan1 = vcatprm(5)
+ wvan2 =vcatprm(6)
+ opt = dx(1)**2+dx(2)**2
+ rsecp = opt+dx(3)**2
+ rs = sqrt(rsecp)
+ rthrp = rsecp*rs
+ rfourp = rthrp*rs
+ rsixp = rfourp*rsecp
+ reight=rsixp*rsecp
+ Ir = 1.0d0/rs
+ Irsecp = 1/rsecp
+ Irthrp = Irsecp/rs
+ Irfourp = Irthrp/rs
+ Irsixp = 1/rsixp
+ Ireight=1/reight
+ Irtw=Irsixp*Irsixp
+ Irthir=Irtw/rs
+ Irfourt=Irthir/rs
+ opt1 = (4*rs*dx(3)*wdip)
+ opt2 = 6*rsecp*wquad1*opt
+ opt3 = wquad1*wquad2p*Irsixp
+ opt4 = (wvan1*wvan2**12)
+ opt5 = opt4*12*Irfourt
+ opt6 = 2*wvan1*wvan2**6
+ opt7 = 6*opt6*Ireight
+ opt8 = wdip/v1m
+ opt10 = wdip/v2m
+ opt11 = (rsecp*v2m)**2
+ opt12 = (rsecp*v1m)**2
+ opt14 = (v1m*v2m*rsecp)**2
+ opt15 = -wquad1/v2m**2
+ opt16 = (rthrp*(v1m*v2m)**2)**2
+ opt17 = (v1m**2*rthrp)**2
+ opt18 = -wquad1/rthrp
+ opt19 = (v1m**2*v2m**2)**2
+ Edip=opt8*(v1dpv2)/(rsecp*v2m)
+ do k=1,3
+ dEdipCat(k)=opt8*(v1(k)*rsecp*v2m-((v2(k)/v2m&
+ *rsecp+2*dx(k)*v2m)*v1dpv2))/opt11
+ dEdipCm(k)=opt10*(v2(k)*rsecp*v1m-((v1(k)/v1m&
+ *rsecp-2*dx(k)*v1m)*v1dpv2))/opt12
+ dEdipCalp(k)=wdip*((-v1(k)-v2(k))*rsecp*v1m&
+ *v2m-(-v1(k)/v1m*v2m*rsecp-v2(k)/v2m*v1m*rsecp)&
+ *v1dpv2)/opt14
+ enddo
+ Equad1=-wquad1*v1dpv2**2/(rthrp*(v1m*v2m)**2)
+ do k=1,3
+ dEquad1Cat(k)=-wquad1*(2*v1(k)*v1dpv2*(rthrp*&
+ (v1m*v2m)**2)-(3*dx(k)*rs*(v1m*v2m)**2+2*v1m*2*&
+ v2(k)*1/2*1/v2m*v1m*v2m*rthrp)*v1dpv2**2)/opt16
+ dEquad1Cm(k)=-wquad1*(2*v2(k)*v1dpv2*(rthrp*&
+ (v1m*v2m)**2)-(-3*dx(k)*rs*(v1m*v2m)**2+2*v2m*2*&
+ v1(k)*1/2*1/v1m*v2m*v1m*rthrp)*v1dpv2**2)/opt16
+ dEquad1Calp(k)=opt18*(2*(-v1(k)-v2(k))*v1dpv2* &
+ v1m**2*v2m**2-(-2*v1(k)*v2m**2-2*v2(k)*v1m**2)*&
+ v1dpv2**2)/opt19
+ enddo
+ Equad2=wquad1*wquad2p*Irthrp
+ do k=1,3
+ dEquad2Cat(k)=-3*dx(k)*rs*opt3
+ dEquad2Cm(k)=3*dx(k)*rs*opt3
+ dEquad2Calp(k)=0.0d0
+ enddo
+ Evan1=opt4*Irtw
+ do k=1,3
+ dEvan1Cat(k)=-dx(k)*opt5
+ dEvan1Cm(k)=dx(k)*opt5
+ dEvan1Calp(k)=0.0d0
+ enddo
+ Evan2=-opt6*Irsixp
+ do k=1,3
+ dEvan2Cat(k)=dx(k)*opt7
+ dEvan2Cm(k)=-dx(k)*opt7
+ dEvan2Calp(k)=0.0d0
+ enddo
+ ecation_prot = ecation_prot+Edip+Equad1+Equad2+Evan1+Evan2
+ do k=1,3
+ dEtotalCat(k)=dEdipCat(k)+dEquad1Cat(k)+ &
+ dEquad2Cat(k)+dEvan1Cat(k)+dEvan2Cat(k)
+ dEtotalCm(k)=dEdipCm(k)+dEquad1Cm(k)+ &
+ dEquad2Cm(k)+dEvan1Cm(k)+dEvan2Cm(k)
+ dEtotalCalp(k)=dEdipCalp(k)+dEquad1Calp(k) &
+ +dEquad2Calp(k)+dEvan1Calp(k)+dEvan2Calp(k)
+ enddo
+ dscmag = 0.0d0
+ do k=1,3
+ dscvec(k) = c(k,i+nres)-c(k,i)
+ dscmag = dscmag+dscvec(k)*dscvec(k)
+ enddo
+ dscmag3 = dscmag
+ dscmag = sqrt(dscmag)
+ dscmag3 = dscmag3*dscmag
+ constA = 1+dASGL/dscmag
+ constB = 0.0d0
+ do k=1,3
+ constB = constB+dscvec(k)*dEtotalCm(k)
+ enddo
+ constB = constB*dASGL/dscmag3
+ do k=1,3
+ gg(k) = dEtotalCm(k)+dEtotalCalp(k)
+ gradpepcatx(k,i)=gradpepcatx(k,i)+ &
+ constA*dEtotalCm(k)-constB*dscvec(k)
+ gradpepcat(k,i)=gradpepcat(k,i)+gg(k)
+ gradpepcat(k,j)=gradpepcat(k,j)+dEtotalCat(k)
+ enddo
+ else
+ rcal = 0.0d0
+ do k=1,3
+ r(k) = c(k,j)-c(k,i+nres)
+ rcal = rcal+r(k)*r(k)
+ enddo
+ ract=sqrt(rcal)
+ rocal=1.5
+ epscalc=0.2
+ r0p=0.5*(rocal+sig0(itype(i,1)))
+ r06 = r0p**6
+ r012 = r06*r06
+ Evan1=epscalc*(r012/rcal**6)
+ Evan2=epscalc*2*(r06/rcal**3)
+ r4 = rcal**4
+ r7 = rcal**7
+ do k=1,3
+ dEvan1Cm(k) = 12*r(k)*epscalc*r012/r7
+ dEvan2Cm(k) = 12*r(k)*epscalc*r06/r4
+ enddo
+ do k=1,3
+ dEtotalCm(k)=dEvan1Cm(k)+dEvan2Cm(k)
+ enddo
+ ecation_prot = ecation_prot+ Evan1+Evan2
+ do k=1,3
+ gradpepcatx(k,i)=gradpepcatx(k,i)+ &
+ dEtotalCm(k)
+ gradpepcat(k,i)=gradpepcat(k,i)+dEtotalCm(k)
+ gradpepcat(k,j)=gradpepcat(k,j)-dEtotalCm(k)
+ enddo
+ endif ! 13-16 residues
+ enddo !j
+ enddo !i
+ return
+ end subroutine ecat_prot
+
+!----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+ subroutine eprot_sc_base(escbase)
+ use calc_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.LOCAL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.NAMES'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CALC'
+! include 'COMMON.CONTROL'
+! include 'COMMON.SBRIDGE'
+ logical :: lprn
+!el local variables
+ integer :: iint,itypi,itypi1,itypj,subchap
+ real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
+ real(kind=8) :: evdw,sig0ij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,aa,bb,ssgradlipi,ssgradlipj, &
+ sslipi,sslipj,faclip
+ integer :: ii
+ real(kind=8) :: fracinbuf
+ real (kind=8) :: escbase
+ real (kind=8),dimension(4):: ener
+ real(kind=8) :: b1,b2,b3,b4,egb,eps_in,eps_inout_fac,eps_out
+ real(kind=8) :: ECL,Elj,Equad,Epol,eheadtail,rhead,dGCLOM2,&
+ sqom1,sqom2,sqom12,c1,c2,c3,pom,Lambf,sparrow,&
+ Chif,ChiLambf,bat,eagle,top,bot,botsq,Fcav,dtop,dFdR,dFdOM1,&
+ dFdOM2,w1,w2,w3,dGCLdR,dFdL,dFdOM12,dbot ,&
+ r1,eps_head,alphapol1,pis,facd2,d2,facd1,d1,erdxj,erdxi,federmaus,&
+ dPOLdR1,dFGBdOM2,dFGBdR1,dPOLdFGB1,RR1,MomoFac1,hawk,d1i,d1j,&
+ sig1,sig2,chis12,chis2,ee1,fgb1,a12sq,chis1
+ real(kind=8),dimension(3,2)::chead,erhead_tail
+ real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead
+ integer troll
+ eps_out=80.0d0
+ escbase=0.0d0
+! do i=1,nres_molec(1)
+ do i=ibond_start,ibond_end
+ if (itype(i,1).eq.ntyp1_molec(1)) cycle
+ itypi = itype(i,1)
+ dxi = dc_norm(1,nres+i)
+ dyi = dc_norm(2,nres+i)
+ dzi = dc_norm(3,nres+i)
+ dsci_inv = vbld_inv(i+nres)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+ do j=nres_molec(1)+1,nres_molec(2)+nres_molec(1)
+ itypj= itype(j,2)
+ if (itype(j,2).eq.ntyp1_molec(2))cycle
+ xj=c(1,j+nres)
+ yj=c(2,j+nres)
+ zj=c(3,j+nres)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+ dxj = dc_norm( 1, nres+j )
+ dyj = dc_norm( 2, nres+j )
+ dzj = dc_norm( 3, nres+j )
+! print *,i,j,itypi,itypj
+ d1i = dhead_scbasei(itypi,itypj) !this is shift of dipole/charge
+ d1j = dhead_scbasej(itypi,itypj) !this is shift of dipole/charge
+! d1i=0.0d0
+! d1j=0.0d0
+! BetaT = 1.0d0 / (298.0d0 * Rb)
+! Gay-berne var's
+ sig0ij = sigma_scbase( itypi,itypj )
+ chi1 = chi_scbase( itypi, itypj,1 )
+ chi2 = chi_scbase( itypi, itypj,2 )
+! chi1=0.0d0
+! chi2=0.0d0
+ chi12 = chi1 * chi2
+ chip1 = chipp_scbase( itypi, itypj,1 )
+ chip2 = chipp_scbase( itypi, itypj,2 )
+! chip1=0.0d0
+! chip2=0.0d0
+ chip12 = chip1 * chip2
+! not used by momo potential, but needed by sc_angular which is shared
+! by all energy_potential subroutines
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+ a12sq = rborn_scbasei(itypi,itypj) * rborn_scbasej(itypi,itypj)
+! a12sq = a12sq * a12sq
+! charge of amino acid itypi is...
+ chis1 = chis_scbase(itypi,itypj,1)
+ chis2 = chis_scbase(itypi,itypj,2)
+ chis12 = chis1 * chis2
+ sig1 = sigmap1_scbase(itypi,itypj)
+ sig2 = sigmap2_scbase(itypi,itypj)
+! write (*,*) "sig1 = ", sig1
+! write (*,*) "sig2 = ", sig2
+! alpha factors from Fcav/Gcav
+ b1 = alphasur_scbase(1,itypi,itypj)
+! b1=0.0d0
+ b2 = alphasur_scbase(2,itypi,itypj)
+ b3 = alphasur_scbase(3,itypi,itypj)
+ b4 = alphasur_scbase(4,itypi,itypj)
+! used to determine whether we want to do quadrupole calculations
+! used by Fgb
+ eps_in = epsintab_scbase(itypi,itypj)
+ if (eps_in.eq.0.0) eps_in=1.0
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+! write (*,*) "eps_inout_fac = ", eps_inout_fac
+!-------------------------------------------------------------------
+! tail location and distance calculations
+ DO k = 1,3
+! location of polar head is computed by taking hydrophobic centre
+! and moving by a d1 * dc_norm vector
+! see unres publications for very informative images
+ chead(k,1) = c(k, i+nres) + d1i * dc_norm(k, i+nres)
+ chead(k,2) = c(k, j+nres) + d1j * dc_norm(k, j+nres)
+! distance
+! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ Rhead_distance(k) = chead(k,2) - chead(k,1)
+ END DO
+! pitagoras (root of sum of squares)
+ Rhead = dsqrt( &
+ (Rhead_distance(1)*Rhead_distance(1)) &
+ + (Rhead_distance(2)*Rhead_distance(2)) &
+ + (Rhead_distance(3)*Rhead_distance(3)))
+!-------------------------------------------------------------------
+! zero everything that should be zero'ed
+ evdwij = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ Fcav=0.0d0
+ eheadtail = 0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ Fcav = 0.0d0
+ dFdR = 0.0d0
+ dCAVdOM1 = 0.0d0
+ dCAVdOM2 = 0.0d0
+ dCAVdOM12 = 0.0d0
+ dscj_inv = vbld_inv(j+nres)
+! print *,i,j,dscj_inv,dsci_inv
+! rij holds 1/(distance of Calpha atoms)
+ rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj)
+ rij = dsqrt(rrij)
+!----------------------------
+ CALL sc_angular
+! this should be in elgrad_init but om's are calculated by sc_angular
+! which in turn is used by older potentials
+! om = omega, sqom = om^2
+ sqom1 = om1 * om1
+ sqom2 = om2 * om2
+ sqom12 = om12 * om12
+
+! now we calculate EGB - Gey-Berne
+! It will be summed up in evdwij and saved in evdw
+ sigsq = 1.0D0 / sigsq
+ sig = sig0ij * dsqrt(sigsq)
+! rij_shift = 1.0D0 / rij - sig + sig0ij
+ rij_shift = 1.0/rij - sig + sig0ij
+ IF (rij_shift.le.0.0D0) THEN
+ evdw = 1.0D20
+ RETURN
+ END IF
+ sigder = -sig * sigsq
+ rij_shift = 1.0D0 / rij_shift
+ fac = rij_shift**expon
+ c1 = fac * fac * aa_scbase(itypi,itypj)
+! c1 = 0.0d0
+ c2 = fac * bb_scbase(itypi,itypj)
+! c2 = 0.0d0
+ evdwij = eps1 * eps2rt * eps3rt * ( c1 + c2 )
+ eps2der = eps3rt * evdwij
+ eps3der = eps2rt * evdwij
+! evdwij = 4.0d0 * eps2rt * eps3rt * evdwij
+ evdwij = eps2rt * eps3rt * evdwij
+ c1 = c1 * eps1 * eps2rt**2 * eps3rt**2
+ fac = -expon * (c1 + evdwij) * rij_shift
+ sigder = fac * sigder
+! fac = rij * fac
+! Calculate distance derivative
+ gg(1) = fac
+ gg(2) = fac
+ gg(3) = fac
+! if (b2.gt.0.0) then
+ fac = chis1 * sqom1 + chis2 * sqom2 &
+ - 2.0d0 * chis12 * om1 * om2 * om12
+! we will use pom later in Gcav, so dont mess with it!
+ pom = 1.0d0 - chis1 * chis2 * sqom12
+ Lambf = (1.0d0 - (fac / pom))
+ Lambf = dsqrt(Lambf)
+ sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0)
+! write (*,*) "sparrow = ", sparrow
+ Chif = 1.0d0/rij * sparrow
+ ChiLambf = Chif * Lambf
+ eagle = dsqrt(ChiLambf)
+ bat = ChiLambf ** 11.0d0
+ top = b1 * ( eagle + b2 * ChiLambf - b3 )
+ bot = 1.0d0 + b4 * (ChiLambf ** 12.0d0)
+ botsq = bot * bot
+ Fcav = top / bot
+! print *,i,j,Fcav
+ dtop = b1 * ((Lambf / (2.0d0 * eagle)) + (b2 * Lambf))
+ dbot = 12.0d0 * b4 * bat * Lambf
+ dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow
+! dFdR = 0.0d0
+! write (*,*) "dFcav/dR = ", dFdR
+ dtop = b1 * ((Chif / (2.0d0 * eagle)) + (b2 * Chif))
+ dbot = 12.0d0 * b4 * bat * Chif
+ eagle = Lambf * pom
+ dFdOM1 = -(chis1 * om1 - chis12 * om2 * om12) / (eagle)
+ dFdOM2 = -(chis2 * om2 - chis12 * om1 * om12) / (eagle)
+ dFdOM12 = chis12 * (chis1 * om1 * om12 - om2) &
+ * (chis2 * om2 * om12 - om1) / (eagle * pom)
+
+ dFdL = ((dtop * bot - top * dbot) / botsq)
+! dFdL = 0.0d0
+ dCAVdOM1 = dFdL * ( dFdOM1 )
+ dCAVdOM2 = dFdL * ( dFdOM2 )
+ dCAVdOM12 = dFdL * ( dFdOM12 )
+
+ ertail(1) = xj*rij
+ ertail(2) = yj*rij
+ ertail(3) = zj*rij
+! eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
+! eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
+! eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 &
+! -2.0D0*alf12*eps3der+sigder*sigsq_om12
+! print *,"EOMY",eom1,eom2,eom12
+! erdxi = scalar( ertail(1), dC_norm(1,i+nres) )
+! erdxj = scalar( ertail(1), dC_norm(1,j+nres) )
+! here dtail=0.0
+! facd1 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+! facd2 = dtail(2,itypi,itypj) * vbld_inv(j+nres)
+ DO k = 1, 3
+! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
+! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
+ pom = ertail(k)
+!-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx_scbase(k,i) = gvdwx_scbase(k,i) &
+ - (( dFdR + gg(k) ) * pom)
+! +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
+! +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+! & - ( dFdR * pom )
+ pom = ertail(k)
+!-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx_scbase(k,j) = gvdwx_scbase(k,j) &
+ + (( dFdR + gg(k) ) * pom)
+! +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
+! +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+!c! & + ( dFdR * pom )
+
+ gvdwc_scbase(k,i) = gvdwc_scbase(k,i) &
+ - (( dFdR + gg(k) ) * ertail(k))
+!c! & - ( dFdR * ertail(k))
+
+ gvdwc_scbase(k,j) = gvdwc_scbase(k,j) &
+ + (( dFdR + gg(k) ) * ertail(k))
+!c! & + ( dFdR * ertail(k))
+
+ gg(k) = 0.0d0
+!c! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
+!c! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
+ END DO
+
+! else
+
+! endif
+!Now dipole-dipole
+ if (wdipdip_scbase(2,itypi,itypj).gt.0.0d0) then
+ w1 = wdipdip_scbase(1,itypi,itypj)
+ w2 = -wdipdip_scbase(3,itypi,itypj)/2.0
+ w3 = wdipdip_scbase(2,itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! ECL
+ fac = (om12 - 3.0d0 * om1 * om2)
+ c1 = (w1 / (Rhead**3.0d0)) * fac
+ c2 = (w2 / Rhead ** 6.0d0) &
+ * (4.0d0 + fac * fac -3.0d0 * (sqom1 + sqom2))
+ c3= (w3/ Rhead ** 6.0d0) &
+ * (2.0d0 - 2.0d0*fac*fac +3.0d0*(sqom1 + sqom2))
+ ECL = c1 - c2 + c3
+!c! write (*,*) "w1 = ", w1
+!c! write (*,*) "w2 = ", w2
+!c! write (*,*) "om1 = ", om1
+!c! write (*,*) "om2 = ", om2
+!c! write (*,*) "om12 = ", om12
+!c! write (*,*) "fac = ", fac
+!c! write (*,*) "c1 = ", c1
+!c! write (*,*) "c2 = ", c2
+!c! write (*,*) "Ecl = ", Ecl
+!c! write (*,*) "c2_1 = ", (w2 / Rhead ** 6.0d0)
+!c! write (*,*) "c2_2 = ",
+!c! & (4.0d0 + fac * fac -3.0d0 * (sqom1 + sqom2))
+!c!-------------------------------------------------------------------
+!c! dervative of ECL is GCL...
+!c! dECL/dr
+ c1 = (-3.0d0 * w1 * fac) / (Rhead ** 4.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead ** 7.0d0) &
+ * (4.0d0 + fac * fac - 3.0d0 * (sqom1 + sqom2))
+ c3= (-6.0d0 * w3) / (Rhead ** 7.0d0) &
+ * (2.0d0 - 2.0d0*fac*fac +3.0d0*(sqom1 + sqom2))
+ dGCLdR = c1 - c2 + c3
+!c! dECL/dom1
+ c1 = (-3.0d0 * w1 * om2 ) / (Rhead**3.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead**6.0d0) &
+ * ( om2 * om12 - 3.0d0 * om1 * sqom2 + om1 )
+ c3 =(6.0d0*w3/ Rhead ** 6.0d0)*(om1-2.0d0*(fac)*(-om2))
+ dGCLdOM1 = c1 - c2 + c3
+!c! dECL/dom2
+ c1 = (-3.0d0 * w1 * om1 ) / (Rhead**3.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead**6.0d0) &
+ * ( om1 * om12 - 3.0d0 * sqom1 * om2 + om2 )
+ c3 =(6.0d0*w3/ Rhead ** 6.0d0)*(om2-2.0d0*(fac)*(-om1))
+ dGCLdOM2 = c1 - c2 + c3
+!c! dECL/dom12
+ c1 = w1 / (Rhead ** 3.0d0)
+ c2 = ( 2.0d0 * w2 * fac ) / Rhead ** 6.0d0
+ c3 = (w3/ Rhead ** 6.0d0)*(-4.0d0*fac)
+ dGCLdOM12 = c1 - c2 + c3
+ DO k= 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ END DO
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j+nres) )
+ facd1 = d1i * vbld_inv(i+nres)
+ facd2 = d1j * vbld_inv(j+nres)
+ DO k = 1, 3
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx_scbase(k,i) = gvdwx_scbase(k,i) &
+ - dGCLdR * pom
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx_scbase(k,j) = gvdwx_scbase(k,j) &
+ + dGCLdR * pom
+
+ gvdwc_scbase(k,i) = gvdwc_scbase(k,i) &
+ - dGCLdR * erhead(k)
+ gvdwc_scbase(k,j) = gvdwc_scbase(k,j) &
+ + dGCLdR * erhead(k)
+ END DO
+ endif
+!now charge with dipole eg. ARG-dG
+ if (wqdip_scbase(2,itypi,itypj).gt.0.0d0) then
+ alphapol1 = alphapol_scbase(itypi,itypj)
+ w1 = wqdip_scbase(1,itypi,itypj)
+ w2 = wqdip_scbase(2,itypi,itypj)
+! w1=0.0d0
+! w2=0.0d0
+! pis = sig0head_scbase(itypi,itypj)
+! eps_head = epshead_scbase(itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! R1 - distance between head of ith side chain and tail of jth sidechain
+ R1 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances tail is center of side-chain
+ R1=R1+(c(k,j+nres)-chead(k,1))**2
+ END DO
+!c! Pitagoras
+ R1 = dsqrt(R1)
+
+!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj)
+!c! & +dhead(1,1,itypi,itypj))**2))
+!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj)
+!c! & +dhead(2,1,itypi,itypj))**2))
+
+!c!-------------------------------------------------------------------
+!c! ecl
+ sparrow = w1 * om1
+ hawk = w2 * (1.0d0 - sqom2)
+ Ecl = sparrow / Rhead**2.0d0 &
+ - hawk / Rhead**4.0d0
+!c!-------------------------------------------------------------------
+!c! derivative of ecl is Gcl
+!c! dF/dr part
+ dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 &
+ + 4.0d0 * hawk / Rhead**5.0d0
+!c! dF/dom1
+ dGCLdOM1 = (w1) / (Rhead**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = (2.0d0 * w2 * om2) / (Rhead ** 4.0d0)
+!c--------------------------------------------------------------------
+!c Polarization energy
+!c Epol
+ MomoFac1 = (1.0d0 - chi1 * sqom2)
+ RR1 = R1 * R1 / MomoFac1
+ ee1 = exp(-( RR1 / (4.0d0 * a12sq) ))
+ fgb1 = sqrt( RR1 + a12sq * ee1)
+! eps_inout_fac=0.0d0
+ epol = 332.0d0 * eps_inout_fac * (( alphapol1 / fgb1 )**4.0d0)
+! derivative of Epol is Gpol...
+ dPOLdFGB1 = -(1328.0d0 * eps_inout_fac * alphapol1 ** 4.0d0) &
+ / (fgb1 ** 5.0d0)
+ dFGBdR1 = ( (R1 / MomoFac1) &
+ * ( 2.0d0 - (0.5d0 * ee1) ) ) &
+ / ( 2.0d0 * fgb1 )
+ dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1)) &
+ * (2.0d0 - 0.5d0 * ee1) ) &
+ / (2.0d0 * fgb1)
+ dPOLdR1 = dPOLdFGB1 * dFGBdR1
+! dPOLdR1 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,1) = ((c(k,j+nres)-chead(k,1))/R1)
+ END DO
+
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j+nres) )
+ bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) )
+! bat=0.0d0
+ federmaus = scalar(erhead_tail(1,1),dC_norm(1,j+nres))
+ facd1 = d1i * vbld_inv(i+nres)
+ facd2 = d1j * vbld_inv(j+nres)
+! facd4 = dtail(2,itypi,itypj) * vbld_inv(j+nres)
+
+ DO k = 1, 3
+ hawk = (erhead_tail(k,1) + &
+ facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres)))
+! facd1=0.0d0
+! facd2=0.0d0
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx_scbase(k,i) = gvdwx_scbase(k,i) &
+ - dGCLdR * pom &
+ - dPOLdR1 * (erhead_tail(k,1))
+! & - dGLJdR * pom
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx_scbase(k,j) = gvdwx_scbase(k,j) &
+ + dGCLdR * pom &
+ + dPOLdR1 * (erhead_tail(k,1))
+! & + dGLJdR * pom
+
+
+ gvdwc_scbase(k,i) = gvdwc_scbase(k,i) &
+ - dGCLdR * erhead(k) &
+ - dPOLdR1 * erhead_tail(k,1)
+! & - dGLJdR * erhead(k)
+
+ gvdwc_scbase(k,j) = gvdwc_scbase(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1)
+! & + dGLJdR * erhead(k)
+
+ END DO
+ endif
+! print *,i,j,evdwij,epol,Fcav,ECL
+ escbase=escbase+evdwij+epol+Fcav+ECL
+ call sc_grad_scbase
+ enddo
+ enddo
+
+ return
+ end subroutine eprot_sc_base
+ SUBROUTINE sc_grad_scbase
+ use calc_data
+
+ real (kind=8) :: dcosom1(3),dcosom2(3)
+ eom1 = &
+ eps2der * eps2rt_om1 &
+ - 2.0D0 * alf1 * eps3der &
+ + sigder * sigsq_om1 &
+ + dCAVdOM1 &
+ + dGCLdOM1 &
+ + dPOLdOM1
+
+ eom2 = &
+ eps2der * eps2rt_om2 &
+ + 2.0D0 * alf2 * eps3der &
+ + sigder * sigsq_om2 &
+ + dCAVdOM2 &
+ + dGCLdOM2 &
+ + dPOLdOM2
+
+ eom12 = &
+ evdwij * eps1_om12 &
+ + eps2der * eps2rt_om12 &
+ - 2.0D0 * alf12 * eps3der &
+ + sigder *sigsq_om12 &
+ + dCAVdOM12 &
+ + dGCLdOM12
+
+! print *,eom1,eom2,eom12,i,j,"eom1,2,12",erij(1),erij(2),erij(3)
+! print *,dsci_inv,dscj_inv,dc_norm(2,nres+j),dc_norm(2,nres+i),&
+! gg(1),gg(2),"rozne"
+ DO k = 1, 3
+ dcosom1(k) = rij * (dc_norm(k,nres+i) - om1 * erij(k))
+ dcosom2(k) = rij * (dc_norm(k,nres+j) - om2 * erij(k))
+ gg(k) = gg(k) + eom1 * dcosom1(k) + eom2 * dcosom2(k)
+ gvdwx_scbase(k,i)= gvdwx_scbase(k,i) - gg(k) &
+ + (eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
+ + eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+ gvdwx_scbase(k,j)= gvdwx_scbase(k,j) + gg(k) &
+ + (eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
+ + eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ gvdwc_scbase(k,i)=gvdwc_scbase(k,i)-gg(k)
+ gvdwc_scbase(k,j)=gvdwc_scbase(k,j)+gg(k)
+ END DO
+ RETURN
+ END SUBROUTINE sc_grad_scbase
+
+
+ subroutine epep_sc_base(epepbase)
+ use calc_data
+ logical :: lprn
+!el local variables
+ integer :: iint,itypi,itypi1,itypj,subchap
+ real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
+ real(kind=8) :: evdw,sig0ij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,aa,bb,ssgradlipi,ssgradlipj, &
+ sslipi,sslipj,faclip
+ integer :: ii
+ real(kind=8) :: fracinbuf
+ real (kind=8) :: epepbase
+ real (kind=8),dimension(4):: ener
+ real(kind=8) :: b1,b2,b3,b4,egb,eps_in,eps_inout_fac,eps_out
+ real(kind=8) :: ECL,Elj,Equad,Epol,eheadtail,rhead,dGCLOM2,&
+ sqom1,sqom2,sqom12,c1,c2,c3,pom,Lambf,sparrow,&
+ Chif,ChiLambf,bat,eagle,top,bot,botsq,Fcav,dtop,dFdR,dFdOM1,&
+ dFdOM2,w1,w2,w3,dGCLdR,dFdL,dFdOM12,dbot ,&
+ r1,eps_head,alphapol1,pis,facd2,d2,facd1,d1,erdxj,erdxi,federmaus,&
+ dPOLdR1,dFGBdOM2,dFGBdR1,dPOLdFGB1,RR1,MomoFac1,hawk,d1i,d1j,&
+ sig1,sig2,chis12,chis2,ee1,fgb1,a12sq,chis1
+ real(kind=8),dimension(3,2)::chead,erhead_tail
+ real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead
+ integer troll
+ eps_out=80.0d0
+ epepbase=0.0d0
+! do i=1,nres_molec(1)-1
+ do i=ibond_start,ibond_end
+ if (itype(i,1).eq.ntyp1_molec(1).or.itype(i+1,1).eq.ntyp1_molec(1)) cycle
+!C itypi = itype(i,1)
+ dxi = dc_norm(1,i)
+ dyi = dc_norm(2,i)
+ dzi = dc_norm(3,i)
+! print *,dxi,(-c(1,i)+c(1,i+1))*vbld_inv(i+1)
+ dsci_inv = vbld_inv(i+1)/2.0
+ xi=(c(1,i)+c(1,i+1))/2.0
+ yi=(c(2,i)+c(2,i+1))/2.0
+ zi=(c(3,i)+c(3,i+1))/2.0
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+ do j=nres_molec(1)+1,nres_molec(2)+nres_molec(1)
+ itypj= itype(j,2)
+ if (itype(j,2).eq.ntyp1_molec(2))cycle
+ xj=c(1,j+nres)
+ yj=c(2,j+nres)
+ zj=c(3,j+nres)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+ dxj = dc_norm( 1, nres+j )
+ dyj = dc_norm( 2, nres+j )
+ dzj = dc_norm( 3, nres+j )
+! d1i = dhead_scbasei(itypi) !this is shift of dipole/charge
+! d1j = dhead_scbasej(itypi) !this is shift of dipole/charge
+
+! Gay-berne var's
+ sig0ij = sigma_pepbase(itypj )
+ chi1 = chi_pepbase(itypj,1 )
+ chi2 = chi_pepbase(itypj,2 )
+! chi1=0.0d0
+! chi2=0.0d0
+ chi12 = chi1 * chi2
+ chip1 = chipp_pepbase(itypj,1 )
+ chip2 = chipp_pepbase(itypj,2 )
+! chip1=0.0d0
+! chip2=0.0d0
+ chip12 = chip1 * chip2
+ chis1 = chis_pepbase(itypj,1)
+ chis2 = chis_pepbase(itypj,2)
+ chis12 = chis1 * chis2
+ sig1 = sigmap1_pepbase(itypj)
+ sig2 = sigmap2_pepbase(itypj)
+! write (*,*) "sig1 = ", sig1
+! write (*,*) "sig2 = ", sig2
+ DO k = 1,3
+! location of polar head is computed by taking hydrophobic centre
+! and moving by a d1 * dc_norm vector
+! see unres publications for very informative images
+ chead(k,1) = (c(k,i)+c(k,i+1))/2.0
+! + d1i * dc_norm(k, i+nres)
+ chead(k,2) = c(k, j+nres)
+! + d1j * dc_norm(k, j+nres)
+! distance
+! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ Rhead_distance(k) = chead(k,2) - chead(k,1)
+! print *,gvdwc_pepbase(k,i)
+
+ END DO
+ Rhead = dsqrt( &
+ (Rhead_distance(1)*Rhead_distance(1)) &
+ + (Rhead_distance(2)*Rhead_distance(2)) &
+ + (Rhead_distance(3)*Rhead_distance(3)))
+
+! alpha factors from Fcav/Gcav
+ b1 = alphasur_pepbase(1,itypj)
+! b1=0.0d0
+ b2 = alphasur_pepbase(2,itypj)
+ b3 = alphasur_pepbase(3,itypj)
+ b4 = alphasur_pepbase(4,itypj)
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+ rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj)
+! print *,i,j,rrij
+ rij = dsqrt(rrij)
+!----------------------------
+ evdwij = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ Fcav=0.0d0
+ eheadtail = 0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ Fcav = 0.0d0
+ dFdR = 0.0d0
+ dCAVdOM1 = 0.0d0
+ dCAVdOM2 = 0.0d0
+ dCAVdOM12 = 0.0d0
+ dscj_inv = vbld_inv(j+nres)
+ CALL sc_angular
+! this should be in elgrad_init but om's are calculated by sc_angular
+! which in turn is used by older potentials
+! om = omega, sqom = om^2
+ sqom1 = om1 * om1
+ sqom2 = om2 * om2
+ sqom12 = om12 * om12
+
+! now we calculate EGB - Gey-Berne
+! It will be summed up in evdwij and saved in evdw
+ sigsq = 1.0D0 / sigsq
+ sig = sig0ij * dsqrt(sigsq)
+ rij_shift = 1.0/rij - sig + sig0ij
+ IF (rij_shift.le.0.0D0) THEN
+ evdw = 1.0D20
+ RETURN
+ END IF
+ sigder = -sig * sigsq
+ rij_shift = 1.0D0 / rij_shift
+ fac = rij_shift**expon
+ c1 = fac * fac * aa_pepbase(itypj)
+! c1 = 0.0d0
+ c2 = fac * bb_pepbase(itypj)
+! c2 = 0.0d0
+ evdwij = eps1 * eps2rt * eps3rt * ( c1 + c2 )
+ eps2der = eps3rt * evdwij
+ eps3der = eps2rt * evdwij
+! evdwij = 4.0d0 * eps2rt * eps3rt * evdwij
+ evdwij = eps2rt * eps3rt * evdwij
+ c1 = c1 * eps1 * eps2rt**2 * eps3rt**2
+ fac = -expon * (c1 + evdwij) * rij_shift
+ sigder = fac * sigder
+! fac = rij * fac
+! Calculate distance derivative
+ gg(1) = fac
+ gg(2) = fac
+ gg(3) = fac
+ fac = chis1 * sqom1 + chis2 * sqom2 &
+ - 2.0d0 * chis12 * om1 * om2 * om12
+! we will use pom later in Gcav, so dont mess with it!
+ pom = 1.0d0 - chis1 * chis2 * sqom12
+ Lambf = (1.0d0 - (fac / pom))
+ Lambf = dsqrt(Lambf)
+ sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0)
+! write (*,*) "sparrow = ", sparrow
+ Chif = 1.0d0/rij * sparrow
+ ChiLambf = Chif * Lambf
+ eagle = dsqrt(ChiLambf)
+ bat = ChiLambf ** 11.0d0
+ top = b1 * ( eagle + b2 * ChiLambf - b3 )
+ bot = 1.0d0 + b4 * (ChiLambf ** 12.0d0)
+ botsq = bot * bot
+ Fcav = top / bot
+! print *,i,j,Fcav
+ dtop = b1 * ((Lambf / (2.0d0 * eagle)) + (b2 * Lambf))
+ dbot = 12.0d0 * b4 * bat * Lambf
+ dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow
+! dFdR = 0.0d0
+! write (*,*) "dFcav/dR = ", dFdR
+ dtop = b1 * ((Chif / (2.0d0 * eagle)) + (b2 * Chif))
+ dbot = 12.0d0 * b4 * bat * Chif
+ eagle = Lambf * pom
+ dFdOM1 = -(chis1 * om1 - chis12 * om2 * om12) / (eagle)
+ dFdOM2 = -(chis2 * om2 - chis12 * om1 * om12) / (eagle)
+ dFdOM12 = chis12 * (chis1 * om1 * om12 - om2) &
+ * (chis2 * om2 * om12 - om1) / (eagle * pom)
+
+ dFdL = ((dtop * bot - top * dbot) / botsq)
+! dFdL = 0.0d0
+ dCAVdOM1 = dFdL * ( dFdOM1 )
+ dCAVdOM2 = dFdL * ( dFdOM2 )
+ dCAVdOM12 = dFdL * ( dFdOM12 )
+
+ ertail(1) = xj*rij
+ ertail(2) = yj*rij
+ ertail(3) = zj*rij
+ DO k = 1, 3
+! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
+! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
+ pom = ertail(k)
+!-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres))
+ gvdwc_pepbase(k,i) = gvdwc_pepbase(k,i) &
+ - (( dFdR + gg(k) ) * pom)/2.0
+! print *,gvdwc_pepbase(k,i),i,(( dFdR + gg(k) ) * pom)/2.0
+! +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
+! +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+! & - ( dFdR * pom )
+ pom = ertail(k)
+!-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx_pepbase(k,j) = gvdwx_pepbase(k,j) &
+ + (( dFdR + gg(k) ) * pom)
+! +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
+! +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+!c! & + ( dFdR * pom )
+
+ gvdwc_pepbase(k,i+1) = gvdwc_pepbase(k,i+1) &
+ - (( dFdR + gg(k) ) * ertail(k))/2.0
+! print *,gvdwc_pepbase(k,i+1),i+1,(( dFdR + gg(k) ) * pom)/2.0
+
+!c! & - ( dFdR * ertail(k))
+
+ gvdwc_pepbase(k,j) = gvdwc_pepbase(k,j) &
+ + (( dFdR + gg(k) ) * ertail(k))
+!c! & + ( dFdR * ertail(k))
+
+ gg(k) = 0.0d0
+!c! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
+!c! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
+ END DO
+
+
+ w1 = wdipdip_pepbase(1,itypj)
+ w2 = -wdipdip_pepbase(3,itypj)/2.0
+ w3 = wdipdip_pepbase(2,itypj)
+! w1=0.0d0
+! w2=0.0d0
+!c!-------------------------------------------------------------------
+!c! ECL
+! w3=0.0d0
+ fac = (om12 - 3.0d0 * om1 * om2)
+ c1 = (w1 / (Rhead**3.0d0)) * fac
+ c2 = (w2 / Rhead ** 6.0d0) &
+ * (4.0d0 + fac * fac -3.0d0 * (sqom1 + sqom2))
+ c3= (w3/ Rhead ** 6.0d0) &
+ * (2.0d0 - 2.0d0*fac*fac +3.0d0*(sqom1 + sqom2))
+
+ ECL = c1 - c2 + c3
+
+ c1 = (-3.0d0 * w1 * fac) / (Rhead ** 4.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead ** 7.0d0) &
+ * (4.0d0 + fac * fac - 3.0d0 * (sqom1 + sqom2))
+ c3= (-6.0d0 * w3) / (Rhead ** 7.0d0) &
+ * (2.0d0 - 2.0d0*fac*fac +3.0d0*(sqom1 + sqom2))
+
+ dGCLdR = c1 - c2 + c3
+!c! dECL/dom1
+ c1 = (-3.0d0 * w1 * om2 ) / (Rhead**3.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead**6.0d0) &
+ * ( om2 * om12 - 3.0d0 * om1 * sqom2 + om1 )
+ c3 =(6.0d0*w3/ Rhead ** 6.0d0)*(om1-2.0d0*(fac)*(-om2))
+ dGCLdOM1 = c1 - c2 + c3
+!c! dECL/dom2
+ c1 = (-3.0d0 * w1 * om1 ) / (Rhead**3.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead**6.0d0) &
+ * ( om1 * om12 - 3.0d0 * sqom1 * om2 + om2 )
+ c3 =(6.0d0*w3/ Rhead ** 6.0d0)*(om2-2.0d0*(fac)*(-om1))
+
+ dGCLdOM2 = c1 - c2 + c3
+!c! dECL/dom12
+ c1 = w1 / (Rhead ** 3.0d0)
+ c2 = ( 2.0d0 * w2 * fac ) / Rhead ** 6.0d0
+ c3 = (w3/ Rhead ** 6.0d0)*(-4.0d0*fac)
+ dGCLdOM12 = c1 - c2 + c3
+ DO k= 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ END DO
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j+nres) )
+! facd1 = d1 * vbld_inv(i+nres)
+! facd2 = d2 * vbld_inv(j+nres)
+ DO k = 1, 3
+
+! pom = erhead(k)
+!+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+! gvdwx_pepbase(k,i) = gvdwx_scbase(k,i) &
+! - dGCLdR * pom
+ pom = erhead(k)
+!+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx_pepbase(k,j) = gvdwx_pepbase(k,j) &
+ + dGCLdR * pom
+
+ gvdwc_pepbase(k,i) = gvdwc_pepbase(k,i) &
+ - dGCLdR * erhead(k)/2.0d0
+! print *,gvdwc_pepbase(k,i+1),i+1,- dGCLdR * erhead(k)/2.0d0
+ gvdwc_pepbase(k,i+1) = gvdwc_pepbase(k,i+1) &
+ - dGCLdR * erhead(k)/2.0d0
+! print *,gvdwc_pepbase(k,i+1),i+1,- dGCLdR * erhead(k)/2.0d0
+ gvdwc_pepbase(k,j) = gvdwc_pepbase(k,j) &
+ + dGCLdR * erhead(k)
+ END DO
+! print *,i,j,evdwij,Fcav,ECL,"vdw,cav,ecl"
+ epepbase=epepbase+evdwij+Fcav+ECL
+ call sc_grad_pepbase
+ enddo
+ enddo
+ END SUBROUTINE epep_sc_base
+ SUBROUTINE sc_grad_pepbase
+ use calc_data
+
+ real (kind=8) :: dcosom1(3),dcosom2(3)
+ eom1 = &
+ eps2der * eps2rt_om1 &
+ - 2.0D0 * alf1 * eps3der &
+ + sigder * sigsq_om1 &
+ + dCAVdOM1 &
+ + dGCLdOM1 &
+ + dPOLdOM1
+
+ eom2 = &
+ eps2der * eps2rt_om2 &
+ + 2.0D0 * alf2 * eps3der &
+ + sigder * sigsq_om2 &
+ + dCAVdOM2 &
+ + dGCLdOM2 &
+ + dPOLdOM2
+
+ eom12 = &
+ evdwij * eps1_om12 &
+ + eps2der * eps2rt_om12 &
+ - 2.0D0 * alf12 * eps3der &
+ + sigder *sigsq_om12 &
+ + dCAVdOM12 &
+ + dGCLdOM12
+! om12=0.0
+! eom12=0.0
+! print *,eom1,eom2,eom12,om12,i,j,"eom1,2,12",erij(1),erij(2),erij(3)
+! if (i.eq.30) print *,gvdwc_pepbase(k,i),- gg(k),&
+! (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i)))&
+! *dsci_inv*2.0
+! print *,dsci_inv,dscj_inv,dc_norm(2,nres+j),dc_norm(2,nres+i),&
+! gg(1),gg(2),"rozne"
+ DO k = 1, 3
+ dcosom1(k) = rij * (dc_norm(k,i) - om1 * erij(k))
+ dcosom2(k) = rij * (dc_norm(k,nres+j) - om2 * erij(k))
+ gg(k) = gg(k) + eom1 * dcosom1(k) + eom2 * dcosom2(k)
+ gvdwc_pepbase(k,i)= gvdwc_pepbase(k,i) +0.5*(- gg(k)) &
+ + (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i)))&
+ *dsci_inv*2.0 &
+ - (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0
+ gvdwc_pepbase(k,i+1)= gvdwc_pepbase(k,i+1) +0.5*(- gg(k)) &
+ - (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i))) &
+ *dsci_inv*2.0 &
+ + (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0
+! print *,eom12,eom2,om12,om2
+!eom12*(-dc_norm(k,i)/2.0-om12*dc_norm(k,nres+j)),&
+! (eom2*(erij(k)-om2*dc_norm(k,nres+j)))
+ gvdwx_pepbase(k,j)= gvdwx_pepbase(k,j) + gg(k) &
+ + (eom12*(dc_norm(k,i)-om12*dc_norm(k,nres+j))&
+ + eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ gvdwc_pepbase(k,j)=gvdwc_pepbase(k,j)+gg(k)
+ END DO
+ RETURN
+ END SUBROUTINE sc_grad_pepbase
+ subroutine eprot_sc_phosphate(escpho)
+ use calc_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.LOCAL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.NAMES'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CALC'
+! include 'COMMON.CONTROL'
+! include 'COMMON.SBRIDGE'
+ logical :: lprn
+!el local variables
+ integer :: iint,itypi,itypi1,itypj,subchap
+ real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
+ real(kind=8) :: evdw,sig0ij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,aa,bb,ssgradlipi,ssgradlipj, &
+ sslipi,sslipj,faclip,alpha_sco
+ integer :: ii
+ real(kind=8) :: fracinbuf
+ real (kind=8) :: escpho
+ real (kind=8),dimension(4):: ener
+ real(kind=8) :: b1,b2,b3,b4,egb,eps_in,eps_inout_fac,eps_out
+ real(kind=8) :: ECL,Elj,Equad,Epol,eheadtail,rhead,dGCLOM2,&
+ sqom1,sqom2,sqom12,c1,c2,pom,Lambf,sparrow,&
+ Chif,ChiLambf,bat,eagle,top,bot,botsq,Fcav,dtop,dFdR,dFdOM1,&
+ dFdOM2,w1,w2,dGCLdR,dFdL,dFdOM12,dbot ,&
+ r1,eps_head,alphapol1,pis,facd2,d2,facd1,d1,erdxj,erdxi,federmaus,&
+ dPOLdR1,dFGBdOM2,dFGBdR1,dPOLdFGB1,RR1,MomoFac1,hawk,d1i,d1j,&
+ sig1,sig2,chis12,chis2,ee1,fgb1,a12sq,chis1,Rhead_sq,Qij,dFGBdOM1
+ real(kind=8),dimension(3,2)::chead,erhead_tail
+ real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead
+ integer troll
+ eps_out=80.0d0
+ escpho=0.0d0
+! do i=1,nres_molec(1)
+ do i=ibond_start,ibond_end
+ if (itype(i,1).eq.ntyp1_molec(1)) cycle
+ itypi = itype(i,1)
+ dxi = dc_norm(1,nres+i)
+ dyi = dc_norm(2,nres+i)
+ dzi = dc_norm(3,nres+i)
+ dsci_inv = vbld_inv(i+nres)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+ do j=nres_molec(1)+1,nres_molec(2)+nres_molec(1)-1
+ itypj= itype(j,2)
+ if ((itype(j,2).eq.ntyp1_molec(2)).or.&
+ (itype(j+1,2).eq.ntyp1_molec(2))) cycle
+ xj=(c(1,j)+c(1,j+1))/2.0
+ yj=(c(2,j)+c(2,j+1))/2.0
+ zj=(c(3,j)+c(3,j+1))/2.0
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+ dxj = dc_norm( 1,j )
+ dyj = dc_norm( 2,j )
+ dzj = dc_norm( 3,j )
+ dscj_inv = vbld_inv(j+1)
+
+! Gay-berne var's
+ sig0ij = sigma_scpho(itypi )
+ chi1 = chi_scpho(itypi,1 )
+ chi2 = chi_scpho(itypi,2 )
+! chi1=0.0d0
+! chi2=0.0d0
+ chi12 = chi1 * chi2
+ chip1 = chipp_scpho(itypi,1 )
+ chip2 = chipp_scpho(itypi,2 )
+! chip1=0.0d0
+! chip2=0.0d0
+ chip12 = chip1 * chip2
+ chis1 = chis_scpho(itypi,1)
+ chis2 = chis_scpho(itypi,2)
+ chis12 = chis1 * chis2
+ sig1 = sigmap1_scpho(itypi)
+ sig2 = sigmap2_scpho(itypi)
+! write (*,*) "sig1 = ", sig1
+! write (*,*) "sig1 = ", sig1
+! write (*,*) "sig2 = ", sig2
+! alpha factors from Fcav/Gcav
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+ a12sq = rborn_scphoi(itypi) * rborn_scphoj(itypi)
+
+ b1 = alphasur_scpho(1,itypi)
+! b1=0.0d0
+ b2 = alphasur_scpho(2,itypi)
+ b3 = alphasur_scpho(3,itypi)
+ b4 = alphasur_scpho(4,itypi)
+! used to determine whether we want to do quadrupole calculations
+! used by Fgb
+ eps_in = epsintab_scpho(itypi)
+ if (eps_in.eq.0.0) eps_in=1.0
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+! write (*,*) "eps_inout_fac = ", eps_inout_fac
+!-------------------------------------------------------------------
+! tail location and distance calculations
+ d1i = dhead_scphoi(itypi) !this is shift of dipole/charge
+ d1j = 0.0
+ DO k = 1,3
+! location of polar head is computed by taking hydrophobic centre
+! and moving by a d1 * dc_norm vector
+! see unres publications for very informative images
+ chead(k,1) = c(k, i+nres) + d1i * dc_norm(k, i+nres)
+ chead(k,2) = (c(k, j) + c(k, j+1))/2.0
+! distance
+! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ Rhead_distance(k) = chead(k,2) - chead(k,1)
+ END DO
+! pitagoras (root of sum of squares)
+ Rhead = dsqrt( &
+ (Rhead_distance(1)*Rhead_distance(1)) &
+ + (Rhead_distance(2)*Rhead_distance(2)) &
+ + (Rhead_distance(3)*Rhead_distance(3)))
+ Rhead_sq=Rhead**2.0
+!-------------------------------------------------------------------
+! zero everything that should be zero'ed
+ evdwij = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ Fcav=0.0d0
+ eheadtail = 0.0d0
+ dGCLdR=0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ Fcav = 0.0d0
+ dFdR = 0.0d0
+ dCAVdOM1 = 0.0d0
+ dCAVdOM2 = 0.0d0
+ dCAVdOM12 = 0.0d0
+ dscj_inv = vbld_inv(j+1)/2.0
+!dhead_scbasej(itypi,itypj)
+! print *,i,j,dscj_inv,dsci_inv
+! rij holds 1/(distance of Calpha atoms)
+ rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj)
+ rij = dsqrt(rrij)
+!----------------------------
+ CALL sc_angular
+! this should be in elgrad_init but om's are calculated by sc_angular
+! which in turn is used by older potentials
+! om = omega, sqom = om^2
+ sqom1 = om1 * om1
+ sqom2 = om2 * om2
+ sqom12 = om12 * om12
+
+! now we calculate EGB - Gey-Berne
+! It will be summed up in evdwij and saved in evdw
+ sigsq = 1.0D0 / sigsq
+ sig = sig0ij * dsqrt(sigsq)
+! rij_shift = 1.0D0 / rij - sig + sig0ij
+ rij_shift = 1.0/rij - sig + sig0ij
+ IF (rij_shift.le.0.0D0) THEN
+ evdw = 1.0D20
+ RETURN
+ END IF
+ sigder = -sig * sigsq
+ rij_shift = 1.0D0 / rij_shift
+ fac = rij_shift**expon
+ c1 = fac * fac * aa_scpho(itypi)
+! c1 = 0.0d0
+ c2 = fac * bb_scpho(itypi)
+! c2 = 0.0d0
+ evdwij = eps1 * eps2rt * eps3rt * ( c1 + c2 )
+ eps2der = eps3rt * evdwij
+ eps3der = eps2rt * evdwij
+! evdwij = 4.0d0 * eps2rt * eps3rt * evdwij
+ evdwij = eps2rt * eps3rt * evdwij
+ c1 = c1 * eps1 * eps2rt**2 * eps3rt**2
+ fac = -expon * (c1 + evdwij) * rij_shift
+ sigder = fac * sigder
+! fac = rij * fac
+! Calculate distance derivative
+ gg(1) = fac
+ gg(2) = fac
+ gg(3) = fac
+ fac = chis1 * sqom1 + chis2 * sqom2 &
+ - 2.0d0 * chis12 * om1 * om2 * om12
+! we will use pom later in Gcav, so dont mess with it!
+ pom = 1.0d0 - chis1 * chis2 * sqom12
+ Lambf = (1.0d0 - (fac / pom))
+ Lambf = dsqrt(Lambf)
+ sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0)
+! write (*,*) "sparrow = ", sparrow
+ Chif = 1.0d0/rij * sparrow
+ ChiLambf = Chif * Lambf
+ eagle = dsqrt(ChiLambf)
+ bat = ChiLambf ** 11.0d0
+ top = b1 * ( eagle + b2 * ChiLambf - b3 )
+ bot = 1.0d0 + b4 * (ChiLambf ** 12.0d0)
+ botsq = bot * bot
+ Fcav = top / bot
+ dtop = b1 * ((Lambf / (2.0d0 * eagle)) + (b2 * Lambf))
+ dbot = 12.0d0 * b4 * bat * Lambf
+ dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow
+! dFdR = 0.0d0
+! write (*,*) "dFcav/dR = ", dFdR
+ dtop = b1 * ((Chif / (2.0d0 * eagle)) + (b2 * Chif))
+ dbot = 12.0d0 * b4 * bat * Chif
+ eagle = Lambf * pom
+ dFdOM1 = -(chis1 * om1 - chis12 * om2 * om12) / (eagle)
+ dFdOM2 = -(chis2 * om2 - chis12 * om1 * om12) / (eagle)
+ dFdOM12 = chis12 * (chis1 * om1 * om12 - om2) &
+ * (chis2 * om2 * om12 - om1) / (eagle * pom)
+
+ dFdL = ((dtop * bot - top * dbot) / botsq)
+! dFdL = 0.0d0
+ dCAVdOM1 = dFdL * ( dFdOM1 )
+ dCAVdOM2 = dFdL * ( dFdOM2 )
+ dCAVdOM12 = dFdL * ( dFdOM12 )
+
+ ertail(1) = xj*rij
+ ertail(2) = yj*rij
+ ertail(3) = zj*rij
+ DO k = 1, 3
+! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
+! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
+! if (i.eq.3) print *,'decl0',gvdwx_scpho(k,i),i
+
+ pom = ertail(k)
+! print *,pom,gg(k),dFdR
+!-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx_scpho(k,i) = gvdwx_scpho(k,i) &
+ - (( dFdR + gg(k) ) * pom)
+! +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
+! +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+! & - ( dFdR * pom )
+! pom = ertail(k)
+!-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres))
+! gvdwx_scpho(k,j) = gvdwx_scpho(k,j) &
+! + (( dFdR + gg(k) ) * pom)
+! +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
+! +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+!c! & + ( dFdR * pom )
+
+ gvdwc_scpho(k,i) = gvdwc_scpho(k,i) &
+ - (( dFdR + gg(k) ) * ertail(k))
+!c! & - ( dFdR * ertail(k))
+
+ gvdwc_scpho(k,j) = gvdwc_scpho(k,j) &
+ + (( dFdR + gg(k) ) * ertail(k))/2.0
+
+ gvdwc_scpho(k,j+1) = gvdwc_scpho(k,j+1) &
+ + (( dFdR + gg(k) ) * ertail(k))/2.0
+
+!c! & + ( dFdR * ertail(k))
+
+ gg(k) = 0.0d0
+ ENDDO
+!c! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
+!c! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
+! alphapol1 = alphapol_scpho(itypi)
+ if (wqq_scpho(itypi).ne.0.0) then
+ Qij=wqq_scpho(itypi)/eps_in
+ alpha_sco=1.d0/alphi_scpho(itypi)
+! Qij=0.0
+ Ecl = (332.0d0 * Qij*dexp(-Rhead*alpha_sco)) / Rhead
+!c! derivative of Ecl is Gcl...
+ dGCLdR = (-332.0d0 * Qij*dexp(-Rhead*alpha_sco)* &
+ (Rhead*alpha_sco+1) ) / Rhead_sq
+ if (energy_dec) write(iout,*) "ECL",ECL,Rhead,1.0/rij
+ else if (wqdip_scpho(2,itypi).gt.0.0d0) then
+ w1 = wqdip_scpho(1,itypi)
+ w2 = wqdip_scpho(2,itypi)
+! w1=0.0d0
+! w2=0.0d0
+! pis = sig0head_scbase(itypi,itypj)
+! eps_head = epshead_scbase(itypi,itypj)
+!c!-------------------------------------------------------------------
+
+!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj)
+!c! & +dhead(1,1,itypi,itypj))**2))
+!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj)
+!c! & +dhead(2,1,itypi,itypj))**2))
+
+!c!-------------------------------------------------------------------
+!c! ecl
+ sparrow = w1 * om1
+ hawk = w2 * (1.0d0 - sqom2)
+ Ecl = sparrow / Rhead**2.0d0 &
+ - hawk / Rhead**4.0d0
+!c!-------------------------------------------------------------------
+ if (energy_dec) write(iout,*) "ECLdipdip",ECL,Rhead,&
+ 1.0/rij,sparrow
+
+!c! derivative of ecl is Gcl
+!c! dF/dr part
+ dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 &
+ + 4.0d0 * hawk / Rhead**5.0d0
+!c! dF/dom1
+ dGCLdOM1 = (w1) / (Rhead**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = (2.0d0 * w2 * om2) / (Rhead ** 4.0d0)
+ endif
+
+!c--------------------------------------------------------------------
+!c Polarization energy
+!c Epol
+ R1 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances tail is center of side-chain
+ R1=R1+((c(k,j)+c(k,j+1))/2.0-chead(k,1))**2
+ END DO
+!c! Pitagoras
+ R1 = dsqrt(R1)
+
+ alphapol1 = alphapol_scpho(itypi)
+! alphapol1=0.0
+ MomoFac1 = (1.0d0 - chi2 * sqom1)
+ RR1 = R1 * R1 / MomoFac1
+ ee1 = exp(-( RR1 / (4.0d0 * a12sq) ))
+! print *,"ee1",ee1,a12sq,alphapol1,eps_inout_fac
+ fgb1 = sqrt( RR1 + a12sq * ee1)
+! eps_inout_fac=0.0d0
+ epol = 332.0d0 * eps_inout_fac * (( alphapol1 / fgb1 )**4.0d0)
+! derivative of Epol is Gpol...
+ dPOLdFGB1 = -(1328.0d0 * eps_inout_fac * alphapol1 ** 4.0d0) &
+ / (fgb1 ** 5.0d0)
+ dFGBdR1 = ( (R1 / MomoFac1) &
+ * ( 2.0d0 - (0.5d0 * ee1) ) ) &
+ / ( 2.0d0 * fgb1 )
+ dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1)) &
+ * (2.0d0 - 0.5d0 * ee1) ) &
+ / (2.0d0 * fgb1)
+ dPOLdR1 = dPOLdFGB1 * dFGBdR1
+! dPOLdR1 = 0.0d0
+! dPOLdOM1 = 0.0d0
+ dFGBdOM1 = (((R1 * R1 * chi2 * om1) / (MomoFac1 * MomoFac1)) &
+ * (2.0d0 - 0.5d0 * ee1) ) &
+ / (2.0d0 * fgb1)
+
+ dPOLdOM1 = dPOLdFGB1 * dFGBdOM1
+ dPOLdOM2 = 0.0
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,1) = (((c(k,j)+c(k,j+1))/2.0-chead(k,1))/R1)
+ END DO
+
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j) )
+ bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) )
+! bat=0.0d0
+ federmaus = scalar(erhead_tail(1,1),dC_norm(1,j))
+ facd1 = d1i * vbld_inv(i+nres)
+ facd2 = d1j * vbld_inv(j)
+! facd4 = dtail(2,itypi,itypj) * vbld_inv(j+nres)
+
+ DO k = 1, 3
+ hawk = (erhead_tail(k,1) + &
+ facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres)))
+! facd1=0.0d0
+! facd2=0.0d0
+! if (i.eq.3) print *,'decl1',dGCLdR,dPOLdR1,gvdwc_scpho(k,i),i,&
+! pom,(erhead_tail(k,1))
+
+! print *,'decl',dGCLdR,dPOLdR1,gvdwc_scpho(k,i)
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx_scpho(k,i) = gvdwx_scpho(k,i) &
+ - dGCLdR * pom &
+ - dPOLdR1 * (erhead_tail(k,1))
+! & - dGLJdR * pom
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
+! gvdwx_scpho(k,j) = gvdwx_scpho(k,j) &
+! + dGCLdR * pom &
+! + dPOLdR1 * (erhead_tail(k,1))
+! & + dGLJdR * pom
+
+
+ gvdwc_scpho(k,i) = gvdwc_scpho(k,i) &
+ - dGCLdR * erhead(k) &
+ - dPOLdR1 * erhead_tail(k,1)
+! & - dGLJdR * erhead(k)
+
+ gvdwc_scpho(k,j) = gvdwc_scpho(k,j) &
+ + (dGCLdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1))/2.0
+ gvdwc_scpho(k,j+1) = gvdwc_scpho(k,j+1) &
+ + (dGCLdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1))/2.0
+
+! & + dGLJdR * erhead(k)
+! if (i.eq.3) print *,'decl2',dGCLdR,dPOLdR1,gvdwc_scpho(k,i),i
+
+ END DO
+! if (i.eq.3) print *,i,j,evdwij,epol,Fcav,ECL
+ if (energy_dec) write (iout,'(a22,2i5,4f8.3,f16.3)'), &
+ "escpho:evdw,pol,cav,CL",i,j,evdwij,epol,Fcav,ECL,escpho
+ escpho=escpho+evdwij+epol+Fcav+ECL
+ call sc_grad_scpho
+ enddo
+
+ enddo
+
+ return
+ end subroutine eprot_sc_phosphate
+ SUBROUTINE sc_grad_scpho
+ use calc_data
+
+ real (kind=8) :: dcosom1(3),dcosom2(3)
+ eom1 = &
+ eps2der * eps2rt_om1 &
+ - 2.0D0 * alf1 * eps3der &
+ + sigder * sigsq_om1 &
+ + dCAVdOM1 &
+ + dGCLdOM1 &
+ + dPOLdOM1
+
+ eom2 = &
+ eps2der * eps2rt_om2 &
+ + 2.0D0 * alf2 * eps3der &
+ + sigder * sigsq_om2 &
+ + dCAVdOM2 &
+ + dGCLdOM2 &
+ + dPOLdOM2
+
+ eom12 = &
+ evdwij * eps1_om12 &
+ + eps2der * eps2rt_om12 &
+ - 2.0D0 * alf12 * eps3der &
+ + sigder *sigsq_om12 &
+ + dCAVdOM12 &
+ + dGCLdOM12
+! om12=0.0
+! eom12=0.0
+! print *,eom1,eom2,eom12,om12,i,j,"eom1,2,12",erij(1),erij(2),erij(3)
+! if (i.eq.30) print *,gvdwc_scpho(k,i),- gg(k),&
+! (-eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,i)))&
+! *dsci_inv*2.0
+! print *,dsci_inv,dscj_inv,dc_norm(2,nres+j),dc_norm(2,nres+i),&
+! gg(1),gg(2),"rozne"
+ DO k = 1, 3
+ dcosom1(k) = rij * (dc_norm(k,nres+i) - om1 * erij(k))
+ dcosom2(k) = rij * (dc_norm(k,j) - om2 * erij(k))
+ gg(k) = gg(k) + eom1 * dcosom1(k) + eom2 * dcosom2(k)
+ gvdwc_scpho(k,j)= gvdwc_scpho(k,j) +0.5*( gg(k)) &
+ + (-eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,j)))&
+ *dscj_inv*2.0 &
+ - (eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv*2.0
+ gvdwc_scpho(k,j+1)= gvdwc_scpho(k,j+1) +0.5*( gg(k)) &
+ - (-eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,j))) &
+ *dscj_inv*2.0 &
+ + (eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv*2.0
+ gvdwx_scpho(k,i)= gvdwx_scpho(k,i) - gg(k) &
+ + (eom12*(dc_norm(k,j)-om12*dc_norm(k,nres+i)) &
+ + eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+
+! print *,eom12,eom2,om12,om2
+!eom12*(-dc_norm(k,i)/2.0-om12*dc_norm(k,nres+j)),&
+! (eom2*(erij(k)-om2*dc_norm(k,nres+j)))
+! gvdwx_scpho(k,j)= gvdwx_scpho(k,j) + gg(k) &
+! + (eom12*(dc_norm(k,i)-om12*dc_norm(k,nres+j))&
+! + eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ gvdwc_scpho(k,i)=gvdwc_scpho(k,i)-gg(k)
+ END DO
+ RETURN
+ END SUBROUTINE sc_grad_scpho
+ subroutine eprot_pep_phosphate(epeppho)
+ use calc_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.LOCAL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.NAMES'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CALC'
+! include 'COMMON.CONTROL'
+! include 'COMMON.SBRIDGE'
+ logical :: lprn
+!el local variables
+ integer :: iint,itypi,itypi1,itypj,subchap
+ real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
+ real(kind=8) :: evdw,sig0ij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,aa,bb,ssgradlipi,ssgradlipj, &
+ sslipi,sslipj,faclip
+ integer :: ii
+ real(kind=8) :: fracinbuf
+ real (kind=8) :: epeppho
+ real (kind=8),dimension(4):: ener
+ real(kind=8) :: b1,b2,b3,b4,egb,eps_in,eps_inout_fac,eps_out
+ real(kind=8) :: ECL,Elj,Equad,Epol,eheadtail,rhead,dGCLOM2,&
+ sqom1,sqom2,sqom12,c1,c2,pom,Lambf,sparrow,&
+ Chif,ChiLambf,bat,eagle,top,bot,botsq,Fcav,dtop,dFdR,dFdOM1,&
+ dFdOM2,w1,w2,dGCLdR,dFdL,dFdOM12,dbot ,&
+ r1,eps_head,alphapol1,pis,facd2,d2,facd1,d1,erdxj,erdxi,federmaus,&
+ dPOLdR1,dFGBdOM2,dFGBdR1,dPOLdFGB1,RR1,MomoFac1,hawk,d1i,d1j,&
+ sig1,sig2,chis12,chis2,ee1,fgb1,a12sq,chis1,Rhead_sq,Qij,dFGBdOM1
+ real(kind=8),dimension(3,2)::chead,erhead_tail
+ real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead
+ integer troll
+ real (kind=8) :: dcosom1(3),dcosom2(3)
+ epeppho=0.0d0
+! do i=1,nres_molec(1)
+ do i=ibond_start,ibond_end
+ if (itype(i,1).eq.ntyp1_molec(1)) cycle
+ itypi = itype(i,1)
+ dsci_inv = vbld_inv(i+1)/2.0
+ dxi = dc_norm(1,i)
+ dyi = dc_norm(2,i)
+ dzi = dc_norm(3,i)
+ xi=(c(1,i)+c(1,i+1))/2.0
+ yi=(c(2,i)+c(2,i+1))/2.0
+ zi=(c(3,i)+c(3,i+1))/2.0
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+ do j=nres_molec(1)+1,nres_molec(2)+nres_molec(1)-1
+ itypj= itype(j,2)
+ if ((itype(j,2).eq.ntyp1_molec(2)).or.&
+ (itype(j+1,2).eq.ntyp1_molec(2))) cycle
+ xj=(c(1,j)+c(1,j+1))/2.0
+ yj=(c(2,j)+c(2,j+1))/2.0
+ zj=(c(3,j)+c(3,j+1))/2.0
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+ rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj)
+ rij = dsqrt(rrij)
+ dxj = dc_norm( 1,j )
+ dyj = dc_norm( 2,j )
+ dzj = dc_norm( 3,j )
+ dscj_inv = vbld_inv(j+1)/2.0
+! Gay-berne var's
+ sig0ij = sigma_peppho
+! chi1=0.0d0
+! chi2=0.0d0
+ chi12 = chi1 * chi2
+! chip1=0.0d0
+! chip2=0.0d0
+ chip12 = chip1 * chip2
+! chis1 = 0.0d0
+! chis2 = 0.0d0
+ chis12 = chis1 * chis2
+ sig1 = sigmap1_peppho
+ sig2 = sigmap2_peppho
+! write (*,*) "sig1 = ", sig1
+! write (*,*) "sig1 = ", sig1
+! write (*,*) "sig2 = ", sig2
+! alpha factors from Fcav/Gcav
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+ b1 = alphasur_peppho(1)
+! b1=0.0d0
+ b2 = alphasur_peppho(2)
+ b3 = alphasur_peppho(3)
+ b4 = alphasur_peppho(4)
+ CALL sc_angular
+ sqom1=om1*om1
+ evdwij = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ Fcav=0.0d0
+ eheadtail = 0.0d0
+ dGCLdR=0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ Fcav = 0.0d0
+ dFdR = 0.0d0
+ dCAVdOM1 = 0.0d0
+ dCAVdOM2 = 0.0d0
+ dCAVdOM12 = 0.0d0
+ rij_shift = rij
+ fac = rij_shift**expon
+ c1 = fac * fac * aa_peppho
+! c1 = 0.0d0
+ c2 = fac * bb_peppho
+! c2 = 0.0d0
+ evdwij = c1 + c2
+! Now cavity....................
+ eagle = dsqrt(1.0/rij_shift)
+ top = b1 * ( eagle + b2 * 1.0/rij_shift - b3 )
+ bot = 1.0d0 + b4 * (1.0/rij_shift ** 12.0d0)
+ botsq = bot * bot
+ Fcav = top / bot
+ dtop = b1 * ((1.0/ (2.0d0 * eagle)) + (b2))
+ dbot = 12.0d0 * b4 * (1.0/rij_shift) ** 11.0d0
+ dFdR = ((dtop * bot - top * dbot) / botsq)
+ w1 = wqdip_peppho(1)
+ w2 = wqdip_peppho(2)
+! w1=0.0d0
+! w2=0.0d0
+! pis = sig0head_scbase(itypi,itypj)
+! eps_head = epshead_scbase(itypi,itypj)
+!c!-------------------------------------------------------------------
+
+!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj)
+!c! & +dhead(1,1,itypi,itypj))**2))
+!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj)
+!c! & +dhead(2,1,itypi,itypj))**2))
+
+!c!-------------------------------------------------------------------
+!c! ecl
+ sparrow = w1 * om1
+ hawk = w2 * (1.0d0 - sqom1)
+ Ecl = sparrow * rij_shift**2.0d0 &
+ - hawk * rij_shift**4.0d0
+!c!-------------------------------------------------------------------
+!c! derivative of ecl is Gcl
+!c! dF/dr part
+! rij_shift=5.0
+ dGCLdR = - 2.0d0 * sparrow * rij_shift**3.0d0 &
+ + 4.0d0 * hawk * rij_shift**5.0d0
+!c! dF/dom1
+ dGCLdOM1 = (w1) * (rij_shift**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = (2.0d0 * w2 * om1) * (rij_shift ** 4.0d0)
+ eom1 = dGCLdOM1+dGCLdOM2
+ eom2 = 0.0
+
+ fac = -expon * (c1 + evdwij) * rij_shift+dFdR+dGCLdR
+! fac=0.0
+ gg(1) = fac*xj*rij
+ gg(2) = fac*yj*rij
+ gg(3) = fac*zj*rij
+ do k=1,3
+ gvdwc_peppho(k,j) = gvdwc_peppho(k,j) +gg(k)/2.0
+ gvdwc_peppho(k,j+1) = gvdwc_peppho(k,j+1) +gg(k)/2.0
+ gvdwc_peppho(k,i) = gvdwc_peppho(k,i) -gg(k)/2.0
+ gvdwc_peppho(k,i+1) = gvdwc_peppho(k,i+1) -gg(k)/2.0
+ gg(k)=0.0
+ enddo
+
+ DO k = 1, 3
+ dcosom1(k) = rij* (dc_norm(k,i) - om1 * erij(k))
+ dcosom2(k) = rij* (dc_norm(k,j) - om2 * erij(k))
+ gg(k) = gg(k) + eom1 * dcosom1(k)! + eom2 * dcosom2(k)
+ gvdwc_peppho(k,j)= gvdwc_peppho(k,j) +0.5*( gg(k)) !&
+! - (eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv*2.0
+ gvdwc_peppho(k,j+1)= gvdwc_peppho(k,j+1) +0.5*( gg(k)) !&
+! + (eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv*2.0
+ gvdwc_peppho(k,i)= gvdwc_peppho(k,i) -0.5*( gg(k)) &
+ - (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0
+ gvdwc_peppho(k,i+1)= gvdwc_peppho(k,i+1) - 0.5*( gg(k)) &
+ + (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0
+ enddo
+ epeppho=epeppho+evdwij+Fcav+ECL
+! print *,i,j,evdwij,Fcav,ECL,rij_shift
+ enddo
+ enddo
+ end subroutine eprot_pep_phosphate
+!!!!!!!!!!!!!!!!-------------------------------------------------------------
+ subroutine emomo(evdw)
+ use calc_data
+ use comm_momo
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.LOCAL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.NAMES'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CALC'
+! include 'COMMON.CONTROL'
+! include 'COMMON.SBRIDGE'
+ logical :: lprn
+!el local variables
+ integer :: iint,itypi1,subchap,isel
+ real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,e1,e2,sigm,epsi
+ real(kind=8) :: evdw
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,ssgradlipi,ssgradlipj, &
+ sslipi,sslipj,faclip,alpha_sco
+ integer :: ii
+ real(kind=8) :: fracinbuf
+ real (kind=8) :: escpho
+ real (kind=8),dimension(4):: ener
+ real(kind=8) :: b1,b2,egb
+ real(kind=8) :: Fisocav,ECL,Elj,Equad,Epol,eheadtail,&
+ Lambf,&
+ Chif,ChiLambf,Fcav,dFdR,dFdOM1,&
+ dFdOM2,dFdL,dFdOM12,&
+ federmaus,&
+ d1i,d1j
+! real(kind=8),dimension(3,2)::erhead_tail
+! real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead,Rtail_distance
+ real(kind=8) :: facd4, adler, Fgb, facd3
+ integer troll,jj,istate
+ real (kind=8) :: dcosom1(3),dcosom2(3)
+ eps_out=80.0d0
+ sss_ele_cut=1.0d0
+! print *,"EVDW KURW",evdw,nres
+ do i=iatsc_s,iatsc_e
+! print *,"I am in EVDW",i
+ itypi=iabs(itype(i,1))
+! if (i.ne.47) cycle
+ if (itypi.eq.ntyp1) cycle
+ itypi1=iabs(itype(i+1,1))
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ xi=dmod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=dmod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=dmod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+
+ if ((zi.gt.bordlipbot) &
+ .and.(zi.lt.bordliptop)) then
+!C the energy transfer exist
+ if (zi.lt.buflipbot) then
+!C what fraction I am in
+ fracinbuf=1.0d0- &
+ ((zi-bordlipbot)/lipbufthick)
+!C 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
+! print *, sslipi,ssgradlipi
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+! dsci_inv=dsc_inv(itypi)
+ dsci_inv=vbld_inv(i+nres)
+! write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
+! write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
+!
+! Calculate SC interaction energy.
+!
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+! print *,"JA PIER",i,j,iint,istart(i,iint),iend(i,iint)
+ IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+ call dyn_ssbond_ene(i,j,evdwij)
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
+ 'evdw',i,j,evdwij,' ss'
+! if (energy_dec) write (iout,*) &
+! 'evdw',i,j,evdwij,' ss'
+ do k=j+1,iend(i,iint)
+!C search over all next residues
+ if (dyn_ss_mask(k)) then
+!C check if they are cysteins
+!C write(iout,*) 'k=',k
+
+!c write(iout,*) "PRZED TRI", evdwij
+! evdwij_przed_tri=evdwij
+ call triple_ssbond_ene(i,j,k,evdwij)
+!c if(evdwij_przed_tri.ne.evdwij) then
+!c write (iout,*) "TRI:", evdwij, evdwij_przed_tri
+!c endif
+
+!c write(iout,*) "PO TRI", evdwij
+!C call the energy function that removes the artifical triple disulfide
+!C bond the soubroutine is located in ssMD.F
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
+ 'evdw',i,j,evdwij,'tss'
+ endif!dyn_ss_mask(k)
+ enddo! k
+ ELSE
+!el ind=ind+1
+ itypj=iabs(itype(j,1))
+ if (itypj.eq.ntyp1) cycle
+ CALL elgrad_init(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+
+! if (j.ne.78) cycle
+! dscj_inv=dsc_inv(itypj)
+ dscj_inv=vbld_inv(j+nres)
+ xj=c(1,j+nres)
+ yj=c(2,j+nres)
+ zj=c(3,j+nres)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+ dxj = dc_norm( 1, nres+j )
+ dyj = dc_norm( 2, nres+j )
+ dzj = dc_norm( 3, nres+j )
+! print *,i,j,itypi,itypj
+! d1i=0.0d0
+! d1j=0.0d0
+! BetaT = 1.0d0 / (298.0d0 * Rb)
+! Gay-berne var's
+!1! sig0ij = sigma_scsc( itypi,itypj )
+! chi1=0.0d0
+! chi2=0.0d0
+! chip1=0.0d0
+! chip2=0.0d0
+! not used by momo potential, but needed by sc_angular which is shared
+! by all energy_potential subroutines
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+ a12sq = rborn(itypi,itypj) * rborn(itypj,itypi)
+! a12sq = a12sq * a12sq
+! charge of amino acid itypi is...
+ chis1 = chis(itypi,itypj)
+ chis2 = chis(itypj,itypi)
+ chis12 = chis1 * chis2
+ sig1 = sigmap1(itypi,itypj)
+ sig2 = sigmap2(itypi,itypj)
+! write (*,*) "sig1 = ", sig1
+! chis1=0.0
+! chis2=0.0
+! chis12 = chis1 * chis2
+! sig1=0.0
+! sig2=0.0
+! write (*,*) "sig2 = ", sig2
+! alpha factors from Fcav/Gcav
+ b1cav = alphasur(1,itypi,itypj)
+! b1cav=0.0d0
+ b2cav = alphasur(2,itypi,itypj)
+ b3cav = alphasur(3,itypi,itypj)
+ b4cav = alphasur(4,itypi,itypj)
+! used to determine whether we want to do quadrupole calculations
+ eps_in = epsintab(itypi,itypj)
+ if (eps_in.eq.0.0) eps_in=1.0
+
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+ Rtail = 0.0d0
+! dtail(1,itypi,itypj)=0.0
+! dtail(2,itypi,itypj)=0.0
+
+ DO k = 1, 3
+ ctail(k,1)=c(k,i+nres)-dtail(1,itypi,itypj)*dc_norm(k,nres+i)
+ ctail(k,2)=c(k,j+nres)-dtail(2,itypi,itypj)*dc_norm(k,nres+j)
+ END DO
+!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)))
+
+! write (*,*) "eps_inout_fac = ", eps_inout_fac
+!-------------------------------------------------------------------
+! tail location and distance calculations
+ d1 = dhead(1, 1, itypi, itypj)
+ d2 = dhead(2, 1, itypi, itypj)
+
+ DO k = 1,3
+! location of polar head is computed by taking hydrophobic centre
+! and moving by a d1 * dc_norm vector
+! see unres publications for very informative images
+ chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres)
+ chead(k,2) = c(k, j+nres) + d2 * dc_norm(k, j+nres)
+! distance
+! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ Rhead_distance(k) = chead(k,2) - chead(k,1)
+ END DO
+! pitagoras (root of sum of squares)
+ Rhead = dsqrt( &
+ (Rhead_distance(1)*Rhead_distance(1)) &
+ + (Rhead_distance(2)*Rhead_distance(2)) &
+ + (Rhead_distance(3)*Rhead_distance(3)))
+!-------------------------------------------------------------------
+! zero everything that should be zero'ed
+ evdwij = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ Fcav=0.0d0
+ eheadtail = 0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ Fcav = 0.0d0
+ dFdR = 0.0d0
+ dCAVdOM1 = 0.0d0
+ dCAVdOM2 = 0.0d0
+ dCAVdOM12 = 0.0d0
+ dscj_inv = vbld_inv(j+nres)
+! print *,i,j,dscj_inv,dsci_inv
+! rij holds 1/(distance of Calpha atoms)
+ rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj)
+ rij = dsqrt(rrij)
+!----------------------------
+ CALL sc_angular
+! this should be in elgrad_init but om's are calculated by sc_angular
+! which in turn is used by older potentials
+! om = omega, sqom = om^2
+ sqom1 = om1 * om1
+ sqom2 = om2 * om2
+ sqom12 = om12 * om12
+
+! now we calculate EGB - Gey-Berne
+! It will be summed up in evdwij and saved in evdw
+ sigsq = 1.0D0 / sigsq
+ sig = sig0ij * dsqrt(sigsq)
+! rij_shift = 1.0D0 / rij - sig + sig0ij
+ rij_shift = Rtail - sig + sig0ij
+ IF (rij_shift.le.0.0D0) THEN
+ evdw = 1.0D20
+ RETURN
+ END IF
+ sigder = -sig * sigsq
+ rij_shift = 1.0D0 / rij_shift
+ fac = rij_shift**expon
+ c1 = fac * fac * aa_aq(itypi,itypj)
+! print *,"ADAM",aa_aq(itypi,itypj)
+
+! c1 = 0.0d0
+ c2 = fac * bb_aq(itypi,itypj)
+! c2 = 0.0d0
+ evdwij = eps1 * eps2rt * eps3rt * ( c1 + c2 )
+ eps2der = eps3rt * evdwij
+ eps3der = eps2rt * evdwij
+! evdwij = 4.0d0 * eps2rt * eps3rt * evdwij
+ evdwij = eps2rt * eps3rt * evdwij
+!#ifdef TSCSC
+! IF (bb_aq(itypi,itypj).gt.0) THEN
+! evdw_p = evdw_p + evdwij
+! ELSE
+! evdw_m = evdw_m + evdwij
+! END IF
+!#else
+ evdw = evdw &
+ + evdwij
+!#endif
+
+ c1 = c1 * eps1 * eps2rt**2 * eps3rt**2
+ fac = -expon * (c1 + evdwij) * rij_shift
+ sigder = fac * sigder
+! fac = rij * fac
+! Calculate distance derivative
+ gg(1) = fac
+ gg(2) = fac
+ gg(3) = fac
+! if (b2.gt.0.0) then
+ fac = chis1 * sqom1 + chis2 * sqom2 &
+ - 2.0d0 * chis12 * om1 * om2 * om12
+! we will use pom later in Gcav, so dont mess with it!
+ pom = 1.0d0 - chis1 * chis2 * sqom12
+ Lambf = (1.0d0 - (fac / pom))
+! print *,"fac,pom",fac,pom,Lambf
+ Lambf = dsqrt(Lambf)
+ sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0)
+! print *,"sig1,sig2",sig1,sig2,itypi,itypj
+! write (*,*) "sparrow = ", sparrow
+ Chif = Rtail * sparrow
+! print *,"rij,sparrow",rij , sparrow
+ ChiLambf = Chif * Lambf
+ eagle = dsqrt(ChiLambf)
+ bat = ChiLambf ** 11.0d0
+ top = b1cav * ( eagle + b2cav * ChiLambf - b3cav )
+ bot = 1.0d0 + b4cav * (ChiLambf ** 12.0d0)
+ botsq = bot * bot
+! print *,top,bot,"bot,top",ChiLambf,Chif
+ Fcav = top / bot
+
+ dtop = b1cav * ((Lambf / (2.0d0 * eagle)) + (b2cav * Lambf))
+ dbot = 12.0d0 * b4cav * bat * Lambf
+ dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow
+
+ dtop = b1cav * ((Chif / (2.0d0 * eagle)) + (b2cav * Chif))
+ dbot = 12.0d0 * b4cav * bat * Chif
+ eagle = Lambf * pom
+ dFdOM1 = -(chis1 * om1 - chis12 * om2 * om12) / (eagle)
+ dFdOM2 = -(chis2 * om2 - chis12 * om1 * om12) / (eagle)
+ dFdOM12 = chis12 * (chis1 * om1 * om12 - om2) &
+ * (chis2 * om2 * om12 - om1) / (eagle * pom)
+
+ dFdL = ((dtop * bot - top * dbot) / botsq)
+! dFdL = 0.0d0
+ dCAVdOM1 = dFdL * ( dFdOM1 )
+ dCAVdOM2 = dFdL * ( dFdOM2 )
+ dCAVdOM12 = dFdL * ( dFdOM12 )
+
+ DO k= 1, 3
+ ertail(k) = Rtail_distance(k)/Rtail
+ END DO
+ erdxi = scalar( ertail(1), dC_norm(1,i+nres) )
+ erdxj = scalar( ertail(1), dC_norm(1,j+nres) )
+ facd1 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+ facd2 = dtail(2,itypi,itypj) * vbld_inv(j+nres)
+ DO k = 1, 3
+!c! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
+!c! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
+ pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx(k,i) = gvdwx(k,i) &
+ - (( dFdR + gg(k) ) * pom)
+!c! & - ( dFdR * pom )
+ pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx(k,j) = gvdwx(k,j) &
+ + (( dFdR + gg(k) ) * pom)
+!c! & + ( dFdR * pom )
+
+ gvdwc(k,i) = gvdwc(k,i) &
+ - (( dFdR + gg(k) ) * ertail(k))
+!c! & - ( dFdR * ertail(k))
+
+ gvdwc(k,j) = gvdwc(k,j) &
+ + (( dFdR + gg(k) ) * ertail(k))
+!c! & + ( dFdR * ertail(k))
+
+ gg(k) = 0.0d0
+! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
+! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
+ END DO
+
+
+!c! Compute head-head and head-tail energies for each state
+
+ isel = iabs(Qi) + iabs(Qj)
+! isel=0
+ IF (isel.eq.0) THEN
+!c! No charges - do nothing
+ eheadtail = 0.0d0
+
+ ELSE IF (isel.eq.4) THEN
+!c! Calculate dipole-dipole interactions
+ CALL edd(ecl)
+ eheadtail = ECL
+! eheadtail = 0.0d0
+
+ ELSE IF (isel.eq.1 .and. iabs(Qi).eq.1) THEN
+!c! Charge-nonpolar interactions
+ CALL eqn(epol)
+ eheadtail = epol
+! eheadtail = 0.0d0
+
+ ELSE IF (isel.eq.1 .and. iabs(Qj).eq.1) THEN
+!c! Nonpolar-charge interactions
+ CALL enq(epol)
+ eheadtail = epol
+! eheadtail = 0.0d0
+
+ ELSE IF (isel.eq.3 .and. icharge(itypj).eq.2) THEN
+!c! Charge-dipole interactions
+ CALL eqd(ecl, elj, epol)
+ eheadtail = ECL + elj + epol
+! eheadtail = 0.0d0
+
+ ELSE IF (isel.eq.3 .and. icharge(itypi).eq.2) THEN
+!c! Dipole-charge interactions
+ CALL edq(ecl, elj, epol)
+ eheadtail = ECL + elj + epol
+! eheadtail = 0.0d0
+
+ ELSE IF ((isel.eq.2.and. &
+ iabs(Qi).eq.1).and. &
+ nstate(itypi,itypj).eq.1) THEN
+!c! Same charge-charge interaction ( +/+ or -/- )
+ CALL eqq(Ecl,Egb,Epol,Fisocav,Elj)
+ eheadtail = ECL + Egb + Epol + Fisocav + Elj
+! eheadtail = 0.0d0
+
+ ELSE IF ((isel.eq.2.and. &
+ iabs(Qi).eq.1).and. &
+ nstate(itypi,itypj).ne.1) THEN
+!c! Different charge-charge interaction ( +/- or -/+ )
+ CALL energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad)
+ END IF
+ END IF ! this endif ends the "catch the gly-gly" at the beggining of Fcav
+ evdw = evdw + Fcav + eheadtail
+
+ IF (energy_dec) write (iout,'(2(1x,a3,i3),3f6.2,10f16.7)') &
+ restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
+ 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,&
+ Equad,evdwij+Fcav+eheadtail,evdw
+! evdw = evdw + Fcav + eheadtail
+
+ iF (nstate(itypi,itypj).eq.1) THEN
+ CALL sc_grad
+ END IF
+!c!-------------------------------------------------------------------
+!c! NAPISY KONCOWE
+ END DO ! j
+ END DO ! iint
+ END DO ! i
+!c write (iout,*) "Number of loop steps in EGB:",ind
+!c energy_dec=.false.
+! print *,"EVDW KURW",evdw,nres
+
+ RETURN
+ END SUBROUTINE emomo
+!C------------------------------------------------------------------------------------
+ SUBROUTINE eqq(Ecl,Egb,Epol,Fisocav,Elj)
+ use calc_data
+ use comm_momo
+ real (kind=8) :: facd3, facd4, federmaus, adler,&
+ Ecl,Egb,Epol,Fisocav,Elj,Fgb
+! integer :: k
+!c! Epol and Gpol analytical parameters
+ alphapol1 = alphapol(itypi,itypj)
+ alphapol2 = alphapol(itypj,itypi)
+!c! Fisocav and Gisocav analytical parameters
+ al1 = alphiso(1,itypi,itypj)
+ al2 = alphiso(2,itypi,itypj)
+ al3 = alphiso(3,itypi,itypj)
+ al4 = alphiso(4,itypi,itypj)
+ csig = (1.0d0 &
+ / dsqrt(sigiso1(itypi, itypj)**2.0d0 &
+ + sigiso2(itypi,itypj)**2.0d0))
+!c!
+ pis = sig0head(itypi,itypj)
+ eps_head = epshead(itypi,itypj)
+ Rhead_sq = Rhead * Rhead
+!c! R1 - distance between head of ith side chain and tail of jth sidechain
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R1 = 0.0d0
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances needed by Epol
+ R1=R1+(ctail(k,2)-chead(k,1))**2
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R1 = dsqrt(R1)
+ R2 = dsqrt(R2)
+
+!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj)
+!c! & +dhead(1,1,itypi,itypj))**2))
+!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj)
+!c! & +dhead(2,1,itypi,itypj))**2))
+
+!c!-------------------------------------------------------------------
+!c! Coulomb electrostatic interaction
+ Ecl = (332.0d0 * Qij) / Rhead
+!c! derivative of Ecl is Gcl...
+ dGCLdR = (-332.0d0 * Qij ) / Rhead_sq
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq))
+ Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0)
+ Egb = -(332.0d0 * Qij * eps_inout_fac) / 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 * eps_inout_fac) / (Fgb * 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
+!c! dPOLdR1 = 0.0d0
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
+!c! dPOLdOM2 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Elj
+!c! Lennard-Jones 6-12 interaction between heads
+ pom = (pis / Rhead)**6.0d0
+ Elj = 4.0d0 * eps_head * pom * (pom-1.0d0)
+!c! derivative of Elj is Glj
+ dGLJdR = 4.0d0 * eps_head*(((-12.0d0*pis**12.0d0)/(Rhead**13.0d0))&
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! These things do the dRdX derivatives, that is
+!c! allow us to change what we see from function that changes with
+!c! distance to function that changes with LOCATION (of the interaction
+!c! site)
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1)
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ END DO
+
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j+nres) )
+ bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) )
+ federmaus = scalar(erhead_tail(1,1),dC_norm(1,j+nres))
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
+ facd1 = d1 * vbld_inv(i+nres)
+ facd2 = d2 * vbld_inv(j+nres)
+ facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+ facd4 = dtail(2,itypi,itypj) * vbld_inv(j+nres)
+
+!c! Now we add appropriate partial derivatives (one in each dimension)
+ DO k = 1, 3
+ hawk = (erhead_tail(k,1) + &
+ facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres)))
+ condor = (erhead_tail(k,2) + &
+ facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j+nres)))
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx(k,i) = gvdwx(k,i) &
+ - dGCLdR * pom&
+ - dGGBdR * pom&
+ - dGCVdR * pom&
+ - dPOLdR1 * hawk&
+ - dPOLdR2 * (erhead_tail(k,2)&
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))&
+ - dGLJdR * pom
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx(k,j) = gvdwx(k,j)+ dGCLdR * pom&
+ + dGGBdR * pom+ dGCVdR * pom&
+ + dPOLdR1 * (erhead_tail(k,1)&
+ -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres)))&
+ + dPOLdR2 * condor + dGLJdR * pom
+
+ gvdwc(k,i) = gvdwc(k,i) &
+ - dGCLdR * erhead(k)&
+ - dGGBdR * erhead(k)&
+ - dGCVdR * erhead(k)&
+ - dPOLdR1 * erhead_tail(k,1)&
+ - dPOLdR2 * erhead_tail(k,2)&
+ - dGLJdR * erhead(k)
+
+ gvdwc(k,j) = gvdwc(k,j) &
+ + dGCLdR * erhead(k) &
+ + dGGBdR * erhead(k) &
+ + dGCVdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1) &
+ + dPOLdR2 * erhead_tail(k,2)&
+ + dGLJdR * erhead(k)
+
+ END DO
+ RETURN
+ END SUBROUTINE eqq
+!c!-------------------------------------------------------------------
+ SUBROUTINE energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad)
+ use comm_momo
+ use calc_data
+
+ double precision eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad
+ double precision ener(4)
+ double precision dcosom1(3),dcosom2(3)
+!c! used in Epol derivatives
+ double precision facd3, facd4
+ double precision federmaus, adler
+ integer istate,ii,jj
+ real (kind=8) :: Fgb
+! print *,"CALLING EQUAD"
+!c! Epol and Gpol analytical parameters
+ alphapol1 = alphapol(itypi,itypj)
+ alphapol2 = alphapol(itypj,itypi)
+!c! Fisocav and Gisocav analytical parameters
+ al1 = alphiso(1,itypi,itypj)
+ al2 = alphiso(2,itypi,itypj)
+ al3 = alphiso(3,itypi,itypj)
+ al4 = alphiso(4,itypi,itypj)
+ csig = (1.0d0 / dsqrt(sigiso1(itypi, itypj)**2.0d0&
+ + sigiso2(itypi,itypj)**2.0d0))
+!c!
+ w1 = wqdip(1,itypi,itypj)
+ w2 = wqdip(2,itypi,itypj)
+ pis = sig0head(itypi,itypj)
+ eps_head = epshead(itypi,itypj)
+!c! First things first:
+!c! We need to do sc_grad's job with GB and Fcav
+ eom1 = eps2der * eps2rt_om1 &
+ - 2.0D0 * alf1 * eps3der&
+ + sigder * sigsq_om1&
+ + dCAVdOM1
+ eom2 = eps2der * eps2rt_om2 &
+ + 2.0D0 * alf2 * eps3der&
+ + sigder * sigsq_om2&
+ + dCAVdOM2
+ eom12 = evdwij * eps1_om12 &
+ + eps2der * eps2rt_om12 &
+ - 2.0D0 * alf12 * eps3der&
+ + sigder *sigsq_om12&
+ + dCAVdOM12
+!c! now some magical transformations to project gradient into
+!c! three cartesian vectors
+ DO k = 1, 3
+ dcosom1(k) = rij * (dc_norm(k,nres+i) - om1 * erij(k))
+ dcosom2(k) = rij * (dc_norm(k,nres+j) - om2 * erij(k))
+ gg(k) = gg(k) + eom1 * dcosom1(k) + eom2 * dcosom2(k)
+!c! this acts on hydrophobic center of interaction
+ gvdwx(k,i)= gvdwx(k,i) - gg(k) &
+ + (eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))&
+ + eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+ gvdwx(k,j)= gvdwx(k,j) + gg(k) &
+ + (eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))&
+ + eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+!c! this acts on Calpha
+ gvdwc(k,i)=gvdwc(k,i)-gg(k)
+ gvdwc(k,j)=gvdwc(k,j)+gg(k)
+ END DO
+!c! sc_grad is done, now we will compute
+ eheadtail = 0.0d0
+ eom1 = 0.0d0
+ eom2 = 0.0d0
+ eom12 = 0.0d0
+ DO istate = 1, nstate(itypi,itypj)
+!c*************************************************************
+ IF (istate.ne.1) THEN
+ IF (istate.lt.3) THEN
+ ii = 1
+ ELSE
+ ii = 2
+ END IF
+ jj = istate/ii
+ d1 = dhead(1,ii,itypi,itypj)
+ d2 = dhead(2,jj,itypi,itypj)
+ DO k = 1,3
+ chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres)
+ chead(k,2) = c(k, j+nres) + d2 * dc_norm(k, j+nres)
+ Rhead_distance(k) = chead(k,2) - chead(k,1)
+ END DO
+!c! pitagoras (root of sum of squares)
+ Rhead = dsqrt( &
+ (Rhead_distance(1)*Rhead_distance(1)) &
+ + (Rhead_distance(2)*Rhead_distance(2)) &
+ + (Rhead_distance(3)*Rhead_distance(3)))
+ END IF
+ Rhead_sq = Rhead * Rhead
+
+!c! R1 - distance between head of ith side chain and tail of jth sidechain
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R1 = 0.0d0
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R1=R1+(ctail(k,2)-chead(k,1))**2
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R1 = dsqrt(R1)
+ R2 = dsqrt(R2)
+ Ecl = (332.0d0 * Qij) / (Rhead * eps_in)
+!c! Ecl = 0.0d0
+!c! write (*,*) "Ecl = ", Ecl
+!c! derivative of Ecl is Gcl...
+ dGCLdR = (-332.0d0 * Qij ) / (Rhead_sq * eps_in)
+!c! dGCLdR = 0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Generalised Born Solvent Polarization
+ ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq))
+ Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0)
+ Egb = -(332.0d0 * Qij * eps_inout_fac) / Fgb
+!c! Egb = 0.0d0
+!c! write (*,*) "a1*a2 = ", a12sq
+!c! write (*,*) "Rhead = ", Rhead
+!c! write (*,*) "Rhead_sq = ", Rhead_sq
+!c! write (*,*) "ee = ", ee
+!c! write (*,*) "Fgb = ", Fgb
+!c! write (*,*) "fac = ", eps_inout_fac
+!c! write (*,*) "Qij = ", Qij
+!c! write (*,*) "Egb = ", Egb
+!c! Derivative of Egb is Ggb...
+!c! dFGBdR is used by Quad's later...
+ dGGBdFGB = -(-332.0d0 * Qij * eps_inout_fac) / (Fgb * Fgb)
+ dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )&
+ / ( 2.0d0 * Fgb )
+ dGGBdR = dGGBdFGB * dFGBdR
+!c! dGGBdR = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Fisocav - isotropic cavity creation term
+ pom = Rhead * csig
+ top = al1 * (dsqrt(pom) + al2 * pom - al3)
+ bot = (1.0d0 + al4 * pom**12.0d0)
+ botsq = bot * bot
+ FisoCav = top / bot
+ dtop = al1 * ((1.0d0 / (2.0d0 * dsqrt(pom))) + al2)
+ dbot = 12.0d0 * al4 * pom ** 11.0d0
+ dGCVdR = ((dtop * bot - top * dbot) / botsq) * csig
+!c! dGCVdR = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Polarization energy
+!c! Epol
+ MomoFac1 = (1.0d0 - chi1 * sqom2)
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR1 = ( R1 * R1 ) / MomoFac1
+ RR2 = ( R2 * R2 ) / MomoFac2
+ ee1 = exp(-( RR1 / (4.0d0 * a12sq) ))
+ ee2 = exp(-( RR2 / (4.0d0 * a12sq) ))
+ fgb1 = sqrt( RR1 + a12sq * ee1 )
+ fgb2 = sqrt( RR2 + a12sq * ee2 )
+ epol = 332.0d0 * eps_inout_fac * (&
+ (( alphapol1 / fgb1 )**4.0d0)+((alphapol2/fgb2) ** 4.0d0 ))
+!c! epol = 0.0d0
+!c! derivative of Epol is Gpol...
+ dPOLdFGB1 = -(1328.0d0 * eps_inout_fac * alphapol1 ** 4.0d0)&
+ / (fgb1 ** 5.0d0)
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0)&
+ / (fgb2 ** 5.0d0)
+ dFGBdR1 = ( (R1 / MomoFac1) &
+ * ( 2.0d0 - (0.5d0 * ee1) ) )&
+ / ( 2.0d0 * fgb1 )
+ dFGBdR2 = ( (R2 / MomoFac2) &
+ * ( 2.0d0 - (0.5d0 * ee2) ) ) &
+ / ( 2.0d0 * fgb2 )
+ dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1)) &
+ * ( 2.0d0 - 0.5d0 * ee1) ) &
+ / ( 2.0d0 * fgb1 )
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) &
+ * ( 2.0d0 - 0.5d0 * ee2) ) &
+ / ( 2.0d0 * fgb2 )
+ dPOLdR1 = dPOLdFGB1 * dFGBdR1
+!c! dPOLdR1 = 0.0d0
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
+ pom = (pis / Rhead)**6.0d0
+ Elj = 4.0d0 * eps_head * pom * (pom-1.0d0)
+!c! Elj = 0.0d0
+!c! derivative of Elj is Glj
+ dGLJdR = 4.0d0 * eps_head &
+ * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) &
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))
+!c! dGLJdR = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Equad
+ IF (Wqd.ne.0.0d0) THEN
+ Beta1 = 5.0d0 + 3.0d0 * (sqom12 - 1.0d0) &
+ - 37.5d0 * ( sqom1 + sqom2 ) &
+ + 157.5d0 * ( sqom1 * sqom2 ) &
+ - 45.0d0 * om1*om2*om12
+ fac = -( Wqd / (2.0d0 * Fgb**5.0d0) )
+ Equad = fac * Beta1
+!c! Equad = 0.0d0
+!c! derivative of Equad...
+ dQUADdR = ((2.5d0 * Wqd * Beta1) / (Fgb**6.0d0)) * dFGBdR
+!c! dQUADdR = 0.0d0
+ dQUADdOM1 = fac* (-75.0d0*om1 + 315.0d0*om1*sqom2 - 45.0d0*om2*om12)
+!c! dQUADdOM1 = 0.0d0
+ dQUADdOM2 = fac* (-75.0d0*om2 + 315.0d0*sqom1*om2 - 45.0d0*om1*om12)
+!c! dQUADdOM2 = 0.0d0
+ dQUADdOM12 = fac * ( 6.0d0*om12 - 45.0d0*om1*om2 )
+ ELSE
+ Beta1 = 0.0d0
+ Equad = 0.0d0
+ END IF
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! Angular stuff
+ eom1 = dPOLdOM1 + dQUADdOM1
+ eom2 = dPOLdOM2 + dQUADdOM2
+ eom12 = dQUADdOM12
+!c! now some magical transformations to project gradient into
+!c! three cartesian vectors
+ DO k = 1, 3
+ dcosom1(k) = rij * (dc_norm(k,nres+i) - om1 * erij(k))
+ dcosom2(k) = rij * (dc_norm(k,nres+j) - om2 * erij(k))
+ tuna(k) = eom1 * dcosom1(k) + eom2 * dcosom2(k)
+ END DO
+!c! Radial stuff
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1)
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ END DO
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j+nres) )
+ bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) )
+ federmaus = scalar(erhead_tail(1,1),dC_norm(1,j+nres))
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
+ facd1 = d1 * vbld_inv(i+nres)
+ facd2 = d2 * vbld_inv(j+nres)
+ facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+ facd4 = dtail(2,itypi,itypj) * vbld_inv(j+nres)
+ DO k = 1, 3
+ hawk = erhead_tail(k,1) + &
+ facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres))
+ condor = erhead_tail(k,2) + &
+ facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j+nres))
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+!c! this acts on hydrophobic center of interaction
+ gheadtail(k,1,1) = gheadtail(k,1,1) &
+ - dGCLdR * pom &
+ - dGGBdR * pom &
+ - dGCVdR * pom &
+ - dPOLdR1 * hawk &
+ - dPOLdR2 * (erhead_tail(k,2) &
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))&
+ - dGLJdR * pom &
+ - dQUADdR * pom&
+ - tuna(k) &
+ + (eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))&
+ + eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+!c! this acts on hydrophobic center of interaction
+ gheadtail(k,2,1) = gheadtail(k,2,1) &
+ + dGCLdR * pom &
+ + dGGBdR * pom &
+ + dGCVdR * pom &
+ + dPOLdR1 * (erhead_tail(k,1) &
+ -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) &
+ + dPOLdR2 * condor &
+ + dGLJdR * pom &
+ + dQUADdR * pom &
+ + tuna(k) &
+ + (eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
+ + eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+
+!c! this acts on Calpha
+ gheadtail(k,3,1) = gheadtail(k,3,1) &
+ - dGCLdR * erhead(k)&
+ - dGGBdR * erhead(k)&
+ - dGCVdR * erhead(k)&
+ - dPOLdR1 * erhead_tail(k,1)&
+ - dPOLdR2 * erhead_tail(k,2)&
+ - dGLJdR * erhead(k) &
+ - dQUADdR * erhead(k)&
+ - tuna(k)
+!c! this acts on Calpha
+ gheadtail(k,4,1) = gheadtail(k,4,1) &
+ + dGCLdR * erhead(k) &
+ + dGGBdR * erhead(k) &
+ + dGCVdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1) &
+ + dPOLdR2 * erhead_tail(k,2) &
+ + dGLJdR * erhead(k) &
+ + dQUADdR * erhead(k)&
+ + tuna(k)
+ END DO
+ ener(istate) = ECL + Egb + Epol + Fisocav + Elj + Equad
+ eheadtail = eheadtail &
+ + wstate(istate, itypi, itypj) &
+ * dexp(-betaT * ener(istate))
+!c! foreach cartesian dimension
+ DO k = 1, 3
+!c! foreach of two gvdwx and gvdwc
+ DO l = 1, 4
+ gheadtail(k,l,2) = gheadtail(k,l,2) &
+ + wstate( istate, itypi, itypj ) &
+ * dexp(-betaT * ener(istate)) &
+ * gheadtail(k,l,1)
+ gheadtail(k,l,1) = 0.0d0
+ END DO
+ END DO
+ END DO
+!c! Here ended the gigantic DO istate = 1, 4, which starts
+!c! at the beggining of the subroutine
+
+ DO k = 1, 3
+ DO l = 1, 4
+ gheadtail(k,l,2) = gheadtail(k,l,2) / eheadtail
+ END DO
+ gvdwx(k,i) = gvdwx(k,i) + gheadtail(k,1,2)
+ gvdwx(k,j) = gvdwx(k,j) + gheadtail(k,2,2)
+ gvdwc(k,i) = gvdwc(k,i) + gheadtail(k,3,2)
+ gvdwc(k,j) = gvdwc(k,j) + gheadtail(k,4,2)
+ DO l = 1, 4
+ gheadtail(k,l,1) = 0.0d0
+ gheadtail(k,l,2) = 0.0d0
+ END DO
+ END DO
+ eheadtail = (-dlog(eheadtail)) / betaT
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ dQUADdOM1 = 0.0d0
+ dQUADdOM2 = 0.0d0
+ dQUADdOM12 = 0.0d0
+ RETURN
+ END SUBROUTINE energy_quad
+!!-----------------------------------------------------------
+ SUBROUTINE eqn(Epol)
+ use comm_momo
+ use calc_data
+
+ double precision facd4, federmaus,epol
+ alphapol1 = alphapol(itypi,itypj)
+!c! R1 - distance between head of ith side chain and tail of jth sidechain
+ R1 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R1=R1+(ctail(k,2)-chead(k,1))**2
+ END DO
+!c! Pitagoras
+ R1 = dsqrt(R1)
+
+!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj)
+!c! & +dhead(1,1,itypi,itypj))**2))
+!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj)
+!c! & +dhead(2,1,itypi,itypj))**2))
+!c--------------------------------------------------------------------
+!c Polarization energy
+!c Epol
+ MomoFac1 = (1.0d0 - chi1 * sqom2)
+ RR1 = R1 * R1 / MomoFac1
+ ee1 = exp(-( RR1 / (4.0d0 * a12sq) ))
+ fgb1 = sqrt( RR1 + a12sq * ee1)
+ epol = 332.0d0 * eps_inout_fac * (( alphapol1 / fgb1 )**4.0d0)
+ dPOLdFGB1 = -(1328.0d0 * eps_inout_fac * alphapol1 ** 4.0d0) &
+ / (fgb1 ** 5.0d0)
+ dFGBdR1 = ( (R1 / MomoFac1) &
+ * ( 2.0d0 - (0.5d0 * ee1) ) ) &
+ / ( 2.0d0 * fgb1 )
+ dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1)) &
+ * (2.0d0 - 0.5d0 * ee1) ) &
+ / (2.0d0 * fgb1)
+ dPOLdR1 = dPOLdFGB1 * dFGBdR1
+!c! dPOLdR1 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
+ DO k = 1, 3
+ erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1)
+ END DO
+ bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) )
+ federmaus = scalar(erhead_tail(1,1),dC_norm(1,j+nres))
+ facd1 = d1 * vbld_inv(i+nres)
+ facd4 = dtail(2,itypi,itypj) * vbld_inv(j+nres)
+
+ DO k = 1, 3
+ hawk = (erhead_tail(k,1) + &
+ facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres)))
+
+ gvdwx(k,i) = gvdwx(k,i) &
+ - dPOLdR1 * hawk
+ gvdwx(k,j) = gvdwx(k,j) &
+ + dPOLdR1 * (erhead_tail(k,1) &
+ -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres)))
+
+ gvdwc(k,i) = gvdwc(k,i) - dPOLdR1 * erhead_tail(k,1)
+ gvdwc(k,j) = gvdwc(k,j) + dPOLdR1 * erhead_tail(k,1)
+
+ END DO
+ RETURN
+ END SUBROUTINE eqn
+ SUBROUTINE enq(Epol)
+ use calc_data
+ use comm_momo
+ double precision facd3, adler,epol
+ alphapol2 = alphapol(itypj,itypi)
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R2 = dsqrt(R2)
+
+!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj)
+!c! & +dhead(1,1,itypi,itypj))**2))
+!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj)
+!c! & +dhead(2,1,itypi,itypj))**2))
+!c------------------------------------------------------------------------
+!c Polarization energy
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR2 = R2 * R2 / MomoFac2
+ ee2 = exp(-(RR2 / (4.0d0 * a12sq)))
+ fgb2 = sqrt(RR2 + a12sq * ee2)
+ epol = 332.0d0 * eps_inout_fac * ((alphapol2/fgb2) ** 4.0d0 )
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0) &
+ / (fgb2 ** 5.0d0)
+ dFGBdR2 = ( (R2 / MomoFac2) &
+ * ( 2.0d0 - (0.5d0 * ee2) ) ) &
+ / (2.0d0 * fgb2)
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) &
+ * (2.0d0 - 0.5d0 * ee2) ) &
+ / (2.0d0 * fgb2)
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! (See comments in Eqq)
+ DO k = 1, 3
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ END DO
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
+ facd2 = d2 * vbld_inv(j+nres)
+ facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+ DO k = 1, 3
+ condor = (erhead_tail(k,2) &
+ + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j+nres)))
+
+ gvdwx(k,i) = gvdwx(k,i) &
+ - dPOLdR2 * (erhead_tail(k,2) &
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))
+ gvdwx(k,j) = gvdwx(k,j) &
+ + dPOLdR2 * condor
+
+ gvdwc(k,i) = gvdwc(k,i) &
+ - dPOLdR2 * erhead_tail(k,2)
+ gvdwc(k,j) = gvdwc(k,j) &
+ + dPOLdR2 * erhead_tail(k,2)
+
+ END DO
+ RETURN
+ END SUBROUTINE enq
+ SUBROUTINE eqd(Ecl,Elj,Epol)
+ use calc_data
+ use comm_momo
+ double precision facd4, federmaus,ecl,elj,epol
+ alphapol1 = alphapol(itypi,itypj)
+ w1 = wqdip(1,itypi,itypj)
+ w2 = wqdip(2,itypi,itypj)
+ pis = sig0head(itypi,itypj)
+ eps_head = epshead(itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! R1 - distance between head of ith side chain and tail of jth sidechain
+ R1 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R1=R1+(ctail(k,2)-chead(k,1))**2
+ END DO
+!c! Pitagoras
+ R1 = dsqrt(R1)
+
+!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj)
+!c! & +dhead(1,1,itypi,itypj))**2))
+!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj)
+!c! & +dhead(2,1,itypi,itypj))**2))
+
+!c!-------------------------------------------------------------------
+!c! ecl
+ sparrow = w1 * Qi * om1
+ hawk = w2 * Qi * Qi * (1.0d0 - sqom2)
+ Ecl = sparrow / Rhead**2.0d0 &
+ - hawk / Rhead**4.0d0
+ dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 &
+ + 4.0d0 * hawk / Rhead**5.0d0
+!c! dF/dom1
+ dGCLdOM1 = (w1 * Qi) / (Rhead**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = (2.0d0 * w2 * Qi * Qi * om2) / (Rhead ** 4.0d0)
+!c--------------------------------------------------------------------
+!c Polarization energy
+!c Epol
+ MomoFac1 = (1.0d0 - chi1 * sqom2)
+ RR1 = R1 * R1 / MomoFac1
+ ee1 = exp(-( RR1 / (4.0d0 * a12sq) ))
+ fgb1 = sqrt( RR1 + a12sq * ee1)
+ epol = 332.0d0 * eps_inout_fac * (( alphapol1 / fgb1 )**4.0d0)
+!c! epol = 0.0d0
+!c!------------------------------------------------------------------
+!c! derivative of Epol is Gpol...
+ dPOLdFGB1 = -(1328.0d0 * eps_inout_fac * alphapol1 ** 4.0d0) &
+ / (fgb1 ** 5.0d0)
+ dFGBdR1 = ( (R1 / MomoFac1) &
+ * ( 2.0d0 - (0.5d0 * ee1) ) ) &
+ / ( 2.0d0 * fgb1 )
+ dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1)) &
+ * (2.0d0 - 0.5d0 * ee1) ) &
+ / (2.0d0 * fgb1)
+ dPOLdR1 = dPOLdFGB1 * dFGBdR1
+!c! dPOLdR1 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
+!c! dPOLdOM2 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Elj
+ pom = (pis / Rhead)**6.0d0
+ Elj = 4.0d0 * eps_head * pom * (pom-1.0d0)
+!c! derivative of Elj is Glj
+ dGLJdR = 4.0d0 * eps_head &
+ * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) &
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1)
+ END DO
+
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j+nres) )
+ bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) )
+ federmaus = scalar(erhead_tail(1,1),dC_norm(1,j+nres))
+ facd1 = d1 * vbld_inv(i+nres)
+ facd2 = d2 * vbld_inv(j+nres)
+ facd4 = dtail(2,itypi,itypj) * vbld_inv(j+nres)
+
+ DO k = 1, 3
+ hawk = (erhead_tail(k,1) + &
+ facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres)))
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx(k,i) = gvdwx(k,i) &
+ - dGCLdR * pom&
+ - dPOLdR1 * hawk &
+ - dGLJdR * pom
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx(k,j) = gvdwx(k,j) &
+ + dGCLdR * pom &
+ + dPOLdR1 * (erhead_tail(k,1) &
+ -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres))) &
+ + dGLJdR * pom
+
+
+ gvdwc(k,i) = gvdwc(k,i) &
+ - dGCLdR * erhead(k) &
+ - dPOLdR1 * erhead_tail(k,1) &
+ - dGLJdR * erhead(k)
+
+ gvdwc(k,j) = gvdwc(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1) &
+ + dGLJdR * erhead(k)
+
+ END DO
+ RETURN
+ END SUBROUTINE eqd
+ SUBROUTINE edq(Ecl,Elj,Epol)
+! IMPLICIT NONE
+ use comm_momo
+ use calc_data
+
+ double precision facd3, adler,ecl,elj,epol
+ alphapol2 = alphapol(itypj,itypi)
+ w1 = wqdip(1,itypi,itypj)
+ w2 = wqdip(2,itypi,itypj)
+ pis = sig0head(itypi,itypj)
+ eps_head = epshead(itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R2 = dsqrt(R2)
+
+!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj)
+!c! & +dhead(1,1,itypi,itypj))**2))
+!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj)
+!c! & +dhead(2,1,itypi,itypj))**2))
+
+
+!c!-------------------------------------------------------------------
+!c! ecl
+ sparrow = w1 * Qi * om1
+ hawk = w2 * Qi * Qi * (1.0d0 - sqom2)
+ ECL = sparrow / Rhead**2.0d0 &
+ - hawk / Rhead**4.0d0
+!c!-------------------------------------------------------------------
+!c! derivative of ecl is Gcl
+!c! dF/dr part
+ dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 &
+ + 4.0d0 * hawk / Rhead**5.0d0
+!c! dF/dom1
+ dGCLdOM1 = (w1 * Qi) / (Rhead**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = (2.0d0 * w2 * Qi * Qi * om2) / (Rhead ** 4.0d0)
+!c--------------------------------------------------------------------
+!c Polarization energy
+!c Epol
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR2 = R2 * R2 / MomoFac2
+ ee2 = exp(-(RR2 / (4.0d0 * a12sq)))
+ fgb2 = sqrt(RR2 + a12sq * ee2)
+ epol = 332.0d0 * eps_inout_fac * ((alphapol2/fgb2) ** 4.0d0 )
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0) &
+ / (fgb2 ** 5.0d0)
+ dFGBdR2 = ( (R2 / MomoFac2) &
+ * ( 2.0d0 - (0.5d0 * ee2) ) ) &
+ / (2.0d0 * fgb2)
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) &
+ * (2.0d0 - 0.5d0 * ee2) ) &
+ / (2.0d0 * fgb2)
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2
+!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)))
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! (see comments in Eqq)
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ END DO
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j+nres) )
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
+ facd1 = d1 * vbld_inv(i+nres)
+ facd2 = d2 * vbld_inv(j+nres)
+ facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+ DO k = 1, 3
+ condor = (erhead_tail(k,2) &
+ + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j+nres)))
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx(k,i) = gvdwx(k,i) &
+ - dGCLdR * pom &
+ - dPOLdR2 * (erhead_tail(k,2) &
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) &
+ - dGLJdR * pom
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx(k,j) = gvdwx(k,j) &
+ + dGCLdR * pom &
+ + dPOLdR2 * condor &
+ + dGLJdR * pom
+
+
+ gvdwc(k,i) = gvdwc(k,i) &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k)
+
+ gvdwc(k,j) = gvdwc(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR2 * erhead_tail(k,2) &
+ + dGLJdR * erhead(k)
+
+ END DO
+ RETURN
+ END SUBROUTINE edq
+ SUBROUTINE edd(ECL)
+! IMPLICIT NONE
+ use comm_momo
+ use calc_data
+
+ double precision ecl
+!c! csig = sigiso(itypi,itypj)
+ w1 = wqdip(1,itypi,itypj)
+ w2 = wqdip(2,itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! ECL
+ fac = (om12 - 3.0d0 * om1 * om2)
+ c1 = (w1 / (Rhead**3.0d0)) * fac
+ c2 = (w2 / Rhead ** 6.0d0) &
+ * (4.0d0 + fac * fac -3.0d0 * (sqom1 + sqom2))
+ ECL = c1 - c2
+!c! write (*,*) "w1 = ", w1
+!c! write (*,*) "w2 = ", w2
+!c! write (*,*) "om1 = ", om1
+!c! write (*,*) "om2 = ", om2
+!c! write (*,*) "om12 = ", om12
+!c! write (*,*) "fac = ", fac
+!c! write (*,*) "c1 = ", c1
+!c! write (*,*) "c2 = ", c2
+!c! write (*,*) "Ecl = ", Ecl
+!c! write (*,*) "c2_1 = ", (w2 / Rhead ** 6.0d0)
+!c! write (*,*) "c2_2 = ",
+!c! & (4.0d0 + fac * fac -3.0d0 * (sqom1 + sqom2))
+!c!-------------------------------------------------------------------
+!c! dervative of ECL is GCL...
+!c! dECL/dr
+ c1 = (-3.0d0 * w1 * fac) / (Rhead ** 4.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead ** 7.0d0) &
+ * (4.0d0 + fac * fac - 3.0d0 * (sqom1 + sqom2))
+ dGCLdR = c1 - c2
+!c! dECL/dom1
+ c1 = (-3.0d0 * w1 * om2 ) / (Rhead**3.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead**6.0d0) &
+ * ( om2 * om12 - 3.0d0 * om1 * sqom2 + om1 )
+ dGCLdOM1 = c1 - c2
+!c! dECL/dom2
+ c1 = (-3.0d0 * w1 * om1 ) / (Rhead**3.0d0)
+ c2 = (-6.0d0 * w2) / (Rhead**6.0d0) &
+ * ( om1 * om12 - 3.0d0 * sqom1 * om2 + om2 )
+ dGCLdOM2 = c1 - c2
+!c! dECL/dom12
+ c1 = w1 / (Rhead ** 3.0d0)
+ c2 = ( 2.0d0 * w2 * fac ) / Rhead ** 6.0d0
+ dGCLdOM12 = c1 - c2
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! (see comments in Eqq)
+ DO k= 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ END DO
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j+nres) )
+ facd1 = d1 * vbld_inv(i+nres)
+ facd2 = d2 * vbld_inv(j+nres)
+ DO k = 1, 3
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx(k,i) = gvdwx(k,i) - dGCLdR * pom
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j+nres))
+ gvdwx(k,j) = gvdwx(k,j) + dGCLdR * pom
+
+ gvdwc(k,i) = gvdwc(k,i) - dGCLdR * erhead(k)
+ gvdwc(k,j) = gvdwc(k,j) + dGCLdR * erhead(k)
+ END DO
+ RETURN
+ END SUBROUTINE edd
+ SUBROUTINE elgrad_init(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+! IMPLICIT NONE
+ use comm_momo
+ use calc_data
+
+ real(kind=8) :: eheadtail,Egb,Ecl,Elj,Equad,Epol,Rb
+ eps_out=80.0d0
+ itypi = itype(i,1)
+ itypj = itype(j,1)
+!c! 1/(Gas Constant * Thermostate temperature) = BetaT
+!c! ENABLE THIS LINE WHEN USING CHECKGRAD!!!
+!c! t_bath = 300
+!c! BetaT = 1.0d0 / (t_bath * Rb)i
+ Rb=0.001986d0
+ BetaT = 1.0d0 / (298.0d0 * Rb)
+!c! Gay-berne var's
+ sig0ij = sigma( itypi,itypj )
+ chi1 = chi( itypi, itypj )
+ chi2 = chi( itypj, itypi )
+ chi12 = chi1 * chi2
+ chip1 = chipp( itypi, itypj )
+ chip2 = chipp( itypj, itypi )
+ chip12 = chip1 * chip2
+! chi1=0.0
+! chi2=0.0
+! chi12=0.0
+! chip1=0.0
+! chip2=0.0
+! chip12=0.0
+!c! not used by momo potential, but needed by sc_angular which is shared
+!c! by all energy_potential subroutines
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+!c! location, location, location
+! xj = c( 1, nres+j ) - xi
+! yj = c( 2, nres+j ) - yi
+! zj = c( 3, nres+j ) - zi
+ dxj = dc_norm( 1, nres+j )
+ dyj = dc_norm( 2, nres+j )
+ dzj = dc_norm( 3, nres+j )
+!c! distance from center of chain(?) to polar/charged head
+!c! write (*,*) "istate = ", 1
+!c! write (*,*) "ii = ", 1
+!c! write (*,*) "jj = ", 1
+ d1 = dhead(1, 1, itypi, itypj)
+ d2 = dhead(2, 1, itypi, itypj)
+!c! ai*aj from Fgb
+ a12sq = rborn(itypi,itypj) * rborn(itypj,itypi)
+!c! a12sq = a12sq * a12sq
+!c! charge of amino acid itypi is...
+ Qi = icharge(itypi)
+ Qj = icharge(itypj)
+ Qij = Qi * Qj
+!c! chis1,2,12
+ chis1 = chis(itypi,itypj)
+ chis2 = chis(itypj,itypi)
+ chis12 = chis1 * chis2
+ sig1 = sigmap1(itypi,itypj)
+ sig2 = sigmap2(itypi,itypj)
+!c! write (*,*) "sig1 = ", sig1
+!c! write (*,*) "sig2 = ", sig2
+!c! alpha factors from Fcav/Gcav
+ b1cav = alphasur(1,itypi,itypj)
+! b1cav=0.0
+ b2cav = alphasur(2,itypi,itypj)
+ b3cav = alphasur(3,itypi,itypj)
+ b4cav = alphasur(4,itypi,itypj)
+ wqd = wquad(itypi, itypj)
+!c! used by Fgb
+ eps_in = epsintab(itypi,itypj)
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+!c! write (*,*) "eps_inout_fac = ", eps_inout_fac
+!c!-------------------------------------------------------------------
+!c! tail location and distance calculations
+ Rtail = 0.0d0
+ DO k = 1, 3
+ ctail(k,1)=c(k,i+nres)-dtail(1,itypi,itypj)*dc_norm(k,nres+i)
+ ctail(k,2)=c(k,j+nres)-dtail(2,itypi,itypj)*dc_norm(k,nres+j)
+ END DO
+!c! tail distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+ Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 )
+ Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 )
+ Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 )
+ Rtail = dsqrt( &
+ (Rtail_distance(1)*Rtail_distance(1)) &
+ + (Rtail_distance(2)*Rtail_distance(2)) &
+ + (Rtail_distance(3)*Rtail_distance(3)))
+!c!-------------------------------------------------------------------
+!c! Calculate location and distance between polar heads
+!c! distance between heads
+!c! for each one of our three dimensional space...
+ d1 = dhead(1, 1, itypi, itypj)
+ d2 = dhead(2, 1, itypi, itypj)
+
+ DO k = 1,3
+!c! location of polar head is computed by taking hydrophobic centre
+!c! and moving by a d1 * dc_norm vector
+!c! see unres publications for very informative images
+ chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres)
+ chead(k,2) = c(k, j+nres) + d2 * dc_norm(k, j+nres)
+!c! distance
+!c! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+!c! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ Rhead_distance(k) = chead(k,2) - chead(k,1)
+ END DO
+!c! pitagoras (root of sum of squares)
+ Rhead = dsqrt( &
+ (Rhead_distance(1)*Rhead_distance(1)) &
+ + (Rhead_distance(2)*Rhead_distance(2)) &
+ + (Rhead_distance(3)*Rhead_distance(3)))
+!c!-------------------------------------------------------------------
+!c! zero everything that should be zero'ed
+ Egb = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ eheadtail = 0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ RETURN
+ END SUBROUTINE elgrad_init
projects / unres4.git / blobdiff