Multiconfiguration Time-dependent Hartree Research Articles

Prediction through quantum dynamics simulations: Photo-excited cyclobutanone.

Quantum dynamics simulations are becoming a standard tool for simulating photo-excited molecular systems involving a manifold of coupled states, known as non-adiabatic dynamics. While these simulations have had many successes in explaining experiments and giving details of non-adiabatic transitions, the question remains as to their predictive power. In this work, we present a set of quantum dynamics simulations on cyclobutanone using both grid-based multi-configuration time-dependent Hartree and direct dynamics variational multi-configuration Gaussian methods. The former used a parameterized vibronic coupling model Hamiltonian, and the latter generated the potential energy surfaces on the fly. The results give a picture of the non-adiabatic behavior of this molecule and were used to calculate the signal from a gas-phase ultrafast electron diffraction (GUED) experiment. Corresponding experimental results will be obtained and presented at a later stage for comparison to test the predictive power of the methods. The results show that over the first 500fs after photo-excitation to the S2 state, cyclobutanone relaxes quickly to the S1 state, but only a small population relaxes further to the S0 state. No significant transfer of population to the triplet manifold is found. It is predicted that the GUED experiments over this time scale will see signals related mostly to the C-O stretch motion and elongation of the molecular ring along the C-C-O axis.

The photochemistry of Rydberg-excited cyclobutanone: Photoinduced processes and ground state dynamics.

Owing to ring strain, cyclic ketones exhibit complex excited state dynamics with multiple competing photochemical channels active on the ultrafast timescale. While the excited state dynamics of cyclobutanone after π* ← n excitation into the lowest-energy excited singlet (S1) state has been extensively studied, the dynamics following 3s ← n excitation into the higher-lying singlet Rydberg (S2) state are less well understood. Herein, we employ fully quantum multiconfigurational time-dependent Hartree (MCTDH) simulations using a model Hamiltonian as well as "on-the-fly" trajectory-based surface-hopping dynamics (TSHD) simulations to study the relaxation dynamics of cyclobutanone following 3s ← n excitation and to predict the ultrafast electron diffraction scattering signature of these relaxation dynamics. Our MCTDH and TSHD simulations indicate that relaxation from the initially-populated singlet Rydberg (S2) state occurs on the timescale of a few hundreds of femtoseconds to a picosecond, consistent with the symmetry-forbidden nature of the state-to-state transition involved. There is no obvious involvement of excited triplet states within the timeframe of our simulations (<2ps). After non-radiative relaxation to the electronic ground state (S0), vibrationally hot cyclobutanone has sufficient internal energy to form multiple fragmented products including C2H4 + CH2CO (C2; 20%) and C3H6 + CO (C3; 2.5%). We discuss the limitations of our MCTDH and TSHD simulations, how these may influence the excited state dynamics we observe, and-ultimately-the predictive power of the simulated experimental observable.

Ultrafast Excited-State Nonadiabatic Dynamics in Pt(II) Donor-Bridge-Acceptor Assemblies: A Quantum Approach for Optical Control.

The ultrafast nonadiabatic excited state dynamics of (PTZ-N-benzyl-acetylide) (trans-bis-trimethylphosphine) Pt(II) (acetylide-NDI-bis-methyl) 1, representative of a series of Pt(II) donor-bridge-acceptor assemblies experimentally studied by the Weinstein group, University of Sheffield, is investigated by means of wavepacket propagations based on the multiconfiguration time-dependent Hartree (MCTDH) method. On the basis of electronic structure data obtained at the time-dependent density functional theory (TD-DFT) level, the subpicosecond decay is simulated by solving an 11 electronic states multimode problem, up to 18 vibrational normal modes, including both spin-orbit coupling (SOC) and vibronic coupling. A careful analysis of the results, within the diabatic representation, provides the key features of the spin-vibronic mechanism at work in this complex, distinguishing between the spin-orbit and vibronically activated ultrafast processes within the excited states manifold. The knowledge of the key active normal modes that promote selectively the population of specific electronic excited states opens a route toward optical control by selectively exciting these modes in order to drive the associated nonadiabatic processes. Relevant simulations, over 2 ps, are proposed to assess the impact of these selective vibrational excitations on the branching ratio between the primary photoproducts, namely, bridge-acceptor charge-transfer (CT) and donor-acceptor charge-separated (CS) electronic states. Whereas the excitation of the localized acetylide bridge C≡C bond stretching does not modify drastically the population of the low-lying electronic states within the first two ps, vibrational excitation of the out-of-plane twisting motion of the N-benzyl group linked to the donor entity favors the population of the 1,3CS states at the expense of the lowest 1,3CT states. This quantum study opens the route to IR optical control experiments based on the specific alteration of vibrational normal modes that activate vibronic couplings between key electronic excited states. However, the presence of critical crossings along the PES channels associated with these normal modes and the role of concurrent SOC driven ultrafast transfers of population should not be underestimated.

