#!/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