1st running version of UNRES HM by FP and AL
[unres.git] / source / unres / src_MD / readrtns.F
index 10cde78..f1cc85b 100644 (file)
@@ -96,6 +96,7 @@ c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
 C Set up the time limit (caution! The time must be input in minutes!)
       read_cart=index(controlcard,'READ_CART').gt.0
       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
+      call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
@@ -372,7 +373,6 @@ C
       endif
       call reada(controlcard,"Q_NP",Q_np,0.1d0)
       usampl = index(controlcard,"USAMPL").gt.0
-
       mdpdb = index(controlcard,"MDPDB").gt.0
       call reada(controlcard,"T_BATH",t_bath,300.0d0)
       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1) 
@@ -749,6 +749,10 @@ C 12/1/95 Added weight for the multi-body term WCORR
        call reada(weightcard,'WTORD',wtor_d,1.0D0)
        call reada(weightcard,'WANG',wang,1.0D0)
        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
+       call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
+       call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
+       call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
+       call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
        call reada(weightcard,'SCAL14',scal14,0.4D0)
        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
@@ -778,11 +782,16 @@ C 12/1/95 Added weight for the multi-body term WCORR
        weights(18)=scal14
        weights(21)=wsccor
       endif
+       weights(25)=wdfa_dist
+       weights(26)=wdfa_tor
+       weights(27)=wdfa_nei
+       weights(28)=wdfa_beta
 
       if(me.eq.king.or..not.out1file)
      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
-     &  wturn4,wturn6
+     &  wturn4,wturn6,
+     &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
    10 format (/'Energy-term weights (unscaled):'//
      & 'WSCC=   ',f10.6,' (SC-SC)'/
      & 'WSCP=   ',f10.6,' (SC-p)'/
@@ -801,7 +810,11 @@ C 12/1/95 Added weight for the multi-body term WCORR
      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
-     & 'WTURN6= ',f10.6,' (turns, 6th order)')
+     & 'WTURN6= ',f10.6,' (turns, 6th order)'/
+     & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
+     & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
+     & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
+     & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
       if(me.eq.king.or..not.out1file)then
        if (wcorr4.gt.0.0d0) then
         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
@@ -829,7 +842,8 @@ C 12/1/95 Added weight for the multi-body term WCORR
       if(me.eq.king.or..not.out1file)
      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
-     &  wturn4,wturn6
+     &  wturn4,wturn6,
+     &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
    22 format (/'Energy-term weights (scaled):'//
      & 'WSCC=   ',f10.6,' (SC-SC)'/
      & 'WSCP=   ',f10.6,' (SC-p)'/
@@ -848,7 +862,11 @@ C 12/1/95 Added weight for the multi-body term WCORR
      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
-     & 'WTURN6= ',f10.6,' (turns, 6th order)')
+     & 'WTURN6= ',f10.6,' (turns, 6th order)'/
+     & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
+     & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
+     & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
+     & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
       if(me.eq.king.or..not.out1file)
      & write (iout,*) "Reference temperature for weights calculation:",
      &  temp0
@@ -940,8 +958,8 @@ C 10/03/12 Adam: Recalculate coordinates with new side chain positions
          call chainbuild
         endif  
 C Following 2 lines for diagnostics; comment out if not needed
-        write (iout,*) "After sideadd"
-        call intout
+c        write (iout,*) "After sideadd"
+c        call intout
       endif
       if (indpdb.eq.0) then
 C Read sequence if not taken from the pdb file.
@@ -1039,6 +1057,21 @@ C 8/13/98 Set limits to generating the dihedral angles
 cd      print *,'NNT=',NNT,' NCT=',NCT
       if (itype(1).eq.21) nnt=2
       if (itype(nres).eq.21) nct=nct-1
+
+C     Bartek:READ init_vars
+C     Initialize variables!
+C     Juyong:READ read_info
+C     READ fragment information!!
+C     both routines should be in dfa.F file!!
+
+      if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
+     &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
+       call init_dfa_vars
+       print*, 'init_dfa_vars finished!'
+       call read_dfa_info
+       print*, 'read_dfa_info finished!'
+      endif
+C
       if (pdbref) then
         if(me.eq.king.or..not.out1file)
      &   write (iout,'(a,i3)') 'nsup=',nsup
@@ -1128,6 +1161,13 @@ c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
       if (constr_dist.gt.0) then
         call read_dist_constr
       endif
