Hello everyone,
sorry for my belated answer.
I figure you want to calculate a temperature for the whole system by calculating a temperature
for the core/shell oscillators and the non core/shell pairs separately with the follwing?
compute CSequ CS temp/cs O_c O_s #CS = group of Core/Shell pairs
compute PCequ PC temp #PC = group of Point/Charge atoms (non core/shell)
variable N_CS equal count(O_c)
variable N_PC equal count(PC)
variable NTOT equal N_CS+N_PC
variable T_equ equal ((CSequN_CS)+(PCequN_PC))/NTOT
thermo_modify temp T_equ
That is in generally not necessary and might lead to the wrong results. The temp/cs command computes the
temperature of the whole system (including the non-core/shell part) you specify in the “group” keyword of
"compute ID group … ". The extra groups O_core and O_shell are merely necessary to specify for the compute
to know which atoms need to be treated as core/shell oscillators. Thus it is possible to treat a system mixed
with polarized ions via the core/shell model and non-polarized ions.
In your example NaAlSi3O8 you need to define your oxygen atoms in a core and a shell group i.e.:
group O_c type atomtype(O_core)
group O_s type atomtype(O_shell)
then the compute
compute CSequ all temp/cs O_c O_s
will compute a temperature according to the physical system NaAlSi3O8 where the movement of the oxygen ions
is defined as the center of mass movement of the core and shell particles belonging to the oxygen.
If you for example want to apply the core/shell polarizability also to silicon you add the corresponding types to
the core and shell groups i.e.:
group cores type atomtype(O_core) atomtype(Si_core)
group shells type atomtype(O_shell) atomtype(Si_shell)
then the compute
compute CSequ all temp/cs cores shells
will compute a temperature according to the physical system NaAlSi3O8 where the movement of the oxygen and silicon
ions are defined as the center of mass movement of the core and shell particles belonging to oxygen and silicon
respectively. (Naturally you would need to alter your data file accordingly in order to include the silicon.)
Now in order to let LAMMPS know that the temperature it gives out in the output corresponds to the
physical temperature you specify:
thermo_modify temp CSequ
Next for the application of a thermostat and setting velocities you also need to use the calculated temperature by the compute CSequ.
For the velocity the bias option removes the introduced (by an initial distribution of velocities) motion between two
core/shell oscillators, in order to tell the command which compute to use to obtain the bias, the compute CSequ needs to be specified.
velocity all create $T 134 dist gaussian mom yes rot no bias yes temp CSequ
The applied gaussian distribution of velocities is scaled to $T before the bias is removed. Therefore it is
necessary to rescale the velocities to the temperature of $T via:
velocity all scale $T temp CSequ
The thermostat works the same way, you need the fix_modify command to let the thermostat know that
the temperature it scales to is the physical temperature computed by CSequ:
fix thermoberendsen all temp/berendsen $T $T 0.4 # choose your thermostat with the physical target temperature of $T
fix nve all nve #integration of the equation of motion
fix_modify thermoberendsen temp CSequ
Hope that clarifies this issue and all your questions from your email.
Hendrik can comment, but I agree with Axel that the compute temp/cs doc page
only indicates it should be used for the core/shell pairs.
Sorry for that unclear part in the documentation I will update it as soon as possible.
Best
Hendrik