From 2e9839b51311dd5257b34cfb222eff5e1180d5be Mon Sep 17 00:00:00 2001 From: Steven T Date: Sat, 8 Sep 2018 01:10:12 -0400 Subject: [PATCH 1/9] Added new tests for punchout as well as for OTF --- src/punchout.py | 0 tests/test_GaussianProcess.py | 3 +- tests/test_OTF.py | 90 ++++++++++++++----- tests/test_files/qe_input_1.in | 10 +-- tests/test_punchout.py | 156 +++++++++++++++++++++++++++++++++ 5 files changed, 230 insertions(+), 29 deletions(-) create mode 100644 src/punchout.py create mode 100644 tests/test_punchout.py diff --git a/src/punchout.py b/src/punchout.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/test_GaussianProcess.py b/tests/test_GaussianProcess.py index 571983d64..a9d931f9e 100644 --- a/tests/test_GaussianProcess.py +++ b/tests/test_GaussianProcess.py @@ -9,7 +9,8 @@ import pytest import numpy as np - +import sys +sys.path.append('../src') from gp import GaussianProcess, minus_like_hyp from env import ChemicalEnvironment from struc import Structure diff --git a/tests/test_OTF.py b/tests/test_OTF.py index 3c59d5095..38ef7097d 100644 --- a/tests/test_OTF.py +++ b/tests/test_OTF.py @@ -15,17 +15,28 @@ import sys import numpy as np sys.path.append("../src") -from otf import parse_qe_input, parse_qe_forces, OTF +from otf import parse_qe_input, parse_qe_forces, OTF, edit_qe_input_positions,\ + run_espresso from gp import GaussianProcess +from struc import Structure -def fake_espresso(noa: int): - """ Returns a list of random forces. +def fake_espresso(qe_input: str, structure: Structure): """ - forces = [np.random.normal(loc=0,scale=1,size=3) for n in range(noa)] - print(forces) - return forces + Returns a list of random forces. Takes the same argument as real ESPRESSO + in order to allow for substitution. + """ + noa = len(structure.positions) + return [np.random.normal(loc=0,scale=1,size=3) for _ in range(noa)] +def fake_predict_on_structure(self, structure): + """ + Substitutes in the predict_on_structure method of GaussianProcess + """ + structure.forces = [np.random.randn(3) for n in structure.positions] + structure.stds = [np.random.randn(3) for n in structure.positions] + + return structure class Fake_GP(GaussianProcess): """ @@ -35,6 +46,10 @@ class Fake_GP(GaussianProcess): def __init__(self, kernel): super(GaussianProcess, self).__init__() + self.sigma_f = 0 + self.length_scale = 0 + self.sigma_n = 0 + pass def train(self): @@ -49,24 +64,15 @@ def update_db(self, structure, forces): """ pass - def predict(self, structure, _): + def predict(self, chemenv, i): """ Substitutes in the predict method of GaussianProcess """ - structure.forces = [np.random.randn(3) for n in range(structure.nat)] - structure.stds = [np.random.randn(3) for n in range(structure.nat)] + #structure.forces = [np.random.randn(3) for n in range(structure.nat)] + #structure.stds = [np.random.randn(3) for n in range(structure.nat)] - return structure - - def predict_on_structure(self, structure): - """ - Substitutes in the predict_on_structure method of GaussianProcess - """ - - structure.forces = [np.random.randn(3) for n in structure.positions] - structure.stds = [np.random.randn(3) for n in structure.positions] + return chemenv - return structure def cleanup_otf_run(): @@ -74,6 +80,7 @@ def cleanup_otf_run(): os.system('rm pwscf.out') os.system('rm pwscf.wfc') os.system('rm pwscf.save') + os.system('rm otf_run.out') # ------------------------------------------------------ # fixtures @@ -155,7 +162,6 @@ def test_cell_parsing(qe_input, exp_cell): ) def test_cell_parsing(qe_input, mass_dict): positions, species, cell, masses = parse_qe_input(qe_input) - print(masses) assert masses == mass_dict @pytest.mark.parametrize('qe_output,exp_forces', @@ -178,7 +184,8 @@ def test_espresso_calling_1(test_otf_engine_1): False), 'PWSCF_COMMAND not found ' \ 'in environment' - forces = test_otf_engine_1.run_espresso() + forces = run_espresso(test_otf_engine_1.qe_input, + test_otf_engine_1.structure) assert isinstance(forces, list) assert len(forces) == len(test_otf_engine_1.structure.forces) @@ -191,6 +198,30 @@ def test_espresso_calling_1(test_otf_engine_1): cleanup_otf_run() +def test_espresso_input_edit(): + """ + Load a structure in from qe_input_1, change the position and cell, + then edit and re-parse + :return: + """ + os.system('cp test_files/qe_input_1.in .') + positions, species, cell, masses = parse_qe_input('./qe_input_1.in') + struc = Structure(cell,species,positions,masses) + + struc.vec1 += np.random.randn(3) + struc.positions[0]+= np.random.randn(3) + + edit_qe_input_positions('./qe_input_1.in', structure=struc) + + positions, species, cell, masses = parse_qe_input('./qe_input_1.in') + + assert np.equal(positions[0],struc.positions[0]).all() + assert np.equal(struc.vec1,cell[0, :]).all() + + os.system('rm ./qe_input_1.in') + + + # ------------------------------------------------------ # test otf methods # ------------------------------------------------------ @@ -215,5 +246,18 @@ def test_update_1(test_otf_engine_1): for i,pos in enumerate(test_otf_engine_1.structure.positions): assert np.isclose(pos,target_positions[i],rtol=1e-6).all() - - +# ------------------------------------------------------ +# test otf runs +# ------------------------------------------------------ +# Under development +""" +def test_otf_1(): + os.system('cp ./test_files/qe_input_1.in ./pwscf.in') + + otf = OTF('./pwscf.in', .1, 2, kernel='two_body', + cutoff=10) + otf.run_espresso = fake_espresso + otf.gp = Fake_GP(kernel='') + otf.run() + cleanup_otf_run() +""" \ No newline at end of file diff --git a/tests/test_files/qe_input_1.in b/tests/test_files/qe_input_1.in index 8c1342c6b..ccd22336c 100644 --- a/tests/test_files/qe_input_1.in +++ b/tests/test_files/qe_input_1.in @@ -21,14 +21,14 @@ &IONS / &CELL -/ -K_POINTS {gamma} +/ ATOMIC_SPECIES H 1.0 H.pbe-kjpaw.UPF CELL_PARAMETERS {angstrom} - 5 0 0 - 0 5 0 - 0 0 5 +5.0 0.0 0.0 +0.0 5.0 0.0 +0.0 0.0 5.0 ATOMIC_POSITIONS {angstrom} H 2.51857 2.5 2.5 H 4.48143 2.5 2.5 +K_POINTS {gamma} diff --git a/tests/test_punchout.py b/tests/test_punchout.py new file mode 100644 index 000000000..a968db37f --- /dev/null +++ b/tests/test_punchout.py @@ -0,0 +1,156 @@ +#!/usr/bin/env python3 +# pylint: disable=redefined-outer-name + +"""" Punchout test suite based on py.test + +Steven Torrisi +""" + +import pytest + +import numpy as np + +from punchout import is_within_d_box, punchout +from test_GaussianProcess import get_random_structure +from struc import Structure + + +# ----------------------- +# Test helper methods +# ----------------------- + +def test_d_within_box(): + """ + Run with a few examples + :return: + """ + + r1 = np.array([0, 0, 0]) + r2 = np.array([0, 0, .5]) + r3 = np.array([.5, .5, .5]) + r4 = np.array([.6, .5, .5]) + r5 = np.array([1.0, 2.0, 3.0]) + + assert is_within_d_box(r1, r2, d=1) + assert is_within_d_box(r1, -r2, d=1) + assert not is_within_d_box(r1, r2, d=.5) + assert not is_within_d_box(r1, r2, d=.1) + + assert is_within_d_box(r1, r3, 1) + assert is_within_d_box(r2, r3, 1) + + assert not is_within_d_box(r1, r4, 1) + assert is_within_d_box(r1, r5, 6.0) + assert not is_within_d_box(r1, r5, 3.0) + + assert is_within_d_box(r3, r4, .2) + assert is_within_d_box(r3, r4, .5) + assert not is_within_d_box(r3, r4, .1) + + # Run for a hundred random pairs + for n in range(100): + rand1 = np.random.randn(3) + rand2 = np.random.randn(3) + + radius = np.linalg.norm(rand1 - rand2) + assert is_within_d_box(rand1, rand2, radius * 2) + assert not is_within_d_box(rand1, rand2, radius) + + +# ----------------------- +# Test punchout +# ----------------------- + +@pytest.fixture(scope='module') +def test_punchout_random_struc(): + cell = 10 * np.eye(3) + species = ['A', 'B', 'C', 'D', 'E'] * 3 + mass_dict = {spec: np.random.randn(1) for spec in species} + noa = 15 + struct, _ = get_random_structure(cell, species, cutoff=5, noa=noa) + struct.mass_dict = mass_dict + + print(struct) + target_atom = np.random.randint(0, noa) + d = np.random.uniform(1, 8) + punchstruc = punchout(struct, target_atom, d=d) + assert isinstance(punchstruc, Structure) + + yield (punchstruc, d) + + del punchstruc + + +def test_punchout_boxed(test_punchout_random_struc): + punchstruc = test_punchout_random_struc[0] + d = test_punchout_random_struc[1] + + for pos1 in punchstruc.positions: + for pos2 in punchstruc.positions: + assert is_within_d_box(pos1, pos2, d) + + +def test_punchout_edges_1(): + """ + Test a few hand-defined edge cases to ensure punched out atoms near the + edges correctly have their periodic neighbors + :return: + """ + species = ['A'] + cell = np.eye(3) + positions = [.5 * np.ones(3)] + struc = Structure(cell, species, positions, cutoff=1, + mass_dict={'A': 1}) + + p1 = punchout(structure=struc, atom=0, d=.25,center=False) + assert len(p1.positions) == 1 + assert np.equal(p1.positions, .5 * np.ones(3)).all() + + p2 = punchout(structure=struc, atom=0, d=2.,center=False) + assert len(p2.positions) == 27 + + p3 = punchout(structure=struc, atom=0, d=1.0,center=True) + assert np.equal(p3.positions,np.zeros(3)).all() + +def test_punchout_edges_2(): + species = ['A'] * 2 + cell = np.eye(3) + positions = [np.zeros(3), .25 * np.ones(3)] + struc = Structure(cell, species, positions, cutoff=1, + mass_dict={'A': 1}) + p1 = punchout(structure=struc, atom=0, d=.01,center=False) + print(p1) + assert len(p1.positions) == 1 + assert np.equal(p1.positions, np.zeros(3)).all() + + p2 = punchout(structure=struc, atom=0, d=.5,center=False) + assert np.equal(p2.positions, positions).all() + + p3 = punchout(structure=struc, atom=1, d=1.,center=False) + assert np.equal(p3.positions, positions).all() + + +def test_punchout_edges_3(): + species = ['A'] * 2 + cell = np.eye(3) + positions = [np.zeros(3), .9 * np.ones(3)] + struc = Structure(cell, species, positions, cutoff=1, + mass_dict={'A': 1}) + p1 = punchout(structure=struc, atom=0, d=.01,center=False) + assert len(p1.positions) == 1 + assert np.equal(p1.positions, np.zeros(3)).all() + + p2 = punchout(structure=struc, atom=0, d=.5,center=False) + assert np.isclose(p2.positions, [np.zeros(3), -.1 * np.ones(3)]).all() + + p3 = punchout(structure=struc, atom=1, d=.5,center=True) + assert np.isclose(p3.positions, [.1 * np.ones(3),np.zeros(3) ]).all() + + p4 = punchout(structure=struc, atom=1, d=1.,center=False) + assert np.isclose(p4.positions, [np.ones(3), .9*np.ones(3)]).all() + + + + + + From 4d71925e721dac63fc5ca27e36e858a47100a509 Mon Sep 17 00:00:00 2001 From: Steven T Date: Sat, 8 Sep 2018 01:10:42 -0400 Subject: [PATCH 2/9] Various refactors to otf; added default argument for punchout_d, moved espresso methods out of OTF (still in otf for now) --- src/otf.py | 170 +++++++++++++++++++++++------------- src/punchout.py | 228 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 336 insertions(+), 62 deletions(-) diff --git a/src/otf.py b/src/otf.py index 3cb48b988..acf833aea 100644 --- a/src/otf.py +++ b/src/otf.py @@ -14,11 +14,11 @@ from math import sqrt from gp import GaussianProcess from env import ChemicalEnvironment - +from punchout import punchout class OTF(object): def __init__(self, qe_input: str, dt: float, number_of_steps: int, - kernel: str, cutoff: float): + kernel: str, cutoff: float, punchout_d: float=None): """ On-the-fly learning engine, containing methods to run OTF calculation @@ -27,12 +27,16 @@ def __init__(self, qe_input: str, dt: float, number_of_steps: int, :param number_of_steps: int, Number of steps :param kernel: Type of kernel for GP regression model :param cutoff: Cutoff radius for kernel + :param punchout_d: Box distance around a high-uncertainty atom to + punch out """ self.qe_input = qe_input self.dt = dt self.Nsteps = number_of_steps self.gp = GaussianProcess(kernel) + self.cutoff = cutoff + self.punchout_d = punchout_d positions, species, cell, masses = parse_qe_input(self.qe_input) self.structure = Structure(lattice=cell, species=species, @@ -40,6 +44,7 @@ def __init__(self, qe_input: str, dt: float, number_of_steps: int, mass_dict=masses) self.curr_step = 0 + self.train_structure = None # Create blank output file with time header and structure information with open('otf_run.out', 'w') as f: @@ -88,15 +93,24 @@ def run_and_train(self): # Run espresso and write out results self.write_to_output('=' * 20 + '\n' + 'Calling QE... ') - forces = self.run_espresso() + + if self.punchout_d == None: + self.train_structure = self.structure + elif self.train_structure is None: + target_atom = np.random.randint(0,len(self.structure.positions)) + self.train_structure = punchout(self.structure,target_atom, + d=self.punchout_d) + + + forces = run_espresso(self.qe_input,self.train_structure) self.write_to_output('Done.\n') # Write input positions and force results qe_strings = 'Resultant Positions and Forces:\n' - for n in range(self.structure.nat): - qe_strings += self.structure.species[n] + ': ' + for n in range(self.train_structure.nat): + qe_strings += self.train_structure.species[n] + ': ' for i in range(3): - qe_strings += '%.3f ' % self.structure.positions[n][i] + qe_strings += '%.3f ' % self.train_structure.positions[n][i] qe_strings += '\t ' for i in range(3): qe_strings += '%.3f ' % forces[n][i] @@ -105,7 +119,7 @@ def run_and_train(self): # Update hyperparameters and write results self.write_to_output('Updating database hyperparameters...\n') - self.gp.update_db(self.structure, forces) + self.gp.update_db(self.train_structure, forces) self.gp.train() self.write_to_output('New GP Hyperparameters:\n' + 'Signal variance: \t' + @@ -133,53 +147,6 @@ def update_positions(self): self.structure.prev_positions[i] = np.array(temp_pos) - def run_espresso(self): - """ - Calls quantum espresso from input located at self.qe_input - - :return: List [nparray] List of forces - """ - - self.write_to_output('') - self.write_qe_input_positions() - - pw_loc = os.environ.get('PWSCF_COMMAND', 'pwscf.x') - - qe_command = '{0} < {1} > {2}'.format(pw_loc, self.qe_input, - 'pwscf.out') - - os.system(qe_command) - - return parse_qe_forces('pwscf.out') - - def write_qe_input_positions(self): - """ - Write the current configuration of the OTF structure to the - qe input file - """ - - with open(self.qe_input, 'r') as f: - lines = f.readlines() - - file_pos_index = None - for i, line in enumerate(lines): - if 'ATOMIC_POSITIONS' in line: - file_pos_index = int(i + 1) - assert file_pos_index is not None, 'Failed to find positions in input' - - for pos_index, line_index in enumerate( - range(file_pos_index, len(lines))): - pos_string = ' '.join([self.structure.species[pos_index], - str(self.structure.positions[pos_index][0]), - str(self.structure.positions[pos_index][1]), - str(self.structure.positions[pos_index][ - 2])]) - lines[line_index] = str(pos_string + '\n') - - with open(self.qe_input, 'w') as f: - for line in lines: - f.write(line) - def write_to_output(self, string: str, output_file: str = 'otf_run.out'): """ Write a string or list of strings to the output file. @@ -220,8 +187,11 @@ def write_config(self): string += str('%.6e' % self.structure.stds[i][j]) + ' ' string += '\n' - self.write_to_output(string) + self.write_to_output(string) + # TODO @Simon and Jon perhaps this should return the integer index of an + # atom which has a variance outside of bounds, and it could return + # something like -1 if all the atoms have acceptable force variance? def is_std_in_bound(self): """ Evaluate if the model error are within predefined criteria @@ -231,9 +201,30 @@ def is_std_in_bound(self): # Some decision making criteria + return True + +def run_espresso(qe_input: str, structure: Structure): + """ + Calls quantum espresso from input located at self.qe_input + + :return: List [nparray] List of forces + """ + + edit_qe_input_positions(qe_input, structure) + + pw_loc = os.environ.get('PWSCF_COMMAND', 'pwscf.x') + + qe_command = '{0} < {1} > {2}'.format(pw_loc, qe_input, + 'pwscf.out') + + os.system(qe_command) + + return parse_qe_forces('pwscf.out') + + def parse_qe_input(qe_input: str): """ Reads the positions, species, and cell in from the qe input file @@ -301,6 +292,62 @@ def parse_qe_input(qe_input: str): return positions, species, cell, masses +def edit_qe_input_positions(qe_input: str, structure: Structure): + """ + Write the current configuration of the OTF structure to the + qe input file + """ + + with open(qe_input, 'r') as f: + lines = f.readlines() + + + file_pos_index = None + cell_index = None + nat = None + for i, line in enumerate(lines): + if 'ATOMIC_POSITIONS' in line: + file_pos_index = int(i + 1) + if 'CELL_PARAMETERS' in line: + cell_index = int(i + 1) + if 'nat' in line: + nat = int(line.split('=')[1]) + + assert file_pos_index is not None, 'Failed to find positions in input' + assert cell_index is not None, 'Failed to find cell in input' + assert nat is not None, 'Failed to find nat in input' + + for pos_index, line_index in enumerate( + range(file_pos_index, file_pos_index+structure.nat)): + pos_string = ' '.join([structure.species[pos_index], + str(structure.positions[pos_index][ + 0]), + str(structure.positions[pos_index][ + 1]), + str(structure.positions[pos_index][ + 2])]) + lines[line_index] = str(pos_string + '\n') + + # TODO current assumption: the new structure has fewer atoms than the + # previous one. If we are always 'editing' a copy of the larger structure + # than this will be okay with the punchout method. + for line_index in range(file_pos_index+structure.nat, + file_pos_index + nat): + lines[line_index]='' + + lines[cell_index] = ' '.join([str(x) for x in structure.vec1]) + '\n' + lines[cell_index + 1] = ' '.join([str(x) for x in structure.vec2]) \ + +'\n' + lines[cell_index + 2] = ' '.join([str(x) for x in structure.vec3]) \ + +'\n' + + + with open(qe_input, 'w') as f: + for line in lines: + f.write(line) + + + def parse_qe_forces(outfile: str): """ Get forces from a pwscf file in Ryd/bohr @@ -376,12 +423,11 @@ def parse_otf_output(outfile: str): results['species'] = species print(results) - if __name__ == '__main__': - # os.system('cp ../tests/test_files/qe_input_2.in pwscf.in') + #os.system('cp ../tests/test_files/qe_input_2.in pwscf.in') - # otf = OTF('pwscf.in', .1, 10, kernel='two_body', - # cutoff=10) - # otf.run() - # parse_output('otf_run.out') - pass + otf = OTF('pwscf.in', .1, 10, kernel='two_body', + cutoff=10) + otf.run() + # parse_output('otf_run.out') + pass \ No newline at end of file diff --git a/src/punchout.py b/src/punchout.py index e69de29bb..c22b40a2e 100644 --- a/src/punchout.py +++ b/src/punchout.py @@ -0,0 +1,228 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Code to facilitate the 'punch-out' scheme of large-scale GP-accelerated +molecular dynamics simulations. +""" + +import numpy as np + +from struc import Structure + + +def is_within_d_box(pos1: np.ndarray,pos2: np.ndarray,d)-> bool: + """ + Return if pos2 is within a cube of side length d centered around pos1, + :param pos1: + :param pos2: + :param d: + :return: + """ + + isnear_x = abs(pos1[0] - pos2[0]) <= d/2 + isnear_y = abs(pos1[1] - pos2[1]) <= d/2 + isnear_z = abs(pos1[2] - pos2[2]) <= d/2 + + return isnear_x and isnear_y and isnear_z + + + +def punchout(structure: Structure, atom: int, d: float, center : bool = True): + """ + + :param structure: a Structure type object + """ + a = structure.vec1 + b = structure.vec2 + c = structure.vec3 + + shift_vectors = [ np.zeros(3), + a,b,c, + -a,-b,-c, + a+b,a-b,a+c,a-c, + b+c,b-c, + -a+b,-a-b, + -a+c,-a-c, + -b+c,-b-c, + a+b+c, + -a+b+c, a-b+c, a+b-c, + -a-b+c, -a+b-c, a-b-c, + -a-b-c + ] + + new_pos = [] + new_prev_pos = [] # Necessary so new structure has non-zero velocity + new_species = [] + + assert 0 <= atom <= len(structure.positions), 'Atom index is greater ' \ + 'than number of atoms ' \ + 'in structure' + target_pos = structure.positions[atom] + + + for i, pos in enumerate(structure.positions): + # Compute all shifted positions to handle edge cases + shifted_positions = [pos + shift for shift in shift_vectors] + + for j, shifted_pos in enumerate(shifted_positions): + if is_within_d_box(target_pos, shifted_pos,d): + new_pos.append(shifted_pos) + new_prev_pos.append(structure.prev_positions + [i]+shift_vectors[j]) + new_species.append(structure.species[i]) + + # Set up other new structural properties + newlatt = d*np.eye(3) + new_mass_dict = {} + + for spec in set(new_species): + new_mass_dict[spec] = structure.mass_dict[spec] + + # Instantiate new structure, and set the previous positions manually + + newstruc = Structure(newlatt,new_species,new_pos, + structure.cutoff,new_mass_dict) + newstruc.prev_positions = list(new_prev_pos) + + # Center new structure at origin + if center: + newstruc.translate_positions(-target_pos) + + return newstruc + + + + + +#OLD DRAFT METHOD: DEPRECATED + +def simple_cube_local_cell(atoms): + """ + Turns a structure in to a simple cubic cell with 1 angstrom of 'padding' on all sides. + Used for 'punch-outs'. + + :param structure: + :return: + """ + + positions = [at.position for at in atoms] + + minx = min([pos[0] for pos in positions]) + miny = min([pos[1] for pos in positions]) + minz = min([pos[2] for pos in positions]) + + maxx = max([pos[0] for pos in positions]) + maxy = max([pos[1] for pos in positions]) + maxz = max([pos[2] for pos in positions]) + + x_extent = maxx-minx + y_extent = maxy-miny + z_extent = maxz-minz + + x_distances=[] + y_distances=[] + z_distances=[] + + for i,pos1 in enumerate(positions): + for j,pos2 in enumerate(positions): + if pos1[0]-pos2[0]!=0: + x_distances.append(np.abs(pos1[0]-pos2[0])/2.) + if pos1[1]-pos2[1]!=0: + y_distances.append(np.abs(pos1[1]-pos2[1])/2.) + if pos1[2]-pos2[2]!=0: + z_distances.append(np.abs(pos1[2]-pos2[2])/2.) + + + #### + # NOT RELEVANT RN: Give each atom a buffer of 1 angstrom on the side of the cubes + #for at in atoms: + # at.position+=np.array([1,1,1]) + + xijbar= np.mean(x_distances) + yijbar= np.mean(y_distances) + zijbar= np.mean(z_distances) + + lattice = np.diag([x_extent+xijbar,y_extent+yijbar,z_extent+zijbar]) + + return Structure(atoms,alat=1,lattice=lattice) + + +if __name__=='__main__': + + + + #from ase.build.bulk import bulk + #from ase.build.supercells import make_supercell + + #ats = bulk('Si','fcc',a=5.431) + + #Sats = make_supercell(ats,7*np.eye(3)) + + #struc=ase_to_structure(Sats,alat=5.14,fractional=False,perturb=.5) + + #print(struc) + #struc = punchout_cell(struc,target_position=(5.0,5.0,5.0),box_dist=5.0) + + #print(struc) + """spec ='H' + + positions=[] + + for i in range(5): + for j in range(5): + for k in range(5): + positions.append((i,j,k)) + print(positions) + print(len(positions)) + + alat =1.0 + lattice = 5*np.eye(3) + + atoms = [Atom(position=pos,element='H') for pos in positions] + + + struc= Structure(atoms,alat=1,lattice=5*np.eye(3),fractional=False) + print(struc[2]) + b = get_box_neighbors(struc,target_atom=struc[2],box_dist=.5) + for x in b: + print(x) + + + c = generate_ultracell(struc) + for x in c: + pass + #print(x) + print(len(struc)) + print(len(c)) + + ultra_pos = [tuple(x.position) for x in c] + print(len(ultra_pos)) + print(len(set(ultra_pos))) + + print(get_box_neighbors(struc,target_pos=[0,0,0],box_dist=1.0)) + print('') + nbr_atoms = get_box_neighbors(c,target_pos=[0,0,0],box_dist=1.0) + + print(isinstance(struc,Structure)) + print(simple_cube_local_cell(center_atoms(nbr_atoms))) + """ + + + + + + + + + + + + + + + + + + + + From 8eaf64eab53c0eb56c75f876d5b16992a0d7316f Mon Sep 17 00:00:00 2001 From: Steven T Date: Sat, 8 Sep 2018 01:11:26 -0400 Subject: [PATCH 3/9] Added new method to struc to translate entire structure (useful for punchout) --- src/struc.py | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/src/struc.py b/src/struc.py index 10827c9f1..b116e1425 100644 --- a/src/struc.py +++ b/src/struc.py @@ -1,6 +1,6 @@ -from numpy import zeros +from numpy import zeros,ndarray from numpy.linalg import inv - +from typing import List class Structure(object): """ @@ -14,7 +14,8 @@ class Structure(object): """ - def __init__(self, lattice, species, positions, cutoff, mass_dict=None): + def __init__(self, lattice: ndarray, species:List[str], + positions: List[ndarray],cutoff: float,mass_dict: dict = None): self.lattice = lattice self.vec1 = lattice[0, :] self.vec2 = lattice[1, :] @@ -44,6 +45,19 @@ def __init__(self, lattice, species, positions, cutoff, mass_dict=None): self.mass_dict = mass_dict + def translate_positions(self, vector : ndarray=zeros(3)): + """ + Translate all positions, and previous positions by vector + :param vector: vector to translate by + :return: + """ + for pos in self.positions: + pos += vector + + for prev_pos in self.prev_positions: + prev_pos += vector + + @staticmethod def calc_bond_list(unique_species): """Converts unique species to a list of bonds. @@ -75,6 +89,10 @@ def calc_triplet_list(unique_species): triplet_list.append([m, n, p]) return triplet_list + + + + if __name__ == '__main__': # create simple structure pass From d810c068415571b47c2e44b9d13cf1561f0a7ac5 Mon Sep 17 00:00:00 2001 From: Steven T Date: Sat, 8 Sep 2018 01:33:16 -0400 Subject: [PATCH 4/9] Punchout added to OTF class/flow --- src/otf.py | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/src/otf.py b/src/otf.py index acf833aea..401403088 100644 --- a/src/otf.py +++ b/src/otf.py @@ -103,6 +103,7 @@ def run_and_train(self): forces = run_espresso(self.qe_input,self.train_structure) + self.write_to_output('Done.\n') # Write input positions and force results @@ -213,15 +214,20 @@ def run_espresso(qe_input: str, structure: Structure): :return: List [nparray] List of forces """ - edit_qe_input_positions(qe_input, structure) + run_qe_path = qe_input+'_run' + + os.system(' '.join(['cp',qe_input,run_qe_path])) + edit_qe_input_positions(run_qe_path, structure) pw_loc = os.environ.get('PWSCF_COMMAND', 'pwscf.x') - qe_command = '{0} < {1} > {2}'.format(pw_loc, qe_input, + qe_command = '{0} < {1} > {2}'.format(pw_loc, run_qe_path, 'pwscf.out') os.system(qe_command) + #os.system(' '.join(['rm', run_qe_path])) + return parse_qe_forces('pwscf.out') @@ -310,8 +316,11 @@ def edit_qe_input_positions(qe_input: str, structure: Structure): file_pos_index = int(i + 1) if 'CELL_PARAMETERS' in line: cell_index = int(i + 1) + # Load nat into variable then overwrite it with new nat if 'nat' in line: nat = int(line.split('=')[1]) + nat_index = int(i) + lines[nat_index]='nat = '+str(structure.nat)+'\n' assert file_pos_index is not None, 'Failed to find positions in input' assert cell_index is not None, 'Failed to find cell in input' @@ -328,9 +337,10 @@ def edit_qe_input_positions(qe_input: str, structure: Structure): 2])]) lines[line_index] = str(pos_string + '\n') - # TODO current assumption: the new structure has fewer atoms than the - # previous one. If we are always 'editing' a copy of the larger structure - # than this will be okay with the punchout method. + # TODO current assumption: if there is a new structure, then the new + # structure has fewer atoms than the previous one. If we are always + # 'editing' a version of the larger structure than this will be okay with + # the punchout method. for line_index in range(file_pos_index+structure.nat, file_pos_index + nat): lines[line_index]='' @@ -427,7 +437,7 @@ def parse_otf_output(outfile: str): #os.system('cp ../tests/test_files/qe_input_2.in pwscf.in') otf = OTF('pwscf.in', .1, 10, kernel='two_body', - cutoff=10) + cutoff=10,punchout_d=None) otf.run() # parse_output('otf_run.out') pass \ No newline at end of file From 1450ed01025f02a506afae23caddf2982f179421 Mon Sep 17 00:00:00 2001 From: Steven T Date: Mon, 10 Sep 2018 14:45:53 -0400 Subject: [PATCH 5/9] Added qe_util, rcut.py, and unit tests for rcut helper methods --- src/qe_util.py | 0 src/rcut.py | 72 ++++++++++++++++++++++++++++++++++++++++++++++ tests/test_rcut.py | 0 3 files changed, 72 insertions(+) create mode 100644 src/qe_util.py create mode 100644 src/rcut.py create mode 100644 tests/test_rcut.py diff --git a/src/qe_util.py b/src/qe_util.py new file mode 100644 index 000000000..e69de29bb diff --git a/src/rcut.py b/src/rcut.py new file mode 100644 index 000000000..fbabd63f5 --- /dev/null +++ b/src/rcut.py @@ -0,0 +1,72 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Self-Contained Code which allows for the +""" + +import numpy as np +import numpy.random as rand + +from math import sin, cos +from struc import Structure +from qe_util import run_espresso + +def perturb_position(pos: np.array, r: float=.1, rscale = .02,): + """ + + :param pos: + :param r: + :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,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 = .02): + + + 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 np.linalg.norm(pos-central_pos) > r_fix: + new_positions.append(perturb_position(pos,mean_pert,pert_sigma)) + else: + new_positions.append(pos) + + newstruc = Structure(lattice=structure.lattice, species=structure.species, + positions=new_positions,cutoff=structure.cutoff) + + return newstruc + + + + + + + + + + + + + + diff --git a/tests/test_rcut.py b/tests/test_rcut.py new file mode 100644 index 000000000..e69de29bb From d47e5e0388c5ee8a94bf8725a060c2ecaad8acab Mon Sep 17 00:00:00 2001 From: Steven T Date: Mon, 10 Sep 2018 14:46:11 -0400 Subject: [PATCH 6/9] Refactored otf for qe_util integration --- src/otf.py | 251 +++++++++-------------------------------------------- 1 file changed, 40 insertions(+), 211 deletions(-) diff --git a/src/otf.py b/src/otf.py index 401403088..819953e09 100644 --- a/src/otf.py +++ b/src/otf.py @@ -8,17 +8,18 @@ """ import numpy as np -import os import datetime from struc import Structure from math import sqrt from gp import GaussianProcess from env import ChemicalEnvironment from punchout import punchout +from qe_util import run_espresso, parse_qe_input + class OTF(object): def __init__(self, qe_input: str, dt: float, number_of_steps: int, - kernel: str, cutoff: float, punchout_d: float=None): + kernel: str, cutoff: float, punchout_d: float = None): """ On-the-fly learning engine, containing methods to run OTF calculation @@ -94,15 +95,16 @@ def run_and_train(self): # Run espresso and write out results self.write_to_output('=' * 20 + '\n' + 'Calling QE... ') - if self.punchout_d == None: + # If not in punchout mode, run QE on the entire structure + if self.punchout_d is None: self.train_structure = self.structure + # If first run, pick a random atom to punch out a structure around elif self.train_structure is None: - target_atom = np.random.randint(0,len(self.structure.positions)) - self.train_structure = punchout(self.structure,target_atom, + target_atom = np.random.randint(0, len(self.structure.positions)) + self.train_structure = punchout(self.structure, target_atom, d=self.punchout_d) - - forces = run_espresso(self.qe_input,self.train_structure) + forces = run_espresso(self.qe_input, self.train_structure) self.write_to_output('Done.\n') @@ -136,19 +138,22 @@ def update_positions(self): Apply a timestep to self.structure based on current structure's forces. """ - # Maintain list of elemental masses in amu to calculate acceleration + # Precompute dt squared for efficiency + dtdt = self.dt ** 2 + + for i, pre_pos in enumerate(self.structure.prev_positions): + temp_pos = np.array(self.structure.positions[i]) + mass = self.structure.mass_dict[self.structure.species[i]] + pos = self.structure.positions[i] + forces = self.structure.forces[i] - for i, prev_pos in enumerate(self.structure.prev_positions): - temp_pos = self.structure.positions[i] - self.structure.positions[i] = 2 * self.structure.positions[i] - \ - prev_pos + self.dt ** 2 * \ - self.structure.forces[i] \ - / self.structure.mass_dict[ - self.structure.species[i]] + self.structure.positions[i] = 2 * pos - pre_pos + dtdt * forces / \ + mass self.structure.prev_positions[i] = np.array(temp_pos) - def write_to_output(self, string: str, output_file: str = 'otf_run.out'): + @staticmethod + def write_to_output(string: str, output_file: str = 'otf_run.out'): """ Write a string or list of strings to the output file. :param string: String to print to output @@ -190,9 +195,7 @@ def write_config(self): self.write_to_output(string) - # TODO @Simon and Jon perhaps this should return the integer index of an - # atom which has a variance outside of bounds, and it could return - # something like -1 if all the atoms have acceptable force variance? + # TODO change this to use the signal variance def is_std_in_bound(self): """ Evaluate if the model error are within predefined criteria @@ -202,196 +205,21 @@ def is_std_in_bound(self): # Some decision making criteria + # SKETCH OF EVENTUAL CODE: + """ + high_error_atom = -1 + + if self.punchout_d is not None: + self.train_structure= punchout(self.structure,atom= + high_error_atom, d= self.punchout_d ) + return False + """ return True - -def run_espresso(qe_input: str, structure: Structure): - """ - Calls quantum espresso from input located at self.qe_input - - :return: List [nparray] List of forces - """ - - run_qe_path = qe_input+'_run' - - os.system(' '.join(['cp',qe_input,run_qe_path])) - edit_qe_input_positions(run_qe_path, structure) - - pw_loc = os.environ.get('PWSCF_COMMAND', 'pwscf.x') - - qe_command = '{0} < {1} > {2}'.format(pw_loc, run_qe_path, - 'pwscf.out') - - os.system(qe_command) - - #os.system(' '.join(['rm', run_qe_path])) - - return parse_qe_forces('pwscf.out') - - -def parse_qe_input(qe_input: str): - """ - Reads the positions, species, and cell in from the qe input file - - :param qe_input: str, Path to PWSCF input file - :return: List[nparray], List[str], nparray, Positions, species, - 3x3 Bravais cell - """ - positions = [] - species = [] - cell = [] - - with open(qe_input) as f: - lines = f.readlines() - - # Find the cell and positions in the output file - cell_index = None - positions_index = None - nat = None - species_index = None - - for i, line in enumerate(lines): - if 'CELL_PARAMETERS' in line: - cell_index = int(i + 1) - if 'ATOMIC_POSITIONS' in line: - positions_index = int(i + 1) - if 'nat' in line: - nat = int(line.split('=')[1]) - if 'ATOMIC_SPECIES' in line: - species_index = int(i + 1) - - assert cell_index is not None, 'Failed to find cell in input' - assert positions_index is not None, 'Failed to find positions in input' - assert nat is not None, 'Failed to find number of atoms in input' - assert species_index is not None, 'Failed to find atomic species in input' - - # Load cell - for i in range(cell_index, cell_index + 3): - cell_line = lines[i].strip() - cell.append(np.fromstring(cell_line, sep=' ')) - cell = np.array(cell) - - # Check cell IO - assert cell != [], 'Cell failed to load' - assert np.shape(cell) == (3, 3), 'Cell failed to load correctly' - - # Load positions - for i in range(positions_index, positions_index + nat): - line_string = lines[i].strip().split() - species.append(line_string[0]) - - pos_string = ' '.join(line_string[1:4]) - - positions.append(np.fromstring(pos_string, sep=' ')) - # Check position IO - assert positions != [], "Positions failed to load" - - # Load masses - masses = {} - for i in range(species_index, species_index + len(set(species))): - # Expects lines of format like: H 1.0 H_pseudo_name.ext - line = lines[i].strip().split() - masses[line[0]] = float(line[1]) - - return positions, species, cell, masses - - -def edit_qe_input_positions(qe_input: str, structure: Structure): - """ - Write the current configuration of the OTF structure to the - qe input file - """ - - with open(qe_input, 'r') as f: - lines = f.readlines() - - - file_pos_index = None - cell_index = None - nat = None - for i, line in enumerate(lines): - if 'ATOMIC_POSITIONS' in line: - file_pos_index = int(i + 1) - if 'CELL_PARAMETERS' in line: - cell_index = int(i + 1) - # Load nat into variable then overwrite it with new nat - if 'nat' in line: - nat = int(line.split('=')[1]) - nat_index = int(i) - lines[nat_index]='nat = '+str(structure.nat)+'\n' - - assert file_pos_index is not None, 'Failed to find positions in input' - assert cell_index is not None, 'Failed to find cell in input' - assert nat is not None, 'Failed to find nat in input' - - for pos_index, line_index in enumerate( - range(file_pos_index, file_pos_index+structure.nat)): - pos_string = ' '.join([structure.species[pos_index], - str(structure.positions[pos_index][ - 0]), - str(structure.positions[pos_index][ - 1]), - str(structure.positions[pos_index][ - 2])]) - lines[line_index] = str(pos_string + '\n') - - # TODO current assumption: if there is a new structure, then the new - # structure has fewer atoms than the previous one. If we are always - # 'editing' a version of the larger structure than this will be okay with - # the punchout method. - for line_index in range(file_pos_index+structure.nat, - file_pos_index + nat): - lines[line_index]='' - - lines[cell_index] = ' '.join([str(x) for x in structure.vec1]) + '\n' - lines[cell_index + 1] = ' '.join([str(x) for x in structure.vec2]) \ - +'\n' - lines[cell_index + 2] = ' '.join([str(x) for x in structure.vec3]) \ - +'\n' - - - with open(qe_input, 'w') as f: - for line in lines: - f.write(line) - - - -def parse_qe_forces(outfile: str): - """ - Get forces from a pwscf file in Ryd/bohr - - :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]) * 13.605698066 - - 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) - - return forces - - -# TODO Currently won't work: needs to be re-done in light of new output -# formatting, but could wait until we finalize an output file format , -# will be useful for data analysis +# TODO Currently won't work: needs to be re-done when we finalize our output +# formatting def parse_otf_output(outfile: str): """ Parse the output of a otf run for analysis @@ -433,11 +261,12 @@ def parse_otf_output(outfile: str): results['species'] = species print(results) + if __name__ == '__main__': - #os.system('cp ../tests/test_files/qe_input_2.in pwscf.in') + # os.system('cp ../tests/test_files/qe_input_2.in pwscf.in') - otf = OTF('pwscf.in', .1, 10, kernel='two_body', - cutoff=10,punchout_d=None) - otf.run() - # parse_output('otf_run.out') - pass \ No newline at end of file + otf = OTF('pwscf.in', .1, 10, kernel='two_body', + cutoff=10, punchout_d=None) + otf.run() + # parse_output('otf_run.out') + pass From 552446c6b9fef096080efbf0618ee8912fcd1256 Mon Sep 17 00:00:00 2001 From: Steven T Date: Mon, 10 Sep 2018 14:46:25 -0400 Subject: [PATCH 7/9] Minor changes to punchout --- src/punchout.py | 232 ++++++++++++++---------------------------------- 1 file changed, 65 insertions(+), 167 deletions(-) diff --git a/src/punchout.py b/src/punchout.py index c22b40a2e..6938d63a8 100644 --- a/src/punchout.py +++ b/src/punchout.py @@ -3,6 +3,8 @@ """ Code to facilitate the 'punch-out' scheme of large-scale GP-accelerated molecular dynamics simulations. + +Steven Torrisi """ import numpy as np @@ -10,69 +12,80 @@ from struc import Structure -def is_within_d_box(pos1: np.ndarray,pos2: np.ndarray,d)-> bool: +def is_within_d_box(pos1: np.array, pos2: np.array, d: float) -> bool: """ Return if pos2 is within a cube of side length d centered around pos1, - :param pos1: - :param pos2: - :param d: - :return: + :param pos1: First position + :type pos1: np.array + :param pos2: Second position + :type pos2: np.array + :param d: Side length of cube + :type d: float + :return: bool """ - isnear_x = abs(pos1[0] - pos2[0]) <= d/2 - isnear_y = abs(pos1[1] - pos2[1]) <= d/2 - isnear_z = abs(pos1[2] - pos2[2]) <= d/2 + isnear_x = abs(pos1[0] - pos2[0]) <= d / 2 + isnear_y = abs(pos1[1] - pos2[1]) <= d / 2 + isnear_z = abs(pos1[2] - pos2[2]) <= d / 2 return isnear_x and isnear_y and isnear_z - -def punchout(structure: Structure, atom: int, d: float, center : bool = True): +def punchout(structure: Structure, atom: int, d: float, center: bool = True): """ - - :param structure: a Structure type object + Punch out a cube of side-length d around the target atom + + :param structure: Structure to punch out a sub-structure from + :type structure: Structure + :param atom: Index of atom to punch out around + :type atom: int + :param d: Side length of cube around target atom + :type d: float + :param center: Return structure centered around 0 or not + :type center: bool """ + + # TODO replace this with the env's method for finding neighbors a = structure.vec1 b = structure.vec2 c = structure.vec3 - shift_vectors = [ np.zeros(3), - a,b,c, - -a,-b,-c, - a+b,a-b,a+c,a-c, - b+c,b-c, - -a+b,-a-b, - -a+c,-a-c, - -b+c,-b-c, - a+b+c, - -a+b+c, a-b+c, a+b-c, - -a-b+c, -a+b-c, a-b-c, - -a-b-c - ] + shift_vectors = [np.zeros(3), + a, b, c, + -a, -b, -c, + a + b, a - b, a + c, a - c, + b + c, b - c, + -a + b, -a - b, + -a + c, -a - c, + -b + c, -b - c, + a + b + c, + -a + b + c, a - b + c, a + b - c, + -a - b + c, -a + b - c, a - b - c, + -a - b - c + ] new_pos = [] - new_prev_pos = [] # Necessary so new structure has non-zero velocity + new_prev_pos = [] # Necessary so new structure has non-zero velocity new_species = [] assert 0 <= atom <= len(structure.positions), 'Atom index is greater ' \ - 'than number of atoms ' \ + 'than number of atoms ' \ 'in structure' target_pos = structure.positions[atom] - for i, pos in enumerate(structure.positions): # Compute all shifted positions to handle edge cases shifted_positions = [pos + shift for shift in shift_vectors] for j, shifted_pos in enumerate(shifted_positions): - if is_within_d_box(target_pos, shifted_pos,d): + if is_within_d_box(target_pos, shifted_pos, d): new_pos.append(shifted_pos) new_prev_pos.append(structure.prev_positions - [i]+shift_vectors[j]) + [i] + shift_vectors[j]) new_species.append(structure.species[i]) # Set up other new structural properties - newlatt = d*np.eye(3) + newlatt = d * np.eye(3) new_mass_dict = {} for spec in set(new_species): @@ -80,149 +93,34 @@ def punchout(structure: Structure, atom: int, d: float, center : bool = True): # Instantiate new structure, and set the previous positions manually - newstruc = Structure(newlatt,new_species,new_pos, - structure.cutoff,new_mass_dict) + newstruc = Structure(newlatt, new_species, new_pos, + structure.cutoff, new_mass_dict) newstruc.prev_positions = list(new_prev_pos) - # Center new structure at origin - if center: - newstruc.translate_positions(-target_pos) - - return newstruc - - - - - -#OLD DRAFT METHOD: DEPRECATED - -def simple_cube_local_cell(atoms): - """ - Turns a structure in to a simple cubic cell with 1 angstrom of 'padding' on all sides. - Used for 'punch-outs'. - - :param structure: - :return: - """ - - positions = [at.position for at in atoms] - - minx = min([pos[0] for pos in positions]) - miny = min([pos[1] for pos in positions]) - minz = min([pos[2] for pos in positions]) - - maxx = max([pos[0] for pos in positions]) - maxy = max([pos[1] for pos in positions]) - maxz = max([pos[2] for pos in positions]) - - x_extent = maxx-minx - y_extent = maxy-miny - z_extent = maxz-minz - - x_distances=[] - y_distances=[] - z_distances=[] - - for i,pos1 in enumerate(positions): - for j,pos2 in enumerate(positions): - if pos1[0]-pos2[0]!=0: - x_distances.append(np.abs(pos1[0]-pos2[0])/2.) - if pos1[1]-pos2[1]!=0: - y_distances.append(np.abs(pos1[1]-pos2[1])/2.) - if pos1[2]-pos2[2]!=0: - z_distances.append(np.abs(pos1[2]-pos2[2])/2.) - - - #### - # NOT RELEVANT RN: Give each atom a buffer of 1 angstrom on the side of the cubes - #for at in atoms: - # at.position+=np.array([1,1,1]) - - xijbar= np.mean(x_distances) - yijbar= np.mean(y_distances) - zijbar= np.mean(z_distances) - - lattice = np.diag([x_extent+xijbar,y_extent+yijbar,z_extent+zijbar]) - - return Structure(atoms,alat=1,lattice=lattice) - - -if __name__=='__main__': - - - - #from ase.build.bulk import bulk - #from ase.build.supercells import make_supercell - - #ats = bulk('Si','fcc',a=5.431) - - #Sats = make_supercell(ats,7*np.eye(3)) - - #struc=ase_to_structure(Sats,alat=5.14,fractional=False,perturb=.5) - - #print(struc) - #struc = punchout_cell(struc,target_position=(5.0,5.0,5.0),box_dist=5.0) - - #print(struc) - """spec ='H' - - positions=[] - - for i in range(5): - for j in range(5): - for k in range(5): - positions.append((i,j,k)) - print(positions) - print(len(positions)) - - alat =1.0 - lattice = 5*np.eye(3) - - atoms = [Atom(position=pos,element='H') for pos in positions] - - - struc= Structure(atoms,alat=1,lattice=5*np.eye(3),fractional=False) - print(struc[2]) - b = get_box_neighbors(struc,target_atom=struc[2],box_dist=.5) - for x in b: - print(x) - - - c = generate_ultracell(struc) - for x in c: - pass - #print(x) - print(len(struc)) - print(len(c)) - - ultra_pos = [tuple(x.position) for x in c] - print(len(ultra_pos)) - print(len(set(ultra_pos))) - - print(get_box_neighbors(struc,target_pos=[0,0,0],box_dist=1.0)) - print('') - nbr_atoms = get_box_neighbors(c,target_pos=[0,0,0],box_dist=1.0) - - print(isinstance(struc,Structure)) - print(simple_cube_local_cell(center_atoms(nbr_atoms))) - """ - - - - - - - - - - - - + # Check if any atoms are unphysically close to each other, add padding + # accordingly + # TODO put in the new method here, except applied on the newstructure + # Sketch of how it will work... + # BELOW CODE CURRENTLY DOES NOTHING; awaiting update to env.py + for i, pos in enumerate(newstruc.positions): + # Compute all shifted positions to handle edge cases + shifted_positions = [pos + shift for shift in shift_vectors if + not np.equal(shift, np.zeros(3)).all()] + for j, shifted_pos in enumerate(shifted_positions): + if is_within_d_box(pos, shifted_pos, 1.0): + # shift box + pass + # Center new structure at origin + if center: + newstruc.translate_positions(-target_pos) + return newstruc +if __name__ == '__main__': + pass From 7ca598776c90f994f978c34156f971f41c96fa1b Mon Sep 17 00:00:00 2001 From: Steven T Date: Mon, 10 Sep 2018 14:46:34 -0400 Subject: [PATCH 8/9] QE Util better documented --- src/qe_util.py | 203 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 203 insertions(+) diff --git a/src/qe_util.py b/src/qe_util.py index e69de29bb..45765d869 100644 --- a/src/qe_util.py +++ b/src/qe_util.py @@ -0,0 +1,203 @@ +#!/usr/bin/env python3 +# pylint: disable=redefined-outer-name +"""" +Assembly of helper functions which call pw.x in the Quantum ESPRESSO suite + +Assumes the environment variable PWSCF_COMMAND is set as a path to the pw.x +binary + +Steven Torrisi +""" + +import os +import numpy as np +from struc import Structure +from typing import List + + +def run_espresso(qe_input: str, structure: Structure, temp: bool = True) -> \ + List[np.array]: + """ + Calls quantum espresso from input located at self.qe_input + + :param qe_input: Path to pwscf.in file + :type qe_input: str + :param structure: Structure object for the edited + :type structure: Structure + :param temp: Run pwscf off of an edited QE input instead of the original + :type temp: Bool + + :return: List [nparray] List of forces + """ + + # If temp, copy the extant file into a new run directory + run_qe_path = qe_input + run_qe_path += '_run ' if temp else '' + + os.system(' '.join(['cp', qe_input, run_qe_path])) + edit_qe_input_positions(run_qe_path, structure) + + pw_loc = os.environ.get('PWSCF_COMMAND', 'pwscf.x') + + qe_command = '{0} < {1} > {2}'.format(pw_loc, run_qe_path, + 'pwscf.out') + + os.system(qe_command) + + return parse_qe_forces('pwscf.out') + + +def parse_qe_input(qe_input: str): + """ + Reads the positions, species, cell, and masses in from the qe input file + + :param qe_input: Path to PWSCF input file + :type qe_input: str + :return: List[nparray], List[str], nparray, Positions, species, + 3x3 Bravais cell + """ + positions = [] + species = [] + cell = [] + + with open(qe_input) as f: + lines = f.readlines() + + # Find the cell and positions in the output file + cell_index = None + positions_index = None + nat = None + species_index = None + + for i, line in enumerate(lines): + if 'CELL_PARAMETERS' in line: + cell_index = int(i + 1) + if 'ATOMIC_POSITIONS' in line: + positions_index = int(i + 1) + if 'nat' in line: + nat = int(line.split('=')[1]) + if 'ATOMIC_SPECIES' in line: + species_index = int(i + 1) + + assert cell_index is not None, 'Failed to find cell in input' + assert positions_index is not None, 'Failed to find positions in input' + assert nat is not None, 'Failed to find number of atoms in input' + assert species_index is not None, 'Failed to find atomic species in input' + + # Load cell + for i in range(cell_index, cell_index + 3): + cell_line = lines[i].strip() + cell.append(np.fromstring(cell_line, sep=' ')) + cell = np.array(cell) + + # Check cell IO + assert cell != [], 'Cell failed to load' + assert np.shape(cell) == (3, 3), 'Cell failed to load correctly' + + # Load positions + for i in range(positions_index, positions_index + nat): + line_string = lines[i].strip().split() + species.append(line_string[0]) + + pos_string = ' '.join(line_string[1:4]) + + positions.append(np.fromstring(pos_string, sep=' ')) + # Check position IO + assert positions != [], "Positions failed to load" + + # Load masses + masses = {} + for i in range(species_index, species_index + len(set(species))): + # Expects lines of format like: H 1.0 H_pseudo_name.ext + line = lines[i].strip().split() + masses[line[0]] = float(line[1]) + + return positions, species, cell, masses + + +def edit_qe_input_positions(qe_input: str, structure: Structure): + """ + Write the current configuration of the OTF structure to the + qe input file + """ + + with open(qe_input, 'r') as f: + lines = f.readlines() + + file_pos_index = None + cell_index = None + nat = None + for i, line in enumerate(lines): + if 'ATOMIC_POSITIONS' in line: + file_pos_index = int(i + 1) + if 'CELL_PARAMETERS' in line: + cell_index = int(i + 1) + # Load nat into variable then overwrite it with new nat + if 'nat' in line: + nat = int(line.split('=')[1]) + nat_index = int(i) + lines[nat_index] = 'nat = ' + str(structure.nat) + '\n' + + assert file_pos_index is not None, 'Failed to find positions in input' + assert cell_index is not None, 'Failed to find cell in input' + assert nat is not None, 'Failed to find nat in input' + + for pos_index, line_index in enumerate( + range(file_pos_index, file_pos_index + structure.nat)): + pos_string = ' '.join([structure.species[pos_index], + str(structure.positions[pos_index][ + 0]), + str(structure.positions[pos_index][ + 1]), + str(structure.positions[pos_index][ + 2])]) + lines[line_index] = str(pos_string + '\n') + + # TODO current assumption: if there is a new structure, then the new + # structure has fewer atoms than the previous one. If we are always + # 'editing' a version of the larger structure than this will be okay with + # the punchout method. + for line_index in range(file_pos_index + structure.nat, + file_pos_index + nat): + lines[line_index] = '' + + lines[cell_index] = ' '.join([str(x) for x in structure.vec1]) + '\n' + lines[cell_index + 1] = ' '.join([str(x) for x in structure.vec2]) \ + + '\n' + lines[cell_index + 2] = ' '.join([str(x) for x in structure.vec3]) \ + + '\n' + + with open(qe_input, 'w') as f: + for line in lines: + f.write(line) + + +def parse_qe_forces(outfile: str): + """ + Get forces from a pwscf file in Ryd/bohr + + :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]) * 13.605698066 + + 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) + + return forces From 8c4aac6f102f2ff21faa82bbb167d6839b8b4726 Mon Sep 17 00:00:00 2001 From: Steven T Date: Mon, 10 Sep 2018 14:47:04 -0400 Subject: [PATCH 9/9] Added more unit tests and modified a few --- tests/test_GaussianProcess.py | 2 +- tests/test_OTF.py | 72 ++++++++++++++++-------------- tests/test_files/qe_input_1.in | 2 +- tests/test_punchout.py | 1 - tests/test_rcut.py | 81 ++++++++++++++++++++++++++++++++++ 5 files changed, 122 insertions(+), 36 deletions(-) diff --git a/tests/test_GaussianProcess.py b/tests/test_GaussianProcess.py index a9d931f9e..6e52619e2 100644 --- a/tests/test_GaussianProcess.py +++ b/tests/test_GaussianProcess.py @@ -27,7 +27,7 @@ def get_random_structure(cell, unique_species, cutoff, noa): for n in range(noa): positions.append(np.random.uniform(-1, 1, 3)) forces.append(np.random.uniform(-1, 1, 3)) - species.append(unique_species[np.random.randint(0, 2)]) + species.append(unique_species[np.random.randint(0, len(unique_species))]) test_structure = Structure(cell, species, positions, cutoff) diff --git a/tests/test_OTF.py b/tests/test_OTF.py index 38ef7097d..2d3a9f9c9 100644 --- a/tests/test_OTF.py +++ b/tests/test_OTF.py @@ -14,11 +14,14 @@ import os import sys import numpy as np + sys.path.append("../src") -from otf import parse_qe_input, parse_qe_forces, OTF, edit_qe_input_positions,\ - run_espresso +from otf import OTF from gp import GaussianProcess from struc import Structure +from qe_util import parse_qe_forces, parse_qe_input, \ + edit_qe_input_positions, run_espresso + def fake_espresso(qe_input: str, structure: Structure): """ @@ -26,17 +29,19 @@ def fake_espresso(qe_input: str, structure: Structure): in order to allow for substitution. """ noa = len(structure.positions) - return [np.random.normal(loc=0,scale=1,size=3) for _ in range(noa)] + return [np.random.normal(loc=0, scale=1, size=3) for _ in range(noa)] + def fake_predict_on_structure(self, structure): - """ + """ Substitutes in the predict_on_structure method of GaussianProcess """ - structure.forces = [np.random.randn(3) for n in structure.positions] - structure.stds = [np.random.randn(3) for n in structure.positions] + structure.forces = [np.random.randn(3) for n in structure.positions] + structure.stds = [np.random.randn(3) for n in structure.positions] + + return structure - return structure class Fake_GP(GaussianProcess): """ @@ -68,13 +73,12 @@ def predict(self, chemenv, i): """ Substitutes in the predict method of GaussianProcess """ - #structure.forces = [np.random.randn(3) for n in range(structure.nat)] - #structure.stds = [np.random.randn(3) for n in range(structure.nat)] + # structure.forces = [np.random.randn(3) for n in range(structure.nat)] + # structure.stds = [np.random.randn(3) for n in range(structure.nat)] return chemenv - def cleanup_otf_run(): os.system('rm otf_run.out') os.system('rm pwscf.out') @@ -82,6 +86,7 @@ def cleanup_otf_run(): os.system('rm pwscf.save') os.system('rm otf_run.out') + # ------------------------------------------------------ # fixtures # ------------------------------------------------------ @@ -152,18 +157,20 @@ def test_cell_parsing(qe_input, exp_cell): positions, species, cell, masses = parse_qe_input(qe_input) assert np.all(exp_cell == cell) + @pytest.mark.parametrize("qe_input,mass_dict", [ ('./test_files/qe_input_1.in', - {'H':1.0}), + {'H': 1.0}), ('./test_files/qe_input_2.in', - {'Al':26.9815385}) + {'Al': 26.9815385}) ] ) def test_cell_parsing(qe_input, mass_dict): positions, species, cell, masses = parse_qe_input(qe_input) assert masses == mass_dict + @pytest.mark.parametrize('qe_output,exp_forces', [ ('./test_files/qe_output_1.out', @@ -185,7 +192,7 @@ def test_espresso_calling_1(test_otf_engine_1): 'in environment' forces = run_espresso(test_otf_engine_1.qe_input, - test_otf_engine_1.structure) + test_otf_engine_1.structure, temp=False) assert isinstance(forces, list) assert len(forces) == len(test_otf_engine_1.structure.forces) @@ -198,6 +205,7 @@ def test_espresso_calling_1(test_otf_engine_1): cleanup_otf_run() + def test_espresso_input_edit(): """ Load a structure in from qe_input_1, change the position and cell, @@ -205,46 +213,44 @@ def test_espresso_input_edit(): :return: """ os.system('cp test_files/qe_input_1.in .') - positions, species, cell, masses = parse_qe_input('./qe_input_1.in') - struc = Structure(cell,species,positions,masses) + positions, species, cell, masses = parse_qe_input('./qe_input_1.in') + struc = Structure(cell, species, positions, masses) struc.vec1 += np.random.randn(3) - struc.positions[0]+= np.random.randn(3) + struc.positions[0] += np.random.randn(3) edit_qe_input_positions('./qe_input_1.in', structure=struc) - positions, species, cell, masses = parse_qe_input('./qe_input_1.in') + positions, species, cell, masses = parse_qe_input('./qe_input_1.in') - assert np.equal(positions[0],struc.positions[0]).all() - assert np.equal(struc.vec1,cell[0, :]).all() + assert np.equal(positions[0], struc.positions[0]).all() + assert np.equal(struc.vec1, cell[0, :]).all() os.system('rm ./qe_input_1.in') - # ------------------------------------------------------ # test otf methods # ------------------------------------------------------ - -#TODO see if there is a better way to to set up the different input runs +# TODO see if there is a better way to to set up the different input runs def test_update_1(test_otf_engine_1): + test_otf_engine_1.structure.prev_positions = [[2.5, 2.5, 2.5], + [4.5, 2.5, 2.5]] - test_otf_engine_1.structure.prev_positions=[[2.5,2.5,2.5], - [4.5,2.5,2.5]] - - test_otf_engine_1.structure.forces=[np.array([.07413986, 0.0, 0.0]), - np.array([-0.07413986, 0.0, 0.0])] + test_otf_engine_1.structure.forces = [np.array([.07413986, 0.0, 0.0]), + np.array([-0.07413986, 0.0, 0.0])] test_otf_engine_1.update_positions() - target_positions=[ np.array([2.53714741,2.5,2.5]), - np.array([4.46285259,2.5,2.5]) - ] + target_positions = [np.array([2.53714741, 2.5, 2.5]), + np.array([4.46285259, 2.5, 2.5]) + ] + + for i, pos in enumerate(test_otf_engine_1.structure.positions): + assert np.isclose(pos, target_positions[i], rtol=1e-6).all() - for i,pos in enumerate(test_otf_engine_1.structure.positions): - assert np.isclose(pos,target_positions[i],rtol=1e-6).all() # ------------------------------------------------------ # test otf runs @@ -260,4 +266,4 @@ def test_otf_1(): otf.gp = Fake_GP(kernel='') otf.run() cleanup_otf_run() -""" \ No newline at end of file +""" diff --git a/tests/test_files/qe_input_1.in b/tests/test_files/qe_input_1.in index ccd22336c..13bc2e6ba 100644 --- a/tests/test_files/qe_input_1.in +++ b/tests/test_files/qe_input_1.in @@ -10,7 +10,7 @@ ecutwfc = 50 ecutrho = 100 ntyp = 1 - nat = 2 +nat = 2 ibrav = 0 / &ELECTRONS diff --git a/tests/test_punchout.py b/tests/test_punchout.py index a968db37f..2c81ef5d2 100644 --- a/tests/test_punchout.py +++ b/tests/test_punchout.py @@ -14,7 +14,6 @@ from test_GaussianProcess import get_random_structure from struc import Structure - # ----------------------- # Test helper methods # ----------------------- diff --git a/tests/test_rcut.py b/tests/test_rcut.py index e69de29bb..5fe7869cf 100644 --- a/tests/test_rcut.py +++ b/tests/test_rcut.py @@ -0,0 +1,81 @@ +#!/usr/bin/env python3 +# pylint: disable=redefined-outer-name + +"""" +Test suite for the rcut.py file, which systematically determines a good cutoff +radius + +Steven Torrisi +""" + +import pytest + +from numpy import arccos, zeros, arctan,isclose,pi,eye,equal,not_equal +from numpy.linalg import norm +from numpy.random import randint,normal,uniform +from test_GaussianProcess import get_random_structure + +import sys +sys.path.append('../src') +from struc import Structure +from qe_util import run_espresso +from rcut import perturb_position,perturb_outside_radius + + +def test_perturb(): + """ + Tests to make sure the randomly distributed positions are actually + uniformly distributed over the surface of a sphere of desired avg radius + """ + trials = 200000 + + running_phi = 0 + running_theta = 0 + running_radius = 0 + for n in range(trials): + pos = zeros(3) + newpos = perturb_position(pos,r=1,rscale=.1) + + rad = norm(newpos) + + running_theta += arccos(newpos[2]/rad) + running_phi += arctan(newpos[1]/newpos[0]) + running_radius+= rad + + avg_phi =running_phi/trials + avg_theta =running_theta/trials + avg_rad = running_radius/trials + assert isclose(avg_phi,0,atol= .2) + assert isclose(avg_theta,pi/2,atol=.2) + assert isclose(avg_rad,1,atol=.2) + + +def test_perturb_outside_radius(): + + # Try for 10 different random structures + + for _ in range(10): + d = abs(normal(10,2)) + cell = d * eye(3) + noa = randint(5,30) + positions = [uniform(0,d,3) for _ in range(noa)] + + struc = Structure(lattice=cell,species=['A']*noa,positions=positions, + cutoff=1) + target_index = randint(0, noa) + target_pos = struc.positions[target_index] + r_fix = abs(normal(5.0,1.0)) + within_radius = [norm(target_pos-pos)