Optimisation of Classical Potential parameters for Charge Equilibration Method from Ab Initio Data

1. Optimisation of Classical Potential parameters for Charge Equilibration Method from Ab Initio Data Emeric Bourasseau Jean-Bernard Maillet Vincent Dubois CEA/DIF Département de Physique Théorique et Appliquée Service de Physique de la Matière Condensée Bruyères-le-Châtel

2. General purpose Development of an optimisation process to obtain transferable classical potentials to simulate complex fluids under high pressure conditions using classical molecular dynamic High pressure conditions: Few experimental data Ab Initio data Transferable parameters: Physical meaning Rational optimisation process

3. Usual optimisation process: Minimisation of an error function that evaluates the accuracy obtained with a set of parameters with respect to a predefined data basis. Three key points: The choice of the potential form Separation of the different physical contributions The building of the data basis Use of reference data that only depend on separated contributions The definition of the error function and the minimisation process Sequential minimisation of different error functions The choice of the potential form The building of the data basis The definition of the error function and the minimisation process

4. Example: Parameter Optimisation of electrostatic potential of H, C, N, O and F atoms Reference data from HF and CH3NO2 systems Key point 1: The potential form: Partial charges must depend on their environment We determine qi and qj with the charge equilibration method Partial charges appear more realistic ≡ more transferable

5. Charge equilibration method: based on the principle of electronegativity equalization + neutrality of each molecule For each molecule i : • System of N linear equations with N unknowns {qia}. Atomic electronegativity self-coulomb integral Coulomb interactions : First parameter (atomic electronegativity)

6. Calculation of Jab (coulomb interaction) • J(r) e2/(40r) if r   Long distance interactions : Overlap at short distances : With Calculation of J°aa (self-coulomb integral) : Second parameter (Slater exponent)

7. Key point 2: The reference data basis: Reference data used to optimise electrostatic parameters should only depend on electrostatic potential contribution Only data calculated on fixed molecular configurations could be used • Suitable data: • Dipole moment of isolated molecule • Atomic charges from an ab initio analysis • Electrostatic potential from electronic density The last two types of data are calculated on several configurations taken from ab initio trajectories of liquid systems.

8. Key point 3: The error function: a is the tolerated error is used to weight the different F functions to obtain comparable values

9. Results:First try: Optimisation of parameters c0 and z of atoms H and F on the basis of HF data 3 parameters (c0H=0 J): c0F, zH and zF. We applied physical considerations to reduce the parameter space. Reference data: Dipole moment of the isolated molecule Electrostatic potential grids taken from 4 configurations of HF in the liquid phase (4 densities between 1.3 and 2.6 g.cm-3 at 500 K)

10. Error function: We performed 50 minimisations from 50 random initial points: we obtained 38 satisfying minima

11. Calculation of the F function in a 3-D space Slice in a plan [zH;zF] for c0F = 4.75 J

12. Best minimum: c0H = 0 J zH = 1.916.1010 m-1 c0F = 5.182 J zF = 4.277.1010 m-1 (Mol. Sim. 31, pp. 705-713, 2005) mcalc = 1.792 D mref = 1.826 D Ab initio electrostatic potential Classical electrostatic potential

13. Results:Second try: Optimisation of parameters c0 and z of atoms H, C, N and O on the basis of CH3NO2 data 7 parameters (c0H=0 J): c0C, c0N,c0O, zH, zC, zN and zO. We applied physical considerations to reduce the parameter space. Reference data: Dipole moment of the isolated molecule Partial charges calculated applying the Bader analysis on 3 configurations of CH3NO2 in the liquid phase (3 densities between 0.9 and 1.6 g.cm-3 at 300 K) Only charges of atoms C, N and O have been used because charges on hydrogen atoms are too close to 0.

