+ enddo
+! write (iout,*) "godl",(godl(k),k=1,constr_homology) ! exponents
+! write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps
+#ifdef DEBUG
+ write (iout,*) "godl",(godl(k),k=1,constr_homology) ! exponents
+ write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps
+#endif
+ if (waga_dist.ge.0.0d0) then
+!
+! For Gaussian-type Urestr
+!
+ odleg=odleg-dLOG(odleg2/constr_homology)+min_odl
+!
+! For Lorentzian-type Urestr
+!
+ else
+ odleg=odleg+odleg2/constr_homology
+ endif
+!
+! write (iout,*) "odleg",odleg ! sum of -ln-s
+! Gradient
+!
+! For Gaussian-type Urestr
+!
+ if (waga_dist.ge.0.0d0) sum_godl=odleg2
+ sum_sgodl=0.0d0
+ do k=1,constr_homology
+! godl=dexp(((-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2))
+! & *waga_dist)+min_odl
+! sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist
+!
+ if(.not.l_homo(k,ii)) cycle
+ if (waga_dist.ge.0.0d0) then
+! For Gaussian-type Urestr
+!
+ sgodl=-godl(k)*distance(k)*sigma_odl(k,ii) ! waga_dist rmvd
+!
+! For Lorentzian-type Urestr
+!
+ else
+ sgodl=-2*sigma_odlir(k,ii)*(distance(k)/(distance(k)**2+ &
+ sigma_odlir(k,ii)**2)**2)
+ endif
+ sum_sgodl=sum_sgodl+sgodl
+
+! sgodl2=sgodl2+sgodl
+! write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE1"
+! write(iout,*) "constr_homology=",constr_homology
+! write(iout,*) i, j, k, "TEST K"
+ enddo
+! print *, "ok",iset
+ if (waga_dist.ge.0.0d0) then
+!
+! For Gaussian-type Urestr
+!
+ grad_odl3=waga_homology(iset)*waga_dist &
+ *sum_sgodl/(sum_godl*dij)
+! print *, "ok"
+!
+! For Lorentzian-type Urestr
+!
+ else
+! Original grad expr modified by analogy w Gaussian-type Urestr grad
+! grad_odl3=-waga_homology(iset)*waga_dist*sum_sgodl
+ grad_odl3=-waga_homology(iset)*waga_dist* &
+ sum_sgodl/(constr_homology*dij)
+! print *, "ok2"
+ endif
+!
+! grad_odl3=sum_sgodl/(sum_godl*dij)
+
+
+! write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE2"
+! write(iout,*) (distance(i,j,k)**2), (2*(sigma_odl(i,j,k))**2),
+! & (-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2))
+
+!cc write(iout,*) godl, sgodl, grad_odl3
+
+! grad_odl=grad_odl+grad_odl3
+
+ do jik=1,3
+ ggodl=grad_odl3*(c(jik,i)-c(jik,j))
+!cc write(iout,*) c(jik,i+1), c(jik,j+1), (c(jik,i+1)-c(jik,j+1))
+!cc write(iout,746) "GRAD_ODL_1", i, j, jik, ggodl,
+!cc & ghpbc(jik,i+1), ghpbc(jik,j+1)
+ ghpbc(jik,i)=ghpbc(jik,i)+ggodl
+ ghpbc(jik,j)=ghpbc(jik,j)-ggodl
+!cc write(iout,746) "GRAD_ODL_2", i, j, jik, ggodl,
+!cc & ghpbc(jik,i+1), ghpbc(jik,j+1)
+! if (i.eq.25.and.j.eq.27) then
+! write(iout,*) "jik",jik,"i",i,"j",j
+! write(iout,*) "sum_sgodl",sum_sgodl,"sgodl",sgodl
+! write(iout,*) "grad_odl3",grad_odl3
+! write(iout,*) "c(",jik,i,")",c(jik,i),"c(",jik,j,")",c(jik,j)
+! write(iout,*) "ggodl",ggodl
+! write(iout,*) "ghpbc(",jik,i,")",
+! & ghpbc(jik,i),"ghpbc(",jik,j,")",
+! & ghpbc(jik,j)
+! endif
+ enddo
+!cc write(iout,778)"TEST: odleg2=", odleg2, "DLOG(odleg2)=",
+!cc & dLOG(odleg2),"-odleg=", -odleg
+
+ enddo ! ii-loop for dist
+#ifdef DEBUG
+ write(iout,*) "------- dist restrs end -------"
+! if (waga_angle.eq.1.0d0 .or. waga_theta.eq.1.0d0 .or.
+! & waga_d.eq.1.0d0) call sum_gradient
+#endif
+! Pseudo-energy and gradient from dihedral-angle restraints from
+! homology templates
+! write (iout,*) "End of distance loop"
+! call flush(iout)
+ kat=0.0d0
+! write (iout,*) idihconstr_start_homo,idihconstr_end_homo
+#ifdef DEBUG
+ write(iout,*) "------- dih restrs start -------"
+ do i=idihconstr_start_homo,idihconstr_end_homo
+ write (iout,*) "gloc_init(",i,icg,")",gloc(i,icg)
+ enddo
+#endif
+ do i=idihconstr_start_homo,idihconstr_end_homo
+ kat2=0.0d0
+! betai=beta(i,i+1,i+2,i+3)
+ betai = phi(i)
+! write (iout,*) "betai =",betai
+ do k=1,constr_homology
+ dih_diff(k)=pinorm(dih(k,i)-betai)
+!d write (iout,'(a8,2i4,2f15.8)') "dih_diff",i,k,dih_diff(k)
+!d & ,sigma_dih(k,i)
+! if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)=
+! & -(6.28318-dih_diff(i,k))
+! if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)=
+! & 6.28318+dih_diff(i,k)
+#ifdef OLD_DIHED
+ kat3=-0.5d0*dih_diff(k)**2*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument
+#else
+ kat3=(dcos(dih_diff(k))-1)*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument
+#endif
+! kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i)
+ gdih(k)=dexp(kat3)
+ kat2=kat2+gdih(k)
+! write(iout,*) "kat2=", kat2, "exp(kat3)=", exp(kat3)
+! write(*,*)""
+ enddo
+! write (iout,*) "gdih",(gdih(k),k=1,constr_homology) ! exps
+! write (iout,*) "i",i," betai",betai," kat2",kat2 ! sum of exps
+#ifdef DEBUG
+ write (iout,*) "i",i," betai",betai," kat2",kat2
+ write (iout,*) "gdih",(gdih(k),k=1,constr_homology)
+#endif
+ if (kat2.le.1.0d-14) cycle
+ kat=kat-dLOG(kat2/constr_homology)
+! write (iout,*) "kat",kat ! sum of -ln-s
+
+!cc write(iout,778)"TEST: kat2=", kat2, "DLOG(kat2)=",
+!cc & dLOG(kat2), "-kat=", -kat
+
+! ----------------------------------------------------------------------
+! Gradient
+! ----------------------------------------------------------------------
+
+ sum_gdih=kat2
+ sum_sgdih=0.0d0
+ do k=1,constr_homology
+#ifdef OLD_DIHED
+ sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i) ! waga_angle rmvd
+#else
+ sgdih=-gdih(k)*dsin(dih_diff(k))*sigma_dih(k,i) ! waga_angle rmvd
+#endif
+! sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle
+ sum_sgdih=sum_sgdih+sgdih
+ enddo
+! grad_dih3=sum_sgdih/sum_gdih
+ grad_dih3=waga_homology(iset)*waga_angle*sum_sgdih/sum_gdih
+! print *, "ok3"
+
+! write(iout,*)i,k,gdih,sgdih,beta(i+1,i+2,i+3,i+4),grad_dih3
+!cc write(iout,747) "GRAD_KAT_1", i, nphi, icg, grad_dih3,
+!cc & gloc(nphi+i-3,icg)
+ gloc(i-3,icg)=gloc(i-3,icg)+grad_dih3
+! if (i.eq.25) then
+! write(iout,*) "i",i,"icg",icg,"gloc(",i,icg,")",gloc(i,icg)
+! endif
+!cc write(iout,747) "GRAD_KAT_2", i, nphi, icg, grad_dih3,
+!cc & gloc(nphi+i-3,icg)
+
+ enddo ! i-loop for dih
+#ifdef DEBUG
+ write(iout,*) "------- dih restrs end -------"
+#endif
+
+! Pseudo-energy and gradient for theta angle restraints from
+! homology templates
+! FP 01/15 - inserted from econstr_local_test.F, loop structure
+! adapted
+
+!
+! For constr_homology reference structures (FP)
+!
+! Uconst_back_tot=0.0d0
+ Eval=0.0d0
+ Erot=0.0d0
+! Econstr_back legacy
+ do i=1,nres
+! do i=ithet_start,ithet_end
+ dutheta(i)=0.0d0
+ enddo
+! do i=loc_start,loc_end
+ do i=-1,nres
+ do j=1,3
+ duscdiff(j,i)=0.0d0
+ duscdiffx(j,i)=0.0d0
+ enddo
+ enddo
+!
+! do iref=1,nref
+! write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
+! write (iout,*) "waga_theta",waga_theta
+ if (waga_theta.gt.0.0d0) then
+#ifdef DEBUG
+ write (iout,*) "usampl",usampl
+ write(iout,*) "------- theta restrs start -------"
+! do i=ithet_start,ithet_end
+! write (iout,*) "gloc_init(",nphi+i,icg,")",gloc(nphi+i,icg)
+! enddo
+#endif
+! write (iout,*) "maxres",maxres,"nres",nres
+
+ do i=ithet_start,ithet_end
+!
+! do i=1,nfrag_back
+! ii = ifrag_back(2,i,iset)-ifrag_back(1,i,iset)
+!
+! Deviation of theta angles wrt constr_homology ref structures
+!
+ utheta_i=0.0d0 ! argument of Gaussian for single k
+ gutheta_i=0.0d0 ! Sum of Gaussians over constr_homology ref structures
+! do j=ifrag_back(1,i,iset)+2,ifrag_back(2,i,iset) ! original loop
+! over residues in a fragment
+! write (iout,*) "theta(",i,")=",theta(i)
+ do k=1,constr_homology
+!
+! dtheta_i=theta(j)-thetaref(j,iref)
+! dtheta_i=thetaref(k,i)-theta(i) ! original form without indexing
+ theta_diff(k)=thetatpl(k,i)-theta(i)
+!d write (iout,'(a8,2i4,2f15.8)') "theta_diff",i,k,theta_diff(k)
+!d & ,sigma_theta(k,i)
+
+!
+ utheta_i=-0.5d0*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta rmvd from Gaussian argument
+! utheta_i=-0.5d0*waga_theta*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta?
+ gtheta(k)=dexp(utheta_i) ! + min_utheta_i?
+ gutheta_i=gutheta_i+gtheta(k) ! Sum of Gaussians (pk)
+! Gradient for single Gaussian restraint in subr Econstr_back
+! dutheta(j-2)=dutheta(j-2)+wfrag_back(1,i,iset)*dtheta_i/(ii-1)
+!
+ enddo
+! write (iout,*) "gtheta",(gtheta(k),k=1,constr_homology) ! exps
+! write (iout,*) "i",i," gutheta_i",gutheta_i ! sum of exps
+
+!
+! Gradient for multiple Gaussian restraint
+ sum_gtheta=gutheta_i
+ sum_sgtheta=0.0d0
+ do k=1,constr_homology
+! New generalized expr for multiple Gaussian from Econstr_back
+ sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i) ! waga_theta rmvd
+!
+! sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i)*waga_theta ! right functional form?
+ sum_sgtheta=sum_sgtheta+sgtheta ! cum variable
+ enddo
+! Final value of gradient using same var as in Econstr_back
+ gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg) &
+ +sum_sgtheta/sum_gtheta*waga_theta &
+ *waga_homology(iset)
+! print *, "ok4"
+
+! dutheta(i-2)=sum_sgtheta/sum_gtheta*waga_theta
+! & *waga_homology(iset)
+! dutheta(i)=sum_sgtheta/sum_gtheta
+!
+! Uconst_back=Uconst_back+waga_theta*utheta(i) ! waga_theta added as weight
+ Eval=Eval-dLOG(gutheta_i/constr_homology)
+! write (iout,*) "utheta(",i,")=",utheta(i) ! -ln of sum of exps
+! write (iout,*) "Uconst_back",Uconst_back ! sum of -ln-s
+! Uconst_back=Uconst_back+utheta(i)
+ enddo ! (i-loop for theta)
+#ifdef DEBUG
+ write(iout,*) "------- theta restrs end -------"
+#endif
+ endif
+!
+! Deviation of local SC geometry
+!
+! Separation of two i-loops (instructed by AL - 11/3/2014)
+!
+! write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
+! write (iout,*) "waga_d",waga_d
+
+#ifdef DEBUG
+ write(iout,*) "------- SC restrs start -------"
+ write (iout,*) "Initial duscdiff,duscdiffx"
+ do i=loc_start,loc_end
+ write (iout,*) i,(duscdiff(jik,i),jik=1,3), &
+ (duscdiffx(jik,i),jik=1,3)
+ enddo
+#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
+! do j=ifrag_back(1,i,iset)+1,ifrag_back(2,i,iset)-1 ! Econstr_back legacy
+! write(iout,*) "xxtab, yytab, zztab"
+! write(iout,'(i5,3f8.2)') i,xxtab(i),yytab(i),zztab(i)
+ do k=1,constr_homology
+!
+ dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str?
+! 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
+! write(iout,*) "dxx, dyy, dzz"
+!d write(iout,'(2i5,4f8.2)') k,i,dxx,dyy,dzz,sigma_d(k,i)
+!
+ usc_diff_i=-0.5d0*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d rmvd from Gaussian argument
+! usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d?
+! uscdiffk(k)=usc_diff(i)
+ guscdiff2(k)=dexp(usc_diff_i) ! without min_scdiff
+! write(iout,*) "i",i," k",k," sigma_d",sigma_d(k,i),
+! & " guscdiff2",guscdiff2(k)
+ guscdiff(i)=guscdiff(i)+guscdiff2(k) !Sum of Gaussians (pk)
+! write (iout,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j),
+! & xxref(j),yyref(j),zzref(j)
+ enddo
+!
+! Gradient
+!
+! Generalized expression for multiple Gaussian acc to that for a single
+! Gaussian in Econstr_back as instructed by AL (FP - 03/11/2014)
+!
+! Original implementation
+! sum_guscdiff=guscdiff(i)
+!
+! sum_sguscdiff=0.0d0
+! do k=1,constr_homology
+! sguscdiff=-guscdiff2(k)*dscdiff(k)*sigma_d(k,i)*waga_d !waga_d?
+! sguscdiff=-guscdiff3(k)*dscdiff(k)*sigma_d(k,i)*waga_d ! w min_uscdiff
+! sum_sguscdiff=sum_sguscdiff+sguscdiff
+! enddo
+!
+! Implementation of new expressions for gradient (Jan. 2015)
+!
+! grad_uscdiff=sum_sguscdiff/(sum_guscdiff*dtab) !?
+ do k=1,constr_homology
+!
+! New calculation of dxx, dyy, and dzz corrected by AL (07/11), was missing and wrong
+! before. Now the drivatives should be correct
+!
+ dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str?
+! 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
+ sum_guscdiff=guscdiff2(k)* &!(dsqrt(dxx*dxx+dyy*dyy+dzz*dzz))* -> wrong!
+ sigma_d(k,i) ! for the grad wrt r'
+! sum_sguscdiff=sum_sguscdiff+sum_guscdiff
+
+!
+! New implementation
+ sum_guscdiff = waga_homology(iset)*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)
+! print *, "ok5"
+!
+#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)
+ write(iout,*) "sum_sguscdiff",sum_guscdiff,waga_homology(iset),waga_d
+ write(iout,*) "dXX_Ctab(",jik,i,")",dXX_Ctab(jik,i)
+ write(iout,*) "dYY_Ctab(",jik,i,")",dYY_Ctab(jik,i)
+ write(iout,*) "dZZ_Ctab(",jik,i,")",dZZ_Ctab(jik,i)
+ write(iout,*) "dXX_C1tab(",jik,i,")",dXX_C1tab(jik,i)
+ write(iout,*) "dYY_C1tab(",jik,i,")",dYY_C1tab(jik,i)
+ write(iout,*) "dZZ_C1tab(",jik,i,")",dZZ_C1tab(jik,i)
+ write(iout,*) "dXX_XYZtab(",jik,i,")",dXX_XYZtab(jik,i)
+ write(iout,*) "dYY_XYZtab(",jik,i,")",dYY_XYZtab(jik,i)
+ write(iout,*) "dZZ_XYZtab(",jik,i,")",dZZ_XYZtab(jik,i)
+ write(iout,*) "duscdiff(",jik,i-1,")",duscdiff(jik,i-1)
+ write(iout,*) "duscdiff(",jik,i,")",duscdiff(jik,i)
+ write(iout,*) "duscdiffx(",jik,i,")",duscdiffx(jik,i)
+! endif
+#endif
+ enddo
+ enddo
+! print *, "ok6"
+!
+! uscdiff(i)=-dLOG(guscdiff(i)/(ii-1)) ! Weighting by (ii-1) required?
+! usc_diff(i)=-dLOG(guscdiff(i)/constr_homology) ! + min_uscdiff ?
+!
+! write (iout,*) i," uscdiff",uscdiff(i)
+!
+! Put together deviations from local geometry
+
+! Uconst_back=Uconst_back+wfrag_back(1,i,iset)*utheta(i)+
+! & wfrag_back(3,i,iset)*uscdiff(i)
+ Erot=Erot-dLOG(guscdiff(i)/constr_homology)
+! write (iout,*) "usc_diff(",i,")=",usc_diff(i) ! -ln of sum of exps
+! write (iout,*) "Uconst_back",Uconst_back ! cum sum of -ln-s
+! Uconst_back=Uconst_back+usc_diff(i)
+!
+! Gradient of multiple Gaussian restraint (FP - 04/11/2014 - right?)
+!
+! New implment: multiplied by sum_sguscdiff
+!
+
+ enddo ! (i-loop for dscdiff)
+
+! endif
+
+#ifdef DEBUG
+ write(iout,*) "------- SC restrs end -------"
+ write (iout,*) "------ After SC loop in e_modeller ------"
+ do i=loc_start,loc_end
+ write (iout,*) "i",i," gradc",(gradc(j,i,icg),j=1,3)
+ write (iout,*) "i",i," gradx",(gradx(j,i,icg),j=1,3)
+ enddo
+ if (waga_theta.eq.1.0d0) then
+ write (iout,*) "in e_modeller after SC restr end: dutheta"
+ do i=ithet_start,ithet_end
+ write (iout,*) i,dutheta(i)
+ enddo
+ endif
+ if (waga_d.eq.1.0d0) then
+ write (iout,*) "e_modeller after SC loop: duscdiff/x"
+ do i=1,nres
+ write (iout,*) i,(duscdiff(j,i),j=1,3)
+ write (iout,*) i,(duscdiffx(j,i),j=1,3)
+ enddo
+ endif
+#endif
+
+! Total energy from homology restraints
+#ifdef DEBUG
+ write (iout,*) "odleg",odleg," kat",kat
+#endif
+!
+! Addition of energy of theta angle and SC local geom over constr_homologs ref strs
+!
+! ehomology_constr=odleg+kat
+!
+! For Lorentzian-type Urestr
+!
+
+ if (waga_dist.ge.0.0d0) then
+!
+! For Gaussian-type Urestr
+!
+ ehomology_constr=(waga_dist*odleg+waga_angle*kat+ &
+ waga_theta*Eval+waga_d*Erot)*waga_homology(iset)
+! write (iout,*) "ehomology_constr=",ehomology_constr
+! print *, "ok7"
+ else
+!
+! For Lorentzian-type Urestr
+!
+ ehomology_constr=(-waga_dist*odleg+waga_angle*kat+ &
+ waga_theta*Eval+waga_d*Erot)*waga_homology(iset)
+! write (iout,*) "ehomology_constr=",ehomology_constr
+ print *, "ok8"
+ endif
+#ifdef DEBUG
+ write (iout,*) "odleg",waga_dist,odleg," kat",waga_angle,kat, &
+ "Eval",waga_theta,eval, &
+ "Erot",waga_d,Erot
+ write (iout,*) "ehomology_constr",ehomology_constr
+#endif
+ return
+!
+! FP 01/15 end
+!
+ 748 format(a8,f12.3,a6,f12.3,a7,f12.3)
+ 747 format(a12,i4,i4,i4,f8.3,f8.3)
+ 746 format(a12,i4,i4,i4,f8.3,f8.3,f8.3)
+ 778 format(a7,1X,f10.3,1X,a4,1X,f10.3,1X,a5,1X,f10.3)
+ 779 format(i3,1X,i3,1X,i2,1X,a7,1X,f7.3,1X,a7,1X,f7.3,1X,a13,1X, &
+ f7.3,1X,a17,1X,f9.3,1X,a10,1X,f8.3,1X,a10,1X,f8.3)
+ end subroutine e_modeller
+
+!----------------------------------------------------------------------------