Problems restarting dynamics

Dear LAMMPS users.
I am new user of LAMMPS. I am doing dynamics that consists in heat a nanoparticle (NP) immersed in liquid. I do dynamics in the next way: firstly I thermalize the system in 100,000 steps using the NPT assemble with a temperature equal to 300 K and an atmosphere and write a restart finishing the dynamic. Then, I do other dynamic, also of 100,000 steps, where I use region commands: one for a sink located far away from the NP and other for the water between the sink and NP. For the NP I use a fix with NVT assemble to raise its temperature from 300 K to 500 K in 100,000 steps, for the region between NP and the sink I use the fix with NPH with an atmosphere, and for sink region I use the fix with NPT assemble with 300 K and an atmosphere, and I write a restart when the dynamic is finished. Until here everything is fine.

The problems comes when I try restart the dynamic. I follow the same way for the regions like above, the unique change that I do is for the fix for the NP. I tell to LAMMPS that keeps the temperature of 500 K reached in the NP, but the dynamic crash in a very few steps (like 24). Water temperature raise too fast and the total pressure fall rapidly. I have try some thing like change the fix style like nve or mvv/dpd for the water in the region between NP and sink but always happens the same. The unique that I have ‘achieved’ is when I decrease the coupling values for the thermostats: the dynamic runs more steps but finally happens the same. I things that is happening some related with technical issues of LAMMPS more than a physical situation, because the firsts dynamics run well and only crash when I change the fix value for the NP. May be is related with the write_command or some like this. I write my input files for you look the unique change that I do

Input to raise the NP temperature after have reached the equilibrium:

variable T index 300
variable T_NP index 500
variable root index GN_tip3pT300

units real
boundary p p p

atom_style full
bond_style harmonic
angle_style harmonic
pair_style lj/cut/coul/long 10.0
kspace_style pppm 1.0e-4

read_restart restart.tip3p_GN.1000
#read_data GN_tip3p.dat

region 1st_shell sphere 0.0 0.0 0.0 9.0 units box side out
region 6th_shell sphere 0.0 0.0 0.0 18.0 units box side in
region cores intersect 2 1st_shell 6th_shell
region outer_shell sphere 0.0 0.0 0.0 18.0 units box side out

pair_coeff 1 1 0.0 0.0
pair_coeff 1 2 0.0 0.0
pair_coeff 1 3 0.0 0.0
pair_coeff 2 2 0.1520721781581 3.15066
pair_coeff 2 3 0.8969153085504 2.878057
pair_coeff 3 3 5.2899687470000 2.62904

bond_coeff 1 600.398055 0.9572
angle_coeff 1 75.0498354 104.52

group water type 1 2
group hydrogen type 1
group oxygen type 2
group gold type 3

group ver region cores
group vver intersect ver water
group shell_out region outer_shell
group sink intersect shell_out water

neighbor 2.0 bin
neigh_modify delay 5 every 1 check yes

compute temp_gold gold temp
compute temp_water water temp
compute temp all temp
compute temp_sink sink temp

compute cc1 all chunk/atom bin/sphere 0.0 0.0 0.0 0.05 20.0 15 discard yes
fix 1 all ave/chunk 100 10 1000 cc1 density/number temp norm sample file radial_profile.dat

fix raise_NP gold nvt temp 300.0 500.0 100.0
fix 34 vver nph iso 1.0 1.0 500.0
fix sink_pp sink npt temp $T $T 100.0 iso 1.0 1.0 500.0
fix recentergold gold recenter 0.0 0.0 0.0 shift all units box

dump dumpxyz all xyz 5000 dump.${root}$T.xyz
dump_modify dumpxyz flush yes

log thermo.${root}T$T

thermo 1000
thermo_modify flush yes

timestep 1.0

thermo_style custom step c_temp_gold c_temp_water c_temp c_temp_sink press vol #c_cc1
thermo 10

run 100000

write_restart restart.tip3p_GN_2

Input to keep the NP temperature and get data:

variable T index 300
variable T_NP index 500
variable root index GN_tip3pT300

units real
boundary p p p

atom_style full
bond_style harmonic
angle_style harmonic
pair_style lj/cut/coul/long 10.0
kspace_style pppm 1.0e-4

