LJ Potential Plotter

released on 11 July 2022 Lin of this site : https://sites.google.com/view/riset-arkundato/lj_pot_plotter

LJ Plot using GNUPLOT

#https://stackoverflow.com/questions/42402853/how-can-i-plot-a-lennard-jones-potential-like-this-one-with-gnuplot

f(r,sigma,epsilon)=4*epsilon*((sigma/r)**12-(sigma/r)**6)

sigma=2.6

epsilon=5

set arrow 1 from first 3,0 to first 3, f(3,sigma,epsilon)

set arrow 2 from first 3,0 to first 2.8, f(2.8,sigma,epsilon)

set label 1 "(r_0, {/Symbol e}_0)" at first 3,0 center offset 0,1

set xlabel "Distance r_0 ({\305})"

set ylabel "Energy E (J/mol)

set xr [2:5.2]

set yr [-6:12]

set key reverse Left at 4,5

plot for [epsi= epsilon-1:epsilon+1:1 ] f(x,sigma,epsi) title sprintf("12-6 LJ {/Symbol s}=%.1f {/Symbol e}=%.1f",sigma,epsi)

set term pngcairo enhanced

set out "LJ.png"

replot

#https://scipython.com/book2/chapter-3-simple-plotting/problems/p33/the-lennard-jones-interaction-potential/

import numpy as np

import matplotlib.pyplot as plt

# The Lennard-Jones parameters:

A = 1.024e-23 # J.nm^6

B = 1.582e-26 # J.nm^12

# Adjust the units of A and B - they have more manageable values

# in K.nm^6 and K.nm^12

kB = 1.381e-23 # Boltzmann's constant, J/K

A, B = A / kB, B / kB

# Interatomic distance, in nm

r = np.linspace(0.3, 1., 1000)

# Interatomic potential

U = B/r**12 - A/r**6

plt.plot(r, U, 'k', lw=2., label='Lennard-Jones')

r0 = (2*B/A)**(1./6)

epsilon = B/r0**12 - A/r0**6

k = 156*B/r0**14 - 42*A/r0**8

V = 0.5 * k * (r-r0)**2 + epsilon

plt.plot(r, V, 'k', ls=':', lw=2., label='Harmonic approximation')

plt.xlim(0.3, 0.6)

plt.ylim(-120, 0)

plt.xlabel('r / nm')

plt.ylabel('Potential energy / K')

plt.legend(loc=4)

plt.show()

Page updated

Google Sites

Report abuse