diff --git a/castep.py b/castep.py index 8908d92..e3df63a 100755 --- a/castep.py +++ b/castep.py @@ -9,7 +9,10 @@ """ import re -import scipy as S +import numpy as np +import logging + +logger = logging.getLogger(__name__) version = 0.1 @@ -32,35 +35,42 @@ dotcastep_poinggroup_RE = re.compile(r"^\s+Point group of crystal =\s+([\+\-]?\d+):") 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 - latticeblock = dotcastep_latt_RE.findall(dotCastep.read())[-1] # Get the last block - handle concat restarts - lattice = [] - lattice.append([float(latticeblock[0]), float(latticeblock[1]), float(latticeblock[2])]) - lattice.append([float(latticeblock[6]), float(latticeblock[7]), float(latticeblock[8])]) - lattice.append([float(latticeblock[12]), float(latticeblock[13]), float(latticeblock[14])]) - # rewind and search for and final atomic positions (these will be absent if e.g. they are all on symmetry positions) - dotCastep.seek(0) - in_atoms = False - pointgroup = None - atoms = [] - for line in dotCastep: - sym_line = dotcastep_poinggroup_RE.search(line) - atom_line = dotcastep_atomline_RE.search(line) - if (in_atoms and atom_line): - atoms.append([atom_line.group(1), float(atom_line.group(2)), \ - float(atom_line.group(3)), float(atom_line.group(4))]) - elif ((not in_atoms) and (dotcastep_infinal_RE.search(line))): - in_atoms = True - elif (sym_line): - pointgroup = int(sym_line.group(1)) - - dotCastep.close() - return (lattice, pointgroup, atoms) + """ + Extract lattice and atom positions from a .castep + file. List of atoms may be empty (e.g. MgO) + """ + with open(seedname+".castep","r") as dotCastep: + # Find the lattice + logger.debug("Find lattice block.") + latticeblock = dotcastep_latt_RE.findall(dotCastep.read())[-1] # Get the last block - handle concat restarts + logger.debug("Found lattice block '%s'.", latticeblock) + lattice = [] + lattice.append([float(latticeblock[0]), float(latticeblock[1]), float(latticeblock[2])]) + lattice.append([float(latticeblock[6]), float(latticeblock[7]), float(latticeblock[8])]) + lattice.append([float(latticeblock[12]), float(latticeblock[13]), float(latticeblock[14])]) + # rewind and search for and final atomic positions (these will be absent if e.g. they are all on symmetry positions) + dotCastep.seek(0) + in_atoms = False + pointgroup = None + atoms = [] + logger.debug("Find atoms block.") + for line in dotCastep: + sym_line = dotcastep_poinggroup_RE.search(line) + atom_line = dotcastep_atomline_RE.search(line) + if (in_atoms and atom_line): + logger.debug("Found atom line '%s'.", line) + logger.debug("Append atom match '%s, %s ,%s, %s'.", atom_line.group(1), atom_line.group(2), atom_line.group(3), atom_line.group(4)) + atoms.append([atom_line.group(1), float(atom_line.group(2)), \ + float(atom_line.group(3)), float(atom_line.group(4))]) + elif ((not in_atoms) and (dotcastep_infinal_RE.search(line))): + logger.debug("Found beginning of atoms block at '%s'.", line) + in_atoms = True + elif (sym_line): + logger.debug("Found symmetry line '%s'.", line.strip()) + logger.debug("Append symmetry match '%s'.", sym_line.group(1)) + pointgroup = int(sym_line.group(1)) + + return (lattice, pointgroup, atoms) # Regular expressions to match a lattice block in a castep .cell file. Note that these @@ -71,49 +81,55 @@ def parse_dotcastep(seedname): dotcell_atoms_end_RE = re.compile(r"^\s*%ENDBLOCK\s+POSITIONS_(?:FRAC|ABS)", re.IGNORECASE) def produce_dotcell(seedname, filename, defcell, atoms): - """ - produce_dotcell: reads .cell (CASTEP cell file - and writes a new .cell file to replacing the - lattice block with a new crystalographic lattice - (which should be supplied as a list of three lists, each with - three elements). Also adds command to fix cell during optimization. - """ - in_lattice = False - in_atoms = False - have_atoms = (atoms != []) # If we have an empty list, no atoms were optimized so just leave them in the .cell file. - inputfile = open(seedname+".cell", "r") - outputfile = open(filename, "w") - for line in inputfile: - if (dotcell_lattice_end_RE.search(line) and in_lattice): - in_lattice = False - elif (dotcell_lattice_start_RE.search(line) and not in_lattice): - outputfile.write("%block LATTICE_CART\n") - outputfile.write(str(defcell[0][0]) + " " + str(defcell[0][1]) + " " + str(defcell[0][2]) + "\n") - outputfile.write(str(defcell[1][0]) + " " + str(defcell[1][1]) + " " + str(defcell[1][2]) + "\n") - outputfile.write(str(defcell[2][0]) + " " + str(defcell[2][1]) + " " + str(defcell[2][2]) + "\n") - outputfile.write("%endblock LATTICE_CART\n") - outputfile.write("FIX_ALL_CELL true\n") - in_lattice = True - elif (dotcell_atoms_end_RE.search(line) and in_atoms and have_atoms): - in_atoms = False - elif ((dotcell_atoms_start_RE.search(line)) and (not in_atoms) and have_atoms): - outputfile.write("%block POSITIONS_FRAC\n") - for atom in atoms: - outputfile.write(" " + atom[0] + " " + str(atom[1]) + " " + str(atom[2]) + " " + str(atom[3]) + "\n") - outputfile.write("%endblock POSITIONS_FRAC\n") - in_atoms = True - elif(not (in_lattice or in_atoms)): - outputfile.write(line) - inputfile.close - outputfile.close - return() + """ + produce_dotcell: reads .cell (CASTEP cell file + and writes a new .cell file to replacing the + lattice block with a new crystalographic lattice + (which should be supplied as a list of three lists, each with + three elements). Also adds command to fix cell during optimization. + """ + in_lattice = False + in_atoms = False + have_atoms = (atoms != []) # If we have an empty list, no atoms were optimized so just leave them in the .cell file. + with open(seedname+".cell", "r") as inputfile, open(filename, "w") as outputfile: + for line in inputfile: + logger.debug("Process line '%s'", line.strip()) + if (dotcell_lattice_end_RE.search(line) and in_lattice): + logger.debug("Reached end of LATTICE_CART block.") + in_lattice = False + elif (dotcell_lattice_start_RE.search(line) and not in_lattice): + logger.debug("Write LATTICE_CART block.") + outputfile.write("%block LATTICE_CART\n") + outputfile.write(str(defcell[0][0]) + " " + str(defcell[0][1]) + " " + str(defcell[0][2]) + "\n") + outputfile.write(str(defcell[1][0]) + " " + str(defcell[1][1]) + " " + str(defcell[1][2]) + "\n") + outputfile.write(str(defcell[2][0]) + " " + str(defcell[2][1]) + " " + str(defcell[2][2]) + "\n") + outputfile.write("%endblock LATTICE_CART\n") + outputfile.write("FIX_ALL_CELL true\n") + in_lattice = True + elif (dotcell_atoms_end_RE.search(line) and in_atoms and have_atoms): + logger.debug("Reached end of POSITIONS_FRAC block.") + in_atoms = False + elif ((dotcell_atoms_start_RE.search(line)) and (not in_atoms) and have_atoms): + logger.debug("Write POSITIONS_FRAC block.") + outputfile.write("%block POSITIONS_FRAC\n") + for atom in atoms: + logger.debug("Write atom '%s'.", atom) + outputfile.write(" " + atom[0] + " " + str(atom[1]) + " " + str(atom[2]) + " " + str(atom[3]) + "\n") + outputfile.write("%endblock POSITIONS_FRAC\n") + in_atoms = True + elif(not (in_lattice or in_atoms)): + logger.debug("Write line back to output file unchanged.") + outputfile.write(line) + else: + logger.debug("Do not write anything to output file.") + return() # regular expression which matches the whole stress tensor block from a .castep file stressRE = re.compile(r"\s\*+\s(?:Symmetrised\s)?Stress\sTensor\s\*+\n.+\n.+?\((\w+)\).+\n.+\n.+\n.+\n\s\*\s+x\s+([\+\-]?\d+.\d+)\s+([\+\-]?\d+.\d+)\s+([\+\-]?\d+.\d+)\s+\*\n\s\*\s+y\s+[\+\-]?\d+.\d+\s+([\+\-]?\d+.\d+)\s+([\+\-]?\d+.\d+)\s+\*\n\s\*\s+z\s+[\+\-]?\d+.\d+\s+[\+\-]?\d+.\d+\s+([\+\-]?\d+.\d+)\s+\*\n") def get_stress_dotcastep(filename): """Extract the stress tensor from a .castep file - + set_facecolor Returns a tuple of (, ) where is a string representing the stress units and is a numpy vector of the elements of the @@ -124,7 +140,7 @@ def get_stress_dotcastep(filename): stressData = stressRE.findall(dotCastep.read())[0] dotCastep.close() units = stressData[0] - stress = S.array([float(stressData[1]),float(stressData[4]), + stress = np.array([float(stressData[1]),float(stressData[4]), float(stressData[6]),float(stressData[5]), float(stressData[3]),float(stressData[2])]) return(units, stress) diff --git a/elastics.py b/elastics.py index 63f548f..de3ec3d 100755 --- a/elastics.py +++ b/elastics.py @@ -8,17 +8,28 @@ """ from __future__ import print_function +import logging import sys import os import re import optparse -import scipy as S +import numpy as np import CijUtil import castep +logger = logging.getLogger(__name__) + version = 1 + +# global P +try: + import matplotlib.pyplot as plt +except ImportError as e: + logger.warning(e) + logger.warning("You need to have matplotlib installed for the --graphics option") + def main(input_options, libmode=False): @@ -28,9 +39,9 @@ def analysePatterns(strain): # these are the IRE conventions, except that integers are 0->5 rather than 1->6 strainDict = {0:"xx",1:"yy",2:"zz",3:"yz", 4:"zx", 5:"xy"} - strainsUsed = S.zeros((6,1)) + strainsUsed = np.zeros((6,1)) - for a in range(0,S.size(strain)): + for a in range(0, np.size(strain)): if strain[a] != 0.0: print(strainDict[a], "component is non-zero") strainsUsed[a] = 1 @@ -43,7 +54,7 @@ def analysePatterns(strain): def cMatrix(symmetryType,TetrHigh): if symmetryType == "Cubic" : - return S.matrix([[1, 7, 7, 0, 0, 0], + return np.matrix([[1, 7, 7, 0, 0, 0], [7, 1, 7, 0, 0, 0], [7, 7, 1, 0, 0, 0], [0, 0, 0, 4, 0, 0], @@ -51,7 +62,7 @@ def cMatrix(symmetryType,TetrHigh): [0, 0, 0, 0, 0, 4]]) elif symmetryType == "Trigonal-high/Hexagonal": - return S.matrix([[1, 7, 8, 9, 0, 0], + return np.matrix([[1, 7, 8, 9, 0, 0], [7, 1, 8, -9, 0, 0], [8, 8, 3, 0, 0, 0], [9, -9, 0, 4, 0, 0], @@ -59,7 +70,7 @@ def cMatrix(symmetryType,TetrHigh): [0, 0, 0, 0, 9, 6]]) elif symmetryType == "Trigonal-low": - return S.matrix([[ 1, 7, 8, 9, 10, 0], + return np.matrix([[ 1, 7, 8, 9, 10, 0], [ 7, 1, 8, -9, -10, 0], [ 8, 8, 3, 0, 0, 0], [ 9, -9, 0, 4, 0, -10], @@ -69,7 +80,7 @@ def cMatrix(symmetryType,TetrHigh): elif symmetryType == "Tetragonal": if TetrHigh == "-1": print("Higher-symmetry tetragonal (422,4mm,4-2m,4/mmm)") - return S.matrix([[1, 7, 8, 0, 0, 0], + return np.matrix([[1, 7, 8, 0, 0, 0], [7, 1, 8, 0, 0, 0], [8, 8, 3, 0, 0, 0], [0, 0, 0, 4, 0, 0], @@ -77,7 +88,7 @@ def cMatrix(symmetryType,TetrHigh): [0, 0, 0, 0, 0, 6]]) else: print("Lower-symmetry tetragonal (4,-4,4/m)") - return S.matrix([[1, 7, 8, 0, 0, 11], + return np.matrix([[1, 7, 8, 0, 0, 11], [7, 1, 8, 0, 0, -11], [8, 8, 3, 0, 0, 0], [0, 0, 0, 4, 0, 0], @@ -85,7 +96,7 @@ def cMatrix(symmetryType,TetrHigh): [11, -11, 0, 0, 0, 6]]) elif symmetryType == "Orthorhombic": - return S.matrix([[ 1, 7, 8, 0, 0, 0], + return np.matrix([[ 1, 7, 8, 0, 0, 0], [ 7, 2, 12, 0, 0, 0], [ 8, 12, 3, 0, 0, 0], [ 0, 0, 0, 4, 0, 0], @@ -93,7 +104,7 @@ def cMatrix(symmetryType,TetrHigh): [ 0, 0, 0, 0, 0, 6]]) elif symmetryType == "Monoclinic": - return S.matrix([[ 1, 7, 8, 0, 10, 0], + return np.matrix([[ 1, 7, 8, 0, 10, 0], [ 7, 2, 12, 0, 14, 0], [ 8, 12, 3, 0, 17, 0], [ 0, 0, 0, 4, 0, 20], @@ -101,7 +112,7 @@ def cMatrix(symmetryType,TetrHigh): [ 0, 0, 0, 20, 0, 6]]) elif symmetryType == "Triclinic": - return S.matrix([[ 1, 7, 8, 9, 10, 11], + return np.matrix([[ 1, 7, 8, 9, 10, 11], [ 7, 2, 12, 13, 14, 15], [ 8, 12, 3, 16, 17, 18], [ 9, 13, 16, 4, 19, 20], @@ -135,14 +146,6 @@ class PotM_Options: global outfile outfile = input_options[0] - if options.graphics: - try: - global P - import pylab as P - except ImportError: - print("You need to have matplotlib installed for the --graphics option", file=sys.stderr) - sys.exit(1) - return options, arguments options, arguments = get_options() @@ -176,28 +179,28 @@ class PotM_Options: # if using graphics, do some initial set-up if options.graphics: - fig = P.figure(num=1, figsize=(9.5,8),facecolor='white') + fig = plt.figure(num=1, figsize=(9.5,8),facecolor='white') fig.subplots_adjust(left=0.07,right=0.97,top=0.97,bottom=0.07,wspace=0.5,hspace=0.5) colourDict = {0: '#BAD0EF', 1:'#FFCECE', 2:'#BDF4CB', 3:'#EEF093',4:'#FFA4FF',5:'#75ECFD'} for index1 in range(6): for index2 in range(6): # position this plot in a 6x6 grid - sp = P.subplot(6,6,6*(index1)+index2+1) + sp = plt.subplot(6,6,6*(index1)+index2+1) sp.set_axis_off() # change the labels on the axes # xlabels = sp.get_xticklabels() - # P.setp(xlabels,'rotation',90,fontsize=7) + # plt.setp(xlabels,'rotation',90,fontsize=7) # ylabels = sp.get_yticklabels() - # P.setp(ylabels,fontsize=7) - P.text(0.4,0.4, "n/a") + # plt.setp(ylabels,fontsize=7) + plt.text(0.4,0.4, "n/a") print("\n<>---------------------------- ANALYSIS ---------------------------------<>") # initialise 1d array to store all 21 unique elastic constants - will be transformed into 6x6 matrix later - finalCijs = S.zeros((21,1)) - errors = S.zeros((21,1)) - + finalCijs = np.zeros((21,1)) + errors = np.zeros((21,1)) + for patt in range(numStrainPatterns//numsteps): print("\nAnalysing pattern", patt+1, ":") @@ -215,9 +218,9 @@ class PotM_Options: # numbering according to IRE conventions (Proc IRE, 1949) if a == 0: - strain = S.array([float(line1[0]),float(line2[1]),float(line3[2]),2*float(line2[2]),2*float(line1[2]),2*float(line1[1])]) + strain = np.array([float(line1[0]),float(line2[1]),float(line3[2]),2*float(line2[2]),2*float(line1[2]),2*float(line1[1])]) else: - strain = S.row_stack((strain,S.array([float(line1[0]),float(line2[1]),float(line3[2]),2*float(line2[2]),2*float(line1[2]),2*float(line1[1])]))) + strain = np.row_stack((strain,np.array([float(line1[0]),float(line2[1]),float(line3[2]),2*float(line2[2]),2*float(line1[2]),2*float(line1[1])]))) # now get corresponding stress data from .castep (units, thisStress) = castep.get_stress_dotcastep(seedname+ @@ -227,7 +230,7 @@ class PotM_Options: if a == 0: stress = thisStress else: - stress = S.row_stack((stress,thisStress)) + stress = np.row_stack((stress,thisStress)) """ @@ -239,12 +242,12 @@ class PotM_Options: """ def __fit(index1, index2): - from scipy import stats, sqrt, square + from scipy import stats # do the fit (cijFitted,intercept,r,tt,stderr) = stats.linregress(strain[:,index2-1],stress[:,index1-1]) - (vmajor,vminor,vmicro) = re.split('\.',S.__version__) + (vmajor,vminor,vmicro) = re.split('\.',np.__version__) if ( int(vmajor) > 0 or int(vminor) >= 7): error = stderr @@ -252,8 +255,8 @@ def __fit(index1, index2): # correct for scipy weirdness - see http://www.scipy.org/scipy/scipy/ticket/8 # This was fixed before 0.7.0 release. Maybe in some versions of 0.6.x too - # will report huge errors if the check is wrong - stderr = S.sqrt((numsteps * stderr**2)/(numsteps-2)) - error = stderr/sqrt(sum(square(strain[:,index2-1]))) + stderr = np.sqrt((numsteps * stderr**2)/(numsteps-2)) + error = stderr/np.sqrt(np.sum(np.square(strain[:,index2-1]))) # print info about the fit print('\n') @@ -269,21 +272,21 @@ def __fit(index1, index2): if options.graphics: # position this plot in a 6x6 grid - sp = P.subplot(6,6,6*(index1-1)+index2) + sp = plt.subplot(6,6,6*(index1-1)+index2) sp.set_axis_on() # change the labels on the axes xlabels = sp.get_xticklabels() - P.setp(xlabels,'rotation',90,fontsize=7) + plt.setp(xlabels,'rotation',90,fontsize=7) ylabels = sp.get_yticklabels() - P.setp(ylabels,fontsize=7) + plt.setp(ylabels,fontsize=7) # colour the plot depending on the strain pattern - sp.set_axis_bgcolor(colourDict[patt]) + sp.set_facecolor(colourDict[patt]) # plot the data - P.plot([strain[0,index2-1],strain[numsteps-1,index2-1]],[cijFitted*strain[0,index2-1]+intercept,cijFitted*strain[numsteps-1,index2-1]+intercept]) - P.plot(strain[:,index2-1],stress[:,index1-1],'ro') + plt.plot([strain[0,index2-1],strain[numsteps-1,index2-1]],[cijFitted*strain[0,index2-1]+intercept,cijFitted*strain[numsteps-1,index2-1]+intercept]) + plt.plot(strain[:,index2-1],stress[:,index1-1],'ro') return cijFitted, error @@ -292,7 +295,7 @@ def __appendOrReplace(valList,erList,val): try: valList.append(val[0]) erList.append(val[1]) - return (sum(valList)/len(valList)), (S.sqrt(sum([x**2 for x in erList])/len(erList)**2)) + return (sum(valList)/len(valList)), (np.sqrt(sum([x**2 for x in erList])/len(erList)**2)) except NameError: return val[0], val[1] @@ -305,7 +308,7 @@ def __createListAndAppend(val): return val[0], newList, val[1], errorList - cij = S.zeros(21) + cij = np.zeros(21) # Analyse the patterns to see which strains were applied strainsUsed = analysePatterns(strain[0,:]) @@ -314,13 +317,13 @@ def __createListAndAppend(val): if symmetryType == "Cubic": - if S.all(strainsUsed.transpose() == S.array([[1.0, 0.0, 0.0, 1.0, 0.0, 0.0]])): # strain pattern e1+e4 + if np.all(strainsUsed.transpose() == np.array([[1.0, 0.0, 0.0, 1.0, 0.0, 0.0]])): # strain pattern e1+e4 finalCijs[0], errors[0] = __fit(1,1) # fit C11 fit_21, fit_21_error = __fit(2,1) fit_31, fit_31_error = __fit(3,1) finalCijs[6] = (fit_21 + fit_31)/2 # fit C21+C31 - errors[6] = S.sqrt((fit_21_error**2)/4 + (fit_31_error**2)/4) + errors[6] = np.sqrt((fit_21_error**2)/4 + (fit_31_error**2)/4) finalCijs[3], errors[3] = __fit(4,4) # fit C44 else: @@ -328,7 +331,7 @@ def __createListAndAppend(val): sys.exit(1) elif symmetryType == "Trigonal-high/Hexagonal": - if S.all(strainsUsed.transpose() == S.array([[0.0, 0.0, 1.0, 0.0, 0.0, 0.0]])): # strain pattern e3 (hexagonal) + if np.all(strainsUsed.transpose() == np.array([[0.0, 0.0, 1.0, 0.0, 0.0, 0.0]])): # strain pattern e3 (hexagonal) # fit C13 + C23, and add to list (more values coming...) finalCijs[7], cij13, errors[7], er13 = __createListAndAppend(__fit(1,3)) @@ -336,14 +339,14 @@ def __createListAndAppend(val): finalCijs[2], errors[2] = __fit(3,3) # fit C33 - elif S.all(strainsUsed.transpose() == S.array([[1.0, 0.0, 0.0, 1.0, 0.0, 0.0]])): # strain pattern e1+e4 (hexagonal) + elif np.all(strainsUsed.transpose() == np.array([[1.0, 0.0, 0.0, 1.0, 0.0, 0.0]])): # strain pattern e1+e4 (hexagonal) finalCijs[0], errors[0] = __fit(1,1) # fit C11 finalCijs[6], errors[6] = __fit(2,1) # fit C21 finalCijs[7], errors[7] = __appendOrReplace(cij13,er13,__fit(3,1)) # fit C31 finalCijs[3], errors[3] = __fit(4,4) # fit C44 - elif S.all(strainsUsed.transpose() == S.array([[1.0, 0.0, 0.0, 0.0, 0.0, 0.0]])): + elif np.all(strainsUsed.transpose() == np.array([[1.0, 0.0, 0.0, 0.0, 0.0, 0.0]])): # strain pattern e1 (trigonal-high) @@ -354,7 +357,7 @@ def __createListAndAppend(val): # Should be zero? finalCijs[9], errors[9] = __fit(5,1) # fit C51 - elif S.all(strainsUsed.transpose() == S.array([[0.0, 0.0, 1.0, 1.0, 0.0, 0.0]])): + elif np.all(strainsUsed.transpose() == np.array([[0.0, 0.0, 1.0, 1.0, 0.0, 0.0]])): # strain pattern e3+e4 (trigonal-high) # could recalculate C13/C14/C23/C24/C46 here, but won't just now @@ -367,7 +370,7 @@ def __createListAndAppend(val): sys.exit(1) elif symmetryType == "Trigonal-low": - if S.all(strainsUsed.transpose() == S.array([[1.0, 0.0, 0.0, 0.0, 0.0, 0.0]])): + if np.all(strainsUsed.transpose() == np.array([[1.0, 0.0, 0.0, 0.0, 0.0, 0.0]])): # strain pattern e1 @@ -377,7 +380,7 @@ def __createListAndAppend(val): finalCijs[8], errors[8] = __fit(4,1) # fit C41 finalCijs[9], errors[9] = __fit(5,1) # fit C51 - elif S.all(strainsUsed.transpose() == S.array([[0.0, 0.0, 1.0, 1.0, 0.0, 0.0]])): + elif np.all(strainsUsed.transpose() == np.array([[0.0, 0.0, 1.0, 1.0, 0.0, 0.0]])): # strain pattern e3+e4 # could recalculate C13/C14/C23/C24/C46 here, but won't just now @@ -390,7 +393,7 @@ def __createListAndAppend(val): sys.exit(1) elif symmetryType == "Tetragonal": - if S.all(strainsUsed.transpose() == S.array([[1.0, 0.0, 0.0, 1.0, 0.0, 0.0]])): # strain pattern e1+e4 + if np.all(strainsUsed.transpose() == np.array([[1.0, 0.0, 0.0, 1.0, 0.0, 0.0]])): # strain pattern e1+e4 finalCijs[0], errors[0] = __fit(1,1) # fit C11 finalCijs[6], errors[6] = __fit(2,1) # fit C21 @@ -399,7 +402,7 @@ def __createListAndAppend(val): finalCijs[3], errors[3] = __fit(4,4) # fit C44 - elif S.all(strainsUsed.transpose() == S.array([[0.0, 0.0, 1.0, 0.0, 0.0, 1.0]])): # strain pattern e3+e6 + elif np.all(strainsUsed.transpose() == np.array([[0.0, 0.0, 1.0, 0.0, 0.0, 1.0]])): # strain pattern e3+e6 finalCijs[2], errors[2] = __fit(3,3) # fit C33 finalCijs[5], errors[5] = __fit(6,6) # fit C66 @@ -409,14 +412,14 @@ def __createListAndAppend(val): sys.exit(1) elif symmetryType == "Orthorhombic": - if S.all(strainsUsed.transpose() == S.array([[1.0, 0.0, 0.0, 1.0, 0.0, 0.0]])): # strain pattern e1+e4 + if np.all(strainsUsed.transpose() == np.array([[1.0, 0.0, 0.0, 1.0, 0.0, 0.0]])): # strain pattern e1+e4 finalCijs[0], errors[0] = __fit(1,1) # fit C11 finalCijs[6], cij12, errors[6], er12 = __createListAndAppend(__fit(2,1)) # fit C21 finalCijs[7], cij13, errors[7], er13 = __createListAndAppend(__fit(3,1)) # fit C31 finalCijs[3], errors[3] = __fit(4,4) # fit C44 - elif S.all(strainsUsed.transpose() == S.array([[0.0, 1.0, 0.0, 0.0, 1.0, 0.0]])): # strain pattern e2+e5 + elif np.all(strainsUsed.transpose() == np.array([[0.0, 1.0, 0.0, 0.0, 1.0, 0.0]])): # strain pattern e2+e5 finalCijs[6], errors[6] = __appendOrReplace(cij12,er12,__fit(1,2)) # fit C12 @@ -424,7 +427,7 @@ def __createListAndAppend(val): finalCijs[11], cij23, errors[11], er23 = __createListAndAppend(__fit(3,2)) # fit C32 finalCijs[4], errors[4] = __fit(5,5) # fit C55 - elif S.all(strainsUsed.transpose() == S.array([[0.0, 0.0, 1.0, 0.0, 0.0, 1.0]])): # strain pattern e3+e6 + elif np.all(strainsUsed.transpose() == np.array([[0.0, 0.0, 1.0, 0.0, 0.0, 1.0]])): # strain pattern e3+e6 finalCijs[7], errors[7] = __appendOrReplace(cij13,er13,__fit(1,3)) # fit C13 finalCijs[11], errors[11] = __appendOrReplace(cij23,er23,__fit(2,3)) # fit C23 @@ -436,7 +439,7 @@ def __createListAndAppend(val): sys.exit(1) elif symmetryType == "Monoclinic": - if S.all(strainsUsed.transpose() == S.array([[1.0, 0.0, 0.0, 1.0, 0.0, 0.0]])): # strain pattern e1+e4 + if np.all(strainsUsed.transpose() == np.array([[1.0, 0.0, 0.0, 1.0, 0.0, 0.0]])): # strain pattern e1+e4 finalCijs[0], errors[0] = __fit(1,1) # fit C11 finalCijs[6], cij12, errors[6], er12 = __createListAndAppend(__fit(2,1)) # fit C21 @@ -446,7 +449,7 @@ def __createListAndAppend(val): finalCijs[19], cij64, errors[19], er64 = __createListAndAppend(__fit(6,4)) # fit C64 - elif S.all(strainsUsed.transpose() == S.array([[0.0, 0.0, 1.0, 0.0, 0.0, 1.0]])): # strain pattern e3+e6 + elif np.all(strainsUsed.transpose() == np.array([[0.0, 0.0, 1.0, 0.0, 0.0, 1.0]])): # strain pattern e3+e6 finalCijs[7], errors[7] = __appendOrReplace(cij13,er13,__fit(1,3)) # fit C13 finalCijs[11], cij23, errors[11], er23 = __createListAndAppend(__fit(2,3)) # fit C23 @@ -455,7 +458,7 @@ def __createListAndAppend(val): finalCijs[19], errors[19] = __appendOrReplace(cij64,er64,__fit(4,6)) # fit C46 finalCijs[5], errors[5] = __fit(6,6) # fit C66 - elif S.all(strainsUsed.transpose() == S.array([[0.0, 1.0, 0.0, 0.0, 0.0, 0.0]])): # strain pattern e2 + elif np.all(strainsUsed.transpose() == np.array([[0.0, 1.0, 0.0, 0.0, 0.0, 0.0]])): # strain pattern e2 finalCijs[6], errors[6] = __appendOrReplace(cij12,er12,__fit(1,2)) # fit C12 finalCijs[1], errors[1] = __fit(2,2) # fit C22 @@ -463,7 +466,7 @@ def __createListAndAppend(val): finalCijs[13], cij52, errors[13], er52 = __createListAndAppend(__fit(5,2)) # fit C52 - elif S.all(strainsUsed.transpose() == S.array([[0.0, 0.0, 0.0, 0.0, 1.0, 0.0]])): # strain pattern e5 + elif np.all(strainsUsed.transpose() == np.array([[0.0, 0.0, 0.0, 0.0, 1.0, 0.0]])): # strain pattern e5 finalCijs[9], errors[9] = __appendOrReplace(cij51,er51,__fit(1,5)) # fit C15 finalCijs[13],errors[13] = __appendOrReplace(cij52,er52,__fit(2,5)) # fit C25 @@ -475,7 +478,7 @@ def __createListAndAppend(val): elif symmetryType == "Triclinic": - if S.all(strainsUsed.transpose() == S.array([[1.0, 0.0, 0.0, 0.0, 0.0, 0.0]])): # strain pattern e1 + if np.all(strainsUsed.transpose() == np.array([[1.0, 0.0, 0.0, 0.0, 0.0, 0.0]])): # strain pattern e1 finalCijs[0], errors[0] = __fit(1,1) # fit C11 finalCijs[6], cij12, errors[6], er12 = __createListAndAppend(__fit(2,1)) # fit C21 @@ -484,7 +487,7 @@ def __createListAndAppend(val): finalCijs[9], cij15, errors[9], er15 = __createListAndAppend(__fit(5,1)) # fit C51 finalCijs[10],cij16, errors[10],er16 = __createListAndAppend(__fit(6,1)) # fit C61 - elif S.all(strainsUsed.transpose() == S.array([[0.0, 1.0, 0.0, 0.0, 0.0, 0.0]])): # strain pattern e2 + elif np.all(strainsUsed.transpose() == np.array([[0.0, 1.0, 0.0, 0.0, 0.0, 0.0]])): # strain pattern e2 finalCijs[6], errors[6] = __appendOrReplace(cij12,er12,__fit(1,2)) # fit C12 finalCijs[1], errors[1] = __fit(2,2) # fit C22 @@ -493,7 +496,7 @@ def __createListAndAppend(val): finalCijs[13], cij25, errors[13], er25 = __createListAndAppend(__fit(5,2)) # fit C52 finalCijs[14], cij26, errors[14], er26 = __createListAndAppend(__fit(6,2)) # fit C62 - elif S.all(strainsUsed.transpose() == S.array([[0.0, 0.0, 1.0, 0.0, 0.0, 0.0]])): # strain pattern e3 + elif np.all(strainsUsed.transpose() == np.array([[0.0, 0.0, 1.0, 0.0, 0.0, 0.0]])): # strain pattern e3 finalCijs[7], errors[7] = __appendOrReplace(cij13,er13,__fit(1,3)) # fit C13 finalCijs[11], errors[11] = __appendOrReplace(cij23,er23,__fit(2,3)) # fit C23 @@ -502,7 +505,7 @@ def __createListAndAppend(val): finalCijs[16], cij35, errors[16], er35 = __createListAndAppend(__fit(5,3)) # fit C53 finalCijs[17], cij36, errors[17], er36 = __createListAndAppend(__fit(6,3)) # fit C63 - elif S.all(strainsUsed.transpose() == S.array([[0.0, 0.0, 0.0, 1.0, 0.0, 0.0]])): # strain pattern e4 + elif np.all(strainsUsed.transpose() == np.array([[0.0, 0.0, 0.0, 1.0, 0.0, 0.0]])): # strain pattern e4 finalCijs[8], errors[8] = __appendOrReplace(cij14,er14,__fit(1,4)) # fit C14 finalCijs[12], errors[12] = __appendOrReplace(cij24,er24,__fit(2,4)) # fit C24 @@ -511,7 +514,7 @@ def __createListAndAppend(val): finalCijs[18], cij45, errors[18], er45 = __createListAndAppend(__fit(5,4)) # fit C54 finalCijs[19], cij46, errors[19], er46 = __createListAndAppend(__fit(6,4)) # fit C64 - elif S.all(strainsUsed.transpose() == S.array([[0.0, 0.0, 0.0, 0.0, 1.0, 0.0]])): # strain pattern e5 + elif np.all(strainsUsed.transpose() == np.array([[0.0, 0.0, 0.0, 0.0, 1.0, 0.0]])): # strain pattern e5 finalCijs[9], errors[9] = __appendOrReplace(cij15,er15,__fit(1,5)) # fit C15 finalCijs[13], errors[13] = __appendOrReplace(cij25,er25,__fit(2,5)) # fit C25 @@ -520,7 +523,7 @@ def __createListAndAppend(val): finalCijs[4], errors[4] = __fit(5,5) # fit C55 finalCijs[20], cij56, errors[20], er56 = __createListAndAppend(__fit(6,5)) # fit C65 - elif S.all(strainsUsed.transpose() == S.array([[0.0, 0.0, 0.0, 0.0, 0.0, 1.0]])): # strain pattern e6 + elif np.all(strainsUsed.transpose() == np.array([[0.0, 0.0, 0.0, 0.0, 0.0, 1.0]])): # strain pattern e6 finalCijs[10], errors[10] = __appendOrReplace(cij16,er16,__fit(1,6)) # fit C16 finalCijs[14], errors[14] = __appendOrReplace(cij26,er26,__fit(2,6)) # fit C26 @@ -537,21 +540,21 @@ def __createListAndAppend(val): sys.exit(1) if options.graphics: - P.savefig(os.path.basename(seedname)+'_fits') + plt.savefig(os.path.basename(seedname)+'_fits') cijdat.close() if symmetryType == "Trigonal-high/Hexagonal" or symmetryType == "Trigonal-low": # for these systems, C66 is calculated as a combination of the other Cijs. finalCijs[5] = 0.5*(finalCijs[0]-finalCijs[6]) - errors[5] = S.sqrt(0.25*(errors[0]**2+errors[6]**2)) + errors[5] = np.sqrt(0.25*(errors[0]**2+errors[6]**2)) c = cMatrix(symmetryType,TetrHigh) # Generate the 6x6 matrix of elastic constants # - negative values signify a symmetry relation - finalCijMatrix = S.zeros((6,6)) - finalErrors = S.zeros((6,6)) + finalCijMatrix = np.zeros((6,6)) + finalErrors = np.zeros((6,6)) for i in range(0,6): for j in range(0,6): index = int(c[i,j]) @@ -574,16 +577,16 @@ def __createListAndAppend(val): print("\n<>---------------------------- RESULTS ----------------------------------<>\n") print("Final Cij matrix ("+units+"):") - print(S.array2string(finalCijMatrix,max_line_width=130,suppress_small=True)) + print(np.array2string(finalCijMatrix,max_line_width=130,suppress_small=True)) print("\nErrors on Cij matrix ("+units+"):") - print(S.array2string(finalErrors,max_line_width=130,suppress_small=True)) + print(np.array2string(finalErrors,max_line_width=130,suppress_small=True)) (sij, esij, covsij) = CijUtil.invertCij(finalCijMatrix,finalErrors) print("\nFinal Sij matrix ("+units+"-1):") - print(S.array2string(sij,max_line_width=130,suppress_small=True)) + print(np.array2string(sij,max_line_width=130,suppress_small=True)) print("\nErrors on Sij matrix ("+units+"-1):") - print(S.array2string(esij,max_line_width=130,suppress_small=True)) + print(np.array2string(esij,max_line_width=130,suppress_small=True)) print("\n<>----------------------------------------------------------------------<>\n") if symmetryType == "Cubic": @@ -615,7 +618,7 @@ def __createListAndAppend(val): print("\n<>-----------------------------------------------------------------------<>\n") - S.savetxt(seedname + '_cij.txt', finalCijMatrix) + np.savetxt(seedname + '_cij.txt', finalCijMatrix) if options.latex: CijUtil.latexCij(finalCijMatrix, finalErrors, seedname + '.tex', options.latex_nt) if options.txt: diff --git a/generate_strain.py b/generate_strain.py index 3725106..ec9db3f 100755 --- a/generate_strain.py +++ b/generate_strain.py @@ -11,6 +11,7 @@ """ from __future__ import print_function +import logging import sys import os import optparse @@ -22,6 +23,7 @@ version = 0.1 +logger = logging.getLogger(__name__) def PointGroup2StrainPat(pointgroup): """ @@ -131,6 +133,8 @@ def get_options(input_options, libmode): version="%prog "+str(version)) p.add_option('--debug', '-d', action='store_true', help="Debug mode (output to stdout rather than file)") + p.add_option('--verbose', '-v', action='store_true', + help="Verbose mode (extensive logging)") p.add_option('--steps', '-n', action='store', type='int', dest="numsteps", help='Number of positive strain magnitudes to impose' + @@ -213,6 +217,21 @@ def main(input_options, libmode=False): options, arguments = get_options(input_options, libmode) seedname = arguments[0] + logformat = "%(levelname)s: %(message)s" + if options.debug: + # detailed log format for debug output + logformat = ( + "[%(asctime)s-%(funcName)s()-%(filename)s:%(lineno)s]" + " %(levelname)s: %(message)s" + ) + loglevel = logging.DEBUG + elif options.verbose: + loglevel = logging.INFO + else: + loglevel = logging.WARNING + + logging.basicConfig(level=loglevel, format=logformat) + (cell,pointgroup,atoms) = castep.parse_dotcastep(seedname) # Re-align lattice vectors on cartisian system a, b, c, al, be, ga = cellCART2cellABC(cell) @@ -262,51 +281,52 @@ def main(input_options, libmode=False): print("Lattice type is ", latticeTypes[latticeCode]) print("Number of patterns: "+ str(numStrainPatterns) +"\n") - cijdat = open(seedname+".cijdat","w") - print("Writing strain data to ", seedname+".cijdat") - cijdat.write(str(latticeCode) + ' ' + str(numsteps*2) + ' 0 ' + '0 \n') - cijdat.write(str(maxstrain)+"\n") - - # The order of these three loops matters for the analysis code. - for patt in range(numStrainPatterns): - this_pat = patterns[patt] - for a in range(0,numsteps): - for neg in range(0,2): - if (neg == 0): - this_mag = (float(a)+1) / (float(numsteps)) * maxstrain - else: - this_mag = -1.0 * (float(a)+1) / (float(numsteps)) * maxstrain - disps = [x * this_mag for x in patterns[patt]] - # Build the strain tensor (IRE convention but 1 -> 0 etc.) - this_strain = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0] - # diagonal elements - strain is displacment / lattice vector length - this_strain[0] = disps[0] / np.sqrt(cell[0][0]**2+cell[0][1]**2+cell[0][2]**2) - this_strain[1] = disps[1] / np.sqrt(cell[1][0]**2+cell[1][1]**2+cell[1][2]**2) - this_strain[2] = disps[2] / np.sqrt(cell[2][0]**2+cell[2][1]**2+cell[2][2]**2) - # off diagonals - we only strain upper right corner of cell matrix, so strain is 1/2*du/dx... - this_strain[3] = 0.5 * (disps[3] / np.sqrt(cell[1][0]**2+cell[1][1]**2+cell[1][2]**2)) - this_strain[4] = 0.5 * (disps[4] / np.sqrt(cell[0][0]**2+cell[0][1]**2+cell[0][2]**2)) - this_strain[5] = 0.5 * (disps[5] / np.sqrt(cell[0][0]**2+cell[0][1]**2+cell[0][2]**2)) - - # Deform cell - only apply deformation to upper right corner - defcell = [[cell[0][0]+disps[0], cell[0][1]+disps[5], cell[0][2]+disps[4]], - [cell[1][0], cell[1][1]+disps[1], cell[1][2]+disps[3]], - [cell[2][0], cell[2][1], cell[2][2]+disps[2]]] - - pattern_name = seedname + "_cij__" + str(patt+1) + "__" + str((a*2)+1+neg) - - print("Pattern Name = ", pattern_name) - print("Pattern = ", this_pat) - print("Magnitude = ", this_mag) - cijdat.write(pattern_name+"\n") - cijdat.write(str(this_strain[0]) + " " + str(this_strain[5]) + " " + str(this_strain[4]) + "\n") - cijdat.write(str(this_strain[5]) + " " + str(this_strain[1]) + " " + str(this_strain[3]) + "\n") - cijdat.write(str(this_strain[4]) + " " + str(this_strain[3]) + " " + str(this_strain[2]) + "\n") - castep.produce_dotcell(seedname, pattern_name+".cell", defcell, atoms) - os.symlink(seedname+".param", pattern_name+".param") - - - + with open(seedname+".cijdat","w") as cijdat: + print("Writing strain data to ", seedname+".cijdat") + cijdat.write(str(latticeCode) + ' ' + str(numsteps*2) + ' 0 ' + '0 \n') + cijdat.write(str(maxstrain)+"\n") + + # The order of these three loops matters for the analysis code. + for patt in range(numStrainPatterns): + this_pat = patterns[patt] + for a in range(0,numsteps): + for neg in range(0,2): + if (neg == 0): + this_mag = (float(a)+1) / (float(numsteps)) * maxstrain + else: + this_mag = -1.0 * (float(a)+1) / (float(numsteps)) * maxstrain + disps = [x * this_mag for x in patterns[patt]] + # Build the strain tensor (IRE convention but 1 -> 0 etc.) + this_strain = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0] + # diagonal elements - strain is displacment / lattice vector length + this_strain[0] = disps[0] / np.sqrt(cell[0][0]**2+cell[0][1]**2+cell[0][2]**2) + this_strain[1] = disps[1] / np.sqrt(cell[1][0]**2+cell[1][1]**2+cell[1][2]**2) + this_strain[2] = disps[2] / np.sqrt(cell[2][0]**2+cell[2][1]**2+cell[2][2]**2) + # off diagonals - we only strain upper right corner of cell matrix, so strain is 1/2*du/dx... + this_strain[3] = 0.5 * (disps[3] / np.sqrt(cell[1][0]**2+cell[1][1]**2+cell[1][2]**2)) + this_strain[4] = 0.5 * (disps[4] / np.sqrt(cell[0][0]**2+cell[0][1]**2+cell[0][2]**2)) + this_strain[5] = 0.5 * (disps[5] / np.sqrt(cell[0][0]**2+cell[0][1]**2+cell[0][2]**2)) + + # Deform cell - only apply deformation to upper right corner + defcell = [[cell[0][0]+disps[0], cell[0][1]+disps[5], cell[0][2]+disps[4]], + [cell[1][0], cell[1][1]+disps[1], cell[1][2]+disps[3]], + [cell[2][0], cell[2][1], cell[2][2]+disps[2]]] + + pattern_name = seedname + "_cij__" + str(patt+1) + "__" + str((a*2)+1+neg) + + print("Pattern Name = ", pattern_name) + print("Pattern = ", this_pat) + print("Magnitude = ", this_mag) + cijdat.write(pattern_name+"\n") + cijdat.write(str(this_strain[0]) + " " + str(this_strain[5]) + " " + str(this_strain[4]) + "\n") + cijdat.write(str(this_strain[5]) + " " + str(this_strain[1]) + " " + str(this_strain[3]) + "\n") + cijdat.write(str(this_strain[4]) + " " + str(this_strain[3]) + " " + str(this_strain[2]) + "\n") + castep.produce_dotcell(seedname, pattern_name+".cell", defcell, atoms) + try: + os.symlink(seedname+".param", pattern_name+".param") + except FileExistsError as e: + logger.warning(e) + if __name__ == "__main__": main(sys.argv[1:])