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!
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.
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,
If there is a point defect in region 1, ions in region 2a should have a harmonic response. Is there any way to output the positions of the displaced region 2a ion at the end of Mott-Littleton calculation? The use of move_2a_to_1 gave the region 2a ions on perfect lattice sites which does not make sense.
An example input, output and dump file are attached. As can be seen in the dump file, atoms in region 2a are sitting at the perfect lattice site 1.res (48.8 KB) input.gin (699 Bytes) input.gout (79.6 KB)
I have also gone into the source code move2a1.F90 and saw the following line in the dev log:
! 3/13 Displacements not applied in region 2 as in previous
! version since this causes problems with restarting.
Is there any ways to turn the displacement back on for region 2a ions?
Apologies for this issue. Looking it to it (since it’s been a while), prior to 2013 the possibility to modify the coordinates of region 2a atoms to include the displacement when transferred to region 1 existed (as you found in the source code). Unfortunately the change to not include the displacements wasn’t reflected in the help text & so I’ll correct this. I can probably add the option to include the displacement back as a modifier to move2a1 with a warning note that this may cause issues, but that experienced users can try it if they like. I’ll try to do this and post as an update to version 6.3 in the next day or so.
Regards,
Julian
Just a quick update. I’ve added the sub-option to include the displacements in version 6.3.2 (also 6.3 main trunk). If you update to this version and then include the line;
move2a1 displace
it should hopefully work for you. Help documentation is also updated.