How to optimize the bulk lattice

[frsy]

Dear all,
  I'm new here. Now, I can build and calculate one simple bulk. But how can I optimize the lattice of the bulk while keep the fractional coordinates of all atoms fixed? I looked up the manual and foud ATK ONLY can relax the coordinates of the atom. Do I miss something?
  Regards,

Frsy

[Nordland]

Hey frsy and welcome too

The best way to optimize the geometry in your case, is to calculate the equation of states - or simply the total energy as a function of the lattice constant. The minimum of this curve is then the optimized geometry,
and if the crystal, you are looking at, is defined by a single lattice constant, then this is the superior way.

However if the crystal is defined by more than one lattice constant, it might require alot of calculations in order to determine the minimum, and for that I have create my own little smart script for doing so.
This you can read more about here on the forum

Good luck !

P.S Just out of curiosity, what system are you calculating on?

[frsy]

Thank you, Nordland!
I'm trying to optimize rutile. It has two variable lattice lengths so your script is helpful to me.
Regards,

Frsy

[frsy]

Dear Nordland,
   My input file for optimizating rutile bulk is :
###############################
from stress_optimizers import *
import ATK; ATK.setVerbosityLevel(1)
from ATK.KohnSham import *
# Select Bravais lattice and define lattice constants
bravais_lattice = BodyCenteredTetragonal(
    a=4.594000 * Angstrom,
    c=2.959000 * Angstrom
    )

# Define bulk elements and unit cell coordinates
elements = [Titanium, Oxygen]
coordinates = [[ 0.        ,  0.        ,  0.        ],
               [ 0.69520000,  0.30480000,  0.        ]]

# Set up bulk configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice,
    elements,
    fractional_coordinates=coordinates
    )
nlPrint(bulk_configuration)

method = KohnShamMethod(brillouin_zone_integration_parameters = brillouinZoneIntegrationParameters((5,5,7)),
                        basis_set_parameters = basisSetParameters())
opt_conf = calculateOptimizedBulkConfiguration(bulk_configuration, method, runtimeParameters(1))
################################

    I noticed from the output, the initially calculated stress is very large:
# Stresses
# Stress (A) = -277.951887706 eV/Ang**3
# Stress (C) = -277.700861096 eV/Ang**3
    I think my lattices (4.954, 2.959) are not too far from the experimental ones. Is this stress reliable? Or do I make something wrong?
    Regards,

Frsy

[Nordland]

Hey Frsy.

The stress is calculated in a very crude way, and I don't think they should be taken for more than a measurement if the cell is in optimal
configuration or not.

I will try to run your script and see what I get for the optimial structure..

[frsy]

After debug "pseudoStress" and "findFractionalCoordinates" I found incorrect fractional_coordinates were used to build strained_configuration. So the energy difference is large and thus the absolute stress is also large.
Maybe the sub "findFractionalCoordiantes" doesn't work in my case. (BCT Rutile)

My input fractional coords:
coordinates = [[ 0.        ,  0.        ,  0.        ],
               [ 0.69520000,  0.30480000,  0.        ]]

The initial coords in output:
# Index  Element  x (Ang)  y (Ang)  z (Ang)
      0       Ti     0.00     0.00     0.00
      1       O     -0.91     0.91     1.48

Recognized coords in the pseudoStress for the first time:
# Index  Element  x (Ang)  y (Ang)  z (Ang)
      0       Ti     0.00     0.00     0.00
      1       O     -0.91     0.91     1.48

Strained coords after using "findFractionalCoordinates":
# Index  Element  x (Ang)  y (Ang)  z (Ang)
      0       Ti     0.00     0.00     0.00
      1       O      2.07     0.25    -0.74

It is far from the original coords, the calculated fractional coords may be questionable.

[Anders Blom]

One should use the primitiveVectors(), not conventionalVectors(), when you compute the fractional coordinates. If you change that single line in the routine findFractionalCoordinates() I believe it will work as intended, that is line 12 in stress_optimizers.py should be:

Code

uc = configuration.bravaisLattice().primitiveVectors()

[Nordland]

I have reupload the script in a new version to support cells where the conventional cell and the primitive is not the same.

[frsy]

Yes, Blom.
The stress is normal now. Thank you.

[Anders Blom]

Apparently there were two places in the script that needed modifying, so I hope you used the updated script posted by Nordland.