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.

Downloads: Notebook | Script ?

Requires: AMS2027 or later

../_images/diffusion_supercell_dependence_7_1_45965171.png

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);

Ar atoms, density = 1.0 g/cm^3
Ar atoms, density = 1.0 g/cm^3
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

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")

Diffusion Coefficient Supercell Dependence¶

Diffusion coefficient supercell dependence¶

See also¶

Related examples¶

Related tutorials¶

Related documentation¶

Python Script¶