Adam's unres update
[unres.git] / source / unres / src-HCD-5D / readrtns_CSA.F
index 4fbc0f1..58b86cc 100644 (file)
@@ -149,6 +149,7 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword
       out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
       out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
       call readi(controlcard,'SYM',symetr,1)
+      call readi(controlcard,'PERMUT',npermut,1)
       call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
@@ -167,6 +168,8 @@ c       write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 c      endif
       call readi(controlcard,'NZ_START',nz_start,0)
       call readi(controlcard,'NZ_END',nz_end,0)
+      call readi(controlcard,'NRAN_START',nran_start,0)
+      write (iout,*) "nran_start",nran_start
 c      call readi(controlcard,'IZ_SC',iz_sc,0)
       timlim=60.0D0*timlim
       safety = 60.0d0*safety
@@ -481,6 +484,7 @@ C
       call readi(controlcard,"NSTEP",n_timestep,1000000)
       call readi(controlcard,"NTWE",ntwe,100)
       call readi(controlcard,"NTWX",ntwx,1000)
+      call readi(controlcard,"REST_FREQ",irest_freq,1000)
       call reada(controlcard,"DT",d_time,1.0d-1)
       call reada(controlcard,"DVMAX",dvmax,2.0d1)
       call reada(controlcard,"DAMAX",damax,1.0d1)
@@ -556,6 +560,7 @@ c  if performing umbrella sampling, fragments constrained are read from the frag
      &   " A"
        write(iout,'(a60,i5)')"Frequency of updating interaction list",
      &   imatupdate
+       write(iout,'(a60,i5)')"Restart writing frequency",irest_freq
        if (RESPA) then
         write (iout,'(2a,i4,a)') 
      &    "A-MTS algorithm used; initial time step for fast-varying",
@@ -741,7 +746,7 @@ C
       integer ilen
       external ilen
       integer iperm,tperm
-      integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp,itemp
+      integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp
       double precision sumv
 C
 C Read PDB structure if applicable
@@ -857,16 +862,19 @@ cd      print *,'NNT=',NNT,' NCT=',NCT
      &    chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
       enddo
       call chain_symmetry(nchain,nres,itype,chain_border,
-     &    chain_length,npermchain,tabpermchain)
+     &    chain_length,npermchain,tabpermchain,nchain_group,nequiv,
+     &    iequiv,mapchain)
 c      do i=1,nres
 c        write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
 c     &    ii=1,npermchain)
 c      enddo
+#ifdef DEBUG
       write(iout,*) "residue permutations"
       do i=1,nres
         write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
       enddo
       call flush(iout)
+#endif
       if (itype(1).eq.ntyp1) nnt=2
       if (itype(nres).eq.ntyp1) nct=nct-1
       write (iout,*) "nnt",nnt," nct",nct
@@ -971,25 +979,34 @@ C     &   write (iout,*) 'FTORS',ftors
         call reada(weightcard,"WDIHC",wdihc,0.591D0)
         write (iout,*) "Weight of dihedral angle restraints",wdihc
         read(inp,'(9x,3f7.3)') 
