Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics

Article pubs.acs.org/jcim

Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator Alain Dequidt,*,† Julien Devémy,† and Agílio A. H. Pádua†,‡ †

Institut de Chimie de Clermont-Ferrand, Université Blaise Pascal, CNRS, 63178 Aubière, France Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, United States

‡

ABSTRACT: LAMMPS is a very customizable molecular dynamics simulation software, which can be used to simulate a large diversity of systems. We introduce a new package for simulation of polarizable systems with LAMMPS using thermalized Drude oscillators. The implemented functionalities are described and are illustrated by examples. The implementation was validated by comparing simulation results with published data and using a reference software. Computational performance is also analyzed.

■

temperature.16 This technique is faster but requires careful selection of force field parameters and simulation conditions to avoid spurious effects. Of the three strategies to include polarization listed above, induced point dipoles are the least used in molecular dynamics (MD) because of the computational cost of evaluating multipole electrostatics at long range; induced dipoles are not presently implemented in the LAMMPS code. This method is implemented for instance in the AMOEBA force field of AMBER20 but is not as fast as the Drude approach.21 Fluctuating charge methods are based on electronegativity equalization11 and rely on the point charges placed on existing sites, therefore they are computationally economic. Charge transfer is allowed between molecules or atom groups (although it can be constrained19). However, in planar molecules or atom groups, fluctuating point charges cannot represent off-plane components of polarization. The fluctuating charge model is implemented in LAMMPS through the QEQ package, both in the matrix-inversion (qeq/ point) and extended Lagrangian (qeq/dynamic) methods.22 Drude oscillators represent induced dipoles associated with interaction sites by pairs of charges: the core atom and a Drude particle, bound by a harmonic spring. Thus, Drude models are able to represent off-plane components of polarization and, since they rely on point charges only, are compatible with long-range electrostatics as implemented in most simulation codes. Drude oscillators are computationally more expensive than fluctuating charges because an additional site is present at each polarizable atom. “Drude oscillator” models are often confused with “core− shell” models,23 the former designation being more common in studies of molecular systems and liquids, whereas the latter is more usual in the scope of crystalline ionic materials. The core−shell

INTRODUCTION One of the main evolutions in classical force fields for molecular simulation is the inclusion of explicit polarization, i.e., electrostatic charge distributions which respond to their environments, introducing a higher degree of realism in the representation of molecules and materials. Force fields that include explicit polarization effects can be considered at the frontier between two of the conventional scales of description in molecular modeling: the quantum level in which electrons are treated explicitly, and the atomistic level where the particles are atomic. Interestingly, the inclusion of dispersion interactions in (quantum) density functional theory through semi-classical terms1,2 is also an active field at present, bridging the gap between the same two levels of description. The objective of the package described here is to enrich the capabilities of LAMMPS,3 a free and open-source molecular dynamics software, toward the simulation of molecular systems and materials using polarizable force fields. LAMMPS is an outstanding code in its versatility and ability to study a wide diversity of physical and chemical systems described by atomistic and coarse-grained interaction models. Several methods have been proposed to include explicit polarization in force fields,4 namely induced point dipoles (or multipoles),5−9 fluctuating charges,10−13 and Drude dipoles.14−16 This field has been reviewed recently.17,18 Force fields with explicit polarization are non-additive, requiring calculation of induced dipoles or fluctuating charges during the simulation trajectory. Since the response of the electrons is much faster than the atomic motions, the Born−Oppenheimer approximation holds and the instantaneous dipoles or charges should be computed in a self-consistent way. Solving for the self-consistent field at each time-step is, however, computationally slow. For all of the above-mentioned methods an approximate solution can be sought through an extended Lagrangian approach9,12 and by thermalizing the additional degrees of freedom at low © XXXX American Chemical Society

Received: October 7, 2015

A

DOI: 10.1021/acs.jcim.5b00612 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling

Ideally, the masses of the Drude particles should be small and the stiffness of the harmonic bonds should be large, ensuring that the Drude particles remain close to their cores. However, this implies very fast oscillations of the Drude particles and would require a very small time step, and it is why the relative coordinates of the DP with respect to its DC should be thermostated at low temperature,16 thus approximating the selfconsistent regime. The criteria for choosing the Drude oscillator parameters are detailed below, and a Python script is provided with the Drude package to apply these modifications to a fixedcharge force field (usage of this script is also documented). Topology. A Drude oscillator consists of a pair of “atoms” (in LAMMPS terminology) which are the Drude core (DC) and the Drude particle (DP). The core and particle bear electric charges and are connected to each other by a harmonic bond. So the atom style atom_style full, which allows for electrostatic and bonded interactions, is appropriate. A new fix drude was created to set up the topology (how the atoms are bonded), letting each atom store the reference of its Drude partner. In addition, an array of “Drude types” stores whether each atom type corresponds to a non-polarizable atom, a DC or a DP. Internally, the fix identifies Drude partners, by scanning the bonded list: if one of the partners is a DP, then the other is its associated DC. This method is called when a data file is read, when a new molecule is created or the fix drude is initialized. Most force fields apply special pair interactions between neighboring, bonded sites (which are connected by covalent bonds, valence angles or torsion potentials). The pair interactions (Coulombic and van der Waals) between 1−2 to 1−4 neighbors are often scaled down. Usually, 1−2 and 1−3 pair interactions are entirely removed, whereas 1−4 pair interactions are scaled by a certain factor according to the specification of the force field. The introduction of Drude dipoles requires a redefinition of these special relations, because a DC−DP pair stands for a single polarizable atom. So the same special factors must be applied to the DPs as to the corresponding DCs. Internally, LAMMPS builds a list of special neighbors storing those special relationships. The fix drude modifies this list of special neighbors. The DPs are first removed from the lists of special neighbors, then the lists are scanned and, for each DC found, its corresponding DP is appended just after, with the same special relationship (1−2 to 1−4) as itself. The partner DP is inserted at the beginning of the special list of each DC. DC and DP partners would ideally be 1−1 partners but this concept does not exist in LAMMPS. So they are defined as 1−2 partners. The 1−2 pair interactions should then be entirely removed. Finally, the special list of the DPs are copied from that of their DC. The modifications in the topology from the introduction of Drude oscillators have the following important consequences: • The 1−2 pair interactions must be totally scaled down by a special factor 0.0, otherwise the cores and respective Drude particles will interact. • If using fix shake to constrain bonds or angles, this fix must be invoked on a group of atoms not including the Drude particles. This is because fix shake makes use of the special neighbors list to assess connectivity, and the neighbor relations of the DPs are particular. Thole Damping Function. If two Drude oscillators are too close, for example by belonging adjacent atoms in the same molecule, DCs strongly attract the neighbors’ DPs leading to excessive correlation of induced dipoles, or even to instability of

