storing forces on adjacent beads in fix_pimd

LAMMPS LAMMPS Mailing List Mirror
1

Dear Users,

I am trying to modify the ‘post_force’ method of fix_pimd.cpp in USER-MISC (only for the case method=PIMD).
What I need is to store the forces that are acting on the atoms of the two adjacent replicas in such a way that they will be available within method ‘spring_force’, exactly as already implemented for the positions using “buf_beads”.

What is the best way of doing it?

I tried following what is already implemented for the positions, and did the two following main attempts.

(1)
I defined a new method ‘void mycomm_exec(double **, double **)’, which receives as second argument a pointer to the stored quantity.
I implemented it as ‘void FixMYPIMD::mycomm_exec(double **ptr, double **buffer)’, just a copy and paste of the original method ‘comm_exec’, changing everywhere “buf_beads” with “buffer”.
I defined a new public pointer “double **buf_force” and initialized following exactly what is done for “buf_beads”, both in the constructor and in method ‘comm_init’.
Finally, in method ‘post_force’ I call twice:

mycomm_exec(atom->x,buf_beads); //Store pos

mycomm_exec(atom->f,buf_force); //Store forces

Doing a single call, either for ‘f’ or ‘x’ gives the correct result, but when calling twice as above I get segmentation fault with “Signal code: Address not mapped”.

(2)

I also tried using the original ‘comm_exec’ twice and memcpy:

comm_exec(atom->f); //Store forces in buf_beads
memcpy(buf_force,buf_beads,sizeof(double)atom->nlocal3*size_plan); //Copy to buf_force
comm_exec(atom->x); //Store pos in buf_beads

In this case I noticed that “buf_beads” and “buf_force” end up pointing to the same block of memory.

I thought it would have been a straightforward modification following what already implemented, but at this point I think I am missing something fundamental, most probably about pointers.
I’m quite new to MPI and c++, any hint or suggestions will be of great help.

Thanks!

Best,
Davide

2

I looked briefly at fix_pimd, but don’t understand its data structures or exactly
what the comm_exec() is doing. I will say that this line in your (2)

memcpy(buf_force,buf_beads,sizeof(double)atom->nlocal3*size_plan); //Copy to buf_force

cannot work b/c buf_force and buf_beads are double ** pointers.
They appear to be allocated in 2 steps as np, then as a series of 3*nlocal mallocs
of doubles. So you would have to copy those Np chunks one at a time
with memcpy. Memcpy will only work with a contiguous chunk of memory.

I’m not sure why they didn’t allocate buf_beads with memory->create() as
a 2d or even 3d (np by nlocal by 3) array. That would guarantee that
the low-level memory is contiguous. And you could do the copy with
a single memcpy, e.g.

memcpy(&yourbuf[0][0][0],&buf_beads[0][0][0],npnlocal3*sizeof(double));

Also, both the authors (top of file) are at ALCF (affiliated with Argonne),
so you could certainly contact them to discuss your model and changes
to their code.

Steve

3

Thank you very much Steve, for your reply. Earlier today I had actually solved this problem doing exactly what you just suggested about copying each chunk one at time. Now it’s more clear, though. Also I have contacted one of the author and will discuss with him in the future.
Best,
Davide