real(kind=8) :: wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,wang,&
wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,&
wturn6,wvdwpp,wliptran,wshield,lipscale,wtube, &
- wbond_nucl,wang_nucl
+ wbond_nucl,wang_nucl,wcorr_nucl,wcorr3_nucl,welpp,wtor_nucl,&
+ wtor_d_nucl,welsb,wsbloc,wvdwsb,welpsb,wvdwpsb
#ifdef CLUSTER
real(kind=8) :: scalscp
#endif
iturn3_start,iturn3_end,iturn4_start,iturn4_end,iint_start,&
iint_end,iphi1_start,iphi1_end,itau_start,itau_end,&
ilip_start,ilip_end,itube_start,itube_end
- integer :: ibond_nucl_start,ibond_nucl_end, &
+ integer :: ibond_nucl_start,ibond_nucl_end,iphi_nucl_start,&
+ iphi_nucl_end,iphid_nucl_start,iphid_nucl_end,&
ibondp_nucl_start,ibondp_nucl_end,ithet_nucl_start,ithet_nucl_end
integer,dimension(:),allocatable :: ibond_displ,ibond_count,&
ithet_displ,ithet_count,iphi_displ,iphi_count,iphi1_displ,&
integer,dimension(:),allocatable :: itortyp !(-ntyp1:ntyp1)
integer,dimension(:,:,:),allocatable :: nterm,nlor !(-maxtor:maxtor,-maxtor:maxtor,2)
integer :: ntortyp,nterm_old
+!------torsion nucleic
+ real(kind=8),dimension(:,:),allocatable :: v0_nucl !(-maxtor:maxtor,-maxtor:maxtor,2)
+ real(kind=8),dimension(:,:,:),allocatable :: v1_nucl,v2_nucl !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
+ real(kind=8),dimension(:,:,:),allocatable :: vlor1_nucl !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
+ real(kind=8),dimension(:,:,:),allocatable :: vlor2_nucl,vlor3_nucl !(maxlor,maxtor,maxtor)
+ integer,dimension(:),allocatable :: itortyp_nucl !(-ntyp1:ntyp1)
+ integer,dimension(:,:),allocatable :: nterm_nucl,nlor_nucl !(-maxtor:maxtor,-maxtor:maxtor,2)
+ integer :: ntortyp_nucl,nterm_old_nucl
+
! 6/23/01 - constants for double torsionals
! common /torsiond/
real(kind=8),dimension(:,:,:,:,:,:),allocatable :: v1c,v1s
!--------------------------------------------------------
call ebond_nucl(estr_nucl)
call ebend_nucl(ebe_nucl)
+ call etor_nucl(etors_nucl)
print *,"after ebend", ebe_nucl
#ifdef TIMING
time_enecalc=time_enecalc+MPI_Wtime()-time00
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 &
+wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d &
+wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube&
+Eafmforce+ethetacnstr &
- +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl
+ +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 &
+wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d &
+wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube&
+Eafmforce+ethetacnstr &
- +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl
+ +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
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,&
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)'/ &
'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,&
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)'/ &
'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
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
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
enddo
endif
#endif
+ allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
+ read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
+ print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
+!el from energy module---------
+ allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
+ allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
+
+ allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
+ allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
+ allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
+ allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
+
+ allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
+ allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
+!el---------------------------
+ nterm_nucl(:,:)=0
+ nlor_nucl(:,:)=0
+!el--------------------
+ read (itorp_nucl,*,end=113,err=113) &
+ (itortyp_nucl(i),i=1,ntyp_molec(2))
+ print *,itortyp_nucl(:)
+!c write (iout,*) 'ntortyp',ntortyp
+ do i=1,ntortyp_nucl
+ do j=1,ntortyp_nucl
+ read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
+ print *,nterm_nucl(i,j),nlor_nucl(i,j)
+ v0ij=0.0d0
+ si=-1.0d0
+ do k=1,nterm_nucl(i,j)
+ read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
+ v0ij=v0ij+si*v1_nucl(k,i,j)
+ si=-si
+ enddo
+ do k=1,nlor_nucl(i,j)
+ read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
+ vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
+ v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
+ enddo
+ v0_nucl(i,j)=v0ij
+ enddo
+ enddo
+
! Read of Side-chain backbone correlation parameters
! Modified 11 May 2012 by Adasko
!CC