3.2. ReaxFF (advanced): ZnS¶
This tutorial illustrates ReaxFF parametrization for ZnS and H₂S on ZnS(110).
The tutorial teaches how to
Combine several different force fields from the AMS ReaxFF library into a good initial guess
Select training jobs to enter the training set
Note
The input files have already been prepared. If you simply want to run the
parametrization without preparing the input, copy the contents of
$AMSHOME/scripting/scm/params/examples/ZnS_ReaxFF
to a new empty
directory. ($AMSHOME
is the AMS installation directory).
To open the example in the graphical user interface (GUI): SCM → ParAMS,
File → Open the job_collection.yaml
file in the new directory.
You can then run the parametrization.
3.2.1. Mix force fields from the AMS ReaxFF library¶
The ReaxFF force fields included with AMS cover a large part of the periodic table, however, far from all interactions are included.
In this tutorial, the goal is to model H₂S adsorption on a ZnS surface. No ZnS ReaxFF force field is included with AMS, but the following force fields might help:
Force field |
Chemistry |
Reference |
LiS.ff |
Sulfur |
|
Mue2016.ff |
Organic S-containing compounds |
|
ZnOH.ff |
Zn and ZnO surfaces |
A ReaxFF force field contains parameters in the following categories: general (GEN), atomic (ATM), bond (BND), angle (ANG), off-diagonal (OFD), torsional (TOR), hydrogen bond (HBD).
Here, we’ll use the following initial parameters from the above force fields:
LiS.ff: ATM:S, BND:S.S, ANG:S.S.S
Mue2016.ff: ATM:H, BND:S.H, ANG:H.S.H
ZnOH.ff: GEN, ATM:Zn, BND:H.H, BND:Zn.H, BND:Zn.Zn, BND:Zn.O*, ANG:Zn.Zn.O*, ANG:Zn.O*.Zn, ANG:Zn.O*.O*, ANG:O*.Zn.O*, HBD:O*.H.O*
For ZnOH.ff, the asterisks in e.g. BND:Zn.O* indicate that we will instead copy those parameters to initialize the BND:Zn.S block. That is, the Zn-O interaction in ZnOH.ff will be used as first guess for the Zn-S interaction in the new force field.
Here, we will demonstrate how to construct an initial ReaxFF force field from scratch. Open ParAMS in ReaxFF mode:
In the Parameters tab, you will see the ReaxFF parameters in the GEN block. Those parameters always exist for any ReaxFF force field.
We want to copy the GEN parameters from ZnOH.ff:
GEN:
click OK.ZnOH.ff
and click OKNow, add the ATM:S, BND:S.S, and ANG:S.S.S parameters from LiS.ff:
BND:S.S ANG:S.S.S
and click OK. This adds the BND:S.S
and ANG:S.S.S
blocks to the table, and also the ATM:S
block because it didn’t exist.ANG:S.S.S
, ATM:S from ATM:S
, BND:S.S from BND:S.S
. Click OKRegularly save your work:
reaxff-zns.params
Next, add the ATM:H, BND:S.H and ANG:H.S.H from Mue2016.ff:
BND:S.H ANG:H.S.H
and click OK. This adds the BND:H.S
and ANG:H.S.H
blocks to the table, and also the ATM:H
block because it didn’t exist. The bond block is called BND:H.S
(sorted alphabetically).ANG:H.S.H
, ATM:H from ATM:H
, BND:H.S from BND:H.S
. Click OKNote
For the BND block the order of atoms do not matter (BND:H.S is the same as BND:S.H). For the ANG block the order does matter - the central atom is given in the middle, but for example ANG:H.S.Zn is equivalent to ANG:Zn.S.H
Now add the remaining blocks ATM:Zn, BND:H.H, BND:Zn.H, BND:Zn.Zn, BND:Zn.S, ANG:Zn.Zn.S, ANG:Zn.S.Zn, ANG:Zn.S.S, ANG:S.Zn.S
BND:H.H BND:Zn.H BND:Zn.Zn BND:Zn.S ANG:Zn.Zn.S ANG:Zn.S.Zn ANG:Zn.S.S ANG:S.Zn.S HBD:S.H.S
and click OK.S
with O
in the “from” fields. Click OK.This copies the values from the BND:Zn.O block in ZnOH.ff to the BND:Zn.S block in your new force field, copies the values from the ANG:Zn.Zn.O block in ZnOH.ff to the ANG:Zn.Zn.S block in your new force field, etc.
This is accomplished in the get_initial_parameters()
method of generate_input_files.py
:
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
3.2.2. Choose the parameters to optimize¶
We will not optimize (deactivate)
Any GEN (general) parameters
Any ATM (atomic) parameters
Any parameters with category
DoNotOptimize
(these should never be optimized)Any BND (bond) parameters for H-Zn or H-H (there is no training data for those interactions)
Any BND (bond) parameters for π-bonds (there is no training data)
Any ANG (angle) parameters for S-S-S (there is no training data)
Any parameters with value 0 (usually means “off”)
We will optimize (activate):
Block |
Atoms |
Parameters |
BND |
S-S, H-S, Zn-Zn, S-Zn |
|
ANG |
H-S-H, S-Zn-Zn, Zn-S-Zn, S-S-Zn, S-Zn-S |
|
HBD |
S-H-S |
|
There are various ways to select the above list of parameters for optimization.
For example, the BND parameters can be selected as follows:
S.S:
D_e^sigma
p_be1
p_ovun1
p_be2
p_bo1
p_bo2
H.S:
, Zn.Zn:
, and S.Zn:
The ANG parameters can be selected as follows:
theta
p_val1
, p_val2
,``p_val7``, p_val4
The HBD parameters can be selected as follows:
Verify that you’ve activated the desired parameters:
You should now see only the active parameters.
By default, the Min and the Max values are set to ±20% of the value. Let’s increase the range for all the now active parameters by another 8.3%:
S.S:D_e^sigma
)S.H.S:-p_hb3
)This increases the range (decreases Min and increases Max) for all selected parameters by 8.3%. In the Adjust Min and Max dialog, you can Scale/Shift/Set the Range/Min/Max, for many selected parameters at the same time.
Save your work before continuing:
The function activate_parameters()
in generate_input_files.py
applies some logic to decide whether to activate a parameter during the
parametrization:
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
For the BND parameters, we choose not to optimize any parameters that should always be constant(with category DoNotOptimize
), the π-bond parameters, any parameters with value 0, or the Zn-H and H-H parameters.
Similar conditions are placed for the ANG parameters.
3.2.3. ZnS training and validation set¶
The training data contains
Polymorphs
Energy-volume scans of the zinc blende and rock salt polymorphs of ZnS (wurtzite energy-volume scan in the validation set). See the GFN1-xTB tutorial for how to set up an energy-volume scan reference calculation in AMS and import it into ParAMS.
Enthalpy of formation of the zinc blende, wurtzite, and rocksalt polymorphs of ZnS
Gaseous H₂S
Optimized S-H bond length in H₂S, and H-S-H angle. See the GFN1-xTB tutorial or Import training data (GUI) tutorial for how to import the bond length from a geometry optimization job.
Bond scan for the S-H bond and angle scan for the H-S-H angle in gaseous H₂S. See the ReaxFF (basic): H₂O bond scan tutorial for how to set up a bond scan reference calculation in AMS.
Clean ZnS(110) surface
Optimized bond lengths and angles at a clean ZnS(110) surface. See also: Building Crystals and Slabs and Crystals and Surfaces
Formation energy (surface energy) of the ZnS(110) surface. See also: Surface energy
Forces on distorted ZnS(110) surfaces (*)
H₂S adsorbed on ZnS(110)
Optimized Zn-S and H…S bond lengths for H₂S/ZnS(110)
Adsorption energy of H₂S/ZnS(110). See also: Adsorption
Other
Atomic charges for a variety of systems
(*): The distorted ZnS(110) surfaces were obtained by geometry-optimizing the clean ZnS(110) surface using a previous iteration of the force field.
Note
This tutorial does not cover how to set up all the reference calculations. In the GUI, if you switch to the Jobs panel and double-click in the Detail column, you can use the button Edit in AMSinput to see how the reference calculation was set up (approximately, some details like the maximum number of iterations for the geometry optimization may differ).
If you have run a reference calculation, in the ParAMS GUI you can choose Jobs → Add Job, and browse to the ams.rkf file from the job. That will let you import many properties (energy, forces, stress tensor, charges, …).
The needed .yaml files can be found in
$AMSHOME/scripting/scm/params/examples/ZnS_ReaxFF
, where
$AMSHOME
is the AMS installation directory.
Note
Using File → Merge you can add jobs, training set entries, and engines into the current project. If you select File → Open, you will start a new project (thus losing the parameters that you set up in the previous section!).
job_collection.yaml
file.training_set.yaml
file.validation_set.yaml
file.job_collection_engines.yaml
file.First download and unzip reference_calculation_results.zip
.
It contains trimmed-down results files from BAND DFT calculations.
Then take a look at the import_results()
function in generate_input_files.py
.
It uses a Results Importer to import results from different types of DFT reference calculations, for example PES Scans and geometry optimizations.
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, binary=False)
3.2.4. Run the ZnS ReaxFF parametrization¶
Run
"$AMSBIN/params"
3.2.5. ZnS parametrization results¶
Note
These results only give an example of a possible result. If you rerun the parametrization, your results will likely differ and could either be better or worse.
3.2.5.1. ZnS energies¶
Reference (eV) |
Prediction (eV) |
|
ZnS enthalpy of formation (zincblende) |
-2.737 |
-2.737 |
ZnS enthalpy of formation (rock salt) |
-2.329 |
-2.360 |
H₂S adsorption energy on ZnS(110) (a) |
-0.742 |
-0.559 |
ZnS(110) slab energy relative to bulk (b) |
+0.124 |
+0.257 |
(a): This adsorption energy was not calculated using the in-cell approach
(b): This energy can be converted into a surface energy