Merge pull request #27 from libAtoms/cell_step_param_scaling

Cell step param scaling

Author: bernstei

README_input_parameters.md

@@ -50,14 +50,14 @@
     - cell_proportion = 0.0
     - type_proportion = 0.0
     - [max_step_size]
-        - pos_gmc_each_atom = -1.0
-        - cell_volume_per_atom = -1.0
-        - cell_shear = 0.2
+        - pos_gmc_each_atom = -0.1
+        - cell_volume_per_atom = -0.05
+        - cell_shear_per_rt3_atom = -1.0
         - cell_stretch = 0.2
     - [step_size]
         - pos_gmc_each_atom = -1.0
         - cell_volume_per_atom = -1.0
-        - cell_shear = -1.0
+        - cell_shear_per_rt3_atom = -1.0
         - cell_stretch = -1.0
     - [cell]
         - min_aspect_ratio = 0.8

README_usage.md

@@ -111,10 +111,21 @@ Hopefully section and key names are self explanatory.
 
 ### Additional notes on run parameters
 
-  - all intervals are in NS iterations, except `stdout_report_interval_s` which is in seconds
+  - All intervals are in NS iterations, except `stdout_report_interval_s` which is in seconds
   - `configs.calculator.type` can be `"ASE"` or `"LAMMPS"`
     - if `"LAMMPS"`, `args` consists of `cmds`, with LAMMPS commands, and `types` dict (one key for each species)
     - if `"ASE"`, `args` consists of `module` key with module that defines a `calc` symbol containing an `ase.calculators.Calculator` object
   - `configs.walk.*_traj_len` controls the number of steps in a walk block of that type
   - `configs.walk.*_proportion` controls the fraction of steps overall that are used for that type of move
-  - if `config.walk.type.sGC = true`, a dict of `mu` values, one per species, is required.
+  - If `config.walk.type.sGC = true`, a dict of `mu` values, one per species, is required.
+  - There are position, cell, and atom-type walks, with associated step size parameters that are auto-tuned. Default values depends on
+    an overall volume scale given by `initial_rand_vol_per_atom` and corresponding length scale given by its cube root.
+    - Maxima are in `[configs.walk.max_step_size]` section. Defaults for first three are negative.
+      - `pos_gmc_each_atom`: distance (typically A) that each atom should typically make in GMC step. If negative, used as multiplier for
+        length scale.
+      - `cell_volume_per_atom`: change in volume (typically A^3), will also be scaled by number of atoms. If negative, used as multiplier
+        for volume scale.
+      - `cell_shear_per_rt3_atom`: cell shear magnitude (typically A) which multiplies _normalized_ cell vectors, will also be scaled by 
+        cube root of number of atoms.  If negative used as multiplier for length scale.
+      - `cell_stretch`: cell stretch, fractional (i.e. strain)
+    - Initial values with same key names are in `[configs.walk.step_size]` section. Any negative values are replaced with half the corresponding maximum.

pymatnext/ns_configs/ase_atoms/__init__.py

@@ -62,7 +62,7 @@ class NSConfig_ASE_Atoms():
     filename_suffix = ".extxyz"
     n_quantities = -1
 
-    _step_size_params = ["pos_gmc_each_atom", "cell_volume_per_atom", "cell_shear", "cell_stretch"]
+    _step_size_params = ["pos_gmc_each_atom", "cell_volume_per_atom", "cell_shear_per_rt3_atom", "cell_stretch"]
     _max_E_hist = collections.deque(maxlen=1000)
     _walk_moves = ["gmc", "cell", "type"]
     _Zs = []
@@ -432,9 +432,11 @@ def _prep_walk(self, params, vol_per_atom=None):
         assert set(list(self.max_step_size.keys())) == set(self._step_size_params)
         # max step size for position GMC and cell volume defaults are scaled to volume per atom
         if self.max_step_size["pos_gmc_each_atom"] < 0.0:
-            self.max_step_size["pos_gmc_each_atom"] = (vol_per_atom ** (1.0/3.0)) / 10.0
+            self.max_step_size["pos_gmc_each_atom"] = (vol_per_atom ** (1.0/3.0)) * np.abs(self.max_step_size["pos_gmc_each_atom"])
         if self.max_step_size["cell_volume_per_atom"] < 0.0:
