SAXS cutoff correct gradiend and cutoff scaling
[unres.git] / source / unres / src_MD-M / readrtns_CSA.F
index 7be445d..70bb647 100644 (file)
@@ -98,6 +98,13 @@ 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)
+      write (iout,*) "constr_dist",constr_dist
+      call readi(controlcard,'NSAXS',nsaxs,0)
+      call readi(controlcard,'SAXS_MODE',saxs_mode,0)
+      call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
+      call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
+      write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
+     &   SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
       call readi(controlcard,'SYM',symetr,1)
       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
@@ -241,10 +248,10 @@ C      endif
       if ((lipbufthick*2.0d0).gt.lipthick)
      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
       endif
-      write(iout,*) "bordliptop=",bordliptop
-      write(iout,*) "bordlipbot=",bordlipbot
-      write(iout,*) "bufliptop=",bufliptop
-      write(iout,*) "buflipbot=",buflipbot
+c      write(iout,*) "bordliptop=",bordliptop
+c      write(iout,*) "bordlipbot=",bordlipbot
+c      write(iout,*) "bufliptop=",bufliptop
+c      write(iout,*) "buflipbot=",buflipbot
 
 
       if (me.eq.king .or. .not.out1file ) 
@@ -619,6 +626,7 @@ C Read weights of the subsequent energy terms.
        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
        call reada(weightcard,'TEMP0',temp0,300.0d0)
        call reada(weightcard,'WLT',wliptran,0.0D0)
+       call reada(weightcard,'WSAXS',wsaxs,1.0D0)
        if (index(weightcard,'SOFT').gt.0) ipot=6
 C 12/1/95 Added weight for the multi-body term WCORR
        call reada(weightcard,'WCORRH',wcorr,1.0D0)
@@ -642,6 +650,7 @@ C 12/1/95 Added weight for the multi-body term WCORR
        weights(17)=wbond
        weights(18)=scal14
        weights(21)=wsccor
+       weights(25)=wsaxs
       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,
@@ -803,12 +812,16 @@ c        print *,'Finished reading pdb data'
               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
               nsi=nsi+1
             enddo
+c AL 7/10/16
+c Calculalte only the coordinates of the current sidechain; no need to rebuild
+c whole chain
+            call locate_side_chain(i)
             if(fail) write(iout,*)'Adding sidechain failed for res ',
      &              i,' after ',nsi,' trials'
           endif
          enddo
 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
-         call chainbuild
+c         call chainbuild
         endif  
       endif
       if (indpdb.eq.0) then
@@ -869,6 +882,7 @@ C 8/13/98 Set limits to generating the dihedral angles
         phibound(2,i)=pi
       enddo
       read (inp,*) ndih_constr
+      write (iout,*) "ndish_constr",ndih_constr
       if (ndih_constr.gt.0) then
         read (inp,*) ftors
         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
@@ -963,7 +977,7 @@ c----------------------
           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
           stop 'Error reading reference structure'
 #endif
-   39     call chainbuild
+   39     call chainbuild_extconf
           call setup_var
 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
           nstart_sup=nnt
@@ -980,7 +994,7 @@ czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
         call flush(iout)
         if (constr_dist.gt.0) call read_dist_constr
-        write (iout,*) "After read_dist_constr nhpb",nhpb
+c        write (iout,*) "After read_dist_constr nhpb",nhpb
         if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
         if(me.eq.king.or..not.out1file)
      &   write (iout,*) 'Contact order:',co
@@ -998,6 +1012,8 @@ c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
         enddo
         endif
       endif
+      write (iout,*) "calling read_saxs_consrtr",nsaxs
+      if (nsaxs.gt.0) call read_saxs_constr
 
 
       if (constr_homology.gt.0) then
@@ -1082,7 +1098,10 @@ C initial geometry.
             enddo
             return
           else
+            write (iout,*) "Calling read_ang"
             call read_angles(inp,*36)
+            write (iout,*) "Calling chainbuild"
+            call chainbuild_extconf
           endif
           goto 37
    36     write (iout,'(a)') 'Error reading angle file.'
@@ -1107,6 +1126,53 @@ C initial geometry.
           omeg(i)=-120d0*deg2rad
           if (itype(i).le.0) omeg(i)=-omeg(i)
          enddo
