Gasphase IR spectrum from Molecular Dynamics
********************************************

GFNFF Spectroscopy IR spectrum Molecular Dynamics

**Requires:** AMS2026 or later

-[ Related examples ]-

* Molecular Dynamics with Python

* Conformer Generation, Rescoring, and Filtering

* Conformer Generation for Multiple Molecules

* PLUMED Biasing in AMS Molecular Dynamics

* Automated Transition-State Workflow for a Diels-Alder Addition

-[ Related tutorials ]-

* IR spectrum of ethane tutorial

-[ Related documentation ]-

* Molecular Dynamics task

* GFNFF engine

[image]


Initial imports
===============

   import scm.plams as plams
   import matplotlib.pyplot as plt
   import numpy as np

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

   PLAMS working folder: /path/to/plams_workdir.002


Molecule
========

   from scm.base import ChemicalSystem
   from scm.plams import view

   mol = ChemicalSystem.from_smiles("NC(CO)OCC=O")
   view(mol, width=200, height=200)

[image]


Engine settings
===============

   s = plams.Settings()
   s.input.GFNFF  # GFN-FF force field
   # s.input.MLPotential.Model = "AIMNet2-wB97MD3"   # new ml model in AMS2024

   print(plams.AMSJob(settings=s).get_input())

   # run in serial
   s.runscript.nproc = 1

   Engine GFNFF
   EndEngine


Some general parameters
=======================

   T = 298  # temperature in K
   max_freq = 4000  # plot spectrum up to 4000 cm^-1
   max_dt_fs = 2000  # maximum correlation in fs for dipole derivative acf


Equilibration
=============

"temperature=(500, T, T)" means that in the first half the simulation
the system is cooled from 500 K to the given temperature, and then
kept constant at that temperature.

The initial temperature of 500 K does some preliminary conformer
search.

   eq_job = plams.AMSNVTJob(
       settings=s,
       name="nvt_eq",
       molecule=mol,
       timestep=0.5,
       nsteps=10000,
       temperature=(500, T, T),
   )
   eq_job.run();

   [11.02|12:04:50] JOB nvt_eq STARTED
   [11.02|12:04:50] JOB nvt_eq RUNNING
   [11.02|12:04:52] JOB nvt_eq FINISHED
   [11.02|12:04:52] JOB nvt_eq SUCCESSFUL


NVE production simulation
=========================

The "binlog_dipolemoment" option stores the dipole moment at every
time step.

   job = plams.AMSNVEJob.restart_from(
       eq_job,
       name="nve_single_prod",
       nsteps=50000,
       binlog_dipolemoment=True,
       binlog_time=True,
       samplingfreq=100,
       timestep=0.5,
   )
   job.run();

   [11.02|12:04:52] JOB nve_single_prod STARTED
   [11.02|12:04:52] JOB nve_single_prod RUNNING
   [11.02|12:05:01] JOB nve_single_prod FINISHED
   [11.02|12:05:01] JOB nve_single_prod SUCCESSFUL


Dipole derivative autocorrelation function
==========================================

   times, dipole_deriv_acf = job.results.get_dipole_derivatives_acf(
       start_fs=0, max_dt_fs=max_dt_fs
   )
   plt.plot(times, dipole_deriv_acf)
   plt.xlabel("Time (fs)")
   plt.ylabel("Dipole deriv. autocorrelation (e bohr / fs)^2")
   plt.title("Raw autocorrelation function");

[image]

Ideally, you should set "max_dt_fs" above to a large enough number so
that the autocorrelation function decreases to a constant value of 0
(**and** have a long enough MD simulation to get enough statistics!)


IR spectrum
===========

The IR spectrum is the Fourier transform of the above autocorrelation
function:

   x_freq, y_intens_raw = job.results.get_ir_spectrum_md(
       times, dipole_deriv_acf, max_freq=max_freq
   )
   plt.plot(x_freq, y_intens_raw)
   plt.xlabel("Frequency (cm^-1)")
   plt.title("IR spectrum (from raw autocorrelation function)")
   plt.xlim(500, max_freq);