QuTree: A tree tensor network package.

We present QuTree, a C++ library for tree tensor network approaches. QuTree provides class structures for tensors, tensor trees, and related linear algebra functions that facilitate the fast development of tree tensor network approaches such as the multilayer multiconfigurational time-dependent Hartree approach or the density matrix renormalization group approach and its various extensions. We investigate the efficiency of relevant tensor and tensor network operations and show that the overhead for managing the network structure is negligible, even in cases with a million leaves and small tensors. QuTree focuses on providing simple, high-level routines while retaining easy access to the backend to facilitate novel developments. We demonstrate the capabilities of the package by computing the eigenstates of coupled harmonic oscillator Hamiltonians and performing random circuit simulations on a virtual quantum computer.

Excitation energy transfer and vibronic relaxation through light-harvesting dendrimer building blocks: A nonadiabatic perspective.

The light-harvesting excitonic properties of poly(phenylene ethynylene) (PPE) extended dendrimers (tree-like π-conjugated macromolecules) involve a directional cascade of local excitation energy transfer (EET) processes occurring from the "leaves" (shortest branches) to the "trunk" (longest branch), which can be viewed from a vibronic perspective as a sequence of internal conversions occurring among a connected graph of nonadiabatically coupled locally excited electronic states via conical intersections. The smallest PPE building block that is able to exhibit EET, the asymmetrically meta-substituted PPE oligomer with one acetylenic bond on one side and two parallel ones on the other side (hence, 2-ring and 3-ring para-substituted pseudo-fragments), is a prototype and the focus of the present work. From linear-response time-dependent density functional theory electronic-structure calculations of the molecule as regards its first two nonadiabatically coupled, optically active, singlet excited states, we built a (1 + 2)-state-8-dimensional vibronic-coupling Hamiltonian model for running subsequent multiconfiguration time-dependent Hartree wavepacket relaxations and propagations, yielding both steady-state absorption and emission spectra as well as real-time dynamics. The EET process from the shortest branch to the longest one occurs quite efficiently (about 80% quantum yield) within the first 25fs after light excitation and is mediated vibrationally through acetylenic and quinoidal bond-stretching modes together with a particular role given to the central-ring anti-quinoidal rock-bending mode. Electronic and vibrational energy relaxations, together with redistributions of quantum populations and coherences, are interpreted herein through the lens of a nonadiabatic perspective, showing some interesting segregation among the foremost photoactive degrees of freedom as regards spectroscopy and reactivity.

A tensor network view of multilayer multiconfiguration time-dependent Hartree methods

The multilayer multiconfiguration time-dependent Hartree (ML-MCTDH) method and the density matrix renormalisation group (DMRG) are powerful workhorses applied mostly in different scientific fields. Although both methods are based on tensor network states, very different mathematical languages are used for describing them. This severely limits knowledge transfer and sometimes leads to re-inventions of ideas well known in the other field. Here, we review ML-MCTDH and DMRG theory using both MCTDH expressions and tensor network diagrams. We derive the ML-MCTDH equations of motions using diagrams and compare them with time-dependent and time-independent DMRG algorithms. We further review two selected recent advancements. The first advancement is related to optimising unoccupied single-particle functions in MCTDH, which corresponds to subspace enrichment in the DMRG. The second advancement is related to finding optimal tree structures and on highlighting similarities and differences of tensor networks used in MCTDH and DMRG theories. We hope that this contribution will foster more fruitful cross-fertilisation of ideas between ML-MCTDH and DMRG.

Coherent Transient Localization Mechanism of Interchain Exciton Transport in Regioregular P3HT: A Quantum-Dynamical Study.

