The directory contains examples which can be used
to test features of GMCT and provide examples of
the necessary program input files and settings.
The examples may also be useful as primers for
setting up own calculations.

For the sake of simplicity, the system is always
named by the single character "d".


DTPA (DTPA_titra, DTPA_coop & DTPA_macro)
=========================================

The DTPA example from the program manual
(diethylen-triamine-pentaacetat). See the
program manual for explanation.
The parameters are taken from:
Onufriev, A; Case, D. A ; Ullmann, G. M.,
Biochemistry, 2001, 40, 3413-3419.

DTPA_titra
----------

The system is set up for the calculation of site
titration curves. A trajectory is written with
the program gmct and analyzed with the program
trajectory-analysis. The occupation probabilities
are already given out by gmct. An error analysis
accounting for autocorrelation of the samples is
carried out by trajectory-analysis. The
probabilities and their statistical uncertainties
are found in the directory Prob. The autocorrelation
functions and autocorrelation lengths are written
to the directory ACorr. The covariances between
the occupancies are written to the directory Corr.

An analytical calculation is started
with the command:
GMCT_PATH/smt d > d.smt.out

The Metropolis Monte Carlo simulation
is started with the command:
GMCT_PATH/gmct d > d.gmct.out

The trajectories written by gmct are
analyzed with the command:
GMCT_PATH/traj-analysis d > d.traj-analysis.out

DTPA_coop
---------

The system is set up for the calculation of
binding free energies and cooperativity free
energies as functions of the pH value. The
preset options select the free energy
perturbation method for the calculation.
You can of course also change the settings
to try one of the other free energy calculation
methods. The reaction free energies are found
in the directory Grct. The cooperativity free
energies between pairs of reactions are found
in the directory Coop/2. The covariances
between pairs of product state occupancies
are written to the directory Corr.

An analytical calculation
is started with the command:
GMCT_PATH/smt-coop d > d.smt-coop.out

In principle, the calculation can also be done
from equilibrium simulation data, but statistically
meaningful results for the free energies can not be
obtained for all pH values from trajectories of
affordable length:

The Metropolis Monte Carlo simulation
is started with the command:
GMCT_PATH/gmct d > d.gmct.out

The trajectories written by gmct are
analyzed with the command:
GMCT_PATH/traj-analysis-coop d > d.traj-analysis-coop.out

traj-analysis-coop and smt-coop calculate also
covariances between pairs between the
product states of pairs of the generated stepwise
binding reactions.

DTPA_macro
----------

An example for an analysis of the macroscopic
binding behavior and the global thermodynamic
functions of the system with the Wang-Landau
Monte Carlo method. The DTPA model is a tough
case for the method. The calculation converges
slowly because of the low number of microstates
available to the system.

Note that the histogram flatness criterion is
not too reliable as convergence indicator.
With a large enough number of MC scans a very
low flatness criterion can be obtained even if
the modification factor is still large. It is
thus advisable to check also convergence of the
relative number of states in the bins with the
largest and smallest number of states printed
as ln(gmax/gmin).

The system is set up for the calculation of the
occupation probabilities of macroscopic
protonation and reduction states, the macroscopic
binding free energies, the partition function of
the system and the global thermodynamic functions
free energy, enthalpy and entropy as functions of
the pH value. The global thermodynamic functions
and the binding free energies are found in the
directory Gmacro. The distribution of the number
of states as function of the microstate energy
is found in the directory DOS. 

An analytical calculation is started with:
GMCT_PATH/macro-smt-coop d > d.macro-smt-coop.out

The Wang-Landau Monte Carlo simulation
is started with the command:
GMCT_PATH/macro-mct-coop d > d.macro-mct-coop.out


Hen-egg white lysozyme (HEWL_titra, HEWL_coop & HEWL_macro)
===========================================================

