Does anyone have any relevant examples (lammps input files, analysis code, etc) for calculating elastic constants at finite T from a constant stress (variable cell) simulation?
I’m trying to test a variable cell Langevin implementation (for ASE), but I’m having a hard time getting consistency between my T=0 finite difference (matscipy.elasticity) and finite T strain fluctuation results. I was hoping to use LAMMPS as a “known good” variable cell MD implementation, but my analysis isn’t working for those results either.
I can think of 3 plausible reasons:
- LAMMPS variable cell algorithms don’t sample the correct ensemble (seems very unlikely)
- algorithms are correct but I’m not getting to equilibrium properly because of poorly chosen parameters (pretty plausible)
- my analysis is incorrect (seems most likely, but that’s why I’d like known good data to apply it to)
I found some very old threads here, but they link to other things on the old sandia-hosted mailing list archives, and nothing concrete enough that I can use to analyze my own trajectories.
I’m getting better agreement now, using fix npt and a sufficiently low temperature, although I still haven’t gotten fix press/langevin to work. Less urgent, but if anyone has tried this, I’d be interested in comparing notes.
Hi @noam.bernstein,
Please have a look at other threads on the topic, for example this one. The paper I link is previous work where we convergence issues in the NPT ensemble to compute the full elastic constant tensor.
I believe there were some corrections in the algorithm since (my PhD advisor told me so) but I never tested the new version and I still think some of our results hold: i.e. there is some impact from the coupling constant of the NH integrator on the sampling. This is something that also stemmed out from discussion with other people in the field. So I don’t think your analyses is incorrect but this requires careful processing to be sure about convergence anyway.
If you have issues with press/langevin please tell or open a PR on Github, I did the implementation and would be eager to correct any problem encountered.
Thanks @Germain . I’ll take a look at the other thread and the linked paper.
Meanwhile, I’ll post my lammps input with the npt and press/langevin setups, and see if you think it’s reasonable, before I open a PR.
Here’s a summary. The fix lines I tried were
variable tdamp equal 0.2
variable pdamp equal 10.0
fix NPT all npt temp ${T} ${T} ${tdamp} tri 0 0 ${pdamp}
and
fix NVT all nvt temp ${T} ${T} ${tdamp}
fix PRESS_LANG all press/langevin tri 0 0 ${pdamp} temp ${T} ${T} 200
With FD T=0 I get
Cij FD [[ 2.533 1.28 1.28 -0. -0. -0. ]
Cij FD [ 1.28 2.533 1.28 -0. 0. 0. ]
Cij FD [ 1.28 1.28 2.533 -0. -0. 0. ]
Cij FD [-0. -0. -0. 1.272 0. 0. ]
Cij FD [-0. 0. -0. 0. 1.272 0. ]
Cij FD [-0. 0. 0. 0. 0. 1.272]]
With with the strain fluctation analysis I get
Cij LAMMPS NPT fluct [[ 2.469 1.444 1.238 -0.071 -0.118 0.054]
Cij LAMMPS NPT fluct [ 1.444 2.427 1.271 0.001 -0.036 -0.034]
Cij LAMMPS NPT fluct [ 1.238 1.271 2.639 0.073 0.159 -0.023]
Cij LAMMPS NPT fluct [-0.071 0.001 0.073 1.231 -0.145 0.048]
Cij LAMMPS NPT fluct [-0.118 -0.036 0.159 -0.145 1.244 -0.184]
Cij LAMMPS NPT fluct [ 0.054 -0.034 -0.023 0.048 -0.184 1.346]]
Cij LAMMPS NVT_PRESS_LANG fluct [[89.123 24.077 23.472 1.718 10.133 -0.683]
Cij LAMMPS NVT_PRESS_LANG fluct [24.077 73.889 27.964 -0.548 -0.068 -0.765]
Cij LAMMPS NVT_PRESS_LANG fluct [23.472 27.964 77.726 -0.082 0.844 0.84 ]
Cij LAMMPS NVT_PRESS_LANG fluct [ 1.718 -0.548 -0.082 56.279 -0.193 1.133]
Cij LAMMPS NVT_PRESS_LANG fluct [10.133 -0.068 0.844 -0.193 16.636 -0.633]
Cij LAMMPS NVT_PRESS_LANG fluct [-0.683 -0.765 0.84 1.133 -0.633 56.817]]
Shall I open a PR? Do you want to see more details?
Is you system an Argon LJ crystal? This is to be expected from finite differences due to symmetries. Note that you have to compare finite T with finite T to get an idea of how the Cij evolve with temperature. There is a contribution from stress fluctuations due to temperature.
Coherent and relatively converged but what is your timestep? Pdamp seems low for either metal or real units. It can affects the systems fluctuations. Did you check that you get similar fluctuations in both cases and that your average pressure is correct?
Is you system an Argon LJ crystal?
No, it’s a Morse potential with random parameters that were originally meant to roughly reproduce Al, but something got screwed up in the unit conversion.
I’m basically satisfied with the agreement of LAMMPS NPT. It’s a 2.5 fs timestep, so 0.0025 in metal units, and tdamp is therefore 80 dt, and pdamp is 4000 dt = 50 tdamp.
It’s just the LAMMPS press/langevin that seems very wrong.
The equations are different from NH barostat in the NPT integrator. Therefore there is no direct comparison possible simply with the Pdamp parameter. This is why I asked about your data and the pressure of your simulations.
The goal of barostat is to keep the average pressure correct. This is what they are made for. The fluctuations depend on the underlying equation and algorithm implementation. The langevin simulation is maybe not converged and the fluctuations are maybe very large because the Pdamp parameter is too high or low to reproduce acceptable dynamics. A way to tell if you’re far would be to compare with the pressure fluctuation in NVT for a large system at equilibrium density. At the thermodynamic limits, they should be close in both ensembles at equilibrium. If there is no Pdamp value where this happens, and we do not see convergence of the Cij, then we might have a problem to investigate.
I’ll try press/langevin with other pdamp parameters and see if that fixes anything. The time scales for the pressure fluctuations seem roughly similar, but I haven’t looked any deeper yet.
Indeed, the cell fluctuations are much much faster than pdamp, which is obviously bad. I’ll increase it a lot, and see if the results look better. The fact that the heuristic that sets the barostat time scale is so far off makes this fix harder to use than I was hoping.
This is what it looks with pdamp of 20, 50 and 200, and NVT (with pdamp 10) for comparison. The Cij values from press/langevin are all completely wrong. The fluctuations are all much too small, and the resulting elastic constants are much too large (values by finite difference and fix npt are of order 1)
Cij LAMMPS NVT_PRESS_LANG 20.0 fluct [[210.115 -20.778 -15.698 2.941 21.39 2.516]
Cij LAMMPS NVT_PRESS_LANG 20.0 fluct [-20.778 211.428 -13.959 10.724 44.381 -1.156]
Cij LAMMPS NVT_PRESS_LANG 20.0 fluct [-15.698 -13.959 186.307 -8.395 12.664 -8.434]
Cij LAMMPS NVT_PRESS_LANG 20.0 fluct [ 2.941 10.724 -8.395 211.31 -4.141 -5.855]
Cij LAMMPS NVT_PRESS_LANG 20.0 fluct [ 21.39 44.381 12.664 -4.141 206.15 7.65 ]
Cij LAMMPS NVT_PRESS_LANG 20.0 fluct [ 2.516 -1.156 -8.434 -5.855 7.65 218.452]]
5000it [00:03, 1328.39it/s]
Cij_fluct traj len 5000
T 9.898106206946093 V 1834.9973814828866
Cij LAMMPS NVT_PRESS_LANG 50.0 fluct [[650.78 90.612 9.541 39.709 -0.56 23.084]
Cij LAMMPS NVT_PRESS_LANG 50.0 fluct [ 90.612 556.1 152.875 -61.335 -41.687 -30.636]
Cij LAMMPS NVT_PRESS_LANG 50.0 fluct [ 9.541 152.875 639.926 39.695 119.552 8.243]
Cij LAMMPS NVT_PRESS_LANG 50.0 fluct [ 39.709 -61.335 39.695 724.939 95.032 -23.284]
Cij LAMMPS NVT_PRESS_LANG 50.0 fluct [ -0.56 -41.687 119.552 95.032 821.678 -26.907]
Cij LAMMPS NVT_PRESS_LANG 50.0 fluct [ 23.084 -30.636 8.243 -23.284 -26.907 790.07 ]]
5000it [00:03, 1334.73it/s]
Cij_fluct traj len 5000
T 9.903119026219077 V 1835.0035241656099
Cij LAMMPS NVT_PRESS_LANG 200.0 fluct [[2688.136 1454.571 1513.519 -182.512 104.012 142.538]
Cij LAMMPS NVT_PRESS_LANG 200.0 fluct [1454.571 2899.704 1563.795 506.997 146.528 -347.574]
Cij LAMMPS NVT_PRESS_LANG 200.0 fluct [1513.519 1563.795 2316.49 -325.281 248.97 -413.78 ]
Cij LAMMPS NVT_PRESS_LANG 200.0 fluct [-182.512 506.997 -325.281 1952.111 -281.512 61.23 ]
Cij LAMMPS NVT_PRESS_LANG 200.0 fluct [ 104.012 146.528 248.97 -281.512 1523.481 344.737]
Cij LAMMPS NVT_PRESS_LANG 200.0 fluct [ 142.538 -347.574 -413.78 61.23 344.737 2239.539]]
traj.10.NVT_PRESS_LANG_20.0.propag_equil.dump.x_cell.pdf (18.9 KB)
traj.10.NVT_PRESS_LANG_50.0.propag_equil.dump.x_cell.pdf (20.7 KB)
traj.10.NVT_PRESS_LANG_200.0.propag_equil.dump.x_cell.pdf (18.7 KB)
traj.10.NPT.propag_equil.dump.x_cell.pdf (16.5 KB)
I didn’t have time to have a look into this specific issue after finishing the implementation.
I’ll try to spend some time on it this week. Three questions:
- How did you compute the Cij? Only with the ASE matscipy module?
- What is you trajectory length? 5000 configuration? This is rather short by the way, but the discrepancy is still clear.
- What is the symmetry of your system? Simple cubic like FCC?
The Cij lines that say “FD” in one of the earlier messages are matscipy elasticity using LAMMPSlib with the same pair style as the ASE calculator. Everything else is fluctuations, with lammps (but I also ran my own variable cell Langevin ASE module that agrees with lammps npt).
The system is fcc Al, relaxed before starting, and I do 5e5 steps as equilibration and then sample for 5e5 steps, one config every 100 (so 5000 configs), for the fluctuation calculations.
1 Like