The Topology Sampling of H2SO4·NH3 with Meta-Dynamics Method

Yugai Huang,Yu Shang,Wenqing LiJiajing HeHairong Zhou and Bin Gu

Abstract:The configurations of molecular clusters have significant impacts on their growth into fine particles in atmosphere.In this paper,we explore the topology space of the structure of H2SO4·NH3 dimer with a novel sampling technique of meta-dynamics(MTD)method and ab initio molecular dynamics simulations.The simulations are carried out at the temperatures of both 50 K and 242 K,which represent the typical high and low latitudes of troposphere.The results show that,compared with only traditional MD simulations,the structure samplings are significantly accelerated with MTD method.Therefore,more isomers of the dimer are discovered within the same simulation time scale.In addition,the results show that MTD is more efficient for circumstances with high temperature.

Keywords:Structure sampling,H2SO4 and NH3,meta-dynamics acceleration,temperature effects.

1 Introduction

It is well known that fine aerosol particles in the atmosphere have very important impacts on the weather and climate.In addition,they also affect human health and safety[Andreae (2013);Lohmann and Feichter (2005)] significantly,especially in urban area with many industrial and anthropogenic emissions.These fine aerosol particles have been proved to be closely related to air pollution events [Guo,Hu,Zamora et al.(2014);Zhang,Wang,Luo et al.(2015)].Observations have demonstrated that sulphuric acid,as the principal secondary product of sulphide emissions,is one of the most important sources of aerosol particles in the troposphere [Merikanto,Zapadinsky,Lauri et al.(2007);Kulmala (2003)].Other organic and inorganic emissions,such as ammonia,amines and amino acids,might play as catalysis for the growth of aerosol particles from nanometersize to ultra-fine seeds of cloud condensation nuclei (CCN).The knowledges of generation and growth of the ultra-fine clusters from different source molecules and ions in the atmosphere are of interdisciplinary interests [Merikanto,Zapadinsky,Lauri et al.(2007);Kulmala (2003);Zhang,Khalizov,Wang et al.(2012)].But there are still great challenges to observe directly the micro-structure changes of the small clusters and their growth in both field and laboratory.Thus,molecular simulations,based on classical force fields and first principle theories,have been widely adopted in the studies of nucleation processes at the atomistic level of view.For all molecular simulation studies,the first step is often to design the clustering structures of the molecules.But in natural environments with thermal fluctuation and atmospheric disturbance,there always exist a huge number of isomer configurations for clusters with even very simple components.In practices,some reasonable configurations are designed based on chemical intuitions [Nadykto,Herb,Yu et al.(2014)].The traditional sampling techniques,such as Monte Carlo [Xu and Zhang (2013);Shanker and Bandyopadhyay (2011);Zhan,Chen and Liu (2006)],are very efficient for local minimization.But in dynamics simulations,they are always not efficient because of no guarantee on ergodicity.For more detailed studies on the generation and evolution of atmospheric particles,an efficient sampling method is required to figure out most of the possible clustering isomers.

The innovative technique of meta-dynamics (MTD)[Laio and Gervasio (2008);Fiorin,Klein and Henin (2013)] has been verified to be quite efficient in speeding up the map searching with multiple local minimizations.Its core idea is filling up local potential wells for all possible reactions through randomly adding small bias potentials (VG)along the searching trajectory.Theoretically,for any specific sampling,the ergodicity can be promised with small enough VG.MTD has been applied widely in the studies of biomolecules structural transitions,phase transitions and even combustion reactions [Laio and Gervasio (2008);Fiorin,Klein and Henin (2013);Zheng and Pfaendtner (2014)].

The H2SO4·NH3dimer,as nascent seed of CCN,can be activated to grow into bigger clusters.As the hydrogen-bonds in the H2SO4·NH3dimer are sensitive to environmental conditions,the structure of H2SO4·NH3dimer might change evidently at different latitude.In this paper,to exploring the configurations of H2SO4·NH3dimer in troposphere,both the traditional AIMD and the MTD accelerated AIMD (MTD-AIMD)simulations were carried out.The initial structure of H2SO4·NH3dimer is shown in Fig.1.The detailed simulation theories and methods are given in section 2.The results along with the discussions on the structure sampling and energy calculations are presented in section 3.The conclusions and further applications of the MTD-AIMD simulation procedure in the studies of aerosol particles are given at the end of the paper.

Figure1:The initial configuration of the H2SO4·NH3 dimer with one hydrogen bond(dash line).a)NH3 and b)H2SO4 (color online)

2 Model and simulation methods

In the application of the general MTD technique,a reaction or transition space should be defined with a set of reasonable collective variables (CVs,noted asS).In chemical structure sampling,the general CVs should be independent of the arrangement of all atoms.Using the switch function：

