A guide for efficiently using LatticeFinder to obtain the optimal lattice constants

In this article, we will describe how to use LatticeFinder to find the optimal lattice constants for a 2D or 3D system. We will describe how to efficiently use LatticeFinder for systems that contain either one lattice constant or two lattice constants.

2D and 3D Systems containing one lattice constant

1) Performing a broad overview of lattice constants

To begin it is often useful to perform a broad overview of lattice constants to determine what lattice constants to focus on. For example, for determining the lattice constant of a Au face-centred cubic (FCC) lattice using the RGL potential, with parameters from Baletto et al. (DOI: 10.1063/1.1448484). Here, we first look at lattice constants between 3.0 Å and 5.0 Å in increments of 0.1 Å. This is performed with the following Run_LatticeFinder.py script:

Run_LatticeFinder.py
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
from LatticeFinder import LatticeFinder_Program

symbol = 'Au'
lattice_type = 'FaceCenteredCubic'

lattice_constant_parameters = (3.0,5.0,0.1)

from asap3.Internal.BuiltinPotentials import Gupta
# Parameter sequence: [p, q, a, xi, r0]
r0 = 4.07/(2.0 ** 0.5)
Au_parameters = {'Au': [10.53, 4.30, 0.2197, 1.855, r0]} # Baletto
cutoff = 8
calculator = Gupta(Au_parameters, cutoff=cutoff, debug=False)

size=(16,16,16)

directions=[]
miller=[]

limits = None

no_of_cpus = 2

LatticeFinder_Program(symbol, lattice_type, lattice_constant_parameters, calculator, size=size, directions=directions, miller=miller, limits=limits, no_of_cpus=no_of_cpus)

This gives an energy per atom vs lattice constant plot as below:

Interpolation Scheme Plot

From this plot, we see that the minimum centres about 4.1 Å. Therefore, we will want to next perform a comprehensive search of the optimum lattice constant between 4.0 Å and 4.2 Å at the least. This is because we set this script to scan the lattic constants in increments of 0.1 Å, therefore the precise lattice constant will lie within the 4.1 ± 0.1 Å range.

2) Performing a comprehensive scan of lattice constants across a small range

We would like to perform a indepth scan of lattice constants between 4.0 Å and 4.2 Å. We will do this by looking between 4.0 Å and 4.2 Å in increments of 0.001 Å. We can add this indepth scan to our original LatticeFinder calculation by changing the Run_LatticeFinder.py script to the following:

Run_LatticeFinder.py
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
from LatticeFinder import LatticeFinder_Program
import numpy as np

symbol = 'Au'
lattice_type = 'FaceCenteredCubic'

lattice_constant_parameters = list(np.arange(3.0,4.0,0.1))+list(np.arange(4.0,4.2,0.001))+list(np.arange(4.2,5.01,0.1))

from asap3.Internal.BuiltinPotentials import Gupta
# Parameter sequence: [p, q, a, xi, r0]
r0 = 4.07/(2.0 ** 0.5)
Au_parameters = {'Au': [10.53, 4.30, 0.2197, 1.855, r0]} # Baletto
cutoff = 8
calculator = Gupta(Au_parameters, cutoff=cutoff, debug=False)

size=(16,16,16)

directions=[]
miller=[]
no_of_cpus=1

LatticeFinder_Program(symbol, lattice_type, lattice_constant_parameters, calculator, size=size, directions=directions, miller=miller, no_of_cpus=no_of_cpus)

You can change the lattice_constant_parameters input from

6
lattice_constant_parameters = (3.0,5.0,0.1)

to

7
lattice_constant_parameters = list(np.arange(3.0,4.0,0.1))+list(np.arange(4.0,4.2,0.001))+list(np.arange(4.2,5.01,0.1))

and rerun LatticeFinder. LatticeFinder will perform calculations with lattice constants that you have not already obtained. Once you rerun LatticeFinder on your script again, you will get the following energy per atom vs lattice constant plot as below:

This gives an energy per atom vs lattice constant plot as below:

