This research was funded by the National Science Foundation through award number CBET-1438325. Fellowship support for CJB through NASA Training Grant NX14AN43H is gratefully acknowledged. Simulations were performed on the Palmetto Supercomputer Cluster, which is maintained by the Cyberinfrastructure Technology Group at Clemson University. We thank Dr. Paul J. Meza-Morales for testing the protocol.

1. Generate the adsorbate structure
<ol>
	<li>Create a VASP POSCAR file comprising a supercell with periodic boundary conditions, as you would if you were performing simulations of adsorbates on metal surfaces under vacuum. The supercell should include an initial guess of the adsorbate structure and the metal surface as well as vacuum space above the adsorbate for adding H2O molecules. Details on how to do this are provided in previous work35,44. 
	NOTE: It is important that the height of the vacuum space be at least 12 &#197; above the top of the adsorbate.</li>
	<li>Relax the structure and minimize its energy using the VASP code. Details on how to do this are provided in previous work35,44. This will produce a file called CONTCAR, which will be used in the next section.</li>
</ol>
2. Add explicit H2O molecules
<ol>
	<li>Add N water molecules to the vacuum space in the CONTCAR created in step 1.2 using the MCPliQ code, where N = &#961;V, &#961; is the density of water, and V is the volume of the vacuum space above the adsorbate. 
	NOTE: &#961; should be taken as the density of water as determined by the TIP3P/CHARMM water model at the simulation temperature. V will be refined in the next step.
	<ol>
		<li>Specify the following information in the MCPliQ master_input.txt file: The number of H2O molecules to add (N) by changing the first argument in line 28, the path to the water.txt file by changing the second argument in line 28, and the minimum and maximum height of the supercell that can be occupied by water molecules by changing the Minimum z-coordinate in line 11 and the Maximum z-coordinate in line 12.</li>
		<li>Execute the MCPliQ code by typing mcpliq from the command-line interface to insert water molecules into the CONTCAR file. The code will output one or more files with the file extension of .POSCAR. 
		NOTE: If more than one POSCAR files are produced, they will be named POSCAR_n.POSCAR. Select the file where n is largest.</li>
	</ol>
	</li>
	<li>Generate the LAMMPS input files for an NPT simulation and equilibrate the cell volume using FFMD in the NPT ensemble in LAMMPS.
	<ol>
		<li>Execute the script lmps_bond_angle.py on the .POSCAR file generated in step 2.1.2 by typing lmps_bond_angle.py $filename.POSCAR at the command-line interface, where $filename is the name of the .POSCAR file generated in step 2.1.2. This script creates a file called $filename.POSCAR.bond_angle_info.txt which lists bonds and angles that will be used in the LAMMPS data file.</li>
		<li>Open VMD and select File &#62; New Molecule in the main window to open the Molecule File Browser window. Select VASP_POSCAR from the Determine File Type dropdown menu. Click Browse and navigate to the $filename.POSCAR file. Click Load to open the $filename.POSCAR file.</li>
		<li>Open the Tk Console within VMD by selecting Extensions &#62; Tk Console from the VMD Main window.</li>
		<li>Execute the following command in the Tk console: topo writelammpsdata $WDPATH/data.myadsorbate full, where $WDPATH is the directory on the computer where VMD will write the LAMMPS data file and data.myadsorbate is the name of the LAMMPS data file.</li>
		<li>Delete the Bonds and Angles section at the bottom of the data.myadsorbate file. Then, append the bond and angle lists in the file $filename.POSCAR.bond_angle_info.text into data.myadsorbate. 
		NOTE: The indexes for the O-H bond type and H-O-H angle type for the water molecules in the $filename.POSCAR.bond_angle_info.txt file are both set to 1. Thus, bond and angle types for adsorbates should start counting at 2.</li>
		<li>Edit the data.myadsorbate file by adding the Lennard-Jones parameters to the Pair Coeffs section and the Coulomb parameters to the Atoms section. Lennard-Jones and Coulomb parameters for the H2O molecules, adsorbate atoms, and metal surface atoms need to be added. 
		NOTE: The Lennard-Jones parameters for Pt atoms, water molecules, and adsorbate atoms in this protocol are obtained from the UFF38, TIP3P/CHARMM37, and OPLS-AA39 force fields, respectively. Coulomb parameters for water molecules and adsorbates atoms are obtained from the TIP3P/CHARMM37 and OPLS-AA39 force fields, respectively. Coulomb parameters for Pt atoms are set to 0 in this protocol. Alternatively, calculated partial charges could be used for the Coulomb parameters for adsorbate atoms and Pt atoms.</li>
		<li>Copy the LAMMPS input file input.equil into the directory $WDPATH. Edit the group variable on line 34 to indicate the atom type indexes for the water oxygen and water hydrogen atoms and the group variable on line 35 to indicate the atom type indexes for the Pt and adsorbate atoms.</li>
		<li>Execute the LAMMPS software by typing mpiexec -np XX lmp_mpi &#60; input.equil at the command-line interface, where XX is the number of CPU cores to use, and lmp_mpi is the name of the LAMMPS executable. Doing this will run an energy minimization to refine the H2O configuration, followed by a FFMD simulation performed at constant number of H2O molecules (N), volume (V), and temperature (T) to bring the water to the simulation temperature, followed by a FFMD simulation run at constant N, pressure (P), and temperature (T) to determine the physically correct height of the simulation box. Output files that will be used in section 3 are called data.myadsorbate_npt and log.myadsorbate. 
		NOTE: The duration of the NPT simulation should be long enough to comprise an &#8220;equilibration&#8221; run, where the volume of the supercell comes to steady state, and a &#8220;production&#8221; run, which is used to sample the ensemble averages (here, the height of the supercell). During the equilibration run, the volume of the supercell when plotted against time should level off to a steady state value. Once this occurs, the NPT simulation can be said to be in its production run. Verify equilibration of the NPT simulation by ensuring that the fluctuations in the height of the supercell (lz) are minimal or have converged to a steady value. If large fluctuations occur, then re-generate a H2O configuration by decreasing the timestep on line 92 in the input.equil file and repeating step 2.2.8 or starting again from step 2.1.1 .</li>
	</ol>
	</li>
