Computation of Atmospheric Concentrations of Molecular Clusters from ab initio Thermochemistry

Our protocol provides a flexible and computationally feasible approach to studying weakly-bound molecular clusters and can be readily applied to gain insights into their structure, formation, and abundance. This technique's main advantage is its efficiency and flexibility in treating molecular clusters at different levels of theory from quick force-fields and semi-empirical methods to rigorous quantum mechanical methods. Atmospheric and aerosol chemistry can benefit the most from this approach leading to better models of climate change.

However, any field involving molecular clusters can exploit this approach. For individuals who have never performed this technique, the most challenging steps are the initial program and script installation and their adaptation to the local computing environment. Students new to computational chemistry can overcome the steep learning curves in using high-performance computing clusters through the visualization of explicit step-by-step instructions.

To obtain a minimum energy structure of isolated glycine molecules for use in a genetic algorithm configurational sampling, open a new session in Avogadro and click Build, Insert, Peptide, Glycine, and Insert Peptide to generate a glycine monomer in the visualization window. Click Extensions and Gaussian and edit the first line in the text box as indicated. Click Generate and save the command file as glycine.com.

To obtain a minimum energy structure of isolated water, open a new session in Avogadro and select Build, Insert, and Fragment. Enter water into the filter text box, select the water file and click Insert. Click Extensions and Gaussian and edit the first line in the text box as indicated.

Click Generate and save the command file as water.com. Then transfer the two com files to the computing cluster and run the Gaussian 09 calculations using the appropriate submit script. Once the calculations have finished on the computing cluster, call open babel to generate xyz files of the minimum energy structures entering the command as indicated.

For genetic algorithm-based configurational sampling, add all of the scripts and templates into a folder and copy the folder to the remote cluster. Make sure that all of the scripts are executable and use the commands as indicated to add the location of the scripts directory to the path environmental variable. To obtain a set of low-energy structures for glycine and water at the inexpensive semi-empirical level of theory, create a directory called gly-h2o-n for which n is the number of water molecules and create a subdirectory called GA under the gly-h2o-n directory to run genetic algorithm calculations.

Copy the ogolem input files, monomers Cartesian coordinates, and PBS batch submission script into the GA directory and run the GA calculation using the appropriately modified run. pbs submit script. Once the calculation is complete, change the directory to gly-h2o-n GA pm7 and run the getRotConsts command as indicated where 13 is the number of atoms in the cluster and zero and nine indicate that there are 10 structures with indices zero through nine.

This will compute the rotational constants of the GA optimized clusters and generate a file called rotConstsData_C that contains a sorted list of all the GA optimized cluster configurations, their energies and their rotational constants. Run the similarityAnalysis. py script with rotConstsData_C file as an input to find and save the unique GA optimized clusters.

Pm7 will be used as a file naming label to generate a file called uniqueStructures-pm7.data. This contains a sorted list of the unique GA optimized configurations. In the gly-h2o-n GA directory, use the combine-GA.

csh script to combine the results for multiple comparable GA runs and to generate a new unique structures list named uniqueStructures-pm7. data in the gly-h2o-n GA directory. The working directory should have the exact organization and structure as illustrated.

To refine the structures of the glycine water clusters from the genetic algorithm based on a semi-empirical method to one using a more accurate quantum mechanical method, create a subdirectory called QM under the gly-h2o-n directory. Under the QM directory, create another subdirectory named pw91-sb and copy the uniqueStructures list from the gly-h2o-n GA directory to the QM pw91-sb directory. Change the directory to gly-h2o-n QM pw91-sb and run the small basis set density functional theory script run-pw91-sb.

csh for which sb is a label for this set of calculations, Q is the preferred queue on the computing cluster, and 10 indicates that 10 calculations will be grouped into one batch job. Once the submitted calculations are complete, use the getRotConsts-dft-sb. csh script to extract the energies and compute the rotational constants of the small basis optimized clusters.

Here, pw91 is the density functional used and n is the number of atoms in the cluster. Use the similaryAnalysis. py script as before to identify the unique structures but use sb as the label.

A list of unique configurations optimized at the pw91 631 plus G star level of theory will be saved in the uniqueStructures-sb. data file. In the gly-h2o-n QM directory, use the combined combine-QM.

csh script to combine the results from multiple comparable QM runs. The combine-QM. csh pw91-sb command will generate a new unique structures list named uniqueStructures-sb.

data in the gly-h2o-n QM directory. To further refine the structures of the glycine and water clusters using a better quantum mechanical description, create a subdirectory called pw91-lb under the QM directory. Copy the unique structures list from the QM pw91-sb directory to the QM pw91-lb directory and change the directory to QM pw91-lb.

Run the large basis density functional theory script run-pw91-lb. csh for which lb is a label for this set of calculations, Q is the preferred queue on the computing cluster, and 10 indicates that 10 calculations are to be grouped into one batch job. Once the submitted calculations are complete, use the getRotConsts-dft-lb.

csh command to compute the rotational constants of the large basis optimized clusters. Here, pw91 is the density functional used and n is the number of atoms in the cluster. Use the similarityAnalysis.

py script as before now with lb as the label to generate a list of unique configurations optimized at the pw91 6311 plus plus G star star theory level and save in the uniqueStructures-lb. data file. To obtain the vibrational structure and energies of the glycine and water clusters necessary to compute the desired thermochemical corrections, copy the unique structures list from the QM pw91-lb directory to the QM pw91-lb ultrafine directory and change the directory to QM/pw91-lb ultrafine.

Run the ultrafine density functional theory script run-pw91-lb-ultrafine. csh for which uf is a label for this set of calculations, Q is the preferred queue on the computing cluster, and 10 indicates that 10 calculations are to be grouped into one batch job. This script will automatically generate the inputs for Gaussian 09 and submit all of the calculations.

Once the submitted calculations are complete, use the getRotConsts-dft-lb-ultrafine. csh command to compute the rotational constants of the ultrafine optimized clusters. Here, pw91 is the density functional used and n is the number of atoms in the cluster.

Use the similarityAnalysis. py script as before now with uf as the label to generate and save a list of unique configurations optimized to ultrafine convergence criteria at the pw91 6311 plus plus G star star theory level in the uniqueStructures-uf. data file.

Then run the run-thermo-pw91. csh script with uniqueStructures-uf. data file as the input to compute the thermodynamic corrections.

Copy and paste the command line output to the attached spreadsheet named gly-h2o-n.xls. As the raw energies of this calculation and the subsequent n equals two, three, four, and five calculations are added to the first sheet of the gly-h2o-n. xls spredsheet, the hydrate distribution sheet which yields the equilibrium concentration of hydrates at different temperatures, relative humidity, and initial concentrations of water and glycine will be updated.

Here, lowest electronic energy isomers of glycine-water clusters can be observed. Note how the hydrogen bond network grows in complexity as the number of water molecules increases moving from a mostly planar network to a three-dimensional cage-like structure at n equals five. In this table, an example of the output of the run-thermo-pw91.

csh script is shown. For each cluster, the energy of the pw91 6311 plus plus G star star corresponds to the gas phase electronic energies at the pw91 6311 plus plus G star star level of theory calculated on ultrafine integration grids in units of hartrees as well as the zero point vibrational energy in units of kilocalorie per mole. At each temperature, the calculated enthalpy formation delta H in a Gibbs-free energy formation delta G are given in units of kilocalories per mole and the calculated enthalpy formation S is given in units of calories per mole.

In this table, representative computations of the total Gibbs-free energy change of hydration and of the sequential hydration are shown. Using these data, the atmospheric concentrations of hydrated glycine can be calculated. One must install the correct software and added the included scripts to reflect one's own computing environment.

Adding the location of the scripts to one's path is crucial. This technique was used to determine the catalytic activity of atmospheric water clusters towards peptide bond formation to contribute to the field of prebiotic chemistry.