[image]

There seems to be quite some “noise” in the IR spectrum. One reason
for this is that there is still some signal (or noise?) in the
autocorrelation function at dt = 2000 fs.

However, it’s also possible to use a tapering (window) function to
make the autocorrelation function smoothly decrease to 0. This will
make the resulting IR spectrum look a bit more tidy. See the next
section.


Tapering function for autocorrelation function
==============================================

   def tapered_cosine(x):
       return 0.5 * (np.cos(np.pi * x / np.max(x)) + 1)


   plt.plot(times, tapered_cosine(times))
   plt.title("Tapering / cutoff / window function");

[image]

Now apply this function to the autocorrelation function:

   dipole_deriv_acf_tapered_cosine = dipole_deriv_acf * tapered_cosine(times)
   plt.plot(times, dipole_deriv_acf_tapered_cosine)
   plt.xlabel("Time (fs)")
   plt.title("Autocorrelation cosine tapering");

[image]

And calculate the IR spectrum:

   x_freq, y_intens_cosine = job.results.get_ir_spectrum_md(
       times, dipole_deriv_acf_tapered_cosine, max_freq=max_freq
   )
   plt.plot(x_freq, y_intens_cosine)
   plt.xlabel("Frequency (cm^-1)")
   plt.title("IR spectrum (from cosine-tapered autocorrelation function)")
   plt.xlim(500, max_freq);

[image]

Above we see that using the cosine-tapered autocorrelation function
gives a smoother IR spectrum without affecting the intensities too
much.


Compare to IR spectrum calculated from harmonic approximation
=============================================================

Let’s compare to an IR spectrum calculated with a geometry
optimization + frequencies job, starting from the final frame of the
MD simulation.

   ams_s = plams.Settings()
   ams_s.input.ams.Task = "GeometryOptimization"
   ams_s.input.ams.Properties.NormalModes = "Yes"
   harmonic_job = plams.AMSJob(settings=ams_s + s, name="harmonic", molecule=mol)
   harmonic_job.run();

   [11.02|12:05:02] JOB harmonic STARTED
   [11.02|12:05:02] JOB harmonic RUNNING
   [11.02|12:05:02] JOB harmonic FINISHED
   [11.02|12:05:02] JOB harmonic SUCCESSFUL

   harmonic_freq, harmonic_intens = harmonic_job.results.get_ir_spectrum(
       broadening_type="lorentzian", broadening_width=20
   )
   plt.plot(harmonic_freq, harmonic_intens)
   rescale_factor = np.sum(harmonic_intens) / np.sum(y_intens_cosine)
   plt.plot(x_freq, y_intens_cosine * rescale_factor)  # rescale
   plt.legend(["Harmonic", "MD NVE"])
   plt.title("IR spectrum")
   plt.xlabel("Frequency (cm^-1)")
   plt.xlim(500, max_freq);

[image]

In this case, the MD simulation samples multiple conformers so there
are more peaks than for the harmonic calculation.


View the trajectory in AMSmovie
===============================

   !amsmovie "{job.results.rkfpath()}"


View the normal modes in AMSspectra
===================================

   !amsspectra "{harmonic_job.results.rkfpath(file='engine')}"


Appendix: Average over multiple short NVE simulations
=====================================================

Best practice is to run multiple NVE simulations starting from
different points of the NVT simulation, assuming that the NVT
simulation samples the correct equilibrium distribution of
structures/conformers.

