diff --git a/.gitignore b/.gitignore deleted file mode 100644 index f3c3970..0000000 --- a/.gitignore +++ /dev/null @@ -1,21 +0,0 @@ -# ignoring run example files -examples/* - -# ignoring built test docs -docs/build/ - -# ignoring VS Code workspace file -*.code-workspace - -# ignoring .pyc files -*.pyc - -# ignoring outputs - -# ignoring intermediate files -*.swp - -# ignore MacOS generated files -*.DS_Store - - diff --git a/GRSlib/GRS.py b/GRSlib/GRS.py index 35dab1b..d09e331 100644 --- a/GRSlib/GRS.py +++ b/GRSlib/GRS.py @@ -1,10 +1,11 @@ from GRSlib.parallel_tools import ParallelTools from GRSlib.io.input import Config from GRSlib.converters.convert_factory import convert +from GRSlib.motion.scoring import Scoring, elems, get_desc_count from GRSlib.motion.scoring_factory import scoring +from GRSlib.motion.motion import Gradient, Genetic from GRSlib.motion.scoring import Scoring from GRSlib.motion.motion import Gradient, Optimize, Create - import random, copy, os, glob, shutil import numpy as np @@ -165,6 +166,38 @@ def get_score(self,data): raise RuntimeError(">>> Found unmatched BASIS for target and current descriptors") return score + def get_ensemble(self,data): + print("get ensemble called.") + print("Using elems:",elems) + print("Ensemble Target") + scores=None + try: + scores = self.score.ensemble_score( + n_totconfig=10, # Example value + data_path='bcc.data', # Example value + cross_weight=1.0, + self_weight=1.0, + randomize_comps=False, + mincellsize=54, + maxcellsize=55, + target_comps={'W': 1.0}, # Ensure this is a dictionary + min_typ_global='box', + soft_strength=0.0, + nelements=len(elems), # Ensure elems is defined + n_descs=get_desc_count('coupling_coefficients.yace'), + mask=None, # Set this as needed + rand_comp=1) + print("Scores returned from ensemble_score:", scores) + except Exception as e: + print(f"An error occurred while calculating the ensemble score: {e}") + return None # Optionally return None or handle the error as needed + # Check if score was calculated + if scores is None or len(scores) == 0: + print("No scores were calculated returning None") + raise RuntimeError("Ensemble score could not be calculated.") + + return scores # Ensure you return the score +#""" self.score = Scoring(data, self.current_desc, self.target_desc, self.prior_desc, self.pt, self.config)score = self.score.ensemble_score() return score""" def propose_structure(self): """ @@ -187,21 +220,17 @@ def genetic_move(self,data): #3) Hybridize, Mutate based on set of rules and probabilities #4) Store socring information with best-of-generation and best-overall isolated #5) Loop until generation limit or scoring residual below threshold - if data == None: - data = self.propose_structure() - self.descriptors['current']= self.convert_to_desc(data) - try: - self.descriptors['target'] = np.load(self.config.sections['TARGET'].target_fdesc) - except: - self.descriptors['target'] = self.convert_to_desc(self.config.sections['TARGET'].target_fname) - - self.score = Scoring(self.pt, self.config, self.loss_func, data, self.descriptors) # Set scoring class to assign scores to moves - self.genmove = Optimize(self.pt, self.config, self.score) #Set desired motion class with scoring attached - - - #self.genmove.tournament_selection() - #for iterations in top_candidates[1], convert.ase_to_lammps - #self.write_output() + print("Called Genetic_Move") + if self.current_desc is None: + self.current_desc = self.propose_structure() + #if data == None: + #data = self.propose_structure() + + self.descriptors['current'] = self.convert_to_desc(data) + self.descriptors['target'] = self.convert_to_desc(self.config.sections['TARGET'].target_fname) + self.genmove = Genetic(self.pt, self.config, data, self.descriptors) + #Dont want to make a func call the default here since the user will define this? + #Need a fallback to provide a good default if a genetic move is called. # @self.pt.single_timeit def gradient_move(self,data): diff --git a/GRSlib/Ver0_Files/opt_tools.py b/GRSlib/Ver0_Files/opt_tools.py deleted file mode 100644 index 625fee5..0000000 --- a/GRSlib/Ver0_Files/opt_tools.py +++ /dev/null @@ -1,612 +0,0 @@ -import numpy as vnp -from ase.io import read,write -#from ase.atoms import Atoms -from ase import Atoms, Atom -from ase.db import connect -from icet.tools import enumerate_structures -import numpy as np -import os -from ase.db import connect -from ase.build import bulk, fcc111, add_adsorbate - -#manually tabulated minimum soft core cutoffs (if not tabulated here ASE-generated default will be used) -min_soft = {'Cr':1.65,'Fe':1.7,'Si':1.2,'V':1.65} -# TODO - make minimum cutoffs an input - -def bound_descs(descs,low_bound=-vnp.inf,up_bound=vnp.inf): - maskb = descs > up_bound - maskb2 = descs < low_bound - descs[maskb] = vnp.nan - descs[maskb2] = -vnp.nan - descs = vnp.nan_to_num(descs,nan=0.0,posinf=0.0,neginf=0.0) - return descs - -def build_target(start,save_all=False): - from ase.io import read,write - from ase import Atoms,Atom - from ase.ga.utilities import closest_distances_generator, CellBounds - from ase.ga.startgenerator import StartGenerator - from ase.data import atomic_numbers, atomic_names, atomic_masses, covalent_radii - #from __future__ import print_function - import sys,os - import ctypes - import numpy as np - from lammps import lammps, LMP_TYPE_ARRAY, LMP_STYLE_GLOBAL - - # get mpi settings from lammps - def run_struct(atoms,fname,maxcut=7.0): - - lmp = lammps() - me = lmp.extract_setting("world_rank") - nprocs = lmp.extract_setting("world_size") - - - cmds = ["-screen", "none", "-log", "none"] - lmp = lammps(cmdargs = cmds) - - #def set_atoms(atoms,atid=0): - # write('iter_%d.data' % atid,atoms,format='lammps-data') - # lmp.command('read_data iter_%d.data' % atid ) - # lmp.command('mass 1 180.94788') - # lmp.command(f"run {nsteps}") - - def run_lammps(dgradflag): - - # simulation settings - fname = file_prefix - lmp.command("clear") - lmp.command("info all out log") - lmp.command('units metal') - lmp.command('atom_style atomic') - lmp.command("boundary p p p") - lmp.command("atom_modify map hash") - lmp.command('neighbor 2.3 bin') - # boundary - lmp.command('boundary p p p') - # read atoms - lmp.command('read_data %s.data' % fname ) - utypes = [] - for atom in atoms: - if atom.symbol not in utypes: - utypes.append(atom.symbol) - for ind,typ in enumerate(utypes): - number = atomic_numbers[typ] - mass = atomic_masses[number] - lmp.command('mass %d %f' % (ind+1,mass)) - - lmp.command("pair_style zero %f" % maxcut) - lmp.command(f"pair_coeff * *") - - if dgradflag: - lmp.command(f"compute pace all pace coupling_coefficients.yace 1 1") - else: - lmp.command(f"compute pace all pace coupling_coefficients.yace 1 0") - - # run - - lmp.command(f"thermo 100") - #lmp.command(f"run {nsteps}") - lmp.command(f"run 0") - - - # declare compute pace variables - - dgradflag = 0 - run_lammps(dgradflag) - lmp_pace = lmp.numpy.extract_compute("pace", LMP_STYLE_GLOBAL, LMP_TYPE_ARRAY) - descriptor_grads = lmp_pace[ : len(atoms), : -1] - #descriptor_grads = bound_descs(lmp_pace[ : len(atoms), : -1]) - return descriptor_grads - #start = 'supercell_target.cif' - file_prefix = 'iter_%d' % 0 - if type(start) == str: - try: - atoms = read(start) - except: - raise TypeError("unrecognized file type %s" % inp) - elif type(start) == Atoms: - # except - atoms = start - else: - atoms = start - - #atoms = read(start) - - write('%s.data' % file_prefix,atoms,format='lammps-data') - start_arr = run_struct(atoms, '%s.data'% file_prefix) - if save_all: - vnp.save('A_target_all.npy',start_arr) - avg_start = vnp.average(start_arr,axis=0) - var_start = vnp.var(start_arr,axis=0) - vnp.save('target_descriptors.npy',avg_start) - vnp.save('target_var_descriptors.npy',var_start) - return avg_start, var_start - -def rand_comp(model_dict): - old_vals = vnp.array(list(model_dict.values())) - vals = vnp.random.rand(len(old_vals)) - newvals = vals/vals.sum() - new_dict = {list(model_dict.keys())[ii]:newvals[ii] for ii in range(len(newvals))} - return new_dict - -def generate_random_integers(sum_value, n): - # Generate a list of n - 1 random integers - random_integers = [vnp.random.randint(0, sum_value) for _ in range(n - 1)] - # Ensure the sum of the random integers is less than or equal to sum_value - random_integers.sort() - # Calculate the Nth integer to ensure the sum is equal to sum_value - random_integers.append(sum_value - sum(random_integers)) - return random_integers - - -class System_Enum: - def __init__(self,base_species): - self.a = None - self.blocks = [] - self.crystal_structures = [] - self.lattice_bases = [] - self.stored_lattice = {'a':{}} - self.base_species = base_species - return None - - - def set_crystal_structures(self,structs = ['bcc']): - self.crystal_structures = structs - - def set_lattice_constant(self,a): - if type(a) == float: - self.a = [a]*len(self.crystal_structures) - elif type(a) == list: - self.a = a - else: - raise TypeError("cannot use type other than float or list for lattice constant(s)") - - return None - - def set_substitutional_blocks(self,blocks=None): - assert len(self.crystal_structures) >= 1, "set crystal structures before defining substitutional blocks - default is 'fcc' " - if blocks == None: - blocks = [ [self.base_species] ] * len(self.crystal_structures) - else: - blocks = blocks - self.blocks = blocks - - def enumerate_structures(self, min_int_mult, max_int_mult): - all_structs = [] - assert self.a != None, "set_lattice_constant first, then do enumeration" - assert self.blocks != None, "set_substitutuional_blocks first, then do structure enumeration" - for icrystal,crystal in enumerate(self.crystal_structures): - primitive = bulk(self.base_species[0] , crystal , a=self.a[icrystal] , cubic=False) - print ('generating structures for crystal: %s' % crystal) - enumerated = enumerate_structures(primitive, range(min_int_mult, max_int_mult), self.blocks[icrystal]) - sublist = [ i for i in enumerated ] - self.stored_lattice['a'][crystal] = self.a[icrystal] - all_structs.append(sublist) - self.all_structs = all_structs - return all_structs - -def starting_generation(pop_size,all_species,cell,typ='ase',nchem = 1,use_template=None): - pop = [] - if typ == 'ase': - from ase.io import read,write - from ase import Atoms,Atom - from ase.ga.utilities import closest_distances_generator, CellBounds - from ase.ga.startgenerator import StartGenerator - from ase.data import atomic_numbers - volume = vnp.dot(cell[2],vnp.cross(cell[0],cell[1])) - # Target cell volume for the initial structures, in angstrom^3 - #volume = 240. # Specify the 'building blocks' from which the initial structures - # will be constructed. Here we take single Ag atoms as building - # blocks, 24 in total. - #blocks = [('Ag', 24)] - ## We may also write: - #blocks = ['Ag'] * 24 - dsize = int(pop_size/nchem) - sort_specs = sorted(all_species) - nats = len(all_species) - uspecs = list(set(sort_specs)) - if len(uspecs) == 1: - block_sets = [ [(uspecs[0],len(sort_specs))] ]* nchem - else: - block_sets = [] - for iii in range(nchem): - natoms_per_type = generate_random_integers(len(sort_specs), len(uspecs)) - block_sets.append( [ tuple([uspecs[ij],natoms_per_type[ij]]) for ij in range(len(uspecs)) ] ) - #blocks = [('Ti', 4), ('O', 8)] - # Generate a dictionary with the closest allowed interatomic distances - Zs = [ atomic_numbers[sym] for sym in list(set(all_species))] - blmin = closest_distances_generator(atom_numbers=Zs, ratio_of_covalent_radii=0.5) - natoms = len(all_species) - slab = Atoms('',cell=cell, pbc=True) - for block in block_sets: - # Initialize the random structure generator - #sg = StartGenerator(slab, all_species, blmin, box_volume=volume, - sg = StartGenerator(slab, block, blmin, number_of_variable_cell_vectors=0) - # and add them to the database - for i in range(dsize): - a = sg.get_new_candidate() - pop.append(a) - if use_template: - pop = [use_template] + pop - return pop - -def at_to_lmp(atoms,index,temperature=10000.0,min_typ='temp', coefftypes=True,soft_strength=10000.0): - from ase.data import atomic_numbers, atomic_names, atomic_masses, covalent_radii - from ase.ga.utilities import closest_distances_generator - #s=generate.format(xx,yy,zz,dx,dz,dz,n_atoms,seed,index,random.uniform(0.0,1.0)) - fname = 'ats_%s.data' % index - write(fname,atoms,format='lammps-data') - generate=\ -""" -units metal -boundary p p p - -read_data {} - -log log_{}.lammps -""" - s1 = generate.format(fname,index) - atutypes =[] - for atom in atoms: - if atom.symbol not in atutypes: - atutypes.append(atom.symbol) - if coefftypes: - utypes = get_desc_count(coefffile='coupling_coefficients.yace',return_elems=True) - else: - utypes = atutypes.copy() - if len(atutypes) < len(utypes): - utypes = atutypes.copy() - typstr = ' ' - try: - Zs = [ atomic_numbers[sym] for sym in utypes] - blmin = closest_distances_generator(atom_numbers=Zs, ratio_of_covalent_radii=0.8) - min_soft_utypes = list(blmin.values()) - except KeyError: - print('missing tabulated minimum bonds, applying minimums from ASE') - for indtyp,utype in enumerate(utypes): - atnum = atomic_numbers[utype] - atmass = atomic_masses[atnum] - s1 += 'mass %d %f\n' % (indtyp+1, atmass) - typstr += ' %s' % utype - #print('typstr final: %s' % typstr) - #print (utypes,atutypes) - if min_typ =='temp': - generate2a=\ -"""pair_style hybrid/overlay soft %2.3f mliap model mliappy LATER descriptor ace coupling_coefficients.yace -pair_coeff * * soft %f\n""" % (max(min_soft_utypes),soft_strength) - generate2b=\ -"""pair_coeff * * mliap %s - -thermo 10 -fix nve all nve -fix lan all langevin %d 100 1.0 48279 - -velocity all create %d 4928459 dist gaussian -""" % (typstr, int(temperature), int(temperature*2)) - generate2 = generate2a+generate2b - elif min_typ == 'box': - generate2a =\ -"""pair_style hybrid/overlay soft %2.3f mliap model mliappy LATER descriptor ace coupling_coefficients.yace -pair_coeff * * soft %f\n""" % (max(min_soft_utypes),soft_strength) - generate2b =\ -"""pair_coeff * * mliap %s - -thermo 10 -velocity all create 1. 4928459 dist gaussian -fix 1b all box/relax iso 0.0 vmax 0.001""" % typstr - generate2 = generate2a+generate2b - else: - generate2a =\ -"""pair_style hybrid/overlay soft %2.3f mliap model mliappy LATER descriptor ace coupling_coefficients.yace -pair_coeff * * soft %f\n""" % (max(min_soft_utypes),soft_strength) - generate2b = \ -"""pair_coeff * * mliap %s - -thermo 10 -velocity all create 1. 4928459 dist gaussian""" % typstr - generate2 = generate2a+generate2b -#minimize 1e-8 1e-8 1000 1000" - #s = generate.format(fname,index,vnp.random.uniform(0.0,1.0)) - s = s1 + generate2 - return s - - -from crystal_enum import * -def prim_crystal(elem_list): - all_prims = get_prim_structs(elem_list, multiatom=False) - myind = np.random.choice(range(len(all_prims))) - return all_prims[myind] - - -from hnf import * -#elems[0],desired_size,volfrac=1.0,cubic=True,override_lat='fcc',override_a=2.98 -def bulk_template(elem,desired_size,volfrac=1.0,cubic=True,override_lat=None,override_a=None): - if not override_a: - prim = bulk(elem,cubic=False) - if cubic: - atoms = bulk(elem,cubic=True) - else: - atoms = prim - - if type(desired_size) == tuple: - atoms = atoms*desired_size - else: - atoms = prim * (desired_size,desired_size,desired_size) - return atoms - else: - if override_lat != None: - prim = bulk(elem,override_lat,a=override_a,cubic=False) - elif override_lat == None: - prim = bulk(elem,a=override_a,cubic=False) - if cubic: - if override_lat != None: - atoms = bulk(elem,override_lat,override_a,cubic=True) - elif override_lat == None: - atoms = bulk(elem,override_a,cubic=True) - else: - atoms = prim - - if type(desired_size) == tuple: - atoms = atoms*desired_size - else: - atoms = prim * (desired_size,desired_size,desired_size) - cell_this = atoms.get_cell() - atoms.set_cell(cell_this*volfrac,scale_atoms=True) - return atoms - -def internal_basic(atoms,index=0,min_typ=None,soft_strength=1.0): - s = at_to_lmp(atoms,index,min_typ=min_typ,soft_strength=soft_strength) - return s - -def bulk_sis_template(base_species,cellmaxmult,crystal_types=['bcc'],lattice_constants=[2.88],cellmin=1): - base_species_pairs = [base_species, - ] - collected_atoms =[] - assert len(crystal_types) == 1, "must have only one crystal type at a time" - chem_lst = ['%s']*len(base_species) - chem_str = '-'.join(b for b in chem_lst) % tuple(base_species) - this_db_prefix = chem_str + '_' + crystal_types[0] + '_' + str(cellmaxmult) - if os.path.isfile('%s.db' % this_db_prefix): - print('has db') - db = connect('%s.db' % this_db_prefix) - for row in db.select(): - atoms = row.toatoms() - collected_atoms.append(atoms) - else: - print('generating db') - for base_species in base_species_pairs: - db = connect('%s.db' % this_db_prefix) - se = System_Enum(base_species) - crystal_structures = crystal_types - se.set_lattice_constant(lattice_constants) - se.set_crystal_structures(crystal_structures) - se.set_substitutional_blocks([base_species,base_species]) - astrcts = se.enumerate_structures(cellmin,cellmaxmult+1) - chem_lst = ['%s']*len(base_species) - chem_str = '-'.join(b for b in chem_lst) % tuple(base_species) - - for icrystal,crystal in enumerate(crystal_structures): - this_db_prefix = chem_str + '_' + crystal + '_' + str(cellmaxmult) - strcts = astrcts[icrystal] - for istrct , strct in enumerate(strcts): - if len(strct) > 1: - db.write(strct) - collected_atoms.append(strct) - else: - db.write(strct*(2,1,1)) - collected_atoms.append(strct*(2,1,1)) - return collected_atoms -#base_species = ['Cr','Fe'] -#cellmaxmult=4 -#bulk_sis_template(base_species,cellmaxmult,crystal_types=['bcc'],lattice_constants=[2.88]) - -def internal_generate_cell(index,desired_size=4,template=None,desired_comps={'Ni':1.0},use_template=None,min_typ='temp',soft_strength=10000,sis_freeze=False): - #from ase.build import bulk - #from ase import Atoms,Atom - if template == None: - #from ase.build import bulk - #from ase import Atoms,Atom - chems = list(desired_comps.keys()) - template = Atoms([chems[0]]*desired_size) - atoms_base = bulk(chems[0]) - vol_base = vnp.dot(vnp.cross(atoms_base.get_cell()[0],atoms_base.get_cell()[1]),atoms_base.get_cell()[2]) - a_simp = vol_base**(1/3) - #cells = get_hnfs(hnf_trs=[desired_size]) - cells_all = get_hnfs(hnf_trs=[desired_size]) - toli = int((desired_size)**(1/3)) - cells= limit_mats_len(cells_all,desired_size,tol=0.16) - #print ('ncells',len(cells)) - try: - cell = a_simp * cells[vnp.random.choice(range(len(cells)))] - except ValueError: - cell = a_simp * cells_all [vnp.random.choice(range(len(cells_all)))] - #print ('hnf',cell) - norms = np.array([np.linalg.norm(k) for k in cell]) - norms /= (desired_size) - #print('hnf norms',norms) - #cell = a_simp * cells[-2] - template.set_cell(cell) - template.set_scaled_positions(vnp.random.uniform(0,1,(desired_size,3))) - else: - tempalte = template - - new_comps = {elem:int(round(len(template)*cmp))/len(template) for elem,cmp in desired_comps.items()} - print ('for structure of size:%d'% len(template),'desired compositions:', desired_comps,'will be replaced with', new_comps) - all_species = get_target_comp_s(desired_size=len(template),desired_comps=new_comps,parent_cell=template.get_cell()) - #all_species = get_target_comp(desired_size=len(template),desired_comps=new_comps,parent_cell=template.get_cell()) - if len(all_species) >= len(template): - print('doing cheap fix for cell size') - all_species=all_species[:len(template)] - elif len(all_species) <= len(template): - print('doing cheap grow for cell size') - diff = len(template)-len(all_species) - all_species = all_species + all_species[:diff] - #print ('all specs vs template',len(all_species),len(template)) - #assert len(all_species)== len(template), "composition list size must match size of template atoms" - cellg = template.get_cell() - scpos = template.get_scaled_positions() - - if use_template: - if not sis_freeze: - #rnd = starting_generation(1,all_species,cellg,typ='ase',use_template=template)[0] - np.random.shuffle(all_species) - rnd = Atoms(all_species) - rnd.set_cell(cellg) - rnd.set_pbc(True) - rnd.set_scaled_positions(scpos) - else: - rnd = template - else: - rnd = starting_generation(1,all_species,cellg,typ='ase')[0] - s = at_to_lmp(rnd,index,min_typ=min_typ,soft_strength=soft_strength) - return s - -def get_comp(atoms,symbols): - comps = {symbol: 0.0 for symbol in symbols} - counts = {symbol: 0 for symbol in symbols} - atsymbols = [atom.symbol for atom in atoms] - for atsymbol in atsymbols: - counts[atsymbol] +=1 - for symbol in symbols: - comps[symbol] = counts[symbol]/len(atoms) - return comps - -def flip_one_atom(atoms,types): - new_atoms = atoms.copy() - flip_ind = vnp.random.randint(0,len(atoms)) - flip_current = new_atoms[flip_ind].symbol - excluded = [typ for typ in types if typ != flip_current] - flip_to_ind = vnp.random.randint(0,len(excluded)) - flip_to_type = excluded[flip_to_ind] - new_atoms[flip_ind].symbol = flip_to_type - return new_atoms - -def flip_N_atoms(atoms,types,fraction=0.25): - pert_inds = vnp.random.choice(range(len(atoms)),size=int(len(atoms)*fraction) ) - new_atoms = atoms.copy() - for pert_ind in pert_inds: - flip_ind = vnp.random.randint(0,len(atoms)) - flip_current = new_atoms[flip_ind].symbol - excluded = [typ for typ in types if typ != flip_current] - flip_to_ind = vnp.random.randint(0,len(excluded)) - flip_to_type = excluded[flip_to_ind] - new_atoms[flip_ind].symbol = flip_to_type - return new_atoms - - -def add_atom(atoms,symbols,tol = 0.5): - from ase import Atom,Atoms - from ase.ga.utilities import closest_distances_generator - from ase.data import atomic_numbers - from ase.neighborlist import primitive_neighbor_list - blmin = closest_distances_generator(atom_numbers=[atomic_numbers[symbol] for symbol in symbols] + [atomic_numbers['Ne']], ratio_of_covalent_radii=0.5) - def readd(): - symbol = vnp.random.choice(symbols) - rnd_pos_scale = vnp.random.rand(1,3) - rnd_pos = vnp.matmul(atoms.get_cell(),rnd_pos_scale.T) - rnd_pos = rnd_pos.T[0] - new_atom = Atom('Ne',rnd_pos) - tst_atoms = atoms.copy() - tst_atoms.append(new_atom) - tst_atoms.wrap() - rc = 5. - - atinds = [atom.index for atom in tst_atoms] - at_dists = {i:[] for i in atinds} - all_dists = [] - nl = primitive_neighbor_list('ijdD',pbc=tst_atoms.pbc,positions=tst_atoms.positions ,cell=atoms.get_cell(),cutoff=rc) - bond_types = {i:[] for i in atinds} - for i,j in zip(nl[0],nl[-1]): - at_dists[i].append(j) - for i,j in zip(nl[0],nl[1]): - bond_types[i].append( (atomic_numbers[tst_atoms[i].symbol] , atomic_numbers[tst_atoms[j].symbol]) ) - return symbol, tst_atoms, at_dists, rnd_pos, bond_types - symbol, tst_atoms , at_dists , rnd_pos, bond_types = readd() - bondtyplst = list(bond_types.keys()) - syms = [tst_atom.symbol for tst_atom in tst_atoms] - tst_id = syms.index('Ne') - tst_dists = at_dists[tst_id] - tst_bonds = bond_types[tst_id] - conds = all([ vnp.linalg.norm(tst_dist) >= blmin[(atomic_numbers[symbol] , tst_bonds[i][1])] for i,tst_dist in enumerate(tst_dists)]) - while not conds: - symbol , tst_atoms, at_dists , rnd_pos, bond_types = readd() - syms = [tst_atom.symbol for tst_atom in tst_atoms] - tst_id = syms.index('Ne') - tst_dists = at_dists[tst_id] - tst_bonds = bond_types[tst_id] - #conds = all([ vnp.linalg.norm(tst_dist) >= tol for tst_dist in tst_dists]) - #conds = all([ vnp.linalg.norm(tst_dist) >= blmin[tst_bonds[i]] for i,tst_dist in enumerate(tst_dists)]) - conds = all([ vnp.linalg.norm(tst_dist) >= blmin[(atomic_numbers[symbol] , tst_bonds[i][1])]-tol for i,tst_dist in enumerate(tst_dists)]) - atoms.append(Atom(symbol,rnd_pos)) - return atoms - -def get_target_comp_s(desired_size,desired_comps,parent_cell): - num_ats = {key:int(desired_comps[key]*desired_size) for key in list(desired_comps.keys())} - symbols =[] - for key,nrepeat in num_ats.items(): - symbols.extend([key]*nrepeat) - return symbols - -def get_target_comp(desired_size,desired_comps,parent_cell): - from ase import Atoms - initial_a = Atoms('',pbc=True) - initial_a.set_cell(parent_cell) - - current_comps = {key:0.0 for key in list(desired_comps.keys())} - current_size = 0 - - chems = list(desired_comps.keys()) - comp_conds = all([current_comps[chem] == desired_comps[chem] for chem in chems]) - these_atoms = initial_a.copy() - max_iter = 100 - this_iter = 0 - toli = 1/(desired_size) - while current_size <= desired_size or not comp_conds and this_iter <= max_iter: - #t = 1000 * time.time() # current time in milliseconds - #vnp.random.seed(int(t) % 2**32) - if current_size < desired_size: - tst_ats = add_atom(these_atoms,chems) - tst_comps = get_comp(tst_ats, chems) - comp_conds = all([tst_comps[chem] == desired_comps[chem] for chem in chems]) - these_atoms = tst_ats.copy() - current_size = len(these_atoms) - Qi = [vnp.abs(tst_comps[chem] - desired_comps[chem]) for chem in chems] - Qi = round(vnp.sum(Qi),8) - elif current_size == desired_size: - #tst_ats = flip_one_atom(these_atoms,chems) - tst_ats = flip_N_atoms(these_atoms,chems,fraction=vnp.random.rand()) - tst_comps = get_comp(tst_ats, chems) - comp_conds = all([tst_comps[chem] == desired_comps[chem] for chem in chems]) - tstQi = [vnp.abs(tst_comps[chem] - desired_comps[chem]) for chem in chems] - tstQi = round(vnp.sum(tstQi),8) - #print (len(these_atoms),tst_comps,tstQi) - if tstQi <= Qi+toli: - these_atoms = tst_ats.copy() - Qi = tstQi - #print ('in composition loop',tst_comps,comp_conds) - if current_size==desired_size and comp_conds: - break - this_iter += 1 - return list(tst_ats.symbols) - -def get_desc_count(coefffile,return_elems=False): - with open(coefffile,'r') as readcoeff: - lines = readcoeff.readlines() - elemline = [line for line in lines if 'elements' in line][0] - elemstr = elemline.split(':')[-1] - elemstr2 = elemstr.replace('[','') - elemstr3 = elemstr2.replace(']','') - elemstr4 = elemstr3.replace(',','') - elems = elemstr4.split() - nelements = len(elems) - desclines = [line for line in lines if 'mu0' in line] - if not return_elems: - return int(len(desclines)/nelements) - #return int(len(desclines)) - else: - return elems - -#rnd = internal_generate_cell(0,desired_comps={'W':0.9,'H':0.1}) -#rnd = internal_generate_cell(index,desired_comps={'W':0.9,'H':0.1}) -#rnd = internal_generate_cell(1,desired_comps={'Cr':1.0}) -#print(rnd) diff --git a/GRSlib/converters/sections/lammps_ace.py b/GRSlib/converters/sections/lammps_ace.py deleted file mode 100644 index 2667aa0..0000000 --- a/GRSlib/converters/sections/lammps_ace.py +++ /dev/null @@ -1,68 +0,0 @@ -from GRSlib.converters.convert import Convert -from GRSlib.converters.sections.lammps_base import Base, _extract_compute_np -import numpy as np - -class Ace(Convert): - - def __init__(self, name, pt, config): - super().__init__(name, pt, config) - self._data = {} - self._i = 0 - self._lmp = None - self._row_index = 0 - self.pt.check_lammps() - - def _prepare_lammps(self): - #Base._set_structure(self) - # this is super clean when there is only one value per key, needs reworking - # self._set_variables(**_lammps_variables(config.sections["ACE"].__dict__)) - - self._lmp.command(f"variable rcutfac equal {max(self.config.sections['BASIS'].rcutfac)}") - self._lmp.command(f"pair_style zero {max(self.config.sections['BASIS'].rcutfac)}") - self._lmp.command("pair_coeff * *") - - compute_name='pace' - numtypes = len(self.config.sections['BASIS'].elements) - base_pace = "compute %s all pace coupling_coefficients.yace %s %s" % (compute_name,self.config.sections['BASIS'].bikflag,self.config.sections['BASIS'].dgradflag) - self._lmp.command(base_pace) - return compute_name - - def run_lammps_single(self,data): - Base._initialize_lammps(self) - Base._set_structure(self) - self._lmp.command("read_data %s" % data) - self._prepare_lammps() - Base._set_neighbor_list(self) - Base._run_lammps(self) - descriptor_vals = self._collect_lammps_single() - self._lmp = self.pt.close_lammps() - return descriptor_vals - - def _collect_lammps_single(self): - num_atoms = self._lmp.extract_global("natoms") - num_types = self.config.sections['BASIS'].numtypes - with open('coupling_coefficients.yace','r') as readcoeff: - lines = readcoeff.readlines() - elemline = [line for line in lines if 'elements' in line][0] - elemstr = elemline.split(':')[-1] - elemstr2 = elemstr.replace('[','') - elemstr3 = elemstr2.replace(']','') - elemstr4 = elemstr3.replace(',','') - elems = elemstr4.split() - nelements = len(elems) - desclines = [line for line in lines if 'mu0' in line] - - #ncols_pace = int(len(desclines)/nelements) - ncols_pace = int(len(desclines)/nelements) + nelements - nrows_pace = num_atoms - lmp_pace = _extract_compute_np(self._lmp, "pace", 0, 2, (nrows_pace, ncols_pace)) - - if (np.isinf(lmp_pace)).any() or (np.isnan(lmp_pace)).any(): - self.pt.single_print('WARNING! Applying np.nan_to_num()') - lmp_pace = np.nan_to_num(lmp_pace) - if (np.isinf(lmp_pace)).any() or (np.isnan(lmp_pace)).any(): - raise ValueError('Nan in computed data of file') -# print("Got descriptors, returning") -# print(np.shape(lmp_pace)) - return lmp_pace - diff --git a/GRSlib/io/sections/basis.py b/GRSlib/io/sections/basis.py deleted file mode 100644 index 1a13eec..0000000 --- a/GRSlib/io/sections/basis.py +++ /dev/null @@ -1,195 +0,0 @@ -import numpy as np -import itertools -from GRSlib.io.sections.sections import Section - -try: - - from fitsnap3lib.lib.sym_ACE.pa_gen import * - from fitsnap3lib.lib.sym_ACE.yamlpace_tools.potential import * - from fitsnap3lib.lib.sym_ACE.wigner_couple import * - from fitsnap3lib.lib.sym_ACE.clebsch_couple import * - - class Basis(Section): - - def __init__(self, name, config, pt, infile, args): - super().__init__(name, config, pt, infile, args) - - allowedkeys = ['descriptor', 'numTypes', 'ranks', 'lmax', 'nmax', 'mumax', 'nmaxbase', 'rcutfac', 'lambda', - 'elements', 'bzeroflag', 'erefs', 'rcinner', 'drcinner', 'RPI_heuristic', 'lmin', - 'bikflag', 'dgradflag','wigner_flag','b_basis','manuallabs'] - for value_name in config['BASIS']: - if value_name in allowedkeys: continue - else: - raise RuntimeError(">>> Found unmatched variable in BASIS section of input: ",value_name) - self.descriptor = self.get_value("BASIS", "descriptor", "ACE") - self.numtypes = len(self.get_value("BASIS", "elements", "H")) - self.ranks = self.get_value("BASIS","ranks","3").split() - self.lmin = self.get_value("BASIS", "lmin", "0").split() - self.lmax = self.get_value("BASIS", "lmax", "2").split() - self.nmax = self.get_value("BASIS", "nmax", "2").split() - #self.mumax = self.get_value("BASIS","mumax", "1") - self.nmaxbase = self.get_value("BASIS", "nmaxbase", "16","int") - self.rcutfac = self.get_value("BASIS", "rcutfac", "4.5").split() - self.lmbda = self.get_value("BASIS","lambda",'1.35').split() - self.rcinner = self.get_value("BASIS","rcinner",'0.0').split() - self.drcinner = self.get_value("BASIS","drcinner",'0.01').split() - self.elements = self.get_value("BASIS", "elements", "H").split() - self.mumax = len(self.elements) - #self.erefs = self.get_value("ACE", "erefs", "0.0").split() - self.erefs = [0.0] * len(self.elements) - self.bikflag = self.get_value("BASIS", "bikflag", "1", "bool") - self.dgradflag = self.get_value("BASIS", "dgradflag", "0", "bool") - self.b_basis = self.get_value("BASIS" , "b_basis" , "pa_tabulated") - self.manuallabs = self.get_value("BASIS", "manuallabs", 'None') - self.type_mapping = {} - for i, atom_type in enumerate(self.elements): - self.type_mapping[atom_type] = i+1 - - self.bzeroflag = self.get_value("BASIS", "bzeroflag", "1", "bool") - self.wigner_flag = self.get_value("BASIS", "wigner_flag", "1", "bool") - - #if self.bikflag: - # self._assert_dependency('bikflag', "CALCULATOR", "per_atom_energy", True) - self.lmax_dct = {int(rnk):int(lmx) for rnk,lmx in zip(self.ranks,self.lmax)} - if self.b_basis != 'pa_tabulated': - self.pt.single_print('WARNING: Only change ACE basis flags if you know what you are doing!') - self._generate_b_list() - self._write_couple() - Section.num_desc = len(self.blist) - #print("Generating ACE lables yields this many cofficients: ",Section.num_desc) - self.delete() - - def _generate_b_list(self): - self.blist = [] - self.nus = [] - self.blank2J = [] - prefac = 1.0 - i = 0 - - if self.manuallabs != 'None': - with open(self.manuallabs,'r') as readjson: - labdata = json.load(readjson) - ranked_chem_nus = [list(ik) for ik in list(labdata.values())] - elif self.manuallabs == 'None' and self.b_basis == 'minsub': - from fitsnap3lib.lib.sym_ACE.rpi_lib import descriptor_labels_YSG - if type(self.lmin) == list: - if len(self.lmin) == 1: - self.lmin = self.lmin * len(self.ranks) - ranked_chem_nus = [descriptor_labels_YSG(int(rnk), int(self.nmax[ind]), int(self.lmax[ind]), int(self.mumax),lmin = int(self.lmin[ind]) ) for ind,rnk in enumerate(self.ranks)] - else: - ranked_chem_nus = [descriptor_labels_YSG(int(rnk), int(self.nmax[ind]), int(self.lmax[ind]), int(self.mumax),lmin = int(self.lmin) ) for ind,rnk in enumerate(self.ranks)] - elif self.manuallabs == 'None' and self.b_basis == 'pa_tabulated': - ranked_chem_nus = [] - if len(self.lmin) == 1: - self.lmin = self.lmin * len(self.ranks) - for ind,rank in enumerate(self.ranks): - rank = int(rank) - PA_lammps, not_compat = pa_labels_raw(rank,int(self.nmax[ind]),int(self.lmax[ind]), int(self.mumax),lmin = int(self.lmin[ind]) ) - ranked_chem_nus.append(PA_lammps) - if len(not_compat) > 0: - self.pt.single_print('Functions incompatible with lammps for rank %d : '% rank, not_compat) - highranks = [int(r) for r in self.ranks if int(r) >= 5] - warnflag = any([ self.lmax_dct[rank] >= 5 and self.lmin[ind] > 1 for ind,rank in enumerate(highranks)]) - if warnflag: - self.pt.single_print('WARNING: lmax and lmin for your current max rank will generate descriptors that cannot be entered into LAMMPS_PACE - try a lower lmax for ranks >= 4' % warnflag[0]) - nus_unsort = [item for sublist in ranked_chem_nus for item in sublist] - nus = nus_unsort.copy() - mu0s = [] - mus =[] - ns = [] - ls = [] - for nu in nus_unsort: - mu0ii,muii,nii,lii = get_mu_n_l(nu) - mu0s.append(mu0ii) - mus.append(tuple(muii)) - ns.append(tuple(nii)) - ls.append(tuple(lii)) - nus.sort(key = lambda x : mus[nus_unsort.index(x)],reverse = False) - nus.sort(key = lambda x : ns[nus_unsort.index(x)],reverse = False) - nus.sort(key = lambda x : ls[nus_unsort.index(x)],reverse = False) - nus.sort(key = lambda x : mu0s[nus_unsort.index(x)],reverse = False) - nus.sort(key = lambda x : len(x),reverse = False) - nus.sort(key = lambda x : mu0s[nus_unsort.index(x)],reverse = False) - byattyp = srt_by_attyp(nus) - #config.nus = [item for sublist in list(byattyp.values()) for item in sublist] - for atype in range(self.numtypes): - nus = byattyp[str(atype)] - for nu in nus: - i += 1 - mu0,mu,n,l,L = get_mu_n_l(nu,return_L=True) - if L != None: - flat_nu = [mu0] + mu + n + l + list(L) - else: - flat_nu = [mu0] + mu + n + l - self.blist.append([i] + flat_nu) - self.nus.append(nu) - self.blank2J.append([prefac]) - self.ncoeff = int(len(self.blist)/self.numtypes) - if not self.bzeroflag: - self.blank2J = np.reshape(self.blank2J, (self.numtypes, int(len(self.blist)/self.numtypes))) - onehot_atoms = np.ones((self.numtypes, 1)) - self.blank2J = np.concatenate((onehot_atoms, self.blank2J), axis=1) - self.blank2J = np.reshape(self.blank2J, (len(self.blist) + self.numtypes)) - else: - self.blank2J = np.reshape(self.blank2J, len(self.blist)) - - def _write_couple(self): - @self.pt.sub_rank_zero - def decorated_write_couple(): - reference_ens = [0.0] * len(self.elements) - bondinds=range(len(self.elements)) - bonds = [b for b in itertools.product(bondinds,bondinds)] - bondstrs = ['[%d, %d]' % b for b in bonds] - assert len(self.lmbda) == len(bondstrs), "must provide rc, lambda, for each BOND type" - assert len(self.rcutfac) == len(bondstrs), "must provide rc, lambda, for each BOND type" - if len(self.lmbda) == 1: - lmbdavals = self.lmbda - rcvals = self.rcutfac - rcinnervals = self.rcinner - drcinnervals = self.drcinner - if len(self.lmbda) > 1: - lmbdavals = {bondstr:lmb for bondstr,lmb in zip(bondstrs,self.lmbda)} - rcvals = {bondstr:lmb for bondstr,lmb in zip(bondstrs,self.rcutfac)} - rcinnervals = {bondstr:lmb for bondstr,lmb in zip(bondstrs,self.rcinner)} - drcinnervals = {bondstr:lmb for bondstr,lmb in zip(bondstrs,self.drcinner)} - - ldict = {int(rank):int(lmax) for rank,lmax in zip(self.ranks,self.lmax)} - L_R = 0 - M_R = 0 - rankstrlst = ['%s']*len(self.ranks) - rankstr = ''.join(rankstrlst) % tuple(self.ranks) - lstrlst = ['%s']*len(self.ranks) - lstr = ''.join(lstrlst) % tuple(self.lmax) - if not self.wigner_flag: - try: - with open('cg_LR_%d_r%s_lmax%s.pickle' %(L_R,rankstr,lstr),'rb') as handle: - ccs = pickle.load(handle) - except FileNotFoundError: - ccs = get_cg_coupling(ldict,L_R=L_R) - #print (ccs) - #store them for later so they don't need to be recalculated - store_generalized(ccs, coupling_type='cg',L_R=L_R) - else: - try: - with open('wig_LR_%d_r%s_lmax%s.pickle' %(L_R,rankstr,lstr),'rb') as handle: - ccs = pickle.load(handle) - except FileNotFoundError: - ccs = get_wig_coupling(ldict,L_R) - #print (ccs) - #store them for later so they don't need to be recalculated - store_generalized(ccs, coupling_type='wig',L_R=L_R) - - apot = AcePot(self.elements, reference_ens, [int(k) for k in self.ranks], [int(k) for k in self.nmax], [int(k) for k in self.lmax], self.nmaxbase, rcvals, lmbdavals, rcinnervals, drcinnervals, [int(k) for k in self.lmin], self.b_basis, **{'ccs':ccs[M_R]}) - apot.write_pot('coupling_coefficients') - - decorated_write_couple() - -except ModuleNotFoundError: - - class Ace(Section): - """ - Dummy class for factory to read if torch is not available for import. - """ - def __init__(self, name, config, pt, infile, args): - super().__init__(name, config, pt, infile, args) - raise ModuleNotFoundError("Missing sympy or pyyaml modules.") diff --git a/GRSlib/io/sections/scoring.py b/GRSlib/io/sections/scoring.py deleted file mode 100644 index 777d36c..0000000 --- a/GRSlib/io/sections/scoring.py +++ /dev/null @@ -1,31 +0,0 @@ -from GRSlib.io.sections.sections import Section - -#pt = ParallelTools() -#output = Output() -class Scoring(Section): - - def __init__(self, name, config, pt,infile, args): - super().__init__(name, config, pt, infile,args) - self.allowedkeys = ['score_type','moments', 'moments_coeff', 'moments_bonus', 'strength_target', - 'strength_prior','exact_distribution','internal_entropy','ensemble_entropy'] - self._check_section() - self.score_type = self.get_value("SCORING", "score_type", None) - self.strength_target = self.get_value("SCORING", "strength_target", 1.0) - self.strength_prior = self.get_value("SCORING", "strength_prior", 0.0) - self.exact_distribution = self.get_value("SCORING", "exact_distribution", False) - if self.score_type == "moments": - self.moments = self.get_value("SCORING", "moments", "mean stdev").split() - #options are : mean stdev skew kurtosis - self.moments_coeff = self.get_value("SCORING", "moments_coeff", "1.0 0.1").split() - #Requires one number per moment - self.moments_bonus = self.get_value("SCORING", "moments_bonus", "0.0 0.0").split() - #Score reduction if exact moment value is matched - elif self.score_type == "entropy": - self.internal_entropy = self.get_value("SCORING", "internal_entropy", 1.0) - # Strength of interactions within a single structure. - # Positive value prefers uniform atom environments, negative value promotes diversity (higher entropy) - self.ensemble_entropy = self.get_value("SCORING", "ensemble_entropy", 1.0) - # Strength of interactions with prior structures. - # Positive value attracts current structure to priors, negative value promotes diversity (higher entropy) - self.delete() - diff --git a/GRSlib/motion/genetic_moves/density.py b/GRSlib/motion/genetic_moves/density.py index f2a23e3..4d1bfee 100644 --- a/GRSlib/motion/genetic_moves/density.py +++ b/GRSlib/motion/genetic_moves/density.py @@ -16,9 +16,10 @@ def readd(): rnd_pos_scale = vnp.random.rand(1,3) rnd_pos = vnp.matmul(atoms.get_cell(),rnd_pos_scale.T) rnd_pos = rnd_pos.T[0] + #new atom being created and added to a list of already existing atoms new_atom = Atom('Ne',rnd_pos) tst_atoms = atoms.copy() - tst_atoms.append(new_atom) + tst_atoms.append(new_atom) #add a new atom to copy tst_atoms.wrap() rc = 5. @@ -32,6 +33,7 @@ def readd(): for i,j in zip(nl[0],nl[1]): bond_types[i].append( (atomic_numbers[tst_atoms[i].symbol] , atomic_numbers[tst_atoms[j].symbol]) ) return symbol, tst_atoms, at_dists, rnd_pos, bond_types + #positioning check for the atoms to make sure the periodic boundary conditions are met symbol, tst_atoms , at_dists , rnd_pos, bond_types = readd() bondtyplst = list(bond_types.keys()) syms = [tst_atom.symbol for tst_atom in tst_atoms] @@ -54,8 +56,54 @@ def readd(): def remove_atom(atoms,symbols,tol = 0.5): #TODO Currently this is a copy/paste of the old code, needs work. blmin = closest_distances_generator(atom_numbers=[atomic_numbers[symbol] for symbol in symbols] + [atomic_numbers['Ne']], ratio_of_covalent_radii=0.5) - return atoms + remove_index = vnp.random.choice(len(atoms)) + remove = atoms[remove_index] -def change_cell(): - #TODO Currently this is a copy/paste of the old code, needs work. - return atoms + tst_atoms = atoms.copy() + tst_atoms.pop(remove_index) #remove the randomly selected atom from copy + + tst_atoms.wrap() + + atinds = [atom.index for atom in tst_atoms] + at_dists = {i: [] for i in atinds} + nl = primitive_neighbor_list('ijdD', pbc=tst_atoms.pbc, positions=tst_atoms.positions, + cell = atoms.get_cell(), cutoff=5.0) + bond_types = {i: [] for i in atinds} + for i,j in zip(nl[0], nl[-1]): + at_dists[i].append(j) + for i,j in zip(nl[0],nl[l]): + bond_types[i].append((atomic_numbers[tst_atoms[i].symbol], + atomic_numbers[tst_atoms[j].symbol])) + + for i, atom in enumerate(tst_atoms): + tst_dists = at_dists[i] + tst_bonds = bond_types[i] + conds = all([vnp.linalg.norm(tst_atoms[i].position - test_atoms[j].position) >= + blmin[(atomic_numbers[atom.symbol], tst_bonds[k][1])] - tol for k,j in enumerate(tst_dists)]) + if not conds: + tst_atoms.append(remove) + return atoms + return tst_atoms + + +def change_cell(): #change cell size through scale variable. + cell = atoms.get_cell() #current cell + new_atoms = atoms.copy() + new_cell = cell * scale + new_atoms.set_cell(new_cell, scale_atoms=True) #scale with wanted density + return new_atoms + +#change cell + +#Change density but keep crystals structure polyhedral template to show what map hit is + +#Change volume + +#Cell is the box that has the atoms + +#Change cell - lets them change the ase atoms object - a variable that could scale the object by + +# ase doc note +# --get_cell(complete=False)[source] +#Get the three unit cell vectors as a :ase.cell.Cell` object. +#The Cell object resembles a 3x3 ndarray, and cell[i, j] is the jth Cartesian coordinate of the ith cell vector. diff --git a/GRSlib/motion/lossfunc/moments.py b/GRSlib/motion/lossfunc/moments.py deleted file mode 100644 index 11b5b83..0000000 --- a/GRSlib/motion/lossfunc/moments.py +++ /dev/null @@ -1,158 +0,0 @@ -from GRSlib.motion.scoring import Scoring -import jax.numpy as jnp -import numpy as np -from jax import grad, jit -from functools import partial - -class Moments(Scoring): - def __init__(self, *args): #pt, config, target_desc, prior_desc): - self.pt, self.config, descriptors = args - self.target_desc = descriptors.get('target',None).copy() - self.prior_desc = descriptors.get('prior',None).copy() - self.n_descriptors = np.shape(self.target_desc)[1] - self.mask = list(range(self.n_descriptors)) - self.n_params = 1 #Variables LAMMPS needs to know about - self.n_elements = self.config.sections['BASIS'].numtypes #Variables LAMMPS needs to know about - self.loss_ff_grad = grad(self.construct_loss) - self.grad_loss = grad(self.construct_loss) - self.mode = "update" - - def __call__(self, *args): - #The arguments that this function brings in are super improtant and are expected by LAMMPS MLIAP package. - #They are (elems, current_desc, beta, energy) - #Integer values of LAMMPS atom types are in elems - #Descriptors as a per-atom array into current_desc. - #Per-atom forces are expected for beta - #Per-atom energy is expected for energy, need to do some testing if per-atom values can be reported back. - if self.mode=="score": - elems, current_desc, beta, energy = args - self.n_atoms = np.shape(current_desc)[0] - - #TODO Explain - score = self.construct_loss(current_desc, self.target_desc) - energy[:] = 0 - energy[0] = float(self.config.sections["SCORING"].strength_target)*score #Scaled score (energy) between current and target - forces = self.grad_loss(current_desc, self.target_desc) #Forces between current and target - beta[:,:]= 0 - beta[:,self.mask] = float(self.config.sections["SCORING"].strength_target)*forces #Scaled forces between current and target - - #TODO Explain - score = self.construct_loss(current_desc, self.prior_desc) -# energy[0] += float(self.config.sections["SCORING"].strength_prior)*score #Scaled score (energy) between current and prior -# print(" Target, Prior Scores: ", energy[0], score) - forces = self.grad_loss(current_desc, self.prior_desc) #Forces between current and prior structures - beta[:,self.mask] += float(self.config.sections["SCORING"].strength_prior)*forces #Scaled forces between current and prior - - elif self.mode=="update": - self.update(args) - beta = self.grad_loss(self.target_desc, self.target_desc) - - - def set_mode_update(self): - self.mode="update" - - def set_mode_score(self): - self.mode="score" - - def update(self,*args): - pt, config, descriptors = args[0] - self.target_desc = descriptors.get('target',None).copy() - self.prior_desc = descriptors.get('prior',None).copy() - self.n_descriptors = np.shape(self.target_desc)[1] - self.mask = list(range(self.n_descriptors)) - - if self.n_elements > 1: - self.current_desc = self.current_desc.flatten() - self.target_desc = self.target_desc.flatten() - self.prior_desc = self.prior_desc.flatten() - self.mode = "score" - - @partial(jit, static_argnums=(0,)) - def construct_loss(self, current_desc, target_desc): - #This needs to be a dynamic call like is done in descriptor calcs - loss_ff = 0 - set_of_moments = [] - - if (any(x == 'mean' for x in self.config.sections['SCORING'].moments)): - first_mom = self.first_moment(current_desc, target_desc) - else: - first_mom = None - - if (any(x == 'stdev' for x in self.config.sections['SCORING'].moments)): - second_mom = self.second_moment(current_desc, target_desc) - else: - second_mom = None - - if (any(x == 'skew' for x in self.config.sections['SCORING'].moments)): - third_mom = self.third_moment(current_desc, target_desc) - else: - third_mom = None - - if (any(x == 'kurt' for x in self.config.sections['SCORING'].moments)): - fourth_mom = self.fourth_moment(current_desc, target_desc) - else: - fourth_mom = None - - for item in [first_mom, second_mom, third_mom, fourth_mom]: - if item != None: - set_of_moments.append(item) - - for item in set_of_moments: - loss_ff += item - return loss_ff - - @partial(jit, static_argnums=(0,)) - def first_moment(self, current_desc, target_desc): - current_avg = jnp.average(current_desc, axis=0) - target_avg = jnp.average(target_desc, axis=0) - tst_residual = jnp.sum(jnp.nan_to_num(jnp.abs(current_avg-target_avg))) - is_zero = jnp.array(jnp.isclose(tst_residual,jnp.zeros(tst_residual.shape)),dtype=int) - bonus = -jnp.sum(is_zero*(float(self.config.sections['SCORING'].moments_bonus[0]))) - tst_residual_final = tst_residual*float(self.config.sections['SCORING'].moments_coeff[0]) + bonus #MAE + bonus - return tst_residual_final - - @partial(jit, static_argnums=(0,)) - def second_moment(self, current_desc, target_desc): - current_std = jnp.std(current_desc, axis=0) - target_std = jnp.std(target_desc, axis=0) - tst_residual = jnp.sum(jnp.nan_to_num(jnp.abs(current_std-target_std))) - is_zero = jnp.array(jnp.isclose(tst_residual,jnp.zeros(tst_residual.shape)),dtype=int) - bonus = -jnp.sum(is_zero*float(self.config.sections['SCORING'].moments_bonus[1])) - tst_residual_final = tst_residual*float(self.config.sections['SCORING'].moments_coeff[1]) + bonus #MAE + bonus - return tst_residual_final - - @partial(jit, static_argnums=(0,)) - def third_moment(self, current_desc, target_desc): - #Showing my work for Pearsons skew = (3(mean-median)/stdev)) - current_avg = jnp.average(current_desc, axis=0) - target_avg = jnp.average(target_desc, axis=0) - current_std = jnp.std(current_desc, axis=0) - target_std = jnp.std(target_desc, axis=0) - current_med = jnp.median(current_desc, axis=0) - target_med = jnp.median(target_desc, axis=0) - - current_skew = 3.0*(current_avg-current_med)/current_std - target_skew = 3.0*(target_avg-target_med)/target_std - - tst_residual = jnp.sum(jnp.nan_to_num(jnp.abs(current_skew-target_skew))) - is_zero = jnp.array(jnp.isclose(tst_residual,jnp.zeros(tst_residual.shape)),dtype=int) - bonus = -jnp.sum(is_zero*float(self.config.sections['SCORING'].moments_bonus[2])) - tst_residual_final = tst_residual*float(self.config.sections['SCORING'].moments_coeff[2]) + bonus #MAE + bonus - return tst_residual_final - - @partial(jit, static_argnums=(0,)) - def fourth_moment(self, current_desc, target_desc): - #Showing my work for Kurtosis = Avg(z^4.0)-3 where z=(x-avg(x))/stdev(x) - current_avg = jnp.average(current_desc, axis=0) - target_avg = jnp.average(target_desc, axis=0) - current_std = jnp.std(current_desc, axis=0) - target_std = jnp.std(target_desc, axis=0) - - current_kurt = jnp.average(((current_desc-current_avg)/current_std)**4.0)-3.0 - target_kurt = jnp.average(((target_desc-target_avg)/target_std)**4.0)-3.0 - - tst_residual = jnp.sum(jnp.nan_to_num(jnp.abs(current_kurt-target_kurt))) - is_zero = jnp.array(jnp.isclose(tst_residual,jnp.zeros(tst_residual.shape)),dtype=int) - bonus = -jnp.sum(is_zero*float(self.config.sections['SCORING'].moments_bonus[3])) - tst_residual_final = tst_residual*float(self.config.sections['SCORING'].moments_coeff[3]) + bonus #MAE + bonus - return tst_residual_final diff --git a/GRSlib/motion/scoring.py b/GRSlib/motion/scoring.py index 4d39f2f..597a8a6 100644 --- a/GRSlib/motion/scoring.py +++ b/GRSlib/motion/scoring.py @@ -1,11 +1,35 @@ +from GRSlib.Ver0_Files.opt_tools import internal_generate_cell +from GRSlib.parallel_tools import ParallelTools +#from GRSlib.motion.lossfunc.moments import LossFunction +from GRSlib.Ver0_Files.opt_tools import get_desc_count #from GRSlib.parallel_tools import ParallelTools #from GRSlib.motion.lossfunc.moments import Moments #from GRSlib.motion.lossfunc import Gradient from GRSlib.converters.sections.lammps_base import Base, _extract_compute_np +from examples.simple_test.GRS_protocol import GRSModel, GRSSampler import lammps, lammps.mliap from lammps.mliap.loader import * from functools import partial -import numpy as np +import numpy as vnp +from GRSlib.Ver0_Files.opt_tools import * +from GRSlib.GRS import * +#from examples.simple_test.GRS_protocol import * +#Scoring has to be a class within motion because we want a consistent reference for scores, ans this +#refrence will be LAMMPS using a constructed potential energy surface from the representation loss function +n_totconfig = 10 +data_path = 'bcc.data' +cross_weight =1.000000 +self_weight = 1.000000 +randomize_comps= False # # flag to use randomized compositions for elements in the dictionary: target_comps = {'Cr':1.0 } +mincellsize = 54 +maxcellsize=55 +target_comps = {'W:1.0'} +min_typ_global='box' #box or min +soft_strength=0.0 +elems=get_desc_count('coupling_coefficients.yace',return_elems=True) +nelements= len(elems) +n_descs= get_desc_count('coupling_coefficients.yace') +rand_comp =1 #Scoring has to be a class within motion because we want a consistent reference for scores, and this #refrence will be LAMMPS using a constructed potential energy surface from the representation loss function. @@ -94,4 +118,111 @@ def _extract_commands(self,string): #Can be given a block of text where it will split them into individual commands add_lmp_lines = [x for x in string.splitlines() if x.strip() != ''] for line in add_lmp_lines: - self.lmp.command(line) + self.lmp.command(line) +#internal generate cell is only used in scoring.py but defined in opt_tools + gsqsmodel only in scoring.py currently +#class ensemble_score(): #will take in the target, and compare it to multiple generated structures -- look at fitsnap sections +# min and maxatoms the same as maxcellsize? +# target_comps is in input +#numelements = num types? + def ensemble_score(self, n_totconfig, data_path, cross_weight, self_weight, randomize_comps, mincellsize, maxcellsize, target_comps, min_typ_global, soft_strength, nelements, n_descs, mask, rand_comp): + if mask == None: + mask = range(n_descs) + self.mask = mask # Generates the multiple structures + scores = [] # Initialize a list to store scores + + print(f"Starting ensemble_score with {n_totconfig} configurations.") # Debugging line + + for i in range(1, n_totconfig + 1): + print(f"Configuration {i}/{n_totconfig} - Using indices: {mask}") # Debugging line + + # Generate the cell + g = internal_generate_cell(i, desired_size=vnp.random.choice(range(mincellsize, maxcellsize)), template=None, desired_comps=target_comps, use_template=None, min_typ=min_typ_global, soft_strength=soft_strength) + print(g) + print(f"Cell generated for configuration {i}: {g}") # Debugging line + + em = GRSModel(nelements, n_descs, mask=mask) + sampler = GRSSampler(em, g) + em.K_cross = cross_weight + em.K_self = self_weight + + print(f"Running minimization for configuration {i}.") # Debugging line + # Run the minimization process + sampler.run("minimize 1e-6 1e-6 1000 10000") + + print(f"Minimization completed for configuration {i}. Writing data.") # Debugging line + # Write the data to a file + sampler.run("write_data %s/sample.%i.dat " % (data_path, i)) + + print(f"Data written for configuration {i}. Updating model.") # Debugging line + sampler.update_model() # Updating the model combines generated structures + + # Calculate the score for the current configuration + score = self.get_score(g) # Pass the generated structure to get_score + print(f"Score for configuration {i}: {score}") # Debugging line + + if score is None: + print(f"Score for configuration {i} is None.") # Debugging line + + scores.append(score) # Append the score to the list + + print("Final scores list:", scores) # Debugging line + return scores # Return the list of scores +"""def ensemble_score(self, n_totconfig, data_path, cross_weight, self_weight, randomize_comps, mincellsize, maxcellsize, target_comps, min_typ_global, soft_strength, nelements, n_descs, mask, rand_comp): + self.mask = mask # Generates the multiple structures + scores = [] # Initialize a list to store scores + print(f"Starting ensemble_score with {n_totconfig} configs.") + for i in range(1, n_totconfig + 1): + print(i, "/", n_totconfig, "Using indices:", mask) + + if not randomize_comps: + g = internal_generate_cell(i, desired_size=vnp.random.choice(range(mincellsize, maxcellsize)), template=None, desired_comps=target_comps, use_template=None, min_typ=min_typ_global, soft_strength=soft_strength) + else: + target_comps_rnd = rand_comp(target_comps) # Randomize compositions + g = internal_generate_cell(i, desired_size=vnp.random.choice(range(mincellsize, maxcellsize)), template=None, desired_comps=target_comps_rnd, use_template=None, min_typ=min_typ_global, soft_strength=soft_strength) + + em = GRSModel(nelements, n_descs, mask=mask) + sampler = GRSSampler(em, g) + em.K_cross = cross_weight + em.K_self = self_weight + + # Run the minimization process + sampler.run("minimize 1e-6 1e-6 1000 10000") + + # Write the data to a file + sampler.run("write_data %s/sample.%i.dat " % (data_path, i)) + + sampler.update_model() # Updating the model combines generated structures + + # Calculate the score for the current configuration + score = self.get_score() # Assuming get_score is defined to return the score for the current state + if score is None: + print(f"Score for configuration {i} is None.") + scores.append(score) # Append the score to the list + print("Final scores list:", scores) + return scores # Return the list of scores +""" +""" def ensemble_score(self, n_totconfig, data_path,cross_weight, self_weight, randomize_comps, mincellsize, maxcellsize, target_comps, min_typ_global,soft_strength, nelements,n_descs,mask,rand_comp): + self.mask=mask#generates the multiple structures -- needs internal generate cell , some of these should be defined in the input file like n_totconfig if they choose multiple and the crossweight and self weigths + i = 1 + while i <= n_totconfig: + print(i,"/",n_totconfig,"Using indicies :",mask) + if not randomize_comps: + g = internal_generate_cell(i,desired_size=vnp.random.choice(range(mincellsize,maxcellsize)),template=None,desired_comps=target_comps,use_template=None,min_typ=min_typ_global,soft_strength=soft_strength) + else: + target_comps_rnd = rand_comp(target_comps) #randomize_comp in input? and target_comps in input? + g = internal_generate_cell(i,desired_size=vnp.random.choice(range(mincellsize,maxcellsize)),template=None,desired_comps=target_comps_rnd,use_template=None,min_typ=min_typ_global,soft_strength=soft_strength) + em=GRSModel(nelements,n_descs,mask=mask) + sampler = GRSSampler(em, g) + em.K_cross = cross_weight + em.K_self = self_weight + #min type + sampler.run("minimize 1e-6 1e-6 1000 10000") + + sampler.run("write_data %s/sample.%i.dat " % (data_path, i)) + + sampler.update_model() #updating the model is how all of the generated structures get combined into one. + +""" +# total number of configurations +#n_totconfig = int(sys.argv[2]) + # self.lmp.command(line) diff --git a/README.md b/README.md deleted file mode 100644 index 20f6d5d..0000000 --- a/README.md +++ /dev/null @@ -1,33 +0,0 @@ -# Structure_Generation - -## Requirements -FitSNAP - (post library update) -ASE - v 3.22.1 -lammps - LAMMPS (post mliap update) - -You must add the fitsnap folder to your PYTHONPATH in order to use -the ACE descriptor setup functionalities. Otherwise, you must bring your own -ACE descriptor file (.yace). - -Add the `/path/to/Structure_Generation/lib` folder to your PYTHONPATH - -`export PYTHONPATH=$PYTHONPATH:/path/to/Structure_Generation/lib` - -*NOTE fitsnap must also be in your pythonpath* - -## Parent folder contents - -`lib` - folder containing helper functions and classes for structure generation & optimization routines -`tools` - some extra helper scripts here - may be removed. - -`examples` folder containing examples for local (lammps) and global optimization (genetic algorithm, simulated annealing) - See further details in the README of the examples folder - - -## Definition of GRS - -For theoretical definitions and methods, see the manuscript at: *link_here* -*TODO* update with software-relevant definitions - -## USAGE - -follow procedure in `examples/simple_test` first. diff --git a/testing/coupling_coefficients.yace b/testing/coupling_coefficients.yace deleted file mode 100644 index d5e8e51..0000000 --- a/testing/coupling_coefficients.yace +++ /dev/null @@ -1,29 +0,0 @@ -elements: [W] -E0: [0.000000] -deltaSplineBins: 0.001000 -embeddings: - 0: {ndensity: 1, FS_parameters: [1.0, 1.0], npoti: FinnisSinclair, rho_core_cutoff: 100000, drho_core_cutoff: 250} -bonds: - [0, 0]: {nradmax: 8, lmax: 3, nradbasemax: 8, radbasename: ChebExpCos, radparameters: [1.4], radcoefficients: [[[1, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0]], [[0, 1, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0]], [[0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0]], [[0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0]], [[0, 0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0]], [[0, 0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 1, 0, 0]], [[0, 0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 0, 1, 0]], [[0, 0, 0, 0, 0, 0, 0, 1], [0, 0, 0, 0, 0, 0, 0, 1], [0, 0, 0, 0, 0, 0, 0, 1], [0, 0, 0, 0, 0, 0, 0, 1]]], prehc: 0, lambdahc: 1.4, rcut: 5.5, dcut: 0.01, rcut_in: 0.0, dcut_in: 0.01, inner_cutoff_type: distance} -functions: - 0: - - {mu0: 0, rank: 1, ndensity: 1, num_ms_combs: 1, mus: [0], ns: [1], ls: [0], ms_combs: [0], ctildes: [1.0]} - - {mu0: 0, rank: 1, ndensity: 1, num_ms_combs: 1, mus: [0], ns: [2], ls: [0], ms_combs: [0], ctildes: [1.0]} - - {mu0: 0, rank: 1, ndensity: 1, num_ms_combs: 1, mus: [0], ns: [3], ls: [0], ms_combs: [0], ctildes: [1.0]} - - {mu0: 0, rank: 1, ndensity: 1, num_ms_combs: 1, mus: [0], ns: [4], ls: [0], ms_combs: [0], ctildes: [1.0]} - - {mu0: 0, rank: 1, ndensity: 1, num_ms_combs: 1, mus: [0], ns: [5], ls: [0], ms_combs: [0], ctildes: [1.0]} - - {mu0: 0, rank: 1, ndensity: 1, num_ms_combs: 1, mus: [0], ns: [6], ls: [0], ms_combs: [0], ctildes: [1.0]} - - {mu0: 0, rank: 1, ndensity: 1, num_ms_combs: 1, mus: [0], ns: [7], ls: [0], ms_combs: [0], ctildes: [1.0]} - - {mu0: 0, rank: 1, ndensity: 1, num_ms_combs: 1, mus: [0], ns: [8], ls: [0], ms_combs: [0], ctildes: [1.0]} - - {mu0: 0, rank: 2, ndensity: 1, num_ms_combs: 1, mus: [0, 0], ns: [1, 1], ls: [0, 0], ms_combs: [0, 0], ctildes: [1.0]} - - {mu0: 0, rank: 2, ndensity: 1, num_ms_combs: 3, mus: [0, 0], ns: [1, 1], ls: [1, 1], ms_combs: [-1, 1, 0, 0, 1, -1], ctildes: [0.5773502588272095, -0.5773502588272095, 0.5773502588272095]} - - {mu0: 0, rank: 2, ndensity: 1, num_ms_combs: 5, mus: [0, 0], ns: [1, 1], ls: [2, 2], ms_combs: [-2, 2, -1, 1, 0, 0, 1, -1, 2, -2], ctildes: [0.4472135901451111, -0.4472135901451111, 0.4472135901451111, -0.4472135901451111, 0.4472135901451111]} - - {mu0: 0, rank: 2, ndensity: 1, num_ms_combs: 7, mus: [0, 0], ns: [1, 1], ls: [3, 3], ms_combs: [-3, 3, -2, 2, -1, 1, 0, 0, 1, -1, 2, -2, 3, -3], ctildes: [0.37796446681022644, -0.37796446681022644, 0.37796446681022644, -0.37796446681022644, 0.37796446681022644, -0.37796446681022644, 0.37796446681022644]} - - {mu0: 0, rank: 3, ndensity: 1, num_ms_combs: 1, mus: [0, 0, 0], ns: [1, 1, 1], ls: [0, 0, 0], ms_combs: [0, 0, 0], ctildes: [1.0]} - - {mu0: 0, rank: 3, ndensity: 1, num_ms_combs: 3, mus: [0, 0, 0], ns: [1, 1, 1], ls: [0, 1, 1], ms_combs: [0, -1, 1, 0, 0, 0, 0, 1, -1], ctildes: [0.33333332136784577, -0.33333332136784577, 0.33333332136784577]} - - {mu0: 0, rank: 3, ndensity: 1, num_ms_combs: 5, mus: [0, 0, 0], ns: [1, 1, 1], ls: [0, 2, 2], ms_combs: [0, -2, 2, 0, -1, 1, 0, 0, 0, 0, 1, -1, 0, 2, -2], ctildes: [0.1999999952104794, -0.1999999952104794, 0.1999999952104794, -0.1999999952104794, 0.1999999952104794]} - - {mu0: 0, rank: 3, ndensity: 1, num_ms_combs: 7, mus: [0, 0, 0], ns: [1, 1, 1], ls: [0, 3, 3], ms_combs: [0, -3, 3, 0, -2, 2, 0, -1, 1, 0, 0, 0, 0, 1, -1, 0, 2, -2, 0, 3, -3], ctildes: [0.14285713817113876, -0.14285713817113876, 0.14285713817113876, -0.14285713817113876, 0.14285713817113876, -0.14285713817113876, 0.14285713817113876]} - - {mu0: 0, rank: 3, ndensity: 1, num_ms_combs: 9, mus: [0, 0, 0], ns: [1, 1, 1], ls: [1, 1, 2], ms_combs: [-1, -1, 2, -1, 0, 1, -1, 1, 0, 0, -1, 1, 0, 0, 0, 0, 1, -1, 1, -1, 0, 1, 0, -1, 1, 1, -2], ctildes: [0.1999999952104794, -0.14142135359335306, 0.08164966239889981, -0.14142135359335306, 0.16329932479779963, -0.14142135359335306, 0.08164966239889981, -0.14142135359335306, 0.1999999952104794]} - - {mu0: 0, rank: 3, ndensity: 1, num_ms_combs: 15, mus: [0, 0, 0], ns: [1, 1, 1], ls: [1, 2, 3], ms_combs: [-1, -2, 3, -1, -1, 2, -1, 0, 1, -1, 1, 0, -1, 2, -1, 0, -2, 2, 0, -1, 1, 0, 0, 0, 0, 1, -1, 0, 2, -2, 1, -2, 1, 1, -1, 0, 1, 0, -1, 1, 1, -2, 1, 2, -3], ctildes: [0.14285713817113876, -0.11664237211493589, 0.09035078971800603, -0.06388765326736667, 0.03688555491939627, -0.08247860448643873, 0.10432810273894777, -0.11065666757424353, 0.10432810273894777, -0.08247860448643873, 0.03688555491939627, -0.06388765326736667, 0.09035078971800603, -0.11664237211493589, 0.14285713817113876]} - - {mu0: 0, rank: 3, ndensity: 1, num_ms_combs: 19, mus: [0, 0, 0], ns: [1, 1, 1], ls: [2, 2, 2], ms_combs: [-2, 0, 2, -2, 1, 1, -2, 2, 0, -1, -1, 2, -1, 0, 1, -1, 1, 0, -1, 2, -1, 0, -2, 2, 0, -1, 1, 0, 0, 0, 0, 1, -1, 0, 2, -2, 1, -2, 1, 1, -1, 0, 1, 0, -1, 1, 1, -2, 2, -2, 0, 2, -1, -1, 2, 0, -2], ctildes: [0.106904496560845, -0.13093073535988964, 0.106904496560845, -0.13093073535988964, 0.0534522482804225, 0.0534522482804225, -0.13093073535988964, 0.106904496560845, 0.0534522482804225, -0.106904496560845, 0.0534522482804225, 0.106904496560845, -0.13093073535988964, 0.0534522482804225, 0.0534522482804225, -0.13093073535988964, 0.106904496560845, -0.13093073535988964, 0.106904496560845]} - - {mu0: 0, rank: 3, ndensity: 1, num_ms_combs: 29, mus: [0, 0, 0], ns: [1, 1, 1], ls: [2, 3, 3], ms_combs: [-2, -1, 3, -2, 0, 2, -2, 1, 1, -2, 2, 0, -2, 3, -1, -1, -2, 3, -1, -1, 2, -1, 0, 1, -1, 1, 0, -1, 2, -1, -1, 3, -2, 0, -3, 3, 0, -2, 2, 0, -1, 1, 0, 0, 0, 0, 1, -1, 0, 2, -2, 0, 3, -3, 1, -3, 2, 1, -2, 1, 1, -1, 0, 1, 0, -1, 1, 1, -2, 1, 2, -3, 2, -3, 1, 2, -2, 0, 2, -1, -1, 2, 0, -2, 2, 1, -3], ctildes: [0.058321186057467944, -0.08247860448643873, 0.09035078971800603, -0.08247860448643873, 0.058321186057467944, -0.09221388589046331, 0.07142856908556938, -0.026082025684736942, -0.026082025684736942, 0.07142856908556938, -0.09221388589046331, 0.09221388589046331, -0.0, -0.055328333787121764, 0.07377110983879254, -0.055328333787121764, -0.0, 0.09221388589046331, -0.09221388589046331, 0.07142856908556938, -0.026082025684736942, -0.026082025684736942, 0.07142856908556938, -0.09221388589046331, 0.058321186057467944, -0.08247860448643873, 0.09035078971800603, -0.08247860448643873, 0.058321186057467944]} diff --git a/testing/log.lammps b/testing/log.lammps deleted file mode 100644 index 57a87cd..0000000 --- a/testing/log.lammps +++ /dev/null @@ -1,100 +0,0 @@ -LAMMPS (12 Jun 2025 - Development - 919ef6d369) -OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. - using 1 OpenMP thread(s) per MPI task -Total wall time: 0:00:00 - atom_style atomic - read_data notbcc.data -Reading data file ... - orthogonal box = (-3.1804 -3.1804 -3.1804) to (6.3608 6.3608 6.3608) - 1 by 1 by 1 MPI processor grid - reading atoms ... - 54 atoms - read_data CPU = 0.001 seconds - pair_style mliap model mliappy LATER descriptor ace coupling_coefficients.yace -Loading python model deferred. - pair_coeff * * W - neighbor 2.3 bin - neigh_modify one 10000 - thermo 10 - thermo_style custom step etotal temp press - -Loading python model complete. -run 0 -WARNING: No fixes with time integration, atoms won't move -For more information see https://docs.lammps.org/err0028 (src/verlet.cpp:60) -Neighbor list info ... - update: every = 1 steps, delay = 0 steps, check = yes - max neighbors/atom: 10000, page size: 100000 - master list distance cutoff = 7.8 - ghost atom cutoff = 7.8 - binsize = 3.9, bins = 3 3 3 - 1 neighbor lists, perpetual/occasional/extra = 1 0 0 - (1) pair mliap, perpetual - attributes: full, newton on - pair build: full/bin/atomonly - stencil: full/bin/3d - bin: standard -Per MPI rank memory allocation (min/avg/max) = 3.793 | 3.793 | 3.793 Mbytes - Step TotEng Temp Press - 0 1.4656356e-308 0 -63832.528 -Loop time of 3.073e-06 on 1 procs for 0 steps with 54 atoms - -97.6% CPU use with 1 MPI tasks x 1 OpenMP threads - -MPI task timing breakdown: -Section | min time | avg time | max time |%varavg| %total ---------------------------------------------------------------- -Pair | 0 | 0 | 0 | 0.0 | 0.00 -Neigh | 0 | 0 | 0 | 0.0 | 0.00 -Comm | 0 | 0 | 0 | 0.0 | 0.00 -Output | 0 | 0 | 0 | 0.0 | 0.00 -Modify | 0 | 0 | 0 | 0.0 | 0.00 -Other | | 3.073e-06 | | |100.00 - -Nlocal: 54 ave 54 max 54 min -Histogram: 1 0 0 0 0 0 0 0 0 0 -Nghost: 807 ave 807 max 807 min -Histogram: 1 0 0 0 0 0 0 0 0 0 -Neighs: 0 ave 0 max 0 min -Histogram: 1 0 0 0 0 0 0 0 0 0 -FullNghs: 7262 ave 7262 max 7262 min -Histogram: 1 0 0 0 0 0 0 0 0 0 - -Total # of neighbors = 7262 -Ave neighs/atom = 134.48148 -Neighbor list builds = 0 -Dangerous builds = 0 -run 0 -WARNING: No fixes with time integration, atoms won't move -For more information see https://docs.lammps.org/err0028 (src/verlet.cpp:60) -Per MPI rank memory allocation (min/avg/max) = 3.793 | 3.793 | 3.793 Mbytes - Step TotEng Temp Press - 0 15.577732 0 -63832.528 -Loop time of 1.0616e-05 on 1 procs for 0 steps with 54 atoms - -75.4% CPU use with 1 MPI tasks x 1 OpenMP threads - -MPI task timing breakdown: -Section | min time | avg time | max time |%varavg| %total ---------------------------------------------------------------- -Pair | 0 | 0 | 0 | 0.0 | 0.00 -Neigh | 0 | 0 | 0 | 0.0 | 0.00 -Comm | 0 | 0 | 0 | 0.0 | 0.00 -Output | 0 | 0 | 0 | 0.0 | 0.00 -Modify | 0 | 0 | 0 | 0.0 | 0.00 -Other | | 1.062e-05 | | |100.00 - -Nlocal: 54 ave 54 max 54 min -Histogram: 1 0 0 0 0 0 0 0 0 0 -Nghost: 807 ave 807 max 807 min -Histogram: 1 0 0 0 0 0 0 0 0 0 -Neighs: 0 ave 0 max 0 min -Histogram: 1 0 0 0 0 0 0 0 0 0 -FullNghs: 7262 ave 7262 max 7262 min -Histogram: 1 0 0 0 0 0 0 0 0 0 - -Total # of neighbors = 7262 -Ave neighs/atom = 134.48148 -Neighbor list builds = 0 -Dangerous builds = 0 -Total wall time: 0:00:00 diff --git a/testing/testing.py b/testing/testing.py deleted file mode 100644 index 0727f08..0000000 --- a/testing/testing.py +++ /dev/null @@ -1,36 +0,0 @@ -from time import time -from mpi4py import MPI -from GRSlib.GRS import GRS -import numpy as np -from ase.build import bulk -from ase.io import read,write - -comm = MPI.COMM_WORLD -rank = comm.Get_rank() -nprocs = comm.Get_size() - -settings = "../GenericInput.in" -grs = GRS(settings, comm=comm) - -#testing of io class -#grs.config.view_state() -#----------------------- - -#testing of convert class -#attributes = [attr for attr in dir(grs.convert) if not attr.startswith('__')] -#print("attr of grs.convert:") -#print(attributes) - -#atoms = grs.convert.lammps_to_ase('bcc.data') -#print(atoms) -#file = grs.convert.ase_to_lammps(atoms) -#print(file) -#grs.convert.run_lammps_single('bcc.data') -#current_desc = grs.convert_to_desc(file) -#----------------------- -#grs.genetic_move.tournament_selection(data=None) - -score = grs.get_score('bcc.data') -print(score) -#print("!") -exit()