APBS tutorial

Nathan Baker

Washington University in St. Louis
Department of Biochemistry and Molecular Biophysics
Center for Computational Biology

This file is part of APBS.

This program is free software; 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 2 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

Linking APBS statically or dynamically with other modules is making a combined work based on APBS. Thus, the terms and conditions of the GNU General Public License cover the whole combination.

SPECIAL GPL EXCEPTION: In addition, as a special exception, the copyright holders of APBS give you permission to combine the APBS program with free software programs and libraries that are released under the GNU LGPL or with code included in releases of ISIM, Ion Simulator Interface, PMV, PyMOL SMOL, VMD, and Vision. Such combined software may be linked with APBS and redistributed together in original or modified form as mere aggregation without requirement that the entire work be under the scope of the GNU General Public License. This special exception permission is also extended to any software listed in the SPECIAL GPL EXCEPTION clauses by the PMG, FEtk, MC, or MALOC libraries.

Note that people who make modified versions of APBS are not obligated to grant this special exception for their modified versions; it is their choice whether to do so. The GNU General Public License gives permission to release a modified version without this exception; this exception also makes it possible to release a modified version which carries forward this exception.


Table of Contents

1. How to read this tutorial
2. How do I get my structure ready for electrostatics calculations?
2.1. PQR format
2.2. XML format
2.3. Generating PQR files from PDB files (PDB2PQR)
3. How do I visualize the electrostatic potential around my biomolecule?
3.1. VMD
3.2. PyMOL
3.3. PMV
4. How do I calculate a solvation energy?
4.1. Polar solvation
4.2. Apolar solvation
5. How do I calculate a binding energy?
5.1. Solvation energy contributions to binding
5.2. Including Coulombic contributions
5.3. Hey! You don't have any parameters for my ligand!
5.4. A ligand binding example
6. How do I calculate solvation forces?
7. How do I calculate a pKa?
7.1. Introduction
7.2. Application to lysozyme
8. My calculation requires too much memory!
8.1. Parallel calculations: an example
8.2. Asynchronous sequential calculations
9. How do I use APBS with my molecular dynamics program (for MM/PBSA, etc.)?
10. Can I run APBS via the web? (Gemstone)
10.1. Getting Gemstone
10.2. Preparing a structure with PDB2PQR
10.3. Perform an electrostatics calculation with APBS
11. More examples...
12. None of this answered my question...

List of Figures

3.1. Arc repressor isocontours in VMD
3.2. Arc repressor surface potential in VMD
3.3. +/- 1 kT/e electrostatic potential isocontours of FAS2 in PyMOL
3.4. +/- 5 kT/e electrostatic surface potential of FAS2 in PyMOL
4.1. Full solvation energy cycle
5.1. Binding free energy cycle
7.1. pKa perturbation free energy cycle schematic
7.2. Active site of HEWL
8.1. Actin dimer isocontours from parallel focusing
10.1. Gemstone PDB2PQR calculation
10.2. Gemstone APBS/Calculation screen
10.3. Gemstone APBS/Grid screen
10.4. Gemstone APBS/Physics screen
10.5. Gemstone APBS/File I/O screen
10.6. Gemstone APBS/Calculation screen (run complete)

List of Tables

7.1. Model amino acid pKa values for common titratable groups; data from Nielsen et al (see footnote)

List of Examples

4.1. Born ion PQR
4.2. Example born input file

List of Equations

4.1. Born ion polar solvation energy
5.1. Binding free energy formula
5.2. Solvation contributions to the binding free energy
5.3. Coulombic contributions to the binding free energy
5.4. Binding free energy
7.1. Acid dissociation free energy
7.2. Transfer free energy

Chapter 1. How to read this tutorial

