Merge pull request #380 from hheenen/core_shell_documentation

updated documentation and examples for coreshell

doc/src/Section_howto.txt

@@ -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,15 +2774,21 @@ 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:
 
-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
+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
+(...) :pre
+
 The additional section in the date file would be formatted like this:
 
 CS-Info         # header of additional section :pre
@@ -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, Fincham, 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)

examples/coreshell/in.coreshell

@@ -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
 

examples/coreshell/in.coreshell.thermostats

+# Testsystem for core-shell model compared to Mitchell and Fincham
+# 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