X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.f90;h=2fc27730ff6238fb5320ac2abb2e3d1ab6f61fed;hb=bbbdc8e18680625d3004f414aec255e9ca7b7353;hp=03a0e344471ecd287e885cda016c5163715dab63;hpb=4367d241fbb2bc284580092d2d177b7c79ac3a42;p=unres4.git diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index 03a0e34..2fc2773 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -45,6 +45,7 @@ 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/ @@ -124,6 +125,8 @@ gshieldc_t3,gshieldc_loc_t3,gshieldx_t4,gshieldc_t4, & gshieldc_loc_t4,gshieldx_ll,gshieldc_ll,gshieldc_loc_ll,& grad_shield,gg_tube,gg_tube_sc,gradafm !(3,maxres) +!-----------------------------NUCLEIC GRADIENT + real(kind=8),dimension(:,:),allocatable ::gradb_nucl,gradbx_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) @@ -225,7 +228,10 @@ 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 @@ -356,6 +362,7 @@ 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 @@ -421,6 +428,7 @@ ! Calculate the bond-stretching energy ! call ebond(estr) + print *,"EBOND",estr ! write(iout,*) "in etotal afer ebond",ipot ! @@ -537,7 +545,10 @@ else etube=0.0d0 endif - +!-------------------------------------------------------- + call ebond_nucl(estr_nucl) + call ebend_nucl(ebe_nucl) + print *,"after ebend", ebe_nucl #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 #endif @@ -583,6 +594,21 @@ 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" @@ -625,6 +651,10 @@ 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, 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 @@ -688,6 +718,9 @@ Eafmforce=energia(23) ethetacnstr=energia(24) etube=energia(25) + estr_nucl=energia(32) + ebe_nucl=energia(33) + #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc & @@ -695,7 +728,8 @@ +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& - +Eafmforce+ethetacnstr + +Eafmforce+ethetacnstr & + +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc & @@ -703,8 +737,8 @@ +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& - +Eafmforce+ethetacnstr - + +Eafmforce+ethetacnstr & + +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl #endif energia(0)=etot ! detecting NaNQ @@ -829,6 +863,9 @@ real(kind=8) :: eello_turn6,eello_turn3,eello_turn4,ebe,escloc real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,Uconst,esccor,eliptran,& 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) @@ -862,6 +899,9 @@ Eafmforce=energia(23) ethetacnstr=energia(24) etube=energia(25) + estr_nucl=energia(32) + ebe_nucl=energia(33) + #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,& estr,wbond,ebe,wang,& @@ -870,7 +910,9 @@ ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,& edihcnstr,ethetacnstr,ebr*nss,& - Uconst,eliptran,wliptran,Eafmforce,etube,wtube,etot + Uconst,eliptran,wliptran,Eafmforce,etube,wtube, & ! till now protein + estr_nucl,wbond_nucl,ebe_nucl,wang_nucl, & + etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ & @@ -898,6 +940,8 @@ '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)'/ & 'ETOT= ',1pE16.6,' (total)') #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,& @@ -907,7 +951,9 @@ ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,& ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, & - etube,wtube,etot + etube,wtube, & + estr_nucl,wbond_nucl, ebe_nucl,wang_nucl,& + etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ & @@ -934,6 +980,8 @@ '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)'/ & 'ETOT= ',1pE16.6,' (total)') #endif return @@ -1006,7 +1054,7 @@ !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 @@ -1161,7 +1209,7 @@ !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) @@ -1302,7 +1350,7 @@ 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), @@ -1573,7 +1621,7 @@ 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 @@ -1596,7 +1644,7 @@ 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, & @@ -1608,6 +1656,7 @@ !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 @@ -1640,6 +1689,7 @@ enddo ! j enddo ! iint enddo ! i +! print *,"ZALAMKA", evdw ! write (iout,*) "Number of loop steps in EGB:",ind !ccc energy_dec=.false. return @@ -1757,7 +1807,7 @@ 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,& @@ -5143,11 +5193,13 @@ ! 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,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 @@ -5156,6 +5208,7 @@ "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 @@ -5184,6 +5237,11 @@ 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 @@ -6514,7 +6572,7 @@ 'etor',i,etors_ii if (lprn) & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') & - restyp(itype(i-2,1)),i-2,restyp(itype(i-1,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) @@ -6623,7 +6681,7 @@ '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,1)),i-2,restyp(itype(i-1,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 @@ -6818,7 +6876,7 @@ 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,1)),i-2,restyp(itype(i-1,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 @@ -10280,10 +10338,8 @@ +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 @@ -10305,9 +10361,8 @@ +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 @@ -10464,7 +10519,9 @@ +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 @@ -10498,7 +10555,9 @@ +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) + @@ -10514,7 +10573,9 @@ +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 @@ -10724,7 +10785,7 @@ ! 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 & @@ -12399,7 +12460,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !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) @@ -12486,7 +12547,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !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) @@ -12612,7 +12673,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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), @@ -12732,7 +12793,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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), @@ -12982,7 +13043,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -13003,7 +13064,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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,& @@ -13273,7 +13334,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -13294,7 +13355,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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,& @@ -13437,7 +13498,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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,& @@ -13566,7 +13627,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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,& @@ -16096,6 +16157,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -18327,7 +18390,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !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 @@ -18488,7 +18551,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !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 @@ -18725,12 +18788,12 @@ write(iout,*) 'Calling CHECK_ECARTIN 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" @@ -18826,7 +18889,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !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 @@ -18929,6 +18992,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -18937,6 +19001,23 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -19424,6 +19505,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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)) !(3,maxres) allocate(grad_shield_side(3,50,nres)) allocate(grad_shield_loc(3,50,nres)) @@ -19552,6 +19635,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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. !---------------------- @@ -19598,6 +19684,283 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 + !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- end module energy