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

From: NAVEEN P MICHAUD AGRAWAL (nmichaud_at_jhu.edu)
Date: Fri Jul 29 2005 - 23:30:23 CDT

I had this same exact problem, and it was apparent when I viewed the merged trajectories in VMD because I
used CHARMM to merge them. This issue was raised in the NAMD wiki: http://www.ks.uiuc.edu/Research/
namd/wiki/index.cgi?NamdAndCHARMM. I hacked the dcdplugin.c file from the newest version of catdcd to
take a trajectory and write out a new one with the correct cell information. Tthis was back in January so I don't
remember exactly what i did, but I've included a tar of the source directory (which includes the files i copied
from catdcd...hopefully I
am allowed to redistribute it). You can see my changes by doing a diff between dcdcellfix.c and dcdplugin.c
(from the catdcd
program). To compile it just type

gcc -I./ -DTEST_DCDPLUGIN dcdcellfix.c -o dcdcellfix -lm

and then

./dcdcellfix [trajectory filename] [output trajectory]

If you don't specify an output file it will default to test.trj. Hopefully this will be helpful to anyone else having this
problem.

Naveen

----- Original Message -----
From: Mauricio Carrillo Tripp <trippm_at_gmail.com>
Date: Friday, July 29, 2005 5:36 pm
Subject: namd-l: Report on a trajectory incompatibility between NAMD and CHARMM: Unfolding and RDFs

> 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<http://chem.acad.wabash.edu/%7Etrippm/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<http://chem.acad.wabash.edu/%7Etrippm/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<http://chem.acad.wabash.edu/%7Etrippm/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
> simulationset 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
> fileopen 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 - 05:18:55 CST