#ifdef MPL
c endif
#endif
+#define DEBUG
+#ifdef DEBUG
+ call enerprint(energia,fact)
+#endif
+#undef DEBUG
if (calc_grad) then
C
C Sum up the components of the Cartesian gradient.
& 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)+
& 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)=fact(1)*wsc*gvdwx(j,i)
& +fact(1)*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)
+
endif
enddo
& 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)+
& 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)=fact(1)*wsc*gvdwx(j,i)+
& fact(1)*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)
+
endif
enddo
#endif
& +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)) then
+
+C write(iout,*) "tu,", i,j,aa_aq(itypi,itypj)-aa,
+C & bb_aq(itypi,itypj)-bb,
+C & sslipi,sslipj
+C endif
C write(iout,*),aa,aa_lip(itypi,itypj),aa_aq(itypi,itypj)
C checking the distance
c if (lprn) then
sigm=dabs(aa/bb)**(1.0D0/6.0D0)
epsi=bb**2/aa
+C#define DEBUG
#ifdef DEBUG
write (iout,'(2(a3,i3,2x),17(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
& evdwij
write (iout,*) "partial sum", evdw, evdw_t
#endif
+C#undef DEBUG
c endif
if (calc_grad) then
C Calculate gradient components.
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)
- if (j.le.1) cycle
+ if (j.lt.1) cycle
C if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
C & .or.itype(j+2).eq.ntyp1
C &) cycle
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
C 12/26/95 - for the evaluation of multi-body H-bonding interactions
ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
if (shield_mode.gt.0) then
+C#define DEBUG
+#ifdef DEBUG
+ write(iout,*) "ees_compon",i,j,el1,el2,
+ & fac_shield(i),fac_shield(j)
+#endif
+C#undef DEBUG
C fac_shield(i)=0.4
C fac_shield(j)=0.6
el1=el1*fac_shield(i)**2*fac_shield(j)**2
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
& (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
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
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
& *fac_shield(i)*fac_shield(j)
gacontp_hb3(k,num_conti,i)=gggp(k)
+ & *fac_shield(i)*fac_shield(j)
+
gacontm_hb1(k,num_conti,i)=ghalfm
& +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
& + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
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
C & .or. itype(i+5).eq.ntyp1
C & .or. itype(i).eq.ntyp1
C & .or. itype(i-1).eq.ntyp1
- & ) goto 178
+ & ) goto 179
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
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)),
call transpose2(auxmat2(1,1),pizda(1,1))
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))
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
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)
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)
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)
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
+ 179 continue
else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1) then
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
C changes suggested by Ana to avoid out of bounds
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))
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))
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
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
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)
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)
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)
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
projects / unres.git / blobdiff