bk-2016-1234.ch010

Downloaded by UNIV OF CALIFORNIA SAN DIEGO on December 12, 2016 | http://pubs.acs.org Publication Date (Web): December 1, 2016 | doi: 10.1021/bk-2016-1234.ch010

Chapter 10

Estimating Ground State Entanglement Entropy Using Path Integral Molecular Dynamics Dmitri Iouchtchenko and Pierre-Nicholas Roy* Department of Chemistry, University of Waterloo, Ontario, N2L 3G1, Canada *E-mail: [email protected].

We provide a perspective on the calculation of ground state entanglement properties of many-body bosonic systems with continuous degrees of freedom using path integral molecular dynamics sampling. To date, such ground state calculations have been performed using quantum Monte Carlo sampling and we show here that the Langevin Equation path integral ground state (LePIGS) dynamics sampling approach can be practical. A key challenge is the computation of reduced density matrices required for the estimation of the so-called Rényi entanglement entropy. We provide proof of concept results using a novel path integral molecular dynamics implementation of the permutation operator. This implementation allows us to efficiently estimate particle entanglement entropy in arbitrary continuum systems, among other quantities.

Introduction In recent years, the replica trick (1) has made it possible to efficiently estimate the entanglement entropy of ground state systems using path integral Monte Carlo (PIMC) (2). This is true both for systems on a lattice (2, 3) and in the continuum (4, 5). While PIMC simulations are extremely useful for various model systems, the need to devise updates for each system under study makes them impractical for molecular systems. The main advantage of the path integral molecular dynamics (PIMD) approach is that updates are automatically determined from classical © 2016 American Chemical Society Tanaka et al.; Recent Progress in Quantum Monte Carlo ACS Symposium Series; American Chemical Society: Washington, DC, 2016.


equations of motion: the paths are guided by a classical force. This allows for the simulation of arbitrary systems without any extra effort once the integrator is implemented. Since every update is a generic global update, there is a cost for this flexibility in the form of a decrease in computational efficiency (6). We show that it is possible to use the replica trick with PIMD to calculate ground state entanglement entropy. We do this by adding the ability to sample arbitrary path integral path connectivities to an existing PIMD integrator. We then benchmark our implementation using a simple harmonic model system.

Derivation of Estimator For a reduced ground state density operator

for subsystem A of total system AB, the replica trick gives access to the second Rényi entropy

We may express the trace in the position representation as

where the numerator has permuted labels, while the denominator does not. This is equivalent to the expectation of the permutation operator , which permutes the coordinates of particles in A, in a replicated system:

We do not have direct access to the ground state of an arbitrary . However, using path integral ground state (PIGS) (7, Hamiltonian 8), we express the ground state as

where is an arbitrary trial function. Following the standard procedure for path integrals, we introduce τ and P (which satisfy β = (P −1)τ with P an odd integer) and break the low-temperature imaginary-time propagator into high-temperature propagators:

Each of the high-temperature propagators may then be approximated using the Trotter factorization by a product of kinetic and potential factors. The resulting “paths” are best represented pictorially, as in Figure 1. 146 Tanaka et al.; Recent Progress in Quantum Monte Carlo ACS Symposium Series; American Chemical Society: Washington, DC, 2016.


Figure 1. Example paths for two interacting particles. Each path is composed of P beads. Wavy segments correspond to the kinetic energy operator and dashed segments correspond to the potential energy operator.

Using this graphical notation, but omitting the irrelevant beads, we may write eq. (3) as

where the action of the permutation operator is clearly evident. It turns out that by simulating an ensemble composed of paths for two modified copies of the system, we may estimate this quantity via a ratio of estimators. The necessary modification is the removal of links and interactions at the middle position for paths in A. We refer to this as “breaking” the paths, and it results in the distribution

Such a distribution corresponds to configurations in a sector different from the Z-sector which is typically employed for properties derived from the partition function (9). With this change to the expanded ensemble, the appropriate estimator becomes the ratio

Implementation Details To perform PIMD simulations, we use the Molecular Modelling Toolkit (MMTK) (10), which contains a Langevin equation path integral ground state (LePIGS) integrator (11). In usual applications of LePIGS, there is no distinction 147 Tanaka et al.; Recent Progress in Quantum Monte Carlo ACS Symposium Series; American Chemical Society: Washington, DC, 2016.


