Hi Eng,
Also please note a subtle difference between pmemd and sander in terms of how each defines molecules by default. Sander just takes what is in the prmtop, in terms of molecule definitions, as the truth, and uses this information in things like the virial (pressure) calculations. PMEMD, on the other hand, actually does a check for covalent bonds in any putative "molecules", and at least in the past, it was common for a "molecule" in a prmtop to include a reduced disulfide (ie., thiols). Thus, if a disulfide is not joined between two subunits, pmemd may well see it as two separate molecules. The reasons for this different behavior in pmemd have to do mostly with complexities introduced in force fields that use extra points which typically are positioned in a frame defined by other atoms. The pmemd default behavior MAY be overridden to match sander behavior exactly by use of the pmemd-specific no_intermolecular_bonds switch. This is generally not something that would have much effect at all on pressure unl
ess such a molecule was numerous in the system (like solvent, which seems sort of unlikely). On the other hand, it is conceivable that there could be a bit of strain introduced across a disulfide by assuming two separate molecules in an ntp calculation, depending on exact geometries, as pressure rescaling is done based on molecular center-of-mass (my suspicion - very small impact, but I have never checked). The discussion of the relevant parameter is in the Amber manuals (starting with Amber 10, when the option was introduced to support a high performance implementation of force fields with extra points). The manual entry is as follows:
no_intermolecular_bonds In &cntrl. New variable controlling molecule definition. If 1, any
molecules (ie., molecules as defined by the prmtop) joined by a covalent bond are
fused to form a single molecule for purposes of pressure and virial-related operations;
if 0 then the old behaviour (use prmtop molecule definitions) pertains.
The default is 1; a value of 0 is not supported with forcefields using extra points.
This option was necessitated in order to efficiently parallelize model systems with
extra points. This redefinition of molecules actually allows for a more correct
treatment of molecules during pressure adjustments and should produce better results
with less strain on covalent bonds joining prmtop-defined molecules, but if
the default value is used for a NTP simulation, results will differ slightly relative
to sander if any intermolecular bonding was applied in forming the prmtop (eg.,
a cyx-cyx bridge was added between two peptides that originated in a PDB file,
with each peptide having its own "TER" card). If consistency with sander is more
important to you, and you are not using extra points, then you may want to set
no_intermolecular_bonds to 0.
Regards - Bob Duke
________________________________________
From: Eng H Yap [eyap.einstein.yu.edu]
Sent: Wednesday, December 05, 2012 8:56 AM
To: amber.ambermd.org
Subject: [AMBER] FW: EGB discrepancy between sander and pmemd
Hi Jason,
Wow, I am glad you found my mistake early, before I produce more MD runs!
It works now with CYS changed to CYX following instructions in AmberTools manual. EGB in sander and pmemd matches.
I will *always* remember to check for disulfide bonds from now on!
Thank you so much,
Enghui
From: Jason Swails [jason.swails.gmail.com]
Sent: Tuesday, December 04, 2012 6:40 PM
To: AMBER Mailing List
Subject: Re: [AMBER] EGB discrepancy between sander and pmemd
You do have strongly overlapping atoms, which is a (fortunate!!) byproduct
of a different mistake. You forgot to add two CYS-CYS disulfide bonds when
you built your topology file. By leaving the residue named CYS (instead of
CYX), tleap automatically added hydrogens on the S atoms that were supposed
to be bonded to each other. These two hydrogens overlap strongly, and it's
these interactions that blew up your system.
This is something important to think about every time you go to create a
topology file of a protein. You have to make sure before-hand that not
only are the CYS residues involved in a disulfide bond renamed to CYX, but
that the necessary bond is formed. Loading a Amber prmtop into VMD will
display these bonds and help you diagnose this problem. I have been burned
by this enough myself (and had to throw away large amounts of useless MD),
which is why I'm so vigilant about it now.
Once you rebuild your topology file with the necessary disulfide bonds,
your errors should go away.
Good luck,
Jason
On Tue, Dec 4, 2012 at 6:29 PM, Jason Swails <jason.swails.gmail.com> wrote:
>
>
> On Tue, Dec 4, 2012 at 6:02 PM, case <case.biomaps.rutgers.edu> wrote:
>
>> On Tue, Dec 04, 2012, Eng H Yap wrote:
>> >
>> > I am still trying to set igb=0 as suggested by Dave, but it seems then
>> > I'll need ntb=1, and I'll have to solvate the box in leap. Is that
>> > right?
>>
>> Sorry: easier path is to set igb=6, which turns off GB, but doesn't
>> require
>> (or allow) ntb=1.
>>
>
> I'll see if your results are reproducible with my version...
>
> pmemd actually prohibits all vacuum calculations -- there's a check
> against igb==0 with ntb==0 and igb==6. igb==0 with ntb==0 and igb==6 are
> sander-only options.
>
> All the best,
> Jason
>
> --
> Jason M. Swails
> Quantum Theory Project,
> University of Florida
> Ph.D. Candidate
> 352-392-4032
>
--
Jason M. Swails
Quantum Theory Project,
University of Florida
Ph.D. Candidate
352-392-4032
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Dec 05 2012 - 12:30:02 PST
