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 'COMMON.CONTROL'
 370       include 'COMMON.CHAIN'
 371       include 'COMMON.IOUNITS'
 372       include 'COMMON.SBRIDGE'
 373       integer ifrag_(2,100),ipair_(2,100)
 374       double precision wfrag_(100),wpair_(100)
 375       character*500 controlcard
 376       logical normalize,next
 377       integer restr_type
 378       double precision xlink(4,0:4) /
 379 c           a          b       c     sigma
 380      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
 381      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
 382      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
 383      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
 384      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
 385       write (iout,*) "Calling read_dist_constr"
 386 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 387 c      call flush(iout)
 388       next=.true.
 389
 390       DO WHILE (next)
 391
 392       call card_concat(controlcard,.true.)
 393       next = index(controlcard,"NEXT").gt.0
 394       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
 395       write (iout,*) "restr_type",restr_type
 396       call readi(controlcard,"NFRAG",nfrag_,0)
 397       call readi(controlcard,"NFRAG",nfrag_,0)
 398       call readi(controlcard,"NPAIR",npair_,0)
 399       call readi(controlcard,"NDIST",ndist_,0)
 400       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 401       if (restr_type.eq.10)
 402      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
 403       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 404       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 405       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 406       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 407       normalize = index(controlcard,"NORMALIZE").gt.0
 408       write (iout,*) "WBOLTZD",wboltzd
 409 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 410 c      write (iout,*) "IFRAG"
 411 c      do i=1,nfrag_
 412 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 413 c      enddo
 414 c      write (iout,*) "IPAIR"
 415 c      do i=1,npair_
 416 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 417 c      enddo
 418       if (nfrag_.gt.0) then
 419         nres0=nres
 420         read(inp,'(a)') pdbfile
 421         write (iout,*)
 422      & "Distance restraints will be constructed from structure ",pdbfile
 423         open(ipdbin,file=pdbfile,status='old',err=11)
 424         call readpdb(.true.)
 425         nres=nres0
 426         close(ipdbin)
 427       endif
 428       do i=1,nfrag_
 429         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 430         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 431      &    ifrag_(2,i)=nstart_sup+nsup-1
 432 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 433 c        call flush(iout)
 434         if (wfrag_(i).eq.0.0d0) cycle
 435         do j=ifrag_(1,i),ifrag_(2,i)-1
 436           do k=j+1,ifrag_(2,i)
 437 c            write (iout,*) "j",j," k",k
 438             ddjk=dist(j,k)
 439             if (restr_type.eq.1) then
 440               nhpb=nhpb+1
 441               irestr_type(nhpb)=1
 442               ihpb(nhpb)=j
 443               jhpb(nhpb)=k
 444               dhpb(nhpb)=ddjk
 445               forcon(nhpb)=wfrag_(i)
 446             else if (constr_dist.eq.2) then
 447               if (ddjk.le.dist_cut) then
 448                 nhpb=nhpb+1
 449                 irestr_type(nhpb)=1
 450                 ihpb(nhpb)=j
 451                 jhpb(nhpb)=k
 452                 dhpb(nhpb)=ddjk
 453                 forcon(nhpb)=wfrag_(i)
 454               endif
 455             else if (restr_type.eq.3) then
 456               nhpb=nhpb+1
 457               irestr_type(nhpb)=1
 458               ihpb(nhpb)=j
 459               jhpb(nhpb)=k
 460               dhpb(nhpb)=ddjk
 461               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 462             endif
 463             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
 464      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 465           enddo
 466         enddo
 467       enddo
 468       do i=1,npair_
 469         if (wpair_(i).eq.0.0d0) cycle
 470         ii = ipair_(1,i)
 471         jj = ipair_(2,i)
 472         if (ii.gt.jj) then
 473           itemp=ii
 474           ii=jj
 475           jj=itemp
 476         endif
 477         do j=ifrag_(1,ii),ifrag_(2,ii)
 478           do k=ifrag_(1,jj),ifrag_(2,jj)
 479             if (restr_type.eq.1) then
 480               nhpb=nhpb+1
 481               irestr_type(nhpb)=1
 482               ihpb(nhpb)=j
 483               jhpb(nhpb)=k
 484               dhpb(nhpb)=ddjk
 485               forcon(nhpb)=wfrag_(i)
 486             else if (constr_dist.eq.2) then
 487               if (ddjk.le.dist_cut) then
 488                 nhpb=nhpb+1
 489                 irestr_type(nhpb)=1
 490                 ihpb(nhpb)=j
 491                 jhpb(nhpb)=k
 492                 dhpb(nhpb)=ddjk
 493                 forcon(nhpb)=wfrag_(i)
 494               endif
 495             else if (restr_type.eq.3) then
 496               nhpb=nhpb+1
 497               irestr_type(nhpb)=1
 498               ihpb(nhpb)=j
 499               jhpb(nhpb)=k
 500               dhpb(nhpb)=ddjk
 501               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 502             endif
 503             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
 504      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 505           enddo
 506         enddo
 507       enddo
 508
 509 c      print *,ndist_
 510       write (iout,*) "Distance restraints as read from input"
 511       do i=1,ndist_
 512         if (restr_type.eq.11) then
 513           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
 514      &     dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
 515 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
 516           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
 517           nhpb=nhpb+1
 518           irestr_type(nhpb)=11
 519           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
 520      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
 521      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
 522           if (ibecarb(nhpb).gt.0) then
 523             ihpb(nhpb)=ihpb(nhpb)+nres
 524             jhpb(nhpb)=jhpb(nhpb)+nres
 525           endif
 526         else if (constr_dist.eq.10) then
 527 c Cross-lonk Markov-like potential
 528           call card_concat(controlcard,.true.)
 529           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
 530           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
 531           ibecarb(nhpb+1)=0
 532           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
 533           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
 534           if (index(controlcard,"ZL").gt.0) then
 535             link_type=1
 536           else if (index(controlcard,"ADH").gt.0) then
 537             link_type=2
 538           else if (index(controlcard,"PDH").gt.0) then
 539             link_type=3
 540           else if (index(controlcard,"DSS").gt.0) then
 541             link_type=4
 542           else
 543             link_type=0
 544           endif
 545           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
 546      &       xlink(1,link_type))
 547           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
 548      &       xlink(2,link_type))
 549           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
 550      &       xlink(3,link_type))
 551           call reada(controlcard,"SIGMA",forcon(nhpb+1),
 552      &       xlink(4,link_type))
 553           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
 554 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
 555 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
 556           if (forcon(nhpb+1).le.0.0d0 .or.
 557      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
 558           nhpb=nhpb+1
 559           irestr_type(nhpb)=10
 560           if (ibecarb(nhpb).gt.0) then
 561             ihpb(nhpb)=ihpb(nhpb)+nres
 562             jhpb(nhpb)=jhpb(nhpb)+nres
 563           endif
 564           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
 565      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
 566      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
 567      &     irestr_type(nhpb)
 568         else
 569 C        print *,"in else"
 570           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
 571      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
 572           if (forcon(nhpb+1).gt.0.0d0) then
 573           nhpb=nhpb+1
 574           if (dhpb1(nhpb).eq.0.0d0) then
 575             irestr_type(nhpb)=1
 576           else
 577             irestr_type(nhpb)=2
 578           endif
 579           if (ibecarb(nhpb).gt.0) then
 580             ihpb(nhpb)=ihpb(nhpb)+nres
 581             jhpb(nhpb)=jhpb(nhpb)+nres
 582           endif
 583           if (dhpb(nhpb).eq.0.0d0)
 584      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 585         endif
 586         write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
 587      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
 588         endif
 589 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
 590 C        if (forcon(nhpb+1).gt.0.0d0) then
 591 C          nhpb=nhpb+1
 592 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 593       enddo
 594
 595       ENDDO ! next
 596
 597       fordepthmax=0.0d0
 598       if (normalize) then
 599         do i=nss+1,nhpb
 600           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
 601      &      fordepthmax=fordepth(i)
 602         enddo
 603         do i=nss+1,nhpb
 604           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
 605         enddo
 606       endif
 607       if (nhpb.gt.nss)  then
 608         write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
 609      &  "The following",nhpb-nss,
 610      &  " distance restraints have been imposed:",
 611      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
 612      &  "  score"," type"
 613         do i=nss+1,nhpb
 614           write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
 615      &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
 616      &  irestr_type(i)
 617         enddo
 618       endif
 619       write (iout,*)
 620       call hpb_partition
 621       call flush(iout)
 622       return
 623    11 write (iout,*)"read_dist_restr: error reading reference structure"
 624       stop
 625       end
 626 c====-------------------------------------------------------------------
 627       subroutine read_constr_homology
 628
 629       include 'DIMENSIONS'
 630       include 'DIMENSIONS.ZSCOPT'
 631       include 'DIMENSIONS.FREE'
 632 #ifdef MPI
 633       include 'mpif.h'
 634 #endif
 635       include 'COMMON.SETUP'
 636       include 'COMMON.CONTROL'
 637       include 'COMMON.CHAIN'
 638       include 'COMMON.IOUNITS'
 639       include 'COMMON.GEO'
 640       include 'COMMON.INTERACT'
 641       include 'COMMON.NAMES'
 642       include 'COMMON.HOMRESTR'
 643 c
 644 c For new homol impl
 645 c
 646       include 'COMMON.VAR'
 647 c     include 'include_unres/COMMON.VAR'
 648 c
 649
 650 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
 651 c    &                 dist_cut
 652 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
 653 c    &    sigma_odl_temp(maxres,maxres,max_template)
 654       character*2 kic2
 655       character*24 model_ki_dist, model_ki_angle
 656       character*500 controlcard
 657       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
 658       integer idomain(max_template,maxres)
 659       logical lprn /.true./
 660       integer ilen
 661       external ilen
 662       logical unres_pdb,liiflag
 663 c
 664 c     FP - Nov. 2014 Temporary specifications for new vars
 665 c
 666       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
 667       double precision, dimension (max_template,maxres) :: rescore
 668       double precision, dimension (max_template,maxres) :: rescore2
 669       character*24 tpl_k_rescore
 670 c -----------------------------------------------------------------
 671 c Reading multiple PDB ref structures and calculation of retraints
 672 c not using pre-computed ones stored in files model_ki_{dist,angle}
 673 c FP (Nov., 2014)
 674 c -----------------------------------------------------------------
 675 c
 676 c
 677 c Alternative: reading from input
 678       call card_concat(controlcard,.true.)
 679       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
 680       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
 681       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
 682       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
 683       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
 684       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
 685       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
 686       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
 687       call readi(controlcard,"IHSET",ihset,1)
 688       if (homol_nset.gt.1)then
 689          call card_concat(controlcard,.true.)
 690          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
 691          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 692           write(iout,*) "iset homology_weight "
 693 c         do i=1,homol_nset
 694 c          write(iout,*) i,waga_homology(i)
 695 c         enddo
 696          endif
 697          iset=mod(kolor,homol_nset)+1
 698       else
 699        iset=1
 700        waga_homology(1)=1.0
 701       endif
 702 c     write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
 703
 704 cd      write (iout,*) "nnt",nnt," nct",nct
 705 cd      call flush(iout)
 706
 707
 708       lim_odl=0
 709       lim_dih=0
 710 c
 711 c  New
 712 c
 713       lim_theta=0
 714       lim_xx=0
 715 c
 716 c  Reading HM global scores (prob not required)
 717 c
 718       do i = nnt,nct
 719         do k=1,constr_homology
 720          idomain(k,i)=0
 721         enddo
 722       enddo
 723 c     open (4,file="HMscore")
 724 c     do k=1,constr_homology
 725 c       read (4,*,end=521) hmscore_tmp
 726 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
 727 c       write(*,*) "Model", k, ":", hmscore(k)
 728 c     enddo
 729 c521  continue
 730
 731       ii=0
 732       do i = nnt,nct-2
 733         do j=i+2,nct
 734         ii=ii+1
 735         ii_in_use(ii)=0
 736         enddo
 737       enddo
 738 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
 739
 740       if (read_homol_frag) then
 741         call read_klapaucjusz
 742       else
 743
 744       do k=1,constr_homology
 745
 746         read(inp,'(a)') pdbfile
 747 c  Next stament causes error upon compilation (?)
 748 c       if(me.eq.king.or. .not. out1file)
 749 c         write (iout,'(2a)') 'PDB data will be read from file ',
 750 c    &   pdbfile(:ilen(pdbfile))
 751          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
 752      &  pdbfile(:ilen(pdbfile))
 753         open(ipdbin,file=pdbfile,status='old',err=33)
 754         goto 34
 755   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
 756      &  pdbfile(:ilen(pdbfile))
 757         stop
 758   34    continue
 759 c        print *,'Begin reading pdb data'
 760 c
 761 c Files containing res sim or local scores (former containing sigmas)
 762 c
 763
 764         write(kic2,'(bz,i2.2)') k
 765
 766         tpl_k_rescore="template"//kic2//".sco"
 767
 768         unres_pdb=.false.
 769         call readpdb
 770         do i=1,2*nres
 771           do j=1,3
 772             crefjlee(j,i)=c(j,i)
 773           enddo
 774         enddo
 775 #ifdef DEBUG
 776         do i=1,nres
 777           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
 778      &      (crefjlee(j,i+nres),j=1,3)
 779         enddo
 780 #endif
 781         write (iout,*) "read_constr_homology: after reading pdb file"
 782         call flush(iout)
 783
 784 c
 785 c     Distance restraints
 786 c
 787 c          ... --> odl(k,ii)
 788 C Copy the coordinates from reference coordinates (?)
 789         do i=1,2*nres
 790           do j=1,3
 791             c(j,i)=cref(j,i)
 792 c           write (iout,*) "c(",j,i,") =",c(j,i)
 793           enddo
 794         enddo
 795 c
 796 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
 797 c
 798 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
 799           open (ientin,file=tpl_k_rescore,status='old')
 800           if (nnt.gt.1) rescore(k,1)=0.0d0
 801           do irec=nnt,nct ! loop for reading res sim
 802             if (read2sigma) then
 803              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
 804      &                                idomain_tmp
 805              i_tmp=i_tmp+nnt-1
 806              idomain(k,i_tmp)=idomain_tmp
 807              rescore(k,i_tmp)=rescore_tmp
 808              rescore2(k,i_tmp)=rescore2_tmp
 809              write(iout,'(a7,i5,2f10.5,i5)') "rescore",
 810      &                      i_tmp,rescore2_tmp,rescore_tmp,
 811      &                                idomain_tmp
 812             else
 813              idomain(k,irec)=1
 814              read (ientin,*,end=1401) rescore_tmp
 815
 816 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
 817              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
 818 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
 819             endif
 820           enddo
 821  1401   continue
 822         close (ientin)
 823         if (waga_dist.ne.0.0d0) then
 824           ii=0
 825           do i = nnt,nct-2
 826             do j=i+2,nct
 827
 828               x12=c(1,i)-c(1,j)
 829               y12=c(2,i)-c(2,j)
 830               z12=c(3,i)-c(3,j)
 831               distal=dsqrt(x12*x12+y12*y12+z12*z12)
 832 c              write (iout,*) k,i,j,distal,dist2_cut
 833
 834             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
 835      &            .and. distal.le.dist2_cut ) then
 836
 837               ii=ii+1
 838               ii_in_use(ii)=1
 839               l_homo(k,ii)=.true.
 840
 841 c             write (iout,*) "k",k
 842 c             write (iout,*) "i",i," j",j," constr_homology",
 843 c    &                       constr_homology
 844               ires_homo(ii)=i
 845               jres_homo(ii)=j
 846               odl(k,ii)=distal
 847               if (read2sigma) then
 848                 sigma_odl(k,ii)=0
 849                 do ik=i,j
 850                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
 851                 enddo
 852                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
 853                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
 854      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 855               else
 856                 if (odl(k,ii).le.dist_cut) then
 857                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
 858                 else
 859 #ifdef OLDSIGMA
 860                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 861      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
 862 #else
 863                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 864      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 865 #endif
 866                 endif
 867               endif
 868               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
 869             else
 870               ii=ii+1
 871               l_homo(k,ii)=.false.
 872             endif
 873             enddo
 874           enddo
 875         lim_odl=ii
 876         endif
 877 c
 878 c     Theta, dihedral and SC retraints
 879 c
 880         if (waga_angle.gt.0.0d0) then
 881 c         open (ientin,file=tpl_k_sigma_dih,status='old')
 882 c         do irec=1,maxres-3 ! loop for reading sigma_dih
 883 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
 884 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
 885 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
 886 c    &                            sigma_dih(k,i+nnt-1)
 887 c         enddo
 888 c1402   continue
 889 c         close (ientin)
 890           do i = nnt+3,nct
 891             if (idomain(k,i).eq.0) then
 892                sigma_dih(k,i)=0.0
 893                cycle
 894             endif
 895             dih(k,i)=phiref(i) ! right?
 896 c           read (ientin,*) sigma_dih(k,i) ! original variant
 897 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
 898 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
 899 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 900 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
 901 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
 902
 903             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
 904      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
 905 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
 906 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
 907 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
 908 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 909 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
 910 c   Instead of res sim other local measure of b/b str reliability possible
 911             if (sigma_dih(k,i).ne.0)
 912      &      sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
 913 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
 914           enddo
 915           lim_dih=nct-nnt-2
 916         endif
 917
 918         if (waga_theta.gt.0.0d0) then
 919 c         open (ientin,file=tpl_k_sigma_theta,status='old')
 920 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
 921 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
 922 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
 923 c    &                              sigma_theta(k,i+nnt-1)
 924 c         enddo
 925 c1403   continue
 926 c         close (ientin)
 927
 928           do i = nnt+2,nct ! right? without parallel.
 929 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
 930 c         do i=ithet_start,ithet_end ! with FG parallel.
 931              if (idomain(k,i).eq.0) then
 932               sigma_theta(k,i)=0.0
 933               cycle
 934              endif
 935              thetatpl(k,i)=thetaref(i)
 936 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
 937 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
 938 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 939 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
 940 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
 941              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
 942      &                        rescore(k,i-2))/3.0
 943 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
 944              if (sigma_theta(k,i).ne.0)
 945      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 946
 947 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
 948 c                             rescore(k,i-2) !  right expression ?
 949 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
 950           enddo
 951         endif
 952
 953         if (waga_d.gt.0.0d0) then
 954 c       open (ientin,file=tpl_k_sigma_d,status='old')
 955 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
 956 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
 957 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
 958 c    &                          sigma_d(k,i+nnt-1)
 959 c         enddo
 960 c1404   continue
 961
 962           do i = nnt,nct ! right? without parallel.
 963 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
 964 c         do i=loc_start,loc_end ! with FG parallel.
 965                if (itype(i).eq.10) cycle
 966                if (idomain(k,i).eq.0 ) then
 967                   sigma_d(k,i)=0.0
 968                   cycle
 969                endif
 970                xxtpl(k,i)=xxref(i)
 971                yytpl(k,i)=yyref(i)
 972                zztpl(k,i)=zzref(i)
 973 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
 974 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
 975 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
 976 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
 977                sigma_d(k,i)=rescore(k,i) !  right expression ?
 978                if (sigma_d(k,i).ne.0)
 979      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
 980
 981 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
 982 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
 983 c              read (ientin,*) sigma_d(k,i) ! 1st variant
 984           enddo
 985         endif
 986       enddo
 987 c
 988 c remove distance restraints not used in any model from the list
 989 c shift data in all arrays
 990 c
 991       if (waga_dist.ne.0.0d0) then
 992         ii=0
 993         liiflag=.true.
 994         do i=nnt,nct-2
 995          do j=i+2,nct
 996           ii=ii+1
 997           if (ii_in_use(ii).eq.0.and.liiflag) then
 998             liiflag=.false.
 999             iistart=ii