Interpolation Scheme Plot

The data from this is shown in results_file.txt

Symbol: Au
Lattice_type: FaceCenteredCubic
calculator: <asap.RGL object at 0x0x54b110>
size: (16, 16, 16)
directions: []
miller: []
Lattice Constant Parameters: ['c']

Properties of System: 

Total energy: -62720.92122649953 eV
No. of atoms: 16384 Atoms (Note the number of atoms along each natural direction of the bulk is (16, 16, 16))
Cohesive energy: -3.8281812272033404 eV/Atom

Total Volume: 276148.8097280001 Angstroms^3
Volume per atom: 16.854785750000005 Angstroms^3/Atom

Stress tensor:
[[1.29870561e-03 1.29377728e-19 1.68652768e-19]
 [1.29377728e-19 1.29870561e-03 1.29962170e-19]
 [1.68652768e-19 1.29962170e-19 1.29870561e-03]]

Bulk Modulus: 218.1717839967987 GPa

How to obtain the Bulk Modulus

In your results_file.txt you will also be given a value for the bulk modulus. If you want the bulk modulus with a good amount of accurancy, you often need to perform an indepth scan of lattice constants across a good range of data points. Furthermore, the range of lattice constants given must be regular, i.e. incremental.

For example, we have obtain the energies of an Au FCC crystal for lattice constants between 3.6 Å and 4.6 Å in increments of 0.001 Å. This is given in the Run_LatticeFinder.py script below:

Run_LatticeFinder.py
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
from LatticeFinder import LatticeFinder_Program

symbol = 'Au'
lattice_type = 'FaceCenteredCubic'

lattice_constant_parameters = (3.6,4.6,0.001)

from asap3.Internal.BuiltinPotentials import Gupta
# Parameter sequence: [p, q, a, xi, r0]
r0 = 4.07/(2.0 ** 0.5)
Au_parameters = {'Au': [10.53, 4.30, 0.2197, 1.855, r0]} # Baletto
cutoff = 8
calculator = Gupta(Au_parameters, cutoff=cutoff, debug=False)

size=(16,16,16)

directions=[]
miller=[]
no_of_cpus=1

LatticeFinder_Program(symbol, lattice_type, lattice_constant_parameters, calculator, size=size, directions=directions, miller=miller, no_of_cpus=no_of_cpus)

This gives the following energy per atom vs lattice constant plot:

Interpolation Scheme Plot

and the data is given in results_file.txt

Symbol: Au
Lattice_type: FaceCenteredCubic
calculator: <asap.RGL object at 0x0x54b110>
size: (16, 16, 16)
directions: []
miller: []
Lattice Constant Parameters: ['c']

Properties of System: 

Total energy: -62720.92122649953 eV
No. of atoms: 16384 Atoms (Note the number of atoms along each natural direction of the bulk is (16, 16, 16))
Cohesive energy: -3.8281812272033404 eV/Atom

Total Volume: 276148.8097280001 Angstroms^3
Volume per atom: 16.854785750000005 Angstroms^3/Atom

Stress tensor:
[[1.29870561e-03 1.29377728e-19 1.68652768e-19]
 [1.29377728e-19 1.29870561e-03 1.29962170e-19]
 [1.68652768e-19 1.29962170e-19 1.29870561e-03]]

Bulk Modulus: 218.1717839967987 GPa

We can see how well the lattice constant and been modelled using a energy vs volume plot given by the EOS module in ASE

Interpolation Scheme Plot

See Equation of state (EOS) in ASE for more information.

2D and 3D Systems containing two lattice constant

1) Performing a broad overview of lattice constants

Like in the one lattice constant case, it is often useful to perform a broad overview of lattice constants to determine what lattice constants to focus on. For example, for determining the lattice constant of a Au hexagonal closed packed (HCP) lattice using the RGL potential, with parameters from Baletto et al. (DOI: 10.1063/1.1448484). Here, we first look at lattice constants for c between 2.0 Å and 5.0 Å in increments of 0.1 Å, and a between 3.0 Å and 6.0 Å in increments of 0.1 Å. This is performed with the following Run_LatticeFinder.py script:

