From da78c7b0bb9631eec7d26a01daf78187c0ab6217 Mon Sep 17 00:00:00 2001
From: kkly1995 <37456935+kkly1995@users.noreply.github.com>
Date: Wed, 27 Apr 2022 15:25:45 -0500
Subject: [PATCH] Create lj_npt

---
 bin/lj_npt | 88 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 88 insertions(+)
 create mode 100755 bin/lj_npt

diff --git a/bin/lj_npt b/bin/lj_npt
new file mode 100755
index 0000000..d65d9a2
--- /dev/null
+++ b/bin/lj_npt
@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+import numpy as np
+import mcase.particle_moves as mcp
+import mcase.cell_moves as mcc
+from ase.io import read, write
+from ase.units import kB
+from ase.calculators.lj import LennardJones
+from mcase.reader import read_input
+from sys import argv
+
+settings = read_input(argv[1])
+rng = np.random.default_rng()
+
+# setup system
+walker = read(settings['xyz'])
+walker.calc = LennardJones(epsilon=float(settings['epsilon']),\
+        sigma=float(settings['sigma']))
+
+# setup simulation
+temperature = float(settings['temperature'])
+pressure = float(settings['pressure'])
+timestep = float(settings['timestep'])
+cellstep = float(settings['cellstep'])
+number_samples = int(settings['samples'])
+beta = 1.0/(kB*temperature)
+scale = np.sqrt(2*timestep)
+
+# successful read, print out stuff
+print('read inputs from ' + argv[1])
+print('particles read from ' + settings['xyz'])
+print('Lennard-jones potential, with epsilon = %s and sigma = %s'\
+        % (settings['epsilon'], settings['sigma']))
+print('temperature: %s K (beta = %s eV^-1)' % (temperature, beta))
+print('pressure: %s eV/A^3' % pressure)
+print('timestep: %s (length scale = %s)' % (timestep, scale))
+print('cell step: %s A' % cellstep)
+print('number of samples to calculate: %d' % number_samples)
+
+# start sampling
+p_accepted = 0 # for particle moves
+c_accepted = 0 # for cell moves
+r = walker.get_positions()
+energy = walker.get_potential_energy()
+force = walker.get_forces()
+cell = walker.cell.cellpar()[:3] # stick to orthorhombic changes
+volume = walker.get_volume()
+N = len(walker)
+
+for _ in range(number_samples):
+    # particle move
+    dr = rng.normal(size=r.shape, scale=scale) + timestep*beta*force
+    r_new = r + dr
+    walker.set_positions(r_new)
+    energy_new = walker.get_potential_energy()
+    force_new = walker.get_forces()
+    # accept or reject
+    if rng.random() < mcp.acceptance_ratio(r, r_new, energy, energy_new,\
+            force, force_new, beta, timestep):
+        energy = energy_new
+        force = np.copy(force_new)
+        r = walker.get_positions(wrap=True) # wrap to avoid blowing up
+        p_accepted += 1
+    else:
+        walker.set_positions(r)
+    # cell move
+    dc = cellstep*(rng.random(size=3) - 0.5)
+    cell_new = cell + dc
+    walker.set_cell(cell_new, scale_atoms=True)
+    volume_new = walker.get_volume()
+    energy_new = walker.get_potential_energy()
+    if rng.random() < mcc.acceptance_ratio(energy, energy_new,\
+            volume, volume_new, N, pressure, beta):
+        energy = energy_new
+        force = walker.get_forces()
+        r = walker.get_positions(wrap=True)
+        cell = np.copy(cell_new)
+        volume = volume_new
+        c_accepted += 1
+    else:
+        walker.set_cell(cell, scale_atoms=True)
+    print(energy, cell[0], cell[1], cell[2])
+
+# finished
+print('particle move acceptance rate: %.2f' % (p_accepted / number_samples))
+print('cell move acceptance rate: %.2f' % (c_accepted / number_samples))
+
+print('writing final configuration to', settings['save_xyz'])
+write(settings['save_xyz'], walker, format='extxyz')