CODEWITHSUNDEEP
pythonmatplotlib
user12298516
user12298516

Reputation:

Plotting a harmonic potential as a function of distance or angles

I'm trying to plot the harmonic potential U = 0.5 k(x-x_0)^2. as a function of distance and the angles. I'm taking the distance between the atoms and the angles between them from a .dat file, which I'll provide below, but I'm getting the strange plotting result. Am I missing something in the code??

Here is my code

import numpy as np
import matplotlib.pyplot as plt
#angles
a_i, a_j, a_k, theta_0, k_theta = np.loadtxt('angles.dat', unpack=True)

#bonds
b_i, b_j, r_0, k_b = np.loadtxt('bonds.dat',unpack=True)

def U_angles(theta, theta_0, k_theta):
    U = 0.5*k_theta*(theta-theta_0)**2
    return U

def U_bonds(r, r_0, k_b):
    U = 0.5*k_b*(r-r_0)**2
    return U
r = b_j - b_i
plt.plot(r, U_bonds(r, r_0, k_b))
plt.show()

Here are the data file, which I'm using Data files

Upvotes: 0

Views: 138

Answers (1)

Mr. T
Mr. T

Reputation: 12410

Your array r contains some values multiple times and is not sorted. Matplotlib plots the values as you provide them. You have to sort your array and find the corresponding U values. Luckily, np.unique() does all this for us - returns unique values, sorts them, and returns their indexes, if you ask nicely:

import numpy as np
import matplotlib.pyplot as plt
#angles
a_i, a_j, a_k, theta_0, k_theta = np.loadtxt('angles.dat', unpack=True)

#bonds
b_i, b_j, r_0, k_b = np.loadtxt('bonds.dat',unpack=True)

def U_angles(theta, theta_0, k_theta):
    U = 0.5*k_theta*(theta-theta_0)**2
    return U

def U_bonds(r, r_0, k_b):
    U = 0.5*k_b*(r-r_0)**2
    return U

r = b_j - b_i
Uplot = U_bonds(r, r_0, k_b)

#find unique, sorted r-values and retrieve their indexes
rplot, rind = np.unique(r, return_index=True)

plt.plot(rplot, Uplot[rind])
plt.show()

Sample output: enter image description here

However, this is a patchwork solution. If U_bonds were a function containing a complicated calculation, you would fare better to only calculate the U-values for unique r-values.

...
r = b_j - b_i

#find unique, sorted r-values and retrieve their indexes
_, rind = np.unique(r, return_index=True)

plt.plot(r[rind], U_bonds(r[rind], r_0[rind], k_b[rind]))
plt.show()

Upvotes: 1

Related Questions