andreww / elastic-constants Goto Github PK

What happens if there if there were pre-existing constrains in the cell (i.e. a cell_constrains block)? Should they should be removed when creating the new .cell files with the conflicting keyword FIX_ALL_CELL or is this meaningless (as the zero strain cell won't be at energy minimum).

Should automate the test cases, and tell it not to check results to arbitrary precision (5-8 dp is enough).

As reported by email - this needs investigation but I'm adding it here so I don't forget

The version check at line 254 of elastics.py is not robust on recent systems (where scipy reports version strings line '0.12.0'. Should probably just calculate the error ourselves everywhere, or check our calculated error against scipy's and issue a warning if they differ.

If the user supplies a lattice type to generate_strain (with the -l option) the code fails with an undefined variable error. This is the "suptype" variable, that should be user-settable. See modifications in Josh's tree as a starting point.

It does not run with python 3! Shouldn't be too hard to fix (just print statements, I think).

Seems we can now get CASTEP to create a cell file - worth investigating if this is useful (we currently parse this information from the .castep output). See:

By the way, the way I would normally recommend to extract an optimized geometry from a castep run is to set the parameter

write_cell_structure : true

which will write a new .cell file with the co-ordinates at high precision. By comparison "castep2cell" loses precision as it can only
extract the number of significant digits printed in the .castep file.

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 = []
...