In MD: Non-zero average forces acting on atoms of system at higher temperature

[alok_shukla1]

Hello,

I am trying to plot the forces acting on a my system. My expectation are zero average force on atoms at higher temperature as well as absolute zero temperature. I calculated avg force by using the following:

Avg force of a iterations= sum(modulus of force on each atom of that iterations)/(no of atoms)

Then plotted the forces with respect to iterations or time [ time in pico seconds are obtained by diving iteration number by 1000 cause one iteration is of one femto second].

From the plots of Si, one can see non-zero average force. Could any one explain why I am getting non-zero average force at higher temperature?

For example, for Silicon, I have performed molecular dynamics (MD) calculation for ~20 ps. And got forces ~1.15 eV/Angstrom at 500K.

Your help will be highly appreciable. Please do provide some references.

Please find the plots in the attachment:

[ferenc_karsai]

Please post your input and output files according to the forum guidelines (INCAR, KPOINTS, POSCAR, POTCAR, OUTCAR).

[alok_shukla1]

Dear ferenc_karsai,

Please find the INCAR, KPOINTS, POSCAR, POTCAR and OUTCAR files.

Since my OUTCAR file size is huge, I have deleted the bottom part of it. I have also attached a script file to extract forces from OUTCAR file.


Looking forward to your suggestions.


Thanks

[toms_bucko1]

In fact, it's the force components rather than their absolute values or the norm of the force vector that should average to zero. One can easily prove this analytically for a 1D harmonic system with potential energy V=C/2(q-q0)^2:

- force: F=-C(q-q_0)

-partition function: Q=\int_{-\infty}^{\infty} dq\,e^{-C(q-q_0)^2/2k_B T} = \sqrt{2\pi k_B T/C}

- mean force: <F> = -\frac{C}{Q} \int_{-\infty}^{\infty} dq\, (q-q_0) \,e^{-C(q-q_0)^2/2k_B T} = \sqrt{\frac{C k_B T}{2\pi}}[e^{-C(q-q_0)^2/2k_B T]_{-\infty}^{\infty} = 0

- mean abs. value of force: <|F|> = \frac{C}{Q} \int_{q_0}^{\infty} dq\, (q-q_0) \,e^{-C(q-q0)^2/2k_B T} + \frac{C}{Q} \int_{-\inft}^{q_0} dq\, (q_0-q) \,e^{-C(q-q0)^2/2k_B T} = \sqrt{2\pi k_B T/C} [-k_B T \,e^{-C(q-q_0)^2/2k_B T]_{q_0}^{\infty} - \sqrt{2\pi k_B T/C} [-k_B T \,e^{-C(q-q_0)^2/2k_B T]_{-\infty}^{q_0} = \sqrt{\frac{2C k_B T}{\pi}} > 0


As a quick numerical test, let's use the OUTCAR file provided by alok_shukla1 in his previous contribution and inspect (say) the x-component of the force on the first atom (feel free to look at other component of other atom - I guess it's pretty obvious how to modify the commands to achieve that):


grep -A 2 TOTAL-FORCE OUTCAR |awk 'BEGIN {a=0.;c=1} {if (NR%4==3){c+=1;a+=$4}} END {print a/c}'
0.0328574

as one can see, this value is already small although the MD was very short (<7000 steps) and one can therefore not expect to have a good convergence of the ensemble average. Visualizing (e.g., using gnuplot, xmgrace, origin,...) the data stored in xxx created as follows:

grep -A 2 TOTAL-FORCE OUTCAR |awk '{if (NR%4==3) {print $4}}' > xxx

reveals that the force component approximately oscillates symmetrically around zero, as it should. From this behavior it is clear that when one takes an absolute value of the same force component, one must necessarily obtain a positive value:

grep -A 2 TOTAL-FORCE OUTCAR |awk 'BEGIN {a=0.;c=1} {if (NR%4==3) {c+=1;a+=($4*$4)**0.5}} END {print a/c}'
0.638877

The behavior reported by alok_shukla1 is therefore correct.