Struggle with NX-Hybrid-Surface Hopping

[Anna U]

Hello everyone,

I am currently working with NewtonX to study the Surface Hopping of a fluorophore in solution. I run some Trajectories, but I am not really happy with the result, as the energy fluctuation is  huge (40 to 150 eV). I attached the energy-time-plots ($NX/plot) as a reference. 

I am using the following parameters to run the calculation (I am happy to share more detailed data if needed):

- QM-part: Turbomole TD-DFT  with def2-TZVPD and wb97x

- MM-part: mm3.prm

- initial conditions: random velocities with 1 eV to start with (the hybrid initial conditions setup is not working, as it gets stuck while calling nxinp)

Besides that, I was working along the NX-Tutorial

My Questions are the following:

Where does this huge energy fluctuation come from?

Is it the random velocities which are not good enough?

Is the plotted potential energy the sum of the QM and the MM part and if so, is it possible to plot just the QM-Energy?

Any help is highly appreciated, as I am the only one in my research group working with NX.

Thank you very much!

Anna

1)  4-aminophthalimide in a DMSO-solvent box

2) 4-aminophthalimide in a DMF-solvent box

[Michal Kochman]

Hello, Anna,

When you start the surface hopping simulation, the energy of the occupied state immediately starts to change very rapidly, but the total energy of the system is conserved. Do I understand correctly?

The problem might have to do with the initial conditions (atomic positions and velocities) for the surface hopping simulation. For example, maybe the system was in a far-from-equilibrium configuration at t = 0. It then relaxed towards thermal equilibrium. How did you generate the initial conditions? The usual procedure is to carry out a thermostatted molecular dynamics simulation in which the geometry of the solute is frozen. The solvent is then in equilibrium with the solute.

On an unrelated note, the def2-TZVPD basis set contains diffuse basis functions. This may be problematic in QM/MM calculations if the MM point charges are included in the QM Hamiltonian. Certainly for test simulations a smaller basis set will be sufficient.

Best, Michał

[Anna U]

Hello Michal,

thanks for your fast reply!

Yes you understood that correct. Also  I am unsure about the quite big energy fluctuations during the Trajectory. I guess this is due to the summation of the energies of the MM part. Is this correct?

I equilibrated the starting geometry with Tinker, running a NVT and  a NPT calculation. The whole system was used for it (the solute was not frozen), so I will try to equilibrate with a frozen solute. The initial conditions for NX were created with random velocities at 1 eV for the whole system. I think this is not ideal as well. But I cannot get the hybrid_mixed_initcond.pl to work. It gets stuck while calling nxinp.

I chose the diffuse function, because they are describing the molecule best, but I did not know it could interfere with the point charges. They are included in my calculation. Therefore, I will give it a try with a smaller basis.

Best regards

Anna

[Michal Kochman]

>> I guess this is due to the summation of the energies of the MM part. Is this correct?

In truth, I don't know what's happening. My thinking is as follows: the total energy is conserved, so from the technical side, Newton-X is probably running correctly. Meaning, it's not a bug in the code. But someone who is more knowledgeable about this particular program should weigh in on this.

On the other hand, the state energies look a bit strange, and that may mean a problem with the initial conditions. Maybe they were set up incorrectly, or in a way that is inconsistent with the surface hopping simulation. Generating initial conditions for a surface hopping simulation with QM/MM solvent is actually a complex procedure. (See below.) Hence, my guess is that the sudden jump in the state energies is due to the initial conditions being unrealistic. For example, the solvent may be too hot, or too cold.

When running a surface hopping simulation with QM/MM solvent, you should be seeing something like this:

Here, the time is in fs, and the total energy is in hartree. The state energies are fluctuating somewhat, but there's no net increase or decrease. This is a simulation for fluorazene, a molecule with a long excited-state lifetime, so the ground state (S0) is separated in energy from the excited states (S1-S4) throughout the entire simulation.

In order to generate initial conditions for a QM/MM simulation, you may need to go beyond the scripts provided with Newton-X, and do some coding on your own. Below I'm pasting an excerpt from the SI for a paper that I wrote which due to be published in JPCA. I'm not necessarily saying this is the best way to do it. But it should give you a general idea of what would need to be done.



The initial conditions for the NAMD simulations were defined in such a way as to
reproduce, insofar as possible, the experimental conditions used by Druzhinin and
co-workers. Setting up the simulations was a fairly complex, multi-stage procedure. This
was in part because the system under study is described using the hybrid QM/MM method,
such that the solute and the solvent are treated on a different footing. I followed the de
facto standard practice in NAMD simulations of photoinduced processes in the solution
phase, in which ZPVE energy is imparted on the solute molecule, but not on the
solvent.

