correction of dihedral homology restraints phi(i+3)

diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F
index 34807fe..4e4a386 100644 (file)
--- a/source/cluster/wham/src-M/energy_p_new.F
+++ b/source/cluster/wham/src-M/energy_p_new.F
@@ -4177,7 +4177,7 @@ c      write (iout,*) idihconstr_start_homo,idihconstr_end_homo
       do i=idihconstr_start_homo,idihconstr_end_homo
         kat2=0.0d0
 c        betai=beta(i,i+1,i+2,i+3)
-        betai = phi(i+3)
+        betai = phi(i)
 c       write (iout,*) "betai =",betai
         do k=1,constr_homology
           dih_diff(k)=pinorm(dih(k,i)-betai)

diff --git a/source/cluster/wham/src-M/initialize_p.F b/source/cluster/wham/src-M/initialize_p.F
index 8b01a36..1c83869 100644 (file)
--- a/source/cluster/wham/src-M/initialize_p.F
+++ b/source/cluster/wham/src-M/initialize_p.F
@@ -620,10 +620,10 @@ c      include 'COMMON.SETUP'
      &   " lim_dih",lim_dih
 #ifdef MPL
       call int_bounds(lim_odl,link_start_homo,link_end_homo)
-      call int_bounds(lim_dih-nnt+1,idihconstr_start_homo,
+      call int_bounds(lim_dih,idihconstr_start_homo,
      &  idihconstr_end_homo)
-      idihconstr_start_homo=idihconstr_start_homo+nnt-1
-      idihconstr_end_homo=idihconstr_end_homo+nnt-1
+      idihconstr_start_homo=idihconstr_start_homo+nnt-1+3
+      idihconstr_end_homo=idihconstr_end_homo+nnt-1+3
       if (me.eq.king .or. .not. out1file)
      &  write (iout,*) 'Processor',fg_rank,' CG group',kolor,
      &  ' absolute rank',MyRank,
@@ -634,8 +634,8 @@ c      include 'COMMON.SETUP'
 #else
       link_start_homo=1
       link_end_homo=lim_odl
-      idihconstr_start_homo=nnt
-      idihconstr_end_homo=lim_dih
+      idihconstr_start_homo=nnt+3
+      idihconstr_end_homo=lim_dih+3
       write (iout,*)
      &  ' lim_odl',lim_odl,' link_start=',link_start_homo,
      &  ' link_end',link_end_homo,' lim_dih',lim_dih,

diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F
index 1680dab..51ac51f 100644 (file)
--- a/source/cluster/wham/src-M/readrtns.F
+++ b/source/cluster/wham/src-M/readrtns.F
@@ -968,6 +968,7 @@ c====-------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.GEO'
       include 'COMMON.INTERACT'
+      include 'COMMON.NAMES'
       include 'COMMON.HOMRESTR'
 c
 c For new homol impl
@@ -992,9 +993,11 @@ c    &    sigma_odl_temp(maxres,maxres,max_template)
 c
 c     FP - Nov. 2014 Temporary specifications for new vars
 c
-      double precision rescore_tmp,x12,y12,z12,rescore2_tmp
+      double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
+     &                 rescore3_tmp
       double precision, dimension (max_template,maxres) :: rescore
       double precision, dimension (max_template,maxres) :: rescore2
+      double precision, dimension (max_template,maxres) :: rescore3
       character*24 tpl_k_rescore
 c -----------------------------------------------------------------
 c Reading multiple PDB ref structures and calculation of retraints
@@ -1048,8 +1051,6 @@ cd      call flush(iout)
 c
 c  New
 c
-      lim_theta=0
-      lim_xx=0
 c
 c  Reading HM global scores (prob not required)
 c
@@ -1145,14 +1146,18 @@ c
 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
           open (ientin,file=tpl_k_rescore,status='old')
           if (nnt.gt.1) rescore(k,1)=0.0d0
-          do irec=nnt,maxdim ! loop for reading res sim 
+          do irec=nnt,nct ! loop for reading res sim 
             if (read2sigma) then
              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
-     &                                idomain_tmp
+     &                                rescore3_tmp,idomain_tmp
              i_tmp=i_tmp+nnt-1
              idomain(k,i_tmp)=idomain_tmp
              rescore(k,i_tmp)=rescore_tmp
              rescore2(k,i_tmp)=rescore2_tmp
+             rescore3(k,i_tmp)=rescore3_tmp
+             write(iout,'(a7,i5,3f10.5,i5)') "rescore",
+     &                      i_tmp,rescore2_tmp,rescore_tmp,
+     &                                rescore3_tmp,idomain_tmp
             else
              idomain(k,irec)=1
              read (ientin,*,end=1401) rescore_tmp
@@ -1252,7 +1257,8 @@ c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
 c   Instead of res sim other local measure of b/b str reliability possible
-            sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
+            if (sigma_dih(k,i).ne.0)
+     &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
           enddo
           lim_dih=nct-nnt-2 
@@ -1284,14 +1290,14 @@ c            read (ientin,*) sigma_theta(k,i) ! 1st variant
              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
      &                        rescore(k,i-2))/3.0
 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