model is implemented in MD codes suitable for materials, such as DL_POLY24 and LAMMPS. Both models are based on the same idea of representing induced dipoles by a pair of point charges connected by a spring. However, what is conventionally called the Drude oscillator model includes two important features: one is a technique to handle thermostating of the additional degrees of freedom, and the other a special treatment of the interactions between neighboring induced dipoles located on the same molecule (Thole damping). The Drude model is therefore more complex to implement (due to the use of specific thermostats, special neighbor relations and short-range damping functions) and to use (in terms of simulation setup) than the core−shell model. The CHARMM25 and NAMD26 molecular dynamics codes implement explicit polarization through Drude dipoles, and the CHARMM force field for biomolecules includes Drude polarization. However, those codes are tailored for organic molecules and biomacromolecules, and are not equipped with the variety of interaction potential models suitable to represent wide classes of molecular, covalent, ionic or metallic materials. Therefore, the present package fills in a missing functionality in LAMMPS, making it an even more complete code in terms of the interaction models that can be simulated and of the diversity of systems that can be studied. As explained above, correct use of Drude dipole polarization is not trivial and requires detailed knowledge of the main implementation principles and its design, as well as how to put all the elements together in LAMMPS for a successful usage leading to rigorous results. We think this intricacy justifies the present article. In the next section, the present implementation of Drude oscillators is documented. In the subsequent section usage examples of the package are illustrated. Then the implementation is validated by comparison with published results. Finally, the performance of the implementation is assessed.

■

IMPLEMENTATION OF DRUDE OSCILLATORS LAMMPS is very customizable and allows choices between several atom styles, bond styles, pair interaction styles, etc. The present implementation makes use of existing infrastructure in the LAMMPS code and its packages as much as possible. For example, Drude cores (DC) and Drude particles (DP) will be described by “normal” atoms, and the core−Drude bonds by standard harmonic bonds. Note that any style of LAMMPS bonds can be used instead of harmonic bonds. For instance, FENE bonds with finite extension could be used to limit the DC-DP distance. On the other hand, no anisotropic bond style exists in LAMMPS yet, so anisotropic polarizability27 cannot be modeled. Due to the special nature of the induced Drude dipoles, some new specific properties were included in the form of a new fix drude. In the specification of the force field and per-atom properties, namely charge and mass, the user has to pay attention to the characteristics of the Drude oscillators: how much mass to assign to the DP, which should be subtracted from that of the initial, non-polarizable atom in order to obtain the mass of the respective DC; also, what value of charge to place on each DP (depending on the polarizability of the atom) which should be subtracted from the charge of the initial, non-polarizable atom to yield the charge of the respective DC. Finally, the force constant of the DC−DP bonds as well as their zero equilibrium length need to be specified. These characteristics of atom and bond types (masses, force constants), or properties of individual atoms (charges) are specified in the LAMMPS data file. B


Article

Journal of Chemical Information and Modeling trajectory or collapse of DPs onto nearby DCs. The behavior of neighboring Drude dipoles within the same molecule backbone is more realistic and stable if a suitable screening or damping function is employed at short range. The Thole’s damping function Sij which was implemented dampens the interaction potential by28,29 qiqj ϕ(rij) = Sij(rij) rij (1) sijrij ⎞ ⎛ ⎟ exp( − s r ) Sij(rij) = 1 − ⎜1 + ij ij ⎝ 2 ⎠

(2)

The Thole scaling coefficient sij is determined by the polarizability of the atoms, αi, and by a Thole damping parameter, ai: sij =

aij 1/3

(αij)

=

(ai + aj)/2 [(αiαj)1/2 ]1/3

Figure 1. Illustration of the contribution of pair style thole. The full line represents the Coulomb potential computed by LAMMPS between the charges ±q corresponding to the dipoles. It can be decomposed into a short-range contribution (direct space, shaded area) and a long-range contribution (reciprocal space) computed using Ewald summation techniques. The dashed line represents the non-diverging Thole potential. Note that the difference between Coulomb and Thole potential is only at short range. The vertical line represents the cutoff distance of the Thole correction and of the short-range part of Coulomb interaction (direct space). The contribution of the pair style thole is the (negative) difference between the dashed and full lines.

(3)

The mixing rule for polarizabilities is the geometric mean between the atom-specific values for Thole damping parameters is the arithmetic mean. A suitable value is a = 2.6 for all atom types, but in certain force fields that value can depend upon the atom type. Thole’s screening is applied to all Drude oscillators whatever their special neighbor relation. It is computed using only the charges ±q corresponding to the induced dipoles, i.e., the charge of a DP and the opposite charge on the respective DC (and not the total charge of the DC). LAMMPS allows for the definition of mixed pair styles which can be added up through the hybrid/overlay option of the pair_style command. We implemented a pair style called thole as a correction to standard Coulomb interaction, so that their sum yields the Thole’s screened interaction. This pair style is to be used as hybrid/overlay with any LAMMPS pair style implementing Coulomb interactions, including long-range styles using Ewald summation, PPPM, etc. Internally, for each pair of interacting “atoms”, if both atoms are DP or DC, then the standard (or screened by 1−2 to 1−4 factors) Coulomb interaction between the ±q charges is subtracted, and Thole’s interaction is added instead. The situation with Ewald summation is illustrated in Figure 1. To summarize, the parameters required for Thole’s interaction are the following: 1. atom polarizability αi in units of length cube (the same used for the calculation of the Drude charge or the Drude force constant) 2. Thole damping parameter ai, by default 2.6 3. cutoff radius, since the Thole and Coulomb interactions only differ at short range This pair style thole is handy because it can be used in combination with any other Coulomb pair style. But it is slightly slower than a direct integration with existing pair styles. We also implemented an integrated version for the pair style lj/cut/ coul/long (short-range Lennard-Jones + long-range Coulomb) called lj/cut/thole/long. The savings in computational time are of the order of 10% (see Performance section). Thermostating the Drude Oscillators. The degrees of freedom corresponding to the Drude oscillators, i.e., the motion of the DPs relative to the respective DCs, should be thermalized at a low temperature TD, of a few K. This avoids fast vibrations of the small reduced masses (allowing a reasonable time step), and

the dynamics of the system approximates the self-consistent regime.16 We implemented two thermostating mehods, using either a Nosé−Hoover16 or a Langevin30 thermostat. Both thermostats act on the reduced degrees of freedom, which are defined for each Drude pair16 by reduced masses, positions, velocities, and forces: M′ = M + m ,

m′ =

Mm M′

(4)

X′ =

MX + mx , M′

x′ = x − X

V′ =

MV + mv , M′

v′ = v − V

(5) (6)

Mf − mF (7) M′ In these equations, upper case denotes DC or center-of-mass values, and lower case denotes DP or dipole values. Primes denote the transformed (reduced) values, while bare letters denote the original values. This transform conserves the total kinetic energy and the virial defined with absolute positions: 1 1 (MV 2 + mv 2) = (M′V ′2 + m′v′2 ) (8) 2 2 F′ = F + f ,

f′ =

XF + xf = X ′F ′ + x′f ′

(9)

Internally, the Langevin thermostat is invoked after the force calculation. It works by adding new contributions to the forces. The Langevin forces for target temperatures TC (centers of mass) and TD (Drude oscillators) are computed as F′ = − C

M′ V ′ + F ′r τC

(10) DOI: 10.1021/acs.jcim.5b00612 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling f′ = −

m′ v′ + f ′r τD

all types of DP, and calculating the Drude charges for individual atom types from the atom polarizabilities using eq 16. This choice is followed in the polarizable CHARMM force field. • Schröder and Steinhauser35 suggest adopting universal charge qD = −1.0 e and mass mD = 0.1 g mol−1 (or u) for all DPs, and calculating the force constant for each type of DC−DP bond from eq 16. The time steps used by these authors are between 0.5 and 2 fs, with the degrees of freedom of the Drude oscillators kept cold at 1 K. In both of these force fields, hydrogen atoms are treated as non-polarizable. We implemented a Python script called polarizer distributed with the USER-DRUDE package to convert a nonpolarizable data file (here data.lmp) into a polarizable data file, data-p.lmp. It can be invoked as

(11)

where F′r and f ′r are Gaussian random forces characterized by zero mean and a variance of 2kBTCM′/(dtτC) and 2kBTDm′/(dtτD), respectively. The real forces acting on the particles are then computed from the inverse transform: M F′ − f ′ M′ m f= F′ + f ′ M′

F=

(12) (13)

The Langevin thermostat for the Drude degrees of freedom is implemented in the fix langevin/drude command. Optionally, the total Langevin force on centers of mass can be canceled by subtracting an equal part of it from all the forces F′. This avoids a random drift of the total center of mass of the system. The Nosé−Hoover thermostat uses the existing thermostats in LAMMPS. Just before the velocity-Verlet half-steps updating positions and velocities, the transform eqs 4−7 is performed. Just after the half-steps, the inverse transform is invoked. So in particular the computation of forces, the building of pair lists, and the trajectory output are performed in the normal representation. The commands for the direct and inverse transforms implemented by the Drude package are called fix drude/ transform/direct and fix drude/transform/ inverse.

