Dear Lammps users,
I am working on coupling lammps with an external fluid solver. The immersed boundary interface was used to exchange the force information between fluid solver and lammps force **f.
The algorithm is like this:
- calculate interface force f_int
- apply f_int to fluid, and update fluid
- apply f_int to lammps through
double **f = atom->**f;
f[i][0] += f_int;
basically it just changes the particle force directly through my interface code.
4) finish one time step and loop again.
It works great if I use “fix 1 all nve” for deformable molecules. However, when I use “fix 1 all rigid single” for rigid molecules, The molecule doesn’t move at all.
I am wondering if particle force **f is initialized to zero or ignored somewhere in “fix rigid” command so that the step 3) is not working.
Any advice would be appreciated. Thank you very much for your help.
Dear Lammps users,
I am working on coupling lammps with an external fluid solver. The immersed
boundary interface was used to exchange the force information between fluid
solver and lammps force **f.
The algorithm is like this:
1) calculate interface force f_int
2) apply f_int to fluid, and update fluid
3) apply f_int to lammps through
double **f = atom->**f;
f[i][0] += f_int;
basically it just changes the particle force directly through my interface
code.
4) finish one time step and loop again.
It works great if I use "fix 1 all nve" for deformable molecules. However,
when I use "fix 1 all rigid single" for rigid molecules, The molecule
doesn't move at all.
I am wondering if particle force **f is initialized to zero or ignored
somewhere in "fix rigid" command so that the step 3) is not working.
no. fix rigid *does* use the atom->f array.
are you sure that you are using fix rigid correctly?
you are speaking of "molecules" (i.e. plural) yet are referring to
using the "single" option of fix rigid (i.e. one large rigid body over
all particles in the group). that doesn't quite fix together.
axel.
Dear Axel,
Thank you for your timely reply. Actually I treat all the atoms in the simulation as one rigid body. So I think "fix rigid single " is correct. Based on the description of the command “fix rigid”, it will sum over all the forces over all the atoms within that rigid body to get the resultant force and torque. And then it will be used to update atom velocity and position. That is exactly what I want. I also print the force value within fix_rigid.cpp, the output is either zero or extremely small(10^-20).
In my code, I did not write a new fix, instead, I use the lammps pointer to get the **f directly, and change it, e.g., double **f = lmp->atom->**f, It should work, as demonstrated with the case where “fix all nve” was used. I will double check it and let you know once I figure it out.
Thank you again for your reply.
Sincerely,
JT
I’m not clear on what you are doing. To apply externally computed forces
to the particles being time-integrated by LAMMPS, you need to
apply the forces at the correct point in the timestep. I don’t think
you can do that thru the Python interface.
What you can do is use fix external, which is designed to allow
an external code (C or Fortran) give external forces to LAMMPS so
that they are applied at the correct point in the timestep.
Steve
Dear Steve and Axel,
Thank you very much for the reply. I have figured it out. The problem is that during each run, the force vector will be cleared once. e.g., through force_clear() in verlet.cpp. This is also confirmed through testing
lammps_gather_atoms(lmp,“f”,1,3,f)
f[1][2] = epsilon;
lammps_scatter_atoms(lmp,“f”,1,3,f)
in simple.cpp in folder COUPLE/simple. No effect at all even f is changed. However, it works for position x and velocity v, as these two are continuous during each time step.
In my case, I followed “fix external” source code and created a new fix style. Instead of setting function pointer to callback, I directly change the public external force vector **fexternal in my external driver code through fix pointer static_cast conversion, as I need to access the public member of a child class. Then the fix will add **fexternal to **f through post_force(). Now the code runs successfully.
Hope it also helps others as well. Thank you very much!