.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples_generated/md_acetonit/acn.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_generated_md_acetonit_acn.py: .. _acn_md_tutorial: Equilibrating an MD box of acetonitrile ======================================= Goals ===== In this tutorial we learn how to perform a thermal equilibration of a box of acetonitrile molecules using ASE. We will: * build a coarse-grained CH₃–C≡N triatomic model. * set up a periodic box at the experimental density (298 K). * apply rigid-molecule constraints. * use the ACN force field. * equilibrate with Langevin dynamics. * scale small box of 27 molecules to 216 molecules. ACN model ========= The acetonitrile force field implemented in ASE (:mod:`ase.calculators.acn`) is an interaction potential between three-site linear molecules. The atoms of the methyl group are treated as a single site centered on the methyl carbon (hydrogens are not explicit). Therefore: * assign the **methyl mass** to the outer carbon (``m_me``), * use the atomic sequence **Me–C–N** repeated for all molecules, * keep molecules **rigid** during MD with :class:`FixLinearTriatomic`. .. GENERATED FROM PYTHON SOURCE LINES 31-43 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np import ase.units as units from ase import Atoms from ase.calculators.acn import ACN, m_me, r_cn, r_mec from ase.constraints import FixLinearTriatomic from ase.io import Trajectory from ase.md import Langevin from ase.visualize.plot import plot_atoms .. GENERATED FROM PYTHON SOURCE LINES 44-49 Step 1: molecule ---------------- Build one CH3–C≡N as a linear triatomic “C–C–N”. The first carbon is the methyl site; we assign it the CH3 mass. Rotate slightly to avoid perfect alignment with the cell axes. .. GENERATED FROM PYTHON SOURCE LINES 49-74 .. code-block:: Python pos = [[0, 0, -r_mec], [0, 0, 0], [0, 0, r_cn]] atoms = Atoms('CCN', positions=pos) atoms.rotate(30, 'x') masses = atoms.get_masses() masses[0] = m_me atoms.set_masses(masses) mol = atoms.copy() mol.set_pbc(False) mol.set_cell([20, 20, 20]) mol.center() fig, ax = plt.subplots(figsize=(4, 4)) plot_atoms( mol, ax=ax, rotation='30x,30y,0z', show_unit_cell=0, radii=0.75, ) ax.set_axis_off() plt.tight_layout() plt.show() .. image-sg:: /examples_generated/md_acetonit/images/sphx_glr_acn_001.png :alt: acn :srcset: /examples_generated/md_acetonit/images/sphx_glr_acn_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 75-79 Step 2: Set up small box of 27-molecules ---------------------------------------- Match density 0.776 g/cm^3 at 298 K. Compute cubic box length L from mass and density. Build 3×3×3 supercell and enable PBC. .. GENERATED FROM PYTHON SOURCE LINES 79-100 .. code-block:: Python density = 0.776 / 1e24 # g / Å^3 L = ((masses.sum() / units.mol) / density) ** (1 / 3.0) atoms.set_cell((L, L, L)) atoms.center() atoms = atoms.repeat((3, 3, 3)) atoms.set_pbc(True) box27 = atoms.copy() fig, ax = plt.subplots(figsize=(5, 5)) plot_atoms( box27, ax=ax, rotation='35x,35y,0z', show_unit_cell=2, radii=0.75, ) ax.set_axis_off() plt.tight_layout() plt.show() .. image-sg:: /examples_generated/md_acetonit/images/sphx_glr_acn_002.png :alt: acn :srcset: /examples_generated/md_acetonit/images/sphx_glr_acn_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 101-104 Step 3: Set constraints ----------------------- Keep each “C–C–N” rigid during MD using FixLinearTriatomic. .. GENERATED FROM PYTHON SOURCE LINES 104-108 .. code-block:: Python nm = 27 triples = [(3 * i, 3 * i + 1, 3 * i + 2) for i in range(nm)] atoms.constraints = FixLinearTriatomic(triples=triples) .. GENERATED FROM PYTHON SOURCE LINES 109-113 Step 4: MD run for 27-molecules system -------------------------------------- Assign ACN with cutoff = half the smallest box edge. Langevin MD at 300 K, 1 fs timestep. Save a frame every step. .. GENERATED FROM PYTHON SOURCE LINES 113-127 .. code-block:: Python atoms.calc = ACN(rc=np.min(np.diag(atoms.cell)) / 2) tag = 'acn_27mol_300K' md = Langevin( atoms, 1 * units.fs, temperature_K=300, friction=0.01, logfile=tag + '.log', ) traj = Trajectory(tag + '.traj', 'w', atoms) md.attach(traj.write, interval=10) md.run(5000) # 5 ps @ 1 fs .. rst-class:: sphx-glr-script-out .. code-block:: none True .. GENERATED FROM PYTHON SOURCE LINES 128-132 Step 5: scale system size to 216 molecules ------------------------------------------ Repeat 2×2×2 to reach 216 molecules. Reapply constraints and update the ACN cutoff for the new cell. .. GENERATED FROM PYTHON SOURCE LINES 132-141 .. code-block:: Python atoms.set_constraint() atoms = atoms.repeat((2, 2, 2)) nm = 216 triples = [(3 * i, 3 * i + 1, 3 * i + 2) for i in range(nm)] atoms.constraints = FixLinearTriatomic(triples=triples) atoms.calc = ACN(rc=np.min(np.diag(atoms.cell)) / 2) .. GENERATED FROM PYTHON SOURCE LINES 142-144 Step 6: MD run for 216-molecules system --------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 144-177 .. code-block:: Python tag = 'acn_216mol_300K' md = Langevin( atoms, 2 * units.fs, temperature_K=300, friction=0.01, logfile=tag + '.log', ) traj = Trajectory(tag + '.traj', 'w', atoms) md.attach(traj.write, interval=10) times_ps, epots, ekins, etots, temps = [], [], [], [], [] sample_interval = 10 # sample every 10 MD steps for lighter plots def sample(): # Time in ps (same as MDLogger: dyn.get_time() / (1000 * units.fs)) t_ps = md.get_time() / (1000.0 * units.fs) ep = atoms.get_potential_energy() # eV total ek = atoms.get_kinetic_energy() # eV total T = atoms.get_temperature() # K times_ps.append(t_ps) epots.append(ep) ekins.append(ek) etots.append(ep + ek) temps.append(T) # initial sample at t=0 sample() md.attach(sample, interval=sample_interval) md.run(1000) # 6 ps @ 2 fs .. rst-class:: sphx-glr-script-out .. code-block:: none True .. GENERATED FROM PYTHON SOURCE LINES 178-181 Plot Instantaneous temperature vs time. Does the system equilibrated well? What is the average temperature? Should we run longer simulations? .. GENERATED FROM PYTHON SOURCE LINES 181-190 .. code-block:: Python fig, ax = plt.subplots(figsize=(6, 4)) ax.plot(times_ps, temps, label='T (K)') ax.set_xlabel('Time (ps)') ax.set_ylabel('Temperature (K)') ax.legend(loc='best') ax.grid(True, linewidth=0.5, alpha=0.5) plt.tight_layout() plt.show() .. image-sg:: /examples_generated/md_acetonit/images/sphx_glr_acn_003.png :alt: acn :srcset: /examples_generated/md_acetonit/images/sphx_glr_acn_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 191-201 Next steps ---------- * View trajectories:: ase gui acn_27mol_300K.traj ase gui acn_216mol_300K.traj * Plot other thermodynamic quantities (p.e., k.e., and total energy). .. rst-class:: sphx-glr-timing **Total running time of the script:** (4 minutes 7.694 seconds) .. _sphx_glr_download_examples_generated_md_acetonit_acn.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: acn.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: acn.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: acn.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_