Romain Dupuis
77 Massachusetts Ave, E19-722, Cambridge, MA 02139

Archives de l’auteur

Séminairé à l’INSP « Autour des effets nucléaires quantiques » (about quantum nuclear effects)

19/01/2016 – Ce séminaire à été donné à l’INSP (Institut des NanoSciences de Paris) afin de présenter l’importance de prendre en compte les effets nucléaires quantiques lorsque l’on s’intéresse à des systèmes composés d’atomes légers (H, Li) ou sous contraintes. La méthode des intégrales des chemins a été présentée puis ses points forts ont été mis en évidence au travers de deux exemples : le fractionnement isotopique du Li entre un minéral et une solution ; le déplacement des atomes d’hydrogène dans la Portlandite (hydroxyde) sous très haute pression.

Fractionation of silicon isotopes in liquids: the importance of configurational disorder

Fractionation of silicon isotopes in liquids: the importance of configurational disorder

Dupuis, Romain ; Benoit, M. ; Nardin, E. ; Méheut, M.



Abstract Silicon isotopes are a promising tool to assess low-temperature geochemical processes such as weathering or chert precipitation. However, their use is hampered by an insufficient understanding of the fractionation associated with elementary processes such as precipitation or dissolution. In particular, the respective contributions of kinetic and equilibrium processes remain to be determined. In this work, equilibrium fractionation factors for silicon isotopes have been calculated using first-principles methods for quartz, kaolinite, and dissolved silicic acid (H4SiO4 and H3SiO4−) at 300 K. The two liquid systems are treated both as realistically as possible, and as consistently with the solids as possible. They are first simulated by ab initio molecular dynamics, then individual snapshots are extracted from the trajectories and relaxed, giving inherent structures (IS) and their fractionation properties are calculated. The fractionation properties of these IS are then calculated. A significant variability of the fractionation properties (σ= 0.4‰) is observed between the independent snapshots, emphasizing the importance of configurational disorder on the fractionation properties of solutions. Furthermore, a correlation is observed between the fractionation properties of these snapshots and the mean Si-O distances, consistent with calculations on minerals. This correlation is used to identify other parameters influencing the fractionation, such as the solvation layer. It is also used to reduce the number of configurations to be computed, and therefore the computational cost. At 300 K, we find a fractionation factor of + 2.1±0.2‰ between quartz and H4SiO4, + 0.4±0.2‰ between kaolinite and H4SiO4, and -1.6±0.3‰ between H3SiO4− and H4SiO4. These calculated solid-solution fractionations show important disagreement with natural observations in low-temperature systems, arguing against isotopic equilibration during silicon precipitation in these environments. On the other hand, the large fractionation associated with the de-protonation of silicic acid suggests the importance of speciation, and in particular pH, for the fractionation of silicon isotopes at low temperature.

Postdoc at DIPC

Invited PhD at the DIPC (San Sebastian, Spain) from september 2014 to december 2014.

Postdoctoral contract at DIPC starting from 10 Dec. 2014 on the study of cementitious materials with ab initio calculations.

Seminar at DIPC

Seminar titled : Realistic calculations of the isotopic fractionation factors of Si and Li for equilibriums involving liquid phases on the 17th of July at the DIPC of Donostia.

Homepage of the DIPC :

J. Chem. Theory Comput., 2014, 10 (4), pp 1440–1453

Efficient calculation of free energy differences associated with isotopic substitution using Path Integral Molecular Dynamics

Marsalek, Ondrej ; Chen, Pei-Yang ; Dupuis, Romain ; Benoit, Magali ; Méheut, Merlin ; Bačić, Zlatko ; Tuckerman, Mark E.


The problem of computing free energy differences due to isotopic substitution in chemical systems is discussed. The shift in the equilibrium properties of a system upon isotopic substitution is a purely quantum mechanical effect that can be quantified using the Feynman path integral approach. In this paper, we explore two developments that lead to a highly efficient path integral scheme. First, we employ a mass switching function inspired by the work of Ceriotti and Markland [ J. Chem. Phys. 2013138, 014112] that is based on the inverse square root of the mass and which leads to a perfectly constant free energy derivative with respect to the switching parameter in the harmonic limit. We show that even for anharmonic systems, this scheme allows a single-point thermodynamic integration approach to be used in the construction of free energy differences. In order to improve the efficiency of the calculations even further, however, we derive a set of free energy derivative estimators based on the fourth-order scheme of Takahashi and Imada [ J. Phys. Soc. Jpn. 198453, 3765]. The Takahashi–Imada procedure generates a primitive fourth-order estimator that allows the number of imaginary time slices in the path-integral approach to be reduced substantially. However, as with all primitive estimators, its convergence is plagued by numerical noise. In order to alleviate this problem, we derive a fourth-order virial estimator based on a transferring of the difference between second- and fourth-order primitive estimators, which remains relatively constant as a function of the number of configuration samples, to the second-order virial estimator. We show that this new estimator converges as smoothly as the second-order virial estimator but requires significantly fewer imaginary time points.

