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
  80 #ifdef AIX
  81       call flush_(iout)
  82 #else
  83       call flush(iout)
  84 #endif
  85       if (min_var) iopt=1
  86       return
  87       end
  88 c--------------------------------------------------------------------------
  89       subroutine molread
  90 C
  91 C Read molecular data.
  92 C
  93       implicit none
  94       include 'DIMENSIONS'
  95       include 'COMMON.IOUNITS'
  96       include 'COMMON.GEO'
  97       include 'COMMON.VAR'
  98       include 'COMMON.INTERACT'
  99       include 'COMMON.LOCAL'
 100       include 'COMMON.NAMES'
 101       include 'COMMON.CHAIN'
 102       include 'COMMON.FFIELD'
 103       include 'COMMON.SBRIDGE'
 104       include 'COMMON.HEADER'
 105       include 'COMMON.CONTROL'
 106       include 'COMMON.CONTACTS'
 107       include 'COMMON.TIME1'
 108       include 'COMMON.TORCNSTR'
 109 #ifdef MPL
 110       include 'COMMON.INFO'
 111 #endif
 112       character*4 sequence(maxres)
 113       character*800 weightcard
 114       integer rescode
 115       double precision x(maxvar)
 116       integer itype_pdb(maxres)
 117       logical seq_comp
 118       integer i,j
 119       write (iout,*) " MOLREAD: NRES",NRES
 120 C
 121 C Body
 122 C
 123 C Read weights of the subsequent energy terms.
 124       call card_concat(weightcard)
 125       call reada(weightcard,'WSC',wsc,1.0d0)
 126       call reada(weightcard,'WLONG',wsc,wsc)
 127       call reada(weightcard,'WSCP',wscp,1.0d0)
 128       call reada(weightcard,'WELEC',welec,1.0D0)
 129       call reada(weightcard,'WVDWPP',wvdwpp,welec)
 130       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 131       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 132       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 133       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 134       call reada(weightcard,'WTURN3',wturn3,1.0D0)
 135       call reada(weightcard,'WTURN4',wturn4,1.0D0)
 136       call reada(weightcard,'WTURN6',wturn6,1.0D0)
 137       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 138       call reada(weightcard,'WBOND',wbond,1.0D0)
 139       call reada(weightcard,'WTOR',wtor,1.0D0)
 140       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 141       call reada(weightcard,'WANG',wang,1.0D0)
 142       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 143       call reada(weightcard,'SCAL14',scal14,0.4D0)
 144       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 145       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 146       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 147       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 148       call reada(weightcard,"D0CM",d0cm,3.78d0)
 149       call reada(weightcard,"AKCM",akcm,15.1d0)
 150       call reada(weightcard,"AKTH",akth,11.0d0)
 151       call reada(weightcard,"AKCT",akct,12.0d0)
 152       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 153       call reada(weightcard,"V2SS",v2ss,7.61d0)
 154       call reada(weightcard,"V3SS",v3ss,13.7d0)
 155       call reada(weightcard,"EBR",ebr,-5.50D0)
 156 C     Bartek
 157       call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
 158       call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
 159       call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
 160       call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
 161       if (index(weightcard,'SOFT').gt.0) ipot=6
 162 C 12/1/95 Added weight for the multi-body term WCORR
 163       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 164       do i=1,maxres-1
 165         do j=i+1,maxres
 166           dyn_ssbond_ij(i,j)=1.0d300
 167         enddo
 168       enddo
 169       call reada(weightcard,"HT",Ht,0.0D0)
 170       if (dyn_ss) then
 171         ss_depth=ebr/wsc-0.25*eps(1,1)
 172         Ht=Ht/wsc-0.25*eps(1,1)
 173         akcm=akcm*wstrain/wsc
 174         akth=akth*wstrain/wsc
 175         akct=akct*wstrain/wsc
 176         v1ss=v1ss*wstrain/wsc
 177         v2ss=v2ss*wstrain/wsc
 178         v3ss=v3ss*wstrain/wsc
 179       else
 180         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 181       endif
 182       write (iout,'(/a)') "Disulfide bridge parameters:"
 183       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 184         write (iout,'(a,f10.2)') 'S-S depth: ',ss_depth
 185       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 186       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 187       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 188      & ' v3ss:',v3ss
 189       write (iout,'(2(a,f10.2))') 'ht:',ht,' eps:', eps(1,1)
 190       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 191       weights(1)=wsc
 192       weights(2)=wscp
 193       weights(3)=welec
 194       weights(4)=wcorr
 195       weights(5)=wcorr5
 196       weights(6)=wcorr6
 197       weights(7)=wel_loc
 198       weights(8)=wturn3
 199       weights(9)=wturn4
 200       weights(10)=wturn6
 201       weights(11)=wang
 202       weights(12)=wscloc
 203       weights(13)=wtor
 204       weights(14)=wtor_d
 205       weights(15)=wstrain
 206       weights(16)=wvdwpp
 207       weights(17)=wbond
 208       weights(18)=scal14
 209       weights(22)=wdfa_dist
 210       weights(23)=wdfa_tor
 211       weights(24)=wdfa_nei
 212       weights(25)=wdfa_beta
 213       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 214      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
 215      &  wturn4,wturn6,wsccor,wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
 216    10 format (/'Energy-term weights (unscaled):'//
 217      & 'WSCC=   ',f10.6,' (SC-SC)'/
 218      & 'WSCP=   ',f10.6,' (SC-p)'/
 219      & 'WELEC=  ',f10.6,' (p-p electr)'/
 220      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 221      & 'WBOND=  ',f10.6,' (stretching)'/
 222      & 'WANG=   ',f10.6,' (bending)'/
 223      & 'WSCLOC= ',f10.6,' (SC local)'/
 224      & 'WTOR=   ',f10.6,' (torsional)'/
 225      & 'WTORD=  ',f10.6,' (double torsional)'/
 226      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 227      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 228      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 229      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 230      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 231      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 232      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 233      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 234      & 'WSCCOR= ',f10.6,' (SC-backbone torsional correalations)'/
 235      & 'WDFAD=  ',f10.6,' (DFA distance)'/
 236      & 'WDFAT=  ',f10.6,' (DFA torsional)'/
 237      & 'WDFAN=  ',f10.6,' (DFA neighbors)'/
 238      & 'WDFAB=  ',f10.6,' (DFA beta)'/)
 239       if (wcorr4.gt.0.0d0) then
 240         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 241      &   'between contact pairs of peptide groups'
 242         write (iout,'(2(a,f5.3/))')
 243      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 244      &  'Range of quenching the correlation terms:',2*delt_corr
 245       else if (wcorr.gt.0.0d0) then
 246         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 247      &   'between contact pairs of peptide groups'
 248       endif
 249       write (iout,'(a,f8.3)')
 250      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 251       write (iout,'(a,f8.3)')
 252      &  'General scaling factor of SC-p interactions:',scalscp
 253       r0_corr=cutoff_corr-delt_corr
 254       do i=1,20
 255         aad(i,1)=scalscp*aad(i,1)
 256         aad(i,2)=scalscp*aad(i,2)
 257         bad(i,1)=scalscp*bad(i,1)
 258         bad(i,2)=scalscp*bad(i,2)
 259       enddo
 260
 261 #ifdef AIX
 262       call flush_(iout)
 263 #else
 264       call flush(iout)
 265 #endif
 266       print *,'indpdb=',indpdb,' pdbref=',pdbref
 267
 268 C Read sequence if not taken from the pdb file.
 269       if (iscode.gt.0) then
 270         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 271       else
 272         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 273       endif
 274 C Convert sequence to numeric code
 275       do i=1,nres
 276         itype(i)=rescode(i,sequence(i),iscode)
 277       enddo
 278       if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
 279         write (iout,*)
 280      &   "Glycine is the first full residue, initial dummy deleted"
 281         do i=1,nres
 282           itype(i)=itype(i+1)
 283         enddo
 284         nres=nres-1
 285       endif
 286       if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
 287         write (iout,*)
 288      &   "Glycine is the last full residue, terminal dummy deleted"
 289         nres=nres-1
 290       endif
 291       print *,nres
 292       print '(20i4)',(itype(i),i=1,nres)
 293
 294       do i=1,nres
 295 #ifdef PROCOR
 296         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 297 #else
 298         if (itype(i).eq.21) then
 299 #endif
 300           itel(i)=0
 301 #ifdef PROCOR
 302         else if (itype(i+1).ne.20) then
 303 #else
 304         else if (itype(i).ne.20) then
 305 #endif
 306           itel(i)=1
 307         else
 308           itel(i)=2
 309         endif
 310       enddo
 311       write (iout,*) "ITEL"
 312       do i=1,nres-1
 313         write (iout,*) i,itype(i),itel(i)
 314       enddo
 315
 316       print *,'Call Read_Bridge.'
 317       call read_bridge
 318       if (with_dihed_constr) then
 319
 320       read (inp,*) ndih_constr
 321       if (ndih_constr.gt.0) then
 322         read (inp,*) ftors
 323         write (iout,*) 'FTORS',ftors
 324         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 325         write (iout,*)
 326      &   'There are',ndih_constr,' constraints on phi angles.'
 327         do i=1,ndih_constr
 328           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 329         enddo
 330         do i=1,ndih_constr
 331           phi0(i)=deg2rad*phi0(i)
 332           drange(i)=deg2rad*drange(i)
 333         enddo
 334       endif
 335
 336       endif
 337
 338       nnt=1
 339       nct=nres
 340       print *,'NNT=',NNT,' NCT=',NCT
 341       if (itype(1).eq.21) nnt=2
 342       if (itype(nres).eq.21) nct=nct-1
 343       if (nstart.lt.nnt) nstart=nnt
 344       if (nend.gt.nct .or. nend.eq.0) nend=nct
 345       write (iout,*) "nstart",nstart," nend",nend
 346       nres0=nres
 347
 348 C     Juyong:READ init_vars
 349 C     Initialize variables!
 350 C     Juyong:READ read_info
 351 C     READ fragment information!!
 352 C     both routines should be in dfa.F file!!
 353
 354       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
 355      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 356 #ifdef DEBUG
 357        write (iout,*) "Calling init_dfa_vars"
 358 #ifdef AIX
 359       call flush_(iout)
 360 #else
 361       call flush(iout)
 362 #endif
 363 #endif
 364        call init_dfa_vars
 365 #ifdef DEBUG
 366        write (iout,*) 'init_dfa_vars finished!'
 367 #ifdef AIX
 368       call flush_(iout)
 369 #else
 370       call flush(iout)
 371 #endif
 372 #endif
 373        call read_dfa_info
 374 #ifdef DEBUG
 375        write (iout,*) 'read_dfa_info finished!'
 376 #ifdef AIX
 377       call flush_(iout)
 378 #else
 379       call flush(iout)
 380 #endif
 381 #endif
 382       endif
 383
 384       if (constr_homology.gt.0) then
 385         call read_constr_homology(print_homology_restraints)
 386       endif
 387
 388 c      if (pdbref) then
 389 c        read(inp,'(a)') pdbfile
 390 c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
 391 c        open(ipdbin,file=pdbfile,status='old',err=33)
 392 c        goto 34
 393 c  33    write (iout,'(a)') 'Error opening PDB file.'
 394 c        stop
 395 c  34    continue
 396 c        print *,'Begin reading pdb data'
 397 c        call readpdb
 398 c        print *,'Finished reading pdb data'
 399 c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
 400 c        do i=1,nres
 401 c          itype_pdb(i)=itype(i)
 402 c        enddo
 403 c        close (ipdbin)
 404 c        write (iout,'(a,i3)') 'nsup=',nsup
 405 c        nstart_seq=nnt
 406 c        if (nsup.le.(nct-nnt+1)) then
 407 c          do i=0,nct-nnt+1-nsup
 408 c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 409 c              nstart_seq=nnt+i
 410 c              goto 111
 411 c            endif
 412 c          enddo
 413 c          write (iout,'(a)')
 414 c     &            'Error - sequences to be superposed do not match.'
 415 c          stop
 416 c        else
 417 c          do i=0,nsup-(nct-nnt+1)
 418 c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 419 c     &      then
 420 c              nstart_sup=nstart_sup+i
 421 c              nsup=nct-nnt+1
 422 c              goto 111
 423 c            endif
 424 c          enddo
 425 c          write (iout,'(a)')
 426 c     &            'Error - sequences to be superposed do not match.'
 427 c        endif
 428 c  111   continue
 429 c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 430 c     &                 ' nstart_seq=',nstart_seq
 431 c      endif
 432       call init_int_table
 433       call setup_var
 434       write (iout,*) "molread: REFSTR",refstr
 435       if (refstr) then
 436         if (.not.pdbref) then
 437           call read_angles(inp,*38)
 438           goto 39
 439    38     write (iout,'(a)') 'Error reading reference structure.'
 440 #ifdef MPL
 441           call mp_stopall(Error_Msg)
 442 #else
 443           stop 'Error reading reference structure'
 444 #endif
 445    39     call chainbuild
 446           nstart_sup=nnt
 447           nstart_seq=nnt
 448           nsup=nct-nnt+1
 449           do i=1,2*nres
 450             do j=1,3
 451               cref(j,i)=c(j,i)
 452             enddo
 453           enddo
 454         endif
 455         call contact(print_contact_map,ncont_ref,icont_ref)
 456       endif
 457 c Read distance restraints
 458       if (constr_dist.gt.0) then
 459         call read_dist_constr
 460         call hpb_partition
 461       endif
 462       return
 463       end
 464 c-----------------------------------------------------------------------------
 465       logical function seq_comp(itypea,itypeb,length)
 466       implicit none
 467       integer length,itypea(length),itypeb(length)
 468       integer i
 469       do i=1,length
 470         if (itypea(i).ne.itypeb(i)) then
 471           seq_comp=.false.
 472           return
 473         endif
 474       enddo
 475       seq_comp=.true.
 476       return
 477       end
 478 c-----------------------------------------------------------------------------
 479       subroutine read_bridge
 480 C Read information about disulfide bridges.
 481       implicit none
 482       include 'DIMENSIONS'
 483       include 'COMMON.IOUNITS'
 484       include 'COMMON.GEO'
 485       include 'COMMON.VAR'
 486       include 'COMMON.INTERACT'
 487       include 'COMMON.LOCAL'
 488       include 'COMMON.NAMES'
 489       include 'COMMON.CHAIN'
 490       include 'COMMON.FFIELD'
 491       include 'COMMON.SBRIDGE'
 492       include 'COMMON.HEADER'
 493       include 'COMMON.CONTROL'
 494       include 'COMMON.TIME1'
 495 #ifdef MPL
 496       include 'COMMON.INFO'
 497 #endif
 498       integer i,j
 499 C Read bridging residues.
 500       read (inp,*) ns,(iss(i),i=1,ns)
 501       write(iout,*)'ns=',ns
 502 C Check whether the specified bridging residues are cystines.
 503       do i=1,ns
 504         if (itype(iss(i)).ne.1) then
 505           write (iout,'(2a,i3,a)')
 506      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 507      &   ' can form a disulfide bridge?!!!'
 508           write (*,'(2a,i3,a)')
 509      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 510      &   ' can form a disulfide bridge?!!!'
 511 #ifdef MPL
 512          call mp_stopall(error_msg)
 513 #else
 514          stop
 515 #endif
 516         endif
 517       enddo
 518 C Read preformed bridges.
 519       if (ns.gt.0) then
 520       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 521       if (nss.gt.0) then
 522         nhpb=nss
 523 C Check if the residues involved in bridges are in the specified list of
 524 C bridging residues.
 525         do i=1,nss
 526           do j=1,i-1
 527             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 528      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 529               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 530      &      ' contains residues present in other pairs.'
 531               write (*,'(a,i3,a)') 'Disulfide pair',i,
 532      &      ' contains residues present in other pairs.'
 533 #ifdef MPL
 534               call mp_stopall(error_msg)
 535 #else
 536               stop
 537 #endif
 538             endif
 539           enddo
 540           do j=1,ns
 541             if (ihpb(i).eq.iss(j)) goto 10
 542           enddo
 543           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 544    10     continue
 545           do j=1,ns
 546             if (jhpb(i).eq.iss(j)) goto 20
 547           enddo
 548           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 549    20     continue
 550 c          dhpb(i)=dbr
 551 c          forcon(i)=fbr
 552         enddo
 553         do i=1,nss
 554           ihpb(i)=ihpb(i)+nres
 555           jhpb(i)=jhpb(i)+nres
 556         enddo
 557       endif
 558       endif
 559       if (ns.gt.0.and.dyn_ss) then
 560           do i=nss+1,nhpb
 561             ihpb(i-nss)=ihpb(i)
 562             jhpb(i-nss)=jhpb(i)
 563             forcon(i-nss)=forcon(i)
 564             dhpb(i-nss)=dhpb(i)
 565           enddo
 566           nhpb=nhpb-nss
 567           nss=0
 568           call hpb_partition
 569           do i=1,ns
 570             dyn_ss_mask(iss(i))=.true.
 571 c          write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
 572           enddo
 573       endif
 574       print *, "Leaving brigde read"
 575       return
 576       end
 577 c----------------------------------------------------------------------------
 578       subroutine read_angles(kanal,*)
 579       implicit none
 580       include 'DIMENSIONS'
 581       include 'COMMON.GEO'
 582       include 'COMMON.VAR'
 583       include 'COMMON.CHAIN'
 584       include 'COMMON.IOUNITS'
 585       integer i,kanal
 586       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 587       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 588       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 589       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 590       do i=1,nres
 591         theta(i)=deg2rad*theta(i)
 592         phi(i)=deg2rad*phi(i)
 593         alph(i)=deg2rad*alph(i)
 594         omeg(i)=deg2rad*omeg(i)
 595       enddo
 596       return
 597    10 return1
 598       end
 599 c----------------------------------------------------------------------------
 600       subroutine reada(rekord,lancuch,wartosc,default)
 601       implicit none
 602       character*(*) rekord,lancuch
 603       double precision wartosc,default
 604       integer ilen,iread
 605       external ilen
 606       iread=index(rekord,lancuch)
 607       if (iread.eq.0) then
 608         wartosc=default
 609         return
 610       endif
 611       iread=iread+ilen(lancuch)+1
 612       read (rekord(iread:),*) wartosc
 613       return
 614       end
 615 c----------------------------------------------------------------------------
 616       subroutine multreada(rekord,lancuch,tablica,dim,default)
 617       implicit none
 618       integer dim,i
 619       double precision tablica(dim),default
 620       character*(*) rekord,lancuch
 621       integer ilen,iread
 622       external ilen
 623       do i=1,dim
 624         tablica(i)=default
 625       enddo
 626       iread=index(rekord,lancuch)
 627       if (iread.eq.0) return
 628       iread=iread+ilen(lancuch)+1
 629       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 630    10 return
 631       end
 632 c----------------------------------------------------------------------------
 633       subroutine readi(rekord,lancuch,wartosc,default)
 634       implicit none
 635       character*(*) rekord,lancuch
 636       integer wartosc,default
 637       integer ilen,iread
 638       external ilen
 639       iread=index(rekord,lancuch)
 640       if (iread.eq.0) then
 641         wartosc=default
 642         return
 643       endif
 644       iread=iread+ilen(lancuch)+1
 645       read (rekord(iread:),*) wartosc
 646       return
 647       end
 648 c----------------------------------------------------------------------------
 649       subroutine multreadi(rekord,lancuch,tablica,dim,default)
 650       implicit none
 651       integer dim,i
 652       integer tablica(dim),default
 653       character*(*) rekord,lancuch
 654       character*80 aux
 655       integer ilen,iread
 656       external ilen
 657       do i=1,dim
 658         tablica(i)=default
 659       enddo
 660       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 661       if (iread.eq.0) return
 662       iread=iread+ilen(lancuch)+1
 663       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 664    10 return
 665       end
 666 c----------------------------------------------------------------------------
 667       subroutine card_concat(card)
 668       include 'DIMENSIONS'
 669       include 'COMMON.IOUNITS'
 670       character*(*) card
 671       character*80 karta,ucase
 672       external ilen
 673       read (inp,'(a)') karta
 674       karta=ucase(karta)
 675       card=' '
 676       do while (karta(80:80).eq.'&')
 677         card=card(:ilen(card)+1)//karta(:79)
 678         read (inp,'(a)') karta
 679         karta=ucase(karta)
 680       enddo
 681       card=card(:ilen(card)+1)//karta
 682       return
 683       end
 684 c----------------------------------------------------------------------------
 685       subroutine openunits
 686       implicit none
 687       include 'DIMENSIONS'
 688 #ifdef MPI
 689       include "mpif.h"
 690       character*3 liczba
 691       include "COMMON.MPI"
 692 #endif
 693       include 'COMMON.IOUNITS'
 694       include 'COMMON.CONTROL'
 695       integer lenpre,lenpot,ilen
 696       external ilen
 697       character*16 cformat,cprint
 698       character*16 ucase
 699       integer lenint,lenout
 700       call getenv('INPUT',prefix)
 701       call getenv('OUTPUT',prefout)
 702       call getenv('INTIN',prefintin)
 703       call getenv('COORD',cformat)
 704       call getenv('PRINTCOOR',cprint)
 705       call getenv('SCRATCHDIR',scratchdir)
 706       from_bx=.true.
 707       from_cx=.false.
 708       if (index(ucase(cformat),'CX').gt.0) then
 709         from_cx=.true.
 710         from_bx=.false.
 711       endif
 712       from_cart=.true.
 713       lenpre=ilen(prefix)
 714       lenout=ilen(prefout)
 715       lenint=ilen(prefintin)
 716 C Get the names and open the input files
 717       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
 718 #ifdef MPI
 719       write (liczba,'(bz,i3.3)') me
 720       outname=prefout(:lenout)//'_clust.out_'//liczba
 721 #else
 722       outname=prefout(:lenout)//'_clust.out'
 723 #endif
 724       if (from_bx) then
 725         intinname=prefintin(:lenint)//'.bx'
 726       else if (from_cx) then
 727         intinname=prefintin(:lenint)//'.cx'
 728       else
 729         intinname=prefintin(:lenint)//'.int'
 730       endif
 731       rmsname=prefintin(:lenint)//'.rms'
 732       open (jplot,file=prefout(:ilen(prefout))//'.tex',
 733      &       status='unknown')
 734       open (jrms,file=rmsname,status='unknown')
 735       open(iout,file=outname,status='unknown')
 736 C Get parameter filenames and open the parameter files.
 737       call getenv('BONDPAR',bondname)
 738       open (ibond,file=bondname,status='old')
 739       call getenv('THETPAR',thetname)
 740       open (ithep,file=thetname,status='old')
 741       call getenv('ROTPAR',rotname)
 742       open (irotam,file=rotname,status='old')
 743       call getenv('TORPAR',torname)
 744       open (itorp,file=torname,status='old')
 745       call getenv('TORDPAR',tordname)
 746       open (itordp,file=tordname,status='old')
 747       call getenv('FOURIER',fouriername)
 748       open (ifourier,file=fouriername,status='old')
 749       call getenv('ELEPAR',elename)
 750       open (ielep,file=elename,status='old')
 751       call getenv('SIDEPAR',sidename)
 752       open (isidep,file=sidename,status='old')
 753       call getenv('SIDEP',sidepname)
 754       open (isidep1,file=sidepname,status="old")
 755       call getenv('SCCORPAR',sccorname)
 756       open (isccor,file=sccorname,status="old")
 757 #ifndef OLDSCP
 758 C
 759 C 8/9/01 In the newest version SCp interaction constants are read from a file
 760 C Use -DOLDSCP to use hard-coded constants instead.
 761 C
 762       call getenv('SCPPAR',scpname)
 763       open (iscpp,file=scpname,status='old')
 764 #endif
 765       return
 766       end
 767 c-------------------------------------------------------------------------------
 768       subroutine read_dist_constr
 769       implicit real*8 (a-h,o-z)
 770       include 'DIMENSIONS'
 771       include 'COMMON.CONTROL'
 772       include 'COMMON.CHAIN'
 773       include 'COMMON.IOUNITS'
 774       include 'COMMON.SBRIDGE'
 775       integer ifrag_(2,100),ipair_(2,100)
 776       double precision wfrag_(100),wpair_(100)
 777       character*500 controlcard
 778 c      write (iout,*) "Calling read_dist_constr"
 779 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 780       call card_concat(controlcard)
 781 c      call flush(iout)
 782 c      write (iout,'(a)') controlcard
 783       call readi(controlcard,"NFRAG",nfrag_,0)
 784       call readi(controlcard,"NPAIR",npair_,0)
 785       call readi(controlcard,"NDIST",ndist_,0)
 786       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 787       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 788       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 789       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 790       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 791       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 792       write (iout,*) "IFRAG"
 793       do i=1,nfrag_
 794         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 795       enddo
 796       write (iout,*) "IPAIR"
 797       do i=1,npair_
 798         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 799       enddo
 800 #ifdef AIX
 801       call flush_(iout)
 802 #else
 803       call flush(iout)
 804 #endif
 805       if (.not.refstr .and. nfrag_.gt.0) then
 806         write (iout,*)
 807      &  "ERROR: no reference structure to compute distance restraints"
 808         write (iout,*)
 809      &  "Restraints must be specified explicitly (NDIST=number)"
 810         stop
 811       endif
 812       if (nfrag_.lt.2 .and. npair_.gt.0) then
 813         write (iout,*) "ERROR: Less than 2 fragments specified",
 814      &   " but distance restraints between pairs requested"
 815         stop
 816       endif
 817 #ifdef AIX
 818       call flush_(iout)
 819 #else
 820       call flush(iout)
 821 #endif
 822       do i=1,nfrag_
 823         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 824         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 825      &    ifrag_(2,i)=nstart_sup+nsup-1
 826 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 827 c        call flush(iout)
 828         if (wfrag_(i).gt.0.0d0) then
 829         do j=ifrag_(1,i),ifrag_(2,i)-1
 830           do k=j+1,ifrag_(2,i)
 831             write (iout,*) "j",j," k",k
 832             ddjk=dist(j,k)
 833             if (constr_dist.eq.1) then
 834             nhpb=nhpb+1
 835             ihpb(nhpb)=j
 836             jhpb(nhpb)=k
 837               dhpb(nhpb)=ddjk
 838             forcon(nhpb)=wfrag_(i)
 839             else if (constr_dist.eq.2) then
 840               if (ddjk.le.dist_cut) then
 841                 nhpb=nhpb+1
 842                 ihpb(nhpb)=j
 843                 jhpb(nhpb)=k
 844                 dhpb(nhpb)=ddjk
 845                 forcon(nhpb)=wfrag_(i)
 846               endif
 847             else
 848               nhpb=nhpb+1
 849               ihpb(nhpb)=j
 850               jhpb(nhpb)=k
 851               dhpb(nhpb)=ddjk
 852               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 853             endif
 854             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
 855      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 856           enddo
 857         enddo
 858         endif
 859       enddo
 860       do i=1,npair_
 861         if (wpair_(i).gt.0.0d0) then
 862         ii = ipair_(1,i)
 863         jj = ipair_(2,i)
 864         if (ii.gt.jj) then
 865           itemp=ii
 866           ii=jj
 867           jj=itemp
 868         endif
 869         do j=ifrag_(1,ii),ifrag_(2,ii)
 870           do k=ifrag_(1,jj),ifrag_(2,jj)
 871             nhpb=nhpb+1
 872             ihpb(nhpb)=j
 873             jhpb(nhpb)=k
 874             forcon(nhpb)=wpair_(i)
 875             dhpb(nhpb)=dist(j,k)
 876             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 877      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 878           enddo
 879         enddo
 880         endif
 881       enddo
 882       do i=1,ndist_
 883        if (constr_dist.eq.11) then
 884          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
 885      &     ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
 886          fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
 887        else
 888         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
 889      &     ibecarb(i),forcon(nhpb+1)
 890        endif
 891         if (forcon(nhpb+1).gt.0.0d0) then
 892           nhpb=nhpb+1
 893           if (ibecarb(i).gt.0) then
 894             ihpb(i)=ihpb(i)+nres
 895             jhpb(i)=jhpb(i)+nres
 896           endif
 897           if (dhpb(nhpb).eq.0.0d0)
 898      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 899         endif
 900       enddo
 901       do i=1,nhpb
 902        if (constr_dist.eq.11) then
 903           write (iout,'(a,3i5,2f8.2,i2,2f10.1)') "+dist.constr11 ",
 904      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i),
 905      &     fordepth(i)
 906        else
 907           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
 908      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
 909        endif
 910       enddo
 911 #ifdef AIX
 912       call flush_(iout)
 913 #else
 914       call flush(iout)
 915 #endif
 916       return
 917       end
 918
 919 c====-------------------------------------------------------------------
 920       subroutine read_constr_homology(lprn)
 921
 922       include 'DIMENSIONS'
 923 #ifdef MPI
 924       include 'mpif.h'
 925 #endif
 926       include 'COMMON.SETUP'
 927       include 'COMMON.CONTROL'
 928       include 'COMMON.CHAIN'
 929       include 'COMMON.IOUNITS'
 930       include 'COMMON.GEO'
 931       include 'COMMON.INTERACT'
 932       include 'COMMON.NAMES'
 933       include 'COMMON.HOMRESTR'
 934 c
 935 c For new homol impl
 936 c
 937       include 'COMMON.VAR'
 938 c     include 'include_unres/COMMON.VAR'
 939 c
 940
 941 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
 942 c    &                 dist_cut
 943 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
 944 c    &    sigma_odl_temp(maxres,maxres,max_template)
 945       character*2 kic2
 946       character*24 model_ki_dist, model_ki_angle
 947       character*500 controlcard
 948       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
 949       integer idomain(max_template,maxres)
 950       integer ilen
 951       external ilen
 952       logical lprn
 953       logical unres_pdb,liiflag
 954 c
 955 c     FP - Nov. 2014 Temporary specifications for new vars
 956 c
 957       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
 958       double precision, dimension (max_template,maxres) :: rescore
 959       double precision, dimension (max_template,maxres) :: rescore2
 960       character*24 tpl_k_rescore
 961 c -----------------------------------------------------------------
 962 c Reading multiple PDB ref structures and calculation of retraints
 963 c not using pre-computed ones stored in files model_ki_{dist,angle}
 964 c FP (Nov., 2014)
 965 c -----------------------------------------------------------------
 966 c
 967 c
 968 c Alternative: reading from input
 969 #ifdef DEBUG
 970       write (iout,*) "BEGIN READ HOMOLOGY INFO"
 971 #ifdef AIX
 972       call flush_(iout)
 973 #else
 974       call flush(iout)
 975 #endif
 976 #endif
 977       call card_concat(controlcard)
 978       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
 979       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
 980       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
 981       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
 982       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
 983       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
 984       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
 985       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
 986       if (homol_nset.gt.1)then
 987          call readi(controlcard,"ISET",iset,1)
 988          call card_concat(controlcard)
 989          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
 990       else
 991         iset=1
 992         waga_homology(1)=1.0
 993       endif
 994 c
 995 #ifdef DEBUG
 996       write(iout,*) "read_constr_homology iset",iset
 997       write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
 998 #ifdef AIX
 999       call flush_(iout)
1000 #else
1001       call flush(iout)
1002 #endif
1003 #endif
1004 cd      write (iout,*) "nnt",nnt," nct",nct
1005 cd      call flush(iout)
1006
1007
1008       lim_odl=0
1009       lim_dih=0
1010 c
1011 c  Reading HM global scores (prob not required)
1012 c
1013       do i = nnt,nct
1014         do k=1,constr_homology
1015          idomain(k,i)=0
1016         enddo
1017       enddo
1018 c     open (4,file="HMscore")
1019 c     do k=1,constr_homology
1020 c       read (4,*,end=521) hmscore_tmp
1021 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
1022 c       write(*,*) "Model", k, ":", hmscore(k)
1023 c     enddo
1024 c521  continue
1025
1026       ii=0
1027       do i = nnt,nct-2
1028         do j=i+2,nct
1029         ii=ii+1
1030         ii_in_use(ii)=0
1031         enddo
1032       enddo
1033 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1034
1035       write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1036       do k=1,constr_homology
1037
1038         read(inp,'(a)') pdbfile
1039 c        write (iout,*) "k ",k," pdbfile ",pdbfile
1040 c  Next stament causes error upon compilation (?)
1041 c       if(me.eq.king.or. .not. out1file)
1042 c         write (iout,'(2a)') 'PDB data will be read from file ',
1043 c    &   pdbfile(:ilen(pdbfile))
1044          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1045      &  pdbfile(:ilen(pdbfile))
1046         open(ipdbin,file=pdbfile,status='old',err=33)
1047         goto 34
1048   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
1049      &  pdbfile(:ilen(pdbfile))
1050         stop
1051   34    continue
1052 c        print *,'Begin reading pdb data'
1053 c
1054 c Files containing res sim or local scores (former containing sigmas)
1055 c
1056
1057         write(kic2,'(bz,i2.2)') k
1058
1059         tpl_k_rescore="template"//kic2//".sco"
1060
1061         unres_pdb=.false.
1062         call readpdb
1063         do i=1,2*nres
1064           do j=1,3
1065             crefjlee(j,i)=c(j,i)
1066           enddo
1067         enddo
1068 #ifdef DEBUG
1069         do i=1,nres
1070           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1071      &      (crefjlee(j,i+nres),j=1,3)
1072         enddo
1073         write (iout,*) "READ HOMOLOGY INFO"
1074         write (iout,*) "read_constr_homology x: after reading pdb file"
1075         write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1076         write (iout,*) "waga_dist",waga_dist
1077         write (iout,*) "waga_angle",waga_angle
1078         write (iout,*) "waga_theta",waga_theta
1079         write (iout,*) "waga_d",waga_d
1080         write (iout,*) "dist_cut",dist_cut
1081 #endif
1082 #ifdef AIX
1083       call flush_(iout)
1084 #else
1085       call flush(iout)
1086 #endif
1087
1088 c
1089 c     Distance restraints
1090 c
1091 c          ... --> odl(k,ii)
1092 C Copy the coordinates from reference coordinates (?)
1093         do i=1,2*nres
1094           do j=1,3
1095             c(j,i)=cref(j,i)
1096 c           write (iout,*) "c(",j,i,") =",c(j,i)
1097           enddo
1098         enddo
1099 c
1100 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1101 c
1102 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1103           open (ientin,file=tpl_k_rescore,status='old')
1104           if (nnt.gt.1) rescore(k,1)=0.0d0
1105           do irec=nnt,nct ! loop for reading res sim
1106             if (read2sigma) then
1107              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1108      &                                idomain_tmp
1109              i_tmp=i_tmp+nnt-1
1110              idomain(k,i_tmp)=idomain_tmp
1111              rescore(k,i_tmp)=rescore_tmp
1112              rescore2(k,i_tmp)=rescore2_tmp
1113              write(iout,'(a7,i5,2f10.5,i5)') "rescore",
1114      &                      i_tmp,rescore2_tmp,rescore_tmp,
1115      &                                idomain_tmp
1116             else
1117              idomain(k,irec)=1
1118              read (ientin,*,end=1401) rescore_tmp
1119
1120 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1121              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1122 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1123             endif
1124           enddo
1125  1401   continue
1126         close (ientin)
1127         if (waga_dist.ne.0.0d0) then
1128           ii=0
1129           do i = nnt,nct-2
1130             do j=i+2,nct
1131
1132               x12=c(1,i)-c(1,j)
1133               y12=c(2,i)-c(2,j)
1134               z12=c(3,i)-c(3,j)
1135               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1136 c              write (iout,*) k,i,j,distal,dist2_cut
1137
1138             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1139      &            .and. distal.le.dist2_cut ) then
1140
1141               ii=ii+1
1142               ii_in_use(ii)=1
1143               l_homo(k,ii)=.true.
1144
1145 c             write (iout,*) "k",k
1146 c             write (iout,*) "i",i," j",j," constr_homology",
1147 c    &                       constr_homology
1148               ires_homo(ii)=i
1149               jres_homo(ii)=j
1150               odl(k,ii)=distal
1151               if (read2sigma) then
1152                 sigma_odl(k,ii)=0
1153                 do ik=i,j
1154                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1155                 enddo
1156                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1157                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1158      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1159               else
1160                 if (odl(k,ii).le.dist_cut) then
1161                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1162                 else
1163 #ifdef OLDSIGMA
1164                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1165      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1166 #else
1167                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1168      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1169 #endif
1170                 endif
1171               endif
1172               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1173             else
1174               ii=ii+1
1175               l_homo(k,ii)=.false.
1176             endif
1177             enddo
1178           enddo
1179         lim_odl=ii
1180         endif
1181 c
1182 c     Theta, dihedral and SC retraints
1183 c
1184         if (waga_angle.gt.0.0d0) then
1185 c         open (ientin,file=tpl_k_sigma_dih,status='old')
1186 c         do irec=1,maxres-3 ! loop for reading sigma_dih
1187 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1188 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1189 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1190 c    &                            sigma_dih(k,i+nnt-1)
1191 c         enddo
1192 c1402   continue
1193 c         close (ientin)
1194           do i = nnt+3,nct
1195             if (idomain(k,i).eq.0) then
1196                sigma_dih(k,i)=0.0
1197                cycle
1198             endif
1199             dih(k,i)=phiref(i) ! right?
1200 c           read (ientin,*) sigma_dih(k,i) ! original variant
1201 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
1202 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1203 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1204 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
1205 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
1206
1207             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1208      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1209 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1210 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
1211 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1212 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1213 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
1214 c   Instead of res sim other local measure of b/b str reliability possible
1215             if (sigma_dih(k,i).ne.0)
1216      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1217 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1218           enddo
1219           lim_dih=nct-nnt-2
1220         endif
1221
1222         if (waga_theta.gt.0.0d0) then
1223 c         open (ientin,file=tpl_k_sigma_theta,status='old')
1224 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1225 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1226 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1227 c    &                              sigma_theta(k,i+nnt-1)
1228 c         enddo
1229 c1403   continue
1230 c         close (ientin)
1231
1232           do i = nnt+2,nct ! right? without parallel.
1233 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
1234 c         do i=ithet_start,ithet_end ! with FG parallel.
1235              if (idomain(k,i).eq.0) then
1236               sigma_theta(k,i)=0.0
1237               cycle
1238              endif
1239              thetatpl(k,i)=thetaref(i)
1240 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1241 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
1242 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1243 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
1244 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
1245              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1246      &                        rescore(k,i-2))/3.0
1247 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1248              if (sigma_theta(k,i).ne.0)
1249      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1250
1251 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1252 c                             rescore(k,i-2) !  right expression ?
1253 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1254           enddo
1255         endif
1256
1257         if (waga_d.gt.0.0d0) then
1258 c       open (ientin,file=tpl_k_sigma_d,status='old')
1259 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1260 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1261 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1262 c    &                          sigma_d(k,i+nnt-1)
1263 c         enddo
1264 c1404   continue
1265
1266           do i = nnt,nct ! right? without parallel.
1267 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
1268 c         do i=loc_start,loc_end ! with FG parallel.
1269                if (itype(i).eq.10) cycle
1270                if (idomain(k,i).eq.0 ) then
1271                   sigma_d(k,i)=0.0
1272                   cycle
1273                endif
1274                xxtpl(k,i)=xxref(i)
1275                yytpl(k,i)=yyref(i)
1276                zztpl(k,i)=zzref(i)
1277 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1278 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1279 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1280 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
1281                sigma_d(k,i)=rescore(k,i) !  right expression ?
1282                if (sigma_d(k,i).ne.0)
1283      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1284
1285 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
1286 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1287 c              read (ientin,*) sigma_d(k,i) ! 1st variant
1288                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
1289           enddo
1290         endif
1291       enddo
1292 c
1293 c remove distance restraints not used in any model from the list
1294 c shift data in all arrays
1295 c
1296       if (waga_dist.ne.0.0d0) then
1297         ii=0
1298         liiflag=.true.
1299         do i=nnt,nct-2
1300          do j=i+2,nct
1301           ii=ii+1
1302           if (ii_in_use(ii).eq.0.and.liiflag) then
1303             liiflag=.false.
1304             iistart=ii
1305           endif
1306           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1307      &                   .not.liiflag.and.ii.eq.lim_odl) then
1308              if (ii.eq.lim_odl) then
1309               iishift=ii-iistart+1
1310              else
1311               iishift=ii-iistart
1312              endif
1313              liiflag=.true.
1314              do ki=iistart,lim_odl-iishift
1315               ires_homo(ki)=ires_homo(ki+iishift)
1316               jres_homo(ki)=jres_homo(ki+iishift)
1317               ii_in_use(ki)=ii_in_use(ki+iishift)
1318               do k=1,constr_homology
1319                odl(k,ki)=odl(k,ki+iishift)
1320                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1321                l_homo(k,ki)=l_homo(k,ki+iishift)
1322               enddo
1323              enddo
1324              ii=ii-iishift
1325              lim_odl=lim_odl-iishift
1326           endif
1327          enddo
1328         enddo
1329       endif
1330       if (constr_homology.gt.0) call homology_partition
1331       if (constr_homology.gt.0) call init_int_table
1332 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1333 cd     & "lim_xx=",lim_xx
1334 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1335 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1336 c
1337 c Print restraints
1338 c
1339       if (.not.lprn) return
1340 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1341       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1342        write (iout,*) "Distance restraints from templates"
1343        do ii=1,lim_odl
1344        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1345      &  ii,ires_homo(ii),jres_homo(ii),
1346      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1347      &  ki=1,constr_homology)
1348        enddo
1349        write (iout,*) "Dihedral angle restraints from templates"
1350        do i=nnt+3,nct
1351         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1352      &      (rad2deg*dih(ki,i),
1353      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1354        enddo
1355        write (iout,*) "Virtual-bond angle restraints from templates"
1356        do i=nnt+2,nct
1357         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1358      &      (rad2deg*thetatpl(ki,i),
1359      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1360        enddo
1361        write (iout,*) "SC restraints from templates"
1362        do i=nnt,nct
1363         write(iout,'(i5,100(4f8.2,4x))') i,
1364      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1365      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1366        enddo
1367       endif
1368 c -----------------------------------------------------------------
1369       return
1370       end
1371 c----------------------------------------------------------------------