added source code

[unres.git] / doc / WHAM.TXT

          WHAM (Weighted Histogram Analysis Method)
          Processing results of UNRES/MREMD simulations
          ---------------------------------------------

TABLE OF CONTENTS
-----------------

1. License terms

2. References

3. Functions of the program

4. Installation

5. Running the program

6. Input and output files 
   6.1. Summary of files
   6.2. The main input file
        6.2.1. General data
        6.2.2 Molecule and energy parameter data
              6.2.2.1. General information
              6.2.2.2. Sequence information
              6.2.2.3. Dihedral angle restraint information
              6.2.2.4. Disulfide-bridge data
        6.2.3. Energy-term weights and parameter files
        6.2.4. (M)REMD/Hamiltonian (M)REMD setting specification
        6.2.5. Information of files from which to read conformations
        6.2.6. Information of reference structure and comparing scheme
   6.3. The structure of the main output file (out)
   6.4. The thermodynamic quantity and ensemble average (stat) files
   6.5. The conformation summary with classification (stat) files
   6.6. The histogram files
   6.7. The rmsd-radius of gyration potential of mean force files
   6.8. The PDB files
   6.8. The compresses Cartesian coordinates (cx) file.

7. Support

1. LICENSE TERMS
----------------

* This software is provided free of charge to academic users, subject to the
  condition that no part of it be sold or used otherwise for commercial
  purposes, including, but not limited to its incorporation into commercial
  software packages, without written consent from the authors. For permission
  contact Prof. H. A. Scheraga, Cornell University.

* This software package is provided on an "as is" basis. We in no way warrant
  either this software or results it may produce.

* Reports or publications using this software package must contain an
  acknowledgment to the authors and the NIH Resource in the form commonly
used
  in academic research.

2. REFERENCES
-------------

[1]  S. Kumar, D. Bouzida, R.H. Swendsen, P.A. Kollman, J.M. Rosenberg.
     The weighted histogram analysis method for free-energy calculations
     on biomolecules. I. The method.
     J. Comput. Chem., 1992, 13, 1011-1021.

[2]  A. Liwo, M. Khalili, C. Czaplewski, S. Kalinowski, S. Oldziej, K. Wachucik,
     H.A. Scheraga.
     Modification and optimization of the united-residue (UNRES) potential
     energy function for canonical simulations. I. Temperature dependence of the
     effective energy function and tests of the optimization method with single
     training proteins.
     J. Phys. Chem. B, 2007, 111, 260-285.

[3]  S. Oldziej, A. Liwo, C. Czaplewski, J. Pillardy, H.A. Scheraga.
     Optimization of the UNRES force field by hierarchical design of the
     potential-energy landscape. 2. Off-lattice tests of the method with single
     proteins.  J. Phys. Chem. B., 2004, 108, 16934-16949.

[4]  S. Oldziej, A. Liwo, C. Czaplewski, J. Pillardy, H.A. Scheraga.
     Optimization of the UNRES force field by hierarchical design of the
     potential-energy landscape. 2. Off-lattice tests of the method with single
     proteins.  J. Phys. Chem. B., 2004, 108, 16934-16949.

3. FUNCTIONS OF THE PROGRAM
---------------------------

The program processes the results of replica exchange (REMD) or multiplexed
replica exchange molecular dynamics (MREMD) simulations with UNRES to compute 
the probabilities of the obtained conformations to occur at particular
temperatures. The program is based on the variant of the weighted histogram
analysis (WHAM) method [1] described in ref [2].

The program outputs the following information:

a) Temperature profiles of thermodynamic and structural ensemble-averaged
   quantities.

b) Histograms of native-likeness measure q (defined by eqs 8-11 of ref [2]).

c) Optionally the most probable conformations at REMD temperatures.

d) Optionally the coordinates with information to compute probabilities
   for the conformations to occur at any temperature.

The program takes usually UNRES compressed coordinate files (cx files) from
MREMD obtained by using the TRAJ1FILE option. The user can request to 
partition the whole run into equal slices (or windows), each starting from,
say, snapshot n (for each trajectory) and ending at snapshot n+1.
Alternatively, the UNRES Cartesian coordinate (x files) can be input; however,
they must contain only the analyzed portion of the trajectories; they
are usually prepared from single trajectories by using xdrf2x.

Two versions of the program are provided:

a) Canonical version which treats single polypeptide chains; the source code
is in WHAM/src directory.

b) Version for oligomeric proteins; multiple chains are handled by inserting
dummy residues in the sequence; the source code is in WHAM/src-M directory.

