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

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