Diffusion Coefficient Temperature Dependence
********************************************

Molecular Dynamics Diffusion Lennard-Jones MSD Parallel jobs

Run Lennard-Jones argon simulations at several temperatures in
parallel, extract diffusion coefficients from the mean-squared
displacement, and fit an Arrhenius relation to estimate the prefactor
and activation barrier.

**Requires:** AMS2026 or later

-[ Related examples ]-

* Molecular Dynamics with Python

* Diffusion Coefficient Supercell Dependence

* ReaxFF Density Equilibration

* Basic MD Trajectory Analysis with PLAMS

* Viscosity from MD with the Green-Kubo Formalism

-[ Related tutorials ]-

* Ionic conductivity tutorial

-[ Related documentation ]-

* Molecular Dynamics task

* Lennard-Jones

* Trajectory analysis utilities

[image]


Diffusion coefficient temperature dependence
============================================

The diffusion coefficient often follows an Arrhenius relation with
temperature. This example runs Lennard-Jones argon simulations at
several temperatures in parallel, extracts the diffusion coefficient
from each trajectory, and fits "ln(D)" versus "1/T" to estimate the
prefactor and activation barrier.

   from scm.plams import Settings


   def get_lennard_jones_settings():
       s = Settings()
       s.input.lennardjones.eps = 3.785e-4
       s.input.lennardjones.rmin = 3.81637
       s.input.lennardjones.cutoff = 8.0
       s.runscript.nproc = 1
       return s

Use a 50-atom argon box and evaluate the diffusion coefficient at 50,
100, 200, and 300 K.

   from scm.plams import Molecule, Atom, packmol, preoptimize, view

   T_list = [50, 100, 200, 300]
   density = 1.0  # g/cm^3
   max_correlation_time_fs = 6000
   start_time_fit_fs = 3000

   Ar = Molecule()
   Ar.add_atom(Atom(symbol="Ar", coords=(0, 0, 0)))

   mol = packmol(Ar, n_molecules=50, density=density)
   mol = preoptimize(mol, settings=get_lennard_jones_settings(), maxiterations=50)
   view(mol, height=200, width=200, direction="tilt_z", padding=-2)

   [image: image generated from notebook][image]

If you have multiple cores available, a parallel "JobRunner" can start
the trajectories concurrently.

   from scm.plams import config, JobRunner

   maxjobs = 4
   config.default_jobrunner = JobRunner(parallel=True, maxjobs=maxjobs)

To keep PLAMS from blocking too early, the production jobs use
"use_prerun=True" and the MSD results are only accessed after all jobs
have been submitted.

   from scm.plams import AMSNVTJob, AMSMSDJob, log

   msd_jobs = []
   DT_list = []
   for T in T_list:
       eq_job = AMSNVTJob(
           name=f"eq_T_{T}",
           settings=get_lennard_jones_settings(),
           molecule=mol,
           temperature=T,
           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_T_{T}",
           nsteps=50000,
           samplingfreq=25,
           temperature=T,
           thermostat="NHC",
           use_prerun=True,
       )
       prod_job.run()
       msd_job = AMSMSDJob(
           prod_job, name=f"msd_T_{T}", max_correlation_time_fs=max_correlation_time_fs
       )
       msd_job.run()
       msd_jobs.append(msd_job)


   for T, msd_job in zip(T_list, msd_jobs):
       D = msd_job.results.get_diffusion_coefficient(start_time_fit_fs=start_time_fit_fs)
       DT_list.append(D)
       log(f"Temperature {T} K: D = {D:.2e} m^2 s^-1")

   [23.03|14:56:29] JOB eq_T_50 STARTED
   [23.03|14:56:29] JOB prod_T_50 STARTED
   [23.03|14:56:29] JOB msd_T_50 STARTED
   [23.03|14:56:29] JOB eq_T_100 STARTED
   [23.03|14:56:29] Waiting for job eq_T_50 to finish
   [23.03|14:56:29] JOB prod_T_100 STARTED
   [23.03|14:56:29] JOB msd_T_100 STARTED
   [23.03|14:56:29] JOB eq_T_200 STARTED
   [23.03|14:56:29] JOB eq_T_50 RUNNING
   [23.03|14:56:29] Waiting for job prod_T_50 to finish
   ... output trimmed ....
   [23.03|14:56:44] JOB msd_T_200 FINISHED
   [23.03|14:56:44] JOB msd_T_50 SUCCESSFUL
   [23.03|14:56:44] JOB msd_T_100 SUCCESSFUL
   [23.03|14:56:44] Temperature 50 K: D = 1.73e-09 m^2 s^-1
   [23.03|14:56:44] Temperature 100 K: D = 6.03e-09 m^2 s^-1
   [23.03|14:56:44] JOB msd_T_200 SUCCESSFUL
   [23.03|14:56:44] JOB msd_T_300 SUCCESSFUL
   [23.03|14:56:44] Waiting for job msd_T_200 to finish
   [23.03|14:56:44] Temperature 200 K: D = 1.19e-08 m^2 s^-1
   [23.03|14:56:44] Temperature 300 K: D = 1.75e-08 m^2 s^-1

   import numpy as np
   import matplotlib.pyplot as plt
   from scm.plams.tools.plot import linear_fit_extrapolate_to_0


   inverted_T_list = 1.0 / np.array(T_list)
   ln_D_list = np.log(DT_list)
   fit_x, fit_y, slope, intercept = linear_fit_extrapolate_to_0(inverted_T_list, ln_D_list)

   D0 = np.exp(intercept)
   boltzmann_eV_per_K = 8.617e-5
   barrier = -boltzmann_eV_per_K * slope
   title = f"D0: {D0:.2e} m²/s. Barrier: {barrier:.3f} eV"
   log(title)

   fig, ax = plt.subplots(figsize=(6, 3))
   ax.plot(inverted_T_list * 1000, ln_D_list, ".", label="Simulation results")
   ax.plot(fit_x * 1000, fit_y, label="Linear fit")
   ax.set_title(title)
   ax.legend()
   ax.set_xlabel("1000/T (K⁻¹)")
   ax.set_ylabel("ln D")
   ax;

   [23.03|14:56:44] D0: 2.48e-08 m²/s. Barrier: 0.012 eV

   [image: image generated from notebook][image]


