Gideon-crash.py

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)