Purpose of MEAD and our extensions

MEAD
consists of a library of C++ objects and some applications that
use these objects for modeling electrostatic properties of molecules.
MEAD
is an acronym for macroscopic electrostatics with atomic detail.
That is, macroscopic continuum electrostatics is applied at a molecular scale,
partitioning the system in regions with different dielectric constants
(molecule, solvent, membrane, ...). Solvent regions can contain mobile ions.
The electrostatic potential is computed according to the
linearized Poisson-Boltzmann equation.
Despite its simplicity, this model can be used to compute molecular properties
with very good accuracy (pKa values, binding constants,
electrostatic solvation energies, ...).

Our extensions allow a detailed modeling of (biological) macromolecules
including the possibility to account for a membrane environment with
an electrostatic trans-membrane potential.
Further extensions allow the visualization of electrostatic potentials,
charge distributions, electrolyte distributions and dielectric distributions.
The data can be given out as three-dimensional volumetric data in OpenDX format,
or within a cut plane or along a line in ASCII format for plotting with your favorite software.
The figure below shows the trans-membrane potential across the ammonia transporter Amt-1
from Archaeoglobus fulgidus as an example application for the visualization extensions.

The trans-membrane potential distribution across the ammonium transporter Amt-1.
The figure was taken from
→ this paper.
The structure of Amt-1 is described in Andrade, Susana L. A.; Dickmanns, Antje; Ficner, Ralph and Einsle, Oliver,
2005, PNAS, 102, 14994-14999.
The extracellular side is shown at the top and the intracellular side at the bottom.
The horizontal black lines indicate the boundaries of the membrane core and headgroup regions.
The dark outer contour of Amt-1 denotes a projection of the solvent accessible surface of
a transporter trimer into a plane perpendicular to the membrane.
The lighter inner contour shows a projection of a slice of Amt-1 of 5 A thickness
into the same projection plane.
The slice shows the trans-membrane pore including putative ammonia positions,
the twin-histidine motive and the Phe-gate.
a) The trans-membrane potential is plotted in a slice plane cutting through the
transporter's trans-membrane pore.
The potential is projected into a plane perpendicular to the membrane, while
the slice plane is slightly tilted relative to the membrane normal to follow
the course of the trans-membrane pore.
The values at the white contours denote the fraction of the trans-membrane potential
at the respective coordinate.
It can be seen that the membrane potential distribution within the protein
does not show a simple linear dependency on the z coordinate.
b) The mobile source charge distribution of the trans-membrane potential
in the same slice and projection planes as in a).
Darker red or blue shades denote higher negative or positive charge
density, respectively.
It can be seen that most of the unbalanced charge is concentrated
close to the membrane and in the depressions of the protein surface.

A lipid membrane can be represented by three dielectric slabs that
model the hydrophobic, ion-inaccessible core region and the hydrophilic headgroup regions
that can be penetrated by mobile ions.
Regions within the membrane boundaries that are not part of the protein or the membrane
can be specified (e.g., water filled protein cavities or a channel through the membrane).

The protein interior can comprise different dielectric regions that can be used,
for example, to account for
protein regions that are modeled classically and such that are modeled with
a quantum chemical approach.
In addition there can be regions that model protein cavities or channels.

The solvent phase(s) are modeled by separate dielectric regions that can also
contain mobile ions.

Documentation

The documentation for the original
MEAD package
is the file README found in the root directory of the distribution.
Below, you can find the information from this README file
and information about our changes and additional application programs.
Example calculations with all necessary input files can be found
in the subdirectory "examples" of the distribution.

The
MEAD
library is best explored by looking at the source code.
The best starting point may be to look at some of the simpler
programs like potential or solvate and then one of the
solver programs (my_xyz_solver).
Don't start with multiflex or
GCEM,
because these programs and the objects used by them
are too complex for a start.
A brief outline of the design of the
MEAD
library can be found in → this paper.

The calculation of binding properties with a continuum electrostatics/molecular mechanics model
is described in a short version on the page about
→
GMCT &
GCEM.
A more detailed description is given in
→ this thesis.
Parts of this description are included in the user manual of
GMCT
(found in the directory doc of the distribution)
and the → corresponding paper.
See → this paper for a review of
titration calculations with continuum electrostatics within a classical two-state model.

General information on MEAD usage

Here, you find information on general options and input files common to all MEAD programs

show more detail

program usage

Program settings are in most cases specified with command line options.
Input files are described below.
The last command line arguments specify the name of the molecule.
The program is called with:

program_dir/program {options} {molname(s)}

where molname(s) are the names of the molecule(s) for which the calculation
will be done.
molname(s) are used as prefix for some input files.

thickness of an ion-exclusion layer (Stern layer) for determining the ion-inaccessible volume

-ionicstr float

0.0

ionic strength in mol/l

-kBolt float

5.984e-06

Boltzmann constant in units of \(\mathrm{e}_{\circ}^{2}/\left(\unicode[Times,serif]{xC5}\mathrm{K}\right)\) (elementary charges squared per Ângström and Kelvin)
This constant must be adjusted if other units of charge, length or temperature are used.

-econv float

332.063202

conversion factor from \(\mathrm{e}_{\circ}^{2}/\left(\unicode[Times,serif]{xC5}\mathrm{K}\right)\) (elementary charges squared per Ângström and Kelvin)
to \(\mathrm{kcal}/\mathrm{mol}\)
This constant must be adjusted if units of energy and electrostatic potentials other than
\(\mathrm{kcal}/\mathrm{mol}\) and \(\mathrm{kcal}/\left(\mathrm{mol}\right)\mathrm{e}_{\circ}\)
are needed in the output.

-conconv float

6.022214e-04

conversion factor from \(\mathrm{mol}/\mathrm{l}\) (moles per liter)
to \(1/\unicode[Times,serif]{xC5}^{3}\) (particles per cubic Å)
This constant must be adjusted if other units of length or concentration are used.

-epssave_oldway

flag that triggers the old style of averaging the dielectric constant between grid points
of differing dielectric constants.
Between two potential lattice points in the
finite difference method. The new way involves
inverse averaging and is similar to a proposal
by McCammon. The old way is a simple mean.
The new way is significantly more accurate
the option to do it the old way is only provided
for the sake of reproducing old experimental
results. It is not recommended otherwise.

-converge_oldway

Revert to the old way of testing for convergence
of the successive over-relaxation (SOR) method of
solving the finite difference representation of the linearized Poisson Boltzmann equation.
The new method gives improved long range accuracy
for large lattices, but at a sometimes substantial
computational cost. See the NEWS file for further
discussion.

-blab1-blab2-blab3

flag that controls the verbosity of the programs while writing to stdout,
specifying no blab flag at all is least verbose.
In the latter case, only essential output is written to stdout.
Writing to output files is not affected.

visualization options (my_xyz_solver):

option andexpected data type

default setting(if any)

meaning

-write_pot

flag that triggers the output of the electrostatic potential

-write_rho

flag that triggers the output of the charge distribution

-write_eps

flag that triggers the output of the dielectric constant distribution

-write_ely

flag that triggers the output of the ionic strength distribution

-x float

x coordinate of the grid center used for the OpenDX volumetric data file.

-y float

y coordinate of the grid center used for the OpenDX volumetric data file.

-z float

z coordinate of the grid center used for the OpenDX volumetric data file.

-space float

grid spacing used for the OpenDX volumetric data file.

-count int {int int}

edge length of the grid used for the OpenDX volumetric data output in grid points.
If one number is specified, all edges have the same length.
If three numbers are specified, they correspond to the edge lengths in
x, y and z direction.
For the output of the electrostatic potential,
the grid must be entirely covered by the grid used in
the calculation of the electrostatic potential
(specified in the .ogm or .mgm file).

-gmt 7 floats 1 int {6 more floats}

Specification of this option triggers the output of the quantities requested by
write_pot, write_rho, write_eps and write_ely
as two-dimensional curves for plotting with the
GMT
program suite or other plotting software.
The arguments are
xcenter, ycenter, zcenter, xnormal, ynormal, znormal, grid spacing, number of grid points
{xcenter2, ycenter2, zcenter2, xnormal2, ynormal2, znormal2}
The first three numbers define the x, y and z coordinates of the grid / plane center.
The next three numbers define components of the normal vector of the plane.
The following integral number defines the edge length of the curve in grid points.
Optionally, six more values can be specified to define the base point and the normal vector
of a second plane.
In this case, the values of the requested functions are
calculated in this second plane then and projected onto the first plane.
This feature can for example be useful for a trans-membrane pore that is not
perfectly normal to the membrane plane as shown in the figure on the top of the page.

input files

