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

The protocol models, thermal and quantum phenomena in liquid-phase heterogeneous catalysis. It is the first, to our knowledge, to incorporate quantum mechanics with complete sampling of an explicit liquid environment. Configurations of liquid molecules generated from this protocol represent those expected under actual reaction conditions and can be used to explore molecular level phenomena that depend on the spatial arrangements of molecules.

The configurations of liquid molecules generated by this protocol provide insight into the roles that solvent plays on liquid phase heterogeneous catalysis. If you're trying this protocol for the first time, I recommend that you first ensure that you have access to and can run the VASP, MCPliQ, VMD, and LAMMPS software. Demonstrating the procedure will be Tianjun Xie, a graduate student from my laboratory.

After generating the adsorbate structure, generate the LAMMPS input files for an NPT simulation and equilibrate the cell volume using FFMD. Copy the LAMMPS input file into the working directory. Edit the group variable on line 34 to indicate the atom type indexes for the water oxygen and water hydrogen atoms.

In the group variable on line 35, to indicate the atom type indexes for the platinum and adsorbate atoms. On line 17 of the input file, edit the run step variable to set the duration of the NPT simulation long enough to comprise an equilibration run and a production run. Execute the LAMMPS software, by typing the command at the command-line interface, which includes the information of the number of CPU cores to use and the name of the LAMMPS executable.

The energy minimization refines the water molecule configuration and is followed by a FFMD simulation performed at constant number of water molecules, volume, and temperature to bring the water to the simulation temperature. Another FFMD simulation is then run at constant number of water molecules, pressure, and temperature to determine the physically correct height of the simulation box. The output files will be used later.

After the NPT simulation, plot the height of the supercell against time. The point where it levels off to a steady state value is the point in the NPT simulation when the production can begin. Verify equilibration of the NPT simulation by ensuring that the fluctuations in the height of the supercell are minimal or have converged to a steady value.

If large fluctuations occur, open the input. equil file and decrease the timestep on line 92 to regenerate the water molecule configuration and execute the LAMMPS software as before. To begin, type at the command-line interface to execute the script.

This script outputs the average supercell height from the production run portion of the NPT simulation to a TXT file. The script reads the length of the cell's z-dimension at 1, 000 femtosecond intervals, which is the default interval for printing information in LAMMPS. If a different printing interval is desired, it can be changed by editing line 20 of the get_npt_lz.

py script and line 16 of LAMMPS input. equil file. The script detects and discards the first two nanoseconds worth of lz values, as they comprise the equilibration portion of the simulation.

The duration of the equilibration run can be changed by editing line 19 of the file. The remaining three nanoseconds comprise the production portion, and are thus used to compute the average z-dimension length. In addition, the script outputs another TXT file which provides values of lz as a function of the timestep as well as a PNG file, which plots the same data.

The plot can be used to verify equilibration of the NPT simulation. To reconstruct the supercell using the average height determined in NPT, copy the previously generated data file into a new working directory and rename it as data.myadsorbate. Then edit the new data file to change zlo to 0.0 and zhi to the lz value from the average value output in the TXT file.

Copy the LAMMPS input file into the new working directory. 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 platinum and adsorbate atoms. Then on line 16, edit the runStep variable so that it is long enough to comprise an equilibration run and a production run.

Type the command to execute LAMMPS into the command-line interface to execute the LAMMPS software. This will run a constant NVT simulation on the water molecules, and the key output file is generated. The NVT simulation comprises an equilibration portion and a production portion.

The production portion starts when the energy of the system plotted against time levels off. To begin, open the hydrogen bond lifetime script. Change the actual start variable on line 22 to set the timestep of the first timeframe.

Change the timestep variable on line 23 to set how often frames are written to the LAMMPS trajectory file. Change N_first and N-last variables on line 24 and 25 to set the first and last timesteps, and change the nevery variable on line 26 to set whether consecutive frames are considered or skipped. Set the number of lines per frame section of the trajectory file by changing the frame line 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 water molecules. The script analyzes the water configurations in the production run and determines if any water 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 picoseconds. LAMMPS writes the configuration of water molecules to the file every 1, 000 femtoseconds, which is the default in the provided LAMMPS input file. It detects and discards the first two nanoseconds worth of configurations in the file as they comprise the equilibration portion of the simulation and uses the remaining three nanoseconds to calculate hydrogen-bond lifetimes.

To execute the script, type at the command-line interface. It then generates a DAT file. Plot the data in the file to view the distribution of hydrogen-bond lifetimes that occurred during the NVT simulation.

To determine the time increment to use for the time sampling interval, use a time increment greater than or equal to the maximum hydrogen-bond lifetime. Determine the number of configurations from the production run of the NVT FFMD trajectory to sample such that the minimum time between configurations is equal to or greater than the time sampling interval identified previously. On the previously written frame extraction script, edit the default value for the number of frames variable on line 21 to specify the number of configurations to extract.

To execute the script, type the script name at the command-line interface. This will output a list of simulation times corresponding to the configurations that should be extracted from the NVT simulation file. These configurations can be used as starting structures in AIMD or QM simulations.

In this procedure, FFMD was used to generate the initial configuration of water molecules. The AIMD simulation shows that a water molecule that is originally hydrogen-bonded to a sugar alcohol adsorbate on a platinum-111 surface abstracts the hydrogen from the alcohol adsorbate and deposits a second hydrogen on the platinum-111 surface. The structures of liquid water molecules are dependent on input settings.

Setting these improperly can have unintended influences on the water structures. In this figure, the left side is the starting structure for a FFMD run. And the right side is within one picosecond of starting the simulation.

The FFMD simulation blows up due to unphysically large force settings, causing water molecules to move far away from the surface. The configurations can be used in the quantum mechanics, or QMM simulation, or they could be used to analyze statistics related to the spatial positions of the molecules. This technique paves the way for researchers to explore the roles that liquid reaction environments have on catalysis by generating actual configurations of liquid molecules at catalytic interfaces.