X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.f90;h=949e6784d880a33d71d284b42875f6a32a978cc4;hb=cfa004813be91a2b6c43e3e5fd6821b31198b81f;hp=fdf457697c7664460b3e4a5db567913e728d0dee;hpb=05fce52032feb59f7e96cd897fb82ca2aa90a888;p=unres4.git diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index fdf4576..949e678 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -180,7 +180,7 @@ 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 @@ -213,6 +213,14 @@ #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 @@ -1269,9 +1277,11 @@ ! 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 @@ -1286,6 +1296,13 @@ 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) @@ -1334,9 +1351,46 @@ ! 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) @@ -1349,6 +1403,11 @@ ! 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 @@ -1380,7 +1439,7 @@ ! 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) @@ -1402,11 +1461,18 @@ 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 @@ -2625,6 +2691,12 @@ 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) @@ -2643,6 +2715,12 @@ 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) & @@ -2663,6 +2741,13 @@ 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) @@ -2704,11 +2789,14 @@ ! 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 @@ -2730,7 +2818,7 @@ 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,& @@ -2761,12 +2849,62 @@ 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 @@ -2786,8 +2924,8 @@ 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, @@ -2804,8 +2942,8 @@ ! 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 @@ -2813,9 +2951,10 @@ ! ! 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 @@ -2834,9 +2973,9 @@ !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 @@ -2856,8 +2995,8 @@ !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 @@ -2866,9 +3005,9 @@ ! ! 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 @@ -2911,7 +3050,7 @@ !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) @@ -2930,13 +3069,16 @@ 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 @@ -3125,26 +3267,48 @@ 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) @@ -3158,14 +3322,24 @@ !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 @@ -3247,8 +3421,12 @@ 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 @@ -3296,12 +3474,22 @@ 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 @@ -3315,18 +3503,30 @@ !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 @@ -3358,6 +3558,7 @@ enddo endif endif + 128 continue ! t_eelecij=t_eelecij+MPI_Wtime()-time00 return end subroutine eelecij @@ -3801,9 +4002,12 @@ ! 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 @@ -3815,6 +4019,12 @@ 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) @@ -3826,20 +4036,67 @@ ! 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)') & @@ -3847,7 +4104,8 @@ ! ! 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 @@ -4119,16 +4377,19 @@ ! 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 @@ -4136,7 +4397,7 @@ 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 ! @@ -9564,21 +9825,27 @@ ! " 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)) & @@ -9989,8 +10256,8 @@ ! 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)') @@ -10172,8 +10439,8 @@ 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) @@ -10228,6 +10495,257 @@ enddo return end subroutine check_ecart +#ifdef CARGRAD +!----------------------------------------------------------------------------- + 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,ddc1,ddcn + real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres) + 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 + 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=1.0D-5 + 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. @@ -10268,8 +10786,8 @@ ! call intcartderiv ! call checkintcartgrad call zerograd - aincr=1.0D-4 - write(iout,*) 'Calling CHECK_ECARTINT.' + aincr=2.0D-5 + write(iout,*) 'Calling CHECK_ECARTINT.',aincr nf=0 icall=0 call geom_to_var(nvar,x) @@ -10466,6 +10984,7 @@ enddo return end subroutine check_ecartint +#endif !----------------------------------------------------------------------------- subroutine check_eint ! Check the gradient of energy in internal coordinates. @@ -10672,6 +11191,49 @@ 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) ! @@ -11290,9 +11852,12 @@ ! 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 @@ -11307,6 +11872,12 @@ 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) @@ -11337,16 +11908,58 @@ 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 @@ -11377,7 +11990,7 @@ ! 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) @@ -11399,6 +12012,9 @@ 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 @@ -11435,9 +12051,11 @@ ! 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 @@ -11452,6 +12070,18 @@ 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) @@ -11482,15 +12112,61 @@ 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 @@ -11522,7 +12198,7 @@ ! 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) @@ -11544,6 +12220,10 @@ 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 @@ -13039,16 +13719,19 @@ 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)) @@ -13072,7 +13755,7 @@ ! ! 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 @@ -13752,7 +14435,7 @@ ! 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 @@ -13839,6 +14522,32 @@ (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 @@ -14119,7 +14828,7 @@ ! 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) @@ -15896,7 +16605,7 @@ !----------------------------------------------------------------------------- subroutine alloc_ener_arrays !EL Allocation of arrays used by module energy - + use MD_data, only: mset !el local variables integer :: i,j