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