Si cluster simulation

Dear lammps community,

I’m new to lammps and wanted to do some first calculations on Si clusters with oxygen on their surface.
I want to use the COMB potential for this, unfortunately the calculation stops after the first step with the segmentation fault error:

Step Temp E_pair E_mol TotEng Press
0 0 -238.48186 0 -238.48186 -109.78007
Segmentation fault

I let the Qeq equilibration info write to a file and I get

Charge equilibration on step 0
iteration: 0, enegtot 18.3786, enegmax 49.7514, fq deviation: 11.1742
iteration: 1, enegtot 18.0347, enegmax 43.6399, fq deviation: 9.8116
iteration: 2, enegtot 17.3506, enegmax 33.5105, fq deviation: 7.5594
iteration: 3, enegtot 16.394, enegmax 21.0967, fq deviation: 5.01134
iteration: 4, enegtot 15.6893, enegmax 9.76827, fq deviation: 2.48154
iteration: 5, enegtot 36.8182, enegmax 242.488, fq deviation: 38.6245
iteration: 6, enegtot 15.9529, enegmax 27.0752, fq deviation: 4.39309
iteration: 7, enegtot -6.93144, enegmax 2340.52, fq deviation: 48.1158
iteration: 8, enegtot 47845.8, enegmax 6.16675e+06, fq deviation: 94401.5
iteration: 9, enegtot -nan, enegmax nan, fq deviation: nan
iteration: 10, enegtot -nan, enegmax nan, fq deviation: nan
iteration: 11, enegtot -nan, enegmax nan, fq deviation: nan

I then came across a thread by K. McKenna ( Date: Wed, 20 Mar 2013 13:07:30 +0000) who faced similar problems, when adding oxygen to a Ti surface.
I therefore simplified my system and now I first want to simulate a pure Si cluster with 71 Si atoms in a quadratic box, but the calculation still stops as described above and I get the same segmentation fault error.
The charge equilibration looks fine now, though, so the oxygen is not causing the problem.

My input file looks like this:

units metal

atom_style charge
pair_style comb
boundary p p p
newton on

read_data Si71_comb.data

pair_coeff * * …/…/potentials/ffield.comb Si

group Si type 1
timestep 0.001
thermo 50
thermo_modify format float %15.8g
dump 1 Si xyz 10 Si71_GEOMETRY.xyz
dump_modify 1 format “%s %15.15g %15.15g %15.15g” element Si

fix 1 Si nvt tchain 4 temp 0.0 1000.0 100
fix 2 Si qeq/comb 50 0.0001 file charge.out
run 500
unfix 1

fix 3 Si nvt tchain 4 temp 1000.0 1000.0 100
run 100
unfix 3

fix 4 Si nvt tchain 4 temp 1000.0 500.0 100
run 100
unfix 4

fix 5 Si nvt tchain 4 temp 500.0 0.0 100
run 100
unfix 5
unfix 2

minimize 1.0e-10 1.0e-14 1000 1000

undump 1

And the data file:

71 atoms
1 atom types

-25.00000 25.00000 xlo xhi
-25.00000 25.00000 ylo yhi
-25.00000 25.00000 zlo zhi

Masses

1 28.0855

Atoms

