Fixing Ionic Positions in Mott-Littleton Calculations

GULP
1

Dear Prof. Gale,

Hope you are well!

TL;DR is it possible to fix region 1 ions in a Mott-Littleton calculation whilst including the response from region 2a in the defect energy? Also, is it possible to output region 2a positions after defect optimisation (or single point ML runs)?

I am currently using GULP to carry out defect calculations using the Mott-Littleton (ML) method. I have a few questions that could use your expertise (or anyone else who is familiar with ML calcs):

I understand GULP can perform a single-point ML run that produces the defect energy of the system without optimising ions in any region. I have confirmed this by generating a dump file and using the option word ‘move2a_to_1’: ions in region 2a sit perfectly on the lattice sites and are not displaced (as they shouldn’t). However, I aim to obtain a defect energy value such that ions in region 1 are not relaxed, and ions in region 2a are displaced in response to how ions are positioned in region 1. Is this possible?

So far, I have tried the following methods to fix ions in region 1:

  • do a single point calc, dump out region 1 ions in restart file, flag them so they are fixed, change the keyword to ‘opti conp defect regi_before’, compare region 1 coordinates before and after the optimisation. I found that inserting flags does not keep region 1 ions fixed.

Sn cor 4.12377 0.00000 0.00000 2.42000 1.0000 0.0000 0 0 0 0 0
Sn she 4.12377 0.00000 0.00000 2.42000 1.0000 0.0000 0 0 0 0 0
(I have also tried 1 1 1 for a quick comparison)

  • Similar to the above, instead of flagging them, I appended ‘fix xyz’ at the end of line describing region 1 atoms from the dump file, and follow the same optimisation protocol described above. I found that this method also did not fix ions in region 1. (The reason why I tried ‘fix xyz’ was I discovered when using serial GULP 6.0, these terms or its variants were added at the end of each line describing ions in region 1, but not the case in parallelised GULP 6.0)

Sn cor 4.12377 0.00000 0.00000 2.42000 1.0000 0.0000 0 0 fix xyz
Sn she 4.12377 0.00000 0.00000 2.42000 1.0000 0.0000 0 0 fix xyz

Am I doing this wrong? Or did I miss something from the theory? Happy to supply input files if needed!

Thank you very much for your time!

Cyril

2

A quick update on the progress after I have posted the thread above:

I have managed to make the flags work by removing ‘conp’ keyword after inspecting the source code. Now starting from an optimised lattice, I managed to fix ionic positions in region 1, but when I use ‘move2a_to_1’ and had a look at the dump file, ions in region 2a do not seem to be displaced, and hence I suspect when a) performing a single point ML calculation, region 2a ions are not displaced and b) performing a defect optimisation with fixed region 1 ions, is equivalent of performing a single point ML calculation where region 2a ions are not displaced.

Consequently, the question remains: is there any method which allow displacements of region 2a ions to be considered in the total defect energy whilst fixing all ions in region 1.

3

Hi Cyril,
The region 2 ions should be relaxed even for a single point calculation of region 1 since the relaxation of region 2 is implicit. For example, if I take example6.gin and run it with “grad conp” in the keywords instead of opti then the output contains:

Components of defect energy :

Region 1 - region 1 = 23.79004558 eV
Region 1 - region 2a (unrelaxed) = -5.28237060 eV
Region 1 - 2a (relaxed - correction) = -1.95513251 eV
Region 1 (Total) = 16.55254247 eV
Region 2a = 0.97706856 eV
Region 2b = -0.63516969 eV

Total defect energy = 16.89444133 eV

Largest displacement in region 2 = 0.0160 Angstroms

As you can see, there is a relaxation correction in the energy and a non-zero displacement for region 2 listed.
Hope this helps,

Julian

4

Hi Julian,

These all make sense now. Thank you very much for the prompt and detailed explanations.

Best
Cyril