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