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)
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
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
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
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"
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
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 &
+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 &
+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
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)
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,&
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)'/ &
'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,&
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)'/ &
'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
!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
enddo ! j
enddo ! iint
enddo ! i
+! print *,"ZALAMKA", evdw
! write (iout,*) "Number of loop steps in EGB:",ind
!ccc energy_dec=.false.
return
! 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
!
"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
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
+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
+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
+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
+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)
+
+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
! 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 &
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
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"
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
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))
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
projects / unres4.git / blobdiff