source/unres/src-HCD-5D/read_constr_homology.F

   1       subroutine read_constr_homology
   2       implicit none
   3       include 'DIMENSIONS'
   4 #ifdef MPI
   5       include 'mpif.h'
   6 #endif
   7       include 'COMMON.SETUP'
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.HOMOLOGY'
  10       include 'COMMON.CHAIN'
  11       include 'COMMON.IOUNITS'
  12       include 'COMMON.MD'
  13       include 'COMMON.QRESTR'
  14       include 'COMMON.GEO'
  15       include 'COMMON.INTERACT'
  16       include 'COMMON.NAMES'
  17 c
  18 c For new homol impl
  19 c
  20       include 'COMMON.VAR'
  21 c
  22
  23 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
  24 c    &                 dist_cut
  25 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
  26 c    &    sigma_odl_temp(maxres,maxres,max_template)
  27       character*2 kic2
  28       character*24 model_ki_dist, model_ki_angle
  29       character*500 controlcard
  30       integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
  31      & ik,iistart,nres_temp
  32       integer ilen
  33       external ilen
  34       logical liiflag,lfirst
  35       integer i01,i10
  36 c
  37 c     FP - Nov. 2014 Temporary specifications for new vars
  38 c
  39       double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
  40      &                 rescore3_tmp
  41       double precision, dimension (max_template,maxres) :: rescore
  42       double precision, dimension (max_template,maxres) :: rescore2
  43       double precision, dimension (max_template,maxres) :: rescore3
  44       double precision distal
  45       character*24 pdbfile,tpl_k_rescore
  46 c -----------------------------------------------------------------
  47 c Reading multiple PDB ref structures and calculation of retraints
  48 c not using pre-computed ones stored in files model_ki_{dist,angle}
  49 c FP (Nov., 2014)
  50 c -----------------------------------------------------------------
  51 c
  52 c
  53 c Alternative: reading from input
  54       call card_concat(controlcard)
  55       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
  56       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
  57       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
  58       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
  59       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
  60       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
  61       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
  62       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
  63       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
  64       if(.not.read2sigma.and.start_from_model) then
  65           if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
  66      &      write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
  67           start_from_model=.false.
  68       endif
  69       if(start_from_model .and. (me.eq.king .or. .not. out1file))
  70      &    write(iout,*) 'START_FROM_MODELS is ON'
  71       if(start_from_model .and. rest) then
  72         if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
  73          write(iout,*) 'START_FROM_MODELS is OFF'
  74          write(iout,*) 'remove restart keyword from input'
  75         endif
  76       endif
  77       if (homol_nset.gt.1)then
  78          call card_concat(controlcard)
  79          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
  80          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
  81           write(iout,*) "iset homology_weight "
  82           do i=1,homol_nset
  83            write(iout,*) i,waga_homology(i)
  84           enddo
  85          endif
  86          iset=mod(kolor,homol_nset)+1
  87       else
  88        iset=1
  89        waga_homology(1)=1.0
  90       endif
  91
  92 cd      write (iout,*) "nnt",nnt," nct",nct
  93 cd      call flush(iout)
  94
  95
  96       lim_odl=0
  97       lim_dih=0
  98 c
  99 c      write(iout,*) 'nnt=',nnt,'nct=',nct
 100 c
 101       do i = nnt,nct
 102         do k=1,constr_homology
 103          idomain(k,i)=0
 104         enddo
 105       enddo
 106
 107       ii=0
 108       do i = nnt,nct-2
 109         do j=i+2,nct
 110         ii=ii+1
 111         ii_in_use(ii)=0
 112         enddo
 113       enddo
 114
 115       if (read_homol_frag) then
 116         call read_klapaucjusz
 117       else
 118
 119       do k=1,constr_homology
 120
 121         read(inp,'(a)') pdbfile
 122         if(me.eq.king .or. .not. out1file)
 123      &    write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
 124      &  pdbfile(:ilen(pdbfile))
 125         open(ipdbin,file=pdbfile,status='old',err=33)
 126         goto 34
 127   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
 128      &  pdbfile(:ilen(pdbfile))
 129         stop
 130   34    continue
 131 c        print *,'Begin reading pdb data'
 132 c
 133 c Files containing res sim or local scores (former containing sigmas)
 134 c
 135
 136         write(kic2,'(bz,i2.2)') k
 137
 138         tpl_k_rescore="template"//kic2//".sco"
 139
 140         unres_pdb=.false.
 141         nres_temp=nres
 142         if (read2sigma) then
 143           call readpdb_template(k)
 144         else
 145           call readpdb
 146         endif
 147         nres_chomo(k)=nres
 148         nres=nres_temp
 149 c
 150 c     Distance restraints
 151 c
 152 c          ... --> odl(k,ii)
 153 C Copy the coordinates from reference coordinates (?)
 154         do i=1,2*nres_chomo(k)
 155           do j=1,3
 156             c(j,i)=cref(j,i)
 157 c           write (iout,*) "c(",j,i,") =",c(j,i)
 158           enddo
 159         enddo
 160 c
 161 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
 162 c
 163 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
 164           open (ientin,file=tpl_k_rescore,status='old')
 165           if (nnt.gt.1) rescore(k,1)=0.0d0
 166           do irec=nnt,nct ! loop for reading res sim
 167             if (read2sigma) then
 168              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
 169      &                                rescore3_tmp,idomain_tmp
 170              i_tmp=i_tmp+nnt-1
 171              idomain(k,i_tmp)=idomain_tmp
 172              rescore(k,i_tmp)=rescore_tmp
 173              rescore2(k,i_tmp)=rescore2_tmp
 174              rescore3(k,i_tmp)=rescore3_tmp
 175              if (.not. out1file .or. me.eq.king)
 176      &        write(iout,'(a7,i5,3f10.5,i5)') "rescore",
 177      &                      i_tmp,rescore2_tmp,rescore_tmp,
 178      &                                rescore3_tmp,idomain_tmp
 179             else
 180              idomain(k,irec)=1
 181              read (ientin,*,end=1401) rescore_tmp
 182
 183 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
 184              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
 185 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
 186             endif
 187           enddo
 188  1401   continue
 189         close (ientin)
 190         if (waga_dist.ne.0.0d0) then
 191           ii=0
 192           do i = nnt,nct-2
 193             do j=i+2,nct
 194
 195               x12=c(1,i)-c(1,j)
 196               y12=c(2,i)-c(2,j)
 197               z12=c(3,i)-c(3,j)
 198               distal=dsqrt(x12*x12+y12*y12+z12*z12)
 199 c              write (iout,*) k,i,j,distal,dist2_cut
 200
 201             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
 202      &            .and. distal.le.dist2_cut ) then
 203
 204               ii=ii+1
 205               ii_in_use(ii)=1
 206               l_homo(k,ii)=.true.
 207
 208 c             write (iout,*) "k",k
 209 c             write (iout,*) "i",i," j",j," constr_homology",
 210 c    &                       constr_homology
 211               ires_homo(ii)=i
 212               jres_homo(ii)=j
 213               odl(k,ii)=distal
 214               if (read2sigma) then
 215                 sigma_odl(k,ii)=0
 216                 do ik=i,j
 217                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
 218                 enddo
 219                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
 220                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
 221      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 222               else
 223                 if (odl(k,ii).le.dist_cut) then
 224                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
 225                 else
 226 #ifdef OLDSIGMA
 227                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 228      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
 229 #else
 230                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 231      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 232 #endif
 233                 endif
 234               endif
 235               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
 236             else
 237               ii=ii+1
 238               l_homo(k,ii)=.false.
 239             endif
 240             enddo
 241           enddo
 242         lim_odl=ii
 243         endif
 244 c        write (iout,*) "Distance restraints set"
 245 c        call flush(iout)
 246 c
 247 c     Theta, dihedral and SC retraints
 248 c
 249         if (waga_angle.gt.0.0d0) then
 250 c         open (ientin,file=tpl_k_sigma_dih,status='old')
 251 c         do irec=1,maxres-3 ! loop for reading sigma_dih
 252 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
 253 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
 254 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
 255 c    &                            sigma_dih(k,i+nnt-1)
 256 c         enddo
 257 c1402   continue
 258 c         close (ientin)
 259           do i = nnt+3,nct
 260             if (idomain(k,i).eq.0) then
 261                sigma_dih(k,i)=0.0
 262                cycle
 263             endif
 264             dih(k,i)=phiref(i) ! right?
 265 c           read (ientin,*) sigma_dih(k,i) ! original variant
 266 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
 267 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
 268 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 269 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
 270 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
 271
 272             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
 273      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
 274 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
 275 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
 276 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
 277 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 278 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
 279 c   Instead of res sim other local measure of b/b str reliability possible
 280             if (sigma_dih(k,i).ne.0)
 281      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
 282 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
 283           enddo
 284           lim_dih=nct-nnt-2
 285         endif
 286 c        write (iout,*) "Dihedral angle restraints set"
 287 c        call flush(iout)
 288
 289         if (waga_theta.gt.0.0d0) then
 290 c         open (ientin,file=tpl_k_sigma_theta,status='old')
 291 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
 292 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
 293 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
 294 c    &                              sigma_theta(k,i+nnt-1)
 295 c         enddo
 296 c1403   continue
 297 c         close (ientin)
 298
 299           do i = nnt+2,nct ! right? without parallel.
 300 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
 301 c         do i=ithet_start,ithet_end ! with FG parallel.
 302              if (idomain(k,i).eq.0) then
 303               sigma_theta(k,i)=0.0
 304               cycle
 305              endif
 306              thetatpl(k,i)=thetaref(i)
 307 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
 308 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
 309 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 310 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
 311 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
 312              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
 313      &                        rescore(k,i-2))/3.0
 314 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
 315              if (sigma_theta(k,i).ne.0)
 316      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 317
 318 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 319 c                             rescore(k,i-2) !  right expression ?
 320 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
 321           enddo
 322         endif
 323 c        write (iout,*) "Angle restraints set"
 324 c        call flush(iout)
 325
 326         if (waga_d.gt.0.0d0) then
 327 c       open (ientin,file=tpl_k_sigma_d,status='old')
 328 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
 329 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
 330 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
 331 c    &                          sigma_d(k,i+nnt-1)
 332 c         enddo
 333 c1404   continue
 334
 335           do i = nnt,nct ! right? without parallel.
 336 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
 337 c         do i=loc_start,loc_end ! with FG parallel.
 338                if (itype(i).eq.10) cycle
 339                if (idomain(k,i).eq.0 ) then
 340                   sigma_d(k,i)=0.0
 341                   cycle
 342                endif
 343                xxtpl(k,i)=xxref(i)
 344                yytpl(k,i)=yyref(i)
 345                zztpl(k,i)=zzref(i)
 346 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
 347 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
 348 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
 349 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
 350                sigma_d(k,i)=rescore3(k,i) !  right expression ?
 351                if (sigma_d(k,i).ne.0)
 352      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
 353
 354 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
 355 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
 356 c              read (ientin,*) sigma_d(k,i) ! 1st variant
 357           enddo
 358         endif
 359       enddo
 360 c      write (iout,*) "SC restraints set"
 361 c      call flush(iout)
 362 c
 363 c remove distance restraints not used in any model from the list
 364 c shift data in all arrays
 365 c
 366 c      write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
 367       if (waga_dist.ne.0.0d0) then
 368         ii=0
 369         liiflag=.true.
 370         lfirst=.true.
 371         do i=nnt,nct-2
 372          do j=i+2,nct
 373           ii=ii+1
 374 c          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
 375 c     &            .and. distal.le.dist2_cut ) then
 376 c          write (iout,*) "i",i," j",j," ii",ii
 377 c          call flush(iout)
 378           if (ii_in_use(ii).eq.0.and.liiflag.or.
 379      &     ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
 380             liiflag=.false.
 381             i10=ii
 382             if (lfirst) then
 383               lfirst=.false.
 384               iistart=ii
 385             else
 386               if(i10.eq.lim_odl) i10=i10+1
 387               do ki=0,i10-i01-1
 388                ires_homo(iistart+ki)=ires_homo(ki+i01)
 389                jres_homo(iistart+ki)=jres_homo(ki+i01)
 390                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
 391                do k=1,constr_homology
 392                 odl(k,iistart+ki)=odl(k,ki+i01)
 393                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
 394                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
 395                enddo
 396               enddo
 397               iistart=iistart+i10-i01
 398             endif
 399           endif
 400           if (ii_in_use(ii).ne.0.and..not.liiflag) then
 401              i01=ii
 402              liiflag=.true.
 403           endif
 404          enddo
 405         enddo
 406         lim_odl=iistart-1
 407       endif
 408 c      write (iout,*) "Removing distances completed"
 409 c      call flush(iout)
 410       endif ! .not. klapaucjusz
 411
 412       if (constr_homology.gt.0) call homology_partition
 413 c      write (iout,*) "After homology_partition"
 414 c      call flush(iout)
 415       if (constr_homology.gt.0) call init_int_table
 416 c      write (iout,*) "After init_int_table"
 417 c      call flush(iout)
 418 c      write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
 419 c      write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
 420 c
 421 c Print restraints
 422 c
 423       if (.not.out_template_restr) return
 424 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 425       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 426        write (iout,*) "Distance restraints from templates"
 427        do ii=1,lim_odl
 428        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
 429      &  ii,ires_homo(ii),jres_homo(ii),
 430      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
 431      &  ki=1,constr_homology)
 432        enddo
 433        write (iout,*) "Dihedral angle restraints from templates"
 434        do i=nnt+3,nct
 435         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
 436      &      (rad2deg*dih(ki,i),
 437      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
 438        enddo
 439        write (iout,*) "Virtual-bond angle restraints from templates"
 440        do i=nnt+2,nct
 441         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
 442      &      (rad2deg*thetatpl(ki,i),
 443      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
 444        enddo
 445        write (iout,*) "SC restraints from templates"
 446        do i=nnt,nct
 447         write(iout,'(i5,100(4f8.2,4x))') i,
 448      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
 449      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
 450        enddo
 451       endif
 452 c -----------------------------------------------------------------
 453       return
 454       end
 455 c----------------------------------------------------------------------
 456       subroutine read_klapaucjusz
 457       implicit none
 458       include 'DIMENSIONS'
 459 #ifdef MPI
 460       include 'mpif.h'
 461 #endif
 462       include 'COMMON.SETUP'
 463       include 'COMMON.CONTROL'
 464       include 'COMMON.HOMOLOGY'
 465       include 'COMMON.CHAIN'
 466       include 'COMMON.IOUNITS'
 467       include 'COMMON.MD'
 468       include 'COMMON.GEO'
 469       include 'COMMON.INTERACT'
 470       include 'COMMON.NAMES'
 471       character*256 fragfile
 472       integer ninclust(maxclust),inclust(max_template,maxclust),
 473      &  nresclust(maxclust),iresclust(maxres,maxclust),nclust
 474
 475       character*2 kic2
 476       character*24 model_ki_dist, model_ki_angle
 477       character*500 controlcard
 478       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
 479      & ik,ll,ii,kk,iistart,iishift,lim_xx
 480       double precision distal
 481       logical lprn /.true./
 482       integer nres_temp
 483       integer ilen
 484       external ilen
 485       logical liiflag
 486 c
 487 c
 488       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
 489       double precision, dimension (max_template,maxres) :: rescore
 490       double precision, dimension (max_template,maxres) :: rescore2
 491       character*24 pdbfile,tpl_k_rescore
 492
 493 c
 494 c For new homol impl
 495 c
 496       include 'COMMON.VAR'
 497 c
 498       call getenv("FRAGFILE",fragfile)
 499       open(ientin,file=fragfile,status="old",err=10)
 500       read(ientin,*) constr_homology,nclust
 501       l_homo = .false.
 502       sigma_theta=0.0
 503       sigma_d=0.0
 504       sigma_dih=0.0
 505 c Read pdb files
 506       do k=1,constr_homology
 507         read(ientin,'(a)') pdbfile
 508         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
 509      &  pdbfile(:ilen(pdbfile))
 510         open(ipdbin,file=pdbfile,status='old',err=33)
 511         goto 34
 512   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
 513      &  pdbfile(:ilen(pdbfile))
 514         stop
 515   34    continue
 516         unres_pdb=.false.
 517         nres_temp=nres
 518         call readpdb_template(k)
 519         nres_chomo(k)=nres
 520         nres=nres_temp
 521         do i=1,nres
 522           rescore(k,i)=0.2d0
 523           rescore2(k,i)=1.0d0
 524         enddo
 525       enddo
 526 c Read clusters
 527       do i=1,nclust
 528         read(ientin,*) ninclust(i),nresclust(i)
 529         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
 530         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
 531       enddo
 532 c
 533 c Loop over clusters
 534 c
 535       do l=1,nclust
 536         do ll = 1,ninclust(l)
 537
 538         k = inclust(ll,l)
 539         do i=1,nres
 540           idomain(k,i)=0
 541         enddo
 542         do i=1,nresclust(l)
 543           if (nnt.gt.1)  then
 544             idomain(k,iresclust(i,l)+1) = 1
 545           else
 546             idomain(k,iresclust(i,l)) = 1
 547           endif
 548         enddo
 549 c
 550 c     Distance restraints
 551 c
 552 c          ... --> odl(k,ii)
 553 C Copy the coordinates from reference coordinates (?)
 554         nres_temp=nres
 555         nres=nres_chomo(k)
 556         do i=1,2*nres
 557           do j=1,3
 558             c(j,i)=chomo(j,i,k)
 559 c           write (iout,*) "c(",j,i,") =",c(j,i)
 560           enddo
 561         enddo
 562         call int_from_cart(.true.,.false.)
 563         call sc_loc_geom(.false.)
 564         do i=1,nres
 565           thetaref(i)=theta(i)
 566           phiref(i)=phi(i)
 567         enddo
 568         nres=nres_temp
 569         if (waga_dist.ne.0.0d0) then
 570           ii=0
 571           do i = nnt,nct-2
 572             do j=i+2,nct
 573
 574               x12=c(1,i)-c(1,j)
 575               y12=c(2,i)-c(2,j)
 576               z12=c(3,i)-c(3,j)
 577               distal=dsqrt(x12*x12+y12*y12+z12*z12)
 578 c              write (iout,*) k,i,j,distal,dist2_cut
 579
 580             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
 581      &            .and. distal.le.dist2_cut ) then
 582
 583               ii=ii+1
 584               ii_in_use(ii)=1
 585               l_homo(k,ii)=.true.
 586
 587 c             write (iout,*) "k",k
 588 c             write (iout,*) "i",i," j",j," constr_homology",
 589 c    &                       constr_homology
 590               ires_homo(ii)=i
 591               jres_homo(ii)=j
 592               odl(k,ii)=distal
 593               if (read2sigma) then
 594                 sigma_odl(k,ii)=0
 595                 do ik=i,j
 596                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
 597                 enddo
 598                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
 599                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
 600      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 601               else
 602                 if (odl(k,ii).le.dist_cut) then
 603                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
 604                 else
 605 #ifdef OLDSIGMA
 606                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 607      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
 608 #else
 609                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 610      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 611 #endif
 612                 endif
 613               endif
 614               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
 615             else
 616               ii=ii+1
 617 c              l_homo(k,ii)=.false.
 618             endif
 619             enddo
 620           enddo
 621         lim_odl=ii
 622         endif
 623 c
 624 c     Theta, dihedral and SC retraints
 625 c
 626         if (waga_angle.gt.0.0d0) then
 627           do i = nnt+3,nct
 628             if (idomain(k,i).eq.0) then
 629 c               sigma_dih(k,i)=0.0
 630                cycle
 631             endif
 632             dih(k,i)=phiref(i)
 633             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
 634      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
 635 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
 636 c     &       " sigma_dihed",sigma_dih(k,i)
 637             if (sigma_dih(k,i).ne.0)
 638      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
 639           enddo
 640           lim_dih=nct-nnt-2
 641         endif
 642
 643         if (waga_theta.gt.0.0d0) then
 644           do i = nnt+2,nct
 645              if (idomain(k,i).eq.0) then
 646 c              sigma_theta(k,i)=0.0
 647               cycle
 648              endif
 649              thetatpl(k,i)=thetaref(i)
 650              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
 651      &                        rescore(k,i-2))/3.0
 652              if (sigma_theta(k,i).ne.0)
 653      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 654           enddo
 655         endif
 656
 657         if (waga_d.gt.0.0d0) then
 658           do i = nnt,nct
 659                if (itype(i).eq.10) cycle
 660                if (idomain(k,i).eq.0 ) then
 661 c                  sigma_d(k,i)=0.0
 662                   cycle
 663                endif
 664                xxtpl(k,i)=xxref(i)
 665                yytpl(k,i)=yyref(i)
 666                zztpl(k,i)=zzref(i)
 667                sigma_d(k,i)=rescore(k,i)
 668                if (sigma_d(k,i).ne.0)
 669      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
 670                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
 671           enddo
 672         endif
 673       enddo ! l
 674       enddo ! ll
 675 c
 676 c remove distance restraints not used in any model from the list
 677 c shift data in all arrays
 678 c
 679       if (waga_dist.ne.0.0d0) then
 680         ii=0
 681         liiflag=.true.
 682         do i=nnt,nct-2
 683          do j=i+2,nct
 684           ii=ii+1
 685           if (ii_in_use(ii).eq.0.and.liiflag) then
 686             liiflag=.false.
 687             iistart=ii
 688           endif
 689           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
 690      &                   .not.liiflag.and.ii.eq.lim_odl) then
 691              if (ii.eq.lim_odl) then
 692               iishift=ii-iistart+1
 693              else
 694               iishift=ii-iistart
 695              endif
 696              liiflag=.true.
 697              do ki=iistart,lim_odl-iishift
 698               ires_homo(ki)=ires_homo(ki+iishift)
 699               jres_homo(ki)=jres_homo(ki+iishift)
 700               ii_in_use(ki)=ii_in_use(ki+iishift)
 701               do k=1,constr_homology
 702                odl(k,ki)=odl(k,ki+iishift)
 703                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
 704                l_homo(k,ki)=l_homo(k,ki+iishift)
 705               enddo
 706              enddo
 707              ii=ii-iishift
 708              lim_odl=lim_odl-iishift
 709           endif
 710          enddo
 711         enddo
 712       endif
 713
 714       return
 715    10 stop "Error in fragment file"
 716       end
 717