Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
import matscipy
import numpy as np
import ase
from ase import Atoms
from ase.calculators.lj import LennardJones
# set up a heptamer
r0 = 1.0
N = 6
phi = np.linspace(0., (N-1)* 2*np.pi/N, N)
x0 = 2.5 # position of center
y0 = 2.5
x = x0 + np.r_[0, r0 * np.cos(phi)]
y = y0 + np.r_[0, r0 * np.sin(phi)]
X = np.c_[x, y, np.zeros(N+1)]
at = ase.Atoms("H%d" % len(X),
positions=X,
cell=np.diag([5, 5, 1.0]),
pbc=[True, True, True])
lj = LennardJones(epsilon=1., sigma=1.)
at.set_calculator(lj);
kBT = 0.1
nsteps = 10**5
dt = 1.e-4
E0 = at.get_potential_energy(at)
print("Initial Energy = ", E0)
Evals = [0.]
#integrate and compute energies
for j in range(nsteps):
x = at.get_positions().reshape(-1)
f = at.get_forces().reshape(-1)
x += f*dt + np.sqrt(2. * kBT *dt) * np.random.normal(size=len(x))
at.set_positions(x.reshape((len(at), 3)))
Evals.append(at.get_potential_energy(at)-E0)
if j % 1000 == 0:
print(j, np.mean(Evals))
Ebar = np.mean(Evals)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.