1 1 0.0000 -5.429356 -2.714678 -2.714000
2 1 0.0000 -5.429356 -2.714678 2.714000
3 1 0.0000 -5.429356 2.714678 -2.714000
4 1 0.0000 -5.429356 2.714678 2.714000
5 1 0.0000 5.429356 -2.714678 -2.714000
6 1 0.0000 5.429356 -2.714678 2.714000
7 1 0.0000 5.429356 2.714678 -2.714000
8 1 0.0000 5.429356 2.714678 2.714000
9 1 0.0000 -5.429356 0.000000 0.000000
10 1 0.0000 5.429356 0.000000 0.000000
11 1 0.0000 -4.072017 -4.072017 1.357000
12 1 0.0000 -4.072017 4.072017 -1.357000
13 1 0.0000 4.072017 -4.072017 -1.357000
14 1 0.0000 4.072017 4.072017 1.357000
15 1 0.0000 -4.072017 -1.357339 4.072000
16 1 0.0000 -4.072017 1.357339 -4.072000
17 1 0.0000 4.072017 -1.357339 -4.072000
18 1 0.0000 4.072017 1.357339 4.072000
19 1 0.0000 -4.072017 -1.357339 -1.357000
20 1 0.0000 -4.072017 1.357339 1.357000
21 1 0.0000 4.072017 -1.357339 1.357000
22 1 0.0000 4.072017 1.357339 -1.357000
23 1 0.0000 -2.714678 -5.429356 -2.714000
24 1 0.0000 -2.714678 -5.429356 2.714000
25 1 0.0000 -2.714678 5.429356 -2.714000
26 1 0.0000 -2.714678 5.429356 2.714000
27 1 0.0000 2.714678 -5.429356 -2.714000
28 1 0.0000 2.714678 -5.429356 2.714000
29 1 0.0000 2.714678 5.429356 -2.714000
30 1 0.0000 2.714678 5.429356 2.714000
31 1 0.0000 -2.714678 -2.714678 -5.429000
32 1 0.0000 -2.714678 -2.714678 5.429000
33 1 0.0000 -2.714678 2.714678 -5.429000
34 1 0.0000 -2.714678 2.714678 5.429000
35 1 0.0000 2.714678 -2.714678 -5.429000
36 1 0.0000 2.714678 -2.714678 5.429000
37 1 0.0000 2.714678 2.714678 -5.429000
38 1 0.0000 2.714678 2.714678 5.429000
39 1 0.0000 -2.714678 -2.714678 0.000000
40 1 0.0000 -2.714678 2.714678 0.000000
41 1 0.0000 2.714678 -2.714678 0.000000
42 1 0.0000 2.714678 2.714678 0.000000
43 1 0.0000 -2.714678 0.000000 -2.714000
44 1 0.0000 -2.714678 0.000000 2.714000
45 1 0.0000 2.714678 0.000000 -2.714000
46 1 0.0000 2.714678 0.000000 2.714000
47 1 0.0000 -1.357339 -4.072017 4.072000
48 1 0.0000 -1.357339 4.072017 -4.072000
49 1 0.0000 1.357339 -4.072017 -4.072000
50 1 0.0000 1.357339 4.072017 4.072000
51 1 0.0000 -1.357339 -4.072017 -1.357000
52 1 0.0000 -1.357339 4.072017 1.357000
53 1 0.0000 1.357339 -4.072017 1.357000
54 1 0.0000 1.357339 4.072017 -1.357000
55 1 0.0000 -1.357339 -1.357339 -4.072000
56 1 0.0000 -1.357339 1.357339 4.072000
57 1 0.0000 1.357339 -1.357339 4.072000
58 1 0.0000 1.357339 1.357339 -4.072000
59 1 0.0000 -1.357339 -1.357339 1.357000
60 1 0.0000 -1.357339 1.357339 -1.357000
61 1 0.0000 1.357339 -1.357339 -1.357000
62 1 0.0000 1.357339 1.357339 1.357000
63 1 0.0000 0.000000 -5.429356 0.000000
64 1 0.0000 0.000000 5.429356 0.000000
65 1 0.0000 0.000000 -2.714678 -2.714000
66 1 0.0000 0.000000 -2.714678 2.714000
67 1 0.0000 0.000000 2.714678 -2.714000
68 1 0.0000 0.000000 2.714678 2.714000
69 1 0.0000 0.000000 0.000000 -5.429000
70 1 0.0000 0.000000 0.000000 5.429000
71 1 0.0000 0.000000 0.000000 0.000000

I first want to set up a simple routine to see if the script works, but even this is not working for me and I don’t know where the error is.
I hope someone can help me with this.
Paul

What version of LAMMPS was used? More comments below, thanks.

Ray

[...]

timestep 0.001

a time step of 1 fs is simply too large for COMB. Try 0.2 or 0.1 fs.

thermo 50
thermo_modify format float %15.8g
dump 1 Si xyz 10 Si71_GEOMETRY.xyz
dump_modify 1 format "%s %15.15g %15.15g %15.15g" element Si

fix 1 Si nvt tchain 4 temp 0.0 1000.0 100
fix 2 Si qeq/comb 50 0.0001 file charge.out

For a starting system, it is better to equilibrate charge every MD
step, i.e., set Nevery to 1 instead of 50.

Hi Ray,

I have lammps version 23 Nov 2013.

I changed the time step and the charge equilibration step, but I still have the same problem.

I also generated a core dump:

(gdb) where
#0 0x00000000005f3fdb in LAMMPS_NS::PairComb::direct(int, int, int, int, double, double, double, double, double, double, double, double, double, double&, double&) ()
#1 0x00000000005f61d0 in LAMMPS_NS::PairComb::compute(int, int) ()
#2 0x00000000004e4709 in LAMMPS_NS::Verlet::run(int) ()
#3 0x00000000006b201e in LAMMPS_NS::Run::command(int, char**) ()
#4 0x0000000000433369 in void LAMMPS_NS::Input::command_creator<LAMMPS_NS::Run>(LAMMPS_NS::LAMMPS*, int, char**) ()
#5 0x0000000000430676 in LAMMPS_NS::Input::execute_command() ()
#6 0x00000000004312d4 in LAMMPS_NS::Input::file() ()
#7 0x000000000061c546 in main ()
(gdb)

but I don’t know what this is supposed to tell me.
Paul

Something else must be going on. You can send me the full input files
and I will take a look.

Ray

Hi Ray,

I have lammps version 23 Nov 2013.

I changed the time step and the charge equilibration step, but I still have the same problem.

The problem is likely in the charge equilibration, so it might be better to start with a simpler model like tersoff and pre-equilibration the system before switching to comb.

Your starting configuration may be too far away from a state where the qeq converges to a meaningful state.

Axel.

Okay, so this has nothing to do with pair_style comb nor qeq/comb.
Your script did not work with either Tersoff or LJ/cut.

The reason is simple. 0 temperature in MD is an ill-defined quantity.
Just replace your 0.0 K with any number other than 0.0, say 0.1 if
you prefer, and it works fine.

Ray

Yes you are right, thanks a lot!
Just as an addition: I left the end temperature at 0.0 K and I get:

ERROR: Target temperature for fix nvt/npt/nph cannot be 0.0 (…/fix_nh.cpp:123)

Maybe such a message could be included for when one sets the starting temperature to 0 K.
Paul