MoRiBS-PIMC: A program to simulate molecular rotors in bosonic solvents using path-integral Monte Carlo

Abstract

We provide the source code of our in-house program MoRiBS-PIMC. This program was developed to simulate rigid molecules rotating in bosonic clusters composed of helium atoms, parahydrogen molecules or any other bosonic point solvent particles. The program can be employed to obtain superfluid response, structural and energetic properties as well as imaginary time correlation functions of dipole operators. These quantities can be used to interpret and predict the results of spectroscopic Andronikashvili experiments. The software is based on the latest advances in the simulation of the quantum rotation of non-linear rigid rotors and in the sampling of bosonic permutations. The program has been parallelized to improve its performance and new techniques have been implemented to obtain symmetry-adapted simulation results. The usage and robustness of the program is demonstrated with some illustrative examples.

Program summary

Program title: MoRiBS-PIMC

Catalogue identifier: AFAD_v1_0

Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AFAD_v1_0.html

Program obtainable from: CPC Program Library, Queen’s University, Belfast, N. Ireland

Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html

No. of lines in distributed program, including test data, etc.: 4161981

No. of bytes in distributed program, including test data, etc.: 25399583

Distribution format: tar.gz

Programming language: C++ and Fortran77.

Computer: Any computers with OpenMP libraries installed.

Operating system: Linux.

RAM: On the order of hundreds of MBytes and system-dependent. See Section 5 for examples.

Supplementary material: The full input data for the examples is available for download.

Classification: 16.02, 16.10, 16.13, 23.

External routines: Lapack routine, dsyev, is needed in a subprogram asymrho.f. Dynamic library for quadruple precision floating number arithmetics is needed for subprograms asymrho.f and symrho.f and linden.f to tabulate rotational propagator, energy and heat capacity estimator. The source codes of these programs are distributed with the main program.

Nature of problem: Many body quantum physics, rigid-body rotation, atomic and molecular clusters, bosonic exchange and superfluidity

Solution method: Path-integral Monte Carlo simulations and worm algorithm for permutation sampling

Restrictions: See Section 1 of the manuscript for details.

Unusual features: This is the first program that combines the worm algorithm for bosonic permutation sampling and general rigid-body rotation for asymmetric tops. The program has extensive applications in quantum molecular dynamics and it has been parallelized via OpenMP.

Additional comments: The program can be used to interpret and predict spectroscopic Andronikashvili experiments conducted on microscopic superfluids.

Running time: System dependent. See Section 5 for examples.

Introduction

Classical molecular dynamics simulations have been a mainstay in the molecular physics community for the study of a variety of systems. The classical approximation certainly produces reliable results at high enough temperatures for which the nuclear quantum effects are either mild or washed out by thermal fluctuations. Methodologies that explicitly account for nuclear quantum effects have become increasingly sought after with the advances in cryogenic physics. These methodologies are especially useful for atoms and molecules with light-massed nuclei such as hydrogen and helium due to their quantum behaviour at temperatures up to a few Kelvin. The Path-Integral Monte Carlo (PIMC) approach is one such methodology and is the framework for the present software package.

The path-integral formalism was first proposed by Feynman  [1], [2] as a general picture to describe quantum statistical mechanics. It represents each particle by a polymer consisting of $P$ beads, with each bead corresponding to a replica of the system. The particle to ring polymer transformation embodies the quantum smearing of an atom in space and is also the basis for the isomorphism with the classical simulation of polymers. PIMC is a method to simulate quantum particles by sampling the configurations of their corresponding “polymers”  [3], [4]. The isomorphism makes it possible to perform quantum simulations at the cost of classical mechanical simulations and allows for PIMC to be computationally feasible for large-scale quantum systems. Of key importance is the fact that the exchange of indistinguishable particles is conveniently represented by the connectivities of the corresponding polymers. These two features represent why PIMC is a highly convenient method for quantum simulations. At least for systems of bosons, PIMC is the only presently known method capable of furnishing exact numerical estimates in principle of physical quantities, including the superfluid density and the condensate fraction   [5].

