diff --git a/AppOutputExtractor/FHIaims/FHIaimsOutputExtractor.py b/AppOutputExtractor/FHIaims/FHIaimsOutputExtractor.py index 633bec2..33055c6 100644 --- a/AppOutputExtractor/FHIaims/FHIaimsOutputExtractor.py +++ b/AppOutputExtractor/FHIaims/FHIaimsOutputExtractor.py @@ -1,5 +1,10 @@ -# +''' +Author: Dr Woongkyu Jee, Dong-Gi Kang +''' + +# +import time from AppOutputExtractor.OutputExtractor import BaseExtractor from AppOutputExtractor.FHIaims.FHIaimsMolecule import molecule as fmol from AppOutputExtractor.FHIaims.FHIaimsMolecule import calculate_rmsd_molecules @@ -8,334 +13,548 @@ import ParsingSupport import os,re +import numpy as np import string,json class extractor(BaseExtractor): - def __init__(self,app_version='22',tag=None): - ''' - ''' - super().__init__(app='FHIaims',version=app_version) - - # set app output patterns - module_path = os.path.dirname(os.path.abspath(__file__)) + '/OutputPattern' # getting this module path, '__file__' - self.patterns = self.load_patterns(module_path) - - # memo - self.tag = tag - - # shellcommand obj - self.shell = shellcommand() - - def set_output_filepath(self,path): - - if os.path.exists(path): - self.output_filepath = path - else: - self.output_filepath = None - print('in {} method "set_output_filepath()", cannot find the file at: "{}" '.format(__file__,path)) - - def set_input_geometry_filepath(self,path): - - if os.path.exists(path): - self.input_geometry_filepath = path - self.input_geometry = fmol(path) - else: - self.input_geometry_filepath = None - print('in {} method "set_input_geometry_filepath()", cannot find the file at: "{}" '.format(__file__,path)) - - def set_output_geometry_filepath(self,path): - - if os.path.exists(path): - self.output_geometry_filepath = path - self.output_geometry = fmol(path) - else: - self.output_geometry_filepath = None - print('in {} method "set_output_geometry_filepath()", cannot find the file at: "{}" '.format(__file__,path)) - - def check_filepaths(self): - #print('App output : {}'.format(self.output_filepath)) - #print('geometry input : {}'.format(self.input_geometry_filepath)) - #print('geometry output: {}'.format(self.output_geometry_filepath)) - ''' - field: (0) AppOutput (1) InputGeometry (2) OutputGeometry - ''' - return [self.output_filepath,self.input_geometry_filepath,self.output_geometry_filepath] - - def get_input_molecule(self): - checker = self.check_filepaths()[1] - if checker: - return self.input_geometry - else: - print('input geometry is not loaded!') - - def get_output_molecule(self): - checker = self.check_filepaths()[2] - if checker: - return self.output_geometry - else: - print('output geometry is not loaded!') - - ''' - Interaction with app output file - ''' - - def check_calculation_success(self): - - self.shell.set_tarfile(self.output_filepath) - cmd = self.shell.grep(self.patterns['SUCCESS']['pattern']) - shell_res = self.shell.execute(cmd) - - if shell_res != None: - self.output_success_tag = True - else: - self.output_success_tag = False - - return self.output_success_tag #!!! - - def check_calculation_runtime(self): - # wall clock time - self.shell.set_tarfile(self.output_filepath) - cmd = self.shell.pipe(self.shell.grep(self.patterns['APP_RUNTIME']['pattern']),self.shell.awk(self.patterns['APP_RUNTIME']['wtime_token'])) - target = self.shell.execute(cmd) #!!! - - try: - target = float(target) - return target - except: - print('failed to get calculation wtime') - return None - - def check_parallel_task(self): - # used cpus - self.shell.set_tarfile(self.output_filepath) - cmd = self.shell.pipe(self.shell.grep(self.patterns['APP_RESOURCE_USED']['pattern']),self.shell.awk(self.patterns['APP_RESOURCE_USED']['token'])) - target = self.shell.execute(cmd) #!!! - - try: - target = int(target) - return target - except: - print('failed to get parallel task number, recheck the app output file') - return None - ''' - Loading SCF converged blocks ... possibly useful for further app output collation - ''' - def set_scf_blocks(self): - ''' - * special blocks: - self.scf_converged_blocks[0] -> first SCF converged blocks [line_start,line_end] - self.scf_converged_blocks[-1]-> final SCF converged blocks - ''' - pattern = self.patterns['BEGIN_SCF']['pattern'].replace("'","") - self.total_lnumber, self.scf_block_lines = ParsingSupport.find_pattern_with_last_word(self.output_filepath,pattern) # GET LINE NUMBERS OF SCF (CONVERGED) BLOCKS - - self.scf_converged_blocklines = [] #!!! - self.scf_converged_blocks = [] #!!! - - # IF ITEM IN ITERABLE SAVE THE LINE NUMBEERS [START,END] - for i, item in enumerate(self.scf_block_lines[:-1]): - curr_tag = int(self.scf_block_lines[i][1]) - next_tag = int(self.scf_block_lines[i+1][1]) - - if next_tag < curr_tag: - - block_start = self.scf_block_lines[i][0] - block_end = self.scf_block_lines[i+1][0] - - self.scf_converged_blocklines.append([block_start,block_end]) - - # FIANL SCF CONVERGED BLOCK (BEFORE APP FINALISATION) - block_start = self.scf_block_lines[-1][0] - block_end = self.total_lnumber - self.scf_converged_blocklines.append([block_start,block_end]) - - # SAVE THE BLOCKS ... 'self.scf_converged_blocks' -> python list - for item in self.scf_converged_blocklines: - self.scf_converged_blocks.append(ParsingSupport.get_lines(self.output_filepath,item[0],item[1])) - - def get_scf_blocks(self): - return self.scf_converged_blocks - - def get_number_of_scf_blocks(self): - return len(self.scf_converged_blocks) - - ''' - AppOutput Collation Methods - ''' - - def get_total_energy(self,block=-1): - - pattern_str = self.patterns['SCF_ENERGY']['pattern'].replace("'","") - token = int(self.patterns['SCF_ENERGY']['token']) - 1 - pattern = re.compile(pattern_str) - - for line in self.scf_converged_blocks[block]: - matching = pattern.search(line) - if matching: - target = float(line.split()[ token ]) - break - return target - - def get_dipole(self,block=-1): - - pattern_str = self.patterns['DIPOLE']['pattern'].replace("'","") - token = int(self.patterns['DIPOLE']['token']) - 1 - pattern = re.compile(pattern_str) - - for line in self.scf_converged_blocks[block]: - matching = pattern.search(line) - if matching: - target = float(line.split()[ token ]) - break - return target - - def get_dipole_moment(self,block=-1): - - pattern_str = self.patterns['DIPOLE_MOMENT']['pattern'].replace("'","") - token_x = int(self.patterns['DIPOLE_MOMENT']['token_x']) - 1 - token_y = int(self.patterns['DIPOLE_MOMENT']['token_y']) - 1 - token_z = int(self.patterns['DIPOLE_MOMENT']['token_z']) - 1 - pattern = re.compile(pattern_str) - - target = [] - - for line in self.scf_converged_blocks[block]: - matching = pattern.search(line) - if matching: - target.append( float(line.split()[ token_x ]) ) - target.append( float(line.split()[ token_y ]) ) - target.append( float(line.split()[ token_z ]) ) - break - return target - - def get_homolumo(self,block=-1): - - ''' - field: (0) HOMO (1) LUMO (2) HOMO-LUMO - ''' - target = [] - - # HOMO - pattern_str = self.patterns['HOMO']['pattern'].replace("'","") - token = int(self.patterns['HOMO']['token']) - 1 - pattern = re.compile(pattern_str) - for line in self.scf_converged_blocks[block]: - matching = pattern.search(line) - if matching: - target.append( float(line.split()[ token ]) ) - break - # LUMO - pattern_str = self.patterns['LUMO']['pattern'].replace("'","") - token = int(self.patterns['LUMO']['token']) - 1 - pattern = re.compile(pattern_str) - for line in self.scf_converged_blocks[block]: - matching = pattern.search(line) - if matching: - target.append( float(line.split()[ token ]) ) - break - # HOMOLUMO GAP - pattern_str = self.patterns['HOMOLUMO']['pattern'].replace("'","") - token = int(self.patterns['HOMOLUMO']['token']) - 1 - pattern = re.compile(pattern_str) - for line in self.scf_converged_blocks[block]: - matching = pattern.search(line) - if matching: - target.append( float(line.split()[ token ]) ) - break - return target - - - # Getters Miscs - - def get_patterns(self): - # return type 'json' - return self.patterns - - def get_tag(self): - return self.tag + def __init__(self,app_version='22',tag=None): + ''' + ''' + super().__init__(app='FHIaims',version=app_version) + + # set app output patterns + module_path = os.path.dirname(os.path.abspath(__file__)) + '/OutputPattern' # getting this module path, '__file__' + self.patterns = self.load_patterns(module_path) + + # memo + self.tag = tag + + # shellcommand obj + self.shell = shellcommand() + self.scf_converged_blocks = [] + + def set_output_filepath(self, path): + if os.path.exists(path): + self.output_filepath = path + else: + self.output_filepath = None + print('in {} method "set_output_filepath()", cannot find the file at: "{}" '.format(__file__,path)) + + def set_input_geometry_filepath(self,path): + if os.path.exists(path): + self.input_geometry_filepath = path + self.input_geometry = fmol(path) + else: + self.input_geometry_filepath = None + print('in {} method "set_input_geometry_filepath()", cannot find the file at: "{}" '.format(__file__,path)) + + def set_output_geometry_filepath(self,path): + if os.path.exists(path): + self.output_geometry_filepath = path + self.output_geometry = fmol(path) + else: + self.output_geometry_filepath = None + print('in {} method "set_output_geometry_filepath()", cannot find the file at: "{}" '.format(__file__,path)) + + def check_filepaths(self): + #print('App output : {}'.format(self.output_filepath)) + #print('geometry input : {}'.format(self.input_geometry_filepath)) + #print('geometry output: {}'.format(self.output_geometry_filepath)) + ''' + field: (0) AppOutput (1) InputGeometry (2) OutputGeometry + ''' + return [self.output_filepath,self.input_geometry_filepath,self.output_geometry_filepath] + + def get_input_molecule(self): + checker = self.check_filepaths()[1] + if checker: + return self.input_geometry + else: + print('input geometry is not loaded!') + + def get_output_molecule(self): + checker = self.check_filepaths()[2] + if checker: + return self.output_geometry + else: + print('output geometry is not loaded!') + + + + ''' + Interaction with app output file + ''' + + def check_calculation_success(self): + self.shell.set_tarfile(self.output_filepath) + cmd = self.shell.grep(self.patterns['SUCCESS']['pattern']) + shell_res = self.shell.execute(cmd) + + if shell_res != None: + self.output_success_tag = True + else: + self.output_success_tag = False + + return self.output_success_tag #!!! + + def check_calculation_runtime(self): + # wall clock time + self.shell.set_tarfile(self.output_filepath) + cmd = self.shell.pipe(\ + self.shell.grep(self.patterns['APP_RUNTIME']['pattern'])\ + ,self.shell.awk(self.patterns['APP_RUNTIME']['wtime_token']) + ) + target = self.shell.execute(cmd) #!!! + + try: + target = float(target) + return target + except: + print('failed to get calculation wtime') + return None + + def check_parallel_task(self): + # used cpus + self.shell.set_tarfile(self.output_filepath) + cmd = self.shell.pipe(\ + self.shell.grep(self.patterns['APP_RESOURCE_USED']['pattern'])\ + ,self.shell.awk(self.patterns['APP_RESOURCE_USED']['token']) + ) + target = self.shell.execute(cmd) #!!! + + try: + target = int(target) + return target + except: + print('failed to get parallel task number, recheck the app output file') + return None + + + + ''' + Loading SCF converged blocks ... possibly useful for further app output collation + ''' + + #@property + def set_scf_blocks(self) -> list: + ''' + * special blocks: + self.scf_converged_blocks[0] -> first SCF converged blocks [line_start,line_end] + self.scf_converged_blocks[-1]-> final SCF converged blocks + ''' + pattern = self.patterns['BEGIN_SCF']['pattern'].replace("'","") + self.total_lnumber, self.scf_block_lines = \ + ParsingSupport.find_pattern_with_last_word(self.output_filepath,pattern) \ + # GET LINE NUMBERS OF SCF (CONVERGED) BLOCKS + + self.scf_converged_blocklines = [] #!!! + self.scf_converged_blocks = [] #!!! + + # IF ITEM IN ITERABLE SAVE THE LINE NUMBEERS [START,END] + for i, item in enumerate(self.scf_block_lines[:-1]): + curr_tag = int(self.scf_block_lines[i][1]) + next_tag = int(self.scf_block_lines[i+1][1]) + + if next_tag < curr_tag: + + block_start = self.scf_block_lines[i][0] + block_end = self.scf_block_lines[i+1][0] + + self.scf_converged_blocklines.append([block_start,block_end]) + + # FIANL SCF CONVERGED BLOCK (BEFORE APP FINALISATION) + block_start = self.scf_block_lines[-1][0] + block_end = self.total_lnumber + self.scf_converged_blocklines.append([block_start,block_end]) + + # SAVE THE BLOCKS ... 'self.scf_converged_blocks' -> python list + for item in self.scf_converged_blocklines: + self.scf_converged_blocks.append(ParsingSupport.get_lines(self.output_filepath,item[0],item[1])) + return self.scf_converged_blocks + + + #@property + def get_species(self, atom_order) -> list: + #get_species = [x for x in self.get_atom_order.tolist()] + get_species = list(set([item for sublist in atom_order for item in sublist])) + get_species = sorted(get_species) + return get_species + + + #@property + def get_no_atoms(self) -> int: + with open(self.output_filepath, 'r') as f: + lines = f.readlines() + for i in lines: + if self.patterns['NO_ATOMS']['pattern'] in i: + no_atoms = int(i.split()[5]) + return no_atoms + + + # REVIEW-delete: property decorator on scf_converged_blocks does same function + #def get_scf_blocks(self): + # return self.scf_converged_blocks + + + @property + def get_number_of_scf_blocks(self) -> int: + return len(self.set_scf_blocks) + + ''' + AppOutput Collation Methods + ''' + + def get_total_energy(self, block=-1): + pattern_str = self.patterns['SCF_ENERGY']['pattern'].replace("'","") + token = int(self.patterns['SCF_ENERGY']['token']) - 1 + pattern = re.compile(pattern_str) + for i in self.set_scf_blocks[block]: + matching = pattern.search(i) + if matching: + target = float(i.split()[ token ]) + break + return target + + + + #@property + def get_atom_order(self, no_atoms, block=-1) -> np.ndarray: + self.set_scf_blocks() + + pattern_str = self.patterns['SCF_GEOMETRY_END']['pattern'].replace("'", "") + pattern = re.compile(pattern_str) + + start_index = None + self.match_atom = np.empty((no_atoms), dtype=object) + for numj, j in enumerate(self.scf_converged_blocks[block]): + matching = pattern.search(j) + if matching: + start_index = numj + 2 + elif start_index is not None and j.strip() == '': + end_index = numj - 1 + atomic_structure = self.scf_converged_blocks[block][start_index: end_index] + for numk, k in enumerate(atomic_structure): + numbers = [x for x in k.split()] + self.match_atom[numk] = numbers[-1] + break + self.match_atom = np.reshape(self.match_atom, (no_atoms, 1)) + return self.match_atom + + #@property + def get_geometries(self, no_atoms, block=-1) -> np.ndarray: + #if not self.scf_converged_blocks: + self.set_scf_blocks() + + pattern_str = self.patterns['SCF_GEOMETRY_BEGIN']['pattern'].replace("'", "") + pattern = re.compile(pattern_str) + + if block == -1 or block == len(self.scf_converged_blocks)-1: + pattern_str = self.patterns['SCF_GEOMETRY_END']['pattern'].replace("'", "") + pattern = re.compile(pattern_str) + + start_index = None + self.geo = np.zeros((no_atoms, 3)) + cnt = 0 + for numj, j in enumerate(self.scf_converged_blocks[block]): + matching = pattern.search(j) + if matching: + start_index = numj + 2 + end_index = numj + no_atoms + 2 + atomic_structure = self.scf_converged_blocks[block][start_index: end_index] + for numk, k in enumerate(atomic_structure): + numbers = [x for x in k.split()] + self.geo[numk] = list(map(float, numbers[1:4])) # Convert the rest to float and store in self.geo + cnt += 1 + start_index = None + + if cnt == 0: + cnt = 1 + else: pass + self.geo = np.reshape(self.geo, (cnt, int(no_atoms), 3)) + return self.geo + + def get_sp_geometries(self, path) -> np.ndarray: + new_path = os.path.join(os.path.dirname(path), 'geometry.in') + with open(new_path, 'r') as f: + lines = f.readlines() + lines = [x.split() for x in lines] + lines = np.array(lines) + shape = np.shape(lines) + #atom_str = np.reshape(lines[:,0], (shape[0], -1)) + self.atom_label = np.reshape(lines[:,-1], (shape[0], -1)) + self.coordinate = lines[:, 1:-1].astype(float) + return self.coordinate, self.atom_label + + def get_sp_forces(self, no_atoms, path) -> np.ndarray: + pattern_str = self.patterns['SCF_FORCE']['pattern'].replace("'","") + pattern = re.compile(pattern_str) + start_index = None + self.forces = np.zeros((no_atoms, 3)) + cnt = 0 + with open(path, 'r') as f: + lines = f.readlines() + + for numi, i in enumerate(lines): + matching = pattern.search(i) + if matching: + start_index = numi + 1 + elif start_index is not None and i.strip() == '': + end_index = numi + force = lines[start_index:end_index] + + for numj, j in enumerate(force): + numbers = list(map(float, j.strip().split()[-3:])) + + self.forces[numj] = numbers + start_index = None + cnt += 1 + self.forces = np.reshape(self.forces, (-1, 3)) + return self.forces + + def get_sp_total_energy(self, path): + pattern_str = self.patterns['SCF_ENERGY']['pattern'].replace("'","") + token = int(self.patterns['SCF_ENERGY']['token']) - 1 + pattern = re.compile(pattern_str) + with open(path, 'r') as f: + lines = f.readlines() + for i in lines: + matching = pattern.search(i) + if matching: + target = float(i.split()[ token ]) + break + return target + + def get_sp_atom_order(self): + return self.atom_label + + def get_sp_species(self): + return list(set(self.atom_label.flatten().tolist())) + + #@property + def get_forces(self, no_atoms=12, block=-1) -> np.ndarray: + pattern_str = self.patterns['SCF_FORCE']['pattern'].replace("'", "") + pattern = re.compile(pattern_str) + start_index = None + self.forces = np.zeros((no_atoms, 3)) + cnt = 0 + #print(self.scf_converged_blocks) + #print(len(self.scf_converged_blocks)) + #print(block) + #for i in self.scf_converged_blocks[block]: + # for j in i: + # print(j, end='') + if not self.scf_converged_blocks: + raise ValueError("scf_converged_blocks is empty!") + + for numj, j in enumerate(self.scf_converged_blocks[block]): + matching = pattern.search(j) + if matching: + start_index = numj + 1 + elif start_index is not None and j.strip() == '': + end_index = numj + force = self.scf_converged_blocks[block][start_index:end_index] + + for numk, k in enumerate(force): + numbers = list(map(float, k.strip().split()[-3:])) + self.forces[numk] = numbers + start_index = None + cnt += 1 + + self.forces = np.reshape(self.forces, (cnt, -1, 3)) + return self.forces + + + #@property + def get_vib_eigvec(self, no_atoms) -> np.ndarray: + ''' + Read whole file contents (not necessary to read in blocks as we need all eigenvector of vibrational mode + ''' + check_vib = [x for x in os.listdir('./') if 'vibration' in x] + if len(check_vib) == 0: + raise FileNotFoundError("Cannot find 'vibration' directory") + else: + check_vib = [x for x in os.listdir('./') if 'vibration' in x][0] + + vib_xyz = [os.path.join(check_vib, x) for x in os.listdir(check_vib) if '_' and '.xyz' in x][0] + + with open(vib_xyz, 'r') as f: + lines = f.readlines() + + self.eigvec = np.zeros((no_atoms*3, no_atoms, 3)) + start_index = None + block_counter = -1 + for numi, i in enumerate(lines): + if 'frequency' in i: + start_index = numi + 1 + block_counter += 1 + elif start_index is not None and (i.strip().split()[0] in ['Al', 'F']): + data = list(map(float, i.strip().split()[-3:])) # convert last three elements to float + atom_index = numi - start_index + self.eigvec[block_counter, atom_index, :] = data + elif i.strip() == str(self.get_no_atoms): + start_index = None + + return self.eigvec + + + def get_dipole(self,block=-1): + pattern_str = self.patterns['DIPOLE']['pattern'].replace("'","") + token = int(self.patterns['DIPOLE']['token']) - 1 + pattern = re.compile(pattern_str) + + for line in self.scf_converged_blocks[block]: + matching = pattern.search(line) + if matching: + target = float(line.split()[ token ]) + break + return target + + def get_dipole_moment(self,block=-1): + pattern_str = self.patterns['DIPOLE_MOMENT']['pattern'].replace("'","") + token_x = int(self.patterns['DIPOLE_MOMENT']['token_x']) - 1 + token_y = int(self.patterns['DIPOLE_MOMENT']['token_y']) - 1 + token_z = int(self.patterns['DIPOLE_MOMENT']['token_z']) - 1 + pattern = re.compile(pattern_str) + + target = [] + + for line in self.scf_converged_blocks[block]: + matching = pattern.search(line) + if matching: + target.append( float(line.split()[ token_x ]) ) + target.append( float(line.split()[ token_y ]) ) + target.append( float(line.split()[ token_z ]) ) + break + return target + + def get_homolumo(self,block=-1): + + ''' + field: (0) HOMO (1) LUMO (2) HOMO-LUMO + ''' + target = [] + + # HOMO + pattern_str = self.patterns['HOMO']['pattern'].replace("'","") + token = int(self.patterns['HOMO']['token']) - 1 + pattern = re.compile(pattern_str) + for line in self.scf_converged_blocks[block]: + matching = pattern.search(line) + if matching: + target.append( float(line.split()[ token ]) ) + break + # LUMO + pattern_str = self.patterns['LUMO']['pattern'].replace("'","") + token = int(self.patterns['LUMO']['token']) - 1 + pattern = re.compile(pattern_str) + for line in self.scf_converged_blocks[block]: + matching = pattern.search(line) + if matching: + target.append( float(line.split()[ token ]) ) + break + # HOMOLUMO GAP + pattern_str = self.patterns['HOMOLUMO']['pattern'].replace("'","") + token = int(self.patterns['HOMOLUMO']['token']) - 1 + pattern = re.compile(pattern_str) + for line in self.scf_converged_blocks[block]: + matching = pattern.search(line) + if matching: + target.append( float(line.split()[ token ]) ) + break + return target + + + # Getters Miscs + + def get_patterns(self): + # return type 'json' + return self.patterns + + def get_tag(self): + return self.tag if __name__ == '__main__': - file_root = '/Users/woongkyujee/Desktop/Python/FHI22_samples/runs/run_1' - main_output = file_root + '/FHIaims.out' - input_geo = file_root + '/geometry.in' - output_geo = file_root + '/geometry.in.next_step' - - ext2 = extractor() - ext2.set_output_filepath(main_output) - ext2.set_input_geometry_filepath(input_geo) - ext2.set_output_geometry_filepath(output_geo) - - print('check filepaths()') - print(ext2.check_filepaths()) # if None in ext2.check_filepaths(): - print('calculation success check: {}'.format(ext2.check_calculation_success())) - - print('calculation runtime') - rtime = ext2.check_calculation_runtime() - print(rtime) - - print('calculation parallel tasks') - ptask = ext2.check_parallel_task() - print(ptask) - - - - - - - - ### EXTRACTION - - ext2.set_scf_blocks() # load scf blocks - - # Energy Check - print('init E') - init_E = ext2.get_total_energy(0) - print(init_E) - print('final E') - final_E = ext2.get_total_energy() - print(final_E) - - # Dipole Check - print('init P') - init_p = ext2.get_dipole(0) - print(init_p) - initial_p_elem = ext2.get_dipole_moment(0) - print(initial_p_elem) - - print('final P') - final_p = ext2.get_dipole() - print(final_p) - final_p_elem = ext2.get_dipole_moment(0) - print(final_p_elem) - - # HOMOLUMO CHECK - print('init homo-lumo, list [homo,lumo,homo-lumo]') - init_hl = ext2.get_homolumo(0) - print(init_hl) - print('final homo-lumo, list [homo,lumo,homo-lumo]') - final_hl = ext2.get_homolumo() - print(final_hl) - - ''' - Unit test with 'ext2' instance - ''' - print('--- input') - ext2.input_geometry.show_info() - print('--- output') - ext2.output_geometry.show_info() - print('in - out geometry rmsd') - rmsd = calculate_rmsd_molecules(ext2.input_geometry,ext2.output_geometry) - print(rmsd) - - #print('output check --') - #ext2.check_output_success() - #print(ext2.output_success_tag) - - ''' - Unit test getting SCF Blocks - ''' + file_root = '/Users/woongkyujee/Desktop/Python/FHI22_samples/runs/run_1' + main_output = file_root + '/FHIaims.out' + input_geo = file_root + '/geometry.in' + output_geo = file_root + '/geometry.in.next_step' + + ext2 = extractor() + ext2.set_output_filepath(main_output) + ext2.set_input_geometry_filepath(input_geo) + ext2.set_output_geometry_filepath(output_geo) + + print('check filepaths()') + print(ext2.check_filepaths()) # if None in ext2.check_filepaths(): + print('calculation success check: {}'.format(ext2.check_calculation_success())) + + print('calculation runtime') + rtime = ext2.check_calculation_runtime() + print(rtime) + + print('calculation parallel tasks') + ptask = ext2.check_parallel_task() + print(ptask) + + + + + + + + ### EXTRACTION + + ext2.set_scf_blocks() # load scf blocks + + # Energy Check + print('init E') + init_E = ext2.get_total_energy(0) + print(init_E) + print('final E') + final_E = ext2.get_total_energy() + print(final_E) + + # Dipole Check + print('init P') + init_p = ext2.get_dipole(0) + print(init_p) + initial_p_elem = ext2.get_dipole_moment(0) + print(initial_p_elem) + + print('final P') + final_p = ext2.get_dipole() + print(final_p) + final_p_elem = ext2.get_dipole_moment(0) + print(final_p_elem) + + # HOMOLUMO CHECK + print('init homo-lumo, list [homo,lumo,homo-lumo]') + init_hl = ext2.get_homolumo(0) + print(init_hl) + print('final homo-lumo, list [homo,lumo,homo-lumo]') + final_hl = ext2.get_homolumo() + print(final_hl) + + ''' + Unit test with 'ext2' instance + ''' + print('--- input') + ext2.input_geometry.show_info() + print('--- output') + ext2.output_geometry.show_info() + print('in - out geometry rmsd') + rmsd = calculate_rmsd_molecules(ext2.input_geometry,ext2.output_geometry) + print(rmsd) + + #print('output check --') + #ext2.check_output_success() + #print(ext2.output_success_tag) + + ''' + Unit test getting SCF Blocks + ''' diff --git a/AppOutputExtractor/FHIaims/FHIaimsVib.py b/AppOutputExtractor/FHIaims/FHIaimsVib.py new file mode 100644 index 0000000..f48cf27 --- /dev/null +++ b/AppOutputExtractor/FHIaims/FHIaimsVib.py @@ -0,0 +1,137 @@ +from AppOutputExtractor.OutputExtractor import BaseExtractor +from AppOutputExtractor.FHIaims.FHIaimsMolecule import molecule as fmol +from AppOutputExtractor.FHIaims.FHIaimsMolecule import calculate_rmsd_molecules +from AppOutputExtractor.FHIaims.FHIaimsOutputExtractor import extractor + +from ShellCommand import shellcommand +import ParsingSupport + +import os +import shutil +import string,json + + +class aimsvibcalc(BaseExtractor): + + def __init__(self, app_version='22', tag=None): + ''' + ''' + app_output = './aims.out' + self.extractor = extractor() + self.extractor.set_output_filepath(app_output) + self.no_atoms = self.extractor.get_no_atoms() + self.geometries = self.extractor.get_geometries(self.no_atoms) + self.order = self.extractor.get_atom_order(self.no_atoms) + self.species = self.extractor.get_species(self.order) + + + self.ucl_id = 'uccatka' + self.job_time = '2:00:00' + #self.job_name = 'test' + self.memory = '2' + self.cpu_core = '40' # for Young 40 core = 1 node + self.payment = 'Gold' + self.budgets = 'UCL_chemM_Woodley' + self.path_binary = '/home/uccatka/software/fhi-aims.221103/build/aims.221103.scalapack.mpi.x' + self.vib_path_binary = '/home/uccatka/software/fhi-aims.221103/build/src/vibrations/numerical_vibrations.pl' + self.path_fhiaims_species = '/home/uccatka/software/fhi-aims.221103/species_defaults/defaults_2020/light' + self.step_size = 0.05 + + super().__init__(app='FHIaims',version=app_version) + + # set app output patterns + module_path = os.path.dirname(os.path.abspath(__file__)) + '/OutputPattern' # getting this module path, '__file__' + self.patterns = self.load_patterns(module_path) + + # memo + self.tag = tag + + # shellcommand obj + self.shell = shellcommand() + + return None + + def make_job_submit(self, job_name, loc='./vibration', step_size='0.0025'): + ''' Write 'submit.sh' job script for SGE system ''' + path = os.path.join(loc, 'submit.sh') + with open(path, 'a') as f: + f.write("#!/bin/bash -l\n") + f.write("\n") + f.write("#$ -S /bin/bash\n") + f.write(f"#$ -l h_rt={self.job_time}\n") + f.write(f"#$ -l mem={self.memory}G\n") + f.write(f"#$ -N {job_name}\n") + f.write(f"#$ -pe mpi {self.cpu_core}\n") + f.write("#$ -cwd\n") + f.write("\n") + f.write(f"#$ -P {self.payment}\n") + f.write(f"#$ -A {self.budgets}\n") + f.write("\n") + f.write("#$ -m e\n") + f.write(f"#$ -M {self.ucl_id}@ucl.ac.uk\n") + f.write("\n") + f.write("module purge\n") + f.write("module load gerun\n") + f.write("module load userscripts\n") + f.write("module load gcc-libs/4.9.2\n") + f.write("module unload -f compilers mpi\n") + f.write("module load beta-modules\n") + f.write("module load gcc-libs/10.2.0\n") + f.write("module load openblas/0.3.7-serial/gnu-4.9.2\n") + f.write("module load compilers/intel/2019/update5\n") + f.write("module load mpi/intel/2018/update3/intel\n") + f.write("module load cmake/3.21.1\n\n") + + f.write(f"{self.vib_path_binary} {job_name}_{step_size} {step_size} > vibres.out\n") + + @property + def vib_calc_prep(self): + #shutil.copy('./control.in', './vibration') + geo_next = 'geometry.in.next_step' + geo = 'geometry.in' + vib_dir = 'vibration' + geometry_files = [x for x in os.listdir('./') if '.in' in x] + if os.path.exists(geo_next): + shutil.copy(geo_next, f'{vib_dir}/{geo}') + shutil.copy('hessian.aims', vib_dir) + else: + shutil.copy(geo, vib_dir) + + basis_set_files = [os.path.join(self.path_fhiaims_species, x) for x in os.listdir(self.path_fhiaims_species)] + basis_set_all = [x.split('_')[1] for x in os.listdir(self.path_fhiaims_species)] + basis_set_index = [basis_set_all.index(x) for x in basis_set_all if x in self.species] + + with open(os.path.join(vib_dir, 'control.in'), 'a') as f: + f.write("xc pbesol\n") + f.write("spin none\n") + f.write("relativistic atomic_zora scalar\n") + f.write("charge 0.\n\n") + f.write("# SCF convergence\n") + f.write("occupation_type gaussian 0.01\n") + f.write("mixer pulay\n") + f.write("n_max_pulay 10\n") + f.write("charge_mix_param 0.5\n") + f.write("sc_accuracy_rho 1E-5\n") + f.write("sc_accuracy_eev 1E-3\n") + f.write("sc_accuracy_etot 1E-6\n") + f.write("sc_accuracy_forces 1E-4\n") + f.write("sc_iter_limit 1500\n\n") + + for i in basis_set_index: + with open(basis_set_files[i], 'r') as ff: + lines = ff.read() + f.write(lines) + f.write('\n') + return None + +if __name__ == "__main__": + vib = aimsvibcalc() + os.mkdir('vibration') + vib.vib_calc_prep + current_dir_name = os.path.basename(os.getcwd()) + vib.make_job_submit(f'n{current_dir_name}') + os.chdir('vibration') + os.system('qsub submit.sh') + + + diff --git a/AppOutputExtractor/FHIaims/MLTTV_spliter.py b/AppOutputExtractor/FHIaims/MLTTV_spliter.py new file mode 100644 index 0000000..0fc84a2 --- /dev/null +++ b/AppOutputExtractor/FHIaims/MLTTV_spliter.py @@ -0,0 +1,37 @@ +from ase import io + +def split_xyz_file(input_file, train_file, valid_file, test_file): + with open(input_file, 'r') as infile: + train_out = open(train_file, 'w') + valid_out = open(valid_file, 'w') + test_out = open(test_file, 'w') + + while True: + # Read the header line containing the number of atoms + header = infile.readline() + if not header: + break # End of file + + num_atoms = int(header.strip()) + block_lines = [infile.readline() for _ in range(num_atoms + 1)] + + # Determine which file to write to based on the current index + i = infile.tell() # Get current position in file + if i % 5 < 3: + output_file = train_out + elif i % 5 == 3: + output_file = valid_out + else: + output_file = test_out + + # Write the block to the chosen file + output_file.write(header) + output_file.writelines(block_lines) + + train_out.close() + valid_out.close() + test_out.close() + +# Usage: +split_xyz_file('Training_set.xyz', 'Training_set_test.xyz', 'Validation_set_test.xyz', 'Testing_set_test.xyz') + diff --git a/AppOutputExtractor/FHIaims/MLTrainingDataGenerator.py b/AppOutputExtractor/FHIaims/MLTrainingDataGenerator.py new file mode 100644 index 0000000..b11fa9b --- /dev/null +++ b/AppOutputExtractor/FHIaims/MLTrainingDataGenerator.py @@ -0,0 +1,604 @@ + +""" +dev note: +work on breathing method""" + +''' +Author: Dong-Gi Kang +Prepare ML-IP data using FHI-aims output +Training data type: vibrational mode of a cluster + +[it retreive the vibrational mode cluster geometry and forces from single point calculation and +generates extended xyz format of Training_set.xyz: contains the total energy, atomic coordination, atomic forces] +The Training_set.xyz will be placed in FIT directory and each vibrational mode ext xyz files are generated in ext_xyz directory (The directories are automatically generated from the code) + +N.B. Change the UCL id, budget code, executable path for FHI-aims and fhi-aims species directory + + + +help: +python MLTrainingDataGenerator.py -h +(execute the file at the directory where the geometry.in, control.in, viration (dir) located) + + +1. python {code.py} --mode run --eigenvector="7 8 9 10" would grab 7th, 8th, 9th 10th (can selectively) then, modify the GM with the step_size (GM geometry + eigenvector * step_size) and prepare the individual directories and submit the single point calculations. + +2. python {code.py} --mode retrieve --eigenvector="7 8 9 10" would grab the generated data from the [1.] and make the ext xyz for each vibrational mode and store into the ext_xyz directory + +3. [python {code.py} --mode make_extxyz] would grab the all data from ext_xyz and make ext xyz format of Training_set.xyz in FIT directory + +4. if you want to trianing MACE or GAP ML-IP use MACE_lib.py or second_GAP.py +''' + +import os +import sys +import random +import numpy as np +import argparse +from itertools import groupby +from AppOutputExtractor.FHIaims.FHIaimsOutputExtractor import extractor + +class ML_train_generator(extractor): + + def __init__(self, app_version='22', tag=None): + + self.breathing_called = False + + #self.extractor = extractor() + #self.extractor.set_output_filepath(app_output) + #self.no_atoms = self.extractor.get_no_atoms() + #self.geometries = self.extractor.get_geometries(self.no_atoms) + #self.order = self.extractor.get_atom_order(self.no_atoms) + #ID = self.extractor.get_species(self.order) + #self.forces = self.extractor.get_forces(self.no_atoms) + #self.vib_eigvecs = self.extractor.get_vib_eigvec(self.no_atoms) + + self.ucl_id = 'uccatka' + self.job_time = '2:00:00' + self.job_name = 'test' + self.memory = '1' + self.cpu_core = '40' # for Young 40 core = 1 node + self.payment = 'Gold' + self.budgets = 'UCL_chemM_Woodley' + self.path_binary = '/home/uccatka/software/fhi-aims.221103/build/aims.221103.scalapack.mpi.x' + self.path_fhiaims_species = '/home/uccatka/software/fhi-aims.221103/species_defaults/defaults_2020/light' + self.step_size = 0.1 ##### STEP SIZE ##### + return None + + + def initiate(self): + app_output = './aims.out' + self.extractor = extractor() + self.extractor.set_output_filepath(app_output) + self.no_atoms = self.extractor.get_no_atoms() + self.geometries = self.extractor.get_geometries(self.no_atoms) + self.order = self.extractor.get_atom_order(self.no_atoms) + self.ID = self.extractor.get_species(self.order) + self.forces = self.extractor.get_forces(self.no_atoms) + self.vib_eigvecs = self.extractor.get_vib_eigvec(self.no_atoms) + + + def mod_xyz_w_vib(self): + ''' Modify LM geometries to the frames of vibrational mode frames ''' + Lambda = len(np.arange(-1, 1+self.step_size, self.step_size)) * self.no_atoms*3 + self.mod_sp = np.zeros((Lambda, self.no_atoms, 3)) + cnt = 0 + for i in range(self.no_atoms * 3): # 3N dimension + for numj, j in enumerate(np.arange(-1, 1+self.step_size, self.step_size)): # -1 to 1 in every step size + j = np.round(j, 2) + frame = self.geometries[-1] + self.vib_eigvecs[i] * j + self.mod_sp[cnt] = np.round(frame, 8) + cnt += 1 + self.mod_sp = np.reshape(self.mod_sp, (self.no_atoms*3, len(np.arange(-1, 1+self.step_size, self.step_size)), self.no_atoms, 3)) + return self.mod_sp + + + def mod_xyz_w_rand_pair_vib(self): + + list_eigvecs = list(range(6, self.no_atoms*3)) + random.shuffle(list_eigvecs) + self.pairs_eigvecs = [[list_eigvecs[i], list_eigvecs[i+1]] for i in range(0, len(list_eigvecs), 2)] + fname_pairs = [f"{x}-{y}" for x, y in self.pairs_eigvecs] + + #Lambda = len(np.arange(-1, 1+self.step_size, self.step_size)) * (self.no_atoms*3-6) # range of steps for all vib. mode, except E(3) + #self.mod_sp_pair = np.zeros((Lambda, self.no_atoms, 3)) + + self.mod_sp_pair = np.zeros((len(self.pairs_eigvecs), len(np.arange(-1, 1+self.step_size, self.step_size)), self.no_atoms, 3)) # range of steps for all vib. mode, except E(3) + + cnt = 0 + for numi, i in enumerate(self.pairs_eigvecs): + for numj, j in enumerate(np.arange(-1, 1+self.step_size, self.step_size)): + j = np.round(j, 2) + frame = self.geometries[-1] + (self.vib_eigvecs[i[0]]+self.vib_eigvecs[i[1]]) * j + #self.mod_sp_pair[cnt] = np.round(frame, 8) + self.mod_sp_pair[numi][numj] = np.round(frame, 8) + cnt += 1 + + self.mod_sp_pair = np.reshape(self.mod_sp_pair, (len(self.pairs_eigvecs), numj+1, self.no_atoms, 3)) + return self.mod_sp_pair, fname_pairs + + + def breathing(self): + scale = np.arange(0.6, 1+self.step_size, self.step_size) + Lambda = len(scale) #* self.no_atoms*3 + self.mod_sp_breath = np.zeros((Lambda, self.no_atoms, 3)) + # shift the centre of mass of the structure to (0, 0, 0) + coord = self.geometries[0] + com = coord.sum(axis=0) + com = com / int(self.no_atoms) + coord_x = np.subtract(coord[:, 0], com[0], out=coord[:, 0]) + coord_y = np.subtract(coord[:, 1], com[1], out=coord[:, 1]) + coord_z = np.subtract(coord[:, 2], com[2], out=coord[:, 2]) + coord = list(zip(coord_x, coord_y, coord_z)) + coord = np.array(coord) + cnt = 0 + + for numj, j in enumerate(scale): + j = np.round(j, 2) + frame = coord * j + self.mod_sp_breath[cnt] = np.round(frame, 8) + cnt += 1 + + self.mod_sp_breath = np.reshape(self.mod_sp_breath, (len(scale), self.no_atoms, 3)) + self.breathing_called = True + return self.mod_sp_breath, scale + + + #@property + def geometry_for_sp(self, mod_sp): + ''' Convert the modified geometry (mod_xyz_w_vib) to {geometry.in} format for FHI-aims ''' + # vibrational modes + if not self.breathing_called: + print("@@@@@@@") + placer = np.full((self.no_atoms, 1), 'atom') + placer_species = np.reshape(self.order, (-1, 1)) + shape = np.shape(mod_sp) + self.for_sp = np.empty((shape[0], shape[1], self.no_atoms, 5), dtype=object) + for i in range(shape[0]): + for j in range(shape[1]): + form = np.concatenate((placer, mod_sp[i][j], placer_species), axis=1) + self.for_sp[i][j] = form + return self.for_sp, self.no_atoms + + # breathing mode + else: + print("*******") + placer_breath = np.full((self.no_atoms, 1), 'atom') + placer_species_breath = np.reshape(self.order, (-1, 1)) + shape_breath = np.shape(mod_sp) + self.for_sp = np.empty((shape_breath[0], shape_breath[1], 5), dtype=object) + for i in range(shape_breath[0]): + form = np.concatenate((placer_breath, mod_sp[i], placer_species_breath), axis=1) + self.for_sp[i] = form + return self.for_sp, self.no_atoms + + + @property + def xyz_from_opti(self): + ''' prepare training data from every SCF converged cycles of a optimisation ''' + train_xyz = 'xyz_from_opti.xyz' + exist = [x for x in os.listdir('./') if train_xyz in x] + if len(exist) != 0: + os.remove(exist[0]) + else: pass + for i in range(len(self.extractor.set_scf_blocks)): + self.energy = self.extractor.get_total_energy(i) + self.geometry = self.geometries[i] + self.force = self.forces[i] + xyz = np.round(np.concatenate((self.geometry, self.force), axis=1), 9) + xyz = np.concatenate((self.order, xyz), axis=1) + + with open('xyz_from_opti.xyz', 'a') as f: + f.write(f'{self.no_atoms}\n') + f.write(f'Properties=species:S:1:pos:R:3:forces:R:3 energy={self.energy} pbc="F F F"\n') + np.savetxt(f, xyz, fmt="%s", delimiter=" ") + print(f"total of {i+1} SCF converged structures are prepared in {train_xyz}") + + + def make_sp_control(self, path): + ''' Write {control.in} file ''' + basis_set_files = [os.path.join(self.path_fhiaims_species, x) for x in os.listdir(self.path_fhiaims_species)] + basis_set_all = [x.split('_')[1] for x in os.listdir(self.path_fhiaims_species)] + basis_set_index = [basis_set_all.index(x) for x in basis_set_all if x in self.ID] + + path = os.path.join(path, 'control.in') + with open(path, 'a') as f: + f.write("#\n") + f.write("xc pbesol\n") + f.write("spin none\n") + f.write("relativistic atomic_zora scalar\n") + f.write("charge 0.\n\n") + f.write("# SCF convergence\n") + f.write("occupation_type gaussian 0.01\n") + f.write("mixer pulay\n") + f.write("n_max_pulay 10\n") + f.write("charge_mix_param 0.5\n") + f.write("sc_accuracy_rho 1E-5\n") + f.write("sc_accuracy_eev 1E-3\n") + f.write("sc_accuracy_etot 1E-6\n") + f.write("sc_accuracy_forces 1E-4\n") + f.write("sc_iter_limit 1500\n") + f.write("# Relaxation\n\n") + #f.write("relax_geometry bfgs 1.e-3\n") + for i in basis_set_index: + with open(basis_set_files[i], 'r') as ff: + lines = ff.read() + f.write(lines) + f.write('\n') + return None + + + def make_job_submit(self, path): + ''' Write 'submit.sh' job script for SGE system ''' + _, last_part = os.path.split(path) + _, second_last_part = os.path.split(os.path.dirname(path)) + + # Combine the last two parts + last_two_parts = f"{second_last_part}_{last_part}" + + path = os.path.join(path, 'submit.sh') + with open(path, 'a') as f: + f.write("#!/bin/bash -l\n") + f.write('\n') + f.write("#$ -S /bin/bash\n") + f.write(f"#$ -l h_rt={self.job_time}\n") + f.write(f"#$ -l mem={self.memory}G\n") + f.write(f"#$ -N p{last_two_parts}\n") + f.write(f"#$ -pe mpi {self.cpu_core}\n") + f.write("#$ -cwd\n") + f.write("\n") + f.write(f"#$ -P {self.payment}\n") + f.write(f"#$ -A {self.budgets}\n") + + f.write("module load gerun\n") + f.write("module load userscripts\n") + f.write("module unload -f compilers mpi gcc-libs\n") + f.write("module load gcc-libs/4.9.2\n") + f.write("module unload -f compilers mpi\n") + f.write("module load beta-modules\n") + f.write("module load openblas/0.3.7-serial/gnu-4.9.2\n") + f.write("module load compilers/intel/2019/update5\n") + f.write("module load mpi/intel/2018/update3/intel\n") + + f.write("\n") + f.write("####$ -m e\n") + f.write(f"####$ -M {self.ucl_id}@ucl.ac.uk\n") + f.write("\n") + f.write(f"gerun {self.path_binary} > aims.out\n") + + + @staticmethod + def sorting_key(path): + parts = path.split('/') + second_key = int(parts[1]) if parts[1] != "breathing" else float('inf') + third_key = float(parts[2].split('_')[1]) # Consider lambda value regardless of the second part + return second_key, third_key + + def retrieve_results_c(self, eigenvectors): + print("---retrieve---") + eigvec_path = [os.path.join('sp', str(eigvec)) for eigvec in eigenvectors] + sp_path = [os.path.join(dirpath, fname) for dirpath in eigvec_path for fname in os.listdir(dirpath)] + lambda_path = [os.path.join(dirpath, fname) for dirpath in sp_path for fname in os.listdir(dirpath) if fname == 'aims.out'] + aims_out_path = sorted(lambda_path, key=self.sorting_key) + aims_out_path = [list(group) for key, group in groupby(aims_out_path, lambda x: x.split('/')[1])] + + if not os.path.exists('ext_xyz'): + os.mkdir('ext_xyz') + cnt = 0 + for numi, i in enumerate(aims_out_path): + filename = f"ext_xyz/ext_{i[0].split('/')[1]}_eigv.xyz" # + with open(filename, 'a') as f: # + for j in i: + ex = extractor() + ex.set_output_filepath(j) + #ex.set_scf_blocks + + no_atoms = ex.get_no_atoms() + geometries, atom_label = ex.get_sp_geometries(j) + forces = ex.get_sp_forces(no_atoms, j) + total_energy = ex.get_sp_total_energy(j) + + coulomb_E, coulomb_F = self.coulomb_E_F(atom_label, geometries) + + # subtract coulomb energy and force + energy = total_energy - coulomb_E + forces = forces - coulomb_F + form = np.concatenate((ex.get_sp_atom_order(), geometries, forces), axis=1) + + f.write(str(no_atoms) + '\n') + f.write(f'Lattice="0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0" Properties=species:S:1:pos:R:3:forces:R:3 energy={energy} pbc="F F F"\n') + np.savetxt(f, form, fmt="%s", delimiter=" ") + + + def retrieve_results(self, eigenvectors): + print("---retrieve---") + eigvec_path = [os.path.join('sp', str(eigvec)) for eigvec in eigenvectors] + sp_path = [os.path.join(dirpath, fname) for dirpath in eigvec_path for fname in os.listdir(dirpath)] + lambda_path = [os.path.join(dirpath, fname) for dirpath in sp_path for fname in os.listdir(dirpath) if fname == 'aims.out'] + aims_out_path = sorted(lambda_path, key=self.sorting_key) + aims_out_path = [list(group) for key, group in groupby(aims_out_path, lambda x: x.split('/')[1])] + + if not os.path.exists('ext_xyz'): + os.mkdir('ext_xyz') + cnt = 0 + for numi, i in enumerate(aims_out_path): + filename = f"ext_xyz/ext_{i[0].split('/')[1]}_eigv.xyz" # + with open(filename, 'a') as f: # + for j in i: + ex = extractor() + ex.set_output_filepath(j) + #ex.set_scf_blocks + + + no_atoms = ex.get_no_atoms() + geometries, atom_label = ex.get_sp_geometries(j) + forces = ex.get_sp_forces(no_atoms, j) + total_energy = ex.get_sp_total_energy(j) + + ## subtract coulomb energy and force + print("Coulomb interactions are eliminated") + coulomb_E, coulomb_F = self.coulomb_E_F(atom_label, geometries) + total_energy = total_energy - coulomb_E + forces = forces - coulomb_F + form = np.concatenate((ex.get_sp_atom_order(), geometries, forces), axis=1) + + f.write(str(no_atoms) + '\n') + f.write(f'Lattice="0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0" Properties=species:S:1:pos:R:3:forces:R:3 energy={total_energy} pbc="F F F"\n') + np.savetxt(f, form, fmt="%s", delimiter=" ") + + + + def make_extxyz(self): + if not os.path.exists('FIT'): + os.mkdir('FIT') + else: pass + if os.path.exists('./FIT/Training_set.xyz'): + os.remove('./FIT/Training_set.xyz') + print("You may want to check the .xyz files in the FIT") + else: + with open('FIT/Training_set.xyz', 'a') as outfile: + filenames = [file for file in os.listdir('ext_xyz') if file.endswith('.xyz')] + sorted_filenames = sorted(filenames, key=lambda x: int(x.split('_')[1]) if x.split('_')[1] != 'breathing' else float('inf')) + for file in sorted_filenames: + print(file) + with open(os.path.join('ext_xyz', file), 'r') as infile: + for line in infile: + outfile.write(line) + + + def split_xyz_file(self, input_file, train_file, valid_file, test_file): + with open(input_file, 'r') as infile: + train_out = open(train_file, 'w') + valid_out = open(valid_file, 'w') + test_out = open(test_file, 'w') + + block_counter = 0 + line = infile.readline() + + while line: + if line.strip().isdigit(): + no_atoms = int(line.strip()) + block_lines = [line] + [infile.readline() for _ in range(no_atoms + 1)] # Read the block + + if block_counter % 5 < 3: + output_file = train_out + elif block_counter % 5 == 3: + output_file = valid_out + else: + output_file = test_out + + output_file.writelines(block_lines) + block_counter += 1 + + line = infile.readline() + + train_out.close() + valid_out.close() + test_out.close() + + + def coulomb_energy(self, r, cat_q, an_q): + return (cat_q * an_q) / r * 14.3996439067522 + + def coulomb_force(self, r, unit_r, cat_q, an_q): + return (cat_q * an_q) / r**2 * unit_r * 14.3996439067522 + + def coulomb_E_F(self, atom_label, structure): + self.charges = {"Al": 3, "F": -1} + coulomb_e = 0.0 + forces = np.zeros_like(structure) + + for i in range(len(structure)): + for j in range(i+1, len(structure)): + coord1 = structure[i] + coord2 = structure[j] + atom1 = atom_label[i][0] + atom2 = atom_label[j][0] + + r_vec = coord2 - coord1 + r = np.linalg.norm(r_vec) + unit_r = r_vec / r # unit vec + + # energy + energy_pair = self.coulomb_energy(r, self.charges[atom1], self.charges[atom2]) + coulomb_e += energy_pair + + # force + force_pair = self.coulomb_force(r, unit_r, self.charges[atom1], self.charges[atom2]) + + # add forces to atoms + forces[i] -= force_pair + forces[j] += force_pair + + return coulomb_e, forces + + + + + + +# executing the code using the class +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument("--eigenvector", type=str, help="A string of space-separated eigenvector indicies. For example, '7 8 9 10'") + parser.add_argument("--mode", type=str, choices=["run", "run_pair", "breath", "retrieve", "retrieve_c", "make_extxyz", "make_extxyz_"], default="run", help="Specify 'run' to execute the first part of the code, 'retrieve' to execute the second part of the code, or 'make_extxyz' to append all .xyz files into Training_set.xyz.") + args = parser.parse_args() + + ml = ML_train_generator() + step_size = ml.step_size ##### STEP SIZE ##### + + + # + # run + # + if args.mode == "run": + ml.initiate() + app_output = './aims.out' + + mod_sp = ml.mod_xyz_w_vib() # for each of vib. mode + sp_frame, no_atoms = ml.geometry_for_sp(mod_sp) + + shape = np.shape(sp_frame) + if args.eigenvector == 'all': + indicies = list(range(7, no_atoms*3+1)) + print("all eigenvectors without rotational and translational\n") + else: + indicies = list(map(int, args.eigenvector.split())) + + if not os.path.exists('sp'): + os.mkdir('sp') + else: pass + + for i in indicies: # Now we only iterate over the specified indicies + if not os.path.exists(os.path.join('sp', str(i))): + os.mkdir(f'sp/{str(i)}') + else: pass + + for numj, j in enumerate(np.arange(-1, 1+step_size, step_size)): + j = str(np.round(j, 2)) + os.mkdir(f'sp/{str(i)}/lambda_{j}') + + with open(f'sp/{i}/lambda_{j}/geometry.in', 'w') as f: + for row in sp_frame[i-1][numj]: + line = ' '.join(str(x) for x in row) + f.write(line + '\n') + ml.make_sp_control(f'sp/{i}/lambda_{j}') + ml.make_job_submit(f'sp/{i}/lambda_{j}') + os.chdir(f'sp/{i}/lambda_{j}') + os.system('qsub submit.sh') # submit jobs + os.chdir('../../../') + + # + # randomly pair up eigenvectors + # + elif args.mode == "run_pair": + ml.initiate() + app_output = './aims.out' + mod_sp, fname_pairs = ml.mod_xyz_w_rand_pair_vib() # randomly paired vib. mode + sp_frame, no_atoms = ml.geometry_for_sp(mod_sp) + shape = np.shape(sp_frame) + + if args.eigenvector == 'all': + #indicies = list(range(15)) + indicies = list(range(3*no_atoms))[6:] + print("all paired eigenvectors without E(3), (rotational and translational)\n") + else: + indicies = list(map(int, args.eigenvector.split())) + + if not os.path.exists('sp'): + os.mkdir('sp') + else: pass + + for i in range(int(len(indicies)/2)): # Now we only iterate over the specified indicies + #i = i+1 + if not os.path.exists(os.path.join('sp', fname_pairs[i])): + os.mkdir(f'sp/{fname_pairs[i]}_pair') + else: pass + + for numj, j in enumerate(np.arange(-1, 1+step_size, step_size)): + j = str(np.round(j, 2)) + os.mkdir(f'sp/{fname_pairs[i]}_pair/lambda_{j}') + with open(f'sp/{fname_pairs[i]}_pair/lambda_{j}/geometry.in', 'w') as f: + for row in sp_frame[i][numj]: + line = ' '.join(str(x) for x in row) + f.write(line + '\n') + ml.make_sp_control(f'sp/{fname_pairs[i]}_pair/lambda_{j}') + ml.make_job_submit(f'sp/{fname_pairs[i]}_pair/lambda_{j}') + os.chdir(f'sp/{fname_pairs[i]}_pair/lambda_{j}') + os.system('qsub submit.sh') # submit jobs + os.chdir('../../../') + + # + # preparen and run breathing mode single point calc + # + if args.mode == "breath": + ml.initiate() + if not os.path.exists('sp'): + os.mkdir('sp') + if not os.path.exists('sp/breathing'): + os.mkdir('sp/breathing') # for breathing mode + + mod_sp_breath, scale = ml.breathing() + #print(mod_sp_breath) + sp_frame, no_atoms = ml.geometry_for_sp(mod_sp_breath) + + # breathing + for numk, k in enumerate(scale): + k = str(np.round(k, 2)) + os.mkdir(f'sp/breathing/lambda_{k}') + + with open(f'sp/breathing/lambda_{k}/geometry.in', 'w') as f: + for row in sp_frame[numk]: + line = ' '.join(str(x) for x in row) + f.write(line + '\n') + ml.make_sp_control(f'sp/breathing/lambda_{k}') + ml.make_job_submit(f'sp/breathing/lambda_{k}') + os.chdir(f'sp/breathing/lambda_{k}') + os.system('qsub submit.sh') # submit job + os.chdir('../../../') + + # + # Collect data from single point calculated data + # + elif args.mode == "retrieve": + ml.initiate() + no_atoms = ml.no_atoms + if args.eigenvector == 'all': + indicies = list(range(7, no_atoms*3+1)) + indicies.append('breathing') + print(indicies) + print("all eigenvectors without rotational and translational\n") + else: + indicies = list(args.eigenvector.split()) + indicies = [int(x) if x.isdigit() else x for x in indicies] + print(indicies) + ml.retrieve_results(indicies) + + # collect coulomb subtracted data + elif args.mode == "retrieve_c": + ml.initiate() + no_atoms = ml.no_atoms + if args.eigenvector == 'all': + indicies = list(range(7, no_atoms*3+1)) + indicies.append('breathing') + print(indicies) + print("all eigenvectors without rotational and translational\n") + else: + indicies = list(args.eigenvector.split()) + indicies = [int(x) if x.isdigit() else x for x in indicies] + print(indicies) + ml.retrieve_results_c(indicies) + + # + # make training data and split to train, test, valid data + # + elif args.mode == "make_extxyz": + ml.make_extxyz() + print("splitting training, test, validation data") + ml.split_xyz_file('./FIT/Training_set.xyz', './FIT/Training_set_test.xyz', './FIT/Validation_set_test.xyz', './FIT/Testing_set_test.xyz') + + # dev + elif args.mode == "make_extxyz_": + ml.split_xyz_file('./Training_set.xyz', './Training_set_test.xyz', './Validation_set_test.xyz', './Testing_set_test.xyz') + + diff --git a/AppOutputExtractor/FHIaims/OutputPattern/FHIaims_22_patterns.json b/AppOutputExtractor/FHIaims/OutputPattern/FHIaims_22_patterns.json index e1068ae..dfcf5ee 100644 --- a/AppOutputExtractor/FHIaims/OutputPattern/FHIaims_22_patterns.json +++ b/AppOutputExtractor/FHIaims/OutputPattern/FHIaims_22_patterns.json @@ -1,54 +1,85 @@ { - "//compatibility": { "FHIAIMS22": "comment" }, - - "SUCCESS": { - "pattern": "'Have a nice day.'" - }, - "SCF_CONVERGED": { - "pattern": "'Self-consistency cycle converged.'", - "token": "None" - }, - "BEGIN_SCF": { - "pattern": "'Begin self-consistency iteration #'", - "token": "None" - }, - - - "SCF_ENERGY": { - "pattern": "'Total energy corrected :'", - "token": "6" - }, - "DIPOLE": { - "pattern": "'Absolute dipole moment'", - "token": "6" - }, - "DIPOLE_MOMENT": { - "pattern": "'Total dipole moment'", - "token_x": "7", "token_y": "8", "token_z": "9" - }, - "HOMO": { - "pattern": "'Highest occupied state'", - "token": "6" - }, - "LUMO": { - "pattern": "'Lowest unoccupied state'", - "token": "6" - }, - "HOMOLUMO": { - "pattern": "'Overall HOMO-LUMO gap'", - "token": "4" - }, - - - - - "APP_RUNTIME": { - "pattern": "'| Total time '", - "ctime_token": "5", - "wtime_token": "7" - }, - "APP_RESOURCE_USED": { - "pattern": "'parallel tasks.'", - "token": "2" - } + "//compatibility": { "FHIAIMS22": "comment" }, + + "SUCCESS": { + "pattern": "'Have a nice day.'" + }, + "SCF_CONVERGED": { + "pattern": "'Self-consistency cycle converged.'", + "token": "None" + }, + "BEGIN_SCF": { + "pattern": "'Begin self-consistency iteration #'", + "token": "None" + }, + + "NO_ATOMS": { + "pattern": "| Number of atoms", + "token": "6" + }, + + "SCF_ENERGY": { + "pattern": "'Total energy corrected :'", + "token": "6" + }, + + "SCF_GEOMETRY_BEGIN": { + "pattern": "'Updated atomic structure:'", + "token": "" + }, + + "SP": { + "pattern": "'Geometry relaxation not requested: no relaxation will be performed.'", + "token": "" + }, + + "SCF_GEOMETRY_END": { + "pattern": "'Final atomic structure:'", + "token": "" + }, + + + "SP_GEOMETRY_END": { + "pattern": "| Atomic structure:''", + "token": "" + }, + + "SCF_FORCE": { + "pattern":"'Total atomic forces '", + "token": "" + }, + + "DIPOLE": { + "pattern": "'Absolute dipole moment'", + "token": "6" + }, + "DIPOLE_MOMENT": { + "pattern": "'Total dipole moment'", + "token_x": "7", "token_y": "8", "token_z": "9" + }, + "HOMO": { + "pattern": "'Highest occupied state'", + "token": "6" + }, + "LUMO": { + "pattern": "'Lowest unoccupied state'", + "token": "6" + }, + "HOMOLUMO": { + "pattern": "'Overall HOMO-LUMO gap'", + "token": "4" + }, + + + + + "APP_RUNTIME": { + "pattern": "'| Total time '", + "ctime_token": "5", + "wtime_token": "7" + }, + "APP_RESOURCE_USED": { + "pattern": "'parallel tasks.'", + "token": "2" + } } diff --git a/AppOutputExtractor/FHIaims/__pycache__/FHIaimsMolecule.cpython-38.pyc b/AppOutputExtractor/FHIaims/__pycache__/FHIaimsMolecule.cpython-38.pyc new file mode 100644 index 0000000..07770ff Binary files /dev/null and b/AppOutputExtractor/FHIaims/__pycache__/FHIaimsMolecule.cpython-38.pyc differ diff --git a/AppOutputExtractor/FHIaims/__pycache__/FHIaimsMolecule.cpython-39.pyc b/AppOutputExtractor/FHIaims/__pycache__/FHIaimsMolecule.cpython-39.pyc new file mode 100644 index 0000000..929a4fc Binary files /dev/null and b/AppOutputExtractor/FHIaims/__pycache__/FHIaimsMolecule.cpython-39.pyc differ diff --git a/AppOutputExtractor/FHIaims/__pycache__/FHIaimsOutputExtractor.cpython-38.pyc b/AppOutputExtractor/FHIaims/__pycache__/FHIaimsOutputExtractor.cpython-38.pyc new file mode 100644 index 0000000..28dfda5 Binary files /dev/null and b/AppOutputExtractor/FHIaims/__pycache__/FHIaimsOutputExtractor.cpython-38.pyc differ diff --git a/AppOutputExtractor/FHIaims/__pycache__/FHIaimsOutputExtractor.cpython-39.pyc b/AppOutputExtractor/FHIaims/__pycache__/FHIaimsOutputExtractor.cpython-39.pyc new file mode 100644 index 0000000..a6cc185 Binary files /dev/null and b/AppOutputExtractor/FHIaims/__pycache__/FHIaimsOutputExtractor.cpython-39.pyc differ diff --git a/AppOutputExtractor/FHIaims/__pycache__/__init__.cpython-38.pyc b/AppOutputExtractor/FHIaims/__pycache__/__init__.cpython-38.pyc new file mode 100644 index 0000000..dbd4f58 Binary files /dev/null and b/AppOutputExtractor/FHIaims/__pycache__/__init__.cpython-38.pyc differ diff --git a/AppOutputExtractor/FHIaims/__pycache__/__init__.cpython-39.pyc b/AppOutputExtractor/FHIaims/__pycache__/__init__.cpython-39.pyc new file mode 100644 index 0000000..5f47396 Binary files /dev/null and b/AppOutputExtractor/FHIaims/__pycache__/__init__.cpython-39.pyc differ diff --git a/AppOutputExtractor/FHIaims/firstblock.txt b/AppOutputExtractor/FHIaims/firstblock.txt deleted file mode 100644 index de806f1..0000000 --- a/AppOutputExtractor/FHIaims/firstblock.txt +++ /dev/null @@ -1,424 +0,0 @@ - Begin self-consistency iteration # 15 - - Date : 20230426, Time : 001902.578 ------------------------------------------------------------- - Pulay mixing of updated and previous charge densities. - Renormalizing the density to the exact electron count on the 3D integration grid. - | Formal number of electrons (from input files) : 180.0000000000 - | Integrated number of electrons on 3D grid : 180.0000000000 - | Charge integration error : -0.0000000000 - | Normalization factor for density and gradient : 1.0000000000 - - Evaluating partitioned Hartree potential by multipole expansion. - | Original multipole sum: apparent total charge = 0.428591E-12 - | Sum of charges compensated after spline to logarithmic grids = 0.142513E-05 - | Analytical far-field extrapolation by fixed multipoles: - | Hartree multipole sum: apparent total charge = 0.427607E-12 - Summing up the Hartree potential. - Time summed over all CPUs for potential: real work 1.124 s, elapsed 1.163 s - | RMS charge density error from multipole expansion : 0.351565E-02 - - Integrating Hamiltonian matrix: batch-based integration. - Time summed over all CPUs for integration: real work 1.671 s, elapsed 1.888 s - - Updating Kohn-Sham eigenvalues and eigenvectors using ELSI and the ELPA eigensolver. - Starting ELPA eigensolver - Finished transformation to standard eigenproblem - | Time : 0.001 s - Finished solving standard eigenproblem - | Time : 0.003 s - Finished back-transformation of eigenvectors - | Time : 0.000 s - - Obtaining occupation numbers and electronic chemical potential using ELSI. - | Note that, for insulating systems, the printed 'chemical potential' value is not uniquely defined. - | It can be anywhere in the energy gap, as long as it correctly separates occupied and unoccupied states. - | In systems with a gap, the physically relevant chemical potential is the VBM or HOMO. - - | Chemical potential (Fermi level): -3.55139282 eV - Highest occupied state (VBM) at -5.40798974 eV - | Occupation number: 2.00000000 - - Lowest unoccupied state (CBM) at -2.81995551 eV - | Occupation number: 0.00000000 - - Overall HOMO-LUMO gap: 2.58803422 eV. - - Total energy components: - | Sum of eigenvalues : -27553.01557732 Ha -749755.70100900 eV - | XC energy correction : -803.04599761 Ha -21851.99341228 eV - | XC potential correction : 1050.09861755 Ha 28574.63724494 eV - | Free-atom electrostatic energy: -16214.56575149 Ha -441220.78316434 eV - | Hartree energy correction : -5.85013205 Ha -159.19019262 eV - | Entropy correction : 0.00000000 Ha 0.00000000 eV - | --------------------------- - | Total energy : -43526.37884093 Ha -1184413.03053330 eV - | Total energy, T -> 0 : -43526.37884093 Ha -1184413.03053330 eV <-- do not rely on this value for anything but (periodic) metals - | Electronic free energy : -43526.37884093 Ha -1184413.03053330 eV - - Derived energy quantities: - | Kinetic energy : 49353.68881641 Ha 1342982.20287661 eV - | Electrostatic energy : -92077.02165973 Ha -2505543.23999763 eV - | Energy correction for multipole - | error in Hartree potential : 0.00685185 Ha 0.18644836 eV - | Sum of eigenvalues per atom : -187438.92525225 eV - | Total energy (T->0) per atom : -296103.25763333 eV <-- do not rely on this value for anything but (periodic) metals - | Electronic free energy per atom : -296103.25763333 eV - Evaluating new KS density and force components. - Integration grid: deviation in total charge ( - N_e) = 4.263256E-13 - - atomic forces [eV/Ang]: - ----------------------- - atom # 1 - Hellmann-Feynman : -0.427134E+01 -0.430813E+01 0.282176E-01 - Ionic forces : 0.000000E+00 0.000000E+00 0.000000E+00 - Multipole : 0.353165E-01 0.357636E-01 -0.239781E-03 - Hartree pot. SCF incomplete : -0.608032E-04 -0.609048E-04 0.136866E-06 - Pulay + GGA : 0.427693E+01 0.431786E+01 -0.326554E-01 - ---------------------------------------------------------------- - Total forces( 1) : 0.408444E-01 0.454293E-01 -0.467746E-02 - atom # 2 - Hellmann-Feynman : -0.477720E+02 0.483083E+02 0.133806E+00 - Ionic forces : 0.000000E+00 0.000000E+00 0.000000E+00 - Multipole : 0.193617E+00 -0.209143E+00 -0.208047E-02 - Hartree pot. SCF incomplete : 0.108976E-03 -0.107322E-03 -0.786876E-07 - Pulay + GGA : 0.474587E+02 -0.481052E+02 -0.127628E+00 - ---------------------------------------------------------------- - Total forces( 2) : -0.119620E+00 -0.608182E-02 0.409783E-02 - atom # 3 - Hellmann-Feynman : 0.434191E+01 0.440368E+01 0.627914E-01 - Ionic forces : 0.000000E+00 0.000000E+00 0.000000E+00 - Multipole : -0.358109E-01 -0.361969E-01 -0.519479E-03 - Hartree pot. SCF incomplete : 0.596425E-04 0.607448E-04 0.614732E-06 - Pulay + GGA : -0.424105E+01 -0.428178E+01 -0.659270E-01 - ---------------------------------------------------------------- - Total forces( 3) : 0.651007E-01 0.857600E-01 -0.365440E-02 - atom # 4 - Hellmann-Feynman : 0.483864E+02 -0.475926E+02 0.157051E+00 - Ionic forces : 0.000000E+00 0.000000E+00 0.000000E+00 - Multipole : -0.209188E+00 0.195108E+00 -0.216362E-02 - Hartree pot. SCF incomplete : -0.108428E-03 0.108169E-03 -0.137507E-06 - Pulay + GGA : -0.481631E+02 0.472726E+02 -0.150675E+00 - ---------------------------------------------------------------- - Total forces( 4) : 0.139753E-01 -0.124847E+00 0.421225E-02 - - - Self-consistency convergence accuracy: - | Change of charge density : 0.4665E-06 - | Change of unmixed KS density : 0.1721E-06 - | Change of sum of eigenvalues : -0.3754E-03 eV - | Change of total energy : 0.1584E-08 eV - | Change of forces : 0.1661E+01 eV/A - - Writing Kohn-Sham eigenvalues. - - State Occupation Eigenvalue [Ha] Eigenvalue [eV] - 1 2.00000 -3583.266931 -97505.65422 - 2 2.00000 -3583.266752 -97505.64936 - 3 2.00000 -591.276714 -16089.45800 - 4 2.00000 -591.276537 -16089.45321 - 5 2.00000 -505.933039 -13767.13846 - 6 2.00000 -505.932863 -13767.13367 - 7 2.00000 -505.932546 -13767.12503 - 8 2.00000 -505.932369 -13767.12022 - 9 2.00000 -505.932251 -13767.11702 - 10 2.00000 -505.932075 -13767.11222 - 11 2.00000 -139.804775 -3804.28148 - 12 2.00000 -139.804602 -3804.27677 - 13 2.00000 -115.928567 -3154.57682 - 14 2.00000 -115.928395 -3154.57214 - 15 2.00000 -115.927537 -3154.54878 - 16 2.00000 -115.927363 -3154.54405 - 17 2.00000 -115.926698 -3154.52596 - 18 2.00000 -115.926525 -3154.52124 - 19 2.00000 -91.712908 -2495.63520 - 20 2.00000 -91.712792 -2495.63205 - 21 2.00000 -91.712735 -2495.63050 - 22 2.00000 -91.712619 -2495.62733 - 23 2.00000 -91.712237 -2495.61695 - 24 2.00000 -91.712064 -2495.61224 - 25 2.00000 -91.711535 -2495.59784 - 26 2.00000 -91.711485 -2495.59648 - 27 2.00000 -91.711361 -2495.59311 - 28 2.00000 -91.711311 -2495.59174 - 29 2.00000 -31.488813 -856.85421 - 30 2.00000 -31.488644 -856.84960 - 31 2.00000 -23.725402 -645.60104 - 32 2.00000 -23.725235 -645.59650 - 33 2.00000 -23.722667 -645.52662 - 34 2.00000 -23.722497 -645.52198 - 35 2.00000 -23.720931 -645.47936 - 36 2.00000 -23.720762 -645.47476 - 37 2.00000 -18.724882 -509.52997 - 38 2.00000 -18.724336 -509.51511 - 39 2.00000 -14.766266 -401.81054 - 40 2.00000 -14.766098 -401.80598 - 41 2.00000 -14.765869 -401.79974 - 42 2.00000 -14.765700 -401.79515 - 43 2.00000 -14.764575 -401.76451 - 44 2.00000 -14.764407 -401.75996 - 45 2.00000 -14.762413 -401.70569 - 46 2.00000 -14.762384 -401.70492 - 47 2.00000 -14.762243 -401.70107 - 48 2.00000 -14.762215 -401.70030 - 49 2.00000 -5.278065 -143.62345 - 50 2.00000 -5.277904 -143.61907 - 51 2.00000 -4.879229 -132.77058 - 52 2.00000 -4.879137 -132.76808 - 53 2.00000 -4.879062 -132.76602 - 54 2.00000 -4.878969 -132.76351 - 55 2.00000 -4.877889 -132.73411 - 56 2.00000 -4.877722 -132.72958 - 57 2.00000 -4.877071 -132.71185 - 58 2.00000 -4.876903 -132.70728 - 59 2.00000 -4.876853 -132.70592 - 60 2.00000 -4.876686 -132.70137 - 61 2.00000 -4.874795 -132.64993 - 62 2.00000 -4.874795 -132.64992 - 63 2.00000 -4.874626 -132.64532 - 64 2.00000 -4.874625 -132.64530 - 65 2.00000 -3.143615 -85.54213 - 66 2.00000 -3.143459 -85.53787 - 67 2.00000 -3.135137 -85.31143 - 68 2.00000 -3.134974 -85.30699 - 69 2.00000 -3.131508 -85.21267 - 70 2.00000 -3.131306 -85.20717 - 71 2.00000 -0.825678 -22.46785 - 72 2.00000 -0.822183 -22.37274 - 73 2.00000 -0.774831 -21.08421 - 74 2.00000 -0.772088 -21.00959 - 75 2.00000 -0.771301 -20.98818 - 76 2.00000 -0.771094 -20.98254 - 77 2.00000 -0.767042 -20.87227 - 78 2.00000 -0.766868 -20.86753 - 79 2.00000 -0.765256 -20.82366 - 80 2.00000 -0.757438 -20.61094 - 81 2.00000 -0.745651 -20.29019 - 82 2.00000 -0.695996 -18.93900 - 83 2.00000 -0.448978 -12.21730 - 84 2.00000 -0.411528 -11.19824 - 85 2.00000 -0.267517 -7.27951 - 86 2.00000 -0.247519 -6.73533 - 87 2.00000 -0.246598 -6.71027 - 88 2.00000 -0.208022 -5.66057 - 89 2.00000 -0.205206 -5.58393 - 90 2.00000 -0.198740 -5.40799 - 91 0.00000 -0.103631 -2.81996 - 92 0.00000 -0.083839 -2.28139 - 93 0.00000 -0.057101 -1.55380 - 94 0.00000 -0.008598 -0.23396 - 95 0.00000 0.002333 0.06349 - 96 0.00000 0.014782 0.40225 - 97 0.00000 0.118000 3.21093 - 98 0.00000 0.187317 5.09715 - 99 0.00000 0.189417 5.15430 - 100 0.00000 0.205704 5.59750 - 101 0.00000 0.218198 5.93747 - 102 0.00000 0.237262 6.45624 - 103 0.00000 0.246051 6.69539 - 104 0.00000 0.247602 6.73760 - 105 0.00000 0.252724 6.87696 - 106 0.00000 0.263359 7.16635 - 107 0.00000 0.275904 7.50772 - 108 0.00000 0.292204 7.95128 - 109 0.00000 0.339497 9.23818 - 110 0.00000 0.360949 9.82193 - 111 0.00000 0.382825 10.41719 - 112 0.00000 0.412553 11.22615 - 113 0.00000 0.486200 13.23018 - 114 0.00000 0.525701 14.30505 - - Highest occupied state (VBM) at -5.40798974 eV - | Occupation number: 2.00000000 - - Lowest unoccupied state (CBM) at -2.81995551 eV - | Occupation number: 0.00000000 - - Overall HOMO-LUMO gap: 2.58803422 eV. - - | Chemical Potential : -3.55139282 eV - | Note that, for insulating systems, the printed 'chemical potential' value is not uniquely defined. - | It can be anywhere in the energy gap, as long as it correctly separates occupied and unoccupied states. - | In systems with a gap, the physically relevant chemical potential is the VBM or HOMO. - - Self-consistency cycle converged. - - ------------------------------------------------------------- - End self-consistency iteration # 15 : max(cpu_time) wall_clock(cpu1) - | Time for this iteration : 2.934 s 5.277 s - | Charge density & force component update : 1.998 s 3.265 s - | Density mixing : 0.017 s 0.030 s - | Hartree multipole update : 0.008 s 0.016 s - | Hartree multipole summation : 0.304 s 0.591 s - | Hartree pot. SCF incomplete forces : 0.211 s 0.407 s - | Integration : 0.517 s 0.961 s - | Solution of K.-S. eqns. : 0.005 s 0.005 s - | Total energy evaluation : 0.001 s 0.000 s - - Partial memory accounting: - | Current value for overall tracked memory usage: - | Minimum: 0.580 MB (on task 0) - | Maximum: 0.580 MB (on task 0) - | Average: 0.580 MB - | Peak value for overall tracked memory usage: - | Minimum: 2.738 MB (on task 0 after allocating grid_partition) - | Maximum: 2.742 MB (on task 1 after allocating grid_partition) - | Average: 2.740 MB - | Largest tracked array allocation so far: - | Minimum: 1.625 MB (all_coords on task 0) - | Maximum: 1.628 MB (all_coords on task 1) - | Average: 1.627 MB - Note: These values currently only include a subset of arrays which are explicitly tracked. - The "true" memory usage will be greater. ------------------------------------------------------------- - ------------------------------------------------------------ - Computing monopole / dipole moments - | Total electronic charge [e] : 0.180000000000000E+03 - | Total ionic charge [e] : 0.180000000000000E+03 - | Total charge [e] : 0.000000000000000E+00 - | Total dipole moment [eAng] : 0.830045719734522E-02 0.106112352433544E-01 -0.592209490176695E-02 - | Absolute dipole moment : 0.147162193208820E-01 eAng / 0.706739976625820E-01 Debye . - Removing unitary transformations (pure translations, rotations) from forces on atoms. - Atomic forces before filtering: - | Net force on center of mass : 0.300206E-03 0.260269E-03 -0.217901E-04 eV/A - | Net torque on center of mass: -0.308482E-04 -0.330117E-04 0.102555E-03 eV - Atomic forces after filtering: - | Net force on center of mass : 0.836279E-17 -0.223008E-16 0.139380E-17 eV/A - | Net torque on center of mass: -0.368783E-18 0.221270E-17 0.236021E-16 eV - - Energy and forces in a compact form: - | Total energy uncorrected : -0.118441303053330E+07 eV - | Total energy corrected : -0.118441303053330E+07 eV <-- do not rely on this value for anything but (periodic) metals - | Electronic free energy : -0.118441303053330E+07 eV - Total atomic forces (unitary forces cleaned) [eV/Ang]: - | 1 0.407581257835650E-01 0.453754197891383E-01 -0.467153716479163E-02 - | 2 -0.119707868441570E+00 -0.615961708595437E-02 0.408928843029504E-02 - | 3 0.650368938195568E-01 0.856837919148820E-01 -0.364943747681396E-02 - | 4 0.139128488384485E-01 -0.124899594618066E+00 0.423168621131056E-02 - - ------------------------------------ - Start decomposition of the XC Energy - ------------------------------------ - X and C from original XC functional choice - Hartree-Fock Energy : 0.000000000 Ha 0.000000000 eV - X Energy : -792.435935834 Ha -21563.278941591 eV - C Energy : -10.610061781 Ha -288.714470688 eV - XC Energy w/o HF : -803.045997615 Ha -21851.993412279 eV - Total XC Energy : -803.045997615 Ha -21851.993412279 eV - ------------------------------------ - LDA X and C from self-consistent density - X Energy LDA : -759.825825678 Ha -20675.912695541 eV - C Energy LDA : -18.362895984 Ha -499.679823159 eV - ------------------------------------ - End decomposition of the XC Energy - ------------------------------------ - ------------------------------------------------------------- - Relaxation / MD: End force evaluation. : max(cpu_time) wall_clock(cpu1) - | Time for this force evaluation : 16.103 s 32.525 s - ------------------------------------------------------------- - Geometry optimization: Attempting to predict improved coordinates. - - Removing unitary transformations (pure translations, rotations) from forces on atoms. - Atomic forces before filtering: - | Net force on center of mass : 0.836279E-17 -0.223008E-16 0.139380E-17 eV/A - | Net torque on center of mass: -0.368783E-18 0.221270E-17 0.236021E-16 eV - Atomic forces after filtering: - | Net force on center of mass : -0.139380E-16 -0.223008E-16 0.139380E-17 eV/A - | Net torque on center of mass: -0.368783E-18 0.221270E-17 0.000000E+00 eV - Net remaining forces (excluding translations, rotations) in present geometry: - || Forces on atoms || = 0.124900E+00 eV/A. - Maximum force component is 0.124900E+00 eV/A. - Present geometry is not yet converged. - - Relaxation step number 1: Predicting new coordinates. - - Advancing geometry using trust radius method. - Allocating 0.098 MB for eigenvec_stored - | Hessian has 0 negative and 6 zero eigenvalues. - | Positive eigenvalues (eV/A^2): 2.01E+00 ... 2.75E+01 - | Use Quasi-Newton step of length |H^-1 F| = 4.19E-02 A. - Finished advancing geometry - | Time : 0.001 s - Updated atomic structure: - x [A] y [A] z [A] - atom -0.01452809 -0.02249747 -0.00488966 O - atom 2.12598529 -0.12169905 -0.00447988 Pb - atom 2.01126266 2.01794109 0.00297922 O - atom -0.12954139 2.11661581 -0.00505601 Pb ------------------------------------------------------------- - Writing the current geometry to file "geometry.in.next_step". - Writing estimated Hessian matrix to file 'hessian.aims' - ------------------------------------------------------------- - Begin self-consistency loop: Re-initialization. - - Date : 20230426, Time : 001908.086 ------------------------------------------------------------- - - Initializing index lists of integration centers etc. from given atomic structure: - | Number of centers in hartree potential : 4 - | Number of centers in hartree multipole : 4 - | Number of centers in electron density summation: 4 - | Number of centers in basis integrals : 4 - | Number of centers in integrals : 4 - | Number of centers in hamiltonian : 4 - Hamiltonian matrix size: - | Size of matrix non-packed: 12880 - | Size of matrix packed: 12392 - Partitioning the integration grid into batches with parallel hashing+maxmin method. - | Number of batches: 1885 - | Maximal batch size: 141 - | Minimal batch size: 67 - | Average batch size: 75.412 - | Standard deviation of batch sizes: 19.544 - - Integration load balanced across 2 MPI tasks. - Work distribution over tasks is as follows: - Initializing partition tables, free-atom densities, potentials, etc. across the integration grid (initialize_grid_storage). - | Species 1: outer_partition_radius set to 6.064445521808699 AA . - | Species 2: outer_partition_radius set to 7.333634274623765 AA . - | Species 3: outer_partition_radius set to 6.071521564059752 AA . - | Species 4: outer_partition_radius set to 7.541002451619759 AA . - | Species 5: outer_partition_radius set to 6.037031146364955 AA . - | The sparse table of interatomic distances needs 0.16 kbyte instead of 0.13 kbyte of memory. - | Using the partition_type stratmann_smoother will reduce your memory usage. - | Net number of integration points: 142152 - | of which are non-zero points : 134924 - Renormalizing the initial density to the exact electron count on the 3D integration grid. - | Initial density: Formal number of electrons (from input files) : 180.0000000000 - | Integrated number of electrons on 3D grid : 179.9999998583 - | Charge integration error : -0.0000001417 - | Normalization factor for density and gradient : 1.0000000008 - Renormalizing the free-atom superposition density to the exact electron count on the 3D integration grid. - | Formal number of electrons (from input files) : 180.0000000000 - | Integrated number of electrons on 3D grid : 179.9999998583 - | Charge integration error : -0.0000001417 - | Normalization factor for density and gradient : 1.0000000008 - Obtaining max. number of non-zero basis functions in each batch (get_n_compute_maxes). - Calculating total energy contributions from superposition of free atom densities. - Initialize hartree_potential_storage - Integrating overlap matrix. - Time summed over all CPUs for integration: real work 0.621 s, elapsed 0.634 s - Normalizing ScaLAPACK eigenvectors - Finished Gram-Schmidt orthonormalization - | Time : 0.018 s - - End scf reinitialization - timings : max(cpu_time) wall_clock(cpu1) - | Time for scf. reinitialization : 0.504 s 1.010 s - | Boundary condition initialization : 0.000 s 0.001 s - | Integration : 0.187 s 0.369 s - | Grid partitioning : 0.132 s 0.258 s - | Preloading free-atom quantities on grid : 0.119 s 0.241 s - | Free-atom superposition energy : 0.061 s 0.120 s - | K.-S. eigenvector reorthonormalization : 0.009 s 0.019 s ------------------------------------------------------------- - Evaluating new KS density. - Integration grid: deviation in total charge ( - N_e) = -3.410605E-13 - - Time for density update prior : max(cpu_time) wall_clock(cpu1) - | self-consistency iterative process : 0.381 s 0.749 s - ------------------------------------------------------------- diff --git a/AppOutputExtractor/FHIaims/geoA.txt b/AppOutputExtractor/FHIaims/geoA.txt deleted file mode 100644 index adde186..0000000 --- a/AppOutputExtractor/FHIaims/geoA.txt +++ /dev/null @@ -1 +0,0 @@ -atom 0 0 1 Pb diff --git a/AppOutputExtractor/FHIaims/geoB.txt b/AppOutputExtractor/FHIaims/geoB.txt deleted file mode 100644 index dac2437..0000000 --- a/AppOutputExtractor/FHIaims/geoB.txt +++ /dev/null @@ -1 +0,0 @@ -atom 0 0 0 Pb diff --git a/AppOutputExtractor/FHIaims/lastblock.txt b/AppOutputExtractor/FHIaims/lastblock.txt deleted file mode 100644 index f3811bd..0000000 --- a/AppOutputExtractor/FHIaims/lastblock.txt +++ /dev/null @@ -1,449 +0,0 @@ - Begin self-consistency iteration # 9 - - Date : 20230426, Time : 002042.420 ------------------------------------------------------------- - Pulay mixing of updated and previous charge densities. - Renormalizing the density to the exact electron count on the 3D integration grid. - | Formal number of electrons (from input files) : 180.0000000000 - | Integrated number of electrons on 3D grid : 180.0000000000 - | Charge integration error : -0.0000000000 - | Normalization factor for density and gradient : 1.0000000000 - - Evaluating partitioned Hartree potential by multipole expansion. - | Original multipole sum: apparent total charge = 0.261720E-13 - | Sum of charges compensated after spline to logarithmic grids = 0.142393E-05 - | Analytical far-field extrapolation by fixed multipoles: - | Hartree multipole sum: apparent total charge = 0.253357E-13 - Summing up the Hartree potential. - Time summed over all CPUs for potential: real work 1.140 s, elapsed 1.145 s - | RMS charge density error from multipole expansion : 0.351863E-02 - - Integrating Hamiltonian matrix: batch-based integration. - Time summed over all CPUs for integration: real work 1.291 s, elapsed 1.498 s - - Updating Kohn-Sham eigenvalues and eigenvectors using ELSI and the ELPA eigensolver. - Starting ELPA eigensolver - Finished transformation to standard eigenproblem - | Time : 0.000 s - Finished solving standard eigenproblem - | Time : 0.015 s - Finished back-transformation of eigenvectors - | Time : 0.000 s - - Obtaining occupation numbers and electronic chemical potential using ELSI. - | Note that, for insulating systems, the printed 'chemical potential' value is not uniquely defined. - | It can be anywhere in the energy gap, as long as it correctly separates occupied and unoccupied states. - | In systems with a gap, the physically relevant chemical potential is the VBM or HOMO. - - | Chemical potential (Fermi level): -3.56862396 eV - Highest occupied state (VBM) at -5.44468173 eV - | Occupation number: 2.00000000 - - Lowest unoccupied state (CBM) at -2.81614171 eV - | Occupation number: 0.00000000 - - Overall HOMO-LUMO gap: 2.62854001 eV. - - Total energy components: - | Sum of eigenvalues : -27553.03389257 Ha -749756.19939224 eV - | XC energy correction : -803.04731297 Ha -21852.02920497 eV - | XC potential correction : 1050.10043214 Ha 28574.68662260 eV - | Free-atom electrostatic energy: -16214.55472269 Ha -441220.48305533 eV - | Hartree energy correction : -5.84339317 Ha -159.00681835 eV - | Entropy correction : 0.00000000 Ha 0.00000000 eV - | --------------------------- - | Total energy : -43526.37888926 Ha -1184413.03184828 eV - | Total energy, T -> 0 : -43526.37888926 Ha -1184413.03184828 eV <-- do not rely on this value for anything but (periodic) metals - | Electronic free energy : -43526.37888926 Ha -1184413.03184828 eV - - Derived energy quantities: - | Kinetic energy : 49353.70099579 Ha 1342982.53429447 eV - | Electrostatic energy : -92077.03257208 Ha -2505543.53693778 eV - | Energy correction for multipole - | error in Hartree potential : 0.00676088 Ha 0.18397294 eV - | Sum of eigenvalues per atom : -187439.04984806 eV - | Total energy (T->0) per atom : -296103.25796207 eV <-- do not rely on this value for anything but (periodic) metals - | Electronic free energy per atom : -296103.25796207 eV - Evaluating new KS density and force components. - Integration grid: deviation in total charge ( - N_e) = -6.536993E-13 - - atomic forces [eV/Ang]: - ----------------------- - atom # 1 - Hellmann-Feynman : -0.432587E+01 -0.435919E+01 -0.133006E-01 - Ionic forces : 0.000000E+00 0.000000E+00 0.000000E+00 - Multipole : 0.338709E-01 0.341236E-01 0.102694E-03 - Hartree pot. SCF incomplete : 0.908767E-06 -0.102766E-06 0.102549E-05 - Pulay + GGA : 0.429244E+01 0.432511E+01 0.128747E-01 - ---------------------------------------------------------------- - Total forces( 1) : 0.441586E-03 0.500880E-04 -0.322210E-03 - atom # 2 - Hellmann-Feynman : -0.475920E+02 0.472541E+02 -0.124252E-02 - Ionic forces : 0.000000E+00 0.000000E+00 0.000000E+00 - Multipole : 0.196036E+00 -0.194719E+00 -0.108403E-03 - Hartree pot. SCF incomplete : 0.310265E-06 -0.579718E-06 0.218626E-05 - Pulay + GGA : 0.473958E+02 -0.470600E+02 0.167459E-02 - ---------------------------------------------------------------- - Total forces( 2) : -0.240986E-03 -0.569363E-03 0.325850E-03 - atom # 3 - Hellmann-Feynman : 0.432610E+01 0.435955E+01 0.203048E-01 - Ionic forces : 0.000000E+00 0.000000E+00 0.000000E+00 - Multipole : -0.338725E-01 -0.341249E-01 -0.161459E-03 - Hartree pot. SCF incomplete : -0.376153E-06 0.531878E-06 0.102378E-05 - Pulay + GGA : -0.429230E+01 -0.432497E+01 -0.204666E-01 - ---------------------------------------------------------------- - Total forces( 3) : -0.707873E-04 0.454105E-03 -0.322278E-03 - atom # 4 - Hellmann-Feynman : 0.475942E+02 -0.472515E+02 0.229773E-01 - Ionic forces : 0.000000E+00 0.000000E+00 0.000000E+00 - Multipole : -0.196094E+00 0.194670E+00 -0.207324E-03 - Hartree pot. SCF incomplete : -0.491079E-06 0.782087E-06 0.218462E-05 - Pulay + GGA : -0.473983E+02 0.470569E+02 -0.224597E-01 - ---------------------------------------------------------------- - Total forces( 4) : -0.130584E-03 0.644953E-04 0.312443E-03 - - - Self-consistency convergence accuracy: - | Change of charge density : 0.1024E-06 - | Change of unmixed KS density : 0.1180E-06 - | Change of sum of eigenvalues : 0.1077E-03 eV - | Change of total energy : 0.0000E+00 eV - | Change of forces : 0.1670E+01 eV/A - - Writing Kohn-Sham eigenvalues. - - State Occupation Eigenvalue [Ha] Eigenvalue [eV] - 1 2.00000 -3583.266836 -97505.65163 - 2 2.00000 -3583.266835 -97505.65161 - 3 2.00000 -591.276628 -16089.45567 - 4 2.00000 -591.276627 -16089.45564 - 5 2.00000 -505.932953 -13767.13613 - 6 2.00000 -505.932953 -13767.13610 - 7 2.00000 -505.932456 -13767.12258 - 8 2.00000 -505.932455 -13767.12255 - 9 2.00000 -505.932168 -13767.11474 - 10 2.00000 -505.932167 -13767.11472 - 11 2.00000 -139.804701 -3804.27947 - 12 2.00000 -139.804700 -3804.27944 - 13 2.00000 -115.928494 -3154.57483 - 14 2.00000 -115.928493 -3154.57481 - 15 2.00000 -115.927456 -3154.54657 - 16 2.00000 -115.927455 -3154.54654 - 17 2.00000 -115.926631 -3154.52413 - 18 2.00000 -115.926630 -3154.52411 - 19 2.00000 -91.712831 -2495.63311 - 20 2.00000 -91.712830 -2495.63308 - 21 2.00000 -91.712713 -2495.62989 - 22 2.00000 -91.712712 -2495.62986 - 23 2.00000 -91.712167 -2495.61504 - 24 2.00000 -91.712166 -2495.61501 - 25 2.00000 -91.711460 -2495.59579 - 26 2.00000 -91.711459 -2495.59576 - 27 2.00000 -91.711412 -2495.59449 - 28 2.00000 -91.711411 -2495.59446 - 29 2.00000 -31.488755 -856.85263 - 30 2.00000 -31.488754 -856.85260 - 31 2.00000 -23.725347 -645.59954 - 32 2.00000 -23.725346 -645.59952 - 33 2.00000 -23.722590 -645.52453 - 34 2.00000 -23.722589 -645.52450 - 35 2.00000 -23.720894 -645.47836 - 36 2.00000 -23.720893 -645.47833 - 37 2.00000 -18.725199 -509.53858 - 38 2.00000 -18.725197 -509.53853 - 39 2.00000 -14.766205 -401.80887 - 40 2.00000 -14.766204 -401.80885 - 41 2.00000 -14.765800 -401.79785 - 42 2.00000 -14.765799 -401.79782 - 43 2.00000 -14.764535 -401.76345 - 44 2.00000 -14.764534 -401.76342 - 45 2.00000 -14.762357 -401.70417 - 46 2.00000 -14.762356 -401.70414 - 47 2.00000 -14.762335 -401.70359 - 48 2.00000 -14.762334 -401.70356 - 49 2.00000 -5.278042 -143.62284 - 50 2.00000 -5.278040 -143.62278 - 51 2.00000 -4.879163 -132.76878 - 52 2.00000 -4.879162 -132.76876 - 53 2.00000 -4.879068 -132.76618 - 54 2.00000 -4.879067 -132.76616 - 55 2.00000 -4.877849 -132.73303 - 56 2.00000 -4.877848 -132.73301 - 57 2.00000 -4.877018 -132.71041 - 58 2.00000 -4.877017 -132.71039 - 59 2.00000 -4.876814 -132.70485 - 60 2.00000 -4.876813 -132.70483 - 61 2.00000 -4.874744 -132.64853 - 62 2.00000 -4.874743 -132.64850 - 63 2.00000 -4.874743 -132.64850 - 64 2.00000 -4.874742 -132.64847 - 65 2.00000 -3.143604 -85.54182 - 66 2.00000 -3.143596 -85.54161 - 67 2.00000 -3.135060 -85.30932 - 68 2.00000 -3.135048 -85.30901 - 69 2.00000 -3.131614 -85.21555 - 70 2.00000 -3.131478 -85.21186 - 71 2.00000 -0.825800 -22.47117 - 72 2.00000 -0.823418 -22.40633 - 73 2.00000 -0.774882 -21.08561 - 74 2.00000 -0.772117 -21.01038 - 75 2.00000 -0.771329 -20.98892 - 76 2.00000 -0.771324 -20.98879 - 77 2.00000 -0.767155 -20.87536 - 78 2.00000 -0.766961 -20.87007 - 79 2.00000 -0.765305 -20.82501 - 80 2.00000 -0.757471 -20.61183 - 81 2.00000 -0.746069 -20.30158 - 82 2.00000 -0.696268 -18.94642 - 83 2.00000 -0.448759 -12.21134 - 84 2.00000 -0.412199 -11.21651 - 85 2.00000 -0.267997 -7.29256 - 86 2.00000 -0.248056 -6.74994 - 87 2.00000 -0.246873 -6.71776 - 88 2.00000 -0.208774 -5.68102 - 89 2.00000 -0.205011 -5.57862 - 90 2.00000 -0.200088 -5.44468 - 91 0.00000 -0.103491 -2.81614 - 92 0.00000 -0.084161 -2.29015 - 93 0.00000 -0.056578 -1.53956 - 94 0.00000 -0.009237 -0.25136 - 95 0.00000 0.002869 0.07808 - 96 0.00000 0.014594 0.39712 - 97 0.00000 0.117962 3.20991 - 98 0.00000 0.187064 5.09026 - 99 0.00000 0.189739 5.16305 - 100 0.00000 0.205703 5.59747 - 101 0.00000 0.218397 5.94288 - 102 0.00000 0.237274 6.45654 - 103 0.00000 0.245847 6.68983 - 104 0.00000 0.246602 6.71039 - 105 0.00000 0.252011 6.85757 - 106 0.00000 0.262894 7.15370 - 107 0.00000 0.275788 7.50457 - 108 0.00000 0.292493 7.95914 - 109 0.00000 0.340114 9.25496 - 110 0.00000 0.361529 9.83770 - 111 0.00000 0.383527 10.43629 - 112 0.00000 0.413601 11.25467 - 113 0.00000 0.487286 13.25971 - 114 0.00000 0.525067 14.28781 - - Highest occupied state (VBM) at -5.44468173 eV - | Occupation number: 2.00000000 - - Lowest unoccupied state (CBM) at -2.81614171 eV - | Occupation number: 0.00000000 - - Overall HOMO-LUMO gap: 2.62854001 eV. - - | Chemical Potential : -3.56862396 eV - | Note that, for insulating systems, the printed 'chemical potential' value is not uniquely defined. - | It can be anywhere in the energy gap, as long as it correctly separates occupied and unoccupied states. - | In systems with a gap, the physically relevant chemical potential is the VBM or HOMO. - - Self-consistency cycle converged. - - ------------------------------------------------------------- - End self-consistency iteration # 9 : max(cpu_time) wall_clock(cpu1) - | Time for this iteration : 2.594 s 5.109 s - | Charge density & force component update : 1.573 s 3.259 s - | Density mixing : 0.021 s 0.042 s - | Hartree multipole update : 0.008 s 0.017 s - | Hartree multipole summation : 0.299 s 0.592 s - | Hartree pot. SCF incomplete forces : 0.201 s 0.388 s - | Integration : 0.481 s 0.762 s - | Solution of K.-S. eqns. : 0.008 s 0.017 s - | Total energy evaluation : 0.001 s 0.000 s - - Partial memory accounting: - | Current value for overall tracked memory usage: - | Minimum: 0.677 MB (on task 0) - | Maximum: 0.677 MB (on task 0) - | Average: 0.677 MB - | Peak value for overall tracked memory usage: - | Minimum: 3.721 MB (on task 0 after allocating d_wave) - | Maximum: 3.745 MB (on task 1 after allocating d_wave) - | Average: 3.733 MB - | Largest tracked array allocation so far: - | Minimum: 1.628 MB (all_coords on task 0) - | Maximum: 1.628 MB (all_coords on task 1) - | Average: 1.628 MB - Note: These values currently only include a subset of arrays which are explicitly tracked. - The "true" memory usage will be greater. ------------------------------------------------------------- - ------------------------------------------------------------ - Computing monopole / dipole moments - | Total electronic charge [e] : 0.180000000000000E+03 - | Total ionic charge [e] : 0.180000000000000E+03 - | Total charge [e] : 0.000000000000000E+00 - | Total dipole moment [eAng] : 0.291817735833790E-04 0.415541790880761E-04 -0.459470451462385E-03 - | Absolute dipole moment : 0.462267694605846E-03 eAng / 0.222002032286256E-02 Debye . - Removing unitary transformations (pure translations, rotations) from forces on atoms. - Atomic forces before filtering: - | Net force on center of mass : -0.771996E-06 -0.675133E-06 -0.619516E-05 eV/A - | Net torque on center of mass: -0.169029E-04 -0.170005E-04 0.762398E-04 eV - Atomic forces after filtering: - | Net force on center of mass : 0.217781E-19 0.435562E-19 0.435562E-19 eV/A - | Net torque on center of mass: -0.161253E-18 -0.460979E-19 0.230489E-18 eV - - Energy and forces in a compact form: - | Total energy uncorrected : -0.118441303184828E+07 eV - | Total energy corrected : -0.118441303184828E+07 eV <-- do not rely on this value for anything but (periodic) metals - | Electronic free energy : -0.118441303184828E+07 eV - Total atomic forces (unitary forces cleaned) [eV/Ang]: - | 1 0.433397563009987E-03 0.585757860002957E-04 -0.320665674486224E-03 - | 2 -0.250155323972936E-03 -0.578627007355561E-03 0.319953271563094E-03 - | 3 -0.622114156351087E-04 0.445952552708970E-03 -0.320724801532806E-03 - | 4 -0.121030823401943E-03 0.740986686462952E-04 0.321437204455936E-03 - - ------------------------------------ - Start decomposition of the XC Energy - ------------------------------------ - X and C from original XC functional choice - Hartree-Fock Energy : 0.000000000 Ha 0.000000000 eV - X Energy : -792.437011372 Ha -21563.308208476 eV - C Energy : -10.610301600 Ha -288.720996495 eV - XC Energy w/o HF : -803.047312972 Ha -21852.029204971 eV - Total XC Energy : -803.047312972 Ha -21852.029204971 eV - ------------------------------------ - LDA X and C from self-consistent density - X Energy LDA : -759.827180764 Ha -20675.949569325 eV - C Energy LDA : -18.362981857 Ha -499.682159873 eV - ------------------------------------ - End decomposition of the XC Energy - ------------------------------------ - ------------------------------------------------------------- - Relaxation / MD: End force evaluation. : max(cpu_time) wall_clock(cpu1) - | Time for this force evaluation : 10.652 s 21.153 s - ------------------------------------------------------------- - Geometry optimization: Attempting to predict improved coordinates. - - Removing unitary transformations (pure translations, rotations) from forces on atoms. - Atomic forces before filtering: - | Net force on center of mass : 0.217781E-19 0.435562E-19 0.435562E-19 eV/A - | Net torque on center of mass: -0.161253E-18 -0.460979E-19 0.230489E-18 eV - Atomic forces after filtering: - | Net force on center of mass : 0.000000E+00 0.435562E-19 0.435562E-19 eV/A - | Net torque on center of mass: -0.161253E-18 -0.460979E-19 0.276587E-18 eV - Net remaining forces (excluding translations, rotations) in present geometry: - || Forces on atoms || = 0.578627E-03 eV/A. - Maximum force component is 0.578627E-03 eV/A. - Present geometry is converged. - ------------------------------------------------------------- - Final atomic structure: - x [A] y [A] z [A] - atom -0.00417371 -0.01235615 -0.00644413 O - atom 2.13657531 -0.13217828 -0.00288119 Pb - atom 2.00074444 2.00751378 0.00134131 O - atom -0.13996757 2.12738103 -0.00346233 Pb ------------------------------------------------------------- - ------------------------------------------------------------------------------- - Final output of selected total energy values: - - The following output summarizes some interesting total energy values - at the end of a run (AFTER all relaxation, molecular dynamics, etc.). - - | Total energy of the DFT / Hartree-Fock s.c.f. calculation : -1184413.031848285 eV - | Final zero-broadening corrected energy (caution - metals only) : -1184413.031848285 eV - | For reference only, the value of 1 Hartree used in FHI-aims is : 27.211384500 eV - - Before relying on these values, please be sure to understand exactly which - total energy value is referred to by a given number. Different objects may - all carry the same name 'total energy'. Definitions: - - Total energy of the DFT / Hartree-Fock s.c.f. calculation: - | Note that this energy does not include ANY quantities calculated after the - | s.c.f. cycle, in particular not ANY RPA, MP2, etc. many-body perturbation terms. - - Final zero-broadening corrected energy: - | For metallic systems only, a broadening of the occupation numbers at the Fermi - | level can be extrapolated back to zero broadening by an electron-gas inspired - | formula. For all systems that are not real metals, this value can be - | meaningless and should be avoided. - ------------------------------------------------------------------------------- - Methods described in the following list of references were used in this FHI-aims run. - If you publish the results, please make sure to cite these reference if they apply. - FHI-aims is an academic code, and for our developers (often, Ph.D. students - and postdocs), scientific credit in the community is essential. - Thank you for helping us! - - For any use of FHI-aims, please cite: - - Volker Blum, Ralf Gehrke, Felix Hanke, Paula Havu, Ville Havu, - Xinguo Ren, Karsten Reuter, and Matthias Scheffler - 'Ab initio molecular simulations with numeric atom-centered orbitals' - Computer Physics Communications 180, 2175-2196 (2009) - http://doi.org/10.1016/j.cpc.2009.06.022 - - - The ELSI infrastructure was used in your run to solve the Kohn-Sham electronic structure. - Please check out http://elsi-interchange.org to learn more. - If scalability is important for your project, please acknowledge ELSI by citing: - - V. W-z. Yu, F. Corsetti, A. Garcia, W. P. Huhn, M. Jacquelin, W. Jia, - B. Lange, L. Lin, J. Lu, W. Mi, A. Seifitokaldani, A. Vazquez-Mayagoitia, - C. Yang, H. Yang, and V. Blum - 'ELSI: A unified software interface for Kohn-Sham electronic structure solvers' - Computer Physics Communications 222, 267-285 (2018). - http://doi.org/10.1016/j.cpc.2017.09.007 - - - For the real-space grid partitioning and parallelization used in this calculation, please cite: - - Ville Havu, Volker Blum, Paula Havu, and Matthias Scheffler, - 'Efficient O(N) integration for all-electron electronic structure calculation' - 'using numerically tabulated basis functions' - Journal of Computational Physics 228, 8367-8379 (2009). - http://doi.org/10.1016/j.jcp.2009.08.008 - - Of course, there are many other important community references, e.g., those cited in the - above references. Our list is limited to references that describe implementations in the - FHI-aims code. The reason is purely practical (length of this list) - please credit others as well. - ------------------------------------------------------------- - Leaving FHI-aims. - Date : 20230426, Time : 002047.673 - - Computational steps: - | Number of self-consistency cycles : 58 - | Number of SCF (re)initializations : 5 - | Number of relaxation steps : 4 - - Detailed time accounting : max(cpu_time) wall_clock(cpu1) - | Total time : 68.455 s 133.429 s - | Preparation time : 0.852 s 1.209 s - | Boundary condition initalization : 0.001 s 0.002 s - | Grid partitioning : 0.721 s 1.205 s - | Preloading free-atom quantities on grid : 0.575 s 0.987 s - | Free-atom superposition energy : 0.299 s 0.560 s - | Total time for integrations : 24.510 s 47.720 s - | Total time for solution of K.-S. equations : 0.650 s 1.156 s - | Total time for EV reorthonormalization : 0.035 s 0.062 s - | Total time for density & force components : 32.514 s 61.670 s - | Total time for mixing : 0.819 s 1.562 s - | Total time for Hartree multipole update : 0.514 s 0.944 s - | Total time for Hartree multipole sum : 6.466 s 12.902 s - | Total time for total energy evaluation : 0.042 s 0.168 s - | Total time NSC force correction : 1.025 s 2.027 s - | Total time for scaled ZORA corrections : 0.000 s 0.000 s - - Partial memory accounting: - | Residual value for overall tracked memory usage across tasks: 0.000000 MB (should be 0.000000 MB) - | Peak values for overall tracked memory usage: - | Minimum: 3.721 MB (on task 0 after allocating d_wave) - | Maximum: 3.745 MB (on task 1 after allocating d_wave) - | Average: 3.733 MB - | Largest tracked array allocation: - | Minimum: 1.628 MB (all_coords on task 0) - | Maximum: 1.628 MB (all_coords on task 1) - | Average: 1.628 MB - Note: These values currently only include a subset of arrays which are explicitly tracked. - The "true" memory usage will be greater. - - Have a nice day. diff --git a/AppOutputExtractor/FHIaims/old_MLTrainingDataGenerator.py b/AppOutputExtractor/FHIaims/old_MLTrainingDataGenerator.py new file mode 100644 index 0000000..1fdfa9f --- /dev/null +++ b/AppOutputExtractor/FHIaims/old_MLTrainingDataGenerator.py @@ -0,0 +1,209 @@ +''' +Author: Dong-Gi Kang +Prepare ML-IP data using FHI-aims output +Training data type: vibrational mode of a cluster +''' +import os +import sys +import numpy as np +from AppOutputExtractor.FHIaims.FHIaimsOutputExtractor import extractor + +class ML_train_generator(extractor): + + def __init__(self, app_version='22', tag=None): + app_output = './aims.out' + + self.extractor = extractor() + self.extractor.set_output_filepath(app_output) + + self.species = self.extractor.get_species + self.no_atoms = self.extractor.get_no_atoms + self.geometries = self.extractor.get_geometries + self.order = self.extractor.get_atom_order + self.forces = self.extractor.get_forces + + try: + self.vib_eigvecs = self.extractor.get_vib_eigvec + except: + pass + + self.ucl_id = 'uccatka' + self.job_time = '2:00:00' + self.job_name = 'test' + self.memory = '1' + self.cpu_core = '40' # for Young 40 core = 1 node + self.payment = 'Gold' + self.budgets = 'UCL_chemM_Woodley' + self.path_binary = '/home/uccatka/software/fhi-aims.221103/build/aims.221103.scalapack.mpi.x' + self.path_fhiaims_species = '/home/uccatka/software/fhi-aims.221103/species_defaults/defaults_2020/light' + self.step_size = 0.05 + return None + + + @property + def mod_xyz_w_vib(self): + ''' Modify LM geometry to array of vibrational mode frames ''' + Lambda = len(np.arange(-1, 1+self.step_size, self.step_size)) * self.no_atoms*3 + self.mod_sp = np.zeros((Lambda, self.no_atoms, 3)) + cnt = 0 + for i in range(self.no_atoms * 3): + for numj, j in enumerate(np.arange(-1, 1+self.step_size, self.step_size)): + j = np.round(j, 2) + frame = self.geometries[-1] + self.vib_eigvecs[i] * j + self.mod_sp[cnt] = np.round(frame, 8) + cnt += 1 + self.mod_sp = np.reshape(self.mod_sp, (self.no_atoms*3, len(np.arange(-1, 1+self.step_size, self.step_size)), self.no_atoms, 3)) + return self.mod_sp + + + @property + def geometry_for_sp(self): + ''' Convert the modified geometry (mod_xyz_w_vib) to {geometry.in} format for FHI-aims ''' + placer = np.full((self.no_atoms, 1), 'atom') + placer_species = np.reshape(self.order, (-1, 1)) + shape = np.shape(self.mod_sp) + self.for_sp = np.empty((shape[0], shape[1], self.no_atoms, 5), dtype=object) + + for i in range(shape[0]): + for j in range(shape[1]): + form = np.concatenate((placer, self.mod_sp[i][j], placer_species), axis=1) + self.for_sp[i][j] = form + return self.for_sp + + + @property + def xyz_from_opti(self): + ''' prepare training data from every SCF converged cycles of a optimisation ''' + train_xyz = 'xyz_from_opti.xyz' + exist = [x for x in os.listdir('./') if train_xyz in x] + if len(exist) != 0: + os.remove(exist[0]) + else: pass + for i in range(len(self.extractor.set_scf_blocks)): + self.energy = self.extractor.get_total_energy(i) + self.geometry = self.geometries[i] + self.force = self.forces[i] + xyz = np.round(np.concatenate((self.geometry, self.force), axis=1), 9) + xyz = np.concatenate((self.order, xyz), axis=1) + + with open('xyz_from_opti.xyz', 'a') as f: + f.write(f'{self.no_atoms}\n') + f.write(f'Properties-species:S:1:pos:R:3:forces:R:3 energy={self.energy} pbc="F F F"\n') + np.savetxt(f, xyz, fmt="%s", delimiter=" ") + print(f"total of {i+1} SCF converged structures are prepared in {train_xyz}") + + + def make_sp_control(self, path): + ''' Write {control.in} file ''' + basis_set_files = [os.path.join(self.path_fhiaims_species, x) for x in os.listdir(self.path_fhiaims_species)] + basis_set_all = [x.split('_')[1] for x in os.listdir(self.path_fhiaims_species)] + basis_set_index = [basis_set_all.index(x) for x in basis_set_all if x in self.species] + + path = os.path.join(path, 'control.in') + with open(path, 'a') as f: + f.write("#\n") + f.write("xc pbesol\n") + f.write("spin none\n") + f.write("relativistic atomic_zora scalar\n") + f.write("charge 0.\n\n") + f.write("# SCF convergence\n") + f.write("occupation_type gaussian 0.01\n") + f.write("mixer pulay\n") + f.write("n_max_pulay 10\n") + f.write("charge_mix_param 0.5\n") + f.write("sc_accuracy_rho 1E-5\n") + f.write("sc_accuracy_eev 1E-3\n") + f.write("sc_accuracy_etot 1E-6\n") + f.write("sc_accuracy_forces 1E-4\n") + f.write("sc_iter_limit 1500\n") + f.write("# Relaxation\n\n") + #f.write("relax_geometry bfgs 1.e-3\n") + for i in basis_set_index: + with open(basis_set_files[i], 'r') as ff: + lines = ff.read() + f.write(lines) + f.write('\n') + return None + + + def make_job_submit(self, path): + ''' Write 'submit.sh' job script for SGE system ''' + _, last_part = os.path.split(path) + _, second_last_part = os.path.split(os.path.dirname(path)) + + # Combine the last two parts + last_two_parts = f"{second_last_part}_{last_part}" + + path = os.path.join(path, 'submit.sh') + with open(path, 'a') as f: + f.write("#!/bin/bash -l\n") + f.write('\n') + f.write("#$ -S /bin/bash\n") + f.write(f"#$ -l h_rt={self.job_time}\n") + f.write(f"#$ -l mem={self.memory}G\n") + f.write(f"#$ -N p{last_two_parts}\n") + f.write(f"#$ -pe mpi {self.cpu_core}\n") + f.write("#$ -cwd\n") + f.write("\n") + f.write(f"#$ -P {self.payment}\n") + f.write(f"#$ -A {self.budgets}\n") + + f.write("module load gerun\n") + f.write("module load userscripts\n") + f.write("module unload -f compilers mpi gcc-libs\n") + f.write("module load gcc-libs/4.9.2\n") + f.write("module unload -f compilers mpi\n") + f.write("module load beta-modules\n") + f.write("module load openblas/0.3.7-serial/gnu-4.9.2\n") + f.write("module load compilers/intel/2019/update5\n") + f.write("module load mpi/intel/2018/update3/intel\n") + + f.write("\n") + f.write("#$ -m e\n") + f.write(f"#$ -M {self.ucl_id}@ucl.ac.uk\n") + f.write("\n") + f.write(f"gerun {self.path_binary} > aims.out\n") + + + + + +import argparse +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument("--eigenvector", type=str, help="A string of space-separated eigenvector indices. For example, '7 8 9 10'") + args = parser.parse_args() + + # Split the string of indices into a list of integers + indices = list(map(int, args.eigenvector.split())) + + app_output = './aims.out' + step_size = 0.05 + ml = ML_train_generator() + + ml.mod_xyz_w_vib + sp_frame = ml.geometry_for_sp + shape = np.shape(sp_frame) + if not os.path.exists('sp'): + os.mkdir('sp') + else: pass + + for i in indices: # Now we only iterate over the specified indices + if not os.path.exists(os.path.join('sp', str(i+1))): + os.mkdir(f'sp/{str(i)}') + else: pass + + for numj, j in enumerate(np.arange(-1, 1+step_size, step_size)): + j = str(np.round(j, 2)) + os.mkdir(f'sp/{str(i)}/lambda_{j}') + with open(f'sp/{i}/lambda_{j}/geometry.in', 'w') as f: + for row in sp_frame[i-1][numj]: + line = ' '.join(str(x) for x in row) + f.write(line + '\n') + ml.make_sp_control(f'sp/{i}/lambda_{j}') + ml.make_job_submit(f'sp/{i}/lambda_{j}') + os.chdir(f'sp/{i}/lambda_{j}') + os.system('qsub submit.sh') + os.chdir('../../../') + + diff --git a/AppOutputExtractor/FHIaims/retrieve_sp_extxyz.py b/AppOutputExtractor/FHIaims/retrieve_sp_extxyz.py new file mode 100644 index 0000000..60625a0 --- /dev/null +++ b/AppOutputExtractor/FHIaims/retrieve_sp_extxyz.py @@ -0,0 +1,57 @@ +from AppOutputExtractor.FHIaims.FHIaimsOutputExtractor import extractor +import os +from itertools import groupby +import numpy as np + +eigvec_path = [os.path.join('sp', x) for x in os.listdir('sp')] +sp_path = [os.path.join(dirpath, fname) for dirpath in eigvec_path for fname in os.listdir(dirpath)] +lambda_path = [os.path.join(dirpath, fname) for dirpath in sp_path for fname in os.listdir(dirpath) if fname == 'aims.out'] +aims_out_path = sorted(lambda_path, key=lambda x: (int(x.split('/')[1]), float(x.split('/')[2].split('_')[1]))) +aims_out_path = [list(group) for key, group in groupby(aims_out_path, lambda x: int(x.split('/')[1]))] +ex = extractor() +if not os.path.exists('ext_xyz'): + os.mkdir('ext_xyz') + +for numi, i in enumerate(aims_out_path): + total_energy = [] + geometry = [] + + for j in i: + print(j) + ex.set_output_filepath(j) + ex.set_scf_blocks + + ex.get_no_atoms + + ex.get_sp_geometries(j) + ex.get_sp_atom_order() + ex.get_sp_species() + ex.get_forces + force_shape = np.shape(ex.get_forces) + get_forces = np.reshape(ex.get_forces, (force_shape[1], force_shape[2])) + + form = np.concatenate((ex.get_sp_atom_order(), ex.get_sp_geometries(j), get_forces), axis=1) + + total_energy.append(ex.get_total_energy()) + geometry.append(form) + + for numk, k in enumerate(total_energy): + with open(f"ext_xyz/ext_{j.split('/')[1]}_eigv.xyz", 'a') as f: + f.write(str(force_shape[1]) + '\n') + f.write(f'Properties-species:S:1:pos:R:3:forces:R:3 energy={total_energy[numk]} pbc="F F F"\n') + np.savetxt(f, geometry[numk], fmt="%s", delimiter=" ") + + + #ex.get_vib_eigvec +# for k in range(len(ex.set_scf_blocks)): +# #print(ex.get_total_energy(k)) +# #print('Geometry') +# print(ex.get_sp_geometries) +# #print('Atomic forces') +# #print(ex.get_forces(k)) +# #print('Eigenvector of vibrational modes') +# #print(ex.get_vib_eigvec) +# #print() +# #print() +# ex.get_total_energy + diff --git a/AppOutputExtractor/FHIaims/test.py b/AppOutputExtractor/FHIaims/test.py deleted file mode 100644 index 73a58d3..0000000 --- a/AppOutputExtractor/FHIaims/test.py +++ /dev/null @@ -1,45 +0,0 @@ - -import re - -def find_pattern_in_file(file_path, pattern): - with open(file_path) as f: - # Read the file line by line - lines = f.readlines() - - # Use regular expression to search for the pattern - regex = re.compile(pattern) - matching_lines = [i+1 for i, line in enumerate(lines) if regex.search(line)] - - # Return the line numbers where the pattern is found - return matching_lines - -def find_pattern_with_last_word(filename, pattern): - - with open(filename, 'r') as file: - lines = file.readlines() - - matches = [] - - for i, line in enumerate(lines): - match = re.search(pattern, line) - if match: - last_word = line.strip().split()[-1] - matches.append((i+1, last_word)) - - return matches - - -if __name__ == '__main__': - - file_path = "/Users/woongkyujee/Desktop/Python/FHI22_samples/runs/run_1/FHIaims.out" - pattern = r"SCF" # The pattern you want to search for - pattern = 'Begin self-consistency iteration #' - - #matching_lines = find_pattern_in_file(file_path, pattern) - #print("Matching line numbers:", matching_lines) - - - print('----- test 2 -----') - - matches = find_pattern_with_last_word(file_path,pattern) - print(matches) diff --git a/AppOutputExtractor/FHIaims/testing_extractor.py b/AppOutputExtractor/FHIaims/testing_extractor.py new file mode 100644 index 0000000..7ae55cc --- /dev/null +++ b/AppOutputExtractor/FHIaims/testing_extractor.py @@ -0,0 +1,14 @@ +from AppOutputExtractor.FHIaims.FHIaimsOutputExtractor import extractor + + +app_output = './aims.out' + +extractor = extractor() +extractor.set_output_filepath(app_output) + +species = extractor.get_species +no_atoms = extractor.get_no_atoms +geometries = extractor.get_geometries +order = extractor.get_atom_order +forces = extractor.get_forces +print(forces) diff --git a/AppOutputExtractor/__pycache__/OutputExtractor.cpython-38.pyc b/AppOutputExtractor/__pycache__/OutputExtractor.cpython-38.pyc new file mode 100644 index 0000000..0428200 Binary files /dev/null and b/AppOutputExtractor/__pycache__/OutputExtractor.cpython-38.pyc differ diff --git a/AppOutputExtractor/__pycache__/OutputExtractor.cpython-39.pyc b/AppOutputExtractor/__pycache__/OutputExtractor.cpython-39.pyc new file mode 100644 index 0000000..c990e1a Binary files /dev/null and b/AppOutputExtractor/__pycache__/OutputExtractor.cpython-39.pyc differ diff --git a/AppOutputExtractor/__pycache__/__init__.cpython-38.pyc b/AppOutputExtractor/__pycache__/__init__.cpython-38.pyc new file mode 100644 index 0000000..fde4232 Binary files /dev/null and b/AppOutputExtractor/__pycache__/__init__.cpython-38.pyc differ diff --git a/AppOutputExtractor/__pycache__/__init__.cpython-39.pyc b/AppOutputExtractor/__pycache__/__init__.cpython-39.pyc new file mode 100644 index 0000000..3bc35b0 Binary files /dev/null and b/AppOutputExtractor/__pycache__/__init__.cpython-39.pyc differ diff --git a/FHIaimsMolecule.py b/FHIaimsMolecule.py new file mode 100644 index 0000000..8f97f1e --- /dev/null +++ b/FHIaimsMolecule.py @@ -0,0 +1,101 @@ +from NonPeriodic.Molecule import BaseMolecule +from NonPeriodic.Atom import BaseAtom + +import os,re +import numpy as np + +''' + +''' + +class atom(BaseAtom): + + def __init__(self,atom_type='atom',atom_attr='default',x=0.,y=0.,z=0.): + + ''' + + ''' + + super().__init__(atom_type=atom_type,atom_attr=atom_attr,x=x,y=y,z=z) + + +class molecule(BaseMolecule): + + def __init__(self,geometry_file): + + ''' + + ''' + + super().__init__() # BaseMolecule NonPeriodic/ (number_of_atoms,atom_list) + + try: + with open(geometry_file,'r') as f: + self.file_exist = True + self.file_path = geometry_file + pattern = re.compile(r"atom") + + for line in f: + if pattern.search(line): + ls = line.split() + new_atom = atom(ls[4],ls[0],ls[1],ls[2],ls[3]) + self.add_atom(new_atom) + + except FileNotFoundError as e: + self.file_exist = False + print(e) # pass + + def is_exist(self): + return self.file_exist + + +''' + +''' + +def calculate_rmsd_molecules(moleculeA,moleculeB): + + if moleculeA.get_number_of_atoms() != moleculeB.get_number_of_atoms(): + return False + else: + rmsd = 0. + + for atomA,atomB in zip(moleculeA.get_atomlist(),moleculeB.get_atomlist()): + cartA = np.array(atomA.get_cart()) + cartB = np.array(atomB.get_cart()) + dev = np.linalg.norm(cartA - cartB) + rmsd = rmsd + dev + + try: + rmsd = rmsd/(float(moleculeA.get_number_of_atoms())*3.) + return rmsd + except ZeroDivisionError as e: + print(e) + return None + + +''' + +''' + +if __name__ == '__main__': + + print('Molecule - 1') + fmol = molecule('/Users/woongkyujee/Desktop/Python/AppOutputAnalysis/unit_tests/run_1/geometry.in') + print(fmol.is_exist()) + fmol.show_info() + + print('Molecule - 2') + fmol_2 = molecule('/Users/woongkyujee/Desktop/Python/FHI22_samples/runs/run_2/geometry.in') + fmol_2.show_info() + #fmol_2.get_atom(2).show_info() + + print('rmsd test 1') + print(calculate_rmsd_molecules(fmol,fmol_2)) + + print('rmsd test 2') + fmolA = molecule('geoA.txt') + fmolB = molecule('geoB.txt') + print(fmolA.show_info()) + print(fmolB.show_info()) + print(calculate_rmsd_molecules(fmolA,fmolB)) diff --git a/FHIaimsOutputExtractor.py b/FHIaimsOutputExtractor.py new file mode 100644 index 0000000..6fd5816 --- /dev/null +++ b/FHIaimsOutputExtractor.py @@ -0,0 +1,510 @@ +''' +Author: Dr Woongkyu Jee, Dong-Gi Kang +''' + + +# +import time +from AppOutputExtractor.OutputExtractor import BaseExtractor +from AppOutputExtractor.FHIaims.FHIaimsMolecule import molecule as fmol +from AppOutputExtractor.FHIaims.FHIaimsMolecule import calculate_rmsd_molecules + +from ShellCommand import shellcommand +import ParsingSupport + +import os,re +import numpy as np +import string,json + +class extractor(BaseExtractor): + + def __init__(self,app_version='22',tag=None): + ''' + ''' + super().__init__(app='FHIaims',version=app_version) + + # set app output patterns + module_path = os.path.dirname(os.path.abspath(__file__)) + '/OutputPattern' # getting this module path, '__file__' + self.patterns = self.load_patterns(module_path) + + # memo + self.tag = tag + + # shellcommand obj + self.shell = shellcommand() + self.scf_converged_blocks = [] + + def set_output_filepath(self,path): + if os.path.exists(path): + self.output_filepath = path + else: + self.output_filepath = None + print('in {} method "set_output_filepath()", cannot find the file at: "{}" '.format(__file__,path)) + + def set_input_geometry_filepath(self,path): + if os.path.exists(path): + self.input_geometry_filepath = path + self.input_geometry = fmol(path) + else: + self.input_geometry_filepath = None + print('in {} method "set_input_geometry_filepath()", cannot find the file at: "{}" '.format(__file__,path)) + + def set_output_geometry_filepath(self,path): + if os.path.exists(path): + self.output_geometry_filepath = path + self.output_geometry = fmol(path) + else: + self.output_geometry_filepath = None + print('in {} method "set_output_geometry_filepath()", cannot find the file at: "{}" '.format(__file__,path)) + + def check_filepaths(self): + #print('App output : {}'.format(self.output_filepath)) + #print('geometry input : {}'.format(self.input_geometry_filepath)) + #print('geometry output: {}'.format(self.output_geometry_filepath)) + ''' + field: (0) AppOutput (1) InputGeometry (2) OutputGeometry + ''' + return [self.output_filepath,self.input_geometry_filepath,self.output_geometry_filepath] + + def get_input_molecule(self): + checker = self.check_filepaths()[1] + if checker: + return self.input_geometry + else: + print('input geometry is not loaded!') + + def get_output_molecule(self): + checker = self.check_filepaths()[2] + if checker: + return self.output_geometry + else: + print('output geometry is not loaded!') + + + + ''' + Interaction with app output file + ''' + + def check_calculation_success(self): + self.shell.set_tarfile(self.output_filepath) + cmd = self.shell.grep(self.patterns['SUCCESS']['pattern']) + shell_res = self.shell.execute(cmd) + + if shell_res != None: + self.output_success_tag = True + else: + self.output_success_tag = False + + return self.output_success_tag #!!! + + def check_calculation_runtime(self): + # wall clock time + self.shell.set_tarfile(self.output_filepath) + cmd = self.shell.pipe(\ + self.shell.grep(self.patterns['APP_RUNTIME']['pattern'])\ + ,self.shell.awk(self.patterns['APP_RUNTIME']['wtime_token']) + ) + target = self.shell.execute(cmd) #!!! + + try: + target = float(target) + return target + except: + print('failed to get calculation wtime') + return None + + def check_parallel_task(self): + # used cpus + self.shell.set_tarfile(self.output_filepath) + cmd = self.shell.pipe(\ + self.shell.grep(self.patterns['APP_RESOURCE_USED']['pattern'])\ + ,self.shell.awk(self.patterns['APP_RESOURCE_USED']['token']) + ) + target = self.shell.execute(cmd) #!!! + + try: + target = int(target) + return target + except: + print('failed to get parallel task number, recheck the app output file') + return None + + + + ''' + Loading SCF converged blocks ... possibly useful for further app output collation + ''' + + #@property + def set_scf_blocks(self) -> list: + ''' + * special blocks: + self.scf_converged_blocks[0] -> first SCF converged blocks [line_start,line_end] + self.scf_converged_blocks[-1]-> final SCF converged blocks + ''' + pattern = self.patterns['BEGIN_SCF']['pattern'].replace("'","") + self.total_lnumber, self.scf_block_lines = \ + ParsingSupport.find_pattern_with_last_word(self.output_filepath,pattern) \ + # GET LINE NUMBERS OF SCF (CONVERGED) BLOCKS + + self.scf_converged_blocklines = [] #!!! + self.scf_converged_blocks = [] #!!! + + # IF ITEM IN ITERABLE SAVE THE LINE NUMBEERS [START,END] + for i, item in enumerate(self.scf_block_lines[:-1]): + curr_tag = int(self.scf_block_lines[i][1]) + next_tag = int(self.scf_block_lines[i+1][1]) + + if next_tag < curr_tag: + + block_start = self.scf_block_lines[i][0] + block_end = self.scf_block_lines[i+1][0] + + self.scf_converged_blocklines.append([block_start,block_end]) + + # FIANL SCF CONVERGED BLOCK (BEFORE APP FINALISATION) + block_start = self.scf_block_lines[-1][0] + block_end = self.total_lnumber + self.scf_converged_blocklines.append([block_start,block_end]) + + # SAVE THE BLOCKS ... 'self.scf_converged_blocks' -> python list + for item in self.scf_converged_blocklines: + self.scf_converged_blocks.append(ParsingSupport.get_lines(self.output_filepath,item[0],item[1])) + return self.scf_converged_blocks + + + #@property + def get_species(self, atom_order) -> list: + #get_species = [x for x in self.get_atom_order.tolist()] + get_species = list(set([item for sublist in atom_order for item in sublist])) + get_species = sorted(get_species) + return get_species + + + #@property + def get_no_atoms(self) -> int: + with open(self.output_filepath, 'r') as f: + lines = f.readlines() + for i in lines: + if self.patterns['NO_ATOMS']['pattern'] in i: + no_atoms = int(i.split()[5]) + return no_atoms + + + # REVIEW-delete: property decorator on scf_converged_blocks does same function + #def get_scf_blocks(self): + # return self.scf_converged_blocks + + + @property + def get_number_of_scf_blocks(self) -> int: + return len(self.set_scf_blocks) + + ''' + AppOutput Collation Methods + ''' + + def get_total_energy(self, block=-1): + pattern_str = self.patterns['SCF_ENERGY']['pattern'].replace("'","") + token = int(self.patterns['SCF_ENERGY']['token']) - 1 + pattern = re.compile(pattern_str) + for i in self.set_scf_blocks[block]: + matching = pattern.search(i) + if matching: + target = float(i.split()[ token ]) + break + return target + + + + #@property + def get_atom_order(self, no_atoms, block=-1) -> np.ndarray: + self.set_scf_blocks() + + pattern_str = self.patterns['SCF_GEOMETRY_END']['pattern'].replace("'", "") + pattern = re.compile(pattern_str) + + start_index = None + self.match_atom = np.empty((no_atoms), dtype=object) + for numj, j in enumerate(self.scf_converged_blocks[block]): + matching = pattern.search(j) + if matching: + start_index = numj + 2 + elif start_index is not None and j.strip() == '': + end_index = numj - 1 + atomic_structure = self.scf_converged_blocks[block][start_index: end_index] + for numk, k in enumerate(atomic_structure): + numbers = [x for x in k.split()] + self.match_atom[numk] = numbers[-1] + break + self.match_atom = np.reshape(self.match_atom, (no_atoms, 1)) + return self.match_atom + + #@property + def get_geometries(self, no_atoms, block=-1) -> np.ndarray: + #if not self.scf_converged_blocks: + self.set_scf_blocks() + + pattern_str = self.patterns['SCF_GEOMETRY_BEGIN']['pattern'].replace("'", "") + pattern = re.compile(pattern_str) + + if block == -1 or block == len(self.scf_converged_blocks)-1: + pattern_str = self.patterns['SCF_GEOMETRY_END']['pattern'].replace("'", "") + pattern = re.compile(pattern_str) + + start_index = None + self.geo = np.zeros((no_atoms, 3)) + cnt = 0 + for numj, j in enumerate(self.scf_converged_blocks[block]): + matching = pattern.search(j) + if matching: + start_index = numj + 2 + end_index = numj + no_atoms + 2 + atomic_structure = self.scf_converged_blocks[block][start_index: end_index] + for numk, k in enumerate(atomic_structure): + numbers = [x for x in k.split()] + self.geo[numk] = list(map(float, numbers[1:4])) # Convert the rest to float and store in self.geo + cnt += 1 + start_index = None + + if cnt == 0: + cnt = 1 + else: pass + self.geo = np.reshape(self.geo, (cnt, int(no_atoms), 3)) + return self.geo + + def get_sp_geometries(self, path) -> np.ndarray: + new_path = os.path.join(os.path.dirname(path), 'geometry.in') + with open(new_path, 'r') as f: + lines = f.readlines() + lines = [x.split() for x in lines] + lines = np.array(lines) + shape = np.shape(lines) + #atom_str = np.reshape(lines[:,0], (shape[0], -1)) + self.atom_label = np.reshape(lines[:,-1], (shape[0], -1)) + self.coordinate = lines[:, 1:-1].astype(float) + return self.coordinate + + def get_sp_atom_order(self): + return self.atom_label + + def get_sp_species(self): + return list(set(self.atom_label.flatten().tolist())) + + #@property + def get_forces(self, no_atoms, block=-1) -> np.ndarray: + pattern_str = self.patterns['SCF_FORCE']['pattern'].replace("'", "") + pattern = re.compile(pattern_str) + start_index = None + self.forces = np.zeros((no_atoms, 3)) + cnt = 0 + for numj, j in enumerate(self.scf_converged_blocks[block]): + matching = pattern.search(j) + if matching: + start_index = numj + 1 + elif start_index is not None and j.strip() == '': + end_index = numj + force = self.scf_converged_blocks[block][start_index:end_index] + for numk, k in enumerate(force): + numbers = list(map(float, k.strip().split()[-3:])) + self.forces[numk] = numbers + start_index = None + cnt += 1 + + self.forces = np.reshape(self.forces, (cnt, 12, 3)) + return self.forces + + + #@property + def get_vib_eigvec(self, no_atoms) -> np.ndarray: + ''' + Read whole file contents (not necessary to read in blocks as we need all eigenvector of vibrational mode + ''' + check_vib = [x for x in os.listdir('./') if 'vibration' in x] + if len(check_vib) == 0: + raise FileNotFoundError("Cannot find 'vibration' directory") + else: + check_vib = [x for x in os.listdir('./') if 'vibration' in x][0] + + vib_xyz = [os.path.join(check_vib, x) for x in os.listdir(check_vib) if '_' and '.xyz' in x][0] + with open(vib_xyz, 'r') as f: + lines = f.readlines() + + self.eigvec = np.zeros((no_atoms*3, no_atoms, 3)) + start_index = None + block_counter = -1 + for numi, i in enumerate(lines): + if 'frequency' in i: + start_index = numi + 1 + block_counter += 1 + elif start_index is not None and (i.strip().split()[0] in ['Al', 'F']): + data = list(map(float, i.strip().split()[-3:])) # convert last three elements to float + atom_index = numi - start_index + self.eigvec[block_counter, atom_index, :] = data + elif i.strip() == str(self.get_no_atoms): + start_index = None + + return self.eigvec + + + def get_dipole(self,block=-1): + pattern_str = self.patterns['DIPOLE']['pattern'].replace("'","") + token = int(self.patterns['DIPOLE']['token']) - 1 + pattern = re.compile(pattern_str) + + for line in self.scf_converged_blocks[block]: + matching = pattern.search(line) + if matching: + target = float(line.split()[ token ]) + break + return target + + def get_dipole_moment(self,block=-1): + pattern_str = self.patterns['DIPOLE_MOMENT']['pattern'].replace("'","") + token_x = int(self.patterns['DIPOLE_MOMENT']['token_x']) - 1 + token_y = int(self.patterns['DIPOLE_MOMENT']['token_y']) - 1 + token_z = int(self.patterns['DIPOLE_MOMENT']['token_z']) - 1 + pattern = re.compile(pattern_str) + + target = [] + + for line in self.scf_converged_blocks[block]: + matching = pattern.search(line) + if matching: + target.append( float(line.split()[ token_x ]) ) + target.append( float(line.split()[ token_y ]) ) + target.append( float(line.split()[ token_z ]) ) + break + return target + + def get_homolumo(self,block=-1): + + ''' + field: (0) HOMO (1) LUMO (2) HOMO-LUMO + ''' + target = [] + + # HOMO + pattern_str = self.patterns['HOMO']['pattern'].replace("'","") + token = int(self.patterns['HOMO']['token']) - 1 + pattern = re.compile(pattern_str) + for line in self.scf_converged_blocks[block]: + matching = pattern.search(line) + if matching: + target.append( float(line.split()[ token ]) ) + break + # LUMO + pattern_str = self.patterns['LUMO']['pattern'].replace("'","") + token = int(self.patterns['LUMO']['token']) - 1 + pattern = re.compile(pattern_str) + for line in self.scf_converged_blocks[block]: + matching = pattern.search(line) + if matching: + target.append( float(line.split()[ token ]) ) + break + # HOMOLUMO GAP + pattern_str = self.patterns['HOMOLUMO']['pattern'].replace("'","") + token = int(self.patterns['HOMOLUMO']['token']) - 1 + pattern = re.compile(pattern_str) + for line in self.scf_converged_blocks[block]: + matching = pattern.search(line) + if matching: + target.append( float(line.split()[ token ]) ) + break + return target + + + # Getters Miscs + + def get_patterns(self): + # return type 'json' + return self.patterns + + def get_tag(self): + return self.tag + + + + +if __name__ == '__main__': + + file_root = '/Users/woongkyujee/Desktop/Python/FHI22_samples/runs/run_1' + main_output = file_root + '/FHIaims.out' + input_geo = file_root + '/geometry.in' + output_geo = file_root + '/geometry.in.next_step' + + ext2 = extractor() + ext2.set_output_filepath(main_output) + ext2.set_input_geometry_filepath(input_geo) + ext2.set_output_geometry_filepath(output_geo) + + print('check filepaths()') + print(ext2.check_filepaths()) # if None in ext2.check_filepaths(): + print('calculation success check: {}'.format(ext2.check_calculation_success())) + + print('calculation runtime') + rtime = ext2.check_calculation_runtime() + print(rtime) + + print('calculation parallel tasks') + ptask = ext2.check_parallel_task() + print(ptask) + + + + + + + + ### EXTRACTION + + ext2.set_scf_blocks() # load scf blocks + + # Energy Check + print('init E') + init_E = ext2.get_total_energy(0) + print(init_E) + print('final E') + final_E = ext2.get_total_energy() + print(final_E) + + # Dipole Check + print('init P') + init_p = ext2.get_dipole(0) + print(init_p) + initial_p_elem = ext2.get_dipole_moment(0) + print(initial_p_elem) + + print('final P') + final_p = ext2.get_dipole() + print(final_p) + final_p_elem = ext2.get_dipole_moment(0) + print(final_p_elem) + + # HOMOLUMO CHECK + print('init homo-lumo, list [homo,lumo,homo-lumo]') + init_hl = ext2.get_homolumo(0) + print(init_hl) + print('final homo-lumo, list [homo,lumo,homo-lumo]') + final_hl = ext2.get_homolumo() + print(final_hl) + + ''' + Unit test with 'ext2' instance + ''' + print('--- input') + ext2.input_geometry.show_info() + print('--- output') + ext2.output_geometry.show_info() + print('in - out geometry rmsd') + rmsd = calculate_rmsd_molecules(ext2.input_geometry,ext2.output_geometry) + print(rmsd) + + #print('output check --') + #ext2.check_output_success() + #print(ext2.output_success_tag) + + ''' + Unit test getting SCF Blocks + ''' diff --git a/FHIaimsVib.py b/FHIaimsVib.py new file mode 100644 index 0000000..6adbd8a --- /dev/null +++ b/FHIaimsVib.py @@ -0,0 +1,133 @@ +from AppOutputExtractor.OutputExtractor import BaseExtractor +from AppOutputExtractor.FHIaims.FHIaimsMolecule import molecule as fmol +from AppOutputExtractor.FHIaims.FHIaimsMolecule import calculate_rmsd_molecules +from AppOutputExtractor.FHIaims.FHIaimsOutputExtractor import extractor + +from ShellCommand import shellcommand +import ParsingSupport + +import os +import shutil +import string,json + + +class aimsvibcalc(BaseExtractor): + + def __init__(self, app_version='22', tag=None): + ''' + ''' + self.extractor = extractor() + app_output = './aims.out' + self.extractor.set_output_filepath(app_output) + self.species = self.extractor.get_species + #print(self.species) + + self.ucl_id = 'uccatka' + self.job_time = '2:00:00' + #self.job_name = 'test' + self.memory = '2' + self.cpu_core = '40' # for Young 40 core = 1 node + self.payment = 'Gold' + self.budgets = 'UCL_chemM_Woodley' + self.path_binary = '/home/uccatka/software/fhi-aims.221103/build/aims.221103.scalapack.mpi.x' + self.vib_path_binary = '/home/uccatka/software/fhi-aims.221103/build/src/vibrations/numerical_vibrations.pl' + self.path_fhiaims_species = '/home/uccatka/software/fhi-aims.221103/species_defaults/defaults_2020/light' + self.step_size = 0.05 + + super().__init__(app='FHIaims',version=app_version) + + # set app output patterns + module_path = os.path.dirname(os.path.abspath(__file__)) + '/OutputPattern' # getting this module path, '__file__' + self.patterns = self.load_patterns(module_path) + + # memo + self.tag = tag + + # shellcommand obj + self.shell = shellcommand() + + return None + + def make_job_submit(self, job_name, loc='./vibration', step_size='0.0025'): + ''' Write 'submit.sh' job script for SGE system ''' + path = os.path.join(loc, 'submit.sh') + with open(path, 'a') as f: + f.write("#!/bin/bash -l\n") + f.write("\n") + f.write("#$ -S /bin/bash\n") + f.write(f"#$ -l h_rt={self.job_time}\n") + f.write(f"#$ -l mem={self.memory}G\n") + f.write(f"#$ -N {job_name}\n") + f.write(f"#$ -pe mpi {self.cpu_core}\n") + f.write("#$ -cwd\n") + f.write("\n") + f.write(f"#$ -P {self.payment}\n") + f.write(f"#$ -A {self.budgets}\n") + f.write("\n") + f.write("#$ -m e\n") + f.write(f"#$ -M {self.ucl_id}@ucl.ac.uk\n") + f.write("\n") + f.write("module purge\n") + f.write("module load gerun\n") + f.write("module load userscripts\n") + f.write("module load gcc-libs/4.9.2\n") + f.write("module unload -f compilers mpi\n") + f.write("module load beta-modules\n") + f.write("module load gcc-libs/10.2.0\n") + f.write("module load openblas/0.3.7-serial/gnu-4.9.2\n") + f.write("module load compilers/intel/2019/update5\n") + f.write("module load mpi/intel/2018/update3/intel\n\n") + + f.write(f"gerun {self.vib_path_binary} {job_name}_{step_size} {step_size} > vibres.out\n") + + @property + def vib_calc_prep(self): + #shutil.copy('./control.in', './vibration') + geo_next = 'geometry.in.next_step' + geo = 'geometry.in' + vib_dir = 'vibration' + geometry_files = [x for x in os.listdir('./') if '.in' in x] + if geo_next in geometry_files: + shutil.copy(geo_next, f'{vib_dir}/{geo}') + shutil.copy('hessian.aims', vib_dir) + else: + shutil.copy(geo, vib_dir) + + basis_set_files = [os.path.join(self.path_fhiaims_species, x) for x in os.listdir(self.path_fhiaims_species)] + basis_set_all = [x.split('_')[1] for x in os.listdir(self.path_fhiaims_species)] + basis_set_index = [basis_set_all.index(x) for x in basis_set_all if x in self.species] + + with open(os.path.join(vib_dir, 'control.in'), 'a') as f: + f.write("xc pbesol\n") + f.write("spin none\n") + f.write("relativistic atomic_zora scalar\n") + f.write("charge 0.\n\n") + f.write("# SCF convergence\n") + f.write("occupation_type gaussian 0.01\n") + f.write("mixer pulay\n") + f.write("n_max_pulay 10\n") + f.write("charge_mix_param 0.5\n") + f.write("sc_accuracy_rho 1E-5\n") + f.write("sc_accuracy_eev 1E-3\n") + f.write("sc_accuracy_etot 1E-6\n") + f.write("sc_accuracy_forces 1E-4\n") + f.write("sc_iter_limit 1500\n\n") + + for i in basis_set_index: + with open(basis_set_files[i], 'r') as ff: + lines = ff.read() + f.write(lines) + f.write('\n') + return None + +if __name__ == "__main__": + vib = aimsvibcalc() + os.mkdir('vibration') + vib.vib_calc_prep + current_dir_name = os.path.basename(os.getcwd()) + vib.make_job_submit(f'n{current_dir_name}') + os.chdir('vibration') + os.system('qsub submit.sh') + + + diff --git a/MLTTV_spliter.py b/MLTTV_spliter.py new file mode 100644 index 0000000..0fc84a2 --- /dev/null +++ b/MLTTV_spliter.py @@ -0,0 +1,37 @@ +from ase import io + +def split_xyz_file(input_file, train_file, valid_file, test_file): + with open(input_file, 'r') as infile: + train_out = open(train_file, 'w') + valid_out = open(valid_file, 'w') + test_out = open(test_file, 'w') + + while True: + # Read the header line containing the number of atoms + header = infile.readline() + if not header: + break # End of file + + num_atoms = int(header.strip()) + block_lines = [infile.readline() for _ in range(num_atoms + 1)] + + # Determine which file to write to based on the current index + i = infile.tell() # Get current position in file + if i % 5 < 3: + output_file = train_out + elif i % 5 == 3: + output_file = valid_out + else: + output_file = test_out + + # Write the block to the chosen file + output_file.write(header) + output_file.writelines(block_lines) + + train_out.close() + valid_out.close() + test_out.close() + +# Usage: +split_xyz_file('Training_set.xyz', 'Training_set_test.xyz', 'Validation_set_test.xyz', 'Testing_set_test.xyz') + diff --git a/MLTrainingDataGenerator.py b/MLTrainingDataGenerator.py new file mode 100644 index 0000000..0051810 --- /dev/null +++ b/MLTrainingDataGenerator.py @@ -0,0 +1,567 @@ + +""" +dev note: +work on breathing method""" + +''' +Author: Dong-Gi Kang +Prepare ML-IP data using FHI-aims output +Training data type: vibrational mode of a cluster + +[it retreive the vibrational mode cluster geometry and forces from single point calculation and +generates extended xyz format of Training_set.xyz: contains the total energy, atomic coordination, atomic forces] +The Training_set.xyz will be placed in FIT directory and each vibrational mode ext xyz files are generated in ext_xyz directory (The directories are automatically generated from the code) + +N.B. Change the UCL id, budget code, executable path for FHI-aims and fhi-aims species directory + + + +help: +python MLTrainingDataGenerator.py -h +(execute the file at the directory where the geometry.in, control.in, viration (dir) located) + + +1. python {code.py} --mode run --eigenvector="7 8 9 10" would grab 7th, 8th, 9th 10th (can selectively) then, modify the GM with the step_size (GM geometry + eigenvector * step_size) and prepare the individual directories and submit the single point calculations. + +2. python {code.py} --mode retrieve --eigenvector="7 8 9 10" would grab the generated data from the [1.] and make the ext xyz for each vibrational mode and store into the ext_xyz directory + +3. [python {code.py} --mode make_extxyz] would grab the all data from ext_xyz and make ext xyz format of Training_set.xyz in FIT directory + +4. if you want to trianing MACE or GAP ML-IP use MACE_lib.py or second_GAP.py +''' + +import os +import sys +import random +import numpy as np +import argparse +from itertools import groupby +from AppOutputExtractor.FHIaims.FHIaimsOutputExtractor import extractor + +class ML_train_generator(extractor): + + def __init__(self, app_version='22', tag=None): + + self.breathing_called = False + + #self.extractor = extractor() + #self.extractor.set_output_filepath(app_output) + #self.no_atoms = self.extractor.get_no_atoms() + #self.geometries = self.extractor.get_geometries(self.no_atoms) + #self.order = self.extractor.get_atom_order(self.no_atoms) + #ID = self.extractor.get_species(self.order) + #self.forces = self.extractor.get_forces(self.no_atoms) + #self.vib_eigvecs = self.extractor.get_vib_eigvec(self.no_atoms) + + self.ucl_id = 'uccatka' + self.job_time = '2:00:00' + self.job_name = 'test' + self.memory = '1' + self.cpu_core = '40' # for Young 40 core = 1 node + self.payment = 'Gold' + self.budgets = 'UCL_chemM_Woodley' + self.path_binary = '/home/uccatka/software/fhi-aims.221103/build/aims.221103.scalapack.mpi.x' + self.path_fhiaims_species = '/home/uccatka/software/fhi-aims.221103/species_defaults/defaults_2020/light' + self.step_size = 0.1 ##### STEP SIZE ##### + return None + + + def initiate(self): + app_output = './aims.out' + self.extractor = extractor() + self.extractor.set_output_filepath(app_output) + self.no_atoms = self.extractor.get_no_atoms() + self.geometries = self.extractor.get_geometries(self.no_atoms) + self.order = self.extractor.get_atom_order(self.no_atoms) + ID = self.extractor.get_species(self.order) + self.forces = self.extractor.get_forces(self.no_atoms) + self.vib_eigvecs = self.extractor.get_vib_eigvec(self.no_atoms) + + + def mod_xyz_w_vib(self): + ''' Modify LM geometries to the frames of vibrational mode frames ''' + Lambda = len(np.arange(-1, 1+self.step_size, self.step_size)) * self.no_atoms*3 + self.mod_sp = np.zeros((Lambda, self.no_atoms, 3)) + cnt = 0 + for i in range(self.no_atoms * 3): # 3N dimension + for numj, j in enumerate(np.arange(-1, 1+self.step_size, self.step_size)): # -1 to 1 in every step size + j = np.round(j, 2) + frame = self.geometries[-1] + self.vib_eigvecs[i] * j + self.mod_sp[cnt] = np.round(frame, 8) + cnt += 1 + self.mod_sp = np.reshape(self.mod_sp, (self.no_atoms*3, len(np.arange(-1, 1+self.step_size, self.step_size)), self.no_atoms, 3)) + return self.mod_sp + + + def mod_xyz_w_rand_pair_vib(self): + + list_eigvecs = list(range(6, self.no_atoms*3)) + random.shuffle(list_eigvecs) + pairs_eigvecs = [[list_eigvecs[i], list_eigvecs[i+1]] for i in range(0, len(list_eigvecs), 2)] + + #Lambda = len(np.arange(-1, 1+self.step_size, self.step_size)) * (self.no_atoms*3-6) # range of steps for all vib. mode, except E(3) + #self.mod_sp_pair = np.zeros((Lambda, self.no_atoms, 3)) + + self.mod_sp_pair = np.zeros((len(pairs_eigvecs), len(np.arange(-1, 1+self.step_size, self.step_size)), self.no_atoms, 3)) # range of steps for all vib. mode, except E(3) + + cnt = 0 + for numi, i in enumerate(pairs_eigvecs): + for numj, j in enumerate(np.arange(-1, 1+self.step_size, self.step_size)): + j = np.round(j, 2) + frame = self.geometries[-1] + (self.vib_eigvecs[i[0]]+self.vib_eigvecs[i[1]]) * j + #self.mod_sp_pair[cnt] = np.round(frame, 8) + self.mod_sp_pair[numi][numj] = np.round(frame, 8) + cnt += 1 + + self.mod_sp_pair = np.reshape(self.mod_sp_pair, (len(pairs_eigvecs), numj+1, self.no_atoms, 3)) + return self.mod_sp_pair + + + def breathing(self): + scale = np.arange(0.6, 1+self.step_size, self.step_size) + Lambda = len(scale) #* self.no_atoms*3 + self.mod_sp_breath = np.zeros((Lambda, self.no_atoms, 3)) + # shift the centre of mass of the structure to (0, 0, 0) + coord = self.geometries[0] + com = coord.sum(axis=0) + com = com / int(self.no_atoms) + coord_x = np.subtract(coord[:, 0], com[0], out=coord[:, 0]) + coord_y = np.subtract(coord[:, 1], com[1], out=coord[:, 1]) + coord_z = np.subtract(coord[:, 2], com[2], out=coord[:, 2]) + coord = list(zip(coord_x, coord_y, coord_z)) + coord = np.array(coord) + cnt = 0 + + for numj, j in enumerate(scale): + j = np.round(j, 2) + frame = coord * j + self.mod_sp_breath[cnt] = np.round(frame, 8) + cnt += 1 + + self.mod_sp_breath = np.reshape(self.mod_sp_breath, (len(scale), self.no_atoms, 3)) + self.breathing_called = True + return self.mod_sp_breath, scale + + + #@property + def geometry_for_sp(self, mod_sp): + ''' Convert the modified geometry (mod_xyz_w_vib) to {geometry.in} format for FHI-aims ''' + # vibrational modes + if not self.breathing_called: + print("@@@@@@@") + placer = np.full((self.no_atoms, 1), 'atom') + placer_species = np.reshape(self.order, (-1, 1)) + shape = np.shape(mod_sp) + self.for_sp = np.empty((shape[0], shape[1], self.no_atoms, 5), dtype=object) + for i in range(shape[0]): + for j in range(shape[1]): + form = np.concatenate((placer, mod_sp[i][j], placer_species), axis=1) + self.for_sp[i][j] = form + return self.for_sp, self.no_atoms + + # breathing mode + else: + print("*******") + placer_breath = np.full((self.no_atoms, 1), 'atom') + placer_species_breath = np.reshape(self.order, (-1, 1)) + shape_breath = np.shape(mod_sp) + self.for_sp = np.empty((shape_breath[0], shape_breath[1], 5), dtype=object) + for i in range(shape_breath[0]): + form = np.concatenate((placer_breath, mod_sp[i], placer_species_breath), axis=1) + self.for_sp[i] = form + return self.for_sp, self.no_atoms + + + @property + def xyz_from_opti(self): + ''' prepare training data from every SCF converged cycles of a optimisation ''' + train_xyz = 'xyz_from_opti.xyz' + exist = [x for x in os.listdir('./') if train_xyz in x] + if len(exist) != 0: + os.remove(exist[0]) + else: pass + for i in range(len(self.extractor.set_scf_blocks)): + self.energy = self.extractor.get_total_energy(i) + self.geometry = self.geometries[i] + self.force = self.forces[i] + xyz = np.round(np.concatenate((self.geometry, self.force), axis=1), 9) + xyz = np.concatenate((self.order, xyz), axis=1) + + with open('xyz_from_opti.xyz', 'a') as f: + f.write(f'{self.no_atoms}\n') + f.write(f'Properties=species:S:1:pos:R:3:forces:R:3 energy={self.energy} pbc="F F F"\n') + np.savetxt(f, xyz, fmt="%s", delimiter=" ") + print(f"total of {i+1} SCF converged structures are prepared in {train_xyz}") + + + def make_sp_control(self, path): + ''' Write {control.in} file ''' + basis_set_files = [os.path.join(self.path_fhiaims_species, x) for x in os.listdir(self.path_fhiaims_species)] + basis_set_all = [x.split('_')[1] for x in os.listdir(self.path_fhiaims_species)] + basis_set_index = [basis_set_all.index(x) for x in basis_set_all if x in ID] + + path = os.path.join(path, 'control.in') + with open(path, 'a') as f: + f.write("#\n") + f.write("xc pbesol\n") + f.write("spin none\n") + f.write("relativistic atomic_zora scalar\n") + f.write("charge 0.\n\n") + f.write("# SCF convergence\n") + f.write("occupation_type gaussian 0.01\n") + f.write("mixer pulay\n") + f.write("n_max_pulay 10\n") + f.write("charge_mix_param 0.5\n") + f.write("sc_accuracy_rho 1E-5\n") + f.write("sc_accuracy_eev 1E-3\n") + f.write("sc_accuracy_etot 1E-6\n") + f.write("sc_accuracy_forces 1E-4\n") + f.write("sc_iter_limit 1500\n") + f.write("# Relaxation\n\n") + #f.write("relax_geometry bfgs 1.e-3\n") + for i in basis_set_index: + with open(basis_set_files[i], 'r') as ff: + lines = ff.read() + f.write(lines) + f.write('\n') + return None + + + def make_job_submit(self, path): + ''' Write 'submit.sh' job script for SGE system ''' + _, last_part = os.path.split(path) + _, second_last_part = os.path.split(os.path.dirname(path)) + + # Combine the last two parts + last_two_parts = f"{second_last_part}_{last_part}" + + path = os.path.join(path, 'submit.sh') + with open(path, 'a') as f: + f.write("#!/bin/bash -l\n") + f.write('\n') + f.write("#$ -S /bin/bash\n") + f.write(f"#$ -l h_rt={self.job_time}\n") + f.write(f"#$ -l mem={self.memory}G\n") + f.write(f"#$ -N p{last_two_parts}\n") + f.write(f"#$ -pe mpi {self.cpu_core}\n") + f.write("#$ -cwd\n") + f.write("\n") + f.write(f"#$ -P {self.payment}\n") + f.write(f"#$ -A {self.budgets}\n") + + f.write("module load gerun\n") + f.write("module load userscripts\n") + f.write("module unload -f compilers mpi gcc-libs\n") + f.write("module load gcc-libs/4.9.2\n") + f.write("module unload -f compilers mpi\n") + f.write("module load beta-modules\n") + f.write("module load openblas/0.3.7-serial/gnu-4.9.2\n") + f.write("module load compilers/intel/2019/update5\n") + f.write("module load mpi/intel/2018/update3/intel\n") + + f.write("\n") + f.write("####$ -m e\n") + f.write(f"####$ -M {self.ucl_id}@ucl.ac.uk\n") + f.write("\n") + f.write(f"gerun {self.path_binary} > aims.out\n") + + + @staticmethod + def sorting_key(path): + parts = path.split('/') + second_key = int(parts[1]) if parts[1] != "breathing" else float('inf') + third_key = float(parts[2].split('_')[1]) # Consider lambda value regardless of the second part + return second_key, third_key + + def retrieve_results(self, eigenvectors): + print("---retrieve---") + eigvec_path = [os.path.join('sp', str(eigvec)) for eigvec in eigenvectors] + sp_path = [os.path.join(dirpath, fname) for dirpath in eigvec_path for fname in os.listdir(dirpath)] + lambda_path = [os.path.join(dirpath, fname) for dirpath in sp_path for fname in os.listdir(dirpath) if fname == 'aims.out'] + aims_out_path = sorted(lambda_path, key=self.sorting_key) + aims_out_path = [list(group) for key, group in groupby(aims_out_path, lambda x: x.split('/')[1])] + + if not os.path.exists('ext_xyz'): + os.mkdir('ext_xyz') + cnt = 0 + for numi, i in enumerate(aims_out_path): + filename = f"ext_xyz/ext_{i[0].split('/')[1]}_eigv.xyz" # + with open(filename, 'a') as f: # + for j in i: + ex = extractor() + ex.set_output_filepath(j) + #ex.set_scf_blocks + + no_atoms = ex.get_no_atoms() + geometries, atom_label = ex.get_sp_geometries(j) + forces = ex.get_sp_forces(no_atoms, j) + total_energy = ex.get_sp_total_energy(j) + #forces = np.round(forces, 8) + + coulomb_E, coulomb_F = self.coulomb_E_F(atom_label, geometries) + #print(j) + #print("structure") + #print(geometries) + #print("atomic force") + #print(forces) + #print("coulomb force") + #print(coulomb_F) + #print("coulomb energy") + #print(coulomb_E) + #print("total energy") + #print(total_energy) + #print(atom_label) + #print() + #print() + + # subtract coulomb energy and force + energy = total_energy - coulomb_E + forces = forces - coulomb_F + form = np.concatenate((ex.get_sp_atom_order(), geometries, forces), axis=1) + + f.write(str(no_atoms) + '\n') + f.write(f'Lattice="0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0" Properties=species:S:1:pos:R:3:forces:R:3 energy={energy} pbc="F F F"\n') + np.savetxt(f, form, fmt="%s", delimiter=" ") + + + def make_extxyz(self): + if not os.path.exists('FIT'): + os.mkdir('FIT') + else: pass + if os.path.exists('./FIT/Training_set.xyz'): + os.remove('./FIT/Training_set.xyz') + print("You may want to check the .xyz files in the FIT") + else: + with open('FIT/Training_set.xyz', 'a') as outfile: + filenames = [file for file in os.listdir('ext_xyz') if file.endswith('.xyz')] + sorted_filenames = sorted(filenames, key=lambda x: int(x.split('_')[1]) if x.split('_')[1] != 'breathing' else float('inf')) + for file in sorted_filenames: + print(file) + with open(os.path.join('ext_xyz', file), 'r') as infile: + for line in infile: + outfile.write(line) + + + def split_xyz_file(self, input_file, train_file, valid_file, test_file): + with open(input_file, 'r') as infile: + train_out = open(train_file, 'w') + valid_out = open(valid_file, 'w') + test_out = open(test_file, 'w') + + block_counter = 0 + line = infile.readline() + + while line: + if line.strip().isdigit(): + no_atoms = int(line.strip()) + block_lines = [line] + [infile.readline() for _ in range(no_atoms + 1)] # Read the block + + if block_counter % 5 < 3: + output_file = train_out + elif block_counter % 5 == 3: + output_file = valid_out + else: + output_file = test_out + + output_file.writelines(block_lines) + block_counter += 1 + + line = infile.readline() + + train_out.close() + valid_out.close() + test_out.close() + + + def coulomb_energy(self, r, cat_q, an_q): + return (cat_q * an_q) / r * 14.3996439067522 + + def coulomb_force(self, r, unit_r, cat_q, an_q): + return (cat_q * an_q) / r**2 * unit_r * 14.3996439067522 + + def coulomb_E_F(self, atom_label, structure): + self.charges = {"Al": 3, "F": -1} + coulomb_e = 0.0 + forces = np.zeros_like(structure) + + for i in range(len(structure)): + for j in range(i+1, len(structure)): + coord1 = structure[i] + coord2 = structure[j] + atom1 = atom_label[i][0] + atom2 = atom_label[j][0] + + r_vec = coord2 - coord1 + r = np.linalg.norm(r_vec) + unit_r = r_vec / r # unit vec + + # energy + energy_pair = self.coulomb_energy(r, self.charges[atom1], self.charges[atom2]) + coulomb_e += energy_pair + + # force + force_pair = self.coulomb_force(r, unit_r, self.charges[atom1], self.charges[atom2]) + + # add forces to atoms + forces[i] -= force_pair + forces[j] += force_pair + + return coulomb_e, forces + + + + + + +# executing the code using the class +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument("--eigenvector", type=str, help="A string of space-separated eigenvector indicies. For example, '7 8 9 10'") + parser.add_argument("--mode", type=str, choices=["run", "run_pair", "breath", "retrieve", "make_extxyz", "make_extxyz_"], default="run", help="Specify 'run' to execute the first part of the code, 'retrieve' to execute the second part of the code, or 'make_extxyz' to append all .xyz files into Training_set.xyz.") + args = parser.parse_args() + + ml = ML_train_generator() + step_size = ml.step_size ##### STEP SIZE ##### + + + # + # run + # + if args.mode == "run": + ml.initiate() + app_output = './aims.out' + + mod_sp = ml.mod_xyz_w_vib() # for each of vib. mode + sp_frame, no_atoms = ml.geometry_for_sp(mod_sp) + + shape = np.shape(sp_frame) + if args.eigenvector == 'all': + indicies = list(range(7, no_atoms*3+1)) + print("all eigenvectors without rotational and translational\n") + else: + indicies = list(map(int, args.eigenvector.split())) + + if not os.path.exists('sp'): + os.mkdir('sp') + else: pass + + for i in indicies: # Now we only iterate over the specified indicies + if not os.path.exists(os.path.join('sp', str(i))): + os.mkdir(f'sp/{str(i)}') + else: pass + + for numj, j in enumerate(np.arange(-1, 1+step_size, step_size)): + j = str(np.round(j, 2)) + os.mkdir(f'sp/{str(i)}/lambda_{j}') + + with open(f'sp/{i}/lambda_{j}/geometry.in', 'w') as f: + for row in sp_frame[i-1][numj]: + line = ' '.join(str(x) for x in row) + f.write(line + '\n') + ml.make_sp_control(f'sp/{i}/lambda_{j}') + ml.make_job_submit(f'sp/{i}/lambda_{j}') + os.chdir(f'sp/{i}/lambda_{j}') + os.system('qsub submit.sh') # submit jobs + os.chdir('../../../') + + # + # randomly pair up eigenvectors + # + elif args.mode == "run_pair": + ml.initiate() + app_output = './aims.out' + + mod_sp = ml.mod_xyz_w_rand_pair_vib() # randomly paired vib. mode + sp_frame, no_atoms = ml.geometry_for_sp(mod_sp) + + shape = np.shape(sp_frame) + if args.eigenvector == 'all': + + indicies = list(range(15)) + print("all paired eigenvectors without E(3), (rotational and translational)\n") + else: + indicies = list(map(int, args.eigenvector.split())) + + if not os.path.exists('sp'): + os.mkdir('sp') + else: pass + + for i in indicies: # Now we only iterate over the specified indicies + i = i+1 + if not os.path.exists(os.path.join('sp', str(i))): + os.mkdir(f'sp/{str(i)}_pair') + else: pass + + for numj, j in enumerate(np.arange(-1, 1+step_size, step_size)): + j = str(np.round(j, 2)) + os.mkdir(f'sp/{str(i)}_pair/lambda_{j}') + + with open(f'sp/{i}_pair/lambda_{j}/geometry.in', 'w') as f: + for row in sp_frame[i-1][numj]: + line = ' '.join(str(x) for x in row) + f.write(line + '\n') + ml.make_sp_control(f'sp/{i}_pair/lambda_{j}') + ml.make_job_submit(f'sp/{i}_pair/lambda_{j}') + os.chdir(f'sp/{i}_pair/lambda_{j}') + os.system('qsub submit.sh') # submit jobs + os.chdir('../../../') + + # + # preparen and run breathing mode single point calc + # + if args.mode == "breath": + ml.initiate() + if not os.path.exists('sp'): + os.mkdir('sp') + if not os.path.exists('sp/breathing'): + os.mkdir('sp/breathing') # for breathing mode + + mod_sp_breath, scale = ml.breathing() + #print(mod_sp_breath) + sp_frame, no_atoms = ml.geometry_for_sp(mod_sp_breath) + + # breathing + for numk, k in enumerate(scale): + k = str(np.round(k, 2)) + os.mkdir(f'sp/breathing/lambda_{k}') + + with open(f'sp/breathing/lambda_{k}/geometry.in', 'w') as f: + for row in sp_frame[numk]: + line = ' '.join(str(x) for x in row) + f.write(line + '\n') + ml.make_sp_control(f'sp/breathing/lambda_{k}') + ml.make_job_submit(f'sp/breathing/lambda_{k}') + os.chdir(f'sp/breathing/lambda_{k}') + os.system('qsub submit.sh') # submit job + os.chdir('../../../') + + # + # Collect data from single point calculated data + # + elif args.mode == "retrieve": + ml.initiate() + no_atoms = ml.no_atoms + if args.eigenvector == 'all': + indicies = list(range(7, no_atoms*3+1)) + indicies.append('breathing') + print(indicies) + print("all eigenvectors without rotational and translational\n") + else: + indicies = list(args.eigenvector.split()) + indicies = [int(x) if x.isdigit() else x for x in indicies] + print(indicies) + ml.retrieve_results(indicies) + + # + # make training data and split to train, test, valid data + # + elif args.mode == "make_extxyz": + ml.make_extxyz() + print("splitting training, test, validation data") + ml.split_xyz_file('./FIT/Training_set.xyz', './FIT/Training_set_test.xyz', './FIT/Validation_set_test.xyz', './FIT/Testing_set_test.xyz') + + # dev + elif args.mode == "make_extxyz_": + ml.split_xyz_file('./Training_set.xyz', './Training_set_test.xyz', './Validation_set_test.xyz', './Testing_set_test.xyz') + + diff --git a/README.md b/README.md index 11070e3..5aa068c 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,5 @@ # MAOA -#### Contributor: Dr. Woongkyu Jee, Dong-Gi Kang +#### Contributors: Dr. Woongkyu Jee (collect data from DFT optimisation output), Dong-Gi Kang, PhD candidate (the vibration, ML-IP workflow)
* * * #### This repository contains scripts that can take essential data from atomic simulation calculations
### Software that can work with: FHI-aims, GULP @@ -9,14 +9,20 @@ It has functions to take data from output file:
✅ Forces
✅ Energies
✅ Eigenvector of vibrational mode
- - +✅ Vibration calculation
+✅ Single point calculation of the frames of vibrational modes
+✅ Single point calculation of breathing mode frames
+✅ Preparing extended xyz file format (total energy, geometry, atomic forces)
+✅ Train, Test, Valid split
+✅ MACE (ML-IP): training
+✅ MACE (ML-IP): optimisation/single point calculation using the trained model
+✅ MACE (ML-IP): vibration calculation using the trained model
+✅ MACE (ML-IP): nudge elastic band calculation using the trained model
+✅ MACE (ML-IP): Plot dimer interaction energy with the interatomic distances of trained data in .html format (plotly)
+
☉ Upcoming function:
- **1.** Preparing extended xyz format - **2.** Take top _n_ many local minima of KLMC output (top_structure) to make FHI-aims calculation directories (geometry.in, control.in, submit.sh)
- **3.** preparing vibration calculation, and submit the job (SGE)
+ **1.** N/A + -## Log
-✅ May 2023: geometry, forces, energy, eigenvector of vibration mode.
diff --git a/SystemInfo.py b/SystemInfo.py new file mode 100644 index 0000000..185fbba --- /dev/null +++ b/SystemInfo.py @@ -0,0 +1,29 @@ +# + +from datetime import datetime + +def time_diff_in_seconds(time_str1, time_str2): + # Convert the time strings to datetime objects + dt1 = datetime.strptime(time_str1, "Date : %Y%m%d, Time : %H%M%S.%f") + dt2 = datetime.strptime(time_str2, "Date : %Y%m%d, Time : %H%M%S.%f") + + # Calculate the time difference in seconds + time_diff = abs((dt2 - dt1).total_seconds()) + + # Return the time difference in seconds + return time_diff + + +if __name__ == '__main__': + + time_str1 = "Date : 20230426, Time : 002047.673" + time_str2 = "Date : 20230426, Time : 001834.244" + + time_diff = time_diff_in_seconds(time_str1, time_str2) + + print("Time difference in seconds:", time_diff) + + + + + diff --git a/for_MACE/MACE_lib.py b/for_MACE/MACE_lib.py new file mode 100755 index 0000000..77bdd30 --- /dev/null +++ b/for_MACE/MACE_lib.py @@ -0,0 +1,447 @@ +import os +import subprocess +import argparse +from ase import io +from ase.optimize import BFGS +from ase.vibrations import Vibrations +from mace.calculators import MACECalculator +import numpy as np +from ase import Atoms +import plotly.graph_objects as go + + +""" +[For training only:] +python script_name.py --model_path "MACE_n6_test_model" --training_data_path "./Training_set_test.xyz" --testing_data_path "./Testing_set_test.xyz" --validation_data_path "./Validation_set_test.xyz" --radius 6.0 --epochs 100 --device_type "cuda" --layers "128x0e+128x1o+128x2e" + + +[For training and then single point calculation, you can use:] +python script_name.py --model_path "MACE_n6_test_model" --training_data_path "./Training_set_test.xyz" --testing_data_path "./Testing_set_test.xyz" --validation_data_path "./Validation_set_test.xyz" --radius 6.0 --epochs 100 --device_type "cuda" --layers "128x0e+128x1o+128x2e" --target_stru "path_to_structure_file" --model_path "path_to_model_file" --opt_struc_path "path_to_optimized_structure_file" + + +[For training and then optimization, you can use:] +python script_name.py --model_path "MACE_n6_test_model" --training_data_path "./Training_set_test.xyz" --testing_data_path "./Testing_set_test.xyz" --validation_data_path "./Validation_set_test.xyz" --radius 6.0 --epochs 100 --device_type "cuda" --layers "128x0e+128x1o+128x2e" --target_stru "path_to_structure_file" --model_path "path_to_model_file" --opt_struc_path "path_to_optimized_structure_file" --optimize + + +[For only single point calculation (without training):] +python script_name.py --target_stru "path_to_structure_file" --model_path "path_to_model_file" --opt_struc_path "path_to_optimized_structure_file" --skip_training + + +[For only optimization (without training):] +python script_name.py --target_stru "path_to_structure_file" --model_path "path_to_model_file" --opt_struc_path "path_to_optimized_structure_file" --optimize --skip_training +""" + + +class MACE: + def __init__(self, args): + self.args = args + + def MACE_training(self): + layers = 'x'.join(self.args.layers.split(' ')) # remove spaces around '+' + subprocess.check_call([ + "python", "/home/uccatka/software/mace/scripts/run_train.py", + "--name", self.args.model_path, + "--train_file", self.args.training_data_path, + "--test_file", self.args.testing_data_path, + "--valid_file", self.args.validation_data_path, + "--config_type_weights", '{"Default":1.0}', + "--model", "MACE", + #"--E0s", "{9:0.000, 13:0.000}", # for the only IP data + #"--E0s", "{9:-2707.428895973, 13:-6596.914328816}", # for the PBEsol data + "--E0s", "{9:-2711.537676517, 13:-6543.933824960}", + "--hidden_irreps", layers, + "--r_max", str(self.args.radius), + "--batch_size", str(self.args.batch_size), + "--max_num_epochs", str(self.args.epochs), + "--swa", + "--start_swa", "10", + "--ema", + "--ema_decay", "0.99", + "--amsgrad", + "--restart_latest", + "--device", self.args.device_type + ]) + + def calculation(self): + if not (self.args.target_stru and self.args.model_path): + print("Target structure and model path must be specified for calculation") + return + + if self.args.optimize: + initial = io.read(self.args.target_stru) + initial.set_calculator(MACECalculator(model_path=f"{self.args.model_path}.model", device=self.args.device_type)) + self.optimize(initial, self.args.opt_struc_path) + else: + self.single_point(self.args.target_stru, f"{self.args.model_path}.model", self.args.device_type) + + if self.args.vibration == 'y': + print("\n" * 2) + opt_struc_vib = io.read(self.args.opt_struc_path) + vib = Vibrations(opt_struc_vib) + vib.run() + vib.summary() + + + def optimize(self, initial, opt_struc): + dyn = BFGS(initial, trajectory=f'{opt_struc}_dummy.traj', logfile=f'{opt_struc}.log') + dyn.run(fmax=0.001) + #io.write(f"{opt_struc}.xyz", initial, format="xyz") + + read_traj = io.read(f'{opt_struc}_dummy.traj', index=":") + io.write(f'{opt_struc}_traj.xyz', read_traj, format='extxyz') + opt_trajectory = f'{opt_struc}_traj.xyz' + + read_opt = io.read(f'{opt_struc}_dummy.traj', index="-1:") + io.write(f'{opt_struc}.xyz', read_opt, format='extxyz') + return opt_struc, opt_trajectory + + + def single_point(self, target_stru, model_path, device_type): + structure = io.read(target_stru) + calculator = MACECalculator(model_path=f"{model_path}.model", device=device_type) + structure.set_calculator(calculator) + energy = structure.get_potential_energy() + forces = structure.get_forces() + return energy, forces + + + def dimer_curve(self, model_path, device_type, atom1='Al', atom2='F', distance_range=(0.0, 5.0), num_points=51): + # List of distances to calculate energies for + distances = np.linspace(*distance_range, num_points) + energies_cat_an = [] + energies_an_an = [] + energies_cat_cat = [] + BM_values = [] + buck4_values = [] + print("single point calculation of dimers at range of distance") + for d in distances: + d = round(d, 2) + print(f'{d} Ang') + # cation-anion interaction (MACE) + cat_an_dimer = Atoms(f'{atom1}{atom2}', positions=[(0, 0, 0), (0, 0, d)]) + calculator = MACECalculator(model_path=f"{model_path}.model", device=device_type) + cat_an_dimer.set_calculator(calculator) + energy_cat_an = cat_an_dimer.get_potential_energy() + energies_cat_an.append(energy_cat_an) + + # anion-anion interaction (MACE) + an_an_dimer = Atoms(f'{atom2}{atom2}', positions=[(0, 0, 0), (0, 0, d)]) + #calculator = MACECalculator(model_path=f"{model_path}.model", device=device_type) + an_an_dimer.set_calculator(calculator) + energy_an_an = an_an_dimer.get_potential_energy() + energies_an_an.append(energy_an_an) + + # cation-cation interaction (MACE) + cat_cat_dimer = Atoms(f'{atom1}{atom1}', positions=[(0, 0, 0), (0, 0, d)]) + #calculator = MACECalculator(model_path=f"{model_path}.model", device=device_type) + cat_cat_dimer.set_calculator(calculator) + energy_cat_cat = cat_cat_dimer.get_potential_energy() + energies_cat_cat.append(energy_cat_cat) + + # Calculate Al-F BM value for the current distance + BM_value = self.BM(d) + BM_values.append(BM_value) + + buck4_value = self.buck4(d) + buck4_values.append(buck4_value) + + # + # Getting interatomic distances + print("Retrieve training data") + with open('Training_set_test.xyz', 'r') as f: + full_lines = f.readlines() + + no_atoms = full_lines[0] + check_continue, cluster_set, clusters, ID, ID_set = [], [], [], [], [] + for numi, i in enumerate(full_lines): + if len(i) > 10 and "Properties" not in i: + check_continue.append(numi) + if numi - check_continue[-1] == 0: + clusters.append(i.split()[1:4]) + ID.append(i.split()[0]) + else: pass + + else: + if len(clusters) != 0: + clusters = np.array(clusters).astype(float) + cluster_set.append(clusters) # make nested list + ID_set.append(ID) # same here + else: pass + ID, clusters = [], [] + + cluster_set = np.array(cluster_set[:-1]) + + cat_cat_dist, an_an_dist, cat_an_dist = [], [], [] + print("Calculating ineteratomic distances") + for i in range(len(cluster_set)): + (npairs_all, npairs_cat_cat, npairs_an_an, npairs_cat_an, all_dist, c_c_dist, a_a_dist, c_a_dist) = \ + self.RDF(no_atoms, cluster_set[i], ID_set[i]) + + cat_cat_dist += c_c_dist # concatenate the list + an_an_dist += a_a_dist + cat_an_dist += c_a_dist + + + # + # PLOT + # + print("Start plotting") + BM_color = "rgb(10, 120, 24)" + + fig = go.FigureWidget() + + energy_trace_cat_an = fig.add_scatter( + x=distances, + y=energies_cat_an, + mode='lines', + name='MACE: cation - anion', + line=dict(shape='linear', color=BM_color) + ) + + energy_trace_an_an = fig.add_scatter( + x=distances, + y=energies_an_an, + mode='lines', + name='MACE: anion - anion', + line=dict(shape="linear", color="firebrick") + ) + + energy_trace_cat_cat = fig.add_scatter( + x=distances, + y=energies_cat_cat, + mode='lines', + name='MACE: cation - cation', + line=dict(shape="linear", color="blue") + ) + + BM_trace = fig.add_scatter( + x=distances, + y=BM_values, + mode='lines', + name='Al-F Born-Mayer', + line=dict(color=BM_color, dash="dot") + ) + + buck4_trace = fig.add_scatter( + x=distances, + y=buck4_values, + mode='lines', + name='F-F Buck4', + line=dict(color="firebrick", dash="dot") + ) + + + #### upper panel (hover) + cat_an_dist_trace = fig.add_histogram( + x=cat_an_dist, #all_het_dist, + xbins=dict(start=0, end=6, size=0.005), + marker_color=BM_color, + name="cat-an dist", + yaxis="y2") + + cat_cat_dist_trace = fig.add_histogram( + x=cat_cat_dist, #all_homo_dist, + xbins=dict(start=0, end=6, size=0.005), + marker_color="blue", + name="cat-cat dist", + yaxis="y2") + + an_an_dist_trace = fig.add_histogram( + x=an_an_dist, + xbins=dict(start=0, end=6, size=0.005), + marker_color="firebrick", + name="an-an dist", + yaxis="y2") + + + # + ### Settings/layout ### + # + fig.layout = dict( + xaxis=dict( + domain=[0, 0.8], + range=[0,6.0], + showgrid=False, + zeroline=False, + title="Interatomic distance / Å"), + yaxis=dict( + domain=[0, 0.8], + range=[-20, 50], + showgrid=False, + zeroline=True, + title="Potential energy / eV"), + legend=dict( + x=0.85, + y=1.0, + ), + margin=dict(l=80, r=80, t=80, b=80), + width=1400, + height=800, + hovermode="closest", + bargap=0.8, + font=dict(size=20), + xaxis2=dict( + domain=[0.85, 1], + showgrid=False, + zeroline=False), + + # hover plot + yaxis2=dict(domain=[0.85, 1], showgrid=False, zeroline=False, title="Count") + ) + + + fig.write_html('./dimer_curve.html') + print("A dimer curve figure is saved") + fig.show() + + + def RDF(self, no_of_atoms, coord, ID): + # Calculate the interatomic disntaces + all_dist = [] + c_a_dist = [] + c_c_dist = [] + a_a_dist = [] + npairs_all = 0 + npairs_cat_an = 0 + npairs_cat_cat = 0 + npairs_an_an = 0 + + all_dup_filter = [] + cat_an_dup_filter = [] + cat_cat_dup_filter = [] + an_an_dup_filter = [] + for i in range(len(ID)): + for j in range(i+1, len(ID)): + npairs_all += 1 + distance = np.linalg.norm(coord[i, :] - coord[j, :]) + all_dist.append(distance) + + # Interatomic distance between hetero species (cat-an) + if ID[i] != ID[j] and (str(i)+str(j) not in cat_an_dup_filter): + npairs_cat_an += 1 + distance = np.round(np.linalg.norm(coord[i, :] - coord[j, :]), 9) + c_a_dist.append(distance) + cat_an_dup_filter.append(str(i)+str(j)) + + # Interatomic distance between homo species + anions_list = ['F', 'Cl', 'Br'] # temporal + if ID[i] in anions_list: + if ID[i] == ID[j] and i != j and (str(j)+str(i) not in an_an_dup_filter): + npairs_an_an += 1 + distance = np.round(np.linalg.norm(coord[i,:] - coord[j, :]), 9) + a_a_dist.append(distance) + an_an_dup_filter.append(str(i)+str(j)) + + else: + if ID[i] == ID[j] and i != j and (str(j) + str(i) not in cat_cat_dup_filter): + npairs_cat_cat += 1 + distance = np.round(np.linalg.norm(coord[i,:] - coord[j, :]), 9) + c_c_dist.append(distance) + cat_cat_dup_filter.append(str(i) + str(j)) + + return (npairs_all, npairs_cat_cat, npairs_an_an, npairs_cat_an, all_dist, c_c_dist, a_a_dist, c_a_dist) + + + # + # analytical potentials + # + def BM(self, x): + return 3760 * np.exp(-x / 0.222) + + + def buck4(self, x): # 2.73154 Å F-F distance + if x.all() < 2.0: + return 1127.7 * np.exp(-x / 0.2753) + elif 2.0 <= x.all() < 2.726: + return ( + -3.976 * x**5 + + 49.0486 * x**4 + - 241.8573 * x**3 + + 597.2668 * x**2 + - 741.117 * x + + 371.2706 + ) + elif 2.726 <= x.all() < 3.031: + return -0.361 * x**3 + 3.2362 * x**2 - 9.6271 * x + 9.4816 + elif x.all() >= 3.031: + return -15.83 / x**6 + + + def coulomb(self, r, cat_q, an_q): + return (cat_q * an_q) / r * 14.3996439067522 + + def coulomb_force(self, r, unit_r, cat_q, an_q): + return (cat_q * an_q) / r**2 * unit_r * 14.3996439067522 + + def coulomb_energy(self, structure): + coulomb_e = 0.0 + forces = np.zeros_like(self.positions) + + for i in range(len(self.positions)): + for j in range(i+1, len(self.positions)): + coord1 = self.positions[i] + coord2 = self.positions[j] + atom1 = self.species[i] + atom2 = self.species[j] + + # Calculate the distance and unit vector between the two atoms + r_vec = coord2 - coord1 + r = np.linalg.norm(r_vec) + unit_r = r_vec / r + + # Calculate the Coulomb potential energy between the pair + energy_pair = self.coulomb_energy(r, self.charges[atom1], self.charges[atom2]) + + # Calculate the Coulomb force between the pair + force_pair = self.coulomb_force(r, unit_r, self.charges[atom1], self.charges[atom2]) + + # Add to the total energy + coulomb_e += energy_pair + + # Add forces to atoms + forces[i] -= force_pair + forces[j] += force_pair + + return coulomb_e, forces + + + + + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description='Training model') + parser.add_argument('--model_path', default="MACE_model", help='Name of the model') + parser.add_argument('--training_data_path', default="./Training_set_test.xyz", help='Path to training data') + parser.add_argument('--testing_data_path', default="./Testing_set_test.xyz", help='Path to testing data') + parser.add_argument('--validation_data_path', default="./Validation_set_test.xyz", help='Path to validation data') + parser.add_argument('--radius', type=float, default=6.0, help='Radius (default = 6.0)') + parser.add_argument('--batch_size', type=int, default=5, help='batch size (default = 5)') + parser.add_argument('--epochs', type=int, default=100, help='Number of epochs') + parser.add_argument('--device_type', default="gpu", help='Type of device') + parser.add_argument('--layers', default="128x0e+128x1o+128x2e", help='Layers') + + parser.add_argument('--target_stru', type=str, help="Path to the structure file") + parser.add_argument('--opt_struc_path', type=str, help="Path to the optimized structure file") + parser.add_argument('--vibration', type=str, default='n', help="Perform vibration? (y/n)") + parser.add_argument('--optimize', action='store_true', help="Perform structure optimization") + + parser.add_argument('--dimer', action='store_true', help="Calculate a dimer curve") + + parser.add_argument('--skip_training', action='store_true', help="Skip training and perform only calculation") + + args = parser.parse_args() + + os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE" # environment variable setting + + mace = MACE(args) + if not args.skip_training: + mace.MACE_training() + if args.target_stru and args.model_path: + mace.calculation() + + if args.dimer: + if args.model_path: + mace.dimer_curve(args.model_path, args.device_type) + else: + print("Model path must be specified to calculate a dimer curve") + diff --git a/for_MACE/MACE_train_test_val_split.py b/for_MACE/MACE_train_test_val_split.py new file mode 100755 index 0000000..6270501 --- /dev/null +++ b/for_MACE/MACE_train_test_val_split.py @@ -0,0 +1,49 @@ +import sys +import os +import subprocess +import numpy as np +from ase import io + +# required module for MACE training + +#module load amd-modules +#module load Python/3.10.4-GCCcore-11.3.0 +#module load foss/2022a CUDA/11.7.1 NCCL/2.12.12-GCCcore-11.3.0-CUDA-11.7.1 +#source ~/venv/mace/bin/activate + +training_data_path = sys.argv[1] + +#orig_train_file = io.read(training_data_path, index=":") +#id_1 = int(len(orig_train_file)*0.8) +#id_2 = int(len(orig_train_file)*0.9) +#train, test, valid = np.split(orig_train_file, [id_1, id_2]) +# +#train = io.read(training_data_path, index=f":{len(train)}") +#test = io.read(training_data_path, index=f"{len(train)}:{len(train)+len(test)}") +#valid = io.read(training_data_path, index=f"{len(train)+len(test)}:") +# +#train = io.write('./Training_set_test.xyz', train, format='extxyz') +#test = io.write( './Testing_set_test.xyz', test, format='extxyz') +#valid = io.write('./Validation_set_test.xyz', valid, format='extxyz') +##test_run_sample = io.write('./FIT/tesing_sample.xyz', sample, format='extxyz') + +orig_train_file = io.read(training_data_path, index=":") + +train = [] +valid = [] +test = [] + +for i, data in enumerate(orig_train_file): + if i % 5 < 3: # First three items in every set of five + train.append(data) + elif i % 5 == 3: # Fourth item in every set of five + valid.append(data) + elif i % 5 == 4: # Fifth item in every set of five + test.append(data) + +# Write the data to the respective files +io.write('./Training_set_test.xyz', train, format='extxyz') +io.write('./Validation_set_test.xyz', valid, format='extxyz') +io.write('./Testing_set_test.xyz', test, format='extxyz') + + diff --git a/for_MACE/geoopt.py b/for_MACE/geoopt.py new file mode 100755 index 0000000..33d9f30 --- /dev/null +++ b/for_MACE/geoopt.py @@ -0,0 +1,155 @@ +import sys +import os +import argparse +import subprocess +from ase import io +from ase.optimize import BFGS +from ase.vibrations import Vibrations +from mace.calculators import MACECalculator + + +class MACE: + def __init__(self): + pass + + def training_MACE(self): + ## Define the command-line arguments + #parser = argparse.ArgumentParser(description="MACE Training Instructions") + #parser.add_argument('out_model_name', type=str, help="Output model name") + #parser.add_argument('training_data_path', type=str, help="Path to the training data file") + #parser.add_argument('testing_data_path', type=str, help="Path to the testing data file") + #parser.add_argument('radius', type=float, help="Radius") + #parser.add_argument('epochs', type=int, help="Number of epochs") + #parser.add_argument('device_type', type=str, default='cpu', help="Device type (default: cpu)") + #parser.add_argument('layers', type=str, nargs='?', default='128x0e + 128x1o + 129x2e', + # help="NN layers (default: '128x0e + 128x1o + 129x2e')") + + ## Parse the command-line arguments + #args = parser.parse_args() + + ## Access the parsed arguments + #out_model_name = args.out_model_name + #training_data_path = args.training_data_path + #testing_data_path = args.testing_data_path + #radius = args.radius + #epochs = args.epochs + #device_type = args.device_type + #layers = args.layers + + # Process the layers argument + if layers is None: + layers = '128x0e + 128x1o + 129x2e' + + # [TEMPORARY] Read the training, testing, and validation data + train = io.read(training_data_path, index=":30") + test = io.read(training_data_path, index="-10:") + valid = io.read(training_data_path, index="30:35") + + # Write the training, testing, and validation data to files + io.write('./FIT/Training_set_test.xyz', train, format='extxyz') + io.write('./FIT/Testing_set_test.xyz', test, format='extxyz') + io.write('./FIT/Validation_set_test.xyz', valid, format='extxyz') + + + def optimize(self, initial_placer, opt_struc): + dyn = BFGS(initial_placer, trajectory=f'{opt_struc}_dummy.traj') + dyn.run(fmax=0.001) + io.write(opt_struc, initial_placer, format="xyz") + + read_traj = io.read(f'{opt_struc}_dummy.traj', index=":") + io.write(f'{opt_struc}_traj.xyz', read_traj) + opt_trajectory = f'{opt_struc}_traj.xyz' + return opt_struc, opt_trajectory + + + def single_point(self, target_stru, model_path, device): + structure = io.read(target_stru) + calculator = MACECalculator(model_path=model_path, device=device) + structure.set_calculator(calculator) + energy = structure.get_potential_energy() + forces = structure.get_forces() + return energy, forces + + + def calculation(self): + #parser = argparse.ArgumentParser(description="Geometric Optimization Instructions") + #parser.add_argument('target_stru', type=str, help="Path to the structure file") + #parser.add_argument('model_path', type=str, help="Path to the model file") + #parser.add_argument('opt_struc_path', type=str, help="Path to the optimized structure file") + #parser.add_argument('vibration', type=str, default='n', help="Perform vibration? (y/n)") + #parser.add_argument('device', type=str, default='gpu', help="Device to use for computation (default: gpu)") + #parser.add_argument('--optimize', action='store_true', help="Perform structure optimization") + + #args = parser.parse_args() + + #target_stru = args.target_stru + #model_path = args.model_path + #opt_struc = args.opt_struc_path + #vibration = args.vibration.lower() + #device = args.device + + if args.optimize: + initial = io.read(target_stru) + initial.set_calculator(MACECalculator(model_path=model_path, device=device)) + self.optimize(initial, opt_struc) + else: + self.single_point(target_stru, model_path, device) + + if vibration == 'y': + print("\n" * 2) + opt_struc_vib = io.read(opt_struc) + vib = Vibrations(opt_struc_vib) + vib.run() + vib.summary() + + + +if __name__ == '__main__': + # Define the command-line arguments + parser = argparse.ArgumentParser(description="MACE Training Instructions") + parser.add_argument('out_model_name', type=str, help="Output model name") + parser.add_argument('training_data_path', type=str, help="Path to the training data file") + parser.add_argument('testing_data_path', type=str, help="Path to the testing data file") + parser.add_argument('radius', type=float, help="Radius") + parser.add_argument('epochs', type=int, help="Number of epochs") + parser.add_argument('device_type', type=str, default='cpu', help="Device type (default: cpu)") + parser.add_argument('layers', type=str, nargs='?', default='128x0e + 128x1o + 129x2e', + help="NN layers (default: '128x0e + 128x1o + 129x2e')") + + # Parse the command-line arguments + args = parser.parse_args() + + # Access the parsed arguments + out_model_name = args.out_model_name + training_data_path = args.training_data_path + testing_data_path = args.testing_data_path + radius = args.radius + epochs = args.epochs + device_type = args.device_type + layers = args.layers + + # Create an instance of the MACE class + mace = MACE() + + # Perform the MACE training + mace.training_MACE() + + + + parser = argparse.ArgumentParser(description="Geometric Optimization Instructions") + parser.add_argument('target_stru', type=str, help="Path to the structure file") + parser.add_argument('model_path', type=str, help="Path to the model file") + parser.add_argument('opt_struc_path', type=str, help="Path to the optimized structure file") + parser.add_argument('vibration', type=str, default='n', help="Perform vibration? (y/n)") + parser.add_argument('device', type=str, default='gpu', help="Device to use for computation (default: gpu)") + parser.add_argument('--optimize', action='store_true', help="Perform structure optimization") + + args = parser.parse_args() + + target_stru = args.target_stru + model_path = args.model_path + opt_struc = args.opt_struc_path + vibration = args.vibration.lower() + device = args.device + # Perform the calculation (optimization or single-point) + mace.calculation() diff --git a/for_MACE/geoopt_2.py b/for_MACE/geoopt_2.py new file mode 100644 index 0000000..edf746b --- /dev/null +++ b/for_MACE/geoopt_2.py @@ -0,0 +1,136 @@ +import os +import sys +import argparse +from ase import io +from ase.optimize import MDMin +from ase.dyneb import DyNEB +from ase.optimize import BFGS +from ase.vibrations import Vibrations +from ase.build import bulk +from ase.calculators.emt import EMT +from ase.phonons import Phonons +from mace.calculators import MACECalculator + + +class MACE: + def __init__(self): + pass + + def training_MACE(self, training_data_path, layers): + # [TEMPORARY] Read the training, testing, and validation data + train = io.read(training_data_path, index=":30") + test = io.read(training_data_path, index="-10:") + valid = io.read(training_data_path, index="30:35") + + # Write the training, testing, and validation data to files + io.write('./FIT/Training_set_test.xyz', train, format='extxyz') + io.write('./FIT/Testing_set_test.xyz', test, format='extxyz') + io.write('./FIT/Validation_set_test.xyz', valid, format='extxyz') + + def optimize(self, initial_placer, opt_struc): + dyn = BFGS(initial_placer, trajectory=f'{opt_struc}_dummy.traj') + dyn.run(fmax=0.001) + io.write(opt_struc, initial_placer, format="xyz") + + read_traj = io.read(f'{opt_struc}_dummy.traj', index=":") + io.write(f'{opt_struc}_traj.xyz', read_traj) + opt_trajectory = f'{opt_struc}_traj.xyz' + return opt_struc, opt_trajectory + + def single_point(self, target_stru, model_path, device): + structure = io.read(target_stru) + calculator = MACECalculator(model_path=model_path, device=device) + structure.set_calculator(calculator) + energy = structure.get_potential_energy() + forces = structure.get_forces() + return energy, forces + + def calculation(self, target_stru, model_path, opt_struc, vibration, device, optimize): + if optimize: + initial = io.read(target_stru) + initial.set_calculator(MACECalculator(model_path=model_path, device=device)) + self.optimize(initial, opt_struc) + else: + self.single_point(target_stru, model_path, device) + + if vibration == 'y': + print("\n" * 2) + opt_struc_vib = io.read(opt_struc) + vib = Vibrations(opt_struc_vib) + vib.run() + vib.summary() + + +if __name__ == '__main__': + + # Create the top-level parser + parser = argparse.ArgumentParser(description="MACE [Training] and [Geometric Optimization Instructions]") + subparsers = parser.add_subparsers(dest='command', help='sub-command help') + + # Create the sub-parser for the training command + parser_train = subparsers.add_parser('training', help='Training-related instructions') + parser_train.add_argument('--out_model_name', type=str, help="Output model name") + parser_train.add_argument('--training_data_path', type=str, default="./FIT/Training_set_test.xyz", help="Path to the training data file") + parser_train.add_argument('--testing_data_path', type=str, default="./FIT/Testing_set_test.xyz", help="Path to the testing data file") + parser_train.add_argument('--validation_data_path', type=str, default="./FIT/Validation_set_test.xyz", help="Path to validation data (default=./FIT/Validation_set_test.xyz)") + parser_train.add_argument('--radius', type=float, default="6.0", help="Radius") + parser_train.add_argument('--epochs', type=int, default="100", help="Number of epochs") + parser_train.add_argument('--device_type', type=str, default='cpu', help="Device type (default: cpu)") + parser_train.add_argument('--layers', type=str, nargs='?', default='128x0e + 128x1o + 129x2e', help="NN layers (default: '128x0e + 128x1o + 129x2e')") + + # Create the sub-parser for the calculation command + parser_calc = subparsers.add_parser('calculation', help='Calculation-related instructions') + parser_calc.add_argument('--optimize', action='store_true', help="Perform structure optimization") + parser_calc.add_argument('--target_stru', type=str, help="Path to the structure file") + parser_calc.add_argument('--model_path', type=str, help="Path to the model file") + parser_calc.add_argument('--opt_struc_path', type=str, help="Path to the optimized structure file") + parser_calc.add_argument('--vibration', type=str, default='n', help="Perform vibration? (y/n)") + parser_calc.add_argument('--device', type=str, default='gpu', help="Device to use for computation (default: gpu)") + + args = parser.parse_args() + + command = args.command + + if command == 'training': + # Access the parsed training arguments + out_model_name = args.out_model_name + training_data_path = args.training_data_path + testing_data_path = args.testing_data_path + radius = args.radius + epochs = args.epochs + device_type = args.device_type + layers = args.layers + + # Perform the training-related instructions... + + print("Instructions for Training:") + print(f"- Output model name: {out_model_name}") + print(f"- Training data path: {training_data_path}") + print(f"- Testing data path: {testing_data_path}") + print(f"- Validation data path: {validation_data_path}") + print(f"- Radius: {radius}") + print(f"- Epochs: {epochs}") + print(f"- Device type: {device_type}") + print(f"- Layers: {layers}") + + elif command == 'calculation': + # Access the parsed calculation arguments + optimize = args.optimize + target_stru = args.target_stru + model_path = args.model_path + opt_struc = args.opt_struc_path + vibration = args.vibration.lower() + device = args.device + + # Perform the calculation-related instructions... + + print("Instructions for Calculation:") + print(f"- Optimize: {optimize}") + print(f"- Target structure path: {target_stru}") + print(f"- Model path: {model_path}") + print(f"- Optimized structure path: {opt_struc}") + print(f"- Vibration: {vibration}") + print(f"- Device: {device}") + + + diff --git a/for_MACE/geoopt_orig.py b/for_MACE/geoopt_orig.py new file mode 100755 index 0000000..34a68df --- /dev/null +++ b/for_MACE/geoopt_orig.py @@ -0,0 +1,68 @@ +import sys +import os +from ase import io +from ase.optimize import BFGS +from ase.vibrations import Vibrations +from mace.calculators import MACECalculator + +def optimize(initial_placer, opt_struc): + dyn = BFGS(initial_placer, trajectory='dummy.traj') + dyn.run(fmax=0.001) + io.write(opt_struc, initial_placer, format="extxyz") + + read_traj = io.read('dummy.traj', index=":") + io.write('traj.xyz', read_traj) + +def single_point(self, target_stru, model_path, device): + structure = io.read(target_stru) + calculator = MACECalculator(model_path=model_path, device=device) + structure.set_calculator(calculator) + energy = structure.get_potential_energy() + forces = structure.get_forces() + return energy, forces + + +if sys.argv[1] == '-h': + print() + print("Instruction:") + print("python geoopt.py {struc_path} {model_path} {opt_struc} {vibration} {DEVICE}") +else: pass + + +#try: +struc_path = sys.argv[1] +model_path = sys.argv[2] +opt_struc = sys.argv[3] +vibration = sys.argv[4].lower() + +if vibration == None: + vibration == 'n' +else: pass + +DEVICE = sys.argv[5] +if DEVICE == None: + DEVCIE = 'gpu' +else: pass + +initial = io.read(struc_path) +initial.set_calculator(MACECalculator(model_path=model_path, device=DEVICE)) +#optimize(initial, opt_struc) +energy, force = single_point(struc_path, model_path, device) +print(energy, forces) + + +if vibration == 'y': + print() + print() + opt_struc_vib = io.read(opt_struc) + vib = Vibrations(opt_struc_vib) + vib.run() + vib.summary() +else: pass + +#except IndexError or NameError: +# print("Error") +# print("python geoopt.py {struc path want to opt} {model path} {name for opt struc} {vibration (default n)} {device (default gpu)}") +# print() + + diff --git a/for_MACE/geoopt_orig_2.py b/for_MACE/geoopt_orig_2.py new file mode 100644 index 0000000..26d394b --- /dev/null +++ b/for_MACE/geoopt_orig_2.py @@ -0,0 +1,48 @@ +import os +import argparse +from ase import io +from ase.optimize import BFGS +from ase.vibrations import Vibrations +from mace.calculators import MACECalculator + +def optimize(initial_placer, opt_struc): + dyn = BFGS(initial_placer, trajectory='dummy.traj') + dyn.run(fmax=0.001) + io.write(opt_struc, initial_placer, format="extxyz") + + read_traj = io.read('dummy.traj', index=":") + io.write('traj.xyz', read_traj) + +def single_point(target_stru, model_path, device): + structure = io.read(target_stru) + calculator = MACECalculator(model_path=model_path, device=device) + structure.set_calculator(calculator) + energy = structure.get_potential_energy() + forces = structure.get_forces() + return energy, forces + +def main(args): + initial = io.read(args.struc_path) + initial.set_calculator(MACECalculator(model_path=args.model_path, device=args.device)) + + #optimize(initial, args.opt_struc) + energy, force = single_point(args.struc_path, args.model_path, args.device) + print(energy, forces) + + if args.vibration: + opt_struc_vib = io.read(args.opt_struc) + vib = Vibrations(opt_struc_vib) + vib.run() + vib.summary() + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description='Optimization and vibration calculation.') + parser.add_argument('struc_path', help='Path of the structure to be optimized') + parser.add_argument('model_path', help='Path of the model') + parser.add_argument('opt_struc', help='Name for the optimized structure') + parser.add_argument('-v', '--vibration', action='store_true', help='Perform vibration (default is False)') + parser.add_argument('-d', '--device', default='gpu', help='Device to use (default is gpu)') + + args = parser.parse_args() + main(args) + diff --git a/for_MACE/neb.py b/for_MACE/neb.py new file mode 100755 index 0000000..fd027eb --- /dev/null +++ b/for_MACE/neb.py @@ -0,0 +1,42 @@ +import sys +from ase import io +from ase.optimize import MDMin +from ase.dyneb import DyNEB +from ase.optimize import BFGS +from mace.calculators import MACECalculator + +# Read initial and final states: +start = sys.argv[1] +final = sys.argv[2] +model_to_path = sys.argv[3] +DEVICE = sys.argv[4] +if DEVICE == None: + DEVICE == 'gpu' +else: pass + +initial = io.read(start) +final = io.read(final) +nimages=13 + +# Make a band consisting of 5 images: +images = [initial] +images += [initial.copy() for i in range(nimages)] +images += [final] + +#neb = NEB(images,climb=True) +neb = DyNEB(images, fmax=0.05, dynamic_relaxation=True) + +# Interpolate linearly the potisions of the three middle images: +neb.interpolate() +calcs=[ MACECalculator(model_path=model_to_path, device=DEVICE) for i in range(nimages) ] + +# Set calculators: +for i in range(nimages): + images[i+1].calc = calcs[i] + +# Optimize: +optimizer = BFGS(neb, trajectory='NEB.traj') +optimizer.run(fmax=0.05) + +io.write("traj.xyz", images, write_info=False) + diff --git a/for_MACE/phonon.py b/for_MACE/phonon.py new file mode 100755 index 0000000..1f4e850 --- /dev/null +++ b/for_MACE/phonon.py @@ -0,0 +1,40 @@ +from ase.io import read, write +from ase.build import bulk +from ase.calculators.emt import EMT +from ase.phonons import Phonons + +# Setup crystal and EMT calculator +atoms = read('start-opt.xyz') + +# Phonon calculator +N = 7 +ph = Phonons(atoms, EMT(), supercell=(N, N, N), delta=0.05) +ph.run() + +# Read forces and assemble the dynamical matrix +ph.read(acoustic=True) +ph.clean() + +path = atoms.cell.bandpath('GXULGK', npoints=100) +bs = ph.get_band_structure(path) + +dos = ph.get_dos(kpts=(20, 20, 20)).sample_grid(npts=100, width=1e-3) + +# Plot the band structure and DOS: +import matplotlib.pyplot as plt # noqa +fig = plt.figure(1, figsize=(7, 4)) +ax = fig.add_axes([.12, .07, .67, .85]) + +emax = 0.035 +bs.plot(ax=ax, emin=0.0, emax=emax) + +dosax = fig.add_axes([.8, .07, .17, .85]) +dosax.fill_between(dos.get_weights(), dos.get_energies(), y2=0, color='grey', + edgecolor='k', lw=1) + +dosax.set_ylim(0, emax) +dosax.set_yticks([]) +dosax.set_xticks([]) +dosax.set_xlabel("DOS", fontsize=18) + +fig.savefig('Al_phonon.png') diff --git a/retrieve_sp_extxyz.py b/retrieve_sp_extxyz.py new file mode 100644 index 0000000..60625a0 --- /dev/null +++ b/retrieve_sp_extxyz.py @@ -0,0 +1,57 @@ +from AppOutputExtractor.FHIaims.FHIaimsOutputExtractor import extractor +import os +from itertools import groupby +import numpy as np + +eigvec_path = [os.path.join('sp', x) for x in os.listdir('sp')] +sp_path = [os.path.join(dirpath, fname) for dirpath in eigvec_path for fname in os.listdir(dirpath)] +lambda_path = [os.path.join(dirpath, fname) for dirpath in sp_path for fname in os.listdir(dirpath) if fname == 'aims.out'] +aims_out_path = sorted(lambda_path, key=lambda x: (int(x.split('/')[1]), float(x.split('/')[2].split('_')[1]))) +aims_out_path = [list(group) for key, group in groupby(aims_out_path, lambda x: int(x.split('/')[1]))] +ex = extractor() +if not os.path.exists('ext_xyz'): + os.mkdir('ext_xyz') + +for numi, i in enumerate(aims_out_path): + total_energy = [] + geometry = [] + + for j in i: + print(j) + ex.set_output_filepath(j) + ex.set_scf_blocks + + ex.get_no_atoms + + ex.get_sp_geometries(j) + ex.get_sp_atom_order() + ex.get_sp_species() + ex.get_forces + force_shape = np.shape(ex.get_forces) + get_forces = np.reshape(ex.get_forces, (force_shape[1], force_shape[2])) + + form = np.concatenate((ex.get_sp_atom_order(), ex.get_sp_geometries(j), get_forces), axis=1) + + total_energy.append(ex.get_total_energy()) + geometry.append(form) + + for numk, k in enumerate(total_energy): + with open(f"ext_xyz/ext_{j.split('/')[1]}_eigv.xyz", 'a') as f: + f.write(str(force_shape[1]) + '\n') + f.write(f'Properties-species:S:1:pos:R:3:forces:R:3 energy={total_energy[numk]} pbc="F F F"\n') + np.savetxt(f, geometry[numk], fmt="%s", delimiter=" ") + + + #ex.get_vib_eigvec +# for k in range(len(ex.set_scf_blocks)): +# #print(ex.get_total_energy(k)) +# #print('Geometry') +# print(ex.get_sp_geometries) +# #print('Atomic forces') +# #print(ex.get_forces(k)) +# #print('Eigenvector of vibrational modes') +# #print(ex.get_vib_eigvec) +# #print() +# #print() +# ex.get_total_energy + diff --git a/testing_extractor.py b/testing_extractor.py new file mode 100644 index 0000000..7ae55cc --- /dev/null +++ b/testing_extractor.py @@ -0,0 +1,14 @@ +from AppOutputExtractor.FHIaims.FHIaimsOutputExtractor import extractor + + +app_output = './aims.out' + +extractor = extractor() +extractor.set_output_filepath(app_output) + +species = extractor.get_species +no_atoms = extractor.get_no_atoms +geometries = extractor.get_geometries +order = extractor.get_atom_order +forces = extractor.get_forces +print(forces)