+c from old chainbuild
+C
+C Define the origin and orientation of the coordinate system and locate the
+C first three CA's and SC(2).
+C
+      call orig_frame
+*
+* Build the alpha-carbon chain.
+*
+      do i=4,nres
+       call locate_next_res(i)
+      enddo     
+C
+C First and last SC must coincide with the corresponding CA.
+C
+      do j=1,3
+       dc(j,nres+1)=0.0D0
+        dc_norm(j,nres+1)=0.0D0
+       dc(j,nres+nres)=0.0D0
+        dc_norm(j,nres+nres)=0.0D0
+        c(j,nres+1)=c(j,1)
+        c(j,nres+nres)=c(j,nres)
+      enddo
+C
+C Define the origin and orientation of the coordinate system and locate the
+C first three CA's and SC(2).
+C
+      call orig_frame
+*
+* Build the alpha-carbon chain.
+*
+      do i=4,nres
+       call locate_next_res(i)
+      enddo     
+C
+C First and last SC must coincide with the corresponding CA.
+C
+      do j=1,3
+       dc(j,nres+1)=0.0D0
+        dc_norm(j,nres+1)=0.0D0
+       dc(j,nres+nres)=0.0D0
+        dc_norm(j,nres+nres)=0.0D0
+        c(j,nres+1)=c(j,1)
+        c(j,nres+nres)=c(j,nres)
+      enddo
+
+c
         else
           if(me.eq.king.or..not.out1file)
      &     write (iout,'(a)') 'Random-generated initial geometry.'
@@ -1782,6 +1848,8 @@ c----------------------------------------------------------------------------
       include 'DIMENSIONS'
       include 'COMMON.MINIM'
       include 'COMMON.IOUNITS'
+      include 'COMMON.CONTROL'
+      include 'COMMON.SETUP'
       character*80 ucase
       character*320 minimcard
       call card_concat(minimcard)
@@ -1794,12 +1862,18 @@ c----------------------------------------------------------------------------
       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
+#ifdef MPI
+      if (.not. out1file .or. me.eq.king) then
+#endif
       write (iout,'(/80(1h*)/20x,a/80(1h*))') 
      &         'Options in energy minimization:'
       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
      & 'MinMin:',MinMin,' MinFun:',MinFun,
      & ' TolF:',TolF,' RTolF:',RTolF
+#ifdef MPI
+      endif
+#endif
       return
       end
 c----------------------------------------------------------------------------
@@ -2062,6 +2136,8 @@ C Get parameter filenames and open the parameter files.
       open (ielep,file=elename,status='old')
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old')
+      call getenv_loc('LIPTRANPAR',liptranname)
+      open (iliptranpar,file=liptranname,status='old')
 #else
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
      &  readonly)
@@ -2093,7 +2169,7 @@ c      print *,ithep_pdb," opened"
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old',readonly)
       call getenv_loc('LIPTRANPAR',liptranname)
-      open (iliptranpar,file=liptranname,status='old',action='read')
+      open (iliptranpar,file=liptranname,status='old',readonly)
 #ifndef CRYST_SC
       call getenv_loc('ROTPARPDB',rotname_pdb)
       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
@@ -2389,6 +2465,67 @@ CCCC NOW PROPERTIES FOR AFM
       end
 
 c-------------------------------------------------------------------------------
+      subroutine read_saxs_constr
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+      include 'COMMON.SETUP'
+      include 'COMMON.CONTROL'
+      include 'COMMON.CHAIN'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.SBRIDGE'
+      double precision cm(3)
+c      read(inp,*) nsaxs
+      write (iout,*) "Calling read_saxs nsaxs",nsaxs
+      call flush(iout)
+      if (saxs_mode.eq.0) then
+c SAXS distance distribution
+      do i=1,nsaxs
+        read(inp,*) distsaxs(i),Psaxs(i)
+      enddo
+      Cnorm = 0.0d0
+      do i=1,nsaxs
+        Cnorm = Cnorm + Psaxs(i)
+      enddo
+      write (iout,*) "Cnorm",Cnorm
+      do i=1,nsaxs
+        Psaxs(i)=Psaxs(i)/Cnorm
+      enddo
+      write (iout,*) "Normalized distance distribution from SAXS"
+      do i=1,nsaxs
+        write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
+      enddo
+      else
+c SAXS "spheres".
+      do i=1,nsaxs
+        read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
+      enddo
+      do j=1,3
+        cm(j)=0.0d0
+      enddo
+      do i=1,nsaxs
+        do j=1,3
+          cm(j)=cm(j)+Csaxs(j,i)
+        enddo
+      enddo
+      do j=1,3
+        cm(j)=cm(j)/nsaxs
+      enddo
+      do i=1,nsaxs
+        do j=1,3
+          Csaxs(j,i)=Csaxs(j,i)-cm(j)
+        enddo
+      enddo
+      write (iout,*) "SAXS sphere coordinates"
+      do i=1,nsaxs
+        write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
+      enddo
+      endif
+      return
+      end
+c-------------------------------------------------------------------------------
       subroutine read_dist_constr
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
@@ -2542,6 +2679,7 @@ c-------------------------------------------------------------------------------
       include 'COMMON.MD'
       include 'COMMON.GEO'
       include 'COMMON.INTERACT'
+      include 'COMMON.NAMES'
 c
 c For new homol impl
 c
@@ -2559,12 +2697,15 @@ c    &    sigma_odl_temp(maxres,maxres,max_template)
       logical lprn /.true./
       integer ilen
       external ilen
+      logical liiflag
 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
@@ -2580,9 +2721,23 @@ c Alternative: reading from input
       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
+      call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
       call readi(controlcard,"HOMOL_NSET",homol_nset,1)       
       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
