RE: Umbrella Sampling input parameters

From: Vermaas, Joshua (Joshua.Vermaas_at_nrel.gov)
Date: Tue Jul 24 2018 - 09:45:08 CDT

Hi jun,

Are the higher force constants being passed along to WHAM? How much simulation time is behind each window? Are your overlaps still ok with the higher force constant? What was the WHAM termination tolerance? You are right to assume that the PMFs should be the same, but for a variety of reasons, WHAM might get you unexpected (and incorrect) results.

Higher force constants tend to reduce the variance within a window, reducing overlap between adjacent umbrellas.

-Josh

On 2018-07-23 22:19:05-06:00 Junwoong Yoon wrote:

Hi Josh,
I'm running two series of umbrella samplings with the same coordinate files but with different force constants.
It seems like their pmf values are different. (higher force constant yields higher pmf values)
I expected that the pmf values would be similar even with different forces constants because WHAM will unbias them.
Also how the force constants affect the overlaps between windows?
Thank you
Jun

On Thu, Jul 19, 2018 at 7:01 PM, Jérôme Hénin <jerome.henin_at_ibpc.fr<mailto:jerome.henin_at_ibpc.fr>> wrote:
Indeed CHARMM's "harmonic" constants are defined differently from those of the rest of the world. NAMD's implementation is faithful to CHARMM, for understandable reasons. Colvars isn't because it is more generic. Unfortunately we must live with these historical oddities.
Jerome

On Thu, 19 Jul 2018 at 22:56, Vermaas, Joshua <Joshua.Vermaas_at_nrel.gov<mailto:Joshua.Vermaas_at_nrel.gov>> wrote:
Hi Jun,

Read the colvar docs again. See http://www.ks.uiuc.edu/Research/namd/2.12/ug/node58.html#SECTION000135400000000000000.Vermaas%40nrel.gov%7C94df613a9ef34b0b82f808d5f11c95f0%7Ca0f29d7e28cd4f5484427885aee7c080%7C0%7C0%7C636680027451421563&sdata=R2QkFqHwvCYRzDtTvrc%2B37GcThFf3%2B3m8k%2Br3bIoBZ0%3D&reserved=0>
This potential form (k/2 (x-x0)^2) matches the form many WHAM codes use (including Grossfields if you read the docs), so the adjustments are minimal.

-Josh

On 2018-07-19 14:34:26-06:00 Junwoong Yoon wrote:

Hi Josh,
I have a confusion about what the force constant in WHAM should be.
In NAMD, the the harmonic has k(x-x0)^2, on the other hand, the harmonic is 1/2*k(x-x0)^2 in WHAM.
So if I set my forceConstant to 2.5 kcal/mol/A^2 with 1A width in NAMD, should I input 5 kcal/mol in WHAM to generate PMF?
Thank you
Jun

On Mon, Jul 16, 2018 at 9:48 PM, Vermaas, Joshua <Joshua.Vermaas_at_nrel.gov<mailto:Joshua.Vermaas_at_nrel.gov>> wrote:
Hi Jun,

I tend to do alot of post-processing, and to be honest I'm confused as to how you got a PMF without knowing what is in your .traj files. Here is my workflow:
Step #1, look at your .traj files. The first row labels the contents. For me, the first column is the timestep, and the second is the value of the colvar, and the third is a printout of the value for the center (I turn on "outputCenters" so I can trace the path of my sampling in replica exchange umbrella sampling. This is not required for conventional umbrella sampling).
Step #2, plot histograms/exchanges based on the columns identified earlier. Use whatever tools you like to make the histograms and identify overlap.
Step #3, determine PMF with a WHAM code. Alan Grossfield's (http://membrane.urmc.rochester.edu/content/whame7c080%7C0%7C0%7C636676292658097288&sdata=6f4B6SX0p1Jc%2FmXRC%2Fq6AxYni3CAPfhZbGnriAbid2Y%3D&reserved=0>) is the easiest to use because of its documentation, but once you understand the math involved, you can roll your own or adapt another one that suits your needs.

Since WHAM is iterative, in my experience an infinity anywhere quickly propagates, so I'm not sure how/why it didn't for you.

-Josh

On 2018-07-16 19:11:18-06:00 Junwoong Yoon wrote:

Hi Josh,
I could do umbrella sampling and WHAM analysis but I still have some things needed to be clear.
1) I could get PMF from WHAM. However, I don't know how I can make histogram plot for each window and see there are sufficient overlaps. Which file should I use for it? or How can I generate that file. From WHAM, I only got the PMF and probability data, and from Umbrella Sampling, I only got .traj file.
2) I have "inf" free energy for one reaction coordinate. (everything else looks reasonable). Is this because of not sufficient overlap?
I really appreciate your taking time to help me and others.
Jun

On Fri, Jul 13, 2018 at 9:30 PM, Vermaas, Joshua <Joshua.Vermaas_at_nrel.gov<mailto:Joshua.Vermaas_at_nrel.gov>> wrote:
Hi Jun,

1.) I was being a bit careless with my nomenclature. This is not SMD as in the SMD module of NAMD, but instead using a moving umbrella to do SMD. To do SMD, your centers should be the initial value, and the targetCenters would be your desired value at the end of the simulation. In this case, if I were pulling a molecule down, my "ref" would be a selection at the center of the membrane, and the "main" selection would be my molecule. So long as targetCenters < centers and the force constant and targetNumSteps are reasonable.

2.) Math it out. I like adjacent umbrellas/windows to cross at a height of ~0.5kcal/mol. If your width is 1, this comes down to solving 0.5=k (dx/2)**2 for your value of dx.

-Josh

On 2018-07-12 16:45:25-06:00 Junwoong Yoon wrote:

Hi Josh,
Thank you so much for your explanation. It was very helpful to prepare my umbrella sampling strategy.
But I still have couple questions that bother me.
1) For SMD prior to umbrella sampling, can I use SMDDir 0 0 -1 if I am trying to pull a molecule down along the z-direction?
Or should I still need to define a normalized direction between fixed atoms and SMD atoms for SMDDir?
How will this be different if I'm calculating distance in x-y plane? For example, pulling two molecules on a surface while their heights(z-dir value) are fixed?
2) To make sure neighboring windows to overlap, what forceConstant should I use if my width is 0.1, and each window is equally spaced by 1A?
    I used 0.025 based on the tutorial, but I want to clarify how does 0.025 scaled forceConstant allow windows to overlap.
My colvars file looks like,
colvar {
    width 0.1
    distance Z { }
    harmonic {
        centers $CENTER
        forceConstant 0.025
Thank you
- Jun

On Tue, Jul 10, 2018 at 6:40 AM, Vermaas, Joshua <Joshua.Vermaas_at_nrel.gov<mailto:Joshua.Vermaas_at_nrel.gov>> wrote:
Hi Jun,

There are two ways of setting up umbrella sampling with colvars. The most common approach is to use an initial SMD pull to cover your reaction coordinate space (center = initial value, targetCenter = final value), and then pick frames from that pull to seed your umbrellas. In that case, each umbrella would have a center, with no defined targetCenter (since the umbrella isn't moving). The other alternative combines your SMD and umbrella sampling together, and sequentially moves umbrellas. In this case, you would use the targetNumStages parameter to move your umbrellas sequentially along your chosen reaction coordinate. This has two drawbacks: 1.) You can't easily make your simulations longer if you are not converged, since the same trajectory is used for all states. 2.) Your umbrella sampling calculation becomes serial, as each window depends on the result of the last one. Since oftentimes umbrella sampling has like a microsecond of MD behind it, this would take prohibitively long for real systems. So if it were me, I'd use the more common approach.

Another piece of advice is to just set your width to 1 unless you are combining dissimilar collective variables together. The applied force constant is unrelated to the upper or lowerboundary parameters. It is instead determined by the forceconstant parameter applied to the harmonic and the width of the collective variable, as described in the documentation of the forceconstant parameter. This also clarifies that you indeed should NOT set upper/lower limits. This is effectively already done by the harmonic potentials, and adding extra biases will skew WHAM analysis.

-Josh

Example umbrella sampling configuration file:

#############################################################
## Collective Variables ##
#############################################################
# Global parameters
colvarsTrajFrequency 50
colvarsTrajAppend off
analysis off
colvar {
        name contactcount
#Here I don't follow my own advice and set the width to be not 1. This was so that the whole reaction coordinate could be mapped to the range (0,1)
        width 977.429
        coordNum {
                expNumer 2
                expDenom 24
                cutoff 10.0
                group1 {
                        atomsFile beta.pdb
                        atomsCol B
                        atomsColValue 1
                }
                group2 {
                        atomsFile beta.pdb
                        atomsCol B
                        atomsColValue 2
                }
        }
}
harmonic {
    name contactpull
    colvars contactcount
    centers 684.20000973
#Technically, I didn't need to specify targetCenters, since the default is for targetCenters=centers
    targetCenters 684.200010
    forceConstant 2000.000000
    targetNumSteps 100
}

On 2018-07-10 01:00:28-06:00 owner-namd-l_at_ks.uiuc.edu<mailto:owner-namd-l_at_ks.uiuc.edu> wrote:

Hi all,
I want to implement Umbrella Sampling and questions about input parameters.
My colvars is a distance between two molecules.
If my initial distance is 11, and I want to make free energy profile along the distance in range of from 2 to 30 with width of 0.1.
Then what should be my "centers" and "targetCenters"?
If I set "centers" to 11 and "targetCenters" to 31, then will I not be able to explore the region between 2 and 11?
And will the force constant be [upperboundary - lowerboundary]*width^2 ?
I heard that I should not define lower/upperBoundary, or lower/uppderWallConstant for Umbrella Sampling, is it correct?
Thank you
Jun

This archive was generated by hypermail 2.1.6 : Sun Sep 15 2019 - 23:19:53 CDT