Porting of electrostatics to the GPU version of DL_MESO_DPD

The porting of DL_MESO_DPD [1,2] to graphic cards (GPUs) was reported in deliverable D4.2 of E-CAM[3] (for a single GPU) and deliverable D4.3 [4] (for multiple GPUs) (Figure 1), and has now been extended to include electrostatics, with two alternative schemes as explained below. This work was recently reported on deliverable D4.4[5].

Figure 1: DL_MESO strong scaling results on PizDaint, obtained using 1.8 billion particles for 256 to 2048 GPUs. Results show very good scaling, with efficiency always above 89% for 2048 GPUs.

To allow Dissipative Particle Dynamics (DPD) methods to treat systems with electrically charged particles, several approaches have been proposed in the literature, mostly based on the Ewald summation method [6]. The DL_MESO_DPD  code includes Standard Ewald and Smooth Particle Mesh Ewald (SPME) methods (in version 2.7, released in December 2018). Accordingly, here the same methods are implemented for the single-GPU version of the code. Continue reading…


CTMQC, a module for excited-state nonadiabatic dynamics


CTMQC is a module for excited-state nonadiabatic dynamics. It is used to simulate the coupled dynamics of electrons and nuclei (ideally in gas phase molecular systems) in response to, for instance, an initial electronic excitation.

The CTMQC module is based on the coupled-trajectory mixed quantum-classical (CT-MQC) algorithm [1,2] that has been derived starting from the evolution equations in the framework the exact factorization of the electron-nuclear wavefunction [3,4,5]. The CTMQC algorithm belongs to the family of quantum-classical methods, as the time evolution of the nuclear degrees of freedom is treated within the classical approximation, whereas electronic dynamics is treated fully quantum mechanically. Basically, the nuclei evolve as point particles, following classical trajectories, while the electrons generate the potential inducing such time evolution.

In its current implementation (used in Refs. [6,7]), the module cannot deal with arbitrary nuclear dimensions, but it is restricted to treat up to 3-dimensional problems, which gives the possibility to compare quantum-classical results easily and directly with quantum wavepacket dynamics. CTMQC has been analyzed and benchmarked against exact propagation results on typical low-dimensional model systems [1,2,6,7], and applied for the simulation of the photo-initiated ring-opening process of Oxirane [8]. For this study, CTMQC has been implemented in a developer version of the CPMD electronic structure package based on time-dependent density functional theory. Concerning electronic input properties, the CTMQC module requires a grid representation of the adiabatic potential energy surfaces and of the nonadiabatic coupling vectors, since the electronic dynamics is represented and solved in the adiabatic basis.

This feature allows the algorithm to be easily adaptable, in the current form, to any quantum chemistry electronic structure package. The number of electronic states to be included is not limited and can be specified as input.

Practical application and exploitation of the code
The purpose of the module is to familiarize the user with a new simulation technique, i.e., the CTMQC method, for treating problems where electronic excited states are populated during the molecular dynamics. Photo-activated ultrafast processes are typical situations in which an approach like CTMQC can be used to predict molecular properties, like structures, quantum yields, or quantum coherence.
The module is designed to apply the CTMQC procedure to one-, two-, and three-dimensional model systems where an arbitrary number of electronic states are coupled via the nuclear dynamics. Tully model systems [9] are within the class of problems that can be treated by the module, as well as a wide class of multidimensional problems involving, for instance, ultrafast radiationless relaxation of photo-excited molecules [10] through conical intersections.


Software documentation can be found in our E-CAM software Library here.


[1] S. K. Min, F. Agostini, E. K. U. Gross Coupled-trajectory quantum-classical approach to electronic decoherence in nonadiabatic processes Phys. Rev. Lett. 115 (2015) 073001
[2] F. Agostini, S. K. Min, A. Abedi, E. K. U. Gross Quantum-classical nonadiabatic dynamics: Coupled- vs independent-trajectory methods J. Chem. Theory Comput. 12 (2016) 2127
[3] A. Abedi, N. T. Maitra, E. K. U. Gross Exact factorization of the time-dependent electron-nuclear wave function Phys. Rev. Lett. 105 (2010) 123002
[4] A. Abedi, F. Agostini, Y. Suzuki, E. K. U. Gross Dynamical steps that bridge piecewise adiabatic shapes in the exact time-dependent potential energy surface Phys. Rev. Lett. 110 (2013) 263001
[5] F. Agostini, B. F. E. Curchod, R. Vuilleumier, I. Tavernelli, E. K. U. Gross, TDDFT and Quantum-Classical Dynamics: A Universal Tool Describing the Dynamics of Matter Springer International Publishing (2018) 1
[7] G. H. Gossel, F. Agostini, N. T. Maitra Coupled-trajectory mixed quantum-classical algorithm: A deconstruction J. Chem. Theory Comput. 14 (2018) 4513
[8] S. K. Min, F. Agostini, I. Tavernelli, E. K. U. Gross Ab initio nonadiabatic dynamics with coupled trajectories: A rigorous approach to quantum (de)coherence J. Phys. Chem. Lett. 8 (2017) 3048
[9] J. C. Tully Molecular dynamics with electronic transitions J. Chem. Phys. 93 (1990) 1061
[10] B. F. E. Curchod, F. Agostini On the dynamics through a conical intersection J. Phys. Chem. Lett. 8 (2017) 831




