- module energy
+ module energy
!-----------------------------------------------------------------------------
use io_units
use names
! print *,"Processor",myrank," BROADCAST iorder"
! FG master sets up the WEIGHTS_ array which will be broadcast to the
! FG slaves as WEIGHTS array.
- ! weights_(1)=wsc
+ weights_(1)=wsc
weights_(2)=wscp
weights_(3)=welec
weights_(4)=wcorr
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
.or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0 &
.or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then
#endif
-! write(iout,*),"just befor eelec call"
+! print *,"just befor eelec call"
call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
-! write (iout,*) "ELEC calc"
+! print *, "ELEC calc"
else
ees=0.0d0
evdw1=0.0d0
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",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
! to calculate the el-loc multibody terms of various order.
!
!AL el mu=0.0d0
+
#ifdef PARMAT
do i=ivec_start+2,ivec_end+2
#else
do i=3,nres+1
#endif
if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ if (itype(i-2,1).eq.0) then
+ iti = nloctyp
+ else
iti = itype2loc(itype(i-2,1))
+ endif
else
iti=nloctyp
endif
#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
else
iti1=nloctyp
endif
+! print *,i,iti
b1(1,i-2)=b(3,iti)
b1(2,i-2)=b(5,iti)
b2(1,i-2)=b(2,iti)
! 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
!d write (iout,*) 'mu2',mu2(:,i-2)
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) &
then
- call matmat2(CC(1,1,i-2),Ugder(1,1,i-2),CUgder(1,1,i-2))
+ call matmat2(CC(1,1,i-1),Ugder(1,1,i-2),CUgder(1,1,i-2))
call matmat2(DD(1,1,i-2),Ugder(1,1,i-2),DUgder(1,1,i-2))
call matmat2(Dtilde(1,1,i-2),Ug2der(1,1,i-2),DtUg2der(1,1,i-2))
- call matvec2(Ctilde(1,1,i-2),obrot_der(1,i-2),Ctobrder(1,i-2))
+ call matvec2(Ctilde(1,1,i-1),obrot_der(1,i-2),Ctobrder(1,i-2))
call matvec2(Dtilde(1,1,i-2),obrot2_der(1,i-2),Dtobr2der(1,i-2))
! Vectors and matrices dependent on a single virtual-bond dihedral.
- call matvec2(DD(1,1,i-2),b1tilde(1,iti1),auxvec(1))
+ call matvec2(DD(1,1,i-2),b1tilde(1,i-1),auxvec(1))
call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2))
call matvec2(Ug2der(1,1,i-2),auxvec(1),Ug2Db1tder(1,i-2))
- call matvec2(CC(1,1,i-2),Ub2(1,i-2),CUgb2(1,i-2))
- call matvec2(CC(1,1,i-2),Ub2der(1,i-2),CUgb2der(1,i-2))
+ call matvec2(CC(1,1,i-1),Ub2(1,i-2),CUgb2(1,i-2))
+ call matvec2(CC(1,1,i-1),Ub2der(1,i-2),CUgb2der(1,i-2))
call matmat2(EUg(1,1,i-2),CC(1,1,iti1),EUgC(1,1,i-2))
call matmat2(EUgder(1,1,i-2),CC(1,1,iti1),EUgCder(1,1,i-2))
call matmat2(EUg(1,1,i-2),DD(1,1,iti1),EUgD(1,1,i-2))
+a23*gmuij1(2)&
+a32*gmuij1(3)&
+a33*gmuij1(4))&
- *fac_shield(i)*fac_shield(j)
+ *fac_shield(i)*fac_shield(j)&
+ *sss_ele_cut
+
!c write(iout,*) "derivative over thatai"
!c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
!c & a33*gmuij1(4)
+a33*gmuij2(4)
gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+&
geel_loc_ij*wel_loc&
- *fac_shield(i)*fac_shield(j)
+ *fac_shield(i)*fac_shield(j)&
+ *sss_ele_cut
+
!c Derivative over j residue
geel_loc_ji=a22*gmuji1(1)&
gloc(nphi+j,icg)=gloc(nphi+j,icg)+&
geel_loc_ji*wel_loc&
- *fac_shield(i)*fac_shield(j)
+ *fac_shield(i)*fac_shield(j)&
+ *sss_ele_cut
+
geel_loc_ji=&
+a22*gmuji2(1)&
!c & a33*gmuji2(4)
gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+&
geel_loc_ji*wel_loc&
- *fac_shield(i)*fac_shield(j)
+ *fac_shield(i)*fac_shield(j)&
+ *sss_ele_cut
#endif
! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
a_temp(1,2)=a23
a_temp(2,1)=a32
a_temp(2,2)=a33
- iti1=itortyp(itype(i+1,1))
- iti2=itortyp(itype(i+2,1))
- iti3=itortyp(itype(i+3,1))
+ iti1=i+1
+ iti2=i+2
+ iti3=i+3
! write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3
call transpose2(EUg(1,1,i+1),e1t(1,1))
call transpose2(Eug(1,1,i+2),e2t(1,1))
call matvec2(ae3(1,1),gUb2(1,i+2),auxgvec(1))
!c auxilary matrix auxgEvec1 of E matix with Ub2 constant
call matvec2(gtae3(1,1),Ub2(1,i+2),auxgEvec3(1))
- s2=scalar2(b1(1,iti1),auxvec(1))
+ s2=scalar2(b1(1,i+1),auxvec(1))
!c derivative of theta i+1 with constant i+3
gs13=scalar2(gtb1(1,i+1),auxvec(1))
!c derivative of theta i+2 with constant i+1
call transpose2(EUgder(1,1,i+1),e1tder(1,1))
call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1))
call matvec2(auxmat(1,1),Ub2(1,i+3),auxvec(1))
- s1=scalar2(b1(1,iti2),auxvec(1))
+ s1=scalar2(b1(1,i+1),auxvec(1))
call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3) &
#ifdef TIMING
time01=MPI_Wtime()
#endif
+!#define DEBUG
#ifdef DEBUG
write (iout,*) "sum_gradient gvdwc, gvdwx"
do i=1,nres
! write (iout,'(i5,3f10.5,2x,f10.5)')
! & i,(gcorr4_turn(j,i),j=1,3),gel_loc_turn4(i)
! enddo
- write (iout,*) "gvdwc gvdwc_scp gvdwc_scpp"
- do i=1,nres
- write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
- i,(gvdwc(j,i),j=1,3),(gvdwc_scp(j,i),j=1,3),&
- (gvdwc_scpp(j,i),j=1,3)
- enddo
- write (iout,*) "gelc_long gvdwpp gel_loc_long"
- do i=1,nres
- write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
- i,(gelc_long(j,i),j=1,3),(gvdwpp(j,i),j=1,3),&
- (gelc_loc_long(j,i),j=1,3)
- enddo
+! write (iout,*) "gvdwc gvdwc_scp gvdwc_scpp"
+! do i=1,nres
+! write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
+! i,(gvdwc(j,i),j=1,3),(gvdwc_scp(j,i),j=1,3),&
+! (gvdwc_scpp(j,i),j=1,3)
+! enddo
+! write (iout,*) "gelc_long gvdwpp gel_loc_long"
+! do i=1,nres
+! write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
+! i,(gelc_long(j,i),j=1,3),(gvdwpp(j,i),j=1,3),&
+! (gelc_loc_long(j,i),j=1,3)
+! enddo
call flush(iout)
#endif
#ifdef SPLITELE
+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) &
! call intcartderiv
! call checkintcartgrad
call zerograd
- aincr=2.0D-5
+ aincr=1.0D-7
write(iout,*) 'Calling CHECK_ECARTINT.',aincr
nf=0
icall=0
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
rs,rthrp,rfourp,rsixp,reight,Irsixp,Ireight,Irtw,Irfourt, &
opt1,opt2,opt3,opt4,opt5,opt6,opt7,opt8,opt9,opt10,opt11,opt12,&
opt13,opt14,opt15,opt16,opt17,opt18,opt19, &
- Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip
+ Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip,&
+ ndiv,ndivi
real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,&
gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, &
dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, &
ecation_prot=0.0d0
! first lets calculate interaction with peptide groups
if (nres_molec(5).eq.0) return
- wconst=78
- wdip =1.092777950857032D2
- wdip=wdip/wconst
- wmodquad=-2.174122713004870D4
- wmodquad=wmodquad/wconst
- wquad1 = 3.901232068562804D1
- wquad1=wquad1/wconst
- wquad2 = 3
- wquad2=wquad2/wconst
- wvan1 = 0.1
- wvan2 = 6
itmp=0
do i=1,4
itmp=itmp+nres_molec(i)
if (zi.lt.0) zi=zi+boxzsize
do j=itmp+1,itmp+nres_molec(5)
+! print *,"WTF",itmp,j,i
+! all parameters were for Ca2+ to approximate single charge divide by two
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+ wconst=78*ndiv
+ wdip =1.092777950857032D2
+ wdip=wdip/wconst
+ wmodquad=-2.174122713004870D4
+ wmodquad=wmodquad/wconst
+ wquad1 = 3.901232068562804D1
+ wquad1=wquad1/wconst
+ wquad2 = 3
+ wquad2=wquad2/wconst
+ wvan1 = 0.1
+ wvan2 = 6
+! itmp=0
+
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
E2 = -wquad1*Irthrp*wquad2+wvan1*(wvan2**12*Irtwelv-&
2*wvan2**6*Irsistp)
ecation_prot = ecation_prot+E1+E2
+! print *,"ecatprot",i,j,ecation_prot,rcpm
dE1dr = -2*costhet*wdip*Irthrp-&
(4*wmodquad*Irfiftp+3*wquad1*Irfourp)*sin2thet
dE2dr = 3*wquad1*wquad2*Irfourp- &
enddo
cm1mag=sqrt(cm1(1)**2+cm1(2)**2+cm1(3)**2)
do j=itmp+1,itmp+nres_molec(5)
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
endif
! enddo
! enddo
- if(itype(i,1).eq.15.or.itype(i,1).eq.16) then
+! 15- Glu 16-Asp
+ if((itype(i,1).eq.15.or.itype(i,1).eq.16).or.&
+ ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.&
+ (itype(i,1).eq.25))) then
if(itype(i,1).eq.16) then
inum=1
else
! The weights of the energy function calculated from
!The quantum mechanical GAMESS simulations of calcium with ASP/GLU
- wh2o=78
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ ndivi=0.5
+ else
+ ndivi=1.0
+ endif
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
+ wh2o=78*ndivi*ndiv
wc = vcatprm(1)
wc=wc/wh2o
wdip =vcatprm(2)
v1dpv2 = v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
! The weights of the energy function calculated from
!The quantum mechanical GAMESS simulations of ASN/GLN with calcium
- wh2o=78
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
+ wh2o=78*ndiv
wdip =vcatprm(2)
wdip=wdip/wh2o
wquad1 =vcatprm(3)
use calc_data
use comm_momo
real (kind=8) :: facd3, facd4, federmaus, adler,&
- Ecl,Egb,Epol,Fisocav,Elj,Fgb
+ Ecl,Egb,Epol,Fisocav,Elj,Fgb,debkap
! integer :: k
!c! Epol and Gpol analytical parameters
alphapol1 = alphapol(itypi,itypj)
dGCLdOM12 = 0.0d0
ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq))
Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0)
- Egb = -(332.0d0 * Qij * eps_inout_fac) / Fgb
+ debkap=debaykap(itypi,itypj)
+ Egb = -(332.0d0 * Qij *&
+ (1.0/eps_in-dexp(-debkap*Fgb)/eps_out)) / Fgb
! print *,"EGB WTF",Qij,eps_inout_fac,Fgb,itypi,itypj,eps_in,eps_out
!c! Derivative of Egb is Ggb...
- dGGBdFGB = -(-332.0d0 * Qij * eps_inout_fac) / (Fgb * Fgb)
+ dGGBdFGB = -(-332.0d0 * Qij * &
+ (1.0/eps_in-dexp(-debkap*Fgb)/eps_out))/(Fgb*Fgb)&
+ -(332.0d0 * Qij *&
+ (dexp(-debkap*Fgb)*debkap/eps_out))/ Fgb
dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )/ ( 2.0d0 * Fgb )
dGGBdR = dGGBdFGB * dFGBdR
!c!-------------------------------------------------------------------
projects / unres4.git / blobdiff