Title

Home

Research

Publications

CV

Software

GMCT

GCEM

ext. MEAD


Purpose of MEAD and our extensions

This page describes our extensions to the molecular electrostatics software suite MEAD. The original version of MEAD was written by Donald Bashford and used to be found → here. Unfortunately, Prof. Bashford retired and the original version is not available anymore. All functionality is however still included in extended MEAD.

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.

membrane potential distribution across the ammonium transporter Amt-1
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.

program options and their defaults

general options:

option and
expected data type
default setting
(if any)
meaning
-epsin float
-epsin1 float
-epsin2 float
-epsin3 float
dielectric constant of a molecular interior region
-epsext float80dielectric constant of an exterior / solvent region (potential, my_xyz_solver)
-epssol float80dielectric constant of an exterior / solvent region (multiflex and solvate)
-epsvac float1.0dielectric constant of vacuum (solvate)
-T float298.15absolute temperature in K
-solrad float1.4solvent probe radius for determining a solvent accessible surface / the solvent inaccessible volume
-sterln float2.0thickness of an ion-exclusion layer (Stern layer) for determining the ion-inaccessible volume
-ionicstr float0.0ionic 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 and
expected 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:

    ATOM/HETATM atnum atname resname resnum x y z charge radius {confid chainid siteid instid}
    ...


    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_1 grid_spacing_1 grid_dimension_1
    ...
    grid_center_N grid_spacing_N grid_dimension_N


    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)

    • traditional format (for command line option -ProteinField/-ReactionField):

      coordinate_1
      ...
      coordinate_N


      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.
    • extended format II (for command line option -fpt):

      site_of_1 instance_of_1 coordinate_1 charge_1
      ...
      site_of_N instance_of_N coordinate_N charge_N


      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 and
expected data type
default setting
(if any)
meaning
-grid_detail int2 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.
valuemeaning
0coarse, not advisable for production
1economic, similar to the settings in the old multiflex examples, can result in significant discretization errors
2high, somewhat finer and larger grids than for 1, solvation energies and interaction energies should be largely converged with respect to grid spacing
3very high, even finer and larger innermost grids, ensures very high-quality solvation and interaction energies also for long-range interactions in membrane proteins
4ultra 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 float1 Floating point number that defines the dielectric constant of the reference environment used in force field parametrization.
-eps1set string String giving the prefix of a .pqr file that define the dielectric region 1 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). 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.
-inside1 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 and
expected 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 stringrotlib 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.
-print_chm Flag that triggers detailed output regarding the molecular mechanics terms to stdout (independent on the blab level).
-write_binary_rotlib Flag that triggers writing of the rotamer library as binary archive (using the Boost serialization library).

other program specific options:
option and
expected data type
default setting
(if any)
meaning
-gcem_dir stringgcem_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 float1e99 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 float1e99 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:
    valuemeaning
    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:
    valuemeaning
    0 The site binds ligands from the outer membrane side.
    1 The site binds ligands from the inner membrane side.
  • sitename.est
    An .est file defines forms of a titratable site. Example:

    label NTP NTD1 NTD2 NTD3
    Gmodel 0 10.881 10.881 10.881
    EpsRef 80 80 80 80
    IstrRef 0.1 0.1 0.1 0.1
    proton 1 0 0 0
    center N
    LYS CA 0.21 0.18 0.18 0.18
    LYS HA 0.10 0.10 0.10 0.10
    LYS N -0.30 -0.96 -0.96 -0.96
    LYS HT1 0.33 NaN 0.34 0.34
    LYS HT2 0.33 0.34 NaN 0.34
    LYS HT3 0.33 0.34 0.34 NaN

    keyword / linemeaning
    label labels for the forms
    Gmodel 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.
  • rotlibname.txt/dat
    Squirell backbone dependent rotamer library in ASCII format or as Boost binary archive. The file format is described on the rotamer library website of the Dunbrack group.
  • 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} sitename backgroundname

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 and
expected 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.
-ionicstr float Ionic strength in the exterior region outside the solvent inaccessible volume of backgroundname (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.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} sitename backgroundname

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 and
expected 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

    extended format II (for command line option -fpt):

    site_of_1 instance_of_1 coordinate_1 charge_1
    ...
    site_of_N instance_of_N coordinate_N charge_N


    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} sitename backgroundname

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

    extended format II (for command line option -fpt):

    site_of_1 instance_of_1 coordinate_1 charge_1
    ...
    site_of_N instance_of_N coordinate_N charge_N

    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} sitename backgroundname

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 and
expected 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

    extended format II (for command line option -fpt):

    site_of_1 instance_of_1 coordinate_1 charge_1
    ...
    site_of_N instance_of_N coordinate_N charge_N


    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 and
expected 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.
-inside1 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

    extended format II (for command line option -fpt):

    site_of_1 instance_of_1 coordinate_1 charge_1
    ...
    site_of_N instance_of_N coordinate_N charge_N


    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.