between paths and particles, since there is always a one-to-one mapping between them. To sample from a distribution with broken paths required modifying the integrator to deal with explicit paths, rather than paths implied by particles. This leads to some subtleties involving the interpretation of β, since the latter typically corresponds to both the simulation temperature and the length of each path. When paths may grow or shrink, their τ values remain unmodified, while their effective β values change. On the other hand, the simulation temperature must not vary during a simulation. Like the path integral ground state (PILE) integrator for finite temperature (12), the LePIGS integrator consists of the following five actions for each time step: 1. 2. 3. 4. 5.

apply the thermostat for half a time step; apply the force fields (interaction potentials and their gradients) for half a time step; propagate free particles for a full time step; apply the force fields for half a time step; and apply the thermostat for half a time step.

Path Connectivity Since force fields are implemented in Cartesian coordinates, the paths must be represented in Cartesian coordinates when the force fields are applied. However, the thermostat and free particle propagation act on the normal mode representation of the paths. This requires the integrator to convert between Cartesian and normal mode coordinates at each step, and this conversion contains one of the main differences between the PILE and LePIGS integrators. For both, we may write the spring energy for a single one-dimensional path with the P-dimensional vector of Cartesian coordinates q as

where

148 Tanaka et al.; Recent Progress in Quantum Monte Carlo ACS Symposium Series; American Chemical Society: Washington, DC, 2016.


for PILE, but

for LePIGS. The difference between these two matrices in the corner elements encodes the boundary conditions: PILE paths are closed (bead P connects to bead 1), while LePIGS paths are open. It turns out that the matrix in eq. (11) is diagonalized by the discrete Fourier transform, while the one in eq. (12) is diagonalized by the type-II discrete cosine transform. Both transforms are available in the FFTW package (13). By treating paths explicitly and recognizing this similarity in structure for the two integrators, we were able to combine them into a single integrator which handles both open and closed paths. This was a natural step along the way to implementing path breaking.

Force Field Scaling Another requirement for path breaking, and which was also a component in merging the integrators, was force field scaling. All the beads making up a single path in finite temperature simulations using the PILE integrator experience the same force due to the potential operator . In ground state simulations using LePIGS, the force acting on the two beads on the ends of a path is only half of the force acting on the other beads, due to the open boundary conditions. In path breaking, we wish to remove some force field terms altogether. Since MMTK handles each force field term internally using a PyFFEnergyTermObject struct, it was fairly straightforward to implement both force field modifications using the unified approach of force field scaling. We added the ability to scale every force field term independently by an arbitrary real value, which could be used to both reduce forces on the end beads of LePIGS paths and to remove some terms by scaling them to zero.

Reconnectors To make the implementation as general as possible, we introduced the possibility to obtain arbitrary path connectivity and force field scaling. This was done by providing the Cython (14) class PIReconnector with the following methods: 149 Tanaka et al.; Recent Progress in Quantum Monte Carlo ACS Symposium Series; American Chemical Society: Washington, DC, 2016.

•

double getScaling(self, Py_ssize_t t) to get the scaling of force field term t; void setScaling(self, Py_ssize_t t, double x) to set the scaling of force field term t to x; void openPath(self, Py_ssize_t p) to open the closed path p; void closePath(self, Py_ssize_t p) to close the open path p; void breakPath(self, Py_ssize_t p) to break the open path p into two paths; void joinPaths(self, Py_ssize_t p1, Py_ssize_t p2) to join two open paths p1 and p2 into a single path; Py_ssize_t[:] disassemblePath(self, Py_ssize_t p) to disassemble the path p into its constituent beads; and void assemblePath(self, Py_ssize_t[:] p) to assemble a path from a collection of beads p.

•


• • • • • •

An instance of a class which inherits from PIReconnector may be passed to the integrator, in which case its bint step(self, int step_num, double[:,:] x) method will be called with the step number and Cartesian coordinates of all beads at the end of each integration step. There exists a PIGSReconnector, which performs the necessary setup for a generic LePIGS simulation, but the end user is free to create a reconnector that suits their purpose. For example, the results presented in the following section were generated with the use of a custom reconnector, which inherits from PIGSReconnector and adds path breaking and force field scaling for particle entanglement entropy.