Module SCDM_WFs implements the selected columns of the density matrix (SCDM) method [1] for building localized Wannier Functions (WFs). Wannier90 [2] is a post-processing tool for the computation of the Maximally Localised Wannier Functions (MLWFs) [3,4,5], which have been increasingly adopted by the electronic structure community for different purposes. The reasons are manifold: MLWFs provide an insightful chemical analysis of the nature of bonding, and its evolution during, say, a chemical reaction. They play for solids a role similar to localized orbitals in molecular systems. In the condensed matter community, they are used in the construction of model Hamiltonians for, e.g., correlated-electron and magnetic systems. Also, they are pivotal 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, and many other applications, e.g. dielectric response and polarization in materials, ballistic transport, analysis of phonons, photonic crystals, cold atom lattices, and the local dielectric responses of insulators, for reference see [3]. This module is a first step towards the automation of MLWFs. In the original Wannier90 framework, automation of MLWFs is hindered by the difficult step of choosing a set of initial localized functions with the correct symmetries and centers to use as an initial guess for the optimization. As a result, high throughput calculations (HTC) and big data analysis with MLWFs have proved to be problematic to implement.

This module is part of the newly developed Wannier90 utilities within the pilot project on Electronic Structure Functionalities for Multi-Thread Workflows. The module is part of the pw2wannier interface between the popular QUANTUM ESPRESSO code link and Wannier90. It will be part of the next version of QUANTUM ESPRESSO v.6.3 and Wannier90. Moreover, it has been successfully added in a developer branch of the AiiDA workflow [6] to perform HTC on large material datasets.

Practical application and exploitation of the code

The SCDM-k method [1] removes the need for an initial guess altogether by using information contained in the single-particle density matrix. In fact, the columns of the density matrix are localized in real space and can be used as a vocabulary to build the localized WFs. The SCDM-k method can be used in isolation to generate well localized WFs. More interestingly is the possibility of coupling the SCDM-k method to Wannier90. The core idea is to use WFs generated by the SCDM-k method as an initial guess in the optimization procedure within Wannier90. This module is a big step towards the automation of WFs and simplification of the use of the Wannier90 program. The module is therefore intended for all the scientists that benefit from the use of WFs in their research. Furthermore, by making the code more accessible and easier to use, this module will certainly increase the popularity of the Wannier90 code.

[1] A. Damle, L. Lin, L. Ying SCDM-k: Localized orbitals for solids via selected columns of the density matrix J.Comp.Phys. 334 (2017) 1
[2] A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, N. Marzari wannier90: A tool for obtaining maximally-localised Wannier functions Com. Phys. Comm. 178 (2008) 685
[3] 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
[4] N. Marzari, D. Vanderbilt Maximally localized generalized Wannier functions for composite energy bands Phys. Rev. B 56 (1997) 12847
[5] I. Souza, N. Marzari, D. Vanderbilt Maximally localized Wannier functions for entangled energy bands Phys. Rev. B 65 (2001) 035109
[6] G. Pizzi, A. Cepellotti, R. Sabatini, N. Marzari, B. Kozinsky AiiDA: automated interactive infrastructure and database for computational science Comp. Mat. Sci. 111 (2016) 218


QQ-Interface (Quantics-QChem-Interface)

