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