[MMTK] Atom ordering and some other problems

Konrad Hinsen hinsen at cnrs-orleans.fr
Mon Aug 6 06:49:31 UTC 2007

On 16.07.2007, at 18:25, Jaroslaw Kalinowski wrote:

> programs to  calculate Poisson-Boltzman electrostatics), second it is
> sometimes easier to use plain arrays rather than objects. For example
> ParticleProperties are fine for representing properties defined for  
> each
> atom in a universe, but how could one represent data for selected  
> objects
> in the universe, or how could one represent two-atom properties  
> like bond
> orders defined for each pair? Arrays seem to be a choice, but for  
> that one

ParticleProperties objects behave like a dictionary, so I'd say the  
first choice for subsets of the whole system would be a dictionary  
mapping atom objects to the associated property value. For atom  
pairs, I'd use a tuple of atoms as the key.

Of course, there is nothing wrong with using arrays and atom indices,  
as long as you keep in mind that every time your program is run, the  
atom indices could be assigned differently in general. Note also that  
atom indices are not defined before the first call to  
universe.configuration(); in case of doubt, call it explicitly.

> 1. Ordering given by u.atomList() is not the same as ordering of
> configuration array, is it?


> 2. If not than how one can iterate over atoms in the same order as  
> is used
> for arrays?

sorted_atoms = copy(universe.atomList())
sorted_atoms.sort(lambda a, b: cmp(a.index, b.index))
for atom in sorted_atoms:
     print atom

> 3. Is it possible to get atom number in pdb file for peptide chains  
> after
> they have been created using PDBPeptideChain.createPeptideChain?

No; in fact, there is no guarantee that every atom even has a PDB  
serial number, e.g. if the MMTK model has hydrogens but the PDB file  

There is one particular situation in which it is possible to access  
the PDB atom numbers: if

1) you use PDBConfiguration.createAll()
2) the PDB file contains all hydrogen atoms of the MMTK model
3) you don't add anything else to the universe

then the internal atom indices are in the same order as the atoms in  
the PDB file. This feature is used for interfacing with other  
simulation programs, e.g. when reading DCD trajectory files.

> 4. Does pickling and unpickling preserve atomList ordering and/or  
> array
> ordering?

Both. Note however that nothing guarantees that atomList ordering  
remains constant even between two successive calls to atomList(). It  
does in the current implementation, but it might not do so in the  

> 5. Does using copy.deepcopy do?

If you copy the whole universe, it does in practice although no  
written rule guarantees that it does. It might not do so in the future.

> 6. Does saving using PDBOutputFile and than reading using  
> PDBConfiguration do?


> 7. And can one know in advance what would be atom ordering in the pdb
> written?


> 8. What is the atom ordering in NetCDF trajectories?

It is given by the system description in the netCDF variable  
"description". Every atom object in that description contains the  
index of that atom in the netCDF arrays.

> 9. Can one assume that each residue has a continuous slice in
> configuration array?

No, except in the special situation described above.

> 10. If so, than is the atom ordering in such slice identical for each
> residue of a given type?

Only of the PDB file has identical atom order for all residues of a  
given type.

> Some more problems.
> Is it possible to store additional data in a trajectory? I have a QM
> forcefield that produces atomic charges, bond orders etc. Is it  
> possible
> to have my function called to provide data for each trajectory frame?

The integrators and minimizers in MMTK don't do this at the moment.  
It should be possible to add such a a feature by writing a suitable  
action object, but I'd better don't promise this without some  

If you write trajectories step by step from Python (using the  
SnapshotGenerator()), you can add any data you wish.

> I would like to try doing QM/MM calculations. It seems that it is not
> possible to apply Amber forcefield to the part of a system, so I was

You can calculate energies and forces for a subset of the system.  
However, there is currently no way to tell the integrators to use  
this feature. If all you need is a simple velocity verlet algorithm,  
you could just rewrite it in a few lines of Python (it's fast enough)  
and calculate energies for subsets.

The "proper" way to do QM/MM in MMTK is to write a new force field.  
That force field would exclude the QM subsystem from the Amber force  
field term list and add the QM energies and forces instead.

> thinking about cloning MM part and adding QM part to a copy. Than I  
> could
> just write a forcefield that would just copy energy and gradients  
> from MM
> universe and add my QM terms. Do you think it is a good approach?

That sounds doable, but you have to make sure that this really gives  
you the QM/MM mixture that you want. The forces on the MM subsystem  
would include interactions with the QM subsystem, for example.

> It is not possible to change protonation state after PeptideChain  
> has been
> created, is it?

No. As far as MMTK is concerned, different protonation states for a  
residue are different residues.

It is possible with moderate effort to replace a residue in a chain  
by another residue, and thus change protonation states. You could  
look at the routine that creates S-S bridges as an example. However,  
any operation that changes the number of atoms invalidates the atom  
indices. They will be reassigned when needed, but this takes time, as  
does the subsequent recalculation of the force field terms etc.  
Changing the system is thus an expensive operation.

If you need changing protonation states in the course of a  
simulation, it is easier and more efficient to start with the fully  
protonated state and modify the force field terms to make unwanted  
hydrogens disappear from all interactions.

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://starship.python.net/pipermail/mmtk/attachments/20070806/febb8fc6/attachment.htm 

More information about the mmtk mailing list