X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Funres%2Fenergy.F90;h=6c60e98dfc180682913be231305b7686444a1ffe;hb=aa10e88041f16954d81953a1ecfbfc956e089f72;hp=52fb0ff35b4211deb75f46c614fb14440f81b3cf;hpb=93156c3542905e5b694a64496717b10b9c6a49a8;p=unres4.git diff --git a/source/unres/energy.F90 b/source/unres/energy.F90 index 52fb0ff..6c60e98 100644 --- a/source/unres/energy.F90 +++ b/source/unres/energy.F90 @@ -1,4 +1,4 @@ - module energy + module energy !----------------------------------------------------------------------------- use io_units use names @@ -378,13 +378,20 @@ 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 chainbuild_cart endif - print *,'Processor',myrank,' calling etotal ipot=',ipot - print *,'Processor',myrank,' nnt=',nnt,' nct=',nct +! print *,'Processor',myrank,' calling etotal ipot=',ipot +! print *,'Processor',myrank,' nnt=',nnt,' nct=',nct #else ! if (modecalc.eq.12.or.modecalc.eq.14) then ! call int_from_cart1(.false.) @@ -793,6 +800,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 +840,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) @@ -852,7 +861,7 @@ 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 @@ -916,8 +925,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 +1051,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) @@ -1187,6 +1196,11 @@ 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 @@ -1255,8 +1269,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 +1345,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 +2027,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 +2090,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 @@ -2837,10 +2851,16 @@ #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 @@ -3013,7 +3033,7 @@ ! 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 @@ -11321,7 +11341,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) @@ -11552,7 +11572,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) & @@ -14020,8 +14040,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 @@ -14077,7 +14097,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 @@ -14311,8 +14331,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 @@ -14368,7 +14388,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 @@ -25514,7 +25534,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) @@ -25559,10 +25579,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!-------------------------------------------------------------------