4. INSTALLATION
---------------

Customize Makefile to your system. See section 7 of the description of UNRES
for compiler flags that are used to created executables for a particular
force field. There are already several Makefiles prepared for various systems
and force fields.

Run make in the WHAM/src directory WHAM/src-M directory for multichain
version. Make sure that MPI is installed on your system; the present program 
runs only in parallel mode.

5. RUNNING THE PROGRAM
----------------------

The program requires a parallel system to run. Depending on system,
either the wham.csh C-shell script (in WHAM/bin directory) can be started
using mpirun or the binary in the C-shell script must be executed through
mpirun. See the wham.csh C-shell script and section 6 for the files 
processed by the program.

6. INPUT AND OUTPUT FILES
-------------------------

6.1. SUMMARY OF THE FILES
-------------------------

The C-shell script wham.csh is used to run the program (see the WHAM/bin
directory). The data files that the script needs are mostly the same as 
for UNRES (see section 6 of UNRES description). In addition, the environmental
variable CONTFUN specifies the method to assess whether two side chains
are at contact; if CONTFUN=GB, the criterion defined by eq 8 of ref 4 is
used to assess whether two side chains are at contact. Also, the parameter
files from the C-shell scripts are overridden if the data from Hamiltonian
MREMD are processed; if so, the parameter files are defined in the main
input file.

The main input file must have inp extension. If it is INPUT.inp, the output
files are as follows:

INPUT.out_POTxxx - output files from different processors (INPUT.out_000 is the
     main output file). POT is the identifier of the sidechain-sidechain
     potential.

INPUT_POT_GB_xxx.stat or INPUT_POT_slice_YYXXX.stat- the summary conformation-
    classification file from processor xxx (each processor handles part of 
    conformations); the second occurs if the run is partitioned into slices.

INPUT.thermal or INPUT_slice_yy.thermal - thermodynamic functions and
    temperature profiles of the ensemble averages (the second form if the
    run is partitioned into slices).

INPUT_T_xxx.pdb or INPUT_slice_yy_T_xxx.pdb - top conformations the number
    of these conformations is selected by the user) in PDB format.

INPUT.cx - the compressed UNRES coordinate file with information to compute
    the probability of a given conformation at any temperature. 

INPUT.hist INPUT_slice_xx.hist INPUT_par_yy.hist INPUT_par_yy_slice_zz.x 
    - histograms of q at MREMD temperatures.

INPUT.ent INPUT_slice_xx.ent INPUT_par_yy.ent INPUT_par_yy_slice_xx.ent
    - the histogram(s) of energy density.

INPUT.rmsrgy INPUT_par_yy.rmsrgy INPUT_slice_xx.rmsrgy or
             INPUT_par_yy_slice_xx.rmsrgy 
    - the 2D histogram(s) of rmsd from the experimental structure and radius 
      of gyration.

6.2. MAIN INPUT FILE
--------------------

This file has the same structure as the UNRES input file; most of the data are 
input in a keyword-based form (see section 7.1 of UNRES description). The data
are grouped into records, referred to as lines. Each record, except for the 
records that are input in non-keyword based form, can be continued by placing
an ampersand (&) in column 80. Such a format is referred to as the data list
format.

In the following description, the default values are given in parentheses.

6.2.1. General data (data list format)
--------------------------------------

N_ENE (N_ENE_MAX) - the number of energy components

SYM (1) - number of chains with same sequence (for oligomeric proteins only),

HAMIL_REP - if present, Hamiltonian process the results of replica exchange runs
            (replicas with different parameters of the energy function)

NPARMSET (1) - number of energy parameter sets (>1 only for Hamiltonian
    replica exchange simulations)

SEPARATE_PARSET - if present, HREMD was run in a mode such that only temperature
    but not energy-function parameters was exchanged 

IPARMPRINT (1) - number of parameter set with which to construct conformational
    ensembles; important only when HREMD runs are processed

ENE_ONLY - if present, only conformational energies will be calculated and
    printed; no WHAM iteration

EINICHECK (2) - > 0 compare the conformational energies against those stored in
    the coordinate file(s); 1: compare but print only a warning message if
    different; 2: compare and terminate the program if different; 0: don't
    compare.

MAXIT (5000) - maximum number of iterations in solving WHAM equations

ISAMPL (1) - input conformation sampling frequency (e.g., if ISAMPL=5, only
    each 5th conformation will be read)

NSLICE (1) - number of "slices" or "windows" into which each trajectory will
    be partitioned; each slice will be analyzed independently

