implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
include "DIMENSIONS.COMPAR"
+ include 'COMMON.CONTROL'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.LOCAL'
#endif
C#undef DEBUG
c endif
+ 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
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
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
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/
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)
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
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
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/
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)
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
* 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
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
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
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)
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
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)
c & write (iout,'(a4,2i4,8f9.5)')'chuj',
c & i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4)
- eel_loc=eel_loc+eel_loc_ij*sss
+ eel_loc=eel_loc+eel_loc_ij
C Now derivative over eel_loc
if (calc_grad) then
if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
& +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)
& +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)
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)
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
& 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
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 enddo
cgrad enddo
C Remaining derivatives of eello
+ 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
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
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
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,
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
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)
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)
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
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)
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
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)
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))
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))
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)
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
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)
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)
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)
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
projects / unres.git / blobdiff