From 6b923476b96c297ec695624edd3f4b32bb5c1a57 Mon Sep 17 00:00:00 2001 From: Hendrik Heenen Date: Tue, 14 Feb 2017 13:14:22 +0100 Subject: [PATCH] updated documentation and examples for coreshell --- doc/src/Section_howto.txt | 81 ++++++++++++------- examples/coreshell/in.coreshell | 13 ++-- examples/coreshell/in.coreshell.thermostats | 86 +++++++++++++++++++++ 3 files changed, 146 insertions(+), 34 deletions(-) create mode 100644 examples/coreshell/in.coreshell.thermostats diff --git a/doc/src/Section_howto.txt b/doc/src/Section_howto.txt index 5e8f676209..cfd75d5704 100644 --- a/doc/src/Section_howto.txt +++ b/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,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) diff --git a/examples/coreshell/in.coreshell b/examples/coreshell/in.coreshell index b03946306b..b436b75a86 100644 --- a/examples/coreshell/in.coreshell +++ b/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 diff --git a/examples/coreshell/in.coreshell.thermostats b/examples/coreshell/in.coreshell.thermostats new file mode 100644 index 0000000000..28b408d3a9 --- /dev/null +++ b/examples/coreshell/in.coreshell.thermostats @@ -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