real(kind=8),dimension(:,:,:),allocatable :: gacont !(3,maxconts,maxres)
integer,dimension(:),allocatable :: ishield_list
integer,dimension(:,:),allocatable :: shield_list
+ real(kind=8),dimension(:),allocatable :: enetube,enecavtube
!
! 12/26/95 - H-bonding contacts
! common /contacts_hb/
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 !(3,maxres)
+ grad_shield,gg_tube,gg_tube_sc,gradafm !(3,maxres)
+!-----------------------------NUCLEIC GRADIENT
+ real(kind=8),dimension(:,:),allocatable ::gradb_nucl,gradbx_nucl, &
+ gvdwpsb1,gelpp,gvdwpsb,gelsbc,gelsbx,gvdwsbx,gvdwsbc,gsbloc,&
+ gsblocx,gradcorr_nucl,gradxorr_nucl,gradcorr3_nucl,gradxorr3_nucl
! 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 :: uygrad,uzgrad !(3,3,2,maxres)
!-----------------------------------------------------------------------------
! common /przechowalnia/
- real(kind=8),dimension(:,:,:),allocatable :: zapas !(max_dim,maxconts,max_fg_procs)
+ real(kind=8),dimension(:,:,:),allocatable :: zapas
+ real(kind=8),dimension(:,:,:,:),allocatable ::zapas2 !(max_dim,maxconts,max_fg_procs)
real(kind=8),dimension(:,:,:),allocatable :: fromto !(3,3,maxdim)(maxdim=(maxres-1)*(maxres-2)/2)
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
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,eliptran,etube
+ real(kind=8) :: eello_turn3,eello_turn4,estr,ebe,eliptran,etube, &
+ Eafmforce,ethetacnstr
real(kind=8) :: ecorr,ecorr5,ecorr6,eturn6
-
+! now energies for nulceic alone parameters
+ real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,&
+ ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
+ ecorr3_nucl
#ifdef MPI
real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw
! shielding effect varibles for MPI
#endif
!
! Compute the side-chain and electrostatic interaction energy
- print *, "Before EVDW"
+! print *, "Before EVDW"
! goto (101,102,103,104,105,106) ipot
select case(ipot)
! Lennard-Jones potential.
if (shield_mode.eq.2) then
call set_shield_fac2
endif
+ print *,"AFTER EGB",ipot,evdw
!mc
!mc Sep-06: egb takes care of dynamic ss bonds too
!mc
! print *,"Processor",myrank," left VEC_AND_DERIV"
if (ipot.lt.6) then
#ifdef SPLITELE
- print *,"after ipot if", ipot
+! 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 &
! Calculate the bond-stretching energy
!
call ebond(estr)
+ print *,"EBOND",estr
! write(iout,*) "in etotal afer ebond",ipot
!
! Calculate the virtual-bond-angle energy.
!
if (wang.gt.0d0) then
- call ebend(ebe)
+ call ebend(ebe,ethetacnstr)
else
ebe=0
endif
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
else
etube=0.0d0
endif
-
+!--------------------------------------------------------
+ print *,"before",ees,evdw1,ecorr
+ call ebond_nucl(estr_nucl)
+ call ebend_nucl(ebe_nucl)
+ call etor_nucl(etors_nucl)
+ call esb_gb(evdwsb,eelsb)
+ call epp_nucl_sub(evdwpp,eespp)
+ call epsb(evdwpsb,eelpsb)
+ call esb(esbloc)
+ call multibody_hb_nucl(ecorr_nucl,ecorr3_nucl,n_corr,n_corr1)
+
+ print *,"after ebend", ebe_nucl
#ifdef TIMING
time_enecalc=time_enecalc+MPI_Wtime()-time00
#endif
energia(20)=Uconst+Uconst_back
energia(21)=esccor
energia(22)=eliptran
+ energia(23)=Eafmforce
+ energia(24)=ethetacnstr
energia(25)=etube
+!---------------------------------------------------------------
+ energia(26)=evdwpp
+ energia(27)=eespp
+ energia(28)=evdwpsb
+ energia(29)=eelpsb
+ energia(30)=evdwsb
+ energia(31)=eelsb
+ energia(32)=estr_nucl
+ energia(33)=ebe_nucl
+ energia(34)=esbloc
+ energia(35)=etors_nucl
+ energia(36)=etors_d_nucl
+ energia(37)=ecorr_nucl
+ energia(38)=ecorr3_nucl
+!----------------------------------------------------------------------
! 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"
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, &
- eliptran,etube
+ eliptran,etube, Eafmforce,ethetacnstr
+ real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,&
+ ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
+ ecorr3_nucl
+
integer :: i
#ifdef MPI
integer :: ierr
Uconst=energia(20)
esccor=energia(21)
eliptran=energia(22)
+ Eafmforce=energia(23)
+ ethetacnstr=energia(24)
etube=energia(25)
+ evdwpp=energia(26)
+ eespp=energia(27)
+ evdwpsb=energia(28)
+ eelpsb=energia(29)
+ evdwsb=energia(30)
+ eelsb=energia(31)
+ estr_nucl=energia(32)
+ ebe_nucl=energia(33)
+ esbloc=energia(34)
+ etors_nucl=energia(35)
+ etors_d_nucl=energia(36)
+ ecorr_nucl=energia(37)
+ ecorr3_nucl=energia(38)
+
+
#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+wliptran*eliptran+wtube*etube
+ +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube&
+ +Eafmforce+ethetacnstr &
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl
#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+wliptran*eliptran+wtube*etube
+ +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube&
+ +Eafmforce+ethetacnstr &
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl
#endif
energia(0)=etot
! detecting NaNQ
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,eliptran,&
- etube
+ etube,ethetacnstr,Eafmforce
+ real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,&
+ ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
+ ecorr3_nucl
etot=energia(0)
evdw=energia(1)
Uconst=energia(20)
esccor=energia(21)
eliptran=energia(22)
+ Eafmforce=energia(23)
+ ethetacnstr=energia(24)
etube=energia(25)
+ evdwpp=energia(26)
+ eespp=energia(27)
+ evdwpsb=energia(28)
+ eelpsb=energia(29)
+ evdwsb=energia(30)
+ eelsb=energia(31)
+ estr_nucl=energia(32)
+ ebe_nucl=energia(33)
+ esbloc=energia(34)
+ etors_nucl=energia(35)
+ etors_d_nucl=energia(36)
+ ecorr_nucl=energia(37)
+ ecorr3_nucl=energia(38)
+
#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,eliptran,wliptran,etube,wtube,etot
+ edihcnstr,ethetacnstr,ebr*nss,&
+ Uconst,eliptran,wliptran,Eafmforce,etube,wtube, & ! till now protein
+ estr_nucl,wbond_nucl,ebe_nucl,wang_nucl, &
+ evdwpp,wvdwpp,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb,&
+ evdwsb,wvdwsb,eelsb,welsb,esbloc,wsbloc,etors_nucl,wtor_nucl,&
+ etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,&
+ ecorr3_nucl,wcorr3_nucl, &
+ 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)'/ &
+ 'ESTR_nucl=',1pE16.6,' WEIGHT=',1pD16.6,' (stretching for nucleic)'/ &
+ 'EBE_nucl=',1pE16.6,' WEIGHT=',1pD16.6,' (bending for nucleic)'/ &
+ 'EVDW_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(phosphate-phosphate VDW)'/ &
+ 'EESPP_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(phosphate-phosphate elec)'/ &
+ 'EVDWPSB_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(phosphate-sugarbase VDW)'/ &
+ 'EESPSB_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(phosphate-sugarbase elec)'/ &
+ 'EVDWSB_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(sugarbase-sugarbase VDW)'/ &
+ 'EESSB_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(sugarbase-sugarbase elec)'/ &
+ 'ESBLOC_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(sugarbase rotamer)'/ &
+ 'ETORS_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(torsional)'/ &
+ 'ETORSD_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(double torsional)'/ &
+ 'ECORR_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(multibody 4th order)'/ &
+ 'ECORR3_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(multibody 3th order)'/ &
'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,eliptran,wliptran,etube,wtube,etot
+ ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, &
+ etube,wtube, &
+ estr_nucl,wbond_nucl, ebe_nucl,wang_nucl,&
+ evdwpp,wvdwpp,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb&
+ evdwsb,wvdwsb,eelsb,welsb,esbloc,wsbloc,etors_nucl,wtor_nucl&
+ etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,&
+ ecorr3_nucl,wcorr3_nucl, &
+ 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)'/ &
+ 'ESTR_nucl= ',1pE16.6,' WEIGHT=',1pD16.6,' (stretching for nucleic)'/ &
+ 'EBE_nucl=',1pE16.6,' WEIGHT=',1pD16.6,' (bending for nucleic)'/ &
+ 'EVDW_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(phosphate-phosphate VDW)'/ &
+ 'EESPP_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(phosphate-phosphate elec)'/ &
+ 'EVDWPSB_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(phosphate-sugarbase VDW)'/ &
+ 'EESPSB_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(phosphate-sugarbase elec)'/ &
+ 'EVDWSB_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(sugarbase-sugarbase VDW)'/ &
+ 'EESSB_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(sugarbase-sugarbase elec)'/ &
+ 'ESBLOC_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(sugarbase rotamer)'/ &
+ 'ETORS_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(torsional)'/ &
+ 'ETORSD_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(double torsional)'/ &
+ 'ECORR_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(multibody 4th order)'/ &
+ 'ECORR3_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(multibody 3th order)'/ &
'ETOT= ',1pE16.6,' (total)')
#endif
return
! allocate(gacont(3,nres/4,iatsc_s:iatsc_e)) !(3,maxconts,maxres)
do i=iatsc_s,iatsc_e
- itypi=iabs(itype(i))
+ itypi=iabs(itype(i,1))
if (itypi.eq.ntyp1) cycle
- itypi1=iabs(itype(i+1))
+ itypi1=iabs(itype(i+1,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
!d write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
!d & 'iend=',iend(i,iint)
do j=istart(i,iint),iend(i,iint)
- itypj=iabs(itype(j))
+ itypj=iabs(itype(j,1))
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
!d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
!d epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
!d write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)')
-!d & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
+!d & restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj),
!d & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm,
!d & (c(k,i),k=1,3),(c(k,j),k=1,3)
evdw=evdw+evdwij
! print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=iabs(itype(i))
+ itypi=iabs(itype(i,1))
if (itypi.eq.ntyp1) cycle
- itypi1=iabs(itype(i+1))
+ itypi1=iabs(itype(i+1,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
!
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
- itypj=iabs(itype(j))
+ itypj=iabs(itype(j,1))
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
!d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
!d epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
!d write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
-!d & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
+!d & restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj),
!d & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
!d & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
!d & (c(k,i),k=1,3),(c(k,j),k=1,3)
! endif
!el ind=0
do i=iatsc_s,iatsc_e
- itypi=iabs(itype(i))
+ itypi=iabs(itype(i,1))
if (itypi.eq.ntyp1) cycle
- itypi1=iabs(itype(i+1))
+ itypi1=iabs(itype(i+1,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
!el ind=ind+1
- itypj=iabs(itype(j))
+ itypj=iabs(itype(j,1))
if (itypj.eq.ntyp1) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
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 & restyp(itypi,1),i,restyp(itypj,1),j,
!d & epsi,sigm,chi1,chi2,chip1,chip2,
!d & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
!d & om1,om2,om12,1.0D0/dsqrt(rrij),
!el ind=0
do i=iatsc_s,iatsc_e
!C print *,"I am in EVDW",i
- itypi=iabs(itype(i))
+ itypi=iabs(itype(i,1))
! if (i.ne.47) cycle
if (itypi.eq.ntyp1) cycle
- itypi1=iabs(itype(i+1))
+ itypi1=iabs(itype(i+1,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
sslipi=0.0d0
ssgradlipi=0.0
endif
- print *, sslipi,ssgradlipi
+! 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))
+ itypj=iabs(itype(j,1))
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,&
! 1.0d0/vbld(j+nres) !d
-! write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
+! write (iout,*) "i",i," j", j," itype",itype(i,1),itype(j,1)
sig0ij=sigma(itypi,itypj)
chi1=chi(itypi,itypj)
chi2=chi(itypj,itypi)
if (rij_shift.le.0.0D0) then
evdw=1.0D20
!d write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-!d & restyp(itypi),i,restyp(itypj),j,
+!d & restyp(itypi,1),i,restyp(itypj,1),j,
!d & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
return
endif
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, &
+ restyp(itypi,1),i,restyp(itypj,1),j, &
epsi,sigm,chi1,chi2,chip1,chip2, &
eps1,eps2rt**2,eps3rt**2,sig,sig0ij, &
om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, &
!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
+! print *,"ZALAMKA", evdw
! Calculate gradient components.
e1=e1*eps1*eps2rt**2*eps3rt**2
enddo ! j
enddo ! iint
enddo ! i
+! print *,"ZALAMKA", evdw
! write (iout,*) "Number of loop steps in EGB:",ind
!ccc energy_dec=.false.
return
! if (icall.eq.0) lprn=.true.
!el ind=0
do i=iatsc_s,iatsc_e
- itypi=iabs(itype(i))
+ itypi=iabs(itype(i,1))
if (itypi.eq.ntyp1) cycle
- itypi1=iabs(itype(i+1))
+ itypi1=iabs(itype(i+1,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
!el ind=ind+1
- itypj=iabs(itype(j))
+ itypj=iabs(itype(j,1))
if (itypj.eq.ntyp1) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
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,&
+ restyp(itypi,1),i,restyp(itypj,1),j,&
epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),&
chi1,chi2,chip1,chip2,&
eps1,eps2rt**2,eps3rt**2,&
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=iabs(itype(i))
+ itypi=iabs(itype(i,1))
if (itypi.eq.ntyp1) cycle
- itypi1=iabs(itype(i+1))
+ itypi1=iabs(itype(i+1,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
!d write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
!d & 'iend=',iend(i,iint)
do j=istart(i,iint),iend(i,iint)
- itypj=iabs(itype(j))
+ itypj=iabs(itype(j,1))
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
eello_turn4=0.0d0
!el ind=0
do i=iatel_s,iatel_e
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
num_conti=0
! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
do j=ielstart(i),ielend(i)
- if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
+ if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) cycle
!el ind=ind+1
iteli=itel(i)
itelj=itel(j)
endif
! if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then
if (i.gt. nnt+2 .and. i.lt.nct+2) then
- iti = itortyp(itype(i-2))
+ iti = itortyp(itype(i-2,1))
else
iti=ntortyp+1
endif
! if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
if (i.gt. nnt+1 .and. i.lt.nct+1) then
- iti1 = itortyp(itype(i-1))
+ iti1 = itortyp(itype(i-1,1))
else
iti1=ntortyp+1
endif
-! print *,iti,i,"iti",iti1,itype(i-1),itype(i-2)
+! print *,iti,i,"iti",iti1,itype(i-1,1),itype(i-2,1)
!d write (iout,*) '*******i',i,' iti1',iti
!d write (iout,*) 'b1',b1(:,iti)
!d write (iout,*) 'b2',b2(:,iti)
enddo
! if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
if (i.gt. nnt+1 .and. i.lt.nct+1) then
- if (itype(i-1).le.ntyp) then
- iti1 = itortyp(itype(i-1))
+ if (itype(i-1,1).le.ntyp) then
+ iti1 = itortyp(itype(i-1,1))
else
iti1=ntortyp+1
endif
#endif
#endif
!d do i=1,nres
-!d iti = itortyp(itype(i))
+!d iti = itortyp(itype(i,1))
!d write (iout,*) i
!d do j=1,2
!d write (iout,'(2f10.5,5x,2f10.5,5x,2f10.5)')
! write (iout,*) 'i',i,' fac',fac
enddo
endif
- print *,wel_loc,"wel_loc",wcorr4,wcorr5,wcorr6,wturn3,wturn4, &
- wturn6
+! 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
! 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
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1 &
+ .or. itype(i+2,1).eq.ntyp1 .or. itype(i+3,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
num_cont_hb(i)=num_conti
enddo
do i=iturn4_start,iturn4_end
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 &
- .or. itype(i+3).eq.ntyp1 &
- .or. itype(i+4).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1 &
+ .or. itype(i+3,1).eq.ntyp1 &
+ .or. itype(i+4,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
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) &
+ if (wturn4.gt.0.0d0 .and. itype(i+2,1).ne.ntyp1) &
call eturn4(i,eello_turn4)
num_cont_hb(i)=num_conti
enddo ! i
!
! Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
!
+ print *,"iatel_s,iatel_e,",iatel_s,iatel_e
do i=iatel_s,iatel_e
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
-! write (iout,*) i,j,itype(i),itype(j)
- if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle
+! write (iout,*) i,j,itype(i,1),itype(j,1)
+ if (itype(j,1).eq.ntyp1.or. itype(j+1,1).eq.ntyp1) cycle
call eelecij(i,j,ees,evdw1,eel_loc)
enddo ! j
num_cont_hb(i)=num_conti
a32=a32*fac
a33=a33*fac
!d write (iout,'(4i5,4f10.5)')
-!d & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
+!d & i,itortyp(itype(i,1)),j,itortyp(itype(j,1)),a22,a23,a32,a33
!d write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
!d write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
!d & uy(:,j),uz(:,j)
a_temp(1,2)=a23
a_temp(2,1)=a32
a_temp(2,2)=a33
- iti1=itortyp(itype(i+1))
- iti2=itortyp(itype(i+2))
- iti3=itortyp(itype(i+3))
+ iti1=itortyp(itype(i+1,1))
+ iti2=itortyp(itype(i+2,1))
+ iti3=itortyp(itype(i+3,1))
! write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3
call transpose2(EUg(1,1,i+1),e1t(1,1))
call transpose2(Eug(1,1,i+2),e2t(1,1))
!d print '(a)','Enter ESCP'
!d write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
do i=iatscp_s,iatscp_e
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
yi=0.5D0*(c(2,i)+c(2,i+1))
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
- if (itype(j).eq.ntyp1) cycle
- itypj=iabs(itype(j))
+ if (itype(j,1).eq.ntyp1) cycle
+ itypj=iabs(itype(j,1))
! Uncomment following three lines for SC-p interactions
! xj=c(1,nres+j)-xi
! yj=c(2,nres+j)-yi
!d print '(a)','Enter ESCP'
!d write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
do i=iatscp_s,iatscp_e
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
yi=0.5D0*(c(2,i)+c(2,i+1))
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
- itypj=iabs(itype(j))
+ itypj=iabs(itype(j,1))
if (itypj.eq.ntyp1) cycle
! Uncomment following three lines for SC-p interactions
! xj=c(1,nres+j)-xi
! 18/07/06 MC: Use the convention that the first nss pairs are SS bonds
if (.not.dyn_ss .and. i.le.nss) then
! 15/02/13 CC dynamic SSbond - additional check
- if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. &
- iabs(itype(jjj)).eq.1) then
+ if (ii.gt.nres .and. iabs(itype(iii,1)).eq.1 .and. &
+ iabs(itype(jjj,1)).eq.1) then
call ssbond_ene(iii,jjj,eij)
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
!-----------------------------------------------------------------------------
deltat1,deltat2,deltat12,ed,pom1,pom2,eom1,eom2,eom12,&
cosphi,ggk
- itypi=iabs(itype(i))
+ itypi=iabs(itype(i,1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
dzi=dc_norm(3,nres+i)
! dsci_inv=dsc_inv(itypi)
dsci_inv=vbld_inv(nres+i)
- itypj=iabs(itype(j))
+ itypj=iabs(itype(j,1))
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(nres+j)
xj=c(1,nres+j)-xi
! if (.not.allocated(gradbx)) allocate(gradbx(3,nres)) !(3,maxres)
do i=ibondp_start,ibondp_end
- if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
- if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
+ if (itype(i-1,1).eq.ntyp1 .and. itype(i,1).eq.ntyp1) cycle
+ if (itype(i-1,1).eq.ntyp1 .or. itype(i,1).eq.ntyp1) then
!C estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
!C do j=1,3
!C gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) &
! endif
enddo
estr=0.5d0*AKP*estr+estr1
+! print *,"estr_bb",estr,AKP
!
! 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
!
do i=ibond_start,ibond_end
- iti=iabs(itype(i))
+ iti=iabs(itype(i,1))
+ if (iti.eq.0) print *,"WARNING WRONG SETTTING",i
if (iti.ne.10 .and. iti.ne.ntyp1) then
nbi=nbondterm(iti)
if (nbi.eq.1) then
"estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,&
AKSC(1,iti),AKSC(1,iti)*diff*diff
estr=estr+0.5d0*AKSC(1,iti)*diff*diff
+! print *,"estr_sc",estr
do j=1,3
gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
enddo
usumsqder=usumsqder+ud(j)*uprod2
enddo
estr=estr+uprod/usum
+! print *,"estr_sc",estr,i
+
+ if (energy_dec) write (iout,*) &
+ "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,&
+ AKSC(1,iti),uprod/usum
do j=1,3
gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres)
enddo
etheta=0.0D0
! write (*,'(a,i2)') 'EBEND ICG=',icg
do i=ithet_start,ithet_end
- if (itype(i-1).eq.ntyp1) cycle
+ if (itype(i-1,1).eq.ntyp1) cycle
! Zero the energy function and its derivative at 0 or pi.
call splinthet(theta(i),0.5d0*delta,ss,ssd)
- it=itype(i-1)
- ichir1=isign(1,itype(i-2))
- ichir2=isign(1,itype(i))
- if (itype(i-2).eq.10) ichir1=isign(1,itype(i-1))
- if (itype(i).eq.10) ichir2=isign(1,itype(i-1))
- if (itype(i-1).eq.10) then
- itype1=isign(10,itype(i-2))
- ichir11=isign(1,itype(i-2))
- ichir12=isign(1,itype(i-2))
- itype2=isign(10,itype(i))
- ichir21=isign(1,itype(i))
- ichir22=isign(1,itype(i))
+ it=itype(i-1,1)
+ ichir1=isign(1,itype(i-2,1))
+ ichir2=isign(1,itype(i,1))
+ if (itype(i-2,1).eq.10) ichir1=isign(1,itype(i-1,1))
+ if (itype(i,1).eq.10) ichir2=isign(1,itype(i-1,1))
+ if (itype(i-1,1).eq.10) then
+ itype1=isign(10,itype(i-2,1))
+ ichir11=isign(1,itype(i-2,1))
+ ichir12=isign(1,itype(i-2,1))
+ itype2=isign(10,itype(i,1))
+ ichir21=isign(1,itype(i,1))
+ ichir22=isign(1,itype(i,1))
endif
- if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+ if (i.gt.3 .and. itype(i-2,1).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
y(1)=0.0D0
y(2)=0.0D0
endif
- if (i.lt.nres .and. itype(i).ne.ntyp1) then
+ if (i.lt.nres .and. itype(i,1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
if (phii1.ne.phii1) phii1=150.0
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 (itype(i-1).eq.ntyp1) cycle
- if (itype(i-2).eq.ntyp1.or.itype(i).eq.ntyp1) cycle
- if (iabs(itype(i+1)).eq.20) iblock=2
- if (iabs(itype(i+1)).ne.20) iblock=1
+ if (itype(i-1,1).eq.ntyp1) cycle
+ if (itype(i-2,1).eq.ntyp1.or.itype(i,1).eq.ntyp1) cycle
+ if (iabs(itype(i+1,1)).eq.20) iblock=2
+ if (iabs(itype(i+1,1)).ne.20) iblock=1
dethetai=0.0d0
dephii=0.0d0
dephii1=0.0d0
theti2=0.5d0*theta(i)
- ityp2=ithetyp((itype(i-1)))
+ ityp2=ithetyp((itype(i-1,1)))
do k=1,nntheterm
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
- if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
+ if (i.gt.3 .and. itype(max0(i-3,1),1).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
#else
phii=phi(i)
#endif
- ityp1=ithetyp((itype(i-2)))
+ ityp1=ithetyp((itype(i-2,1)))
! propagation of chirality for glycine type
do k=1,nsingle
cosph1(k)=dcos(k*phii)
enddo
else
phii=0.0d0
- ityp1=ithetyp(itype(i-2))
+ ityp1=ithetyp(itype(i-2,1))
do k=1,nsingle
cosph1(k)=0.0d0
sinph1(k)=0.0d0
enddo
endif
- if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
+ if (i.lt.nres .and. itype(i+1,1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
if (phii1.ne.phii1) phii1=150.0
#else
phii1=phi(i+1)
#endif
- ityp3=ithetyp((itype(i)))
+ ityp3=ithetyp((itype(i,1)))
do k=1,nsingle
cosph2(k)=dcos(k*phii1)
sinph2(k)=dsin(k*phii1)
enddo
else
phii1=0.0d0
- ityp3=ithetyp(itype(i))
+ ityp3=ithetyp(itype(i,1))
do k=1,nsingle
cosph2(k)=0.0d0
sinph2(k)=0.0d0
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
escloc=0.0D0
! write (iout,'(a)') 'ESC'
do i=loc_start,loc_end
- it=itype(i)
+ it=itype(i,1)
if (it.eq.ntyp1) cycle
if (it.eq.10) goto 1
nlobit=nlob(iabs(it))
delta=0.02d0*pi
escloc=0.0D0
do i=loc_start,loc_end
- if (itype(i).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1) cycle
costtab(i+1) =dcos(theta(i+1))
sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
cosfac=dsqrt(cosfac2)
sinfac2=0.5d0/(1.0d0-costtab(i+1))
sinfac=dsqrt(sinfac2)
- it=iabs(itype(i))
+ it=iabs(itype(i,1))
if (it.eq.10) goto 1
!
! Compute the axes of tghe local cartesian coordinates system; store in
y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
enddo
do j = 1,3
- z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i)))
+ z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i,1)))
enddo
! write (2,*) "i",i
! write (2,*) "x_prime",(x_prime(j),j=1,3)
! Compute the energy of the ith side cbain
!
! write (2,*) "xx",xx," yy",yy," zz",zz
- it=iabs(itype(i))
+ it=iabs(itype(i,1))
do j = 1,65
x(j) = sc_parmin(j,it)
enddo
!c diagnostics - remove later
xx1 = dcos(alph(2))
yy1 = dsin(alph(2))*dcos(omeg(2))
- zz1 = -dsign(1.0,dfloat(itype(i)))*dsin(alph(2))*dsin(omeg(2))
+ zz1 = -dsign(1.0,dfloat(itype(i,1)))*dsin(alph(2))*dsin(omeg(2))
write(2,'(3f8.1,3f9.3,1x,3f9.3)') &
alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,&
xx1,yy1,zz1
! & dscp1,dscp2,sumene
! sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
escloc = escloc + sumene
-! write (2,*) "i",i," escloc",sumene,escloc,it,itype(i)
+! write (2,*) "i",i," escloc",sumene,escloc,it,itype(i,1)
! & ,zz,xx,yy
!#define DEBUG
#ifdef DEBUG
!
! Compute the gradient of esc
!
-! zz=zz*dsign(1.0,dfloat(itype(i)))
+! zz=zz*dsign(1.0,dfloat(itype(i,1)))
pom_s1=(1.0d0+x(63))/(0.1d0 + dscp1)**2
pom_s16=6*(1.0d0+x(64))/(0.1d0 + dscp1**6)**2
pom_s2=(1.0d0+x(65))/(0.1d0 + dscp2)**2
+(sumene2x+sumene4x*cost2tab(i+1))*(s2+s2_6) &
+(pom1+pom2)*pom_dx
#ifdef DEBUG
- write(2,*), "de_dxx = ", de_dxx,de_dxx_num,itype(i)
+ write(2,*), "de_dxx = ", de_dxx,de_dxx_num,itype(i,1)
#endif
!
sumene1y=x(3) + 2*x(6)*yy + x(9)*xx + x(10)*zz
+(sumene2y+sumene4y*cost2tab(i+1))*(s2+s2_6) &
+(pom1-pom2)*pom_dy
#ifdef DEBUG
- write(2,*), "de_dyy = ", de_dyy,de_dyy_num,itype(i)
+ write(2,*), "de_dyy = ", de_dyy,de_dyy_num,itype(i,1)
#endif
!
de_dzz =(x(24) +2*x(27)*zz +x(28)*xx +x(30)*yy &
+x(60)*xx*yy)*cost2tab(i+1)*(s2+s2_6) &
+ ( x(14) + 2*x(17)*zz+ x(18)*xx + x(20)*yy)*(s2+s2_6)
#ifdef DEBUG
- write(2,*), "de_dzz = ", de_dzz,de_dzz_num,itype(i)
+ write(2,*), "de_dzz = ", de_dzz,de_dzz_num,itype(i,1)
#endif
!
de_dt = 0.5d0*sumene3*cost2tab(i+1)*(s1+s1_6) &
-0.5d0*sumene4*sint2tab(i+1)*(s2+s2_6) &
+pom1*pom_dt1+pom2*pom_dt2
#ifdef DEBUG
- write(2,*), "de_dt = ", de_dt,de_dt_num,itype(i)
+ write(2,*), "de_dt = ", de_dt,de_dt_num,itype(i,1)
#endif
!
!
dZZ_Ci(k)=0.0d0
do j=1,3
dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1) &
- *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
+ *dsign(1.0d0,dfloat(itype(i,1)))*dC_norm(j,i+nres)
dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1) &
- *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
+ *dsign(1.0d0,dfloat(itype(i,1)))*dC_norm(j,i+nres)
enddo
dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
etors=0.0D0
do i=iphi_start,iphi_end
etors_ii=0.0D0
- if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 &
- .or. itype(i).eq.ntyp1) cycle
- itori=itortyp(itype(i-2))
- itori1=itortyp(itype(i-1))
+ if (itype(i-2,1).eq.ntyp1.or. itype(i-1,1).eq.ntyp1 &
+ .or. itype(i,1).eq.ntyp1) cycle
+ itori=itortyp(itype(i-2,1))
+ itori1=itortyp(itype(i-1,1))
phii=phi(i)
gloci=0.0D0
! Proline-Proline pair is a special case...
'etor',i,etors_ii
if (lprn) &
write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') &
- restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,&
+ restyp(itype(i-2,1),1),i-2,restyp(itype(i-1,1),1),i-1,itori,itori1,&
(v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6)
gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
! write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
! lprn=.true.
etors=0.0D0
do i=iphi_start,iphi_end
- if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 &
- .or. itype(i-3).eq.ntyp1 &
- .or. itype(i).eq.ntyp1) cycle
+ if (itype(i-2,1).eq.ntyp1 .or. itype(i-1,1).eq.ntyp1 &
+ .or. itype(i-3,1).eq.ntyp1 &
+ .or. itype(i,1).eq.ntyp1) cycle
etors_ii=0.0D0
- if (iabs(itype(i)).eq.20) then
+ if (iabs(itype(i,1)).eq.20) then
iblock=2
else
iblock=1
endif
- itori=itortyp(itype(i-2))
- itori1=itortyp(itype(i-1))
+ itori=itortyp(itype(i-2,1))
+ itori1=itortyp(itype(i-1,1))
phii=phi(i)
gloci=0.0D0
! Regular cosine and sine terms
'etor',i,etors_ii-v0(itori,itori1,iblock)
if (lprn) &
write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') &
- restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,&
+ restyp(itype(i-2,1),1),i-2,restyp(itype(i-1,1),1),i-1,itori,itori1,&
(v1(j,itori,itori1,iblock),j=1,6),&
(v2(j,itori,itori1,iblock),j=1,6)
gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
! write(iout,*) "a tu??"
do i=iphid_start,iphid_end
etors_d_ii=0.0D0
- if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 &
- .or. itype(i-3).eq.ntyp1 &
- .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
- itori=itortyp(itype(i-2))
- itori1=itortyp(itype(i-1))
- itori2=itortyp(itype(i))
+ if (itype(i-2,1).eq.ntyp1 .or. itype(i-1,1).eq.ntyp1 &
+ .or. itype(i-3,1).eq.ntyp1 &
+ .or. itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
+ itori=itortyp(itype(i-2,1))
+ itori1=itortyp(itype(i-1,1))
+ itori2=itortyp(itype(i,1))
phii=phi(i)
phii1=phi(i+1)
gloci1=0.0D0
gloci2=0.0D0
iblock=1
- if (iabs(itype(i+1)).eq.20) iblock=2
+ if (iabs(itype(i+1,1)).eq.20) iblock=2
! Regular cosine and sine terms
do j=1,ntermd_1(itori,itori1,itori2,iblock)
! write (iout,*) "EBACK_SC_COR",itau_start,itau_end
esccor=0.0D0
do i=itau_start,itau_end
- if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
+ if ((itype(i-2,1).eq.ntyp1).or.(itype(i-1,1).eq.ntyp1)) cycle
esccor_ii=0.0D0
- isccori=isccortyp(itype(i-2))
- isccori1=isccortyp(itype(i-1))
+ isccori=isccortyp(itype(i-2,1))
+ isccori1=isccortyp(itype(i-1,1))
! write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1)
phii=phi(i)
! 2 = Ca...Ca...Ca...SC
! 3 = SC...Ca...Ca...SCi
gloci=0.0D0
- if (((intertyp.eq.3).and.((itype(i-2).eq.10).or. &
- (itype(i-1).eq.10).or.(itype(i-2).eq.ntyp1).or. &
- (itype(i-1).eq.ntyp1))) &
- .or. ((intertyp.eq.1).and.((itype(i-2).eq.10) &
- .or.(itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1) &
- .or.(itype(i).eq.ntyp1))) &
- .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or. &
- (itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. &
- (itype(i-3).eq.ntyp1)))) cycle
- if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.ntyp1)) cycle
- if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.ntyp1)) &
+ if (((intertyp.eq.3).and.((itype(i-2,1).eq.10).or. &
+ (itype(i-1,1).eq.10).or.(itype(i-2,1).eq.ntyp1).or. &
+ (itype(i-1,1).eq.ntyp1))) &
+ .or. ((intertyp.eq.1).and.((itype(i-2,1).eq.10) &
+ .or.(itype(i-2,1).eq.ntyp1).or.(itype(i-1,1).eq.ntyp1) &
+ .or.(itype(i,1).eq.ntyp1))) &
+ .or.((intertyp.eq.2).and.((itype(i-1,1).eq.10).or. &
+ (itype(i-1,1).eq.ntyp1).or.(itype(i-2,1).eq.ntyp1).or. &
+ (itype(i-3,1).eq.ntyp1)))) cycle
+ if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1,1).eq.ntyp1)) cycle
+ if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres,1).eq.ntyp1)) &
cycle
do j=1,nterm_sccor(isccori,isccori1)
v1ij=v1sccor(j,intertyp,isccori,isccori1)
gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
if (lprn) &
write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') &
- restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,isccori,isccori1,&
+ restyp(itype(i-2,1),1),i-2,restyp(itype(i-1,1),1),i-1,isccori,isccori1,&
(v1sccor(j,intertyp,isccori,isccori1),j=1,6),&
(v2sccor(j,intertyp,isccori,isccori1),j=1,6)
gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
allocate(dipderx(3,5,4,maxconts,nres))
!
- iti1 = itortyp(itype(i+1))
+ iti1 = itortyp(itype(i+1,1))
if (j.lt.nres-1) then
- itj1 = itortyp(itype(j+1))
+ itj1 = itortyp(itype(j+1,1))
else
itj1=ntortyp+1
endif
if (l.eq.j+1) then
! parallel orientation of the two CA-CA-CA frames.
if (i.gt.1) then
- iti=itortyp(itype(i))
+ iti=itortyp(itype(i,1))
else
iti=ntortyp+1
endif
- itk1=itortyp(itype(k+1))
- itj=itortyp(itype(j))
+ itk1=itortyp(itype(k+1,1))
+ itj=itortyp(itype(j,1))
if (l.lt.nres-1) then
- itl1=itortyp(itype(l+1))
+ itl1=itortyp(itype(l+1,1))
else
itl1=ntortyp+1
endif
else
! Antiparallel orientation of the two CA-CA-CA frames.
if (i.gt.1) then
- iti=itortyp(itype(i))
+ iti=itortyp(itype(i,1))
else
iti=ntortyp+1
endif
- itk1=itortyp(itype(k+1))
- itl=itortyp(itype(l))
- itj=itortyp(itype(j))
+ itk1=itortyp(itype(k+1,1))
+ itl=itortyp(itype(l,1))
+ itj=itortyp(itype(j,1))
if (j.lt.nres-1) then
- itj1=itortyp(itype(j+1))
+ itj1=itortyp(itype(j+1,1))
else
itj1=ntortyp+1
endif
!d write (iout,*)
!d & 'EELLO5: Contacts have occurred for peptide groups',i,j,
!d & ' and',k,l
- itk=itortyp(itype(k))
- itl=itortyp(itype(l))
- itj=itortyp(itype(j))
+ itk=itortyp(itype(k,1))
+ itl=itortyp(itype(l,1))
+ itj=itortyp(itype(j,1))
eello5_1=0.0d0
eello5_2=0.0d0
eello5_3=0.0d0
! i i C
! C
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
- itk=itortyp(itype(k))
+ itk=itortyp(itype(k,1))
s1= scalar2(AEAb1(1,2,imat),CUgb2(1,i))
s2=-scalar2(AEAb2(1,1,imat),Ug2Db1t(1,k))
s3= scalar2(AEAb2(1,1,imat),CUgb2(1,k))
!
! 4/7/01 AL Component s1 was removed, because it pertains to the respective
! energy moment and not to the cluster cumulant.
- iti=itortyp(itype(i))
+ iti=itortyp(itype(i,1))
if (j.lt.nres-1) then
- itj1=itortyp(itype(j+1))
+ itj1=itortyp(itype(j+1,1))
else
itj1=ntortyp+1
endif
- itk=itortyp(itype(k))
- itk1=itortyp(itype(k+1))
+ itk=itortyp(itype(k,1))
+ itk1=itortyp(itype(k+1,1))
if (l.lt.nres-1) then
- itl1=itortyp(itype(l+1))
+ itl1=itortyp(itype(l+1,1))
else
itl1=ntortyp+1
endif
! 4/7/01 AL Component s1 was removed, because it pertains to the respective
! energy moment and not to the cluster cumulant.
!d write (2,*) 'eello_graph4: wturn6',wturn6
- iti=itortyp(itype(i))
- itj=itortyp(itype(j))
+ iti=itortyp(itype(i,1))
+ itj=itortyp(itype(j,1))
if (j.lt.nres-1) then
- itj1=itortyp(itype(j+1))
+ itj1=itortyp(itype(j+1,1))
else
itj1=ntortyp+1
endif
- itk=itortyp(itype(k))
+ itk=itortyp(itype(k,1))
if (k.lt.nres-1) then
- itk1=itortyp(itype(k+1))
+ itk1=itortyp(itype(k+1,1))
else
itk1=ntortyp+1
endif
- itl=itortyp(itype(l))
+ itl=itortyp(itype(l,1))
if (l.lt.nres-1) then
- itl1=itortyp(itype(l+1))
+ itl1=itortyp(itype(l+1,1))
else
itl1=ntortyp+1
endif
j=i+4
k=i+1
l=i+3
- iti=itortyp(itype(i))
- itk=itortyp(itype(k))
- itk1=itortyp(itype(k+1))
- itl=itortyp(itype(l))
- itj=itortyp(itype(j))
+ iti=itortyp(itype(i,1))
+ itk=itortyp(itype(k,1))
+ itk1=itortyp(itype(k+1,1))
+ itl=itortyp(itype(l,1))
+ itj=itortyp(itype(j,1))
!d write (2,*) 'itk',itk,' itk1',itk1,' itl',itl,' itj',itj
!d write (2,*) 'i',i,' k',k,' j',j,' l',l
!d if (i.ne.1 .or. j.ne.3 .or. k.ne.2 .or. l.ne.4) then
wturn6*gcorr6_turn_long(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)
-
-
-
+ +wtube*gg_tube(j,i) &
+ +wbond_nucl*gradb_nucl(j,i)
enddo
enddo
#else
wturn6*gcorr6_turn_long(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)
-
-
+ +wtube*gg_tube(j,i) &
+ +wbond_nucl*gradb_nucl(j,i)
enddo
enddo
wsccor*gsccorc(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) &
+wturn4*gshieldc_loc_t4(j,i) &
+wel_loc*gshieldc_ll(j,i) &
+wel_loc*gshieldc_loc_ll(j,i) &
- +wtube*gg_tube(j,i)
+ +wtube*gg_tube(j,i) &
+ +wbond_nucl*gradb_nucl(j,i)
+
#else
wturn6*gcorr6_turn(j,i)+ &
wsccor*gsccorc(j,i) &
+wscloc*gscloc(j,i) &
+ +gradafm(j,i) &
+wliptran*gliptranc(j,i) &
+welec*gshieldc(j,i) &
+welec*gshieldc_loc(j,) &
+wturn4*gshieldc_loc_t4(j,i) &
+wel_loc*gshieldc_ll(j,i) &
+wel_loc*gshieldc_loc_ll(j,i) &
- +wtube*gg_tube(j,i)
+ +wtube*gg_tube(j,i) &
+ +wbond_nucl*gradb_nucl(j,i)
+
+wturn3*gshieldx_t3(j,i) &
+wturn4*gshieldx_t4(j,i) &
+wel_loc*gshieldx_ll(j,i)&
- +wtube*gg_tube_sc(j,i)
+ +wtube*gg_tube_sc(j,i) &
+ +wbond_nucl*gradbx_nucl(j,i)
+
enddo
! include 'COMMON.CALC'
! include 'COMMON.IOUNITS'
real(kind=8), dimension(3) :: dcosom1,dcosom2
-
+! print *,"wchodze"
eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 &
! Derivatives in alpha and omega:
!
do i=2,nres-1
-! dsci=dsc(itype(i))
+! dsci=dsc(itype(i,1))
dsci=vbld(i+nres)
#ifdef OSF
alphi=alph(i)
! write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
!d write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
!d & 'iend=',iend(i,iint)
do j=istart(i,iint),iend(i,iint)
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
! write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
!d write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
!d & 'iend=',iend(i,iint)
do j=istart(i,iint),iend(i,iint)
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
! print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
!
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
!d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
!d epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
!d write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
-!d & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
+!d & restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj),
!d & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
!d & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
!d & (c(k,i),k=1,3),(c(k,j),k=1,3)
! print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
!
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
!d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
!d epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
!d write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)')
-!d & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
+!d & restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj),
!d & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
!d & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
!d & (c(k,i),k=1,3),(c(k,j),k=1,3)
! endif
!el ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
!el ind=ind+1
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
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 & restyp(itypi,1),i,restyp(itypj,1),j,
!d & epsi,sigm,chi1,chi2,chip1,chip2,
!d & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
!d & om1,om2,om12,1.0D0/dsqrt(rrij),
! endif
!el ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
!el ind=ind+1
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
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 & restyp(itypi,1),i,restyp(itypj,1),j,
!d & epsi,sigm,chi1,chi2,chip1,chip2,
!d & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq),
!d & om1,om2,om12,1.0D0/dsqrt(rrij),
! if (icall.eq.0) lprn=.false.
!el ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
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'
+! 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)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
! write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
! & 1.0d0/vbld(j+nres)
-! write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
+! write (iout,*) "i",i," j", j," itype",itype(i,1),itype(j,1)
sig0ij=sigma(itypi,itypj)
chi1=chi(itypi,itypj)
chi2=chi(itypj,itypi)
if (rij_shift.le.0.0D0) then
evdw=1.0D20
!d write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-!d & restyp(itypi),i,restyp(itypj),j,
+!d & restyp(itypi,1),i,restyp(itypj,1),j,
!d & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
return
endif
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,&
+ restyp(itypi,1),i,restyp(itypj,1),j,&
epsi,sigm,chi1,chi2,chip1,chip2,&
eps1,eps2rt**2,eps3rt**2,sig,sig0ij,&
om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,&
! if (icall.eq.0) lprn=.false.
!el ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
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)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
! write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
! & 1.0d0/vbld(j+nres)
-! write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
+! write (iout,*) "i",i," j", j," itype",itype(i,1),itype(j,1)
sig0ij=sigma(itypi,itypj)
chi1=chi(itypi,itypj)
chi2=chi(itypj,itypi)
if (rij_shift.le.0.0D0) then
evdw=1.0D20
!d write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-!d & restyp(itypi),i,restyp(itypj),j,
+!d & restyp(itypi,1),i,restyp(itypj,1),j,
!d & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
return
endif
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,&
+ restyp(itypi,1),i,restyp(itypj,1),j,&
epsi,sigm,chi1,chi2,chip1,chip2,&
eps1,eps2rt**2,eps3rt**2,sig,sig0ij,&
om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,&
! if (icall.eq.0) lprn=.true.
!el ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
!el ind=ind+1
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
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,&
+ restyp(itypi,1),i,restyp(itypj,1),j,&
epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),&
chi1,chi2,chip1,chip2,&
eps1,eps2rt**2,eps3rt**2,&
! if (icall.eq.0) lprn=.true.
!el ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=itype(i,1)
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=itype(i+1,1)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
!el ind=ind+1
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
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,&
+ restyp(itypi,1),i,restyp(itypj,1),j,&
epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),&
chi1,chi2,chip1,chip2,&
eps1,eps2rt**2,eps3rt**2,&
! Loop over i,i+2 and i,i+3 pairs of the peptide groups
!
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
+ if (itype(i,1).eq.ntyp1.or. itype(i+1,1).eq.ntyp1 &
+ .or. itype(i+2,1).eq.ntyp1 .or. itype(i+3,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
num_cont_hb(i)=num_conti
enddo
do i=iturn4_start,iturn4_end
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 &
- .or. itype(i+3).eq.ntyp1 &
- .or. itype(i+4).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1 &
+ .or. itype(i+3,1).eq.ntyp1 &
+ .or. itype(i+4,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
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) &
+ if (wturn4.gt.0.0d0 .and. itype(i+2,1).ne.ntyp1) &
call eturn4(i,eello_turn4)
num_cont_hb(i)=num_conti
enddo ! i
! Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
!
do i=iatel_s,iatel_e
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
- if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
+ if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) cycle
call eelecij_scale(i,j,ees,evdw1,eel_loc)
enddo ! j
num_cont_hb(i)=num_conti
a32=a32*fac
a33=a33*fac
!d write (iout,'(4i5,4f10.5)')
-!d & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
+!d & i,itortyp(itype(i,1)),j,itortyp(itype(j,1)),a22,a23,a32,a33
!d write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
!d write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i),
!d & uy(:,j),uz(:,j)
! & " iatel_e_vdw",iatel_e_vdw
call flush(iout)
do i=iatel_s_vdw,iatel_e_vdw
- if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1.or. itype(i+1,1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
! & ' ielend',ielend_vdw(i)
call flush(iout)
do j=ielstart_vdw(i),ielend_vdw(i)
- if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
+ if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) cycle
!el ind=ind+1
iteli=itel(i)
itelj=itel(j)
!d print '(a)','Enter ESCP'
!d write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
do i=iatscp_s,iatscp_e
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
yi=0.5D0*(c(2,i)+c(2,i+1))
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
! Uncomment following three lines for SC-p interactions
! xj=c(1,nres+j)-xi
!d print '(a)','Enter ESCP'
!d write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
do i=iatscp_s,iatscp_e
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
yi=0.5D0*(c(2,i)+c(2,i+1))
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
- itypj=itype(j)
+ itypj=itype(j,1)
if (itypj.eq.ntyp1) cycle
! Uncomment following three lines for SC-p interactions
! xj=c(1,nres+j)-xi
!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
enddo
if (n.le.nphi+ntheta) goto 10
do i=2,nres-1
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
galphai=0.0D0
gomegai=0.0D0
do k=1,3
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
+ gradb_nucl(j,i)=0.0d0
+ gradbx_nucl(j,i)=0.0d0
do intertyp=1,3
gloc_sc(intertyp,i,icg)=0.0d0
enddo
do j=1,3
dcostheta(j,1,i)=-(dc_norm(j,i-1)+cost*dc_norm(j,i-2))/&
vbld(i-1)
- if (itype(i-1).ne.ntyp1) dtheta(j,1,i)=-dcostheta(j,1,i)/sint
+ if (itype(i-1,1).ne.ntyp1) dtheta(j,1,i)=-dcostheta(j,1,i)/sint
dcostheta(j,2,i)=-(dc_norm(j,i-2)+cost*dc_norm(j,i-1))/&
vbld(i)
- if (itype(i-1).ne.ntyp1) dtheta(j,2,i)=-dcostheta(j,2,i)/sint
+ if (itype(i-1,1).ne.ntyp1) dtheta(j,2,i)=-dcostheta(j,2,i)/sint
enddo
enddo
#if defined(MPI) && defined(PARINTDER)
#else
do i=3,nres
#endif
- if ((itype(i-1).ne.10).and.(itype(i-1).ne.ntyp1)) then
+ if ((itype(i-1,1).ne.10).and.(itype(i-1,1).ne.ntyp1)) then
cost1=dcos(omicron(1,i))
sint1=sqrt(1-cost1*cost1)
cost2=dcos(omicron(2,i))
dcosomicron(j,2,2,i)=-(dc_norm(j,i-1) &
+cost2*(-dc_norm(j,i-1+nres)))/ &
vbld(i-1+nres)
-! write(iout,*) "vbld", i,itype(i),vbld(i-1+nres)
+! write(iout,*) "vbld", i,itype(i,1),vbld(i-1+nres)
domicron(j,2,2,i)=-1/sint2*dcosomicron(j,2,2,i)
enddo
endif
#else
do i=4,nres
#endif
-! if (itype(i-1).eq.21 .or. itype(i-2).eq.21 ) cycle
+! if (itype(i-1,1).eq.21 .or. itype(i-2,1).eq.21 ) cycle
! the conventional case
sint=dsin(theta(i))
sint1=dsin(theta(i-1))
ctgt=cost/sint
ctgt1=cost1/sint1
cosg_inv=1.0d0/cosg
- if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then
+ if (itype(i-1,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1) then
dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1) &
-(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2)
dphi(j,1,i)=cosg_inv*dsinphi(j,1,i)
! Obtaining the gamma derivatives from cosine derivative
else
do j=1,3
- if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then
+ if (itype(i-1,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1) then
dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3* &
dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp* &
dc_norm(j,i-3))/vbld(i-2)
do i=3,nres
!elwrite(iout,*) " vecpr",i,nres
#endif
- if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle
-! if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10).or.
-! & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1)) cycle
+ if ((itype(i-2,1).eq.ntyp1).or.(itype(i-2,1).eq.10)) cycle
+! if ((itype(i-2,1).eq.ntyp1).or.(itype(i-2,1).eq.10).or.
+! & (itype(i-1,1).eq.ntyp1).or.(itype(i,1).eq.ntyp1)) cycle
!c dtauangle(j,intertyp,dervityp,residue number)
!c INTERTYP=1 SC...Ca...Ca..Ca
! the conventional case
#else
do i=4,nres
#endif
- if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or. &
- (itype(i-2).eq.ntyp1).or.(itype(i-3).eq.ntyp1)) cycle
+ if ((itype(i-1,1).eq.ntyp1).or.(itype(i-1,1).eq.10).or. &
+ (itype(i-2,1).eq.ntyp1).or.(itype(i-3,1).eq.ntyp1)) cycle
! the conventional case
sint=dsin(omicron(1,i))
sint1=dsin(theta(i-1))
do i=3,nres
#endif
! the conventional case
- if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or. &
- (itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle
+ if ((itype(i-1,1).eq.ntyp1).or.(itype(i-1,1).eq.10).or. &
+ (itype(i-2,1).eq.ntyp1).or.(itype(i-2,1).eq.10)) cycle
sint=dsin(omicron(1,i))
sint1=dsin(omicron(2,i-1))
sing=dsin(tauangle(3,i))
#else
do i=2,nres-1
#endif
- if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
fac5=1.0d0/dsqrt(2*(1+dcos(theta(i+1))))
fac6=fac5/vbld(i)
fac7=fac5*fac5
write (iout,*) &
"Analytical (upper) and numerical (lower) gradient of alpha"
do i=2,nres-1
- if(itype(i).ne.10) then
+ if(itype(i,1).ne.10) then
do j=1,3
dcji=dc(j,i-1)
dc(j,i-1)=dcji+aincr
write (iout,*) &
"Analytical (upper) and numerical (lower) gradient of omega"
do i=2,nres-1
- if(itype(i).ne.10) then
+ if(itype(i,1).ne.10) then
do j=1,3
dcji=dc(j,i-1)
dc(j,i-1)=dcji+aincr
(cref(3,jl,kkk)-cref(3,il,kkk))**2)
dij=dist(il,jl)
qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
- if (itype(il).ne.10 .or. itype(jl).ne.10) then
+ if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then
nl=nl+1
d0ijCM=dsqrt( &
(cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
(cref(3,jl,kkk)-cref(3,il,kkk))**2)
dij=dist(il,jl)
qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
- if (itype(il).ne.10 .or. itype(jl).ne.10) then
+ if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then
nl=nl+1
d0ijCM=dsqrt( &
(cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
dqwol(k,jl)=dqwol(k,jl)-ddqij
enddo
- if (itype(il).ne.10 .or. itype(jl).ne.10) then
+ if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then
nl=nl+1
d0ijCM=dsqrt( &
(cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
dqwol(k,il)=dqwol(k,il)+ddqij
dqwol(k,jl)=dqwol(k,jl)-ddqij
enddo
- if (itype(il).ne.10 .or. itype(jl).ne.10) then
+ if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then
nl=nl+1
d0ijCM=dsqrt( &
(cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
!el allocate(dyn_ssbond_ij(iatsc_s:iatsc_e,nres))
!el allocate(dyn_ssbond_ij(0:nres+4,nres))
- itypi=itype(i)
+ itypi=itype(i,1)
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
dsci_inv=vbld_inv(i+nres)
- itypj=itype(j)
+ itypj=itype(j,1)
xj=c(1,nres+j)-c(1,nres+i)
yj=c(2,nres+j)-c(2,nres+i)
zj=c(3,nres+j)-c(3,nres+i)
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,1)
+ 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,1)
+ 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,1)
+ 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)
real(kind=8) :: fracinbuf,eliptran,sslip,positi,ssgradlip
integer :: i
eliptran=0.0
- print *, "I am in eliptran"
+! 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))&
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1).or.(i.eq.nres))&
cycle
positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize))
enddo
! here starts the side chain transfer
do i=ilip_start,ilip_end
- if (itype(i).eq.ntyp1) cycle
+ if (itype(i,1).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 lipbufthick is thickenes of lipid buffore
sslip=sscalelip(fracinbuf)
ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
- eliptran=eliptran+sslip*liptranene(itype(i))
+ eliptran=eliptran+sslip*liptranene(itype(i,1))
gliptranx(3,i)=gliptranx(3,i) &
- +ssgradlip*liptranene(itype(i))
+ +ssgradlip*liptranene(itype(i,1))
gliptranc(3,i-1)= gliptranc(3,i-1) &
- +ssgradlip*liptranene(itype(i))
+ +ssgradlip*liptranene(itype(i,1))
!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))
+ eliptran=eliptran+sslip*liptranene(itype(i,1))
gliptranx(3,i)=gliptranx(3,i) &
- +ssgradlip*liptranene(itype(i))
+ +ssgradlip*liptranene(itype(i,1))
gliptranc(3,i-1)= gliptranc(3,i-1) &
- +ssgradlip*liptranene(itype(i))
+ +ssgradlip*liptranene(itype(i,1))
!C print *, "doing sscalefor top part",sslip,fracinbuf
else
- eliptran=eliptran+liptranene(itype(i))
+ eliptran=eliptran+liptranene(itype(i,1))
!C print *,"I am in true lipid"
endif
endif ! if in lipid or buffor
!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),dimension(3) :: vectube
real(kind=8) :: Etube,xtemp,xminact,yminact,&
ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi, &
sc_aa_tube,sc_bb_tube
!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
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle
!C now calculate distance from center of tube and direction vectors
xmin=boxxsize
ymin=boxysize
do i=itube_start,itube_end
!C Lets not jump over memory as we use many times iti
- iti=itype(i)
+ iti=itype(i,1)
!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 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),dimension(3) :: vectube
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
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
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,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 lipbufthick is thickenes of lipid buffore
sstube=sscalelip(fracinbuf)
ssgradtube=-sscagradlip(fracinbuf)/tubebufthick
-!C print *,ssgradtube, sstube,tubetranene(itype(i))
+!C print *,ssgradtube, sstube,tubetranene(itype(i,1))
enetube(i)=enetube(i)+sstube*tubetranenepep
!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
-!C &+ssgradtube*tubetranene(itype(i))
+!C &+ssgradtube*tubetranene(itype(i,1))
!C gg_tube(3,i-1)= gg_tube(3,i-1)
-!C &+ssgradtube*tubetranene(itype(i))
+!C &+ssgradtube*tubetranene(itype(i,1))
!C print *,"doing sccale for lower part"
elseif (positi.gt.buftubetop) then
fracinbuf=1.0d0- &
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 &+ssgradtube*tubetranene(itype(i,1))
!C gg_tube(3,i-1)= gg_tube(3,i-1)
-!C &+ssgradtube*tubetranene(itype(i))
+!C &+ssgradtube*tubetranene(itype(i,1))
!C print *, "doing sscalefor top part",sslip,fracinbuf
else
sstube=1.0d0
!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)
+ iti=itype(i,1)
!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 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 print *,ssgradtube, sstube,tubetranene(itype(i,1))
+ enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i,1))
!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
-!C &+ssgradtube*tubetranene(itype(i))
+!C &+ssgradtube*tubetranene(itype(i,1))
!C gg_tube(3,i-1)= gg_tube(3,i-1)
-!C &+ssgradtube*tubetranene(itype(i))
+!C &+ssgradtube*tubetranene(itype(i,1))
!C print *,"doing sccale for lower part"
elseif (positi.gt.buftubetop) then
fracinbuf=1.0d0- &
sstube=sscalelip(fracinbuf)
ssgradtube=sscagradlip(fracinbuf)/tubebufthick
- enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i))
+ enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i,1))
!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
-!C &+ssgradtube*tubetranene(itype(i))
+!C &+ssgradtube*tubetranene(itype(i,1))
!C gg_tube(3,i-1)= gg_tube(3,i-1)
-!C &+ssgradtube*tubetranene(itype(i))
+!C &+ssgradtube*tubetranene(itype(i,1))
!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))
+ enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i,1))
!C print *,"I am in true lipid"
endif
else
end subroutine calctube2
!=====================================================================================================================================
subroutine calcnano(Etube)
- real(kind=8) :: vectube(3),enetube(nres*2), &
- enecavtube(nres*2)
+ real(kind=8),dimension(3) :: vectube
+
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
+ integer:: i,j,iti,r
Etube=0.0d0
- print *,itube_start,itube_end,"poczatek"
+! print *,itube_start,itube_end,"poczatek"
do i=itube_start,itube_end
enetube(i)=0.0d0
enetube(i+nres)=0.0d0
!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
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle
!C now calculate distance from center of tube and direction vectors
xmin=boxxsize
ymin=boxysize
!C fac=fac+faccav
!C 667 continue
endif
-
+ if (energy_dec) write(iout,*),i,rdiff,enetube(i),enecavtube(i)
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
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)
+ iti=itype(i,1)
!C lets ommit dummy atoms for now
if ((iti.eq.ntyp1) &
!C in UNRES uncomment the line below as GLY has no side-chain...
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
+ if (energy_dec) write(iout,*),i,rdiff,enetube(i+nres),enecavtube(i+nres)
enddo
Etube=Etube+enetube(i)+enetube(i+nres)+enecavtube(i) &
+enecavtube(i+nres)
enddo
+! do i=1,20
+! print *,"begin", i,"a"
+! do r=1,10000
+! rdiff=r/100.0d0
+! rdiff6=rdiff**6.0d0
+! sc_aa_tube=sc_aa_tube_par(i)
+! sc_bb_tube=sc_bb_tube_par(i)
+! enetube(i)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6
+! denominator=(1.0d0+dcavtub(i)*rdiff6*rdiff6)
+! enecavtube(i)= &
+! (bcavtub(i)*rdiff+acavtub(i)*dsqrt(rdiff)+ccavtub(i)) &
+! /denominator
+
+! print '(5(f10.3,1x))',rdiff,enetube(i),enecavtube(i),enecavtube(i)+enetube(i)
+! enddo
+! print *,"end",i,"a"
+! enddo
!C print *,"ETUBE", etube
return
end subroutine calcnano
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
+!C if ((itype(i,1).eq.ntyp1).and.itype(i+1,1).eq.ntyp1) cycle
ishield_list(i)=0
- if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+ if ((itype(i,1).eq.ntyp1).and.itype(i+1,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
+ if ((itype(k,1).eq.ntyp1).or.(itype(k,1).eq.10)) cycle
dist_pep_side=0.0
dist_side_calf=0.0
do j=1,3
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))
+ short=short_r_sidechain(itype(k,1))
+ long=long_r_sidechain(itype(k,1))
costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2)
sinthet=short/dist_pep_side*costhet
!C now costhet_grad
enddo
fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield)
-!C write(2,*) "TOTAL VOLUME",i,itype(i),fac_shield(i)
+!C write(2,*) "TOTAL VOLUME",i,itype(i,1),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
allocate(ielstart_vdw(nres))
allocate(ielend_vdw(nres))
!(maxres)
+ allocate(nint_gr_nucl(nres))
+ allocate(nscp_gr_nucl(nres))
+ allocate(ielstart_nucl(nres))
+ allocate(ielend_nucl(nres))
+!(maxres)
+ allocate(istart_nucl(nres,maxint_gr))
+ allocate(iend_nucl(nres,maxint_gr))
+!(maxres,maxint_gr)
+ allocate(iscpstart_nucl(nres,maxint_gr))
+ allocate(iscpend_nucl(nres,maxint_gr))
+!(maxres,maxint_gr)
+ allocate(ielstart_vdw_nucl(nres))
+ allocate(ielend_vdw_nucl(nres))
allocate(lentyp(0:nfgtasks-1))
!(0:maxprocs-1)
allocate(grad_shield(3,-1:nres))
allocate(gg_tube_sc(3,-1:nres))
allocate(gg_tube(3,-1:nres))
+ allocate(gradafm(3,-1:nres))
+ allocate(gradb_nucl(3,-1:nres))
+ allocate(gradbx_nucl(3,-1:nres))
+ allocate(gvdwpsb1(3,-1:nres))
+ allocate(gelpp(3,-1:nres))
+ allocate(gvdwpsb(3,-1:nres))
+ allocate(gelsbc(3,-1:nres))
+ allocate(gelsbx(3,-1:nres))
+ allocate(gvdwsbx(3,-1:nres))
+ allocate(gvdwsbc(3,-1:nres))
+ allocate(gsbloc(3,-1:nres))
+ allocate(gsblocx(3,-1:nres))
+ allocate(gradcorr_nucl(3,-1:nres))
+ allocate(gradxorr_nucl(3,-1:nres))
+ allocate(gradcorr3_nucl(3,-1:nres))
+ allocate(gradxorr3_nucl(3,-1:nres))
+
!(3,maxres)
allocate(grad_shield_side(3,50,nres))
allocate(grad_shield_loc(3,50,nres))
! 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))
+ allocate(enetube(nres*2))
+ allocate(enecavtube(nres*2))
+
!(maxres)
dyn_ss_mask(:)=.false.
!----------------------
return
end subroutine alloc_ener_arrays
+!-----------------------------------------------------------------
+ subroutine ebond_nucl(estr_nucl)
+!c
+!c Evaluate the energy of stretching of the CA-CA and CA-SC virtual bonds
+!c
+
+ real(kind=8),dimension(3) :: u,ud
+ real(kind=8) :: usum,uprod,uprod1,uprod2,usumsqder
+ real(kind=8) :: estr_nucl,diff
+ integer :: iti,i,j,k,nbi
+ estr_nucl=0.0d0
+!C print *,"I enter ebond"
+ if (energy_dec) &
+ write (iout,*) "ibondp_start,ibondp_end",&
+ ibondp_nucl_start,ibondp_nucl_end
+ do i=ibondp_nucl_start,ibondp_nucl_end
+ if (itype(i-1,2).eq.ntyp1_molec(2) .or. &
+ itype(i,2).eq.ntyp1_molec(2)) cycle
+! estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
+! do j=1,3
+! gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax)
+! & *dc(j,i-1)/vbld(i)
+! enddo
+! if (energy_dec) write(iout,*)
+! & "estr1",i,vbld(i),distchainmax,
+! & gnmr1(vbld(i),-1.0d0,distchainmax)
+
+ diff = vbld(i)-vbldp0_nucl
+ if(energy_dec)write(iout,*) "estr_nucl_bb" , i,vbld(i),&
+ vbldp0_nucl,diff,AKP_nucl*diff*diff
+ estr_nucl=estr_nucl+diff*diff
+ print *,estr_nucl
+ do j=1,3
+ gradb_nucl(j,i-1)=AKP_nucl*diff*dc(j,i-1)/vbld(i)
+ enddo
+!c write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3)
+ enddo
+ estr_nucl=0.5d0*AKP_nucl*estr_nucl
+ print *,"partial sum", estr_nucl,AKP_nucl
+
+ if (energy_dec) &
+ write (iout,*) "ibondp_start,ibondp_end",&
+ ibond_nucl_start,ibond_nucl_end
+
+ do i=ibond_nucl_start,ibond_nucl_end
+!C print *, "I am stuck",i
+ iti=itype(i,2)
+ if (iti.eq.ntyp1_molec(2)) cycle
+ nbi=nbondterm_nucl(iti)
+!C print *,iti,nbi
+ if (nbi.eq.1) then
+ diff=vbld(i+nres)-vbldsc0_nucl(1,iti)
+
+ if (energy_dec) &
+ write (iout,*) "estr_nucl_sc", i,iti,vbld(i+nres),vbldsc0_nucl(1,iti),diff, &
+ AKSC_nucl(1,iti),AKSC_nucl(1,iti)*diff*diff
+ estr_nucl=estr_nucl+0.5d0*AKSC_nucl(1,iti)*diff*diff
+ print *,estr_nucl
+ do j=1,3
+ gradbx_nucl(j,i)=AKSC_nucl(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
+ enddo
+ else
+ do j=1,nbi
+ diff=vbld(i+nres)-vbldsc0_nucl(j,iti)
+ ud(j)=aksc_nucl(j,iti)*diff
+ u(j)=abond0_nucl(j,iti)+0.5d0*ud(j)*diff
+ enddo
+ uprod=u(1)
+ do j=2,nbi
+ uprod=uprod*u(j)
+ enddo
+ usum=0.0d0
+ usumsqder=0.0d0
+ do j=1,nbi
+ uprod1=1.0d0
+ uprod2=1.0d0
+ do k=1,nbi
+ if (k.ne.j) then
+ uprod1=uprod1*u(k)
+ uprod2=uprod2*u(k)*u(k)
+ endif
+ enddo
+ usum=usum+uprod1
+ usumsqder=usumsqder+ud(j)*uprod2
+ enddo
+ estr_nucl=estr_nucl+uprod/usum
+ do j=1,3
+ gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres)
+ enddo
+ endif
+ enddo
+!C print *,"I am about to leave ebond"
+ return
+ end subroutine ebond_nucl
+
+!-----------------------------------------------------------------------------
+ subroutine ebend_nucl(etheta_nucl)
+ real(kind=8),dimension(nntheterm_nucl+1) :: coskt,sinkt !mmaxtheterm
+ real(kind=8),dimension(nsingle_nucl+1) :: cosph1,sinph1,cosph2,sinph2 !maxsingle
+ real(kind=8),dimension(ndouble_nucl+1,ndouble_nucl+1) :: cosph1ph2,sinph1ph2 !maxdouble,maxdouble
+ logical :: lprn=.true., lprn1=.false.
+!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_nucl,ccl,ssl,scl,csl,ethetacnstr
+! local variables for constrains
+ real(kind=8) :: difi,thetiii
+ integer itheta
+ etheta_nucl=0.0D0
+! print *,"ithet_start",ithet_nucl_start," ithet_end",ithet_nucl_end,nres
+ do i=ithet_nucl_start,ithet_nucl_end
+ if ((itype(i-1,2).eq.ntyp1_molec(2)).or.&
+ (itype(i-2,2).eq.ntyp1_molec(2)).or. &
+ (itype(i,2).eq.ntyp1_molec(2))) cycle
+ dethetai=0.0d0
+ dephii=0.0d0
+ dephii1=0.0d0
+ theti2=0.5d0*theta(i)
+ ityp2=ithetyp_nucl(itype(i-1,2))
+ do k=1,nntheterm_nucl
+ coskt(k)=dcos(k*theti2)
+ sinkt(k)=dsin(k*theti2)
+ enddo
+ if (i.gt.3 .and. itype(i-2,2).ne.ntyp1_molec(2)) then
+#ifdef OSF
+ phii=phi(i)
+ if (phii.ne.phii) phii=150.0
+#else
+ phii=phi(i)
+#endif
+ ityp1=ithetyp_nucl(itype(i-2,2))
+ do k=1,nsingle_nucl
+ cosph1(k)=dcos(k*phii)
+ sinph1(k)=dsin(k*phii)
+ enddo
+ else
+ phii=0.0d0
+ ityp1=nthetyp_nucl+1
+ do k=1,nsingle_nucl
+ cosph1(k)=0.0d0
+ sinph1(k)=0.0d0
+ enddo
+ endif
+
+ if (i.lt.nres .and. itype(i,2).ne.ntyp1_molec(2)) then
+#ifdef OSF
+ phii1=phi(i+1)
+ if (phii1.ne.phii1) phii1=150.0
+ phii1=pinorm(phii1)
+#else
+ phii1=phi(i+1)
+#endif
+ ityp3=ithetyp_nucl(itype(i,2))
+ do k=1,nsingle_nucl
+ cosph2(k)=dcos(k*phii1)
+ sinph2(k)=dsin(k*phii1)
+ enddo
+ else
+ phii1=0.0d0
+ ityp3=nthetyp_nucl+1
+ do k=1,nsingle_nucl
+ cosph2(k)=0.0d0
+ sinph2(k)=0.0d0
+ enddo
+ endif
+ ethetai=aa0thet_nucl(ityp1,ityp2,ityp3)
+ do k=1,ndouble_nucl
+ do l=1,k-1
+ ccl=cosph1(l)*cosph2(k-l)
+ ssl=sinph1(l)*sinph2(k-l)
+ scl=sinph1(l)*cosph2(k-l)
+ csl=cosph1(l)*sinph2(k-l)
+ cosph1ph2(l,k)=ccl-ssl
+ cosph1ph2(k,l)=ccl+ssl
+ sinph1ph2(l,k)=scl+csl
+ sinph1ph2(k,l)=scl-csl
+ enddo
+ enddo
+ if (lprn) then
+ write (iout,*) "i",i," ityp1",ityp1," ityp2",ityp2,&
+ " ityp3",ityp3," theti2",theti2," phii",phii," phii1",phii1
+ write (iout,*) "coskt and sinkt",nntheterm_nucl
+ do k=1,nntheterm_nucl
+ write (iout,*) k,coskt(k),sinkt(k)
+ enddo
+ endif
+ do k=1,ntheterm_nucl
+ ethetai=ethetai+aathet_nucl(k,ityp1,ityp2,ityp3)*sinkt(k)
+ dethetai=dethetai+0.5d0*k*aathet_nucl(k,ityp1,ityp2,ityp3)&
+ *coskt(k)
+ if (lprn)&
+ write (iout,*) "k",k," aathet",aathet_nucl(k,ityp1,ityp2,ityp3),&
+ " ethetai",ethetai
+ enddo
+ if (lprn) then
+ write (iout,*) "cosph and sinph"
+ do k=1,nsingle_nucl
+ write (iout,*) k,cosph1(k),sinph1(k),cosph2(k),sinph2(k)
+ enddo
+ write (iout,*) "cosph1ph2 and sinph2ph2"
+ do k=2,ndouble_nucl
+ do l=1,k-1
+ write (iout,*) l,k,cosph1ph2(l,k),cosph1ph2(k,l),&
+ sinph1ph2(l,k),sinph1ph2(k,l)
+ enddo
+ enddo
+ write(iout,*) "ethetai",ethetai
+ endif
+ do m=1,ntheterm2_nucl
+ do k=1,nsingle_nucl
+ aux=bbthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph1(k)&
+ +ccthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph1(k)&
+ +ddthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph2(k)&
+ +eethet_nucl(k,m,ityp1,ityp2,ityp3)*sinph2(k)
+ ethetai=ethetai+sinkt(m)*aux
+ dethetai=dethetai+0.5d0*m*aux*coskt(m)
+ dephii=dephii+k*sinkt(m)*(&
+ ccthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph1(k)-&
+ bbthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph1(k))
+ dephii1=dephii1+k*sinkt(m)*(&
+ eethet_nucl(k,m,ityp1,ityp2,ityp3)*cosph2(k)-&
+ ddthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph2(k))
+ if (lprn) &
+ write (iout,*) "m",m," k",k," bbthet",&
+ bbthet_nucl(k,m,ityp1,ityp2,ityp3)," ccthet",&
+ ccthet_nucl(k,m,ityp1,ityp2,ityp3)," ddthet",&
+ ddthet_nucl(k,m,ityp1,ityp2,ityp3)," eethet",&
+ eethet_nucl(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+ enddo
+ enddo
+ if (lprn) &
+ write(iout,*) "ethetai",ethetai
+ do m=1,ntheterm3_nucl
+ do k=2,ndouble_nucl
+ do l=1,k-1
+ aux=ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+&
+ ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+&
+ ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+&
+ ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)
+ ethetai=ethetai+sinkt(m)*aux
+ dethetai=dethetai+0.5d0*m*coskt(m)*aux
+ dephii=dephii+l*sinkt(m)*(&
+ -ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)-&
+ ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+&
+ ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+&
+ ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+ dephii1=dephii1+(k-l)*sinkt(m)*( &
+ -ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+&
+ ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+&
+ ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)-&
+ ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+ if (lprn) then
+ write (iout,*) "m",m," k",k," l",l," ffthet", &
+ ffthet_nucl(l,k,m,ityp1,ityp2,ityp3), &
+ ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)," ggthet",&
+ ggthet_nucl(l,k,m,ityp1,ityp2,ityp3),&
+ ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+ write (iout,*) cosph1ph2(l,k)*sinkt(m), &
+ cosph1ph2(k,l)*sinkt(m),&
+ sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m)
+ endif
+ enddo
+ enddo
+ enddo
+10 continue
+ if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') &
+ i,theta(i)*rad2deg,phii*rad2deg, &
+ phii1*rad2deg,ethetai
+ etheta_nucl=etheta_nucl+ethetai
+! print *,i,"partial sum",etheta_nucl
+ if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang_nucl*dephii
+ if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang_nucl*dephii1
+ gloc(nphi+i-2,icg)=wang_nucl*dethetai
+ enddo
+ return
+ end subroutine ebend_nucl
+!----------------------------------------------------
+ subroutine etor_nucl(etors_nucl)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.VAR'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCAL'
+! include 'COMMON.TORSION'
+! include 'COMMON.INTERACT'
+! include 'COMMON.DERIV'
+! include 'COMMON.CHAIN'
+! include 'COMMON.NAMES'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.FFIELD'
+! include 'COMMON.TORCNSTR'
+! include 'COMMON.CONTROL'
+ real(kind=8) :: etors_nucl,edihcnstr
+ logical :: lprn
+!el local variables
+ integer :: i,j,iblock,itori,itori1
+ real(kind=8) :: phii,gloci,v1ij,v2ij,cosphi,sinphi,&
+ vl1ij,vl2ij,vl3ij,pom1,difi,etors_ii,pom
+! Set lprn=.true. for debugging
+ lprn=.false.
+! lprn=.true.
+ etors_nucl=0.0D0
+! print *,"iphi_nucl_start/end", iphi_nucl_start,iphi_nucl_end
+ do i=iphi_nucl_start,iphi_nucl_end
+ if (itype(i-2,2).eq.ntyp1_molec(2) .or. itype(i-1,2).eq.ntyp1_molec(2) &
+ .or. itype(i-3,2).eq.ntyp1_molec(2) &
+ .or. itype(i,2).eq.ntyp1_molec(2)) cycle
+ etors_ii=0.0D0
+ itori=itortyp_nucl(itype(i-2,2))
+ itori1=itortyp_nucl(itype(i-1,2))
+ phii=phi(i)
+! print *,i,itori,itori1
+ gloci=0.0D0
+!C Regular cosine and sine terms
+ do j=1,nterm_nucl(itori,itori1)
+ v1ij=v1_nucl(j,itori,itori1)
+ v2ij=v2_nucl(j,itori,itori1)
+ cosphi=dcos(j*phii)
+ sinphi=dsin(j*phii)
+ etors_nucl=etors_nucl+v1ij*cosphi+v2ij*sinphi
+ if (energy_dec) etors_ii=etors_ii+&
+ v1ij*cosphi+v2ij*sinphi
+ gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
+ enddo
+!C Lorentz terms
+!C v1
+!C E = SUM ----------------------------------- - v1
+!C [v2 cos(phi/2)+v3 sin(phi/2)]^2 + 1
+!C
+ cosphi=dcos(0.5d0*phii)
+ sinphi=dsin(0.5d0*phii)
+ do j=1,nlor_nucl(itori,itori1)
+ vl1ij=vlor1_nucl(j,itori,itori1)
+ vl2ij=vlor2_nucl(j,itori,itori1)
+ vl3ij=vlor3_nucl(j,itori,itori1)
+ pom=vl2ij*cosphi+vl3ij*sinphi
+ pom1=1.0d0/(pom*pom+1.0d0)
+ etors_nucl=etors_nucl+vl1ij*pom1
+ if (energy_dec) etors_ii=etors_ii+ &
+ vl1ij*pom1
+ pom=-pom*pom1*pom1
+ gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom
+ enddo
+!C Subtract the constant term
+ etors_nucl=etors_nucl-v0_nucl(itori,itori1)
+ if (energy_dec) write (iout,'(a6,i5,0pf7.3)') &
+ 'etor',i,etors_ii-v0_nucl(itori,itori1)
+ if (lprn) &
+ write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') &
+ restyp(itype(i-2,2),2),i-2,restyp(itype(i-1,2),2),i-1,itori,itori1, &
+ (v1_nucl(j,itori,itori1),j=1,6),(v2_nucl(j,itori,itori1),j=1,6)
+ gloc(i-3,icg)=gloc(i-3,icg)+wtor_nucl*gloci
+!c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
+ enddo
+ return
+ end subroutine etor_nucl
+!------------------------------------------------------------
+ subroutine epp_nucl_sub(evdw1,ees)
+!C
+!C This subroutine calculates the average interaction energy and its gradient
+!C in the virtual-bond vectors between non-adjacent peptide groups, based on
+!C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
+!C The potential depends both on the distance of peptide-group centers and on
+!C the orientation of the CA-CA virtual bonds.
+!C
+ integer :: i,j,k,iteli,itelj,num_conti,isubchap,ind
+ 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,fac,evdw1ij
+ integer xshift,yshift,zshift
+ real(kind=8),dimension(3):: ggg,gggp,gggm,erij
+ real(kind=8) :: ees,eesij
+!c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
+ real(kind=8) scal_el /0.5d0/
+ t_eelecij=0.0d0
+ ees=0.0D0
+ evdw1=0.0D0
+ ind=0
+!c
+!c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
+!c
+ print *,"iatel_s_nucl,iatel_e_nucl",iatel_s_nucl,iatel_e_nucl
+ do i=iatel_s_nucl,iatel_e_nucl
+ if (itype(i,2).eq.ntyp1_molec(2) .or. itype(i+1,2).eq.ntyp1_molec(2)) cycle
+ dxi=dc(1,i)
+ dyi=dc(2,i)
+ dzi=dc(3,i)
+ dx_normi=dc_norm(1,i)
+ dy_normi=dc_norm(2,i)
+ dz_normi=dc_norm(3,i)
+ 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
+
+ do j=ielstart_nucl(i),ielend_nucl(i)
+ if (itype(j,2).eq.ntyp1_molec(2) .or. itype(j+1,2).eq.ntyp1_molec(2)) cycle
+ ind=ind+1
+ dxj=dc(1,j)
+ dyj=dc(2,j)
+ dzj=dc(3,j)
+! xj=c(1,j)+0.5D0*dxj-xmedi
+! yj=c(2,j)+0.5D0*dyj-ymedi
+! zj=c(3,j)+0.5D0*dzj-zmedi
+ xj=c(1,j)+0.5D0*dxj
+ yj=c(2,j)+0.5D0*dyj
+ zj=c(3,j)+0.5D0*dzj
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ 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
+!c write (2,*)"ij",i,j," r0pp",r0pp," rij",rij," epspp",epspp
+ fac=(r0pp**2/rij)**3
+ ev1=epspp*fac*fac
+ ev2=epspp*fac
+ evdw1ij=ev1-2*ev2
+ fac=(-ev1-evdw1ij)/rij
+! write (2,*)"fac",fac," ev1",ev1," ev2",ev2," evdw1ij",evdw1ij
+ if (energy_dec) write(iout,'(2i5,a9,f10.4)') i,j,"evdw1ij",evdw1ij
+ evdw1=evdw1+evdw1ij
+!C
+!C Calculate contributions to the Cartesian gradient.
+!C
+ ggg(1)=fac*xj
+ ggg(2)=fac*yj
+ ggg(3)=fac*zj
+ do k=1,3
+ gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
+ gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
+ enddo
+!c phoshate-phosphate electrostatic interactions
+ rij=dsqrt(rij)
+ fac=1.0d0/rij
+ eesij=dexp(-BEES*rij)*fac
+! write (2,*)"fac",fac," eesijpp",eesij
+ if (energy_dec) write(iout,'(2i5,a9,f10.4)') i,j,"eesijpp",eesij
+ ees=ees+eesij
+!c fac=-eesij*fac
+ fac=-(fac+BEES)*eesij*fac
+ ggg(1)=fac*xj
+ ggg(2)=fac*yj
+ ggg(3)=fac*zj
+!c write(2,*) "ggg",i,j,ggg(1),ggg(2),ggg(3)
+!c write(2,*) "gelpp",i,(gelpp(k,i),k=1,3)
+!c write(2,*) "gelpp",j,(gelpp(k,j),k=1,3)
+ do k=1,3
+ gelpp(k,i)=gelpp(k,i)-ggg(k)
+ gelpp(k,j)=gelpp(k,j)+ggg(k)
+ enddo
+ enddo ! j
+ enddo ! i
+!c ees=332.0d0*ees
+ ees=AEES*ees
+ do i=nnt,nct
+!c write (2,*) "i",i," gelpp",(gelpp(k,i),k=1,3)
+ do k=1,3
+ gvdwpp(k,i)=6*gvdwpp(k,i)
+!c gelpp(k,i)=332.0d0*gelpp(k,i)
+ gelpp(k,i)=AEES*gelpp(k,i)
+ enddo
+!c write (2,*) "i",i," gelpp",(gelpp(k,i),k=1,3)
+ enddo
+!c write (2,*) "total EES",ees
+ return
+ end subroutine epp_nucl_sub
+!---------------------------------------------------------------------
+ subroutine epsb(evdwpsb,eelpsb)
+! use comm_locel
+!C
+!C This subroutine calculates the excluded-volume interaction energy between
+!C peptide-group centers and side chains and its gradient in virtual-bond and
+!C side-chain vectors.
+!C
+ real(kind=8),dimension(3):: ggg
+ 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,rij,evdwpsb,eelpsb
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
+ integer xshift,yshift,zshift
+
+!cd print '(a)','Enter ESCP'
+!cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
+ eelpsb=0.0d0
+ evdwpsb=0.0d0
+ print *,"iatscp_s_nucl,iatscp_e_nucl",iatscp_s_nucl,iatscp_e_nucl
+ do i=iatscp_s_nucl,iatscp_e_nucl
+ if (itype(i,2).eq.ntyp1_molec(2) &
+ .or. itype(i+1,2).eq.ntyp1_molec(2)) cycle
+ 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_nucl(i)
+
+ do j=iscpstart_nucl(i,iint),iscpend_nucl(i,iint)
+ itypj=itype(j,2)
+ if (itypj.eq.ntyp1_molec(2)) cycle
+!C Uncomment following three lines for SC-p interactions
+!c xj=c(1,nres+j)-xi
+!c yj=c(2,nres+j)-yi
+!c zj=c(3,nres+j)-zi
+!C Uncomment following three lines for Ca-p interactions
+! xj=c(1,j)-xi
+! yj=c(2,j)-yi
+! zj=c(3,j)-zi
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ 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)
+ fac=rrij**expon2
+ e1=fac*fac*aad_nucl(itypj)
+ e2=fac*bad_nucl(itypj)
+ if (iabs(j-i) .le. 2) then
+ e1=scal14*e1
+ e2=scal14*e2
+ endif
+ evdwij=e1+e2
+ evdwpsb=evdwpsb+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a4)') &
+ 'evdw2',i,j,evdwij,"tu4"
+!C
+!C Calculate contributions to the gradient in the virtual-bond and SC vectors.
+!C
+ fac=-(evdwij+e1)*rrij
+ ggg(1)=xj*fac
+ ggg(2)=yj*fac
+ ggg(3)=zj*fac
+ do k=1,3
+ gvdwpsb1(k,i)=gvdwpsb1(k,i)-ggg(k)
+ gvdwpsb(k,j)=gvdwpsb(k,j)+ggg(k)
+ enddo
+ enddo
+
+ enddo ! iint
+ enddo ! i
+ do i=1,nct
+ do j=1,3
+ gvdwpsb(j,i)=expon*gvdwpsb(j,i)
+ gvdwpsb1(j,i)=expon*gvdwpsb1(j,i)
+ enddo
+ enddo
+ return
+ end subroutine epsb
+
+!------------------------------------------------------
+ subroutine esb_gb(evdwsb,eelsb)
+ use comm_locel
+ use calc_data_nucl
+ integer :: iint,itypi,itypi1,itypj,subchap,num_conti2
+ real(kind=8) :: xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
+ real(kind=8) :: evdw,sig0iji,evdwsb,eelsb,ecorr,eelij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,aa,bb,faclip,sig0ij
+ integer :: ii
+ logical lprn
+ evdw=0.0D0
+ eelsb=0.0d0
+ ecorr=0.0d0
+ evdwsb=0.0D0
+ lprn=.false.
+ ind=0
+! print *,"iastsc_nucl",iatsc_s_nucl,iatsc_e_nucl
+ do i=iatsc_s_nucl,iatsc_e_nucl
+ num_conti=0
+ num_conti2=0
+ itypi=itype(i,2)
+! PRINT *,"I=",i,itypi
+ if (itypi.eq.ntyp1_molec(2)) cycle
+ itypi1=itype(i+1,2)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ xi=dmod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=dmod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=dmod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+ dsci_inv=vbld_inv(i+nres)
+!C
+!C Calculate SC interaction energy.
+!C
+ do iint=1,nint_gr_nucl(i)
+! print *,"tu?",i,istart_nucl(i,iint),iend_nucl(i,iint)
+ do j=istart_nucl(i,iint),iend_nucl(i,iint)
+ ind=ind+1
+! print *,"JESTEM"
+ itypj=itype(j,2)
+ if (itypj.eq.ntyp1_molec(2)) cycle
+ dscj_inv=vbld_inv(j+nres)
+ sig0ij=sigma_nucl(itypi,itypj)
+ chi1=chi_nucl(itypi,itypj)
+ chi2=chi_nucl(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip_nucl(itypi,itypj)
+ chip2=chip_nucl(itypj,itypi)
+ chip12=chip1*chip2
+! xj=c(1,nres+j)-xi
+! yj=c(2,nres+j)-yi
+! zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(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
+
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+!C Calculate angle-dependent terms of energy and contributions to their
+!C derivatives.
+ erij(1)=xj*rij
+ erij(2)=yj*rij
+ erij(3)=zj*rij
+ om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
+ om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
+ om12=dxi*dxj+dyi*dyj+dzi*dzj
+ call sc_angular_nucl
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+sig0ij
+! print *,rij_shift,"rij_shift"
+!c write (2,*) " rij",1.0D0/rij," sig",sig," sig0ij",sig0ij,
+!c & " rij_shift",rij_shift
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D20
+ return
+ endif
+ sigder=-sig*sigsq
+!c---------------------------------------------------------------
+ rij_shift=1.0D0/rij_shift
+ fac=rij_shift**expon
+ e1=fac*fac*aa_nucl(itypi,itypj)
+ e2=fac*bb_nucl(itypi,itypj)
+ evdwij=eps1*eps2rt*(e1+e2)
+!c write (2,*) "eps1",eps1," eps2rt",eps2rt,
+!c & " e1",e1," e2",e2," evdwij",evdwij
+ eps2der=evdwij
+ evdwij=evdwij*eps2rt
+ evdwsb=evdwsb+evdwij
+ if (lprn) then
+ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
+ restyp(itypi,2),i,restyp(itypj,2),j, &
+ epsi,sigm,chi1,chi2,chip1,chip2, &
+ eps1,eps2rt**2,sig,sig0ij, &
+ om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,&
+ evdwij
+ write (iout,*) "aa",aa_nucl(itypi,itypj)," bb",bb_nucl(itypi,itypj)
+ endif
+
+ if (energy_dec) write (iout,'(a6,2i5,e15.3,a4)') &
+ 'evdw',i,j,evdwij,"tu3"
+
+
+!C Calculate gradient components.
+ e1=e1*eps1*eps2rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ fac=rij*fac
+!c fac=0.0d0
+!C Calculate the radial part of the gradient
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+!C Calculate angular part of the gradient.
+ call sc_grad_nucl
+ call eelsbij(eelij,num_conti2)
+ if (energy_dec .and. &
+ (j.eq.i+1.or.j.eq.nres-i+1.or.j.eq.nres-i.or.j.eq.nres-i+2)) &
+ write (istat,'(e14.5)') evdwij
+ eelsb=eelsb+eelij
+ enddo ! j
+ enddo ! iint
+ num_cont_hb(i)=num_conti2
+ enddo ! i
+!c write (iout,*) "Number of loop steps in EGB:",ind
+!cccc energy_dec=.false.
+ return
+ end subroutine esb_gb
+!-------------------------------------------------------------------------------
+ subroutine eelsbij(eesij,num_conti2)
+ use comm_locel
+ use calc_data_nucl
+ real(kind=8),dimension(3) :: ggg,gggp,gggm,dcosb,dcosg
+ real(kind=8),dimension(3,3) :: erder,uryg,urzg,vryg,vrzg
+ 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,num_conti2
+
+!c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
+ real(kind=8) scal_el /0.5d0/
+ integer :: iteli,itelj,kkk,kkll,m,isubchap
+ real(kind=8) :: ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp,facfac
+ real(kind=8) :: ees,evdw1,eel_loc,aaa,bbb,ael3i,ael63i,ael32i
+ real(kind=8) :: dx_normj,dy_normj,dz_normj,&
+ r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,fac5,fac6,&
+ el1,el2,el3,el4,eesij,ees0ij,facvdw,facel,fac1,ecosa,&
+ ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der,&
+ a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij,&
+ ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp,&
+ ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,&
+ ecosgp,ecosam,ecosbm,ecosgm,ghalf,itypi,itypj
+ ind=ind+1
+ itypi=itype(i,2)
+ itypj=itype(j,2)
+! print *,i,j,itypi,itypj,istype(i),istype(j),"????"
+ ael6i=ael6_nucl(itypi,itypj)
+ ael3i=ael3_nucl(itypi,itypj)
+ ael63i=ael63_nucl(itypi,itypj)
+ ael32i=ael32_nucl(itypi,itypj)
+!c write (iout,*) "eelecij",i,j,itype(i),itype(j),
+!c & ael6i,ael3i,ael63i,al32i,rij,rrij
+ dxj=dc(1,j+nres)
+ dyj=dc(2,j+nres)
+ dzj=dc(3,j+nres)
+ dx_normi=dc_norm(1,i+nres)
+ dy_normi=dc_norm(2,i+nres)
+ dz_normi=dc_norm(3,i+nres)
+ dx_normj=dc_norm(1,j+nres)
+ dy_normj=dc_norm(2,j+nres)
+ dz_normj=dc_norm(3,j+nres)
+!c xj=c(1,j)+0.5D0*dxj-xmedi
+!c yj=c(2,j)+0.5D0*dyj-ymedi
+!c zj=c(3,j)+0.5D0*dzj-zmedi
+ if (ipot_nucl.ne.2) then
+ cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
+ cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
+ cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
+ else
+ cosa=om12
+ cosb=om1
+ cosg=om2
+ endif
+ r3ij=rij*rrij
+ r6ij=r3ij*r3ij
+ fac=cosa-3.0D0*cosb*cosg
+ facfac=fac*fac
+ fac1=3.0d0*(cosb*cosb+cosg*cosg)
+ fac3=ael6i*r6ij
+ fac4=ael3i*r3ij
+ fac5=ael63i*r6ij
+ fac6=ael32i*r6ij
+!c write (iout,*) "r3ij",r3ij," r6ij",r6ij," fac",fac," fac1",fac1,
+!c & " fac2",fac2," fac3",fac3," fac4",fac4," fac5",fac5," fac6",fac6
+ el1=fac3*(4.0D0+facfac-fac1)
+ el2=fac4*fac
+ el3=fac5*(2.0d0-2.0d0*facfac+fac1)
+ el4=fac6*facfac
+ eesij=el1+el2+el3+el4
+!C 12/26/95 - for the evaluation of multi-body H-bonding interactions
+ ees0ij=4.0D0+facfac-fac1
+
+ if (energy_dec) then
+ if(j.eq.i+1.or.j.eq.nres-i+1.or.j.eq.nres-i.or.j.eq.nres-i+2) &
+ write (istat,'(2a1,i4,1x,2a1,i4,4f10.5,3e12.5,$)') &
+ sugartyp(istype(i)),restyp(itypi,2),i,sugartyp(istype(j)),&
+ restyp(itypj,2),j,1.0d0/rij,cosa,cosb,cosg,fac*r3ij, &
+ (4.0D0+facfac-fac1)*r6ij,(2.0d0-2.0d0*facfac+fac1)*r6ij
+ write (iout,'(a6,2i5,e15.3)') 'ees',i,j,eesij
+ endif
+
+!C
+!C Calculate contributions to the Cartesian gradient.
+!C
+ facel=-3.0d0*rrij*(eesij+el1+el3+el4)
+ fac1=fac
+!c erij(1)=xj*rmij
+!c erij(2)=yj*rmij
+!c erij(3)=zj*rmij
+!*
+!* Radial derivatives. First process both termini of the fragment (i,j)
+!*
+ ggg(1)=facel*xj
+ ggg(2)=facel*yj
+ ggg(3)=facel*zj
+ do k=1,3
+ gelsbc(k,j)=gelsbc(k,j)+ggg(k)
+ gelsbc(k,i)=gelsbc(k,i)-ggg(k)
+ gelsbx(k,j)=gelsbx(k,j)+ggg(k)
+ gelsbx(k,i)=gelsbx(k,i)-ggg(k)
+ enddo
+!*
+!* Angular part
+!*
+ ecosa=2.0D0*fac3*fac1+fac4+(-4.0d0*fac5+2.0d0*fac6)*fac1
+ fac4=-3.0D0*fac4
+ fac3=-6.0D0*fac3
+ fac5= 6.0d0*fac5
+ fac6=-6.0d0*fac6
+ ecosb=fac3*(fac1*cosg+cosb)+cosg*fac4+(cosb+2*fac1*cosg)*fac5+&
+ fac6*fac1*cosg
+ ecosg=fac3*(fac1*cosb+cosg)+cosb*fac4+(cosg+2*fac1*cosb)*fac5+&
+ fac6*fac1*cosb
+ do k=1,3
+ dcosb(k)=rij*(dc_norm(k,i+nres)-erij(k)*cosb)
+ dcosg(k)=rij*(dc_norm(k,j+nres)-erij(k)*cosg)
+ enddo
+ do k=1,3
+ ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
+ enddo
+ do k=1,3
+ gelsbx(k,i)=gelsbx(k,i)-ggg(k) &
+ +(ecosa*(dc_norm(k,j+nres)-cosa*dc_norm(k,i+nres))&
+ + ecosb*(erij(k)-cosb*dc_norm(k,i+nres)))*vbld_inv(i+nres)
+ gelsbx(k,j)=gelsbx(k,j)+ggg(k) &
+ +(ecosa*(dc_norm(k,i+nres)-cosa*dc_norm(k,j+nres))&
+ + ecosg*(erij(k)-cosg*dc_norm(k,j+nres)))*vbld_inv(j+nres)
+ gelsbc(k,j)=gelsbc(k,j)+ggg(k)
+ gelsbc(k,i)=gelsbc(k,i)-ggg(k)
+ enddo
+! IF ( (wcorr_nucl.gt.0.0d0.or.wcorr3_nucl.gt.0.0d0) .and.
+ IF ( j.gt.i+1 .and.&
+ num_conti.le.maxconts) THEN
+!C
+!C Calculate the contact function. The ith column of the array JCONT will
+!C contain the numbers of atoms that make contacts with the atom I (of numbers
+!C greater than I). The arrays FACONT and GACONT will contain the values of
+!C the contact function and its derivative.
+ r0ij=2.20D0*sigma(itypi,itypj)
+!c write (2,*) "ij",i,j," rij",1.0d0/rij," r0ij",r0ij
+ call gcont(rij,r0ij,1.0D0,0.2d0/r0ij,fcont,fprimcont)
+!c write (2,*) "fcont",fcont
+ if (fcont.gt.0.0D0) then
+ num_conti=num_conti+1
+ num_conti2=num_conti2+1
+
+ if (num_conti.gt.maxconts) then
+ write (iout,*) 'WARNING - max. # of contacts exceeded;',&
+ ' will skip next contacts for this conf.'
+ else
+ jcont_hb(num_conti,i)=j
+!c write (iout,*) "num_conti",num_conti,
+!c & " jcont_hb",jcont_hb(num_conti,i)
+!C Calculate contact energies
+ cosa4=4.0D0*cosa
+ wij=cosa-3.0D0*cosb*cosg
+ cosbg1=cosb+cosg
+ cosbg2=cosb-cosg
+ fac3=dsqrt(-ael6i)*r3ij
+!c write (2,*) "ael6i",ael6i," r3ij",r3ij," fac3",fac3
+ ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
+ if (ees0tmp.gt.0) then
+ ees0pij=dsqrt(ees0tmp)
+ else
+ ees0pij=0
+ endif
+ ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
+ if (ees0tmp.gt.0) then
+ ees0mij=dsqrt(ees0tmp)
+ else
+ ees0mij=0
+ endif
+ ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
+ ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
+!c write (iout,*) "i",i," j",j,
+!c & " ees0m",ees0m(num_conti,i)," ees0p",ees0p(num_conti,i)
+ ees0pij1=fac3/ees0pij
+ ees0mij1=fac3/ees0mij
+ fac3p=-3.0D0*fac3*rrij
+ ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
+ ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
+ ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
+ ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
+ ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
+ ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
+ ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
+ ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
+ ecosap=ecosa1+ecosa2
+ ecosbp=ecosb1+ecosb2
+ ecosgp=ecosg1+ecosg2
+ ecosam=ecosa1-ecosa2
+ ecosbm=ecosb1-ecosb2
+ ecosgm=ecosg1-ecosg2
+!C End diagnostics
+ facont_hb(num_conti,i)=fcont
+ fprimcont=fprimcont/rij
+ do k=1,3
+ 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
+!C Derivatives due to the contact function
+ gacont_hbr(1,num_conti,i)=fprimcont*xj
+ gacont_hbr(2,num_conti,i)=fprimcont*yj
+ gacont_hbr(3,num_conti,i)=fprimcont*zj
+ do k=1,3
+!c
+!c Gradient of the correlation terms
+!c
+ gacontp_hb1(k,num_conti,i)= &
+ (ecosap*(dc_norm(k,j+nres)-cosa*dc_norm(k,i+nres)) &
+ + ecosbp*(erij(k)-cosb*dc_norm(k,i+nres)))*vbld_inv(i+nres)
+ gacontp_hb2(k,num_conti,i)= &
+ (ecosap*(dc_norm(k,i+nres)-cosa*dc_norm(k,j+nres)) &
+ + ecosgp*(erij(k)-cosg*dc_norm(k,j+nres)))*vbld_inv(j+nres)
+ gacontp_hb3(k,num_conti,i)=gggp(k)
+ gacontm_hb1(k,num_conti,i)= &
+ (ecosam*(dc_norm(k,j+nres)-cosa*dc_norm(k,i+nres)) &
+ + ecosbm*(erij(k)-cosb*dc_norm(k,i+nres)))*vbld_inv(i+nres)
+ gacontm_hb2(k,num_conti,i)= &
+ (ecosam*(dc_norm(k,i+nres)-cosa*dc_norm(k,j+nres))&
+ + ecosgm*(erij(k)-cosg*dc_norm(k,j+nres)))*vbld_inv(j+nres)
+ gacontm_hb3(k,num_conti,i)=gggm(k)
+ enddo
+ endif
+ endif
+ ENDIF
+ return
+ end subroutine eelsbij
+!------------------------------------------------------------------
+ subroutine sc_grad_nucl
+ use comm_locel
+ use calc_data_nucl
+ real(kind=8),dimension(3) :: dcosom1,dcosom2
+ eom1=eps2der*eps2rt_om1+sigder*sigsq_om1
+ eom2=eps2der*eps2rt_om2+sigder*sigsq_om2
+ eom12=evdwij*eps1_om12+eps2der*eps2rt_om12+sigder*sigsq_om12
+ do k=1,3
+ dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
+ dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
+ enddo
+ do k=1,3
+ gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
+ enddo
+ do k=1,3
+ gvdwsbx(k,i)=gvdwsbx(k,i)-gg(k) &
+ +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))&
+ +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+ gvdwsbx(k,j)=gvdwsbx(k,j)+gg(k) &
+ +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
+ +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ enddo
+!C
+!C Calculate the components of the gradient in DC and X
+!C
+ do l=1,3
+ gvdwsbc(l,i)=gvdwsbc(l,i)-gg(l)
+ gvdwsbc(l,j)=gvdwsbc(l,j)+gg(l)
+ enddo
+ return
+ end subroutine sc_grad_nucl
+!-----------------------------------------------------------------------
+ subroutine esb(esbloc)
+!C Calculate the local energy of a side chain and its derivatives in the
+!C corresponding virtual-bond valence angles THETA and the spherical angles
+!C ALPHA and OMEGA derived from AM1 all-atom calculations.
+!C added by Urszula Kozlowska. 07/11/2007
+!C
+ real(kind=8),dimension(3):: x_prime,y_prime,z_prime
+ real(kind=8),dimension(9):: x
+ real(kind=8) :: sumene,dsc_i,dp2_i,xx,yy,zz,sumene1, &
+ sumene2,sumene3,sumene4,s1,s1_6,s2,s2_6,&
+ de_dxx,de_dyy,de_dzz,de_dt,s1_t,s1_6_t,s2_t,s2_6_t
+ real(kind=8),dimension(3):: dXX_Ci1,dYY_Ci1,dZZ_Ci1,dXX_Ci,&
+ dYY_Ci,dZZ_Ci,dXX_XYZ,dYY_XYZ,dZZ_XYZ,dt_dCi,dt_dCi1
+ real(kind=8) :: esbloc,delta,cosfac2,cosfac,sinfac2,sinfac,de_dtt,&
+ cossc,cossc1,cosfac2xx,sinfac2yy,pom1,pom
+ integer::it,nlobit,i,j,k
+! common /sccalc/ time11,time12,time112,theti,it,nlobit
+ delta=0.02d0*pi
+ esbloc=0.0D0
+ do i=loc_start_nucl,loc_end_nucl
+ if (itype(i,2).eq.ntyp1_molec(2)) cycle
+ costtab(i+1) =dcos(theta(i+1))
+ sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
+ cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
+ sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
+ cosfac2=0.5d0/(1.0d0+costtab(i+1))
+ cosfac=dsqrt(cosfac2)
+ sinfac2=0.5d0/(1.0d0-costtab(i+1))
+ sinfac=dsqrt(sinfac2)
+ it=itype(i,2)
+ if (it.eq.10) goto 1
+
+!c
+!C Compute the axes of tghe local cartesian coordinates system; store in
+!c x_prime, y_prime and z_prime
+!c
+ do j=1,3
+ x_prime(j) = 0.00
+ y_prime(j) = 0.00
+ z_prime(j) = 0.00
+ enddo
+!C write(2,*) "dc_norm", dc_norm(1,i+nres),dc_norm(2,i+nres),
+!C & dc_norm(3,i+nres)
+ do j = 1,3
+ x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
+ y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
+ enddo
+ do j = 1,3
+ z_prime(j) = -uz(j,i-1)
+ enddo
+
+ xx=0.0d0
+ yy=0.0d0
+ zz=0.0d0
+ do j = 1,3
+ xx = xx + x_prime(j)*dc_norm(j,i+nres)
+ yy = yy + y_prime(j)*dc_norm(j,i+nres)
+ zz = zz + z_prime(j)*dc_norm(j,i+nres)
+ enddo
+
+ xxtab(i)=xx
+ yytab(i)=yy
+ zztab(i)=zz
+ it=itype(i,2)
+ do j = 1,9
+ x(j) = sc_parmin_nucl(j,it)
+ enddo
+#ifdef CHECK_COORD
+!Cc diagnostics - remove later
+ xx1 = dcos(alph(2))
+ yy1 = dsin(alph(2))*dcos(omeg(2))
+ zz1 = -dsin(alph(2))*dsin(omeg(2))
+ write(2,'(3f8.1,3f9.3,1x,3f9.3)') &
+ alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,&
+ xx1,yy1,zz1
+!C," --- ", xx_w,yy_w,zz_w
+!c end diagnostics
+#endif
+ sumene = enesc_nucl(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ esbloc = esbloc + sumene
+ if (energy_dec) write(iout,*) "i",i," esbloc",sumene,esbloc,xx,yy,zz
+ if (energy_dec) write(iout,*) "x",(x(k),k=1,9)
+#ifdef DEBUG
+ write (2,*) "x",(x(k),k=1,9)
+!C
+!C This section to check the numerical derivatives of the energy of ith side
+!C chain in xx, yy, zz, and theta. Use the -DDEBUG compiler option or insert
+!C #define DEBUG in the code to turn it on.
+!C
+ write (2,*) "sumene =",sumene
+ aincr=1.0d-7
+ xxsave=xx
+ xx=xx+aincr
+ write (2,*) xx,yy,zz
+ sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ de_dxx_num=(sumenep-sumene)/aincr
+ xx=xxsave
+ write (2,*) "xx+ sumene from enesc=",sumenep,sumene
+ yysave=yy
+ yy=yy+aincr
+ write (2,*) xx,yy,zz
+ sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ de_dyy_num=(sumenep-sumene)/aincr
+ yy=yysave
+ write (2,*) "yy+ sumene from enesc=",sumenep,sumene
+ zzsave=zz
+ zz=zz+aincr
+ write (2,*) xx,yy,zz
+ sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ de_dzz_num=(sumenep-sumene)/aincr
+ zz=zzsave
+ write (2,*) "zz+ sumene from enesc=",sumenep,sumene
+ costsave=cost2tab(i+1)
+ sintsave=sint2tab(i+1)
+ cost2tab(i+1)=dcos(0.5d0*(theta(i+1)+aincr))
+ sint2tab(i+1)=dsin(0.5d0*(theta(i+1)+aincr))
+ sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ de_dt_num=(sumenep-sumene)/aincr
+ write (2,*) " t+ sumene from enesc=",sumenep,sumene
+ cost2tab(i+1)=costsave
+ sint2tab(i+1)=sintsave
+!C End of diagnostics section.
+#endif
+!C
+!C Compute the gradient of esc
+!C
+ de_dxx=x(1)+2*x(4)*xx+x(7)*zz+x(8)*yy
+ de_dyy=x(2)+2*x(5)*yy+x(8)*xx+x(9)*zz
+ de_dzz=x(3)+2*x(6)*zz+x(7)*xx+x(9)*yy
+ de_dtt=0.0d0
+#ifdef DEBUG
+ write (2,*) "x",(x(k),k=1,9)
+ write (2,*) "xx",xx," yy",yy," zz",zz
+ write (2,*) "de_xx ",de_xx," de_yy ",de_yy,&
+ " de_zz ",de_zz," de_tt ",de_tt
+ write (2,*) "de_xx_num",de_dxx_num," de_yy_num",de_dyy_num,&
+ " de_zz_num",de_dzz_num," de_dt_num",de_dt_num
+#endif
+!C
+ cossc=scalar(dc_norm(1,i),dc_norm(1,i+nres))
+ cossc1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres))
+ cosfac2xx=cosfac2*xx
+ sinfac2yy=sinfac2*yy
+ do k = 1,3
+ dt_dCi(k) = -(dc_norm(k,i-1)+costtab(i+1)*dc_norm(k,i))*&
+ vbld_inv(i+1)
+ dt_dCi1(k)= -(dc_norm(k,i)+costtab(i+1)*dc_norm(k,i-1))*&
+ vbld_inv(i)
+ pom=(dC_norm(k,i+nres)-cossc*dC_norm(k,i))*vbld_inv(i+1)
+ pom1=(dC_norm(k,i+nres)-cossc1*dC_norm(k,i-1))*vbld_inv(i)
+!c write (iout,*) "i",i," k",k," pom",pom," pom1",pom1,
+!c & " dt_dCi",dt_dCi(k)," dt_dCi1",dt_dCi1(k)
+!c write (iout,*) "dC_norm",(dC_norm(j,i),j=1,3),
+!c & (dC_norm(j,i-1),j=1,3)," vbld_inv",vbld_inv(i+1),vbld_inv(i)
+ dXX_Ci(k)=pom*cosfac-dt_dCi(k)*cosfac2xx
+ dXX_Ci1(k)=-pom1*cosfac-dt_dCi1(k)*cosfac2xx
+ dYY_Ci(k)=pom*sinfac+dt_dCi(k)*sinfac2yy
+ dYY_Ci1(k)=pom1*sinfac+dt_dCi1(k)*sinfac2yy
+ dZZ_Ci1(k)=0.0d0
+ dZZ_Ci(k)=0.0d0
+ do j=1,3
+ dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres)
+ dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres)
+ enddo
+
+ dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
+ dYY_XYZ(k)=vbld_inv(i+nres)*(y_prime(k)-yy*dC_norm(k,i+nres))
+ dZZ_XYZ(k)=vbld_inv(i+nres)*(z_prime(k)-zz*dC_norm(k,i+nres))
+!c
+ dt_dCi(k) = -dt_dCi(k)/sinttab(i+1)
+ dt_dCi1(k)= -dt_dCi1(k)/sinttab(i+1)
+ enddo
+
+ do k=1,3
+ dXX_Ctab(k,i)=dXX_Ci(k)
+ dXX_C1tab(k,i)=dXX_Ci1(k)
+ dYY_Ctab(k,i)=dYY_Ci(k)
+ dYY_C1tab(k,i)=dYY_Ci1(k)
+ dZZ_Ctab(k,i)=dZZ_Ci(k)
+ dZZ_C1tab(k,i)=dZZ_Ci1(k)
+ dXX_XYZtab(k,i)=dXX_XYZ(k)
+ dYY_XYZtab(k,i)=dYY_XYZ(k)
+ dZZ_XYZtab(k,i)=dZZ_XYZ(k)
+ enddo
+ do k = 1,3
+!c write (iout,*) "k",k," dxx_ci1",dxx_ci1(k)," dyy_ci1",
+!c & dyy_ci1(k)," dzz_ci1",dzz_ci1(k)
+!c write (iout,*) "k",k," dxx_ci",dxx_ci(k)," dyy_ci",
+!c & dyy_ci(k)," dzz_ci",dzz_ci(k)
+!c write (iout,*) "k",k," dt_dci",dt_dci(k)," dt_dci",
+!c & dt_dci(k)
+!c write (iout,*) "k",k," dxx_XYZ",dxx_XYZ(k)," dyy_XYZ",
+!c & dyy_XYZ(k)," dzz_XYZ",dzz_XYZ(k)
+ gsbloc(k,i-1)=gsbloc(k,i-1)+de_dxx*dxx_ci1(k) &
+ +de_dyy*dyy_ci1(k)+de_dzz*dzz_ci1(k)+de_dt*dt_dCi1(k)
+ gsbloc(k,i)=gsbloc(k,i)+de_dxx*dxx_Ci(k) &
+ +de_dyy*dyy_Ci(k)+de_dzz*dzz_Ci(k)+de_dt*dt_dCi(k)
+ gsblocx(k,i)= de_dxx*dxx_XYZ(k)&
+ +de_dyy*dyy_XYZ(k)+de_dzz*dzz_XYZ(k)
+ enddo
+!c write(iout,*) "ENERGY GRAD = ", (gsbloc(k,i-1),k=1,3),
+!c & (gsbloc(k,i),k=1,3),(gsblocx(k,i),k=1,3)
+
+!C to check gradient call subroutine check_grad
+
+ 1 continue
+ enddo
+ return
+ end subroutine esb
+!=-------------------------------------------------------
+ real(kind=8) function enesc_nucl(x,xx,yy,zz,cost2,sint2)
+! implicit none
+ real(kind=8),dimension(9):: x(9)
+ real(kind=8) :: xx,yy,zz,cost2,sint2,sumene1,sumene2, &
+ sumene3,sumene4,sumene,dsc_i,dp2_i,dscp1,dscp2,s1,s1_6,s2,s2_6
+ integer i
+!c write (2,*) "enesc"
+!c write (2,*) "x",(x(i),i=1,9)
+!c write(2,*)"xx",xx," yy",yy," zz",zz," cost2",cost2," sint2",sint2
+ sumene=x(1)*xx+x(2)*yy+x(3)*zz+x(4)*xx**2 &
+ + x(5)*yy**2+x(6)*zz**2+x(7)*xx*zz+x(8)*xx*yy &
+ + x(9)*yy*zz
+ enesc_nucl=sumene
+ return
+ end function enesc_nucl
+!-----------------------------------------------------------------------------
+ subroutine multibody_hb_nucl(ecorr,ecorr3,n_corr,n_corr1)
+#ifdef MPI
+ include 'mpif.h'
+ integer,parameter :: max_cont=2000
+ integer,parameter:: max_dim=2*(8*3+6)
+ integer, parameter :: msglen1=max_cont*max_dim
+ integer,parameter :: msglen2=2*msglen1
+ integer source,CorrelType,CorrelID,Error
+ real(kind=8) :: buffer(max_cont,max_dim)
+ integer status(MPI_STATUS_SIZE)
+ integer :: ierror,nbytes
+#endif
+ real(kind=8),dimension(3):: gx(3),gx1(3)
+ real(kind=8) :: time00
+ logical lprn,ldone
+ integer i,j,i1,j1,jj,kk,num_conti,num_conti1,nn
+ real(kind=8) ecorr,ecorr3
+ integer :: n_corr,n_corr1,mm,msglen
+!C Set lprn=.true. for debugging
+ lprn=.false.
+ n_corr=0
+ n_corr1=0
+#ifdef MPI
+ if(.not.allocated(zapas2)) allocate(zapas2(3,maxconts,nres,8))
+
+ if (nfgtasks.le.1) goto 30
+ if (lprn) then
+ write (iout,'(a)') 'Contact function values:'
+ do i=nnt,nct-1
+ write (iout,'(2i3,50(1x,i2,f5.2))') &
+ i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), &
+ j=1,num_cont_hb(i))
+ enddo
+ endif
+!C Caution! Following code assumes that electrostatic interactions concerning
+!C a given atom are split among at most two processors!
+ CorrelType=477
+ CorrelID=fg_rank+1
+ ldone=.false.
+ do i=1,max_cont
+ do j=1,max_dim
+ buffer(i,j)=0.0D0
+ enddo
+ enddo
+ mm=mod(fg_rank,2)
+!c write (*,*) 'MyRank',MyRank,' mm',mm
+ if (mm) 20,20,10
+ 10 continue
+!c write (*,*) 'Sending: MyRank',MyRank,' mm',mm,' ldone',ldone
+ if (fg_rank.gt.0) then
+!C Send correlation contributions to the preceding processor
+ msglen=msglen1
+ nn=num_cont_hb(iatel_s_nucl)
+ call pack_buffer(max_cont,max_dim,iatel_s,0,buffer)
+!c write (*,*) 'The BUFFER array:'
+!c do i=1,nn
+!c write (*,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,30)
+!c enddo
+ if (ielstart_nucl(iatel_s_nucl).gt.iatel_s_nucl+ispp) then
+ msglen=msglen2
+ call pack_buffer(max_cont,max_dim,iatel_s+1,30,buffer)
+!C Clear the contacts of the atom passed to the neighboring processor
+ nn=num_cont_hb(iatel_s_nucl+1)
+!c do i=1,nn
+!c write (*,'(i2,9(3f8.3,2x))') i,(buffer(i,j+30),j=1,30)
+!c enddo
+ num_cont_hb(iatel_s_nucl)=0
+ endif
+!cd write (iout,*) 'Processor ',fg_rank,MyRank,
+!cd & ' is sending correlation contribution to processor',fg_rank-1,
+!cd & ' msglen=',msglen
+!c write (*,*) 'Processor ',fg_rank,MyRank,
+!c & ' is sending correlation contribution to processor',fg_rank-1,
+!c & ' msglen=',msglen,' CorrelType=',CorrelType
+ time00=MPI_Wtime()
+ call MPI_Send(buffer,msglen,MPI_DOUBLE_PRECISION,fg_rank-1, &
+ CorrelType,FG_COMM,IERROR)
+ time_sendrecv=time_sendrecv+MPI_Wtime()-time00
+!cd write (iout,*) 'Processor ',fg_rank,
+!cd & ' has sent correlation contribution to processor',fg_rank-1,
+!cd & ' msglen=',msglen,' CorrelID=',CorrelID
+!c write (*,*) 'Processor ',fg_rank,
+!c & ' has sent correlation contribution to processor',fg_rank-1,
+!c & ' msglen=',msglen,' CorrelID=',CorrelID
+!c msglen=msglen1
+ endif ! (fg_rank.gt.0)
+ if (ldone) goto 30
+ ldone=.true.
+ 20 continue
+!c write (*,*) 'Receiving: MyRank',MyRank,' mm',mm,' ldone',ldone
+ if (fg_rank.lt.nfgtasks-1) then
+!C Receive correlation contributions from the next processor
+ msglen=msglen1
+ if (ielend_nucl(iatel_e_nucl).lt.nct_molec(2)-1) msglen=msglen2
+!cd write (iout,*) 'Processor',fg_rank,
+!cd & ' is receiving correlation contribution from processor',fg_rank+1,
+!cd & ' msglen=',msglen,' CorrelType=',CorrelType
+!c write (*,*) 'Processor',fg_rank,
+!c &' is receiving correlation contribution from processor',fg_rank+1,
+!c & ' msglen=',msglen,' CorrelType=',CorrelType
+ time00=MPI_Wtime()
+ nbytes=-1
+ do while (nbytes.le.0)
+ call MPI_Probe(fg_rank+1,CorrelType,FG_COMM,status,IERROR)
+ call MPI_Get_count(status,MPI_DOUBLE_PRECISION,nbytes,IERROR)
+ enddo
+!c print *,'Processor',myrank,' msglen',msglen,' nbytes',nbytes
+ call MPI_Recv(buffer,nbytes,MPI_DOUBLE_PRECISION, &
+ fg_rank+1,CorrelType,FG_COMM,status,IERROR)
+ time_sendrecv=time_sendrecv+MPI_Wtime()-time00
+!c write (*,*) 'Processor',fg_rank,
+!c &' has received correlation contribution from processor',fg_rank+1,
+!c & ' msglen=',msglen,' nbytes=',nbytes
+!c write (*,*) 'The received BUFFER array:'
+!c do i=1,max_cont
+!c write (*,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,60)
+!c enddo
+ if (msglen.eq.msglen1) then
+ call unpack_buffer(max_cont,max_dim,iatel_e_nucl+1,0,buffer)
+ else if (msglen.eq.msglen2) then
+ call unpack_buffer(max_cont,max_dim,iatel_e_nucl,0,buffer)
+ call unpack_buffer(max_cont,max_dim,iatel_e_nucl+1,30,buffer)
+ else
+ write (iout,*) &
+ 'ERROR!!!! message length changed while processing correlations.'
+ write (*,*) &
+ 'ERROR!!!! message length changed while processing correlations.'
+ call MPI_Abort(MPI_COMM_WORLD,Error,IERROR)
+ endif ! msglen.eq.msglen1
+ endif ! fg_rank.lt.nfgtasks-1
+ if (ldone) goto 30
+ ldone=.true.
+ goto 10
+ 30 continue
+#endif
+ if (lprn) then
+ write (iout,'(a)') 'Contact function values:'
+ do i=nnt_molec(2),nct_molec(2)-1
+ write (iout,'(2i3,50(1x,i2,f5.2))') &
+ i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), &
+ j=1,num_cont_hb(i))
+ enddo
+ endif
+ ecorr=0.0D0
+ ecorr3=0.0d0
+!C Remove the loop below after debugging !!!
+ do i=nnt_molec(2),nct_molec(2)
+ do j=1,3
+ gradcorr_nucl(j,i)=0.0D0
+ gradxorr_nucl(j,i)=0.0D0
+ gradcorr3_nucl(j,i)=0.0D0
+ gradxorr3_nucl(j,i)=0.0D0
+ enddo
+ enddo
+ print *,"iatsc_s_nucl,iatsc_e_nucl",iatsc_s_nucl,iatsc_e_nucl
+!C Calculate the local-electrostatic correlation terms
+ do i=iatsc_s_nucl,iatsc_e_nucl
+ i1=i+1
+ num_conti=num_cont_hb(i)
+ num_conti1=num_cont_hb(i+1)
+ print *,i,num_conti,num_conti1
+ do jj=1,num_conti
+ j=jcont_hb(jj,i)
+ do kk=1,num_conti1
+ j1=jcont_hb(kk,i1)
+!c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
+!c & ' jj=',jj,' kk=',kk
+ if (j1.eq.j+1 .or. j1.eq.j-1) then
+!C
+!C Contacts I-J and (I+1)-(J+1) or (I+1)-(J-1) occur simultaneously.
+!C The system gains extra energy.
+!C Tentative expression & coefficients; assumed d(stacking)=4.5 A,
+!C parallel dipoles of stacknig bases and sin(mui)sin(muj)/eps/d^3=0.7
+!C Need to implement full formulas 34 and 35 from Liwo et al., 1998.
+!C
+ ecorr=ecorr+ehbcorr_nucl(i,j,i+1,j1,jj,kk,0.528D0,0.132D0)
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
+ 'ecorrh',i,j,ehbcorr_nucl(i,j,i+1,j1,jj,kk,0.528D0,0.132D0)
+ n_corr=n_corr+1
+ else if (j1.eq.j) then
+!C
+!C Contacts I-J and I-(J+1) occur simultaneously.
+!C The system loses extra energy.
+!C Tentative expression & c?oefficients; assumed d(stacking)=4.5 A,
+!C parallel dipoles of stacknig bases and sin(mui)sin(muj)/eps/d^3=0.7
+!C Need to implement full formulas 32 from Liwo et al., 1998.
+!C
+!c write (iout,*) 'ecorr3: i=',i,' j=',j,' i1=',i1,' j1=',j1,
+!c & ' jj=',jj,' kk=',kk
+ ecorr3=ecorr3+ehbcorr3_nucl(i,j,i+1,j,jj,kk,0.310D0,-0.155D0)
+ endif
+ enddo ! kk
+ do kk=1,num_conti
+ j1=jcont_hb(kk,i)
+!c write (iout,*) 'ecorr3: i=',i,' j=',j,' i1=',i1,' j1=',j1,
+!c & ' jj=',jj,' kk=',kk
+ if (j1.eq.j+1) then
+!C Contacts I-J and (I+1)-J occur simultaneously.
+!C The system loses extra energy.
+ ecorr3=ecorr3+ehbcorr3_nucl(i,j,i,j+1,jj,kk,0.310D0,-0.155D0)
+ endif ! j1==j+1
+ enddo ! kk
+ enddo ! jj
+ enddo ! i
+ return
+ end subroutine multibody_hb_nucl
+!-----------------------------------------------------------
+ real(kind=8) function ehbcorr_nucl(i,j,k,l,jj,kk,coeffp,coeffm)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.DERIV'
+! include 'COMMON.INTERACT'
+! include 'COMMON.CONTACTS'
+ real(kind=8),dimension(3) :: gx,gx1
+ logical :: lprn
+!el local variables
+ 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, &
+ rlocshield
+
+ lprn=.false.
+ eij=facont_hb(jj,i)
+ ekl=facont_hb(kk,k)
+ ees0pij=ees0p(jj,i)
+ ees0pkl=ees0p(kk,k)
+ ees0mij=ees0m(jj,i)
+ ees0mkl=ees0m(kk,k)
+ ekont=eij*ekl
+ ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl)
+! print *,"ehbcorr_nucl",ekont,ees
+!cd ees=-(coeffp*ees0pkl+coeffm*ees0mkl)
+!C Following 4 lines for diagnostics.
+!cd ees0pkl=0.0D0
+!cd ees0pij=1.0D0
+!cd ees0mkl=0.0D0
+!cd ees0mij=1.0D0
+!cd write (iout,*)'Contacts have occurred for nucleic bases',
+!cd & i,j,' fcont:',eij,' eij',' eesij',ees0pij,ees0mij,' and ',k,l
+!cd & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' ees=',ees
+!C Calculate the multi-body contribution to energy.
+! ecorr_nucl=ecorr_nucl+ekont*ees
+!C Calculate multi-body contributions to the gradient.
+ coeffpees0pij=coeffp*ees0pij
+ coeffmees0mij=coeffm*ees0mij
+ coeffpees0pkl=coeffp*ees0pkl
+ coeffmees0mkl=coeffm*ees0mkl
+ do ll=1,3
+ gradxorr_nucl(ll,i)=gradxorr_nucl(ll,i) &
+ -ekont*(coeffpees0pkl*gacontp_hb1(ll,jj,i)+&
+ coeffmees0mkl*gacontm_hb1(ll,jj,i))
+ gradxorr_nucl(ll,j)=gradxorr_nucl(ll,j) &
+ -ekont*(coeffpees0pkl*gacontp_hb2(ll,jj,i)+&
+ coeffmees0mkl*gacontm_hb2(ll,jj,i))
+ gradxorr_nucl(ll,k)=gradxorr_nucl(ll,k) &
+ -ekont*(coeffpees0pij*gacontp_hb1(ll,kk,k)+&
+ coeffmees0mij*gacontm_hb1(ll,kk,k))
+ gradxorr_nucl(ll,l)=gradxorr_nucl(ll,l) &
+ -ekont*(coeffpees0pij*gacontp_hb2(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb2(ll,kk,k))
+ gradlongij=ees*ekl*gacont_hbr(ll,jj,i)- &
+ ekont*(coeffpees0pkl*gacontp_hb3(ll,jj,i)+ &
+ coeffmees0mkl*gacontm_hb3(ll,jj,i))
+ gradcorr_nucl(ll,j)=gradcorr_nucl(ll,j)+gradlongij
+ gradcorr_nucl(ll,i)=gradcorr_nucl(ll,i)-gradlongij
+ gradlongkl=ees*eij*gacont_hbr(ll,kk,k)- &
+ ekont*(coeffpees0pij*gacontp_hb3(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb3(ll,kk,k))
+ gradcorr_nucl(ll,l)=gradcorr_nucl(ll,l)+gradlongkl
+ gradcorr_nucl(ll,k)=gradcorr_nucl(ll,k)-gradlongkl
+ gradxorr_nucl(ll,i)=gradxorr_nucl(ll,i)-gradlongij
+ gradxorr_nucl(ll,j)=gradxorr_nucl(ll,j)+gradlongij
+ gradxorr_nucl(ll,k)=gradxorr_nucl(ll,k)-gradlongkl
+ gradxorr_nucl(ll,l)=gradxorr_nucl(ll,l)+gradlongkl
+ enddo
+ ehbcorr_nucl=ekont*ees
+ return
+ end function ehbcorr_nucl
+!-------------------------------------------------------------------------
+
+ real(kind=8) function ehbcorr3_nucl(i,j,k,l,jj,kk,coeffp,coeffm)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.DERIV'
+! include 'COMMON.INTERACT'
+! include 'COMMON.CONTACTS'
+ real(kind=8),dimension(3) :: gx,gx1
+ logical :: lprn
+!el local variables
+ 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, &
+ rlocshield
+
+ lprn=.false.
+ eij=facont_hb(jj,i)
+ ekl=facont_hb(kk,k)
+ ees0pij=ees0p(jj,i)
+ ees0pkl=ees0p(kk,k)
+ ees0mij=ees0m(jj,i)
+ ees0mkl=ees0m(kk,k)
+ ekont=eij*ekl
+ ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl)
+!cd ees=-(coeffp*ees0pkl+coeffm*ees0mkl)
+!C Following 4 lines for diagnostics.
+!cd ees0pkl=0.0D0
+!cd ees0pij=1.0D0
+!cd ees0mkl=0.0D0
+!cd ees0mij=1.0D0
+!cd write (iout,*)'Contacts have occurred for nucleic bases',
+!cd & i,j,' fcont:',eij,' eij',' eesij',ees0pij,ees0mij,' and ',k,l
+!cd & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' ees=',ees
+!C Calculate the multi-body contribution to energy.
+! ecorr=ecorr+ekont*ees
+!C Calculate multi-body contributions to the gradient.
+ coeffpees0pij=coeffp*ees0pij
+ coeffmees0mij=coeffm*ees0mij
+ coeffpees0pkl=coeffp*ees0pkl
+ coeffmees0mkl=coeffm*ees0mkl
+ do ll=1,3
+ gradxorr3_nucl(ll,i)=gradxorr3_nucl(ll,i) &
+ -ekont*(coeffpees0pkl*gacontp_hb1(ll,jj,i)+&
+ coeffmees0mkl*gacontm_hb1(ll,jj,i))
+ gradxorr3_nucl(ll,j)=gradxorr3_nucl(ll,j) &
+ -ekont*(coeffpees0pkl*gacontp_hb2(ll,jj,i)+ &
+ coeffmees0mkl*gacontm_hb2(ll,jj,i))
+ gradxorr3_nucl(ll,k)=gradxorr3_nucl(ll,k) &
+ -ekont*(coeffpees0pij*gacontp_hb1(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb1(ll,kk,k))
+ gradxorr3_nucl(ll,l)=gradxorr3_nucl(ll,l) &
+ -ekont*(coeffpees0pij*gacontp_hb2(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb2(ll,kk,k))
+ gradlongij=ees*ekl*gacont_hbr(ll,jj,i)- &
+ ekont*(coeffpees0pkl*gacontp_hb3(ll,jj,i)+ &
+ coeffmees0mkl*gacontm_hb3(ll,jj,i))
+ gradcorr3_nucl(ll,j)=gradcorr3_nucl(ll,j)+gradlongij
+ gradcorr3_nucl(ll,i)=gradcorr3_nucl(ll,i)-gradlongij
+ gradlongkl=ees*eij*gacont_hbr(ll,kk,k)- &
+ ekont*(coeffpees0pij*gacontp_hb3(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb3(ll,kk,k))
+ gradcorr3_nucl(ll,l)=gradcorr3_nucl(ll,l)+gradlongkl
+ gradcorr3_nucl(ll,k)=gradcorr3_nucl(ll,k)-gradlongkl
+ gradxorr3_nucl(ll,i)=gradxorr3_nucl(ll,i)-gradlongij
+ gradxorr3_nucl(ll,j)=gradxorr3_nucl(ll,j)+gradlongij
+ gradxorr3_nucl(ll,k)=gradxorr3_nucl(ll,k)-gradlongkl
+ gradxorr3_nucl(ll,l)=gradxorr3_nucl(ll,l)+gradlongkl
+ enddo
+ ehbcorr3_nucl=ekont*ees
+ return
+ end function ehbcorr3_nucl
+#ifdef MPI
+ subroutine pack_buffer(dimen1,dimen2,atom,indx,buffer)
+ integer dimen1,dimen2,atom,indx,numcont,i,ii,k,j,num_kont,num_kont_old
+ real(kind=8):: buffer(dimen1,dimen2)
+ num_kont=num_cont_hb(atom)
+ do i=1,num_kont
+ do k=1,8
+ do j=1,3
+ buffer(i,indx+(k-1)*3+j)=zapas2(j,i,atom,k)
+ enddo ! j
+ enddo ! k
+ buffer(i,indx+25)=facont_hb(i,atom)
+ buffer(i,indx+26)=ees0p(i,atom)
+ buffer(i,indx+27)=ees0m(i,atom)
+ buffer(i,indx+28)=d_cont(i,atom)
+ buffer(i,indx+29)=dfloat(jcont_hb(i,atom))
+ enddo ! i
+ buffer(1,indx+30)=dfloat(num_kont)
+ return
+ end subroutine pack_buffer
+!c------------------------------------------------------------------------------
+ subroutine unpack_buffer(dimen1,dimen2,atom,indx,buffer)
+ integer dimen1,dimen2,atom,indx,numcont,i,ii,k,j,num_kont,num_kont_old
+ real(kind=8):: buffer(dimen1,dimen2)
+! double precision zapas
+! common /contacts_hb/ zapas(3,maxconts,maxres,8),
+! & facont_hb(maxconts,maxres),ees0p(maxconts,maxres),
+! & ees0m(maxconts,maxres),d_cont(maxconts,maxres),
+! & num_cont_hb(maxres),jcont_hb(maxconts,maxres)
+ num_kont=buffer(1,indx+30)
+ num_kont_old=num_cont_hb(atom)
+ num_cont_hb(atom)=num_kont+num_kont_old
+ do i=1,num_kont
+ ii=i+num_kont_old
+ do k=1,8
+ do j=1,3
+ zapas2(j,ii,atom,k)=buffer(i,indx+(k-1)*3+j)
+ enddo ! j
+ enddo ! k
+ facont_hb(ii,atom)=buffer(i,indx+25)
+ ees0p(ii,atom)=buffer(i,indx+26)
+ ees0m(ii,atom)=buffer(i,indx+27)
+ d_cont(i,atom)=buffer(i,indx+28)
+ jcont_hb(ii,atom)=buffer(i,indx+29)
+ enddo ! i
+ return
+ end subroutine unpack_buffer
+!c------------------------------------------------------------------------------
+#endif
+
+!----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
end module energy