[gmx-users] Accuracy of force calculation and doubts on PME parameter choice

[Berk Hess]

Hi,

A test case like with two atoms is an unrealistic system.
For two atoms one could use a cut-off of half the box length and accurate Ewald summation
and gets things right. But you get things right for a periodic replication of two atoms
which is never a realistic system.

For many realistic systems the LJ error might be more problematic than the Coulomb error.
All LJ forces beyond the cut-off are attractive and therefore you miss a large amount
of attractive energy in a condensed system.

The Coulomb forces on the other hand will have an inaccuracy on the order of 1 kJ/mol nm,
but these are local inaccuracies. The potential will not deviate systematically, but fluctuate
around the correct value when looking across the whole box volume. Moreover, the Coulomb
forces will have differing signs, meaning that for a condensed system there will be some
cancellation of errors.

What accuracy you need for the Coulomb forces can differ strongly depending on the system
and the questions you are asking. Point charges in bio-molecular force fields are just approximations
and the error of the approximation is much higher that the force error.
The real question to answer is if for the particular property you are interested in
a certain PME force error would give you a systematic deviation in that property.

Berk.


Date: Mon, 7 Jul 2008 23:25:37 -0400
From: nicklefan at gmail.com
To: gmx-users at gromacs.org
Subject: [gmx-users] Accuracy of force calculation and doubts on PME	parameter choice

Dear all:
 
In most studies, we focus on the accuracy of potential energy. How about the force? Consider two systems:
 
1. Two LJ particles (sigma=0.33nm, epsilon=0.625kJ/mol): If the force is cutoff at 2.5*sigma, so the maximum error in force is at the cutoff length and equal to ~0.08 kJ/mol nm.
 
2. Two point charges in a 3*3*3.3nm box (charge separation: 0.5 nm, q1=e, q2=-e): I tried the following calculations:
a. Ewald summation, rcutoff=1.0nm, ewald_tol=1e-5, FFT_spacing=0.1nm: This produces: f1=-f2=542.2 kJ/mol nm;
b. PME, rcutoff=1.0 nm, ewald_tol=1e-5; 4th order, FFT_spacing=0.1nm, this produces: f1=-f2=543.9 kJ/mol nm;
c. PME, rcutoff=1.0 nm, ewald_tol=1e-5; 4th order, FFT_spacing=0.11nm, this produces: f1=543.97 kJ/mol nm, f2=-544.96 kJ/mol nm;
d. PME, rcutoff=1.0 nm, ewald_tol=1e-5; 4th order, FFT_spacing=0.12nm, this produces: f1=543.43 kJ/mol nm, f2=-545.49 kJ/mol nm;
 
Clearly, for LJ particles, the maximum accuracy is ~0.08kJ/mol nm - which sounds small. But for electrostatics, the error is quite large: if we take the Ewald results to be the exact solution, then for choice d (which seems to be a popular choice among gmx-users), the error amounts to be ~ 2.0 kJ/mol nm. Is this error acceptable?

 
Thanks for your time.
 
Nick

_________________________________________________________________
Express yourself instantly with MSN Messenger! Download today it's FREE!
http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20080708/077f5e24/attachment.html>