Title

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 options and their defaults

general options:

visualization options (my_xyz_solver):

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 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:

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 int	2	This option influences the detail of the auto-generated grids for the computation of the electrostatic potentials with the finite difference method. A higher value will result in a larger and finer grids. The default value of 2 will in most cases suffice for high quality results, sometimes a value of 3 can be advantageous (especially for membrane proteins, where the apolar environment leads to more far-reaching electrostatic interactions). If user supplied grids exist, the setting affects the automatic adjustment of the innermost grid's size to the size of each site. A higher value will result in a larger spacing of the site to the grid boundary. value meaning 0 coarse, not advisable for production 1 economic, similar to the settings in the old multiflex examples, can result in significant discretization errors 2 high, somewhat finer and larger grids than for 1, solvation energies and interaction energies should be largely converged with respect to grid spacing 3 very high, even finer and larger innermost grids, ensures very high-quality solvation and interaction energies also for long-range interactions in membrane proteins 4 ultra high, extremely fine grids, should not result in significant numeric differences of the result relative to 3, needs very much memory and computation time
-epsext float		Floating point number that defines the dielectric constant of the solvent region and the region specified by elycavset (if applicable).
-epshead float		Floating point number that defines the dielectric constant of the membrane head group region (if applicable / membz is defined).
-epscore float		Floating point number that defines the dielectric constant of the membrane head group region (if applicable / membz is defined).
-epsin1 float		Floating point number that defines the dielectric constant of the region specified by eps1set (if applicable).
-epsin2 float		Floating point number that defines the dielectric constant of the region specified by eps2set (if applicable).
-epsin3 float		Floating point number that defines the dielectric constant of the region specified by eps3set (if applicable).
-eps_ff float	1	Floating point number that defines the dielectric constant of the reference environment used in force field parametrization.
-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.

input files

• molname.pqr
  The receptor structure. The file format is described above under general options and input files.
• molname.ogm and molname.mgm (optional)
  defines the electrostatics grids for the site in the receptor environment and their the model compound in the reference environment (bulk solution), respectively. The file format is described above under general options and input files.
• molname.rot
  defines the sites of molname ans their types For each site, the file contains a line of the format:
  
  site_label site_type
  
  The site label is constructed from sitename-chainid-resid, where chainid and resid correspond to the data found in molname.pqr (chainid set to 0 if absent in molname.pqr). site_type is one of:
  

  value	meaning
  0	ignored as site
  1	flexible and titratable (site_name.est required, see below) user defined conformers are read from molname.pqr flexible parts must be present for each conformer with unique instance IDs differing from 0, For amino acid residues, additional sidechain rotamers are generated for amino acid residues (currently not implemented for Pro), using dihedral angles from the Squirell backbone dependent rotamer library. For each conformer, the different forms defined in site_name.est are generated.
  2	flexible and non-titratable (site_name.est not required) user defined conformers are read from molname.pqr flexible parts must be present for each conformer with unique instance IDs differing from 0, For amino acid residues, additional sidechain rotamers are generated for amino acid residues (currently not implemented for Pro), using dihedral angles from the Squirell backbone dependent rotamer library.
  3	non-flexible and titratable (site_name.est required, see below) For the conformer found in molname.pqr, the different forms defined in site_name.est are generated
  4	quantum mechanically treated site (QM site, site_label.qst required, see below) Conformational flexibility / orientational polarization effects are thought to be considered within the QM treatment. The model energy should also contain any energy terms due to bound ligands, and any other energy contribution that apart from the interactions with the rest of the protein. GCEM The model energy of a QM site is expected to be completely defined in site_label.qst, where site_label is given by sitename-chainid-resid. The coordinates and atomic partial charges of the QM site are taken from molname.pqr or separate .pqr files as described below for site_label.qst.

• molname.con
  defines the connectivity of the sites to the membrane sides. The file is only required for membrane proteins. For each site, the file contains a line of the format:
  
  site_label connectivity
  
  The site label is constructed from sitename-chainid-resid, where chainid and resid correspond to the data found in molname.pqr (chainid set to 0 if absent in molname.pqr). connectivity is one of:
  

  value	meaning
  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 / 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

• 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.

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:

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.

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:

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.

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:

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.

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:

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.

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:

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.
-inside	1	This option defines which membrane side is "inside" in the electrophysiological sense, where the membrane potential is measured at the inside relative to the outside. A value of 0 states that the lower membrane side is inside, while a value of 1 states that the upper side is inside.
-capacitance		This flag triggers the calculation of the capacitance of the system. The capacitance is the accumulated charge in Ampere seconds = Coulomb per membrane potential in Volt. Hence, the capacitance is given by default in Farad (F = As/V).

input files

• backgroundname.pqr and additional .pqr files used to define the dielectric regions
  are structure files in .pqr format. The file format is described above under general options and input files.
• backgroundname.ogm
  defines the electrostatics grids for the site in the receptor environment. The file format is described above under general options and input files.

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

  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.

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:

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 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: nformat 2 format -0.733004200213 3 format -0.043364101728026 4
ipmf int	0	Flag whether or not (1 = yes, 0 = no) to output of a proton motive force (pmf \(= \mu_{\mathrm{H}^{+}}^{\mathrm{out}} - \mu_{\mathrm{H}^{+}}^{\mathrm{in}} - \mathrm{F}\Delta \psi = \bar{\mu}_{\mathrm{H}^{+}}^{\mathrm{out}} - \bar{\mu}_{\mathrm{H}^{+}}^{\mathrm{in}} \)) in the curves instead separate proton chemical potentials and an electric membrane potential. If ipmf is set to 1, the membrane potential is substituted with the pmf and the proton chemical potential on the outside is omitted.
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.

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.

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.

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.

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:

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 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.

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.

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 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.

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 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.

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)