Dear LAMMPS Users,

I am trying to replicate the system described in http://dx.doi.org/10.1103/PhysRevB.86.045213 however I am having trouble obtaining the correct elastic constants.

It is PbTe with NaCl-type structure.

Potential includes Coulomb part, Morse part between Pb-Te and Te-Te pairs, and cosine/squared interaction for Te-Pb-Te angle.

pair_style hybrid/overlay morse/smooth/linear 8 coul/cut 20.0

pair_coeff * * coul/cut 20.0

pair_coeff 1 2 morse/smooth/linear 0.465 0.863 3.68 7

pair_coeff 2 2 morse/smooth/linear 0.494 1.51 4.22 8

angle_style cosine/squared

angle_coeff 1 0.340 90.0

I am setting the system via read_restart command from custom made file (PbTe with NaCl-type structure and angles: twelve 90 degree angles Te-Pb-Te per every Pb atom).

When I try to calculate the elastic constants, the minimization (fix 3 all box/relax aniso 0.0) does not reach pressure close to 0 (i.e. it is usually lower than -100 000 bars and does not depend on initial lattice parameter), despite zeroing forces. Also, the minimization itself sometimes converges to unphysical energy values (e.g. 1e+200 eV).

I am using the ELASTIC example.

I have tried different cutoff distances, initial lattice parameters, types of minimization (cg, sd) however without success.

I would be very gratefull if someone could tell me the possible reasons for such behaviour.

Rafal