Amber masthead
Filler image AmberTools21 Amber20 Manuals Tutorials Force Fields Contacts History
Filler image

Useful links:

Amber Home
Download Amber
Installation
Amber Citations
GPU Support
Updates
Mailing Lists
For Educators
File Formats
Materials for Educators → Particle Mesh Ewald Electrostatics

Download the library of Matlab or GNU Octave functions for this demonstration.

Unwrapping the Library

Unpack the tarball in its own directory, or directly in your MATLAB documents directory. Open Matlab or Octave; from here onward, unless otherwise noted, all commands will be given on the interpreter command line, not the Linux or Windows shell command line. If you use Matlab and install in a unique directory, type path(path, '⋖your directory⋗') in the interpreter or in any script that will call the functions to put them in the path. For octave, use the addpath() function.

The first thing to do is to make sure that things are working. Run the TestAll script on the interpreter command line. The results should appear roughly as follows:

This test concerns a diffuse system of ions in a cubic box.

Error in PME direct space calculation:  0.000110 kcal/mol-A
Error in PME recip. space calculation:  0.000004 kcal/mol-A
Error in PME complete calculation:      0.000110 kcal/mol-A

All errors should be less than 0.001 kcal/mol-A in this example.
For condensed-phase systems, errors of up to 0.01 kcal/mol-A are acceptable.

This test concerns a diffuse system of ions in a truncated octahedron.

Error in PME direct space calculation:  0.000320 kcal/mol-A
Error in PME recip. space calculation:  0.000011 kcal/mol-A
Error in PME complete calculation:      0.000321 kcal/mol-A

All errors should be less than 0.005 kcal/mol-A in this example.
For condensed-phase systems, errors of up to 0.01 kcal/mol-A are acceptable.

The examples driven by TestAll, drawing upon coordinate and charge information in the Tests directory, are examples of the forces that emerge from a PME calculation on these systems of ions. Presently, the code does not handle exclusions, but the only way that PME handles such things is to calculate everything in the manner that this library does and then subtract off the contributions of the excluded interactions. A more detailed implementation, based on this Matlab library, can be found in the mdgx simulations engine distributed with AmberTools.

Another nice feature of the library is the extensive comments added to every function. Open any of them in the Matlab or Octave editor (or vi or emacs) and the purpose and usage should become clear. One can also type help ⋖function name⋗ in the interpreter to get the critical information. Note that many functions return more than one argument if prompted. This is another useful feature of the Matlab code to allow users to reach in at any point in the workflow and pull out information of interest.

A central quantity in Ewald calculations is the box geometry, specified by six independent pieces of information. Intuitively, we speak in terms of the box edge lengths and three major angles, as described for crystal lattices. These values in turn determine the six unique elements of 3 x 3 upper-triangular matrices U and invU to take coordinates into and back out of expressions based on units of box stacking (fractional coordinates). The CompXfrm function will create transformation matrices based on the six box dimensions, such as those found on the last line of a periodic inpcrd file (remember to scale the final three numbers, box angles, by π/180.0). Example:

gdim = [70.0 60.0 60.0 pi/2.0 2.0/3.0*pi pi/1.8];
[U, invU] = CompXfrm(gdim);

U =  0.0143    0.0025    0.0085      invU = 70.0000  -10.4189  -30.0000
          0    0.0169    0.0017                   0   59.0885   -5.2898
          0         0    0.0193                   0         0   51.6916

The transformation matrices are much more convenient mathematical quantities for manipulating coordinates. Because the regular grid tilts and stretches to fit the box geometry, expressing all coordinates in fractional coordinates is a useful step before finding the alignment of each particle to the mesh.

Go to the next section

Return to the Introduction

"How's that for maxed out?"

Last modified: Aug 17, 2018