Transient localization has been proposed as a transport mechanism in organic materials, for both charge carriers and excitons. Here, we characterize a quantum coherent transient localization mechanism using full quantum simulations of an H-aggregated model system representative of regioregular poly(3-hexylthiophene) (rrP3HT). A Frenkel-Holstein Hamiltonian parametrized from first principles is considered, including local high-frequency modes and anharmonic, site-correlated interchain modes. Quantum-dynamical calculations are carried out using the Multi-Layer Multi-Configuration Time-Dependent Hartree (ML-MCTDH) method for a 13-site system with 195 vibrational modes, under periodic boundary conditions. It is shown that temporary localization of exciton polarons alternates with resonant transfer driven by interchain modes. While the transport process is mainly determined by exciton-polarons at the low-energy band edge, persistent coupling with the excitonic manifold is observed, giving rise to a nonadiabatic excitonic flux. This elementary transport mechanism remains preserved for limited static disorder and gives way to Anderson localization when the static disorder becomes dominant.

Quantum dynamics studies on reactive scattering of Cl with rotationally excited H[formula omitted]

In this work, extensive multi-configurational time-dependent Hartree (MCTDH) calculations in conjugation with the neural-network (NN) fitted potential energy surface (Science, 347 (2015), 60) were launched to study dynamics of Cl reacting with rotationally excited H2. Due to its direct abstraction mechanism, Cl + H2→ HCl + H in the ground state has no clear reaction resonances, which is well reproduced Here. For the reaction of Cl with H2 in rotationally excited states, the reactive probabilities are enhanced due to a Polanyi-like rule for rotationally excited states. Such Polanyi-like rule can be interpreted through the similar way as for the cases of vibrationally excited states. Moreover, distinctive resonances peaks are investigated for these cases implying a model for generating a quasi-bound rotation state. Similar to the chemical bond softening model for the vibrationally excited cases, a probable model is proposed in this work.

Eigenstate calculation in the state-averaged (multi-layer) multi-configurational time-dependent Hartree approach.

A new approach for the calculation of eigenstates with the state-averaged (multi-layer) multi-configurational time-dependent Hartree (MCTDH) approach is presented. The approach is inspired by the recent work of Larsson [J. Chem. Phys. 151, 204102 (2019)]. It employs local optimization of the basis sets at each node of the multi-layer MCTDH tree and successive downward and upward sweeps to obtain a globally converged result. At the top node, the Hamiltonian represented in the basis of the single-particle functions (SPFs) of the first layer is diagonalized. Here p wavefunctions corresponding to the p lowest eigenvalues are computed by a block Lanczos approach. At all other nodes, a non-linear operator consisting of the respective mean-field Hamiltonian matrix and a projector onto the space spanned by the respective SPFs is considered. Here, the eigenstate corresponding to the lowest eigenvalue is computed using a short iterative Lanczos scheme. Two different examples are studied to illustrate the new approach: the calculation of the vibrational states of methyl and acetonitrile. The calculations for methyl employ the single-layer MCTDH approach, a general potential energy surface, and the correlation discrete variable representation. A five-layer MCTDH representation and a sum of product-type Hamiltonian are used in the acetonitrile calculations. Very fast convergence and order of magnitude reductions in the numerical effort compared to the previously used block relaxation scheme are found. Furthermore, a detailed comparison with the results of Avila and Carrington [J. Chem. Phys. 134, 054126 (2011)] for acetonitrile highlights the potential problems of convergence tests for high-dimensional systems.

Quasi-direct Quantum Molecular Dynamics: The Time-Dependent Adaptive Density-Guided Approach for Potential Energy Surface Construction.

