3.3 Parallelization levels

In QUANTUM ESPRESSO several MPI parallelization levels are implemented, in which both calculations and data structures are distributed across processors. Processors are organized in a hierarchy of groups, which are identified by different MPI communicators level. The groups hierarchy is as follow:

• world: is the group of all processors (MPI_COMM_WORLD).
• images: Processors can then be divided into different "images", each corresponding to a different self-consistent or linear-response calculation, loosely coupled to others.
• pools: each image can be subpartitioned into "pools", each taking care of a group of k-points.
• bands: each pool is subpartitioned into "band groups", each taking care of a group of Kohn-Sham orbitals (also called bands, or wavefunctions) (still experimental)
• PW: orbitals in the PW basis set, as well as charges and density in either reciprocal or real space, are distributed across processors. This is usually referred to as "PW parallelization". All linear-algebra operations on array of PW / real-space grids are automatically and effectively parallelized. 3D FFT is used to transform electronic wave functions from reciprocal to real space and vice versa. The 3D FFT is parallelized by distributing planes of the 3D grid in real space to processors (in reciprocal space, it is columns of G-vectors that are distributed to processors).
• tasks: In order to allow good parallelization of the 3D FFT when the number of processors exceeds the number of FFT planes, FFTs on Kohn-Sham states are redistributed to "task" groups so that each group can process several wavefunctions at the same time.
• linear-algebra group: A further level of parallelization, independent on PW or k-point parallelization, is the parallelization of subspace diagonalization / iterative orthonormalization. Both operations required the diagonalization of arrays whose dimension is the number of Kohn-Sham states (or a small multiple of it). All such arrays are distributed block-like across the ``linear-algebra group'', a subgroup of the pool of processors, organized in a square 2D grid. As a consequence the number of processors in the linear-algebra group is given by n², where n is an integer; n² must be smaller than the number of processors in the PW group. The diagonalization is then performed in parallel using standard linear algebra operations. (This diagonalization is used by, but should not be confused with, the iterative Davidson algorithm). The preferred option is to use ScaLAPACK; alternative built-in algorithms are anyway available.

Note however that not all parallelization levels are implemented in all codes!

3.3.0.1 About communications

Images and pools are loosely coupled and processors communicate between different images and pools only once in a while, whereas processors within each pool are tightly coupled and communications are significant. This means that Gigabit ethernet (typical for cheap PC clusters) is ok up to 4-8 processors per pool, but fast communication hardware (e.g. Mirynet or comparable) is absolutely needed beyond 8 processors per pool.

3.3.0.2 Choosing parameters

: To control the number of processors in each group, command line switches: -nimage, -npools, -nband, -ntg, -ndiag or -northo (shorthands, respectively: -ni, -nk, -nb, -nt, -nd) are used. As an example consider the following command line:

mpirun -np 4096 ./neb.x -ni 8 -nk 2 -nt 4 -nd 144 -i my.input

This executes a NEB calculation on 4096 processors, 8 images (points in the configuration space in this case) at the same time, each of which is distributed across 512 processors. k-points are distributed across 2 pools of 256 processors each, 3D FFT is performed using 4 task groups (64 processors each, so the 3D real-space grid is cut into 64 slices), and the diagonalization of the subspace Hamiltonian is distributed to a square grid of 144 processors (12x12).

Default values are: -ni 1 -nk 1 -nt 1 ; nd is set to 1 if ScaLAPACK is not compiled, it is set to the square integer smaller than or equal to half the number of processors of each pool.

3.3.0.3 Massively parallel calculations

For very large jobs (i.e. O(1000) atoms or more) or for very long jobs, to be run on massively parallel machines (e.g. IBM BlueGene) it is crucial to use in an effective way all available parallelization levels. Without a judicious choice of parameters, large jobs will find a stumbling block in either memory or CPU requirements. Note that I/O may also become a limiting factor.

Since v.4.1, ScaLAPACK can be used to diagonalize block distributed matrices, yielding better speed-up than the internal algorithms for large ( > 1000 x 1000) matrices, when using a large number of processors (> 512). You need to have -D__SCALAPACK added to DFLAGS in make.sys, LAPACK_LIBS set to something like:

    LAPACK_LIBS = -lscalapack -lblacs -lblacsF77init -lblacs -llapack

The repeated -lblacs is not an error, it is needed! configure tries to find a ScaLAPACK library, unless configure -with-scalapack=no is specified. If it doesn't, inquire with your system manager on the correct way to link it.

A further possibility to expand scalability, especially on machines like IBM BlueGene, is to use mixed MPI-OpenMP. The idea is to have one (or more) MPI process(es) per multicore node, with OpenMP parallelization inside a same node. This option is activated by configure -with-openmp, which adds preprocessing flag -D__OPENMP and one of the following compiler options:

ifort	-openmp
xlf	-qsmp=omp
PGI	-mp
ftn	-mp=nonuma

OpenMP parallelization is currently implemented and tested for the following combinations of FFTs and libraries:

internal FFTW copy	requires -D__FFTW
ESSL	requires -D__ESSL or -D__LINUX_ESSL, link with -lesslsmp

Currently, ESSL (when available) are faster than internal FFTW.

3.3.1 Understanding parallel I/O

In parallel execution, each processor has its own slice of data (Kohn-Sham orbitals, charge density, etc), that have to be written to temporary files during the calculation, or to data files at the end of the calculation. This can be done in two different ways:

• ``distributed'': each processor writes its own slice to disk in its internal format to a different file.
• ``collected'': all slices are collected by the code to a single processor that writes them to disk, in a single file, using a format that doesn't depend upon the number of processors or their distribution.

The ``distributed'' format is fast and simple, but the data so produced is readable only by a job running on the same number of processors, with the same type of parallelization, as the job who wrote the data, and if all files are on a file system that is visible to all processors (i.e., you cannot use local scratch directories: there is presently no way to ensure that the distribution of processes across processors will follow the same pattern for different jobs).

Currently, CP uses the ``collected'' format; PWscf uses the ``distributed'' format, but has the option to write the final data file in ``collected'' format (input variable wf_collect) so that it can be easily read by CP and by other codes running on a different number of processors.

In addition to the above, other restrictions to file interoperability apply: e.g., CP can read only files produced by PWscf for the k = 0 case.

The directory for data is specified in input variables outdir and prefix (the former can be specified as well in environment variable ESPRESSO_TMPDIR): outdir/prefix.save. A copy of pseudopotential files is also written there. If some processor cannot access the data directory, the pseudopotential files are read instead from the pseudopotential directory specified in input data. Unpredictable results may follow if those files are not the same as those in the data directory!

IMPORTANT: Avoid I/O to network-mounted disks (via NFS) as much as you can! Ideally the scratch directory outdir should be a modern Parallel File System. If you do not have any, you can use local scratch disks (i.e. each node is physically connected to a disk and writes to it) but you may run into trouble anyway if you need to access your files that are scattered in an unpredictable way across disks residing on different nodes.

You can use input variable disk_io to reduce the the amount of I/O done by pw.x. Since v.5.1, the dafault value is disk_io='low', so the code will store wavefunctions into RAM and not on disk during the calculation. Specify disk_io='medium' only if you have too many k-points and you run into trouble with memory; choose disk_io='none' if you do not need to keep final data files.

For very large cp.x runs, you may consider using wf_collect=.false., memory='small' and saverho=.false. to reduce I/O to the strict minimum.