GridFF - ProkopHapala/FireCore GitHub Wiki

Main Idea

In typically simulations involving molecules on substrate, the number of surface atoms is substantially higher than the number of atoms in molecule. Therefore the computational effort spend to evaluate pair-wise non-covalent interactions between atoms of the molecule and the substrate are typically a bottleneck in classical force-field simulations, as they are typically more numerous than molecule-to-molecule pairwise-atomic interactions. At the same time, for practical reason the many such simulations keep the substrate rigid. This makes the simulation behavior more controlable, and eliminates necessity to find or develop reliable classical froce-field parameters for the substrate, which are generally less developed than classical force-fields for organic molecules.

With these assumptions, we can achieve substantial performance advantage and even improved accuracy by replacing pair-wise interactins of molecule atoms with substrate atoms, by interaction of the molecule atoms with a single potential $V(\vec{r}_i,T_i)$, where $\vec{r}_i$ is position of an atom and $T_i$ is atomic type. As long as the inter-atomic potentials are additive we can calculate such potential $V(\vec{r}_i,T_i)$ simply by projecting the interaction energy on a 3D grid. This means to evaluate sum of interaction energy between all atoms of the substrate and the test atom $i$ with type $T_i$ placed to set of positions $[ \vec{r}_i ]$.

$V(\vec r_i,T_i)=\sum_j{ E_{i,j}(\vec r_i- \vec r_j|,T_i,T_j) }$

Factorized Potentials

However, evaluating, storing and dispatching of large number of such potentials for each atomic type can be rather cumbesome, inconvenient and inefficient. Assuming grid spacing 0.1A and single-point precission we can estimate this grid to consume 4 MB (100x100x100x4) per cubic nanometer for each atomic type. To eliminate this problem we can exploit the fact that some popular non-covalent inter-atomic potentials, such as Coulomb and Morse can be easily factorized in such a way, that the parameters of the test atom $i$ (corresponding to molecule) can be etracted in front of the sum, which allows to store projected potentials of surface atoms $V_j(r)$ independent of the molecule atom types $T_i$.

This is trivial for Coulomb:

$$ V_{Coulomb} = \sum_j \frac{Q_i Q_j}{|\vec{r} - \vec{r}_j|} = Q_i \sum_j \frac{Q_j}{|\vec{r} - \vec{r}_j|} $$

But it can be rather easily done also for Morse assuming Lorentz-Berthelot mixing rules:

$$ V_{Morse} = V_{Pauli} + V_{Lond} = \sum_j e_i e_j ( \exp(2k( |\vec{r} - \vec{r}_j| - R_i - R_j ) ) - 2 \exp(k(|\vec{r} - \vec{r}_j|- R_i - R_j ) ) ) \newline $$

$$ V_{Pauli} = e_i \sum_j e_j ( \exp(2k( |\vec{r} - \vec{r}_j| - R_i - R_j ) )
= e_i \exp(2kR_i) \sum_j e_j ( \exp(2k( |\vec{r} - \vec{r}_j| - R_j ) ) $$

$$ V_{Lond} = -2 e_i \sum_j e_j \exp(k( |\vec{r} - \vec{r}_j| - R_i - R_j ) )
= -2 e_i \exp(kR_i) \sum_j e_j ( \exp(k( |\vec{r} - \vec{r}_j| - R_j ) )
$$

Hydrogen Bond Correction

We want to implement hydrogen bond correction reflecting the fact that crystal field of strongly polarized atoms can deform electron cloud of neighboring particles. We found that this can be rather well modeled by decreasing repulsive part of Lennarad-Jones or Morse potential. Nevertheless the mixing rules for this correction deviate from Lorentz-Berthelot. In particular we modify the repulsion by factor $1+H_{ij}$ where $H_{ij}=\max(0, H_i H_j)$ (Notice that $H_{ij}$ is always negative and it is non-zero only if $H_i$ and $H_j$ have negative sign. We call $H_i$ and $H_j$ prameters a "hydrogen bond charge" as they typically correlate very well with atomic charge and electronegativity of the particles (more acidic hydrogens has higher $H_i>1$, more electronegative particles with rich free electron pairs like $N,O,F$ has negative $H_i<1$). This mixing rules require us to store the relevant repulsive part of the surface potential separately for each atomic type $j$.

Efficient implementation

Both GPU and CPU architectures are optimized for memory access aligned by 16-64 bytes. In particular OpenCL architecture prefer to use float4 (16 bytes) and float8 (32 bytes), and x86 CPU architecture is optimized for SIMD vector operations such as AVX256 (32 bytes) and AVX512 (64 bytes). Therefore we aim to fit our grid projected surface potential in float4 or double8 data-types.

Fortunately usage of binary substances as a substrate (such as NaCl, MgO, MoS2 ...) allows us to desing a strategy which perfectly fits these harware constrains. We use following memory layout

float4 Vgrid_i{
   V_Paul_A,              // Pauli potential of surface type A
   V_Paul_B,              // Pauli potential of surface type B
   V_Lond_A + V_Lond_B,   // Combined London potential of both A,B (all atom of the substrate)  
   V_Coul                 // Electrostatic poential of all atoms of the substrate
};

This layout allows us to easily apply relevant scaling factor $f_{ij} = e_i (1+\max(0,H_i H_j))$ separately to Pauli repulsion summed over all surface atoms of given type $j$.