On Thu, Jan 22, 2015 at 3:37 PM, Joseph Baker <bakerj.tcnj.edu> wrote:
> Gerald, the molecules are the ones in the amber ionic liquid tutorial, so
> bigger than water. I just tried to set skinnb = 4.0 under &ewald in my
> input script, and am getting basically the same behavior as with the
> default skinnb (total energy vs time stays roughly flat for ~ 1 ns or so,
> and then starts to really shoot up, giving an overall drift of about 6 x
> 10^-4 kT/ns/dof).
>
I'm not sure if this was always true, but in the current versions of the
code skinnb has no effect at all on the results (just in efficiency). [1]
Jason, I was not using SHAKE in the 1 fs timestep runs, however in the 2 fs
> runs I did. I have two runs (that started from different annealing
> protocols, so this is not really comparing apples to apples I suppose) and
> those used the following parameters
>
> dt = 0.002
> ntc=ntf=2
> And default values for tol and dsum_tol
> Results in drift = 1.0 x 10^-2 kT/ns/dof
>
> dt = 0.002
> ntc=ntf=2
> tol = 0.000001
> dsum_tol = 0.000001
> Results in drift = 2.6 x 10^-3 kT/ns/dof
>
I suspect this is largely due to the tolerance changes. You could always
compute the forces for each starting position and compare them, but I
suspect the distribution of force magnitudes will be roughly the same
(which I think is really what matters most to an integrator wrt time-step
dependence of energy drift).
This makes it appear that there is almost an order of magnitude improvement
> in the drift when those tolerance values are changed, however, the second
> case was the result of a longer (and probably better) annealing scheme. If
> you don't think that this is telling, I can run a more controlled test of
> the tolerance parameter by using the same starting configuration in each
> case. Is it worth is to make the tolerance even smaller than 10^-6?
I don't run much NVE myself, but I *think* I recall seeing that my
colleagues that _do_ set tol=1e-8. The default SHAKE tolerance is 1e-5, so
you've only reduced it by 1 order of magnitude. At some point that will
stop being a significant contribution to the integration error, but I don't
know for sure what point that is.
To be honest, I'm not sure what dsum_tol and rsum_tol *really* do, but
based on a brief code inspection here it is: I think dsum_tol is used to
determine the width of neutralizing charge distributions (i.e., the Ewald
coefficient). It is chosen so that the error in erfc at the cutoff falls
below the tolerance.
This Ewald coefficient, though, has a direct impact on the error in the
reciprocal space sum, too. The narrower you make the Gaussian charge
distributions, the faster the direct space sum decays, but the more
frequencies you need in the reciprocal lattice vectors to fit to the (more
jagged) charge distribution. So even though you made dsum_tol small, the
error in the reciprocal space sum could still be the primary contributor to
your energy drift.
I would suggest *also* decreasing rsum_tol by an order of magnitude (or
more) and see if that helps clean up the drift for you. I *think* that
making both dsum_tol and rsum_tol small will force pmemd to use a small
Ewald coefficient with a large number of reciprocal vectors (enough recip.
vectors to force the reciprocal sum error to fall below rsum_tol). If
that's true, then this should help. If it's not, then I'm outta ideas :).
Good luck,
Jason
--
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Jan 22 2015 - 13:30:02 PST
