Diffusion Coefficient Supercell Dependence¶
Run Lennard-Jones argon simulations for several box sizes, calculate diffusion coefficients from the mean-squared displacement, and extrapolate the results to infinite system size with a linear fit in 1/L.
Diffusion coefficient supercell dependence¶
Finite-size effects make the diffusion coefficient depend on the size of the simulation cell. This example runs Lennard-Jones argon boxes of different sizes, computes the diffusion coefficient from the mean-squared displacement, and linearly extrapolates the results to the infinite-system limit.
from scm.inputs import AMS, LennardJones
from scm.plams import Settings
def get_lennard_jones_settings():
s = Settings()
s.input = AMS(Engine=LennardJones(Eps=3.785e-4, RMin=3.81637, Cutoff=8.0))
s.runscript.nproc = 1
return s
Set up three argon boxes at the same density and visualize the resulting supercells.
from scm.base import ChemicalSystem
from scm.plams import packmol, preoptimize, plot_image_grid, view
sizes = [50, 150, 450]
inverted_L_list = []
density = 1.0 # g/cm^3
Ar = ChemicalSystem()
Ar.add_atom("Ar", coords=(0, 0, 0))
boxes = []
imgs = {}
for size in sizes:
mol = packmol(Ar, n_molecules=size, density=density)
mol = preoptimize(mol, settings=get_lennard_jones_settings(), maxiterations=50)
boxes.append(mol)
L = mol.lattice[0][0]
inverted_L_list.append(1.0 / L)
print(f"{size} Ar atoms, density = {density} g/cm^3")
imgs[f"{size} Ar atoms"] = view(mol, height=200, width=200, direction="tilt_z")
plot_image_grid(imgs, rows=1);
50 Ar atoms, density = 1.0 g/cm^3
150 Ar atoms, density = 1.0 g/cm^3
450 Ar atoms, density = 1.0 g/cm^3
For each system, run an NVT equilibration, an NVT production simulation, and an MSD analysis. The diffusion coefficients are then fitted as a function of 1/L.
D_list = []
temperature = 400
max_correlation_time_fs = 6000
start_time_fit_fs = 3000
from scm.plams import AMSNVTJob, AMSMSDJob, log
for size, mol in zip(sizes, boxes):
eq_job = AMSNVTJob(
name=f"eq_{size}",
settings=get_lennard_jones_settings(),
molecule=mol,
temperature=temperature,
nsteps=10000,
samplingfreq=500,
thermostat="Berendsen",
tau=100,
timestep=1.0,
writebonds=False,
writecharges=False,
writemolecules=False,
writevelocities=False,
)
eq_job.run()
prod_job = AMSNVTJob.restart_from(
eq_job,
name=f"prod_{size}",
nsteps=60000,
samplingfreq=25,
temperature=temperature,
thermostat="NHC",
)
prod_job.run()
msd_job = AMSMSDJob(
prod_job, name=f"msd_{size}", max_correlation_time_fs=max_correlation_time_fs
)
msd_job.run()
D = msd_job.results.get_diffusion_coefficient(start_time_fit_fs=start_time_fit_fs)
D_list.append(D)
log(f"{size} atoms: D = {D:.2e} m^2 s^-1")
[16.10|09:25:24] JOB eq_50 STARTED
[16.10|09:25:24] JOB eq_50 RUNNING
[16.10|09:25:26] JOB eq_50 FINISHED
[16.10|09:25:26] JOB eq_50 SUCCESSFUL
[16.10|09:25:27] JOB prod_50 STARTED
[16.10|09:25:27] JOB prod_50 RUNNING
[16.10|09:25:38] JOB prod_50 FINISHED
[16.10|09:25:38] JOB prod_50 SUCCESSFUL
[16.10|09:25:38] JOB msd_50 STARTED
[16.10|09:25:38] JOB msd_50 RUNNING
... output trimmed ....
[16.10|09:26:07] JOB eq_450 SUCCESSFUL
[16.10|09:26:07] JOB prod_450 STARTED
[16.10|09:26:07] JOB prod_450 RUNNING
[16.10|09:27:19] JOB prod_450 FINISHED
[16.10|09:27:19] JOB prod_450 SUCCESSFUL
[16.10|09:27:19] JOB msd_450 STARTED
[16.10|09:27:19] JOB msd_450 RUNNING
[16.10|09:27:20] JOB msd_450 FINISHED
[16.10|09:27:21] JOB msd_450 SUCCESSFUL
[16.10|09:27:21] 450 atoms: D = 2.40e-08 m^2 s^-1
import matplotlib.pyplot as plt
from scm.plams import log
from scm.plams.tools.plot import linear_fit_extrapolate_to_0
fit_x, fit_y, slope, intercept = linear_fit_extrapolate_to_0(inverted_L_list, D_list)
dilute_D = intercept
title = f"Extrapolated (dilute) D: {dilute_D:.2e} m²/s"
log(title)
fig, ax = plt.subplots(figsize=(6, 3))
ax.plot(inverted_L_list, D_list, ".", label="Simulation results")
ax.plot(fit_x, fit_y, label="Linear fit")
ax.set_title(title)
ax.legend()
ax.set_xlabel("1/L (Å⁻¹)")
ax.set_ylabel("Diffusion coefficient D (m²/s)")
ax;
[16.10|09:27:21] Extrapolated (dilute) D: 2.61e-08 m²/s
See also¶
Python Script¶
#!/usr/bin/env python
# coding: utf-8
# ## Diffusion coefficient supercell dependence
#
# Finite-size effects make the diffusion coefficient depend on the size of the simulation cell. This example runs Lennard-Jones argon boxes of different sizes, computes the diffusion coefficient from the mean-squared displacement, and linearly extrapolates the results to the infinite-system limit.
#
from scm.inputs import AMS, LennardJones
from scm.plams import Settings
def get_lennard_jones_settings():
s = Settings()
s.input = AMS(Engine=LennardJones(Eps=3.785e-4, RMin=3.81637, Cutoff=8.0))
s.runscript.nproc = 1
return s
# Set up three argon boxes at the same density and visualize the resulting supercells.
from scm.base import ChemicalSystem
from scm.plams import packmol, preoptimize, plot_image_grid, view
sizes = [50, 150, 450]
inverted_L_list = []
density = 1.0 # g/cm^3
Ar = ChemicalSystem()
Ar.add_atom("Ar", coords=(0, 0, 0))
boxes = []
imgs = {}
for size in sizes:
mol = packmol(Ar, n_molecules=size, density=density)
mol = preoptimize(mol, settings=get_lennard_jones_settings(), maxiterations=50)
boxes.append(mol)
L = mol.lattice[0][0]
inverted_L_list.append(1.0 / L)
print(f"{size} Ar atoms, density = {density} g/cm^3")
imgs[f"{size} Ar atoms"] = view(mol, height=200, width=200, direction="tilt_z")
plot_image_grid(imgs, rows=1, save_path="picture1.png");
# For each system, run an NVT equilibration, an NVT production simulation, and an MSD analysis. The diffusion coefficients are then fitted as a function of ``1/L``.
#
D_list = []
temperature = 400
max_correlation_time_fs = 6000
start_time_fit_fs = 3000
from scm.plams import AMSNVTJob, AMSMSDJob, log
for size, mol in zip(sizes, boxes):
eq_job = AMSNVTJob(
name=f"eq_{size}",
settings=get_lennard_jones_settings(),
molecule=mol,
temperature=temperature,
nsteps=10000,
samplingfreq=500,
thermostat="Berendsen",
tau=100,
timestep=1.0,
writebonds=False,
writecharges=False,
writemolecules=False,
writevelocities=False,
)
eq_job.run()
prod_job = AMSNVTJob.restart_from(
eq_job,
name=f"prod_{size}",
nsteps=60000,
samplingfreq=25,
temperature=temperature,
thermostat="NHC",
)
prod_job.run()
msd_job = AMSMSDJob(
prod_job, name=f"msd_{size}", max_correlation_time_fs=max_correlation_time_fs
)
msd_job.run()
D = msd_job.results.get_diffusion_coefficient(start_time_fit_fs=start_time_fit_fs)
D_list.append(D)
log(f"{size} atoms: D = {D:.2e} m^2 s^-1")
import matplotlib.pyplot as plt
from scm.plams import log
from scm.plams.tools.plot import linear_fit_extrapolate_to_0
fit_x, fit_y, slope, intercept = linear_fit_extrapolate_to_0(inverted_L_list, D_list)
dilute_D = intercept
title = f"Extrapolated (dilute) D: {dilute_D:.2e} m²/s"
log(title)
fig, ax = plt.subplots(figsize=(6, 3))
ax.plot(inverted_L_list, D_list, ".", label="Simulation results")
ax.plot(fit_x, fit_y, label="Linear fit")
ax.set_title(title)
ax.legend()
ax.set_xlabel("1/L (Å⁻¹)")
ax.set_ylabel("Diffusion coefficient D (m²/s)")
ax;
ax.figure.savefig("picture2.png")