whereRijis the distance between any atom pair ofiandj.R0is the standard bonding criteria,to describe the bond strength of the pair of atoms (0 ≤aij≤1).All topology information of the cluster can be found in the symmetric contact matrix ofaij.Setting the largest eigenvalue of the matrix asand its corresponding eigenvector asa social permutation invariant (SPRINT)coordinates of a cluster can be defined as[Pietrucci and Andreoni (2011)]：

Here,“sorted” indicates the components of the eigenvector are sorted in ascending order.The SPRINT coordinates are invariant for the permutation of all atoms in a species.They represent the bonding complexity of each atomic species in the cluster [Pietrucci and Andreoni (2011)].In the subspace of CVs,a time-dependent bias potential (VG)as the summation of randomly placed small Gaussian hills along the simulation trajectory is defined as：

wheretis the simulation time,ωis the energy rate of the Gaussian hill,siis the hill width in each dimension of the subspace,anddis the total number of the CVs.With the randomly placed bias-potential,all the local barriers in the reaction space should be overcome sooner or later.As a result the sampling efficiency can be dramatically increased.Furthermore,the free energy profile of the reaction can be reconstructed as a reverse of the summing of the bias-potential [Fiorin,Klein and Henin (2013)].

In this work,the structures of the H2SO4·NH3dimer are sampled in the subspace of the SPRINT coordinates of five protons.As the configuration of H2SO4·NH3dimer is determined by the hydrogen-bonding pattern,only the hydrogen related contact matrix elements are none zero.The bonding criteria was set as 2.4 Å.The height of the Gaussian potential hill was set as 2.0 kcal/mol,which is small enough compared with the hydrogen bond strength.The width of the hills was set as 1.5 a.u.The time step of the AIMD was set as 1 fs.The bias potential hills were added every 200 steps along the MD trajectory.Density functional theories (DFT)have been widely used inab initiosimulation studies on atmospheric particles [Temelso,Phan and Shields (2012);Nadykto,Yu,Jakovleva et al.(2011);Chon,Lee and Lin (2014);Elm,Bilde and Mikkelsen (2013);DePalma,Doren and Johnston (2014);Nadykto,Herb,Yu et al.(2014);Herb,Xu,Yu et al.(2013)].In addition to DFT,several higher level of theories,such as Hartree-Fock (HF),Moller-Plesset perturbation (MP),coupled cluster with singles and doubles (CCSD)et al.,have also been used in geometry optimization,spectrum calculation and reaction energy analysis of aerosol related clusters [Zhang,Khalizov,Wang et al.(2012);Zhao,Khalizov,Zhang et al.(2009);Xu and Zhang (2012)].In this paper,the all electron hybrid-DFT(PBE0)level of theory [Adamo and Barone (1999)] was employed.By partially including the Hartree-Fock exchange term,hybrid-DFT (PBE0)is better than DFT in eliminating the fake self-interactions of electrons [Perdew and Zunger (1981)].Thus,hybrid-DFT is promising in deal with ions and radicals with unpaired electrons and holes [Cora,Alfredsson,Mallia et al.(2004);Morisanchez,Cohen and Yang (2006);Sai,Barbara and Leung (2011);Schmidt and Kummel (2016)].The accuracy of the hybrid-DFT scheme was verified by our former work [Wang,Huang,Gu et al 2016].The aug-cc-pV(T+d)Z basis sets,which including diffuse orbitals,were employed for sulphur atoms [Jr,Peterson and Wilson (2001)],while the aug-cc-pVTZ basis sets were used for all other atomic species.