Surface Science 618 (2013) 72–77

Development of Monte-Carlo simulations for nano-patterning surfaces associated with MM-EPES analysis: Application to different Si(111) nanoporous surfaces


EPES (elastic peak electron spectroscopy) allows measuring the percentage of elastic backscattered electrons ηe from a surface excited by primary electrons. However, this method must be combined with Monte Carlo simulations to get quantitative information. After a brief description of the algorithm used in this work (named MC2), we focused on the adaptation of this simulation for nanoporous surfaces (named MC2-NP). The theoretical results obtained put in evidence the dependence of ηe versus pore diameter (d), depth of the pores (h) and covering rate (CR) of the pores on the surface. Results obtained on surfaces having cylinder-shaped and cone-shaped holes with nanometer dimensions are presented too. To validate theoretical results obtained with MC2-NP, silicon(111) nanoporous surfaces have been prepared with an anodized aluminum oxide (AAO) template and by argon ion bombardment in an UHV chamber. Uniform nanohole arrays were formed as a replica of ordered lattice pattern of the template. Then EPES experimental measurements have been performed on planar and nanoporous Si(111) surfaces using a retarding field analyzer (RFA). The experimental results put in evidence that the percentage of the elastically backscattered electrons is influenced by the patterning of the surface. Then comparing values of ηe obtained experimentally with those obtained with MC2-NP simulations, we show the sensitivity of the EPES method for studying nanoporous surfaces. In this way, we expect fast estimation of nanohole’s dimensions by in-situ MM-EPES (Multi-Mode EPES) without other techniques such as for example scanning electron microscopy.


  • EPES;
  • MM-EPES;
  • Nanohole;
  • Nanoporous alumina mask;
  • Surface plasmon;
  • Monte Carlo simulation

Research Projects


The fractionation of isotopes is widely used in geology to interpret global processes such as CO2 cycle and weathering[1]. The calculation of the isotopic fractionation factor at the equilibrium in minerals is now faithfully performed ([2], [3]) and gives insights on the mechanisms of fractionation.

However, in natural systems, most of equilibriums involve a liquid phase. But the fast growth of measurement methods (in complicated systems -biotic, liquid) has not been catched up by calculations. This is challenging because, for consistency, solutions and liquids have to be treated with the same method. The harmonic approx- imation that is commonly accepted to treat minerals was straightforwardly extended to liquids. Results on iron isotopes for the equilibrium of two species in a solution are promising[4] but the values of the fractionation at the equilibrium between minerals and solution are unfortunately in disagreement with the experiment. It suggests that the error (due to the anharmonicity in liquids) cancels out for similar phases but can not be neglected for mineral/solution equilibriums.

PhD project

The objective of my PhD project was to calculate realistic and consistent fractionation factors, which requires the use of advanced ab initio methods performed on HPC centers. Two approaches were chosen, depending on the kind of equilibrium. Concerning species in liquids, following the approach proposed in [4], we investigated within the harmonic approximation the effect of the configurational disorder in liquids on the fractionation factor estimation . Then, we studied solution/mineral equilibriums with a consistent method that takes into account the effect of anharmonicity (Path Integral Molecular Dynamics coupled with Thermodynamics Integration). Efforts have been done to reduce the important computational resources required by this method.
Configurational disorder in solutions The fractionation factor can be computed as an average over indepen- dent configurations extracted from an ergodic trajectory. In former work on the fractionation factor, only one configuration was considered. In this study, we investigated the viability of this approximation and proposed a more efficient sampling that is based on the study of the liquids structural properties and on its impact of the fractionation properties.

Trajectories of H4SiO4, H3SiO− 4 or H2SiO2− 4 in water were simulated within the framework of Car-Parrinello MD (with the CPMD package with Troullier-Martins pseudopotentials and the BLYP functional). Analysis codes were written to study the structural liquid properties.

The method that was developed is based on the work of Rustad and Bylaska[5]: Several configurations are relaxed (with the pw-scf code Quantum Espresso); The vibrational spectra is computed using Density Func- tional Perturbation Theory (for liquids the Γ-point approximation has been tested); Finally the fractionation factor is calculated at the temperature of the dynamics as detailed in Méheut et al.[3].

A fine study of the vibrational properties of different configurations has been done in order to point out the modes that contribute the most to the fractionation properties. We also investigated the correlation between the structural properties of liquids and the calculated value of the fractionation factor with analytical in-house codes (i.e. statistics on the number of water molecules in the first solvation layer, bond-length distribution or size of the water rings around the silicon atom).

