Custom fix for bond orientation leads to center-of-mass drift in parallel runs

LAMMPS LAMMPS Development
1

Dear LAMMPS developer:

I am currently working on repeat and extend the procedure of a paper in polymer crystallization as attached. The idea of the paper is to apply an external orientational force on bonds inside a given region to a specific direction, which intends to accelerate the regular packing of the polymer chain. The external potential given in the original paper is:

E = k * (theta)^2 for theta < 0.5π
E = k * (pi - theta)^2 for theta ≥ 0.5π

Where theta is the angle between the bond and the orientation direction.

To the best of my knowledge, LAMMPS does not currently offer a built-in fix with this functionality (the closest might be angle_style dipole, but that is molecule-specific). Therefore, I have written a custom fix (based on version 29Aug2024 Update 1, compiles without errors). With respect to the original paper and to avoid the costly computation of arccos, I use the following equation to represent the potential energy:

E = k*(sin(theta-theta0))^2

Where theta is still the angle between the bond and the orientation direction. Theta0 is the equilibrium angle between the chemical bond and the orientation direction. The reason theta0 was added because the original paper was simulating a polymer chain with equilibrium bond angle at 180 degree.

The syntax of this fix is simple:

fix ID group-ID region “type” bondtype k theta0 rx ry rz “type” bondtype k theta0 rx ry rz …

  • region : region in which the bond midpoint must lie to apply the force.
  • bondtype : bond type to which this potential applies.
  • k : force constant.
  • theta0 : reference angle (in degrees) between the bond and the crystal direction vector (rx,ry,rz).

When testing this fix on a simple OPLS-UA HDPE system, I could observed “flying-ice” draft (center-of-mass velocity develops, not the random thermal motion) right after the system was compacted from extreme low density state. I suspect this may be due to a parallel communication issue that I did not currently correctly set. The simulation setting, including in, data, potential file, was attached. The drafting is significant after timestep 40000.

I also tried (and attached) some simple two atoms with one bond (point from 1 to 2) system in size of 20A*20A*20A, with explicit command “newton on on”:

  • case 1: running with -np 24, no bond length constrain, atom 2 lost force when bond length > 3.65 at step 146, already set comm_modify cutoff 5.0
  • case 2: running with -np 24, no bond length constrain, atom 2 lost force when bond length > ~13, which probably due to bond length > half of
  • cell case 3: running with -np 24, has bond length constrain, significant “flying-ice” draft after step 6000
  • case 3_nofix: running with -np 24, has bond length constrain, new fix disabled, no “flying-ice” draft.

I am not sure where I did wrong in the fix’s coding or the setting of the simulation. I would appreciate any guidance.

Best regards

Jiawei Zhao [email protected]
SINOPEC (Beijing) Research Institute of Chemical Industry

original paper :
deformation-and-fracture-processes-of-a-lamellar-structure-in-polyethylene-at-the-molecular-level-by-a-coarse-grained.pdf (6.5 MB)

fix’s src:
fix_orient_crystal.h (938 Bytes)
fix_orient_crystal.cpp (9.0 KB)

small trial systems and results:
case_3_nofix.zip (874.0 KB)
case_3.zip (992.3 KB)
case_2.zip (920.1 KB)
case_1.zip (944.9 KB)

polymer simulation files:
polymer_simulation.zip (483.2 KB)

2

Please understand that your chances that somebody will read the paper, then read and test your source code, figure out what the problem is, and then explain it to you are extremely slim. What you are asking for is a lot of work. In fact, I would claim that I would not even trust the findings of somebody that does this, because a person that would voluntarily do somebody else’s work without getting some kind of compensation for it does not seem like somebody making smart choices.

So basically, you need to recognize that any post of the kind “here is my code and what I want it to do, now please tell me what I am doing wrong” is not likely to get any kind of response except for the kind of what you are getting from me now (if any at all). :disappointed_face:

Your chances will be much better, if you ask very specific questions (or you need to figure out a way to hire somebody with the debugging skills you need to resolve your problem), which means that you have to do more debugging and check out in detail the data you get from your computation versus what you can manually compute. In the simplest approach that means adding a lot of printf() statements to your code and compare the output against your expectations.

If you are desperate, you might try using AI coding tools to debug your issue. From my personal experience, I would not expect the AI to actually suggest a correct solution (in most of the cases where I tried using an AI to resolve a problem, it would eliminate the problem rather by disabling the added functionality), it may provide some useful hints for where to look for the issue (it did in most of my attempts) and thus help to understand it. Understanding a problem is usually the most important step in finding a solution.

3

Dear Axel,
Thanks for the reply.

4

Here is a specific problem I would like to asked:

1st
I just found I have a typo in the previous post, which should be:

case 1: running with -np 24, no bond length constrain, atom 2’s force become 0 when bond length > 3.65 at step 146, already set comm_modify cutoff 5.0
case 2: running with -np 1, no bond length constrain, atom 2’s force become 0 when bond length > ~13, which probably due to bond length > half of cell

Those are two identical simulations run by serial and parallel. In my code, what I did (or believe have done) is for loop the bond in each processor at post_force, for each bond do the logic and apply force for both atoms, then relay on the reverse communication to add the force for ghost atoms if needed, I wonder why the reverse communication failed in the parallel run much earlier then the serial one, which may be the cause.

Jiawei Zhao

5

The “automatic” reverse communication has already happened when Fix::post_force() is called. You would have to use Fix::pre_reverse() instead. See 4.7. How a timestep works — LAMMPS documentation

If you want to do your own reverse communication, you have to either first clear the forces on ghost atoms, since the reverse communication does not do that, or allocate your own array of size atom->nmax and then clear it and add forces to it and then call your own reverse communication, for which you need to implement the corresponding pack/unpack functions, and finally add the collected forces on local atoms to the force arrays of the atoms.

2 Likes
6

Thanks! I relocate the main logic from Fix::post_force() to Fix::pre_reverse() instead. It make sense now and the behavior is what I expected.

2 Likes
7

@apm0zjw Thanks for reporting back. This will hopefully help others that have a similar problem and are searching for a solution.

2 Likes