-            self.max_step_size["cell_volume_per_atom"] = vol_per_atom / 20.0
+            self.max_step_size["cell_volume_per_atom"] = vol_per_atom * np.abs(self.max_step_size["cell_volume_per_atom"])
+        if self.max_step_size["cell_shear_per_rt3_atom"] < 0.0:
+            self.max_step_size["cell_shear_per_rt3_atom"] = (vol_per_atom ** (1.0/3.0)) * np.abs(self.max_step_size["cell_shear_per_rt3_atom"])
 
         # actual step sizes
         self.step_size = params["step_size"].copy()

pymatnext/ns_configs/ase_atoms/ase_atoms_params.py

@@ -27,16 +27,16 @@
     "type_proportion": 0.0,
 
     "max_step_size": {
-        "pos_gmc_each_atom": -1.0,
-        "cell_volume_per_atom": -1.0,
-        "cell_shear": 0.2,
+        "pos_gmc_each_atom": -0.1,
+        "cell_volume_per_atom": -0.05,
+        "cell_shear_per_rt3_atom": -1.0,
         "cell_stretch": 0.2
     },
 
     "step_size": {
         "pos_gmc_each_atom": -1.0,
         "cell_volume_per_atom": -1.0,
-        "cell_shear": -1.0,
+        "cell_shear_per_rt3_atom": -1.0,
         "cell_stretch": -1.0
     },
 

pymatnext/ns_configs/ase_atoms/walks_ase_calc.py

@@ -126,12 +126,15 @@ def walk_cell(ns_atoms, Emax, rng):
 
     Returns
     -------
-    [("cell_volume_per_atom", int n_attempt, int n_success), ("cell_shear", ...), ("cell_stretch", ...)] info on move params and attempts/successes
+    [("cell_volume_per_atom", int n_attempt, int n_success),
+     ("cell_shear_per_rt3_atom", int n_attempt, n_success),
+     ("cell_stretch", int n_attempt, n_success)] with n_att and n_acc number of attempted and accepted
+     move for each submove type
     """
     atoms = ns_atoms.atoms
     N_atoms = len(atoms)
-    step_size_volume = N_atoms * ns_atoms.step_size["cell_volume_per_atom"]
-    step_size_shear = ns_atoms.step_size["cell_shear"]
+    step_size_volume = ns_atoms.step_size["cell_volume_per_atom"] * N_atoms
+    step_size_shear = ns_atoms.step_size["cell_shear_per_rt3_atom"] * (N_atoms ** (1.0 / 3.0))
     step_size_stretch = ns_atoms.step_size["cell_stretch"]
     min_aspect_ratio = ns_atoms.move_params["cell"]["min_aspect_ratio"]
     flat_V_prior = ns_atoms.move_params["cell"]["flat_V_prior"]
@@ -197,7 +200,7 @@ def walk_cell(ns_atoms, Emax, rng):
             n_acc[move] += 1
 
     return [("cell_volume_per_atom", n_att["volume"], n_acc["volume"]),
-            ("cell_shear", n_att["shear"], n_acc["shear"]),
+            ("cell_shear_per_rt3_atom", n_att["shear"], n_acc["shear"]),
             ("cell_stretch", n_att["stretch"], n_acc["stretch"])]
 
 

pymatnext/ns_configs/ase_atoms/walks_lammps.py

@@ -187,18 +187,19 @@ def walk_cell(ns_atoms, Emax, rng):
 
     Returns
     -------
-    [("cell_volume_per_atom", n_att, n_acc),
-     ("cell_shear", n_att, n_acc),
-     ("cell_stretch", n_att, n_acc)] with n_att and n_acc number of attempted and accepted
+    [("cell_volume_per_atom", int n_attempt, int n_success),
+     ("cell_shear_per_rt3_atom", int n_attempt, n_success),
+     ("cell_stretch", int n_attempt, n_success)] with n_att and n_acc number of attempted and accepted
      move for each submove type
     """
     ns_atoms.calc.command("unfix NS")
     atoms = ns_atoms.atoms
     # for LAMMPS RanMars RNG
     lammps_seed = rng.integers(1, 900000000)
 
-    step_size_volume = len(atoms) * ns_atoms.step_size["cell_volume_per_atom"]
-    step_size_shear = ns_atoms.step_size["cell_shear"]
+    N_atoms = len(atoms)
+    step_size_volume = ns_atoms.step_size["cell_volume_per_atom"] * N_atoms
+    step_size_shear = ns_atoms.step_size["cell_shear_per_rt3_atom"] * (N_atoms ** (1.0 / 3.0))
     step_size_stretch = ns_atoms.step_size["cell_stretch"]
 
     types, pos, vel = set_lammps_from_atoms(ns_atoms)
