working torsional

diff --git a/source/unres/control.F90 b/source/unres/control.F90
index a45ba29..86f23f9 100644 (file)
--- a/source/unres/control.F90
+++ b/source/unres/control.F90
@@ -832,6 +832,10 @@
       iphi_end=iturn3_end+2
       iturn3_start=iturn3_start-1
       iturn3_end=iturn3_end-1
+      call int_bounds(nct_molec(2)-nnt_molec(2)-2,iphi_nucl_start,iphi_nucl_end)
+      iphi_nucl_start=iphi_nucl_start+nnt_molec(2)+2
+      iphi_nucl_end=iphi_nucl_end+nnt_molec(2)+2
+      print *,"KURDE",iphi_nucl_start,iphi_nucl_end
       call int_bounds(nres_molec(1)-3,itau_start,itau_end)
       itau_start=itau_start+3
       itau_end=itau_end+3
@@ -1394,6 +1398,8 @@
       iphi_end=nct_molec(1)
       iphi1_start=4
       iphi1_end=nres_molec(1)
+      iphi_nucl_start=4+nres_molec(1)
+      iphi_nucl_end=nres_molec(1)+nres_molec(2)
       idihconstr_start=1
       idihconstr_end=ndih_constr
       ithetaconstr_start=1

diff --git a/source/unres/data/energy_data.f90 b/source/unres/data/energy_data.f90
index c41fc11..844aa85 100644 (file)
--- a/source/unres/data/energy_data.f90
+++ b/source/unres/data/energy_data.f90
@@ -65,7 +65,8 @@
       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
@@ -179,7 +180,8 @@
        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,&
@@ -261,6 +263,15 @@
       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 

diff --git a/source/unres/energy.f90 b/source/unres/energy.f90
index 2fc2773..fb4caac 100644 (file)
--- a/source/unres/energy.f90
+++ b/source/unres/energy.f90
@@ -548,6 +548,7 @@
 !--------------------------------------------------------
       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
@@ -718,8 +719,20 @@
       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 &
@@ -729,7 +742,10 @@
        +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 &
@@ -738,7 +754,10 @@
        +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
@@ -899,8 +918,19 @@
       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,&
@@ -912,6 +942,10 @@
         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)'/ &
@@ -942,6 +976,17 @@
        '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,&
@@ -953,6 +998,10 @@
         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)'/ &
@@ -982,6 +1031,17 @@
        '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
@@ -19960,6 +20020,86 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
       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
 
 !-----------------------------------------------------------------------------
 !-----------------------------------------------------------------------------

diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90
index 25a5e24..4c7f3a5 100644 (file)
--- a/source/unres/io_config.f90
+++ b/source/unres/io_config.f90
@@ -1681,6 +1681,48 @@
       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