PaPIM: A code for Quantum Time Correlation Functions

 

PaPIM code is a package to study the (quantum) properties of materials, and in particular time correlation functions, via the so-called mixed quantum-classical methods. In these schemes, quantum evolution is approximated by appropriately combining a set of classical trajectories for the system. Several quantum effects, for example, the possibility to find atoms in classically forbidden regions (tunneling), are reproduced at a manageable fraction of the cost of exact solutions.

The PaPIM module is a high-performance Fortran 90/95 MPI parallelized package for calculating system’s time-dependent observables. The code represents the current optimized assembly of the following modules:

  • PIM_wd and PIM_qcfmodules (described in deliverable D3.3) for exact quantum sampling of the Wigner phase space probability distribution function and the corresponding calculation of specific quantum correlation functions, respectively;
  • ClassMC module (described in D3.1) for Monte Carlo sampling of classical Maxwell-Boltzmann distribution and calculation of corresponding correlation-functions;
  • PotMod module (described in D3.1), a library for model potentials and interfaces to external codes for potential energy calculations used by the sampling modules. This module is currently being enhanced with an interface to couple PaPIM with the CP2K package for electronic structure calculations;
  • AuxMod module (described in D3.1) which provides a tailored set of MPI commands used for code parallelisation as well as input handling subroutines.

Practical application and exploitation of the code

The code has been extensively used for the calculation of the infrared absorption spectrum of CH5+ in the gas phase. [1] This highly flexible molecule is considered a standard benchmark of approximate quantum methods, and has experimental interest, for example, in the context of green chemistry. The calculations performed with PaPIM were used to benchmark both the PIM method for time-correlation functions [2] and to realize the code performance analysis.