Apart from liquids, we also computed the Quartz and Kaolinite minerals within the harmonic approxima- tion with the same pseudopotentials as for the dynamics. Quartz is used as a reference to compare the results to experiments and to other calculations[3] performed with the PBE functional.

Two publications on this study will be submitted for publication soon: One regarding the employed method- ology and the second one on the geological implications of the obtained results.

Anharmonicity effects in solution

Since the fractionation factor between a solution and a mineral is impor- tant in geology (i.e. for precipitations, dissolutions and adsorptions), our objective was to consider a consistent method for both phases which takes into account quantum anharmonic effects. A method of choice is the Path Integral Molecular Dynamics coupled with Thermodynamics Integration[6]. The fractionation factor between Li+ in water and Li2O was studied. For comparison with experiments[7], clays have also been computed within the harmonic approximation.

For this study, we closely collaborated with the group of research of Mark. E. Tuckerman, which is a spe- cialist of PIMD, in order to develop this method and optimize our calculations on geological systems. This collaboration has resulted in the publication of the developed methodology in J. Chem. Theory Comput. [8].

The test of this method on ”real” systems has first been done with empirical potentials in order to have an estimation of the best parameters for the path-integrals dynamics and to observe the differences between anharmonic (performed with PINY package) and harmonic calculations (with Gulp and PINY). At the same time, we studied the reliability of the potentials. We compared the vibrational properties, which drive the fractionation, with experimental data[9] and with ab initio calculations. We observed an important effect due to anharmonicity in solution whereas it is negligible in Li2O, and the obtained fractionation factors are in very good agreement with the experiments.

The ongoing part of this study consists in ab initio calculations on similar systems (with anharmonic methods) and on more complex Li-containing minerals (with harmonic methods) in order to compare the results more precisely with experiments. The objective is to finally have a consistent and realistic calculation of the fraction- ation factor between a mineral and a solution which would be a world first.

This work will result in the publication of two additional articles which are currently in preparation.


During this project, I worked with specialists of several methods for applications to realistic systems in the field of geology. Some of these methods are widely used in large fields such as physics or biology (namely DFT, MD and CPMD) and PIMD is an attractive method that allows the calculation of very fine quantum effects. The main objective of my PhD thesis, which was to develop a method to compute the fractionation factor in complicated systems involving liquid phases, was successfully achieved. The results were presented at the Goldschmidt conference in 2013 and grew interest for geologists, for instance in order to understand the fractionation properties of silicon at extreme pH. More developments are now foreseen to improve the methods in order to have even more realistic results.


[1] Opfergelt S. and Delmelle P. (2012) silicon isotopes and continental weathering processes: Assessing controls on Si transfer to the ocean. C R Geosci. 344, 723–738.

[2] Schauble E. A. (2011) First-principles estimates of equilibrium magnesium isotope fractionation in silicate, oxide, carbonate and hexaaquamagnesium(2+) crystals. Geochim. Cosmochim. Acta 75, 844 – 869.

[3] M´eheut M., Lazzeri M., Balan E. and Mauri F. (2007) Equilibrium isotopic fractionation in the kaolinite, quartz, water system: Prediction from first-principles density-functional theory. Geochim. Cosmochim. Acta 71, 3170 – 3181.

[4] Beard B. L., Handler R. M., Scherer M. M., Wu L., Czaja A. D., Heimann A. and Johnson C. M. (2010) Iron isotope fractionation between aqueous ferrous iron and goethite. Earth Planet. Sci. Lett. 295, 241 – 250.

[5] Rustad J. R. and Bylaska E. J. (2007) Ab Initio Calculation of Isotopic Fractionation in B(OH)3(aq) and BOH4-(aq). J. Am. Chem. Soc. 129, 2222–2223.

[6] Ceriotti M. and Markland T. E. (2013) Efficient methods and practical guidelines for simulating isotope effects. J. Chem. Phys. 138, 014112

[7] Vigier N., Decarreau A., Millot R., Carignan J., Petit S. and France-Lanord C. (2008) Quantifying Li isotope fractionation during smectite formation and implications for the Li cycle. Geochim. Cosmochim. Acta 72, 780 – 792.

[8] Marsalek O., Chen P-Y, Dupuis R., Benoit M., M´eheut M., Baˇci´c Z. and Mark E. Tuckerman (2014) Efficient Calculation of Free Energy Differences Associated with Isotopic Substitution Using Path-Integral Molecular Dynamics. J. Chem. Theory Comput. 10 (4), 1440-1453

[9] Prabhatasree G., Choudhury N. and Chaplot S. L. (2004) Lattice dynamics of lithium oxide. Pramana journal of physics 63, 409-412