+c     &     (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
      &     (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
         write (iout,*) "The secprob array"
         do i=nnt,nct
           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
         enddo
         ndih_constr=0
+        iconstr_dih=0
         do i=nnt+3,nct
           if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1 
      &    .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
             ndih_constr=ndih_constr+1
             idih_constr(ndih_constr)=i
+            iconstr_dih(i)=ndih_constr
             sumv=0.0d0
             do j=1,3
               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
               sumv=sumv+vpsipred(j,ndih_constr)
             enddo 
-            do j=1,3
-              vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
-            enddo
+            if (sumv.gt.0.0d0) then
+              do j=1,3
+                vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
+              enddo
+            else
+              vpsipred(1,ndih_constr)=1.0d0
+              vpsipred(2,ndih_constr)=0.0d0
+              vpsipred(3,ndih_constr)=0.0d0
+            endif
             phibound(1,ndih_constr)=phihel*deg2rad
             phibound(2,ndih_constr)=phibet*deg2rad
             sdihed(1,ndih_constr)=sigmahel*deg2rad
@@ -1125,8 +1142,11 @@ c      write (iout,*) "After read_dist_constr nhpb",nhpb
       endif
       write (iout,*) "calling read_saxs_consrtr",nsaxs
       if (nsaxs.gt.0) call read_saxs_constr
-
+c      write (iout,*) "After read_saxs_constr"
+c      call flush(iout) 
       if (constr_homology.gt.0) then
+c        write (iout,*) "Calling read_constr_homology"
+c        call flush(iout)
         call read_constr_homology
         if (indpdb.gt.0 .or. pdbref) then
           do i=1,2*nres
@@ -1184,7 +1204,8 @@ c      write (iout,*) "After read_dist_constr nhpb",nhpb
  335        continue
             unres_pdb=.false.
             nres_temp=nres
-            call readpdb
+c            call readpdb
+            call readpdb_template(nmodel_start+1)
             close(ipdbin)
             if (nres.ge.nres_temp) then
               nmodel_start=nmodel_start+1
@@ -1195,22 +1216,11 @@ c      write (iout,*) "After read_dist_constr nhpb",nhpb
                 enddo
               enddo
             else
-c              itemp=nres
-c              nres=nres_temp
-c              call gen_rand_conf(itemp,*115)
-c              nmodel_start=nmodel_start+1
-c              do i=1,2*nres
-c                do j=1,3
-c                  chomo(j,i,nmodel_start)=c(j,i)
-c                enddo
-c              enddo
-c              goto 116
-  115         if (me.eq.king .or. .not. out1file)
-     &          write (iout,'(a,2i5,1x,a)') 
+              if (me.eq.king .or. .not. out1file)
+     &          write (iout,'(a,2i7,1x,a)') 
      &           "Different number of residues",nres_temp,nres,
      &           " model skipped."
             endif
-  116       continue
             nres=nres_temp
           enddo
   332     continue
@@ -2341,6 +2351,7 @@ c1out       open(iout,file=outname,status='unknown')
 #else
       if (me.eq.king .or. .not.out1file)
      &    open(iout,file=outname,status='unknown')
+#define DEBUG
 #ifdef DEBUG
       if (fg_rank.gt.0) then
         write (liczba,'(i3.3)') myrank/nfgtasks
@@ -2349,6 +2360,7 @@ c1out       open(iout,file=outname,status='unknown')
      &   status='unknown')
       endif
 #endif
+#undef DEBUG
       if(me.eq.king) then
        open(igeom,file=intname,status='unknown',access='append')
        open(ipdb,file=pdbname,status='unknown')
@@ -2734,10 +2746,10 @@ c            write (iout,*) "j",j," k",k
             endif
 #ifdef MPI
             if (.not.out1file .or. me.eq.king) 
-     &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
+     &      write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 #else
-            write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
+            write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 #endif
           enddo
@@ -2781,10 +2793,10 @@ c            write (iout,*) "j",j," k",k
             endif
 #ifdef MPI
             if (.not.out1file .or. me.eq.king)
-     &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
+     &      write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 #else
-            write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
+            write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 #endif
           enddo
@@ -2811,13 +2823,13 @@ c     &    fordepth_peak(nhpb_peak+1),npeak
           ipeak(2,npeak)=i
 #ifdef MPI
           if (.not.out1file .or. me.eq.king)
-     &    write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+     &    write (iout,'(a,5i6,2f8.2,2f10.5,i5)') "+dist.restr ",
      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
 #else
-          write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+          write (iout,'(a,5i6,2f8.2,2f10.5,i5)') "+dist.restr ",
      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
@@ -2840,11 +2852,11 @@ c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
           irestr_type(nhpb)=11
 #ifdef MPI
           if (.not.out1file .or. me.eq.king)
-     &    write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+     &    write (iout,'(a,4i6,2f8.2,2f10.5,i5)') "+dist.restr ",
      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
 #else
-          write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+          write (iout,'(a,4i6,2f8.2,2f10.5,i5)') "+dist.restr ",
      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
 #endif
@@ -2904,12 +2916,12 @@ c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
           endif
 #ifdef MPI
           if (.not.out1file .or. me.eq.king)
-     &    write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+     &    write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
      &     irestr_type(nhpb)
 #else
-          write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+          write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
      &     irestr_type(nhpb)
@@ -2938,10 +2950,10 @@ C        print *,"in else"
         endif
 #ifdef MPI
           if (.not.out1file .or. me.eq.king)
-     &    write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
+     &    write (iout,'(a,4i6,f8.2,f10.1)') "+dist.restr ",
      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
 #else
-          write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
+          write (iout,'(a,4i6,f8.2,f10.1)') "+dist.restr ",
      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
 #endif
         endif
@@ -2984,13 +2996,13 @@ C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
             dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
 #ifdef MPI
             if (.not.out1file .or. me.eq.king) then
-            write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+            write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
      &       irestr_type(nhpb)
             endif
 #else
-            write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+            write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
      &       irestr_type(nhpb)
@@ -3054,8 +3066,8 @@ 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,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
-     & ik,iistart,nres_temp
+      integer ki,i,ii,j,k,l,ii_in_use(maxdim_cont),i_tmp,idomain_tmp,
+     & irec,ik,iistart,nres_temp
       integer ilen
       external ilen
       logical liiflag,lfirst
@@ -3108,7 +3120,7 @@ c      endif
          call card_concat(controlcard)
          read(controlcard,*) (waga_homology(i),i=1,homol_nset) 
          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
-          write(iout,*) "iset homology_weight "
+c          write(iout,*) "iset homology_weight "
           do i=1,homol_nset
            write(iout,*) i,waga_homology(i)
           enddo
@@ -3122,30 +3134,30 @@ c      endif
 cd      write (iout,*) "nnt",nnt," nct",nct
 cd      call flush(iout)
 
+      if (read_homol_frag) then
+        call read_klapaucjusz
+      else
 
       lim_odl=0
       lim_dih=0
 c
 c      write(iout,*) 'nnt=',nnt,'nct=',nct
 c
-      do i = nnt,nct
-        do k=1,constr_homology
-         idomain(k,i)=0
-        enddo
-      enddo
-
-      ii=0
-      do i = nnt,nct-2 
-        do j=i+2,nct 
-        ii=ii+1
-        ii_in_use(ii)=0
-        enddo
-      enddo
-
-      if (read_homol_frag) then
-        call read_klapaucjusz
-      else
+c      do i = nnt,nct
+c        do k=1,constr_homology
+c         idomain(k,i)=0
+c        enddo
+c      enddo
+       idomain=0
 
+c      ii=0
+c      do i = nnt,nct-2 
+c        do j=i+2,nct 
+c        ii=ii+1
+c        ii_in_use(ii)=0
+c        enddo
+c      enddo
+      ii_in_use=0
       do k=1,constr_homology
 
         read(inp,'(a)') pdbfile
@@ -3265,8 +3277,8 @@ c    &                       constr_homology
               endif
               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) 
             else
