On Mon, 2014-02-24 at 10:31 +0100, FyD wrote:
> Dear Robert,
>
> > The cutoff value of "0.d0" in pmemd for a pme run is equivalent to
> > requesting that the default values be used - probably 8.0 angstrom,
> > if memory serves, assuming igb 0 for a pme simulation.
>
> Yes - this corresponds to what I get...
>
> I compare MD simulations (with PME) with cut = 12 vs cut = 0
> the obtained results are very similar.
>
> This means 'cut = O' does not do what I want; I wonder if crashing is
> not better than generating data in this case...
>
> > PME does strange things at small-ish cutoff, but does not allow 0
> > (in the sense that it says "oh, you wanted 8"), and I presume later
> > catches and disallows cut < 0. In my experience, reducing the
> > direct space cut below roughly 7.0 angstrom causes unacceptable
> > error limits unless you increase grid density for reciprocal space
> > (there is basically a balancing act between reciprocal and direct
> > space in terms of what is necessary to keep error acceptable). And
> > if you think long cutoffs are expensive, play around a bit with
> > greatly increasing nfft1,2,3 in a large system :-}. So all I am
> > really answering here is why pmemd ran - it basically ignored you...
>
> ok
>
> > I would have to read the paper done with amber 4.1 to have a clue
> > what they were really up to,
>
> See http://pubs.acs.org/doi/full/10.1021/jp9717655 - It is written:
>
> "The heat of vaporization is computed from the average intermolecular
> interaction energy Eint via dHvap = -Eint + RT
> For molecules that do not have an internal nonbonded interaction
> beyond 1?4 (e.g., the chloromethanes), Eint is calculated
> straightforwardly as the sum of EELEC and ENONB in the SANDER output,
> divided by the number of molecules in the system. For all other
> molecules, EELEC and ENONB also contain contributions from
> intramolecular nonbonded interactions. To correct for these, we
> performed a short simulation with the nonbonded cutoff radius set to
> zero. Due to the residue-based cutoff in AMBER, EELEC and ENONB now
> accumulate only the intraresidue nonbonded interactions that can be
> subtracted from the total to yield the intermolecular portion needed
> for the calculation of dHvap."
>
> > but pmemd is designed to always use either
> > pme or generalized Born, so resurrecting an old copy of sander may
> > be your best bet..
>
> I use PME & would like to continue to use PME.
> Resurrecting an old copy of sander? uff...
I would actually suggest that you really _don't_ want to use PME for
this particular calculation. What you're effectively trying to do is
eliminate the intra-molecular nonbonded interactions so that you can
compute the inter-molecular nonbonded interactions exclusively. It
seems that the strategy is to compute the full electrostatics (using
PME) and then subtract out the intramolecular nonbonded energies (using
no PME and a group-based cutoff of 0 to ensure that atoms interact only
with other atoms in the same residue).
This is a trick that takes advantage of a key assumption: each solvent
molecule is made up of a single residue only.
> Does it means this is not possible to carefully study/check solvent
> boxes with Amber nowadays?
I've never tried pure group-based cutoffs with no Ewald sum. It may not
work anymore since it's a somewhat unusual simulation (nor is it
particularly general since it relies on the above assumption). If
group-based cutoffs still work and the assumption above holds for your
system, that's the easiest solution. There are alternatives, though,
that depend on the nature of your 'molecule'.
1) Does the molecule have any atoms that are no farther than 2 bonds
away? If not (like with water and chloromethane), then there are no
intramolecular nonbonded interactions and this step is unnecessary.
2) Does the molecule have any atoms that are no farther than 3 bonds
away? If not (like with ethane and ethanol), then you can simply ignore
the 1-4 VDW and 1-4 EEL contributions (again, making this step
unnecessary).
3) The most tedious scenario if there are intramolecular energies you
need to get rid of. You can, using ParmEd, zero out the charges and van
der Waals parameters for every molecule *except* the molecule whose
intramolecular NB interactions you want and simply compute a
(non-periodic) energy in which you know the EEL and VDW contributions
are strictly intramolecular. You can then do this for every molecule in
your system and sum them up to get the total intramolecular nonbonded
energy (or, alternatively, if your molecules are rigid due to
constraints you can simply multiply by the number of molecules you
have).
If you are comfortable with Python programming, you can do step 3 in a
single Python script using the ParmEd API like so:
from chemistry.amber.readparm import AmberParm
from ParmedTools.ParmedActions import change, addljtype
for i in range(num_molecules):
parm = AmberParm('path/to/prmtop')
change(parm, 'CHARGE', '!:%d' % i, 0.0)
addljtype(parm, '!:%d' % i, radius=0.0, epsilon=0.0)
parm.writeParm('tmp.parm7')
# System call to compute energy
Alternatively you can write a quick NAB program to do something
equivalent. In both cases, help can be found in the AmberTools manual.
HTH,
Jason
--
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Mon Feb 24 2014 - 06:00:06 PST
