From fb0d16ed54cdd257c7d54688088988dcd02c76ab Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 21 Dec 2022 10:28:04 -0800 Subject: [PATCH 01/50] allow a0 to be specified as an input in some cases a seems to be far from the ASE reference --- pynta/calculator.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pynta/calculator.py b/pynta/calculator.py index 651effa6..94bf9273 100644 --- a/pynta/calculator.py +++ b/pynta/calculator.py @@ -135,15 +135,15 @@ def add_sella_constraint(cons,d): constructor(**constraint_dict) return -def get_lattice_parameter(metal,surface_type,software,software_kwargs,da=0.1,options={"xatol":1e-4}): +def get_lattice_parameter(metal,surface_type,software,software_kwargs,da=0.1,options={"xatol":1e-4},a0=None): soft = name_to_ase_software(software)(**software_kwargs) def f(a): slab = bulk(metal,surface_type[:3],a=a) slab.calc = soft slab.pbc = (True, True, True) return slab.get_potential_energy() - - a0 = reference_states[chemical_symbols.index(metal)]['a'] + if a0 is None: + a0 = reference_states[chemical_symbols.index(metal)]['a'] avals = np.arange(a0-da,a0+da,0.01) outavals = [] Evals = [] From 5282c58369e4631ae5a769875c0333fb9078d102 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 21 Dec 2022 11:42:16 -0800 Subject: [PATCH 02/50] enable freezing of the first n atoms this should allow slab freezing that is agnostic to what atoms are present --- pynta/tasks.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/pynta/tasks.py b/pynta/tasks.py index 5aa63b65..8aa8f60f 100644 --- a/pynta/tasks.py +++ b/pynta/tasks.py @@ -151,6 +151,11 @@ def run_task(self, fw_spec): sp.set_constraint(FixAtoms( indices=[atom.index for atom in sp if atom.symbol == sym] )) + elif c.split()[0] == "freeze" and c.split()[1] == "up" and c.split()[2] == "to": + n = int(c.split()[3]) + sp.set_constraint(FixAtoms( + indices=list(range(n)) + )) opt_kwargs["trajectory"] = label+".traj" From 7084f22423aefc126da2c93be448fc7cb5475ad6 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 21 Dec 2022 11:42:44 -0800 Subject: [PATCH 03/50] use the element agnostic slab freeze --- pynta/main.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/pynta/main.py b/pynta/main.py index 30b8e50f..44b16c25 100644 --- a/pynta/main.py +++ b/pynta/main.py @@ -331,7 +331,7 @@ def setup_adsorbates(self,initial_guess_finished=False): if self.software != "XTB": constraints = ["freeze up to {}".format(self.freeze_ind)] else: - constraints = ["freeze all "+self.metal] + constraints = ["freeze up to "+str(self.nslab)] else: #gas phase big_slab_ads = structure target_site_num = None #no slab so can't run site analysis @@ -365,7 +365,7 @@ def setup_adsorbates(self,initial_guess_finished=False): optfws2.append(fwopt2) vib_obj_dict = {"software": self.software, "label": adsname, "software_kwargs": software_kwargs, - "constraints": ["freeze all "+self.metal]} + "constraints": ["freeze up to "+str(self.nslab)]} cfw = collect_firework(xyzs,True,["vibrations_firework"],[vib_obj_dict],["vib.json"],[],parents=optfws2,label=adsname) self.adsorbate_fw_dict[adsname] = optfws2 @@ -404,7 +404,7 @@ def setup_adsorbates(self,initial_guess_finished=False): if self.software != "XTB": constraints = ["freeze up to {}".format(self.freeze_ind)] else: - constraints = ["freeze all "+self.metal] + constraints = ["freeze up to "+str(self.nslab)] xyz = os.path.join(prefix_path,str(prefix)+".xyz") init_path = os.path.join(prefix_path,prefix+"_init.xyz") assert os.path.exists(init_path), init_path @@ -424,7 +424,7 @@ def setup_adsorbates(self,initial_guess_finished=False): optfws2.append(fwopt2) vib_obj_dict = {"software": self.software, "label": ad, "software_kwargs": software_kwargs, - "constraints": ["freeze all "+self.metal]} + "constraints": ["freeze up to "+str(self.nslab)]} cfw = collect_firework(xyzs,True,["vibrations_firework"],[vib_obj_dict],["vib.json"],[True,False],parents=optfws2,label=ad) self.adsorbate_fw_dict[ad] = optfws2 @@ -442,11 +442,11 @@ def setup_transition_states(self,adsorbates_finished=False): "run_kwargs": {"fmax" : 0.02, "steps" : 70},"constraints": ["freeze up to {}".format(self.freeze_ind)],"sella":True,"order":1,} else: opt_obj_dict = {"software":self.software,"label":"prefix","socket":self.socket,"software_kwargs":self.software_kwargs_TS, - "run_kwargs": {"fmax" : 0.02, "steps" : 70},"constraints": ["freeze all "+self.metal],"sella":True,"order":1,} + "run_kwargs": {"fmax" : 0.02, "steps" : 70},"constraints": ["freeze up to "+str(self.nslab)],"sella":True,"order":1,} vib_obj_dict = {"software":self.software,"label":"prefix","socket":self.socket,"software_kwargs":self.software_kwargs, - "constraints": ["freeze all "+self.metal]} + "constraints": ["freeze up to "+str(self.nslab)]} IRC_obj_dict = {"software":self.software,"label":"prefix","socket":self.socket,"software_kwargs":self.software_kwargs, - "run_kwargs": {"fmax" : 0.1, "steps" : 70},"constraints":["freeze all "+self.metal]} + "run_kwargs": {"fmax" : 0.1, "steps" : 70},"constraints":["freeze up to "+str(self.nslab)]} for i,rxn in enumerate(self.rxns_dict): ts_path = os.path.join(self.path,"TS"+str(i)) os.makedirs(ts_path) From 56ff0596ada288867f6dfca10019be16e1eec4da Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 21 Dec 2022 12:35:58 -0800 Subject: [PATCH 04/50] extend freeze up to option to all optimizations, vib and irc --- pynta/tasks.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/pynta/tasks.py b/pynta/tasks.py index 8aa8f60f..415c9861 100644 --- a/pynta/tasks.py +++ b/pynta/tasks.py @@ -199,6 +199,11 @@ def run_task(self, fw_spec): sp.set_constraint(FixAtoms( indices=[atom.index for atom in sp if atom.symbol == sym] )) + elif c.split()[0] == "freeze" and c.split()[1] == "up" and c.split()[2] == "to": + n = int(c.split()[3]) + sp.set_constraint(FixAtoms( + indices=list(range(n)) + )) opt = Sella(sp,trajectory=label+".traj",order=order) try: @@ -410,6 +415,12 @@ def run_task(self, fw_spec): indices=[atom.index for atom in sp if atom.symbol == sym] )) indices = [atom.index for atom in sp if atom.symbol != sym] + elif c.split()[0] == "freeze" and c.split()[1] == "up" and c.split()[2] == "to": + n = int(c.split()[3]) + sp.set_constraint(FixAtoms( + indices=list(range(n)) + )) + indices = list(range(n)) vib = Vibrations(sp,indices=indices) vib.run() @@ -759,6 +770,11 @@ def run_task(self, fw_spec): sp.set_constraint(FixAtoms( indices=[atom.index for atom in sp if atom.symbol == sym] )) + elif c.split()[0] == "freeze" and c.split()[1] == "up" and c.split()[2] == "to": + n = int(c.split()[3]) + sp.set_constraint(FixAtoms( + indices=list(range(n)) + )) else: raise ValueError("Could not interpret constraint: {}".format(c)) From ddf433fee18181508b15e148717726cfa6eab551 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 21 Dec 2022 14:50:34 -0800 Subject: [PATCH 05/50] give slab calculations full priority --- pynta/main.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pynta/main.py b/pynta/main.py index 44b16c25..80e19199 100644 --- a/pynta/main.py +++ b/pynta/main.py @@ -130,7 +130,7 @@ def generate_slab(self): if self.software != "XTB": fwslab = optimize_firework(os.path.join(self.path,"slab_init.xyz"),self.software,"slab", opt_method="BFGSLineSearch",socket=self.socket,software_kwargs=self.software_kwargs, - run_kwargs={"fmax" : 0.01},out_path=os.path.join(self.path,"slab.xyz"),constraints=["freeze up to {}".format(self.freeze_ind)]) + run_kwargs={"fmax" : 0.01},out_path=os.path.join(self.path,"slab.xyz"),constraints=["freeze up to {}".format(self.freeze_ind)],priority=1000) wfslab = Workflow([fwslab], name=self.label+"_slab") self.launchpad.add_wf(wfslab) self.launch() From f2a0037d6549e8d99f8c56c69172792b4987b2c3 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 21 Dec 2022 17:45:35 -0800 Subject: [PATCH 06/50] fix bug in constrained vibrational frequency calculation --- pynta/tasks.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pynta/tasks.py b/pynta/tasks.py index 415c9861..e9061bde 100644 --- a/pynta/tasks.py +++ b/pynta/tasks.py @@ -420,7 +420,9 @@ def run_task(self, fw_spec): sp.set_constraint(FixAtoms( indices=list(range(n)) )) - indices = list(range(n)) + indices = [ i for i in range(len(sp)) if i >= n] + else: + raise ValueError vib = Vibrations(sp,indices=indices) vib.run() From 718c3d85cacce281797ed51e40c092e74fc9d873 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 17 Jan 2023 16:22:11 -0800 Subject: [PATCH 07/50] add a non-pbc fallback algorithm for xtb optimizations --- pynta/calculator.py | 167 ++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 163 insertions(+), 4 deletions(-) diff --git a/pynta/calculator.py b/pynta/calculator.py index 94bf9273..65af5a12 100644 --- a/pynta/calculator.py +++ b/pynta/calculator.py @@ -11,6 +11,8 @@ from sella import Sella, Constraints import scipy.optimize as opt import copy +from copy import deepcopy +import itertools def get_energy_atom_bond(atoms,ind1,ind2,k,deq): bd,d = get_distances([atoms.positions[ind1]], [atoms.positions[ind2]], cell=atoms.cell, pbc=atoms.pbc) @@ -87,8 +89,8 @@ def calculate(self, atoms=None, properties=None, system_changes=calculator.all_c self.results["free_energy"] += energy self.results["forces"] += forces -def run_harmonically_forced_xtb(atoms,atom_bond_potentials,site_bond_potentials,nslab,method="GFN1-xTB", - constraints=[]): +def run_harmonically_forced_xtb(atoms,atom_bond_potentials,site_bond_potentials,nslab, + molecule_to_atom_maps=None,ase_to_mol_num=None,method="GFN1-xTB",constraints=[]): """ Optimize TS guess using xTB + harmonic forcing terms determined by atom_bond_potentials and site_bond_potentials """ @@ -115,13 +117,170 @@ def run_harmonically_forced_xtb(atoms,atom_bond_potentials,site_bond_potentials, try: opt.run(fmax=0.02,steps=150) - except Exception as e: - return None,None,None + except Exception as e: #no pbc fallback + return run_harmonically_forced_xtb_no_pbc(atoms,atom_bond_potentials,site_bond_potentials,nslab, + molecule_to_atom_maps=molecule_to_atom_maps,ase_to_mol_num=ase_to_mol_num, + constraints=constraints,method=method,dthresh=4.0) Eharm,Fharm = atoms.calc.get_energy_forces() return atoms,Eharm,Fharm +def run_harmonically_forced_xtb_no_pbc(atoms,atom_bond_potentials,site_bond_potentials,nslab, + molecule_to_atom_maps,ase_to_mol_num=None, + constraints=[],method="GFN1-xTB",dthresh=4.0): + """ + This algorithm extends the slab and selects a section of the slab in which the adsorbates are closest + together, truncates the slab around it and optimizes without pbc based on the truncated slab before + translating the result back to the original periodic slab + + dthresh is the threshold x-y distance a slab atom must be away from an adsorbate atom to be truncated + """ + if ase_to_mol_num is None: #assume only one adsorbate and translate that as a unit + ase_to_mol_num = {i+nslab:0 for i in range(len(atoms) - nslab)} + + if not isinstance(molecule_to_atom_maps[0],list): + molecule_to_atom_maps = [molecule_to_atom_maps] + + target = atoms.cell[0][:2] + atoms.cell[1][:2] + + site_poss = [] + site_inds = [] + site_mols = [] + for site_bond_potential in site_bond_potentials: + ind = site_bond_potential["ind"] + site_pos = site_bond_potential["site_pos"] + mol_ind = ase_to_mol_num[ind] + + if mol_ind not in site_mols: + site_poss.append(site_pos) + site_inds.append(ind) + site_mols.append(mol_ind) + else: + i = site_mols.index(mol_ind) + if np.linalg.norm(target-site_pos[:2]) > np.linalg.norm(target-site_poss[i][:2]): + site_poss[i] = site_pos + site_inds[i] = ind + + trans = get_best_translation(site_poss,atoms.cell) + mol_to_trans = {site_mols[i]: trans[i] for i in range(len(site_mols))} + + new_site_potentials = [] + for i,site_potential in enumerate(site_bond_potentials): + sitep = deepcopy(site_potential) + mol_ind = ase_to_mol_num[site_potential["ind"]] + sitep["site_pos"] = (np.array(sitep["site_pos"]) + mol_to_trans[mol_ind]).tolist() + new_site_potentials.append(sitep) + + + slab = atoms[:nslab] + ad = atoms[nslab:] + bigslab = slab * (2,2,1) + + for ind,molind in ase_to_mol_num.items(): + ad.positions[ind-nslab] += mol_to_trans[molind] + + + bigad = bigslab+ad + + delinds = [] + for i,atom in enumerate(bigslab): + dist = np.inf + for j,a in enumerate(ad): + if np.linalg.norm(atom.position[:2]-a.position[:2]) < dist: + dist = np.linalg.norm(atom.position[:2]-a.position[:2]) + if dist > dthresh: + delinds.append(i) + + for ind in reversed(sorted(delinds)): + del bigad[ind] + + bigad.pbc = (False,False,False) + new_nslab = len(bigad) - len(ad) + + for site_potential in new_site_potentials: + site_potential["ind"] += new_nslab-nslab + + new_atom_bond_potentials = [] + for atom_bond_potential in atom_bond_potentials: + abpot = deepcopy(atom_bond_potential) + abpot["ind1"] += new_nslab-nslab + abpot["ind2"] += new_nslab-nslab + new_atom_bond_potentials.append(abpot) + + new_constraints = [] + for constraint in constraints: + if constraint == "freeze slab": + new_constraints.append(constraint) + continue + else: + c = deepcopy(constraint) + for i in range(len(c["indices"])): + c["indices"][i] += new_nslab-nslab + new_constraints.append(c) + + cons = Constraints(bigad) + for c in new_constraints: + if isinstance(c,dict): + add_sella_constraint(cons,c) + elif c == "freeze slab": + for i,atom in enumerate(bigad): #freeze the slab + if i < new_nslab: + cons.fix_translation(atom.index) + else: + raise ValueError("Constraint {} not understood".format(c)) + + + hfxtb = HarmonicallyForcedXTB(method="GFN1-xTB", + atom_bond_potentials=new_atom_bond_potentials, + site_bond_potentials=new_site_potentials) + + bigad.calc = hfxtb + + opt = Sella(bigad,constraints=cons,trajectory="xtbharm.traj",order=0) + + try: + opt.run(fmax=0.02,steps=150) + except: + return None,None,None + + Eharm,Fharm = bigad.calc.get_energy_forces() + + newad = bigad[new_nslab:] + for ind,molind in ase_to_mol_num.items(): + newad.positions[ind-nslab] -= mol_to_trans[molind] + + outadslab = slab + newad + + outadslab.pbc = (True,True,False) + + return outadslab,Eharm,Fharm + +def get_best_translation(poss,cell): + target = cell[0][:2] + cell[1][:2] + pos2ds = [np.array(pos[:2]) for pos in poss] + translations = [np.zeros(2),cell[0][:2],cell[1][:2],cell[0][:2] + cell[1][:2]] + mindist = np.inf + minpos2ddist = np.inf + tranout = None + for transs in itertools.product(translations,repeat=len(pos2ds)): + dist = 0.0 + pos2ddist = 0.0 + for i in range(len(pos2ds)): + dist += np.linalg.norm(pos2ds[i]+transs[i]-target) + for j in range(i): + pos2ddist += np.linalg.norm(pos2ds[i]+transs[i]-(pos2ds[j]+transs[j])) + if abs(pos2ddist - minpos2ddist) > 0.1 and pos2ddist < minpos2ddist: + tranout = transs + mindist = dist + minpos2ddist = pos2ddist + elif abs(pos2ddist - minpos2ddist) < 0.1 and dist < mindist: + tranout = transs + mindist = dist + minpos2ddist = pos2ddist + + return [np.array([x[0],x[1],0.0]) for x in tranout] + def add_sella_constraint(cons,d): """ construct a constraint from a dictionary that is the input to the constraint From 711988d6a3f44e00a709df4b9d4d95cc2768555e Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 17 Jan 2023 16:23:05 -0800 Subject: [PATCH 08/50] handle weird mono and di atomic species separately if rdkit fails --- pynta/mol.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/pynta/mol.py b/pynta/mol.py index 25dbd60d..420ee8fe 100644 --- a/pynta/mol.py +++ b/pynta/mol.py @@ -37,7 +37,20 @@ def get_desorbed_with_map(mol): return molcopy,out_map def get_conformer(desorbed): - rdmol,rdmap = desorbed.to_rdkit_mol(remove_h=False,return_mapping=True) + try: + rdmol,rdmap = desorbed.to_rdkit_mol(remove_h=False,return_mapping=True) + except Exception as e: + syms = [a.symbol for a in desorbed.atoms] + indmap = {i:i for i in range(len(desorbed.atoms))} + if len(desorbed.atoms) == 1: + atoms = Atoms(syms[0],positions=[(0,0,0)]) + return atoms,indmap + elif len(desorbed.atoms) == 2: + atoms = Atoms(syms[0]+syms[1],positions=[(0,0,0),(1.3,0,0)]) + return atoms,indmap + else: + raise e + indmap = {i:rdmap[a] for i,a in enumerate(desorbed.atoms)} Chem.AllChem.EmbedMultipleConfs(rdmol,numConfs=1,randomSeed=1) conf = rdmol.GetConformer() From bb81b579ef8686174ed807d999b92dfd028f12ad Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 17 Jan 2023 16:23:45 -0800 Subject: [PATCH 09/50] fix xtb optimization calls --- pynta/mol.py | 3 ++- pynta/tasks.py | 21 ++++++++++++++++----- 2 files changed, 18 insertions(+), 6 deletions(-) diff --git a/pynta/mol.py b/pynta/mol.py index 420ee8fe..ce20c4c3 100644 --- a/pynta/mol.py +++ b/pynta/mol.py @@ -115,7 +115,8 @@ def generate_adsorbate_guesses(mol,ads,full_slab,cas,mol_to_atoms_map,metal, for i,geo in enumerate(geos): #freeze bonds for messier first opt geo_out,Eharm,Fharm = run_harmonically_forced_xtb(geo,[],site_bond_params_lists[i],len(full_slab), - method="GFN1-xTB",constraints=constraint_list) + molecule_to_atom_maps=mol_to_atoms_map,ase_to_mol_num=None, + method="GFN1-xTB",constraints=constraint_list) if geo_out: geo_out.calc = None geos_out.append(geo_out) diff --git a/pynta/tasks.py b/pynta/tasks.py index e9061bde..d6efb0aa 100644 --- a/pynta/tasks.py +++ b/pynta/tasks.py @@ -533,13 +533,23 @@ def run_task(self, fw_spec): mols = [mol_dict[name] for name in species_names] ts_dict["mols"] = [mol.to_adjacency_list() for mol in mols] ts_dict["ads_sizes"] = [ads_size(mol) for mol in mols] - ts_dict["template_mol_map"] = get_template_mol_map(reactants,mols) + template_mol_map = get_template_mol_map(reactants,mols) + ts_dict["template_mol_map"] = template_mol_map ts_dict["reverse_names"] = reverse_names ts_dict["molecule_to_atom_maps"] = [{value:key for key,value in gratom_to_molecule_atom_maps[name].items()} for name in species_names] with open(os.path.join(ts_path,"info.json"),'w') as f: json.dump(ts_dict,f) + template_to_ase = {i:get_ase_index(i,template_mol_map,mol_atom_maps,nslab,ads_sizes) for i in range(len(reactants.atoms))} + ase_to_mol_num = {} + for tind,aind in template_to_ase.items(): + if aind: + for i,mol_map in enumerate(template_mol_map): + if tind in mol_map.keys(): + ase_to_mol_num[aind] = i + break + tsstructs = get_unique_TS_structs(adsorbates,species_names,cas,nslab,num_surf_sites,mol_dict, gratom_to_molecule_atom_maps,gratom_to_molecule_surface_atom_maps, facet,metal) @@ -571,7 +581,7 @@ def run_task(self, fw_spec): print("number of TS guesses after filtering by max distance between sites") print(len(out_tsstructs)) - inputs = [ (out_tsstructs[j],new_atom_bond_potential_lists[j],new_site_bond_potential_lists[j],nslab,new_constraint_lists[j],ts_path,j) for j in range(len(out_tsstructs))] + inputs = [ (out_tsstructs[j],new_atom_bond_potential_lists[j],new_site_bond_potential_lists[j],nslab,new_constraint_lists[j],ts_path,j,molecule_to_atom_maps,ase_to_mol_num) for j in range(len(out_tsstructs))] #with mp.Pool(nprocs) as pool: # outputs = pool.map(map_harmonically_forced_xtb,inputs) @@ -821,10 +831,11 @@ def run_task(self, fw_spec): return FWAction() def map_harmonically_forced_xtb(input): - tsstruct,atom_bond_potentials,site_bond_potentials,nslab,constraints,ts_path,j = input + tsstruct,atom_bond_potentials,site_bond_potentials,nslab,constraints,ts_path,j,molecule_to_atom_maps,ase_to_mol_num = input os.makedirs(os.path.join(ts_path,str(j))) - sp,Eharm,Fharm = run_harmonically_forced_xtb(tsstruct,atom_bond_potentials,site_bond_potentials,nslab,method="GFN1-xTB", - constraints=constraints) + sp,Eharm,Fharm = run_harmonically_forced_xtb(tsstruct,atom_bond_potentials,site_bond_potentials,nslab, + molecule_to_atom_maps=molecule_to_atom_maps,ase_to_mol_num=ase_to_mol_num, + method="GFN1-xTB",constraints=constraints) if sp: if "initial_charges" in sp.arrays.keys(): #avoid bug in ase From de2baac385be0841d7a870a5586cc32c43910698 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 17 Jan 2023 16:23:57 -0800 Subject: [PATCH 10/50] import numpy in utils.py --- pynta/utils.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pynta/utils.py b/pynta/utils.py index 4610b799..065471ea 100644 --- a/pynta/utils.py +++ b/pynta/utils.py @@ -5,6 +5,7 @@ from ase.io import write, read from copy import deepcopy from importlib import import_module +import numpy as np def get_unique_sym(geoms): ''' Check for the symmetry equivalent structures in the given files From ee8619c67e6534c194b219e8a4f4a2c985579386 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 17 Jan 2023 16:25:12 -0800 Subject: [PATCH 11/50] remove the bond relaxation xtb optimization for bidentates on Li with the non-pbc fallback this tends to ruin lots of formerly not bad guesses --- pynta/mol.py | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/pynta/mol.py b/pynta/mol.py index ce20c4c3..85fbaaa0 100644 --- a/pynta/mol.py +++ b/pynta/mol.py @@ -152,14 +152,7 @@ def generate_adsorbate_guesses(mol,ads,full_slab,cas,mol_to_atoms_map,metal, xyzsout.append(geos_out[Eind]) site_bond_params_lists_final.append(site_bond_params_lists_out[Eind]) - xyzfinals = [] - for i,xyz in enumerate(xyzsout): - #unfreeze bonds (should be better than rdkit's guess) - geo_out,Eharm,Fharm = run_harmonically_forced_xtb(xyz,[],site_bond_params_lists_final[i],len(full_slab), - method="GFN1-xTB",constraints=["freeze slab"]) - xyzfinals.append(geo_out) - - return xyzfinals + return xyzsout site_bond_length_dict = { ("ontop",None,None): 1.826370311, From 02256b81b85c77a56dfe2b6d56fe1710d3ef69a1 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 17 Jan 2023 16:25:48 -0800 Subject: [PATCH 12/50] handle fold and longbridge site types --- pynta/mol.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pynta/mol.py b/pynta/mol.py index 85fbaaa0..4722c9e5 100644 --- a/pynta/mol.py +++ b/pynta/mol.py @@ -190,6 +190,10 @@ def generate_adsorbate_guesses(mol,ads,full_slab,cas,mol_to_atoms_map,metal, } def get_site_bond_length(sitetype,atomtype=None,metal=None): + if "fold" in sitetype: + sitetype = "fcc" + if sitetype == "longbridge" or sitetype == "shortbridge": + sitetype = "bridge" if (sitetype,atomtype,metal) in site_bond_length_dict.keys(): return site_bond_length_dict[(sitetype,atomtype,metal)] elif (sitetype,atomtype,None) in site_bond_length_dict.keys(): From 7610862b94ac6129b0e0ad8e89593f60d44e7964 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 1 Feb 2023 17:26:20 -0800 Subject: [PATCH 13/50] add advanced get_unique_sites function mostly copied from latest acat --- pynta/mol.py | 68 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) diff --git a/pynta/mol.py b/pynta/mol.py index 4722c9e5..8e6c1360 100644 --- a/pynta/mol.py +++ b/pynta/mol.py @@ -325,6 +325,74 @@ def place_adsorbate(ads,slab,atom_surf_inds,sites,metal): else: raise ValueError +def get_unique_sites(cas, unique_composition=False, + unique_subsurf=False, + return_signatures=False, + return_site_indices=False, + about=None, site_list=None): + """Function mostly copied from the ACAT software + Get all symmetry-inequivalent adsorption sites (one + site for each type). + + Parameters + ---------- + unique_composition : bool, default False + Take site composition into consideration when + checking uniqueness. + + unique_subsurf : bool, default False + Take subsurface element into consideration when + checking uniqueness. + + return_signatures : bool, default False + Whether to return the unique signatures of the + sites instead. + + return_site_indices: bool, default False + Whether to return the indices of each unique + site (in the site list). + + about: numpy.array, default None + If specified, returns unique sites closest to + this reference position. + + """ + + if site_list is None: + sl = cas.site_list + else: + sl = site_list + key_list = ['site', 'morphology'] + if unique_composition: + if not cas.composition_effect: + raise ValueError('the site list does not include ' + + 'information of composition') + key_list.append('composition') + if unique_subsurf: + key_list.append('subsurf_element') + else: + if unique_subsurf: + raise ValueError('to include the subsurface element, ' + + 'unique_composition also need to be set to True') + if return_signatures: + sklist = sorted([[s[k] for k in key_list] for s in sl]) + return sorted(list(sklist for sklist, _ in groupby(sklist))) + else: + seen_tuple = [] + uni_sites = [] + if about is not None: + sl = sorted(sl, key=lambda x: get_mic(x['position'], + about, cas.cell, return_squared_distance=True)) + for i, s in enumerate(sl): + sig = tuple(s[k] for k in key_list) + if sig not in seen_tuple: + seen_tuple.append(sig) + if return_site_indices: + s = i + uni_sites.append(s) + + return uni_sites + def generate_unique_site_additions(geo,cas,nslab,site_bond_params_list=[],sites_list=[]): nads = len(geo) - nslab #label sites with unique noble gas atoms From 98a638c11cb1e8d851e8030805d98ce3ec690e7e Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 1 Feb 2023 17:27:34 -0800 Subject: [PATCH 14/50] add function for selecting unique sites and site pairs based on fingerprinting uses the site identity and morphology for each site and the xy distance and the signed z distance --- pynta/mol.py | 45 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/pynta/mol.py b/pynta/mol.py index 8e6c1360..e503757b 100644 --- a/pynta/mol.py +++ b/pynta/mol.py @@ -2,6 +2,7 @@ from ase.io import read, write from ase.data import covalent_radii from ase import Atoms +from ase.geometry import get_distances from acat.adsorption_sites import SlabAdsorptionSites from acat.adsorbate_coverage import SlabAdsorbateCoverage from acat.settings import site_heights @@ -393,6 +394,50 @@ def get_unique_sites(cas, unique_composition=False, return uni_sites +def generate_unique_placements(slab,cas): + nslab = len(slab) + middle = sum(slab.cell)/2.0 + + unique_single_sites = get_unique_sites(cas,about=middle) + + unique_site_pairs = dict() # (site1site,site1morph,),(site2site,site2morph),xydist,zdist + for unique_site in unique_single_sites: + uni_site_fingerprint = (unique_site["site"],unique_site["morphology"]) + for site in cas.get_sites(): + site_fingerprint = (site["site"],site["morphology"]) + bd,d = get_distances([unique_site["position"]], [site["position"]], cell=slab.cell, pbc=(True,True,False)) + xydist = np.linalg.norm(bd[0][0][:1]) + zdist = bd[0][0][2] + + fingerprint = (uni_site_fingerprint,site_fingerprint,round(xydist,3),round(zdist,3)) + + if fingerprint in unique_site_pairs.keys(): + current_sites = unique_site_pairs[fingerprint] + current_dist = np.linalg.norm(sum([s["position"][:1] for s in current_sites])/2-middle[:1]) + possible_dist = np.linalg.norm((unique_site["position"][:1]+site["position"][:1])/2-middle[:1]) + if possible_dist < current_dist: + unique_site_pairs[fingerprint] = [unique_site,site] + else: + unique_site_pairs[fingerprint] = [unique_site,site] + + unique_site_pairs_lists = list(unique_site_pairs.values()) + unique_site_lists = [[unique_site] for unique_site in unique_single_sites] + + single_site_bond_params_lists = [] + for unique_site_list in unique_site_lists: + pos = unique_site_list[0]["position"] + single_site_bond_params_lists.append([{"site_pos": pos,"ind": None, "k": 100.0, "deq": 0.0}]) + + double_site_bond_params_lists = [] + for unique_site_pair_list in unique_site_pairs_lists: + bond_params_list = [] + for site in unique_site_pair_list: + pos = site["position"] + bond_params_list.append({"site_pos": pos,"ind": None, "k": 100.0, "deq": 0.0}) + double_site_bond_params_lists.append(bond_params_list) + + return unique_site_lists,unique_site_pairs_lists,single_site_bond_params_lists,double_site_bond_params_lists + def generate_unique_site_additions(geo,cas,nslab,site_bond_params_list=[],sites_list=[]): nads = len(geo) - nslab #label sites with unique noble gas atoms From 2125f032ca4ebf32f837ad5763536ad9c2f2ebc6 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 1 Feb 2023 17:28:03 -0800 Subject: [PATCH 15/50] switch to using generate_unique_placements for slab analysis --- pynta/main.py | 18 ++++-------------- 1 file changed, 4 insertions(+), 14 deletions(-) diff --git a/pynta/main.py b/pynta/main.py index 80e19199..10ec7532 100644 --- a/pynta/main.py +++ b/pynta/main.py @@ -149,22 +149,12 @@ def analyze_slab(self): self.cas = cas - single_geoms,single_site_bond_params_lists,single_sites_lists = generate_unique_site_additions(full_slab,cas,len(full_slab),site_bond_params_list=[],sites_list=[]) - - double_geoms_full = [] - double_site_bond_params_lists_full = [] - double_sites_lists_full = [] - for i in range(len(single_geoms)): - double_geoms,double_site_bond_params_lists,double_sites_lists = generate_unique_site_additions(single_geoms[i], - cas,len(full_slab),single_site_bond_params_lists[i],single_sites_lists[i]) - double_geoms_full.extend(double_geoms) - double_site_bond_params_lists_full.extend(double_site_bond_params_lists) - double_sites_lists_full.extend(double_sites_lists) + unique_site_lists,unique_site_pairs_lists,single_site_bond_params_lists,double_site_bond_params_lists = generate_unique_placements(full_slab,cas) self.single_site_bond_params_lists = single_site_bond_params_lists - self.single_sites_lists = single_sites_lists - self.double_site_bond_params_lists = double_site_bond_params_lists_full - self.double_sites_lists = double_sites_lists_full + self.single_sites_lists = unique_site_lists + self.double_site_bond_params_lists = double_site_bond_params_lists + self.double_sites_lists = unique_site_pairs_lists def generate_mol_dict(self): """ From 658f618e9820a1a31c9f6366a5cd84dd9b0dbe18 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 1 Feb 2023 17:28:29 -0800 Subject: [PATCH 16/50] enable fmaxopt and fmaxirc specification in Pynta --- pynta/main.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/pynta/main.py b/pynta/main.py index 10ec7532..a75c7155 100644 --- a/pynta/main.py +++ b/pynta/main.py @@ -1,5 +1,5 @@ from pynta.tasks import * -from pynta.mol import get_adsorbate, generate_unique_site_additions, generate_adsorbate_guesses, get_name +from pynta.mol import get_adsorbate, generate_unique_site_additions, generate_adsorbate_guesses, get_name,generate_unique_placements from molecule.molecule import Molecule import ase.build from ase.io import read, write @@ -35,7 +35,7 @@ def __init__(self,path,rxns_file,surface_type,metal,label,launchpad_path=None,fw TS_opt_software_kwargs=None, lattice_opt_software_kwargs={'kpts': (25,25,25), 'ecutwfc': 70, 'degauss':0.02, 'mixing_mode': 'plain'}, reset_launchpad=False,queue_adapter_path=None,num_jobs=25,max_num_hfsp_opts=None,#max_num_hfsp_opts is mostly for fast testing - Eharmtol=3.0,Eharmfiltertol=30.0,Ntsmin=5,frozen_layers=2): + Eharmtol=3.0,Eharmfiltertol=30.0,Ntsmin=5,frozen_layers=2,fmaxopt=0.05,fmaxirc=0.1): self.surface_type = surface_type if launchpad_path: @@ -108,6 +108,8 @@ def __init__(self,path,rxns_file,surface_type,metal,label,launchpad_path=None,fw self.Eharmfiltertol = Eharmfiltertol self.Ntsmin = Ntsmin self.max_num_hfsp_opts = max_num_hfsp_opts + self.fmaxopt = fmaxopt + self.fmaxirc = fmaxirc def generate_slab(self): """ @@ -130,7 +132,7 @@ def generate_slab(self): if self.software != "XTB": fwslab = optimize_firework(os.path.join(self.path,"slab_init.xyz"),self.software,"slab", opt_method="BFGSLineSearch",socket=self.socket,software_kwargs=self.software_kwargs, - run_kwargs={"fmax" : 0.01},out_path=os.path.join(self.path,"slab.xyz"),constraints=["freeze up to {}".format(self.freeze_ind)],priority=1000) + run_kwargs={"fmax" : self.fmaxopt},out_path=os.path.join(self.path,"slab.xyz"),constraints=["freeze up to {}".format(self.freeze_ind)],priority=1000) wfslab = Workflow([fwslab], name=self.label+"_slab") self.launchpad.add_wf(wfslab) self.launch() @@ -348,7 +350,7 @@ def setup_adsorbates(self,initial_guess_finished=False): fwopt2 = optimize_firework(os.path.join(self.path,"Adsorbates",adsname,str(prefix),"weakopt_"+str(prefix)+".xyz"), self.software,str(prefix), opt_method="QuasiNewton",socket=self.socket,software_kwargs=software_kwargs, - run_kwargs={"fmax" : 0.02, "steps" : 70},parents=[fwopt],constraints=constraints, + run_kwargs={"fmax" : self.fmaxopt, "steps" : 70},parents=[fwopt],constraints=constraints, ignore_errors=True, metal=self.metal, facet=self.surface_type, target_site_num=target_site_num, priority=3) optfws.append(fwopt) optfws.append(fwopt2) @@ -407,7 +409,7 @@ def setup_adsorbates(self,initial_guess_finished=False): fwopt2 = optimize_firework(os.path.join(self.path,"Adsorbates",ad,str(prefix),"weakopt_"+str(prefix)+".xyz"), self.software,str(prefix), opt_method="QuasiNewton",socket=self.socket,software_kwargs=software_kwargs, - run_kwargs={"fmax" : 0.02, "steps" : 70},parents=[fwopt],constraints=constraints, + run_kwargs={"fmax" : self.fmaxopt, "steps" : 70},parents=[fwopt],constraints=constraints, ignore_errors=True, metal=self.metal, facet=self.surface_type, target_site_num=target_site_num, priority=3) optfws.append(fwopt) optfws.append(fwopt2) @@ -429,14 +431,14 @@ def setup_transition_states(self,adsorbates_finished=False): """ if self.software != "XTB": opt_obj_dict = {"software":self.software,"label":"prefix","socket":self.socket,"software_kwargs":self.software_kwargs_TS, - "run_kwargs": {"fmax" : 0.02, "steps" : 70},"constraints": ["freeze up to {}".format(self.freeze_ind)],"sella":True,"order":1,} + "run_kwargs": {"fmax" : self.fmaxopt, "steps" : 70},"constraints": ["freeze up to {}".format(self.freeze_ind)],"sella":True,"order":1,} else: opt_obj_dict = {"software":self.software,"label":"prefix","socket":self.socket,"software_kwargs":self.software_kwargs_TS, "run_kwargs": {"fmax" : 0.02, "steps" : 70},"constraints": ["freeze up to "+str(self.nslab)],"sella":True,"order":1,} vib_obj_dict = {"software":self.software,"label":"prefix","socket":self.socket,"software_kwargs":self.software_kwargs, "constraints": ["freeze up to "+str(self.nslab)]} IRC_obj_dict = {"software":self.software,"label":"prefix","socket":self.socket,"software_kwargs":self.software_kwargs, - "run_kwargs": {"fmax" : 0.1, "steps" : 70},"constraints":["freeze up to "+str(self.nslab)]} + "run_kwargs": {"fmax" : self.fmaxirc, "steps" : 70},"constraints":["freeze up to "+str(self.nslab)]} for i,rxn in enumerate(self.rxns_dict): ts_path = os.path.join(self.path,"TS"+str(i)) os.makedirs(ts_path) From 0a6d41f2127c6a0901f1dbc8c4bc3a5ab5f21975 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 1 Feb 2023 17:30:16 -0800 Subject: [PATCH 17/50] add fmaxopt and fmaxirc to documentation --- docsrc/README.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docsrc/README.rst b/docsrc/README.rst index 1d084a6f..1aa6538f 100644 --- a/docsrc/README.rst +++ b/docsrc/README.rst @@ -311,6 +311,10 @@ specified in terms of what `ASE `_ expects. **lattice_opt_software_kwargs**: this is a dictionary of keyword arguments that should be different from software_kwargs when optimizing the lattice constant +**fmaxopt**: this is the fmax used for all slab, adsorbate and TS geometry optimization jobs + +**fmaxirc**: this is the fmax used for IRC calculations + Lastly there are the Pynta energy filter criteria: **Eharmtol**: a tolerance such that all TS/adsorbate guesses are always calculated if they have energies less than Emin * Eharmtol From b0763424545292e53f59097e9d2c027ed8ee4233 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Thu, 2 Feb 2023 16:23:04 -0800 Subject: [PATCH 18/50] fix variable name --- pynta/tasks.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pynta/tasks.py b/pynta/tasks.py index d6efb0aa..a877514e 100644 --- a/pynta/tasks.py +++ b/pynta/tasks.py @@ -541,7 +541,7 @@ def run_task(self, fw_spec): with open(os.path.join(ts_path,"info.json"),'w') as f: json.dump(ts_dict,f) - template_to_ase = {i:get_ase_index(i,template_mol_map,mol_atom_maps,nslab,ads_sizes) for i in range(len(reactants.atoms))} + template_to_ase = {i:get_ase_index(i,template_mol_map,gratom_to_molecule_atom_maps,nslab,ads_sizes) for i in range(len(reactants.atoms))} ase_to_mol_num = {} for tind,aind in template_to_ase.items(): if aind: From 98d557589a1b71dbcc446f5f69efe274c2f2c157 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Thu, 2 Feb 2023 16:37:35 -0800 Subject: [PATCH 19/50] other definition fix --- pynta/tasks.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pynta/tasks.py b/pynta/tasks.py index a877514e..b063e61e 100644 --- a/pynta/tasks.py +++ b/pynta/tasks.py @@ -541,7 +541,7 @@ def run_task(self, fw_spec): with open(os.path.join(ts_path,"info.json"),'w') as f: json.dump(ts_dict,f) - template_to_ase = {i:get_ase_index(i,template_mol_map,gratom_to_molecule_atom_maps,nslab,ads_sizes) for i in range(len(reactants.atoms))} + template_to_ase = {i:get_ase_index(i,template_mol_map,gratom_to_molecule_atom_maps,nslab,[ads_size(mol) for mol in mols]) for i in range(len(reactants.atoms))} ase_to_mol_num = {} for tind,aind in template_to_ase.items(): if aind: From fded9f00af6bb2cd216957b18ad511507a163331 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sun, 12 Feb 2023 09:36:03 -0800 Subject: [PATCH 20/50] fix launching slab job --- pynta/main.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pynta/main.py b/pynta/main.py index a75c7155..cba2f59a 100644 --- a/pynta/main.py +++ b/pynta/main.py @@ -135,7 +135,7 @@ def generate_slab(self): run_kwargs={"fmax" : self.fmaxopt},out_path=os.path.join(self.path,"slab.xyz"),constraints=["freeze up to {}".format(self.freeze_ind)],priority=1000) wfslab = Workflow([fwslab], name=self.label+"_slab") self.launchpad.add_wf(wfslab) - self.launch() + self.launch(single_job=True) while not os.path.exists(self.slab_path): #wait until slab optimizes, this is required anyway and makes the rest of the code simpler time.sleep(1) self.slab = read(self.slab_path) @@ -459,13 +459,13 @@ def setup_transition_states(self,adsorbates_finished=False): fw = Firework([ts_task],parents=parents,name="TS"+str(i)+"est",spec={"_priority": 10}) self.fws.append(fw) - def launch(self): + def launch(self,single_job=False): """ Call appropriate rapidfire function """ if self.queue: rapidfirequeue(self.launchpad,self.fworker,self.qadapter,njobs_queue=self.njobs_queue,nlaunches="infinite") - elif not self.queue and self.num_jobs == 1: + elif not self.queue and (self.num_jobs == 1 or single_job): rapidfire(self.launchpad,self.fworker,nlaunches="infinite") else: launch_multiprocess(self.launchpad,self.fworker,"INFO","infinite",self.num_jobs,5) From 194eb773e458dd14e8837390c6b672ee16b020ce Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sun, 12 Feb 2023 09:43:13 -0800 Subject: [PATCH 21/50] allow generate slab to be run independently --- pynta/main.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pynta/main.py b/pynta/main.py index cba2f59a..47d1f840 100644 --- a/pynta/main.py +++ b/pynta/main.py @@ -111,7 +111,7 @@ def __init__(self,path,rxns_file,surface_type,metal,label,launchpad_path=None,fw self.fmaxopt = fmaxopt self.fmaxirc = fmaxirc - def generate_slab(self): + def generate_slab(self,async=False): """ generates and optimizes a small scale slab that can be scaled to a large slab as needed optimization occurs through fireworks and this process waits until the optimization is completed @@ -136,6 +136,8 @@ def generate_slab(self): wfslab = Workflow([fwslab], name=self.label+"_slab") self.launchpad.add_wf(wfslab) self.launch(single_job=True) + if async: + return while not os.path.exists(self.slab_path): #wait until slab optimizes, this is required anyway and makes the rest of the code simpler time.sleep(1) self.slab = read(self.slab_path) From 824a06e15fa7abf2026233b261d9df71f05830b7 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sun, 12 Feb 2023 09:46:16 -0800 Subject: [PATCH 22/50] rename async --- pynta/main.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pynta/main.py b/pynta/main.py index 47d1f840..79dd7f6c 100644 --- a/pynta/main.py +++ b/pynta/main.py @@ -111,7 +111,7 @@ def __init__(self,path,rxns_file,surface_type,metal,label,launchpad_path=None,fw self.fmaxopt = fmaxopt self.fmaxirc = fmaxirc - def generate_slab(self,async=False): + def generate_slab(self,skip_sync=False): """ generates and optimizes a small scale slab that can be scaled to a large slab as needed optimization occurs through fireworks and this process waits until the optimization is completed @@ -136,7 +136,7 @@ def generate_slab(self,async=False): wfslab = Workflow([fwslab], name=self.label+"_slab") self.launchpad.add_wf(wfslab) self.launch(single_job=True) - if async: + if skip_sync: return while not os.path.exists(self.slab_path): #wait until slab optimizes, this is required anyway and makes the rest of the code simpler time.sleep(1) From 08ca551244e4fdc779d57cbfab81610573ddfb0f Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sun, 12 Feb 2023 09:50:24 -0800 Subject: [PATCH 23/50] let it skip the slab launch --- pynta/main.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pynta/main.py b/pynta/main.py index 79dd7f6c..eab2a53c 100644 --- a/pynta/main.py +++ b/pynta/main.py @@ -111,7 +111,7 @@ def __init__(self,path,rxns_file,surface_type,metal,label,launchpad_path=None,fw self.fmaxopt = fmaxopt self.fmaxirc = fmaxirc - def generate_slab(self,skip_sync=False): + def generate_slab(self,skip_launch=False): """ generates and optimizes a small scale slab that can be scaled to a large slab as needed optimization occurs through fireworks and this process waits until the optimization is completed @@ -135,9 +135,9 @@ def generate_slab(self,skip_sync=False): run_kwargs={"fmax" : self.fmaxopt},out_path=os.path.join(self.path,"slab.xyz"),constraints=["freeze up to {}".format(self.freeze_ind)],priority=1000) wfslab = Workflow([fwslab], name=self.label+"_slab") self.launchpad.add_wf(wfslab) - self.launch(single_job=True) - if skip_sync: + if skip_launch: return + self.launch(single_job=True) while not os.path.exists(self.slab_path): #wait until slab optimizes, this is required anyway and makes the rest of the code simpler time.sleep(1) self.slab = read(self.slab_path) From 69ef7b190c28811848a233e717b61b98690f6665 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 21 Feb 2023 16:17:40 -0800 Subject: [PATCH 24/50] deepcopy position vectors to prevent bug --- pynta/mol.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pynta/mol.py b/pynta/mol.py index e503757b..b80d2bc1 100644 --- a/pynta/mol.py +++ b/pynta/mol.py @@ -425,14 +425,14 @@ def generate_unique_placements(slab,cas): single_site_bond_params_lists = [] for unique_site_list in unique_site_lists: - pos = unique_site_list[0]["position"] + pos = deepcopy(unique_site_list[0]["position"]) single_site_bond_params_lists.append([{"site_pos": pos,"ind": None, "k": 100.0, "deq": 0.0}]) double_site_bond_params_lists = [] for unique_site_pair_list in unique_site_pairs_lists: bond_params_list = [] for site in unique_site_pair_list: - pos = site["position"] + pos = deepcopy(site["position"]) bond_params_list.append({"site_pos": pos,"ind": None, "k": 100.0, "deq": 0.0}) double_site_bond_params_lists.append(bond_params_list) From 8e8d0c96410cd4bd616cc5295b40e7cf48fb15dd Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Thu, 2 Mar 2023 12:20:07 -0800 Subject: [PATCH 25/50] revive fhardmax and allow fizzled parents for the second opt job and collect jobs --- pynta/main.py | 6 ++++-- pynta/tasks.py | 12 ++++++------ 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/pynta/main.py b/pynta/main.py index eab2a53c..5243c724 100644 --- a/pynta/main.py +++ b/pynta/main.py @@ -35,7 +35,7 @@ def __init__(self,path,rxns_file,surface_type,metal,label,launchpad_path=None,fw TS_opt_software_kwargs=None, lattice_opt_software_kwargs={'kpts': (25,25,25), 'ecutwfc': 70, 'degauss':0.02, 'mixing_mode': 'plain'}, reset_launchpad=False,queue_adapter_path=None,num_jobs=25,max_num_hfsp_opts=None,#max_num_hfsp_opts is mostly for fast testing - Eharmtol=3.0,Eharmfiltertol=30.0,Ntsmin=5,frozen_layers=2,fmaxopt=0.05,fmaxirc=0.1): + Eharmtol=3.0,Eharmfiltertol=30.0,Ntsmin=5,frozen_layers=2,fmaxopt=0.05,fmaxirc=0.1,fmaxopthard=0.05): self.surface_type = surface_type if launchpad_path: @@ -110,6 +110,7 @@ def __init__(self,path,rxns_file,surface_type,metal,label,launchpad_path=None,fw self.max_num_hfsp_opts = max_num_hfsp_opts self.fmaxopt = fmaxopt self.fmaxirc = fmaxirc + self.fmaxopthard = fmaxopthard def generate_slab(self,skip_launch=False): """ @@ -353,7 +354,8 @@ def setup_adsorbates(self,initial_guess_finished=False): self.software,str(prefix), opt_method="QuasiNewton",socket=self.socket,software_kwargs=software_kwargs, run_kwargs={"fmax" : self.fmaxopt, "steps" : 70},parents=[fwopt],constraints=constraints, - ignore_errors=True, metal=self.metal, facet=self.surface_type, target_site_num=target_site_num, priority=3) + ignore_errors=True, metal=self.metal, facet=self.surface_type, target_site_num=target_site_num, priority=3, fmaxhard=self.fmaxopthard, + allow_fizzled_parents=True) optfws.append(fwopt) optfws.append(fwopt2) optfws2.append(fwopt2) diff --git a/pynta/tasks.py b/pynta/tasks.py index b063e61e..94e45d97 100644 --- a/pynta/tasks.py +++ b/pynta/tasks.py @@ -58,7 +58,7 @@ def run_task(self, fw_spec): def optimize_firework(xyz,software,label,opt_method=None,sella=None,socket=False,order=0,software_kwargs={},opt_kwargs={}, run_kwargs={},constraints=[],parents=[],out_path=None,time_limit_hrs=np.inf,fmaxhard=0.0,ignore_errors=False, - target_site_num=None,metal=None,facet=None,priority=1): + target_site_num=None,metal=None,facet=None,priority=1,allow_fizzled_parents=False): d = {"xyz" : xyz, "software" : software,"label" : label} if opt_method: d["opt_method"] = opt_method if software_kwargs: d["software_kwargs"] = software_kwargs @@ -81,7 +81,7 @@ def optimize_firework(xyz,software,label,opt_method=None,sella=None,socket=False if out_path is None: out_path = os.path.join(directory,label+".xyz") t2 = FileTransferTask({'files': [{'src': label+'.xyz', 'dest': out_path}, {'src': label+'.traj', 'dest': os.path.join(directory,label+".traj")}], 'mode': 'copy', 'ignore_errors' : ignore_errors}) - return Firework([t1,t2],parents=parents,name=label+"opt",spec={"_priority": priority}) + return Firework([t1,t2],parents=parents,name=label+"opt",spec={"_allow_fizzled_parents": allow_fizzled_parents,"_priority": priority}) @explicit_serialize class MolecularOptimizationTask(OptimizationTask): @@ -612,16 +612,16 @@ def run_task(self, fw_spec): ctask = MolecularCollect({"xyzs":xyzsout,"check_symm":True,"fw_generators": ["optimize_firework",["vibrations_firework","IRC_firework","IRC_firework"]], "fw_generator_dicts": [self["opt_obj_dict"],[self["vib_obj_dict"],irc_obj_dict_forward,irc_obj_dict_reverse]], "out_names": ["opt.xyz",["vib.json","irc_forward.traj","irc_reverse.traj"]],"future_check_symms": [True,False], "label": "TS"+str(index)+"_"+rxn_name}) - cfw = Firework([ctask],name="TS"+str(index)+"_"+rxn_name+"_collect",spec={"_priority": 5}) + cfw = Firework([ctask],name="TS"+str(index)+"_"+rxn_name+"_collect",spec={"_allow_fizzled_parents": True, "_priority": 5}) newwf = Workflow([cfw],name='rxn_'+str(index)+str(rxn_name)) return FWAction(detours=newwf) #using detour allows us to inherit children from the original collect to the subsequent collects else: return FWAction() -def collect_firework(xyzs,check_symm,fw_generators,fw_generator_dicts,out_names,future_check_symms,parents=[],label=""): +def collect_firework(xyzs,check_symm,fw_generators,fw_generator_dicts,out_names,future_check_symms,parents=[],label="",allow_fizzled_parents=False): task = MolecularCollect({"xyzs": xyzs, "check_symm": check_symm, "fw_generators": fw_generators, "fw_generator_dicts": fw_generator_dicts, "out_names": out_names, "future_check_symms": future_check_symms, "label": label}) - return Firework([task],parents=parents,name=label+"collect",spec={"_priority": 5}) + return Firework([task],parents=parents,name=label+"collect",spec={"_allow_fizzled_parents": allow_fizzled_parents,"_priority": 5}) @explicit_serialize class MolecularCollect(CollectTask): @@ -667,7 +667,7 @@ def run_task(self, fw_spec): task = MolecularCollect({"xyzs": out_xyzs,"check_symm": future_check_symms[0], "fw_generators": fw_generators[1:],"fw_generator_dicts": fw_generator_dicts[1:], "out_names": out_names[1:],"future_check_symms": future_check_symms[1:],"label": self["label"]}) - cfw = Firework([task],parents=fws,name=self["label"]+"collect") + cfw = Firework([task],parents=fws,name=self["label"]+"collect",spec={"_allow_fizzled_parents":True,"_priority": 4}) newwf = Workflow(fws+[cfw],name=self["label"]+"collect"+str(-len(self["fw_generators"]))) return FWAction(detours=newwf) #using detour allows us to inherit children from the original collect to the subsequent collects else: From ec5397545019a318a408be6371c30b28399ab69a Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 7 Mar 2023 16:34:08 -0800 Subject: [PATCH 26/50] allow fizzled parents for vib collect and TSest --- pynta/main.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pynta/main.py b/pynta/main.py index 5243c724..e2ace25d 100644 --- a/pynta/main.py +++ b/pynta/main.py @@ -422,7 +422,7 @@ def setup_adsorbates(self,initial_guess_finished=False): vib_obj_dict = {"software": self.software, "label": ad, "software_kwargs": software_kwargs, "constraints": ["freeze up to "+str(self.nslab)]} - cfw = collect_firework(xyzs,True,["vibrations_firework"],[vib_obj_dict],["vib.json"],[True,False],parents=optfws2,label=ad) + cfw = collect_firework(xyzs,True,["vibrations_firework"],[vib_obj_dict],["vib.json"],[True,False],parents=optfws2,label=ad,allow_fizzled_parents=False) self.adsorbate_fw_dict[ad] = optfws2 logging.error(self.adsorbate_fw_dict.keys()) self.fws.extend(optfws+[cfw]) @@ -460,7 +460,7 @@ def setup_transition_states(self,adsorbates_finished=False): if not adsorbates_finished: for m in reactants+products: parents.extend(self.adsorbate_fw_dict[m]) - fw = Firework([ts_task],parents=parents,name="TS"+str(i)+"est",spec={"_priority": 10}) + fw = Firework([ts_task],parents=parents,name="TS"+str(i)+"est",spec={"_allow_fizzled_parents": True,"_priority": 10}) self.fws.append(fw) def launch(self,single_job=False): From b2abd0da90ebababe597936661f59e2121d9f718 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 7 Mar 2023 17:56:12 -0800 Subject: [PATCH 27/50] fix bug in generating the template_to_ase map --- pynta/tasks.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pynta/tasks.py b/pynta/tasks.py index 94e45d97..6372945a 100644 --- a/pynta/tasks.py +++ b/pynta/tasks.py @@ -541,7 +541,9 @@ def run_task(self, fw_spec): with open(os.path.join(ts_path,"info.json"),'w') as f: json.dump(ts_dict,f) - template_to_ase = {i:get_ase_index(i,template_mol_map,gratom_to_molecule_atom_maps,nslab,[ads_size(mol) for mol in mols]) for i in range(len(reactants.atoms))} + molecule_to_atom_maps = [{value:key for key,value in gratom_to_molecule_atom_maps[name].items()} for name in species_names] + template_to_ase = {i:get_ase_index(i,template_mol_map,molecule_to_atom_maps, + nslab,[ads_size(mol) for mol in mols]) for i in range(len(reactants.atoms))} ase_to_mol_num = {} for tind,aind in template_to_ase.items(): if aind: From 978d210ea185044aa97aa40427382c6826882e4d Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 7 Mar 2023 18:23:03 -0800 Subject: [PATCH 28/50] allow failed xtb runs --- pynta/tasks.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pynta/tasks.py b/pynta/tasks.py index 6372945a..423d4a59 100644 --- a/pynta/tasks.py +++ b/pynta/tasks.py @@ -847,7 +847,9 @@ def map_harmonically_forced_xtb(input): json.dump(d,f) write(os.path.join(ts_path,str(j),"xtb.xyz"),sp) xyz = os.path.join(ts_path,str(j),"xtb.xyz") - return (sp,Eharm,xyz) + return (sp,Eharm,xyz) + else: + return (None,None,None) class StructureError(Exception): pass class TimeLimitError(Exception): pass From 202b764efc8fa94736fb4a93fb75a712758aa7f9 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 8 Mar 2023 09:09:33 -0800 Subject: [PATCH 29/50] fix bug placing gas phase species --- pynta/transitionstate.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pynta/transitionstate.py b/pynta/transitionstate.py index 23d11be7..ada9a91f 100644 --- a/pynta/transitionstate.py +++ b/pynta/transitionstate.py @@ -226,7 +226,7 @@ def get_unique_TS_structs(adsorbates,species_names,cas,nslab,num_surf_sites,mol_ c += 1 site = sites[c] - add_adsorbate_to_site(adslab,adsorbate=adss[1],site=site,height=gas_height) + add_adsorbate_to_site(adslab,adsorbate=adss[1],surf_ind=0,site=site,height=gas_height) if len(adss) == 2: tsstructs.append(adslab) else: @@ -235,7 +235,7 @@ def get_unique_TS_structs(adsorbates,species_names,cas,nslab,num_surf_sites,mol_ while site2["occupied"] == True and site2 != site: c += 1 site2 = sites[c] - add_adsorbate_to_site(adslab,adsorbate=adss[2],site=site) + add_adsorbate_to_site(adslab,adsorbate=adss[2],surf_ind=0,site=site2,height=gas_height) if len(adss) == 3: tsstructs.append(adslab) else: From 5d83eca6aaf7d24dac760c9fed223c8c0764037a Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 8 Mar 2023 09:19:26 -0800 Subject: [PATCH 30/50] handle case where the first species is gas phase --- pynta/tasks.py | 2 +- pynta/transitionstate.py | 11 ++++++++--- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/pynta/tasks.py b/pynta/tasks.py index 423d4a59..0310d59f 100644 --- a/pynta/tasks.py +++ b/pynta/tasks.py @@ -552,7 +552,7 @@ def run_task(self, fw_spec): ase_to_mol_num[aind] = i break - tsstructs = get_unique_TS_structs(adsorbates,species_names,cas,nslab,num_surf_sites,mol_dict, + tsstructs = get_unique_TS_structs(adsorbates,species_names,slab,cas,nslab,num_surf_sites,mol_dict, gratom_to_molecule_atom_maps,gratom_to_molecule_surface_atom_maps, facet,metal) diff --git a/pynta/transitionstate.py b/pynta/transitionstate.py index ada9a91f..887d3079 100644 --- a/pynta/transitionstate.py +++ b/pynta/transitionstate.py @@ -159,7 +159,7 @@ def determine_TS_construction(reactant_names,reactant_mols,product_names,product return forward,ordered_reacting_species -def get_unique_TS_structs(adsorbates,species_names,cas,nslab,num_surf_sites,mol_dict, +def get_unique_TS_structs(adsorbates,species_names,slab,cas,nslab,num_surf_sites,mol_dict, gratom_to_molecule_atom_maps,gratom_to_molecule_surface_atom_maps, facet,metal,gas_height=5.0): """ @@ -168,9 +168,14 @@ def get_unique_TS_structs(adsorbates,species_names,cas,nslab,num_surf_sites,mol_ tsstructs = [] ordered_adsorbates = [adsorbates[name] for name in species_names] for adss in itertools.product(*ordered_adsorbates): - adslab = adss[0] if len(adss) == 1: - tsstructs.append(adslab) + if num_surf_sites[0] > 0: + tsstructs.append(adss[0]) + else: + adslab = deepcopy(slab) + site = cas.get_sites()[0] + add_adsorbate_to_site(adslab,adsorbate=adss[0],surf_ind=0,site=site,height=gas_height) + tsstructs.append(adslab) else: if num_surf_sites[1] == 1: name = species_names[1] From 8edb0c0e762c109b5ff80511c8bcbff75069e922 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 8 Mar 2023 18:31:53 -0800 Subject: [PATCH 31/50] ensure surf_ind specified --- pynta/transitionstate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pynta/transitionstate.py b/pynta/transitionstate.py index 887d3079..3b0c22d0 100644 --- a/pynta/transitionstate.py +++ b/pynta/transitionstate.py @@ -220,7 +220,7 @@ def get_unique_TS_structs(adsorbates,species_names,slab,cas,nslab,num_surf_sites while site["occupied"] == True: c += 1 site = sites[c] - add_adsorbate_to_site(adslab,adsorbate=adss[2],site=site,height=gas_height) + add_adsorbate_to_site(adslab,adsorbate=adss[2],surf_ind=0,site=site,height=gas_height) elif num_surf_sites[1] == 0: adcovl1 = SlabAdsorbateCoverage(adslab,adsorption_sites=cas) From f93762a00e326956ed9616dad387c15d2d998f5c Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Fri, 10 Mar 2023 09:59:37 -0800 Subject: [PATCH 32/50] fix bug in adslab definition --- pynta/transitionstate.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/pynta/transitionstate.py b/pynta/transitionstate.py index 3b0c22d0..8b095dee 100644 --- a/pynta/transitionstate.py +++ b/pynta/transitionstate.py @@ -168,14 +168,14 @@ def get_unique_TS_structs(adsorbates,species_names,slab,cas,nslab,num_surf_sites tsstructs = [] ordered_adsorbates = [adsorbates[name] for name in species_names] for adss in itertools.product(*ordered_adsorbates): + if num_surf_sites[0] > 0: + adslab = adss[0].copy() + else: + adslab = slab.copy() + site = cas.get_sites()[0] + add_adsorbate_to_site(adslab,adsorbate=adss[0],surf_ind=0,site=site,height=gas_height) if len(adss) == 1: - if num_surf_sites[0] > 0: - tsstructs.append(adss[0]) - else: - adslab = deepcopy(slab) - site = cas.get_sites()[0] - add_adsorbate_to_site(adslab,adsorbate=adss[0],surf_ind=0,site=site,height=gas_height) - tsstructs.append(adslab) + tsstructs.append(adslab) else: if num_surf_sites[1] == 1: name = species_names[1] From b567fb8495861223e7e67692a950d11b22053127 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Fri, 10 Mar 2023 16:55:56 -0800 Subject: [PATCH 33/50] handle gas phase species in non-pbc fallback --- pynta/calculator.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/pynta/calculator.py b/pynta/calculator.py index 65af5a12..634c2e4c 100644 --- a/pynta/calculator.py +++ b/pynta/calculator.py @@ -162,6 +162,14 @@ def run_harmonically_forced_xtb_no_pbc(atoms,atom_bond_potentials,site_bond_pote site_poss[i] = site_pos site_inds[i] = ind + for mol_ind in set(ase_to_mol_num.values()): #handle gas phase species that don't have site_bond_potentials + if mol_ind not in site_mols: + for key,val in ase_to_mol_num.items(): + if val == mol_ind: + site_mols.append(mol_ind) + site_poss.append(atoms.positions[key]) + break + trans = get_best_translation(site_poss,atoms.cell) mol_to_trans = {site_mols[i]: trans[i] for i in range(len(site_mols))} From 7829d3365cf9066245ac88810b424e381b70c305 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 14 Mar 2023 17:22:56 -0700 Subject: [PATCH 34/50] adjust get_best_translation before the translation was chosen first to keep the sites close to each other and then to keep them close to the center now for the second criteria instead of looking at the sites we look at the distances of the vertices of a box surrounding the adsorbate in xy fro the center --- pynta/calculator.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/pynta/calculator.py b/pynta/calculator.py index 634c2e4c..a9937eb4 100644 --- a/pynta/calculator.py +++ b/pynta/calculator.py @@ -170,7 +170,14 @@ def run_harmonically_forced_xtb_no_pbc(atoms,atom_bond_potentials,site_bond_pote site_poss.append(atoms.positions[key]) break - trans = get_best_translation(site_poss,atoms.cell) + aposx = [a.position[0] for a in atoms[nslab:]] + aposy = [a.position[1] for a in atoms[nslab:]] + apos = [np.array([min(aposx),min(aposy)]), + np.array([min(aposx),max(aposy)]), + np.array([max(aposx),min(aposy)]), + np.array([max(aposx),max(aposy)])] + trans = get_best_translation(site_poss,apos,atoms.cell) + mol_to_trans = {site_mols[i]: trans[i] for i in range(len(site_mols))} new_site_potentials = [] @@ -264,7 +271,7 @@ def run_harmonically_forced_xtb_no_pbc(atoms,atom_bond_potentials,site_bond_pote return outadslab,Eharm,Fharm -def get_best_translation(poss,cell): +def get_best_translation(poss,apos,cell): target = cell[0][:2] + cell[1][:2] pos2ds = [np.array(pos[:2]) for pos in poss] translations = [np.zeros(2),cell[0][:2],cell[1][:2],cell[0][:2] + cell[1][:2]] @@ -275,7 +282,8 @@ def get_best_translation(poss,cell): dist = 0.0 pos2ddist = 0.0 for i in range(len(pos2ds)): - dist += np.linalg.norm(pos2ds[i]+transs[i]-target) + for pos in apos: + dist += np.linalg.norm(pos+transs-target) for j in range(i): pos2ddist += np.linalg.norm(pos2ds[i]+transs[i]-(pos2ds[j]+transs[j])) if abs(pos2ddist - minpos2ddist) > 0.1 and pos2ddist < minpos2ddist: From 7f82c4d63fe3c0e8d1659efe939b1388a222f38f Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 15 Mar 2023 17:04:35 -0700 Subject: [PATCH 35/50] searching a 3x3 extension for non-pbc placement instead of a 2x2 extension this greatly improves our ability to keep the adsorbates away from the edges --- pynta/calculator.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pynta/calculator.py b/pynta/calculator.py index a9937eb4..da04b189 100644 --- a/pynta/calculator.py +++ b/pynta/calculator.py @@ -190,7 +190,7 @@ def run_harmonically_forced_xtb_no_pbc(atoms,atom_bond_potentials,site_bond_pote slab = atoms[:nslab] ad = atoms[nslab:] - bigslab = slab * (2,2,1) + bigslab = slab * (3,3,1) for ind,molind in ase_to_mol_num.items(): ad.positions[ind-nslab] += mol_to_trans[molind] @@ -272,9 +272,10 @@ def run_harmonically_forced_xtb_no_pbc(atoms,atom_bond_potentials,site_bond_pote return outadslab,Eharm,Fharm def get_best_translation(poss,apos,cell): - target = cell[0][:2] + cell[1][:2] + target = (cell[0][:2] + cell[1][:2])*1.5 pos2ds = [np.array(pos[:2]) for pos in poss] - translations = [np.zeros(2),cell[0][:2],cell[1][:2],cell[0][:2] + cell[1][:2]] + translations = [np.zeros(2),cell[0][:2],cell[1][:2],cell[0][:2] + cell[1][:2], + 2*cell[0][:2], 2*cell[1][:2], 2*cell[0][:2] + cell[1][:2],cell[0][:2] + 2*cell[1][:2],2*cell[0][:2] + 2*cell[1][:2]] mindist = np.inf minpos2ddist = np.inf tranout = None From c66b229e27c1d145010381d4eb172eda54a0d51d Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 15 Mar 2023 17:05:23 -0700 Subject: [PATCH 36/50] account for site position in site selection and atom removal --- pynta/calculator.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/pynta/calculator.py b/pynta/calculator.py index da04b189..2515cc1a 100644 --- a/pynta/calculator.py +++ b/pynta/calculator.py @@ -170,8 +170,8 @@ def run_harmonically_forced_xtb_no_pbc(atoms,atom_bond_potentials,site_bond_pote site_poss.append(atoms.positions[key]) break - aposx = [a.position[0] for a in atoms[nslab:]] - aposy = [a.position[1] for a in atoms[nslab:]] + aposx = [a.position[0] for a in atoms[nslab:]] + [a[0] for a in site_poss] + aposy = [a.position[1] for a in atoms[nslab:]] + [a[1] for a in site_poss] apos = [np.array([min(aposx),min(aposy)]), np.array([min(aposx),max(aposy)]), np.array([max(aposx),min(aposy)]), @@ -204,6 +204,9 @@ def run_harmonically_forced_xtb_no_pbc(atoms,atom_bond_potentials,site_bond_pote for j,a in enumerate(ad): if np.linalg.norm(atom.position[:2]-a.position[:2]) < dist: dist = np.linalg.norm(atom.position[:2]-a.position[:2]) + for sitep in new_site_potentials: + if np.linalg.norm(atom.position[:2]-sitep["site_pos"][:2]) < dist: + dist = np.linalg.norm(atom.position[:2]-sitep["site_pos"][:2]) if dist > dthresh: delinds.append(i) From 95dcb15e890305339fe6e53f8b341d213fb372bd Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Mon, 15 May 2023 15:59:37 -0700 Subject: [PATCH 37/50] add HFSP firework generation --- pynta/tasks.py | 60 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) diff --git a/pynta/tasks.py b/pynta/tasks.py index 0310d59f..e3969488 100644 --- a/pynta/tasks.py +++ b/pynta/tasks.py @@ -832,6 +832,66 @@ def run_task(self, fw_spec): return FWAction() +def HFSP_firework(xyz,atom_bond_potentials,site_bond_potentials,nslab,constraints,molecule_to_atom_maps,ase_to_mol_num, + out_path=None,label="",parents=[],ignore_errors=False): + d = {"xyz": xyz, "atom_bond_potentials": atom_bond_potentials, "site_bond_potentials": site_bond_potentials, + "nslab": nslab, "constraints": constraints, "molecule_to_atom_maps": molecule_to_atom_maps, "ase_to_mol_num": ase_to_mol_num, + "label": label, "ignore_errors": ignore_errors} + t1 = MolecularHFSP(d) + directory = os.path.dirname(xyz) + if out_path is None: out_path = os.path.join(directory,label+".traj") + t2 = FileTransferTask({'files': [{'src': label+'.xyz', 'dest': out_path}, {'src': "xtbharm.traj", 'dest': os.path.join(directory,label+".traj")}], + 'mode': 'copy', 'ignore_errors' : ignore_errors}) + return Firework([t1,t2],parents=parents,name=label+"HFSP") + +@explicit_serialize +class MolecularHFSP(OptimizationTask): + required_params = ["xyz","atom_bond_potentials","site_bond_potentials","nslab","constraints","molecule_to_atom_maps","ase_to_mol_num"] + optional_params = ["label","ignore_errors","method"] + + def run_task(self, fw_spec): + label = self["label"] if "label" in self.keys() else "xtb" + ignore_errors = self["ignore_errors"] if "ignore_errors" in self.keys() else False + method = self["method"] if "method" in self.keys() else "GFN1-xTB" + + atom_bond_potentials = self["atom_bond_potentials"] + site_bond_potentials = self["site_bond_potentials"] + nslab = self["nslab"] + molecule_to_atom_maps = self["molecule_to_atom_maps"] + ase_to_mol_num = self["ase_to_mol_num"] + constraints = self["constraints"] + xyz = self['xyz'] + + errors = [] + + suffix = os.path.split(xyz)[-1].split(".")[-1] + + try: + if suffix == "xyz": + sp = read(xyz) + elif suffix == "traj": #take last point on trajectory + sp = Trajectory(xyz)[-1] + else: #assume xyz + sp = read(xyz) + except Exception as e: + if not ignore_errors: + raise e + else: + errors.append(e) + + spout,Eharm,Fharm = run_harmonically_forced_xtb(sp,atom_bond_potentials,site_bond_potentials,nslab, + molecule_to_atom_maps=molecule_to_atom_maps,ase_to_mol_num=ase_to_mol_num, + method="GFN1-xTB",constraints=constraints) + if spout: + if "initial_charges" in sp.arrays.keys(): #avoid bug in ase + del sp.arrays["initial_charges"] + write(label+".xyz", spout) + converged = True + else: + converged = False + + return FWAction(stored_data={"error": errors,"converged": converged}) + def map_harmonically_forced_xtb(input): tsstruct,atom_bond_potentials,site_bond_potentials,nslab,constraints,ts_path,j,molecule_to_atom_maps,ase_to_mol_num = input os.makedirs(os.path.join(ts_path,str(j))) From 8afac2c22ba6998d7218d249e45f984555758883 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 17 May 2023 17:18:38 -0700 Subject: [PATCH 38/50] update HFSP optimizations to avoid defining constraints using Sella --- pynta/calculator.py | 55 +++++++++++++++++++++++++++++++-------------- 1 file changed, 38 insertions(+), 17 deletions(-) diff --git a/pynta/calculator.py b/pynta/calculator.py index 2515cc1a..fed8e3d0 100644 --- a/pynta/calculator.py +++ b/pynta/calculator.py @@ -94,18 +94,28 @@ def run_harmonically_forced_xtb(atoms,atom_bond_potentials,site_bond_potentials, """ Optimize TS guess using xTB + harmonic forcing terms determined by atom_bond_potentials and site_bond_potentials """ - cons = Constraints(atoms) - for c in constraints: if isinstance(c,dict): - add_sella_constraint(cons,c) + constraint = construct_constraint(c) + sp.set_constraint(constraint) + elif c == "freeze half slab": + sp.set_constraint(FixAtoms([ + atom.index for atom in sp if atom.position[2] < sp.cell[2, 2] / 2. + ])) elif c == "freeze slab": - for i,atom in enumerate(atoms): #freeze the slab - if i < nslab: - cons.fix_translation(atom.index) - else: - raise ValueError("Constraint {} not understood".format(c)) - + sp.set_constraint(FixAtoms( + indices=list(range(nslab)) + )) + elif c.split()[0] == "freeze" and c.split()[1] == "all": #ex: "freeze all Cu" + sym = c.split()[2] + sp.set_constraint(FixAtoms( + indices=[atom.index for atom in sp if atom.symbol == sym] + )) + elif c.split()[0] == "freeze" and c.split()[1] == "up" and c.split()[2] == "to": + n = int(c.split()[3]) + sp.set_constraint(FixAtoms( + indices=list(range(n)) + )) hfxtb = HarmonicallyForcedXTB(method="GFN1-xTB", atom_bond_potentials=atom_bond_potentials, @@ -237,17 +247,28 @@ def run_harmonically_forced_xtb_no_pbc(atoms,atom_bond_potentials,site_bond_pote c["indices"][i] += new_nslab-nslab new_constraints.append(c) - cons = Constraints(bigad) for c in new_constraints: if isinstance(c,dict): - add_sella_constraint(cons,c) + constraint = construct_constraint(c) + sp.set_constraint(constraint) + elif c == "freeze half slab": + sp.set_constraint(FixAtoms([ + atom.index for atom in sp if atom.position[2] < sp.cell[2, 2] / 2. + ])) elif c == "freeze slab": - for i,atom in enumerate(bigad): #freeze the slab - if i < new_nslab: - cons.fix_translation(atom.index) - else: - raise ValueError("Constraint {} not understood".format(c)) - + sp.set_constraint(FixAtoms( + indices=list(range(nslab)) + )) + elif c.split()[0] == "freeze" and c.split()[1] == "all": #ex: "freeze all Cu" + sym = c.split()[2] + sp.set_constraint(FixAtoms( + indices=[atom.index for atom in sp if atom.symbol == sym] + )) + elif c.split()[0] == "freeze" and c.split()[1] == "up" and c.split()[2] == "to": + n = int(c.split()[3]) + sp.set_constraint(FixAtoms( + indices=list(range(n)) + )) hfxtb = HarmonicallyForcedXTB(method="GFN1-xTB", atom_bond_potentials=new_atom_bond_potentials, From 3d734747e0fc92d6b8abec14c8e4bb15bb293c99 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 17 May 2023 17:30:35 -0700 Subject: [PATCH 39/50] unwrap dictionaries properly in HFSP firework --- pynta/tasks.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pynta/tasks.py b/pynta/tasks.py index e3969488..48455387 100644 --- a/pynta/tasks.py +++ b/pynta/tasks.py @@ -857,8 +857,8 @@ def run_task(self, fw_spec): atom_bond_potentials = self["atom_bond_potentials"] site_bond_potentials = self["site_bond_potentials"] nslab = self["nslab"] - molecule_to_atom_maps = self["molecule_to_atom_maps"] - ase_to_mol_num = self["ase_to_mol_num"] + molecule_to_atom_maps = {int(k):v for k,v in self["molecule_to_atom_maps"].items()} + ase_to_mol_num = {int(k):v for k,v in self["ase_to_mol_num"].items()} constraints = self["constraints"] xyz = self['xyz'] From 3425133c9b2a11817283b702780063f437a7a3a0 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 17 May 2023 17:32:55 -0700 Subject: [PATCH 40/50] fix unwrapping in HFSP firework --- pynta/tasks.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pynta/tasks.py b/pynta/tasks.py index 48455387..b21e39e9 100644 --- a/pynta/tasks.py +++ b/pynta/tasks.py @@ -857,7 +857,7 @@ def run_task(self, fw_spec): atom_bond_potentials = self["atom_bond_potentials"] site_bond_potentials = self["site_bond_potentials"] nslab = self["nslab"] - molecule_to_atom_maps = {int(k):v for k,v in self["molecule_to_atom_maps"].items()} + molecule_to_atom_maps = [{int(k):v for k,v in x.items()} for x in self["molecule_to_atom_maps"]] ase_to_mol_num = {int(k):v for k,v in self["ase_to_mol_num"].items()} constraints = self["constraints"] xyz = self['xyz'] From 3cbba1adb8114ad6f334b6620870868a613c6ae8 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 17 May 2023 23:36:44 -0700 Subject: [PATCH 41/50] fix bug in HFSP constraints --- pynta/calculator.py | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/pynta/calculator.py b/pynta/calculator.py index fed8e3d0..082cf850 100644 --- a/pynta/calculator.py +++ b/pynta/calculator.py @@ -97,23 +97,23 @@ def run_harmonically_forced_xtb(atoms,atom_bond_potentials,site_bond_potentials, for c in constraints: if isinstance(c,dict): constraint = construct_constraint(c) - sp.set_constraint(constraint) + atoms.set_constraint(constraint) elif c == "freeze half slab": - sp.set_constraint(FixAtoms([ - atom.index for atom in sp if atom.position[2] < sp.cell[2, 2] / 2. + atoms.set_constraint(FixAtoms([ + atom.index for atom in atoms if atom.position[2] < atoms.cell[2, 2] / 2. ])) elif c == "freeze slab": - sp.set_constraint(FixAtoms( + atoms.set_constraint(FixAtoms( indices=list(range(nslab)) )) elif c.split()[0] == "freeze" and c.split()[1] == "all": #ex: "freeze all Cu" sym = c.split()[2] - sp.set_constraint(FixAtoms( - indices=[atom.index for atom in sp if atom.symbol == sym] + atoms.set_constraint(FixAtoms( + indices=[atom.index for atom in atoms if atom.symbol == sym] )) elif c.split()[0] == "freeze" and c.split()[1] == "up" and c.split()[2] == "to": n = int(c.split()[3]) - sp.set_constraint(FixAtoms( + atoms.set_constraint(FixAtoms( indices=list(range(n)) )) @@ -250,23 +250,23 @@ def run_harmonically_forced_xtb_no_pbc(atoms,atom_bond_potentials,site_bond_pote for c in new_constraints: if isinstance(c,dict): constraint = construct_constraint(c) - sp.set_constraint(constraint) + atoms.set_constraint(constraint) elif c == "freeze half slab": - sp.set_constraint(FixAtoms([ - atom.index for atom in sp if atom.position[2] < sp.cell[2, 2] / 2. + bigad.set_constraint(FixAtoms([ + atom.index for atom in bigad if atom.position[2] < bigad.cell[2, 2] / 2. ])) elif c == "freeze slab": - sp.set_constraint(FixAtoms( - indices=list(range(nslab)) + bigad.set_constraint(FixAtoms( + indices=list(range(new_nslab)) )) elif c.split()[0] == "freeze" and c.split()[1] == "all": #ex: "freeze all Cu" sym = c.split()[2] - sp.set_constraint(FixAtoms( - indices=[atom.index for atom in sp if atom.symbol == sym] + bigad.set_constraint(FixAtoms( + indices=[atom.index for atom in bigad if atom.symbol == sym] )) elif c.split()[0] == "freeze" and c.split()[1] == "up" and c.split()[2] == "to": n = int(c.split()[3]) - sp.set_constraint(FixAtoms( + atoms.set_constraint(FixAtoms( indices=list(range(n)) )) From 3e6df604473fafd1d8f2fe093e173058c5009ad6 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Thu, 18 May 2023 09:28:15 -0700 Subject: [PATCH 42/50] import constraints in calculator.py --- pynta/calculator.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pynta/calculator.py b/pynta/calculator.py index 082cf850..7ad7a98d 100644 --- a/pynta/calculator.py +++ b/pynta/calculator.py @@ -7,6 +7,7 @@ from ase.calculators import calculator from ase.data import reference_states,chemical_symbols from ase.build import bulk +from ase.constraints import * from pynta.utils import name_to_ase_software from sella import Sella, Constraints import scipy.optimize as opt From fa8391084e80a8e0156863a1a1de3fe66de269b9 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Thu, 18 May 2023 10:38:27 -0700 Subject: [PATCH 43/50] move construct_constraint to utils --- pynta/calculator.py | 1 + pynta/tasks.py | 11 ----------- pynta/utils.py | 14 ++++++++++++++ 3 files changed, 15 insertions(+), 11 deletions(-) diff --git a/pynta/calculator.py b/pynta/calculator.py index 7ad7a98d..f2273d75 100644 --- a/pynta/calculator.py +++ b/pynta/calculator.py @@ -14,6 +14,7 @@ import copy from copy import deepcopy import itertools +from pynta.utils import * def get_energy_atom_bond(atoms,ind1,ind2,k,deq): bd,d = get_distances([atoms.positions[ind1]], [atoms.positions[ind2]], cell=atoms.cell, pbc=atoms.pbc) diff --git a/pynta/tasks.py b/pynta/tasks.py index b21e39e9..f5ccf3b7 100644 --- a/pynta/tasks.py +++ b/pynta/tasks.py @@ -994,17 +994,6 @@ def reconstruct_firework(new_task,old_task,task_list,full=True): tasks.append(reconstruct_task(d)) return Firework(tasks) -def construct_constraint(d): - """ - construct a constrain from a dictionary that is the input to the constraint - constructor plus an additional "type" key that indices the name of the constraint - returns the constraint - """ - constraint_dict = copy.deepcopy(d) - constructor = getattr("ase.constraints",constraint_dict["type"]) - del constraint_dict["type"] - return constructor(**constraint_dict) - def get_fizzled_fws(lpad): return [lpad.get_fw_by_id(i) for i in lpad.get_fw_ids_in_wfs() if lpad.get_fw_by_id(i).state == "FIZZLED"] diff --git a/pynta/utils.py b/pynta/utils.py index 065471ea..e901e330 100644 --- a/pynta/utils.py +++ b/pynta/utils.py @@ -3,9 +3,12 @@ import ase from ase.utils.structure_comparator import SymmetryEquivalenceCheck from ase.io import write, read +import ase.constraints from copy import deepcopy from importlib import import_module import numpy as np +import copy + def get_unique_sym(geoms): ''' Check for the symmetry equivalent structures in the given files @@ -236,3 +239,14 @@ def clean_pynta_path(path,save_initial_guess=True): os.remove(pa) else: os.remove(os.path.join(path,p,ad,ind)) + +def construct_constraint(d): + """ + construct a constrain from a dictionary that is the input to the constraint + constructor plus an additional "type" key that indices the name of the constraint + returns the constraint + """ + constraint_dict = copy.deepcopy(d) + constructor = getattr(ase.constraints,constraint_dict["type"]) + del constraint_dict["type"] + return constructor(**constraint_dict) \ No newline at end of file From c425aba0d91a2b3d03ba80b85b3af6ace54d5bb1 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Thu, 18 May 2023 10:38:51 -0700 Subject: [PATCH 44/50] remove constraints input from Sella --- pynta/calculator.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pynta/calculator.py b/pynta/calculator.py index f2273d75..70dcfa28 100644 --- a/pynta/calculator.py +++ b/pynta/calculator.py @@ -125,7 +125,7 @@ def run_harmonically_forced_xtb(atoms,atom_bond_potentials,site_bond_potentials, atoms.calc = hfxtb - opt = Sella(atoms,constraints=cons,trajectory="xtbharm.traj",order=0) + opt = Sella(atoms,trajectory="xtbharm.traj",order=0) try: opt.run(fmax=0.02,steps=150) @@ -278,7 +278,7 @@ def run_harmonically_forced_xtb_no_pbc(atoms,atom_bond_potentials,site_bond_pote bigad.calc = hfxtb - opt = Sella(bigad,constraints=cons,trajectory="xtbharm.traj",order=0) + opt = Sella(bigad,trajectory="xtbharm.traj",order=0) try: opt.run(fmax=0.02,steps=150) From 977bb8bd74054c92acb439b125173face4d525d9 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Thu, 18 May 2023 10:39:28 -0700 Subject: [PATCH 45/50] switch remaining xtb constraints handling from Sella to ASE --- pynta/calculator.py | 4 ++++ pynta/mol.py | 2 +- pynta/transitionstate.py | 2 +- 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/pynta/calculator.py b/pynta/calculator.py index 70dcfa28..11cf82c8 100644 --- a/pynta/calculator.py +++ b/pynta/calculator.py @@ -247,6 +247,10 @@ def run_harmonically_forced_xtb_no_pbc(atoms,atom_bond_potentials,site_bond_pote c = deepcopy(constraint) for i in range(len(c["indices"])): c["indices"][i] += new_nslab-nslab + if "a1" in c: + c["a1"] += new_nslab-nslab + if "a2" in c: + c["a2"] += new_nslab-nslab new_constraints.append(c) for c in new_constraints: diff --git a/pynta/mol.py b/pynta/mol.py index b80d2bc1..5461b9d2 100644 --- a/pynta/mol.py +++ b/pynta/mol.py @@ -97,7 +97,7 @@ def generate_adsorbate_guesses(mol,ads,full_slab,cas,mol_to_atoms_map,metal, mol_fixed_bond_pairs = [[mol.atoms.index(bd.atom1),mol.atoms.index(bd.atom2)] for bd in mol.get_all_edges() if (not bd.atom1.is_surface_site()) and (not bd.atom2.is_surface_site())] atom_fixed_bond_pairs = [[mol_to_atoms_map[pair[0]]+len(full_slab),mol_to_atoms_map[pair[1]]+len(full_slab)]for pair in mol_fixed_bond_pairs] - constraint_list = [{"type": "fix_bond", "indices": pair} for pair in atom_fixed_bond_pairs]+["freeze slab"] + constraint_list = [{"type": "FixBondLength", "a1": pair[0], "a2": pair[1]} for pair in atom_fixed_bond_pairs]+["freeze slab"] geos = [] for i,sites_list in enumerate(sites_lists): diff --git a/pynta/transitionstate.py b/pynta/transitionstate.py index 8b095dee..044b05cd 100644 --- a/pynta/transitionstate.py +++ b/pynta/transitionstate.py @@ -519,7 +519,7 @@ def generate_constraints_harmonic_parameters(tsstructs,adsorbates,slab,forward_t if tsstruct_valid: if fixed_bond_pairs: - constraint_list = [{"type": "fix_bond", "indices": pair} for pair in fixed_bond_pairs]+["freeze slab"] + constraint_list = [{"type": "FixBondLength", "a1": pair[0], "a2": pair[1]} for pair in fixed_bond_pairs]+["freeze slab"] constraint_lists.append(constraint_list) else: constraint_lists.append(["freeze slab"]) From 5d18c2c98d13db6cad4a2cccbbb5d712262b4e64 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sun, 21 May 2023 10:58:41 -0700 Subject: [PATCH 46/50] fix output directory --- pynta/tasks.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pynta/tasks.py b/pynta/tasks.py index f5ccf3b7..2061235a 100644 --- a/pynta/tasks.py +++ b/pynta/tasks.py @@ -839,7 +839,7 @@ def HFSP_firework(xyz,atom_bond_potentials,site_bond_potentials,nslab,constraint "label": label, "ignore_errors": ignore_errors} t1 = MolecularHFSP(d) directory = os.path.dirname(xyz) - if out_path is None: out_path = os.path.join(directory,label+".traj") + if out_path is None: out_path = os.path.join(directory,label+".xyz") t2 = FileTransferTask({'files': [{'src': label+'.xyz', 'dest': out_path}, {'src': "xtbharm.traj", 'dest': os.path.join(directory,label+".traj")}], 'mode': 'copy', 'ignore_errors' : ignore_errors}) return Firework([t1,t2],parents=parents,name=label+"HFSP") From 663d9f922717363c5b3dee6edcd578f2f1a96498 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 23 May 2023 20:41:32 -0700 Subject: [PATCH 47/50] fix handling of multiple constraints set_constraint clears prior constraints so it is important to only call it once with all desired constraints --- pynta/calculator.py | 31 ++++++++++++++++------------ pynta/tasks.py | 49 +++++++++++++++++++++++++++------------------ 2 files changed, 48 insertions(+), 32 deletions(-) diff --git a/pynta/calculator.py b/pynta/calculator.py index 11cf82c8..73ed1416 100644 --- a/pynta/calculator.py +++ b/pynta/calculator.py @@ -96,29 +96,32 @@ def run_harmonically_forced_xtb(atoms,atom_bond_potentials,site_bond_potentials, """ Optimize TS guess using xTB + harmonic forcing terms determined by atom_bond_potentials and site_bond_potentials """ + out_constraints = [] for c in constraints: if isinstance(c,dict): constraint = construct_constraint(c) - atoms.set_constraint(constraint) + out_constraints.append(constraint) elif c == "freeze half slab": - atoms.set_constraint(FixAtoms([ + out_constraints.append(FixAtoms([ atom.index for atom in atoms if atom.position[2] < atoms.cell[2, 2] / 2. ])) elif c == "freeze slab": - atoms.set_constraint(FixAtoms( + out_constraints.append(FixAtoms( indices=list(range(nslab)) )) elif c.split()[0] == "freeze" and c.split()[1] == "all": #ex: "freeze all Cu" sym = c.split()[2] - atoms.set_constraint(FixAtoms( + out_constraints.append(FixAtoms( indices=[atom.index for atom in atoms if atom.symbol == sym] )) elif c.split()[0] == "freeze" and c.split()[1] == "up" and c.split()[2] == "to": n = int(c.split()[3]) - atoms.set_constraint(FixAtoms( + out_constraints.append(FixAtoms( indices=list(range(n)) )) - + + atoms.set_constraint(out_constraints) + hfxtb = HarmonicallyForcedXTB(method="GFN1-xTB", atom_bond_potentials=atom_bond_potentials, site_bond_potentials=site_bond_potentials) @@ -252,30 +255,32 @@ def run_harmonically_forced_xtb_no_pbc(atoms,atom_bond_potentials,site_bond_pote if "a2" in c: c["a2"] += new_nslab-nslab new_constraints.append(c) - + + out_constraints = [] for c in new_constraints: if isinstance(c,dict): constraint = construct_constraint(c) - atoms.set_constraint(constraint) + out_constraints.append(constraint) elif c == "freeze half slab": - bigad.set_constraint(FixAtoms([ + out_constraints.append(FixAtoms([ atom.index for atom in bigad if atom.position[2] < bigad.cell[2, 2] / 2. ])) elif c == "freeze slab": - bigad.set_constraint(FixAtoms( + out_constraints.append(FixAtoms( indices=list(range(new_nslab)) )) elif c.split()[0] == "freeze" and c.split()[1] == "all": #ex: "freeze all Cu" sym = c.split()[2] - bigad.set_constraint(FixAtoms( + out_constraints.append(FixAtoms( indices=[atom.index for atom in bigad if atom.symbol == sym] )) elif c.split()[0] == "freeze" and c.split()[1] == "up" and c.split()[2] == "to": n = int(c.split()[3]) - atoms.set_constraint(FixAtoms( + out_constraints.append(FixAtoms( indices=list(range(n)) )) - + atoms.set_constraint(out_constraints) + hfxtb = HarmonicallyForcedXTB(method="GFN1-xTB", atom_bond_potentials=new_atom_bond_potentials, site_bond_potentials=new_site_potentials) diff --git a/pynta/tasks.py b/pynta/tasks.py index 2061235a..9e5f3172 100644 --- a/pynta/tasks.py +++ b/pynta/tasks.py @@ -138,25 +138,27 @@ def run_task(self, fw_spec): if not sella: assert order == 0 + out_constraints = [] for c in constraints: if isinstance(c,dict): constraint = construct_constraint(c) - sp.set_constraint(constraint) + out_constraints.append(constraint) elif c == "freeze half slab": - sp.set_constraint(FixAtoms([ + out_constraints.append(FixAtoms([ atom.index for atom in sp if atom.position[2] < sp.cell[2, 2] / 2. ])) elif c.split()[0] == "freeze" and c.split()[1] == "all": #ex: "freeze all Cu" sym = c.split()[2] - sp.set_constraint(FixAtoms( + out_constraints.append(FixAtoms( indices=[atom.index for atom in sp if atom.symbol == sym] )) elif c.split()[0] == "freeze" and c.split()[1] == "up" and c.split()[2] == "to": n = int(c.split()[3]) - sp.set_constraint(FixAtoms( + out_constraints.append(FixAtoms( indices=list(range(n)) )) - + + sp.set_constraint(out_constraints) opt_kwargs["trajectory"] = label+".traj" opt = opt_method(sp,**opt_kwargs) @@ -186,25 +188,28 @@ def run_task(self, fw_spec): # errors.append(e) else: + out_constraints = [] for c in constraints: if isinstance(c,dict): constraint = construct_constraint(c) - sp.set_constraint(constraint) + out_constraints.append(constraint) elif c == "freeze half slab": - sp.set_constraint(FixAtoms([ + out_constraints.append(FixAtoms([ atom.index for atom in sp if atom.position[2] < sp.cell[2, 2] / 2. ])) elif c.split()[0] == "freeze" and c.split()[1] == "all": #ex: "freeze all Cu" sym = c.split()[2] - sp.set_constraint(FixAtoms( + out_constraints.append(FixAtoms( indices=[atom.index for atom in sp if atom.symbol == sym] )) elif c.split()[0] == "freeze" and c.split()[1] == "up" and c.split()[2] == "to": n = int(c.split()[3]) - sp.set_constraint(FixAtoms( + out_constraints.append(FixAtoms( indices=list(range(n)) )) + sp.set_constraint(out_constraints) + opt = Sella(sp,trajectory=label+".traj",order=order) try: if np.isinf(time_limit_hrs): @@ -400,30 +405,33 @@ def run_task(self, fw_spec): sp.calc = SocketIOCalculator(software,log=sys.stdout,unixsocket=unixsocket) if socket else software constraints = deepcopy(self["constraints"]) if "constraints" in self.keys() else [] + out_constraints = [] for c in constraints: if isinstance(c,dict): constraint = construct_constraint(c) - sp.set_constraint(constraint) + out_constraints.append(constraint) elif c == "freeze half slab": - sp.set_constraint(FixAtoms([ + out_constraints.append(FixAtoms([ atom.index for atom in sp if atom.position[2] < sp.cell[2, 2] / 2. ])) indices = [atom.index for atom in sp if atom.position[2] > sp.cell[2, 2] / 2.] elif c.split()[0] == "freeze" and c.split()[1] == "all": #ex: "freeze all Cu" sym = c.split()[2] - sp.set_constraint(FixAtoms( + out_constraints.append(FixAtoms( indices=[atom.index for atom in sp if atom.symbol == sym] )) indices = [atom.index for atom in sp if atom.symbol != sym] elif c.split()[0] == "freeze" and c.split()[1] == "up" and c.split()[2] == "to": n = int(c.split()[3]) - sp.set_constraint(FixAtoms( + out_constraints.append(FixAtoms( indices=list(range(n)) )) indices = [ i for i in range(len(sp)) if i >= n] else: raise ValueError - + + sp.set_constraint(out_constraints) + vib = Vibrations(sp,indices=indices) vib.run() @@ -773,25 +781,28 @@ def run_task(self, fw_spec): sp.calc = SocketIOCalculator(software,log=sys.stdout,unixsocket=unixsocket) if socket else software constraints = deepcopy(self["constraints"]) if "constraints" in self.keys() else [] - + + out_constraints = [] for c in constraints: if c == "freeze half slab": - sp.set_constraint(FixAtoms([ + out_constraints.append(FixAtoms([ atom.index for atom in sp if atom.position[2] < sp.cell[2, 2] / 2. ])) elif c.split()[0] == "freeze" and c.split()[1] == "all": #ex: "freeze all Cu" sym = c.split()[2] - sp.set_constraint(FixAtoms( + out_constraints.append(FixAtoms( indices=[atom.index for atom in sp if atom.symbol == sym] )) elif c.split()[0] == "freeze" and c.split()[1] == "up" and c.split()[2] == "to": n = int(c.split()[3]) - sp.set_constraint(FixAtoms( + out_constraints.append(FixAtoms( indices=list(range(n)) )) else: raise ValueError("Could not interpret constraint: {}".format(c)) - + + sp.set_constraint(out_constraints) + opt = IRC(sp,trajectory=label+"_irc.traj",dx=0.1,eta=1e-4,gamma=0.4) try: if forward: From e7a50c49842a2d75baf6b52125ed80e2510671ea Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Fri, 30 Jun 2023 11:39:49 -0700 Subject: [PATCH 48/50] skip fixing of indices if indices not present --- pynta/calculator.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pynta/calculator.py b/pynta/calculator.py index 73ed1416..38e51761 100644 --- a/pynta/calculator.py +++ b/pynta/calculator.py @@ -248,8 +248,9 @@ def run_harmonically_forced_xtb_no_pbc(atoms,atom_bond_potentials,site_bond_pote continue else: c = deepcopy(constraint) - for i in range(len(c["indices"])): - c["indices"][i] += new_nslab-nslab + if "indices" in c: + for i in range(len(c["indices"])): + c["indices"][i] += new_nslab-nslab if "a1" in c: c["a1"] += new_nslab-nslab if "a2" in c: From 60d33f4349da799f2f4e50982cbd0333c0e3b3d1 Mon Sep 17 00:00:00 2001 From: tdprice-858 <76229513+tdprice-858@users.noreply.github.com> Date: Mon, 10 Jul 2023 23:16:22 -0700 Subject: [PATCH 49/50] Update calculator.py Adding a fix to run_harmonically_forced_xtb_no_pbc function. Needed to apply the new_constraints to the bigad atom object. Test shows adsorbate is now properly constrained. --- pynta/calculator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pynta/calculator.py b/pynta/calculator.py index 38e51761..bcee94c0 100644 --- a/pynta/calculator.py +++ b/pynta/calculator.py @@ -285,7 +285,7 @@ def run_harmonically_forced_xtb_no_pbc(atoms,atom_bond_potentials,site_bond_pote hfxtb = HarmonicallyForcedXTB(method="GFN1-xTB", atom_bond_potentials=new_atom_bond_potentials, site_bond_potentials=new_site_potentials) - + bigad.set_constraint(out_constraints) bigad.calc = hfxtb opt = Sella(bigad,trajectory="xtbharm.traj",order=0) From dbcc227ffeb33a3becaabd2b09fb4dbee36db1ab Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Thu, 13 Jul 2023 14:59:38 -0700 Subject: [PATCH 50/50] delete improper constraint setting --- pynta/calculator.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pynta/calculator.py b/pynta/calculator.py index bcee94c0..81d4e351 100644 --- a/pynta/calculator.py +++ b/pynta/calculator.py @@ -280,7 +280,6 @@ def run_harmonically_forced_xtb_no_pbc(atoms,atom_bond_potentials,site_bond_pote out_constraints.append(FixAtoms( indices=list(range(n)) )) - atoms.set_constraint(out_constraints) hfxtb = HarmonicallyForcedXTB(method="GFN1-xTB", atom_bond_potentials=new_atom_bond_potentials,