Name: computation of atmospheric concentrations of molecular clusters from ab initio thermochemistry
Uploaded: 2020-04-08T15:00:00
Duration: 12 min 11 s
Description: Watch this Scientific Journal Video about Computation of Atmospheric Concentrations of Molecular Clusters from ab initio Thermochemistry at JoVE.com

This project was supported by grants CHE-1229354, CHE-1662030, CHE-1721511, and CHE-1903871 from the National Science Foundation (GCS), the Arnold and Mabel Beckman Foundation Beckman Scholar Award (AGG), and the Barry M. Goldwater Scholarship (AGG). High-performance computing resources of the MERCURY Consortium (http://www.mercuryconsortium.org) were used.

1. Finding the minimum energy structure of isolated glycine and water
NOTE: The goal here is two-fold: (i) to obtain minimum energy structures of isolated water and glycine molecules for use in the genetic algorithm configurational sampling, (ii) and to compute the thermodynamic corrections to the gas phase energies of these molecules for use in the calculation of atmospheric concentrations.
<ol>
	<li>On your local computer, open a new session of Avogadro.
	<ol>
		<li>Click Build &#62; Insert &#62; Peptide and select Gly from the Insert Peptide window to generate a glycine monomer in the visualization window.</li>
		<li>Click Extensions &#62; Gaussian and edit the first line in the text box to read &#39;# pw91pw91/6-311++G** int(Acc2E=12,UltraFine) scf(conver=12) opt(tight,maxcyc=300) freq&#39;. Click Generate and save the input file as glycine.com.</li>
		<li>Please note that if the molecule has significant conformational flexibility, as glycine does55, it is critical to perform conformational analysis to identify the global minimum structure and other low-lying conformers. OpenBabel54 provides robust conformational search tools utilizing different algorithms and quick force fields. While conformers are allowed to relax and interconvert during GA and subsequent calculations, it is sometimes necessary to run multiple GA calculations, each starting with a different conformer.</li>
	</ol>
	</li>
	<li>On your local computer, open a new session in Avogadro.
	<ol>
		<li>Click Build &#62; Insert &#62; Fragment and search for &#34;water&#34; from the Insert Fragment window to get the coordinates of water.</li>
		<li>Click Extensions &#62; Gaussian and edit the first line in the text box to read &#39;# pw91pw91/6-311++G** int(Acc2E=12,UltraFine) scf(conver=12) opt(tight,maxcyc=300) freq&#39;. Click Generate and save the input file as water.com.</li>
	</ol>
	</li>
	<li>Transfer the two .com files to the remote cluster. Once you log into the remote cluster, call Gaussian 09 in a batch submission script to start the calculation. When the calculations finish, extract the Cartesian coordinates (.xyz files) of the minimum energy structures by calling OpenBabel. For glycine, the command to execute is: 
	obabel -ig09 glycine.log -oxyz &#62; glycine.xyz 
	These two .xyz files will be used by the GA configurational sampling in the next step.</li>
</ol>
2. Genetic-algorithm-based configurational sampling of Gly(H2O)n=1-5 clusters
NOTE: The goal here is to obtain a set of low-energy structures for Gly(H2O)n=1-5 at the inexpensive semi-empirical level of theory, using the PM746 model implemented in MOPAC47. It is imperative that the working directory has the exact organization and structure as shown in Figure 2. This is to ensure that the custom shell and Python scripts work without failures.
<ol>
	<li>Copy all necessary scripts to the remote cluster and add their location to $PATH
	<ol>
		<li>Put all scripts and template files to a folder (Eg. scripts) and copy it to the remote cluster</li>
		<li>Make sure all the scripts are executable</li>
		<li>Add the location of the scripts directory to the $PATH environmental variable by entering the following commands in a terminal. The default location of the scripts is set to $HOME/JoVE-demo/scripts, however, one can define an environmental variable called $SCRIPTS_HOME pointing to the directory containing the scripts and add $SCRIPTS_HOME to one&#39;s path
		<ol>
			<li>Bash shell: 
			export SCRIPTS_HOME=/path/to/scripts 
			export PATH=${SCRIPTS_HOME}:${PATH}</li>
			<li>Tcsh/Csh shell: 
			setenv SCRIPTS_HOME /path/to/scripts 
			setenv PATH ${SCRIPTS_HOME}:${PATH}</li>
		</ol>
		</li>
	</ol>
	</li>
	<li>On the remote cluster, set up and run a GA calculation:
	<ol>
		<li>Create a directory called gly-h2o-n where n is the number of water molecules.</li>
		<li>Create a subdirectory called GA under the gly-h2o-n directory to run genetic algorithm calculations.</li>
		<li>Copy the OGOLEM input files (Eg. pm7.ogo), monomers Cartesian coordinates (Eg. glycine.xyz, water.xyz) and PBS batch submission script (Eg. run.pbs) into the GA directory.</li>
		<li>Make the necessary changes to the OGOLEM input file and batch submission file.</li>
		<li>Submit the calculation. When the calculation starts, OGOLEM will create a new directory named as the prefix of the OGOLEM input file (Eg. pm7) in the GA directory and store newly generated coordinates there.</li>
	</ol>
	</li>
	<li>Once the calculation is complete, compile the energies and rotational constants, and use that information to determine which are the unique low-energy structures:
	<ol>
		<li>Change directory to gly-h2o-n/GA/pm7 and</li>
		<li>Extract the energies and compute the rotational constants of the GA-optimized clusters with the command: 
		getRotConsts-GA.csh N 0 99 
		where N is the number of atoms in the molecular cluster and &#39;0 99&#39; indicates that the GA pool size is 100, with indices running from 0 through 99. This will generate a file called rotConstsData_C which contains a sorted list of all the GA-optimized cluster configurations, their energies, and their rotational constants.</li>
		<li>Execute the command: 
		similarityAnalysis.py pm7 rotConstsData_C 
		where pm7 will be used as a file-naming label, to find and save the unique GA-optimized clusters. This will generate a file called uniqueStructures-pm7.data which contains a sorted list of the unique GA-optimized configurations. This is a list of unique local minimum structures for the Gly(H2O)n cluster optimized at the PM7 level of theory, and these structures are now ready to be refined using DFT.</li>
	</ol>
	</li>
	<li>Go up to the gly-h2o-n/GA directory and combine the results from multiple comparable GA runs using the combine-GA.csh script. The syntax is: 
	combine-GA.csh &#60;label&#62; &#60;list of directories with GA runs&#62; 
	In this particular case, the command: 
	combine-GA.csh pm7 pm7 
	will generate a new unique structures list named &#39;uniqueStructures-pm7.data&#39; in the gly-h2o-n/GA directory.</li>
</ol>
3. Refinement using QM method with a small basis set
NOTE: The goal here is to refine the configurational sampling of the Gly(H2O)n=1-5 clusters using a better quantum-mechanical description to obtain a smaller but more accurate set of Gly(H2O)n=1-5 cluster structures. The starting structures for this step are the outputs of Step 2.
<ol>
	<li>Prepare and run the small basis set DFT calculation:
	<ol>
		<li>Create a subdirectory called QM under the gly-h2o-n directory. Under the QM directory, create another subdirectory named pw91-sb.</li>
		<li>Copy the unique structures list (uniqueStructures-pm7.data) from the gly-h2o-n/GA directory to the QM/pw91-sb directory.</li>
		<li>Change directory to that gly-h2o-n/QM/pw91-sb.</li>
		<li>Run the small basis set DFT configurational sampling script using the command: 
		run-pw91-sb.csh uniqueStructures-pm7.data sb QUEUE 10 
		where sb is a label for this set of calculations, QUEUE 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 the calculations. Enter &#39;test&#39; for the &#39;QUEUE&#39; to do a dry run.</li>
	</ol>
	</li>
	<li>Once the submitted calculations are complete, extract and analyze the results.
	<ol>
		<li>Extract the energies and compute the rotational constants of the small-basis-optimized clusters using the command: 
		getRotConsts-dft-sb.csh pw91 N 
		where pw91 indicates that the PW91 density functional was used, and N is the number of atoms in the cluster. That will create a file named rotConstsData_C.</li>
		<li>Now identify the unique structures with the command: 
		similarityAnalysis.py sb rotConstsData_C 
		where sb is used as a file-naming label. There will now be a list of unique configurations optimized at the PW91/6-31+G* level of theory saved in the file uniqueStructures-sb.data.</li>
	</ol>
	</li>
	<li>Go up to the gly-h2o-n/QM directory and combine the results from multiple comparable QM runs using the combine-QM.csh script. The syntax is: 
	combine-QM.csh &#60;label&#62; &#60;list of directories with QM calcs&#62; 
	In this particular case, the command: 
	combine-QM.csh sb pw91-sb 
	will generate a new unique structures list named &#39;uniqueStructures-sb.data&#39; in the gly-h2o-n/QM directory.</li>
</ol>
4. Further refinement using QM method with a large basis set
NOTE: The goal here is to further refine the configurational sampling of the Gly(H2O)n=1-5 clusters using a better quantum-mechanical description. The starting structures for this step are the outputs of Step 3.
<ol>
	<li>Submit more reliable calculations using a larger basis set.
	<ol>
		<li>Create a subdirectory called pw91-lb under the QM directory.</li>
		<li>Copy the unique structures list (uniqueStructures-sb.data) from the gly-h2o-n/QM directory to the gly-h2o-n/QM/pw91-lb directory and change to that directory.</li>
		<li>Run the large-basis DFT configurational sampling script with the command: 
		run-pw91-lb.csh uniqueStructures-sb.data lb QUEUE 10 
		where lb is a label for this set of calculations, QUEUE 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 the calculations. Enter &#39;test&#39; for the &#39;QUEUE&#39; to do a dry run testing.</li>
	</ol>
	</li>
	<li>Once the submitted calculations are complete, extract and analyze the data
	<ol>
		<li>Compute the rotational constants of the large-basis-optimized clusters with the command: 
		getRotConsts-dft-lb.csh pw91 N 
		where pw91 indicates that the PW91 density functional was used, and N is the number of atoms in the cluster.</li>
		<li>Now identify the unique structures with the command: 
		similarityAnalysis.py lb rotConstsData_C 
		where lb is used as a file-naming label. You now have a list of unique configurations optimized at the PW91/6-311++G** level of theory saved in the file uniqueStructures-lb.data.</li>
	</ol>
	</li>
</ol>
5. Final Energy and Thermodynamic Correction Calculations
NOTE: The goal here is to obtain the vibrational structure and energies of the Gly(H2O)n=1-5 clusters using a large basis set and an ultrafine integration grid in order to compute the desired thermochemical corrections.
<ol>
	<li>Starting with results from the previous step, submit more reliable calculations.
	<ol>
		<li>Create a subdirectory called ultrafine under the QM/pw91-lb directory. Then copy the unique structures list (uniqueStructures-lb.data) from the QM/pw91-lb directory to the QM/pw91-lb/ultrafine directory and change to that directory.</li>
		<li>Submit the ultrafine large-basis DFT script with the command: 
		run-pw91-lb-ultrafine.csh uniqueStructures-lb.data uf QUEUE 10 
		where uf is a label for this set of calculations, QUEUE 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 the calculations. Enter &#39;test&#39; for the &#39;QUEUE&#39; to do a dry run testing.</li>
	</ol>
	</li>
	<li>Once the submitted calculations are complete, extract and analyze the data
	<ol>
		<li>Extract the energies and compute the rotational constants of the large-basis-optimized clusters with the command: 
		getRotConsts-dft-lb-ultrafine.csh pw91 N 
		where pw91 indicates that the PW91 density functional was used, and N is the number of atoms in the cluster.</li>
		<li>Now identify the unique structures with the command: 
		similarityAnalysis.py uf rotConstsData_C 
		where uf is used as a file-naming label. You now have a list of unique configurations optimized at the PW91/6-311++G** level of theory saved in the file uniqueStructures-uf.data.</li>
	</ol>
	</li>
	<li>Perform a final extraction of information needed to calculate thermodynamic corrections. Use that information to compute the thermodynamic corrections.
	<ol>
		<li>Extract the final electronic energies, rotational constants and vibrational frequencies, and use them to calculate thermodynamic corrections using the command: 
		run-thermo-pw91.csh uniqueStructures-uf.data</li>
		<li>Copy/paste the command-line output to the &#39;Raw_Energies&#39; sheet of the Excel spreadsheet named &#39;gly-h2o-n.xlsx&#39;. You would need to do this for the monomers (glycine and water) as well as the lowest energy member of each hydrate (gly-h2o-n, where n=1,2, &#8230;).</li>
		<li>As the raw energies are added to the first sheet of the &#39;gly-h2o-n.xlsx&#39; spreadsheet, the subsequent &#39;Binding_Energies&#39; and &#39;Hydrate_Distribution&#39; sheets are automatically updated. In particular, the &#39;Hydrate_Distribution&#39; sheet yields the equilibrium concentration of hydrates at different temperatures (Eg. 298.15K), relative humidity (20%, 50%, 100%) and initial concentrations of water ([H2O]) and glycine ([Glycine]). The theory behind these calculations is described in the next step.</li>
	</ol>
	</li>
</ol>
6. Computing atmospheric concentrations of Gly(H2O)n=0-5 clusters at room temperature at sea-level
NOTE: This is accomplished by first copying the thermodynamic data generated in the previous step into a spreadsheet and calculating the Gibbs free energies of sequential hydration. Then, the Gibbs free energies are used to calculate equilibrium constants for each sequential hydration. Finally, a set of linear equations are solved to get the equilibrium concentration of the hydrates for a given concentration of monomers, temperature and pressure.
<ol>
	<li>Start by setting up a system of chemical equilibria for the sequential hydration of glycine as shown below: 
	<img alt="Equation 1" src="/files/ftp_upload/60964/60964equ01.jpg" /></li>
	<li>Compute the equilibrium constants Kn using Kn = e-&#916;Gn/(kBT), where n is the level of hydration, &#916;Gn is the Gibbs free energy change of the nth hydration reaction, kB is Boltzmann&#39;s constant, and T is temperature. 
	<img alt="Equation 2" src="/files/ftp_upload/60964/60964equ02.jpg" /></li>
	<li>Set up the equation for the conservation of mass, using the assumption that the sum of the equilibrium concentrations of the hydrated and un-hydrated glycine clusters equals the initial concentration of isolated glycine [Gly]0. Rewrite this system of six simultaneous equations, using some algebraic rearrangement of the equilibrium constant expressions, as 
	<img alt="Equation 3" src="/files/ftp_upload/60964/60964equ03.jpg" /></li>
	<li>Solve the system of equations shown above to obtain the equilibrium concentrations of Gly(H2O)n = 0-5 using an experimental value56,57,58 for the concentration of glycine in the atmosphere, [Gly]0 = 2.9 x 106 cm-3, and the concentration of water in the atmosphere at 100% relative humidity and a temperature of 298.15 K59, [H2O] = 7.7 x 1017 cm-3.</li>
</ol>

The accuracy of the data generated by this protocol depends mainly on three things: (i) the variety of configurations sampled by Step 2, (ii) the accuracy of the electronic structure of the system, (iii) and the accuracy of the thermodynamic corrections. Each of these factors can be addressed by modifying the method by editing the included scripts. The first factor is easily overcome with the use of a larger initial pool of randomly generated structures, more numerous iterations of the GA, and a looser definition of the criteria involved in the GA. In addition, one may use a different semi-empirical method such as the self-consistent charge density-functional tight-binding (SCC-DFTB)62 model and the effective fragment potential (EFP)63 model in order to explore the effects of different physical descriptions. The main limitation here is the inability of the method to form or break covalent bonds, meaning that the monomers are frozen. The GA procedure only finds the most stable relative positions of these frozen monomers according to the semi-empirical description.
The accuracy of the electronic structure of the system can be improved in a variety of ways, each with its computational cost. One may choose a better density functional, such as M06-2X64 and wB97X-V65, or quantum mechanical (QM) method such as the M&#248;ller-Plesset66,67,68 (MPn) perturbation theories and coupled-cluster69 (CC) methods in order to improve the physical description of the system. In the hierarchy of functionals, the performance generally improves upon going from generalized-gradient approximation (GGA) functionals like PW91 to range-separated hybrid functionals like wB97X-D and meta-GGA hybrid functionals like M06-2X.
The disadvantage of DFT methods is that a systematic convergence towards an accurate value is not possible; however, DFT methods are computationally inexpensive and there is a wide variety of functionals for a wide variety of applications.
Energies calculated using wavefunction methods like MP2 and CCSD(T) in conjunction with correlation consistent basis sets of increasing cardinal number ([aug-]cc-pV[D,T,Q,...]Z) converge towards their complete basis set limit systematically, but the computational cost of each calculation becomes prohibitive as the system size grows. Further refinement of the electronic structure can be accomplished by using explicitly correlated basis sets70 and by extrapolating to the complete basis set (CBS)71 limit. Our recent work suggests that a density-fitted explicitly correlated second-order M&#248;ller-Plesset (DF-MP2-F12) perturbative approach yields energies approaching that of MP2/CBS computations32. Modification of the current protocol to use different electronic structure methods involves two steps: (i) prepare a template input file following the syntax given by the software, (ii) and edit the run-pw91-sb.csh, run-pw91-lb.csh, and run-pw91-lb-ultrafine.csh scripts to generate the correct input file syntax as well as the correct submit script for the software.
Lastly, the accuracy of the thermodynamic corrections depends on the electronic structure method as well as the description of the PES around the global minimum. An accurate description of the PES requires the computation of third- and higher-order derivatives of the PES with respect to displacements in the nuclear degrees of freedom, such as the quartic force field72,73 (QFF), which is an exceptionally costly task. The current protocol uses the harmonic oscillator approximation to the vibrational frequencies, resulting in the need to compute only up to second derivatives of the PES. This approach becomes problematic in systems with high anharmonicity, such as very floppy molecules and symmetric double-well potentials due to the large difference in the true PES and the harmonic PES. Furthermore, the cost of having a high-quality PES from a computationally demanding electronic structure method only compounds the problem of cost for vibrational frequency calculations. One approach to overcome this is to use the electronic energies from a high-quality electronic structure calculation along with vibrational frequencies computed on a lower quality PES, resulting in a balance between cost and accuracy. The current protocol can be modified to use different PES descriptions as described in the previous paragraph; however, one may also edit the vibrational frequency keywords in the scripts and templates to compute anharmonic vibrational frequencies.
Two crucial issues for any configurational sampling protocol are the initial method for sampling the potential energy surface and the criteria used to identify each cluster. We have made extensive use of a variety of methods in our previous work. For the first issue, the initial method for sampling the potential energy surface, we have made the choice of using GA with semi-empirical methods based on these factors. Configurational sampling using chemical intuition26, random sampling, and molecular dynamics (MD)29,30, fail to find putative global minima regularly for clusters larger than 10 monomers, as we observed in our studies of water clusters18. We have successfully used basin hopping (BH) to study the complex PES of (H2O)1174, but it required the manual inclusion of some potential low energy isomers the BH algorithm did not find. A comparison of the performance of BH and GA in finding the global minimum of water clusters, (H2O)n=10-20 demonstrated that GA consistently found the global minimum faster than BH75. GA as implemented in OGOLEM and CLUSTER is very versatile because it can be applied to any molecular cluster and it can interface with a vast number of packages with classical force field, semi-empirical, density functional, and ab initio capabilities. The choice of PM7 is driven by its speed and reasonable accuracy. Virtually any other semi-empirical method would have significantly higher computational cost.
As for the second issue, we have explored using different criteria to identify unique structures ranging from electronic energies, dipole moments, overlap RMSDs and rotational constants. Using dipole moments proved difficult because both the dipole moment components were dependent on the molecule&#39;s orientation and the total dipole moment was very sensitive to geometry differences in such a way that it was difficult to set thresholds determining is structures are the same or unique. A combination of electronic energies and rotational constants proved to be most useful.
The current criteria for deeming two structures unique is based on an energy difference threshold of 0.10 kcal mol-1 and rotational constant difference of 1%. Therefore, two structures are considered different if their energies differ by more than 0.10 kcal mol-1 (~0.00015 a.u.) AND any of their three rotational constants (A, B, C) differ by more than 1%. Substantial internal benchmarks over the years found these thresholds to be reasonable choices. Our configurational sampling approach and screening methodology has been applied to very weakly bound clusters such as polyaromatic hydrocarbons complexed with water76,77 as well as strongly bound ternary sulfate hydrates containing ammonia and amines32. For clusters where there are different protonation states to be considered, the best approach is to run various GA calculations, each starting with monomers in different protonation states. This ensures that structures with different protonation states are carefully considered. However, the low-level DFT calculations often allow protonation states to change during the course of the geometry optimization, thereby yielding the most stable protonation state regardless of the starting geometry.
Our GA configurational sampling methods should work well even for floppy molecules as long as the GA codes are interfaced with general, non-parameterized methods that allow the monomers to adopt different configurations during the course of the GA run. For example, interfacing GA with PM7 would allow monomers&#39; structures to change, but if their bonds break as would happen when protonation states change, the structures may get discarded as unacceptable candidates.
We have considered different ways of correcting the shortcomings of the harmonic approximation, especially those arising from low vibrational frequencies. Incorporating the quasi-harmonic approximation into the current methodology is not difficult. However, there are still questions about the quasi-harmonic method, especially when it comes to the cutoff frequency below which it will be applied. Also, there are no rigorous benchmarking works examining the reliability of the quasi-RRHO approximation even though conventional wisdom suggests it should be an improvement over RRHO approximation.
The protocol thus presented may be generalized to any system of noncovalently-bound gas phase molecular clusters. It may also be generalized to use any semi-empirical method, electronic structure method and software, and vibrational analysis method and software by editing the scripts and templates. This assumes that the user is comfortable with the Linux command-line interface, Python scripting, and high-performance computing. The unfamiliar syntax and look of the Linux operating system and lack of scripting experience is the largest pitfall in this protocol and is where new students struggle the most. This protocol has been used successfully in a variety of implementations for years in our group, mostly focusing on the effects of sulfuric acid and ammonia on aerosol formation. Further improvements to this protocol will involve a more robust interface to more electronic structure software, alternative implementations of the genetic algorithm, and possibly the use of newer methods for faster computations of electronic and vibrational energies. Our current applications of this protocol are exploring the importance of amino acids in the early stages of aerosol formation in the current atmosphere and in the formation of larger biological molecules in prebiotic environments.

The most uncertain parameter in atmospheric studies of climate change is the exact extent to which cloud particles reflect incoming solar radiation. Aerosols, which are particulate matter suspended in a gas, form cloud particles called cloud condensation nuclei (CCN) that scatter incoming radiation, thus preventing its absorption and the subsequent heating of the atmosphere1. A detailed understanding of this net cooling effect requires an understanding of the growth of aerosols into CCNs, which in turn requires an understanding of the growth of small molecular clusters into aerosol particles. Recent work has suggested that aerosol formation is initiated by molecular clusters of 3 nm in diameter or less2; however, this size regime is difficult to access using experimental techniques3,4. Therefore, a computational modeling approach is desired in order to overcome this experimental limitation.
Using our modeling approach described below, we can analyze the growth of any hydrated cluster. Because we are interested in the role of water in the formation of large biological molecules from smaller constituents in pre-biotic environments, we illustrate our approach with glycine. The challenges encountered and tools needed to address those research questions are very similar to those involved in the study of atmospheric aerosols and prenucleation clusters5,6,7,8,9,10,11,12,13,14,15. Here, we examine hydrated glycine clusters starting from an isolated glycine molecule followed by a series of stepwise additions of up to five water molecules. The final goal is to calculate the equilibrium concentrations of Gly(H2O)n=0-5 clusters in the atmosphere at room temperature at sea-level and a relative humidity (RH) of 100 %.
A small number of these sub-nanometer molecular clusters grow into a metastable critical cluster (1-3 nm in diameter) either by adding other vapor molecules or coagulating on existing clusters. These critical clusters have a favorable growth profile leading to the formation of much larger (up to 50-100 nm) cloud condensation nuclei (CCN), which directly affect the precipitation efficiency of clouds as well their ability to reflect incident light. Therefore, having a good understanding of the thermodynamics of molecular clusters and their equilibrium distributions should lead to more accurate predictions of the impact of aerosols on the global climate.
A descriptive model of aerosol formation requires accurate thermodynamics of molecular cluster formation. The computation of accurate thermodynamics of molecular cluster formation requires the identification of the most stable configurations, which involves finding the global and local minima on the cluster&#39;s potential energy surface (PES)16. This process is called configurational sampling and can be achieved through a variety of techniques, including those based on molecular dynamics (MD)17,18,19,20, Monte Carlo (MC)21,22, and genetic algorithms (GA)23,24,25.
Different protocols have been developed over the years to obtain the structure and thermodynamics of atmospheric hydrates at a high level of theory. These protocols differed in the choice of (i) configurational sampling method, (ii) nature of low-level method used in the configurational sampling, and (iii) the hierarchy of higher-level methods used to refine the results in the subsequent steps.
The configurational sampling methods included chemical intuition26, random sampling27,28, molecular dynamics (MD)29,30, basin hopping (BH)31, and genetic algorithm (GA)24,25,32. The most common low-level methods employed with these sampling methods are force fields or semi-empirical models such as PM6, PM7 and SCC-DFTB. These are often followed by DFT calculations with increasingly larger basis sets and more reliable functionals from the higher rungs of Jacob&#39;s ladder33. In some cases, these are followed by higher level wavefunction methods such as MP2, CCSD(T), and the cost efficient DLPNO-CCSD(T)34,35.
Kildgaard et al.36 developed a systematic method where water molecules are added at points on the Fibonacci spheres37 around smaller hydrated or unhydrated clusters to generate candidates for larger clusters. Unphysical and redundant candidates are removed based on close contact thresholds and root-mean-square distance between different conformers. Subsequent optimizations using the PM6 semi-empirical method and a hierarchy of DFT and wavefunction methods are used to get a set of low energy conformers at a high level of theory.
The artificial bee colony (ABC) algorithm38 is a new configurational sampling approach that has recently been implemented by Zhang et al. to study molecular clusters in a program called ABCluster39. Kubecka et al.40 used ABCluster for configurational sampling followed by low-level reoptimizations using the tight-binding GFN-xTB semi-empirical method41. They further refined the structures and energies using DFT methods followed by final energies using DLPNO-CCSD(T).
Regardless of the method, configurational sampling starts with a randomly- or nonrandomly-generated distribution of points on the PES. Each point corresponds to a specific geometry of the molecular cluster in question and is generated by the sampling method. Then the closest local minimum is found for each point by following the &#34;downhill&#34; direction on the PES. The set of minima thus found correspond to those geometries of the molecular cluster that are stable, at least for some time. Here, the shape of the PES and the evaluation of the energy at each point on the surface will be sensitive to the physical description of the system where a more accurate physical description results in a more computationally expensive energy calculation. We will specifically use the GA method implemented in the OGOLEM25 program, which has been successfully applied to a variety of global optimization and configurational sampling problems42,43,44,45, to generate the initial set of sampling points. The PES will be described by the PM7 model46 implemented in the MOPAC2016 program47. This combination is employed because it generates a larger variety of points compared to the MD and MC methods and finds the local minima faster than more-detailed descriptions of the PES.
The set of GA-optimized local minima are taken as the starting geometries for a series of screening steps, which lead to a set of low-lying minimum energy. This part of the protocol begins by optimizing the set of unique GA-optimized structures using density-functional theory (DFT) with a small basis set. This set of optimizations will generally give a smaller set of unique local minimum structures which are modeled in more detail compared to the GA-optimized semi-empirical structures. Then another round of DFT optimizations are performed on this smaller set of structures using a larger basis set. Again, this step will generally give a smaller set of unique structures which are modeled in more detail compared to the small basis DFT step. The final set of unique structures are then optimized to a tighter convergence and the harmonic vibrational frequencies are calculated. After this step we have everything we need to compute the equilibrium concentrations of the clusters in the atmosphere. The overall approach is summarized diagrammatically in Figure 1. We will use the PW9148 generalized-gradient approximation (GGA) exchange-correlation functional in the Gaussian0949 implementation of DFT along with two variations of the Pople50 basis set (6-31+G* for the small basis step and 6-311++G** for the large basis step). This particular combination of exchange-correlation functional and basis set was chosen due to its previous success in computing accurate Gibbs free energies of formation for atmospheric clusters51,52.
This protocol assumes that the user has access to a high-performance computing (HPC) cluster with the portable batch system53 (PBS), MOPAC2016 (http://openmopac.net/MOPAC2016.html)47, OGOLEM (https://www.ogolem.org)25, Gaussian 09 (https://gaussian.com)49, and OpenBabel54 (http://openbabel.org/wiki/Main_Page) software installed following their specific installation instructions. Each step in this protocol also uses a set of in-house shell and Python 2.7 scripts which must be saved to a directory that is included in the user&#39;s $PATH environmental variable. All necessary environmental modules and execution permissions to run all of the above programs must also be loaded into the user&#39;s session. The disk and memory usage by the GA code (OGOLEM) and semi-empirical codes (MOPAC) are very small by modern computer resource standards. The overall memory and disk usage for OGOLEM/MOPAC depends on how many threads one wants to use, and even then, the resource usage will be small compared to the capabilities of most HPC systems. The resource needs of the QM methods depend on the size of the clusters and the level of theory used. The advantage of using this protocol is that one can vary the level of theory to be able to calculate the final set of low energy structures, keeping in mind that usually faster calculations lead to more uncertainty in accuracy of the results.
For the sake of clarity, the user&#39;s local computer will be referred to as &#34;local computer&#34; while the HPC cluster they have access to will be referred to as &#34;remote cluster&#34;.

<table border="0" cellpadding="0" cellspacing="0" style="width:1639px;" width="1638">
	
	<tbody>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">Avogadro</td>
			<td style="width:483px;">https://avogadro.cc</td>
			<td style="width:123px;"></td>
			<td style="width:617px;">Open-source molecular visualization program</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">Gaussian [09/16] Software</td>
			<td style="width:483px;">http://www.gaussian.com/</td>
			<td style="width:123px;"></td>
			<td>Commercial ab initio electronic structure program</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">MOPAC 2016</td>
			<td style="width:483px;">http://openmopac.net/MOPAC2016.html</td>
			<td style="width:123px;"></td>
			<td>Open-source semi-empirical program</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">OGOLEM Software</td>
			<td style="width:483px;">https://www.ogolem.org</td>
			<td style="width:123px;"></td>
			<td>Genetic algorithm-based global optimization program</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">OpenBabel</td>
			<td style="width:483px;">http://openbabel.org/wiki/Main_Page</td>
			<td style="width:123px;"></td>
			<td>Open-source cheminformatics library</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">calcRotConsts.py</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Python script to compute rotational constants</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">calcSymmetry.csh</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Shell script to calculate symmetry number of a molecule given Cartesian coordinates</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">combine-GA.csh</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Shell script to combine energy and rotational constants from different GA directories</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">combine-QM.csh</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Shell script to combine energy and rotational constants from different QM directories</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">gaussianE.csh</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Shell script to extract Gaussian 09 energies</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">gaussianFreqs.csh</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Shell script to extract Gaussian 09 vibrational frequencies</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">getrotconsts</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Executable to calculate rotational constants given a molecule&#39;s Cartesian coordinates</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">getRotConsts-dft-lb.csh</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Shell script to compute rotational constants for a batch of large basis DFT optimized structures</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">getRotConsts-dft-lb-ultrafine.csh</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Shell script to compute rotational constants for a batch of ultrafine DFT optimized structures</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">getRotConsts-dft-sb.csh</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Shell script to compute rotational constants for a batch of small basis DFT optimized structures</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">getRotConsts-GA.csh</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Shell script to compute rotational constants for a batch of genetic algorithm optimized structures</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">global-minimum-coords.xyz</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Cartesian coordinates of global minimum structures of gly-(h2o)n, where n=0-5</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">make-thermo-gaussian.csh</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Shell script to extract data from Gaussian output files and make input files for the thermo.pl script</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">ogolem-input-file.ogo</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Ogolem sample input file</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">ogolem-submit-script.pbs</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>PBS batch submission file for Ogolem calculations</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">README.docx</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Clarifications to help readers use the scripts effectively</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">runogolem.csh</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Shell script to run OGOLEM</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">run-pw91-lb.csh</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Shell script to run a batch of large basis DFT optimization calculations</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">run-pw91-lb-ultrafine.csh</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Shell script to run a batch of ultrafine DFT optimization calculations</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">run-pw91-sb.csh</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Shell script to run a batch of small basis DFT optimization calculations</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">run-thermo-pw91.csh</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Shell script to compute the thermodynamic corrections for a batch of DFT optimized structures</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">similarityAnalysis.py</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Python script to determine unique structures based on rotational constants and energies</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">symmetry</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Executable to calculate molecular symmetry given Cartesian coordinates</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">symmetry.c</td>
			<td style="width:483px;">(C) 1996, 2003 S. Patchkovskii, Serguei.Patchkovskii@sympatico.ca</td>
			<td style="width:123px;"></td>
			<td style="width:617px;">C code to determine the molecular symmstry of a molecule given Cartesian coordinates</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">template-marcy.pbs</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Template for a PBS submit script which uses OGOLEM</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">template-pw91.com</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Template Gaussian 09 input</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">template-pw91-HL.com</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Template Gaussian 09 input for ultrafine DFT optimization</td>
		</tr>
		<tr height="42">
			<td height="42" style="height:42px;width:416px;">thermo.pl</td>
			<td style="width:483px;">https://www.nist.gov/mml/csd/chemical-informatics-research-group/products-and-services/program-computing-ideal-gas</td>
			<td style="width:123px;"></td>
			<td>Perl open-source script to compute ideal gas thermodynamic corrections</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">gly-h2o-n.xlsx</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Excel spreadsheet for the complete protocol</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">table-1.xlsx</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Excel spreadsheet</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">table-2.xlsx</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Excel spreadsheet</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">table-3.xlsx</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Excel spreadsheet</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">table-4.xlsx</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Excel spreadsheet</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">water.xyz</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Cartesian coordinates of water</td>
		</tr>
		<tr height="21">
			<td height="21" style="height:21px;width:416px;">glycine.xyz</td>
			<td style="width:483px;">Shields Group, Department of Chemistry, Furman University</td>
			<td style="width:123px;"></td>
			<td>Cartesian coordinates of glycine</td>
		</tr>
	</tbody>
</table>

computation of atmospheric concentrations of molecular clusters from ab initio thermochemistry

The computational study of the formation and growth of atmospheric aerosols requires an accurate Gibbs free energy surface, which can be obtained from gas phase electronic structure and vibrational frequency calculations. These quantities are valid for those atmospheric clusters whose geometries correspond to a minimum on their potential energy surfaces. The Gibbs free energy of the minimum energy structure can be used to predict atmospheric concentrations of the cluster under a variety of conditions such as temperature and pressure. We present a computationally inexpensive procedure built on a genetic algorithm-based configurational sampling followed by a series of increasingly accurate screening calculations. The procedure starts by generating and evolving the geometries of a large set of configurations using semi-empirical models then refines the resulting unique structures at a series of high-level ab initio levels of theory. Finally, thermodynamic corrections are computed for the resulting set of minimum-energy structures and used to compute the Gibbs free energies of formation, equilibrium constants, and atmospheric concentrations. We present the application of this procedure to the study of hydrated glycine clusters under ambient conditions.

Watch this Scientific Journal Video about Computation of Atmospheric Concentrations of Molecular Clusters from ab initio Thermochemistry at JoVE.com

Computation of Atmospheric Concentrations of Molecular Clusters from ab initio Thermochemistry

The atmospheric concentrations of weakly bound molecular clusters can be computed from the thermochemical properties of low energy structures found through a multi-step configurational sampling methodology utilizing a genetic algorithm and semi-empirical and ab initio quantum chemistry.

Computation of Atmospheric Concentrations of Molecular Clusters from ab initio Thermochemistry

The computational study of the formation and growth of atmospheric aerosols requires an accurate Gibbs free energy surface, which can be obtained from gas ...

Education

Research

JoVE Core

JoVE Journal

Chemistry

Thermochemistry

Chemical Equilibrium

SCIENTISTS IN ACTION

Isolation and Chemical Characterization of Lipid A from Gram-negative Bacteria

Lipopolysaccharide (LPS) is the major cell surface molecule of gram-negative bacteria, deposited on the outer leaflet of the outer membrane bilayer. LPS ...

A Simple Method for the Size Controlled Synthesis of Stable Oligomeric Clusters of Gold Nanoparticles under Ambient Conditions

Reducing dilute aqueous HAuCl4 with sodium thiocyanate (NaSCN) under alkaline conditions produces 2 to 3 nm diameter nanoparticles. Stable grape-like ...

Construction and Systematical Symmetric Studies of a Series of Supramolecular Clusters with Binary or Ternary Ammonium Triphenylacetates

Functions of clusters in nano or sub-nano scale significantly depend on not only kinds of their components but also arrangements, or symmetry, of their ...

Improved Heterojunction Quality in Cu2O-based Solar Cells Through the Optimization of Atmospheric Pressure Spatial Atomic Layer Deposited Zn1-xMgxO

Improved Heterojunction Quality in Cu2O-based Solar Cells Through the Optimization of Atmospheric Pressure Spatial Atomic Layer Deposited Zn1-xMgxO

Atmospheric pressure spatial atomic layer deposition (AP-SALD) was used to deposit n-type ZnO and Zn1-xMgxO thin films onto p-type thermally oxidized Cu2O ...

Analyzing the Photo-oxidation of 2-propanol at Indoor Air Level Concentrations Using Field Asymmetric Ion Mobility Spectrometry

We demonstrate a versatile protocol to be used for determining the effectiveness of photocatalysts in degrading indoor air concentration (ppb) volatile ...

Direct Imaging of Laser-driven Ultrafast Molecular Rotation

We present a method for visualizing laser-induced, ultrafast molecular rotational wave packet dynamics. We have developed a new 2-dimensional Coulomb ...

Fizzy Extraction of Volatile Organic Compounds Combined with Atmospheric Pressure Chemical Ionization Quadrupole Mass Spectrometry

Chemical analysis of volatile and semivolatile compounds dissolved in liquid samples can be challenging. The dissolved components need to be brought to ...

Combustion Chemistry of Fuels: Quantitative Speciation Data Obtained from an Atmospheric High-temperature Flow Reactor with Coupled Molecular-beam Mass Spectrometer

This manuscript describes a high-temperature flow reactor experiment coupled to the powerful molecular beam mass spectrometry (MBMS) technique. This ...

Light-driven Molecular Motors on Surfaces for Single Molecular Imaging

The design and synthesis of a synthetic system that aims for the direct visualization of a synthetic rotary motor at the single molecule level on surfaces ...

Synthesis of In37P20(O2CR)51 Clusters and Their Conversion to InP Quantum Dots

Synthesis of In37P20O2CR51 Clusters and Their Conversion to InP Quantum Dots

This text presents a method for the synthesis of In37P20(O2C14H27)51 clusters and their conversion to indium phosphide quantum dots. The In37P20(O2CR)51 ...