AMBER Archive (2008)

Subject: Re: AMBER: pmemd iwrap trouble

From: Robert Duke (
Date: Tue Aug 26 2008 - 22:17:29 CDT

Okay gang, I have reviewed the pmemd vs. sander code as it applies to
wrapping in general, and wrapping the truncated octahedron in particular.
The code, as I said, is basically derived from code that existed in amber 6
sander. Since then, there has been 1 bugfix, applied to both sander and
pmemd. In addition, in sander 6 there was conditionally compiled code that
opted to use the fortran 90 sign() intrinsic instead of a conditional (foo
.gt. bar) to determine whether or not to add or subtract a value to
coordinates in wrapping molecules for the truncated octahedron. Since amber
6 the code that uses sign() has become the default for sander but not pmemd.
This is probably a bit faster, but insignificantly so, given the frequency
of wrapping. The other "feature" you get from using sign() is slightly
enhanced precision, whereby this function distinguishes between +0.d0
and -0.d0. Using the conditional as pmemd does, all values that evaluate
to 0.d0 to within the precision of storage wrap in one direction; for sander
using sign(), half the 0.d0 values go one way and the other half go the
other way. This is basically an insignificant difference in a rare event,
but would have potential to cause sander and pmemd to do a slightly
different wrap. I don't think you can do something equivalent in C, and the
f90 spec states that sign() is the only place where this distinction is made
in f90. So for all intents and purposes, sander and pmemd do the same
thing, and it is correct... I do not know what the correct machinations to
get what you want in ptraj are, but Tom C. is obviously the authority, so I
would do what he says. I tend to just do single step wraps in pmemd to
prevent restart overflow, and I never worry about imaging truncated
octahedra (I rarely use them). There was an indication from Tom C. that a
recent release of ptraj was broken, so that is potentially a source of some
grief here. I think we have beat this to death by now, and I hope everyone
has what they need (if not, keep squawking); I am reasonably confident that
pmemd is doing the "right" thing.
Regards - Bob Duke

----- Original Message -----
From: "Robert Duke" <>
To: <>
Sent: Tuesday, August 26, 2008 3:49 PM
Subject: Re: AMBER: pmemd iwrap trouble

> You might try wrapping in 1 step of 0.01 fsec (ie., dt = .00001). That
> way the before and after delta in all atom positions should be reduced
> 100-fold, and you should see a smaller energy delta. You would expect
> some difference, as you are really minimizing two structures that differ
> by more than just the wrapping operation - they also differ by the
> distance moved by all atoms in the dynamics step. Also, using the 32
> angstom water buffer sounds pretty excessive. Is your protein complex
> perhaps not a good candidate for fitting in a truncated octahedron? (so a
> globular protein approximating a sphere is the best candidate; elongated
> structures may fit better in an orthogonal unit cell). Also note my
> earlier comments about no_intermolecular_bonds.
> Regards - Bob
> ----- Original Message -----
> From: "Lars SkjŠrven" <>
> To: <>
> Sent: Tuesday, August 26, 2008 2:36 PM
> Subject: Re: AMBER: pmemd iwrap trouble
>> Hi Tom,
>> I performed a 1 step energy minimization of the restart files before
>> and after the wrapping in pmemd10. The energies are indeed similar,
>> but not identical: -2.1748E+06 and -2.1750E+06.
>> When I do "image center familiar" in ptraj of the restart files before
>> and after pmemd wrapping, and then an rmsd in ptraj of the two new
>> restart files, it gives an rmsd of 1.6┼. visualizing in vmd shows that
>> most of the protein aligns, but two subunits are translated slightly
>> with respect to the 10 other subunits.
>> Ptraj scripts used:
>> ---
>> trajin b4_wrapping.rst
>> trajin after_wrapping.rst
>> center :1-6456
>> image center familiar
>> trajout reimaged.rst restart
>> ----
>> and then,
>> ---
>> trajin reimaged.rst.1
>> trajin reimaged.rst.2
>> rms first mass out rmsd.dat @N,CA,C
>> ----
>> Yes, I am still concerned about the validity of my simulation.. :-)
>> Do I have any other choice to tackle this overflow issue? Can I use
>> ptraj to reimage the water back in the box (which works fine), and
>> make a new restart file (which I assume do not contain velocities,
>> right?), and copy the velocities from the old restart file (which are
>> about to overflow)? then, use this new and edited restart file to
>> continue my simulation?
>> When I built my system I had trouble with "solvateoct 10┼" completely
>> solvating my protein, instead I had to use "solvateoct 32┼"
>> (
>> using only 10┼ left big chunks of my protein out of the waterbox. for
>> any other proteins besides this big monster protein I have worked
>> with, "solvateoct 10┼" works fine. can my problem with wrapping have
>> anything to do with this issue? probably not.. ?
>> Thanks for all help on this issue.
>> Best regards,
>> Lars Skjaerven
>> University of Bergen, Norway
>> On Tue, Aug 26, 2008 at 6:59 PM, Robert Duke <> wrote:
>>> Hi Tom,
>>> Thanks, that was my understanding that the imaging in sander and pmemd
>>> was
>>> the same, and "equivalent" to a truncated octahedron, but that the image
>>> would not look correct in something like vmd. I was hoping to get Tom
>>> D. to
>>> comment on this since he did the code. I don't use truncated
>>> octahedrons
>>> that much myself, but have had 0 problems with other imaging, aside from
>>> problems with not being quite sure what to do in ptraj sometimes, and
>>> also
>>> ptraj understandably can't keep up with box size changes associated with
>>> a
>>> constant pressure simulation (so wrapping a whole trajectory rather than
>>> a
>>> frame really blows up; one might expect that...). Anyway, thanks for
>>> the
>>> input. I am going to confirm that my code is identical to sander in
>>> amber
>>> 10 with regard to wrapping, but that is what I expect. There IS a known
>>> problem with the amber 9 release of pmemd, for which a patch was
>>> released,
>>> so some folks with an unpatched version of 9 may see a few molecules not
>>> wrap correctly (so do the patch...).
>>> Regards - Bob
>>> ----- Original Message ----- From: "Thomas Cheatham III" <>
>>> To: <>
>>> Sent: Tuesday, August 26, 2008 12:48 PM
>>> Subject: Re: AMBER: pmemd iwrap trouble
>>>>> I get exactly the same results from pmemd9, pmemd10, and sander10 in
>>>>> regard the wrapping (iwrap=1); 2 of my subunits are displaced
>>>>> (regardless of ioutfm=1 or 0). However, by using the following ptraj
>>>>> script (instead of iwrap) it results in a truncated octahedral fitted
>>>>> nicely around the protein complex (as expected):
>>>>> center :1-5332
>>>>> image center familiar
>>>>> trajout reimaged.rst restart
>>>>> However, I see no trouble when I use solvatebox instead of solvateoct
>>>>> (for both sander and pmemd). Unfortunately I cant go back using
>>>>> solvatebox at this time, after 2 months of simulations..
>>>> I've been following this thread and think that it is just a
>>>> misunderstanding and that nothing is wrong with the imaging, it just
>>>> doesn't look like you would expect. sander/pmemd image in triclinic
>>>> space
>>>> so for a protein this looks like a slanted box with the protein
>>>> potentially sticking out...
>>>> __PPP__
>>>> / PPP /
>>>> /___PPP/
>>>> ...rather than the more "familiar" truncated octahedron shape. To
>>>> verify
>>>> this, look at the ptraj imaged trajectory without the familiar keyword;
>>>> this will look like what you are seeing from sander/pmemd and they are
>>>> equivalent, albeit different ways of looking at things.
>>>> There is also a way in pmemd to generate images in alternative unit
>>>> cells,
>>>> image xoffset 1.0 yoffset 0.0 zoffset 0.0
>>>> which you can do for each direction and then "see" the whole packing in
>>>> the unit cell.
>>>> If you are still concerned, you could also try a single point energy
>>>> minimization on both the "familiar" and alternative restrt file to
>>>> verify
>>>> the energies are similar.
>>>> imaging in ptraj is not broken (now) to the best of my knowledge; there
>>>> was a brief problem with AmberTools 1.0 that lead to errors in the box
>>>> size, but this was working in earlier versions and is working in
>>>> current
>>>> versions.
>>>> --tom
>>>> -----------------------------------------------------------------------
>>>> The AMBER Mail Reflector
>>>> To post, send mail to
>>>> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
>>>> to
>>> -----------------------------------------------------------------------
>>> The AMBER Mail Reflector
>>> To post, send mail to
>>> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
>>> to
>> -----------------------------------------------------------------------
>> The AMBER Mail Reflector
>> To post, send mail to
>> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
>> to
> -----------------------------------------------------------------------
> The AMBER Mail Reflector
> To post, send mail to
> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
> to

The AMBER Mail Reflector
To post, send mail to
To unsubscribe, send "unsubscribe amber" (in the *body* of the email)