X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new-sep_barrier.F;h=44384a6048fa76830a8b441ddbfe1971a8b7c641;hb=4d4be87fb43102a6a5be442699a953ef327ecbde;hp=e0b1f44e7d761743517e3450d65fc282b50329ee;hpb=aedd2f72960d9cec69fcc7f0b6fd6555c7e22312;p=unres.git 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 e0b1f44..44384a6 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 @@ -6,7 +6,7 @@ C if(r.lt.r_cut-rlamb) then C sscale=1.0d0 C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then C gamm=(r-(r_cut-rlamb))/rlamb - sscalelip=1.0d0+r*r*(2*r-3.0d0) + sscalelip=1.0d0+r*r*(2.0d0*r-3.0d0) C else C sscale=0d0 C endif @@ -20,7 +20,7 @@ C if(r.lt.r_cut-rlamb) then C sscagrad=0.0d0 C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then C gamm=(r-(r_cut-rlamb))/rlamb - sscagradlip=r*(6*r-6.0d0) + sscagradlip=r*(6.0d0*r-6.0d0) C else C sscagrad=0.0d0 C endif @@ -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 @@ -1455,8 +1515,8 @@ C do i=iturn3_start,iturn3_end if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1 & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 - & .or. itype(i+4).eq.ntyp1 +C & .or. itype(i-1).eq.ntyp1 +C & .or. itype(i+4).eq.ntyp1 & ) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -1482,8 +1542,8 @@ C if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 & .or. itype(i+3).eq.ntyp1 & .or. itype(i+4).eq.ntyp1 - & .or. itype(i+5).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 +C & .or. itype(i+5).eq.ntyp1 +C & .or. itype(i-1).eq.ntyp1 & ) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -1511,8 +1571,8 @@ c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c do i=iatel_s,iatel_e if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 - & .or. itype(i+2).eq.ntyp1 - & .or. itype(i-1).eq.ntyp1 +C & .or. itype(i+2).eq.ntyp1 +C & .or. itype(i-1).eq.ntyp1 &) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -1533,8 +1593,8 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) do j=ielstart(i),ielend(i) if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 - & .or.itype(j+2).eq.ntyp1 - & .or.itype(j-1).eq.ntyp1 +C & .or.itype(j+2).eq.ntyp1 +C & .or.itype(j-1).eq.ntyp1 &) cycle call eelecij_scale(i,j,ees,evdw1,eel_loc) enddo ! j @@ -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 @@ -2614,9 +2683,8 @@ C Uncomment following three lines for Ca-p interactions 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.lt.1.0d0) then fac=rrij**expon2 e1=fac*fac*aad(itypj,iteli) @@ -2635,7 +2703,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)/expon ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac @@ -2691,6 +2759,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 +2790,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(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 @@ -2753,9 +2828,9 @@ 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))) - if (sss.gt.0.0d0) then + sss=sscale(1.0d0/(dsqrt(rrij))) + sssgrad=sscagrad(1.0d0/(dsqrt(rrij))) + if (sss.gt.0.0d0) then fac=rrij**expon2 e1=fac*fac*aad(itypj,iteli) @@ -2769,11 +2844,12 @@ C Uncomment following three lines for Ca-p interactions evdw2=evdw2+evdwij*sss if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))') & 'evdw2',i,j,sss,evdwij + 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)/expon ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac