From c6c44b5fe065780d2d1dcaf6fffb9e786528d152 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 24 Jul 2015 17:57:41 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13713 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- doc/Manual.html | 18 +- doc/Manual.txt | 5 + doc/Section_commands.html | 69 +++--- doc/Section_commands.txt | 9 + doc/compute_reduce.html | 9 +- doc/compute_reduce.txt | 9 +- doc/compute_temp_drude.html | 62 +++++ doc/compute_temp_drude.txt | 57 +++++ doc/fix_drude.html | 60 +++++ doc/fix_drude.txt | 56 +++++ doc/fix_drude_transform.html | 195 +++++++++++++++ doc/fix_drude_transform.txt | 183 ++++++++++++++ doc/fix_langevin_drude.html | 272 ++++++++++++++++++++ doc/fix_langevin_drude.txt | 253 +++++++++++++++++++ doc/pair_bop.html | 4 +- doc/pair_bop.txt | 4 +- doc/pair_born.html | 15 +- doc/pair_born.txt | 14 +- doc/pair_buck.html | 20 +- doc/pair_buck.txt | 19 +- doc/pair_coul.html | 21 +- doc/pair_coul.txt | 20 +- doc/tutorial_drude.html | 464 +++++++++++++++++++++++++++++++++++ doc/tutorial_drude.txt | 456 ++++++++++++++++++++++++++++++++++ 24 files changed, 2221 insertions(+), 73 deletions(-) create mode 100644 doc/compute_temp_drude.html create mode 100644 doc/compute_temp_drude.txt create mode 100644 doc/fix_drude.html create mode 100644 doc/fix_drude.txt create mode 100644 doc/fix_drude_transform.html create mode 100644 doc/fix_drude_transform.txt create mode 100644 doc/fix_langevin_drude.html create mode 100644 doc/fix_langevin_drude.txt create mode 100644 doc/tutorial_drude.html create mode 100644 doc/tutorial_drude.txt diff --git a/doc/Manual.html b/doc/Manual.html index fe2a778867..71715fa4d5 100644 --- a/doc/Manual.html +++ b/doc/Manual.html @@ -199,13 +199,13 @@ it gives quick access to documentation for all LAMMPS commands.
6.21 Calculating viscosity
- 6.22 Calculating a diffusion coefficient + 6.22 Calculating a diffusion coefficient
- 6.23 Using chunks to calculate system properties + 6.23 Using chunks to calculate system properties
- 6.24 Setting parameters for pppm/disp + 6.24 Setting parameters for pppm/disp
- 6.25 Adiabatic core/shell model + 6.25 Adiabatic core/shell model
  • Example problems @@ -412,6 +412,16 @@ it gives quick access to documentation for all LAMMPS commands. + + + + + + + + + + diff --git a/doc/Manual.txt b/doc/Manual.txt index 86f93af678..616bf4c5ec 100644 --- a/doc/Manual.txt +++ b/doc/Manual.txt @@ -228,6 +228,11 @@ it gives quick access to documentation for all LAMMPS commands. :link(howto_19,Section_howto.html#howto_19) :link(howto_20,Section_howto.html#howto_20) :link(howto_21,Section_howto.html#howto_21) +:link(howto_22,Section_howto.html#howto_22) +:link(howto_23,Section_howto.html#howto_23) +:link(howto_24,Section_howto.html#howto_24) +:link(howto_25,Section_howto.html#howto_25) +:link(howto_26,Section_howto.html#howto_26) :link(mod_1,Section_modify.html#mod_1) :link(mod_2,Section_modify.html#mod_2) diff --git a/doc/Section_commands.html b/doc/Section_commands.html index cad8cdfffa..a1c750c741 100644 --- a/doc/Section_commands.html +++ b/doc/Section_commands.html @@ -433,12 +433,13 @@ g = GPU, i = USER-INTEL, k = KOKKOS, o = USER-OMP, t = OPT. package.

    - - - - - - + + + + + +
    adapt/fepaddtorqueatcave/spatial/spherecolvarsgle
    imdipilangevin/efflb/fluidlb/momentumlb/pc
    lb/rigid/pc/spherelb/viscousmesomeso/stationarynph/effnpt/eff
    nve/effnvt/effnvt/sllod/effphononpimdqbmsst
    qeq/reaxqmmmqtbreax/c/bondsreax/c/speciessaed/vtk
    smdtemp/rescale/effti/rsti/springttm/mod +
    adapt/fepaddtorqueatcave/spatial/spheredrudedrude/transform/direct
    drude/transform/reversecolvarsgleimdipilangevin/drude
    langevin/efflb/fluidlb/momentumlb/pclb/rigid/pc/spherelb/viscous
    mesomeso/stationarynph/effnpt/effnve/effnvt/eff
    nvt/sllod/effphononpimdqbmsstqeq/reaxqmmm
    qtbreax/c/bondsreax/c/speciessaed/vtksmdtemp/rescale/eff
    ti/rsti/springttm/mod
    @@ -473,8 +474,8 @@ package.

    - - +
    ackland/atombasal/atomfepke/effke/atom/effmeso_e/atom
    meso_rho/atommeso_t/atomsaedtemp/efftemp/deform/efftemp/region/eff
    temp/rotatexrd +
    meso_rho/atommeso_t/atomsaedtemp/drudetemp/efftemp/deform/eff
    temp/region/efftemp/rotatexrd
    @@ -492,30 +493,31 @@ KOKKOS, o = USER-OMP, t = OPT.
    - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + +
    nonehybridhybrid/overlayadp (o)
    airebo (o)beck (go)bodybop
    born (go)born/coul/long (cgo)born/coul/msm (o)born/coul/wolf (go)
    brownian (o)brownian/poly (o)buck (cgko)buck/coul/cut (cgko)
    buck/coul/long (cgko)buck/coul/msm (o)buck/long/coul/long (o)colloid (go)
    comb (o)comb3coul/cut (gko)coul/debye (gko)
    coul/dsf (gko)coul/long (gko)coul/msmcoul/streitz
    coul/wolf (ko)dpd (o)dpd/tstat (o)dsmc
    eam (cgkot)eam/alloy (cgkot)eam/fs (cgkot)eim (o)
    gauss (go)gayberne (gio)gran/hertz/history (o)gran/hooke (co)
    gran/hooke/history (o)hbond/dreiding/lj (o)hbond/dreiding/morse (o)kim
    lcbopline/lj (o)lj/charmm/coul/charmm (cko)lj/charmm/coul/charmm/implicit (cko)
    lj/charmm/coul/long (cgiko)lj/charmm/coul/msmlj/class2 (cgko)lj/class2/coul/cut (cko)
    lj/class2/coul/long (cgko)lj/cut (cgikot)lj/cut/coul/cut (cgko)lj/cut/coul/debye (cgko)
    lj/cut/coul/dsf (gko)lj/cut/coul/long (cgikot)lj/cut/coul/msm (go)lj/cut/dipole/cut (go)
    lj/cut/dipole/longlj/cut/tip4p/cut (o)lj/cut/tip4p/long (ot)lj/expand (cgko)
    lj/gromacs (cgko)lj/gromacs/coul/gromacs (cko)lj/long/coul/long (o)lj/long/dipole/long
    lj/long/tip4p/longlj/smooth (co)lj/smooth/linear (o)lj96/cut (cgo)
    lubricate (o)lubricate/poly (o)lubricateUlubricateU/poly
    meam (o)mie/cut (o)morse (cgot)nb3b/harmonic (o)
    nm/cut (o)nm/cut/coul/cut (o)nm/cut/coul/long (o)peri/eps
    peri/lps (o)peri/pmb (o)peri/vesreax
    rebo (o)resquared (go)snapsoft (go)
    sw (cgkio)table (gko)tersoff (cko)tersoff/mod (ko)
    tersoff/zbl (ko)tip4p/cut (o)tip4p/long (o)tri/lj (o)
    yukawa (go)yukawa/colloid (go)zbl (o) +
    born (go)born/coul/long (cgo)born/coul/long/csborn/coul/msm (o)
    born/coul/wolf (go)brownian (o)brownian/poly (o)buck (cgko)
    buck/coul/cut (cgko)buck/coul/long (cgko)buck/coul/long/csbuck/coul/msm (o)
    buck/long/coul/long (o)colloid (go)comb (o)comb3
    coul/cut (gko)coul/debye (gko)coul/dsf (gko)coul/long (gko)
    coul/long/cscoul/msmcoul/streitzcoul/wolf (ko)
    dpd (o)dpd/tstat (o)dsmceam (cgkot)
    eam/alloy (cgkot)eam/fs (cgkot)eim (o)gauss (go)
    gayberne (gio)gran/hertz/history (o)gran/hooke (co)gran/hooke/history (o)
    hbond/dreiding/lj (o)hbond/dreiding/morse (o)kimlcbop
    line/lj (o)lj/charmm/coul/charmm (cko)lj/charmm/coul/charmm/implicit (cko)lj/charmm/coul/long (cgiko)
    lj/charmm/coul/msmlj/class2 (cgko)lj/class2/coul/cut (cko)lj/class2/coul/long (cgko)
    lj/cut (cgikot)lj/cut/coul/cut (cgko)lj/cut/coul/debye (cgko)lj/cut/coul/dsf (gko)
    lj/cut/coul/long (cgikot)lj/cut/coul/msm (go)lj/cut/dipole/cut (go)lj/cut/dipole/long
    lj/cut/tip4p/cut (o)lj/cut/tip4p/long (ot)lj/expand (cgko)lj/gromacs (cgko)
    lj/gromacs/coul/gromacs (cko)lj/long/coul/long (o)lj/long/dipole/longlj/long/tip4p/long
    lj/smooth (co)lj/smooth/linear (o)lj96/cut (cgo)lubricate (o)
    lubricate/poly (o)lubricateUlubricateU/polymeam (o)
    mie/cut (o)morse (cgot)nb3b/harmonic (o)nm/cut (o)
    nm/cut/coul/cut (o)nm/cut/coul/long (o)peri/epsperi/lps (o)
    peri/pmb (o)peri/vesreaxrebo (o)
    resquared (go)snapsoft (go)sw (cgkio)
    table (gko)tersoff (cko)tersoff/mod (ko)tersoff/zbl (ko)
    tip4p/cut (o)tip4p/long (o)tri/lj (o)yukawa (go)
    yukawa/colloid (go)zbl (o)

    These are additional pair styles in USER packages, which can be used @@ -530,7 +532,8 @@ package. lj/sdk/coul/long (go)lj/sdk/coul/msm (o)lj/sf (o)meam/spline meam/sw/splinequipreax/csph/heatconduction sph/idealgassph/ljsph/rhosumsph/taitwater -sph/taitwater/morrissrptersoff/table (o)tip4p/long/soft (o) +sph/taitwater/morrissrptersoff/table (o)thole +tip4p/long/soft (o)

    diff --git a/doc/Section_commands.txt b/doc/Section_commands.txt index ad18e51b3e..600ef8038e 100644 --- a/doc/Section_commands.txt +++ b/doc/Section_commands.txt @@ -602,10 +602,14 @@ package"_Section_start.html#start_3. "addtorque"_fix_addtorque.html, "atc"_fix_atc.html, "ave/spatial/sphere"_fix_ave_spatial_sphere.html, +"drude"_fix_drude.html, +"drude/transform/direct"_fix_drude_transform.html, +"drude/transform/reverse"_fix_drude_transform.html, "colvars"_fix_colvars.html, "gle"_fix_gle.html, "imd"_fix_imd.html, "ipi"_fix_ipi.html, +"langevin/drude"_fix_langevin_drude.html, "langevin/eff"_fix_langevin_eff.html, "lb/fluid"_fix_lb_fluid.html, "lb/momentum"_fix_lb_momentum.html, @@ -726,6 +730,7 @@ package"_Section_start.html#start_3. "meso_rho/atom"_compute_meso_rho_atom.html, "meso_t/atom"_compute_meso_t_atom.html, "saed"_compute_saed.html, +"temp/drude"_compute_temp_drude.html, "temp/eff"_compute_temp_eff.html, "temp/deform/eff"_compute_temp_deform_eff.html, "temp/region/eff"_compute_temp_region_eff.html, @@ -754,6 +759,7 @@ KOKKOS, o = USER-OMP, t = OPT. "bop"_pair_bop.html, "born (go)"_pair_born.html, "born/coul/long (cgo)"_pair_born.html, +"born/coul/long/cs"_pair_born.html, "born/coul/msm (o)"_pair_born.html, "born/coul/wolf (go)"_pair_born.html, "brownian (o)"_pair_brownian.html, @@ -761,6 +767,7 @@ KOKKOS, o = USER-OMP, t = OPT. "buck (cgko)"_pair_buck.html, "buck/coul/cut (cgko)"_pair_buck.html, "buck/coul/long (cgko)"_pair_buck.html, +"buck/coul/long/cs"_pair_buck.html, "buck/coul/msm (o)"_pair_buck.html, "buck/long/coul/long (o)"_pair_buck_long.html, "colloid (go)"_pair_colloid.html, @@ -770,6 +777,7 @@ KOKKOS, o = USER-OMP, t = OPT. "coul/debye (gko)"_pair_coul.html, "coul/dsf (gko)"_pair_coul.html, "coul/long (gko)"_pair_coul.html, +"coul/long/cs"_pair_coul.html, "coul/msm"_pair_coul.html, "coul/streitz"_pair_coul.html, "coul/wolf (ko)"_pair_coul.html, @@ -883,6 +891,7 @@ package"_Section_start.html#start_3. "sph/taitwater/morris"_pair_sph_taitwater_morris.html, "srp"_pair_srp.html, "tersoff/table (o)"_pair_tersoff.html, +"thole"_pair_thole.html, "tip4p/long/soft (o)"_pair_lj_soft.html :tb(c=4,ea=c) :line diff --git a/doc/compute_reduce.html b/doc/compute_reduce.html index 6a02a01eea..e5eb69ea9c 100644 --- a/doc/compute_reduce.html +++ b/doc/compute_reduce.html @@ -25,7 +25,7 @@ reduce/region arg = region-ID region-ID = ID of region to use for choosing atoms -
  • mode = sum or min or max or ave +
  • mode = sum or min or max or ave or sumsq or avesq
  • one or more inputs can be listed @@ -71,7 +71,12 @@ vectors. 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_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 @@ + +

    LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Commands +
    + + + + + + +
    + +

    compute temp/drude command +

    +

    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 @@ + +
    LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Commands +
    + + + + + + +
    + +

    fix drude command +

    +

    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 @@ + + + + +
    LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Commands +
    + + + + + + +
    + +

    fix drude/transform/direct command +

    +

    fix drude/transform/inverse command +

    +

    Syntax: +

    +
    fix ID group-ID style keyword value ... 
+
    + +

    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 +

    +
    1. temperature of the centers of mass (temperature units) +
    2. temperature of the dipoles (temperature units) +
    3. number of degrees of freedom of the centers of mass +
    4. number of degrees of freedom of the dipoles +
    5. kinetic energy of the centers of mass (energy units) +
    6. kinetic energy of the dipoles (energy units) +
    +

    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 @@ + + + + +
    LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Commands +
    + + + + + + +
    + +

    fix langevin/drude command +

    +

    Syntax: +

    +
    fix ID group-ID langevin/drude Tcom damp_com seed_com Tdrude damp_drude seed_drude keyword values ... 
+
    + +

    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} } + } \). +
    +\(f_r'\) is a random force proportional to +\(\sqrt { \frac {2\, k_B \mathtt{Tdrude}\, m'} + {\mathrm dt\, \mathtt{damp\_drude} } + } \). +
    +

    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: +

    +
    1. temperature of the centers of mass (temperature units) +
    2. temperature of the dipoles (temperature units) +
    3. number of degrees of freedom of the centers of mass +
    4. number of degrees of freedom of the dipoles +
    5. kinetic energy of the centers of mass (energy units) +
    6. kinetic energy of the dipoles (energy units) +
    +

    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 born/coul/long command

    +

    pair_style born/coul/long/cs command +

    pair_style born/coul/long/cuda command

    pair_style born/coul/long/gpu command @@ -37,12 +39,12 @@ 

    pair_style style args
    -
    • 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 
        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 buck/coul/long command

    +

    pair_style buck/coul/long/cs command +

    pair_style buck/coul/long/cuda command

    pair_style buck/coul/long/gpu command @@ -47,7 +49,7 @@ 

    pair_style style args
    -
    • 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 
        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 command

    +

    pair_style coul/long/cs command +

    pair_style coul/long/omp command

    pair_style coul/long/gpu command @@ -67,6 +69,7 @@ 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 @@ -91,6 +94,7 @@ pair_coeff 2 2 3.5 pair_coeff * * 
    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.
+
    
+
    • Lamoureux and Roux 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. 
+
+
    • Schroeder and Steinhauser 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. 
+
    
+
    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:
+
    
+
    1. 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. 
+
+
    2. 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 functions. 
+
    
+
    
+
+
    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 
+
    
+
    1. the atom polarizability in units of cubic length
+
    2. the screening factor of the Thole function (optional, default value specified by the pair_style command)
+
    3. the cutoff (optional, default value defined by the pair_style command) 
+
    
+
    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)
+