-              ii=ii+1
-              l_homo(k,ii)=.false.
+c              ii=ii+1
+c              l_homo(k,ii)=.false.
             endif
             enddo
           enddo
@@ -3456,26 +3468,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,100(2f8.2,1x,l1,4x))') 
+       write(iout,'(3i7,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,nct
-        write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
+        write (iout,'(i7,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,nct
-        write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
+        write (iout,'(i7,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,nct
-        write(iout,'(i5,100(4f8.2,4x))') i,
+        write(iout,'(i7,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
@@ -3641,14 +3653,17 @@ c----------------------------------------------------------------------
       character*2 kic2
       character*24 model_ki_dist, model_ki_angle
       character*500 controlcard
-      integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
-     & ik,ll,ii,kk,iistart,iishift,lim_xx
+      integer ki, i, j, jj,k, l, ii_in_use(maxdim_cont),i_tmp,
+     & idomain_tmp,
+     & ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,
+     & i01,i10,nnt_chain,nct_chain
+      integer itype_temp(maxres)
       double precision distal
       logical lprn /.true./
       integer nres_temp
       integer ilen
       external ilen
-      logical liiflag
+      logical liiflag,lfirst
 c
 c
       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
@@ -3662,14 +3677,48 @@ c For new homol impl
 c
       include 'COMMON.VAR'
 c
+c      write (iout,*) "READ_KLAPAUCJUSZ"
+c      print *,"READ_KLAPAUCJUSZ"
+c      call flush(iout)
       call getenv("FRAGFILE",fragfile) 
+      write (iout,*) "Opening", fragfile
+      call flush(iout)
       open(ientin,file=fragfile,status="old",err=10)
-      read(ientin,*) constr_homology,nclust
-      nmodel_start=constr_homology
-      l_homo = .false.
+c      write (iout,*) " opened"
+c      call flush(iout)
+
       sigma_theta=0.0
       sigma_d=0.0
       sigma_dih=0.0
+      l_homo = .false.
+
+      nres_temp=nres
+      itype_temp=itype
+      ii=0
+      lim_odl=0
+
+c      write (iout,*) "Entering loop"
+c      call flush(iout)
+
+      DO IGR = 1,NCHAIN_GROUP
+
+c      write (iout,*) "igr",igr
+      call flush(iout)
+      read(ientin,*) constr_homology,nclust
+
+      if (start_from_model) then
+        nmodel_start=constr_homology
+      else
+        nmodel_start=0
+      endif
+
+      ii_old=lim_odl
+
+      ichain=iequiv(1,igr)
+      nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1
+      nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1
+c      write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain
+
 c Read pdb files 
       do k=1,constr_homology 
         read(ientin,'(a)') pdbfile
@@ -3683,10 +3732,10 @@ c Read pdb files
         stop
   34    continue
         unres_pdb=.false.
-        nres_temp=nres
+c        nres_temp=nres
         call readpdb_template(k)
         nres_chomo(k)=nres
-        nres=nres_temp
+c        nres=nres_temp
         do i=1,nres
           rescore(k,i)=0.2d0
           rescore2(k,i)=1.0d0
@@ -3705,6 +3754,7 @@ c
         do ll = 1,ninclust(l)
         
         k = inclust(ll,l)
+c        write (iout,*) "l",l," ll",ll," k",k
         do i=1,nres
           idomain(k,i)=0
         enddo
@@ -3720,7 +3770,7 @@ c     Distance restraints
 c
 c          ... --> odl(k,ii)
 C Copy the coordinates from reference coordinates (?)
-        nres_temp=nres
+c        nres_temp=nres
         nres=nres_chomo(k)
         do i=1,2*nres
           do j=1,3
@@ -3734,11 +3784,13 @@ c           write (iout,*) "c(",j,i,") =",c(j,i)
           thetaref(i)=theta(i)
           phiref(i)=phi(i)
         enddo
-        nres=nres_temp
+c        nres=nres_temp
         if (waga_dist.ne.0.0d0) then
-          ii=0
-          do i = nnt,nct-2 
-            do j=i+2,nct 
+          ii=ii_old
+c          do i = nnt,nct-2 
+          do i = nnt_chain,nct_chain-2 
+c            do j=i+2,nct 
+            do j=i+2,nct_chain
 
               x12=c(1,i)-c(1,j)
               y12=c(2,i)-c(2,j)
@@ -3755,9 +3807,9 @@ c              write (iout,*) k,i,j,distal,dist2_cut
 
 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     &                       constr_homology
+              ires_homo(ii)=i+chain_border1(1,igr)-1
+              jres_homo(ii)=j+chain_border1(1,igr)-1
               odl(k,ii)=distal
               if (read2sigma) then
                 sigma_odl(k,ii)=0
 c     Theta, dihedral and SC retraints
 c
         if (waga_angle.gt.0.0d0) then
-          do i = nnt+3,nct 
+          do i = nnt_chain+3,nct_chain
+            iii=i+chain_border1(1,igr)-1
             if (idomain(k,i).eq.0) then 
 c               sigma_dih(k,i)=0.0
                cycle
             endif
-            dih(k,i)=phiref(i)
-            sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
+            dih(k,iii)=phiref(i)
+            sigma_dih(k,iii)=
+     &         (rescore(k,i)+rescore(k,i-1)+
      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
 c     &       " sigma_dihed",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))
+            if (sigma_dih(k,iii).ne.0)
+     &       sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii))
           enddo
-          lim_dih=nct-nnt-2 
+c          lim_dih=nct-nnt-2 
         endif
 
         if (waga_theta.gt.0.0d0) then
-          do i = nnt+2,nct
+          do i = nnt_chain+2,nct_chain
+             iii=i+chain_border1(1,igr)-1
              if (idomain(k,i).eq.0) then  
 c              sigma_theta(k,i)=0.0
               cycle
              endif
-             thetatpl(k,i)=thetaref(i)
-             sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
+             thetatpl(k,iii)=thetaref(i)
+             sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+
      &                        rescore(k,i-2))/3.0
