Python Potentials and Ewald Sums

Test programs for potentials, forces and torques

Two program files are provided: test_pot_atom.py and test_pot_linear.py, for pair potentials between, respectively, atoms and linear molecules. These load, at runtime, a module containing a function to calculate the necessary potential, forces and torques. The aim is to demonstrate the numerical testing of the analytical derivatives which go into the forces and torques: small displacements and rotations are applied in order to do this. The test is performed for a randomly selected configuration. Some parameters are used to prevent serious overlap, which might produce numerical overflow, while keeping the particles close enough together to give non-zero results. The values of these parameters may be adjusted via the input file in individual cases. To run the programs without any tweaking, simply give (through standard input in the usual way) a record containing the bare minimum information, namely a string which identifies the model of interest, for example {"model":"at"}. The supplied examples (with their identifying strings) are, for test_pot_atom.py:

• "at" test_pot_at.py the Axilrod-Teller three-body potential
• "bend" test_pot_bend.py the angle-bending part of a polymer chain potential
• "twist" test_pot_twist.py the angle-torsion part of a polymer chain potential

and for test_pot_linear.py

• "dd" test_pot_dd.py the dipole-dipole potential
• "dq" test_pot_dq.py the dipole-quadrupole and quadrupole-dipole potential
• "gb" test_pot_gb.py the Gay-Berne potential
• "qq" test_pot_qq.py the quadrupole-quadrupole potential

T-tensor program

The program t_tensor.py compares the calculation of multipole energies by two methods: using explicit formulae based on trigonometric functions of the Euler angles, and via the Cartesian T-tensors. Two linear molecules are placed in random positions and orientations, within a specified range of separations, and some of the contributions to the electrostatic energies and forces are calculated. The program may be run using an empty input record {}, so as to take the program defaults, or various parameters may be specified using JSON format.

Several of the tensor manipulations are neatly expressed using NumPy library functions such as outer (outer product) and einsum (Einstein summation).

The force between the molecules is calculated from the analytical derivative of the T-tensor with respect to the separation vector. This naturally leads to formulae where the original T-tensor of rank n is replaced by one of rank n+1.

The torque on each molecule is calculated by formulae similar to those used for torques on multipoles in an external field, field gradient, etc., but in which the field terms are replaced by tensors based on T and the multipoles on the other molecule. This naturally leads to formulae involving the Levi-Civita (antisymmetric) symbol.

In practical applications, the formulae would usually be incorporated in a scheme for handling long-range forces in periodic boundaries (e.g. Ewald sum).

Ewald program

The k-space and r-space contributions to the Ewald sum are illustrated in ewald_module.py and we provide a program ewald.py to test these. The program reads in a configuration file cnf.inp in the usual format: any of the Lennard-Jones or hard-sphere configurations would be suitable. Charges are assigned to the atoms in an arbitrary way. The program itself adds the surface term (the self term is included in the k-space routine). Then, a comparison is made with the brute force summation over all pairs in shells of periodic boxes surrounding the central box. For default parameters, and test configurations with N=256, reasonable convergence is obtained within 8-10 shells. One can adjust the screening parameter kappa within reason (and the number of k-vectors may need changing as well): the contributions of r-space and k-space terms will change, but their sum should remain approximately constant.

There is also a comparison with a simplified particle-mesh Ewald method. As discussed in the text, the charge distribution is assigned to a cubic mesh, Fourier transformed by FFT, and used to calculate the total potential energy, using the solution of Poisson's equation in Fourier space. In doing so, accuracy is improved by optimizing the so-called influence function G. In this example, we use a simple sharpening function discussed by

• V Ballenegger, JJ Cerda, C Holm, J Chem Theo Comp, 8, 936 (2012),

but more sophisticated optimized functions are possible. It is easy to comment out this sharpening function, to see the extent of the correction; it is reasonably significant for the default parameter values.

See below for more discussion of the mesh function, provided in mesh_module.py, and elsewhere for the FFT routine which is illustrated in fft3dwrap.py.

Mesh program

The program mesh.py generates a random configuration of a small number of charges and illustrates the way this may be assigned to a regular cubic mesh using the triangular-shaped cloud distribution described in

• RW Hockney, JW Eastwood, Computer simulation using particles (Adam Hilger, Bristol, 1988).

The function generating the charge density is provided in mesh_module.py. The mesh dimension is, by default, kept small enough to print out the whole array for inspection afterwards. The number of charges and mesh dimension may be adjusted by the user, via JSON parameters.