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.HOMRESTR'
 512 c
 513 c For new homol impl
 514 c
 515       include 'COMMON.VAR'
 516 c     include 'include_unres/COMMON.VAR'
 517 c
 518
 519 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
 520 c    &                 dist_cut
 521 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
 522 c    &    sigma_odl_temp(maxres,maxres,max_template)
 523       character*2 kic2
 524       character*24 model_ki_dist, model_ki_angle
 525       character*500 controlcard
 526       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
 527       integer idomain(max_template,maxres)
 528       logical lprn /.true./
 529       integer ilen
 530       external ilen
 531       logical unres_pdb
 532 c
 533 c     FP - Nov. 2014 Temporary specifications for new vars
 534 c
 535       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
 536       double precision, dimension (max_template,maxres) :: rescore
 537       double precision, dimension (max_template,maxres) :: rescore2
 538       character*24 tpl_k_rescore
 539 c -----------------------------------------------------------------
 540 c Reading multiple PDB ref structures and calculation of retraints
 541 c not using pre-computed ones stored in files model_ki_{dist,angle}
 542 c FP (Nov., 2014)
 543 c -----------------------------------------------------------------
 544 c
 545 c
 546 c Alternative: reading from input
 547       call card_concat(controlcard,.true.)
 548       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
 549       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
 550       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
 551       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
 552       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
 553       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
 554       call readi(controlcard,"HOMOL_SET",homol_nset,1)
 555       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
 556       call readi(controlcard,"IHSET",ihset,1)
 557       if (homol_nset.gt.1)then
 558          call card_concat(controlcard,.true.)
 559          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
 560          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 561           write(iout,*) "iset homology_weight "
 562 c         do i=1,homol_nset
 563 c          write(iout,*) i,waga_homology(i)
 564 c         enddo
 565          endif
 566          iset=mod(kolor,homol_nset)+1
 567       else
 568        iset=1
 569        waga_homology(1)=1.0
 570       endif
 571 c     write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
 572
 573 cd      write (iout,*) "nnt",nnt," nct",nct
 574 cd      call flush(iout)
 575
 576
 577       lim_odl=0
 578       lim_dih=0
 579 c
 580 c  New
 581 c
 582       lim_theta=0
 583       lim_xx=0
 584 c
 585 c  Reading HM global scores (prob not required)
 586 c
 587       do i = nnt,nct
 588         do k=1,constr_homology
 589          idomain(k,i)=0
 590         enddo
 591       enddo
 592 c     open (4,file="HMscore")
 593 c     do k=1,constr_homology
 594 c       read (4,*,end=521) hmscore_tmp
 595 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
 596 c       write(*,*) "Model", k, ":", hmscore(k)
 597 c     enddo
 598 c521  continue
 599
 600       ii=0
 601       do i = nnt,nct-2
 602         do j=i+2,nct
 603         ii=ii+1
 604         ii_in_use(ii)=0
 605         enddo
 606       enddo
 607 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 608
 609       do k=1,constr_homology
 610
 611         read(inp,'(a)') pdbfile
 612 c  Next stament causes error upon compilation (?)
 613 c       if(me.eq.king.or. .not. out1file)
 614 c         write (iout,'(2a)') 'PDB data will be read from file ',
 615 c    &   pdbfile(:ilen(pdbfile))
 616          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
 617      &  pdbfile(:ilen(pdbfile))
 618         open(ipdbin,file=pdbfile,status='old',err=33)
 619         goto 34
 620   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
 621      &  pdbfile(:ilen(pdbfile))
 622         stop
 623   34    continue
 624 c        print *,'Begin reading pdb data'
 625 c
 626 c Files containing res sim or local scores (former containing sigmas)
 627 c
 628
 629         write(kic2,'(bz,i2.2)') k
 630
 631         tpl_k_rescore="template"//kic2//".sco"
 632
 633         unres_pdb=.false.
 634         call readpdb
 635         do i=1,2*nres
 636           do j=1,3
 637             crefjlee(j,i)=c(j,i)
 638           enddo
 639         enddo
 640 #ifdef DEBUG
 641         do i=1,nres
 642           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
 643      &      (crefjlee(j,i+nres),j=1,3)
 644         enddo
 645 #endif
 646         write (iout,*) "read_constr_homology: after reading pdb file"
 647         call flush(iout)
 648
 649 c
 650 c     Distance restraints
 651 c
 652 c          ... --> odl(k,ii)
 653 C Copy the coordinates from reference coordinates (?)
 654         do i=1,2*nres
 655           do j=1,3
 656             c(j,i)=cref(j,i)
 657 c           write (iout,*) "c(",j,i,") =",c(j,i)
 658           enddo
 659         enddo
 660 c
 661 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
 662 c
 663 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
 664           open (ientin,file=tpl_k_rescore,status='old')
 665           if (nnt.gt.1) rescore(k,1)=0.0d0
 666           do irec=nnt,maxdim ! loop for reading res sim
 667             if (read2sigma) then
 668              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
 669      &                                idomain_tmp
 670              i_tmp=i_tmp+nnt-1
 671              idomain(k,i_tmp)=idomain_tmp
 672              rescore(k,i_tmp)=rescore_tmp
 673              rescore2(k,i_tmp)=rescore2_tmp
 674             else
 675              idomain(k,irec)=1
 676              read (ientin,*,end=1401) rescore_tmp
 677
 678 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
 679              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
 680 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
 681             endif
 682           enddo
 683  1401   continue
 684         close (ientin)
 685         if (waga_dist.ne.0.0d0) then
 686           ii=0
 687           do i = nnt,nct-2
 688             do j=i+2,nct
 689
 690               x12=c(1,i)-c(1,j)
 691               y12=c(2,i)-c(2,j)
 692               z12=c(3,i)-c(3,j)
 693               distal=dsqrt(x12*x12+y12*y12+z12*z12)
 694 c              write (iout,*) k,i,j,distal,dist2_cut
 695
 696             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
 697      &            .and. distal.le.dist2_cut ) then
 698
 699               ii=ii+1
 700               ii_in_use(ii)=1
 701               l_homo(k,ii)=.true.
 702
 703 c             write (iout,*) "k",k
 704 c             write (iout,*) "i",i," j",j," constr_homology",
 705 c    &                       constr_homology
 706               ires_homo(ii)=i
 707               jres_homo(ii)=j
 708               odl(k,ii)=distal
 709               if (read2sigma) then
 710                 sigma_odl(k,ii)=0
 711                 do ik=i,j
 712                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
 713                 enddo
 714                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
 715                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
 716      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 717               else
 718                 if (odl(k,ii).le.dist_cut) then
 719                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
 720                 else
 721 #ifdef OLDSIGMA
 722                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 723      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
 724 #else
 725                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 726      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 727 #endif
 728                 endif
 729               endif
 730               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
 731             else
 732               ii=ii+1
 733               l_homo(k,ii)=.false.
 734             endif
 735             enddo
 736           enddo
 737         lim_odl=ii
 738         endif
 739 c
 740 c     Theta, dihedral and SC retraints
 741 c
 742         if (waga_angle.gt.0.0d0) then
 743 c         open (ientin,file=tpl_k_sigma_dih,status='old')
 744 c         do irec=1,maxres-3 ! loop for reading sigma_dih
 745 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
 746 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
 747 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
 748 c    &                            sigma_dih(k,i+nnt-1)
 749 c         enddo
 750 c1402   continue
 751 c         close (ientin)
 752           do i = nnt+3,nct
 753             if (idomain(k,i).eq.0) then
 754                sigma_dih(k,i)=0.0
 755                cycle
 756             endif
 757             dih(k,i)=phiref(i) ! right?
 758 c           read (ientin,*) sigma_dih(k,i) ! original variant
 759 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
 760 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
 761 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 762 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
 763 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
 764
 765             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
 766      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
 767 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
 768 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
 769 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
 770 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 771 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
 772 c   Instead of res sim other local measure of b/b str reliability possible
 773             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
 774 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
 775           enddo
 776           lim_dih=nct-nnt-2
 777         endif
 778
 779         if (waga_theta.gt.0.0d0) then
 780 c         open (ientin,file=tpl_k_sigma_theta,status='old')
 781 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
 782 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
 783 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
 784 c    &                              sigma_theta(k,i+nnt-1)
 785 c         enddo
 786 c1403   continue
 787 c         close (ientin)
 788
 789           do i = nnt+2,nct ! right? without parallel.
 790 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
 791 c         do i=ithet_start,ithet_end ! with FG parallel.
 792              if (idomain(k,i).eq.0) then
 793               sigma_theta(k,i)=0.0
 794               cycle
 795              endif
 796              thetatpl(k,i)=thetaref(i)
 797 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
 798 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
 799 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 800 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
 801 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
 802              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
 803      &                        rescore(k,i-2))/3.0
 804 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
 805              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 806
 807 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 808 c                             rescore(k,i-2) !  right expression ?
 809 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
 810           enddo
 811         endif
 812         lim_theta=nct-nnt-1
 813
 814         if (waga_d.gt.0.0d0) then
 815 c       open (ientin,file=tpl_k_sigma_d,status='old')
 816 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
 817 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
 818 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
 819 c    &                          sigma_d(k,i+nnt-1)
 820 c         enddo
 821 c1404   continue
 822
 823           do i = nnt,nct ! right? without parallel.
 824 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
 825 c         do i=loc_start,loc_end ! with FG parallel.
 826                if (itype(i).eq.10) cycle
 827                if (idomain(k,i).eq.0 ) then
 828                   sigma_d(k,i)=0.0
 829                   cycle
 830                endif
 831                xxtpl(k,i)=xxref(i)
 832                yytpl(k,i)=yyref(i)
 833                zztpl(k,i)=zzref(i)
 834 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
 835 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
 836 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
 837 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
 838                sigma_d(k,i)=rescore(k,i) !  right expression ?
 839                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
 840
 841 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
 842 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
 843 c              read (ientin,*) sigma_d(k,i) ! 1st variant
 844                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
 845           enddo
 846           lim_xx=nct-nnt+1
 847         endif
 848       enddo
 849 c
 850 c remove distance restraints not used in any model from the list
 851 c shift data in all arrays
 852 c
 853       if (waga_dist.ne.0.0d0) then
 854         ii=0
 855         do i=nnt,nct-2
 856          do j=i+2,nct
 857           ii=ii+1
 858           if (ii_in_use(ii).eq.0) then
 859              do ki=ii,lim_odl-1
 860               ires_homo(ki)=ires_homo(ki+1)
 861               jres_homo(ki)=jres_homo(ki+1)
 862               ii_in_use(ki)=ii_in_use(ki+1)
 863               do k=1,constr_homology
 864                odl(k,ki)=odl(k,ki+1)
 865                sigma_odl(k,ki)=sigma_odl(k,ki+1)
 866                l_homo(k,ki)=l_homo(k,ki+1)
 867               enddo
 868              enddo
 869              ii=ii-1
 870              lim_odl=lim_odl-1
 871           endif
 872          enddo
 873         enddo
 874       endif
 875       if (constr_homology.gt.0) call homology_partition
 876       if (constr_homology.gt.0) call init_int_table
 877 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
 878 cd     & "lim_xx=",lim_xx
 879 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
 880 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
 881 c
 882 c Print restraints
 883 c
 884       if (.not.lprn) return
 885 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 886       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 887        write (iout,*) "Distance restraints from templates"
 888        do ii=1,lim_odl
 889        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
 890      &  ii,ires_homo(ii),jres_homo(ii),
 891      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
 892      &  ki=1,constr_homology)
 893        enddo
 894        write (iout,*) "Dihedral angle restraints from templates"
 895        do i=nnt+3,lim_dih
 896         write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*dih(ki,i),
 897      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
 898        enddo
 899        write (iout,*) "Virtual-bond angle restraints from templates"
 900        do i=nnt+2,lim_theta
 901         write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
 902      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
 903        enddo
 904        write (iout,*) "SC restraints from templates"
 905        do i=nnt,lim_xx
 906         write(iout,'(i5,100(4f8.2,4x))') i,
 907      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
 908      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
 909        enddo
 910       endif
 911 c -----------------------------------------------------------------
 912       return
 913       end
 914 c----------------------------------------------------------------------