From: Oleksii Zdorevskyi (zdorevskyi_at_bitp.kiev.ua)
Date: Fri Oct 21 2022 - 05:29:37 CDT
Dear Marcelo,
Thank you very much for your comprehensive reply.
According to your suggestions, I have changed my algorithm. Now, I am
doing it in the following way:
1. MM energy - E(MM). To calculate it, I am using the NAMDEnergy VMD
plugin. I create a selection of all the atoms in the system excluding the
QM region and calculate the non-bonded interactions inside this selection.
According to the NAMD Manual, these interactions include van der Waals and
Electrostatics. My simulation is performed on a protein in vacuum with
only structurally resolved water, so no PME is used.
2. E(QM/MM)=E(QM/MM:namd)+E(QM/MM:orca). To derive the part of this energy
treated by NAMD - E(QM/MM:namd), I first apply the same procedure as in
step 1, but for the whole system. In this way, I derive E(tot:namd). Then,
the desired QM/MM:namd energy will be:
E(QM/MM:namd) = E(tot:namd) - E(MM)
3. The remaining parts I need, are E(QM/MM:orca) - the part of QM/MM
energy calculated by ORCA, and E(QM) - the pure energy of the QM region.
I have a question regarding these parts. The ORCA output is being written
to the "base dir", however, these are temporary files which are
overwritten every step. Therefore, after my simulation is done, I can see
the ORCA output files only for the last frame.
The only thing which is being recorded is the QMENERGY in the NAMD log
file, which, as you mentioned, includes both terms:
QMENERGY = E(QM) + E(QM/MM:orca).
Is there any way to record the ORCA output at every frame, to enable
separating those two (E(QM), and E(QM/MM:orca)) after the simulation?
Thanks.
Best regards,
Oleksii
> 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