I try to compress a bonded network structure in one direction with two plates. For this I use two plane regions with a repulsive harmonic potential. This planes move with constant velocity and uniaxial compress the network structure. To obtain Young’s modulus of the network I have to calculate the change of the load on the plates per displacement.
Here is the latest version of my input script:
units lj
atom_style angle
boundary p p p
angle_style ...
bond_style ...
read_data network.dat
group net type 1
variable dt equal 0.0001
variable t equal step*${dt}
variable lux equal 169-v_t
variable uux equal 170-v_t
variable dx equal ramp(0,-84)
variable mdx equal ramp(0,84)
region unet block ${lux} ${uux} 0 170 0 170
group unet dynamic net region unet every 10
dump dum5 unet custom 10 ufx.dump fx
region uplane plane 170.0 5.0 5.0 1.0 0.0 0.0 side out move v_dx NULL NULL
region lplane plane 0.0 5.0 5.0 1.0 0.0 0.0 side in move v_mdx NULL NULL
pair_style ....
pair_coeff ...
fix 1 all nve
fix upot net wall/region uplane harmonic 250.0 0.0 1.0
fix lpot net wall/region lplane harmonic 250.0 0.0 1.0
compute ffx unet reduce sum fx
variable ff equal c_ffx
#fix uload unet print 10 "${mdx} ${ff} " file uload.dat screen no
timestep 0.01
run 840000
Using fix uload results in:
ERROR: Fix uload does not allow use with a dynamic group (…/modify.cpp:251)
In this input script I tried to define a dynamic (block) region (only for the upper plane), that is moving with each plate and reaches from the plate to the cutoff of the harmonic potential. With this region I wanted to define a group and calculate the sum of forces (or change of forces) of network particles inside this group.
But by looking at the ufx.dump it seems the group unet is always empty, for each time step the output is:
ITEM: TIMESTEP
…
ITEM: NUMBER OF ATOMS
0
ITEM: BOX BOUNDS pp pp pp
0.0000000000000000e+00 1.7000000000000000e+02
0.0000000000000000e+00 1.7000000000000000e+02
0.0000000000000000e+00 1.7000000000000000e+02
ITEM: ATOMS fx
I’m pretty sure the movement of the plates works as intended. Can someone please help with the calculation of the force/ load on the plates?
Thank you.
The quoted input is incomplete and missing the data file, which takes away many possible options to figure out what is happening by recreating your calculation and making modifications and observations.
That said, it is a needlessly convoluted and complex input. Why not use fix indent or fix wall/harmonic? Both provide direct access to the energy and force due to the interaction of the wall/plate with atoms.
Or even simpler, why not just set up a simple periodic bulk system and use fix deform? That is a much more conventional and simpler (and probably also more accurate) approach that does not have to deal with any interfaces.
Angle, bond and pair style differ from simulation to simulation. I use systems with DPD/SRP potentials, mixed with FENE bons/ harmonic bonds, WCA potential and cosine angle potential. The systems themselfs are always bead- spring models of polymernetworks in explicit or implicit solvent. The explicit solvent model uses beads with DPD potential and the implicit solvent model is a modified WCA potential (- without explicit solvent particles). With the explicit solvent model the “fix deform” approach did not yield the right behavior of the system during simulation. The explicit solvent has to be able to evade the moving walls.
I’m sure it is a needlessly convoluted and complex input
Without the need to compute Young’s modulus the needlessly convoluted and complex input got the job done. How would the input look with fix wall/harmonic? There would not be any need for all the additional regions and groups and the fix wall region would be replaced by fix wall harmonic:
fix uwall net wall/harmonic xhi v_uux 250.0 0.0 1.0
fix lwall net wall/harmonic xlo v_t 250.0 0.0 1.0
the doc page says " A lo face interacts with particles near the lower side of the simulation box in that dimension. A hi face interacts with particles near the upper side of the simulation box in that dimension." So if I use “xlo” with a wall at the lower side of the simulation box does the wall interact with the particles between the wall in direction of the center of the simulationbox up to the cutoff distance to the wall? How can the force due to the interaction of the wall with atoms be accessed directly? Do I use the fix ID in conjunction with a fix print command?:
fix uload net print 10 "${t} f_uwall[2] " file uload.dat screen no
The doc page of “fix_wall” states : " This fix computes a global scalar energy and a global vector of forces, which can be accessed by various output commands. Note that the scalar energy is the sum of interactions with all defined walls. If you want the energy on a per-wall basis, you need to use multiple fix wall commands. The length of the vector is equal to the number of walls defined by the fix. Each vector value is the normal force on a specific wall. Note that an outward force on a wall will be a negative value for lo walls and a positive value for hi walls. The scalar and vector values calculated by this fix are “extensive”."
results in “ERROR: Thermo fix vector is accessed out-of-range (…/thermo.cpp:949)”
It seems only f_XYZ[1] is in range.
Like the doc pages states: “The length of the vector is equal to the number of walls defined by the fix. Each vector value is the normal force on a specific wall”
Another problem is that fix wall/harmonic doesn’t allow a periodic boundry in normal direction to the wall. I fixed this problem with two additional walls at the borders of the simulation box in order to keep all solvent particles inside the simulation box and I hope this wont change the mean density of the system too much.
Since the the wall’s normal is parallel to the x-, y-, or z-axis, there is no point in storing the other force components. They would be zero. The wall only acts in the corresponding x-, y-, or z-direction.
By definition, if you replace a periodic direction with two walls (that keep the cell dimensions the same), the mean density of the system remains constant.
If you are unlucky, the mean density of the system will be the only thing that remains constant.
Near a planar wall, the local environment of a particle is fundamentally anisotropic, whereas that same particle’s local environment is isotropic (at a large enough length scale) in a bulk simulation under periodic boundaries. I would expect pronounced transverse layering and altered parallel correlations near such a wall, and I would always assume that molecular behaviour will be different from in the bulk, unless proven otherwise. In your particular situation, I would not be confident that there is an elastic modulus for your system, rather than a spatially varying elastic modulus and hence inhomogeneous (and even possibly nonlinear!) strain responses to the stress throughout the material.
This is just my opinion, of course, but I hope you will take the time to verify that your simulation is physically relevant, not just that it runs without crashing.
To answer the question of your error, using anything except the “all” group with fix print doesn’t make sense. Just change it to fix uload all print 10 ... and your system should run, but you should take everyone else’s advice into account concerning the practicality of walls.
In terms of the region processing, it might be more easily done in post by dumping the forces of all atoms as well as the positions, then write a script to select the atoms within that simple cutoff of the wall location, which you’re outputting with the fix print command. Then you won’t have to deal with dynamic groups or regions at all.
My only other concern is that you already have a free surface at the location of the wall, and you’re not trying to separate bonded atoms? The best way to get us to help you is to provide a simple test input with the issues you’re facing, the necessary data files to run that input, and (if relevant) the command you executed the simulation with.
