These are additional pair styles in USER packages, which can be used @@ -530,7 +532,8 @@ package.
Each listed input is operated on independently. For per-atom inputs, the group specified with this command means only atoms within the diff --git a/doc/compute_reduce.txt b/doc/compute_reduce.txt index c4cf88356a..0074e990f3 100644 --- a/doc/compute_reduce.txt +++ b/doc/compute_reduce.txt @@ -18,7 +18,7 @@ style = {reduce} or {reduce/region} :l {reduce} arg = none {reduce/region} arg = region-ID region-ID = ID of region to use for choosing atoms :pre -mode = {sum} or {min} or {max} or {ave} :l +mode = {sum} or {min} or {max} or {ave} or {sumsq} or {avesq} :l one or more inputs can be listed :l input = x, y, z, vx, vy, vz, fx, fy, fz, c_ID, c_ID\[N\], f_ID, f_ID\[N\], v_name :l x,y,z,vx,vy,vz,fx,fy,fz = atom attribute (position, velocity, force component) @@ -58,7 +58,12 @@ The reduction operation is specified by the {mode} setting. The {sum} option adds the values in the vector into a global total. The {min} or {max} options find the minimum or maximum value across all vector values. The {ave} setting adds the vector values into a global total, -then divides by the number of values in the vector. +then divides by the number of values in the vector. The {sumsq} +option sums the square of the values in the vector into a global +total. The {avesq} setting does the same as {sumsq}, then divdes the +sum of squares by the number of values. The last two options can be +useful for calculating the variance of some quantity, e.g. variance = +sumsq - ave^2. Each listed input is operated on independently. For per-atom inputs, the group specified with this command means only atoms within the diff --git a/doc/compute_temp_drude.html b/doc/compute_temp_drude.html new file mode 100644 index 0000000000..6e50495d5a --- /dev/null +++ b/doc/compute_temp_drude.html @@ -0,0 +1,62 @@ + +
Syntax: +
+compute ID group-ID temp/drude ++
Examples: +
+compute TDRUDE all temp/drude ++
Description: +
+Define a computation that calculates the temperature based on the +center-of-mass velocities of pairs of Drude cores and Drude particles, +bonded by springs. This compute is designed to be used with the +thermalized Drude oscillator model. Polarizable +models in LAMMPS are described in this +Section. +
+Specifically, this compute enables calculation of the temperature of +the Drude particles in relative coordinates with respect to their +cores. +
+Output info: +
+This compute calculates a global scalar (the temperature) and a global +vector of length 6 (KE tensor), which can be accessed by indices 1-6. +These values can be used by any command that uses global scalar or +vector values from a compute as input. See this +section for an overview of LAMMPS output +options. +
+The scalar value calculated by this compute is "intensive". The +vector are "extensive". +
+The scalar value will be in temperature units. The +vector values will be in energy units. +
+Restrictions: +
+The number of core-Drude pairs contributing to the temperature is +assumed to be constant for the duration of the run. +
+Related commands: none +
+Default: none +
+ diff --git a/doc/compute_temp_drude.txt b/doc/compute_temp_drude.txt new file mode 100644 index 0000000000..0b36fc5431 --- /dev/null +++ b/doc/compute_temp_drude.txt @@ -0,0 +1,57 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +compute temp/drude command :h3 + +[Syntax:] + +compute ID group-ID temp/drude :pre + +ID, group-ID are documented in "compute"_compute.html command +temp/drude = style name of this compute command :ul + +[Examples:] + +compute TDRUDE all temp/drude :pre + +[Description:] + +Define a computation that calculates the temperature based on the +center-of-mass velocities of pairs of Drude cores and Drude particles, +bonded by springs. This compute is designed to be used with the +"thermalized Drude oscillator model"_tutorial_drude.html. Polarizable +models in LAMMPS are described in "this +Section"_Section_howto.html#howto_25. + +Specifically, this compute enables calculation of the temperature of +the Drude particles in relative coordinates with respect to their +cores. + +[Output info:] + +This compute calculates a global scalar (the temperature) and a global +vector of length 6 (KE tensor), which can be accessed by indices 1-6. +These values can be used by any command that uses global scalar or +vector values from a compute as input. See "this +section"_Section_howto.html#howto_15 for an overview of LAMMPS output +options. + +The scalar value calculated by this compute is "intensive". The +vector are "extensive". + +The scalar value will be in temperature "units"_units.html. The +vector values will be in energy "units"_units.html. + +[Restrictions:] + +The number of core-Drude pairs contributing to the temperature is +assumed to be constant for the duration of the run. + +[Related commands:] none + +[Default:] none diff --git a/doc/fix_drude.html b/doc/fix_drude.html new file mode 100644 index 0000000000..398abd701e --- /dev/null +++ b/doc/fix_drude.html @@ -0,0 +1,60 @@ + +Syntax: +
+fix ID group-ID drude flag1 flag2 ... flagN ++
Examples: +
+fix 1 all drude 1 1 0 1 0 2 2 2 +fix 1 all drude C C N C N D D D ++
Description: +
+Assign each atom type in the system to be one of 3 kinds of atoms +within the Drude polarization model. This compute is designed to be +used with the thermalized Drude oscillator +model. Polarizable models in LAMMPS +are described in this Section. +
+The three possible types can be designated with an integer (0,1,2) +or capital letter (N,C,D): +
+Restrictions: +
+This fix should be invoked before any other commands that implement +the Drude oscillator model, such as fix +langevin_drude, fix +drude/transform, compute +temp/drude, pair_style +thole. +
+Related commands: +
+fix langevin_drude, fix +drude/transform, compute +temp/drude, pair_style +thole +
+Default: None +
+ diff --git a/doc/fix_drude.txt b/doc/fix_drude.txt new file mode 100644 index 0000000000..8c78fb0772 --- /dev/null +++ b/doc/fix_drude.txt @@ -0,0 +1,56 @@ + +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +fix drude command :h3 + +[Syntax:] + +fix ID group-ID drude flag1 flag2 ... flagN :pre + +ID, group-ID are documented in "fix"_fix.html command +drude = style name of this fix command +tag = Drude flag for each atom type (1 to N) in the system :ul + +[Examples:] + +fix 1 all drude 1 1 0 1 0 2 2 2 +fix 1 all drude C C N C N D D D :pre + +[Description:] + +Assign each atom type in the system to be one of 3 kinds of atoms +within the Drude polarization model. This compute is designed to be +used with the "thermalized Drude oscillator +model"_tutorial_drude.html. Polarizable models in LAMMPS +are described in "this Section"_Section_howto.html#howto_25. + +The three possible types can be designated with an integer (0,1,2) +or capital letter (N,C,D): + +0 or N = non-polarizable atom (not part of Drude model) +1 or C = Drude core +2 or D = Drude electron :ul + +[Restrictions:] + +This fix should be invoked before any other commands that implement +the Drude oscillator model, such as "fix +langevin_drude"_fix_langevin_drude.html, "fix +drude/transform"_fix_drude_transform.html, "compute +temp/drude"_compute_temp_drude.html, "pair_style +thole"_pair_thole.html. + +[Related commands:] + +"fix langevin_drude"_fix_langevin_drude.html, "fix +drude/transform"_fix_drude_transform.html, "compute +temp/drude"_compute_temp_drude.html, "pair_style +thole"_pair_thole.html + +[Default:] None diff --git a/doc/fix_drude_transform.html b/doc/fix_drude_transform.html new file mode 100644 index 0000000000..ee62f4726b --- /dev/null +++ b/doc/fix_drude_transform.html @@ -0,0 +1,195 @@ + + + + +Syntax: +
+fix ID group-ID style keyword value ... ++
temp value = yes or no = do or do not calculate reduced temperatures of core and Drude particles ++ +
Examples: +
+fix 3 all drude/transform/direct +fix 1 all drude/transform/inverse temp yes ++
Description: +
+Transform the coordinates of Drude oscillators from real to reduced +and back for thermalizing the Drude oscillators as described in +(Lamoureux) using a Nose-Hoover thermostat. This fix is +designed to be used with the thermalized Drude oscillator +model. Polarizable models in LAMMPS are +described in this Section. +
+Drude oscillators are a pair of atoms representing a single +polarizable atom. Ideally, the mass of Drude particles would vanish +and their positions would be determined self-consistently by iterative +minimization of the energy, the cores' positions being fixed. It is +however more efficient and it yields comparable results, if the Drude +oscillators (the motion of the Drude particle relative to the core) +are thermalized at a low temperature. In that case, the Drude +particles need a small mass. +
+The thermostats act on the reduced degrees of freedom, which are +defined by the following equations. Note that in these equations +upper case denotes atomic or center of mass values and lower case +denotes Drude particle or dipole values. Primes denote the transformed +(reduced) values, while bare letters denote the original values. +
+Masses: \begin{equation} M' = M + m \end{equation} +\begin{equation} m' = \frac {M\, m } {M'} \end{equation} +Positions: \begin{equation} X' = \frac {M\, X + m\, x} {M'} +\end{equation} \begin{equation} x' = x - X \end{equation} +Velocities: \begin{equation} V' = \frac {M\, V + m\, v} {M'} +\end{equation} \begin{equation} v' = v - V \end{equation} +Forces: \begin{equation} F' = F + f \end{equation} +\begin{equation} f' = \frac { M\, f - m\, F} {M'} +\end{equation} +
+This transform conserves the total kinetic energy +\begin{equation} \frac 1 2 \, (M\, V^2\ + m\, v^2) += \frac 1 2 \, (M'\, V'^2\ + m'\, v'^2) \end{equation} +and the virial defined with absolute positions +\begin{equation} X\, F + x\, f = X'\, F' + x'\, f' \end{equation} +
+The temp keyword specifies wheterh temperatures in reduced units for +cores and Drude particles are calculated. If the temp option is set +to yes the reduced temperatures, degrees of freedom, and kinetic +energies are calculated and can be accessed as explained below; +otherwise the reduced values are not calculated. +
+This fix requires each atom know whether it is a Drude particle or +not. You must therefore use the fix drude command to +specify the Drude status of each atom type. +
+IMPORTANT NOTE: only the Drude core atoms need to be in the group +specified for this fix. A Drude electron will be transformed together +with its cores even if it is not itself in the group. It is safe to +include Drude electrons or non-polarizable atoms in the group. The +non-polarizable atoms will simply not be transformed. +
+This fix does NOT perform time integration. It only transform masses, +coordinates, velocities and forces. Thus you must use separate time +integration fixes, like fix nve or fix +npt to actually update the velocities and positions of +atoms. In order to thermalize the reduced degrees of freedom at +different temperatures, two Nose-Hoover thermostats must be defined, +acting on two distinct groups. +
+IMPORTANT NOTE: The fix drude/transform/direct command must appear +before any Nose-Hoover thermostating fixes. The fix +drude/transform/inverse command must appear after any Nose-Hoover +thermostating fixes. +
+Example: +
+fix fDIRECT all drude/transform/direct +fix fNVT gCORES nvt temp 300.0 300.0 100.0 +fix fNVT gDRUDES nvt temp 1.0 1.0 100.0 +fix fINVERSE all drude/transform/inverse ++
In this example, gCORES is the group of the atom cores and gDRUDES +is the group of the Drude particles (electrons). The centers of mass +of the Drude oscillators will be thermostated at 300.0 and the +internal degrees of freedom will be thermostated at 1.0. +
+In addition, if you want to use a barostat to simulate a system at +constant pressure, only one of the Nose-Hoover fixes must be npt, +the other one should be nvt. You must add a compute temp/com and a +fix_modify command so that the temperature of the npt fix be just +that of its group but the pressure be the overall pressure +thermo_press. +
+Example: +compute cTEMP_CORE gCORES temp/com +fix fDIRECT all drude/transform/direct +fix fNPT gCORES npt temp 298.0 298.0 100.0 iso 1.0 1.0 500.0 +fix_modify fNPT temp cTEMP_CORE press thermo_press +fix fNVT gDRUDES nvt temp 5.0 5.0 100.0 +fix fINVERSE all drude/transform/inverse ++
In this example, gCORES is the group of the atom cores and gDRUDES +is the group of the Drude particles. The centers of mass of the Drude +oscillators will be thermostated at 298.0 and the internal degrees of +freedom will be thermostated at 5.0. The whole system will be +barostated at 1.0. +
+In order to avoid the flying ice cube problem (irreversible transfer +of linear momentum to the center of mass of the system), you may need +to add a fix momentum command like such as +
+fix fMOMENTUM all momentum 100 linear 1 1 1 ++
Restart, fix_modify, output, run start/stop, minimize info: +
+If the temp yes keyword is used in fix drude/transform/inverse +this fix computes a global vector with 6 components which can be +accessed by various output commands. The +meaning of the components are +
+No information about this fix is written to binary restart +files. +
+Restrictions: none +
+Related commands: +
+fix drude, +fix langevin/drude, +compute temp/drude, +pair_style thole +
+Default: +
+The option defaults are temp = no. +
+(Lamoureux) Lamoureux and Roux, J Chem Phys, 119, 3025-3039 (2003). +
+ diff --git a/doc/fix_drude_transform.txt b/doc/fix_drude_transform.txt new file mode 100644 index 0000000000..c795a96b3b --- /dev/null +++ b/doc/fix_drude_transform.txt @@ -0,0 +1,183 @@ + + + +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +fix drude/transform/direct command :h3 +fix drude/transform/inverse command :h3 + +[Syntax:] + +fix ID group-ID style keyword value ... :pre + +ID, group-ID are documented in "fix"_fix.html command :ulb,l +style = {drude/transform/direct} or {drude/transform/inverse} :l +zero or more keywords may be appended to drude/transform/inverse :l +keyword = {temp} :l + {temp} value = {yes} or {no} = do or do not calculate reduced temperatures of core and Drude particles :pre +:ule + +[Examples:] + +fix 3 all drude/transform/direct +fix 1 all drude/transform/inverse temp yes :pre + +[Description:] + +Transform the coordinates of Drude oscillators from real to reduced +and back for thermalizing the Drude oscillators as described in +"(Lamoureux)"_#Lamoureux using a Nose-Hoover thermostat. This fix is +designed to be used with the "thermalized Drude oscillator +model"_tutorial_drude.html. Polarizable models in LAMMPS are +described in "this Section"_Section_howto.html#howto_25. + +Drude oscillators are a pair of atoms representing a single +polarizable atom. Ideally, the mass of Drude particles would vanish +and their positions would be determined self-consistently by iterative +minimization of the energy, the cores' positions being fixed. It is +however more efficient and it yields comparable results, if the Drude +oscillators (the motion of the Drude particle relative to the core) +are thermalized at a low temperature. In that case, the Drude +particles need a small mass. + +The thermostats act on the reduced degrees of freedom, which are +defined by the following equations. Note that in these equations +upper case denotes atomic or center of mass values and lower case +denotes Drude particle or dipole values. Primes denote the transformed +(reduced) values, while bare letters denote the original values. + +Masses: \begin\{equation\} M' = M + m \end\{equation\} +\begin\{equation\} m' = \frac \{M\, m \} \{M'\} \end\{equation\} +Positions: \begin\{equation\} X' = \frac \{M\, X + m\, x\} \{M'\} +\end\{equation\} \begin\{equation\} x' = x - X \end\{equation\} +Velocities: \begin\{equation\} V' = \frac \{M\, V + m\, v\} \{M'\} +\end\{equation\} \begin\{equation\} v' = v - V \end\{equation\} +Forces: \begin\{equation\} F' = F + f \end\{equation\} +\begin\{equation\} f' = \frac \{ M\, f - m\, F\} \{M'\} +\end\{equation\} + +This transform conserves the total kinetic energy +\begin\{equation\} \frac 1 2 \, (M\, V^2\ + m\, v^2) += \frac 1 2 \, (M'\, V'^2\ + m'\, v'^2) \end\{equation\} +and the virial defined with absolute positions +\begin\{equation\} X\, F + x\, f = X'\, F' + x'\, f' \end\{equation\} + +The {temp} keyword specifies wheterh temperatures in reduced units for +cores and Drude particles are calculated. If the {temp} option is set +to {yes} the reduced temperatures, degrees of freedom, and kinetic +energies are calculated and can be accessed as explained below; +otherwise the reduced values are not calculated. + +:line + +This fix requires each atom know whether it is a Drude particle or +not. You must therefore use the "fix drude"_fix_drude.html command to +specify the Drude status of each atom type. + +IMPORTANT NOTE: only the Drude core atoms need to be in the group +specified for this fix. A Drude electron will be transformed together +with its cores even if it is not itself in the group. It is safe to +include Drude electrons or non-polarizable atoms in the group. The +non-polarizable atoms will simply not be transformed. + +:line + +This fix does NOT perform time integration. It only transform masses, +coordinates, velocities and forces. Thus you must use separate time +integration fixes, like "fix nve"_fix_nve.html or "fix +npt"_fix_nh.html to actually update the velocities and positions of +atoms. In order to thermalize the reduced degrees of freedom at +different temperatures, two Nose-Hoover thermostats must be defined, +acting on two distinct groups. + +IMPORTANT NOTE: The {fix drude/transform/direct} command must appear +before any Nose-Hoover thermostating fixes. The {fix +drude/transform/inverse} command must appear after any Nose-Hoover +thermostating fixes. + +Example: + +fix fDIRECT all drude/transform/direct +fix fNVT gCORES nvt temp 300.0 300.0 100.0 +fix fNVT gDRUDES nvt temp 1.0 1.0 100.0 +fix fINVERSE all drude/transform/inverse :pre + +In this example, {gCORES} is the group of the atom cores and {gDRUDES} +is the group of the Drude particles (electrons). The centers of mass +of the Drude oscillators will be thermostated at 300.0 and the +internal degrees of freedom will be thermostated at 1.0. + +In addition, if you want to use a barostat to simulate a system at +constant pressure, only one of the Nose-Hoover fixes must be {npt}, +the other one should be {nvt}. You must add a {compute temp/com} and a +{fix_modify} command so that the temperature of the {npt} fix be just +that of its group but the pressure be the overall pressure +{thermo_press}. + +Example: :b + +compute cTEMP_CORE gCORES temp/com +fix fDIRECT all drude/transform/direct +fix fNPT gCORES npt temp 298.0 298.0 100.0 iso 1.0 1.0 500.0 +fix_modify fNPT temp cTEMP_CORE press thermo_press +fix fNVT gDRUDES nvt temp 5.0 5.0 100.0 +fix fINVERSE all drude/transform/inverse :pre + +In this example, {gCORES} is the group of the atom cores and {gDRUDES} +is the group of the Drude particles. The centers of mass of the Drude +oscillators will be thermostated at 298.0 and the internal degrees of +freedom will be thermostated at 5.0. The whole system will be +barostated at 1.0. + +In order to avoid the flying ice cube problem (irreversible transfer +of linear momentum to the center of mass of the system), you may need +to add a {fix momentum} command like such as + +fix fMOMENTUM all momentum 100 linear 1 1 1 :pre + +:line + +[Restart, fix_modify, output, run start/stop, minimize info:] + +If the {temp yes} keyword is used in {fix drude/transform/inverse} +this fix computes a global vector with 6 components which can be +accessed by various "output commands"_Section_howto.html#howto_15. The +meaning of the components are + +temperature of the centers of mass (temperature units) +temperature of the dipoles (temperature units) +number of degrees of freedom of the centers of mass +number of degrees of freedom of the dipoles +kinetic energy of the centers of mass (energy units) +kinetic energy of the dipoles (energy units) :ol + +No information about this fix is written to "binary restart +files"_restart.html. + +[Restrictions:] none + +[Related commands:] + +"fix drude"_fix_drude.html, +"fix langevin/drude"_fix_langevin_drude.html, +"compute temp/drude"_compute_temp_drude.html, +"pair_style thole"_pair_thole.html + +[Default:] + +The option defaults are temp = no. + +:line + +:link(Lamoureux) +[(Lamoureux)] Lamoureux and Roux, J Chem Phys, 119, 3025-3039 (2003). diff --git a/doc/fix_langevin_drude.html b/doc/fix_langevin_drude.html new file mode 100644 index 0000000000..c4669040d9 --- /dev/null +++ b/doc/fix_langevin_drude.html @@ -0,0 +1,272 @@ + + + + +Syntax: +
+fix ID group-ID langevin/drude Tcom damp_com seed_com Tdrude damp_drude seed_drude keyword values ... ++
zero value = no or yes + no = do not set total random force on centers of mass to zero + yes = set total random force on centers of mass to zero ++ +
Examples: +
+fix 3 all langevin/drude 300.0 100.0 19377 1.0 20.0 83451 +fix 1 all langevin/drude 298.15 100.0 19377 5.0 10.0 83451 zero yes ++
Description: +
+Apply two Langevin thermostats as described in (Jiang) for +thermalizing the reduced degrees of freedom of Drude oscillators. +This link describes how to use the thermalized Drude oscillator +model in LAMMPS and polarizable models in LAMMPS +are discussed in this Section. +
+Drude oscillators are a way to simulate polarizables atoms, by +splitting them into a core and a Drude particle bound by a harmonic +bond. The thermalization works by transforming the particles degrees +of freedom by these equations. In these equations upper case denotes +atomic or center of mass values and lower case denotes Drude particle +or dipole values. Primes denote the transformed (reduced) values, +while bare letters denote the original values. +
+Velocities: +\begin{equation} V' = \frac {M\, V + m\, v} {M'} \end{equation} +\begin{equation} v' = v - V \end{equation} +Masses: +\begin{equation} M' = M + m \end{equation} +\begin{equation} m' = \frac {M\, m } {M'} \end{equation} +The Langevin forces are computed as +\begin{equation} F' = - \frac {M'} {\mathtt{damp\_com}}\, V' + F_r' \end{equation} +\begin{equation} f' = - \frac {m'} {\mathtt{damp\_drude}}\, v' + f_r' \end{equation} +\(F_r'\) is a random force proportional to +\(\sqrt { \frac {2\, k_B \mathtt{Tcom}\, m'} + {\mathrm dt\, \mathtt{damp\_com} } + } \). +Then the real forces acting on the particles are computed from the inverse +transform: +\begin{equation} F = \frac M {M'}\, F' - f' \end{equation} +\begin{equation} f = \frac m {M'}\, F' + f' \end{equation} +
+For Drude pairs (core + electron), the center of mass and the dipole +are thermostated if (and only if) the core atom is in the specified +group. +
+Note that the thermostat effect of this fix is applied to only the +translational degrees of freedom of the particles, which is an +important consideration if finite-size particles, which have +rotational degrees of freedom, are being thermostated. The +translational degrees of freedom can also have a bias velocity removed +from them before thermostating takes place; see the description below. +
+IMPORTANT NOTE: Like the fix langevin command, +this fix does NOT perform time integration. It only modifies forces to +effect thermostating. Thus you must use a separate time integration +fix, like fix nve or fix nph to actually +update the velocities and positions of atoms using the modified +forces. Likewise, this fix should not normally be used on atoms that +also have their temperature controlled by another fix - e.g. by fix +nvt or fix temp/rescale commands. +
+See this howto section of the manual for +a discussion of different ways to compute temperature and perform +thermostating. +
+This fix requires each atom know whether it is a Drude particle or +not. You must therefore use the fix drude command to +specify the Drude status of each atom type. +
+IMPORTANT NOTE: only the Drude core atoms need to be in the group +specified for this fix. A Drude electron will be transformed together +with its cores even if it is not itself in the group. It is safe to +include Drude electrons or non-polarizable atoms in the group. The +non-polarizable atoms will simply not be thermostatted as if they had +a massless Drude partner (electron). +
+IMPORTANT NOTE: Ghost atoms need to know their velocity for this fix +to act correctly. You must use the comm_modify +command to enable this, e.g. +
+comm_modify vel yes ++
Tcom is the target temperature of the centers of mass, which would +be used to thermostate the non-polarizable atoms. Tdrude is the +(normally low) target temperature of the core-Drude particle pairs +(dipoles). Tcom and Tdrude can be specified as an equal-style +variable. If the value is a variable, it should be +specified as v_name, where name is the variable name. In this case, +the variable will be evaluated each timestep, and its value used to +determine the target temperature. +
+Equal-style variables can specify formulas with various mathematical +functions, and include thermo_style command +keywords for the simulation box parameters and timestep and elapsed +time. Thus it is easy to specify a time-dependent temperature. +
+Like other fixes that perform thermostating, this fix can be used with +compute commands that remove a "bias" from the atom +velocities. E.g. removing the center-of-mass velocity from a group of +atoms. This is not done by default, but only if the +fix_modify command is used to assign a temperature +compute to this fix that includes such a bias term. See the doc pages +for individual compute commands to determine which ones +include a bias. In this case, the thermostat works in the following +manner: bias is removed from each atom, thermostating is performed on +the remaining thermal degrees of freedom, and the bias is added back +in. NOTE: this feature has not been tested. +
+Note: The temperature thermostating the core-Drude particle pairs +should be chosen low enough, so as to mimic as closely as possible the +self-consistent minimization. It must however be high enough, so that +the dipoles can follow the local electric field exerted by the +neighbouring atoms. The optimal value probably depends on the +temperature of the centers of mass and on the mass of the Drude +particles. +
+damp_com is the characteristic time for reaching thermal equilibrium +of the centers of mass. For example, a value of 100.0 means to relax +the temperature of the centers of mass in a timespan of (roughly) 100 +time units (tau or fmsec or psec - see the units +command). damp_drude is the characteristic time for reaching +thermal equilibrium of the dipoles. It is typically a few timesteps. +
+The number seed_com and seed_drude are positive integers. They set +the seeds of the Marsaglia random number generators used for +generating the random forces on centers of mass and on the +dipoles. Each processor uses the input seed to generate its own unique +seed and its own stream of random numbers. Thus the dynamics of the +system will not be identical on two runs on different numbers of +processors. +
+The keyword zero can be used to eliminate drift due to the +thermostat on centers of mass. Because the random forces on different +centers of mass are independent, they do not sum exactly to zero. As +a result, this fix applies a small random force to the entire system, +and the momentum of the total center of mass of the system undergoes a +slow random walk. If the keyword zero is set to yes, the total +random force on the centers of mass is set exactly to zero by +subtracting off an equal part of it from each center of mass in the +group. As a result, the total center of mass of a system with zero +initial momentum will not drift over time. +
+Example for rigid bodies in the NPT ensemble: +
+comm_modify vel yes +fix TEMP all langevin/drude 300. 100. 1256 1. 20. 13977 zero yes +fix NPH ATOMS rigid/nph/small molecule iso 1. 1. 500. +fix NVE DRUDES nve +thermo_style custom step cpu etotal ke pe ebond ecoul elong press vol temp c_TATOM c_TDRUDE f_TEMP[1] f_TEMP[2] ++
Comments: +
+Restart, fix_modify, output, run start/stop, minimize info: +
+No information about this fix is written to binary restart +files. Because the state of the random number generator +is not saved in restart files, this means you cannot do "exact" +restarts with this fix, where the simulation continues on the same as +if no restart had taken place. However, in a statistical sense, a +restarted simulation should produce the same behavior. +
+The fix_modify temp option is supported by this +fix. You can use it to assign a temperature compute +you have defined to this fix which will be used in its thermostating +procedure, as described above. For consistency, the group used by the +compute should include the group of this fix and the Drude particles. +
+This fix computes a global vector with 6 components which can be +accessed by various output commands. The +meaning of the components are as follows: +
+This fix is not invoked during energy minimization. +
+Restrictions: none +
+Related commands: +
+fix langevin, +fix drude/transform, +compute temp/drude, +pair_style thole +
+Default: +
+The option defaults are zero = no. +
+(Jiang) Jiang, Hardy, Phillips, MacKerell, Schulten, and Roux, J +Phys Chem Lett, 2, 87-92 (2011). +
+ diff --git a/doc/fix_langevin_drude.txt b/doc/fix_langevin_drude.txt new file mode 100644 index 0000000000..f936490b95 --- /dev/null +++ b/doc/fix_langevin_drude.txt @@ -0,0 +1,253 @@ + + + +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +fix langevin/drude command :h3 + +[Syntax:] + +fix ID group-ID langevin/drude Tcom damp_com seed_com Tdrude damp_drude seed_drude keyword values ... :pre + +ID, group-ID are documented in "fix"_fix.html command :ulb,l +langevin/drude = style name of this fix command :l +Tcom = desired temperature of the centers of mass (temperature units) :l +damp_com = damping parameter for the thermostat on centers of mass (time units) :l +seed_com = random number seed to use for white noise of the thermostat on centers of mass (positive integer) :l +Tdrude = desired temperature of the Drude oscillators (temperature units) :l +damp_drude = damping parameter for the thermostat on Drude oscillators (time units) :l +seed_drude = random number seed to use for white noise of the thermostat on Drude oscillators (positive integer) :l +zero or more keyword/value pairs may be appended :l +keyword = {zero} :l + {zero} value = {no} or {yes} + {no} = do not set total random force on centers of mass to zero + {yes} = set total random force on centers of mass to zero :pre +:ule + +[Examples:] + +fix 3 all langevin/drude 300.0 100.0 19377 1.0 20.0 83451 +fix 1 all langevin/drude 298.15 100.0 19377 5.0 10.0 83451 zero yes :pre + +[Description:] + +Apply two Langevin thermostats as described in "(Jiang)"_#Jiang for +thermalizing the reduced degrees of freedom of Drude oscillators. +This link describes how to use the "thermalized Drude oscillator +model"_tutorial_drude.html in LAMMPS and polarizable models in LAMMPS +are discussed in "this Section"_Section_howto.html#howto_25. + +Drude oscillators are a way to simulate polarizables atoms, by +splitting them into a core and a Drude particle bound by a harmonic +bond. The thermalization works by transforming the particles degrees +of freedom by these equations. In these equations upper case denotes +atomic or center of mass values and lower case denotes Drude particle +or dipole values. Primes denote the transformed (reduced) values, +while bare letters denote the original values. + +Velocities: +\begin\{equation\} V' = \frac \{M\, V + m\, v\} \{M'\} \end\{equation\} +\begin\{equation\} v' = v - V \end\{equation\} +Masses: +\begin\{equation\} M' = M + m \end\{equation\} +\begin\{equation\} m' = \frac \{M\, m \} \{M'\} \end\{equation\} +The Langevin forces are computed as +\begin\{equation\} F' = - \frac \{M'\} \{\mathtt\{damp\_com\}\}\, V' + F_r' \end\{equation\} +\begin\{equation\} f' = - \frac \{m'\} \{\mathtt\{damp\_drude\}\}\, v' + f_r' \end\{equation\} +\(F_r'\) is a random force proportional to +\(\sqrt \{ \frac \{2\, k_B \mathtt\{Tcom\}\, m'\} + \{\mathrm dt\, \mathtt\{damp\_com\} \} + \} \). :b +\(f_r'\) is a random force proportional to +\(\sqrt \{ \frac \{2\, k_B \mathtt\{Tdrude\}\, m'\} + \{\mathrm dt\, \mathtt\{damp\_drude\} \} + \} \). :b +Then the real forces acting on the particles are computed from the inverse +transform: +\begin\{equation\} F = \frac M \{M'\}\, F' - f' \end\{equation\} +\begin\{equation\} f = \frac m \{M'\}\, F' + f' \end\{equation\} + +For Drude pairs (core + electron), the center of mass and the dipole +are thermostated if (and only if) the core atom is in the specified +group. + +Note that the thermostat effect of this fix is applied to only the +translational degrees of freedom of the particles, which is an +important consideration if finite-size particles, which have +rotational degrees of freedom, are being thermostated. The +translational degrees of freedom can also have a bias velocity removed +from them before thermostating takes place; see the description below. + +IMPORTANT NOTE: Like the "fix langevin"_fix_langevin.html command, +this fix does NOT perform time integration. It only modifies forces to +effect thermostating. Thus you must use a separate time integration +fix, like "fix nve"_fix_nve.html or "fix nph"_fix_nh.html to actually +update the velocities and positions of atoms using the modified +forces. Likewise, this fix should not normally be used on atoms that +also have their temperature controlled by another fix - e.g. by "fix +nvt"_fix_nh.html or "fix temp/rescale"_fix_temp_rescale.html commands. + +See "this howto section"_Section_howto.html#howto_16 of the manual for +a discussion of different ways to compute temperature and perform +thermostating. + +:line + +This fix requires each atom know whether it is a Drude particle or +not. You must therefore use the "fix drude"_fix_drude.html command to +specify the Drude status of each atom type. + +IMPORTANT NOTE: only the Drude core atoms need to be in the group +specified for this fix. A Drude electron will be transformed together +with its cores even if it is not itself in the group. It is safe to +include Drude electrons or non-polarizable atoms in the group. The +non-polarizable atoms will simply not be thermostatted as if they had +a massless Drude partner (electron). + +IMPORTANT NOTE: Ghost atoms need to know their velocity for this fix +to act correctly. You must use the "comm_modify"_comm_modify.html +command to enable this, e.g. + +comm_modify vel yes :pre + +:line + +{Tcom} is the target temperature of the centers of mass, which would +be used to thermostate the non-polarizable atoms. {Tdrude} is the +(normally low) target temperature of the core-Drude particle pairs +(dipoles). {Tcom} and {Tdrude} can be specified as an equal-style +"variable"_variable.html. If the value is a variable, it should be +specified as v_name, where name is the variable name. In this case, +the variable will be evaluated each timestep, and its value used to +determine the target temperature. + +Equal-style variables can specify formulas with various mathematical +functions, and include "thermo_style"_thermo_style.html command +keywords for the simulation box parameters and timestep and elapsed +time. Thus it is easy to specify a time-dependent temperature. + +Like other fixes that perform thermostating, this fix can be used with +"compute commands"_compute.html that remove a "bias" from the atom +velocities. E.g. removing the center-of-mass velocity from a group of +atoms. This is not done by default, but only if the +"fix_modify"_fix_modify.html command is used to assign a temperature +compute to this fix that includes such a bias term. See the doc pages +for individual "compute commands"_compute.html to determine which ones +include a bias. In this case, the thermostat works in the following +manner: bias is removed from each atom, thermostating is performed on +the remaining thermal degrees of freedom, and the bias is added back +in. NOTE: this feature has not been tested. + +Note: The temperature thermostating the core-Drude particle pairs +should be chosen low enough, so as to mimic as closely as possible the +self-consistent minimization. It must however be high enough, so that +the dipoles can follow the local electric field exerted by the +neighbouring atoms. The optimal value probably depends on the +temperature of the centers of mass and on the mass of the Drude +particles. + +{damp_com} is the characteristic time for reaching thermal equilibrium +of the centers of mass. For example, a value of 100.0 means to relax +the temperature of the centers of mass in a timespan of (roughly) 100 +time units (tau or fmsec or psec - see the "units"_units.html +command). {damp_drude} is the characteristic time for reaching +thermal equilibrium of the dipoles. It is typically a few timesteps. + +The number {seed_com} and {seed_drude} are positive integers. They set +the seeds of the Marsaglia random number generators used for +generating the random forces on centers of mass and on the +dipoles. Each processor uses the input seed to generate its own unique +seed and its own stream of random numbers. Thus the dynamics of the +system will not be identical on two runs on different numbers of +processors. + +The keyword {zero} can be used to eliminate drift due to the +thermostat on centers of mass. Because the random forces on different +centers of mass are independent, they do not sum exactly to zero. As +a result, this fix applies a small random force to the entire system, +and the momentum of the total center of mass of the system undergoes a +slow random walk. If the keyword {zero} is set to {yes}, the total +random force on the centers of mass is set exactly to zero by +subtracting off an equal part of it from each center of mass in the +group. As a result, the total center of mass of a system with zero +initial momentum will not drift over time. + +:line + +Example for rigid bodies in the NPT ensemble: + +comm_modify vel yes +fix TEMP all langevin/drude 300. 100. 1256 1. 20. 13977 zero yes +fix NPH ATOMS rigid/nph/small molecule iso 1. 1. 500. +fix NVE DRUDES nve +thermo_style custom step cpu etotal ke pe ebond ecoul elong press vol temp c_TATOM c_TDRUDE f_TEMP\[1\] f_TEMP\[2\] :pre + +Comments: + +Drude particles should not be in the rigid group, otherwise the Drude oscillators will be frozen and the system will lose its polarizability. +{zero yes} avoids a drift of the center of mass of the system, but is a bit slower. +use two different random seeds to avoid unphysical correlations. +temperature is controlled by the fix {langevin/drude}, so the time-integration fixes do not thermostate. +don't forget to time-integrate both cores and Drude particles. +pressure is time-integrated only once by using {nve} for Drude particles and {nph} for atoms/cores (or vice versa). Do not use {nph} for both. +contrary to the alternative thermostating using Nose-Hoover thermostat fix {npt} and fix {drude/transform}, the {fix_modify} command is not required here, because the fix {nph} computes the global pressure even if its group is {ATOMS}. This is what we want. If we thermostated {ATOMS} using {npt}, the pressure should be the global one, but the temperature should be only that of the cores. That's why the command {fix_modify} should be called in that case. +f_TEMP\[1\] and f_TEMP\[2\] contain the reduced temperatures of the cores/atoms and of the Drude particles (see below). They should be 300. and 1. on average here. :ul + +:line + +[Restart, fix_modify, output, run start/stop, minimize info:] + +No information about this fix is written to "binary restart +files"_restart.html. Because the state of the random number generator +is not saved in restart files, this means you cannot do "exact" +restarts with this fix, where the simulation continues on the same as +if no restart had taken place. However, in a statistical sense, a +restarted simulation should produce the same behavior. + +The "fix_modify"_fix_modify.html {temp} option is supported by this +fix. You can use it to assign a temperature "compute"_compute.html +you have defined to this fix which will be used in its thermostating +procedure, as described above. For consistency, the group used by the +compute should include the group of this fix and the Drude particles. + +This fix computes a global vector with 6 components which can be +accessed by various "output commands"_Section_howto.html#howto_15. The +meaning of the components are as follows: + +temperature of the centers of mass (temperature units) +temperature of the dipoles (temperature units) +number of degrees of freedom of the centers of mass +number of degrees of freedom of the dipoles +kinetic energy of the centers of mass (energy units) +kinetic energy of the dipoles (energy units) :ol + +This fix is not invoked during "energy minimization"_minimize.html. + +[Restrictions:] none + +[Related commands:] + +"fix langevin"_fix_langevin.html, +"fix drude/transform"_fix_drude_transform.html, +"compute temp/drude"_compute_temp_drude.html, +"pair_style thole"_pair_thole.html + +[Default:] + +The option defaults are zero = no. + +:line + +:link(Jiang) +[(Jiang)] Jiang, Hardy, Phillips, MacKerell, Schulten, and Roux, J +Phys Chem Lett, 2, 87-92 (2011). diff --git a/doc/pair_bop.html b/doc/pair_bop.html index a8ed6a2f4d..73479a3916 100644 --- a/doc/pair_bop.html +++ b/doc/pair_bop.html @@ -40,7 +40,9 @@ transferability to different phases can approach that of quantum mechanical methods. This potential is similar to the original BOP developed by Pettifor (Pettifor_1, Pettifor_2, Pettifor_3) and later updated -by Murdick, Zhou, and Ward (Murdick, Ward). +by Murdick, Zhou, and Ward (Murdick, Ward). As of +summer 2015, BOP potential files for these systems are provide with +LAMMPS: AlCu, CCu, CdTe, CuH, GaAs.The BOP potential consists of three terms:
diff --git a/doc/pair_bop.txt b/doc/pair_bop.txt index d1c7e77392..9c071aaa9a 100644 --- a/doc/pair_bop.txt +++ b/doc/pair_bop.txt @@ -34,7 +34,9 @@ transferability to different phases can approach that of quantum mechanical methods. This potential is similar to the original BOP developed by Pettifor ("Pettifor_1"_#Pettifor_1, "Pettifor_2"_#Pettifor_2, "Pettifor_3"_#Pettifor_3) and later updated -by Murdick, Zhou, and Ward ("Murdick"_#Murdick, "Ward"_#Ward). +by Murdick, Zhou, and Ward ("Murdick"_#Murdick, "Ward"_#Ward). As of +summer 2015, BOP potential files for these systems are provide with +LAMMPS: AlCu, CCu, CdTe, CuH, GaAs. The BOP potential consists of three terms: diff --git a/doc/pair_born.html b/doc/pair_born.html index f0b34ed097..bd2ba588fd 100644 --- a/doc/pair_born.html +++ b/doc/pair_born.html @@ -17,6 +17,8 @@pair_style style args-
born args = cutoff
cutoff = global cutoff for non-Coulombic interactions (distance units)
- born/coul/long args = cutoff (cutoff2)
+ born/coul/long or born/coul/long/cs args = cutoff (cutoff2)
cutoff = global cutoff for non-Coulombic (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units)
born/coul/msm args = cutoff (cutoff2)
@@ -60,7 +62,9 @@ pair_coeff * * 6.08 0.317 2.340 24.18 11.51
pair_coeff 1 1 6.08 0.317 2.340 24.18 11.51
pair_style born/coul/long 10.0 +pair_style born/coul/long/cs 10.0 pair_style born/coul/long 10.0 8.0 +pair_style born/coul/long/cs 10.0 8.0 pair_coeff * * 6.08 0.317 2.340 24.18 11.51 pair_coeff 1 1 6.08 0.317 2.340 24.18 11.51@@ -102,7 +106,12 @@ term.
The born/coul/wolf style adds a Coulombic term as described for the Wolf potential in the coul/wolf pair style.
-Note that these potentials are related to the Buckingham
+ Style born/coul/long/cs is identical to born/coul/long except that
+a term is added for the core/shell model
+to allow charges on core and shell particles to be separated by r =
+0.0.
+ Note that these potentials are related to the Buckingham
potential.
The following coefficients must be defined for each pair of atoms
diff --git a/doc/pair_born.txt b/doc/pair_born.txt
index 676b87882a..8da606f9da 100644
--- a/doc/pair_born.txt
+++ b/doc/pair_born.txt
@@ -11,6 +11,7 @@ pair_style born command :h3
pair_style born/omp command :h3
pair_style born/gpu command :h3
pair_style born/coul/long command :h3
+pair_style born/coul/long/cs command :h3
pair_style born/coul/long/cuda command :h3
pair_style born/coul/long/gpu command :h3
pair_style born/coul/long/omp command :h3
@@ -24,11 +25,11 @@ pair_style born/coul/wolf/omp command :h3
pair_style style args :pre
-style = {born} or {born/coul/long} or {born/coul/msm} or {born/coul/wolf}
+style = {born} or {born/coul/long} or {born/coul/long/cs} or {born/coul/msm} or {born/coul/wolf}
args = list of arguments for a particular style :ul
{born} args = cutoff
cutoff = global cutoff for non-Coulombic interactions (distance units)
- {born/coul/long} args = cutoff (cutoff2)
+ {born/coul/long} or {born/coul/long/cs} args = cutoff (cutoff2)
cutoff = global cutoff for non-Coulombic (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units)
{born/coul/msm} args = cutoff (cutoff2)
@@ -46,7 +47,9 @@ pair_coeff * * 6.08 0.317 2.340 24.18 11.51
pair_coeff 1 1 6.08 0.317 2.340 24.18 11.51 :pre
pair_style born/coul/long 10.0
+pair_style born/coul/long/cs 10.0
pair_style born/coul/long 10.0 8.0
+pair_style born/coul/long/cs 10.0 8.0
pair_coeff * * 6.08 0.317 2.340 24.18 11.51
pair_coeff 1 1 6.08 0.317 2.340 24.18 11.51 :pre
@@ -88,8 +91,13 @@ term.
The {born/coul/wolf} style adds a Coulombic term as described for the
Wolf potential in the "coul/wolf"_pair_coul.html pair style.
+Style {born/coul/long/cs} is identical to {born/coul/long} except that
+a term is added for the "core/shell model"_Section_howto.html#howto_25
+to allow charges on core and shell particles to be separated by r =
+0.0.
+
Note that these potentials are related to the "Buckingham
-potential"_pair_born.html.
+potential"_pair_buck.html.
The following coefficients must be defined for each pair of atoms
types via the "pair_coeff"_pair_coeff.html command as in the examples
diff --git a/doc/pair_buck.html b/doc/pair_buck.html
index 1e14667251..046d2f7121 100644
--- a/doc/pair_buck.html
+++ b/doc/pair_buck.html
@@ -31,6 +31,8 @@
pair_style style args-
buck args = cutoff
@@ -55,7 +57,7 @@
buck/coul/cut args = cutoff (cutoff2)
cutoff = global cutoff for Buckingham (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units)
- buck/coul/long args = cutoff (cutoff2)
+ buck/coul/long or buck/coul/long/cs args = cutoff (cutoff2)
cutoff = global cutoff for Buckingham (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units)
buck/coul/msm args = cutoff (cutoff2)
@@ -75,7 +77,9 @@ pair_coeff 1 1 100.0 1.5 200.0 9.0
pair_coeff 1 1 100.0 1.5 200.0 9.0 8.0
pair_style buck/coul/long 10.0 +pair_style buck/coul/long/cs 10.0 pair_style buck/coul/long 10.0 8.0 +pair_style buck/coul/long/cs 10.0 8.0 pair_coeff * * 100.0 1.5 200.0 pair_coeff 1 1 100.0 1.5 200.0 9.0@@ -109,6 +113,11 @@ A,C and Coulombic terms. If two cutoffs are specified, the first is used as the cutoff for the A,C terms, and the second is the cutoff for the Coulombic term. +
Style buck/coul/long/cs is identical to buck/coul/long except that +a term is added for the core/shell model +to allow charges on core and shell particles to be separated by r = +0.0. +
Note that these potentials are related to the Born-Mayer-Huggins potential.
@@ -179,7 +188,7 @@ I,J pairs must be specified explicitly. for the energy of the exp() and 1/r^6 portion of the pair interaction.The buck/coul/long pair style supports the -pair_modify table option ti tabulate the +pair_modify table option to tabulate the short-range portion of the long-range Coulombic interaction.
These styles support the pair_modify tail option for adding long-range @@ -196,8 +205,9 @@ respa command. They do not support the inner,
Restrictions:
-The buck/coul/long style is part of the KSPACE package. It is only -enabled if LAMMPS was built with that package (which it is by +
The buck/coul/long style is part of the KSPACE package. The +buck/coul/long/cs style is part of the CORESHELL package. They are +only enabled if LAMMPS was built with that package (which it is by default). See the Making LAMMPS section for more info.
diff --git a/doc/pair_buck.txt b/doc/pair_buck.txt index 2409c24ee9..2452f10e07 100644 --- a/doc/pair_buck.txt +++ b/doc/pair_buck.txt @@ -17,6 +17,7 @@ pair_style buck/coul/cut/gpu command :h3 pair_style buck/coul/cut/kk command :h3 pair_style buck/coul/cut/omp command :h3 pair_style buck/coul/long command :h3 +pair_style buck/coul/long/cs command :h3 pair_style buck/coul/long/cuda command :h3 pair_style buck/coul/long/gpu command :h3 pair_style buck/coul/long/kk command :h3 @@ -28,14 +29,14 @@ pair_style buck/coul/msm/omp command :h3 pair_style style args :pre -style = {buck} or {buck/coul/cut} or {buck/coul/long} or {buck/coul/msm} +style = {buck} or {buck/coul/cut} or {buck/coul/long} or {buck/coul/long/cs} or {buck/coul/msm} args = list of arguments for a particular style :ul {buck} args = cutoff cutoff = global cutoff for Buckingham interactions (distance units) {buck/coul/cut} args = cutoff (cutoff2) cutoff = global cutoff for Buckingham (and Coulombic if only 1 arg) (distance units) cutoff2 = global cutoff for Coulombic (optional) (distance units) - {buck/coul/long} args = cutoff (cutoff2) + {buck/coul/long} or {buck/coul/long/cs} args = cutoff (cutoff2) cutoff = global cutoff for Buckingham (and Coulombic if only 1 arg) (distance units) cutoff2 = global cutoff for Coulombic (optional) (distance units) {buck/coul/msm} args = cutoff (cutoff2) @@ -55,7 +56,9 @@ pair_coeff 1 1 100.0 1.5 200.0 9.0 pair_coeff 1 1 100.0 1.5 200.0 9.0 8.0 :pre pair_style buck/coul/long 10.0 +pair_style buck/coul/long/cs 10.0 pair_style buck/coul/long 10.0 8.0 +pair_style buck/coul/long/cs 10.0 8.0 pair_coeff * * 100.0 1.5 200.0 pair_coeff 1 1 100.0 1.5 200.0 9.0 :pre @@ -89,6 +92,11 @@ A,C and Coulombic terms. If two cutoffs are specified, the first is used as the cutoff for the A,C terms, and the second is the cutoff for the Coulombic term. +Style {buck/coul/long/cs} is identical to {buck/coul/long} except that +a term is added for the "core/shell model"_Section_howto.html#howto_25 +to allow charges on core and shell particles to be separated by r = +0.0. + Note that these potentials are related to the "Born-Mayer-Huggins potential"_pair_born.html. @@ -159,7 +167,7 @@ These styles support the "pair_modify"_pair_modify.html shift option for the energy of the exp() and 1/r^6 portion of the pair interaction. The {buck/coul/long} pair style supports the -"pair_modify"_pair_modify.html table option ti tabulate the +"pair_modify"_pair_modify.html table option to tabulate the short-range portion of the long-range Coulombic interaction. These styles support the pair_modify tail option for adding long-range @@ -176,8 +184,9 @@ respa"_run_style.html command. They do not support the {inner}, [Restrictions:] -The {buck/coul/long} style is part of the KSPACE package. It is only -enabled if LAMMPS was built with that package (which it is by +The {buck/coul/long} style is part of the KSPACE package. The +{buck/coul/long/cs} style is part of the CORESHELL package. They are +only enabled if LAMMPS was built with that package (which it is by default). See the "Making LAMMPS"_Section_start.html#start_3 section for more info. diff --git a/doc/pair_coul.html b/doc/pair_coul.html index 176e331f89..db95f8848a 100644 --- a/doc/pair_coul.html +++ b/doc/pair_coul.html @@ -35,6 +35,8 @@pair_style coul/long 10.0 +pair_style coul/long/cs 10.0 pair_coeff * *
pair_style coul/msm 10.0 @@ -225,6 +229,10 @@ option. The Coulombic cutoff specified for this style means that pairwise interactions within this distance are computed directly; interactions outside that distance are computed in reciprocal space. +Style coul/long/cs is identical to coul/long except that a term is +added for the core/shell model to allow +charges on core and shell particles to be separated by r = 0.0. +
Styles tip4p/cut and tip4p/long implement the coulomb part of the TIP4P water model of (Jorgensen), which introduces a massless site located a short distance away from the oxygen atom @@ -333,16 +341,15 @@ to be specified in an input script that reads a restart file.
Restrictions:
The coul/long, coul/msm and tip4p/long styles are part of the -KSPACE package. They are only enabled if LAMMPS was built with that -package (which it is by default). - See the Making LAMMPS section -for more info. +KSPACE package. The coul/long/cs style is part of the CORESHELL +package. They are only enabled if LAMMPS was built with that package +(which it is by default). See the Making +LAMMPS section for more info.
Related commands:
-pair_coeff, pair_style -hybrid/overlay -kspace_style +
pair_coeff, pair_style, +hybrid/overlay, kspace_style
Default: none
diff --git a/doc/pair_coul.txt b/doc/pair_coul.txt index a6fbda3785..746878fe4d 100644 --- a/doc/pair_coul.txt +++ b/doc/pair_coul.txt @@ -19,6 +19,7 @@ pair_style coul/dsf/gpu command :h3 pair_style coul/dsf/kk command :h3 pair_style coul/dsf/omp command :h3 pair_style coul/long command :h3 +pair_style coul/long/cs command :h3 pair_style coul/long/omp command :h3 pair_style coul/long/gpu command :h3 pair_style coul/long/kk command :h3 @@ -39,6 +40,7 @@ pair_style coul/cut cutoff pair_style coul/debye kappa cutoff pair_style coul/dsf alpha cutoff pair_style coul/long cutoff +pair_style coul/long/cs cutoff pair_style coul/long/gpu cutoff pair_style coul/wolf alpha cutoff pair_style coul/streitz cutoff keyword alpha @@ -63,6 +65,7 @@ pair_style coul/dsf 0.05 10.0 pair_coeff * * :pre pair_style coul/long 10.0 +pair_style coul/long/cs 10.0 pair_coeff * * :pre pair_style coul/msm 10.0 @@ -197,6 +200,10 @@ option. The Coulombic cutoff specified for this style means that pairwise interactions within this distance are computed directly; interactions outside that distance are computed in reciprocal space. +Style {coul/long/cs} is identical to {coul/long} except that a term is +added for the "core/shell model"_Section_howto.html#howto_25 to allow +charges on core and shell particles to be separated by r = 0.0. + Styles {tip4p/cut} and {tip4p/long} implement the coulomb part of the TIP4P water model of "(Jorgensen)"_#Jorgensen, which introduces a massless site located a short distance away from the oxygen atom @@ -305,16 +312,15 @@ This pair style can only be used via the {pair} keyword of the [Restrictions:] The {coul/long}, {coul/msm} and {tip4p/long} styles are part of the -KSPACE package. They are only enabled if LAMMPS was built with that -package (which it is by default). - See the "Making LAMMPS"_Section_start.html#start_3 section -for more info. +KSPACE package. The {coul/long/cs} style is part of the CORESHELL +package. They are only enabled if LAMMPS was built with that package +(which it is by default). See the "Making +LAMMPS"_Section_start.html#start_3 section for more info. [Related commands:] -"pair_coeff"_pair_coeff.html, "pair_style -hybrid/overlay"_pair_hybrid.html -"kspace_style"_kspace_style.html +"pair_coeff"_pair_coeff.html, "pair_style, +hybrid/overlay"_pair_hybrid.html, "kspace_style"_kspace_style.html [Default:] none diff --git a/doc/tutorial_drude.html b/doc/tutorial_drude.html new file mode 100644 index 0000000000..62509b0f90 --- /dev/null +++ b/doc/tutorial_drude.html @@ -0,0 +1,464 @@ + + + + +Tutorial for Thermalized Drude oscillators in LAMMPS +
+This tutorial explains how to use Drude +oscillators in LAMMPS to simulate polarizable +systems. As an illustration, the input files for a simulation of 250 +phenol molecules are documented. First of all, LAMMPS has to be +compiled with the USER-DRUDE package activated. Then, the data file +and input scripts have to be modified to include the Drude dipoles and +how to handle them. +
+
+ +Overview of Drude induced dipoles +
+Polarizable atoms acquire an induced electric dipole moment under the +action of an external electric field, for example the electric field +created by the surrounding particles. Drude oscillators represent +these dipoles by two fixed charges: the core (DC) and the Drude +particle (DP) bound by a harmonic potential. The Drude particle can be +thought of as the electron cloud whose center can be displaced from +the position of the the corresponding nucleus. +
+The sum of the masses of a core-Drude pair should be the mass of the +initial (unsplit) atom, \(m_C + m_D = m\). The sum of their charges +should be the charge of the initial (unsplit) atom, \(q_C + q_D = q\). +A harmonic potential between the core and Drude partners should be +present, with force constant \(k_D\) and an equilibrium distance of +zero. The (half-)stiffness of the harmonic bond +\(K_D = k_D/2\) and the Drude charge \(q_D\) are related to the atom +polarizability \(\alpha\) by +
+\begin{equation} K_D = \frac 1 2\, \frac {q_D^2} \alpha +\end{equation} +
+Ideally, the mass of the Drude particle should be small, and the +stiffness of the harmonic bond should be large, so that the Drude +particle remains close ot the core. 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. +
+
The motion of of the Drude particles can be calculated by minimizing +the energy of the induced dipoles at each timestep, by an interative, +self-consistent procedure. The Drude particles can be massless and +therefore do not contribute to the kinetic energy. However, the +relaxed method is computationall slow. An extended-lagrangian method +can be used to calculate the positions of the Drude particles, but +this requires them to have mass. It is important in this case to +decouple the degrees of freedom associated with the Drude oscillators +from those of the normal atoms. Thermalizing the Drude dipoles at +temperatures comparable to the rest of the simulation leads to several +problems (kinetic energy transfer, very short timestep, etc.), which +can be remediated by the "cold Drude" technique (Lamoureux and +Roux). +
+Two closely related models are used to represent polarization through +"charges on a spring": the core-shell model and the Drude +model. Although the basic idea is the same, the core-shell model is +normally used for ionic/crystalline materials, whereas the Drude model +is normally used for molecular systems and fluid states. In ionic +crystals the symmetry around each ion and the distance between them +are such that the core-shell model is sufficiently stable. But to be +applicable to molecular/covalent systems the Drude model includes two +important features: +
+Preparation of the data file +
+The data file is similar to a standard LAMMPS data file for +atom_style full. The DPs and the harmonic bonds connecting them +to their DC should appear in the data file as normal atoms and bonds. +
+You can use the polarizer tool (Python script distributed with the +USER-DRUDE package) to convert a non-polarizable data file (here +data.102494.lmp) to a polarizable data file (data-p.lmp) +
+polarizer -q -f phenol.dff data.102494.lmp data-p.lmp ++
This will automatically insert the new atoms and bonds. +The masses and charges of DCs and DPs are computed +from phenol.dff, as well as the DC-DP bond constants. The file +phenol.dff contains the polarizabilities of the atom types +and the mass of the Drude particles, for instance: +
+# units: kJ/mol, A, deg +# kforce is in the form k/2 r_D^2 +# type m_D/u q_D/e k_D alpha/A3 thole +OH 0.4 -1.0 4184.0 0.63 0.67 +CA 0.4 -1.0 4184.0 1.36 2.51 +CAI 0.4 -1.0 4184.0 1.09 2.51 ++
The hydrogen atoms are absent from this file, so they will be treated +as non-polarizable atoms. In the non-polarizable data file +data.102494.lmp, atom names corresponding to the atom type numbers +have to be specified as comments at the end of lines of the Masses +section. You probably need to edit it to add these names. It should +look like +
+Masses ++
1 12.011 # CAI +2 12.011 # CA +3 15.999 # OH +4 1.008 # HA +5 1.008 # HO ++
Basic input file +
+The atom style should be set to (or derive from) full, so that you +can define atomic charges and molecular bonds, angles, dihedrals... +
+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 recognize that you are using Drude oscillators, you should +use the fix drude. The command is +
+fix DRUDE all drude 1 1 1 0 0 2 2 2 ++
or, equivalently +
+fix DRUDE all drude C C C N N D D D ++
The 0, 1, 2 (or N, C, D) following the drude keyword have the +following meaning: There is one tag for each atom type. This tag is 1 +(or C) for DCs, 2 (or D) for DPs and 0 (or N) for non-polarizable +atoms. Here the atom types 1 to 3 (C and O atoms) are DC, atom types +4 and 5 (H atoms) are non-polarizable and the atom types 6 to 8 are +the newly created DPs. +
+By recognizing the fix drude, LAMMPS will find and store matching +DC-DP pairs and will treat DP as equivalent to their DC in the +special bonds relations. It may be necessary to extend the space +for storing such special relations. In this case extra space should +be reserved by using the extra keyword of the special_bonds +command. With our phenol, there is 1 more special neighbor for which +space is required. Otherwise LAMMPS crashes and gives the required +value. +
+special_bonds lj/coul 0.0 0.0 0.5 extra 1 ++
Let us assume we want to run a simple NVT simulation at 300 K. Note +that Drude oscillators need to be thermalized at a low temperature in +order to approximate a self-consistent field (SCF), therefore it is not +possible to simulate an NVE ensemble with this package. Since dipoles +are approximated by a charged DC-DP pair, the pair_style must +include Coulomb interactions, for instance lj/cut/coul/long with +kspace_style pppm. For example, with a cutoff of 10. and a precision +1.e-4: +
+pair_style lj/cut/coul/long 10.0 +kspace_style pppm 1.0e-4 ++
As compared to the non-polarizable input file, pair_coeff lines need +to be added for the DPs. Since the DPs have no Lennard-Jones +interactions, their epsilon is 0. so the only pair_coeff line +that needs to be added is +
+pair_coeff * 6* 0.0 0.0 # All-DPs ++
Now for the thermalization, the simplest choice is to use the fix +langevin/drude. +
+fix LANG all langevin/drude 300. 100 12435 1. 20 13977 ++
This applies a Langevin thermostat at temperature 300. to the centers +of mass of the DC-DP pairs, with relaxation time 100 and with random +seed 12345. This fix applies also a Langevin thermostat at temperature +1. to the relative motion of the DPs around their DCs, with relaxation +time 20 and random seed 13977. Only the DCs and non-polarizable +atoms need to be in this fix's group. LAMMPS will thermostate the DPs +together with their DC. For this, ghost atoms need to know their +velocities. Thus you need to add the following command: +
+comm_modify vel yes ++
In order to avoid that the center of mass of the whole system +drifts due to the random forces of the Langevin thermostat on DCs, you +can add the zero yes option at the end of the fix line. +
+If the fix shake is used to constrain the C-H bonds, it should be +invoked after the fix langevin/drude for more accuracy. +
+fix SHAKE ATOMS shake 0.0001 20 0 t 4 5 ++
IMPORTANT NOTE: The group of the fix shake must not include the DPs. +If the group ATOMS is defined by non-DPs atom types, you could use +
+Since the fix langevin/drude does not perform time integration (just +modification of forces but no position/velocity updates), the fix +nve should be used in conjunction. +
+fix NVE all nve ++
Finally, do not forget to update the atom type elements if you use +them in a dump_modify ... element ... command, by adding the element +type of the DPs. Here for instance +
+dump DUMP all custom 10 dump.lammpstrj id mol type element x y z ix iy iz +dump_modify DUMP element C C O H H D D D ++
The input file should now be ready for use! +
+You will notice that the global temperature thermo_temp computed by +LAMMPS is not 300. K as wanted. This is because LAMMPS treats DPs as +standard atoms in his default compute. If you want to output the +temperatures of the DC-DP pair centers of mass and of the DPs relative +to their DCs, you should use thermo_style custom with respectively +f_LANG[1] and f_LANG[2]. These should be close to 300. and +1. on average. +
+thermo_style custom step temp f_LANG[1] f_LANG[2] ++
Thole screening +
+Dipolar interactions represented by point charges on springs may not +be stable, for example if the atomic polarizability is too high for +instance, 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 exagerated. Often, +special bond relations prevent bonded neighboring atoms to see the +charge of each other's DP, so that the problem does not always appear. +It is possible to use screened dipole dipole interactions by using 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. In our example, we would use +
+pair_style hybrid/overlay lj/cut/coul/long 10.0 thole 2.6 10.0 ++
This tells LAMMPS that we are using two pair_styles. The first one is +as above (lj/cut/coul/long 10.0). The second one is a thole +pair_style with default screening factor 2.6 (Noskov) and +cutoff 10.0. +
+Since hybrid/overlay does not support mixing rules, 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 pair_coeff section of the input file. In our +example, this will look like: +
+pair_coeff 1 1 lj/cut/coul/long 0.0700 3.550 +pair_coeff 1 2 lj/cut/coul/long 0.0700 3.550 +pair_coeff 1 3 lj/cut/coul/long 0.1091 3.310 +pair_coeff 1 4 lj/cut/coul/long 0.0458 2.985 +pair_coeff 2 2 lj/cut/coul/long 0.0700 3.550 +pair_coeff 2 3 lj/cut/coul/long 0.1091 3.310 +pair_coeff 2 4 lj/cut/coul/long 0.0458 2.985 +pair_coeff 3 3 lj/cut/coul/long 0.1700 3.070 +pair_coeff 3 4 lj/cut/coul/long 0.0714 2.745 +pair_coeff 4 4 lj/cut/coul/long 0.0300 2.420 +pair_coeff * 5 lj/cut/coul/long 0.0000 0.000 +pair_coeff * 6* lj/cut/coul/long 0.0000 0.000 +pair_coeff 1 1 thole 1.090 2.510 +pair_coeff 1 2 thole 1.218 2.510 +pair_coeff 1 3 thole 0.829 1.590 +pair_coeff 1 6 thole 1.090 2.510 +pair_coeff 1 7 thole 1.218 2.510 +pair_coeff 1 8 thole 0.829 1.590 +pair_coeff 2 2 thole 1.360 2.510 +pair_coeff 2 3 thole 0.926 1.590 +pair_coeff 2 6 thole 1.218 2.510 +pair_coeff 2 7 thole 1.360 2.510 +pair_coeff 2 8 thole 0.926 1.590 +pair_coeff 3 3 thole 0.630 0.670 +pair_coeff 3 6 thole 0.829 1.590 +pair_coeff 3 7 thole 0.926 1.590 +pair_coeff 3 8 thole 0.630 0.670 +pair_coeff 6 6 thole 1.090 2.510 +pair_coeff 6 7 thole 1.218 2.510 +pair_coeff 6 8 thole 0.829 1.590 +pair_coeff 7 7 thole 1.360 2.510 +pair_coeff 7 8 thole 0.926 1.590 +pair_coeff 8 8 thole 0.630 0.670 ++
For the thole pair style the coefficients are +
+The special neighbors have charge-charge and charge-dipole +interactions screened by the coul factors of the special_bonds +command (0., 0., and 0.5 in the example above). Without using the +pair_style thole, dipole-dipole interactions are screened by the +same factor. By using the pair_style thole, dipole-dipole +interactions are screened by Thole's function, whatever their special +relationship (except within each DC-DP pair of course). Consider for +example 1-2 neighbors: using the pair_style thole, their dipoles +will see each other (despite the coul factor being 0.) and the +interactions between these dipoles will be damped by Thole's +function. +
+Thermostats and barostats +
+Using a Nose-Hoover barostat with the langevin/drude thermostat is +straightforward using fix nph instead of nve. For example: +
+fix NPH all nph iso 1. 1. 500 ++
It is also possible to use a Nose-Hoover instead of a Langevin +thermostat. This requires to use fix +drude/transform just before and after the +time intergation 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. 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 would use +
+fix DIRECT all drude/transform/direct +fix NVT1 ATOMS nvt temp 300. 300. 100 +fix NVT2 DRUDES nvt temp 1. 1. 20 +fix INVERSE all drude/transform/inverse ++
For our phenol example, the groups would be defined as +
+group ATOMS type 1 2 3 4 5 # DCs and non-polarizable atoms +group CORES type 1 2 3 # DCs +group DRUDES type 6 7 8 # DPs ++
Note that with the fixes drude/transform, it is not required to +specify comm_modify vel yes because the fixes do it anyway (several +times and for the forces also). To avoid the flying ice cube artifact +(Lamoureux), where the atoms progressively freeze and the +center of mass of the whole system drifts faster and faster, the fix +momentum can be used. For instance: +
+fix MOMENTUM all momentum 100 linear 1 1 1 ++
It is a bit more tricky to run a NPT simulation with Nose-Hoover +barostat and thermostat. First, the volume should be integrated only +once. So the fix for DCs and atoms should be npt while the fix for +DPs should be 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 then use +the fix_modify command for this. In the end, the block of +instructions for thermostating and barostating will look like +
+compute TATOMS ATOMS temp +fix DIRECT all drude/transform/direct +fix NPT ATOMS npt temp 300. 300. 100 iso 1. 1. 500 +fix_modify NPT temp TATOMS press thermo_press +fix NVT DRUDES nvt temp 1. 1. 20 +fix INVERSE all drude/transform/inverse ++
Rigid bodies +
+You may want to simulate molecules as rigid bodies (but polarizable). +Common cases are water models such as SWM4-NDP, which is a +kind of polarizable TIP4P water. The rigid bodies and the DPs should +be integrated separately, even with the Langevin thermostat. Let us +review the different thermostats and ensemble combinations. +
+NVT ensemble using Langevin thermostat: +
+comm_modify vel yes +fix LANG all langevin/drude 300. 100 12435 1. 20 13977 +fix RIGID ATOMS rigid/nve/small molecule +fix NVE DRUDES nve ++
NVT ensemble using Nose-Hoover thermostat: +
+fix DIRECT all drude/transform/direct +fix RIGID ATOMS rigid/nvt/small molecule temp 300. 300. 100 +fix NVT DRUDES nvt temp 1. 1. 20 +fix INVERSE all drude/transform/inverse ++
NPT ensemble with Langevin thermostat: +
+comm_modify vel yes +fix LANG all langevin/drude 300. 100 12435 1. 20 13977 +fix RIGID ATOMS rigid/nph/small molecule iso 1. 1. 500 +fix NVE DRUDES nve ++
NPT ensemble using Nose-Hoover thermostat: +
+compute TATOM ATOMS temp +fix DIRECT all drude/transform/direct +fix RIGID ATOMS rigid/npt/small molecule temp 300. 300. 100 iso 1. 1. 500 +fix_modify RIGID temp TATOM press thermo_press +fix NVT DRUDES nvt temp 1. 1. 20 +fix INVERSE all drude/transform/inverse ++
(Lamoureux) Lamoureux and Roux, J Chem Phys, 119, 3025-3039 (2003) +
+ + +(Schroeder) Schröder and Steinhauser, J Chem Phys, 133, +154511 (2010). +
+ + +(Jiang) Jiang, Hardy, Phillips, MacKerell, Schulten, and Roux, + J Phys Chem Lett, 2, 87-92 (2011). +
+ + +(Thole) Chem Phys, 59, 341 (1981). +
+ + +(Noskov) Noskov, Lamoureux and Roux, J Phys Chem B, 109, 6705 (2005). +
+ + +(SWM4-NDP) Lamoureux, Harder, Vorobyov, Roux, MacKerell, Chem Phys +Let, 418, 245-249 (2006) +
+ diff --git a/doc/tutorial_drude.txt b/doc/tutorial_drude.txt new file mode 100644 index 0000000000..c574395c35 --- /dev/null +++ b/doc/tutorial_drude.txt @@ -0,0 +1,456 @@ + + + +Tutorial for Thermalized Drude oscillators in LAMMPS :h3 + +This tutorial explains how to use "Drude +oscillators"_drude_oscillators.html in LAMMPS to simulate polarizable +systems. As an illustration, the input files for a simulation of 250 +phenol molecules are documented. First of all, LAMMPS has to be +compiled with the USER-DRUDE package activated. Then, the data file +and input scripts have to be modified to include the Drude dipoles and +how to handle them. + +:line + +[Overview of Drude induced dipoles] + +Polarizable atoms acquire an induced electric dipole moment under the +action of an external electric field, for example the electric field +created by the surrounding particles. Drude oscillators represent +these dipoles by two fixed charges: the core (DC) and the Drude +particle (DP) bound by a harmonic potential. The Drude particle can be +thought of as the electron cloud whose center can be displaced from +the position of the the corresponding nucleus. + +The sum of the masses of a core-Drude pair should be the mass of the +initial (unsplit) atom, \(m_C + m_D = m\). The sum of their charges +should be the charge of the initial (unsplit) atom, \(q_C + q_D = q\). +A harmonic potential between the core and Drude partners should be +present, with force constant \(k_D\) and an equilibrium distance of +zero. The (half-)stiffness of the "harmonic bond"_bond_harmonic.html +\(K_D = k_D/2\) and the Drude charge \(q_D\) are related to the atom +polarizability \(\alpha\) by + +\begin\{equation\} K_D = \frac 1 2\, \frac \{q_D^2\} \alpha +\end\{equation\} + +Ideally, the mass of the Drude particle should be small, and the +stiffness of the harmonic bond should be large, so that the Drude +particle remains close ot the core. 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 Roux"_#Lamoureux suggest adopting a global +half-stiffness, \(K_D\) = 500 kcal/(mol Å2) — +which corresponds to a force constant \(k_D\) = 4184 kJ/(mol +Å2) — for all types of core-Drude bond, a +global mass \(m_D\) = 0.4 g/mol (or u) for all types of Drude +particle, and to calculate the Drude charges for individual atom types +from the atom polarizabilities using equation (1). This choice is +followed in the polarizable CHARMM force field. :ulb,l + +"Schroeder and Steinhauser"_#Schroeder suggest adopting a global +charge \(q_D\) = -1.0e and a global mass \(m_D\) = 0.1 g/mol (or u) +for all Drude particles, and to calculate the force constant for each +type of core-Drude bond from equation (1). The timesteps 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 these force fields +hydrogen atoms are treated as non-polarizable. :ule,l + +The motion of of the Drude particles can be calculated by minimizing +the energy of the induced dipoles at each timestep, by an interative, +self-consistent procedure. The Drude particles can be massless and +therefore do not contribute to the kinetic energy. However, the +relaxed method is computationall slow. An extended-lagrangian method +can be used to calculate the positions of the Drude particles, but +this requires them to have mass. It is important in this case to +decouple the degrees of freedom associated with the Drude oscillators +from those of the normal atoms. Thermalizing the Drude dipoles at +temperatures comparable to the rest of the simulation leads to several +problems (kinetic energy transfer, very short timestep, etc.), which +can be remediated by the "cold Drude" technique ("Lamoureux and +Roux"_#Lamoureux). + +Two closely related models are used to represent polarization through +"charges on a spring": the core-shell model and the Drude +model. Although the basic idea is the same, the core-shell model is +normally used for ionic/crystalline materials, whereas the Drude model +is normally used for molecular systems and fluid states. In ionic +crystals the symmetry around each ion and the distance between them +are such that the core-shell model is sufficiently stable. But to be +applicable to molecular/covalent systems the Drude model includes two +important features: + +The possibility to thermostat the additional degrees of freedom +associated with the induced dipoles at very low temperature, in terms +of the reduced coordinates of the Drude particles with respect to +their cores. This makes the trajectory close to that of relaxed +induced dipoles. :olb,l + +The Drude dipoles on covalently bonded atoms interact too strongly +due to the short distances, so an atom may capture the Drude particle +(shell) of a neighbor, or the induced dipoles within the same molecule +may align too much. To avoid this, damping at short of the +interactions between the point charges composing the induced dipole +can be done by "Thole"_#Thole functions. :ole,l + +:line + +[Preparation of the data file] + +The data file is similar to a standard LAMMPS data file for +{atom_style full}. The DPs and the {harmonic bonds} connecting them +to their DC should appear in the data file as normal atoms and bonds. + +You can use the {polarizer} tool (Python script distributed with the +USER-DRUDE package) to convert a non-polarizable data file (here +{data.102494.lmp}) to a polarizable data file ({data-p.lmp}) + +polarizer -q -f phenol.dff data.102494.lmp data-p.lmp :pre + +This will automatically insert the new atoms and bonds. +The masses and charges of DCs and DPs are computed +from {phenol.dff}, as well as the DC-DP bond constants. The file +{phenol.dff} contains the polarizabilities of the atom types +and the mass of the Drude particles, for instance: + +# units: kJ/mol, A, deg +# kforce is in the form k/2 r_D^2 +# type m_D/u q_D/e k_D alpha/A3 thole +OH 0.4 -1.0 4184.0 0.63 0.67 +CA 0.4 -1.0 4184.0 1.36 2.51 +CAI 0.4 -1.0 4184.0 1.09 2.51 :pre + +The hydrogen atoms are absent from this file, so they will be treated +as non-polarizable atoms. In the non-polarizable data file +{data.102494.lmp}, atom names corresponding to the atom type numbers +have to be specified as comments at the end of lines of the {Masses} +section. You probably need to edit it to add these names. It should +look like + +Masses :pre + +1 12.011 # CAI +2 12.011 # CA +3 15.999 # OH +4 1.008 # HA +5 1.008 # HO :pre + +:line + +[Basic input file] + +The atom style should be set to (or derive from) {full}, so that you +can define atomic charges and molecular bonds, angles, dihedrals... + +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 recognize that you are using Drude oscillators, you should +use the fix {drude}. The command is + +fix DRUDE all drude 1 1 1 0 0 2 2 2 :pre + +or, equivalently + +fix DRUDE all drude C C C N N D D D :pre + +The 0, 1, 2 (or N, C, D) following the {drude} keyword have the +following meaning: There is one tag for each atom type. This tag is 1 +(or C) for DCs, 2 (or D) for DPs and 0 (or N) for non-polarizable +atoms. Here the atom types 1 to 3 (C and O atoms) are DC, atom types +4 and 5 (H atoms) are non-polarizable and the atom types 6 to 8 are +the newly created DPs. + +By recognizing the fix {drude}, LAMMPS will find and store matching +DC-DP pairs and will treat DP as equivalent to their DC in the +{special bonds} relations. It may be necessary to extend the space +for storing such special relations. In this case extra space should +be reserved by using the {extra} keyword of the {special_bonds} +command. With our phenol, there is 1 more special neighbor for which +space is required. Otherwise LAMMPS crashes and gives the required +value. + +special_bonds lj/coul 0.0 0.0 0.5 extra 1 :pre + +Let us assume we want to run a simple NVT simulation at 300 K. Note +that Drude oscillators need to be thermalized at a low temperature in +order to approximate a self-consistent field (SCF), therefore it is not +possible to simulate an NVE ensemble with this package. Since dipoles +are approximated by a charged DC-DP pair, the {pair_style} must +include Coulomb interactions, for instance {lj/cut/coul/long} with +{kspace_style pppm}. For example, with a cutoff of 10. and a precision +1.e-4: + +pair_style lj/cut/coul/long 10.0 +kspace_style pppm 1.0e-4 :pre + +As compared to the non-polarizable input file, {pair_coeff} lines need +to be added for the DPs. Since the DPs have no Lennard-Jones +interactions, their {epsilon} is 0. so the only {pair_coeff} line +that needs to be added is + +pair_coeff * 6* 0.0 0.0 # All-DPs :pre + +Now for the thermalization, the simplest choice is to use the "fix +langevin/drude"_fix_langevin_drude.html. + +fix LANG all langevin/drude 300. 100 12435 1. 20 13977 :pre + +This applies a Langevin thermostat at temperature 300. to the centers +of mass of the DC-DP pairs, with relaxation time 100 and with random +seed 12345. This fix applies also a Langevin thermostat at temperature +1. to the relative motion of the DPs around their DCs, with relaxation +time 20 and random seed 13977. Only the DCs and non-polarizable +atoms need to be in this fix's group. LAMMPS will thermostate the DPs +together with their DC. For this, ghost atoms need to know their +velocities. Thus you need to add the following command: + +comm_modify vel yes :pre + +In order to avoid that the center of mass of the whole system +drifts due to the random forces of the Langevin thermostat on DCs, you +can add the {zero yes} option at the end of the fix line. + +If the fix {shake} is used to constrain the C-H bonds, it should be +invoked after the fix {langevin/drude} for more accuracy. + +fix SHAKE ATOMS shake 0.0001 20 0 t 4 5 :pre + +IMPORTANT NOTE: The group of the fix {shake} must not include the DPs. +If the group {ATOMS} is defined by non-DPs atom types, you could use + +Since the fix {langevin/drude} does not perform time integration (just +modification of forces but no position/velocity updates), the fix +{nve} should be used in conjunction. + +fix NVE all nve :pre + +Finally, do not forget to update the atom type elements if you use +them in a {dump_modify ... element ...} command, by adding the element +type of the DPs. Here for instance + +dump DUMP all custom 10 dump.lammpstrj id mol type element x y z ix iy iz +dump_modify DUMP element C C O H H D D D :pre + +The input file should now be ready for use! + +You will notice that the global temperature {thermo_temp} computed by +LAMMPS is not 300. K as wanted. This is because LAMMPS treats DPs as +standard atoms in his default compute. If you want to output the +temperatures of the DC-DP pair centers of mass and of the DPs relative +to their DCs, you should use {thermo_style custom} with respectively +{f_LANG\[1\]} and {f_LANG\[2\]}. These should be close to 300. and +1. on average. + +thermo_style custom step temp f_LANG\[1\] f_LANG\[2\] :pre + +:line + +[Thole screening] + +Dipolar interactions represented by point charges on springs may not +be stable, for example if the atomic polarizability is too high for +instance, 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 exagerated. Often, +special bond relations prevent bonded neighboring atoms to see the +charge of each other's DP, so that the problem does not always appear. +It is possible to use screened dipole dipole interactions by using the +"{pair_style thole}"_pair_thole.html. 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. In our example, we would use + +pair_style hybrid/overlay lj/cut/coul/long 10.0 thole 2.6 10.0 :pre + +This tells LAMMPS that we are using two pair_styles. The first one is +as above ({lj/cut/coul/long 10.0}). The second one is a {thole} +pair_style with default screening factor 2.6 ("Noskov"_#Noskov) and +cutoff 10.0. + +Since {hybrid/overlay} does not support mixing rules, 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 {pair_coeff} section of the input file. In our +example, this will look like: + +pair_coeff 1 1 lj/cut/coul/long 0.0700 3.550 +pair_coeff 1 2 lj/cut/coul/long 0.0700 3.550 +pair_coeff 1 3 lj/cut/coul/long 0.1091 3.310 +pair_coeff 1 4 lj/cut/coul/long 0.0458 2.985 +pair_coeff 2 2 lj/cut/coul/long 0.0700 3.550 +pair_coeff 2 3 lj/cut/coul/long 0.1091 3.310 +pair_coeff 2 4 lj/cut/coul/long 0.0458 2.985 +pair_coeff 3 3 lj/cut/coul/long 0.1700 3.070 +pair_coeff 3 4 lj/cut/coul/long 0.0714 2.745 +pair_coeff 4 4 lj/cut/coul/long 0.0300 2.420 +pair_coeff * 5 lj/cut/coul/long 0.0000 0.000 +pair_coeff * 6* lj/cut/coul/long 0.0000 0.000 +pair_coeff 1 1 thole 1.090 2.510 +pair_coeff 1 2 thole 1.218 2.510 +pair_coeff 1 3 thole 0.829 1.590 +pair_coeff 1 6 thole 1.090 2.510 +pair_coeff 1 7 thole 1.218 2.510 +pair_coeff 1 8 thole 0.829 1.590 +pair_coeff 2 2 thole 1.360 2.510 +pair_coeff 2 3 thole 0.926 1.590 +pair_coeff 2 6 thole 1.218 2.510 +pair_coeff 2 7 thole 1.360 2.510 +pair_coeff 2 8 thole 0.926 1.590 +pair_coeff 3 3 thole 0.630 0.670 +pair_coeff 3 6 thole 0.829 1.590 +pair_coeff 3 7 thole 0.926 1.590 +pair_coeff 3 8 thole 0.630 0.670 +pair_coeff 6 6 thole 1.090 2.510 +pair_coeff 6 7 thole 1.218 2.510 +pair_coeff 6 8 thole 0.829 1.590 +pair_coeff 7 7 thole 1.360 2.510 +pair_coeff 7 8 thole 0.926 1.590 +pair_coeff 8 8 thole 0.630 0.670 :pre + +For the {thole} pair style the coefficients are + +the atom polarizability in units of cubic length +the screening factor of the Thole function (optional, default value specified by the pair_style command) +the cutoff (optional, default value defined by the pair_style command) :ol + +The special neighbors have charge-charge and charge-dipole +interactions screened by the {coul} factors of the {special_bonds} +command (0., 0., and 0.5 in the example above). Without using the +pair_style {thole}, dipole-dipole interactions are screened by the +same factor. By using the pair_style {thole}, dipole-dipole +interactions are screened by Thole's function, whatever their special +relationship (except within each DC-DP pair of course). Consider for +example 1-2 neighbors: using the pair_style {thole}, their dipoles +will see each other (despite the {coul} factor being 0.) and the +interactions between these dipoles will be damped by Thole's +function. + +:line + +[Thermostats and barostats] + +Using a Nose-Hoover barostat with the {langevin/drude} thermostat is +straightforward using fix {nph} instead of {nve}. For example: + +fix NPH all nph iso 1. 1. 500 :pre + +It is also possible to use a Nose-Hoover instead of a Langevin +thermostat. This requires to use "{fix +drude/transform}"_fix_drude_transform.html just before and after the +time intergation 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. 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 would use + +fix DIRECT all drude/transform/direct +fix NVT1 ATOMS nvt temp 300. 300. 100 +fix NVT2 DRUDES nvt temp 1. 1. 20 +fix INVERSE all drude/transform/inverse :pre + +For our phenol example, the groups would be defined as + +group ATOMS type 1 2 3 4 5 # DCs and non-polarizable atoms +group CORES type 1 2 3 # DCs +group DRUDES type 6 7 8 # DPs :pre + +Note that with the fixes {drude/transform}, it is not required to +specify {comm_modify vel yes} because the fixes do it anyway (several +times and for the forces also). To avoid the flying ice cube artifact +"(Lamoureux)"_#Lamoureux, where the atoms progressively freeze and the +center of mass of the whole system drifts faster and faster, the {fix +momentum} can be used. For instance: + +fix MOMENTUM all momentum 100 linear 1 1 1 :pre + +It is a bit more tricky to run a NPT simulation with Nose-Hoover +barostat and thermostat. First, the volume should be integrated only +once. So the fix for DCs and atoms should be {npt} while the fix for +DPs should be {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 then use +the {fix_modify} command for this. In the end, the block of +instructions for thermostating and barostating will look like + +compute TATOMS ATOMS temp +fix DIRECT all drude/transform/direct +fix NPT ATOMS npt temp 300. 300. 100 iso 1. 1. 500 +fix_modify NPT temp TATOMS press thermo_press +fix NVT DRUDES nvt temp 1. 1. 20 +fix INVERSE all drude/transform/inverse :pre + +:line + +[Rigid bodies] + +You may want to simulate molecules as rigid bodies (but polarizable). +Common cases are water models such as "SWM4-NDP"_#SWM4-NDP, which is a +kind of polarizable TIP4P water. The rigid bodies and the DPs should +be integrated separately, even with the Langevin thermostat. Let us +review the different thermostats and ensemble combinations. + +NVT ensemble using Langevin thermostat: + +comm_modify vel yes +fix LANG all langevin/drude 300. 100 12435 1. 20 13977 +fix RIGID ATOMS rigid/nve/small molecule +fix NVE DRUDES nve :pre + +NVT ensemble using Nose-Hoover thermostat: + +fix DIRECT all drude/transform/direct +fix RIGID ATOMS rigid/nvt/small molecule temp 300. 300. 100 +fix NVT DRUDES nvt temp 1. 1. 20 +fix INVERSE all drude/transform/inverse :pre + +NPT ensemble with Langevin thermostat: + +comm_modify vel yes +fix LANG all langevin/drude 300. 100 12435 1. 20 13977 +fix RIGID ATOMS rigid/nph/small molecule iso 1. 1. 500 +fix NVE DRUDES nve :pre + +NPT ensemble using Nose-Hoover thermostat: + +compute TATOM ATOMS temp +fix DIRECT all drude/transform/direct +fix RIGID ATOMS rigid/npt/small molecule temp 300. 300. 100 iso 1. 1. 500 +fix_modify RIGID temp TATOM press thermo_press +fix NVT DRUDES nvt temp 1. 1. 20 +fix INVERSE all drude/transform/inverse :pre + +:line + +:link(Lamoureux) +[(Lamoureux)] Lamoureux and Roux, J Chem Phys, 119, 3025-3039 (2003) + +:link(Schroeder) +[(Schroeder)] Schröder and Steinhauser, J Chem Phys, 133, +154511 (2010). + +:link(Jiang) +[(Jiang)] Jiang, Hardy, Phillips, MacKerell, Schulten, and Roux, + J Phys Chem Lett, 2, 87-92 (2011). + +:link(Thole) +[(Thole)] Chem Phys, 59, 341 (1981). + +:link(Noskov) +[(Noskov)] Noskov, Lamoureux and Roux, J Phys Chem B, 109, 6705 (2005). + +:link(SWM4-NDP) +[(SWM4-NDP)] Lamoureux, Harder, Vorobyov, Roux, MacKerell, Chem Phys +Let, 418, 245-249 (2006) +