source/cluster/wham/src/readrtns.F

   1       subroutine read_control
   2 C
   3 C Read molecular data
   4 C
   5       implicit none
   6       include 'DIMENSIONS'
   7       include 'sizesclu.dat'
   8       include 'COMMON.IOUNITS'
   9       include 'COMMON.TIME1'
  10       include 'COMMON.SBRIDGE'
  11       include 'COMMON.CONTROL'
  12       include 'COMMON.CLUSTER'
  13       include 'COMMON.CHAIN'
  14       include 'COMMON.HEADER'
  15       include 'COMMON.FFIELD'
  16       include 'COMMON.FREE'
  17       character*320 controlcard,ucase
  18 #ifdef MPL
  19       include 'COMMON.INFO'
  20 #endif
  21       integer i
  22
  23       read (INP,'(a80)') titel
  24       call card_concat(controlcard)
  25
  26       call readi(controlcard,'NRES',nres,0)
  27       write (iout,*) "NRES",NRES
  28       call readi(controlcard,'RESCALE',rescale_mode,2)
  29       call readi(controlcard,'PDBOUT',outpdb,0)
  30       call readi(controlcard,'MOL2OUT',outmol2,0)
  31       refstr=(index(controlcard,'REFSTR').gt.0)
  32       write (iout,*) "REFSTR",refstr
  33       pdbref=(index(controlcard,'PDBREF').gt.0)
  34       dyn_ss=(index(controlcard,'DYN_SS').gt.0)
  35       iscode=index(controlcard,'ONE_LETTER')
  36       tree=(index(controlcard,'MAKE_TREE').gt.0)
  37       min_var=(index(controlcard,'MINVAR').gt.0)
  38       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
  39       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
  40       call readi(controlcard,'NCUT',ncut,0)
  41       call readi(controlcard,'NCLUST',nclust,5)
  42       call readi(controlcard,'NSTART',nstart,0)
  43       call readi(controlcard,'NEND',nend,0)
  44       call reada(controlcard,'ECUT',ecut,10.0d0)
  45       call reada(controlcard,'PROB',prob_limit,0.99d0)
  46       write (iout,*) "Probability limit",prob_limit
  47       lgrp=(index(controlcard,'LGRP').gt.0)
  48       caonly=(index(controlcard,'CA_ONLY').gt.0)
  49       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
  50       if (ncut.gt.0)
  51      & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
  52       call readi(controlcard,'IOPT',iopt,2)
  53       lside = index(controlcard,"SIDE").gt.0
  54       efree = index(controlcard,"EFREE").gt.0
  55       call readi(controlcard,'NTEMP',nT,1)
  56       write (iout,*) "nT",nT
  57       call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
  58       write (iout,*) "nT",nT
  59       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
  60       do i=1,nT
  61         beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
  62       enddo
  63       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
  64       lprint_cart=index(controlcard,"PRINT_CART") .gt.0
  65       lprint_int=index(controlcard,"PRINT_INT") .gt.0
  66       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
  67       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  68       write (iout,*) "with_dihed_constr ",with_dihed_constr,
  69      & " CONSTR_DIST",constr_dist
  70
  71       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
  72       write (iout,*) "with_homology_constr ",with_dihed_constr,
  73      & " CONSTR_HOMOLOGY",constr_homology
  74       print_homology_restraints=
  75      & index(controlcard,"PRINT_HOMOLOGY_RESTRAINTS").gt.0
  76       print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0
  77       print_homology_models=
  78      & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0
  79       read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
  80
  81 #ifdef AIX
  82       call flush_(iout)
  83 #else
  84       call flush(iout)
  85 #endif
  86       if (min_var) iopt=1
  87       return
  88       end
  89 c--------------------------------------------------------------------------
  90       subroutine molread
  91 C
  92 C Read molecular data.
  93 C
  94       implicit none
  95       include 'DIMENSIONS'
  96       include 'COMMON.IOUNITS'
  97       include 'COMMON.GEO'
  98       include 'COMMON.VAR'
  99       include 'COMMON.INTERACT'
 100       include 'COMMON.LOCAL'
 101       include 'COMMON.NAMES'
 102       include 'COMMON.CHAIN'
 103       include 'COMMON.FFIELD'
 104       include 'COMMON.SBRIDGE'
 105       include 'COMMON.HEADER'
 106       include 'COMMON.CONTROL'
 107       include 'COMMON.CONTACTS'
 108       include 'COMMON.TIME1'
 109       include 'COMMON.TORCNSTR'
 110 #ifdef MPL
 111       include 'COMMON.INFO'
 112 #endif
 113       character*4 sequence(maxres)
 114       character*800 weightcard
 115       integer rescode
 116       double precision x(maxvar)
 117       integer itype_pdb(maxres)
 118       logical seq_comp
 119       integer i,j
 120       write (iout,*) " MOLREAD: NRES",NRES
 121 C
 122 C Body
 123 C
 124 C Read weights of the subsequent energy terms.
 125       call card_concat(weightcard)
 126       call reada(weightcard,'WSC',wsc,1.0d0)
 127       call reada(weightcard,'WLONG',wsc,wsc)
 128       call reada(weightcard,'WSCP',wscp,1.0d0)
 129       call reada(weightcard,'WELEC',welec,1.0D0)
 130       call reada(weightcard,'WVDWPP',wvdwpp,welec)
 131       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 132       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 133       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 134       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 135       call reada(weightcard,'WTURN3',wturn3,1.0D0)
 136       call reada(weightcard,'WTURN4',wturn4,1.0D0)
 137       call reada(weightcard,'WTURN6',wturn6,1.0D0)
 138       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 139       call reada(weightcard,'WBOND',wbond,1.0D0)
 140       call reada(weightcard,'WTOR',wtor,1.0D0)
 141       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 142       call reada(weightcard,'WANG',wang,1.0D0)
 143       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 144       call reada(weightcard,'SCAL14',scal14,0.4D0)
 145       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 146       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 147       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 148       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 149       call reada(weightcard,"D0CM",d0cm,3.78d0)
 150       call reada(weightcard,"AKCM",akcm,15.1d0)
 151       call reada(weightcard,"AKTH",akth,11.0d0)
 152       call reada(weightcard,"AKCT",akct,12.0d0)
 153       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 154       call reada(weightcard,"V2SS",v2ss,7.61d0)
 155       call reada(weightcard,"V3SS",v3ss,13.7d0)
 156       call reada(weightcard,"EBR",ebr,-5.50D0)
 157 C     Bartek
 158       call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
 159       call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
 160       call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
 161       call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
 162       if (index(weightcard,'SOFT').gt.0) ipot=6
 163 C 12/1/95 Added weight for the multi-body term WCORR
 164       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 165       do i=1,maxres-1
 166         do j=i+1,maxres
 167           dyn_ssbond_ij(i,j)=1.0d300
 168         enddo
 169       enddo
 170       call reada(weightcard,"HT",Ht,0.0D0)
 171       if (dyn_ss) then
 172         ss_depth=ebr/wsc-0.25*eps(1,1)
 173         Ht=Ht/wsc-0.25*eps(1,1)
 174         akcm=akcm*wstrain/wsc
 175         akth=akth*wstrain/wsc
 176         akct=akct*wstrain/wsc
 177         v1ss=v1ss*wstrain/wsc
 178         v2ss=v2ss*wstrain/wsc
 179         v3ss=v3ss*wstrain/wsc
 180       else
 181         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 182       endif
 183       write (iout,'(/a)') "Disulfide bridge parameters:"
 184       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 185         write (iout,'(a,f10.2)') 'S-S depth: ',ss_depth
 186       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 187       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 188       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 189      & ' v3ss:',v3ss
 190       write (iout,'(2(a,f10.2))') 'ht:',ht,' eps:', eps(1,1)
 191       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 192       weights(1)=wsc
 193       weights(2)=wscp
 194       weights(3)=welec
 195       weights(4)=wcorr
 196       weights(5)=wcorr5
 197       weights(6)=wcorr6
 198       weights(7)=wel_loc
 199       weights(8)=wturn3
 200       weights(9)=wturn4
 201       weights(10)=wturn6
 202       weights(11)=wang
 203       weights(12)=wscloc
 204       weights(13)=wtor
 205       weights(14)=wtor_d
 206       weights(15)=wstrain
 207       weights(16)=wvdwpp
 208       weights(17)=wbond
 209       weights(18)=scal14
 210       weights(22)=wdfa_dist
 211       weights(23)=wdfa_tor
 212       weights(24)=wdfa_nei
 213       weights(25)=wdfa_beta
 214       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 215      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
 216      &  wturn4,wturn6,wsccor,wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
 217    10 format (/'Energy-term weights (unscaled):'//
 218      & 'WSCC=   ',f10.6,' (SC-SC)'/
 219      & 'WSCP=   ',f10.6,' (SC-p)'/
 220      & 'WELEC=  ',f10.6,' (p-p electr)'/
 221      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 222      & 'WBOND=  ',f10.6,' (stretching)'/
 223      & 'WANG=   ',f10.6,' (bending)'/
 224      & 'WSCLOC= ',f10.6,' (SC local)'/
 225      & 'WTOR=   ',f10.6,' (torsional)'/
 226      & 'WTORD=  ',f10.6,' (double torsional)'/
 227      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 228      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 229      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 230      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 231      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 232      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 233      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 234      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 235      & 'WSCCOR= ',f10.6,' (SC-backbone torsional correalations)'/
 236      & 'WDFAD=  ',f10.6,' (DFA distance)'/
 237      & 'WDFAT=  ',f10.6,' (DFA torsional)'/
 238      & 'WDFAN=  ',f10.6,' (DFA neighbors)'/
 239      & 'WDFAB=  ',f10.6,' (DFA beta)'/)
 240       if (wcorr4.gt.0.0d0) then
 241         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 242      &   'between contact pairs of peptide groups'
 243         write (iout,'(2(a,f5.3/))')
 244      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 245      &  'Range of quenching the correlation terms:',2*delt_corr
 246       else if (wcorr.gt.0.0d0) then
 247         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 248      &   'between contact pairs of peptide groups'
 249       endif
 250       write (iout,'(a,f8.3)')
 251      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 252       write (iout,'(a,f8.3)')
 253      &  'General scaling factor of SC-p interactions:',scalscp
 254       r0_corr=cutoff_corr-delt_corr
 255       do i=1,20
 256         aad(i,1)=scalscp*aad(i,1)
 257         aad(i,2)=scalscp*aad(i,2)
 258         bad(i,1)=scalscp*bad(i,1)
 259         bad(i,2)=scalscp*bad(i,2)
 260       enddo
 261
 262 #ifdef AIX
 263       call flush_(iout)
 264 #else
 265       call flush(iout)
 266 #endif
 267       print *,'indpdb=',indpdb,' pdbref=',pdbref
 268
 269 C Read sequence if not taken from the pdb file.
 270       if (iscode.gt.0) then
 271         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 272       else
 273         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 274       endif
 275 C Convert sequence to numeric code
 276       do i=1,nres
 277         itype(i)=rescode(i,sequence(i),iscode)
 278       enddo
 279       if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
 280         write (iout,*)
 281      &   "Glycine is the first full residue, initial dummy deleted"
 282         do i=1,nres
 283           itype(i)=itype(i+1)
 284         enddo
 285         nres=nres-1
 286       endif
 287       if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
 288         write (iout,*)
 289      &   "Glycine is the last full residue, terminal dummy deleted"
 290         nres=nres-1
 291       endif
 292       print *,nres
 293       print '(20i4)',(itype(i),i=1,nres)
 294
 295       do i=1,nres
 296 #ifdef PROCOR
 297         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 298 #else
 299         if (itype(i).eq.21) then
 300 #endif
 301           itel(i)=0
 302 #ifdef PROCOR
 303         else if (itype(i+1).ne.20) then
 304 #else
 305         else if (itype(i).ne.20) then
 306 #endif
 307           itel(i)=1
 308         else
 309           itel(i)=2
 310         endif
 311       enddo
 312       write (iout,*) "ITEL"
 313       do i=1,nres-1
 314         write (iout,*) i,itype(i),itel(i)
 315       enddo
 316
 317       print *,'Call Read_Bridge.'
 318       call read_bridge
 319       if (with_dihed_constr) then
 320
 321       read (inp,*) ndih_constr
 322       if (ndih_constr.gt.0) then
 323         read (inp,*) ftors
 324         write (iout,*) 'FTORS',ftors
 325         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 326         write (iout,*)
 327      &   'There are',ndih_constr,' constraints on phi angles.'
 328         do i=1,ndih_constr
 329           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 330         enddo
 331         do i=1,ndih_constr
 332           phi0(i)=deg2rad*phi0(i)
 333           drange(i)=deg2rad*drange(i)
 334         enddo
 335       endif
 336
 337       endif
 338
 339       nnt=1
 340       nct=nres
 341       print *,'NNT=',NNT,' NCT=',NCT
 342       if (itype(1).eq.21) nnt=2
 343       if (itype(nres).eq.21) nct=nct-1
 344       if (nstart.lt.nnt) nstart=nnt
 345       if (nend.gt.nct .or. nend.eq.0) nend=nct
 346       write (iout,*) "nstart",nstart," nend",nend
 347       nres0=nres
 348
 349 C     Juyong:READ init_vars
 350 C     Initialize variables!
 351 C     Juyong:READ read_info
 352 C     READ fragment information!!
 353 C     both routines should be in dfa.F file!!
 354
 355       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
 356      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 357 #ifdef DEBUG
 358        write (iout,*) "Calling init_dfa_vars"
 359 #ifdef AIX
 360       call flush_(iout)
 361 #else
 362       call flush(iout)
 363 #endif
 364 #endif
 365        call init_dfa_vars
 366 #ifdef DEBUG
 367        write (iout,*) 'init_dfa_vars finished!'
 368 #ifdef AIX
 369       call flush_(iout)
 370 #else
 371       call flush(iout)
 372 #endif
 373 #endif
 374        call read_dfa_info
 375 #ifdef DEBUG
 376        write (iout,*) 'read_dfa_info finished!'
 377 #ifdef AIX
 378       call flush_(iout)
 379 #else
 380       call flush(iout)
 381 #endif
 382 #endif
 383       endif
 384
 385       if (constr_homology.gt.0) then
 386         call read_constr_homology(print_homology_restraints)
 387       endif
 388
 389 c      if (pdbref) then
 390 c        read(inp,'(a)') pdbfile
 391 c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
 392 c        open(ipdbin,file=pdbfile,status='old',err=33)
 393 c        goto 34
 394 c  33    write (iout,'(a)') 'Error opening PDB file.'
 395 c        stop
 396 c  34    continue
 397 c        print *,'Begin reading pdb data'
 398 c        call readpdb
 399 c        print *,'Finished reading pdb data'
 400 c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
 401 c        do i=1,nres
 402 c          itype_pdb(i)=itype(i)
 403 c        enddo
 404 c        close (ipdbin)
 405 c        write (iout,'(a,i3)') 'nsup=',nsup
 406 c        nstart_seq=nnt
 407 c        if (nsup.le.(nct-nnt+1)) then
 408 c          do i=0,nct-nnt+1-nsup
 409 c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 410 c              nstart_seq=nnt+i
 411 c              goto 111
 412 c            endif
 413 c          enddo
 414 c          write (iout,'(a)')
 415 c     &            'Error - sequences to be superposed do not match.'
 416 c          stop
 417 c        else
 418 c          do i=0,nsup-(nct-nnt+1)
 419 c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 420 c     &      then
 421 c              nstart_sup=nstart_sup+i
 422 c              nsup=nct-nnt+1
 423 c              goto 111
 424 c            endif
 425 c          enddo
 426 c          write (iout,'(a)')
 427 c     &            'Error - sequences to be superposed do not match.'
 428 c        endif
 429 c  111   continue
 430 c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 431 c     &                 ' nstart_seq=',nstart_seq
 432 c      endif
 433       call init_int_table
 434       call setup_var
 435       write (iout,*) "molread: REFSTR",refstr
 436       if (refstr) then
 437         if (.not.pdbref) then
 438           call read_angles(inp,*38)
 439           goto 39
 440    38     write (iout,'(a)') 'Error reading reference structure.'
 441 #ifdef MPL
 442           call mp_stopall(Error_Msg)
 443 #else
 444           stop 'Error reading reference structure'
 445 #endif
 446    39     call chainbuild
 447           nstart_sup=nnt
 448           nstart_seq=nnt
 449           nsup=nct-nnt+1
 450           do i=1,2*nres
 451             do j=1,3
 452               cref(j,i)=c(j,i)
 453             enddo
 454           enddo
 455         endif
 456         call contact(print_contact_map,ncont_ref,icont_ref)
 457       endif
 458 c Read distance restraints
 459       if (constr_dist.gt.0) then
 460         call read_dist_constr
 461         call hpb_partition
 462       endif
 463       return
 464       end
 465 c-----------------------------------------------------------------------------
 466       logical function seq_comp(itypea,itypeb,length)
 467       implicit none
 468       integer length,itypea(length),itypeb(length)
 469       integer i
 470       do i=1,length
 471         if (itypea(i).ne.itypeb(i)) then
 472           seq_comp=.false.
 473           return
 474         endif
 475       enddo
 476       seq_comp=.true.
 477       return
 478       end
 479 c-----------------------------------------------------------------------------
 480       subroutine read_bridge
 481 C Read information about disulfide bridges.
 482       implicit none
 483       include 'DIMENSIONS'
 484       include 'COMMON.IOUNITS'
 485       include 'COMMON.GEO'
 486       include 'COMMON.VAR'
 487       include 'COMMON.INTERACT'
 488       include 'COMMON.LOCAL'
 489       include 'COMMON.NAMES'
 490       include 'COMMON.CHAIN'
 491       include 'COMMON.FFIELD'
 492       include 'COMMON.SBRIDGE'
 493       include 'COMMON.HEADER'
 494       include 'COMMON.CONTROL'
 495       include 'COMMON.TIME1'
 496 #ifdef MPL
 497       include 'COMMON.INFO'
 498 #endif
 499       integer i,j
 500 C Read bridging residues.
 501       read (inp,*) ns,(iss(i),i=1,ns)
 502       write(iout,*)'ns=',ns
 503 C Check whether the specified bridging residues are cystines.
 504       do i=1,ns
 505         if (itype(iss(i)).ne.1) then
 506           write (iout,'(2a,i3,a)')
 507      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 508      &   ' can form a disulfide bridge?!!!'
 509           write (*,'(2a,i3,a)')
 510      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 511      &   ' can form a disulfide bridge?!!!'
 512 #ifdef MPL
 513          call mp_stopall(error_msg)
 514 #else
 515          stop
 516 #endif
 517         endif
 518       enddo
 519 C Read preformed bridges.
 520       if (ns.gt.0) then
 521       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 522       if (nss.gt.0) then
 523         nhpb=nss
 524 C Check if the residues involved in bridges are in the specified list of
 525 C bridging residues.
 526         do i=1,nss
 527           do j=1,i-1
 528             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 529      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 530               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 531      &      ' contains residues present in other pairs.'
 532               write (*,'(a,i3,a)') 'Disulfide pair',i,
 533      &      ' contains residues present in other pairs.'
 534 #ifdef MPL
 535               call mp_stopall(error_msg)
 536 #else
 537               stop
 538 #endif
 539             endif
 540           enddo
 541           do j=1,ns
 542             if (ihpb(i).eq.iss(j)) goto 10
 543           enddo
 544           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 545    10     continue
 546           do j=1,ns
 547             if (jhpb(i).eq.iss(j)) goto 20
 548           enddo
 549           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 550    20     continue
 551 c          dhpb(i)=dbr
 552 c          forcon(i)=fbr
 553         enddo
 554         do i=1,nss
 555           ihpb(i)=ihpb(i)+nres
 556           jhpb(i)=jhpb(i)+nres
 557         enddo
 558       endif
 559       endif
 560       if (ns.gt.0.and.dyn_ss) then
 561           do i=nss+1,nhpb
 562             ihpb(i-nss)=ihpb(i)
 563             jhpb(i-nss)=jhpb(i)
 564             forcon(i-nss)=forcon(i)
 565             dhpb(i-nss)=dhpb(i)
 566           enddo
 567           nhpb=nhpb-nss
 568           nss=0
 569           call hpb_partition
 570           do i=1,ns
 571             dyn_ss_mask(iss(i))=.true.
 572 c          write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
 573           enddo
 574       endif
 575       print *, "Leaving brigde read"
 576       return
 577       end
 578 c----------------------------------------------------------------------------
 579       subroutine read_angles(kanal,*)
 580       implicit none
 581       include 'DIMENSIONS'
 582       include 'COMMON.GEO'
 583       include 'COMMON.VAR'
 584       include 'COMMON.CHAIN'
 585       include 'COMMON.IOUNITS'
 586       integer i,kanal
 587       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 588       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 589       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 590       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 591       do i=1,nres
 592         theta(i)=deg2rad*theta(i)
 593         phi(i)=deg2rad*phi(i)
 594         alph(i)=deg2rad*alph(i)
 595         omeg(i)=deg2rad*omeg(i)
 596       enddo
 597       return
 598    10 return1
 599       end
 600 c----------------------------------------------------------------------------
 601       subroutine reada(rekord,lancuch,wartosc,default)
 602       implicit none
 603       character*(*) rekord,lancuch
 604       double precision wartosc,default
 605       integer ilen,iread
 606       external ilen
 607       iread=index(rekord,lancuch)
 608       if (iread.eq.0) then
 609         wartosc=default
 610         return
 611       endif
 612       iread=iread+ilen(lancuch)+1
 613       read (rekord(iread:),*) wartosc
 614       return
 615       end
 616 c----------------------------------------------------------------------------
 617       subroutine multreada(rekord,lancuch,tablica,dim,default)
 618       implicit none
 619       integer dim,i
 620       double precision tablica(dim),default
 621       character*(*) rekord,lancuch
 622       integer ilen,iread
 623       external ilen
 624       do i=1,dim
 625         tablica(i)=default
 626       enddo
 627       iread=index(rekord,lancuch)
 628       if (iread.eq.0) return
 629       iread=iread+ilen(lancuch)+1
 630       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 631    10 return
 632       end
 633 c----------------------------------------------------------------------------
 634       subroutine readi(rekord,lancuch,wartosc,default)
 635       implicit none
 636       character*(*) rekord,lancuch
 637       integer wartosc,default
 638       integer ilen,iread
 639       external ilen
 640       iread=index(rekord,lancuch)
 641       if (iread.eq.0) then
 642         wartosc=default
 643         return
 644       endif
 645       iread=iread+ilen(lancuch)+1
 646       read (rekord(iread:),*) wartosc
 647       return
 648       end
 649 c----------------------------------------------------------------------------
 650       subroutine multreadi(rekord,lancuch,tablica,dim,default)
 651       implicit none
 652       integer dim,i
 653       integer tablica(dim),default
 654       character*(*) rekord,lancuch
 655       character*80 aux
 656       integer ilen,iread
 657       external ilen
 658       do i=1,dim
 659         tablica(i)=default
 660       enddo
 661       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 662       if (iread.eq.0) return
 663       iread=iread+ilen(lancuch)+1
 664       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 665    10 return
 666       end
 667 c----------------------------------------------------------------------------
 668       subroutine card_concat(card)
 669       include 'DIMENSIONS'
 670       include 'COMMON.IOUNITS'
 671       character*(*) card
 672       character*80 karta,ucase
 673       external ilen
 674       read (inp,'(a)') karta
 675       karta=ucase(karta)
 676       card=' '
 677       do while (karta(80:80).eq.'&')
 678         card=card(:ilen(card)+1)//karta(:79)
 679         read (inp,'(a)') karta
 680         karta=ucase(karta)
 681       enddo
 682       card=card(:ilen(card)+1)//karta
 683       return
 684       end
 685 c----------------------------------------------------------------------------
 686       subroutine openunits
 687       implicit none
 688       include 'DIMENSIONS'
 689 #ifdef MPI
 690       include "mpif.h"
 691       character*3 liczba
 692       include "COMMON.MPI"
 693 #endif
 694       include 'COMMON.IOUNITS'
 695       include 'COMMON.CONTROL'
 696       integer lenpre,lenpot,ilen
 697       external ilen
 698       character*16 cformat,cprint
 699       character*16 ucase
 700       integer lenint,lenout
 701       call getenv('INPUT',prefix)
 702       call getenv('OUTPUT',prefout)
 703       call getenv('INTIN',prefintin)
 704       call getenv('COORD',cformat)
 705       call getenv('PRINTCOOR',cprint)
 706       call getenv('SCRATCHDIR',scratchdir)
 707       from_bx=.true.
 708       from_cx=.false.
 709       if (index(ucase(cformat),'CX').gt.0) then
 710         from_cx=.true.
 711         from_bx=.false.
 712       endif
 713       from_cart=.true.
 714       lenpre=ilen(prefix)
 715       lenout=ilen(prefout)
 716       lenint=ilen(prefintin)
 717 C Get the names and open the input files
 718       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
 719 #ifdef MPI
 720       write (liczba,'(bz,i3.3)') me
 721       outname=prefout(:lenout)//'_clust.out_'//liczba
 722 #else
 723       outname=prefout(:lenout)//'_clust.out'
 724 #endif
 725       if (from_bx) then
 726         intinname=prefintin(:lenint)//'.bx'
 727       else if (from_cx) then
 728         intinname=prefintin(:lenint)//'.cx'
 729       else
 730         intinname=prefintin(:lenint)//'.int'
 731       endif
 732       rmsname=prefintin(:lenint)//'.rms'
 733       open (jplot,file=prefout(:ilen(prefout))//'.tex',
 734      &       status='unknown')
 735       open (jrms,file=rmsname,status='unknown')
 736       open(iout,file=outname,status='unknown')
 737 C Get parameter filenames and open the parameter files.
 738       call getenv('BONDPAR',bondname)
 739       open (ibond,file=bondname,status='old')
 740       call getenv('THETPAR',thetname)
 741       open (ithep,file=thetname,status='old')
 742       call getenv('ROTPAR',rotname)
 743       open (irotam,file=rotname,status='old')
 744       call getenv('TORPAR',torname)
 745       open (itorp,file=torname,status='old')
 746       call getenv('TORDPAR',tordname)
 747       open (itordp,file=tordname,status='old')
 748       call getenv('FOURIER',fouriername)
 749       open (ifourier,file=fouriername,status='old')
 750       call getenv('ELEPAR',elename)
 751       open (ielep,file=elename,status='old')
 752       call getenv('SIDEPAR',sidename)
 753       open (isidep,file=sidename,status='old')
 754       call getenv('SIDEP',sidepname)
 755       open (isidep1,file=sidepname,status="old")
 756       call getenv('SCCORPAR',sccorname)
 757       open (isccor,file=sccorname,status="old")
 758 #ifndef OLDSCP
 759 C
 760 C 8/9/01 In the newest version SCp interaction constants are read from a file
 761 C Use -DOLDSCP to use hard-coded constants instead.
 762 C
 763       call getenv('SCPPAR',scpname)
 764       open (iscpp,file=scpname,status='old')
 765 #endif
 766       return
 767       end
 768 c-------------------------------------------------------------------------------
 769       subroutine read_dist_constr
 770       implicit real*8 (a-h,o-z)
 771       include 'DIMENSIONS'
 772 #ifdef MPI
 773       include 'mpif.h'
 774 #endif
 775       include 'COMMON.CONTROL'
 776       include 'COMMON.CHAIN'
 777       include 'COMMON.IOUNITS'
 778       include 'COMMON.SBRIDGE'
 779       integer ifrag_(2,100),ipair_(2,100)
 780       double precision wfrag_(100),wpair_(100)
 781       character*500 controlcard
 782       logical lprn /.true./
 783       logical normalize,next
 784       integer restr_type
 785       double precision xlink(4,0:4) /
 786 c           a          b       c     sigma
 787      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
 788      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
 789      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
 790      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
 791      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
 792       write (iout,*) "Calling read_dist_constr"
 793 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 794 c      call flush(iout)
 795       next=.true.
 796
 797       DO WHILE (next)
 798
 799       call card_concat(controlcard)
 800       next = index(controlcard,"NEXT").gt.0
 801       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
 802       write (iout,*) "restr_type",restr_type
 803       call readi(controlcard,"NFRAG",nfrag_,0)
 804       call readi(controlcard,"NFRAG",nfrag_,0)
 805       call readi(controlcard,"NPAIR",npair_,0)
 806       call readi(controlcard,"NDIST",ndist_,0)
 807       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 808       if (restr_type.eq.10)
 809      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
 810       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 811       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 812       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 813       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 814       normalize = index(controlcard,"NORMALIZE").gt.0
 815       write (iout,*) "WBOLTZD",wboltzd
 816 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 817 c      write (iout,*) "IFRAG"
 818 c      do i=1,nfrag_
 819 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 820 c      enddo
 821 c      write (iout,*) "IPAIR"
 822 c      do i=1,npair_
 823 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 824 c      enddo
 825       if (nfrag_.gt.0) then
 826         nres0=nres
 827         read(inp,'(a)') pdbfile
 828         write (iout,*)
 829      & "Distance restraints will be constructed from structure ",pdbfile
 830         open(ipdbin,file=pdbfile,status='old',err=11)
 831         call readpdb(.true.)
 832         nres=nres0
 833         close(ipdbin)
 834       endif
 835       do i=1,nfrag_
 836         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 837         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 838      &    ifrag_(2,i)=nstart_sup+nsup-1
 839 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 840 c        call flush(iout)
 841         if (wfrag_(i).eq.0.0d0) cycle
 842         do j=ifrag_(1,i),ifrag_(2,i)-1
 843           do k=j+1,ifrag_(2,i)
 844 c            write (iout,*) "j",j," k",k
 845             ddjk=dist(j,k)
 846             if (restr_type.eq.1) then
 847               nhpb=nhpb+1
 848               irestr_type(nhpb)=1
 849               ihpb(nhpb)=j
 850               jhpb(nhpb)=k
 851               dhpb(nhpb)=ddjk
 852               forcon(nhpb)=wfrag_(i)
 853             else if (constr_dist.eq.2) then
 854               if (ddjk.le.dist_cut) then
 855                 nhpb=nhpb+1
 856                 irestr_type(nhpb)=1
 857                 ihpb(nhpb)=j
 858                 jhpb(nhpb)=k
 859                 dhpb(nhpb)=ddjk
 860                 forcon(nhpb)=wfrag_(i)
 861               endif
 862             else if (restr_type.eq.3) then
 863               nhpb=nhpb+1
 864               irestr_type(nhpb)=1
 865               ihpb(nhpb)=j
 866               jhpb(nhpb)=k
 867               dhpb(nhpb)=ddjk
 868               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 869             endif
 870             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
 871      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 872           enddo
 873         enddo
 874       enddo
 875       do i=1,npair_
 876         if (wpair_(i).eq.0.0d0) cycle
 877         ii = ipair_(1,i)
 878         jj = ipair_(2,i)
 879         if (ii.gt.jj) then
 880           itemp=ii
 881           ii=jj
 882           jj=itemp
 883         endif
 884         do j=ifrag_(1,ii),ifrag_(2,ii)
 885           do k=ifrag_(1,jj),ifrag_(2,jj)
 886             if (restr_type.eq.1) then
 887               nhpb=nhpb+1
 888               irestr_type(nhpb)=1
 889               ihpb(nhpb)=j
 890               jhpb(nhpb)=k
 891               dhpb(nhpb)=ddjk
 892               forcon(nhpb)=wfrag_(i)
 893             else if (constr_dist.eq.2) then
 894               if (ddjk.le.dist_cut) then
 895                 nhpb=nhpb+1
 896                 irestr_type(nhpb)=1
 897                 ihpb(nhpb)=j
 898                 jhpb(nhpb)=k
 899                 dhpb(nhpb)=ddjk
 900                 forcon(nhpb)=wfrag_(i)
 901               endif
 902             else if (restr_type.eq.3) then
 903               nhpb=nhpb+1
 904               irestr_type(nhpb)=1
 905               ihpb(nhpb)=j
 906               jhpb(nhpb)=k
 907               dhpb(nhpb)=ddjk
 908               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 909             endif
 910             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
 911      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 912           enddo
 913         enddo
 914       enddo
 915
 916 c      print *,ndist_
 917       write (iout,*) "Distance restraints as read from input"
 918       do i=1,ndist_
 919         if (restr_type.eq.11) then
 920           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
 921      &     dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
 922 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
 923           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
 924           nhpb=nhpb+1
 925           irestr_type(nhpb)=11
 926           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
 927      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
 928      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
 929           if (ibecarb(nhpb).gt.0) then
 930             ihpb(nhpb)=ihpb(nhpb)+nres
 931             jhpb(nhpb)=jhpb(nhpb)+nres
 932           endif
 933         else if (constr_dist.eq.10) then
 934 c Cross-lonk Markov-like potential
 935           call card_concat(controlcard)
 936           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
 937           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
 938           ibecarb(nhpb+1)=0
 939           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
 940           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
 941           if (index(controlcard,"ZL").gt.0) then
 942             link_type=1
 943           else if (index(controlcard,"ADH").gt.0) then
 944             link_type=2
 945           else if (index(controlcard,"PDH").gt.0) then
 946             link_type=3
 947           else if (index(controlcard,"DSS").gt.0) then
 948             link_type=4
 949           else
 950             link_type=0
 951           endif
 952           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
 953      &       xlink(1,link_type))
 954           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
 955      &       xlink(2,link_type))
 956           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
 957      &       xlink(3,link_type))
 958           call reada(controlcard,"SIGMA",forcon(nhpb+1),
 959      &       xlink(4,link_type))
 960           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
 961 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
 962 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
 963           if (forcon(nhpb+1).le.0.0d0 .or.
 964      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
 965           nhpb=nhpb+1
 966           irestr_type(nhpb)=10
 967           if (ibecarb(nhpb).gt.0) then
 968             ihpb(nhpb)=ihpb(nhpb)+nres
 969             jhpb(nhpb)=jhpb(nhpb)+nres
 970           endif
 971           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
 972      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
 973      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
 974      &     irestr_type(nhpb)
 975         else
 976 C        print *,"in else"
 977           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
 978      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
 979           if (forcon(nhpb+1).gt.0.0d0) then
 980           nhpb=nhpb+1
 981           if (dhpb1(nhpb).eq.0.0d0) then
 982             irestr_type(nhpb)=1
 983           else
 984             irestr_type(nhpb)=2
 985           endif
 986           if (ibecarb(nhpb).gt.0) then
 987             ihpb(nhpb)=ihpb(nhpb)+nres
 988             jhpb(nhpb)=jhpb(nhpb)+nres
 989           endif
 990           if (dhpb(nhpb).eq.0.0d0)
 991      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 992         endif
 993         write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
 994      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
 995         endif
 996 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
 997 C        if (forcon(nhpb+1).gt.0.0d0) then
 998 C          nhpb=nhpb+1
 999 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1000       enddo
1001
1002       ENDDO ! next
1003
1004       fordepthmax=0.0d0
1005       if (normalize) then
1006         do i=nss+1,nhpb
1007           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
1008      &      fordepthmax=fordepth(i)
1009         enddo
1010         do i=nss+1,nhpb
1011           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
1012         enddo
1013       endif
1014       if (nhpb.gt.nss)  then
1015         write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
1016      &  "The following",nhpb-nss,
1017      &  " distance restraints have been imposed:",
1018      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
1019      &  "  score"," type"
1020         do i=nss+1,nhpb
1021           write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
1022      &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
1023      &  irestr_type(i)
1024         enddo
1025       endif
1026       call hpb_partition
1027       call flush(iout)
1028       return
1029    11 write (iout,*)"read_dist_restr: error reading reference structure"
1030       stop
1031       end
1032 c====-------------------------------------------------------------------
1033       subroutine read_constr_homology(lprn)
1034
1035       include 'DIMENSIONS'
1036 #ifdef MPI
1037       include 'mpif.h'
1038 #endif
1039       include 'COMMON.SETUP'
1040       include 'COMMON.CONTROL'
1041       include 'COMMON.CHAIN'
1042       include 'COMMON.IOUNITS'
1043       include 'COMMON.GEO'
1044       include 'COMMON.INTERACT'
1045       include 'COMMON.NAMES'
1046       include 'COMMON.HOMRESTR'
1047 c
1048 c For new homol impl
1049 c
1050       include 'COMMON.VAR'
1051 c     include 'include_unres/COMMON.VAR'
1052 c
1053
1054 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1055 c    &                 dist_cut
1056 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1057 c    &    sigma_odl_temp(maxres,maxres,max_template)
1058       character*2 kic2
1059       character*24 model_ki_dist, model_ki_angle
1060       character*500 controlcard
1061       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1062       integer idomain(max_template,maxres)
1063       integer ilen
1064       external ilen
1065       logical lprn
1066       logical unres_pdb,liiflag
1067 c
1068 c     FP - Nov. 2014 Temporary specifications for new vars
1069 c
1070       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
1071       double precision, dimension (max_template,maxres) :: rescore
1072       double precision, dimension (max_template,maxres) :: rescore2
1073       character*24 tpl_k_rescore
1074 c -----------------------------------------------------------------
1075 c Reading multiple PDB ref structures and calculation of retraints
1076 c not using pre-computed ones stored in files model_ki_{dist,angle}
1077 c FP (Nov., 2014)
1078 c -----------------------------------------------------------------
1079 c
1080 c
1081 c Alternative: reading from input
1082 #ifdef DEBUG
1083       write (iout,*) "BEGIN READ HOMOLOGY INFO"
1084 #ifdef AIX
1085       call flush_(iout)
1086 #else
1087       call flush(iout)
1088 #endif
1089 #endif
1090       call card_concat(controlcard)
1091       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1092       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1093       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1094       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1095       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1096       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1097       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1098       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1099       if (homol_nset.gt.1)then
1100          call readi(controlcard,"ISET",iset,1)
1101          call card_concat(controlcard)
1102          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1103       else
1104         iset=1
1105         waga_homology(1)=1.0
1106       endif
1107 c
1108 #ifdef DEBUG
1109       write(iout,*) "read_constr_homology iset",iset
1110       write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
1111 #ifdef AIX
1112       call flush_(iout)
1113 #else
1114       call flush(iout)
1115 #endif
1116 #endif
1117 cd      write (iout,*) "nnt",nnt," nct",nct
1118 cd      call flush(iout)
1119
1120
1121       lim_odl=0
1122       lim_dih=0
1123 c
1124 c  Reading HM global scores (prob not required)
1125 c
1126       do i = nnt,nct
1127         do k=1,constr_homology
1128          idomain(k,i)=0
1129         enddo
1130       enddo
1131 c     open (4,file="HMscore")
1132 c     do k=1,constr_homology
1133 c       read (4,*,end=521) hmscore_tmp
1134 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
1135 c       write(*,*) "Model", k, ":", hmscore(k)
1136 c     enddo
1137 c521  continue
1138
1139       ii=0
1140       do i = nnt,nct-2
1141         do j=i+2,nct
1142         ii=ii+1
1143         ii_in_use(ii)=0
1144         enddo
1145       enddo
1146 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1147
1148       if (read_homol_frag) then
1149         call read_klapaucjusz
1150       else
1151       write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1152       do k=1,constr_homology
1153
1154         read(inp,'(a)') pdbfile
1155 c        write (iout,*) "k ",k," pdbfile ",pdbfile
1156 c  Next stament causes error upon compilation (?)
1157 c       if(me.eq.king.or. .not. out1file)
1158 c         write (iout,'(2a)') 'PDB data will be read from file ',
1159 c    &   pdbfile(:ilen(pdbfile))
1160          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1161      &  pdbfile(:ilen(pdbfile))
1162         open(ipdbin,file=pdbfile,status='old',err=33)
1163         goto 34
1164   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
1165      &  pdbfile(:ilen(pdbfile))
1166         stop
1167   34    continue
1168 c        print *,'Begin reading pdb data'
1169 c
1170 c Files containing res sim or local scores (former containing sigmas)
1171 c
1172
1173         write(kic2,'(bz,i2.2)') k
1174
1175         tpl_k_rescore="template"//kic2//".sco"
1176
1177         unres_pdb=.false.
1178         call readpdb
1179         do i=1,2*nres
1180           do j=1,3
1181             crefjlee(j,i)=c(j,i)
1182           enddo
1183         enddo
1184 #ifdef DEBUG
1185         do i=1,nres
1186           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1187      &      (crefjlee(j,i+nres),j=1,3)
1188         enddo
1189         write (iout,*) "READ HOMOLOGY INFO"
1190         write (iout,*) "read_constr_homology x: after reading pdb file"
1191         write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1192         write (iout,*) "waga_dist",waga_dist
1193         write (iout,*) "waga_angle",waga_angle
1194         write (iout,*) "waga_theta",waga_theta
1195         write (iout,*) "waga_d",waga_d
1196         write (iout,*) "dist_cut",dist_cut
1197 #endif
1198 #ifdef AIX
1199       call flush_(iout)
1200 #else
1201       call flush(iout)
1202 #endif
1203
1204 c
1205 c     Distance restraints
1206 c
1207 c          ... --> odl(k,ii)
1208 C Copy the coordinates from reference coordinates (?)
1209         do i=1,2*nres
1210           do j=1,3
1211             c(j,i)=cref(j,i)
1212 c           write (iout,*) "c(",j,i,") =",c(j,i)
1213           enddo
1214         enddo
1215 c
1216 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1217 c
1218 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1219           open (ientin,file=tpl_k_rescore,status='old')
1220           if (nnt.gt.1) rescore(k,1)=0.0d0
1221           do irec=nnt,nct ! loop for reading res sim
1222             if (read2sigma) then
1223              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1224      &                                idomain_tmp
1225              i_tmp=i_tmp+nnt-1
1226              idomain(k,i_tmp)=idomain_tmp
1227              rescore(k,i_tmp)=rescore_tmp
1228              rescore2(k,i_tmp)=rescore2_tmp
1229              write(iout,'(a7,i5,2f10.5,i5)') "rescore",
1230      &                      i_tmp,rescore2_tmp,rescore_tmp,
1231      &                                idomain_tmp
1232             else
1233              idomain(k,irec)=1
1234              read (ientin,*,end=1401) rescore_tmp
1235
1236 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1237              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1238 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1239             endif
1240           enddo
1241  1401   continue
1242         close (ientin)
1243         if (waga_dist.ne.0.0d0) then
1244           ii=0
1245           do i = nnt,nct-2
1246             do j=i+2,nct
1247
1248               x12=c(1,i)-c(1,j)
1249               y12=c(2,i)-c(2,j)
1250               z12=c(3,i)-c(3,j)
1251               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1252 c              write (iout,*) k,i,j,distal,dist2_cut
1253
1254             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1255      &            .and. distal.le.dist2_cut ) then
1256
1257               ii=ii+1
1258               ii_in_use(ii)=1
1259               l_homo(k,ii)=.true.
1260
1261 c             write (iout,*) "k",k
1262 c             write (iout,*) "i",i," j",j," constr_homology",
1263 c    &                       constr_homology
1264               ires_homo(ii)=i
1265               jres_homo(ii)=j
1266               odl(k,ii)=distal
1267               if (read2sigma) then
1268                 sigma_odl(k,ii)=0
1269                 do ik=i,j
1270                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1271                 enddo
1272                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1273                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1274      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1275               else
1276                 if (odl(k,ii).le.dist_cut) then
1277                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1278                 else
1279 #ifdef OLDSIGMA
1280                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1281      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1282 #else
1283                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1284      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1285 #endif
1286                 endif
1287               endif
1288               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1289             else
1290               ii=ii+1
1291               l_homo(k,ii)=.false.
1292             endif
1293             enddo
1294           enddo
1295         lim_odl=ii
1296         endif
1297 c
1298 c     Theta, dihedral and SC retraints
1299 c
1300         if (waga_angle.gt.0.0d0) then
1301 c         open (ientin,file=tpl_k_sigma_dih,status='old')
1302 c         do irec=1,maxres-3 ! loop for reading sigma_dih
1303 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1304 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1305 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1306 c    &                            sigma_dih(k,i+nnt-1)
1307 c         enddo
1308 c1402   continue
1309 c         close (ientin)
1310           do i = nnt+3,nct
1311             if (idomain(k,i).eq.0) then
1312                sigma_dih(k,i)=0.0
1313                cycle
1314             endif
1315             dih(k,i)=phiref(i) ! right?
1316 c           read (ientin,*) sigma_dih(k,i) ! original variant
1317 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
1318 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1319 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1320 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
1321 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
1322
1323             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1324      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1325 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1326 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
1327 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1328 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1329 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
1330 c   Instead of res sim other local measure of b/b str reliability possible
1331             if (sigma_dih(k,i).ne.0)
1332      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1333 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1334           enddo
1335           lim_dih=nct-nnt-2
1336         endif
1337
1338         if (waga_theta.gt.0.0d0) then
1339 c         open (ientin,file=tpl_k_sigma_theta,status='old')
1340 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1341 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1342 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1343 c    &                              sigma_theta(k,i+nnt-1)
1344 c         enddo
1345 c1403   continue
1346 c         close (ientin)
1347
1348           do i = nnt+2,nct ! right? without parallel.
1349 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
1350 c         do i=ithet_start,ithet_end ! with FG parallel.
1351              if (idomain(k,i).eq.0) then
1352               sigma_theta(k,i)=0.0
1353               cycle
1354              endif
1355              thetatpl(k,i)=thetaref(i)
1356 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1357 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
1358 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1359 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
1360 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
1361              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1362      &                        rescore(k,i-2))/3.0
1363 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1364              if (sigma_theta(k,i).ne.0)
1365      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1366
1367 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1368 c                             rescore(k,i-2) !  right expression ?
1369 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1370           enddo
1371         endif
1372
1373         if (waga_d.gt.0.0d0) then
1374 c       open (ientin,file=tpl_k_sigma_d,status='old')
1375 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1376 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1377 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1378 c    &                          sigma_d(k,i+nnt-1)
1379 c         enddo
1380 c1404   continue
1381
1382           do i = nnt,nct ! right? without parallel.
1383 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
1384 c         do i=loc_start,loc_end ! with FG parallel.
1385                if (itype(i).eq.10) cycle
1386                if (idomain(k,i).eq.0 ) then
1387                   sigma_d(k,i)=0.0
1388                   cycle
1389                endif
1390                xxtpl(k,i)=xxref(i)
1391                yytpl(k,i)=yyref(i)
1392                zztpl(k,i)=zzref(i)
1393 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1394 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1395 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1396 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
1397                sigma_d(k,i)=rescore(k,i) !  right expression ?
1398                if (sigma_d(k,i).ne.0)
1399      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1400
1401 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
1402 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1403 c              read (ientin,*) sigma_d(k,i) ! 1st variant
1404                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
1405           enddo
1406         endif
1407       enddo
1408 c
1409 c remove distance restraints not used in any model from the list
1410 c shift data in all arrays
1411 c
1412       if (waga_dist.ne.0.0d0) then
1413         ii=0
1414         liiflag=.true.
1415         do i=nnt,nct-2
1416          do j=i+2,nct
1417           ii=ii+1
1418           if (ii_in_use(ii).eq.0.and.liiflag) then
1419             liiflag=.false.
1420             iistart=ii
1421           endif
1422           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1423      &                   .not.liiflag.and.ii.eq.lim_odl) then
1424              if (ii.eq.lim_odl) then
1425               iishift=ii-iistart+1
1426              else
1427               iishift=ii-iistart
1428              endif
1429              liiflag=.true.
1430              do ki=iistart,lim_odl-iishift
1431               ires_homo(ki)=ires_homo(ki+iishift)
1432               jres_homo(ki)=jres_homo(ki+iishift)
1433               ii_in_use(ki)=ii_in_use(ki+iishift)
1434               do k=1,constr_homology
1435                odl(k,ki)=odl(k,ki+iishift)
1436                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1437                l_homo(k,ki)=l_homo(k,ki+iishift)
1438               enddo
1439              enddo
1440              ii=ii-iishift
1441              lim_odl=lim_odl-iishift
1442           endif
1443          enddo
1444         enddo
1445       endif
1446
1447       endif ! .not. klapaucjusz
1448
1449       if (constr_homology.gt.0) call homology_partition
1450       if (constr_homology.gt.0) call init_int_table
1451 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1452 cd     & "lim_xx=",lim_xx
1453 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1454 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1455 c
1456 c Print restraints
1457 c
1458       if (.not.lprn) return
1459 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1460       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1461        write (iout,*) "Distance restraints from templates"
1462        do ii=1,lim_odl
1463        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1464      &  ii,ires_homo(ii),jres_homo(ii),
1465      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1466      &  ki=1,constr_homology)
1467        enddo
1468        write (iout,*) "Dihedral angle restraints from templates"
1469        do i=nnt+3,nct
1470         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1471      &      (rad2deg*dih(ki,i),
1472      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1473        enddo
1474        write (iout,*) "Virtual-bond angle restraints from templates"
1475        do i=nnt+2,nct
1476         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1477      &      (rad2deg*thetatpl(ki,i),
1478      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1479        enddo
1480        write (iout,*) "SC restraints from templates"
1481        do i=nnt,nct
1482         write(iout,'(i5,100(4f8.2,4x))') i,
1483      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1484      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1485        enddo
1486       endif
1487 c -----------------------------------------------------------------
1488       return
1489       end
1490 c----------------------------------------------------------------------
1491       subroutine read_klapaucjusz
1492
1493       include 'DIMENSIONS'
1494 #ifdef MPI
1495       include 'mpif.h'
1496 #endif
1497       include 'COMMON.SETUP'
1498       include 'COMMON.CONTROL'
1499       include 'COMMON.CHAIN'
1500       include 'COMMON.IOUNITS'
1501       include 'COMMON.GEO'
1502       include 'COMMON.INTERACT'
1503       include 'COMMON.NAMES'
1504       include 'COMMON.HOMRESTR'
1505       character*256 fragfile
1506       integer ninclust(maxclust),inclust(max_template,maxclust),
1507      &  nresclust(maxclust),iresclust(maxres,maxclust)
1508
1509       character*2 kic2
1510       character*24 model_ki_dist, model_ki_angle
1511       character*500 controlcard
1512       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1513       integer idomain(max_template,maxres)
1514       logical lprn /.true./
1515       integer ilen
1516       external ilen
1517       logical unres_pdb,liiflag
1518 c
1519 c
1520       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
1521       double precision, dimension (max_template,maxres) :: rescore
1522       double precision, dimension (max_template,maxres) :: rescore2
1523       character*24 tpl_k_rescore
1524
1525 c
1526 c For new homol impl
1527 c
1528       include 'COMMON.VAR'
1529 c
1530       double precision chomo(3,maxres2+2,max_template)
1531       call getenv("FRAGFILE",fragfile)
1532       open(ientin,file=fragfile,status="old",err=10)
1533       read(ientin,*) constr_homology,nclust
1534       l_homo = .false.
1535       sigma_theta=0.0
1536       sigma_d=0.0
1537       sigma_dih=0.0
1538 c Read pdb files
1539       do k=1,constr_homology
1540         read(ientin,'(a)') pdbfile
1541         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
1542      &  pdbfile(:ilen(pdbfile))
1543         open(ipdbin,file=pdbfile,status='old',err=33)
1544         goto 34
1545   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
1546      &  pdbfile(:ilen(pdbfile))
1547         stop
1548   34    continue
1549         unres_pdb=.false.
1550         call readpdb
1551         do i=1,2*nres
1552           do j=1,3
1553             chomo(j,i,k)=c(j,i)
1554           enddo
1555         enddo
1556         do i=1,nres
1557           rescore(k,i)=0.2d0
1558           rescore2(k,i)=1.0d0
1559         enddo
1560       enddo
1561 c Read clusters
1562       do i=1,nclust
1563         read(ientin,*) ninclust(i),nresclust(i)
1564         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
1565         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
1566       enddo
1567 c
1568 c Loop over clusters
1569 c
1570       do l=1,nclust
1571         do ll = 1,ninclust(l)
1572
1573         k = inclust(ll,l)
1574         do i=1,nres
1575           idomain(k,i)=0
1576         enddo
1577         do i=1,nresclust(l)
1578           if (nnt.gt.1)  then
1579             idomain(k,iresclust(i,l)+1) = 1
1580           else
1581             idomain(k,iresclust(i,l)) = 1
1582           endif
1583         enddo
1584 c
1585 c     Distance restraints
1586 c
1587 c          ... --> odl(k,ii)
1588 C Copy the coordinates from reference coordinates (?)
1589         do i=1,2*nres
1590           do j=1,3
1591             c(j,i)=chomo(j,i,k)
1592 c           write (iout,*) "c(",j,i,") =",c(j,i)
1593           enddo
1594         enddo
1595         call int_from_cart(.true.,.false.)
1596         call sc_loc_geom(.false.)
1597         do i=1,nres
1598           thetaref(i)=theta(i)
1599           phiref(i)=phi(i)
1600         enddo
1601         if (waga_dist.ne.0.0d0) then
1602           ii=0
1603           do i = nnt,nct-2
1604             do j=i+2,nct
1605
1606               x12=c(1,i)-c(1,j)
1607               y12=c(2,i)-c(2,j)
1608               z12=c(3,i)-c(3,j)
1609               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1610 c              write (iout,*) k,i,j,distal,dist2_cut
1611
1612             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1613      &            .and. distal.le.dist2_cut ) then
1614
1615               ii=ii+1
1616               ii_in_use(ii)=1
1617               l_homo(k,ii)=.true.
1618
1619 c             write (iout,*) "k",k
1620 c             write (iout,*) "i",i," j",j," constr_homology",
1621 c    &                       constr_homology
1622               ires_homo(ii)=i
1623               jres_homo(ii)=j
1624               odl(k,ii)=distal
1625               if (read2sigma) then
1626                 sigma_odl(k,ii)=0
1627                 do ik=i,j
1628                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1629                 enddo
1630                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1631                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1632      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1633               else
1634                 if (odl(k,ii).le.dist_cut) then
1635                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1636                 else
1637 #ifdef OLDSIGMA
1638                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1639      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1640 #else
1641                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1642      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1643 #endif
1644                 endif
1645               endif
1646               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1647             else
1648               ii=ii+1
1649 c              l_homo(k,ii)=.false.
1650             endif
1651             enddo
1652           enddo
1653         lim_odl=ii
1654         endif
1655 c
1656 c     Theta, dihedral and SC retraints
1657 c
1658         if (waga_angle.gt.0.0d0) then
1659           do i = nnt+3,nct
1660             if (idomain(k,i).eq.0) then
1661 c               sigma_dih(k,i)=0.0
1662                cycle
1663             endif
1664             dih(k,i)=phiref(i)
1665             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1666      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1667 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
1668 c     &       " sigma_dihed",sigma_dih(k,i)
1669             if (sigma_dih(k,i).ne.0)
1670      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1671           enddo
1672           lim_dih=nct-nnt-2
1673         endif
1674
1675         if (waga_theta.gt.0.0d0) then
1676           do i = nnt+2,nct
1677              if (idomain(k,i).eq.0) then
1678 c              sigma_theta(k,i)=0.0
1679               cycle
1680              endif
1681              thetatpl(k,i)=thetaref(i)
1682              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1683      &                        rescore(k,i-2))/3.0
1684              if (sigma_theta(k,i).ne.0)
1685      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1686           enddo
1687         endif
1688
1689         if (waga_d.gt.0.0d0) then
1690           do i = nnt,nct
1691                if (itype(i).eq.10) cycle
1692                if (idomain(k,i).eq.0 ) then
1693 c                  sigma_d(k,i)=0.0
1694                   cycle
1695                endif
1696                xxtpl(k,i)=xxref(i)
1697                yytpl(k,i)=yyref(i)
1698                zztpl(k,i)=zzref(i)
1699                sigma_d(k,i)=rescore(k,i)
1700                if (sigma_d(k,i).ne.0)
1701      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1702                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
1703           enddo
1704         endif
1705       enddo ! l
1706       enddo ! ll
1707 c
1708 c remove distance restraints not used in any model from the list
1709 c shift data in all arrays
1710 c
1711       if (waga_dist.ne.0.0d0) then
1712         ii=0
1713         liiflag=.true.
1714         do i=nnt,nct-2
1715          do j=i+2,nct
1716           ii=ii+1
1717           if (ii_in_use(ii).eq.0.and.liiflag) then
1718             liiflag=.false.
1719             iistart=ii
1720           endif
1721           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1722      &                   .not.liiflag.and.ii.eq.lim_odl) then
1723              if (ii.eq.lim_odl) then
1724               iishift=ii-iistart+1
1725              else
1726               iishift=ii-iistart
1727              endif
1728              liiflag=.true.
1729              do ki=iistart,lim_odl-iishift
1730               ires_homo(ki)=ires_homo(ki+iishift)
1731               jres_homo(ki)=jres_homo(ki+iishift)
1732               ii_in_use(ki)=ii_in_use(ki+iishift)
1733               do k=1,constr_homology
1734                odl(k,ki)=odl(k,ki+iishift)
1735                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1736                l_homo(k,ki)=l_homo(k,ki+iishift)
1737               enddo
1738              enddo
1739              ii=ii-iishift
1740              lim_odl=lim_odl-iishift
1741           endif
1742          enddo
1743         enddo
1744       endif
1745
1746       return
1747    10 stop "Error infragment file"
1748       end
1749 c----------------------------------------------------------------------