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         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
 884      &     ibecarb(i),forcon(nhpb+1)
 885         if (forcon(nhpb+1).gt.0.0d0) then
 886           nhpb=nhpb+1
 887           if (ibecarb(i).gt.0) then
 888             ihpb(i)=ihpb(i)+nres
 889             jhpb(i)=jhpb(i)+nres
 890           endif
 891           if (dhpb(nhpb).eq.0.0d0)
 892      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 893         endif
 894       enddo
 895       do i=1,nhpb
 896           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
 897      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
 898       enddo
 899 #ifdef AIX
 900       call flush_(iout)
 901 #else
 902       call flush(iout)
 903 #endif
 904       return
 905       end
 906
 907 c====-------------------------------------------------------------------
 908       subroutine read_constr_homology(lprn)
 909
 910       include 'DIMENSIONS'
 911 #ifdef MPI
 912       include 'mpif.h'
 913 #endif
 914       include 'COMMON.SETUP'
 915       include 'COMMON.CONTROL'
 916       include 'COMMON.CHAIN'
 917       include 'COMMON.IOUNITS'
 918       include 'COMMON.GEO'
 919       include 'COMMON.INTERACT'
 920       include 'COMMON.NAMES'
 921       include 'COMMON.HOMRESTR'
 922 c
 923 c For new homol impl
 924 c
 925       include 'COMMON.VAR'
 926 c     include 'include_unres/COMMON.VAR'
 927 c
 928
 929 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
 930 c    &                 dist_cut
 931 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
 932 c    &    sigma_odl_temp(maxres,maxres,max_template)
 933       character*2 kic2
 934       character*24 model_ki_dist, model_ki_angle
 935       character*500 controlcard
 936       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
 937       integer idomain(max_template,maxres)
 938       integer ilen
 939       external ilen
 940       logical lprn
 941       logical unres_pdb
 942 c
 943 c     FP - Nov. 2014 Temporary specifications for new vars
 944 c
 945       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
 946       double precision, dimension (max_template,maxres) :: rescore
 947       double precision, dimension (max_template,maxres) :: rescore2
 948       character*24 tpl_k_rescore
 949 c -----------------------------------------------------------------
 950 c Reading multiple PDB ref structures and calculation of retraints
 951 c not using pre-computed ones stored in files model_ki_{dist,angle}
 952 c FP (Nov., 2014)
 953 c -----------------------------------------------------------------
 954 c
 955 c
 956 c Alternative: reading from input
 957 #ifdef DEBUG
 958       write (iout,*) "BEGIN READ HOMOLOGY INFO"
 959 #ifdef AIX
 960       call flush_(iout)
 961 #else
 962       call flush(iout)
 963 #endif
 964 #endif
 965       call card_concat(controlcard)
 966       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
 967       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
 968       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
 969       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
 970       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
 971       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
 972       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
 973       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
 974       if (homol_nset.gt.1)then
 975          call readi(controlcard,"ISET",iset,1)
 976          call card_concat(controlcard)
 977          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
 978       else
 979         iset=1
 980         waga_homology(1)=1.0
 981       endif
 982 c
 983 #ifdef DEBUG
 984       write(iout,*) "read_constr_homology iset",iset
 985       write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
 986 #ifdef AIX
 987       call flush_(iout)
 988 #else
 989       call flush(iout)
 990 #endif
 991 #endif
 992 cd      write (iout,*) "nnt",nnt," nct",nct
 993 cd      call flush(iout)
 994
 995
 996       lim_odl=0
 997       lim_dih=0
 998 c
 999 c  Reading HM global scores (prob not required)
