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.
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.Here, you find information on general options and input files common to all MEAD programs
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:
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 float | 80 | dielectric constant of an exterior / solvent region (potential, my_xyz_solver) |
-epssol float | 80 | dielectric constant of an exterior / solvent region (multiflex and solvate) |
-epsvac float | 1.0 | dielectric constant of vacuum (solvate) |
-T float | 298.15 | absolute temperature in K |
-solrad float | 1.4 | solvent probe radius for determining a solvent accessible surface / the solvent inaccessible volume |
-sterln float | 2.0 | thickness of an ion-exclusion layer (Stern layer) for determining the ion-inaccessible volume |
-ionicstr float | 0.0 | ionic strength in mol/l |
-kBolt float | 5.984e-06 | Boltzmann constant in units of \(\mathrm{e}_{\circ}^{2}/\left(\unicode[Times,serif]{xC5}\mathrm{K}\right)\) (elementary charges squared per Ângström and Kelvin) This constant must be adjusted if other units of charge, length or temperature are used. |
-econv float | 332.063202 | conversion factor from \(\mathrm{e}_{\circ}^{2}/\left(\unicode[Times,serif]{xC5}\mathrm{K}\right)\) (elementary charges squared per Ângström and Kelvin) to \(\mathrm{kcal}/\mathrm{mol}\) This constant must be adjusted if units of energy and electrostatic potentials other than \(\mathrm{kcal}/\mathrm{mol}\) and \(\mathrm{kcal}/\left(\mathrm{mol}\right)\mathrm{e}_{\circ}\) are needed in the output. |
-conconv float | 6.022214e-04 | conversion factor from \(\mathrm{mol}/\mathrm{l}\) (moles per liter) to \(1/\unicode[Times,serif]{xC5}^{3}\) (particles per cubic Å) This constant must be adjusted if other units of length or concentration are used. |
-epssave_oldway | flag that triggers the old style of averaging the dielectric constant between grid points of differing dielectric constants. Between two potential lattice points in the finite difference method. The new way involves inverse averaging and is similar to a proposal by McCammon. The old way is a simple mean. The new way is significantly more accurate the option to do it the old way is only provided for the sake of reproducing old experimental results. It is not recommended otherwise. | |
-converge_oldway | Revert to the old way of testing for convergence of the successive over-relaxation (SOR) method of solving the finite difference representation of the linearized Poisson Boltzmann equation. The new method gives improved long range accuracy for large lattices, but at a sometimes substantial computational cost. See the NEWS file for further discussion. | |
-blab1 -blab2 -blab3 | flag that controls the verbosity of the programs while writing to stdout, specifying no blab flag at all is least verbose. In the latter case, only essential output is written to stdout. Writing to output files is not affected. |
visualization options (my_xyz_solver):
option 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. |
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)
→
GCEM
is a program for the automated preparation of the necessary input
for
→
GMCT
from a continuum electrostatics / molecular mechanics model.
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.
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:
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.
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 int | 2 |
This option influences the detail of the auto-generated grids for the
computation of the electrostatic potentials with the finite difference method.
A higher value will result in a larger and finer grids.
The default value of 2 will in most cases suffice for high quality results,
sometimes a value of 3 can be advantageous (especially for membrane proteins,
where the apolar environment leads to more far-reaching electrostatic interactions).
If user supplied grids exist, the setting affects the automatic adjustment of the
innermost grid's size to the size of each site.
A higher value will result in a larger spacing of the site to the grid boundary.
|
||||||||||||
-epsext float | Floating point number that defines the dielectric constant of the solvent region and the region specified by elycavset (if applicable). | |||||||||||||
-epshead float | Floating point number that defines the dielectric constant of the membrane head group region (if applicable / membz is defined). | |||||||||||||
-epscore float | Floating point number that defines the dielectric constant of the membrane head group region (if applicable / membz is defined). | |||||||||||||
-epsin1 float | Floating point number that defines the dielectric constant of the region specified by eps1set (if applicable). | |||||||||||||
-epsin2 float | Floating point number that defines the dielectric constant of the region specified by eps2set (if applicable). | |||||||||||||
-epsin3 float | Floating point number that defines the dielectric constant of the region specified by eps3set (if applicable). | |||||||||||||
-eps_ff float | 1 | Floating point number that defines the dielectric constant of the reference environment used in force field parametrization. | ||||||||||||
-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. | |||||||||||||
-inside | 1 | This option defines which membrane side is "inside" in the electrophysiological sense, where the membrane potential is measured at the inside relative to the outside. A value of 0 states that the lower membrane side is inside, while a value of 1 states that the upper side is inside, where lower and upper refer to lower and higher values of the z coordinate, respectively. | ||||||||||||
-capacitance | This flag triggers the calculation of the capacitance of the system. The capacitance is the accumulated charge in Ampere seconds = Coulomb per membrane potential in Volt. Hence, the capacitance is given by default in Farad (F = As/V). The calculation of the capacitance requires tighter convergence criteria for the electrostatic potentials, which increases the required computation time. The capacitive energy is only needed if there are multiple global conformations. The difference between the the capacitive energy terms of different global conformations is negligibly small under normal conditions of naturally occurring values of the membrane potential, but depends quadratically on membrane potential. Caution: The capacitance and hence the capacitive energy depend on the system size (mainly on the area of the covered membrane region). Therefore, the same grid size must be used for the innermost grid in the calculation of the capacitance for all global conformations (normally automatically taken care of). |
program specific options for the molecular mechanics part:
option 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 string | rotlib | Prefix of the name of a file containing the Squirell backbone dependent rotamer library (tested with 2002 and 2010 versions). First, it is attempted to read rotlibname.dat as binary boost archive. If the binary file is not found, it is attempted to read rotlibname.txt. The file is not read if there are only non-flexible sites (according to the definitions in molname.rot). |
-skip_stat_mmint | Flag that triggers the omission of molecular mechanics contributions to the intrinsic energies of each non-flexible site (according to the specification of the site type in molname.rot). | |
-skip_stat_mmg | Flag that triggers the omission of molecular mechanics contributions to each interaction energy that involves a non-flexible site (according to the specification of the site type in molname.rot). | |
-empirical_energy | Flag that triggers the use of an empirical energy function for the part of the intrinsic energy that involves atoms of the site itself and nearby backbone residues that are thought to determine the rotamer propensity in the backbone dependent rotamer library. The empirical energy is defined as \(-\beta^{-1} \ln\left[ p\right]\), where p is the propensity for the corresponding sidechain rotamer in the rotamer library. | |
-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 string | gcem_input | Name of the subdirectory to which gcem will write input files for the continuum electrostatics solvers and the corresponding job script(s). In addition, the directory will contain some deprecated files that are no longer used but may be helpful in debugging while making changes to the source code. |
-mead_path string | Path to the root directory of the MEAD distribution used for the calls to the continuum electrostatics solvers in the job script(s). The option is useful if the solver jobs are run on a different system, for example a remote computer cluster. | |
-filter_ctof float | 1e99 | Intrinsic energy cutoff (in kcal/mol) for eliminating conformers with excessively high intrinsic energy, which makes them unlikely to be populated significantly in equilibrium. A conformer (maybe a sidechain rotamer) consisting of N instances (may be different binding forms of the sidechain) is eliminated if the minimum intrinsic energy of any instance of the conformer is larger than the minimum intrinsic energy of all instances of the site plus filter_ctof. |
-filter_ctof_gs float | 1e99 | Energy cutoff (in kcal/mol) for eliminating conformers with excessively high energy, which makes them unlikely to be populated significantly in equilibrium. The elimination criterion is almost identical to the Goldstein criterion of dead end elimination, with the exception that the hypothetical minimum energy must be larger by filter_ctof_gs rather than by 0. A conformer (maybe a sidechain rotamer) consisting of N instances (may be different binding forms of the sidechain) is eliminated if all of its instances fulfill the modified Goldstein criterion. |
-cap_max_nin int | Flag that triggers the reduction of the number of instances of any site to a number ≤ cap_max_nin. Keep no more than cap_max_nin instances of a site with lowest sum of intrinsic energy and sum of minimum interaction energies with all other sites. The method is equivalent to repeated application of the modified Goldstein criterion to all conformers with successively decreasing filter_ctof_gs until the number of instances is equal to or smaller than cap_max_nin. | |
-print_debug | Flag that triggers detailed output regarding the read and generated data to stdout (independent on the blab level). The flag is mainly intended for verification purposes while debugging new program features. |
value | meaning |
0 | ignored as site |
1 |
flexible and titratable (site_name.est required, see below)
user defined conformers are read from molname.pqr flexible parts must be present for each conformer with unique instance IDs differing from 0, For amino acid residues, additional sidechain rotamers are generated for amino acid residues (currently not implemented for Pro), using dihedral angles from the Squirell backbone dependent rotamer library. For each conformer, the different forms defined in site_name.est are generated. |
2 |
flexible and non-titratable (site_name.est not required)
user defined conformers are read from molname.pqr flexible parts must be present for each conformer with unique instance IDs differing from 0, For amino acid residues, additional sidechain rotamers are generated for amino acid residues (currently not implemented for Pro), using dihedral angles from the Squirell backbone dependent rotamer library. |
3 |
non-flexible and titratable (site_name.est required, see below)
For the conformer found in molname.pqr, the different forms defined in site_name.est are generated |
4 |
quantum mechanically treated site (QM site, site_label.qst required, see below)
Conformational flexibility / orientational polarization effects are thought to be considered within the QM treatment. The model energy should also contain any energy terms due to bound ligands, and any other energy contribution that apart from the interactions with the rest of the protein. GCEM The model energy of a QM site is expected to be completely defined in site_label.qst, where site_label is given by sitename-chainid-resid. The coordinates and atomic partial charges of the QM site are taken from molname.pqr or separate .pqr files as described below for site_label.qst. |
value | meaning |
0 | The site binds ligands from the outer membrane side. |
1 | The site binds ligands from the inner membrane side. |
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 / line | meaning |
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 |
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. |
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.
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 settings are specified with command line options. Input files are described below. The program is called with:
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.
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). |
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.
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 settings are specified with command line options. Input files are described below. The program is called with:
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.
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). |
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):
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.
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 settings are specified with command line options. Input files are described below. The program is called with:
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.
general options are described above
backgroundname.diel
is a file defining the dielectric regions
format of a single line corresponding to a dielectric region:
backgroundname.ely
is a file defining the electrolyte regions
format of a single line corresponding to an electrolyte region:
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):
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.
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 settings are specified with command line options. Input files are described below. The program is called with:
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.
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. |
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):
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.
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 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:
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.
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. | |
-inside | 1 | This option defines which membrane side is "inside" in the electrophysiological sense, where the membrane potential is measured at the inside relative to the outside. A value of 0 states that the lower membrane side is inside, while a value of 1 states that the upper side is inside. |
-capacitance | This flag triggers the calculation of the capacitance of the system. The capacitance is the accumulated charge in Ampere seconds = Coulomb per membrane potential in Volt. Hence, the capacitance is given by default in Farad (F = As/V). |
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):
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.
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 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:
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:
option and expected data type | default setting (if any) | meaning | |||||||||
temp float | 298.15 | absolute temperature in K. | |||||||||
gmct_outfile string | molname.gmct.out | Name of the GMCT output file. | |||||||||
curve_dir string | my_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 int | 0 | Program verbosity as described under general program options. | |||||||||
nformat int |
Keyword that triggers the reading of a block of format descriptions
for the chemical potentials (same order as in gmct_outfile.
The number gives the number of format lines in the block.
Each line consists of the keyword format followed by a conversion factor
and the number of decimal digits after the comma to which the converted
chemical potential, pmf or membrane potential should be rounded.
If there is a proton chemical potential and an electron chemical potential
which should appear in the form of a pH value and a reduction potential in V in the
program output,
the block could look like this:
|
||||||||||
ipmf int | 0 | Flag whether or not (1 = yes, 0 = no) to output of a proton motive force (pmf \(= \mu_{\mathrm{H}^{+}}^{\mathrm{out}} - \mu_{\mathrm{H}^{+}}^{\mathrm{in}} - \mathrm{F}\Delta \psi = \bar{\mu}_{\mathrm{H}^{+}}^{\mathrm{out}} - \bar{\mu}_{\mathrm{H}^{+}}^{\mathrm{in}} \)) in the curves instead separate proton chemical potentials and an electric membrane potential. If ipmf is set to 1, the membrane potential is substituted with the pmf and the proton chemical potential on the outside is omitted. | |||||||||
iprob int | 0 | Flag whether or not (1 = yes, 0 = no) to output occupation probability curves for each instance. | |||||||||
ibprob int | 0 | Flag whether or not (1 = yes, 0 = no) to output occupation probability curves for each possible number of ligands bound to a site (permutation of the possible numbers of bound ligands of each type). | |||||||||
ifprob int | 0 | Flag whether or not (1 = yes, 0 = no) to output occupation probability curves for each (binding) form of a site. | |||||||||
icprob int | 0 | Flag whether or not (1 = yes, 0 = no) to output occupation probability curves for each global conformation and each conformer of a site. Conformers are expected be ordered sequentially in gmct_outfile (automatically done if the input for GMCT was generated with GCEM ) , where each conformer possesses the same (binding) forms. | |||||||||
ibmean int | 0 | Flag whether or not (1 = yes, 0 = no) to output the average number of bound ligands of each type bound by a site. | |||||||||
imaxp int | 0 | Flag whether or not (1 = yes, 0 = no) to output the most highly populated instance of a site. | |||||||||
ibmaxp int | 0 | Flag whether or not (1 = yes, 0 = no) to output the most highly populated number of bound ligands of each ligand type for a site. | |||||||||
ifmaxp int | 0 | Flag whether or not (1 = yes, 0 = no) to output the most highly populated (binding) form of a site. | |||||||||
icmaxp int | 0 | Flag whether or not (1 = yes, 0 = no) to output the most highly populated global conformation and the most highly populated conformer of a site. | |||||||||
ihalf int | 0 | Flag whether or not (1 = yes, 0 = no) to output the midpoints for each transition between the possible pairs of instances of a site. One chemical potential is varied, while the others have fixed values. The midpoint is the value of the varied chemical potential at which the population of the two instances is equal at fixed values of all other chemical potentials. | |||||||||
ibhalf int | 0 | Flag whether or not (1 = yes, 0 = no) to output the midpoints for each transition between consecutive binding numbers of each ligand type and each site. The midpoint is the chemical potential of the corresponding ligand type at which the population of the consecutive binding numbers is equal at fixed values of all other chemical potentials, the pmfs and the membrane potential. An example application of this option is the determination of \(\mathrm{p}K_{1/2}\) values, which in general depend on the chemical potentials of all other ligands present in the system and a possibly present membrane potential. | |||||||||
ifhalf int | 0 | Flag whether or not (1 = yes, 0 = no) to output the midpoints for each transition between the possible pairs of (binding) forms of a site. One chemical potential is varied, while the others have fixed values. The midpoint is the value of the varied chemical potential at which the population of the two (binding) forms is equal at fixed values of all other chemical potentials. | |||||||||
ichalf int | 0 | Flag whether or not (1 = yes, 0 = no) to output the midpoints for each transition between the possible pairs of global conformations and of conformers of a site. One chemical potential is varied, while the others have fixed values. The midpoint is the value of the varied chemical potential at which the population of the two global conformations or site conformers is equal at fixed values of all other chemical potentials. | |||||||||
iginst int | 0 | Flag whether or not (1 = yes, 0 = no) to output the free energy of the instances of a site calculated from \(G = -\beta^{-1} \ln\left[p\right]\), where \(p\) is the occupation probability of the instance. | |||||||||
igbind int | 0 | Flag whether or not (1 = yes, 0 = no) to output the free energy of each possible number of ligands bound to a site (permutation of the possible numbers of bound ligands of each type) calculated from \(G = -\beta^{-1} \ln\left[p\right]\), where \(p\) is the occupation probability of the binding number. | |||||||||
igconf int | 0 | Flag whether or not (1 = yes, 0 = no) to output the free energy of each global conformation and each site conformer calculated from \(G = -\beta^{-1} \ln\left[p\right]\), where \(p\) is the occupation probability of the global conformation or site conformer. | |||||||||
itraj int | 0 | Flag whether or not (1 = yes, 0 = no) to write assembled structures of molname with all sites in their most highly populated instances for each permutation of chemical potential values. |
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.
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:
where molname is the name of the site used as prefix for the corresponding .pqr file.
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.
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:
where molname is the name of the site used as prefix for the corresponding .pqr file.
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.
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:
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.
This program calculates the electrostatic potential due to the charge distribution of molname
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:
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. |
This program calculates the electrostatic solvation energy of a molecule in a solvent.
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:
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 calculates the electrostatic part of the transfer energy for bringing a compound from the bulk solvent to a molecular environment.
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:
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 is a program for the automated preparation of the necessary input for continuum electrostatic calculations on binding equilibria within a traditional two-state model.
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:
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 float | dielectric constant of a molecular interior region | |
-epssol float | 80 | dielectric constant of an exterior / solvent region |
-ionicstr float | 0.0 | ionic 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. |
This program calculates titration curves with an approximate analytical method (reduced site method).
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 specific options:
-cutoff float | 20.5 | Cutoff for the reduced site method. |
-pHrange pH_min pH_max | -20 30.1 | pH range for the titration calculation |
-dry | Flag that causes redti to do a "dry run" in which it prints the number of sites to be included in the reduced site calculation at each pH point. This is useful for checking whether a calculation will require a prohibitive amount of CPU time since CPU time will go exponentially in the reduced site number. |
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 |
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.
© 2011-2023 RTU