A bonded network in semi-priodic system (y-fixed and x-periodic)

Dear all,
I am using LAMMPS to study a bonded network which is semi-periodic. My x-directions are periodic, and y-directions are fixed, and no atoms form any bond across the y-boundaries. However, I have bonds that do cross the x-boundaries. And I am using a simple harmonic bond style. I simulate in the microcanonical ensemble, i.e. NVE, from this simple unstressed/unperturbed network where the distance between two atoms is equal to their equilibrium length, so no force is acting between any two atoms, and no potential energy is expected.
However, I am observing a shrinkage in the x-direction due to all the bonds that crossed the boundaries.

Could anyone please help me?
LAMMPS version use: LAMMPS (29 Aug 2024 - Update 2)
The input file, data file and log files are attached. Also, the images of before and after are attached.
in.simple (2.2 KB)
log.lammps (15.0 KB)
STATE-Z5.0771 (333.5 KB)
t0
t1000

That’s a lot of bond types. Are you certain the bond lengths are calculated correctly across the x periodic boundary?

Alternatively, you could swap a BPM bond style which automatically sets the equilibrium length to the initial length: 8.5.4. Bonded particle models — LAMMPS documentation. Then you only need one bond type.

Dear @jtclemm ,
Thank you so much for your prompt response. Yes, it did work for the BPM bond style.

This led me to a follow-up question: Is the harmonic bond style not designed to incorporate the semi-periodic systems?

Once again, thank you @jtclemm.
Cheers,
Harish

Those two bond styles represent different models, so you are comparing apples to oranges.

Rather than choosing a bond style based on whether it results in the behavior you expect, you have to answer the question of what should the model be and then decide on whether your expectation has to be revised or whether there might be a problem with your (initial) geometry or a problem with the topology.

All bonds in LAMMPS work across periodic boundary conditions. I just guessed your most likely error was defining the correct coefficients since there’s a lot of room for mistakes (did you correctly remap across the periodic box in your preprocessing?). The bpm bond styles, like bpm/spring, are just useful for stress free networks since they automatically calculate equilibrium lengths and remove this potential source of error.

Actually, bpm/spring has the same functional form as harmonic if you turn off damping. It’s just a matter of what is used for an equilibrium length (and maybe a factor of 2 in the stiffness).

That still makes it a different model in my book. The equilibrium distance is part of the model.
That said, there is also a variant of bond style harmonic that takes the initial distance as input in bond style harmonic/restrain.

What if you can exactly reproduce the behavior of bpm/spring with harmonic and harmonic with bpm/spring? :wink: You could say they only differ in how you input an equilibrium length.

I have seen that bond style before. If I recall correctly, instead of saving a reference length it stores initial atom coordinates and recalculates the reference length (which creates bugs if your simulation geometry ever changes). Would be interesting to benchmark against bpm/spring.

Thank you @jtclemm and @akohlmey for your replies.
Here is my presence for the entire conversation:
(1) Yes, they are two different types of bonds. In “bond stype bpm/spring”, equilibrium bond lengths (ro) are at the time = 0, and all the particles are of point size. Which is “same” as defining “bond style harmonic” and in “bond coeff” provide “ro” by calculating it at time = 0 based on the initial topology.

(2) I have taken care of bond lengths and wraparound effects carefully. Yes, the damping off and stiffness factor 2 comparison is giving me the expected results.

(3) I have checked with " harmonic/restrain" with only “bonid atom1 atom2” in the data file, and it gives me the SAME results as “harmonic”

(4) So, theoretically, “bond style harmonic/restrain” and “bond style bpm/spring” should give identical results except for the factor of “2” in the force calculation. This is not happening at the moment. It is only giving identical results for the “double periodic” system and not for the “semi-periodic” systems in my case.

I am sorry if I am missing something here.

Thank you again. :slight_smile:

Yes, although the bpm bond style has a lot of optional features which will change results if turned on. For instance, you probably want break no smooth no unless you want to model fracture.

Looking at your data file, it looks like this may not be exactly the case for your harmonic implementation. E.g. bond #2150 has an equilibrium length of 140 which does not reflect the distance to the nearest periodic image (~2).

Given periodic issues in your data file for harmonic, I’m surprised harmonic/restrain gives similar errors. I assumed it correct PBCs, but I can’t say what might be wrong without looking into it further. Maybe there’s something with image flags.

Theoretically yes, but (as stated above) there are some conditions.

I’m biased, but personally I think the bpm/spring bond style may be the most appropriate for your system since it was designed to do exactly what you are trying to accomplish.

Hope this helps.