- module energy
+ module energy
!-----------------------------------------------------------------------------
use io_units
use names
wscbase=weights(46)
wscpho=weights(47)
wpeppho=weights(48)
+! welpsb=weights(28)*fact(1)
+!
+! wcorr_nucl= weights(37)*fact(1)
+! wcorr3_nucl=weights(38)*fact(2)
+! wtor_nucl= weights(35)*fact(1)
+! wtor_d_nucl=weights(36)*fact(2)
+
endif
time_Bcast=time_Bcast+MPI_Wtime()-time00
time_Bcastw=time_Bcastw+MPI_Wtime()-time00
call AFMforce(Eafmforce)
else if (selfguide.gt.0) then
call AFMvel(Eafmforce)
+ else
+ Eafmforce=0.0d0
endif
endif
if (tubemode.eq.1) then
eespp=0.0d0
endif
! write(iout,*) ecorr_nucl,"ecorr_nucl",nres_molec(2)
-! print *,"before ecatcat"
+! print *,"before ecatcat",wcatcat
if (nfgtasks.gt.1) then
if (fg_rank.eq.0) then
call ecatcat(ecationcation)
epeppho=0.0
endif
! call ecatcat(ecationcation)
-! print *,"after ebend", ebe_nucl
+! print *,"after ebend", wtor_nucl
#ifdef TIMING
time_enecalc=time_enecalc+MPI_Wtime()-time00
#endif
! 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"
- energia(41)=ecation_prot
- energia(42)=ecationcation
+ energia(42)=ecation_prot
+ energia(41)=ecationcation
energia(46)=escbase
energia(47)=epepbase
energia(48)=escpho
etors_d_nucl=energia(36)
ecorr_nucl=energia(37)
ecorr3_nucl=energia(38)
- ecation_prot=energia(41)
- ecationcation=energia(42)
+ ecation_prot=energia(42)
+ ecationcation=energia(41)
escbase=energia(46)
epepbase=energia(47)
escpho=energia(48)
wtor=weights(13)*fact(1)
wtor_d=weights(14)*fact(2)
wsccor=weights(21)*fact(1)
+ welpsb=weights(28)*fact(1)
+ wcorr_nucl= weights(37)*fact(1)
+ wcorr3_nucl=weights(38)*fact(2)
+ wtor_nucl= weights(35)*fact(1)
+ wtor_d_nucl=weights(36)*fact(2)
return
end subroutine rescale_weights
etors_d_nucl=energia(36)
ecorr_nucl=energia(37)
ecorr3_nucl=energia(38)
- ecation_prot=energia(41)
- ecationcation=energia(42)
+ ecation_prot=energia(42)
+ ecationcation=energia(41)
escbase=energia(46)
epepbase=energia(47)
escpho=energia(48)
ecorr,wcorr,&
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, &
+ ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforce, &
etube,wtube, &
estr_nucl,wbond_nucl, ebe_nucl,wang_nucl,&
- evdwpp,wvdwpp_nucl,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb&
- evdwsb,wvdwsb,eelsb,welsb,esbloc,wsbloc,etors_nucl,wtor_nucl&
+ evdwpp,wvdwpp_nucl,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,ecation_prot,wcatprot,ecationcation,wcatcat, &
escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,&
! write(iout,*)"c ", c(1,:), c(2,:), c(3,:)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
- sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
- sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ sss_ele_cut=sscale_ele(1.0d0/(rij))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij))
! print *,sss_ele_cut,sss_ele_grad,&
! 1.0d0/(rij),r_cut_ele,rlamb_ele
if (sss_ele_cut.le.0.0) cycle
fac=rij*fac
! print *,'before fac',fac,rij,evdwij
fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
- /sigma(itypi,itypj)*rij
+ *rij
! print *,'grad part scale',fac, &
! evdwij*sss_ele_grad/sss_ele_cut &
! /sigma(itypi,itypj)*rij
#endif
#else
if (i.gt. nnt+2 .and. i.lt.nct+2) then
+! write(iout,*) "i,",molnum(i)
+! print *, "i,",molnum(i),i,itype(i-2,1)
+ if (molnum(i).eq.1) then
iti = itype2loc(itype(i-2,1))
else
iti=nloctyp
endif
+ else
+ iti=nloctyp
+ endif
!c write (iout,*) "i",i-1," itype",itype(i-2)," iti",iti
!c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
if (i.gt. nnt+1 .and. i.lt.nct+1) then
! if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
if (i.gt. nnt+1 .and. i.lt.nct+1) then
if (itype(i-1,1).eq.0) then
- iti1=ntortyp+1
+ iti1=nloctyp
elseif (itype(i-1,1).le.ntyp) then
iti1 = itype2loc(itype(i-1,1))
else
+wcorr3_nucl*gradcorr3_nucl(j,i) +&
wcatprot* gradpepcat(j,i)+ &
wcatcat*gradcatcat(j,i)+ &
- wscbase*gvdwc_scbase(j,i) &
+ wscbase*gvdwc_scbase(j,i)+ &
wpepbase*gvdwc_pepbase(j,i)+&
wscpho*gvdwc_scpho(j,i)+&
wpeppho*gvdwc_peppho(j,i)
+gradafm(j,i) &
+wliptran*gliptranc(j,i) &
+welec*gshieldc(j,i) &
- +welec*gshieldc_loc(j,) &
+ +welec*gshieldc_loc(j,i) &
+wcorr*gshieldc_ec(j,i) &
+wcorr*gshieldc_loc_ec(j,i) &
+wturn3*gshieldc_t3(j,i) &
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
- sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
- sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ sss_ele_cut=sscale_ele(1.0d0/(rij))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij))
sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
if (sss_ele_cut.le.0.0) cycle
if (sss.lt.1.0d0) then
sigder=fac*sigder
fac=rij*fac
fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
- /sigma(itypi,itypj)*rij-sss_grad/(1.0-sss)*rij &
+ *rij-sss_grad/(1.0-sss)*rij &
/sigmaii(itypi,itypj))
! fac=0.0d0
! Calculate the radial part of the gradient
rij=dsqrt(rrij)
sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
- sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
- sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ sss_ele_cut=sscale_ele(1.0d0/(rij))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij))
if (sss_ele_cut.le.0.0) cycle
if (sss.gt.0.0d0) then
sigder=fac*sigder
fac=rij*fac
fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
- /sigma(itypi,itypj)*rij+sss_grad/sss*rij &
+ *rij+sss_grad/sss*rij &
/sigmaii(itypi,itypj))
! fac=0.0d0
projects / unres4.git / blobdiff