Fitting a Coulomb-Buckingham force field for solid electrolytes

GULP
1

Hi,
I am trying to fit an interatomic potential of the Coulomb-Buckingham form for Li3PS4, a typical solid electrolyte, by targeting the elastic properties of one of its known stable phases: the beta phase. The reference structure provided in the input corresponds to the minimum energy structure derived with the quantum espresso code and the full elastic tensor Cij (measured in the experiments) has been provided as target property. The Coulomb charges have been set equal to the formal charges, i.e. + 1 for Li, +5 for P and -2 for S and the initial guesses for the Buckingham potential correspond to the known parametrization of a different electrolyte (Li3ClO).

I have therefore produced the following input file, for a Li12P4S16 structure. The fitting converges to a parametrization with a very large sum of squares (~60000) and the agreement with the observables given as reference is very poor: for instance C11 is predicted to be 337 GPa, as compared to the given 34 GPa. Modifying the initial guess for the Buckingham parameters within a quite wide range does not solve the problem and it seems that the BFGS method always converges to this “very-far off” minimum. I am new to the use of GULP, so:

  1. do you see specific issues with the procedure that I have followed and/or with the current input file that I have created, such as, e.g., too many fitting parameters and too few observables?
  2. I know that a Coulomb-Buckingham potential should provide a correct description of this material as it is described very well in the literature for other electrolytes with similar features; so, what initial guess for the Buckingham parameters would one start with in the absence of any previous knowledge?

Any advice on this would be very helpful.
Thanks.
Best regards,
Lorenzo Gigli

fit conp simul opti
cell
6.18608204081008 8.09698016944516 13.1261169247833 90.000000 90.000000 90.0000 0 0 0 0 0 0
fractional
Li 0.389032 0.0318183 0.667829
Li 0.610967 0.96818 0.33217
Li 0.889031 0.96818 0.832169
Li 0.110968 0.0318184 0.167829
Li 0.610967 0.531818 0.33217
Li 0.389032 0.468181 0.667829
Li 0.110968 0.468181 0.16783
Li 0.889031 0.531818 0.832169
Li 0.499999 4.90626e-09 -6.1819e-08
Li -2.80902e-07 -1.59742e-07 0.499999
Li 0.499999 0.499999 -7.70879e-08
Li -1.68722e-07 0.499999 0.499999
P 0.153742 0.25 0.910995
P 0.846256 0.749999 0.0890036
P 0.653742 0.749999 0.589004
P 0.346257 0.25 0.410995
S 0.826102 0.249999 0.898712
S 0.173897 0.749999 0.101287
S 0.326101 0.749999 0.601286
S 0.673897 0.249999 0.398712
S 0.27273 0.249999 0.0591669
S 0.727269 0.749999 0.940832
S 0.77273 0.749999 0.440831
S 0.227269 0.249999 0.559167
S 0.273039 0.037691 0.847452
S 0.726959 0.962308 0.152547
S 0.773039 0.962308 0.652546
S 0.22696 0.0376906 0.347452
S 0.726959 0.53769 0.152547
S 0.273039 0.462308 0.847452
S 0.22696 0.462309 0.347452
S 0.773039 0.53769 0.652546
species
Li core 1.0
P core 5.0
S core -2.0

observables
elastic 21
1 1 34.0
1 2 16.0
1 3 11.0
2 2 50.0
2 3 17.0
3 3 31.0
4 4 10.0
5 5 12.0
6 6 12.0
1 4 0.0
1 5 0.0
1 6 0.0
2 4 0.0
2 5 0.0
2 6 0.0
3 4 0.0
3 5 0.0
3 6 0.0
4 5 0.0
4 6 0.0
4 7 0.0
end

buck
S core S core 10000.0 0.20 15.0 1 1 1
buck
P core P core 10000.0 0.20 15.0 1 1 1
buck
Li core Li core 360.5269 0.16098 0.001 1 1 0
buck
Li core S core 400.0 0.30 0.001 1 1 0
buck
Li core P core 400.0 0.30 0.001 1 1 0
buck
P core S core 0.0 1.0 0.001 0 0 0

2

Hi Lorenzo,
There are a couple of comments I can make here:

  1. Formal charge models will nearly always overestimate elastic constants since the reality is that the atoms aren’t formally charged and the electrostatics is a major contribution to the energy. This is especially true when you try to have P with +5. Using a shell model would partially offset this (for lower symmetry environments) since the polarisation can to counter the excessive charges. Alternative you need to look at partial charges.
  2. You’re using “simul” as a keyword but this only applies to shell model fits and you don’t have shells. Once you have a guess at the potentials that gives a stable optimisation then you should try to use “relax” as a keyword since this will give the best fit as the properties will be computed at the optimised structure rather than the input one (where they may not be meaningful). See the Phil. Mag. B GULP paper on fitting for details of these points and general information on fitting potentials.
    Regards,
    Julian
3

Hi Julian,
thanks, this is very helpful.

Best regards,
Lorenzo