Let’s make this more explicit with another NVT simulation, followed by
multiple NVE simulations from different points of the NVT simulation.
See also the “Molecular Dynamics with Python” tutorial.

   nvt_prod_job = plams.AMSNVTJob.restart_from(
       eq_job,
       name="nvt_prod",
       nsteps=50000,
       samplingfreq=100,
       timestep=0.5,
       thermostat="NHC",
       tau=100,
       temperature=T,
   )
   nvt_prod_job.run(watch=True);

   [11.02|12:05:10] JOB nvt_prod STARTED
   [11.02|12:05:10] JOB nvt_prod RUNNING
   [11.02|12:05:10] nvt_prod: AMS 2025.206  RunTime: Feb11-2026 12:05:10  ShM Nodes: 1  Procs: 1
   [11.02|12:05:10] nvt_prod: Starting MD calculation:
   [11.02|12:05:10] nvt_prod: --------------------
   ... output trimmed ....
   [11.02|12:05:19] nvt_prod:         50000   25000.00     252.      -2.55492            n/a            0.1
   [11.02|12:05:19] nvt_prod: MD calculation finished.
   [11.02|12:05:19] nvt_prod: NORMAL TERMINATION
   [11.02|12:05:19] JOB nvt_prod FINISHED
   [11.02|12:05:19] JOB nvt_prod SUCCESSFUL

   nvespawner_job = plams.AMSNVESpawnerJob(
       nvt_prod_job,
       name="nvespawner-" + nvt_prod_job.name,
       n_nve=10,  # the number of NVE simulations to run
       timestep=0.5,
       binlog_time=True,
       binlog_dipolemoment=True,
       nsteps=20000,
   )
   nvespawner_job.run();

   [11.02|12:05:19] JOB nvespawner-nvt_prod STARTED
   [11.02|12:05:19] JOB nvespawner-nvt_prod RUNNING
   [11.02|12:05:19] JOB nvespawner-nvt_prod/nve1 STARTED
   [11.02|12:05:19] JOB nvespawner-nvt_prod/nve1 RUNNING
   [11.02|12:05:23] JOB nvespawner-nvt_prod/nve1 FINISHED
   ... output trimmed ....
   [11.02|12:05:55] JOB nvespawner-nvt_prod/nve10 RUNNING
   [11.02|12:05:58] JOB nvespawner-nvt_prod/nve10 FINISHED
   [11.02|12:05:59] JOB nvespawner-nvt_prod/nve10 SUCCESSFUL
   [11.02|12:05:59] JOB nvespawner-nvt_prod FINISHED
   [11.02|12:05:59] JOB nvespawner-nvt_prod SUCCESSFUL

Let’s check that the temperature during the NVE is not too far from
the requested temperature.

   avg_T = nvespawner_job.results.get_mean_temperature()

   print(f"Set temperature during NVT: {T:.1f} K")
   print(f"Mean temperature during NVE: {avg_T:.1f}")

   Set temperature during NVT: 298.0 K
   Mean temperature during NVE: 289.4

Calculate the average dipole derivative autocorrelation function.

To calculate the IR spectrum from our custom set of averaged data, we
directly call the "power_spectrum" function from PLAMS:

   avg_x, avg_y = nvespawner_job.results.get_dipole_derivatives_acf(
       start_fs=0, max_dt_fs=max_dt_fs
   )
   avg_y *= tapered_cosine(avg_x)

   x_freq_multiple, y_intens_cosine_multiple = plams.trajectories.analysis.power_spectrum(
       times, avg_y, max_freq=max_freq
   )

   plt.plot(x_freq, y_intens_cosine)
   plt.plot(x_freq_multiple, y_intens_cosine_multiple)
   plt.xlabel("Frequency (cm^-1)")
   plt.title("IR spectrum from multiple NVE simulations")
   plt.legend(["single", "multiple"])
   plt.xlim(500, 4000);

[image]


See also
========


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

* Molecular Dynamics with Python

* Conformer Generation, Rescoring, and Filtering

* Conformer Generation for Multiple Molecules

* PLUMED Biasing in AMS Molecular Dynamics

* Automated Transition-State Workflow for a Diels-Alder Addition


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

* IR spectrum of ethane tutorial


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

* Molecular Dynamics task

