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       call readi(controlcard,"HOMOL_SET",homol_nset,1)
 459       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
 460       call readi(controlcard,"IHSET",ihset,1)
 461       if (homol_nset.gt.1)then
 462          call card_concat(controlcard,.true.)
 463          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
 464          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 465           write(iout,*) "iset homology_weight "
 466 c         do i=1,homol_nset
 467 c          write(iout,*) i,waga_homology(i)
 468 c         enddo
 469          endif
 470          iset=mod(kolor,homol_nset)+1
 471       else
 472        iset=1
 473        waga_homology(1)=1.0
 474       endif
 475 c     write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
 476
 477 cd      write (iout,*) "nnt",nnt," nct",nct
 478 cd      call flush(iout)
 479
 480
 481       lim_odl=0
 482       lim_dih=0
 483 c
 484 c  New
 485 c
 486 c
 487 c  Reading HM global scores (prob not required)
 488 c
 489       do i = nnt,nct
 490         do k=1,constr_homology
 491          idomain(k,i)=0
 492         enddo
 493       enddo
 494 c     open (4,file="HMscore")
 495 c     do k=1,constr_homology
 496 c       read (4,*,end=521) hmscore_tmp
 497 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
 498 c       write(*,*) "Model", k, ":", hmscore(k)
 499 c     enddo
 500 c521  continue
 501
 502       ii=0
 503       do i = nnt,nct-2
 504         do j=i+2,nct
 505         ii=ii+1
 506         ii_in_use(ii)=0
 507         enddo
 508       enddo
 509 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 510
 511       if (read_homol_frag) then
 512         call read_klapaucjusz
 513       else
 514
 515       do k=1,constr_homology
 516
 517         read(inp,'(a)') pdbfile
 518 c  Next stament causes error upon compilation (?)
 519 c       if(me.eq.king.or. .not. out1file)
 520 c         write (iout,'(2a)') 'PDB data will be read from file ',
 521 c    &   pdbfile(:ilen(pdbfile))
 522          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
 523      &  pdbfile(:ilen(pdbfile))
 524         open(ipdbin,file=pdbfile,status='old',err=33)
 525         goto 34
 526   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
 527      &  pdbfile(:ilen(pdbfile))
 528         stop
 529   34    continue
 530 c        print *,'Begin reading pdb data'
 531 c
 532 c Files containing res sim or local scores (former containing sigmas)
 533 c
 534
 535         write(kic2,'(bz,i2.2)') k
 536
 537         tpl_k_rescore="template"//kic2//".sco"
 538
 539         unres_pdb=.false.
 540         call readpdb_template(k)
 541 cref        do i=1,2*nres
 542 cref          do j=1,3
 543 cref            crefjlee(j,i)=c(j,i)
 544 cref          enddo
 545 cref        enddo
 546 #ifdef DEBUG
 547         do i=1,nres
 548           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
 549      &      (crefjlee(j,i+nres),j=1,3)
 550         enddo
 551 #endif
 552         write (iout,*) "read_constr_homology: after reading pdb file"
 553         call flush(iout)
 554
 555 c
 556 c     Distance restraints
 557 c
 558 c          ... --> odl(k,ii)
 559 C Copy the coordinates from reference coordinates (?)
 560         do i=1,2*nres
 561           do j=1,3
 562 c            c(j,i)=cref(j,i)
 563 c           write (iout,*) "c(",j,i,") =",c(j,i)
 564           enddo
 565         enddo
 566 c
 567 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
 568 c
 569 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
 570           open (ientin,file=tpl_k_rescore,status='old')
 571           if (nnt.gt.1) rescore(k,1)=0.0d0
 572           do irec=nnt,nct ! loop for reading res sim
 573             if (read2sigma) then
 574              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
 575      &                                rescore3_tmp,idomain_tmp
 576              i_tmp=i_tmp+nnt-1
 577              idomain(k,i_tmp)=idomain_tmp
 578              rescore(k,i_tmp)=rescore_tmp
 579              rescore2(k,i_tmp)=rescore2_tmp
 580              rescore3(k,i_tmp)=rescore3_tmp
 581              write(iout,'(a7,i5,3f10.5,i5)') "rescore",
 582      &                      i_tmp,rescore2_tmp,rescore_tmp,
 583      &                                rescore3_tmp,idomain_tmp
 584             else
 585              idomain(k,irec)=1
 586              read (ientin,*,end=1401) rescore_tmp
 587
 588 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
 589              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
 590 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
 591             endif
 592           enddo
 593  1401   continue
 594         close (ientin)
 595         if (waga_dist.ne.0.0d0) then
 596           ii=0
 597           do i = nnt,nct-2
 598             do j=i+2,nct
 599
 600               x12=c(1,i)-c(1,j)
 601               y12=c(2,i)-c(2,j)
 602               z12=c(3,i)-c(3,j)
 603               distal=dsqrt(x12*x12+y12*y12+z12*z12)
 604 c              write (iout,*) k,i,j,distal,dist2_cut
 605
 606             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
 607      &            .and. distal.le.dist2_cut ) then
 608
 609               ii=ii+1
 610               ii_in_use(ii)=1
 611               l_homo(k,ii)=.true.
 612
 613 c             write (iout,*) "k",k
 614 c             write (iout,*) "i",i," j",j," constr_homology",
 615 c    &                       constr_homology
 616               ires_homo(ii)=i
 617               jres_homo(ii)=j
 618               odl(k,ii)=distal
 619               if (read2sigma) then
 620                 sigma_odl(k,ii)=0
 621                 do ik=i,j
 622                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
 623                 enddo
 624                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
 625                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
 626      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 627               else
 628                 if (odl(k,ii).le.dist_cut) then
 629                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
 630                 else
 631 #ifdef OLDSIGMA
 632                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 633      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
 634 #else
 635                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 636      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 637 #endif
 638                 endif
 639               endif
 640               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
 641             else
 642               ii=ii+1
 643               l_homo(k,ii)=.false.
 644             endif
 645             enddo
 646           enddo
 647         lim_odl=ii
 648         endif
 649 c
 650 c     Theta, dihedral and SC retraints
 651 c
 652         if (waga_angle.gt.0.0d0) then
 653 c         open (ientin,file=tpl_k_sigma_dih,status='old')
 654 c         do irec=1,maxres-3 ! loop for reading sigma_dih
 655 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
 656 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
 657 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
 658 c    &                            sigma_dih(k,i+nnt-1)
 659 c         enddo
 660 c1402   continue
 661 c         close (ientin)
 662           do i = nnt+3,nct
 663             if (idomain(k,i).eq.0) then
 664                sigma_dih(k,i)=0.0
 665                cycle
 666             endif
 667             dih(k,i)=phiref(i) ! right?
 668 c           read (ientin,*) sigma_dih(k,i) ! original variant
 669 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
 670 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
 671 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 672 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
 673 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
 674
 675             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
 676      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
 677 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
 678 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
 679 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
 680 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 681 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
 682 c   Instead of res sim other local measure of b/b str reliability possible
 683             if (sigma_dih(k,i).ne.0)
 684      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
 685 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
 686           enddo
 687           lim_dih=nct-nnt-2
 688         endif
 689
 690         if (waga_theta.gt.0.0d0) then
 691 c         open (ientin,file=tpl_k_sigma_theta,status='old')
 692 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
 693 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
 694 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
 695 c    &                              sigma_theta(k,i+nnt-1)
 696 c         enddo
 697 c1403   continue
 698 c         close (ientin)
 699
 700           do i = nnt+2,nct ! right? without parallel.
 701 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
 702 c         do i=ithet_start,ithet_end ! with FG parallel.
 703              if (idomain(k,i).eq.0) then
 704               sigma_theta(k,i)=0.0
 705               cycle
 706              endif
 707              thetatpl(k,i)=thetaref(i)
 708 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
 709 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
 710 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 711 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
 712 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
 713              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
 714      &                        rescore(k,i-2))/3.0
 715 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
 716              if (sigma_theta(k,i).ne.0)
 717      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 718
 719 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 720 c                             rescore(k,i-2) !  right expression ?
 721 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
 722           enddo
 723         endif
 724
 725         if (waga_d.gt.0.0d0) then
 726 c       open (ientin,file=tpl_k_sigma_d,status='old')
 727 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
 728 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
 729 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
 730 c    &                          sigma_d(k,i+nnt-1)
 731 c         enddo
 732 c1404   continue
 733
 734           do i = nnt,nct ! right? without parallel.
 735 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
 736 c         do i=loc_start,loc_end ! with FG parallel.
 737                if (itype(i).eq.10) cycle
 738                if (idomain(k,i).eq.0 ) then
 739                   sigma_d(k,i)=0.0
 740                   cycle
 741                endif
 742                xxtpl(k,i)=xxref(i)
 743                yytpl(k,i)=yyref(i)
 744                zztpl(k,i)=zzref(i)
 745 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
 746 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
 747 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
 748 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
 749                sigma_d(k,i)=rescore3(k,i) !  right expression ?
 750                if (sigma_d(k,i).ne.0)
 751      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
 752
 753 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
 754 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
 755 c              read (ientin,*) sigma_d(k,i) ! 1st variant
 756           enddo
 757         endif
 758       enddo
 759 c
 760 c remove distance restraints not used in any model from the list
 761 c shift data in all arrays
 762 c
 763       if (waga_dist.ne.0.0d0) then
 764         ii=0
 765         liiflag=.true.
 766         do i=nnt,nct-2
 767          do j=i+2,nct
 768           ii=ii+1
 769           if (ii_in_use(ii).eq.0.and.liiflag) then
 770             liiflag=.false.
 771             iistart=ii
 772           endif
 773           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
 774      &                   .not.liiflag.and.ii.eq.lim_odl) then
 775              if (ii.eq.lim_odl) then
 776               iishift=ii-iistart+1
 777              else
 778               iishift=ii-iistart
 779              endif
 780              liiflag=.true.
 781              do ki=iistart,lim_odl-iishift
 782               ires_homo(ki)=ires_homo(ki+iishift)
 783               jres_homo(ki)=jres_homo(ki+iishift)
 784               ii_in_use(ki)=ii_in_use(ki+iishift)
 785               do k=1,constr_homology
 786                odl(k,ki)=odl(k,ki+iishift)
 787                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
 788                l_homo(k,ki)=l_homo(k,ki+iishift)
 789               enddo
 790              enddo
 791              ii=ii-iishift
 792              lim_odl=lim_odl-iishift
 793           endif
 794          enddo
 795         enddo
 796       endif
 797
 798       endif ! .not. klapaucjusz
 799
 800       if (constr_homology.gt.0) call homology_partition
 801       if (constr_homology.gt.0) call init_int_table
 802 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
 803 cd     & "lim_xx=",lim_xx
 804 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
 805 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
 806 c
 807 c Print restraints
 808 c
 809       if (.not.lprn) return
 810 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 811       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 812        write (iout,*) "Distance restraints from templates"
 813        do ii=1,lim_odl
 814        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
 815      &  ii,ires_homo(ii),jres_homo(ii),
 816      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
 817      &  ki=1,constr_homology)
 818        enddo
 819        write (iout,*) "Dihedral angle restraints from templates"
 820        do i=nnt+3,nct
 821         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
 822      &      (rad2deg*dih(ki,i),
 823      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
 824        enddo
 825        write (iout,*) "Virtual-bond angle restraints from templates"
 826        do i=nnt+2,nct
 827         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
 828      &      (rad2deg*thetatpl(ki,i),
 829      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
 830        enddo
 831        write (iout,*) "SC restraints from templates"
 832        do i=nnt,nct
 833         write(iout,'(i5,100(4f8.2,4x))') i,
 834      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
 835      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
 836        enddo
 837       endif
 838 c -----------------------------------------------------------------
 839       return
 840       end
 841 c----------------------------------------------------------------------
 842       subroutine read_klapaucjusz
 843
 844       include 'DIMENSIONS'
 845       include 'DIMENSIONS.ZSCOPT'
 846       include 'DIMENSIONS.FREE'
 847 #ifdef MPI
 848       include 'mpif.h'
 849 #endif
 850       include 'COMMON.SETUP'
 851       include 'COMMON.CONTROL'
 852       include 'COMMON.CHAIN'
 853       include 'COMMON.IOUNITS'
 854       include 'COMMON.GEO'
 855       include 'COMMON.INTERACT'
 856       include 'COMMON.NAMES'
 857       include 'COMMON.HOMRESTR'
 858       character*256 fragfile
 859       integer ninclust(maxclust),inclust(max_template,maxclust),
 860      &  nresclust(maxclust),iresclust(maxres,maxclust)
 861
 862       character*2 kic2
 863       character*24 model_ki_dist, model_ki_angle
 864       character*500 controlcard
 865       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
 866       integer idomain(max_template,maxres)
 867       logical lprn /.true./
 868       integer ilen
 869       external ilen
 870       logical unres_pdb,liiflag
 871 c
 872 c
 873       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
 874       double precision, dimension (max_template,maxres) :: rescore
 875       double precision, dimension (max_template,maxres) :: rescore2
 876       character*24 tpl_k_rescore
 877
 878 c
 879 c For new homol impl
 880 c
 881       include 'COMMON.VAR'
 882 c
 883       double precision chomo(3,maxres2+2,max_template)
 884       call getenv("FRAGFILE",fragfile)
 885       open(ientin,file=fragfile,status="old",err=10)
 886       read(ientin,*) constr_homology,nclust
 887       l_homo = .false.
 888       sigma_theta=0.0
 889       sigma_d=0.0
 890       sigma_dih=0.0
 891 c Read pdb files
 892       do k=1,constr_homology
 893         read(ientin,'(a)') pdbfile
 894         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
 895      &  pdbfile(:ilen(pdbfile))
 896         open(ipdbin,file=pdbfile,status='old',err=33)
 897         goto 34
 898   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
 899      &  pdbfile(:ilen(pdbfile))
 900         stop
 901   34    continue
 902         unres_pdb=.false.
 903         call readpdb_template(k)
 904         do i=1,2*nres
 905           do j=1,3
 906             chomo(j,i,k)=c(j,i)
 907           enddo
 908         enddo
 909         do i=1,nres
 910           rescore(k,i)=0.2d0
 911           rescore2(k,i)=1.0d0
 912         enddo
 913       enddo
 914 c Read clusters
 915       do i=1,nclust
 916         read(ientin,*) ninclust(i),nresclust(i)
 917         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
 918         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
 919       enddo
 920 c
 921 c Loop over clusters
 922 c
 923       do l=1,nclust
 924         do ll = 1,ninclust(l)
 925
 926         k = inclust(ll,l)
 927         do i=1,nres
 928           idomain(k,i)=0
 929         enddo
 930         do i=1,nresclust(l)
 931           if (nnt.gt.1)  then
 932             idomain(k,iresclust(i,l)+1) = 1
 933           else
 934             idomain(k,iresclust(i,l)) = 1
 935           endif
 936         enddo
 937 c
 938 c     Distance restraints
 939 c
 940 c          ... --> odl(k,ii)
 941 C Copy the coordinates from reference coordinates (?)
 942         do i=1,2*nres
 943           do j=1,3
 944             c(j,i)=chomo(j,i,k)
 945 c           write (iout,*) "c(",j,i,") =",c(j,i)
 946           enddo
 947         enddo
 948         call int_from_cart1(.false.)
 949         call int_from_cart(.true.,.false.)
 950         call sc_loc_geom(.false.)
 951         do i=1,nres
 952           thetaref(i)=theta(i)
 953           phiref(i)=phi(i)
 954         enddo
 955         if (waga_dist.ne.0.0d0) then
 956           ii=0
 957           do i = nnt,nct-2
 958             do j=i+2,nct
 959
 960               x12=c(1,i)-c(1,j)
 961               y12=c(2,i)-c(2,j)
 962               z12=c(3,i)-c(3,j)
 963               distal=dsqrt(x12*x12+y12*y12+z12*z12)
 964 c              write (iout,*) k,i,j,distal,dist2_cut
 965
 966             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
 967      &            .and. distal.le.dist2_cut ) then
 968
 969               ii=ii+1
 970               ii_in_use(ii)=1
 971               l_homo(k,ii)=.true.
 972
 973 c             write (iout,*) "k",k
 974 c             write (iout,*) "i",i," j",j," constr_homology",
 975 c    &                       constr_homology
 976               ires_homo(ii)=i
 977               jres_homo(ii)=j
 978               odl(k,ii)=distal
 979               if (read2sigma) then
 980                 sigma_odl(k,ii)=0
 981                 do ik=i,j
 982                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
 983                 enddo
 984                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
 985                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
 986      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 987               else
 988                 if (odl(k,ii).le.dist_cut) then
 989                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
 990                 else
 991 #ifdef OLDSIGMA
 992                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 993      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
 994 #else
 995                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 996      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 997 #endif
 998                 endif
 999               endif
1000               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1001             else
1002               ii=ii+1
1003 c              l_homo(k,ii)=.false.
1004             endif
1005             enddo
1006           enddo
1007         lim_odl=ii
1008         endif
1009 c
1010 c     Theta, dihedral and SC retraints
1011 c
1012         if (waga_angle.gt.0.0d0) then
1013           do i = nnt+3,nct
1014             if (idomain(k,i).eq.0) then
1015 c               sigma_dih(k,i)=0.0
1016                cycle
1017             endif
1018             dih(k,i)=phiref(i)
1019             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1020      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1021 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
1022 c     &       " sigma_dihed",sigma_dih(k,i)
1023             if (sigma_dih(k,i).ne.0)
1024      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1025           enddo
1026           lim_dih=nct-nnt-2
1027         endif
1028
1029         if (waga_theta.gt.0.0d0) then
1030           do i = nnt+2,nct
1031              if (idomain(k,i).eq.0) then
1032 c              sigma_theta(k,i)=0.0
1033               cycle
1034              endif
1035              thetatpl(k,i)=thetaref(i)
1036              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1037      &                        rescore(k,i-2))/3.0
1038              if (sigma_theta(k,i).ne.0)
1039      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1040           enddo
1041         endif
1042
1043         if (waga_d.gt.0.0d0) then
1044           do i = nnt,nct
1045                if (itype(i).eq.10) cycle
1046                if (idomain(k,i).eq.0 ) then
1047 c                  sigma_d(k,i)=0.0
1048                   cycle
1049                endif
1050                xxtpl(k,i)=xxref(i)
1051                yytpl(k,i)=yyref(i)
1052                zztpl(k,i)=zzref(i)
1053                sigma_d(k,i)=rescore(k,i)
1054                if (sigma_d(k,i).ne.0)
1055      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1056                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
1057           enddo
1058         endif
1059       enddo ! l
1060       enddo ! ll
1061 c
1062 c remove distance restraints not used in any model from the list
1063 c shift data in all arrays
1064 c
1065       if (waga_dist.ne.0.0d0) then
1066         ii=0
1067         liiflag=.true.
1068         do i=nnt,nct-2
1069          do j=i+2,nct
1070           ii=ii+1
1071           if (ii_in_use(ii).eq.0.and.liiflag) then
1072             liiflag=.false.
1073             iistart=ii
1074           endif
1075           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1076      &                   .not.liiflag.and.ii.eq.lim_odl) then
1077              if (ii.eq.lim_odl) then
1078               iishift=ii-iistart+1
1079              else
1080               iishift=ii-iistart
1081              endif
1082              liiflag=.true.
1083              do ki=iistart,lim_odl-iishift
1084               ires_homo(ki)=ires_homo(ki+iishift)
1085               jres_homo(ki)=jres_homo(ki+iishift)
1086               ii_in_use(ki)=ii_in_use(ki+iishift)
1087               do k=1,constr_homology
1088                odl(k,ki)=odl(k,ki+iishift)
1089                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1090                l_homo(k,ki)=l_homo(k,ki+iishift)
1091               enddo
1092              enddo
1093              ii=ii-iishift
1094              lim_odl=lim_odl-iishift
1095           endif
1096          enddo
1097         enddo
1098       endif
1099
1100       return
1101    10 stop "Error infragment file"
1102       end
1103 c----------------------------------------------------------------------