Dear all:
I’m trying to simulate the shear loading of a screw dislocation system with the flexible boundary conditions [PHYSICAL REVIEW MATERIALS 6, 013608 (2022)], which requires users to control the center of mass (COM) of boundary atoms. During the simulation, the boundary atoms are first given an initial velocity, and then the velocity of the COM is kept by cancelling out the average force acting on it.
In the following script, we mimic such an operation in a pristine Fe bulk system. The initial velocity is given a large value (i.e., variable v equal 1), and only the x component of the velocity of the COM of the upper layer is controlled. However, it is found that although the fx acting on the COM is almost zero (in the order of 10E-17), the x coordinate of the COM does not change linearly with respect of time. This means that the velocity of the COM of the upper layer is not controlled.
variable T equal 300
variable v equal 1 # The initial velocity of the upper layer
units metal
boundary p p s
atom_style atomic
timestep 0.001
read_data init.lmp
pair_style eam/alloy
pair_coeff * * /home/wsyang/W_screw/CRSS/test_Fe_screw/FeC.eam Fe
region upper block INF INF INF INF 131 INF units box # Define the geometry of upper layer
group upper region upper
# 1. Thermal equilibrium at ${T}
velocity all create ${T} 123456
fix 1 all nvt temp ${T} ${T} $(100*dt)
thermo_style custom step temp press pxz time lz
thermo 100
run 2000
unfix 1
# 2. The following code snippet intends to mimic the shear loading of a dislocation system
reset_timestep 0
compute uppercom upper com # The COM of upper layer
variable uppercom_x equal c_uppercom[1] # x coordinate of the COM
variable uppercom_y equal c_uppercom[2] # y coordinate of the COM
variable uppercom_z equal c_uppercom[3] # z coordinate of the COM
compute force_upper upper property/atom fx fy fz # The interatomic forces of the atoms in group “upper”
compute fx_ave_upper upper reduce ave c_force_upper[1] # Averaged fx
compute fy_ave_upper upper reduce ave c_force_upper[2] # Averaged fy
compute fz_ave_upper upper reduce ave c_force_upper[3] # Averaged fz
# Take the negative value of the averaged forces
variable com_fx_upper equal -1*c_fx_ave_upper
variable com_fy_upper equal -1*c_fy_ave_upper
variable com_fz_upper equal -1*c_fz_ave_upper
velocity upper set $v NULL NULL units box # Give the atoms in the upper layer an initial velocity
fix 2 upper addforce v_com_fx_upper 0.0 0.0 # Add the negative of averaged forces, in order to make sure that the x component of the force acting on COM is cancelled out.
fix 3 all nve
fix 4 all recenter INIT INIT INIT
# Export the positions of COM and the forces acting on COM
dump 2 upper custom 1 upper_layer.traj id x y z vx vy vz
fix save_com_upper all print 1 “${uppercom_x} ${uppercom_y} ${uppercom_z} ${com_fx_upper} ${com_fy_upper} ${com_fz_upper}” file upper_com.txt screen no
run 500
This is the excerpt of upper_com.txt, which contains the position (uppercom_x) and the fx of the COM (com_fx_upper). The velocity of COM is not kept at initial velocity.
# Fix print output for fix save_com_upper
41.5225401516348 41.5268708876333 136.349736033633 -9.86864910777917e-19 0.000736533811731725 0.00193306135662619
41.5234804061713 41.5269397440762 136.350020042822 -7.23700934570472e-18 0.000868785888786798 0.00175421337732491
41.5244206399892 41.5270084504159 136.35030374893 4.60536958363028e-18 0.000995305422799789 0.00157587520005436
41.5253608430891 41.5270769847932 136.350587182769 -4.60536958363028e-18 0.0011135987438461 0.00139840497668016
41.5263010054822 41.5271453267704 136.350870375 -9.86864910777917e-18 0.00122436134735561 0.00122330821090209
41.5272411172262 41.5272134572105 136.351153355876 1.71056584534839e-17 0.00132996234097077 0.00104999413924816
41.5281811684377 41.5272813578682 136.351436155341 2.63163976207445e-18 0.00143207032705165 0.000876834371772904
41.5291211493286 41.5273490111023 136.351718803312 1.39805862360205e-17 0.00153045938285516 0.000702937428621962
41.5300610502483 41.5274163999135 136.352001329834 9.86864910777917e-18 0.00162557960045334 0.000528788875884722
41.5310008617159 41.5274835078677 136.352283764996 1.38983474934557e-17 0.00171680720931052 0.000354438667139515
# Omitted