* GFNFF engine


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

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

   # ## Initial imports

   import scm.plams as plams
   import matplotlib.pyplot as plt
   import numpy as np

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


   # ## Molecule

   from scm.base import ChemicalSystem
   from scm.plams import view

   mol = ChemicalSystem.from_smiles("NC(CO)OCC=O")
   view(mol, width=200, height=200, picture_path="picture1.png")


   # ## Engine settings

   s = plams.Settings()
   s.input.GFNFF  # GFN-FF force field
   # s.input.MLPotential.Model = "AIMNet2-wB97MD3"   # new ml model in AMS2024

   print(plams.AMSJob(settings=s).get_input())

   # run in serial
   s.runscript.nproc = 1


   # ## Some general parameters

   T = 298  # temperature in K
   max_freq = 4000  # plot spectrum up to 4000 cm^-1
   max_dt_fs = 2000  # maximum correlation in fs for dipole derivative acf


   # ## Equilibration
   #
   # ``temperature=(500, T, T)`` means that in the first half the simulation the system is cooled from 500 K to the given temperature, and then kept constant at that temperature.
   #
   # The initial temperature of 500 K does some preliminary conformer search.

   eq_job = plams.AMSNVTJob(
       settings=s,
       name="nvt_eq",
       molecule=mol,
       timestep=0.5,
       nsteps=10000,
       temperature=(500, T, T),
   )
   eq_job.run()
   # ## NVE production simulation
   #
   # The ``binlog_dipolemoment`` option stores the dipole moment at every time step.

   job = plams.AMSNVEJob.restart_from(
       eq_job,
       name="nve_single_prod",
       nsteps=50000,
       binlog_dipolemoment=True,
       binlog_time=True,
       samplingfreq=100,
       timestep=0.5,
   )
   job.run()
   # ## Dipole derivative autocorrelation function

   times, dipole_deriv_acf = job.results.get_dipole_derivatives_acf(start_fs=0, max_dt_fs=max_dt_fs)
   plt.plot(times, dipole_deriv_acf)
   plt.xlabel("Time (fs)")
   plt.ylabel("Dipole deriv. autocorrelation (e bohr / fs)^2")
   plt.title("Raw autocorrelation function")
   plt.gcf().savefig("picture2.png")
   plt.clf()


   # Ideally, you should set ``max_dt_fs`` above to a large enough number so that the autocorrelation function decreases to a constant value of 0 (**and** have a long enough MD simulation to get enough statistics!)

   # ## IR spectrum
   #
   # The IR spectrum is the Fourier transform of the above autocorrelation function:

   x_freq, y_intens_raw = job.results.get_ir_spectrum_md(times, dipole_deriv_acf, max_freq=max_freq)
   plt.plot(x_freq, y_intens_raw)
   plt.xlabel("Frequency (cm^-1)")
   plt.title("IR spectrum (from raw autocorrelation function)")
   plt.xlim(500, max_freq)
   plt.gcf().savefig("picture3.png")
   plt.clf()


   # There seems to be quite some "noise" in the IR spectrum. One reason for this is that there is still some signal (or noise?) in the autocorrelation function at dt = 2000 fs.
   #
   # However, it's also possible to use a tapering (window) function to make the autocorrelation function smoothly decrease to 0. This will make the resulting IR spectrum look a bit more tidy. See the next section.

   # ## Tapering function for autocorrelation function


   def tapered_cosine(x):
       return 0.5 * (np.cos(np.pi * x / np.max(x)) + 1)


   plt.plot(times, tapered_cosine(times))
   plt.title("Tapering / cutoff / window function")
   plt.gcf().savefig("picture4.png")
   plt.clf()


   # Now apply this function to the autocorrelation function:

   dipole_deriv_acf_tapered_cosine = dipole_deriv_acf * tapered_cosine(times)
   plt.plot(times, dipole_deriv_acf_tapered_cosine)
   plt.xlabel("Time (fs)")
   plt.title("Autocorrelation cosine tapering")
   plt.gcf().savefig("picture5.png")
   plt.clf()


   # And calculate the IR spectrum:

   x_freq, y_intens_cosine = job.results.get_ir_spectrum_md(times, dipole_deriv_acf_tapered_cosine, max_freq=max_freq)
   plt.plot(x_freq, y_intens_cosine)
   plt.xlabel("Frequency (cm^-1)")
   plt.title("IR spectrum (from cosine-tapered autocorrelation function)")
   plt.xlim(500, max_freq)
   plt.gcf().savefig("picture6.png")
   plt.clf()


   # Above we see that using the cosine-tapered autocorrelation function gives a smoother IR spectrum without affecting the intensities too much.

   # ## Compare to IR spectrum calculated from harmonic approximation

   # Let's compare to an IR spectrum calculated with a geometry optimization + frequencies job, starting from the final frame of the MD simulation.

   ams_s = plams.Settings()
   ams_s.input.ams.Task = "GeometryOptimization"
   ams_s.input.ams.Properties.NormalModes = "Yes"
   harmonic_job = plams.AMSJob(settings=ams_s + s, name="harmonic", molecule=mol)
   harmonic_job.run()
   harmonic_freq, harmonic_intens = harmonic_job.results.get_ir_spectrum(broadening_type="lorentzian", broadening_width=20)
   plt.plot(harmonic_freq, harmonic_intens)
   rescale_factor = np.sum(harmonic_intens) / np.sum(y_intens_cosine)
   plt.plot(x_freq, y_intens_cosine * rescale_factor)  # rescale
   plt.legend(["Harmonic", "MD NVE"])
   plt.title("IR spectrum")
   plt.xlabel("Frequency (cm^-1)")
   plt.xlim(500, max_freq)
   plt.gcf().savefig("picture7.png")
   plt.clf()


   # In this case, the MD simulation samples multiple conformers so there are more peaks than for the harmonic calculation.

   # ## View the trajectory in AMSmovie

   # get_ipython().system('amsmovie "{job.results.rkfpath()}"')


   # ## View the normal modes in AMSspectra

   # get_ipython().system('amsspectra "{harmonic_job.results.rkfpath(file=\'engine\')}"')


   # ## Appendix: Average over multiple short NVE simulations
   #
   # Best practice is to run multiple NVE simulations starting from different points of the NVT simulation, assuming that the NVT simulation samples the correct equilibrium distribution of structures/conformers.
   #
   # Let's make this more explicit with another NVT simulation, followed by multiple NVE simulations from different points of the NVT simulation. See also the "Molecular Dynamics with Python" tutorial.

   nvt_prod_job = plams.AMSNVTJob.restart_from(
       eq_job,
       name="nvt_prod",
       nsteps=50000,
       samplingfreq=100,
       timestep=0.5,
       thermostat="NHC",
       tau=100,
       temperature=T,
   )
   nvt_prod_job.run(watch=True)
   nvespawner_job = plams.AMSNVESpawnerJob(
       nvt_prod_job,
       name="nvespawner-" + nvt_prod_job.name,
       n_nve=10,  # the number of NVE simulations to run
       timestep=0.5,
       binlog_time=True,
       binlog_dipolemoment=True,
       nsteps=20000,
   )
   nvespawner_job.run()
   # Let's check that the temperature during the NVE is not too far from the requested temperature.

   avg_T = nvespawner_job.results.get_mean_temperature()

   print(f"Set temperature during NVT: {T:.1f} K")
   print(f"Mean temperature during NVE: {avg_T:.1f}")


   # Calculate the average dipole derivative autocorrelation function.
   #
   # To calculate the IR spectrum from our custom set of averaged data, we directly call the ``power_spectrum`` function from PLAMS:

   avg_x, avg_y = nvespawner_job.results.get_dipole_derivatives_acf(start_fs=0, max_dt_fs=max_dt_fs)
   avg_y *= tapered_cosine(avg_x)

   x_freq_multiple, y_intens_cosine_multiple = plams.trajectories.analysis.power_spectrum(times, avg_y, max_freq=max_freq)


   plt.plot(x_freq, y_intens_cosine)
   plt.plot(x_freq_multiple, y_intens_cosine_multiple)
   plt.xlabel("Frequency (cm^-1)")
   plt.title("IR spectrum from multiple NVE simulations")
   plt.legend(["single", "multiple"])
   plt.xlim(500, 4000)
   plt.gcf().savefig("picture8.png")
   plt.clf()
