X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.F90;h=094c72491f660f876d62349d3ae9b044fed4e3e5;hb=9e219e6555133b785d50b29273d99bdbea99f9f6;hp=3870bb0aaeef6d24551f54ef0c7b9693d9426464;hpb=da0ffdc0802840f3e36a2f41166ad72ced8f7845;p=unres4.git diff --git a/source/unres/energy.F90 b/source/unres/energy.F90 index 3870bb0..094c724 100644 --- a/source/unres/energy.F90 +++ b/source/unres/energy.F90 @@ -1,4 +1,4 @@ - module energy + module energy !----------------------------------------------------------------------------- use io_units use names @@ -793,6 +793,8 @@ call AFMforce(Eafmforce) else if (selfguide.gt.0) then call AFMvel(Eafmforce) + else + Eafmforce=0.0d0 endif endif if (tubemode.eq.1) then @@ -831,7 +833,7 @@ 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) @@ -916,8 +918,8 @@ ! 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 @@ -1042,8 +1044,8 @@ 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) @@ -1255,8 +1257,8 @@ 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) @@ -1331,11 +1333,11 @@ 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,& @@ -2013,8 +2015,8 @@ ! 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 @@ -2076,7 +2078,7 @@ 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 @@ -2848,6 +2850,7 @@ 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) @@ -3037,7 +3040,7 @@ 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-1),Ub2(1,i-2),CUgb2(1,i-2)) @@ -5053,9 +5056,9 @@ 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)) @@ -5086,7 +5089,7 @@ 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 @@ -5186,7 +5189,7 @@ 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) & @@ -11320,7 +11323,7 @@ +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) @@ -11551,7 +11554,7 @@ +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) & @@ -14019,8 +14022,8 @@ 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 @@ -14076,7 +14079,7 @@ 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 @@ -14310,8 +14313,8 @@ 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 @@ -14367,7 +14370,7 @@ 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 @@ -25513,7 +25516,7 @@ 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) @@ -25558,10 +25561,15 @@ 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!-------------------------------------------------------------------