Through collaborations the code is also currently employed by several groups in their study of: properties of H2 molecules in clathrates (materials for capture and storage of hydrogen and CO2 in energy applications (University College Dublin); infrared characterisation of molecules, and from it understand the effect that the environment has on their chemical properties, in the atmosphere (Université Pierre et Marie Curie); hydrogen at extreme pressures in the context of geophysical applications (Ecole Normale Supérieure Paris); new potentials to efficiently characterise the chemical reactivity of small water clusters, again with possible applications on the physics of the atmosphere in reactions related to greenhouse effect (University of Bochum).

More description of the code and its systematic tests are reported in the E-CAM deliverable D3.3.

 

[1] O. Asvany, P. K. P, B. Redlich, I. Hegemann, S. Schlemmer, D. Marx Understanding the infrared spectrum of bare CH5+ Science 309 (2005) 1219

[2] M. Monteferrante, S. Bonella, G. Ciccotti Quantum dynamical structure factor of liquid neon via a quasiclassical symmetrized method J. Chem. Phys. 138 (2013) 054118


Share

Spring shooting – A module for improving efficiency of transition path sampling

 

Transition path sampling is most efficient when paths are generated from the top of the free energy barrier. However, complex (biomolecular) activated processes, such as nucleation or protein binding/unbinding, can have asymmetric and peaked barriers. Using uniform selection on these type of processes will not be efficient, as it, on average, results in selected points that are not on the top of the barrier. Paths generated from these points have a low acceptance probability and accepted transition paths decorrelate slowly, resulting in a low overall efficiency. The Spring shooting module was developed to increase the efficiency of path sampling of these types of barriers, without any prior knowledge of the barrier shape. The spring shooting algorithm uses a shooting point selector that is biased with a spring potential. This bias pulls the selection of points towards the transition state at the top of the barrier. The paths that are generated from points selected by this biased selector therefore have an increased acceptance probability and the decorrelation between accepted transition paths is also increased. This results in a higher overall efficiency. The spring shooting algorithm is described in more detail in a paper by Brotzakis and Bolhuis. [1]  This module was developed during the ESDW on classical molecular dynamics held in Amsterdam.

 

[1] Z. F. Brotzakis, P. G. Bolhuis A one-way shooting algorithm for transition path sampling of asymmetric barriers J. Chem. Phys. 145 (2016) 164112

Share

Symmetry Adapted Wannier Functions – a Component of the Wannier90

 

Symmetry Adapted Wannier Functions is a module within Wannier90 which is devoted to the construction of Wannier function (WF) with a given symmetry. The procedure implemented in this module enables one to control the symmetry and center of the WFs and also simplifies the minimisation of the spread functional under these symmetry constraints.

This module is part of the nine modules reported in Deliverable D2.3 which together deal with the implementation of symmetry adapted WFs, to improve the symmetery of the WFs and related electronic-structure quantities, such as band structure and density of states; improvements in the interpolation of band structures, developments in the selection of the k-point mesh to increase accuracy, ability of performing non-collinear spin calculations as well as interface layer modules to tight-binding codes.

Starting from an E-CAM ESDW3 in San Sebastian organised by the Wannier90 developers, a set of nine modules were produced to meet the desire of the electronic-structure community to extend the use of WFs, and in particular of Maximally Localised Wannier Functions (MLWFs), to a broader class of physical and chemical problems by adding new functionality to the Wannier90 code.

All modules are accessible through the Wannier90 code, which in turn is interfaced with the all the most popular DFT codes. Wannier90 is used as a postprocessing tool. Therefore, the end users of electronic-structure codes, such as DFT, Tight Binding and Quantum Monte Carlo codes, that are interfaced with these modules via Wannier90, will benefit from the functionalities they provide, e.g. WFs with improved symmetry, spin-orbit calculations etc., and they can focus on developing new ideas, and new science without needing to rewrite functionalities that are already established.

Practical application and exploitation of the code

Wannier functions are an important class of functions which enable one to obtain a real-space picture of the electronic structure of a system. They provide an insightful chemical analysis of the nature of bonding, and chemical reaction in condensed-matter physics, similar to the role played by localised molecular orbitals in chemistry. They are also a powerful tool in the study of dielectric properties via the modern theory of polarisation. In the condensed-matter community WFs are employed in the construction of model Hamiltonians for, e.g., correlated-electron and magnetic systems (to study new quantum phases of matter) and are used as building blocks in first-principles Tight Binding Hamiltonians, where chemically accurate Hamiltonians are constructed directly on the Wannier basis, rather than fitted or inferred from macroscopic considerations. [1]

Wannier90 [2] is a program that, for a given system, generates the Wannier functions with minimum spatial spreads, known as MLWFs, among the class of all possible WFs. The locality of MLWFs can be exploited to compute, among other things, band-structure, density of states and Fermi surfaces at modest computational cost.

The developed modules have been used to study the properties of strongly correlated materials and to assess the quality of high-level quantum methods. [3]

 

[1] A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, N. Marzari wannier90: A tool for obtaining maximally-localised wannier functions Comput. Phys. Commun 178 (2008) 685

[2] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, D. Vanderbilt Maximally localized wannier functions: Theory and applications Rev. Mod. Phys. 84 (2012) 1419

[3] L. Boehnke, F. Nilsson, F. Aryasetiawan, P. Werner When strong correlations become weak: Consistent merging of GW and DMFT Phys. Rev. B 94 (2016) 201106

Share

PIM_wd: Module for sampling of the quantum Wigner distribution

 

The PIM_wd module implements the exact quantum Wigner probability distribution function sampling algorithm of the Phase Integration Method [1], and is the main subroutine for the quantum correlation function calculations in the PaPIM code. The module samples the thermal Wigner density using a generalised Monte Carlo scheme for sampling phase space points. The scheme combines the Penalty [2] and Kennedy [3] algorithms to sample noisy probability densities. This is necessary because the estimator of the quantum thermal density is not known analytically but must be computed via a statistical average affected by uncertainty. The sampled points are the basis for the calculation of time-independent and time-dependent system observables.

The module was developed as the main component of the PaPIM code, but also as a standalone subroutine that can be easily implemented in other methods (e.g. the whole family of so-called linearised approximations of quantum dynamics) for which phase space sampling of the exact quantum Wigner probability distribution is required. Because the Phase Integration Method samples a set of independent phase space points, independent instances of the PIM_wd module can be run in parallel in order to parallelise the phase space sampling. In the PaPIM package, the parallelisation is accomplished using MPI, which has proved to provide good scalability of the PaPIM code. The module will also be adapted for HTC capabilities.

Practical application and exploitation of the code

The code has been extensively used for the calculation of the infrared absorption spectrum of CH5+ in the gas phase. [4] This highly flexible molecule is considered a standard benchmark of approximate quantum methods, and has experimental interest, for example, in the context of green chemistry.

This module is part of the modules in deliverable D3.3 which were developed during the E-CAM ESDW7.

 

[1] M. Monteferrante, S. Bonella, G. Ciccotti Quantum dynamical structure factor of liquid neon via a quasiclassical symmetrized method J. Chem. Phys. 138 (2013) 054118

[2] D. M. Ceperley, M. Dewing The penalty method for random walks with uncertain energies J. Chem. Phys. 110 (1999) 9812

[3] A. D. Kennedy, J. Kuti Noise without Noise: A New Monte Carlo Method Phys. Rev. Lett. 54 (1985) 2473

[4] O. Asvany, P. K. P, B. Redlich, I. Hegemann, S. Schlemmer, D. Marx Understanding the infrared spectrum of bare CH5+ Science 309 (2005) 1219

Share

Coarse-Graining module, a Component of the Hierarchical Equilibration Strategy for Polymer Melts

To study the properties of polymer melts by numerical simulations, equilibrated configurations must be prepared. However, the relaxation time for high molecular weight polymer melts is huge and increases, according to reptation theory, with the third power of the molecular weight. Hence, an effective method for decreasing the equilibration time is required. The hierarchical strategy pioneered in Ref. [1] is a particularly suitable way to do this. The present module provides a part of that method.

To decrease the relaxation time, microscopic monomers are coarse-grained (CG) by mapping each subchain with N_{b} monomers onto a soft blob. The CG system is then characterized by a much lower molecular weight and thus is equilibrated quickly. The present module provides a python script which performs this coarse-graining procedure. The implementation details can be seen in the module’s documentation on our software Library here. This module is part of a set of codes that together implement the Hierarchical Equilibration strategy of Ref. [1], in the ESPResSO++ [2] (for the complete list of modules, see here under ESPResSO++).

 

Practical application and exploitation of the code

The development of a multiscale method for polymer blends and block copolymers is fundamentally new and needs to be based on first-principles theory. This is therefore an intellectual challenge in its own right. Furthermore, this paves the way to analyze the physical properties of novel composite materials that attract the attention of industrial companies. Such materials may be promising ingredients of new products like e.g. efficient and environment-friendly car tires. The implementation of the Hierarchical Equilibration strategy in the ESPResSO++ package is a step towards achieving this goal. In particular,  the practical application of this strategy is the E-CAM pilot project in collaboration with Michelin aimed at studying the Rheological Properties of New Composite Materials.

E-CAM deliverables D4.2 and D4.3 contain more information on the suite of programs developed under this pilot project.

 

[1] Zhang, G., Moreira, L. A., Stuehn, T., Daoulas, K. C., and Kremer, K., Equilibration of High Molecular Weight Polymer Melts: A Hierarchical Strategy, ACS Macro Lett., 3, 198-203 (2014)

[2] ESPResSo++ is the “Extensible Software Package for Research in Soft Matter based upon C++”, a general-purpose simulation package for soft-matter research, mainly developed at the Max Planck Institute for Polymer Research Mainz. It is freely available under the GNU Public License. http://www.espresso-pp.de/

Share

Contact Map – a package for analyzing and exploring contacts, from a trajectory generated by MD

 

Contacts can be an important tool for defining (meta)stable states in processes involving biomolecules. For example, an analysis of contacts can be particularly useful when defining bound states during a binding processes between proteins, DNA, and small molecules (such as potential drugs).

The contacts analyzed by the contact_map package can be either intermolecular or intramolecular, and can be analyzed on a residue-residue basis or an atom-atom basis.

This package makes it very easy to answer questions like:

  • What contacts are present in a trajectory?
  • Which contacts are most common in a trajectory?
  • What is the difference between the frequency of contacts in one trajectory and another? (Or with a specific frame, such as a PDB entry.)
  • For a particular residue-residue contact pair of interest, which atoms are most frequently in contact?

It also facilitates visualization of the contact matrix, with colors representing the fraction of trajectory time that the contact was present. Full documentation available at http://contact-map.readthedocs.io/.

Information about software installation, testing and a link to the source code, can be found in our E-CAM software Library here.

Practical application and exploitation of the code

The practical application of this software module is the pilot project in collaboration with BiKi Technologies on “Binding Kinetics“, sustained by an E-CAM postdoctoral researcher at University of Amsterdam.  The project aims at investigating the binding/unbinding of a selective reversible inhibitor for protein GSK3β.

Contacts between a ligand and a protein are an excellent way to characterize “hotspots” – states where the ligand stays for a significant amount of time, but not nearly as long as in the final binding pocket. These hotspots are metastable states in path sampling, and should be treated with a multiple state approach. Therefore, attempting to identify those states would be a necessarily preliminary step to prepare the path sampling simulation.

Other more general applications to this module include protein-protein aggregation or DNA-protein binding, as well as large scale conformational changes in biomolecules, such as protein folding.

 

Share

LocConQubit, a module for the construction of controlled pulses on isolated qubit systems using the Local Control Theory

 

The LocConQubit module implements the Local Control Theory[1,2], an algorithm for on-the-fly construction of a time-dependent potential that drives the evolution of a Hamiltonian towards one of its eigenstates. The algorithm is applicable to any Hamiltonian that is separable into a time-dependent and into a time-independent part, where the first part is directly incorporated into the algorithm, while the latter defines the basis of system states from which a designated target state is selected. States with vanishing interaction elements cannot be treated with the aforementioned algorithm. The algorithm is fine-tuned by the user with a single parameter in order to assure physical range of the generated time-dependent potential. This free parameter can be time-dependent while certain constrains in pulse generation can be directly incorporated into the algorithm. The module is accompanied with subroutines for pulse frequency analysis, post-processing, fidelity calculation and visualization of pulses and system evolution. The module is written in Python 3 programming language and is an addition to the open source QuTiP software package. The module uses the OpenMP functionalities available in QuTiP to parallelize the calculation of the pulse fidelity in order to search more efficiently for an optimal control pulse.

Additional module documentation, which includes background information on the Local Control Theory, information about software installation and testing and a link to the source code, can be found in our E-CAM software Library here

Practical application and exploitation of the code

The practical application of this software module is the pilot project with IBM on “Quantum Computing” sustained by an E-CAM postdoctoral researcher at École Polytechnique Fédérale de Lausanne (EPFL).

This module enables to construct more efficient control pulses for superconducting transmon qubits coupled to a single tunable coupler whose energy is controlled with an external electromagnetic pulse. By properly modulating the energy of the tunable coupler with an external control pulse, the coupler operates as a quantum logic gate between coupled qubits. To improve gate performance and thus overall performance of quantum computers, pulses are tailored to make gate operations faster while maintaining at the same time the highest possible fidelity. The Local Control Theory was applied to these systems to generate efficient state preparation pulses which transfer populations completely from one qubit state to the other, as well as pulses for the SWAP gates which completely exchange quantum states between two qubits. A set of pulses capable of transferring populations with a full fidelity to designated target states was generated and, by post-processing this set, an optimal set of pulses for experimental implementation was obtained. This set is currently being tested at IBM. In parallel, capabilities as well as limits of the Local Control theory to manipulate such systems have been investigated in detail. Results of this work are going to be published in two scientific papers. In addition, the current OpenMP parallelization will be upgraded with a more advance parallelization scheme that will enable more efficient utilization of \acs{HPC} resources and an easier implementation of parallelized optimization techniques.

 

[1] B. F. E. Curchod, T. J. Penfold, U. Rothlisberger and I. Tavernelli, Local control theory in trajectory-based nonadiabatic dynamics, Phys. Rev. A, vol. 84, p. 042507, 2011. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevA.84.042507

[2] V. Engel, C. Meier, and D. J. Tannor, Local Control Theory: Recent Applications to Energy and Particle Transfer Processes in Molecules, John Wiley Sons, Inc., 2009, pp. 29–101. [Online]. Available: http: //dx.doi.org/10.1002/9780470431917.ch2

Share

GRASP Sampling – a module to build a representative data set for a fitting procedure

GRASP_sampling performs a stratified sampling of the configurations, described by vectors, of a system to build a representative training set in a fitting procedure. Given a list of candidate configurations, and selected the size (N) of the training set required, the module executes the combinatorial optimization that maximizes the following dissimilarity score (DS) among the elements of the training set:

../../../_images/dissimilarity_score.png

In this formula, the j-th configuration in the sum is the j-th nearest one to the l-th configuration and dij is the Euclidean distance between the l-th and j-th configurations. M is the number of the nearest configurations considered in the score. The exponential weight makes the score near independent from the particular value of M, if it is larger than 4-6.

The combinatorial optimization that maximizes the dissimilarity score is performed using the greedy randomized adaptive search procedure[1]  (GRASP) algorithm. A stratified sampling can be performed without a combinatorial optimization using classical statistical techniques (for example Latin hypercube sampling), the GRASP sampling becomes useful when the selection is restricted to a predeterminated set of configurations, generated or sampled with specific internal constrains. This is the case of the molecular configurations generated in a molecular dynamics simulation.

The complete module documentation, including a link to the source code, can be found in our repository here

Motivation and exploitation

The application of the GRASP algorithm to perform a stratified sampling is described in a recent publication [2] by the E-CAM partners at Scuola Normale Superiore (SNS), that we previously reported here.

The motivation behind this software module is the pilot project with industry “Quantum Mechanical Parameterisation of Metal Ions in Proteins” sustained by an E-CAM postdoctoral researcher from SNS.

 

[1] Feo, T. A.; Resende, M. G. Greedy randomized adaptive search procedures. J. Glob. Optim. 1995, 6, 109−133

[2] Francesco Fracchia, Gianluca Del Frate, Giordano Mancini, Walter Rocchia, and Vincenzo Barone, Force Field Parametrization of Metal Ions from Statistical Learning Techniques, J. Chem. Theory Comput. 2018, 14, 255−273

Share

Geomoltools: A set of software modules to easily manipulate molecular geometries

Geomoltools is a set of eight pre- and post-treatment Fortran codes that can be used to easily manipulate molecular geometries, allowing to minimize the average energy obtained for a range of internuclear distances for the dimers of each element, and decrease the computational cost of a DFT calculation.

The set of codes are:

  • mol2xyz: converts a .mol file into an ordered .xyz file
  • pastemol: joins two .xyz files
  • movemol: translates and aligns the molecule with some predefined axes
  • stackmol: generates (manually or randomly) different stacking arrangements between two molecules
  • geodiff: compares the internal coordinates of two molecules
  • xyz2zmt_s: converts the cartesian coordinates contained in a .xyz file into Z-matrix (2 possible formats)
  • zmt2xyz_s: converts a Z-matrix (from 2 possible formats) into cartesian coordinates
  • ucubcellgen: calculates the vectors of a unit cell given some atomic coordinates.

Modules source codes can be found here.  For a detailed explanation of the main programs, please have a look to this file. A complete tutorial on how to use the different codes from the package Geomoltools in order to manipulate (rotate, translate, join, pack, convert, etc.) molecular geometries, can be found at this address.

Motivation and exploitation

These modules have been used to study the stacking arrangements of acceptor:donor molecules for organic photovolatics polymers by high-throughput computation with the SIESTA code. This set of codes are available under the GNU General Public License (GPL) version 2.

Share

Path density for OpenPathSampling

Module path density implements path density calculations for the OpenPathSampling (OPS) package, including a generic multidimensional sparse histogram, and plotting functions for the two-dimensional case. Path density plots provide a way to visualize kinetic information obtained from path sampling, such as the mechanism of a rare event. In addition, the code in this module can also be used to visualize thermodynamic information such as free energy landscapes.

This module has been incorporated into the core of OPS, an open-source Python package for path sampling that wraps around other classical Molecular Dynamics (MD) codes [1]. An easy-to-read article on the use of path sampling methods to study rare events, and the role of the OPS package to performing these simulations can be found here.

At first glance, a typical path density plot may appear similar to a two-dimensional free energy landscape plot. They are both “heatmap”-type plots, plotting a two-dimensional histogram in some pair of collective variables. However, path density differs from free energy in several important respects:

  • A path density plot is histogrammed according to the number of paths, not the number of configurations. So if a cell is visited more than once during a path, it still only gets counted once.
  • A path density plot may interpolate across cells that the path jumps over. This is because it is assumed that the input must actually be continuous.

These differences can prevent metastable regions from overwhelming the transition regions in the plot. When looking at mechanisms, the path density is a more useful tool than the raw configurational probability.

Module documentation can be found here, including a link to the source code. This and other software modules for studying the thermodynamics and kinetics of rare events where recently documented in deliverable D1.2.: Classical MD E-CAM modules I, available here.

Motivation and exploitation

The path density is one of the most important tools for visualizing mechanisms, and is often one of the first things to analyze in order to draw scientific conclusions about the mechanism from transition path sampling simulations. This module was used to illustrate the differences between dynamics of the wild-type and oncogenic mutant forms of KRas, as part of one student’s master’s thesis and another student’s bachelor’s thesis at the University of Amsterdam. Results from those projects are currently in preparation for publication [2].

 

[1] Jan-Hendrik Prinz, David W.H. Swenson, Peter G. Bolhuis, and John D. Chodera. OpenPathSampling: A Python framework for path sampling simulations. I. Introduction and usage. In prep.
[2] Sander Roet, Ferry Hooft, Peter G. Bolhuis, David W.H. Swenson, and Jocelyne Vreede. Simulating the dynamics of oncogenic and wild-type KRas. In prep.

Share