+
+
+      if (constr_homology.gt.0) then
+        call read_constr_homology
+      endif
+
+
       if (nhpb.gt.0) call hpb_partition
 c      write (iout,*) "After read_dist_constr nhpb",nhpb
 c      call flush(iout)
@@ -1251,18 +1291,6 @@ C Generate distance constraints, if the PDB structure is to be regularized.
         write (iout,'(20i4)') (iss(i),i=1,ns)
        if (dyn_ss) then
           write(iout,*)"Running with dynamic disulfide-bond formation"
-          do i=nss+1,nhpb
-            ihpb(i-nss)=ihpb(i)
-            jhpb(i-nss)=jhpb(i)
-            forcon(i-nss)=forcon(i)
-            dhpb(i-nss)=dhpb(i)
-          enddo
-          nhpb=nhpb-nss
-          nss=0
-          call hpb_partition
-          do i=1,ns
-            dyn_ss_mask(iss(i))=.true.
-          enddo
        else
         write (iout,'(/a/)') 'Pre-formed links are:' 
        do i=1,nss
@@ -1270,14 +1298,27 @@ C Generate distance constraints, if the PDB structure is to be regularized.
          i2=jhpb(i)-nres
          it1=itype(i1)
          it2=itype(i2)
-         if (me.eq.king.or..not.out1file)
-     &    write (iout,'(2a,i3,3a,i3,a,3f10.3)')
+          write (iout,'(2a,i3,3a,i3,a,3f10.3)')
      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
      &    ebr,forcon(i)
        enddo
        write (iout,'(a)')
        endif
       endif
+      if (ns.gt.0.and.dyn_ss) then
+          do i=nss+1,nhpb
+            ihpb(i-nss)=ihpb(i)
+            jhpb(i-nss)=jhpb(i)
+            forcon(i-nss)=forcon(i)
+            dhpb(i-nss)=dhpb(i)
+          enddo
+          nhpb=nhpb-nss
+          nss=0
+          call hpb_partition
+          do i=1,ns
+            dyn_ss_mask(iss(i))=.true.
+          enddo
+      endif
       if (i2ndstr.gt.0) call secstrp2dihc
 c      call geom_to_var(nvar,x)
 c      call etotal(energia(0))
@@ -1356,7 +1397,8 @@ C Check whether the specified bridging residues are cystines.
 C Read preformed bridges.
       if (ns.gt.0) then
       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
-      write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
+      if(fg_rank.eq.0)
+     & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
       if (nss.gt.0) then
         nhpb=nss
 C Check if the residues involved in bridges are in the specified list of
@@ -2484,7 +2526,7 @@ c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
         if (wfrag_(i).gt.0.0d0) then
         do j=ifrag_(1,i),ifrag_(2,i)-1
           do k=j+1,ifrag_(2,i)
-            write (iout,*) "j",j," k",k
+c            write (iout,*) "j",j," k",k
             ddjk=dist(j,k)
             if (constr_dist.eq.1) then
             nhpb=nhpb+1
@@ -2574,6 +2616,343 @@ c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
       return
       end
 c-------------------------------------------------------------------------------