The first step was to generate the nuclear positions and velocities of the fluorazene
molecule. To this end, the ground-state minimum-energy geometry of fluorazene was
optimized at the CAM-B3LYP level, and its normal modes were calculated numerically.
Afterwards, 1000 phase space points (which is to say, 1000 sets of nuclear positions and
velocities) were sampled from the Wigner distribution. As pointed out in Refs. 47−49,
semiclassical simulations are prone to an artificial leakage of zero-point vibrational energy
from the stretching modes of hydrogen-heavy atom bonds to other vibrational modes. In
order to alleviate this problem, when generating the Wigner distribution, I froze the nine
vibrational modes with the highest frequencies, which correspond to C–H stretching modes.
At the second stage of the procedure, the solvent was equilibrated with the solute. This
was achieved by placing the fluorazene molecule (the ground-state minimum-energy
geometry) at the center of a 500-molecule ACN nanodroplet. Afterwards, the solvent was
equilibrated with the solute by propagating a molecular dynamics trajectory in the
canonical (N V T ) ensemble.

During the equilibration, the geometry of the fluorazene molecule was kept frozen, and
the potential energy surface of the entire system was calculated at the MM level of theory.
The charge distribution of the solute was represented by a set of point charges placed at the
nuclear positions. These charges were obtained by fitting them to the electrostatic potential
generated by the fluorazene molecule. The fit was performed with the Merz-Singh-Kollman
(MKS) scheme implemented in Gaussian 16. The temperature of the solvent was
maintained at 298 K by imposing the Langevin thermostat with a friction coefficient of
γ = 1 ps−1 .

The Langevin equation of motion was propagated with the use of the
Brünger-Brooks-Karplus (BBK) integrator 51 with a time step of 0.5 fs. The equilibration
period lasted for 100 ps.

Following the equilibration period, 1000 solvent configurations were sampled from the
molecular dynamics simulation at intervals of 2 ps, and were combined with the solute
geometries generated previously by sampling from the Wigner distribution. In effect, this
means that the solvent was equilibrated with the minimum-energy geometry of fluorazene,
and not with the individual geometries sampled from the Wigner distribution, which
include vibrational displacements from the minimum-energy geometry. This approximation
is acceptable for a relatively rigid molecule such as fluorazene.

The 1000 combined solvent-solute geometries were used as the basis for the simulation
of the photoabsorption spectrum of fluorazene with the nuclear ensemble method. 52,53 In
the calculation of the spectrum, a Gaussian line shape function was used with a variance
of σ 2 = 0.04 eV2 . In Figure S9, the resulting spectrum is compared to the experimental
spectrum of fluorazene in ACN solution. It can be seen that, in the simulated spectrum,
the first photoabsorption band is blue-shifted by roughly 0.8 eV with respect to the observed
spectrum. This is in part because the calculation was performed with the relatively small
def2-SVP basis set.

The fact that the calculated spectrum is blue-shifted with respect to experiment has
implications for the choice of initial conditions. The TA experiments reported in Ref. 27
were performed with a pump pulse wavelength of 290 nm, which corresponds to a photon
energy of 4.28 eV. (Some experiments were also performed with photoexcitation at 266 nm,
but their results were reported to be similar to those with the longer wavelength. 27 ) Because
the TDDFT method overestimates the energy range of the first photoabsorption band of
fluorazene, I must also assume a higher photon energy of the pump pulse. Accordingly,
I decided to set the pump pulse photon energy to 5.1±0.1 eV.

Phase space points (sets of nuclear positions and velocities along with an initial adiabatic
state from S1 to S4 ) with excitation energies in the range 5.1±0.1 eV were sampled with
probabilities proportional to their oscillator strengths. A total of 60 phase space points were
selected in this way, and each was used as the starting point for a single NAMD trajectory.
The simulations were propagated for a time period of 1.5 ps.

[Michal Kochman]

>> I chose the diffuse function, because they are describing the molecule best, but I did not know it could interfere with the point charges. They are included in my calculation.

Please see the paper "Avoiding Electron Spill-Out in QM/MM Calculations on Excited States with Simple Pseudopotentials" by Khah et al., and "Benchmark Studies on the Building Blocks of DNA. 2. Effect of Biological Environment on the Electronic Excitation Spectrum of Nucleobases" by Szalay et al. The first paper deals directly with the QM/MM method. The second is relevant because it shows that solvation causes a destabilization of the diffuse, Rydberg-type states relative to compact, valence-type states. Often you don't absolutely need diffuse basis functions when modeling photochemistry in the condensed phase.

[Anna U]

Thanks a lot for your detailed answer. I now know what to work on.

I got a single question left: In the plot you provided, is there the energy of the QM system plottet or the energy of the whole system?

If I plot the QM energy only, it looks much better (energies fluctuation is just some eV). This shows as well, that the solvent is not well equilibrated...