Why Meadionize?

Meadionize is an improved version of the Autoionize plugin for VMD. Autoionize randomly places sodium and chlorine counterions around a solvated molecule. Autoionize has been designed for molecules with small net charges, where the electrostatic interactions between the molecule and the counterions are relatively weak, and the molecular structure and function are not sensitive to the counterion distribution. That is not the case for highly charged systems, most notably nucleic acids, which are surrounded by a clowd of counterions; correct simulations of these systems require placing ions according to the electrostatic potential of the molecule. Meadionize addresses this problem by placing the ions into the minima of the electrostatic potential map generated by the 'potential' utility of the MEAD program package by Don Bashford by solving the Poisson-Boltzmann equation. Meadionize accepts all MEAD configuration parameters, generates all necessary input files for the potential program, executes it, and uses its output to place the counterions. New: Meadionize version 1.1 can now add divalent ions.

Installation

Create a directory for Meadionize: e.g., under your home directory

    mkdir $HOME/meadionize

Download and save files meadionize.tcl (the main script), ions.top (the ion topology), and pkgIndex.tcl (the package index that tells VMD where to look for the package) into that directory. Add the following two lines to your $HOME/.vmdrc file (create one if you do not have it):

    global env

    lappend auto_path $env(HOME)/meadionize

Then, download and install the MEAD program package following the included installation instructions.

Usage

Like all other VMD plugins, you first need to load Meadionize into VMD. Type the following command in the VMD console (better yet, Tk console):

    package require meadionize

This command actually loads Meadionize into VMD and displays the installed plugin version. Run 'meadionize' with no parameters to get brief help on Meadionize syntax. Meadionize requires the following mandatory parameters: the PSF (structure) file name, the PDB (coordinate) file name, the Charmm parameter file name, positive and negative ion types (as of this point, $\ensuremath{\operatorname{Na}}^+$ , $K^+$ , $\ensuremath{\operatorname{Ca}}^{2 +}$ , $\ensuremath{\operatorname{Mg}}^{2 +}$ , and $\ensuremath{\operatorname{Cl}}^-$ ions are supported), and either ionic strength, or explicit numbers of the positive and negative ions to add. Other parameters are optional; their default values rarely need to be changed, except the fg option (fine grid resolution), which may need to be increased to $0.5 - 1.0$A when running Meadionize on machines with low RAM and/or adding ions to an exceptionally large system. This is an example of Meadionize command line:

    meadionize -psf solvated.psf -pdb solvated.pdb -par par_all27_prot_lipid.prm -ipos na -ineg cl -is 0.1

Since calculating the electrostatic potential map for large molecules can take hours, it is recommended to run Meadionize using VMD in the text-only mode. To do that, one needs to copy the above two command lines into a file (e.g., do_ionize.tcl), add a command 'quit' to the end of the file to tell VMD to stop after executing Meadionize, and run the following command in the UNIX shell (csh or tcsh):

    vmd -dispdev text < do_ionize.tcl >& do_ionize.log &

For bash or ksh, the shell command syntax is slightly different:

    vmd -dispdev text < do_ionize.tcl > do_ionize.log 2>&1 &

To monitor Meadionize in real time, use the following command (any shell):

    tail -f do_ionize.log

Output

Meadionize prints out diagnostic messages about performed steps as well as repeats messages from the potential utility. Often, these messages include the following warnings:

WARNING: SAVanal_calc: vertex found with count = 2

WARNING: SAVanal_calc: vertex found with count = 2

WARNING: SAVanal_calc: vertex found with count = 1

These warnings are harmless, and arise from numerical degeneracies in the calculation of the molecular surface (see You and Bashford, J. Comp. Chem. 16, 743 (1995)). Don Bashford explains these warnings in more details in his post on the Computational Chemistry List.

Algorithm

Meadionize performs the following basic steps:

1. If an ionic strength is requested, Meadionize finds the number of water molecules and calculates the numbers of the positive and negative ions (alternatively, these numbers are given as command line parameters).

2. Prepares necessary input files and calls the 'potential' utility of the MEAD package to solve the Poisson-Boltzmann equation and calculate the electrostatic potential map.

3. In a random order, replaces the water molecules at the electrostatic potential minima (for positive ions) or maxima (for negative ions) with the corresponding ions. Each time, the new ion is placed so that a minimum distance between any two ions, as well as a minimum distance between any ion and themolecule, are maintained.

At the first step, Meadionize computes the numbers of positive and negative ions from two conditions: zero net charge of the system and the ionic strength:

• $n_- z_- - n_+ z_+ = Q_{\ensuremath{\operatorname{mol}}}$ , where $n_-$ ($n_+$ ) and $z_-$ ($z_+$ ) are numbers and charges of the negative (positive) ions, respectively, and $Q_{\ensuremath{\operatorname{mol}}}$ is the net charge of the molecule before adding ions.
• $n_- z^{^{} 2}_- + n_+ z^2_+ = 2 N_{\ensuremath{\operatorname{is}}}$ , where $N_{\ensuremath{\operatorname{is}}}$ is a quantity proportional to the total number of ions, a function of the requested ionic strength.

In the second condition, $N_{\ensuremath{\operatorname{is}}} = C N_A V_{H_2 O}$ , where $C$ is the ionic strength (mol), $N_A = 6.022 \cdot 10^{23}$ is the Avogadro number, and $V_{H_2 O}$is the water volume (L). The latter can be represented as $V_{H_2 O} = N_{H_2 O} / \rho = N_{H_2 O} m_{H_2 O} / \rho_m$ , where $N_{H_2 O}$ is the number of water molecules in the system, $\rho$ is the volume water density, $\rho_m = 0.982 g / \ensuremath{\operatorname{cm}}^3$ is the mass water density (TIP3 water model, Jorgensen et al, J. Chem. Phys. 79, 926 (1983)), and $m_{H_2 O} = 18 \cdot 1.67262 \cdot 10^{- 24} g$ is the water molecule mass. Combining the above, one obtains the following formulas:

• $N_{\ensuremath{\operatorname{is}}} = 0.01846 C N_{H_2 O}$
• $\left. n_+ = [ 2 N_{\ensuremath{\operatorname{is}}} - Q_{\ensuremath{\operatorname{mol}}} z_- \right] / [ z_+ \left( z_+ + z_- \right)]$
• $n_- = \left[ Q_{\ensuremath{\operatorname{mol}}} + n_+ z_+] / z_- \right.$

License

Meadionize is released under the GNU public license. I thank Don Bashford for developing MEAD, which made Meadionize possible. All comments, suggestions, and bug reports are very welcome.