Basic molecular dynamics analysis

Note: This example requires AMS2023 or later.

This example illustrates how to calculate the basic

  • velocity autocorrelation function (VACF)

  • diffusion coefficient from the integral of the VACF

  • power spectrum (Fourier transform of the VACF)

  • viscosity from the Green-Kubo relation (integral of the off-diagonal pressure tensor autocorrelation function)

  • density along the z axis

For details about the functions, see the AMSResults API.

Technical

The example only shows how to technically calculate the numbers. For real simulations, run longer MD simulations and carefully converge any calculated quantities.

The viscosity requires especially long MD simulations.

Note

More advanced analysis is possible by setting up an AMSAnalysisJob job.

See also: Molecular Dynamics with Python

To follow along, either

Worked Example

Initial Imports

import os
import numpy as np
import matplotlib.pyplot as plt
from scm.plams import init, packmol, Settings, AMSJob, from_smiles, plot_molecule

Run Simple MD Simulation of Water

Run a short MD simulation of 16 water molecules in a box.

# this line is not required in AMS2025+
init()

mol = packmol(from_smiles("O"), n_molecules=16, density=1.0)
s = Settings()
s.input.ams.Task = "MolecularDynamics"
s.input.ReaxFF.ForceField = "Water2017.ff"
s.input.ams.MolecularDynamics.CalcPressure = "Yes"
s.input.ams.MolecularDynamics.InitialVelocities.Temperature = 300
s.input.ams.MolecularDynamics.Trajectory.SamplingFreq = 1
s.input.ams.MolecularDynamics.TimeStep = 0.5
s.input.ams.MolecularDynamics.NSteps = 2000
s.runscript.nproc = 1
os.environ["OMP_NUM_THREADS"] = "1"
job = AMSJob(settings=s, molecule=mol, name="md")
job.run()
PLAMS working folder: /path/plams/examples/BasicMDAnalysis/plams_workdir
[03.03|11:50:24] JOB md STARTED
[03.03|11:50:24] JOB md RUNNING
[03.03|11:50:28] JOB md FINISHED
[03.03|11:50:28] JOB md SUCCESSFUL





<scm.plams.interfaces.adfsuite.ams.AMSResults at 0x12f4185b0>
plot_molecule(mol, rotation=("80x,10y,0z"));
../../_images/BasicMDPostanalysis_5_0.png

Or alternatively, load a previously run MD job:

# job = AMSJob.load_external('/path/to/ams.rkf')
results = job.results

Velocity Autocorrelation Function

times, vacf = results.get_velocity_acf(start_fs=0, max_dt_fs=250, normalize=False)
normalized_vacf = vacf / vacf[0]
plt.plot(times, normalized_vacf)
plt.xlabel("Time (fs)")
plt.ylabel("Velocity autocorrelation function")
plt.title("Velocity autocorrelation function")
plt.savefig("plams_vacf.png")
A = np.stack((times, normalized_vacf), axis=1)
np.savetxt("plams_vacf.txt", A, header="Time(fs) VACF")
../../_images/BasicMDPostanalysis_9_0.png

Diffusion Coefficient

t_D, D = results.get_diffusion_coefficient_from_velocity_acf(times, vacf)
plt.plot(t_D, D)
plt.xlabel("Time (fs)")
plt.ylabel("D (m²s⁻¹)")
plt.title("Diffusion coefficient")
plt.savefig("plams_vacf_D.png")
A = np.stack((t_D, D), axis=1)
np.savetxt("plams_vacf_D.txt", A, header="time(fs) D(m^2*s^-1)")
../../_images/BasicMDPostanalysis_11_0.png

Power Spectrum

freq, intensities = results.get_power_spectrum(times, vacf, number_of_points=1000)
plt.plot(freq, intensities)
plt.xlabel("Frequency (cm⁻¹)")
plt.ylabel("Power spectrum (arbitrary units)")
plt.title("Power spectrum")
plt.savefig("plams_power_spectrum.png")
A = np.stack((freq, intensities), axis=1)
np.savetxt("plams_power_spectrum.txt", A, header="Frequency(cm^-1) PowerSpectrum")
../../_images/BasicMDPostanalysis_13_0.png

Green-Kubo Viscosity

t, viscosity = results.get_green_kubo_viscosity(start_fs=0, max_dt_fs=250)  # do not do this for NPT simulations
plt.plot(t, viscosity)
plt.xlabel("Time (fs)")
plt.ylabel("Viscosity (mPa*s)")
plt.title("Viscosity")
plt.savefig("plams_green_kubo_viscosity.png")
A = np.stack((t, viscosity), axis=1)
np.savetxt("plams_green_kubo_viscosity.txt", A, header="Time(fs) Viscosity(mPa*s)")
../../_images/BasicMDPostanalysis_15_0.png

