|
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
| 
