Cluster formation, Single-point energies & Optimizations
Obtaining crystals
Different siliceous/zeolite structures are found at IZA.
The structure that is used to create the initial cluster is: MFI.
1. Setting up the QM/MM embedded-cluster
The first thing that must be done is cutting a large siliceous framework into a smaller siliceous cluster. Siliceous simply means that the structure contains only Si and O atoms. Each Si atom is bound to 4 O atoms, and each O atom is bound to 2 Si atoms. As such, the structure contains this repeating unit:
The siliceous cluster is constructed by running the input file construct_mfi.py, located in the directory zeol_covalent/1_construct_silicate.
from chemsh import *
# Creation of mfi unit cell fragment
mfi = Fragment(coords='mfi.c', connect_mode='covalent')
mfi.save('mfi.xyz','xyz')
# Creation of mfi cluster
cluster = mfi.construct_cluster(radius_cluster=25.0, origin_atom=4, adjust_charge='coordination_scaled',
radius_active=10.0, bq_margin=5.0, bq_density=2)
# Saving of mfi cluster
cluster.save('mfi.pun', 'pun')
cluster.save('mfi.pun.xyz', 'xyz')
The periodic input structure mfi is created as a fragment with coordinates read from mfi.c.
By using the construct_cluster() function on the fragment, the large mfi framework is cut to produce an mfi cluster. The variables included in this function are:
radius_cluster: Determines the size of the cluster. This is in units Bohr (\(a_{0}\))origin_atom: The central atom that the cluster is based around.adjust_charge: Adjusts the charges of the covalently bonded atoms at the periphery surface of the cluster to account for the bonds that were cleaved as a result of the cluster generationradius_active: Sets the size of the sampling region. This should be at least as large as the largest optimisation region you intend to use.bq_marginandbq_density: Set the number and position of the point charges around the edge of the cluster to approximate the missing periodic electrostatic interactions. (‘F’ is arbitrarily used to denote this inside the cluster).
Any structure can be saved by running the save() function, giving a filename and filetype. In our example, we have used .pun files, due to our method of defining the qm/mm regions in the structure, but alternatives include .cjson and .xyz.
The resulting cluster from the construct_cluster() function is saved as mfi.pun. This file contains the coordinates, charge and connectivity of all of the atoms in the cluster. We will continue to use this file in the production of a zeolite cluster as well as further siliceous cluster calculations. Therefore, it has been copied across to succeeding directories.
2. Conversion to zeolite
Now that a silicate base has been established, we can make the necessary changes to produce a zeolite by running the construct_zsm5.py file, located in the directory zeol_covalent/2_construct_zeol. Below is an explanation of the code that carries out this procedure.
Before any ammendments can be made, the mfi cluster fragment is produced using coordinates read from mfi.pun. After this, the mfi cluster fragment is displayed using the print_frag() function:
# Creation of mfi cluster fragment
mfi = Fragment(coords='mfi.pun', connect_mode='covalent')
# Displays mfi cluster before amendments
print("The structure of the silicate is:")
print_frag(mfi)
The first step in converting the silicate to a zeolite is exchanging an Si atom for an Al atom. All of the atom labels and atomic numbers in the cluster can be accessed/altered with names and znums, specifying the index. As python counts from 0, the first atom has index 0, and thus the 551st atom has index 550.
mfi.names[0] = 'Al'
mfi.znums[0] = 13
However, if we compare the Al centre of the cluster now with the Si centre of the cluster before the substitution…
As you can see, the charges of the hypothetical complex are not balanced. Si has a charge of 4+ in the fragment, and each O has a charge of 2-. Thus, with a Si atom bound to 4 O atoms, the overall charge is 4-. Al has a charge of 3+ in the fragment. Bound to 4 O atoms, the overall charge is 5-.
Therefore, we must add a H atom the structure to compensate for the charge. This H addition is carried out via the insert() function using positional arguments.
mfi.insert('H', 550)
The coordinates and charge of the newly inserted atom are updated respectively to replicate that of the Si atom that was deleted.
mfi.coords[550] = -1.65333239913411e+01, -6.10680708115438e+00, 8.97525903483508e+00
mfi.charges[550] = 0.00000
The updated mfi cluster fragment is displayed using the print_frag() function, and the contents saved in a new file zsm5.pun using the save() function. As this is needed in further zeolite scripts, it has been copied across to succeeding directories.
3. Single-point energy calculation of the silicate
A single-point energy calculation can be carried out on the silicate by running the input file mfi_sp.py found in the directory zeol_covalent/3_silicate_sp.
Like in section 2, the fragment mfi is created using coords read from mfi.pun.
After we cut the cluster, a quantum mechanical (qm) region and molecular mechanical (mm) region is formed. Each region is processed using a different calculator. We specify how each of the regions are dealt with using qm_theory and mm_theory.
qm_theory = NWChem(method='scf', basis='3-21g', maxiter=100)
mm_theory = GULP(ff=ff)
qm_theory uses NWChem as the calculator and thus uses NWChem arguments. The qm_theory method scf and the basis file 3-21g have been chosen as the default due to their speed.
mm_theory uses GULP as the calculator and thus uses GULP arguments; the forcefield used for gulp is zeolite_gulp, which is introduced first.
If you would like to obtain a more accurate value, at the cost of a more lengthy procedure, changes can be made to the qm_theory and mm_theory .i.e. changing the parameters to include dft in the qm_theory method with conjugate=True in mm_theory.
# Method of dft used instead of scf
qm_theory = NWChem(method='dft',
functional='b3lyp',
basis='3-21g',
charge=0,
unique_listing=True,
path="mpirun -np 4 nwchem")
# Additional conjugate parameter
mm_theory = GULP(ff=ff,
conjugate=True)
Now, we must define the atoms that in the qm region that will be used in the NWChem calculator.
qm_region = mfi.getRegion(1)
All atoms in the QM region are defined by a ‘1’ suffix on their labels in the zsm5.pun file. The getRegion(1) function places these selected atoms in a list called qm_region.
Note
If you wish to change the atoms selected to be part of the QM region, open the zsm5.pun file using a text editor. This can be done by entering the command vi zsm5.pun. To edit the file, enter ‘i’ and make the necessary changes. When the changes are complete, exit the editor mode by entering esc. Leave the visualiser and save the changes by entering ‘:wq’.
We introduce a dipole-adjust which adjusts the charge of the qm/mm boundary atoms for a more accurate energy value. Each cluster will have it’s own specific dipole-adjustment depending on the atoms present, which is set using bond_modifiers.
silicate_modifiers = {("si","o1"):0.3}
The above parameters all come together to form the qmmm object which is necessary to run the energy calculation. There are a few additional arguments defined within the object. The use of a dipole adjust can be initiated if the parameter dipole_adjust is set to true.
qmmm = QMMM(qm=qm_theory,
mm=mm_theory,
frag=mfi,
coupling='covalent',
qm_region=qm_region,
bond_modifiers=silicate_modifiers,
embedding='electrostatic',
dipole_adjust=True)
Finally, the sp object is prepared using the theory qmmm as an argument. This is what is used to calculate the single-point (sp) energy of the cluster when we apply the run() function.
sp = SP(theory=qmmm, gradients=False)
sp.run()
4. Single-point energy calculation of the zeolite
A single-point energy calculation can be carried out on a zeolite cluster by running the file zsm5_sp.py found in the directory zeol_covalent/4_zeolite_sp.
Instead of using an mfi cluster fragment created from mfi coordinates, the zeolite cluster fragment zsm5 is created using coords read from zsm5.pun.
zsm5 = Fragment(coords='zsm5.pun', connect_mode='covalent')
The rest of the python script is exactly the same as that used in calculating a silicate cluster sp energy, with zsm5 replacing mfi every time the fragment is called.
5. Optimization of silicate
We can optimize the clusters that we have produced to find their lowest energy state. During the optimization process, the atoms in the cluster will change their coordinates slightly with each successive cycle that runs until the lowest energy state is found, or the energy has ‘converged’.
Optimization of the silicate can be done by running the input file mfi_opt.py. This is found in the directory zeol_covalent/5_silicate_opt.
The layout replicates that of the single point energy: the qm_theory, mm_theory, qm_region, silicate_modifiers and qmmm object are all defined. After, there are a few minor differences.
To begin, we must define the atoms that will be frozen during the optimization process. These atoms are any that lie outside of the qm_region. As such, we use the getNotQm() function to define these atoms.
frozen = getNotQM(mfi.natoms, qm_region)
This function uses the number of atoms in the cluster fragment (accessed via natoms) as well as the current qm_region as the parameters. The resulting atoms are all placed into a list, called frozen. To access this function, we must import it from the file getFrozen.
from chemsh.cluster.getFrozen import getNotQM
Instead of using the sp object to find the single-point energy of the cluster, we use the opt object to optimize the cluster.
opt = Opt(theory=qmmm, algorithm="lbfgs", tolerance=0.0045, trust_radius="energy", maxcycle=500, maxene=500, frozen=frozen) opt.run()
The opt object includes the following parameters:
theory: Specifies where the object gets its atomic data from. Set to qmmm to use data from theqmmmobject.algorithm: Specifies the optimisation algorithm. Here it is defaulted as lbfgs.tolerance: Determines how narrow the convergence criteria is. The smaller the tolerance, the more cycles the optimization will require. A wide tolerance may only require 1 cycle, and there may be no change in the structure/energy at all.maxcycle: Determines the number of optimisation cycles. When the energy has converged the cycles will stop automatically.maxene: Determines the maximum number of energy and gradient calculations per cycle.frozen: Specifies which atoms will be frozen during the calculation. Uses list produced bygetNotQm()function.
The optimization calculation is carried out using the run() function on the opt object.
After the structure has been optimized, the reverse dipole adjust is applied to the fragment using the opposite adjustments.
qmmm.boundary_charge_adjust(reverse_silicate_modifiers)silicate_modifiers = {("si","o1"):0.3} reverse_silicate_modifiers = {("si","o1"):-0.3}
As you can see, the python script includes not only the silicate_modifiers, but the reverse_silicate_modifiers too. This is so we do not have a double dipole-adjustment applied to the cluster in further calculations.
The optimized structure is saved in the file optimized_mfi.pun using the save() function. We will use this file from now on in following procedures to produce more precise mfi fragments, and thus it has been copied across to succeeding directories.
6. Optimization of zeolite
Optimization of the zeolite can be executed by running the input file zsm5_opt.py in the directory zeol_covalent/6_zeol_opt.
Instead of using an mfi cluster fragment created from mfi coordinates, the zeolite fragment zsm5 is created using coords read from zsm5.pun.
zsm5 = Fragment(coords='zsm5.pun', connect_mode='covalent')
The rest of the python script is exactly the same as that used in calculating a silicate cluster sp energy, with zsm5 replacing mfi every time the fragment is called.
The resulting optimized structure is saved in the file, optimised_zsm5.pun, which will be used in all following zsm5 calculations for more precise fragments.