From e1aeda758c35ccdf5a5b88e825b6ca391c2dc1ae Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Sun, 6 Aug 2017 18:56:28 +0200 Subject: [PATCH] working torsional --- source/unres/control.F90 | 6 ++ source/unres/data/energy_data.f90 | 15 +++- source/unres/energy.f90 | 144 ++++++++++++++++++++++++++++++++++++- source/unres/io_config.f90 | 42 +++++++++++ 4 files changed, 203 insertions(+), 4 deletions(-) diff --git a/source/unres/control.F90 b/source/unres/control.F90 index a45ba29..86f23f9 100644 --- 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 --- 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 --- 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 --- 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 -- 1.7.9.5