Towards Efficient Direct Semiclassical Molecular Dynamics for Complex Molecular Systems

Chemical processes are intrinsically quantum mechanical and quantum effects cannot be excluded a priori. Classical dynamics that use fitted force fields have been routinely applied to complex molecular systems. But since the force fields used in classical dynamics are tuned to fit experimental and/or electronic structure data, the harmonic potential approximation and the negligibility of quantum effects are artificially and ad hoc compensated. Also, fitting atomic forces is usually a trade-off between the desired accuracy and the human and computational effort required to construct them, and it is often biased by the functional forms chosen. Thus, it can happen that the force field is not transferable, i.e. it cannot be applied a priori to other molecular systems. In addition, force fields do not account for bond dissociation or excited vibrational processes, due to the harmonic approximation.

To bypass these force field limitations, an alternative is the direct dynamics (on-the-fly) approach, with which the nuclear classical dynamics is coupled with atomic forces calculated from quantum mechanical electronic structure theory. Direct semiclassical molecular dynamics employs thousands of direct dynamics trajectories to calculate the Feynman Path Integral propagator, and reproduces quantitative quantum effects with errors often smaller than 1%, making it a very promising tool for including quantum effects for complex molecular systems. Direct semiclassical dynamics incurs much lower computation cost than purely quantum dynamics, but still calls for substantial reduction of computation cost for application to complex and interesting molecular systems on large HPC machines.

The high computation cost of direct semiclassical dynamics comes from two sources. One is the large number of trajectories needed. The other is the enormous computation cost to calculate a single trajectory. In this talk, we present our efforts in containing computation costs from these two sources in order to make direct semi-classical dynamics feasible on modern HPC systems.

A single trajectory of a direct semiclassical dynamics simulation may take days to weeks on a powerful multi-core processor. For instance, our on-going study of 10-atom glycine with the B3LYP/6-31G** electronic structure theory takes about 11.5 days on two quad-core Intel Xeon 2.26GHz processors (8 cores total) for a trajectory of 5000 time steps. To reduce the single trajectory calculation time, we developed a mathematical method to utilize directional data buried in previously calculated quantum data for future time steps, thereby reducing the expensive quantum mechanical electronic structure calculations. With the new method, we are able to reduce the computation time of a 5000-step trajectory to about 2 days with almost the same accuracy.

A simulation study for glycine requires hundreds of thousands to even millions of trajectories when a usual semiclassical method is used. To reduce this requirement, we have developed a semiclassical algorithm that can calculate molecular spectra and vibrational eigenfunctions with only a few trajectories and faithfully reproduce the result from a calculation with thousands of trajectories. This has allowed us to perform first direct ab initio semiclassical dynamics, and then determine the power spectrum and anharmonic vibrational frequencies, for a molecule as large as glycine. Along this line, we have developed a method to investigate more and more complex systems. Recently we reproduced quantum mechanical tunneling in the resonating umbrella inversion for ammonia using substantially fewer trajectories. For the on-going investigation for 10-atom glycine, we estimate that our developed code needs only hundreds of trajectories to reach fidelity comparable to hundreds of thousands of trajectories that a usual semiclassical method would require.

Since the computations of different trajectories are embarrassingly parallel, our trajectory-reduction semiclassical method combined with our directional information utilization technique for reducing the single-trajectory computation time, promises a feasible HPC solution for addressing complex molecular problems. Only hundreds of multi-core nodes for two days, or tens of multi-core nodes for a month, are required to complete the direct semiclassical dynamics simulation of glycine.