From 688de817003e3e5f38e18eee73f13c25419c7f00 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Mon, 9 May 2016 18:23:55 +0200 Subject: [PATCH] bugfix in shield FGPROC>1 --- source/unres/src_MD-M/COMMON.SHIELD | 8 +- source/unres/src_MD-M/energy_p_new-sep_barrier.F | 14 +- source/unres/src_MD-M/energy_p_new_barrier.F | 294 +++++++++++++++++++++- source/unres/src_MD-M/energy_split-sep.F | 4 +- source/unres/src_MD-M/readrtns_CSA.F | 6 + 5 files changed, 307 insertions(+), 19 deletions(-) diff --git a/source/unres/src_MD-M/COMMON.SHIELD b/source/unres/src_MD-M/COMMON.SHIELD index 6cbf5a8..2ead23e 100644 --- a/source/unres/src_MD-M/COMMON.SHIELD +++ b/source/unres/src_MD-M/COMMON.SHIELD @@ -1,13 +1,15 @@ double precision VSolvSphere,VSolvSphere_div,long_r_sidechain, & short_r_sidechain,fac_shield,grad_shield_side,grad_shield, - & buff_shield,wshield + & buff_shield,wshield,lipscale,sslipi,sslipj, + & ssgradlipi,ssgradlipj integer ishield_list,shield_list,ees0plist common /shield/ VSolvSphere,VSolvSphere_div,buff_shield, - & long_r_sidechain(ntyp), + & long_r_sidechain(ntyp),sslipi,sslipj, + & ssgradlipi,ssgradlipj, & short_r_sidechain(ntyp),fac_shield(maxres),wshield, & grad_shield_side(3,maxcontsshi,-1:maxres), & grad_shield(3,-1:maxres), - & grad_shield_loc(3,maxcontsshi,-1:maxres), + & grad_shield_loc(3,maxcontsshi,-1:maxres),lipscale, & ishield_list(maxres),shield_list(maxcontsshi,maxres), & ees0plist(maxconts,maxres) 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 0a638ef..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 @@ -2683,7 +2683,6 @@ 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))) sssgrad=sscagrad(1.0d0/(dsqrt(rrij))) if (sss.lt.1.0d0) then @@ -2704,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) + fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/expon ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac @@ -2795,7 +2794,7 @@ C Uncomment following three lines for Ca-p interactions if (xj.lt.0) xj=xj+boxxsize yj=mod(yj,boxysize) if (yj.lt.0) yj=yj+boxysize - zj=mod(zi,boxzsize) + 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 @@ -2831,7 +2830,7 @@ C Uncomment following three lines for Ca-p interactions rrij=1.0D0/(xj*xj+yj*yj+zj*zj) sss=sscale(1.0d0/(dsqrt(rrij))) sssgrad=sscagrad(1.0d0/(dsqrt(rrij))) - if (sss.gt.0.0d0) then + if (sss.gt.0.0d0) then fac=rrij**expon2 e1=fac*fac*aad(itypj,iteli) @@ -2845,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) + fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/expon ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 03bc765..c0670a7 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -3421,6 +3421,7 @@ C include 'COMMON.FFIELD' include 'COMMON.TIME1' include 'COMMON.SPLITELE' + include 'COMMON.SHIELD' dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), @@ -3537,6 +3538,30 @@ c end if 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 call eelecij(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) @@ -3589,13 +3614,35 @@ c if ((zmedi.gt.((0.5d0)*boxzsize)).or. c & (zmedi.lt.((-0.5d0)*boxzsize))) then c go to 196 c endif - xmedi=mod(xmedi,boxxsize) if (xmedi.lt.0) xmedi=xmedi+boxxsize ymedi=mod(ymedi,boxysize) 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.0 + endif + else + sslipi=0.0d0 + ssgradlipi=0.0 + endif num_conti=num_cont_hb(i) c write(iout,*) "JESTEM W PETLI" call eelecij(i,i+3,ees,evdw1,eel_loc) @@ -3637,6 +3684,29 @@ c & .or. itype(i-1).eq.ntyp1 if (ymedi.lt.0) ymedi=ymedi+boxysize zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize + if ((zmedi.gt.bordlipbot) + &.and.(zmedi.lt.bordliptop)) then +C the energy transfer exist + if (zmedi.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((zmedi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipi=sscalelip(fracinbuf) + ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick + elseif (zmedi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zmedi)/lipbufthick) + sslipi=sscalelip(fracinbuf) + ssgradlipi=sscagradlip(fracinbuf)/lipbufthick + else + sslipi=1.0d0 + ssgradlipi=0.0 + endif + else + sslipi=0.0d0 + ssgradlipi=0.0 + endif +C print *,sslipi,"TU?!" C xmedi=xmedi+xshift*boxxsize C ymedi=ymedi+yshift*boxysize C zmedi=zmedi+zshift*boxzsize @@ -3771,6 +3841,28 @@ C zj=c(3,j)+0.5D0*dzj-zmedi zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize 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 dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 xj_safe=xj yj_safe=yj @@ -3865,13 +3957,20 @@ C fac_shield(j)=0.6 el2=el2*fac_shield(i)**2*fac_shield(j)**2 eesij=(el1+el2) ees=ees+eesij +C FOR NOW SHIELD IS NOT USED WITH LIPSCALE +C & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) else fac_shield(i)=1.0 fac_shield(j)=1.0 eesij=(el1+el2) ees=ees+eesij + &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) +C print *,"TUCC",(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0 endif evdw1=evdw1+evdwij*sss + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) +C print *,sslipi,sslipj,lipscale**2, +C & (sslipi+sslipj)/2.0d0*lipscale**2+1.0d0 cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, cd & 1.0D0/dsqrt(rrmij),evdwij,eesij, @@ -3891,7 +3990,9 @@ C Calculate contributions to the Cartesian gradient. C #ifdef SPLITELE facvdw=-6*rrmij*(ev1+evdwij)*sss + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) facel=-3*rrmij*(el1+eesij) + &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) fac1=fac erij(1)=xj*rmij erij(2)=yj*rmij @@ -3990,6 +4091,12 @@ C gelc_long(k,j-1)=gelc_long(k,j-1) C & +grad_shield(k,j)*eesij/fac_shield(j) enddo C print *,"bafter", gelc_long(1,i), gelc_long(1,j) +C Lipidic part for lipscale + gelc_long(3,j)=gelc_long(3,j)+ + & ssgradlipj*eesij/2.0d0*lipscale**2 + + gelc_long(3,i)=gelc_long(3,i)+ + & ssgradlipi*eesij/2.0d0*lipscale**2 * * Loop over residues i+1 thru j-1. @@ -4001,8 +4108,13 @@ cgrad enddo cgrad enddo if (sss.gt.0.0) then ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) else ggg(1)=0.0 ggg(2)=0.0 @@ -4018,6 +4130,12 @@ c 9/28/08 AL Gradient compotents will be summed only at the end gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo +C Lipidic part for scaling weight + gvdwpp(3,j)=gvdwpp(3,j)+ + & sss*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ + & sss*ssgradlipi*evdwij/2.0d0*lipscale**2 + * * Loop over residues i+1 thru j-1. * @@ -4029,6 +4147,7 @@ cgrad enddo #else C MARYSIA facvdw=(ev1+evdwij)*sss + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) facel=(el1+eesij) fac1=fac fac=-3*rrmij*(facvdw+facvdw+facel) @@ -4064,12 +4183,22 @@ cgrad enddo cgrad enddo c 9/28/08 AL Gradient compotents will be summed only at the end ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) + ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) do k=1,3 gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) enddo + gvdwpp(3,j)=gvdwpp(3,j)+ + & sss*ssgradlipj*evdwij/2.0d0*lipscale**2 + gvdwpp(3,i)=gvdwpp(3,i)+ + & sss*ssgradlipi*evdwij/2.0d0*lipscale**2 + #endif * * Angular part @@ -4088,6 +4217,7 @@ cd & (dcosg(k),k=1,3) do k=1,3 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))* & fac_shield(i)**2*fac_shield(j)**2 + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) enddo c do k=1,3 c ghalf=0.5D0*ggg(k) @@ -4109,10 +4239,12 @@ C print *,"before22", gelc_long(1,i), gelc_long(1,j) & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)) & *fac_shield(i)**2*fac_shield(j)**2 + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) 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)) & *fac_shield(i)**2*fac_shield(j)**2 + & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0) gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo @@ -4327,6 +4459,8 @@ 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 Now derivative over eel_loc if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & (shield_mode.gt.0)) then @@ -4383,6 +4517,8 @@ C Calculate patrial derivative for theta angle & +a32*gmuij1(3) & +a33*gmuij1(4)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + c write(iout,*) "derivative over thatai" c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), c & a33*gmuij1(4) @@ -4399,6 +4535,8 @@ c & a33*gmuij2(4) gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & geel_loc_ij*wel_loc & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + c Derivative over j residue geel_loc_ji=a22*gmuji1(1) @@ -4412,6 +4550,7 @@ c & a33*gmuji1(4) gloc(nphi+j,icg)=gloc(nphi+j,icg)+ & geel_loc_ji*wel_loc & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) geel_loc_ji= & +a22*gmuji2(1) @@ -4424,6 +4563,8 @@ c & a33*gmuji2(4) gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+ & geel_loc_ji*wel_loc & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + #endif cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij @@ -4440,22 +4581,35 @@ 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) + C 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)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l) gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l) cgrad ghalf=0.5d0*ggg(l) cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf enddo + gel_loc_long(3,j)=gel_loc_long(3,j)+ + & ssgradlipj*eel_loc_ij/2.0d0*lipscale/ + & ((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + + gel_loc_long(3,i)=gel_loc_long(3,i)+ + & ssgradlipi*eel_loc_ij/2.0d0*lipscale/ + & ((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + cgrad do k=i+1,j2 cgrad do l=1,3 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l) @@ -4466,18 +4620,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 @@ -4723,7 +4881,42 @@ C Third- and fourth-order contributions from turns & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 j=i+2 -c write (iout,*) "eturn3",i,j,j1,j2 +C xj=(c(1,j)+c(1,j+1))/2.0d0 +C yj=(c(2,j)+c(2,j+1))/2.0d0 + 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 + 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 +C sslipj=0.0 +C ssgradlipj=0.0d0 + +C write (iout,*) "eturn3",i,j,j1,j2 a_temp(1,1)=a22 a_temp(1,2)=a23 a_temp(2,1)=a32 @@ -4752,26 +4945,34 @@ c auxalary matrix for i+2 and constant i+1 call matmat2(a_temp(1,1),auxgmatt1(1,1),gpizda1(1,1)) call matmat2(a_temp(1,1),auxgmatt2(1,1),gpizda2(1,1)) if (shield_mode.eq.0) then - fac_shield(i)=1.0 - fac_shield(j)=1.0 + fac_shield(i)=1.0d0 + fac_shield(j)=1.0d0 C else C fac_shield(i)=0.4 C fac_shield(j)=0.6 endif C if (j.eq.78) C & write(iout,*) i,j,fac_shield(i),fac_shield(j) - eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2)) + eello_turn3=eello_turn3+ +C & 1.0*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + &0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) - eello_t3=0.5d0*(pizda(1,1)+pizda(2,2)) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + eello_t3= + &0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) #ifdef NEWCORR C Derivatives in theta gloc(nphi+i,icg)=gloc(nphi+i,icg) & +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3 & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg) & +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3 & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + #endif C if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') @@ -4828,6 +5029,8 @@ C Derivatives in gamma(i) call matmat2(a_temp(1,1),auxmat3(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),auxmat3(1,1)) @@ -4835,6 +5038,8 @@ 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 c ghalf1=0.5d0*agg(l,1) @@ -4849,6 +5054,7 @@ c ghalf4=0.5d0*agg(l,4) 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)!+agg(l,1) a_temp(1,2)=aggi1(l,2)!+agg(l,2) @@ -4858,6 +5064,7 @@ c ghalf4=0.5d0*agg(l,4) 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)!+ghalf1 a_temp(1,2)=aggj(l,2)!+ghalf2 a_temp(2,1)=aggj(l,3)!+ghalf3 @@ -4866,6 +5073,8 @@ c ghalf4=0.5d0*agg(l,4) 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) a_temp(2,1)=aggj1(l,3) @@ -4874,7 +5083,18 @@ c ghalf4=0.5d0*agg(l,4) 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 + gshieldc_t3(3,i)=gshieldc_t3(3,i)+ + & ssgradlipi*eello_t3/4.0d0*lipscale + gshieldc_t3(3,j)=gshieldc_t3(3,j)+ + & ssgradlipj*eello_t3/4.0d0*lipscale + gshieldc_t3(3,i-1)=gshieldc_t3(3,i-1)+ + & ssgradlipi*eello_t3/4.0d0*lipscale + gshieldc_t3(3,j-1)=gshieldc_t3(3,j-1)+ + & ssgradlipj*eello_t3/4.0d0*lipscale + +C print *,ssgradlipi,ssgradlipj,eello_t3,sslipi,sslipj return end C------------------------------------------------------------------------------- @@ -4923,6 +5143,37 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC cd call checkint_turn4(i,a_temp,eello_turn4_num) c write (iout,*) "eturn4 i",i," j",j," j1",j1," j2",j2 c write(iout,*)"WCHODZE W PROGRAM" + 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 + a_temp(1,1)=a22 a_temp(1,2)=a23 a_temp(2,1)=a32 @@ -5004,6 +5255,8 @@ C fac_shield(j)=0.4 endif eello_turn4=eello_turn4-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + eello_t4=-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) c write(iout,*)'chujOWO', auxvec(1),b1(1,iti2) @@ -5063,13 +5316,17 @@ cd & ' eello_turn4_num',8*eello_turn4_num gloc(nphi+i,icg)=gloc(nphi+i,icg) & -(gs13+gsE13+gsEE1)*wturn4 & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg) & -(gs23+gs21+gsEE2)*wturn4 & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg) & -(gs32+gsE31+gsEE3)*wturn4 & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) c gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)- c & gs2 @@ -5087,6 +5344,8 @@ C Derivatives in gamma(i) 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)) call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1)) @@ -5096,6 +5355,8 @@ 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)) call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1)) @@ -5108,6 +5369,8 @@ 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 C Derivatives of this turn contributions in DC(i+2) if (j.lt.nres-1) then @@ -5128,6 +5391,8 @@ 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 C Remaining derivatives of this turn contribution @@ -5147,6 +5412,8 @@ 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) a_temp(2,1)=aggi1(l,3) @@ -5162,6 +5429,8 @@ 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) a_temp(2,1)=aggj(l,3) @@ -5177,6 +5446,8 @@ 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) a_temp(2,1)=aggj1(l,3) @@ -5193,7 +5464,16 @@ C Remaining derivatives of this turn contribution c write (iout,*) "s1",s1," s2",s2," s3",s3," s1+s2+s3",s1+s2+s3 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 return end C----------------------------------------------------------------------------- diff --git a/source/unres/src_MD-M/energy_split-sep.F b/source/unres/src_MD-M/energy_split-sep.F index b3f924b..f79deea 100644 --- a/source/unres/src_MD-M/energy_split-sep.F +++ b/source/unres/src_MD-M/energy_split-sep.F @@ -293,7 +293,7 @@ C call sum_energy(energia,.true.) C call enerprint C write (iout,*) "Exit ETOTAL_LONG" - call flush(iout) +C call flush(iout) return end c------------------------------------------------------------------------------ @@ -627,6 +627,6 @@ C write (iout,*) "ETOTAL_SHORT before SUM_ENERGY" call flush(iout) call sum_energy(energia,.true.) C write (iout,*) "Exit ETOTAL_SHORT" - call flush(iout) +C call flush(iout) return end diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index e40717f..d0c10c7 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -766,10 +766,16 @@ C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,"BTRISS",btriss,0.021D0) call reada(weightcard,"CTRISS",ctriss,1.001D0) call reada(weightcard,"DTRISS",dtriss,1.001D0) + call reada(weightcard,"LIPSCALE",lipscale,1.3D0) write (iout,*) "ATRISS=", atriss write (iout,*) "BTRISS=", btriss write (iout,*) "CTRISS=", ctriss write (iout,*) "DTRISS=", dtriss + if (shield_mode.gt.0) then + lipscale=0.0d0 + write (iout,*) "liscale not used in shield mode" + endif + write (iout,*) "lipid scaling factor", lipscale dyn_ss=(index(weightcard,'DYN_SS').gt.0) do i=1,maxres dyn_ss_mask(i)=.false. -- 1.7.9.5