1
$\begingroup$

First, let just me stress that I'm not a an expert in computation chemistry, so now the problem:

We have GCMC molecular simulation, in the Grand Canonical ensemble, using the standard metropolis algorithm.

We would like to port the program to CUDA and run it in parallel. There are several techniques to parallelize such simulation to system with local interaction. Unfortunately in our simulation we compute a total energy for the system taking into account all the particles in the system. We use a primitive model, and given the elementary charge $e$ and the absolute permittivity of vacuum $\epsilon_0$, the interaction between two charges species $i$ and $j$ separated by a distance $r_{ij}$ is calculated acconrding to: $$\begin{cases} u^{\text{el}}(r_{ij}) = \frac{q_iq_je^2}{4\pi\epsilon_0\epsilon_sr_{ij}}\\ u^{\text{hs}}(r_{ij}) = \infty,& \quad \text{if } r_{ij} \leq R_i + R_j\\ u^{\text{hs}}(r_{ij}) = 0, &\quad \text{if } r_{ij} > R_i + R_j\\ \end{cases} $$

The hard radii sphere is set to $2\mathring{A}$. A one-body external field was use to used to correct for the long range electrostatic interaction and to account for the cell boundary contraints along the z-direction acting as hard walls.: $$ \begin{cases} v^{\text{ext}}(z_i) = \infty,& \quad \text{when } z_i \geq L_z||\;z_i \leq <0\\ v^{\text{ext}}(z_i) = q_ie\varphi^{\text{ext}}(z_i), &\quad \text{when } 0 < z_i < L_z\\ \end{cases} $$ where $\varphi^{\text{ext}}(z_i)$ is the potential produced by the external charge distribution.

The full configurational energy of the system is then: $$ U^\text{conf}= \sum _i^{N_0+N_s} v^\text{ext}(z_i) + \sum_{1<j}^{N_0+N_s} u^\text{el}(r_{ij}) + u^\text{hs}(r_{ij})$$

The heaviest part of the simulation the the calculation of $U^\text{conf}$.

The main iteration is like this:

for i<-ion_0 to last_ion do:
   actual_energy <- compute_u_conf(i)
   propose_a_move(i)
   new_energy <- compute_u_conf(i)
   accept_or_reject_move(i)

Here is where domain decomposition does not fit. Even if we can compute in parallel a move, then the computation of the energy would not be consistent because others particles are moving. This would break the Markov chain and make the system not convergent because the acceptance of the proposed move is based also on the total energy calculated.

It is possible to calculate the energy only locally and using some sort of mean for the rest of the space? Is there a way to parallelize the calculation of the energy being consistent? Have someone any idea on who can it be done? Any help is welcome.

Improve this question
$\endgroup$
3
  • $\begingroup$ I am not sure, but do you really need to go towards using domain decomposition? You could use CUDA and run it in on the GPU, but still do everything in a "shared memory context". To me, domain decomposition only becomes necessary if you want to go CUDA + MPI on multi-GPU... then your problem occurs... $\endgroup$
    – BlaB
    Jun 15, 2017 at 12:32
  • $\begingroup$ @BlaisB: Your approach of doing everything in a "shared memory context" would allow you to run only one CUDA block, because thread synchronization applies only to a block of threads, not the whole grid. This is not enough for an efficient implementation. $\endgroup$
    – Jakub Klinkovský
    Jun 15, 2017 at 13:54
  • $\begingroup$ @Blais: The problem is not the degree of parallelization, nor a shared memory problem. The problem is that in this simple Metropolis algorithm the energy is computed in sequence for each ion after the acceptance/rejection of the preview move. This ensure that system will converge. This is clearly not possible in parallel, we can only compute the energy based on a prior global state. So where my question come from: will the system still converge if we try to move x ions at a time based on a prevew global state and then make a global acceptance/rejection? $\endgroup$
    – dual0
    Jun 16, 2017 at 8:58

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.