source/wham/src-M/molread_zs.F

   1       subroutine molread(*)
   2 C
   3 C Read molecular data.
   4 C
   5       implicit real*8 (a-h,o-z)
   6       include 'DIMENSIONS'
   7       include 'DIMENSIONS.ZSCOPT'
   8       include 'DIMENSIONS.FREE'
   9       include 'COMMON.IOUNITS'
  10       include 'COMMON.GEO'
  11       include 'COMMON.VAR'
  12       include 'COMMON.INTERACT'
  13       include 'COMMON.LOCAL'
  14       include 'COMMON.NAMES'
  15       include 'COMMON.CHAIN'
  16       include 'COMMON.FFIELD'
  17       include 'COMMON.SBRIDGE'
  18       include 'COMMON.TORCNSTR'
  19       include 'COMMON.CONTROL'
  20       character*4 sequence(maxres)
  21       integer rescode
  22       double precision x(maxvar)
  23       character*320 controlcard,ucase
  24       dimension itype_pdb(maxres)
  25       logical seq_comp
  26       call card_concat(controlcard,.true.)
  27       call reada(controlcard,'SCAL14',scal14,0.4d0)
  28       call reada(controlcard,'SCALSCP',scalscp,1.0d0)
  29       call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
  30       call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
  31       r0_corr=cutoff_corr-delt_corr
  32       call readi(controlcard,"NRES",nres,0)
  33       iscode=index(controlcard,"ONE_LETTER")
  34       if (nres.le.0) then
  35         write (iout,*) "Error: no residues in molecule"
  36         return1
  37       endif
  38       if (nres.gt.maxres) then
  39         write (iout,*) "Error: too many residues",nres,maxres
  40       endif
  41       write(iout,*) 'nres=',nres
  42 C Read sequence of the protein
  43       if (iscode.gt.0) then
  44         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
  45       else
  46         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
  47       endif
  48 C Convert sequence to numeric code
  49       do i=1,nres
  50         itype(i)=rescode(i,sequence(i),iscode)
  51       enddo
  52       write (iout,*) "Numeric code:"
  53       write (iout,'(20i4)') (itype(i),i=1,nres)
  54       do i=1,nres-1
  55 #ifdef PROCOR
  56         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
  57 #else
  58         if (itype(i).eq.ntyp1) then
  59 #endif
  60           itel(i)=0
  61 #ifdef PROCOR
  62         else if (iabs(itype(i+1)).ne.20) then
  63 #else
  64         else if (iabs(itype(i)).ne.20) then
  65 #endif
  66           itel(i)=1
  67         else
  68           itel(i)=2
  69         endif
  70       enddo
  71        write (iout,*) "ITEL"
  72        do i=1,nres-1
  73          write (iout,*) i,itype(i),itel(i)
  74        enddo
  75       call read_bridge
  76
  77       if (with_dihed_constr) then
  78
  79       read (inp,*) ndih_constr
  80       if (ndih_constr.gt.0) then
  81 C        read (inp,*) ftors
  82 C        write (iout,*) 'FTORS',ftors
  83         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
  84      &   i=1,ndih_constr)
  85         write (iout,*)
  86      &   'There are',ndih_constr,' constraints on phi angles.'
  87         do i=1,ndih_constr
  88           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
  89      &    ftors(i)
  90         enddo
  91         do i=1,ndih_constr
  92           phi0(i)=deg2rad*phi0(i)
  93           drange(i)=deg2rad*drange(i)
  94         enddo
  95       endif
  96
  97       endif
  98       if (with_theta_constr) then
  99 C with_theta_constr is keyword allowing for occurance of theta constrains
 100       read (inp,*) ntheta_constr
 101 C ntheta_constr is the number of theta constrains
 102       if (ntheta_constr.gt.0) then
 103 C        read (inp,*) ftors
 104         read (inp,*) (itheta_constr(i),theta_constr0(i),
 105      &  theta_drange(i),for_thet_constr(i),
 106      &  i=1,ntheta_constr)
 107 C the above code reads from 1 to ntheta_constr
 108 C itheta_constr(i) residue i for which is theta_constr
 109 C theta_constr0 the global minimum value
 110 C theta_drange is range for which there is no energy penalty
 111 C for_thet_constr is the force constant for quartic energy penalty
 112 C E=k*x**4
 113 C        if(me.eq.king.or..not.out1file)then
 114          write (iout,*)
 115      &   'There are',ntheta_constr,' constraints on phi angles.'
 116          do i=1,ntheta_constr
 117           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
 118      &    theta_drange(i),
 119      &    for_thet_constr(i)
 120          enddo
 121 C        endif
 122         do i=1,ntheta_constr
 123           theta_constr0(i)=deg2rad*theta_constr0(i)
 124           theta_drange(i)=deg2rad*theta_drange(i)
 125         enddo
 126 C        if(me.eq.king.or..not.out1file)
 127 C     &   write (iout,*) 'FTORS',ftors
 128 C        do i=1,ntheta_constr
 129 C          ii = itheta_constr(i)
 130 C          thetabound(1,ii) = phi0(i)-drange(i)
 131 C          thetabound(2,ii) = phi0(i)+drange(i)
 132 C        enddo
 133       endif ! ntheta_constr.gt.0
 134       endif! with_theta_constr
 135       nnt=1
 136       nct=nres
 137       if (itype(1).eq.ntyp1) nnt=2
 138       if (itype(nres).eq.ntyp1) nct=nct-1
 139       write(iout,*) 'NNT=',NNT,' NCT=',NCT
 140       write (iout,*) "calling read_saxs_consrtr",nsaxs
 141       if (nsaxs.gt.0) call read_saxs_constr
 142
 143       if (constr_homology.gt.0) then
 144 c       write (iout,*) "About to call read_constr_homology"
 145 c       call flush(iout)
 146         call read_constr_homology
 147        write (iout,*) "Exit read_constr_homology"
 148        call flush(iout)
 149 cref        if (indpdb.gt.0 .or. pdbref) then
 150 cref          do i=1,2*nres
 151 cref            do j=1,3
 152 cref              c(j,i)=crefjlee(j,i)
 153 cref              cref(j,i)=crefjlee(j,i)
 154 cref            enddo
 155 cref          enddo
 156 cref        endif
 157       else
 158         homol_nset=0
 159       endif
 160
 161
 162       call setup_var
 163       call init_int_table
 164       if (ns.gt.0) then
 165         write (iout,'(/a,i3,a)') 'The chain contains',ns,
 166      &  ' disulfide-bridging cysteines.'
 167         write (iout,'(20i4)') (iss(i),i=1,ns)
 168        if (dyn_ss) then
 169           write(iout,*)"Running with dynamic disulfide-bond formation"
 170        else
 171         write (iout,'(/a/)') 'Pre-formed links are:'
 172         do i=1,nss
 173           i1=ihpb(i)-nres
 174           i2=jhpb(i)-nres
 175           it1=itype(i1)
 176           it2=itype(i2)
 177          write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 178      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
 179      &    dhpb(i),ebr,forcon(i)
 180         enddo
 181       endif
 182       endif
 183       write (iout,'(a)')
 184       if (ns.gt.0.and.dyn_ss) then
 185           do i=nss+1,nhpb
 186             ihpb(i-nss)=ihpb(i)
 187             jhpb(i-nss)=jhpb(i)
 188             forcon(i-nss)=forcon(i)
 189             dhpb(i-nss)=dhpb(i)
 190           enddo
 191           nhpb=nhpb-nss
 192           nss=0
 193           call hpb_partition
 194           do i=1,ns
 195             dyn_ss_mask(iss(i))=.true.
 196           enddo
 197       endif
 198       return
 199       end
 200 c-----------------------------------------------------------------------------
 201       logical function seq_comp(itypea,itypeb,length)
 202       implicit none
 203       integer length,itypea(length),itypeb(length)
 204       integer i
 205       do i=1,length
 206         if (itypea(i).ne.itypeb(i)) then
 207           seq_comp=.false.
 208           return
 209         endif
 210       enddo
 211       seq_comp=.true.
 212       return
 213       end
 214 c-----------------------------------------------------------------------------
 215       subroutine read_bridge
 216 C Read information about disulfide bridges.
 217       implicit real*8 (a-h,o-z)
 218       include 'DIMENSIONS'
 219       include 'DIMENSIONS.ZSCOPT'
 220       include 'COMMON.IOUNITS'
 221       include 'COMMON.GEO'
 222       include 'COMMON.VAR'
 223       include 'COMMON.INTERACT'
 224       include 'COMMON.NAMES'
 225       include 'COMMON.CHAIN'
 226       include 'COMMON.FFIELD'
 227       include 'COMMON.SBRIDGE'
 228 C Read bridging residues.
 229       read (inp,*) ns,(iss(i),i=1,ns)
 230       print *,'ns=',ns
 231       write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 232 C Check whether the specified bridging residues are cystines.
 233       do i=1,ns
 234         if (itype(iss(i)).ne.1) then
 235           write (iout,'(2a,i3,a)')
 236      &   'Do you REALLY think that the residue ',
 237      &    restyp(itype(iss(i))),i,
 238      &   ' can form a disulfide bridge?!!!'
 239           write (*,'(2a,i3,a)')
 240      &   'Do you REALLY think that the residue ',
 241      &    restyp(itype(iss(i))),i,
 242      &   ' can form a disulfide bridge?!!!'
 243          stop
 244         endif
 245       enddo
 246 C Read preformed bridges.
 247       if (ns.gt.0) then
 248       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 249       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 250       if (nss.gt.0) then
 251         nhpb=nss
 252 C Check if the residues involved in bridges are in the specified list of
 253 C bridging residues.
 254         do i=1,nss
 255           do j=1,i-1
 256             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 257      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 258               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 259      &      ' contains residues present in other pairs.'
 260               write (*,'(a,i3,a)') 'Disulfide pair',i,
 261      &      ' contains residues present in other pairs.'
 262               stop
 263             endif
 264           enddo
 265           do j=1,ns
 266             if (ihpb(i).eq.iss(j)) goto 10
 267           enddo
 268           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 269    10     continue
 270           do j=1,ns
 271             if (jhpb(i).eq.iss(j)) goto 20
 272           enddo
 273           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 274    20     continue
 275           dhpb(i)=dbr
 276           forcon(i)=fbr
 277         enddo
 278         do i=1,nss
 279           ihpb(i)=ihpb(i)+nres
 280           jhpb(i)=jhpb(i)+nres
 281         enddo
 282       endif
 283       endif
 284       return
 285       end
 286 c------------------------------------------------------------------------------
 287       subroutine read_angles(kanal,iscor,energ,iprot,*)
 288       implicit real*8 (a-h,o-z)
 289       include 'DIMENSIONS'
 290       include 'DIMENSIONS.ZSCOPT'
 291       include 'COMMON.INTERACT'
 292       include 'COMMON.SBRIDGE'
 293       include 'COMMON.GEO'
 294       include 'COMMON.VAR'
 295       include 'COMMON.CHAIN'
 296       include 'COMMON.IOUNITS'
 297       character*80 lineh
 298       read(kanal,'(a80)',end=10,err=10) lineh
 299       read(lineh(:5),*,err=8) ic
 300       read(lineh(6:),*,err=8) energ
 301       goto 9
 302     8 ic=1
 303       print *,'error, assuming e=1d10',lineh
 304       energ=1d10
 305       nss=0
 306     9 continue
 307       read(lineh(18:),*,end=10,err=10) nss
 308       IF (NSS.LT.9) THEN
 309         read (lineh(20:),*,end=10,err=10)
 310      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
 311       ELSE
 312         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
 313         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
 314      &    I=9,NSS),iscor
 315       ENDIF
 316 c      print *,"energy",energ," iscor",iscor
 317       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 318       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 319       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 320       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 321       do i=1,nres
 322         theta(i)=deg2rad*theta(i)
 323         phi(i)=deg2rad*phi(i)
 324         alph(i)=deg2rad*alph(i)
 325         omeg(i)=deg2rad*omeg(i)
 326       enddo
 327       return
 328    10 return1
 329       end
 330 c-------------------------------------------------------------------------------
 331       subroutine read_saxs_constr
 332       implicit real*8 (a-h,o-z)
 333       include 'DIMENSIONS'
 334       include 'DIMENSIONS.ZSCOPT'
 335       include 'DIMENSIONS.FREE'
 336 #ifdef MPI
 337       include 'mpif.h'
 338 #endif
 339       include 'COMMON.CONTROL'
 340       include 'COMMON.CHAIN'
 341       include 'COMMON.IOUNITS'
 342       include 'COMMON.SBRIDGE'
 343       double precision cm(3)
 344 c      read(inp,*) nsaxs
 345       write (iout,*) "Calling read_saxs nsaxs",nsaxs
 346       call flush(iout)
 347       if (saxs_mode.eq.0) then
 348 c SAXS distance distribution
 349       do i=1,nsaxs
 350         read(inp,*) distsaxs(i),Psaxs(i)
 351       enddo
 352       Cnorm = 0.0d0
 353       do i=1,nsaxs
 354         Cnorm = Cnorm + Psaxs(i)
 355       enddo
 356       write (iout,*) "Cnorm",Cnorm
 357       do i=1,nsaxs
 358         Psaxs(i)=Psaxs(i)/Cnorm
 359       enddo
 360       write (iout,*) "Normalized distance distribution from SAXS"
 361       do i=1,nsaxs
 362         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
 363       enddo
 364       Wsaxs0=0.0d0
 365       do i=1,nsaxs
 366         Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
 367       enddo
 368       write (iout,*) "Wsaxs0",Wsaxs0
 369       else
 370 c SAXS "spheres".
 371       do i=1,nsaxs
 372         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
 373       enddo
 374       do j=1,3
 375         cm(j)=0.0d0
 376       enddo
 377       do i=1,nsaxs
 378         do j=1,3
 379           cm(j)=cm(j)+Csaxs(j,i)
 380         enddo
 381       enddo
 382       do j=1,3
 383         cm(j)=cm(j)/nsaxs
 384       enddo
 385       do i=1,nsaxs
 386         do j=1,3
 387           Csaxs(j,i)=Csaxs(j,i)-cm(j)
 388         enddo
 389       enddo
 390       write (iout,*) "SAXS sphere coordinates"
 391       do i=1,nsaxs
 392         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
 393       enddo
 394       endif
 395       return
 396       end
 397 c====-------------------------------------------------------------------
 398       subroutine read_constr_homology
 399
 400       include 'DIMENSIONS'
 401       include 'DIMENSIONS.ZSCOPT'
 402       include 'DIMENSIONS.FREE'
 403 #ifdef MPI
 404       include 'mpif.h'
 405 #endif
 406       include 'COMMON.SETUP'
 407       include 'COMMON.CONTROL'
 408       include 'COMMON.CHAIN'
 409       include 'COMMON.IOUNITS'
 410       include 'COMMON.GEO'
 411       include 'COMMON.INTERACT'
 412       include 'include_unres/COMMON.NAMES'
 413       include 'COMMON.HOMRESTR'
 414 c
 415 c For new homol impl
 416 c
 417       include 'COMMON.VAR'
 418 c     include 'include_unres/COMMON.VAR'
 419 c
 420
 421 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
 422 c    &                 dist_cut
 423 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
 424 c    &    sigma_odl_temp(maxres,maxres,max_template)
 425       character*2 kic2
 426       character*24 model_ki_dist, model_ki_angle
 427       character*500 controlcard
 428       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
 429       integer idomain(max_template,maxres)
 430       logical lprn /.true./
 431       integer ilen
 432       external ilen
 433       logical unres_pdb,liiflag
 434 c
 435 c     FP - Nov. 2014 Temporary specifications for new vars
 436 c
 437       double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
 438      &                 rescore3_tmp
 439       double precision, dimension (max_template,maxres) :: rescore
 440       double precision, dimension (max_template,maxres) :: rescore2
 441       double precision, dimension (max_template,maxres) :: rescore3
 442       character*24 tpl_k_rescore
 443 c -----------------------------------------------------------------
 444 c Reading multiple PDB ref structures and calculation of retraints
 445 c not using pre-computed ones stored in files model_ki_{dist,angle}
 446 c FP (Nov., 2014)
 447 c -----------------------------------------------------------------
 448 c
 449 c
 450 c Alternative: reading from input
 451       call card_concat(controlcard,.true.)
 452       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
 453       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
 454       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
 455       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
 456       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
 457       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
 458       dist1cut=(index(controlcard,'DIST1CUT').gt.0)
 459       call readi(controlcard,"HOMOL_SET",homol_nset,1)
 460       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
 461       call readi(controlcard,"IHSET",ihset,1)
 462       if (homol_nset.gt.1)then
 463          call card_concat(controlcard,.true.)
 464          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
 465          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 466           write(iout,*) "iset homology_weight "
 467 c         do i=1,homol_nset
 468 c          write(iout,*) i,waga_homology(i)
 469 c         enddo
 470          endif
 471          iset=mod(kolor,homol_nset)+1
 472       else
 473        iset=1
 474        waga_homology(1)=1.0
 475       endif
 476 c     write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
 477
 478 cd      write (iout,*) "nnt",nnt," nct",nct
 479 cd      call flush(iout)
 480
 481
 482       lim_odl=0
 483       lim_dih=0
 484 c
 485 c  New
 486 c
 487 c
 488 c  Reading HM global scores (prob not required)
 489 c
 490       do i = nnt,nct
 491         do k=1,constr_homology
 492          idomain(k,i)=0
 493         enddo
 494       enddo
 495 c     open (4,file="HMscore")
 496 c     do k=1,constr_homology
 497 c       read (4,*,end=521) hmscore_tmp
 498 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
 499 c       write(*,*) "Model", k, ":", hmscore(k)
 500 c     enddo
 501 c521  continue
 502
 503       ii=0
 504       do i = nnt,nct-2
 505         do j=i+2,nct
 506         ii=ii+1
 507         ii_in_use(ii)=0
 508         enddo
 509       enddo
 510 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 511
 512       if (read_homol_frag) then
 513         call read_klapaucjusz
 514       else
 515
 516       do k=1,constr_homology
 517
 518         read(inp,'(a)') pdbfile
 519 c  Next stament causes error upon compilation (?)
 520 c       if(me.eq.king.or. .not. out1file)
 521 c         write (iout,'(2a)') 'PDB data will be read from file ',
 522 c    &   pdbfile(:ilen(pdbfile))
 523          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
 524      &  pdbfile(:ilen(pdbfile))
 525         open(ipdbin,file=pdbfile,status='old',err=33)
 526         goto 34
 527   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
 528      &  pdbfile(:ilen(pdbfile))
 529         stop
 530   34    continue
 531 c        print *,'Begin reading pdb data'
 532 c
 533 c Files containing res sim or local scores (former containing sigmas)
 534 c
 535
 536         write(kic2,'(bz,i2.2)') k
 537
 538         tpl_k_rescore="template"//kic2//".sco"
 539
 540         unres_pdb=.false.
 541         call readpdb_template(k)
 542 cref        do i=1,2*nres
 543 cref          do j=1,3
 544 cref            crefjlee(j,i)=c(j,i)
 545 cref          enddo
 546 cref        enddo
 547 #ifdef DEBUG
 548         do i=1,nres
 549           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
 550      &      (crefjlee(j,i+nres),j=1,3)
 551         enddo
 552 #endif
 553         write (iout,*) "read_constr_homology: after reading pdb file"
 554         call flush(iout)
 555
 556 c
 557 c     Distance restraints
 558 c
 559 c          ... --> odl(k,ii)
 560 C Copy the coordinates from reference coordinates (?)
 561         do i=1,2*nres
 562           do j=1,3
 563 c            c(j,i)=cref(j,i)
 564 c           write (iout,*) "c(",j,i,") =",c(j,i)
 565           enddo
 566         enddo
 567 c
 568 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
 569 c
 570 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
 571           open (ientin,file=tpl_k_rescore,status='old')
 572           if (nnt.gt.1) rescore(k,1)=0.0d0
 573           do irec=nnt,nct ! loop for reading res sim
 574             if (read2sigma) then
 575              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
 576      &                                rescore3_tmp,idomain_tmp
 577              i_tmp=i_tmp+nnt-1
 578              idomain(k,i_tmp)=idomain_tmp
 579              rescore(k,i_tmp)=rescore_tmp
 580              rescore2(k,i_tmp)=rescore2_tmp
 581              rescore3(k,i_tmp)=rescore3_tmp
 582              write(iout,'(a7,i5,3f10.5,i5)') "rescore",
 583      &                      i_tmp,rescore2_tmp,rescore_tmp,
 584      &                                rescore3_tmp,idomain_tmp
 585             else
 586              idomain(k,irec)=1
 587              read (ientin,*,end=1401) rescore_tmp
 588
 589 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
 590              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
 591 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
 592             endif
 593           enddo
 594  1401   continue
 595         close (ientin)
 596         if (waga_dist.ne.0.0d0) then
 597           ii=0
 598           do i = nnt,nct-2
 599             do j=i+2,nct
 600
 601               x12=c(1,i)-c(1,j)
 602               y12=c(2,i)-c(2,j)
 603               z12=c(3,i)-c(3,j)
 604               distal=dsqrt(x12*x12+y12*y12+z12*z12)
 605 c              write (iout,*) k,i,j,distal,dist2_cut
 606            if (dist1cut .and. k.gt.1) then
 607               ii=ii+1
 608               if (l_homo(1,ii)) then
 609                 ii_in_use(ii)=1
 610                 l_homo(k,ii)=.true.
 611                 ires_homo(ii)=i
 612                 jres_homo(ii)=j
 613                 odl(k,ii)=distal
 614                 sigma_odl(k,ii)=sigma_odl(1,ii)
 615               else
 616                 l_homo(k,ii)=.false.
 617               endif
 618            else
 619             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
 620      &            .and. distal.le.dist2_cut ) then
 621
 622               ii=ii+1
 623               ii_in_use(ii)=1
 624               l_homo(k,ii)=.true.
 625
 626 c             write (iout,*) "k",k
 627 c             write (iout,*) "i",i," j",j," constr_homology",
 628 c    &                       constr_homology
 629               ires_homo(ii)=i
 630               jres_homo(ii)=j
 631               odl(k,ii)=distal
 632               if (read2sigma) then
 633                 sigma_odl(k,ii)=0
 634                 do ik=i,j
 635                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
 636                 enddo
 637                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
 638                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
 639      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 640               else
 641                 if (odl(k,ii).le.dist_cut) then
 642                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
 643                 else
 644 #ifdef OLDSIGMA
 645                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 646      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
 647 #else
 648                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 649      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 650 #endif
 651                 endif
 652               endif
 653               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
 654             else
 655               ii=ii+1
 656               l_homo(k,ii)=.false.
 657             endif
 658            endif
 659             enddo
 660           enddo
 661         lim_odl=ii
 662         endif
 663 c
 664 c     Theta, dihedral and SC retraints
 665 c
 666         if (waga_angle.gt.0.0d0) then
 667 c         open (ientin,file=tpl_k_sigma_dih,status='old')
 668 c         do irec=1,maxres-3 ! loop for reading sigma_dih
 669 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
 670 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
 671 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
 672 c    &                            sigma_dih(k,i+nnt-1)
 673 c         enddo
 674 c1402   continue
 675 c         close (ientin)
 676           do i = nnt+3,nct
 677             if (idomain(k,i).eq.0) then
 678                sigma_dih(k,i)=0.0
 679                cycle
 680             endif
 681             dih(k,i)=phiref(i) ! right?
 682 c           read (ientin,*) sigma_dih(k,i) ! original variant
 683 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
 684 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
 685 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 686 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
 687 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
 688
 689             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
 690      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
 691 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
 692 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
 693 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
 694 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 695 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
 696 c   Instead of res sim other local measure of b/b str reliability possible
 697             if (sigma_dih(k,i).ne.0)
 698      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
 699 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
 700           enddo
 701           lim_dih=nct-nnt-2
 702         endif
 703
 704         if (waga_theta.gt.0.0d0) then
 705 c         open (ientin,file=tpl_k_sigma_theta,status='old')
 706 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
 707 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
 708 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
 709 c    &                              sigma_theta(k,i+nnt-1)
 710 c         enddo
 711 c1403   continue
 712 c         close (ientin)
 713
 714           do i = nnt+2,nct ! right? without parallel.
 715 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
 716 c         do i=ithet_start,ithet_end ! with FG parallel.
 717              if (idomain(k,i).eq.0) then
 718               sigma_theta(k,i)=0.0
 719               cycle
 720              endif
 721              thetatpl(k,i)=thetaref(i)
 722 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
 723 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
 724 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 725 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
 726 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
 727              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
 728      &                        rescore(k,i-2))/3.0
 729 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
 730              if (sigma_theta(k,i).ne.0)
 731      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 732
 733 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 734 c                             rescore(k,i-2) !  right expression ?
 735 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
 736           enddo
 737         endif
 738
 739         if (waga_d.gt.0.0d0) then
 740 c       open (ientin,file=tpl_k_sigma_d,status='old')
 741 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
 742 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
 743 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
 744 c    &                          sigma_d(k,i+nnt-1)
 745 c         enddo
 746 c1404   continue
 747
 748           do i = nnt,nct ! right? without parallel.
 749 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
 750 c         do i=loc_start,loc_end ! with FG parallel.
 751                if (itype(i).eq.10) cycle
 752                if (idomain(k,i).eq.0 ) then
 753                   sigma_d(k,i)=0.0
 754                   cycle
 755                endif
 756                xxtpl(k,i)=xxref(i)
 757                yytpl(k,i)=yyref(i)
 758                zztpl(k,i)=zzref(i)
 759 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
 760 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
 761 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
 762 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
 763                sigma_d(k,i)=rescore3(k,i) !  right expression ?
 764                if (sigma_d(k,i).ne.0)
 765      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
 766
 767 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
 768 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
 769 c              read (ientin,*) sigma_d(k,i) ! 1st variant
 770           enddo
 771         endif
 772       enddo
 773 c
 774 c remove distance restraints not used in any model from the list
 775 c shift data in all arrays
 776 c
 777       if (waga_dist.ne.0.0d0) then
 778         ii=0
 779         liiflag=.true.
 780         do i=nnt,nct-2
 781          do j=i+2,nct
 782           ii=ii+1
 783           if (ii_in_use(ii).eq.0.and.liiflag) then
 784             liiflag=.false.
 785             iistart=ii
 786           endif
 787           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
 788      &                   .not.liiflag.and.ii.eq.lim_odl) then
 789              if (ii.eq.lim_odl) then
 790               iishift=ii-iistart+1
 791              else
 792               iishift=ii-iistart
 793              endif
 794              liiflag=.true.
 795              do ki=iistart,lim_odl-iishift
 796               ires_homo(ki)=ires_homo(ki+iishift)
 797               jres_homo(ki)=jres_homo(ki+iishift)
 798               ii_in_use(ki)=ii_in_use(ki+iishift)
 799               do k=1,constr_homology
 800                odl(k,ki)=odl(k,ki+iishift)
 801                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
 802                l_homo(k,ki)=l_homo(k,ki+iishift)
 803               enddo
 804              enddo
 805              ii=ii-iishift
 806              lim_odl=lim_odl-iishift
 807           endif
 808          enddo
 809         enddo
 810       endif
 811
 812       endif ! .not. klapaucjusz
 813
 814       if (constr_homology.gt.0) call homology_partition
 815       if (constr_homology.gt.0) call init_int_table
 816 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
 817 cd     & "lim_xx=",lim_xx
 818 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
 819 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
 820 c
 821 c Print restraints
 822 c
 823       if (.not.lprn) return
 824 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 825       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 826        write (iout,*) "Distance restraints from templates"
 827        do ii=1,lim_odl
 828        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
 829      &  ii,ires_homo(ii),jres_homo(ii),
 830      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
 831      &  ki=1,constr_homology)
 832        enddo
 833        write (iout,*) "Dihedral angle restraints from templates"
 834        do i=nnt+3,nct
 835         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
 836      &      (rad2deg*dih(ki,i),
 837      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
 838        enddo
 839        write (iout,*) "Virtual-bond angle restraints from templates"
 840        do i=nnt+2,nct
 841         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
 842      &      (rad2deg*thetatpl(ki,i),
 843      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
 844        enddo
 845        write (iout,*) "SC restraints from templates"
 846        do i=nnt,nct
 847         write(iout,'(i5,100(4f8.2,4x))') i,
 848      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
 849      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
 850        enddo
 851       endif
 852 c -----------------------------------------------------------------
 853       return
 854       end
 855 c----------------------------------------------------------------------
 856       subroutine read_klapaucjusz
 857
 858       include 'DIMENSIONS'
 859       include 'DIMENSIONS.ZSCOPT'
 860       include 'DIMENSIONS.FREE'
 861 #ifdef MPI
 862       include 'mpif.h'
 863 #endif
 864       include 'COMMON.SETUP'
 865       include 'COMMON.CONTROL'
 866       include 'COMMON.CHAIN'
 867       include 'COMMON.IOUNITS'
 868       include 'COMMON.GEO'
 869       include 'COMMON.INTERACT'
 870       include 'COMMON.NAMES'
 871       include 'COMMON.HOMRESTR'
 872       character*256 fragfile
 873       integer ninclust(maxclust),inclust(max_template,maxclust),
 874      &  nresclust(maxclust),iresclust(maxres,maxclust)
 875
 876       character*2 kic2
 877       character*24 model_ki_dist, model_ki_angle
 878       character*500 controlcard
 879       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
 880       integer idomain(max_template,maxres)
 881       logical lprn /.true./
 882       integer ilen
 883       external ilen
 884       logical unres_pdb,liiflag
 885 c
 886 c
 887       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
 888       double precision, dimension (max_template,maxres) :: rescore
 889       double precision, dimension (max_template,maxres) :: rescore2
 890       character*24 tpl_k_rescore
 891
 892 c
 893 c For new homol impl
 894 c
 895       include 'COMMON.VAR'
 896 c
 897       double precision chomo(3,maxres2+2,max_template)
 898       call getenv("FRAGFILE",fragfile)
 899       open(ientin,file=fragfile,status="old",err=10)
 900       read(ientin,*) constr_homology,nclust
 901       l_homo = .false.
 902       sigma_theta=0.0
 903       sigma_d=0.0
 904       sigma_dih=0.0
 905 c Read pdb files
 906       do k=1,constr_homology
 907         read(ientin,'(a)') pdbfile
 908         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
 909      &  pdbfile(:ilen(pdbfile))
 910         open(ipdbin,file=pdbfile,status='old',err=33)
 911         goto 34
 912   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
 913      &  pdbfile(:ilen(pdbfile))
 914         stop
 915   34    continue
 916         unres_pdb=.false.
 917         call readpdb_template(k)
 918         do i=1,2*nres
 919           do j=1,3
 920             chomo(j,i,k)=c(j,i)
 921           enddo
 922         enddo
 923         do i=1,nres
 924           rescore(k,i)=0.2d0
 925           rescore2(k,i)=1.0d0
 926         enddo
 927       enddo
 928 c Read clusters
 929       do i=1,nclust
 930         read(ientin,*) ninclust(i),nresclust(i)
 931         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
 932         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
 933       enddo
 934 c
 935 c Loop over clusters
 936 c
 937       do l=1,nclust
 938         do ll = 1,ninclust(l)
 939
 940         k = inclust(ll,l)
 941         do i=1,nres
 942           idomain(k,i)=0
 943         enddo
 944         do i=1,nresclust(l)
 945           if (nnt.gt.1)  then
 946             idomain(k,iresclust(i,l)+1) = 1
 947           else
 948             idomain(k,iresclust(i,l)) = 1
 949           endif
 950         enddo
 951 c
 952 c     Distance restraints
 953 c
 954 c          ... --> odl(k,ii)
 955 C Copy the coordinates from reference coordinates (?)
 956         do i=1,2*nres
 957           do j=1,3
 958             c(j,i)=chomo(j,i,k)
 959 c           write (iout,*) "c(",j,i,") =",c(j,i)
 960           enddo
 961         enddo
 962         call int_from_cart1(.false.)
 963         call int_from_cart(.true.,.false.)
 964         call sc_loc_geom(.false.)
 965         do i=1,nres
 966           thetaref(i)=theta(i)
 967           phiref(i)=phi(i)
 968         enddo
 969         if (waga_dist.ne.0.0d0) then
 970           ii=0
 971           do i = nnt,nct-2
 972             do j=i+2,nct
 973
 974               x12=c(1,i)-c(1,j)
 975               y12=c(2,i)-c(2,j)
 976               z12=c(3,i)-c(3,j)
 977               distal=dsqrt(x12*x12+y12*y12+z12*z12)
 978 c              write (iout,*) k,i,j,distal,dist2_cut
 979
 980             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
 981      &            .and. distal.le.dist2_cut ) then
 982
 983               ii=ii+1
 984               ii_in_use(ii)=1
 985               l_homo(k,ii)=.true.
 986
 987 c             write (iout,*) "k",k
 988 c             write (iout,*) "i",i," j",j," constr_homology",
 989 c    &                       constr_homology
 990               ires_homo(ii)=i
 991               jres_homo(ii)=j
 992               odl(k,ii)=distal
 993               if (read2sigma) then
 994                 sigma_odl(k,ii)=0
 995                 do ik=i,j
 996                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
 997                 enddo
 998                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
 999                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1000      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1001               else
1002                 if (odl(k,ii).le.dist_cut) then
1003                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1004                 else
1005 #ifdef OLDSIGMA
1006                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1007      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1008 #else
1009                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1010      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1011 #endif
1012                 endif
1013               endif
1014               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1015             else
1016               ii=ii+1
1017 c              l_homo(k,ii)=.false.
1018             endif
1019             enddo
1020           enddo
1021         lim_odl=ii
1022         endif
1023 c
1024 c     Theta, dihedral and SC retraints
1025 c
1026         if (waga_angle.gt.0.0d0) then
1027           do i = nnt+3,nct
1028             if (idomain(k,i).eq.0) then
1029 c               sigma_dih(k,i)=0.0
1030                cycle
1031             endif
1032             dih(k,i)=phiref(i)
1033             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1034      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1035 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
1036 c     &       " sigma_dihed",sigma_dih(k,i)
1037             if (sigma_dih(k,i).ne.0)
1038      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1039           enddo
1040           lim_dih=nct-nnt-2
1041         endif
1042
1043         if (waga_theta.gt.0.0d0) then
1044           do i = nnt+2,nct
1045              if (idomain(k,i).eq.0) then
1046 c              sigma_theta(k,i)=0.0
1047               cycle
1048              endif
1049              thetatpl(k,i)=thetaref(i)
1050              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1051      &                        rescore(k,i-2))/3.0
1052              if (sigma_theta(k,i).ne.0)
1053      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1054           enddo
1055         endif
1056
1057         if (waga_d.gt.0.0d0) then
1058           do i = nnt,nct
1059                if (itype(i).eq.10) cycle
1060                if (idomain(k,i).eq.0 ) then
1061 c                  sigma_d(k,i)=0.0
1062                   cycle
1063                endif
1064                xxtpl(k,i)=xxref(i)
1065                yytpl(k,i)=yyref(i)
1066                zztpl(k,i)=zzref(i)
1067                sigma_d(k,i)=rescore(k,i)
1068                if (sigma_d(k,i).ne.0)
1069      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1070                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
1071           enddo
1072         endif
1073       enddo ! l
1074       enddo ! ll
1075 c
1076 c remove distance restraints not used in any model from the list
1077 c shift data in all arrays
1078 c
1079       if (waga_dist.ne.0.0d0) then
1080         ii=0
1081         liiflag=.true.
1082         do i=nnt,nct-2
1083          do j=i+2,nct
1084           ii=ii+1
1085           if (ii_in_use(ii).eq.0.and.liiflag) then
1086             liiflag=.false.
1087             iistart=ii
1088           endif
1089           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1090      &                   .not.liiflag.and.ii.eq.lim_odl) then
1091              if (ii.eq.lim_odl) then
1092               iishift=ii-iistart+1
1093              else
1094               iishift=ii-iistart
1095              endif
1096              liiflag=.true.
1097              do ki=iistart,lim_odl-iishift
1098               ires_homo(ki)=ires_homo(ki+iishift)
1099               jres_homo(ki)=jres_homo(ki+iishift)
1100               ii_in_use(ki)=ii_in_use(ki+iishift)
1101               do k=1,constr_homology
1102                odl(k,ki)=odl(k,ki+iishift)
1103                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1104                l_homo(k,ki)=l_homo(k,ki+iishift)
1105               enddo
1106              enddo
1107              ii=ii-iishift
1108              lim_odl=lim_odl-iishift
1109           endif
1110          enddo
1111         enddo
1112       endif
1113
1114       return
1115    10 stop "Error infragment file"
1116       end
1117 c----------------------------------------------------------------------