-             sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
+             if (sigma_theta(k,i).ne.0)
+     &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 
 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 c                             rescore(k,i-2) !  right expression ?
 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
           enddo
         endif
-        lim_theta=nct-nnt-1 
 
         if (waga_d.gt.0.0d0) then
 c       open (ientin,file=tpl_k_sigma_d,status='old')
@@ -1317,15 +1323,14 @@ c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
-               sigma_d(k,i)=rescore(k,i) !  right expression ?
-               sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
+               sigma_d(k,i)=rescore3(k,i) !  right expression ?
+               if (sigma_d(k,i).ne.0)
+     &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
 
 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
 c              read (ientin,*) sigma_d(k,i) ! 1st variant
-               if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
           enddo
-          lim_xx=nct-nnt+1 
         endif
       enddo
 c
@@ -1374,17 +1379,19 @@ cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
      &  ki=1,constr_homology)
        enddo
        write (iout,*) "Dihedral angle restraints from templates"
-       do i=nnt+3,lim_dih
-        write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*dih(ki,i),
+       do i=nnt+3,nct
+        write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
+     &      (rad2deg*dih(ki,i),
      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
        enddo
        write (iout,*) "Virtual-bond angle restraints from templates"
-       do i=nnt+2,lim_theta
-        write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
+       do i=nnt+2,nct
+        write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
+     &      (rad2deg*thetatpl(ki,i),
      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
        enddo
        write (iout,*) "SC restraints from templates"
-       do i=nnt,lim_xx
+       do i=nnt,nct
         write(iout,'(i5,100(4f8.2,4x))') i,
      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index b3c4c22..5eca704 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -7060,11 +7060,12 @@ c      write (iout,*) idihconstr_start_homo,idihconstr_end_homo
       do i=idihconstr_start_homo,idihconstr_end_homo
         kat2=0.0d0
 c        betai=beta(i,i+1,i+2,i+3)
-        betai = phi(i+3)
+        betai = phi(i)
 c       write (iout,*) "betai =",betai
         do k=1,constr_homology
           dih_diff(k)=pinorm(dih(k,i)-betai)
-c         write (iout,*) "dih_diff(",k,") =",dih_diff(k)
+cd          write (iout,'(a8,2i4,2f15.8)') "dih_diff",i,k,dih_diff(k)
+cd     &                  ,sigma_dih(k,i)
 c          if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)=
 c     &                                   -(6.28318-dih_diff(i,k))
 c          if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)=
@@ -7172,6 +7173,9 @@ c
 c         dtheta_i=theta(j)-thetaref(j,iref)
 c         dtheta_i=thetaref(k,i)-theta(i) ! original form without indexing
           theta_diff(k)=thetatpl(k,i)-theta(i)
+cd          write (iout,'(a8,2i4,2f15.8)') "theta_diff",i,k,theta_diff(k)
+cd     &                  ,sigma_theta(k,i)
+
 c
           utheta_i=-0.5d0*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta rmvd from Gaussian argument
 c         utheta_i=-0.5d0*waga_theta*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta?
@@ -7242,7 +7246,7 @@ c                                    Original sign inverted for calc of gradient
           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
+cd          write(iout,'(2i5,4f8.2)') k,i,dxx,dyy,dzz,sigma_d(k,i)
 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?

diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F
index 4d92752..570600a 100644 (file)
--- a/source/unres/src_MD-M/initialize_p.F
+++ b/source/unres/src_MD-M/initialize_p.F
@@ -1487,10 +1487,10 @@ cd     &   " lim_dih",lim_dih
 #ifdef MPI
       if (me.eq.king .or. .not. out1file) write (iout,*) "MPI"
       call int_bounds(lim_odl,link_start_homo,link_end_homo)
-      call int_bounds(lim_dih-nnt+1,idihconstr_start_homo,
+      call int_bounds(lim_dih,idihconstr_start_homo,
      &  idihconstr_end_homo)
-      idihconstr_start_homo=idihconstr_start_homo+nnt-1
-      idihconstr_end_homo=idihconstr_end_homo+nnt-1
+      idihconstr_start_homo=idihconstr_start_homo+nnt-1+3
+      idihconstr_end_homo=idihconstr_end_homo+nnt-1+3
       if (me.eq.king .or. .not. out1file) 
      &  write (iout,*) 'Processor',fg_rank,' CG group',kolor,
      &  ' absolute rank',MyRank,
@@ -1502,8 +1502,8 @@ cd     &   " lim_dih",lim_dih
       write (iout,*) "Not MPI"
       link_start_homo=1
       link_end_homo=lim_odl
-      idihconstr_start_homo=nnt
-      idihconstr_end_homo=lim_dih
+      idihconstr_start_homo=nnt+3
+      idihconstr_end_homo=lim_dih+3
       write (iout,*) 
      &  ' lim_odl',lim_odl,' link_start=',link_start_homo,
      &  ' link_end',link_end_homo,' lim_dih',lim_dih,

diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F
index f944233..2a591cf 100644 (file)
--- a/source/unres/src_MD-M/readrtns_CSA.F
+++ b/source/unres/src_MD-M/readrtns_CSA.F
@@ -2542,6 +2542,7 @@ c-------------------------------------------------------------------------------
       include 'COMMON.MD'
       include 'COMMON.GEO'
       include 'COMMON.INTERACT'
+      include 'COMMON.NAMES'
 c
 c For new homol impl
 c
@@ -2562,9 +2563,11 @@ c    &    sigma_odl_temp(maxres,maxres,max_template)
 c
 c     FP - Nov. 2014 Temporary specifications for new vars
 c
-      double precision rescore_tmp,x12,y12,z12,rescore2_tmp
+      double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
+     &                 rescore3_tmp
       double precision, dimension (max_template,maxres) :: rescore
       double precision, dimension (max_template,maxres) :: rescore2
+      double precision, dimension (max_template,maxres) :: rescore3
       character*24 pdbfile,tpl_k_rescore
 c -----------------------------------------------------------------
 c Reading multiple PDB ref structures and calculation of retraints
@@ -2615,11 +2618,6 @@ cd      call flush(iout)
       lim_odl=0
       lim_dih=0
 c
-c  New
-c
-      lim_theta=0
-      lim_xx=0
-c
       write(iout,*) 'nnt=',nnt,'nct=',nct
 c
       do i = nnt,nct
@@ -2683,14 +2681,18 @@ c
 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
           open (ientin,file=tpl_k_rescore,status='old')
           if (nnt.gt.1) rescore(k,1)=0.0d0
-          do irec=nnt,maxdim ! loop for reading res sim 
+          do irec=nnt,nct ! loop for reading res sim 
             if (read2sigma) then
              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
-     &                                idomain_tmp
+     &                                rescore3_tmp,idomain_tmp
              i_tmp=i_tmp+nnt-1
              idomain(k,i_tmp)=idomain_tmp
              rescore(k,i_tmp)=rescore_tmp
              rescore2(k,i_tmp)=rescore2_tmp
+             rescore3(k,i_tmp)=rescore3_tmp
+             write(iout,'(a7,i5,3f10.5,i5)') "rescore",
+     &                      i_tmp,rescore2_tmp,rescore_tmp,
+     &                                rescore3_tmp,idomain_tmp
             else
              idomain(k,irec)=1
              read (ientin,*,end=1401) rescore_tmp
@@ -2790,7 +2792,8 @@ c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
 c   Instead of res sim other local measure of b/b str reliability possible
-            sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
+            if (sigma_dih(k,i).ne.0)
+     &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
           enddo
           lim_dih=nct-nnt-2 
@@ -2822,14 +2825,14 @@ c            read (ientin,*) sigma_theta(k,i) ! 1st variant
              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
      &                        rescore(k,i-2))/3.0
 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