FIMIN (0.001) - maximum average difference between window free energies
    between the current and the previous iteration

ENSEMBLES (0) - number of conformations (ranked according to probabilities) to
    be output to PDB file at each MREMD temperature; 0 means that no
    conformations will be output. Non-zero values should not be used when NSLICE>1

CLASSIFY - if present, each conformation will be assigned a class, according
to the scheme described in ref [3]

DELTA (0.01) - one dimension bin size of the histogram in q

DELTRMS (0.05) - rms dimension bin size in rms-radius of gyration histograms

DELTRGY (0.05) - radius of gyration bin size in rms-radius of gyration histograms

NQ (1) - number of q's (can be for entire molecule, fragments, and pairs of
    fragments) 

CXFILE - produce the compressed coordinate file with information necessary to
    compute the probabilities of conformations at any temperature

HISTOUT - if present, the histograms of q at MREMD temperatures are
    constructed and printed to main output file

HISTFILE - if present, the histograms are also printed to separate files

ENTFILE - if present, histogram of density of states (entropy) is constructed
    and printed

RMSRGYMAP - if present, 2D histograms of radius of rmsd and radius of gyration at MREMD
    temperatures are constructed and printed

WITH_DIHED_CONSTR - if present, dihedral-angle restraints were imposed in the
    processed MREMD simulations

RESCALE (1) - Choice of the type of temperature dependence of the force field.
0  - no temperature dependence
1  - homographic dependence (not implemented yet with any force field)
2  - hyperbolic tangent dependence [18].

6.2.2 Molecule and energy parameter data
----------------------------------------

6.2.2.1. General information
----------------------------

SCAL14 (0.4) - scale factor of backbone-electrostatic 1,4-interactions

SCALSCP (1.0) - scale factor of SC-p interactions

CUTOFF (7.0) - cut-off on backbone-electrostatic interactions to compute 4-
    and higher-order correlations

DELT_CORR (0.5) - thickness of the distance range in which the energy is
decreased to zero 

ONE_LETTER - if present, the sequence is to be read in 1-letter code,
    otherwise 3-letter code

6.2.2.2. Sequence information
-----------------------------

1st record (keyword-based input):

NRES - number of residues, including the UNRES dummy terminal residues, if present

Next records: amino-acid sequence

3-letter code: Sequence is input in format 20(1X,A3) 

1-letter code: Sequence is input in format 80A1

6.2.2.3. Dihedral angle restraint information
---------------------------------------------

This is the information about dihedral-angle restraints, if any are present.
It is specified only when WITH_DIHED_CONSTR is present in the first record.

1st line: ndih_constr - number of restraints (free format)

2nd line: ftors - force constant (free format)

Each of the following ndih_constr lines:

idih_constr(i),phi0(i),drange(i)  (free format)

idih_constr(i) - the number of the dihedral angle gamma corresponding to the
ith restraint

phi0(i) - center of dihedral-angle restraint

drange(i) - range of flat well (no restraints for phi0(i) +/- drange(i))

6.2.2.4. Disulfide-bridge data
------------------------------

1st line: NS, (ISS(I),I=1,NS)    (free format)

NS - number of cystine residues forming disulfide bridges

ISS(I) - the number of the Ith disulfide-bonding cystine in the sequence

2nd line: NSS, (IHPB(I),JHPB(I),I=1,NSS) (free format)

NSS - number of disulfide bridges

IHPB(I),JHPB(I) - the first and the second residue of ith disulfide link

Because the input is in free format, each line can be split

6.2.3. Energy-term weights and parameter files
----------------------------------------------

There are NPARMSET records specified below.

All items described in this section are input in keyword-based mode.

1st record: Weights for the following energy terms:

WSC (1.0) -    side-chain-side-chain interaction energy

WSCP (1.0) -   side chain-peptide group interaction energy

WELEC (1.0) -  peptide-group-peptide group interaction energy

WEL_LOC (1.0)- third-order backbone-local correlation energy

WCORR (1.0) -  fourth-order backbone-local correlation energy

WCORR5 (1.0) - fifth-order backbone-local correlation energy

WCORR6 (1.0) - sixth-order backbone-local correlation energy

WTURN3 (1.0) - third-order backbone-local correlation energy of pairs of 
               peptide groups separated by a single peptide group

WTURN4 (1.0) - fourth-order backbone-local correlation energy of pairs of 
               peptide groups separated by two peptide groups

WTURN6 (1.0) - sixth-order backbone-local correlation energy for pairs of 
               peptide groups separated by four peptide groups

