Merge pull request lammps#4405 from tylercollins5737/develop

vcm example script and log added

doc/src/Examples.rst

@@ -146,6 +146,8 @@ Lowercase directories
 +-------------+------------------------------------------------------------------+
 | streitz     | use of Streitz/Mintmire potential with charge equilibration      |
 +-------------+------------------------------------------------------------------+
+| stress_vcm  | removing binned rigid body motion from binned stress profile     |
++-------------+------------------------------------------------------------------+
 | tad         | temperature-accelerated dynamics of vacancy diffusion in bulk Si |
 +-------------+------------------------------------------------------------------+
 | threebody   | regression test input for a variety of manybody potentials       |

doc/src/compute_temp_chunk.rst

@@ -184,11 +184,24 @@ temp/chunk calculation to a file is to use the
 The keyword/value option pairs are used in the following ways.
 
 The *com* keyword can be used with a value of *yes* to subtract the
-velocity of the center-of-mass for each chunk from the velocity of the
-atoms in that chunk, before calculating either the global or per-chunk
-temperature.  This can be useful if the atoms are streaming or
+velocity of the center-of-mass (VCM) for each chunk from the velocity of
+the atoms in that chunk, before calculating either the global or per-chunk
+temperature. This can be useful if the atoms are streaming or
 otherwise moving collectively, and you wish to calculate only the
-thermal temperature.
+thermal temperature. This per-chunk VCM bias can be used in other fixes and
+computes that can incorporate a temperature bias. If this compute is used
+as a temperature bias in other commands then this bias is subtracted from
+each atom, the command runs with the remaining thermal velocities, and
+then the bias is added back in. This includes thermostatting
+fixes like :doc:`fix nvt <fix_nh>`,
+:doc:`fix temp/rescale <fix_temp_rescale>`,
+:doc:`fix temp/berendsen <fix_temp_berendsen>`, and
+:doc:`fix langevin <fix_langevin>`, and computes like
+:doc:`compute stress/atom <compute_stress_atom>` and
+:doc:`compute pressure <compute_pressure>`. See the input script in
+examples/stress_vcm for an example of how to use the *com* keyword in
+conjunction with compute stress/atom to create a stress profile of a rigid
+body while removing the overall motion of the rigid body.
 
 For the *bias* keyword, *bias-ID* refers to the ID of a temperature
 compute that removes a "bias" velocity from each atom.  This also

examples/README

@@ -112,6 +112,7 @@ snap:     examples for using several bundled SNAP potentials
 srd:      stochastic rotation dynamics (SRD) particles as solvent
 steinhardt: Steinhardt-Nelson Q_l and W_l parameters usng orientorder/atom
 streitz:  Streitz-Mintmire potential for Al2O3
+stress_vcm: removing binned rigid body motion from binned stress profile
 tad:      temperature-accelerated dynamics of vacancy diffusion in bulk Si
 template: examples for using atom_style template and comparing to atom style molecular
 tersoff:  regression test input for Tersoff variants

examples/stress_vcm/README

@@ -0,0 +1,32 @@
+README stress_vcm 
+=================
+
+Contents:
+
+-   in.stress_vcm: Example script showing how to remove binned 
+    velocities of center of mass (VCM) from stress calculations.
+-   stress_comparison.19Nov24.png: Plot shows the stress 
+    calculated in bars on the y axis for each positional bin on 
+    the x axis. Plotted are three different time steps from 
+    stress profiles with and without the VCM removed. Plot 
+    generated using Python.
+-   stress_xx.19Nov24.out: Output file generated by fix ave/time.
+-   log.19Nov24.stress_vcm.g++.1: LAMMPS log file with 1 proc.
+-   log.19Nov24.stress_vcm.g++.4: LAMMPS log file with 4 procs.
+
+Notes:
+
+-   Running this script as-is will generate two files. A log 
+    file with thermodynamic data and a stress_xx.out file 
+    containing the binned stress profile with the VCM removed.
+-   To generate the binned stress profile without removing the 
+    VCM then the compute stress/atom command at step three 
+    needs the last keyword "ch_temp_vcm" to be replaced with 
+    "NULL". 
+-   Uncommenting the line under "Atom dump" will generate an 
+    all atom dump file every 50 time steps containing atom ID,
+    type, and xyz coordinates.
+-   Uncommenting the lines under "Image dumps" will generate 
+    .jpg image files every 250 timesteps.
+-   Uncommenting lines under "Movie dump" will generate a .avi 
+    movie file showing timesteps every 125 timesteps. 