14. Error function: We performed 500 minimisations from 500 random initial points: we obtained 12 satisfying minima The complete investigation of the error function in the 7-D space is too long : 1019 evaluations would be needed. We used physical considerations to determine the best minimum. c0H = 0 J zH = 1.877.1010 m-1 c0C = -0.765 J zC = 2.313.1010 m-1 c0N = 0.429 J zN = 3.019.1010 m-1 c0O = 4.256 J zO = 3.194.1010 m-1 Physical tendencies: c0C<c0H<c0N<c0O and zC<zN<zO

15. Those parameters show interesting realism: c0H = 0 J zH = 1.877.1010 m-1 c0C = -0.765 J zC = 2.313.1010 m-1 c0N = 0.429 J zN = 3.019.1010 m-1 c0O = 4.256 J zO = 3.194.1010 m-1 c0H = 0 J zH = 1.877.1010 m-1 c0H = 0 J zH = 1.916.1010 m-1 c0C = -0.765 J zC = 2.313.1010 m-1 c0N = 0.429 J zN = 3.019.1010 m-1 c0O = 4.256 J zO = 3.194.1010 m-1 c0F = 5.182 J zF = 4.277.1010 m-1 Physical tendencies: c0C<c0H<c0N<c0O<c0F and zC<zN<zO<zF mcalc = 3.72 D mref = 3.56 D RRMS : Configuration 1 9.57 % (7.8 % on CNO) Configuration 2 10.3 % (8.5 % on CNO) Configuration 3 7.9 % (6.7 % on CNO)

16. -0.003 -0.028 -0.412 -0.398 H H O O 0.423 0.427 C N C N H H 0.358 0.362 O O H H 0.0 -0.04 -0.404 -0.418 0.05 0.08 0.197 -0.34 H O -0.29 C N H 0.375 O H 0.188 -0.353 0.267 Comparison of partial charges obtained on a molecule taken from a liquid configuration of nitromethane at 300 K and 1.125 g.cm-3 Ab Initio (Bader Analysis) Our model Reax FF(1) qCH3 = 0.354 qCH3 = 0.439 qCH3 = 0.471 (1) A. C. van Duin, J. Phys. Chem. A 105, pp. 9396-9409 (2001)

17. Results:Third try: Optimisation of parameters c0 and z of atoms H, C, N, O and F on the basis of HF and CH3NO2 data 9 parameters (c0H=0 J): c0C, c0N,c0O, c0F, zH, zC, zN, zO and zF. We applied physical considerations to reduce the parameter space. Reference data: Dipole moment of the isolated HF molecule Electrostatic potential grids taken from 4 configurations of HF in the liquid phase (4 densities between 1.3 and 2.6 g.cm-3 at 500 K) Dipole moment of the isolated CH3NO2 molecule Partial charges calculated applying the Bader analysis on 3 configurations of CH3NO2 in the liquid phase (3 densities between 0.9 and 1.6 g.cm-3 at 300 K)

18. Error function: We performed 150 minimisations from 150 random initial points: we obtained 20 satisfying minima The complete investigation of the error function in the 9-D space is impossible: 1023 evaluations would be needed. We used physical considerations to determine the best minimum.

19. Those parameters show interesting realism: c0H = 0 J zH = 1.888.1010 m-1 c0C = -0.769 J zC = 2.312.1010 m-1 c0N = 0.324 J zN = 3.060.1010 m-1 c0O = 4.433 J zO = 3.314.1010 m-1 c0F = 5.879 J zF = 4.706.1010 m-1 mcalc(CH3NO2) = 3.73 D mref(CH3NO2) = 3.56 D mcalc(HF) = 1.783 D mref(HF) = 1.826 D RRMS (CH3NO2) : Configuration 1 9.5 % (7.7 % on CNO) Configuration 2 10.4 % (8.5 % on CNO) Configuration 3 8.0 % (6.7 % on CNO)

20. Conclusions We have developed a global method to obtain transferable parameters used in charge equilibration method. The optimised parameters reproduce accurately ab initio electrostatic properties and ab initio partial charges in the liquid phase. The process is not completely automatic: man hand is always needed to determine the best minimum Perspectives: To extend the reference data basis To perform the same work for other potential contributions. This is in progress (see Vincent Dubois' poster)