<jupyter_start><jupyter_text>Magnetized PlasmaThis codes tests the magnetic integrator by reproducing results of Bernu Journal de Physique, **42** L253 (1981). This notebook can be found in the `sarkas/notebook` folder of the Github repo.We simulate a One Component Plasma composed of $N = 1000$ Hydrogen atoms at a density of $n = 1.62 \times 10^{30}$ $N/m^3$, and temperature $T = 0.3$ $eV$, which leads to $\Gamma = 100$. The system is in a constant magnetic field $B = B_0 \hat{\mathbf z}$ with $B_0 = 17 \times 10^6$ Gauss. The ratio of the cyclotron frequency to the plasma frequency is $\beta \sim 1.0$. The YAML input file can be found at [input_file](https://raw.githubusercontent.com/murillo-group/sarkas/master/docs/examples/Magnetized_OCP/input_files/mag_ocp_mks.yaml) and this notebook at [notebook](https://raw.githubusercontent.com/murillo-group/sarkas/master/docs/examples/Magnetized_OCP/Magnetized_Plasma.ipynb).<jupyter_code># Import the usual libraries
%pylab
%matplotlib inline
import os

plt.style.use('MSUstyle')
# Import sarkas
from sarkas.processes import Simulation, PostProcess, PreProcess

# Create the file path to the YAML input file
input_file_name = os.path.join('input_files',  'mag_ocp_mks.yaml')<jupyter_output>Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib<jupyter_text>PreProcessing<jupyter_code>preproc = PreProcess(input_file_name)
preproc.setup(read_yaml=True)
preproc.run(timing=True, loops=100, pppm_estimate = True)<jupyter_output>[38;2;83;80;84m






     _______.     ___      .______       __  ___      ___           _______.
    /       |    /   \     |   _  \     |  |/  /     /   \         /       |
   |   (----`   /  ^  \    |  |_)  |    |  '  /     /  ^  \       |   (----`
    \   \      /  /_\  \   |      /     |    <     /  /_\  \       \   \    
.----)   |    /  _____  \  |  |\  \----.|  .  \   /  _____  \  .----)   |   
|_______/    /__/     \__\ | _| `._____||__|\__\ /__/     \__\ |_______/    
                                                                            

[0m
An open-source pure-python molecular dynamics code for non-ideal plasmas.



* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
                                 Preprocessing                                  
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

Job ID:  ocp_mag
Job directory:  Simulations/ocp_mag
Equilibration dumps [...]<jupyter_text>PreProcessing shows that we have chosen good parameters for the pppm algorithm. The black dot in the first plot shows our choice of parameters. This leads to a Total Force Error $\sim 9 \times 10^8$. The horizontal and vertical lines in the above plots represent the value $r = L/2 a_{\rm ws}$. A common technique used in magnetized MD is to equilibrate the system in two steps. First without the presence of the magnetic field and then with the magnetic field turned on. This is because, the magnetic field can greatly delay the relaxation to equilibrium see [Ott et al PRE 2017](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.95.013209). Hence the magnetization phase before the production phase.We confirm that our choice of parameters leads to the same physical system as the one simulated by Bernu. $\Gamma =95$, $\beta = 1.0$. Contrarily to Bernu we choose to run our simulation for a much longer time.     Production:     No. of production steps = 66190     Total production time = 7.2809e-13 [s] ~ 1219 w_p T_prod     snapshot interval step = 10     snapshot interval time = 1.1000e-16 [s] = 0.1843 w_p T_snap    Total number of snapshots = 6619 Bernu ran a simulation for 244 plasma periods. We decided to run for a time 5 times longer. This is because we want to show the time slicing option in the Postprocessing phase.The estimated times, on this machine, are     ----------------------- Total Estimated Times ------------------------     Equilibration Time: 0 hrs 1 min 13 sec    Magnetization Time: 0 hrs 1 min 6 sec    Production Time: 0 hrs 33 min 35 sec    Total Run Time: 0 hrs 35 min 55 secand the entire run requires a minimum of 132 MB of space on the hard drive.     Total minimum needed space: 0 GB 151 MB 219 KB 26 bytesWe are ready to start the simulation.  Simulation<jupyter_code>sim = Simulation(input_file_name)
sim.setup(read_yaml=True)
sim.run()<jupyter_output>[38;2;24;69;49m






     _______.     ___      .______       __  ___      ___           _______.
    /       |    /   \     |   _  \     |  |/  /     /   \         /       |
   |   (----`   /  ^  \    |  |_)  |    |  '  /     /  ^  \       |   (----`
    \   \      /  /_\  \   |      /     |    <     /  /_\  \       \   \    
.----)   |    /  _____  \  |  |\  \----.|  .  \   /  _____  \  .----)   |   
|_______/    /__/     \__\ | _| `._____||__|\__\ /__/     \__\ |_______/    
                                                                            

[0m
An open-source pure-python molecular dynamics code for non-ideal plasmas.



* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
                                   Simulation                                   
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

Job ID:  ocp_mag
Job directory:  Simulations/ocp_mag
Equilibration dumps [...]<jupyter_text>The run completed in a time faster than the predicted one. Woohoo! PostProcessing<jupyter_code>postproc = PostProcess(input_file_name)
postproc.setup(read_yaml=True)

postproc.therm.setup(postproc.parameters)<jupyter_output>* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
                                 Postprocessing                                 
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

Job ID:  ocp_mag
Job directory:  Simulations/ocp_mag
Equilibration dumps directory:  Simulations/ocp_mag/Simulation/Equilibration/dumps
Production dumps directory:  Simulations/ocp_mag/Simulation/Production/dumps
PostProcessing directory:  Simulations/ocp_mag/PostProcessing<jupyter_text>Before proceeding to the calculation of the diffusion coefficient we check for conservation of energy during the production phase.<jupyter_code>postproc.therm.temp_energy_plot(postproc, 'production')<jupyter_output><empty_output><jupyter_text>The energy remained fairly constant over the entire production phase. The Temperature similarly oscillated widly ($\sim 10 \%$) but its average is exactly the desired value. Both the Temperature and the Total Energy show a Gaussian distribution around their means.We are ready to calculate the diffusion coefficient.  TransportTransport coefficient are calculated using the `TransportCoefficient` class in `sarkas.tools.transport`. Once the class is imported the diffusion coefficient is calculated by calling the method `diffusion`. This method requires the parameters of the simulation passed as `postproc.parameters`. The `phase` argument in this case is not needed since `production` is the default phase. As mentioned above we decided to run our simulation 4 times longer than Bernu's. This is because we wanted to divide our long simulation into 4 small simulations of length ~ 276 $\omega_p T$ (Bernu's was 244 $\omega_p T$), calculate the diffusion coefficients for each of them, and average them. Sarkas was created with this workflow in mind and the option ``no_slices`` allows to do exactly this.The ``diffusion`` method will calculate 4 velocity auto-correlation functions, their diffusion coefficients, and plot their averages with a shaded area.In the case of a magnetized plasma we have two diffusion coefficients:$$ D_{\parallel}(\tau) = \frac{k_B T}{m} \int_0^\tau dt Z_{\parallel}(t), \quad D_{\perp}(\tau) = \frac{k_B T}{m} \int_0^\tau dt Z_{\perp}(t)$$where, for the case of a constant magnetic field in the $z$-direction$$ Z_{\parallel}(t) = \langle v_z(t)v_z(0) \rangle, \quad Z_{\perp}(t) = \frac{1}{2} \left [ \langle v_x(t)v_x(0) \rangle + \langle v_y(t)v_y(0) \rangle \right ],$$are the velocity autocorrelation functions. These functions are calculated by the class `VelocityAutocorrelationFunction` in `sarkas.tools.observables`. This class calculates and saves $Z(t)$ for each axis and the total $Z(t) = \frac{1}{3} \left [ Z_x(t) + Z_y(t) + Z_z(t) \right ]$. The `diffusion` method takes care of dividing it into the perpendicular and parallel. It is important to note that the `diffusion` method assume a constant magnetic field in the $z$-direction. Arbitrary directions of the magnetic field are not supported at the moment.<jupyter_code>from sarkas.tools.transport import TransportCoefficient
diffusion = TransportCoefficient.diffusion(postproc.parameters, phase = 'production', no_slices = 5)<jupyter_output>======================= Diffusion Coefficient ======================== 

Data saved in: 
 Simulations/ocp_mag/PostProcessing/VelocityAutoCorrelationFunction/Production/Diffusion_ocp_mag.csv

No. of slices = 5
No. dumps per slice = 1324
Time interval of autocorrelation function = 1.4564e-13 [s] ~ 244 w_p T

Calculating vacf for slice 0/5.

Parsing particles' velocities.<jupyter_text>Once the calculation is complete `Sarkas` produces a figure containing two ``subplots``: a plot of $Z_{\perp}(t)$ and $Z_{\parallel}(t)$ on the left and a plot of $D_{\parallel}(\tau)$ and $D_{\perp}(\tau)$ on the right $y$-axis.The above plots . Therefore, we re-plot the diffusion only by rescaling it by $\omega_p a_{\rm ws}^2$.<jupyter_code>rescaling = postproc.parameters.total_plasma_frequency * postproc.parameters.a_ws**2
fig, ax = plt.subplots(1,1, figsize=(10,7))
ax2 = ax.twiny()
ax.plot(diffusion["Time"]*1e15, diffusion["H Parallel Diffusion avg"]/rescaling,
       label = r'$D_{\parallel}$')
ax.fill_between(
    diffusion["Time"]*1e15, 
    diffusion["H Parallel Diffusion avg"]/rescaling - diffusion["H Parallel Diffusion std"]/rescaling,
    diffusion["H Parallel Diffusion avg"]/rescaling + diffusion["H Parallel Diffusion std"]/rescaling,
    alpha = 0.2)

# Perpendicular
ax.plot(diffusion["Time"]*1e15, diffusion["H Perpendicular Diffusion avg"]/rescaling,
        label = r'$D_{\perp}$')
ax.fill_between(
    diffusion["Time"]*1e15, 
    diffusion["H Perpendicular Diffusion avg"]/rescaling - diffusion["H Perpendicular Diffusion std"]/rescaling,
    diffusion["H Perpendicular Diffusion avg"]/rescaling + diffusion["H Perpendicular Diffusion std"]/rescaling,
    alpha = 0.2)

ax.plot(diffusion["Time"]*1e15, diffusion["H Parallel Diffusion avg"].expanding().mean()/rescaling,
       ls = '--', label = 'Mov Avg')

ax.plot(diffusion["Time"]*1e15, diffusion["H Perpendicular Diffusion avg"].expanding().mean()/rescaling,
        ls = '--', label = r'Mov Avg')
# ax.axhline(coefficient["H Parallel Diffusion avg"].iloc[10]/rescaling)
# ax2.plot(diffusion["H Parallel Diffusion avg"].expanding().mean()/rescaling, ls = '')
# ax.grid(False)
ax.legend(ncol = 2)
ax.set_xlim(diffusion["Time"].iloc[1]*1e15, 1.5*diffusion["Time"].iloc[-1]*1e15)
_ = ax.set( xlabel = r'$\tau$ [fs]', ylabel =r'$D_{\parallel,\perp}$ [$\omega_p a_{ws}^2$]', xscale = 'log')

ax2.set_xlim([1,1500*1.5])
ax2.grid(alpha = 0.1)
ax2.set( xscale = 'log', xlabel = 'Index')<jupyter_output><empty_output><jupyter_text>This gives a clearer picture of the diffusion coefficient. Note that the diffusion is a constant and not a function of time. This plot is meant to show how the diffusion changes as the upper limit of the time integration is varied. The final value to be published would be an average over the last ~ 800 diffusion values.<jupyter_code>print('The perpendicular diffusion is D_x = {:.4f} w_p a_ws^2'.format(
    1/rescaling * diffusion["H Perpendicular Diffusion avg"].iloc[-800:].mean()) ) 
print('The parallel diffusion is D_z = {:.4f} w_p a_ws^2'.format(
    1/rescaling * diffusion["H Parallel Diffusion avg"].iloc[-800:].mean()) )<jupyter_output>The perpendicular diffusion is D_x = 0.0024 w_p a_ws^2
The parallel diffusion is D_z = 0.0026 w_p a_ws^2
{'fraction_non_alphanumeric': 0.1352482560525236, 'fraction_numerical': 0.02256873204759951, 'mean_word_length': 2.946243523316062, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 27, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 8, 'word_list_counts': {'cursed_substrings.json': 61, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '616305', 'license': 'MIT', 'n_tokens_mistral': 4087, 'n_tokens_neox': 3799, 'n_words': 1137, 'path': 'docs/examples/Magnetized_OCP/Magnetized_Plasma.ipynb', 'repo_name': 'pwessels-uhh/sarkas', 'size': '12185'}