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