name.pqr
contains a molecule structure in a format similar to that of a PDB file but
with atomic partial charges and radii in the occupancy and
B factor columns, respectively.
More specifically, lines beginning with either "ATOM" or "HETATM" (no leading spaces)
are interpreted as a set of tokens separated by one or more spaces or TAB characters.
Note that the .pqr format does not support some PDB features such
as a altLoc fields, and a one character chainID between
resName and resSeq. Doing so would break the whitespace
separated tokens convention that allows for easy processing
with perl scripts, etc.
Instead we optionally allow for additional digits at the end of the line
to specify global conformation, chain, residue and instance,
where instance is a certain form of a site or residue
(see →the page on
GMCT
and
GCEM.
If you have a PDB file and need to generate a PQR file, this
implies making some choices about the charges and radii. This
is similar to making a choice about what force field to use in
an MD simulation. MEAD per se, doesn't make the choice for
you. However, in the utilities subdir are some tools that may
be useful if you want to use CHARMM or PARSE parameters. The
Amber program suite comes with a program, ambpdb, that allows
you to generate a PDB or PQR file given Amber format files.
Another option is the program
pdb2pqr.
format:

The contents of the first two columns are ignored.
The last four columns are optional
(used, e.g., by GCEM).
atname is the atom name.
resname is the residue name.
resnum is the residue number.
x,y,z are the x, y, and z coordinates of the atom (floating point numbers).
charge is the atomic partial charge of the atom.
radius is the atom radius.
confid, chainid siteid and instid are integer numbers that identify the global conformation,
the polymer chain, the site and the instance to which the atom belongs, respectively.

name.ogm and name.mgm
definine the cubic grids for the computation of electrostatic potentials
with the finite difference method.
format:

grid_center is the grid center and can be given as three floating point numbers
specifying the coordinates directly or as centering style.
Optionally, the coordinate values can be enclosed by parenthesis and separated by commas.
For the centering style, three options are available.
ON_ORIGIN specifies that the grid center is placed on the origin of the coordinate system
(0.0, 0.0, 0.0).
ON_GEOM_CENT specifies that the grid center is to be placed on the geometric center
of the molecule / receptor or site denoted by name.
ON_CENT_OF_INTEREST specifies that the grid center is to be placed on the geometric center
of a site of interest.
grid_spacing is the spacing between two grid points in Å.
grid_dimension is the edge length of the cubic grid in grid points
and must be an odd integer.
Each grid must be smaller than the previous grid, and normally it will also
have a smaller grid spacing.
As a rule of thumb, the outermost grid should cover at least the
charge distribution for which the electrostatic potential is to be
computed plus 15 Å in each direction.
The outermost grid also has to include all
points at which the electrostatic potential is to be computed.

name.fpt
is a file containing the coordinates at which the electrostatic potential shall be calculated
(traditional format)
or a file containing the atomic partial charges and their coordinates
(extended format I for the the application programs named my_xyz_solver)
or a file containing the atomic partial charges and their coordinates for each instance
of each site found in a receptor
(extended format II for the the application programs named my_xyz_solver)

where coordinate consists of three floating point values that denote the x, y and z
coordinates of a point respectively.
Optionally, the coordinate values can be enclosed by parenthesis and separated by commas.
The line breaks can be substituted with any number of whitespace characters.

extended format I (for command line option -pf):

coordinate_1 charge_1
...
coordinate_N charge_N

coordinate consists of three floating point values that denote the x, y and z
coordinates of a point respectively.
Optionally, the coordinate values can be enclosed by parenthesis and separated by commas.
charge is the atomic partial charge at the preceding coordinate.
The line breaks can be substituted with any number of whitespace characters.

site and instance denote the site and instance to whose
charge distribution the following atomic partial charge belongs.
coordinate consists of three floating point values that denote the x, y and z
coordinates of a point respectively.
Optionally, the coordinate values can be enclosed by parenthesis and separated by commas.
charge is the atomic partial charge at the preceding coordinate.
The line breaks can be substituted with any number of whitespace characters.

name.potat
This is a binary file produced by some programs,
to avoid costly recalculation of electrostatic potentials.
It contains the potential at each atom of name
generated by some set of charges. Variations of
name may denote charge states, sites, conformers,
solvated or vacuum or uniform dielectric environments,
depending on the application. Atomic coordinates and radii
and the generating charges are also included for the sake of
consistency checking.
These files allow multiflex, etc. to avoid unnecessary recalculations
when all you want to do is add or alter some site,
but you must be careful about which .potat files you keep.
Specify the -blab2 flag for a blow-by-blow account of attempts
to read and write .potat files in multiflex.

GCEM

→
GCEM
is a program for the automated preparation of the necessary input
for
→
GMCT
from a continuum electrostatics / molecular mechanics model.

show more detail

GCEM
is a program for the automated preparation of the necessary input
for
GMCT
from a continuum electrostatics / molecular mechanics model.
Details about the applications and the underlying model of
GCEM
can be found on a
→ separate page.
Examples for the usage of the program can be found in the directory examples/gcem.

program usage

GCEM
computes the energy terms used in the microstate energy function of
→
GMCT
using a continuum electrostatics / molecular mechanics model.
A schematic view of the model is shown on
→ this page
(where not all features need to be present).
GCEM
has a number of new features as compared to multiflex
and is based on a generalized formulation of binding theory,
which offers a wider application range and more transparency.

GCEM
considers one global conformation at a time.
That is, a separate calculation needs to be set up for each
global conformation.

Program settings are specified with command line options.
Input files are described below.
The last command line argument is the name of the molecule.
The program is called with:

program_dir/gcem {options} molname

where molname is the name of the molecule / receptor.
molname is used as prefix for some input files.

Each
GCEM
calculation consist of at least two program runs.
The preprocessing run generates sidechain rotamers,
calculates molecular mechanics energy terms, writes
a restart file for the postprocessing run and
sets up the necessary input for the continuum
electrostatics calculations that are done by
independent
MEAD
programs (my_xyz_solver).
This approach allows a very efficient and simple
parallelization without communication overhead
(see the PERL scripts rst-gcem-xyz.pl provided
with the examples for an example).
The postprocessing can also be run several times
alternating with recomputation of the electrostatic
energies to eliminate energetically unfavorable conformers.
Thereby, the continuum electrostatics calculations
are refined resulting in a decrease of the error introduced
by the inflation of the low dielectric molecule interior
by the high-energy conformers.

The intended purpose of the interior dielectric regions
eps1set, eps2set and eps3set
was for modeling a quantum chemically treated region,
a classically treated region and
a region of solvent filled cavities or pores inside the molecular structure that
are to be excluded from the membrane dielectric, but they might as well be used
for other purposes.
The interior regions are inaccessible to mobile ions.
The intended purpose of the region elycavset was to model
ion accessible cavities or depressions in the protein surface
reaching into the membrane boundaries (e.g. gorges leading to a channel entrance
and trans-membrane pores with large diameter).

A membrane can be modeled by three-layer dielectric slab
representing the polar, ion-accessible headgroup regions
(optional) and the apolar ion-inaccessible core region of the membrane.
The membrane is requested with membz.

molname can contain multiple sites which can occur
in multiple forms called instances.
Details can be found on a
→ separate page.

program options and their defaults

general options are described above

program specific options for the continuum electrostatics part:

option andexpected data type

default setting(if any)

meaning

-grid_detail int

2

This option influences the detail of the auto-generated grids for the
computation of the electrostatic potentials with the finite difference method.
A higher value will result in a larger and finer grids.
The default value of 2 will in most cases suffice for high quality results,
sometimes a value of 3 can be advantageous (especially for membrane proteins,
where the apolar environment leads to more far-reaching electrostatic interactions).
If user supplied grids exist, the setting affects the automatic adjustment of the
innermost grid's size to the size of each site.
A higher value will result in a larger spacing of the site to the grid boundary.

value

meaning

0

coarse, not advisable for production

1

economic, similar to the settings in the old multiflex examples, can result in significant discretization errors

2

high, somewhat finer and larger grids than for 1, solvation energies and interaction energies should be largely converged with respect to grid spacing

3

very high, even finer and larger innermost grids, ensures very high-quality solvation and interaction energies also for long-range interactions in membrane proteins

4

ultra high, extremely fine grids, should not result in significant numeric differences of the result relative to 3, needs very much memory and computation time

-epsext float

Floating point number that defines the dielectric constant of the
solvent region and the region specified by elycavset (if applicable).

-epshead float

Floating point number that defines the dielectric constant of the
membrane head group region (if applicable / membz is defined).

-epscore float

Floating point number that defines the dielectric constant of the
membrane head group region (if applicable / membz is defined).

-epsin1 float

Floating point number that defines the dielectric constant of the
region specified by eps1set (if applicable).

-epsin2 float

Floating point number that defines the dielectric constant of the
region specified by eps2set (if applicable).

-epsin3 float

Floating point number that defines the dielectric constant of the
region specified by eps3set (if applicable).

-eps_ff float

1

Floating point number that defines the dielectric constant of the
reference environment used in force field parametrization.

String giving the prefix of a .pqr file that define the
dielectric region 2 with the dielectric constant epsin2 (if applicable).

-eps3set string

String giving the prefix of a .pqr file that define the
dielectric region 3 with the dielectric constant epsin3 (if applicable).

-elycavset string

String giving the prefix of a .pqr file that defines the
an exterior region with the dielectric constant epsext
and ionic strength defined by ionicstr or
ionicstr1 and ionicstr2 according to the membrane side (if applicable).

-membz z_lower_core z_upper_core z_lower_head z_upper_head

Boundaries of a three-layer dielectric slab
representing the polar, ion-accessible headgroup regions
and the apolar ion-inaccessible core region of the membrane (if applicable).
The membrane is perpendicular to the z-direction, hence the
components of its normal vector are given by (0,0,1).
The headgroup region extends from
z=z_lower_head to z=z_upper_head excluding the core region.
The core region extends from
z=z_lower_core to z=z_upper_core excluding the core region.
The headgroup region can be omitted by setting
z_lower_core = z_lower_head and
z_upper_core = z_upper_head.
The ionic strength on the membrane sides can be set to equal values
by specifying ionicstr or set to distinct values
by specifying ionicstr1 and ionicstr2.
The ionic strength within the boundaries of the core region
(e.g., in a channel pore defined by membhole or elycavset)
is linearly interpolated between ionicstr1 and ionicstr2.

-ionicstr float

Ionic strength in the exterior region(s) (by default in mol/l).
In receptor environment overridden by ionicstr1, ionicstr2 if specified.
If ionicstr is not specified, it is set to the same value as ionicstr1
or ionicstr2 if those values were specified (in that order of precedence).
Also used as default ionic strength in exterior region for the reference environment
of the model compounds if not specified otherwise in the corresponding
sitename.est or sitename.qst files.

-ionicstr1 float

Ionic strength in the upper exterior region(s) (z > z_upper_core) (by default in mol/l).
Overrides ionicstr. If not specified, set to ionicstr if that value was specified.

-ionicstr2 float

Ionic strength in the lower exterior region(s) (z < z_lower_core) (by default in mol/l).
Overrides ionicstr. If not specified, set to ionicstr if that value was specified.

-membhole radius {cent_x cent_y}

Exclude a cylindrical hole from the membrane.
The radius of the hole is given by float_1.
The optional arguments cent_x and cent_y specify the x and y
coordinates of the cylinder center, respectively (default values 0.0 and 0.0).
This option works seldom satisfactory for real proteins with non-cylindrical shape.
The use of elycavset and/or one of the interior dielectric regions, for example
eps3set, is to be preferred.

-inside

1

This option defines which membrane side is "inside" in the electrophysiological sense,
where the membrane potential is measured at the inside relative to the outside.
A value of 0 states that the lower membrane side is inside, while a value of 1 states
that the upper side is inside, where lower and upper refer to lower and higher values
of the z coordinate, respectively.

-capacitance

This flag triggers the calculation of the capacitance of the
system.
The capacitance is the accumulated charge in Ampere seconds = Coulomb per membrane potential in Volt.
Hence, the capacitance is given by default in Farad (F = As/V).
The calculation of the capacitance requires tighter convergence criteria for the
electrostatic potentials, which increases the required computation time.
The capacitive energy is only needed if there are multiple global conformations.
The difference between the the capacitive energy terms
of different global conformations is negligibly small under normal conditions
of naturally occurring values of the membrane potential, but depends quadratically on membrane potential.
Caution: The capacitance and hence the capacitive energy depend on the system size
(mainly on the area of the covered membrane region).
Therefore, the same grid size must be used for the innermost grid in the calculation
of the capacitance for all global conformations (normally automatically taken care of).

program specific options for the molecular mechanics part:

option andexpected data type

default setting(if any)

meaning

-charmm_par string

Name of a CHARMM parameter file
containing the CHARMM force field parameters.
Examples are provided with the examples in the subdirectories of examples/gcem.
The file is not read if there are only non-flexible sites
(according to the definitions in molname.rot)
and skip_stat_mmg and skip_stat_mmint
are given on the command line.

-charmm_top string

Name of a CHARMM topology file
containing the residue topology definitions.
Double bonds must be specified as such for correct automatic identification of the
sidechain torsions.
Examples are provided with in the subdirectories of examples/gcem.
The file is not read if there are only non-flexible sites
(according to the definitions in molname.rot)
and skip_stat_mmg and skip_stat_mmint
are given on the command line.

-rotlib string

rotlib

Prefix of the name of a file
containing the
Squirell backbone dependent rotamer library
(tested with 2002 and 2010 versions).
First, it is attempted to read rotlibname.dat as binary boost archive.
If the binary file is not found, it is attempted to read
rotlibname.txt.
The file is not read if there are only non-flexible sites
(according to the definitions in molname.rot).

-skip_stat_mmint

Flag that triggers the omission of molecular mechanics contributions
to the intrinsic energies of each non-flexible site
(according to the specification of the site type in molname.rot).

-skip_stat_mmg

Flag that triggers the omission of molecular mechanics contributions
to each interaction energy that involves a non-flexible site
(according to the specification of the site type in molname.rot).

-empirical_energy

Flag that triggers the use of an empirical energy function for the
part of the intrinsic energy that involves atoms of the site itself
and nearby backbone residues that are thought to determine the rotamer
propensity in the backbone dependent rotamer library.
The empirical energy is defined as
\(-\beta^{-1} \ln\left[ p\right]\), where p is the propensity for the
corresponding sidechain rotamer in the rotamer library.

Flag that triggers writing of the rotamer library as binary archive
(using the Boost serialization library).

other program specific options:

option andexpected data type

default setting(if any)

meaning

-gcem_dir string

gcem_input

Name of the subdirectory to which gcem will write input files for
the continuum electrostatics solvers and the corresponding job script(s).
In addition, the directory will contain some deprecated files that are
no longer used but may be helpful in debugging while making changes to the
source code.

-mead_path string

Path to the root directory of the
MEAD
distribution used for the calls to the continuum electrostatics solvers
in the job script(s).
The option is useful if the solver jobs are run on a different system,
for example a remote computer cluster.

-filter_ctof float

1e99

Intrinsic energy cutoff (in kcal/mol) for eliminating conformers with excessively
high intrinsic energy, which makes them unlikely to be populated significantly in equilibrium.
A conformer (maybe a sidechain rotamer) consisting of N instances
(may be different binding forms of the sidechain) is eliminated if
the minimum intrinsic energy of any instance of the conformer is
larger than the minimum intrinsic energy of all instances of the site
plus filter_ctof.

-filter_ctof_gs float

1e99

Energy cutoff (in kcal/mol) for eliminating conformers with excessively high energy,
which makes them unlikely to be populated significantly in equilibrium.
The elimination criterion is almost identical to the Goldstein criterion of
dead end elimination, with the exception that the hypothetical minimum energy
must be larger by filter_ctof_gs rather than by 0.
A conformer (maybe a sidechain rotamer) consisting of N instances
(may be different binding forms of the sidechain) is eliminated if
all of its instances fulfill the modified Goldstein criterion.

-cap_max_nin int

Flag that triggers the reduction of the number of instances of any site
to a number ≤ cap_max_nin.
Keep no more than cap_max_nin instances of a site with lowest sum of intrinsic energy
and sum of minimum interaction energies with all other sites.
The method is equivalent to repeated application of the modified Goldstein criterion
to all conformers with successively decreasing filter_ctof_gs until
the number of instances is equal to or smaller than cap_max_nin.

-print_debug

Flag that triggers detailed output regarding the read and generated data
to stdout (independent on the blab level). The flag is mainly intended
for verification purposes while debugging new program features.

input files

molname.pqr
The receptor structure.
The file format is described above under general options and input files.

molname.ogm and molname.mgm (optional)
defines the electrostatics grids for the site in the receptor environment
and their the model compound in the reference environment (bulk solution), respectively.
The file format is described above under general options and input files.

molname.rot
defines the sites of molname ans their types
For each site, the file contains a line of the format:

site_label site_type

The site label is constructed from
sitename-chainid-resid,
where chainid and resid correspond to the data found in molname.pqr
(chainid set to 0 if absent in molname.pqr).
site_type is one of:

value

meaning

0

ignored as site

1

flexible and titratable (site_name.est required, see below)
user defined conformers are read from molname.pqr
flexible parts must be present for each conformer
with unique instance IDs differing from 0,
For amino acid residues, additional sidechain rotamers are
generated for amino acid residues (currently not implemented for Pro),
using dihedral angles from the Squirell backbone dependent rotamer library.
For each conformer, the different forms defined in
site_name.est are generated.

2

flexible and non-titratable (site_name.est not required)
user defined conformers are read from molname.pqr
flexible parts must be present for each conformer
with unique instance IDs differing from 0,
For amino acid residues, additional sidechain rotamers are
generated for amino acid residues (currently not implemented for Pro),
using dihedral angles from the Squirell backbone dependent rotamer library.

3

non-flexible and titratable (site_name.est required, see below)
For the conformer found in molname.pqr,
the different forms defined in site_name.est are generated

4

quantum mechanically treated site (QM site, site_label.qst required, see below)
Conformational flexibility / orientational polarization effects are thought to be
considered within the QM treatment.
The model energy should also contain any energy terms due to bound ligands,
and any other energy contribution that apart from the interactions with the rest of the protein.
GCEM
The model energy of a QM site is expected to be completely defined in
site_label.qst,
where site_label is given by sitename-chainid-resid.
The coordinates and atomic partial charges of the QM site
are taken from molname.pqr or separate .pqr files
as described below for site_label.qst.

molname.con
defines the connectivity of the sites to the membrane sides.
The file is only required for membrane proteins.
For each site, the file contains a line of the format:

site_label connectivity

The site label is constructed from
sitename-chainid-resid,
where chainid and resid correspond to
the data found in molname.pqr
(chainid set to 0 if absent in molname.pqr).
connectivity is one of:

model energy, (relative) chemical potentials of the forms
(see also the program manual of
GMCT)

EpsRef

relative dielectric constant of the reference environment
to which the Gmodel value refers, specification is optional,
set to the same value as epsext if not specified

IstrRef

specification is optional,
ionic strength of the reference environment
to which the Gmodel value refers

proton, electron ...

Any keyword or line patterns not matching one of the other entries
of this table names a ligand type.
The following numbers specify the numbers of the ligand type bound by each form.

center

ignored, formerly used to define the atom of the model compound
to be used as center of the finite difference grids,
GCEM uses
the geometric center of the model compound.

remaining lines

residue name, atom name, atomic partial charges of the atoms in each form

site_label.qst
A .qst file defines forms of a QM site.
Example:

label

FES1010

FES1111

Gmodel

-1.4232

0

EpsRef

1

1

IstrRef

0

0

proton

1

2

electron

1

2

keyword / line

meaning

label

labels for the forms

Gmodel

model energy, (relative) chemical potentials of the forms
(see also the program manual of
GMCT)

EpsRef

specification is optional,
relative dielectric constant of the reference environment
to which the Gmodel value refers

IstrRef

specification is optional,
ionic strength of the reference environment
to which the Gmodel value refers

proton, electron ...

Any keyword not found among among the preceding keywords names a ligand type.
The following numbers specify the numbers of the ligand type bound by each form.

Conformational flexibility / orientational polarization effects are thought to be
considered within the QM treatment.
The model energy should also contain any energy terms due to bound ligands,
and any other energy contribution that apart from the interactions with the rest of the protein.
GCEM
For each instance, a structure file named
site_label_instance_label.pqr is expected to exist.
The model energy of a QM site is expected to be completely defined in
extended_site_label.qst,
where extended_site_label is constructed from the site_label
as found in molname.rot by inserting the conformer id confid.
The extended site label is given by
sitename-confid-chainid-resid, where the confid.
Atomic coordinates and partial charges of atoms of the site found
in any structure of the QM site are ignored if also found in molname.pqr.
Capping atoms and groups are assumed not to be present in the structures —
the user has to remove them during structure preparation.
Bonded molecular mechanics terms involving link-atoms
are assumed to be the equal for all instances of the QM site, and thus neglected.
If you wish to include such energy contributions, you can add them to the model
energy.
There are unavoidable ambiguities and inconsistencies in the treatment of
the link regions.
Therefore, it is advisable to choose the extent of the QM site such that
any significant conformational flexibility occurs well within the QM site.
In this way, the above assumption (equal energy contributions for
all instances of the QM site due to the link region) is justified
just as for the other site types.
The program functionality for the QM sites was not yet tested much
and may thus change in the future.

CHARMM parameter and residue topology files
The file format is described on www.charmm.org.
The automatic determination of the amino acid sidechain torsion angles requires
that double bonds are specified as such in the residue topology file.

my_2diel_solver

This program computes the electrostatic potential
and the corresponding electrostatic energy terms
of a site in a two-dielectric environment.
Additional features enable the use of this program for
visualization purposes and as helper program for
GCEM.

show more detail

The program calculates the electrostatic interaction energy of
sitename with backgroundname and the
Born solvation energy of sitename.
The calculation of the Born solvation energy requires a second
calculation with one of the programs my_xyz_solver as
reference point and to cancel the grid artifacts.

All quantities that are discretized on the finite difference grid
can be written as OpenDX three-dimensional volume data, two-dimensional
curves in slice planes or one-dimensional curves along lines.
The generated data can be used for visualization and verification
purposes. One might, for example, want to verify the correct assignment
of dielectric regions in a convenient way using a molecular structure viewer
that can read OpenDX volume data, as for example
VMD or PyMol.
The data can also be written as two-dimensional curves for plotting with the
GMT
program suite or other plotting software.

program usage

Program settings are specified with command line options.
Input files are described below.
The program is called with:

program_dir/my_2diel_solver {options} sitenamebackgroundname

where sitename is the name of the site used as prefix for
the corresponding .pqr and .ogm files.
backgroundname is the name of the background
(structure containing all atoms not belonging to any site)
and used as prefix for the corresponding .pqr file.
Outside of the use with GCEM,
sitename can refer to any molecular structure
whose electrostatic potential, solvation energy or
interaction energy with a second molecule
backgroundname within
the two-dielectric environment are of interest.

The intended purpose of the interior dielectric regions
eps1set and eps2set
was for modeling a quantum chemically treated region and a
a classically treated region, but they might as well be used
for other purposes.
The interior regions are inaccessible to mobile ions.
The intended purpose of the region elycavset was to model
ion accessible cavities or depressions in the protein surface
reaching into the membrane boundaries (e.g. gorges leading to a channel entrance
and trans-membrane pores with large diameter).

Examples for the program usage can be found in the subdirectories
located in examples/gcem.
GCEM
will create the input files and a job script job.sh using the my_xyz_solvers
in the directory gcem_dir/gcem.

program options and their defaults

general options are described above

program specific options for the continuum electrostatics calculation:

option andexpected data type

default setting(if any)

meaning

-epsext float

Floating point number that defines the dielectric constant of the
solvent region outside the solvent inaccessible volume of
backgroundname.

-epsin float

Floating point number that defines the dielectric constant of
the interior of backgroundname.

input files

sitename.pqr,
backgroundname.pqr and
additional .pqr files used to define the dielectric regions
are structure files in .pqr format.
The file format is described above under general options and input files.

sitename.mgm
defines the electrostatics grids for the site in the receptor environment.
The file format is described above under general options and input files.

my_3diel_solver

This program computes the electrostatic potential of a site in a three-dielectric environment
and the corresponding electrostatic energy terms.
Additional features enable the use of this program for
visualization purposes and as helper program for
GCEM.

show more detail

The program calculates the electrostatic interaction energy of
sitename with backgroundname and the
Born solvation energy of sitename.
The calculation of the Born solvation energy requires a second
calculation with one of the programs my_xyz_solver as
reference point and to cancel the grid artifacts.
If the option -fpt is specified, the electrostatic
interaction energies with other sites can be calculated.

All quantities that are discretized on the finite difference grid
can be written as OpenDX three-dimensional volume data, two-dimensional
curves in slice planes or one-dimensional curves along lines.
The generated data can be used for visualization and verification
purposes. One might, for example, want to verify the correct assignment
of dielectric regions in a convenient way using a molecular structure viewer
that can read OpenDX volume data, as for example
VMD or PyMol.
The data can also be written as two-dimensional curves for plotting with the
GMT
program suite or other plotting software.

program usage

Program settings are specified with command line options.
Input files are described below.
The program is called with:

program_dir/my_3diel_solver {options} sitenamebackgroundname

where sitename is the name of the site used as prefix for
the corresponding .pqr and .ogm files.
backgroundname is the name of the background
(structure containing all atoms not belonging to any site)
and used as prefix for the corresponding .pqr file.
Outside of the use with GCEM,
sitename can refer to any molecular structure
whose electrostatic potential, solvation energy or
interaction energy with a second molecule
backgroundname within
the three-dielectric environment are of interest.

The intended purpose of the interior dielectric regions
eps1set and eps2set
was for modeling a quantum chemically treated region and a
a classically treated region, but they might as well be used
for other purposes.
The interior regions are inaccessible to mobile ions.
The intended purpose of the region elycavset was to model
ion accessible cavities or depressions in the protein surface
reaching into the membrane boundaries (e.g. gorges leading to a channel entrance
and trans-membrane pores with large diameter).

Examples for the program usage can be found in the subdirectories
located in examples/gcem.
GCEM
will create the input files and a job script job.sh using the my_xyz_solvers
in the directory gcem_dir/gcem.

program options and their defaults

general options are described above

program specific options for the continuum electrostatics calculation:

option andexpected data type

default setting(if any)

meaning

-fpt string

Prefix of a file fptname.fpt for the calculation of site-site interaction energies

-epsext float

Floating point number that defines the dielectric constant of the
solvent region and the region specified by elycavset (if applicable).

-epsin1 float

Floating point number that defines the dielectric constant of the
region specified by eps1set (if applicable).

-epsin2 float

Floating point number that defines the dielectric constant of the
region specified by eps2set (if applicable).

-epshomo float

Floating point number that defines the dielectric constant of an
additional homogeneous dielectric (optional).

-eps1set string

String giving the prefix of a .pqr file that define the
dielectric region 1 with the dielectric constant epsin1 (if applicable).

-eps2set string

String giving the prefix of a .pqr file that define the
dielectric region 2 with the dielectric constant epsin2 (if applicable).

-ionicstr float

Ionic strength in the exterior region(s) (by default in mol/l).

input files

sitename.pqr,
backgroundname.pqr and
additional .pqr files used to define the dielectric regions
are structure files in .pqr format.
The file format is described above under general options and input files.

sitename.ogm
defines the electrostatics grids for the site in the receptor environment.
The file format is described above under general options and input files.

fptname.fpt
is a file containing the atomic partial charges and their coordinates for each instance
of each site found in a receptor

site and instance denote the site and instance to whose
charge distribution the following atomic partial charge belongs.
coordinate consists of three floating point values that denote the x, y and z
coordinates of a point respectively.
Optionally, the coordinate values can be enclosed by parenthesis and separated by commas.
charge is the atomic partial charge at the preceding coordinate.
The line breaks can be substituted with any number of whitespace characters.

my_Ndiel_solver

This program computes the electrostatic potential of a site in an
environment with an arbitrary number of dielectric and electrolyte regions
and the corresponding electrostatic energy terms.
Additional features enable the use of this program for
visualization purposes and as helper program for
GCEM.

show more detail

The program calculates the electrostatic interaction energy of
sitename with backgroundname and the
Born solvation energy of sitename.
The calculation of the Born solvation energy requires a second
calculation with one of the programs my_xyz_solver as
reference point and to cancel the grid artifacts.
If the option -fpt is specified, the electrostatic
interaction energies with other sites can be calculated.

All quantities that are discretized on the finite difference grid
can be written as OpenDX three-dimensional volume data, two-dimensional
curves in slice planes or one-dimensional curves along lines.
The generated data can be used for visualization and verification
purposes. One might, for example, want to verify the correct assignment
of dielectric regions in a convenient way using a molecular structure viewer
that can read OpenDX volume data, as for example
VMD or PyMol.
The data can also be written as two-dimensional curves for plotting with the
GMT
program suite or other plotting software.

program usage

Program settings are specified with command line options.
Input files are described below.
The program is called with:

program_dir/my_Ndiel_solver {options} sitenamebackgroundname

where sitename is the name of the site used as prefix for
the corresponding .pqr and .ogm files.
backgroundname is the name of the background
(structure containing all atoms not belonging to any site)
and used as prefix for the corresponding .pqr file.
Outside of the use with GCEM,
sitename can refer to any molecular structure
whose electrostatic potential, solvation energy or
interaction energy with a second molecule
backgroundname within
the N-dielectric environment are of interest.

An example of the program usage can be found in the directory
examples/my_Ndiel_solver.

program options and their defaults

general options are described above

input files

sitename.pqr,
backgroundname.pqr and
additional .pqr files used to define the dielectric regions
are structure files in .pqr format.
The file format is described above under general options and input files.

sitename.ogm
defines the electrostatics grids for the site in the receptor environment.
The file format is described above under general options and input files.

backgroundname.diel
is a file defining the dielectric regions

format of a single line corresponding to a dielectric region:

pqrname eps solrad

pqrname is the prefix of a structure file that determines the
ion-inaccessible volume of the region.eps is the relative dielectric constant of the dielectric region.solrad is the solvent probe sphere radius used to define the
dielectric region.
The priority of the regions increases from the top to the bottom of the
list. Solvent inaccessible regions of previous entries are overridden by the
solvent inaccessible regions of following entries.

backgroundname.ely
is a file defining the electrolyte regions

format of a single line corresponding to an electrolyte region:

pqrname istr ionrad

pqrname is the prefix of a structure file that determines the
ion-inaccessible volume of the region.istr is the ionic strength in the ion-accessible volume
of the electrolyte region.ionrad is the ion radius (Stern layer radius) used to define
the electrolyte region.
The priority of the regions increases from the top to the bottom of the
list. Ion accessible regions of previous entries are overridden by the
ion-accessible regions of following entries.

fptname.fpt
is a file containing the atomic partial charges and their coordinates for each instance
of each site found in a receptor

site and instance denote the site and instance to whose
charge distribution the following atomic partial charge belongs.
coordinate consists of three floating point values that denote the x, y and z
coordinates of a point respectively.
Optionally, the coordinate values can be enclosed by parenthesis and separated by commas.
charge is the atomic partial charge at the preceding coordinate.
The line breaks can be substituted with any number of whitespace characters.

my_memb_solver

This program computes the electrostatic potential
and the corresponding electrostatic energy terms
of a site in a
environment that models a lipid membrane, the receptor and the
solvent phases above and below the membrane
with up to 6 dielectric and 2 electrolyte regions.
Additional features enable the use of this program for
visualization purposes and as helper program for
GCEM.

show more detail

The program calculates the electrostatic interaction energy of
sitename with backgroundname and the
Born solvation energy of sitename.
The calculation of the Born solvation energy requires a second
calculation with one of the programs my_xyz_solver as
reference point and to cancel the grid artifacts.
If the option -fpt is specified, the electrostatic
interaction energies with other sites can be calculated.

All quantities that are discretized on the finite difference grid
can be written as OpenDX three-dimensional volume data, two-dimensional
curves in slice planes or one-dimensional curves along lines.
The generated data can be used for visualization and verification
purposes. One might, for example, want to verify the correct assignment
of dielectric regions in a convenient way using a molecular structure viewer
that can read OpenDX volume data, as for example
VMD or PyMol.
The data can also be written as two-dimensional curves for plotting with the
GMT
program suite or other plotting software.

program usage

Program settings are specified with command line options.
Input files are described below.
The program is called with:

program_dir/my_memb_solver {options} sitenamebackgroundname

where sitename is the name of the site used as prefix for
the corresponding .pqr and .ogm files.
backgroundname is the name of the background
(structure containing all atoms not belonging to any site)
and used as prefix for the corresponding .pqr file.
Outside of the use with GCEM,
sitename can refer to any molecular structure
whose electrostatic potential, solvation energy or
interaction energy with a second molecule
backgroundname within
the membrane environment are of interest.

The intended purpose of the interior dielectric regions
eps1set, eps2set and eps3set
was for modeling a quantum chemically treated region,
a classically treated region and
a region of solvent filled cavities or pores inside the molecular structure that
are to be excluded from the membrane dielectric, but they might as well be used
for other purposes.
The interior regions are inaccessible to mobile ions.
The intended purpose of the region elycavset was to model
ion accessible cavities or depressions in the protein surface
reaching into the membrane boundaries (e.g. gorges leading to a channel entrance
and trans-membrane pores with large diameter).

A membrane can be modeled by three-layer dielectric slab
representing the polar, ion-accessible headgroup regions
(optional) and the apolar ion-inaccessible core region of the membrane.
The membrane is requested with membz.

Examples for the program usage can be found in the directory
examples/gcem/bR of the distribution.
GCEM
will create the input files and a job script job.sh using the my_xyz_solvers
in the directory gcem_dir/gcem.

program options and their defaults

general options are described above

program specific options for the continuum electrostatics calculation:

option andexpected data type

default setting(if any)

meaning

-fpt string

Prefix of a file fptname.fpt for the calculation of site-site interaction energies

-epsext float

Floating point number that defines the dielectric constant of the
solvent region and the region specified by elycavset (if applicable).

-epshead float

Floating point number that defines the dielectric constant of the
membrane head group region (if applicable / membz is defined).

-epscore float

Floating point number that defines the dielectric constant of the
membrane head group region (if applicable / membz is defined).

-epsin1 float

Floating point number that defines the dielectric constant of the
region specified by eps1set (if applicable).

-epsin2 float

Floating point number that defines the dielectric constant of the
region specified by eps2set (if applicable).

-epsin3 float

Floating point number that defines the dielectric constant of the
region specified by eps3set (if applicable).

-epshomo float

Floating point number that defines the dielectric constant of an
additional homogeneous dielectric (optional).

-eps1set string

String giving the prefix of a .pqr file that define the
dielectric region 1 with the dielectric constant epsin1 (if applicable).

-eps2set string

String giving the prefix of a .pqr file that define the
dielectric region 2 with the dielectric constant epsin2 (if applicable).

-eps3set string

String giving the prefix of a .pqr file that define the
dielectric region 3 with the dielectric constant epsin3 (if applicable).

-elycavset string

String giving the prefix of a .pqr file that defines the
an exterior region with the dielectric constant epsext
and ionic strength defined by ionicstr or
ionicstr1 and ionicstr2 according to the membrane side (if applicable).

-membz z_lower_core z_upper_core z_lower_head z_upper_head

Boundaries of a three-layer dielectric slab
representing the polar, ion-accessible headgroup regions
and the apolar ion-inaccessible core region of the membrane (if applicable).
The membrane is perpendicular to the z-direction, hence the
components of its normal vector are given by (0,0,1).
The headgroup region extends from
z=z_lower_head to z=z_upper_head excluding the core region.
The core region extends from
z=z_lower_core to z=z_upper_core excluding the core region.
The headgroup region can be omitted by setting
z_lower_core = z_lower_head and
z_upper_core = z_upper_head.
The ionic strength on the membrane sides can be set to equal values
by specifying ionicstr or set to distinct values
by specifying ionicstr1 and ionicstr2.
The ionic strength within the boundaries of the core region
(e.g., in a channel pore defined by membhole or elycavset)
is linearly interpolated between ionicstr1 and ionicstr2.

-ionicstr float

Ionic strength in the exterior region(s) (by default in mol/l).

-ionicstr1 float

Ionic strength in the upper exterior region(s) (z > z_upper_core) (by default in mol/l).
Overrides ionicstr. If specified, also ionicstr2 must be specified.

-ionicstr2 float

Ionic strength in the lower exterior region(s) (z < z_lower_core) (by default in mol/l).
Overrides ionicstr. If specified, also ionicstr1 must be specified.

-membhole radius {cent_x cent_y}

Exclude a cylindrical hole from the membrane.
The radius of the hole is given by float_1.
The optional arguments cent_x and cent_y specify the x and y
coordinates of the cylinder center, respectively (default values 0.0 and 0.0).
This option works seldom satisfactory for real proteins with non-cylindrical shape.
The use of elycavset and/or one of the interior dielectric regions, for example
eps3set, is to be preferred.

input files

sitename.pqr,
backgroundname.pqr and
additional .pqr files used to define the dielectric regions
are structure files in .pqr format.
The file format is described above under general options and input files.

sitename.ogm
defines the electrostatics grids for the site in the receptor environment.
The file format is described above under general options and input files.

fptname.fpt
is a file containing the atomic partial charges and their coordinates for each instance
of each site found in a receptor

site and instance denote the site and instance to whose
charge distribution the following atomic partial charge belongs.
coordinate consists of three floating point values that denote the x, y and z
coordinates of a point respectively.
Optionally, the coordinate values can be enclosed by parenthesis and separated by commas.
charge is the atomic partial charge at the preceding coordinate.
The line breaks can be substituted with any number of whitespace characters.

my_membpot_solver

This program computes the electrostatic trans-membrane potential
and the corresponding electrostatic energy terms in a
environment that models a lipid membrane, the protein and the
solvent phases above and below the membrane
with up to 6 dielectric and 2 electrolyte regions
(same as for my_memb_solver).
Additional features enable the use of this program for
visualization purposes and as helper program for
GCEM.

show more detail

The program calculates the electrostatic interaction energy of
backgroundname with the charge distribution causing the electrostatic trans-membrane
potential.
If the option -fpt is specified, the electrostatic
interaction energies of sites with the charge distribution causing the electrostatic trans-membrane
potential are calculated.
If the option -capacitance is specified, the program calculates
the capacitance of the receptor-membrane system.

All quantities that are discretized on the finite difference grid
can be written as OpenDX three-dimensional volume data, two-dimensional
curves in slice planes or one-dimensional curves along lines.
The generated data can be used for visualization and verification
purposes. One might, for example, want to verify the correct assignment
of dielectric regions in a convenient way using a molecular structure viewer
that can read OpenDX volume data, as for example
VMD or PyMol.
The data can also be written as two-dimensional curves for plotting with the
GMT
program suite or other plotting software.

The theoretical basis of computing the electrostatic trans-membrane potential
is described in
→ this paper (Roux, 1997).
A constant offset potential of \(\pm 0.5 \Psi\)
is added to all ion-accessible grid points on the inner and outer side of the membrane, respectively.
Equivalently, an effective uniform charge distribution can be assigned to
the ion accessible volume on either membrane side
(see → here).
For the finite difference representation, an effective charge of
\(q^{\mathrm{eff}} = \pm \frac{\overline{\kappa}^{2} h^{3} \Psi}{8 \pi}\),
is assigned to each ion accessible grid point
of the inner or outer membrane side, respectively.
Here, \(\overline{\kappa}^{2}\) is the inverse Debye length and
\(h\) is the grid spacing.
As it turned out, this alternative way of adding the offset potential
was also implemented by Roux and coworkers for the
PBEQ module of CHARMM.

program usage

Program settings are specified with command line options.
Input files are described below.
The last command line argument is the name of the molecule.
The program is called with:

program_dir/my_membpot_solver {options} backgroundname

backgroundname is the name of the background
(structure containing all atoms not belonging to any site)
and used as prefix for the corresponding .pqr and .ogm files.
Outside of the use with GCEM,
backgroundname can refer to any molecular structure
whose energy inside the membrane potential is of interest.

The intended purpose of the interior dielectric regions
eps1set, eps2set and eps3set
was for modeling a quantum chemically treated region,
a classically treated region and
a region of solvent filled cavities or pores inside the molecular structure that
are to be excluded from the membrane dielectric, but they might as well be used
for other purposes.
The interior regions are inaccessible to mobile ions.
The intended purpose of the region elycavset was to model
ion accessible cavities or depressions in the protein surface
reaching into the membrane boundaries (e.g. gorges leading to a channel entrance
and trans-membrane pores with large diameter).

A membrane can be modeled by three-layer dielectric slab
representing the polar, ion-accessible headgroup regions
(optional) and the apolar ion-inaccessible core region of the membrane.
The membrane is requested with membz.

Examples for the program usage can be found in the directories
examples/my_membpot_solver and examples/gcem/bR of the distribution.
GCEM
will create the input files and a job script job.sh using the my_membpot_solver
in the directory gcem_dir/membpot.

program options and their defaults

general options are described above

program specific options for the continuum electrostatics calculation:

option andexpected data type

default setting(if any)

meaning

-fpt string

Prefix of a file fptname.fpt for the calculation of site-site interaction energies

-epsext float

Floating point number that defines the dielectric constant of the
solvent region and the region specified by elycavset (if applicable).

-epshead float

Floating point number that defines the dielectric constant of the
membrane head group region (if applicable / membz is defined).

-epscore float

Floating point number that defines the dielectric constant of the
membrane head group region (if applicable / membz is defined).

-epsin1 float

Floating point number that defines the dielectric constant of the
region specified by eps1set (if applicable).

-epsin2 float

Floating point number that defines the dielectric constant of the
region specified by eps2set (if applicable).

-epsin3 float

Floating point number that defines the dielectric constant of the
region specified by eps3set (if applicable).

-epshomo float

Floating point number that defines the dielectric constant of an
additional homogeneous dielectric (optional).

-eps1set string

String giving the prefix of a .pqr file that define the
dielectric region 1 with the dielectric constant epsin1 (if applicable).

-eps2set string

String giving the prefix of a .pqr file that define the
dielectric region 2 with the dielectric constant epsin2 (if applicable).

-eps3set string

String giving the prefix of a .pqr file that define the
dielectric region 3 with the dielectric constant epsin3 (if applicable).

-elycavset string

String giving the prefix of a .pqr file that defines the
an exterior region with the dielectric constant epsext
and ionic strength defined by ionicstr or
ionicstr1 and ionicstr2 according to the membrane side (if applicable).

-membz z_lower_core z_upper_core z_lower_head z_upper_head

Boundaries of a three-layer dielectric slab
representing the polar, ion-accessible headgroup regions
and the apolar ion-inaccessible core region of the membrane (if applicable).
The membrane is perpendicular to the z-direction, hence the
components of its normal vector are given by (0,0,1).
The headgroup region extends from
z=z_lower_head to z=z_upper_head excluding the core region.
The core region extends from
z=z_lower_core to z=z_upper_core excluding the core region.
The headgroup region can be omitted by setting
z_lower_core = z_lower_head and
z_upper_core = z_upper_head.
The ionic strength on the membrane sides can be set to equal values
by specifying ionicstr or set to distinct values
by specifying ionicstr1 and ionicstr2.
The ionic strength within the boundaries of the core region
(e.g., in a channel pore defined by membhole or elycavset)
is linearly interpolated between ionicstr1 and ionicstr2.

-ionicstr float

Ionic strength in the exterior region(s) (by default in mol/l).
Overridden by ionicstr1, ionicstr2 if specified.

-ionicstr1 float

Ionic strength in the upper exterior region(s) (z > z_upper_core) (by default in mol/l).
Overrides ionicstr. If specified, also ionicstr2 must be specified.

-ionicstr2 float

Ionic strength in the lower exterior region(s) (z < z_lower_core) (by default in mol/l).
Overrides ionicstr. If specified, also ionicstr1 must be specified.

-membhole radius {cent_x cent_y}

Exclude a cylindrical hole from the membrane.
The radius of the hole is given by float_1.
The optional arguments cent_x and cent_y specify the x and y
coordinates of the cylinder center, respectively (default values 0.0 and 0.0).
This option works seldom satisfactory for real proteins with non-cylindrical shape.
The use of elycavset and/or one of the interior dielectric regions, for example
eps3set, is to be preferred.

-inside

1

This option defines which membrane side is "inside" in the electrophysiological sense,
where the membrane potential is measured at the inside relative to the outside.
A value of 0 states that the lower membrane side is inside, while a value of 1 states
that the upper side is inside.

-capacitance

This flag triggers the calculation of the capacitance of the
system.
The capacitance is the accumulated charge in Ampere seconds = Coulomb per membrane potential in Volt.
Hence, the capacitance is given by default in Farad (F = As/V).

input files

backgroundname.pqr and
additional .pqr files used to define the dielectric regions
are structure files in .pqr format.
The file format is described above under general options and input files.

backgroundname.ogm
defines the electrostatics grids for the site in the receptor environment.
The file format is described above under general options and input files.

fptname.fpt
is a file containing the atomic partial charges and their coordinates for each instance
of each site found in a receptor

site and instance denote the site and instance to whose
charge distribution the following atomic partial charge belongs.
coordinate consists of three floating point values that denote the x, y and z
coordinates of a point respectively.
Optionally, the coordinate values can be enclosed by parenthesis and separated by commas.
charge is the atomic partial charge at the preceding coordinate.
The line breaks can be substituted with any number of whitespace characters.

get-curves-gmct

This program analyzes the output of
→
GMCT
and extracts information as for example
net binding probabilities of sites summed over all instances,
and \(\mathrm{p}K_{1/2}\) values.
The MEAD
library is only used for reading and writing of instance structures from the
GCEM
input.

show more detail

The program reads a
GMCT
output file, generates derived data and writes the results
to the directory curve_dir.

get-curves-gmct writes N-dimensional curves that describe the dependence
of a quantity (e.g., a binding probability) on the chemical potentials
of the ligands and the membrane potential.
One- two- and three-dimensional curves can, for example, be used for plotting with the
GMT
program suite or other plotting software.
The curves are written to subdirectories of curve_dir whose names
correspond to the respective keywords in molname.gc.setup
(example: the output of the midpoints of binding titrations is triggered with the
input option ibindhalf and the corresponding curves are written to curve_dir/bind_half).

The site and background structures in .pqr format needed to
construct assembled structures of molname with all
sites in their most highly populated instances are read from
gcem_dir (only needed if the input option itraj is set to 1).
Molecular structures in .pqr format can be visualized with molecular structure viewers
as for example
VMD or PyMol.

program usage

Program settings are specified within a setup file.
Input files are described below.
The last command line argument is the name of the molecule.
The program is called with:

program_dir/get-curves-gmct {options} molname

where molname is the name of the site used as prefix for
the corresponding setup file molname.gc.setup.

Examples for the program usage can be found in the directories
examples/DTPA_titra and examples/HEWL_titra of the
GMCT
distribution.

Name of the subdirectory that will contain the curves and structure files created by get-curves-gmct
(created automatically if not existent).

gcem_dir string

..

Name of the directory containing the gcem input.
The .pqr files are expected to be found in subdirectories with the name gcem_input or qmpb_input
If there are multiple conformations, the program expects to find the corresponding .pqr files
in the subdirectories (gcem|qmpb)_dir/confname.

blab int

0

Program verbosity as described under general program options.

nformat int

Keyword that triggers the reading of a block of format descriptions
for the chemical potentials (same order as in gmct_outfile.
The number gives the number of format lines in the block.
Each line consists of the keyword format followed by a conversion factor
and the number of decimal digits after the comma to which the converted
chemical potential, pmf or membrane potential should be rounded.
If there is a proton chemical potential and an electron chemical potential
which should appear in the form of a pH value and a reduction potential in V in the
program output,
the block could look like this:

nformat

2

format

-0.733004200213

3

format

-0.043364101728026

4

ipmf int

0

Flag whether or not (1 = yes, 0 = no)
to output of a proton motive force
(pmf \(= \mu_{\mathrm{H}^{+}}^{\mathrm{out}}
- \mu_{\mathrm{H}^{+}}^{\mathrm{in}}
- \mathrm{F}\Delta \psi =
\bar{\mu}_{\mathrm{H}^{+}}^{\mathrm{out}}
- \bar{\mu}_{\mathrm{H}^{+}}^{\mathrm{in}}
\))
in the curves instead separate
proton chemical potentials and an electric membrane potential.
If ipmf is set to 1, the membrane potential is substituted with the pmf and the
proton chemical potential on the outside is omitted.

Flag whether or not (1 = yes, 0 = no)
to output occupation probability curves for each possible number of
ligands bound to a site
(permutation of the possible numbers of bound ligands of each type).

ifprob int

0

Flag whether or not (1 = yes, 0 = no)
to output occupation probability curves for each (binding) form of a site.

icprob int

0

Flag whether or not (1 = yes, 0 = no)
to output occupation probability curves for each global conformation and
each conformer of a site.
Conformers are expected be ordered sequentially in gmct_outfile
(automatically done if the input for
GMCT was generated with
GCEM)
, where each conformer possesses the same (binding) forms.

ibmean int

0

Flag whether or not (1 = yes, 0 = no)
to output the average number of bound ligands of each type bound by a site.

imaxp int

0

Flag whether or not (1 = yes, 0 = no)
to output the most highly populated instance of a site.

ibmaxp int

0

Flag whether or not (1 = yes, 0 = no)
to output the most highly populated number of bound ligands of each ligand type for a site.

ifmaxp int

0

Flag whether or not (1 = yes, 0 = no)
to output the most highly populated (binding) form of a site.

icmaxp int

0

Flag whether or not (1 = yes, 0 = no)
to output the most highly populated global conformation and
the most highly populated conformer of a site.

ihalf int

0

Flag whether or not (1 = yes, 0 = no)
to output the midpoints for each transition between the possible pairs
of instances of a site.
One chemical potential is varied, while the others have fixed values.
The midpoint is the value of the varied chemical potential
at which the population of the two instances is equal
at fixed values of all other chemical potentials.

ibhalf int

0

Flag whether or not (1 = yes, 0 = no)
to output the midpoints for each transition between consecutive
binding numbers of each ligand type and each site.
The midpoint is the chemical potential of the corresponding ligand type
at which the population of the consecutive binding numbers is equal
at fixed values of all other chemical potentials, the pmfs and the membrane potential.
An example application of this option is the determination of
\(\mathrm{p}K_{1/2}\) values, which in general depend on the
chemical potentials of all other ligands present in the system
and a possibly present membrane potential.

ifhalf int

0

Flag whether or not (1 = yes, 0 = no)
to output the midpoints for each transition between the possible pairs
of (binding) forms of a site.
One chemical potential is varied, while the others have fixed values.
The midpoint is the value of the varied chemical potential
at which the population of the two (binding) forms is equal
at fixed values of all other chemical potentials.

ichalf int

0

Flag whether or not (1 = yes, 0 = no)
to output the midpoints for each transition between the possible pairs
of global conformations and of conformers of a site.
One chemical potential is varied, while the others have fixed values.
The midpoint is the value of the varied chemical potential
at which the population of the two global conformations or site conformers is equal
at fixed values of all other chemical potentials.

iginst int

0

Flag whether or not (1 = yes, 0 = no)
to output the free energy of the instances of a site
calculated from \(G = -\beta^{-1} \ln\left[p\right]\),
where \(p\) is the occupation probability of the instance.

igbind int

0

Flag whether or not (1 = yes, 0 = no)
to output the free energy of each possible number of
ligands bound to a site
(permutation of the possible numbers of bound ligands of each type)
calculated from \(G = -\beta^{-1} \ln\left[p\right]\),
where \(p\) is the occupation probability of the binding number.

igconf int

0

Flag whether or not (1 = yes, 0 = no)
to output the free energy of each global conformation and each site conformer
calculated from \(G = -\beta^{-1} \ln\left[p\right]\),
where \(p\) is the occupation probability of the global conformation or site conformer.

itraj int

0

Flag whether or not (1 = yes, 0 = no)
to write assembled structures of molname with all
sites in their most highly populated instances for each permutation of
chemical potential values.

input files

molname.gc.setup
contains the program options as described above.

molname.gmct.out (or user specified name)
contains the calculation results written by
GMCT to stdout.
The results comprise the
populations of the instances and the global conformations
for each permutation of chemical potential values as computed from a
Monte Carlo simulation with the program gmct or an
analytical calculation with smt.
The format is described in the user manual of the
program suite
GMCT.

.pqr files for the instances of all sites and the background.
(only needed if the input option itraj is set to 1)
The file format is described above under general options and input files.

pqr2SolvAccVol

This program computes an analytical representation of the
solvent (in)accessible volume of a molecular structure and
writes it to a binary or ASCII file which can be read by
the programs my_xyz_solver to avoid time-consuming
recomputation of the volume.

show more detail

Program settings are specified with command line options.
Input files are described below.
The program uses a .pqr file as input (file format described above).
The command line option -solrad specifies the solvent probe sphere radius.
(described under general options above).
The last command line argument is the name of the molecule.
The program is called with:

program_dir/pqr2SolvAccVol {options} molname

where molname is the name of the site used as prefix for
the corresponding .pqr file.

pqr2IonAccVol

This program computes an analytical representation of the
ion inaccessible volume (solvent inaccessible volume with
Stern layer radius as probe radius and extended by the Stern layer radius)
of a molecular structure and writes it to a binary or ASCII file which
can be read by the programs my_xyz_solver to avoid time-consuming
recomputation of the volume.

show more detail

Program settings are specified with command line options.
Input files are described below.
The program uses a .pqr file as input (file format described above).
The program uses a .pqr file as input (file format described above).
The command line option -sterln specifies the ion probe sphere radius
that defines the thickness of the Stern layer.
The command line option -gmt specifies the dimensions of
the curve, the projection plane and optionally a second slice plane
(described under general options above).
The last command line argument is the name of the molecule.
The program is called with:

program_dir/pqr2IonAccVol {options} molname

where molname is the name of the site used as prefix for
the corresponding .pqr file.

pqr2crv

This program computes a projection of the solvent inaccessible volume of a
molecular structure or a slice thereof into a user defined plane for visualization purposes.
An example is shown in the above figure showing the trans-membrane potential distribution across Amt-1.

show more detail

program usage

Program settings are specified with command line options.
Input files are described below.
The program uses a .pqr file as input (file format described above).
The command line option -solrad specifies the solvent probe sphere radius.
The command line option -gmt specifies the dimensions of
the curve, the projection plane and optionally a second slice plane
(described under general options above).
The last command line argument is the name of the molecule.
The program is called with:

program_dir/pqr2crv {options} molname

where molname is the name of the site used as prefix for
the corresponding .pqr file.

If the flag -blab3 is set, .pqr file(s) for visualization
of the projection and slice planes will be written.
These files can be used to check the correctness of the
settings for the projection and slice planes in a molecule viewer
like VMD or
PyMol.

potential

This program calculates the electrostatic potential due to the charge distribution of molname

show more detail

Potential calculates the potential due to the molecule, molname
whose coordinates radii, and charges are specified in a file molname.pqr,
and writes to standard output its
values at points specified by the molname.fpt file. The units of the
output potentials are input charge units divided by input length
units. For typical PDB derived input files, this would be
elementary charges/Ångström. In that case,
to get (kcal/mol)/elementary charge, multiply by 332.063202.

The CoarseFieldOut and CoarsefieldInit options will cause the
coarsest potential lattice to be written to, or initialized from, an
AVS (Advanced Visualization System) "field" file. By default, the
units are as above for the the output potentials, but if you give the
option " AvsScaleFactor f", where f is a floating point number, the
field will be scaled by that factor.

The program is called with:

program_dir/potential {options} {molname}

Potential requires the files molname.pqr and molname.ogm. molname.fpt
is not strictly required, but if neither the .fpt file nor the
CoarseFieldOut option are used, the program produces no output of the
calculated potentials. The epsin option is mandatory.

prefix of input file name.fld.
This option triggers reading of initial values for the electrostatic potential
for the coarsest (first) grid level in AVS "field" format from name.fld.

-AvsScaleFactor float

Scale factor for unit conversion of the the electrostatic potential
for the coarsest (first) grid level in AVS "field" format.

solvate

This program calculates the electrostatic solvation energy of a molecule
in a solvent.

show more detail

Solvate calculates the Born solvation energy of a molecule — that is,
the difference in the electrostatic work required to bring its atom
charges from zero to their full values in solvent versus vacuum.

The program is called with:

program_dir/solvate {options} {molname}

molname is the molecule for
which the the calculation is to be done and whose coordinates radii
and charges are specified in a file molname.pqr. Solvate requires a
molname.pqr file and a molname.ogm file (see below) as inputs.
The -epsin option is mandatory.

The Born solvation energy is written to standard output in kcal/mol.
Physical conditions and units for I/O can be set by flags on the
command line (see general options above). By default solvate assumes we are going
from vacuum (eps=1) to water (eps=80). Note that solvate uses the
epssol and epsvac flags rather than the epsext options to control
solvent conditions. You can try it out on a sphere to check agreement
with the Born formula. See the example in the directory example/solvate/born.

An example for the program usage can be found in the directory
examples/solvate of the distribution.

general options are described above

program specific options:

option andexpected data type

default setting(if any)

meaning

-ReactionField

flag that triggers the output of the electrostatic reaction potential
(difference between the electrostatic potential in protein + solvent environment vs. vacuum)
at the coordinates specified in molname.fpt (solvate and solinprot only).
The corresponding potential values will be written to molname.rf.

solinprot

Solinprot calculates the electrostatic part of the transfer energy for bringing a compound from the
bulk solvent to a molecular environment.

show more detail

Calculates the "solvation" energy of the molecule named by solute
(for which there must be solute.pqr and solute.ogm files) inside the molecule
named by protein (for which the must be a protein.pqr file) which is,
in turn, in some solvent (water, by default). The interior of the
solute is presumed to have the dielectric constant, epsin1 (given by
the epsin1 flag) and regions interior to the protein but exterior to
the solute are presumed to have the dielectric constant, epsin2 (given
by the epsin2 flag) and regions exterior to both protein and solute
are presumed to have dielectric constant, epsext (80.0 by default).
The protein may contain charges and their interaction with the solute
will contribute to "solvation energy." The solvated energy is
calculated relative to a vacuum calculation in which the dielectric
constant has a value of epsin1 inside the solute and epsvac (1.0 by
default) outside.
The calculation works like this First the potential due to the solute
charges (call them rho_solute) in the above described dielectric
environment is calculated. Call this potential phi_sol. Next the
potential due to rho_solute is calculated in the vacuum dielectric
environment described above. Call this potential phi_vac. The
reaction field component of the solvation energy is then,
(rho_solute*phi_sol - rho_solute*phi_vac)/2, where "*" indicates a
suitable sum or integral of charge times potential. So far, this is
the same as the solvate program except for the three dielectric
environment on the solvated side. We also need the contribution due
to protein charges, rho_protein, interacting with the solute
rho_protein*phi_sol.

The program is called with:

program_dir/solinprot {options} solute protein

The program expects to find the file solute.pqr
containing the solute's coordinates charges and radii
and the file solute.ogm containg the grid definitions.
In addition, the file protein.pqr is needed
which contains the coordinates, charges and radii of the
protein
The options are similar to those for the solvate program except that the
epsin1 and epsin2 options are mandatory and the
epsin option is forbidden. The epsext option is used instead of epssol for
specifying the solvent dielectric.
The ProteinField flag is also available.

general options are described above
program specific options:

option andexpected data type

default setting(if any)

meaning

-ReactionField

flag that triggers the output of the electrostatic reaction potential
(difference between the electrostatic potential in protein + solvent environment vs. vacuum)
at the coordinates specified in molname.fpt.
The corresponding potential values will be written to molname.rf.

-ProteinField

flag that triggers the output of the electrostatic potential due to the
protein charge distribution at the coordinates specified in
molname.fpt.
The corresponding potential values will be written to molname.pf.

multiflex

Multiflex is a program for the automated preparation of the necessary input
for continuum electrostatic calculations on binding equilibria
within a traditional two-state model.

show more detail

Multiflex does the electrostatic part of a titration calculation for a
multi site titrating molecule. It can do single conformer calculations
based on the methods described in Karplus and Bashford (1990) and
Bashford and Gerwert (1992)
which assumes a rigid molecule, or
it can included limited conformational flexibility by the method of
You and Bashford (1995).
In the latter case, the user must
supply the coordinates for the conformational variants and the
corresponding non electrostatic energies. This can be done with
a program like CHARMM.

See Ullmann1999 for a review of
titration calculations with continuum electrostatics within a classical two-state model.
→
GCEM
(see program description above)
provides an alternative to multiflex with additional features
and a more detailed receptor model.

For single conformer calculation, multiflex works like the old
multimead program. It takes molname.pqr, molname.ogm, molname.mgm,
molname.sites and site_type.st files as inputs (see below) and as its
main outputs, produces a molname.pkint file, which contains the
calculated intrinsic pKas, a molname.g file, which contains site site
interactions in units of charge squared per length and a molname.summ
file which summarizes the self and background contributions to the
intrinsic pK of each site.
The molname.pkint file and the molname.g
file can be used directly as input to redti, the program for
calculating titration curves.
Alternatively, these files can be used with Monte Carlo simulation
programs
(e.g., XMCTI)
to compute titration curves for systems, in which the
computational cost of redti is prohibitive.
Multiflex also produces a file,
molname.sitename.summ for each titrating site which contains some
summary information that is mainly interesting for multi conformer
(flexible) calculations and it produces a large number of .potat
files, which are binary files that are useful if a job is interrupted
and restarted (see below). The epsin flag is mandatory. Other flags
can be used to change units and/or physical conditions or include a
membrane.

For the "flexible" calculations which involve multiple conformers,
multiflex needs the input files described above but it also needs
additional files For each flexible site there must be a file,
molname.sitename.confs (see below).
These .confs files tell about the possible conformers of that site
and their non-electrostatic energies (see below under input files).
For each conformer named in a .confs file, there must be a .pqr file
having the coordinates, charges and radii for the whole protein with
the site sitename in conformer confname
(see below under input files).
You might need a lot of of .pqr files, for example, the
You and Bashford calculations on lysozyme had 36 conformers for each of 12
sites and one molname.pqr file is always needed so 433 .pqr files were
needed for each titration calculation.
Flexible and single conformer sites can co exist in the same molecule.
The way multiflex tells the difference is that flexible sites have
confs files and single conformer files don't.

Program settings are specified with command line options.
Special input files are described below.
Example applications are found in the directory examples/multiflex
of the
MEAD
distribution.
The last command line arguments specify the name of the molecule.
The program is called with:

program_dir/multiflex {options} molname

where molname is the name of the molecule for which the calculation
will be done.
molname is used as prefix for some input files.

Examples for the program usage can be found in the directory
examples/multiflex of the distribution.

Set up a membrane parallel to the x-y plane (membrane normal in z direction).
The membrane is modeled as low-dielectric slab whose lower and upper boundaries
are given by z_lower and z_upper, respectively.

-membhole radius {cent_x cent_y}

Exclude a cylindrical hole from the membrane.
The radius of the hole is given by float_1.
The optional arguments cent_x and cent_y specify the x and y coordinates of the cylinder
center, respectively (default values 0.0 and 0.0).

-site int

Integral number N whose specification causes multiflex
to do the electrostatics calculations only for the Nth site specified in the .sites file.

input files

molname.pqr
is the structure of molname.
The file format is described above under general options and input files.

molname.sites
contains a list of the titratable sites of molname.
For each site, the file contains a line of the following format:

res_num site_type {chain_id}

res_num is the residue number of the site which will be matched to the
residue numbers of the atoms in molname.pqr.
site_type is the type of the site (usually similar to the residue name).
A file named site_type.st is needed for each site_type (see below).
The optional argument chain_id is the chain id of the site.
chain_id is matched to the corresponding chain_id in
molname.pqr and must consequently also be specified in molname.pqr also if specified
here.

site_type.st
multiflex expects a file with the name site_type.st to specify
details concerning each site_type that appears in the
molname.sites file.
The first line contains one floating point number which
is the model compound pK for that type of site.
All remaining lines are of the format:

resname atomname prot_charge deprot_charge

where resname and atomname, along with the
res_nums given in the molname.sites file match an
atom in the molname.pqr file that is part of
a titrating site.
prot_charge is the charge of this atom in
the protonated state and deprot_charge is its charge in the
deprotonated state.
It is expected that the sum of all
prot_charges subtracted from the sum of all deprot_charges
will equal one.
An extension of this file's syntax is planned to allow a regular
expression to be given that will specify how a model compound
is to be made from the atom coordinates in molname.pqr.

molname.sitename.conf
multiflex expects a file with the name molname.sitename.conf
to exist for each flexible site (otherwise the site is treated as non-flexible).
Here, "sitename" is constructed from the molname.sites file by the procedure,

"7 GLU" → "GLU-7"

These .confs files tell about the possible conformers of that site.
They contain lines of the form

confname mac_non_elstat_energy mod_non_elstat_energy

where confname is the name of the conformer and the next two entries
are non electrostatic energies of the conformer in the macromolecule
and in the model compound corresponding to this site. The energies
must be given in kcal/mol unless the econv flag (see below)
has been given to change the energy units.

molname.sitename.confname.pqr
For each conformer named in a .confs file, there must be a .pqr file
having the coordinates, charges and radii for the whole protein with
the site sitename in conformer confname
The file format is described above under general options and input files.

molname.ogm and molname.mgm
defines the electrostatics grids for the site in the receptor environment
and for the model compound in the reference environment (bulk solution), respectively.
The file format is described above under general options and input files.

output files

molname.pkint
contains the calculated intrinsic pKa values of the sites in pK units.
Each site has no more than two protonation forms (protonated and deprotonated).
The protonation forms are called instances in the nomenclature of
GCEM.
The intrinsic pKa \(\mathrm{p}K_{\mathrm{a},i}^{\mathrm{int}}\) of a site \(i\)
is equal to the relative intrinsic energy of the two protonation forms in
pK units
\[
\mathrm{p}K_{\mathrm{a},i}^{\mathrm{int}} = \frac{1}{\mathrm{RT}\ln 10}
\left(\GintNoConf{i}{\mathrm{deprot}} - \GintNoConf{i}{\mathrm{prot}}\right)
\]

molname.g
contains the (relative) interaction energies between
each pair of sites in units of
\( \mathrm{e}_{\circ}^{2} / \unicode[Times,serif]{xC5} \)
(charge squared per Ångström).
The relative interaction energy between the sites \(i\) and \(j\) in terms of
instance-instance interaction energies as used by
GCEM.
is given by
\[
g_{i,j} = \GinterNoConf{i}{ \mathrm{prot}}{j}{ \mathrm{prot}} - \GinterNoConf{i}{\mathrm{deprot}}{j}{\mathrm{deprot}}
- \left(
\GinterNoConf{i}{ \mathrm{prot}}{j}{\mathrm{deprot}} - \GinterNoConf{i}{\mathrm{deprot}}{j}{\mathrm{deprot}}
+\GinterNoConf{i}{\mathrm{deprot}}{j}{ \mathrm{prot}} - \GinterNoConf{i}{\mathrm{deprot}}{j}{\mathrm{deprot}}
\right)
\]
which becomes after collection of terms
\[
g_{i,j} = \GinterNoConf{i}{ \mathrm{prot}}{j}{ \mathrm{prot}}
+\GinterNoConf{i}{\mathrm{deprot}}{j}{\mathrm{deprot}}
-\GinterNoConf{i}{ \mathrm{prot}}{j}{\mathrm{deprot}}
-\GinterNoConf{i}{\mathrm{deprot}}{j}{ \mathrm{prot}}
\]
For each site, the file contains a line of the following format:

\(i\) \(j\) \(g_{i,j}\)

where \(i\) and \(j\) are the site ids in the order of the
molname.sites file and \(g_{i,j}\) is the interaction energy between the sites

molname.summ
contains a summary of the background and self energy contributions
to the intrinsic pKa values which summarizes the self (Born solvation)
and background contributions to the
intrinsic pKa value of each site.
See Ullmann1999 for a review of
titration calculations with continuum electrostatics within a classical two-state model.

molname.sitename.summ
contains a summary of the energy contributions
to the intrinsic pKa values file which summarizes the energy contributions to the
intrinsic pKa value of each conformer confname of the flexible site sitename.

redti

This program calculates titration curves with an approximate analytical method
(reduced site method).

show more detail

Redti solves the multiple site titration curve problem given a set of
intrinsic pKas (molname.pkint) and a site-site interaction matrix
(molname.g) using the reduced site method described by
Bashford & Karplus (1991) J. Phys. Chem. vol. 95, pp. 9556-61.
The input files can be obtained from multiflex (see above for a description of the files).
Redti is written in C rather than C++.
Its command line syntax is:

program_dir/redti {options} molname

where molname is the name of the molecule for which the calculation
will be done.
molname is used as prefix for some input files.

program specific options:

-cutoff float

20.5

Cutoff for the reduced site method.

-pHrange pH_min pH_max

-20 30.1

pH range for the titration calculation

-dry

Flag that causes redti to do a "dry run" in which it prints the number of sites
to be included in the reduced site calculation at each pH point. This
is useful for checking whether a calculation will require a prohibitive
amount of CPU time since CPU time will go exponentially in the reduced
site number.

input files

molname.pkint
contains the calculated intrinsic pKa values of the sites in pK units.

The MEAD distribution was updated to version 2.3.5.
This version fixes compilation issues with newer compilers
(GCC >= 7, Clang 8, ICC 2018).
Version 2.3.4 fixed compilation issues of the tool get-curves-gmct
with newer versions of the Boost library.
Changes in version 2.3.3 were restricted to a bugfix for the application
GCEM.
The bug occurred only in calculations with QM sites but without explicit
ligand atoms.
The bug resulted in a job script with calls to the electrostatics solvers
for the receptor environment that missed the QM-site related options
-epsin1 and -eps1set.
Version 2.3.2 fixed a critical bug that caused
rare segfaults in calculations with very big and/or fine grids
(see the file ChangeLog for details).
The bug was inherited from the original MEAD distribution
which is thus also affected.
All users are encouraged to update.

Feedback and bug reports

Your opinion and hints are very welcome.
Please provide detailed information about the program input and output when reporting a bug.
Often, running a program with a higher verbosity level (command line option -blab3)
helps to clarify the source of an error.