source/wham/src-HCD/read_constr_homology.F

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