Hi,
I am trying to find the forces acting on a given atom by all other
atoms of type k. To do that I have done the following things:
I defined a new variable per atom that is going to store this quantity
so I did the usual things, declare, create, destroy, zero it when the
force is zeroed.
atom.h
class Atom : protected Pointers { public:
double **ffromtype;
atom.cpp
Atom::~Atom()
memory->destroy_2d_double_array(ffromtype);
atom_vec_atomic.cpp
void AtomVecAtomic::grow(int n)
ffromtype = atom->ffromtype =
memory->grow_2d_double_array(atom->ffromtype,nmax,3,"atom:ffromtype");
verlet.cpp
void Verlet::force_clear(int vflag)
double **ffromtype = atom->ffromtype;
for (i = 0; i < nall; i++) {
ffromtype[i][0] = 0.0;
ffromtype[i][1] = 0.0;
ffomrtype[i][2] = 0.0;
}
I only have atoms that are LJ with a cut potential, so:
pair_lj_cut.cpp
void PairLJCut::compute(int eflag, int vflag)
right after the force is added:
ffromtype[i][0] += delx*fforce;
ffromtype[i][1] += dely*fforce;
ffromtype[i][2] += delz*fforce;
if (newton_pair || j < nlocal) {
ffromtype[j][0] -= delx*fforce;
ffromtype[j][1] -= dely*fforce;
ffromtype[j][2] -= delz*fforce;
for testing purposes first I want to get the same forces as the program
has in the atom->f so that I would know that I accounted for everything
that modifies the force in the program.
So I am writing out the forces acting o a given group of atoms (with
the custom fix sumforce) using the atom->f and the atom->ffromtype
variables and comparing them. I found that for a longer simulation
there are time intervals during which they are the same and intervals
when they are not. Meaning there is something in the program that
changes the forces that I am not aware of.
I have tried the following:
- modify all the fixes that act on my system to act on ffromtype as well
this did not change anything confirming that the fix wrote out the
forces before any other fixes acted on it (fixes act in the order specified in the inputfile ?)
- summing up ffromtype for vflag=2 where the atom->f is summed up (in
virial_compute), this did not change anything (as it should not)
Is there something else that could modify the forces that I am not
taking into account ?
Thanks,
Andras