1000 c
1001       do i = nnt,nct
1002         do k=1,constr_homology
1003          idomain(k,i)=0
1004         enddo
1005       enddo
1006 c     open (4,file="HMscore")
1007 c     do k=1,constr_homology
1008 c       read (4,*,end=521) hmscore_tmp
1009 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
1010 c       write(*,*) "Model", k, ":", hmscore(k)
1011 c     enddo
1012 c521  continue
1013
1014       ii=0
1015       do i = nnt,nct-2
1016         do j=i+2,nct
1017         ii=ii+1
1018         ii_in_use(ii)=0
1019         enddo
1020       enddo
1021 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1022
1023       write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1024       do k=1,constr_homology
1025
1026         read(inp,'(a)') pdbfile
1027 c        write (iout,*) "k ",k," pdbfile ",pdbfile
1028 c  Next stament causes error upon compilation (?)
1029 c       if(me.eq.king.or. .not. out1file)
1030 c         write (iout,'(2a)') 'PDB data will be read from file ',
1031 c    &   pdbfile(:ilen(pdbfile))
1032          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1033      &  pdbfile(:ilen(pdbfile))
1034         open(ipdbin,file=pdbfile,status='old',err=33)
1035         goto 34
1036   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
1037      &  pdbfile(:ilen(pdbfile))
1038         stop
1039   34    continue
1040 c        print *,'Begin reading pdb data'
1041 c
1042 c Files containing res sim or local scores (former containing sigmas)
1043 c
1044
1045         write(kic2,'(bz,i2.2)') k
1046
1047         tpl_k_rescore="template"//kic2//".sco"
1048
1049         unres_pdb=.false.
1050         call readpdb
1051         do i=1,2*nres
1052           do j=1,3
1053             crefjlee(j,i)=c(j,i)
1054           enddo
1055         enddo
1056 #ifdef DEBUG
1057         do i=1,nres
1058           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1059      &      (crefjlee(j,i+nres),j=1,3)
1060         enddo
1061         write (iout,*) "READ HOMOLOGY INFO"
1062         write (iout,*) "read_constr_homology x: after reading pdb file"
1063         write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1064         write (iout,*) "waga_dist",waga_dist
1065         write (iout,*) "waga_angle",waga_angle
1066         write (iout,*) "waga_theta",waga_theta
1067         write (iout,*) "waga_d",waga_d
1068         write (iout,*) "dist_cut",dist_cut
1069 #endif
1070 #ifdef AIX
1071       call flush_(iout)
1072 #else
1073       call flush(iout)
1074 #endif
1075
1076 c
1077 c     Distance restraints
1078 c
1079 c          ... --> odl(k,ii)
1080 C Copy the coordinates from reference coordinates (?)
1081         do i=1,2*nres
1082           do j=1,3
1083             c(j,i)=cref(j,i)
1084 c           write (iout,*) "c(",j,i,") =",c(j,i)
1085           enddo
1086         enddo
1087 c
1088 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1089 c
1090 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1091           open (ientin,file=tpl_k_rescore,status='old')
1092           if (nnt.gt.1) rescore(k,1)=0.0d0
1093           do irec=nnt,nct ! loop for reading res sim
1094             if (read2sigma) then
1095              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1096      &                                idomain_tmp
1097              i_tmp=i_tmp+nnt-1
1098              idomain(k,i_tmp)=idomain_tmp
1099              rescore(k,i_tmp)=rescore_tmp
1100              rescore2(k,i_tmp)=rescore2_tmp
1101              write(iout,'(a7,i5,2f10.5,i5)') "rescore",
1102      &                      i_tmp,rescore2_tmp,rescore_tmp,
1103      &                                idomain_tmp
1104             else
1105              idomain(k,irec)=1
1106              read (ientin,*,end=1401) rescore_tmp
1107
1108 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1109              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1110 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1111             endif
1112           enddo
1113  1401   continue
1114         close (ientin)
1115         if (waga_dist.ne.0.0d0) then
1116           ii=0
1117           do i = nnt,nct-2
1118             do j=i+2,nct
1119
1120               x12=c(1,i)-c(1,j)
1121               y12=c(2,i)-c(2,j)
1122               z12=c(3,i)-c(3,j)
1123               distal=dsqrt(x12*x12+y12*y12+z12*z12)
1124 c              write (iout,*) k,i,j,distal,dist2_cut
1125
1126             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1127      &            .and. distal.le.dist2_cut ) then
1128
1129               ii=ii+1
1130               ii_in_use(ii)=1
1131               l_homo(k,ii)=.true.
1132
1133 c             write (iout,*) "k",k
1134 c             write (iout,*) "i",i," j",j," constr_homology",
1135 c    &                       constr_homology
1136               ires_homo(ii)=i
1137               jres_homo(ii)=j
1138               odl(k,ii)=distal
1139               if (read2sigma) then
1140                 sigma_odl(k,ii)=0
1141                 do ik=i,j
1142                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1143                 enddo
1144                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1145                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1146      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1147               else
1148                 if (odl(k,ii).le.dist_cut) then
1149                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1150                 else
1151 #ifdef OLDSIGMA
1152                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1153      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1154 #else
1155                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1156      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1157 #endif
1158                 endif
1159               endif
1160               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1161             else
1162               ii=ii+1
1163               l_homo(k,ii)=.false.
1164             endif
1165             enddo
1166           enddo
1167         lim_odl=ii
1168         endif
1169 c
1170 c     Theta, dihedral and SC retraints
1171 c
1172         if (waga_angle.gt.0.0d0) then
1173 c         open (ientin,file=tpl_k_sigma_dih,status='old')
1174 c         do irec=1,maxres-3 ! loop for reading sigma_dih
1175 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1176 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1177 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1178 c    &                            sigma_dih(k,i+nnt-1)
1179 c         enddo
1180 c1402   continue
1181 c         close (ientin)
1182           do i = nnt+3,nct
1183             if (idomain(k,i).eq.0) then
1184                sigma_dih(k,i)=0.0
1185                cycle
1186             endif
1187             dih(k,i)=phiref(i) ! right?
1188 c           read (ientin,*) sigma_dih(k,i) ! original variant
1189 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
1190 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1191 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1192 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
1193 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
1194
1195             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1196      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1197 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1198 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
1199 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1200 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1201 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
1202 c   Instead of res sim other local measure of b/b str reliability possible
1203             if (sigma_dih(k,i).ne.0)
1204      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1205 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1206           enddo
1207           lim_dih=nct-nnt-2
1208         endif
1209
1210         if (waga_theta.gt.0.0d0) then
1211 c         open (ientin,file=tpl_k_sigma_theta,status='old')
1212 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1213 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1214 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1215 c    &                              sigma_theta(k,i+nnt-1)
1216 c         enddo
1217 c1403   continue
1218 c         close (ientin)
1219
1220           do i = nnt+2,nct ! right? without parallel.
1221 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
1222 c         do i=ithet_start,ithet_end ! with FG parallel.
1223              if (idomain(k,i).eq.0) then
1224               sigma_theta(k,i)=0.0
1225               cycle
1226              endif
1227              thetatpl(k,i)=thetaref(i)
1228 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1229 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
1230 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1231 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
1232 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
1233              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1234      &                        rescore(k,i-2))/3.0
1235 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1236              if (sigma_theta(k,i).ne.0)
1237      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1238
1239 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1240 c                             rescore(k,i-2) !  right expression ?
1241 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1242           enddo
1243         endif
1244
1245         if (waga_d.gt.0.0d0) then
1246 c       open (ientin,file=tpl_k_sigma_d,status='old')
1247 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1248 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1249 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1250 c    &                          sigma_d(k,i+nnt-1)
1251 c         enddo
1252 c1404   continue
1253
1254           do i = nnt,nct ! right? without parallel.
1255 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
1256 c         do i=loc_start,loc_end ! with FG parallel.
1257                if (itype(i).eq.10) cycle
1258                if (idomain(k,i).eq.0 ) then
1259                   sigma_d(k,i)=0.0
1260                   cycle
1261                endif
1262                xxtpl(k,i)=xxref(i)
1263                yytpl(k,i)=yyref(i)
1264                zztpl(k,i)=zzref(i)
1265 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1266 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1267 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1268 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
1269                sigma_d(k,i)=rescore(k,i) !  right expression ?
1270                if (sigma_d(k,i).ne.0)
1271      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1272
1273 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
1274 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1275 c              read (ientin,*) sigma_d(k,i) ! 1st variant
1276                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
1277           enddo
1278         endif
1279       enddo
1280 c
1281 c remove distance restraints not used in any model from the list
1282 c shift data in all arrays
1283 c
1284       if (waga_dist.ne.0.0d0) then
1285         ii=0
1286         do i=nnt,nct-2
1287          do j=i+2,nct
1288           ii=ii+1
1289           if (ii_in_use(ii).eq.0) then
1290              do ki=ii,lim_odl-1
1291               ires_homo(ki)=ires_homo(ki+1)
1292               jres_homo(ki)=jres_homo(ki+1)
1293               ii_in_use(ki)=ii_in_use(ki+1)
1294               do k=1,constr_homology
1295                odl(k,ki)=odl(k,ki+1)
1296                sigma_odl(k,ki)=sigma_odl(k,ki+1)
1297                l_homo(k,ki)=l_homo(k,ki+1)
1298               enddo
1299              enddo
1300              ii=ii-1
1301              lim_odl=lim_odl-1
1302           endif
1303          enddo
1304         enddo
1305       endif
1306       if (constr_homology.gt.0) call homology_partition
1307       if (constr_homology.gt.0) call init_int_table
1308 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1309 cd     & "lim_xx=",lim_xx
1310 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1311 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1312 c
1313 c Print restraints
1314 c
1315       if (.not.lprn) return
1316 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1317       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1318        write (iout,*) "Distance restraints from templates"
1319        do ii=1,lim_odl
1320        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
1321      &  ii,ires_homo(ii),jres_homo(ii),
1322      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1323      &  ki=1,constr_homology)
1324        enddo
1325        write (iout,*) "Dihedral angle restraints from templates"
1326        do i=nnt+3,nct
1327         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1328      &      (rad2deg*dih(ki,i),
1329      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1330        enddo
1331        write (iout,*) "Virtual-bond angle restraints from templates"
1332        do i=nnt+2,nct
1333         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1334      &      (rad2deg*thetatpl(ki,i),
1335      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1336        enddo
1337        write (iout,*) "SC restraints from templates"
1338        do i=nnt,nct
1339         write(iout,'(i5,100(4f8.2,4x))') i,
1340      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1341      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1342        enddo
1343       endif
1344 c -----------------------------------------------------------------
1345       return
1346       end
1347 c----------------------------------------------------------------------