Plot a well known double well surface

Question

I am trying to plot a well known energy landscape from physics, the Muller Brown potential. Taken from the literature (https://arxiv.org/abs/1701.01241):

With my contour plot, however, I cannot see the two wells, it looks as if it was just a single Gaussian. Am I doing something wrong?

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mp

# Parameters of Muller Brown Potential (cf. Bonfati & Cob 2017)
A = np.array([-200, -100, -170, 15])
a = np.array([-1, -1, -6.5, 0.7])
b = np.array([0, 0, 11, 0.6])
c = np.array([-10, -10, -6.5, 0.7])
x_m = np.array([1, 0, -0.5, -1])
y_m = np.array([0, 0.5, 1.5, 1])

x = np.linspace(-1.5, 1, 1000)
y = np.linspace(-0.5, 2, 1000)
XX, YY = np.meshgrid(x, y)


Z =A[0]*np.exp( a[0]*(XX-x_m[0])**2 + b[0]*(XX-x_m[0])*(YY-y_m[0]) + c[0]*(YY-y_m[0])**2 )
+  A[1]*np.exp( a[1]*(XX-x_m[1])**2 + b[1]*(XX-x_m[1])*(YY-y_m[1]) + c[1]*(YY-y_m[1])**2 )
+  A[2]*np.exp( a[2]*(XX-x_m[2])**2 + b[2]*(XX-x_m[2])*(YY-y_m[2]) + c[2]*(YY-y_m[2])**2 )
+  A[3]*np.exp( a[3]*(XX-x_m[3])**2 + b[3]*(XX-x_m[3])*(YY-y_m[3]) + c[3]*(YY-y_m[3])**2 )



fig, ax = plt.subplots()

c=ax.contourf(XX, YY, Z)
plt.colorbar(c)
ax.set_xlabel('x')
ax.set_ylabel('y')

edit: Setting the limits of the plot region to the ones the authors does not seem to help. If I use

x = np.linspace(-2, 0, 1000)
y = np.linspace(0, 2, 1000)

I see this:

user_3944322 · Answer

You forgot the parantheses around the expression for Z so that you effectively evaluated only the first summand.

x = np.linspace(-1.5, 0.2, 1000)
y = np.linspace(0, 1.9, 1000)
XX, YY = np.meshgrid(x, y)

Z = (A[0]*np.exp( a[0]*(XX-x_m[0])**2 + b[0]*(XX-x_m[0])*(YY-y_m[0]) + c[0]*(YY-y_m[0])**2 )
    +A[1]*np.exp( a[1]*(XX-x_m[1])**2 + b[1]*(XX-x_m[1])*(YY-y_m[1]) + c[1]*(YY-y_m[1])**2 )
    +A[2]*np.exp( a[2]*(XX-x_m[2])**2 + b[2]*(XX-x_m[2])*(YY-y_m[2]) + c[2]*(YY-y_m[2])**2 )
    +A[3]*np.exp( a[3]*(XX-x_m[3])**2 + b[3]*(XX-x_m[3])*(YY-y_m[3]) + c[3]*(YY-y_m[3])**2 ))

Bonus:

m1 = (-0.558223634633024, 1.441725841804669)
m2 = (-0.050010822998206, 0.466694104871972)
s1 = (-0.822001558732732, 0.624312802814871)
plt.plot(*m1, 'm*'), plt.text(*m1, "  Min 1")
plt.plot(*m2, 'm*'), plt.text(*m2, "  Min 2")
plt.plot(*s1, 'bo'), plt.text(*s1, "  Saddle 1")

Plot a well known double well surface

1 Answers