We present a new quasi-direct quantum molecular dynamics computational method which offers a compromise between quantum dynamics using a precomputed potential energy surface (PES) and fully direct quantum dynamics. This method is termed the time-dependent adaptive density-guided approach (TD-ADGA) and is a method for constructing a PES on the fly during a dynamics simulation. This is achieved by acquisition of new single-point (SP) calculations and refitting of the PES, depending on the need of the dynamics. The TD-ADGA is a further development of the adaptive density-guided approach (ADGA) for PES construction where the placement of SPs is guided by the density of the nuclear wave function. In TD-ADGA, the ADGA framework has been integrated into the time propagation of the time-dependent nuclear wave function and we use the reduced one-mode density of this wave function to guide when and where new SPs are placed. The PES is thus extended or updated if the wave function moves into new areas or if a certain area becomes more important. Here, we derive equations for the reduced one-mode density for the time-dependent Hartree (TDH) method and for multiconfiguration time-dependent Hartree (MCTDH) methods, but the TD-ADGA can be used with any time-dependent wave function method as long as a density is available. The TD-ADGA method has been investigated on molecular systems containing single- and double-minimum potentials and on single-mode and multi-mode systems. We explore different approaches to handle the fact that the TD-ADGA involves a PES that changes during the computation and show how results can be obtained that are in very good agreement with results obtained by using an accurate reference PES. Dynamics with TD-ADGA is essentially a black box procedure, where only the initialization of the system and how to compute SPs must be provided. The TD-ADGA thus makes it easier to carry out quantum molecular dynamics and the quasi-direct framework opens up the possibility to compute quantum dynamics accurately for larger molecular systems.

Quantum and classical molecular dynamics for H atom scattering from graphene.

This work presents systematic comparisons between classical molecular dynamics (cMD) and quantum dynamics (QD) simulations of 15-dimensional and 75-dimensional models in their description of H atom scattering from graphene. We use an experimentally validated full-dimensional neural network potential energy surface of a hydrogen atom interacting with a large cell of graphene containing 24 carbon atoms. For quantum dynamics simulations, we apply Monte Carlo canonical polyadic decomposition to transform the original potential energy surface (PES) into a sum of products form and use the multi-layer multi-configuration time-dependent Hartree method to simulate the quantum scattering of a hydrogen or deuterium atom with an initial kinetic energy of 1.96 or 0.96 eV and an incident angle of 0°, i.e., perpendicular to the graphene surface. The cMD and QD initial conditions have been carefully chosen in order to be as close as possible. Our results show little differences between cMD and QD simulations when the incident energy of the H atom is equal to 1.96 eV. However, a large difference in sticking probability is observed when the incident energy of the H atom is equal to 0.96 eV, indicating the predominance of quantum effects. To the best of our knowledge, our work provides the first benchmark of quantum against classical simulations for a system of this size with a realistic PES. Additionally, new projectors are implemented in the Heidelberg multi-configuration time-dependent Hartree package for the calculation of the atom scattering energy transfer distribution as a function of outgoing angles.

Accurate Potential Energy Surfaces Using Atom-Centered Potentials and Minimal High-Level Data.

We demonstrate that a Δ-density functional theory (Δ-DFT) approach based on atom-centered potentials (ACPs) represents a computationally inexpensive and accurate method for representing potential energy surfaces (PESs) for the HONO and HFCO molecules and vibrational frequencies derived therefrom. Using as few as 100 CCSD(T)-F12a reference energies, ACPs developed for use with B3LYP/def2-TZVPP are shown to produce PESs for HONO and HFCO with mean absolute errors of 27.7 and 5.8 cm-1, respectively. Application of the multiconfigurational time-dependent Hartree (MCTDH) method with ACP-corrected B3LYP/def2-TZVPP PESs produces vibrational frequencies for cis- and trans-HONO with mean absolute percent errors (MAPEs) of 0.8 and 1.1, compared to 0.8 obtained for the two isomers with CCSD(T)-F12a/cc-pVTZ-F12/MCTDH. For HFCO, the vibrational frequencies obtained using the present (Δ-DFT)/MCTDH approach give a MAPE of 0.1, which is the error obtained with CCSD(T)-F12a/cc-pVTZ-F12/MCTDH. The ACP approach is therefore successful in representing a PES calculated at a high level of theory (CCSD(T)-F12a) and a promising method for the development of a general protocol for the representation of accurate molecular PESs and the calculation of molecular properties from them.

Temperature Dependence of the Electronic Absorption Spectrum of NO2.