examples/stress_vcm/in.stress_vcm

@@ -0,0 +1,113 @@
+# Removing Binned Velocities of Center of Mass (VCM) from Stress 
+
+# This example shows how to remove rigid body motion from 
+# binned stress calculations. This uses a combination of commands
+# from compute chunk/atom, compute temp/chunk, compute 
+# stress/atom and fix ave/time. We'll show how these commands 
+# work in the context of a shockwave experiment on a cube of 
+# atoms. To shock the cube, a rectangular region of atoms is 
+# frozen, moved into the cube with a constant velocity along the
+# x direction, and then unfrozen. As the shockwave begins 
+# propagating, the body of the cube also moves along the x 
+# direction. To better understand the stress dynamics of the 
+# cube we remove the velocity component belonging to the overall 
+# motion of each bin.
+
+units		    metal
+boundary        p p p
+atom_style	    atomic
+lattice		    fcc 5.3589
+processors 		1 * *
+
+# Defining regions for box and atoms.
+# In this experiment an elongated simulation cell is 
+# defined in the x direction to allow for non-periodic 
+# motion of the atoms.
+
+region		    box1 block -3 24 0 12 0 12 units lattice
+region		    box2 block 0 12 0 12 0 12 units lattice
+
+# Creating box and atoms
+
+create_box	    1 box1
+create_atoms    1 region box2
+
+mass		    1 40.00
+
+# Adding energy to the system
+
+velocity        all create 600.0 9999
+
+pair_style	    lj/cut 10 
+pair_coeff	    1 1 0.04 3.405
+
+# Begin time integration
+
+timestep        2e-3
+
+fix	            fix_nve all nve
+
+thermo		    100
+
+run             500
+
+#--------------------------------------#
+# Chunk, Stress, and VCM removal steps #
+#--------------------------------------#
+
+# 1. Create 20 equispaced bins sliced along the x direction.
+# -"units reduced" normalizes the distance from 0.0 to 1.0
+variable        nbins index 20
+variable        fraction equal 1.0/v_nbins
+variable        volfrac equal 1/(vol*${fraction}) 
+compute         ch_id all chunk/atom bin/1d x lower ${fraction} units reduced 
+
+# 2. Calculate temperature bins with VCM aka COM velocities removed.
+compute         ch_temp_vcm all temp/chunk ch_id com yes
+
+# 3. Compute per atom stress with VCM removed via temp-ID.
+# -The velocities from specified temp-ID are used to compute stress.
+# -Stress/atom units are pressure*volume! Optionally handled next step.
+compute         atom_stress_vcm all stress/atom ch_temp_vcm
+
+# 4. Divide out bin volume from xx stress component.
+variable        stress atom -(c_atom_stress_vcm[1])/(vol*${fraction})
+
+# 5. Sum the per atom stresses in each bin.
+compute         ch_stress_vcm all reduce/chunk ch_id sum v_stress
+
+# 6. Average and output to file.
+# -The average output is every 100 steps with samples collected 20 times with 5 step intervals.
+fix             ave_stress_vcm all ave/time 5 20 100 c_ch_stress_vcm mode vector file stress_xx.out
+
+#--------------------------------------#
+
+# Piston compressing along x direction
+
+region 			piston block -1 1 INF INF INF INF units lattice
+group          	piston region piston
+fix 			fix_piston piston move linear 5 0 0 units box # strain rate ~ 8e10 1/s 
+
+thermo_style    custom step temp ke pe lx ly lz pxx pyy pzz econserve 
+
+# Atom dump
+
+# dump		    atom_dump all atom 50 dump.vcm
+
+# # Image dumps
+
+# dump		    2 all image 250 image.*.jpg type type &
+# 		        axes yes 0.8 0.02 view 60 -30
+# dump_modify	    2 pad 1
+
+# # Movie dump
+
+# dump		    3 all movie 125 movie.avi type type &
+# 		        axes yes 0.8 0.02 view 60 -30
+# dump_modify	    3 pad 1
+
+run	            500
+
+unfix           fix_piston
+
+run             1500