</ol>
3. Extract the proper height of the supercell
<ol>
	<li>Execute the script get_npt_lz.py on the log.myadsorbate file by typing get_npt_lz.py log.myadsorbate at the command-line interface. This script outputs the average supercell height from the &#8220;production&#8221; run portion of the NPT simulation in the avg_lz.txt file. 
	NOTE: The get_npt_lz.py script assumes that LAMMPS writes the length of the cell z-dimension (lz) to the log.myadsorbate file every 1000 fs (customizable at line 20 of the get_npt_lz.py script), which is the default in the provided input.equil LAMMPS input file. The get_npt_lz.py script detects and discards the first 2 ns (customizable at line 19 of the get_npt_lz.py script) worth of lz values in the log.myadsorbate file, as they comprise the equilibration portion of the simulation, while the remaining 3 ns comprise the &#8220;production&#8221; portion and are thus used by the get_npt_lz.py script to compute the average z-dimension length. In addition to the avg_lz.txt file, the get_npt_lz.py script outputs a file called npt_data.txt, which provides values of lz as a function of the timestep, as well as a file called npt_plot.png, which plots the same data. The plot can be used to verify equilibration of the NPT simulation.</li>
	<li>Reconstruct the supercell using the average height determined in NPT.
	<ol>
		<li>Copy the data.myadsorbate_npt file into a new directory, referred to here as $WD2PATH, and rename it data.myadsorbate.</li>
		<li>Edit the new data.myadsorbate file so that the lz height is equal to the average value output from the get_npt_lz.py script by changing the zlo and zhi arguments in the data.myadsorbate file such that zlo is 0.0 and zhi is the lz value from the avg_lz.txt file produced in step 3.1.</li>
	</ol>
	</li>
</ol>
4. Generate configurations of H2O molecules
<ol>
	<li>Copy the LAMMPS input file input.prod into $WD2PATH. Edit the group variable on line 32 to indicate the atom type indexes for the water oxygen and water hydrogen atoms and the group variable on line 33 to indicate the atom type indexes for the Pt and adsorbate atoms.</li>
	<li>Execute the LAMMPS software by typing mpiexec -np XX lmp_mpi &#60; input.prod into the command-line interface, where XX is the number of CPU cores to use, and lmp_mpi is the name of the LAMMPS executable. Doing this will run a constant NVT simulation on the H2O molecules. The key output file from this simulation is the dump.myadsorbate.lammpstrj file. 
	NOTE: The duration of the NVT simulation should be long enough to comprise an equilibration run, where the energy of the system comes to steady state, and a production run, from which the ensemble averages (here, the spatial positions of the water molecules) are sampled. During the equilibration run, the energy of the system when plotted against time should level off to a steady state value. Once this occurs, the NVT simulation can be said to be in its production run.</li>
</ol>
5. Determine the hydrogen bond lifetime for proper time sampling
<ol>
	<li>Edit the hb_lifetime_dist.py script to specify: the timestep of the first frame of the dump.myadsorbate.lammpstrj file by changing the actualStart variable on line 22, how often frames are written to the LAMMPS trajectory file by changing the timestep variable on line 23, the first and last timesteps the script should consider (i.e., the production portion of the trajectory) by changing the N_first and N_last variables on lines 24 and 25, whether consecutive frames are considered or frames are skipped by changing the nevery variable on line 26, and the number of lines per frame section of the trajectory file by changing the frameLine variable on line 27. Additionally, edit lines 31 through 35 to specify which atom types within the data.myadsorbate file belong to the adsorbate and which atom types belong to the H2O molecules. 
	NOTE: The hb_lifetime_dist.py script analyzes the H2O configurations in the production run and determines if any H2O molecules are hydrogen bonded to the adsorbate. It then counts the simulation time that each hydrogen bond remains intact and reports this information as a distribution of hydrogen bond lifetimes in units of ps. The specific version of the script that is provided with this protocol assumes that LAMMPS writes the configuration of H2O molecules to the dump.myadsorbate.lammpstrj file every 1000 fs, which is the default in the provided input.prod LAMMPS input file. It detects and discards the first 2 ns worth of configurations in the the dump.myadsorbate.lammpstrj file, as they comprise the equilibration portion of the simulation, and uses the remaining 3 ns to calculate hydrogen bond lifetimes.</li>
	<li>Execute the script hb_lifetime_dist.py on the dump.myadsorbate.lammpstrj file by typing hb_lifetime_dist.py at the command-line interface. Doing this will produce a file called distribution_HB_lifetime.dat.</li>
	<li>Plot the data in the distribution_HB_lifetime.dat file to view the distribution of hydrogen bond lifetimes that occurred during the NVT simulation.</li>
	<li>Determine the time increment to use for the time sampling interval based on the calculated hydrogen bond lifetimes. The best choice is the maximum hydrogen bond lifetime; alternatively, a value that would capture the 95% confidence interval can be used.</li>
</ol>
6. Sample configurations of liquid H2O molecules
<ol>
	<li>Determine the number of configurations from the production run of the NVT FFMD trajectory for further calculations. The number of configurations should be selected so that the minimum time between configurations is equal to or greater than the time sampling interval identified in section 5.
	<ol>
		<li>Edit the default value for the num_frames variable on line 21 of the lammps_frames.py script to specify the number of configurations to extract.</li>
		<li>Execute the script lammps_frames.py on the file dump.myadsorbate.lammpstrj by typing lammps_frames.py at the command-line interface. Doing this will output a list of simulation times corresponding to the configurations that should be extracted from the dump.myadsorbate.lammpstrj file. These configurations can be used as starting structures in AIMD or QM simulations. 
		NOTE: 1) The lammps_frames.py script automatically detects the LAMMPS log and dump files as well as the production portion of the trajectory within the dump file and divides the number of configurations within the dump file into 10 groups. Alternatively, the user can specify the log file, the dump file, and the number of configurations from the command-line interface using the -l, -d, and -n options, respectively. To do so, the user should type lammps_frames.py -n XX -l $logfilename -d $dumpfilename at the command-line interface, where XX is the desired number of configurations, $logfilename is the name of the LAMMPS logfile, and $dumpfilename is the name of the LAMMPS trajectory (dump) file. The simulation times that are output refer to the median times in each group. 2) If the configurations will be calculated in VASP with the LDIPOLE flag turned on, a small layer of vacuum space should be added to the top of the supercell above the water layer. This will facilitate convergence of the electronic structure in the VASP calculation. Adding an additional 3 &#197; of vacuum space above the H2O molecules has been successful in the simulations discussed below.</li>
	</ol>
	</li>