read_restart restart.tip3p_GN_2
#read_data GN_tip3p.dat

region 1st_shell sphere 0.0 0.0 0.0 9.0 units box side out
region 6th_shell sphere 0.0 0.0 0.0 18.0 units box side in
region cores intersect 2 1st_shell 6th_shell
region outer_shell sphere 0.0 0.0 0.0 18.0 units box side out

pair_coeff 1 1 0.0 0.0
pair_coeff 1 2 0.0 0.0
pair_coeff 1 3 0.0 0.0
pair_coeff 2 2 0.1520721781581 3.15066
pair_coeff 2 3 0.8969153085504 2.878057
pair_coeff 3 3 5.2899687470000 2.62904

bond_coeff 1 600.398055 0.9572
angle_coeff 1 75.0498354 104.52

group water type 1 2
group hydrogen type 1
group oxygen type 2
group gold type 3

group ver region cores
group vver intersect ver water
group shell_out region outer_shell
group sink intersect shell_out water

neighbor 2.0 bin
neigh_modify delay 5 every 1 check yes

compute temp_gold gold temp
compute temp_water water temp
compute temp all temp
compute temp_sink sink temp

compute cc1 all chunk/atom bin/sphere 0.0 0.0 0.0 0.05 20.0 15 discard yes
fix 1 all ave/chunk 100 10 1000 cc1 density/number temp norm sample file radial_profile.dat

fix raise_NP gold nvt temp 500.0 500.0 100.0
fix 34 vver nph iso 1.0 1.0 500.0
fix sink_pp sink npt temp $T $T 100.0 iso 1.0 1.0 500.0
fix recentergold gold recenter 0.0 0.0 0.0 shift all units box

dump dumpxyz all xyz 5000 dump.${root}$T.xyz
dump_modify dumpxyz flush yes

log thermo.${root}T$T

thermo 1000
thermo_modify flush yes

timestep 1.0

thermo_style custom step c_temp_gold c_temp_water c_temp c_temp_sink press vol #c_cc1
thermo 10

run 100000

write_restart restart.tip3p_GN_3

Thank you so much to all for your help.

Best regards,

Oscar Gutiérrez.

It looks like you have both a fix nph and fix not command
defined. That’s bad for dynamics as both will try to adjust
the box size (to control the pressure). The pressure
in LAMMPS is global (regardless of the group used
in the fix nph or fix not command) and it should
be controlled exactly once (or not at all). All other
integration fixes should only control temperature (e.g.
for other groups of atoms).

Axel - should we add an error check that a script
does not define multiple pressure-control fixes?

Steve

Axel - should we add an error check that a script
does not define multiple pressure-control fixes?

Steve,

we already have a check in fix_nh.cpp that checks for collisions with fix deform on individual box change parameters.
to make this more general, would be a bit more complex. let me think about it.

axel.

I was thinking something simpler, independent of fix deform.
Just a check that multiple npt or nph fixes are not defined.
Maybe it is OK to have 2 npt fixes is they do not control
the same dimension. Though that’s an odd thing to do, so

maybe a warning instead of an error is sufficient.

Also checking there is no conflict with fix deform in the dimensions
controlled by npt vs deform would be a bonus. I think you’re
saying that is already in place.

Steve

I was thinking something simpler, independent of fix deform.
Just a check that multiple npt or nph fixes are not defined.

this may look simple, but the number of possible permutations are growing significantly with each new feature that would modify box information.
depending on string matches of the fix name is not a very safe thing to do in the first place, and it gets pretty complex, if you consider that there will be accelerated versions and rigid versions of everything.

so doing a more generic check (and getting rid of multiple variables in the fix class by replacing them with a single bitmask) is the approach that looks much better to me in the long term. please see PR #1874 https://github.com/lammps/lammps/pull/1874
as indicated, this opens the door to additional checks in a generic place and thus we could get rid of the more specific checks that will only test for subsets.

Maybe it is OK to have 2 npt fixes is they do not control
the same dimension. Though that’s an odd thing to do, so

maybe a warning instead of an error is sufficient.

Also checking there is no conflict with fix deform in the dimensions
controlled by npt vs deform would be a bonus. I think you’re
saying that is already in place.

there are some, but not everywhere.

axel.