subroutine etotal(energia)
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
- use MD_data, only: totT
+ use MD_data
#ifndef ISNAN
external proc_proc
#ifdef WINPGI
#ifdef MPI
real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw
+! shielding effect varibles for MPI
+! real(kind=8) fac_shieldbuf(maxres),
+! & grad_shield_locbuf(3,maxcontsshi,-1:maxres),
+! & grad_shield_sidebuf(3,maxcontsshi,-1:maxres),
+! & grad_shieldbuf(3,-1:maxres)
+! integer ishield_listbuf(maxres),
+! &shield_listbuf(maxcontsshi,maxres)
+
! print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
! & " nfgtasks",nfgtasks
if (nfgtasks.gt.1) then
! include 'COMMON.SBRIDGE'
logical :: lprn
!el local variables
- integer :: iint,itypi,itypi1,itypj
+ 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
integer :: ii
!cccc energy_dec=.false.
! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
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)
! alf1=0.0D0
! alf2=0.0D0
! alf12=0.0D0
- 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)
! write(iout,*)"c ", c(1,:), c(2,:), c(3,:)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
+ sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+! print *,sss_ele_cut,sss_ele_grad,&
+! 1.0d0/(rij),r_cut_ele,rlamb_ele
+ if (sss_ele_cut.le.0.0) cycle
! Calculate angle-dependent terms of energy and contributions to their
! derivatives.
call sc_angular
! write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,& !d
! " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2," fac",fac !d
evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij
+ evdw=evdw+evdwij*sss_ele_cut
if (lprn) then
sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac
+! print *,'before fac',fac,rij,evdwij
+ fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
+ /sigma(itypi,itypj)*rij
+! print *,'grad part scale',fac, &
+! evdwij*sss_ele_grad/sss_ele_cut &
+! /sigma(itypi,itypj)*rij
! fac=0.0d0
! Calculate the radial part of the gradient
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
+! print *,'before sc_grad', gg(1),gg(2),gg(3)
! Calculate angular part of the gradient.
call sc_grad
ENDIF ! dyn_ss
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
num_conti=0
call eelecij(i,i+2,ees,evdw1,eel_loc)
if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
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
num_conti=num_cont_hb(i)
call eelecij(i,i+3,ees,evdw1,eel_loc)
if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) &
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
+
! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
! include 'COMMON.VECTORS'
! include 'COMMON.FFIELD'
! include 'COMMON.TIME1'
- real(kind=8),dimension(3) :: ggg,gggp,gggm,erij,dcosb,dcosg
+ real(kind=8),dimension(3) :: ggg,gggp,gggm,erij,dcosb,dcosg,xtemp
real(kind=8),dimension(3,3) :: erder,uryg,urzg,vryg,vrzg
real(kind=8),dimension(2,2) :: acipa !el,a_temp
!el real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1
real(kind=8),dimension(4) :: muij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
+ integer xshift,yshift,zshift
!el integer :: num_conti,j1,j2
!el real(kind=8) :: a22,a23,a32,a33,dxi,dyi,dzi,dx_normi,dy_normi,&
!el dz_normi,xmedi,ymedi,zmedi
0.0d0,0.0d0,1.0d0/),shape(unmat))
! integer :: maxconts=nres/4
!el local variables
- integer :: k,i,j,iteli,itelj,kkk,l,kkll,m
+ integer :: k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap
real(kind=8) :: ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp
real(kind=8) :: ees,evdw1,eel_loc,aaa,bbb,ael3i
real(kind=8) :: dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj,&
dx_normj=dc_norm(1,j)
dy_normj=dc_norm(2,j)
dz_normj=dc_norm(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-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
+ 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
rrmij=1.0D0/rij
rij=dsqrt(rij)
+!C print *,xmedi,ymedi,zmedi,xj,yj,zj,boxxsize,rij
+ sss_ele_cut=sscale_ele(rij)
+ sss_ele_grad=sscagrad_ele(rij)
+! sss_ele_cut=1.0d0
+! sss_ele_grad=0.0d0
+! print *,sss_ele_cut,sss_ele_grad,&
+! (rij),r_cut_ele,rlamb_ele
+! if (sss_ele_cut.le.0.0) go to 128
+
rmij=1.0D0/rij
r3ij=rrmij*rmij
r6ij=r3ij*r3ij
eesij=el1+el2
! 12/26/95 - for the evaluation of multi-body H-bonding interactions
ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
- ees=ees+eesij
- evdw1=evdw1+evdwij
+ ees=ees+eesij*sss_ele_cut
+ evdw1=evdw1+evdwij*sss_ele_cut
!d write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
!d & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
!d & 1.0D0/dsqrt(rrmij),evdwij,eesij,
! Calculate contributions to the Cartesian gradient.
!
#ifdef SPLITELE
- facvdw=-6*rrmij*(ev1+evdwij)
- facel=-3*rrmij*(el1+eesij)
+ facvdw=-6*rrmij*(ev1+evdwij)*sss_ele_cut
+ facel=-3*rrmij*(el1+eesij)*sss_ele_cut
fac1=fac
erij(1)=xj*rmij
erij(2)=yj*rmij
!
! Radial derivatives. First process both termini of the fragment (i,j)
!
- ggg(1)=facel*xj
- ggg(2)=facel*yj
- ggg(3)=facel*zj
+ ggg(1)=facel*xj+sss_ele_grad*rmij*eesij*xj
+ ggg(2)=facel*yj+sss_ele_grad*rmij*eesij*yj
+ ggg(3)=facel*zj+sss_ele_grad*rmij*eesij*zj
+
! do k=1,3
! ghalf=0.5D0*ggg(k)
! gelc(k,i)=gelc(k,i)+ghalf
!grad gelc(l,k)=gelc(l,k)+ggg(l)
!grad enddo
!grad enddo
- ggg(1)=facvdw*xj
- ggg(2)=facvdw*yj
- ggg(3)=facvdw*zj
+ ggg(1)=facvdw*xj+sss_ele_grad*rmij*evdwij*xj
+ ggg(2)=facvdw*yj+sss_ele_grad*rmij*evdwij*yj
+ ggg(3)=facvdw*zj+sss_ele_grad*rmij*evdwij*zj
! do k=1,3
! ghalf=0.5D0*ggg(k)
! gvdwpp(k,i)=gvdwpp(k,i)+ghalf
!grad enddo
!grad enddo
#else
- facvdw=ev1+evdwij
- facel=el1+eesij
+ facvdw=(ev1+evdwij)*sss_ele_cut
+ facel=(el1+eesij)*sss_ele_cut
fac1=fac
fac=-3*rrmij*(facvdw+facvdw+facel)
erij(1)=xj*rmij
!
! Radial derivatives. First process both termini of the fragment (i,j)
!
- ggg(1)=fac*xj
- ggg(2)=fac*yj
- ggg(3)=fac*zj
+ ggg(1)=fac*xj+sss_ele_grad*rmij*(eesij+evdwij)*xj
+ ggg(2)=fac*yj+sss_ele_grad*rmij*(eesij+evdwij)*yj
+ ggg(3)=fac*zj+sss_ele_grad*rmij*(eesij+evdwij)*zj
! do k=1,3
! ghalf=0.5D0*ggg(k)
! gelc(k,i)=gelc(k,i)+ghalf
!d print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
!d & (dcosg(k),k=1,3)
do k=1,3
- ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
+ ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*sss_ele_cut
enddo
! do k=1,3
! ghalf=0.5D0*ggg(k)
do k=1,3
gelc(k,i)=gelc(k,i) &
+(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) &
- + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+ + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)&
+ *sss_ele_cut
gelc(k,j)=gelc(k,j) &
+(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) &
- + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+ + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)&
+ *sss_ele_cut
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
enddo
+
IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 &
.or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 &
.or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) &
+a33*muij(4)
! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
-
+! eel_loc_ij=0.0
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
'eelloc',i,j,eel_loc_ij
! if (energy_dec) write (iout,*) "a22",a22," a23",a23," a32",a32," a33",a33
! if (energy_dec) write (iout,*) "muij",muij
! write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
-
- eel_loc=eel_loc+eel_loc_ij
+
+ eel_loc=eel_loc+eel_loc_ij*sss_ele_cut
! Partial derivatives in virtual-bond dihedral angles gamma
if (i.gt.1) &
gel_loc_loc(i-1)=gel_loc_loc(i-1)+ &
- a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j) &
- +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)
+ (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j) &
+ +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)) &
+ *sss_ele_cut
gel_loc_loc(j-1)=gel_loc_loc(j-1)+ &
- a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) &
- +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
+ (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) &
+ +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)) &
+ *sss_ele_cut
! Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
- do l=1,3
- ggg(l)=agg(l,1)*muij(1)+ &
- agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
+! do l=1,3
+! ggg(1)=(agg(1,1)*muij(1)+ &
+! agg(1,2)*muij(2)+agg(1,3)*muij(3)+agg(1,4)*muij(4)) &
+! *sss_ele_cut &
+! +eel_loc_ij*sss_ele_grad*rmij*xj
+! ggg(2)=(agg(2,1)*muij(1)+ &
+! agg(2,2)*muij(2)+agg(2,3)*muij(3)+agg(2,4)*muij(4)) &
+! *sss_ele_cut &
+! +eel_loc_ij*sss_ele_grad*rmij*yj
+! ggg(3)=(agg(3,1)*muij(1)+ &
+! agg(3,2)*muij(2)+agg(3,3)*muij(3)+agg(3,4)*muij(4)) &
+! *sss_ele_cut &
+! +eel_loc_ij*sss_ele_grad*rmij*zj
+ xtemp(1)=xj
+ xtemp(2)=yj
+ xtemp(3)=zj
+
+ do l=1,3
+ ggg(l)=(agg(l,1)*muij(1)+ &
+ agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))&
+ *sss_ele_cut &
+ +eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
+
gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
!grad ghalf=0.5d0*ggg(l)
!grad enddo
! Remaining derivatives of eello
do l=1,3
- gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+ &
- aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
- gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+ &
- aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
- gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+ &
- aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
- gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+ &
- aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
+ gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+ &
+ aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))&
+ *sss_ele_cut
+!+eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
+ gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+ &
+ aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3) &
+ +aggi1(l,4)*muij(4))&
+ *sss_ele_cut
+!+eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
+ gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+ &
+ aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))&
+ *sss_ele_cut
+!+eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
+ gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+ &
+ aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3) &
+ +aggj1(l,4)*muij(4))&
+ *sss_ele_cut
+!+eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
enddo
ENDIF
! Change 12/26/95 to calculate four-body contributions to H-bonding energy
ees0mij=0
endif
! ees0mij=0.0D0
- ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
- ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
+ ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij) &
+ *sss_ele_cut
+
+ ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) &
+ *sss_ele_cut
+
! Diagnostics. Comment out or remove after debugging!
! ees0p(num_conti,i)=0.5D0*fac3*ees0pij
! ees0m(num_conti,i)=0.5D0*fac3*ees0mij
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
+ gggp(1)=gggp(1)+ees0pijp*xj &
+ +ees0p(num_conti,i)/sss_ele_cut*rmij*xj*sss_ele_grad
+ gggp(2)=gggp(2)+ees0pijp*yj &
+ +ees0p(num_conti,i)/sss_ele_cut*rmij*yj*sss_ele_grad
+ gggp(3)=gggp(3)+ees0pijp*zj &
+ +ees0p(num_conti,i)/sss_ele_cut*rmij*zj*sss_ele_grad
+
+ gggm(1)=gggm(1)+ees0mijp*xj &
+ +ees0m(num_conti,i)/sss_ele_cut*rmij*xj*sss_ele_grad
+
+ gggm(2)=gggm(2)+ees0mijp*yj &
+ +ees0m(num_conti,i)/sss_ele_cut*rmij*yj*sss_ele_grad
+
+ gggm(3)=gggm(3)+ees0mijp*zj &
+ +ees0m(num_conti,i)/sss_ele_cut*rmij*zj*sss_ele_grad
+
! Derivatives due to the contact function
gacont_hbr(1,num_conti,i)=fprimcont*xj
gacont_hbr(2,num_conti,i)=fprimcont*yj
!grad ghalfm=0.5D0*gggm(k)
gacontp_hb1(k,num_conti,i)= & !ghalfp+
(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i)) &
- + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+ + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) &
+ *sss_ele_cut
+
gacontp_hb2(k,num_conti,i)= & !ghalfp+
(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j)) &
- + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
- gacontp_hb3(k,num_conti,i)=gggp(k)
+ + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)&
+ *sss_ele_cut
+
+ gacontp_hb3(k,num_conti,i)=gggp(k) &
+ *sss_ele_cut
+
gacontm_hb1(k,num_conti,i)= & !ghalfm+
(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i)) &
- + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+ + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) &
+ *sss_ele_cut
+
gacontm_hb2(k,num_conti,i)= & !ghalfm+
(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j)) &
- + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
- gacontm_hb3(k,num_conti,i)=gggm(k)
+ + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) &
+ *sss_ele_cut
+
+ gacontm_hb3(k,num_conti,i)=gggm(k) &
+ *sss_ele_cut
+
enddo
! Diagnostics. Comment out or remove after debugging!
!diag do k=1,3
enddo
endif
endif
+ 128 continue
! t_eelecij=t_eelecij+MPI_Wtime()-time00
return
end subroutine eelecij
! include 'COMMON.CONTROL'
real(kind=8),dimension(3) :: ggg
!el local variables
- integer :: i,iint,j,k,iteli,itypj
+ 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
+ e1,e2,evdwij,rij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
+ integer xshift,yshift,zshift
evdw2=0.0D0
evdw2_14=0.0d0
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(i)
! yj=c(2,nres+j)-yi
! zj=c(3,nres+j)-zi
! 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)-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)
+ rij=dsqrt(1.0d0/rrij)
+ sss_ele_cut=sscale_ele(rij)
+ sss_ele_grad=sscagrad_ele(rij)
+! print *,sss_ele_cut,sss_ele_grad,&
+! (rij),r_cut_ele,rlamb_ele
+ if (sss_ele_cut.le.0.0) cycle
fac=rrij**expon2
e1=fac*fac*aad(itypj,iteli)
e2=fac*bad(itypj,iteli)
if (iabs(j-i) .le. 2) then
e1=scal14*e1
e2=scal14*e2
- evdw2_14=evdw2_14+e1+e2
+ evdw2_14=evdw2_14+(e1+e2)*sss_ele_cut
endif
evdwij=e1+e2
- evdw2=evdw2+evdwij
+ evdw2=evdw2+evdwij*sss_ele_cut
! if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)') &
! 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),&
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
!
! Calculate contributions to the gradient in the virtual-bond and SC vectors.
!
- fac=-(evdwij+e1)*rrij
+ fac=-(evdwij+e1)*rrij*sss_ele_cut
+ fac=fac+evdwij*sss_ele_grad/rij/expon
ggg(1)=xj*fac
ggg(2)=yj*fac
ggg(3)=zj*fac
! if (.not.allocated(gradbx)) allocate(gradbx(3,nres)) !(3,maxres)
do i=ibondp_start,ibondp_end
+ if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
- estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
- do j=1,3
- gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) &
- *dc(j,i-1)/vbld(i)
- enddo
- if (energy_dec) write(iout,*) &
- "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
+!C estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
+!C do j=1,3
+!C gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) &
+!C *dc(j,i-1)/vbld(i)
+!C enddo
+!C if (energy_dec) write(iout,*) &
+!C "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
+ diff = vbld(i)-vbldpDUM
else
diff = vbld(i)-vbldp0
+ endif
if (energy_dec) write (iout,'(a7,i5,4f7.3)') &
"estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
estr=estr+diff*diff
gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
enddo
! write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3)
- endif
+! endif
enddo
estr=0.5d0*AKP*estr+estr1
!
etors=0.0D0
do i=iphi_start,iphi_end
if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 &
+ .or. itype(i-3).eq.ntyp1 &
.or. itype(i).eq.ntyp1) cycle
etors_ii=0.0D0
if (iabs(itype(i)).eq.20) then
do i=iphid_start,iphid_end
etors_d_ii=0.0D0
if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 &
+ .or. itype(i-3).eq.ntyp1 &
.or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
itori=itortyp(itype(i-2))
itori1=itortyp(itype(i-1))
! " sigder",sigder
! write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
! write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
+!C print *,sss_ele_cut,'in sc_grad'
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)
+ gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*sss_ele_cut
+!C print *,'gg',k,gg(k)
enddo
! write (iout,*) "gg",(gg(k),k=1,3)
do k=1,3
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
+ +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv &
+ *sss_ele_cut
+
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
+ +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv &
+ *sss_ele_cut
+
! write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
! +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
! write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
! Check the gradient of the virtual-bond and SC vectors in the internal
! coordinates.
!
- aincr=1.0d-7
- aincr2=5.0d-8
+ aincr=1.0d-6
+ aincr2=5.0d-7
call cartder
write (iout,'(a)') '**************** dx/dalpha'
write (iout,'(a)')
nf=0
nfl=0
call zerograd
- aincr=1.0D-7
- print '(a)','CG processor',me,' calling CHECK_CART.'
+ aincr=1.0D-5
+ print '(a)','CG processor',me,' calling CHECK_CART.',aincr
nf=0
icall=0
call geom_to_var(nvar,x)
enddo
return
end subroutine check_ecart
+#ifdef CARGRAD
!-----------------------------------------------------------------------------
subroutine check_ecartint
! Check the gradient of the energy in Cartesian coordinates.
!el integer :: icall
!el common /srutu/ icall
real(kind=8),dimension(6) :: ggg,ggg1
- real(kind=8),dimension(3) :: cc,xx,ddc,ddx
+ real(kind=8),dimension(3) :: cc,xx,ddc,ddx,ddc1,ddcn
real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
- real(kind=8),dimension(3) :: dcnorm_safe,dxnorm_safe
+ real(kind=8),dimension(3) :: dcnorm_safe1,dcnorm_safe2,dxnorm_safe
real(kind=8),dimension(6,0:nres) :: grad_s,grad_s1 !(6,0:maxres)
real(kind=8),dimension(nres) :: phi_temp,theta_temp,alph_temp,omeg_temp !(maxres)
real(kind=8),dimension(0:n_ene) :: energia,energia1
write(iout,*) 'Calling CHECK_ECARTINT.'
nf=0
icall=0
+ write (iout,*) "Before geom_to_var"
+ call geom_to_var(nvar,x)
+ write (iout,*) "after geom_to_var"
+ write (iout,*) "split_ene ",split_ene
+ call flush(iout)
+ if (.not.split_ene) then
+ write(iout,*) 'Calling CHECK_ECARTINT if'
+ call etotal(energia)
+!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
+ etot=energia(0)
+ write (iout,*) "etot",etot
+ call flush(iout)
+!el call enerprint(energia)
+!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
+ call flush(iout)
+ write (iout,*) "enter cartgrad"
+ call flush(iout)
+ call cartgrad
+!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
+ write (iout,*) "exit cartgrad"
+ call flush(iout)
+ icall =1
+ do i=1,nres
+ write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
+ enddo
+ do j=1,3
+ grad_s(j,0)=gcart(j,0)
+ enddo
+!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
+ do i=1,nres
+ do j=1,3
+ grad_s(j,i)=gcart(j,i)
+ grad_s(j+3,i)=gxcart(j,i)
+ enddo
+ enddo
+ else
+write(iout,*) 'Calling CHECK_ECARTIN else.'
+!- split gradient check
+ call zerograd
+ call etotal_long(energia)
+!el call enerprint(energia)
+ call flush(iout)
+ write (iout,*) "enter cartgrad"
+ call flush(iout)
+ call cartgrad
+ write (iout,*) "exit cartgrad"
+ call flush(iout)
+ icall =1
+ write (iout,*) "longrange grad"
+ do i=1,nres
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
+ (gxcart(j,i),j=1,3)
+ enddo
+ do j=1,3
+ grad_s(j,0)=gcart(j,0)
+ enddo
+ do i=1,nres
+ do j=1,3
+ grad_s(j,i)=gcart(j,i)
+ grad_s(j+3,i)=gxcart(j,i)
+ enddo
+ enddo
+ call zerograd
+ call etotal_short(energia)
+!el call enerprint(energia)
+ call flush(iout)
+ write (iout,*) "enter cartgrad"
+ call flush(iout)
+ call cartgrad
+ write (iout,*) "exit cartgrad"
+ call flush(iout)
+ icall =1
+ write (iout,*) "shortrange grad"
+ do i=1,nres
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
+ (gxcart(j,i),j=1,3)
+ enddo
+ do j=1,3
+ grad_s1(j,0)=gcart(j,0)
+ enddo
+ do i=1,nres
+ do j=1,3
+ grad_s1(j,i)=gcart(j,i)
+ grad_s1(j+3,i)=gxcart(j,i)
+ enddo
+ enddo
+ endif
+ write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
+! do i=1,nres
+ do i=nnt,nct
+ do j=1,3
+ if (nnt.gt.1 .and. i.eq.nnt) ddc1(j)=c(j,1)
+ if (nct.lt.nres .and. i.eq.nct) ddcn(j)=c(j,nres)
+ ddc(j)=c(j,i)
+ ddx(j)=c(j,i+nres)
+ dcnorm_safe1(j)=dc_norm(j,i-1)
+ dcnorm_safe2(j)=dc_norm(j,i)
+ dxnorm_safe(j)=dc_norm(j,i+nres)
+ enddo
+ do j=1,3
+ c(j,i)=ddc(j)+aincr
+ if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=c(j,1)+aincr
+ if (nct.lt.nres .and. i.eq.nct) c(j,nres)=c(j,nres)+aincr
+ if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ call int_from_cart1(.false.)
+ if (.not.split_ene) then
+ call etotal(energia1)
+ etot1=energia1(0)
+ write (iout,*) "ij",i,j," etot1",etot1
+ else
+!- split gradient
+ call etotal_long(energia1)
+ etot11=energia1(0)
+ call etotal_short(energia1)
+ etot12=energia1(0)
+ endif
+!- end split gradient
+! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
+ c(j,i)=ddc(j)-aincr
+ if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)-aincr
+ if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)-aincr
+ if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ call int_from_cart1(.false.)
+ if (.not.split_ene) then
+ call etotal(energia1)
+ etot2=energia1(0)
+ write (iout,*) "ij",i,j," etot2",etot2
+ ggg(j)=(etot1-etot2)/(2*aincr)
+ else
+!- split gradient
+ call etotal_long(energia1)
+ etot21=energia1(0)
+ ggg(j)=(etot11-etot21)/(2*aincr)
+ call etotal_short(energia1)
+ etot22=energia1(0)
+ ggg1(j)=(etot12-etot22)/(2*aincr)
+!- end split gradient
+! write (iout,*) "etot21",etot21," etot22",etot22
+ endif
+! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
+ c(j,i)=ddc(j)
+ if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)
+ if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)
+ if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i-1)=dcnorm_safe1(j)
+ dc_norm(j,i)=dcnorm_safe2(j)
+ dc_norm(j,i+nres)=dxnorm_safe(j)
+ enddo
+ do j=1,3
+ c(j,i+nres)=ddx(j)+aincr
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ call int_from_cart1(.false.)
+ if (.not.split_ene) then
+ call etotal(energia1)
+ etot1=energia1(0)
+ else
+!- split gradient
+ call etotal_long(energia1)
+ etot11=energia1(0)
+ call etotal_short(energia1)
+ etot12=energia1(0)
+ endif
+!- end split gradient
+ c(j,i+nres)=ddx(j)-aincr
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ call int_from_cart1(.false.)
+ if (.not.split_ene) then
+ call etotal(energia1)
+ etot2=energia1(0)
+ ggg(j+3)=(etot1-etot2)/(2*aincr)
+ else
+!- split gradient
+ call etotal_long(energia1)
+ etot21=energia1(0)
+ ggg(j+3)=(etot11-etot21)/(2*aincr)
+ call etotal_short(energia1)
+ etot22=energia1(0)
+ ggg1(j+3)=(etot12-etot22)/(2*aincr)
+!- end split gradient
+ endif
+! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
+ c(j,i+nres)=ddx(j)
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dxnorm_safe(j)
+ call int_from_cart1(.false.)
+ enddo
+ write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+ i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
+ if (split_ene) then
+ write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+ i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),&
+ k=1,6)
+ write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+ i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),&
+ ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
+ endif
+ enddo
+ return
+ end subroutine check_ecartint
+#else
+!-----------------------------------------------------------------------------
+ subroutine check_ecartint
+! Check the gradient of the energy in Cartesian coordinates.
+ use io_base, only: intout
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CONTROL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.VAR'
+! include 'COMMON.CONTACTS'
+! include 'COMMON.MD'
+! include 'COMMON.LOCAL'
+! include 'COMMON.SPLITELE'
+ use comm_srutu
+!el integer :: icall
+!el common /srutu/ icall
+ real(kind=8),dimension(6) :: ggg,ggg1
+ real(kind=8),dimension(3) :: cc,xx,ddc,ddx
+ real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
+ real(kind=8),dimension(3) :: dcnorm_safe,dxnorm_safe
+ real(kind=8),dimension(6,0:nres) :: grad_s,grad_s1 !(6,0:maxres)
+ real(kind=8),dimension(nres) :: phi_temp,theta_temp,alph_temp,omeg_temp !(maxres)
+ real(kind=8),dimension(0:n_ene) :: energia,energia1
+ integer :: uiparm(1)
+ real(kind=8) :: urparm(1)
+!EL external fdum
+ integer :: i,j,k,nf
+ real(kind=8) :: rlambd,aincr,etot,etot1,etot11,etot12,etot2,&
+ etot21,etot22
+ r_cut=2.0d0
+ rlambd=0.3d0
+ icg=1
+ nf=0
+ nfl=0
+ call intout
+! call intcartderiv
+! call checkintcartgrad
+ call zerograd
+ aincr=2.0D-5
+ write(iout,*) 'Calling CHECK_ECARTINT.',aincr
+ nf=0
+ icall=0
call geom_to_var(nvar,x)
if (.not.split_ene) then
call etotal(energia)
enddo
return
end subroutine check_ecartint
+#endif
!-----------------------------------------------------------------------------
subroutine check_eint
! Check the gradient of energy in internal coordinates.
endif
return
end function sscale
+ real(kind=8) function sscale_grad(r)
+! include "COMMON.SPLITELE"
+ real(kind=8) :: r,gamm
+ if(r.lt.r_cut-rlamb) then
+ sscale_grad=0.0d0
+ else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+ gamm=(r-(r_cut-rlamb))/rlamb
+ sscale_grad=gamm*(6*gamm-6.0d0)/rlamb
+ else
+ sscale_grad=0d0
+ endif
+ return
+ end function sscale_grad
+
+!!!!!!!!!! PBCSCALE
+ real(kind=8) function sscale_ele(r)
+! include "COMMON.SPLITELE"
+ real(kind=8) :: r,gamm
+ if(r.lt.r_cut_ele-rlamb_ele) then
+ sscale_ele=1.0d0
+ else if(r.le.r_cut_ele.and.r.ge.r_cut_ele-rlamb_ele) then
+ gamm=(r-(r_cut_ele-rlamb_ele))/rlamb_ele
+ sscale_ele=1.0d0+gamm*gamm*(2*gamm-3.0d0)
+ else
+ sscale_ele=0d0
+ endif
+ return
+ end function sscale_ele
+
+ real(kind=8) function sscagrad_ele(r)
+ real(kind=8) :: r,gamm
+! include "COMMON.SPLITELE"
+ if(r.lt.r_cut_ele-rlamb_ele) then
+ sscagrad_ele=0.0d0
+ else if(r.le.r_cut_ele.and.r.ge.r_cut_ele-rlamb_ele) then
+ gamm=(r-(r_cut_ele-rlamb_ele))/rlamb_ele
+ sscagrad_ele=gamm*(6*gamm-6.0d0)/rlamb_ele
+ else
+ sscagrad_ele=0.0d0
+ endif
+ return
+ end function sscagrad_ele
+!!!!!!!!!!!!!!!
!-----------------------------------------------------------------------------
subroutine elj_long(evdw)
!
! include 'COMMON.CONTROL'
logical :: lprn
!el local variables
- integer :: iint,itypi,itypi1,itypj
+ integer :: iint,itypi,itypi1,itypj,subchap
real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig,sig0ij,rij_shift
- real(kind=8) :: sss,e1,e2,evdw
+ real(kind=8) :: sss,e1,e2,evdw,sss_grad
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
+
evdw=0.0D0
!cccc energy_dec=.false.
! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
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
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
alf1=alp(itypi)
alf2=alp(itypj)
alf12=0.5D0*(alf1+alf2)
- 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)
+! Searching for nearest neighbour
+ 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
+
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)
sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
-
+ sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
+ if (sss_ele_cut.le.0.0) cycle
if (sss.lt.1.0d0) then
! Calculate angle-dependent terms of energy and contributions to their
! write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
! & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij*(1.0d0-sss)
+ evdw=evdw+evdwij*(1.0d0-sss)*sss_ele_cut
if (lprn) then
sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac
+ fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
+ /sigma(itypi,itypj)*rij-sss_grad/(1.0-sss)*rij &
+ /sigmaii(itypi,itypj))
! fac=0.0d0
! Calculate the radial part of the gradient
gg(1)=xj*fac
! include 'COMMON.CONTROL'
logical :: lprn
!el local variables
- integer :: iint,itypi,itypi1,itypj
+ integer :: iint,itypi,itypi1,itypj,subchap
real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig0ij,sig
- real(kind=8) :: sss,e1,e2,evdw,rij_shift
+ real(kind=8) :: sss,e1,e2,evdw,rij_shift,sss_grad
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
evdw=0.0D0
!cccc energy_dec=.false.
! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
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
+ 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)
+
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
alf1=alp(itypi)
alf2=alp(itypj)
alf12=0.5D0*(alf1+alf2)
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+! 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)
+! Searching for nearest neighbour
+ 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
+
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)
sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+ sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
+ sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ if (sss_ele_cut.le.0.0) cycle
if (sss.gt.0.0d0) then
! write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
! & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij*sss
+ evdw=evdw+evdwij*sss*sss_ele_cut
if (lprn) then
sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac
+ fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
+ /sigma(itypi,itypj)*rij+sss_grad/sss*rij &
+ /sigmaii(itypi,itypj))
+
! fac=0.0d0
! Calculate the radial part of the gradient
gg(1)=xj*fac
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))*scalfac
+ gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac&
+ *sss_ele_cut
enddo
! write (iout,*) "gg",(gg(k),k=1,3)
do k=1,3
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*scalfac
+ +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac&
+ *sss_ele_cut
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*scalfac
+ +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac&
+ *sss_ele_cut
! write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
! & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
! write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
!
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
- use MD_data, only: totT
+ use MD_data, only: totT,usampl,eq_time
#ifndef ISNAN
external proc_proc
#ifdef WINPGI
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
use energy_data
- use MD_data, only: totT
+ use MD_data, only: totT,usampl,eq_time
#ifdef MPI
include 'mpif.h'
#endif
(gxcart(j,i),j=1,3)
enddo
#endif
+#ifdef CARGRAD
+#ifdef DEBUG
+ write (iout,*) "CARGRAD"
+#endif
+ do i=nres,1,-1
+ do j=1,3
+ gcart(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
+! gcart_new(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
+ enddo
+! write (iout,'(i5,3f10.5,5x,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), &
+! (gcart_new(j,i),j=1,3),(gxcart(j,i),j=1,3)
+ enddo
+! Correction: dummy residues
+ if (nnt.gt.1) then
+ do j=1,3
+! gcart_new(j,nnt)=gcart_new(j,nnt)+gcart_new(j,1)
+ gcart(j,nnt)=gcart(j,nnt)+gcart(j,1)
+ enddo
+ endif
+ if (nct.lt.nres) then
+ do j=1,3
+! gcart_new(j,nct)=gcart_new(j,nct)+gcart_new(j,nres)
+ gcart(j,nct)=gcart(j,nct)+gcart(j,nres)
+ enddo
+ endif
+#endif
#ifdef TIMING
time_cartgrad=time_cartgrad+MPI_Wtime()-time00
#endif
! Obtaining the gamma derivatives from sine derivative
if (phi(i).gt.-pi4.and.phi(i).le.pi4.or. &
phi(i).gt.pi34.and.phi(i).le.pi.or. &
- phi(i).gt.-pi.and.phi(i).le.-pi34) then
+ phi(i).ge.-pi.and.phi(i).le.-pi34) then
call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
call vecpr(dc_norm(1,i-3),dc_norm(1,i-1),vp2)
call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
!-----------------------------------------------------------------------------
subroutine alloc_ener_arrays
!EL Allocation of arrays used by module energy
-
+ use MD_data, only: mset
!el local variables
integer :: i,j
allocate(mu(2,nres))
allocate(muder(2,nres))
allocate(Ub2(2,nres))
- do i=1,nres
- Ub2(1,i)=0.0d0
- Ub2(2,i)=0.0d0
- enddo
+ Ub2(1,:)=0.0d0
+ Ub2(2,:)=0.0d0
allocate(Ub2der(2,nres))
allocate(Ctobr(2,nres))
allocate(Ctobrder(2,nres))
! allocate(qinfrag(50,nprocs/20),wfrag(50,nprocs/20)) !(50,maxprocs/20)
! allocate(qinpair(100,nprocs/20),wpair(100,nprocs/20)) !(100,maxprocs/20)
allocate(mset(0:nprocs)) !(maxprocs/20)
- do i=0,nprocs
- mset(i)=0
- enddo
+ mset(:)=0
! allocate(ifrag(2,50,nprocs/20)) !(2,50,maxprocs/20)
! allocate(ipair(2,100,nprocs/20)) !(2,100,maxprocs/20)
allocate(dUdconst(3,0:nres))
! and side-chain vectors in theta or phi.
allocate(dyn_ssbond_ij(0:nres+4,0:nres+4))
!(maxres,maxres)
- do i=1,nres
- do j=i+1,nres
- dyn_ssbond_ij(i,j)=1.0d300
- enddo
- enddo
+! do i=1,nres
+! do j=i+1,nres
+ dyn_ssbond_ij(:,:)=1.0d300
+! enddo
+! enddo
if (nss.gt.0) then
allocate(idssb(nss),jdssb(nss))
endif
allocate(dyn_ss_mask(nres))
!(maxres)
- do i=1,nres
- dyn_ss_mask(i)=.false.
- enddo
+ dyn_ss_mask(:)=.false.
!----------------------
! common.sccor
! Parameters of the SCCOR term
projects / unres4.git / blobdiff