Strategies for deciding simulation parameters for 2D alchemical metadynamics

[wehs7661]

This issue is for the discussion of the second problem, which is mainly about the strategies for deciding simulation parameters for 2D alchemical metadynamics for the CB8-G3 system.

Briefly, the questions described in the corresponding jupyter notebook (any commit after 545046e) can be summarized as follows:

• Are the strategies I used for deciding metadynamics parameters appropriate?
• Can the bias be regarded as quasi-stationary if the Gaussian height in the unphysical region is still very high?
• What are the general suggestions you would give about deciding the parameters for a wall potential? How large should the energy penalty typically be?
• What should I modify in my protocol for free energy calculations to consider the reweighting of the wall potential?
• Is there a way to prohibit the exploration of regions such as $\lambda=40$ and $n=0$?

[wehs7661]

Regarding the methods for restricting the sampling region in the CV space, here is one idea I proposed when brainstorming the solution:

Say that we are interested in the sampling problem of a host-guest binding complex characteristic of two CVs: lambda and the number of water molecules in the binding cavity (N). We could first run a vanilla simulation for the coupled and the uncoupled states, respectively, and use plumed driver to figure out the values that are physically meaningful at a certain lambda value. For example, say that in the coupled state (lambda = 0), the range of N is from 2 to 5, whereas for the uncoupled state, N ranges from 10 to 15. Then, we could consider restricting the sampling space as shown in the figure below, where the shaded region is the only region that the system is allowed to sample.

With UPPER_WALLS and LOWER_WALLS , we could define walls like V_b1 and V_b2. To express V_b3 and V_b4, which define the boundary of the allowed region, we could use CUSTOM to define CV1 = N - 8 * lambda and CV2 = 4N + 5 * lambda, so instead of having V_1 = k(x-2)^{2} and V_2=(x-15)^{2} that defines a rectangular region, we could use V_1 =k(CV1 + 2)^2 and V2=k(CV2-7)^2 to define a trapezoidal region. I haven't tried out this approach though.

Certainly, in many cases, we don't know the right pathway and a trapezoidal region might oversimplify the problem, but I was trying to figure out at least a starting point of methods that use exiting PLUMED tools to sort of address our needs.