[lammps-users] Fix addforce

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
variable f_z1 atom -0.6667

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!!


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

That looks like the correct way to apply a positition-dependent force.
Why don't you print out the forces in dump file for a run 0
command and convince yourself that the resulting
forces are different than if you had applied a constant force.


2011/1/18 Juan J. Meléndez <[email protected]>: