updated documentation and examples for coreshell

This commit is contained in:
Hendrik Heenen 2017-02-14 13:14:22 +01:00
parent cb982f2f28
commit 6b923476b9
3 changed files with 146 additions and 34 deletions

View File

@ -2573,7 +2573,7 @@ well.
6.26 Adiabatic core/shell model :link(howto_26),h4
The adiabatic core-shell model by "Mitchell and
Finchham"_#MitchellFinchham is a simple method for adding
Fincham"_#MitchellFincham is a simple method for adding
polarizability to a system. In order to mimic the electron shell of
an ion, a satellite particle is attached to it. This way the ions are
split into a core and a shell where the latter is meant to react to
@ -2667,13 +2667,16 @@ bond_coeff 1 63.014 0.0
bond_coeff 2 25.724 0.0 :pre
When running dynamics with the adiabatic core/shell model, the
following issues should be considered. Since the relative motion of
the core and shell particles corresponds to the polarization, typical
thermostats can alter the polarization behaviour, meaning the shell
will not react freely to its electrostatic environment. This is
critical during the equilibration of the system. Therefore
it's typically desirable to decouple the relative motion of the
core/shell pair, which is an imaginary degree of freedom, from the
following issues should be considered. The relative motion of
the core and shell particles corresponds to the polarization,
hereby an instantaneous relaxation of the shells is approximated
and a fast core/shell spring frequency ensures a nearly constant
internal kinetic energy during the simulation.
Thermostats can alter this polarization behaviour, by scaling the
internal kinetic energy, meaning the shell will not react freely to
its electrostatic environment.
Therefore it is typically desirable to decouple the relative motion of
the core/shell pair, which is an imaginary degree of freedom, from the
real physical system. To do that, the "compute
temp/cs"_compute_temp_cs.html command can be used, in conjunction with
any of the thermostat fixes, such as "fix nvt"_fix_nh.html or "fix
@ -2704,6 +2707,22 @@ fix thermostatequ all nve # integrator as needed f
fix_modify thermoberendsen temp CSequ
thermo_modify temp CSequ # output of center-of-mass derived temperature :pre
The pressure for the core/shell system is computed via the regular
LAMMPS convention by "treating the cores and shells as individual
particles"_#MitchellFincham2. For the thermo output of the pressure
as well as for the application of a barostat, it is necessary to
use an additional "pressure"_compute_pressure compute based on the
default "temperature"_compute_temp and specifying it as a second
argument in "fix modify"_fix_modify.html and
"thermo_modify"_thermo_modify.html resulting in:
(...)
compute CSequ all temp/cs cores shells
compute thermo_press_lmp all pressure thermo_temp # pressure for individual particles
thermo_modify temp CSequ press thermo_press_lmp # modify thermo to regular pressure
fix press_bar all npt temp 300 300 0.04 iso 0 0 0.4
fix_modify press_bar temp CSequ press thermo_press_lmp # pressure modification for correct kinetic scalar :pre
If "compute temp/cs"_compute_temp_cs.html is used, the decoupled
relative motion of the core and the shell should in theory be
stable. However numerical fluctuation can introduce a small
@ -2724,24 +2743,18 @@ temp/cs"_compute_temp_cs.html command to the {temp} keyword of the
velocity all create 1427 134 bias yes temp CSequ
velocity all scale 1427 temp CSequ :pre
It is important to note that the polarizability of the core/shell
pairs is based on their relative motion. Therefore the choice of
spring force and mass ratio need to ensure much faster relative motion
of the 2 atoms within the core/shell pair than their center-of-mass
velocity. This allow the shells to effectively react instantaneously
to the electrostatic environment. This fast movement also limits the
timestep size that can be used.
To maintain the correct polarizability of the core/shell pairs, the
kinetic energy of the internal motion shall remain nearly constant.
Therefore the choice of spring force and mass ratio need to ensure
much faster relative motion of the 2 atoms within the core/shell pair
than their center-of-mass velocity. This allows the shells to
effectively react instantaneously to the electrostatic environment and
limits energy transfer to or from the core/shell oscillators.
This fast movement also dictates the timestep that can be used.
The primary literature of the adiabatic core/shell model suggests that
the fast relative motion of the core/shell pairs only allows negligible
energy transfer to the environment. Therefore it is not intended to
decouple the core/shell degree of freedom from the physical system
during production runs. In other words, the "compute
temp/cs"_compute_temp_cs.html command should not be used during
production runs and is only required during equilibration. This way one
is consistent with literature (based on the code packages DL_POLY or
GULP for instance).
energy transfer to the environment.
The mentioned energy transfer will typically lead to a small drift
in total energy over time. This internal energy can be monitored
using the "compute chunk/atom"_compute_chunk_atom.html and "compute
@ -2761,14 +2774,20 @@ command, to use as input to the "compute
chunk/atom"_compute_chunk_atom.html command to define the core/shell
pairs as chunks.
For example,
For example if core/shell pairs are the only molecules:
read_data NaCl_CS_x0.1_prop.data
compute prop all property/atom molecule
compute cs_chunk all chunk/atom c_prop
compute cstherm all temp/chunk cs_chunk temp internal com yes cdof 3.0 # note the chosen degrees of freedom for the core/shell pairs
fix ave_chunk all ave/time 10 1 10 c_cstherm file chunk.dump mode vector :pre
For example if core/shell pairs and other molecules are present:
fix csinfo all property/atom i_CSID # property/atom command
read_data NaCl_CS_x0.1_prop.data fix csinfo NULL CS-Info # atom property added in the data-file
compute prop all property/atom i_CSID
compute cs_chunk all chunk/atom c_prop
compute cstherm all temp/chunk cs_chunk temp internal com yes cdof 3.0 # note the chosen degrees of freedom for the core/shell pairs
fix ave_chunk all ave/time 10 1 10 c_cstherm file chunk.dump mode vector :pre
(...) :pre
The additional section in the date file would be formatted like this:
@ -2890,9 +2909,13 @@ Phys, 79, 926 (1983).
:link(Shinoda)
[(Shinoda)] Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004).
:link(MitchellFinchham)
[(Mitchell and Finchham)] Mitchell, Finchham, J Phys Condensed Matter,
:link(MitchellFincham)
[(Mitchell and Fincham)] Mitchell, Finchham, J Phys Condensed Matter,
5, 1031-1038 (1993).
:link(MitchellFincham2)
[(Fincham)] Fincham, Mackrodt and Mitchell, J Phys Condensed Matter,
6, 393-404 (1994).
:link(howto-Lamoureux)
[(Lamoureux and Roux)] G. Lamoureux, B. Roux, J. Chem. Phys 119, 3025 (2003)

View File

@ -40,7 +40,8 @@ thermo 50
thermo_style custom step etotal pe ke temp press &
epair evdwl ecoul elong ebond fnorm fmax vol
compute CSequ all temp/cs cores shells
compute CStemp all temp/cs cores shells
compute thermo_press_lmp all pressure thermo_temp # press for correct kinetic scalar
# output via chunk method
@ -49,16 +50,18 @@ compute CSequ all temp/cs cores shells
#compute cstherm all temp/chunk cs_chunk temp internal com yes cdof 3.0
#fix ave_chunk all ave/time 100 1 100 c_cstherm file chunk.dump mode vector
thermo_modify temp CSequ
thermo_modify temp CStemp press thermo_press_lmp
# velocity bias option
velocity all create 1427 134 dist gaussian mom yes rot no bias yes temp CSequ
velocity all scale 1427 temp CSequ
velocity all create 1427 134 dist gaussian mom yes rot no bias yes temp CStemp
velocity all scale 1427 temp CStemp
# thermostating using the core/shell decoupling
fix thermoberendsen all temp/berendsen 1427 1427 0.4
fix nve all nve
fix_modify thermoberendsen temp CSequ
fix_modify thermoberendsen temp CStemp
# 2 fmsec timestep

View File

@ -0,0 +1,86 @@
# Testsystem for core-shell model compared to Mitchel and Finchham
# Hendrik Heenen, June 2014
# ------------------------ INITIALIZATION ----------------------------
units metal
dimension 3
boundary p p p
atom_style full
# ----------------------- ATOM DEFINITION ----------------------------
fix csinfo all property/atom i_CSID
read_data data.coreshell fix csinfo NULL CS-Info
group cores type 1 2
group shells type 3 4
neighbor 2.0 bin
comm_modify vel yes
# ------------------------ FORCE FIELDS ------------------------------
kspace_style ewald 1.0e-6
pair_style born/coul/long/cs 20.0 20.0 # A, rho, sigma=0, C, D
pair_coeff * * 0.0 1.000 0.00 0.00 0.00
pair_coeff 3 3 487.0 0.23768 0.00 1.05 0.50 #Na-Na
pair_coeff 3 4 145134.0 0.23768 0.00 6.99 8.70 #Na-Cl
pair_coeff 4 4 405774.0 0.23768 0.00 72.40 145.40 #Cl-Cl
bond_style harmonic
bond_coeff 1 63.014 0.0
bond_coeff 2 25.724 0.0
# ------------------------ Equilibration Run -------------------------------
reset_timestep 0
thermo 50
thermo_style custom step etotal pe ke temp press &
epair evdwl ecoul elong ebond fnorm fmax vol
compute CStemp all temp/cs cores shells
compute thermo_press_lmp all pressure thermo_temp # press for correct kinetic scalar
# output via chunk method
#compute prop all property/atom i_CSID
#compute cs_chunk all chunk/atom c_prop
#compute cstherm all temp/chunk cs_chunk temp internal com yes cdof 3.0
#fix ave_chunk all ave/time 100 1 100 c_cstherm file chunk.dump mode vector
thermo_modify temp CStemp press thermo_press_lmp
# 2 fmsec timestep
timestep 0.002
# velocity bias option
velocity all create 1427 134 dist gaussian mom yes rot no bias yes temp CStemp
velocity all scale 1427 temp CStemp
# thermostating using the core/shell decoupling
fix thermoberendsen all temp/berendsen 1427 1427 0.4
fix nve all nve
fix_modify thermoberendsen temp CStemp
run 500
unfix thermoberendsen
unfix nve
fix npt_equ all npt temp 1427 1427 0.04 iso 0 0 0.4
fix_modify npt_equ temp CStemp press thermo_press_lmp # pressure for correct kinetic scalar
run 500
unfix npt_equ
# ------------------------ Dynamic Run -------------------------------
fix npt_dyn all npt temp 1427 1427 0.04 iso 0 0 0.4
fix_modify npt_dyn temp CStemp press thermo_press_lmp # pressure for correct kinetic scalar
run 1000