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

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