source/cluster/wham/src-HCD-5D/read_constr_homology.F

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