Elastic constants and box sizes

Dear LAMMPS users,
I am developing a MEAM potential but I have a problem obtaining the same/similar results among different box sizes, in theory using periodic boundaries should solve this situation but it doesn’t, I can understand if there is a small fluctuation between different box sizes but the results that I obtain gives totally different results, I am using esencially the same script as in the example, I have to mention that the position of atoms is defined by a bulk file.
I have searched for an answer in the forums including lowering the strain (“up” variable) until it converges but none seems to solve the problem.
Although making tests I have noticed that with some combinations of the strain (“up” variable), and the energy tolerance (“Etol”) and the force tolerance (“Ftol”) for the minimization gives, sometimes, constant results along different box sizes for some structures but not for all.

Any hint on how to solve this is much appreciated,

Mario Muralles

There is not enough info in your email to guess if the problem

is with your potential or your input script or if variation in the results

is normal.


One thing you could try is to switch from your meam potential to another
"stable" potential that have been used to model the same system while
leaving the rest of your input unchanged. If the problems remain at least
you know that something else is wrong beyond your new potential
implementation. Note as well that the elastic constant example script
provided with the lammps distro was updated very recently after a mistake
was identified in such a script. Note sure this could be triggering your
issues but make sure you start your work employing the corrected version of
the script.

Thank Steve and Carlos you for the answers,

Sorry to not specify about my inputs, I am developing a potential with the objective of reproduce a phase transition with the correct elastic constants in each phase, right now I am setting the potential to each individual phase in it low state energies aproaching to the right parameters by searching the screening parameters, this meaning that I am testing many potentials, I am working with a cubic and a monoclinic structure representing each phase, at the beginning I started to calculate the elastic constants on a unit cell and since the boundary is set to periodic I thought that there shouldn’t be any problem, then as a proof-check I tested bigger boxes having different results.
As Carlos mentioned about the issues with the old version of the script to get the elastic constants I decide on checking on the new version of the script and indeed it solved the problems relating to the cubic structure but I still observe the size problems with the monoclinic structure, I have to say that some of the potentials show good agreement about the elastic constants but others present an abismal differences.

Thank you for your help and comments,


Upgrading to the new script was just a step. Do as I suggested and run your tests using an already developed/tested potential. Debugging is easier with fewer variables. Theoretically one cannot expect to see big differences between Cijs by just changing the box size. You can always send your input files although there is no guarantee people will spend time looking at them. So far from your description the problem seems to be located on your side of the fence thus little motivation for the developers to chime in. Yet, users sometimes step up to certain tasks (me being one of them) when properly motivated.

"it solved the problems relating to the cubic structure." I doubt that
this is true, since the error in the old script only affected
non-orthorhombic cells.

Based on your description, you are clearly asking too much of the
elastic constant script. It will only give correct (size-independent)
results when all the relaxation steps are adequately converged. The
only way to tell if this is true is by carefully looking at the
output. If you are sampling a space of arbitrary potential
parameters, it is likely that in many cases the script is failing to
converge the initial cell relaxation, before it even gets to apply

Before you adopt this shotgun (seen Craig Ventner) approach, you first
need to develop a robust protocol for detecting unconverged


Thank you for your comments and help,
As Aidan mention the problem with the cubic structure was related with the relaxation, the old script presented some problems with lower force tolerances, giving me different results, now it shows me the same results converging after a certain force tolerance as expected.
Also as Carlos mentioned I have tried other potentials (made of the info given in the papers) and still the problem continued.

At the end I found that for the monoclinic structure something happens by the time of reading the data file, I used to generate the coordinates of the atoms in MATLAB defining the unit cell and making larger boxes in the same way as the command replicate does, but after trying everything else I decided to try to define them in the input script, I used an input found in the forum and it didnt work so well, then I decided to use my parameters and it totally worked, problem solved, still don’t understand why it didn’t work with the data file and specially with the other input.

Anyway thank you so much for your help,