Added a script to calculate elastic compliance tensor

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14819 f3b2605a-c512-4ea7-a41b-209d697bcdaa

examples/ELASTIC/compliance.py

@@ -0,0 +1,95 @@
+#!/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
+