-             if (sigma_theta(k,i).ne.0)
-     &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
+             if (sigma_theta(k,iii).ne.0)
+     &       sigma_theta(k,iii)=1.0d0/
+     &       (sigma_theta(k,iii)*sigma_theta(k,iii))
           enddo
         endif
 
         if (waga_d.gt.0.0d0) then
-          do i = nnt,nct
+          do i = nnt_chain,nct_chain
+             iii=i+chain_border1(1,igr)-1
                if (itype(i).eq.10) cycle 
                if (idomain(k,i).eq.0 ) then 
 c                  sigma_d(k,i)=0.0
                   cycle
                endif
-               xxtpl(k,i)=xxref(i)
-               yytpl(k,i)=yyref(i)
-               zztpl(k,i)=zzref(i)
-               sigma_d(k,i)=rescore(k,i)
-               if (sigma_d(k,i).ne.0)
-     &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
-               if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
+               xxtpl(k,iii)=xxref(i)
+               yytpl(k,iii)=yyref(i)
+               zztpl(k,iii)=zzref(i)
+               sigma_d(k,iii)=rescore(k,i)
+               if (sigma_d(k,iii).ne.0)
+     &          sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii))
+c               if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
           enddo
         endif
       enddo ! l
 c remove distance restraints not used in any model from the list
 c shift data in all arrays
 c