This will automatically insert the new atoms and bonds and the Drude Types section. The masses and charges of DCs and DPs are computed from information in the file phenol.dff, as well as the DC−DP bond constants. The file phenol.dff contains the polarizabilities of the atom types and the masses of the DPs, for instance:

■

LAMMPS INPUT FILES In this section we explain how to run a simulation with the Drude package, using phenol as the example. We assume that the user has available the input files for the corresponding non-polarizable system. Several molecular builder tools are available to prepare LAMMPS input files for systems composed of molecules or materials, such as VMD TopoTools,31 Moltemplate,32 or fftool.33 First, LAMMPS has to be compiled with the USERDRUDE package activated. Then, the data file and input scripts have to be modified to include the Drude dipoles and how to handle them, tasks that can be simplified by using the polarizer tool supplied with this package. Data File. Each polarizable atom should be split into a pair of atoms DC and DP representing the Drude oscillator. The sum of the masses of each DC and DP should equal the mass of the initial atom, M + m = Ma (14)

Hydrogen atoms are absent from this file so they would be treated as non-polarizable atoms. LAMMPS does not store atom names internally, but rather numerical atom type IDs. Therefore, to establish the correspondence, in the non-polarizable data file data.lmp atom names corresponding to the atom type numbers have to be specified as comments at the end of lines of the Masses section, as in

Those comments with the atom type names are read by the polarizer script but are ignored by LAMMPS. In this example, atom types 1, 2, and 3 will be polarizable atoms, whereas atom types 4 and 5 will be non-polarizable. Therefore, three types of Drude particle will be added, and the polarizable system will have eight atom types. The masses and charges of all the atoms of the three polarizable types will be modified to the values of the corresponding Drude cores; a new Drude particle of the required type will be added for each core, and new bonds will be created between cores and Drude particles. Basic Input Script. The polarizer tool also outputs certain lines related to the input script (the use of these lines will be explained below). In order for LAMMPS to handle atoms which are bonded and carry electrostatic charge, the atom style full should be used, with harmonic bonds.

and similarly, the sum of the charges of DC and DP should equal the charge of the initial atom, Q + q = Qa

(15)

DC−DP bonds should appear as standard bonds in the Bonds section of the data file. They should have zero equilibrium length. The stiffness of the harmonic bond kD34 and the Drude charge q are related to the atom polarizability α by4 kD = q 2 / α

(16)

The values of Drude mass, Drude charge, and force constant can be chosen following different strategies, as in the following examples of polarizable force fields: • Lamoureux and Roux16 suggest adopting a global halfstiffness, kD = 4184 kJ mol−1 Å−2, for all types of DC−DP bond, and a universal mass, mD = 0.4 g mol−1 (or u), for

Electrostatic interactions are of course required for modeling polarizability. So the pair style must include Coulomb interactions, for instance lj/cut/coul/long with kspace D


Article


modification of forces but no position/velocity updates), fix nve should be used in conjunction:

style pppm. For example, with a cutoff of 12.0 and a precision 10−4:

Note that the global temperature thermo_temp computed by LAMMPS is not 300.0 K as desired. This is because LAMMPS treats DPs as standard atoms in its default compute. In order to output the temperatures of the DC−DP pair centers of mass and of the DPs relative to their DCs, the compute temp/drude should be used:

The scaling factor for 1−2 interactions should be set to zero using the special_bonds command. Also, it may be necessary to extend the space for storing such special relations, in virtue of the additional Drude particles. In this case extra space should be reserved by using the extra keyword of the special_bonds command. For phenol, there is one more special neighbor for which space is required. If omitted, LAMMPS will exit with a message indicating the required value.

Then the thermo_style custom can be invoked with respectively c_TDRUDE[1] and c_TDRUDE[2]. These should be close to 300.0 and 1.0 on average, respectively, in this example:

Pair interactions involving the non-polarizable atoms, cores, and Drudes can be conveniently supplied in the input script. As compared to the non-polarizable input file, pair_coeff lines need to be added for the DPs. Since the DPs (supposedly atom type 6 and above in our example) have no Lennard-Jones interactions, their epsilon is 0.0, so the only pair_coeff line that needs to be added is

In the phenol example, fix shake is used to constrain C−H and O−H bonds. As explained above, it must be applied to a group-ID not including the Drude particles; for example, it can be invoked on the ATOMS group-ID:

For several purposes, such as applying fixes for thermostating and barostating the system, among other uses, it is convenient to define atom groups for atoms + DCs, DCs, and DPs. For phenol, such atom groups would have been defined as

To summarize, the block of specific commands in the LAMMPS input script (except pair style infomation) to simulate phenol using cold Drude dipoles, with the Langevin thermostat, could be

The new fix drude command specifies if an atom is • tag N, a non-polarizable atom • tag C, a polarizable core • tag D, a Drude particle In the example of phenol, this information would be given as Thole Screening. As explained above, dipolar interactions represented by point charges on springs may not be stable. For example if the atomic polarizability is too high a DP can escape from its DC and be captured by another DC, which makes the force and energy diverge and the simulation crash. Even without reaching this extreme case, the correlation between nearby dipoles on the same molecule may be exaggerated. Often, special bond relations prevent bonded neighboring atoms to interact with the charge of each other’s DP, so that the problem does not always appear. It is possible to use screened dipole−dipole interactions through the pair_style thole. This is implemented as a correction to the Coulomb pair styles, which dampens at short distance the interactions between the charges representing the induced dipoles. It is to be used as hybrid/overlay with any standard coul pair style. For phenol, we could use

Let us consider first how to run a simple NVT simulation at 300 K using real units. Since Drude oscillators need to be thermalized at a low temperature in order to approximate selfconsistent dynamics, it is not possible to simulate an NVE ensemble with this package. For the thermalization, the simplest choice is to use the fix langevin/drude, This applies a Langevin thermostat at temperature TC = 300.0 to the centers of mass of the DC−DP pairs, with relaxation time τC = 100 and with random seed 12 345. This fix applies also a Langevin thermostat at temperature TD = 1.0 to the relative motion of the DPs around their DCs, with relaxation time τD = 20 and random seed 13 977. Only the DCs and non-polarizable atoms need to be in this fix’s group. LAMMPS will thermostat the DPs together with their DC. For this, ghost atoms need to know their velocities. Thus, the following command has to be included:

This tells LAMMPS that two pair styles are combined. The first is as above (lj/cut/coul/long 12.0), and the second is a thole pair style with default screening factor a = 2.6 and cutoff 12.0. Since hybrid/overlay does not support mixing rules in LAMMPS, the interaction coefficients of all the pairs of atom types with I < J should be explicitly defined. The output of the polarizer script can be used to complete the

The center of mass of the whole system may drift due to the random forces of the Langevin thermostat on DCs, but this drift can be avoided by adding the zero yes option at the end of the fix langevin/drude line. If fix shake is used to constrain bonds or angles, it should preferably be added after the fix langevin/drude for more accuracy. Since the fix langevin/drude does not perform time integration (just E


Article


The fix drude/transform/inverse performs the reverse transformation. For a NVT simulation with the DCs and atoms at 300 K and the DPs at 1 K relative to their DC one could use the following block of commands:

pair_coeff section of the input file. For phenol, this could look like

To avoid the flying ice cube artifact,16 where the atoms progressively lose kinetic energy and the center of mass of the whole system drifts, the fix momentum can be invoked. For instance, Running an NpT simulation with Nosé−Hoover barostat and thermostat is slightly more involved. First, the volume should be integrated only once. Therefore, DCs and non-polarizable atoms should be regulated by fix npt while DPs should be regulated by fix nvt (or vice versa). Second, the fix npt computes a global pressure and thus a global temperature whatever the fix group. We do want the pressure to correspond to the whole system, but we want the temperature to correspond to the fix group only. We must thus use the fix modify command for this. In the end, the block of instructions for a NpT simulation would look like

For the thole pair style the coefficients are the atom polarizability α, the Thole damping parameter a (optional, default value specified by the pair_style command), and the cutoff (optional, default value defined by the pair_style command). The special neighbors have Coulomb interactions screened by the coul factors of the special_bonds command (in the example above 0.0, 0.0, and 0.5 for 1−2, 1−3, and 1−4 special neighbors, respectively), so without using the pair_style thole, dipole−dipole interactions are screened by the same factor. By using pair_style thole, dipole−dipole interactions are screened by Thole functions, whatever their special relationship (except within each DC−DP pair of course). Consider for example 1−2 neighbors: using pair_style thole the dipoles will interact (despite the coul special factor being 0.0), and the interactions between these dipoles will be damped by a Thole function. Thermostats and Barostats. The use of a Langevin thermostat for the Drude oscillators coupled with a constantvolume integrator has been documented in the section above describing a basic input script. To perform a simulation at constant NpT, a Nosé−Hoover barostat can be combined with the langevin/drude thermostat in a straightforward manner, by using fix nph instead of fix nve. For example, a trajectory at a constant pressure of 1.0 bar can be generated using

Similarly to the preceding case where the Langevin thermostat was used, the correct temperatures of cores and Drude particles, using center-of-mass velocities and relative velocities, respectively, can be calculated using the compute temp/drude command. Rigid Bodies. It is possible to simulate molecules as rigid bodies (but polarizable). Common cases handled as rigid bodies are water models with structures similar to that of TIP4P (which has a mass-less charged site at a distance from the O atom bisecting the H−O−H angle). Examples are SWM4NDP,36 which has a Drude dipole on the O atom, and COS/ G3,37 which has a Drude dipole on the mass-less charge site. The rigid bodies and the DPs should be time-integrated separately, even with the Langevin thermostat. Let us review the different thermostats and ensemble combinations. NVT ensemble using Langevin thermostat:

NVT ensemble using Nosé−Hoover thermostat:

NpT ensemble with Langevin thermostat:

It is also possible to use a Nosé−Hoover thermostat instead of a Langevin thermostat. This requires the command fix drude/ transform to be invoked just before and after the time integration fixes. The fix drude/transform/direct converts the atomic masses, positions, velocities and forces into a reduced representation, where the DCs transform into the centers of mass of the DC−DP pairs and the DPs transform into their relative position with respect to their DC.

NpT ensemble with Nosé−Hoover thermostat:

F


Article


■

VALIDATION In this section the results of using the USER-DRUDE package are compared with those of literature sources for different systems which are well documented. We chose the following: • PSPC water,16 a three-site water model with one Drude dipole on the O atom, in which the Drude particle has a positive charge. PSPC water can be simulated with constrained bonds and angle in LAMMPS, using fix shake. • SWM4-NDP water,36 a four-site water model with one Drude dipole on the O atom, in which the Drude particle has a negative charge. Due to the mass-less fourth site carrying a partial charge, this water model is simulated as a rigid body. • Ethanol29 with three Drude dipoles on the C and O atoms. Thole functions are used in this example. The simulation conditions for the three test systems, whose force field parameters were taken from the sources listed above, are given in Table 1. The DC−DP force constant was fixed at

parameters and electrostatic terms, including Drude parameters, from Noskov et al.,29 combined with intramolecular terms taken from the OPLS force field (for simplicity). Since NAMD only allows Langevin thermostating for Drude oscillators, we only compare the results obtained with the Langevin thermostat, using relaxation times of 100 fs for DCs and 20 fs for DPs. The properties deduced from NAMD and LAMMPS trajectories are given in Table 3. In all cases, the quantities calculated include Table 3. Comparison of the Results Obtained for Ethanol Using the Drude Dipole Implementation in LAMMPS and in NAMD

Table 1. Simulation Conditions for the Test Systems Used in the Validation of the Drude Dipole Implementation in LAMMPS N molecules cutoff radius (Å) PPPM tolerance temperature (K) Drude temp (K) pressure (bar) time step (fs) trajectory (ns)

PSPC

SWM4-NDP

Ethanol

500 12.0 1.0 × 10−4 300.0 1.0 1.0 1.0 1.0

500 12.0 1.0 × 10−4 300.0 1.0 1.0 1.0 1.0

250 12.0 1.0 × 10−4 300.0 1.0 1.0 1.0 3.0

Nosé−Hoover

literature

v (Å3) Δu (kJ mol−1) D (×10−9 m2 s−1) ϵ

31.7 ± 0.4 33.8 ± 0.3 3.97 ± 0.01 183 ± 30

32.12 33.22 3.97 ± 0.08 185 ± 30

v (Å3) Δu (kJ mol−1) D (×10−9 m2 s−1) ϵ

SWM4-NDP 30.2 ± 0.3 30.2 ± 0.3 41.2 ± 0.3 41.2 ± 0.3 1.21 ± 0.01 2.34 ± 0.01 76 ± 11 81 ± 13

29.91 41.52 2.33 ± 0.02 79 ± 3

v (Å ) Δu (kJ mol−1) rO‑O (Å) rO‑HO (Å)

93.7 ± 1.2 40.0 ± 1.1 2.67 ± 0.01 1.75 ± 0.01

94.1 ± 1.3 40.4 ± 1.1 2.66 ± 0.01 1.75 ± 0.01

D (×10−9 m2 s−1) ⟨μ⟩ (D) ϵ

0.27 ± 0.01 2.38 ± 0.02 22.5 ± 5.7

0.27 ± 0.01 2.41 ± 0.02 21.4 ± 4.6

■

PERFORMANCE Thermalized Drude oscillators allow for as long an integration time step as non-polarizable systems. However, the polarizable simulations remain more time-consuming than their nonpolarizable equivalents, for the following reasons: 1 New degrees of freedom and interactions are introduced (DP’s positions, velocities, DC−DP bonds). At least, pair interactions scale with the number density of “atoms” squared. Here, DPs are introduced at constant volume, so the number density increases (even though the mass density is constant). 2 Extra computations are required, notably for Thole screening and real/reduced coordinate transforms. Typically, without the OMP package the simulation CPU time is 2−3 times longer for polarizable systems using Drude dipoles than for the non-polarizable ones. In Table 4 and Figure 2 we compare the simulation times in the case of ethanol, with 250 ethanol molecules simulated during 500 000 steps in the NVT ensemble on 12 processors using MPI communication. The time for computing non-bonded interactions (Pair+Kspce) is largely increased with polarizable models. Using the combined Thole screening allows for a small acceleration with respect to using a specific pair style for Thole screening, which is combined via hybrid/overlay with a Coulomb term. Bonded interactions,

Table 2. Comparison of the Results Obtained Using the Drude Dipole Implementation in LAMMPS, with Both Langevin and Nosé−Hoover Thermostats, with Values from the Literature16,36 PSPC 31.6 ± 0.4 33.8 ± 0.3 1.59 ± 0.01 192 ± 36

NAMD

molecular volume v, molecular nearest-neighbor distances rA‑B, diffusion coefficient D, relative permittivity ϵ, and molecular dipole moment μ. The diffusion coefficients were calculated from the mean-squared displacement and the relative permittivity from the fluctuations of the total dipole moment of the simulation box. The comparisons show a good agreement for all quantities except for the diffusion coefficient of water using a Langevin thermostat. This is a known issue of Langevin thermostats, which tend to slow dynamics down due to artificial friction, so this is not related to the implementation of Drude oscillators. Therefore, these comparisons validate the implementation in LAMMPS for molecules containing one or more Drude oscillators (use of Thole damping), and also for rigid (using constrains or rigidbody integrators) and flexible molecules.

kD = 4184.0 kJ mol−1 Å−2, and the mass of the DPs was m = 0.4 g mol−1. Simulations were performed in periodic cubic boxes, with long-range electrostatics calculated using the PPPM Ewald summation. The properties calculated for the test systems are compared to values from literature in Table 2. For ethanol the level of details

Langevin

LAMMPS 3

reported in the literature29 is not sufficient for a thorough comparison, therefore we chose to test the LAMMPS implementation against reference simulations on the same force field model performed with NAMD.26 We used the Lennard-Jones G


Article

Journal of Chemical Information and Modeling Table 4. Comparison of Computation Times in Secondsa NP−NH NP−L P−NH−HO P−L−HO P−NH−C P−L−C a

Total

Pair

Bond

Kspce

Neigh

Comm

Other

2149.0 2112.0 7265.0 5919.0 6215.0 4975.0

1280.0 1268.2 3748.6 3748.6 3033.9 2993.4

177.34 175.18 183.12 180.16 182.62 182.62

342.81 333.49 1214.14 1113.78 933.10 976.28

76.03 71.00 154.92 145.21 123.28 117.72

163.10 151.96 333.84 389.32 333.48 373.68

109.25 111.75 1597.36 341.84 1608.71 331.87

NP/P, non-polarizable/polarizable; NH/L, Nosé−Hoover/Langevin; HO/C, hybrid overlay/combined. See Figure 2 for the other notations.

throughout, allowing users to understand the function of the different pieces and how they need to be articulated to make the model work. Comparisons with published results and with simulations made with another molecular dynamics software validate the implementation. Performance analysis shows that the cost of adding explicit polarization through Drude dipoles is about 3-fold in terms of CPU time when compared to a nonpolarizable model, for a typical molecular fluid, mainly because of the computation of pair and electrostatic interactions involving the additional particles of the Drude oscillators. The package has been included in the standard LAMMPS version. It requires the command make YES-DRUDE when compiling LAMMPS. Some of the examples used in the paper including the polarizability script can be obtained with the LAMMPS source code.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected].

