Hi all:

I’m trying to simulate impurity segregation to a grain boundary in ZrO2. To this end, an elastic force must be applied to each impurity atom, which is one of the driving forces for atomic movement.

In particular, I am applying some force to yttrium atoms by means of the “fix addforce” command. The applied force on each atom is not constant, but depends on the atom’s position. I have written it as follows:

kspace_style ewald 1.0e-4
velocity all create 1500.0 82986 dist gaussian mom no

# Calculating the force over each yttrium atom

variable f_x1 atom -0.6667*(x-10.26)/sqrt((x-10.26)(x-10.26)+(y-12.17)(y-12.17)+(z-16.2325)(z-16.2325))
variable f_y1 atom -0.6667
(y-12.17)/sqrt((x-10.26)(x-10.26)+(y-12.17)(y-12.17)+(z-16.2325)(z-16.2325))
variable f_z1 atom -0.6667
(z-16.2325)/sqrt((x-10.26)(x-10.26)+(y-12.17)(y-12.17)+(z-16.2325)*(z-16.2325))

fix fuerza_1 yttrium addforce v_f_x1 v_f_y1 v_f_z1

# Calculating msd within the grains
compute dcmoxi_bulk oxygen_bulk msd com yes
compute dcmzir_bulk zirconium_bulk msd com yes
compute dcmytr_bulk yttrium_bulk msd com yes

# Same in the GB
compute dcmoxi_gb oxygen_gb msd com yes
compute dcmzir_gb zirconium_gb msd com yes
compute dcmytr_gb yttrium_gb msd com yes

#Output settings
thermo_style custom step temp etotal press lx ly lz vol c_dcmoxi_bulk[1] c_dcmoxi_bulk[2] c_dcmoxi_bulk[3] c_dcmoxi_bulk[4] c_dcmzir_bulk[1] c_dcmzir_bulk[2] c_dcmzir_bulk[3] c_dcmzir_bulk[4] c_dcmytr_bulk[1] c_dcmytr_bulk[2] c_dcmytr_bulk[3] c_dcmytr_bulk[4] c_dcmoxi_gb[1] c_dcmoxi_gb[2] c_dcmoxi_gb[3] c_dcmoxi_gb[4] c_dcmzir_gb[1] c_dcmzir_gb[2] c_dcmzir_gb[3] c_dcmzir_gb[4] c_dcmytr_gb[1] c_dcmytr_gb[2] c_dcmytr_gb[3] c_dcmytr_gb[4]
thermo 100
thermo_modify lost warn

#Molecular Dynamics settings
fix 4 all temp/berendsen 1500.0 1500.0 0.5
fix 5 all nph iso 0 0 0.5 nreset 1

run 100000

The issue is that I get the same result when I use a constant force (i.e., when I set f_x1 atom 0.5 and so on). So is there anything wrong in my script? If so, how can I apply an atom-dependent force?