WBOND (1.0)  - virtual-bond-stretching energy

WANG (1.0) -   virtual-bond-angle-bending energy

WTOR (1.0) -   virtual-bond-torsional energy

WTORD (1.0) -  virtual-bond-double-torsional energy

WSCCOR (1.0) - sequence-specific virtual-bond-torsional energy

WDIHC (0.0)  - dihedral-angle-restraint energy

WHPB (1.0)   - distance-restraint energy

2nd record: Parameter files. If filename is not specified that corresponds to
particular parameters, the respective name from the C-shell script will be
assigned. If no files are to be specified, an empty line must be inserted.

BONDPAR - bond-stretching parameters

THETPAR - backbone virtual-bond-angle-bending parameters

ROTPAR  - side-chain-rotamer parameters

TORPAR  - backbone-torsional parameters

TORDPAR - backbone-double-torsional parameters

FOURIER - backbone-local - backbone-electrostatic correlation parameters

SCCORAR - sequence-specific backbone-torsional parameters (not used at
          present)

SIDEPAR - side-chain-side-chain-interaction parameters

ELEPAR  - backbone-electrostatic-interaction parameters

SCPPAR  - backbone-side-chain-interaction parameters

6.2.4. (M)REMD/Hamiltonian (M)REMD setting specification
--------------------------------------------------------

If HAMIL_REP is present in general data, read the following group of records
only once; otherwise, read for each parameter set (NPARSET times total)

NT (1) - number of temperatures

REPLICA - if present, replicas in temperatures were specified with this parameter set

UMBRELLA - if present, umbrella-sampling was run with this parameter set

READ_ISET - if present, umbrella-sampling-window number is read from the compressed Cartesian
            coordinate (cx) file even if the data are not from umbrella-sampling run(s). 
            ISET is present in the cx files from the present version of UNRES.

Following NT records are for consecutive temperature replicas; each record is
organized as keyword-based input:

TEMP (298.0) - initial temperature of this replica (replicas in MREMD)

FI (0.0) - initial values of the dimensionless free energies for all q-restraint
       windows for this replica (NR values)

KH (100.0) - force constants of q restraints (NR values)

Q0 (0.0d0) - q-restraint centers (NR values)

6.2.5. Information of files from which to read conformations
------------------------------------------------------------

If HAMIL_REP is present in general data, read the following two records
only once; otherwise, read for each parameter set (NPARSET times total)

1st record (keyword-based input):

For temperature replica only ONE record is read; for non-(M)REMD runs, NT
records must be supplied. The records are in keyword-based format.

NFILE_ASC - number of files in ASCII format (UNRES Cartesian coordinate (x)
            files) for current parameter set

NFILE_CX - number of compressed coordinate files (cx files) for current
           parameter set.

NFILE_BIN - number of binary coordinate files (now obsolete because it
           requires initial conversion of ASCII format trajectories into binary format)

It is strongly recommended to use cx files from (M)REMD runs with TRAJ1FILE
option. Multitude of trajectory files which are opened and closed by different
processors might impair file system accessibility. Should you wish to process
trajectories each one of which is stored in a separate file, better collate
the required slices of them first to an x file by using the xdrf2x program
piped to the UNIX cat command.

2nd record:

coordinate file name(s) without extension

6.2.6. Information of reference structure and comparing scheme
-----------------------------------------------------------------

The following records pertain to setting up the classification of conformation
aimed ultimately at obtaining a class numbers. Fragments and pairs of
fragments are specified and compared against those of reference structure in
terms of secondary structure, number of contacts, rmsd, virtual-bond-valence
and dihedral angles, etc. Then the class number is constructed as described in
ref 3. A brief description of comparison procedure is as follows:

1. Elementary fragments usually corresponding to elements of secondary
or supersecondary structure are selected. Based on division into fragments,
levels of structural hierarchy are defined.

2. At level 1, each fragment is checked for agreement with the corresponding
fragment in the native structure. Comparison is carried out at two levels:
the secondary structure agreement and the contact-pattern agreement level. 

At the secondary structure level the secondary structure (helix, strand 
or undefined) in the fragment is compared with that in the native fragment 
in a residue-wise manner. Score 0 is assigned if the structure is different 
in more than 1/3 of the fragment, 1 is assigned otherwise.

The contact-pattern agreement level compares the contacts between the peptide
groups of the backbone of the fragment and the native fragment and also
compares their virtual-bond dihedral angles gamma. It is allowed to shift
the sequence by up to 3 residues to obtain contact pattern match. A score
of 0 is assigned if more than 1/3 of native contacts do not occur or
there is more than 60 deg (usually, but this cutoff can be changed) maximum
difference in gamma. Otherwise score 1 is assigned. 

