I introduced a new member to atom_vec_full called double **efield. efield holds the value of the electric field at points in space where atoms are located, thus efield looks like
efield[nlocal][3 coordinates]. The values for efield are taken from the force calculation in pair_lj_cut_coul_long.cpp and are assigned in the compute() function. (I modeled efield on the existing double **f, using the same function calls, such as memory->grow_2d_double_array(), etc.)
If I do a simulation where I print out the efield values for a particular atom at each time step, then I get two different results depending on the input script.
Inputs 1 and 2 generate different values for efield for timesteps larger than 10 timesteps. But the thermodynamic quantities from thermo are the same.
I've gone over the code in write_restart.cpp and output.cpp but I don't see how this makes sense. Is there some part of atom_vec_full.cpp that is affected by writing restarts? Although the differences after 10 timesteps are minor, the differences grow with time until one of the systems blow up due to poor energy conservation.
Any help on this would be great! Thanks in advance.