Skip to content

Commit

Permalink
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Small improvements to documentation for thermostats. Bug fixes for pa…
Browse files Browse the repository at this point in the history
…rticles initialization.
lucianogsilvestri committed May 11, 2024
1 parent 2b0b974 commit 5262f36
Showing 13 changed files with 1,502 additions and 221 deletions.
380 changes: 274 additions & 106 deletions docs/documentation/Features_files/Berendsen_NB/Berendsen_Thermostat.ipynb

Large diffs are not rendered by default.

100 changes: 100 additions & 0 deletions docs/documentation/Features_files/Berendsen_NB/berendsen_script.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,100 @@
# Import the usual libraries
import numpy as np
import matplotlib.pyplot as plt

import os
plt.style.use('MSUstyle')

from multiprocessing import Process

# Import sarkas
from sarkas.processes import PreProcess, Simulation

def run_simulation_berendsen(input_file_name, tau, eq_steps, cycles):

args = {
'Integrator': {'thermalization_timestep': 0, # Timesteps before turning on the thermostat
'berendsen_tau': tau}, # Change tau for each simulation
"Parameters":{ "rand_seed": 123456,
"equilibration_steps": eq_steps},
"IO": # Store all simulations' data in simulations_dir,
# but save the dumps in different subfolders (job_dir)
{
"job_dir": f"tau_cycles{cycles:.2f}",
"verbose": False # This is so not to print to screen for every run
},
}

sim = Simulation(input_file_name)
sim.setup(read_yaml=True, other_inputs=args)
sim.run()

print(f'tau cycle = {cycles:.2f} Done')

def run_simulation_langevin(input_file_name, gamma, eq_steps, cycles):

args = {
'Integrator': {'langevin_gamma': gamma}, # Change tau for each simulation
"Parameters":{ "rand_seed": 123456,
"equilibration_steps": eq_steps},
"IO": # Store all simulations' data in simulations_dir,
# but save the dumps in different subfolders (job_dir)
{
"job_dir": f"gamma_cycles{cycles:.2f}",
"verbose": False # This is so not to print to screen for every run
},
}

sim = Simulation(input_file_name)
sim.setup(read_yaml=True, other_inputs=args)
sim.run()

print(f'gamma cycle = {cycles:.2f} Done')


if __name__ == '__main__':
# Create the file path to the YAML input file
input_file_name = os.path.join('input_files', 'yocp_quickstart.yaml' )

cycles = np.array([0.01, 0.1, 0.2, 1.0, 5.0, 10.0])

pre = PreProcess(input_file_name)
pre.setup(read_yaml=True)

tau_p = 2.0 * np.pi /pre.parameters.total_plasma_frequency

taus = - tau_p / np.log(0.01) * cycles / pre.parameters.dt
eq_steps = np.rint(2.0 * cycles * tau_p/ pre.parameters.dt).astype(int)

processes = []
# Create the file path to the YAML input file
input_file_name = os.path.join('input_files', 'yocp_quickstart.yaml' )

for i, tau in enumerate(taus):

p0 = Process(target = run_simulation_berendsen, args = (input_file_name, tau, eq_steps[i], cycles[i],) )
processes.append(p0)
p0.start()

# Wait for all processes to finish
for p in processes:
p.join()

processes = []
# Create the file path to the YAML input file
input_file_name = os.path.join('input_files', 'yocp_quickstart_langevin.yaml' )

pre = PreProcess(input_file_name)
pre.setup(read_yaml=True)

tau_p = 2.0 * np.pi /pre.parameters.total_plasma_frequency

gammas = - np.log(0.01)/(2.0 * tau_p * cycles)
for i, gamma in enumerate(gammas):
p0 = Process(target = run_simulation_langevin, args = (input_file_name, gamma, eq_steps[i], cycles[i],) )
processes.append(p0)
p0.start()

# Wait for all processes to finish
for p in processes:
p.join()
Original file line number Diff line number Diff line change
@@ -6,7 +6,7 @@ Particles:
mass: 2.0089e-23 # g
num: 1000 # total number of particles
Z: 1.976599 # degree of ionization
temperature: 5.e+03 # T = 0.5 eV
temperature_eV: 0.5 # T = 0.5 eV

Potential:
type: Yukawa # potential type
@@ -16,29 +16,23 @@ Potential:
Integrator:
type: Verlet
dt: 5.0e-17 # sec
equilibration_steps: 5000 # number of timesteps for the equilibrium
production_steps: 0 # number of timesteps after the equilibrium
eq_dump_step: 2
prod_dump_step: 10
boundary_conditions: periodic # REQUIRED
thermalization: yes # OPTIONAL. Default = yes
thermostat_type: berendsen # REQUIRED if thermalization is yes
thermalization_timestep: 50 # REQUIRED if thermalization is yes
berendsen_tau: 1.0 # REQUIRED if thermostat: berendsen
thermostate_temperatures_eV: 0.5 # OPTIONAL Default = Species.temperature_eV
thermalization_timestep: 500 # REQUIRED if thermalization is yes
berendsen_tau: 2.0 # REQUIRED if thermostat: berendsen

Parameters:
units: cgs # units
load_method: random_no_reject
equilibration_steps: 5000 # REQUIRED
production_steps: 0 # REQUIRED
equilibration_steps: 1500 # REQUIRED
production_steps: 2500 # REQUIRED
eq_dump_step: 2 # REQUIRED
prod_dump_step: 10 # REQUIRED

IO:
verbose: yes
md_simulations_dir: Berendsen_runs # dir name to save data. The default is "Checkpoint"
job_dir: tau_2
verbose: False
job_dir: berendsen_test

Observables:
- Thermodynamics:
Original file line number Diff line number Diff line change
@@ -0,0 +1,42 @@
# keywords: YOCP, Yukawa, PP, mks, random_no_reject
Particles:
- Species:
name: C
number_density: 1.13666931822e+23 # N/cc
mass: 2.0089e-23 # g
num: 1000 # total number of particles
Z: 1.976599 # degree of ionization
temperature_eV: 0.5 # T = 0.5 eV

Potential:
type: Yukawa # potential type
method: PP
rc: 6.0e-8 # cm, cut-off radius

Integrator:
equilibration_type: langevin
production_type: verlet
dt: 5.0e-17 # sec
boundary_conditions: periodic # REQUIRED
thermalization: no # OPTIONAL. Default = yes
langevin_type: BBK
langevin_gamma: 9.2774e+13 # [Hertz]

Parameters:
units: cgs # units
load_method: random_no_reject
equilibration_steps: 1500 # REQUIRED
production_steps: 2500 # REQUIRED
eq_dump_step: 2 # REQUIRED
prod_dump_step: 10 # REQUIRED

IO:
verbose: False
job_dir: langevin_test

Observables:
- Thermodynamics:
phase: production

- RadialDistributionFunction:
no_bins: 250
8 changes: 7 additions & 1 deletion docs/documentation/why_sarkas.rst
Original file line number Diff line number Diff line change
@@ -155,9 +155,15 @@ screening parameters and measure their diffusion coefficient. An example script
rng = np.random.default_rng()
# Create arrays of screening parameters
kappas = np.linspace(1, 5)
# Run 10 simulations
# Run 10 parallel simulations
processes = []
for i, kappa in enumerate(kappas):
p0 = Process(target = launch, args = ( rng.integers(low = 2**32), kappa) )
processes.append(p0)
p0.start()
# Wait for all simulations to finish
for p in processes:
p.join()
Notice how both the simulation and the postprocessing can be done all in one script.
Loading

0 comments on commit 5262f36

Please sign in to comment.