Computing strain energy of crosslinked polymer networks

Hello.

I’m running Brownian dynamics simulations of semiflexible polymer chains in LAMMPS to analyze viscoelastic behavior. I’m applying periodic displacements via fix move to a group of atoms at the boundary of the simulation box to mimic rheological forcing. I’ve been trying to compute the strain energy of the entire network, but I’m running into some problems.

At first I was summing up the inner product of the forces/atom with the displacements/atom over the course of the simulation. But this incorporated the energy input from the fix langevin thermostat and lead to a linear rise in the computed quantity that dwarfs the input from the periodic applied strains. I have 3 concerns/questions:

  1. I’ve been trying to come up with a fix for this. I’d prefer to simply discard/subtract the the random force vector used by fix langevin from the total force on each atom. This requires modifying fix langevin to create a custom fix of my own that stores the computed per-atom force vectors for the drag and random forces in fix langevin, but I’m having trouble seeing how a LAMMPS fix interfaces with the dump command in the C++ classes to properly add a new function into a custom fix C++ class. Is there a clear example of a fix command that stores a dumpable per-atom vector I can model this off of? And, per the documentation section 11.7, should I be using the method, “compute_vector”? Can this get me per-atom vector?

  2. A second option is to recompute the drag force in fix langevin as a per atom vector variable in my LAMMPS script and track the individual, per-atom force contributions of every defined bond, angle, dihedral, and pair potential I have in my simulations, but somehow I think this would be more complicated and require more of an interface with the underlying C++ code.

  3. Is there an equivalent of the strain energy as a thermodynamic variable in LAMMPS that I’m missing? Something that captures the work done by bond forces, angle forces, dihedral forces, pair potential forces, as well as the viscous drag force from a Brownian dynamics model that I’m overlooking such as ‘pe’ (though I doubt pe stores any information about the viscous drag)?

I’ve using a post-processing work-around to use Fourier analysis to fit the transformed work computation to the transform of a linear ramp, and inverting the difference to obtain the waveform atop the monotonic thermal input, but something more direct like the custom fix mentioned in (1) above is more accurate and faithful to the data. I’ve been working on this problem for a few days, and the code is a little complex for me to adeptly modify it, so any help/tips is appreciated.

Hello.

I’m running Brownian dynamics simulations of semiflexible polymer chains in LAMMPS to analyze viscoelastic behavior. I’m applying periodic displacements via fix move to a group of atoms at the boundary of the simulation box to mimic rheological forcing. I’ve been trying to compute the strain energy of the entire network, but I’m running into some problems.

At first I was summing up the inner product of the forces/atom with the displacements/atom over the course of the simulation. But this incorporated the energy input from the fix langevin thermostat and lead to a linear rise in the computed quantity that dwarfs the input from the periodic applied strains. I have 3 concerns/questions:

  1. I’ve been trying to come up with a fix for this. I’d prefer to simply discard/subtract the the random force vector used by fix langevin from the total force on each atom. This requires modifying fix langevin to create a custom fix of my own that stores the computed per-atom force vectors for the drag and random forces in fix langevin, but I’m having trouble seeing how a LAMMPS fix interfaces with the dump command in the C++ classes to properly add a new function into a custom fix C++ class. Is there a clear example of a fix command that stores a dumpable per-atom vector I can model this off of? And, per the documentation section 11.7, should I be using the method, “compute_vector”? Can this get me per-atom vector?

no. compute_vector we provide a global vector. to have a fix store per-atom data that can be accessed from dumps and elsewhere, you need to do more steps:
there are a few flags that need to be set in the constructor. in your case you would have to have

int peratom_flag; // 0/1 if per-atom data is stored
int size_peratom_cols; // 0 = vector, N = columns in peratom array
int peratom_freq; // frequency per-atom data is available at

those you be set to: 1, 3 and 1, respectively.

you also need to define a callback and adjust the per MPI rank array sizes to accommodate that atoms migrate between MPI ranks.

you can look at, for example, fix ava/atom as an example for a fix with per-atom properties. the data is then stored in the class member array.

please also keep in mind that adding per-atom data causes an access ambiguity between the global vector and the per atom data. that should not be a problem if you only access the data from a dump, but can be a problem with other accesses and properties that are already provided.

  1. A second option is to recompute the drag force in fix langevin as a per atom vector variable in my LAMMPS script and track the individual, per-atom force contributions of every defined bond, angle, dihedral, and pair potential I have in my simulations, but somehow I think this would be more complicated and require more of an interface with the underlying C++ code.

there is a mechanism for that kind of data collection only for pair-wise additive pair styles, not for bonded interactions.

  1. Is there an equivalent of the strain energy as a thermodynamic variable in LAMMPS that I’m missing? Something that captures the work done by bond forces, angle forces, dihedral forces, pair potential forces, as well as the viscous drag force from a Brownian dynamics model that I’m overlooking such as ‘pe’ (though I doubt pe stores any information about the viscous drag)?

no.