The nitrogen dioxide (NO2) radical is composed of the two most abundant elements in the atmosphere, where it can be formed in a variety of ways including combustion, detonation of energetic materials, and lightning. Relevant also to smog and ozone cycles, together these processes span a wide range of temperatures. Remarkably, high-resolution NO2 electronic absorption spectra have only been reported in a narrow range below about 300 K. Previously, we reported [ J. Phys. Chem. A 2021, 125, 5519-5533] the construction of quasi-diabatic potential energy surfaces (PESs) for the lowest four electronic states (X̃, Ã, B̃, and C̃) of NO2. In addition to three-dimensional PESs based on explicitly correlated MRCI(Q)-F12/VTZ-F12 ab initio data, the geometry dependence of each component of the dipoles and transition dipoles was also mapped into fitted surfaces. The multiconfigurational time-dependent Hartree (MCTDH) method was then used to compute the 0 K electronic absorption spectrum (from the ground rovibrational initial state) employing those energy and transition dipole surfaces. Here, in an extension of that work, we report an investigation into the effects of elevated temperature on the spectrum, considering the effects of the population of rotationally and vibrationally excited initial states. The calculations are complemented by new experimental measurements. Spectral contributions from hundreds of rotational states up to N = 20 and from 200 individually-characterized vibrational states were computed. A spectral simulation tool was developed that enables modeling the spectrum at various temperatures─by weighting individual spectral contributions via the partition function, or for pure excited initial states, which can be probed via transient absorption spectroscopy. We validate these results against experimental absorption spectroscopy data at high temperatures, as well as via a new measurement from the (1,0,1) initial vibrational state.

A numerically exact correlation discrete variable representation for multi-configurational time-dependent Hartree calculations.

The correlation discrete variable representation (CDVR) enables (multilayer) multi-configurational time-dependent Hartree (MCTDH) calculations with general potentials. The CDVR employs a set of grids corresponding to single-particle functions to efficiently evaluate all potential matrix elements appearing in the MCTDH equations of motion. In standard CDVR approaches, the number of grid points employed is tied to the number of corresponding single-particle functions. This limits the accuracy of the quadrature, which can be achieved for a given single-particle function basis. In this work, an extended CDVR approach that facilitates a numerically exact quadrature of all potential matrix elements is introduced. The number of grid points employed can be increased independent of the number of corresponding single-particle function to achieve any desired quadrature accuracy. The properties of the new scheme are illustrated by numerical calculations studying the photodissociation of NOCl and the vibrational states of CH3. Fast convergence with respect to the number of additional quadrature points is observed: Employing a single additional point in each physical or logical coordinate already ensures negligible quadrature errors.

Exploring disordered quantum spin models with a multilayer multiconfigurational approach

Numerical simulations of quantum spin models are crucial for a profound understanding of many-body phenomena in a variety of research areas in physics. An outstanding problem is the availability of methods to tackle systems that violate area-laws of entanglement entropy. Such scenarios cover a wide range of compelling physical situations including disordered quantum spin systems among others. In this work, we employ a numerical technique referred to as multi-layer multi-configuration time-dependent Hartree (ML-MCTDH) to evaluate the ground state of several disordered spin models. ML-MCTDH has previously been used to study problems of high-dimensional quantum dynamics in molecular and ultracold physics but is here applied to study spin systems for the first time. We exploit the inherent flexibility of the method to present results in one and two spatial dimensions and treat challenging setups that incorporate long-range interactions as well as disorder. Our results suggest that the hierarchical multi-layering inherent to ML-MCTDH allows to tackle a wide range of quantum many-body problems such as spin dynamics of varying dimensionality.

Efficient simulation of multielectron dynamics in molecules under intense laser pulses: implementation of the multiconfiguration time-dependent Hartree–Fock method based on the adaptive finite element method

We present an implementation of the multiconfiguration time-dependent Hartree–Fock method based on the adaptive finite element method for molecules under intense laser pulses. For efficient simulations, orbital functions are propagated by a stable propagator using the short iterative Arnoldi scheme and our implementation is parallelized for distributed memory computing. This is demonstrated by simulating high-harmonic generation from a water molecule and achieves a simulation of multielectron dynamics with overwhelmingly less computational time, compared to our previous work [Sawada, R. et al. Phys. Rev. A, 2016, 93, 023434].

Confining CO2 inside sI clathrate‐hydrates: The impact of the CO2–water interaction on quantized dynamics

