correction of dihedral homology restraints phi(i+3)

[unres.git] / source / unres / src_MD-M / readrtns_CSA.F

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)