Run_LatticeFinder.py
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
from LatticeFinder import LatticeFinder_Program

symbol = 'Au'
lattice_type = 'HexagonalClosedPacked'

lattice_constant_parameters = {'a': (2.0,5.0,0.1), 'c': (3.0,6.0,0.1)}

from asap3.Internal.BuiltinPotentials import Gupta
# Parameter sequence: [p, q, a, xi, r0]
r0 = 4.07/(2.0 ** 0.5)
Au_parameters = {'Au': [10.53, 4.30, 0.2197, 1.855, r0]} # Baletto
cutoff = 8
calculator = Gupta(Au_parameters, cutoff=cutoff, debug=False)

size_single = 28
size=(size_single,size_single,size_single)

directions=[]
miller=[]

limits = {'a': (2.6,3.2), 'c': (4.4,5.0)}
no_of_cpus = 8


LatticeFinder_Program(symbol, lattice_type, lattice_constant_parameters, calculator, size=size, directions=directions, miller=miller, limits=limits, no_of_cpus=no_of_cpus)

This gives an energy per atom vs lattice constant plots as below. We have use the limits section to zoom in on the important area of the potential energy surface to make it easier to see the bottom of the well.

Interpolation Scheme Plot Interpolation Scheme Plot

From this plot, we see that the minimum centres about a = 2.9 Å and c = 4.7 Å. Therefore, we will want to next perform a comprehensive search of the optimum lattice constant for a between 2.8 Å and 3.0 Å and c between 4.6 Å and 4.8 Å, at the least. This is because we set this script to scan the lattic constants in increments of 0.1 Å, therefore the precise lattice constant will lie within the a = 2.9 ± 0.1 Å and c = 4.7 ± 0.1 Å ranges.

2) Performing a comprehensive scan of lattice constants across a small range

We would like to perform a indepth scan of lattice constants between for a between 2.8 Å and 3.0 Å and c between 4.6 Å and 4.8 Å in increments of 0.001 Å. Unlike for LatticeFinder calculation for systems with one lattice constant, you can not add data easily to the plots made for systems with two lattice constants and therefore it is best to perform a new LatticeFinder calculation for your system. The Run_LatticeFinder.py file for running this is shown below:

Run_LatticeFinder.py
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
from LatticeFinder import LatticeFinder_Program

symbol = 'Au'
lattice_type = 'HexagonalClosedPacked'

lattice_constant_parameters = {'a': (2.8,3.0,0.001), 'c': (4.6,4.8,0.001)}

from asap3.Internal.BuiltinPotentials import Gupta
# Parameter sequence: [p, q, a, xi, r0]
r0 = 4.07/(2.0 ** 0.5)
Au_parameters = {'Au': [10.53, 4.30, 0.2197, 1.855, r0]} # Baletto
cutoff = 8
calculator = Gupta(Au_parameters, cutoff=cutoff, debug=False)

size_single = 28
size=(size_single,size_single,size_single)

directions=[]
miller=[]

limits = {'a': (2.6,3.2), 'c': (4.4,5.0)}
no_of_cpus = 8


LatticeFinder_Program(symbol, lattice_type, lattice_constant_parameters, calculator, size=size, directions=directions, miller=miller, limits=limits, no_of_cpus=no_of_cpus)

You can change the lattice_constant_parameters input from

6
lattice_constant_parameters = {'a': (2.0,5.0,0.1), 'c': (3.0,6.0,0.1)}

to

6
lattice_constant_parameters = {'a': (2.8,3.0,0.001), 'c': (4.6,4.8,0.001)}

and rerun LatticeFinder. LatticeFinder will perform calculations with lattice constants that you have not already obtained. Once you rerun LatticeFinder on your script again, you will get the following energy per atom vs lattice constant plot as below:

This gives an energy per atom vs lattice constant plots as below:

Interpolation Scheme Plot Interpolation Scheme Plot

The data from this is shown in results_file.txt