The total score of a fragment is an octal number consisting of bits
hereafter referred to S (secondary structure) C (contact match) and H
(sHift) (they are in the order HCS). Their values are as follows:

S - 1 native secondary structure; 0 otherwise;
C - 1 native contact pattern; 0 otherwise;
H - 1 contact match obtained without sequence shift 0 otherwise.

For example, octal 7 (111) corresponds to native secondary structure, native
contact pattern, and no need to shift the sequence for contact match;
octal 1 (001) corresponds to native secondary structure only (i.e., nonnative
contact pattern).

3. At level 2, contacts between (i) the peptide groups or (ii) the side chains 
within pairs of fragments are compared. Case (i) holds when we seek contacts
between the strands of a larger beta-sheet formed by two fragments, case (ii)
when we seek the interhelix or helix-beta sheet contacts. Additionally,
the pairs of fragments are compared with their native counterparts by rmsd.
Score 0 is assigned to a pair of fragments, if it has less than 2/3 native 
contacts and too large rmsd (a cut-off of 0.1 A/residue is set), score 1 if
it has enough native contacts and sufficiently low rmsd, but the sequence 
has to be shifted to obtain a match, and score 2, if sufficient match is 
obtained without shift.

4. At level 3 and higher, triads, quadruplets,..., etc. of fragments are
compared in terms of rmsd from their native counterparts (the last level
corresponds to comparing whole molecules). The score (0, 1, or 2) is assigned 
to each composite fragment as in the case of level 2.

5. The TOTAL class number of a structure is a binary number composed of
parts of scores of fragments, fragment pairs, etc. It is illustrated 
on the following example; it is assumed that the molecule has three fragment 
as in the case of 1igd.

level 1      level 2                   level 3
123 123 123||1-2 1-3 2-3 1-2 1-3 2-3 || 1-2-3 | 1-2-3 ||
sss|ccc|hhh|| c   c   c | h   h   h  ||   r   |   h   ||

Bits s, c, and h of level 1 are explained in point 2; bits c and h of level
2 pertain to contact-pattern match and shift; bits r and h of level 3 pertain 
to rmsd match and shift for level 3.

The input is specified as follows:


Program to classify structures

1st record (keyword-based input):

VERBOSE : if present, detailed output in classification (use if you want to 
       fill up the disk)

PDBREF : if present, the reference structure is read from the pdb

BINARY : if present, the class will be output in octal/quaternary/binary format
      for levels 1, 2, and 3, respectively

DONT_MERGE_HELICES : if present, the pieces of helices that contain only
      small breaks of hydrogen-bonding contacts (e.g., a kink) are not merged
      in a larger helix

NLEVEL=n : number of classification levels

n>0 - the fragments for n levels will be defined manually
n<0 - the number of levels is -n and the fragments will be detected automatically

START=n : the number of conformation at which to start

END=n : the number of conformation at which to end

FREQ=n (1) : sampling frequency of conformations; e.g. FREQ=2 means that every
      second conformation will be considered

CUTOFF_UP=x : upper boundary of rmsd cutoff (the value is per 50 residues)

CUTOFF_LOW=x : lower boundary of rmsd cutoff (per 50 residues)

RMSUP_LIM=x : lower absolute boundary of rmsd cutoff (regardless of fragment 
    length)

RMSUPUP_LIM=x : upper absolute boundary of rmsd cutoff (regardless of fragment 
    length)

FRAC_SEC=x (0.66666) the fraction of native secondary structure 
     to consider a fragment native in secondary structure

2nd record:

For nlevel < 0 (automatic fragment assignment):

SPLIT_BET=n (0) : if 1, the hairpins are split into strands and strands are
     considered elementary fragments

ANGCUT_HEL=x (50): cutoff on gamma angle differences from the native for a helical
     fragment

MAXANG_HEL=x (60) : as above but maximum cutoff

ANGCUT_BET=x (90), MAXANG_BET=x (360), ANGCUT_STRAND=xi (90), MAXANG_STRAND=x (360) 
     same but for a hairpin or sheet fragment.

FRAC_MIN=x (0.6666) : minimum fraction of native secondary structure

NC_FRAC_HEL=x (0.5) : fraction of native contacts for a helical fragment

NC_REQ_HEL=x (0) : minimum required number of contacts

