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