+      start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
+      if(.not.read2sigma.and.start_from_model) then
+          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) 
+     &      write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
+          start_from_model=.false.
+      endif
+      if(start_from_model .and. (me.eq.king .or. .not. out1file))
+     &    write(iout,*) 'START_FROM_MODELS is ON'
+      if(start_from_model .and. rest) then 
+        if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+         write(iout,*) 'START_FROM_MODELS is OFF'
+         write(iout,*) 'remove restart keyword from input'
+        endif
+      endif
       if (homol_nset.gt.1)then
          call card_concat(controlcard)
          read(controlcard,*) (waga_homology(i),i=1,homol_nset) 
@@ -2605,12 +2760,7 @@ 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      write(iout,*) 'nnt=',nnt,'nct=',nct
 c
       do i = nnt,nct
         do k=1,constr_homology
@@ -2629,11 +2779,8 @@ c
       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))
-         write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
+        if(me.eq.king .or. .not. out1file)
+     &    write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
      &  pdbfile(:ilen(pdbfile))
         open(ipdbin,file=pdbfile,status='old',err=33)
         goto 34
 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)=0.5d0*(rescore_tmp+0.5d0)
-             rescore2(k,i_tmp)=0.5d0*(rescore2_tmp+0.5d0)
+             rescore(k,i_tmp)=rescore_tmp
+             rescore2(k,i_tmp)=rescore2_tmp
+             rescore3(k,i_tmp)=rescore3_tmp
+             if (.not. out1file .or. me.eq.king)
+     &        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
@@ -2697,8 +2849,15 @@ c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
           do i = nnt,nct-2 
             do j=i+2,nct 
 
-            if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0) then
-                  
+              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) 
+c              write (iout,*) k,i,j,distal,dist2_cut
+
+            if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
+     &            .and. distal.le.dist2_cut ) then
+
               ii=ii+1
               ii_in_use(ii)=1
               l_homo(k,ii)=.true.
@@ -2708,13 +2867,6 @@ 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
               if (read2sigma) then
                 sigma_odl(k,ii)=0
      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
               else
                 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)
+                 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) 
                 else
 #ifdef OLDSIGMA
-                 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error 
+                 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* 
      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
 #else
-                 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error 
+                 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* 
      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 #endif
                 endif
               endif
-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)
-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
+              sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) 
             else
               ii=ii+1
               l_homo(k,ii)=.false.
@@ -2784,15 +2924,16 @@ 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)
-            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
+            sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
+     &                     rescore(k,i-2)+rescore(k,i-3))/4.0
+c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
 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))
+            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 
@@ -2821,17 +2962,17 @@ 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)
-             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))
+             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
+             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')
@@ -2857,15 +2998,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
@@ -2874,30 +3014,40 @@ c shift data in all arrays
 c
       if (waga_dist.ne.0.0d0) then
         ii=0
+        liiflag=.true.
         do i=nnt,nct-2 
          do j=i+2,nct 
           ii=ii+1
-          if (ii_in_use(ii).eq.0) then 
-             do ki=ii,lim_odl-1
-              ires_homo(ki)=ires_homo(ki+1)
-              jres_homo(ki)=jres_homo(ki+1)
-              ii_in_use(ki)=ii_in_use(ki+1)
+          if (ii_in_use(ii).eq.0.and.liiflag) then
+            liiflag=.false.
+            iistart=ii
+          endif
+          if (ii_in_use(ii).ne.0.and..not.liiflag.or.
+     &                   .not.liiflag.and.ii.eq.lim_odl) then
+             if (ii.eq.lim_odl) then
+              iishift=ii-iistart+1
+             else
+              iishift=ii-iistart
+             endif
+             liiflag=.true.
+             do ki=iistart,lim_odl-iishift
+              ires_homo(ki)=ires_homo(ki+iishift)
+              jres_homo(ki)=jres_homo(ki+iishift)
+              ii_in_use(ki)=ii_in_use(ki+iishift)
               do k=1,constr_homology
-               odl(k,ki)=odl(k,ki+1)
-               sigma_odl(k,ki)=sigma_odl(k,ki+1)
-               l_homo(k,ki)=l_homo(k,ki+1)
+               odl(k,ki)=odl(k,ki+iishift)
+               sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
+               l_homo(k,ki)=l_homo(k,ki+iishift)
               enddo
              enddo
-             ii=ii-1
-             lim_odl=lim_odl-1
+             ii=ii-iishift
+             lim_odl=lim_odl-iishift
           endif
          enddo
         enddo
       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
@@ -2908,24 +3058,26 @@ cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
        write (iout,*) "Distance restraints from templates"
        do ii=1,lim_odl
-       write(iout,'(3i5,10(2f8.2,1x,l1,4x))') 
+       write(iout,'(3i5,100(2f8.2,1x,l1,4x))') 
      &  ii,ires_homo(ii),jres_homo(ii),
      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(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),
+       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,10(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
-        write(iout,'(i5,10(4f8.2,4x))') i,
+       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)
        enddo