From dccd769920b8e5df46970ba1cf59e9ba8d69a4bd Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Tue, 25 Feb 2020 22:19:53 -0500 Subject: [PATCH 01/21] fix pep8 violations in otf --- flare/otf.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/flare/otf.py b/flare/otf.py index 43727b8b4..a6f5a0e60 100644 --- a/flare/otf.py +++ b/flare/otf.py @@ -86,7 +86,7 @@ def __init__(self, dft_input: str, dt: float, number_of_steps: int, dft_softwarename: str = "qe", no_cpus: int = 1, npool: int = None, mpi: str = "srun", dft_kwargs=None, - store_dft_output: Tuple[Union[str,List[str]],str] = None): + store_dft_output: Tuple[Union[str, List[str]], str] = None): self.dft_input = dft_input self.dt = dt @@ -103,7 +103,6 @@ def __init__(self, dft_input: str, dt: float, number_of_steps: int, positions, species, cell, masses = \ self.dft_module.parse_dft_input(self.dft_input) - self.structure = struc.Structure(cell=cell, species=species, positions=positions, mass_dict=masses, @@ -228,7 +227,6 @@ def run(self): if (self.dft_count-1) < self.freeze_hyps: self.train_gp() - # write gp forces if counter >= self.skip and not self.dft_step: self.update_temperature(new_pos) @@ -243,7 +241,7 @@ def run(self): def run_dft(self): """Calculates DFT forces on atoms in the current structure. - + If OTF has store_dft_output set, then the specified DFT files will be copied with the current date and time prepended in the format 'Year.Month.Day:Hour:Minute:Second:'. @@ -266,7 +264,7 @@ def run_dft(self): time_curr = time.time() - self.start_time self.output.write_to_log('number of DFT calls: %i \n' % self.dft_count) self.output.write_to_log('wall time from start: %.2f s \n' % time_curr) - + # Store DFT outputs in another folder if desired # specified in self.store_dft_output if self.store_dft_output is not None: @@ -332,8 +330,8 @@ def update_temperature(self, new_pos: 'ndarray'): new_pos (np.ndarray): Positions of atoms in the next MD frame. """ KE, temperature, velocities = \ - md.calculate_temperature(new_pos, self.structure, self.dt, - self.noa) + md.calculate_temperature(new_pos, self.structure, self.dt, + self.noa) self.KE = KE self.temperature = temperature self.velocities = velocities From fa63946ab32d97ea710d040edeb631f5c66980e4 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Tue, 25 Feb 2020 22:33:17 -0500 Subject: [PATCH 02/21] add custom module option and documentation --- flare/otf.py | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/flare/otf.py b/flare/otf.py index a6f5a0e60..feb6d9d04 100644 --- a/flare/otf.py +++ b/flare/otf.py @@ -74,6 +74,12 @@ class OTF: single file name, or a list of several. Copied files will be prepended with the date and time with the format 'Year.Month.Day:Hour:Minute:Second:'. + custom_module: A custom module that can be used in place of the DFT + modules available in the FLARE package. The module must contain + two functions: parse_dft_input, which takes a file name (in string + format) as input and returns the positions, species, cell, and + masses of a structure of atoms; and run_dft_par, which takes a + number of DFT related inputs and returns the forces on all atoms. """ def __init__(self, dft_input: str, dt: float, number_of_steps: int, gp: gp.GaussianProcess, dft_loc: str, @@ -86,7 +92,8 @@ def __init__(self, dft_input: str, dt: float, number_of_steps: int, dft_softwarename: str = "qe", no_cpus: int = 1, npool: int = None, mpi: str = "srun", dft_kwargs=None, - store_dft_output: Tuple[Union[str, List[str]], str] = None): + store_dft_output: Tuple[Union[str, List[str]], str] = None, + custom_module=None): self.dft_input = dft_input self.dt = dt @@ -97,7 +104,11 @@ def __init__(self, dft_input: str, dt: float, number_of_steps: int, self.skip = skip self.dft_step = True self.freeze_hyps = freeze_hyps - self.dft_module = dft_software[dft_softwarename] + + if custom_module is not None: + self.dft_module = custom_module + else: + self.dft_module = dft_software[dft_softwarename] # parse input file positions, species, cell, masses = \ From d353d49f0829e7690ebcfd625f18000c4d54420a Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Wed, 26 Feb 2020 11:20:52 -0500 Subject: [PATCH 03/21] add velocity initializationn functions --- flare/md.py | 40 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/flare/md.py b/flare/md.py index 33184141b..8378251fa 100644 --- a/flare/md.py +++ b/flare/md.py @@ -34,3 +34,43 @@ def calculate_temperature(new_pos, structure, dt, noa): temperature = 2 * KE / ((3 * noa - 3) * kb) return KE, temperature, velocities + + +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 + + +def multicomponent_velocities(temperature, masses): + """Draw velocities from the Maxwell-Boltzmann distribution for particles of varying mass. Assumes mass is given in amu.""" + + 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 From d2a45bea5ff28706437a006971e7423d6237b21e Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Wed, 26 Feb 2020 11:29:23 -0500 Subject: [PATCH 04/21] add helper functions for supercell construction --- flare/util.py | 58 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/flare/util.py b/flare/util.py index b768ce1b6..30ef285ab 100644 --- a/flare/util.py +++ b/flare/util.py @@ -8,6 +8,64 @@ import numpy as np + +def get_supercell_positions(sc_size: int, cell: np.ndarray, + positions: np.ndarray): + """Returns the positions of a supercell of atoms, with the number of cells in each direction fixed. + + Args: + sc_size (int): Size of the supercell. + cell (np.ndarray): 3x3 array of cell vectors. + positions (np.ndarray): Positions of atoms in the unit cell. + + Returns: + np.ndarray: Positions of atoms in the supercell. + """ + + 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 np.array(sc_positions) + + +def supercell_custom(cell: np.ndarray, positions: np.ndarray, + size1: int, size2: int, size3: int): + """Returns the positions of a supercell of atoms with a chosen number of cells in each direction. + + Args: + cell (np.ndarray): 3x3 array of cell vectors. + positions (np.ndarray): Positions of atoms in the unit cell. + size1 (int): Number of cells along the first cell vector. + size2 (int): Number of cells along the second cell vector. + size3 (int): Number of cells along the third cell vector. + + Returns: + np.ndarray: Positions of atoms in the supercell. + """ + + 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) + # Dictionary mapping elements to their atomic number (Z) _element_to_Z = {'H': 1, 'He': 2, From 91fce54080c17c659e04854e882adb82793c8efe Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Wed, 26 Feb 2020 11:39:24 -0500 Subject: [PATCH 05/21] add docstrings --- flare/md.py | 27 +++++++++++++++++++++++---- flare/util.py | 10 ++++++---- 2 files changed, 29 insertions(+), 8 deletions(-) diff --git a/flare/md.py b/flare/md.py index 8378251fa..8025194c2 100644 --- a/flare/md.py +++ b/flare/md.py @@ -1,4 +1,5 @@ import numpy as np +from typing import List def update_positions(dt, noa, structure): @@ -37,8 +38,17 @@ def calculate_temperature(new_pos, structure, dt, noa): def get_random_velocities(noa: int, temperature: float, mass: float): - """Draw velocities from Maxwell-Boltzmann distribution. Assumes mass - is given in amu.""" + """Draw velocities from the Maxwell-Boltzmann distribution, assuming a + fixed mass for all particles in amu. + + Args: + noa (int): Number of atoms in the system. + temperature (float): Temperature of the system. + mass (float): Mass of each particle in amu. + + Returns: + np.ndarray: Particle velocities, corrected to give zero center of mass motion. + """ mass_md = mass * 0.000103642695727 kb = 0.0000861733034 @@ -51,8 +61,17 @@ def get_random_velocities(noa: int, temperature: float, mass: float): return corrected_velocities -def multicomponent_velocities(temperature, masses): - """Draw velocities from the Maxwell-Boltzmann distribution for particles of varying mass. Assumes mass is given in amu.""" +def multicomponent_velocities(temperature: float, masses: List[float]): + """Draw velocities from the Maxwell-Boltzmann distribution for particles of + varying mass. + + Args: + temperature (float): Temperature of the system. + masses (List[float]): Particle masses in amu. + + Returns: + np.ndarray: Particle velocities, corrected to give zero center of mass motion. + """ noa = len(masses) velocities = np.zeros((noa, 3)) diff --git a/flare/util.py b/flare/util.py index 30ef285ab..87a276051 100644 --- a/flare/util.py +++ b/flare/util.py @@ -11,7 +11,8 @@ def get_supercell_positions(sc_size: int, cell: np.ndarray, positions: np.ndarray): - """Returns the positions of a supercell of atoms, with the number of cells in each direction fixed. + """Returns the positions of a supercell of atoms, with the number of cells + in each direction fixed. Args: sc_size (int): Size of the supercell. @@ -39,15 +40,16 @@ def get_supercell_positions(sc_size: int, cell: np.ndarray, def supercell_custom(cell: np.ndarray, positions: np.ndarray, size1: int, size2: int, size3: int): - """Returns the positions of a supercell of atoms with a chosen number of cells in each direction. - + """Returns the positions of a supercell of atoms with a chosen number of + cells in each direction. + Args: cell (np.ndarray): 3x3 array of cell vectors. positions (np.ndarray): Positions of atoms in the unit cell. size1 (int): Number of cells along the first cell vector. size2 (int): Number of cells along the second cell vector. size3 (int): Number of cells along the third cell vector. - + Returns: np.ndarray: Positions of atoms in the supercell. """ From 83d0d58af88e19edccac4cfd7bcd9efd03db735c Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Wed, 26 Feb 2020 12:06:16 -0500 Subject: [PATCH 06/21] minor rearrangement of otf inputs --- flare/otf.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/flare/otf.py b/flare/otf.py index feb6d9d04..0e36af654 100644 --- a/flare/otf.py +++ b/flare/otf.py @@ -89,9 +89,8 @@ def __init__(self, dft_input: str, dt: float, number_of_steps: int, calculate_energy: bool = False, output_name: str = 'otf_run', max_atoms_added: int = 1, freeze_hyps: int = 10, rescale_steps: List[int] = [], rescale_temps: List[int] = [], - dft_softwarename: str = "qe", - no_cpus: int = 1, npool: int = None, mpi: str = "srun", - dft_kwargs=None, + dft_softwarename: str = "qe", no_cpus: int = 1, + npool: int = None, mpi: str = "srun", dft_kwargs=None, store_dft_output: Tuple[Union[str, List[str]], str] = None, custom_module=None): From 3b9b542ac640dfce7e308c84313b1c40d4c6ffed Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Wed, 26 Feb 2020 12:13:29 -0500 Subject: [PATCH 07/21] remove random print statements from otf.py and output.py --- flare/otf.py | 1 - flare/output.py | 2 -- 2 files changed, 3 deletions(-) diff --git a/flare/otf.py b/flare/otf.py index 0e36af654..9377ea852 100644 --- a/flare/otf.py +++ b/flare/otf.py @@ -178,7 +178,6 @@ def run(self): self.start_time = time.time() while self.curr_step < self.number_of_steps: - print('curr_step:', self.curr_step) # run DFT and train initial model if first step and DFT is on if self.curr_step == 0 and self.std_tolerance != 0: # call dft and update positions diff --git a/flare/output.py b/flare/output.py index 2a043a1f3..467ba2f5f 100644 --- a/flare/output.py +++ b/flare/output.py @@ -219,8 +219,6 @@ def write_md_config(self, dt, curr_step, structure, string += f'{velocities[i][j]:8.4} ' string += '\n' - print(curr_step) - print(structure.species_labels) self.write_xyz_config(curr_step, structure, dft_step) self.write_xyz(curr_step, structure.stds, structure.species_labels, "std", header) From a9000c6b541c3706a3ddd1d12d5a03f053a54022 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Wed, 26 Feb 2020 14:52:49 -0500 Subject: [PATCH 08/21] report a consistent number of significant figures in otf output files --- flare/output.py | 11 +++++------ flare/util.py | 1 + 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/flare/output.py b/flare/output.py index 467ba2f5f..17c5c0d58 100644 --- a/flare/output.py +++ b/flare/output.py @@ -162,7 +162,6 @@ def write_header(self, cutoffs, kernel_name: str, if self.always_flush: f.flush() - # TO DO: this module should be removed in the future def write_md_config(self, dt, curr_step, structure, temperature, KE, local_energies, start_time, dft_step, velocities): @@ -207,16 +206,16 @@ def write_md_config(self, dt, curr_step, structure, for i in range(len(structure.positions)): string += f'{structure.species_labels[i]} ' for j in range(3): - string += f"{structure.positions[i][j]:8.4} " + string += f"{structure.positions[i][j]:8.4f} " string += '\t' for j in range(3): - string += f"{structure.forces[i][j]:8.4} " + string += f"{structure.forces[i][j]:8.4f} " string += '\t' for j in range(3): - string += f'{structure.stds[i][j]:8.4} ' + string += f'{structure.stds[i][j]:8.4f} ' string += '\t' for j in range(3): - string += f'{velocities[i][j]:8.4} ' + string += f'{velocities[i][j]:8.4f} ' string += '\n' self.write_xyz_config(curr_step, structure, dft_step) @@ -341,7 +340,7 @@ def write_hyps(self, hyp_labels, hyps, start_time, like, like_grad, f.write(f'likelihood gradient: {like_grad}\n') if start_time: time_curr = time.time() - start_time - f.write(f'wall time from start: {time_curr:.2} s \n') + f.write(f'wall time from start: {time_curr:.2f} s \n') if self.always_flush: f.flush() diff --git a/flare/util.py b/flare/util.py index 87a276051..759aed3ef 100644 --- a/flare/util.py +++ b/flare/util.py @@ -68,6 +68,7 @@ def supercell_custom(cell: np.ndarray, positions: np.ndarray, return np.array(sc_positions) + # Dictionary mapping elements to their atomic number (Z) _element_to_Z = {'H': 1, 'He': 2, From 124d6704ab9b41baf825a06e4885ae2d87a6f374 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Wed, 26 Feb 2020 15:41:30 -0500 Subject: [PATCH 09/21] improve output appearance --- flare/output.py | 42 ++++++++++++++++++++++++------------------ 1 file changed, 24 insertions(+), 18 deletions(-) diff --git a/flare/output.py b/flare/output.py index 17c5c0d58..d222bcee2 100644 --- a/flare/output.py +++ b/flare/output.py @@ -151,9 +151,9 @@ def write_header(self, cutoffs, kernel_name: str, # report previous positions headerstring += '\nprevious positions (A):\n' for i in range(len(structure.positions)): - headerstring += f'{structure.species_labels[i]} ' + headerstring += f'{structure.species_labels[i]:5}' for j in range(3): - headerstring += f'{structure.prev_positions[i][j]:10.4} ' + headerstring += f'{structure.prev_positions[i][j]:8.4f}' headerstring += '\n' headerstring += '-' * 80 + '\n' @@ -194,28 +194,34 @@ def write_md_config(self, dt, curr_step, structure, string += f'\nSimulation Time: {(dt * curr_step):.3} ps \n' # Construct Header line - string += 'El Position (A) \t\t\t\t ' + n_space = 24 + string += str.ljust('El', 5) + string += str.center('Position (A)', n_space) + string += ' ' * 4 if not dft_step: - string += 'GP Force (ev/A) ' + string += str.center('GP Force (ev/A)', n_space) + string += ' ' * 4 else: - string += 'DFT Force (ev/A) ' - string += '\t\t\t\t Std. Dev (ev/A) \t' - string += '\t\t\t\t Velocities (A/ps) \n' + string += str.center('DFT Force (ev/A)', n_space) + string += ' ' * 4 + string += str.center('Std. Dev (ev/A)', n_space) + ' ' * 4 + string += str.center('Velocities (A/ps)', n_space) + '\n' # Construct atom-by-atom description for i in range(len(structure.positions)): - string += f'{structure.species_labels[i]} ' + string += f'{structure.species_labels[i]:5}' + # string += '\t' for j in range(3): - string += f"{structure.positions[i][j]:8.4f} " - string += '\t' + string += f'{structure.positions[i][j]:8.4f}' + string += ' ' * 4 for j in range(3): - string += f"{structure.forces[i][j]:8.4f} " - string += '\t' + string += f'{structure.forces[i][j]:8.4f}' + string += ' ' * 4 for j in range(3): - string += f'{structure.stds[i][j]:8.4f} ' - string += '\t' + string += f'{structure.stds[i][j]:8.4f}' + string += ' ' * 4 for j in range(3): - string += f'{velocities[i][j]:8.4f} ' + string += f'{velocities[i][j]:8.4f}' string += '\n' self.write_xyz_config(curr_step, structure, dft_step) @@ -331,12 +337,12 @@ def write_hyps(self, hyp_labels, hyps, start_time, like, like_grad, if hyp_labels is not None: for i, label in enumerate(hyp_labels): - f.write(f'Hyp{i} : {label} = {hyps[i]}\n') + f.write(f'Hyp{i} : {label} = {hyps[i]:.4f}\n') else: for i, hyp in enumerate(hyps): - f.write(f'Hyp{i} : {hyp}\n') + f.write(f'Hyp{i} : {hyp:.4f}\n') - f.write(f'likelihood: {like}\n') + f.write(f'likelihood: {like:.4f}\n') f.write(f'likelihood gradient: {like_grad}\n') if start_time: time_curr = time.time() - start_time From ecf766465781623011868af6b78bc1df81a03494 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Wed, 26 Feb 2020 15:58:48 -0500 Subject: [PATCH 10/21] remove xyz files --- flare/otf.py | 5 ----- flare/output.py | 8 +------- 2 files changed, 1 insertion(+), 12 deletions(-) diff --git a/flare/otf.py b/flare/otf.py index 9377ea852..0a0755cbd 100644 --- a/flare/otf.py +++ b/flare/otf.py @@ -311,9 +311,6 @@ def train_gp(self): """Optimizes the hyperparameters of the current GP model.""" self.gp.train(self.output) - self.output.write_hyps(self.gp.hyp_labels, self.gp.hyps, - self.start_time, - self.gp.likelihood, self.gp.likelihood_gradient) def update_positions(self, new_pos: 'ndarray'): """Performs a Verlet update of the atomic positions. @@ -351,5 +348,3 @@ def record_state(self): self.local_energies, self.start_time, self.dft_step, self.velocities) - self.output.write_xyz_config(self.curr_step, self.structure, - self.dft_step) diff --git a/flare/output.py b/flare/output.py index d222bcee2..16020b182 100644 --- a/flare/output.py +++ b/flare/output.py @@ -36,9 +36,7 @@ def __init__(self, basename: str = 'otf_run', """ self.basename = f"{basename}" self.outfiles = {} - filesuffix = {'log': '.out', 'xyz': '.xyz', - 'fxyz': '-f.xyz', 'hyps': "-hyps.dat", - 'std': '-std.xyz', 'stat': '-stat.dat'} + filesuffix = {'log': '.out', 'hyps': '-hyps.dat'} for filetype in filesuffix.keys(): self.open_new_log(filetype, filesuffix[filetype]) @@ -224,10 +222,6 @@ def write_md_config(self, dt, curr_step, structure, string += f'{velocities[i][j]:8.4f}' string += '\n' - self.write_xyz_config(curr_step, structure, dft_step) - self.write_xyz(curr_step, structure.stds, structure.species_labels, - "std", header) - string += '\n' string += f'temperature: {temperature:.2f} K \n' string += f'kinetic energy: {KE:.6f} eV \n' From 6dac359b14677eb6f740c6e585bde9536030e5bb Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Wed, 26 Feb 2020 16:18:53 -0500 Subject: [PATCH 11/21] comment out lines causing gpfa tests to fail --- flare/output.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/flare/output.py b/flare/output.py index 16020b182..cf01cf5fc 100644 --- a/flare/output.py +++ b/flare/output.py @@ -390,9 +390,9 @@ def write_gp_dft_comparison(self, curr_step, frame, string += '\n' - self.write_xyz_config(curr_step, frame, forces=frame.forces, - stds=frame.stds, forces_2=dft_forces, - dft_step=True) + # self.write_xyz_config(curr_step, frame, forces=frame.forces, + # stds=frame.stds, forces_2=dft_forces, + # dft_step=True) mae = np.mean(error) * 1000 mac = np.mean(np.abs(dft_forces)) * 1000 @@ -433,7 +433,7 @@ def write_gp_dft_comparison(self, curr_step, frame, stat += f' {dt}\n' self.outfiles['log'].write(string) - self.outfiles['stat'].write(stat) + # self.outfiles['stat'].write(stat) if self.always_flush: self.outfiles['log'].flush() From 424c46d9953aba33b491b01ed1b7c940fab3a7ed Mon Sep 17 00:00:00 2001 From: Steven Torrisi Date: Wed, 26 Feb 2020 16:23:32 -0500 Subject: [PATCH 12/21] Add comments to md.py --- flare/md.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/flare/md.py b/flare/md.py index 8025194c2..1422ee908 100644 --- a/flare/md.py +++ b/flare/md.py @@ -49,12 +49,14 @@ def get_random_velocities(noa: int, temperature: float, mass: float): Returns: np.ndarray: Particle velocities, corrected to give zero center of mass motion. """ - + + # Use FLARE mass units (time = ps, length = A, energy = eV) mass_md = mass * 0.000103642695727 kb = 0.0000861733034 std = np.sqrt(kb * temperature / mass_md) velocities = np.random.normal(scale=std, size=(noa, 3)) - + + # Remove center-of-mass motion vel_sum = np.sum(velocities, axis=0) corrected_velocities = velocities - vel_sum / noa @@ -79,6 +81,7 @@ def multicomponent_velocities(temperature: float, masses: List[float]): mom_tot = np.array([0., 0., 0.]) for n, mass in enumerate(masses): + # Convert to FLARE mass units (time = ps, length = A, energy = eV) mass_md = mass * 0.000103642695727 std = np.sqrt(kb * temperature / mass_md) rand_vel = np.random.normal(scale=std, size=(3)) @@ -86,7 +89,7 @@ def multicomponent_velocities(temperature: float, masses: List[float]): mom_curr = rand_vel * mass_md mom_tot += mom_curr - # correct momentum + # Correct momentum, remove center of mass motion mom_corr = mom_tot / noa for n, mass in enumerate(masses): mass_md = mass * 0.000103642695727 From 4fc3bb727fc475cb3038529b9d82c9a001d0f9d1 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Wed, 26 Feb 2020 16:25:37 -0500 Subject: [PATCH 13/21] change dft_softwarename to force_source --- flare/otf.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/flare/otf.py b/flare/otf.py index 0a0755cbd..581927750 100644 --- a/flare/otf.py +++ b/flare/otf.py @@ -89,10 +89,9 @@ def __init__(self, dft_input: str, dt: float, number_of_steps: int, calculate_energy: bool = False, output_name: str = 'otf_run', max_atoms_added: int = 1, freeze_hyps: int = 10, rescale_steps: List[int] = [], rescale_temps: List[int] = [], - dft_softwarename: str = "qe", no_cpus: int = 1, + force_source: Union[str, object] = "qe", no_cpus: int = 1, npool: int = None, mpi: str = "srun", dft_kwargs=None, - store_dft_output: Tuple[Union[str, List[str]], str] = None, - custom_module=None): + store_dft_output: Tuple[Union[str, List[str]], str] = None): self.dft_input = dft_input self.dt = dt @@ -104,10 +103,10 @@ def __init__(self, dft_input: str, dt: float, number_of_steps: int, self.dft_step = True self.freeze_hyps = freeze_hyps - if custom_module is not None: - self.dft_module = custom_module + if isinstance(force_source, str): + self.dft_module = dft_software[force_source] else: - self.dft_module = dft_software[dft_softwarename] + self.dft_module = force_source # parse input file positions, species, cell, masses = \ From b641411c471e0b1d8f6bd5d5b8e51cfff0a89fa2 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Wed, 26 Feb 2020 16:31:06 -0500 Subject: [PATCH 14/21] fix otf unit tests --- tests/test_OTF_cp2k.py | 4 ++-- tests/test_OTF_cp2k_par.py | 2 +- tests/test_OTF_vasp.py | 14 +++++++------- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/tests/test_OTF_cp2k.py b/tests/test_OTF_cp2k.py index 8e5ee617e..79cb26d79 100644 --- a/tests/test_OTF_cp2k.py +++ b/tests/test_OTF_cp2k.py @@ -62,7 +62,7 @@ def test_otf_h2(): otf = OTF(cp2k_input, dt, number_of_steps, gp, dft_loc, std_tolerance_factor, init_atoms=[0], calculate_energy=True, max_atoms_added=1, - dft_softwarename="cp2k", + force_source="cp2k", output_name='h2_otf_cp2k') otf.run() @@ -111,7 +111,7 @@ def test_otf_al(): std_tolerance_factor, init_atoms=[0], calculate_energy=True, output_name='al_otf_cp2k', freeze_hyps=freeze_hyps, skip=5, - dft_softwarename="cp2k", + force_source="cp2k", max_atoms_added=max_atoms_added) otf.run() diff --git a/tests/test_OTF_cp2k_par.py b/tests/test_OTF_cp2k_par.py index a2514b089..a8e38bb47 100644 --- a/tests/test_OTF_cp2k_par.py +++ b/tests/test_OTF_cp2k_par.py @@ -66,7 +66,7 @@ def test_otf_h2_par(): otf = OTF(cp2k_input, dt, number_of_steps, gp, dft_loc, std_tolerance_factor, init_atoms=[0], calculate_energy=True, max_atoms_added=1, - dft_softwarename="cp2k", + force_source="cp2k", no_cpus=2, par=True, mpi="mpi", output_name='h2_otf_cp2k_par') diff --git a/tests/test_OTF_vasp.py b/tests/test_OTF_vasp.py index aafa0c2c8..618786e90 100644 --- a/tests/test_OTF_vasp.py +++ b/tests/test_OTF_vasp.py @@ -35,17 +35,17 @@ def test_otf_h2(): energy_force_kernel = en.two_body_force_en gp = GaussianProcess(kernel=kernel, - kernel_grad=kernel_grad, - hyps=hyps, - cutoffs=cutoffs, - hyp_labels=hyp_labels, - energy_force_kernel=energy_force_kernel, - maxiter=50) + kernel_grad=kernel_grad, + hyps=hyps, + cutoffs=cutoffs, + hyp_labels=hyp_labels, + energy_force_kernel=energy_force_kernel, + maxiter=50) otf = OTF(vasp_input, dt, number_of_steps, gp, dft_loc, std_tolerance_factor, init_atoms=[0], calculate_energy=True, max_atoms_added=1, - no_cpus=1, dft_softwarename='vasp', + no_cpus=1, force_source='vasp', output_name='h2_otf_vasp') otf.run() From 6eccb622d562a04ba901e1589d1dd6a237a0f335 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Wed, 26 Feb 2020 16:36:17 -0500 Subject: [PATCH 15/21] fix otf docstring; remove offending mgp unit test for now --- flare/otf.py | 16 ++++----- tests/test_mgp.py | 91 +++++++++++++++++++++++------------------------ 2 files changed, 53 insertions(+), 54 deletions(-) diff --git a/flare/otf.py b/flare/otf.py index 581927750..0d1be9b01 100644 --- a/flare/otf.py +++ b/flare/otf.py @@ -54,8 +54,14 @@ class OTF: velocities of the atoms are rescaled. Defaults to []. rescale_temps (List[int], optional): List of rescaled temperatures. Defaults to []. - dft_softwarename (str, optional): DFT code used to calculate - ab initio forces during training. Defaults to "qe". + force_source (Union[str, object], optional): DFT code used to calculate + ab initio forces during training. A custom module can be used here + in place of the DFT modules available in the FLARE package. The + module must contain two functions: parse_dft_input, which takes a + file name (in string format) as input and returns the positions, + species, cell, and masses of a structure of atoms; and run_dft_par, + which takes a number of DFT related inputs and returns the forces + on all atoms. Defaults to "qe". no_cpus (int, optional): Number of cpus used during training. Defaults to 1. npool (int, optional): Number of k-point pools for DFT @@ -74,12 +80,6 @@ class OTF: single file name, or a list of several. Copied files will be prepended with the date and time with the format 'Year.Month.Day:Hour:Minute:Second:'. - custom_module: A custom module that can be used in place of the DFT - modules available in the FLARE package. The module must contain - two functions: parse_dft_input, which takes a file name (in string - format) as input and returns the positions, species, cell, and - masses of a structure of atoms; and run_dft_par, which takes a - number of DFT related inputs and returns the forces on all atoms. """ def __init__(self, dft_input: str, dt: float, number_of_steps: int, gp: gp.GaussianProcess, dft_loc: str, diff --git a/tests/test_mgp.py b/tests/test_mgp.py index 6ff38c3b1..552b88cc2 100644 --- a/tests/test_mgp.py +++ b/tests/test_mgp.py @@ -162,49 +162,48 @@ def test_2_plus_3_body(otf_object, structure, params): -# ASSUMPTION: You have a Lammps executable with the mgp pair style with $lmp -# as the corresponding environment variable. -@pytest.mark.skipif(not os.environ.get('lmp', - False), reason='lmp not found ' - 'in environment: Please install LAMMPS ' - 'and set the $lmp env. ' - 'variable to point to the executatble.') -def test_lammps(otf_object, structure): - atom_types = [1, 2] - atom_masses = [108, 127] - atom_species = [1, 2] * 27 - - # create data file - lammps_location = 'AgI_Molten_15.txt' - data_file_name = 'tmp.data' - data_text = lammps_calculator.lammps_dat(structure, atom_types, - atom_masses, atom_species) - lammps_calculator.write_text(data_file_name, data_text) - - # create lammps input - style_string = 'mgp' - coeff_string = '* * {} Ag I yes yes'.format(lammps_location) - lammps_executable = '$lmp' - dump_file_name = 'tmp.dump' - input_file_name = 'tmp.in' - output_file_name = 'tmp.out' - input_text = \ - lammps_calculator.generic_lammps_input(data_file_name, style_string, - coeff_string, dump_file_name) - lammps_calculator.write_text(input_file_name, input_text) - - lammps_calculator.run_lammps(lammps_executable, input_file_name, - output_file_name) - - lammps_forces = lammps_calculator.lammps_parser(dump_file_name) - - # check that lammps agrees with gp to within 1 meV/A - forces = otf_object.force_list[-1] - assert(np.abs(lammps_forces[0, 1] - forces[0, 1]) < 1e-3) - - os.system('rm tmp.in tmp.out tmp.dump tmp.data' - ' log.lammps') - os.system('rm '+lammps_location) - os.system('rm grid3*.npy') - os.system('rm -r kv3') - +# # ASSUMPTION: You have a Lammps executable with the mgp pair style with $lmp +# # as the corresponding environment variable. +# @pytest.mark.skipif(not os.environ.get('lmp', +# False), reason='lmp not found ' +# 'in environment: Please install LAMMPS ' +# 'and set the $lmp env. ' +# 'variable to point to the executatble.') +# def test_lammps(otf_object, structure): +# atom_types = [1, 2] +# atom_masses = [108, 127] +# atom_species = [1, 2] * 27 + +# # create data file +# lammps_location = 'AgI_Molten_15.txt' +# data_file_name = 'tmp.data' +# data_text = lammps_calculator.lammps_dat(structure, atom_types, +# atom_masses, atom_species) +# lammps_calculator.write_text(data_file_name, data_text) + +# # create lammps input +# style_string = 'mgp' +# coeff_string = '* * {} Ag I yes yes'.format(lammps_location) +# lammps_executable = '$lmp' +# dump_file_name = 'tmp.dump' +# input_file_name = 'tmp.in' +# output_file_name = 'tmp.out' +# input_text = \ +# lammps_calculator.generic_lammps_input(data_file_name, style_string, +# coeff_string, dump_file_name) +# lammps_calculator.write_text(input_file_name, input_text) + +# lammps_calculator.run_lammps(lammps_executable, input_file_name, +# output_file_name) + +# lammps_forces = lammps_calculator.lammps_parser(dump_file_name) + +# # check that lammps agrees with gp to within 1 meV/A +# forces = otf_object.force_list[-1] +# assert(np.abs(lammps_forces[0, 1] - forces[0, 1]) < 1e-3) + +# os.system('rm tmp.in tmp.out tmp.dump tmp.data' +# ' log.lammps') +# os.system('rm '+lammps_location) +# os.system('rm grid3*.npy') +# os.system('rm -r kv3') From 7a2183e4493a5f1c0b38ec93e97c35c6b3b89487 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Wed, 26 Feb 2020 16:37:54 -0500 Subject: [PATCH 16/21] move velocity helper functions to util.py --- flare/md.py | 61 --------------------------------------------------- flare/util.py | 61 +++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 61 insertions(+), 61 deletions(-) diff --git a/flare/md.py b/flare/md.py index 1422ee908..cab5451a9 100644 --- a/flare/md.py +++ b/flare/md.py @@ -35,64 +35,3 @@ def calculate_temperature(new_pos, structure, dt, noa): temperature = 2 * KE / ((3 * noa - 3) * kb) return KE, temperature, velocities - - -def get_random_velocities(noa: int, temperature: float, mass: float): - """Draw velocities from the Maxwell-Boltzmann distribution, assuming a - fixed mass for all particles in amu. - - Args: - noa (int): Number of atoms in the system. - temperature (float): Temperature of the system. - mass (float): Mass of each particle in amu. - - Returns: - np.ndarray: Particle velocities, corrected to give zero center of mass motion. - """ - - # Use FLARE mass units (time = ps, length = A, energy = eV) - mass_md = mass * 0.000103642695727 - kb = 0.0000861733034 - std = np.sqrt(kb * temperature / mass_md) - velocities = np.random.normal(scale=std, size=(noa, 3)) - - # Remove center-of-mass motion - vel_sum = np.sum(velocities, axis=0) - corrected_velocities = velocities - vel_sum / noa - - return corrected_velocities - - -def multicomponent_velocities(temperature: float, masses: List[float]): - """Draw velocities from the Maxwell-Boltzmann distribution for particles of - varying mass. - - Args: - temperature (float): Temperature of the system. - masses (List[float]): Particle masses in amu. - - Returns: - np.ndarray: Particle velocities, corrected to give zero center of mass motion. - """ - - noa = len(masses) - velocities = np.zeros((noa, 3)) - kb = 0.0000861733034 - mom_tot = np.array([0., 0., 0.]) - - for n, mass in enumerate(masses): - # Convert to FLARE mass units (time = ps, length = A, energy = eV) - 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, remove center of mass motion - 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/util.py b/flare/util.py index 759aed3ef..96aaec224 100644 --- a/flare/util.py +++ b/flare/util.py @@ -9,6 +9,67 @@ import numpy as np +def get_random_velocities(noa: int, temperature: float, mass: float): + """Draw velocities from the Maxwell-Boltzmann distribution, assuming a + fixed mass for all particles in amu. + + Args: + noa (int): Number of atoms in the system. + temperature (float): Temperature of the system. + mass (float): Mass of each particle in amu. + + Returns: + np.ndarray: Particle velocities, corrected to give zero center of mass motion. + """ + + # Use FLARE mass units (time = ps, length = A, energy = eV) + mass_md = mass * 0.000103642695727 + kb = 0.0000861733034 + std = np.sqrt(kb * temperature / mass_md) + velocities = np.random.normal(scale=std, size=(noa, 3)) + + # Remove center-of-mass motion + vel_sum = np.sum(velocities, axis=0) + corrected_velocities = velocities - vel_sum / noa + + return corrected_velocities + + +def multicomponent_velocities(temperature: float, masses: List[float]): + """Draw velocities from the Maxwell-Boltzmann distribution for particles of + varying mass. + + Args: + temperature (float): Temperature of the system. + masses (List[float]): Particle masses in amu. + + Returns: + np.ndarray: Particle velocities, corrected to give zero center of mass motion. + """ + + noa = len(masses) + velocities = np.zeros((noa, 3)) + kb = 0.0000861733034 + mom_tot = np.array([0., 0., 0.]) + + for n, mass in enumerate(masses): + # Convert to FLARE mass units (time = ps, length = A, energy = eV) + 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, remove center of mass motion + mom_corr = mom_tot / noa + for n, mass in enumerate(masses): + mass_md = mass * 0.000103642695727 + velocities[n] -= mom_corr / mass_md + + return velocities + + def get_supercell_positions(sc_size: int, cell: np.ndarray, positions: np.ndarray): """Returns the positions of a supercell of atoms, with the number of cells From 1087a738b5d87db6862f7a17b5e7e89040ddf16d Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Wed, 26 Feb 2020 16:50:05 -0500 Subject: [PATCH 17/21] fix david's pep8 catch --- flare/md.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/flare/md.py b/flare/md.py index cab5451a9..0c327ac9c 100644 --- a/flare/md.py +++ b/flare/md.py @@ -1,6 +1,7 @@ -import numpy as np from typing import List +import numpy as np + def update_positions(dt, noa, structure): dtdt = dt ** 2 From 1345b4f71c80f0d19146b5409fcb6fd78e33a900 Mon Sep 17 00:00:00 2001 From: Jonathan Vandermause Date: Wed, 26 Feb 2020 17:03:06 -0500 Subject: [PATCH 18/21] report hyps in otf output file --- flare/otf.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/flare/otf.py b/flare/otf.py index 0d1be9b01..d050bd801 100644 --- a/flare/otf.py +++ b/flare/otf.py @@ -310,6 +310,9 @@ def train_gp(self): """Optimizes the hyperparameters of the current GP model.""" self.gp.train(self.output) + self.output.write_hyps(self.gp.hyp_labels, self.gp.hyps, + self.start_time, + self.gp.likelihood, self.gp.likelihood_gradient) def update_positions(self, new_pos: 'ndarray'): """Performs a Verlet update of the atomic positions. From e1aa89ebb34332b4ca0c6804d6f388a86e19f754 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Wed, 26 Feb 2020 17:13:58 -0500 Subject: [PATCH 19/21] change mgp -> mgpf to adapt to new lmp exe --- tests/test_mgp.py | 90 +++++++++++++++++++++++------------------------ 1 file changed, 45 insertions(+), 45 deletions(-) diff --git a/tests/test_mgp.py b/tests/test_mgp.py index 552b88cc2..609422a1f 100644 --- a/tests/test_mgp.py +++ b/tests/test_mgp.py @@ -162,48 +162,48 @@ def test_2_plus_3_body(otf_object, structure, params): -# # ASSUMPTION: You have a Lammps executable with the mgp pair style with $lmp -# # as the corresponding environment variable. -# @pytest.mark.skipif(not os.environ.get('lmp', -# False), reason='lmp not found ' -# 'in environment: Please install LAMMPS ' -# 'and set the $lmp env. ' -# 'variable to point to the executatble.') -# def test_lammps(otf_object, structure): -# atom_types = [1, 2] -# atom_masses = [108, 127] -# atom_species = [1, 2] * 27 - -# # create data file -# lammps_location = 'AgI_Molten_15.txt' -# data_file_name = 'tmp.data' -# data_text = lammps_calculator.lammps_dat(structure, atom_types, -# atom_masses, atom_species) -# lammps_calculator.write_text(data_file_name, data_text) - -# # create lammps input -# style_string = 'mgp' -# coeff_string = '* * {} Ag I yes yes'.format(lammps_location) -# lammps_executable = '$lmp' -# dump_file_name = 'tmp.dump' -# input_file_name = 'tmp.in' -# output_file_name = 'tmp.out' -# input_text = \ -# lammps_calculator.generic_lammps_input(data_file_name, style_string, -# coeff_string, dump_file_name) -# lammps_calculator.write_text(input_file_name, input_text) - -# lammps_calculator.run_lammps(lammps_executable, input_file_name, -# output_file_name) - -# lammps_forces = lammps_calculator.lammps_parser(dump_file_name) - -# # check that lammps agrees with gp to within 1 meV/A -# forces = otf_object.force_list[-1] -# assert(np.abs(lammps_forces[0, 1] - forces[0, 1]) < 1e-3) - -# os.system('rm tmp.in tmp.out tmp.dump tmp.data' -# ' log.lammps') -# os.system('rm '+lammps_location) -# os.system('rm grid3*.npy') -# os.system('rm -r kv3') +# ASSUMPTION: You have a Lammps executable with the mgp pair style with $lmp +# as the corresponding environment variable. +@pytest.mark.skipif(not os.environ.get('lmp', + False), reason='lmp not found ' + 'in environment: Please install LAMMPS ' + 'and set the $lmp env. ' + 'variable to point to the executatble.') +def test_lammps(otf_object, structure): + atom_types = [1, 2] + atom_masses = [108, 127] + atom_species = [1, 2] * 27 + + # create data file + lammps_location = 'AgI_Molten_15.txt' + data_file_name = 'tmp.data' + data_text = lammps_calculator.lammps_dat(structure, atom_types, + atom_masses, atom_species) + lammps_calculator.write_text(data_file_name, data_text) + + # create lammps input + style_string = 'mgpf' + coeff_string = '* * {} Ag I yes yes'.format(lammps_location) + lammps_executable = '$lmp' + dump_file_name = 'tmp.dump' + input_file_name = 'tmp.in' + output_file_name = 'tmp.out' + input_text = \ + lammps_calculator.generic_lammps_input(data_file_name, style_string, + coeff_string, dump_file_name) + lammps_calculator.write_text(input_file_name, input_text) + + lammps_calculator.run_lammps(lammps_executable, input_file_name, + output_file_name) + + lammps_forces = lammps_calculator.lammps_parser(dump_file_name) + + # check that lammps agrees with gp to within 1 meV/A + forces = otf_object.force_list[-1] + assert(np.abs(lammps_forces[0, 1] - forces[0, 1]) < 1e-3) + + os.system('rm tmp.in tmp.out tmp.dump tmp.data' + ' log.lammps') + os.system('rm '+lammps_location) + os.system('rm grid3*.npy') + os.system('rm -r kv3') From ef94823034612324771ea3b0b1fd340242c5dd58 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Wed, 26 Feb 2020 17:15:32 -0500 Subject: [PATCH 20/21] close issue #154 --- tests/test_mgp.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/test_mgp.py b/tests/test_mgp.py index 609422a1f..ee6ad9020 100644 --- a/tests/test_mgp.py +++ b/tests/test_mgp.py @@ -182,7 +182,8 @@ def test_lammps(otf_object, structure): lammps_calculator.write_text(data_file_name, data_text) # create lammps input - style_string = 'mgpf' + style_string = 'mgpf' # use mgpf for the new lmp executable pair_style, + # because mgp is the energy version coeff_string = '* * {} Ag I yes yes'.format(lammps_location) lammps_executable = '$lmp' dump_file_name = 'tmp.dump' From 42483435858e13f364da92c778c23a6870513952 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Wed, 26 Feb 2020 17:54:01 -0500 Subject: [PATCH 21/21] change lmp --- tests/test_mgp.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_mgp.py b/tests/test_mgp.py index ee6ad9020..934e35fd9 100644 --- a/tests/test_mgp.py +++ b/tests/test_mgp.py @@ -182,7 +182,7 @@ def test_lammps(otf_object, structure): lammps_calculator.write_text(data_file_name, data_text) # create lammps input - style_string = 'mgpf' # use mgpf for the new lmp executable pair_style, + style_string = 'mgpf' # use mgpf for the new lammps executable pair_style, # because mgp is the energy version coeff_string = '* * {} Ag I yes yes'.format(lammps_location) lammps_executable = '$lmp'