</ol>

The authors disclose no conflicts of interest.

The method as presented was selected for its ease of implementation, but multiple customizations could be made. For one, the force fields used in the FFMD simulations can be modified. Changing the force field parameters and/or potentials can be done by editing the LAMMPS input and data files. Similarly, solvents other than H2O could be employed. To make this modification, the desired solvent molecule would need to be inserted starting from step 2.1.1, and the LAMMPS input files would need to be edited to incorporate the appropriate potentials and parameters. Inserting the new solvent molecule would also require supplying the internal coordinates of the solvent molecule in a .txt file analogous to the water.txt file.
Another modification that could be made is to modify the area of the surface slab. The results discussed in this manuscript employed 3 Pt x 3 Pt or 4 Pt x 4 Pt surface slabs, which have surface areas less than 120 &#197;2. As the slab surface area increases, the computational expense also increases. Computational expense has the largest impact on section 5 of this protocol. If the data processing steps in section 5 become computationally prohibitive, big data post processing strategies such as those discussed in Li et al. 201845 can be employed.
Possible sources of uncertainty for this procedure include the force field employed, the sampling method, and the sampling frequency. The water structure is determined by the force field that is used, meaning that the choice of force field could influence the specific configurations of H2O molecules. Our group has assessed how the choices of force field for H2O molecules and Pt atoms influence the interaction energies calculated in FFMD and found that the choice of force field contributes less than 0.1 eV to this interaction energy. Another source of uncertainty is the sampling method, which influences the specific configurations that are used to calculate a quantity of interest. Our group has compared the performance of the &#8220;time sampling&#8221; method presented in this protocol with an &#8220;energy sampling&#8221; method, which is biased to lower energy configurations of H2O molecules, on the interaction energies calculated in DFT and found both of these sampling methods give statistically equal values35,46. The sampling frequency can also influence the results. We have assessed how increasing the number of configurations from 10 to 30,000 influences the average interaction energies calculated in FFMD for 40 different C3HxO3 adsorbates and found that the sampling frequency contributes less than 0.1 eV to the average interaction energy44.
The main limitation to this method is that the adsorbates are approximated by the structures under vacuum during the FFMD simulations. In reality, the adsorbates would exhibit conformational changes (bond stretches, angle bends, torsional motions, etc.) due to normal thermal movements, including interactions with solvent molecules. Attempts to include conformational changes of adsorbates into the FFMD simulations would require detailed development of force fields for catalytic surface adsorbates, i.e., which comprise terms that describe bond stretches, angle bends, and torsional terms, amongst others. As a future direction of this protocol, we are developing such force fields for adsorbates at solid surfaces, which we will use to determine the extent to which using rigid adsorbates influences the results.

