This code goal is to develop a simple microscopic simulation based on the Lennard-Jones potential Default build :
$ cmake -S . -B Build -G Ninja
$ cd Build
$ ninja In order to enable optimizations :
$ cmake -S . -B Build -G Ninja -DCMAKE_BUILD_TYPE=Release
$ cd Build
$ ninja To run the default configuration
Other usefull targets are :
$ ninja doc
$ ninja plot
$ ninja vmd A file containing cartesian coordinates of a particle system.
The format is the following :
> 0 1 > type x y z > type x y z First line is the number of input types
Type is ignored here
When the input is parsed we need to be checking for duplicates particles
as it's physically wrong to have two particles at the same place
These are the main formulas used in the simulation
$$\tilde{U}^{LJ} = \frac{4 \epsilon^*}{2}
\sum_{i=1}^{N}
\sum_{j > j}
\left[
\left( \frac{\sigma}{r} \right)^{12}
- 2 \times \left( \frac{\sigma}{r} \right)^{6}
\right]
= \sum_{i=1}^{N}
\sum_{j > j}
u_{ij}$$
Lennard Jones Potential : $$U^{LJ}(r) = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12}
- 2 \times \left( \frac{\sigma}{r} \right)^{6} \right]
+ \delta U$$
Forces computations (gradient of Lennard Jones) : $$\frac{\partial u_{ij}}{\partial x_i} =
-48 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12}
- \left( \frac{\sigma}{r} \right)^{6} \right]
\times \frac{x_i - x_j}{r^2_{ij}}$$
$$\left\{
\begin{array}{l}
P_5(r) = 1
- 10 \left( \frac{r - r_{min}}{r_{max} - r_{min}} \right)^{3}
+ 15 \left( \frac{r - r_{min}}{r_{max} - r_{min}} \right)^{4}
- 6 \left( \frac{r - r_{min}}{r_{max} - r_{min}} \right)^{5} \\\
\frac{\partial P_5}{\partial x_i} =
- \frac{30}{r_{max} - r_{min}} \left( \frac{r - r_{min}}{r_{max} - r_{min}} \right)^{2}
+ \frac{60}{r_{max} - r_{min}} \left( \frac{r - r_{min}}{r_{max} - r_{min}} \right)^{3}
- \frac{30}{r_{max} - r_{min}} \left( \frac{r - x_{min}}{x_{max} - x{min}} \right)^{4}
\times \frac{x_j - x_i}{r} \\\
x_{min} \leq x \leq x_{max} \\\
(x_{min} = R_c - \delta_r, x_{max} = R_c + \delta_r)
\end{array}
\right .$$
Then the correction Lennard Jones and Forces with sigmoid for points between
R_c - dr < r_ij < R_c + dr
to keep the equation continuous and differentiable (necessary as we are
considering a cut off radius)
$$U^{LJ}(r) = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12}
- 2 \times \left( \frac{\sigma}{r} \right)^{6} \right] \times P_5(r)$$
Forces : (gradient of the new L_J term) $$\frac{\partial u^s_{ij}}{\partial x_i} =
\frac{\partial u_{ij}}{\partial x_i} \times P_5(r)
+ \frac{\partial P_5}{\partial x_i} \times U_{LJ}$$
Periodic Conditions and main formula : $$\tilde{U}^{LJ}(R_{cut}) = \frac{4 \epsilon^*}{2}
\sum_{\vec{\eta}}^{N_{sym}}
\sum_{i=1}^{N}
\sum_{j, r_{ij} \leq R_{cut}}
\left[
\left( \frac{\sigma}{r} \right)^{12}
- 2 \times \left( \frac{\sigma}{r} \right)^{6}
\right]$$
Velocity Verlet
$$\left\{
\begin{array}{l}
m_i \frac{\partial \vec{r_i}}{\partial t}(t + \frac{1}{2} \partial t) =
m_i \frac{\partial \vec{r_i}}{\partial t}(t)
- \frac{1}{2} \vec{\nabla}_i U(t)
\cdot \partial t \\\
\vec{r}_i (t + \partial(t)) =
\vec{r}_i (t) + \frac{\partial \vec{r}_i}{\partial t} (t + \frac{1}{2} \partial(t)) \cdot \partial t \\\
m_i \frac{\partial \vec{r_i}}{\partial t} (t + \partial t) =
m_i \frac{\partial \vec{r_i}}{\partial t}(t + \frac{1}{2} \partial t)
- \frac{1}{2} \vec{\nabla}_i U(t + \partial t)
\cdot \partial t \\\
\end{array}
\right .$$
$$\frac{dT}{dt} = \frac{T_0 - T}{\tau}$$
PDB formatted file with the simulation state at different time steps,
by default results are dumped every iterations
CSV for the Temperatures, Kinetic Energy, Potential Energy, Sum of the energies
by default results are also dumped every iterations
https://fortranwiki.org/fortran/show/Keywords https://fortran-lang.org/en/learn/ https://en.wikipedia.org/wiki/Lennard-Jones_potential https://en.wikipedia.org/wiki/Periodic_boundary_conditions https://en.wikipedia.org/wiki/Verlet_integration https://www.biostat.jhsph.edu/~iruczins/teaching/260.655/links/pdbformat.pdf https://curc.readthedocs.io/en/latest/programming/OpenMP-Fortran.html
