From: Marcelo C. R. Melo (melomcr_at_gmail.com)
Date: Mon Oct 17 2022 - 11:42:47 CDT
Hi Oleksii,
You understood it correctly, but I would caution you regarding the use of
different QM packages. ORCA receives a cloud of classical point charges
surrounding the QM region in order to run its calculations in electrostatic
embedding, and it calculates electrostatic interactions between nearby
classical charges internally (NAMD does not see that calculation).
Other packages (such as MOPAC, which we also used in that manuscript you
cited) do not calculate and return the results of electrostatic interaction
calculations. Therefore, the determination of E(QM/MM) will take different
forms depending on the QM package used.
That said, since you are using ORCA, you will get the energy from
electrostatic interactions between classical charges and QM charges from
the ORCA output. In this case, QMENERGY will include E(QM) and a portion
of E(QM/MM). NAMD will calculate long range electrostatic interactions (if
you have PME turned on) and (as you noted) van der Waals between the
classical and quantum part, which is the rest of E(QM/MM). Therefore, your
job is a bit more complicated because you will have to extract some of the
electrostatic energy from ORCA's output. I am not an ORCA expert, but I
believe you can get this QM-MM energy due to electrostatic embedding from
ORCA's output.
As for separating E(MM) and E(QM/MM), you would have to select an arbitrary
collection of atoms (the ones composing the QM region) and calculate
the van der Waals between those and the surrounding classical atoms, up to
the non-bonded cutoff. I believe VMD has a plugin that can help with that,
but others will be of more help than me here. NAMD does not treat QM-MM van
der Waals interactions differently from MM-MM van der Waals interactions,
so any traditional way to compute this would work.
Finally, the last (and usually complicated) issue with QM/MM simulations is
long range electrostatics. Using ORCA, part of the long range
electrostatics will be calculated through its electrostatic embedding, and
should be available from its output. However, NAMD's PME does calculate
long range interactions between QM atoms and MM atoms that are much more
distant than the electrostatic embedding cutoff. I am not sure how to
calculate that particular fraction of E(QM/MM) in isolation after the
simulation. That is not a fraction of energy we tend to save. To determine
this after the simulation, you would have to keep the partial charges
placed on QM atoms throughout the simulation (something you can save using
NAMD keywords when initiating the simulation), and then back calculate the
interactions using timesteps of the DCD. This takes advantage of the fact
that classical point charges do not change during a simulation, only QM
partial charges do.
Let me know if this helps,
Marcelo
On Fri, Oct 14, 2022 at 6:04 AM Oleksii Zdorevskyi <zdorevskyi_at_bitp.kiev.ua>
wrote:
> Dear NAMD community,
>
> I have a question regarding the calculation of interaction energies in
> QM/MM simulations. I have browsed through previous NAMD-l discussions,
> however, did not find any similar topics.
>
> According to the NAMD QM/MM paper (Melo, Marcelo CR, et al. "NAMD goes
> quantum: an integrative suite for hybrid simulations." Nature methods 15.5
> (2018): 351-354.), the total *potential* energy of the system is computed
> in the following way:
>
> E(total) = E(MM) + E(QM) + E(QM/MM),
>
> where E(MM) is the energy of the classical part, E(QM) - is the
> single-point energy calculated by the QC software, and E(QM/MM) is the
> interaction energy between the classical and quantum regions.
>
> Now, I need to plot these 3 components as separate time series.
>
> When I browse through the log file, I can see the "QMENERGY" keyword,
> which corresponds to the E(QM) (actually, I checked it by comparing this
> value with the corresponding single-point energy from ORCA output).
>
> Then, I can extract the energy called "POTENTIAL" from the same log file.
> If I understand correctly, this energy contains all the interactions in
> the system, including van der Waals between the classical and quantum
> part, as well as Coulomb interaction between classical charges, and
> single-point charges derived from the previous QM step (when we use
> "electrostatic embedding"). Does this energy describe E(tot), listed in
> the formula above, correct?
>
> If I subtract the "QMENERGY" (E(QM)) from the "POTENTIAL" (E(tot)), it
> will give me the following:
> E(MM)+E(QM/MM),
>
> so E(MM) and E(QM/MM) can be only obtained as a sum.
>
> However, I need to plot separate contributions for E(MM), E(QM), and
> E(QM/MM).
>
> Are there any means to extract them from NAMD output?
>
> I will appreciate your help.
>
> Best regards,
> Oleksii
>
>
>
This archive was generated by hypermail 2.1.6 : Tue Dec 13 2022 - 14:32:44 CST