Modeling molecular phenomena involved in heterogeneous catalysis under liquid conditions is necessary for understanding catalytic function; however, this remains challenging because it requires a fine balance between chemical accuracy and computational expense. In general, since catalysis involves the breaking and forming of chemical bonds, quantum mechanics must be used to at least some degree; however, long simulations are challenging in quantum mechanics, as they require significant computer resources. Since molecules in the liquid phase are under constant thermal motion, simulations must also include configurational sampling, i.e., they must incorporate multiple spatial arrangements of the liquid molecules, as each different spatial arrangement (i.e., each configuration) has a different energy. This means that multiple configurations of liquid molecules must be simulated for each catalytic species of interest. These needs &#8211; to use quantum mechanics and to perform multiple calculations per catalytic species &#8211; can render modeling in heterogeneous catalysis under liquid phase computationally intractable. The purpose of the method described herein is to enable computationally tractable simulations of phenomena in heterogeneous catalysis under liquid phase.
We are particularly interested in heterogeneously catalyzed reactions that are carried out under liquid water. Water molecules have significant influence on catalytic phenomena, such as interacting with catalytic species (e.g., via dispersion forces and hydrogen bonding)1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23, participating in catalytic reactions1,7,8,9,15,21,22,24,25,26,27, and influencing reaction pathways and/or catalytic rates1,11,12,15,18,23,25,27,28,29,30,31. Modeling of these phenomena has been performed using QM and/or ab initio molecular dynamics (AIMD)1,2,6,7,14,22,25,27,28,32,33,34, force field molecular dynamics (FFMD)35, and quantum mechanics/molecular mechanics (QM/MM)10. In AIMD and FFMD, the atoms in the system are moved pursuant to Newton&#8217;s equations of motion according to the forces acting upon them. In AIMD, the system energy and forces are calculated with quantum mechanics, whereas in FFMD, the system energy and forces are calculated using force fields, which are algebraic expressions that are parameterized based on experimental or QM data. In QM/MM, the portion of the system where the bond breaking and forming occurs is calculated with QM, and the remainder of the system is calculated with MM, which employs force fields. Because they directly employ QM, AIMD and QM/MM are better suited for capturing the bond breaking and forming that occurs in aqueous phase heterogeneous catalysis; however, FFMD is significantly more computationally tractable and thus better suited for generating the configurations of liquid H2O molecules. The method presented in this protocol balances chemical accuracy and computational expense by employing a combination of QM and FFMD.
Specifically, this method uses FFMD simulations for generating configurations of liquid H2O and QM to calculate system energies. FFMD is carried out using LAMMPS.36 The force fields used in FFMD in this work employ Lennard-Jones + Coulomb (LJ+C) potentials, where the LJ parameters have been taken from the TIP3P/CHARMM model37 for H2O, the universal force field38 (UFF) for Pt, and the OPLS-AA force field39 for catalytic species, and the Coulomb parameters have been taken from the TIP3P/CHARMM37 model for H2O and the OPLS-AA force field39 for catalytic species. The Coulomb parameters for Pt atoms have been set to 0. QM calculations are performed using the VASP code40,41,42, which is a density functional theory (DFT) code. Water molecule insertions are performed with a code developed in-house called Monte Carlo Plug-in for Quantum Methods (MCPliQ). File conversions from VASP to LAMMPS in this protocol are performed with the Visual Molecular Dynamics (VMD) software43.
The protocol is intended to generate configurations of liquid water molecules around catalytic species on flat transition metal surfaces at low coverage. Coverage is denoted &#952;&#160;and defined as the number of adsorbates per surface metal atom (i.e., the number of surface adsorbates normalized by the number of metal atoms in the topmost layer of the metal slab in the catalyst model). In this manuscript, low coverage is defined as &#952; &#8804; 1/9 monolayer (ML), where 1 ML means one catalytic species per surface metal atom. The catalyst models should be placed in periodic simulation boxes. The simulation boxes do not have to be cubes. This manuscript demonstrates the use of the protocol for generating configurations of liquid H2O that can be used to calculate quantities of interest in aqueous phase heterogeneous catalysis.
This protocol requires that the user has access to installed and working versions of the VASP, MCPliQ, LAMMPS, and VMD software. More information about VASP (https://www.vasp.at/), LAMMPS (https://Lammps.sandia.gov/), and VMD (https://www.ks.uiuc.edu/Research/vmd/) are available on their respective websites. The MCPliQ software is documented at https://github.com/getman-research-group/JoVE_article, along with all input files and Python scripts mentioned in this protocol. This protocol assumes that the executables and scripts mentioned within will be run on a high-performance research computer and are installed in a directory that is in the user&#8217;s $PATH variable. If an executable or script is placed in a location that is not in the user&#8217;s $PATH, then the path to the executable must be included to execute it. Executables and scripts are executed in steps 2.1.2, 2.2.1, 2.2.8, 3.1, 4.2, 5.2, and 6.1.2. For example, to execute the MCPliQ code in step 2.1.2 from a directory that is not in the user&#8217;s $PATH, the user would type $PATHTOMCPLIQ/mcpliq at the command-line interface instead of mcpliq, where $PATHTOMCPLIQ is the location where the mcpliq executable has been stored (e.g., $PATHTOMCPLIQ might be ~/bin). Before starting this protocol, all executables and scripts should be given executable permissions (e.g., in Linux, this could be done by typing chmod +x mcpliq at the command-line interface from the directory where the mcpliq executable is stored). Further, any modules required by any of the software or scripts should be loaded (these dependencies will be specific to individual installations of the various software and the computer where the simulations will be run).

<table>
	<tbody>
		<tr>
			<td>VASP software</td>
			<td>Computational Materials Physics, Dept. of Physics, University of Vienna</td>
			<td>vasp.5.4.4</td>
			<td>Standard parallel VASP executable in the newest version.</td>
		</tr>
		<tr>
			<td>LAMMPS software</td>
			<td>Sandia National Laboratory</td>
			<td>31Mar17-dp</td>
			<td>Double-precision, parallel LAMMPS executable from 31 March 2017.</td>
		</tr>
		<tr>
			<td>VMD software</td>
			<td>Theoretical and Computational Biophysics Group, University of Illinois at Urbana-Champaign</td>
			<td>1.9.3</td>
			<td>Standard VMD executable in the newest version.</td>
		</tr>
		<tr>
			<td>MCPliQ software</td>
			<td>Getman Research Group, Dept. of Chemical and Biomolecular Engineering, Clemson University</td>
			<td></td>
			<td>Executable and input files for the MCPliQ software availabe from the Getman Research Group GitHub page.</td>
		</tr>
		<tr>
			<td>JoVE article scripts</td>
			<td>Getman Research Group, Dept. of Chemical and Biomolecular Engineering, Clemson University</td>
			<td></td>
			<td>Python scripts for this JoVE manuscript available from the Getman Research Group GitHub page.</td>
		</tr>
		<tr>
			<td>H2O PDB file</td>
			<td>Getman Research Group, Dept. of Chemical and Biomolecular Engineering, Clemson University or RCSB Protein Data Bank</td>
			<td></td>
			<td>PDB file for a water molecule, available from the Getman Research Group GitHub page or at http://www.rcsb.org/ligand/HOH.</td>
		</tr>
	</tbody>
</table>

multiscale sampling of a heterogeneous water/metal catalyst interface using density functional theory and force-field molecular dynamics

A significant number of heterogeneously-catalyzed chemical processes occur under liquid conditions, but simulating catalyst function under such conditions is challenging when it is necessary to include the solvent molecules. The bond breaking and forming processes modeled in these systems necessitate the use of quantum chemical methods. Since molecules in the liquid phase are under constant thermal motion, simulations must also include configurational sampling. This means that multiple configurations of liquid molecules must be simulated for each catalytic species of interest. The goal of the protocol presented here is to generate and sample trajectories of configurations of liquid water molecules around catalytic species on flat transition metal surfaces in a way that balances chemical accuracy with computational expense. Specifically, force field molecular dynamics (FFMD) simulations are used to generate configurations of liquid molecules that can subsequently be used in quantum mechanics-based methods such as density functional theory or ab initio molecular dynamics. To illustrate this, in this manuscript, the protocol is used for catalytic intermediates that could be involved in the pathway for the decomposition of glycerol (C3H8O3). The structures that are generated using FFMD are modeled in DFT in order to estimate the enthalpies of solvation of the catalytic species and to identify how H2O molecules participate in catalytic decompositions.

One use of this protocol is to calculate the energies of interaction between the liquid water and catalytic species, i.e., &#916;Eint35:
&#8710;Eint=ECatalytic species+H2O+EClean catalyst surface-ECatalytic species-EClean catalyst surface+H2O
where ECatalytic species+H2O is the energy of a configuration of H2O molecules around a catalytic species on a metal surface, EClean catalyst surface is the energy of the clean catalyst surface in vacuum, ECatalytic species is the energy of the catalytic species on a metal surface in vacuum, and EClean catalyst surface+H2O is the energy of the configuration of H2O over the catalyst surface with the catalytic species removed. The positions of the H2O molecules used to calculate ECatalytic species+H2O and EClean catalyst surface+H2O should be identical. All values of E are calculated using the VASP code. The quantity &#916;Eint includes all of the physical and chemical interactions between all of the molecules in the liquid water structure and the catalytic species and gives a reasonable estimate of the enthalpy of solvation of the catalytic species, which is needed to calculate its free energy of solvation and total free energy. Table 1 provides values for &#916;Eint calculated for species on a Pt(111) catalyst surface with chemical formulas equal to CxHyOz in units of eV (1 eV = 96.485 kJ/mol). The values were calculated at coverages &#8804;1/9 ML.35,46&#160;The values reported are the averages taken over 10 configurations of liquid H2O, and the uncertainties are reported as standard deviations. All the values are negative, indicating favorable interactions with water.
Another application of this protocol is to generate starting structures for AIMD. Movie 1 is a movie of an AIMD trajectory that was started from a configuration generated by this protocol. At the start of this movie, a COH adsorbate is shown on a Pt(111) surface under a structure of liquid H2O. One H2O molecule is emphasized, which formed a hydrogen bond with COH. Over the course of the movie, this H2O molecule abstracts the proton from the COH adsorbate and deposits a second hydrogen atom on the Pt(111) surface. The H2O molecule thus helps to catalyze the reaction COH* + * &#8594;&#160;CO* + H*, where the *s indicate catalytic sites. This simulation highlights the main strength and the main purpose of the multiscale sampling method described herein. Numerous configurations of H2O molecules are generated with FFMD, due to its strength in computational tractability. However, a limitation of FFMD is that it cannot capture bond breaking and forming unless a reactive force field is implemented. AIMD uses quantum mechanics to calculate energies and thus can capture bond breaking and forming. However, AIMD is too computationally demanding to generate all of the configurations of H2O molecules necessary to ensure sufficient sampling has been achieved. Thus, this protocol combines the two methods.
The structures of liquid H2O molecules generated by this procedure are dependent on input settings. Setting these improperly can have unintended influences on the water structures. For example, when the intermolecular distances become too small or when other parameters in the molecular dynamics input files are set improperly or take on unphysical values, the water structure can become unreasonable. Under these circumstances, the structure of water will &#8220;blow up&#8221; unintendedly during the FFMD trajectory. Figure 1 shows an example of this. The snapshot on the left-hand side is the starting structure for a FFMD run, and the snapshot on the right-hand side is a snapshot taken within 1 ps of starting the simulation. As can be seen, the H2O molecules have moved far away from the surface. This is caused by improper settings made in the simulation input files and is not a structure that is likely to occur in reality.
<img alt="Figure 1" class="xfigimg" src="/files/ftp_upload/59284/59284fig1v2.jpg" /> 
Figure 1: Example of a negative result. The force field molecular dynamics simulation &#8220;blows up&#8221; due to an unphysical setting or value. Left hand image: The starting geometry of the Pt(111) surface, adsorbate, and liquid water structure. Right hand image: The geometry of the Pt(111) surface, adsorbate, and liquid water structure less than 1 ps later. In the right-hand image, the H2O molecules have separated from the surface due to unphysically large forces. <a href="https://www.jove.com/files/ftp_upload/59284/59284fig1v2large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>
<img alt="Movie 1" src="/files/ftp_upload/59284/59284mov1.jpg" /> 
Movie 1: Ab initio molecular dynamics (AIMD) simulation initiated from a configuration generated in multiscale sampling. A H2O molecule that is originally hydrogen bonded to a COH adsorbate on a Pt(111) surface abstracts the proton from COH and deposits a second hydrogen on the Pt(111) surface. This bond breaking and forming event can be captured by AIMD but not with force field molecular dynamics (FFMD) unless a reactive force field is used. The initial configuration of H2O molecules used in this AIMD simulation was generated using FFMD as described in this manuscript. <a href="https://www.jove.com/files/ftp_upload/59284/movie1.mp4" target="_blank">Please click here to view this video. (Right-click to download.)</a>
<table border="1" fo:keep-together.within-page="1" fo:keep-with-next.within-page="always">
	<tbody>
		<tr>
			<td>Catalytic Species</td>
			<td>&#8710;Eint (eV)</td>
		</tr>
		<tr>
			<td>COH</td>
			<td>-0.70 &#177; 0.07</td>
		</tr>
		<tr>
			<td>CO&#160;</td>
			<td>-0.03 &#177; 0.03</td>
		</tr>
		<tr>
			<td>CH2OH</td>
			<td>-0.64 &#177; 0.12</td>
		</tr>
		<tr>
			<td>CHO-CHOH-CH2OH</td>
			<td>-0.93 &#177; 0.22</td>
		</tr>
		<tr>
			<td>COH-COH-CH2OH</td>
			<td>-0.87 &#177; 0.23</td>
		</tr>
		<tr>
			<td>COH-CHOH-COH</td>
			<td>-1.72 &#177; 0.26</td>
		</tr>
		<tr>
			<td>CHOH-COH-CO</td>
			<td>-1.57 &#177; 0.25</td>
		</tr>
		<tr>
			<td>CHO-CO-CO</td>
			<td>-0.31 &#177; 0.19</td>
		</tr>
	</tbody>
</table>
Table 1: Water-catalytic species interaction energy results. Interaction energies in eV calculated for eight CxHyOz adsorbates on Pt(111). Values reported are the averages taken over multiple configurations of liquid H2O. The uncertainties are the standard deviations of the averages. 1 eV = 96.485 kJ/mol.

Watch this Scientific Journal Video about Multiscale Sampling of a Heterogeneous Water/Metal Catalyst Interface using Density Functional Theory and Force-Field Molecular Dynamics at JoVE.com

Multiscale Sampling of a Heterogeneous Water/Metal Catalyst Interface using Density Functional Theory and Force-Field Molecular Dynamics

The goal of the protocol presented here is to generate and sample trajectories of configurations of liquid water molecules around catalytic species on a flat transition metal surface. The sampled configurations can be used as starting structures in quantum mechanics-based methods.

A significant number of heterogeneously-catalyzed chemical processes occur under liquid conditions, but simulating catalyst function under such conditions ...

multiscale-sampling-heterogeneous-watermetal-catalyst-interface-using

Research

JoVE Journal

Chemistry

12.6K Views.  Clemson University. The goal of the protocol presented here is to generate and sample trajectories of configurations of liquid water molecules around catalytic species on a flat transition metal surface. The sampled configurations can be used as starting structures in quantum mechanics-based methods.

Video: Multiscale Sampling of a Heterogeneous Water/Metal Catalyst Interface using Density Functional Theory and Force-Field Molecular Dynamics

HKUST-1 as a Heterogeneous Catalyst for the Synthesis of Vanillin

Vanillin (4-hydoxy-3-methoxybenzaldehyde) is the main component of the extract of vanilla bean. The natural vanilla scent is a mixture of approximately 200 different odorant compounds in addition to vanillin. The natural extraction of vanillin (from the orchid Vanilla planifolia, Vanilla tahitiensis and Vanilla pompon) represents only 1% of the worldwide production and since this process is expensive and very long, the rest of the production of vanillin is synthesized. Many biotechnological approaches can be used for the synthesis of vanillin from lignin, phenolic stilbenes, isoeugenol, eugenol, guaicol, etc., with the disadvantage of harming the environment since these processes use strong oxidizing agents and toxic solvents. Thus, eco-friendly alternatives on the production of vanillin are very desirable and thus, under current investigation. Porous coordination polymers (PCPs) are a new class of highly crystalline materials that recently have been used for catalysis. HKUST-1 (Cu3(BTC)2(H2O)3, BTC = 1,3,5-benzene-tricarboxylate) is a very well known PCP which has been extensively studied as a heterogeneous catalyst. Here, we report a synthetic strategy for the production of vanillin by the oxidation of trans-ferulic acid using HKUST-1 as a catalyst.

The conversion of trans-ferulic acid to vanillin was achieved by heterogeneous catalysis. HKUST-1 was employed in this synthesis and the essential step in the catalytic process was the generation of unsaturated metal sites. Thus, when the catalyst was activated under vacuum, full vanillin conversion (yield of 95%) was obtained.

Vanillin (4-hydoxy-3-methoxybenzaldehyde) is the main component of the extract of vanilla bean. The natural vanilla scent is a mixture of approximately ...

Facile Synthesis of Worm-like Micelles by Visible Light Mediated Dispersion Polymerization Using Photoredox Catalyst

Presented herein is a protocol for the facile synthesis of worm-like micelles by visible light mediated dispersion polymerization. This approach begins with the synthesis of a hydrophilic poly(oligo(ethylene glycol) methyl ether methacrylate) (POEGMA) homopolymer using reversible addition-fragmentation chain-transfer (RAFT) polymerization. Under mild visible light irradiation (&#955; = 460 nm, 0.7 mW/cm2), this macro-chain transfer agent (macro-CTA) in the presence of a ruthenium based photoredox catalyst, Ru(bpy)3Cl2 can be chain extended with a second monomer to form a well-defined block copolymer in a process known as Photoinduced Electron Transfer RAFT (PET-RAFT). When PET-RAFT is used to chain extend POEGMA with benzyl methacrylate (BzMA) in ethanol (EtOH), polymeric nanoparticles with different morphologies are formed in situ according to a polymerization-induced self-assembly (PISA) mechanism. Self-assembly into nanoparticles presenting POEGMA chains at the corona and poly(benzyl methacrylate) (PBzMA) chains in the core occurs in situ due to the growing insolubility of the PBzMA block in ethanol. Interestingly, the formation of highly pure worm-like micelles can be readily monitored by observing the onset of a highly viscous gel in situ due to nanoparticle entanglements occurring during the polymerization. This process thereby allows for a more reproducible synthesis of worm-like micelles simply by monitoring the solution viscosity during the course of the polymerization. In addition, the light stimulus can be intermittently applied in an ON/OFF manner demonstrating temporal control over the nanoparticle morphology.

This article describes a process for producing polymeric self-assembled nanoparticles using visible light mediated dispersion polymerization. Using low energy visible light to control the polymerization allows for the reproducible formation of self-assembled worm-like micelles at high solids content.

Presented herein is a protocol for the facile synthesis of worm-like micelles by visible light mediated dispersion polymerization. This approach begins ...

Synthesis and Performance Characterizations of Transition Metal Single Atom Catalyst for Electrochemical CO2 Reduction

This protocol presents both the synthesis method of the Ni single atom catalyst, and the electrochemical testing of its catalytic activity and selectivity in aqueous CO2 reduction. Different from traditional metal nanocrystals, the synthesis of metal single atoms involves a matrix material that can confine those single atoms and prevent them from aggregation. We report an electrospinning and thermal annealing method to prepare Ni single atoms dispersed and coordinated in a graphene shell, as active centers for CO2 reduction to CO. During the synthesis, N dopants play a critical role in generating graphene vacancies to trap Ni atoms. Aberration-corrected scanning transmission electron microscopy and three-dimensional atom probe tomography were employed to identify the single Ni atomic sites in graphene vacancies. Detailed setup of electrochemical CO2 reduction apparatus coupled with an on-line gas chromatography is also demonstrated. Compared to metallic Ni, Ni single atom catalyst exhibit dramatically improved CO2 reduction and suppressed H2 evolution side reaction.

Synthesis and Performance Characterizations of Transition Metal Single Atom Catalyst for Electrochemical CO2 Reduction

Here, we present a protocol for the synthesis and electrochemical testing of transition metal single atoms coordinated in graphene vacancies as active centers for selective carbon dioxide reduction to carbon monoxide in aqueous solutions.

This protocol presents both the synthesis method of the Ni single atom catalyst, and the electrochemical testing of its catalytic activity and selectivity ...

Preparation of Functional Silica Using a Bioinspired Method

The goal of the protocols described herein is to synthesize bioinspired silica materials, perform enzyme encapsulation therein, and partially or totally purify the same by acid elution. By combining sodium silicate with a polyfunctional bioinspired additive, silica is rapidly formed at ambient conditions upon neutralization. 
The effect of neutralization rate and biomolecule addition point on silica yield are investigated, and biomolecule immobilization efficiency is reported for varying addition point. In contrast to other porous silica synthesis methods, it is shown that the mild conditions required for bioinspired silica synthesis are fully compatible with the encapsulation of delicate biomolecules. Additionally, mild conditions are used across all synthesis and modification steps, making bioinspired silica a promising target for the scale-up and commercialization as both a bare material and active support medium.
The synthesis is shown to be highly sensitive to conditions, i.e., the neutralization rate and final synthesis pH, however tight control over these parameters is demonstrated through the use of auto titration methods, leading to high reproducibility in reaction progression pathway and yield. 
Therefore, bioinspired silica is an excellent active material support choice, showing versatility towards many current applications, not limited to those demonstrated here, and potency in future applications.

Here, we present a protocol to synthesize bioinspired silica materials and immobilize enzymes therein. Silica is synthesized by combining sodium silicate and an amine &#39;additive&#39;, which neutralize at a controlled rate. Material properties and function can be altered either by in situ enzyme immobilization or post-synthetic acid elution of encapsulated additives.

The goal of the protocols described herein is to synthesize bioinspired silica materials, perform enzyme encapsulation therein, and partially or totally ...

Heterogeneous Removal of Water-Soluble Ruthenium Olefin Metathesis Catalyst from Aqueous Media Via Host-Guest Interaction

A highly efficient transition metal catalyst removal method is developed. The water-soluble catalyst contains a newly-designed NHC ligand for the catalyst removal via host-guest interactions. The new NHC ligand possesses an adamantyl (guest) tethered linear ethylene glycol units for hydrophobic inclusion into the cavity of a &#946;-cyclodextrin (&#946;-CD) host compound. The new NHC ligand was applied to a Ru-based olefin metathesis catalyst. The Ru catalyst demonstrated excellent activity in representative ring-closing metathesis (RCM) and ring-opening metathesis polymerization (ROMP) reactions in aqueous media as well as organic solvent, CH2Cl2. After the reaction was complete, the lingering Ru residue was removed from the aqueous solution with the efficiency of more than 99% (53 ppm of Ru residue) by simple filtration utilizing a host-guest interaction between insoluble silica-grafted &#946;-CD (host) and the adamantyl moiety (guest) on the catalyst. The new Ru catalyst also demonstrated high removal efficiency via extraction when the reaction is run in organic solvent by partitioning the crude reaction mixture between layers of diethyl ether and water. In this way, the catalyst stays in aqueous layer only. In organic layer, the residual Ru amount was only 0.14 ppm in the RCM reactions of diallyl compounds.

A removable water-soluble N-heterocyclic carbene (NHC) ligand in aqueous media via host-guest interaction has been developed. We demonstrated representative olefin metathesis reactions in water as well as in dichloromethane. Via either host-guest interaction or extraction, the residual ruthenium (Ru) catalyst was as low as 0.14 ppm after the reaction.

A highly efficient transition metal catalyst removal method is developed. The water-soluble catalyst contains a newly-designed NHC ligand for the catalyst ...

Imaging Corrosion at the Metal-Paint Interface Using Time-of-Flight Secondary Ion Mass Spectrometry

Corrosion developed at the paint and aluminum (Al) metal-paint interface of an aluminum alloy is analyzed using time-of-flight secondary ion mass spectrometry (ToF-SIMS), illustrating that SIMS is a suitable technique to study the chemical distribution at a metal-paint interface. The painted Al alloy coupons are immersed in a salt solution or exposed to air only. SIMS provides chemical mapping and 2D molecular imaging of the interface, allowing direct visualization of the morphology of the corrosion products formed at the metal-paint interface and mapping of the chemical after corrosion occurs. The experimental procedure of this method is presented to provide technical details to facilitate similar research and highlight pitfalls that may be encountered during such experiments.

Time-of-flight secondary ion mass spectrometry is applied to demonstrate the chemical mapping and corrosion morphology at the metal-paint interface of an aluminum alloy after being exposed to a salt solution compared with a specimen exposed to air.

Corrosion developed at the paint and aluminum (Al) metal-paint interface of an aluminum alloy is analyzed using time-of-flight secondary ion mass ...

[(DPEPhos)(bcp)Cu]PF6: A General and Broadly Applicable Copper-Based Photoredox Catalyst

Our group recently reported the use of [(DPEPhos)(bcp)Cu]PF6 as a general copper-based photoredox catalyst which proved efficient to promote the activation of a broad variety of organic halides, including unactivated ones. These can then participate in various radical transformations such as reduction and cyclization reactions, as well as in the direct arylation of several (hetero)arenes. These transformations provide a straightforward access to a range of small molecules of interest in synthetic chemistry, as well as to biologically active natural products. Altogether, [(DPEPhos)(bcp)Cu]PF6 acts as a convenient photoredox catalyst which appears to be an attractive, cheap and complementary alternative to the state-of-the-art iridium- and ruthenium-based photoredox catalysts. Here, we report a detailed protocol for the synthesis of [(DPEPhos)(bcp)Cu]PF6, as well as NMR and spectroscopic characterizations, and we illustrate its use in synthetic chemistry for the direct arylation of (hetero)arenes and radical cyclization of organic halides. In particular, the direct arylation of N-methylpyrrole with 4-iodobenzonitrile to afford&#160;4-(1-methyl-1H-pyrrol-2-yl)benzonitrile and the radical cyclization of N-benzoyl-N-[(2-iodoquinolin-3-yl)methyl]cyanamide to afford natural product luotonin A are detailed. The scope and limitations of this copper-based photoredox catalyst are also briefly discussed.

[DPEPhosbcpCu]PF6: A General and Broadly Applicable Copper-Based Photoredox Catalyst

Detailed and general protocols are presented for the synthesis of [(DPEPhos)(bcp)Cu]PF6, a general copper-based photoredox catalyst, and for its use in synthetic chemistry for the direct arylation of C-H bonds in (hetero)arenes and radical cyclization of organic halides.

Our group recently reported the use of [(DPEPhos)(bcp)Cu]PF6 as a general copper-based photoredox catalyst which proved efficient to promote the ...

Vibrational Spectra of a N719-Chromophore/Titania Interface from Empirical-Potential Molecular-Dynamics Simulation, Solvated by a Room Temperature Ionic Liquid

The accurate molecular simulation prediction of vibrational spectra, and other structural, energetic and spectral characteristics, of photo-active metal-oxide surfaces in contact with light-absorbing dyes is an ongoing thorny and elusive challenge in physical chemistry. With this in mind, a molecular-dynamics (MD) simulation were performed by using optimized empirical potentials for a well-representative and prototypical dye-sensitized solar cell (DSC) solvated by a widely-studied room temperature ionic liquid (RTIL), in the guise of a [bmim]+[NTf2]- RTIL solvating an N719-sensitizing dye adsorbed onto 101 anatase-titania. In doing so, important insights were gleaned into how using a RTIL as the electrolytic hole acceptor modulates the dynamical and vibrational properties of a N719 dye, estimating the spectra for the DSC photo-active interface via Fourier transformation of mass-weighted velocity autocorrelation functions from MD. The acquired vibrational spectra were compared with the experiment spectra and those sampled from ab initio molecular dynamics (AIMD); in particular, various empirical-potential spectra generated from MD provide insight into how partial-charge charge parameterization of the ionic liquid affects vibrational spectra prediction. In any event, careful fitting of empirical force-field models has been shown to be an effective tool in handling DSC vibrational properties, when validated by AIMD and an experiment.

A dye-sensitized solar cell was solvated by RTILs; using optimized empirical potentials, a molecular dynamics simulation was applied to compute vibrational properties. The obtained vibrational spectra were compared with experiment and ab initio molecular dynamics; various empirical potential spectra show how partial-charge charge parameterization of the ionic liquid affects vibrational spectra prediction.

The accurate molecular simulation prediction of vibrational spectra, and other structural, energetic and spectral characteristics, of photo-active ...

Development of Heterogeneous Enantioselective Catalysts using Chiral Metal-Organic Frameworks (MOFs)

Substrate size discrimination by the pore size and homogeneity of the chiral environment at the reaction sites are important issues in the validation of the reaction site in metal-organic framework (MOF)&#8211;based catalysts in an enantioselective catalytic reaction system. Therefore, a method of validating the reaction site of MOF-based catalysts is necessary to investigate this issue. Substrate size discrimination by pore size was accomplished by comparing the substrate size versus the reaction rate in two different types of carbonyl-ene reactions with two kinds of MOFs. The MOF catalysts were used to compare the performance of the two reaction types (Zn-mediated stoichiometric and Ti-catalyzed carbonyl-ene reactions) in two different media. Using the proposed method, it was observed that the entire MOF crystal participated in the reaction, and the interior of the crystal pore played an important role in exerting chiral control when the reaction was stoichiometric. Homogeneity of the chiral environment of MOF catalysts was established by the size control method for a particle used in the Zn-mediated stoichiometric reaction system. The protocol proposed for the catalytic reaction revealed that the reaction mainly occurred on the catalyst surface regardless of the substrate size, which reveals the actual reaction sites in MOF-based heterogeneous catalysts. This method for reaction site validation of MOF catalysts suggests various considerations for developing heterogeneous enantioselective MOF catalysts.

Development of Heterogeneous Enantioselective Catalysts using Chiral Metal-Organic Frameworks MOFs

Here, we present a protocol for active site validation of metal-organic framework catalysts by comparing stoichiometric and catalytic carbonyl-ene reactions to find out whether a reaction takes place on the inner or outer surface of metal-organic frameworks.

Substrate size discrimination by the pore size and homogeneity of the chiral environment at the reaction sites are important issues in the validation of ...

Analyzing Melts and Fluids from Ab Initio Molecular Dynamics Simulations with the UMD Package

We have developed a Python-based open-source package to analyze the results stemming from ab initio molecular-dynamics simulations of fluids. The package is best suited for applications on natural systems, like silicate and oxide melts, water-based fluids, and various supercritical fluids. The package is a collection of Python scripts that include two major libraries dealing with file formats and with crystallography. All the scripts are run at the command line. We propose a simplified format to store the atomic trajectories and relevant thermodynamic information of the simulations, which is saved in UMD files, standing for Universal Molecular Dynamics. The UMD package allows the computation of a series of structural, transport and thermodynamic properties. Starting with the pair-distribution function it defines bond lengths, builds an interatomic connectivity matrix, and eventually determines the chemical speciation. Determining the lifetime of the chemical species allows running a full statistical analysis. Then dedicated scripts compute the mean-square displacements for the atoms as well as for the chemical species. The implemented self-correlation analysis of the atomic velocities yields the diffusion coefficients and the vibrational spectrum. The same analysis applied on the stresses yields the viscosity. The package is available via the GitHub website and via its own dedicated page of the ERC IMPACT project as open-access package.

Melts and fluids are ubiquitous vectors of mass transport in natural systems. We have developed an open-source package to analyze ab initio molecular-dynamics simulations of such systems. We compute structural (bonding, clusterization, chemical speciation), transport (diffusion, viscosity) and thermodynamic properties (vibrational spectrum).

We have developed a Python-based open-source package to analyze the results stemming from ab initio molecular-dynamics simulations of fluids. The package ...