Fortran Potentials and Ewald Sums

Test programs for potentials, forces and torques

Two program files are provided: test_pot_atom.f90 and test_pot_linear.f90, for pair potentials between, respectively, atoms and linear molecules. These are combined, as appropriate, with modules which contain a subroutine 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 namelist in individual cases; to run the programs without any tweaking, simply give an empty namelist &nml / to standard input in the usual way. The supplied examples are:

• test_pot_at the Axilrod-Teller three-body potential
• test_pot_bend the angle-bending part of a polymer chain potential
• test_pot_dd the dipole-dipole potential
• test_pot_dq the dipole-quadrupole and quadrupole-dipole potential
• test_pot_gb the Gay-Berne potential
• test_pot_qq the quadrupole-quadrupole potential
• test_pot_twist the angle-torsion part of a polymer chain potential

In all cases, the SConstruct file builds these programs in a directory whose name is taken from the module name above, but the executable file is named test_pot_atom or test_pot_linear as appropriate.

T-tensor program

The program t_tensor 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 namelist &nml /, so as to take the program defaults, or various parameters may be specified in the namelist.

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 and we provide a program ewald 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, and elsewhere for the FFT routine which is illustrated in fft3dwrap.

Mesh program

The program mesh 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. 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 namelist parameters.