Hen-egg white lysozyme (HEWL is often used as
model system for protonation state calculations.
See for example:

Bashford, D.; and Karplus, M.; Biochemistry,
1990, 29, 10219-10225

Webb, H.; Tynan-Connolly, B. M.; Lee, G. M.;
Farrell, D.; O'Meara, F.; Sndergaard, C. R.;
Teilum, K.; Hewage, C.; McIntosh, L. P. and
Nielsen, J. E., Proteins 2010, 79, 685-702.

The example system was prepared from PDB 4LZT
with our in-house program GCEM based on the MEAD
library (Ullmann, R. T., to be published, see
also the program manual for explanation of the
energy terms and the HEWL example for GCEM that
is included in our MEAD distribution). A single
conformer was used for each site as found in the
crystal structure. A detailed model was chosen
for all titratable sites with explicit positions
for each of the dissociable hydrogen atoms. All
aspartate, glutamate, and the C-terminus are
modeled with four tautomers for the protonated
form (cis and trans configuration of the
titratable proton at each carboxyl oxygen) and
a single deprotonated form. All lysine residues
and the N-terminus are modeled with a single
instance for the protonated form and three
rotamers for the deprotonated form (loss of
each of the three hydrogen atoms of the amino
group). The arginine residues are modeled with
a single protonated instance and five tautomers
of the deprotonated form (loss of each hydrogen
atom of the guanidino group). The tyrosine
residues are modeled with a single instance for
the deprotonated form and three rotamers for
the dissociable hydrogen atom of the hydroxyl
group. The histidine residue is modeled with
four instances corresponding to the fully
deprotonated imidazolate form, the singly
protonated delta and epsilon tautomers and the
fully protonated imidazolium form. This detailed
model yields pK1/2 values that compare very
well to the experimental ones (see for example
Nielsen2010). The increased detail relative
to traditional two-state models for proto-
nation state calculation is of course paid
with a somewhat higher computational cost.

The same set of calculations as for the DTPA
example was prepared (see above and comments
in the setup file located in each directory).
Analytical calculations are no longer possible
at this system size. The Monte Carlo simulation
methods provide highly efficient ways to
obtain accurate and precise results for larger
models. The Hybrid method is an approximate,
semianalytic method to calculate observables.
This method can be much more efficient than
the Monte Carlo simulation methods if the system
does not contain too many strongly interacting
sites, as is the case for HEWL. In strongly
coupled systems, it can be necessary to use
large clusters for fairly accurate results.
In such cases, the Monte Carlo method is much
more efficient than the Hybrid method. The
same finding applies for systems with several
strongly interacting sites which are modeled
with many instances.

HEWL_titra
----------

The system is set up for the calculation of
site titration curves.

The Metropolis Monte Carlo simulation
is started with the command:
GMCT_PATH/gmct d > d.gmct.out

The trajectories written by gmct are
analyzed with the command:
GMCT_PATH/traj-analysis d > d.traj-analysis.out

pK(1/2) values can be extracted from the gmct output
and written to my_curves/bind_half/ with
MEAD_PATH/apps/get-curves-gmct/get-curves-gmct d

The calculation can take a while
(25 min on an Intel Core i5 430M).

Alternatively, we can use the approximate
Hybrid method to compute the titration curves
and compare the results to those of the rigor-
ous Monte Carlo simulation. The maximum cluster
site (clustnsitemax) is set to 5. For HEWL,
even the Tanford-Roxby limit of a single site
per cluster works quite well. This is because
most sites are highly solvent exposed, resulting
in a reduction of the electrostatic site-site
interactions due to solvent screening.

In strongly coupled systems, it can be necessary
to use large clusters for fairly accurate results.
In such cases, the Monte Carlo method is much
more efficient than the Hybrid method.
The same finding applies for systems, with multiple
interacting sites which are modeled with more than
a few instances.

The approximate Hybrid calculation is started
with the command:
GMCT_PATH/hybrid d > d.hybrid.out

HEWL_grct
---------

The system is set up for the calculation of the
protonation free energies of the two catalytic
residues in the active site of the enzyme
(Glu-35 and Asp-52) as functions of the pH
value. The preset options select the free
energy perturbation method for the calculations.
You can of course also change the settings to
try one of the other free energy calculation
methods (see the manual).

The protonation free energies for the other
sites can be added by removing the comment
signs in front of the corresponding definitions
in the file d.rnames.The calculation of the
cooperativity free energies between all pairs
of protonation reactions can be triggered by
removing the comment sign in front of the option
icoop in the setup file (takes a long time if
many sites are included due to the combinatorics
giving a maximum of 465 pairs of protonation
reactions for the 31 protonatable sites). For
the calculation of cooperativity free energies
you might want to increase the number of
chimeric intermediates (nifep) and the number
of MC scans (nmcfull). Also the use of a lower
statistical error tolerance (tolg_stderr) might
be reasonable for this purpose.

It can be seen that the protonation free energy
of the active site residues deviates markedly
from that of isolated sites in solution. The
protonation free energies of Asp-52 and Glu-35
are not simple linear functions of the pH
value and do not exhibit the constant slope of
RT ln[10] / pH unit that would be expected for
a single isolated site. The nearly constant
difference between the protonation free energies
of the catalytic residues might be important
for allowing the enzyme to perform its function
over a broad range of pH values. Studying the
protonation free energy of the residues thus
provides information that is not immediately 
obvious from their titration curves. Similar
findings apply to other enzymes.
See for example:

Bombarda, E.; Ullmann, G. M., J. Phys. Chem. B,
2010, 114, 1994-2003.

The free energy calculations are started
with the command:
GMCT_PATH/mct-coop d > d.mct-coop.out

The calculation can take some time
(1 hour 40 min on an Intel Core i5 430M).

HEWL_macro
----------

An example for an analysis of the macroscopic
protonation probabilities and the global
thermodynamic functions of the system with the
Wang-Landau Monte Carlo method.

The global thermodynamic functions are found
in the directory Gmacro. The macroscopic
binding probabilities are printed to stdout.

The Wang-Landau Monte Carlo simulation
is started with the command:
GMCT_PATH/macro-mct-coop d > d.macro-mct-coop.out

The calculation can take a while
(1 hour 40 min on an Intel Core i7 860).


Pseudomonas aeruginosa azurin (PaAz_coop & PaAz_macro)
======================================================

The Pseudomonas aeruginosa azurin example
from the program manual. The system provides
a motivating example for the utility of GMCT
in studying the coupling of the binding
equilibria of different ligands and of the
conformational equilibria of the receptor.

A detailed study and description of the
system can be found in:

Ullmann, R. T. and Ullmann, G. M., 2011
J. Phys. Chem. B., 115, 10346-10359,
dx.doi.org/10.1021/jp204644h.

This example system is modeled in greater
detail than the HEWL example and makes more
use of the capabilities of GMCT. Sidechain
rotamers are included for surface exposed
flexible residues. The protonation forms of
the titratable sites and hydroxyl rotamers
are modeled in the same way as for HEWL.
The conformations of the peptide flip region
are modeled explicitly. Two ligand types
with distinct chemical potentials appear in
this example (proton and electron).

PaAz_coop
---------

The system is set up for the calculation of
the reduction free energy of the copper
center, the protonation free energy of His-35
and the cooperativity free energy between the
two reactions at a pH value of 7.0. The
preset options select the free energy
perturbation method for the calculation.
You can of course also change the settings to
try one of the other free energy calculation
methods. The reaction free energies are found
in the directory Grct. The cooperativity free
energies are found in the directory Coop/2.

The free energy calculations are started
with the command:
GMCT_PATH/mct-coop d > d.mct-coop.out

The calculation can take a while
(20 min on an Intel Core i5 430M).

An analytical calculation is not possible
due to the relatively large system size.

The calculation can in principle also be done
using equilibrium simulation data, but very
long simulations would be required to obtain
statistically meaningful results:

The Metropolis Monte Carlo simulation
is started with the command:
GMCT_PATH/gmct d > d.gmct.out

The trajectories written by gmct are
analyzed with the command:
GMCT_PATH/traj-analysis-coop d > d.traj-analysis-coop.out

PaAz_macro
----------

An example for an analysis of macroscopic
binding behavior and global thermodynamic
functions of the system with the Wang-Landau
Monte Carlo method. The system is set up for
the calculation of the occupation
probabilities of macroscopic protonation and
reduction states, the macroscopic binding
free energies, the partition function of the
system and the global thermodynamic functions
free energy, enthalpy and entropy at a pH
value of 7.0.

An analytical calculation is not possible
due to the relatively large system size.

The Wang-Landau Monte Carlo simulation
is started with the command:
GMCT_PATH/macro-mct-coop d > d.macro-mct-coop.out

The calculation can take a while
(5 hours 20 min on an Intel Core i7 860).


Bennett's example system (bennett_setup_A & bennett_setup_B)
============================================================

The example system from Charles Bennett's paper that
introduced the Bennett acceptance ratio (BAR) method:
Bennett, C. H.; J. Comp. Phys., 1976, 22, 245-268,
dx.doi.org/10.1016/0021-9991(76)90078-4

The two different, equivalent setups are described
in our paper on a generalized formulation of the
free energy perturbation (FEP) theory:
Ullmann, R. T.; Ullmann, G. M., J. Phys. Chem. B,
2011, 115, 507-521, dx.doi.org/10.1021/jp1093838

A free energy calculation with a FEP simulation
is started with the command:
GMCT_PATH/mct-coop d > d.mct-coop.out

An analytical free energy calculation
is started with the command:
GMCT_PATH/smt-coop d > d.smt-coop.out


Alchemical FEP example (relative_binding_free_energy)
=====================================================

A greatly simplified, discretized representation
of an alchemical binding free energy calculation.
The example appeared also in our FEP paper.

The system is set up for the calculation of the
relative binding free energy of the two ligands,
the protonation free energy of the titratable
site and of the cooperativity free energy between
the exchange of the ligands and the protonation
of the titratable site.

This example system is very challenging for the
free energy simulation methods, although it is
very simple to calculate the free energies and
cooperativity free energy analytically.

A free energy calculation with a FEP simulation
is started with the command:
GMCT_PATH/mct-coop d > d.mct-coop.out

An analytical free energy calculation
is started with the command:
GMCT_PATH/smt-coop d > d.smt-coop.out


Accounting for electrostatic and electrochemical
trans-membrane potential differences (2_site_bR)
================================================

A reduced model of the proposed "back-pressure valve"
function of Asp-85 and Asp-115 in bacteriorhodopsin
(bR) from Halobacterium salinarum. The residues
Asp-85 and Asp-115 are included for one monomer.
The intrinsic energies, membrane potential
interaction energies and the site-site interaction
energy are taken from a full continuum electrostatic
model of bR (see also the program manual for
explanation of the energy terms and the bR example
for GCEM included in our MEAD distribution).

See:
Calimet, N.; Ullmann, G. M., J. Mol. Biol.,
2004, 339, 571-589
and
Bombarda, E.; Becker, T.; Ullmann, G. M.,
J. Am. Chem. Soc., 2006, 128, 12129-12139.

The example illustrates the necessary settings for
the treatment of trans-membrane electrochemical
potential differences and of electrostatic trans-
membrane potentials. The electrostatic transmembrane
potential calculations are implemented in the MEAD
program my_membpot_solver. It turned out, that our
implementation is very similar to that in the CHARMM
module PBEQ, although the corresponding publication
contains an alternative, equivalent formulation.
See the MEAD documentation on gwdg.de/~tullman
for the underlying theory of our implementation and
the original paper of Benoit Roux on the topic:

Roux, Benoit; Biophys. J., 1997, 73, 2980-2989.

An analytical calculation of equilibrium
occupation probabilities is started with
the command:
GMCT_PATH/smt d > d.smt.out

A Metropolis Monte Carlo simulation for
the calculation of equilibrium occupation 
probabilities is started with the command:
GMCT_PATH/gmct d

An analytical free energy calculation
is started with the command:
GMCT_PATH/smt-coop d > d.smt-coop.out

A free energy calculation with a FEP simulation
is started with the command:
GMCT_PATH/mct-coop d > d.mct-coop.out