+
+      subroutine read_constr_homology
+
+      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+      include 'COMMON.SETUP'
+      include 'COMMON.CONTROL'
+      include 'COMMON.CHAIN'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.MD'
+      include 'COMMON.GEO'
+      include 'COMMON.INTERACT'
+c
+c For new homol impl
+c
+      include 'COMMON.VAR'
+c
+
+c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
+c    &                 dist_cut
+c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
+c    &    sigma_odl_temp(maxres,maxres,max_template)
+      character*2 kic2
+      character*24 model_ki_dist, model_ki_angle
+      character*500 controlcard
+      integer ki, i, j, k, l
+      logical lprn /.true./
+c
+c     FP - Nov. 2014 Temporary specifications for new vars
+c
+      double precision rescore_tmp,x12,y12,z12
+      double precision, dimension (max_template,maxres) :: rescore
+      character*24 pdbfile,tpl_k_rescore
+c -----------------------------------------------------------------
+c Reading multiple PDB ref structures and calculation of retraints
+c not using pre-computed ones stored in files model_ki_{dist,angle}
+c FP (Nov., 2014)
+c -----------------------------------------------------------------
+c
+c
+c Alternative: reading from input
+      call card_concat(controlcard)
+      call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
+      call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
+      call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
+      call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
+      call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
+      lim_odl=0
+      lim_dih=0
+c
+c  New
+c
+      lim_theta=0
+      lim_xx=0
+c
+c  Reading HM global scores (prob not required)
+c
+c     open (4,file="HMscore")
+c     do k=1,constr_homology
+c       read (4,*,end=521) hmscore_tmp
+c       hmscore(k)=hmscore_tmp ! Another transformation can be used 
+c       write(*,*) "Model", k, ":", hmscore(k)
+c     enddo
+c521  continue
+
+c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
+      do k=1,constr_homology
+
+        read(inp,'(a)') pdbfile
+c  Next stament causes error upon compilation (?)
+c       if(me.eq.king.or. .not. out1file)
+c         write (iout,'(2a)') 'PDB data will be read from file ',
+c    &   pdbfile(:ilen(pdbfile))
+        open(ipdbin,file=pdbfile,status='old',err=33)
+        goto 34
+  33    write (iout,'(a)') 'Error opening PDB file.'
+        stop
+  34    continue
+c        print *,'Begin reading pdb data'
+c
+c Files containing res sim or local scores (former containing sigmas)
+c
+
+        write(kic2,'(bz,i2.2)') k
+
+        tpl_k_rescore="template"//kic2//".sco"
+c       tpl_k_sigma_odl="template"//kic2//".sigma_odl"
+c       tpl_k_sigma_dih="template"//kic2//".sigma_dih"
+c       tpl_k_sigma_theta="template"//kic2//".sigma_theta"
+c       tpl_k_sigma_d="template"//kic2//".sigma_d"
+
+        call readpdb
+c
+c     Distance restraints
+c
+c          ... --> odl(k,ii)
+C Copy the coordinates from reference coordinates (?)
+        do i=1,2*nres
+          do j=1,3
+            c(j,i)=cref(j,i)
+c           write (iout,*) "c(",j,i,") =",c(j,i)
+          enddo
+        enddo
+c
+c From read_dist_constr (commented out 25/11/2014 <-> res sim)
+c
+c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
+          open (ientin,file=tpl_k_rescore,status='old')
+          do irec=1,maxdim ! loop for reading res sim 
+            if (irec.eq.1) then
+               rescore(k,irec)=0.0d0
+               goto 1301 
+            endif
+            read (ientin,*,end=1401) rescore_tmp
+c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
+            rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
+c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
+ 1301     continue
+          enddo  
+ 1401   continue
+          close (ientin)        
+c         open (ientin,file=tpl_k_sigma_odl,status='old')
+c         do irec=1,maxdim ! loop for reading sigma_odl
+c            read (ientin,*,end=1401) i, j, 
+c    &                                sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
+c            sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
+c    &       sigma_odl_temp(i+nnt-1,j+nnt-1,k) 
+c         enddo
+c 1401   continue
+c         close (ientin)
+        if (waga_dist.gt.0.0d0) then
+          ii=0
+          do i = nnt,nct-2 ! right? without parallel.
+            do j=i+2,nct ! right?
+c         do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology 
+c           do j=i+2,nres ! ibid
+c         do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
+c           do j=i+2,nct ! ibid
+              ii=ii+1
+c             write (iout,*) "k",k
+c             write (iout,*) "i",i," j",j," constr_homology",
+c    &                       constr_homology
+              ires_homo(ii)=i
+              jres_homo(ii)=j
+c
+c Attempt to replace dist(i,j) by its definition in ...
+c
+              x12=c(1,i)-c(1,j)
+              y12=c(2,i)-c(2,j)
+              z12=c(3,i)-c(3,j)
+              distal=dsqrt(x12*x12+y12*y12+z12*z12)
+              odl(k,ii)=distal
+c
+c             odl(k,ii)=dist(i,j)
+c             write (iout,*) "dist(",i,j,") =",dist(i,j)
+c             write (iout,*) "distal = ",distal
+c             write (iout,*) "odl(",k,ii,") =",odl(k,ii)
+c            write(iout,*) "rescore(",k,i,") =",rescore(k,i),
+c    &                      "rescore(",k,j,") =",rescore(k,j)
+c
+c  Calculation of sigma from res sim
+c
+c             if (odl(k,ii).le.6.0d0) then
+c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
+c    Other functional forms possible depending on odl(k,ii), eg.
+c
+            if (odl(k,ii).le.dist_cut) then
+              sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
+c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
+            else
+              sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error 
+     &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
+
+c   Following expr replaced by a positive exp argument
+c             sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
+c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
+
+c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
+c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
+            endif
+c
+              sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
+c             sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
+c
+c             sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
+c    &                        sigma_odl_temp(i,j,k)  ! not inverse because of use of res. similarity
+            enddo
+c           read (ientin,*) sigma_odl(k,ii) ! 1st variant
+          enddo
+c         lim_odl=ii
+c         if (constr_homology.gt.0) call homology_partition
+        endif
+c
+c     Theta, dihedral and SC retraints
+c
+        if (waga_angle.gt.0.0d0) then
+c         open (ientin,file=tpl_k_sigma_dih,status='old')
+c         do irec=1,maxres-3 ! loop for reading sigma_dih
+c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
+c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
+c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
+c    &                            sigma_dih(k,i+nnt-1)
+c         enddo
+c1402   continue
+c         close (ientin)
+          do i = nnt+3,nct ! right? without parallel.
+c         do i=1,nres ! alternative for bounds acc to readpdb?
+c         do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
+c         do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
+            dih(k,i)=phiref(i) ! right?
+c           read (ientin,*) sigma_dih(k,i) ! original variant
+c             write (iout,*) "dih(",k,i,") =",dih(k,i)
+c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
+c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
+c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
+c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
+
+            sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
+     &                     rescore(k,i-2)+rescore(k,i-3)  !  right expression ?
+c
+c           write (iout,*) "Raw sigmas for dihedral angle restraints"
+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))
+c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
+            if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
+c           if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
+          enddo
+        endif
+
+        if (waga_theta.gt.0.0d0) then
+c         open (ientin,file=tpl_k_sigma_theta,status='old')
+c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
+c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
+c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
+c    &                              sigma_theta(k,i+nnt-1)
+c         enddo
+c1403   continue
+c         close (ientin)
+
+          do i = nnt+2,nct ! right? without parallel.
+c         do i = i=1,nres ! alternative for bounds acc to readpdb?
+c         do i=ithet_start,ithet_end ! with FG parallel.
+             thetatpl(k,i)=thetaref(i)
+c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
+c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
+c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
+c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
+c            read (ientin,*) sigma_theta(k,i) ! 1st variant
+             sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
+     &                        rescore(k,i-2) !  right expression ?
+             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)
+             if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
+          enddo
+        endif
+
+        if (waga_d.gt.0.0d0) then
+c       open (ientin,file=tpl_k_sigma_d,status='old')
+c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
+c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
+c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
+c    &                          sigma_d(k,i+nnt-1)
+c         enddo
+c1404   continue
+          close (ientin)
+
+          do i = nnt,nct ! right? without parallel.
+c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
+c         do i=loc_start,loc_end ! with FG parallel.
+             if (itype(i).eq.10) goto 1 ! right?
+               xxtpl(k,i)=xxref(i)
+               yytpl(k,i)=yyref(i)
+               zztpl(k,i)=zzref(i)
+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))
+
+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?
+    1     continue
+          enddo
+        endif
+        close(ientin)
+      enddo
+      if (waga_dist.gt.0.0d0) lim_odl=ii
+      if (constr_homology.gt.0) call homology_partition
+      if (constr_homology.gt.0) call init_int_table
+      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
+     & "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
+c Print restraints
+c
+      if (.not.lprn) return
+      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
+      write (iout,*) "Distance restraints from templates"
+      do ii=1,lim_odl
+       write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
+     &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
+      enddo
+      write (iout,*) "Dihedral angle restraints from templates"
+      do i=nnt+3,lim_dih
+        write (iout,'(i5,10(2f8.2,4x))') 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,10(2f8.2,4x))') 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
+        write(iout,'(i5,10(4f8.2,4x))') i,
+     &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
+     &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
+      enddo
+
+c -----------------------------------------------------------------
+      return
+      end
+c----------------------------------------------------------------------
+
 #ifdef WINIFL
       subroutine flush(iu)
       return
@@ -2585,6 +2964,7 @@ c-------------------------------------------------------------------------------
       return
       end
 #endif
+
 c------------------------------------------------------------------------------
       subroutine copy_to_tmp(source)
       include "DIMENSIONS"