Ionic and Molecular Liquids: Long-Range Forces and Rigid-Body Constraints

1. Charusita Chakravarty Indian Institute of Technology Delhi Ionic and Molecular Liquids:Long-Range Forces and Rigid-Body Constraints

2. Electrostatic Interactions • Ionic liquids: NaCl, SiO2, NH4Cl • mobile charge carriers which are atomic or molecular entities • Simple ionic melts: model ions as point charges with Coulombic interactions. Short-range repulsions control ionic radii. • Molecular Liquids • Electronic charge distributions show significant deviations from spherical symmetry • Can be modeled by: (a) multipole moment expansions or (b) arrays of partial charges water benzene

3. Multipole Expansion for Electrostatic Potential Represent the electronic charge distribution of a molecule by a set of multipoles:

4. Range of Electrostatic Interactions Long-range interactions: Tail correction will diverge for 1/rn interactions with n greater than or equal to 3; therefore minimum image convention cannot be applied

5. O 0.9572 Å 109.47 2d- (-1.04e) H H 2d- (-0.8472e) 1.0 Å d+ (+0.52e) d+ (+0.52e) O O H H 109.47 d+ (+0.4238e) d+ (+0.4238e) Partial Charge Distributions of Some Typical Molecules Dipole + Quadrupole TIP4P Water SPC/E Water Quadrupole Molecular Nitrogen

6. Array of Point Charges Coulombic contribution to the potential energy for an array of N charges that form a charge-neutral system: Electrostatic potential • Particle i interacts with all other charges and their mirror images but not with itself • Gaussian units • Cannot apply minimum image convention because sum converges very slowly

7. Poisson Equation Potential Energy and Forces Electrostatic Potential Charge Distribution Linear differential equation:

8. Ewald Summation for Point Charges Screened charge distribution: Converges fast in real space Co Point Charge Distribution: Converges slowly Gaussian compensating charge distribution: can be analytically evaluated in real space

9. Ewald Summation Screening a point charge to convert the long-range Coulombic interaction into a short-range interaction Evaluating the real-space contribution due to the screened charges The Poisson equation in reciprocal space for the compensating screened charge distribution Evaluating the reciprocal space contribution Self-interaction correction

10. Screening a point charge

11. Electrostatic potential due to a Gaussian charge distribution

13. Real Space Contribution The value of a must be chosen so that the range of interaction of the screened charges is small enough that a real space cutoff of rc < L/2 can be used and the minimum image convention can be applied

14. Poisson Equation in Reciprocal Space Fourier series representation of a function in a cubic box with edge length L and volume V under periodic boundary conditions: Fourier coefficients Poisson equation in real space Poisson equation in reciprocal space Reminder: Fourier transforms of derivatives 

15. Unit Point Charge at origin: Array of point charges Array of Gaussian charges Convolution of Point Charge distribution and smearing function FT of Point charge Array Green’s Function Fourier transform of smearing Function

16. Fourier Part of Ewald Sum Array of N Gaussian point charges with periodic images: Corresponding Electrostatic potential in reciprocal space Electrostatic potential in real space can be obtained using:

17. Fourier Part of Ewald Sum (contd.) Electrostatic potential in real space Reciprocal Space contribution to potential energy System is embedded in a medium with an inifinite dielectric constant

18. Correction for Self-Interaction Must remove potential energy contribution due to a continuous Gaussian cloudof charge qi and a point charge qi located at the centre of the cloud. Electrostatic potential due to Gaussian centred at origin

19. Coulombic Interaction expressed as an Ewald Sum Reciprocal space Self-Interaction Real Space Important: For molecules, the self-interaction correction must be modified because partial charges on the same molecule will not interact

20. Accuracy of Ewald Summation • Convergence parameters: • Width of Gaussian in real space, a • Real space cut-off, rc • Cutoff in Fourier or reciprocal space, kc=2p/Lnc

21. Calculating Ewald Sums for NaCl Na+ Cl- Crystal Liquid

22. Hands-on Exercise:Calculating the Madelung Constant for NaCl The electrostatic energy of a structure of 2N ions of charge +/- q is where a is the Madelung constant and rnn is the distance between the nearest neighbours.

23. Symmetric Stretch 3657cm-1 Bend 1595cm-1 Asymmetric Strech 3756cm-1 Rotational and Vibrational Modes of Water http://chsfpc5.chem.ncsu.edu/~franzen/CH795N/lecture/XIV/XIV.html http://www1.lsbu.ac.uk/water/vibrat.htm

