Dear developers,
In a Metallic model, the applied force on a rigid group is read from “3 global vector of forces” computed by command FIX SETFORCE 0 0 0 .
The problem is the fluctuation of force values are more than ten time of real values. It other words the force plot is too rough (attached image).
What is wrong here? How can I acquire more smooth output?
To make the script easier to look, I’ve simplified my model. It contains two blocks of Si each includes Rigid/Thermostat/Free layers. The force applied on the top and the resultant force is read from below.
Regards.
Rouhollah,
########## Initialization ########
units metal
dimension 3
boundary p p s
atom_style atomic
variable Gpa equal 0.1 #external pressure in Giga pascal
######### create geometry ###########
lattice diamond 5.431 #Si
region s_box block 0.0 30.0 0.0 10.0 0.0 20.0
create_box 3 s_box
############### Lower Block ##############
region rigid_down block 0.0 30.0 0.0 10.0 0.0 2.0
region thermos_down block 0.0 30.0 0.0 10.0 2.0 4.0
region newt_down block 0.0 30.0 0.0 10.0 4.0 8.0
region block_down union 3 newt_down thermos_down rigid_down # Total block_down volume
create_atoms 1 region block_down
############## Upper Block ##############
region newt_up block 0.0 30.0 0.0 10.0 8.5 12.0
region thermos_up block 0.0 30.0 0.0 10.0 12.0 14.0
region rigid_up block 0.0 30.0 0.0 10.0 14.0 16.0
region block_up union 3 newt_up rigid_up thermos_up
region mobile_up union 2 newt_up thermos_up
create_atoms 2 region block_up
################### Assign Groups to various regions ############
group RIGID_UP region rigid_up
group t_up region thermos_up
group NEWT_UP region newt_up
group mobile_up region mobile_up
group block_up region block_up
group THERMO_UP subtract t_up RIGID_UP NEWT_UP
group RIGID_DOWN region rigid_down
group t_down region thermos_down
group NEWT_DOWN region newt_down
group THERMOS_DOWN subtract t_down RIGID_DOWN NEWT_DOWN
group rigid_up_down union RIGID_UP RIGID_DOWN
group THERMO union THERMO_UP THERMOS_DOWN
group NEWTON union NEWT_UP NEWT_DOWN # Newtonian layer for block_up and block_down
group FREE subtract all rigid_up_down # (= thermostat+newtonian)
set group as particular type
set group RIGID_DOWN type 1
set group THERMOS_DOWN type 2
set group NEWT_DOWN type 3
set group RIGID_UP type 1
set group THERMO_UP type 2
set group NEWT_UP type 3
############### Interatomic potential ##################
pair_style tersoff
pair_coeff * * SiC.tersoff Si Si Si
mass 1 28.0855 ## Si
mass 2 28.0855 ## Si
mass 3 28.0855 ## Si
############# initial velocities relaxation #############
velocity all create 300 1235 mom yes rot yes dist gaussian
minimize 1.0e-25 1.0e-25 5000 10000
RELAXATION
fix 1 FREE nvt temp 300.0 300.0 0.2
fix 2 RIGID_DOWN setforce 0.0 0.0 0.0
fix 3 RIGID_UP rigid/nvt single temp 300.0 300.0 0.2
thermo 500
thermo_style custom step time cpu cpuremain timeremain temp epair emol etotal press
############ read forces from fix set 0 0 0 ##########
variable tt equal time
variable FX equal f_2[1]
variable FY equal f_2[2]
variable FZ equal f_2[3]
fix myprint all print 100 “{tt} {FX} {FY} {FZ}” &
file OUTPUT.txt screen no title “TIME Fx Fy Fz”