Third_order from PHONON pkg does not respect symmetry constraints

I’ve been trying to verify the output from the third_order command and ran across a weird issue: the force constants do not seem to follow the correct constraints. The most basic constraint is that force constants are invariant to the order of the index (tensor should be symmetric):
image

Below I pasted the 6 lines dumped third_order whose first terms should be equal.
The format of the output is i / alpha / j / beta / k gamma=1 gamma=2 gamma=3. With this I am just looking at gamma=1 to keep things simple. I expected all 6 of these to have the first gamma1 term but they do not. I am not sure if I am just interpreting the format wrong or if the output is wrong. Any help would be appreciated!

1 1 2 1 3       20.81668171        0.86736174       -5.20417043
1 1 3 1 2       10.40834086        0.00000000      -14.7451495
2 1 3 1 1       27.75557562      -13.87778781       -3.46944695
2 1 1 1 3       20.81668171        0.86736174       -5.20417043
3 1 1 1 2       10.40834086        0.00000000      -14.74514955
3 1 2 1 1       27.75557562      -13.87778781       -3.46944695

It is very difficult to provide a meaningful response without knowing more details.
For example, it would matter what kind of system you are studying and what kind of potentials you are using and whether you have properly minimized the system before running the third_order command. Since LAMMPS uses floating point math and the third_order command is based on finite differences, it is quite easily possible to have deviations from the ideal case, specifically of not everything else is “perfect” and you are using a potential that uses spline interpolation somewhere.

Since this is a very special command that not many of the LAMMPS developers (and users) have experience with, your best bet for assistance would be to contact the author that contributed the code. Contact info should be in the source code.

1 Like

Yeah I’ll try e-mailing.

The system is FCC at its equilibrium position with 4 unit cells and to be honest it should not matter the symmetry constraints are true regardless. It is effectively a statement of Clairaut’s theorem that derivatives should be interchangeable so the structure should not matter. It also does not matter if it is at equilibrium. And this isn’t really a small deviation it is completely different, hence I think my interpretation is wrong (or the output formatting is wrong).

Hello Ethan,

Axel is quite right about deviations arising due to finite differences.
Sometimes the differences arise from delta-x and/or potential cutoffs.

Please post your input file.

Best,
Charlie

Here’s the input file. I’ll try playing around with the inputs more to see if anything changes.

in.lammps_dynmat_test (1.3 KB)

Thanks!

I’ll see if I get some time to check it out this week, but first thing I would start off with is to change the delta x to 1e-6. Usually the sweet spot for delta x is somewhere around there. Also, you don’t need a skin of 1. Your skin just needs to be bigger than your delta x. Finally, one last thing is do not run the dynamical matrix command before hand, I need to fix a bug where the last movement gets undone.

1 Like

I tried a couple things but could not resolve the issues I was having. Changing the value of dx does have rather drastic effects on the results, but of the 3 cases (1e-6,1e-8, 1e-9) I tried none solved the issue. I’m also a little surprised that the results are so sensitive to the choice of dx. I would have thought if you have a stable configuration that overly small dx would not hurt your results.

Tried another code that uses finite differences & LAMMPS as a force calculator see here. There are no issues with symmetry here and it is really not that sensitive to the choice of dx. I have no clue what is going on (sorry not that helpful I know) but I’m pretty confident something is not right. Whether that is in the code or my input script I don’t know.

Something that’s always been a little suspicious to me is how fast the LAMMPS code runs. I’m not an MPI expert but the implementation in ModeCode and LAMMPS should not differ than much but the runtime is about 100x different.