Dear all,
I was trying to develop charges on a crystal by using Qeq. I found the ionization potential and electron affinity parameters to get the Mulliken electronegativity, which is used by Qeq. But it gave NaN on charges and those energies. The crystal structure is from the standard database. It seems there is something wrong blowing up the system. Can anyone know how to deal with it? I was using the same method successfully once before, but I have no idea why this time it is not working.
Here is the input.
units metal
dimension 3
#boundary p p p
atom_style charge
read_data bi2o3_qeq.data
pair_style buck/coul/long 12.0
kspace_style ewald 1.0e-6
pair_coeff 1 1 6452.86 0.5089 0.0000
pair_coeff 1 2 6048.95 0.5162 0.0000
pair_coeff 2 2 5670.33 0.5235 8.1163
neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes
group type1 type 1
compute charge1 type1 property/atom q
compute q1 type1 reduce ave c_charge1
group type2 type 2
compute charge2 type2 property/atom q
compute q2 type2 reduce ave c_charge2
variable qtot equal count(type1)*c_q1+count(type2)*c_q2
thermo_style custom step pe c_q1 c_q2 v_qtot spcpu
thermo 10
timestep 0.0001
velocity all create 300.0 1281937
fix 1 all nve
fix 2 all qeq/dynamic 1 10 1.0e-3 100 bi2o3.parm
run 100
write_data bi2o3_qeq_o.data
Here is the error:
WARNING: Charges did not converge at step 0: nan (…/fix_qeq_dynamic.cpp:165)
ERROR on proc 0: Non-numeric atom coords - simulation unstable (…/domain.cpp:521)
Best wishes
Shaw
Your charges did not converge so atoms probably ended up with spurious charge values, causing Coulomb potential to misbehave.
-
Start with more meaningful initial charges, i.e., charges that are close to equilibrium charges or just some fractional values instead of zeros would help a lot.
-
Check your QEq parameters in “bi2o3.parm” as they might not make sense.
Dear all,
I was trying to develop charges on a crystal by using Qeq. I found the
ionization potential and electron affinity parameters to get the Mulliken
electronegativity, which is used by Qeq. But it gave NaN on charges and
those energies. The crystal structure is from the standard database. It
seems there is something wrong blowing up the system. Can anyone know how
to deal with it? I was using the same method successfully once before, but
I have no idea why this time it is not working.
what happens, if you run that input withOUT fix qeq and just assign fixed
charges to some small initial guess value? do you still get NaNs?
axel.
I removed the Qeq part. But it still cannot calculate the potential energy and left as -nan in the output.
I removed the Qeq part. But it still cannot calculate the potential energy
and left as -nan in the output.
that means, that most likely your input geometry in the data file is
bogus, e.g. due to atoms overlapping atoms.
axel.
I have used the topo to check with the data file. But I did not find overlapping in the geometry. What I am wondering is that no matter how bad the potential is, there should be an energy output. But here, only nan left.
The following is the data file.
83 atoms
2 atom types
-1.363 6.813 xlo xhi
-1.363 6.813 ylo yhi
-1.363 6.813 zlo zhi
Masses
1 209
2 16
Atoms
1 2 -2 1.363 4.088 4.088
2 1 2.486486486 0 0 0
3 2 -2 1.363 4.088 6.813
4 1 2.486486486 5.45 2.725 2.725
5 2 -2 1.363 -1.363 4.088
6 1 2.486486486 2.725 2.725 0
7 2 -2 6.813 4.088 6.813
8 1 2.486486486 2.725 5.45 2.725
9 2 -2 1.363 -1.363 -1.363
10 1 2.486486486 0 2.725 2.725
11 2 -2 1.363 4.088 1.363
12 1 2.486486486 5.45 5.45 5.45
13 2 -2 1.363 4.088 -1.363
14 1 2.486486486 2.725 0 2.725
15 2 -2 6.813 4.088 1.363
16 1 2.486486486 2.725 2.725 5.45
17 2 -2 4.088 1.363 4.088
18 2 -2 6.813 1.363 4.088
19 2 -2 4.088 1.363 -1.363
20 2 -2 6.813 6.813 4.088
21 2 -2 -1.363 1.363 -1.363
22 2 -2 1.363 1.363 4.088
23 2 -2 -1.363 1.363 4.088
24 2 -2 1.363 6.813 4.088
25 2 -2 4.088 4.088 1.363
26 2 -2 4.088 6.813 1.363
27 2 -2 -1.363 4.088 1.363
28 2 -2 4.088 6.813 6.813
29 2 -2 -1.363 -1.363 1.363
30 2 -2 4.088 1.363 1.363
31 2 -2 4.088 -1.363 1.363
32 2 -2 4.088 1.363 6.813
33 1 2.486486486 2.725 2.725 2.725
34 2 -2 1.363 4.088 4.088
35 1 2.486486486 2.725 5.45 0
36 2 -2 1.363 4.088 1.363
37 1 2.486486486 0 5.45 2.725
38 2 -2 1.363 4.088 4.088
39 1 2.486486486 0 2.725 0
40 2 -2 1.363 4.088 1.363
41 1 2.486486486 2.725 5.45 5.45
42 2 -2 1.363 4.088 4.088
43 1 2.486486486 2.725 2.725 2.725
44 2 -2 1.363 4.088 1.363
45 1 2.486486486 0 2.725 5.45
46 2 -2 1.363 4.088 4.088
47 1 2.486486486 0 5.45 2.725
48 2 -2 1.363 4.088 1.363
49 1 2.486486486 2.725 2.725 2.725
50 2 -2 4.088 1.363 4.088
51 1 2.486486486 0 2.725 5.45
52 2 -2 1.363 1.363 4.088
53 1 2.486486486 2.725 0 5.45
54 2 -2 4.088 1.363 4.088
55 1 2.486486486 0 0 2.725
56 2 -2 1.363 1.363 4.088
57 1 2.486486486 5.45 2.725 5.45
58 2 -2 4.088 1.363 4.088
59 1 2.486486486 2.725 2.725 2.725
60 2 -2 1.363 1.363 4.088
61 1 2.486486486 5.45 0 2.725
62 2 -2 4.088 1.363 4.088
63 1 2.486486486 2.725 0 5.45
64 2 -2 1.363 1.363 4.088
65 1 2.486486486 2.725 2.725 2.725
66 2 -2 4.088 4.088 1.363
67 1 2.486486486 5.45 0 2.725
68 2 -2 4.088 1.363 1.363
69 1 2.486486486 5.45 2.725 0
70 1 2.486486486 2.725 0 0
71 1 2.486486486 5.45 5.45 2.725
72 1 2.486486486 2.725 5.45 0
73 1 2.486486486 2.725 0 5.45
74 2 -2 4.088 1.363 4.088
75 1 2.486486486 5.45 2.725 5.45
76 2 -2 1.363 4.088 1.363
77 1 2.486486486 0 2.725 0
78 1 2.486486486 0 5.45 2.725
79 2 -2 1.363 1.363 4.088
80 2 -2 1.363 4.088 4.088
81 1 2.486486486 0 2.725 5.45
82 1 2.486486486 2.725 5.45 5.45
83 1 2.486486486 0 0 2.725
I have used the topo to check with the data file. But I did not find
overlapping in the geometry. What I am wondering is that no matter how bad
the potential is, there should be an energy output. But here, only nan left.
look at the number of entries: you have 83 atoms, but your compound is
Bi_2O_3 (i.e. 5 atoms per units). that *cannot* be correct, when assuming
that you want to model a bulk crystal with no defects. also, have you
checked overlaps through periodic boundaries?
this simple input will show you how many overlaps there are:
units real
atom_style charge
boundary p p p
read_data bi2o3.data
pair_style zero 5.0
pair_coeff * *
delete_atoms overlap 0.01 all all
write_data bi2o3-new.data
if you look at the output, you'll see:
Deleted 44 atoms, new total = 39
that means the delete_atoms command identifies not just a few but 44(!)
atoms that are within less than 0.01 angstrom of each other.
it also looks as if one position is missing, since the remaining 39 atoms
is as unreasonable as 83 (not divisible by 5) for a bulk system.
axel.
Thanks for your replying. Now I think I might check with the cif file which I used to build up a crystal. There might be a mistake in the original cif file so the output structure might be wrong.