In contrast to many quantum mechanical phenomena that occur in the atomic scale, superfluidity is famous for being a macroscopic quantum effect. This phenomenon can for instance be observed in a container filled with liquid ⁴He. For this reason, the study of superfluidity has occupied a special place in the physics community. How large does a quantum system have to be to start exhibiting some of the signatures of superfluid behaviour? A breakthrough in this direction was the Gottingen experiment conducted by Grebenev et al.  [6]. They carried out high-resolution infrared (IR) spectroscopic measurements on helium droplets doped with an OCS molecule as a chromophore. They found that when OCS is surrounded by as few as 60 ⁴He atoms, the IR spectrum exhibits narrow rotational lines with a rotational constant greater than that of the bare OCS molecule. They attributed the reduced rotational constant (increased moment of inertia) to the adiabatic following of the normal component of the ⁴He to the OCS rotation. The sharp rotational lines were attributed to the coherent decoupling (zero viscosity) between the normal and superfluid components. Evidently, nano-scale systems can still display superfluid behaviour. At the macroscopic scale, helium is the only known superfluid because other substances solidify at the low temperatures required for this phenomenon to take place. Nano-clusters of other substances may also exhibit superfluidity due to the smaller number of neighbours (weaker interaction) from dimension reduction. Clusters composed of para-hydrogen ($p{H}_2$) molecules are good candidates for the observation of superfluidity. In 2000, the same research group published a similar study on OCS${\left(p{H}_2\right)}_{14\text{–}16}$ clusters and they observed evidence of $p{H}_2$  superfluidity  [7] as manifested by the absence of a Q-branch in the rovibrational spectrum. So far, $p{H}_2$  is the only microscopic superfluid known to us other than helium.

These two groundbreaking experiments have motivated a number of subsequent studies of microscopic superfluidity. These experiments have been referred to as microscopic (or spectroscopic) Andronikashvili experiment due to their connection to the historical Andronikashvili experiment where one measures the superfluid fraction of macroscopic systems  [8], [9]. The essence of these experiments is to dope a molecular chromophore in a nano-scale cluster or droplet and study the rotational structure in the spectrum of the complex. It is often desirable or even necessary to conduct a high-level theoretical simulation for the complex due to the complexity of the spectrum of a weakly bound system and employ the simulation results to interpret or even predict their superfluid behaviour.

In this contribution, we provide the source code of our in-house PIMC program that has been successfully used to simulate molecular rotors in bosonic solvents in the past decade. The italic letters form the acronym of the first part of the name, MoRiBS, while the second half (PIMC) indicates the simulation methodology. The main structure of the program was developed to handle rotation of one linear rigid rotor in a bath of one type of bosons and was written in C++ by one of us (NB)  [10], [11], [12], [13], [14], [15], [16], [17], [18]. Other authors (TZ, GG and HL) made significant contributions to parallelize the program and to extend its applicability to handle an arbitrary number of one type of generic asymmetric rigid rotors. Some of these latter contributions are programmed as Fortran77 subroutines and they are called by the C++ main program. Since the beginning, the development of MoRiBS-PIMC has been under the supervision of PNR. The functionality of the program is to perform a quantum mechanical simulation for an $NVT$ ensemble using the PIMC method. Equilibrium properties such as energies, structures, superfluid fractions, and effective moments of inertia (effective rotational constants) may be obtained via the analysis of simulation runs. The program is applicable to the following types of systems:

• An arbitrary number of one type of point-like particles such as ⁴He clusters that may be treated as either boltzmannons or bosons (including bosonic exchange symmetry). Fermions are beyond the scope of the present program;

• An arbitrary number of one type of non-linear rigid rotors such as H₂O clusters where the rotors are treated as boltzmannons and the exchange effect between rotors is ignored;

• One linear rigid rotor such as CO₂, surrounded by an arbitrary number of one type of point-like boltzmannons or bosons;

• An arbitrary number of one type of non-linear rigid rotors surrounded by an arbitrary number of one type of point-like boltzmannons or bosons.

The present version of our program does not support periodic boundary conditions but we aim to support this feature in future developments.

The connection between our program and the spectroscopic Andronikashvili experiment is clear. Compared to other existent PIMC programs, MoRiBS-PIMC features the combination of bosonic exchange of point-like particles and asymmetric top rotation. The rest of the paper is organized as follows. Section  2 provides a brief account of the PIMC methodology that is focused on the treatment of rigid rotors, bosonic exchange, and the calculation of superfluid fractions. Section  3 details how the methodology is implemented in the computer program along with the implementations of bosonic permutation sampling and symmetry-adapted configuration sampling in simulations. Section  4 gives a brief introduction to the input and output files and Section  5 provides simulation results of several examples. Finally, Section  6 provides some concluding remarks.