!-----------------------------------------------------------------------------
! Maximum number of SC local term fitting function coefficiants
integer,parameter :: maxsccoef=65
+! Maximum number of local shielding effectors
+ integer,parameter :: maxcontsshi=50
!-----------------------------------------------------------------------------
! commom.calc common/calc/
!-----------------------------------------------------------------------------
! common /contacts/
! Change 12/1/95 - common block CONTACTS1 included.
! common /contacts1/
+
integer,dimension(:),allocatable :: num_cont !(maxres)
integer,dimension(:,:),allocatable :: jcont !(maxconts,maxres)
- real(kind=8),dimension(:,:),allocatable :: facont !(maxconts,maxres)
+ real(kind=8),dimension(:,:),allocatable :: facont,ees0plist !(maxconts,maxres)
real(kind=8),dimension(:,:,:),allocatable :: gacont !(3,maxconts,maxres)
+ integer,dimension(:),allocatable :: ishield_list
+ integer,dimension(:,:),allocatable :: shield_list
!
! 12/26/95 - H-bonding contacts
! common /contacts_hb/
real(kind=8),dimension(:,:),allocatable :: gvdwc,gelc,gelc_long,&
gvdwpp,gvdwc_scpp,gradx_scp,gvdwc_scp,ghpbx,ghpbc,&
gradcorr,gradcorr_long,gradcorr5_long,gradcorr6_long,&
- gcorr6_turn_long,gradxorr,gradcorr5,gradcorr6 !(3,maxres)
+ gcorr6_turn_long,gradxorr,gradcorr5,gradcorr6,gliptran,gliptranc,&
+ gliptranx, &
+ gshieldx,gshieldc,gshieldc_loc,gshieldx_ec,&
+ gshieldc_ec,gshieldc_loc_ec,gshieldx_t3, &
+ gshieldc_t3,gshieldc_loc_t3,gshieldx_t4,gshieldc_t4, &
+ gshieldc_loc_t4,gshieldx_ll,gshieldc_ll,gshieldc_loc_ll,&
+ grad_shield,gg_tube,gg_tube_sc,gradafm !(3,maxres)
! real(kind=8),dimension(:,:),allocatable :: gloc,gloc_x !(maxvar,2)
real(kind=8),dimension(:,:),allocatable :: gel_loc,gel_loc_long,&
gcorr3_turn,gcorr4_turn,gcorr6_turn,gradb,gradbx !(3,maxres)
! real(kind=8),dimension(:,:,:),allocatable :: dtheta !(3,2,maxres)
real(kind=8),dimension(:,:),allocatable :: gscloc,gsclocx !(3,maxres)
! real(kind=8),dimension(:,:,:),allocatable :: dphi,dalpha,domega !(3,3,maxres)
+ real(kind=8),dimension(:,:,:),allocatable :: grad_shield_side, &
+ grad_shield_loc ! (3,maxcontsshileding,maxnres)
! integer :: nfl,icg
! common /deriv_loc/
+ real(kind=8), dimension(:),allocatable :: fac_shield
real(kind=8),dimension(3,5,2) :: derx,derx_turn
! common /deriv_scloc/
real(kind=8),dimension(:,:),allocatable :: dXX_C1tab,dYY_C1tab,&
integer :: n_corr,n_corr1,ierror
real(kind=8) :: etors,edihcnstr,etors_d,esccor,ehpb
real(kind=8) :: evdw,evdw1,evdw2,evdw2_14,escloc,ees,eel_loc
- real(kind=8) :: eello_turn3,eello_turn4,estr,ebe
+ real(kind=8) :: eello_turn3,eello_turn4,estr,ebe,eliptran,etube, &
+ Eafmforce,ethetacnstr
real(kind=8) :: ecorr,ecorr5,ecorr6,eturn6
#ifdef MPI
#endif
!
! Compute the side-chain and electrostatic interaction energy
-!
+! print *, "Before EVDW"
! goto (101,102,103,104,105,106) ipot
select case(ipot)
! Lennard-Jones potential.
! 50 continue
end select
! continue
-
+! print *,"after EGB"
+! shielding effect
+ if (shield_mode.eq.2) then
+ call set_shield_fac2
+ endif
!mc
!mc Sep-06: egb takes care of dynamic ss bonds too
!mc
#ifdef TIMING
time_vec=time_vec+MPI_Wtime()-time01
#endif
-! print *,"Processor",myrank," left VEC_AND_DERIV"
+! print *,"Processor",myrank," left VEC_AND_DERIV"
if (ipot.lt.6) then
#ifdef SPLITELE
+! print *,"after ipot if", ipot
if (welec.gt.0d0.or.wvdwpp.gt.0d0.or.wel_loc.gt.0d0.or. &
wturn3.gt.0d0.or.wturn4.gt.0d0 .or. wcorr.gt.0.0d0 &
.or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0 &
.or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0 &
.or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then
#endif
+! print *,"just befor eelec call"
call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
-! write (iout,*) "ELEC calc"
+! write (iout,*) "ELEC calc"
else
ees=0.0d0
evdw1=0.0d0
! write (iout,*) "Soft-sphere SCP potential"
call escp_soft_sphere(evdw2,evdw2_14)
endif
-!elwrite(iout,*) "in etotal before ebond",ipot
+! write(iout,*) "in etotal before ebond",ipot
!
! Calculate the bond-stretching energy
!
call ebond(estr)
-!elwrite(iout,*) "in etotal afer ebond",ipot
+! write(iout,*) "in etotal afer ebond",ipot
!
! Calculate the disulfide-bridge and other energy and the contributions
! Calculate the virtual-bond-angle energy.
!
if (wang.gt.0d0) then
- call ebend(ebe)
+ call ebend(ebe,ethetacnstr)
else
ebe=0
endif
Uconst=0.0d0
Uconst_back=0.0d0
endif
-!elwrite(iout,*) "after Econstr"
+ call flush(iout)
+! write(iout,*) "after Econstr"
+
+ if (wliptran.gt.0) then
+! print *,"PRZED WYWOLANIEM"
+ call Eliptransfer(eliptran)
+ else
+ eliptran=0.0d0
+ endif
+ if (fg_rank.eq.0) then
+ if (AFMlog.gt.0) then
+ call AFMforce(Eafmforce)
+ else if (selfguide.gt.0) then
+ call AFMvel(Eafmforce)
+ endif
+ endif
+ if (tubemode.eq.1) then
+ call calctube(etube)
+ else if (tubemode.eq.2) then
+ call calctube2(etube)
+ elseif (tubemode.eq.3) then
+ call calcnano(etube)
+ else
+ etube=0.0d0
+ endif
#ifdef TIMING
time_enecalc=time_enecalc+MPI_Wtime()-time00
energia(17)=estr
energia(20)=Uconst+Uconst_back
energia(21)=esccor
+ energia(22)=eliptran
+ energia(23)=Eafmforce
+ energia(24)=ethetacnstr
+ energia(25)=etube
! Here are the energies showed per procesor if the are more processors
! per molecule then we sum it up in sum_energy subroutine
! print *," Processor",myrank," calls SUM_ENERGY"
logical :: reduce
real(kind=8) :: evdw,evdw2,evdw2_14,ees,evdw1,ecorr,ecorr5,ecorr6
real(kind=8) :: eel_loc,eello_turn3,eello_turn4,eturn6,ebe,escloc
- real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,esccor,etot
+ real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,esccor,etot, &
+ eliptran,etube, Eafmforce,ethetacnstr
integer :: i
#ifdef MPI
integer :: ierr
estr=energia(17)
Uconst=energia(20)
esccor=energia(21)
+ eliptran=energia(22)
+ Eafmforce=energia(23)
+ ethetacnstr=energia(24)
+ etube=energia(25)
#ifdef SPLITELE
etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 &
+wang*ebe+wtor*etors+wscloc*escloc &
+wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
+wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 &
+wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d &
- +wbond*estr+Uconst+wsccor*esccor
+ +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube&
+ +Eafmforce+ethetacnstr
#else
etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) &
+wang*ebe+wtor*etors+wscloc*escloc &
+wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
+wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 &
+wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d &
- +wbond*estr+Uconst+wsccor*esccor
+ +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube&
+ +Eafmforce+ethetacnstr
+
#endif
energia(0)=etot
! detecting NaNQ
!el local variables
real(kind=8) :: etot,evdw,evdw2,ees,evdw1,ecorr,ecorr5,ecorr6,eel_loc
real(kind=8) :: eello_turn6,eello_turn3,eello_turn4,ebe,escloc
- real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,Uconst,esccor
+ real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,Uconst,esccor,eliptran,&
+ etube,ethetacnstr,Eafmforce
etot=energia(0)
evdw=energia(1)
estr=energia(17)
Uconst=energia(20)
esccor=energia(21)
+ eliptran=energia(22)
+ Eafmforce=energia(23)
+ ethetacnstr=energia(24)
+ etube=energia(25)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,&
estr,wbond,ebe,wang,&
ecorr,wcorr,&
ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,&
eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,&
- edihcnstr,ebr*nss,&
- Uconst,etot
+ edihcnstr,ethetacnstr,ebr*nss,&
+ Uconst,eliptran,wliptran,Eafmforce,etube,wtube,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)'/ &
+ 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ &
'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ &
'UCONST= ',1pE16.6,' (Constraint energy)'/ &
+ 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/&
+ 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ &
+ 'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/ &
'ETOT= ',1pE16.6,' (total)')
#else
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,&
ecorr,wcorr,&
ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,&
eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,&
- ebr*nss,Uconst,etot
+ ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, &
+ etube,wtube,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)'/ &
+ 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ &
'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ &
'UCONST=',1pE16.6,' (Constraint energy)'/ &
+ 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ &
+ 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ &
+ 'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/ &
'ETOT= ',1pE16.6,' (total)')
#endif
return
! include 'COMMON.NAMES'
! include 'COMMON.IOUNITS'
! include 'COMMON.CONTACTS'
- real(kind=8),dimension(3) :: gg
+ real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj
integer :: num_conti
!el local variables
integer :: i,itypi,iint,j,itypi1,itypj,k
! 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_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
evdwij=e1+e2
!d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
!d epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
! include 'COMMON.INTERACT'
! include 'COMMON.IOUNITS'
! include 'COMMON.NAMES'
- real(kind=8),dimension(3) :: gg
+ real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj
logical :: scheck
!el local variables
integer :: i,iint,j,itypi,itypi1,k,itypj
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_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
evdwij=e_augm+e1+e2
!d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
!d epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
! Calculate whole angle-dependent part of epsilon and contributions
! to its derivatives
fac=(rrij*sigsq)**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij
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_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
!d write (iout,'(2(a3,i3,2x),15(0pf7.3))')
!d & restyp(itypi),i,restyp(itypj),j,
!d & epsi,sigm,chi1,chi2,chip1,chip2,
real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
real(kind=8) :: evdw,sig0ij
real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
- dist_temp, dist_init
+ dist_temp, dist_init,aa,bb,ssgradlipi,ssgradlipj, &
+ sslipi,sslipj,faclip
integer :: ii
+ real(kind=8) :: fracinbuf
+
!cccc energy_dec=.false.
! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
! if (icall.eq.0) lprn=.false.
!el ind=0
do i=iatsc_s,iatsc_e
+!C print *,"I am in EVDW",i
itypi=iabs(itype(i))
+! if (i.ne.47) cycle
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
zi=dmod(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
+! print *, sslipi,ssgradlipi
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
'evdw',i,j,evdwij,' ss'
! if (energy_dec) write (iout,*) &
! 'evdw',i,j,evdwij,' ss'
+ 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
!el ind=ind+1
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
+! if (j.ne.78) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
! write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,&
if (yj.lt.0) yj=yj+boxysize
zj=dmod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
+! print *,"tu",xi,yi,zi,xj,yj,zj
+! print *,"tu2",j,j+nres,c(1,j),c(1,j+nres)
+! this fragment set correct epsilon for lipid phase
+ 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
+!------------------------------------------------
dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
xj_safe=xj
yj_safe=yj
!---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ faclip=fac
+ e1=fac*fac*aa!(itypi,itypj)
+ e2=fac*bb!(itypi,itypj)
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij*sss_ele_cut
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!(itypi,itypj)
write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
restyp(itypi),i,restyp(itypj),j, &
epsi,sigm,chi1,chi2,chip1,chip2, &
evdwij
endif
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
- 'evdw',i,j,evdwij !,"egb"
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2e10.2,e11.3)')&
+ 'evdw',i,j,evdwij,xi,xj,rij !,"egb"
+!C print *,i,j,c(1,i),c(1,j),c(2,i),c(2,j),c(3,i),c(3,j)
! if (energy_dec) write (iout,*) &
! 'evdw',i,j,evdwij
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
+!C Calculate the radial part of the gradient
+ gg_lipi(3)=eps1*(eps2rt*eps2rt)&
+ *(eps3rt*eps3rt)*sss_ele_cut/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
+
! print *,'before sc_grad', gg(1),gg(2),gg(3)
! Calculate angular part of the gradient.
call sc_grad
!---------------------------------------------------------------
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_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij+e_augm
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_aq(itypi,itypj)/&
+ bb_aq(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
restyp(itypi),i,restyp(itypj),j,&
epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),&
! include 'COMMON.NAMES'
! include 'COMMON.IOUNITS'
! include 'COMMON.CONTACTS'
- real(kind=8),dimension(3) :: gg
+ real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj
!d print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct
!el local variables
integer :: i,iint,j,itypi,itypi1,itypj,k
real(kind=8) :: auxvec(2),auxmat(2,2)
integer :: i,iti1,iti,k,l
real(kind=8) :: sin1,cos1,sin2,cos2,dwacos2,dwasin2
-
+! print *,"in set matrices"
!
! Compute the virtual-bond-torsional-angle dependent quantities needed
! to calculate the el-loc multibody terms of various order.
#else
do i=3,nres+1
#endif
+! print *,i,"i"
if (i .lt. nres+1) then
sin1=dsin(phi(i))
cos1=dcos(phi(i))
else
iti1=ntortyp+1
endif
+! print *,iti,i,"iti",iti1,itype(i-1),itype(i-2)
!d write (iout,*) '*******i',i,' iti1',iti
!d write (iout,*) 'b1',b1(:,iti)
!d write (iout,*) 'b2',b2(:,iti)
!el local variables
integer :: i,k,j
real(kind=8) :: ees,evdw1,eel_loc,eello_turn3,eello_turn4
- real(kind=8) :: fac,t_eelecij
+ real(kind=8) :: fac,t_eelecij,fracinbuf
!d write(iout,*) 'In EELEC'
+! print *,"IN EELEC"
!d do i=1,nloctyp
!d write(iout,*) 'Type',i
!d write(iout,*) 'B1',B1(:,i)
! write (iout,*) 'i',i,' fac',fac
enddo
endif
+! print *,wel_loc,"wel_loc",wcorr4,wcorr5,wcorr6,wturn3,wturn4, &
+! wturn6
if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 &
.or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or. &
wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
#ifdef TIMING
time01=MPI_Wtime()
#endif
+! print *, "before set matrices"
call set_matrices
+! print *, "after set matrices"
+
#ifdef TIMING
time_mat=time_mat+MPI_Wtime()-time01
#endif
endif
+! print *, "after set matrices"
!d do i=1,nres-1
!d write (iout,*) 'i=',i
!d do k=1,3
!
-
+! print *,"before iturn3 loop"
do i=iturn3_start,iturn3_end
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 &
.or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
+ xmedi=dmod(xmedi,boxxsize)
+ if (xmedi.lt.0) xmedi=xmedi+boxxsize
+ ymedi=dmod(ymedi,boxysize)
+ if (ymedi.lt.0) ymedi=ymedi+boxysize
+ zmedi=dmod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
num_conti=0
- call eelecij(i,i+2,ees,evdw1,eel_loc)
+ if ((zmedi.gt.bordlipbot) &
+ .and.(zmedi.lt.bordliptop)) then
+!C the energy transfer exist
+ if (zmedi.lt.buflipbot) then
+!C what fraction I am in
+ fracinbuf=1.0d0- &
+ ((zmedi-bordlipbot)/lipbufthick)
+!C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zmedi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zmedi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
+! print *,i,sslipi,ssgradlipi
+ call eelecij(i,i+2,ees,evdw1,eel_loc)
if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
num_cont_hb(i)=num_conti
enddo
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
+ xmedi=dmod(xmedi,boxxsize)
+ if (xmedi.lt.0) xmedi=xmedi+boxxsize
+ ymedi=dmod(ymedi,boxysize)
+ if (ymedi.lt.0) ymedi=ymedi+boxysize
+ zmedi=dmod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ if ((zmedi.gt.bordlipbot) &
+ .and.(zmedi.lt.bordliptop)) then
+!C the energy transfer exist
+ if (zmedi.lt.buflipbot) then
+!C what fraction I am in
+ fracinbuf=1.0d0- &
+ ((zmedi-bordlipbot)/lipbufthick)
+!C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zmedi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zmedi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
+
num_conti=num_cont_hb(i)
call eelecij(i,i+3,ees,evdw1,eel_loc)
if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) &
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
+ xmedi=dmod(xmedi,boxxsize)
+ if (xmedi.lt.0) xmedi=xmedi+boxxsize
+ ymedi=dmod(ymedi,boxysize)
+ if (ymedi.lt.0) ymedi=ymedi+boxysize
+ zmedi=dmod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ if ((zmedi.gt.bordlipbot) &
+ .and.(zmedi.lt.bordliptop)) then
+!C the energy transfer exist
+ if (zmedi.lt.buflipbot) then
+!C what fraction I am in
+ fracinbuf=1.0d0- &
+ ((zmedi-bordlipbot)/lipbufthick)
+!C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zmedi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zmedi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
+
! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
! include 'COMMON.VECTORS'
! include 'COMMON.FFIELD'
! include 'COMMON.TIME1'
- real(kind=8),dimension(3) :: ggg,gggp,gggm,erij,dcosb,dcosg
+ real(kind=8),dimension(3) :: ggg,gggp,gggm,erij,dcosb,dcosg,xtemp
real(kind=8),dimension(3,3) :: erder,uryg,urzg,vryg,vrzg
real(kind=8),dimension(2,2) :: acipa !el,a_temp
!el real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1
real(kind=8),dimension(4) :: muij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,rlocshield,fracinbuf
+ integer xshift,yshift,zshift,ilist,iresshield
!el integer :: num_conti,j1,j2
!el real(kind=8) :: a22,a23,a32,a33,dxi,dyi,dzi,dx_normi,dy_normi,&
!el dz_normi,xmedi,ymedi,zmedi
0.0d0,0.0d0,1.0d0/),shape(unmat))
! integer :: maxconts=nres/4
!el local variables
- integer :: k,i,j,iteli,itelj,kkk,l,kkll,m
+ integer :: k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap
real(kind=8) :: ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp
real(kind=8) :: ees,evdw1,eel_loc,aaa,bbb,ael3i
real(kind=8) :: dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj,&
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-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
+ 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
+
+ isubchap=0
+ dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ 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
+!C print *,i,j
+ 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
rrmij=1.0D0/rij
rij=dsqrt(rij)
+!C print *,xmedi,ymedi,zmedi,xj,yj,zj,boxxsize,rij
+ sss_ele_cut=sscale_ele(rij)
+ sss_ele_grad=sscagrad_ele(rij)
+! sss_ele_cut=1.0d0
+! sss_ele_grad=0.0d0
+! print *,sss_ele_cut,sss_ele_grad,&
+! (rij),r_cut_ele,rlamb_ele
+! if (sss_ele_cut.le.0.0) go to 128
+
rmij=1.0D0/rij
r3ij=rrmij*rmij
r6ij=r3ij*r3ij
evdwij=ev1+ev2
el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
el2=fac4*fac
- eesij=el1+el2
+! eesij=el1+el2
+ if (shield_mode.gt.0) then
+!C fac_shield(i)=0.4
+!C fac_shield(j)=0.6
+ el1=el1*fac_shield(i)**2*fac_shield(j)**2
+ el2=el2*fac_shield(i)**2*fac_shield(j)**2
+ eesij=(el1+el2)
+ ees=ees+eesij*sss_ele_cut
+!C FOR NOW SHIELD IS NOT USED WITH LIPSCALE
+!C & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+ else
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+ eesij=(el1+el2)
+ ees=ees+eesij &
+ *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)*sss_ele_cut
+!C print *,"TUCC",(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0
+ endif
+
! 12/26/95 - for the evaluation of multi-body H-bonding interactions
ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
- ees=ees+eesij
- evdw1=evdw1+evdwij
+! ees=ees+eesij*sss_ele_cut
+ evdw1=evdw1+evdwij*sss_ele_cut &
+ *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
!d write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
!d & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
!d & 1.0D0/dsqrt(rrmij),evdwij,eesij,
! Calculate contributions to the Cartesian gradient.
!
#ifdef SPLITELE
- facvdw=-6*rrmij*(ev1+evdwij)
- facel=-3*rrmij*(el1+eesij)
+ facvdw=-6*rrmij*(ev1+evdwij)*sss_ele_cut &
+ *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+ facel=-3*rrmij*(el1+eesij)*sss_ele_cut &
+ *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
fac1=fac
erij(1)=xj*rmij
erij(2)=yj*rmij
!
! Radial derivatives. First process both termini of the fragment (i,j)
!
- ggg(1)=facel*xj
- ggg(2)=facel*yj
- ggg(3)=facel*zj
+ ggg(1)=facel*xj+sss_ele_grad*rmij*eesij*xj* &
+ ((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+ ggg(2)=facel*yj+sss_ele_grad*rmij*eesij*yj* &
+ ((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+ ggg(3)=facel*zj+sss_ele_grad*rmij*eesij*zj* &
+ ((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
+ 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*sss_ele_cut
+ gshieldx(k,iresshield)=gshieldx(k,iresshield)+ &
+ rlocshield &
+ +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0 &
+ *sss_ele_cut
+ gshieldc(k,iresshield-1)=gshieldc(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)*eesij/fac_shield(j) &
+ *2.0*sss_ele_cut
+ gshieldx(k,iresshield)=gshieldx(k,iresshield)+ &
+ rlocshield &
+ +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0 &
+ *sss_ele_cut
+ 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 &
+ *sss_ele_cut
+
+ gshieldc(k,j)=gshieldc(k,j)+ &
+ grad_shield(k,j)*eesij/fac_shield(j)*2.0 &
+ *sss_ele_cut
+
+ gshieldc(k,i-1)=gshieldc(k,i-1)+ &
+ grad_shield(k,i)*eesij/fac_shield(i)*2.0 &
+ *sss_ele_cut
+
+ gshieldc(k,j-1)=gshieldc(k,j-1)+ &
+ grad_shield(k,j)*eesij/fac_shield(j)*2.0 &
+ *sss_ele_cut
+
+ enddo
+ endif
+
+
! do k=1,3
! ghalf=0.5D0*ggg(k)
! gelc(k,i)=gelc(k,i)+ghalf
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
enddo
+ gelc_long(3,j)=gelc_long(3,j)+ &
+ ssgradlipj*eesij/2.0d0*lipscale**2&
+ *sss_ele_cut
+
+ gelc_long(3,i)=gelc_long(3,i)+ &
+ ssgradlipi*eesij/2.0d0*lipscale**2&
+ *sss_ele_cut
+
+
!
! Loop over residues i+1 thru j-1.
!
!grad gelc(l,k)=gelc(l,k)+ggg(l)
!grad enddo
!grad enddo
- ggg(1)=facvdw*xj
- ggg(2)=facvdw*yj
- ggg(3)=facvdw*zj
+ ggg(1)=facvdw*xj+sss_ele_grad*rmij*evdwij*xj &
+ *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+ ggg(2)=facvdw*yj+sss_ele_grad*rmij*evdwij*yj &
+ *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+ ggg(3)=facvdw*zj+sss_ele_grad*rmij*evdwij*zj &
+ *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
! do k=1,3
! ghalf=0.5D0*ggg(k)
! gvdwpp(k,i)=gvdwpp(k,i)+ghalf
gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
enddo
-!
-! Loop over residues i+1 thru j-1.
+
+!C Lipidic part for scaling weight
+ gvdwpp(3,j)=gvdwpp(3,j)+ &
+ sss_ele_cut*ssgradlipj*evdwij/2.0d0*lipscale**2
+ gvdwpp(3,i)=gvdwpp(3,i)+ &
+ sss_ele_cut*ssgradlipi*evdwij/2.0d0*lipscale**2
+!! Loop over residues i+1 thru j-1.
!
!grad do k=i+1,j-1
!grad do l=1,3
!grad enddo
!grad enddo
#else
- facvdw=ev1+evdwij
- facel=el1+eesij
+ facvdw=(ev1+evdwij)*sss_ele_cut &
+ *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
+ facel=(el1+eesij)*sss_ele_cut
fac1=fac
fac=-3*rrmij*(facvdw+facvdw+facel)
erij(1)=xj*rmij
!
! Radial derivatives. First process both termini of the fragment (i,j)
!
- ggg(1)=fac*xj
- ggg(2)=fac*yj
- ggg(3)=fac*zj
+ ggg(1)=fac*xj+sss_ele_grad*rmij*(eesij+evdwij)*xj
+ ggg(2)=fac*yj+sss_ele_grad*rmij*(eesij+evdwij)*yj
+ ggg(3)=fac*zj+sss_ele_grad*rmij*(eesij+evdwij)*zj
! do k=1,3
! ghalf=0.5D0*ggg(k)
! gelc(k,i)=gelc(k,i)+ghalf
!grad enddo
!grad enddo
! 9/28/08 AL Gradient compotents will be summed only at the end
- ggg(1)=facvdw*xj
- ggg(2)=facvdw*yj
- ggg(3)=facvdw*zj
+ ggg(1)=facvdw*xj &
+ *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+ ggg(2)=facvdw*yj &
+ *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+ ggg(3)=facvdw*zj &
+ *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
do k=1,3
gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
enddo
+ gvdwpp(3,j)=gvdwpp(3,j)+ &
+ sss_ele_cut*ssgradlipj*evdwij/2.0d0*lipscale**2
+ gvdwpp(3,i)=gvdwpp(3,i)+ &
+ sss_ele_cut*ssgradlipi*evdwij/2.0d0*lipscale**2
+
#endif
!
! Angular part
!d print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
!d & (dcosg(k),k=1,3)
do k=1,3
- ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
+ ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*sss_ele_cut &
+ *fac_shield(i)**2*fac_shield(j)**2 &
+ *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
enddo
! do k=1,3
! ghalf=0.5D0*ggg(k)
do k=1,3
gelc(k,i)=gelc(k,i) &
+(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) &
- + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+ + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)&
+ *sss_ele_cut &
+ *fac_shield(i)**2*fac_shield(j)**2 &
+ *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
gelc(k,j)=gelc(k,j) &
+(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) &
- + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+ + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)&
+ *sss_ele_cut &
+ *fac_shield(i)**2*fac_shield(j)**2 &
+ *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
enddo
+
IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 &
.or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 &
.or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
! Contribution to the local-electrostatic energy coming from the i-j pair
eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) &
+a33*muij(4)
-! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+ endif
+ eel_loc_ij=eel_loc_ij &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+!C Now derivative over eel_loc
+ 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)&
+ *sss_ele_cut
+ gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+ &
+ rlocshield &
+ +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i) &
+ *sss_ele_cut
+
+ 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) &
+ *sss_ele_cut
+ gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+ &
+ rlocshield &
+ +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j) &
+ *sss_ele_cut
+
+ 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) &
+ *sss_ele_cut
+ gshieldc_ll(k,j)=gshieldc_ll(k,j)+ &
+ grad_shield(k,j)*eel_loc_ij/fac_shield(j) &
+ *sss_ele_cut
+ gshieldc_ll(k,i-1)=gshieldc_ll(k,i-1)+ &
+ grad_shield(k,i)*eel_loc_ij/fac_shield(i) &
+ *sss_ele_cut
+ gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+ &
+ grad_shield(k,j)*eel_loc_ij/fac_shield(j) &
+ *sss_ele_cut
+ enddo
+ endif
+
+
+! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
+! eel_loc_ij=0.0
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
'eelloc',i,j,eel_loc_ij
! if (energy_dec) write (iout,*) "a22",a22," a23",a23," a32",a32," a33",a33
! if (energy_dec) write (iout,*) "muij",muij
! write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
-
- eel_loc=eel_loc+eel_loc_ij
+
+ eel_loc=eel_loc+eel_loc_ij*sss_ele_cut
! Partial derivatives in virtual-bond dihedral angles gamma
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)
+ (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)) &
+ *sss_ele_cut &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
gel_loc_loc(j-1)=gel_loc_loc(j-1)+ &
- a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) &
- +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)
+ (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)) &
+ *sss_ele_cut &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
! Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
- 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)
+! do l=1,3
+! ggg(1)=(agg(1,1)*muij(1)+ &
+! agg(1,2)*muij(2)+agg(1,3)*muij(3)+agg(1,4)*muij(4)) &
+! *sss_ele_cut &
+! +eel_loc_ij*sss_ele_grad*rmij*xj
+! ggg(2)=(agg(2,1)*muij(1)+ &
+! agg(2,2)*muij(2)+agg(2,3)*muij(3)+agg(2,4)*muij(4)) &
+! *sss_ele_cut &
+! +eel_loc_ij*sss_ele_grad*rmij*yj
+! ggg(3)=(agg(3,1)*muij(1)+ &
+! agg(3,2)*muij(2)+agg(3,3)*muij(3)+agg(3,4)*muij(4)) &
+! *sss_ele_cut &
+! +eel_loc_ij*sss_ele_grad*rmij*zj
+ xtemp(1)=xj
+ xtemp(2)=yj
+ xtemp(3)=zj
+
+ 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))&
+ *sss_ele_cut &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) &
+ +eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
+
+
gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
!grad ghalf=0.5d0*ggg(l)
!grad gel_loc(l,i)=gel_loc(l,i)+ghalf
!grad gel_loc(l,j)=gel_loc(l,j)+ghalf
enddo
+ gel_loc_long(3,j)=gel_loc_long(3,j)+ &
+ ssgradlipj*eel_loc_ij/2.0d0*lipscale/ &
+ ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
+
+ gel_loc_long(3,i)=gel_loc_long(3,i)+ &
+ ssgradlipi*eel_loc_ij/2.0d0*lipscale/ &
+ ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
+
!grad do k=i+1,j2
!grad do l=1,3
!grad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
!grad enddo
! Remaining derivatives of eello
do l=1,3
- gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+ &
- aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
- 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)
- 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)
- 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)
+ 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))&
+ *sss_ele_cut &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+!+eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
+ 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))&
+ *sss_ele_cut &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+!+eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
+ 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))&
+ *sss_ele_cut &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+!+eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
+ 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))&
+ *sss_ele_cut &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+!+eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
enddo
ENDIF
! Change 12/26/95 to calculate four-body contributions to H-bonding energy
else
ees0pij=0
endif
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0d0
+ fac_shield(j)=1.0d0
+ else
+ ees0plist(num_conti,i)=j
+ endif
! ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
if (ees0tmp.gt.0) then
ees0mij=0
endif
! ees0mij=0.0D0
- ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
- ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
+ ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij) &
+ *sss_ele_cut &
+ *fac_shield(i)*fac_shield(j)
+
+ ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) &
+ *sss_ele_cut &
+ *fac_shield(i)*fac_shield(j)
+
! Diagnostics. Comment out or remove after debugging!
! ees0p(num_conti,i)=0.5D0*fac3*ees0pij
! ees0m(num_conti,i)=0.5D0*fac3*ees0mij
gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
enddo
- gggp(1)=gggp(1)+ees0pijp*xj
- gggp(2)=gggp(2)+ees0pijp*yj
- gggp(3)=gggp(3)+ees0pijp*zj
- gggm(1)=gggm(1)+ees0mijp*xj
- gggm(2)=gggm(2)+ees0mijp*yj
- gggm(3)=gggm(3)+ees0mijp*zj
+ gggp(1)=gggp(1)+ees0pijp*xj &
+ +ees0p(num_conti,i)/sss_ele_cut*rmij*xj*sss_ele_grad
+ gggp(2)=gggp(2)+ees0pijp*yj &
+ +ees0p(num_conti,i)/sss_ele_cut*rmij*yj*sss_ele_grad
+ gggp(3)=gggp(3)+ees0pijp*zj &
+ +ees0p(num_conti,i)/sss_ele_cut*rmij*zj*sss_ele_grad
+
+ gggm(1)=gggm(1)+ees0mijp*xj &
+ +ees0m(num_conti,i)/sss_ele_cut*rmij*xj*sss_ele_grad
+
+ gggm(2)=gggm(2)+ees0mijp*yj &
+ +ees0m(num_conti,i)/sss_ele_cut*rmij*yj*sss_ele_grad
+
+ gggm(3)=gggm(3)+ees0mijp*zj &
+ +ees0m(num_conti,i)/sss_ele_cut*rmij*zj*sss_ele_grad
+
! Derivatives due to the contact function
gacont_hbr(1,num_conti,i)=fprimcont*xj
gacont_hbr(2,num_conti,i)=fprimcont*yj
!grad ghalfm=0.5D0*gggm(k)
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)
+ + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) &
+ *sss_ele_cut*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)
- gacontp_hb3(k,num_conti,i)=gggp(k)
+ + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)&
+ *sss_ele_cut*fac_shield(i)*fac_shield(j)
+
+ gacontp_hb3(k,num_conti,i)=gggp(k) &
+ *sss_ele_cut*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)
+ + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) &
+ *sss_ele_cut*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)
- gacontm_hb3(k,num_conti,i)=gggm(k)
+ + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) &
+ *sss_ele_cut*fac_shield(i)*fac_shield(j)
+
+ gacontm_hb3(k,num_conti,i)=gggm(k) &
+ *sss_ele_cut*fac_shield(i)*fac_shield(j)
+
enddo
! Diagnostics. Comment out or remove after debugging!
!diag do k=1,3
enddo
endif
endif
+ 128 continue
! t_eelecij=t_eelecij+MPI_Wtime()-time00
return
end subroutine eelecij
!el dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,&
!el num_conti,j1,j2
!el local variables
- integer :: i,j,l
- real(kind=8) :: eello_turn3
+ integer :: i,j,l,k,ilist,iresshield
+ real(kind=8) :: eello_turn3,zj,fracinbuf,eello_t3, rlocshield
j=i+2
! write (iout,*) "eturn3",i,j,j1,j2
+ zj=(c(3,j)+c(3,j+1))/2.0d0
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ if ((zj.lt.0)) write (*,*) "CHUJ"
+ if ((zj.gt.bordlipbot) &
+ .and.(zj.lt.bordliptop)) then
+!C the energy transfer exist
+ if (zj.lt.buflipbot) then
+!C what fraction I am in
+ fracinbuf=1.0d0- &
+ ((zj-bordlipbot)/lipbufthick)
+!C lipbufthick is thickenes of lipid buffore
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zj.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipj=1.0d0
+ ssgradlipj=0.0
+ endif
+ else
+ sslipj=0.0d0
+ ssgradlipj=0.0
+ endif
+
a_temp(1,1)=a22
a_temp(1,2)=a23
a_temp(2,1)=a32
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))
- eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2))
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0d0
+ fac_shield(j)=1.0d0
+ endif
+
+ eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2)) &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+ eello_t3= &
+ 0.5d0*(pizda(1,1)+pizda(2,2)) &
+ *fac_shield(i)*fac_shield(j)
+
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2))
+ 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)
+ 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)
+ 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
+
!d write (2,*) 'i,',i,' j',j,'eello_turn3',
!d & 0.5d0*(pizda(1,1)+pizda(2,2)),
!d & ' eello_turn3_num',4*eello_turn3_num
call matmat2(EUgder(1,1,i+1),EUg(1,1,i+2),auxmat2(1,1))
call transpose2(auxmat2(1,1),auxmat3(1,1))
call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1))
- gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2))
+ gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2))&
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
! Derivatives in gamma(i+1)
call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1))
call transpose2(auxmat2(1,1),auxmat3(1,1))
call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1))
gel_loc_turn3(i+1)=gel_loc_turn3(i+1) &
- +0.5d0*(pizda(1,1)+pizda(2,2))
+ +0.5d0*(pizda(1,1)+pizda(2,2)) &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
! Cartesian derivatives
do l=1,3
! ghalf1=0.5d0*agg(l,1)
a_temp(2,2)=aggi(l,4)!+ghalf4
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))
+ +0.5d0*(pizda(1,1)+pizda(2,2)) &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
a_temp(1,1)=aggi1(l,1)!+agg(l,1)
a_temp(1,2)=aggi1(l,2)!+agg(l,2)
a_temp(2,1)=aggi1(l,3)!+agg(l,3)
a_temp(2,2)=aggi1(l,4)!+agg(l,4)
call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1) &
- +0.5d0*(pizda(1,1)+pizda(2,2))
+ +0.5d0*(pizda(1,1)+pizda(2,2)) &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
a_temp(1,1)=aggj(l,1)!+ghalf1
a_temp(1,2)=aggj(l,2)!+ghalf2
a_temp(2,1)=aggj(l,3)!+ghalf3
a_temp(2,2)=aggj(l,4)!+ghalf4
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))
+ +0.5d0*(pizda(1,1)+pizda(2,2)) &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
a_temp(1,1)=aggj1(l,1)
a_temp(1,2)=aggj1(l,2)
a_temp(2,1)=aggj1(l,3)
a_temp(2,2)=aggj1(l,4)
call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
gcorr3_turn(l,j1)=gcorr3_turn(l,j1) &
- +0.5d0*(pizda(1,1)+pizda(2,2))
+ +0.5d0*(pizda(1,1)+pizda(2,2)) &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
enddo
+ gshieldc_t3(3,i)=gshieldc_t3(3,i)+ &
+ ssgradlipi*eello_t3/4.0d0*lipscale
+ gshieldc_t3(3,j)=gshieldc_t3(3,j)+ &
+ ssgradlipj*eello_t3/4.0d0*lipscale
+ gshieldc_t3(3,i-1)=gshieldc_t3(3,i-1)+ &
+ ssgradlipi*eello_t3/4.0d0*lipscale
+ gshieldc_t3(3,j-1)=gshieldc_t3(3,j-1)+ &
+ ssgradlipj*eello_t3/4.0d0*lipscale
+
return
end subroutine eturn3
!-----------------------------------------------------------------------------
!el dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,&
!el num_conti,j1,j2
!el local variables
- integer :: i,j,iti1,iti2,iti3,l
- real(kind=8) :: eello_turn4,s1,s2,s3
+ integer :: i,j,iti1,iti2,iti3,l,k,ilist,iresshield
+ real(kind=8) :: eello_turn4,s1,s2,s3,zj,fracinbuf,eello_t4,&
+ rlocshield
j=i+3
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!d call checkint_turn4(i,a_temp,eello_turn4_num)
! write (iout,*) "eturn4 i",i," j",j," j1",j1," j2",j2
+ zj=(c(3,j)+c(3,j+1))/2.0d0
+ 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
+
a_temp(1,1)=a22
a_temp(1,2)=a23
a_temp(2,1)=a32
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))
- eello_turn4=eello_turn4-(s1+s2+s3)
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+ endif
+
+ eello_turn4=eello_turn4-(s1+s2+s3) &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+ eello_t4=-(s1+s2+s3) &
+ *fac_shield(i)*fac_shield(j)
+!C Now derivative over shield:
+ 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)
+ 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)
+ 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
+
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
'eturn4',i,j,-(s1+s2+s3)
!d write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
s1=scalar2(b1(1,iti2),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)
+ gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3) &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
! 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(ae3(1,1),e2tder(1,1),auxmat(1,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)
+ gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3) &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
! 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),e2t(1,1),auxmat3(1,1))
call matmat2(auxmat3(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
- gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3)
+ gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3) &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
! Cartesian derivatives
! Derivatives of this turn contributions in DC(i+2)
if (j.lt.nres-1) then
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
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)
+ gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3)&
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
enddo
endif
! Remaining derivatives of this turn contribution
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))
- gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3)
+ gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3) &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+
a_temp(1,1)=aggi1(l,1)
a_temp(1,2)=aggi1(l,2)
a_temp(2,1)=aggi1(l,3)
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))
- gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3)
+ gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3) &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+
a_temp(1,1)=aggj(l,1)
a_temp(1,2)=aggj(l,2)
a_temp(2,1)=aggj(l,3)
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))
- gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3)
+ gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3) &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+
a_temp(1,1)=aggj1(l,1)
a_temp(1,2)=aggj1(l,2)
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))
! write (iout,*) "s1",s1," s2",s2," s3",s3," s1+s2+s3",s1+s2+s3
- gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3)
+ gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) &
+ *fac_shield(i)*fac_shield(j) &
+ *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
enddo
+ gshieldc_t4(3,i)=gshieldc_t4(3,i)+ &
+ ssgradlipi*eello_t4/4.0d0*lipscale
+ gshieldc_t4(3,j)=gshieldc_t4(3,j)+ &
+ ssgradlipj*eello_t4/4.0d0*lipscale
+ gshieldc_t4(3,i-1)=gshieldc_t4(3,i-1)+ &
+ ssgradlipi*eello_t4/4.0d0*lipscale
+ gshieldc_t4(3,j-1)=gshieldc_t4(3,j-1)+ &
+ ssgradlipj*eello_t4/4.0d0*lipscale
+
return
end subroutine eturn4
!-----------------------------------------------------------------------------
! include 'COMMON.CONTROL'
real(kind=8),dimension(3) :: ggg
!el local variables
- integer :: i,iint,j,k,iteli,itypj
+ integer :: i,iint,j,k,iteli,itypj,subchap
real(kind=8) :: evdw2,evdw2_14,xi,yi,zi,xj,yj,zj,rrij,fac,&
- e1,e2,evdwij
+ e1,e2,evdwij,rij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
+ integer xshift,yshift,zshift
evdw2=0.0D0
evdw2_14=0.0d0
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))
+ 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)
! yj=c(2,nres+j)-yi
! zj=c(3,nres+j)-zi
! Uncomment following three lines for Ca-p interactions
- xj=c(1,j)-xi
- yj=c(2,j)-yi
- zj=c(3,j)-zi
- rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
- 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
- endif
- evdwij=e1+e2
- evdw2=evdw2+evdwij
-! if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)') &
-! 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),&
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
- 'evdw2',i,j,evdwij
-!
-! Calculate contributions to the gradient in the virtual-bond and SC vectors.
-!
- fac=-(evdwij+e1)*rrij
- ggg(1)=xj*fac
- ggg(2)=yj*fac
- ggg(3)=zj*fac
-!grad if (j.lt.i) then
+! 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)
+ 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
+ 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)
+ rij=dsqrt(1.0d0/rrij)
+ sss_ele_cut=sscale_ele(rij)
+ sss_ele_grad=sscagrad_ele(rij)
+! print *,sss_ele_cut,sss_ele_grad,&
+! (rij),r_cut_ele,rlamb_ele
+ if (sss_ele_cut.le.0.0) cycle
+ 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)*sss_ele_cut
+ endif
+ evdwij=e1+e2
+ evdw2=evdw2+evdwij*sss_ele_cut
+! if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)') &
+! 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),&
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
+ 'evdw2',i,j,evdwij
+!
+! Calculate contributions to the gradient in the virtual-bond and SC vectors.
+!
+ fac=-(evdwij+e1)*rrij*sss_ele_cut
+ fac=fac+evdwij*sss_ele_grad/rij/expon
+ ggg(1)=xj*fac
+ ggg(2)=yj*fac
+ ggg(3)=zj*fac
+!grad if (j.lt.i) then
!d write (iout,*) 'j<i'
! Uncomment following three lines for SC-p interactions
! do k=1,3
ehpb=ehpb+2*eij
!d write (iout,*) "eij",eij
endif
+ 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
+ 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
+ if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, &
+ 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,*) "beta nmr",
+!c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ else
+ dd=dist(ii,jj)
+ 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,*) "beta 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
+ do j=1,3
+ ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
+ ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
+ enddo
+ do k=1,3
+ ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+ ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+ enddo
else
-! Calculate the distance between the two points and its difference from the
-! target distance.
- dd=dist(ii,jj)
- rdis=dd-dhpb(i)
-! Get the force constant corresponding to this distance.
- waga=forcon(i)
-! Calculate the contribution to energy.
- ehpb=ehpb+waga*rdis*rdis
-!
-! Evaluate gradient.
-!
- fac=waga*rdis/dd
-!d print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
-!d & ' waga=',waga,' fac=',fac
- do j=1,3
- ggg(j)=fac*(c(j,jj)-c(j,ii))
- enddo
-!d print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
-! If this is a SC-SC distance, we need to calculate the contributions to the
-! Cartesian gradient in the SC vectors (ghpbx).
- if (iii.lt.ii) then
+ dd=dist(ii,jj)
+ 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
+ if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, &
+ 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
-!grad do j=iii,jjj-1
-!grad do k=1,3
-!grad ghpbc(k,j)=ghpbc(k,j)+ggg(k)
-!grad enddo
-!grad enddo
- do k=1,3
- ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
- ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
- enddo
+ endif
+!cgrad do j=iii,jjj-1
+!cgrad do k=1,3
+!cgrad ghpbc(k,j)=ghpbc(k,j)+ggg(k)
+!cgrad enddo
+!cgrad enddo
+ do k=1,3
+ ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+ ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+ enddo
endif
enddo
- ehpb=0.5D0*ehpb
+ if (constr_dist.ne.11) ehpb=0.5D0*ehpb
+
return
end subroutine edis
!-----------------------------------------------------------------------------
end subroutine theteng
#else
!-----------------------------------------------------------------------------
- subroutine ebend(etheta)
+ subroutine ebend(etheta,ethetacnstr)
!
! Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
! angles gamma and its derivatives in consecutive thetas and gammas.
!el local variables
integer :: i,k,iblock,ityp1,ityp2,ityp3,l,m
real(kind=8) :: dethetai,dephii,dephii1,theti2,phii,phii1,ethetai
- real(kind=8) :: aux,etheta,ccl,ssl,scl,csl
+ real(kind=8) :: aux,etheta,ccl,ssl,scl,csl,ethetacnstr
+! local variables for constrains
+ real(kind=8) :: difi,thetiii
+ integer itheta
etheta=0.0D0
do i=ithet_start,ithet_end
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
gloc(nphi+i-2,icg)=wang*dethetai
enddo
+!-----------thete constrains
+! if (tor_mode.ne.2) then
+ ethetacnstr=0.0d0
+!C print *,ithetaconstr_start,ithetaconstr_end,"TU"
+ do i=ithetaconstr_start,ithetaconstr_end
+ 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
+ if (energy_dec) then
+ write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", &
+ i,itheta,rad2deg*thetiii, &
+ rad2deg*theta_constr0(i), rad2deg*theta_drange(i), &
+ rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, &
+ gloc(itheta+nphi-2,icg)
+ endif
+ enddo
+! endif
+
return
end subroutine ebend
#endif
real(kind=8),dimension(3) :: gx,gx1
logical :: lprn
!el local variables
- integer :: i,j,k,l,jj,kk,ll
+ integer :: i,j,k,l,jj,kk,ll,ilist,m, iresshield
real(kind=8) :: coeffp,coeffm,eij,ekl,ees0pij,ees0pkl,ees0mij,&
ees0mkl,ees,coeffpees0pij,coeffmees0mij,&
- coeffpees0pkl,coeffmees0mkl,gradlongij,gradlongkl
+ coeffpees0pkl,coeffmees0mkl,gradlongij,gradlongkl, &
+ rlocshield
lprn=.false.
eij=facont_hb(jj,i)
!grad enddo
! write (iout,*) "ehbcorr",ekont*ees
ehbcorr=ekont*ees
+ 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)
+ 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)
+ 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)
+ 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)
+ 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
+ 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
return
end function ehbcorr
#ifdef MOMENT
#ifdef MPI
include 'mpif.h'
#endif
- real(kind=8),dimension(3,nres) :: gradbufc,gradbufx,gradbufc_sum,&
+ real(kind=8),dimension(3,-1:nres) :: gradbufc,gradbufx,gradbufc_sum,&
gloc_scbuf !(3,maxres)
real(kind=8),dimension(4*nres) :: glocbuf !(4*maxres)
call flush(iout)
#endif
#ifdef SPLITELE
- do i=1,nct
+ do i=0,nct
do j=1,3
gradbufc(j,i)=wsc*gvdwc(j,i)+ &
wscp*(gvdwc_scp(j,i)+gvdwc_scpp(j,i))+ &
wcorr5*gradcorr5_long(j,i)+ &
wcorr6*gradcorr6_long(j,i)+ &
wturn6*gcorr6_turn_long(j,i)+ &
- wstrain*ghpbc(j,i)
+ wstrain*ghpbc(j,i) &
+ +wliptran*gliptranc(j,i) &
+ +gradafm(j,i) &
+ +welec*gshieldc(j,i) &
+ +wcorr*gshieldc_ec(j,i) &
+ +wturn3*gshieldc_t3(j,i)&
+ +wturn4*gshieldc_t4(j,i)&
+ +wel_loc*gshieldc_ll(j,i)&
+ +wtube*gg_tube(j,i)
+
+
+
enddo
enddo
#else
- do i=1,nct
+ do i=0,nct
do j=1,3
gradbufc(j,i)=wsc*gvdwc(j,i)+ &
wscp*(gvdwc_scp(j,i)+gvdwc_scpp(j,i))+ &
wcorr5*gradcorr5_long(j,i)+ &
wcorr6*gradcorr6_long(j,i)+ &
wturn6*gcorr6_turn_long(j,i)+ &
- wstrain*ghpbc(j,i)
+ wstrain*ghpbc(j,i) &
+ +wliptran*gliptranc(j,i) &
+ +gradafm(j,i) &
+ +welec*gshieldc(j,i)&
+ +wcorr*gshieldc_ec(j,i) &
+ +wturn4*gshieldc_t4(j,i) &
+ +wel_loc*gshieldc_ll(j,i)&
+ +wtube*gg_tube(j,i)
+
+
+
enddo
enddo
#endif
enddo
call flush(iout)
#endif
- do i=1,nres
+ do i=0,nres
do j=1,3
gradbufc_sum(j,i)=gradbufc(j,i)
enddo
#ifdef TIMING
! time_allreduce=time_allreduce+MPI_Wtime()-time00
#endif
- do i=nnt,nres
+ do i=0,nres
do k=1,3
gradbufc(k,i)=0.0d0
enddo
do j=1,3
gradbufc(j,nres-1)=gradbufc_sum(j,nres)
enddo
- do i=nres-2,nnt,-1
+ do i=nres-2,-1,-1
do j=1,3
gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1)
enddo
call flush(iout)
#endif
!el#undef DEBUG
- do i=1,nres
+ do i=-1,nres
do j=1,3
gradbufc_sum(j,i)=gradbufc(j,i)
gradbufc(j,i)=0.0d0
do j=1,3
gradbufc(j,nres-1)=gradbufc_sum(j,nres)
enddo
- do i=nres-2,nnt,-1
+ do i=nres-2,-1,-1
do j=1,3
gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1)
enddo
!el if (.not.allocated(gradx)) allocate(gradx(3,nres,2)) !(3,maxres,2)
!el if (.not.allocated(gradc)) allocate(gradc(3,nres,2)) !(3,maxres,2)
!el-----------------
- do i=1,nct
+ do i=-1,nct
do j=1,3
#ifdef SPLITELE
gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ &
wcorr6*gradcorr6(j,i)+ &
wturn6*gcorr6_turn(j,i)+ &
wsccor*gsccorc(j,i) &
- +wscloc*gscloc(j,i)
+ +wscloc*gscloc(j,i) &
+ +wliptran*gliptranc(j,i) &
+ +gradafm(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) &
+ +wtube*gg_tube(j,i)
+
+
#else
gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ &
wel_loc*gel_loc(j,i)+ &
wcorr6*gradcorr6(j,i)+ &
wturn6*gcorr6_turn(j,i)+ &
wsccor*gsccorc(j,i) &
- +wscloc*gscloc(j,i)
+ +wscloc*gscloc(j,i) &
+ +gradafm(j,i) &
+ +wliptran*gliptranc(j,i) &
+ +welec*gshieldc(j,i) &
+ +welec*gshieldc_loc(j,) &
+ +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) &
+ +wtube*gg_tube(j,i)
+
+
+
#endif
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*gsccorx(j,i) &
- +wscloc*gsclocx(j,i)
+ +wscloc*gsclocx(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)&
+ +wtube*gg_tube_sc(j,i)
+
+
enddo
enddo
#ifdef DEBUG
call MPI_Barrier(FG_COMM,IERR)
time_barrier_g=time_barrier_g+MPI_Wtime()-time00
time00=MPI_Wtime()
- call MPI_Reduce(gradbufc(1,1),gradc(1,1,icg),3*nres,&
+ call MPI_Reduce(gradbufc(1,0),gradc(1,0,icg),3*nres+3,&
MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
call MPI_Reduce(gradbufx(1,1),gradx(1,1,icg),3*nres,&
MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
do k=1,3
gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*sss_ele_cut
!C print *,'gg',k,gg(k)
- enddo
+ enddo
+! print *,i,j,gg_lipi(3),gg_lipj(3),sss_ele_cut
! write (iout,*) "gg",(gg(k),k=1,3)
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 &
*sss_ele_cut
- gvdwx(k,j)=gvdwx(k,j)+gg(k) &
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)&
+(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
+eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv &
*sss_ele_cut
!grad enddo
!grad enddo
do l=1,3
- gvdwc(l,i)=gvdwc(l,i)-gg(l)
- gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
enddo
return
end subroutine sc_grad
! Check the gradient of the virtual-bond and SC vectors in the internal
! coordinates.
!
- aincr=1.0d-7
- aincr2=5.0d-8
+ aincr=1.0d-6
+ aincr2=5.0d-7
call cartder
write (iout,'(a)') '**************** dx/dalpha'
write (iout,'(a)')
nf=0
nfl=0
call zerograd
- aincr=1.0D-7
- print '(a)','CG processor',me,' calling CHECK_CART.'
+ aincr=1.0D-5
+ print '(a)','CG processor',me,' calling CHECK_CART.',aincr
nf=0
icall=0
call geom_to_var(nvar,x)
! call intcartderiv
! call checkintcartgrad
call zerograd
- aincr=1.0D-6
- write(iout,*) 'Calling CHECK_ECARTINT.'
+ aincr=2.0D-5
+ write(iout,*) 'Calling CHECK_ECARTINT.',aincr
nf=0
icall=0
call geom_to_var(nvar,x)
endif
return
end function sscale
+ real(kind=8) function sscale_grad(r)
+! include "COMMON.SPLITELE"
+ real(kind=8) :: r,gamm
+ if(r.lt.r_cut-rlamb) then
+ sscale_grad=0.0d0
+ else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+ gamm=(r-(r_cut-rlamb))/rlamb
+ sscale_grad=gamm*(6*gamm-6.0d0)/rlamb
+ else
+ sscale_grad=0d0
+ endif
+ return
+ end function sscale_grad
+
!!!!!!!!!! PBCSCALE
real(kind=8) function sscale_ele(r)
! include "COMMON.SPLITELE"
endif
return
end function sscagrad_ele
+ real(kind=8) function sscalelip(r)
+ real(kind=8) r,gamm
+ sscalelip=1.0d0+r*r*(2.0d0*r-3.0d0)
+ return
+ end function sscalelip
+!C-----------------------------------------------------------------------
+ real(kind=8) function sscagradlip(r)
+ real(kind=8) r,gamm
+ sscagradlip=r*(6.0d0*r-6.0d0)
+ return
+ end function sscagradlip
+
!!!!!!!!!!!!!!!
!-----------------------------------------------------------------------------
subroutine elj_long(evdw)
! include 'COMMON.IOUNITS'
! include 'COMMON.CONTACTS'
real(kind=8),parameter :: accur=1.0d-10
- real(kind=8),dimension(3) :: gg
+ real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj
!el local variables
integer :: i,iint,j,k,itypi,itypi1,itypj
real(kind=8) :: xi,yi,zi,xj,yj,zj,rij,sss,rrij,fac,eps0ij
rrij=1.0D0/rij
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
evdwij=e1+e2
evdw=evdw+(1.0d0-sss)*evdwij
!
! include 'COMMON.IOUNITS'
! include 'COMMON.CONTACTS'
real(kind=8),parameter :: accur=1.0d-10
- real(kind=8),dimension(3) :: gg
+ real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj
!el local variables
integer :: i,iint,j,k,itypi,itypi1,itypj,num_conti
real(kind=8) :: xi,yi,zi,xj,yj,zj,rij,sss,rrij,fac,eps0ij
rrij=1.0D0/rij
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
evdwij=e1+e2
evdw=evdw+sss*evdwij
!
! include 'COMMON.INTERACT'
! include 'COMMON.IOUNITS'
! include 'COMMON.NAMES'
- real(kind=8),dimension(3) :: gg
+ real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj
logical :: scheck
!el local variables
integer :: i,iint,j,k,itypi,itypi1,itypj
if (sss.lt.1.0d0) then
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_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
evdwij=e_augm+e1+e2
!d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
!d epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
! include 'COMMON.INTERACT'
! include 'COMMON.IOUNITS'
! include 'COMMON.NAMES'
- real(kind=8),dimension(3) :: gg
+ real(kind=8),dimension(3) :: gg,gg_lipi,gg_lipj
logical :: scheck
!el local variables
integer :: i,iint,j,k,itypi,itypi1,itypj
if (sss.gt.0.0d0) then
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_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
evdwij=e_augm+e1+e2
!d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
!d epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
! Calculate whole angle-dependent part of epsilon and contributions
! to its derivatives
fac=(rrij*sigsq)**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij*(1.0d0-sss)
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_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
!d write (iout,'(2(a3,i3,2x),15(0pf7.3))')
!d & restyp(itypi),i,restyp(itypj),j,
!d & epsi,sigm,chi1,chi2,chip1,chip2,
! Calculate whole angle-dependent part of epsilon and contributions
! to its derivatives
fac=(rrij*sigsq)**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij*sss
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_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
!d write (iout,'(2(a3,i3,2x),15(0pf7.3))')
!d & restyp(itypi),i,restyp(itypj),j,
!d & epsi,sigm,chi1,chi2,chip1,chip2,
!el local variables
integer :: iint,itypi,itypi1,itypj,subchap
real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig,sig0ij,rij_shift
- real(kind=8) :: sss,e1,e2,evdw
+ real(kind=8) :: sss,e1,e2,evdw,sss_grad
real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
- dist_temp, dist_init
+ dist_temp, dist_init,aa,bb,fracinbuf,sslipi,sslipj,&
+ ssgradlipi,ssgradlipj
+
evdw=0.0D0
!cccc energy_dec=.false.
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)
!
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
+! call dyn_ssbond_ene(i,j,evdwij)
+! evdw=evdw+evdwij
+! if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
+! 'evdw',i,j,evdwij,' ss'
+! if (energy_dec) write (iout,*) &
+! 'evdw',i,j,evdwij,' ss'
+! 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
!el ind=ind+1
itypj=itype(j)
if (itypj.eq.ntyp1) cycle
if (yj.lt.0) yj=yj+boxysize
zj=mod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
+ if ((zj.gt.bordlipbot) &
+ .and.(zj.lt.bordliptop)) then
+!C the energy transfer exist
+ if (zj.lt.buflipbot) then
+!C what fraction I am in
+ fracinbuf=1.0d0- &
+ ((zj-bordlipbot)/lipbufthick)
+!C lipbufthick is thickenes of lipid buffore
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zj.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipj=1.0d0
+ ssgradlipj=0.0
+ endif
+ else
+ sslipj=0.0d0
+ ssgradlipj=0.0
+ endif
+ 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
+
dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
xj_safe=xj
yj_safe=yj
sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
if (sss_ele_cut.le.0.0) cycle
if (sss.lt.1.0d0) then
!---------------------------------------------------------------
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
evdw=evdw+evdwij*(1.0d0-sss)*sss_ele_cut
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_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
restyp(itypi),i,restyp(itypj),j,&
epsi,sigm,chi1,chi2,chip1,chip2,&
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac
- fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
- /sigma(itypi,itypj)*rij
+ fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
+ /sigma(itypi,itypj)*rij-sss_grad/(1.0-sss)*rij &
+ /sigmaii(itypi,itypj))
! fac=0.0d0
! Calculate the radial part of the gradient
gg(1)=xj*fac
gg(3)=zj*fac
! Calculate angular part of the gradient.
call sc_grad_scale(1.0d0-sss)
+ ENDIF !mask_dyn_ss
endif
enddo ! j
enddo ! iint
!el local variables
integer :: iint,itypi,itypi1,itypj,subchap
real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig0ij,sig
- real(kind=8) :: sss,e1,e2,evdw,rij_shift
+ real(kind=8) :: sss,e1,e2,evdw,rij_shift,sss_grad
real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
- dist_temp, dist_init
+ dist_temp, dist_init,aa,bb,fracinbuf,sslipi,sslipj,&
+ ssgradlipi,ssgradlipj
evdw=0.0D0
!cccc energy_dec=.false.
! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
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)
!
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
+ call dyn_ssbond_ene(i,j,evdwij)
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
+ 'evdw',i,j,evdwij,' ss'
+ 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
+
+! if (energy_dec) write (iout,*) &
+! 'evdw',i,j,evdwij,' ss'
+ ELSE
!el ind=ind+1
itypj=itype(j)
if (itypj.eq.ntyp1) cycle
if (yj.lt.0) yj=yj+boxysize
zj=mod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
+ if ((zj.gt.bordlipbot) &
+ .and.(zj.lt.bordliptop)) then
+!C the energy transfer exist
+ if (zj.lt.buflipbot) then
+!C what fraction I am in
+ fracinbuf=1.0d0- &
+ ((zj-bordlipbot)/lipbufthick)
+!C lipbufthick is thickenes of lipid buffore
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zj.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipj=1.0d0
+ ssgradlipj=0.0
+ endif
+ else
+ sslipj=0.0d0
+ ssgradlipj=0.0
+ endif
+ 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
+
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
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+ sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
if (sss_ele_cut.le.0.0) cycle
!---------------------------------------------------------------
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
! write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
! & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij*sss
+ evdw=evdw+evdwij*sss*sss_ele_cut
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_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
restyp(itypi),i,restyp(itypj),j,&
epsi,sigm,chi1,chi2,chip1,chip2,&
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac
- fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
- /sigma(itypi,itypj)*rij
+ fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
+ /sigma(itypi,itypj)*rij+sss_grad/sss*rij &
+ /sigmaii(itypi,itypj))
! fac=0.0d0
! Calculate the radial part of the gradient
! Calculate angular part of the gradient.
call sc_grad_scale(sss)
endif
+ ENDIF !mask_dyn_ss
enddo ! j
enddo ! iint
enddo ! i
!---------------------------------------------------------------
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_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
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_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
restyp(itypi),i,restyp(itypj),j,&
epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),&
!---------------------------------------------------------------
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_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+(evdwij+e_augm)*sss
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_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj)
write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
restyp(itypi),i,restyp(itypj),j,&
epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),&
#ifdef TIMING
time01=MPI_Wtime()
#endif
+! print *, "before set matrices"
call set_matrices
+! print *,"after set martices"
#ifdef TIMING
time_mat=time_mat+MPI_Wtime()-time01
#endif
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
+ xmedi=dmod(xmedi,boxxsize)
+ if (xmedi.lt.0) xmedi=xmedi+boxxsize
+ ymedi=dmod(ymedi,boxysize)
+ if (ymedi.lt.0) ymedi=ymedi+boxysize
+ zmedi=dmod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
num_conti=0
call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
+ xmedi=dmod(xmedi,boxxsize)
+ if (xmedi.lt.0) xmedi=xmedi+boxxsize
+ ymedi=dmod(ymedi,boxysize)
+ if (ymedi.lt.0) ymedi=ymedi+boxysize
+ zmedi=dmod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
num_conti=num_cont_hb(i)
call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) &
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
+ xmedi=dmod(xmedi,boxxsize)
+ if (xmedi.lt.0) xmedi=xmedi+boxxsize
+ ymedi=dmod(ymedi,boxysize)
+ if (ymedi.lt.0) ymedi=ymedi+boxysize
+ zmedi=dmod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
! include 'COMMON.VECTORS'
! include 'COMMON.FFIELD'
! include 'COMMON.TIME1'
- real(kind=8),dimension(3) :: ggg,gggp,gggm,erij,dcosb,dcosg
+ real(kind=8),dimension(3) :: ggg,gggp,gggm,erij,dcosb,dcosg,xtemp
real(kind=8),dimension(3,3) :: erder,uryg,urzg,vryg,vrzg
real(kind=8),dimension(2,2) :: acipa !el,a_temp
!el real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1
real(kind=8),dimension(4) :: muij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,sss_grad
+ integer xshift,yshift,zshift
+
!el integer :: num_conti,j1,j2
!el real(kind=8) :: a22,a23,a32,a33,dxi,dyi,dzi,dx_normi,dy_normi,&
!el dz_normi,xmedi,ymedi,zmedi
0.0d0,1.0d0,0.0d0,&
0.0d0,0.0d0,1.0d0/),shape(unmat))
!el local variables
- integer :: i,j,k,l,iteli,itelj,kkk,kkll,m
+ integer :: i,j,k,l,iteli,itelj,kkk,kkll,m,isubchap
real(kind=8) :: aaa,bbb,ael6i,ael3i,dxj,dyj,dzj
real(kind=8) :: xj,yj,zj,rij,rrmij,rmij,sss,r3ij,r6ij,fac
real(kind=8) :: cosa,cosb,cosg,ev1,ev2,fac3,fac4,evdwij
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-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
+ isubchap=0
+ dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ 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
+!C print *,i,j
+ 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
rrmij=1.0D0/rij
rij=dsqrt(rij)
rmij=1.0D0/rij
! For extracting the short-range part of Evdwpp
sss=sscale(rij/rpp(iteli,itelj))
+ sss_ele_cut=sscale_ele(rij)
+ sss_ele_grad=sscagrad_ele(rij)
+ sss_grad=sscale_grad((rij/rpp(iteli,itelj)))
+! sss_ele_cut=1.0d0
+! sss_ele_grad=0.0d0
+ if (sss_ele_cut.le.0.0) go to 128
r3ij=rrmij*rmij
r6ij=r3ij*r3ij
eesij=el1+el2
! 12/26/95 - for the evaluation of multi-body H-bonding interactions
ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
- ees=ees+eesij
- evdw1=evdw1+evdwij*(1.0d0-sss)
+ ees=ees+eesij*sss_ele_cut
+ evdw1=evdw1+evdwij*(1.0d0-sss)*sss_ele_cut
!d write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
!d & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
!d & 1.0D0/dsqrt(rrmij),evdwij,eesij,
! Calculate contributions to the Cartesian gradient.
!
#ifdef SPLITELE
- facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
- facel=-3*rrmij*(el1+eesij)
+ facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss_ele_cut
+ facel=-3*rrmij*(el1+eesij)*sss_ele_cut
fac1=fac
erij(1)=xj*rmij
erij(2)=yj*rmij
!
! Radial derivatives. First process both termini of the fragment (i,j)
!
- ggg(1)=facel*xj
- ggg(2)=facel*yj
- ggg(3)=facel*zj
+ ggg(1)=facel*xj+sss_ele_grad*rmij*eesij*xj
+ ggg(2)=facel*yj+sss_ele_grad*rmij*eesij*yj
+ ggg(3)=facel*zj+sss_ele_grad*rmij*eesij*zj
! do k=1,3
! ghalf=0.5D0*ggg(k)
! gelc(k,i)=gelc(k,i)+ghalf
!grad gelc(l,k)=gelc(l,k)+ggg(l)
!grad enddo
!grad enddo
- ggg(1)=facvdw*xj
- ggg(2)=facvdw*yj
- ggg(3)=facvdw*zj
+ ggg(1)=facvdw*xj+sss_ele_grad*rmij*evdwij*xj*(1.0d0-sss) &
+ -evdwij*sss_ele_cut/rij*sss_grad/rpp(iteli,itelj)*xj
+ ggg(2)=facvdw*yj+sss_ele_grad*rmij*evdwij*yj*(1.0d0-sss) &
+ -evdwij*sss_ele_cut/rij*sss_grad/rpp(iteli,itelj)*yj
+ ggg(3)=facvdw*zj+sss_ele_grad*rmij*evdwij*zj*(1.0d0-sss) &
+ -evdwij*sss_ele_cut/rij*sss_grad/rpp(iteli,itelj)*zj
! do k=1,3
! ghalf=0.5D0*ggg(k)
! gvdwpp(k,i)=gvdwpp(k,i)+ghalf
!grad enddo
!grad enddo
#else
- facvdw=ev1+evdwij*(1.0d0-sss)
- facel=el1+eesij
+ facvdw=(ev1+evdwij)*(1.0d0-sss)*sss_ele_cut
+ facel=(el1+eesij)*sss_ele_cut
fac1=fac
fac=-3*rrmij*(facvdw+facvdw+facel)
erij(1)=xj*rmij
!d print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
!d & (dcosg(k),k=1,3)
do k=1,3
- ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
+ ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k) )*sss_ele_cut
enddo
! do k=1,3
! ghalf=0.5D0*ggg(k)
do k=1,3
gelc(k,i)=gelc(k,i) &
+(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) &
- + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+ + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)&
+ *sss_ele_cut
gelc(k,j)=gelc(k,j) &
+(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) &
- + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+ + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)&
+ *sss_ele_cut
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
enddo
'eelloc',i,j,eel_loc_ij
! write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) !d
- eel_loc=eel_loc+eel_loc_ij
+ eel_loc=eel_loc+eel_loc_ij*sss_ele_cut
! Partial derivatives in virtual-bond dihedral angles gamma
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)
+ (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)) &
+ *sss_ele_cut
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)
+ (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)) &
+ *sss_ele_cut
+ xtemp(1)=xj
+ xtemp(2)=yj
+ xtemp(3)=zj
+
! Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
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)
+ ggg(l)=(agg(l,1)*muij(1)+ &
+ agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))&
+ *sss_ele_cut &
+ +eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
+
gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
!grad ghalf=0.5d0*ggg(l)
!grad enddo
! Remaining derivatives of eello
do l=1,3
- gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+ &
- aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
- 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)
- 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)
- 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)
+ 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))&
+ *sss_ele_cut
+
+ 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))&
+ *sss_ele_cut
+
+ 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))&
+ *sss_ele_cut
+
+ 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))&
+ *sss_ele_cut
+
enddo
ENDIF
! Change 12/26/95 to calculate four-body contributions to H-bonding energy
ees0mij=0
endif
! ees0mij=0.0D0
- ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
- ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
+ ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij) &
+ *sss_ele_cut
+
+ ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) &
+ *sss_ele_cut
+
! Diagnostics. Comment out or remove after debugging!
! ees0p(num_conti,i)=0.5D0*fac3*ees0pij
! ees0m(num_conti,i)=0.5D0*fac3*ees0mij
gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
enddo
- gggp(1)=gggp(1)+ees0pijp*xj
- gggp(2)=gggp(2)+ees0pijp*yj
- gggp(3)=gggp(3)+ees0pijp*zj
- gggm(1)=gggm(1)+ees0mijp*xj
- gggm(2)=gggm(2)+ees0mijp*yj
- gggm(3)=gggm(3)+ees0mijp*zj
+! gggp(1)=gggp(1)+ees0pijp*xj
+! gggp(2)=gggp(2)+ees0pijp*yj
+! gggp(3)=gggp(3)+ees0pijp*zj
+! gggm(1)=gggm(1)+ees0mijp*xj
+! gggm(2)=gggm(2)+ees0mijp*yj
+! gggm(3)=gggm(3)+ees0mijp*zj
+ gggp(1)=gggp(1)+ees0pijp*xj &
+ +ees0p(num_conti,i)/sss_ele_cut*rmij*xj*sss_ele_grad
+ gggp(2)=gggp(2)+ees0pijp*yj &
+ +ees0p(num_conti,i)/sss_ele_cut*rmij*yj*sss_ele_grad
+ gggp(3)=gggp(3)+ees0pijp*zj &
+ +ees0p(num_conti,i)/sss_ele_cut*rmij*zj*sss_ele_grad
+
+ gggm(1)=gggm(1)+ees0mijp*xj &
+ +ees0m(num_conti,i)/sss_ele_cut*rmij*xj*sss_ele_grad
+
+ gggm(2)=gggm(2)+ees0mijp*yj &
+ +ees0m(num_conti,i)/sss_ele_cut*rmij*yj*sss_ele_grad
+
+ gggm(3)=gggm(3)+ees0mijp*zj &
+ +ees0m(num_conti,i)/sss_ele_cut*rmij*zj*sss_ele_grad
+
! Derivatives due to the contact function
gacont_hbr(1,num_conti,i)=fprimcont*xj
gacont_hbr(2,num_conti,i)=fprimcont*yj
!
!grad ghalfp=0.5D0*gggp(k)
!grad ghalfm=0.5D0*gggm(k)
- 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)
- 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)
- gacontp_hb3(k,num_conti,i)=gggp(k)
- 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)
- 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)
- gacontm_hb3(k,num_conti,i)=gggm(k)
- enddo
- ENDIF ! wcorr
- endif ! num_conti.le.maxconts
- endif ! fcont.gt.0
- endif ! j.gt.i+1
- if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
- do k=1,4
- do l=1,3
- ghalf=0.5d0*agg(l,k)
- aggi(l,k)=aggi(l,k)+ghalf
+! 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)
+! 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)
+! gacontp_hb3(k,num_conti,i)=gggp(k)
+! 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)
+! 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)
+! gacontm_hb3(k,num_conti,i)=gggm(k)
+ 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) &
+ *sss_ele_cut
+
+ 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)&
+ *sss_ele_cut
+
+ gacontp_hb3(k,num_conti,i)=gggp(k) &
+ *sss_ele_cut
+
+ 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) &
+ *sss_ele_cut
+
+ 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) &
+ *sss_ele_cut
+
+ gacontm_hb3(k,num_conti,i)=gggm(k) &
+ *sss_ele_cut
+
+ enddo
+ ENDIF ! wcorr
+ endif ! num_conti.le.maxconts
+ endif ! fcont.gt.0
+ endif ! j.gt.i+1
+ if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
+ do k=1,4
+ do l=1,3
+ ghalf=0.5d0*agg(l,k)
+ aggi(l,k)=aggi(l,k)+ghalf
aggi1(l,k)=aggi1(l,k)+agg(l,k)
aggj(l,k)=aggj(l,k)+ghalf
enddo
enddo
endif
endif
+ 128 continue
! t_eelecij=t_eelecij+MPI_Wtime()-time00
return
end subroutine eelecij_scale
real(kind=8) :: scal_el=0.5d0
#endif
!el local variables
- integer :: i,j,k,iteli,itelj,num_conti
+ integer :: i,j,k,iteli,itelj,num_conti,isubchap
real(kind=8) :: dxi,dyi,dzi,dxj,dyj,dzj,aaa,bbb
real(kind=8) :: xj,yj,zj,rij,rrmij,sss,r3ij,r6ij,evdw1,&
dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,&
dx_normj,dy_normj,dz_normj,rmij,ev1,ev2,evdwij,facvdw
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,sss_grad
+ integer xshift,yshift,zshift
+
evdw1=0.0D0
! write (iout,*) "iatel_s_vdw",iatel_s_vdw,
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
+ xmedi=dmod(xmedi,boxxsize)
+ if (xmedi.lt.0) xmedi=xmedi+boxxsize
+ ymedi=dmod(ymedi,boxysize)
+ if (ymedi.lt.0) ymedi=ymedi+boxysize
+ zmedi=dmod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
num_conti=0
! write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
! & ' ielend',ielend_vdw(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-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
+ isubchap=0
+ dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ 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
+!C print *,i,j
+ 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
rrmij=1.0D0/rij
rij=dsqrt(rij)
sss=sscale(rij/rpp(iteli,itelj))
+ sss_ele_cut=sscale_ele(rij)
+ sss_ele_grad=sscagrad_ele(rij)
+ sss_grad=sscale_grad((rij/rpp(iteli,itelj)))
+ if (sss_ele_cut.le.0.0) cycle
if (sss.gt.0.0d0) then
rmij=1.0D0/rij
r3ij=rrmij*rmij
if (energy_dec) then
write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
endif
- evdw1=evdw1+evdwij*sss
+ evdw1=evdw1+evdwij*sss*sss_ele_cut
!
! Calculate contributions to the Cartesian gradient.
!
- facvdw=-6*rrmij*(ev1+evdwij)*sss
- ggg(1)=facvdw*xj
- ggg(2)=facvdw*yj
- ggg(3)=facvdw*zj
+ facvdw=-6*rrmij*(ev1+evdwij)*sss*sss_ele_cut
+! ggg(1)=facvdw*xj
+! ggg(2)=facvdw*yj
+! ggg(3)=facvdw*zj
+ ggg(1)=facvdw*xj+sss_ele_grad*rmij*evdwij*xj*sss &
+ +evdwij*sss_ele_cut/rij*sss_grad/rpp(iteli,itelj)*xj
+ ggg(2)=facvdw*yj+sss_ele_grad*rmij*evdwij*yj*sss &
+ +evdwij*sss_ele_cut/rij*sss_grad/rpp(iteli,itelj)*yj
+ ggg(3)=facvdw*zj+sss_ele_grad*rmij*evdwij*zj*sss &
+ +evdwij*sss_ele_cut/rij*sss_grad/rpp(iteli,itelj)*zj
+
do k=1,3
gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
! include 'COMMON.CONTROL'
real(kind=8),dimension(3) :: ggg
!el local variables
- integer :: i,iint,j,k,iteli,itypj
- real(kind=8) :: xi,yi,zi,xj,yj,zj,rrij,sss,fac,e1,e2
+ integer :: i,iint,j,k,iteli,itypj,subchap
+ real(kind=8) :: xi,yi,zi,xj,yj,zj,rrij,sss,fac,e1,e2,sss_grad,rij
real(kind=8) :: evdw2,evdw2_14,evdwij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
+
evdw2=0.0D0
evdw2_14=0.0d0
!d print '(a)','Enter ESCP'
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))
+ 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)
! yj=c(2,nres+j)-yi
! zj=c(3,nres+j)-zi
! 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)
+ 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
+ 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)
- sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
-
+ rij=dsqrt(1.0d0/rrij)
+ sss_ele_cut=sscale_ele(rij)
+ sss_ele_grad=sscagrad_ele(rij)
+! print *,sss_ele_cut,sss_ele_grad,&
+! (rij),r_cut_ele,rlamb_ele
+ if (sss_ele_cut.le.0.0) cycle
+ sss=sscale((rij/rscp(itypj,iteli)))
+ sss_grad=sscale_grad(rij/rscp(itypj,iteli))
if (sss.lt.1.0d0) then
fac=rrij**expon2
if (iabs(j-i) .le. 2) then
e1=scal14*e1
e2=scal14*e2
- evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
+ evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)*sss_ele_cut
endif
evdwij=e1+e2
- evdw2=evdw2+evdwij*(1.0d0-sss)
+ evdw2=evdw2+evdwij*(1.0d0-sss)*sss_ele_cut
if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))') &
'evdw2',i,j,sss,evdwij
!
! Calculate contributions to the gradient in the virtual-bond and SC vectors.
!
- fac=-(evdwij+e1)*rrij*(1.0d0-sss)
+ fac=-(evdwij+e1)*rrij*(1.0d0-sss)*sss_ele_cut
+ fac=fac+evdwij*sss_ele_grad/rij/expon*(1.0d0-sss)&
+ -evdwij*sss_ele_cut/rij/expon*sss_grad/rscp(itypj,iteli)
ggg(1)=xj*fac
ggg(2)=yj*fac
ggg(3)=zj*fac
! include 'COMMON.CONTROL'
real(kind=8),dimension(3) :: ggg
!el local variables
- integer :: i,iint,j,k,iteli,itypj
- real(kind=8) :: xi,yi,zi,xj,yj,zj,rrij,sss,fac,e1,e2
+ integer :: i,iint,j,k,iteli,itypj,subchap
+ real(kind=8) :: xi,yi,zi,xj,yj,zj,rrij,sss,fac,e1,e2,sss_grad,rij
real(kind=8) :: evdw2,evdw2_14,evdwij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
+
evdw2=0.0D0
evdw2_14=0.0d0
!d print '(a)','Enter ESCP'
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))
+ 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)
! yj=c(2,nres+j)-yi
! zj=c(3,nres+j)-zi
! Uncomment following three lines for Ca-p interactions
- xj=c(1,j)-xi
- yj=c(2,j)-yi
- zj=c(3,j)-zi
- rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
-
- sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
+! 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)
+ 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
+ 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)
+ rij=dsqrt(1.0d0/rrij)
+ sss_ele_cut=sscale_ele(rij)
+ sss_ele_grad=sscagrad_ele(rij)
+! print *,sss_ele_cut,sss_ele_grad,&
+! (rij),r_cut_ele,rlamb_ele
+ if (sss_ele_cut.le.0.0) cycle
+ sss=sscale(rij/rscp(itypj,iteli))
+ sss_grad=sscale_grad(rij/rscp(itypj,iteli))
if (sss.gt.0.0d0) then
fac=rrij**expon2
if (iabs(j-i) .le. 2) then
e1=scal14*e1
e2=scal14*e2
- evdw2_14=evdw2_14+(e1+e2)*sss
+ evdw2_14=evdw2_14+(e1+e2)*sss*sss_ele_cut
endif
evdwij=e1+e2
- evdw2=evdw2+evdwij*sss
+ evdw2=evdw2+evdwij*sss*sss_ele_cut
if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))') &
'evdw2',i,j,sss,evdwij
!
! Calculate contributions to the gradient in the virtual-bond and SC vectors.
!
- fac=-(evdwij+e1)*rrij*sss
+ fac=-(evdwij+e1)*rrij*sss*sss_ele_cut
+ fac=fac+evdwij*sss_ele_grad/rij/expon*sss &
+ +evdwij*sss_ele_cut/rij/expon*sss_grad/rscp(itypj,iteli)
+
ggg(1)=xj*fac
ggg(2)=yj*fac
ggg(3)=zj*fac
!el local variables
integer :: i,nres6
real(kind=8) :: evdw,evdw1,evdw2,evdw2_14,esccor,etors_d,etors
- real(kind=8) :: ehpb,escloc,estr,ebe,edihcnstr
+ real(kind=8) :: ehpb,escloc,estr,ebe,edihcnstr,ethetacnstr
nres6=6*nres
! write(iout,'(a,i2)')'Calling etotal_short ipot=',ipot
!
! Calculate the virtual-bond-angle energy.
!
- call ebend(ebe)
+ call ebend(ebe,ethetacnstr)
!
! Calculate the SC local energy.
!
endif
return
end function gnmr1prim
-!-----------------------------------------------------------------------------
+!----------------------------------------------------------------------------
+ real(kind=8) function rlornmr1(y,ymin,ymax,sigma)
+ real(kind=8) y,ymin,ymax,sigma
+ real(kind=8) wykl /4.0d0/
+ if (y.lt.ymin) then
+ rlornmr1=(ymin-y)**wykl/((ymin-y)**wykl+sigma**wykl)
+ else if (y.gt.ymax) then
+ rlornmr1=(y-ymax)**wykl/((y-ymax)**wykl+sigma**wykl)
+ else
+ rlornmr1=0.0d0
+ endif
+ return
+ end function rlornmr1
+!------------------------------------------------------------------------------
+ real(kind=8) function rlornmr1prim(y,ymin,ymax,sigma)
+ real(kind=8) y,ymin,ymax,sigma
+ real(kind=8) wykl /4.0d0/
+ if (y.lt.ymin) then
+ rlornmr1prim=-(ymin-y)**(wykl-1)*sigma**wykl*wykl/ &
+ ((ymin-y)**wykl+sigma**wykl)**2
+ else if (y.gt.ymax) then
+ rlornmr1prim=(y-ymax)**(wykl-1)*sigma**wykl*wykl/ &
+ ((y-ymax)**wykl+sigma**wykl)**2
+ else
+ rlornmr1prim=0.0d0
+ endif
+ return
+ end function rlornmr1prim
+
real(kind=8) function harmonic(y,ymax)
! implicit none
real(kind=8) :: y,ymax
#ifdef DEBUG
write (iout,*) "gcart, gxcart, gloc before int_to_cart"
#endif
- do i=1,nct
+ do i=0,nct
do j=1,3
gcart(j,i)=gradc(j,i,icg)
gxcart(j,i)=gradx(j,i,icg)
#ifdef DEBUG
write (iout,*) "CARGRAD"
#endif
- do i=nres,1,-1
+ do i=nres,0,-1
do j=1,3
gcart(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
! gcart_new(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
! include 'COMMON.SCCOR'
!
!el local variables
- integer :: i,j,intertyp
+ integer :: i,j,intertyp,k
! Initialize Cartesian-coordinate gradient
!
! if (.not.allocated(gradx)) allocate(gradx(3,nres,2)) !(3,maxres,2)
! allocate(gloc_sc(3,nres,10)) !(3,0:maxres2,10)maxres2=2*maxres
!elwrite(iout,*) "icg",icg
- do i=1,nres
+ do i=-1,nres
do j=1,3
gvdwx(j,i)=0.0D0
gradx_scp(j,i)=0.0D0
gradx(j,i,icg)=0.0d0
gscloc(j,i)=0.0d0
gsclocx(j,i)=0.0d0
+ gliptran(j,i)=0.0d0
+ gliptranx(j,i)=0.0d0
+ gliptranc(j,i)=0.0d0
+ gshieldx(j,i)=0.0d0
+ gshieldc(j,i)=0.0d0
+ gshieldc_loc(j,i)=0.0d0
+ gshieldx_ec(j,i)=0.0d0
+ gshieldc_ec(j,i)=0.0d0
+ gshieldc_loc_ec(j,i)=0.0d0
+ gshieldx_t3(j,i)=0.0d0
+ gshieldc_t3(j,i)=0.0d0
+ gshieldc_loc_t3(j,i)=0.0d0
+ gshieldx_t4(j,i)=0.0d0
+ gshieldc_t4(j,i)=0.0d0
+ gshieldc_loc_t4(j,i)=0.0d0
+ gshieldx_ll(j,i)=0.0d0
+ gshieldc_ll(j,i)=0.0d0
+ gshieldc_loc_ll(j,i)=0.0d0
+ gg_tube(j,i)=0.0d0
+ gg_tube_sc(j,i)=0.0d0
+ gradafm(j,i)=0.0d0
do intertyp=1,3
gloc_sc(intertyp,i,icg)=0.0d0
enddo
enddo
enddo
+ do i=1,nres
+ do j=1,maxcontsshi
+ shield_list(j,i)=0
+ do k=1,3
+!C print *,i,j,k
+ grad_shield_side(k,j,i)=0.0d0
+ grad_shield_loc(k,j,i)=0.0d0
+ enddo
+ enddo
+ ishield_list(i)=0
+ enddo
+
!
! Initialize the gradient of local energy terms.
!
! Obtaining the gamma derivatives from sine derivative
if (phi(i).gt.-pi4.and.phi(i).le.pi4.or. &
phi(i).gt.pi34.and.phi(i).le.pi.or. &
- phi(i).gt.-pi.and.phi(i).le.-pi34) then
+ phi(i).ge.-pi.and.phi(i).le.-pi34) then
call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
call vecpr(dc_norm(1,i-3),dc_norm(1,i-1),vp2)
call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
ljXs=sig-sig0ij
ljA=eps1*eps2rt**2*eps3rt**2
- ljB=ljA*bb(itypi,itypj)
- ljA=ljA*aa(itypi,itypj)
- ljxm=ljXs+(-2.0D0*aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+ ljB=ljA*bb_aq(itypi,itypj)
+ ljA=ljA*aa_aq(itypi,itypj)
+ ljxm=ljXs+(-2.0D0*aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
ssXs=d0cm
deltat1=1.0d0-om1
! Stop and plot energy and derivative as a function of distance
if (checkstop) then
ssm=ssC-0.25D0*ssB*ssB/ssA
- ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
+ ljm=-0.25D0*ljB*bb_aq(itypi,itypj)/aa_aq(itypi,itypj)
if (ssm.lt.ljm .and. &
dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
nicheck=1000
havebond=.false.
ljd=rij-ljXs
fac=(1.0D0/ljd)**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
eij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=eij*eps3rt
eps3der=eij*eps2rt
eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
else
havebond=.false.
- ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
- d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB
+ ljm=-0.25D0*ljB*bb_aq(itypi,itypj)/aa_aq(itypi,itypj)
+ d_ljm(1)=-0.5D0*bb_aq(itypi,itypj)/aa_aq(itypi,itypj)*ljB
d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt- &
alf12/eps3rt)
return
end subroutine dyn_ssbond_ene
+!--------------------------------------------------------------------------
+ subroutine triple_ssbond_ene(resi,resj,resk,eij)
+! implicit none
+! Includes
+ use calc_data
+ use comm_sschecks
+! include 'DIMENSIONS'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.VAR'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CALC'
+#ifndef CLUST
+#ifndef WHAM
+ use MD_data
+! include 'COMMON.MD'
+! use MD, only: totT,t_bath
+#endif
+#endif
+ double precision h_base
+ external h_base
+
+!c Input arguments
+ integer resi,resj,resk,m,itypi,itypj,itypk
+
+!c Output arguments
+ double precision eij,eij1,eij2,eij3
+
+!c Local variables
+ logical havebond
+!c integer itypi,itypj,k,l
+ double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
+ double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
+ double precision xik,yik,zik,xjk,yjk,zjk,dxk,dyk,dzk
+ double precision sig0ij,ljd,sig,fac,e1,e2
+ double precision dcosom1(3),dcosom2(3),ed
+ double precision pom1,pom2
+ double precision ljA,ljB,ljXs
+ double precision d_ljB(1:3)
+ double precision ssA,ssB,ssC,ssXs
+ double precision ssxm,ljxm,ssm,ljm
+ double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
+ eij=0.0
+ if (dtriss.eq.0) return
+ i=resi
+ j=resj
+ k=resk
+!C write(iout,*) resi,resj,resk
+ itypi=itype(i)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+ dsci_inv=vbld_inv(i+nres)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ itypj=itype(j)
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ dscj_inv=vbld_inv(j+nres)
+ itypk=itype(k)
+ xk=c(1,nres+k)
+ yk=c(2,nres+k)
+ zk=c(3,nres+k)
+
+ dxk=dc_norm(1,nres+k)
+ dyk=dc_norm(2,nres+k)
+ dzk=dc_norm(3,nres+k)
+ dscj_inv=vbld_inv(k+nres)
+ xij=xj-xi
+ xik=xk-xi
+ xjk=xk-xj
+ yij=yj-yi
+ yik=yk-yi
+ yjk=yk-yj
+ zij=zj-zi
+ zik=zk-zi
+ zjk=zk-zj
+ rrij=(xij*xij+yij*yij+zij*zij)
+ rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
+ rrik=(xik*xik+yik*yik+zik*zik)
+ rik=dsqrt(rrik)
+ rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
+ rjk=dsqrt(rrjk)
+!C there are three combination of distances for each trisulfide bonds
+!C The first case the ith atom is the center
+!C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
+!C distance y is second distance the a,b,c,d are parameters derived for
+!C this problem d parameter was set as a penalty currenlty set to 1.
+ if ((iabs(j-i).le.2).or.(iabs(i-k).le.2)) then
+ eij1=0.0d0
+ else
+ eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**6+ctriss)
+ endif
+!C second case jth atom is center
+ if ((iabs(j-i).le.2).or.(iabs(j-k).le.2)) then
+ eij2=0.0d0
+ else
+ eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**6+ctriss)
+ endif
+!C the third case kth atom is the center
+ if ((iabs(i-k).le.2).or.(iabs(j-k).le.2)) then
+ eij3=0.0d0
+ else
+ eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**6+ctriss)
+ endif
+!C eij2=0.0
+!C eij3=0.0
+!C eij1=0.0
+ eij=eij1+eij2+eij3
+!C write(iout,*)i,j,k,eij
+!C The energy penalty calculated now time for the gradient part
+!C derivative over rij
+ fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5) &
+ -eij2**2/dtriss*(2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)
+ gg(1)=xij*fac/rij
+ gg(2)=yij*fac/rij
+ gg(3)=zij*fac/rij
+ do m=1,3
+ gvdwx(m,i)=gvdwx(m,i)-gg(m)
+ gvdwx(m,j)=gvdwx(m,j)+gg(m)
+ enddo
+
+ do l=1,3
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ enddo
+!C now derivative over rik
+ fac=-eij1**2/dtriss* &
+ (-2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5) &
+ -eij3**2/dtriss*(2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
+ gg(1)=xik*fac/rik
+ gg(2)=yik*fac/rik
+ gg(3)=zik*fac/rik
+ do m=1,3
+ gvdwx(m,i)=gvdwx(m,i)-gg(m)
+ gvdwx(m,k)=gvdwx(m,k)+gg(m)
+ enddo
+ do l=1,3
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)
+ gvdwc(l,k)=gvdwc(l,k)+gg(l)
+ enddo
+!C now derivative over rjk
+ fac=-eij2**2/dtriss* &
+ (-2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)- &
+ eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
+ gg(1)=xjk*fac/rjk
+ gg(2)=yjk*fac/rjk
+ gg(3)=zjk*fac/rjk
+ do m=1,3
+ gvdwx(m,j)=gvdwx(m,j)-gg(m)
+ gvdwx(m,k)=gvdwx(m,k)+gg(m)
+ enddo
+ do l=1,3
+ gvdwc(l,j)=gvdwc(l,j)-gg(l)
+ gvdwc(l,k)=gvdwc(l,k)+gg(l)
+ enddo
+ return
+ end subroutine triple_ssbond_ene
+
+
+
!-----------------------------------------------------------------------------
real(kind=8) function h_base(x,deriv)
! A smooth function going 0->1 in range [0,1]
diff=newnss-nss
!mc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
-
+! print *,newnss,nss,maxdim
do i=1,nss
found=.false.
+! print *,newnss
do j=1,newnss
+!! print *,j
if (idssb(i).eq.newihpb(j) .and. &
jdssb(i).eq.newjhpb(j)) found=.true.
enddo
#ifndef CLUST
#ifndef WHAM
+! write(iout,*) "found",found,i,j
if (.not.found.and.fg_rank.eq.0) &
write(iout,'(a15,f12.2,f8.1,2i5)') &
"SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
do i=1,newnss
found=.false.
do j=1,nss
+! print *,i,j
if (newihpb(i).eq.idssb(j) .and. &
newjhpb(i).eq.jdssb(j)) found=.true.
enddo
#ifndef CLUST
#ifndef WHAM
+! write(iout,*) "found",found,i,j
if (.not.found.and.fg_rank.eq.0) &
write(iout,'(a15,f12.2,f8.1,2i5)') &
"SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
return
end subroutine dyn_set_nss
+! Lipid transfer energy function
+ subroutine Eliptransfer(eliptran)
+!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
+ real(kind=8) :: fracinbuf,eliptran,sslip,positi,ssgradlip
+ integer :: i
+ eliptran=0.0
+! print *, "I am in eliptran"
+ do i=ilip_start,ilip_end
+!C do i=1,1
+ if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1).or.(i.eq.nres))&
+ cycle
+
+ positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize))
+ if (positi.le.0.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
+
+!C print *,"doing sccale for lower part"
+!C print *,i,sslip,fracinbuf,ssgradlip
+ 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
+ if (energy_dec) write(iout,*) i,"eliptran=",eliptran,positi,sslip
+ enddo
+! here starts the side chain transfer
+ do i=ilip_start,ilip_end
+ 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
+ if (energy_dec) write(iout,*) i,"eliptran=",eliptran
+ enddo
+ return
+ end subroutine Eliptransfer
+!----------------------------------NANO FUNCTIONS
+!C-----------------------------------------------------------------------
+!C-----------------------------------------------------------
+!C This subroutine is to mimic the histone like structure but as well can be
+!C utilizet to nanostructures (infinit) small modification has to be used to
+!C make it finite (z gradient at the ends has to be changes as well as the x,y
+!C gradient has to be modified at the ends
+!C The energy function is Kihara potential
+!C E=4esp*((sigma/(r-r0))^12 - (sigma/(r-r0))^6)
+!C 4eps is depth of well sigma is r_minimum r is distance from center of tube
+!C and r0 is the excluded size of nanotube (can be set to 0 if we want just a
+!C simple Kihara potential
+ subroutine calctube(Etube)
+ real(kind=8) :: vectube(3),enetube(nres*2)
+ real(kind=8) :: Etube,xtemp,xminact,yminact,&
+ ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi, &
+ sc_aa_tube,sc_bb_tube
+ integer :: i,j,iti
+ Etube=0.0d0
+ do i=itube_start,itube_end
+ enetube(i)=0.0d0
+ enetube(i+nres)=0.0d0
+ enddo
+!C first we calculate the distance from tube center
+!C for UNRES
+ do i=itube_start,itube_end
+!C lets ommit dummy atoms for now
+ if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle
+!C now calculate distance from center of tube and direction vectors
+ xmin=boxxsize
+ ymin=boxysize
+! Find minimum distance in periodic box
+ do j=-1,1
+ vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+ vectube(1)=vectube(1)+boxxsize*j
+ vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+ vectube(2)=vectube(2)+boxysize*j
+ xminact=abs(vectube(1)-tubecenter(1))
+ yminact=abs(vectube(2)-tubecenter(2))
+ if (xmin.gt.xminact) then
+ xmin=xminact
+ xtemp=vectube(1)
+ endif
+ if (ymin.gt.yminact) then
+ ymin=yminact
+ ytemp=vectube(2)
+ endif
+ enddo
+ vectube(1)=xtemp
+ vectube(2)=ytemp
+ vectube(1)=vectube(1)-tubecenter(1)
+ vectube(2)=vectube(2)-tubecenter(2)
+
+!C print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1)
+!C print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2)
+
+!C as the tube is infinity we do not calculate the Z-vector use of Z
+!C as chosen axis
+ vectube(3)=0.0d0
+!C now calculte the distance
+ tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+!C now normalize vector
+ vectube(1)=vectube(1)/tub_r
+ vectube(2)=vectube(2)/tub_r
+!C calculte rdiffrence between r and r0
+ rdiff=tub_r-tubeR0
+!C and its 6 power
+ rdiff6=rdiff**6.0d0
+!C for vectorization reasons we will sumup at the end to avoid depenence of previous
+ enetube(i)=pep_aa_tube/rdiff6**2.0d0+pep_bb_tube/rdiff6
+!C write(iout,*) "TU13",i,rdiff6,enetube(i)
+!C print *,rdiff,rdiff6,pep_aa_tube
+!C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
+!C now we calculate gradient
+ fac=(-12.0d0*pep_aa_tube/rdiff6- &
+ 6.0d0*pep_bb_tube)/rdiff6/rdiff
+!C write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i),
+!C &rdiff,fac
+!C now direction of gg_tube vector
+ do j=1,3
+ gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0
+ gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0
+ enddo
+ enddo
+!C basically thats all code now we split for side-chains (REMEMBER to sum up at the END)
+!C print *,gg_tube(1,0),"TU"
+
+
+ do i=itube_start,itube_end
+!C Lets not jump over memory as we use many times iti
+ iti=itype(i)
+!C lets ommit dummy atoms for now
+ if ((iti.eq.ntyp1) &
+!C in UNRES uncomment the line below as GLY has no side-chain...
+!C .or.(iti.eq.10)
+ ) cycle
+ xmin=boxxsize
+ ymin=boxysize
+ do j=-1,1
+ vectube(1)=mod((c(1,i+nres)),boxxsize)
+ vectube(1)=vectube(1)+boxxsize*j
+ vectube(2)=mod((c(2,i+nres)),boxysize)
+ vectube(2)=vectube(2)+boxysize*j
+
+ xminact=abs(vectube(1)-tubecenter(1))
+ yminact=abs(vectube(2)-tubecenter(2))
+ if (xmin.gt.xminact) then
+ xmin=xminact
+ xtemp=vectube(1)
+ endif
+ if (ymin.gt.yminact) then
+ ymin=yminact
+ ytemp=vectube(2)
+ endif
+ enddo
+ vectube(1)=xtemp
+ vectube(2)=ytemp
+!C write(iout,*), "tututu", vectube(1),tubecenter(1),vectube(2),
+!C & tubecenter(2)
+ vectube(1)=vectube(1)-tubecenter(1)
+ vectube(2)=vectube(2)-tubecenter(2)
+
+!C as the tube is infinity we do not calculate the Z-vector use of Z
+!C as chosen axis
+ vectube(3)=0.0d0
+!C now calculte the distance
+ tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+!C now normalize vector
+ vectube(1)=vectube(1)/tub_r
+ vectube(2)=vectube(2)/tub_r
+
+!C calculte rdiffrence between r and r0
+ rdiff=tub_r-tubeR0
+!C and its 6 power
+ rdiff6=rdiff**6.0d0
+!C for vectorization reasons we will sumup at the end to avoid depenence of previous
+ sc_aa_tube=sc_aa_tube_par(iti)
+ sc_bb_tube=sc_bb_tube_par(iti)
+ enetube(i+nres)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6
+ fac=-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff- &
+ 6.0d0*sc_bb_tube/rdiff6/rdiff
+!C now direction of gg_tube vector
+ do j=1,3
+ gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac
+ gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
+ enddo
+ enddo
+ do i=itube_start,itube_end
+ Etube=Etube+enetube(i)+enetube(i+nres)
+ enddo
+!C print *,"ETUBE", etube
+ return
+ end subroutine calctube
+!C TO DO 1) add to total energy
+!C 2) add to gradient summation
+!C 3) add reading parameters (AND of course oppening of PARAM file)
+!C 4) add reading the center of tube
+!C 5) add COMMONs
+!C 6) add to zerograd
+!C 7) allocate matrices
+
+
+!C-----------------------------------------------------------------------
+!C-----------------------------------------------------------
+!C This subroutine is to mimic the histone like structure but as well can be
+!C utilizet to nanostructures (infinit) small modification has to be used to
+!C make it finite (z gradient at the ends has to be changes as well as the x,y
+!C gradient has to be modified at the ends
+!C The energy function is Kihara potential
+!C E=4esp*((sigma/(r-r0))^12 - (sigma/(r-r0))^6)
+!C 4eps is depth of well sigma is r_minimum r is distance from center of tube
+!C and r0 is the excluded size of nanotube (can be set to 0 if we want just a
+!C simple Kihara potential
+ subroutine calctube2(Etube)
+ real(kind=8) :: vectube(3),enetube(nres*2)
+ real(kind=8) :: Etube,xtemp,xminact,yminact,&
+ ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi,fracinbuf,&
+ sstube,ssgradtube,sc_aa_tube,sc_bb_tube
+ integer:: i,j,iti
+ Etube=0.0d0
+ do i=itube_start,itube_end
+ enetube(i)=0.0d0
+ enetube(i+nres)=0.0d0
+ enddo
+!C first we calculate the distance from tube center
+!C first sugare-phosphate group for NARES this would be peptide group
+!C for UNRES
+ do i=itube_start,itube_end
+!C lets ommit dummy atoms for now
+
+ if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle
+!C now calculate distance from center of tube and direction vectors
+!C vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+!C if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
+!C vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+!C if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize
+ xmin=boxxsize
+ ymin=boxysize
+ do j=-1,1
+ vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+ vectube(1)=vectube(1)+boxxsize*j
+ vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+ vectube(2)=vectube(2)+boxysize*j
+
+ xminact=abs(vectube(1)-tubecenter(1))
+ yminact=abs(vectube(2)-tubecenter(2))
+ if (xmin.gt.xminact) then
+ xmin=xminact
+ xtemp=vectube(1)
+ endif
+ if (ymin.gt.yminact) then
+ ymin=yminact
+ ytemp=vectube(2)
+ endif
+ enddo
+ vectube(1)=xtemp
+ vectube(2)=ytemp
+ vectube(1)=vectube(1)-tubecenter(1)
+ vectube(2)=vectube(2)-tubecenter(2)
+
+!C print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1)
+!C print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2)
+
+!C as the tube is infinity we do not calculate the Z-vector use of Z
+!C as chosen axis
+ vectube(3)=0.0d0
+!C now calculte the distance
+ tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+!C now normalize vector
+ vectube(1)=vectube(1)/tub_r
+ vectube(2)=vectube(2)/tub_r
+!C calculte rdiffrence between r and r0
+ rdiff=tub_r-tubeR0
+!C and its 6 power
+ rdiff6=rdiff**6.0d0
+!C THIS FRAGMENT MAKES TUBE FINITE
+ positi=mod((c(3,i)+c(3,i+1))/2.0d0,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,bordtubebot,buftubebot,bordtubetop
+ if ((positi.gt.bordtubebot) &
+ .and.(positi.lt.bordtubetop)) then
+!C the energy transfer exist
+ if (positi.lt.buftubebot) then
+ fracinbuf=1.0d0- &
+ ((positi-bordtubebot)/tubebufthick)
+!C lipbufthick is thickenes of lipid buffore
+ sstube=sscalelip(fracinbuf)
+ ssgradtube=-sscagradlip(fracinbuf)/tubebufthick
+!C print *,ssgradtube, sstube,tubetranene(itype(i))
+ enetube(i)=enetube(i)+sstube*tubetranenepep
+!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
+!C &+ssgradtube*tubetranene(itype(i))
+!C gg_tube(3,i-1)= gg_tube(3,i-1)
+!C &+ssgradtube*tubetranene(itype(i))
+!C print *,"doing sccale for lower part"
+ elseif (positi.gt.buftubetop) then
+ fracinbuf=1.0d0- &
+ ((bordtubetop-positi)/tubebufthick)
+ sstube=sscalelip(fracinbuf)
+ ssgradtube=sscagradlip(fracinbuf)/tubebufthick
+ enetube(i)=enetube(i)+sstube*tubetranenepep
+!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
+!C &+ssgradtube*tubetranene(itype(i))
+!C gg_tube(3,i-1)= gg_tube(3,i-1)
+!C &+ssgradtube*tubetranene(itype(i))
+!C print *, "doing sscalefor top part",sslip,fracinbuf
+ else
+ sstube=1.0d0
+ ssgradtube=0.0d0
+ enetube(i)=enetube(i)+sstube*tubetranenepep
+!C print *,"I am in true lipid"
+ endif
+ else
+!C sstube=0.0d0
+!C ssgradtube=0.0d0
+ cycle
+ endif ! if in lipid or buffor
+
+!C for vectorization reasons we will sumup at the end to avoid depenence of previous
+ enetube(i)=enetube(i)+sstube* &
+ (pep_aa_tube/rdiff6**2.0d0+pep_bb_tube/rdiff6)
+!C write(iout,*) "TU13",i,rdiff6,enetube(i)
+!C print *,rdiff,rdiff6,pep_aa_tube
+!C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
+!C now we calculate gradient
+ fac=(-12.0d0*pep_aa_tube/rdiff6- &
+ 6.0d0*pep_bb_tube)/rdiff6/rdiff*sstube
+!C write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i),
+!C &rdiff,fac
+
+!C now direction of gg_tube vector
+ do j=1,3
+ gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0
+ gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0
+ enddo
+ gg_tube(3,i)=gg_tube(3,i) &
+ +ssgradtube*enetube(i)/sstube/2.0d0
+ gg_tube(3,i-1)= gg_tube(3,i-1) &
+ +ssgradtube*enetube(i)/sstube/2.0d0
+
+ enddo
+!C basically thats all code now we split for side-chains (REMEMBER to sum up at the END)
+!C print *,gg_tube(1,0),"TU"
+ do i=itube_start,itube_end
+!C Lets not jump over memory as we use many times iti
+ iti=itype(i)
+!C lets ommit dummy atoms for now
+ if ((iti.eq.ntyp1) &
+!!C in UNRES uncomment the line below as GLY has no side-chain...
+ .or.(iti.eq.10) &
+ ) cycle
+ vectube(1)=c(1,i+nres)
+ vectube(1)=mod(vectube(1),boxxsize)
+ if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
+ vectube(2)=c(2,i+nres)
+ vectube(2)=mod(vectube(2),boxysize)
+ if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize
+
+ vectube(1)=vectube(1)-tubecenter(1)
+ vectube(2)=vectube(2)-tubecenter(2)
+!C THIS FRAGMENT MAKES TUBE FINITE
+ 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,bordtubebot,buftubebot,bordtubetop
+
+ if ((positi.gt.bordtubebot) &
+ .and.(positi.lt.bordtubetop)) then
+!C the energy transfer exist
+ if (positi.lt.buftubebot) then
+ fracinbuf=1.0d0- &
+ ((positi-bordtubebot)/tubebufthick)
+!C lipbufthick is thickenes of lipid buffore
+ sstube=sscalelip(fracinbuf)
+ ssgradtube=-sscagradlip(fracinbuf)/tubebufthick
+!C print *,ssgradtube, sstube,tubetranene(itype(i))
+ enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i))
+!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
+!C &+ssgradtube*tubetranene(itype(i))
+!C gg_tube(3,i-1)= gg_tube(3,i-1)
+!C &+ssgradtube*tubetranene(itype(i))
+!C print *,"doing sccale for lower part"
+ elseif (positi.gt.buftubetop) then
+ fracinbuf=1.0d0- &
+ ((bordtubetop-positi)/tubebufthick)
+
+ sstube=sscalelip(fracinbuf)
+ ssgradtube=sscagradlip(fracinbuf)/tubebufthick
+ enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i))
+!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
+!C &+ssgradtube*tubetranene(itype(i))
+!C gg_tube(3,i-1)= gg_tube(3,i-1)
+!C &+ssgradtube*tubetranene(itype(i))
+!C print *, "doing sscalefor top part",sslip,fracinbuf
+ else
+ sstube=1.0d0
+ ssgradtube=0.0d0
+ enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i))
+!C print *,"I am in true lipid"
+ endif
+ else
+!C sstube=0.0d0
+!C ssgradtube=0.0d0
+ cycle
+ endif ! if in lipid or buffor
+!CEND OF FINITE FRAGMENT
+!C as the tube is infinity we do not calculate the Z-vector use of Z
+!C as chosen axis
+ vectube(3)=0.0d0
+!C now calculte the distance
+ tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+!C now normalize vector
+ vectube(1)=vectube(1)/tub_r
+ vectube(2)=vectube(2)/tub_r
+!C calculte rdiffrence between r and r0
+ rdiff=tub_r-tubeR0
+!C and its 6 power
+ rdiff6=rdiff**6.0d0
+!C for vectorization reasons we will sumup at the end to avoid depenence of previous
+ sc_aa_tube=sc_aa_tube_par(iti)
+ sc_bb_tube=sc_bb_tube_par(iti)
+ enetube(i+nres)=(sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6)&
+ *sstube+enetube(i+nres)
+!C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
+!C now we calculate gradient
+ fac=(-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff-&
+ 6.0d0*sc_bb_tube/rdiff6/rdiff)*sstube
+!C now direction of gg_tube vector
+ do j=1,3
+ gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac
+ gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
+ enddo
+ gg_tube_SC(3,i)=gg_tube_SC(3,i) &
+ +ssgradtube*enetube(i+nres)/sstube
+ gg_tube(3,i-1)= gg_tube(3,i-1) &
+ +ssgradtube*enetube(i+nres)/sstube
+
+ enddo
+ do i=itube_start,itube_end
+ Etube=Etube+enetube(i)+enetube(i+nres)
+ enddo
+!C print *,"ETUBE", etube
+ return
+ end subroutine calctube2
+!=====================================================================================================================================
+ subroutine calcnano(Etube)
+ real(kind=8) :: vectube(3),enetube(nres*2), &
+ enecavtube(nres*2)
+ real(kind=8) :: Etube,xtemp,xminact,yminact,&
+ ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,denominator,faccav,&
+ sc_aa_tube,sc_bb_tube,zmin,ztemp,zminact
+ integer:: i,j,iti
+
+ Etube=0.0d0
+! print *,itube_start,itube_end,"poczatek"
+ do i=itube_start,itube_end
+ enetube(i)=0.0d0
+ enetube(i+nres)=0.0d0
+ enddo
+!C first we calculate the distance from tube center
+!C first sugare-phosphate group for NARES this would be peptide group
+!C for UNRES
+ do i=itube_start,itube_end
+!C lets ommit dummy atoms for now
+ if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle
+!C now calculate distance from center of tube and direction vectors
+ xmin=boxxsize
+ ymin=boxysize
+ zmin=boxzsize
+
+ do j=-1,1
+ vectube(1)=dmod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+ vectube(1)=vectube(1)+boxxsize*j
+ vectube(2)=dmod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+ vectube(2)=vectube(2)+boxysize*j
+ vectube(3)=dmod((c(3,i)+c(3,i+1))/2.0d0,boxzsize)
+ vectube(3)=vectube(3)+boxzsize*j
+
+
+ xminact=dabs(vectube(1)-tubecenter(1))
+ yminact=dabs(vectube(2)-tubecenter(2))
+ zminact=dabs(vectube(3)-tubecenter(3))
+
+ if (xmin.gt.xminact) then
+ xmin=xminact
+ xtemp=vectube(1)
+ endif
+ if (ymin.gt.yminact) then
+ ymin=yminact
+ ytemp=vectube(2)
+ endif
+ if (zmin.gt.zminact) then
+ zmin=zminact
+ ztemp=vectube(3)
+ endif
+ enddo
+ vectube(1)=xtemp
+ vectube(2)=ytemp
+ vectube(3)=ztemp
+
+ vectube(1)=vectube(1)-tubecenter(1)
+ vectube(2)=vectube(2)-tubecenter(2)
+ vectube(3)=vectube(3)-tubecenter(3)
+
+!C print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1)
+!C print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2)
+!C as the tube is infinity we do not calculate the Z-vector use of Z
+!C as chosen axis
+!C vectube(3)=0.0d0
+!C now calculte the distance
+ tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+!C now normalize vector
+ vectube(1)=vectube(1)/tub_r
+ vectube(2)=vectube(2)/tub_r
+ vectube(3)=vectube(3)/tub_r
+!C calculte rdiffrence between r and r0
+ rdiff=tub_r-tubeR0
+!C and its 6 power
+ rdiff6=rdiff**6.0d0
+!C for vectorization reasons we will sumup at the end to avoid depenence of previous
+ enetube(i)=pep_aa_tube/rdiff6**2.0d0+pep_bb_tube/rdiff6
+!C write(iout,*) "TU13",i,rdiff6,enetube(i)
+!C print *,rdiff,rdiff6,pep_aa_tube
+!C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
+!C now we calculate gradient
+ fac=(-12.0d0*pep_aa_tube/rdiff6- &
+ 6.0d0*pep_bb_tube)/rdiff6/rdiff
+!C write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i),
+!C &rdiff,fac
+ if (acavtubpep.eq.0.0d0) then
+!C go to 667
+ enecavtube(i)=0.0
+ faccav=0.0
+ else
+ denominator=(1.0d0+dcavtubpep*rdiff6*rdiff6)
+ enecavtube(i)= &
+ (bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff)+ccavtubpep) &
+ /denominator
+ enecavtube(i)=0.0
+ faccav=((bcavtubpep*1.0d0+acavtubpep/2.0d0/dsqrt(rdiff)) &
+ *denominator-(bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff) &
+ +ccavtubpep)*rdiff6**2.0d0/rdiff*dcavtubpep*12.0d0) &
+ /denominator**2.0d0
+!C faccav=0.0
+!C fac=fac+faccav
+!C 667 continue
+ endif
+
+ do j=1,3
+ gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0
+ gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0
+ enddo
+ enddo
+
+ do i=itube_start,itube_end
+ enecavtube(i)=0.0d0
+!C Lets not jump over memory as we use many times iti
+ iti=itype(i)
+!C lets ommit dummy atoms for now
+ if ((iti.eq.ntyp1) &
+!C in UNRES uncomment the line below as GLY has no side-chain...
+!C .or.(iti.eq.10)
+ ) cycle
+ xmin=boxxsize
+ ymin=boxysize
+ zmin=boxzsize
+ do j=-1,1
+ vectube(1)=dmod((c(1,i+nres)),boxxsize)
+ vectube(1)=vectube(1)+boxxsize*j
+ vectube(2)=dmod((c(2,i+nres)),boxysize)
+ vectube(2)=vectube(2)+boxysize*j
+ vectube(3)=dmod((c(3,i+nres)),boxzsize)
+ vectube(3)=vectube(3)+boxzsize*j
+
+
+ xminact=dabs(vectube(1)-tubecenter(1))
+ yminact=dabs(vectube(2)-tubecenter(2))
+ zminact=dabs(vectube(3)-tubecenter(3))
+
+ if (xmin.gt.xminact) then
+ xmin=xminact
+ xtemp=vectube(1)
+ endif
+ if (ymin.gt.yminact) then
+ ymin=yminact
+ ytemp=vectube(2)
+ endif
+ if (zmin.gt.zminact) then
+ zmin=zminact
+ ztemp=vectube(3)
+ endif
+ enddo
+ vectube(1)=xtemp
+ vectube(2)=ytemp
+ vectube(3)=ztemp
+
+!C write(iout,*), "tututu", vectube(1),tubecenter(1),vectube(2),
+!C & tubecenter(2)
+ vectube(1)=vectube(1)-tubecenter(1)
+ vectube(2)=vectube(2)-tubecenter(2)
+ vectube(3)=vectube(3)-tubecenter(3)
+!C now calculte the distance
+ tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+!C now normalize vector
+ vectube(1)=vectube(1)/tub_r
+ vectube(2)=vectube(2)/tub_r
+ vectube(3)=vectube(3)/tub_r
+
+!C calculte rdiffrence between r and r0
+ rdiff=tub_r-tubeR0
+!C and its 6 power
+ rdiff6=rdiff**6.0d0
+ sc_aa_tube=sc_aa_tube_par(iti)
+ sc_bb_tube=sc_bb_tube_par(iti)
+ enetube(i+nres)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6
+!C enetube(i+nres)=0.0d0
+!C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
+!C now we calculate gradient
+ fac=-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff- &
+ 6.0d0*sc_bb_tube/rdiff6/rdiff
+!C fac=0.0
+!C now direction of gg_tube vector
+!C Now cavity term E=a(x+bsqrt(x)+c)/(1+dx^12)
+ if (acavtub(iti).eq.0.0d0) then
+!C go to 667
+ enecavtube(i+nres)=0.0d0
+ faccav=0.0d0
+ else
+ denominator=(1.0d0+dcavtub(iti)*rdiff6*rdiff6)
+ enecavtube(i+nres)= &
+ (bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(rdiff)+ccavtub(iti)) &
+ /denominator
+!C enecavtube(i)=0.0
+ faccav=((bcavtub(iti)*1.0d0+acavtub(iti)/2.0d0/dsqrt(rdiff)) &
+ *denominator-(bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(rdiff) &
+ +ccavtub(iti))*rdiff6**2.0d0/rdiff*dcavtub(iti)*12.0d0) &
+ /denominator**2.0d0
+!C faccav=0.0
+ fac=fac+faccav
+!C 667 continue
+ endif
+!C print *,"TUT",i,iti,rdiff,rdiff6,acavtub(iti),denominator,
+!C & enecavtube(i),faccav
+!C print *,"licz=",
+!C & (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti))
+!C print *,"finene=",enetube(i+nres)+enecavtube(i)
+ do j=1,3
+ gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac
+ gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
+ enddo
+ enddo
+
+
+
+ do i=itube_start,itube_end
+ Etube=Etube+enetube(i)+enetube(i+nres)+enecavtube(i) &
+ +enecavtube(i+nres)
+ enddo
+!C print *,"ETUBE", etube
+ return
+ end subroutine calcnano
+
+!===============================================
+!--------------------------------------------------------------------------------
+!C first for shielding is setting of function of side-chains
+
+ subroutine set_shield_fac2
+ real(kind=8) :: div77_81=0.974996043d0, &
+ div4_81=0.2222222222d0
+ real (kind=8) :: dist_pep_side,dist_side_calf,dist_pept_group, &
+ scale_fac_dist,fac_help_scale,VofOverlap,VolumeTotal,costhet,&
+ short,long,sinthet,costhet_fac,sh_frac_dist,rkprim,cosphi, &
+ sinphi,cosphi_fac,pep_side0pept_group,cosalfa,fac_alfa_sin
+!C the vector between center of side_chain and peptide group
+ real(kind=8),dimension(3) :: pep_side_long,side_calf, &
+ pept_group,costhet_grad,cosphi_grad_long, &
+ cosphi_grad_loc,pep_side_norm,side_calf_norm, &
+ sh_frac_dist_grad,pep_side
+ integer i,j,k
+!C write(2,*) "ivec",ivec_start,ivec_end
+ do i=1,nres
+ fac_shield(i)=0.0d0
+ do j=1,3
+ grad_shield(j,i)=0.0d0
+ enddo
+ enddo
+ do i=ivec_start,ivec_end
+!C do i=1,nres-1
+!C if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+ ishield_list(i)=0
+ if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+!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=sqrt(dist_pep_side)
+ dist_pept_group=sqrt(dist_pept_group)
+ dist_side_calf=sqrt(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 print *,ishield_list(i),i
+!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
+ 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,itype(i),fac_shield(i)
+ enddo
+ return
+ end subroutine set_shield_fac2
+!----------------------------------------------------------------------------
+! SOUBROUTINE FOR AFM
+ subroutine AFMvel(Eafmforce)
+ use MD_data, only:totTafm
+ real(kind=8),dimension(3) :: diffafm
+ real(kind=8) :: afmdist,Eafmforce
+ integer :: i
+!C Only for check grad COMMENT if not used for checkgrad
+!C totT=3.0d0
+!C--------------------------------------------------------
+!C print *,"wchodze"
+ afmdist=0.0d0
+ Eafmforce=0.0d0
+ do i=1,3
+ diffafm(i)=c(i,afmend)-c(i,afmbeg)
+ afmdist=afmdist+diffafm(i)**2
+ enddo
+ afmdist=dsqrt(afmdist)
+! totTafm=3.0
+ Eafmforce=0.5d0*forceAFMconst &
+ *(distafminit+totTafm*velAFMconst-afmdist)**2
+!C Eafmforce=-forceAFMconst*(dist-distafminit)
+ do i=1,3
+ gradafm(i,afmend-1)=-forceAFMconst* &
+ (distafminit+totTafm*velAFMconst-afmdist) &
+ *diffafm(i)/afmdist
+ gradafm(i,afmbeg-1)=forceAFMconst* &
+ (distafminit+totTafm*velAFMconst-afmdist) &
+ *diffafm(i)/afmdist
+ enddo
+! print *,'AFM',Eafmforce,totTafm*velAFMconst,afmdist
+ return
+ end subroutine AFMvel
+!---------------------------------------------------------
+ subroutine AFMforce(Eafmforce)
+
+ real(kind=8),dimension(3) :: diffafm
+! real(kind=8) ::afmdist
+ real(kind=8) :: afmdist,Eafmforce
+ integer :: i
+ afmdist=0.0d0
+ Eafmforce=0.0d0
+ do i=1,3
+ diffafm(i)=c(i,afmend)-c(i,afmbeg)
+ afmdist=afmdist+diffafm(i)**2
+ enddo
+ afmdist=dsqrt(afmdist)
+! print *,afmdist,distafminit
+ Eafmforce=-forceAFMconst*(afmdist-distafminit)
+ do i=1,3
+ gradafm(i,afmend-1)=-forceAFMconst*diffafm(i)/afmdist
+ gradafm(i,afmbeg-1)=forceAFMconst*diffafm(i)/afmdist
+ enddo
+!C print *,'AFM',Eafmforce
+ return
+ end subroutine AFMforce
+
!-----------------------------------------------------------------------------
#ifdef WHAM
subroutine read_ssHist
allocate(grij_hb_cont(3,maxconts,nres))
!(3,maxconts,maxres)
allocate(facont_hb(maxconts,nres))
+
allocate(ees0p(maxconts,nres))
allocate(ees0m(maxconts,nres))
allocate(d_cont(maxconts,nres))
+ allocate(ees0plist(maxconts,nres))
+
!(maxconts,maxres)
allocate(num_cont_hb(nres))
!(maxres)
!(6,maxdim)
allocate(dxds(6,nres))
!(6,maxres)
- allocate(gradx(3,nres,0:2))
- allocate(gradc(3,nres,0:2))
+ allocate(gradx(3,-1:nres,0:2))
+ allocate(gradc(3,-1:nres,0:2))
!(3,maxres,2)
- allocate(gvdwx(3,nres))
- allocate(gvdwc(3,nres))
- allocate(gelc(3,nres))
- allocate(gelc_long(3,nres))
- allocate(gvdwpp(3,nres))
- allocate(gvdwc_scpp(3,nres))
- allocate(gradx_scp(3,nres))
- allocate(gvdwc_scp(3,nres))
- allocate(ghpbx(3,nres))
- allocate(ghpbc(3,nres))
- allocate(gradcorr(3,nres))
- allocate(gradcorr_long(3,nres))
- allocate(gradcorr5_long(3,nres))
- allocate(gradcorr6_long(3,nres))
- allocate(gcorr6_turn_long(3,nres))
- allocate(gradxorr(3,nres))
- allocate(gradcorr5(3,nres))
- allocate(gradcorr6(3,nres))
+ allocate(gvdwx(3,-1:nres))
+ allocate(gvdwc(3,-1:nres))
+ allocate(gelc(3,-1:nres))
+ allocate(gelc_long(3,-1:nres))
+ allocate(gvdwpp(3,-1:nres))
+ allocate(gvdwc_scpp(3,-1:nres))
+ allocate(gradx_scp(3,-1:nres))
+ allocate(gvdwc_scp(3,-1:nres))
+ allocate(ghpbx(3,-1:nres))
+ allocate(ghpbc(3,-1:nres))
+ allocate(gradcorr(3,-1:nres))
+ allocate(gradcorr_long(3,-1:nres))
+ allocate(gradcorr5_long(3,-1:nres))
+ allocate(gradcorr6_long(3,-1:nres))
+ allocate(gcorr6_turn_long(3,-1:nres))
+ allocate(gradxorr(3,-1:nres))
+ allocate(gradcorr5(3,-1:nres))
+ allocate(gradcorr6(3,-1:nres))
+ allocate(gliptran(3,-1:nres))
+ allocate(gliptranc(3,-1:nres))
+ allocate(gliptranx(3,-1:nres))
+ allocate(gshieldx(3,-1:nres))
+ allocate(gshieldc(3,-1:nres))
+ allocate(gshieldc_loc(3,-1:nres))
+ allocate(gshieldx_ec(3,-1:nres))
+ allocate(gshieldc_ec(3,-1:nres))
+ allocate(gshieldc_loc_ec(3,-1:nres))
+ allocate(gshieldx_t3(3,-1:nres))
+ allocate(gshieldc_t3(3,-1:nres))
+ allocate(gshieldc_loc_t3(3,-1:nres))
+ allocate(gshieldx_t4(3,-1:nres))
+ allocate(gshieldc_t4(3,-1:nres))
+ allocate(gshieldc_loc_t4(3,-1:nres))
+ allocate(gshieldx_ll(3,-1:nres))
+ allocate(gshieldc_ll(3,-1:nres))
+ allocate(gshieldc_loc_ll(3,-1:nres))
+ allocate(grad_shield(3,-1:nres))
+ allocate(gg_tube_sc(3,-1:nres))
+ allocate(gg_tube(3,-1:nres))
+ allocate(gradafm(3,-1:nres))
!(3,maxres)
+ allocate(grad_shield_side(3,50,nres))
+ allocate(grad_shield_loc(3,50,nres))
+! grad for shielding surroing
allocate(gloc(0:maxvar,0:2))
allocate(gloc_x(0:maxvar,2))
!(maxvar,2)
- allocate(gel_loc(3,nres))
- allocate(gel_loc_long(3,nres))
- allocate(gcorr3_turn(3,nres))
- allocate(gcorr4_turn(3,nres))
- allocate(gcorr6_turn(3,nres))
- allocate(gradb(3,nres))
- allocate(gradbx(3,nres))
+ allocate(gel_loc(3,-1:nres))
+ allocate(gel_loc_long(3,-1:nres))
+ allocate(gcorr3_turn(3,-1:nres))
+ allocate(gcorr4_turn(3,-1:nres))
+ allocate(gcorr6_turn(3,-1:nres))
+ allocate(gradb(3,-1:nres))
+ allocate(gradbx(3,-1:nres))
!(3,maxres)
allocate(gel_loc_loc(maxvar))
allocate(gel_loc_turn3(maxvar))
allocate(g_corr5_loc(maxvar))
allocate(g_corr6_loc(maxvar))
!(maxvar)
- allocate(gsccorc(3,nres))
- allocate(gsccorx(3,nres))
+ allocate(gsccorc(3,-1:nres))
+ allocate(gsccorx(3,-1:nres))
!(3,maxres)
- allocate(gsccor_loc(nres))
+ allocate(gsccor_loc(-1:nres))
!(maxres)
- allocate(dtheta(3,2,nres))
+ allocate(dtheta(3,2,-1:nres))
!(3,2,maxres)
- allocate(gscloc(3,nres))
- allocate(gsclocx(3,nres))
+ allocate(gscloc(3,-1:nres))
+ allocate(gsclocx(3,-1:nres))
!(3,maxres)
- allocate(dphi(3,3,nres))
- allocate(dalpha(3,3,nres))
- allocate(domega(3,3,nres))
+ allocate(dphi(3,3,-1:nres))
+ allocate(dalpha(3,3,-1:nres))
+ allocate(domega(3,3,-1:nres))
!(3,3,maxres)
! common /deriv_scloc/
allocate(dXX_C1tab(3,nres))
!----------------------
! common.MD
! common /mdgrad/
- allocate(gcart(3,0:nres))
- allocate(gxcart(3,0:nres))
+ allocate(gcart(3,-1:nres))
+ allocate(gxcart(3,-1:nres))
!(3,0:MAXRES)
- allocate(gradcag(3,nres))
- allocate(gradxag(3,nres))
+ allocate(gradcag(3,-1:nres))
+ allocate(gradxag(3,-1:nres))
!(3,MAXRES)
! common /back_constr/
!el in energy:Econstr_back allocate((:),allocatable :: utheta,ugamma,uscdiff !(maxfrag_back)
! enddo
! enddo
- if (nss.gt.0) then
- allocate(idssb(nss),jdssb(nss))
+! if (nss.gt.0) then
+ allocate(idssb(maxdim),jdssb(maxdim))
+! allocate(newihpb(nss),newjhpb(nss))
!(maxdim)
- endif
+! endif
+ allocate(ishield_list(nres))
+ allocate(shield_list(50,nres))
allocate(dyn_ss_mask(nres))
+ allocate(fac_shield(nres))
!(maxres)
dyn_ss_mask(:)=.false.
!----------------------
projects / unres4.git / blobdiff