Non-numeric positions - simulation unstable in fix gcmc command

[email protected] <[email protected]>

Dear all,

I’m now trying to use fix gcmc command in my simulation.

My simulation box is cubic (with 555nm^3) and contains a water droplet in the center. Besides the droplet, there is some void space. Before running this control file, I’ve run my simulation in NVT ensemble for 0.5 ns to make the droplet configuration stable and everything was well.

However, when I started to use fix gcmc command in this control file, I get the error message: non-numeric positions - simulation unstable. I’ve dumped my atoms coordinates, but I didn’t see atoms out of box there.

My Lammps version is 7Aug2019 and I’ve attached my input file here. Could anyone tell me why I get this error and how to resolve this problem?

Thank you so much!!!


H2O.txt (549 Bytes)

control01 (1.95 KB) (462 KB)

You are getting too creative too soon. Whenever you see behavior like this, you need to take a couple of steps back. I isolated the problem to your SHAKE command. I then reverted your script to more closely resemble examples/in.gcmc.h2o, and it appears to work now:

#fix wshake water shake 0.0001 20 0 t 1 2 mol h2omol

#fix mynvt all nvt temp {temp} {temp} 0.1

fix wshake water shake 0.0001 20 0 b 1 a 1 mol h2omol

fix mynvt all nvt temp {temp} {temp} 100

Note also that your NVT constant was very short.


Dear Aidan,

Thank you so much for your email. I’ve tried your suggestions as control01 file attached here.

(1) However, I still get the error: non-numeric positions - simulation unstable. Am I missing something?

(2) And also, my unit in the control file is metal, so the time unit is picosecond, which is different from fs in the example input file. Hence, for your suggestions
fix mynvt all nvt temp {temp} {temp} 100

Isn’t the relaxation time so long here?

(3) btw, I also tried smaller relaxation time, 1 and 10 ps for fix nvt. But I still get the same error message.

My lammps version is 7Aug2019.

Looking forward to hearing from you.

Best regards,

control01 (1.95 KB)

1 Like

Dear Aiden,

I think I may find the reason why the system crashes. I use GCMC command for 10 molecule exchange trials every 100 timesteps (1 timestep is 1 fs). I got the result that the system’s temperature increases crazy!!! My input file and output file are as follows. And I have a question: why the temperature would increase crazy for my system? Could you please give me some sense of that?

Best regards,

#input file
#SPCE_water input file

#variable command
variable mu index -0.25
variable disp index 0.5
variable temp index 298
variable tfac equal 3.0/3.0

#global setting
units metal
dimension 3
atom_style full
pair_style lj/cut/coul/long 14.0
kspace_style ewald 0.0001
bond_style harmonic
angle_style harmonic

pair_coeff 1 1 6.73853e-3 3.16556
pair_coeff 1 2 0.0 0.0
pair_coeff 2 2 0.0 0.0
bond_coeff 1 24.0291 1.0
angle_coeff 1 1.98472 109.47

molecule h2omol H2O.txt

#MD setting
group water type 1 2
neighbor 2.0 bin
neigh_modify every 1 delay 1 check yes
velocity all create ${temp} 54654
timestep 0.0001
run_style verlet

#minimize 1.0e-4 1.0e-6 100 1000
reset_timestep 0

#rigid constrints with thermostat

fix wshake water shake 0.0001 20 0 b 1 a 1 mol h2omol
fix mynvt all nvt temp {temp} {temp} 10

#define dynamic temperature

compute_modify thermo_temp dynamic/dof yes
compute_modify mynvt_temp dynamic/dof yes
timestep 0.0001

#gcmc calculation
fix mygcmc water gcmc 100 10 0 0 54341 {temp} {mu} {disp} mol & h2omol tfac_insert {tfac} group water shake wshake

#atom counts

variable oxygen atom “type==1”
variable hydrogen atom “type==2”
group oxygen dynamic all var oxygen
group hydrogen dynamic all var hydrogen
variable nO equal count(oxygen)
variable nH equal count(hydrogen)


variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
variable racc equal f_mygcmc[8]/(f_mygcmc[7]+0.1)
thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_racc v_nO v_nH

thermo 100
#dump mydump all atom 100 dump.atom
#dump 1 all dcd 1000 SPCE_water_01.dcd
dump 2 all xyz 1000
dump_modify 2 element O H
run 10000


#output file

LAMMPS (7 Aug 2019)
using 1 OpenMP thread(s) per MPI task
Reading restart file …
restart file = 7 Aug 2019, LAMMPS = 7 Aug 2019
WARNING: Restart file used different # of processors (src/read_restart.cpp:738)
restoring atom style full from restart
orthogonal box = (-10 -10 -10) to (40 40 40)
1 by 1 by 1 MPI processor grid
restoring pair style lj/cut/coul/long/intel from restart
restoring bond style harmonic/intel from restart
restoring angle style harmonic/intel from restart
3000 atoms
2000 bonds
1000 angles
Finding 1-2 1-3 1-4 neighbors …
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
2 = max # of 1-2 neighbors
1 = max # of 1-3 neighbors
1 = max # of 1-4 neighbors
2 = max # of special neighbors
special bonds CPU = 0.000854969 secs
read_restart CPU = 0.00404406 secs
Read molecule h2omol:
3 atoms with max type 2
2 bonds with max type 1
1 angles with max type 1
0 dihedrals with max type 0
0 impropers with max type 0
3000 atoms in group water
Finding SHAKE clusters …
0 = # of size 2 clusters
0 = # of size 3 clusters
0 = # of size 4 clusters
1000 = # of frozen angles
find clusters CPU = 0.000793934 secs
WARNING: Molecule attributes do not match system attributes (src/molecule.cpp:1386)
dynamic group oxygen defined
dynamic group hydrogen defined
Ewald initialization …
using 12-bit tables for long-range coulomb (src/kspace.cpp:323)
G vector (1/distance) = 0.200214
estimated absolute RMS force accuracy = 0.00166133
estimated relative force accuracy = 0.000115373
KSpace vectors: actual max1d max3d = 1054 8 2456
kxmax kymax kzmax = 8 8 8

Let me repreat my previous advice. You are getting too creative too soon. Whenever you see behavior like this, you need to take a couple of steps back. It is not possible for me or anyone else to exhaustively search for every possible mistake you may have made, in the hope of finding actual mistakes. Instead the onus is on you to start with a known working example and incrementally add desired changes one step at a time. It is a lot easier to find a mistake if there is only one possible place to look for it.

1 Like

I have already responded to this issue twice.

Dear Aidan,

Thank you so much for your kind reply. You’re right! I run my simulation on my laptop successfully. And I finally found what was going wrong.

Here, I just want to report an issue about gcmc on intel-mpi. Before, I run my simulation in the cluster with mpi-intel. On the cluster, it used default “lj/cut/coul/long/intel style” to run my simulation.

With the same input file, data file and LAMMPS-7AUG2019 version, my simulation crashes on the cluster but ran successfully on my laptop. I think there maybe something wrong when lammps run gcmc command on intel-mpi. Just let you know.

Thank you again and have a great night!

Best wishes,