Skip to content

Instantly share code, notes, and snippets.

@cortner
Created December 4, 2017 13:31
Show Gist options
  • Save cortner/9e2bb6f5cec9994b5d6d1342cc450573 to your computer and use it in GitHub Desktop.
Save cortner/9e2bb6f5cec9994b5d6d1342cc450573 to your computer and use it in GitHub Desktop.
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