Computing
After Parallel briey discussing the often neglected, but in praxis frequently encountered, issue of trivially parallel computing, we turn to parallel computing with information exchange. Our illustration is the replica exchange method, also called parallel tempering. Closely related (but less good) is the Jump-Walker (J-Walker) method. Parallel tempering has a low, but non-trivial, communication overhead between the processes. Information about the temperatures, energies and nearest neighbor locations of n systems is communicated, and logical decisions have to be made. This algorithm gives some kind of ideal start into the world of parallel processing and is also of major practical importance. In some applications it competes with the multicanonical algorithm, in others both algorithms supplement one another. For instance, combinations of them have been explored in simulations of peptides (small proteins).
1
Trivially parallel computing
An example of a trivially parallel large scale endeavor is the Folding@Home distributed computing project, which achieved sustained speeds of about thirty teraops in a molecular dynamics simulation There are many tasks which simply can be done in parallel, for instance canonical MC simulations at dierent temperatures and/or distinct lattices. For such applications the gain due to having many (identical) processors available scales linearly with the number of processors. In practice trivially parallel applications are frequently encountered. Another example are simulations of spin glass replica. It is far more cost ecient to run trivially parallel applications on PC clusters than on dedicated parallel machines (supercomputers), which feature expensive communication networks. But, it appears always easier to get CPU time funding on supercomputers than to nd large-scale trivially parallel clusters. For MC simulations of single systems, which last (very) long in real time,
2
parallelization is still trivial, but there is some overhead due to multiple equilibration runs. Typically the CPU time consumption of such a simulation is divided into ttot = teq + tmeas . For parallelization of such a computation on n independent nodes the single processor CPU time scales as t1 = teq + tmeas/n . tot How many processors should one choose (assuming a suciently large supply)? A good rule of thumb appears to be tmeas n = nint teq .
The total processor time used becomes tn = n teq + tmeas 2 tmeas, which limits tot
3
the CPU time loss due to multiple equilibrations to less than two Rloss tn 2 tmeas = tot <2, ttot teq + tmeas
and the gain is a reduction of the run time by the factor
1 Rgain
t1 2 teq 2 tot = = . ttot teq + n teq 1 + n
The worst case is n = 2, for which the improvement factor in CPU time is 2/3. teq . Obviously this parallelization makes most sense for tmeas
4
Replica exchange method (parallel tempering)
The method of multiple Markov chains was introduced by Geyer and, independently, by Hukusima and Nemoto under the name replica exchange method. The latter work was inuenced by the simulated tempering, which coined the name parallel tempering for an exchange of temperatures in multiple Markov chains. But not that simulated tempering is distinct method and actually a special case of the method of expanded ensembles. A precursor of such ideas is a 1986 PRL by Swendsen and Wang, where a general replica exchange method is formulated (apparently its generality prevented one from identifying the important special case). In practice one does not exchange conguration but temperatures or similar parameters, so that the amount of exchanged information is small. Therefore, the method is particularly well suited for parallel processing on clusters of workstations, which lack the expensive fast communication hardware of the
5
dedicated parallel computers. The Jump Walker (J-Walker) approach does not exchange congurations, but feeds replica from a higher temperature into a simulation at a lower temperature, and has consequently problems with balance. Parallel tempering performs n canonical MC simulations at dierent -values with Boltzmann weight factors wB,i(E
(k)
)=e
i E (k)
= eH , i = 0, . . . , n 1 ,
where 0 < 1 < ... < n2 < n1, and allows to exchange neighboring -values: i1 i Their joint weight is ei1Ei1iEi = eH and the i1 i exchange leads to
6
for i = 1, . . . , n 1 .
H
= =
i1Ei
(k)
iEi1 iEi
(k)
(k )
(k)
i1Ei1
(k )
(i i1) Ei
Ei1
(k )
which is accepted or rejected according to the Metropolis algorithm, i.e., with probability one for H 0 and with probability exp( H) for H > 0. The CPU time spent on the exchange should be less than 50% of the total updating time. In practice it is normally much less (a few percent). The i spacing should to be determined so that a reasonably large acceptance rate is obtained for each pair. This can be done by a recursion, which is a slight modication of one by Kerler and Rehberg. Let us denote by apt the acceptance rate of the exchange. Assume i pt m ai > 0 for all i = 1, . . . , n 1. Iteration m is performed with the values i ,
7
m+1 i = 0, . . . , n 1. We dene the i values of the next iteration by m+1 m = 0 0 m+1 m+1 m m and i = i1 + am (i i1) for i = 1, . . . , n 1 , i
where am i =
m
apt,m i
with
m
=
For large acceptance rates apt the distance whereas it i shrinks for small acceptance rates. The denition of the ai coecients guarantees m+1 m n1 = n1, i.e. the initial n1 0 distance is kept. A problem of the recursion is still that the statistical stays noise similar in each iteration step. It may be preferable to combine all the obtained estimates with suitable weight factors wk . The nal estimate after m iterations reads then i =
m k wk i k=1 m k k=1 w
. n1 pt,m m m ) (i i1 i=1 ai m m i i1 is increased,
m m n1 0
.
8
Two reasonable choices for the weights are wk = min {apt,k } , i
i
which determines the weight from the worst exchange acceptance rate and suppresses the statistics of runs with a low wk , and 1 wk = with 2 = 2
n
1 apt,k i
2
,
i=1
which relies on what would be the correct statistical weighting of the acceptance rates, if they were statistically independent (what they are not).
9
Computer implementation
We like to program the exchange using MPI. Initially the indices of the values are chosen to agree with the rank MY ID of the process. The exchanges lead then to the situation for which the indices become an arbitrary permutation of the numbers 0, . . . , n 1. Two approaches are possible: 1. All the n processes send the needed information (the energies of their congurations) to a master processor. The master handles the exchange and sends each process its new value back (which can agree with the old value). 2. The involved processes handle the exchanges through point to point communication. We consider only the second approach, because for it the number of sweeps per time unit has a chance to scale linearly with the number of processors, whereas for
10
the rst approach the master processor will become the bottleneck of the simulation when many processors are employed (when only few processors are used, it may be more ecient to let a master processor handle the exchanges). Each process has to know the identities of the processes where the neighboring values reside. They are stored in a neighbor array NEIGH(2). The initial values are those dened in the phbwrite mpi.f test program of 1MPI. With each exchange the nearest neighbor and next nearest neighbor connections get modied and a correct bookkeeping is the main diculty of implementing the parallel tempering algorithm. We cannot deal with all pairs simultaneously, but need one buer process between each pair. This leads to three distinct interaction groups, starting either with 0, 1 or 2, which we label by INDEX = 0, 1, 2 as indicated for eight processes in the following gure (e.g., for INDEX = 0 the buer processes are i = 2 and i = 5):
11
2 INDEX
1
0 0 1 2 3 i of i 4 5 6 7
Our MPI Fortran subroutine p pt mpi.f implements the exchange (as in the gure the variable INDEX takes on the values 0, 1 or 2). After the lowest value is determined, the higher values follow according to
12
the scheme of the gure. The thus identied processes send their iact, IACPT PT (for apt) and NEIGH(1) values to their NEIGH(2) neighbors, i.e. to the processes i which accommodate their upper, neighboring values. These neighbor processes decide whether the exchange is accepted or rejected, relying on the logical-valued function LPT EX IA(BETA1, BETA2, IACT1, IACT2) , which implements the Metropolis criterion. For either outcome the processes send two integers back to their lower neighbors (the originally initiating processes) and one integer to their upper neighbors. If the exchange is rejected, the rst element of each MPI SEND is set to 2, thus signaling to the receiver that nothing changed. If accepted, the receivers store the addresses of their new neighbors in the NEIGH array. In addition, two more actions are performed: First, the initiating processes inform their left neighbors about the change by sending a new neighbor address to NDEST3R. This send is performed in any case, where the message 2 informs NDEST3R that nothing changed. Second, the processes whose values changed re-calculate their heatbath weights and interchange also their parallel tempering acceptance rates IACPT PT.
13
After a certain number of exchanges, the values are updated by either of the subroutines pt rec0 mpi.f, pt rec1 mpi.f or pt rec2 mpi.f . In practice it may be unavoidable that one (or several) of the acceptance rates apt i are zero. Applying the Kerler-Rehberg recursion blindly gives the undesirable result i1 = i. To avoid this, we replace all apt = 0 values by apt = 0.5/Nupdate i min where Nupdate is the number of updates done. This works when there are enough i values to bridge the [1, n1] interval. Otherwise, the only solution to the problem of apt = 0 values is to increase the number of processes. i We could gather the needed permutation of values, as well as their corresponding exchange acceptance rates, using point to point MPI communication.
14
However, it is shorter to rely on the MPI ALLGATHER (SENDBUF, ISENDCOUNT, SENDTYPE, RECVBUF, IRECVCOUNT, RECVTYPE, COMM, IERROR) instruction. The arguments of MPI ALLGATHER are: SENDBUF, the starting element of the send buer; ISENDCOUNT (integer), the number of elements in the send buer; SENDTYPE, the data type of the send buer elements; RECVBUF, the rst element of the receive buer; IRECVCOUNT (integer), the number of elements in the receive buer; RECVTYPE, the data type of the receive buer; COMM, the communicator; IERROR, the MPI error code.
15
Finally a usefull instruction, in particular for debugging purposes, is MPI BCAST(BUF, ICOUNT, DATATYPE, ISOURCE, COMM, IERROR) which broadcasts the buer BUF from the process with rank ISOURCE to all processes specied by the communicator COMM. The number of elements in the buer is given by ICOUNT, the datatype by DATATYPE, and IERROR is the MPI error code.
16