I am currently working on modifying the “fix_widom” source code, and I am encountering an issue related to forces update. Specifically, I want to calculate the difference between the forces before and after insertion of an atom. Here is the relevant code snippet from “fix_widom.cpp”:
void FixWidom::attempt_atomic_insertion_full() {
double ** previousF = atom->f;
// insertion of an atom
// …
// update forces, including the inserted atom
force->setup();
double ** nowF = atom->f;
// compute force difference between before and after insertion
for (int i = 0; i < natoms; i++) {
resultF[i][0] = nowF[i][0] - previousF[i][0];
resultF[i][1] = nowF[i][1] - previousF[i][1];
resultF[i][2] = nowF[i][2] - previousF[i][2];
}
}
However, I have noticed that the forces cannot be changed back after using “force->setup()”. I am unsure whether this will affect the subsequent simulation steps, such as velocity update.
Could someone please advise me on how to properly calculate the force difference while ensuring the correct functioning of the simulation?
This is mostly a question of understanding of C++ programming and data structures.
You have made a copy of the pointer to the storage of the forces, not of the actual force data.
Thank you for your reply. But I may not have made myself clear. When I invoke the “force->setup()”, I think I will change the forces of the system. After I get the “resultF”, I want the “atom->f” to recover to “previousF”. For example, there are 10 atoms, and the I insetion an atom to get the 10 atoms’ forces difference. After that the simulation is still about the 10 atoms. So, the atom->f must be “previousF”.
There are multiple issues with the piece of code you quoted and what you are saying, but - I repeat - what you are copying is a pointer and not the data itself. That is a fundamental problem.
If you look at the source code of the Force class, the call to Force::setup() is effectively a NOP.
To reset energy and forces, you have to do all the steps that are done in the FixWidom::energy_full() function plus the reverse communication (see the comments there).
In addition, you have to keep in mind that you are reneighboring, which means that the local order of atoms can change, so you must not compare by index directly, but have to loop over atom IDs and also make a copy of the atom->map() data from the before state, to be able to compare same with same.
“My understanding of what you said about deep copy and shallow copy is that changing the force in the fix_widom instance does not change the force in other instances. I have another question: why is it that when the atomic coordinates are the same, the forces calculated by energy_full() in the fix_widom.cpp file are different from the “actual” forces (where “actual” forces refer to the forces outputted by dump)?”
I don’t see why they would be different. You have to back up claims like this with hard evidence.
Given that you have now repeatedly displayed difficulties grasping rather simple C/C++ data structure concepts, it is not too much of a stretch to assume that you are mixing up things. Like the fact that indexing in local arrays and dump output (by atom->tag) is different.
In the attempt_atomic_insertion_full() function of fix_widom , I output the coordinates and forces before and after insertion an atom. I test and I probably knew that the forces hasn’t been updated when calling “fix_widom”, but the “atom->x” has been updated. It is said that the “atom->f” correspond to last frame “atom->x”. So when execute the statement “energy_stored = energy_full()”, the “atom->f” updates. After that the “atom-f” correspond to the“atom->x". But in my out.lammpstrj and log file, the forces are different.
I am so sorry if I caused any confusion in our communication as I am a beginner in C++. Also, I want to clarify that all my tests were conducted without using MPI, i.e., by running the command ‘lmp_mpi -i in.xxx > log’. testForce.tgz (234.2 KB)