reading lattice parameters from .castep

[AndyDuff1]

I discovered an issue when using a castep.castep file for which lattice parameters have been held fixed. Normally of course when using CASTEP's internal geometry optimization this wouldn't be the case. However if one performed the geometry optimization instead using a grid-based approach (a grid of different lattice parameters), and then did a final run relaxing atom positions but keeping lattice vectors fixed, then one might have a .castep file like this.

Currently the code assumes the .castep corresponds to a geometry optimization run. I.e., 'Unit Cell' is looked for after 'Final Configuration:'. However if a grid-based approach is used there is no 'Unit Cell' after 'Final Configuration'. Of course in this case any 'Unit Cell' will do, as they are all the same (unchanged over the simulation).

I didn't see any guidance on contributions to this code so I've included a fix for this here. If you agree with the change perhaps you could incorporate it into your codebase?:

Change to castep.py:

After:
dotcastep_latt_RE = ...
add:
#New regular expression targeting just the 'Unit Cell' section
dotcastep_latt_fixed_RE = re.compile(r"""\s+Unit\sCell\s*\n\s+-+\s*\n
\s+Real\sLattice(A)\s+Reciprocal\sLattice(1/A)\s*\n
\s+([+-]?\d+.\d+)\s+([+-]?\d+.\d+)\s+([+-]?\d+.\d+)\s+([+-]?\d+.\d+)\s+([+-]?\d+.\d+)\s+([+-]?\d+.\d+)\s*\n
\s+([+-]?\d+.\d+)\s+([+-]?\d+.\d+)\s+([+-]?\d+.\d+)\s+([+-]?\d+.\d+)\s+([+-]?\d+.\d+)\s+([+-]?\d+.\d+)\s*\n
\s+([+-]?\d+.\d+)\s+([+-]?\d+.\d+)\s+([+-]?\d+.\d+)\s+([+-]?\d+.\d+)\s+([+-]?\d+.\d+)\s+([+-]?\d+.\d+)\s*\n""", re.VERBOSE)

The start of parse_dotcastep should then read:

def parse_dotcastep(seedname):
"""
Extract lattice and atom positions from a .castep
file. List of atoms may be empty (e.g. MgO)
"""
dotCastep = open(seedname+".castep","r")
# Find the lattice. First attempt to locate 'Final Configuration' followed by 'Unit Cell'
latticeblock = dotcastep_latt_RE.findall(dotCastep.read()) # Get the last block - handle concat restarts
if latticeblock:
print("Lattice block identified after 'Final Configuration'")
latticeblock = latticeblock[-1]
else:
print("No 'Unit Cell' after 'Final Configuration' - assuming fixed lattice relaxation.")
# If 'Unit Cell' does not proceed 'Final Configuration', unit cell parameters may have been fixed (using
# FIX_ALL_CELL in .cell file). This may occur e.g. if a grid method was used to optimize geometry followed
# by a fixed lattice parameter relaxation. In this case rewind and look instead just for 'Unit Cell'
dotCastep.seek(0)
latticeblock = dotcastep_latt_fixed_RE.findall(dotCastep.read())[-1]
lattice = []
...

[andreww]

I've not thought it through very carefully, but I'm happy to look at a patch for this (fork the repo, commit a change on a branch on your fork, and submit a pull request). However, I think you would need to be careful from the point of view of the science being valid. I think there needs to be a check at the analysis stage that the fitted stress-strain relationship (for all Cijs) gives zero stress for zero strain (probably for the case of a grid of fixed lattice constants case and for the cell relaxation case) akin to the check in castep that the atoms are at a local minimum before it'll do a phonon calculation. We could also check that the relationship is linear I suppose.

[AndyDuff1]

Okay thanks for the reply - I've made a pull request.