program options and their defaults

program options:

option and
expected data type
default setting
(if any)
meaning
temp float298.15 absolute temperature in K.
gmct_outfile stringmolname.gmct.out Name of the GMCT output file.
curve_dir stringmy_curves 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 int0 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 int0 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.
iprob int0 Flag whether or not (1 = yes, 0 = no) to output occupation probability curves for each instance.
ibprob int0 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 int0 Flag whether or not (1 = yes, 0 = no) to output occupation probability curves for each (binding) form of a site.
icprob int0 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 int0 Flag whether or not (1 = yes, 0 = no) to output the average number of bound ligands of each type bound by a site.
imaxp int0 Flag whether or not (1 = yes, 0 = no) to output the most highly populated instance of a site.
ibmaxp int0 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 int0 Flag whether or not (1 = yes, 0 = no) to output the most highly populated (binding) form of a site.
icmaxp int0 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 int0 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 int0 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 int0 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 int0 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 int0 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 int0 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 int0 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 int0 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.

general options are described above
program specific options:
option and
expected data type
default setting
(if any)
meaning
-CoarseFieldOut string prefix of output file name.fld. This option triggers writing of electrostatic potential for the coarsest (first) grid level in AVS "field" format to name.fld.
-CoarseFieldInit string 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 and
expected 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 and
expected 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.

program specific options:
option and
expected data type
default setting
(if any)
meaning
-epsin floatdielectric constant of a molecular interior region
-epssol float80dielectric constant of an exterior / solvent region
-ionicstr float0.0ionic strength in mol/l
-nopotat setting this flag prevents multiflex from writing .potat files for saving computed electrostatic potentials as restart files for future runs.
-membz z_lower z_upper 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 float20.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.
  • molname.g
    contains the (relative) interaction energies in units of \( \mathrm{e}_{\circ}^{2} / \unicode[Times,serif]{xC5} \) (charge squared per A)

Downloads

MEAD is free software in the sense of the definition given by the Free Software Foundation. You can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 1, or (at your option) any later version. For more details, see the files README, INSTALLATION and COPYING of the distribution. Examples for the program usage with all necessary input files can be found in the directory examples of the MEAD distribution.

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.

extended version of the program suite MEAD download my_mead-2.3.5.tar.bz2

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.

  • electrolyte
  • mobile ions
  • ionic strength
  • ionic distribution
  • Debye-Huckel
  • Debye-Huckel length
  • ionic concentration
  • charge distribution
  • atomic partial charge
  • atomic radius
  • atomic radii
  • solvent accessible surface
  • analytical solvent accessible surface
  • solvent inaccessible volume
  • analytical solvent inaccessible volume
  • Stern layer
  • solvent radius
  • probe sphere radius
  • ion exclusion layer
  • ion exclusion radius
  • ion accessible surface
  • analytical ion accessible surface
  • ion inaccessible volume
  • analytical ion inaccessible volume
  • dielectric constant
  • dielectric permittivity
  • pH value
  • pH
  • pKa
  • pKa value
  • pK(1/2)
  • pK(1/2)
  • pK1/2
  • reduction potential
  • redox potential
  • chemical potential
  • electrochemical potential
  • electrochemical potential gradient
  • chemical potential gradient
  • electrostatic potential
  • electrostatic potential solver
  • electrostatic potential calculation
  • electrostatic potential visualization
  • electrostatic potential distribution
  • electrostatic potential map
  • surface electrostatic potential
  • membrane potential
  • trans-membrane potential
  • transmembrane potential
  • electric membrane potential
  • membrane potential calculation
  • trans-membrane potential calculation
  • transmembrane potential calculation
  • electric membrane potential calculation
  • membrane potential visualization
  • trans-membrane potential visualization
  • transmembrane potential visualization
  • electric membrane potential visualization
  • membrane potential distribution
  • trans-membrane potential distribution
  • transmembrane potential distribution
  • electric membrane potential distribution
  • membrane potential map
  • trans-membrane potential map
  • transmembrane potential map
  • electric membrane potential map
  • membrane potential solver
  • trans-membrane potential solver
  • transmembrane potential solver
  • electric membrane potential solver
  • capacitance
  • capacitance calculation
  • electric distance
  • electric distance calculation
  • membrane
  • lipid membrane
  • cell membrane
  • compartment
  • continuum electrostatics
  • finite-difference method
  • Poisson-Boltzmann equation
  • linearized Poisson-Boltzmann equation
  • PBE
  • LPBE
  • FDPBE
  • Poisson-Boltzmann equation solver
  • linearized Poisson-Boltzmann equation solver
  • PBE solver
  • LPBE solver
  • FDPBE solver
  • OpenDX
  • volume data
  • volumetric data
© 2011-2021 RTU