The QQ-Interface module connects the full quantum nonadiabatic wavefunction propagation code Quantics to the time-dependent density functional theory (TDDFT) module of the electronic structure program Q-Chem. Q-Chem provides analytic gradients, Hessians and derivative couplings at TDDFT level. With this module, it is possible to use the Q-Chem TDDFT module for excited state direct dynamics calculations. Quantics will start Q-Chem calculations whenever needed, prepare the input file from a template and will read the Q-Chem output file. The Q-Chem results are stored in the Quantics database and can be used in dynamics simulations. Due to the modular design of Quantics the TDDFT module of Q-Chem can be used for all dynamics simulations, e.g. direct dynamics variational multi-configurational Gaussian (dd-vMCG) or surface hopping simulations.

This module is part of a set of new functionalities developed for the Quantics program package during the E-CAM Extended Software Development Worksop: Quantum MD held at the University College Dublin.

Practical application and exploitation of the code

The module will be used to examine the nonadiabatic excited state dynamics of small to medium-sized molecules. The TDDFT module of Q-Chem allows treating systems that are too large for efficient multireference, such as CASSCF calculations. Until now photoinduced dynamics simulations of such molecules were only possible using trajectory-based algorithms. With Quantics a full quantum-mechanical description of the nuclear motion is possible.


PLUMED wrapper for OpenPathSampling


PLUMED is a widely used and versatile rare-event sampling and analysis code that can be used with various Molecular Dynamics (MD) engines. It has a very intuitive and versatile syntax for the definition of Collective Variables (CVs), and a wide variety of sampling methods, which accounts for its widespread use. The present module allows PLUMED and OPS to be used together. More details on the module can be found here.

Practical application and exploitation of the code

Transition path sampling simulations and analysis rely on accurate state definitions. Such states are typically defined as volumes in a Collective Variables space. OPS already supports a number of CVs, including the ones defined in the MDTraj python library. PLUMED offers a wide variety of extra CVs, which are enabled in OPS by this module. Many of PLUMED’s dozens of CVs have a biomolecular focus, but they are also general enough for other applications. PLUMED’s popularity (over 500 citations in 4 years after the release of PLUMED2 [1]) is greatly based on the fact that it works with many MD codes. OPS is now added to that list. The PLUMED code is well-maintained and documented for both users and developers. Several tutorials and a mailing list are available to address FAQs. More information about PLUMED is available here.


[1] G. Tribello, M. Bonomi, D. Branduardi, C. Camilloni, G. Bussi PLUMED 2: New feathers for an old bird Comput. Phys. Commun. 185 (2014) 604




Improving I/O of DL_MESO_DPD files using SIONlib


This module implements the SIONlib library to optimize the I/O (writing/reading) of the trajectory files generated by DL_MESO_DPD, the Dissipative Particle Dynamics (DPD) code from the DL_MESO package. SIONlib is a library for writing and reading binary data to/from several thousands of processors into one or a small number of physical files. For parallel access to files, only the open and close functions are collective, while the writing and reading of files can be done asynchronously. [1] In DL_MESO_DPD’s last release (version 2.6), the MPI version of DL_MESO_DPD generates multiple trajectory files, one for each MPI task. The interface with SIONlib optimizes the data writing so that just one physical file is produced from several MPI tasks. This drastic reduction in the number of output files is a benefit for the I/O of the code, and simplifies the maintenance of the output, especially for a large number of MPI tasks.

This module is part of the newly developed utilities for the DL_MESO_DPD code within the pilot project on Polarizable Mesoscale Models.

Practical application and exploitation of the code

The implementation of this module generates a single trajectory file (history.sion) in a parallel run of DL_MESO_DPD, instead of multiple (HISTORY) ones. Accordingly, analogous modifications have to be implemented in the post-processing utilities that read the HISTORY files. As an example, the changes were implemented in a formatting utility. Besides showing how to adapt the reading, this allows a robust check of the implementation, since the output is human readable, contains the full trajectories, and can be readily compared with outputs obtained using the standard version of DL_MESO_DPD.

The next released version of DL_MESO_DPD (in development) will tackle the writing of files differently, producing a single trajectory file from the start. However, the interface proposed here provides this feature to the users of version 2.6, and represents an alternative solution for the handling of the trajectories.

It should be noted that this implementation is meant to show the feasibility of the interfacing, not to deal with all the possible cases. Thus, the module’s functionality is restricted to the relevant case in which: i) the simulation is run in parallel using MPI, ii) a single SIONlib-type physical file is produced, and iii) the post-processing is done by a single process.

While SIONlib is optimized for a large number of MPI tasks, even the reduction from several output files to just one represents a benefit, for example when it comes to the maintenance of the simulation output.


[1] http://www.fz-juelich.de/ias/jsc/EN/Expertise/Support/Software/SIONlib/_node.html


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


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


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


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