-             sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
+             if (sigma_theta(k,i).ne.0)
+     &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 
 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 c                             rescore(k,i-2) !  right expression ?
 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
           enddo
         endif
-        lim_theta=nct-nnt-1 
 
         if (waga_d.gt.0.0d0) then
 c       open (ientin,file=tpl_k_sigma_d,status='old')
@@ -2855,15 +2858,14 @@ c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
-               sigma_d(k,i)=rescore(k,i) !  right expression ?
-               sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
+               sigma_d(k,i)=rescore3(k,i) !  right expression ?
+               if (sigma_d(k,i).ne.0)
+     &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
 
 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
 c              read (ientin,*) sigma_d(k,i) ! 1st variant
-               if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
           enddo
-          lim_xx=nct-nnt+1 
         endif
       enddo
 c
@@ -2894,8 +2896,6 @@ c
       endif
       if (constr_homology.gt.0) call homology_partition
       if (constr_homology.gt.0) call init_int_table
-cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
-cd     & "lim_xx=",lim_xx
 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
 c
@@ -2912,17 +2912,19 @@ cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
      &  ki=1,constr_homology)
        enddo
        write (iout,*) "Dihedral angle restraints from templates"
-       do i=nnt+3,lim_dih
-        write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*dih(ki,i),
+       do i=nnt+3,nct
+        write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
+     &      (rad2deg*dih(ki,i),
      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
        enddo
        write (iout,*) "Virtual-bond angle restraints from templates"
-       do i=nnt+2,lim_theta
-        write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
+       do i=nnt+2,nct
+        write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
+     &      (rad2deg*thetatpl(ki,i),
      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
        enddo
        write (iout,*) "SC restraints from templates"
-       do i=nnt,lim_xx
+       do i=nnt,nct
         write(iout,'(i5,100(4f8.2,4x))') i,
      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)

diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F
index d0775ad..a98fe50 100644 (file)
--- a/source/wham/src-M/energy_p_new.F
+++ b/source/wham/src-M/energy_p_new.F
@@ -3673,7 +3673,7 @@ c      write (iout,*) idihconstr_start_homo,idihconstr_end_homo
       do i=idihconstr_start_homo,idihconstr_end_homo
         kat2=0.0d0
 c        betai=beta(i,i+1,i+2,i+3)
-        betai = phi(i+3)
+        betai = phi(i)
 c       write (iout,*) "betai =",betai
         do k=1,constr_homology
           dih_diff(k)=pinorm(dih(k,i)-betai)

diff --git a/source/wham/src-M/initialize_p.F b/source/wham/src-M/initialize_p.F
index 459c312..3f5f299 100644 (file)
--- a/source/wham/src-M/initialize_p.F
+++ b/source/wham/src-M/initialize_p.F
@@ -632,10 +632,10 @@ c      include 'COMMON.SETUP'
      &   " lim_dih",lim_dih
 #ifdef MPL
       call int_bounds(lim_odl,link_start_homo,link_end_homo)
-      call int_bounds(lim_dih-nnt+1,idihconstr_start_homo,
+      call int_bounds(lim_dih,idihconstr_start_homo,
      &  idihconstr_end_homo)
-      idihconstr_start_homo=idihconstr_start_homo+nnt-1
-      idihconstr_end_homo=idihconstr_end_homo+nnt-1
+      idihconstr_start_homo=idihconstr_start_homo+nnt-1+3
+      idihconstr_end_homo=idihconstr_end_homo+nnt-1+3
       if (me.eq.king .or. .not. out1file)
      &  write (iout,*) 'Processor',fg_rank,' CG group',kolor,
      &  ' absolute rank',MyRank,
@@ -646,8 +646,8 @@ c      include 'COMMON.SETUP'
 #else
       link_start_homo=1
       link_end_homo=lim_odl
-      idihconstr_start_homo=nnt
-      idihconstr_end_homo=lim_dih
+      idihconstr_start_homo=nnt+3
+      idihconstr_end_homo=lim_dih+3
       write (iout,*)
      &  ' lim_odl',lim_odl,' link_start=',link_start_homo,
      &  ' link_end',link_end_homo,' lim_dih',lim_dih,

diff --git a/source/wham/src-M/molread_zs.F b/source/wham/src-M/molread_zs.F
index dee1afc..0c1a211 100644 (file)
--- a/source/wham/src-M/molread_zs.F
+++ b/source/wham/src-M/molread_zs.F
@@ -339,6 +339,7 @@ c====-------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.GEO'
       include 'COMMON.INTERACT'
+      include 'include_unres/COMMON.NAMES'
       include 'COMMON.HOMRESTR'
 c
 c For new homol impl
@@ -363,9 +364,11 @@ c    &    sigma_odl_temp(maxres,maxres,max_template)
 c
 c     FP - Nov. 2014 Temporary specifications for new vars
 c
-      double precision rescore_tmp,x12,y12,z12,rescore2_tmp
+      double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
+     &                 rescore3_tmp
       double precision, dimension (max_template,maxres) :: rescore
       double precision, dimension (max_template,maxres) :: rescore2
+      double precision, dimension (max_template,maxres) :: rescore3
       character*24 tpl_k_rescore
 c -----------------------------------------------------------------
 c Reading multiple PDB ref structures and calculation of retraints
@@ -410,8 +413,6 @@ cd      call flush(iout)
 c
 c  New
 c
-      lim_theta=0
-      lim_xx=0
 c
 c  Reading HM global scores (prob not required)
 c
@@ -494,14 +495,18 @@ c
 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
           open (ientin,file=tpl_k_rescore,status='old')
           if (nnt.gt.1) rescore(k,1)=0.0d0
-          do irec=nnt,maxdim ! loop for reading res sim 
+          do irec=nnt,nct ! loop for reading res sim 
             if (read2sigma) then
              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
-     &                                idomain_tmp
+     &                                rescore3_tmp,idomain_tmp
              i_tmp=i_tmp+nnt-1
              idomain(k,i_tmp)=idomain_tmp
              rescore(k,i_tmp)=rescore_tmp
              rescore2(k,i_tmp)=rescore2_tmp
+             rescore3(k,i_tmp)=rescore3_tmp
+             write(iout,'(a7,i5,3f10.5,i5)') "rescore",
+     &                      i_tmp,rescore2_tmp,rescore_tmp,
+     &                                rescore3_tmp,idomain_tmp
             else
              idomain(k,irec)=1
              read (ientin,*,end=1401) rescore_tmp
@@ -601,7 +606,8 @@ c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
 c   Instead of res sim other local measure of b/b str reliability possible
-            sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
+            if (sigma_dih(k,i).ne.0)
+     &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
           enddo
           lim_dih=nct-nnt-2 
@@ -633,14 +639,14 @@ c            read (ientin,*) sigma_theta(k,i) ! 1st variant
              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
      &                        rescore(k,i-2))/3.0
 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
-             sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
+             if (sigma_theta(k,i).ne.0)
+     &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 
 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 c                             rescore(k,i-2) !  right expression ?
 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
           enddo
         endif
-        lim_theta=nct-nnt-1 
 
         if (waga_d.gt.0.0d0) then
 c       open (ientin,file=tpl_k_sigma_d,status='old')
@@ -666,15 +672,14 @@ c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
-               sigma_d(k,i)=rescore(k,i) !  right expression ?
-               sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
+               sigma_d(k,i)=rescore3(k,i) !  right expression ?
+               if (sigma_d(k,i).ne.0)
+     &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
 
 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
 c              read (ientin,*) sigma_d(k,i) ! 1st variant
-               if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
           enddo
-          lim_xx=nct-nnt+1 
         endif
       enddo
 c
@@ -723,17 +728,19 @@ cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
      &  ki=1,constr_homology)
        enddo
        write (iout,*) "Dihedral angle restraints from templates"
-       do i=nnt+3,lim_dih
-        write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*dih(ki,i),
+       do i=nnt+3,nct
+        write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
+     &      (rad2deg*dih(ki,i),
      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
        enddo
        write (iout,*) "Virtual-bond angle restraints from templates"
-       do i=nnt+2,lim_theta
-        write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
+       do i=nnt+2,nct
+        write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
+     &      (rad2deg*thetatpl(ki,i),
      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
        enddo
        write (iout,*) "SC restraints from templates"
-       do i=nnt,lim_xx
+       do i=nnt,nct
         write(iout,'(i5,100(4f8.2,4x))') i,
      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)