How to get 2D free energy landscape with umbrella sampling ?

Dear lammps users

I am trying to get a 2D free energy landscape by umbrella sampling in lammps colvars.For example the free energy change along with (x=0,1,2…10) and (y=0,1,2…10).
In the simulation, two colvars components and harmonic biases potentials are defined:

colvar { # x

name one


}

colvar{ #y
name two

}

harmonic {

name h_1
colvars one

}

harmonic {

name h_2
colvars two

}

But I found that the two components change simultaneously, which means that only the free energies at points like (0,0), (1,1), …(10,10) could be computed. While the ones like (0,1),(0,2)…(0,10) could not. And My question is how could I define the biases potentials and get the 2D free energy landscape with umbrella sampling?
Thanks in advance for any suggestion.

Best,
Zhiqiang

Hi Jhon, first of all, please do realize all other free energy methods involve 1) a bias to enhance statistical sampling and 2) a free-energy estimator. For umbrella sampling the bias is in this case is always the harmonic potential, but the free energy estimator varies (WHAM, M-BAR, etc) and is typically a post-processing phase (i.e. you won’t get a PMF from LAMMPS + Colvars directly). In any case, you can use the colvars.traj files and the harmonic potential parameters (defined by you) as input for the post-processing step.

You are sampling a single “point” in the 2D space, defined by the fixed centers of the two harmonic restraints. You are responsible for changing the centers of the two restraints in the input of each umbrella sampling window/replica, and process the aggregated data using WHAM or other tool. You will need multiple input files for this, of course.

Lastly, please consider using a single 2D potential:

harmonic {

name umb

colvars one two

centers x_0 y_0

forceConstant K

}

which is of course equivalent to two independent harmonic restraints, but is easier to maintain when you change x_0 and y_0.

Giacomo

Dear Fiorin,

Thanks for your kindly suggestions.
Thus, my understanding is that if the potential in one simulation:

harmonic {

name umb

colvars one two

centers x_0 y_0

forceConstant K

}

used and x_0, y_0 change from 0 to 10 respectively. What I could get after processed by WHAM is Series of points in the reaction path defined by x_0, y_0 in the simulation. To get the 2D free energy landscape, series of simulations are needed for WHAM to post processing.

For another thing, apart from umbrella sampling, to get a 2D free energy landscape in problems like the rotation of transmembrane protein(all-atom simulations), which method do you recommend(like ABF or metadynamics)?

Have a nice day.
Thanks
Zhiqiang

Dear Fiorin,

Thanks for your kindly suggestions.
Thus, my understanding is that if the potential in one simulation:

harmonic {
    name umb
    colvars one two
    centers x_0 y_0
    forceConstant K
}
used and x_0, y_0 change from 0 to 10 respectively. What I could get
after processed by WHAM is Series of points in the reaction path defined
by x_0, y_0 in the simulation. To get the 2D free energy landscape,
series of simulations are needed for WHAM to post processing.

For another thing, apart from umbrella sampling, to get a 2D free energy
landscape in problems like the rotation of transmembrane protein(all-atom
simulations), which method do you recommend(like ABF or metadynamics)?

I can't answer on that without knowing the specifics and/or having done
some test runs on the TM protein of interest. I can only say that each
method has either (a) a different way of biasing the sampling or (b) a
different way of measuring the free energy changes.

Since it's by all means a complex problem, I would look first at the
literature on similar systems, and pick a collective variable that will
fit. The choice of method/algorithm is secondary, particularly with
Colvars that allows you to test many ones out.

Giacomo

Dear Fiorin,

Thanks for your kindly suggestions.
Thus, my understanding is that if the potential in one simulation:

harmonic {
    name umb
    colvars one two
    centers x_0 y_0
    forceConstant K
}
used and x_0, y_0 change from 0 to 10 respectively. What I could get after
processed by WHAM is Series of points in the reaction path defined by x_0,
y_0 in the simulation. To get the 2D free energy landscape, series of
simulations are needed for WHAM to post processing.

For another thing, apart from umbrella sampling, to get a 2D free energy
landscape in problems like the rotation of transmembrane protein(all-atom
simulations), which method do you recommend(like ABF or metadynamics)?

I can't answer on that without knowing the specifics and/or having done some
test runs on the TM protein of interest. I can only say that each method
has either (a) a different way of biasing the sampling or (b) a different
way of measuring the free energy changes.

Since it's by all means a complex problem, I would look first at the
literature on similar systems, and pick a collective variable that will fit.
The choice of method/algorithm is secondary, particularly with Colvars that
allows you to test many ones out.

in addition to that, it is always a very good idea to devise a (much)
simpler test system to practice and explore how to effectively and
correctly compute free energies with that. there are numerous examples
discussed in the relevant literature. it is good practice to start by
reproducing some of those and then work your way up. computing free
energy profiles has two main challenges: managing the technology well,
and understanding the dynamics of your system well to choose the right
methodology and settings. it'll be much easier to address these two
challenges one after the other rather than both at the same time.

axel.