I am getting some conflicting results on this topic while searching for literature online, so I am asking this question here with respect to a simplified experimental setup. If we are trying to theoretically model an experimental setup for excitation and relaxation of a molecule in solution due to external radiation, there should be a sink for the energy dissipation, which is observed in the experiment as the heat dissipation from the solution. So in my understanding there should be an external "sink" attached to the MD setup. On the other hand the energy dissipation in the solution then becomes dependent on the coupling strength between the MD system and the thermostat, which can bias the results. So how should I approach this issue?

P.S.: It would be very helpful if there are some relevant literature pointed in the answers.

  • 2
    $\begingroup$ I think the dependence on the heat sink should vanish in the limit that you use an infinite number of solvent molecules. So the NVT ensemble should be usable when you include a sufficient number of solvent molecules. However, I don't know what the required number is, and whether it is beyond what is affordable for a typical user. $\endgroup$
    – wzkchem5
    Mar 13, 2022 at 8:19
  • $\begingroup$ Hi, Thanks for your response. My thinking is that the sufficiently large number of solvent molecules would act as a heat sink compared to a smaller system. In that aspect, for a small system you might need the thermostat to model the dissipation rate correctly, while for a larger system, the available solvent DOF would be available for heat transfer, so we might not need a thermostat (NVE). Does this makes sense? $\endgroup$
    – mykd
    Mar 13, 2022 at 8:40
  • $\begingroup$ Yes, if the system is enormously large, a heat bath is not needed at all. For some large but not excessively large system size, you may have a scenario where a heat bath is necessary, but the coupling constant of the bath with the system is not important as long as it is not unrealistically small and not unrealistically large. The problem is I have no idea of the exact range of the system size for which you encounter the latter scenario. $\endgroup$
    – wzkchem5
    Mar 13, 2022 at 10:40
  • 1
    $\begingroup$ I don't do NEMD, but, this might be useful? or a start? youtube.com/watch?v=DgAFUfspntY&t=802s $\endgroup$
    – B. Kelly
    Mar 15, 2022 at 23:41

1 Answer 1



Thermostatting / barostatting in ensembles such as NVT, NPT, etc, introduces an exchange of energy between the simulated system and a virtual "bath". You should ensure that whatever you are observing in simulation is a property of the system, not of its coupling to the thermostat / barostat. For a first pass, I would recommend equilibrating the system in a relevant ensemble, then switching to NVE for the excitation and subsequent dynamics.

The Physics

It looks like you are interested in exciting some species in a simulation volume, and tracking where that energy subsequently goes. Most likely, (1) you have some species or degrees of freedom to which energy can be transferred, and (2) the excitation is with an electric field / photon. Quantum-mechanically, the absorption of a photon is fast. In a classical potentials-based molecular-dynamics simulation, a classical interaction of charged particles with a field can also be simulated - ignoring polarizability or including it with a shell model.

Most likely, absorption / excitation is the fastest process. The next slower one is the coupling of the initially excited degree of freedom (electronic / phononic / translation / rotation) to others within the system. Most likely (see "nonlinearity" in limitations), the excited degree of freedom will transfer energy to lower-energy ones, which are definitionally lower frequency and work on longer timescales. If excitations are sparse and/or localized, there will ultimately be some heat diffusion happening as well. While all of this is happening, the system has more energy than it should, and is out of equilibrium.

Here is the most important part: in simulation, the thermal coupling to "bath" is usually much faster than in experiment, because simulation is a much smaller system than in reality. Heat transfer from something like an excited molecule to the substrate on which it sits takes nanoseconds in real life, see for example Fig 7 in this paper. But we cannot tolerate a thermostat that takes so long in simulation, which means that thermostatting pulls energy from the simulation more rapidly than is really physical. It does so by coupling to higher-energy degrees of freedom in the system and slowing them down.

Practically, this means that thermostatting will affect more dynamics than is physical in a non-equilibrium simulation with excitations. The easiest solution to this is to equilibrate the system first, then simulate the non-equilibrium dynamics in NVE, and note the temperature rise. The next easiest is to simulate a range of thermostatting parameters, and only trust the results in the range where changing the thermostat does not change the observed dynamics.


There is a rigorous body of work simulating the excitations of molecular motions and vibrations in liquids. Notably, they have used both classical MD and ab initio MD. The ab initio lets them study polarizability and second-order (induced-dipole) effects. In approximate chronological order, there are a few pure simulation papers, one and two and three, and a few experiment + simulation papers: one and two and three. I am not sure what to say or highlight about these works other than they have looked at a variety of aspects of this type of simulation, and benchmarked them with experiments. A keen reader could pay attention to the analyses used to track the partitioning of energy between various degrees of freedom, and how the simulated excitation strengths correspond to the experimental ones.

Absolutely shameless plug of my own work: I've used classical MD to look at the excitation of ionic translations with light, pre-print here. I took five steps to ensure I could trust the results: (1) used a system with the fewest possible confounding degrees of freedom (2) used NVE and noted the temperature rise (3) made a very large simulation to get enough rare-event statistics (4) used an excitation field that is equal to the experimental one to minimize confounding contributions from nonlinearity (5) interpreted only the most salient, qualitative trends rather than the quantitative rates. The simulation is probably quantitatively wrong, but the trends are probably qualitatively right.


First, this answer deals mostly with classical, nuclear-dynamics stuff. It may be possible to just put a field into an ab initio simulation and simulate electronic excitations. But I don't know for sure. If you want to do a proper electronic excitation, someone else should answer with how best to do it in ab initio, TD-DFT, etc.

Second, the excitation and dynamics of interest may be rare; beware averages. This happens most at higher frequencies / energies as quantum effect become prominent. See for example this paper about water (disclosure: I worked with some of these people), where to match the simulated effect to experiment the excitation fraction had to be taken into account. Experimental excitation strength / density may be much smaller than what is practical to simulate while still getting a simulated signal - but the dynamics at simulate-able strength may not be the same as at an experimental one. Separately, the dynamics of interest could also be a small shoulder superimposed on some larger effect you may not care about, as in my work linked above.

Third, nonlinearity happens. Nonlinearity and anharmonicity of vibrations are integral to the transfer of energy between degrees of freedom, and to heat transport in the first place. Polarizability and induced-dipole effects may happen. If electrons are involved, weird stuff like impact ionization may happen. Nonlinear processes such as phonon upconversion turn on with excitation strength. In general, I recommend being very careful about the excitation strength used in simulation as that will affect subsequent dynamics: try a few, or go as low as possible, or match an experiment.

  • $\begingroup$ Really great answer! Specifically "in simulation, the thermal coupling to "bath" is usually much faster than in experiment, because simulation is a much smaller system than in reality." I guess I completely forgot about the heat transfer rate and how the incident field intensity might affect it. $\endgroup$
    – mykd
    Jul 25, 2022 at 10:45

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .