storing forces on adjacent beads in fix_pimd

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.

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”.


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.



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.


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.


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.