[MMTK] Atom ordering and some other problems
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
> atom in a universe, but how could one represent data for selected
> 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:
> 3. Is it possible to get atom number in pdb file for peptide chains
> 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
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
> 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
> 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
> 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
> 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.
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...
More information about the mmtk