Bash script to extract lattice parameters from CONTCAR

Question

I am using VASP to calculate material properties. I want to know how to extract the lattice parameters from my cells automatically with a script. I know about the method with VESTA (described here How to extract lattice constants a and c from relaxed hexagonal structure using VASP?), but doing this manually for every composition takes way too much time.

Since I am familiar with bash from running my calculations I would like to use a bash-script that calculates a, b, c, α, β, and γ from the relaxed CONTCAR file that I can loop over my directories. Is this possible to do in bash or do I need to resort to python or something similar?

Edit:

as requested, an example would be this CONTCAR file. Since there are 64 individual atoms and the the relevant lines are 2-5 in all cases, I am omitting the lines at the bottom. This is a supercell, which means the lattice parameters are multiplied, but the script should not need to address that.

2x2x2 supercell ternary nitride         
1.00000000000000     
 8.3651760123584129    0.0104817210801190   -0.0093895223039569
 0.0104823789019837    8.3504089378499629   -0.0085602185310777
-0.0093916727121912   -0.0085617113679278    8.3470735253222390
Al   Ti   N 
16    16    32
Direct
0.0003223292397248  0.9972239351012413  0.0011915008883630
0.4991700057920037  0.9984995833496634  0.5000934141482651
...

A preferred output for this cell would then be

 8.36519 8.35042 8.34708 90.1176 90.1289 89.8562

for a, b, c, α, β, and γ. Separation by tabs, commas, or semicolons would also be acceptable.

Answer 1

I'd recommend ASE, a ready solution not to re-invent the wheel and ensure all the boundary cases are properly handled. Also this code will be easier to maintain in future:

#!/usr/bin/env python

import sys

from ase.io.vasp import read_vasp
from ase.geometry import cell_to_cellpar


try:
    ase_obj = read_vasp(sys.argv[1])
except IndexError:
    sys.exit('A CIF path must be given!')

print("%.5f %.5f %.5f %.4f %.4f %.4f" % tuple(
    cell_to_cellpar(ase_obj.cell)
))

If you want to extract rather the conventional cell, not primitive, try to apply spglib in between:

#!/usr/bin/env python

import sys

from ase.io.vasp import read_vasp
from ase import Atoms
from ase.geometry import cell_to_cellpar
import spglib


try:
    ase_obj = read_vasp(sys.argv[1])
except IndexError:
    sys.exit('A CIF path must be given!')

cell, scaled_pos, numbers = spglib.standardize_cell(ase_obj, to_primitive=False)

ase_obj = Atoms(
    scaled_positions=scaled_pos,
    numbers=numbers,
    cell=cell,
    pbc=True
)
print("%.5f %.5f %.5f %.4f %.4f %.4f" % tuple(
    cell_to_cellpar(ase_obj.cell)
))

Answer 2

Of course, you can use python. For example:

• CONTCAR template:

  GeTe
  1.00000000000000
  4.2642961037835692   -0.0182827012387123   -0.0127716218381237
  2.2215279123437899    3.6399683038796313   -0.0127716218406244
  2.2215279123442269    1.2405016904345516    3.4220882400997064
  Ge   Te
  1     1
  Direct
  -0.0008620876194062 -0.0008620876194055 -0.0008620876194060
   0.5218621726194035  0.5218621076194079  0.5218621136194089
  
   0.00000000E+00  0.00000000E+00  0.00000000E+00
   0.00000000E+00  0.00000000E+00  0.00000000E+00
  

• Lattice parameters from VESTA:

  a =  4.26435 Å      α = 58.8490° 
  b =  4.26435 Å      β = 58.8490°
  c =  4.26435 Å      γ = 58.8490°
  

• Python script:

  import numpy as np
  
  with open("CONTCAR","r") as f:
      lines=f.readlines()
  
  lat_mat=np.zeros((3,3))
  for i in range(3):
      lat_mat[i,:]=list(map(float,lines[2+i].strip().split()))
  
  a=np.sqrt(np.sum(lat_mat[0,:]**2))
  b=np.sqrt(np.sum(lat_mat[1,:]**2))
  c=np.sqrt(np.sum(lat_mat[2,:]**2))
  print(a,b,c)
  
  alpha = np.arccos(np.dot(lat_mat[1,:],lat_mat[2,:])/b/c)*360/2/np.pi
  beta  = np.arccos(np.dot(lat_mat[0,:],lat_mat[2,:])/a/c)*360/2/np.pi
  gamma = np.arccos(np.dot(lat_mat[0,:],lat_mat[1,:])/a/b)*360/2/np.pi
  print(alpha,beta,gamma)
  

• Python output:

  4.26435442150775   4.264354421585447  4.264354421518477
  
  58.84900754415152  58.84900754404418  58.849007544603154
  

You may realize this with a bash script. Also, the atomic simulation environment (ASE) can do this.

Hope it helps.