+#endif
+ do i=loc_start,loc_end
+ usc_diff_i=0.0d0 ! argument of Gaussian for single k
+ guscdiff(i)=0.0d0 ! Sum of Gaussians over constr_homology ref structures
+c do j=ifrag_back(1,i,iset)+1,ifrag_back(2,i,iset)-1 ! Econstr_back legacy
+c write(iout,*) "xxtab, yytab, zztab"
+c write(iout,'(i5,3f8.2)') i,xxtab(i),yytab(i),zztab(i)
+ do k=1,constr_homology
+c
+ dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str?
+c Original sign inverted for calc of gradients (s. Econstr_back)
+ dyy=-yytpl(k,i)+yytab(i) ! ibid y
+ dzz=-zztpl(k,i)+zztab(i) ! ibid z
+c write(iout,*) "dxx, dyy, dzz"
+c write(iout,'(2i5,3f8.2)') k,i,dxx,dyy,dzz
+c
+ usc_diff_i=-0.5d0*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d rmvd from Gaussian argument
+c usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d?
+c uscdiffk(k)=usc_diff(i)
+ guscdiff2(k)=dexp(usc_diff_i) ! without min_scdiff
+ guscdiff(i)=guscdiff(i)+dexp(usc_diff_i) !Sum of Gaussians (pk)
+c write (iout,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j),
+c & xxref(j),yyref(j),zzref(j)
+ enddo
+c
+c Gradient
+c
+c Generalized expression for multiple Gaussian acc to that for a single
+c Gaussian in Econstr_back as instructed by AL (FP - 03/11/2014)
+c
+c Original implementation
+c sum_guscdiff=guscdiff(i)
+c
+c sum_sguscdiff=0.0d0
+c do k=1,constr_homology
+c sguscdiff=-guscdiff2(k)*dscdiff(k)*sigma_d(k,i)*waga_d !waga_d?
+c sguscdiff=-guscdiff3(k)*dscdiff(k)*sigma_d(k,i)*waga_d ! w min_uscdiff
+c sum_sguscdiff=sum_sguscdiff+sguscdiff
+c enddo
+c
+c Implementation of new expressions for gradient (Jan. 2015)
+c
+c grad_uscdiff=sum_sguscdiff/(sum_guscdiff*dtab) !?
+ do k=1,constr_homology
+c
+c New calculation of dxx, dyy, and dzz corrected by AL (07/11), was missing and wrong
+c before. Now the drivatives should be correct
+c
+ dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str?
+c Original sign inverted for calc of gradients (s. Econstr_back)
+ dyy=-yytpl(k,i)+yytab(i) ! ibid y
+ dzz=-zztpl(k,i)+zztab(i) ! ibid z
+c
+c New implementation
+c
+ sum_guscdiff=guscdiff2(k)*!(dsqrt(dxx*dxx+dyy*dyy+dzz*dzz))* -> wrong!
+ & sigma_d(k,i) ! for the grad wrt r'
+c sum_sguscdiff=sum_sguscdiff+sum_guscdiff
+c
+c
+c New implementation
+ sum_guscdiff = waga_d*sum_guscdiff
+ do jik=1,3
+ duscdiff(jik,i-1)=duscdiff(jik,i-1)+
+ & sum_guscdiff*(dXX_C1tab(jik,i)*dxx+
+ & dYY_C1tab(jik,i)*dyy+dZZ_C1tab(jik,i)*dzz)/guscdiff(i)
+ duscdiff(jik,i)=duscdiff(jik,i)+
+ & sum_guscdiff*(dXX_Ctab(jik,i)*dxx+
+ & dYY_Ctab(jik,i)*dyy+dZZ_Ctab(jik,i)*dzz)/guscdiff(i)
+ duscdiffx(jik,i)=duscdiffx(jik,i)+
+ & sum_guscdiff*(dXX_XYZtab(jik,i)*dxx+
+ & dYY_XYZtab(jik,i)*dyy+dZZ_XYZtab(jik,i)*dzz)/guscdiff(i)
+c
+#ifdef DEBUG
+ write(iout,*) "jik",jik,"i",i
+ write(iout,*) "dxx, dyy, dzz"
+ write(iout,'(2i5,3f8.2)') k,i,dxx,dyy,dzz
+ write(iout,*) "guscdiff2(",k,")",guscdiff2(k)
+c write(iout,*) "sum_sguscdiff",sum_sguscdiff
+cc write(iout,*) "dXX_Ctab(",jik,i,")",dXX_Ctab(jik,i)
+c write(iout,*) "dYY_Ctab(",jik,i,")",dYY_Ctab(jik,i)
+c write(iout,*) "dZZ_Ctab(",jik,i,")",dZZ_Ctab(jik,i)
+c write(iout,*) "dXX_C1tab(",jik,i,")",dXX_C1tab(jik,i)
+c write(iout,*) "dYY_C1tab(",jik,i,")",dYY_C1tab(jik,i)
+c write(iout,*) "dZZ_C1tab(",jik,i,")",dZZ_C1tab(jik,i)
+c write(iout,*) "dXX_XYZtab(",jik,i,")",dXX_XYZtab(jik,i)
+c write(iout,*) "dYY_XYZtab(",jik,i,")",dYY_XYZtab(jik,i)
+c write(iout,*) "dZZ_XYZtab(",jik,i,")",dZZ_XYZtab(jik,i)
+c write(iout,*) "duscdiff(",jik,i-1,")",duscdiff(jik,i-1)
+c write(iout,*) "duscdiff(",jik,i,")",duscdiff(jik,i)
+c write(iout,*) "duscdiffx(",jik,i,")",duscdiffx(jik,i)
+c endif
+#endif
+ enddo
+ enddo
+c
+c uscdiff(i)=-dLOG(guscdiff(i)/(ii-1)) ! Weighting by (ii-1) required?
+c usc_diff(i)=-dLOG(guscdiff(i)/constr_homology) ! + min_uscdiff ?
+c
+c write (iout,*) i," uscdiff",uscdiff(i)
+c
+c Put together deviations from local geometry