include 'COMMON.INTERACT'
include 'COMMON.SBRIDGE'
include 'COMMON.CHAIN'
+ include 'COMMON.SHIELD'
+ include 'COMMON.CONTROL'
double precision fact(6)
cd write(iout, '(a,i2)')'Calling etotal ipot=',ipot
cd print *,'nnt=',nnt,' nct=',nct
C
C Calculate electrostatic (H-bonding) energy of the main chain.
C
- 106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
+ 106 continue
+C write(iout,*) "shield_mode",shield_mode,ethetacnstr
+ if (shield_mode.eq.1) then
+ call set_shield_fac
+ else if (shield_mode.eq.2) then
+ call set_shield_fac2
+ endif
+ call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
C
C Calculate excluded-volume interaction energy between peptide groups
C and side chains.
C
C Calculate the virtual-bond-angle energy.
C
- call ebend(ebe)
+ call ebend(ebe,ethetacnstr)
cd print *,'Bend energy finished.'
C
C Calculate the SC local energy.
C 21/5/07 Calculate local sicdechain correlation energy
C
call eback_sc_corr(esccor)
+
+ if (wliptran.gt.0) then
+ call Eliptransfer(eliptran)
+ endif
+
C
C 12/1/95 Multi-body terms
C
call multibody_eello(ecorr,ecorr5,ecorr6,eturn6,n_corr,n_corr1)
c write (*,*) 'n_corr=',n_corr,' n_corr1=',n_corr1
c print *,ecorr,ecorr5,ecorr6,eturn6
+ else
+ ecorr=0.0d0
+ ecorr5=0.0d0
+ ecorr6=0.0d0
+ eturn6=0.0d0
endif
if (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) then
call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
endif
-c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
+
+c write(iout,*) "TEST_ENE",constr_homology
+ if (constr_homology.ge.1) then
+ call e_modeller(ehomology_constr)
+ else
+ ehomology_constr=0.0d0
+ endif
+c write(iout,*) "TEST_ENE",ehomology_constr
+
+
+ write (iout,*) "ft(6)",fact(6),wliptran,eliptran
#ifdef SPLITELE
+ if (shield_mode.gt.0) then
+ etot=fact(1)*wsc*(evdw+fact(6)*evdw_t)+fact(1)*wscp*evdw2
+ & +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
+ & +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
+ & +wliptran*eliptran
+ else
etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2+welec*fact(1)*ees
& +wvdwpp*evdw1
& +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+ & +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
+ & +wbond*estr+wsccor*fact(1)*esccor+ehomology_constr
+ & +wliptran*eliptran
+ endif
#else
+ if (shield_mode.gt.0) then
+ etot=fact(1)wsc*(evdw+fact(6)*evdw_t)+fact(1)*wscp*evdw2
+ & +welec*fact(1)*(ees+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
+ & +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+ehomology_constr
+ & +wliptran*eliptran
+ else
etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2
& +welec*fact(1)*(ees+evdw1)
& +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+ & +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
+ & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
+ & +wliptran*eliptran
+ endif
#endif
+
energia(0)=etot
energia(1)=evdw
-c call enerprint(energia(0),frac)
#ifdef SCP14
energia(2)=evdw2-evdw2_14
energia(17)=evdw2_14
energia(18)=estr
energia(19)=esccor
energia(20)=edihcnstr
+ energia(24)=ehomology_constr
energia(21)=evdw_t
+c energia(24)=ethetacnstr
+ energia(22)=eliptran
c detecting NaNQ
#ifdef ISNAN
#ifdef AIX
#ifdef SPLITELE
do i=1,nct
do j=1,3
+ if (shield_mode.eq.0) then
gradc(j,i,icg)=wsc*gvdwc(j,i)+wscp*gvdwc_scp(j,i)+
& welec*fact(1)*gelc(j,i)+wvdwpp*gvdwpp(j,i)+
& wbond*gradb(j,i)+
& wcorr6*fact(5)*gradcorr6(j,i)+
& wturn6*fact(5)*gcorr6_turn(j,i)+
& wsccor*fact(2)*gsccorc(j,i)
+ & +wliptran*gliptranc(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)
+ else
+ gradc(j,i,icg)=fact(1)*wsc*gvdwc(j,i)
+ & +fact(1)*wscp*gvdwc_scp(j,i)+
+ & welec*fact(1)*gelc(j,i)+fact(1)*wvdwpp*gvdwpp(j,i)+
+ & wbond*gradb(j,i)+
+ & wstrain*ghpbc(j,i)+
+ & wcorr*fact(3)*gradcorr(j,i)+
+ & wel_loc*fact(2)*gel_loc(j,i)+
+ & wturn3*fact(2)*gcorr3_turn(j,i)+
+ & wturn4*fact(3)*gcorr4_turn(j,i)+
+ & wcorr5*fact(4)*gradcorr5(j,i)+
+ & wcorr6*fact(5)*gradcorr6(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
#else
- do i=1,nct
+ do i=1,nct
do j=1,3
+ if (shield_mode.eq.0) then
gradc(j,i,icg)=wsc*gvdwc(j,i)+wscp*gvdwc_scp(j,i)+
& welec*fact(1)*gelc(j,i)+wstrain*ghpbc(j,i)+
& wbond*gradb(j,i)+
& wcorr6*fact(5)*gradcorr6(j,i)+
& wturn6*fact(5)*gcorr6_turn(j,i)+
& wsccor*fact(2)*gsccorc(j,i)
+ & +wliptran*gliptranc(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)
- enddo
+ & +wliptran*gliptranx(j,i)
+ else
+ gradc(j,i,icg)=fact(1)*wsc*gvdwc(j,i)+
+ & fact(1)*wscp*gvdwc_scp(j,i)+
+ & welec*fact(1)*gelc(j,i)+wstrain*ghpbc(j,i)+
+ & wbond*gradb(j,i)+
+ & wcorr*fact(3)*gradcorr(j,i)+
+ & wel_loc*fact(2)*gel_loc(j,i)+
+ & wturn3*fact(2)*gcorr3_turn(j,i)+
+ & wturn4*fact(3)*gcorr4_turn(j,i)+
+ & wcorr5*fact(4)*gradcorr5(j,i)+
+ & wcorr6*fact(5)*gradcorr6(j,i)+
+ & wturn6*fact(5)*gcorr6_turn(j,i)+
+ & wsccor*fact(2)*gsccorc(j,i)
+ & +wliptran*gliptranc(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)
+ endif
+ enddo
#endif
enddo
& +wturn3*fact(2)*gel_loc_turn3(i)
& +wturn6*fact(5)*gel_loc_turn6(i)
& +wel_loc*fact(2)*gel_loc_loc(i)
- & +wsccor*fact(1)*gsccor_loc(i)
+c & +wsccor*fact(1)*gsccor_loc(i)
+c ROZNICA Z WHAMem
enddo
endif
+ if (dyn_ss) call dyn_set_nss
return
end
C------------------------------------------------------------------------
esccor=energia(19)
edihcnstr=energia(20)
estr=energia(18)
+ ehomology_constr=energia(24)
+c ethetacnstr=energia(24)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1,
& wvdwpp,
& ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5),
& eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2),
& eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5),
- & esccor,wsccor*fact(1),edihcnstr,ebr*nss,etot
+ & esccor,wsccor*fact(1),edihcnstr,ehomology_constr,ebr*nss,
+ & etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
+ & 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#else
& ecorr6,wcorr6*fact(5),eel_loc,wel_loc*fact(2),
& eello_turn3,wturn3*fact(2),eello_turn4,wturn4*fact(3),
& eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor,
- & edihcnstr,ebr*nss,etot
+ & edihcnstr,ehomology_constr,ebr*nss,
+ & etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
+ & 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#endif
integer icant
external icant
cd print *,'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
+c ROZNICA DODANE Z WHAM
+c do i=1,210
+c do j=1,2
+c eneps_temp(j,i)=0.0d0
+c enddo
+c enddo
+cROZNICA
+
evdw=0.0D0
evdw_t=0.0d0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
- if (itypi.eq.21) cycle
- itypi1=itype(i+1)
+ itypi=iabs(itype(i))
+ if (itypi.eq.ntyp1) cycle
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
cd & 'iend=',iend(i,iint)
do j=istart(i,iint),iend(i,iint)
- itypj=itype(j)
- if (itypj.eq.21) cycle
+ itypj=iabs(itype(j))
+ if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
zj=c(3,nres+j)-zi
c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e1+e2
ij=icant(itypi,itypj)
+c ROZNICA z WHAM
+c eneps_temp(1,ij)=eneps_temp(1,ij)+e1/dabs(eps0ij)
+c eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps0ij
+c
+
cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
cd write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)')
cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
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(itypi,itypj).gt.0.0d0) then
+ if (bb.gt.0.0d0) then
evdw=evdw+evdwij
else
evdw_t=evdw_t+evdwij
evdw=0.0D0
evdw_t=0.0d0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
- if (itypi.eq.21) cycle
- itypi1=itype(i+1)
+ itypi=iabs(itype(i))
+ if (itypi.eq.ntyp1) cycle
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
C
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
- itypj=itype(j)
- if (itypj.eq.21) cycle
+ itypj=iabs(itype(j))
+ if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
zj=c(3,nres+j)-zi
rij=1.0D0/r_inv_ij
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e_augm+e1+e2
ij=icant(itypi,itypj)
cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
- if (bb(itypi,itypj).gt.0.0d0) then
+ if (bb.gt.0.0d0) then
evdw=evdw+evdwij
else
evdw_t=evdw_t+evdwij
c endif
ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
- if (itypi.eq.21) cycle
- itypi1=itype(i+1)
+ itypi=iabs(itype(i))
+ if (itypi.eq.ntyp1) cycle
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
ind=ind+1
- itypj=itype(j)
- if (itypj.eq.21) cycle
+ itypj=iabs(itype(j))
+ if (itypj.eq.ntyp1) cycle
dscj_inv=vbld_inv(j+nres)
chi1=chi(itypi,itypj)
chi2=chi(itypj,itypi)
C Calculate whole angle-dependent part of epsilon and contributions
C to its derivatives
fac=(rrij*sigsq)**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
ij=icant(itypi,itypj)
aux=eps1*eps2rt**2*eps3rt**2
- if (bb(itypi,itypj).gt.0.0d0) then
+ if (bb.gt.0.0d0) then
evdw=evdw+evdwij
else
evdw_t=evdw_t+evdwij
endif
if (calc_grad) then
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
cd & restyp(itypi),i,restyp(itypj),j,
cd & epsi,sigm,chi1,chi2,chip1,chip2,
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
+ include 'COMMON.SBRIDGE'
logical lprn
common /srutu/icall
integer icant
external icant
+ integer xshift,yshift,zshift
+ logical energy_dec /.false./
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
evdw_t=0.0d0
c if (icall.gt.0) lprn=.true.
ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
- if (itypi.eq.21) cycle
- itypi1=itype(i+1)
+ itypi=iabs(itype(i))
+ if (itypi.eq.ntyp1) cycle
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ 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
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
C
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
+ IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+
+c write(iout,*) "PRZED ZWYKLE", evdwij
+ call dyn_ssbond_ene(i,j,evdwij)
+c write(iout,*) "PO ZWYKLE", evdwij
+
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+ & 'evdw',i,j,evdwij,' ss'
+C triple bond artifac removal
+ do k=j+1,iend(i,iint)
+C search over all next residues
+ if (dyn_ss_mask(k)) then
+C check if they are cysteins
+C write(iout,*) 'k=',k
+
+c write(iout,*) "PRZED TRI", evdwij
+ evdwij_przed_tri=evdwij
+ call triple_ssbond_ene(i,j,k,evdwij)
+c if(evdwij_przed_tri.ne.evdwij) then
+c write (iout,*) "TRI:", evdwij, evdwij_przed_tri
+c endif
+
+c write(iout,*) "PO TRI", evdwij
+C call the energy function that removes the artifical triple disulfide
+C bond the soubroutine is located in ssMD.F
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+ & 'evdw',i,j,evdwij,'tss'
+ endif!dyn_ss_mask(k)
+ enddo! k
+ ELSE
ind=ind+1
- itypj=itype(j)
- if (itypj.eq.21) cycle
+ itypj=iabs(itype(j))
+ if (itypj.eq.ntyp1) cycle
dscj_inv=vbld_inv(j+nres)
sig0ij=sigma(itypi,itypj)
chi1=chi(itypi,itypj)
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)
+ 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 write(iout,*) "czy jest 0", bb-bb_lip(itypi,itypj),
+C & bb-bb_aq(itypi,itypj)
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=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-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
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))
+ if (sss.le.0.0d0) cycle
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
call sc_angular
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
- if (bb(itypi,itypj).gt.0) then
- evdw=evdw+evdwij
+ if (bb.gt.0) then
+ evdw=evdw+evdwij*sss
else
- evdw_t=evdw_t+evdwij
+ evdw_t=evdw_t+evdwij*sss
endif
ij=icant(itypi,itypj)
aux=eps1*eps2rt**2*eps3rt**2
c & " ij",ij," eneps",aux*e1/dabs(eps(itypi,itypj)),
c & aux*e2/eps(itypi,itypj)
c if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
-c write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-c & restyp(itypi),i,restyp(itypj),j,
-c & epsi,sigm,chi1,chi2,chip1,chip2,
-c & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
-c & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
-c & evdwij
-c write (iout,*) "pratial sum", evdw,evdw_t
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
+C#define DEBUG
+#ifdef DEBUG
+C write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+C & restyp(itypi),i,restyp(itypj),j,
+C & epsi,sigm,chi1,chi2,chip1,chip2,
+C & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
+C & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
+C & evdwij
+ write (iout,*) "pratial sum", evdw,evdw_t,e1,e2,fac,aa
+#endif
+C#undef DEBUG
c endif
if (calc_grad) then
C Calculate gradient components.
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac
+ fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
+ 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
C Calculate the radial part of the gradient
gg(1)=xj*fac
gg(2)=yj*fac
C Calculate angular part of the gradient.
call sc_grad
endif
+ ENDIF ! dyn_ss
enddo ! j
enddo ! iint
enddo ! i
c if (icall.gt.0) lprn=.true.
ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
- if (itypi.eq.21) cycle
- itypi1=itype(i+1)
+ itypi=iabs(itype(i))
+ if (itypi.eq.ntyp1) cycle
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
dsci_inv=vbld_inv(i+nres)
+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
C
C Calculate SC interaction energy.
C
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
ind=ind+1
- itypj=itype(j)
- if (itypj.eq.21) cycle
+ itypj=iabs(itype(j))
+ if (itypj.eq.ntyp1) cycle
dscj_inv=vbld_inv(j+nres)
sig0ij=sigma(itypi,itypj)
r0ij=r0(itypi,itypj)
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)
+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 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
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
fac_augm=rrij**expon
e_augm=augm(itypi,itypj)*fac_augm
evdwij=evdwij*eps2rt*eps3rt
- if (bb(itypi,itypj).gt.0.0d0) then
+ if (bb.gt.0.0d0) then
evdw=evdw+evdwij+e_augm
else
evdw_t=evdw_t+evdwij+e_augm
gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
enddo
do k=1,3
- gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
& +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
& +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
- gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipi(k)
& +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
& +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
enddo
C
do k=i,j-1
do l=1,3
- gvdwc(l,k)=gvdwc(l,k)+gg(l)
+ gvdwc(l,k)=gvdwc(l,k)+gg(l)+gg_lipi(l)
enddo
enddo
+ do l=1,3
+ gvdwc(l,j)=gvdwc(l,j)+gg_lipj(l)
+ enddo
return
end
c------------------------------------------------------------------------------
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
+ 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),
gcorr_loc(i)=0.0d0
enddo
do i=iatel_s,iatel_e
- if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle
+ if (i.le.1) cycle
+ if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
+ & .or. ((i+2).gt.nres)
+ & .or. ((i-1).le.0)
+ & .or. itype(i+2).eq.ntyp1
+ & .or. itype(i-1).eq.ntyp1
+ &) cycle
+C endif
if (itel(i).eq.0) goto 1215
dxi=dc(1,i)
dyi=dc(2,i)
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
num_conti=0
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
do j=ielstart(i),ielend(i)
- if (itype(j).eq.21 .or. itype(j+1).eq.21) cycle
+ if (j.le.1) cycle
+ if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
+ & .or.((j+2).gt.nres)
+ & .or.((j-1).le.0)
+ & .or.itype(j+2).eq.ntyp1
+ & .or.itype(j-1).eq.ntyp1
+ &) cycle
+C endif
if (itel(j).eq.0) goto 1216
ind=ind+1
iteli=itel(i)
dx_normj=dc_norm(1,j)
dy_normj=dc_norm(2,j)
dz_normj=dc_norm(3,j)
- xj=c(1,j)+0.5D0*dxj-xmedi
- yj=c(2,j)+0.5D0*dyj-ymedi
- zj=c(3,j)+0.5D0*dzj-zmedi
+ 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
+ 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
+
rij=xj*xj+yj*yj+zj*zj
+ sss=sscale(sqrt(rij))
+ sssgrad=sscagrad(sqrt(rij))
rrmij=1.0D0/rij
rij=dsqrt(rij)
rmij=1.0D0/rij
c write (iout,*) "i",i,iteli," j",j,itelj," eesij",eesij
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 fac_shield(i)=0.4
+C fac_shield(j)=0.6
+C#define DEBUG
+#ifdef DEBUG
+ write(iout,*) "ees_compon",i,j,el1,el2,
+ & fac_shield(i),fac_shield(j)
+#endif
+C#undef DEBUG
+ 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
+ else
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+ eesij=(el1+el2)
ees=ees+eesij
- evdw1=evdw1+evdwij
+ endif
+C ees=ees+eesij
+ evdw1=evdw1+evdwij*sss
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,
C Calculate contributions to the Cartesian gradient.
C
#ifdef SPLITELE
- facvdw=-6*rrmij*(ev1+evdwij)
+ facvdw=-6*rrmij*(ev1+evdwij)*sss
facel=-3*rrmij*(el1+eesij)
fac1=fac
erij(1)=xj*rmij
ggg(1)=facel*xj
ggg(2)=facel*yj
ggg(3)=facel*zj
+
+ if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+ & (shield_mode.gt.0)) then
+C print *,i,j
+ do ilist=1,ishield_list(i)
+ iresshield=shield_list(ilist,i)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
+ & *2.0
+ gshieldx(k,iresshield)=gshieldx(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0
+ gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
+C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
+C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C if (iresshield.gt.i) then
+C do ishi=i+1,iresshield-1
+C gshieldc(k,ishi)=gshieldc(k,ishi)+rlocshield
+C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C
+C enddo
+C else
+C do ishi=iresshield,i
+C gshieldc(k,ishi)=gshieldc(k,ishi)-rlocshield
+C & -grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
+C
+C enddo
+C endif
+C enddo
+C enddo
+ enddo
+ enddo
+ do ilist=1,ishield_list(j)
+ iresshield=shield_list(ilist,j)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
+ & *2.0
+ gshieldx(k,iresshield)=gshieldx(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
+ gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
+ enddo
+ enddo
+
+ do k=1,3
+ gshieldc(k,i)=gshieldc(k,i)+
+ & grad_shield(k,i)*eesij/fac_shield(i)*2.0
+ gshieldc(k,j)=gshieldc(k,j)+
+ & grad_shield(k,j)*eesij/fac_shield(j)*2.0
+ gshieldc(k,i-1)=gshieldc(k,i-1)+
+ & grad_shield(k,i)*eesij/fac_shield(i)*2.0
+ gshieldc(k,j-1)=gshieldc(k,j-1)+
+ & grad_shield(k,j)*eesij/fac_shield(j)*2.0
+
+ enddo
+ endif
+
do k=1,3
ghalf=0.5D0*ggg(k)
gelc(k,i)=gelc(k,i)+ghalf
gelc(l,k)=gelc(l,k)+ggg(l)
enddo
enddo
- ggg(1)=facvdw*xj
- ggg(2)=facvdw*yj
- ggg(3)=facvdw*zj
+C ggg(1)=facvdw*xj
+C ggg(2)=facvdw*yj
+C ggg(3)=facvdw*zj
+ 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
+ else
+ ggg(1)=0.0
+ ggg(2)=0.0
+ ggg(3)=0.0
+ endif
do k=1,3
ghalf=0.5D0*ggg(k)
gvdwpp(k,i)=gvdwpp(k,i)+ghalf
enddo
enddo
#else
- facvdw=ev1+evdwij
+ facvdw=(ev1+evdwij)*sss
facel=el1+eesij
fac1=fac
fac=-3*rrmij*(facvdw+facvdw+facel)
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
enddo
do k=1,3
ghalf=0.5D0*ggg(k)
gelc(k,i)=gelc(k,i)+ghalf
& +(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
+
gelc(k,j)=gelc(k,j)+ghalf
& +(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
enddo
do k=i+1,j-1
do l=1,3
& +a33*muij(4)
cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
cd write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+C else
+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)
eel_loc=eel_loc+eel_loc_ij
C Partial derivatives in virtual-bond dihedral angles gamma
if (calc_grad) then
+ if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+ & (shield_mode.gt.0)) then
+C print *,i,j
+
+ do ilist=1,ishield_list(i)
+ iresshield=shield_list(ilist,i)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,i)*eel_loc_ij
+ & /fac_shield(i)
+C & *2.0
+ gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
+ gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
+ & +rlocshield
+ enddo
+ enddo
+ do ilist=1,ishield_list(j)
+ iresshield=shield_list(ilist,j)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,j)*eel_loc_ij
+ & /fac_shield(j)
+C & *2.0
+ gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
+ gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
+ & +rlocshield
+
+ enddo
+ enddo
+ do k=1,3
+ gshieldc_ll(k,i)=gshieldc_ll(k,i)+
+ & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
+ gshieldc_ll(k,j)=gshieldc_ll(k,j)+
+ & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
+ gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+
+ & grad_shield(k,i)*eel_loc_ij/fac_shield(i)
+ gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
+ & grad_shield(k,j)*eel_loc_ij/fac_shield(j)
+ enddo
+ endif
if (i.gt.1)
& 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)
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)
+
cd call checkint3(i,j,mu1,mu2,a22,a23,a32,a33,acipa,eel_loc_ij)
cd write(iout,*) 'agg ',agg
cd write(iout,*) 'aggi ',aggi
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)
+
enddo
do k=i+2,j2
do l=1,3
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)
+
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)
+
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)
+
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)
+
enddo
endif
ENDIF
fac3=dsqrt(-ael6i)*r3ij
ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0d0
+ fac_shield(j)=1.0d0
+ else
+ ees0plist(num_conti,i)=j
+C fac_shield(i)=0.4d0
+C fac_shield(j)=0.6d0
+ endif
c ees0mij=0.0D0
ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
+ & *fac_shield(i)*fac_shield(j)
+
ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
+ & *fac_shield(i)*fac_shield(j)
+
C Diagnostics. Comment out or remove after debugging!
c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
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)
+
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)
+
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)
+ & *fac_shield(i)*fac_shield(j)
+
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)
+
gacontm_hb3(k,num_conti,i)=gggm(k)
+ & *fac_shield(i)*fac_shield(j)
+
enddo
endif
C Diagnostics. Comment out or remove after debugging!
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
+ include 'COMMON.SHIELD'
+ include 'COMMON.CONTROL'
+
dimension ggg(3)
double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2),
& e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2),
& aggj(3,4),aggj1(3,4),a_temp(2,2)
common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2
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.((i+5).gt.nres)
+C & .or.((i-1).le.0)
+C end of changes suggested by Ana
+ & .or. itype(i+2).eq.ntyp1
+ & .or. itype(i+3).eq.ntyp1
+C & .or. itype(i+5).eq.ntyp1
+C & .or. itype(i).eq.ntyp1
+C & .or. itype(i-1).eq.ntyp1
+ & ) goto 179
+
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Third-order contributions
call matmat2(EUg(1,1,i+1),EUg(1,1,i+2),auxmat(1,1))
call transpose2(auxmat(1,1),auxmat1(1,1))
call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+C else
+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)
+ eello_t3=0.5d0*(pizda(1,1)+pizda(2,2))
+ & *fac_shield(i)*fac_shield(j)
+
cd write (2,*) 'i,',i,' j',j,'eello_turn3',
cd & 0.5d0*(pizda(1,1)+pizda(2,2)),
cd & ' eello_turn3_num',4*eello_turn3_num
if (calc_grad) then
+C Derivatives in shield mode
+ if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+ & (shield_mode.gt.0)) then
+C print *,i,j
+
+ do ilist=1,ishield_list(i)
+ iresshield=shield_list(ilist,i)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,i)*eello_t3/fac_shield(i)
+C & *2.0
+ gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,i)*eello_t3/fac_shield(i)
+ gshieldc_t3(k,iresshield-1)=gshieldc_t3(k,iresshield-1)
+ & +rlocshield
+ enddo
+ enddo
+ do ilist=1,ishield_list(j)
+ iresshield=shield_list(ilist,j)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,j)*eello_t3/fac_shield(j)
+C & *2.0
+ gshieldx_t3(k,iresshield)=gshieldx_t3(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,j)*eello_t3/fac_shield(j)
+ gshieldc_t3(k,iresshield-1)=gshieldc_t3(k,iresshield-1)
+ & +rlocshield
+
+ enddo
+ enddo
+
+ do k=1,3
+ gshieldc_t3(k,i)=gshieldc_t3(k,i)+
+ & grad_shield(k,i)*eello_t3/fac_shield(i)
+ gshieldc_t3(k,j)=gshieldc_t3(k,j)+
+ & grad_shield(k,j)*eello_t3/fac_shield(j)
+ gshieldc_t3(k,i-1)=gshieldc_t3(k,i-1)+
+ & grad_shield(k,i)*eello_t3/fac_shield(i)
+ gshieldc_t3(k,j-1)=gshieldc_t3(k,j-1)+
+ & grad_shield(k,j)*eello_t3/fac_shield(j)
+ enddo
+ endif
+
C Derivatives in gamma(i)
call matmat2(EUgder(1,1,i+1),EUg(1,1,i+2),auxmat2(1,1))
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)
+
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))
call matmat2(a_temp(1,1),pizda(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)
+
C Cartesian derivatives
do l=1,3
a_temp(1,1)=aggi(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)
+
a_temp(1,1)=aggi1(l,1)
a_temp(1,2)=aggi1(l,2)
a_temp(2,1)=aggi1(l,3)
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)
+
a_temp(1,1)=aggj(l,1)
a_temp(1,2)=aggj(l,2)
a_temp(2,1)=aggj(l,3)
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)
+
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)
+
enddo
endif
- else if (j.eq.i+3 .and. itype(i+2).ne.21) then
+ 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
+C & .or.((i+5).gt.nres)
+C & .or.((i-1).le.0)
+C end of changes suggested by Ana
+ & .or. itype(i+3).eq.ntyp1
+ & .or. itype(i+4).eq.ntyp1
+C & .or. itype(i+5).eq.ntyp1
+ & .or. itype(i).eq.ntyp1
+C & .or. itype(i-1).eq.ntyp1
+ & ) goto 178
+
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Fourth-order contributions
call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+C else
+C fac_shield(i)=0.4
+C fac_shield(j)=0.6
+ endif
eello_turn4=eello_turn4-(s1+s2+s3)
+ & *fac_shield(i)*fac_shield(j)
+ eello_t4=-(s1+s2+s3)
+ & *fac_shield(i)*fac_shield(j)
+
cd write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
cd & ' eello_turn4_num',8*eello_turn4_num
C Derivatives in gamma(i)
if (calc_grad) then
+ if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+ & (shield_mode.gt.0)) then
+C print *,i,j
+
+ do ilist=1,ishield_list(i)
+ iresshield=shield_list(ilist,i)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,i)*eello_t4/fac_shield(i)
+C & *2.0
+ gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,i)*eello_t4/fac_shield(i)
+ gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1)
+ & +rlocshield
+ enddo
+ enddo
+ do ilist=1,ishield_list(j)
+ iresshield=shield_list(ilist,j)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,j)*eello_t4/fac_shield(j)
+C & *2.0
+ gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(k,ilist,j)*eello_t4/fac_shield(j)
+ gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1)
+ & +rlocshield
+
+ enddo
+ enddo
+
+ do k=1,3
+ gshieldc_t4(k,i)=gshieldc_t4(k,i)+
+ & grad_shield(k,i)*eello_t4/fac_shield(i)
+ gshieldc_t4(k,j)=gshieldc_t4(k,j)+
+ & grad_shield(k,j)*eello_t4/fac_shield(j)
+ gshieldc_t4(k,i-1)=gshieldc_t4(k,i-1)+
+ & grad_shield(k,i)*eello_t4/fac_shield(i)
+ gshieldc_t4(k,j-1)=gshieldc_t4(k,j-1)+
+ & grad_shield(k,j)*eello_t4/fac_shield(j)
+ enddo
+ endif
+
call transpose2(EUgder(1,1,i+1),e1tder(1,1))
call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1))
call matvec2(auxmat(1,1),Ub2(1,i+3),auxvec(1))
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)
+
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)
+
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(auxmat(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)
+
C Cartesian derivatives
C Derivatives of this turn contributions in DC(i+2)
if (j.lt.nres-1) then
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)
+
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)
+
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)
+
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)
+
a_temp(1,1)=aggj1(l,1)
a_temp(1,2)=aggj1(l,2)
a_temp(2,1)=aggj1(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,j1)=gcorr4_turn(l,j1)-(s1+s2+s3)
+ & *fac_shield(i)*fac_shield(j)
+
enddo
endif
+ 178 continue
endif
return
end
c write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e,
c & ' scal14',scal14
do i=iatscp_s,iatscp_e
- if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
iteli=itel(i)
c write (iout,*) "i",i," iteli",iteli," nscp_gr",nscp_gr(i),
c & " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
xi=0.5D0*(c(1,i)+c(1,i+1))
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
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
- itypj=itype(j)
- if (itypj.eq.21) cycle
+ itypj=iabs(itype(j))
+ if (itypj.eq.ntyp1) cycle
C Uncomment following three lines for SC-p interactions
c xj=c(1,nres+j)-xi
c yj=c(2,nres+j)-yi
c zj=c(3,nres+j)-zi
C Uncomment following three lines for Ca-p interactions
- xj=c(1,j)-xi
- yj=c(2,j)-yi
- zj=c(3,j)-zi
+ xj=c(1,j)
+ 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
+
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
+ sss=sscale(1.0d0/(dsqrt(rrij)))
+ if (sss.le.0.0d0) cycle
+ sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
fac=rrij**expon2
e1=fac*fac*aad(itypj,iteli)
e2=fac*bad(itypj,iteli)
if (iabs(j-i) .le. 2) then
e1=scal14*e1
e2=scal14*e2
- evdw2_14=evdw2_14+e1+e2
+ evdw2_14=evdw2_14+(e1+e2)*sss
endif
evdwij=e1+e2
c write (iout,*) i,j,evdwij
- evdw2=evdw2+evdwij
+ evdw2=evdw2+evdwij*sss
if (calc_grad) then
C
C Calculate contributions to the gradient in the virtual-bond and SC vectors.
C
- fac=-(evdwij+e1)*rrij
+ fac=-(evdwij+e1)*rrij*sss
+ fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/expon
ggg(1)=xj*fac
ggg(2)=yj*fac
ggg(3)=zj*fac
include 'COMMON.DERIV'
include 'COMMON.VAR'
include 'COMMON.INTERACT'
+ include 'COMMON.CONTROL'
dimension ggg(3)
ehpb=0.0D0
cd print *,'edis: nhpb=',nhpb,' fbr=',fbr
endif
C 24/11/03 AL: SS bridges handled separately because of introducing a specific
C distance and angle dependent SS bond potential.
- if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+C if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
+C & iabs(itype(jjj)).eq.1) then
+C call ssbond_ene(iii,jjj,eij)
+C ehpb=ehpb+2*eij
+C else
+ if (.not.dyn_ss .and. i.le.nss) then
+ if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
+ & iabs(itype(jjj)).eq.1) then
call ssbond_ene(iii,jjj,eij)
ehpb=ehpb+2*eij
- else
+ endif !ii.gt.neres
+ else if (ii.gt.nres .and. jj.gt.nres) then
+c Restraints from contact prediction
+ dd=dist(ii,jj)
+ if (constr_dist.eq.11) then
+C ehpb=ehpb+fordepth(i)**4.0d0
+C & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+ ehpb=ehpb+fordepth(i)**4.0d0
+ & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+ fac=fordepth(i)**4.0d0
+ & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj,
+C & ehpb,fordepth(i),dd
+C print *,"TUTU"
+C write(iout,*) ehpb,"atu?"
+C ehpb,"tu?"
+C fac=fordepth(i)**4.0d0
+C & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+ else !constr_dist.eq.11
+ if (dhpb1(i).gt.0.0d0) then
+ ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c write (iout,*) "beta nmr",
+c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ else !dhpb(i).gt.0.00
+
C Calculate the distance between the two points and its difference from the
C target distance.
dd=dist(ii,jj)
C Evaluate gradient.
C
fac=waga*rdis/dd
+ endif !dhpb(i).gt.0
+ endif
cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
cd & ' waga=',waga,' fac=',fac
do j=1,3
ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
enddo
endif
+ else !ii.gt.nres
+C write(iout,*) "before"
+ dd=dist(ii,jj)
+C write(iout,*) "after",dd
+ if (constr_dist.eq.11) then
+ ehpb=ehpb+fordepth(i)**4.0d0
+ & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+ fac=fordepth(i)**4.0d0
+ & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+C ehpb=ehpb+fordepth(i)**4*rlornmr1(dd,dhpb(i),dhpb1(i))
+C fac=fordepth(i)**4*rlornmr1prim(dd,dhpb(i),dhpb1(i))/dd
+C print *,ehpb,"tu?"
+C write(iout,*) ehpb,"btu?",
+C & dd,dhpb(i),dhpb1(i),fordepth(i),forcon(i)
+C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj,
+C & ehpb,fordepth(i),dd
+ else
+ if (dhpb1(i).gt.0.0d0) then
+ ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c write (iout,*) "alph nmr",
+c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ else
+ rdis=dd-dhpb(i)
+C Get the force constant corresponding to this distance.
+ waga=forcon(i)
+C Calculate the contribution to energy.
+ ehpb=ehpb+waga*rdis*rdis
+c write (iout,*) "alpha reg",dd,waga*rdis*rdis
+C
+C Evaluate gradient.
+C
+ fac=waga*rdis/dd
+ endif
+ endif
+ do j=1,3
+ ggg(j)=fac*(c(j,jj)-c(j,ii))
+ enddo
+cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
+C If this is a SC-SC distance, we need to calculate the contributions to the
+C Cartesian gradient in the SC vectors (ghpbx).
+ if (iii.lt.ii) then
+ do j=1,3
+ ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
+ ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
+ enddo
+ endif
do j=iii,jjj-1
do k=1,3
ghpbc(k,j)=ghpbc(k,j)+ggg(k)
enddo
endif
enddo
- ehpb=0.5D0*ehpb
+ if (constr_dist.ne.11) ehpb=0.5D0*ehpb
return
end
C--------------------------------------------------------------------------
include 'COMMON.VAR'
include 'COMMON.IOUNITS'
double precision erij(3),dcosom1(3),dcosom2(3),gg(3)
- itypi=itype(i)
+ itypi=iabs(itype(i))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
dsci_inv=dsc_inv(itypi)
- itypj=itype(j)
+ itypj=iabs(itype(j))
dscj_inv=dsc_inv(itypj)
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
do k=1,3
gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
enddo
- do k=1,3
- ghpbx(k,i)=ghpbx(k,i)-gg(k)
- & +(eom12*dc_norm(k,nres+j)+eom1*erij(k))*dsci_inv
- ghpbx(k,j)=ghpbx(k,j)+gg(k)
- & +(eom12*dc_norm(k,nres+i)+eom2*erij(k))*dscj_inv
+ do k=1,3
+ ghpbx(k,i)=ghpbx(k,i)-gg(k)
+ & +(eom12*dc_norm(k,nres+j)+eom1*erij(k))*dsci_inv
+ ghpbx(k,j)=ghpbx(k,j)+gg(k)
+ & +(eom12*dc_norm(k,nres+i)+eom2*erij(k))*dscj_inv
+ enddo
+C
+C Calculate the components of the gradient in DC and X
+C
+ do k=i,j-1
+ do l=1,3
+ ghpbc(l,k)=ghpbc(l,k)+gg(l)
+ enddo
+ enddo
+ return
+ end
+C--------------------------------------------------------------------------
+
+
+c LICZENIE WIEZOW Z ROWNANIA ENERGII MODELLERA
+ subroutine e_modeller(ehomology_constr)
+ implicit real*8 (a-h,o-z)
+
+ include 'DIMENSIONS'
+
+ integer nnn, i, j, k, ki, irec, l
+ integer katy, odleglosci, test7
+ real*8 odleg, odleg2, odleg3, kat, kat2, kat3, gdih(max_template)
+ real*8 distance(max_template),distancek(max_template),
+ & min_odl,godl(max_template),dih_diff(max_template)
+
+c
+c FP - 30/10/2014 Temporary specifications for homology restraints
+c
+ double precision utheta_i,gutheta_i,sum_gtheta,sum_sgtheta,
+ & sgtheta
+ double precision, dimension (maxres) :: guscdiff,usc_diff
+ double precision, dimension (max_template) ::
+ & gtheta,dscdiff,uscdiffk,guscdiff2,guscdiff3,
+ & theta_diff
+
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.CHAIN'
+ include 'COMMON.GEO'
+ include 'COMMON.DERIV'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.HOMRESTR'
+c
+ include 'COMMON.SETUP'
+ include 'COMMON.NAMES'
+
+ do i=1,max_template
+ distancek(i)=9999999.9
+ enddo
+
+ odleg=0.0d0
+
+c Pseudo-energy and gradient from homology restraints (MODELLER-like
+c function)
+C AL 5/2/14 - Introduce list of restraints
+c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
+#ifdef DEBUG
+ write(iout,*) "------- dist restrs start -------"
+ write (iout,*) "link_start_homo",link_start_homo,
+ & " link_end_homo",link_end_homo
+#endif
+ do ii = link_start_homo,link_end_homo
+ i = ires_homo(ii)
+ j = jres_homo(ii)
+ dij=dist(i,j)
+c write (iout,*) "dij(",i,j,") =",dij
+ do k=1,constr_homology
+ if(.not.l_homo(k,ii)) cycle
+ distance(k)=odl(k,ii)-dij
+c write (iout,*) "distance(",k,") =",distance(k)
+c
+c For Gaussian-type Urestr
+c
+ distancek(k)=0.5d0*distance(k)**2*sigma_odl(k,ii) ! waga_dist rmvd from Gaussian argument
+c write (iout,*) "sigma_odl(",k,ii,") =",sigma_odl(k,ii)
+c write (iout,*) "distancek(",k,") =",distancek(k)
+c distancek(k)=0.5d0*waga_dist*distance(k)**2*sigma_odl(k,ii)
+c
+c For Lorentzian-type Urestr
+c
+ if (waga_dist.lt.0.0d0) then
+ sigma_odlir(k,ii)=dsqrt(1/sigma_odl(k,ii))
+ distancek(k)=distance(k)**2/(sigma_odlir(k,ii)*
+ & (distance(k)**2+sigma_odlir(k,ii)**2))
+ endif
+ 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)
+ & min_odl=distancek(kk)
+ enddo
+c write (iout,* )"min_odl",min_odl
+#ifdef DEBUG
+ write (iout,*) "ij dij",i,j,dij
+ write (iout,*) "distance",(distance(k),k=1,constr_homology)
+ write (iout,*) "distancek",(distancek(k),k=1,constr_homology)
+ write (iout,* )"min_odl",min_odl
+#endif
+ odleg2=0.0d0
+ do k=1,constr_homology
+c Nie wiem po co to liczycie jeszcze raz!
+c odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/
+c & (2*(sigma_odl(i,j,k))**2))
+ if(.not.l_homo(k,ii)) cycle
+ if (waga_dist.ge.0.0d0) then
+c
+c For Gaussian-type Urestr
+c
+ godl(k)=dexp(-distancek(k)+min_odl)
+ odleg2=odleg2+godl(k)
+c
+c For Lorentzian-type Urestr
+c
+ else
+ odleg2=odleg2+distancek(k)
+ endif
+
+ccc write(iout,779) i,j,k, "odleg2=",odleg2, "odleg3=", odleg3,
+ccc & "dEXP(odleg3)=", dEXP(odleg3),"distance(i,j,k)^2=",
+ccc & distance(i,j,k)**2, "dist(i+1,j+1)=", dist(i+1,j+1),
+ccc & "sigma_odl(i,j,k)=", sigma_odl(i,j,k)
+
+ enddo
+c write (iout,*) "godl",(godl(k),k=1,constr_homology) ! exponents
+c write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps
+#ifdef DEBUG
+ write (iout,*) "godl",(godl(k),k=1,constr_homology) ! exponents
+ write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps
+#endif
+ if (waga_dist.ge.0.0d0) then
+c
+c For Gaussian-type Urestr
+c
+ odleg=odleg-dLOG(odleg2/constr_homology)+min_odl
+c
+c For Lorentzian-type Urestr
+c
+ else
+ odleg=odleg+odleg2/constr_homology
+ endif
+c
+#ifdef GRAD
+c write (iout,*) "odleg",odleg ! sum of -ln-s
+c Gradient
+c
+c For Gaussian-type Urestr
+c
+ if (waga_dist.ge.0.0d0) sum_godl=odleg2
+ sum_sgodl=0.0d0
+ do k=1,constr_homology
+c godl=dexp(((-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2))
+c & *waga_dist)+min_odl
+c sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist
+c
+ if(.not.l_homo(k,ii)) cycle
+ if (waga_dist.ge.0.0d0) then
+c For Gaussian-type Urestr
+c
+ sgodl=-godl(k)*distance(k)*sigma_odl(k,ii) ! waga_dist rmvd
+c
+c For Lorentzian-type Urestr
+c
+ else
+ sgodl=-2*sigma_odlir(k,ii)*(distance(k)/(distance(k)**2+
+ & sigma_odlir(k,ii)**2)**2)
+ endif
+ sum_sgodl=sum_sgodl+sgodl
+
+c sgodl2=sgodl2+sgodl
+c write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE1"
+c write(iout,*) "constr_homology=",constr_homology
+c write(iout,*) i, j, k, "TEST K"
+ enddo
+ if (waga_dist.ge.0.0d0) then
+c
+c For Gaussian-type Urestr
+c
+ grad_odl3=waga_homology(iset)*waga_dist
+ & *sum_sgodl/(sum_godl*dij)
+c
+c For Lorentzian-type Urestr
+c
+ else
+c Original grad expr modified by analogy w Gaussian-type Urestr grad
+c grad_odl3=-waga_homology(iset)*waga_dist*sum_sgodl
+ grad_odl3=-waga_homology(iset)*waga_dist*
+ & sum_sgodl/(constr_homology*dij)
+ endif
+c
+c grad_odl3=sum_sgodl/(sum_godl*dij)
+
+
+c write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE2"
+c write(iout,*) (distance(i,j,k)**2), (2*(sigma_odl(i,j,k))**2),
+c & (-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2))
+
+ccc write(iout,*) godl, sgodl, grad_odl3
+
+c grad_odl=grad_odl+grad_odl3
+
+ do jik=1,3
+ ggodl=grad_odl3*(c(jik,i)-c(jik,j))
+ccc write(iout,*) c(jik,i+1), c(jik,j+1), (c(jik,i+1)-c(jik,j+1))
+ccc write(iout,746) "GRAD_ODL_1", i, j, jik, ggodl,
+ccc & ghpbc(jik,i+1), ghpbc(jik,j+1)
+ ghpbc(jik,i)=ghpbc(jik,i)+ggodl
+ ghpbc(jik,j)=ghpbc(jik,j)-ggodl
+ccc write(iout,746) "GRAD_ODL_2", i, j, jik, ggodl,
+ccc & ghpbc(jik,i+1), ghpbc(jik,j+1)
+c if (i.eq.25.and.j.eq.27) then
+c write(iout,*) "jik",jik,"i",i,"j",j
+c write(iout,*) "sum_sgodl",sum_sgodl,"sgodl",sgodl
+c write(iout,*) "grad_odl3",grad_odl3
+c write(iout,*) "c(",jik,i,")",c(jik,i),"c(",jik,j,")",c(jik,j)
+c write(iout,*) "ggodl",ggodl
+c write(iout,*) "ghpbc(",jik,i,")",
+c & ghpbc(jik,i),"ghpbc(",jik,j,")",
+c & ghpbc(jik,j)
+c endif
+ enddo
+#endif
+ccc write(iout,778)"TEST: odleg2=", odleg2, "DLOG(odleg2)=",
+ccc & dLOG(odleg2),"-odleg=", -odleg
+
+ enddo ! ii-loop for dist
+#ifdef DEBUG
+ write(iout,*) "------- dist restrs end -------"
+c if (waga_angle.eq.1.0d0 .or. waga_theta.eq.1.0d0 .or.
+c & waga_d.eq.1.0d0) call sum_gradient
+#endif
+c Pseudo-energy and gradient from dihedral-angle restraints from
+c homology templates
+c write (iout,*) "End of distance loop"
+c call flush(iout)
+ kat=0.0d0
+c write (iout,*) idihconstr_start_homo,idihconstr_end_homo
+#ifdef DEBUG
+ write(iout,*) "------- dih restrs start -------"
+ do i=idihconstr_start_homo,idihconstr_end_homo
+ write (iout,*) "gloc_init(",i,icg,")",gloc(i,icg)
+ enddo
+#endif
+ do i=idihconstr_start_homo,idihconstr_end_homo
+ kat2=0.0d0
+c betai=beta(i,i+1,i+2,i+3)
+ betai = phi(i+3)
+c write (iout,*) "betai =",betai
+ do k=1,constr_homology
+ dih_diff(k)=pinorm(dih(k,i)-betai)
+c write (iout,*) "dih_diff(",k,") =",dih_diff(k)
+c if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)=
+c & -(6.28318-dih_diff(i,k))
+c if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)=
+c & 6.28318+dih_diff(i,k)
+
+ kat3=-0.5d0*dih_diff(k)**2*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument
+c kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i)
+ gdih(k)=dexp(kat3)
+ kat2=kat2+gdih(k)
+c write(iout,*) "kat2=", kat2, "exp(kat3)=", exp(kat3)
+c write(*,*)""
+ enddo
+c write (iout,*) "gdih",(gdih(k),k=1,constr_homology) ! exps
+c write (iout,*) "i",i," betai",betai," kat2",kat2 ! sum of exps
+#ifdef DEBUG
+ write (iout,*) "i",i," betai",betai," kat2",kat2
+ write (iout,*) "gdih",(gdih(k),k=1,constr_homology)
+#endif
+ if (kat2.le.1.0d-14) cycle
+ kat=kat-dLOG(kat2/constr_homology)
+c write (iout,*) "kat",kat ! sum of -ln-s
+
+ccc write(iout,778)"TEST: kat2=", kat2, "DLOG(kat2)=",
+ccc & dLOG(kat2), "-kat=", -kat
+
+#ifdef GRAD
+c ----------------------------------------------------------------------
+c Gradient
+c ----------------------------------------------------------------------
+
+ sum_gdih=kat2
+ sum_sgdih=0.0d0
+ do k=1,constr_homology
+ sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i) ! waga_angle rmvd
+c sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle
+ sum_sgdih=sum_sgdih+sgdih
+ enddo
+c grad_dih3=sum_sgdih/sum_gdih
+ grad_dih3=waga_homology(iset)*waga_angle*sum_sgdih/sum_gdih
+
+c write(iout,*)i,k,gdih,sgdih,beta(i+1,i+2,i+3,i+4),grad_dih3
+ccc write(iout,747) "GRAD_KAT_1", i, nphi, icg, grad_dih3,
+ccc & gloc(nphi+i-3,icg)
+ gloc(i,icg)=gloc(i,icg)+grad_dih3
+c if (i.eq.25) then
+c write(iout,*) "i",i,"icg",icg,"gloc(",i,icg,")",gloc(i,icg)
+c endif
+ccc write(iout,747) "GRAD_KAT_2", i, nphi, icg, grad_dih3,
+ccc & gloc(nphi+i-3,icg)
+#endif
+ enddo ! i-loop for dih
+#ifdef DEBUG
+ write(iout,*) "------- dih restrs end -------"
+#endif
+
+c Pseudo-energy and gradient for theta angle restraints from
+c homology templates
+c FP 01/15 - inserted from econstr_local_test.F, loop structure
+c adapted
+
+c
+c For constr_homology reference structures (FP)
+c
+c Uconst_back_tot=0.0d0
+ Eval=0.0d0
+ Erot=0.0d0
+c Econstr_back legacy
+#ifdef GRAD
+ do i=1,nres
+c do i=ithet_start,ithet_end
+ dutheta(i)=0.0d0
+c enddo
+c do i=loc_start,loc_end
+ do j=1,3
+ duscdiff(j,i)=0.0d0
+ duscdiffx(j,i)=0.0d0
+ enddo
+ enddo
+#endif
+c
+c do iref=1,nref
+c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
+c write (iout,*) "waga_theta",waga_theta
+ if (waga_theta.gt.0.0d0) then
+#ifdef DEBUG
+ write (iout,*) "usampl",usampl
+ write(iout,*) "------- theta restrs start -------"
+c do i=ithet_start,ithet_end
+c write (iout,*) "gloc_init(",nphi+i,icg,")",gloc(nphi+i,icg)
+c enddo
+#endif
+c write (iout,*) "maxres",maxres,"nres",nres
+
+ do i=ithet_start,ithet_end
+c
+c do i=1,nfrag_back
+c ii = ifrag_back(2,i,iset)-ifrag_back(1,i,iset)
+c
+c Deviation of theta angles wrt constr_homology ref structures
+c
+ utheta_i=0.0d0 ! argument of Gaussian for single k
+ gutheta_i=0.0d0 ! Sum of Gaussians over constr_homology ref structures
+c do j=ifrag_back(1,i,iset)+2,ifrag_back(2,i,iset) ! original loop
+c over residues in a fragment
+c write (iout,*) "theta(",i,")=",theta(i)
+ do k=1,constr_homology
+c
+c dtheta_i=theta(j)-thetaref(j,iref)
+c dtheta_i=thetaref(k,i)-theta(i) ! original form without indexing
+ theta_diff(k)=thetatpl(k,i)-theta(i)
+c
+ utheta_i=-0.5d0*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta rmvd from Gaussian argument
+c utheta_i=-0.5d0*waga_theta*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta?
+ gtheta(k)=dexp(utheta_i) ! + min_utheta_i?
+ gutheta_i=gutheta_i+dexp(utheta_i) ! Sum of Gaussians (pk)
+c Gradient for single Gaussian restraint in subr Econstr_back
+c dutheta(j-2)=dutheta(j-2)+wfrag_back(1,i,iset)*dtheta_i/(ii-1)
+c
+ enddo
+c write (iout,*) "gtheta",(gtheta(k),k=1,constr_homology) ! exps
+c write (iout,*) "i",i," gutheta_i",gutheta_i ! sum of exps
+
+c
+#ifdef GRAD
+c Gradient for multiple Gaussian restraint
+ sum_gtheta=gutheta_i
+ sum_sgtheta=0.0d0
+ do k=1,constr_homology
+c New generalized expr for multiple Gaussian from Econstr_back
+ sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i) ! waga_theta rmvd
+c
+c sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i)*waga_theta ! right functional form?
+ sum_sgtheta=sum_sgtheta+sgtheta ! cum variable
+ enddo
+c
+c Final value of gradient using same var as in Econstr_back
+ dutheta(i-2)=sum_sgtheta/sum_gtheta*waga_theta
+ & *waga_homology(iset)
+c dutheta(i)=sum_sgtheta/sum_gtheta
+c
+c Uconst_back=Uconst_back+waga_theta*utheta(i) ! waga_theta added as weight
+#endif
+ Eval=Eval-dLOG(gutheta_i/constr_homology)
+c write (iout,*) "utheta(",i,")=",utheta(i) ! -ln of sum of exps
+c write (iout,*) "Uconst_back",Uconst_back ! sum of -ln-s
+c Uconst_back=Uconst_back+utheta(i)
+ enddo ! (i-loop for theta)
+#ifdef DEBUG
+ write(iout,*) "------- theta restrs end -------"
+#endif
+ endif
+c
+c Deviation of local SC geometry
+c
+c Separation of two i-loops (instructed by AL - 11/3/2014)
+c
+c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
+c write (iout,*) "waga_d",waga_d
+
+#ifdef DEBUG
+ write(iout,*) "------- SC restrs start -------"
+ write (iout,*) "Initial duscdiff,duscdiffx"
+ do i=loc_start,loc_end
+ write (iout,*) i,(duscdiff(jik,i),jik=1,3),
+ & (duscdiffx(jik,i),jik=1,3)
enddo
-C
-C Calculate the components of the gradient in DC and X
-C
- do k=i,j-1
- do l=1,3
- ghpbc(l,k)=ghpbc(l,k)+gg(l)
+#endif
+ do i=loc_start,loc_end
+ usc_diff_i=0.0d0 ! argument of Gaussian for single k
+ guscdiff(i)=0.0d0 ! Sum of Gaussians over constr_homology ref structures
+c do j=ifrag_back(1,i,iset)+1,ifrag_back(2,i,iset)-1 ! Econstr_back legacy
+c write(iout,*) "xxtab, yytab, zztab"
+c write(iout,'(i5,3f8.2)') i,xxtab(i),yytab(i),zztab(i)
+ do k=1,constr_homology
+c
+ dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str?
+c Original sign inverted for calc of gradients (s. Econstr_back)
+ dyy=-yytpl(k,i)+yytab(i) ! ibid y
+ dzz=-zztpl(k,i)+zztab(i) ! ibid z
+c write(iout,*) "dxx, dyy, dzz"
+c write(iout,'(2i5,3f8.2)') k,i,dxx,dyy,dzz
+c
+ usc_diff_i=-0.5d0*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d rmvd from Gaussian argument
+c usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d?
+c uscdiffk(k)=usc_diff(i)
+ guscdiff2(k)=dexp(usc_diff_i) ! without min_scdiff
+ guscdiff(i)=guscdiff(i)+dexp(usc_diff_i) !Sum of Gaussians (pk)
+c write (iout,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j),
+c & xxref(j),yyref(j),zzref(j)
+ enddo
+c
+c Gradient
+c
+c Generalized expression for multiple Gaussian acc to that for a single
+c Gaussian in Econstr_back as instructed by AL (FP - 03/11/2014)
+c
+c Original implementation
+c sum_guscdiff=guscdiff(i)
+c
+c sum_sguscdiff=0.0d0
+c do k=1,constr_homology
+c sguscdiff=-guscdiff2(k)*dscdiff(k)*sigma_d(k,i)*waga_d !waga_d?
+c sguscdiff=-guscdiff3(k)*dscdiff(k)*sigma_d(k,i)*waga_d ! w min_uscdiff
+c sum_sguscdiff=sum_sguscdiff+sguscdiff
+c enddo
+c
+c Implementation of new expressions for gradient (Jan. 2015)
+c
+c grad_uscdiff=sum_sguscdiff/(sum_guscdiff*dtab) !?
+#ifdef GRAD
+ do k=1,constr_homology
+c
+c New calculation of dxx, dyy, and dzz corrected by AL (07/11), was missing and wrong
+c before. Now the drivatives should be correct
+c
+ dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str?
+c Original sign inverted for calc of gradients (s. Econstr_back)
+ dyy=-yytpl(k,i)+yytab(i) ! ibid y
+ dzz=-zztpl(k,i)+zztab(i) ! ibid z
+c
+c New implementation
+c
+ sum_guscdiff=guscdiff2(k)*!(dsqrt(dxx*dxx+dyy*dyy+dzz*dzz))* -> wrong!
+ & sigma_d(k,i) ! for the grad wrt r'
+c sum_sguscdiff=sum_sguscdiff+sum_guscdiff
+c
+c
+c New implementation
+ sum_guscdiff = waga_homology(iset)*waga_d*sum_guscdiff
+ do jik=1,3
+ duscdiff(jik,i-1)=duscdiff(jik,i-1)+
+ & sum_guscdiff*(dXX_C1tab(jik,i)*dxx+
+ & dYY_C1tab(jik,i)*dyy+dZZ_C1tab(jik,i)*dzz)/guscdiff(i)
+ duscdiff(jik,i)=duscdiff(jik,i)+
+ & sum_guscdiff*(dXX_Ctab(jik,i)*dxx+
+ & dYY_Ctab(jik,i)*dyy+dZZ_Ctab(jik,i)*dzz)/guscdiff(i)
+ duscdiffx(jik,i)=duscdiffx(jik,i)+
+ & sum_guscdiff*(dXX_XYZtab(jik,i)*dxx+
+ & dYY_XYZtab(jik,i)*dyy+dZZ_XYZtab(jik,i)*dzz)/guscdiff(i)
+c
+#ifdef DEBUG
+ write(iout,*) "jik",jik,"i",i
+ write(iout,*) "dxx, dyy, dzz"
+ write(iout,'(2i5,3f8.2)') k,i,dxx,dyy,dzz
+ write(iout,*) "guscdiff2(",k,")",guscdiff2(k)
+c write(iout,*) "sum_sguscdiff",sum_sguscdiff
+cc write(iout,*) "dXX_Ctab(",jik,i,")",dXX_Ctab(jik,i)
+c write(iout,*) "dYY_Ctab(",jik,i,")",dYY_Ctab(jik,i)
+c write(iout,*) "dZZ_Ctab(",jik,i,")",dZZ_Ctab(jik,i)
+c write(iout,*) "dXX_C1tab(",jik,i,")",dXX_C1tab(jik,i)
+c write(iout,*) "dYY_C1tab(",jik,i,")",dYY_C1tab(jik,i)
+c write(iout,*) "dZZ_C1tab(",jik,i,")",dZZ_C1tab(jik,i)
+c write(iout,*) "dXX_XYZtab(",jik,i,")",dXX_XYZtab(jik,i)
+c write(iout,*) "dYY_XYZtab(",jik,i,")",dYY_XYZtab(jik,i)
+c write(iout,*) "dZZ_XYZtab(",jik,i,")",dZZ_XYZtab(jik,i)
+c write(iout,*) "duscdiff(",jik,i-1,")",duscdiff(jik,i-1)
+c write(iout,*) "duscdiff(",jik,i,")",duscdiff(jik,i)
+c write(iout,*) "duscdiffx(",jik,i,")",duscdiffx(jik,i)
+c endif
+#endif
+ enddo
+ enddo
+#endif
+c
+c uscdiff(i)=-dLOG(guscdiff(i)/(ii-1)) ! Weighting by (ii-1) required?
+c usc_diff(i)=-dLOG(guscdiff(i)/constr_homology) ! + min_uscdiff ?
+c
+c write (iout,*) i," uscdiff",uscdiff(i)
+c
+c Put together deviations from local geometry
+
+c Uconst_back=Uconst_back+wfrag_back(1,i,iset)*utheta(i)+
+c & wfrag_back(3,i,iset)*uscdiff(i)
+ Erot=Erot-dLOG(guscdiff(i)/constr_homology)
+c write (iout,*) "usc_diff(",i,")=",usc_diff(i) ! -ln of sum of exps
+c write (iout,*) "Uconst_back",Uconst_back ! cum sum of -ln-s
+c Uconst_back=Uconst_back+usc_diff(i)
+c
+c Gradient of multiple Gaussian restraint (FP - 04/11/2014 - right?)
+c
+c New implment: multiplied by sum_sguscdiff
+c
+
+ enddo ! (i-loop for dscdiff)
+
+c endif
+
+#ifdef DEBUG
+ write(iout,*) "------- SC restrs end -------"
+ write (iout,*) "------ After SC loop in e_modeller ------"
+ do i=loc_start,loc_end
+ write (iout,*) "i",i," gradc",(gradc(j,i,icg),j=1,3)
+ write (iout,*) "i",i," gradx",(gradx(j,i,icg),j=1,3)
enddo
+ if (waga_theta.eq.1.0d0) then
+ write (iout,*) "in e_modeller after SC restr end: dutheta"
+ do i=ithet_start,ithet_end
+ write (iout,*) i,dutheta(i)
+ enddo
+ endif
+ if (waga_d.eq.1.0d0) then
+ write (iout,*) "e_modeller after SC loop: duscdiff/x"
+ do i=1,nres
+ write (iout,*) i,(duscdiff(j,i),j=1,3)
+ write (iout,*) i,(duscdiffx(j,i),j=1,3)
enddo
+ endif
+#endif
+
+c Total energy from homology restraints
+#ifdef DEBUG
+ write (iout,*) "odleg",odleg," kat",kat
+ write (iout,*) "odleg",odleg," kat",kat
+ write (iout,*) "Eval",Eval," Erot",Erot
+ write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
+ write (iout,*) "waga_dist ",waga_dist,"waga_angle ",waga_angle
+ write (iout,*) "waga_theta ",waga_theta,"waga_d ",waga_d
+ write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
+#endif
+c
+c Addition of energy of theta angle and SC local geom over constr_homologs ref strs
+c
+c ehomology_constr=odleg+kat
+c
+c For Lorentzian-type Urestr
+c
+
+ if (waga_dist.ge.0.0d0) then
+c
+c For Gaussian-type Urestr
+c
+ ehomology_constr=(waga_dist*odleg+waga_angle*kat+
+ & waga_theta*Eval+waga_d*Erot)*waga_homology(iset)
+c write (iout,*) "ehomology_constr=",ehomology_constr
+ else
+c
+c For Lorentzian-type Urestr
+c
+ ehomology_constr=(-waga_dist*odleg+waga_angle*kat+
+ & waga_theta*Eval+waga_d*Erot)*waga_homology(iset)
+c write (iout,*) "ehomology_constr=",ehomology_constr
+ endif
+#ifdef DEBUG
+ write (iout,*) "iset",iset," waga_homology",waga_homology(iset)
+ write (iout,*) "odleg",waga_dist,odleg," kat",waga_angle,kat,
+ & " Eval",waga_theta,Eval," Erot",waga_d,Erot
+ write (iout,*) "ehomology_constr",ehomology_constr
+#endif
return
+
+ 748 format(a8,f12.3,a6,f12.3,a7,f12.3)
+ 747 format(a12,i4,i4,i4,f8.3,f8.3)
+ 746 format(a12,i4,i4,i4,f8.3,f8.3,f8.3)
+ 778 format(a7,1X,f10.3,1X,a4,1X,f10.3,1X,a5,1X,f10.3)
+ 779 format(i3,1X,i3,1X,i2,1X,a7,1X,f7.3,1X,a7,1X,f7.3,1X,a13,1X,
+ & f7.3,1X,a17,1X,f9.3,1X,a10,1X,f8.3,1X,a10,1X,f8.3)
end
C--------------------------------------------------------------------------
+
+C--------------------------------------------------------------------------
subroutine ebond(estr)
c
c Evaluate the energy of stretching of the CA-CA and CA-SC virtual bonds
logical energy_dec /.false./
double precision u(3),ud(3)
estr=0.0d0
+ estr1=0.0d0
do i=nnt+1,nct
- if (itype(i-1).eq.21 .or. itype(i).eq.21) then
- estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
- do j=1,3
- gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax)
- & *dc(j,i-1)/vbld(i)
- enddo
- if (energy_dec) write(iout,*)
- & "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
- 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
+C gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax)
+C & *dc(j,i-1)/vbld(i)
+C enddo
+C if (energy_dec) write(iout,*)
+C & "estr1",i,vbld(i),distchainmax,
+C & gnmr1(vbld(i),-1.0d0,distchainmax)
+C else
+ if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
+ diff = vbld(i)-vbldpDUM
+ else
diff = vbld(i)-vbldp0
c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff
+ endif
estr=estr+diff*diff
do j=1,3
gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
enddo
- endif
-
+C endif
+C write (iout,'(a7,i5,4f7.3)')
+C & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
enddo
estr=0.5d0*AKP*estr+estr1
c
c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
c
do i=nnt,nct
- iti=itype(i)
- if (iti.ne.10 .and. iti.ne.21) then
+ iti=iabs(itype(i))
+ if (iti.ne.10 .and. iti.ne.ntyp1) then
nbi=nbondterm(iti)
if (nbi.eq.1) then
diff=vbld(i+nres)-vbldsc0(1,iti)
end
#ifdef CRYST_THETA
C--------------------------------------------------------------------------
- subroutine ebend(etheta)
+ subroutine ebend(etheta,ethetacnstr)
C
C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
C angles gamma and its derivatives in consecutive thetas and gammas.
include 'COMMON.IOUNITS'
include 'COMMON.NAMES'
include 'COMMON.FFIELD'
+ include 'COMMON.TORCNSTR'
common /calcthet/ term1,term2,termm,diffak,ratak,
& ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq,
& delthe0,sig0inv,sigtc,sigsqtc,delthec,it
double precision y(2),z(2)
delta=0.02d0*pi
- time11=dexp(-2*time)
- time12=1.0d0
+c time11=dexp(-2*time)
+c time12=1.0d0
etheta=0.0D0
c write (iout,*) "nres",nres
c write (*,'(a,i2)') 'EBEND ICG=',icg
c write (iout,*) ithet_start,ithet_end
do i=ithet_start,ithet_end
- if (itype(i-1).eq.21) cycle
+ if (i.le.2) cycle
+ if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
+ & .or.itype(i).eq.ntyp1) cycle
C Zero the energy function and its derivative at 0 or pi.
call splinthet(theta(i),0.5d0*delta,ss,ssd)
it=itype(i-1)
- if (i.gt.3 .and. itype(i-2).ne.21) then
+ ichir1=isign(1,itype(i-2))
+ ichir2=isign(1,itype(i))
+ if (itype(i-2).eq.10) ichir1=isign(1,itype(i-1))
+ if (itype(i).eq.10) ichir2=isign(1,itype(i-1))
+ if (itype(i-1).eq.10) then
+ itype1=isign(10,itype(i-2))
+ ichir11=isign(1,itype(i-2))
+ ichir12=isign(1,itype(i-2))
+ itype2=isign(10,itype(i))
+ ichir21=isign(1,itype(i))
+ ichir22=isign(1,itype(i))
+ endif
+ if (i.eq.3) then
+ y(1)=0.0D0
+ y(2)=0.0D0
+ else
+ if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
- icrc=0
- call proc_proc(phii,icrc)
+c icrc=0
+c call proc_proc(phii,icrc)
if (icrc.eq.1) phii=150.0
#else
phii=phi(i)
y(1)=0.0D0
y(2)=0.0D0
endif
- if (i.lt.nres .and. itype(i).ne.21) then
+ endif
+ if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
- icrc=0
- call proc_proc(phii1,icrc)
+c icrc=0
+c call proc_proc(phii1,icrc)
if (icrc.eq.1) phii1=150.0
phii1=pinorm(phii1)
z(1)=cos(phii1)
C In following comments this theta will be referred to as t_c.
thet_pred_mean=0.0d0
do k=1,2
- athetk=athet(k,it)
- bthetk=bthet(k,it)
+ athetk=athet(k,it,ichir1,ichir2)
+ bthetk=bthet(k,it,ichir1,ichir2)
+ if (it.eq.10) then
+ athetk=athet(k,itype1,ichir11,ichir12)
+ bthetk=bthet(k,itype2,ichir21,ichir22)
+ endif
thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k)
enddo
c write (iout,*) "thet_pred_mean",thet_pred_mean
thet_pred_mean=thet_pred_mean*ss+a0thet(it)
c write (iout,*) "thet_pred_mean",thet_pred_mean
C Derivatives of the "mean" values in gamma1 and gamma2.
- dthetg1=(-athet(1,it)*y(2)+athet(2,it)*y(1))*ss
- dthetg2=(-bthet(1,it)*z(2)+bthet(2,it)*z(1))*ss
+ dthetg1=(-athet(1,it,ichir1,ichir2)*y(2)
+ &+athet(2,it,ichir1,ichir2)*y(1))*ss
+ dthetg2=(-bthet(1,it,ichir1,ichir2)*z(2)
+ & +bthet(2,it,ichir1,ichir2)*z(1))*ss
+ if (it.eq.10) then
+ dthetg1=(-athet(1,itype1,ichir11,ichir12)*y(2)
+ &+athet(2,itype1,ichir11,ichir12)*y(1))*ss
+ dthetg2=(-bthet(1,itype2,ichir21,ichir22)*z(2)
+ & +bthet(2,itype2,ichir21,ichir22)*z(1))*ss
+ endif
if (theta(i).gt.pi-delta) then
call theteng(pi-delta,thet_pred_mean,theta0(it),f0,fprim0,
& E_tc0)
if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)
- 1215 continue
+c 1215 continue
enddo
C Ufff.... We've done all this!!!
+C now constrains
+ ethetacnstr=0.0d0
+C print *,ithetaconstr_start,ithetaconstr_end,"TU"
+ do i=1,ntheta_constr
+ itheta=itheta_constr(i)
+ thetiii=theta(itheta)
+ difi=pinorm(thetiii-theta_constr0(i))
+ if (difi.gt.theta_drange(i)) then
+ difi=difi-theta_drange(i)
+ ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+ gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
+ & +for_thet_constr(i)*difi**3
+ else if (difi.lt.-drange(i)) then
+ difi=difi+drange(i)
+ ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+ gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
+ & +for_thet_constr(i)*difi**3
+ else
+ difi=0.0
+ endif
+C if (energy_dec) then
+C write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",
+C & i,itheta,rad2deg*thetiii,
+C & rad2deg*theta_constr0(i), rad2deg*theta_drange(i),
+C & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,
+C & gloc(itheta+nphi-2,icg)
+C endif
+ enddo
return
end
C---------------------------------------------------------------------------
end
#else
C--------------------------------------------------------------------------
- subroutine ebend(etheta)
+ subroutine ebend(etheta,ethetacnstr)
C
C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
C angles gamma and its derivatives in consecutive thetas and gammas.
include 'COMMON.NAMES'
include 'COMMON.FFIELD'
include 'COMMON.CONTROL'
+ include 'COMMON.TORCNSTR'
double precision coskt(mmaxtheterm),sinkt(mmaxtheterm),
& cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle),
& sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble),
etheta=0.0D0
c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
do i=ithet_start,ithet_end
- if (itype(i-1).eq.21) cycle
+ if (i.eq.2) cycle
+c print *,i,itype(i-1),itype(i),itype(i-2)
+ if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1)
+ & .or.(itype(i).eq.ntyp1)) cycle
+C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF
+
+ if (iabs(itype(i+1)).eq.20) iblock=2
+ if (iabs(itype(i+1)).ne.20) iblock=1
dethetai=0.0d0
dephii=0.0d0
dephii1=0.0d0
theti2=0.5d0*theta(i)
- ityp2=ithetyp(itype(i-1))
+ ityp2=ithetyp((itype(i-1)))
do k=1,nntheterm
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
- if (i.gt.3 .and. itype(i-2).ne.21) then
+ if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
#else
phii=phi(i)
#endif
- ityp1=ithetyp(itype(i-2))
+ ityp1=ithetyp((itype(i-2)))
do k=1,nsingle
cosph1(k)=dcos(k*phii)
sinph1(k)=dsin(k*phii)
sinph1(k)=0.0d0
enddo
endif
- if (i.lt.nres .and. itype(i).ne.21) then
+ if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
if (phii1.ne.phii1) phii1=150.0
#else
phii1=phi(i+1)
#endif
- ityp3=ithetyp(itype(i))
+ ityp3=ithetyp((itype(i)))
do k=1,nsingle
cosph2(k)=dcos(k*phii1)
sinph2(k)=dsin(k*phii1)
c write (iout,*) "i",i," ityp1",itype(i-2),ityp1,
c & " ityp2",itype(i-1),ityp2," ityp3",itype(i),ityp3
c call flush(iout)
- ethetai=aa0thet(ityp1,ityp2,ityp3)
+ ethetai=aa0thet(ityp1,ityp2,ityp3,iblock)
do k=1,ndouble
do l=1,k-1
ccl=cosph1(l)*cosph2(k-l)
enddo
endif
do k=1,ntheterm
- ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3)*sinkt(k)
- dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3)
+ ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3,iblock)*sinkt(k)
+ dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3,iblock)
& *coskt(k)
if (lprn)
- & write (iout,*) "k",k," aathet",aathet(k,ityp1,ityp2,ityp3),
+ & write (iout,*) "k",k," aathet",
+ & aathet(k,ityp1,ityp2,ityp3,iblock),
& " ethetai",ethetai
enddo
if (lprn) then
endif
do m=1,ntheterm2
do k=1,nsingle
- aux=bbthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)
- & +ccthet(k,m,ityp1,ityp2,ityp3)*sinph1(k)
- & +ddthet(k,m,ityp1,ityp2,ityp3)*cosph2(k)
- & +eethet(k,m,ityp1,ityp2,ityp3)*sinph2(k)
+ aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)
+ & +ccthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k)
+ & +ddthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)
+ & +eethet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k)
ethetai=ethetai+sinkt(m)*aux
dethetai=dethetai+0.5d0*m*aux*coskt(m)
dephii=dephii+k*sinkt(m)*(
- & ccthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)-
- & bbthet(k,m,ityp1,ityp2,ityp3)*sinph1(k))
+ & ccthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)-
+ & bbthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k))
dephii1=dephii1+k*sinkt(m)*(
- & eethet(k,m,ityp1,ityp2,ityp3)*cosph2(k)-
- & ddthet(k,m,ityp1,ityp2,ityp3)*sinph2(k))
+ & eethet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)-
+ & ddthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k))
if (lprn)
& write (iout,*) "m",m," k",k," bbthet",
- & bbthet(k,m,ityp1,ityp2,ityp3)," ccthet",
- & ccthet(k,m,ityp1,ityp2,ityp3)," ddthet",
- & ddthet(k,m,ityp1,ityp2,ityp3)," eethet",
- & eethet(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+ & bbthet(k,m,ityp1,ityp2,ityp3,iblock)," ccthet",
+ & ccthet(k,m,ityp1,ityp2,ityp3,iblock)," ddthet",
+ & ddthet(k,m,ityp1,ityp2,ityp3,iblock)," eethet",
+ & eethet(k,m,ityp1,ityp2,ityp3,iblock)," ethetai",ethetai
enddo
enddo
if (lprn)
do m=1,ntheterm3
do k=2,ndouble
do l=1,k-1
- aux=ffthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
- & ffthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+
- & ggthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
- & ggthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)
+ aux=ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+
+ & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l)+
+ & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+
+ & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)
ethetai=ethetai+sinkt(m)*aux
dethetai=dethetai+0.5d0*m*coskt(m)*aux
dephii=dephii+l*sinkt(m)*(
- & -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)-
- & ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
- & ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
- & ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+ & -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)-
+ & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+
+ & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+
+ & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l))
dephii1=dephii1+(k-l)*sinkt(m)*(
- & -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
- & ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
- & ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)-
- & ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+ & -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+
+ & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+
+ & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)-
+ & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l))
if (lprn) then
write (iout,*) "m",m," k",k," l",l," ffthet",
- & ffthet(l,k,m,ityp1,ityp2,ityp3),
- & ffthet(k,l,m,ityp1,ityp2,ityp3)," ggthet",
- & ggthet(l,k,m,ityp1,ityp2,ityp3),
- & ggthet(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+ & ffthet(l,k,m,ityp1,ityp2,ityp3,iblock),
+ & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)," ggthet",
+ & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock),
+ & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock),
+ & " ethetai",ethetai
write (iout,*) cosph1ph2(l,k)*sinkt(m),
& cosph1ph2(k,l)*sinkt(m),
& sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m)
etheta=etheta+ethetai
if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
- gloc(nphi+i-2,icg)=wang*dethetai
+c gloc(nphi+i-2,icg)=wang*dethetai
+ gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai
+ enddo
+C now constrains
+ ethetacnstr=0.0d0
+C print *,ithetaconstr_start,ithetaconstr_end,"TU"
+ do i=1,ntheta_constr
+ itheta=itheta_constr(i)
+ thetiii=theta(itheta)
+ difi=pinorm(thetiii-theta_constr0(i))
+ if (difi.gt.theta_drange(i)) then
+ difi=difi-theta_drange(i)
+ ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+ gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
+ & +for_thet_constr(i)*difi**3
+ else if (difi.lt.-drange(i)) then
+ difi=difi+drange(i)
+ ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+ gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg)
+ & +for_thet_constr(i)*difi**3
+ else
+ difi=0.0
+ endif
+C if (energy_dec) then
+C write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",
+C & i,itheta,rad2deg*thetiii,
+C & rad2deg*theta_constr0(i), rad2deg*theta_drange(i),
+C & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,
+C & gloc(itheta+nphi-2,icg)
+C endif
enddo
return
end
c write (iout,'(a)') 'ESC'
do i=loc_start,loc_end
it=itype(i)
- if (it.eq.21) cycle
+ if (it.eq.ntyp1) cycle
if (it.eq.10) goto 1
- nlobit=nlob(it)
+ nlobit=nlob(iabs(it))
c print *,'i=',i,' it=',it,' nlobit=',nlobit
c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
theti=theta(i+1)-pipol
do iii=-1,1
do j=1,nlobit
- expfac=dexp(bsc(j,it)-0.5D0*contr(j,iii)+emin)
+ expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin)
cd print *,'j=',j,' expfac=',expfac
escloc_i=escloc_i+expfac
do k=1,3
dersc12=0.0d0
do j=1,nlobit
- expfac=dexp(bsc(j,it)-0.5D0*contr(j)+emin)
+ expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j)+emin)
escloc_i=escloc_i+expfac
do k=1,2
dersc(k)=dersc(k)+Ax(k,j)*expfac
delta=0.02d0*pi
escloc=0.0D0
do i=loc_start,loc_end
- if (itype(i).eq.21) cycle
+ if (itype(i).eq.ntyp1) cycle
costtab(i+1) =dcos(theta(i+1))
sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
cosfac=dsqrt(cosfac2)
sinfac2=0.5d0/(1.0d0-costtab(i+1))
sinfac=dsqrt(sinfac2)
- it=itype(i)
+ it=iabs(itype(i))
if (it.eq.10) goto 1
c
C Compute the axes of tghe local cartesian coordinates system; store in
y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
enddo
do j = 1,3
- z_prime(j) = -uz(j,i-1)
+ z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i)))
enddo
c write (2,*) "i",i
c write (2,*) "x_prime",(x_prime(j),j=1,3)
C Compute the energy of the ith side cbain
C
c write (2,*) "xx",xx," yy",yy," zz",zz
- it=itype(i)
+ it=iabs(itype(i))
do j = 1,65
x(j) = sc_parmin(j,it)
enddo
Cc diagnostics - remove later
xx1 = dcos(alph(2))
yy1 = dsin(alph(2))*dcos(omeg(2))
- zz1 = -dsin(alph(2))*dsin(omeg(2))
+c zz1 = -dsin(alph(2))*dsin(omeg(2))
+ zz1 = -dsign(1.0d0,itype(i))*dsin(alph(2))*dsin(omeg(2))
write(2,'(3f8.1,3f9.3,1x,3f9.3)')
& alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,
& xx1,yy1,zz1
dZZ_Ci1(k)=0.0d0
dZZ_Ci(k)=0.0d0
do j=1,3
- dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres)
- dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres)
+ dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)
+ & *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
+ dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)
+ & *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
enddo
dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
c lprn=.true.
etors=0.0D0
do i=iphi_start,iphi_end
- if (itype(i-2).eq.21 .or. itype(i-1).eq.21
- & .or. itype(i).eq.21) cycle
+ if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
+ & .or. itype(i).eq.ntyp1) cycle
itori=itortyp(itype(i-2))
itori1=itortyp(itype(i-1))
phii=phi(i)
difi=phii-phi0(i)
if (difi.gt.drange(i)) then
difi=difi-drange(i)
- edihcnstr=edihcnstr+0.25d0*ftors*difi**4
- gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
else if (difi.lt.-drange(i)) then
difi=difi+drange(i)
- edihcnstr=edihcnstr+0.25d0*ftors*difi**4
- gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
endif
! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
c lprn=.true.
etors=0.0D0
do i=iphi_start,iphi_end
- if (itype(i-2).eq.21 .or. itype(i-1).eq.21
- & .or. itype(i).eq.21) cycle
+ if (i.le.2) cycle
+ if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
+ & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
if (itel(i-2).eq.0 .or. itel(i-1).eq.0) goto 1215
+ if (iabs(itype(i)).eq.20) then
+ iblock=2
+ else
+ iblock=1
+ endif
itori=itortyp(itype(i-2))
itori1=itortyp(itype(i-1))
phii=phi(i)
gloci=0.0D0
C Regular cosine and sine terms
- do j=1,nterm(itori,itori1)
- v1ij=v1(j,itori,itori1)
- v2ij=v2(j,itori,itori1)
+ do j=1,nterm(itori,itori1,iblock)
+ v1ij=v1(j,itori,itori1,iblock)
+ v2ij=v2(j,itori,itori1,iblock)
cosphi=dcos(j*phii)
sinphi=dsin(j*phii)
etors=etors+v1ij*cosphi+v2ij*sinphi
C
cosphi=dcos(0.5d0*phii)
sinphi=dsin(0.5d0*phii)
- do j=1,nlor(itori,itori1)
+ do j=1,nlor(itori,itori1,iblock)
vl1ij=vlor1(j,itori,itori1)
vl2ij=vlor2(j,itori,itori1)
vl3ij=vlor3(j,itori,itori1)
gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom
enddo
C Subtract the constant term
- etors=etors-v0(itori,itori1)
+ etors=etors-v0(itori,itori1,iblock)
if (lprn)
& write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
& restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
- & (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6)
+ & (v1(j,itori,itori1,1),j=1,6),(v2(j,itori,itori1,1),j=1,6)
gloc(i-3,icg)=gloc(i-3,icg)+wtor*fact*gloci
c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
1215 continue
edihi=0.0d0
if (difi.gt.drange(i)) then
difi=difi-drange(i)
- edihcnstr=edihcnstr+0.25d0*ftors*difi**4
- gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
- edihi=0.25d0*ftors*difi**4
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
+ edihi=0.25d0*ftors(i)*difi**4
else if (difi.lt.-drange(i)) then
difi=difi+drange(i)
- edihcnstr=edihcnstr+0.25d0*ftors*difi**4
- gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
- edihi=0.25d0*ftors*difi**4
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
+ edihi=0.25d0*ftors(i)*difi**4
else
difi=0.0d0
endif
c lprn=.true.
etors_d=0.0D0
do i=iphi_start,iphi_end-1
- if (itype(i-2).eq.21 .or. itype(i-1).eq.21
- & .or. itype(i).eq.21 .or. itype(i+1).eq.21) cycle
+ if (i.le.3) cycle
+ if ((itype(i-2).eq.ntyp1).or.itype(i-3).eq.ntyp1.or.
+ & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1).or.
+ & (itype(i+1).eq.ntyp1)) cycle
if (itel(i-2).eq.0 .or. itel(i-1).eq.0 .or. itel(i).eq.0)
& goto 1215
itori=itortyp(itype(i-2))
phii1=phi(i+1)
gloci1=0.0D0
gloci2=0.0D0
+ iblock=1
+ if (iabs(itype(i+1)).eq.20) iblock=2
C Regular cosine and sine terms
- do j=1,ntermd_1(itori,itori1,itori2)
- v1cij=v1c(1,j,itori,itori1,itori2)
- v1sij=v1s(1,j,itori,itori1,itori2)
- v2cij=v1c(2,j,itori,itori1,itori2)
- v2sij=v1s(2,j,itori,itori1,itori2)
+ do j=1,ntermd_1(itori,itori1,itori2,iblock)
+ v1cij=v1c(1,j,itori,itori1,itori2,iblock)
+ v1sij=v1s(1,j,itori,itori1,itori2,iblock)
+ v2cij=v1c(2,j,itori,itori1,itori2,iblock)
+ v2sij=v1s(2,j,itori,itori1,itori2,iblock)
cosphi1=dcos(j*phii)
sinphi1=dsin(j*phii)
cosphi2=dcos(j*phii1)
gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
enddo
- do k=2,ntermd_2(itori,itori1,itori2)
+ do k=2,ntermd_2(itori,itori1,itori2,iblock)
do l=1,k-1
- v1cdij = v2c(k,l,itori,itori1,itori2)
- v2cdij = v2c(l,k,itori,itori1,itori2)
- v1sdij = v2s(k,l,itori,itori1,itori2)
- v2sdij = v2s(l,k,itori,itori1,itori2)
+ v1cdij = v2c(k,l,itori,itori1,itori2,iblock)
+ v2cdij = v2c(l,k,itori,itori1,itori2,iblock)
+ v1sdij = v2s(k,l,itori,itori1,itori2,iblock)
+ v2sdij = v2s(l,k,itori,itori1,itori2,iblock)
cosphi1p2=dcos(l*phii+(k-l)*phii1)
cosphi1m2=dcos(l*phii-(k-l)*phii1)
sinphi1p2=dsin(l*phii+(k-l)*phii1)
c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
esccor=0.0D0
do i=itau_start,itau_end
- if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1) cycle
+ if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
esccor_ii=0.0D0
isccori=isccortyp(itype(i-2))
isccori1=isccortyp(itype(i-1))
integer dimen1,dimen2,atom,indx
double precision buffer(dimen1,dimen2)
double precision zapas
- common /contacts_hb/ zapas(3,20,maxres,7),
- & facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres),
- & num_cont_hb(maxres),jcont_hb(20,maxres)
+ common /contacts_hb/ zapas(3,ntyp,maxres,7),
+ & facont_hb(ntyp,maxres),ees0p(ntyp,maxres),ees0m(ntyp,maxres),
+ & num_cont_hb(maxres),jcont_hb(ntyp,maxres)
num_kont=buffer(1,indx+26)
num_kont_old=num_cont_hb(atom)
num_cont_hb(atom)=num_kont+num_kont_old
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
include 'COMMON.CONTACTS'
+ include 'COMMON.SHIELD'
+
double precision gx(3),gx1(3)
logical lprn
lprn=.false.
& ekont*(coeffp*ees0pij*gacontp_hb3(ll,kk,k)+
& coeffm*ees0mij*gacontm_hb3(ll,kk,k))
enddo
- enddo
+ enddo
+ if (shield_mode.gt.0) then
+ j=ees0plist(jj,i)
+ l=ees0plist(kk,k)
+C print *,i,j,fac_shield(i),fac_shield(j),
+C &fac_shield(k),fac_shield(l)
+ if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+ & (fac_shield(k).gt.0).and.(fac_shield(l).gt.0)) then
+ do ilist=1,ishield_list(i)
+ iresshield=shield_list(ilist,i)
+ do m=1,3
+ rlocshield=grad_shield_side(m,ilist,i)*ehbcorr/fac_shield(i)
+C & *2.0
+ gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(m,ilist,i)*ehbcorr/fac_shield(i)
+ gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
+ &+rlocshield
+ enddo
+ enddo
+ do ilist=1,ishield_list(j)
+ iresshield=shield_list(ilist,j)
+ do m=1,3
+ rlocshield=grad_shield_side(m,ilist,j)*ehbcorr/fac_shield(j)
+C & *2.0
+ gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(m,ilist,j)*ehbcorr/fac_shield(j)
+ gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
+ & +rlocshield
+ enddo
+ enddo
+ do ilist=1,ishield_list(k)
+ iresshield=shield_list(ilist,k)
+ do m=1,3
+ rlocshield=grad_shield_side(m,ilist,k)*ehbcorr/fac_shield(k)
+C & *2.0
+ gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(m,ilist,k)*ehbcorr/fac_shield(k)
+ gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
+ & +rlocshield
+ enddo
+ enddo
+ do ilist=1,ishield_list(l)
+ iresshield=shield_list(ilist,l)
+ do m=1,3
+ rlocshield=grad_shield_side(m,ilist,l)*ehbcorr/fac_shield(l)
+C & *2.0
+ gshieldx_ec(m,iresshield)=gshieldx_ec(m,iresshield)+
+ & rlocshield
+ & +grad_shield_loc(m,ilist,l)*ehbcorr/fac_shield(l)
+ gshieldc_ec(m,iresshield-1)=gshieldc_ec(m,iresshield-1)
+ & +rlocshield
+ enddo
+ enddo
+C print *,gshieldx(m,iresshield)
+ do m=1,3
+ gshieldc_ec(m,i)=gshieldc_ec(m,i)+
+ & grad_shield(m,i)*ehbcorr/fac_shield(i)
+ gshieldc_ec(m,j)=gshieldc_ec(m,j)+
+ & grad_shield(m,j)*ehbcorr/fac_shield(j)
+ gshieldc_ec(m,i-1)=gshieldc_ec(m,i-1)+
+ & grad_shield(m,i)*ehbcorr/fac_shield(i)
+ gshieldc_ec(m,j-1)=gshieldc_ec(m,j-1)+
+ & grad_shield(m,j)*ehbcorr/fac_shield(j)
+
+ gshieldc_ec(m,k)=gshieldc_ec(m,k)+
+ & grad_shield(m,k)*ehbcorr/fac_shield(k)
+ gshieldc_ec(m,l)=gshieldc_ec(m,l)+
+ & grad_shield(m,l)*ehbcorr/fac_shield(l)
+ gshieldc_ec(m,k-1)=gshieldc_ec(m,k-1)+
+ & grad_shield(m,k)*ehbcorr/fac_shield(k)
+ gshieldc_ec(m,l-1)=gshieldc_ec(m,l-1)+
+ & grad_shield(m,l)*ehbcorr/fac_shield(l)
+
+ enddo
+ endif
+ endif
endif
ehbcorr=ekont*ees
return
& auxmat(2,2)
iti1 = itortyp(itype(i+1))
if (j.lt.nres-1) then
- itj1 = itortyp(itype(j+1))
+ if (itype(j).le.ntyp) then
+ itj1 = itortyp(itype(j+1))
+ else
+ itj1=ntortyp+1
+ endif
else
itj1=ntortyp+1
endif
enddo
if (l.eq.j+1) then
C parallel orientation of the two CA-CA-CA frames.
- if (i.gt.1) then
+c if (i.gt.1) then
+ if (i.gt.1 .and. itype(i).le.ntyp) then
iti=itortyp(itype(i))
else
iti=ntortyp+1
endif
itk1=itortyp(itype(k+1))
itj=itortyp(itype(j))
- if (l.lt.nres-1) then
+c if (l.lt.nres-1) then
+ if (l.lt.nres-1 .and. itype(l+1).le.ntyp) then
itl1=itortyp(itype(l+1))
else
itl1=ntortyp+1
C End vectors
else
C Antiparallel orientation of the two CA-CA-CA frames.
- if (i.gt.1) then
+c if (i.gt.1) then
+ if (i.gt.1 .and. itype(i).le.ntyp) then
iti=itortyp(itype(i))
else
iti=ntortyp+1
itk1=itortyp(itype(k+1))
itl=itortyp(itype(l))
itj=itortyp(itype(j))
- if (j.lt.nres-1) then
+c if (j.lt.nres-1) then
+ if (j.lt.nres-1 .and. itype(j+1).le.ntyp) then
itj1=itortyp(itype(j+1))
else
itj1=ntortyp+1
C 4/7/01 AL Component s1 was removed, because it pertains to the respective
C energy moment and not to the cluster cumulant.
iti=itortyp(itype(i))
- if (j.lt.nres-1) then
+c if (j.lt.nres-1) then
+ if (j.lt.nres-1 .and. itype(j+1).le.ntyp) then
itj1=itortyp(itype(j+1))
else
itj1=ntortyp+1
endif
itk=itortyp(itype(k))
itk1=itortyp(itype(k+1))
- if (l.lt.nres-1) then
+c if (l.lt.nres-1) then
+ if (l.lt.nres-1 .and. itype(l+1).le.ntyp) then
itl1=itortyp(itype(l+1))
else
itl1=ntortyp+1
cd write (2,*) 'eello_graph4: wturn6',wturn6
iti=itortyp(itype(i))
itj=itortyp(itype(j))
- if (j.lt.nres-1) then
+c if (j.lt.nres-1) then
+ if (j.lt.nres-1 .and. itype(j+1).le.ntyp) then
itj1=itortyp(itype(j+1))
else
itj1=ntortyp+1
endif
itk=itortyp(itype(k))
- if (k.lt.nres-1) then
+c if (k.lt.nres-1) then
+ if (k.lt.nres-1 .and. itype(k+1).le.ntyp) then
itk1=itortyp(itype(k+1))
else
itk1=ntortyp+1
scalar=sc
return
end
+C-----------------------------------------------------------------------
+ double precision function sscale(r)
+ double precision r,gamm
+ include "COMMON.SPLITELE"
+ if(r.lt.r_cut-rlamb) then
+ sscale=1.0d0
+ else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+ gamm=(r-(r_cut-rlamb))/rlamb
+ sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0)
+ else
+ sscale=0d0
+ endif
+ return
+ end
+C-----------------------------------------------------------------------
+C-----------------------------------------------------------------------
+ double precision function sscagrad(r)
+ double precision r,gamm
+ include "COMMON.SPLITELE"
+ if(r.lt.r_cut-rlamb) then
+ sscagrad=0.0d0
+ else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+ gamm=(r-(r_cut-rlamb))/rlamb
+ sscagrad=gamm*(6*gamm-6.0d0)/rlamb
+ else
+ sscagrad=0.0d0
+ endif
+ return
+ end
+C-----------------------------------------------------------------------
+C first for shielding is setting of function of side-chains
+ subroutine set_shield_fac2
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SHIELD'
+ include 'COMMON.INTERACT'
+C this is the squar root 77 devided by 81 the epislion in lipid (in protein)
+ double precision div77_81/0.974996043d0/,
+ &div4_81/0.2222222222d0/,sh_frac_dist_grad(3)
+
+C the vector between center of side_chain and peptide group
+ double precision pep_side(3),long,side_calf(3),
+ &pept_group(3),costhet_grad(3),cosphi_grad_long(3),
+ &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3)
+C the line belowe needs to be changed for FGPROC>1
+ do i=1,nres-1
+ if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+ ishield_list(i)=0
+Cif there two consequtive dummy atoms there is no peptide group between them
+C the line below has to be changed for FGPROC>1
+ VolumeTotal=0.0
+ do k=1,nres
+ if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle
+ dist_pep_side=0.0
+ dist_side_calf=0.0
+ do j=1,3
+C first lets set vector conecting the ithe side-chain with kth side-chain
+ pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0
+C pep_side(j)=2.0d0
+C and vector conecting the side-chain with its proper calfa
+ side_calf(j)=c(j,k+nres)-c(j,k)
+C side_calf(j)=2.0d0
+ pept_group(j)=c(j,i)-c(j,i+1)
+C lets have their lenght
+ dist_pep_side=pep_side(j)**2+dist_pep_side
+ dist_side_calf=dist_side_calf+side_calf(j)**2
+ dist_pept_group=dist_pept_group+pept_group(j)**2
+ enddo
+ dist_pep_side=dsqrt(dist_pep_side)
+ dist_pept_group=dsqrt(dist_pept_group)
+ dist_side_calf=dsqrt(dist_side_calf)
+ do j=1,3
+ pep_side_norm(j)=pep_side(j)/dist_pep_side
+ side_calf_norm(j)=dist_side_calf
+ enddo
+C now sscale fraction
+ sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield
+C print *,buff_shield,"buff"
+C now sscale
+ if (sh_frac_dist.le.0.0) cycle
+C If we reach here it means that this side chain reaches the shielding sphere
+C Lets add him to the list for gradient
+ ishield_list(i)=ishield_list(i)+1
+C ishield_list is a list of non 0 side-chain that contribute to factor gradient
+C this list is essential otherwise problem would be O3
+ shield_list(ishield_list(i),i)=k
+C Lets have the sscale value
+ if (sh_frac_dist.gt.1.0) then
+ scale_fac_dist=1.0d0
+ do j=1,3
+ sh_frac_dist_grad(j)=0.0d0
+ enddo
+ else
+ scale_fac_dist=-sh_frac_dist*sh_frac_dist
+ & *(2.0d0*sh_frac_dist-3.0d0)
+ fac_help_scale=6.0d0*(sh_frac_dist-sh_frac_dist**2)
+ & /dist_pep_side/buff_shield*0.5d0
+C remember for the final gradient multiply sh_frac_dist_grad(j)
+C for side_chain by factor -2 !
+ do j=1,3
+ sh_frac_dist_grad(j)=fac_help_scale*pep_side(j)
+C sh_frac_dist_grad(j)=0.0d0
+C scale_fac_dist=1.0d0
+C print *,"jestem",scale_fac_dist,fac_help_scale,
+C & sh_frac_dist_grad(j)
+ enddo
+ endif
+C this is what is now we have the distance scaling now volume...
+ short=short_r_sidechain(itype(k))
+ long=long_r_sidechain(itype(k))
+ costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2)
+ sinthet=short/dist_pep_side*costhet
+C now costhet_grad
+C costhet=0.6d0
+C sinthet=0.8
+ costhet_fac=costhet**3*short**2*(-0.5d0)/dist_pep_side**4
+C sinthet_fac=costhet**2*0.5d0*(short**3/dist_pep_side**4*costhet
+C & -short/dist_pep_side**2/costhet)
+C costhet_fac=0.0d0
+ do j=1,3
+ costhet_grad(j)=costhet_fac*pep_side(j)
+ enddo
+C remember for the final gradient multiply costhet_grad(j)
+C for side_chain by factor -2 !
+C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1
+C pep_side0pept_group is vector multiplication
+ pep_side0pept_group=0.0d0
+ do j=1,3
+ pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j)
+ enddo
+ cosalfa=(pep_side0pept_group/
+ & (dist_pep_side*dist_side_calf))
+ fac_alfa_sin=1.0d0-cosalfa**2
+ fac_alfa_sin=dsqrt(fac_alfa_sin)
+ rkprim=fac_alfa_sin*(long-short)+short
+C rkprim=short
+
+C now costhet_grad
+ cosphi=1.0d0/dsqrt(1.0d0+rkprim**2/dist_pep_side**2)
+C cosphi=0.6
+ cosphi_fac=cosphi**3*rkprim**2*(-0.5d0)/dist_pep_side**4
+ sinphi=rkprim/dist_pep_side/dsqrt(1.0d0+rkprim**2/
+ & dist_pep_side**2)
+C sinphi=0.8
+ do j=1,3
+ cosphi_grad_long(j)=cosphi_fac*pep_side(j)
+ &+cosphi**3*0.5d0/dist_pep_side**2*(-rkprim)
+ &*(long-short)/fac_alfa_sin*cosalfa/
+ &((dist_pep_side*dist_side_calf))*
+ &((side_calf(j))-cosalfa*
+ &((pep_side(j)/dist_pep_side)*dist_side_calf))
+C cosphi_grad_long(j)=0.0d0
+ cosphi_grad_loc(j)=cosphi**3*0.5d0/dist_pep_side**2*(-rkprim)
+ &*(long-short)/fac_alfa_sin*cosalfa
+ &/((dist_pep_side*dist_side_calf))*
+ &(pep_side(j)-
+ &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side)
+C cosphi_grad_loc(j)=0.0d0
+ enddo
+C print *,sinphi,sinthet
+ VofOverlap=VSolvSphere/2.0d0*(1.0d0-dsqrt(1.0d0-sinphi*sinthet))
+ & /VSolvSphere_div
+C & *wshield
+C now the gradient...
+ do j=1,3
+ grad_shield(j,i)=grad_shield(j,i)
+C gradient po skalowaniu
+ & +(sh_frac_dist_grad(j)*VofOverlap
+C gradient po costhet
+ & +scale_fac_dist*VSolvSphere/VSolvSphere_div/4.0d0*
+ &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*(
+ & sinphi/sinthet*costhet*costhet_grad(j)
+ & +sinthet/sinphi*cosphi*cosphi_grad_long(j)))
+ & )*wshield
+C grad_shield_side is Cbeta sidechain gradient
+ grad_shield_side(j,ishield_list(i),i)=
+ & (sh_frac_dist_grad(j)*-2.0d0
+ & *VofOverlap
+ & -scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*
+ &(1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*(
+ & sinphi/sinthet*costhet*costhet_grad(j)
+ & +sinthet/sinphi*cosphi*cosphi_grad_long(j)))
+ & )*wshield
+
+ grad_shield_loc(j,ishield_list(i),i)=
+ & scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*
+ &(1.0d0/(dsqrt(1.0d0-sinphi*sinthet))*(
+ & sinthet/sinphi*cosphi*cosphi_grad_loc(j)
+ & ))
+ & *wshield
+ enddo
+ VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
+ enddo
+ fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield)
+C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
+ enddo
+ return
+ end
+C first for shielding is setting of function of side-chains
+ subroutine set_shield_fac
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SHIELD'
+ include 'COMMON.INTERACT'
+C this is the squar root 77 devided by 81 the epislion in lipid (in protein)
+ double precision div77_81/0.974996043d0/,
+ &div4_81/0.2222222222d0/,sh_frac_dist_grad(3)
+
+C the vector between center of side_chain and peptide group
+ double precision pep_side(3),long,side_calf(3),
+ &pept_group(3),costhet_grad(3),cosphi_grad_long(3),
+ &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3)
+C the line belowe needs to be changed for FGPROC>1
+ do i=1,nres-1
+ if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+ ishield_list(i)=0
+Cif there two consequtive dummy atoms there is no peptide group between them
+C the line below has to be changed for FGPROC>1
+ VolumeTotal=0.0
+ do k=1,nres
+ if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle
+ dist_pep_side=0.0
+ dist_side_calf=0.0
+ do j=1,3
+C first lets set vector conecting the ithe side-chain with kth side-chain
+ pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0
+C pep_side(j)=2.0d0
+C and vector conecting the side-chain with its proper calfa
+ side_calf(j)=c(j,k+nres)-c(j,k)
+C side_calf(j)=2.0d0
+ pept_group(j)=c(j,i)-c(j,i+1)
+C lets have their lenght
+ dist_pep_side=pep_side(j)**2+dist_pep_side
+ dist_side_calf=dist_side_calf+side_calf(j)**2
+ dist_pept_group=dist_pept_group+pept_group(j)**2
+ enddo
+ dist_pep_side=dsqrt(dist_pep_side)
+ dist_pept_group=dsqrt(dist_pept_group)
+ dist_side_calf=dsqrt(dist_side_calf)
+ do j=1,3
+ pep_side_norm(j)=pep_side(j)/dist_pep_side
+ side_calf_norm(j)=dist_side_calf
+ enddo
+C now sscale fraction
+ sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield
+C print *,buff_shield,"buff"
+C now sscale
+ if (sh_frac_dist.le.0.0) cycle
+C If we reach here it means that this side chain reaches the shielding sphere
+C Lets add him to the list for gradient
+ ishield_list(i)=ishield_list(i)+1
+C ishield_list is a list of non 0 side-chain that contribute to factor gradient
+C this list is essential otherwise problem would be O3
+ shield_list(ishield_list(i),i)=k
+C Lets have the sscale value
+ if (sh_frac_dist.gt.1.0) then
+ scale_fac_dist=1.0d0
+ do j=1,3
+ sh_frac_dist_grad(j)=0.0d0
+ enddo
+ else
+ scale_fac_dist=-sh_frac_dist*sh_frac_dist
+ & *(2.0*sh_frac_dist-3.0d0)
+ fac_help_scale=6.0*(sh_frac_dist-sh_frac_dist**2)
+ & /dist_pep_side/buff_shield*0.5
+C remember for the final gradient multiply sh_frac_dist_grad(j)
+C for side_chain by factor -2 !
+ do j=1,3
+ sh_frac_dist_grad(j)=fac_help_scale*pep_side(j)
+C print *,"jestem",scale_fac_dist,fac_help_scale,
+C & sh_frac_dist_grad(j)
+ enddo
+ endif
+C if ((i.eq.3).and.(k.eq.2)) then
+C print *,i,sh_frac_dist,dist_pep,fac_help_scale,scale_fac_dist
+C & ,"TU"
+C endif
+
+C this is what is now we have the distance scaling now volume...
+ short=short_r_sidechain(itype(k))
+ long=long_r_sidechain(itype(k))
+ costhet=1.0d0/dsqrt(1.0+short**2/dist_pep_side**2)
+C now costhet_grad
+C costhet=0.0d0
+ costhet_fac=costhet**3*short**2*(-0.5)/dist_pep_side**4
+C costhet_fac=0.0d0
+ do j=1,3
+ costhet_grad(j)=costhet_fac*pep_side(j)
+ enddo
+C remember for the final gradient multiply costhet_grad(j)
+C for side_chain by factor -2 !
+C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1
+C pep_side0pept_group is vector multiplication
+ pep_side0pept_group=0.0
+ do j=1,3
+ pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j)
+ enddo
+ cosalfa=(pep_side0pept_group/
+ & (dist_pep_side*dist_side_calf))
+ fac_alfa_sin=1.0-cosalfa**2
+ fac_alfa_sin=dsqrt(fac_alfa_sin)
+ rkprim=fac_alfa_sin*(long-short)+short
+C now costhet_grad
+ cosphi=1.0d0/dsqrt(1.0+rkprim**2/dist_pep_side**2)
+ cosphi_fac=cosphi**3*rkprim**2*(-0.5)/dist_pep_side**4
+
+ do j=1,3
+ cosphi_grad_long(j)=cosphi_fac*pep_side(j)
+ &+cosphi**3*0.5/dist_pep_side**2*(-rkprim)
+ &*(long-short)/fac_alfa_sin*cosalfa/
+ &((dist_pep_side*dist_side_calf))*
+ &((side_calf(j))-cosalfa*
+ &((pep_side(j)/dist_pep_side)*dist_side_calf))
+
+ cosphi_grad_loc(j)=cosphi**3*0.5/dist_pep_side**2*(-rkprim)
+ &*(long-short)/fac_alfa_sin*cosalfa
+ &/((dist_pep_side*dist_side_calf))*
+ &(pep_side(j)-
+ &cosalfa*side_calf(j)/dist_side_calf*dist_pep_side)
+ enddo
+
+ VofOverlap=VSolvSphere/2.0d0*(1.0-costhet)*(1.0-cosphi)
+ & /VSolvSphere_div
+ & *wshield
+C now the gradient...
+C grad_shield is gradient of Calfa for peptide groups
+C write(iout,*) "shield_compon",i,k,VSolvSphere,scale_fac_dist,
+C & costhet,cosphi
+C write(iout,*) "cosphi_compon",i,k,pep_side0pept_group,
+C & dist_pep_side,dist_side_calf,c(1,k+nres),c(1,k),itype(k)
+ do j=1,3
+ grad_shield(j,i)=grad_shield(j,i)
+C gradient po skalowaniu
+ & +(sh_frac_dist_grad(j)
+C gradient po costhet
+ &-scale_fac_dist*costhet_grad(j)/(1.0-costhet)
+ &-scale_fac_dist*(cosphi_grad_long(j))
+ &/(1.0-cosphi) )*div77_81
+ &*VofOverlap
+C grad_shield_side is Cbeta sidechain gradient
+ grad_shield_side(j,ishield_list(i),i)=
+ & (sh_frac_dist_grad(j)*-2.0d0
+ & +scale_fac_dist*costhet_grad(j)*2.0d0/(1.0-costhet)
+ & +scale_fac_dist*(cosphi_grad_long(j))
+ & *2.0d0/(1.0-cosphi))
+ & *div77_81*VofOverlap
+
+ grad_shield_loc(j,ishield_list(i),i)=
+ & scale_fac_dist*cosphi_grad_loc(j)
+ & *2.0d0/(1.0-cosphi)
+ & *div77_81*VofOverlap
+ enddo
+ VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
+ enddo
+ fac_shield(i)=VolumeTotal*div77_81+div4_81
+C write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
+ enddo
+ return
+ end
+C--------------------------------------------------------------------------
+C-----------------------------------------------------------------------
+ double precision function sscalelip(r)
+ double precision r,gamm
+ include "COMMON.SPLITELE"
+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)
+C else
+C sscale=0d0
+C endif
+ return
+ end
+C-----------------------------------------------------------------------
+ double precision function sscagradlip(r)
+ double precision r,gamm
+ include "COMMON.SPLITELE"
+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)
+C else
+C sscagrad=0.0d0
+C endif
+ return
+ end
+
+C-----------------------------------------------------------------------
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+ subroutine Eliptransfer(eliptran)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CALC'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SPLITELE'
+ include 'COMMON.SBRIDGE'
+C this is done by Adasko
+C print *,"wchodze"
+C structure of box:
+C water
+C--bordliptop-- buffore starts
+C--bufliptop--- here true lipid starts
+C lipid
+C--buflipbot--- lipid ends buffore starts
+C--bordlipbot--buffore ends
+ eliptran=0.0
+ write(iout,*) "I am in?"
+ do i=1,nres
+C do i=1,1
+ if (itype(i).eq.ntyp1) cycle
+ positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize))
+ if (positi.le.0) positi=positi+boxzsize
+C print *,i
+C first for peptide groups
+c for each residue check if it is in lipid or lipid water border area
+ if ((positi.gt.bordlipbot)
+ &.and.(positi.lt.bordliptop)) then
+C the energy transfer exist
+ if (positi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*pepliptran
+ gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0
+ gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
+C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
+ elseif (positi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*pepliptran
+ gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0
+ gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
+C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
+C print *, "doing sscalefor top part"
+C print *,i,sslip,fracinbuf,ssgradlip
+ else
+ eliptran=eliptran+pepliptran
+C print *,"I am in true lipid"
+ endif
+C else
+C eliptran=elpitran+0.0 ! I am in water
+ endif
+ enddo
+C print *, "nic nie bylo w lipidzie?"
+C now multiply all by the peptide group transfer factor
+C eliptran=eliptran*pepliptran
+C now the same for side chains
+CV do i=1,1
+ do i=1,nres
+ if (itype(i).eq.ntyp1) cycle
+ positi=(mod(c(3,i+nres),boxzsize))
+ if (positi.le.0) positi=positi+boxzsize
+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)
+C print *,positi,bordlipbot,buflipbot
+ if ((positi.gt.bordlipbot)
+ & .and.(positi.lt.bordliptop)) then
+C the energy transfer exist
+ if (positi.lt.buflipbot) then
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*liptranene(itype(i))
+ gliptranx(3,i)=gliptranx(3,i)
+ &+ssgradlip*liptranene(itype(i))
+ gliptranc(3,i-1)= gliptranc(3,i-1)
+ &+ssgradlip*liptranene(itype(i))
+C print *,"doing sccale for lower part"
+ elseif (positi.gt.bufliptop) then
+ fracinbuf=1.0d0-
+ &((bordliptop-positi)/lipbufthick)
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*liptranene(itype(i))
+ gliptranx(3,i)=gliptranx(3,i)
+ &+ssgradlip*liptranene(itype(i))
+ gliptranc(3,i-1)= gliptranc(3,i-1)
+ &+ssgradlip*liptranene(itype(i))
+C print *, "doing sscalefor top part",sslip,fracinbuf
+ else
+ eliptran=eliptran+liptranene(itype(i))
+C print *,"I am in true lipid"
+ endif
+ endif ! if in lipid or buffor
+C else
+C eliptran=elpitran+0.0 ! I am in water
+ enddo
+ return
+ end
+C-------------------------------------------------------------------------------------