Dear All,
I have a question on the accuracy of calculating force per atom (Class2 force field) in Lammps.
I am doing my own calculation for force per atom, and I keep getting the difference between my result and Lammps force per atom around 10^(-2) force unit.
I have checked carefully the Lammps source code on bond, angle, dihedral and improper force calculations. The equations are the same as mine which are derived from the potential energy.
I can’t find the source of such difference, can anyone give me some comments?
Is there any assumption used in Lammps to calculate the force per atom? How much accuracy we can expect?
Thank you for any comments,
Lili Zhang
Lammps uses double precision floating point math, which has about 15 digits relative accuracy. When summing multiple force terms, there can be truncation errors, and those depend on the order of summation and differences in magnitude.
You have to provide more details about what you compare, before one can speculate what is causing the difference.
Axel
Dear Dr. Kohlmeyer,
Thanks for your comment. Please see the attached file for the error analysis.
My procedure for calculating the force per atom is very straightforward:
(1) I used Lammps output of the atoms positions
(2) I derived the force equation from the Class2 potential functions: LJ, bond, angle, dihedral, improper and coulombic energies.
(3) sum all the forces on the atom I am interested.
My result is not the same as Lammps force per atom output.
So I am wondering if Lammps considers something more than I included.
Thank you very much,
Lili
Lili,
You have been at this for a long time already wondering why your code and lammps don’t yield the same results. I see two ways to put an end to your agony. One, you find a simple system with analytical solution to compare the numbers to. Or, you find another code that has the same potential implemented and use it as a third referee to rule one way or the other.
Carlos
My procedure for calculating the force per atom is very straightforward:
(1) I used Lammps output of the atoms positions
(2) I derived the force equation from the Class2 potential functions: LJ, bond, angle, dihedral, improper and coulombic energies.
(3) sum all the forces on the atom I am interested.
My result is not the same as Lammps force per atom output.
So I am wondering if Lammps considers something more than I included.
Ghost atoms? Neighbor list length? Cutoff function? These are also something to check.
Ray
Dear Ray,
For the purpose of checking our code, I didn’t use period boundary condition, so there shouldn’t be “ghost atoms,” am I right?
We used neighbor modify to make sure the building neighbor list in Lammps and ours is the same.
And we use the same cutoff distance too.
The disturbing thing is the magnitude of the difference with Lammps for the force per atom is not big, but there is a difference.
The reason we calculate the force by our own is we are developing continuum level stress calculation, which is not the Virial stress Lammps calculates.
Thanks,
Lili
Lili,
Since you seem to be familiar enough with the LAMMPS implementation of the class2 potential, why not turn everything off and then turn bond, angle, dihedral, and etc back on one by one while doing the same with your in-house code? This way you will at least know what part is causing you the trouble.
Ray