Re: [AMBER] "closest" command in cpptraj

From: Daniel Roe <>
Date: Mon, 16 Apr 2012 17:25:54 -0400


On Mon, Apr 16, 2012 at 3:57 PM, Alfredo Quevedo
<> wrote:
> (LIG) with the corresponding 10 WAT molecules. How I am supossed to
> indicate cpptraj to conserve the rest of the receptor and calculate WAT
> distances just based in LIG residue?

Currently it is not possible in ptraj or cpptraj to specify separate
masks for stripping and distance calculation in the closest command,
although this would be a nice feature and may be included in an
upcoming release, maybe a patch.

Currently the only way I can think of to do this is a bit involved but
could probably be made somewhat tractable through scripting. You could
generate PDB files of the ligand with the 10 closest molecules,
generate corresponding PDB files of the receptor, then concatenate
them outside of cpptraj and convert them back into a trajectory. So
for example to generate the ligand/water and receptor PDBs in one

trajin Solvated-trajectory.mdcrd
closest 10 :LIG first closestout archive.dat
outtraj Ligand-10closest.pdb pdb multi
strip :WAT,LIG outprefix strip
outtraj Receptor.pdb pdb multi

This will generate two pdb files per frame, one named
Ligand-10closest.pdb.X and the other named Receptor.pdb.X, where X is
the frame number. Note the unstrip command, which undoes the effect of
the closest command so that the receptor is present for 'strip'. You
can concatenate each of these PDBs together for each frame (removing
the END cards) to generate frames with receptor/ligand/10 waters using
a script like (assuming 10 frames):

for (( frame=1; frame <= 10; frame++)) ; do
  cat Receptor.pdb.$frame Ligand-10closest.pdb.$frame | grep -v "END
" > Output.pdb.$frame

Then you can generate a parm containing the receptor, ligand, and 10
waters either with tleap or you can do it with cpptraj. For the
following example I have a system which the receptor is residues 1 -
268, the ligand is residue 269, and the first 10 water molecules are

parm Solvated.parm7
parmstrip !(:1-279)
parmwrite out closest.Solvated.parm7

Substitute your actual parm name here. For this example a parm
'closest.Solvated.parm7' would be generated containing residues 1 to
279. Adjust the mask accordingly for your system.

Last, you would combine all of the Receptor/ligand/10 water PDBs back
into a trajectory.

echo "parm closest.Solvated.parm7" >
for (( frame=1; frame <= 10; frame++)) ; do
  echo "trajin Output.pdb.$frame" >>
echo "trajout netcdf nobox" >>

cpptraj -i

At this point you should have a topology (closest.Solvated.parm7) and
netcdf trajectory ( containing the receptor, ligand, and
10 closest waters. Phew!

Thanks for the good feature suggestion. Let me know if you have any
more troubles or if I can be of any more help.


AMBER mailing list
Received on Mon Apr 16 2012 - 14:30:02 PDT
Custom Search