1000           endif
1001           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1002      &                   .not.liiflag.and.ii.eq.lim_odl) then
1003              if (ii.eq.lim_odl) then
1004               iishift=ii-iistart+1
1005              else
1006               iishift=ii-iistart
1007              endif
1008              liiflag=.true.
1009              do ki=iistart,lim_odl-iishift
1010               ires_homo(ki)=ires_homo(ki+iishift)
1011               jres_homo(ki)=jres_homo(ki+iishift)
1012               ii_in_use(ki)=ii_in_use(ki+iishift)
1013               do k=1,constr_homology
1014                odl(k,ki)=odl(k,ki+iishift)
1015                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1016                l_homo(k,ki)=l_homo(k,ki+iishift)
1017               enddo
1018              enddo
1019              ii=ii-iishift
1020              lim_odl=lim_odl-iishift
1021           endif
1022          enddo
1023         enddo
1024       endif
1025
1026       endif ! .not. klapaucjusz
1027
1028       if (constr_homology.gt.0) call homology_partition
1029       if (constr_homology.gt.0) call init_int_table
1030 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1031 cd     & "lim_xx=",lim_xx
1032 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1033 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1034 c
1035 c Print restraints
1036 c
1037       if (.not.lprn) return
1038 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1039       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1040        write (iout,*) "Distance restraints from templates"
1041        do ii=1,lim_odl
1042        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1043      &  ii,ires_homo(ii),jres_homo(ii),
1044      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1045      &  ki=1,constr_homology)
1046        enddo
1047        write (iout,*) "Dihedral angle restraints from templates"
1048        do i=nnt+3,nct
1049         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1050      &      (rad2deg*dih(ki,i),
1051      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1052        enddo
1053        write (iout,*) "Virtual-bond angle restraints from templates"
1054        do i=nnt+2,nct
1055         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1056      &      (rad2deg*thetatpl(ki,i),
1057      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1058        enddo
1059        write (iout,*) "SC restraints from templates"
1060        do i=nnt,nct
1061         write(iout,'(i5,100(4f8.2,4x))') i,
1062      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1063      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1064        enddo
1065       endif
1066 c -----------------------------------------------------------------
1067       return
1068       end
1069 c----------------------------------------------------------------------
1070       subroutine read_klapaucjusz
1071
1072       include 'DIMENSIONS'
1073       include 'DIMENSIONS.ZSCOPT'
1074       include 'DIMENSIONS.FREE'
1075 #ifdef MPI
1076       include 'mpif.h'
1077 #endif
1078       include 'COMMON.SETUP'
1079       include 'COMMON.CONTROL'
1080       include 'COMMON.CHAIN'
1081       include 'COMMON.IOUNITS'
1082       include 'COMMON.GEO'
1083       include 'COMMON.INTERACT'
1084       include 'COMMON.NAMES'
1085       include 'COMMON.HOMRESTR'
1086       character*256 fragfile
1087       integer ninclust(maxclust),inclust(max_template,maxclust),
1088      &  nresclust(maxclust),iresclust(maxres,maxclust)
1089
1090       character*2 kic2
1091       character*24 model_ki_dist, model_ki_angle
1092       character*500 controlcard
1093       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1094       integer idomain(max_template,maxres)
1095       logical lprn /.true./
1096       integer ilen
1097       external ilen
1098       logical unres_pdb,liiflag
1099 c
1100 c
1101       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
1102       double precision, dimension (max_template,maxres) :: rescore
1103       double precision, dimension (max_template,maxres) :: rescore2
1104       character*24 tpl_k_rescore
1105
1106 c
1107 c For new homol impl
1108 c
1109       include 'COMMON.VAR'
1110 c
1111       double precision chomo(3,maxres2+2,max_template)
1112       call getenv("FRAGFILE",fragfile)
1113       open(ientin,file=fragfile,status="old",err=10)
1114       read(ientin,*) constr_homology,nclust
1115       l_homo = .false.
1116       sigma_theta=0.0
1117       sigma_d=0.0
1118       sigma_dih=0.0
1119 c Read pdb files
1120       do k=1,constr_homology
1121         read(ientin,'(a)') pdbfile
1122         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
1123      &  pdbfile(:ilen(pdbfile))
1124         open(ipdbin,file=pdbfile,status='old',err=33)
1125         goto 34
1126   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
1127      &  pdbfile(:ilen(pdbfile))
1128         stop
1129   34    continue
1130         unres_pdb=.false.
1131         call readpdb
1132         do i=1,2*nres
1133           do j=1,3
1134             chomo(j,i,k)=c(j,i)
1135           enddo
1136         enddo
1137         do i=1,nres
1138           rescore(k,i)=0.2d0
1139           rescore2(k,i)=1.0d0
1140         enddo
1141       enddo
1142 c Read clusters
1143       do i=1,nclust
1144         read(ientin,*) ninclust(i),nresclust(i)
1145         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
1146         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
1147       enddo
1148 c
1149 c Loop over clusters
1150 c
1151       do l=1,nclust
1152         do ll = 1,ninclust(l)
1153
1154         k = inclust(ll,l)
1155         do i=1,nres
1156           idomain(k,i)=0
1157         enddo
1158         do i=1,nresclust(l)
1159           if (nnt.gt.1)  then
1160             idomain(k,iresclust(i,l)+1) = 1
1161           else
1162             idomain(k,iresclust(i,l)) = 1
1163           endif
1164         enddo
1165 c
1166 c     Distance restraints
1167 c
1168 c          ... --> odl(k,ii)
1169 C Copy the coordinates from reference coordinates (?)
1170         do i=1,2*nres
1171           do j=1,3
1172             c(j,i)=chomo(j,i,k)
1173 c           write (iout,*) "c(",j,i,") =",c(j,i)
1174           enddo
1175         enddo
1176         call int_from_cart(.true.,.false.)
1177         call sc_loc_geom(.false.)
1178         do i=1,nres
1179           thetaref(i)=theta(i)
1180           phiref(i)=phi(i)
1181         enddo
1182         if (waga_dist.ne.0.0d0) then
1183           ii=0
1184           do i = nnt,nct-2
1185             do j=i+2,nct
1186
1187               x12=c(1,i)-c(1,j)
1188               y12=c(2,i)-c(2,j)
1189               z12=c(3,i)-c(3,j)
1190               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1191 c              write (iout,*) k,i,j,distal,dist2_cut
1192
1193             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1194      &            .and. distal.le.dist2_cut ) then
1195
1196               ii=ii+1
1197               ii_in_use(ii)=1
1198               l_homo(k,ii)=.true.
1199
1200 c             write (iout,*) "k",k
1201 c             write (iout,*) "i",i," j",j," constr_homology",
1202 c    &                       constr_homology
1203               ires_homo(ii)=i
1204               jres_homo(ii)=j
1205               odl(k,ii)=distal
1206               if (read2sigma) then
1207                 sigma_odl(k,ii)=0
1208                 do ik=i,j
1209                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1210                 enddo
1211                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1212                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1213      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1214               else
1215                 if (odl(k,ii).le.dist_cut) then
1216                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1217                 else
1218 #ifdef OLDSIGMA
1219                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1220      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1221 #else
1222                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1223      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1224 #endif
1225                 endif
1226               endif
1227               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1228             else
1229               ii=ii+1
1230 c              l_homo(k,ii)=.false.
1231             endif
1232             enddo
1233           enddo
1234         lim_odl=ii
1235         endif
1236 c
1237 c     Theta, dihedral and SC retraints
1238 c
1239         if (waga_angle.gt.0.0d0) then
1240           do i = nnt+3,nct
1241             if (idomain(k,i).eq.0) then
1242 c               sigma_dih(k,i)=0.0
1243                cycle
1244             endif
1245             dih(k,i)=phiref(i)
1246             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1247      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1248 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
1249 c     &       " sigma_dihed",sigma_dih(k,i)
1250             if (sigma_dih(k,i).ne.0)
1251      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1252           enddo
1253           lim_dih=nct-nnt-2
1254         endif
1255
1256         if (waga_theta.gt.0.0d0) then
1257           do i = nnt+2,nct
1258              if (idomain(k,i).eq.0) then
1259 c              sigma_theta(k,i)=0.0
1260               cycle
1261              endif
1262              thetatpl(k,i)=thetaref(i)
1263              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1264      &                        rescore(k,i-2))/3.0
1265              if (sigma_theta(k,i).ne.0)
1266      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1267           enddo
1268         endif
1269
1270         if (waga_d.gt.0.0d0) then
1271           do i = nnt,nct
1272                if (itype(i).eq.10) cycle
1273                if (idomain(k,i).eq.0 ) then
1274 c                  sigma_d(k,i)=0.0
1275                   cycle
1276                endif
1277                xxtpl(k,i)=xxref(i)
1278                yytpl(k,i)=yyref(i)
1279                zztpl(k,i)=zzref(i)
1280                sigma_d(k,i)=rescore(k,i)
1281                if (sigma_d(k,i).ne.0)
1282      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1283                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
1284           enddo
1285         endif
1286       enddo ! l
1287       enddo ! ll
1288 c
1289 c remove distance restraints not used in any model from the list
1290 c shift data in all arrays
1291 c
1292       if (waga_dist.ne.0.0d0) then
1293         ii=0
1294         liiflag=.true.
1295         do i=nnt,nct-2
1296          do j=i+2,nct
1297           ii=ii+1
1298           if (ii_in_use(ii).eq.0.and.liiflag) then
1299             liiflag=.false.
1300             iistart=ii
1301           endif
1302           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1303      &                   .not.liiflag.and.ii.eq.lim_odl) then
1304              if (ii.eq.lim_odl) then
1305               iishift=ii-iistart+1
1306              else
1307               iishift=ii-iistart
1308              endif
1309              liiflag=.true.
1310              do ki=iistart,lim_odl-iishift
1311               ires_homo(ki)=ires_homo(ki+iishift)
1312               jres_homo(ki)=jres_homo(ki+iishift)
1313               ii_in_use(ki)=ii_in_use(ki+iishift)
1314               do k=1,constr_homology
1315                odl(k,ki)=odl(k,ki+iishift)
1316                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1317                l_homo(k,ki)=l_homo(k,ki+iishift)
1318               enddo
1319              enddo
1320              ii=ii-iishift
1321              lim_odl=lim_odl-iishift
1322           endif
1323          enddo
1324         enddo
1325       endif
1326
1327       return
1328    10 stop "Error infragment file"
1329       end
1330 c----------------------------------------------------------------------