forked from lijiext/lammps
95 lines
2.0 KiB
Python
Executable File
95 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 __future__ import print_function
|
|
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
|
|
|
|
with open("log.lammps",'r') as logfile:
|
|
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], end="")
|
|
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], end="")
|
|
print()
|