From bda47a25a295b42b189045617a3da7fde6184a55 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Tue, 3 Sep 2019 13:32:18 -0400 Subject: [PATCH] remove untested modules --- flare/modules/__init__.py | 0 flare/modules/__init__.pyc | Bin 121 -> 0 bytes flare/modules/ab_initio_activation.py | 77 ---- flare/modules/activation_parser.py | 51 --- flare/modules/analyze_gp.py | 160 ------- flare/modules/analyze_md.py | 100 ----- flare/modules/analyze_otf.py | 212 ---------- flare/modules/ase_calculator.py | 109 ----- flare/modules/ase_otf.py | 213 ---------- flare/modules/ase_otf_logger.py | 105 ----- flare/modules/ase_otf_md.py | 61 --- flare/modules/atomeye_helper.py | 297 ------------- flare/modules/calculate_rdf.py | 89 ---- flare/modules/convergence.py | 177 -------- flare/modules/crystals.py | 125 ------ flare/modules/cutoff_sweep.py | 143 ------- flare/modules/data_sweep.py | 155 ------- flare/modules/flare_parsers/__init__.py | 0 flare/modules/flare_parsers/add_all_parser.py | 397 ------------------ .../flare_parsers/gp_from_aimd_parser.py | 296 ------------- flare/modules/flare_parsers/old_vac_parser.py | 218 ---------- flare/modules/flare_parsers/otf_parser_v0.py | 267 ------------ .../flare_parsers/parse_output_temp.py | 106 ----- flare/modules/gp_calculator.py | 59 --- flare/modules/gp_from_aimd.py | 392 ----------------- flare/modules/gp_vs_deepmd.py | 102 ----- flare/modules/hyp_opt.py | 267 ------------ flare/modules/initialize_velocities.py | 40 -- flare/modules/lammps/__init__.py | 0 flare/modules/lammps/coulomb_util.py | 132 ------ flare/modules/lammps/eam.py | 129 ------ flare/modules/lammps/mc.py | 171 -------- flare/modules/qe_input.py | 277 ------------ flare/modules/qe_parsers.py | 147 ------- flare/modules/rcut.py | 153 ------- flare/modules/test_otf.py | 78 ---- flare/modules/vacancy_diffusion.py | 97 ----- flare/modules/validate.py | 200 --------- .../{modules/flare_parsers => }/otf_parser.py | 2 +- tests/test_parse_otf.py | 2 +- 40 files changed, 2 insertions(+), 5604 deletions(-) delete mode 100644 flare/modules/__init__.py delete mode 100644 flare/modules/__init__.pyc delete mode 100644 flare/modules/ab_initio_activation.py delete mode 100644 flare/modules/activation_parser.py delete mode 100644 flare/modules/analyze_gp.py delete mode 100644 flare/modules/analyze_md.py delete mode 100644 flare/modules/analyze_otf.py delete mode 100644 flare/modules/ase_calculator.py delete mode 100644 flare/modules/ase_otf.py delete mode 100644 flare/modules/ase_otf_logger.py delete mode 100644 flare/modules/ase_otf_md.py delete mode 100644 flare/modules/atomeye_helper.py delete mode 100644 flare/modules/calculate_rdf.py delete mode 100644 flare/modules/convergence.py delete mode 100644 flare/modules/crystals.py delete mode 100644 flare/modules/cutoff_sweep.py delete mode 100644 flare/modules/data_sweep.py delete mode 100644 flare/modules/flare_parsers/__init__.py delete mode 100644 flare/modules/flare_parsers/add_all_parser.py delete mode 100644 flare/modules/flare_parsers/gp_from_aimd_parser.py delete mode 100644 flare/modules/flare_parsers/old_vac_parser.py delete mode 100644 flare/modules/flare_parsers/otf_parser_v0.py delete mode 100644 flare/modules/flare_parsers/parse_output_temp.py delete mode 100644 flare/modules/gp_calculator.py delete mode 100644 flare/modules/gp_from_aimd.py delete mode 100644 flare/modules/gp_vs_deepmd.py delete mode 100644 flare/modules/hyp_opt.py delete mode 100644 flare/modules/initialize_velocities.py delete mode 100644 flare/modules/lammps/__init__.py delete mode 100644 flare/modules/lammps/coulomb_util.py delete mode 100644 flare/modules/lammps/eam.py delete mode 100644 flare/modules/lammps/mc.py delete mode 100644 flare/modules/qe_input.py delete mode 100644 flare/modules/qe_parsers.py delete mode 100644 flare/modules/rcut.py delete mode 100644 flare/modules/test_otf.py delete mode 100644 flare/modules/vacancy_diffusion.py delete mode 100644 flare/modules/validate.py rename flare/{modules/flare_parsers => }/otf_parser.py (99%) diff --git a/flare/modules/__init__.py b/flare/modules/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/flare/modules/__init__.pyc b/flare/modules/__init__.pyc deleted file mode 100644 index c1f2108f6c17ae3701ee962ac0fb1f4ff8f9abf7..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 121 zcmZSn%*&-(dN(GS0SXv_v;z dict: - """ - Get information about the run from the header of the file - :param outfile: - :return: - """ - with open(outfile, 'r') as f: - lines = f.readlines() - - header_info = {} - - stopreading = None - - for line in lines: - if '---' in line or '=' in line: - stopreading = lines.index(line) - break - - if stopreading is None: - raise Exception("OTF File is malformed") - - for line in lines[:stopreading]: - # TODO Update this in full - if 'frames' in line: - header_info['frames'] = int(line.split(':')[1]) - if 'kernel' in line: - header_info['kernel'] = line.split(':')[1].strip() - if 'number of hyperparameters:' in line: - header_info['n_hyps'] = int(line.split(':')[1]) - if 'optimization algorithm' in line: - header_info['algo'] = line.split(':')[1].strip() - if 'number of atoms' in line: - header_info['atoms'] = int(line.split(':')[1]) - if 'timestep' in line: - header_info['dt'] = float(line.split(':')[1]) - if 'system species' in line: - line = line.split(':')[1] - line = line.split("'") - - species = [item for item in line if item.isalpha()] - - header_info['species'] = set(species) - - return header_info - - -def parse_md_information(outfile: str = 'otf_run.out') -> (List[str], np.array, - np.array, np.array): - """ - Parse information about the OTF MD run - :param outfile: - :return: - """ - with open(outfile, 'r') as f: - lines = f.readlines() - - frame_indices = [lines.index(line) + 3 for line in lines if - line[:2] == '-F'] - - # Parse number of atoms and number of timesteps - md_nat_line = [line.split(':')[-1] for line in lines[0:20] if - 'number of atoms' in line] - md_nat = int(md_nat_line[0]) - - md_steps_line = [line.split(':')[-1] for line in lines[0:20] if - 'number of frames' in line] - md_steps = int(md_steps_line[0]) - - # Number of atoms is constant between frames and number of frames is known - run_positions = np.empty(shape=(md_steps, md_nat, 3)) - run_forces = np.empty(shape=(md_steps, md_nat, 3)) - run_stds = np.empty(shape=(md_steps, md_nat, 3)) - run_vels = np.empty(shape=(md_steps,md_nat,3)) - - # Assume species order is constant for MD run - # Get species from first frame - species = [] - first_frame = frame_indices[0] - for n in range(md_nat): - specie = strip_and_split(lines[first_frame + n])[0] - species.append(specie) - - # Get frame information and populate arrays - for frame_n, frame_index in enumerate(frame_indices): - - for at in range(md_nat): - data = strip_and_split(lines[frame_index + at]) - run_positions[frame_n, at, :] = data[1:4] - run_forces[frame_n, at, :] = data[4:7] - run_stds[frame_n, at, :] = data[7:10] - run_vels[frame_n,at,:] = data[10:13] - - return species, run_positions, run_forces, run_stds, run_vels - - -def parse_dft_information(outfile: str) -> (List[str], List[np.array], \ - List[np.array]): - """ - Parse the output of a otf run for analysis - :param outfile: str, Path to file - :return: dict{int:value,'species':list}, Dict of positions, forces, - vars indexed by frame and of species - """ - - with open(outfile, 'r') as f: - lines = f.readlines() - - # DFT indices start with ~ and are 2 lines above the data - dft_indices = [lines.index(line) + 3 for line in lines if line[:2] == '*-'] - # Frame indices start with - and are 2 lines above the data - - # Species and positions may vary from DFT run to DFT run - - # TODO turn into numpy array and extend array with each new data point; - # would be more efficient for large files - - dft_species = [] - dft_positions = [] - dft_forces = [] - dft_velocities = [] - - for frame_n, frame_index in enumerate(dft_indices): - - dft_run_nat = 0 - - dft_run_species = [] - dft_run_positions = [] - dft_run_forces = [] - dft_run_velocities = [] - - # Read atoms, positions, forces in - - while len(lines[frame_index + dft_run_nat])>=2: - dft_run_nat += 1 - - for at in range(dft_run_nat): - data = strip_and_split(lines[frame_index + at]) - - dft_run_species.append(data[0].strip(':')) - - curr_position = np.array([data[1], data[2], data[3]], - dtype='float') - curr_force = np.array([data[4], data[5], data[6]], dtype='float') - curr_vel = np.array([data[7], data[8], data[9]],dtype='float') - - dft_run_positions.append(curr_position) - dft_run_forces.append(curr_force) - dft_run_velocities.append(curr_vel) - - dft_species.append(dft_run_species) - dft_positions.append(np.array(dft_run_positions)) - dft_forces.append(np.array(dft_run_forces)) - dft_velocities.append(np.array(dft_run_velocities)) - - return dft_species, dft_positions, dft_forces, dft_velocities - - -def parse_hyperparameter_information(outfile: str = 'otf_run.out')-> List[ - np.array]: - - with open(outfile, 'r') as f: - lines = f.readlines() - - header_info = parse_header_information(outfile) - n_hyps = header_info['n_hyps'] - - hyp_indices = [i+1 for i,text in enumerate(lines) if 'GP hyper' in - text] - - hyperparameter_set = [np.zeros(n_hyps)]*len(hyp_indices) - - for set_index,line_index in enumerate(hyp_indices): - - for j in range(n_hyps): - - cur_hyp = lines[line_index+j].strip().split('=')[1] - - hyperparameter_set[set_index][j] = float(cur_hyp) - - return hyperparameter_set - - -def parse_time_information(outfile: str = 'otf_run.out'): - # with open(outfile, 'r') as f: - # lines = f.readlines() - - raise NotImplementedError - - -if __name__ == '__main__': - pass diff --git a/flare/modules/ase_calculator.py b/flare/modules/ase_calculator.py deleted file mode 100644 index c598850ce..000000000 --- a/flare/modules/ase_calculator.py +++ /dev/null @@ -1,109 +0,0 @@ -import numpy as np -import multiprocessing as mp -import concurrent.futures -import copy -import sys -sys.path.append('..') -from flare.env import AtomicEnvironment -from flare.struc import Structure -from flare.mff.mff_new import MappedForceField -from ase.calculators.calculator import Calculator - -class FLARE_Calculator(Calculator): - def __init__(self, gp_model, mff_model, use_mapping=False): - super().__init__() # all set to default values,TODO: change - self.mff_model = mff_model - self.gp_model = gp_model - self.use_mapping = use_mapping - self.results = {} - - def get_potential_energy(self, atoms=None, force_consistent=False): - # TODO: to be implemented - return 1 - - nat = len(atoms) - struc_curr = Structure(atoms.cell, ['A']*nat, - atoms.positions) - local_energies = np.zeros(nat) - - for n in range(nat): - chemenv = AtomicEnvironment(struc_curr, n, - self.gp_model.cutoffs) - local_energies[n] = self.gp_model.predict_local_energy(chemenv) - - return np.sum(local_energies) - - def get_forces(self, atoms): - if self.use_mapping: - return self.get_forces_mff(atoms) - else: - return self.get_forces_gp(atoms) - - def get_forces_gp(self, atoms): - nat = len(atoms) - struc_curr = Structure(atoms.cell, ['A']*nat, - atoms.positions) - - forces = np.zeros((nat, 3)) - stds = np.zeros((nat, 3)) - for n in range(nat): - chemenv = AtomicEnvironment(struc_curr, n, - self.gp_model.cutoffs) - for i in range(3): - force, std = self.gp_model.predict(chemenv, i + 1) - forces[n][i] = float(force) - stds[n][i] = np.sqrt(np.absolute(std)) - - self.results['stds'] = stds - atoms.get_uncertainties = self.get_uncertainties - - return forces - - def get_forces_mff(self, atoms): - nat = len(atoms) - struc_curr = Structure(atoms.cell, ['A']*nat, - atoms.positions) - - forces = np.zeros((nat, 3)) - stds = np.zeros((nat, 3)) - for n in range(nat): - chemenv = AtomicEnvironment(struc_curr, n, - self.mff_model.GP.cutoffs) - f, v = self.mff_model.predict(chemenv, mean_only=False) - forces[n] = f - stds[n] = np.sqrt(np.absolute(v)) - - self.results['stds'] = stds - atoms.get_uncertainties = self.get_uncertainties - return forces - - def get_stress(self, atoms): - return np.eye(3) - - def calculation_required(self, atoms, quantities): - return True - - def get_uncertainties(self): - return self.results['stds'] - - def train_gp(self, monitor=True): - self.gp_model.train(monitor) - - def build_mff(self, skip=True): - # l_bound not implemented - - if skip and (self.curr_step in self.non_mapping_steps): - return 1 - - # set svd rank based on the training set, grid number and threshold 1000 - grid_params = self.mff_model.grid_params - struc_params = self.mff_model.struc_params - - train_size = len(self.gp_model.training_data) - rank_2 = np.min([1000, grid_params['grid_num_2'], train_size*3]) - rank_3 = np.min([1000, grid_params['grid_num_3'][0]**3, train_size*3]) - grid_params['svd_rank_2'] = rank_2 - grid_params['svd_rank_3'] = rank_3 - - self.mff_model = MappedForceField(self.gp_model, grid_params, struc_params) - diff --git a/flare/modules/ase_otf.py b/flare/modules/ase_otf.py deleted file mode 100644 index 7086bcab8..000000000 --- a/flare/modules/ase_otf.py +++ /dev/null @@ -1,213 +0,0 @@ -import os -import sys -sys.path.append('..') -from copy import deepcopy - -from flare.struc import Structure - -import numpy as np -from ase.md.md import MolecularDynamics -from ase import units - - -class OTF(MolecularDynamics): - def __init__(self, atoms, timestep, - trajectory=None, - # on-the-fly parameters - dft_calc=None, dft_count=None, - std_tolerance_factor: float=1, - prev_pos_init: np.ndarray=None, par:bool=False, - skip: int=0, init_atoms: list=[], - calculate_energy=False, - max_atoms_added=1, freeze_hyps=1, - rescale_steps=[], rescale_temps=[], add_all=False, - no_cpus=1, - # mff parameters - use_mapping: bool=False, - non_mapping_steps: list=[], - l_bound: float=None, two_d: bool=False): - - ''' - output_name => logfile - npool: added - prev_pos_init: relaunch mode not implemented - par: consider merging into gp_calculator instead of here - skip: not implemented - calculate_energy: not implemented - l_bound: mff update l_bound, not implemented - - dft calculator is set outside of the otf module, and input as dft_calc, - so that different calculators can be used - - TODO: - 1. replace Structure with the ase internal Atoms - ''' - - MolecularDynamics.__init__(self, atoms, timestep, trajectory) - - self.std_tolerance = std_tolerance_factor - self.noa = len(atoms.positions) - self.max_atoms_added = max_atoms_added - self.freeze_hyps = freeze_hyps - self.dft_calc = dft_calc - if dft_count is None: - self.dft_count = 0 - else: - self.dft_count = dft_count - - # params for mapped force field - self.use_mapping = use_mapping - self.non_mapping_steps = non_mapping_steps - self.two_d = two_d - - # initialize local energies - if calculate_energy: - self.local_energies = np.zeros(self.noa) - else: - self.local_energies = None - - # initialize otf - if init_atoms is None: - self.init_atoms = [int(n) for n in range(self.noa)] - else: - self.init_atoms = init_atoms - - # rescale temperatures - self.rescale_steps = rescale_steps - self.rescale_temps = rescale_temps - - - def is_std_in_bound(self, atom_list): - # set uncertainty threshold - if self.std_tolerance == 0: - self.std_in_bound = True - self.target_atom = -1 - elif self.std_tolerance > 0: - noise = np.abs(self.atoms.calc.gp_model.hyps[-1]) - threshold = self.std_tolerance * noise - else: - threshold = np.abs(self.std_tolerance) - - # find max std - max_std = 0 - for atom, std in enumerate(self.stds): - std_curr = np.max(std) - - if std_curr > max_std: - max_std = std_curr - target_atom = atom - - # if above threshold, return atom - if max_std > threshold and atom not in atom_list: - self.std_in_bound = False - self.target_atom = target_atom - else: - self.std_in_bound = True - self.target_atom = -1 - - def otf_run(self, steps): - """Perform a number of time steps.""" - # initialize gp by a dft calculation - if not self.atoms.calc.gp_model.training_data: - self.dft_count = 0 - self.std_in_bound = False - self.target_atom = 0 - self.stds = [] - dft_forces = self.call_DFT() - - # update gp model - atom_struc = Structure(self.atoms.cell, - ['A']*len(self.atoms.positions), - self.atoms.positions) - self.atoms.calc.gp_model.update_db(atom_struc, dft_forces, - custom_range=self.init_atoms) - - # train calculator - self.train() - print(self.atoms.calc.mff_model) - - if not self.initialized: - self.initialize() - else: - if self.have_the_atoms_been_changed(): - raise NotImplementedError( - "You have modified the atoms since the last timestep.") - - for i in range(steps): - print('step:', i) - self.step() - self.nsteps += 1 - self.stds = self.atoms.get_uncertainties() - - # figure out if std above the threshold - self.call_observers() - self.is_std_in_bound([]) - - if not self.std_in_bound: - # call dft/eam - print('calling dft') - dft_forces = self.call_DFT() - - # update gp - print('updating gp') - self.update_GP(dft_forces) - - - def call_DFT(self): - prev_calc = self.atoms.calc - calc = deepcopy(self.dft_calc) - self.atoms.set_calculator(calc) - forces = self.atoms.get_forces() - self.call_observers() - self.atoms.set_calculator(prev_calc) - self.dft_count += 1 - return forces - - def update_GP(self, dft_forces): - atom_count = 0 - atom_list = [] - gp_model = self.atoms.calc.gp_model - while (not self.std_in_bound and atom_count < - self.max_atoms_added): - # build gp structure from atoms - atom_struc = Structure(self.atoms.cell, - ['A']*len(self.atoms.positions), - self.atoms.positions) - - # update gp model - gp_model.update_db(atom_struc, dft_forces, - custom_range=[self.target_atom]) - - if gp_model.alpha is None: - gp_model.set_L_alpha() - else: - gp_model.update_L_alpha() - - atom_list.append(self.target_atom) - forces = self.atoms.calc.get_forces_gp(self.atoms) # needed before get_uncertainties - self.stds = self.atoms.get_uncertainties() - - # write added atom to the log file, - # refer to ase.optimize.optimize.Dynamics - self.observers[0][0].add_atom_info(self.target_atom, - self.stds[self.target_atom]) - - self.is_std_in_bound(atom_list) - atom_count += 1 - - self.train() - - def train(self, monitor=True, skip=False): - calc = self.atoms.calc - if (self.dft_count-1) < self.freeze_hyps: - calc.gp_model.train(monitor) - self.observers[0][0].write_hyps(calc.gp_model.hyp_labels, - calc.gp_model.hyps, calc.gp_model.like, - calc.gp_model.like_grad) - else: - calc.gp_model.set_L_alpha() - - # build mapped force field - if self.use_mapping: - calc.build_mff(skip) - diff --git a/flare/modules/ase_otf_logger.py b/flare/modules/ase_otf_logger.py deleted file mode 100644 index 994b25202..000000000 --- a/flare/modules/ase_otf_logger.py +++ /dev/null @@ -1,105 +0,0 @@ -import numpy as np -import datetime -import time - -from ase.md import MDLogger -from ase import units - -from flare.modules.gp_calculator import GPCalculator - -class OTFLogger(MDLogger): - - def __init__(self, dyn, atoms, logfile, header=False, stress=False, - peratom=False, mode="w"): - - super().__init__(dyn, atoms, logfile, header, stress, - peratom, mode) - - self.write_header_info() - self.start_time = time.time() - - - def write_header_info(self): - gp_model = self.atoms.calc.gp_model - self.logfile.write(str(datetime.datetime.now())) - self.logfile.write('\nnumber of cpu cores: ') # TODO - self.logfile.write('\ncutoffs: '+str(gp_model.cutoffs)) - self.logfile.write('\nkernel: '+gp_model.kernel.__name__) - self.logfile.write('\nnumber of parameters: '+str(len(gp_model.hyps))) - self.logfile.write('\nhyperparameters: '+str(gp_model.hyps)) - self.logfile.write('\nhyperparameter optimization algorithm: '+gp_model.algo) - self.logfile.write('\nuncertainty tolerance: '+str(self.dyn.std_tolerance)) - self.logfile.write('\nperiodic cell:\n'+str(self.atoms.cell)) - self.logfile.write('\n') - - def write_hyps(self, hyp_labels, hyps, like, like_grad): - self.logfile.write('\nGP hyperparameters: \n') - for i, label in enumerate(jyp_labels): - self.logfile.write('Hyp{} : {} = {}\n'.format(i, label, hyps[i])) - - self.logfile.write('likelihood: '+str(like)+'\n') - self.logfile.write('likelihood gradient: '+str(like_grad)+'\n') - self.logfile.write('wall time from start: %.2f s \n' % time.time()-self.start_time) - - def __call__(self): - self.logfile.write('\n'+50*'-') - if self.dyn is not None: - steps = self.dyn.get_time() - t = steps / (1000*units.fs) - if type(self.atoms.calc) != GPCalculator: # TODO - self.logfile.write('\n-*Frame: '+str(steps)) - else: - self.logfile.write('\n-Frame: '+str(steps)) - self.logfile.write('\nSimulation time: '+str(t)) - self.logfile.write('\n') - - # add positions, forces and stds to be written - positions = self.atoms.get_positions() - forces = self.atoms.get_forces() - velocities = self.atoms.get_velocities() - if type(self.atoms.calc) == GPCalculator: # TODO: make it compatible with mff - stds = self.atoms.get_uncertainties() - force_str = 'GP Forces' - else: - stds = np.zeros(positions.shape) - force_str = 'DFT Forces' - self.logfile.write((10*' ')+'Positions'+(11*' '+'|'+11*' ')+ - force_str+(11*' '+'|'+11*' ')+ - 'Velocities'+(10*' '+'|'+11*' ')+ - 'Uncertainties\n') - - template = '{:9f} {:9f} {:9f} | {:9f} {:9f} {:9f} | {:9f} {:9f} {:9f} | {:9f} {:9f} {:9f}' - for atom in range(len(positions)): - dat = template.format(\ - positions[atom][0], positions[atom][1], positions[atom][2], - forces[atom][0], forces[atom][1], forces[atom][2], - velocities[atom][0],velocities[atom][1],velocities[atom][2], - stds[atom][0], stds[atom][1], stds[atom][2]) - self.logfile.write(dat+'\n') - - # write energy, temperature info - epot = self.atoms.get_potential_energy() - ekin = self.atoms.get_kinetic_energy() - temp = ekin / (1.5 * units.kB * self.natoms) - if self.peratom: - epot /= self.natoms - ekin /= self.natoms - self.logfile.write('\ntotal energy: '+str(epot+ekin)) - self.logfile.write('\nkinetic energy: '+str(ekin)) - self.logfile.write('\ntemperature: '+str(temp)) - self.logfile.write('\nwall time from start: '+str(time.time()-self.start_time)) - - self.logfile.write('\n') - self.logfile.flush() - - - def add_atom_info(self, target_atom, uncertainty): - self.logfile.write('\nAdding atom '+str(target_atom)+' to the training set with uncertainty: '+str(uncertainty)) - - def write_mff_train(self, mff_model, train_time): - train_size = len(mff_model.GP.training_data) - self.logfile.write('\ntraining set size: {}\n'.format(train_size)) - self.logfile.write('lower bound: {}\n'.format(mff_model.l_bound)) - self.logfile.write('mff l_bound: {}\n'.format(mff_model.grid_params['bounds_2'][0,0])) - self.logfile.write('building mapping time: {}'.format(train_time)) - diff --git a/flare/modules/ase_otf_md.py b/flare/modules/ase_otf_md.py deleted file mode 100644 index 5cc7d0778..000000000 --- a/flare/modules/ase_otf_md.py +++ /dev/null @@ -1,61 +0,0 @@ -import os -import sys -sys.path.append('../..') -from flare.struc import Structure -from flare.modules.ase_otf import OTF - -import numpy as np -from ase.calculators.espresso import Espresso -from ase.calculators.eam import EAM -from ase.md.npt import NPT -from ase.md.nvtberendsen import NVTBerendsen -from ase.md.nptberendsen import NPTBerendsen -from ase.md.verlet import VelocityVerlet -from ase.md.md import MolecularDynamics -from ase import units - - -class OTF_VelocityVerlet(VelocityVerlet, OTF): - def __init__(self, atoms, timestep=None, trajectory=None, dt=None, - **kwargs): - - VelocityVerlet.__init__(self, atoms, timestep, trajectory, - dt=dt) - - OTF.__init__(self, atoms, timestep, trajectory, **kwargs) - - -class OTF_NVTBerendsen(NVTBerendsen, OTF): - def __init__(self, atoms, timestep, temperature, taut, fixcm=True, - trajectory=None, **kwargs): - - NVTBerendsen.__init__(self, atoms, timestep, temperature, taut, - fixcm, trajectory) - - OTF.__init__(self, atoms, timestep, trajectory, **kwargs) - - -class OTF_NPTBerendsen(NPTBerendsen, OTF): - def __init__(self, atoms, timestep, temperature, taut=0.5e3 * - units.fs, pressure=1.01325, taup=1e3 * units.fs, - compressibility=4.57e-5, fixcm=True, trajectory=None, - **kwargs): - - NPTBerendsen.__init__(self, atoms, timestep, temperature, taut, - pressure, taup, - compressibility, fixcm, trajectory) - - OTF.__init__(self, atoms, timestep, trajectory, **kwargs) - - -class OTF_NPT(NPT, OTF): - def __init__(self, atoms, timestep, temperature, externalstress, - ttime, pfactor, mask=None, trajectory=None, **kwargs): - - NPT.__init__(self, atoms, timestep, temperature, - externalstress, ttime, pfactor, mask, - trajectory) - - OTF.__init__(self, atoms, timestep, trajectory, **kwargs) - - diff --git a/flare/modules/atomeye_helper.py b/flare/modules/atomeye_helper.py deleted file mode 100644 index e2ee7bcec..000000000 --- a/flare/modules/atomeye_helper.py +++ /dev/null @@ -1,297 +0,0 @@ -import numpy as np -import math -from copy import copy -from typing import List -import subprocess -from flare import struc - - -def write_cfgs_from_otf(otf_run, cell, species, start_frame, no_frames, - folder_name, image_quality, scr_dest, trans_vec, - skip=1): - - # make folder for cfgs - subprocess.call('mkdir %s' % folder_name, shell=True) - subprocess.call('mkdir %s/Pic' % folder_name, shell=True) - - scr_anim_text = '%s\n' % image_quality - - # make cfgs - no_digits = int(np.ceil(math.log(no_frames, 10))) - - for n in range(no_frames): - frame_no = n - frame_no_padded = str(frame_no).zfill(no_digits) - frame_string = frame_no_padded+'.cfg' - frame_dest = folder_name+'/'+frame_string - - scr_anim_text += '%s %s/Pic/%s.jpg\n' % \ - (frame_dest, folder_name, frame_no_padded) - - positions = otf_run.position_list[start_frame + n * skip] - - write_cfg_file(frame_dest, positions, species, cell) - - # write animation directions for AtomEye - write_file(scr_dest, scr_anim_text) - - -def write_cfgs_from_md(md_trajectory, start_frame, no_frames, folder_name, - image_quality, scr_dest, trans_vec, skip=1): - # make folder for cfgs - subprocess.call('mkdir %s' % folder_name, shell=True) - subprocess.call('mkdir %s/Pic' % folder_name, shell=True) - - scr_anim_text = '%s\n' % image_quality - - # make cfgs - no_digits = int(np.ceil(math.log(no_frames, 10))) - cell = md_trajectory.cell - for n in range(no_frames): - frame_no = n - frame_no_padded = str(frame_no).zfill(no_digits) - frame_string = frame_no_padded+'.cfg' - frame_dest = folder_name+'/'+frame_string - - scr_anim_text += '%s %s/Pic/%s.jpg\n' % \ - (frame_dest, folder_name, frame_no_padded) - - positions = \ - np.array(md_trajectory.MD_data[start_frame+n*skip] - ['positions']) + \ - trans_vec - species = md_trajectory.MD_data[start_frame+n*skip]['elements'] - - write_cfg_file(frame_dest, positions, species, cell) - - # write animation directions for AtomEye - write_file(scr_dest, scr_anim_text) - - -def write_cfgs_from_pos(pos_list, cell, folder_name, image_quality, scr_dest, - species): - # make folder for cfgs - subprocess.call('mkdir %s' % folder_name, shell=True) - subprocess.call('mkdir %s/Pic' % folder_name, shell=True) - - scr_anim_text = '%s\n' % image_quality - - # make cfgs - no_frames = len(pos_list) - no_digits = int(np.ceil(math.log(no_frames, 10))) - - for n, pos in enumerate(pos_list): - # pos = np.load(pos_file) - frame_no = n - frame_no_padded = str(frame_no).zfill(no_digits) - frame_string = frame_no_padded+'.cfg' - frame_dest = folder_name+'/'+frame_string - - scr_anim_text += '%s %s/Pic/%s.jpg\n' % \ - (frame_dest, folder_name, frame_no_padded) - - write_cfg_file(frame_dest, pos, species, cell) - - # write animation directions for AtomEye - write_file(scr_dest, scr_anim_text) - - -def write_cfg_file(file_name: str, positions: np.ndarray, species: List[str], - cell: np.ndarray) -> None: - """write cfg file that can be interpreted by AtomEye. - assumes orthorombic unit cell. - - :param file_name: destination of cfg file - :type file_name: str - :param positions: positions of atoms - :type positions: np.ndarray - :param species: atom species - :type species: List[str] - :param cell: unit cell vectors - :type cell: np.ndarray - :return: creates the cfg file - :rtype: None - """ - - cfg_text = get_cfg_text(positions, species, cell) - write_file(file_name, cfg_text) - - -def get_cfg_text(positions: np.ndarray, species: List[str], - cell: np.ndarray) -> str: - """returns cfg text - - :param positions: Cartesian coordinates of atomic positions - :type positions: np.ndarray - :param species: list of atomic species (determines atom size) - :param cell: cell of unit vectors - :type cell: np.ndarray - :return: cfg text - :rtype: str - """ - - cfg_header = get_cfg_header(positions.shape[0], cell) - reduced_coordinates = calculate_reduced_coordinates(positions, cell) - position_text = get_reduced_coordinate_text(reduced_coordinates, species) - cfg_text = cfg_header + position_text - - return cfg_text - - -def get_cfg_header(number_of_particles: int, cell: np.ndarray) -> str: - """creates cfg header from atom positions and unit cell. - assumes unit cell is orthorombic. - - :param positions: Nx3 array of atomic positions - :type positions: np.ndarray - :param cell: 3x3 array of cell vectors (cell vectors are rows) - :type cell: np.ndarray - """ - - cfg_text = """Number of particles = %i -# (required) this must be the first line - -A = 1.0 Angstrom (basic length-scale) -# (optional) basic length-scale: default A = 1.0 [Angstrom] - -H0(1,1) = %f A -H0(1,2) = %f A -H0(1,3) = %f A -# (required) this is the supercell's 1st edge, in A - -H0(2,1) = %f A -H0(2,2) = %f A -H0(2,3) = %f A -# (required) this is the supercell's 2nd edge, in A - -H0(3,1) = %f A -H0(3,2) = %f A -H0(3,3) = %f A -# (required) this is the supercell's 3rd edge, in A - -Transform(1,1) = 1 -Transform(1,2) = 0 -Transform(1,3) = 0 -Transform(2,1) = 0 -Transform(2,2) = 1 -Transform(2,3) = 0 -Transform(3,1) = 0 -Transform(3,2) = 0 -Transform(3,3) = 1 -# (optional) apply additional transformation on H0: H = H0 * Transform; -# default = Identity matrix. - -eta(1,1) = 0 -eta(1,2) = 0 -eta(1,3) = 0 -eta(2,2) = 0 -eta(2,3) = 0 -eta(3,3) = 0 -# (optional) apply additional Lagrangian strain on H0: -# H = H0 * sqrt(Identity_matrix + 2 * eta); -# default = zero matrix. - -# ENSUING ARE THE ATOMS, EACH ATOM DESCRIBED BY A ROW -# 1st entry is atomic mass in a.m.u. -# 2nd entry is the chemical symbol (max 2 chars) - -# 3rd entry is reduced coordinate s1 (dimensionless) -# 4th entry is reduced coordinate s2 (dimensionless) -# 5th entry is reduced coordinate s3 (dimensionless) -# real coordinates x = s * H, x, s are 1x3 row vectors - -# 6th entry is d(s1)/dt in basic rate-scale R -# 7th entry is d(s2)/dt in basic rate-scale R -# 8th entry is d(s3)/dt in basic rate-scale R -R = 1.0 [ns^-1] -# (optional) basic rate-scale: default R = 1.0 [ns^-1] -""" % (number_of_particles, - cell[0, 0], cell[0, 1], cell[0, 2], - cell[1, 0], cell[1, 1], cell[1, 2], - cell[2, 0], cell[2, 1], cell[2, 2]) - - return cfg_text - - -def calculate_reduced_coordinates(positions: np.ndarray, - cell: np.ndarray) -> np.ndarray: - """convert raw cartesian coordinates to reduced coordinates with each atom - wrapped back into the unit cell. assumes unit cell is orthorombic. - - :param positions: Nx3 array of atomic positions - :type positions: np.ndarray - :param cell: 3x3 array of cell vectors (cell vectors are rows) - :type cell: np.ndarray - :return: Nx3 array of reduced coordinates - :rtype: np.ndarray - """ - - species = ['Co'] * len(positions) - - structure = struc.Structure(cell, species, positions) - rel_pos = \ - structure.raw_to_relative(structure.positions, - structure.cell_transpose, - structure.cell_dot_inverse) - - rel_wrap = rel_pos - np.floor(rel_pos) - - # reduced_coordinates = np.zeros((positions.shape[0], 3)) - # for m in range(positions.shape[0]): - # for n in range(3): - # trial_coord = positions[m, n] / cell[n, n] - - # # reduced coordinates must be between 0 and 1 - # trans = np.floor(trial_coord) - # reduced_coordinates[m, n] = trial_coord - trans - - return rel_wrap - - -def get_reduced_coordinate_text(reduced_coordinates: np.ndarray, - species: List[str]) -> str: - """records reduced coordinates in cfg format. - - :param reduced_coordinates: array of reduced coordinates - :type reduced_coordinates: np.ndarray - :param species: list of atomic species, which determines atom size - :type species: List[str] - :return: cfg string of reduced coordinates. - :rtype: str - """ - - reduced_text = """ -# ENSUING ARE THE ATOMS, EACH ATOM DESCRIBED BY A ROW -# 1st entry is atomic mass in a.m.u. -# 2nd entry is the chemical symbol (max 2 chars) - -# 3rd entry is reduced coordinate s1 (dimensionless) -# 4th entry is reduced coordinate s2 (dimensionless) -# 5th entry is reduced coordinate s3 (dimensionless) -# real coordinates x = s * H, x, s are 1x3 row vectors - -# 6th entry is d(s1)/dt in basic rate-scale R -# 7th entry is d(s2)/dt in basic rate-scale R -# 8th entry is d(s3)/dt in basic rate-scale R -R = 1.0 [ns^-1] -# (optional) basic rate-scale: default R = 1.0 [ns^-1] -""" - - for spec, coord in zip(species, reduced_coordinates): - # use arbitrary mass, label - reduced_text += \ - '1.0 %s %f %f %f 0 0 0 \n' % (spec, coord[0], coord[1], coord[2]) - - return reduced_text - - -def write_file(file_name: str, text: str): - with open(file_name, 'w') as fin: - fin.write(text) - - -if __name__ == '__main__': - reduced_coordinates = np.array([[1, 2, 3], [4, 5, 6]]) - species = ['A', 'B'] - test = get_reduced_coordinate_text(reduced_coordinates, species) - print(test) diff --git a/flare/modules/calculate_rdf.py b/flare/modules/calculate_rdf.py deleted file mode 100644 index 05507036f..000000000 --- a/flare/modules/calculate_rdf.py +++ /dev/null @@ -1,89 +0,0 @@ -import numpy as np -import sys -from flare import struc, env - - -def calculate_rdf(position_list, cell, species, snaps, cutoff, bins, - cell_vol=None): - - # assume cubic cell by default - if cell_vol is None: - cell_vol = cell[0, 0] ** 3 - - # collect interatomic distances - r_list = [] - delta_r = cutoff / bins - atom_count = 0 - nat = position_list[0].shape[0] - cutoffs = np.array([cutoff]) - for snap in snaps: - positions = position_list[snap] - structure = struc.Structure(cell, species, positions) - - for n in range(len(positions)): - env_curr = env.AtomicEnvironment(structure, n, cutoffs) - atom_count += 1 - for bond in env_curr.bond_array_2: - r_list.append(bond[0]) - r_list = np.array(r_list) - radial_hist, _ = \ - np.histogram(r_list, bins=bins, range=(0, cutoff)) - - # weight the histogram - rs = np.linspace(delta_r/2, cutoff-delta_r/2, bins) - rho = nat / cell_vol - weights = (4 * np.pi * rho / 3) * ((rs+delta_r)**3 - rs**3) - rad_dist = radial_hist / (atom_count * weights) - - return rs, rad_dist, atom_count - - -def calculate_species_rdf(position_list, spec1, spec2, cell, species, snaps, - cutoff, bins, cell_vol=None): - - # assume cubic cell by default - if cell_vol is None: - cell_vol = cell[0, 0] ** 3 - - # collect interatomic distances - r_list = [] - delta_r = cutoff / bins - atom_count = 0 - nat = position_list[0].shape[0] - cutoffs = np.array([cutoff]) - - # compute concentration of species 2 - positions = position_list[snaps[0]] - struc_ex = struc.Structure(cell, species, positions) - spec2_count = 0 - for spec in struc_ex.coded_species: - if spec == spec2: - spec2_count += 1 - spec2_conc = spec2_count / nat - - for snap in snaps: - positions = position_list[snap] - structure = struc.Structure(cell, species, positions) - - for n in range(len(positions)): - env_curr = env.AtomicEnvironment(structure, n, cutoffs) - ctype = env_curr.ctype - - if ctype == spec1: - atom_count += 1 - - for bond, spec in zip(env_curr.bond_array_2, env_curr.etypes): - if spec == spec2: - r_list.append(bond[0]) - - r_list = np.array(r_list) - radial_hist, _ = \ - np.histogram(r_list, bins=bins, range=(0, cutoff)) - - # weight the histogram - rs = np.linspace(delta_r/2, cutoff-delta_r/2, bins) - rho = (nat * spec2_conc) / cell_vol - weights = (4 * np.pi * rho / 3) * ((rs+delta_r)**3 - rs**3) - rad_dist = radial_hist / (atom_count * weights) - - return rs, rad_dist, atom_count diff --git a/flare/modules/convergence.py b/flare/modules/convergence.py deleted file mode 100644 index 74774e7c1..000000000 --- a/flare/modules/convergence.py +++ /dev/null @@ -1,177 +0,0 @@ -import numpy as np -import sys -import os -import qe_input -import qe_parsers -import datetime -import time - - -def convergence(txt_name, input_file_string, output_file_string, pw_loc, - calculation, scf_inputs, nks, ecutwfcs, rho_facs, - electron_maxstep=100, metal=False): - # initialize update text - txt = update_init() - write_file(txt_name, txt) - - # converged run - initial_input_name = input_file_string + '.in' - initial_output_name = output_file_string + '.out' - - scf = qe_input.QEInput(initial_input_name, initial_output_name, pw_loc, - calculation, scf_inputs, - electron_maxstep=electron_maxstep, - metal=metal) - - # perform and time converged scf run - time0 = time.time() - scf.run_espresso(npool=False) - time1 = time.time() - scf_time = time1 - time0 - - conv_energy, conv_forces = qe_parsers.parse_scf(initial_output_name) - - txt += record_results(scf_inputs['kvec'][0], scf_inputs['ecutwfc'], - scf_inputs['ecutrho'], conv_energy, conv_forces, - conv_energy, conv_forces, scf_time) - write_file(txt_name, txt) - - # loop over parameters - for nk in nks: - for ecutwfc in ecutwfcs: - for rho_fac in rho_facs: - ecutrho = ecutwfc * rho_fac - scf_inputs['kvec'] = np.array([nk, nk, nk]) - scf_inputs['ecutrho'] = ecutrho - scf_inputs['ecutwfc'] = ecutwfc - - input_name = input_file_string + \ - '_nk%i_e%i_rho%i.in' % (nk, ecutwfc, rho_fac) - output_name = output_file_string + \ - '_nk%i_e%i_rho%i.out' % (nk, ecutwfc, rho_fac) - - scf = qe_input.QEInput(input_name, output_name, - pw_loc, calculation, scf_inputs, - electron_maxstep=electron_maxstep, - metal=metal) - - # perform and time scf run - time0 = time.time() - scf.run_espresso(npool=False) - time1 = time.time() - scf_time = time1 - time0 - - energy, forces = qe_parsers.parse_scf(output_name) - - txt += record_results(nk, ecutwfc, ecutrho, energy, forces, - conv_energy, conv_forces, scf_time) - write_file(txt_name, txt) - - txt += update_fin() - write_file(txt_name, txt) - - # remove output directory - if os.path.isdir('output'): - os.system('rm -r output') - if os.path.isdir('__pycache__'): - os.system('rm -r __pycache__') - - -def reshape_forces(forces: list) -> np.ndarray: - forces_array = np.array(forces) - forces_array = forces_array.reshape(forces_array.size) - return forces_array - - -def print_results(nk, ecutwfc, ecutrho, energy, forces, conv_energy, - conv_forces): - print('\n') - print('nk: %i' % nk) - print('ecutwfc: %.2f' % ecutwfc) - print('ecutrho: %.2f' % ecutrho) - print('energy: %f' % energy) - - en_diff = energy - conv_energy - print('energy difference from converged value: %.2e eV' % en_diff) - - # reshape converged force - conv_forces_array = reshape_forces(conv_forces) - forces_array = reshape_forces(forces) - force_diff = conv_forces_array - forces_array - force_MAE = np.mean(np.abs(force_diff)) - print('force MAE: %.2e eV/A' % force_MAE) - - max_err = np.max(np.abs(force_diff)) - print('max force error: %.2e eV/A' % max_err) - print('\n') - - -def record_results(nk, ecutwfc, ecutrho, energy, forces, conv_energy, - conv_forces, scf_time): - txt = """ - -Inputs: -nk: %i -ecutwfc: %.2f -ecutrho: %.2f - -Convergence results: -energy: %f eV.""" % (nk, ecutwfc, ecutrho, energy) - - en_diff = energy - conv_energy - txt += """ -energy difference from converged value: %.2e eV.""" % en_diff - - # reshape converged force - conv_forces_array = reshape_forces(conv_forces) - forces_array = reshape_forces(forces) - force_diff = conv_forces_array - forces_array - force_MAE = np.mean(np.abs(force_diff)) - txt += """ -force MAE: %.2e eV/A.""" % force_MAE - - max_err = np.max(np.abs(force_diff)) - txt += """ -max force error: %.2e eV/A. -scf time: %.2f s. -""" % (max_err, scf_time) - - return txt - - -# ----------------------------------------------------------------------------- -# monitor code progress -# ----------------------------------------------------------------------------- - -def write_file(fname, text): - with open(fname, 'w') as fin: - fin.write(text) - - -def update_init(): - init_text = """Quantum Espresso convergence test: %s. -Author: Jonathan Vandermause. -""" % str(datetime.datetime.now()) - - return init_text - - -def update_fin(): - fin_text = """ -------------------------------------------------------------------------------- - JOB DONE. -------------------------------------------------------------------------------- -""" - return fin_text - -if __name__ == '__main__': - nk = 5 - ecutwfc = 50 - ecutrho = 100 - energy = 24.5 - forces = [np.array([1, 2, 3])] - conv_energy = 25 - conv_forces = [np.array([4, 5, 6])] - record_test = record_results(nk, ecutwfc, ecutrho, energy, forces, - conv_energy, conv_forces, 4) - print(record_test) diff --git a/flare/modules/crystals.py b/flare/modules/crystals.py deleted file mode 100644 index c96e81bad..000000000 --- a/flare/modules/crystals.py +++ /dev/null @@ -1,125 +0,0 @@ -import numpy as np -from ase.build import fcc111, add_adsorbate -from ase.visualize import view -from ase.io import write - - -def get_supercell_positions(sc_size, cell, positions): - sc_positions = [] - for m in range(sc_size): - vec1 = m * cell[0] - for n in range(sc_size): - vec2 = n * cell[1] - for p in range(sc_size): - vec3 = p * cell[2] - - # append translated positions - for pos in positions: - sc_positions.append(pos+vec1+vec2+vec3) - - return sc_positions - - -def supercell_custom(cell, positions, size1, size2, size3): - sc_positions = [] - for m in range(size1): - vec1 = m * cell[0] - for n in range(size2): - vec2 = n * cell[1] - for p in range(size3): - vec3 = p * cell[2] - - # append translated positions - for pos in positions: - sc_positions.append(pos+vec1+vec2+vec3) - - return np.array(sc_positions) - - -# ----------------------------------------------------------------------------- -# fcc helper functions -# ----------------------------------------------------------------------------- - -def fcc_positions(cube_lat): - positions = [np.array([0, 0, 0]), - np.array([cube_lat/2, cube_lat/2, 0]), - np.array([0, cube_lat/2, cube_lat/2]), - np.array([cube_lat/2, 0, cube_lat/2])] - return positions - -# ----------------------------------------------------------------------------- -# diamond helper functions -# ----------------------------------------------------------------------------- - - -def cubic_diamond_positions(cube_lat): - positions = [np.array([0, 0, 0]), - np.array([cube_lat/2, cube_lat/2, 0]), - np.array([0, cube_lat/2, cube_lat/2]), - np.array([cube_lat/2, 0, cube_lat/2]), - np.array([cube_lat/4, cube_lat/4, cube_lat/4]), - np.array([3*cube_lat/4, 3*cube_lat/4, cube_lat/4]), - np.array([cube_lat/4, 3*cube_lat/4, 3*cube_lat/4]), - np.array([3*cube_lat/4, cube_lat/4, 3*cube_lat/4])] - return positions - - -def cubic_rocksalt_positions(cube_lat): - positions = np.array([[0, 0, 0], - [cube_lat/2, 0, 0], - [cube_lat/2, cube_lat/2, 0], - [cube_lat, cube_lat/2, 0], - [cube_lat/2, 0, cube_lat/2], - [cube_lat, 0, cube_lat/2], - [0, cube_lat/2, cube_lat/2], - [cube_lat/2, cube_lat/2, cube_lat/2]]) - - return positions - - -def primitive_diamond_positions(prim_lat): - positions = [np.array([0, 0, 0]), - np.array([prim_lat/2, prim_lat/2, prim_lat/2])] - return positions - - -# ----------------------------------------------------------------------------- -# slab helper functions -# ----------------------------------------------------------------------------- - -def get_fcc111_slab(layers, size, element, vacuum): - slab = fcc111(element, size=(size, size, layers), vacuum=vacuum) - return slab - - -def fcc111_and_adsorbate(layers, size, element, vacuum, height, position): - slab = fcc111(element, size=(size, size, layers)) - add_adsorbate(slab, element, height, position) - slab.center(vacuum=vacuum, axis=2) - return slab - - -# ----------------------------------------------------------------------------- -# water helper functions -# ----------------------------------------------------------------------------- - - -def water_coordinates(ox_pos: np.ndarray, - theta: float, phi: float) -> list: - H_angle = 104.45 * (2*np.pi / 360) - OH_len = 95.84e-12 - pass - - -if __name__ == '__main__': - layers = 2 - size = 2 - element = 'Pd' - vacuum = 10 - height = 1 - position = 'hcp' - - slab_test = fcc111_and_adsorbate(layers, size, element, vacuum, height, - position) - print(slab_test.positions) - print(slab_test.cell) diff --git a/flare/modules/cutoff_sweep.py b/flare/modules/cutoff_sweep.py deleted file mode 100644 index edace0a2a..000000000 --- a/flare/modules/cutoff_sweep.py +++ /dev/null @@ -1,143 +0,0 @@ -import numpy as np -import sys -from flare import gp, env, struc, kernels -from flare.modules import analyze_gp, qe_parsers -from mc_kernels import mc_simple -from scipy.optimize import minimize -import time -import datetime - - -def sweep(txt_name, data_file, cell, training_snaps, cutoffs, kernel, - kernel_grad, initial_hyps, par): - - # set up text file - txt = update_init() - write_file(txt_name, txt) - - # define md_trajectory object - md_trajectory = analyze_gp.MDAnalysis(data_file, cell) - - # set up training run - hyps = initial_hyps - - for cutoff in cutoffs: - gp_test = \ - analyze_gp.get_gp_from_snaps(md_trajectory, training_snaps, - kernel, kernel_grad, hyps, cutoff, - par=par) - gp_test.algo = 'BFGS' - gp_test.hyps = hyps - - # train gp model - time0 = time.time() - gp_test.train(monitor=True) - time1 = time.time() - training_time = time1 - time0 - - likelihood = gp_test.like - - hyps = gp_test.hyps - - txt += """\n - cutoff: {} - optimized hyperparameters: - """.format(cutoff) - txt += str(hyps) - txt += """ - likelihood: %.5f - training time: %.2f s""" % (likelihood, training_time) - - write_file(txt_name, txt) - - -def sweep_and_test(txt_name, data_file, cell, training_snaps, cutoffs, kernel, - kernel_grad, initial_hyps, par, test_snaps): - # set up text file - txt = update_init() - write_file(txt_name, txt) - - # define md_trajectory object - md_trajectory = analyze_gp.MDAnalysis(data_file, cell) - - # set up training run - hyps = initial_hyps - - for cutoff in cutoffs: - gp_test = \ - analyze_gp.get_gp_from_snaps(md_trajectory, training_snaps, - kernel, kernel_grad, hyps, cutoff, - par=par) - - gp_test.algo = 'BFGS' - gp_test.hyps = hyps - - # train gp model - time0 = time.time() - gp_test.train(monitor=True) - time1 = time.time() - training_time = time1 - time0 - - likelihood = gp_test.like - - hyps = gp_test.hyps - - txt += """\n - cutoff: {} - optimized hyperparameters: - """.format(cutoff) - txt += str(hyps) - txt += """ - likelihood: %.5f - training time: %.2f s""" % (likelihood, training_time) - - write_file(txt_name, txt) - - # test model - all_predictions, all_variances, all_forces = \ - analyze_gp.predict_forces_on_test_set(gp_test, md_trajectory, - test_snaps, cutoff) - - training_set_size = len(training_snaps) * 32 * 3 - avg_force = np.mean(np.abs(all_forces)) - max_force = np.max(np.abs(all_forces)) - mae = np.mean(np.abs(all_predictions - all_forces)) - max_err = np.max(np.abs(all_predictions - all_forces)) - avg_std = np.mean(np.sqrt(all_variances)) - max_std = np.max(np.sqrt(all_variances)) - - txt += """\n -training_set_size = %i -average force = %.4f -max force = %.4f -mean absolute error = %.4f -max error = %.4f -average std = %.4f -max std = %.4f -\n""" % (training_set_size, avg_force, max_force, mae, max_err, avg_std, - max_std) - - write_file(txt_name, txt) - - -def write_file(fname, text): - with open(fname, 'w') as fin: - fin.write(text) - - -def update_init(): - init_text = """Cutoff test. -Date and time: %s. -Author: Jonathan Vandermause. -""" % str(datetime.datetime.now()) - - return init_text - - -def update_fin(): - fin_text = """ -------------------------------------------------------------------------------- - JOB DONE. -------------------------------------------------------------------------------- -""" - return fin_text diff --git a/flare/modules/data_sweep.py b/flare/modules/data_sweep.py deleted file mode 100644 index eb9704145..000000000 --- a/flare/modules/data_sweep.py +++ /dev/null @@ -1,155 +0,0 @@ -import numpy as np -import sys -from flare import gp, struc -from flare.modules import analyze_gp -import time -import datetime - - -def data_sweep(data_file, cell, training_snaps, cutoffs, kernel, kernel_grad, - hyps, test_snaps, txt_name): - """Uses coordinates and forces from a Quantum Espresso AIMD output file to - test the performance of a GP model as a function of the training set size. - """ - - md_trajectory = analyze_gp.MDAnalysis(data_file, cell) - - # set up text file - txt = update_init() - write_file(txt_name, txt) - - # initialize gp model - gp_model = gp.GaussianProcess(kernel, kernel_grad, hyps, - cutoffs, opt_algorithm='BFGS') - - for n, snap in enumerate(training_snaps): - # update gp - time0 = time.time() - structure = md_trajectory.get_structure_from_snap(snap) - forces = md_trajectory.get_forces_from_snap(snap) - gp_model.update_db(structure, forces) - gp_model.set_L_alpha() - time1 = time.time() - - # predict on test set - all_predictions, all_variances, all_forces = \ - analyze_gp.predict_forces_on_test_set(gp_model, md_trajectory, - test_snaps, cutoffs) - time2 = time.time() - - update_time = time1 - time0 - prediction_time = time2 - time1 - - # record predictions for parity plots - if n == 0: - np.save('aimd_forces', all_forces) - pred_name = 'preds_' + str(snap) - vars_name = 'vars_' + str(snap) - np.save(pred_name, all_predictions) - np.save(vars_name, all_variances) - - # compute and report error - training_set_size = len(gp_model.training_labels_np) - avg_force = np.mean(np.abs(all_forces)) - max_force = np.max(np.abs(all_forces)) - mae = np.mean(np.abs(all_predictions - all_forces)) - max_err = np.max(np.abs(all_predictions - all_forces)) - avg_std = np.mean(np.sqrt(all_variances)) - max_std = np.max(np.sqrt(all_variances)) - - txt += """\n -training_set_size = %i -average force = %.4f -max force = %.4f -mean absolute error = %.4f -max error = %.4f -average std = %.4f -max std = %.4f -update time (s) = %.4f -prediction time (s) = %.4f -\n""" % (training_set_size, avg_force, max_force, mae, max_err, avg_std, - max_std, update_time, prediction_time) - - write_file(txt_name, txt) - - -def data_sweep_update(data_file, cell, training_snaps, cutoffs, kernel, - kernel_grad, hyps, test_snaps, txt_name): - """Uses coordinates and forces from a Quantum Espresso AIMD output file to - test the performance of a GP model as a function of the training set size. - """ - - init_time = time.time() - md_trajectory = analyze_gp.MDAnalysis(data_file, cell) - - # set up text file - txt = update_init() - write_file(txt_name, txt) - - # initialize gp model - gp_model = gp.GaussianProcess(kernel, kernel_grad, hyps, - cutoffs, opt_algorithm='BFGS') - - for n, snap in enumerate(training_snaps): - # update gp - structure = md_trajectory.get_structure_from_snap(snap) - forces = md_trajectory.get_forces_from_snap(snap) - gp_model.update_db(structure, forces) - - if n == 0: - gp_model.set_L_alpha() - else: - gp_model.update_L_alpha_v1() - - # predict on test set - all_predictions, all_variances, all_forces = \ - analyze_gp.predict_forces_on_test_set(gp_model, md_trajectory, - test_snaps, cutoffs) - - # compute and report error - training_set_size = len(gp_model.training_labels_np) - avg_force = np.mean(np.abs(all_forces)) - max_force = np.max(np.abs(all_forces)) - mae = np.mean(np.abs(all_predictions - all_forces)) - max_err = np.max(np.abs(all_predictions - all_forces)) - avg_std = np.mean(np.sqrt(all_variances)) - max_std = np.max(np.sqrt(all_variances)) - time_curr = time.time() - time_from_start = time_curr - init_time - - txt += """\n -training_set_size = %i -average force = %.4f -max force = %.4f -mean absolute error = %.4f -max error = %.4f -average std = %.4f -max std = %.4f -time from start (s) = %.4f -\n""" % (training_set_size, avg_force, max_force, mae, max_err, avg_std, - max_std, time_from_start) - - write_file(txt_name, txt) - - -def write_file(fname, text): - with open(fname, 'w') as fin: - fin.write(text) - - -def update_init(): - init_text = """Data sweep. -Date and time: %s. -Author: Jonathan Vandermause. -""" % str(datetime.datetime.now()) - - return init_text - - -def update_fin(): - fin_text = """ -------------------------------------------------------------------------------- - JOB DONE. -------------------------------------------------------------------------------- -""" - return fin_text diff --git a/flare/modules/flare_parsers/__init__.py b/flare/modules/flare_parsers/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/flare/modules/flare_parsers/add_all_parser.py b/flare/modules/flare_parsers/add_all_parser.py deleted file mode 100644 index db2ee3908..000000000 --- a/flare/modules/flare_parsers/add_all_parser.py +++ /dev/null @@ -1,397 +0,0 @@ -import sys -import numpy as np -from typing import List, Tuple -sys.path.append('../../otf/otf_engine') -import gp, env, struc, kernels, otf - - -class OtfAnalysis: - def __init__(self, filename, calculate_energy=False): - self.filename = filename - - self.calculate_energy = calculate_energy - - self.header = parse_header_information(filename) - - position_list, force_list, uncertainty_list, velocity_list,\ - dft_frames, temperatures, times, msds, dft_times, energies = \ - self.parse_pos_otf(filename) - - self.position_list = position_list - self.force_list = force_list - self.uncertainty_list = uncertainty_list - self.velocity_list = velocity_list - self.dft_frames = dft_frames - self.temperatures = temperatures - self.times = times - self.msds = msds - self.dft_times = dft_times - - if self.calculate_energy: - self.energies = energies - - gp_position_list, gp_force_list, gp_uncertainty_list,\ - gp_velocity_list, gp_hyp_list, \ - gp_species_list = self.extract_gp_info(filename) - - self.gp_position_list = gp_position_list - self.gp_force_list = gp_force_list - self.gp_uncertainty_list = gp_uncertainty_list - self.gp_velocity_list = gp_velocity_list - self.gp_hyp_list = gp_hyp_list - self.gp_species_list = gp_species_list - - def make_gp(self, cell=None, kernel=None, kernel_grad=None, algo=None, - call_no=None, cutoffs=None, hyps=None, init_gp=None): - - if init_gp is None: - # Use run's values as extracted from header - # TODO Allow for kernel gradient in header - if cell is None: - cell = self.header['cell'] - if kernel is None: - kernel = self.header['kernel'] - if kernel_grad is None: - raise Exception('Kernel gradient not supplied') - if algo is None: - algo = self.header['algo'] - if cutoffs is None: - cutoffs = self.header['cutoffs'] - if call_no is None: - call_no = len(self.gp_position_list) - if hyps is None: - gp_hyps = self.gp_hyp_list[call_no-1][-1] - else: - gp_hyps = hyps - - gp_model = gp.GaussianProcess(kernel, kernel_grad, gp_hyps, - cutoffs, opt_algorithm=algo) - else: - gp_model = init_gp - call_no = len(self.gp_position_list) - gp_hyps = self.gp_hyp_list[call_no-1][-1] - gp_model.hyps = gp_hyps - - for (positions, forces, atoms, _, species) in \ - zip(self.gp_position_list[:call_no], - self.gp_force_list[:call_no], - self.gp_atom_list[:call_no], self.gp_hyp_list[:call_no], - self.gp_species_list[:call_no]): - - struc_curr = struc.Structure(cell, species, positions) - - gp_model.update_db(struc_curr, forces, custom_range=atoms) - - gp_model.set_L_alpha() - - return gp_model - - @staticmethod - def get_gp_activation(gp_model): - pass - - def parse_pos_otf(self, filename): - """ - Exclusively parses MD run information - :param filename: - :return: - """ - position_list = [] - force_list = [] - uncertainty_list = [] - velocity_list = [] - temperatures = [] - dft_frames = [] - dft_times = [] - times = [] - msds = [] - energies = [] - - with open(filename, 'r') as f: - lines = f.readlines() - - n_steps = 0 - - for index, line in enumerate(lines): - if line.startswith("number of atoms"): - at_line = line.split() - noa = int(at_line[3]) - - # number of hyperparameters - if line.startswith("number of hyperparameters"): - line_curr = line.split(':') - noh = int(line_curr[-1]) - - # DFT frame - if line.startswith("*-Frame"): - dft_frame_line = line.split() - dft_frames.append(int(dft_frame_line[1])) - dft_time_line = lines[index+1].split() - dft_times.append(float(dft_time_line[-2])) - - # MD frame - if line.startswith("-Frame"): - n_steps += 1 - time_line = lines[index+1].split() - sim_time = float(time_line[2]) - times.append(sim_time) - - _, positions, forces, uncertainties, velocities = \ - parse_snapshot(lines, index, noa, False, noh) - - position_list.append(positions) - force_list.append(forces) - uncertainty_list.append(uncertainties) - velocity_list.append(velocities) - - temp_line = lines[index+4+noa].split() - temperatures.append(float(temp_line[-2])) - - if self.calculate_energy: - en_line = lines[index+5+noa].split() - energies.append(float(en_line[2])) - - msds.append(np.mean((positions - position_list[0])**2)) - - return position_list, force_list, uncertainty_list, velocity_list,\ - dft_frames, temperatures, times, msds, dft_times, energies - - def extract_gp_info(self, filename): - """ - Exclusively parses DFT run information - :param filename: - :return: - """ - species_list = [] - position_list = [] - force_list = [] - uncertainty_list = [] - velocity_list = [] - hyp_list = [] - - with open(filename, 'r') as f: - lines = f.readlines() - - for index, line in enumerate(lines): - - # number of atoms - if line.startswith("number of atoms"): - line_curr = line.split() - noa = int(line_curr[3]) - - # number of hyperparameters - if line.startswith("number of hyperparameters"): - line_curr = line.split() - noh = int(line_curr[3]) - - if line.startswith("*-"): - # keep track of atoms added to training set - ind_count = index+1 - line_check = lines[ind_count] - hyps_added = [] - while not line_check.startswith("-"): - - # keep track of hyperparameters - if line_check.startswith("GP hyperparameters:"): - hyps = [] - for hyp_line in lines[(ind_count+1):(ind_count+1+noh)]: - hyp_line = hyp_line.split() - hyps.append(float(hyp_line[-1])) - hyps = np.array(hyps) - hyps_added.append(hyps) - - ind_count += 1 - try: - line_check = lines[ind_count] - except: - break - # catch case where final frame is a DFT call - if line_check.startswith('Run complete'): - break - - hyp_list.append(hyps_added) - - # add DFT positions and forces - line_curr = line.split() - - # TODO: generalize this to account for arbitrary starting list - append_atom_lists(species_list, position_list, force_list, - uncertainty_list, velocity_list, - lines, index, noa, True, noh) - - return position_list, force_list, uncertainty_list, velocity_list,\ - hyp_list, species_list - - def output_md_structures(self): - """ - Returns structure objects corresponding to the MD frames of an OTF run. - :param filename: - :return: - """ - - positions = self.position_list - structures = [] - cell = self.header['cell'] - species = self.header['species'] - forces = self.force_list - stds = self.uncertainty_list - for i in range(len(positions)): - cur_struc = struc.Structure(cell=cell, species=species, - positions=positions[i]) - cur_struc.forces = forces[i] - cur_struc.stds = stds[i] - structures.append(cur_struc) - return structures - - -def append_atom_lists(species_list: List[str], - position_list: List[np.ndarray], - force_list: List[np.ndarray], - uncertainty_list: List[np.ndarray], - velocity_list: List[np.ndarray], - lines: List[str], index: int, noa: int, - dft_call: bool, noh: int) -> None: - - """Update lists containing atom information at each snapshot.""" - - species, positions, forces, uncertainties, velocities = \ - parse_snapshot(lines, index, noa, dft_call, noh) - - species_list.append(species) - position_list.append(positions) - force_list.append(forces) - uncertainty_list.append(uncertainties) - velocity_list.append(velocities) - - -def parse_snapshot(lines, index, noa, dft_call, noh): - """Parses snapshot of otf output file.""" - - # initialize values - species = [] - positions = np.zeros((noa, 3)) - forces = np.zeros((noa, 3)) - uncertainties = np.zeros((noa, 3)) - velocities = np.zeros((noa, 3)) - - # Current setting for # of lines to skip after Frame marker - skip = 3 - - for count, frame_line in enumerate(lines[(index+skip):(index+skip+noa)]): - # parse frame line - spec, position, force, uncertainty, velocity = \ - parse_frame_line(frame_line) - - # update values - species.append(spec) - positions[count] = position - forces[count] = force - uncertainties[count] = uncertainty - velocities[count] = velocity - - return species, positions, forces, uncertainties, velocities - - -def strip_and_split(line): - """ - Helper function which saves a few lines of code elsewhere - :param line: - :return: - """ - - line = line.strip().split() - stripped_line = [subline.strip() for subline in line] - - return stripped_line - - -def parse_header_information(outfile: str = 'otf_run.out') -> dict: - """ - Get information about the run from the header of the file - :param outfile: - :return: - """ - with open(outfile, 'r') as f: - lines = f.readlines() - - header_info = {} - - stopreading = None - - for line in lines: - if '---' in line or '=' in line: - stopreading = lines.index(line) - break - - if stopreading is None: - raise Exception("OTF File is malformed") - - for i, line in enumerate(lines[:stopreading]): - # TODO Update this in full - if 'cutoffs' in line: - line = line.split(':')[1].strip() - line = line.strip('[').strip(']') - line = line.split() - cutoffs = [] - for val in line: - cutoffs.append(float(val)) - header_info['cutoffs'] = cutoffs - if 'frames' in line: - header_info['frames'] = int(line.split(':')[1]) - if 'kernel' in line: - header_info['kernel'] = line.split(':')[1].strip() - if 'number of hyperparameters:' in line: - header_info['n_hyps'] = int(line.split(':')[1]) - if 'optimization algorithm' in line: - header_info['algo'] = line.split(':')[1].strip() - if 'number of atoms' in line: - header_info['atoms'] = int(line.split(':')[1]) - if 'timestep' in line: - header_info['dt'] = float(line.split(':')[1]) - if 'system species' in line: - line = line.split(':')[1] - line = line.split("'") - - species = [item for item in line if item.isalpha()] - - header_info['species_set'] = set(species) - if 'periodic cell' in line: - vectors = [] - for cell_line in lines[i+1:i+4]: - cell_line = cell_line.strip().replace('[', '').replace(']', '') - vec = cell_line.split() - vector = [float(vec[0]), float(vec[1]), float(vec[2])] - vectors.append(vector) - header_info['cell'] = np.array(vectors) - if 'previous positions' in line: - struc_spec = [] - prev_positions = [] - for pos_line in lines[i+1:i+1+header_info.get('atoms', 0)]: - pos = pos_line.split() - struc_spec.append(pos[0]) - prev_positions.append((float(pos[1]), float(pos[2]), - float(pos[3]))) - header_info['species'] = struc_spec - header_info['prev_positions'] = np.array(prev_positions) - - return header_info - - -def parse_frame_line(frame_line): - """parse a line in otf output. - - :param frame_line: frame line to be parsed - :type frame_line: string - :return: species, position, force, uncertainty, and velocity of atom - :rtype: list, np.arrays - """ - - frame_line = frame_line.split() - - spec = str(frame_line[0]) - position = np.array([float(n) for n in frame_line[1:4]]) - force = np.array([float(n) for n in frame_line[4:7]]) - uncertainty = np.array([float(n) for n in frame_line[7:10]]) - velocity = np.array([float(n) for n in frame_line[10:13]]) - - return spec, position, force, uncertainty, velocity diff --git a/flare/modules/flare_parsers/gp_from_aimd_parser.py b/flare/modules/flare_parsers/gp_from_aimd_parser.py deleted file mode 100644 index 8d8638ce2..000000000 --- a/flare/modules/flare_parsers/gp_from_aimd_parser.py +++ /dev/null @@ -1,296 +0,0 @@ -import sys -import numpy as np -from typing import List, Tuple -from flare import gp, env, struc, kernels, otf - - -class OtfAnalysis: - def __init__(self, filename, calculate_energy=False): - self.filename = filename - self.calculate_energy = calculate_energy - self.header = parse_header_information(filename) - - position_list, species_list, gp_force_list, dft_force_list, \ - gp_uncertainty_list, update_frames, msds, energies, \ - atom_list, hyp_list = \ - self.parse_file(filename) - - self.position_list = position_list - self.species_list = species_list - self.gp_force_list = gp_force_list - self.dft_force_list = dft_force_list - self.gp_uncertainty_list = gp_uncertainty_list - self.update_frames = update_frames - self.msds = msds - self.atom_list = atom_list - self.hyp_list = hyp_list - - # assumption: species list is the same across AIMD simulation - self.species_labels = species_list[0] - _, coded_species = struc.get_unique_species(self.species_labels) - self.species = coded_species - - if self.calculate_energy: - self.energies = energies - - def make_gp(self, cell=None, kernel=None, kernel_grad=None, algo=None, - call_no=None, cutoffs=None, hyps=None, init_gp=None, - energy_force_kernel=None): - - if init_gp is None: - # Use run's values as extracted from header - # TODO Allow for kernel gradient in header - if cell is None: - cell = self.header['cell'] - if kernel is None: - kernel = self.header['kernel'] - if kernel_grad is None: - raise Exception('Kernel gradient not supplied') - if algo is None: - algo = self.header['algo'] - if cutoffs is None: - cutoffs = self.header['cutoffs'] - if call_no is None: - call_no = len(self.update_frames) - if hyps is None: - gp_hyps = self.hyp_list[call_no-1] - else: - gp_hyps = hyps - - gp_model = \ - gp.GaussianProcess(kernel, kernel_grad, gp_hyps, - cutoffs, opt_algorithm=algo, - energy_force_kernel=energy_force_kernel) - else: - gp_model = init_gp - call_no = len(self.update_frames) - gp_hyps = self.hyp_list[call_no-1][-1] - gp_model.hyps = gp_hyps - - for n, frame in enumerate(self.update_frames): - positions = self.position_list[frame] - struc_curr = \ - struc.Structure(cell, self.species, positions, - species_labels=self.species_labels) - atoms = self.atom_list[n] - forces = self.dft_force_list[frame] - gp_model.update_db(struc_curr, forces, custom_range=atoms) - - gp_model.set_L_alpha() - - return gp_model - - def parse_file(self, filename): - """ - Exclusively parses MD run information - :param filename: - :return: - """ - position_list = [] - species_list = [] - gp_force_list = [] - dft_force_list = [] - gp_uncertainty_list = [] - update_frames = [] - msds = [] - energies = [] - atom_list = [] - hyp_list = [] - - with open(filename, 'r') as f: - lines = f.readlines() - - for index, line in enumerate(lines): - if line.startswith("number of atoms"): - at_line = line.split() - noa = int(at_line[3]) - - # number of hyperparameters - if line.startswith("number of hyperparameters"): - line_curr = line.split(':') - noh = int(line_curr[-1]) - - # DFT forces - if line.startswith("*-Frame"): - species, positions, forces, _ = \ - parse_snapshot(lines, index, noa) - - position_list.append(positions) - species_list.append(species) - msds.append(np.mean((positions - position_list[0])**2)) - dft_force_list.append(forces) - - frame_line = line.split() - frame_curr = (int(frame_line[1])) - - # GP forces - if line.startswith("-Frame"): - _, _, forces, uncertainties = \ - parse_snapshot(lines, index, noa) - - gp_force_list.append(forces) - gp_uncertainty_list.append(uncertainties) - - if self.calculate_energy: - en_line = lines[index+5+noa].split() - energies.append(float(en_line[2])) - - # Update frame - if line.startswith("Adding atom"): - # keep track of atoms added to training set - atoms_added = [] - line_split = line.split() - atom_strings = line_split[2:-4] - for n, atom_string in enumerate(atom_strings): - if n == 0: - atoms_added.append(int(atom_string[1:-1])) - else: - atoms_added.append(int(atom_string[0:-1])) - atom_list.append(atoms_added) - update_frames.append(frame_curr) - - # Hyperparameter block - if line.startswith("GP hyperparameters:"): - hyps = [] - for hyp_line in lines[(index+1):(index+1+noh)]: - hyp_line = hyp_line.split() - hyps.append(float(hyp_line[-1])) - hyps = np.array(hyps) - hyp_list.append(hyps) - - return position_list, species_list, gp_force_list, dft_force_list, \ - gp_uncertainty_list, update_frames, msds, energies, \ - atom_list, hyp_list - - -def parse_snapshot(lines, index, noa): - """Parses snapshot of otf output file.""" - - # initialize values - species = [] - positions = np.zeros((noa, 3)) - forces = np.zeros((noa, 3)) - uncertainties = np.zeros((noa, 3)) - - # Current setting for # of lines to skip after Frame marker - skip = 3 - - for count, frame_line in enumerate(lines[(index+skip):(index+skip+noa)]): - # parse frame line - spec, position, force, uncertainty = \ - parse_frame_line(frame_line) - - # update values - species.append(spec) - positions[count] = position - forces[count] = force - uncertainties[count] = uncertainty - - return species, positions, forces, uncertainties - - -def strip_and_split(line): - """ - Helper function which saves a few lines of code elsewhere - :param line: - :return: - """ - - line = line.strip().split() - stripped_line = [subline.strip() for subline in line] - - return stripped_line - - -def parse_header_information(outfile: str = 'otf_run.out') -> dict: - """ - Get information about the run from the header of the file - :param outfile: - :return: - """ - with open(outfile, 'r') as f: - lines = f.readlines() - - header_info = {} - - stopreading = None - - for line in lines: - if '---' in line or '=' in line: - stopreading = lines.index(line) - break - - if stopreading is None: - raise Exception("OTF File is malformed") - - for i, line in enumerate(lines[:stopreading]): - # TODO Update this in full - if 'cutoffs' in line: - line = line.split(':')[1].strip() - line = line.strip('[').strip(']') - line = line.split() - cutoffs = [] - for val in line: - try: - cutoffs.append(float(val)) - except: - cutoffs.append(float(val[:-1])) - header_info['cutoffs'] = cutoffs - if 'frames' in line: - header_info['frames'] = int(line.split(':')[1]) - if 'kernel' in line: - header_info['kernel'] = line.split(':')[1].strip() - if 'number of hyperparameters:' in line: - header_info['n_hyps'] = int(line.split(':')[1]) - if 'optimization algorithm' in line: - header_info['algo'] = line.split(':')[1].strip() - if 'number of atoms' in line: - header_info['atoms'] = int(line.split(':')[1]) - if 'timestep' in line: - header_info['dt'] = float(line.split(':')[1]) - if 'system species' in line: - line = line.split(':')[1] - line = line.split("'") - - species = [item for item in line if item.isalpha()] - - header_info['species_set'] = set(species) - if 'periodic cell' in line: - vectors = [] - for cell_line in lines[i+1:i+4]: - cell_line = \ - cell_line.strip().replace('[', '').replace(']', '') - vec = cell_line.split() - vector = [float(vec[0]), float(vec[1]), float(vec[2])] - vectors.append(vector) - header_info['cell'] = np.array(vectors) - if 'previous positions' in line: - struc_spec = [] - prev_positions = [] - for pos_line in lines[i+1:i+1+header_info.get('atoms', 0)]: - pos = pos_line.split() - struc_spec.append(pos[0]) - prev_positions.append((float(pos[1]), float(pos[2]), - float(pos[3]))) - header_info['species'] = struc_spec - header_info['prev_positions'] = np.array(prev_positions) - - return header_info - - -def parse_frame_line(frame_line): - """parse a line in otf output. - :param frame_line: frame line to be parsed - :type frame_line: string - :return: species, position, force, uncertainty, and velocity of atom - :rtype: list, np.arrays - """ - - frame_line = frame_line.split() - - spec = str(frame_line[0]) - position = np.array([float(n) for n in frame_line[1:4]]) - force = np.array([float(n) for n in frame_line[4:7]]) - uncertainty = np.array([float(n) for n in frame_line[7:10]]) - - return spec, position, force, uncertainty diff --git a/flare/modules/flare_parsers/old_vac_parser.py b/flare/modules/flare_parsers/old_vac_parser.py deleted file mode 100644 index d3291c029..000000000 --- a/flare/modules/flare_parsers/old_vac_parser.py +++ /dev/null @@ -1,218 +0,0 @@ -import sys -import numpy as np -import crystals -sys.path.append('../otf_engine') -import gp, env, struc, kernels, otf - - -class OtfAnalysis: - def __init__(self, filename): - self.filename = filename - - positions, dft_frames, temperatures, times, msds, dft_times = \ - self.parse_pos_otf(filename) - - self.positions = positions - self.dft_frames = dft_frames - self.temperatures = temperatures - self.times = times - self.msds = msds - self.dft_times = dft_times - - gp_pos_list, gp_force_list, gp_atom_list, gp_hyp_list, \ - gp_cutoff_radius, gp_species_list = self.extract_gp_info(filename) - - self.gp_pos_list = gp_pos_list - self.gp_force_list = gp_force_list - self.gp_atom_list = gp_atom_list - self.gp_hyp_list = gp_hyp_list - self.gp_cutoff_radius = gp_cutoff_radius - self.gp_species_list = gp_species_list - - def make_gp(self, cell, kernel, bodies, algo, call_no, cutoffs=None): - gp_model = gp.GaussianProcess(kernel, bodies, algo, cutoffs=cutoffs) - gp_hyps = self.gp_hyp_list[call_no-1] - gp_model.hyps = gp_hyps - - for count, (positions, forces, atom, _, species) in \ - enumerate(zip(self.gp_pos_list, self.gp_force_list, - self.gp_atom_list, self.gp_hyp_list, - self.gp_species_list)): - if count < call_no: - struc_curr = struc.Structure(cell, species, positions, - self.gp_cutoff_radius) - - gp_model.update_db(struc_curr, forces, custom_range=[atom]) - - gp_model.set_L_alpha() - - return gp_model - - @staticmethod - def get_gp_activation(gp_model): - pass - - def parse_pos_otf(self, filename): - positions = [] - temperatures = [] - dft_frames = [] - dft_times = [] - times = [] - msds = [] - - with open(filename, 'r') as f: - lines = f.readlines() - - n_steps = 0 - - for index, line in enumerate(lines): - if line.startswith("Number of Atoms"): - at_line = line.split() - noa = int(at_line[3]) - - if line.startswith("-*"): - dft_line = line.split() - dft_frames.append(int(dft_line[1])) - dft_times.append(float(dft_line[4])) - - if line.startswith("- Frame"): - n_steps += 1 - pos = [] - - frame_line = line.split() - sim_time = float(frame_line[5]) - times.append(sim_time) - - for frame_line in lines[(index+2):(index+2+noa)]: - frame_line = frame_line.split() - curr_pos = np.zeros(shape=(3,)) - curr_pos[0] = str(frame_line[1]) - curr_pos[1] = str(frame_line[2]) - curr_pos[2] = str(frame_line[3]) - - pos.append(curr_pos) - - temp_line = lines[index+2+noa].split() - temperatures.append(float(temp_line[1])) - - positions.append(np.array(pos)) - - msds.append(np.mean((pos - positions[0])**2)) - - return positions, dft_frames, temperatures, times, msds, dft_times - - def extract_gp_info(self, filename): - pos_list = [] - species_list = [] - atom_list = [] - force_list = [] - hyp_list = [] - - with open(filename, 'r') as f: - lines = f.readlines() - - for index, line in enumerate(lines): - if line.startswith("Structure Cutoff Radius:"): - line_curr = line.split() - cutoff_radius = float(line_curr[3]) - - if line.startswith("Number of Atoms"): - line_curr = line.split() - noa = int(line_curr[3]) - - if line.startswith("# Hyperparameters"): - line_curr = line.split() - noh = int(line_curr[2]) - - if line.startswith("Calling DFT with training atoms"): - line_curr = line.split() - - pos = [] - fcs = [] - hyps = [] - species = [] - - for frame_line in lines[(index+4):(index+4+noh)]: - frame_line = frame_line.split() - hyps.append(float(frame_line[5])) - hyps = np.array(hyps) - hyp_list.append(hyps) - hyp_list.append(hyps) - - for frame_line in lines[(index+12):(index+12+noa)]: - frame_line = frame_line.split() - curr_pos = np.zeros(shape=(3,)) - curr_pos[0] = str(frame_line[1]) - curr_pos[1] = str(frame_line[2]) - curr_pos[2] = str(frame_line[3]) - - dft_fc = np.zeros(shape=(3,)) - dft_fc[0] = str(frame_line[4]) - dft_fc[1] = str(frame_line[5]) - dft_fc[2] = str(frame_line[6]) - - species.append(str(frame_line[0])) - pos.append(curr_pos) - fcs.append(dft_fc) - - pos = np.array(pos) - pos_list.append(pos) - pos_list.append(pos) - - species_list.append(species) - species_list.append(species) - - fcs = np.array(fcs) - force_list.append(fcs) - force_list.append(fcs) - - atom_list.append(0) - atom_list.append(30) - - if line.startswith("Calling DFT due to"): - line_curr = line.split() - - pos = [] - fcs = [] - hyps = [] - species = [] - - for frame_line in lines[(index+4):(index+4+noh)]: - frame_line = frame_line.split() - hyps.append(float(frame_line[5])) - hyps = np.array(hyps) - hyp_list.append(hyps) - - for frame_line in lines[(index+12):(index+12+noa)]: - frame_line = frame_line.split() - curr_pos = np.zeros(shape=(3,)) - curr_pos[0] = str(frame_line[1]) - curr_pos[1] = str(frame_line[2]) - curr_pos[2] = str(frame_line[3]) - - dft_fc = np.zeros(shape=(3,)) - dft_fc[0] = str(frame_line[4]) - dft_fc[1] = str(frame_line[5]) - dft_fc[2] = str(frame_line[6]) - - pos.append(curr_pos) - fcs.append(dft_fc) - - species.append(str(frame_line[0])) - - pos = np.array(pos) - pos_list.append(pos) - - fcs = np.array(fcs) - force_list.append(fcs) - - species_list.append(species) - - atom_list.append(int(line_curr[5])) - - return pos_list, force_list, atom_list, hyp_list, cutoff_radius,\ - species_list - - -def get_gp_from_file(filename, cell): - pass diff --git a/flare/modules/flare_parsers/otf_parser_v0.py b/flare/modules/flare_parsers/otf_parser_v0.py deleted file mode 100644 index 7c44afa05..000000000 --- a/flare/modules/flare_parsers/otf_parser_v0.py +++ /dev/null @@ -1,267 +0,0 @@ -import sys -import numpy as np -from flare.modules import crystals -from typing import List, Tuple -from flare import gp, env, struc, kernels, otf - - -class OtfAnalysis: - def __init__(self, filename, calculate_energy=False): - self.filename = filename - - self.calculate_energy = calculate_energy - - position_list, force_list, uncertainty_list, velocity_list,\ - dft_frames, temperatures, times, msds, dft_times, energies = \ - self.parse_pos_otf(filename) - - self.position_list = position_list - self.force_list = force_list - self.uncertainty_list = uncertainty_list - self.velocity_list = velocity_list - self.dft_frames = dft_frames - self.temperatures = temperatures - self.times = times - self.msds = msds - self.dft_times = dft_times - - if self.calculate_energy: - self.energies = energies - - gp_position_list, gp_force_list, gp_uncertainty_list,\ - gp_velocity_list, gp_atom_list, gp_hyp_list, \ - gp_species_list = self.extract_gp_info(filename) - - self.gp_position_list = gp_position_list - self.gp_force_list = gp_force_list - self.gp_uncertainty_list = gp_uncertainty_list - self.gp_velocity_list = gp_velocity_list - self.gp_atom_list = gp_atom_list - self.gp_hyp_list = gp_hyp_list - self.gp_species_list = gp_species_list - - def make_gp(self, cell, kernel, kernel_grad, algo, call_no, - start_list, cutoffs): - gp_hyps = self.gp_hyp_list[call_no-1] - gp_model = gp.GaussianProcess(kernel, kernel_grad, gp_hyps, - cutoffs) - - for count, (positions, forces, atom, _, species) in \ - enumerate(zip(self.gp_position_list, self.gp_force_list, - self.gp_atom_list, self.gp_hyp_list, - self.gp_species_list)): - - struc_curr = struc.Structure(cell, species, positions) - - # add atoms in start list - if count == 0: - gp_model.update_db(struc_curr, forces, - custom_range=start_list) - - # add one atom - elif count < call_no: - gp_model.update_db(struc_curr, forces, custom_range=[atom]) - - gp_model.set_L_alpha() - - return gp_model - - @staticmethod - def get_gp_activation(gp_model): - pass - - def parse_pos_otf(self, filename): - position_list = [] - force_list = [] - uncertainty_list = [] - velocity_list = [] - temperatures = [] - dft_frames = [] - dft_times = [] - times = [] - msds = [] - energies = [] - - with open(filename, 'r') as f: - lines = f.readlines() - - n_steps = 0 - - for index, line in enumerate(lines): - if line.startswith("Number of Atoms"): - at_line = line.split() - noa = int(at_line[3]) - - # number of hyperparameters - if line.startswith("# Hyperparameters"): - line_curr = line.split() - noh = int(line_curr[2]) - - if line.startswith("-*"): - dft_line = line.split() - dft_frames.append(int(dft_line[1])) - dft_times.append(float(dft_line[4])) - - if line.startswith("- Frame"): - n_steps += 1 - frame_line = line.split() - sim_time = float(frame_line[5]) - times.append(sim_time) - - _, positions, forces, uncertainties, velocities = \ - parse_snapshot(lines, index, noa, False, noh) - - position_list.append(positions) - force_list.append(forces) - uncertainty_list.append(uncertainties) - velocity_list.append(velocities) - - temp_line = lines[index+2+noa].split() - temperatures.append(float(temp_line[1])) - - if self.calculate_energy: - en_line = lines[index+5+noa].split() - energies.append(float(en_line[2])) - - msds.append(np.mean((positions - position_list[0])**2)) - - return position_list, force_list, uncertainty_list, velocity_list,\ - dft_frames, temperatures, times, msds, dft_times, energies - - def extract_gp_info(self, filename): - species_list = [] - position_list = [] - force_list = [] - uncertainty_list = [] - velocity_list = [] - atom_list = [] - hyp_list = [] - - with open(filename, 'r') as f: - lines = f.readlines() - - for index, line in enumerate(lines): - # cutoff radius - # if line.startswith("Structure Cutoff Radius:"): - # line_curr = line.split() - # cutoff_radius = float(line_curr[3]) - - # number of atoms - if line.startswith("Number of Atoms"): - line_curr = line.split() - noa = int(line_curr[3]) - - # number of hyperparameters - if line.startswith("# Hyperparameters"): - line_curr = line.split() - noh = int(line_curr[2]) - - # TODO: change otf program so that this if statement is unnecessary - if line.startswith("Calling DFT with training atoms"): - line_curr = line.split() - - # TODO: write function for updating hyps list - hyps = [] - for frame_line in lines[(index+4):(index+4+noh)]: - frame_line = frame_line.split() - hyps.append(float(frame_line[4])) - hyps = np.array(hyps) - hyp_list.append(hyps) - - # TODO: generalize this to account for arbitrary starting list - append_atom_lists(species_list, position_list, force_list, - uncertainty_list, velocity_list, - lines, index, noa, True, noh) - - atom_list.append(0) - - if line.startswith("Calling DFT due to"): - line_curr = line.split() - - # TODO: write function for updating hyps list - hyps = [] - for frame_line in lines[(index+4):(index+4+noh)]: - frame_line = frame_line.split() - hyps.append(float(frame_line[4])) - hyps = np.array(hyps) - hyp_list.append(hyps) - - # TODO: generalize atom list update to account for arbitrary - # list - append_atom_lists(species_list, position_list, force_list, - uncertainty_list, velocity_list, - lines, index, noa, True, noh) - atom_list.append(int(line_curr[5])) - - return position_list, force_list, uncertainty_list, velocity_list,\ - atom_list, hyp_list, species_list - - -def append_atom_lists(species_list: List[str], - position_list: List[np.ndarray], - force_list: List[np.ndarray], - uncertainty_list: List[np.ndarray], - velocity_list: List[np.ndarray], - lines: List[str], index: int, noa: int, - dft_call: bool, noh: int) -> None: - - """Update lists containing atom information at each snapshot.""" - - species, positions, forces, uncertainties, velocities = \ - parse_snapshot(lines, index, noa, dft_call, noh) - - species_list.append(species) - position_list.append(positions) - force_list.append(forces) - uncertainty_list.append(uncertainties) - velocity_list.append(velocities) - - -def parse_snapshot(lines, index, noa, dft_call, noh): - """Parses snapshot of otf output file.""" - - # initialize values - species = [] - positions = np.zeros((noa, 3)) - forces = np.zeros((noa, 3)) - uncertainties = np.zeros((noa, 3)) - velocities = np.zeros((noa, 3)) - - # number of lines to skip after frame line - if dft_call is True: - skip = 9 + noh - else: - skip = 2 - - for count, frame_line in enumerate(lines[(index+skip):(index+skip+noa)]): - # parse frame line - spec, position, force, uncertainty, velocity = \ - parse_frame_line(frame_line) - - # update values - species.append(spec) - positions[count] = position - forces[count] = force - uncertainties[count] = uncertainty - velocities[count] = velocity - - return species, positions, forces, uncertainties, velocities - - -def parse_frame_line(frame_line): - """parse a line in otf output. - :param frame_line: frame line to be parsed - :type frame_line: string - :return: species, position, force, uncertainty, and velocity of atom - :rtype: list, np.arrays - """ - - frame_line = frame_line.split() - - spec = str(frame_line[0]) - position = np.array([float(n) for n in frame_line[1:4]]) - force = np.array([float(n) for n in frame_line[4:7]]) - uncertainty = np.array([float(n) for n in frame_line[7:10]]) - velocity = np.array([float(n) for n in frame_line[10:13]]) - - return spec, position, force, uncertainty, velocity \ No newline at end of file diff --git a/flare/modules/flare_parsers/parse_output_temp.py b/flare/modules/flare_parsers/parse_output_temp.py deleted file mode 100644 index baef652e1..000000000 --- a/flare/modules/flare_parsers/parse_output_temp.py +++ /dev/null @@ -1,106 +0,0 @@ -import sys -import numpy as np - -from matplotlib import pyplot as plt - - -def read_pos(filename): - pos = [] - - with open(filename, 'r') as outf: - lines = outf.readlines() - - for line in lines: - line = line.split() - - curr_pos = np.zeros(shape=(3,)) - curr_pos[0] = str(line[1]) - curr_pos[1] = str(line[2]) - curr_pos[2] = str(line[3]) - pos.append(curr_pos) - - return np.array(pos) - - -def get_v(pos, prev_pos, next_pos): - v = np.zeros(shape=pos.shape) - - for i in range(pos.shape[0]): - v[i, :] = (next_pos[i, :] - prev_pos[i, :]) / (2 * dt) - - return v - - -def get_temp(pos, prev_pos, next_pos): - prev_positions_m = prev_pos * angst2m - positions_m = pos * angst2m - next_positions_m = next_pos * angst2m - - v = get_v(pos=positions_m, prev_pos=prev_positions_m, next_pos=next_positions_m) - v_magn = np.empty(v.shape[0]) - - for i in range(v_magn.shape[0]): - v_magn[i] = np.sqrt(v[i, 0]**2 + v[i, 1]**2 + v[i, 2]**2) - - T = (1/((3 * n_atoms - 3) * k_b)) * np.sum((v_magn**2) * m_al) - - return T - - -def parse_pos_otf(filename): - traj = [] - - with open(filename, 'r') as f: - lines = f.readlines() - - n_steps = 0 - - for index, line in enumerate(lines): - if line.startswith("Number of Atoms"): - at_line = line.split() - noa = int(at_line[3]) - - if line.startswith("- Frame"): - n_steps += 1 - pos = [] - - for frame_line in lines[(index+2):(index+2+noa)]: - frame_line = frame_line.split() - curr_pos = np.zeros(shape=(3,)) - curr_pos[0] = str(frame_line[1]) - curr_pos[1] = str(frame_line[2]) - curr_pos[2] = str(frame_line[3]) - - pos.append(curr_pos) - - traj.append(np.array(pos)) - - return traj - -if __name__ == "__main__": - - otf_filename = sys.argv[1] - - # params - n_atoms = 31 - n_comp = 5000 - dt = 1e-15 # sec - m_al = (26.98154 * 0.001) / (6.022140857 * 1e23) # kg/atom - k_b = 1.38064852 * 1e-23 # J/K - angst2m = 1e-10 # Angstrom to m - - # get otf temperature - temps_otf = [] - pos_traj = parse_pos_otf(filename=otf_filename) - - for i in range(len(pos_traj)-2): - temps_otf.append(get_temp(pos=pos_traj[i+1], prev_pos=pos_traj[i], next_pos=pos_traj[i+2])) - - # plot - plt.plot(list(range(len(pos_traj)-2)), temps_otf) - plt.xlabel('Frame') - plt.ylabel('Instantaneous Temperature [K]') - plt.show() - - np.savetxt('temps_otf.txt', np.array(temps_otf)) - diff --git a/flare/modules/gp_calculator.py b/flare/modules/gp_calculator.py deleted file mode 100644 index 007640ce3..000000000 --- a/flare/modules/gp_calculator.py +++ /dev/null @@ -1,59 +0,0 @@ -import numpy as np -import multiprocessing as mp -import concurrent.futures -import copy -import sys -from flare import gp, env, struc, kernels, otf -from ase.calculators.calculator import Calculator - -class GPCalculator(Calculator): - def __init__(self, gp_model): - super().__init__() # all set to default values,TODO: change - self.gp_model = gp_model - - def get_potential_energy(self, atoms=None, force_consistent=False): - return 0 - nat = len(atoms) - species = atoms.get_chemical_symbols() - - struc_curr = struc.Structure(atoms.cell[:], species, - atoms.positions) - local_energies = np.zeros(nat) - - for n in range(nat): - chemenv = env.AtomicEnvironment(struc_curr, n, - self.gp_model.cutoffs) - local_energies[n] = self.gp_model.predict_local_energy(chemenv) - - return np.sum(local_energies) - - def get_forces(self, atoms): - nat = len(atoms) - species = atoms.get_chemical_symbols() - - struc_curr = struc.Structure(atoms.cell[:], species, - atoms.positions) - - forces = np.zeros((nat, 3)) - stds = np.zeros((nat, 3)) - for n in range(nat): - chemenv = env.AtomicEnvironment(struc_curr, n, - self.gp_model.cutoffs) - for i in range(3): - force, std = self.gp_model.predict(chemenv, i + 1) - forces[n][i] = float(force) - stds[n][i] = np.sqrt(np.absolute(std)) - - self.results['stds'] = stds - atoms.get_uncertainties = self.get_uncertainties - - return forces - - def get_stress(self, atoms): - return np.eye(3) - - def calculation_required(self, atoms, quantities): - return True - - def get_uncertainties(self): - return self.results['stds'] diff --git a/flare/modules/gp_from_aimd.py b/flare/modules/gp_from_aimd.py deleted file mode 100644 index 614504c26..000000000 --- a/flare/modules/gp_from_aimd.py +++ /dev/null @@ -1,392 +0,0 @@ -import sys -import numpy as np -import datetime -import time -from typing import List -import copy -import multiprocessing as mp -import concurrent.futures -from flare import struc, gp, env, qe_util, md - - -class ActiveGp: - def __init__(self, position_list, force_list, species, cell, - gp: gp.GaussianProcess, - std_tolerance_factor: float = 1, par: bool=False, - init_atoms: List[int]=None, - calculate_energy=False, output_name='otf_run.out', - max_atoms_added=1, freeze_hyps=10, no_cpus=1, - zeroed_atoms=False, no_zeroed=0): - - self.position_list = position_list - self.force_list = force_list - self.species = species - self.cell = cell - self.gp = gp - self.std_tolerance = std_tolerance_factor - self.dft_step = True - self.freeze_hyps = freeze_hyps - - self.zeroed_atoms = zeroed_atoms - self.no_zeroed = no_zeroed - - positions = position_list[0] - self.structure = struc.Structure(cell=cell, species=species, - positions=positions) - - self.noa = self.structure.positions.shape[0] - self.atom_list = list(range(self.noa)) - self.curr_step = 0 - - self.dt = 1.0 - self.number_of_steps = 100 - self.temperature = 1.0 - self.KE = 1.0 - self.velocities = np.zeros((self.noa, 3)) - - self.max_atoms_added = max_atoms_added - - # initialize local energies - if calculate_energy: - self.local_energies = np.zeros(self.noa) - else: - self.local_energies = None - - # set atom list for initial dft run - if init_atoms is None: - self.init_atoms = [int(n) for n in range(self.noa)] - else: - self.init_atoms = init_atoms - - self.dft_count = 0 - - # set pred function - if not par and not calculate_energy: - self.pred_func = self.predict_on_structure - elif par and not calculate_energy: - self.pred_func = self.predict_on_structure_par - elif not par and calculate_energy: - self.pred_func = self.predict_on_structure_en - elif par and calculate_energy: - self.pred_func = self.predict_on_structure_par_en - self.par = par - - self.output_name = output_name - - # set number of cpus for qe runs - self.no_cpus = no_cpus - - def run(self): - write_header(self.gp.cutoffs, self.gp.kernel_name, self.gp.hyps, - self.gp.algo, self.dt, self.number_of_steps, - self.structure, self.output_name, - self.std_tolerance) - - self.start_time = time.time() - - for positions, forces in zip(self.position_list, self.force_list): - - # run DFT and train initial model if first step and DFT is on - if self.curr_step == 0 and self.std_tolerance != 0: - # update positions and forces - self.structure.positions = positions - self.structure.wrap_positions() - self.structure.forces = copy.deepcopy(forces) - self.dft_step = True - self.record_state() - - # make initial gp model and predict forces - self.update_gp(self.init_atoms, forces) - self.dft_count += 1 - if (self.dft_count-1) < self.freeze_hyps: - self.train_gp() - - # after step 1, try predicting with GP model - else: - # update positions and forces - self.structure.positions = positions - self.structure.wrap_positions() - self.structure.forces = copy.deepcopy(forces) - self.dft_step = True - self.record_state() - - # predict with gp - self.pred_func() - - # get max uncertainty atoms - if self.zeroed_atoms: - std_in_bound, target_atoms = self.is_std_in_bound_zeroed() - mae = np.mean(np.abs(self.structure.forces - [0:-self.zeroed_atoms] - - forces[0:-self.zeroed_atoms])) - mac = np.mean(np.abs(forces[0:-self.zeroed_atoms])) - else: - std_in_bound, target_atoms = self.is_std_in_bound() - mae = np.mean(np.abs(self.structure.forces - forces)) - mac = np.mean(np.abs(forces)) - - # record GP forces - self.dft_step = False - self.record_state() - - write_to_output('\nmae: %.4f \n' % mae, self.output_name) - write_to_output('\nmac: %.4f \n' % mac, self.output_name) - - # add max uncertainty atoms to training set - if not std_in_bound: - self.update_gp(target_atoms, forces) - self.dft_count += 1 - if (self.dft_count-1) < self.freeze_hyps: - self.train_gp() - - self.curr_step += 1 - - conclude_run(self.output_name) - - def predict_on_atom(self, atom): - chemenv = env.AtomicEnvironment(self.structure, atom, self.gp.cutoffs) - comps = [] - stds = [] - # predict force components and standard deviations - for i in range(3): - force, var = self.gp.predict(chemenv, i+1) - comps.append(float(force)) - stds.append(np.sqrt(np.abs(var))) - - return comps, stds - - def predict_on_atom_en(self, atom): - chemenv = env.AtomicEnvironment(self.structure, atom, self.gp.cutoffs) - comps = [] - stds = [] - # predict force components and standard deviations - for i in range(3): - force, var = self.gp.predict(chemenv, i+1) - comps.append(float(force)) - stds.append(np.sqrt(np.abs(var))) - - # predict local energy - local_energy = self.gp.predict_local_energy(chemenv) - return comps, stds, local_energy - - def predict_on_structure_par(self): - n = 0 - with concurrent.futures.ProcessPoolExecutor() as executor: - for res in executor.map(self.predict_on_atom, self.atom_list): - for i in range(3): - self.structure.forces[n][i] = res[0][i] - self.structure.stds[n][i] = res[1][i] - n += 1 - - def predict_on_structure_par_en(self): - n = 0 - with concurrent.futures.ProcessPoolExecutor() as executor: - for res in executor.map(self.predict_on_atom_en, self.atom_list): - for i in range(3): - self.structure.forces[n][i] = res[0][i] - self.structure.stds[n][i] = res[1][i] - self.local_energies[n] = res[2] - n += 1 - - def predict_on_structure(self): - for n in range(self.structure.nat): - chemenv = env.AtomicEnvironment(self.structure, n, self.gp.cutoffs) - for i in range(3): - force, var = self.gp.predict(chemenv, i + 1) - self.structure.forces[n][i] = float(force) - self.structure.stds[n][i] = np.sqrt(np.abs(var)) - - def predict_on_structure_en(self): - for n in range(self.structure.nat): - chemenv = env.AtomicEnvironment(self.structure, n, self.gp.cutoffs) - for i in range(3): - force, var = self.gp.predict(chemenv, i + 1) - self.structure.forces[n][i] = float(force) - self.structure.stds[n][i] = np.sqrt(np.abs(var)) - self.local_energies[n] = self.gp.predict_local_energy(chemenv) - - def update_gp(self, train_atoms, dft_frcs): - write_to_output('\nAdding atom {} to the training set.\n' - .format(train_atoms), - self.output_name) - write_to_output('Uncertainty: {}.\n' - .format(self.structure.stds[train_atoms[0]]), - self.output_name) - - # update gp model - self.gp.update_db(self.structure, dft_frcs, - custom_range=train_atoms) - - self.gp.set_L_alpha() - - def train_gp(self): - self.gp.train(True) - write_hyps(self.gp.hyp_labels, self.gp.hyps, - self.start_time, self.output_name, - self.gp.like, self.gp.like_grad) - - def is_std_in_bound(self): - # set uncertainty threshold - if self.std_tolerance == 0: - return True, -1 - elif self.std_tolerance > 0: - threshold = self.std_tolerance * np.abs(self.gp.hyps[-1]) - else: - threshold = np.abs(self.std_tolerance) - - # sort max stds - max_stds = np.zeros((self.noa)) - for atom, std in enumerate(self.structure.stds): - max_stds[atom] = np.max(std) - stds_sorted = np.argsort(max_stds) - target_atoms = list(stds_sorted[-self.max_atoms_added:]) - - # if above threshold, return atom - if max_stds[stds_sorted[-1]] > threshold: - return False, target_atoms - else: - return True, [-1] - - # zeroed atoms must appear last - def is_std_in_bound_zeroed(self): - # set uncertainty threshold - if self.std_tolerance == 0: - return True, -1 - elif self.std_tolerance > 0: - threshold = self.std_tolerance * np.abs(self.gp.hyps[-1]) - else: - threshold = np.abs(self.std_tolerance) - - # sort max stds - max_stds = np.zeros((self.noa - self.no_zeroed)) - for atom, std in enumerate(self.structure.stds[0:-self.no_zeroed]): - max_stds[atom] = np.max(std) - stds_sorted = np.argsort(max_stds) - target_atoms = list(stds_sorted[-self.max_atoms_added:]) - - # if above threshold, return atom - if max_stds[stds_sorted[-1]] > threshold: - return False, target_atoms - else: - return True, [-1] - - def record_state(self): - write_md_config(self.dt, self.curr_step, self.structure, - self.temperature, self.KE, - self.local_energies, self.start_time, - self.output_name, self.dft_step, - self.velocities) - - -def write_to_output(string: str, output_file: str = 'otf_run.out'): - with open(output_file, 'a') as f: - f.write(string) - - -def write_header(cutoffs, kernel_name, hyps, algo, dt, Nsteps, structure, - output_name, std_tolerance): - - with open(output_name, 'w') as f: - f.write(str(datetime.datetime.now()) + '\n') - - if std_tolerance < 0: - std_string = \ - 'uncertainty tolerance: {} eV/A\n'.format(np.abs(std_tolerance)) - elif std_tolerance > 0: - std_string = \ - 'uncertainty tolerance: {} times noise \n'\ - .format(np.abs(std_tolerance)) - else: - std_string = '' - - headerstring = '' - headerstring += 'cutoffs: {}\n'.format(cutoffs) - headerstring += 'kernel: {}\n'.format(kernel_name) - headerstring += 'number of hyperparameters: {}\n'.format(len(hyps)) - headerstring += 'hyperparameters: {}' \ - '\n'.format(hyps) - headerstring += 'hyperparameter optimization algorithm: {}' \ - '\n'.format(algo) - headerstring += std_string - headerstring += 'timestep (ps): {}\n'.format(dt) - headerstring += 'number of atoms: {}\n'.format(structure.nat) - headerstring += 'system species: {}\n'.format(set(structure.species)) - headerstring += 'periodic cell: \n' - headerstring += str(structure.cell) - - headerstring += '\n' - headerstring += '-' * 80 + '\n' - - write_to_output(headerstring, output_name) - - -def write_md_config(dt, curr_step, structure, temperature, KE, - local_energies, start_time, output_name, dft_step, - velocities): - string = '' - - # Mark if a frame had DFT forces with an asterisk - if not dft_step: - string += "\n-Frame: " + str(curr_step) - else: - string += '-' * 80 + '\n' - string += "\n*-Frame: " + str(curr_step) - - string += '\nSimulation Time: %.3f ps \n' % (dt * curr_step) - - # Construct Header line - string += 'El Position (A) \t\t\t\t ' - if not dft_step: - string += 'GP Force (ev/A) ' - else: - string += 'DFT Force (ev/A) ' - string += '\t\t\t\t Std. Dev (ev/A) \n' - - # Construct atom-by-atom description - for i in range(len(structure.positions)): - string += structure.species[i] + ' ' - for j in range(3): - string += str("%.8f" % structure.positions[i][j]) + ' ' - string += '\t' - for j in range(3): - string += str("%.8f" % structure.forces[i][j]) + ' ' - string += '\t' - for j in range(3): - string += str("%.8e" % structure.stds[i][j]) + ' ' - string += '\n' - - string += '\n' - - # calculate potential and total energy - if local_energies is not None: - pot_en = np.sum(local_energies) - tot_en = KE + pot_en - string += \ - 'potential energy: %.6f eV \n' % pot_en - string += 'total energy: %.6f eV \n' % tot_en - - string += 'wall time from start: %.2f s \n' % \ - (time.time() - start_time) - - write_to_output(string, output_name) - - -def write_hyps(hyp_labels, hyps, start_time, output_name, like, like_grad): - write_to_output('\nGP hyperparameters: \n', output_name) - - for i, label in enumerate(hyp_labels): - write_to_output('Hyp{} : {} = {}\n'.format(i, label, hyps[i]), - output_name) - - write_to_output('likelihood: '+str(like)+'\n', output_name) - write_to_output('likelihood gradient: '+str(like_grad)+'\n', output_name) - time_curr = time.time() - start_time - write_to_output('wall time from start: %.2f s \n' % time_curr, - output_name) - - -def conclude_run(output_name): - footer = '-' * 20 + '\n' - footer += 'Run complete. \n' - - write_to_output(footer, output_name) diff --git a/flare/modules/gp_vs_deepmd.py b/flare/modules/gp_vs_deepmd.py deleted file mode 100644 index 6d9369e7e..000000000 --- a/flare/modules/gp_vs_deepmd.py +++ /dev/null @@ -1,102 +0,0 @@ -import numpy as np -import sys -import qe_parsers, analyze_gp, eam, analyze_md, otf_parser_v2 -from scipy.optimize import minimize -import time -import copy -sys.path.append('../otf_engine') -import gp, env, struc, kernels, otf, md, md_run - - -def get_structure_from_dataset(filename, snap): - box_test = np.load(file_name + 'box.npy') - coord_test = np.load(file_name + 'coord.npy') - force_test = np.load(file_name + 'force.npy') - energy_test = np.load(file_name + 'energy.npy') - - positions = coord_test[snap].reshape(-1, 3) - forces = force_test[snap].reshape(-1, 3) - cell = box_test[snap].reshape(3, 3) - species = ['Cu'] * positions.shape[0] - - struc_curr = struc.Structure(cell, species, positions) - - return struc_curr, forces - - -def predict_forces_on_structure(gp_model, structure, forces, cutoffs): - noa = len(structure.positions) - - predictions = np.zeros([noa, 3]) - variances = np.zeros([noa, 3]) - - for atom in range(noa): - env_curr = env.AtomicEnvironment(structure, atom, cutoffs) - - for n in range(3): - d = n+1 - comp_pred, comp_var = gp_model.predict(env_curr, d) - predictions[atom, n] = comp_pred - variances[atom, n] = comp_var - - # mean absolute error - MAE = np.mean(np.abs(forces-predictions)) - - # mean absolute std - stds = np.sqrt(variances) - MAS = np.mean(stds) - - return predictions, variances, forces, MAE, MAS - - -def predict_forces_on_snaps(gp_model, filename, snaps, cutoffs): - preds = np.array([]) - frcs = np.array([]) - for snap in snaps: - struc_curr, forces = get_structure_from_dataset(filename, snap) - predictions, variances, forces, MAE, MAS = \ - predict_forces_on_structure(gp_model, struc_curr, forces, cutoffs) - preds = np.append(preds, predictions) - frcs = np.append(frcs, forces) - return preds, frcs - -# load copper training set (1/9) -box_test = np.load('/Users/jonpvandermause/Downloads/Cu/set.000/box.npy') -coord_test = np.load('/Users/jonpvandermause/Downloads/Cu/set.000/coord.npy') -force_test = np.load('/Users/jonpvandermause/Downloads/Cu/set.000/force.npy') -energy_test = np.load('/Users/jonpvandermause/Downloads/Cu/set.000/energy.npy') - -# initialize gp -kernel = kernels.three_body -kernel_grad = kernels.three_body_grad -hyps = np.array([1, 1, 1]) -cutoffs = np.array([4.2, 4.2]) -opt_algorithm = 'BFGS' -maxiter = 50 - -gp_model = gp.GaussianProcess(kernel, kernel_grad, hyps, cutoffs, - opt_algorithm=opt_algorithm, maxiter=maxiter) - -snaps = [100, 125, 150, 175, 200, 225, 250, 275, 300] -atoms = [1, 25, 50, 75, 100, 125, 150, 175, 200, 225, 250] - -for snap in snaps: - positions = coord_test[snap].reshape(-1, 3) - forces = force_test[snap].reshape(-1, 3) - cell = box_test[snap].reshape(3, 3) - species = ['Cu'] * positions.shape[0] - - struc_curr = struc.Structure(cell, species, positions) - - gp_model.update_db(struc_curr, forces, atoms) - -gp_model.hyps = np.array([0.00658963, 1.16094213, 0.10436963]) -gp_model.set_L_alpha() - -snaps = [100, 200, 300] -file_name = '/Users/jonpvandermause/Downloads/Cu/set.009/' - -predictions, forces = \ - predict_forces_on_snaps(gp_model, file_name, snaps, cutoffs) - -print(np.sqrt(np.mean((predictions - forces)**2))) diff --git a/flare/modules/hyp_opt.py b/flare/modules/hyp_opt.py deleted file mode 100644 index ec1579475..000000000 --- a/flare/modules/hyp_opt.py +++ /dev/null @@ -1,267 +0,0 @@ -import math -from math import exp -import sys -import numpy as np -from scipy.linalg import solve_triangular -import time -sys.path.append('otf/otf_engine') -import env, gp, struc - - -def two_body_grad(env1, env2, d1, d2, hyps): - return two_body_grad_from_env(env1.bond_array, - env1.bond_types, - env2.bond_array, - env2.bond_types, - d1, d2, hyps) - - -def two_body_grad_from_env(bond_array_1, bond_types_1, bond_array_2, - bond_types_2, d1, d2, hyps): - sig = hyps[0] - ls = hyps[1] - S = sig * sig - L = 1 / (ls * ls) - sig_conv = 2 * sig - ls_conv = -2 / (ls * ls * ls) - - kern = 0 - sig_derv = 0 - ls_derv = 0 - - x1_len = len(bond_types_1) - x2_len = len(bond_types_2) - - for m in range(x1_len): - r1 = bond_array_1[m, 0] - coord1 = bond_array_1[m, d1] - typ1 = bond_types_1[m] - - for n in range(x2_len): - r2 = bond_array_2[n, 0] - coord2 = bond_array_2[n, d2] - typ2 = bond_types_2[n] - - # check that bonds match - if typ1 == typ2: - rr = (r1-r2)*(r1-r2) - kern += S*L*exp(-0.5*L*rr)*coord1*coord2*(1-L*rr) - sig_derv += L*exp(-0.5*L*rr)*coord1*coord2*(1-L*rr) * sig_conv - ls_derv += 0.5*coord1*coord2*S*exp(-L*rr/2) * \ - (2+L*rr*(-5+L*rr))*ls_conv - - kern_grad = np.array([sig_derv, ls_derv]) - - return kern, kern_grad - - -def get_likelihood_and_gradients(training_data, training_labels_np, - kernel_grad, hyps, sigma_n): - - number_of_hyps = len(hyps) - - # initialize matrices - size = len(training_data)*3 - k_mat = np.zeros([size, size]) - - # add a matrix to include noise variance: - hyp_mat = np.zeros([size, size, number_of_hyps+1]) - - ds = [1, 2, 3] - - # calculate elements - for m_index in range(size): - x_1 = training_data[int(math.floor(m_index / 3))] - d_1 = ds[m_index % 3] - - for n_index in range(m_index, size): - x_2 = training_data[int(math.floor(n_index / 3))] - d_2 = ds[n_index % 3] - - # calculate kernel and gradient - cov = kernel_grad(x_1, x_2, d_1, d_2, hyps) - - # store kernel value - k_mat[m_index, n_index] = cov[0] - k_mat[n_index, m_index] = cov[0] - - # store gradients (excluding noise variance) - for p_index in range(number_of_hyps): - hyp_mat[m_index, n_index, p_index] = cov[1][p_index] - hyp_mat[n_index, m_index, p_index] = cov[1][p_index] - - # add gradient of noise variance - hyp_mat[:, :, number_of_hyps] = np.eye(size) * 2 * sigma_n - - # matrix manipulation - ky_mat = k_mat + sigma_n ** 2 * np.eye(size) - ky_mat_inv = np.linalg.inv(ky_mat) - alpha = np.matmul(ky_mat_inv, training_labels_np) - alpha_mat = np.matmul(alpha.reshape(alpha.shape[0], 1), - alpha.reshape(1, alpha.shape[0])) - like_mat = alpha_mat - ky_mat_inv - - # calculate likelihood - like = (-0.5*np.matmul(training_labels_np, alpha) - - 0.5*math.log(np.linalg.det(ky_mat)) - - math.log(2 * np.pi) * k_mat.shape[1] / 2) - - # calculate likelihood gradient - like_grad = np.zeros(number_of_hyps + 1) - for n in range(number_of_hyps + 1): - like_grad[n] = 0.5 * np.trace(np.matmul(like_mat, hyp_mat[:, :, n])) - - return like, like_grad - - -def likelihood_ascent(training_data, training_labels_np, kernel_grad, - old_hyps, old_sigma_n, step_factor, max_steps, tol, - verbose=False): - number_of_hyps = len(old_hyps) # excluding noise - new_hyps = old_hyps - new_sigma_n = old_sigma_n - old_like = 1e16 - - for _ in range(max_steps): - like, like_grad = get_likelihood_and_gradients(training_data, - training_labels_np, - kernel_grad, new_hyps, - new_sigma_n) - new_hyps = new_hyps + like_grad[0:number_of_hyps] * step_factor - new_sigma_n = new_sigma_n + like_grad[number_of_hyps] * step_factor - - # check tolerance - if (like > old_like) and (like - old_like < tol): - print(like-old_like) - break - - old_like = like - - if verbose: - print(like) - print(_) - - return new_hyps, new_sigma_n, like - - -# testing ground -if __name__ == '__main__': - # make env1 - cell = np.eye(3) - cutoff = np.linalg.norm(np.array([0.5, 0.5, 0.5])) + 0.001 - - positions_1 = [np.array([0, 0, 0]), np.array([0.1, 0.2, 0.3])] - species_1 = ['B', 'A'] - atom_1 = 0 - test_structure_1 = struc.Structure(cell, species_1, positions_1, cutoff) - env1 = env.ChemicalEnvironment(test_structure_1, atom_1) - - # make env 2 - positions_2 = [np.array([0, 0, 0]), np.array([0.25, 0.3, 0.4])] - species_2 = ['B', 'A'] - atom_2 = 0 - test_structure_2 = struc.Structure(cell, species_2, positions_2, cutoff) - env2 = env.ChemicalEnvironment(test_structure_2, atom_2) - - # set hyperparameters - d1 = 1 - d2 = 1 - sig = 3.1415 - ls = 1 - hyp = np.array([sig, ls]) - - # check if kernel matches old kernel - assert(np.isclose(env.two_body(env1, env2, d1, d2, sig, ls), - two_body_grad(env1, env2, d1, d2, hyp)[0])) - - # check if two body derivative formulas are correct - kern, kern_grad = two_body_grad(env1, env2, d1, d2, hyp) - S_derv = kern_grad[0] - L_derv = kern_grad[1] - - delta = 1e-8 - tol = 1e-5 - new_sig = sig + delta - new_ls = ls + delta - - S_derv_brute = (env.two_body(env1, env2, d1, d2, new_sig, ls) - - env.two_body(env1, env2, d1, d2, sig, ls)) / delta - - L_derv_brute = (env.two_body(env1, env2, d1, d2, sig, new_ls) - - env.two_body(env1, env2, d1, d2, sig, ls)) / delta - - assert(np.isclose(S_derv, S_derv_brute, tol)) - assert(np.isclose(L_derv, L_derv_brute, tol)) - - # check that new likelihood matches old likelihood - sig_n = 5 - forces = [np.array([1, 2, 3]), np.array([4, 5, 6])] - - gp_test = gp.GaussianProcess('two_body') - hyp_gp = [sig, ls, sig_n] - gp_test.update_db(test_structure_1, forces) - gp_like = gp_test.like_hyp(hyp_gp) - - tb = gp_test.training_data - training_labels_np = gp_test.training_labels_np - like, like_grad = get_likelihood_and_gradients(tb, training_labels_np, - two_body_grad, hyp, - sig_n) - new_like = like - - assert(gp_like == new_like) - - # check that likelihood gradient is correct - delta = 1e-7 - new_sig = np.array([sig + delta, ls]) - new_ls = np.array([sig, ls + delta]) - new_n = sig_n + delta - - sig_grad_brute = (get_likelihood_and_gradients(tb, training_labels_np, - two_body_grad, new_sig, - sig_n)[0] - - get_likelihood_and_gradients(tb, training_labels_np, - two_body_grad, hyp, - sig_n)[0])\ - / delta - - ls_grad_brute = (get_likelihood_and_gradients(tb, training_labels_np, - two_body_grad, new_ls, - sig_n)[0] - - get_likelihood_and_gradients(tb, training_labels_np, - two_body_grad, hyp, - sig_n)[0])\ - / delta - - n_grad_brute = (get_likelihood_and_gradients(tb, training_labels_np, - two_body_grad, hyp, - new_n)[0] - - get_likelihood_and_gradients(tb, training_labels_np, - two_body_grad, hyp, - sig_n)[0])\ - / delta - - tol = 1e-3 - assert(np.isclose(like_grad[0], sig_grad_brute, tol)) - assert(np.isclose(like_grad[1], ls_grad_brute)) - assert(np.isclose(like_grad[2], n_grad_brute)) - - # test likelihood ascent - hyp = np.array([10, 10]) - sig_n = 3.89 - step_factor = 0.1 - max_steps = 1000 - tol = 1e-8 - new_hyps, new_sigma_n, like = likelihood_ascent(tb, training_labels_np, - two_body_grad, - hyp, sig_n, step_factor, - max_steps, tol, True) - - gp_test.train() - print(new_hyps) - print(new_sigma_n) - print(like) - print(gp_test.sigma_f) - print(gp_test.length_scale) - print(gp_test.sigma_n) - print(gp_test.get_likelihood()) diff --git a/flare/modules/initialize_velocities.py b/flare/modules/initialize_velocities.py deleted file mode 100644 index 3c22052dc..000000000 --- a/flare/modules/initialize_velocities.py +++ /dev/null @@ -1,40 +0,0 @@ -import numpy as np - - -def get_random_velocities(noa: int, temperature: float, mass: float): - """Draw velocities from Maxwell-Boltzmann distribution. Assumes mass - is given in amu.""" - - mass_md = mass * 0.000103642695727 - kb = 0.0000861733034 - std = np.sqrt(kb * temperature / mass_md) - velocities = np.random.normal(scale=std, size=(noa, 3)) - - vel_sum = np.sum(velocities, axis=0) - corrected_velocities = velocities - vel_sum / noa - - return corrected_velocities - - -# TODO: implement velocity function for multicomponent systems -def multicomponent_velocities(temperature, masses): - noa = len(masses) - velocities = np.zeros((noa, 3)) - kb = 0.0000861733034 - mom_tot = np.array([0., 0., 0.]) - - for n, mass in enumerate(masses): - mass_md = mass * 0.000103642695727 - std = np.sqrt(kb * temperature / mass_md) - rand_vel = np.random.normal(scale=std, size=(3)) - velocities[n] = rand_vel - mom_curr = rand_vel * mass_md - mom_tot += mom_curr - - # correct momentum - mom_corr = mom_tot / noa - for n, mass in enumerate(masses): - mass_md = mass * 0.000103642695727 - velocities[n] -= mom_corr / mass_md - - return velocities diff --git a/flare/modules/lammps/__init__.py b/flare/modules/lammps/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/flare/modules/lammps/coulomb_util.py b/flare/modules/lammps/coulomb_util.py deleted file mode 100644 index d04f2f66b..000000000 --- a/flare/modules/lammps/coulomb_util.py +++ /dev/null @@ -1,132 +0,0 @@ -import numpy as np -import subprocess -import os - - -class EwaldCalculator: - def __init__(self, struc, charges, cutoff, accuracy, lammps_folder, - lammps_executable): - self.struc = struc - self.charges = charges - self.cutoff = cutoff - self.accuracy = accuracy - self.lammps_folder = lammps_folder - self.lammps_executable = lammps_executable - - self.input_file = lammps_folder + '/tmp.in' - self.output_file = lammps_folder + '/tmp.out' - self.dat_file = lammps_folder + '/tmp.data' - self.dump_file = lammps_folder + '/tmp.dump' - - self.input_text = self.lammps_ewald_input() - self.dat_text = self.lammps_ewald_dat() - - def get_ewald_forces(self): - self.lammps_ewald_generator() - self.run_ewald() - forces = self.lammps_ewald_parser() - return forces - - def run_ewald(self): - # create input and data files - self.lammps_ewald_generator() - - # run lammps - lammps_command = '%s < %s > %s' % (self.lammps_executable, - self.input_file, - self.output_file) - os.system(lammps_command) - - def lammps_ewald_generator(self): - self.write_file(self.input_file, self.input_text) - self.write_file(self.dat_file, self.dat_text) - - def lammps_ewald_parser(self): - forces = [] - - with open(self.dump_file, 'r') as outf: - lines = outf.readlines() - - for count, line in enumerate(lines): - if line.startswith('ITEM: ATOMS'): - force_start = count - - for line in lines[force_start+1:]: - fline = line.split() - forces.append([float(fline[-3]), - float(fline[-2]), - float(fline[-1])]) - - return np.array(forces) - - def lammps_ewald_input(self): - input_text = """# long range coulomb calculation -units metal -atom_style charge -dimension 3 -boundary p p p -read_data %s - -pair_style coul/long %f -pair_coeff * * -kspace_style ewald %e - -thermo_style one -dump 1 all custom 1 %s id type x y z q fx fy fz -dump_modify 1 sort id -run 0 -""" % (self.dat_file, self.cutoff, self.accuracy, self.dump_file) - - return input_text - - def lammps_ewald_dat(self): - dat_text = """Header of the LAMMPS data file - -%i atoms -1 atom types -""" % (self.struc.nat) - - dat_text += self.lammps_cell_text() - - dat_text += """ -Masses - -1 1 - -Atoms -""" - - dat_text += self.lammps_pos_text() - - return dat_text - - def lammps_cell_text(self): - """ Write cell from structure object. Assumes orthorombic periodic - cell.""" - - cell_text = """ -0.0 %f xlo xhi -0.0 %f ylo yhi -0.0 %f zlo zhi -%f %f %f xy xz yz -""" % (self.struc.cell[0, 0], - self.struc.cell[1, 1], - self.struc.cell[2, 2], - self.struc.cell[1, 0], - self.struc.cell[2, 0], - self.struc.cell[2, 1]) - - return cell_text - - def lammps_pos_text(self): - pos_text = '\n' - for count, pos in enumerate(self.struc.positions): - pos_text += '%i 1 %f %f %f %f \n' % \ - (count+1, self.charges[self.struc.coded_species[count]], - pos[0], pos[1], pos[2]) - return pos_text - - @staticmethod - def write_file(file, text): - with open(file, 'w') as fin: - fin.write(text) diff --git a/flare/modules/lammps/eam.py b/flare/modules/lammps/eam.py deleted file mode 100644 index b0b70275a..000000000 --- a/flare/modules/lammps/eam.py +++ /dev/null @@ -1,129 +0,0 @@ -import numpy as np -import subprocess -import os - - -class EAM_Force_Calculator: - def __init__(self, struc, style_string, coeff_string, lammps_folder, - lammps_executable): - self.struc = struc - self.style_string = style_string - self.coeff_string = coeff_string - self.lammps_folder = lammps_folder - self.lammps_executable = lammps_executable - - self.input_file = lammps_folder + '/tmp.in' - self.output_file = lammps_folder + '/tmp.out' - self.dat_file = lammps_folder + '/tmp.data' - self.dump_file = lammps_folder + '/tmp.dump' - - self.input_text = self.lammps_input() - self.dat_text = self.lammps_dat() - - def get_forces(self): - self.lammps_generator() - self.run_ewald() - forces = self.lammps_parser() - return forces - - def run_ewald(self): - # create input and data files - self.lammps_generator() - - # run lammps - lammps_command = '%s < %s > %s' % (self.lammps_executable, - self.input_file, - self.output_file) - os.system(lammps_command) - - def lammps_generator(self): - self.write_file(self.input_file, self.input_text) - self.write_file(self.dat_file, self.dat_text) - - def lammps_parser(self): - forces = [] - - with open(self.dump_file, 'r') as outf: - lines = outf.readlines() - - for count, line in enumerate(lines): - if line.startswith('ITEM: ATOMS'): - force_start = count - - for line in lines[force_start+1:]: - fline = line.split() - forces.append([float(fline[-3]), - float(fline[-2]), - float(fline[-1])]) - - return np.array(forces) - - def lammps_input(self): - input_text = """# lammps input file created with eam.py. -units metal -atom_style atomic -dimension 3 -boundary p p p -read_data %s - -pair_style %s -pair_coeff %s - -thermo_style one -dump 1 all custom 1 %s id type x y z fx fy fz -dump_modify 1 sort id -run 0 -""" % (self.dat_file, self.style_string, self.coeff_string, self.dump_file) - - return input_text - - def lammps_dat(self): - dat_text = """Header of the LAMMPS data file - -%i atoms -1 atom types -""" % (self.struc.nat) - - dat_text += self.lammps_cell_text() - - dat_text += """ -Masses - -1 1 - -Atoms -""" - - dat_text += self.lammps_pos_text() - - return dat_text - - def lammps_cell_text(self): - """ Write cell from structure object. Assumes orthorombic periodic - cell.""" - - cell_text = """ -0.0 %f xlo xhi -0.0 %f ylo yhi -0.0 %f zlo zhi -%f %f %f xy xz yz -""" % (self.struc.cell[0, 0], - self.struc.cell[1, 1], - self.struc.cell[2, 2], - self.struc.cell[1, 0], - self.struc.cell[2, 0], - self.struc.cell[2, 1]) - - return cell_text - - def lammps_pos_text(self): - pos_text = '\n' - for count, pos in enumerate(self.struc.positions): - pos_text += '%i 1 %f %f %f \n' % \ - (count+1, pos[0], pos[1], pos[2]) - return pos_text - - @staticmethod - def write_file(file, text): - with open(file, 'w') as fin: - fin.write(text) diff --git a/flare/modules/lammps/mc.py b/flare/modules/lammps/mc.py deleted file mode 100644 index bc951e297..000000000 --- a/flare/modules/lammps/mc.py +++ /dev/null @@ -1,171 +0,0 @@ -import numpy as np -import subprocess -import os - - -class EAM_Force_Calculator: - def __init__(self, struc, style_string, coeff_string, lammps_folder, - lammps_executable, atom_types, atom_masses, species, - charges=None): - self.struc = struc - self.style_string = style_string - self.coeff_string = coeff_string - self.lammps_folder = lammps_folder - self.lammps_executable = lammps_executable - self.atom_types = atom_types - self.atom_masses = atom_masses - self.species = species - self.charges = charges - - self.input_file = lammps_folder + '/tmp.in' - self.output_file = lammps_folder + '/tmp.out' - self.dat_file = lammps_folder + '/tmp.data' - self.dump_file = lammps_folder + '/tmp.dump' - - self.input_text = self.lammps_input() - self.dat_text = self.lammps_dat() - - def get_forces(self): - self.lammps_generator() - self.run_ewald() - forces = self.lammps_parser() - return forces - - def run_ewald(self): - # create input and data files - self.lammps_generator() - - # run lammps - lammps_command = '%s < %s > %s' % (self.lammps_executable, - self.input_file, - self.output_file) - os.system(lammps_command) - - def lammps_generator(self): - self.write_file(self.input_file, self.input_text) - self.write_file(self.dat_file, self.dat_text) - - def lammps_parser(self): - forces = [] - - with open(self.dump_file, 'r') as outf: - lines = outf.readlines() - - for count, line in enumerate(lines): - if line.startswith('ITEM: ATOMS'): - force_start = count - - for line in lines[force_start+1:]: - fline = line.split() - forces.append([float(fline[-3]), - float(fline[-2]), - float(fline[-1])]) - - return np.array(forces) - - def lammps_input(self): - input_text = """# lammps input file created with eam.py. -units metal -atom_style atomic -dimension 3 -boundary p p p -newton off -read_data %s - -pair_style %s -pair_coeff %s - -thermo_style one -dump 1 all custom 1 %s id type x y z fx fy fz -dump_modify 1 sort id -run 0 -""" % (self.dat_file, self.style_string, self.coeff_string, self.dump_file) - - return input_text - - def lammps_dat(self): - dat_text = """Header of the LAMMPS data file - -%i atoms -%i atom types -""" % (self.struc.nat, len(self.atom_types)) - - dat_text += self.lammps_cell_text() - dat_text += """ -Masses - -""" - mass_text = '' - for atom_type, atom_mass in zip(self.atom_types, self.atom_masses): - mass_text += '%i %i\n' % (atom_type, atom_mass) - dat_text += mass_text - dat_text += """ -Atoms -""" - dat_text += self.lammps_pos_text() - - return dat_text - - def lammps_dat_charged(self): - dat_text = """Header of the LAMMPS data file - -%i atoms -%i atom types -""" % (self.struc.nat, len(self.atom_types)) - - dat_text += self.lammps_cell_text() - dat_text += """ -Masses - -""" - mass_text = '' - for atom_type, atom_mass in zip(self.atom_types, - self.atom_masses): - mass_text += '%i %i\n' % (atom_type, atom_mass) - dat_text += mass_text - dat_text += """ -Atoms -""" - dat_text += self.lammps_pos_text_charged() - - return dat_text - - def lammps_cell_text(self): - """ Write cell from structure object. Assumes orthorombic periodic - cell.""" - - cell_text = """ -0.0 %f xlo xhi -0.0 %f ylo yhi -0.0 %f zlo zhi -%f %f %f xy xz yz -""" % (self.struc.cell[0, 0], - self.struc.cell[1, 1], - self.struc.cell[2, 2], - self.struc.cell[1, 0], - self.struc.cell[2, 0], - self.struc.cell[2, 1]) - - return cell_text - - def lammps_pos_text(self): - pos_text = '\n' - for count, (pos, spec) in enumerate(zip(self.struc.positions, - self.species)): - pos_text += '%i %i %f %f %f \n' % \ - (count+1, spec, pos[0], pos[1], pos[2]) - return pos_text - - def lammps_pos_text_charged(self): - pos_text = '\n' - for count, (pos, chrg, spec) in enumerate(zip(self.struc.positions, - self.charges, - self.species)): - pos_text += '%i %i %f %f %f %f \n' % \ - (count+1, spec, chrg, pos[0], pos[1], pos[2]) - return pos_text - - @staticmethod - def write_file(file, text): - with open(file, 'w') as fin: - fin.write(text) diff --git a/flare/modules/qe_input.py b/flare/modules/qe_input.py deleted file mode 100644 index de48fc908..000000000 --- a/flare/modules/qe_input.py +++ /dev/null @@ -1,277 +0,0 @@ -import numpy as np -import subprocess -import os - - -class BashInput: - def __init__(self, bash_name, bash_inputs): - self.bash_name = bash_name - - self.n = bash_inputs['n'] - self.N = bash_inputs['N'] - self.t = bash_inputs['t'] - self.e = bash_inputs['e'] - self.p = bash_inputs['p'] - self.o = bash_inputs['o'] - self.mem_per_cpu = bash_inputs['mem_per_cpu'] - self.mail_user = bash_inputs['mail_user'] - self.command = bash_inputs['command'] - - self.bash_text = self.get_bash_text() - - def get_bash_text(self): - sh_text = """#!/bin/sh -#SBATCH -n {} -#SBATCH -N {} -#SBATCH -t {}-00:00 -#SBATCH -e {} -#SBATCH -p {} -#SBATCH -o {} -#SBATCH --mem-per-cpu={} -#SBATCH --mail-type=ALL -#SBATCH --mail-user={} - -module load gcc/4.9.3-fasrc01 openmpi/2.1.0-fasrc01 -module load python/3.6.3-fasrc01 - -{}""".format(self.n, self.N, self.t, self.e, self.p, self.o, - self.mem_per_cpu, self.mail_user, self.command) - - return sh_text - - def write_bash_text(self): - with open(self.bash_name, 'w') as fin: - fin.write(self.bash_text) - - -class QEInput: - def __init__(self, input_file_name: str, output_file_name: str, - pw_loc: str, calculation: str, - scf_inputs: dict, md_inputs: dict = None, - press_conv_thr='0.5', electron_maxstep=100, - metal=False, mixing_beta=0.7, smearing='mp', - degauss=0.02): - - self.input_file_name = input_file_name - self.output_file_name = output_file_name - self.pw_loc = pw_loc - self.calculation = calculation - - self.pseudo_dir = scf_inputs['pseudo_dir'] - self.outdir = scf_inputs['outdir'] - self.nat = scf_inputs['nat'] - self.ntyp = scf_inputs['ntyp'] - self.ecutwfc = scf_inputs['ecutwfc'] - self.ecutrho = scf_inputs['ecutrho'] - self.cell = scf_inputs['cell'] - self.species = scf_inputs['species'] - self.positions = scf_inputs['positions'] - self.kvec = scf_inputs['kvec'] - self.ion_names = scf_inputs['ion_names'] - self.ion_masses = scf_inputs['ion_masses'] - self.ion_pseudo = scf_inputs['ion_pseudo'] - self.electron_maxstep = electron_maxstep - self.metal = metal - self.mixing_beta = mixing_beta - self.smearing = smearing - self.degauss = degauss - - # get text blocks - self.species_txt = self.get_species_txt() - self.position_txt = self.get_position_txt() - self.cell_txt = self.get_cell_txt() - self.kpt_txt = self.get_kpt_txt() - - # get md parameters - if self.calculation == 'md': - self.dt = md_inputs['dt'] - self.nstep = md_inputs['nstep'] - self.ion_temperature = md_inputs['ion_temperature'] - self.tempw = md_inputs['tempw'] - - # if vc, get pressure convergence threshold - if self.calculation == 'vc-relax': - self.press_conv_thr = press_conv_thr - - self.input_text = self.get_input_text() - - # write input file - self.write_file() - - def get_species_txt(self): - spectxt = '' - spectxt += 'ATOMIC_SPECIES' - for name, mass, pseudo in zip(self.ion_names, - self.ion_masses, - self.ion_pseudo): - spectxt += '\n {} {} {}'.format(name, mass, pseudo) - - return spectxt - - def get_position_txt(self): - postxt = '' - postxt += 'ATOMIC_POSITIONS {angstrom}' - for spec, pos in zip(self.species, self.positions): - postxt += '\n {} {:1.8f} {:1.8f} {:1.8f}'.format(spec, *pos) - - return postxt - - def get_cell_txt(self): - celltxt = '' - celltxt += 'CELL_PARAMETERS {angstrom}' - for vector in self.cell: - celltxt += '\n {:1.5f} {:1.5f} {:1.5f}'.format(*vector) - - return celltxt - - def get_kpt_txt(self): - ktxt = '' - ktxt += 'K_POINTS automatic' - ktxt += '\n {} {} {} 0 0 0'.format(*self.kvec) - - return ktxt - - def get_input_text(self): - - input_text = """ &control - calculation = '{}' - pseudo_dir = '{}' - outdir = '{}' - tstress = .true. - tprnfor = .true.""".format(self.calculation, self.pseudo_dir, self.outdir) - - # if MD, add time step and number of steps - if self.calculation == 'md': - input_text += """ - dt = {} - nstep = {}""".format(self.dt, self.nstep) - - input_text += """ - / - &system - ibrav= 0 - nat= {} - ntyp= {} - ecutwfc ={} - ecutrho = {}""".format(self.nat, self.ntyp, self.ecutwfc, - self.ecutrho) - - # if metal is true, add smearing - if self.metal is True: - input_text += """ - occupations = 'smearing' - smearing = '{}' - degauss = {}""".format(self.smearing, self.degauss) - - # if MD or relax, don't reduce number of k points based on symmetry, - # since the symmetry might change throughout the calculation - if (self.calculation == 'md') or (self.calculation == 'relax'): - input_text += """ - nosym = .true.""" - - input_text += """ - / - &electrons - conv_thr = 1.0d-6 - mixing_beta = {} - electron_maxstep = {}""".format(self.mixing_beta, self.electron_maxstep) - - # if MD or relax, need to add an &IONS block - if (self.calculation == 'md') or (self.calculation == 'relax') or \ - (self.calculation == 'vc-relax'): - input_text += """ - / - &ions""" - - # if MD, add additional details about ions - if self.calculation == 'md': - input_text += """ - pot_extrapolation = 'second-order' - wfc_extrapolation = 'second-order' - ion_temperature = '{}' - tempw = {}""".format(self.ion_temperature, self.tempw) - - # if vc-relax, add cell block - if self.calculation == 'vc-relax': - input_text += """ - / - &cell - press_conv_thr = {}""".format(self.press_conv_thr) - - # insert species, cell, position and k-point textblocks - input_text += """ - / -{} -{} -{} -{} -""".format(self.species_txt, self.cell_txt, self.position_txt, self.kpt_txt) - - return input_text - - def write_file(self): - with open(self.input_file_name, 'w') as fin: - fin.write(self.input_text) - - def run_espresso(self, npool): - if npool is False: - qe_command = \ - 'mpirun {0} < {1} > {2}'.format( - self.pw_loc, self.input_file_name, self.output_file_name) - else: - qe_command = \ - 'mpirun -npool {0} {1} < {2} > {3}'.format( - npool, self.pw_loc, self.input_file_name, - self.output_file_name) - os.system(qe_command) - # subprocess.call(qe_command, shell=True) - - -def create_structure(el: str, alat: float, size: int, perturb: bool=False, - pass_relax: bool=False, - pass_pos: bool= None): - """ - Create bulk structure with vacancy, return cell, position, number of atoms - - :param el: - :param alat: - :param size: - :param perturb: - :param pass_relax: - :param pass_pos: - :return: - """ - - # create bulk cell - unit_cell = crystal(el, [(0, 0, 0)], spacegroup=225, cellpar=[alat, alat, alat, 90, 90, 90]) - - # size of initial perturbation - pert_size = 0.1 * alat - print("Pert Size", pert_size) - - # make supercell - multiplier = np.identity(3) * size - supercell = make_supercell(unit_cell, multiplier) - - # remove atom - supercell.pop(supercell.get_number_of_atoms() // 2) - - # get unpertubed positions - al_pos = np.asarray(supercell.positions) - - if pass_relax: - al_pos = pass_pos - - if perturb: - for atom in range(al_pos.shape[0]): - for coord in range(3): - al_pos[atom][coord] += np.random.uniform(-pert_size, pert_size) - - cell = supercell.get_cell() - nat = supercell.get_number_of_atoms() - - return cell, al_pos, nat - - -if __name__ == '__main__': - pass diff --git a/flare/modules/qe_parsers.py b/flare/modules/qe_parsers.py deleted file mode 100644 index 086ee93c0..000000000 --- a/flare/modules/qe_parsers.py +++ /dev/null @@ -1,147 +0,0 @@ -import numpy as np - -# unit conversion from http://greif.geo.berkeley.edu/~driver/conversions.html -en_conv = 13.6056925330 # bohr to eV -force_conv = 25.71104309541616 # bohr/Ry to eV/A -stress_conv = 25.71104309541616 / (0.529177208**2) - - -# get energy from scf run -def parse_scf_energy(outfile): - with open(outfile, 'r') as outf: - lines = outf.readlines() - - for line in lines: - if '!' in line: - energy = float(line.split()[-2]) - - return energy - - -def parse_scf(outfile: str): - """ - Get forces from a pwscf file in eV/A - - :param outfile: str, Path to pwscf output file - :return: list[nparray] , List of forces acting on atoms - """ - forces = [] - total_energy = np.nan - - with open(outfile, 'r') as outf: - for line in outf: - if line.lower().startswith('! total energy'): - total_energy = float(line.split()[-2]) - - if line.find('force') != -1 and line.find('atom') != -1: - line = line.split('force =')[-1] - line = line.strip() - line = line.split(' ') - line = [x for x in line if x != ''] - temp_forces = [] - for x in line: - temp_forces.append(float(x)) - forces.append(np.array(list(temp_forces))) - - assert total_energy != np.nan, "Quantum ESPRESSO parser failed to read " \ - "the file {}. Run failed.".format(outfile) - - total_energy = total_energy * en_conv - forces = [force_conv * force for force in forces] - - return total_energy, forces - - -def parse_md_output(outfile): - - steps = {} - - with open(outfile, 'r') as outf: - lines = outf.readlines() - - split_indexes = [N for N in range(len(lines)) if '!' == lines[N][0]] - - step_chunks = [] - for n in range(len(split_indexes)): - step_chunks.append(lines[split_indexes[n]:split_indexes[n+1] if - n != len(split_indexes)-1 else len(lines)]) - - for current_chunk in step_chunks: - - force_start_line = [line for line in current_chunk if - 'Forces acting on atoms' in line][0] - force_end_line = [line for line in current_chunk if - 'Total force' in line][0] - force_start_index = current_chunk.index(force_start_line)+2 - force_end_index = current_chunk.index(force_end_line)-2 - - atoms_start_line = [line for line in current_chunk if - 'ATOMIC_POSITIONS' in line][0] - atoms_end_line = [line for line in current_chunk if - 'kinetic energy' in line][0] - atoms_start_index = current_chunk.index(atoms_start_line)+1 - atoms_end_index = current_chunk.index(atoms_end_line)-3 - - stress_start_line = [line for line in current_chunk if - 'Computing stress' in line][0] - stress_start_index = current_chunk.index(stress_start_line) + 3 - - temperature_line = [line for line in current_chunk if - 'temperature' in line][0] - dyn_line = [line for line in current_chunk if - 'Entering Dynamics' in line][0] - dyn_index = current_chunk.index(dyn_line) - time_index = dyn_index+1 - - forces = [] - for line in current_chunk[force_start_index:force_end_index+1]: - forceline = line.split('=')[-1].split() - forces.append(force_conv * np.array([float(forceline[0]), - float(forceline[1]), - float(forceline[2])])) - forces = np.array(forces) - total_force = float(force_end_line.split('=')[1].strip().split()[0]) - SCF_corr = float(force_end_line.split('=')[2].strip()[0]) - - stress = [] - for line in current_chunk[stress_start_index:stress_start_index+3]: - stress_line = line.split() - stress.append(np.array([float(stress_line[0]), - float(stress_line[1]), - float(stress_line[2])])) - stress = np.array(stress) - stress = stress * stress_conv # convert to eV/A^3 - - positions = [] - elements = [] - for line in current_chunk[atoms_start_index:atoms_end_index+1]: - atomline = line.split() - elements.append(atomline[0]) - positions.append(np.array([float(atomline[1]), float(atomline[2]), - float(atomline[3])])) - positions = np.array(positions) - - toten = \ - float(current_chunk[0].split('=')[-1].strip().split()[0]) * en_conv - temperature_line = temperature_line.split('=')[-1] - temperature = float(temperature_line.split()[0]) - iteration = int(dyn_line.split('=')[-1]) - timeline = current_chunk[time_index].split('=')[-1].strip().split()[0] - time = float(timeline) - Ekin = float(atoms_end_line.split('=')[1].strip().split()[0]) - - steps[iteration] = {'iteration': iteration, - 'forces': forces, - 'stress': stress, - 'positions': positions, - 'elements': elements, - 'temperature': temperature, - 'time': time, - 'energy': toten, - 'ekin': Ekin, - 'kinetic energy': Ekin, - 'total energy': toten, - 'total force': total_force, - 'SCF correction': SCF_corr} - - return(steps) diff --git a/flare/modules/rcut.py b/flare/modules/rcut.py deleted file mode 100644 index 57df9e104..000000000 --- a/flare/modules/rcut.py +++ /dev/null @@ -1,153 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Self-Contained Code which allows for the -""" - -import sys -import os - -import numpy as np -import numpy.random as rand - -sys.path.append('../otf_engine') - -from math import sin, cos -from struc import Structure -from qe_util import run_espresso, parse_qe_input - - -def perturb_position(pos, r_pert: float, rscale): - """ - - :param pos: - :param r_pert: - :param rscale: - :return: - """ - - theta = rand.uniform(-180, 180) - x = rand.uniform(0, 1) - phi = np.arccos(2 * x - 1) - - r = np.abs(rand.normal(loc=r_pert, scale=rscale)) - - newpos = np.zeros(3) - newpos[0] = pos[0] + r * sin(phi) * cos(theta) - newpos[1] = pos[1] + r * sin(phi) * sin(theta) - newpos[2] = pos[2] + r * cos(phi) - - return newpos - - -def perturb_outside_radius(structure: Structure, atom: int, r_fix: float, - mean_pert: float, pert_sigma: float): - """ - - :param structure: - :param atom: - :param r_fix: - :param mean_pert: - :param pert_sigma: - :return: - """ - assert r_fix > 0, "Radius must be positive" - assert atom >= 0, "Atomic index must be non-negative" - assert mean_pert >= 0, "Perturbation must be non-negative" - - new_positions = [] - central_pos = structure.positions[atom] - - for pos in structure.positions: - if is_within_r_periodic(structure, central_pos, pos, r_fix): - new_positions.append(pos) - else: - new_positions.append(perturb_position(pos, mean_pert, pert_sigma)) - - newstruc = Structure(lattice=structure.lattice, species=structure.species, - positions=new_positions, cutoff=structure.cutoff) - - return newstruc - - -def is_within_r_periodic(structure, central_pos, neighbor_pos, radius): - """ - - :param structure: - :param central_pos: - :param neighbor_pos: - :param radius: - :return: - """ - images = structure.get_periodic_images(neighbor_pos, super_check=2) - # Check to see if images are within the radius - for image in images: - if np.linalg.norm(image - central_pos) < radius: - return True - - return False - - -def gauge_force_variance(qe_input: str, trials: int, atom: int, r_fix: - float, mean_pert: float = .1, pert_sigma: float = .02, write_output: - bool = True): - """ - - :param qe_input: - :param trials: - :param atom: - :param r_fix: - :param mean_pert: - :param pert_sigma: - :param write_output: - :return: - """ - positions, species, cell, masses = parse_qe_input(qe_input) - struc = Structure(positions=positions, species=species, lattice=cell, - cutoff=1, mass_dict=masses) - - total_forces = np.empty(shape=(trials, struc.nat, 3)) - pw_loc = os.environ.get('PWSCF_COMMAND') - - for i in range(trials): - pert_struct = perturb_outside_radius(struc, atom, r_fix, mean_pert, - pert_sigma) - pert_forces = run_espresso(qe_input=qe_input, structure=pert_struct, - pw_loc=pw_loc) - - for j in range(struc.nat): - for k in range(3): - total_forces[i, j, k] = float(pert_forces[j][k]) - - if write_output: - with open("rcut.out", 'a') as f: - f.write('Perturbed structure about atom {} and radius {}' - ': \n'.format(atom, r_fix)) - for pos in pert_struct.positions: - f.write(str(pos)+'\n') - f.write("All Forces (normed force on atom {}: {}) : " - "\n".format(atom, - np.linalg.norm(total_forces[i,atom,:]))) - for force in total_forces[i,:]: - f.write(str(force)+'\n') - #f.write(str(total_forces[i, :, :])) - - all_forces = [total_forces[i][atom][0:3] for i in range(trials)] - - report_str = 'Std dev of force on atom {} : {} Ry/Au \n'.format(atom, - np.std([np.linalg.norm(force) for force in all_forces])) - - print(report_str) - - if write_output: - with open("rcut.out", 'a') as f: - f.write(report_str) - - return total_forces - - -if __name__ == '__main__': - # gauge_force_variance(qe_input='qe_input_1.in', trials=5, atom=0, - # r_fix=3.0, - # mean_pert=.1) - pass diff --git a/flare/modules/test_otf.py b/flare/modules/test_otf.py deleted file mode 100644 index 2a716fcc1..000000000 --- a/flare/modules/test_otf.py +++ /dev/null @@ -1,78 +0,0 @@ -import numpy as np -from flare.modules import otf_parser, analyze_md -from flare import kernels, otf, env -import concurrent.futures -import time - - -class TestOTF: - def __init__(self, otf_file, gp_cell, kernel, - kernel_grad, cutoffs, start_list, test_snaps): - self.test_snaps = test_snaps - self.struc = None # populate later - - # make gp and aimd run - otf_run = otf_parser.OtfAnalysis(otf_file) - call_no = len(otf_run.gp_position_list) - self.gp_model = otf_run.make_gp(gp_cell, kernel, kernel_grad, 'BFGS', - call_no, start_list, cutoffs) - - # adapted from analyze_gp.py - def predict_forces_on_test_set(self, aimd_file, aimd_cell): - # get aimd trajectory - aimd_run = analyze_md.MDAnalysis(aimd_file, aimd_cell) - - all_predictions = np.array([]) - all_stds = np.array([]) - all_forces = np.array([]) - - for snap in self.test_snaps: - # get structure and forces from AIMD trajectory - self.struc = aimd_run.get_structure_from_snap(snap) - forces_curr = aimd_run.get_forces_from_snap(snap) - - # predict forces and stds - self.predict_on_structure_par() - print(time.time()) - predictions = self.struc.forces - stds = self.struc.stds - - # append results - all_predictions = np.append(all_predictions, predictions) - all_stds = np.append(all_stds, stds) - all_forces = np.append(all_forces, forces_curr) - - return all_predictions, all_stds, all_forces - - def predict_on_structure_par(self): - atom_list = list(range(self.struc.positions.shape[0])) - n = 0 - with concurrent.futures.ProcessPoolExecutor() as executor: - for res in executor.map(self.predict_on_atom, atom_list): - for i in range(3): - self.struc.forces[n][i] = res[0][i] - self.struc.stds[n][i] = res[1][i] - n += 1 - - def predict_on_structure(self): - for n in range(self.struc.nat): - chemenv = \ - env.AtomicEnvironment(self.struc, n, self.gp_model.cutoffs) - for i in range(3): - force, var = self.gp_model.predict(chemenv, i + 1) - self.struc.forces[n][i] = float(force) - self.struc.stds[n][i] = np.sqrt(np.absolute(var)) - - # take prediction functions from otf - def predict_on_atom(self, atom): - chemenv = \ - env.AtomicEnvironment(self.struc, atom, self.gp_model.cutoffs) - comps = [] - stds = [] - # predict force components and standard deviations - for i in range(3): - force, var = self.gp_model.predict(chemenv, i+1) - comps.append(float(force)) - stds.append(np.sqrt(np.absolute(var))) - - return comps, stds diff --git a/flare/modules/vacancy_diffusion.py b/flare/modules/vacancy_diffusion.py deleted file mode 100644 index 96d67ab5e..000000000 --- a/flare/modules/vacancy_diffusion.py +++ /dev/null @@ -1,97 +0,0 @@ -import numpy as np -import sys -import crystals -import copy -import eam -sys.path.append('../../otf/otf_engine') -import gp, struc, qe_parsers, env - - -# given a gp model, predict vacancy diffusion activation profile -def vac_diff_fcc(gp_model, gp_cell, fcc_cell, species, cutoff, nop): - - # create 2x2x2 fcc supercell with atom 0 removed - alat = fcc_cell[0, 0] - fcc_unit = crystals.fcc_positions(alat) - fcc_super = crystals.get_supercell_positions(2, fcc_cell, fcc_unit) - vac_super = copy.deepcopy(fcc_super) - vac_super.pop(0) - vac_super = np.array(vac_super) - - # create list of positions for the migrating atom - start_pos = vac_super[0] - end_pos = np.array([0, 0, 0]) - diff_vec = end_pos - start_pos - test_list = [] - step = diff_vec / (nop - 1) - for n in range(nop): - test_list.append(start_pos + n*step) - - # predict force on migrating atom - vac_copy = copy.deepcopy(vac_super) - store_res = np.zeros((6, nop)) - - cutoffs = np.array([cutoff, cutoff]) - - for count, pos_test in enumerate(test_list): - vac_copy[0] = pos_test - - struc_curr = struc.Structure(gp_cell, species, vac_copy) - env_curr = env.AtomicEnvironment(struc_curr, 0, cutoffs) - - fc_comp_x, vr_x = gp_model.predict(env_curr, 1) - fc_comp_y, vr_y = gp_model.predict(env_curr, 2) - fc_comp_z, vr_z = gp_model.predict(env_curr, 3) - - store_res[0, count] = fc_comp_x - store_res[1, count] = fc_comp_y - store_res[2, count] = fc_comp_z - store_res[3, count] = vr_x - store_res[4, count] = vr_y - store_res[5, count] = vr_z - - return store_res - - -# given a gp model, predict vacancy diffusion activation profile -def vac_diff_lammps(style_string, coeff_string, lammps_folder, - lammps_executable, gp_cell, fcc_cell, species, - nop): - - # create 2x2x2 fcc supercell with atom 0 removed - alat = fcc_cell[0, 0] - fcc_unit = crystals.fcc_positions(alat) - fcc_super = crystals.get_supercell_positions(2, fcc_cell, fcc_unit) - vac_super = copy.deepcopy(fcc_super) - vac_super.pop(0) - vac_super = np.array(vac_super) - - # create list of positions for the migrating atom - start_pos = vac_super[0] - end_pos = np.array([0, 0, 0]) - diff_vec = end_pos - start_pos - test_list = [] - step = diff_vec / (nop - 1) - for n in range(nop): - test_list.append(start_pos + n*step) - - # predict force on migrating atom - vac_copy = copy.deepcopy(vac_super) - store_res = np.zeros((3, nop)) - - for count, pos_test in enumerate(test_list): - vac_copy[0] = pos_test - - struc_curr = struc.Structure(gp_cell, species, vac_copy) - - lammps_calculator = \ - eam.EAM_Force_Calculator(struc_curr, style_string, coeff_string, - lammps_folder, lammps_executable) - - lammps_forces = lammps_calculator.get_forces() - - store_res[0, count] = lammps_forces[0, 0] - store_res[1, count] = lammps_forces[0, 1] - store_res[2, count] = lammps_forces[0, 2] - - return store_res diff --git a/flare/modules/validate.py b/flare/modules/validate.py deleted file mode 100644 index 333bedb5f..000000000 --- a/flare/modules/validate.py +++ /dev/null @@ -1,200 +0,0 @@ -import sys -import numpy as np -import datetime -import time -from typing import List -import copy -import multiprocessing as mp -import concurrent.futures -from flare import struc, gp, env, qe_util, md, output - - -class Validate: - def __init__(self, qe_input: str, dt: float, number_of_steps: int, - dft_steps, gp: gp.GaussianProcess, pw_loc: str, - prev_pos_init: np.ndarray=None, par: bool=False, skip: int=0, - calculate_energy=False, output_name='otf_run.out', no_cpus=1): - - self.qe_input = qe_input - self.dt = dt - self.number_of_steps = number_of_steps - self.gp = gp - self.pw_loc = pw_loc - self.skip = skip - self.dft_step = True - self.dft_steps = dft_steps - - # parse input file - positions, species, cell, masses = \ - qe_util.parse_qe_input(self.qe_input) - - self.structure = struc.Structure(cell=cell, species=species, - positions=positions, - mass_dict=masses, - prev_positions=prev_pos_init) - - self.noa = self.structure.positions.shape[0] - self.atom_list = list(range(self.noa)) - self.curr_step = 0 - - # initialize local energies - if calculate_energy: - self.local_energies = np.zeros(self.noa) - else: - self.local_energies = None - - self.dft_count = 0 - - # set pred function - if not par and not calculate_energy: - self.pred_func = self.predict_on_structure - elif par and not calculate_energy: - self.pred_func = self.predict_on_structure_par - elif not par and calculate_energy: - self.pred_func = self.predict_on_structure_en - elif par and calculate_energy: - self.pred_func = self.predict_on_structure_par_en - self.par = par - - self.output_name = output_name - - # set number of cpus for qe runs - self.no_cpus = no_cpus - - def run(self): - output.write_header(self.gp.cutoffs, self.gp.kernel_name, self.gp.hyps, - self.gp.algo, self.dt, self.number_of_steps, - self.structure, self.output_name, 1.) - counter = 0 - self.start_time = time.time() - - while self.curr_step < self.number_of_steps: - self.pred_func() - self.dft_step = False - new_pos = md.update_positions(self.dt, self.noa, - self.structure) - - if self.curr_step in self.dft_steps: - # record GP forces - self.update_temperature(new_pos) - self.record_state() - - # run DFT and record forces - self.dft_step = True - self.run_dft() - new_pos = md.update_positions(self.dt, self.noa, - self.structure) - self.update_temperature(new_pos) - self.record_state() - - # write gp forces - if counter >= self.skip and not self.dft_step: - self.update_temperature(new_pos) - self.record_state() - counter = 0 - - counter += 1 - self.update_positions(new_pos) - self.curr_step += 1 - - output.conclude_run(self.output_name) - - def predict_on_atom(self, atom): - chemenv = env.AtomicEnvironment(self.structure, atom, self.gp.cutoffs) - comps = [] - stds = [] - # predict force components and standard deviations - for i in range(3): - force, var = self.gp.predict(chemenv, i+1) - comps.append(float(force)) - stds.append(np.sqrt(np.abs(var))) - - return comps, stds - - def predict_on_atom_en(self, atom): - chemenv = env.AtomicEnvironment(self.structure, atom, self.gp.cutoffs) - comps = [] - stds = [] - # predict force components and standard deviations - for i in range(3): - force, var = self.gp.predict(chemenv, i+1) - comps.append(float(force)) - stds.append(np.sqrt(np.abs(var))) - - # predict local energy - local_energy = self.gp.predict_local_energy(chemenv) - return comps, stds, local_energy - - def predict_on_structure_par(self): - n = 0 - with concurrent.futures.ProcessPoolExecutor() as executor: - for res in executor.map(self.predict_on_atom, self.atom_list): - for i in range(3): - self.structure.forces[n][i] = res[0][i] - self.structure.stds[n][i] = res[1][i] - n += 1 - - def predict_on_structure_par_en(self): - n = 0 - with concurrent.futures.ProcessPoolExecutor() as executor: - for res in executor.map(self.predict_on_atom_en, self.atom_list): - for i in range(3): - self.structure.forces[n][i] = res[0][i] - self.structure.stds[n][i] = res[1][i] - self.local_energies[n] = res[2] - n += 1 - - def predict_on_structure(self): - for n in range(self.structure.nat): - chemenv = env.AtomicEnvironment(self.structure, n, self.gp.cutoffs) - for i in range(3): - force, var = self.gp.predict(chemenv, i + 1) - self.structure.forces[n][i] = float(force) - self.structure.stds[n][i] = np.sqrt(np.abs(var)) - - def predict_on_structure_en(self): - for n in range(self.structure.nat): - chemenv = env.AtomicEnvironment(self.structure, n, self.gp.cutoffs) - for i in range(3): - force, var = self.gp.predict(chemenv, i + 1) - self.structure.forces[n][i] = float(force) - self.structure.stds[n][i] = np.sqrt(np.abs(var)) - self.local_energies[n] = self.gp.predict_local_energy(chemenv) - - def run_dft(self): - output.write_to_output('\nCalling Quantum Espresso...\n', - self.output_name) - - # calculate DFT forces - forces = qe_util.run_espresso_par(self.qe_input, self.structure, - self.pw_loc, self.no_cpus) - self.structure.forces = forces - - # write wall time of DFT calculation - self.dft_count += 1 - output.write_to_output('QE run complete.\n', self.output_name) - time_curr = time.time() - self.start_time - output.write_to_output('number of DFT calls: %i \n' % self.dft_count, - self.output_name) - output.write_to_output('wall time from start: %.2f s \n' % time_curr, - self.output_name) - - def update_positions(self, new_pos): - self.structure.prev_positions = self.structure.positions - self.structure.positions = new_pos - self.structure.wrap_positions() - - def update_temperature(self, new_pos): - KE, temperature, velocities = \ - md.calculate_temperature(new_pos, self.structure, self.dt, - self.noa) - self.KE = KE - self.temperature = temperature - self.velocities = velocities - - def record_state(self): - output.write_md_config(self.dt, self.curr_step, self.structure, - self.temperature, self.KE, - self.local_energies, self.start_time, - self.output_name, self.dft_step, - self.velocities) diff --git a/flare/modules/flare_parsers/otf_parser.py b/flare/otf_parser.py similarity index 99% rename from flare/modules/flare_parsers/otf_parser.py rename to flare/otf_parser.py index c55ced0f6..d5e67bdbf 100644 --- a/flare/modules/flare_parsers/otf_parser.py +++ b/flare/otf_parser.py @@ -2,7 +2,7 @@ import numpy as np from typing import List, Tuple from flare import gp, env, struc, kernels, otf -from flare.mff.utils import get_l_bound + class OtfAnalysis: def __init__(self, filename, calculate_energy=False): diff --git a/tests/test_parse_otf.py b/tests/test_parse_otf.py index b5af8dfe1..2c2f85916 100644 --- a/tests/test_parse_otf.py +++ b/tests/test_parse_otf.py @@ -1,7 +1,7 @@ import os import sys import numpy as np -from flare.modules.flare_parsers.otf_parser import OtfAnalysis +from flare.otf_parser import OtfAnalysis from flare.kernels import two_plus_three_body, two_plus_three_body_grad from flare.env import AtomicEnvironment