Model System Results To verify our implementation, we have applied it to the simplest non-trivial entangled continuum system: two harmonically coupled harmonic oscillators. The Hamiltonian is simply

for which we may obtain all relevant quantities analytically (4). Specifically, the second Rényi entropy for particle entanglement is

in which

and



The trial function employed in the simulations was the ground state of a nearly identical system, but with half the mass. This allowed for shorter simulations than would have been possible with a uniform trial function. However, it was sufficiently different from the exact ground state that β convergence still played a role in the analysis of the entanglement entropy. Convergence of the computed result for one choice of parameters is shown in Figure 2. This convergence proceeds as predicted by a separate numerically exact code. The computed and exact results for several values of the parameters are shown in Figure 3. These appear to match quite well, which suggests that we are capable of computing particle entanglement entropy in the ground state using path integral molecular dynamics.

Concluding Remarks We have presented a PIMD implementation of general path connectivity. Algorithms that allow for Feynman paths to be broken reconnected have been implemented in the MMTK software package. The implementation has been tested for LePIGS simulations of particles in a harmonic trap. Our proof-of-principle results indicate that we have a working and practical implementation for estimating particle entanglement entropy using PIMD. The logical next step is to apply this approach to more complex physical systems. One area of immediate interest is the study of particle entanglement within a low-temperature helium cluster doped with a CO2 molecule. Such clusters have been studied in detail (6, 15), and it would be interesting to see any correlations between the entanglement entropy and other properties, such as the superfluid response to a quantum mechanical probe (16), the so-called microscopic Andronikashvili experiment pioneered by Grebenev, Toennies, and Vilesov (17, 18). In this context, it will also be of great interest to study the relation between particle entanglement entropy and supefluidity in parahydrogen systems where freezing appears to begin at some cluster size consistent with the closing of a solvation shell (19), or where fast rotors such as CO lead to a non-trivial superfluid response (20). These applications to parahydrogen systems will also benefit from the use of tailored trial functions developed for LePIGS simulations of pure hydrogen clusters (21). The flexible implementation described above also opens the door for many other techniques which require sampling from different configuration sectors. For example, by dynamically changing the topology of the paths in a simulation, it will be possible to estimate the entanglement entropy for a continuum molecular system under a spatial bipartitioning (5). Additionally, we have already begun work on a PIMD implementation of the worm algorithm for efficient sampling of bosonic exchange (9). We finally note that our approach should be readily adaptable for other PIGS or variational path integral molecular dynamics implementations such as the Nosé–Hoover thermostat method proposed by Miura (22).



Figure 2. Convergence of entanglement entropy of coupled harmonic oscillators with τ for ωint/ω0 = 4. Dotted curve is the numerically exact result. Solid line is the exact result in the τ → 0 limit.

Figure 3. Entanglement entropy of coupled harmonic oscillators. Markers are the computed results with standard errors. The solid curve corresponds to the exact result. 152 Tanaka et al.; Recent Progress in Quantum Monte Carlo ACS Symposium Series; American Chemical Society: Washington, DC, 2016.

Acknowledgments We thank C. Herdman, A. Del Maestro, and R. Melko for stimulating discussions. We thank the Natural Sciences and Engineering Research Council of Canada and the Canada Foundation for Innovation for financial support.


References 1. 2.

3.

4.

5.

6.

7.

8. 9.

10. 11. 12.

13.

14.