Density Along Axis

z, density = results.get_density_along_axis(axis="z", density_type="mass", bin_width=0.2, atom_indices=None)
plt.plot(z, density)
plt.xlabel("z coordinate (Å)")
plt.ylabel("Density (g/cm³)")
plt.title("Density along z")
plt.savefig("plams_density_along_z.png")
A = np.stack((z, density), axis=1)
np.savetxt("plams_density_along_z.txt", A, header="z(angstrom) density(g/cm^3)")
../../_images/BasicMDPostanalysis_17_0.png

Complete Python code

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

# ## Initial Imports

import os
import numpy as np
import matplotlib.pyplot as plt
from scm.plams import init, packmol, Settings, AMSJob, from_smiles, plot_molecule


# ## Run Simple MD Simulation of Water

# Run a short MD simulation of 16 water molecules in a box.

# this line is not required in AMS2025+
init()

mol = packmol(from_smiles("O"), n_molecules=16, density=1.0)
s = Settings()
s.input.ams.Task = "MolecularDynamics"
s.input.ReaxFF.ForceField = "Water2017.ff"
s.input.ams.MolecularDynamics.CalcPressure = "Yes"
s.input.ams.MolecularDynamics.InitialVelocities.Temperature = 300
s.input.ams.MolecularDynamics.Trajectory.SamplingFreq = 1
s.input.ams.MolecularDynamics.TimeStep = 0.5
s.input.ams.MolecularDynamics.NSteps = 2000
s.runscript.nproc = 1
os.environ["OMP_NUM_THREADS"] = "1"
job = AMSJob(settings=s, molecule=mol, name="md")
job.run()


plot_molecule(mol, rotation=("80x,10y,0z"))


# Or alternatively, load a previously run MD job:

# job = AMSJob.load_external('/path/to/ams.rkf')
results = job.results


# ## Velocity Autocorrelation Function

times, vacf = results.get_velocity_acf(start_fs=0, max_dt_fs=250, normalize=False)
normalized_vacf = vacf / vacf[0]
plt.plot(times, normalized_vacf)
plt.xlabel("Time (fs)")
plt.ylabel("Velocity autocorrelation function")
plt.title("Velocity autocorrelation function")
plt.savefig("plams_vacf.png")
A = np.stack((times, normalized_vacf), axis=1)
np.savetxt("plams_vacf.txt", A, header="Time(fs) VACF")


# ## Diffusion Coefficient

t_D, D = results.get_diffusion_coefficient_from_velocity_acf(times, vacf)
plt.plot(t_D, D)
plt.xlabel("Time (fs)")
plt.ylabel("D (m²s⁻¹)")
plt.title("Diffusion coefficient")
plt.savefig("plams_vacf_D.png")
A = np.stack((t_D, D), axis=1)
np.savetxt("plams_vacf_D.txt", A, header="time(fs) D(m^2*s^-1)")


# ## Power Spectrum

freq, intensities = results.get_power_spectrum(times, vacf, number_of_points=1000)
plt.plot(freq, intensities)
plt.xlabel("Frequency (cm⁻¹)")
plt.ylabel("Power spectrum (arbitrary units)")
plt.title("Power spectrum")
plt.savefig("plams_power_spectrum.png")
A = np.stack((freq, intensities), axis=1)
np.savetxt("plams_power_spectrum.txt", A, header="Frequency(cm^-1) PowerSpectrum")


# ## Green-Kubo Viscosity

t, viscosity = results.get_green_kubo_viscosity(start_fs=0, max_dt_fs=250)  # do not do this for NPT simulations
plt.plot(t, viscosity)
plt.xlabel("Time (fs)")
plt.ylabel("Viscosity (mPa*s)")
plt.title("Viscosity")
plt.savefig("plams_green_kubo_viscosity.png")
A = np.stack((t, viscosity), axis=1)
np.savetxt("plams_green_kubo_viscosity.txt", A, header="Time(fs) Viscosity(mPa*s)")


# ## Density Along Axis

z, density = results.get_density_along_axis(axis="z", density_type="mass", bin_width=0.2, atom_indices=None)
plt.plot(z, density)
plt.xlabel("z coordinate (Å)")
plt.ylabel("Density (g/cm³)")
plt.title("Density along z")
plt.savefig("plams_density_along_z.png")
A = np.stack((z, density), axis=1)
np.savetxt("plams_density_along_z.txt", A, header="z(angstrom) density(g/cm^3)")