source/wham/src/molread_zs.F

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