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