Dear all,
I was trying to do shock simulation study of water by applying MSST method in the lammps frame work. Here is the details of my system:

sample: Water

Forcefield  Reaxc/ ffield.reax.Fe_O_C_H

Equlibration (NPT) have done with reaxc/qeq

Simulation box p p p with size ~252525 A

Total no of atoms: 1500

Condition: 300K, 1 atm

Input script for msst shock:
REAX potential for H / O system
Using for bulk water
general setup 
units real
atom_style charge
read_restart 300000.restart
pair_style reax/c NULL
pair_coeff * * ffield.reax.Fe_O_C_H H O
neighbor 2 bin
neigh_modify every 1 delay 0 check yes
#########outputs
variable sysvol equal vol
variable sysmass equal mass(all)/6.0221367e+23
variable sysdensity equal v_sysmass/v_sysvol/1.0e24
variable time equal step*dt+0.000001
fix 2 all qeq/reax 1 0.0 10.0 1.0e6 reax/c
fix 7 all ave/time 100 10 1000 c_thermo_temp c_thermo_press v_sysdensity file thermo.output.dat
MSST fix
1km/s = .01A/fs
fix msst all msst z 0.15 q 5 tscale 0.2
this is needed to make etotal equal the MSST conserved quantity
fix_modify msst energy yes
variable dhug equal f_msst[1]
variable dray equal f_msst[2]
variable lgr_vel equal f_msst[3]
variable lgr_pos equal f_msst[4]
thermo 10
thermo_style custom step temp ke pe lx ly lz pxx pyy pzz etotal v_dhug v_dray v_lgr_vel v_lgr_pos f_msst
thermo_modify flush yes
timestep 0.05
restart 100000 *.restart
write_data data.final
run 1000000
1 Like
You should use the smallest possible tscale that results in compression.
Masslike parameter q should be chosen to give a rise time on the
timescale of 100fs to 1 ps  meaning compression occurs after 100fs 
1ps.
If your system compress almost immediately, then your q and tscale are
not well chosen.
Your compression rate (Lz/Lz0 = 9/25) is almost 65%, which is
definitely too high. I would choose a smaller shock velocity (why Us
= 15 km/s by the way? That is most likely overdriven and will be
difficult to reach steady state). When you see departures from the
Hugoniot and the Rayleigh line, most definitely they won't come back 
meaning your parameters (tscale, q, Us) are not good.
Ray
Hi Ray,
After a sets of MSST simulation by adjusting “q” and “tscale”, I got some results regarding shock compression of water modelled by reax/c. I have attached the plots of the variables. I have read your PETN study by MSST and reaxFF recently. I have some specific Qs related to MSST framework. Would be highly obliged if you kindly make some suggestion and comment regarding the results (attached) and confusions,

From my results I got that the steady state have reached after ~5 ps but as per the existing literature ( JCP 130 124517 2009, attached ) is concerned i saw that the steady state can reach after ~1 ps. So my confusion is what the physical significance of this “time” required to attain equilibrium state, when we take the averaged data of the variables at the equilibrated portion of the curves ?

I used q=50 Us=5km/s tscale=0.004, if I want to reduce the time of damping which parameter is the most sensitive to be adjusted ?

As I’ve no experience with MSST, so I don’t know whether the shifts of the dhugoniot and drayleigh line is look normal in my plots or not .

Using abinitio + MSST, the reported value of P = 8.3 GPa and T= 434 K but I got 10 GPa and 760 K respectively. Is this due to the wrong parameters ?

what is the f_msst in the thermo_style, is there any documentation about this f_msst ?
Thank you for your help.
shock_msst.pdf (461 KB)
jcp.pdf (501 KB)
Hi Ray,
After a sets of MSST simulation by adjusting "q" and "tscale", I got some
results regarding shock compression of water modelled by reax/c. I have
attached the plots of the variables. I have read your PETN study by MSST and
reaxFF recently. I have some specific Qs related to MSST framework. Would be
highly obliged if you kindly make some suggestion and comment regarding the
results (attached) and confusions,
1. From my results I got that the steady state have reached after ~5 ps but
as per the existing literature ( JCP 130 124517 2009, attached ) is
concerned i saw that the steady state can reach after ~1 ps. So my confusion
is what the physical significance of this "time" required to attain
equilibrium state, when we take the averaged data of the variables at the
equilibrated portion of the curves ?
Yes.
2. I used q=50 Us=5km/s tscale=0.004, if I want to reduce the time of
damping which parameter is the most sensitive to be adjusted ?
Not sure, may consider the mu parameter. But why reduce oscillation
time? From the doc page:
"Under some special highsymmetry conditions, the pressure (volume)
and/or temperature of the system may oscillate for many cycles even
with an appropriate choice of masslike parameter q. Such oscillations
have physical significance in some cases. The optional mu keyword adds
an artificial viscosity that helps break the system symmetry to
equilibrate to the shock Hugoniot and Rayleigh line more rapidly in
such cases."
3. As I've no experience with MSST, so I don't know whether the shifts of
the dhugoniot and drayleigh line is look normal in my plots or not .
Are they around zero?
4. Using abinitio + MSST, the reported value of P = 8.3 GPa and T= 434 K
but I got 10 GPa and 760 K respectively. Is this due to the wrong parameters
?
One is DFT and one is ReaxFF  why should they be the same?
5. what is the f_msst in the thermo_style, is there any documentation about
this f_msst ?
It is explaines in the doc page: look out for the word scalar.
Ray