X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-HCD%2Fenergy_p_new.F;h=e72d5584c13047539fce853ee48eb7aa5bbb7722;hb=f9783c27de185ff3e1d6e874a37a1a8c4d3c93a3;hp=ce7a6a7ad1cdb402f08724cc05283149fecbd3c7;hpb=3d04f1babd4ae90f9a0b8f256a2ce0a2c5751f33;p=unres.git diff --git a/source/wham/src-HCD/energy_p_new.F b/source/wham/src-HCD/energy_p_new.F index ce7a6a7..e72d558 100644 --- a/source/wham/src-HCD/energy_p_new.F +++ b/source/wham/src-HCD/energy_p_new.F @@ -1216,6 +1216,9 @@ C returning jth atom to box & +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 +c write(iout,*) "tu,", i,j,aa_lip(itypi,itypj),bb_lip(itypi,itypj) +c if (aa.ne.aa_aq(itypi,itypj)) write(iout,'(2e15.5)') +c &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj)) xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) @@ -1278,8 +1281,8 @@ c#define DEBUG #endif c#undef DEBUG c endif - if (energy_dec) write (iout,'(a,2i5,3f10.5)') - & 'r sss evdw',i,j,1.0d0/rij,sss,evdwij + if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)') + & 'r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij if (calc_grad) then C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 @@ -1288,6 +1291,12 @@ C Calculate gradient components. fac=rij*fac fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij C Calculate the radial part of the gradient + 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 gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac @@ -1457,6 +1466,12 @@ C Calculate gradient components. fac=rij*fac-2*expon*rrij*e_augm fac=fac+(evdwij+e_augm)*sssgrad/sss*rij C Calculate the radial part of the gradient + 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 gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac @@ -2110,6 +2125,8 @@ C common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + double precision sslipi,sslipj,ssgradlipi,ssgradlipj + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ @@ -2217,6 +2234,7 @@ c end if ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) num_conti=0 call eelecij(i,i+2,ees,evdw1,eel_loc) if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) @@ -2247,6 +2265,7 @@ c & .or. itype(i-1).eq.ntyp1 ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) #ifdef FOURBODY num_conti=num_cont_hb(i) #endif @@ -2287,6 +2306,7 @@ c & .or. itype(i-1).eq.ntyp1 ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi call to_box(xmedi,ymedi,zmedi) + call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi) #ifdef FOURBODY num_conti=num_cont_hb(i) #endif @@ -2358,6 +2378,9 @@ C------------------------------------------------------------------------------- common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij, + & faclipij2 + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ @@ -2393,6 +2416,9 @@ C zj=c(3,j)+0.5D0*dzj-zmedi yj=c(2,j)+0.5D0*dyj zj=c(3,j)+0.5D0*dzj call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) + faclipij=(sslipi+sslipj)/2.0d0*lipscale+1.0d0 + faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0 xj=boxshift(xj-xmedi,boxxsize) yj=boxshift(yj-ymedi,boxysize) zj=boxshift(zj-zmedi,boxzsize) @@ -2432,25 +2458,25 @@ C fac_shield(j)=0.6 el1=el1*fac_shield(i)**2*fac_shield(j)**2 el2=el2*fac_shield(i)**2*fac_shield(j)**2 eesij=(el1+el2) - ees=ees+eesij + ees=ees+eesij*faclipij2 else fac_shield(i)=1.0 fac_shield(j)=1.0 eesij=(el1+el2) - ees=ees+eesij + ees=ees+eesij*faclipij2 endif - evdw1=evdw1+evdwij*sss + evdw1=evdw1+evdwij*sss*faclipij2 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, cd & xmedi,ymedi,zmedi,xj,yj,zj if (energy_dec) then - write (iout,'(a6,2i5,0pf7.3,2i5,3e11.3)') - &'evdw1',i,j,evdwij - &,iteli,itelj,aaa,evdw1,sss - write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij, - &fac_shield(i),fac_shield(j) + write (iout,'(a6,2i5,0pf7.3,2i5,3e11.3)') + &' evdw1',i,j,evdwij,iteli,itelj,aaa,evdw1,sss + write (iout,'(a6,2i5,0pf7.3,6f8.5)') 'ees',i,j,eesij, + & fac_shield(i),fac_shield(j),sslipi,sslipj,faclipij, + & faclipij2 endif C @@ -2468,9 +2494,10 @@ C * Radial derivatives. First process both termini of the fragment (i,j) * if (calc_grad) then - ggg(1)=facel*xj - ggg(2)=facel*yj - ggg(3)=facel*zj + aux=(facel*sss+rmij*sssgrad*eesij)*faclipij2 + ggg(1)=aux*xj + ggg(2)=aux*yj + ggg(3)=aux*zj if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & (shield_mode.gt.0)) then C print *,i,j @@ -2556,6 +2583,11 @@ C gelc_long(k,i-1)=gelc_long(k,i-1) C & +grad_shield(k,i)*eesij/fac_shield(i) C gelc_long(k,j-1)=gelc_long(k,j-1) C & +grad_shield(k,j)*eesij/fac_shield(j) + gelc_long(3,j)=gelc_long(3,j)+ + & ssgradlipj*eesij/2.0d0*lipscale**2*sss + + gelc_long(3,i)=gelc_long(3,i)+ + & ssgradlipi*eesij/2.0d0*lipscale**2*sss enddo C print *,"bafter", gelc_long(1,i), gelc_long(1,j) @@ -2568,7 +2600,7 @@ cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo if (sss.gt.0.0) then - facvdw=facvdw+sssgrad*rmij*evdwij + facvdw=(facvdw+sssgrad*rmij*evdwij)*faclipij2 ggg(1)=facvdw*xj ggg(2)=facvdw*yj ggg(3)=facvdw*zj @@ -2587,6 +2619,11 @@ 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. * @@ -2598,7 +2635,7 @@ cgrad enddo endif ! calc_grad #else C MARYSIA - facvdw=(ev1+evdwij) + facvdw=(ev1+evdwij)*faclipij2 facel=(el1+eesij) fac1=fac fac=-3*rrmij*(facvdw+facvdw+facel)*sss @@ -2642,6 +2679,10 @@ 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 + 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 ! calc_grad #endif * @@ -2661,7 +2702,7 @@ cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3), 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 + & fac_shield(i)**2*fac_shield(j)**2*sss*faclipij2 enddo c do k=1,3 c ghalf=0.5D0*ggg(k) @@ -2682,11 +2723,11 @@ C print *,"before22", gelc_long(1,i), gelc_long(1,j) 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)) - & *fac_shield(i)**2*fac_shield(j)**2 + & *fac_shield(i)**2*fac_shield(j)**2*faclipij2 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 + & *fac_shield(i)**2*fac_shield(j)**2*faclipij2 gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo @@ -2918,7 +2959,7 @@ C fac_shield(i)=0.4 C fac_shield(j)=0.6 endif eel_loc_ij=eel_loc_ij - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eelloc',i,j,eel_loc_ij c if (eel_loc_ij.ne.0) @@ -2982,7 +3023,7 @@ C Calculate patrial derivative for theta angle & +a23*gmuij1(2) & +a32*gmuij1(3) & +a33*gmuij1(4)) - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij c write(iout,*) "derivative over thatai" c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), c & a33*gmuij1(4) @@ -2998,7 +3039,7 @@ c & a33*gmuij2(4) & +a33*gmuij2(4) gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & geel_loc_ij*wel_loc - & *fac_shield(i)*fac_shield(j)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij c Derivative over j residue geel_loc_ji=a22*gmuji1(1) @@ -3011,7 +3052,7 @@ c & a33*gmuji1(4) gloc(nphi+j,icg)=gloc(nphi+j,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij geel_loc_ji= & +a22*gmuji2(1) @@ -3023,7 +3064,7 @@ c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), 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)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij #endif cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij @@ -3032,12 +3073,12 @@ C Partial derivatives in virtual-bond dihedral angles gamma & 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)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij 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) + & *fac_shield(i)*fac_shield(j)*sss*faclipij C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) aux=eel_loc_ij/sss*sssgrad*rmij ggg(1)=aux*xj @@ -3046,13 +3087,18 @@ C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) do l=1,3 ggg(l)=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)*sss + & *fac_shield(i)*fac_shield(j)*sss*faclipij 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/faclipij + + gel_loc_long(3,i)=gel_loc_long(3,i)+ + & ssgradlipi*eel_loc_ij/2.0d0*lipscale/faclipij cgrad do k=i+1,j2 cgrad do l=1,3 cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l) @@ -3062,19 +3108,19 @@ C 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)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss*faclipij 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) + & *fac_shield(i)*fac_shield(j)*sss*faclipij 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) + & *fac_shield(i)*fac_shield(j)*sss*faclipij 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) + & *fac_shield(i)*fac_shield(j)*sss*faclipij enddo endif ! calc_grad @@ -3342,6 +3388,8 @@ C Third- and fourth-order contributions from turns common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij + common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij j=i+2 c write (iout,*) "eturn3",i,j,j1,j2 a_temp(1,1)=a22 @@ -3379,7 +3427,7 @@ C fac_shield(i)=0.4 C fac_shield(j)=0.6 endif eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij eello_t3=0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) if (energy_dec) write (iout,'(6heturn3,2i5,0pf7.3)') i,i+2, @@ -3389,10 +3437,10 @@ C#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) + & *fac_shield(i)*fac_shield(j)*faclipij 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) + & *fac_shield(i)*fac_shield(j)*faclipij C#endif C Derivatives in shield mode @@ -3447,14 +3495,14 @@ C Derivatives in gamma(i) call transpose2(auxmat2(1,1),auxmat3(1,1)) 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) + & *fac_shield(i)*fac_shield(j)*faclipij 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)) call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,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) + & *fac_shield(i)*fac_shield(j)*faclipij C Cartesian derivatives do l=1,3 c ghalf1=0.5d0*agg(l,1) @@ -3468,7 +3516,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i)=gcorr3_turn(l,i) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggi1(l,1)!+agg(l,1) a_temp(1,2)=aggi1(l,2)!+agg(l,2) @@ -3477,7 +3525,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj(l,1)!+ghalf1 a_temp(1,2)=aggj(l,2)!+ghalf2 a_temp(2,1)=aggj(l,3)!+ghalf3 @@ -3485,7 +3533,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j)=gcorr3_turn(l,j) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -3493,7 +3541,7 @@ c ghalf4=0.5d0*agg(l,4) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) gcorr3_turn(l,j1)=gcorr3_turn(l,j1) & +0.5d0*(pizda(1,1)+pizda(2,2)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij enddo endif ! calc_grad @@ -3628,7 +3676,7 @@ C fac_shield(i)=0.6 C fac_shield(j)=0.4 endif eello_turn4=eello_turn4-(s1+s2+s3) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*faclipij eello_t4=-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) c write(iout,*)'chujOWO', auxvec(1),b1(1,iti2) @@ -3705,7 +3753,7 @@ C Derivatives in gamma(i) call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1)) 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) + & *fac_shield(i)*fac_shield(j)*faclipij 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)) @@ -3714,7 +3762,7 @@ C Derivatives in gamma(i+1) call matmat2(auxmat(1,1),e1t(1,1),pizda(1,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) + & *fac_shield(i)*fac_shield(j)*faclipij 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)) @@ -3726,7 +3774,7 @@ C Derivatives in gamma(i+2) call matmat2(auxmat3(1,1),e1t(1,1),pizda(1,1)) 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) + & *fac_shield(i)*fac_shield(j)*faclipij if (calc_grad) then C Cartesian derivatives C Derivatives of this turn contributions in DC(i+2) @@ -3747,7 +3795,7 @@ C Derivatives of this turn contributions in DC(i+2) s3=0.5d0*(pizda(1,1)+pizda(2,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) + & *fac_shield(i)*fac_shield(j)*faclipij enddo endif C Remaining derivatives of this turn contribution @@ -3766,7 +3814,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) 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) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggi1(l,1) a_temp(1,2)=aggi1(l,2) a_temp(2,1)=aggi1(l,3) @@ -3781,7 +3829,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) 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) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj(l,1) a_temp(1,2)=aggj(l,2) a_temp(2,1)=aggj(l,3) @@ -3796,7 +3844,7 @@ C Remaining derivatives of this turn contribution call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) 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) + & *fac_shield(i)*fac_shield(j)*faclipij a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) a_temp(2,1)=aggj1(l,3) @@ -3812,7 +3860,7 @@ C Remaining derivatives of this turn contribution s3=0.5d0*(pizda(1,1)+pizda(2,2)) 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) + & *fac_shield(i)*fac_shield(j)*faclipij enddo endif ! calc_grad @@ -9492,6 +9540,8 @@ C--bufliptop--- here true lipid starts C lipid C--buflipbot--- lipid ends buffore starts C--bordlipbot--buffore ends +c call cartprint +c write (iout,*) "Eliptransfer peplipran",pepliptran eliptran=0.0 do i=1,nres C do i=1,1 @@ -9543,6 +9593,8 @@ CV do i=1,1 if (itype(i).eq.ntyp1) cycle positi=(mod(c(3,i+nres),boxzsize)) if (positi.le.0) positi=positi+boxzsize +c write(iout,*) "i",i," positi",positi,bordlipbot,buflipbot, +c & bordliptop C print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop c for each residue check if it is in lipid or lipid water border area C respos=mod(c(3,i+nres),boxzsize) @@ -9553,8 +9605,11 @@ C the energy transfer exist if (positi.lt.buflipbot) then fracinbuf=1.0d0- & ((positi-bordlipbot)/lipbufthick) +c write (iout,*) "i",i,itype(i)," fracinbuf",fracinbuf +c write (iout,*) "i",i," liptranene",liptranene(itype(i)) C lipbufthick is thickenes of lipid buffore sslip=sscalelip(fracinbuf) +c write (iout,*) "sslip",sslip ssgradlip=-sscagradlip(fracinbuf)/lipbufthick eliptran=eliptran+sslip*liptranene(itype(i)) gliptranx(3,i)=gliptranx(3,i) @@ -9580,6 +9635,7 @@ C print *,"I am in true lipid" endif ! if in lipid or buffor C else C eliptran=elpitran+0.0 ! I am in water +c write (iout,*) "eliptran",eliptran enddo return end