From: Adam Sieradzan Date: Tue, 26 Apr 2016 07:23:07 +0000 (+0200) Subject: bug fix in respa and readpdb X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=eab68ce2fb626e693cd1367fe2bcb717b98c587f bug fix in respa and readpdb --- diff --git a/source/unres/src_MD-M/energy_p_new-sep_barrier.F b/source/unres/src_MD-M/energy_p_new-sep_barrier.F index 3815c13..0a638ef 100644 --- a/source/unres/src_MD-M/energy_p_new-sep_barrier.F +++ b/source/unres/src_MD-M/energy_p_new-sep_barrier.F @@ -662,6 +662,29 @@ c if (icall.eq.0) lprn=.false. if (yi.lt.0) yi=yi+boxysize zi=mod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize + if ((zi.gt.bordlipbot) + &.and.(zi.lt.bordliptop)) then +C the energy transfer exist + if (zi.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((zi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipi=sscalelip(fracinbuf) + ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick + elseif (zi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick) + sslipi=sscalelip(fracinbuf) + ssgradlipi=sscagradlip(fracinbuf)/lipbufthick + else + sslipi=1.0d0 + ssgradlipi=0.0 + endif + else + sslipi=0.0d0 + ssgradlipi=0.0 + endif + dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -707,12 +730,12 @@ C the energy transfer exist if (zj.lt.buflipbot) then C what fraction I am in fracinbuf=1.0d0- - & ((positi-bordlipbot)/lipbufthick) + & ((zj-bordlipbot)/lipbufthick) C lipbufthick is thickenes of lipid buffore sslipj=sscalelip(fracinbuf) ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick - elseif (zi.gt.bufliptop) then - fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) + elseif (zj.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) sslipj=sscalelip(fracinbuf) ssgradlipj=sscagradlip(fracinbuf)/lipbufthick else @@ -723,10 +746,11 @@ C lipbufthick is thickenes of lipid buffore sslipj=0.0d0 ssgradlipj=0.0 endif + aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 - & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 + & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 - & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 xj_safe=xj @@ -788,6 +812,7 @@ cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -822,8 +847,12 @@ C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac - gg_lipi(3)=ssgradlipi*evdwij - gg_lipj(3)=ssgradlipj*evdwij + gg_lipi(3)=eps1*(eps2rt*eps2rt) + &*(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + &+faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi C Calculate angular part of the gradient. call sc_grad_scale(1.0d0-sss) endif @@ -853,6 +882,7 @@ C include 'COMMON.CALC' include 'COMMON.CONTROL' logical lprn + integer xshift,yshift,zshift evdw=0.0D0 ccccc energy_dec=.false. c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -873,6 +903,29 @@ c if (icall.eq.0) lprn=.false. if (yi.lt.0) yi=yi+boxysize zi=mod(zi,boxzsize) if (zi.lt.0) zi=zi+boxzsize + if ((zi.gt.bordlipbot) + &.and.(zi.lt.bordliptop)) then +C the energy transfer exist + if (zi.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((zi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipi=sscalelip(fracinbuf) + ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick + elseif (zi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick) + sslipi=sscalelip(fracinbuf) + ssgradlipi=sscagradlip(fracinbuf)/lipbufthick + else + sslipi=1.0d0 + ssgradlipi=0.0 + endif + else + sslipi=0.0d0 + ssgradlipi=0.0 + endif + dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -918,12 +971,12 @@ C the energy transfer exist if (zj.lt.buflipbot) then C what fraction I am in fracinbuf=1.0d0- - & ((positi-bordlipbot)/lipbufthick) + & ((zj-bordlipbot)/lipbufthick) C lipbufthick is thickenes of lipid buffore sslipj=sscalelip(fracinbuf) ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick - elseif (zi.gt.bufliptop) then - fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) + elseif (zj.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) sslipj=sscalelip(fracinbuf) ssgradlipj=sscagradlip(fracinbuf)/lipbufthick else @@ -935,9 +988,9 @@ C lipbufthick is thickenes of lipid buffore ssgradlipj=0.0 endif aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 - & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 + & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 - & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 + & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 xj_safe=xj yj_safe=yj @@ -998,6 +1051,7 @@ cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon + faclip=fac e1=fac*fac*aa e2=fac*bb evdwij=eps1*eps2rt*eps3rt*(e1+e2) @@ -1032,15 +1086,19 @@ C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac - gg_lipi(3)=ssgradlipi*evdwij - gg_lipj(3)=ssgradlipj*evdwij + gg_lipi(3)=eps1*(eps2rt*eps2rt) + &*(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + &+faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi C Calculate angular part of the gradient. call sc_grad_scale(sss) endif enddo ! j enddo ! iint enddo ! i -c write (iout,*) "Number of loop steps in EGB:",ind +C write (iout,*) "Number of loop steps in EGB:",ind cccc energy_dec=.false. return end @@ -1317,6 +1375,8 @@ c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12 enddo do k=1,3 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac + gg_lipi(k)=gg_lipi(k)*scalfac + gg_lipj(k)=gg_lipj(k)*scalfac enddo c write (iout,*) "gg",(gg(k),k=1,3) do k=1,3 @@ -1589,6 +1649,7 @@ C 13-go grudnia roku pamietnego... double precision unmat(3,3) /1.0d0,0.0d0,0.0d0, & 0.0d0,1.0d0,0.0d0, & 0.0d0,0.0d0,1.0d0/ + integer xhift,yshift,zshift c time00=MPI_Wtime() cd write (iout,*) "eelecij",i,j C print *,"WCHODZE2" @@ -1652,8 +1713,8 @@ C print *,"WCHODZE2" rij=dsqrt(rij) rmij=1.0D0/rij c For extracting the short-range part of Evdwpp - sss=sscale(rij/rpp(iteli,itelj)) - sssgrad=sscagrad(rij/rpp(iteli,itelj)) + sss=sscale(rij) + sssgrad=sscagrad(rij) r3ij=rrmij*rmij r6ij=r3ij*r3ij cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj @@ -1778,9 +1839,9 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo - ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj) - ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj) - ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj) + ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj + ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj + ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj c do k=1,3 c ghalf=0.5D0*ggg(k) c gvdwpp(k,i)=gvdwpp(k,i)+ghalf @@ -1835,9 +1896,9 @@ c 9/28/08 AL Gradient compotents will be summed only at the end C ggg(1)=facvdw*xj C ggg(2)=facvdw*yj C ggg(3)=facvdw*zj - ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj) - ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj) - ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj) + ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj + ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj + ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj do k=1,3 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) @@ -2499,8 +2560,8 @@ c & ' ielend',ielend_vdw(i) rij=xj*xj+yj*yj+zj*zj rrmij=1.0D0/rij rij=dsqrt(rij) - sss=sscale(rij/rpp(iteli,itelj)) - sssgrad=sscagrad(rij/rpp(iteli,itelj)) + sss=sscale(rij) + sssgrad=sscagrad(rij) if (sss.gt.0.0d0) then rmij=1.0D0/rij r3ij=rrmij*rmij @@ -2518,9 +2579,9 @@ C C Calculate contributions to the Cartesian gradient. C facvdw=-6*rrmij*(ev1+evdwij)*sss - ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj) - ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj) - ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj) + ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj + ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj C ggg(1)=facvdw*xj C ggg(2)=facvdw*yj C ggg(3)=facvdw*zj @@ -2552,6 +2613,7 @@ C include 'COMMON.IOUNITS' include 'COMMON.CONTROL' dimension ggg(3) + integer xshift,yshift,zshift evdw2=0.0D0 evdw2_14=0.0d0 CD print '(a)','Enter ESCP KURWA' @@ -2581,6 +2643,13 @@ C Uncomment following three lines for Ca-p interactions 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(zi,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 @@ -2615,8 +2684,8 @@ C Uncomment following three lines for Ca-p interactions rrij=1.0D0/(xj*xj+yj*yj+zj*zj) - sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) - sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) + sss=sscale(1.0d0/(dsqrt(rrij))) + sssgrad=sscagrad(1.0d0/(dsqrt(rrij))) if (sss.lt.1.0d0) then fac=rrij**expon2 e1=fac*fac*aad(itypj,iteli) @@ -2635,7 +2704,7 @@ C Calculate contributions to the gradient in the virtual-bond and SC vectors. C fac=-(evdwij+e1)*rrij*(1.0d0-sss) - fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli) + fac=fac-(evdwij)*sssgrad*dsqrt(rrij) ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac @@ -2691,6 +2760,7 @@ C include 'COMMON.IOUNITS' include 'COMMON.CONTROL' dimension ggg(3) + integer xshift,yshift,zshift evdw2=0.0D0 evdw2_14=0.0d0 cd print '(a)','Enter ESCP' @@ -2721,6 +2791,12 @@ C Uncomment following three lines for Ca-p interactions 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(zi,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 @@ -2753,8 +2829,8 @@ C Uncomment following three lines for Ca-p interactions zj=zj_safe-zi endif rrij=1.0D0/(xj*xj+yj*yj+zj*zj) - sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) - sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli))) + sss=sscale(1.0d0/(dsqrt(rrij))) + sssgrad=sscagrad(1.0d0/(dsqrt(rrij))) if (sss.gt.0.0d0) then fac=rrij**expon2 @@ -2773,7 +2849,7 @@ C C Calculate contributions to the gradient in the virtual-bond and SC vectors. C fac=-(evdwij+e1)*rrij*sss - fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli) + fac=fac+(evdwij)*sssgrad*dsqrt(rrij) ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac diff --git a/source/unres/src_MD-M/energy_split-sep.F b/source/unres/src_MD-M/energy_split-sep.F index e15ab61..b3f924b 100644 --- a/source/unres/src_MD-M/energy_split-sep.F +++ b/source/unres/src_MD-M/energy_split-sep.F @@ -246,7 +246,12 @@ c &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6 if ((wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) .and. ipot.lt.6) then call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1) endif -C +C +C if (wliptran.gt.0) then +C call Eliptransfer(eliptran) +C else +C eliptran=0.0d0 +C endif C If performing constraint dynamics, call the constraint energy C after the equilibration time if(usampl.and.totT.gt.eq_time) then @@ -287,7 +292,7 @@ C energia(20)=Uconst+Uconst_back call sum_energy(energia,.true.) C call enerprint - write (iout,*) "Exit ETOTAL_LONG" +C write (iout,*) "Exit ETOTAL_LONG" call flush(iout) return end @@ -575,6 +580,12 @@ C else esccor=0.0d0 endif + if (wliptran.gt.0) then + call Eliptransfer(eliptran) + else + eliptran=0.0d0 + endif +C print *,eliptran,wliptran C C Put energy components into an array C @@ -611,10 +622,11 @@ C energia(17)=estr energia(19)=edihcnstr energia(21)=esccor - write (iout,*) "ETOTAL_SHORT before SUM_ENERGY" + energia(22)=eliptran +C write (iout,*) "ETOTAL_SHORT before SUM_ENERGY" call flush(iout) call sum_energy(energia,.true.) - write (iout,*) "Exit ETOTAL_SHORT" +C write (iout,*) "Exit ETOTAL_SHORT" call flush(iout) return end diff --git a/source/unres/src_MD-M/readpdb.F b/source/unres/src_MD-M/readpdb.F index 7aa8fd4..48bb75a 100644 --- a/source/unres/src_MD-M/readpdb.F +++ b/source/unres/src_MD-M/readpdb.F @@ -310,7 +310,9 @@ c enddo c enddiagnostic C makes copy of chains write (iout,*) "symetr", symetr - + do j=1,3 + dc(j,0)=c(j,1) + enddo if (symetr.gt.1) then call permut(symetr) nperm=1