Qbox List

Dear Francois,

We recently applied the SCAN functional to simulate water, and found that the sum of all atomic forces is not so close to zero.
Please see the attachment.
Using PBE, the sum of forces is [0.305, -0.509, -0.304] meV/angstrom, while it is [11.678, -23.4177,5.2982] meV/angstrom with SCAN.
Did we do anything wrong? Thanks!

Best,
Ding

Dear Ding,

I have looked at the issue, and it appears that this is the kind of numerical noise that results from using a finite plane wave cutoff. The total force can be reduced somewhat by increasing the cutoff, but that makes the calculations very costly. Note that the value of the total force you observed is still small if we consider the error in forces that affects individual atoms on average. The numerical error is larger for SCAN than for PBE likely because the SCAN xc energy depends more strongly on gradients of the density and wave functions (the kinetic energy density) and this introduces high Fourier components that converge slowly with large ecut.

I assume that you are concerned about the residual total force causing a drift during a long MD simulation (the "flying ice cube" problem). For this reason, I have added a new variable called lock_cm to Qbox. It can be used to lock the center of mass of the system during a simulation. You can do that by using set lock_cm ON before the run command. The position and velocity of the center of mass can be monitored using the status command at regular intervals (it should be exactly fixed and the velocity should be zero if lock_cm is ON).

The lock_cm variable is available in release 1.74.1 which was just published.

Please let me know about your experience if you decide to use this feature, especially if you encounter problems. I did not have time to test it extensively.

Best wishes,
Francois

Dear Francois,

Thank you for your prompt reply!

Sorry that "the sum of forces" in my original post is in fact the sum of atomic forces divided by the number of atoms.

We also compared the atomic forces calculated by Qbox and QE.
We used the same plane wave cut off of 85 Ry and the same pps. Please see the attachments.
Using PBE, the mean absolute error (MAE) between Qbox and QE is as little as 0.38 meV/Angstrom, while the MAE is 13.46 meV/Angstrom with SCAN.
It seems the SCAN calculation with Qbox may give a small residual force.

The new variable lock_cm is very helpful. We will read the source code and test it and let you know.
We are currently using reset_vcm every certain period of time, which also seems working well.

Best,
Ding

Dear Ding,

Thanks for the update. I looked at the QE results, and at first sight it would seem that QE has a much smaller residual total force on atoms. However, this is simply due to the fact that QE arbitrarily subtracts the sum of forces in the force calculation. See the following lines in the QE forces.f90 file:
https://github.com/QEF/q-e/blob/28725df ... .f90#L298

Therefore the sum of forces computed from the QE output only reflects the randomness of rounding of the eighth printed digit, of the order of 1.e-8. It is not clear what the residual force is in QE before subtracting the sum of forces.

I checked that the Qbox SCAN mean error (sum / natoms) can be reduced by increasing ecut. At ecut=300 Ry, the mean error is smaller than 3.e-5 Ha/bohr (1.6 meV/Å). At ecut=500 Ry, the mean error is smaller than 1.3e-5 Ha/bohr (0.66 meV/Å). At ecut=800 Ry, it is 5.2e-6 Ha/bohr (0.26 meV/Å). Of course using such large values of ecut is not practical, but it demonstrates that the implementation is correct. In practice, using the lock_cm variable is probably the best approach to avoid drift.

It has been noted that SCAN is numerically difficult, and two "regularized" SCAN functionals have been proposed (rSCAN and r2SCAN) See https://pubs.acs.org/doi/10.1021/acs.jpclett.0c02405. I have not considered implementing them yet, and I don't know if they would improve the accuracy of the sum of forces for moderate ecut.

Best,
Francois

Thank you very much Francois.