NC_FRAC_BET=x (0.5), NC_REQ_BET=x (0) : same for beta sheet fragments

NC_FRAC_PAIR=x (0.3), NC_REQ_PAIR=x (0) : same for pairs of segments

NSHIFT_HEL=n (3), NSHIFT_BET=n (3), NSHIFT_STRAND=n (3), NSHIFT_PAIR=n (3) : 
      allowed sequence shift to match native and compared structure for the 
      respective types of secondary structure
   
RMS_SINGLE=n (0), CONT_SINGLE=n (1), LOCAL_SINGLE=n (1), RMS_PAIR=n (0), 

CONT_PAIR=n (1) : types of criteria in considering the geometry of a fragment
      or pair native; 1 means that the criterion is turned on

For nlevel > 0 (manual assignment):

Level 1:

1st line:

NFRAG=n : number of elementary fragments

Next lines (one group of lines per each fragment):

1st line:

NPIECE=n : number of segments constituting the fragment

ANGCUT, MAXANG, FRAC_MIN, NC_FRAC, NC_REQ : criterial numbers of native-likeness 
      as for automatic classification

LOCAL, ELCONT, SCCONT, RMS : types of criteria implemented, as for automatic 
      classification except that ELECONT and SCCONT mean that electrostatic or
      side-chain contacts are considered, respectively

NPIECE following lines:

IFRAG1=n, IFRAG2=n : the start and end residue of a continuous segment constituting 
      a fragment

Level 2 and higher:

1st line:

NFRAG=n : number of fragments considered at this level

For each fragment the following line is read:

NPIECE=n : number of elementary fragments (as defined at level 1) constituting this
      composite fragment

IPIECE=i1 i2 ... in: the numbers of these fragments 

NC_FRAC, NC_REQ : contact criteria (valid only for level 2)

ELCONT, SCCONT, RMS : as for level 1; note, that for level 3 and higher the only
       criterion of nativelikeness is rms

3rd (for nlevel<0) or following (for n>0) line:

Name of the file with reference structure (e.g., the pdb file with the 
      experimental structure) 

6.3. The structure of the main output file (out)
------------------------------------------------

The initial portion of the main output file, named INPUT.out_POT_000 
contains information of parameter files specified in the C-shell script, 
compilation info, and the UNRES numeric code of the amino-acid sequence.
Subsequently, actual energy-term weights and parameter files are printed.
If lprint was set at .true. in parmread.F, all energy-function
parameters are printed. If REFSTR was specified in the control-data list,
the program then outputs the read reference-structure coordinates and
partition of structure into fragments. 

Subsequently, the information about the number of structures read in and
those that were rejected is printed followed by succinct information form
the iteration process. Finally, the histograms (also output separately to
specific histogram files; see section 6.6) and the data of the dependence of
free energy, energy, heat capacity, and conformational averages on temperature
are printed (these are also output separately to file described in section
6.6).

