Re: [NAMD] How to decompose REUS results into electrostatic and van der Waals contributions

From: Jeff Comer (jeffcomer_at_gmail.com)
Date: Wed Apr 03 2019 - 15:10:03 CDT

Here's the approach that I use, which is similar to what JC suggested,
except that I don't use the pairInteraction command. However, there
are multiple reasons why you might not be able to do this with the
simulations that you have. Also, you may find that these energies
aren't that useful because they are (1) extremely noisy for large
explicit-solvent systems and (2) hard to interpret owing to
compensation between different energies.

You can calculate U_elec(ξ) and U_vdW(ξ), where ξ is the transition
coordinate for your umbrella sampling. I've done what you're
suggesting from ABF simulations in the following paper:

https://doi.org/10.1021/acs.jctc.8b00830

Here's what I do:
(1) Extract U_elec and U_vdW from the log file. I use awk for this:
awk -v OFMT=%.15g '/^ENERGY: / {print 2e-6*$2,$7}' ../sim.log > sim.elec.log
awk -v OFMT=%.15g '/^ENERGY: / {print 2e-6*$2,$8}' ../sim.log > sim.vdw.log
(2) Extract the transition coordinate ξ at the same simulation steps
from the .colvars.traj file.
(3) Paste together the transition coordinate ξ(step) and the energy
U_elec(step), and so forth.
(4) Calculate the average energies U_elec and U_vdW in bins along ξ to
estimate U_elec(ξ) and U_vdW(ξ).

The noisiness of the energies U_elec and U_vdW will limit your
statistics a lot. To get the best statistics, I would suggest that you
run with "outputEnergies 500" and "colvarsTrajFrequency 500". Smaller
values don't help you much owing to correlation between frames.

If you haven't done this in your simulations, you can recalculate the
energies and transition coordinates from the dcd file, assuming you
are writing a lot of frames (probably dcdFreq 500 or dcdFreq 1000).
You can run the dcd files through NAMD to obtain the .log and
.colvars.traj files:

outputEnergies 1
set ts 0
foreach dcd $dcdList {
    coorfile open dcd $dcd
    while { ![coorfile read] } {
    firstTimestep $ts
    run 0
    incr ts 1
    }
    coorfile close
}

***
I would guess that doing all this might be impossible for the
simulations that you have. However, I would recommend running new
simulations in the end-point states (usually bound and unbound) with
"outputEnergies 500" and "colvarsTrajFrequency 500". This will allow
you to separate different components at least for the bound and
unbound states. We used this approach to separate the ΔG_bind into
ΔU_bind and –TΔS_bind contributions (
https://doi.org/10.3390/toxins10100397 ). From the umbrella sampling,
you know where the free energy minimum is (ξ_min) and the free energy
there ΔG(ξ_min), so you can extract frames near the free energy
minimum and begin your bound simulations from them. You should check
that different starting coordinates give the same results. You can use
flat bottom restraints (Colvars harmonicWalls) to keep the system near
ξ_min. You can do something similar for the unbound state (ξ_far). You
can use the approach above to get <U(ξ_min)> (or even <U_elec(ξ_min)>
and <U_vdw(ξ_min)>) and <U(ξ_far)>. Even in this case, the statistics
will be difficult, so make sure you calculate the error bars.

Jeff

–––––––––––––––––––––––––––––––––––———————
Jeffrey Comer, PhD
Assistant Professor
Institute of Computational Comparative Medicine
Nanotechnology Innovation Center of Kansas State
Kansas State University
Office: P-213 Mosier Hall
Phone: 785-532-6311
Website: http://jeffcomer.us

On Tue, Apr 2, 2019 at 7:26 PM Canal de Sebassen
<thecromicusproductions_at_gmail.com> wrote:
>
> Hello,
>
> I'm post-processing some long REUS simulations and would like to decompose
> the free energy I obtain from WHAM into its vdW and electrostatic components. Is there any way to do this?
>
> Thanks a lot,
>
> Sebastian

This archive was generated by hypermail 2.1.6 : Thu Dec 05 2019 - 23:20:38 CST