[MMTK] changing charge of one atom

[Konrad Hinsen]

On Mar 19, 2007, at 11:27, Davide Ricci wrote:

> I'd like to be able to change the charge of a molecule in my  
> Universe. For example, suppose I created an universe of 128  
> molecules of water, equilibrated and at a certain density. Now I  
> want to perform a MD and, every 200 steps, take a snapshot of the  
> systems, select one or more molecules, change the charge parameter  
> of the FF and calculate the potential energy. Is it possible ?
> How can I change the charge of a molecule, once I've created the  
> Universe ?

This is possible, though a bit tricky. First of all, "charge" is a  
force field parameter, not a physical property of universal meaning.  
There are force fields that have no parameter called "charge", and  
those who do have it may interpret it differently. The way in which  
charges are assigned to each atom also depends on the force field -  
which means that any approach to modifying the charge is also force  
field dependent.

Here is a simple example for water:

from MMTK import *
from MMTK.ForceFields import SPCEForceField

universe = InfiniteUniverse(SPCEForceField())
w1 = Molecule('water', position = Vector(0., 0., 0.))
w2 = Molecule('water', position = Vector(1., 0., 0.))
universe.addObject([w1, w2])

print universe.energyTerms()['electrostatic']

w1.spce_charge = copy(w1.spce_charge)
for i in range(10):
     for atom in w1.atomList():
         w1.spce_charge[w1.getReference(atom)] = 1.1*atom.charge()
     universe.setForceField(SPCEForceField())
     print universe.energyTerms()['electrostatic']


This approach works for all force fields that are implemented as  
subclasses of MMForceField. This includes the SPCEForceField used  
here, but also all variants of Amber (just replace spce_charge by  
amber_charge).

In these force fields, the charges are defined at the level of the  
molecule or, for chain molecules, at the residue level. The  
definition takes the form of a dictionary. There is only one  
dictionary for every molecule/group *definition*, whose use is shared  
by all the molecules/group. If you want to change the charges in only  
one water molecule, you must first make a copy, as shown in the line  
before the loop.

The other tricky aspect is that MMTK caches force field parameters  
for a given system for efficiency reasons. If you want to change a  
force field parameters, you have to delete the cached parameter set.  
The easiest way to do that is to redefine the force field for the  
universe, as is done inside the loop.

Konrad.
--
---------------------------------------------------------------------
Konrad Hinsen
Centre de Biophysique Moléculaire, CNRS Orléans
Synchrotron Soleil - Division Expériences
Saint Aubin - BP 48
91192 Gif sur Yvette Cedex, France
Tel. +33-1 69 35 97 15
E-Mail: hinsen at cnrs-orleans.fr
---------------------------------------------------------------------