diff --git a/nexus/nexus/rmg_analyzer.py b/nexus/nexus/rmg_analyzer.py index f57494ac58..6f79db1183 100644 --- a/nexus/nexus/rmg_analyzer.py +++ b/nexus/nexus/rmg_analyzer.py @@ -11,6 +11,8 @@ import numpy as np from .developer import DevBase, dotdict, obj +from .generic import NexusError +from .rmg_input import RmgInput from .simulation import Simulation, SimulationAnalyzer from .structure import generate_structure from .unit_converter import UnitConverter, convert @@ -32,6 +34,17 @@ def normalize_line(line): #end def normalize_line +def line_numbers(line): + """Return all finite, whitespace-delimited RMG numbers in a line.""" + values = [] + for token in line.replace(',',' ').split(): + value = as_float(token.strip('()[]{}')) + if value is not None: + values.append(value) + return np.array(values,dtype=float) +#end def line_numbers + + class RmgOutData(DevBase): """Read an RMG output file and collect results appropriate to its run mode. @@ -48,27 +61,64 @@ class RmgOutData(DevBase): Absolute path to the output directory. outfile_name : str Name of the RMG output file. + input : RmgInput or None + Parsed control input when the referenced input file is available. setup_info : obj - Run mode and the initial structure, when available. + Parsed setup sections, run mode, structure, and k-point information. run_mode : str or None Short RMG calculation mode: ``"scf"``, ``"nscf"``, or ``"relax"``. geometry : obj or None - Cartesian k-points and k-point weights. + Cell volume and crystal/Cartesian k-point information. + convergence : obj or None + Electronic and ionic convergence indicators and event counts. + timing : obj or None + Total, per-step, and section-resolved timing data in seconds. energy : float or numpy.floating or None Last total energy obtained from the eigenvalue sum. energy_units : str or None Units associated with ``energy``. + energies : numpy.ndarray or None + History of total energies obtained from eigenvalue sums. + energy_units_history : numpy.ndarray or None + Units corresponding to ``energies``. + direct_energies : numpy.ndarray or None + History of directly evaluated total energies. + direct_energy_units : numpy.ndarray or None + Units corresponding to ``direct_energies``. electronic : obj or None Fermi energies, band edges, gaps, k-points, eigenvalues, and - occupations when reported. + occupations, charge, magnetization, force, volume, and per-atom energy + data when reported. + scf : obj or None + SCF energy components, iteration indices, residuals, and timing data. + ionic_steps : obj or None + Detailed per-step ionic records. + position_units : str or None + Units associated with ionic positions. + force_units : str or None + Units associated with ionic forces. + positions : numpy.ndarray or None + Ionic positions with shape ``(nsteps, natoms, 3)``. forces : numpy.ndarray or None Ionic forces with shape ``(nsteps, natoms, 3)``. + charges : numpy.ndarray or None + Ionic charges with shape ``(nsteps, natoms)``. + magnetizations : numpy.ndarray or None + Ionic magnetizations with shape ``(nsteps, natoms)``. + max_forces : numpy.ndarray or None + Maximum ionic force magnitude at each ionic step. structures : obj or None Mapping from ionic-step index to a :class:`Structure` instance. stress : numpy.ndarray or None Stress tensors with shape ``(nsteps, 3, 3)``. + stress_units : str or None + Units associated with stress and pressure values. + pressures : numpy.ndarray or None + Hydrostatic pressure at each reported stress step. pressure : float or numpy.floating or None Last hydrostatic pressure. + produced_files : obj or None + Paths to recognized files produced by an SCF run. Notes ----- @@ -86,6 +136,9 @@ class RmgOutData(DevBase): If ``filepath`` does not identify a regular file. """ + # This pattern represents RMG numbers embedded in records whose structural + # punctuation must also be recognized. Splitting those records on whitespace + # would not reliably separate signed values, brackets, and optional exponents. # Match a signed integer or decimal with an optional E- or D-exponent. # Example: -1.2345D+02 number_pattern = r'[-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[EeDd][-+]?\d+)?' @@ -113,11 +166,12 @@ def __init__(self,filepath): f'Path provided: {filepath}' ) raise IsADirectoryError(msg) - path,outfile_name = os.path.split(filepath) - self.path = path - self.abspath = os.path.abspath(path) - self.outfile_name = outfile_name - self.setup_info = None + path,outfile_name = os.path.split(filepath) + self.path = path + self.abspath = os.path.abspath(path) + self.outfile_name = outfile_name + self.input = None + self.setup_info = None with open(filepath,'r') as input_file: lines = input_file.read().splitlines() @@ -125,30 +179,157 @@ def __init__(self,filepath): # modes: scf, nscf, relax if self.run_mode in {'scf','nscf','relax'}: - self.geometry = None - self.energy = None - self.energy_units = None - self.electronic = None - self.forces = None - self.structures = None - self.stress = None - self.pressure = None + self.geometry = None + self.convergence = None + self.timing = None + self.energy = None + self.energy_units = None + self.energies = None + self.energy_units_history = None + self.direct_energies = None + self.direct_energy_units = None + self.electronic = None + self.scf = None + self.ionic_steps = None + self.position_units = None + self.force_units = None + self.positions = None + self.forces = None + self.charges = None + self.magnetizations = None + self.max_forces = None + self.structures = None + self.stress = None + self.stress_units = None + self.pressures = None + self.pressure = None self.read_geometry() + self.read_convergence(lines) + self.read_timing(lines) self.read_energies(lines) + self.read_scf(lines) self.read_ions(lines) self.read_stress(lines) self.read_electronic(lines) + + # modes: scf + if self.run_mode=='scf': + self.produced_files = None + self.read_produced_files() #end def __init__ def read_setup_info(self,lines): - """Read the run mode and initial structure from the setup report.""" + """Read setup sections, the run mode, and the initial structure. + + Binds ``setup_info`` to an ``obj`` containing normalized setup blocks + and derived grid, lattice, ion, k-point, and structure data. When the + referenced control file is available, also binds ``input``. + + Parameters + ---------- + lines : list of str + Complete RMG log split into lines. + """ setup_info = obj( run_mode = None, structure = None, k_points = None, + files = None, ) + position_heading = 'initial ionic positions and displacements' + + def process_name(text): + """Convert an RMG setup label to a normalized member name.""" + # Parenthetical annotations can occur inside a label with meaningful + # text on either side. Delimiter splitting would discard one side, + # while tokenization cannot identify the full annotation reliably. + # Remove parenthetical annotations from setup labels. + # Example: Grid Points (Linear Anisotropy: 1.000) + text = re.sub(r'\([^)]*\)','',text) + name = '_'.join(text.strip().lower().split()) + return name.replace('/','_').replace('-','_') + #end def process_name + + # Parse the indented setup report into named persistent sections. + sections = obj() + current = None + files = None + grid_points = None + lattice_setup = None + in_setup = False + section_added = False + for raw_line in lines: + stripped = raw_line.strip() + if not in_setup: + if stripped.lower()!='files': + continue + in_setup = True + if normalize_line(stripped).lower().startswith(position_heading): + break + if len(stripped)==0: + continue + if not raw_line[0].isspace(): + section_name = process_name(stripped.rstrip(':')) + section_added = False + if section_name=='files': + current = obj( + control_input_file = None, + data_output_file = None, + ) + files = current + elif section_name=='grid_points': + current = obj( + equivalent_energy_cutoffs = None, + units = None, + ) + grid_points = current + elif section_name=='lattice_setup': + current = obj() + lattice_setup = current + else: + current = obj() + continue + if current is None or ':' not in stripped: + continue + if not section_added and section_name!='k_points': + sections[section_name] = current + section_added = True + label,value = stripped.split(':',1) + name = process_name(label) + value = value.strip() + units = None + + # Convert simple Boolean, integer, and floating-point fields. + upper_value = value.upper() + number = as_float(value) + if upper_value in {'ON','OFF'}: + value = upper_value=='ON' + elif number is not None: + is_integer = ( + number.is_integer() + and not any(c in value.lower() for c in '.e') + ) + value = int(number) if is_integer else number + else: + # Convert numeric sequences, separating a trailing unit label. + tokens = value.replace(',',' ').split() + if len(tokens)>1 and as_float(tokens[-1]) is None: + numeric = [as_float(token) for token in tokens[:-1]] + if None not in numeric: + units = tokens[-1] + tokens = tokens[:-1] + values = [as_float(token) for token in tokens] + if ( + (len(values)>1 or units is not None and len(values)>0) + and None not in values + ): + value = np.array(values,dtype=float) + current[name] = value + if units is not None: + current.units = units + setup_info.update(sections) # Read the calculation mode from the run setup block. run_mode = None @@ -199,7 +380,6 @@ def read_setup_info(self,lines): if len(tokens)>=4: axis_unit = tokens[3].strip(',;') - position_heading = 'initial ionic positions and displacements' position_tables = [] i = 0 # Read each initial-position table, retaining its reported units. @@ -220,8 +400,8 @@ def read_setup_info(self,lines): line = lines[i] if len(line.strip())==0 and len(atoms)>0: break - tokens = line.split() - atom = tokens[0] if len(tokens)>0 else '' + tokens = line.split() + atom = tokens[0] if len(tokens)>0 else '' valid_atom = ( len(atom)>0 and atom[0].isalpha() @@ -237,21 +417,64 @@ def read_setup_info(self,lines): positions.append(values) i += 1 if len(atoms)>0: - position_tables.append(obj( - units = units, - atoms = np.array(atoms,dtype=object), - positions = np.array(positions,dtype=float), - )) + position_tables.append( + obj( + units = units, + atoms = np.array(atoms,dtype=object), + positions = np.array(positions,dtype=float), + ), + ) + + if grid_points is not None: + grid = [] + grid_pe = [] + grid_spacing = [] + for direction in ('x','y','z'): + if direction not in grid_points: + break + text = normalize_line(str(grid_points[direction])) + lower = text.lower() + values = [] + for label in ('total','per pe','spacing'): + start = lower.find(label) + if start<0: + break + remainder = text[start+len(label):].lstrip(' :') + tokens = remainder.split() + value = as_float(tokens[0]) if len(tokens)>0 else None + if value is None: + break + values.append(value) + if len(values)!=3: + break + grid.append(values[0]) + grid_pe.append(values[1]) + grid_spacing.append(values[2]) + if len(grid)==3: + grid_points.grid = np.array(grid,dtype=int) + grid_points.grid_pe = np.array(grid_pe,dtype=int) + grid_points.grid_spacing = np.array(grid_spacing,dtype=float) + grid_points.grid_units = 'a0' + cutoffs = grid_points.equivalent_energy_cutoffs + if cutoffs is not None: + cutoff_values = line_numbers(str(cutoffs)) + if len(cutoff_values)>=2: + grid_points.ecut = cutoff_values[0] + grid_points.ecut_charge = cutoff_values[1] + grid_points.ecut_units = grid_points.units if set(axes)=={'x','y','z'} and len(position_tables)>0: ion_positions = next( (table for table in position_tables if table.units=='B'), position_tables[0], ) - aunits = 'B' if axis_unit in {None,'a0','B','bohr'} else 'A' - axes_array = np.array( - [axes[c] for c in ('x','y','z')],dtype=float) - axes_array = convert(axes_array,aunits,'B') + setup_info.ion_positions = ion_positions + aunits = 'B' if axis_unit in {None,'a0','B','bohr'} else 'A' + reported_axes = np.array( + [axes[c] for c in ('x','y','z')], + dtype=float, + ) + axes_array = convert(reported_axes,aunits,'B') positions = convert(ion_positions.positions,ion_positions.units,'B') valid = ( axes_array.shape==(3,3) @@ -266,6 +489,8 @@ def read_setup_info(self,lines): elem = ion_positions.atoms, pos = positions, ) + if lattice_setup is not None: + lattice_setup.axes = reported_axes kpoints = [] kweights = [] @@ -300,19 +525,57 @@ def read_setup_info(self,lines): setup_info.structure.add_kpoints( kpoints, kweights, - recenter = False, + recenter = False, cell_unit = True, ) + + if files is not None and files.control_input_file is not None: + control_file = str(files.control_input_file) + filepaths = ( + os.path.join(self.path,control_file), + os.path.join(self.path,os.path.basename(control_file)), + os.path.join(os.path.dirname(self.path),control_file), + ) + filepath = next( + (path for path in filepaths if os.path.isfile(path)),None) + if filepath is not None: + try: + self.input = RmgInput(filepath) + except (NexusError,OSError,TypeError,ValueError): + pass self.setup_info = setup_info #end def read_setup_info def read_energies(self,lines): - """Read the final total energy obtained from the eigenvalue sum.""" - label = 'final total energy from eig sum' + """Read eigenvalue-sum and direct total-energy histories. + + Binds the history arrays and their unit arrays, as well as the final + eigenvalue-sum energy and its units. + + Parameters + ---------- + lines : list of str + Complete RMG log split into lines. + """ + energies = [] + energy_units = [] + direct_energies = [] + direct_energy_units = [] for line in lines: - text = normalize_line(line) - lower = text.lower() - if label not in lower: + text = normalize_line(line) + lower = text.lower() + label = None + target_values = None + target_units = None + if 'final total energy from eig sum' in lower: + label = 'final total energy from eig sum' + target_values = energies + target_units = energy_units + elif 'final total energy from direct' in lower: + label = 'final total energy from direct' + target_values = direct_energies + target_units = direct_energy_units + if label is None: continue remainder = text[lower.index(label)+len(label):].lstrip() if len(remainder)==0 or remainder[0] not in {':','='}: @@ -323,8 +586,21 @@ def read_energies(self,lines): value = as_float(tokens[0]) if value is None: continue - self.energy = value - self.energy_units = tokens[1].strip(',;') if len(tokens)>=2 else 'Ha' + target_values.append(value) + target_units.append( + tokens[1].strip(',;') if len(tokens)>=2 else None, + ) + if len(energies)>0: + self.energies = np.array(energies,dtype=float) + self.energy_units_history = np.array(energy_units,dtype=object) + self.energy = self.energies[-1] + self.energy_units = self.energy_units_history[-1] or 'Ha' + if len(direct_energies)>0: + self.direct_energies = np.array(direct_energies,dtype=float) + self.direct_energy_units = np.array( + direct_energy_units, + dtype=object, + ) #end def read_energies @@ -332,7 +608,9 @@ def read_electronic(self,lines): """Parse electronic quantities exposed by ``RmgAnalyzer``. Binds ``electronic`` to an ``obj`` containing Fermi energies, band - edges, gaps, k-point-major eigenvalues, occupations, and k-points. + edges, gaps, charge and magnetization values, summed forces, per-atom + volume and energy, k-point-major eigenvalues and occupations, and + k-points. Parameters ---------- @@ -358,22 +636,31 @@ def assigned_value(text,lower,*labels): #end def assigned_value data = obj( - fermi_energies = [], - valence_band_maxima = [], - conduction_band_minima = [], - band_gaps = [], - kpoints_crystal = None, - kpoints = None, - eigenvalues = None, - occupations = None, + fermi_energies = [], + valence_band_maxima = [], + conduction_band_minima = [], + band_gaps = [], + total_charges = [], + total_magnetizations = [], + absolute_magnetizations = [], + sum_forces = [], + volume_per_atom = [], + energy_per_atom = [], + kpoints_crystal = None, + kpoints = None, + eigenvalues = None, + occupations = None, ) + # Each row can contain several value/occupation pairs with arbitrary + # whitespace inside and outside the brackets. Plain token positions are + # therefore unstable and cannot preserve the pair boundaries safely. # Match one eigenvalue followed by its bracketed occupation. # Example: -6.4238 [2.000] - npat = self.number_pattern + npat = self.number_pattern pair_pattern = re.compile( r'('+npat+r')\s*\[\s*('+npat+r')\s*\]', - re.IGNORECASE + re.IGNORECASE, ) datasets = [] dataset = dotdict() @@ -391,7 +678,22 @@ def assigned_value(text,lower,*labels): 'conduction band minimum', 'conduction band minumm', ) - gap = assigned_value(text,lower,'band gap') + gap = assigned_value(text,lower,'band gap') + total_charge = assigned_value( + text, + lower, + 'total charge in supercell', + ) + total_magnetization = assigned_value( + text, + lower, + 'total magnetization', + ) + absolute_magnetization = assigned_value( + text, + lower, + 'absolute magnetization', + ) if fermi is not None: data.fermi_energies.append(fermi) elif vbm is not None and cbm is not None: @@ -399,10 +701,25 @@ def assigned_value(text,lower,*labels): data.conduction_band_minima.append(cbm) elif gap is not None: data.band_gaps.append(gap) + elif total_charge is not None: + data.total_charges.append(total_charge) + elif total_magnetization is not None: + data.total_magnetizations.append(total_magnetization) + elif absolute_magnetization is not None: + data.absolute_magnetizations.append(absolute_magnetization) + elif lower.startswith('sum force'): + values = line_numbers(text.partition('=')[2]) + if len(values)>=3: + data.sum_forces.append(values[:3]) + elif 'volume and energy per atom' in lower: + values = line_numbers(text.partition('=')[2]) + if len(values)>=2: + data.volume_per_atom.append(values[0]) + data.energy_per_atom.append(values[1]) if 'kohn sham eigenvalues' in lower and 'k-point' in lower: kpoint_start = lower.rfind('k-point')+len('k-point') - kpoint_text = text[kpoint_start:] + kpoint_text = text[kpoint_start:] index_text,separator,coordinates_text = kpoint_text.partition(']') if len(separator)==0 or '[' not in index_text: continue @@ -411,7 +728,7 @@ def assigned_value(text,lower,*labels): except ValueError: continue coordinate_tokens = coordinates_text.lstrip(' :').split() - coordinates = [as_float(v) for v in coordinate_tokens[:3]] + coordinates = [as_float(v) for v in coordinate_tokens[:3]] if len(coordinates)!=3 or None in coordinates: continue if index in dataset: @@ -453,6 +770,10 @@ def assigned_value(text,lower,*labels): for name,values in data.items(): if values is not None: data[name] = np.array(values,dtype=float) + data.energy_units = 'eV' + data.magnetization_units = 'Bohr mag/cell' + data.sum_force_units = 'Ha/a0' + data.energy_per_atom_units = 'eV' if len(dataset)>0: datasets.append(dataset) @@ -478,7 +799,7 @@ def assigned_value(text,lower,*labels): channels = [ candidate[index].channels.get(spin) for index in indices for spin in spins - ] + ] if any( channel is None or len(channel[0])==0 @@ -489,21 +810,27 @@ def assigned_value(text,lower,*labels): if len({len(channel[0]) for channel in channels})!=1: continue data.kpoints_crystal = np.array( - [candidate[i].kpoint for i in indices],dtype=float) - data.eigenvalues = np.array([ - [candidate[i].channels[spin][0] for spin in spins] - for i in indices - ],dtype=float) - data.occupations = np.array([ - [candidate[i].channels[spin][1] for spin in spins] - for i in indices - ],dtype=float) + [candidate[i].kpoint for i in indices], + dtype=float, + ) + data.eigenvalues = np.array( + [[candidate[i].channels[spin][0] for spin in spins] + for i in indices], + dtype=float, + ) + data.occupations = np.array( + [[candidate[i].channels[spin][1] for spin in spins] + for i in indices], + dtype=float, + ) if spins==['none']: data.eigenvalues = data.eigenvalues[:,0,:] data.occupations = data.occupations[:,0,:] if self.setup_info.structure is not None: data.kpoints = np.dot( - data.kpoints_crystal,self.setup_info.structure.kaxes) + data.kpoints_crystal, + self.setup_info.structure.kaxes, + ) break nfound = sum(v.size for v in data.values() if isinstance(v,np.ndarray)) @@ -512,22 +839,130 @@ def assigned_value(text,lower,*labels): #end def read_electronic + def read_scf(self,lines): + """Read SCF energy components, iteration indices, residuals, and times. + + Binds ``scf`` to an ``obj`` containing NumPy histories. Energies are + in Hartree and times are in seconds. + + Parameters + ---------- + lines : list of str + Complete RMG log split into lines. + """ + component_names = { + 'eigenvalue sum' : 'eigenvalue_sum', + 'ion_ion' : 'ion_ion', + 'electrostatic' : 'electrostatic', + 'vxc' : 'vxc', + 'exc' : 'exc', + 'total energy' : 'total_energy', + 'estimated error' : 'estimated_error', + } + values = dotdict( + eigenvalue_sum = [], + ion_ion = [], + electrostatic = [], + vxc = [], + exc = [], + total_energy = [], + estimated_error = [], + ) + + # Detailed summaries contain an optional subset of multiword fields in + # varying order, and some labels include brackets. Positional splitting + # would couple parsing to the current order and fail when fields are absent. + # Match fields within a detailed SCF-iteration summary. + # Example: quench: [md: 0/2 scf: 3/20 step time: 0.20 RMS[dV]: 2e-5] + detail_pattern = re.compile( + r'\bmd\s*:\s*(?P\d+)\s*/|' + r'\bscf\s*:\s*(?P\d+)\s*/|' + r'\bstep\s+time\s*:\s*(?P'+self.number_pattern+r')|' + r'\bscf\s+time\s*:\s*(?P