endif
c print *,"Processor",myrank," computed Utord"
C
- call eback_sc_corr(esccor)
+ if (wsccor.gt.0.0d0) then
+ call eback_sc_corr(esccor)
+ else
+ esccor=0.0d0
+ endif
if (wliptran.gt.0) then
call Eliptransfer(eliptran)
+ else
+ eliptran=0.0d0
endif
-
+#ifdef FOURBODY
C
C 12/1/95 Multi-body terms
C
c write (iout,*) "Calling multibody_hbond"
call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
endif
+#endif
c write (iout,*) "NSAXS",nsaxs
if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
call e_saxs(Esaxs_constr)
c write(iout,*) "TEST_ENE1 ehomology_constr=",ehomology_constr
#ifdef DFA
C BARTEK for dfa test!
+ edfadis=0.0d0
if (wdfa_dist.gt.0) call edfad(edfadis)
c write(iout,*)'edfad is finished!', wdfa_dist,edfadis
+ edfator=0.0d0
if (wdfa_tor.gt.0) call edfat(edfator)
c write(iout,*)'edfat is finished!', wdfa_tor,edfator
+ edfanei=0.0d0
if (wdfa_nei.gt.0) call edfan(edfanei)
c write(iout,*)'edfan is finished!', wdfa_nei,edfanei
+ edfabet=0.0d0
if (wdfa_beta.gt.0) call edfab(edfabet)
c write(iout,*)'edfab is finished!', wdfa_beta,edfabet
+#else
+ edfadis=0.0d0
+ edfator=0.0d0
+ edfanei=0.0d0
+ edfabet=0.0d0
#endif
#ifdef SPLITELE
& +welec*fact(1)*ees
& +fact(1)*wvdwpp*evdw1
& +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
- & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
- & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
+ & +wstrain*ehpb
+ & +wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+ & +wcorr6*fact(5)*ecorr6
+ & +wturn4*fact(3)*eello_turn4
& +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
& +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
edfator = energia(29)
edfanei = energia(30)
edfabet = energia(31)
+ Eafmforc=0.0d0
+ etube=0.0d0
+ Uconst=0.0d0
#ifdef SPLITELE
write(iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1,wvdwpp,
& estr,wbond,ebe,wang,escloc,wscloc,etors,wtor*fact(1),
- & etors_d,wtor_d*fact(2),ehpb,wstrain,ecorr,wcorr*fact(3),
- & ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5),eel_loc,
+ & etors_d,wtor_d*fact(2),ehpb,wstrain,
+#ifdef FOURBODY
+ & ecorr,wcorr*fact(3),
+ & ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5),
+#endif
+ & eel_loc,
& wel_loc*fact(2),eello_turn3,wturn3*fact(2),
- & eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5),
+ & eello_turn4,wturn4*fact(3),
+#ifdef FOURBODY
+ & eello_turn6,wturn6*fact(5),
+#endif
& esccor,wsccor*fact(1),edihcnstr,
& ethetacnstr,ebr*nss,Uconst,wumb,eliptran,wliptran,Eafmforc,
& etube,wtube,esaxs,wsaxs,ehomology_constr,
& 'ETORSD=',1pE16.6,' WEIGHT=',1pE16.6,' (double torsional)'/
& 'EHBP= ',1pE16.6,' WEIGHT=',1pE16.6,
& ' (SS bridges & dist. cnstr.)'/
+#ifdef FOURBODY
& 'ECORR4=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
& 'ECORR5=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
& 'ECORR6=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
+#endif
& 'EELLO= ',1pE16.6,' WEIGHT=',1pE16.6,' (electrostatic-local)'/
& 'ETURN3=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 3rd order)'/
& 'ETURN4=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 4th order)'/
+#ifdef FOURBODY
& 'ETURN6=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 6th order)'/
+#endif
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pE16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (virtual-bond dihedral angle restraints)'/
& 'ETHETC=',1pE16.6,' (virtual-bond angle restraints)'/
#else
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),
& estr,wbond,ebe,wang,escloc,wscloc,etors,wtor*fact(1),
- & etors_d,wtor_d*fact(2),ehpb,wstrain,ecorr,wcorr*fact(3),
+ & etors_d,wtor_d*fact(2),ehpb,
+#ifdef FOURBODY
+ & wstrain,ecorr,wcorr*fact(3),
& ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5),
+#endif
& eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2),
- & eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5),
+ & eello_turn4,wturn4*fact(3),
+#ifdef FOURBODY
+ & eello_turn6,wturn6*fact(5),
+#endif
& esccor,wsccor*fact(1),edihcnstr,
& ethetacnstr,ebr*nss,Uconst,wumb,eliptran,wliptran,Eafmforc,
& etube,wtube,esaxs,wsaxs,ehomology_constr,
& 'ETORSD=',1pE16.6,' WEIGHT=',1pE16.6,' (double torsional)'/
& 'EHBP= ',1pE16.6,' WEIGHT=',1pE16.6,
& ' (SS bridges & dist. restr.)'/
+#ifdef FOURBODY
& 'ECORR4=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
& 'ECORR5=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
& 'ECORR6=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
+#endif
& 'EELLO= ',1pE16.6,' WEIGHT=',1pE16.6,' (electrostatic-local)'/
& 'ETURN3=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 3rd order)'/
& 'ETURN4=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 4th order)'/
+#ifdef FOURBODY
& 'ETURN6=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 6th order)'/
+#endif
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pE16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (virtual-bond dihedral angle restraints)'/
& 'ETHETC=',1pE16.6,' (virtual-bond angle restraints)'/
include 'COMMON.SBRIDGE'
include 'COMMON.NAMES'
include 'COMMON.IOUNITS'
+#ifdef FOURBODY
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+#endif
dimension gg(3)
integer icant
external icant
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
C Change 12/1/95
num_conti=0
C
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
zj=c(3,nres+j)-zi
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
C Change 12/1/95 to calculate four-body interactions
rij=xj*xj+yj*yj+zj*zj
rrij=1.0D0/rij
+ sqrij=dsqrt(rij)
+ sss1=sscale(sqrij)
+ if (sss1.eq.0.0d0) cycle
+ sssgrad1=sscagrad(sqrij)
c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
cd & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm,
cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
if (bb.gt.0.0d0) then
- evdw=evdw+evdwij
+ evdw=evdw+sss1*evdwij
else
- evdw_t=evdw_t+evdwij
+ evdw_t=evdw_t+sss1*evdwij
endif
if (calc_grad) then
C
C Calculate the components of the gradient in DC and X
C
- fac=-rrij*(e1+evdwij)
+ fac=-rrij*(e1+evdwij)*sss1
+ & +evdwij*sssgrad1/sqrij/expon
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
enddo
enddo
endif
+#ifdef FOURBODY
C
C 12/1/95, revised on 5/20/97
C
cd & i,j,(gacont(kk,num_conti,i),kk=1,3)
endif
endif
+#endif
enddo ! j
enddo ! iint
+#ifdef FOURBODY
C Change 12/1/95
num_cont(i)=num_conti
+#endif
enddo ! i
if (calc_grad) then
do i=1,nct
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
C
C Calculate SC interaction energy.
C
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
zj=c(3,nres+j)-zi
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
fac_augm=rrij**expon
e_augm=augm(itypi,itypj)*fac_augm
r_inv_ij=dsqrt(rrij)
rij=1.0D0/r_inv_ij
+ sss1=sscale(rij)
+ if (sss1.eq.0.0d0) cycle
+ sssgrad1=sscagrad(rij)
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
e1=fac*fac*aa
cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
if (bb.gt.0.0d0) then
- evdw=evdw+evdwij
+ evdw=evdw+evdwij*sss1
else
- evdw_t=evdw_t+evdwij
+ evdw_t=evdw_t+evdwij*sss1
endif
if (calc_grad) then
C
C Calculate the components of the gradient in DC and X
C
- fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
+ fac=(-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2))*sss1
+ & +evdwij*sssgrad1*r_inv_ij/expon
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
c alf1=0.0D0
c alf2=0.0D0
c alf12=0.0D0
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
include "DIMENSIONS.COMPAR"
+ include 'COMMON.CONTROL'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.LOCAL'
yi=c(2,nres+i)
zi=c(3,nres+i)
C returning the ith atom to box
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
- if ((zi.gt.bordlipbot)
- &.and.(zi.lt.bordliptop)) then
-C the energy transfer exist
- if (zi.lt.buflipbot) then
-C what fraction I am in
- fracinbuf=1.0d0-
- & ((zi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
- sslipi=sscalelip(fracinbuf)
- ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
- elseif (zi.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
- sslipi=sscalelip(fracinbuf)
- ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
- else
- sslipi=1.0d0
- ssgradlipi=0.0
- endif
- else
- sslipi=0.0d0
- ssgradlipi=0.0
- endif
-
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
yj=c(2,nres+j)
zj=c(3,nres+j)
C returning jth atom to box
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- 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
- aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
- & +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 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
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
-C finding the closest
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- subchap=1
- endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
-
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +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
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
c write (iout,*) i,j,xj,yj,zj
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
- sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
- sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
+ sss=sscale(1.0d0/rij)
+ sssgrad=sscagrad(1.0d0/rij)
if (sss.le.0.0) cycle
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
#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
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
c alf1=0.0D0
c alf2=0.0D0
c alf12=0.0D0
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +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
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
+ sss=sscale(1.0d0/rij)
+ if (sss.eq.0.0d0) cycle
+ sssgrad=sscagrad(1.0d0/rij)
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
call sc_angular
e_augm=augm(itypi,itypj)*fac_augm
evdwij=evdwij*eps2rt*eps3rt
if (bb.gt.0.0d0) then
- evdw=evdw+evdwij+e_augm
+ evdw=evdw+(evdwij+e_augm)*sss
else
- evdw_t=evdw_t+evdwij+e_augm
+ evdw_t=evdw_t+(evdwij+e_augm)*sss
endif
ij=icant(itypi,itypj)
aux=eps1*eps2rt**2*eps3rt**2
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
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
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
+ include 'COMMON.CORRMAT'
double precision auxvec(2),auxmat(2,2)
C
C Compute the virtual-bond-torsional-angle dependent quantities needed
C
c write(iout,*) 'SET_MATRICES nphi=',nphi,nres
do i=3,nres+1
- if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ ii=ireschain(i-2)
+ if (ii.eq.0) cycle
+ innt=chain_border(1,ii)
+ inct=chain_border(2,ii)
+ if (i.gt. innt+2 .and. i.lt.inct+2) then
iti = itype2loc(itype(i-2))
else
iti=nloctyp
endif
c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
- if (i.gt. nnt+1 .and. i.lt.nct+1) then
+ if (i.gt. innt+1 .and. i.lt.inct+1) then
iti1 = itype2loc(itype(i-1))
else
iti1=nloctyp
write (iout,*) 'theta=', theta(i-1)
#endif
#else
+ if (i.gt. innt+2 .and. i.lt.inct+2) then
+c if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ iti = itype2loc(itype(i-2))
+ else
+ iti=nloctyp
+ endif
+c write (iout,*) "i",i-1," itype",itype(i-2)," iti",iti
+c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
+ if (i.gt. nnt+1 .and. i.lt.nct+1) then
+ iti1 = itype2loc(itype(i-1))
+ else
+ iti1=nloctyp
+ endif
c if (i.gt. nnt+2 .and. i.lt.nct+2) then
c iti = itype2loc(itype(i-2))
c else
c write(iout,*) "Macierz EUG",
c & eug(1,1,i-2),eug(1,2,i-2),eug(2,1,i-2),
c & eug(2,2,i-2)
+#ifdef FOURBODY
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0)
& then
call matmat2(CC(1,1,i-2),Ug(1,1,i-2),CUg(1,1,i-2))
call matvec2(Ctilde(1,1,i-1),obrot(1,i-2),Ctobr(1,i-2))
call matvec2(Dtilde(1,1,i-2),obrot2(1,i-2),Dtobr2(1,i-2))
endif
+#endif
else
do k=1,2
Ub2(k,i-2)=0.0d0
#endif
cd write (iout,*) 'mu1',mu1(:,i-2)
cd write (iout,*) 'mu2',mu2(:,i-2)
+#ifdef FOURBODY
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0)
& then
if (calc_grad) then
call matmat2(EUgder(1,1,i-2),DD(1,1,i-1),EUgDder(1,1,i-2))
endif
endif
+#endif
enddo
+#ifdef FOURBODY
C Matrices dependent on two consecutive virtual-bond dihedrals.
C The order of matrices is from left to right.
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0)
endif
enddo
endif
+#endif
return
end
C--------------------------------------------------------------------------
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
+#ifdef FOURBODY
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+#endif
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
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/
eello_turn3=0.0d0
eello_turn4=0.0d0
ind=0
+#ifdef FOURBODY
do i=1,nres
num_cont_hb(i)=0
enddo
+#endif
cd print '(a)','Enter EELEC'
c write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
c call flush(iout)
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
- 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
+ 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)
+#ifdef FOURBODY
num_cont_hb(i)=num_conti
+#endif
enddo
do i=iturn4_start,iturn4_end
if (i.lt.1) cycle
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
-C Return atom into box, boxxsize is size of box in x dimension
-c 194 continue
-c if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
-c if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
-C Condition for being inside the proper box
-c if ((xmedi.gt.((0.5d0)*boxxsize)).or.
-c & (xmedi.lt.((-0.5d0)*boxxsize))) then
-c go to 194
-c endif
-c 195 continue
-c if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize
-c if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize
-C Condition for being inside the proper box
-c if ((ymedi.gt.((0.5d0)*boxysize)).or.
-c & (ymedi.lt.((-0.5d0)*boxysize))) then
-c go to 195
-c endif
-c 196 continue
-c if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
-c if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
-C Condition for being inside the proper box
-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
-
+ call to_box(xmedi,ymedi,zmedi)
+ call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
+#ifdef FOURBODY
num_conti=num_cont_hb(i)
+#endif
c write(iout,*) "JESTEM W PETLI"
call eelecij(i,i+3,ees,evdw1,eel_loc)
if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
& call eturn4(i,eello_turn4)
+#ifdef FOURBODY
num_cont_hb(i)=num_conti
+#endif
enddo ! i
C Loop over all neighbouring boxes
C do xshift=-1,1
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
- 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
-C xmedi=xmedi+xshift*boxxsize
-C ymedi=ymedi+yshift*boxysize
-C zmedi=zmedi+zshift*boxzsize
-
-C Return tom into box, boxxsize is size of box in x dimension
-c 164 continue
-c if (xmedi.gt.((xshift+0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
-c if (xmedi.lt.((xshift-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
-C Condition for being inside the proper box
-c if ((xmedi.gt.((xshift+0.5d0)*boxxsize)).or.
-c & (xmedi.lt.((xshift-0.5d0)*boxxsize))) then
-c go to 164
-c endif
-c 165 continue
-c if (ymedi.gt.((yshift+0.5d0)*boxysize)) ymedi=ymedi-boxysize
-c if (ymedi.lt.((yshift-0.5d0)*boxysize)) ymedi=ymedi+boxysize
-C Condition for being inside the proper box
-c if ((ymedi.gt.((yshift+0.5d0)*boxysize)).or.
-c & (ymedi.lt.((yshift-0.5d0)*boxysize))) then
-c go to 165
-c endif
-c 166 continue
-c if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
-c if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
-cC Condition for being inside the proper box
-c if ((zmedi.gt.((zshift+0.5d0)*boxzsize)).or.
-c & (zmedi.lt.((zshift-0.5d0)*boxzsize))) then
-c go to 166
-c endif
-
-c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
+ call to_box(xmedi,ymedi,zmedi)
+ call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
+#ifdef FOURBODY
num_conti=num_cont_hb(i)
+#endif
C I TU KURWA
do j=ielstart(i),ielend(i)
C do j=16,17
&) cycle
call eelecij(i,j,ees,evdw1,eel_loc)
enddo ! j
+#ifdef FOURBODY
num_cont_hb(i)=num_conti
+#endif
enddo ! i
C enddo ! zshift
C enddo ! yshift
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
+#ifdef FOURBODY
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+#endif
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
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/
xj=c(1,j)+0.5D0*dxj
yj=c(2,j)+0.5D0*dyj
zj=c(3,j)+0.5D0*dzj
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- 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"
- dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- isubchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- isubchap=1
- endif
- enddo
- enddo
- enddo
- if (isubchap.eq.1) then
- xj=xj_temp-xmedi
- yj=yj_temp-ymedi
- zj=zj_temp-zmedi
- else
- xj=xj_safe-xmedi
- yj=yj_safe-ymedi
- zj=zj_safe-zmedi
- endif
-C if ((i+3).lt.j) then !this condition keeps for turn3 and turn4 not subject to PBC
-c 174 continue
-c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
-c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
-C Condition for being inside the proper box
-c if ((xj.gt.((0.5d0)*boxxsize)).or.
-c & (xj.lt.((-0.5d0)*boxxsize))) then
-c go to 174
-c endif
-c 175 continue
-c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
-c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
-C Condition for being inside the proper box
-c if ((yj.gt.((0.5d0)*boxysize)).or.
-c & (yj.lt.((-0.5d0)*boxysize))) then
-c go to 175
-c endif
-c 176 continue
-c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
-c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
-C Condition for being inside the proper box
-c if ((zj.gt.((0.5d0)*boxzsize)).or.
-c & (zj.lt.((-0.5d0)*boxzsize))) then
-c go to 176
-c endif
-C endif !endPBC condintion
-C xj=xj-xmedi
-C yj=yj-ymedi
-C zj=zj-zmedi
+ 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)
rij=xj*xj+yj*yj+zj*zj
-
- sss=sscale(sqrt(rij))
- sssgrad=sscagrad(sqrt(rij))
+ sss=sscale(sqrt(rij))
+ if (sss.eq.0.0d0) return
+ sssgrad=sscagrad(sqrt(rij))
c write (iout,*) "ij",i,j," rij",sqrt(rij)," r_cut",r_cut,
c & " rlamb",rlamb," sss",sss
c if (sss.gt.0.0d0) then
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*sss*faclipij2
else
fac_shield(i)=1.0
fac_shield(j)=1.0
eesij=(el1+el2)
- ees=ees+eesij
+ ees=ees+eesij*sss*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
- ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
- ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
- ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
+ facvdw=facvdw+sssgrad*rmij*evdwij*faclipij2
+ ggg(1)=facvdw*xj
+ ggg(2)=facvdw*yj
+ ggg(3)=facvdw*zj
else
ggg(1)=0.0
ggg(2)=0.0
endif ! calc_grad
#else
C MARYSIA
- facvdw=(ev1+evdwij)*sss
+ facvdw=(ev1+evdwij)*faclipij2
facel=(el1+eesij)
fac1=fac
- fac=-3*rrmij*(facvdw+facvdw+facel)
+ fac=-3*rrmij*(facvdw+facvdw+facel)*sss
+ & +(evdwij+eesij)*sssgrad*rrmij
erij(1)=xj*rmij
erij(2)=yj*rmij
erij(3)=zj*rmij
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)
+ & *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)
& +a23*gmuij1(2)
& +a32*gmuij1(3)
& +a33*gmuij1(4))
- & *fac_shield(i)*fac_shield(j)
+ & *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)
+ & *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)
+ & *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
+ ggg(2)=aux*yj
+ ggg(3)=aux*zj
do l=1,3
- ggg(l)=(agg(l,1)*muij(1)+
+ 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)
+ & *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
C Change 12/26/95 to calculate four-body contributions to H-bonding energy
c if (j.gt.i+1 .and. num_conti.le.maxconts) then
+#ifdef FOURBODY
if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
& .and. num_conti.le.maxconts) then
c write (iout,*) i,j," entered corr"
gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
enddo
gggp(1)=gggp(1)+ees0pijp*xj
+ & +ees0p(num_conti,i)/sss*rmij*xj*sssgrad
gggp(2)=gggp(2)+ees0pijp*yj
+ & +ees0p(num_conti,i)/sss*rmij*yj*sssgrad
gggp(3)=gggp(3)+ees0pijp*zj
+ & +ees0p(num_conti,i)/sss*rmij*zj*sssgrad
gggm(1)=gggm(1)+ees0mijp*xj
+ & +ees0m(num_conti,i)/sss*rmij*xj*sssgrad
gggm(2)=gggm(2)+ees0mijp*yj
+ & +ees0m(num_conti,i)/sss*rmij*yj*sssgrad
gggm(3)=gggm(3)+ees0mijp*zj
+ & +ees0m(num_conti,i)/sss*rmij*zj*sssgrad
C Derivatives due to the contact function
gacont_hbr(1,num_conti,i)=fprimcont*xj
gacont_hbr(2,num_conti,i)=fprimcont*yj
gacontp_hb1(k,num_conti,i)=!ghalfp
& +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
& + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gacontp_hb2(k,num_conti,i)=!ghalfp
& +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
& + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gacontp_hb3(k,num_conti,i)=gggp(k)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
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)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gacontm_hb2(k,num_conti,i)=!ghalfm
& +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
& + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gacontm_hb3(k,num_conti,i)=gggm(k)
& *fac_shield(i)*fac_shield(j)
-
+*sss
enddo
C Diagnostics. Comment out or remove after debugging!
cdiag do k=1,3
endif ! num_conti.le.maxconts
endif ! fcont.gt.0
endif ! j.gt.i+1
+#endif
if (calc_grad) then
if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
do k=1,4
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
- include 'COMMON.CONTACTS'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
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
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
- include 'COMMON.CONTACTS'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
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+3
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
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
yi=0.5D0*(c(2,i)+c(2,i+1))
zi=0.5D0*(c(3,i)+c(3,i+1))
C Returning the ith atom to box
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
+ call to_box(xi,yi,zi)
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
yj=c(2,j)
zj=c(3,j)
C returning the jth atom to box
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- 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
- yj_safe=yj
- zj_safe=zj
- subchap=0
-C Finding the closest jth atom
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- subchap=1
- endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
C sss is scaling function for smoothing the cutoff gradient otherwise
C the gradient would not be continuouse
estr1=0.0d0
c write (iout,*) "distchainmax",distchainmax
do i=nnt+1,nct
+#ifdef FIVEDIAG
+ if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) cycle
+ diff = vbld(i)-vbldp0
+#else
if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
C estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
C do j=1,3
diff = vbld(i)-vbldp0
c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff
endif
+#endif
+ if (energy_dec) write (iout,'(a7,i5,4f7.3)')
+ & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
estr=estr+diff*diff
do j=1,3
gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
nbi=nbondterm(iti)
if (nbi.eq.1) then
diff=vbld(i+nres)-vbldsc0(1,iti)
-C write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
-C & AKSC(1,iti),AKSC(1,iti)*diff*diff
+ if (energy_dec) write (iout,*)
+ & i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
+ & AKSC(1,iti),AKSC(1,iti)*diff*diff
estr=estr+0.5d0*AKSC(1,iti)*diff*diff
do j=1,3
gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
- if (i.eq.3) then
- phii=0.0d0
- ityp1=nthetyp+1
- do k=1,nsingle
- cosph1(k)=0.0d0
- sinph1(k)=0.0d0
- enddo
- else
+cu if (i.eq.3) then
+cu phii=0.0d0
+cu ityp1=nthetyp+1
+cu do k=1,nsingle
+cu cosph1(k)=0.0d0
+cu sinph1(k)=0.0d0
+cu enddo
+cu else
if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
sinph1(k)=0.0d0
enddo
endif
- endif
if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
enddo
return
end
+#ifdef FOURBODY
c------------------------------------------------------------------------------
subroutine multibody(ecorr)
C This subroutine calculates multi-body contributions to energy following
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
double precision gx(3),gx1(3)
logical lprn
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
double precision gx(3),gx1(3)
logical lprn
lprn=.false.
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
double precision gx(3),gx1(3)
logical lprn,ldone
include 'COMMON.LOCAL'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.CHAIN'
include 'COMMON.CONTROL'
include 'COMMON.SHIELD'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.SHIELD'
include 'COMMON.CONTROL'
double precision gx(3),gx1(3)
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VAR'
include 'COMMON.GEO'
cd write (2,*) 'eel_turn6',ekont*eel_turn6
return
end
-
+#endif
crc-------------------------------------------------
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine Eliptransfer(eliptran)
enddo
c min_odl=minval(distancek)
- do kk=1,constr_homology
- if(l_homo(kk,ii)) then
- min_odl=distancek(kk)
- exit
- endif
- enddo
- do kk=1,constr_homology
- if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl)
+ if (nexl.gt.0) then
+ min_odl=0.0d0
+ else
+ do kk=1,constr_homology
+ if(l_homo(kk,ii)) then
+ min_odl=distancek(kk)
+ exit
+ endif
+ enddo
+ do kk=1,constr_homology
+ if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl)
& min_odl=distancek(kk)
- enddo
+ enddo
+ endif
c write (iout,* )"min_odl",min_odl
#ifdef DEBUG
write (iout,*) "ij dij",i,j,dij
projects / unres.git / blobdiff