Improved relaxation of bulk systems using approximative stress.

[Nordland]

Hey everyone.

I have created a script for calculating the approximative stress on the unitcell for bulk crystal. I have combined this
with the relaxation module that comes with ATK, and tries to use the same interface as the one official included.

Given a bulk configuration and method to be used, it will first optimized the internal structure of the unitcell according to
the forces on the atoms, and then it will calculate the approximative stress on the unitcell, and scale the unitcell according to
the stress. It will then calculate the optimized internal geometry again, and calculate the stress again and repeat this cycle until
all the forces and stress are below a certain threshold.

Normally the fastest way to calculate the optimized geometry for bulk configuration in ATK is using the equation of state,
however for more complex systems like wurtzite ZnO, it requires alot of calculations! In these cases I think this script is
superior, however for simple systems like Si(FCC) the equation of states is properly better.

However when the goal is to determine the effect of pressure and strain on the unitcell on the band gap etc, it is important that the unitcell
is in equilibrium, and therefore I find this script quite usefull.

The usage is pretty straight forward - download this script and use it in the following way:

from stress_optimizers import *
....
optimized_conf = calculateOptimizedBulkConfiguration(configuration,
                                                                                         method,
                                                                                         runtime_parameters,
                                                                                         geometric_optimization_parameters,
                                                                                         stress_tolerance = 5.0e-3*Units.eV/Units.Ang**3,
                                                                                         maximum_iterations = 100)

• The configuration is the bulk configuration to be optimized.
• The method is the KohnShamMethod to be used for this optimization.
• The runtime parameters for control the verbosity and checkpoint file.
• The geometric_optimization_parameters is the ATK geometric optimization parameters for the forces.
• The stress tolerance thresshold for the optimization - default 5.0e-3 ev/Ang**3
• Maximum iterations is the maximum number of iterations of the unitcell optimization.

Feel free to use.

Updated 1:
Updated the script to correctly support cells where the conventional cell and the primitive cell is not identical.

Updated 2:
Updated with version that are MPI safe and allows the user to constrain atoms.

[Nordland]

I have applied this functionality on Si(fcc) to see how good agreement, you can get by using only the default parameters
of ATK, between the calculated lattice constant and the experimental value.

# ----------------------------------------
# Lattices parameters (FaceCenteredCubic)
# ----------------------------------------
# A = 5.41146687312 Ang

The experimental value is 5.4306 Ang, and therefore the errors is only 0.3%! and that is using only default parameters!
I have attached the script for reproducing these results ( as well the results )

[Nordland]

It is my experience that one of the lattice types that is notorious hard to perform good calculations of the lattices constants on
is the semi-conducting wurtzite crystals. Again I have tryed only with the default parameters to see how good the agreement is between experimental values
and calculation.

# ----------------------------------------
# Lattices parameters (Hexagonal)
# ----------------------------------------
# A = 3.20426776283 Ang
# C = 5.21758120909 Ang
# ----------------------------------------
# Fractional coordinates
# -----------------------------------------
# Index  Element  x (frac)   y (frac)   z (frac)
#    0      Ga       0.333378051061      0.666729233831       0.0040187886993
#    1      Ga       0.666698321157      0.333359730274       0.504110187473
#    2       N        0.333371154227      0.66669569163         0.380881963248
#    3       N        0.66671438049        0.333401955497       0.880886063456

The experimential values are A = 3.18 Ang, C = 5.1660 Ang and U = 0.385.
The errors for A is 0.7%, and the errors for C is 1%,  and for U is also 1%

Again this is done using the default parameters, it can properly be done even better if the parameters are tuned to
give better results for GaN.

[rosen]

Thank Nordland, it is very kind of you. The script is useful.

[ipsecog]

The pseudo-strain method is very powerful! Very nice work, Nordland. I wish this functionality could be included in VNL soon, it would be very valuable!

However, it would be even better to compute the stress analytically via the Hellmann-Feynman relations, and use some standard relaxation tool like QuasiNewton also for the lattice constants (right?). Perhaps this can be introduced in ATK soon?

[Nordland]

I have updated the script in a new version - please feel free to try it out.

[frsy]

Dear Nordland,
    I suggest you add constaints parameter to the calculateOptimizedBulkConfiguration() so that one can fix atoms in the bulk.
Another suggestion is to add "if processIsMaster():" before print(). The current script messes up the screen in a parallel calculation.

    Regards,

Frsy

[Nordland]

Hey Frsy!

That are two good ideas!

I will update the script accordingly and feel free to test it....

[zdhlover]

 
This method is powerful and beautiful, but I have a problem about the relaxation of bulk system : I have a supercell bulk system, and I want to optimize only the C lattice constant and also the atom coordinates(may be we also need constrain on atom in some case) in series as this: optimize C then used the optimized C optimize the atom coordinates, and repeated the above two steps , ...... , until the C lattice constant and the atom coordinates are all well optimized ,but the a lattice constant and b lattice constant are all fixed in the whole process.

Could someone help me? Thanks a lot!

[Nordland]

Okay. Right now it optimizes the coordinates of the atoms at each steps, so that point is already covered.

Let me see if I can find a way to limit to only one of the lattice constants at the time