[lammps-users] Possible error in velocities of ghost atoms in dpd style

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

I’m doing a calculation which requires the velocities of ghost atoms.
Following Steve Plimpton’s suggestion I’m using atom_style dpd, but
my tests appear to show that the velocities of the ghost atoms are
not identical to those of the parent atoms.

I attach the compute, input and output files which test this. For
simplicity I consider a single FCC cell with lattice parameter 5.0, and
only print out the first few ghosts. At the zeroth step I get the following
results and everything appears to be fine:

atom 0 x= 0.0000000000 y= 0.0000000000 z= 0.0000000000
vx= 0.0004940417 vy= 0.0009573655 vz= 0.0013279029
atom 1 x= 2.5000000000 y= 2.5000000000 z= 0.0000000000
vx= -0.0006975682 vy= 0.0009922590 vz= -0.0000905373
atom 2 x= 2.5000000000 y= 0.0000000000 z= 2.5000000000
vx= -0.0004465096 vy= -0.0010397911 vz= 0.0002940638
atom 3 x= 0.0000000000 y= 2.5000000000 z= 2.5000000000
vx= 0.0006500360 vy= -0.0009098333 vz= -0.0015314294
atom 4 x= 5.0000000000 y= 0.0000000000 z= 0.0000000000
vx= 0.0004940417 vy= 0.0009573655 vz= 0.0013279029
atom 5 x= 7.5000000000 y= 2.5000000000 z= 0.0000000000
vx= -0.0006975682 vy= 0.0009922590 vz= -0.0000905373
atom 6 x= 7.5000000000 y= 0.0000000000 z= 2.5000000000
vx= -0.0004465096 vy= -0.0010397911 vz= 0.0002940638
atom 7 x= 5.0000000000 y= 2.5000000000 z= 2.5000000000
vx= 0.0006500360 vy= -0.0009098333 vz= -0.0015314294
atom 8 x= -5.0000000000 y= 0.0000000000 z= 0.0000000000
vx= 0.0004940417 vy= 0.0009573655 vz= 0.0013279029
atom 9 x= -2.5000000000 y= 2.5000000000 z= 0.0000000000
vx= -0.0006975682 vy= 0.0009922590 vz= -0.0000905373
atom 10 x= -2.5000000000 y= 0.0000000000 z= 2.5000000000
vx= -0.0004465096 vy= -0.0010397911 vz= 0.0002940638
atom 11 x= -5.0000000000 y= 2.5000000000 z= 2.5000000000
vx= 0.0006500360 vy= -0.0009098333 vz= -0.0015314294

From this it’s clear that
4 and 8 are ghosts of atom 0,
5 and 9 are ghosts of atom 1,
6 and 10 are ghosts of atom 2,
7 and 11 are ghosts of atom 3,
and the velocities of the ghosts are identical to those of their parents.

However, after the first step the velocities of the ghosts differ from those
of their parents, and the disagreement gets worse with every timestep.
After just one iteration I obtain:
atom 0 x= 0.0019761668 y= 0.0038294618 z= 0.0053116117
vx= 0.0004923162 vy= 0.0009560366 vz= 0.0013248753
atom 1 x= 2.4972097272 y= 2.5039690359 z= -0.0003621493
vx= -0.0006955450 vy= 0.0009907359 vz= -0.0000914899
atom 2 x= 2.4982139618 y= -0.0041591645 z= 2.5011762552
vx= -0.0004448536 vy= -0.0010381985 vz= 0.0002947187
atom 3 x= 0.0026001442 y= 2.4963606668 z= 2.4938742824
vx= 0.0006480824 vy= -0.0009085740 vz= -0.0015281041
atom 4 x= 5.0019761668 y= 0.0038294618 z= 0.0053116117
vx= 0.0004940417 vy= 0.0009573655 vz= 0.0013279029
atom 5 x= 7.4972097272 y= 2.5039690359 z= -0.0003621493
vx= -0.0006975682 vy= 0.0009922590 vz= -0.0000905373
atom 6 x= 7.4982139618 y= -0.0041591645 z= 2.5011762552
vx= -0.0004465096 vy= -0.0010397911 vz= 0.0002940638
atom 7 x= 5.0026001442 y= 2.4963606668 z= 2.4938742824
vx= 0.0006500360 vy= -0.0009098333 vz= -0.0015314294
atom 8 x= -4.9980238332 y= 0.0038294618 z= 0.0053116117
vx= 0.0004940417 vy= 0.0009573655 vz= 0.0013279029
atom 9 x= -2.5027902728 y= 2.5039690359 z= -0.0003621493
vx= -0.0006975682 vy= 0.0009922590 vz= -0.0000905373
atom 10 x= -2.5017860382 y= -0.0041591645 z= 2.5011762552
vx= -0.0004465096 vy= -0.0010397911 vz= 0.0002940638
atom 11 x= -4.9973998558 y= 2.4963606668 z= 2.4938742824
vx= 0.0006500360 vy= -0.0009098333 vz= -0.0015314294

Is this a bug? Is there a mistake in my test? Or am I simply
misunderstanding something? I’d be very grateful for any suggestions.

Comments for anyone who looks at my program in detail:
In the compute I calculate the total KE and PE and output them
together with the KE and PE as calculated by LAMMPs; these agree
with each other. At every step I also write x and v for the
first 12 atoms to the screen (a copy of the screen output is attached
in check_vel.out). The velocities of the real atoms (0,1,2,3) agree with
those dumped by LAMMPS.

Thanks in advance,

Philip

<<check_vel.out>> <<check_vel.param>> <<compute_heatflux3.cpp>> <<compute_heatflux3.h>>

check_vel.out (9.6 KB)

check_vel.param (694 Bytes)

compute_heatflux3.cpp (5.88 KB)

compute_heatflux3.h (1.06 KB)

2

When in the timestep are you printing the v of ghost atoms? If
it's at the end, after final integration, then they won't be updated.
They are communicated in the middle of the step, so they are
correct for using them in force computations.

Steve

3

Steve,

I print them out at the end of the timestep, so now I understand the discrepancy. However, this suggests that I have a problem with the parent atoms when I evaluate the scalar product v.F, since the velocities and the forces are separated by half a timestep. Of course I can (and will) check how my final results depend on the size of the timestep, which should indicate whether this time shift is significant. Is there any (easy) way to access the velocities of the parent atoms in the middle of the integration, so that I have them at the same time as the forces and as the velocities of the ghosts?

Thanks for your help,

Philip

4

The sequence of the velocity Verlet timestep is:

update v by 1/2 step
update x by full step
communicate x (and v if atom_style dpd)
compute forces
update v by 1/2 step

If your fix is a diagnostic (end of step), and you need ghost v, then
you have no choice but to communicate them. See other fixes
that communicate quantities via comm.cpp, e.g. fix shake

Steve