lammps/examples/ELASTIC/compliance.py

96 lines
2.0 KiB
Python
Executable File

#!/usr/bin/env python
# This file reads in the file log.lammps generated by the script ELASTIC/in.elastic
# It prints out the 6x6 tensor of elastic constants Cij
# followed by the 6x6 tensor of compliance constants Sij
# It uses the same conventions as described in:
# Sprik, Impey and Klein PRB (1984).
# The units of Cij are whatever was used in log.lammps (usually GPa)
# The units of Sij are the inverse of that (usually 1/GPa)
from numpy import zeros
from numpy.linalg import inv
# define logfile layout
nvals = 21
valpos = 4
valstr = '\nElastic Constant C'
# define order of Cij in logfile
cindices = [0]*nvals
cindices[0] = (0,0)
cindices[1] = (1,1)
cindices[2] = (2,2)
cindices[3] = (0,1)
cindices[4] = (0,2)
cindices[5] = (1,2)
cindices[6] = (3,3)
cindices[7] = (4,4)
cindices[8] = (5,5)
cindices[9] = (0,3)
cindices[10] = (0,4)
cindices[11] = (0,5)
cindices[12] = (1,3)
cindices[13] = (1,4)
cindices[14] = (1,5)
cindices[15] = (2,3)
cindices[16] = (2,4)
cindices[17] = (2,5)
cindices[18] = (3,4)
cindices[19] = (3,5)
cindices[20] = (4,5)
# open logfile
logfile = open("log.lammps",'r')
txt = logfile.read()
# search for 21 elastic constants
c = zeros((6,6))
s2 = 0
for ival in range(nvals):
s1 = txt.find(valstr,s2)
if (s1 == -1):
print "Failed to find elastic constants in log file"
exit(1)
s1 += 1
s2 = txt.find("\n",s1)
line = txt[s1:s2]
# print line
words = line.split()
(i1,i2) = cindices[ival]
c[i1,i2] = float(words[valpos])
c[i2,i1] = c[i1,i2]
print "C tensor [GPa]"
for i in range(6):
for j in range(6):
print "%10.8g " % c[i][j],
print
# apply factor of 2 to columns of off-diagonal elements
for i in range(6):
for j in range(3,6):
c[i][j] *= 2.0
s = inv(c)
# apply factor of 1/2 to columns of off-diagonal elements
for i in range(6):
for j in range(3,6):
s[i][j] *= 0.5
print "S tensor [1/GPa]"
for i in range(6):
for j in range(6):
print "%10.8g " % s[i][j],
print