@@ -256,7 +257,7 @@ def walk_cell(ns_atoms, Emax, rng):
                 n_acc[submove_type] = 0
 
     return [("cell_volume_per_atom", n_att["volume"], n_acc["volume"]),
-            ("cell_shear", n_att["shear"], n_acc["shear"]),
+            ("cell_shear_per_rt3_atom", n_att["shear"], n_acc["shear"]),
             ("cell_stretch", n_att["stretch"], n_acc["stretch"])]
 
 def walk_type(ns_atoms, Emax, rng):

tests/test_sample.py

@@ -158,7 +158,7 @@ def do_Morse_ASE(tmp_path, monkeypatch, using_mpi, max_iter=None):
     assert len(list(tmp_path.glob('Morse_ASE.test.traj.*xyz'))) == 1
     assert len(list(ase.io.read(tmp_path / 'Morse_ASE.test.traj.extxyz', ':'))) == int(np.ceil(max_iter_use / traj_interval))
 
-    # from test run 6/13/2025, when max iter was extended to 110 to catch deadlock in old buggy snapshot step_size writing
+    # from test run 6/20/2025, when cell_shear_per_rt3_atom default was increased
     if using_mpi:
         samples_fields_ref = np.asarray([1.09000000e+02, 8.02719236e+00, 1.28609799e+04, 1.60000000e+01])
     else:
@@ -172,9 +172,7 @@ def do_Morse_ASE(tmp_path, monkeypatch, using_mpi, max_iter=None):
     # tolerance loosened so that restart, which isn't perfect due to finite precision in
     # saved cofig file, still passes
     if not np.allclose(samples_fields, samples_fields_ref, rtol=0.02):
-        print("final samples line test", samples_fields)
-        print("final samples line ref ", samples_fields_ref)
-        assert False
+        assert False, f"test {samples_fields} ref {samples_fields_ref}"
 
     # this will fail if number of steps is not divisible by N_samples interval, because
     # clone_hist always saves every line. Also, clone_hist is a hack that's truncated by
@@ -243,21 +241,19 @@ def do_EAM_LAMMPS(tmp_path, monkeypatch, using_mpi, max_iter=None):
     assert len(list(tmp_path.glob('EAM_LAMMPS.test.traj.*xyz'))) == 1
     assert len(list(ase.io.read(tmp_path / 'EAM_LAMMPS.test.traj.extxyz', ':'))) == max_iter_use // traj_interval
 
-    # from test run 12/8/2022
+    # from test run 6/20/2025 when cell shear default was increased
     if using_mpi:
-        fields_ref = np.asarray([2.99000000e+02, -3.91163426e+02,  1.08674253e+04,  1.60000000e+01, 1.87500000e-01,  8.12500000e-01])
+        fields_ref = np.asarray([2.99000000e+02, -3.90328931e+02,  1.99365520e+04,  1.60000000e+01, 1.87500000e-01,  8.12500000e-01])
     else:
-        fields_ref = np.asarray([299, -366.1823642208, 6004.3693892916, 16.0000000000, 0.2500000000, 0.7500000000])
+        fields_ref = np.asarray([2.99000000e+02, -3.90994404e+02,  1.45640509e+04,  1.60000000e+01, 1.87500000e-01,  8.12500000e-01])
 
     with open(tmp_path / 'EAM_LAMMPS.test.NS_samples') as fin:
         for l in fin:
             pass
     fields = np.asarray([float(f) for f in l.strip().split()])
 
     if not np.allclose(fields, fields_ref):
-        print("final line test", fields)
-        print("final line ref ", fields_ref)
-        assert False
+        assert False, f"test {fields} ref {fields_ref}"
 
 
 def do_pressure(tmp_path, monkeypatch, using_mpi):
@@ -308,7 +304,7 @@ def do_pressure(tmp_path, monkeypatch, using_mpi):
 
         Vfirst = float(lfirst.strip().split()[2])
         Vlast = float(llast.strip().split()[2])
-        assert Vfirst / Vlast > 20
+        assert Vfirst / Vlast > 10
 
 
 def do_sGC(tmp_path, monkeypatch, using_mpi):
@@ -366,4 +362,4 @@ def do_sGC(tmp_path, monkeypatch, using_mpi):
         f_13_last, f_29_last = [float(f) for f in llast.strip().split()[4:6]]
 
         assert f_29_first / f_13_first == 1.0
-        assert f_29_last / f_13_last > 5
+        assert f_29_last / f_13_last > 4