AMBER Archive (2008)

Subject: Re: AMBER: pmemd iwrap trouble

From: Carlos Simmerling (carlos.simmerling_at_gmail.com)
Date: Wed Aug 27 2008 - 08:26:54 CDT

Hi Bob,
thanks for looking into this so thoroughly!
carlos

On Tue, Aug 26, 2008 at 11:17 PM, Robert Duke <rduke_at_email.unc.edu> wrote:
> 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" <rduke_at_email.unc.edu>
> To: <amber_at_scripps.edu>
> 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"
>> <lars.skjarven_at_biomed.uib.no>
>> To: <amber_at_scripps.edu>
>> 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Å"
>>> (http://structbio.vanderbilt.edu/archives/amber-archive/2008/1109.php).
>>> 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 <rduke_at_email.unc.edu> 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" <tec3_at_utah.edu>
>>>> To: <amber_at_scripps.edu>
>>>> 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 amber_at_scripps.edu
>>>>> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
>>>>> to majordomo_at_scripps.edu
>>>>>
>>>>
>>>> -----------------------------------------------------------------------
>>>> The AMBER Mail Reflector
>>>> To post, send mail to amber_at_scripps.edu
>>>> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
>>>> to majordomo_at_scripps.edu
>>>>
>>> -----------------------------------------------------------------------
>>> The AMBER Mail Reflector
>>> To post, send mail to amber_at_scripps.edu
>>> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
>>> to majordomo_at_scripps.edu
>>>
>>
>> -----------------------------------------------------------------------
>> The AMBER Mail Reflector
>> To post, send mail to amber_at_scripps.edu
>> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
>> to majordomo_at_scripps.edu
>>
>
> -----------------------------------------------------------------------
> The AMBER Mail Reflector
> To post, send mail to amber_at_scripps.edu
> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
> to majordomo_at_scripps.edu
> 

-- 
===================================================================
Carlos L. Simmerling, Ph.D.
Associate Professor Phone: (631) 632-1336
Center for Structural Biology Fax: (631) 632-1555
CMM Bldg, Room G80
Stony Brook University E-mail: carlos.simmerling_at_gmail.com
Stony Brook, NY 11794-5115 Web: http://comp.chem.sunysb.edu
===================================================================
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
      to majordomo_at_scripps.edu