3
$\begingroup$

How can we check whether or not the power spectra (vibrational density of states) calculated by taking the Fourier transform of the velocity auto-correlation function (VACF) from a molecular dynamics trajectory is converged? How to estimate the statistical error in it?

What is the meaning of the word "convergence" in this context?

$\endgroup$
2
$\begingroup$

There are multiple ways in which the VACF, along with other properties derived from molecular dynamics trajectories like the dipole ACF to give IR spectra, can be considered "converged".

  1. The molecular dynamics trajectory itself must be converged by some time-dependent statistical metric, such as measuring fluctuations in the energy or temperature over a number of steps and ensuring the average is below some threshold. This will be related to the ensemble chosen, the thermostat and barostat and associated parameters (if relevant), the integrator used for time evolution, and the size of the simulation, among other things.

  2. Molecular simulations are layers of approximations piled on top of each other. In the case of purely classical dynamics, this usually refers to the force field parameters for the solute and solvent. Force fields are designed for different systems and may or may not be transferable, but the goal is usually to accurately predict dynamics. In the case of the VACF, if the dynamics of a system are described correctly, then the resulting power spectrum should be as well. However, in the case of the dipole ACF, the force field parameters enter directly as the nuclear charges (rather than just indirectly through nuclear velocities), which places stricter requirements on the quality of the force field. One hopes that newer force fields are better, but improvement over time isn't necessarily systematic and there is no quality hierarchy.

  3. As an extension of 2, if the simulation is quantum mechanical in any way (such as in QM/MM, Born-Oppenheimer MD, or Car–Parrinello MD), then there are approximations in how the electronic structure is treated. For spectroscopic properties, it is more clear how to systematically improve results with regard to electronic structure than force fields due to well-established quantum chemical method hierarchies. Benchmarks of this hierarchy have been performed on many properties, including geometric parameters, dipole moments, and vibrational frequencies. In theory, the error in both dynamics and spectra should decrease up going from Hartree-Fock to MP2 to CCSD to CCSD(T), or by using more advanced density functionals (say, BP86 to PBE0 to M06-2X to B2PLYP). In practice, this would require benchmarking to draw any conclusions.

  4. As an extension of 3, there are approximations in QM/MM of how the classically-treated and QM-treated regions interact with each other. Rather than go into more detail, here is a review of different QM/MM methods covering different definitions of the boundary and how one region may interact with the other.

Regarding estimation of the statistical error in resulting spectra, I am only aware of direct comparison (absolute error) between experimental and computed spectra, because it is assumed that the MD trajectory is converged (point 1). Performing multiple independent MD simulations and calculating the VACF from each should give an error measure closer to true statistical error, though I am unsure of how this would be computed.

Here is good almost-review (with good examples) of using Born-Oppenheimer molecular dynamics to calculate power spectra, IR spectra, and Raman spectra. I am not aware of any comparisons between classical MD and BOMD for computing spectra.

On the surface, it may appear as though the introduction of quantum mechanics makes things more difficult; certainly things become more time-consuming, however it is much clearer how results will change by treating electronic structure explicitly and with less severe approximations, rather than with a force field.

$\endgroup$
  • $\begingroup$ whatever you are suggesting is related to the quality of MD trajectory, the way to check convergence is still not clear to me! As far as my understanding, we can check whether VACF is converging to zero or not within the simulated time, since if the simulation time is small, it may not get converged! And, if we have long enough trajectory, we might consider spectra calculated from different portion of it, ie taking different initial starting point. $\endgroup$ – Sangkha Borah Aug 19 '16 at 3:56
  • $\begingroup$ To your first comment, "whatever you are suggesting is related to the quality of MD trajectory, the way to check convergence is still not clear to me", you should read the peer-reviewed literature for systems related to yours in order to understand this. This is both fundamental and critical to running any kind of MD simulation. $\endgroup$ – pentavalentcarbon Aug 19 '16 at 4:53
  • $\begingroup$ To your second point, "we can check whether VACF is converging to zero or not within the simulated time", that might be valid, but what does the literature say? It sounds like you have some knowledge of this already. $\endgroup$ – pentavalentcarbon Aug 19 '16 at 4:54
  • $\begingroup$ To your third point, "if we have long enough trajectory, we might consider spectra calculated from different portion of it, ie taking different initial starting point", the simplest case would be to take the entirety of your production run and none of your equilibration run, assuming that your system is properly equilibrated. This would ensure your VACF captures all the points (hopefully unbiased) that your trajectory samples. $\endgroup$ – pentavalentcarbon Aug 19 '16 at 4:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.