How gaussians are added to periodic CVs?

[psn417]

Hello everyone, I am trying to learn the Well-Tempered MetaD. I use the alanine dipeptide example. And I write a python script to add the gaussians:

def wrapped_gaussian(x, y, phi, psi, sigma_phi, sigma_psi, height):
    ret = 0
    x_phi = []
    order = [-2, -1, 0, 1, 2]
    for ii in order:
        x_phi.append((phi-x).copy()-ii*np.pi*2)
    x_phi = np.array(x_phi)
    y_psi = []
    for ii in order:
        y_psi.append((psi-y).copy()-ii*np.pi*2)
    y_psi = np.array(y_psi)
    for a in x_phi:
        for b in y_psi:
           ret += np.sum(height * np.exp(-a**2/(2*sigma_phi**2) - b**2/(2*sigma_psi**2)))
    return ret

# read hills
def get_bias():
    hills = np.loadtxt("HILLS").T
    phi = hills[1]
    psi = hills[2]
    sigma_phi = hills[3]
    sigma_psi = hills[4]
    height = hills[5]
    def bias(positions):
        ret = []
        """
        the bias is periodic !!!
        """
        x = positions[0] / 180 * np.pi
        y = positions[1] / 180 * np.pi
        for xx, yy in zip(x, y):
            ret.append(wrapped_gaussian(xx, yy, phi, psi, sigma_phi, sigma_psi, height))
        return np.array(ret)
    return bias

bias = get_bias()

phi, psi = np.mgrid[-180:180:3, -180:180:3]
positions = np.vstack([phi.ravel(), psi.ravel()])
bias_potential = np.reshape(bias(positions).T, phi.shape)

plt.contourf(phi, psi, -bias_potential, levels=50, cmap="Spectral")
plt.colorbar()
plt.title("bias potential")

But I find the result is different (almost identical, but slightly different) from what I get from plumed sum_hills --hills HILLS --bin 120,120. I wonder how the sum_hill add periodic CVs? Does it uses wrapped gaussians?

[psn417]

I read the source code and find the differeces:

1. The gaussians are stretched. They are reduced to zero at a certain cutoff.
2. Only the nearest distance are calculated in periodic CV.

[psn417]

The following code can produce the exactly same result of sum_hills:

def stretched_gaussian(x, y, phi, psi, sigma_phi, sigma_psi, height):
    d1 = (phi-x)/(2*np.pi)
    d1[d1<-0.5] += 1
    d1[d1>0.5] -= 1
    d1 *= 2*np.pi

    d2 = (psi-y)/(2*np.pi)
    d2[d2<-0.5] += 1
    d2[d2>0.5] -= 1
    d2 *= 2*np.pi

    dp2 = d1**2 / (2*sigma_phi**2) + d2**2 / (2*sigma_psi**2)
    tmp = np.zeros(phi.shape)
    # stretch it at cutoff=6.25
    tmp[dp2<6.25] = height[dp2<6.25] * (np.exp(-dp2[dp2<6.25]) * 1.00193418799744762399 - 0.00193418799744762399)
    return tmp.sum()

# read hills
def get_bias():
    hills = np.loadtxt("HILLS").T
    phi = hills[1]
    psi = hills[2]
    sigma_phi = hills[3]
    sigma_psi = hills[4]
    height = hills[5]
    def bias(positions):
        ret = []
        x = positions[0] / 180 * np.pi
        y = positions[1] / 180 * np.pi
        for xx, yy in zip(x, y):
            ret.append(stretched_gaussian(xx, yy, phi, psi, sigma_phi, sigma_psi, height))
        return np.array(ret)
    return bias

bias = get_bias()

phi, psi = np.mgrid[-180:180:3, -180:180:3]
positions = np.vstack([phi.ravel(), psi.ravel()])
bias_potential = np.reshape(bias(positions).T, phi.shape)

plt.contourf(phi, psi, -bias_potential, levels=50, cmap="Spectral")
plt.colorbar()
plt.title("bias potential")

[Cloudac7]

Hello, I am also learning about how to do the sum_hill and what made me curious about, as a result, is why dp2cutoff is chosen to be 6.25? Is there any mathematical reason about the value?

[psn417]

Hello, I am also learning about how to do the sum_hill and what made me curious about, as a result, is why dp2cutoff is chosen to be 6.25? Is there any mathematical reason about the value?

Hi, I think that why they choose 6.25 is not important. The bias doesn't need to be very pericise. We can use reweighting to get the free energy surface. If we use reweighting, sum_hill is completely useless. We just need the metad.bias value of each moment.

You can read the paper "Metadynamics with Adaptive Gaussians" for more detail.

[GiovanniBussi]

Notice that sum_hills uses the same convention. Using reweighting might help if the features of the landscape are narrower than the Gaussians, or if you want to analyze another variable. Besides this, the choice of a Gaussian or of another (similar) kernel is pretty irrelevant, but it's important that you use the same convention in sum_hills and in METAD.

[Cloudac7]

Hello, I am also learning about how to do the sum_hill and what made me curious about, as a result, is why dp2cutoff is chosen to be 6.25? Is there any mathematical reason about the value?

Hi, I think that why they choose 6.25 is not important. The bias doesn't need to be very pericise. We can use reweighting to get the free energy surface. If we use reweighting, sum_hill is completely useless. We just need the metad.bias value of each moment.

You can read the paper "Metadynamics with Adaptive Gaussians" for more detail.

Notice that sum_hills uses the same convention. Using reweighting might help if the features of the landscape are narrower than the Gaussians, or if you want to analyze another variable. Besides this, the choice of a Gaussian or of another (similar) kernel is pretty irrelevant, but it's important that you use the same convention in sum_hills and in METAD.

Thank you kindly for reply.