Moliere Potential

Simulation data from Figure 8 of Vorberger et al

The YAML input file can be found at input_file and this notebook at notebook.

# Import the usual packages.
%pylab
%matplotlib inline
import os
import pandas as pd

from sarkas.processes import PreProcess, Simulation, PostProcess

# Set the plotting style
plt.style.use('MSUstyle')

# Link to the input file.
input_file = os.path.join('input_files', 'moliere_mks.yaml')

Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib

Simulation

sim = Simulation(input_file)
sim.setup(read_yaml = True)
sim.run()

 __            _
/ _\ __ _ _ __| | ____ _ ___
\ \ / _` | '__| |/ / _` / __|
_\ \ (_| | |  |   < (_| \__ \
\__/\__,_|_|  |_|\_\__,_|___/



An open-source pure-python molecular dynamics suite for non-ideal plasmas.





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

Job ID:  Moliere_mks
Job directory:  Simulations/Moliere_mks

Equilibration dumps directory:
 Simulations/Moliere_mks/Simulation/Equilibration/dumps
Production dumps directory:
 Simulations/Moliere_mks/Simulation/Production/dumps

Equilibration Thermodynamics file:
 Simulations/Moliere_mks/Simulation/Equilibration/EquilibrationEnergy_Moliere_mks.csv
Production Thermodynamics file:
 Simulations/Moliere_mks/Simulation/Production/ProductionEnergy_Moliere_mks.csv
Random Seed =  654984647

PARTICLES:
Total No. of particles =  10000
No. of species =  1
Species ID: 0
        Name: Be
        No. of particles = 10000
        Number density = 1.234875e+29 [N/m^3]
        Atomic weight = 9.0122 [a.u.]
        Mass = 1.507397e-26 [kg]
        Mass density = 1.848000e+06 [kg/m^3]
        Charge number/ionization degree = 4.0000
        Charge = 6.408707e-19 [C]
        Temperature = 1.392542e+05 [K] = 1.200000e+01 [eV]
        Debye Length = 1.832054e-11 [1/m^3]
        Plasma Frequency = 6.164440e+14 [rad/s]

SIMULATION AND INITIAL PARTICLE BOX:
Units:  mks
Wigner-Seitz radius = 1.245746e-10 [m]
No. of non-zero box dimensions =  3
Box side along x axis = 3.472931e+01 a_ws = 4.326390e-09 [m]
Box side along y axis = 3.472931e+01 a_ws = 4.326390e-09 [m]
Box side along z axis = 3.472931e+01 a_ws = 4.326390e-09 [m]
Box Volume = 8.097987e-26 [m^3]
Initial particle box side along x axis = 3.472931e+01 a_ws = 4.326390e-09 [m]
Initial particle box side along y axis = 3.472931e+01 a_ws = 4.326390e-09 [m]
Initial particle box side along z axis = 3.472931e+01 a_ws = 4.326390e-09 [m]
Initial particle box Volume = 8.097987e-26 [m^3]
Boundary conditions: periodic

ELECTRON PROPERTIES:
Number density: n_e = 4.939499e+29 [N/m^3]
Wigner-Seitz radius: a_e = 7.847708e-11 [m]
Temperature: T_e = 1.392542e+05 [K] = 1.200000e+01 [eV]
de Broglie wavelength: lambda_deB = 1.997449e-10 [m]
Thomas-Fermi length: lambda_TF = 4.721828e-11 [m]
Fermi wave number: k_F = 2.445502e+10 [1/m]
Fermi Energy: E_F = 2.278552e+01 [eV]
Relativistic parameter: x_F = 9.443531e-03 --> E_F = 2.278501e+01 [eV]
Degeneracy parameter: Theta = 5.266503e-01
Coupling: r_s = 1.483002,  Gamma_e = 1.529071
Warm Dense Matter Parameter: W = 7.5544e-01
Chemical potential: mu = 1.3549e+00 k_B T_e = 7.1357e-01 E_F

POTENTIAL:  moliere

ALGORITHM:  pp
rcut = 5.3219 a_ws = 6.629755e-10 [m]
No. of PP cells per dimension =  6,  6,  6
No. of particles in PP loop =    971
No. of PP neighbors per particle =    150
Tot Force Error = 5.713963e-20

THERMOSTAT:
Type: berendsen
First thermostating timestep, i.e. relaxation_timestep = 50
Berendsen parameter tau: 2.000 [timesteps]
Berendsen relaxation rate: 0.500 [1/timesteps]
Thermostating temperatures:
Species ID 0: T_eq = 1.392542e+05 [K] = 1.200000e+01 [eV]

INTEGRATOR:
Type: verlet
Time step = 3.244415e-17 [s]
Total plasma frequency = 6.164440e+14 [rad/s]
w_p dt = 0.0200 ~ 1/49

Equilibration:
No. of equilibration steps = 1000
Total equilibration time = 3.2444e-14 [s] ~ 20 w_p T_eq
snapshot interval step = 10
snapshot interval time = 3.2444e-16 [s] = 0.2000 w_p T_snap
Total number of snapshots = 100

Production:
No. of production steps = 5000
Total production time = 1.6222e-13 [s] ~ 100 w_p T_prod
snapshot interval step = 10
snapshot interval time = 3.2444e-16 [s] = 0.2000 w_p T_snap
Total number of snapshots = 500


------------------------ Initialization Times ------------------------


Potential Initialization Time: 1 sec 171 msec 420 usec 854 nsec

Particles Initialization Time: 0 sec 1 msec 914 usec 956 nsec

Total Simulation Initialization Time: 1 sec 171 msec 420 usec 854 nsec


--------------------------- Equilibration ----------------------------

Equilibration Time: 0 hrs 3 min 35 sec


----------------------------- Production -----------------------------

Production Time: 0 hrs 24 min 6 sec

Total Time: 0 hrs 27 min 43 sec

Post Processing

postproc = PostProcess(input_file)
postproc.setup(read_yaml = True)
postproc.run()

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
                                 Postprocessing
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

Job ID:  Moliere_mks
Job directory:  Simulations/Moliere_mks
PostProcessing directory:
 Simulations/Moliere_mks/PostProcessing

Equilibration dumps directory:  Simulations/Moliere_mks/Simulation/Equilibration/dumps
Production dumps directory:
 Simulations/Moliere_mks/Simulation/Production/dumps

Equilibration Thermodynamics file:
 Simulations/Moliere_mks/Simulation/Equilibration/EquilibrationEnergy_Moliere_mks.csv
Production Thermodynamics file:
 Simulations/Moliere_mks/Simulation/Production/ProductionEnergy_Moliere_mks.csv


==================== Radial Distribution Function ====================

Data saved in:
 Simulations/Moliere_mks/PostProcessing/RadialDistributionFunction/Production/RadialDistributionFunction_Moliere_mks.h5
Data accessible at: self.ra_values, self.dataframe

No. bins = 500
dr = 0.0106 a_ws = 1.3260e-12 [m]
Maximum Distance (i.e. potential.rc)= 5.3219 a_ws = 6.6298e-10 [m]

5000

Radial Distribution Function Calculation Time: 0 sec 82 msec 224 usec 470 nsec

postproc.rdf.plot(scaling= postproc.rdf.a_ws,
                 xlabel = r'$r/a_{ws}$',
                 title = 'Pair Distribution Function')

<AxesSubplot:title={'center':'Pair Distribution Function'}, xlabel='$r/a_{ws}$'>

from sarkas.potentials import moliere as mol
from sarkas.potentials import coulomb as coul

r_array = postproc.parameters.box_lengths[0] * np.linspace(0.001, 1, 1000)
mol_potential = np.zeros(len(r_array))
coul_potential = np.zeros(len(r_array))

# I need to do this, because Numba compiled the function to take in a float and an array not two arrays
for (ir, r) in enumerate(r_array):
    mol_potential[ir], _ = mol.moliere_force(r, postproc.potential.matrix[:,0, 0])
    coul_potential[ir], _ = coul.coulomb_force(r, postproc.potential.matrix[:,0, 0])

pot_const = postproc.potential.matrix[0,0,0]
plt.plot(r_array/postproc.parameters.a_ws, mol_potential/pot_const, label = 'Moliere')
plt.plot(r_array/postproc.parameters.a_ws, coul_potential/pot_const, label = 'Coulomb')
plt.legend()
plt.yscale('log')
plt.ylabel(r'$U(r)$   [$Q^2/(4\pi \epsilon_0)$]')
plt.xlabel(r'$r/a_{\rm ws}$')
plt.xlim(0,10)
plt.ylim(1e5,1e11)
plt.title('Potential Plot')

Text(0.5, 1.0, 'Potential Plot')