The output files corresponding to non-master processors
(INPUT.out_POT_xxx where xxx>0 contain only the information up to the
iteration protocol. These files can be deleted right after the run.

6.4. The thermodynamic quantity and ensemble average (thermal) files
-----------------------------------------------------------------

The files INPUT.thermal or INPUT_slice_yy.thermal contain thermodynamic,
ensemble-averaged conformation-dependent quantities and their temperature
derivatives. The structure of a record is as follows:

   T            F              E       q_1...q_n   rmsd     Rgy      Cv       var(q_1)...var(q_n)  var(rmsd)    var(Rgy)  cov(q_1,E)...cov(q_n,E) cov(rmsd,E) cov(Rgy,E)
  298.0      -83.91454     -305.28112   0.30647   6.28347  11.61204 0.70886E+01    0.35393E-02    0.51539E+01    0.57012E+00    0.43802E+00 0.62384E+01    0.33912E+01

where:

T: absolute temperature (in K),

F: free energy at T,

E: average energy at T,

q_1..q_n: ensemble-averaged q values at T (usually only the total q corresponding to whole
           molecule is requested, as in the example above, but the user can specify 
           more than one fragment or pair of fragments for which the q's are 
           calculated, If there's no reference structure, this entry contains
           a 0,

rmsd: ensemble-averaged root mean square deviation at T,

Rgy: ensemble-averaged radius of gyration computed from Calpha coordinates at T,

Cv: heat capacity at T,

var(q_1)...var(q_n): variances of q's at T,

var(rmsd): variance of rmsd at T,

var(Rgy): variance of radius of gyration at T, 

cov(q_1,E)...cov(q_n,E): covariances of q's and energy at T,

cov(rmsd,E): covariance of rmsd and energy at T,

cov(Rgy,E): covariance of radius of gyration and energy at T.

According to Camacho and Thirumalali (Europhys. Lett., 35, 627, 1996), the
maximum of the variance of the radius of gyration corresponds to the collapse
point of a polypeptide chain and the maximum variance of q or rmsd corresponds to 
the midpoint of the transition to the native structure. More precisely, these
points are inflection points in the plots of the respective quantities which,
with temperature-independent force field, are proportional to their covariances 
with energy.

6.5. The conformation summary with classification (stat) files
--------------------------------------------------------------

The stat files (with names INPUT_POT_xxx.stat or
INPUT_POT_sliceyyxxx.stat; where yy is the number of a slice and xxx
is the rank of a processor) contain the output of the classification
of subsequent conformations (equally partitioned between processors). The
files can be concatenated by processor rank to get a summary file. Each line
has the following structure (example values are also provided):

                                      |                           level 1                           |       level 2                |    level3 |
                                      |                                                             |                              |           |
                      whole mol       |            frag1           frag2            frag3       cl1 |         level3               |           |
No      energy    rmsd  q      ang dif|n1n2 n3   rms  q    ang   rms  q    ang   rms  q    ang      | nc1nc2 rms q     rms q    cl2|    rms cl3|class
 9999   -122.42   4.285 0.3751  47.8  |4 10 21   0.6 0.33  16.7  3.6 0.42  56.3  0.7 0.12  16.5 737 | 9  0   1.6 0.20  4.3 0.20 20 | 0   4.0 2 |737.20.2

No - number of conformation

whole mol denotes the characteristics of the whole molecule
q - 1-(Wolynes' q)

level 1, 2, and 3 denote the characteristics computed for the respective fragments
as these levels.

n1, n2, n3 - number of native contacts for a given segment

cl1, cl2, cl3 - group of segment classes for segments at level 1, 2, and 3, respectively

class - total class of the conformation

The octal/quaternary/binary numbers denoting the class for a fragment at level 1, 2, 
and 3, respectively, are described in ref. 3

6.6. The histogram files
------------------------

The histogram file with names INPUT_[par_yy][_slice_xx].hist where xx denotes
the number of the slice and yy denotes the number of the parameter if 
SEPARATE_PARSET was specified in input contain histograms of q at replica
temperatures and energy-parameter sets; with SEPARATE_PARSET histograms 
corresponding to subsequent parameter sets are saved in files with par_yy
infixes. The histograms are multidimensional if q is a vector (usually,
however, q corresponds to the entire molecule and, consequently, the
histograms are one-dimensional). The histogram files are printed if histfile
and histout was specified in the control data record.

Each line of a histogram file corresponds to a given (multidimensional) bin in 
q contains the following:

q_1,...,q_n at a given bin (format f6.3 for each)

histogram values for subsequent replica temperatures (format e20.10 for each)

iparm (the number of parameter set; format i5)

If SEPARATE_PARSET was not specified, the entries corresponding to each
parameter follow one another.

The state density (microcanonical entropy) is printed to file(s)
INPUT[_slice_xx].ent. Each line contains the left boundary of the energy
bin and ln(state density) followed by " ent" string. At present, the state
density is calculated correctly only if one energy-parameter set is used.

6.7. The rmsd-radius of gyration potential of mean force files
------------------------------------------

These files with names INPUT[_par_yy][_slice_xx].rmsrgy contain the
two-dimensional potentials of mean force in rmsd and radius of gyration
at all replica-exchange temperatures and for all energy-parameter sets. 
A line contains the left boundaries of the radius of gyration - rmsd bin
(radius of gyration first) (format 2f8.2) and the PMF values at all
replica-exchange temperatures (e14.5), followed by the number of the parameter
set. With SEPARATE_PARSET, the PMFs corresponding to different parameter sets
are printed to separate files.

6.8. The PDB files
------------------

The PDB files with names INPUT_[slice_xx_]Tyyy.pdb, where Tyyy specifies 
a given replica temperature contain the conformations whose probabilities at
replica temperature T sum to 0.99, after sorting the conformations by
probabilities in descending order. The PDB files follow the standard format; 
see ftp://ftp.wwpdb.org/pub/pdb/doc/format_descriptions/Format_v33_Letter.pdf. 
For single-chain proteins, an example is as follows:

REMARK CONF    9059 TEMPERATURE  330.0 RMS   8.86
REMARK DIMENSIONLESS FREE ENERGY   -1.12726E+02
REMARK ENERGY   -2.22574E+01 ENTROPY   -7.87818E+01
ATOM      1  CA  VAL     1       8.480   5.714 -34.044
ATOM      2  CB  VAL     1       9.803   5.201 -33.968
ATOM      3  CA  ASP     2       8.284   2.028 -34.925
ATOM      4  CB  ASP     2       7.460   0.983 -33.832
.
.
.
ATOM    115  CA  LYS    58      28.446  -3.448 -12.936
ATOM    116  CB  LYS    58      26.613  -4.175 -14.514
TER
CONECT    1    3    2
.
.
.
CONECT  113  115  114
CONECT  115  116

where

CONF is the number of the conformation from the processed slice of MREMD
trajectories

TEMPERATURE is the replica temperature

RMS is the Calpha rmsd from the reference (experimental) structure.

DIMENSIONLESS FREE ENERGY is -log(probability) (equation 14 of ref 2)
for the conformation at this replica temperature calculated by WHAM.

ENERGY is the UNRES energy of the conformation at the replica temperature
(note that UNRES energy is in general temperature dependent).

ENTROPY is the omega of equation 15 of ref 2 of the conformation 

In the ATOM entries, CA denotes a Calpha  atom and CB denotes UNRES side-chain
atom. The CONECT entries specify the Calpha(i)-Calpha(i-1),
Calpha(i)-Calpha(i+1) and Calpha(i)-SC(i) links.

The PDB files generated for oligomeric proteins are similar except that
chains are separated with TER and molecules with ENDMDL records and chain
identifiers are included. An example is as follows:

REMARK CONF     765 TEMPERATURE  301.0 RMS  11.89
REMARK DIMENSIONLESS FREE ENERGY   -4.48514E+02
REMARK ENERGY   -3.58633E+02 ENTROPY    1.51120E+02
ATOM      1  CA  GLY A   1      -0.736  11.305  24.600
ATOM      2  CA  TYR A   2      -3.184   9.928  21.998
ATOM      3  CB  TYR A   2      -1.474  10.815  20.433
.
.
.
ATOM     40  CB  MET A  21      -4.033  -2.913  27.189
ATOM     41  CA  GLY A  22      -5.795 -10.240  27.249
TER
ATOM     42  CA  GLY B   1       6.750  -6.905  19.263
ATOM     43  CA  TYR B   2       5.667  -4.681  16.362
.
.
.
ATOM    163  CB  MET D  21       4.439  12.326  -4.950
ATOM    164  CA  GLY D  22      10.096  14.370  -9.301
TER
CONECT    1    2
CONECT    2    4    3
.
.
.
CONECT   39   41   40
CONECT   42   43
.
.
.
CONECT  162  164  163
ENDMDL

6.8. The compressed Cartesian coordinates (cx) files
----------------------------------------------------

These files contain compressed data in the Europort Data Compression XDRF
library format written by Dr. F. van Hoesel, Groeningen University 
(http://hpcv100.rc.rug.nl/xdrfman.html). The files are written
by the cxwrite subroutine. The resulting cx file contains the omega
factors to compute probabilities of conformations at any temperature
and any energy-function parameters if Hamiltonian replica exchange was
performed in the preceding UNRES run. The files have general names 
INPUT[_par_yy][_slice_xx].cx where xx is slice number and yy is parameter-set
number.

The items written to the cx file are as follows (the precision is 5
significant digits):

1) Cartesian coordinates of Calpha and SC sites
2) nss (number of disulfide bonds)
3) if nss > 0: 
   a) ihpb  (first residue of a disulfide link)
   b) jhpb  (second residue of a disulfide link)
4) UNRES energy at that replica temperature that the conformation was at
   snapshot-recording time,
5) ln(omega) of eq 15 of ref 2,
6) Calpha rmsd
7) conformation class number (0 if CLASSIFY was not specified).

7. SUPPORT
----------

   Dr. Adam Liwo
   Faculty of Chemistry, University of Gdansk
   ul. Sobieskiego 18, 80-952 Gdansk Poland.
   phone: +48 58 523 5430
   fax: +48 58 523 5472
   e-mail: adam@chem.univ.gda.pl

   Dr. Cezary Czaplewski
   Faculty of Chemistry, University of Gdansk
   ul. Sobieskiego 18, 80-952 Gdansk Poland.
   phone: +48 58 523 5430
   fax: +48 58 523 5472
   e-mail: czarek@chem.univ.gda.pl

Prepared by Adam Liwo, 02/19/12