X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2Fenergy_p_new_barrier.F;fp=source%2Funres%2Fsrc-HCD-5D%2Fenergy_p_new_barrier.F;h=9ec8107a38505730ffcacf204a23584d50c4c37e;hb=34afc6bc75a1158d5a5a7f08c5a7d84c8c1fb7c4;hp=2d94dc0f8e0fece55b82e226e7e53ad8336843d9;hpb=4f0aefffb3277492cac65561d39aa45e7dc34aa5;p=unres.git diff --git a/source/unres/src-HCD-5D/energy_p_new_barrier.F b/source/unres/src-HCD-5D/energy_p_new_barrier.F index 2d94dc0..9ec8107 100644 --- a/source/unres/src-HCD-5D/energy_p_new_barrier.F +++ b/source/unres/src-HCD-5D/energy_p_new_barrier.F @@ -2013,6 +2013,11 @@ C & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip double precision dist,sscale,sscagrad,sscagradlip,sscalelip double precision boxshift + double precision x0,y0,r012,rij12,facx0, + & facx02,afacx0,bfacx0,abfacx0,Afac,BBfac,Afacsig,Bfacsig +c alpha_GB=0.5d0 +c alpha_GB=0.01d0 +c alpha_GB1=1.0d0+1.0d0/alpha_GB evdw=0.0D0 ccccc energy_dec=.false. C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -2173,67 +2178,150 @@ c & " sig",sig," sig0ij",sig0ij c for diagnostics; uncomment c rij_shift=1.2*sig0ij C I hate to put IF's in the loops, but here don't have another choice!!!! - if (rij_shift.le.0.0D0) then - evdw=1.0D20 +c if (rij_shift.le.0.0D0) then + x0=alpha_GB*(sig-sig0ij) + if (energy_dec) write (iout,*) i,j," x0",x0 + if (rij_shift.le.x0) then +c sig=2.0d0*sig0ij + sigder=-sig*sigsq +c sigder=0.0d0 + fac=rij**expon + rij12=fac*fac +c rij12=1.0d0 + x0=alpha_GB*(sig-sig0ij) + facx0=1.0d0/x0**expon + facx02=facx0*facx0 + r012=((1.0d0+alpha_GB)*(sig-sig0ij))**(2*expon) + afacx0=aa*facx02 + bfacx0=bb*facx0 + abfacx0=afacx0+0.5d0*bfacx0 + Afac=alpha_GB1*abfacx0 + Afacsig=0.5d0*alpha_GB1*bfacx0/(sig-sig0ij) + BBfac=Afac-(afacx0+bfacx0) +c BBfac=0.0d0 + Bfacsig=(-alpha_GB1*(abfacx0+afacx0)+ + & (afacx0+afacx0+bfacx0))/(sig-sig0ij) +c Bfacsig=0.0d0 + Afac=Afac*r012 + Afacsig=Afacsig*r012 +c Afac=1.0d0 +c Afacsig=0.0d0 +c w(x)=4*eps*((1.0+1.0/alpha_GB)*(y0**12-0.5*y0**6)*(r0/x)**12-(1+1/alpha)*(y0**12-0.5*y0**6)+y0**12-y0**6) +c eps1=1.0d0 +c eps2rt=1.0d0 +c eps3rt=1.0d0 + e1 = eps1*eps2rt*eps3rt*Afac*rij12 + e2 = -eps1*eps2rt*eps3rt*BBfac + evdwij = e1+e2 + eps2der=evdwij*eps3rt + eps3der=evdwij*eps2rt +c eps2der=0.0d0 +c eps3der=0.0d0 +c eps1_om12=0.0d0 + evdwij=evdwij*eps2rt*eps3rt +c Afacsig=0.d0 +c Bfacsig=0.0d0 + if (lprn) then + write (iout,*) "aa",aa," bb",bb," sig0ij",sig0ij + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa + write (iout,'(2(a3,i3,2x),18(0pf9.5))') + & restyp(itypi),i,restyp(itypj),j, + & epsi,sigm,chi1,chi2,chip1,chip2, + & eps1,eps2rt**2,eps3rt**2,sig,sig0ij, + & eps1*eps2rt**2*eps3rt**2,om1,om2,om12, + & 1.0D0/rij,rij_shift, + & evdwij + endif + if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)') + & 'RE r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij + evdw=evdw+evdwij*sss +C Calculate gradient components. + e1=e1*eps2rt*eps3rt + sigder=-expon*eps1*eps2rt*eps2rt*eps3rt*eps3rt + & *(Afacsig*rij12-Bfacsig)*sigder + fac=-2.0d0*expon*e1*rij*rij +c print '(2i4,6f8.4)',i,j,sss,sssgrad* +c & evdwij,fac,sigma(itypi,itypj),expon + fac=fac+evdwij*sssgrad/sss*rij +c fac=0.0d0 +c write (iout,*) "sigder",sigder," fac",fac," e1",e1, +c & " e2",e2," sss",sss," sssgrad",sssgrad,"esp123", +c & eps1*eps2rt**2*eps3rt**2 +C Calculate the radial part of the gradient + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi +C gg_lipi(3)=0.0d0 +C gg_lipj(3)=0.0d0 + gg(1)=xj*fac + gg(2)=yj*fac + gg(3)=zj*fac +c evdw=1.0D20 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))') cd & restyp(itypi),i,restyp(itypj),j, cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) - return - endif - sigder=-sig*sigsq +c return + else + rij_shift=1.0D0/rij_shift + sigder=-sig*sigsq c--------------------------------------------------------------- - rij_shift=1.0D0/rij_shift - fac=rij_shift**expon + fac=rij_shift**expon C here to start with C if (c(i,3).gt. - faclip=fac - e1=fac*fac*aa - e2=fac*bb - evdwij=eps1*eps2rt*eps3rt*(e1+e2) - eps2der=evdwij*eps3rt - eps3der=evdwij*eps2rt + faclip=fac + e1=fac*fac*aa + e2=fac*bb + evdwij=eps1*eps2rt*eps3rt*(e1+e2) + eps2der=evdwij*eps3rt + eps3der=evdwij*eps2rt C write(63,'(2i3,2e10.3,2f10.5)') i,j,aa,bb, evdwij, C &((sslipi+sslipj)/2.0d0+ C &(2.0d0-sslipi-sslipj)/2.0d0) c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt, c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 - evdwij=evdwij*eps2rt*eps3rt - evdw=evdw+evdwij*sss - if (lprn) then - sigm=dabs(aa/bb)**(1.0D0/6.0D0) - epsi=bb**2/aa - write (iout,'(2(a3,i3,2x),17(0pf7.3))') - & restyp(itypi),i,restyp(itypj),j, - & epsi,sigm,chi1,chi2,chip1,chip2, - & eps1,eps2rt**2,eps3rt**2,sig,sig0ij, - & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, - & evdwij - endif + evdwij=evdwij*eps2rt*eps3rt + evdw=evdw+evdwij*sss + if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)') + & 'GB r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij + if (lprn) then + write (iout,*) "aa",aa," bb",bb," sig0ij",sig0ij + sigm=dabs(aa/bb)**(1.0D0/6.0D0) + epsi=bb**2/aa + write (iout,'(2(a3,i3,2x),17(0pf7.3))') + & restyp(itypi),i,restyp(itypj),j, + & epsi,sigm,chi1,chi2,chip1,chip2, + & eps1,eps2rt**2,eps3rt**2,sig,sig0ij, + & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, + & evdwij + endif - if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)') - & 'r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij C Calculate gradient components. - e1=e1*eps1*eps2rt**2*eps3rt**2 - fac=-expon*(e1+evdwij)*rij_shift - sigder=fac*sigder - fac=rij*fac -c print '(2i4,6f8.4)',i,j,sss,sssgrad* -c & evdwij,fac,sigma(itypi,itypj),expon - fac=fac+evdwij*sssgrad/sss*rij -c fac=0.0d0 + e1=e1*eps1*eps2rt**2*eps3rt**2 + fac=-expon*(e1+evdwij)*rij_shift + sigder=fac*sigder + fac=rij*fac +c print '(2i4,6f8.4)',i,j,sss,sssgrad* +c & evdwij,fac,sigma(itypi,itypj),expon + fac=fac+evdwij*sssgrad/sss*rij +c fac=0.0d0 C Calculate the radial part of the gradient - gg_lipi(3)=eps1*(eps2rt*eps2rt) + gg_lipi(3)=eps1*(eps2rt*eps2rt) & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) - gg_lipj(3)=ssgradlipj*gg_lipi(3) - gg_lipi(3)=gg_lipi(3)*ssgradlipi -C gg_lipi(3)=0.0d0 -C gg_lipj(3)=0.0d0 - gg(1)=xj*fac - gg(2)=yj*fac - gg(3)=zj*fac + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi +C gg_lipi(3)=0.0d0 +C gg_lipj(3)=0.0d0 + gg(1)=xj*fac + gg(2)=yj*fac + gg(3)=zj*fac + endif C Calculate angular part of the gradient. c call sc_grad_scale(sss) call sc_grad