From 43c3da3901d111457ddbca78492bb2c2ea02edec Mon Sep 17 00:00:00 2001 From: VGFletcher <125897513+VGFletcher@users.noreply.github.com> Date: Wed, 4 Sep 2024 13:54:20 +0100 Subject: [PATCH 1/6] Added bond length restriction keyword --- ns_run.py | 137 ++++++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 113 insertions(+), 24 deletions(-) diff --git a/ns_run.py b/ns_run.py index 7863337..c059ae9 100755 --- a/ns_run.py +++ b/ns_run.py @@ -1,3 +1,9 @@ +#This is a modified version of the pymatnest ns_run.py file, with an additional criteria to reject +#configurations if they contain bond lengths less than a given bond length defined by the +#'min_bond_len' keyword + +#Modification made by @V.G.Fletcher + import re, math, time, os import pprint import numpy as np, ase, ase.io @@ -1083,8 +1089,29 @@ def do_MD_atom_walk(at, movement_args, Emax, KEmax): reject_Emax = (final_E >= Emax) reject_KEmax = (KEmax > 0.0 and final_KE >= KEmax) + #################################VGF_MODIFIED######################################### + if ns_args['calc_distances']: + distances = ase.geometry.get_distances(at.get_positions(), pbc=True, cell=at.get_cell()) + bonds = np.sort(distances[1], axis=None) + config_size = len(at) + min_bond = bonds[config_size] + + reject_len = min_bond < ns_args['min_bond_len'] + + num_bonds= bonds[config_size:] < (ns_args['min_bond_len'] + 0.1) + perc_holes = sum(num_bonds)/len(num_bonds) + + reject_pot_hole = perc_holes > ns_args['hole_tolerance'] + + if reject_pot_hole and not (reject_fuzz or reject_Emax or reject_KEmax or reject_len): + ase.io.write(f"potential_holes.traj.{rank}.extxyz", at, parallel=False, format=ns_args['config_file_format'], append=True) + else: + reject_len = False + reject_pot_hole = False + ##################################################################################### + #DOC \item if reject - if reject_fuzz or reject_Emax or reject_KEmax: # reject + if reject_fuzz or reject_Emax or reject_KEmax or reject_len or reject_pot_hole: # reject #DOC \item set positions, velocities, energy back to value before perturbation (maybe should be after?) # print(print_prefix, ": WARNING: reject MD traj Emax ", Emax, " initial E ", orig_E, " velo perturbed E ", pre_MD_E, " final E ",final_E, " KEmax ", KEmax, " KE ", final_KE) at.set_positions(pre_MD_pos) @@ -1444,24 +1471,50 @@ def do_cell_step(at, Emax, p_accept, transform): # set new positions and velocities at.set_cell(new_cell, scale_atoms=True) - + ################################################VGF-MODIFIED####################################### + if ns_args['calc_distances']: + distances = ase.geometry.get_distances(at.get_positions(), pbc=True, cell=at.get_cell()) + bonds = np.sort(distances[1], axis=None) + config_size = len(at) + min_bond = bonds[config_size] + + accept_len = min_bond > ns_args['min_bond_len'] + + num_bonds= bonds[config_size:] < (ns_args['min_bond_len'] + 0.1) + perc_holes = sum(num_bonds)/len(num_bonds) + + accept_pot_hole = perc_holes < ns_args['hole_tolerance'] + else: + accept_len = True + accept_pot_hole = True + ################################################################################################### if Emax is None: return - # calculate new energy - try: - new_energy = eval_energy(at) - except Exception as err: - if ns_args['debug'] >= 4: - print( "eval_energy got exception ", err) - new_energy = 2.0*abs(Emax) + if accept_len: + # calculate new energy + try: + new_energy = eval_energy(at) + except Exception as err: + if ns_args['debug'] >= 4: + print( "eval_energy got exception ", err) + new_energy = 2.0*abs(Emax) #print("error in eval_energy setting new_energy = 2*abs(Emax)=" , new_energy) - # accept or reject - if new_energy < Emax: # accept - at.info['ns_energy'] = new_energy - return True - else: # reject and revert + # accept or reject + if new_energy < Emax and accept_pot_hole: # accept + at.info['ns_energy'] = new_energy + return True + else: # reject and revert + if new_energy < Emax: + ase.io.write(f"potential_holes.traj.{rank}.extxyz", at, parallel=False, format=ns_args['config_file_format'], append=True) + at.set_cell(orig_cell,scale_atoms=False) + at.set_positions(orig_pos) + if ns_args['n_extra_data'] > 0: + at.arrays['ns_extra_data'][...] = extra_data + return False + + else: at.set_cell(orig_cell,scale_atoms=False) at.set_positions(orig_pos) if ns_args['n_extra_data'] > 0: @@ -2088,10 +2141,10 @@ def full_auto_set_stepsizes(walkers, walk_stats, movement_args, comm, Emax, KEma #update step length dir = None if rate < min_rate: - exp = -1.0 + exp = -1 dir = "down" elif rate >= max_rate: - exp = 1.0 + exp = 1 dir = "up" else: exp = None @@ -2176,10 +2229,10 @@ def adjust_step_sizes(walk_stats, movement_args, comm, do_print_rate=True, monit dir = None exp = 0.0 if rate < min_rate: - exp = -1.0 + exp = -1 dir = "down" elif rate >= max_rate: - exp = 1.0 + exp = 1 dir = "up" orig_value = movement_args[key+"_"+suffix] @@ -2551,6 +2604,7 @@ def do_ns_loop(): walker_copy.info['volume'] = walker_copy.get_volume() walker_copy.info['ns_P'] = movement_args['MC_cell_P'] walker_copy.info['iter'] = i_ns_step + walker_copy.info['config_type'] = "config_{}".format(i_ns_step) #VF ADDED for ACEFIT test-train walker_copy.info['config_n_global'] = global_n if walker_copy.has('masses') and walker_copy.has('momenta'): walker_copy.info['ns_KE'] = walker_copy.get_kinetic_energy() @@ -3106,11 +3160,11 @@ def main(): global Z_list import sys - + sys.excepthook = excepthook_mpi_abort - + stacktrace.listen() - + if len(sys.argv) != 1 and len(sys.argv) != 2: usage() sys.exit(1) @@ -3238,6 +3292,12 @@ def main(): ns_args['random_energy_perturbation'] = float(args.pop('random_energy_perturbation', 1.0e-12)) ns_args['n_extra_data'] = int(args.pop('n_extra_data', 0)) ns_args['Z_cell_axis'] = float(args.pop('Z_cell_axis', 10.0)) + #VGF parameters + ns_args['min_bond_len'] = float(args.pop('min_bond_len', 0.0)) + ns_args['hole_tolerance'] = float(args.pop('hole_tolerance',1.0)) + + ns_args['calc_distances'] = (ns_args['min_bond_len'] != 0.0) or (ns_args['hole_tolerance'] != 1.0) + #End VGF parameters # surely there's a cleaner way of doing this? try: @@ -3813,15 +3873,37 @@ def main(): # random initial positions energy = float('nan') n_try = 0 - while n_try < ns_args['random_init_max_n_tries'] and (math.isnan(energy) or energy > ns_args['start_energy_ceiling']): + ###################################VGF-MODIFIED############################### + reject_len = True + reject_pot_hole = True + while (n_try < ns_args['random_init_max_n_tries']) and (((math.isnan(energy) or energy > ns_args['start_energy_ceiling'])) or reject_len or reject_pot_hole): at.set_scaled_positions( rng.float_uniform(0.0, 1.0, (len(at), 3) ) ) if movement_args['2D']: # zero the Z coordiates in a 2D simulation temp_at=at.get_positions() temp_at[:,2]=0.0 at.set_positions(temp_at) + ################################################################## + if ns_args['calc_distances']: + distances = ase.geometry.get_distances(at.get_positions(), pbc=True, cell=at.get_cell()) + bonds = np.sort(distances[1], axis=None) + config_size = len(at) + min_bond = bonds[config_size] + #print('min_bond_here', min_bond) + + reject_len = min_bond < ns_args['min_bond_len'] + + num_bonds= bonds[config_size:] < (ns_args['min_bond_len'] + 0.1) + perc_holes = sum(num_bonds)/len(num_bonds) + + reject_pot_hole = perc_holes > ns_args['hole_tolerance'] + else: + reject_len = False + reject_pot_hole = False + ##################################################################### energy = eval_energy(at) n_try += 1 - if math.isnan(energy) or energy > ns_args['start_energy_ceiling']: + + if math.isnan(energy) or energy > ns_args['start_energy_ceiling'] or reject_len or reject_pot_hole: sys.stderr.write("WARNING: rank %d failed to generate initial config by random positions under max energy %f in %d tries\n" % (rank, ns_args['start_energy_ceiling'], ns_args['random_init_max_n_tries'])) # try FORTRAN config initializer @@ -3847,7 +3929,14 @@ def main(): energy = eval_energy(at) at.info['ns_energy'] = rand_perturb_energy(energy, ns_args['random_energy_perturbation']) at.info['volume'] = at.get_volume() - + ###VGF modified### + if ns_args['calc_distances']: + distances = ase.geometry.get_distances(at.get_positions()) + bonds = np.sort(distances[1], axis=None) + min_bond = bonds[len(at)] + print('min_bond_final', rank, min_bond) + ################### + # Done initialising atomic positions. Now initialise momenta # set KEmax from P and Vmax From 5d3a20a64f97f89b1a772c50a8d6c0310b4f46a0 Mon Sep 17 00:00:00 2001 From: Vincent Fletcher Date: Mon, 9 Sep 2024 17:17:07 +0100 Subject: [PATCH 2/6] edited docs --- ns_run.py | 146 ++++++++++++++++++++++++++++++------------------------ 1 file changed, 80 insertions(+), 66 deletions(-) diff --git a/ns_run.py b/ns_run.py index c059ae9..66a99a0 100755 --- a/ns_run.py +++ b/ns_run.py @@ -902,7 +902,7 @@ def rotate_dir_3N_2D(vec, max_ang): def rej_free_perturb_velo(at, Emax, KEmax, rotate=True): #DOC -#DOC ``rej_free_perturb_velo`` +#DOC ``rej_free_perturb_velo(at, Emax, KEmax, rotate=True)`` if not at.has('momenta') or not at.has('masses'): return @@ -920,22 +920,25 @@ def rej_free_perturb_velo(at, Emax, KEmax, rotate=True): exit_error("rej_free_perturb_velo() called with Emax and KEmax both None\n", 9) orig_KE = at.get_kinetic_energy() - #DOC \item if atom\_velo\_rej\_free\_fully\_randomize, pick random velocities consistent with Emax + #DOC **if** atom\_velo\_rej\_free\_fully\_randomize: + #DOC \item pick random velocities consistent with Emax if movement_args['atom_velo_rej_free_fully_randomize']: # randomize completely, fixed atoms are taken care of at.set_velocities(gen_random_velo(at, KEmax_use)) - #DOC \item else perturb velocities + #DOC **else**: (perturb velocities) else: # perturb velo = at.get_velocities() velo_mag = np.linalg.norm(velo) - #DOC \item if current velocity=0, can't rescale, so pick random velocities consistent with Emax + #DOC **if** current velocity=0: + #DOC \item can't rescale, so pick random velocities consistent with Emax if velo_mag == 0.0: # generate random velocities, fixed atoms are taken care of at.set_velocities(gen_random_velo(at, KEmax_use)) - #DOC \item else, pick new random magnitude consistent with Emax, random rotation of current direction with angle uniform in +/- atom\_velo\_rej\_free\_perturb\_angle + #DOC **else**: + #DOC \item pick new random magnitude consistent with Emax, random rotation of current direction with angle uniform in +/- atom\_velo\_rej\_free\_perturb\_angle else: # using default settings we will do this. # pick new random magnitude - count on dimensionality to make change small @@ -967,14 +970,15 @@ def rej_free_perturb_velo(at, Emax, KEmax, rotate=True): def do_MC_atom_velo_walk(at, movement_args, Emax, nD, KEmax): #DOC -#DOC ``do\_MC\_atom\_velo\_walk`` +#DOC ``do_MC_atom_velo_walk`` - #DOC \item Else if MC\_atom\_velo\_walk\_rej\_free + #DOC **if** MC\_atom\_velo\_walk\_rej\_free: if movement_args['MC_atom_velo_walk_rej_free']: #DOC \item call rej\_free\_perturb\_velo() rej_free_perturb_velo(at, Emax, KEmax) return {} - #DOC \item else do MC pertubation to velocities + #DOC **else**: + #DOC \item do MC pertubation to velocities n_steps = movement_args['velo_traj_len'] step_size = movement_args['MC_atom_velo_step_size'] @@ -1025,8 +1029,9 @@ def do_MD_atom_walk(at, movement_args, Emax, KEmax): if orig_E >= Emax: print(print_prefix, ": WARNING: orig_E =", orig_E, " >= Emax =", Emax) - #DOC \item if MD\_atom\_velo\_pre\_perturb, call do\_MC\_atom\_velo\_walk() for magnitude and rotation + #DOC **if** MD\_atom\_velo\_pre\_perturb: if movement_args['MD_atom_velo_pre_perturb']: + #DOC \item call do_MC_atom_velo_walk() for magnitude and rotation do_MC_atom_velo_walk(at, movement_args, Emax, nD, KEmax) pre_MD_pos = at.get_positions() @@ -1035,8 +1040,8 @@ def do_MD_atom_walk(at, movement_args, Emax, KEmax): pre_MD_extra_data = at.arrays['ns_extra_data'].copy() pre_MD_E = at.info['ns_energy'] - - #DOC \item propagate in time atom\_traj\_len time steps of length MD\_atom\_timestep + #DOC propagate in time atom_traj_len time steps of length MD_atom_timestep + #DOC if movement_args['python_MD']: forces = eval_forces(at) final_E = None @@ -1078,14 +1083,16 @@ def do_MD_atom_walk(at, movement_args, Emax, KEmax): reject_fuzz = False final_KE = eval_energy_KE(at) - #DOC \item If MD\_atom\_reject\_energy\_violation is set, accept/reject entire move on E deviating by less than MD\_atom\_energy\_fuzz times kinetic energy + #DOC **if** MD_atom_reject_energy_violation is set: + #DOC \item accept/reject entire move on E deviating by less than MD\_atom\_energy\_fuzz times kinetic energy if abs(final_E-pre_MD_E) > movement_args['MD_atom_energy_fuzz'] * final_KE: if movement_args['MD_atom_reject_energy_violation']: reject_fuzz = True # else: # print(print_prefix, ": WARNING: MD energy deviation > fuzz*final_KE. Pre-MD, post-MD, difference, final_KE ", pre_MD_E, final_E, final_E-pre_MD_E, final_KE) - #DOC \item accept/reject entire move on E < Emax and KE < KEmax + #DOC accept/reject entire move on E < Emax and KE < KEmax + #DOC reject_Emax = (final_E >= Emax) reject_KEmax = (KEmax > 0.0 and final_KE >= KEmax) @@ -1123,9 +1130,10 @@ def do_MD_atom_walk(at, movement_args, Emax, KEmax): at.arrays['ns_extra_data'][...] = pre_MD_extra_data at.info['ns_energy'] = pre_MD_E n_accept = 0 - #DOC \item else + #DOC **else**: (accepted) else: # accept - #DOC \item flip velocities if MD\_atom\_velo\_flip\_accept + #DOC **if** MD\_atom\_velo\_flip\_accept: + #DOC \item flip velocities # remember to reverse velocities on acceptance to preserve detailed # balance, since velocity is (sometimes) being perturbed, not completely # randomized @@ -1135,7 +1143,8 @@ def do_MD_atom_walk(at, movement_args, Emax, KEmax): at.info['ns_energy'] = final_E n_accept = 1 - #DOC \item if MD\_atom\_velo\_post\_perturb, call do\_MC\_atom\_velo\_walk() for magnitude and rotation + #DOC **if** MD\_atom\_velo\_post\_perturb: + #DOC \item call do\_MC\_atom\_velo\_walk() for magnitude and rotation if movement_args['MD_atom_velo_post_perturb']: do_MC_atom_velo_walk(at, movement_args, Emax, nD, KEmax) @@ -1157,7 +1166,8 @@ def do_MC_atom_walk(at, movement_args, Emax, KEmax): if movement_args['2D']: nD=2 - #DOC \item if MC\_atom\_velocities and MC\_atom\_velocities\_pre\_perturb, call do\_MC\_atom\_velo\_walk() to perturb velocities, magnitude and and rotation + #DOC **if** MC\_atom\_velocities **and** MC\_atom\_velocities\_pre\_perturb: + #DOC \item call do\_MC\_atom\_velo\_walk() to perturb velocities, magnitude and rotation if movement_args['MC_atom_velocities'] and movement_args['MC_atom_velocities_pre_perturb']: do_MC_atom_velo_walk(at, movement_args, Emax, nD, KEmax) @@ -1177,7 +1187,7 @@ def do_MC_atom_walk(at, movement_args, Emax, KEmax): else: rotate_dir_3N(at.arrays['GMC_direction'], movement_args['GMC_dir_perturb_angle']) - #DOC \item if using fortran calculator and not reproducible + #DOC **if** using fortran calculator and not reproducible: if do_calc_fortran and not ns_args['reproducible']: #DOC \item call fortran MC code f\_MC\_MD.MC\_atom\_walk at.wrap() @@ -1218,13 +1228,12 @@ def do_MC_atom_walk(at, movement_args, Emax, KEmax): at.set_positions(orig_pos) n_accept = 0 - #DOC \item else + #DOC **else**: (do python MC) else: - #DOC \item do python MC if movement_args['MC_atom_velocities']: exit_error("MC_atom_velocities only supported for FORTRAN calculator\n", 8) dz=0.0 - #DOC \item if MC_atom_Galilean + #DOC **if** MC_atom_Galilean: if movement_args['MC_atom_Galilean']: #DOC \item go Galilean MC in python @@ -1297,11 +1306,11 @@ def do_MC_atom_walk(at, movement_args, Emax, KEmax): n_try = n_reflect + n_reverse n_accept = n_reflect - #DOC \item else + #DOC **else**: else: - #DOC \item loop atom\_traj\_len times + #DOC \item **loop** atom\_traj\_len times for i_MC_step in range(n_steps): - #DOC \item loop over atoms in random order + #DOC \item **loop** over atoms in random order at_list=list(range(len(at))) rng.shuffle_in_place(at_list) for i_at in at_list: @@ -1564,20 +1573,22 @@ def do_MC_swap_step(at, movement_args, Emax, KEmax): #DOC ``do_MC_swap_step`` Z = at.get_atomic_numbers() - #DOC \item return if all atomic numbers are identical + #DOC **if** all atomic numbers are identical: + #DOC **return** if (Z[:] == Z[0]).all(): # don't try to swap when all atoms are the same return (0, {}) r_cut = movement_args['swap_r_cut'] - #DOC \item randomly pick a desired cluster size + #DOC randomly pick a desired cluster size + #DOC cluster_size = np.where(rng.float_uniform(0,1) < movement_args['swap_probs'])[0][0]+1 if cluster_size > 1: (i_list, j_list) = matscipy.neighbours.neighbour_list('ij', at, r_cut) else: i_list = None j_list = None - #DOC \item pick two clusters with distinct atomic numbers, backing off to smaller clusters on failure to find appropriate larger clusters, but always pick at least a pair of atoms to swap + #DOC pick two clusters with distinct atomic numbers, backing off to smaller clusters on failure to find appropriate larger clusters, but always pick at least a pair of atoms to swap c1 = None c2 = None while (c1 is None or c2 is None or np.all(Z[c1] == Z[c2])): @@ -1625,8 +1636,9 @@ def do_MC_swap_step(at, movement_args, Emax, KEmax): extra_data_2_orig = at.arrays['ns_extra_data'][c2,...].copy() at.arrays['ns_extra_data'][c1,...] = extra_data_2_orig at.arrays['ns_extra_data'][c2,...] = extra_data_1_orig - - #DOC \item accept swap if energy < Emax + #DOC + #DOC **if** energy < Emax: + #DOC \item accept swap new_energy = eval_energy(at) new_KE = eval_energy_KE(at) @@ -1694,7 +1706,7 @@ def do_atom_walk(at, movement_args, Emax, KEmax): #DOC ``do_atom_walk`` n_reps = movement_args['n_atom_steps_per_call'] out = {} - #DOC \item loop n\_atom\_steps\_per\_call times, calling do\_MC\_atom\_walk() or do\_MD\_atom\_walk() + #DOC **loop** n\_atom\_steps\_per\_call times, calling do\_MC\_atom\_walk() or do\_MD\_atom\_walk() for i in range(n_reps): if movement_args['atom_algorithm'] == 'MC': accumulate_stats(out, do_MC_atom_walk(at, movement_args, Emax, KEmax)) @@ -1801,33 +1813,32 @@ def walk_single_walker(at, movement_args, Emax, KEmax): accumulate_stats(out, t_out) else: - #DOC \item create list - #DOC \item do\_atom\_walk :math:`*` n\_atom\_step\_n\_calls possible_moves = ([do_atom_walk] * movement_args['n_atom_steps_n_calls'] + - #DOC \item do\_cell\_volume\_step :math:`*` n\_cell\_volume\_steps [do_MC_cell_volume_step] * movement_args['n_cell_volume_steps'] + - #DOC \item do\_cell\_shear\_step :math:`*` n\_cell\_shear\_steps [do_MC_cell_shear_step] * movement_args['n_cell_shear_steps'] + - #DOC \item do\_cell\_stretch\_step :math:`*` n\_cell\_stretch\_steps [do_MC_cell_stretch_step] * movement_args['n_cell_stretch_steps'] + - #DOC \item do\_swap\_step :math:`*` n\_swap\_steps [do_MC_swap_step] * movement_args['n_swap_steps'] + - #DOC \item do\_semi\_grand\_step :math:`*` n\_semi\_grand\_steps [do_MC_semi_grand_step] * movement_args['n_semi_grand_steps']) out = {} n_model_calls_used = 0 - #DOC \item loop while n\_model\_calls\_used < n\_model\_calls + #DOC **loop** **while** n_model_calls_used < n_model_calls: while n_model_calls_used < movement_args['n_model_calls']: - #DOC \item pick random item from list + #DOC \item pick random item from list: + #DOC \item do\_atom\_walk :math:`*` n\_atom\_step\_n\_calls + #DOC \item do\_cell\_volume\_step :math:`*` n\_cell\_volume\_steps + #DOC \item do\_cell\_shear\_step :math:`*` n\_cell\_shear\_steps + #DOC \item do\_cell\_stretch\_step :math:`*` n\_cell\_stretch\_steps + #DOC \item do\_swap\_step :math:`*` n\_swap\_steps + #DOC \item do\_semi\_grand\_step :math:`*` n\_semi\_grand\_steps move = possible_moves[rng.int_uniform(0, len(possible_moves))] - #DOC \item do move + #DOC \item do move (t_n_model_calls, t_out) = move(at, movement_args, Emax, KEmax) n_model_calls_used += t_n_model_calls accumulate_stats(out, t_out) - #DOC \item perturb final energy by random\_energy\_perturbation + #DOC perturb final energy by random\_energy\_perturbation # perturb final energy at.info['ns_energy'] = rand_perturb_energy( at.info['ns_energy'], ns_args['random_energy_perturbation'], Emax) @@ -1876,8 +1887,8 @@ def full_auto_set_stepsizes(walkers, walk_stats, movement_args, comm, Emax, KEma """Automatically set all step sizes. Returns the time (in seconds) taken for the routine to run.""" #DOC #DOC ``full_auto_set_stepsizes`` - #DOC \item Step sizes for each (H)MC move are set via a loop which performs additional exploration moves, calibrating each step size to obtain an acceptance rate inside a specified range. - + #DOC Step sizes for each (H)MC move are set via a loop which performs additional exploration moves, calibrating each step size to obtain an acceptance rate inside a specified range. + #DOC global print_prefix orig_prefix = print_prefix @@ -1893,8 +1904,8 @@ def full_auto_set_stepsizes(walkers, walk_stats, movement_args, comm, Emax, KEma # in each trial we will evolve walk_n_walkers configurations else: walk_n_walkers = n_samples_per_move_type - #DOC \item The routine is MPI parallelised, so that the wall time goes as 1/num\_of\_processes - + #DOC The routine is MPI parallelised, so that the wall time goes as 1/num_of_processes + #DOC key_list=[] for key, value in walk_stats.items(): key_list.append(key) @@ -1945,10 +1956,11 @@ def full_auto_set_stepsizes(walkers, walk_stats, movement_args, comm, Emax, KEma if (i==6): key_list.append("MC_swap_") - #DOC \item For each (H)MC move type the following is performed + #DOC For each (H)MC move type the following is performed: + #DOC for key in key_list: - #DOC \item Set ''movement\_args'' parameters so that only one (H)MC call is made at a time + #DOC \item Set ``movement_args`` parameters so that only one (H)MC call is made at a time # reprogram n_atom_steps_n_calls, n_cell_volume_steps, n_cell_shear_steps, n_cell_stretch_steps, n_swap_steps according to key # possible key values: # MD_atom # atoms HMC @@ -2032,7 +2044,7 @@ def full_auto_set_stepsizes(walkers, walk_stats, movement_args, comm, Emax, KEma else: exit_error("full_auto_set_stepsizes got key '%s', unkown to this routine\n" % key, 5) - #DOC \item Min and max acceptance rates are copied from parameters MC\_adjust\_min\_rate / MD\_adjust\_min\_rate and MC\_adjust\_max\_rate / MD\_adjust\_max\_rate + #DOC \item Min and max acceptance rates are copied from parameters MC_adjust_min_rate / MD_adjust_min_rate and MC_adjust_max_rate / MD_adjust_max_rate if key.find("MC") == 0: if key.find("MC_atom_step_size") and movement_args['MC_atom_Galilean']: @@ -2056,27 +2068,26 @@ def full_auto_set_stepsizes(walkers, walk_stats, movement_args, comm, Emax, KEma # protects against possible future bugs that would be hard to detect - #DOC \item Step size calibration loop: + #DOC \item **loop** for step size calibration, repeat the following 200/num\_of\_MPI\_processes times: dir = None while True: stats = {} # clear acceptance / trial stats for new step size stats_cumul = {} # clear acceptance / trial stats for new step size - #DOC \item Repeat the following 200/num\_of\_MPI\_processes times: - #DOC \item Copy a configuration from the live set (each MPI process chooses a different configuration) + #DOC \item Copy a configuration from the live set (each MPI process chooses a different configuration) for i in range(first_walker,first_walker + walk_n_walkers): k = i%len(walkers) # cycle through walkers array buf = walkers[k].copy() # copy config k into buffer "buf" for walking (walkers array unchanged) buf.calc = walkers[k].calc # copy calculator - #DOC \item Each MPI processes performs one (H)MC move on its cloned configuration + #DOC \item Each MPI processes performs one (H)MC move on its cloned configuration # build up stats from walkers if (not key=="MC_atom_velo"): stats = walk_single_walker(buf, exploration_movement_args, Emax, KEmax) else: stats = do_MC_atom_velo_walk(buf, exploration_movement_args, Emax, nD, KEmax) - #DOC running statistics for the number of accepted/rejected moves on each process are recorded + #DOC \item running statistics for the number of accepted/rejected moves on each process are recorded accumulate_stats(stats_cumul, stats) first_walker = first_walker + walk_n_walkers # cycle through local samples @@ -2094,23 +2105,26 @@ def full_auto_set_stepsizes(walkers, walk_stats, movement_args, comm, Emax, KEma n_accept = n_accept_g[0] rate = float(n_accept)/float(n_try) - #DOC \item The total number of accepted/rejected moves for this step size (summed across all MPI processes) are estabilshed + #DOC \item The total number of accepted/rejected moves for this step size (summed across all MPI processes) are estabilshed if ((comm is None or comm.rank == 0) and (ns_args['debug'] >= 1)): print( print_prefix, "trial stepsize and accept rate for %s = %e , %f (%d)" % (key, movement_args[key+"_"+suffix], rate, n_try)) if (rate>min_rate and ratemax_rate)): - #DOC \item If this is NOT the first time the loop has been performed for this (H)MC move AND we previously obtained an acceptance rate on one side of the desired range, and now find an acceptance rate on the other side of the desired range - #DOC \item Return the step size that gave an acceptance rate closest to the middle of the desired range. + #DOC \item **if** this is **not** the first time the loop has been performed for this (H)MC move **and** we previously obtained an acceptance rate on one side of the desired range: + #DOC \item now find an acceptance rate on the other side of the desired range + #DOC \item **return** the step size that gave an acceptance rate closest to the middle of the desired range. # check which gives a accept_rate closer to centre of acceptable window # and take that @@ -2128,7 +2142,7 @@ def full_auto_set_stepsizes(walkers, walk_stats, movement_args, comm, Emax, KEma if (comm is None or comm.rank == 0): print( print_prefix, "full_auto_set_stepsizes adjusted %s to %f" % (key+"_"+suffix, movement_args[key+"_"+suffix])) break - + #DOC \item **else**: (this is the first time) else: # this is the first time first_time = False @@ -2172,7 +2186,7 @@ def full_auto_set_stepsizes(walkers, walk_stats, movement_args, comm, Emax, KEma print( print_prefix, "full_auto_set_stepsizes adjusted %s to %f" % (key+"_"+suffix, movement_args[key+"_"+suffix])) break - #DOC \item Return step sizes and time taken for routine to run + #DOC \item **return** step sizes and time taken for routine to run print_prefix = orig_prefix @@ -3431,13 +3445,13 @@ def main(): movement_args['n_model_calls'] = int(args.pop('n_model_calls', 0)) movement_args['do_good_load_balance'] = str_to_logical(args.pop('do_good_load_balance', "F")) - #DOC \item process n\_atom\_steps - #DOC \item If break\_up\_atom\_traj - #DOC \item n\_atom\_steps\_per\_call = 1 - #DOC \item n\_atom\_steps\_n\_calls = n\_atom\_steps - #DOC \item else - #DOC \item n\_atom\_steps\_per\_call = n\_atom\_steps - #DOC \item n\_atom\_steps\_n\_calls = 1 + #DOC \item process n\_atom\_steps + #DOC \item **if** break\_up\_atom\_traj: + #DOC \item n\_atom\_steps\_per\_call = 1 + #DOC \item n\_atom\_steps\_n\_calls = n\_atom\_steps + #DOC \item **else**: + #DOC \item n\_atom\_steps\_per\_call = n\_atom\_steps + #DOC \item n\_atom\_steps\_n\_calls = 1 movement_args['n_atom_steps'] = int(args.pop('n_atom_steps', 1)) movement_args['atom_traj_len'] = int(args.pop('atom_traj_len', 8)) From 611791f8a2b85221756d22b68ab228f91240ac38 Mon Sep 17 00:00:00 2001 From: Vincent Fletcher Date: Mon, 9 Sep 2024 17:39:22 +0100 Subject: [PATCH 3/6] tidied up nearest neighbour implementation --- ns_run.py | 56 ++++++++++++++++--------------------------------------- 1 file changed, 16 insertions(+), 40 deletions(-) diff --git a/ns_run.py b/ns_run.py index 66a99a0..db53f97 100755 --- a/ns_run.py +++ b/ns_run.py @@ -1,6 +1,6 @@ #This is a modified version of the pymatnest ns_run.py file, with an additional criteria to reject #configurations if they contain bond lengths less than a given bond length defined by the -#'min_bond_len' keyword +#'min_nn_dis' keyword #Modification made by @V.G.Fletcher @@ -1097,28 +1097,19 @@ def do_MD_atom_walk(at, movement_args, Emax, KEmax): reject_KEmax = (KEmax > 0.0 and final_KE >= KEmax) #################################VGF_MODIFIED######################################### - if ns_args['calc_distances']: + if ns_args['calc_nn_dis']: distances = ase.geometry.get_distances(at.get_positions(), pbc=True, cell=at.get_cell()) bonds = np.sort(distances[1], axis=None) config_size = len(at) min_bond = bonds[config_size] - reject_len = min_bond < ns_args['min_bond_len'] - - num_bonds= bonds[config_size:] < (ns_args['min_bond_len'] + 0.1) - perc_holes = sum(num_bonds)/len(num_bonds) - - reject_pot_hole = perc_holes > ns_args['hole_tolerance'] - - if reject_pot_hole and not (reject_fuzz or reject_Emax or reject_KEmax or reject_len): - ase.io.write(f"potential_holes.traj.{rank}.extxyz", at, parallel=False, format=ns_args['config_file_format'], append=True) + reject_len = min_bond < ns_args['min_nn_dis'] else: reject_len = False - reject_pot_hole = False ##################################################################################### #DOC \item if reject - if reject_fuzz or reject_Emax or reject_KEmax or reject_len or reject_pot_hole: # reject + if reject_fuzz or reject_Emax or reject_KEmax or reject_len: # reject #DOC \item set positions, velocities, energy back to value before perturbation (maybe should be after?) # print(print_prefix, ": WARNING: reject MD traj Emax ", Emax, " initial E ", orig_E, " velo perturbed E ", pre_MD_E, " final E ",final_E, " KEmax ", KEmax, " KE ", final_KE) at.set_positions(pre_MD_pos) @@ -1481,21 +1472,15 @@ def do_cell_step(at, Emax, p_accept, transform): # set new positions and velocities at.set_cell(new_cell, scale_atoms=True) ################################################VGF-MODIFIED####################################### - if ns_args['calc_distances']: + if ns_args['calc_nn_dis']: distances = ase.geometry.get_distances(at.get_positions(), pbc=True, cell=at.get_cell()) bonds = np.sort(distances[1], axis=None) config_size = len(at) min_bond = bonds[config_size] - accept_len = min_bond > ns_args['min_bond_len'] - - num_bonds= bonds[config_size:] < (ns_args['min_bond_len'] + 0.1) - perc_holes = sum(num_bonds)/len(num_bonds) - - accept_pot_hole = perc_holes < ns_args['hole_tolerance'] + accept_len = min_bond > ns_args['min_nn_dis'] else: accept_len = True - accept_pot_hole = True ################################################################################################### if Emax is None: return @@ -1511,12 +1496,10 @@ def do_cell_step(at, Emax, p_accept, transform): #print("error in eval_energy setting new_energy = 2*abs(Emax)=" , new_energy) # accept or reject - if new_energy < Emax and accept_pot_hole: # accept + if new_energy < Emax: # accept at.info['ns_energy'] = new_energy return True else: # reject and revert - if new_energy < Emax: - ase.io.write(f"potential_holes.traj.{rank}.extxyz", at, parallel=False, format=ns_args['config_file_format'], append=True) at.set_cell(orig_cell,scale_atoms=False) at.set_positions(orig_pos) if ns_args['n_extra_data'] > 0: @@ -3307,10 +3290,10 @@ def main(): ns_args['n_extra_data'] = int(args.pop('n_extra_data', 0)) ns_args['Z_cell_axis'] = float(args.pop('Z_cell_axis', 10.0)) #VGF parameters - ns_args['min_bond_len'] = float(args.pop('min_bond_len', 0.0)) - ns_args['hole_tolerance'] = float(args.pop('hole_tolerance',1.0)) - - ns_args['calc_distances'] = (ns_args['min_bond_len'] != 0.0) or (ns_args['hole_tolerance'] != 1.0) + ns_args['min_nn_dis'] = float(args.pop('min_nn_dis', 0.0)) + ns_args['calc_nn_dis_init'] = ns_args['min_nn_dis'] > 0.0 + + ns_args['calc_nn_dis'] = str_to_logical(args.pop('calc_nn_dis', 'F')) #End VGF parameters # surely there's a cleaner way of doing this? @@ -3889,35 +3872,28 @@ def main(): n_try = 0 ###################################VGF-MODIFIED############################### reject_len = True - reject_pot_hole = True - while (n_try < ns_args['random_init_max_n_tries']) and (((math.isnan(energy) or energy > ns_args['start_energy_ceiling'])) or reject_len or reject_pot_hole): + while (n_try < ns_args['random_init_max_n_tries']) and (((math.isnan(energy) or energy > ns_args['start_energy_ceiling'])) or reject_len): at.set_scaled_positions( rng.float_uniform(0.0, 1.0, (len(at), 3) ) ) if movement_args['2D']: # zero the Z coordiates in a 2D simulation temp_at=at.get_positions() temp_at[:,2]=0.0 at.set_positions(temp_at) ################################################################## - if ns_args['calc_distances']: + if ns_args['calc_nn_dis_init']: distances = ase.geometry.get_distances(at.get_positions(), pbc=True, cell=at.get_cell()) bonds = np.sort(distances[1], axis=None) config_size = len(at) min_bond = bonds[config_size] #print('min_bond_here', min_bond) - reject_len = min_bond < ns_args['min_bond_len'] - - num_bonds= bonds[config_size:] < (ns_args['min_bond_len'] + 0.1) - perc_holes = sum(num_bonds)/len(num_bonds) - - reject_pot_hole = perc_holes > ns_args['hole_tolerance'] + reject_len = min_bond < ns_args['min_nn_dis'] else: reject_len = False - reject_pot_hole = False ##################################################################### energy = eval_energy(at) n_try += 1 - if math.isnan(energy) or energy > ns_args['start_energy_ceiling'] or reject_len or reject_pot_hole: + if math.isnan(energy) or energy > ns_args['start_energy_ceiling'] or reject_len: sys.stderr.write("WARNING: rank %d failed to generate initial config by random positions under max energy %f in %d tries\n" % (rank, ns_args['start_energy_ceiling'], ns_args['random_init_max_n_tries'])) # try FORTRAN config initializer @@ -3944,7 +3920,7 @@ def main(): at.info['ns_energy'] = rand_perturb_energy(energy, ns_args['random_energy_perturbation']) at.info['volume'] = at.get_volume() ###VGF modified### - if ns_args['calc_distances']: + if ns_args['calc_nn_dis_init']: distances = ase.geometry.get_distances(at.get_positions()) bonds = np.sort(distances[1], axis=None) min_bond = bonds[len(at)] From e1ed8d4906e9460834f9485a3e35e7916442142b Mon Sep 17 00:00:00 2001 From: Vincent Fletcher Date: Mon, 9 Sep 2024 17:42:55 +0100 Subject: [PATCH 4/6] modified tag --- ns_run.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/ns_run.py b/ns_run.py index db53f97..99f6f07 100755 --- a/ns_run.py +++ b/ns_run.py @@ -1,9 +1,3 @@ -#This is a modified version of the pymatnest ns_run.py file, with an additional criteria to reject -#configurations if they contain bond lengths less than a given bond length defined by the -#'min_nn_dis' keyword - -#Modification made by @V.G.Fletcher - import re, math, time, os import pprint import numpy as np, ase, ase.io From 536a3b4c859539403c892664f1655b33ffdeed25 Mon Sep 17 00:00:00 2001 From: Vincent Fletcher Date: Mon, 9 Sep 2024 18:32:47 +0100 Subject: [PATCH 5/6] added comments --- ns_run.py | 35 +++++++++++++++++++++++------------ 1 file changed, 23 insertions(+), 12 deletions(-) diff --git a/ns_run.py b/ns_run.py index 99f6f07..2b80185 100755 --- a/ns_run.py +++ b/ns_run.py @@ -479,7 +479,14 @@ def usage(): | Analyzers to apply during run. String consists of semicolon separated pairs of module name and intervals. Module names correspond to analysis modules in NS_PATH/ns_run_analyzers (see there for examples) or elsewhere in PYTHONPATH | Positive interval refers to NS loop, negative to initial walks | default:'' - + + ``min_nn_dis=float`` + | Optional input to add additional rejection criteria when initially generating walkers. If any nearest neighbour distances are less than this length, given in Angstroms, the walker is rejected. + | default: 0.0 + + ``calc_nn_dis=[ T | F ]`` + | While you can set the min_nn_dis to add an additional rejection criteria to the initial walker generation, you can use this keyword to keep the rejection criteria on throughout the sampling. This does increase the calculation time of a run and scales exponentially with the number of atoms. + | default: F """ sys.stderr.write("Usage: %s [ -no_mpi ] < input\n" % sys.argv[0]) sys.stderr.write("input:\n") @@ -619,6 +626,8 @@ def usage(): sys.stderr.write("no_extra_walks_at_all=[ T | F ] (F)\n") sys.stderr.write("track_configs=[ T | F ] (F)\n") sys.stderr.write("track_configs_write=[ T | F ] (F)\n") + sys.stderr.write("min_nn_dis=float, (0.0, reject initial walkers with a nearest neighbour distance less than this value in Angstroms)") + sys.stderr.write("calc_nn_dis=[ T | F ], (F, reject walkers with a nearest neighbour distance less than min_nn_dis throughout the entire sampling procedure)") def excepthook_mpi_abort(exctype, value, tb): print( print_prefix,'Uncaught Exception Type:', exctype) @@ -1090,7 +1099,7 @@ def do_MD_atom_walk(at, movement_args, Emax, KEmax): reject_Emax = (final_E >= Emax) reject_KEmax = (KEmax > 0.0 and final_KE >= KEmax) - #################################VGF_MODIFIED######################################### + #VGF calculate nn_distances and reject if enabled, otherwise skip calculation and dont reject based on nn distance if ns_args['calc_nn_dis']: distances = ase.geometry.get_distances(at.get_positions(), pbc=True, cell=at.get_cell()) bonds = np.sort(distances[1], axis=None) @@ -1100,7 +1109,7 @@ def do_MD_atom_walk(at, movement_args, Emax, KEmax): reject_len = min_bond < ns_args['min_nn_dis'] else: reject_len = False - ##################################################################################### + #VGF end #DOC \item if reject if reject_fuzz or reject_Emax or reject_KEmax or reject_len: # reject @@ -1465,7 +1474,8 @@ def do_cell_step(at, Emax, p_accept, transform): # set new positions and velocities at.set_cell(new_cell, scale_atoms=True) - ################################################VGF-MODIFIED####################################### + + #VGF calculate nn_distances and reject if enabled, otherwise skip calculation and dont reject based on nn distance if ns_args['calc_nn_dis']: distances = ase.geometry.get_distances(at.get_positions(), pbc=True, cell=at.get_cell()) bonds = np.sort(distances[1], axis=None) @@ -1475,7 +1485,8 @@ def do_cell_step(at, Emax, p_accept, transform): accept_len = min_bond > ns_args['min_nn_dis'] else: accept_len = True - ################################################################################################### + #VGF end + if Emax is None: return @@ -3283,12 +3294,12 @@ def main(): ns_args['random_energy_perturbation'] = float(args.pop('random_energy_perturbation', 1.0e-12)) ns_args['n_extra_data'] = int(args.pop('n_extra_data', 0)) ns_args['Z_cell_axis'] = float(args.pop('Z_cell_axis', 10.0)) - #VGF parameters + #VGF nn parameters ns_args['min_nn_dis'] = float(args.pop('min_nn_dis', 0.0)) ns_args['calc_nn_dis_init'] = ns_args['min_nn_dis'] > 0.0 ns_args['calc_nn_dis'] = str_to_logical(args.pop('calc_nn_dis', 'F')) - #End VGF parameters + #VGF end nn parameters # surely there's a cleaner way of doing this? try: @@ -3864,7 +3875,7 @@ def main(): # random initial positions energy = float('nan') n_try = 0 - ###################################VGF-MODIFIED############################### + #VGF added optional nn rejection criteria reject_len = True while (n_try < ns_args['random_init_max_n_tries']) and (((math.isnan(energy) or energy > ns_args['start_energy_ceiling'])) or reject_len): at.set_scaled_positions( rng.float_uniform(0.0, 1.0, (len(at), 3) ) ) @@ -3872,7 +3883,7 @@ def main(): temp_at=at.get_positions() temp_at[:,2]=0.0 at.set_positions(temp_at) - ################################################################## + #VGF calculate nn_distances and reject if enabled, otherwise skip calculation and dont reject based on nn distance if ns_args['calc_nn_dis_init']: distances = ase.geometry.get_distances(at.get_positions(), pbc=True, cell=at.get_cell()) bonds = np.sort(distances[1], axis=None) @@ -3883,7 +3894,7 @@ def main(): reject_len = min_bond < ns_args['min_nn_dis'] else: reject_len = False - ##################################################################### + #VGF end energy = eval_energy(at) n_try += 1 @@ -3913,13 +3924,13 @@ def main(): energy = eval_energy(at) at.info['ns_energy'] = rand_perturb_energy(energy, ns_args['random_energy_perturbation']) at.info['volume'] = at.get_volume() - ###VGF modified### + #VGF if nn rejection criteria enabled, print the final smallest nearest neighbour distance if ns_args['calc_nn_dis_init']: distances = ase.geometry.get_distances(at.get_positions()) bonds = np.sort(distances[1], axis=None) min_bond = bonds[len(at)] print('min_bond_final', rank, min_bond) - ################### + #VGF end # Done initialising atomic positions. Now initialise momenta From 4f120befe9680ae16ab81e3bf905d0824b71dab1 Mon Sep 17 00:00:00 2001 From: VGFletcher <125897513+VGFletcher@users.noreply.github.com> Date: Mon, 9 Sep 2024 19:35:23 +0100 Subject: [PATCH 6/6] updated comment --- ns_run.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ns_run.py b/ns_run.py index 2b80185..2dfcadc 100755 --- a/ns_run.py +++ b/ns_run.py @@ -905,7 +905,7 @@ def rotate_dir_3N_2D(vec, max_ang): def rej_free_perturb_velo(at, Emax, KEmax, rotate=True): #DOC -#DOC ``rej_free_perturb_velo(at, Emax, KEmax, rotate=True)`` +#DOC ``rej_free_perturb_velo`` if not at.has('momenta') or not at.has('masses'): return