Re: [AMBER] PME with cutoff = 0

From: FyD <fyd.q4md-forcefieldtools.org>
Date: Wed, 26 Feb 2014 07:19:35 +0100

Dear Bob,

> I think Jason's reply to this is the direction you need to look. I
> saw your email and was mostly trying to explain why you could
> misunderstand and think pmemd ran a 0 cutoff for you when sander
> wouldn't. The further story here - the code initializes the cut to
> 0.0 to check and see if you set a value for cut in the namelist, and
> assumes that after you read the namelist, if the cut value is still
> 0.0, then the defaults should be set for pme electrostatics cut and
> vdw cut instead (and it is even more complicated really, in that
> there are other variables specific for electrostatics vs. vdw
> cutoffs you could have specified). We could use some weird negative
> value for cut I suppose, but ultimately someone can surprise you
> with input, and the only way to be sure that they provided input for
> a value, at least as far as I know, is to read the namelist twice
> against different initial values for a variable. Anyway, I defer to
> Jason regarding the solution to your actual problem here
> - my intent was merely to try to be sure you understood what pmemd
> had done. By the way, there should be clear indications in the
> mdout output as to the cutoff actually used

Yes...

[fyd.master0 JOB8-pmemd]$ grep -i cut md-cstpres2b.out
  cst pressure MD, nstlim*dt/1000 nsec, PME with cutoff
   cut = 0,
      dielc = 1.00000, cut = 8.00000, intdiel = 1.00000
      Cutoff= 8.000 Tol =0.100E-04
  TESTING RELATIVE ERROR over r ranging from 0.0 to cutoff

One clearly sees, that I requested cut = 0 and Cutoff= 8.000 was used...

- in the "ewald parameters"
> subsection (and of course in an ideal world, you would get the same
> result using an infinitesimal cut near 0 for the direct space calc
> because the reciprocal space calc would compensate - but in the real
> world, if you start messing with the actual pme run parameters, the
> accuracy of your results varies; messing with the parameters is not
> turning off pme, so you are always attempting a full pbc system
> electrostatics evaluation).

Thank you for your detailed explanation.

regards, Francois


> ________________________________________
> From: FyD [fyd.q4md-forcefieldtools.org]
> Sent: Monday, February 24, 2014 1:31 AM
> To: AMBER Mailing List
> Subject: Re: [AMBER] PME with cutoff = 0
>
> 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...
>
> Does it means this is not possible to carefully study/check solvent
> boxes with Amber nowadays?
>
> thank you,
> regards, Francois
>
>
> Dear Ross,
>
>> Both sander and pmemd almost certainly make assumptions about the cutoff.
>> Since the days of AMBER 4.1 there have been huge changes in the way the
>> pairlist is built etc.
>
> I can understand that ;-)
>
>> For example the code use to just rebuild the
>> pairlist at set intervals. These days it is smart and only rebuilds it
>> when an atom has moved more than the distance skinnb such that the
>> pairlist cutoff is really cut+skinnb. Thus reducing cut to zero probably
>> messes up logic in this part of the code. There are likely other issues as
>> well, like divisions by zero inside the code that determines the Ewald
>> coefficient etc.
>
> My understanding is that division by zero should always be checked.
>
>> The fact one code crashes and the other doesn't is likely
>> just dependent on where they hit the first problem. From what you show it
>> looks like pmemd just hangs while sander crashes with an error - both
>> could be considered a crash.
>
> pmemd does not hang; see above.
>
>> Anyway, to cut a long story short neither code officially supports a value
>> of cut=0 for PME runs and it is something that isn't tested. It may be
>> possible to run the simulation in GB in sander by setting igb=6 - this
>> will run a gas phase simulation. I doubt there is anyway to run a periodic
>> simulation with a cutoff of 0. If you really want to then your options are
>> likely see if you can find a PRE-PME copy of sander (such as AMBER 4.1) or
>> (sander-classic from AMBER 6 perhaps?) or crack open the debugger and see
>> if you can trace through the use of the cutoff value and find out where
>> the problems are occurring. This may just be a case of going down the
>> rabbit hole though.
>
> sander-classic from AMBER 6: do you think this is possible to get
> sander-classic from AMBER 6 compiled using 4.4 GNU compiler?
>
> Does it means this is not possible to carefully study/check solvent
> boxes with Amber nowadays?
> See http://pubs.acs.org/doi/full/10.1021/jp9717655
>
>> Sorry I can't offer much more help than that. - Note you could also try
>> cut=0.001 and see if that works. There is almost certainly a 1.0d0/cut^2
>> somewhere in the code and this might work around that and still give you
>> what you want from the simulation.
>
> I am testing...
>
> thank you,
> regards, Francois
>
>
>> On 2/23/14, 9:19 AM, "FyD" <fyd.q4md-forcefieldtools.org> wrote:
>>
>>> Dear All,
>>>
>>> I ran amber10/sander.MPI with PME + cutoff = 0 -> sander crashes
>>> amber12/sander.MPI with PME + cutoff = 0 -> sander crashes
>>> (compiled with intel13 + mkl)
>>> amber12/pmemd.MPI with PME + cutoff = 0 -> pmemd does not crash...
>>>
>>> Printed output: printing stops in section 4:
>>>
>>> --------------------------------------------------------------------------
>>> 4. RESULTS
>>> --------------------------------------------------------------------------
>>>
>>> | # of SOLUTE degrees of freedom (RNDFP): 24936.
>>> | # of SOLVENT degrees of freedom (RNDFS): 0.
>>> | NDFMIN = 24933. NUM_NOSHAKE = 0 CORRECTED RNDFP = 24933.
>>> | TOTAL # of degrees of freedom (RNDF) = 24933.
>>>
>>> -> nothing else is printed...
>>>
>>> In Fox & Kollman J.Phys.Chem.B 1998, 102, 8070, MD simulation with
>>> cutoff = 0 is reported & Sander from Amber 4.1 was used.
>>>
>>> Why does sander crash, while pmemd does not crash in this case?
>>> (with the same input & cut = 12, sander does not crash)
>>>
>>> thank you, regards,
>>> Francois
>>>
>>>
>>>
>>> # Cst pressure simulation with cutoff = 0 Angs; 300 K short
>>> $PARAL -np $NPROC $EXE -O -i md-cstpres0.in \
>>> -o md-cstpres0b.out \
>>> -p $mol.parm7 \
>>> -c md-cstvol.rst7 \
>>> -r md-cstpres0b.rst7 \
>>> -x md-cstpres0b.mdcrd
>>>
>>> my test input below:
>>>
>>> cst pressure MD, nstlim*dt psec, PME with cutoff
>>> &cntrl
>>> nmropt = 0, ioutfm = 1, iwrap = 0,
>>> ntx = 1, irest = 0, ntrx = 1, ntxo = 1,
>>> ntpr = 100, ntwr = 1000, ntwx = 500,
>>> ntave = 0, ntwv = 0, ntwe = 0,
>>> ntf = 2, ntb = 2,
>>> ntc = 2, tol = 0.00001,
>>> dielc = 1, ipol = 0,
>>> cut = 0, ig = 71277, comp = 52.5,
>>> ibelly = 0, ntr = 0, igb = 0,
>>> imin = 0, nstlim = 5000,
>>> nscm = 1000, t = 0.0, dt = 0.002,
>>> temp0 = 300, tempi = 300,
>>> ntt = 1, tautp = 0.2, vlimit = 20,
>>> ntp = 1, pres0 = 1.0, taup = 0.2,
>>> &end
>



_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Feb 25 2014 - 22:30:02 PST
Custom Search