PAPER

View Journal | View Issue

Published on 14 February 2014. Downloaded by Arizona State University on 30/10/2014 17:34:20.

Floating orbital molecular dynamics simulations Cite this: Phys. Chem. Chem. Phys., 2014, 16, 6997

a b Eva Perlt,a Marc Bru ¨ ssel and Barbara Kirchner*

We introduce an alternative ab initio molecular dynamics simulation as a unification of Hartree–Fock molecular dynamics and the floating orbital approach. The general scheme of the floating orbital molecular dynamics method is presented. Moreover, a simple but sophisticated guess for the orbital centers Received 13th November 2013, Accepted 14th February 2014

is provided to reduce the number of electronic structure optimization steps at each molecular dynamics

DOI: 10.1039/c3cp54797c

floating orbital molecular dynamics approach with and without application of the initial guess. Finally,

www.rsc.org/pccp

a water monomer and a water dimer are simulated, and the influence of the orbital floating on certain properties like the dipole moment is investigated.

step. The conservation of total energy and angular momentum is investigated in order to validate the

1 Introduction Molecular dynamics (MD) simulations constitute an established method for the investigation of condensed matter.1–5 With the introduction of ab initio molecular dynamics (AIMD)6 the accuracy of results could be improved, and quantum chemical methods could be applied to the dynamics of fluid systems. A very popular approach in this field is the Car–Parrinello molecular dynamics method (CPMD).7 Therein, explicit electronic degrees of freedom are introduced. The integration time step and thereby the simulation time are, however, limited by the fictitious dynamics. In another recently presented approach this is circumvented by propagating the density matrix by a predictor–corrector algorithm, which is called Second Generation CPMD.8,9 Other integrators and implementations are also possible.10–12 In the electronic structure theory, antisymmetric products of molecular orbitals provide the foundation for the construction of the total wave function. The former ones in turn are generated by expansion in a finite set of atomic basis functions centered on the corresponding nuclei. The choice of the basis functions for calculation of an electronic structure is very important because it determines the quality of the resulting wave functions.13 The idea of releasing the basis functions from the nuclei in molecular structure optimizations arose in 1954 when Hurley published the first study on floating orbitals (FO).14 Shull and Ebbing applied these floating wave functions to H2 and H2+.15

a

¨r Physikalische und Theoretische Chemie, Wilhelm-Ostwald-Institut fu ¨t Leipzig, Linne´str. 2, D-04103 Leipzig, Germany Universita b ¨r Physikalische und Mulliken Center for Theoretical Chemistry, Institut fu ¨t Bonn, Beringstraße 4+6, D-53115 Bonn, Theoretische Chemie, Universita Germany. E-mail: [email protected]

This journal is © the Owner Societies 2014

Later, Frost introduced his floating spherical Gaussian orbital model. He chose the exponents and centers of the basis functions as additional variational parameters and applied them to lithium hydride as well as other molecular compounds.16–18 ´sza ´r have taken up the model recently. They could Tasi and Csa show that it is possible to approach the Hartree–Fock limit energy with floating-centered and floating-exponent s-type Gaussian basis functions.19 In 1970 Ford et al. proposed a similar approach. They also applied 1s functions together with spherical Gaussian functions, which were initially located on the bonds or at the lone pair positions.20 In 1979 and subsequent years, Huber developed the floating orbital geometry optimization (FOGO), in which he employed Hellmann–Feynman forces as well as energy gradients.21–24 Approximately one decade later in 1988, Hurley investigated the water molecule using floating Gaussian functions.25 Around this time, Helgaker ¨f provided a comprehensive study on floating orbitals and Almo in which they investigated a number of properties of several molecules.26 Furthermore, the FO approach was combined with density functional theory and applied to atomic systems by Pakiari and Mohajeri in 2002.27 Another method concerning the displacements of basis functions apart from atom-centered or bond-centered ones deals with distributed Gaussian basis functions.28,29 In that approach, the displacement of the basis functions is dependent on a parametrized function which may, for example, be deduced from an anharmonic potential.30 Because of the bond formation in molecular systems, the electronic shell may no longer be focused on the nuclei. Hence, this is the physical reason for allowing the positions of the basis functions to deviate from the nuclear positions. Furthermore, orbital floating might be advantageous if an external field is ¨f demonstrated the improvement of applied. Helgaker and Almlo the accuracy of the results for properties like dipole moments, static polarizabilities, and infrared intensities.26 Finally, a possible

Phys. Chem. Chem. Phys., 2014, 16, 6997--7005 | 6997

View Article Online

Published on 14 February 2014. Downloaded by Arizona State University on 30/10/2014 17:34:20.

Paper

reduction of the basis set size was hinted at whereas the loss of accuracy is compensated by floating the orbitals. Consequently, higher angular momentum functions might be neglected. Another idea was the expression of higher angular momentum functions in terms of linear combinations of Gaussian functions that are displaced from their respective nucleus.31 A further problem which requires basis functions to be located away from the nuclei is the solvated electron. Marsalek et al. investigated the solvated electron in bulk water with AIMD using a grid of fixed Gaussians, in addition to the conventional atom-centered basis functions and plane waves.32 The positions of the orbitals and their optimization can be treated in diﬀerent ways. One way is to introduce—figuratively spoken—two diﬀerent kinds of atoms: bare nuclei which possess a charge and mass but no basis functions and so-called ghost atoms which have neither nuclear charge nor mass but carry basis functions, as examined by Hurley in his later study in 1988.25 In such an approach, the number of ghost atoms amounts to at least the number of bare nuclei, but it may also be larger, so that the separation of the atoms’ shells is possible. This method allows us to optimize the centers of the basis functions and the molecular geometry simultaneously by applying the usual optimization scheme. Another method employs two diﬀerent optimization schemes, alternatingly: the molecular geometry including nuclei and basis functions is optimized according to the energy gradient. The Hellmann–Feynman force serves to move only the nuclei while the basis functions are fixed. This approach has been successfully applied by Huber in his studies on hydrogen ion clusters.21–24 The great advantage of this method is the easy evaluation of the Hellmann–Feynman forces and the comparably fast performance. For the geometry optimization procedure, there are diﬀerent approaches. Huber applied the Murtagh–Sargent variant of the quasi-Newton algorithm and calculated the gradient explicitly ¨f whereas the Hessian was approximated.33 Helgaker and Almlo calculated both the gradient and the Hessian analytically using the trust-region method.26 In this article, the employment of floating orbitals is extended to molecular dynamics instead of only applying it to the calculation of the electronic energy or the geometry optimization. The displacement of the atoms is conducted based on conventional energy gradients applied to the conventional equations of motion in molecular dynamics (MD) simulations. Within this study, a very simple scheme is applied for the molecular dynamics procedure, i.e., the forces are evaluated from the Hartree–Fock (HF) gradient. After each MD step the centers of the basis functions are optimized according to the derivative of the total energy with respect to the basis function centers. As an optimization algorithm the Broyden–Fletcher–Goldfarb–Shanno algorithm (BFGS) is applied in conjunction with a line search scheme.34–37 However, the additional computational cost for the orbital optimization has to be considered. Therefore, the main goal is the investigation of the basis function displacements with regard to their respective nuclei along the trajectory for covalent as

6998 | Phys. Chem. Chem. Phys., 2014, 16, 6997--7005

PCCP

well as hydrogen bonded substances rather than the increase of performance. The article is organized as follows: in the first part the developed methods are described. The implementation of the quantum chemical methods are detailed, and the optimization scheme as well as the AIMD approach are presented. Furthermore, the embedment of the floating orbital scheme into the molecular dynamics procedure is laid out. In the second part the floating orbital molecular dynamics scheme (FOMD) is validated in terms of energy and angular momentum conservation and finally applied to a single water molecule and a water dimer.

2 Methods 2.1

General scheme

The central issue of this work is the introduction and demonstration of the FOMD approach. Therefore, only the restricted Hartree–Fock approximation is used to determine the electronic structure. However, the extension of the FO approach to other and more sophisticated electronic structure methods has been demonstrated and thus, the extension of the FOMD approach will also be straightforward.27 The Born–Oppenheimer approximation is assumed for all simulations. In this section the floating orbital approach and its implementation in the MD scheme are explained. The term orbital floating covers the optimization of the basis function centers at a stationary molecular geometry. In the floating orbital molecular dynamics approach the orbital floating is performed for each MD step. The force on the atoms is then evaluated from the gradient with the basis functions displaced to their optimized position. Consequently, the trajectory with floated orbitals is not equal to the usual Hartree–Fock molecular dynamics (HFMD) simulations. An overview of the implementation is given in Fig. 1. 2.2

Optimization of the basis function center coordinates

In order to optimize the positions of the basis functions independent of the nuclear positions, another gradient apart from the conventional molecular gradient is required. The total energy E of the system is given by 1XX core E ¼ Vcc þ (1) Pnm Hmn þ Fmn 2 m n where Vcc denotes the nuclear–nuclear repulsion term, Pnm the elements of the density matrix, H core mn the elements of the oneelectron part of the Fock matrix and Fmn the elements of the Fock matrix with the latter ones being calculated according to eqn (2) and (3). H core mn = Tmn + Vmn

core þ Fmn ¼ Hmn

XX l

s

Pls

(2)

1 fm fn jfs fl fm fl jfs fn 2 (3)

This journal is © the Owner Societies 2014

View Article Online

Published on 14 February 2014. Downloaded by Arizona State University on 30/10/2014 17:34:20.

PCCP

Paper

Fig. 1 Implementation of the FOMD scheme. E, C and ea represent the electronic energy, the coeﬃcient matrix and the orbital energies, respectively. Dw and Df represent the convergence criteria for the positions of the basis centers wi and the forces fi. A detailed explanation is given in the subsequent sections.

Tmn stands for kinetic contributions, Vmn for the Coulombic terms and fi represents the basis functions. The ordinary gradient fA as the first derivative of E with respect to nuclear coordinates Ai with the basis functions attached to the atoms reads:38 core @ f f k f f X X @Hmn m n s l 1 fAi ¼ Pnm þ Pnm Pls 2 @A @A i i nm mnls (4) X @Smn @Vcc Qnm þ @Ai @Ai mn

in eqn (4) has been decomposed into its kinetic and Coulombic contributions leading to eqn (7)

fAi ¼

X

Pnm

nm

@Vmn X @Tmn @Vmn þ þ @Ai @wj @wj j/A

!!

X @ fm fn k fs fl 1X Pnm Pls þ 2 mnls @wi j/A

X mn

Qnm

X @Smn j/A

@wj

þ

(7)

@Vcc : @Ai

with S being the overlap matrix and:

1 fm fn k fs fl ¼ fm fn jfs fl fm fl jfs fn 2

Qnm ¼ 2

N=2 X

ea Cma Cna :

(5)

(6)

a

Here, ea denotes the SCF eigenvalues, and C is the matrix of SCF coeﬃcients in a system containing N electrons. The nuclear coordinates are consecutively indexed, i.e., A1 denotes the x-coordinate of nucleus one, A2 its y-coordinate and so on. The same holds for the coordinates of the basis function centers (w j) introduced below. In order to perform FOMD, it is necessary to formally distinguish between nuclear coordinates (Ai) and the coordinates of the basis function centers (w j). To illustrate this fact, the one-electron part of the Fock matrix

This journal is © the Owner Societies 2014

Care has to be taken for the index j running over all indices of basis functions assigned to the nucleus A ( j p A). The vector containing all nuclear coordinates Ai is denoted by A and w represents the respective vector of basis function positions w j. The vector fA with the entries fAi will be denoted in the following as molecular gradient. It can be read from eqn (7) that there are two distinct derivatives of the Coulomb term Vmn to be [email protected] sidered: with respect to nuclear coordinates and with @Ai respect to the centers of those basis functions which are @Vmn . If the positions of the attached to the respective atom @wj basis functions are optimized, another gradient needs to be evaluated, namely the first derivative of the total energy with respect to solely the coordinates w j of the basis functions fm.

Phys. Chem. Chem. Phys., 2014, 16, 6997--7005 | 6999

View Article Online

Paper

PCCP

The vector fw with the entries fwj will hereafter be termed electronic gradient and is calculated in the following way: ! X @Tmn @Vmn fwj ¼ Pnm þ @wj @wj nm

Published on 14 February 2014. Downloaded by Arizona State University on 30/10/2014 17:34:20.

@ fm fn k fs fl 1X þ Pnm Pls @wi 2 mnls

X mn

Qnm

(8)

@Smn : @wj

The electronic gradient holds another dimension compared to the molecular gradient because distinct subshells contribute to the entry of one atom into the molecular gradient while they have an own number in the electronic gradient vector. Furthermore, basis functions of higher angular momentum are floated together as one function independent of the magnetic quantum number, e.g., all three p functions are assumed to have the same position. It should be noted that in the case of optimized orbital positions, i.e., fw = 0, the molecular gradient reduces to the Hellman–Feynman force. The derivatives of the integrals are computed as integrals of derived Gaussian functions using Obara–Saika recurrence relations except for the first derivative of the Coulomb integral with respect to nuclear coordinates, which is evaluated employing Hermite expansion.39 2.3 Unification of floating orbitals with molecular dynamics simulations In this section, the implementation of the general FOMD approach is described. It is important to notice, that the orbital centers are not moving according to a fictitious dynamics as for example in Car–Parrinello molecular dynamics simulations.6 Instead, the positions of the basis functions play the role of additional variables in order to minimize the total energy. For this purpose, the BFGS algorithm was selected to minimize the total energy of the system with respect to the centers of the basis functions. Besides, there are different possible approaches for floating basis functions depending on which functions float individually. In this study, the most general approach is chosen, i.e., each set of contracted basis functions is treated individually. For higher angular momentum functions these are floated together independent of the magnetic quantum number. Hence, it can be assumed that the px, py and pz functions of one contraction have a common center. In the following the general scheme of a FOMD step is described. First, the initial system is subject to ordinary SCF calculations. In their initial positions w(0) the basis functions are located on their corresponding nuclei. To obtain new positions of the basis functions, the corresponding Newton equation (eqn (9)) is solved employing the inverse of an approximated Hessian (H1) and the electronic gradient (eqn (8)). w(n) w(n1) = (H1)(n1) fw(n1)

7000 | Phys. Chem. Chem. Phys., 2014, 16, 6997--7005

(9)

The superscripts in parentheses in eqn (9) denote the step number in the BFGS iteration with the initial quantities indexed by 0. In this work the unit matrix is employed as guess for the Hessian (H(0)). The displacement is scaled with a usual line search procedure using a constant maximal relative step length. At least one SCF calculation has to be carried out in the line search procedure. Afterwards, the Hessian is updated with the usual BFGS update, and a new step direction is computed using eqn (9). There are two criteria to judge the convergence of the BFGS algorithm. Therefore, two vectors are introduced (eqn (10) and (11)): ðnÞ ðn1Þ wj wj xcrit ½ j ¼ (10) max wj ; 1

f crit ½ j ¼

ðnÞ fwj max wj ; 1 maxðjE j; 1Þ

:

(11)

If the value of the uniform norm (JxcritJN and JfcritJN) for both vectors drops below the corresponding threshold, convergence is reached. If this is the case, the orbital positions are printed out, and the molecular gradient (eqn (4)) is evaluated to obtain the forces on the nuclei. The latter are then displaced according to the Velocity Verlet formalism.40 The orbitals are placed on their associate nuclei to be optimized again. The newly obtained nuclear positions are printed out, and the next FOMD step is started with a further SCF calculation. 2.4

Initial guess for positions of the basis function centers

A disadvantage of using the BFGS algorithm to optimize the basis function centers is the need for repeated SCF calculations and gradient evaluations cutting the profits of a smaller basis set. This is one of the reasons that the floating orbital approach is not superior to a larger basis set in normal geometry optimizations. As a solution to this problem, the application of a guess for the orbital positions has already been suggested by Huber.21 Thereby, the number of necessary steps to reach convergence in an optimization procedure may be reduced. For a molecular dynamics simulation, information about the previous step is available. The first idea is the utilization of this information to obtain a reasonable guess for the positions of the basis functions needed in the subsequent MD step. In this work a simple guess is used for the positions employing information about the centers from the previous FOMD step. A vector w(n) rel[k] is introduced describing the distance between all basis function positions and their corresponding nuclei in the kth MD step after the successful BFGS optimization. The elements of the vector are calculated according to eqn (12). Ai( j )[k] denotes the corresponding nuclear coordinate for w(n) [k] j w(n) [k] = w(n) [k] Ai( j )[k] rel j

(12)

(0) In the (k + 1) step w(n) rel[k] is used to obtain a guess for wj [k + 1] (eqn (13)).

w(0) [k + 1] = w(n) [k] Ai( j )[k] j rel

(13)

This journal is © the Owner Societies 2014

View Article Online

Published on 14 February 2014. Downloaded by Arizona State University on 30/10/2014 17:34:20.

PCCP

Paper

3 Computational details

4 Validation of the method

For the evaluation of molecular properties, diﬀerent molecular dynamics simulations have been carried out. Therefore, our own Fortran 2003 code FORPLEASURE has been used.41 The code contains a molecular dynamics procedure which can be applied with or without floating orbitals. So far no periodic boundary conditions and no thermostats are implemented, but to do this is straightforward and will be done in future versions. Reference calculations have been performed using the ORCA program package.42 All simulations were carried out employing Born–Oppenheimer molecular dynamics with Gaussian type orbitals. The basis sets applied are 3-21g (b1),43 4-31g (b2),44 and 6-311g* (b3).45 The temperature was not controlled by a thermostat. The simulations were conducted in the microcanonical (NVE) ensemble. The start atomic velocities are either initialized according to the Boltzmann distribution at a temperature of 400 K or set to zero, whereas the latter ones are labelled with a superscript zero. Simulations with zero initial velocity have been carried out for the following reason. In order to demonstrate the correctness of the implementation, results have been compared to the well-established ORCA program.42 In this code the initial velocities can be either set to zero or read from a file. To avoid numerical inaccuracies and provide comparability, simulations with zero initial velocity have been performed. The nuclei are propagated according to the Velocity Verlet algorithm40 for which the forces are evaluated from the molecular gradient. The integration timestep is set to 0.3 fs. Simulation times are given with the respective results. Simulations in which the guess has been applied are additionally designated with wG. Please note that the simulations are not equilibrated and are run for a short total simulation time. Therefore, all runs suﬀer from a non-equilibrated and shortly sampled trajectory. However, the present study serves to demonstrate the feasibility of the floating orbital method in conjunction with molecular dynamics simulations. Furthermore, the investigated systems are small and possess only few degrees of freedom. Therefore, these error sources are of minor importance here. In order to reduce the penalty due to insuﬃcient equilibration, the first 500 MD steps are not considered for the analyses in Section 5. Molecular integrals are evaluated applying the Obara–Saika recurrence relations, unless otherwise stated.46 SCF convergence is accelerated from the optimal damping algorithm.47 With regard to the reduction of computational cost, electron repulsion integrals are screened according to the Schwarz criterion with a threshold value of 106. The SCF convergence criterion is chosen to be 108. The maximal relative step length for the line search used in the BFGS method is set to 3 104. This value is sufficiently small for the investigated systems in the case studies. As explained in the Methods section, there are two criteria to judge convergence in the BFGS method. The threshold for JxcritJN is set to 108 for all performed calculations. In this work the threshold value for JfcritJN is set to 106. Power spectra are computed from TRAVIS.48 Two-dimensional plots are produced using XMGRACE.49

As introduced in Section 2.1, the force on the particles is evaluated with the basis functions being displaced from the respective nuclei. Therefore, the dynamics is aﬀected by the floating, and it needs to be verified that the drift of the conserved quantities is suﬃciently low. The drifts in the total energy as well as total angular momentum have been evaluated for the water monomer and the water dimer, respectively.

This journal is © the Owner Societies 2014

4.1

Conservation of the total energy

The drift of the total energy is evaluated as the slope of the linear regression and is given in Table 1. The drift of the total energy is lower than 1 kJ mol1 per picosecond in all cases. Especially for the dimer, there is a significant influence of the floating. Regarding the short simulation times of the floating orbital simulations, the value below 1 kJ (mol1 ps1) is considered to be suﬃciently low. 4.2

Conservation of total angular momentum

Also for the total angular momentum the drift DL has been calculated as the slope of the linear regression. The drifts are given in Table 2. If floating of the orbitals is performed, the drift of the total angular momentum is slightly larger as compared to the HFMD simulations. However, the observed drifts are suﬃciently small so that the total angular momentum is conserved.

5 Evaluation of molecular properties 5.1

Dipole moment

With the intention of demonstrating the influence of floating orbitals on molecular properties, the dipole moment has been investigated. The following simulations have been carried out employing initial velocities equal to zero. The drifts of the total energy and the total angular momentum of these simulations were found to be even smaller than those obtained with an Table 1

Drift in total energy DE in kJ (mol1 ps1), simulation time t in ps

Monomer

Dimer

Simulation

t

DE

FO/b1 FO/b1/wG HF/b1 HF/b2 HF/b3

6.0 6.0 15.0 15.0 15.0

1.5 1.3 2.9 4.8 1.9

Table 2

104 104 105 106 105

t

DE

5.4 6.0 15.0 15.0 11.4

1.9 1.0 9.4 3.4 4.0

103 101 104 104 104

Drift in total angular momentum DL in h/ps, simulation time t in ps

Monomer Simulation FO/b1 FO/b1/wG HF/b1 HF/b2 HF/b3

t 6.0 6.0 15.0 15.0 15.0

Dimer DL 8.9 4.8 5.8 2.7 1.1

9

10 108 1010 1010 109

t

DL

5.4 6.0 15.0 15.0 11.4

5.8 1.9 1.2 1.4 1.2

106 106 108 106 106

Phys. Chem. Chem. Phys., 2014, 16, 6997--7005 | 7001

View Article Online

Published on 14 February 2014. Downloaded by Arizona State University on 30/10/2014 17:34:20.

Paper

Fig. 2

PCCP

Total dipole moment of the water dimer (upper panel) and the water monomer (lower panel).

initial velocity according to the Boltzmann distribution at 400 K. In order to assess the correctness and the accuracy of the present implementation, results of HFMD simulations performed using both ForPleasure and the ORCA program have been compared.42 An excellent agreement has been observed for both the trajectories and the dipole moment with a deviation of less than 0.2% for the latter property. In Fig. 2 the dipole moment is plotted for 1.5 ps of the simulation time. Here, the HFMD reference calculations as obtained using ORCA are shown for the b1 basis (red line) and the b3 basis (blue line). The ForPleasure floating orbital simulation without the guess is represented by the black line. If the guess for the basis function positions is applied, this is indicated by additional symbols. Fig. 2 shows that the floating orbital simulations with and without the guess produce the same results. It needs to be pointed out that by applying the guess the computational cost could be reduced by half in the case of the shown simulations without loss of accuracy. During the simulation the molecule changes its geometry and thus, the values obtained along the trajectory should not be compared to the reference data for relaxed geometries. However, the molecular dynamics can be regarded as a fluctuation around the equilibrium geometry and may be considered for a rough comparison. There is a significant diﬀerence in both the magnitude and the fluctuation of the dipole moments. Values obtained by floating orbital simulations are in all cases lower than those from conventional HFMD simulations. For the water monomer the experimental value is 1.855 D50 while it is 2.643 D51 for the dimer. Experimental values are given in Table 3 together with the arithmetic mean value along the observed trajectories. It can be clearly seen that the Hartree–Fock method in conjunction with incomplete basis sets overestimates the

7002 | Phys. Chem. Chem. Phys., 2014, 16, 6997--7005

Table 3 Arithmetic mean of the dipole moment in D along the trajectories for the water monomer and the water dimer, simulation time t in ps. Initial velocity has been set to zero for all simulations

Simulation 0

FO /b1 FO0/b1/wG HF0/b1 HF0/b3 Experimental50,51

Monomer

Dimer

t

Dipole moment

t

Dipole moment

6.0 6.0 15.0 15.0

1.696 1.696 2.386 2.317 1.855

3.0 6.0 6.0 3.0

2.688 2.689 3.538 3.347 2.643

dipole moment. The results obtained by floating orbital simulations agree much better with experimental values, especially in the case of the dimer, although only the 3-21g basis set has been applied. In the case of the water monomer the dipole moment is underestimated by the floating orbital approach, but still closer to the experimental value. 5.2

Power spectra

It has already been mentioned that the dynamics is aﬀected by the floating orbital approach. In order to investigate the influence, power spectra have been obtained from the trajectories using the TRAVIS program.48 The purpose of this section is not the quantitative analysis of the signals and their intensity, but rather an investigation of trends in the peak positions depending on the basis set and the eﬀect of floating. The power spectrum of the water monomer is shown in Fig. 3 for the relevant wavenumbers. Again, the FOMD simulations with and without guess provide identical results indicating suﬃcient convergence during the optimizations. Regarding the peaks below 2000 cm1 a good agreement with the HFMD simulation applying the b3 basis set can be observed. The peak

This journal is © the Owner Societies 2014

View Article Online

Published on 14 February 2014. Downloaded by Arizona State University on 30/10/2014 17:34:20.

PCCP

Fig. 3

Paper

Power spectra of the water monomer as obtained from the diﬀerent trajectories.

positions of the smaller basis sets occur at smaller wavenumbers. A similar trend can be seen for the signals above 4000 cm1. Both b1 and b2 basis sets show smaller wavenumbers in the HFMD simulations, whereas the position of the FOMD/b1 simulations is shifted towards the HFMD/b3 result. However, the agreement of both curves is slightly worse in this case. Overall, the positions of the signals are shifted towards the larger basis set, if the floating is applied. The relevant regions of the power spectrum of the water dimer are shown in Fig. 4. In this case a diﬀerence between the

FOMD simulations with and without guess is observed in the intensity of the peak at around 4100 cm1. This is, however, only due to the diﬀerent simulation times. Again a good agreement for the peak positions of the HFMD/b3 approach and the FOMD simulation in the region below 2000 cm1 is observed. On the other hand, also the HFMD/b1 simulation reproduces the peak position of the larger basis quite well. Only the HFMD/b2 simulation is significantly oﬀ so that the influence of the floating cannot be judged here. The positions of the peaks above 4000 cm1 show again a positive eﬀect of the floating.

Fig. 4 Power spectra of the water dimer as obtained from the diﬀerent trajectories.

This journal is © the Owner Societies 2014

Phys. Chem. Chem. Phys., 2014, 16, 6997--7005 | 7003

View Article Online

Published on 14 February 2014. Downloaded by Arizona State University on 30/10/2014 17:34:20.

Paper

HFMD simulations at smaller basis sets show the peaks at wavenumbers that are too low. Compared to that, the HFMD/b3 peak position is shifted to larger wavenumbers which can be reproduced by the FOMD simulations. The peak shape, in contrast, does not agree that well because the FOMD simulations produce a smaller peak width. Overall, the floating orbital approach is capable of influencing the dynamics in a way that the power spectra show shifts of the peaks towards the positions as obtained by larger basis sets.

6 Conclusions The well-known approach of floating orbitals, which up to now has only been applied to geometry optimizations and the calculation of molecular properties, has been included in a molecular dynamics scheme. The basic theory and necessary algorithms have been introduced, and a guess for the basis function centers has been suggested. The FOMD approach was found to conserve both total energy and angular momentum of the test systems water monomer and dimer to obtain suﬃcient accuracy given the low equilibration time. Concerning code performance no runtimes are given, because up to now the code has not been suﬃciently optimized. Due to repeated SCF and gradient calculations the floating increases the computational cost by two orders of magnitude. However, applying the initial guess to the basis function centers is capable of reducing the runtime by 50%. In order to judge the influence of the proposed method, the dipole moment has been calculated along the trajectory. For this property the eﬀect of the floating orbitals is advantageous, which is due to the fact that this property benefits from the possibility of the centers of charges to move. The eﬀect of polarization and the quality of properties related to it are covered by polarization functions in state-of-the-art calculations. The floating is capable of replacing these computationally demanding atomic orbitals. However, due to the need for repeated energy and gradient evaluations during the optimization of the basis function centers the FOMD approach is computationally very demanding. This problem can be treated in future versions by applying diﬀerent algorithms and implementation of parallelization. A similar trend is observed for the power spectra as representatives for the dynamics of the diﬀerent trajectories. If the 3-21g basis set is applied the observed peak positions are shifted by the FOMD approach towards HFMD positions as obtained with the 6-311g* basis set. The calculation of accurate vibrational spectra and dipole moments is still a challenging task in simulations. This work demonstrates that FOMD will be a suitable approach to increase the accuracy of these properties. In forthcoming studies the FOMD approach will be applied in order to investigate further properties on the fly.

Acknowledgements This work has been financially supported by DFG projects KI-768/5-3 and KI-768/8-1 as well as the SFB 813 ‘‘spin centers’’.

7004 | Phys. Chem. Chem. Phys., 2014, 16, 6997--7005

PCCP

This work is funded by the European Union and the Free State of Saxony. Computer time from the RZ Leipzig is gratefully acknowledged. The authors thank B. Intemann for careful reading of the manuscript.

References 1 K. Laasonen, M. Sprik, M. Parrinello and R. Car, J. Chem. Phys., 1993, 99, 9080–9089. 2 H. Huber, A. J. Dyson and B. Kirchner, Chem. Soc. Rev., 1999, 28, 121–133. 3 I. Bak, J. Hutter and G. Plinks, J. Phys. Chem. A, 2006, 110, 2188–2194. 4 M. Tuckerman, K. Laasonen, M. Sprik and M. Parrinello, J. Phys. Chem., 1995, 99, 5749–5752. 5 B. L. Bhargava, Y. Yasaka and M. L. Klein, Chem. Commun., 2011, 47, 6228–6241. 6 B. Kirchner, P. J. di Dio and J. Hutter, Top. Curr. Chem., 2012, 307, 109–154. 7 R. Car and M. Parrinello, Phys. Rev. Lett., 1985, 55, 2471–2474. ¨hne, M. Krack, F. R. Mohamed and M. Parrinello, 8 T. D. Ku Phys. Rev. Lett., 2007, 98, 066401. ¨hne, M. Krack and M. Parrinello, J. Chem. Theory 9 T. D. Ku Comput., 2009, 5, 235–241. 10 A. M. N. Niklasson, Phys. Rev. Lett., 2008, 100, 123004. 11 J. Jakowski and K. Morokuma, J. Chem. Phys., 2009, 130, 224106. 12 G. Zheng, A. M. N. Niklasson and M. Karplus, J. Chem. Phys., 2011, 135, 044122. 13 T. Helgaker, P. Jørgensen and J. Olsen, Molecular ElectronicStructure Theory, John Wiley & Sons, Chichester, West-Sussex, United Kingdom, 1st edn, 2000, ch. 6. 14 A. C. Hurley, Proc. R. Soc. London, Ser. A, 1954, 226, 170–178. 15 H. Shull and D. D. Ebbing, J. Chem. Phys., 1958, 28, 866–870. 16 A. A. Frost, J. Chem. Phys., 1967, 47, 3707–3713. 17 A. A. Frost, J. Chem. Phys., 1967, 47, 3714–3716. 18 A. A. Frost, J. Phys. Chem., 1968, 72, 1289–1293. ´sza ´r, Chem. Phys. Lett., 2007, 438, 19 G. Tasi and A. G. Csa 139–143. 20 B. Ford, G. G. Hall and J. C. Packer, Int. J. Quantum Chem., 1970, 4, 533–539. 21 H. Huber, Chem. Phys. Lett., 1979, 62, 95–99. 22 H. Huber, Chem. Phys. Lett., 1980, 70, 353–357. 23 H. Huber, Theor. Chim. Acta, 1980, 55, 117–126. 24 H. Huber, THEOCHEM, 1981, 1, 277–284. 25 A. C. Hurley, J. Comput. Chem., 1988, 9, 75–79. ¨f, J. Chem. Phys., 1988, 89, 26 T. Helgaker and J. Almlo 4889–4902. 27 A. H. Pakiari and A. Mohajeri, Int. J. Mod. Phys. C, 2002, 13, 1095–1103. 28 S. Wilson and D. Moncrieﬀ, Mol. Phys., 1993, 80, 461–467. 29 D. Moncrieﬀ and S. Wilson, Mol. Phys., 1994, 82, 523–530. 30 V. N. Glushkov and S. Wilson, Int. J. Quantum Chem., 2007, 107, 2632–2642.

This journal is © the Owner Societies 2014

View Article Online

Published on 14 February 2014. Downloaded by Arizona State University on 30/10/2014 17:34:20.

PCCP

31 J. L. Whitten, J. Chem. Phys., 1966, 44, 359–364. 32 O. Marsalek, F. Uhlig, J. VandeVondele and P. Jungwirth, Acc. Chem. Res., 2012, 45, 23–32. 33 B. A. Murtagh and R. W. H. Sargent, Comput. J., 1970, 13, 185–194. 34 C. G. Broyden, Math. Comput., 1967, 21, 368–381. 35 R. Fletcher, Comput. J., 1970, 13, 317–322. 36 D. Goldfarb, Math. Comput., 1970, 24, 23–26. 37 D. F. Shanno, Math. Comput., 1970, 24, 647–656. 38 A. Szabo and N. S. Ostlund, Modern Quantum Chemistry, Dover Publications Inc., New York, USA, 1st edn, 1996, pp. 440–442. 39 L. E. McMurchie and E. R. Davidson, J. Comput. Phys., 1978, 26, 218–231. 40 L. Verlet, Phys. Rev., 1967, 159, 98–103. ¨ssel 41 ForPleasure V1.1 Copyright B. Kirchner, written by M. Bru and E. Perlt 2011-2012, University of Leipzig, WilhelmOstwald-Institute of Physical and Theoretical Chemistry, 2012.

This journal is © the Owner Societies 2014

Paper

42 F. Neese, Wiley Interdiscip. Rev. Comput. Mol. Sci., 2012, 2, 73–78. 43 J. S. Binkley, J. A. Pople and W. J. Hehre, J. Am. Chem. Soc., 1980, 102, 939–947. 44 R. Ditchfield, W. J. Hehre and J. A. Pople, J. Chem. Phys., 1971, 54, 724–728. 45 R. Krishnan, J. S. Binkley, R. Seeger and J. A. Pople, J. Chem. Phys., 1980, 72, 650–654. 46 S. Obara and A. Saika, J. Chem. Phys., 1986, 84, 3963–3974. `s and C. Le Bris, Int. J. Quantum Chem., 2000, 79, 47 E. Cance 82–90. 48 M. Brehm and B. Kirchner, J. Chem. Inf. Model., 2011, 51, 2007–2023. 49 Grace, (c) 1996-2008 Grace Development Team, see http:// plasma-gate.weizmann.ac.il/Grace, (accessed January 08, 2013). 50 F. J. Lovas, J. Phys. Chem. Ref. Data, 1978, 7, 1445–1750. 51 T. R. Dyke, K. M. Mack and J. S. Muenter, J. Chem. Phys., 1977, 66, 498–510.

Phys. Chem. Chem. Phys., 2014, 16, 6997--7005 | 7005