24. Molecules: Multiple Time-scales • Bonded interactions are much stronger than non-bonded interactions • Intramolecular vibrations have frequencies that are typically an order magnitude greater than those of intermolecular vibrations • MD/MC: time step will be dictated by fastest vibrational mode • Fast, intramolecular vibrational modes do not explore much of configuration space- rapid, essentially harmonic, small amplitude motion about equilibrium geometry • Require efficient sampling of orientational and intermolecular vibrational motions

25. Simulation Methods for Molecules • Freeze out all or some intramolecular modes: • Serve to define vibrationally averaged molecular structure and are completely decoupled from intermolecular vibrations, librations or rotations • Rigid-Body Rotations: • Characterize the mass distribution of the molecule by its moment of inertia tensor • MC: Use orientational moves • MD: Propagate rigid-body equations of motion • Will not work if there are low-frequency vibrational modes • Apply Geometrical Constraints • MD: SHAKE • MD: RATTLE • Multiple Time-Scale Algorithms

26. N2 H2 SF6 CH4 NH3 lCl5 H2O Ix Iy  Iz Rotations of Rigid-Bodies Space-fixed (SF) and Body-fixed (BF) axes (Goldstein, Classical Mechanics) Moments of Inertia of Molecules: Ia < Ib < Ic Linear: Spherical polar angles: (q,f) Non-linear: Euler angles: (f,q,y) (Atkins, Physical Chemistry)

27. Monte Carlo Orientational Moves for Linear Molecules Orientation of a linear molecule is specified by a unit vector u. To change it by a small amount: • Generate a unit vector v with a random orientation. See algorithm to generate random vector on unit sphere • Multiply vector v with a scale factor g, which determines acceptance probability of trial orientational move • Create a sum vector: t = u + gv • Normalise t to obtain trial orientation

28. Euler Angles (f,q,y) Euler’s rotation theorem: Any rotation of a rigid-body may be described by a set of three angles • Rotation, A: Initial orientation of body-fixed axes (X,Y,Z) to final orientation (X’’,Y’’’, Z’) http://stackoverflow.com/questions http://mathworld.wolfram.com/EulerAngles.html

29. Euler angles (f,q,y)(contd.) • Rotate the X-, Y-, and Z-axes about the Z-axis by f (p < f < p ), resulting in the X'-, Y'-, and Z-axes. • Rotate the X'-, Y'-, and Z-axes about the X'-axis by q (0<q<p), resulting in the X'-, Y''-, and Z'-axes. • Rotate the X'-, Y''-, and Z'-axes about the Z'-axis by y ( -p <y<p ), resulting in the X''-, Y'''-, and Z'-axes. Rotation A = BCD, therefore new coordinates are:

30. Monte Carlo Orientational Moves for Non-linear Molecules • Specify the orientation of a rigid body by a quaternion of unit norm Q that may be thought of as a unit vector in four-dimensional space.

31. Applying Geometrical Constraints Lagrangian Equations of Motion Kinetic Energy Potential Energy Cartesian coordinates Geometrical Constraints Define constraint equations and require that system moves tangential to the constraint plane.

32. Introduce a new Lagrangian that contains all the constraints: The equations of motion of the constrained system are: Constraint force acting along coordinate qi

33. Bond Length Constraint for Diatomic Molecule m1 m2 Bond constraint • Constraint forces: • lie along bond direction • are equal in magnitude • opposite in direction • do no work Verlet algorithm: New position with constraint Constraint forces on atom i Unconstrained position

34. Diatomic molecule (contd) Three unknowns In the case of a diatomic molecule, can obtain quadratic equation in l. However, cannot be conveniently generalised for larger molecules with more constraints.

35. Multiple Constraints Verlet algorithm for N atoms in the presence of l constraints: Taylor expansion of the constraints with respect to unconstrained positions. Impose the requirement that after one time step, the constraints must be satisfied. Constraint force acting on i-th atom due to k-th constraint Substitute from Taylor expansion above

36. Multiple Constraints (contd.) Need to find the unknown Lagrange multipliers: Solution by Matrix Inversion Since the Taylor expansion was truncated at first-order iterative scheme will be required to obtain convergence.

37. SHAKE Algorithm • Iterative scheme is applied to each constraint in turn: RATTLE: Similar approach within a velocity Verlet scheme

38. References • D. Frenkel and B. Smit, Understanding Molecular Simulations: From Algorithms to Applications • A. R. Leach, Molecular Modelling: Principles and Applications • D. C. Rapaport, The Art of Molecular Dynamics Simulation (Details of how to implement algorithms for molecular systems) • M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (SHAKE, RATTLE, Ewald subroutines)