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:
read_restart restart.sigma5.dwivedi
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?
Thanks in advance!!
Juanjo