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