This tutorial is designed as a "how to" document to get users acquainted with computational electrostatics calculations using APBS. In order to perform the provided examples, you will need the latest version of APBS (http://apbs.sf.net). Other requirements are listed in the individual sections.

[Important]Important

Note that many of the examples in this tutorial can also be run through the web-based Gemstone service without the need for local installations of the APBS software. More information is available in Chapter 10, Can I run APBS via the web? (Gemstone) .

[Warning]Warning

This document is still under construction and should be completed before the next APBS release. Many of the topics covered in the incomplete sections are demonstrated in the APBS examples directory.

Chapter 2. How do I get my structure ready for electrostatics calculations?

In order to perform electrostatics calculations on your biomolecular structure of interest, you need to provide atomic charge and radius information to APBS. Charges are used to form the biomolecular charge distribution for the Poisson-Boltzmann (PB) equation while the radii are used to construct the dielectric and ionic accessibility functions. This charge and radius information can be provided to APBS in a few different formats which are described in more detail below.

2.1. PQR format

The PQR format provides a very simple way to include parameter information by by replacing the occupancy and temperature columns of a PDB-format structure file with charge ("Q") and radius ("R") information. Unfortunately, the simplicity of this format also limits its extensibility: it can be very difficult to add new atom types and parameters in a generic format without the use of external software such as PDB2PQR. The XML parameter format described below is much easier to modify.

2.2. XML format

The XML format provides a somewhat more complicated format for including parameter information with increased flexibility in formatting, extension, and other modifications. As in the PQR format, atom coordinates are supplemented with charge and radius information. Please see the APBS user guide for complete format specifications.

2.3. Generating PQR files from PDB files (PDB2PQR)

[Note]Note

This tutorial was prepared using PDB2PQR 1.2.1

The PDB2PQR web service and software will convert most PDB files into PQR format with some caveats, in particular the inability to "repair" large regions of missing coordinates. Although PDB2PQR can fix some missing heavy atoms in sidechains, it does not currently have the (nontrivial) capability to model in large regions of missing backbone and sidechain coordinates.

PDB2PQR will also perform hydrogen bond optimization, sidechain rotamer search, limited titration state assignment, ligand parameterization, and APBS input file preparation. Please see http://pdb2pqr.sf.net/ for more details.

As mentioned above, PDB2PQR is discussed in more detailed on the PDB2PQR homepage . Therefore, we will review the minimal steps required to produce a PQR file from a PDB file here. To begin, choose a server from http://pdb2pqr.sf.net/#availability.

2.3.1. Choose the PDB file to convert

Either enter the 4-character PDB ID or accession number (e.g., 1FAS, 1MAH, 1LYS, etc.) or upload your own PDB file.

[Note]Note

Note that, if you choose to enter a 4-character PDB file, PDB2PQR will process all recognizable chains of PDB file as it was deposited in the PDB (e.g., not the biological unit, any related transformations, etc.).

2.3.2. Pick a forcefield

For most applications, the choice is easy: PARSE. This forcefield has been optimized for implicit solvent calculation and is probably the best choice for visualization of protein electrostatics and many common types of energetic calculations for proteins. However, AMBER and CHARMM may be more appropriate if you are attempting to compare directly to simulations performed with those force fields, require nucleic acid support, are simulating ligands parameterized with those force fields, etc.

It is also possible to upload a user-defined forcefield (e.g., to define radii and charges for ligands or unusual residues). Please see the PDB2PQR documentation for more information.

2.3.3. Output naming scheme

This is largely irrelevant to electrostatics calculations but may be important for visualization. When in doubt, choose the "Internal naming scheme" which attempts to conform to IUPAC standards.

2.3.4. Other options

These options fall into two categories: how to build missing atoms (including hydrogens) onto the structure and additional output configuration. Please see http://pdb2pqr.sf.net/ for more details.

Chapter 3. How do I visualize the electrostatic potential around my biomolecule?

There are a number of excellent molecular graphics programs which can be used to visualize electrostatic properties in the context of the biomolecule. Since there are far too many to review in any depth, we will focus on the 3 packages with which we are most familiar.

3.1. VMD

The VMD molecular graphics software package provides support for both the execution of APBS and the visualization of the resulting electrostatic potentials. Documentation on the APBS interface has been provided by the VMD developers at http://www.ks.uiuc.edu/Research/vmd/plugins/apbsrun/. We will supplement this with a basic demonstration of how go from a PDB entry to a plot of structure & potential in VMD using APBS.

[Note]Note

This tutorial was written using VMD 1.8.5.

3.1.1. Generating the PQR

We'll perform this example with arc repressor (PDB ID 1MYK), a DNA binding protein. The first step in the visualization process is generating a PQR file for use with VMD and APBS. Please follow the directions in Section 2.3, “Generating PQR files from PDB files (PDB2PQR)” to generate a PQR file for this structure.

3.1.2. Loading the PQR

Load the PQR file you just created into VMD (FileNew molecule... in the VMD Main window). Adjust the molecule to the desired view. [1] Depending on the force field[2] you used in PDB2PQR, you may notice that VMD is displaying some strange bonds. When a PQR file is loaded into VMD, bond distances are inferred from the PQR radii. These radii were chosen for continuum electrostatics calculations, not visualization. Don't worry if some of the bonds appear a bit strange (hide the hydrogens if it really bothers you).

3.1.3. Running the electrostatics calculation

You're now ready for electrostatics calculations.

  1. Choose ExtensionsAnalysisAPBS electrostatics in the VMD Main window. A new APBS Tool window should appear.

  2. Choose EditSettings... in the APBS Tool window and set the path to your local copy of the APBS executable. Hit OK to close this window.

  3. Select the "0" calculation from the "Individual PB calculations (ELEC):" window. Hit the Edit button. Adjust the APBS settings here as you desire. The defaults are usually fine, although it may be useful to edit the ions to adjust the ionic strength used in the calculation. Hit OK when you're done.

  4. You're now ready to run the calculation. Hit Run APBS in the APBS Tool window. Check the VMD Console window for information about the job while it's running. When the job is finished, a "APBSRun: Load APBS Maps" window should appear. Select "Load files into top molecule" and hit OK.

3.1.4. Electrostatics visualization

3.1.4.1. Isocontour visualization

One of the most popular visualization methods is the isocontour.

  1. First, choose GraphicsRepresentations... from the VMD Main window.

  2. Hit the Create Rep button and change Drawing Method to "Isosurface".

  3. Change Draw from "Points" to "Solid Surface" and Material to "Transparent".

  4. Note that the current isovalue is 0; change this to 1 (or 5 or 10 or whatever) depending on your personal preferences.

  5. To continue the lonstanding tradition of electrostatic potential coloring, choose "ColorID 0" for the Coloring Method.

  6. For the negative isocontour, hit "Create Rep" and select the newly created representation. Change the isocontour value to -1 (or 5 or 10...) and the ColorID to 1.

At this point, you probably have an image that looks something like this (note that I changed the surface to points to make the figure a bit cleaner):

Figure 3.1. Arc repressor isocontours in VMD

Arc repressor isocontours in VMD


3.1.4.2. Visualizing surface potentials

Another popular method of electrostatic potential visualization maps the electrostatic potential to the biomolecular surface. Before proceeding, you may want to delete the two isocontours you just created using the Delete Rep in the Graphical Representations window.

  1. Go to the Graphical Representations window (GraphicsRepresentations... from the VMD Main window) and create a new representation of the molecule with the Create Rep button.

  2. Go to the "Draw style" tab of the Graphical Representations window and change Drawing Method to "Surf"[3] and Coloring Method to "Volume".

  3. Go to the "Trajectory tab" of the Graphical Representations window and change the "Color Scale Data Range:" to -10 to 10 (or whatever your favorite values are).

  4. Based on your version of VMD and your personal preferences, you might want to change the color scale for this image. Go to GraphicsColors... in the VMD Main window and select the "Color Scale" tab from the new "Color Controls" window. The traditional coloring scheme for electrostatics is "RWB" (in the Method menu).

Your molecule now probably looks like this:

Figure 3.2. Arc repressor surface potential in VMD

Arc repressor surface potential in VMD


3.2. PyMOL

The PyMOL molecular graphics software package provides support for both the execution of APBS and the visualization of the resulting electrostatic potentials. We will provide a basic demonstration of how go from a PDB entry to a plot of structure & potential in PyMOL using APBS.

[Note]Note

This tutorial was written using PyMOLX11Hybrid 0.99 on a Mac.

3.2.1. Generating the PQR

We'll perform this example with fasciculin-2 (PDB ID 1FAS), a snake neurotoxin which binds the negatively-charged acetylcholinesterase. Please generate the PQR file using the steps outlined in Section 2.3, “Generating PQR files from PDB files (PDB2PQR)”.

Load the PQR file you created into PyMOL (FileOpen...) and choose your favorite graphical representation.

3.2.2. Performing the electrostatics calcuation

Go to the PluginAPBS Tools... to open the APBS calculation plugin.

  1. Under the "Main" tab of the PyMOL APBS Tools window, select Use another PQR and either browse to (via the Choose Externally Generated PQR: button) or input the path to your PQR file. This step is necessary to ensure you use the radii and charges assigned by PDB2PQR.

  2. Under the "APBS Location" tab of the PyMOL APBS Tools window, either browse to (via the APBS binary location: : button) or input the path to your local APBS binary. It is not necessary to provide a path to the APBS psize.py binary for most biomolecules.

  3. Under the "Temporary File Locations" tab of the PyMOL APBS Tools window, customize the locations of the various temporary files created during the run. This can be useful if you want to save the generated files for later use.

  4. Under the "Configuration" tab of the PyMOL APBS Tools window, hit the Set grid to set the grid spacings. The default values are usually sufficient for all but the most highly charged biomolecules.

  5. Under the "Configuration" tab of the PyMOL APBS Tools window, customize the remaining parameters; the defaults are usually OK[4].

  6. Under the "Configuration" tab of the PyMOL APBS Tools window, hit the Run APBS button to start the APBS calculation. Depending on the speed of your computer, this could take a few minutes. The Run APBS button will become unselected when the calcaultion is finished.

3.2.3. Visualize the electrostatic potential

Before proceeding with the remaining steps, you must load the electrostatic potential data into PyMOL. Under the "Visualization" tab of the PyMOL APBS Tools window, hit the Update button.

3.2.3.1. Electrostatic isocontours

PyMOL makes this step very easy: adjust the positive and negative "Contour" fields to the desired values (usually +/-1, +/-5, or +/-10 kT/e[5]) and hit the Positive Isosurface and Negative Isosurface and Show buttons.

At this point, you probably have a figure that looks something like:

Figure 3.3. +/- 1 kT/e electrostatic potential isocontours of FAS2 in PyMOL

+/- 1 kT/e electrostatic potential isocontours of FAS2 in PyMOL


3.2.3.2. Surface potentials

If you haven't already, hide the isocontours by hitting Positive Isosurface and Negative Isosurface and Hide buttons.

The surface potential is also straightforward to visualize. Set the "Low" and "High" values to the desired values (usually +/-1, +/-5, or +/-10 kT/e) at which the surface colors are clamped at red (-) or blue (+). Check the "Solvent accessible surface" and "Color by potential on sol. acc. surf." buttons to plot the potential on the solvent-accessible (probe-inflated or Lee-Richards) surface[6]. Hit the "Molecular Surface" Show button to load the surface potential.

Figure 3.4. +/- 5 kT/e electrostatic surface potential of FAS2 in PyMOL

+/- 5 kT/e electrostatic surface potential of FAS2 in PyMOL


3.3. PMV

[Warning]Warning

Under construction -- surfaces are currently broken with some PQR structure/force field combinations.



[1] At this point you may notice (to many users' chagrin) that PDB2PQR removed the chain IDs from the structure. Sorry about this; we're trying to fix it in future releases of APBS and PDB2PQR.

[2] Visualization seems to work best with PARSE

[3] FYI, this appears to generate a molecular or Connolly type of surface rather than the solvent-accessible or Lee-Richards surface I personally prefer for this type of visualization...

[4] Although I usually prefer 0.150 concentrations for the +1 and -1 ion species so electrostatic properties aren't overly exaggerated.

[5] Note that PyMOL has a slight units error on the contour values.

[6] In my opinion, the solvent-accessible surface tends to reveal more global features of the surface potential. Tighter surfaces (e.g., van der Waals and molecular or Connolly surfaces) tend to simply map atomic surface onto the biomolecular surface.

Chapter 4. How do I calculate a solvation energy?

Solvation energies are usually decomposed into a free energy cycle as shown in Figure 4.1, “Full solvation energy cycle”. Note that such solvation energies often performed on fixed conformations; as such, they are more correctly called "potentials of mean force". More details on using APBS for the polar and nonpolar portions of such a cycle are given in the following sections.

Figure 4.1. Full solvation energy cycle. This cycle incorporates several processes into the solvation energy (step 1). Step 2 indicates charging of the solute in solution (e.g., inhomogeneous dielectric, ions present). Step 3 indicates the introduction of attractive solute-solvent dispersive interaction interactions (e.g., an integral of Weeks-Chandler-Andersen interactions over the solvent-accessible volume). Step 4 indicates the introduction of repulsive solute-solvent interaction (e.g., cavity formation). Steps 5 and 6 are basically null steps although they could be used to offset unwanted energies added in Steps 3 and 4 above. Finally, Step 6 represents the charging of the solute in a vacuum or homogeneous dielectric environment in the absence of mobile ions.

Full solvation energy cycle. This cycle incorporates several processes into the solvation energy (step 1). Step 2 indicates charging of the solute in solution (e.g., inhomogeneous dielectric, ions present). Step 3 indicates the introduction of attractive solute-solvent dispersive interaction interactions (e.g., an integral of Weeks-Chandler-Andersen interactions over the solvent-accessible volume). Step 4 indicates the introduction of repulsive solute-solvent interaction (e.g., cavity formation). Steps 5 and 6 are basically null steps although they could be used to offset unwanted energies added in Steps 3 and 4 above. Finally, Step 6 represents the charging of the solute in a vacuum or homogeneous dielectric environment in the absence of mobile ions.


4.1. Polar solvation

The full free energy cycle in Figure 4.1, “Full solvation energy cycle” is usually decomposed into polar and nonpolar parts. The polar portion is usually represented by the charging energies in Steps 2 and 6:

Energies returned from APBS electrostatics calculations are charging free energies. Therefore, to calculate the polar contribution to the solvation free energy, we simply need to setup two calculations corresponding to Steps 2 and 6 in Figure 4.1, “Full solvation energy cycle”. However, the electrostatic charging free energies returned by APBS include self-interaction terms. These are the energies of a charge distribution interacting with itself. Such self-interaction energies are typically very large and extremely sensitive to the problem discretization (grid spacing, location, etc.). Therefore, it is very important that the two calculations are performed with identical grid spacings, lengths, and centers, in order to ensure appropriate matching (or "cancellation") of self-energy terms.

One of the canonical examples for polar solvation is the Born ion: a nonpolarizable sphere with a single charge at its center surrounded by an aqueous medium. In the absence of external ions, the polar solvation energy for this sytem is given by Equation 4.1, “Born ion polar solvation energy”.

Equation 4.1. Born ion polar solvation energy. q is the ion charge, a is the ion radius, and the two epsilons denote the internal and external (solution) dielectric constants. This model assumes zero ionic strength.


We can setup a PQR file for the Born ion for use with APBS as shown in Example 4.1, “Born ion PQR” and downloaded here.

Example 4.1. Born ion PQR

REMARK  This is an ion with a 3 A radius and a +1 e charge
ATOM      1   I  ION     1 0.000   0.000   0.000  1.00 3.00
					


We're interested in performing two APBS calculations for the charging free energies in homogeneous and heterogeneous dielectric coefficients. We'll assume the internal dielectric coefficient is 1 (e.g., a vacuum) and the external dielectric coefficient is 78.54 (e.g., water). for these settings, the polar solvation energy expression has the form

where R is the ion radius in Angstroms. One possible APBS input file for this calculation is shown in Example 4.2, “Example born input file” and can be downloaded here.

Example 4.2. Example born input file

# READ IN MOLECULES
read
  mol pqr born.pqr
end


elec name solv # Electrostatics calculation on the solvated state
	mg-manual # Specify the mode for APBS to run
	dime 97 97 97 # The grid dimensions
	nlev 4 # Multigrid level parameter
	grid 0.33 0.33 0.33 # Grid spacing
	gcent mol 1 # Center the grid on molecule 1
	mol 1 # Perform the calculation on molecule 1
	lpbe # Solve the linearized Poisson-Boltzmann equation
	bcfl mdh # Use all multipole moments when calculating the potential
	pdie 1.0 # Solute dielectric
	sdie 78.54 # Solvent dielectric
	chgm spl2 # Spline-based discretization of the delta functions
	srfm mol # Molecular surface definition
	srad 1.4 # Solvent probe radius (for molecular surface)
	swin 0.3 # Solvent surface spline window (not used here)
	sdens 10.0 # Sphere density for accessibility object
	temp 298.15 # Temperature
	calcenergy total # Calculate energies
	calcforce no # Do not calculate forces
end

elec name ref # Calculate potential for reference (vacuum) state
	mg-manual
	dime 97 97 97
	nlev 4
	grid 0.33 0.33 0.33
	gcent mol 1
	mol 1
	lpbe
	bcfl mdh
	pdie 1.0
	sdie 1.0
	chgm spl2
	srfm mol
	srad 1.4
	swin 0.3
	sdens 10.0
	temp 298.15
	calcenergy total
	calcforce no
end

# Calculate solvation energy
print energy solv - ref end

quit
					


Running this example with a recent version of APBS should give an answer of -229.59 kJ/mol which is in good agreement with the -230.62 kJ/mol predicted by the analytic formula above.

Note that the Born example above can be easily generalized to other polar solvation energy calculations. For example, ions could be added to the solv ELEC section, dielectric constants could be modified, surface definitions could be changed (in both ELEC sections!), or more complicated molecules could be examined. Many of the examples included with APBS (e.g., solv and ionize) also demonstrate solvation energy calculations.

[Note]Note

Note that, as molecules get larger, it is important to examine the sensitivity of the calculated solvation energies with respect to grid spacings and dimensions.

4.2. Apolar solvation

Referring back to the full free energy cycle in Figure 4.1, “Full solvation energy cycle”, the nonpolar solvation free energy is usually represented by the energy changes in Steps 3 through 5:

where Step 4 represents the energy of creating a cavity in solution and Steps 3-5 is the energy associated with dispersive interactions between the solute and solvent. There are many possible choices for modeling this nonpolar solvation process. APBS implements a relatively general model described by Wagoner and Baker (PNAS 2006) and references therein. The implementation and invocation of this model is described in more detail in the APBS user guide.

Our basic model for the cavity creation term (Step 4) is motivated by scaled particle theory and has the form

where p is the solvent pressure (press keyword in APBS), Δ V is the solvent-accessible volume of the solute, γ is the solvent surface tension (gamma keyword in APBS), and Δ A is the solvent-accessible surface area of the solute.

Our basic model for the dispersion terms (Steps 3 and 5) follow a Weeks-Chandler-Anderson framework as proposed by Levy, Zhang, and Gallicchio (J Comput Chem, 2002):

where ρ is the bulk solvent density (bconc keyword in APBS), Ω is the problem domain, u(att)(y) is the attractive dispersion interaction between the solute and the solvent at point y with dispersive Lennard-Jones parameters specified in APBS parameter files, and θ(y) describes the solvent accessibility of point y.

[Note]Note

The ability to independently adjust press, gamma, and bconc means that the general nonpolar solvation model presented above can be easily adapted to other popular nonpolar solvation models. For example, setting press and bconc to zero yields a typical solvent-accessible surface area model.

Unlike the Born ion, there are no good simple examples for demonstrating these types of nonpolar calculations. APBS includes several examples of calculations using the apolar model above. Interested readers should examine the alkanes apolar example provided with APBS.

Chapter 5. How do I calculate a binding energy?

In general, implicit solvent models are used to calculation the contribution of solvation to binding free energies. Additional binding free energy contributions (molecular mechanics energies, entropic changes, etc.) must be calculated separately and are not discussed in this tutorial. One exception is the inclusion of intermolecular Coulombic interations; we'll discuss how these can be calculated in APBS below.

Our framework for calculating solvation contributions to binding free energies is shown in Figure 5.1, “ Binding free energy cycle ”.

Figure 5.1.  Binding free energy cycle illustrating binding in terms of transfer free energies from a homogeneous dielectric environment (where interactions are described by Coulomb's law) to an inhomogeneous dielectric environment with differing internal (green) and external (cyan) dielectric constants. The binding (dissociation) free energy is depicted in Step 3.

Binding free energy cycle illustrating binding in terms of transfer free energies from a homogeneous dielectric environment (where interactions are described by Coulomb's law) to an inhomogeneous dielectric environment with differing internal (green) and external (cyan) dielectric constants. The binding (dissociation) free energy is depicted in Step 3.


The binding free energy is given by Equation 5.1, “ Binding free energy formula ”.

Equation 5.1.  Binding free energy formula based on Figure 5.1, “ Binding free energy cycle ”.


The following sections provide more detail on calculating individual terms of this equation.

5.1.  Solvation energy contributions to binding

If we're just interested in calculating the solvation contributions to binding (steps 4 and 2 in Figure 5.1, “ Binding free energy cycle ”), then we simply need to follow the instructions from Chapter 4, How do I calculate a solvation energy? for the complex and isolated components. The solvation energy contribution to the binding is then given by Equation 5.2, “ Solvation contributions to the binding free energy ”.

Equation 5.2.  Solvation contributions to the binding free energy for a two-component complex of molecules "mol 1" and "mol 2". Refer to Figure 5.1, “ Binding free energy cycle ” for numbering.


This solvation energy contribution can be decomposed into polar and nonpolar parts as outlined in Chapter 4, How do I calculate a solvation energy?.

5.2. Including Coulombic contributions

To complete the binding free energy cycle (see Figure 5.1, “ Binding free energy cycle ”), we need to add intermolecular Coulombic contributions to the solvation energy change upon binding to get the total electrostatic/solvent contribution to the binding free energy. In particular, we're interested in the change in Coulombic electrostatic energy upon binding, as given by Equation 5.3, “ Coulombic contributions to the binding free energy ”.

Equation 5.3.  Coulombic contributions to the binding free energy for a two-component complex of molecules "mol 1" and "mol 2". See Equation 5.1, “ Binding free energy formula ” for numbering.


Each of the ΔGcoul quantities in Equation 5.3, “ Coulombic contributions to the binding free energy ” is the sum of pairwise Coulombic interactions between all atoms in the molecule (or complex) for a particular uniform dielectric. In order to combine these Coulombic binding energies with the solvation energies described above, we need to make sure consistent dielectric constants are used. In particular, Coulombic interactions should be calculated using the same uniform dielectric constant as the reference state of the solvation energy above. For example, if solvation energies are calculated for transfering a protein from a medium with uniform dielectric of εin to a medium with internal dielectric εin and external dielectric εout, then Coulombic energies should be calculated using a dielectric of εin. The APBS accessory program tools/manip/coulomb was created to help with the calculation of these individual per-molecule Coulombic energies. Given a PQR file as input, the tools/manip/coulomb program calculates Coulombic energies for a vacuum dielectric (e.g., a uniform dielectric of 1). If the reference dielectric is εin, then all energies returned by tools/manip/coulomb need to be divided by εin.

With the appropriate dielectric constant of εin used for the Coulombic binding energy, a total electrostatic/solvation free energy can be calculated from Equation 5.4, “ Binding free energy ”. Equation 5.1, “ Binding free energy formula ”, Equation 5.2, “ Solvation contributions to the binding free energy ”, and Equation 5.3, “ Coulombic contributions to the binding free energy ” as

Equation 5.4.  Binding free energy


5.3. Hey! You don't have any parameters for my ligand!

PDB2PQR can now (as of version 1.2.0) parameterize ligands, thanks to new features developed in collaboration with the Jens Nielsen group. Please visit the PDB2PQR homepage for more details.

5.4. A ligand binding example

[Warning]Warning

Under construction

Chapter 6. How do I calculate solvation forces?

APBS provides force calculations for both polar and nonpolar solvation following the same procedures used in Chapter 4, How do I calculate a solvation energy?. In general, forces can be obtained by modifying input files used for solvation energy calculations to include calcforce total (for total forces on the entire solute) or calcforce comps to obtain detailed force information for each atom. It is important to note that, like solvation energy calculations, "self-interaction" terms must be removed (cf. the Example 4.2, “Example born input file” example).

Chapter 7. How do I calculate a pKa?

[Important]Important

This tutorial was prepared with Dave Sept as part of a lab for a biomolecular simulation class.

[Note]Note

This tutorial covers Poisson-Boltzmann methods for determining biomolecule pKa values. Other methods for pKa and titration state determination are given in the PDB2PQR examples.

7.1. Introduction

This is a very brief introduction to the concepts behind biomolecular pKas and titration states. More information can be obtained from most biochemistry or biophysics textbooks as well as some of the original articles on pKa evaluation[7].

Recall that the acid dissociation constant Ka describes the dissocation of an acid into its components:

by way of the activities

Under “conditions of ideality,[8] activities can be replaced with concentrations to give

You should also recall that chemical equilibrium constants can be related to free energies by

However, chemists found it easier to use base-10 logarithms than natural logarithms for measuring pH, therefore the pKa is defined as

7.1.1. Amino acid model pKa values

In many calculations, pKa values are assigned based on model values for amino acid side chains to mimic

for an isolated amino acid in solution. Some model pKa values are given in Table 7.1, “Model amino acid pKa values for common titratable groups; data from Nielsen et al (see footnote) ”.

Table 7.1. Model amino acid pKa values for common titratable groups; data from Nielsen et al (see footnote)

Amino acidModel pKa
Arginine13.0
Aspartic acid4.0
Cysteine8.7
C-terminus3.8
Glutamic acid4.4
Histidine6.3
Lysine10.4
N-terminus8.0
Tyrosine9.6


As we'll see in the next section, these model values provide the basis for calculating pKa values in proteins.

7.1.2.  pKa values in proteins

The role of the model pKa values introduced in the previous sections is to move all the chemical (bond making and breaking) complexity of protonation into the model values. In particular, pKa values in proteins are calculated as pertubations of the model compounds according to the free energy cycles shown below.

The pKa for the amino acid in the context of the protein is given by the free energy cycle:

We are interested in determining the unknown ΔaGHA from the known model ΔaGHA,model and the unknown ΔxferGHA and ΔxferGA- according to:

Equation 7.1. Acid dissociation free energy


In general, the quantities ΔxferGHA and ΔxferGA- are obainted from a computational approach. Nearly any free energy calculation method could be used to obtain these energies according to a scheme where the (de)solvation energies of the charged and uncharged amino acids are calculated according to:

Figure 7.1. pKa perturbation free energy cycle schematic

pKa perturbation free energy cycle schematic


Note that these energies have all assumed the same background state of the protein for the calculations; in other words, the same states are used for other titratable groups in the protein for charged and uncharged amino acid states. We'll discuss the implications of this assumption a bit later in this tutorial.

7.1.3. Continuum electrostatics methods for pKa calculations in proteins

Although nearly any free energy method could be used to evaluate the energies of transferring the protonated and unprotonated amino acids from solution into the protein environment, continuum electrostatics offer a (usually) satisfying compromise between accuracy and computational efficiency.

The transfer free energies to be calculated, ΔxferGHA and ΔxferGA-, can be determined from Poisson-Boltzmann (PB) energies. In particular, these energies can be calculated as effective “binding energy calculations” similar to those covered in the APBS examples and tutorial:

Equation 7.2. Transfer free energy


where

  • Gprotein with charged X is the electrostatic energy of the protein with group X bound and all charges on X set to their normal values.

  • Gprotein with uncharged X is the electrostatic energy of the protein with group X bound and all charges on X set to zero.

  • Gcharged X in solution is the electrostatic energy of group X in solution with all charges set to their normal values.

Note that, as with binding energies, Equation 7.2, “Transfer free energy” can be evaluated two ways:

Ways to evaluate transfer free energies

  • Directly in a PB equation by computing the difference of total electrostatic energies (including self energies) from PB calculations for each of the states. In order for this to work, all conformations/grid positions/charge states must be the same in each PB calculation.

  • Indirectly through PB calculations of solvation free energies and Coulomb's law calculations (in a dielectric εp) of electrostatic interaction energies. This method is much more stable.

Both of these methods can be combined, via free energy cycles, to give the desired ΔxferGX. However, when care is taken to use the same grids and conformations for all calculations, the direct method using total electrostatic energies is usually the most stable.

Note that none of the methods discussed above have explicitly allowed for changes in titration state of other groups in the protein during protonation/deprotonation of the acid group of interest. Additionally, none of these methods explicitly provide for conformational changes in the protein coupled to protonation/deprotonoation. As such, we cannot claim to be computing true pKa values with this method. Instead, we are calculating so-called intrinsic pKa values.

7.2. Application to lysozyme

7.2.1. Background

Hen egg white lysozyme (HEWL) is a very popular system for pKa calculations as it has a number of interesting values for its titratable residues. Early pKa work on this enzyme is presented in Tanford C, Roxby R. Interpretation of protein titration curves. Application to lysozyme. Biochemistry. 11 (11), 2192-8, 1972 which also contains the pKa values used in this lab. More recent pKa calculations and a review of some of the methodology can be found in Nielsen JE, Vriend G. Optimizing the hydrogen-bond network in Poisson-Boltzmann equation-based pk(a) calculations. Proteins-Structure Function and Genetics. 43 (4), 403-12, 2001. Finally, the biological relevance of lysozyme is briefly reviewed at the PDB Molecule of the Month site.

HEWL has two active site residues GLU 35 and ASP 52 whose titration states determine the catalytic competency of the enzyme:

Figure 7.2. Active site of HEWL

Active site of HEWL


In particular, the enzyme is only active when ASP 52 is ionized (pKa ≈ 1.2) and GLU 35 is neutral (pKa = 6.3).

Additionally, lysozyme has ASP 66 with a pKa of 6.6 and HIS 15 with a pKa of 5.7.

7.2.2. Preparing the PDB file

Download the PDB entry 2LZT from the PDB; save it as 2LZT-GLU35.pdb. If you have time, you should also visit the PDBSum analysis page as well for additional information about the structure.

[Warning]Problems with explicit water!

It is very important that you remove all explicit water from the PDB file before proceeding. (Why?)

We're going to need to generate a protonated form of GLU 35 to perform our pKa calculations. We will do this with the PDB2PQR web service. PDB2PQR “cleans up” PDB files for electrostatics calculations by performing a number of operations described on the PDB2PQR homepage.

Unless otherwise directed, PDB2PQR adds hydrogens to residues based on model pKa values. Therefore, we will need to specify the titration state of GLU 35 for our pKa calculation by changing the residue name from GLU to GLH using your favorite text editor. Save the result as 2LZT-GLH35.pdb

We're now ready to run PDB2PQR to generate protonated versions of our PDB files. Use the command line version of PDB2PQR or one of the web servers listed on the PDB2PQR homepage to generate protonoated PQR files for for 2LZT-GLU35.pdb and 2LZT-GLH35.pdb. Name your results 2LZT-GLU35.pqr and 2LZT-GLH35.pqr, respectively. Although it is always important to test sensitivity to various force fields, I'd recommend starting with PARSE.

[Note]Note

You can use PDB2PQR to assign titration states with PROPKA but don't do it for the above steps since we need to set the titration states explicit for our calculations.

Recall that we're also going to need the isolated residue for our electrostatics calculations of intrinsic pKa. Use your favorite text editor to extract the entire GLU 35 or GLH 35 residue from 2LZT-GLU35.pqr and 2LZT-GLH35.pqr. Save the results as GLU35.pqr and GLH35.pqr, respectively.

Finally, we'll also need to perform electrostatics calculations on HEWL with the uncharged sidechains. Use your favorite text editor to zero out the charges in 2LZT-GLU35.pqr and 2LZT-GLH35.pqr to create 2LZT-noGLU35.pqr and 2LZT-noGLH35.pqr. This can be done by setting the second-to-last column in the PQR file to zero; e.g.

ATOM    534  N   GLU    35       6.456  16.408  22.487 -0.5163 1.8240
ATOM    535  CA  GLU    35       6.354  15.205  23.349  0.0397 1.9080
ATOM    536  C   GLU    35       6.889  13.941  22.695  0.5366 1.9080
ATOM    537  O   GLU    35       7.705  13.192  23.277 -0.5819 1.6612
ATOM    538  CB  GLU    35       4.906  14.973  23.796  0.0560 1.9080
ATOM    539  CG  GLU    35       4.476  15.932  24.948  0.0136 1.9080
ATOM    540  CD  GLU    35       5.171  15.736  26.255  0.8054 1.9080
ATOM    541  OE1 GLU    35       5.844  14.786  26.570 -0.8188 1.6612
ATOM    542  OE2 GLU    35       5.038  16.682  27.056 -0.8188 1.6612
ATOM    543  H   GLU    35       5.671  16.799  22.056  0.2936 0.6000
ATOM    544  HA  GLU    35       6.856  15.382  24.215  0.1105 1.3870
ATOM    545  HB2 GLU    35       4.318  15.158  23.042 -0.0173 1.4870
ATOM    546  HB3 GLU    35       4.832  14.066  24.142 -0.0173 1.4870
ATOM    547  HG2 GLU    35       4.654  16.857  24.629 -0.0425 1.4870
ATOM    548  HG3 GLU    35       3.500  15.796  25.084 -0.0425 1.4870
				

might become

ATOM    534  N   GLU    35       6.456  16.408  22.487  0.0000 1.8240
ATOM    535  CA  GLU    35       6.354  15.205  23.349  0.0000 1.9080
ATOM    536  C   GLU    35       6.889  13.941  22.695  0.0000 1.9080
ATOM    537  O   GLU    35       7.705  13.192  23.277  0.0000 1.6612
ATOM    538  CB  GLU    35       4.906  14.973  23.796  0.0000 1.9080
ATOM    539  CG  GLU    35       4.476  15.932  24.948  0.0000 1.9080
ATOM    540  CD  GLU    35       5.171  15.736  26.255  0.0000 1.9080
ATOM    541  OE1 GLU    35       5.844  14.786  26.570  0.0000 1.6612
ATOM    542  OE2 GLU    35       5.038  16.682  27.056  0.0000 1.6612
ATOM    543  H   GLU    35       5.671  16.799  22.056  0.0000 0.6000
ATOM    544  HA  GLU    35       6.856  15.382  24.215  0.0000 1.3870
ATOM    545  HB2 GLU    35       4.318  15.158  23.042  0.0000 1.4870
ATOM    546  HB3 GLU    35       4.832  14.066  24.142  0.0000 1.4870
ATOM    547  HG2 GLU    35       4.654  16.857  24.629  0.0000 1.4870
ATOM    548  HG3 GLU    35       3.500  15.796  25.084  0.0000 1.4870
				

7.2.3. Setting up the total electrostatic energy calculations

We will be using focusing calculations to calculate the electrostatic potential and free energies for the systems of interest.

[Warning]Warning

In what follows, we are evaluating total electrostatic free energies -- e.g., energies which contain charge self-interaction terms. We will cancel these self-interaction terms in subsequent steps when we calculate solvation or transfer free energies. Therefore, it is very important that you use the same grid parameters (grid centers, dimensions, spacings, etc.) for every calculation.

Here is a template input that we will use for each of the solvation energy calculations (download it here):

read
	mol pqr compound.pqr # This is the compound for which we will calculate solvation energies
	mol pqr ref.pqr      # This is a compound used as a reference for grid centering
end

elec name inhom          
	mg-auto              # Focusing calculations
	dime 129 129 129     # This is a good grid spacing for this system
	cglen 52.0 66.0 79.0 # These are reasonable coarse grid settings for this system (PDB2PQR-recommended)
	fglen 51.0 59.0 67.0 # These are reasonable fine grid settings for this system (PDB2PQR-recommended)
	cgcent mol 2         # Center the grid on the reference molecule
	fgcent mol 2         # Center the grid on the reference molecule
	mol 1
	lpbe
	bcfl sdh
	pdie 20.00
	sdie 78.54
	srfm smol
	sdens 20.0
	chgm spl2
	srad 1.40
	swin 0.30
	temp 298.15
	calcenergy total
	calcforce no
end

# Print the final energy 
print energy inhom end

quit
				

There are a number of aspects to this input file which are worth noting:

  • In general, compound.pqr will change for each calculation but ref.pqr will not. Choose one molecule to be ref.pqr (2LZT-GLU35.pqr is a good choice) and use it in every calculation.

  • We are using a solute dielectric constant (εp) of 20 (see pdie). This is a common choice for pKa since[9] it is thought to implicitly represent internal relaxation and rearrangement of the solute[10].

You now have enough information to calculate total electrostatic energies for all of the relevant molecules so far: 2LZT-GLU35.pqr, 2LZT-GLH35.pqr, 2LZT-noGLU35.pqr, 2LZT-noGLH35.pqr, GLU35.pqr, and GLH35.pqr. You should be able to construct APBS input files for each of these systems by modifying the template above. Once these input files are constructed, you can run the PB calculation by

$ apbs foo.in | tee foo.out
				

where foo.in is the input file of interest and the output is saved in foo.out.

7.2.4. Setting up the transfer free energy calculations

Recall from our earlier discussion that the transfer free energies can be evaluated in two ways. It is up to you to decide which method to use; both should give comparable results for a sufficiently fine grid. As mentioned above, the direct subtraction of total electrostatic energies is usually the most stable route, assuming identical grids and solute conformations are used for all calculations.

Any differences in results between these two transfer free energy methods are typically due to a lack of convergence in the calculations and can be resolved by decreasing the grid spacing (e.g., increasing the number of grid points).

If you choose to decompose the transfer free energies into solvation and Coulombic energy changes, then you will need to perform two additional electrostatics calculations beyond the APBS calculations outlined in the previous section:

  • First, you will need to re-run all of your APBS calculations with a homogeneous dielectric εs = εp to calculate solvation energies. This can be accomplished by appropriate modification of the input files from the previous section.

  • Second, you will need to calculate Coulombic interaction energies. The APBS program coulomb will evaluate the Coulomb's law electrostatic energy of a set of charges in a vacuum. The program's usage is

    $ coulomb foo.pqr
    							

    where foo.pqr is the particular PQR file you are interested in. Note that this program can also be run with no arguments to obtain an expanded list of options. Don't forget to scale by εp!

7.2.5. Putting it all together

At this point, you should have everything you need to calculate the intrinsic pKa of interest. However, if you get stuck, I've included some example files here that might be helpful:

GLU35.pqr, GLH35.pqr, 2LZT-GLH35.pqr, 2LZT-GLU35.pqr, 2LZT-noGLH35.pqr, 2LZT-noGLU35.pqr

PQR files described above

2LZT-noGLH35.in, 2LZT-GLH35.in, GLH35.in

APBS input files needed to calculate the transfer free energy for GLH directly using total electrostatic energies.

2LZT-noGLH35-vac.in, 2LZT-GLH35-vac.in, GLH35-vac.in

Additional APBS input files needed to calculate the transfer free energy for GLH using polar solvation energies.

2LZT-noGLU35.in, 2LZT-GLU35.in, GLU35.in

APBS input files needed to calculate the transfer free energy for GLU directly using total electrostatic energies.

2LZT-noGLU35-vac.in, 2LZT-GLU35-vac.in, GLU35-vac.in

Additional APBS input files needed to calculate the transfer free energy for GLU using polar solvation energies.

run-apbs.sh

Bash script to run all the APBS input files listed above. Assumes the apbs binary is in your PATH.

run-coul.sh

Bash script to run supporting Coulombic energy calculations for solvation energy-based evaluation of the transfer free energies. Assumes the APBS coulomb binary (usually found in apbs/tools/manip) is in your PATH.

process.sh

Bash script to process APBS output and calculate pKa shifts directly using total electrostatic energies. Assumes results were obtained by running run-apbs.sh, linked above.



[7] The following references provide a good introduction to biomolecular pKa calculations:

[8] Such conditions might happen at extremely low ion and coion concentrations in the absence of other interacting species, etc. In other words, for real biological systems, these conditions are rarely encountered.

[9] Start vigorous waving of hands...

[10] Stop vigorous waving of hands...

Chapter 8. My calculation requires too much memory!

APBS requires approximately 200 B memory per grid point. Memory usage can be predicted before performing the calculations using the tools/manip/psize.py Python script provided with APBS.

If it appears your calculation is going to require more memory than is currently available on your system, you have a few options:

  • APBS calculations can be performed in parallel across multiple processors (hopefully, sharing distributed memory!). This functionality is provided by the mg-para keyword which is described in more detail in the APBS user manual and below.

  • APBS calculations can be broken into a series of smaller, asynchronous runs which (individually) require less memory. This functionality is provided by the mg-para async keyword which is described in more detail in the APBS user manual and below.

  • Submit your calculations through Gemstone to external computational resources.

8.1. Parallel calculations: an example

The actin dimer example provided with the APBS source code (examples/actin-dimer/complex.pqr) is a fairly large system and a good example for parallel focusing calculations. This example will use the actin dimer complex PQR file complex.pqr.

We're going to use an 8-processor parallel calculation to write out the electrostatic potential map for this complex. Each processor will solve a portion of the overall problem using the parallel focusing method (see Baker et al, PNAS, 2001) on a 973 mesh with 20% overlap between meshes for neighboring processors. An example input file for this calculation might look like (download here):

read 
	mol pqr complex.pqr
end

elec name complex
	mg-para
	ofrac 0.1
	pdime 2 2 2
	dime  97  97  97
	fglen 150 115 160
	cglen 156 121 162
	cgcent mol 1
	fgcent mol 1
	mol 1
	npbe
	bcfl sdh
	ion 1 0.150 2.0
	ion -1 0.150 2.0 
	pdie 2.0
	sdie 78.54
	srfm mol
	chgm spl0
	srad 1.4
	swin 0.3
	sdens 10.0
	temp 298.15
	calcenergy total
	calcforce no
	write pot dx pot
end

quit
			

where the pdime 2 2 2 specifies the 8-processor array dimensions, the ofrac 0.1 specifies the 20% overlap between processor calculations, and the dime 97 97 97 specifies the size of each processor's calculation. The write pot dx pot instructs APBS to write out OpenDX-format maps of the potential to 8 files pot#.dx, where # is the number of the particular processor.

Running this input file with an MPI-compiled version of APBS runs 8 parallel focusing calculations, with each calculation generating fine-scale solutions on a different region of the (fglen) problem domain. Note that 8 separate OpenDX files are written by the 8 processors used to perform the calculation. Writing separate OpenDX files allows us to avoid communication in the parallel run and keeps individual file sizes (relatively) small. Additionally, if a user is interested in a specific portion of the problem domain, only a few files are needed to get local potential information.

However, most users are interested in global potentials. For some programs (OpenDX, DataTank), the individual potential files can simply be read into the program separately and the program will reconstruct the global map. Most other programs will require the user to reassemble the global map first; APBS provides the mergedx program for this purpose. mergedx is a simple program that allows users to combine several OpenDX files from a parallel focusing calculation into a single map. This map can be downsampled from the original resolution to provide coarser datasets for fast visualization, etc. For example, the command

$ mergedx 65 65 65 pot0.dx pot1.dx pot2.dx pot3.dx pot4.dx pot5.dx pot6.dx pot7.dx
			

will generate a file gridmerged.dx which has downsampled the much larger dataset contained in the 8 OpenDX files into a 653 file which would be suitable for rough visualization. An example of the mergedx output (see Section 3.1.4.1, “Isocontour visualization” for methods) is shown in

Figure 8.1. Actin dimer isocontours from parallel focusing

Actin dimer isocontours from parallel focusing


Note that downsampling isn't necessary -- and often isn't desirable for high quality visualization or quantitative analysis.

8.2. Asynchronous sequential calculations

The steps described in the previous section can also be performed for systems or binaries which are not equipped for parallel calculations via MPI. In particular, you can add

async n
				

to the ELEC mg-para APBS input file to make the single-processor calculation masquerade as processor n of a parallel calculation.

Scalar maps from asynchronous APBS calculations can be combined using the mergedx program as described above. Currently, energies and forces from asynchronous APBS calculations need to merged manually (e.g., summed) from the individual asynchronous calculation output. This can be accomplished by simple shell scripts.

As a specific example, we can modify the input file above to include a async 0 command in the ELEC statement and thus cause APBS to perform the operations of the first processor in the parallel focusing calculation. The modified input file (available for download here) should look like:

read 
	mol pqr complex.pqr
end

elec name complex
	mg-para
	ofrac 0.1
	pdime 2 2 2
	async 0
	dime  97  97  97
	fglen 150 115 160
	cglen 156 121 162
	cgcent mol 1
	fgcent mol 1
	mol 1
	npbe
	bcfl sdh
	ion 1 0.150 2.0
	ion -1 0.150 2.0 
	pdie 2.0
	sdie 78.54
	srfm mol
	chgm spl0
	srad 1.4
	swin 0.3
	sdens 10.0
	temp 298.15
	calcenergy total
	calcforce no
	write pot dx pot
end

quit
				

This should create an OpenDX-format potential map called pot.dx, corresponding to the output from processor 0 in a parallel focusing calculation. Performing additional APBS calculations with async 1, async 2, ..., async 7 will generate the corresponding OpenDX maps for the remaining processors of the parallel focusing calculations. These can then be reassembled with mergedx as described above.

Chapter 9. How do I use APBS with my molecular dynamics program (for MM/PBSA, etc.)?

Robert Konecny (McCammon group) has developed the iAPBS package which provides a C/C++/FORTRAN interface to APBS for use with AMBER, NAMD, and CHARMM. More information is available from the iAPBS homepage.

APBS also links against developmental versions of the TINKER software package. Public versions of TINKER with APBS support should be available in Fall 2007.

Chapter 10. Can I run APBS via the web? (Gemstone)

Yes, you can run APBS via the web, thanks to the Gemstone service. This section demonstrates how Gemstone can be used to perform most APBS calculations commonly done through the command-line interface. Better yet, it allows you to use NBCR's resources instead of your own!

10.1. Getting Gemstone

You need to download and install the Gemstone extension for the Firefox web browser from http://gemstone.mozdev.org before starting this tutorial.

10.2. Preparing a structure with PDB2PQR

We'll start by preparing a structure with PDB2PQR. First, you need to download a PDB file from the Protein Data Bank. For the purposes of this example, we'll use 1FAS in PDB format (download directly here).

Open the Gemstone interface through the Firefox menu ToolsGemstone. Select the PDB2PQR application through the Gemstone menu UtilsPDB2PQR on the righthand side of the Gemstone screen. Select the PDB file you just downloaded and any other options for PDB2PQR of interest. Your screen should look something like:

Figure 10.1. Gemstone PDB2PQR calculation

Gemstone PDB2PQR calculation


Go ahead and "Run" the Gemstone calculation, check the stdOut and stdErr for any warning messages, and download the resulting PQR file.

10.3. Perform an electrostatics calculation with APBS

We're now ready to perform an APBS calculation of the 1FAS electrostatic potential through Gemstone. Open the ChemistryApbs menu on the righthand side of the Gemstone screen. We'll step through each of the input tabs in the middle pane in the following sections.

10.3.1. Calculation

Move to the Calculation tab in the middle pane.

Most of the "Mathematics" settings can be left as-is; however, feel free to play with the discretization and surface definitions to see how it affects your results.

Under "Molecules", load the 1FAS PQR file you just prepared in the previous section.

Don't select anything under "Energy and Force" for this calculation.

At this point, your screen should look something like this:

Figure 10.2. Gemstone APBS/Calculation screen

Gemstone APBS/Calculation screen


10.3.2. Grid

Next, move to the Grid tab in the middle pane.

Under "General", enter 97 for X, Y, and Z "Grid points".

Under "Coarse Grid", choose 100 for all X, Y, and Z "Lengths". Center the mesh on the molecule.

Under "Fine Grid", choose 80 for all X, Y, and Z "Lengths". Center the mesh on the molecule.

Leave the "Parallel Execution" alone for this calculation.

At this point, your screen should look something like this:

Figure 10.3. Gemstone APBS/Grid screen

Gemstone APBS/Grid screen


10.3.3. Physics

Next, move to the Physics tab in the middle pane.

Choose any values you like here; the defaults should be fine for most visualization calculations.

At this point, your screen should look something like this:

Figure 10.4. Gemstone APBS/Physics screen

Gemstone APBS/Physics screen


10.3.4. File I/O

Next, move to the File I/O tab in the middle pane.

Unless you have other properties you wan to visualize, you can leave these settings as-is.

At this point, your screen should look something like this:

Figure 10.5. Gemstone APBS/File I/O screen

Gemstone APBS/File I/O screen


10.3.5. Running the calculation

Finally, return to the Calculation tab and select "Run APBS". After a while, the results will appear under the Output tab in the middle pane. When the run is complete, you should see something like this:

Figure 10.6. Gemstone APBS/Calculation screen (run complete)

Gemstone APBS/Calculation screen (run complete)


In particular, you should have the option to save the potential file calculated in this APBS run.

10.3.6. Other calculations

Nearly all of the examples in this tutorial can also be performed with Gemstone. Furthermore, input files can be uploaded to use Gemstone as a "portal" to web services at other institutions.

Chapter 11. More examples...

In addition to the relatively simple applications presented here, there are a number of additional examples distributed with the APBS software package (http://apbs.sf.net) which include ionization energies, protein-protein interactions, parallel calculations, etc.

Chapter 12. None of this answered my question...

Please visit the APBS-users mailing list. After you've looked for the answer to your question in the archive, please post it to the mailing list.