+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