Re: setting the polarizabilities for the Drude force field

From: Kenno Vanommeslaeghe (
Date: Fri Aug 02 2013 - 11:01:26 CDT

Hi Frank and Hadi,

I just had a talk with the lead developer of the CHARMM Drude polarizable
force field to settle this.

On 08/01/2013 07:48 AM, Frank Haverkort wrote:
> In the PSF file, the drude charge q_D and the polarizability alpha are
> specified. On the other hand, the force constant for the drude particle
> k_D is specified in the force field parameter file. But since these
> quantities are related via alpha = q_D^2 / k_D, why is it not sufficient
> to define only two of them in the input files? Now it seems that I
> manually have to calculate the force constant k_D from a given alpha and q_D.

The force constant from the parameter file is used, and is currently the
same for all Drude particles. When CHARMM builds an XPLOR psf based on the
CHARMM topology file, it uses the alpha to determine the *charge* on the
drude particle. So it would indeed appear that having both the charge and
the alpha in the psf would be redundant. Probably it's just the format
specifications. Perhaps the more important question is which of the two
NAMD actually uses. If I were to make a guess, I would say the charge, but
maybe a NAMD developer can answer that question for sure.

> This statement seems to contradict the statement I found on this mailing
> list, that the Drude force constants are specified in the force field
> parameter file.

The force constants *are* specified in the parameter file, but in the
CHARMM Drude polarizable FF, we happen to constrain ourselves to using
just one force constant for all Drude particles.

On 08/01/2013 02:39 PM, Hadi Dinpajooh wrote:
> The PSF file contains information about q_D, polarizability (alpha),
> lone-pairs, thole parameter, etc., but it does not contain any information
> about force constant. You should add force constant information to the
> parameter file when you run NAMD simulations.

If we're talking about the CHARMM Drude polarizable FF, all necessary the
force constant information is already in the parameter file in the from of
a wildcard:
X DRUD 500.00 0.000
Versions of NAMD that are new enough to run Drude simulations correctly do
recognize this syntax.

> You do need to calculate the
> force constant, k_D, for given values of alpha and q_D.

Nope, the second option is the correct one:
> However, you can
> calculate q_D for given values of alpha and k_D instead (I believe this is
> how CHARMM generates PSF files).

> If you use a very stiff force constant,
> then the displacement of the Drude particle remains much smaller than any
> interatomic distance. Also note that if you want to use a force constant
> of 1000 kcal/(mol Ang^2), you should set it to 500 kcal/(mol Ang^2) in the
> parameter file due to the potential function used for bond stretches, i.e.
> V(bond) = Kb(b - b0)**2.

You should use the force constant that came with the force field. Choosing
something else may lead to problems.



This archive was generated by hypermail 2.1.6 : Wed Dec 31 2014 - 23:21:29 CST