We report new results on the translational-rotational (T-R) states of the CO2 molecule inside the sI clathrate-hydrate cages. We adopted the multiconfiguration time-dependent Hartree methodology to solve the nuclear molecular Hamiltonian, and to address issues on the T-R couplings. Motivated by experimental X-ray observations on the CO2 orientation in the D and T sI cages, we aim to evaluate the effect of the CO2 -water interaction on quantum dynamics. Thus, we first compared semiempirical and ab initio-based pair interaction model potentials against first-principles DFT-D calculations for ascertaining the importance of nonadditive many-body effects on such guest-host interactions. Our results reveal that the rotational and translational excited states quantum dynamics is remarkably different, with the pattern and density of states clearly affected by the underlying potential model. By analyzing the corresponding the probability density distributions of the calculated T-R eigenstates on both semiempirical and ab initio pair CO2 -water nanocage potentials, we have extracted information on the altered CO2 guest local structure, and we discussed it in connection with experimental data on the orientation of the CO2 molecule in the D and T sI clathrate cages available from neutron diffraction and solid-state NMR studies, as well as in comparison with previous molecular dynamics simulations. Our calculations provide a very sensitive test of the potential quality by predicting the low-lying T-R states and corresponding transitions for the encapsulated CO2 molecule. As such spectroscopic observables have not been measured so far, our results could trigger further detailed experimental and theoretical investigations leading to a quantitative description of the present guest-host interactions.

Nonlinear Molecular Electronic Spectroscopy via MCTDH Quantum Dynamics: From Exact to Approximate Expressions.

We present an accurate and efficient approach to computing the linear and nonlinear optical spectroscopy of a closed quantum system subject to impulsive interactions with an incident electromagnetic field. It incorporates the effect of ultrafast nonadiabatic dynamics by means of explicit numerical propagation of the nuclear wave packet. The fundamental expressions for the evaluation of first- and higher-order response functions are recast in a general form that can be used with any quantum dynamics code capable of computing the overlap of nuclear wave packets evolving in different states. Here we present the evaluation of these expressions with the multiconfiguration time-dependent Hartree (MCTDH) method. Application is made to pyrene, excited to its lowest bright excited state S2 which exhibits a sub-100-fs nonadiabatic decay to a dark state S1. The system is described by a linear vibronic coupling Hamiltonian, parametrized with multiconfiguration electronic structure methods. We show that the ultrafast nonadiabatic dynamics can have a remarkable effect on the spectral line shapes that goes beyond simple lifetime broadening. Furthermore, a widely employed approximate expression based on the time scale separation of dephasing and population relaxation is recast in the same theoretical framework. Application to pyrene shows the range of validity of such approximations.

Charged ultralong-range Rydberg trimers

We show that the recently observed class of long-range ion-Rydberg molecules can be divided into two families of states, which are characterised by their unique electronic structures resulting from the ion-induced admixture of quantum defect-split Rydberg $n$P states with different low-field seeking high-$l$ states. We predict that in both cases these diatomic molecular states can bind additional ground state atoms lying within the orbit of the Rydberg electron, thereby forming charged ultralong-range Rydberg molecules (ULRM) with binding energies similar to that of conventional non-polar ULRM. To demonstrate this, we consider a Rydberg atom interacting with a single ground state atom and an ion. The additional atom breaks the system's cylindrical symmetry, which leads to mixing between states that would otherwise be decoupled. The electronic structure is obtained using exact diagonalisation over a finite basis and the vibrational structure is determined using the Multi-Configuration Time-Dependent Hartree method. Due to the lobe-like structure of the electronic density, bound trimers with both linear and nonlinear geometrical configurations of the three nuclei are possible. The predicted trimer binding energies and excitation series are distinct enough from those of the ion-Rydberg dimer to be observed using current experimental techniques.

Optimal Mode Combination in the Multiconfiguration Time-Dependent Hartree Method through Multivariate Statistics: Factor Analysis and Hierarchical Clustering.

The multiconfiguration time-dependent Hartree (MCTDH) method and its multilayer extension (ML-MCTDH) are powerful algorithms for the efficient computation of nuclear quantum dynamics in high-dimensional systems. By providing time-dependent variational orbitals and an optimal choice of layered effective degrees of freedom, one is able to reduce the computational cost to an amenable number of configurations. However, choices related to selecting properly the mode grouping and tensor tree are strongly system dependent and, thus far, subjectively based on intuition and/or experience. Therefore, herein we detail a new protocol based on multivariate statistics─more specifically, factor analysis and hierarchical clustering─for a reliable and convenient guiding in the optimal design of such complex "system-of-systems" tensor-network decompositions. The advantages of employing the new algorithm and its applicability are tested on water and two floppy protonated water clusters with large amplitude motions.