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

Question on input files/tags, interpreting output, etc.

Please check whether the answer to your question is given in the VASP online manual or has been discussed in this forum previously!

Moderators: Global Moderator, Moderator

Post Reply
Message
Author
alok_shukla1
Newbie
Newbie
Posts: 8
Joined: Sun Nov 17, 2019 3:21 am

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

#1 Post by alok_shukla1 » Sun May 15, 2022 4:23 am

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:
You do not have the required permissions to view the files attached to this post.

ferenc_karsai
Global Moderator
Global Moderator
Posts: 217
Joined: Mon Nov 04, 2019 12:44 pm

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

#2 Post by ferenc_karsai » Mon May 16, 2022 5:59 am

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

alok_shukla1
Newbie
Newbie
Posts: 8
Joined: Sun Nov 17, 2019 3:21 am

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

#3 Post by alok_shukla1 » Mon May 16, 2022 12:03 pm

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
You do not have the required permissions to view the files attached to this post.

toms_bucko1
Newbie
Newbie
Posts: 1
Joined: Wed Feb 19, 2020 5:51 pm

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

#4 Post by toms_bucko1 » Thu Jun 23, 2022 1:39 pm

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.

Post Reply