X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=blobdiff_plain;f=source%2Fwham%2Fsrc-M%2Fenergy_p_new.F;h=b5f5705f17804e35a23c3be28055804723415063;hp=a57f604214e2ce995f838bafc90dbe9330bd9d5e;hb=bde71ffff7527add147c1bee4bad23f02b4d1ab7;hpb=6f352fb19a065386a573d715b09235bbdcfdd744 diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index a57f604..b5f5705 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -250,11 +250,28 @@ C & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) & +wliptran*gliptranc(j,i) + & +welec*gshieldc(j,i) + & +welec*gshieldc_loc(j,i) + & +wcorr*gshieldc_ec(j,i) + & +wcorr*gshieldc_loc_ec(j,i) + & +wturn3*gshieldc_t3(j,i) + & +wturn3*gshieldc_loc_t3(j,i) + & +wturn4*gshieldc_t4(j,i) + & +wturn4*gshieldc_loc_t4(j,i) + & +wel_loc*gshieldc_ll(j,i) + & +wel_loc*gshieldc_loc_ll(j,i) + gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*fact(2)*gsccorx(j,i) & +wliptran*gliptranx(j,i) + & +welec*gshieldx(j,i) + & +wcorr*gshieldx_ec(j,i) + & +wturn3*gshieldx_t3(j,i) + & +wturn4*gshieldx_t4(j,i) + & +wel_loc*gshieldx_ll(j,i) + else gradc(j,i,icg)=fact(1)*wsc*gvdwc(j,i) & +fact(1)*wscp*gvdwc_scp(j,i)+ @@ -312,11 +329,28 @@ C & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) & +wliptran*gliptranc(j,i) + & +welec*gshieldc(j,i) + & +welec*gshieldc_loc(j,i) + & +wcorr*gshieldc_ec(j,i) + & +wcorr*gshieldc_loc_ec(j,i) + & +wturn3*gshieldc_t3(j,i) + & +wturn3*gshieldc_loc_t3(j,i) + & +wturn4*gshieldc_t4(j,i) + & +wturn4*gshieldc_loc_t4(j,i) + & +wel_loc*gshieldc_ll(j,i) + & +wel_loc*gshieldc_loc_ll(j,i) + gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*fact(1)*gsccorx(j,i) & +wliptran*gliptranx(j,i) + & +welec*gshieldx(j,i) + & +wcorr*gshieldx_ec(j,i) + & +wturn3*gshieldx_t3(j,i) + & +wturn4*gshieldx_t4(j,i) + & +wel_loc*gshieldx_ll(j,i) + else gradc(j,i,icg)=fact(1)*wsc*gvdwc(j,i)+ & fact(1)*wscp*gvdwc_scp(j,i)+ @@ -2171,6 +2205,31 @@ C endif if (ymedi.lt.0) ymedi=ymedi+boxysize zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize + zmedi2=mod(zmedi,boxzsize) + if (zmedi2.lt.0) zmedi2=zmedi2+boxzsize + if ((zmedi2.gt.bordlipbot) + &.and.(zmedi2.lt.bordliptop)) then +C the energy transfer exist + if (zmedi2.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((zmedi2-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipi=sscalelip(fracinbuf) + ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick + elseif (zmedi2.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zmedi2)/lipbufthick) + sslipi=sscalelip(fracinbuf) + ssgradlipi=sscagradlip(fracinbuf)/lipbufthick + else + sslipi=1.0d0 + ssgradlipi=0.0d0 + endif + else + sslipi=0.0d0 + ssgradlipi=0.0d0 + endif + num_conti=0 C write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) do j=ielstart(i),ielend(i) @@ -2216,6 +2275,28 @@ C End diagnostics if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize + if ((zj.gt.bordlipbot) + &.and.(zj.lt.bordliptop)) then +C the energy transfer exist + if (zj.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((zj-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipj=sscalelip(fracinbuf) + ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick + elseif (zj.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) + sslipj=sscalelip(fracinbuf) + ssgradlipj=sscagradlip(fracinbuf)/lipbufthick + else + sslipj=1.0d0 + ssgradlipj=0.0 + endif + else + sslipj=0.0d0 + ssgradlipj=0.0 + endif dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 xj_safe=xj yj_safe=yj @@ -2722,6 +2803,7 @@ C fac_shield(j)=0.6 endif eel_loc_ij=eel_loc_ij & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) c write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij C write (iout,'(a6,2i5,0pf7.3)') C & 'eelloc',i,j,eel_loc_ij @@ -2779,11 +2861,13 @@ C Partial derivatives in virtual-bond dihedral angles gamma & (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)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) 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)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) cd call checkint3(i,j,mu1,mu2,a22,a23,a32,a33,acipa,eel_loc_ij) cd write(iout,*) 'agg ',agg @@ -2797,6 +2881,7 @@ C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) ggg(l)=(agg(l,1)*muij(1)+ & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) enddo do k=i+2,j2 @@ -2809,18 +2894,22 @@ C Remaining derivatives of eello 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)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) 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)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) 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)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) 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)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) enddo endif @@ -3083,6 +3172,37 @@ C Third- and fourth-order contributions from turns double precision agg(3,4),aggi(3,4),aggi1(3,4), & aggj(3,4),aggj1(3,4),a_temp(2,2) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2 + zj=(c(3,j)+c(3,j+1))/2.0d0 +C xj=mod(xj,boxxsize) +C if (xj.lt.0) xj=xj+boxxsize +C yj=mod(yj,boxysize) +C if (yj.lt.0) yj=yj+boxysize + zj=mod(zj,boxzsize) + if (zj.lt.0) zj=zj+boxzsize +C if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ" + if ((zj.gt.bordlipbot) + &.and.(zj.lt.bordliptop)) then +C the energy transfer exist + if (zj.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((zj-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipj=sscalelip(fracinbuf) + ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick + elseif (zj.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) + sslipj=sscalelip(fracinbuf) + ssgradlipj=sscagradlip(fracinbuf)/lipbufthick + else + sslipj=1.0d0 + ssgradlipj=0.0 + endif + else + sslipj=0.0d0 + ssgradlipj=0.0 + endif + if (j.eq.i+2) then if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds @@ -3120,8 +3240,11 @@ C fac_shield(j)=0.6 eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + eello_t3=0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) cd write (2,*) 'i,',i,' j',j,'eello_turn3', cd & 0.5d0*(pizda(1,1)+pizda(2,2)), @@ -3176,6 +3299,8 @@ C Derivatives in gamma(i) call matmat2(a_temp(1,1),pizda(1,1),pizda(1,1)) gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + C Derivatives in gamma(i+1) call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1)) call transpose2(auxmat2(1,1),pizda(1,1)) @@ -3183,6 +3308,7 @@ C Derivatives in gamma(i+1) gel_loc_turn3(i+1)=gel_loc_turn3(i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) C Cartesian derivatives do l=1,3 @@ -3194,6 +3320,7 @@ C Cartesian derivatives gcorr3_turn(l,i)=gcorr3_turn(l,i) & +0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) a_temp(1,1)=aggi1(l,1) a_temp(1,2)=aggi1(l,2) @@ -3203,6 +3330,7 @@ C Cartesian derivatives gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) a_temp(1,1)=aggj(l,1) a_temp(1,2)=aggj(l,2) @@ -3212,6 +3340,7 @@ C Cartesian derivatives gcorr3_turn(l,j)=gcorr3_turn(l,j) & +0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) @@ -3221,6 +3350,7 @@ C Cartesian derivatives gcorr3_turn(l,j1)=gcorr3_turn(l,j1) & +0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) enddo endif @@ -3330,6 +3460,7 @@ C & *2.0 s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) C Derivatives in gamma(i+1) call transpose2(EUgder(1,1,i+2),e2tder(1,1)) @@ -3340,6 +3471,7 @@ C Derivatives in gamma(i+1) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) C Derivatives in gamma(i+2) call transpose2(EUgder(1,1,i+3),e3tder(1,1)) @@ -3353,6 +3485,7 @@ C Derivatives in gamma(i+2) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) C Cartesian derivatives @@ -3375,6 +3508,7 @@ C Derivatives of this turn contributions in DC(i+2) ggg(l)=-(s1+s2+s3) gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) enddo endif @@ -3395,6 +3529,7 @@ C Remaining derivatives of this turn contribution s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) a_temp(1,1)=aggi1(l,1) a_temp(1,2)=aggi1(l,2) @@ -3411,6 +3546,7 @@ C Remaining derivatives of this turn contribution s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) a_temp(1,1)=aggj(l,1) a_temp(1,2)=aggj(l,2) @@ -3427,6 +3563,7 @@ C Remaining derivatives of this turn contribution s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) @@ -3443,8 +3580,17 @@ C Remaining derivatives of this turn contribution s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) enddo + gshieldc_t4(3,i)=gshieldc_t4(3,i)+ + & ssgradlipi*eello_t4/4.0d0*lipscale + gshieldc_t4(3,j)=gshieldc_t4(3,j)+ + & ssgradlipj*eello_t4/4.0d0*lipscale + gshieldc_t4(3,i-1)=gshieldc_t4(3,i-1)+ + & ssgradlipi*eello_t4/4.0d0*lipscale + gshieldc_t4(3,j-1)=gshieldc_t4(3,j-1)+ + & ssgradlipj*eello_t4/4.0d0*lipscale endif 178 continue endif