In a previous tutorial, we have presented how to determine the structural ground state of the compound CaO between the Crystal structures NaCl (B1) and CsCl (B2) (How to calculate the structural ground state of Calcium Oxide CaO). In this tutorial, we will present how to calculate the transition pressure from The phase B1 to the phase B2.

We will use a shell script for calculating the Enthalpy versus the Pressure for each phase and another shell script for plotting.

We need to use the following input files:

CaO.param

task     GEOMETRYOPTIMISATION ! The TASK keyword instructs CASTEP what to do
xc_functional   PBE           ! Which exchange-correlation functional to use.
cut_off_energy  500 eV        !
opt_strategy    speed         ! Choose algorithms for best speed
calculate_stress : true

 

CaO-B1_init.cell

%BLOCK LATTICE_ABC
ang    # angstrom units
   3.42 3.42 3.42
   60 60 60
%ENDBLOCK LATTICE_ABC

%BLOCK POSITIONS_FRAC
Ca   0.000000000000000   0.000000000000000   0.000000000000000 
O    0.500000000000000   0.500000000000000   0.500000000000000 
%ENDBLOCK POSITIONS_FRAC

kpoint_mp_grid 10 10 10

symmetry_generate

 

CaO-B2_init.cell

%BLOCK LATTICE_ABC
ang    # angstrom units
   2.90 2.90 2.90
   90 90 90
%ENDBLOCK LATTICE_ABC

%BLOCK POSITIONS_FRAC
Ca   0.000000000000000   0.000000000000000   0.000000000000000 
O    0.500000000000000   0.500000000000000   0.500000000000000 
%ENDBLOCK POSITIONS_FRAC

kpoint_mp_grid 10 10 10

symmetry_generate

 

tail.cell


%BLOCK EXTERNAL_PRESSURE
GPa
    @P@    0.0000000000    0.0000000000
                    @P@    0.0000000000
                                    @P@
%ENDBLOCK EXTERNAL_PRESSURE

 

Script cal-enthalpy.sh

#!/bin/sh
# Developped by Dr. Abderrahmane Reggad
# Email: abde.reggad@gmail.com

seed=CaO

rm -f ${seed}.cell
rm -f ${seed}.castep
rm -f ${seed}_tmplt.cell

var='B1 B2'
for stsy in $var; do

cat ${seed}-${stsy}_init.cell tail.cell >> ${seed}_tmplt.cell

for PR in `seq -w 55 1 70`
do
  
  sed -e "s/@P@/$PR/g" ${seed}_tmplt.cell > $seed.cell
  
  echo " running structure ${stsy} "
  echo " running Pressure= $PR GPa "
  mpirun -np 2 castep.mpi ${seed}
  
#analyse the runs and extract the values we need into a single results file
grep 'Pressure' ${seed}.castep| tail -n 1  | awk '{print $3}' > P
grep 'Final Enthalpy' ${seed}.castep| tail -n 1  | awk '{print $5}' > H
paste P H >> ${seed}-${stsy}_PH.dat
rm P H
done

echo 'finished with results in '${seed}-${stsy}_PH2.dat
echo

done

 

Script plot-HP.sh

#!/bin/sh
# Developped by Dr. Abderrahmane Reggad
# Email: abde.reggad@gmail.com

seed=CaO
rm -f fit_Burch.py

var='B1 B2'
for stsy in $var; do

cat >> fit_Burch.py <<EOF
#!/usr/bin/python3
# Matt Probert 21/02/2017

# simple script to load a text file of P (GPa) vs H (eV)
from scipy.optimize import leastsq
import numpy as np
import matplotlib.pyplot as plt
import sys

#get HP filename and load data - file must be 2 columns containing pressure and enthalpy
pressures, enthalpies = np.loadtxt('${seed}-${stsy}_PH.dat',unpack=True)
col = (np.random.random(), np.random.random(), np.random.random())
#plot the raw data and fitted curve on top to file as a cross-check
plt.plot(pressures,enthalpies, 'r.')

x = np.linspace(min(pressures), max(pressures), 50)
y = np.linspace(min(enthalpies), max(enthalpies), 50)
plt.plot(x, y, c=col, label='${seed}-${stsy}')

plt.xlabel('Presuure')
plt.ylabel('Enthalpy')
plt.title('Enthalpy VS Pressure ')
plt.legend(loc='best')
plt.savefig('BM_curve.png')

EOF

# run python script
python3 fit_Burch.py 
done

 

Execution

$ sh cal-enthalphy.sh

We got the files CaO-B1_PH.dat    CaO-B2_PH.dat

 

$ sh plot-HP.sh

We got the following plot


Now We need to do another calculation for pressure values between 55 GPa and 70 GPa and for that we need to modify in the shell script the following line

for PR in `seq -w 0 10 100`

by the following line

for PR in `seq -w 55 1 70`

Before doing execution, we need to rename the dat files for avoiding their overwriting.

And to re-plotting the results to get the following picture


From the picture we can see that the transition picture equals to 67 GPa. This value is close to the value found in the following paper


Reference: https://www.researchgate.net/figure/The-enthalpy-difference-between-the-B1-phase-and-B2-phase-as-a-function-of-pressure_fig6_230962333