X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.f90;h=2fc27730ff6238fb5320ac2abb2e3d1ab6f61fed;hb=bbbdc8e18680625d3004f414aec255e9ca7b7353;hp=3975993d0115b0bd876d9563fd3cac4ae59a9f11;hpb=7c0faf2ccb6f94ed9a77aa527d3885ba054d3fb2;p=unres4.git diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index 3975993..2fc2773 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -125,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) @@ -226,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 @@ -357,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 @@ -539,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 @@ -585,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" @@ -627,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 @@ -690,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 & @@ -697,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 & @@ -705,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 @@ -831,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) @@ -864,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,& @@ -872,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)'/ & @@ -900,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,& @@ -909,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)'/ & @@ -936,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 @@ -1610,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 @@ -1642,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 @@ -5145,7 +5193,7 @@ ! endif enddo estr=0.5d0*AKP*estr+estr1 - print *,"estr_bb",estr,AKP +! print *,"estr_bb",estr,AKP ! ! 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included ! @@ -5160,7 +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 +! print *,"estr_sc",estr do j=1,3 gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres) enddo @@ -5189,11 +5237,11 @@ usumsqder=usumsqder+ud(j)*uprod2 enddo estr=estr+uprod/usum - print *,"estr_sc",estr,i +! print *,"estr_sc",estr,i if (energy_dec) write (iout,*) & "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,& - AKSC(1,iti),AKSC(1,iti)*diff*diff + AKSC(1,iti),uprod/usum do j=1,3 gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres) enddo @@ -10290,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 @@ -10315,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 @@ -10474,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 @@ -10508,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) + @@ -10524,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 @@ -10734,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 & @@ -16106,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 @@ -18740,7 +18793,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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" @@ -18948,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 @@ -19435,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)) @@ -19612,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