Report on a trajectory incompatibility between NAMD and CHARMM: Unfolding and RDFs

From: Mauricio Carrillo Tripp (trippm_at_gmail.com)
Date: Fri Jul 29 2005 - 16:36:07 CDT

Report on a trajectory incompatibility between NAMD and CHARMM.

I describe a couple of problems we've had when applying Charmm's analysis
tool set to a
trajectory obtained from a NAMD simulation (using Charmm's FF) and the way
we came
around them. I apologize for the length of this posting but I hope that it
will
help other people who, either now or in the future, try to follow the same
analysis procedure.
What I describe might even be a bug in NAMD, and hopefully the developers
will fix it in
future versions. We currently use NAMD 2.5 for Linux-i686-TCP.
We usually build the system using Charmm, creating the psf file both in
Charmm and
XPlore formats. We use the latter to run the simulation on NAMD, obtaining a
DCD trajectory
file. Finally, we go back to Charmm and analyze the DCD.
The systems we work with are lipid membranes on different conditions and
different sizes;
this might be relevant for the problems described below since they are
anisotropic.

Unfolding:

Some time ago I reported to this list a problem I had while trying to unfold
a trajectory obtained from a
NAMD simulation using Charmm's FF. The problem was that the coordinates got
messed up
during the process of unfolding the system using the Charmm's tool set. This
was clearly seen when
looking at the trajectory with VMD, the output would look something like
this:
http://chem.acad.wabash.edu/~trippm/tmp/unfold-namd.png>

We tried a little trick (read below), and that fixed the problem:
http://chem.acad.wabash.edu/~trippm/tmp/unfold-charmm.png>

RDFs:

Recently we ran into some confusing results while trying to calculate a
radial distribution
function, also using the Charmm's tool set, on a trajectory file produced by
NAMD. We used
a couple of different commands, different options, different system size and
different image
transformations. All of them should have given the same result. However,
they didn't.
Most confusing was the fact that any of those results were the same to the
result obtained from
the same analysis made to a trajectory file produced by Charmm (which seemed
to be the right answer).
We didn't understand why would this be happening, until we remembered the
'small' issue we had before
when trying to unfold the system. We applied the same 'trick' and everything
worked just perfect!
In the next figure, panel (a), are shown the results of the analysis done on
the trajectory file produced by
NAMD. It can be clearly seen that there is a disagreement between all of the
methods used.
Panel (b) shows all the results of the same methods applied to the 'fixed'
trajectory file:
http://chem.acad.wabash.edu/~trippm/tmp/rdf-comparison.png>
As shown, only when 'fixing' the trajectory file all the methods used to
calculate the RDF give the
same result. Moreover, they are the same as the result obtained when
analysing the trajectory
produced by Charmm.

The 'fix':

There seems to be a difference between the way NAMD writes the cell's info
to the DCD file and
the way Charmm reads it (and writes it). This is why we tried to reset this
info in the trajectory file.
The first step is to read in the trajectory produced by NAMD into Charmm.
Then, we dump only the
coordinates into a temporary file. After defining the cell's lengths and
geometry, we read in the
coordinates from the temporary file and write the final 'fixed' DCD
trajectory file. After doing this,
all the analisys made with Charmm that involve periodic boundaries
operations and the creation of
images should work. This might not be the best or optimized solution, but it
has worked for us so far.

As I mentioned before, this problem might be only apparent for systems that
are not isotropic, but for consistency's
sake the trajectory file should be corrected anyway. At least until this bug
is fixed in the code, shall it be one.

For those who would like to try it for themselves, here's the Charmm script:
**************************************************************************************************************************************************************

! Charmm's configuration file to RESET the Cell's info in the trajectory
file
! to overcome an incompatibility between NAMD and Charmm
! by trippm_at_wabash.edu and fellers_at_wabash.edu
! Wabash College 290705

set OPT 0 ! set this to [0] if just resetting the cell's info, [1] if
unfolding the system

set NS 500000 ! set this to the total number of steps in the simulation
set NF 1000 ! set this to the actual number of frames in the trajectory file

set boxx 49.00 ! set this to the cell's X length
set boxy 49.00 ! set this to the cell's Y length
set boxz 54.00 ! set this to the cell's Z length

! Read in the PSF file
open unit 50 read card name my_system.psf
read psf card unit 50
close unit 50

! Read in the NAMD DCD trajectory file, and write it out to a temp file
open unit 11 read file name my_system-namd.dcd
open unit 1 write file name tmp.dcd

! Calculate the skip value
set SK @NS
divide SK by @NF

merge coor firstu 11 nunit 1 skip @SK output 1 nfile @NF

close unit 1
close unit 11

! Define the cell's info, Charmm's format...
crystal define tetragonal @boxx @boxy @boxz 90.0 90.0 90.0
open unit 1 read card name image.cry
crystal read unit 1 card
close unit 1

! Read in the trajectory temp file, and write it out to a Charmm's formated
file
open unit 11 read file name tmp.dcd
open unit 1 write file name my_system-charmm.dcd

if OPT .eq. 0 then merge coor firstu 11 nunit 1 skip @SK output 1 nfile @NF
if OPT .eq. 1 then merge coor firstu 11 nunit 1 skip @SK output 1 nfile @NF
unfold

stop

************************************************************************************************************************************************************** 

-- 
Mauricio Carrillo Tripp, PhD
Department of Chemistry
Wabash College
trippm_at_wabash.edu
http://trippm.bajacast.com/

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:41:00 CST