import re, math, time, os
import pprint
import numpy as np, ase, ase.io
import ns_rng
import stacktrace
from copy import deepcopy
import pick_interconnected_clump
try:
import matscipy.neighbours
except:
pass
import collections
from traceback import print_exception
import check_memory
from ase.md.verlet import VelocityVerlet
import importlib
print_prefix=""
[docs]
def usage():
""" Print help to the standard output about the usage of the code and input parameters. The current list of parameters is the following:
.. glossary::
``max_volume_per_atom=float``
| Maximum volume per atom allowed during the run.
| default: 1.0e3
``min_volume_per_atom=float``
| Minimum volume per atom allowed during the run.
| default: 1.0
``start_species=int int [ float ] [, int int [ float ] ... ]``
| MANDATORY (or start_config_file)
| Atomic number; multiplicity; [ not recommended: mass (amu) ]. Info repeated for each species, separated by commas, mass is optional and not recommended.
``start_config_file=str``
| MANDATORY (or start_species)
| Name of file to read in for atom information
| default: ''
``keep_atoms_fixed=int ``
| Number of atoms to fix in surface simulation.
``apply_Z_wall=[T|F] ``
| Whether to have a boundary in the Z direction to keep free particles reaching the other side of "surface" layer due to periodic boundary conditions.
| This functionality is not fully tested! Recommended use with MC evaluator and fortran. If constructing a surface layer, make it parallel to the XY plane and set the X dimension with the "wall" taken into account, no atoms should violate the wall restriction initially!!!
| If True, the wall is set at 10 angstroms below the top of the simulation cell.
``restart_file=path_to_file``
| File for restart configs. Mutually exclusive with ``start_*``, one is required. The file should contain the state of the walkers to continue from along with the restart iteration number. Normally such a file can be the concatenated snapshot files.
``n_walkers=int``
| MANDATORY
| Total number of walkers, i.e. the size of the live set used. It has to be the multiple of the number of processors used and it also has to be larger than the number of processors.
``n_cull=int``
| Number of walkers to kill at each NS iteration. Use of default 1 is strongly recommended.
| deafult: 1
``n_extra_walk_per_task=int``
| default: 0
``n_iter_times_fraction_killed=int``
| MANDATORY
| Number of nested sampling iteration cycles performed per walker. Thus the total number of iterations will be ``n_iter_times_fraction_killed`` / ``(n_cull/n_walkers)``. Either this or ``converge_down_to_T`` is required.
``converge_down_to_T=flot``
| MANDATORY
| temperature down to which Z(T) should be converged. Either this or ``n_iter_times_fraction_killed`` is required.
``T_estimate_finite_diff_lag=int``
| default: 1000
| Lag (in iterations) for doing d log Omega / d E finite difference derivative
``min_Emax=float``
| Termination condition based on Emax: if this value is reached, the iteration will stop.
| No default.
``out_file_prefix=str``
| String used as prefix for the different output files.
| No default.
``energy_calculator= ( ASE | lammps | internal | fortran)``
| Energy calculator.
| default: fortran
``n_extra_data=int``
| Amount of extra data per atom to pass around.
| default: 0
``KEmax_max_T=float``
| If > 0, maximum temperature for estimating KEmax. For constant *P* ensemble 3/2 P Vmax will be used if KEmax_max_T is < 0.
| default: -1
``kB=float``
| Boltzmann constant
| default: 8.6173324e-5 (eV/A)
``start_energy_ceiling_per_atom=float``
| Maximum potential energy per atom for initial configurations. P*Vmax is added to this automatically in case of NpT runs.
| default: 1.0e9
``start_energy_ceiling=float``
| DEPRECATED: use ``start_energy_ceiling_per_atom``. Maximum potential energy for initial configurations. P*Vmax is added to this automatically in case of NpT runs.
| default: 1.0e9
``random_init_max_n_tries=int``
| Maximum number of tries to create initial random atomic configuration
| default 100
``n_model_calls_expected=int``
| Total number of model calls performed during the random walk, expected by the user. The actual number of model calls performed during runtime can be different: this number is divided up among the processors in case of running parallel, and in order to make sure different processors perform about the same number of calls, it can be increased for some. (It is recommended to use n_model_calls_expected/[number of processors] > 20) Either this or the keyword ``n_model_calls`` is mandatory, but the use of this keyword is strongly recommended. If < 0, value will be set to _sum_ of each type of steps (accounting for atom_traj_len), so ``n_*_steps`` will be setting the _number_ of calls, rather than just the ratios.
| default: 0
``n_model_calls=int``
| Number of model calls performed during the random walk. This is the actual number of calls performed, ideally the program sets its value depending on n_model_calls_expected and the number of processors used. Either this or the keyword ``n_model_calls_expected`` is mandatory, though the usage of the latter is strongly recommended.
| default: 0
``do_good_load_balance=[T | F]``
| Whether to do steps in groups with good load balance
| default: F
``n_atom_steps=int``
| Ratio of atomic trajectories will be determined by ``n_atom_steps``/SUM(``n_*_steps``).
| default: 1
``atom_traj_len=int``
| Length of atomic trajectory (MD steps or MC sweeps) in each atomic type step.
| default: 8
``atom_traj_len_cost_multiplier=int``
| Multiplier for cost of an atomic step (set to 1 for MD, TE-HMC, and SP-MC with O(1) cost moves, N for naive SP-MC)
| default: 1
``break_up_atom_traj=[T | F]``
| Whether to intersperse ``n_atom_steps`` atomic sub-trajectories with other types of steps.
| default: F
``n_cell_volume_steps=int``
| Ratio of cell volume steps will be determined by ``n_cell_volume_steps``/SUM(``n_*_steps``).
| default: 1
``n_cell_shear_steps=int``
| Ratio of cell shear steps will be determined by ``n_cell_shear_steps``/SUM(``n_*_steps``).
| default: 1
``n_cell_stretch_steps=int``
| Ratio of cell stretch steps will be determined by ``n_cell_stretch_steps``/SUM(``n_*_steps``).
| default: 1
``n_swap_steps=int``
| Ratio of species swap steps will be determined by ``n_swap_steps``/SUM(``n_*_steps``). It has to be set other than zero for a multicomponent system.
| default: 0
``swap_max_cluster=int``
| Maximum size of interconnected cluster to try to swap.
| default: 1
``swap_r_cut=float``
| Cutoff radius for defining connected atoms for cluster.
| default: 2.5
``swap_cluster_probability_increment=float``
| Factor between prob. of picking increasing larger clusters.
| default: 0.75
``swap_velo=[T | F]``
| If true, swap velocities when swapping atoms, breaking coherence a bit.
| default: F
``no_swap_velo_fix_mag_alt=[T | F]``
| If true, use alternate method for correcting velocity magnitudes when not swapping velocities.
| default: F
``n_semi_grand_steps=int``
| Ratio of species type-change steps will be determined by ``n_semi_grand_steps``/SUM(``n_*_steps``). Only makes sense for a multicomponent system.
| default: 0
``semi_grand_potentials=Z1 : mu1 [, Z2 : mu2 ...]``
| Chemical potentials for semi-grand canonical species-type transmutations
| default: ''
``velo_traj_len=int``
| Number of MC steps in (optional) explicit velocity MC trajectory.
| default: 0
``random_energy_perturbation=float``
| default: 1.0e-12
``atom_algorithm=[MC | MD | GMC]``
| MANDATORY
| Use either Monte Carlo or Molecular dynamics to explore. GMC is an alias for atom_algorithm=MC, MC_atom_Galilean=True.
``MC_atom_velocities=[T | F]``
| This keyword is supported only for energy_calculator=fortran.
| default: F
``MC_atom_velocities_pre_perturb=[T | F]``
| Perturb velocities (rejection free) before MC + velocities walk.
| default: F
``MC_atom_step_size=float``
| Initial atom move step size in units of (max_volume_per_atom * N_atoms)^(1/3)
| default: 1.0
``MC_atom_step_size_max=float``
| Maximum atom step size in units of (max_volume_per_atom * N_atoms)^(1/3).
| default: 1.0
``MC_atom_uniform_rv=[T | F]``
| default: F
``MC_atom_Galilean=[T | F]``
| default: F
``GMC_no_reverse=[T | F]``
| default: T
``MD_atom_velo_pre_perturb=[T | F]``
| Perturb velocities before MD trajectory
| default: F
``MD_atom_velo_post_perturb=[T | F]``
| Perturb velocities after MD trajectory
| default: T
``MD_atom_velo_flip_accept=[T | F]``
| default: F
``MD_atom_timestep=float``
| default: 0.1 (ASE time units)
``MD_atom_timestep_max=float``
| default: 2.0 (ASE time units)
``MD_atom_energy_fuzz=float``
| Tolerance for rejecting non-energy conserving trajectories, as fraction of Kinetic Energy
| default: 1.0e-2
``MD_atom_reject_energy_violation=[ T | F ]``
| Use energy conservation violation (exceeding MD_atom_energy_fuzz * KE) to reject MD trajectories.
| default: F
``python_MD=[ T | F ]``
| Do MD using python code rather than underlying driver
| default: F
``atom_velo_rej_free_fully_randomize=[T | F]``
| If true, randomize velocities completely rather than just perturbing.
| default: F
``atom_velo_rej_free_perturb_angle=float``
| Max angle in radians for random rotations.
| default: 0.3
``MC_atom_velo_step_size=float``
| default: 50.0
``MC_atom_velo_step_size_max=float``
| default: 10000.0
``MC_atom_velo_walk_rej_free=[T | F]``
| default: T. If true, use rejection free algorithm for MC_atom_walk
``MC_cell_P=float``
| Pressure value to be used. The unit of pressure depends on both the energy calculator and on the potential model used. Note that ASE uses
| eV as energy and Angstrom as distance units everywhere, thus this case the pressure in the input have to be eV/Angstrom^3 (in case of using LAMMPS
| the lammpslib module will convert the units for lammps to whatever units are set by the potential type)
| default: 0.0
``MC_cell_flat_V_prior=[T | F]``
| Flat prior for V (use with MC_cell_P=0 and cell moves, requires reweighting configurations when analyzing)
| POORLY TESTED, DO NOT TRUST (YET).
| default: F
``MC_cell_volume_per_atom_step_size=float``
| Initial volume stepsize for volume change.
| default: 5% of the maximum allowed volume
``MC_cell_volume_per_atom_step_size_max=float``
| Maximum allowed volume step size.
| default: 50% of the maximum allowed volume
``MC_cell_volume_per_atom_prob=float``
| default: 1.0
``MC_cell_stretch_step_size=float``
| default: 0.1
``MC_cell_stretch_step_size_max=float``
| default: 1.0
``MC_cell_stretch_prob=float``
| default: 1.0
``MC_cell_shear_step_size=float``
| default: 0.1, in units of (max_volume_per_atom * N_atoms)^(1/3)
``MC_cell_shear_step_size_max=float``
| default: 1.0, in units of (max_volume_per_atom * N_atoms)^(1/3)
``MC_cell_shear_prob=float``
| default: 1.0
``MC_cell_min_aspect_ratio=float``
| Smallest allowed distance between parallel faces for cell normalised to unit volume. A higher value of ``MC_cell_min_aspect_ratio`` restricts the system to more cube-like cell shapes, while a low value allows the system to become essentially flat. In case of 64 atoms the use of ``MC_cell_min_aspect_ratio`` < 0.65 *does* effect the melting transition.
| default: 0.8
``cell_shape_equil_steps=int``
| default: 1000
``full_auto_step_sizes=[T | F]``
| If true (T), automatically calibrate all sizes by performing additional short explorations, including at the start of run. If false (F), use initial input step sizes and make small adjustments to these during the run.
| default: T
``monitor_step_interval_times_fraction_killed=float``
| Divided by ``n_cull/n_walkers`` to get actual monitoring interval in iterations, negative for only using last iteration, 0 for no monitoring
| default: 1
``adjust_step_interval_times_fraction_killed=float``
| Divided by ``n_cull/n_walkers`` to get actual step adjustment interval in iterations, negative for only using last iteration, 0 for no adjust.
| default: 1
``MC_adjust_step_factor=float``
| default: 1.1
``MC_adjust_min_rate=float``
| default: 0.25
``MC_adjust_max_rate=float``
| default: 0.75
``GMC_adjust_min_rate=float``
| default: 0.25
``GMC_adjust_max_rate=float``
| default: 0.75
``GMC_dir_perturb_angle=float``
| default: -1.0
``GMC_dir_perturb_angle_during=float``
| default: 0.0
``MD_adjust_step_factor=float``
| default: 1.1
``MD_adjust_min_rate=float``
| default: 0.5
``MD_adjust_max_rate=float``
| default: 0.95
``ASE_calc_module=str``
| MANDATORY if energy_calculator=ASE
``FORTRAN_model=str``
| MANDATORY if energy_calculator=fortran
``FORTRAN_model_params=str``
| parameters (optional) for fortran model
``LAMMPS_fix_gmc=[T | F]``
| default: F
``LAMMPS_init_cmds=str``
| MANDATORY if energy_calculator=lammps
``LAMMPS_name=str``
| '', arch name for lammps shared object file
``LAMMPS_serial=[T | F]``
| default True, lammps library is serial so do not pass COMM_SELF as a communicator
``LAMMPS_header=str``
| lammpslib.py value default, override lammpslib.py header commands for energy_calculator=lammps
``LAMMPS_header_extra=str``
| '', extra lammpslib.py header commands for energy_calculator=lammps
``LAMMPS_atom_types=str``
| MANDATORY if energy_calculator=lammps
| atomic_symbol1 lammps_type1 [, atomic_symbol2 lammps_type2, ...]
| mapping between atom species and lammps potential types
``config_file_format=str``
| File format in which configurations are printed, e.g. for trajectory and snapshot files.
| default: extxyz
``rng=( numpy | internal | rngstream )``
| Random number generator.
| default: numpy
``profile=rank_to_profile``
| default: -1
``2D=[ T | F ]``
| Perform 2D simulation. Some functionality may be limited with this option!
| The 2D area is in the xy plane, the z axis stays constant during the simulation (set its value using the keyword ``Z_cell_axis``).
| The maximum allowed volume is still defined as volume in this case!
| If constant pressure is set, the p*V term is now p*A*Z_cell_axis, thus set the MC_cell_P parameter multiplied by 1/Z_cell_axis.
| default: F
``Z_cell_axis=float``
| Only used if 2D=T. Value of the Z cell axis, perpendicular to the area in case of a 2D simulation.
| This stays constant during the simulation.
| default: 10.0
``debug=int``
| Verbosity level used in the output file. The larger its value the more info is printed.
| default: 0
``snapshot_interval=int``
| Iteration interval at which a snapshot is created: every process prints out its current walkers in extended xyz format. If it is set <=0, no snapshots will be printed except the final positions at the end of the nested sampling run. Note that when new snapshots are printed, the previous set is deleted. The snapshot files are convenient source to see how the sampling progresses, but these are also the basis to restart a sampling! When using restart, the walkers will be read from these files.
| default: 10000
``snapshot_time=int``
| Max time between snapshots in seconds
| default: 3600
``snapshot_per_parallel_task=[ T | F ]``
| Save a separate snapshot file for each parallel task. Faster I/O (in parallel), but less convenient in some cases.
| default: T
``snapshot_clean=[ T | F ]``
| If true, delete previous iteration snapshot files
| default: T
``random_initialise_pos=[ T | F ]``
| If true, randomize the initial positions
| default: T
``random_initialise_cell=[ T | F ]``
| If true, randomize the initial cell (currently by random walk)
| default: T
``LAMMPS_molecular_info=[ T | F ]``
| If true, create lammps bonds and other molecular MM features from initial atoms config (e.g. can be read from a file)
| default: T
``initial_walk_N_walks=int``
| number of initial random walks to apply to all walkers
| default: ''
``initial_walk_adjust_interval=int``
| in initial walks, interval between walks for adjusting step size
| default: ''
``initial_walk_Emax_offset_per_atom=float``
| offset (per atom) to increase Emax during initial random walk applied to all walkers
| default: 1.0
``initial_walk_only=[T | F]``
| do initial walk only, then quit
| default: F
``traj_interval=int``
| Iteration interval at which the currently culled configuration(s) is/are printed to the trajectory output, in the set format. If it is set <=0, no trajectory files will be printed at all. Useful option for larger runs as the trajectory files can become huge.
| default: 100
``delta_random_seed=int``
| Random number seed to be used in the run. If smaller than 0, a seed from /dev/urandom is used.
| default: -1
``no_extra_walks_at_all=[ T | F ]``
| default: F
``track_configs=[ T | F ]``
| Track configrations across all walks/clones.
| default: F
``track_configs_write=[ T | F ]``
| Write tracked configrations to a file (if false, tracking info will still be in snapshot and traj files)
| default: F
``ns_run_analyzers=string``
| 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")
sys.stderr.write("max_volume_per_atom=float (1e3)\n")
sys.stderr.write("min_volume_per_atom=float (1)\n")
sys.stderr.write("start_species=int int [ float ] [, int int [ float ] ... ] (MANDATORY, this or start_config_file required. atomic_number multiplicity [not recomended: mass (amu)]. Info repeated for each species, separated by commas, mass is optional and not recommended.\n")
sys.stderr.write("start_config_file=str (MANDATORY, this or start_species required. if set filename to read initial atom information from (instead of creating them)\n")
sys.stderr.write("keep_atoms_fixed=int (0, no atoms to be fixed)\n")
sys.stderr.write("apply_Z_wall=[T | F] (F, supported only for fortran MC)\n")
sys.stderr.write("restart_file=path_to_file (file for restart configs. Mutually exclusive with start_*, one is required)\n")
sys.stderr.write("n_walkers=int (MANDATORY)\n")
sys.stderr.write("n_cull=int (1, number of walkers to kill at each NS iteration)\n")
sys.stderr.write("n_extra_walk_per_task=int (0)\n")
sys.stderr.write("n_iter_times_fraction_killed=int (MANDATORY, this or converge_down_to_T required)\n")
sys.stderr.write("converge_down_to_T=float (MANDATORY, this or n_iter_times_fraction_killed required)\n")
sys.stderr.write("T_estimate_finite_diff_lag=int (1000, lag for doing finite difference in current T estimate\n")
sys.stderr.write("min_Emax=float (None. Termination condition based on Emax)\n")
sys.stderr.write("out_file_prefix=str (None)\n")
sys.stderr.write("energy_calculator= ( ASE | lammps | internal | fortran) (fortran)\n")
sys.stderr.write("n_extra_data=int (0, amount of extra data per atom to pass around)\n")
sys.stderr.write("ns_run_analyzers=str ('', analyzers to apply during run\n")
sys.stderr.write("\n")
sys.stderr.write("KEmax_max_T=float (-1, if > 0 maximum temperature for estimating KEmax)\n")
sys.stderr.write("kB=float (8.6173324e-5, Boltzmann constant)\n")
sys.stderr.write("start_energy_ceiling_per_atom=float (1.0e9, max potential energy per atom for initial configs. P*Vmax is added to this automatically)\n")
sys.stderr.write("start_energy_ceiling=float (DEPRECATED: use start_energy_ceiling_per_atom. 1.0e9, max potential energy for initial configs. P*Vmax is added to this automatically)\n")
sys.stderr.write("random_init_max_n_tries=int (100, maximum number of tries for random initial positions under energy ceiling\n")
sys.stderr.write("\n")
sys.stderr.write("n_model_calls_expected=int (0, one of these is required)\n")
sys.stderr.write("n_model_calls=int (0, one of these is required)\n")
sys.stderr.write("do_good_load_balance=[T | F] (F, whether to do steps in groups with good load balance\n")
sys.stderr.write("\n")
sys.stderr.write("n_atom_steps=int (1, number of atomic trajectories \n")
sys.stderr.write("atom_traj_len=int (8, length of atomic trajectory (MD steps or MC sweeps) in each step\n")
sys.stderr.write("atom_traj_len_cost_multiplier=int (1, multiplier for cost of an atomic step (set to 1 for MD, TE-HMC, and SP-MC with O(1) cost moves, N for naive SP-MC)\n")
sys.stderr.write("break_up_atom_traj=[T | F] (F, whether to intersperse n_atom_steps atomic sub-trajectories with other types of steps\n")
sys.stderr.write("\n")
sys.stderr.write("n_cell_volume_steps=int (1, number of cell MC volume steps )\n")
sys.stderr.write("n_cell_shear_steps=int (1, number of cell MC shear steps )\n")
sys.stderr.write("n_cell_stretch_steps=int (1, number of cell MC stretch steps )\n")
sys.stderr.write("n_swap_steps=int (0, number of atom swaps)\n")
sys.stderr.write("swap_max_cluster=int (1, maximum size of interconnected cluster to try to swap)\n")
sys.stderr.write("swap_r_cut=float (2.5, cutoff radius for defining connected atoms for cluster)\n")
sys.stderr.write("swap_cluster_probability_increment=float (0.75, factor between prob. of picking increasing larger clusters)\n")
sys.stderr.write("swap_velo=[T | F] (F, if true, swap velocities when swapping atoms, breaking coherence a bit\n")
sys.stderr.write("swap_velo_fix_mag_alt=[T | F] (F, if true, use alternate correction for velocity magnitudes when not swapping velocities\n")
sys.stderr.write("\n")
sys.stderr.write("n_semi_grand_steps=int (0, number of species type changes)\n")
sys.stderr.write("semi_grand_potentials={ int : float [ , int : float ... ] } (None, chemical potentials for species type changes)\n")
sys.stderr.write("\n")
sys.stderr.write("velo_traj_len=int (0, number of MC sweeps in each velocity MC segement)\n")
sys.stderr.write("\n")
sys.stderr.write("random_energy_perturbation=float (1.0e-12)\n")
sys.stderr.write("atom_algorithm=[MC | MD | GMC] (MANDATORY)\n")
sys.stderr.write("\n")
sys.stderr.write("MC_atom_velocities=[T | F] (F, supported only for energy_calculator=fortran)\n")
sys.stderr.write("MC_atom_velocities_pre_perturb=[T | F] (F, Perturb velocities (rejection free) before MC + velocities walk)\n")
sys.stderr.write("MC_atom_step_size=float (1.0, in units of (max_volume_per_atom * N_atoms)^(1/3) )\n")
sys.stderr.write("MC_atom_step_size_max=float (1.0, in units of (max_volume_per_atom * N_atoms)^(1/3) )\n")
sys.stderr.write("MC_atom_uniform_rv=[T | F] (F)\n")
sys.stderr.write("MC_atom_Galilean=[T | F] (F)\n")
sys.stderr.write("GMC_no_reverse=[T | F] (T)\n")
sys.stderr.write("\n")
sys.stderr.write("MD_atom_velo_pre_perturb=[T | F] (F. Perturb velocities before MD trajectory\n")
sys.stderr.write("MD_atom_velo_post_perturb=[T | F] (T. Perturb velocities after MD trajectory\n")
sys.stderr.write("MD_atom_velo_flip_accept=[T | F] (F)\n")
sys.stderr.write("MD_atom_timestep=float (0.1, in ASE time units)\n")
sys.stderr.write("MD_atom_timestep_max=float (2.0, in ASE time units)\n")
sys.stderr.write("MD_atom_energy_fuzz=float (1.0e-2. Tolerance for rejecting non-energy conserving trajectories, as fraction of KE)\n")
sys.stderr.write("MD_atom_reject_energy_violation=[ T | F ] (F, use energy conservation violation (exceeding MD_atom_energy_fuzz * KE) to reject MD trajectories)\n")
sys.stderr.write("\n")
sys.stderr.write("atom_velo_rej_free_fully_randomize=[T | F] (F. If true, randomize velocities completely rather than just perturbing.\n")
sys.stderr.write("atom_velo_rej_free_perturb_angle=float (0.3. Max angle in radians for random rotations.)\n")
sys.stderr.write("MC_atom_velo_step_size=float (50.0)\n")
sys.stderr.write("MC_atom_velo_step_size_max=float (10000.0)\n")
sys.stderr.write("MC_atom_velo_walk_rej_free=[T | F] (T. If true, use rejection free algorithm for MC_atom_walk\n")
sys.stderr.write("\n")
sys.stderr.write("\n")
sys.stderr.write("MC_cell_P=float (0.0)\n")
sys.stderr.write("MC_cell_flat_V_prior=[T | F] (F). POORLY TESTED, DO NOT TRUST (YET)\n")
sys.stderr.write("MC_cell_volume_per_atom_step_size=float (5% of the maximum allowed volume)\n")
sys.stderr.write("MC_cell_volume_per_atom_step_size_max=float (50% of the maximum allowed volume)\n")
sys.stderr.write("MC_cell_volume_per_atom_prob=float (1.0)\n")
sys.stderr.write("MC_cell_stretch_step_size=float (0.1)\n")
sys.stderr.write("MC_cell_stretch_step_size_max=float (1.0)\n")
sys.stderr.write("MC_cell_stretch_prob=float (1.0)\n")
sys.stderr.write("MC_cell_shear_step_size=float (0.1 in units of (max_volume_per_atom * N_atoms)^(1/3))\n")
sys.stderr.write("MC_cell_shear_step_size_max=float (1.0 in units of (max_volume_per_atom * N_atoms)^(1/3))\n")
sys.stderr.write("MC_cell_shear_prob=float (1.0)\n")
sys.stderr.write("MC_cell_min_aspect_ratio=float (0.8)\n")
sys.stderr.write("cell_shape_equil_steps=int (1000)\n")
sys.stderr.write("\n")
sys.stderr.write("monitor_step_interval_times_fraction_killed=float (1, divided by n_cull/n_walkers to get actual monitoring interval in iterations, negative for only using last iteration, 0 for no monitoring)\n")
sys.stderr.write("adjust_step_interval_times_fraction_killed=float (5, divided by n_cull/n_walkers to get actual adjustment interval in iterations, negative for only using last iteration, 0 for no adjust)\n")
sys.stderr.write("full_auto_step_sizes=[T | F] (T) (T. Automatically calibrate all sizes by performing additional short explorations, including at start of run. F. Use initial input step sizes and make small adjustments to step sizes during run.)\n")
sys.stderr.write("MC_adjust_step_factor=float (1.1)\n")
sys.stderr.write("MC_adjust_min_rate=float (0.25)\n")
sys.stderr.write("MC_adjust_max_rate=float (0.75)\n")
sys.stderr.write("GMC_adjust_min_rate=float (0.25)\n")
sys.stderr.write("GMC_adjust_max_rate=float (0.75)\n")
sys.stderr.write("GMC_dir_perturb_angle=float (-1.0)\n")
sys.stderr.write("GMC_dir_perturb_angle_during=float (0.0)\n")
sys.stderr.write("MD_adjust_step_factor=float (1.1)\n")
sys.stderr.write("MD_adjust_min_rate=float (0.5)\n")
sys.stderr.write("MD_adjust_max_rate=float (0.95)\n")
sys.stderr.write("\n")
sys.stderr.write("ASE_calc_module=str (MANDATORY if energy_calculator=ASE)\n")
sys.stderr.write("FORTRAN_model=str (MANDATORY if energy_calculator=fortran)\n")
sys.stderr.write("FORTRAN_model_params=str (parameters for energy_calculator=fortran)\n")
sys.stderr.write("LAMMPS_fix_gmc=[T | F]\n")
sys.stderr.write("LAMMPS_init_cmds=str (MANDATORY if energy_calculator=lammps)\n")
sys.stderr.write("LAMMPS_name=str ('', arch name for lammps shared object file)\n")
sys.stderr.write("LAMMPS_serial=[T | F] (T, is lammps library serial)\n")
sys.stderr.write("LAMMPS_header=str (lammpslib.py value default, override lammpslib.py header commands for energy_calculator=lammps)\n")
sys.stderr.write("LAMMPS_header_extra=str ('', extra lammpslib.py header commands for energy_calculator=lammps)\n")
sys.stderr.write("LAMMPS_atom_types=symbol int [, symbol int ] ... ('', mapping from atomic symbols to type numbers for LAMMPS ASE interface)\n")
sys.stderr.write("\n")
sys.stderr.write("config_file_format=str (extxyz)\n") # julia
sys.stderr.write("rng=( numpy | internal | rngstream ) (numpy)\n") # julia
sys.stderr.write("profile=rank_to_profile (-1)\n")
sys.stderr.write("2D=[ T | F ] (F)\n")
sys.stderr.write("Z_cell_axis=float (10.0)\n")
sys.stderr.write("debug=debug_level (0, <= 0 for no debugging tests/prints)\n")
sys.stderr.write("snapshot_interval=int (10000, <=0 for no snapshots except final positions)\n")
sys.stderr.write("snapshot_time=int (3600, <=0 for no time based snapshots)\n")
sys.stderr.write("snapshot_per_parallel_task=[T | F] (T, if true save separate snapshot file for each parallel task)\n")
sys.stderr.write("snapshot_clean=[T | F] (T, if true clean previous iter snapshots\n")
sys.stderr.write("random_initialise_pos=[T | F] (T, if true randomize the initial positions\n")
sys.stderr.write("random_initialise_cell=[T | F] (T, if true randomize the initial cell\n")
sys.stderr.write("LAMMPS_molecular_info=[T | F] (T, if true create lammps molecular info\n")
sys.stderr.write("initial_walk_N_walks=int (0 number of rounds for initial walk) \n")
sys.stderr.write("initial_walk_adjust_interval=int (10 interval (in walks) between adjustments of steps during initial walk) \n")
sys.stderr.write("initial_walk_Emax_offset_per_atom=float (1.0, offset (per atom) to add to Emax for initial walks) \n")
sys.stderr.write("initial_walk_only=[T | F] (F, quit after doing initial walk) \n")
sys.stderr.write("traj_interval=int (100, <=0 for no trajectory)\n")
sys.stderr.write("delta_random_seed=int (-1, < 0 for seed from /dev/urandom)\n")
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)")
[docs]
def excepthook_mpi_abort(exctype, value, tb):
print( print_prefix,'Uncaught Exception Type:', exctype)
print( print_prefix,'Value:', value)
print( print_prefix,'Traceback:', tb)
# print_tb(tb)
print_exception(exctype, value, tb)
print( print_prefix, "Aborting")
sys.stdout.flush()
try:
MPI.COMM_WORLD.Abort(1)
except:
pass
sys.exit(1)
[docs]
def exit_error(message, stat):
sys.stderr.write(message)
sys.stderr.flush()
try:
MPI.COMM_WORLD.Abort(stat)
except:
pass
sys.exit(stat)
[docs]
def str_to_logical(x):
if x in [ "t", "T", "true", "True", "TRUE", "1" ]:
return True
if x in [ "f", "F", "false", "False", "FALSE", "0" ]:
return False
exit_error("Failed to parse logical\n",5)
[docs]
def internal_Vlj(dr_vec, Eshift=0.0):
return (np.sum(dr_vec**-12 - dr_vec**-6) - Eshift*len(dr_vec))
# pos is Nx3, l is a scalar box size
[docs]
def energy_internal_pos(pos, l):
n_at=np.size(pos,0)
dr_list = np.zeros( (n_at*(n_at-1)/2,3) )
offset=0
# calc relative pos vector
for i in range(n_at):
dr_list[offset:offset+(n_at-i-1),0] = pos[i,0] - pos[i+1:,0]
dr_list[offset:offset+(n_at-i-1),1] = pos[i,1] - pos[i+1:,1]
dr_list[offset:offset+(n_at-i-1),2] = pos[i,2] - pos[i+1:,2]
offset += n_at-i-1
# apply min image
dr_list[:,:] -= l*np.floor(dr_list[:,:]/l+0.5)
# magnitude of vectors
dr_list_mag = np.sqrt(np.sum(dr_list*dr_list,axis=1))
return np.sum(internal_Vlj(dr_list_mag[np.where(dr_list_mag < internal_cutoff)],Eshift))
[docs]
def energy_internal(at):
return energy_internal_pos(at.get_positions(), at.get_cell()[0,0])
[docs]
def eval_energy_PE(at):
if do_calc_ASE or do_calc_lammps:
if do_calc_lammps:
#NB only MD can make crazy positions, so maybe just do this after MD propagation?
at.wrap()
energy = at.get_potential_energy()
elif do_calc_internal:
energy = energy_internal(at)
elif do_calc_fortran:
energy = f_MC_MD.eval_energy(at)
else:
sys.stderr.write("No way to eval_energy_PE()\n", 5)
return energy
[docs]
def eval_energy_KE(at):
if at.has('momenta') and at.has('masses'):
energy = at.get_kinetic_energy()
else:
energy = 0.0
return energy
[docs]
def eval_energy_PV(at):
return movement_args['MC_cell_P']*at.get_volume()
[docs]
def eval_energy_mu(at):
if movement_args['semi_grand_potentials'] is not None:
mu = movement_args['semi_grand_potentials']
energy = np.sum([ mu[Z] for Z in at.get_atomic_numbers() ])
else:
energy = 0.0
return energy
[docs]
def eval_energy(at, do_KE=True, do_PE=True):
"""Calls appropriate functions to calculate the potential energy, kinetic energy and the p*V term.
"""
energy = eval_energy_PV(at) + eval_energy_mu(at)
if do_PE:
energy += eval_energy_PE(at)
if do_KE:
energy += eval_energy_KE(at)
return energy
[docs]
def eval_forces(at):
if do_calc_ASE or do_calc_lammps:
if do_calc_lammps:
#NB only MD can make crazy positions, so maybe just do this after MD propagation?
at.wrap()
forces = at.get_forces()
elif do_calc_internal:
exit_error('no forces for do_calc_internal', 10)
elif do_calc_fortran:
forces = np.zeros( at.get_positions().shape )
f_MC_MD.eval_forces(at, forces)
else:
sys.stderr.write("No way to eval_forces()\n", 5)
return forces
[docs]
def propagate_NVE_ASE(at, dt, n_steps):
at.calc = pot
# dt is being converted from ASE units to fs
if ns_args['debug'] >= 10:
if movement_args['keep_atoms_fixed'] > 0: #LIVIA
from ase.constraints import FixAtoms
c = FixAtoms(indices=[atom.index for atom in at if atom.index < movement_args['keep_atoms_fixed']])
at.set_constraint(c)
vv = VelocityVerlet(at, timestep=dt, logfile='-')
else:
if movement_args['keep_atoms_fixed'] > 0: #LIVIA
from ase.constraints import FixAtoms
c = FixAtoms(indices=[atom.index for atom in at if atom.index < movement_args['keep_atoms_fixed']])
at.set_constraint(c)
vv = VelocityVerlet(at, timestep=dt)
vv.run(n_steps)
[docs]
def propagate_lammps(at, dt, n_steps, algo, Emax=None):
if pot.first_propagate:
pot.first_propagate = False
else:
pot.lmp.command('unfix 1 ')
# set appropriate fix
if algo == 'NVE':
if movement_args['keep_atoms_fixed'] > 0: #LIVIA
#pot.lmp.command('group mobile id > ' + str(movement_args['keep_atoms_fixed']))
pot.lmp.command('fix 1 mobile nve')
else:
pot.lmp.command('fix 1 all nve')
elif algo == 'GMC':
# GMC does not work with fixed atoms # Ray Yang
if movement_args['keep_atoms_fixed'] > 0:
exit_error("propagate_lammps got algo '%s', which does not work with fixed atoms\n" % algo)
# Hard coded value needed for LAMMPS. Let's hope our RNG maximum is at
# least this large.
pot.lmp.command('fix 1 all ns/gmc {} {} '.format(rng.int_uniform(1, 900000000), Emax))
else:
exit_error(
"propagate_lammps got algo '%s', neither NVE nor GMC\n" % algo)
# NB maybe not do this every time? Just _after_ MD, since that's the only
# way position can become crazy?
at.wrap()
# do propagation, recover from exception if needed
try:
if algo == 'NVE':
pot.propagate(at, properties=['energy', 'forces'],
system_changes=['positions'], n_steps=n_steps, dt=dt)
else:
pot.propagate(at, properties=['energy', 'forces'],
system_changes=['positions'], n_steps=n_steps, dt=dt,
dt_not_real_time=True)
except Exception as err:
# clean up and return failure
if ns_args['debug'] >= 4:
print("propagate_lammps got exception ", err)
pot.restart_lammps(at)
pot.first_propagate = True
return False
return True
[docs]
def velo_rv_mag(n):
if movement_args['2D']:
#unit_rv[:,2] = 0.0 # commented out from here, as unit_rv is unknown to this function at this point
nDOF = 2.0
else:
nDOF = 3.0
# In 3D rv_mag should have prob distrib p(r) = r^(3N-1).
# Using transformation rule p(y) = p(x) |dx/dy|, with p(y) = y^{3N-1} and p(x) = 1,
# one gets dx/dy = y^{3N-1}
# x = y^{3N}
# y = x^{1/3N}
return rng.float_uniform(0.0, 1.0)**(1.0/(nDOF*n))
[docs]
def velo_unit_rv(n):
unit_rv = rng.normal(1.0, (n, 3))
unit_rv /= np.linalg.norm(unit_rv)
return unit_rv
[docs]
def gen_random_velo(at, KEmax, unit_rv=None):
if unit_rv is None:
unit_rv = velo_unit_rv(len(at))
if movement_args['2D']:
unit_rv[:, 2] = 0.0
# propose right magnitude of randomness consistent with mobile atoms
if movement_args['keep_atoms_fixed'] > 0:
rv_mag = velo_rv_mag(len(at)-movement_args['keep_atoms_fixed'])
else:
rv_mag = velo_rv_mag(len(at))
#rv_mag = velo_rv_mag(len(at))
# from Baldock thesis Eq. 11.10
# p^{**} = r \mathbf{S} \hat{\mathbf{r}}
# and 11.11
# S_{ij} = \delta_{ij} (2 m_i [ E_{lim} - U(q))])^{1/2}
# p_i = r (2 m_i)^{1/2} (E-U)^{1/2} \hat{r}_i
# v_i = p_i / m_i = r (2/m)^{1/2} (E-U)^{1/2} \hat{r}_i
masses = at.get_masses()
# the "np.sqrt(2.0/np.array([masses,]*3).transpose()) * np.sqrt(KEmax)" is filling it up with the same term
velocities = rv_mag * np.sqrt(2.0/np.array([masses,]*3).transpose()) * np.sqrt(KEmax) * unit_rv
# zero velocities of fixed atoms (these were not taken into account in the nDOF)
if movement_args['keep_atoms_fixed'] > 0:
velocities[:movement_args['keep_atoms_fixed']]=0.0
return velocities
[docs]
def pairwise(iterable):
a = iter(iterable)
return zip(a, a)
[docs]
def rotate_dir_3N(vec, max_ang):
if max_ang <= 0.0:
return
# apply random rotations
indices = [ (int(i/3), i%3) for i in range(3*vec.shape[0]) ]
rng.shuffle_in_place(indices)
for ((ind_1_a,ind_1_c), (ind_2_a,ind_2_c)) in pairwise(indices):
ang = rng.float_uniform(-max_ang,max_ang)
c_ang = np.cos(ang)
s_ang = np.sin(ang)
v_1 = vec[ind_1_a,ind_1_c] * c_ang + vec[ind_2_a,ind_2_c] * s_ang
v_2 = -vec[ind_1_a,ind_1_c] * s_ang + vec[ind_2_a,ind_2_c] * c_ang
vec[ind_1_a,ind_1_c] = v_1
vec[ind_2_a,ind_2_c] = v_2
[docs]
def rotate_dir_3N_2D(vec, max_ang):
if max_ang <= 0.0:
return
# apply random rotations
indices = [ (int(i/2), i%2) for i in range(2*vec.shape[0]) ]
rng.shuffle_in_place(indices)
for ((ind_1_a,ind_1_c), (ind_2_a,ind_2_c)) in pairwise(indices):
ang = rng.float_uniform(-max_ang,max_ang)
c_ang = np.cos(ang)
s_ang = np.sin(ang)
v_1 = vec[ind_1_a,ind_1_c] * c_ang + vec[ind_2_a,ind_2_c] * s_ang
v_2 = -vec[ind_1_a,ind_1_c] * s_ang + vec[ind_2_a,ind_2_c] * c_ang
vec[ind_1_a,ind_1_c] = v_1
vec[ind_2_a,ind_2_c] = v_2
[docs]
def rej_free_perturb_velo(at, Emax, KEmax, rotate=True):
#DOC
#DOC ``rej_free_perturb_velo``
if not at.has('momenta') or not at.has('masses'):
return
KEmax_use = None
if Emax is not None:
initial_KE = eval_energy_KE(at)
KEmax_use = Emax - (at.info['ns_energy'] - initial_KE)
if KEmax > 0.0:
if Emax is None:
KEmax_use = KEmax
else:
KEmax_use = min(KEmax_use, KEmax)
if KEmax_use is None:
exit_error("rej_free_perturb_velo() called with Emax and KEmax both None\n", 9)
orig_KE = at.get_kinetic_energy()
#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 **else**: (perturb velocities)
else: # perturb
velo = at.get_velocities()
velo_mag = np.linalg.norm(velo)
#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 **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
# WARNING: check this for variable masses
sqrt_masses_2D = np.sqrt(at.get_masses().reshape((len(at), 1)))
scaled_vel = gen_random_velo(at, KEmax_use, velo/velo_mag) * sqrt_masses_2D
if rotate:
if movement_args['2D']:
rotate_dir_3N_2D(scaled_vel, movement_args['atom_velo_rej_free_perturb_angle'])
else:
if movement_args['keep_atoms_fixed'] > 0: # this is the new line
rotate_dir_3N(scaled_vel[movement_args['keep_atoms_fixed']:, :], movement_args['atom_velo_rej_free_perturb_angle'])
else:
rotate_dir_3N(scaled_vel, movement_args['atom_velo_rej_free_perturb_angle'])
at.set_velocities(scaled_vel / sqrt_masses_2D)
if movement_args['keep_atoms_fixed'] > 0: # for surface simulations
random_velo = at.get_velocities()
random_velo[:movement_args['keep_atoms_fixed'], :] = 0
at.set_velocities(random_velo)
new_KE = at.get_kinetic_energy()
# rej_free_perturb_velo expects at.info['ns_energy'] to be set correctly initially
at.info['ns_energy'] += new_KE-orig_KE
[docs]
def do_MC_atom_velo_walk(at, movement_args, Emax, nD, KEmax):
#DOC
#DOC ``do_MC_atom_velo_walk``
#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 **else**:
#DOC \item do MC pertubation to velocities
n_steps = movement_args['velo_traj_len']
step_size = movement_args['MC_atom_velo_step_size']
n_try = n_steps*len(at)
initial_KE = eval_energy_KE(at)
KEmax_use = Emax - (at.info['ns_energy'] - initial_KE) # expecting ns_energy = KE + PE (+PV)
if KEmax > 0.0 and KEmax < KEmax_use:
KEmax_use = KEmax
if do_calc_fortran:
(n_try, n_accept, final_KE) = f_MC_MD.MC_atom_walk_velo(at, n_steps, step_size, nD, KEmax_use)
at.info['ns_energy'] += final_KE-initial_KE
else:
masses = at.get_masses()
velocities = at.get_velocities()
KE = initial_KE
n_accept = 0
for i_step in range(n_steps):
at_list = list(range(len(at)))
rng.shuffle_in_place(at_list)
d_vel = rng.normal(step_size, (len(at), 3))
if movement_args['keep_atoms_fixed'] > 0: # for surface simulations
d_vel[:movement_args['keep_atoms_fixed'], :] = 0
for i_at in at_list:
d_KE = 0.5*masses[i_at]*(np.sum((velocities[i_at,:]+d_vel[i_at,:])**2) - np.sum(velocities[i_at,:]**2))
if KE + d_KE < KEmax_use:
velocities[i_at,:] += d_vel[i_at,:]
KE += d_KE
n_accept += 1
at.set_velocities(velocities)
at.info['ns_energy'] += KE-initial_KE
return {'MC_atom_velo': (n_try, n_accept)}
[docs]
def do_MD_atom_walk(at, movement_args, Emax, KEmax):
#DOC
#DOC ``do_MD_atom_walk``
""" perform MD walk on the configuration """
orig_E = at.info['ns_energy']
nD = 3
if movement_args['2D']:
nD = 2
if orig_E >= Emax:
print(print_prefix, ": WARNING: orig_E =", orig_E, " >= Emax =", Emax)
#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()
pre_MD_velo = at.get_velocities()
if ns_args['n_extra_data'] > 0:
pre_MD_extra_data = at.arrays['ns_extra_data'].copy()
pre_MD_E = at.info['ns_energy']
#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
timestep = movement_args['MD_atom_timestep']
for i in range(movement_args['atom_traj_len']):
at.set_momenta(at.get_momenta()+forces*0.5*timestep)
at.set_positions(at.get_positions()+at.get_momenta()*timestep)
try:
forces = eval_forces(at)
except:
final_E = 2.0*abs(Emax)
break
at.set_momenta(at.get_momenta()+forces*0.5*timestep)
if final_E is None: # didn't abort due to exception in eval_forces()
final_E = eval_energy(at)
else:
if do_calc_ASE:
propagate_NVE_ASE(at, dt=movement_args['MD_atom_timestep'], n_steps=movement_args['atom_traj_len'])
final_E = eval_energy(at)
elif do_calc_lammps:
if propagate_lammps(at, dt=movement_args['MD_atom_timestep'], n_steps=movement_args['atom_traj_len'], algo='NVE'):
final_E = pot.results['energy'] + eval_energy(at, do_PE=False)
else: # propagate returned success == False
final_E = 2.0*abs(Emax)
## print("error in propagate_lammps NVE, setting final_E = 2*abs(Emax) =" , final_E)
elif do_calc_fortran:
# Fortran MD code does not support fixed atoms # Ray Yang
if movement_args['keep_atoms_fixed'] > 0:
exit_error("do_MD_atom_walk() called with keep_atoms_fixed > 0, but no way to do that with fortran code\n", 3)
final_E = f_MC_MD.MD_atom_NVE_walk(
at, n_steps=movement_args['atom_traj_len'],
timestep=movement_args['MD_atom_timestep'],
debug=ns_args['debug'])
final_E += eval_energy(at, do_PE=False, do_KE=False)
else:
exit_error("Need some non-ASE, non-fortran, non-lammps way of doing MD\n",3)
reject_fuzz = False
final_KE = eval_energy_KE(at)
#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 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)
#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)
config_size = len(at)
min_bond = bonds[config_size]
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
#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)
if movement_args['MD_atom_velo_flip_accept']:
at.set_velocities(pre_MD_velo) # TODO: should we be redundant in zeroing?
else:
at.set_velocities(-pre_MD_velo)
if ns_args['n_extra_data'] > 0:
at.arrays['ns_extra_data'][...] = pre_MD_extra_data
at.info['ns_energy'] = pre_MD_E
n_accept = 0
#DOC **else**: (accepted)
else: # 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
if movement_args['MD_atom_velo_flip_accept']:
# TODO: negative of 0 is 0
at.set_velocities(-at.get_velocities()) # is there a faster way of doing this in ASE? Can you do at.velocities?
at.info['ns_energy'] = final_E
n_accept = 1
#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)
return {'MD_atom': (1, n_accept)}
[docs]
def do_MC_atom_walk(at, movement_args, Emax, KEmax):
#DOC
#DOC ``do_MC_atom_walk``
n_steps = movement_args['atom_traj_len']
step_size = movement_args['MC_atom_step_size']
step_size_velo = movement_args['MC_atom_velo_step_size']
n_try = n_steps*(len(at)-movement_args['keep_atoms_fixed']) # does this do anything? actual n_try is set within fortran code
n_accept=0
n_accept_velo = None
nD=3
if movement_args['2D']:
nD=2
#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)
if movement_args['MC_atom_Galilean']:
if movement_args['GMC_dir_perturb_angle'] < 0.0 or np.linalg.norm(at.arrays['GMC_direction']) == 0.0:
# completely random orientation, magnitude 1
at.arrays['GMC_direction'][:,:] = rng.normal(1.0, (len(at), 3))
if (movement_args['keep_atoms_fixed'] > 0):
at.arrays['GMC_direction'][:movement_args['keep_atoms_fixed'],:] = 0.0
at.arrays['GMC_direction'] /= np.linalg.norm(at.arrays['GMC_direction']) # what is this line doing?
elif movement_args['GMC_dir_perturb_angle'] > 0.0:
# apply random rotations
if movement_args['keep_atoms_fixed'] > 0:
rotate_dir_3N(at.arrays['GMC_direction'][movement_args['keep_atoms_fixed']:,:], movement_args['GMC_dir_perturb_angle'])
else:
rotate_dir_3N(at.arrays['GMC_direction'], movement_args['GMC_dir_perturb_angle'])
#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()
if (ns_args['debug'] >= 10 and movement_args['2D']):
pp=at.get_positions()
if (any(pp[:,2])>0.00001 and nD==2): #LIVIA
exit_error("Not a 2D system anymore\n",3)
if movement_args['MC_atom_velocities']:
(n_try, n_accept, n_accept_velo, final_E) = f_MC_MD.MC_atom_walk(at, n_steps, step_size, Emax-eval_energy(at, do_PE=False, do_KE=False), nD, movement_args['keep_atoms_fixed'], KEmax, step_size_velo)
at.info['ns_energy'] = final_E + eval_energy(at, do_PE=False, do_KE=False)
else:
if movement_args['MC_atom_Galilean']:
(n_try, n_accept, final_E) = f_MC_MD.GMC_atom_walk(at, n_steps, step_size, Emax-eval_energy(at, do_PE=False), no_reverse=movement_args['GMC_no_reverse'], pert_ang=movement_args['GMC_dir_perturb_angle_during'])
else:
# fixed atoms MC works by including the number of atoms to be kept fixed - LIVIA
(n_try, n_accept, final_E) = f_MC_MD.MC_atom_walk(at, n_steps, step_size, Emax-eval_energy(at, do_PE=False), nD, movement_args['keep_atoms_fixed'], movement_args['wall_dist'])
at.info['ns_energy'] = final_E + eval_energy(at, do_PE=False, do_KE=True)
elif (do_calc_lammps and movement_args['MC_atom_Galilean'] and ns_args['LAMMPS_fix_gmc']):
orig_pos = at.get_positions()
if (ns_args['debug'] >= 10 and movement_args['2D']):
if (any(orig_pos[:,2])>0.00001): #LIVIA
exit_error("Not a 2D system anymore\n",3)
orig_energy = at.info['ns_energy']
extra_term = eval_energy(at, do_PE=False, do_KE=False)
if propagate_lammps(at, step_size, n_steps, algo='GMC', Emax=Emax-extra_term ):
energy1 = pot.results['energy'] + extra_term
if energy1 < Emax:
at.info['ns_energy'] = energy1
n_accept = n_steps*len(at)
else:
at.info['ns_energy'] = orig_energy
at.set_positions(orig_pos)
n_accept = 0
else: # got an exception, reject traj
at.info['ns_energy'] = orig_energy
at.set_positions(orig_pos)
n_accept = 0
#DOC **else**: (do python MC)
else:
if movement_args['MC_atom_velocities']:
exit_error("MC_atom_velocities only supported for FORTRAN calculator\n", 8)
dz=0.0
#DOC **if** MC_atom_Galilean:
if movement_args['MC_atom_Galilean']:
#DOC \item go Galilean MC in python
do_no_reverse = movement_args['GMC_no_reverse']
if do_no_reverse:
last_good_pos = at.get_positions()
else:
n_reverse = 0
n_reflect = 0
pos = at.get_positions()
d_pos = step_size*at.arrays['GMC_direction']
for i_MC_step in range(n_steps):
if not do_no_reverse:
last_good_pos = at.get_positions()
last_good_d_pos = d_pos.copy()
# perturb direction if needed
rotate_dir_3N(d_pos, movement_args['GMC_dir_perturb_angle_during'])
# step and evaluate
pos += d_pos
at.set_positions(pos)
E = eval_energy(at)
cur_E_is_correct = True
if E >= Emax: # reflect
Fhat = eval_forces(at)
Fhat /= np.linalg.norm(Fhat)
d_pos -= 2.0*Fhat*np.sum(Fhat*d_pos)
n_reflect += 1
if not do_no_reverse: # do reflection step and check for reverse
# step on reflected traj and evaluate
pos += d_pos
at.set_positions(pos)
E = eval_energy(at)
cur_E_is_correct = True
if E >= Emax: # reverse
pos[:,:] = last_good_pos
at.set_positions(pos)
cur_E_is_correct = False
d_pos[:,:] = -last_good_d_pos
n_reflect -= 1
n_reverse += 1
if do_no_reverse: # accept/reject
n_try = 1
if E < Emax: # accept
# save new E, pos, direction
at.info['ns_energy'] = E
at.set_positions(pos)
at.arrays['GMC_direction'][:,:] = d_pos/np.linalg.norm(d_pos)
n_accept = 1
else: # reject
# E and GMC_direction in at were never overwritten, no need to restore, but do need to reverse dir
at.set_positions(last_good_pos)
at.arrays['GMC_direction'] *= -1.0
n_accept = 0
else:
if n_reverse > 0 and not cur_E_is_correct:
E = eval_energy(at)
at.info['ns_energy'] = E
at.set_positions(pos)
at.arrays['GMC_direction'][:,:] = d_pos/np.linalg.norm(d_pos)
n_try = n_reflect + n_reverse
n_accept = n_reflect
#DOC **else**:
else:
#DOC \item **loop** atom\_traj\_len times
for i_MC_step in range(n_steps):
#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:
#DOC \item propose single atom move
if movement_args['MC_atom_uniform_rv']:
dx = rng.float_uniform(-step_size,step_size)
dy = rng.float_uniform(-step_size,step_size)
dz = 0.0
if not movement_args['2D']:
dz = rng.float_uniform(-step_size,step_size)
else:
dx = rng.normal(step_size)
dy = rng.normal(step_size)
dz = 0.0
if not movement_args['2D']:
dz = rng.normal(step_size)
orig_energy = at.info['ns_energy']
orig_pos = at.get_positions()
if ns_args['n_extra_data'] > 0:
orig_extra_data = at.arrays['ns_extra_data'].copy()
new_pos = orig_pos.copy()
new_pos[i_at,:] += (dx, dy,dz)
at.set_positions(new_pos)
#DOC \item accept/reject on E < Emax
energy = eval_energy(at)
if energy >= Emax:
# reject move
at.set_positions(orig_pos)
if ns_args['n_extra_data'] > 0:
at.arrays['ns_extra_data'][...] = orig_extra_data
energy = orig_energy
else:
n_accept += 1
at.info['ns_energy'] = energy
out = {}
if n_accept_velo is not None:
out['MC_atom_velo'] = (n_try, n_accept_velo)
out['MC_atom'] = (n_try, n_accept)
return out
[docs]
def propose_volume_step(at, step_size, flat_V_prior):
dV = rng.normal(step_size*len(at))
orig_V = at.get_volume()
new_V = orig_V+dV
if new_V < 0: # negative number cannot be raised to fractional power, so this is only to avoid fatal error during the run
new_V=abs(new_V)
# print("Warning, the step_size for volume change might be too big, resulted in negative new volume", step_size, dV, orig_V+dV)
#print("TRANSFORM", new_V, orig_V, dV, step_size, len(at))
transform = np.identity(3)*(new_V/orig_V)**(1.0/3.0)
if flat_V_prior:
p_accept = 1.0
else:
p_accept = min(1.0, (new_V/orig_V)**len(at))
return (p_accept, transform)
[docs]
def propose_area_step(at, step_size, flat_V_prior):
dA = rng.normal(step_size*len(at))
orig_cell = at.get_cell()
orig_A = at.get_volume() / orig_cell[2,2]
new_A = orig_A+dA
if new_A < 0: # negative number cannot be raised to fractional power, so this is only to avoid fatal error during the run
new_A=abs(new_A)
# print("Warning, the step_size for volume change might be too big, resulted in negative new volume", step_size, dV, orig_V+dV)
#print("TRANSFORM", new_V, orig_V, dV, step_size, len(at))
transform = np.identity(3)*np.sqrt(new_A/orig_A)
transform[2,2] = 1.0 # no change in the Z direction
if flat_V_prior:
p_accept = 1.0
else:
p_accept = min(1.0, (new_A/orig_A)**len(at))
return (p_accept, transform)
[docs]
def propose_shear_step(at, step_size):
if movement_args['2D']:
# pick random vector
rnd_vec_ind = rng.int_uniform(0, 2)
# turn other two into orthonormal pair
other_vec_ind = list(range(3))
other_vec_ind.remove(rnd_vec_ind)
orig_cell = at.get_cell()
v1 = orig_cell[other_vec_ind[0],:].copy()
v1 /= np.sqrt(np.dot(v1,v1))
#pick random magnitude
rv1 = rng.normal(step_size)
# create new cell and transformation matrix
new_cell = orig_cell.copy()
new_cell[rnd_vec_ind,:] += rv1*v1
else:
# pick random vector
rnd_vec_ind = rng.int_uniform(0, 3)
# turn other two into orthonormal pair
other_vec_ind = list(range(3))
other_vec_ind.remove(rnd_vec_ind)
orig_cell = at.get_cell()
v1 = orig_cell[other_vec_ind[0],:].copy()
v2 = orig_cell[other_vec_ind[1],:].copy()
v1 /= np.sqrt(np.dot(v1,v1))
v2 -= v1*np.dot(v1,v2)
v2 /= np.sqrt(np.dot(v2,v2))
# pick random magnitudes
rv1 = rng.normal(step_size)
rv2 = rng.normal(step_size)
# create new cell and transformation matrix
new_cell = orig_cell.copy()
new_cell[rnd_vec_ind,:] += rv1*v1 + rv2*v2
transform = np.dot(np.linalg.inv(orig_cell), new_cell)
return (1.0, transform)
[docs]
def propose_stretch_step(at, step_size):
if movement_args['2D']:
rnd_v1_ind = 0
rnd_v2_ind = 1
else:
rnd_v1_ind = rng.int_uniform(0, 3)
rnd_v2_ind = rng.int_uniform(0, 3)
if rnd_v1_ind == rnd_v2_ind:
rnd_v2_ind = (rnd_v2_ind+1) % 3
rv = rng.normal(step_size)
transform = np.identity(3)
transform[rnd_v1_ind,rnd_v1_ind] = np.exp(rv)
transform[rnd_v2_ind,rnd_v2_ind] = np.exp(-rv)
return (1.0, transform)
[docs]
def min_aspect_ratio(vol, cell):
min_aspect_ratio = sys.float_info.max
nD=3
if movement_args['2D']:
# in 2D, do not check the Z direction
nD=2
for i in range(nD):
vi = cell[i,:]
vnorm_hat = np.cross(cell[(i+1)%3,:],cell[(i+2)%3,:])
vnorm_hat /= np.sqrt(np.dot(vnorm_hat,vnorm_hat))
min_aspect_ratio = min(min_aspect_ratio, abs(np.dot(vnorm_hat,vi)))
if movement_args['2D']:
# in 2D divide by the square root of area
min_aspect_ratio /= (vol/np.sqrt(np.dot(cell[2,:],cell[2,:])))**(1.0/2.0)
else:
min_aspect_ratio /= vol**(1.0/3.0)
return min_aspect_ratio
[docs]
def do_cell_step(at, Emax, p_accept, transform):
if p_accept < 1.0 and rng.float_uniform(0.0,1.0) > p_accept:
return False
# save old config and apply transformation
orig_cell = at.get_cell()
new_cell = np.dot(orig_cell,transform)
new_vol = abs(np.dot(new_cell[0,:],np.cross(new_cell[1,:],new_cell[2,:])))
# check size and shape constraints
if new_vol > ns_args['max_volume_per_atom']*len(at) or new_vol < ns_args['min_volume_per_atom']*len(at):
return False
if min_aspect_ratio(new_vol, new_cell) < movement_args['MC_cell_min_aspect_ratio']:
return False
# orig_cell already saved previous value
orig_pos = at.get_positions()
if ns_args['n_extra_data'] > 0:
extra_data = at.arrays['ns_extra_data'].copy()
# set new positions and velocities
at.set_cell(new_cell, scale_atoms=True)
#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)
config_size = len(at)
min_bond = bonds[config_size]
accept_len = min_bond > ns_args['min_nn_dis']
else:
accept_len = True
#VGF end
if Emax is None:
return
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
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:
at.arrays['ns_extra_data'][...] = extra_data
return False
[docs]
def do_cell_shape_walk(at, movement_args):
possible_moves = {
'MC_cell_shear': propose_shear_step,
'MC_cell_stretch': propose_stretch_step }
items = list(possible_moves.items())
for i in range(movement_args['cell_shape_equil_steps']):
rng.shuffle_in_place(items)
for key, propose_step_func in items:
(p_accept, transform) = propose_step_func(at, movement_args[key+"_step_size"])
do_cell_step(at, None, p_accept, transform)
[docs]
def do_MC_semi_grand_step(at, movement_args, Emax, KEmax):
global Z_list
Z = at.get_atomic_numbers()
n_types = len(movement_args['semi_grand_potentials'])
at_i = rng.int_uniform(0,len(at))
type_i = rng.int_uniform(0,n_types)
while Z_list[type_i] == Z[at_i]:
type_i = rng.int_uniform(0,n_types)
Z_new = Z.copy()
Z_new[at_i] = Z_list[type_i]
at.set_atomic_numbers(Z_new)
new_energy = eval_energy(at)
new_KE = eval_energy_KE(at)
if new_energy < Emax and (KEmax < 0.0 or new_KE < KEmax): # accept
# update energy
at.info['ns_energy'] = new_energy
n_accept = 1
else:
# undo
at.set_atomic_numbers(Z)
n_accept = 0
return (1, {'MC_semi_grand' : (1, n_accept) })
[docs]
def do_MC_swap_step(at, movement_args, Emax, KEmax):
#DOC
#DOC ``do_MC_swap_step``
Z = at.get_atomic_numbers()
#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 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 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])):
# print(print_prefix, ": do_MC_swap try to find cluster ", cluster_size)
c1 = pick_interconnected_clump.pick_interconnected(rng, len(at), i_list, j_list, cluster_size, r_cut)
c2 = pick_interconnected_clump.pick_interconnected(rng, len(at), i_list, j_list, cluster_size, r_cut)
# print(print_prefix, ": do_MC_swap got ", c1, c2)
# decrement cluster size
cluster_size -= 1
if cluster_size < 1:
cluster_size = 1
# failed to find appropriate
if c1 is None or c2 is None or np.all(Z[c1] == Z[c2]):
# print(print_prefix, ": do_MC_swap giving up on cluster ",c1,c2)
# return 1 for number of model calls so walk will finish, even if no clusters are ever found
return (1, {})
# print(print_prefix, ": do_MC_swap trying cluster", c1, c2)
p_1_orig = at.positions[c1,:].copy()
p_2_orig = at.positions[c2,:].copy()
at.positions[c1,:] = p_2_orig
at.positions[c2,:] = p_1_orig
if not movement_args['swap_velo']:
# TODO: don't do this for the time being, but we should be able to modify this
# If we just swap particles by swapping positions, velocities will follow to the new position, and so will effectively be swapped
# Therefore, if we _don't_ want to actually swap velocities, we have to (apparently) swap them
velocities = at.get_velocities()
v_1_orig = velocities[c1,:].copy()
v_2_orig = velocities[c2,:].copy()
if movement_args['no_swap_velo_fix_mag_alt']:
# make velocity have the same direction, but magnitude from the previous velocity of the other particle so local KE and momentum are preserved
n1 = np.linalg.norm(v_1_orig)
n2 = np.linalg.norm(v_2_orig)
velocities[c1,:] = (v_2_orig.T * n1/n2).T
velocities[c2,:] = (v_1_orig.T * n2/n1).T
else: # this will be executed (unnecessarily) even if all masses are the same
# make velocity have same direction, and magnitude rescaled to keep local kinetic energy the same
masses = at.get_masses()
velocities[c1,:] = (v_2_orig.T * np.sqrt(masses[c2]/masses[c1])).T
velocities[c2,:] = (v_1_orig.T * np.sqrt(masses[c1]/masses[c2])).T
at.set_velocities(velocities)
if ns_args['n_extra_data'] > 0:
extra_data_1_orig = at.arrays['ns_extra_data'][c1,...].copy()
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
#DOC **if** energy < Emax:
#DOC \item accept swap
new_energy = eval_energy(at)
new_KE = eval_energy_KE(at)
if new_energy < Emax and (KEmax < 0.0 or new_KE < KEmax): # accept
at.info['ns_energy'] = new_energy
accept_n = 1
else: # reject
at.positions[c1,:] = p_1_orig
at.positions[c2,:] = p_2_orig
if not movement_args['swap_velo']:
velocities[c1,:] = v_1_orig
velocities[c2,:] = v_2_orig
at.set_velocities(velocities)
if ns_args['n_extra_data'] > 0:
at.arrays['ns_extra_data'][c1,...] = extra_data_1_orig
at.arrays['ns_extra_data'][c2,...] = extra_data_2_orig
accept_n = 0
return (1, {('MC_swap_%d' % len(c1)) : (1, accept_n) })
[docs]
def do_MC_cell_volume_step(at, movement_args, Emax, KEmax):
#DOC
#DOC ``do_MC_cell_volume_step``
step_rv = rng.float_uniform(0.0, 1.0)
if step_rv > movement_args['MC_cell_volume_per_atom_prob']:
return (0, {})
if movement_args['2D']:
(p_accept, transform) = propose_area_step(at, movement_args['MC_cell_volume_per_atom_step_size'], movement_args['MC_cell_flat_V_prior'])
else:
(p_accept, transform) = propose_volume_step(at, movement_args['MC_cell_volume_per_atom_step_size'], movement_args['MC_cell_flat_V_prior'])
if do_cell_step(at, Emax, p_accept, transform):
return (1, {'MC_cell_volume_per_atom' : (1, 1) })
else:
return (1, {'MC_cell_volume_per_atom' : (1, 0) })
[docs]
def do_MC_cell_shear_step(at, movement_args, Emax, KEmax):
#DOC
#DOC ``do_MC_cell_shear_step``
step_rv = rng.float_uniform(0.0, 1.0)
if step_rv > movement_args['MC_cell_shear_prob']:
return (0, {})
(p_accept, transform) = propose_shear_step(at, movement_args['MC_cell_shear_step_size'])
if do_cell_step(at, Emax, p_accept, transform):
return (1, {'MC_cell_shear' : (1, 1) })
else:
return (1, {'MC_cell_shear' : (1, 0) })
[docs]
def do_MC_cell_stretch_step(at, movement_args, Emax, KEmax):
#DOC
#DOC ``do_MC_cell_stretch_step``
step_rv = rng.float_uniform(0.0, 1.0)
if step_rv > movement_args['MC_cell_stretch_prob']:
return (0, {})
(p_accept, transform) = propose_stretch_step(at, movement_args['MC_cell_stretch_step_size'])
if do_cell_step(at, Emax, p_accept, transform):
return (1, {'MC_cell_stretch' : (1, 1) })
else:
return (1, {'MC_cell_stretch' : (1, 0) })
[docs]
def do_atom_walk(at, movement_args, Emax, KEmax):
#DOC
#DOC ``do_atom_walk``
n_reps = movement_args['n_atom_steps_per_call']
out = {}
#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))
elif movement_args['atom_algorithm'] == 'MD':
accumulate_stats(out, do_MD_atom_walk(at, movement_args, Emax, KEmax))
else:
exit_error("do_atom_walk got unknown 'atom_algorithm' = '%s'\n" % movement_args['atom_algorithm'], 5)
return (n_reps*movement_args['atom_traj_len']*movement_args['atom_traj_len_cost_multiplier'], out)
[docs]
def rand_perturb_energy(energy, perturbation, Emax=None):
if Emax is None:
if abs(energy) <= 1.0:
energy += rng.float_uniform(-1.0, 0.0)*perturbation
else:
energy *= (1.0+rng.float_uniform(-1.0, 0.0)*perturbation)
else:
if abs(energy) <= 1.0:
pert = rng.float_uniform(-1.0, 0.0)*perturbation
n_tries = 0
while energy+pert >= Emax and n_tries < 100:
pert = rng.float_uniform(-1.0, 0.0)*perturbation
n_tries += 1
if energy+pert >= Emax:
print(print_prefix, "WARNING: failed to do random energy perturbation below Emax ", energy, Emax)
energy += pert
else:
pert = 1.0 + rng.float_uniform(-1.0, 0.0)*perturbation
n_tries = 0
while energy*pert >= Emax and n_tries < 100:
pert = 1.0 + rng.float_uniform(-1.0, 0.0)*perturbation
n_tries += 1
if energy*pert >= Emax:
print(print_prefix, "WARNING: failed to do random energy perturbation below Emax ", energy, Emax)
energy *= pert
return energy
####################################################################################################
[docs]
def add_to_list(list, ind, costs, nums, walk_len_avail, first_free_slot):
if nums[0] == 0:
return (walk_len_avail, first_free_slot)
n_steps_r = float(walk_len_avail) * (float(costs[0]*nums[0]) / float(sum(costs*nums))) / float(costs[0])
n_steps_r_fract = n_steps_r - math.floor(n_steps_r)
rv = np.random.rand()
if rv < n_steps_r_fract:
n_steps = int(math.floor(n_steps_r)+1)
else:
n_steps = int(math.floor(n_steps_r))
list[first_free_slot:first_free_slot+n_steps] = [ind]*n_steps
first_free_slot += n_steps
walk_len_avail -= n_steps*costs[0]
return (walk_len_avail, first_free_slot)
[docs]
def create_list(costs, nums, walk_len):
list = []
walk_len_avail = walk_len
first_free_slot = 1
for i in range(len(costs)):
(walk_len_avail, first_free_slot) = add_to_list(list, i, costs[i:], nums[i:], walk_len_avail, first_free_slot)
return list
####################################################################################################
[docs]
def walk_single_walker(at, movement_args, Emax, KEmax):
"""Do random walk on a single atoms object."""
#DOC
#DOC ``walk_single_walker``
out = {}
if movement_args['do_good_load_balance']:
possible_moves = np.array([do_atom_walk,
do_MC_cell_volume_step,
do_MC_cell_shear_step,
do_MC_cell_stretch_step,
do_MC_swap_step,
do_MC_semi_grand_step])
nums = np.array([movement_args['n_atom_steps_n_calls'],
movement_args['n_cell_volume_steps'],
movement_args['n_cell_shear_steps'],
movement_args['n_cell_stretch_steps'],
movement_args['n_swap_steps'],
movement_args['n_semi_grand_steps']])
costs = np.array([movement_args['atom_traj_len']*movement_args['atom_traj_len_cost_multiplier'],
1,
1,
1,
1,
1])
list = create_list(costs, nums, movement_args['n_model_calls'])
for move_i in list:
(t_n_model_calls, t_out) = possible_moves[move_i](at, movement_args,
Emax, KEmax)
accumulate_stats(out, t_out)
else:
possible_moves = ([do_atom_walk] * movement_args['n_atom_steps_n_calls'] +
[do_MC_cell_volume_step] * movement_args['n_cell_volume_steps'] +
[do_MC_cell_shear_step] * movement_args['n_cell_shear_steps'] +
[do_MC_cell_stretch_step] * movement_args['n_cell_stretch_steps'] +
[do_MC_swap_step] * movement_args['n_swap_steps'] +
[do_MC_semi_grand_step] * movement_args['n_semi_grand_steps'])
out = {}
n_model_calls_used = 0
#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 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
(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 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)
# DEBUG print("walk_single_walker end ", eval_energy(at, do_PE=False),
# eval_energy(at) #DEBUG)
return out
[docs]
def max_energy(walkers, n, kinetic_only=False):
"""Collect the current energies of the walkers from all the processes and
chooses the right number of highest energies to be culled"""
# do local max
if kinetic_only:
energies_loc = np.array([eval_energy_KE(at) for at in walkers])
else:
energies_loc = np.array([ at.info['ns_energy'] for at in walkers])
volumes_loc = np.array([ at.get_volume() for at in walkers])
if comm is not None:
energies = np.zeros( (comm.size*len(energies_loc)) )
volumes = np.zeros( (comm.size*len(volumes_loc)) )
# comm.barrier() #BARRIER
comm.Allgather( [ energies_loc, MPI.DOUBLE ], [ energies, MPI.DOUBLE ] )
energies = energies.flatten()
comm.Allgather( [ volumes_loc, MPI.DOUBLE ], [ volumes, MPI.DOUBLE ] )
volumes = volumes.flatten()
else:
energies = energies_loc
volumes = volumes_loc
# n is n_cull
Emax_ind = energies.argsort()[-1:-n-1:-1]
Emax = energies[Emax_ind]
Vmax = volumes[Emax_ind]
# WARNING: assumes that each node has equal number of walkers
rank_of_max = np.floor(Emax_ind/len(walkers)).astype(int)
ind_of_max = np.mod(Emax_ind,len(walkers))
return (Emax, Vmax, rank_of_max, ind_of_max)
# WARNING: full_auto_set_stepsizes shouldn't really depend on walk_stats from a real walk, since it does all its own pilot walks
# right now it uses this data structure to figure out what kind of steps actually occur, but it should probably get this
# information right from movement_args, rather than relying on what happened to have happened in the real walks
[docs]
def full_auto_set_stepsizes(walkers, walk_stats, movement_args, comm, Emax, KEmax, size_n_proc):
"""Automatically set all step sizes. Returns the time (in seconds) taken for the routine to run."""
#DOC
#DOC ``full_auto_set_stepsizes``
#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
print_prefix = "full_auto"
if comm is not None:
print_prefix = "%d %s" % (comm.rank, print_prefix)
full_auto_start_time = time.time()
n_samples_per_move_type=200 # total number of (H)MC moves used to set each step length
if comm is not None:
walk_n_walkers = int(np.ceil(float(n_samples_per_move_type)/size_n_proc))
# in each trial we will evolve walk_n_walkers configurations
else:
walk_n_walkers = n_samples_per_move_type
#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)
if (comm is not None):
# Identify move types that are being used
# Not all processes may have used the same move types
key_ints=[]
for key in key_list:
if (key == "MD_atom"):
i = 0
elif (key == "MC_atom"):
i = 1
elif (key == "MC_atom_velo"):
i = 2
elif (key == "MC_cell_shear"):
i = 3
elif (key == "MC_cell_stretch"):
i = 4
elif (key == "MC_cell_volume_per_atom"):
i = 5
elif (key[:7] == "MC_swap"):
i = 6
else:
i = -1
key_ints.append(i)
key_flags = 1*np.asarray([ i in key_ints for i in range(7)])
totalkeys = np.zeros( (len(key_flags)), dtype=int)
comm.Allreduce([key_flags, MPI.INT], [totalkeys, MPI.INT], MPI.SUM)
key_list = []
for i in range(7):
if (totalkeys[i]>0):
if (i==0):
key_list.append("MD_atom")
if (i==1):
key_list.append("MC_atom")
if (i==2):
key_list.append("MC_atom_velo")
if (i==3):
key_list.append("MC_cell_shear")
if (i==4):
key_list.append("MC_cell_stretch")
if (i==5):
key_list.append("MC_cell_volume_per_atom")
if (i==6):
key_list.append("MC_swap_")
#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
# 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
# MC_atom # MC atom sweeps
# MC_atom_velo # MC velocity sweeps
# MC_cell_shear # cell shear move
# MC_cell_stretch # cell stretch move
# MC_cell_volume_per_atom # volume move
exploration_movement_args = deepcopy(movement_args)
# turn all step types off to begin with
exploration_movement_args['n_atom_steps_n_calls'] = 0
exploration_movement_args['n_cell_volume_steps'] = 0
exploration_movement_args['n_cell_shear_steps'] = 0
exploration_movement_args['n_cell_stretch_steps'] = 0
exploration_movement_args['n_swap_steps'] = 0
exploration_movement_args['n_semi_grand_steps'] = 0
exploration_movement_args['MC_atom_velocities']=False
nD=3
if exploration_movement_args['2D']:
nD=2
# check that the total number of attempts for this key is not zero
if key in walk_stats:
(n_try, n_accept) = walk_stats[key]
else:
n_try=0
n_try_g = np.zeros( (1), dtype=int)
if (comm is not None):
n_try_s = np.array( [n_try], dtype = int)
comm.Allreduce([n_try_s, MPI.INT], [n_try_g, MPI.INT], MPI.SUM)
else:
n_try_g[0] = n_try
if (n_try_g[0]==0):
continue # skip this key - it is not used
if (key=="MD_atom" or key=="MC_atom"):
exploration_movement_args['n_atom_steps_n_calls'] = 1
# one call to do_atom_walk per walk_single_walker call
exploration_movement_args['n_atom_steps_per_call'] = 1
# do_atom_walk makes one do_MC_walk/do_MD_walk per call
if key == "MC_atom" and not movement_args['MC_atom_Galilean']:
exploration_movement_args['atom_traj_len'] = 1
# 1 atom sweep per do_MC_walk call
exploration_movement_args['n_model_calls'] = 1 * exploration_movement_args['atom_traj_len_cost_multiplier']
# The do_atom_walk routine reports that it has performed
# #model_calls = the number of complete MC sweeps performed,
# rather than single point evaluations.
else:
exploration_movement_args['n_model_calls'] = exploration_movement_args['atom_traj_len']*exploration_movement_args['atom_traj_len_cost_multiplier']
# The do_atom_walk routine reports that it has performed
# #model_calls = number of single point evaluations (time steps)
elif (key=="MC_atom_velo"):
exploration_movement_args['velo_traj_len']=1
exploration_movement_args['MC_atom_velo_walk_rej_free']=False
# one velocities sweep
elif (key=="MC_cell_shear"):
exploration_movement_args['n_cell_shear_steps'] = 1
exploration_movement_args['n_model_calls'] = 1
# one call to do_MC_cell_shear_step per walk_single_walker call
elif (key=="MC_cell_stretch"):
exploration_movement_args['n_cell_stretch_steps'] = 1
exploration_movement_args['n_model_calls'] = 1
# one call to do_MC_cell_stretch_step per walk_single_walker call
elif (key=="MC_cell_volume_per_atom"):
exploration_movement_args['n_cell_volume_steps'] = 1
exploration_movement_args['n_model_calls'] = 1
# one call to do_MC_cell_volume_step per walk_single_walker call
elif (key[:7] == "MC_swap"):
# skip swap moves, since they have no step size
break
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
if key.find("MC") == 0:
if key.find("MC_atom_step_size") and movement_args['MC_atom_Galilean']:
min_rate = movement_args['GMC_adjust_min_rate']
max_rate = movement_args['GMC_adjust_max_rate']
else:
min_rate = movement_args['MC_adjust_min_rate']
max_rate = movement_args['MC_adjust_max_rate']
suffix="step_size"
elif key.find("MD") == 0:
min_rate = movement_args['MD_adjust_min_rate']
max_rate = movement_args['MD_adjust_max_rate']
suffix="timestep"
else:
exit_error("full_auto_set_stepsizes got key '%s', neither MC nor MD\n" % key, 5)
first_walker = rng.int_uniform(0, len(walkers)) # random starting point for config selection
first_time = True # we will make at least two tries. Logical flag ensures this.
steplength_store = movement_args[key+"_"+suffix]
# protects against possible future bugs that would be hard to detect
#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 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
# 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 \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
(n_try, n_accept) = stats_cumul[key]
if comm is not None:
n_try_s = np.array( [n_try], dtype = int)
n_accept_s = np.array( [n_accept], dtype = int)
n_try_g = np.zeros( (1), dtype=int)
n_accept_g = np.zeros( (1), dtype=int)
# comm.barrier() #BARRIER
comm.Allreduce([n_try_s, MPI.INT], [n_try_g, MPI.INT], MPI.SUM)
comm.Allreduce([n_accept_s, MPI.INT], [n_accept_g, MPI.INT], MPI.SUM)
n_try = n_try_g[0]
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
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 rate<max_rate):
#DOC \item **if** the total acceptance rate is within the desired range:
#DOC \item **return** this stepsize
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**:
else:
if( not first_time ): # dodge this the first time round
# Check whether rate and rate_store are on different sides
# of interval
if ((min(rate,rate_store) < min_rate) and (max(rate,rate_store)>max_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:
#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
target = 0.5*(min_rate+max_rate)
if (abs(rate-target)<abs(rate_store-target)):
# take current step length
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
else:
# take saved step length
movement_args[key+"_"+suffix] = steplength_store
exploration_movement_args[key+"_"+suffix] = steplength_store
rate = rate_store
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
#DOC \item Store step length and acceptance rate
# save current step length and acceptance rate
steplength_store = movement_args[key+"_"+suffix]
rate_store = rate
#DOC \item update step length, by :math:`*` or :math:`/` by MC\_adjust\_step\_factor, to improve acceptance rate
#update step length
dir = None
if rate < min_rate:
exp = -1
dir = "down"
elif rate >= max_rate:
exp = 1
dir = "up"
else:
exp = None
# try to adjust
if dir is not None:
movement_args[key+"_"+suffix] *= movement_args['MC_adjust_step_factor']**exp
exploration_movement_args[key+"_"+suffix] *= exploration_movement_args['MC_adjust_step_factor']**exp
if (comm is None or comm.rank == 0) and (ns_args['debug'] >= 1):
print( print_prefix, "new trial step size for %s = %e" % (key, movement_args[key+"_"+suffix]))
#DOC \item Check that step size is not larger than max allowed value (specified by user), and also that step size is not smaller than 10\ :sup:`-20`\ (useful for detecting errors).
# if exceeded maximum, cap change
if movement_args[key+"_"+suffix] > movement_args[key+"_"+suffix+"_max"]:
movement_args[key+"_"+suffix] = movement_args[key+"_"+suffix+"_max"]
exploration_movement_args[key+"_"+suffix] = exploration_movement_args[key+"_"+suffix+"_max"]
# Error check:
if (movement_args[key+"_"+suffix] < 1.0e-20):
exit_error("full_auto_set_stepsizes got '%s'_'%s' = '%e': too small. Is everything correct?\n" % (key, suffix, movement_args[key+"_"+suffix]), 25)
# if was already at max and didn't really change, break
if (movement_args[key+"_"+suffix] == steplength_store):
dir = None
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 **return** step sizes and time taken for routine to run
print_prefix = orig_prefix
full_auto_end_time = time.time()
duration = full_auto_end_time - full_auto_start_time
return duration
[docs]
def adjust_step_sizes(walk_stats, movement_args, comm, do_print_rate=True, monitor_only=False):
"""
Adjust step size to keep the acceptance ratio at the desired level.
"""
for key in walk_stats:
(n_try, n_accept) = walk_stats[key]
if comm is not None:
n_try_s = np.array( [n_try], dtype = int)
n_accept_s = np.array( [n_accept], dtype = int)
n_try_g = np.zeros( (1), dtype=int)
n_accept_g = np.zeros( (1), dtype=int)
# comm.barrier() #BARRIER
comm.Allreduce([n_try_s, MPI.INT], [n_try_g, MPI.INT], MPI.SUM)
comm.Allreduce([n_accept_s, MPI.INT], [n_accept_g, MPI.INT], MPI.SUM)
n_try = n_try_g[0]
n_accept = n_accept_g[0]
if n_try > 0:
rate = float(n_accept)/float(n_try)
if do_print_rate and comm is None or comm.rank == 0:
print( print_prefix, "accept rate for %s = %f (%d)" % (key, rate, n_try))
if monitor_only:
continue
if key.find("MC") == 0:
if key.find("MC_atom_step_size") and movement_args['MC_atom_Galilean']:
min_rate = movement_args['GMC_adjust_min_rate']
max_rate = movement_args['GMC_adjust_max_rate']
else:
min_rate = movement_args['MC_adjust_min_rate']
max_rate = movement_args['MC_adjust_max_rate']
suffix="step_size"
elif key.find("MD") == 0:
min_rate = movement_args['MD_adjust_min_rate']
max_rate = movement_args['MD_adjust_max_rate']
suffix="timestep"
else:
if (comm is None or comm.rank == 0):
print( "WARNING: adjust_step_size got key '%s', neither MC nor MD\n" % key)
continue
if key+"_"+suffix not in movement_args:
continue
dir = None
exp = 0.0
if rate < min_rate:
exp = -1
dir = "down"
elif rate >= max_rate:
exp = 1
dir = "up"
orig_value = movement_args[key+"_"+suffix]
if dir is not None:
# try to adjust
movement_args[key+"_"+suffix] *= movement_args['MC_adjust_step_factor']**exp
# if exceeded maximum, cap change
if movement_args[key+"_"+suffix] > movement_args[key+"_"+suffix+"_max"]:
movement_args[key+"_"+suffix] = movement_args[key+"_"+suffix+"_max"]
# if was already at max and didn't really change, unset dir
if movement_args[key+"_"+suffix] == orig_value:
dir=None
if dir is not None and (comm is None or comm.rank == 0):
print( print_prefix, "adjust_step_sizes adjusted %s %s to %f" % (key, dir, movement_args[key+"_"+suffix]))
[docs]
def zero_stats(d, movement_args):
for key in movement_args:
m = re.search('(.*)_(step_size|timestep)$', key)
if m is not None:
d[m.group(1)] = (0, 0)
[docs]
def accumulate_stats(d_cumul, d):
for key in d:
if key in d_cumul:
d_cumul[key] = tuple([i1+i2 for i1, i2 in zip(d[key], d_cumul[key])])
else:
d_cumul[key] = d[key]
# figure out n_steps to walk on each iteration to get correct expected number
[docs]
def set_n_from_expected(prop):
if movement_args[prop+'_expected'] > 0:
if movement_args[prop] > 0:
exit_error("Got both "+prop+" and "+prop+"_expected, conflict\n", 5)
if rank == 0:
print( "Calculating %s from %s_expected=%d" % (prop, prop, movement_args[prop+'_expected']))
if max_n_cull_per_task*size == n_cull and n_extra_walk_per_task == 0: # no extra walkers
movement_args[prop] = movement_args[prop+'_expected']
if rank == 0:
print( "No extra walkers (n_cull mod n_tasks == 0), trivial, so average n_walks at kill is 1, and %s=%s_expected" % (prop, prop))
else:
# f_c = n_c/n_t [ fraction of total that are culled (and walked once)]
# f_e = n_e/(n_t-n_c) [fraction of ones that aren't culled that are also walked ]
# n(1) <- f_c n_t + (1-f_c)*(1-f_e) n(1)
# n(l >= 2) <- (1-f_c)(1-f_e) n(l) + (1-f_c) f_e n(l-1)
# f = (1-f_c)f_e / (f_c+f_e - f_c f_e)
# n(1)/n_t = f_c/(1-(1-f_c)(1-f_e)) = f_c/(f_c+f_e-f_c f_e)
# n(l >= 2)/n_t = f n(l-1) = f^(l-1) n(1)/n_t
# n_walks = n(1)/n_t + sum_{l=2}^\infty l n(l)/n_t
# = n(1)/n_t (1 + sum_{l=2}^\infty l f^{l-1})
# = n(1)/n_t (1 + \sum{l=1}^\infty (l+1) f^l)
# = n(1)/n_t (1 + \sum f^l + sum l f^l)
# = n(1)/n_t (1 + f/(1-f) + f/(1-f)**2
f_cull = float(n_cull)/float(ns_args['n_walkers'])
f_extra = float( (max_n_cull_per_task*size-n_cull) + n_extra_walk_per_task*size ) / float(ns_args['n_walkers']-n_cull)
f = (1.0-f_cull)*f_extra/(f_cull+f_extra-f_cull*f_extra)
n_walks = f_cull/(f_cull+f_extra-f_cull*f_extra) * (1.0 + f/(1.0-f) + f/(1.0-f)**2)
if rank == 0:
print("f_cull" ,f_cull,"f_extra", f_extra, "f", f)
print( "Calculated average n_walks at kill = ",n_walks, " from f_cull ",f_cull," f_extra (due to otherwise idle processors doing walks and explicitly requested extra walks) ",f_extra)
if n_walks < 1:
print( "WARNING: got average n_walks < 1, but will always do at least 1 walk, so effective %s_expected will be higher than requested" % prop)
print( "Setting %s = ceiling(%s_expected/n_walks)" % (prop, prop))
movement_args[prop] = int(math.ceil(movement_args[prop+'_expected']/n_walks))
else:
if movement_args[prop] > 0 and rank == 0:
print( "WARNING: using absolute number of "+prop)
if rank == 0:
print( "Final value of %s=%d" % (prop, movement_args[prop]))
[docs]
def additive_init_config(at, Emax):
if do_calc_lammps:
exit_error("python additive_init_config doesn't work with LAMMPS, since it varies list of atoms\n", 10)
pos = at.get_positions()
for i_at in range(1,len(at)):
at_new = at[0:i_at+1]
if do_calc_ASE or do_calc_lammps:
at_new.calc = at.calc
at_new.set_positions(pos[0:i_at+1,:])
success = False
for i_try in range(10):
pos[i_at,:] = np.dot(at_new.get_cell(), rng.float_uniform(0.0, 1.0, (3) ))
at_new.set_positions(pos[0:i_at+1,:])
if movement_args['2D']: # zero the Z coordiates in a 2D simulation
temp_at=at_new.get_positions()
temp_at[:,2]=0.0
at_new.set_positions(temp_at)
energy = eval_energy(at_new)
if energy < Emax:
success = True
break
if not success:
exit_error("Failed 10 times to insert atom %d with Emax %f" % (i_at, Emax), 7)
at.set_positions(pos)
return energy
[docs]
def save_snapshot(snapshot_id):
"""
Save the current walker configurations' as a snapshot in the file ``out_file_prefix.iter.rank.config_file_format``
"""
if comm is not None:
comm.barrier() # if parallel, ensure that we are always in sync, so snapshots are always a consistent set
if ns_args['snapshot_per_parallel_task']:
rank_id = "%d" % rank
else:
rank_id = "ALL"
if ns_args['snapshot_per_parallel_task'] or rank == 0:
try:
snapshot_io = open(ns_args['out_file_prefix']+'snapshot.%s.%s.%s' % (snapshot_id, rank_id, ns_args['config_file_format']), "w")
except:
snapshot_io = open(ns_args['out_file_prefix']+'snapshot.%d.%s.%s' % (snapshot_id, rank_id, ns_args['config_file_format']), "w")
root_walkers_write_t0 = time.time()
for at in walkers:
at.info['volume'] = at.get_volume()
at.info['iter']=snapshot_id
ase.io.write(snapshot_io, at, parallel=False, format=ns_args['config_file_format'])
print( "root walkers write time ", time.time() - root_walkers_write_t0)
if not ns_args['snapshot_per_parallel_task']:
if comm is not None: # gather other walkers to do I/O on master task
if rank == 0: # I/O on master task
for r in range(1,size):
remote_walkers_recv_t0 = time.time()
remote_walkers = comm.recv(source=r, tag=2)
print( "save_snapshot remote walkers recv time ", r, time.time() - remote_walkers_recv_t0)
remote_walkers_write_t0 = time.time()
for at in remote_walkers:
at.info['volume'] = at.get_volume()
at.info['iter']=snapshot_id
ase.io.write(snapshot_io, at, format=ns_args['config_file_format'])
print( "save_snapshot remote walkers write time ", r, time.time() - remote_walkers_write_t0)
else: # not master task
comm.send(walkers, dest=0, tag=2)
if ns_args['snapshot_per_parallel_task'] or rank == 0:
snapshot_io.close()
[docs]
def clean_prev_snapshot(iter):
if iter is not None and ns_args['snapshot_clean']:
snapshot_file=ns_args['out_file_prefix']+'snapshot.%d.%d.%s' % (iter, rank, ns_args['config_file_format'])
try:
os.remove(snapshot_file)
except:
print( print_prefix, ": WARNING: Failed to delete '%s'" % snapshot_file)
[docs]
def do_ns_loop():
"""
This is the main nested sampling loop, doing the iterations.
"""
global print_prefix
global cur_config_ind
# cant keep atoms fixed and change the simulation cell at the moment
if (movement_args['keep_atoms_fixed'] > 0
and movement_args['n_cell_volume_steps'] > 0):
exit_error("cant keep atoms fixed and change the simulation cell\n", 11)
nD = 3 # dimensionality of the system
if movement_args['2D']: # simulate a 2D system only
nD = 2
if rank == 0:
if energy_io.tell() == 0: # tell returns the current stream position
if movement_args['do_velocities']:
nExtraDOF = 0
else:
if movement_args['keep_atoms_fixed'] > 0:
nExtraDOF = (n_atoms-movement_args['keep_atoms_fixed'])*nD
else:
nExtraDOF = n_atoms*nD
energy_io.write("%d %d %d %s %d\n" % (ns_args['n_walkers'], ns_args['n_cull'], nExtraDOF, movement_args['MC_cell_flat_V_prior'], n_atoms))
## print(print_prefix, ": random state ", np.random.get_state())
## if rank == 0:
## print(print_prefix, ": common random state ", common_random_state)
if ns_args['debug'] >= 10 and size <= 1:
for at in walkers:
at.info['n_walks'] = 0
for at in walkers:
at.info['KEmax'] = KEmax
if movement_args['MC_cell_P'] > 0:
print(rank, ": initial enthalpy ", at.info['ns_energy'], " PE ", eval_energy_PE(at), " KE ", eval_energy_KE(at), " PV ", eval_energy_PV(at), " mu ", eval_energy_mu(at), " vol ", at.get_volume())
else:
print(rank, ": initial enthalpy ", at.info['ns_energy'], " PE ", eval_energy_PE(at), " KE ", eval_energy_KE(at), " mu ", eval_energy_mu(at), " vol ",at.get_volume())
sys.stdout.flush()
# stats for purpose of adjusting step size
walk_stats_adjust = {}
# stats for purpose of monitoring acceptance rates
walk_stats_monitor = {}
zero_stats(walk_stats_adjust, movement_args)
zero_stats(walk_stats_monitor, movement_args)
initial_time = time.time()
prev_time = initial_time
step_size_setting_duration = 0.0
total_step_size_setting_duration = 0.0
Emax_of_step = None
Emax_save = []
i_ns_step_save = []
traj_walker_list = []
E_dump_list = []
E_dump_list_times = []
verbose = False
# to avoid errors of unassigned values, if in case of a restart the final number of iter is the same as the starting, stop.
if start_first_iter == ns_args['n_iter']:
print("WARNING: Increase the n_iter_times_fraction_killed variable in the input if you want NS cycles to be performed.")
exit_error("starting iteration and the total number of required iterations are the same,hence no NS cycles will be performed\n", 11)
last_log_X_n = 0.0
i_range_mod_n_cull = np.array(range(ns_args['n_cull']))
i_range_plus_1_mod_n_cull = np.mod(np.array(range(ns_args['n_cull']))+1, ns_args['n_cull'])
log_X_n_term = np.log(ns_args['n_walkers']-i_range_mod_n_cull) - np.log(ns_args['n_walkers']+1-i_range_mod_n_cull)
log_X_n_term_cumsum = np.cumsum(log_X_n_term)
log_X_n_term_cumsum_modified = log_X_n_term_cumsum - np.log(ns_args['n_walkers']+1-i_range_plus_1_mod_n_cull)
log_X_n_term_sum = log_X_n_term_cumsum[-1]
if ns_args['converge_down_to_T'] > 0:
converge_down_to_beta = 1.0/(ns_args['kB']*ns_args['converge_down_to_T'])
log_Z_term_max = np.NINF
prev_snapshot_iter = None
pprev_snapshot_iter = None
last_snapshot_time = time.time()
# for estimating current temperature from d log Omega / d E
if ns_args['T_estimate_finite_diff_lag'] > 0:
log_alpha = np.log(float(ns_args['n_walkers']+1-ns_args['n_cull'])/float(ns_args['n_walkers']+1))
Emax_history = collections.deque(maxlen=ns_args['T_estimate_finite_diff_lag'])
if ns_analyzers is not None:
for (ns_analyzer, ns_analyzer_interval) in ns_analyzers:
ns_analyzer.analyze(walkers, -1, "NS_loop start")
# START MAIN LOOP
i_ns_step = start_first_iter
while ns_args['n_iter'] < 0 or i_ns_step < ns_args['n_iter']:
if ns_args['debug'] == -5:
print(i_ns_step, rank, " ".join(["{:.2f}".format(eval_energy(x))
for x in walkers]))
check_memory.check_memory("start_ns_main_loop")
print_prefix = "%d NS %d" % (rank, i_ns_step)
if ns_args['debug'] >= 4 and ns_args['track_configs']:
for at in walkers:
print(print_prefix, "INFO: 10 config_ind ", at.info['config_ind'], " from ", at.info['from_config_ind'], " at ", at.info['config_ind_time'])
if movement_args['adjust_step_interval'] < 0:
zero_stats(walk_stats_adjust, movement_args)
if movement_args['monitor_step_interval'] < 0:
zero_stats(walk_stats_monitor, movement_args)
if ns_args['debug'] >= 20:
print(print_prefix, "%30s" % ": LOOP_TE START 00 ", i_ns_step, ["%.10f" % eval_energy(at) for at in walkers])
print(print_prefix, "%30s" % ": LOOP_PE START 01 ", i_ns_step, ["%.10f" % eval_energy(at, do_KE=False) for at in walkers])
print(print_prefix, "%30s" % ": LOOP_X START 02 ", i_ns_step, ["%.10f" % at.positions[0, 0] for at in walkers])
# get list of highest energy configs
(Emax, Vmax, cull_rank, cull_ind) = max_energy(walkers, n_cull)
Emax_next = Emax[-1]
if rank == 0 and Emax_of_step is not None and Emax[0] > Emax_of_step:
print(print_prefix, ": WARNING: energy above Emax ", Emax_of_step, " bad energies: ", Emax[np.where(Emax > Emax_of_step)], cull_rank[np.where(Emax > Emax_of_step)], cull_ind[np.where(Emax > Emax_of_step)])
# comm.barrier()
# exit_error("Energy above Emax\n", 5)
if rank == 0 and (i_ns_step > start_first_iter and Emax_next >= Emax_of_step):
print("WARNING: Emax not decreasing ", Emax_of_step, Emax_next)
Emax_of_step = Emax_next
if ns_args['min_Emax'] is not None and Emax_of_step < ns_args['min_Emax']:
if rank == 0:
# if the termination was set by a minimum energy, and it is reached, stop.
print("Leaving loop because Emax=", Emax_of_step, " < min_Emax =", ns_args['min_Emax'])
i_ns_step += 1 # add one so outside loop when one is subtracted to get real last iteration it's still correct
break
if rank == 0:
cur_time = time.time()
output_this_iter = (cur_time > prev_time+60 or i_ns_step == 0 or i_ns_step == ns_args['n_iter'] or (ns_args['n_iter'] > 0 and i_ns_step % max(int(ns_args['n_iter']/1000), 1) == 0))
else:
output_this_iter = False
if ns_args['converge_down_to_T'] > 0:
# see ns_analyse.py calc_log_a() for math
log_a = log_X_n_term_sum*i_ns_step + log_X_n_term_cumsum_modified
# DEBUG if rank == 0:
# DEBUG for ii in range(len(log_a)):
# DEBUG print i_ns_step, "log_a, beta, Es, beta*Es ", log_a[ii], beta, Emax[ii], beta*Emax[ii]
log_Z_term_max = max(log_Z_term_max, np.amax(log_a - converge_down_to_beta * Emax))
log_Z_term_last = log_a[-1]-converge_down_to_beta*Emax[-1]
if output_this_iter:
print("log_Z_term max ", log_Z_term_max, "last ", log_Z_term_last, "diff ", log_Z_term_max-log_Z_term_last)
if log_Z_term_last < log_Z_term_max - 10.0:
if rank == 0:
print(print_prefix, "Leaving loop because Z(%f) is converged" % ns_args['converge_down_to_T'])
i_ns_step += 1 # add one so outside loop when one is subtracted to get real last iteration it's still correct
break
if ns_args['T_estimate_finite_diff_lag'] > 0:
Emax_history.append(Emax_of_step)
if output_this_iter:
if ns_args['T_estimate_finite_diff_lag'] > 0 and len(Emax_history) > 1:
beta_estimate = (len(Emax_history)-1)*log_alpha/(Emax_history[-1]-Emax_history[0])
T_estimate = 1.0/(ns_args['kB']*beta_estimate)
else:
T_estimate = -1
print(i_ns_step, "Emax_of_step ", Emax_of_step, "T_estimate ", T_estimate, " loop time ", cur_time-prev_time-step_size_setting_duration, " time spent setting step sizes: ", step_size_setting_duration)
sys.stdout.flush()
prev_time = cur_time
step_size_setting_duration = 0.0
entries_for_this_rank = np.where(cull_rank == rank)[0]
cull_list = cull_ind[entries_for_this_rank]
if rank == 0 and ns_args['debug'] >= 4 and len(cull_ind[entries_for_this_rank]) > 0:
print(print_prefix, "INFO: 20 cull ", cull_ind[entries_for_this_rank], " on ", rank)
# record Emax walkers energies
if rank == 0:
for (E, V) in zip(Emax, Vmax):
energy_io.write("%d %.60f %.60f\n" % (i_ns_step, E, V))
energy_io.flush()
## Save the energies and corresponding iteration numbers in a list then print(them out only when printing a snapshot)
#Emax_save.extend(Emax)
#i_ns_step_save.extend(n_cull*[i_ns_step])
## if it is time to print((i.e. at the same iteration when a snapshot is written, or at every iter if no snapshots - for smooth restarts))
#if ns_args['snapshot_interval'] < 0 or i_ns_step % ns_args['snapshot_interval'] == ns_args['snapshot_interval']-1:
# for istep,E in zip(i_ns_step_save,Emax_save):
# energy_io.write("%d %.60f\n" % (istep, E))
# energy_io.flush()
# #empty the save lists, so they are ready for the next bunch of saved energies
# Emax_save[:]=[]
# i_ns_step_save[:]=[]
# record Emax walkers configurations
if cull_list is not None:
for (i, global_n_offset) in zip(cull_list, entries_for_this_rank):
if ns_args['debug'] >= 10 and size <= 1:
print(print_prefix, "walker killed at age ", walkers[i].info['n_walks'])
# store culled config in list to be written (when
# snapshot_interval has passed) every traj_interval steps
global_n = i_ns_step*n_cull + global_n_offset
if ns_args['traj_interval'] > 0 and global_n % ns_args['traj_interval'] == ns_args['traj_interval']-1:
walker_copy = walkers[i].copy()
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()
traj_walker_list.append(walker_copy)
# if tracking all configs, save this one that has been culled
if track_traj_io is not None:
at = walkers[i].copy()
at.info['culled'] = True
ase.io.write(track_traj_io, at, parallel=False, format=ns_args['config_file_format'])
if ns_args['E_dump_interval'] > 0 and i_ns_step % ns_args['E_dump_interval'] == 0: # ns_args['E_dump_interval']-1:
if walkers[0].has('masses') and walkers[0].has('momenta'):
E_dump_list.append([w.info['ns_energy'] - w.get_kinetic_energy() for w in walkers])
else:
E_dump_list.append([w.info['ns_energy'] for w in walkers])
E_dump_list_times.append(i_ns_step)
if ns_args['traj_interval'] > 0:
for at in traj_walker_list:
ase.io.write(traj_io, at, parallel=False, format=ns_args['config_file_format'])
traj_io.flush()
traj_walker_list=[]
# print(the recorded Emax walkers configurations to output file)
if (ns_args['snapshot_interval'] < 0 or i_ns_step % ns_args['snapshot_interval'] == ns_args['snapshot_interval']-1 or
(ns_args['snapshot_seq_pairs'] and i_ns_step > 0 and i_ns_step%ns_args['snapshot_interval'] == 0) ) :
##NB if ns_args['traj_interval'] > 0:
##NB for at in traj_walker_list:
##NB ase.io.write(traj_io, at, parallel=False, format=ns_args['config_file_format'])
##NB traj_io.flush()
##NB traj_walker_list=[]
if ns_args['E_dump_interval'] > 0:
if comm is not None:
E_dump_list_all = np.array(comm.allgather(E_dump_list))
else:
E_dump_list_all = np.array(E_dump_list)
if rank == 0:
for i in range(E_dump_list_all.shape[1]):
E_dump_io.write("step %d\n" % E_dump_list_times[i])
if len(E_dump_list_all.shape) == 3:
np.savetxt(E_dump_io, E_dump_list_all[:,i,:])
else:
np.savetxt(E_dump_io, E_dump_list_all[i,:])
E_dump_io.flush()
E_dump_list = []
E_dump_list_all = None
E_dump_list_times = []
# calculate how many will be culled on each rank
n_cull_of_rank = np.array([sum(cull_rank == r) for r in range(size)])
# label configs to be culled
status = np.empty((size, n_walkers), np.object_)
status[:, :] = ''
for r in range(size):
status[r, cull_ind[np.where(cull_rank == r)[0]]] = 'c_t'
if ns_args['debug'] >= 10:
initial_PE_loc = [eval_energy(at, do_KE=False) for at in walkers]
initial_E_loc = [eval_energy(at) for at in walkers]
if comm is not None:
initial_PE = np.array(comm.allgather(initial_PE_loc)).flatten()
initial_E = np.array(comm.allgather(initial_E_loc)).flatten()
else:
initial_PE = np.array(initial_PE_loc)
initial_E = np.array(initial_E_loc)
initial_changed = initial_PE[np.where(status.flatten() == 'c_t')]
initial_unchanged = initial_PE[np.where(status.flatten() == '')]
if ns_args['debug'] >= 30:
for r in range(len(status)):
print(print_prefix, ": initial status ", r,
[s for s in status[r, :]])
# find load balance by cloning on top of excess maxima
recv_ind = []
recv_rank = []
send_ind = []
send_rank = []
cull_inds_to_remove = []
if n_cull > 1: # send/recv for fixing load balance
# CHECK FOR RANDOMNESS ISSUES AND WHICH NODES ARE USED FOR CLONES
for r in range(size):
# maybe remote_r should be chosen completely randomly, rather
# than close to task of extra culled configs
for dr in np.array(list(zip(np.array(range(1, size)), -np.array(range(1, size))))).flatten():
if n_cull_of_rank[r] <= max_n_cull_per_task: # not too many that need to be culled on this rank
break
# this rank has too many to cull, must receive replacement from another node
remote_r = (r+dr) % size
if n_cull_of_rank[remote_r] < max_n_cull_per_task: # send from r+dr to r
n_transfer = min(n_cull_of_rank[r]-max_n_cull_per_task, max_n_cull_per_task-n_cull_of_rank[remote_r])
recv_rank.extend([r]*n_transfer)
send_rank.extend([remote_r]*n_transfer)
local_ind = np.where(status[r, :] == 'c_t')[0][0:n_transfer]
recv_ind.extend(local_ind)
remote_ind = np.where(status[remote_r, :] == '')[0][0:n_transfer]
send_ind.extend(remote_ind)
status[r, local_ind] = 'c_s'
status[remote_r, remote_ind] = 'c_t_a'
n_cull_of_rank[r] -= n_transfer
n_cull_of_rank[remote_r] += n_transfer
# save local random state, and switch to common one
rng.switch_to_common()
# select clones
for r in range(size):
list_clone_target = np.where(status[r, :] == 'c_t')[0]
# assign clones
n_remaining_clones = len(list_clone_target)
while n_remaining_clones > 0:
remote_r = rng.int_uniform(0, size)
n_avail_remote = sum(status[remote_r, :] == '')
if n_avail_remote > 0: # something is available on remote_r
# send from random avail walker on remote_r to clone_target on r
n_transfer = min(n_remaining_clones, n_avail_remote)
# set ranks
send_rank.extend([remote_r]*n_transfer)
recv_rank.extend([r]*n_transfer)
# set indices
r_is = []
for ii in range(n_transfer):
r_i = rng.int_uniform(0, n_walkers)
while status[remote_r, r_i] != '':
r_i = rng.int_uniform(0, n_walkers)
# now r_i should be something with status ''
status[remote_r, r_i] = 'c_s'
r_is.append(r_i)
send_ind.extend(r_is)
status[r, list_clone_target[0:n_transfer]] = 'c_t_a'
recv_ind.extend(list_clone_target[0:n_transfer])
if n_transfer < len(list_clone_target):
list_clone_target = list_clone_target[n_transfer:]
n_remaining_clones -= n_transfer
if ns_args['debug'] >= 20:
print(print_prefix, "%30s" % ": LOOP_TE POST_LOC_CLONE 15 ", i_ns_step, ["%.10f" % eval_energy(at) for at in walkers])
print(print_prefix, "%30s" % ": LOOP_PE POST_LOC_CLONE 16 ", i_ns_step, ["%.10f" % eval_energy(at, do_KE=False) for at in walkers])
print(print_prefix, "%30s" % ": LOOP_X POST_LOC_CLONE 17 ", i_ns_step, ["%.10f" % at.positions[0, 0] for at in walkers])
# make into numpy arrays so that mathematical operations will work
send_rank = np.array(send_rank)
send_ind = np.array(send_ind)
recv_rank = np.array(recv_rank)
recv_ind = np.array(recv_ind)
if ns_args['debug'] >= 10:
if rank == 0:
for i in range(len(send_rank)):
print(print_prefix, "send from ", send_rank[i], send_ind[i],
" to ", recv_rank[i], recv_ind[i])
# save new common state, and restore to local state
rng.switch_to_local()
if n_cull == 1:
if send_rank[0] == recv_rank[0] and send_rank[0] == rank: # local copy
walkers[recv_ind[0]].set_positions(walkers[send_ind[0]].get_positions())
walkers[recv_ind[0]].set_cell(walkers[send_ind[0]].get_cell())
if movement_args['do_velocities']:
walkers[recv_ind[0]].set_velocities(walkers[send_ind[0]].get_velocities())
if movement_args['do_GMC']:
walkers[recv_ind[0]].arrays['GMC_direction'][:, :] = walkers[send_ind[0]].arrays['GMC_direction']
if ns_args['n_extra_data'] > 0:
walkers[recv_ind[0]].arrays['ns_extra_data'][...] = walkers[send_ind[0]].arrays['ns_extra_data']
if ns_args['swap_atomic_numbers']:
walkers[recv_ind[0]].set_atomic_numbers(walkers[send_ind[0]].get_atomic_numbers())
if movement_args['do_velocities']:
walkers[recv_ind[0]].set_masses(walkers[send_ind[0]].get_masses())
if ns_args['track_configs']:
walkers[recv_ind[0]].info['config_ind'] = walkers[send_ind[0]].info['config_ind']
walkers[recv_ind[0]].info['from_config_ind'] = walkers[send_ind[0]].info['from_config_ind']
walkers[recv_ind[0]].info['config_ind_time'] = walkers[send_ind[0]].info['config_ind_time']
walkers[recv_ind[0]].info['ns_energy'] = eval_energy(walkers[recv_ind[0]])
if ns_args['debug'] >= 10 and size <= 1:
walkers[recv_ind[0]].info['n_walks'] = 0
else: # need send/recv
n_send = 3*(n_atoms + 3)
if movement_args['do_velocities']:
n_send += 3*n_atoms
if movement_args['do_GMC']:
n_send += 3*n_atoms
if ns_args['n_extra_data'] > 0:
n_send += ns_args['n_extra_data']*n_atoms
if ns_args['swap_atomic_numbers']:
n_send += n_atoms # Z
if movement_args['do_velocities']:
n_send += n_atoms # mass
if ns_args['track_configs']:
n_send += 3
buf = np.zeros(n_send)
if send_rank[0] == rank: # only one config is sent/received
buf_o = 0
buf[buf_o:buf_o+3*n_atoms] = walkers[send_ind[0]].get_positions().reshape((3*n_atoms)); buf_o += 3*n_atoms
buf[buf_o:buf_o+3*3] = walkers[send_ind[0]].get_cell().reshape((3*3)); buf_o += 3*3
if movement_args['do_velocities']:
buf[buf_o:buf_o+3*n_atoms] = walkers[send_ind[0]].get_velocities().reshape((3*n_atoms)); buf_o += 3*n_atoms
if movement_args['do_GMC']:
buf[buf_o:buf_o+3*n_atoms] = walkers[send_ind[0]].arrays['GMC_direction'].reshape((3*n_atoms)); buf_o += 3*n_atoms
if ns_args['n_extra_data'] > 0:
buf[buf_o:buf_o+ns_args['n_extra_data']*n_atoms] = walkers[send_ind[0]].arrays['ns_extra_data'].reshape((ns_args['n_extra_data']*n_atoms)); buf_o += ns_args['n_extra_data']*n_atoms
if ns_args['swap_atomic_numbers']:
buf[buf_o:buf_o+n_atoms] = walkers[send_ind[0]].get_atomic_numbers(); buf_o += n_atoms
if movement_args['do_velocities']:
buf[buf_o:buf_o+n_atoms] = walkers[send_ind[0]].get_masses(); buf_o += n_atoms
if ns_args['track_configs']:
buf[buf_o] = walkers[send_ind[0]].info['config_ind']; buf_o += 1
buf[buf_o] = walkers[send_ind[0]].info['from_config_ind']; buf_o += 1
buf[buf_o] = walkers[send_ind[0]].info['config_ind_time']; buf_o += 1
comm.Send([buf, MPI.DOUBLE], dest=recv_rank[0], tag=100)
elif recv_rank[0] == rank:
comm.Recv([buf, MPI.DOUBLE], source=send_rank[0], tag=100)
buf_o = 0
walkers[recv_ind[0]].set_positions(buf[buf_o:buf_o+3*n_atoms].reshape((n_atoms, 3))); buf_o += 3*n_atoms
walkers[recv_ind[0]].set_cell(buf[buf_o:buf_o+3*3].reshape((3, 3))); buf_o += 3*3
if movement_args['do_velocities']:
walkers[recv_ind[0]].set_velocities(buf[buf_o:buf_o+3*n_atoms].reshape((n_atoms, 3))); buf_o += 3*n_atoms
if movement_args['do_GMC']:
walkers[recv_ind[0]].arrays['GMC_direction'][:, :] = buf[buf_o:buf_o+3*n_atoms].reshape((n_atoms, 3)); buf_o += 3*n_atoms
if ns_args['n_extra_data'] > 0:
walkers[recv_ind[0]].arrays['ns_extra_data'][...] = buf[buf_o:buf_o+3*n_atoms].reshape(walkers[recv_ind[0]].arrays['ns_extra_data'].shape); buf_o += ns_args['n_extra_data']*n_atoms
if ns_args['swap_atomic_numbers']:
walkers[recv_ind[0]].set_atomic_numbers(buf[buf_o:buf_o+n_atoms].astype(int)); buf_o += n_atoms
if movement_args['do_velocities']:
walkers[recv_ind[0]].set_masses(buf[buf_o:buf_o+n_atoms]); buf_o += n_atoms
if ns_args['track_configs']:
walkers[recv_ind[0]].info['config_ind'] = int(buf[buf_o]); buf_o += 1
walkers[recv_ind[0]].info['from_config_ind'] = int(buf[buf_o]); buf_o += 1
walkers[recv_ind[0]].info['config_ind_time'] = int(buf[buf_o]); buf_o += 1
walkers[recv_ind[0]].info['ns_energy'] = eval_energy(walkers[recv_ind[0]])
else: # complicated construction of sending/receiving buffers
# figure out how much is sent per config
n_data_per_config = 1+3*(n_atoms + 3)
if movement_args['do_velocities']:
n_data_per_config += 3*n_atoms
if movement_args['do_GMC']:
n_data_per_config += 3*n_atoms
if ns_args['n_extra_data'] > 0:
n_data_per_config += ns_args['n_extra_data']*n_atoms
if ns_args['swap_atomic_numbers']:
n_send += n_atoms # Z
if movement_args['do_velocities']:
n_send += n_atoms # mass
if ns_args['track_configs']:
n_data_per_config += 3
# figure out send counts
send_count = [0] * size
for i in np.where(send_rank == rank)[0]:
r_recv = recv_rank[i]
send_count[r_recv] += n_data_per_config
# figure out send displacements
send_displ = [0] * size
send_displ[0] = 0
for i in range(1,size):
send_displ[i] = send_displ[i-1] + send_count[i-1]
# create empty buffer for sending
send_count_tot = sum(send_count)
send_data = np.zeros(send_count_tot)
# copy data to be sent to buffer
send_displ_t = list(send_displ)
for i in np.where(send_rank == rank)[0]:
r_recv = recv_rank[i]
i_send = send_ind[i]
data_o = send_displ_t[r_recv]
send_data[data_o] = walkers[i_send].info['ns_energy']; data_o += 1
send_data[data_o:data_o+3*n_atoms] = walkers[i_send].get_positions().reshape( (3*n_atoms) ); data_o += 3*n_atoms
send_data[data_o:data_o+3*3] = walkers[i_send].get_cell().reshape( (3*3) ); data_o += 3*3
if movement_args['do_velocities']:
send_data[data_o:data_o+3*n_atoms] = walkers[i_send].get_velocities().reshape( (3*n_atoms) ); data_o += 3*n_atoms
if movement_args['do_GMC']:
send_data[data_o:data_o+3*n_atoms] = walkers[i_send].arrays['GMC_direction'].reshape( (3*n_atoms) ); data_o += 3*n_atoms
if ns_args['n_extra_data'] > 0:
send_data[data_o:data_o+ns_args['n_extra_data']*n_atoms] = walkers[i_send].arrays['ns_extra_data'].reshape( (ns_args['n_extra_data']*n_atoms) ); data_o += ns_args['n_extra_data']*n_atoms
if ns_args['swap_atomic_numbers']:
send_data[data_o:data_o+n_atoms] = walkers[i_send].get_atomic_numbers(); data_o += n_atoms
if movement_args['do_velocities']:
send_data[data_o:data_o+n_atoms] = walkers[i_send].get_masses(); data_o += n_atoms
if ns_args['track_configs']:
send_data[data_o] = walkers[i_send].info['config_ind']; data_o += 1
send_data[data_o] = walkers[i_send].info['from_config_ind']; data_o += 1
send_data[data_o] = walkers[i_send].info['config_ind_time']; data_o += 1
send_displ_t[r_recv] = data_o
# figure out recv counts
recv_count = [0] * size
for i in np.where(recv_rank == rank)[0]:
r_send = send_rank[i]
recv_count[r_send] += n_data_per_config
# figure out recv displacements
recv_displ = [0] * size
recv_displ[0] = 0
for i in range(1,size):
recv_displ[i] = recv_displ[i-1] + recv_count[i-1]
# create empty buffer for receiving
recv_count_tot = sum(recv_count)
recv_data = np.zeros(recv_count_tot)
# do communications
if comm is not None:
send_buf = [send_data, send_count, send_displ, MPI.DOUBLE]
recv_buf = (recv_data, recv_count, recv_displ, MPI.DOUBLE)
comm.Alltoallv(send_buf, recv_buf)
else:
recv_data = send_data.copy()
# copy data from recv buffer to walkers
recv_displ_t = list(recv_displ)
for i in np.where(recv_rank == rank)[0]:
r_send = send_rank[i]
i_recv = recv_ind[i]
data_o = recv_displ_t[r_send]
walkers[i_recv].info['ns_energy'] = recv_data[data_o]; data_o += 1
walkers[i_recv].set_positions( recv_data[data_o:data_o+3*n_atoms].reshape( (n_atoms, 3) )); data_o += 3*n_atoms
walkers[i_recv].set_cell( recv_data[data_o:data_o+3*3].reshape( (3, 3) )); data_o += 3*3
if movement_args['do_velocities']:
walkers[i_recv].set_velocities( recv_data[data_o:data_o+3*n_atoms].reshape( (n_atoms, 3) )); data_o += 3*n_atoms
if movement_args['do_GMC']:
walkers[i_recv].arrays['GMC_direction'][:,:] = recv_data[data_o:data_o+3*n_atoms].reshape( (n_atoms, 3) ); data_o += 3*n_atoms
if ns_args['n_extra_data'] > 0:
walkers[i_recv].arrays['ns_extra_data'][...] = recv_data[data_o:data_o+ns_args['n_extra_data']*n_atoms].reshape( walkers[i_recv].arrays['ns_extra_data'].shape ); data_o += ns_args['n_extra_data']*n_atoms
if ns_args['swap_atomic_numbers']:
walkers[i_recv].set_atomic_numbers(recv_data[data_o:data_o+n_atoms].astype(int)); data_o += n_atoms
if movement_args['do_velocities']:
walkers[i_recv].set_masses(recv_data[data_o:data_o+n_atoms]); data_o += n_masses
if ns_args['track_configs']:
walkers[i_recv].info['config_ind'] = int(recv_data[data_o]); data_o += 1
walkers[i_recv].info['from_config_ind'] = int(recv_data[data_o]); data_o += 1
walkers[i_recv].info['config_ind_time'] = int(recv_data[data_o]); data_o += 1
recv_displ_t[r_send] = data_o
if ns_args['debug'] >= 20:
print(print_prefix, "%30s" % ": LOOP_TE POST_CLONE 20 ", i_ns_step, ["%.10f" % eval_energy(at) for at in walkers])
print(print_prefix, "%30s" % ": LOOP_PE POST_CLONE 21 ", i_ns_step, ["%.10f" % eval_energy(at, do_KE=False) for at in walkers])
print(print_prefix, "%30s" % ": LOOP_X POST_CLONE 22 ", i_ns_step, ["%.10f" % at.positions[0, 0] for at in walkers])
if ns_args['track_configs']:
# loop over _all_ clone targets and increment cur_config_ind,
# setting appropriate configs' new config_ind as needed
for r in range(size):
clone_walk_ind = np.where(status[r, :] == 'c_t_a')[0]
for i_at in clone_walk_ind:
if r == rank:
walkers[i_at].info['from_config_ind'] = walkers[i_at].info['config_ind']
walkers[i_at].info['config_ind'] = cur_config_ind
walkers[i_at].info['config_ind_time'] = i_ns_step
cur_config_ind += 1
# move cloned walkers
if (i_ns_step == start_first_iter and movement_args['full_auto_step_sizes']):
# set initial step sizes. Performed here since this is the first time all the arrays are in place
conf_pre=walkers[0].copy()
conf_pre.calc = walkers[0].calc
move_args_pre=deepcopy(movement_args)
walk_stats_pre=walk_single_walker(conf_pre, move_args_pre, Emax_of_step, KEmax)
delta_step_size_setting_duration = full_auto_set_stepsizes(walkers, walk_stats_pre, movement_args, comm, Emax_of_step, KEmax, size)
total_step_size_setting_duration += delta_step_size_setting_duration
step_size_setting_duration += delta_step_size_setting_duration
del(walk_stats_pre)
del(move_args_pre)
del(conf_pre)
sys.stdout.flush()
# walk clone targets
if ns_args['debug'] >= 4:
for i in np.where(status[rank, :] == 'c_s')[0]:
print(print_prefix, "INFO: 30 clone source ", rank, i)
clone_walk_ind = np.where(status[rank, :] == 'c_t_a')[0]
for i_at in clone_walk_ind:
if ns_args['debug'] >= 4:
print(print_prefix, "INFO: 40 WALK clone_target ", rank, i_at)
walk_stats = walk_single_walker(walkers[i_at], movement_args,
Emax_of_step, KEmax)
walkers[i_at].info['last_walked_iter_clone'] = i_ns_step
# if tracking all configs, save this one that has been walked
if track_traj_io is not None:
walkers[i_at].info['iter'] = i_ns_step
ase.io.write(track_traj_io, walkers[i_at], format=ns_args['config_file_format'])
#print("WALK on rank ", rank, "at iteration ", i_ns_step, " walker ", i_at )
if ns_args['debug'] >= 10 and size <= 1:
walkers[i_at].info['n_walks'] += movement_args['n_model_calls']
accumulate_stats(walk_stats_adjust, walk_stats)
accumulate_stats(walk_stats_monitor, walk_stats)
sys.stdout.flush()
if ns_args['debug'] >= 20:
print(print_prefix, "%30s" % ": LOOP_TE POST_CLONE_WALK 25 ", i_ns_step, ["%.10f" % eval_energy(at) for at in walkers])
print(print_prefix, "%30s" % ": LOOP_PE POST_CLONE_WALK 26 ", i_ns_step, ["%.10f" % eval_energy(at, do_KE=False) for at in walkers])
print(print_prefix, "%30s" % ": LOOP_X POST_CLONE_WALK 27 ", i_ns_step, ["%.10f" % at.positions[0, 0] for at in walkers])
# check that everything that should have been changed has, and things
# that shouldn't have, haven't
if ns_args['debug'] >= 10:
final_PE_loc = [eval_energy(at, do_KE=False) for at in walkers]
final_E_loc = [eval_energy(at) for at in walkers]
if comm is not None:
final_PE = np.array(comm.allgather(final_PE_loc)).flatten()
final_E = np.array(comm.allgather(final_E_loc)).flatten()
else:
final_PE = final_PE_loc
final_E = final_E_loc
if rank == 0:
final_status = status.flatten()
for e in initial_unchanged:
if e not in final_PE:
print("initial_PE ", initial_PE)
print("final_PE ", final_PE)
print("initial_E ", initial_E)
print("final_E ", final_E)
print("final_status ", final_status)
print("WARNING: energy that should have been unchanged ", e," missing from final energies")
for e in initial_changed:
if e in final_PE:
print("initial_PE ", initial_PE)
print("final_PE ", final_PE)
print("initial_E ", initial_E)
print("final_E ", final_E)
print("final_status ", final_status)
print("WARNING: energy that should have been changed ", e," still there in final energies")
# walk extras
if not ns_args['no_extra_walks_at_all']:
r_i = rng.int_uniform(0, n_walkers)
# WARNING: this may select walkers for extra walks multiple
# times, yet never re-walk ones that were walked as clone
# targets
while status[rank, r_i] != '' and status[rank, r_i] != 'c_s':
r_i = rng.int_uniform(0, n_walkers)
if ns_args['debug'] >= 4:
print(print_prefix, "INFO: 50 WALK extra ", rank, r_i)
walk_stats = walk_single_walker(walkers[r_i], movement_args,
Emax_of_step, KEmax)
walkers[r_i].info['last_walked_iter_extra'] = i_ns_step
# if tracking all configs, save this one that has been walked
if track_traj_io is not None:
walkers[i_at].info['iter'] = i_ns_step
ase.io.write(track_traj_io, walkers[i_at],
format=ns_args['config_file_format'])
# print("WALK EXTRA on rank ", rank, "at iteration ", i_ns_step,
# " walker ", r_i)
if ns_args['debug'] >= 10 and size <= 1:
#walkers[r_i].info['n_walks'] += movement_args['n_steps'] # LIVIA-this gives error, does not exist
walkers[r_i].info['n_walks'] += movement_args['atom_traj_len']
accumulate_stats(walk_stats_adjust, walk_stats)
accumulate_stats(walk_stats_monitor, walk_stats)
monitored_this_step = False
if movement_args['monitor_step_interval'] != 0 and i_ns_step % abs(movement_args['monitor_step_interval']) == abs(movement_args['monitor_step_interval'])-1:
adjust_step_sizes(walk_stats_monitor, movement_args, comm, monitor_only=True)
zero_stats(walk_stats_monitor, movement_args)
monitored_this_step = True
if movement_args['adjust_step_interval'] != 0 and i_ns_step % abs(movement_args['adjust_step_interval']) == abs(movement_args['adjust_step_interval'])-1:
if (not movement_args['full_auto_step_sizes']):
adjust_step_sizes(walk_stats_adjust, movement_args, comm, do_print_rate=(not monitored_this_step))
else:
delta_step_size_setting_duration = full_auto_set_stepsizes(walkers, walk_stats_adjust, movement_args, comm, Emax_of_step, KEmax, size)
total_step_size_setting_duration += delta_step_size_setting_duration
step_size_setting_duration += delta_step_size_setting_duration
zero_stats(walk_stats_adjust, movement_args)
if ns_args['debug'] >= 20:
print(print_prefix, "%30s" % ": LOOP_TE END 30 ", i_ns_step, ["%.10f" % eval_energy(at) for at in walkers])
print(print_prefix, "%30s" % ": LOOP_PE END 31 ", i_ns_step, ["%.10f" % eval_energy(at, do_KE=False) for at in walkers])
print(print_prefix, "%30s" % ": LOOP_X END 32 ", i_ns_step, ["%.10f" % at.positions[0, 0] for at in walkers])
if ns_args['debug'] >= 30:
for r in range(len(status)):
print(print_prefix, ": final status ", r, [s for s in status[r, :]])
if (rank == 0) and ((ns_args['snapshot_interval'] > 0 and i_ns_step > 0 and i_ns_step % ns_args['snapshot_interval'] == 0) or
(ns_args['snapshot_seq_pairs'] and i_ns_step > 1 and i_ns_step%ns_args['snapshot_interval'] == 1) or
(ns_args['snapshot_time'] > 0 and time.time()-last_snapshot_time > ns_args['snapshot_time'])):
do_snapshot = True
else:
do_snapshot = False
if comm is not None:
do_snapshot = comm.bcast(do_snapshot, root=0)
if do_snapshot:
save_snapshot(i_ns_step)
last_snapshot_time = time.time()
clean_prev_snapshot(pprev_snapshot_iter)
pprev_snapshot_iter = prev_snapshot_iter
prev_snapshot_iter = i_ns_step
if ns_analyzers is not None:
for (ns_analyzer, ns_analyzer_interval) in ns_analyzers:
if ns_analyzer_interval > 0 and (i_ns_step+1)%ns_analyzer_interval == 0:
ns_analyzer.analyze(walkers, i_ns_step, "NS_loop %d" % i_ns_step)
i_ns_step += 1
### END OF MAIN LOOP
# flush remaining traj configs
for at in traj_walker_list:
ase.io.write(traj_io, at, parallel=False, format=ns_args['config_file_format'])
traj_io.flush()
traj_walker_list = []
if ns_args['E_dump_interval'] > 0:
if comm is not None:
E_dump_list_all = np.array(comm.allgather(E_dump_list))
else:
E_dump_list_all = np.array(E_dump_list)
if rank == 0:
for i in range(E_dump_list_all.shape[1]):
E_dump_io.write("step %d\n" % E_dump_list_times[i])
if len(E_dump_list_all.shape) == 3:
np.savetxt(E_dump_io, E_dump_list_all[:,i,:])
else:
np.savetxt(E_dump_io, E_dump_list_all[i,:])
E_dump_io.flush()
cur_time = time.time()
if rank == 0:
print( "LOOP TIME total ",cur_time-initial_time-total_step_size_setting_duration, " per iter ", (cur_time-initial_time-total_step_size_setting_duration)/(i_ns_step+1))
print( "TIME SPENT SETTING STEP SIZES total ",total_step_size_setting_duration)
return i_ns_step-1
[docs]
def main():
""" Main function """
global movement_args
global ns_args, start_first_iter
global max_n_cull_per_task
global size, rank, comm, rng, np, sys, ns_analyzers
global n_cull, n_walkers, n_walkers_per_task
global n_extra_walk_per_task
global do_calc_ASE, do_calc_lammps, do_calc_internal, do_calc_fortran
global energy_io, traj_io, walkers
global n_atoms, KEmax, pot
global MPI, f_MC_MD
global track_traj_io, cur_config_ind
global E_dump_io
global print_prefix
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)
use_mpi = True
if len(sys.argv) == 2:
if sys.argv[1] == "-no_mpi":
use_mpi = False
else:
usage()
sys.exit(1)
# initialize mpi
comm = None
calculator_comm = None
rank = 0
size = 1
if use_mpi:
print("INFO: use_mpi true, importing mpi4py module")
try:
from mpi4py import MPI
except:
sys.stderr.write("Failed to import mpi4py\n")
sys.exit(10)
comm = MPI.COMM_WORLD
calculator_comm = MPI.COMM_SELF
if use_mpi:
try:
rank = comm.Get_rank() # id of process
size = comm.Get_size() # number of processes
except:
exit_error("Failed to get rank or size\n", 10)
if comm is not None:
print("comm ", comm, " size ", size, " rank ", rank)
# read inputs on root, then bcast
if rank == 0:
lines = sys.stdin.readlines()
if len(lines) == 0:
try:
infile = open("ns_inputs", "r")
except:
exit_error("Failed to read ns_inputs file\n", 1)
lines = infile.readlines()
args = {}
if rank == 0:
for line in lines:
if re.match("\s*(#.*)?$", line):
continue
matches = re.match("\s*(\S+)\s*=\s*(.*\S)", line)
if matches is None:
exit_error("Failed to parse line '%s'" % line, 1)
args[matches.group(1)] = matches.group(2)
else:
args = None
if comm is not None:
args = comm.bcast(args, root=0) # send args to other processes
#DOC ``main``: parse arguments
# parse args
ns_args = {}
ns_args['check_memory'] = str_to_logical(args.pop('check_memory', 'F'))
if ns_args['check_memory'] and (comm is None or comm.rank == 0):
check_memory.active = True
# convert from strings to actual args
try:
ns_args['n_walkers'] = int(args.pop('n_walkers'))
except:
exit_error("need number of walkers n_walkers\n", 1)
ns_args['n_cull'] = int(args.pop('n_cull', 1))
ns_args['n_iter_times_fraction_killed'] = float(args.pop('n_iter_times_fraction_killed', -1))
if ns_args['n_iter_times_fraction_killed'] > 0:
ns_args['n_iter'] = int(round(ns_args['n_iter_times_fraction_killed']/(float(ns_args['n_cull'])/float(ns_args['n_walkers']))))
else:
ns_args['n_iter'] = -1
ns_args['converge_down_to_T'] = float(args.pop('converge_down_to_T', -1))
if ns_args['n_iter'] <= 0 and ns_args['converge_down_to_T'] <= 0:
exit_error("need either n_iter_times_fraction_killed or converge_down_to_T", 1)
ns_args['T_estimate_finite_diff_lag'] = int(args.pop('T_estimate_finite_diff_lag', 1000))+1
try:
ns_args['min_Emax'] = float(args.pop('min_Emax'))
except:
ns_args['min_Emax'] = None
ns_args['start_species'] = args.pop('start_species', None)
ns_args['start_config_file'] = args.pop('start_config_file', None)
if ns_args['start_species'] is None and ns_args['start_config_file'] is None:
exit_error("always need start_species or start_config_file, even if restart_file is specified\n",1)
if ns_args['start_species'] is not None and ns_args['start_config_file'] is not None:
exit_error("can't specify both start_species and start_config_file\n",1)
ns_args['restart_file'] = args.pop('restart_file', 'AUTO')
ns_args['max_volume_per_atom'] = float(args.pop('max_volume_per_atom', 1.0e3))
ns_args['min_volume_per_atom'] = float(args.pop('min_volume_per_atom', 1.0))
ns_args['out_file_prefix'] = args.pop('out_file_prefix', '')
if ns_args['out_file_prefix'] != '':
ns_args['out_file_prefix'] += '.'
ns_args['profile'] = int(args.pop('profile', -1))
ns_args['debug'] = int(args.pop('debug', -1))
ns_args['snapshot_interval'] = int(args.pop('snapshot_interval', -1))
ns_args['snapshot_time'] = int(args.pop('snapshot_time', 3600))
ns_args['snapshot_per_parallel_task'] = str_to_logical(args.pop('snapshot_per_parallel_task', 'T'))
ns_args['snapshot_seq_pairs'] = str_to_logical(args.pop('snapshot_seq_pairs', "F"))
ns_args['snapshot_clean'] = str_to_logical(args.pop('snapshot_clean', "T"))
ns_args['random_initialise_pos'] = str_to_logical(args.pop('random_initialise_pos', "T"))
ns_args['random_initialise_cell'] = str_to_logical(args.pop('random_initialise_cell', "T"))
ns_args['LAMMPS_molecular_info'] = str_to_logical(args.pop('LAMMPS_molecular_info', "T"))
ns_args['initial_walk_N_walks'] = int(args.pop('initial_walk_N_walks', 0))
ns_args['initial_walk_adjust_interval'] = int(args.pop('initial_walk_adjust_interval', 1))
ns_args['initial_walk_Emax_offset_per_atom'] = float(args.pop('initial_walk_Emax_offset_per_atom', 1))
ns_args['initial_walk_only'] = str_to_logical(args.pop('initial_walk_only', 'F'))
ns_args['traj_interval'] = int(args.pop('traj_interval', 100))
ns_args['E_dump_interval'] = int(args.pop('E_dump_interval', -1))
ns_args['delta_random_seed'] = int(args.pop('delta_random_seed', -1))
ns_args['n_extra_walk_per_task'] = int(args.pop('n_extra_walk_per_task', 0))
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 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'))
#VGF end nn parameters
# surely there's a cleaner way of doing this?
try:
ns_args['start_energy_ceiling_per_atom'] = float(args.pop('start_energy_ceiling_per_atom'))
except:
ns_args['start_energy_ceiling_per_atom'] = None
try:
ns_args['start_energy_ceiling'] = float(args.pop('start_energy_ceiling'))
except:
ns_args['start_energy_ceiling'] = None
if ns_args['start_energy_ceiling_per_atom'] is not None and ns_args['start_energy_ceiling'] is not None:
# conflict
exit_error("got both start_energy_ceiling and start_energy_ceiling_per_atom\n", 1)
elif ns_args['start_energy_ceiling_per_atom'] is None and ns_args['start_energy_ceiling'] is None:
# neither specified, use default
ns_args['start_energy_ceiling_per_atom'] = 1.0e9
elif ns_args['start_energy_ceiling'] is not None and rank == 0:
# warn on deprecated feature
sys.stderr.write("WARNING: got DEPRECATED start_energy_ceiling\n")
ns_args['random_init_max_n_tries'] = int(args.pop('random_init_max_n_tries', 100))
ns_args['KEmax_max_T'] = float(args.pop('KEmax_max_T', -1))
ns_args['kB'] = float(args.pop('kB', 8.6173324e-5)) # eV/K
# parse energy_calculator
ns_args['energy_calculator'] = args.pop('energy_calculator', 'fortran')
do_calc_ASE = False
do_calc_lammps = False
do_calc_internal = False
do_calc_fortran = False
if ns_args['energy_calculator'] == 'ASE':
try:
ns_args['ASE_calc_module'] = args.pop('ASE_calc_module')
except:
exit_error("need ASE calculator module name ASE_calc_module\n",1)
do_calc_ASE=True
elif ns_args['energy_calculator'] == 'lammps':
try:
from lammpslib import LAMMPSlib
except:
exit_error("energy_calculator=lammps and failed to import lammpslib module\n", 1)
do_calc_lammps=True
try:
ns_args['LAMMPS_init_cmds'] = args.pop('LAMMPS_init_cmds')
except:
exit_error("need LAMMPS initialization commands LAMMPS_init_cmds\n", 1)
ns_args['LAMMPS_name'] = args.pop('LAMMPS_name', os.environ.get('LAMMPS_name', ''))
ns_args['LAMMPS_serial'] = str_to_logical(args.pop('LAMMPS_serial', 'T'))
ns_args['LAMMPS_header'] = args.pop('LAMMPS_header', 'units metal; atom_style atomic; atom_modify map array sort 0 0')
ns_args['LAMMPS_header_extra'] = args.pop('LAMMPS_header_extra', '')
ns_args['LAMMPS_atom_types'] = None
ns_args['LAMMPS_fix_gmc'] = str_to_logical(args.pop('LAMMPS_fix_gmc', "F"))
LAMMPS_atom_types = args.pop('LAMMPS_atom_types', '')
if len(LAMMPS_atom_types) > 0:
if LAMMPS_atom_types == 'TYPE_EQUALS_Z':
ns_args['LAMMPS_atom_types'] = LAMMPS_atom_types
else:
ns_args['LAMMPS_atom_types'] = {}
for type_pair in [s.strip() for s in LAMMPS_atom_types.split(',')]:
f = type_pair.split()
ns_args['LAMMPS_atom_types'][f[0]] = int(f[1])
else:
exit_error("LAMMPS_atom_types is mandatory if calculator type is LAMMPS\n", 1)
elif ns_args['energy_calculator'] == 'internal':
do_calc_internal=True
elif ns_args['energy_calculator'] == 'fortran':
import fortranMCMDpy
do_calc_fortran=True
try:
ns_args['FORTRAN_model'] = args.pop('FORTRAN_model')
except:
exit_error("need FORTRAN model FORTRAN_model\n", 1)
ns_args['FORTRAN_model_params'] = args.pop('FORTRAN_model_params', '0')
f_MC_MD = fortranMCMDpy.fortran_MC_MD(ns_args['FORTRAN_model'])
params = np.array([float(x) for x in ns_args['FORTRAN_model_params'].split()])
f_MC_MD.init_model(params)
else:
exit_error("energy_calculator=%s unknown\n" % ns_args['energy_calculator'], 3)
ns_args['no_extra_walks_at_all'] = str_to_logical(args.pop('no_extra_walks_at_all', "F"))
ns_args['track_configs'] = str_to_logical(args.pop('track_configs', "F"))
ns_args['track_configs_write'] = str_to_logical(args.pop('track_configs_write', "F"))
ns_args['config_file_format'] = args.pop('config_file_format', 'extxyz')
ns_args['rng'] = args.pop('rng', 'numpy')
ns_args['ns_run_analyzers'] = args.pop('ns_run_analyzers', '')
if ns_args['rng'] == 'numpy':
rng = ns_rng.NsRngNumpy(ns_args['delta_random_seed'], comm)
# elif ns_args['rng'] == 'julia':
# import julia
# j = julia.Julia()
# rng = ns_rng.NsRngJulia(j)
elif ns_args['rng'] == 'rngstream':
import rngstream
rng = ns_rng.NsRngStream(ns_args['delta_random_seed'], comm)
elif ns_args['rng'] == 'internal':
rng = ns_rng.NsRngInternal(ns_args['delta_random_seed'], comm)
else:
exit_error("rng=%s unknown\n" % ns_args['rng'], 3)
if do_calc_fortran:
l_seed = f_MC_MD.seed_size()
seed = np.array([0] * l_seed, dtype=np.int32)
for i in range(l_seed):
# maybe we need a better way of getting max possible int
seed[i] = rng.int_uniform(1, sys.maxsize)
f_MC_MD.set_seed(seed)
ns_args['reproducible'] = str_to_logical(args.pop('reproducible', "F"))
if ns_args['reproducible']:
# reset seed after using some random numbers to generate fortran
# seed, so that fortran and non-fortran have the same seed
if ns_args['rng'] == 'numpy':
rng = ns_rng.NsRngNumpy(ns_args['delta_random_seed'], comm)
elif ns_args['rng'] == 'rngstream':
rng = ns_rng.NsRngStream(ns_args['delta_random_seed'], comm)
elif ns_args['rng'] == 'internal':
rng = ns_rng.NsRngInternal(ns_args['delta_random_seed'], comm)
else:
exit_error("rng=%s unknown\n" % ns_args['rng'], 3)
movement_args = {}
movement_args['n_model_calls_expected'] = int(args.pop('n_model_calls_expected', 0))
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
movement_args['n_atom_steps'] = int(args.pop('n_atom_steps', 1))
movement_args['atom_traj_len'] = int(args.pop('atom_traj_len', 8))
movement_args['atom_traj_len_cost_multiplier'] = int(args.pop('atom_traj_len_cost_multiplier', 1))
movement_args['break_up_atom_traj'] = str_to_logical(args.pop('break_up_atom_traj', "T"))
if movement_args['n_atom_steps'] == 0:
movement_args['n_atom_steps_per_call'] = 0
movement_args['n_atom_steps_n_calls'] = 0
else:
if movement_args['break_up_atom_traj']:
movement_args['n_atom_steps_per_call'] = 1
movement_args['n_atom_steps_n_calls'] = movement_args['n_atom_steps']
else:
movement_args['n_atom_steps_per_call'] = movement_args['n_atom_steps']
movement_args['n_atom_steps_n_calls'] = 1
movement_args['n_cell_volume_steps'] = int(args.pop('n_cell_volume_steps', 1))
movement_args['n_cell_shear_steps'] = int(args.pop('n_cell_shear_steps', 1))
movement_args['n_cell_stretch_steps'] = int(args.pop('n_cell_stretch_steps', 1))
movement_args['n_swap_steps'] = int(args.pop('n_swap_steps', 0))
movement_args['swap_max_cluster'] = int(args.pop('swap_max_cluster', 1))
movement_args['swap_r_cut'] = float(args.pop('swap_r_cut', 2.5))
movement_args['swap_cluster_probability_increment'] = float(args.pop('cluster_probability_increment', 0.75))
movement_args['swap_velo'] = str_to_logical(args.pop('swap_velo', "F"))
movement_args['no_swap_velo_fix_mag_alt'] = str_to_logical(args.pop('no_swap_velo_fix_mag_alt', "F"))
movement_args['n_semi_grand_steps'] = int(args.pop('n_semi_grand_steps', 0))
if movement_args['n_semi_grand_steps'] > 0:
try:
m = re.match('\s*{([^}]*)}\s*$', args.pop('semi_grand_potentials'))
semi_grand_potentials = m.group(1).split(',')
except:
exit_error('Got n_semi_grand_steps > 0 but no valid semi_grand_potentials\n', 3)
movement_args['semi_grand_potentials'] = {}
for semi_grand_item in semi_grand_potentials:
(Z_s, mu_s) = semi_grand_item.strip().split(":")
movement_args['semi_grand_potentials'][int(Z_s)] = float(mu_s)
else:
movement_args['semi_grand_potentials'] = None
if movement_args['n_model_calls_expected'] < 0:
movement_args['n_model_calls_expected'] = (movement_args['n_atom_steps']*movement_args['atom_traj_len']*movement_args['atom_traj_len_cost_multiplier'] +
movement_args['n_cell_volume_steps'] +
movement_args['n_cell_shear_steps'] +
movement_args['n_cell_stretch_steps'] +
movement_args['n_swap_steps'])
# initialize swap cluster size probabilities
movement_args['swap_probs'] = np.zeros((movement_args['swap_max_cluster']))
movement_args['swap_probs'][0] = 1.0
for i in range(1, movement_args['swap_max_cluster']):
movement_args['swap_probs'][i] = movement_args['swap_probs'][i-1] * movement_args['swap_cluster_probability_increment']
movement_args['swap_probs'] /= np.sum(movement_args['swap_probs'])
for i in range(1, movement_args['swap_max_cluster']):
movement_args['swap_probs'][i] = movement_args['swap_probs'][i] + movement_args['swap_probs'][i-1]
if (movement_args['n_model_calls_expected'] <= 0 and
movement_args['n_model_calls'] <= 0):
exit_error("Got all of n_model_calls* == 0\n", 3)
if (movement_args['n_atom_steps'] <= 0 and
movement_args['n_cell_volume_steps'] <= 0 and
movement_args['n_cell_shear_steps'] <= 0 and
movement_args['n_cell_stretch_steps'] <= 0 and
movement_args['n_swap_steps'] <= 0 and
movement_args['n_semi_grand_steps'] <= 0):
exit_error("Got all of n_steps_* == 0\n", 3)
movement_args['velo_traj_len'] = int(args.pop('velo_traj_len', 8))
try:
movement_args['atom_algorithm'] = args.pop('atom_algorithm')
except:
exit_error("Failed to read algorithm for atom motion atom_algorithm", 1)
if movement_args['atom_algorithm'] != 'MC' and movement_args['atom_algorithm'] != 'MD' and movement_args['atom_algorithm'] != 'GMC':
exit_error("Got unknown atom_algorithm '%s'\n" % movement_args['atom_algorithm'], 3)
movement_args['MC_atom_velocities'] = str_to_logical(args.pop('MC_atom_velocities', "F"))
movement_args['MC_atom_velocities_pre_perturb'] = str_to_logical(args.pop('MC_atom_velocities_pre_perturb', "F"))
movement_args['MC_atom_step_size'] = float(args.pop('MC_atom_step_size', 1.0))
movement_args['MC_atom_step_size_max'] = float(args.pop('MC_atom_step_size_max', 1.0))
movement_args['MC_atom_velo_step_size'] = float(args.pop('MC_atom_velo_step_size', 50.0))
movement_args['MC_atom_velo_step_size_max'] = float(args.pop('MC_atom_velo_step_size_max', 10000.0))
movement_args['MC_atom_uniform_rv'] = str_to_logical(args.pop('MC_atom_uniform_rv', "F"))
movement_args['MC_atom_Galilean'] = str_to_logical(args.pop('MC_atom_Galilean', "F"))
movement_args['GMC_no_reverse'] = str_to_logical(args.pop('GMC_no_reverse', "T"))
movement_args['do_velocities'] = (movement_args['atom_algorithm'] == 'MD' or movement_args['MC_atom_velocities'])
# atom_algorithm == GMC is just an alias for atom_algorithm = MC, MC_atom_Galilean = True
if movement_args['atom_algorithm'] == 'GMC':
movement_args['atom_algorithm'] = 'MC'
movement_args['MC_atom_Galilean'] = True
movement_args['python_MD'] = str_to_logical(args.pop('python_MD', "F"))
movement_args['MD_atom_velo_pre_perturb'] = str_to_logical(args.pop('MD_atom_velo_pre_perturb', "F"))
movement_args['MD_atom_velo_post_perturb'] = str_to_logical(args.pop('MD_atom_velo_post_perturb', "T"))
movement_args['MD_atom_velo_flip_accept'] = str_to_logical(args.pop('MD_atom_velo_flip_accept', "F"))
movement_args['atom_velo_rej_free_fully_randomize'] = str_to_logical(args.pop('atom_velo_rej_free_fully_randomize', "F"))
movement_args['atom_velo_rej_free_perturb_angle'] = float(args.pop('atom_velo_rej_free_perturb_angle', 0.3))
movement_args['MC_atom_velo_walk_rej_free'] = str_to_logical(args.pop('MC_atom_velo_walk_rej_free', "T"))
movement_args['MD_atom_timestep'] = float(args.pop('MD_atom_timestep', 0.1))
movement_args['MD_atom_timestep_max'] = float(args.pop('MD_atom_timestep_max', 2.0))
movement_args['MD_atom_energy_fuzz'] = float(args.pop('MD_atom_energy_fuzz', 1.0e-2))
movement_args['MD_atom_reject_energy_violation'] = str_to_logical(args.pop('MD_atom_reject_energy_violation', "F"))
movement_args['MC_cell_P'] = float(args.pop('MC_cell_P', 0.0))
movement_args['MC_cell_flat_V_prior'] = str_to_logical(args.pop('MC_cell_flat_V_prior', "F"))
default_value = ns_args['max_volume_per_atom']/20.0 # 5% of maximum allowed volume per atom
movement_args['MC_cell_volume_per_atom_step_size'] = float(args.pop('MC_cell_volume_per_atom_step_size', default_value))
movement_args['MC_cell_volume_per_atom_step_size_max'] = float(args.pop('MC_cell_volume_per_atom_step_size_max', 10.0*default_value)) # 50% of maximum allowed volume per atom
movement_args['MC_cell_volume_per_atom_prob'] = float(args.pop('MC_cell_volume_per_atom_prob', 1.0))
movement_args['MC_cell_stretch_step_size'] = float(args.pop('MC_cell_stretch_step_size', 0.1))
movement_args['MC_cell_stretch_step_size_max'] = float(args.pop('MC_cell_stretch_step_size_max', 1.0))
movement_args['MC_cell_stretch_prob'] = float(args.pop('MC_cell_stretch_prob', 1.0))
movement_args['MC_cell_shear_step_size'] = float(args.pop('MC_cell_shear_step_size', 0.1))
movement_args['MC_cell_shear_step_size_max'] = float(args.pop('MC_cell_shear_step_size_max', 1.0))
movement_args['MC_cell_shear_prob'] = float(args.pop('MC_cell_shear_prob', 1.0))
movement_args['MC_cell_min_aspect_ratio'] = float(args.pop('MC_cell_min_aspect_ratio', 0.8))
movement_args['cell_shape_equil_steps'] = int(args.pop('cell_shape_equil_steps', 1000))
try:
movement_args['monitor_step_interval'] = float(args.pop('monitor_step_interval'))
except:
movement_args['monitor_step_interval_times_fraction_killed'] = float(args.pop('monitor_step_interval_times_fraction_killed', 1))
movement_args['monitor_step_interval'] = int(round(movement_args['monitor_step_interval_times_fraction_killed']/(float(ns_args['n_cull'])/float(ns_args['n_walkers']))))
try:
movement_args['adjust_step_interval'] = float(args.pop('adjust_step_interval'))
except:
movement_args['adjust_step_interval_times_fraction_killed'] = float(args.pop('adjust_step_interval_times_fraction_killed', 1))
movement_args['adjust_step_interval'] = int(round(movement_args['adjust_step_interval_times_fraction_killed']/(float(ns_args['n_cull'])/float(ns_args['n_walkers']))))
# if movement_args['adjust_step_interval'] < 20:
# print("WARNING: step size adjustment would be done too often, at every ", movement_args['adjust_step_interval'], " iteration")
# print("WARNING: adjust_step_interval is increased to 20")
# movement_args['adjust_step_interval'] = 20
movement_args['full_auto_step_sizes'] = str_to_logical(args.pop('full_auto_step_sizes', "T"))
movement_args['MC_adjust_step_factor'] = float(args.pop('MC_adjust_step_factor', 1.5))
movement_args['MC_adjust_min_rate'] = float(args.pop('MC_adjust_min_rate', 0.25))
movement_args['MC_adjust_max_rate'] = float(args.pop('MC_adjust_max_rate', 0.75))
movement_args['GMC_adjust_min_rate'] = float(args.pop('GMC_adjust_min_rate', 0.25))
movement_args['GMC_adjust_max_rate'] = float(args.pop('GMC_adjust_max_rate', 0.75))
movement_args['GMC_dir_perturb_angle'] = float(args.pop('GMC_dir_perturb_angle', -1.0))
movement_args['GMC_dir_perturb_angle_during'] = float(args.pop('GMC_dir_perturb_angle_during', 0.0))
movement_args['MD_adjust_step_factor'] = float(args.pop('MD_adjust_step_factor', 1.1))
movement_args['MD_adjust_min_rate'] = float(args.pop('MD_adjust_min_rate', 0.50))
movement_args['MD_adjust_max_rate'] = float(args.pop('MD_adjust_max_rate', 0.95))
movement_args['2D'] = str_to_logical(args.pop('2D', "F"))
movement_args['keep_atoms_fixed'] = int(args.pop('keep_atoms_fixed', 0))
movement_args['apply_Z_wall'] = str_to_logical(args.pop('apply_Z_wall', "F"))
if movement_args['apply_Z_wall']:
# TODO: RY - the wall_dist should be a parameter and should allow user to change its value
movement_args['wall_dist'] = 10.00 # LIVIA - review this hard coded value
else:
movement_args['wall_dist'] = 0.0
if len(args) > 0:
exit_error(str(args)+"\nUnknown arguments read in\n", 2)
if rank == 0:
print("ns_args ", pprint.pformat(ns_args))
print("movement_args ", pprint.pformat(movement_args))
# initialize in-situ analyzers
try:
ns_analyzers=[]
for analyzer_str in ns_args['ns_run_analyzers'].split(";"):
try:
(analyzer_name, ns_analyzer_interval) = analyzer_str.split()
ns_analyzer_interval = int(ns_analyzer_interval)
try:
analyzer_module = importlib.import_module(analyzer_name.strip())
except:
analyzer_module = importlib.import_module("ns_run_analyzers."+analyzer_name.strip())
ns_analyzers.append((analyzer_module.NSAnalyzer(comm), ns_analyzer_interval))
if rank == 0:
print("Got NSAnalyzer from", analyzer_name.strip(), "ns_analyzer_interval", ns_analyzer_interval)
except:
if rank == 0:
print("Failed to get NSAnalyzer from", analyzer_name.strip())
except:
if rank == 0:
print("no ns_run_analyzers set")
ns_analyzers = None
# initialise potential
if do_calc_ASE:
pot = importlib.import_module(ns_args['ASE_calc_module']).calc
elif do_calc_internal or do_calc_fortran:
pass
elif do_calc_lammps:
init_cmds = [s.strip() for s in ns_args['LAMMPS_init_cmds'].split(';')]
header_cmds = [s.strip() for s in ns_args['LAMMPS_header'].split(';')]
header_extra_cmds = [s.strip() for s in ns_args['LAMMPS_header_extra'].split(';')]
if ns_args['LAMMPS_serial']:
lammps_comm = None
else:
lammps_comm = calculator_comm
if ns_args['debug'] >= 5:
pot = LAMMPSlib(lmpcmds=init_cmds, atom_types=ns_args['LAMMPS_atom_types'], log_file='lammps.%d.log' % rank, keep_alive=True, lammps_name=ns_args['LAMMPS_name'],
lammps_header=header_cmds, lammps_header_extra=header_extra_cmds, comm=lammps_comm, read_molecular_info=ns_args['LAMMPS_molecular_info'])
else:
pot = LAMMPSlib(lmpcmds=init_cmds, atom_types=ns_args['LAMMPS_atom_types'], keep_alive=True, lammps_name=ns_args['LAMMPS_name'],
lammps_header=header_cmds, lammps_header_extra=header_extra_cmds, comm=lammps_comm, read_molecular_info=ns_args['LAMMPS_molecular_info'])
if rank == 0:
print("PRE START_LAMMPS")
sys.stdout.flush()
pot.start_lammps() # so top level things like units will be set
if rank == 0:
print("POST START_LAMMPS")
sys.stdout.flush()
pot.first_propagate = True
else:
exit_error("Need some way of initializing calculator\n", 3)
# figure out numbers of local walkers
rank_of_walker = [0] * ns_args['n_walkers']
if size <= 1:
n_walkers = ns_args['n_walkers']
else:
n_walkers_per_task = ns_args['n_walkers']//size # using // ensures division gives an integer value
if n_walkers_per_task*size != ns_args['n_walkers']:
exit_error("number of walkers %d not divisible by number of MPI processes %d\n" % (ns_args['n_walkers'], size), 5)
last_walker = 0
for i_rank in range(size):
first_walker = last_walker
last_walker = last_walker + n_walkers_per_task
if last_walker > ns_args['n_walkers']:
last_walker = ns_args['n_walkers']
if i_rank == rank:
n_walkers = last_walker-first_walker
my_first_walker = first_walker
my_last_walker = last_walker
if last_walker > first_walker:
rank_of_walker[first_walker:last_walker] = [i_rank]*(last_walker-first_walker)
# figure out number of configs that will be culled on each task
n_cull = ns_args['n_cull']
n_extra_walk_per_task = ns_args['n_extra_walk_per_task']
max_n_cull_per_task = int(n_cull/size)
if max_n_cull_per_task * size != n_cull:
max_n_cull_per_task += 1
# internal model, LJ eps=1, sigma=1, cutoff=3, with PBC cube l=pbc[0,0]
internal_cutoff = 3.0
Eshift = internal_cutoff**-12 - internal_cutoff**-6
set_n_from_expected('n_model_calls')
if rank == 0:
print("Using n_model_calls = ", movement_args['n_model_calls'])
# create list of species, and check for possible problems
try:
species_list = ns_args['start_species'].split(',')
except:
species_list = []
if rank == 0:
init_atoms = ase.io.read(ns_args['start_config_file'], parallel=False)
atomic_numbers = init_atoms.get_atomic_numbers()
for Z in set(atomic_numbers):
n_of_Z = sum(atomic_numbers == Z)
if 'masses' in init_atoms.arrays:
mass_of_Z = init_atoms.get_masses()[np.where(atomic_numbers == Z)[0][0]]
species_list.append("%d %d %f" % (Z, n_of_Z, mass_of_Z))
else:
species_list.append("%d %d" % (Z, n_of_Z))
if comm is not None:
species_list = comm.bcast(species_list, root=0)
if do_calc_lammps:
if not ns_args['LAMMPS_atom_types'] == 'TYPE_EQUALS_Z':
used_chem_symbols = {ase.data.chemical_symbols[int(species.split()[0])] for species in species_list}
if not used_chem_symbols == set(ns_args['LAMMPS_atom_types'].keys()):
exit_error("species in start_species must correspond to those in LAMMPS_atom_types\n", 1)
mass_list = []
Z_list = []
warned_explicit_mass = False
for species in species_list:
species_fields = species.split()
Z_list.append(int(species_fields[0]))
if len(species_fields) == 3:
if not warned_explicit_mass:
if rank == 0:
sys.stderr.write("WARNING: setting masses explicitly. Not recommended, do only if you're sure it's necessary\n")
warned_explicit_mass = True
type_mass = float(species_fields[2])
mass_list.append(type_mass)
if len(mass_list) > 0:
mass_list = np.array(mass_list)
if np.any(mass_list != mass_list[0]) and not movement_args['atom_velo_rej_free_fully_randomize']:
exit_error("ERROR: Masses are not all equal, and atom_velo_rej_free_fully_randomize is false. Refusing to produce incorrect results\n", 1)
created_temp_restart_file = False
if ns_args['restart_file'] == "AUTO":
if rank == 0:
print("DOING restart_file=AUTO")
import glob
sfx = ns_args['config_file_format']
# try to match .0. or .ALL., but using sloppy regexp
print("checking snapshots glob", glob.iglob('%ssnapshot.[0-9]*.[0A]*.%s' % (ns_args['out_file_prefix'], sfx)))
try:
newest_snapshot = max(glob.iglob('%ssnapshot.[0-9]*.[0A]*.%s' % (ns_args['out_file_prefix'], sfx)), key=os.path.getmtime)
print("restarting from ", newest_snapshot)
if re.search('\.ALL\.%s$' % sfx, newest_snapshot) is not None: # latest snapshot is a combined one for all nodes
print("snapshot is combined for all nodes")
ns_args['restart_file'] = newest_snapshot
else:
snapshot_root = re.sub('\.0\.%s' % sfx, "", newest_snapshot)
restart_file = snapshot_root+".ALL."+sfx
print("creating combined snapshot file", restart_file, "from", [tf for tf in glob.iglob('%s.[0-9]*.%s' % (snapshot_root, sfx))])
created_temp_restart_file = True
with open(restart_file, "w") as snapshot_out:
for snapshot_file in glob.iglob('%s.[0-9]*.%s' % (snapshot_root, sfx)):
with open(snapshot_file, "r") as snapshot_in:
for line in snapshot_in:
snapshot_out.write(line)
ns_args['restart_file'] = restart_file
except:
print("no snapshot files found")
ns_args['restart_file'] = ''
if comm is not None:
ns_args['restart_file'] = comm.bcast(ns_args['restart_file'], root=0)
sys.stdout.flush()
# set up walkers
walkers = []
if ns_args['restart_file'] == '': # start from scratch
start_first_iter = 0
# create initial config
if rank == 0:
if ns_args['start_config_file'] is not None:
init_atoms = ase.io.read(ns_args['start_config_file'], parallel=False)
if not 'masses' in init_atoms.arrays:
init_atoms.set_masses([1.0] * len(init_atoms))
else:
# create atoms structs from a list of atomic numbers and numbers of atoms
# always create it slightly smaller than the max to avoid numerical instability with nearly identical volumes
if movement_args['2D']:
lc = 0.999*(ns_args['max_volume_per_atom']/ns_args['Z_cell_axis'])**(1.0/2.0)
init_atoms = ase.Atoms(cell=(lc, lc, ns_args['Z_cell_axis']), pbc=(1, 1, 1))
else:
lc = 0.999*ns_args['max_volume_per_atom']**(1.0/3.0)
init_atoms = ase.Atoms(cell=(lc, lc, lc), pbc=(1, 1, 1))
for species in species_list:
species_fields = species.split()
type_Z = int(species_fields[0])
type_n = int(species_fields[1])
if len(species_fields) == 2:
init_atoms += ase.Atoms([type_Z] * type_n, masses=[1.0] * type_n)
elif len(species_fields) == 3:
type_mass = float(species_fields[2])
init_atoms += ase.Atoms([type_Z] * type_n, masses=[type_mass] * type_n)
else:
exit_error("Each entry in start_species must include atomic number, multiplicity, and optionally mass", 5)
init_atoms.set_cell(init_atoms.get_cell()*float(len(init_atoms))**(1.0/3.0), scale_atoms=True)
ase.io.write(sys.stdout, init_atoms, parallel=False, format=ns_args['config_file_format'])
# ase.io.write(sys.stdout, init_atoms, format=ns_args['config_file_format'])
else: # rank != 0
init_atoms = None
# bcast atoms created on rank == 0
if comm is not None:
init_atoms = comm.bcast(init_atoms, root=0)
# create extra data arrays if needed
if ns_args['n_extra_data'] > 0:
init_atoms.arrays['ns_extra_data'] = np.zeros( (len(init_atoms), ns_args['n_extra_data']) )
# clone initial config into array of walkers
for i_walker in range(n_walkers):
walkers.append(init_atoms.copy())
# set up data structures to track configs as they are cloned and evolve
if ns_args['track_configs']:
if comm is None:
config_ind = 0
else:
config_ind = comm.rank*n_walkers
for at in walkers:
at.set_velocities(np.zeros((len(walkers[0]), 3)))
if ns_args['track_configs']:
at.info['config_ind'] = config_ind
at.info['from_config_ind'] = -1
at.info['config_ind_time'] = -1
config_ind += 1
if do_calc_ASE or do_calc_lammps:
at.calc = pot
if ns_args['track_configs']:
if comm is not None:
cur_config_ind = comm.size*n_walkers
else:
cur_config_ind = n_walkers
# V should have prob distrib p(V) = V^N.
# Using transformation rule p(y) = p(x) |dx/dy|, with p(y) = y^N and p(x) = 1,
# one gets dx/dy = y^N
# x = y^{N+1}
# y = x^{1/(N+1)}
if ns_args['start_energy_ceiling_per_atom'] is not None:
ns_args['start_energy_ceiling'] = ns_args['start_energy_ceiling_per_atom'] * len(init_atoms)
ns_args['start_energy_ceiling'] += movement_args['MC_cell_P']*ns_args['max_volume_per_atom']*len(init_atoms)
# initial positions are just random, up to an energy ceiling
for (i_at, at) in enumerate(walkers):
# randomize cell if P is set, both volume (from appropriate distribution) and shape (from uniform distribution with min aspect ratio limit)
if ns_args['random_initialise_cell']:
energy = float('nan')
if movement_args['MC_cell_P'] > 0.0:
if movement_args['2D']:
lc = (len(at)*ns_args['max_volume_per_atom']/ns_args['Z_cell_axis']*rng.float_uniform(0.0,1.0)**(1.0/float(len(at)+1)))**(1.0/2.0)
temp_cell = np.identity(3) * lc
temp_cell[2,2] = ns_args['Z_cell_axis']
at.set_cell( temp_cell )
do_cell_shape_walk(at, movement_args)
else:
lc = (len(at)*ns_args['max_volume_per_atom']*rng.float_uniform(0.0,1.0)**(1.0/float(len(at)+1)))**(1.0/3.0)
at.set_cell( np.identity(3) * lc )
do_cell_shape_walk(at, movement_args)
if ns_args['random_initialise_pos']:
# random initial positions
energy = float('nan')
n_try = 0
#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) ) )
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)
#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)
config_size = len(at)
min_bond = bonds[config_size]
#print('min_bond_here', min_bond)
reject_len = min_bond < ns_args['min_nn_dis']
else:
reject_len = False
#VGF end
energy = eval_energy(at)
n_try += 1
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
n_try = 0
if do_calc_fortran:
while n_try < ns_args['random_init_max_n_tries'] and (math.isnan(energy) or energy > ns_args['start_energy_ceiling']):
f_MC_MD.init_config(at, ns_args['start_energy_ceiling']-movement_args['MC_cell_P']*ns_args['max_volume_per_atom']*len(init_atoms))
energy = eval_energy(at)
n_try += 1
if math.isnan(energy) or energy > ns_args['start_energy_ceiling']:
sys.stderr.write("WARNING: rank %d failed to generate initial config by fortran config initializer under max energy %f in %d tries\n" % (rank, ns_args['start_energy_ceiling'], ns_args['random_init_max_n_tries']))
# try python config initializer
n_try = 0
while n_try < ns_args['random_init_max_n_tries'] and (math.isnan(energy) or energy > ns_args['start_energy_ceiling']):
energy = additive_init_config(at, ns_args['start_energy_ceiling'])
n_try += 1
# quit if failed to generate acceptable config
if math.isnan(energy) or energy > ns_args['start_energy_ceiling']:
exit_error("Rank %d failed to generate initial config by random, fortran, or (atom by atom addition) python initializer under max energy %f in %d tries each\n" % (rank, ns_args['start_energy_ceiling'], ns_args['random_init_max_n_tries']), 4)
energy = eval_energy(at)
at.info['ns_energy'] = rand_perturb_energy(energy, ns_args['random_energy_perturbation'])
at.info['volume'] = at.get_volume()
#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
# set KEmax from P and Vmax
if (movement_args['do_velocities']):
if ns_args['KEmax_max_T'] > 0.0:
KEmax = 1.5*len(walkers[0])*ns_args['kB']*ns_args['KEmax_max_T']
elif movement_args['MC_cell_P'] > 0.0:
KEmax = 1.5*movement_args['MC_cell_P']*len(walkers[0])*ns_args['max_volume_per_atom']
else:
exit_error("do_velocities is set, but neither KEmax_max_T nor MC_cell_P are > 0, so no heuristic for setting KEmax",4)
for at in walkers:
at.info['KEmax']=KEmax
else:
KEmax = -1.0
# set initial velocities, rejection free
if movement_args['do_velocities']:
for at in walkers:
rej_free_perturb_velo(at, None, KEmax) # adds KE to at.info['ns_energy']
# swap atomic numbers if doing semi-grand canonical ensemble
ns_args['swap_atomic_numbers'] = (movement_args['n_semi_grand_steps'] > 0)
else: # set up walkers with a restart
# swap atomic numbers if doing semi-grand canonical ensemble
ns_args['swap_atomic_numbers'] = (
movement_args['n_semi_grand_steps'] > 0)
at_list = ase.io.read(ns_args['restart_file'], index=":") # LIVIA
for r in range(size):
if rank == r:
walkers = at_list[r*n_walkers:(r+1)*n_walkers] # TODO: RBW – split walkers on different processes? maybe we need to set things up (energies?) before splitting?
print(rank, r, walkers)
for at in walkers:
if np.any(at.get_atomic_numbers()
!= walkers[0].get_atomic_numbers()):
ns_args['swap_atomic_numbers'] = True
# broadcast swap_atomic_numbers in case it was overridden to True
# by presence of configurations with different atomic number lists
if comm is not None:
ns_args['swap_atomic_numbers'] = comm.bcast(
ns_args['swap_atomic_numbers'], root=0)
if ns_args['track_configs']:
if comm is not None:
cur_config_ind = comm.size*n_walkers
else:
cur_config_ind = n_walkers
if 'config_ind' in walkers[0].info:
max_config_ind = max(
[at.info['config_ind'] for at in walkers])
if comm is not None:
comm.allreduce(max_config_ind, op=MPI.MAX)
cur_config_ind = max_config_ind+1
# Calculate ns_energy for manual surface configurations from restart
if movement_args['keep_atoms_fixed'] > 0:
# print("RBW: calc ns_energy for man surf configs from restart") # debug
for (i_at, at) in enumerate(walkers):
if do_calc_ASE or do_calc_lammps:
at.set_calculator(pot)
energy = eval_energy(at)
at.info['ns_energy'] = rand_perturb_energy(
energy, ns_args['random_energy_perturbation'])
at.info['volume'] = at.get_volume()
if movement_args['do_velocities']:
KEmax = walkers[0].info['KEmax']
else:
KEmax = -1.0
# Set init velocities for manual surface configurations from restart
# RBW checked this and seed-averaged initial total energies are ~
# the same, which means velocities *seem* to be set properly
# LBP: only perturb if this is a "fake" restart to start a surface sim.
# proper restart does not need perturbation as that messes up KE
if movement_args['keep_atoms_fixed'] > 0 and walkers[0].info['iter'] < 1:
# print("RBW: set init vels for man surf configs from restart") # debug
if movement_args['do_velocities']:
for at in walkers:
# TODO: RBW – make sure surface restart files contain KEmax calculated from number of free atoms, otherwise, their kinetic energy can be inefficiently high
energy = eval_energy(at)
rej_free_perturb_velo(at, None, KEmax)
# adds KE to at.info['ns_energy']
energy = eval_energy(at)
for at in walkers:
if ns_args['n_extra_data'] > 0 and (not 'ns_extra_data' in at.arrays or at.arrays['ns_extra_data'].size/len(at) != ns_args['n_extra_data']):
at.arrays['ns_extra_data'] = np.zeros( (len(at), ns_args['n_extra_data']) )
if do_calc_ASE or do_calc_lammps:
at.calc = pot
at.info['ns_energy'] = rand_perturb_energy(eval_energy(at), ns_args['random_energy_perturbation'])
if 'iter' in at.info:
start_first_iter = at.info['iter'] + 1
else:
print("ERROR: no iteration number information was found "
"in the restart file")
exit_error("no iteration number information was found in "
"the restart file\n", 5)
key_found = False
for key in at.info:
# check if 'volume=' info is present in the file used for
# restart
if key == 'volume':
movement_args['MC_cell_volume_per_atom_step_size'] = (
at.info['volume']/10.0/len(at))
key_found = True
if not key_found:
print( "WARNING: no volume information was found in the restart file. If volume changes will be done, the starting stepsize will be the default")
sys.stdout.flush()
# add GMC direction if needed
movement_args['do_GMC'] = ((movement_args['atom_algorithm'] == 'MC') and (movement_args['MC_atom_Galilean']))
if movement_args['do_GMC']:
for at in walkers:
if 'GMC_direction' not in at.arrays:
at.arrays['GMC_direction'] = np.zeros((len(at), 3))
# scale initial MC_atom_step_size by max_vol^(1/3)
max_lc = (ns_args['max_volume_per_atom']*len(walkers[0]))**(1.0/3.0)
movement_args['MC_atom_step_size'] *= max_lc
movement_args['MC_atom_step_size_max'] *= max_lc
# scale MC_cell_shear_step_size by max_vol^1.0)
movement_args['MC_cell_shear_step_size'] *= max_lc
movement_args['MC_cell_shear_step_size_max'] *= max_lc
if ns_analyzers is not None:
for (ns_analyzer, ns_analyzer_interval) in ns_analyzers:
if ns_analyzer_interval < 0:
ns_analyzer.analyze(walkers, -1, "initial_walk start")
# do initial walks if needed
if ns_args['initial_walk_N_walks'] > 0 and ns_args['restart_file'] == '':
if rank == 0:
print("doing initial_walks", ns_args['initial_walk_N_walks'])
t0 = time.time()
(Emax, Vmax, cull_rank, cull_ind) = max_energy(walkers, 1)
# WARNING: this assumes that all walkers have same numbers of atoms
Emax = Emax[0] + ns_args['initial_walk_Emax_offset_per_atom']*len(walkers[0])
# do full walks, not shortened walks that account for NS algorithm
# parallelization
save_n_model_calls = movement_args['n_model_calls']
movement_args['n_model_calls'] = movement_args['n_model_calls_expected']
walk_stats_adjust = {}
zero_stats(walk_stats_adjust, movement_args)
for i_initial_walk in range(ns_args['initial_walk_N_walks']):
if rank == 0:
print("initial walk start iter ", i_initial_walk, "time",
time.time()-t0)
print_prefix = "%d initial_walk %d" % (rank, i_initial_walk)
if i_initial_walk > 0 and (i_initial_walk-1) % ns_args['initial_walk_adjust_interval'] == 0: # first adjust is after first walk
# could be done before first walk if full_auto_set_stepsize
# didn't need walk_stats
full_auto_set_stepsizes(walkers, walk_stats_adjust,
movement_args, comm, Emax, -1, size)
walk_stats_adjust = {}
zero_stats(walk_stats_adjust, movement_args)
for (i_at, at) in enumerate(walkers):
print_prefix = "%d initial_walk %d at %d" % (rank, i_initial_walk, i_at)
if ns_args['profile'] == rank:
import cProfile
pr = cProfile.Profile()
walk_stats = pr.runcall(walk_single_walker, at=at, movement_args=movement_args, Emax=Emax, KEmax=-1)
pr.dump_stats(ns_args['out_file_prefix']+'initial_walk.profile.stats')
else:
walk_stats = walk_single_walker(at, movement_args, Emax, -1)
accumulate_stats(walk_stats_adjust, walk_stats)
if ns_analyzers is not None:
for (ns_analyzer, ns_analyzer_interval) in ns_analyzers:
if ns_analyzer_interval < 0 and (ns_analyzer_interval == -1 or (i_initial_walk+1) % (-ns_analyzer_interval) == 0):
ns_analyzer.analyze(walkers, -1, "initial_walk %d" % i_initial_walk)
if ns_args['snapshot_interval'] > 0 and (i_initial_walk+1) % ns_args['snapshot_interval'] == 0:
save_snapshot(i_initial_walk-ns_args['initial_walk_N_walks'])
# restore walk lengths for rest of NS run
movement_args['n_model_calls'] = save_n_model_calls
(KEmax, _, _, _) = max_energy(walkers, 1, kinetic_only=True)
if ns_args['initial_walk_only']:
if comm is not None:
MPI.Finalize()
sys.exit(0)
sys.stdout.flush()
n_atoms = len(walkers[0])
# do NS
# open the file where the trajectory will be printed
if ns_args['restart_file'] == '': # start from scratch, so if this file exists, overwrite it
traj_io = open(ns_args['out_file_prefix']+'traj.%d.%s' % (rank, ns_args['config_file_format']), "w")
if ns_args['track_configs'] and ns_args['track_configs_write']:
track_traj_io = open(ns_args['out_file_prefix']+'track_traj.%d.%s' % (rank, ns_args['config_file_format']), "w")
else:
track_traj_io = None
else: # restart, so the existing file should be appended
# concatenate existing traj file to before restart
print(rank, "truncating traj file to start_first_iter",
start_first_iter)
mode = "r+"
if movement_args['keep_atoms_fixed'] > 0:
mode = "a+"
with open(ns_args['out_file_prefix']+'traj.%d.%s' % (rank, ns_args['config_file_format']), mode) as f:
prev_pos = None
# loop this way with "while True" and "f.readline()" because directly looping over f does not
# set position reported by f.tell() to end of line
# make sure there's somplace to truncate to if traj file has only one config and it's already too late
prev_pos = 0
prev_prev_pos = 0
while True:
l = f.readline()
if not l:
break
m = re.search(r"\biter=(\d+)\b", l)
if m is not None:
cur_iter = int(m.group(1))
if cur_iter >= start_first_iter:
# rewind back to two lines back, before start of this config
f.seek(prev_prev_pos)
f.truncate()
break
prev_prev_pos = prev_pos
prev_pos = f.tell()
traj_io = open(ns_args['out_file_prefix']+'traj.%d.%s' % (rank, ns_args['config_file_format']), "a")
if ns_args['track_configs'] and ns_args['track_configs_write']:
track_traj_io = open(ns_args['out_file_prefix']+'track_traj.%d.%s' % (rank, ns_args['config_file_format']), "a")
else:
track_traj_io = None
# Read the existing traj file and look for the point where we restart from. Truncate the rest.
# This part is not used because the ASE.io.read takes soooo long, that it makes a restart impossible.
#traj_io = open(ns_args['out_file_prefix']+'traj.%d.%s' % (rank, ns_args['config_file_format']), "r+")
#i = 0
#while True:
# at=(ase.io.read(traj_io, format=ns_args['config_file_format'],index=i))
# print("ASE.io.read trajectory", rank, i, at.info['iter'])
# if at.info['iter'] >= start_first_iter:
# at=(ase.io.read(traj_io, format=ns_args['config_file_format'],index=i-1))
# traj_io.truncate()
# break
# i += 1
sys.stdout.flush()
if ns_args['E_dump_interval'] > 0 and rank == 0:
E_dump_io = open(ns_args['out_file_prefix']+'E_dump', "w")
E_dump_io.write("n_walkers %d\n" % ns_args['n_walkers'])
else:
E_dump_io = None
# open the file where the energies will be printed
if rank == 0:
if ns_args['restart_file'] == '': # start from scratch, so if this file exists, overwrite it
energy_io = open(ns_args['out_file_prefix']+'energies', 'w')
else: # restart, so the existing file should be appended
try:
energy_io = open(ns_args['out_file_prefix']+'energies', 'r+')
line_pos = 0
tmp_iter = 0
line = energy_io.readline() # read the first line of nwalker,ncull..etc information
i = 0
line_pos += len(line)
while True: # we do create an infinite loop here :(
line = energy_io.readline() # read lines one by one
if not line: # something went wrong, exit the infinit loop
print("WARNING: end of .energies file reached without finding the iteration number", start_first_iter)
break
i = i + 1
line_pos += len(line) #Track the position of the end of each line
if i % 10000 == 0:
print(rank, "reading .energies file line %d" % i)
if i % n_cull == 0: # if this is n_cull-th line, examine the stored iteration
tmp_split = line.split()
tmp_iter = int(tmp_split[0]) # tmp_iter contains the iteration number of the line as an integer number
if tmp_iter == start_first_iter - 1: # if this is the iteration same as in the snapshot,
print(rank, "truncating energy file at line ", i)
energy_io.seek(line_pos, 0) #Move the file pointer to the end of the line
energy_io.truncate() #delete the rest of the file, as we are restarting from here
break
except FileNotFoundError:
print("WARNING: got restart file, but no corresponding energies file, so creating new one from scratch")
energy_io = open(ns_args['out_file_prefix']+'energies', 'w')
sys.stdout.flush()
if ns_args['profile'] == rank:
print("if")
import cProfile
pr = cProfile.Profile()
final_iter = pr.runcall(do_ns_loop)
pr.dump_stats(ns_args['out_file_prefix']+'profile.stats')
else:
print("else")
final_iter = do_ns_loop()
# cleanup post loop
save_snapshot(final_iter) # this is the final configuration
for at in walkers:
print(rank, ": final energy ", at.info['ns_energy'])
if rank == 0:
energy_io.close()
traj_io.close()
if track_traj_io is not None:
track_traj_io.close()
if E_dump_io is not None:
E_dump_io.close()
if comm is not None:
MPI.Finalize()
sys.exit(0)