subroutine read_constr_homology
      implicit none
      include 'DIMENSIONS'
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.SETUP'
      include 'COMMON.CONTROL'
      include 'COMMON.CHAIN'
      include 'COMMON.IOUNITS'
      include 'COMMON.GEO'
      include 'COMMON.INTERACT'
      include 'COMMON.NAMES'
      include 'COMMON.HOMRESTR'
      include 'COMMON.HOMOLOGY'
c
c For new homol impl
c
      include 'COMMON.VAR'
c     include 'include_unres/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,ii,ik,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,
     & lim_theta,lim_xx,irec,iistart,iishift,i10,i01
      double precision distal
      integer idomain(max_template,maxres)
      logical lfirst
      integer ilen
      external ilen
      logical liiflag
      integer nres_temp
c
c     FP - Nov. 2014 Temporary specifications for new vars
c
      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
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
      call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
      call readi(controlcard,"HOMOL_NSET",homol_nset,1)
      read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
      call readi(controlcard,"IHSET",ihset,1)       
      write (iout,*) "homol_nset ",homol_nset
      if (homol_nset.gt.1)then
         call card_concat(controlcard)
         read(controlcard,*) (waga_homology(i),i=1,homol_nset) 
c         if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
c          write(iout,*) "iset homology_weight "
c         do i=1,homol_nset
c          write(iout,*) i,waga_homology(i)
c         enddo
c         endif
         iset=mod(kolor,homol_nset)+1
      else
       iset=1
       waga_homology(1)=1.0
      endif
c     write(iout,*) "waga_homology(",iset,")",waga_homology(iset)

cd      write (iout,*) "nnt",nnt," nct",nct
cd      call flush(iout)


      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
      do i = nnt,nct
        do k=1,constr_homology
         idomain(k,i)=0
        enddo
      enddo
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

      ii=0
      do i = nnt,nct-2 
        do j=i+2,nct 
        ii=ii+1
        ii_in_use(ii)=0
        enddo
      enddo
c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d

      if (read_homol_frag) then
        call read_klapaucjusz
      else

      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',
     &  pdbfile(:ilen(pdbfile))
        open(ipdbin,file=pdbfile,status='old',err=33)
        goto 34
  33    write (iout,'(a,5x,a)') 'Error opening PDB file',
     &  pdbfile(:ilen(pdbfile))
        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"

        unres_pdb=.false.
        nres_temp=nres
        if (read2sigma) then
          call readpdb_template(k)
          close(ipdbin)
        else
          call readpdb(out_template_coord)
          close(ipdbin)
        endif
        nres_chomo(k)=nres
        nres=nres_temp

        do i=1,2*nres
          do j=1,3
            crefjlee(j,i)=c(j,i)
          enddo
        enddo
#ifdef DEBUG
        do i=1,nres_chomo(k)
          write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
     &      (crefjlee(j,i+nres),j=1,3)
        enddo
        write (iout,*) "read_constr_homology: after reading pdb file"
        call flush(iout)
#endif

c
c     Distance restraints
c
c          ... --> odl(k,ii)
C Copy the coordinates from reference coordinates (?)
        do i=1,2*nres_chomo(k)
          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')
          if (nnt.gt.1) rescore(k,1)=0.0d0
          do irec=nnt,nct ! loop for reading res sim 
            if (read2sigma) then
             read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_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

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)
            endif
          enddo  
 1401   continue
        close (ientin)        
        if (waga_dist.ne.0.0d0) then
          ii=0
          do i = nnt,nct-2 
            do j=i+2,nct 

              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.

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
              odl(k,ii)=distal
              if (read2sigma) then
                sigma_odl(k,ii)=0
                do ik=i,j
                 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
                enddo
                sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
                if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = 
     &        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) 
                else
#ifdef OLDSIGMA
                 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))* 
     &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
#endif
                endif
              endif
              sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) 
            else
              ii=ii+1
              l_homo(k,ii)=.false.
            endif
            enddo
          enddo
        lim_odl=ii
        endif
c        write (iout,*) "Distance restraints set"
c        call flush(iout)
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 
            if (idomain(k,i).eq.0) then 
               sigma_dih(k,i)=0.0
               cycle
            endif
            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))/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
            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 
        endif
c        write (iout,*) "Dihedral angle restraints set"
c        call flush(iout)

        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.
             if (idomain(k,i).eq.0) then  
              sigma_theta(k,i)=0.0
              cycle
             endif
             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))/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
c        write (iout,*) "Angle restraints set"
c        call flush(iout)

        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

          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) cycle 
               if (idomain(k,i).eq.0 ) then 
                  sigma_d(k,i)=0.0
                  cycle
               endif
               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)=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
          enddo
        endif
      enddo
c      write (iout,*) "SC restraints set"
c      call flush(iout)
c
c remove distance restraints not used in any model from the list
c shift data in all arrays
c
c      write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
      if (waga_dist.ne.0.0d0) then
        ii=0
        liiflag=.true.
        lfirst=.true.
        do i=nnt,nct-2 
         do j=i+2,nct 
          ii=ii+1
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
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.
            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) then
             i01=ii
             liiflag=.true.
          endif
         enddo
        enddo
        lim_odl=iistart-1
      endif
c      write (iout,*) "Removing distances completed"
c      call flush(iout)
      endif ! .not. klapaucjusz

      if (constr_homology.gt.0) call homology_partition
