AMBER Archive (2009)

Subject: Re: [AMBER] The Right Neutralization Procedure and Ionic Strength Calculation ?

From: Thomas Cheatham III (tec3_at_utah.edu)
Date: Fri Jul 17 2009 - 11:50:15 CDT

> My first question is related to neutralization procedure.
...
> Let's assume that N1>N2 (overal charge is positive).
>
> Here I have two possibilities how to neutralize given system:
>
> #1
> Add just n=N1-N2 Cl- counterions.

"net-neutralization"

> #2
> Add N1 Cl- counterions and N2 Na+ counterions.

"excess salt"

The third possibility is no salt (which although unphysical due to the
overall net-charge works in simulation since the Ewald forces do not
diverge and then there is also this effective net-neutralizing plasma; in
the past AMBER had two options whereby you could smear out the excess
charge or use the implicit plasma).

In 1998, with DNA, I did simulations on a ns time scale and saw little
difference with each scheme. However, recent posts to the reflector
pointed out some anomalies with RNA equilibration where it extended
drastically if no salt was added... I am able to reproduce this and
recommend you do either #1 or #2. Generally I personally try to match
experiment so I try to find an added salt concentration close to
experimental conditions (usually ~200 mM).

The argument for not adding excess salt with polyanions like nucleic
acids relates to sampling (since you would have to better sample the
ionic distribution) and force field parameters (which were flawed and
led to artificial crystallization)... However, as RNA can bind anions,
it is probably wise to have Cl- around too. Therefore I would vote for
#2.
 
> My second question is related to ionic strength (IS).

Ionic strength is complicated since you have the solute charges and also
ion charges that contribute. The way I logically get around this is to
specify that I add net-neutralizing salt to balance the overall charge and
then add a certain concentration of excess salt (like 200 mM). The
overall ionic strength should include all charges; however you could also
report the ionic strength of the added salt.

> the reference volume (for calculation of NaCl concentration) is the overal
> volume of water molecules as it is reported here:
>
> http://www.ks.uiuc.edu/Research/vmd/mailing_list/vmd-l/3943.html

I have a perl script that I plug in the box size and molarity I want and
it, based on volume (**not correcting for the volume of solute**), tells
me how many ions to add. Remember also that by default LEaP produces
water boxes that are slightly less dense than expected. I append the perl
script below.

--tec3

molarity.perl

#!/usr/sbin/perl

$molar = $ARGV[0];
$box_x = $ARGV[1];
$box_y = $ARGV[2];
$box_z = $ARGV[3];
$box_a = $ARGV[4];
$box_b = $ARGV[5];
$box_g = $ARGV[6];

if ($box_x == 0) {
    printf("Usage: molarity box-x box-y box-z alpha beta gamma\n");
    die;
}

if ($box_y == 0) {
    $box_y = $box_x;
}

if ($box_z == 0) {
    $box_z = $box_x;
}

if ($box_a == 0) {
    $box_a = 90.0;
}

if ($box_b == 0) {
    $box_b = $box_a;
}

if ($box_g == 0) {
    $box_g = $box_a;
}

$torad = 2 * 3.141592654 / 360.0;

$rad_a = $box_a * $torad;
$rad_b = $box_b * $torad;
$rad_g = $box_g * $torad;

$angles = 1 - cos($rad_a)*cos($rad_a) - cos($rad_b)*cos($rad_b) - cos($rad_g)*cos($rad_g);
$angles += 2 * cos($rad_a)*cos($rad_b)*cos($rad_g);
$angles = sqrt($angles);

$volume = $box_x * $box_y * $box_z * $angles;

$molecules = 6.022 * $volume * $molar / 10000;

printf(" MOLARITY = %8.3f\n", $molar);
printf(" Box size = %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n", $box_x, $box_y, $box_z, $box_a, $box_b, $box_g);
printf(" Volume = %8.3f\n", $volume);
printf("\n %8.3f molecules are necessary to make a molarity of %6.2f M\n\n", $molecules, $molar);

_______________________________________________
AMBER mailing list
AMBER_at_ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber