From: Junwoong Yoon (junwoony_at_andrew.cmu.edu)
Date: Mon Jul 23 2018 - 23:18:58 CDT
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> 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>
> wrote:
>
>> Hi Jun,
>>
>> Read the colvar docs again. See http://www.ks.uiuc.edu/
>> Research/namd/2.12/ug/node58.html#SECTION000135400000000000000
>> 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
>> > 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/wham
>>> <https://na01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fmembrane.urmc.rochester.edu%2Fcontent%2Fwham_&data=02%7C01%7CJoshua.Vermaas%40nrel.gov%7Ce3890a31c6464102ad1308d5edb70297%7Ca0f29d7e28cd4f5484427885aee7c080%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> 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> 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 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 : Tue Dec 31 2019 - 23:20:06 CST