Personal tools

DelPhi

How to automate the calculation of electrostatic potentials with the Poisson-Boltzmann solver DelPhi

Biskit contains a wrapper for the DelPHI Poisson-Boltzmann solver from the Honig lab, with the unimaginative name Biskit.Delphi. It's actually more than just a wrapper and can (optionally) automate a whole pipeline from structure repair to hydrogen bond-optimization and the assignment of complete Amber all atom partial charges. An application of this wrapper is the automatic calculation of electrostatic free energies of binding with the Biskit.Dock.DelphiBindingEnergy class (see also the end-user script Dock/delphiBinding.py) .

NEW in Biskit 2.4

Background

Two major inconveniences of DelPHI are addressed by this pipeline:

  1. the default charge files bundled with DelPhi are outdated and always have to be manually adapted to account for the proper handling of terminal residues in proteins or DNA. The Biskit wrapper maps state-of-the-art partial atomic charges from Amber topology files into the structure, automatically discovers residues with deviating atom content (e.g. terminal residues) and creates a custom DelPhi charge file.
  2. Independent DelPhi calculations for different conformations or structures can only be compared if the grid is exactly identical; Delphi can determine grid dimensions automatically but these dimensions cannot be reproduced for another system. So we are back to a manual specification of grid dimensions (center, number of points, distance between points) for most but the most simple applications. For the Biskit wrapper, I have "reverse-engineered" the DelPhi way of determining automatic grid dimensions (with the perfil option), this grid is then specified explicitely in the DelPhi input file and can be copied between different instances of the Delphi wrapper.

The current default workflow of this wrapper is:

  1. take input model/structure, remove hydrogens
  2. optionally: cap chain breaks and probable false N- and C-termini with NME and ACE residues to cancel artificial charges (see Biskit.PDBCleaner)
  3. add and optimize hydrogens with the reduce program (see Biskit.Reduce )
  4. adapt residue and atom names to Amber conventions
  5. match each residue (by atom content) to a residue from a list of Amber residue topology files and assign partial charge to each atom (see Biskit.AtomCharger)
  6. Create custom delphi charge file with Amber partial charges
  7. Determine center and dimensions of the grid used in Delphi calculation
  8. Run Delphi in temporary folder
  9. parse result energies into result dictionary

Usage

>>> inputmodel = PDBModel( 'mystructure.pdb' )
>>> D = Delphi( inputmodel )
>>> result = D.run()
>>> result
{'scharge' :  1.4266   # surface charge
 'egrid' :  9105.51    # total grid energy
 'ecoul' :  -9849.664  # couloumb energy
 'erxn'  :  -664.7469   # corrected reaction field energy
 'erxnt' :  -21048.13  # total reaction field energy
 'eself' :  -20383.39 }  # self reaction field energy

The delphi potential map can be rescued with the parameter f_map:

>>> D =  Delphi( inputmodel, f_map='delphi_mapout.phi' )

This will save the map to the file delphi_mapout.phi. Otherwise, the map is deleted together with the temporary folder.

Note

All energy values are in units of kT. The important terms for the calculation of (free) energy differences are 'egrid', 'erxn' and 'ecoul'.

Grid dimensions

If no further options are given, the Delphi wrapper will determine grid dimensions that are centered on the molecule and lead to at most 60% "filling" of the grid in any of the x, y or z dimensions. In other words, the longest dimension of the solute along x, y or z axis will still have a 20% margin to the grid boundary. This fill-factor can be overriden with the 'perfil' parameter.

The density of the grid is, by default, 2.3 points per Angstroem and can be adjusted with the 'scale' parameter. The number of grid points (gsize) will be calculated accordingly.

For a given combination of scale and perfil, this calculation will give approximately (but not exactly) the same grid dimensions as the one the delphi program would determine by itself. Minor differences may lead to two grid points more or less. This can also lead to some changes in resulting energies.

For more control over grid dimensions, you should call the method Delphi.setGrid() before Delphi.run(). For example:

>>> D = Delphi( inputmodel )
>>> D.setGrid( scale=1.2, perfil=80 )

... will fix a grid centered on the molecule, with 1/1.2 A grid spacing and at least 10% distance margin to the boundary. Another example:

>>> D = Delphi( inputmodel )
>>> D.setGrid( acenter=(0,0,0), scale=1.0, gsize=100 )

... circumvents any calculation and fixes a 100 x 100 x 100 grid centered on 0 and with 1/1.0 A grid spacing -- a box of 100 A x 100 A x 100 A. If the scale parameter is not given, it will default to whatever was specified at the Delphi() constructor (and from there default to 2.3).

Delphi.setGrid() returns the three grid parameters as a dictionary for later re-use. So if you want to use the same grid dimensions for two different structures you can first calculate the dimensions based on one structure and then apply the same dimensions to another Delphi run:

>>> D1 = Delphi( complex )
>>> grid1 = D1.setGrid( scale=2.4, perfil=60 )
>>> print grid1
{'acenter': (0.1, 10., -0.5), 'scale': 2.400, 'gsize': 99 }
>>> ligand = complex.takeChains( [1] ) ## extract part of structure
>>> D2 = Delphi( ligand )
>>> D2.setGrid( **grid1 )

Note: The '**' syntax unpacks the dictionary as keyword parameters into the method.

Customization

The default delphi parameter file can be replaced (parameter template). Note though that you should keep the place holders for input and output files. The default parameter file is taken from:

Biskit/data/delphi/delphi_simple.prm

As always, the place holders (e.g. %(salt)f ) in this file are replaced by the value of a variable of the same name (e.g. salt) within the name space of the Executor instance. In other words:

>>> D = Delphi( inputmodel, salt=0.2, ionrad=2.5 )

...will override the default values for salt and ionradius for which there are the place hoders %(salt)f and %(ionrad)f in the template file.

The default handling and matching of atomic partial charges can be modified by:

  • providing a ready-made Delphi charge file (parameter: f_charges)
  • or providing an alternative list of Amber topology files from which residues are looked up by their atom content (parameter: topologies)
  • or circumvent charge matching (parameter: addcharge=False) and provide a input PDBModel with an atom profile 'partial_charge' that is used instead

Note:

Command configuration: biskit/Biskit/data/defaults/exe_delphi.dat

As always this file can be copied to ~/.biskit/ in your home directory and then be adapted to your local environment

Dependencies:

Two external programs are required for this wrapper:

Test:

Run the built-in unit test to verify your installation:

cd biskit/Biskit
python delphi.py

Should (after lots of other output) yield:

----------------------------------------------------------------------
Ran 2 tests in 12.987s

OK
Help:

The complete list of parameters is documented in the source code (at the moment, the online API documentation is lagging behind, you can generate your own with epydoc though). All methods are documented python-style and can be interrogated at the python prompt:

>>> d = Delphi( PDBModel() )  # create Delphi instance with an empty structure
>>> print d.__init__.__doc__

Which will yield something like that:

@param model: structure for which potential should be calculated
@type  model: PDBModel
@param template: delphi command file template [None=use default]
@type  template: str
@param f_radii: alternative delphi atom radii file [None=use default]
@type  f_radii: str
@param topologies: alternative list of residue charge/topology files
                   [default: amber/residues/all*]
@type  topologies: [ str ]
@param f_charges: alternative delphi charge file
                  [default: create custom]
@type  f_charges: str
@param f_map   : output file name for potential map [None= discard]
@type  f_map   : str
@param addcharge: build atomic partial charges with AtomCharger
                  [default: True]
@type  addcharge: bool

@param protonate: (re-)build hydrogen atoms with reduce program (True)
                  see L{Biskit.Reduce}
@type  protonate: bool
@param autocap: add capping NME and ACE residues to any (auto-detected)
                false N- or C-terminal and chain breaks (default: False)
                see L{Biskit.Reduce} and L{Biskit.PDBCleaner}
@type  autocap: bool

@param indi: interior dilectric (4.0)
@param exdi: exterior dielectric (80.0)
@param salt: salt conc. in M (0.15)
@param ionrad: ion radius (2)
@param prbrad: probe radius (1.4)
@param bndcon: boundary condition (4, delphi default is 2)
@param scale:  grid spacing (2.3)
@param perfil: grid fill factor in % (for automatic grid, 60)

@param kw: additional key=value parameters for Executor:
@type  kw: key=value pairs
::
  debug    - 0|1, keep all temporary files (default: 0)
  verbose  - 0|1, print progress messages to log (log != STDOUT)
  node     - str, host for calculation (None->local) NOT TESTED
                  (default: None)
  nice     - int, nice level (default: 0)
  log      - Biskit.LogFile, program log (None->STOUT) (default: None)