Figure 2. Comparison of computation times for simulations of 250 ethanol molecules for 500 000 steps on 12 processors (MPI) using non-polarizable model, polarizable model with hybrid/overlay Thole screening, and polarizable model with integrated Thole screening. Langevin and Nosé−Hoover thermostats are tested in each case. Pair = Lennard-Jones + short-range Coulomb including Thole screening; Bond = bonds, angles, and dihedrals; Kspce = long-range Coulomb (PPPM); Neigh = pair list constructions; Comm = communications, Outpt = output; and Other = time integration, transforms, force modifications, etc.

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Steve Plimpton and Axel Kohlmeyer for their help regarding implementation and integration in the official LAMMPS version.

■

neighbor list constructions and output have a minimal contribution. The Nosé−Hoover thermostat is significantly slower than the Langevin thermostat.

REFERENCES

(1) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (2) Tkatchenko, A.; Scheffler, M. Accurate Molecular Van Der Waals Interactions from Ground-State Electron Density and Free-Atom Reference Data. Phys. Rev. Lett. 2009, 102, 073005. (3) Plimpton, S. J. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. See also http://lammps. sandia.gov/. (4) Rick, S. W.; Stuart, S. J. Potentials and algorithms for incorporating polarizability in computer simulations; Reviews in Computational Chemistry 18; John Wiley & Sons Inc.: Hoboken, NJ, 2002; pp 89−146. (5) Onsager, L. Electric Moments of Molecules in Liquids. J. Am. Chem. Soc. 1936, 58, 1486−1493. (6) Kirkwood, J. G. The Dielectric Polarization of Polar Liquids. J. Chem. Phys. 1939, 7, 911. (7) Stillinger, F. H.; David, C. W. Polarization Model for Water and its Ionic Dissociation Products. J. Chem. Phys. 1978, 69, 1473−1484. (8) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (9) Sprik, M.; Klein, M. L. A Polarizable Model for Water Using Distributed Charge Sites. J. Chem. Phys. 1988, 89, 7556. (10) Nalewajski, R. F. A. Study of Electronegativity Equalization. J. Phys. Chem. 1985, 89, 2831−2837.