+c      write (iout,*) "ii_old",ii_old
       if (waga_dist.ne.0.0d0) then
-        ii=0
+#ifdef DEBUG
+       write (iout,*) "Distance restraints from templates"
+       do iii=1,lim_odl
+       write(iout,'(4i5,100(2f8.2,1x,l1,4x))') 
+     &  iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii),
+     &  (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii),
+     &  ki=1,constr_homology)
+       enddo
+#endif
+        ii=ii_old
         liiflag=.true.
-        do i=nnt,nct-2 
-         do j=i+2,nct 
+        lfirst=.true.
+        do i=nnt_chain,nct_chain-2
+         do j=i+2,nct_chain
           ii=ii+1
-          if (ii_in_use(ii).eq.0.and.liiflag) then
+c          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
+c     &            .and. distal.le.dist2_cut ) then
+c          write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii)
+c          call flush(iout)
+          if (ii_in_use(ii).eq.0.and.liiflag.or.
+     &     ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
             liiflag=.false.
-            iistart=ii
+            i10=ii
+            if (lfirst) then
+              lfirst=.false.
+              iistart=ii
+            else
+              if(i10.eq.lim_odl) i10=i10+1
+              do ki=0,i10-i01-1
+               ires_homo(iistart+ki)=ires_homo(ki+i01)
+               jres_homo(iistart+ki)=jres_homo(ki+i01)
+               ii_in_use(iistart+ki)=ii_in_use(ki+i01)
+               do k=1,constr_homology
+                odl(k,iistart+ki)=odl(k,ki+i01)
+                sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
+                l_homo(k,iistart+ki)=l_homo(k,ki+i01)
+               enddo
+              enddo
+              iistart=iistart+i10-i01
+            endif
           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
+          if (ii_in_use(ii).ne.0.and..not.liiflag) then
+             i01=ii
              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+iishift)
-               sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
-               l_homo(k,ki)=l_homo(k,ki+iishift)
-              enddo
-             enddo
-             ii=ii-iishift
-             lim_odl=lim_odl-iishift
           endif
          enddo
         enddo
+        lim_odl=iistart-1
       endif
 
+      lll=lim_odl-ii_old
+
+      do i=2,nequiv(igr) 
+
+        ichain=iequiv(i,igr)
+
+        do j=nnt_chain,nct_chain
+          jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr))
+          do k=1,constr_homology
+            dih(k,jj)=dih(k,j)
+            sigma_dih(k,jj)=sigma_dih(k,j)
+            thetatpl(k,jj)=thetatpl(k,j)
+            sigma_theta(k,jj)=sigma_theta(k,j)
+            xxtpl(k,jj)=xxtpl(k,j)
+            yytpl(k,jj)=yytpl(k,j)
+            zztpl(k,jj)=zztpl(k,j)
+            sigma_d(k,jj)=sigma_d(k,j)
+          enddo
+        enddo
+
+        jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr)) 
+c        write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj
+        do j=ii_old+1,lim_odl
+          ires_homo(j+lll)=ires_homo(j)+jj
+          jres_homo(j+lll)=jres_homo(j)+jj
+          do k=1,constr_homology
+            odl(k,j+lll)=odl(k,j)
+            sigma_odl(k,j+lll)=sigma_odl(k,j)
+            l_homo(k,j+lll)=l_homo(k,j)
+          enddo
+        enddo
+
+        ii_old=ii_old+lll
+        lim_odl=lim_odl+lll
+
+      enddo
+
+      ENDDO ! IGR
+
+      if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2 
+      nres=nres_temp
+      itype=itype_temp
+
       return
    10 stop "Error in fragment file"
       end