source/wham/src/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 c     include 'include_unres/COMMON.VAR'
  13       include 'COMMON.INTERACT'
  14       include 'COMMON.LOCAL'
  15       include 'COMMON.NAMES'
  16       include 'COMMON.CHAIN'
  17       include 'COMMON.FFIELD'
  18       include 'COMMON.SBRIDGE'
  19       include 'COMMON.TORCNSTR'
  20       include 'COMMON.CONTROL'
  21       character*4 sequence(maxres)
  22       integer rescode
  23       double precision x(maxvar)
  24       character*320 controlcard,ucase
  25       dimension itype_pdb(maxres)
  26       logical seq_comp
  27       call card_concat(controlcard,.true.)
  28       call reada(controlcard,'SCAL14',scal14,0.4d0)
  29       call reada(controlcard,'SCALSCP',scalscp,1.0d0)
  30       call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
  31       call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
  32 C     Bartek
  33       call reada(controlcard,'WDFAD',wdfa_dist,0.0d0)
  34       call reada(controlcard,'WDFAT',wdfa_tor,0.0d0)
  35       call reada(controlcard,'WDFAN',wdfa_nei,0.0d0)
  36       call reada(controlcard,'WDFAB',wdfa_beta,0.0d0)
  37       write (iout,*) "wdfa_dist",wdfa_dist," wdfa_tor",wdfa_tor,
  38      &  " wdfa_nei",wdfa_nei," wdfa_beta",wdfa_beta
  39       r0_corr=cutoff_corr-delt_corr
  40       call readi(controlcard,"NRES",nres,0)
  41       iscode=index(controlcard,"ONE_LETTER")
  42       if (nres.le.0) then
  43         write (iout,*) "Error: no residues in molecule"
  44         return1
  45       endif
  46       if (nres.gt.maxres) then
  47         write (iout,*) "Error: too many residues",nres,maxres
  48       endif
  49       write(iout,*) 'nres=',nres
  50 C Read sequence of the protein
  51       if (iscode.gt.0) then
  52         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
  53       else
  54         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
  55       endif
  56 C Convert sequence to numeric code
  57       do i=1,nres
  58         itype(i)=rescode(i,sequence(i),iscode)
  59       enddo
  60       if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
  61         write (iout,*)
  62      &   "Glycine is the first full residue, initial dummy deleted"
  63         do i=1,nres
  64           itype(i)=itype(i+1)
  65         enddo
  66         nres=nres-1
  67       endif
  68       if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
  69         write (iout,*)
  70      &   "Glycine is the last full residue, terminal dummy deleted"
  71         nres=nres-1
  72       endif
  73       write (iout,*) "Numeric code:"
  74       write (iout,'(20i4)') (itype(i),i=1,nres)
  75       do i=1,nres-1
  76 #ifdef PROCOR
  77         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
  78 #else
  79         if (itype(i).eq.21) then
  80 #endif
  81           itel(i)=0
  82 #ifdef PROCOR
  83         else if (itype(i+1).ne.20) then
  84 #else
  85         else if (itype(i).ne.20) then
  86 #endif
  87           itel(i)=1
  88         else
  89           itel(i)=2
  90         endif
  91       enddo
  92       call read_bridge
  93
  94       if (with_dihed_constr) then
  95
  96       read (inp,*) ndih_constr
  97       if (ndih_constr.gt.0) then
  98         read (inp,*) ftors
  99         write (iout,*) 'FTORS',ftors
 100         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 101         write (iout,*)
 102      &   'There are',ndih_constr,' constraints on phi angles.'
 103         do i=1,ndih_constr
 104           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 105         enddo
 106         do i=1,ndih_constr
 107           phi0(i)=deg2rad*phi0(i)
 108           drange(i)=deg2rad*drange(i)
 109         enddo
 110       endif
 111
 112       endif
 113
 114       nnt=1
 115       nct=nres
 116       if (itype(1).eq.21) nnt=2
 117       if (itype(nres).eq.21) nct=nct-1
 118       write(iout,*) 'NNT=',NNT,' NCT=',NCT
 119
 120 C     Juyong:READ init_vars
 121 C     Initialize variables!
 122 C     Juyong:READ read_info
 123 C     READ fragment information!!
 124 C     both routines should be in dfa.F file!!
 125
 126       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
 127      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 128        write (iout,*) "Calling init_dfa_vars"
 129        call flush(iout)
 130        call init_dfa_vars
 131        write (iout,*) 'init_dfa_vars finished!'
 132        call flush(iout)
 133        call read_dfa_info
 134        write (iout,*) 'read_dfa_info finished!'
 135        call flush(iout)
 136       endif
 137
 138 c Read distance restraints
 139       if (constr_dist.gt.0) then
 140         if (refstr) call read_ref_structure(*11)
 141         call read_dist_constr
 142         call hpb_partition
 143       endif
 144
 145       if (constr_homology.gt.0) then
 146 c       write (iout,*) "About to call read_constr_homology"
 147 c       call flush(iout)
 148         call read_constr_homology
 149 c       write (iout,*) "Exit read_constr_homology"
 150 c       call flush(iout)
 151         if (indpdb.gt.0 .or. pdbref) then
 152           do i=1,2*nres
 153             do j=1,3
 154               c(j,i)=crefjlee(j,i)
 155               cref(j,i)=crefjlee(j,i)
 156             enddo
 157           enddo
 158         endif
 159 #ifdef DEBUG
 160         write (iout,*) "Array C"
 161         do i=1,nres
 162           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
 163      &      (c(j,i+nres),j=1,3)
 164         enddo
 165         write (iout,*) "Array Cref"
 166         do i=1,nres
 167           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
 168      &      (cref(j,i+nres),j=1,3)
 169         enddo
 170 #endif
 171 #ifdef DEBUG
 172        call int_from_cart1(.false.)
 173        call sc_loc_geom(.false.)
 174        do i=1,nres
 175          thetaref(i)=theta(i)
 176          phiref(i)=phi(i)
 177          write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i)
 178        enddo
 179        do i=1,nres-1
 180          do j=1,3
 181            dc(j,i)=c(j,i+1)-c(j,i)
 182            dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 183          enddo
 184        enddo
 185        do i=2,nres-1
 186          do j=1,3
 187            dc(j,i+nres)=c(j,i+nres)-c(j,i)
 188            dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 189          enddo
 190        enddo
 191 #endif
 192       else
 193         homol_nset=0
 194       endif
 195
 196
 197       call setup_var
 198       call init_int_table
 199       if (ns.gt.0) then
 200         write (iout,'(/a,i3,a)') 'The chain contains',ns,
 201      &  ' disulfide-bridging cysteines.'
 202         write (iout,'(20i4)') (iss(i),i=1,ns)
 203         write (iout,'(/a/)') 'Pre-formed links are:'
 204         do i=1,nss
 205           i1=ihpb(i)-nres
 206           i2=jhpb(i)-nres
 207           it1=itype(i1)
 208           it2=itype(i2)
 209          write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 210      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
 211      &    dhpb(i),ebr,forcon(i)
 212         enddo
 213       endif
 214       write (iout,'(a)')
 215       return
 216    11 stop "Error reading reference structure"
 217       end
 218 c-----------------------------------------------------------------------------
 219       logical function seq_comp(itypea,itypeb,length)
 220       implicit none
 221       integer length,itypea(length),itypeb(length)
 222       integer i
 223       do i=1,length
 224         if (itypea(i).ne.itypeb(i)) then
 225           seq_comp=.false.
 226           return
 227         endif
 228       enddo
 229       seq_comp=.true.
 230       return
 231       end
 232 c-----------------------------------------------------------------------------
 233       subroutine read_bridge
 234 C Read information about disulfide bridges.
 235       implicit real*8 (a-h,o-z)
 236       include 'DIMENSIONS'
 237       include 'DIMENSIONS.ZSCOPT'
 238       include 'COMMON.IOUNITS'
 239       include 'COMMON.GEO'
 240       include 'COMMON.VAR'
 241       include 'COMMON.INTERACT'
 242       include 'COMMON.NAMES'
 243       include 'COMMON.CHAIN'
 244       include 'COMMON.FFIELD'
 245       include 'COMMON.SBRIDGE'
 246 C Read bridging residues.
 247       read (inp,*) ns,(iss(i),i=1,ns)
 248       print *,'ns=',ns
 249       write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 250 C Check whether the specified bridging residues are cystines.
 251       do i=1,ns
 252         if (itype(iss(i)).ne.1) then
 253           write (iout,'(2a,i3,a)')
 254      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 255      &   ' can form a disulfide bridge?!!!'
 256           write (*,'(2a,i3,a)')
 257      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 258      &   ' can form a disulfide bridge?!!!'
 259          stop
 260         endif
 261       enddo
 262 C Read preformed bridges.
 263       if (ns.gt.0) then
 264       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 265       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 266       if (nss.gt.0) then
 267         nhpb=nss
 268 C Check if the residues involved in bridges are in the specified list of
 269 C bridging residues.
 270         do i=1,nss
 271           do j=1,i-1
 272             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 273      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 274               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 275      &      ' contains residues present in other pairs.'
 276               write (*,'(a,i3,a)') 'Disulfide pair',i,
 277      &      ' contains residues present in other pairs.'
 278               stop
 279             endif
 280           enddo
 281           do j=1,ns
 282             if (ihpb(i).eq.iss(j)) goto 10
 283           enddo
 284           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 285    10     continue
 286           do j=1,ns
 287             if (jhpb(i).eq.iss(j)) goto 20
 288           enddo
 289           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 290    20     continue
 291           dhpb(i)=dbr
 292           forcon(i)=fbr
 293         enddo
 294         do i=1,nss
 295           ihpb(i)=ihpb(i)+nres
 296           jhpb(i)=jhpb(i)+nres
 297         enddo
 298       endif
 299       endif
 300       if (ns.gt.0.and.dyn_ss) then
 301 C /06/28/2013 Adasko:ns is number of Cysteins bonded also called half of
 302 C the bond
 303           do i=nss+1,nhpb
 304 C /06/28/2013 Adasko: nss number of full SS bonds
 305             ihpb(i-nss)=ihpb(i)
 306             jhpb(i-nss)=jhpb(i)
 307             forcon(i-nss)=forcon(i)
 308             dhpb(i-nss)=dhpb(i)
 309           enddo
 310           nhpb=nhpb-nss
 311           nss=0
 312           call hpb_partition
 313           do i=1,ns
 314             dyn_ss_mask(iss(i))=.true.
 315 C /06/28/2013 Adasko: dyn_ss_mask which Cysteins can form disulfidebond
 316 c          write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
 317           enddo
 318       endif
 319       return
 320       end
 321 c------------------------------------------------------------------------------
 322       subroutine read_angles(kanal,iscor,energ,iprot,*)
 323       implicit real*8 (a-h,o-z)
 324       include 'DIMENSIONS'
 325       include 'DIMENSIONS.ZSCOPT'
 326       include 'COMMON.INTERACT'
 327       include 'COMMON.SBRIDGE'
 328       include 'COMMON.GEO'
 329       include 'COMMON.VAR'
 330       include 'COMMON.CHAIN'
 331       include 'COMMON.IOUNITS'
 332       character*80 lineh
 333       read(kanal,'(a80)',end=10,err=10) lineh
 334       read(lineh(:5),*,err=8) ic
 335       read(lineh(6:),*,err=8) energ
 336       goto 9
 337     8 ic=1
 338       print *,'error, assuming e=1d10',lineh
 339       energ=1d10
 340       nss=0
 341     9 continue
 342       read(lineh(18:),*,end=10,err=10) nss
 343       IF (NSS.LT.9) THEN
 344         read (lineh(20:),*,end=10,err=10)
 345      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
 346       ELSE
 347         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
 348         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
 349      &    I=9,NSS),iscor
 350       ENDIF
 351 c      print *,"energy",energ," iscor",iscor
 352       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 353       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 354       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 355       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 356       do i=1,nres
 357         theta(i)=deg2rad*theta(i)
 358         phi(i)=deg2rad*phi(i)
 359         alph(i)=deg2rad*alph(i)
 360         omeg(i)=deg2rad*omeg(i)
 361       enddo
 362       return
 363    10 return1
 364       end
 365 c-------------------------------------------------------------------------------
 366       subroutine read_dist_constr
 367       implicit real*8 (a-h,o-z)
 368       include 'DIMENSIONS'
 369       include 'DIMENSIONS.ZSCOPT'
 370       include 'DIMENSIONS.FREE'
 371       include 'COMMON.CONTROL'
 372       include 'COMMON.CHAIN'
 373       include 'COMMON.IOUNITS'
 374       include 'COMMON.SBRIDGE'
 375       integer ifrag_(2,100),ipair_(2,100)
 376       double precision wfrag_(100),wpair_(100)
 377       character*500 controlcard
 378 c      write (iout,*) "Calling read_dist_constr"
 379 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 380 c      call flush(iout)
 381       call card_concat(controlcard,.true.)
 382       call readi(controlcard,"NFRAG",nfrag_,0)
 383       call readi(controlcard,"NPAIR",npair_,0)
 384       call readi(controlcard,"NDIST",ndist_,0)
 385       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 386       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 387       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 388       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 389       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 390       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 391       write (iout,*) "IFRAG"
 392       do i=1,nfrag_
 393         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 394       enddo
 395       write (iout,*) "IPAIR"
 396       do i=1,npair_
 397         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 398       enddo
 399       call flush(iout)
 400       if (.not.refstr .and. nfrag_.gt.0) then
 401         write (iout,*)
 402      &  "ERROR: no reference structure to compute distance restraints"
 403         write (iout,*)
 404      &  "Restraints must be specified explicitly (NDIST=number)"
 405         stop
 406       endif
 407       if (nfrag_.lt.2 .and. npair_.gt.0) then
 408         write (iout,*) "ERROR: Less than 2 fragments specified",
 409      &   " but distance restraints between pairs requested"
 410         stop
 411       endif
 412       call flush(iout)
 413       do i=1,nfrag_
 414         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 415         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 416      &    ifrag_(2,i)=nstart_sup+nsup-1
 417 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 418         call flush(iout)
 419         if (wfrag_(i).gt.0.0d0) then
 420         do j=ifrag_(1,i),ifrag_(2,i)-1
 421           do k=j+1,ifrag_(2,i)
 422             write (iout,*) "j",j," k",k
 423             ddjk=dist(j,k)
 424             if (constr_dist.eq.1) then
 425             nhpb=nhpb+1
 426             ihpb(nhpb)=j
 427             jhpb(nhpb)=k
 428               dhpb(nhpb)=ddjk
 429             forcon(nhpb)=wfrag_(i)
 430             else if (constr_dist.eq.2) then
 431               if (ddjk.le.dist_cut) then
 432                 nhpb=nhpb+1
 433                 ihpb(nhpb)=j
 434                 jhpb(nhpb)=k
 435                 dhpb(nhpb)=ddjk
 436                 forcon(nhpb)=wfrag_(i)
 437               endif
 438             else
 439               nhpb=nhpb+1
 440               ihpb(nhpb)=j
 441               jhpb(nhpb)=k
 442               dhpb(nhpb)=ddjk
 443               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 444             endif
 445             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
 446      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 447           enddo
 448         enddo
 449         endif
 450       enddo
 451       do i=1,npair_
 452         if (wpair_(i).gt.0.0d0) then
 453         ii = ipair_(1,i)
 454         jj = ipair_(2,i)
 455         if (ii.gt.jj) then
 456           itemp=ii
 457           ii=jj
 458           jj=itemp
 459         endif
 460         do j=ifrag_(1,ii),ifrag_(2,ii)
 461           do k=ifrag_(1,jj),ifrag_(2,jj)
 462             nhpb=nhpb+1
 463             ihpb(nhpb)=j
 464             jhpb(nhpb)=k
 465             forcon(nhpb)=wpair_(i)
 466             dhpb(nhpb)=dist(j,k)
 467             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 468      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 469           enddo
 470         enddo
 471         endif
 472       enddo
 473       do i=1,ndist_
 474         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
 475      &     ibecarb(i),forcon(nhpb+1)
 476         if (forcon(nhpb+1).gt.0.0d0) then
 477           nhpb=nhpb+1
 478           if (ibecarb(i).gt.0) then
 479             ihpb(i)=ihpb(i)+nres
 480             jhpb(i)=jhpb(i)+nres
 481           endif
 482           if (dhpb(nhpb).eq.0.0d0)
 483      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 484         endif
 485       enddo
 486       do i=1,nhpb
 487           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
 488      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
 489       enddo
 490       call flush(iout)
 491       return
 492       end
 493
 494
 495
 496 c====-------------------------------------------------------------------
 497       subroutine read_constr_homology
 498
 499       include 'DIMENSIONS'
 500       include 'DIMENSIONS.ZSCOPT'
 501       include 'DIMENSIONS.FREE'
 502 #ifdef MPI
 503       include 'mpif.h'
 504 #endif
 505       include 'COMMON.SETUP'
 506       include 'COMMON.CONTROL'
 507       include 'COMMON.CHAIN'
 508       include 'COMMON.IOUNITS'
 509       include 'COMMON.GEO'
 510       include 'COMMON.INTERACT'
 511       include 'COMMON.NAMES'
 512       include 'COMMON.HOMRESTR'
 513 c
 514 c For new homol impl
 515 c
 516       include 'COMMON.VAR'
 517 c     include 'include_unres/COMMON.VAR'
 518 c
 519
 520 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
 521 c    &                 dist_cut
 522 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
 523 c    &    sigma_odl_temp(maxres,maxres,max_template)
 524       character*2 kic2
 525       character*24 model_ki_dist, model_ki_angle
 526       character*500 controlcard
 527       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
 528       integer idomain(max_template,maxres)
 529       logical lprn /.true./
 530       integer ilen
 531       external ilen
 532       logical unres_pdb,liiflag
 533 c
 534 c     FP - Nov. 2014 Temporary specifications for new vars
 535 c
 536       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
 537       double precision, dimension (max_template,maxres) :: rescore
 538       double precision, dimension (max_template,maxres) :: rescore2
 539       character*24 tpl_k_rescore
 540 c -----------------------------------------------------------------
 541 c Reading multiple PDB ref structures and calculation of retraints
 542 c not using pre-computed ones stored in files model_ki_{dist,angle}
 543 c FP (Nov., 2014)
 544 c -----------------------------------------------------------------
 545 c
 546 c
 547 c Alternative: reading from input
 548       call card_concat(controlcard,.true.)
 549       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
 550       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
 551       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
 552       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
 553       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
 554       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
 555       call readi(controlcard,"HOMOL_SET",homol_nset,1)
 556       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
 557       call readi(controlcard,"IHSET",ihset,1)
 558       if (homol_nset.gt.1)then
 559          call card_concat(controlcard,.true.)
 560          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
 561          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 562           write(iout,*) "iset homology_weight "
 563 c         do i=1,homol_nset
 564 c          write(iout,*) i,waga_homology(i)
 565 c         enddo
 566          endif
 567          iset=mod(kolor,homol_nset)+1
 568       else
 569        iset=1
 570        waga_homology(1)=1.0
 571       endif
 572 c     write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
 573
 574 cd      write (iout,*) "nnt",nnt," nct",nct
 575 cd      call flush(iout)
 576
 577
 578       lim_odl=0
 579       lim_dih=0
 580 c
 581 c  New
 582 c
 583       lim_theta=0
 584       lim_xx=0
 585 c
 586 c  Reading HM global scores (prob not required)
 587 c
 588       do i = nnt,nct
 589         do k=1,constr_homology
 590          idomain(k,i)=0
 591         enddo
 592       enddo
 593 c     open (4,file="HMscore")
 594 c     do k=1,constr_homology
 595 c       read (4,*,end=521) hmscore_tmp
 596 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
 597 c       write(*,*) "Model", k, ":", hmscore(k)
 598 c     enddo
 599 c521  continue
 600
 601       ii=0
 602       do i = nnt,nct-2
 603         do j=i+2,nct
 604         ii=ii+1
 605         ii_in_use(ii)=0
 606         enddo
 607       enddo
 608 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 609
 610       do k=1,constr_homology
 611
 612         read(inp,'(a)') pdbfile
 613 c  Next stament causes error upon compilation (?)
 614 c       if(me.eq.king.or. .not. out1file)
 615 c         write (iout,'(2a)') 'PDB data will be read from file ',
 616 c    &   pdbfile(:ilen(pdbfile))
 617          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
 618      &  pdbfile(:ilen(pdbfile))
 619         open(ipdbin,file=pdbfile,status='old',err=33)
 620         goto 34
 621   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
 622      &  pdbfile(:ilen(pdbfile))
 623         stop
 624   34    continue
 625 c        print *,'Begin reading pdb data'
 626 c
 627 c Files containing res sim or local scores (former containing sigmas)
 628 c
 629
 630         write(kic2,'(bz,i2.2)') k
 631
 632         tpl_k_rescore="template"//kic2//".sco"
 633
 634         unres_pdb=.false.
 635         call readpdb
 636         do i=1,2*nres
 637           do j=1,3
 638             crefjlee(j,i)=c(j,i)
 639           enddo
 640         enddo
 641 #ifdef DEBUG
 642         do i=1,nres
 643           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
 644      &      (crefjlee(j,i+nres),j=1,3)
 645         enddo
 646 #endif
 647         write (iout,*) "read_constr_homology: after reading pdb file"
 648         call flush(iout)
 649
 650 c
 651 c     Distance restraints
 652 c
 653 c          ... --> odl(k,ii)
 654 C Copy the coordinates from reference coordinates (?)
 655         do i=1,2*nres
 656           do j=1,3
 657             c(j,i)=cref(j,i)
 658 c           write (iout,*) "c(",j,i,") =",c(j,i)
 659           enddo
 660         enddo
 661 c
 662 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
 663 c
 664 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
 665           open (ientin,file=tpl_k_rescore,status='old')
 666           if (nnt.gt.1) rescore(k,1)=0.0d0
 667           do irec=nnt,nct ! loop for reading res sim
 668             if (read2sigma) then
 669              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
 670      &                                idomain_tmp
 671              i_tmp=i_tmp+nnt-1
 672              idomain(k,i_tmp)=idomain_tmp
 673              rescore(k,i_tmp)=rescore_tmp
 674              rescore2(k,i_tmp)=rescore2_tmp
 675              write(iout,'(a7,i5,2f10.5,i5)') "rescore",
 676      &                      i_tmp,rescore2_tmp,rescore_tmp,
 677      &                                idomain_tmp
 678             else
 679              idomain(k,irec)=1
 680              read (ientin,*,end=1401) rescore_tmp
 681
 682 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
 683              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
 684 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
 685             endif
 686           enddo
 687  1401   continue
 688         close (ientin)
 689         if (waga_dist.ne.0.0d0) then
 690           ii=0
 691           do i = nnt,nct-2
 692             do j=i+2,nct
 693
 694               x12=c(1,i)-c(1,j)
 695               y12=c(2,i)-c(2,j)
 696               z12=c(3,i)-c(3,j)
 697               distal=dsqrt(x12*x12+y12*y12+z12*z12)
 698 c              write (iout,*) k,i,j,distal,dist2_cut
 699
 700             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
 701      &            .and. distal.le.dist2_cut ) then
 702
 703               ii=ii+1
 704               ii_in_use(ii)=1
 705               l_homo(k,ii)=.true.
 706
 707 c             write (iout,*) "k",k
 708 c             write (iout,*) "i",i," j",j," constr_homology",
 709 c    &                       constr_homology
 710               ires_homo(ii)=i
 711               jres_homo(ii)=j
 712               odl(k,ii)=distal
 713               if (read2sigma) then
 714                 sigma_odl(k,ii)=0
 715                 do ik=i,j
 716                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
 717                 enddo
 718                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
 719                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
 720      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 721               else
 722                 if (odl(k,ii).le.dist_cut) then
 723                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
 724                 else
 725 #ifdef OLDSIGMA
 726                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 727      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
 728 #else
 729                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 730      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 731 #endif
 732                 endif
 733               endif
 734               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
 735             else
 736               ii=ii+1
 737               l_homo(k,ii)=.false.
 738             endif
 739             enddo
 740           enddo
 741         lim_odl=ii
 742         endif
 743 c
 744 c     Theta, dihedral and SC retraints
 745 c
 746         if (waga_angle.gt.0.0d0) then
 747 c         open (ientin,file=tpl_k_sigma_dih,status='old')
 748 c         do irec=1,maxres-3 ! loop for reading sigma_dih
 749 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
 750 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
 751 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
 752 c    &                            sigma_dih(k,i+nnt-1)
 753 c         enddo
 754 c1402   continue
 755 c         close (ientin)
 756           do i = nnt+3,nct
 757             if (idomain(k,i).eq.0) then
 758                sigma_dih(k,i)=0.0
 759                cycle
 760             endif
 761             dih(k,i)=phiref(i) ! right?
 762 c           read (ientin,*) sigma_dih(k,i) ! original variant
 763 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
 764 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
 765 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 766 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
 767 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
 768
 769             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
 770      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
 771 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
 772 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
 773 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
 774 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 775 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
 776 c   Instead of res sim other local measure of b/b str reliability possible
 777             if (sigma_dih(k,i).ne.0)
 778      &      sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
 779 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
 780           enddo
 781           lim_dih=nct-nnt-2
 782         endif
 783
 784         if (waga_theta.gt.0.0d0) then
 785 c         open (ientin,file=tpl_k_sigma_theta,status='old')
 786 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
 787 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
 788 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
 789 c    &                              sigma_theta(k,i+nnt-1)
 790 c         enddo
 791 c1403   continue
 792 c         close (ientin)
 793
 794           do i = nnt+2,nct ! right? without parallel.
 795 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
 796 c         do i=ithet_start,ithet_end ! with FG parallel.
 797              if (idomain(k,i).eq.0) then
 798               sigma_theta(k,i)=0.0
 799               cycle
 800              endif
 801              thetatpl(k,i)=thetaref(i)
 802 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
 803 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
 804 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 805 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
 806 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
 807              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
 808      &                        rescore(k,i-2))/3.0
 809 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
 810              if (sigma_theta(k,i).ne.0)
 811      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 812
 813 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 814 c                             rescore(k,i-2) !  right expression ?
 815 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
 816           enddo
 817         endif
 818
 819         if (waga_d.gt.0.0d0) then
 820 c       open (ientin,file=tpl_k_sigma_d,status='old')
 821 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
 822 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
 823 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
 824 c    &                          sigma_d(k,i+nnt-1)
 825 c         enddo
 826 c1404   continue
 827
 828           do i = nnt,nct ! right? without parallel.
 829 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
 830 c         do i=loc_start,loc_end ! with FG parallel.
 831                if (itype(i).eq.10) cycle
 832                if (idomain(k,i).eq.0 ) then
 833                   sigma_d(k,i)=0.0
 834                   cycle
 835                endif
 836                xxtpl(k,i)=xxref(i)
 837                yytpl(k,i)=yyref(i)
 838                zztpl(k,i)=zzref(i)
 839 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
 840 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
 841 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
 842 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
 843                sigma_d(k,i)=rescore(k,i) !  right expression ?
 844                if (sigma_d(k,i).ne.0)
 845      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
 846
 847 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
 848 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
 849 c              read (ientin,*) sigma_d(k,i) ! 1st variant
 850           enddo
 851         endif
 852       enddo
 853 c
 854 c remove distance restraints not used in any model from the list
 855 c shift data in all arrays
 856 c
 857       if (waga_dist.ne.0.0d0) then
 858         ii=0
 859         liiflag=.true.
 860         do i=nnt,nct-2
 861          do j=i+2,nct
 862           ii=ii+1
 863           if (ii_in_use(ii).eq.0.and.liiflag) then
 864             liiflag=.false.
 865             iistart=ii
 866           endif
 867           if (ii_in_use(ii).ne.0.and..not.liiflag) then
 868              iishift=ii-iistart
 869              liiflag=.true.
 870              do ki=iistart,lim_odl-iishift
 871               ires_homo(ki)=ires_homo(ki+iishift)
 872               jres_homo(ki)=jres_homo(ki+iishift)
 873               ii_in_use(ki)=ii_in_use(ki+iishift)
 874               do k=1,constr_homology
 875                odl(k,ki)=odl(k,ki+iishift)
 876                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
 877                l_homo(k,ki)=l_homo(k,ki+iishift)
 878               enddo
 879              enddo
 880              ii=ii-iishift
 881              lim_odl=lim_odl-iishift
 882           endif
 883          enddo
 884         enddo
 885       endif
 886       if (constr_homology.gt.0) call homology_partition
 887       if (constr_homology.gt.0) call init_int_table
 888 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
 889 cd     & "lim_xx=",lim_xx
 890 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
 891 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
 892 c
 893 c Print restraints
 894 c
 895       if (.not.lprn) return
 896 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 897       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 898        write (iout,*) "Distance restraints from templates"
 899        do ii=1,lim_odl
 900        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
 901      &  ii,ires_homo(ii),jres_homo(ii),
 902      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
 903      &  ki=1,constr_homology)
 904        enddo
 905        write (iout,*) "Dihedral angle restraints from templates"
 906        do i=nnt+3,nct
 907         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
 908      &      (rad2deg*dih(ki,i),
 909      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
 910        enddo
 911        write (iout,*) "Virtual-bond angle restraints from templates"
 912        do i=nnt+2,nct
 913         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
 914      &      (rad2deg*thetatpl(ki,i),
 915      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
 916        enddo
 917        write (iout,*) "SC restraints from templates"
 918        do i=nnt,nct
 919         write(iout,'(i5,100(4f8.2,4x))') i,
 920      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
 921      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
 922        enddo
 923       endif
 924 c -----------------------------------------------------------------
 925       return
 926       end
 927 c----------------------------------------------------------------------