+ vbeta=vbeta+vbet(i,j)*e_gcont
+
+
+ if (dis(i,j) .ge. dfa_cutoff-2*dfa_cutoff_delta) then
+c gradient correction from gcont
+ de_gcont=vbet(i,j)*fprim_gcont/dis(i,j)
+ shetfx(i)=shetfx(i) + de_gcont*rx(i,j)
+ shetfy(i)=shetfy(i) + de_gcont*ry(i,j)
+ shetfz(i)=shetfz(i) + de_gcont*rz(i,j)
+
+ shetfx(j)=shetfx(j) - de_gcont*rx(i,j)
+ shetfy(j)=shetfy(j) - de_gcont*ry(i,j)
+ shetfz(j)=shetfz(j) - de_gcont*rz(i,j)
+
+c energy correction from gcont
+ vbet(i,j)=vbet(i,j)*e_gcont
+ vbetap(i,j)=vbetap(i,j)*e_gcont
+ vbetap1(i,j)=vbetap1(i,j)*e_gcont
+ vbetap2(i,j)=vbetap2(i,j)*e_gcont
+ vbetam(i,j)=vbetam(i,j)*e_gcont
+ vbetam1(i,j)=vbetam1(i,j)*e_gcont
+ vbetam2(i,j)=vbetam2(i,j)*e_gcont
+ endif
+