Parallel simulations on clusters using Message Passing

Introduction

Asap can run on parallel computers and Linux clusters, using the Message Passing Interface to communicate between the processors. Best results are obtained on identical machines, since no load balancing is provided.

The mechanism behind parallel simulations is “Domain Decomposition”: space is divided into regions, and each processor handles its own part of space. Processors are placed in a regular 3D grid, as shown on the figure.

Atoms, that move from one region to another, are transferred between the processors. Each processor also needs to communicate with neighboring processors, so atoms near the boundaries between processors can know where all their neighbors are. To keep communication overhead down to an acceptable level, each processor should handle at least 25000 atoms, preferably at least 100000 atoms.

Each processor is responsible for a part of space, but needs to know the positions of atoms in a slightly larger part of space. In reality, a three-dimensional grid is used instead of the two-dimensional grid shown here.

A processor needs to know the positions of nearby atoms on neighboring processors.

Setting up a parallel simulation

Initially, the atoms must be distributed among the processors. This is handled by the MakeParallelAtoms function, which takes an Atoms object on one or more processors, and returns a special ParallelAtoms object on all processors. The atoms will be distributed correctly among the processors.

MakeParallelAtoms(atoms, cpus, cell=None, pbc=None, distribute=True):

Creates a ParallelAtoms from one or more Atoms ojects. On each processor, atoms is either an Atoms object or None, if more than one processor contributes with atoms, they must contain the same type of information (i.e. if one contains momenta, they all must), and the final ParallelAtoms will contain the atoms from all the original Atoms objects. Duplicated atoms will appear more than once in the results, that will normally be an error.

All processors must specify the same cpus, it is a 3-tuple specifying the size of the grid of processors.

Some rarely-used optional arguments (if used they must be identical on all processors): cell and pbc can be used to override the supercell and the periodicity of the atoms. distribute can be used to prevent the distribution of atoms, it can (rarely!) be useful if you have some extra data to be associated with the atoms, and want to associate that first.

Often, a parallel simulation is started from a file. The Trajectory class has been extended to return a ParallelAtoms object if requested:

Trajectory(filename).get_atoms(n, cpus):

Return Atoms object number n as a ParallelAtoms object, the cpus argument has the same meaning as in MakeParallelAtoms. If cpus is left out, an ordinary Atoms object is returned on the master (cpu number 0) and None on all other cpus.

The ParallelAtoms object

The ParallelAtoms object is derived from the ASE Atoms object, and has some extra functionality required to make parallel simulations work. The following methods are special:

atoms.get_number_of_atoms():

Returns the total number of atoms in the simulation.

len(atoms):

Returns the number of atoms on this processor.

atoms.distribute():

Redistribute the atoms among the cpus, e.g. after adding new atoms.

atoms.get_ids():

Get the ID array. Each atom has a unique ID number between 0 and N-1. This number remains the same as atoms move between processors.

atoms.is_master():

Return True on CPU number 0 and False on all other CPUs.

atoms.get_comm():

Return the communicator. See the Communication using message passing.

Dynamics and filters

Most filters and dynamics objects can directly be used with parallel simulations, provided you import them from Asap instead of importing them directly from the ASE. In many cases, special versions containing small modifications are made to handle atoms that migrate from one processor to another.

Filters, that change the number of atoms may be problematic, in particular if they are combined with a dynamics object which stores information about the atoms. If you need such a combination, and run into problems, please contact the mailing list, most of these problems are solvable, but I do not want to spend time on it until there is a need.

If you need a filter or a dynamics that is not in Asap, please inform the mailing list. Maybe it is trivial to add it. In the meanwhile, try importing it directly from ASE, but please inform the mailing list even if the ASE version works.

File output

File output from parallel simulations require that all information about the atoms is collected at the master node. This is done automatically by the Trajectory object: attaching a Trajectory to a Dynamics object will result in the data being collected and saved correctly, provided the Trajectory object is imported from asap3. The ASE methods for reading and writing files will otherwise only work in serial simulations.

Of course, any file can be read on the master node alone with any ASE method, and then promoted to a parallel atoms object with the MakeParallelAtoms() method described above.

Please see the example Running on parallel computers.

Common Neighbor Analysis and Polyhedral Template Matching

These analysis methods work out of the box, provided that the cutoff given to the method is shorter than the cutoff of the potential, as the necessary information is otherwise not available on this processor.

In particular for PTM this typically means that a cutoff must be specified, since the default (10 Å) will typically be too long.

Special considerations in your own scripts

Most of your script should be almost the same whether you run a serial or a parallel Asap, except that you need to set up the parallel simulation, as described above. There are, however, a few pitfalls to avoid. In particular, it is important that all processors always execute the same commands, as communication between processors require them to be at the same place in the code.

Migration occurs while calculating forces or energies

Atoms migrate between processors when forces or energies are calculated (but not necessarily every time they are calculated). It means that the number of atoms on a given processor may change when get_forces(), get_potential_energies() or get_potential_energy() are called.

This expression will fail with a mysterious error message from Python (frames not aligned or something cryptic like that):

new_p = atoms.get_momenta() + dt * atoms.get_forces()

