Density Functional and Molecular Dynamics Method

The electronic structures were computed using density functional theory (DFT) (see e.g. ref DFT). For the exchange-correlation functional we used the Perdew-Zunger [43] parameterization of the local density approximation (LDA) for this functional, which is based on the free-electron MC simulations by Ceperley and Alder [35]. Density gradient corrections were added, namely the Becke-88 gradient correction for exchange[67] and the Perdew-86 gradient correction for correlation.[38]

The Car-Parrinello (CP) method[49] was applied to perform both dynamic and static calculations. The method performs the classical MD and simultaneously applies DFT to describe the electronic structure, using an extended Lagrangian formulation. The characteristic feature of the Car-Parrinello approach is that the electronic wave function, i.e. the coefficients of the plane wave basis set, are dynamically optimized to be consistent with the changing positions of the atomic nuclei. The actual implementation involves the numerical integration of the equations of motion of second-order Newtonian dynamics. A crucial parameter in this scheme is the fictitious mass $m_{\rm e}$ associated with the dynamics of the electronic degrees of freedom. In practice, $m_{\rm e}$ has to be chosen small enough to ensure fast wave function adaptation to the changing nuclear positions on one hand and sufficiently large to have a workable large time step on the other. In the present work, the equations of motion were integrated using the Verlet [68] algorithm with a mass $m_{\rm e}=1000 ~\mathrm{a.u.}=9.10939 10^{-28}$ kg, which limits the time step to $\delta t=0.19$ fs. To maintain a constant temperature of $T=300$ K a Nosé thermostat [69] was applied with a period of 100 fs.

The CP simulations have been performed using the CP-PAW code package developed by Blöchl. It implements the ab-initio molecular dynamics together with the projector augmented wave (PAW) method.[24] The PAW method uses an augmented plane wave basis for the electronic valence wave functions, and, in the current implementation, frozen atomic wave functions for the core states. Thus it is able to produce the correct wave function and densities also close to the nucleus, including the correct nodal structure of the wave functions. The advantages compared to the pseudopotential approach are that transferability problems are largely avoided, that quantities such as hyperfine parameters and electric field gradients are obtained with high accuracy [33,34] and, most important for the present study, that a smaller basis set as compared to traditional norm-conserving pseudopotentials is required. The frozen core approximation was applied for the 1s electrons of C and O, and up to 2p for Cl. For H, C and O, one projector function per angular-momentum quantum number was used for $s$- and $p$-angular momenta. For Cl, two projector functions were used for $s$- and one for $p$-angular momenta. The Kohn-Sham orbitals of the valence electrons were expanded in plane waves up to a kinetic energy cutoff of 30 Ry.

For reference, we also performed static DFT calculations using the atomic-orbital based ADF package.[70] In these calculations, the Kohn-Sham orbitals were expanded in an uncontracted triple-$\zeta$ Slater-type basis set augmented with one 2p and one 3d polarization function for H, 3d and 4f polarization functions for C, O, and Cl.