Calabrese, P.; Cardy, J. Entanglement entropy and conformal field theory. J. Phys. A 2009, 42, 504005. Hastings, M. B.; González, I.; Kallin, A. B.; Melko, R. G. Measuring Renyi entanglement entropy in quantum Monte Carlo simulations. Phys. Rev. Lett. 2010, 104, 157201. Stéphan, J.-M.; Misguich, G.; Pasquier, V. Rényi entanglement entropies in quantum dimer models: from criticality to topological order. J. Stat. Mech.: Theory Exp. 2012, 2012, P02003. Herdman, C. M.; Roy, P.-N.; Melko, R. G.; Del Maestro, A. Particle entanglement in continuum many-body systems via quantum Monte Carlo. Phys. Rev. B 2014, 89, 140501. Herdman, C. M.; Inglis, S.; Roy, P.-N.; Melko, R.; Del Maestro, A. Pathintegral Monte Carlo method for Rényi entanglement entropies. Phys. Rev. E 2014, 90, 013308. Ing, C.; Hinsen, K.; Yang, J.; Zeng, T.; Li, H.; Roy, P.-N. A path-integral Langevin equation treatment of low-temperature doped helium clusters. J. Chem. Phys. 2012, 136, 224309. Hetenyi, B.; Rabani, E.; Berne, B. J. Path-integral diffusion Monte Carlo: Calculation of observables of many-body systems in the ground state. J. Chem. Phys. 1999, 110, 6143–6153. Sarsa, A.; Schmidt, K. E.; Magro, W. R. A path integral ground state method. J. Chem. Phys. 2000, 113, 1366–1371. Boninsegni, M.; Prokof’ev, N.; Svistunov, B. Worm algorithm and diagrammatic Monte Carlo: A new approach to continuous-space path integral Monte Carlo simulations. Phys. Rev. E 2006, 74, 036701. Hinsen, K. The molecular modeling toolkit: a new approach to molecular simulations. J. Comput. Chem. 2000, 21, 79–85. Constable, S.; Schmidt, M.; Ing, C.; Zeng, T.; Roy, P.-N. Langevin Equation Path Integral Ground State. J. Phys. Chem. A 2013, 117, 7461–7467. Ceriotti, M.; Parrinello, M.; Markland, T. E.; Manolopoulos, D. E. Efficient stochastic thermostatting of path integral molecular dynamics. J. Chem. Phys. 2010, 133, 124104. Frigo, M.; Johnson, S. The Design and Implementation of FFTW3. Proc. IEEE (Special issue on “Program Generation, Optimization, and Platform Adaptation”) 2005, 93, 216–231. Behnel, S.; Bradshaw, R.; Citro, C.; Dalcin, L.; Seljebotn, D.; Smith, K. Cython: The Best of Both Worlds. Comput. Sci. Eng. 2011, 13, 31–39. 153 Tanaka et al.; Recent Progress in Quantum Monte Carlo ACS Symposium Series; American Chemical Society: Washington, DC, 2016.


15. Li, H.; Blinov, N.; Roy, P.-N.; Le Roy, R. J. Path-integral Monte Carlo simulation of ν3 vibrational shifts for CO2 in (He)n clusters critically tests the He–CO2 potential energy surface. J. Chem. Phys. 2009, 130, 144305. 16. Zeng, T.; Roy, P.-N. Microscopic molecular superfluid response: theory and simulations. Rep. Prog. Phys. 2014, 77, 046601. 17. Grebenev, S.; Toennies, J. P.; Vilesov, A. F. Superfluidity within a small Helium-4 cluster: the microscopic Andronikashvili experiment. Science 1998, 279, 2083–2085. 18. Toennies, J. P.; Vilesov, A. F. Spectroscopy of atoms and molecules in liquid helium. Annu. Rev. Phys. Chem. 1998, 49, 1–41. 19. Li, H.; Le Roy, R. J.; Roy, P.-N.; McKellar, A. R. W. Molecular Superfluid: Nonclassical Rotations in Doped Para-Hydrogen Clusters. Phys. Rev. Lett. 2010, 105, 133401. 20. Raston, P. L.; Jäger, W.; Li, H.; Le Roy, R. J.; Roy, P.-N. Persistent molecular superfluid response in doped para-hydrogen clusters. Phys. Rev. Lett. 2012, 108, 253402. 21. Schmidt, M.; Constable, S.; Ing, C.; Roy, P.-N. Inclusion of trial functions in the Langevin equation path integral ground state method: Application to parahydrogen clusters and their isotopologues. J. Chem. Phys. 2014, 140, 234101. 22. Miura, S. A variational path integral molecular dynamics method applied to molecular vibrational fluctuations. Mol. Sim. 2012, 38, 378–383.


bk-2016-1234.ch010

Recommend Documents