■

CONCLUSION We have introduced a new user-contributed package for LAMMPS, called USER-DRUDE, intended for simulation of polarizable systems with thermalized Drude oscillators. In our view this is an interesting addition to the LAMMPS code, allowing a description of explicit polarization in molecular or ionic fluids, and also materials. This package is complementary to other existing modules in LAMMPS, such as the QEQ package (charge equilibration) or the CORESHELL package (adiabatic core−shell model suitable for ionic solids). With the USERDRUDE implementation, Drude oscillators can be thermalized using a Nosé−Hoover or a Langevin thermostat, and shortrange screening functions are provided for covalent systems. The Drude oscillators can be used in flexible or rigid molecules (the latter using constrains or rigid-body integrators). The use of Drude dipoles in LAMMPS requires the combination of a number of elements (thermostats, screening functions, parametrization of the Drude dipoles themselves) and a successful use demands attention. In order to provide clear and detailed guidance, examples have been provided and commented H


Article

Journal of Chemical Information and Modeling (11) Rappe, A. K.; Goddard, W. A. Charge Equilibration for Molecular Dynamics Simulations. J. Phys. Chem. 1991, 95, 3358−3363. (12) Rick, S. W.; Stuart, S. J.; Berne, B. J. Dynamical Fluctuating Charge Force Fields: Application to Liquid Water. J. Chem. Phys. 1994, 101, 6141−6156. (13) Smirnov, K. S.; Graaf, B. v. d. Consistent Implementation of the Electronegativity Equalization Method in Molecular Mechanics and Molecular Dynamics. J. Chem. Soc., Faraday Trans. 1996, 92, 2469− 2474. (14) Drude, P. The theory of optics; Lehrbuch der Optik; Longmans, Green, and Co.: London/New York, 1902. (15) van Maaren, P. J.; van der Spoel, D. Molecular Dynamics Simulations of Water with Novel Shell-Model Potentials. J. Phys. Chem. B 2001, 105, 2618−2626. (16) Lamoureux, G.; Roux, B. Modeling Induced Polarization with Classical Drude Oscillators: Theory and Molecular Dynamics Simulation Algorithm. J. Chem. Phys. 2003, 119, 3025−3039. (17) Lopes, P. E. M.; Roux, B.; Mackerell, A. D. Molecular Modeling and Dynamics Studies with Explicit Inclusion of Electronic Polarizability: Theory and Applications. Theor. Chem. Acc. 2009, 124, 11−28. (18) Huang, J.; Lopes, P. E. M.; Roux, B.; MacKerell, A. D., Jr. Recent Advances in Polarizable Force Fields for Macromolecules: Microsecond Simulations of Proteins Using the Classical Drude Oscillator Model. J. Phys. Chem. Lett. 2014, 5, 3144−3150. (19) Bauer, B. A.; Patel, S. Recent Applications and Developments of Charge Equilibration Force Fields for Modeling Dynamical Charges in Classical Molecular Dynamics Simulations. Theor. Chem. Acc. 2012, 131, 1153−1155. (20) Salomon-Ferrer, R.; Case, D. A.; Walker, R. C. An Overview of the Amber Biomolecular Simulation Package. WIRES: Comput. Mol. Sci. 2013, 3, 198−210. (21) Baker, C. M.; Best, R. B. Matching of Additive and Polarizable Force Fields for Multiscale Condensed Phase Simulations. J. Chem. Theory Comput. 2013, 9, 2826−2837. (22) In their present implementation in LAMMPS, the QEQ methods are meant to treat all of the atoms in the system and are not prepared to mix polarizable and fixed-charged groups of atoms in the same simulation. (23) Mitchell, P. J.; Fincham, D. Shell Model Simulations by Adiabatic Dynamics. J. Phys.: Condens. Matter 1993, 5, 1031−1038. (24) Todorov, I. T.; Smith, W.; Trachenko, K.; Dove, M. T. DL_POLY_3: New Dimensions in Molecular Dynamics Simulations via Massive Parallelism. J. Mater. Chem. 2006, 16, 1911. See also http:// www.scd.stfc.ac.uk/SCD/research/app/44516.aspx. (25) Brooks, B. R.; Brooks, C. L.; Mackerell, A. D.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; Caflisch, A.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer, M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.; York, D. M.; Karplus, M. CHARMM: The Biomolecular Simulation Program. J. Comput. Chem. 2009, 30, 1545−1614. See also http://www.charmm. org/. (26) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. See also http://www.ks.uiuc.edu/Research/namd/. (27) Harder, E.; Anisimov, V. M.; Vorobyov, I. V.; Lopes, P. E. M.; Noskov, S. Y.; MacKerell, A. D.; Roux, B. Atomic Level Anisotropy in the Electrostatic Modeling of Lone Pairs for a Polarizable Force Field Based on the Classical Drude Oscillator. J. Chem. Theory Comput. 2006, 2, 1587−1597. (28) Thole, B. T. Molecular Polarizabilities Calculated with a Modified Dipole Interaction. Chem. Phys. 1981, 59, 341−350. (29) Noskov, S. Y.; Lamoureux, G.; Roux, B. Molecular Dynamics Study of Hydration in Ethanol-Water Mixtures Using a Polarizable Force Field. J. Phys. Chem. B 2005, 109, 6705−6713. (30) Jiang, W.; Hardy, D. J.; Phillips, J. C.; MacKerell, A. D.; Schulten, K.; Roux, B. High-Performance Scalable Molecular Dynamics