See also
========


Related examples
----------------

* Molecular Dynamics with Python

* Diffusion Coefficient Supercell Dependence

* ReaxFF Density Equilibration

* Basic MD Trajectory Analysis with PLAMS

* Viscosity from MD with the Green-Kubo Formalism


Related tutorials
-----------------

* Ionic conductivity tutorial


Related documentation
---------------------

* Molecular Dynamics task

* Lennard-Jones

* Trajectory analysis utilities


Python Script
=============

   #!/usr/bin/env python
   # coding: utf-8

   # ## Diffusion coefficient temperature dependence
   #
   # The diffusion coefficient often follows an Arrhenius relation with temperature. This example runs Lennard-Jones argon simulations at several temperatures in parallel, extracts the diffusion coefficient from each trajectory, and fits ``ln(D)`` versus ``1/T`` to estimate the prefactor and activation barrier.
   #

   from scm.plams import Settings


   def get_lennard_jones_settings():
       s = Settings()
       s.input.lennardjones.eps = 3.785e-4
       s.input.lennardjones.rmin = 3.81637
       s.input.lennardjones.cutoff = 8.0
       s.runscript.nproc = 1
       return s


   # Use a 50-atom argon box and evaluate the diffusion coefficient at 50, 100, 200, and 300 K.
   #

   from scm.plams import Molecule, Atom, packmol, preoptimize, view

   T_list = [50, 100, 200, 300]
   density = 1.0  # g/cm^3
   max_correlation_time_fs = 6000
   start_time_fit_fs = 3000

   Ar = Molecule()
   Ar.add_atom(Atom(symbol="Ar", coords=(0, 0, 0)))

   mol = packmol(Ar, n_molecules=50, density=density)
   mol = preoptimize(mol, settings=get_lennard_jones_settings(), maxiterations=50)
   view(mol, height=200, width=200, direction="tilt_z", padding=-2, picture_path="picture1.png")


   # If you have multiple cores available, a parallel ``JobRunner`` can start the trajectories concurrently.
   #

   from scm.plams import config, JobRunner

   maxjobs = 4
   config.default_jobrunner = JobRunner(parallel=True, maxjobs=maxjobs)


   # To keep PLAMS from blocking too early, the production jobs use ``use_prerun=True`` and the MSD results are only accessed after all jobs have been submitted.
   #

   from scm.plams import AMSNVTJob, AMSMSDJob, log

   msd_jobs = []
   DT_list = []
   for T in T_list:
       eq_job = AMSNVTJob(
           name=f"eq_T_{T}",
           settings=get_lennard_jones_settings(),
           molecule=mol,
           temperature=T,
           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_T_{T}",
           nsteps=50000,
           samplingfreq=25,
           temperature=T,
           thermostat="NHC",
           use_prerun=True,
       )
       prod_job.run()
       msd_job = AMSMSDJob(prod_job, name=f"msd_T_{T}", max_correlation_time_fs=max_correlation_time_fs)
       msd_job.run()
       msd_jobs.append(msd_job)


   for T, msd_job in zip(T_list, msd_jobs):
       D = msd_job.results.get_diffusion_coefficient(start_time_fit_fs=start_time_fit_fs)
       DT_list.append(D)
       log(f"Temperature {T} K: D = {D:.2e} m^2 s^-1")


   import numpy as np
   import matplotlib.pyplot as plt
   from scm.plams.tools.plot import linear_fit_extrapolate_to_0


   inverted_T_list = 1.0 / np.array(T_list)
   ln_D_list = np.log(DT_list)
   fit_x, fit_y, slope, intercept = linear_fit_extrapolate_to_0(inverted_T_list, ln_D_list)

   D0 = np.exp(intercept)
   boltzmann_eV_per_K = 8.617e-5
   barrier = -boltzmann_eV_per_K * slope
   title = f"D0: {D0:.2e} m²/s. Barrier: {barrier:.3f} eV"
   log(title)

   fig, ax = plt.subplots(figsize=(6, 3))
   ax.plot(inverted_T_list * 1000, ln_D_list, ".", label="Simulation results")
   ax.plot(fit_x * 1000, fit_y, label="Linear fit")
   ax.set_title(title)
   ax.legend()
   ax.set_xlabel("1000/T (K⁻¹)")
   ax.set_ylabel("ln D")
   ax
   ax.figure.savefig("picture2.png")