The simulation engine of this work was the quantum module named Quickstep (QS)of CP2K [Vandevondele,Krack,Mohamed et al.(2005)].The MTD plugin was provided by the PLUMED package [Bonomi,Branduardi,Bussi et al.(2009);Tribello,Bonomi,Branduardi et al.(2013)].The Gaussian augmented-plane waves method (GAPW)[Lippert,Hutter and Parrinello (1999)],in which the Hartree energy and potential were calculated by Fourier trans-form techniques,and the Kohn-Sham orbitals were expanded on a Gaussian basis set.The Poisson solver was the Martyna Tuckerman method[Martyna and Tucker-man (1999),which decouples periodic images of the system.The SCF convergence criterion was set to be 10-7for structure sampling.

3 Results and discussions

For the situation of high altitude,the temperature was set as 50 K.Both the AIMD simulations with and without the MTD acceleration for the dimer were carried for 40 ps.As all possible topologies were sampled within 9 ps,the first 9.4 ps of the 40 ps simulation was presented in Fig.2.As shown in the top panel of Fig.2,three kinds of dimer structures of type I,II and III were sampled with MTD method.In the traditional MD simulation,only the structures of type I and II were generated.The SPRINT coordinates of the five hydrogen atoms are shown in the middle panel of Fig.2.The potential energyEpof the system is shown in the bottom panel of Fig.2.The first emergency time of the three structures are indicated with black arrows.

Figure2:Sampled configurations of H2SO4·NH3 dimer (top panel)and the SPRINT CVs of the 5 H atoms (middle panel)in MTD-AIMD simulation;the potential energy Ep (bottom panel)of both MTD-AIMD (black line)and traditional AIMD (green line)simulations at 50 K (color online)

The SPRINT coordinatesSiare very sensitive to the bonding around the protons.For type I configuration,one hydrogen atom is combined with both one nitrogen and one oxygen atoms,forming a typical hydrogen bond.In the cases of both type II and III,there are two hydrogen bonds in the dimer.At 3.8 ps the proton transferring takes place.As a result,one new hydrogen bond was formed.MeanwhileSichanges to be around 5.5.In addition to the CVs,the structure changes can also be reflected in the fluctuation of the energyEp.In the first 3 ps,the potential energy evolution of MTD and MD are basically the same.Both the isomers of type I and II appear at the same point for MTD and MD.In the following simulations,the energy changes differently from each other.The potential energy with MTD keeps increasing as a result of the accumulation of the Gaussian potential wave packets.Hence,the dimer can escape easily from the “valley” of the potential energy surface.In this way,the type III isomer was successfully sampled at 3.86 ps.As a contrast,the traditional MD cannot escape from the local minimum until the end of the simulation.

Figure3:Sampled configurations of H2SO4·NH3 dimer (top)and the SPRINT CVs of the five H atoms (middle)in MTD-AIMD simulation;the potential energy Ep (bottom)of the first 3 ps MTD (black line)and traditional AIMD (green line)simulations at 242 K (color online)

In the troposphere,the environmental temperature increases with the altitude decreases.For low troposphere,the typical temperature of 242 K was chosen for the study of the configurations of the H2SO4·NH3dimer.The possible structures and corresponding evolutions of the SPRINT coordinates of the hydrogen atoms of the dimer are shown in the top and middle panels of Fig.3.The potential energy variation is also shown at the bottom of Fig.3.

As shown in the top panel of Fig.3,not only the type I,II and III structures,but also new transition states of α and β,are sampled in MTD-AIMD simulation.The time of the first emergency of the isomers are pointed with arrows in Fig.3.It can be seen that the generation and breaking of hydrogen bonds,and the transfer of protons are clearly recorded along the MTD trajectory.The number of hydrogen bonds keeps varying between the type I and II with higher frequency in both low than in high troposphere.As a contrast,in traditional AIMD simulation,the type I,II and III isomers are found with much longer simulation time,and the transition states are not spotted.The sampling efficiency is improved evidently with MTD acceleration.

With traditional AIMD,the type III isomer appears much later than that with MTDAIMD.The potential energy barrier between the type I and II is very small as around 0.2 eV,as the corresponding hydrogen-bonding strength is weak.But the energy barrier between (I,II)and III is relatively higher as around 0.8 eV.It makes the structure transition much frequent between type I and II,especially at high temperature.

The more evident fluctuation of the CVs at 242 K compared with at 50 K is easy to understand,as high temperature will induce intense thermal movements and vibrations of the molecules.For the same reason,the fluctuation of the potential energy at higher temperature is also more evident.

Although the fluctuation of the SPRINT coordinates at low altitude with high temperature are much more evident than those at high altitude with low temperature,the structural variation patterns in the H2SO4·NH3dimers are the same in both cases.For example,when the type II structure appears,the correspondingS4andS5is about 4.2 and 5.3 individually.When the type III structure was spotted,theS4andS5is about 4.8 and 5.3 individually.The structure information is consistent for both 50 and 242 K.Hence,the capability of the SPRINT coordinates to represent the structural transformation of the cluster keeps unified for different environments.

4 Conclusions

The micro-structures of the nanoscale aerosol clusters in the troposphere are very important for both environment and climate changes [Zhang,Khalizov,Wang et al.(2012);Vereecken,Glowacki and Pilling (2015)].Although molecular simulations have contributed substantially to the understanding of the formation of sulfur acid related aerosol particles,for further understanding of the physical and chemical reactions,specific configurations of the particles with real atmospheric conditions should be taken into account.

In this work,based onab initiomolecular dynamics (AIMD)method,the configurations of the H2SO4·NH3dimer in atmosphere with both low and high temperatures were sampled with meta-dynamics (MTD)technique and traditional AIMD simulations.It is found that MTD method is effective in overcoming local potential wells,and then sampling structures.For H2SO4·NH3,the type I,II and III structures and the transition states are founded with MTD-AIMD within reasonable simulation time scales.Comparisons between the simulations for high latitude with low temperature and for near ground with high temperature show that,MTD is more efficient for high temperature circumstance with more thermal fluctuations.

The MTD-AIMD simulations might be expected in revealing new insights in the open challenges about atmospheric clusters.We will use it in our studies of the sulfuric acid particles with water vapor and other atmospheric emissions.

Acknowledgement:This work was partially supported by the program of the Key Laboratory for Aerosol Cloud Precipitation of CMA NUIST (No.KDW1304)and the innovation project of NUIST (2017).