From aa456a8235be939a2f24ba02c7df7dfee1ecebd7 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Sat, 29 Apr 2017 17:36:14 +0200 Subject: [PATCH] energy fr PBC workin ene,grad,symp_chek OK, for respa also OK --- source/unres/energy.f90 | 489 ++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 402 insertions(+), 87 deletions(-) diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index 949e678..1a69db0 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -2861,6 +2861,7 @@ 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 @@ -12634,6 +12635,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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_scale(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) @@ -12652,6 +12659,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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_scale(i,i+3,ees,evdw1,eel_loc) if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) & @@ -12672,6 +12685,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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) @@ -12712,11 +12731,15 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! 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,sss_grad + 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 @@ -12735,7 +12758,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 0.0d0,1.0d0,0.0d0,& 0.0d0,0.0d0,1.0d0/),shape(unmat)) !el local variables - integer :: i,j,k,l,iteli,itelj,kkk,kkll,m + integer :: i,j,k,l,iteli,itelj,kkk,kkll,m,isubchap real(kind=8) :: aaa,bbb,ael6i,ael3i,dxj,dyj,dzj real(kind=8) :: xj,yj,zj,rij,rrmij,rmij,sss,r3ij,r6ij,fac real(kind=8) :: cosa,cosb,cosg,ev1,ev2,fac3,fac4,evdwij @@ -12784,15 +12807,63 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 + 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 rrmij=1.0D0/rij rij=dsqrt(rij) rmij=1.0D0/rij ! For extracting the short-range part of Evdwpp sss=sscale(rij/rpp(iteli,itelj)) + sss_ele_cut=sscale_ele(rij) + sss_ele_grad=sscagrad_ele(rij) + sss_grad=sscale_grad((rij/rpp(iteli,itelj))) +! sss_ele_cut=1.0d0 +! sss_ele_grad=0.0d0 + if (sss_ele_cut.le.0.0) go to 128 r3ij=rrmij*rmij r6ij=r3ij*r3ij @@ -12812,8 +12883,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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*(1.0d0-sss) + ees=ees+eesij*sss_ele_cut + evdw1=evdw1+evdwij*(1.0d0-sss)*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, @@ -12828,8 +12899,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! Calculate contributions to the Cartesian gradient. ! #ifdef SPLITELE - facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss) - facel=-3*rrmij*(el1+eesij) + facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss_ele_cut + facel=-3*rrmij*(el1+eesij)*sss_ele_cut fac1=fac erij(1)=xj*rmij erij(2)=yj*rmij @@ -12837,9 +12908,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! ! 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 @@ -12858,9 +12929,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !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*(1.0d0-sss) & + -evdwij*sss_ele_cut/rij*sss_grad/rpp(iteli,itelj)*xj + ggg(2)=facvdw*yj+sss_ele_grad*rmij*evdwij*yj*(1.0d0-sss) & + -evdwij*sss_ele_cut/rij*sss_grad/rpp(iteli,itelj)*yj + ggg(3)=facvdw*zj+sss_ele_grad*rmij*evdwij*zj*(1.0d0-sss) & + -evdwij*sss_ele_cut/rij*sss_grad/rpp(iteli,itelj)*zj ! do k=1,3 ! ghalf=0.5D0*ggg(k) ! gvdwpp(k,i)=gvdwpp(k,i)+ghalf @@ -12880,8 +12954,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !grad enddo !grad enddo #else - facvdw=ev1+evdwij*(1.0d0-sss) - facel=el1+eesij + facvdw=(ev1+evdwij)*(1.0d0-sss)*sss_ele_cut + facel=(el1+eesij)*sss_ele_cut fac1=fac fac=-3*rrmij*(facvdw+facvdw+facel) erij(1)=xj*rmij @@ -12935,7 +13009,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !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) @@ -12954,10 +13028,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -13154,19 +13230,28 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 'eelloc',i,j,eel_loc_ij ! write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) !d - 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 + xtemp(1)=xj + xtemp(2)=yj + xtemp(3)=zj + ! 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) + 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) @@ -13180,14 +13265,22 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !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 + + 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 + + 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 + + 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 + enddo ENDIF ! Change 12/26/95 to calculate four-body contributions to H-bonding energy @@ -13268,8 +13361,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -13317,12 +13414,28 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 +! 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 @@ -13334,20 +13447,46 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! !grad ghalfp=0.5D0*gggp(k) !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) - 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) - 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) - 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) +! 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) +! 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) +! 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) +! 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) + 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) & + *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)& + *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) & + *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) & + *sss_ele_cut + + gacontm_hb3(k,num_conti,i)=gggm(k) & + *sss_ele_cut + enddo ENDIF ! wcorr endif ! num_conti.le.maxconts @@ -13370,6 +13509,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' enddo endif endif + 128 continue ! t_eelecij=t_eelecij+MPI_Wtime()-time00 return end subroutine eelecij_scale @@ -13400,11 +13540,15 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' real(kind=8) :: scal_el=0.5d0 #endif !el local variables - integer :: i,j,k,iteli,itelj,num_conti + integer :: i,j,k,iteli,itelj,num_conti,isubchap 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 + integer xshift,yshift,zshift + evdw1=0.0D0 ! write (iout,*) "iatel_s_vdw",iatel_s_vdw, @@ -13421,6 +13565,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 ! write (iout,*) 'i',i,' ielstart',ielstart_vdw(i), ! & ' ielend',ielend_vdw(i) @@ -13439,13 +13589,59 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 + 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 rrmij=1.0D0/rij rij=dsqrt(rij) sss=sscale(rij/rpp(iteli,itelj)) + sss_ele_cut=sscale_ele(rij) + sss_ele_grad=sscagrad_ele(rij) + sss_grad=sscale_grad((rij/rpp(iteli,itelj))) + if (sss_ele_cut.le.0.0) cycle if (sss.gt.0.0d0) then rmij=1.0D0/rij r3ij=rrmij*rmij @@ -13458,14 +13654,21 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' if (energy_dec) then write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss endif - evdw1=evdw1+evdwij*sss + evdw1=evdw1+evdwij*sss*sss_ele_cut ! ! Calculate contributions to the Cartesian gradient. ! - facvdw=-6*rrmij*(ev1+evdwij)*sss - ggg(1)=facvdw*xj - ggg(2)=facvdw*yj - ggg(3)=facvdw*zj + facvdw=-6*rrmij*(ev1+evdwij)*sss*sss_ele_cut +! ggg(1)=facvdw*xj +! ggg(2)=facvdw*yj +! ggg(3)=facvdw*zj + ggg(1)=facvdw*xj+sss_ele_grad*rmij*evdwij*xj*sss & + +evdwij*sss_ele_cut/rij*sss_grad/rpp(iteli,itelj)*xj + ggg(2)=facvdw*yj+sss_ele_grad*rmij*evdwij*yj*sss & + +evdwij*sss_ele_cut/rij*sss_grad/rpp(iteli,itelj)*yj + ggg(3)=facvdw*zj+sss_ele_grad*rmij*evdwij*zj*sss & + +evdwij*sss_ele_cut/rij*sss_grad/rpp(iteli,itelj)*zj + do k=1,3 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) @@ -13495,9 +13698,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! include 'COMMON.CONTROL' real(kind=8),dimension(3) :: ggg !el local variables - integer :: i,iint,j,k,iteli,itypj - real(kind=8) :: xi,yi,zi,xj,yj,zj,rrij,sss,fac,e1,e2 + integer :: i,iint,j,k,iteli,itypj,subchap + real(kind=8) :: xi,yi,zi,xj,yj,zj,rrij,sss,fac,e1,e2,sss_grad,rij real(kind=8) :: evdw2,evdw2_14,evdwij + real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,& + dist_temp, dist_init + evdw2=0.0D0 evdw2_14=0.0d0 !d print '(a)','Enter ESCP' @@ -13508,6 +13714,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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) @@ -13519,13 +13731,56 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! 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) + 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) - sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) - + 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 + sss=sscale((rij/rscp(itypj,iteli))) + sss_grad=sscale_grad(rij/rscp(itypj,iteli)) if (sss.lt.1.0d0) then fac=rrij**expon2 @@ -13534,16 +13789,18 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' if (iabs(j-i) .le. 2) then e1=scal14*e1 e2=scal14*e2 - evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss) + evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)*sss_ele_cut endif evdwij=e1+e2 - evdw2=evdw2+evdwij*(1.0d0-sss) + evdw2=evdw2+evdwij*(1.0d0-sss)*sss_ele_cut if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))') & 'evdw2',i,j,sss,evdwij ! ! Calculate contributions to the gradient in the virtual-bond and SC vectors. ! - fac=-(evdwij+e1)*rrij*(1.0d0-sss) + fac=-(evdwij+e1)*rrij*(1.0d0-sss)*sss_ele_cut + fac=fac+evdwij*sss_ele_grad/rij/expon*(1.0d0-sss)& + -evdwij*sss_ele_cut/rij/expon*sss_grad/rscp(itypj,iteli) ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac @@ -13600,9 +13857,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! include 'COMMON.CONTROL' real(kind=8),dimension(3) :: ggg !el local variables - integer :: i,iint,j,k,iteli,itypj - real(kind=8) :: xi,yi,zi,xj,yj,zj,rrij,sss,fac,e1,e2 + integer :: i,iint,j,k,iteli,itypj,subchap + real(kind=8) :: xi,yi,zi,xj,yj,zj,rrij,sss,fac,e1,e2,sss_grad,rij real(kind=8) :: evdw2,evdw2_14,evdwij + real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,& + dist_temp, dist_init + evdw2=0.0D0 evdw2_14=0.0d0 !d print '(a)','Enter ESCP' @@ -13613,6 +13873,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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) @@ -13624,13 +13890,59 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! 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 - rrij=1.0D0/(xj*xj+yj*yj+zj*zj) - - sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) +! 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 + sss=sscale(rij/rscp(itypj,iteli)) + sss_grad=sscale_grad(rij/rscp(itypj,iteli)) if (sss.gt.0.0d0) then fac=rrij**expon2 @@ -13639,16 +13951,19 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' if (iabs(j-i) .le. 2) then e1=scal14*e1 e2=scal14*e2 - evdw2_14=evdw2_14+(e1+e2)*sss + evdw2_14=evdw2_14+(e1+e2)*sss*sss_ele_cut endif evdwij=e1+e2 - evdw2=evdw2+evdwij*sss + evdw2=evdw2+evdwij*sss*sss_ele_cut if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))') & 'evdw2',i,j,sss,evdwij ! ! Calculate contributions to the gradient in the virtual-bond and SC vectors. ! - fac=-(evdwij+e1)*rrij*sss + fac=-(evdwij+e1)*rrij*sss*sss_ele_cut + fac=fac+evdwij*sss_ele_grad/rij/expon*sss & + +evdwij*sss_ele_cut/rij/expon*sss_grad/rscp(itypj,iteli) + ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac -- 1.7.9.5