nan in fix/qeq

I used also the parameters in the folder you specified, but the problem insists. When you run the simulation, it ends without error but when you check your trajectory file, the charges are nan. Moreover, the E_coul term that thermo writes is zero. Have you check the charges of each atom after the run? If yes, what are they?

And a segmentation fault error arises if you run with ewald. Have you turned on ewald in your run?

My traj file is

ITEM: TIMESTEP
0
ITEM: NUMBER OF ATOMS
3
ITEM: BOX BOUNDS pp pp pp
-10 10
-10 10
-10 10
ITEM: ATOMS id type q xs ys zs ix iy iz vx vy vz
1 1 nan 0.5 0.5 0.5 0 0 0 0 0 0
2 2 nan 0.544 0.5 0.5 0 0 0 0 0 0
3 1 nan 0.566 0.538105 0.5 0 0 0 0 0 0
ITEM: TIMESTEP
1
ITEM: NUMBER OF ATOMS
3
ITEM: BOX BOUNDS pp pp pp
-10 10
-10 10
-10 10
ITEM: ATOMS id type q xs ys zs ix iy iz vx vy vz
1 1 nan 0.5 0.5 0.5 0 0 0 0 0 0
2 2 nan 0.544 0.5 0.5 0 0 0 0 0 0
3 1 nan 0.566 0.538105 0.5 0 0 0 0 0 0

baris.

I used also the parameters in the folder you specified, but the problem insists. When you run the simulation, it ends without error but when you check your trajectory file, the charges are nan. Moreover, the E_coul term that thermo writes is zero. Have you check the charges of each atom after the run? If yes, what are they?

And a segmentation fault error arises if you run with ewald. Have you turned on ewald in your run?

My traj file is

ITEM: TIMESTEP
0
ITEM: NUMBER OF ATOMS
3
ITEM: BOX BOUNDS pp pp pp
-10 10
-10 10
-10 10
ITEM: ATOMS id type q xs ys zs ix iy iz vx vy vz
1 1 nan 0.5 0.5 0.5 0 0 0 0 0 0
2 2 nan 0.544 0.5 0.5 0 0 0 0 0 0
3 1 nan 0.566 0.538105 0.5 0 0 0 0 0 0
ITEM: TIMESTEP
1
ITEM: NUMBER OF ATOMS
3
ITEM: BOX BOUNDS pp pp pp
-10 10
-10 10
-10 10
ITEM: ATOMS id type q xs ys zs ix iy iz vx vy vz
1 1 nan 0.5 0.5 0.5 0 0 0 0 0 0
2 2 nan 0.544 0.5 0.5 0 0 0 0 0 0
3 1 nan 0.566 0.538105 0.5 0 0 0 0 0 0

baris.

For some reason your system/setup does work with the JPC parameters, but it does work with the QEq parameters from the examples/reax/CHO/param.qeq file.

This is not a bug but related to the sensitivity of the current QEq implementation to its parameters.

Ray

I used also the parameters in the folder you specified, but the problem insists. When you run the simulation, it ends without error but when you check your trajectory file, the charges are nan. Moreover, the E_coul term that thermo writes is zero. Have you check the charges of each atom after the run? If yes, what are they?

Of course. q_H = 0.0625, q_O = -0.125, E_coul = 0.0. But your E_coul is also zero without fix qeq/reax. You might want to check your pair_style, bond_style, angle_style commands.

What version of LAMMPS are you using? Did you use $LMP_DIR/examples/reax/CHO/param.qeq as is? Or did you remove the second line?

And a segmentation fault error arises if you run with ewald. Have you turned on ewald in your run?

Of course. It is not the problem with Ewald, but with bond_style and angle_style. If you remove those and use a pure pair_style, charges are more reasonable with both coul/cut and coul/long.

Ray

I tried both the latest stable and development versions of lammps.

I deleted the line for Carbon, so two lines remained one for H the other for O.

Did you turn off the bonded terms in the in file?

My E_coul term is zero, because I did not assign initial charges to the atoms. Did you assign them?

But Ray, if you use pure pair_style (without bonded terms) how can you evaluate the charge distribution with time (I mean, for example the use of qeq/reax with NVT )?

baris.

I used also the parameters in the folder you specified, but the problem insists. When you run the simulation, it ends without error but when you check your trajectory file, the charges are nan. Moreover, the E_coul term that thermo writes is zero. Have you check the charges of each atom after the run? If yes, what are they?

Of course. q_H = 0.0625, q_O = -0.125, E_coul = 0.0. But your E_coul is also zero without fix qeq/reax. You might want to check your pair_style, bond_style, angle_style commands.

What version of LAMMPS are you using? Did you use $LMP_DIR/examples/reax/CHO/param.qeq as is? Or did you remove the second line?

And a segmentation fault error arises if you run with ewald. Have you turned on ewald in your run?

Of course. It is not the problem with Ewald, but with bond_style and angle_style. If you remove those and use a pure pair_style, charges are more reasonable with both coul/cut and coul/long.

Ray

I tried both the latest stable and development versions of lammps.

I deleted the line for Carbon, so two lines remained one for H the other for O.

Did you turn off the bonded terms in the in file?

Yes, and removed bond_style settings in the structure file.

My E_coul term is zero, because I did not assign initial charges to the atoms. Did you assign them?

This does not make sense. Your E_coul is zero with and without charges. You of course assigned initial charges to the atoms: 0.4 and -0.8 for H and O, respectively.

But Ray, if you use pure pair_style (without bonded terms) how can you evaluate the charge distribution with time (I mean, for example the use of qeq/reax with NVT )?

This does not make sense either. Fix qeq/reax is originally designed for ReaxFF, which is a pair_style without pre-defined bonds and angles.

So, as far as I understand

  1. fix qeq/reax should be used only with pure pair_style (no bonded styles in in.lammps file)
  2. charge evolution with time is not possible with fix qeq/reax command (I mean use fix qeq/reax with NVT), because nonbonded terms deleted.

Am I right?

baris.

I tried both the latest stable and development versions of lammps.

I deleted the line for Carbon, so two lines remained one for H the other for O.

Did you turn off the bonded terms in the in file?

Yes, and removed bond_style settings in the structure file.

My E_coul term is zero, because I did not assign initial charges to the atoms. Did you assign them?

This does not make sense. Your E_coul is zero with and without charges. You of course assigned initial charges to the atoms: 0.4 and -0.8 for H and O, respectively.

But Ray, if you use pure pair_style (without bonded terms) how can you evaluate the charge distribution with time (I mean, for example the use of qeq/reax with NVT )?

This does not make sense either. Fix qeq/reax is originally designed for ReaxFF, which is a pair_style without pre-defined bonds and angles.

So, as far as I understand

  1. fix qeq/reax should be used only with pure pair_style (no bonded styles in in.lammps file)

Yes.

  1. charge evolution with time is not possible with fix qeq/reax command (I mean use fix qeq/reax with NVT), because nonbonded terms deleted.

No. You can sure use fix qeq/reax with fix nvt and get charge evolution.

My excuses for spooning in with a short piece of advice. Baris, you should go back and read the original paper that covers the development of the QEq method. Your questions display a lack of understanding of the working principles of the method, that’s why your continuous back and forth asking Ray to pretty much guide you along your trial and error baby steps.

Carlos