[MMTK] how to calculate density of state for a specific atom

Konrad Hinsen hinsen at cnrs-orleans.fr
Thu Mar 15 16:37:48 UTC 2007

On Mar 15, 2007, at 16:59, Nitin wrote:

>    I have a trajectory file (.nc file) for a collection of atoms  
> e.g. 20 Hydrogen, 10 Oxygen etc. Using xMoldyn one can calculate  
> the density of state (d.o.s.) for the whole set i.e. 20 H + 10 O ,  
> or for a group of atoms e.g. 20H or 10 O only; but how can we  
> select a particular Hydrogen let say the 7th Hydrogen and calculate  
> it's d.o.s. ?

You can make arbitrary atom selections using Python code. The pMoldyn  
input file shown below is an example.


#** Example for Python-coded atom selection in a protein.
#** This file was initially created by saving the pMoldyn
#** input file for a VACF calculation in xMoldyn. The
#** modifications are explained in the comments that
#** start with **.

# pMoldyn data
from MMTK import *
title = 'Velocity Autocorrelation Function (from coordinates)'
trajectory = ['lysozyme_vacuum_nvt_1ns.nc']
log_file = 'vacf-xyz_lysozyme_vacuum_nvt_1ns.log'
time_info = (0, 25000, 1)
#** Comment out the line with the xMoldyn atom selection
#atoms = {'Protein.0': ['*']}
projection_vector = None
weights = 'mass'
differentiation = 'fast'
units_length = Units.nm
time_steps = 2000
output_files = {'vacf': 'VACF_lysozyme_vacuum_nvt_1ns.plot'}

#** Write a function "atoms_code" that takes the trajectory as its
#** argument and that returns an MMTK object whose atoms represent
#** the selection.

#** The example selection function given above selects the sidechains
#** of all alanine residues in all peptide chains of all proteins
#** in the simulated system. All the functions used here, as well as
#** other useful functions for atom selection, are described in the
#** MMTK user manual.

#** To run this calculation, type
#**    pMoldyn --vacf-xyz --input atom_selection_example.inp
#** Watch nMOLDYN print the number of selected atoms (48 for this case),
#** this allows you to check if everything works well.

def atoms_code(trajectory):
     import MMTK
     from MMTK.Proteins import Protein
     selection = MMTK.Collection()
     universe = trajectory.universe
     proteins = universe.objectList(Protein)
     for protein in proteins:
         for chain in protein:
             for residue in chain:
                 # Residue names consist of the three-letter code
                 # followed by the residue number, e.g. Ala11.
                 if residue.name[:3] == 'Ala':
     return selection

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

More information about the mmtk mailing list