First, Python calls get_momenta, which returns an array containing the momenta. Then it calls get_forces, but before it returns, atoms may have migrated between processors; so when it returns the number of atoms may have changed, and the size of the force array will be different from the size of the momentum array. The solution is to get the forces first, forcing migration to happen then:

f = atoms.get_forces()
new_p = atoms.get_momenta() + dt * f

The number of atoms is ambiguous

The number of atoms in the simulation and the number of atoms on this processor are different quantities in a parallel simulation.

len(atoms) returns the number of atoms on this processor. This will vary as atoms migrate.

atoms.get_number_of_atoms() returns the total number of atoms on all processors. This is a constant.

Be careful which one you use. Use len(atoms) when looping over atoms, and when dimensioning arrays, but calculate the energy per atom like this: atoms.get_potential_energy() / atoms.get_number_of_atoms().

Beware of file corruption

If multiple processors open the same file for writing on an NFS mounted filesystem, file corruption is almost certain, even if they write the same data. So avoid that. The Trajectory objects take care of only opening the output file on the master processor, so this is most likely to happen with files you open yourself. Se the example in the next section.

If you write data directly from your own script, consider opening the file with ASE’s paropen function. In that way only one processor writes to the file.

On the CAMD cluster Niflheim, files in /scratch/username are local to each processor. It may be a useful place to dump debugging information.

Beware of conditionals

Beware of conditionals such as if statements in your script. If the different processors take different branches, your simulation will most likely lock up and hang. A typical error is to attempt to print something on the “master” processor, but not on the rest:

master = (asap3.mpi.world.rank == 0)
  ...
if master:
    output = open("filename.txt", "w")

  ...
if master:
    output.write("Energy = " + str(atoms.get_potential_energy()))
    # This will FAIL !

The problem is that the master processor calls get_potential_energy(), but the other processors do not. When calculating the total potential energy, communication between the processors is necessary, and the master processor will wait forever for answers from the other processors. A better way of doing this is:

master = (asap3.mpi.world.rank == 0)
  ...
if master:
    output = open("filename.txt", "w")
else:
    output = open("/dev/null", "w")

  ...
output.write("Energy = " + str(atoms.get_potential_energy()))

Since this is such a common construction, the ASE and Asap contain a helper function for opening files in parallel simulations, imported into asap3 in parallel simulations:

output = paropen("filename.txt", "w")

  ...
output.write("Energy = " + str(atoms.get_potential_energy()))

Keeping data between time steps

Often, you want to compare data at two different times, for example to calculate how far atoms have moved. Keeping data from previous times will usually result in a \(frames not aligned\) error message when you try to use it. It is because the number of atoms, and the indexes of atom, keep changing every time atoms move between processors. If you get (or calculate) data concerning the atoms at one time, and try to use it at another time, the atoms will have migrated between processors, and the one-to-one correspondence between your data and the atoms is lost. The solution is to store the data on the atoms, letting it be migrated along with the atoms.

Most kinds of numerical data can be stored on the atoms in the form of a numpy array. Store the data with atoms.set_array(name, data) and receive it again with atoms.set_array(name).

Be sure to store all data before a call to get_forces(), get_potential_energies() or get_potential_energy(), and receive it again afterwards. Any arrays in your script will become invalid when calling these functions, so although it may look strange to store the data with set_array and retrieve it almost immediately after with get_array, it is necessesary to preserve data during a migration.

The \(registered arrays\) available in Asap version 2 are no longer available. It was a more complicated solution to the same problem.

Uncaught Python exceptions and debugging information

If an exception occurs on a processor causing the script to stop, Asap attempts to stop the script on all processors. If the real error did not occur on the master processor, there is a risk that the error message is lost, and the only error message is from stopping the scripts. To make it possible to debug such cases, a function is provided which redirects output on all processors into files.

DebugOutput(filename, stdout=1, nomaster=False):

Redirect stderr to a different file on each node. The filename should contain %d, which is replaced by the node number. The file is opened with minimal buffering, and stderr (standard error) is redirected to it (also for C/C++ extensions). If the optional argument stdout is true (the default), Python’s sys.stdout is also redirected to the same file. Standard output for C/C++ extensions cannot be redirected. If the optional argument nomaster is set to true, redicection does not occur on the master processor.

Deleting atoms

In a parallel simulation, you cannot delete individual atoms with the del operator, the way you can in serial simulations. The reason is that deletion needs coordination amongst the tasks on different processors.

Starting with version 3.12.4, atoms can be deleted with the delete_atoms_globally method of the atoms.

delete_atoms_globally(self, local_indices=None, global_indices=None):

Delete atoms from a parallel simulation.

local_indices: Remove atoms with these indices on the local MPI task.

global_indices: Remove atoms with these global IDs. The atoms may reside on this or on another MPI task.

Both parameters should be sequences, or a single int. Both can be specified simultaneously, if need be, although that may be confusing.

This method must be called simultaneously by all MPI tasks, but the lists of indices may be diffent on the different tasks. An atom is deleted if it is specified in at least one task. Duplicates are allowed.

Example

Please see the example Running on parallel computers in the Examples section.