#!/usr/bin/env amspython from scm.params import * """ Script to generate a training and validation set for ZnS and H2S adsorption on ZnS(110). create_parameters() illustrates how to copy and paste from different force fields in the AMS ReaxFF library to get a reasonable initial guess of parameter values. NOTE: To use the paths in import_results(), first unzip the file reference_calculation_results.zip. It contains (simplified) ams.rkf and band.rkf results files from DFT reference calculations. The settings for the reference jobs can be found in job_collection.yaml and job_collection_engines.yaml. """ def main(): create_parameters() import_results() def create_parameters(): interf = get_initial_parameters() activate_parameters(interf) # save ffield and yaml format interf.yaml_store("parameter_interface.yaml") interf.write("custom_forcefield.ff") def get_initial_parameters(bounds_scale=1.3): LiS = ReaxFFParameters("LiS.ff", bounds_scale=bounds_scale) # http://dx.doi.org/10.1039/C4CP04532G Mue2016 = ReaxFFParameters("Mue2016.ff", bounds_scale=bounds_scale) # http://dx.doi.org/10.1021/acs.jctc.6b00461 ZnOH = ReaxFFParameters("ZnOH.ff", bounds_scale=bounds_scale) # http://dx.doi.org/10.1016/j.susc.2009.12.012 # another_one = ReaxFFParameters("/absolute/path/to/some_ff.ff", bounds_scale=bounds_scale) # use custom force field not in AMS database # interf is the new parameter interface interf = ReaxFFParameters(None) # copy S parameters from LiS interf.copy_block(LiS, "ATM:S") interf.copy_block(LiS, "BND:S.S") interf.copy_block(LiS, "ANG:S.S.S") # copy H, S.H, and H.S.H parameters from Mue2016 interf.copy_block(Mue2016, "ATM:H") interf.copy_block(Mue2016, "BND:S.H") interf.copy_block(Mue2016, "ANG:H.S.H") # copy the GEN block and many other blocks from ZnOH, replacing O with S interf.copy_block(ZnOH, "GEN") interf.copy_block(ZnOH, "ATM:Zn") interf.copy_block(ZnOH, "BND:H.H", "BND:H.H") interf.copy_block(ZnOH, "BND:Zn.H", "BND:Zn.H") interf.copy_block(ZnOH, "BND:Zn.Zn", "BND:Zn.Zn") interf.copy_block(ZnOH, "BND:Zn.O", "BND:Zn.S") interf.copy_block(ZnOH, "ANG:Zn.Zn.O", "ANG:Zn.Zn.S") interf.copy_block(ZnOH, "ANG:Zn.O.Zn", "ANG:Zn.S.Zn") interf.copy_block(ZnOH, "ANG:Zn.O.O", "ANG:Zn.S.S") interf.copy_block(ZnOH, "ANG:O.Zn.O", "ANG:S.Zn.S") interf.copy_block(ZnOH, "HBD:O.H.O", "HBD:S.H.S") # give a title interf.header["head"] = "ZnS parametrization" return interf def activate_parameters(interf: ReaxFFParameters): # only optimize bond and angle parameters # do not optimize H-Zn and H-H bond parameters # do not optimize pi-bonding parameters # parameters with a value of 0 are usually best kept at 0 for p in interf: if p.metadata["category"] == "DoNotOptimize": p.is_active = False elif p.metadata["block"] == "BND": p.is_active = ( "n/a" not in p.name.lower() and not p.name.endswith("pi") and p.value != 0 and p.metadata["atoms"] != ["H", "Zn"] and p.metadata["atoms"] != ["Zn", "H"] and p.metadata["atoms"] != ["H", "H"] and not "pi bond" in p.metadata["description"].lower() ) elif p.metadata["block"] == "ANG": p.is_active = "n/a" not in p.name.lower() and p.value != 0 and p.metadata["atoms"] != ["S", "S", "S"] elif p.metadata["block"] == "HBD": p.is_active = True else: p.is_active = False def import_results(): ri = ResultsImporter(settings={"units": {"energy": "eV", "forces": "eV/angstrom"}, "remove_bonds": True}) # energy-volume scans ri.add_singlejob("dft/ZnS_pesscans/zincblende/ams.rkf", task="PESScan", properties={"pes": {"weight": 5.0}}) ri.add_singlejob("dft/ZnS_pesscans/rocksalt/ams.rkf", task="PESScan", properties={"pes": {"weight": 4.0}}) ri.add_singlejob("dft/ZnS_pesscans/wurtzite/ams.rkf", task="PESScan", properties="pes", data_set="validation_set") # single-points ri.add_singlejob("dft/ZnS_pesscans/wurtzite_sp/ams.rkf", properties="charges") ri.add_singlejob("dft/ZnS_pesscans/rocksalt_sp/ams.rkf", properties="charges") ri.add_singlejob("dft/ZnS_pesscans/zincblende_sp/ams.rkf", properties="charges") # optimized bond lengths, charges, gasphase H2S ri.add_singlejob( "dft/band_h2s.results/ams.rkf", task="GeometryOptimization", properties=["distance(0,1)", "angle(1,0,2)", "charges"], ) # bond and angle scans, gasphase H2S ri.add_singlejob("dft/bondscan_h2s_pbesol.results/ams.rkf", task="PESScan", properties={"pes": {"weight": 3.0}}) ri.add_singlejob("dft/anglescan_h2s_pbesol.results/ams.rkf", task="PESScan", properties="pes") # H2S adsorbed on ZnS(110), bond lengths and charges ri.add_singlejob( "dft/band_h2s_110.results/ams.rkf", task="GeometryOptimization", properties=["distance(34,19)", "distance(14, 37)", "charges"], ) # enthalpy of formations ri.add_reaction_energy( reactants=["dft/sulfur.results", "dft/Zn.results"], products=["dft/ZnS_pesscans/wurtzite_sp"], normalization="p0", normalization_value=0.5, ) ri.add_reaction_energy( reactants=["dft/sulfur.results", "dft/Zn.results"], products=["dft/ZnS_pesscans/rocksalt_sp"], normalization="p0", normalization_value=1.0, ) ri.add_reaction_energy( reactants=["dft/sulfur.results", "dft/Zn.results"], products=["dft/ZnS_pesscans/zincblende_sp"], normalization="p0", normalization_value=1.0, ) # relative polymorph energies ri.add_reaction_energy( reactants=["dft/ZnS_pesscans/wurtzite_sp"], products=["dft/ZnS_pesscans/zincblende_sp"], normalization="p0", normalization_value=1.0, ) ri.add_reaction_energy( reactants=["dft/ZnS_pesscans/rocksalt_sp"], products=["dft/ZnS_pesscans/zincblende_sp"], normalization="p0", normalization_value=1.0, ) # adsorption energy ri.add_reaction_energy( reactants=["dft/band_110.results", "dft/band_h2s.results"], products=["dft/band_h2s_110.results/ams.rkf"], normalization="p0", normalization_value=1.0, task="GeometryOptimization", ) # "surface energy", in a way ri.add_reaction_energy( reactants=["wurtzite_sp"], products=["dft/band_110_noconstraints.results"], normalization="r0", normalization_value=0.5, task="GeometryOptimization", ) ri.add_singlejob( "dft/band_110_noconstraints.results", task="GeometryOptimization", properties=[ "distance(39,42)", "distance(19, 42)", "distance(42,47)", "distance(27,47)", "distance(4,47)", "distance(4,12)", "distance(3,12)", "distance(12,15)", "distance(5,28)", "angle(42,39,19)", "angle(39,19,31)", "angle(39,19,24)", "angle(39,42,47)", "angle(23,27,47)", "angle(27,47,4)", "dihedral(9,28,27,23)", "charges", ], ) # some singlepoints ri.add_singlejob("dft/band_distorted_ads_110.results/", properties="forces") ri.add_singlejob("dft/band_distorted_clean_110.results", properties="forces") ri.add_reaction_energy(reactants=["dft/band_distorted_ads_110.results"], products=["band_h2s_110"], sigma=2.0) ri.add_reaction_energy(reactants=["dft/band_distorted_clean_110.results"], products=["band_110"], sigma=2.0) # save training_set.yaml, validation_set.yaml, job_collection.yaml ri.store(backup=False) if __name__ == "__main__": main()