c      write (iout,*) "After homology_partition"
c      call flush(iout)
      if (constr_homology.gt.0) call init_int_table
c      write (iout,*) "After init_int_table"
c      call flush(iout)
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.out_template_restr) return
cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
c      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))') 
     &  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)),
     &      (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)),
     &      (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,
     &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
     &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
       enddo
c      endif
c -----------------------------------------------------------------
      return
      end
c----------------------------------------------------------------------
      subroutine read_klapaucjusz

      include 'DIMENSIONS'
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.SETUP'
      include 'COMMON.CONTROL'
      include 'COMMON.HOMOLOGY'
      include 'COMMON.CHAIN'
      include 'COMMON.IOUNITS'
      include 'COMMON.GEO'
      include 'COMMON.INTERACT'
      include 'COMMON.NAMES'
      include 'COMMON.HOMRESTR'
      character*256 fragfile
      integer ninclust(maxclust),inclust(max_template,maxclust),
     &  nresclust(maxclust),iresclust(maxres,maxclust)

      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
      integer idomain(max_template,maxres)
      integer nres_temp
      logical lprn /.true./
      logical lfirst
      integer ilen
      external ilen
      logical liiflag
c
c
      double precision rescore_tmp,x12,y12,z12,rescore2_tmp
      double precision, dimension (max_template,maxres) :: rescore
      double precision, dimension (max_template,maxres) :: rescore2
      character*24 tpl_k_rescore

c
c For new homol impl
c
      include 'COMMON.VAR'
c
      call getenv("FRAGFILE",fragfile) 
      write (iout,*) "read_klapaucjusz ",fragfile
      open(ientin,file=fragfile,status="old",err=10)
      read(ientin,*) constr_homology,nclust
      l_homo = .false.
      sigma_theta=0.0
      sigma_d=0.0
      sigma_dih=0.0
c Read pdb files 
      do k=1,constr_homology 
        read(ientin,'(a)') pdbfile
        write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
     &  pdbfile(:ilen(pdbfile))
        open(ipdbin,file=pdbfile,status='old',err=33)
        goto 34
  33    write (iout,'(a,5x,a)') 'Error opening PDB file',
     &  pdbfile(:ilen(pdbfile))
        stop
  34    continue
        unres_pdb=.false.
        nres_temp=nres
        call readpdb_template(k)
        nres_chomo(k)=nres
        nres=nres_temp
c        do i=1,2*nres
c          do j=1,3
c            chomo(j,i,k)=c(j,i)
c          enddo
c        enddo
        do i=1,nres
          rescore(k,i)=0.2d0
          rescore2(k,i)=1.0d0
        enddo
      enddo
c Read clusters
      do i=1,nclust
        read(ientin,*) ninclust(i),nresclust(i)
        read(ientin,*) (inclust(k,i),k=1,ninclust(i))
        read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
      enddo
c
c Loop over clusters
c
      do l=1,nclust
        do ll = 1,ninclust(l)
        
        k = inclust(ll,l)
        do i=1,nres
          idomain(k,i)=0
        enddo
        do i=1,nresclust(l)
          if (nnt.gt.1)  then
            idomain(k,iresclust(i,l)+1) = 1
          else
            idomain(k,iresclust(i,l)) = 1
          endif
        enddo
c
c     Distance restraints
c
c          ... --> odl(k,ii)
C Copy the coordinates from reference coordinates (?)
        nres_temp=nres
        nres=nres_chomo(k)
        do i=1,2*nres
          do j=1,3
            c(j,i)=chomo(j,i,k)
c           write (iout,*) "c(",j,i,") =",c(j,i)
          enddo
        enddo
        call int_from_cart(.true.,.false.)
        call sc_loc_geom(.false.)
        do i=1,nres
          thetaref(i)=theta(i)
          phiref(i)=phi(i)
        enddo
        nres=nres_temp
        if (waga_dist.ne.0.0d0) then
          ii=0
          do i = nnt,nct-2 
            do j=i+2,nct 

              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.

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
              odl(k,ii)=distal
              if (read2sigma) then
                sigma_odl(k,ii)=0
                do ik=i,j
                 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
                enddo
                sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
                if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = 
     &        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) 
                else
#ifdef OLDSIGMA
                 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))* 
     &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
#endif
                endif
              endif
              sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) 
            else
              ii=ii+1
c              l_homo(k,ii)=.false.
            endif
            enddo
          enddo
        lim_odl=ii
        endif
c
c     Theta, dihedral and SC retraints
c
        if (waga_angle.gt.0.0d0) then
          do i = nnt+3,nct 
            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)+
     &                     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))
          enddo
          lim_dih=nct-nnt-2 
        endif

        if (waga_theta.gt.0.0d0) then
          do i = nnt+2,nct
             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)+
     &                        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))
          enddo
        endif

        if (waga_d.gt.0.0d0) then
          do i = nnt,nct
               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
          enddo
        endif
      enddo ! l
      enddo ! ll
c
c remove distance restraints not used in any model from the list
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.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+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
      endif
#ifdef DEBUG
      write (iout,*) "ires_homo and jres_homo arrays, lim_odl",lim_odl
      do i=1,lim_odl
        write (iout,*) i,ires_homo(i),jres_homo(i)
      enddo
#endif
      return
   10 stop "Error in fragment file"
      end