Simulations of a Polarizable Force Field Based on Classical Drude Oscillators in NAMD. J. Phys. Chem. Lett. 2011, 2, 87−92. (31) Kohlmeyer, A. TopoTools, 2014; https://sites.google.com/site/ akohlmey/software/topotools. (32) Jewett, A.; Lambert, J. Moltemplate, 2014; http://www. moltemplate.org/. (33) Pádua, A. A. H. fftool: a tool to build force field input files for molecular dynamics, 2015; https://github.com/agiliopadua/fftool. (34) Many molecular dynamics codes use K = k/2 as force constants in harmonic bonds and angles. (35) Schröder, C.; Steinhauser, O. Simulating Polarizable Molecular Ionic Liquids with Drude Oscillators. J. Chem. Phys. 2010, 133, 154511. (36) Lamoureux, G.; Harder, E.; Vorobyov, I. V.; Roux, B.; MacKerell, A. D. A Polarizable Model of Water for Molecular Dynamics Simulations of Biomolecules. Chem. Phys. Lett. 2006, 418, 245−249. (37) Yu, H.; van Gunsteren, W. F. Charge-on-Spring Polarizable Water Models Revisited: From Water Clusters to Liquid Water to Ice. J. Chem. Phys. 2004, 121, 9549.

I


Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics

Recommend Documents