forked from sr76/excitingscripts
-
Notifications
You must be signed in to change notification settings - Fork 0
/
THERMAL-cubic.py
executable file
·100 lines (77 loc) · 3.56 KB
/
THERMAL-cubic.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
#!/usr/bin/python
# -*- coding: utf-8 -*-
#_______________________________________________________________________________
import sys
import os
from lxml import etree
import numpy
from pylab import *
from scipy import interpolate
#-------------------------------------------------------------------------------
def checkfile(file):
if (str(os.path.exists(file))=='False'):
sys.exit("ERROR: File "+file+" not found!\n")
return open(file,"r")
#-------------------------------------------------------------------------------
def readfvib(f,t,a):
return t.xpath('/thermodynamicproperties/'+a+'/map/@'+f)
#-------------------------------------------------------------------------------
narg = len(sys.argv)-1
if (narg<3):
print "\nIncorrect number of arguments. **Usage**:\n\n",
print "THERMAL-cubic.py TMIN TMAX NTSTEPS\n"
print "Temperatures should be given in Kelvin\n"
sys.exit()
tmin = float(sys.argv[1]) ; tmax = float(sys.argv[2]) ; ntpt = int(sys.argv[3])
tstep = (tmax-tmin)/float(ntpt) ; temp = []
for i in range(ntpt+1): temp.append(tmin+i*tstep)
#-------------------------------------------------------------------------------
print
print '------------------------------------------------------------------------'
print 'Linear thermal-expansion coefficient for cubic systems'
print '------------------------------------------------------------------------'
bzero = input("\nEnter value for the bulk modulus [GPa] >>>> ")
print
#-------------------------------------------------------------------------------
iname = [] ; tname = []
iname.append("input.xml.-") ; tname.append("thermo.xml.-")
iname.append("input.xml.0") ; tname.append("thermo.xml.0")
iname.append("input.xml.+") ; tname.append("thermo.xml.+")
ifile = [] ; tfile = [] ; itree = [] ; ttree = [] ; alat = [] ; vol = []
tttt = [] ; ffff = [] ; ssss = [] ; cccc = [] ; vvvv = [] ; fvib = []
dvib = []
for i in range(len(iname)):
ifile.append(checkfile(iname[i]))
tfile.append(checkfile(iname[i]))
itree.append(etree.parse(iname[i]))
ttree.append(etree.parse(tname[i]))
alat.append(map(float,itree[i].xpath('/input/structure/crystal/@scale'))[0])
vol.append(alat[i]**3/4.)
tttt.append(readfvib("variable1",ttree[i],'vibrationalfreeenergy'))
ffff.append(readfvib("function1",ttree[i],'vibrationalfreeenergy'))
ssss.append(readfvib("function1",ttree[i],'vibrationalentropy'))
cccc.append(readfvib("function1",ttree[i],'heatcapacity'))
vvvv.append(readfvib("function1",ttree[i],'vibrationalenergy'))
tttt[i].insert(0,"0.") ; ffff[i].insert(0,ffff[i][0])
splf = interpolate.splrep(tttt[i],ffff[i],s=0)
fvib.append(interpolate.splev(temp,splf,der=0))
dvib.append(interpolate.splev(temp,splf,der=1))
outf = open('linear-thermal-expansion','w')
GPa2au = 3.398827e-5
bzero = bzero*GPa2au
for i in range(len(tttt[0])): tttt[0][i]=float(tttt[0][i])
for i in range(len(tttt[1])): tttt[1][i]=float(tttt[1][i])
for i in range(len(tttt[2])): tttt[2][i]=float(tttt[2][i])
tcut= min(max(tttt[0]),max(tttt[1]),max(tttt[2]))
for i in range(len(temp)):
alpha=(dvib[0][i]-dvib[2][i])/(vol[2]-vol[0])/bzero/3/2
if (temp[i] < tcut+1e-10): print>>outf, temp[i], alpha
outf.close()
#-------------------------------------------------------------------------------
style = " o-"
if (ntpt > 40): style = "-"
xlabel = " \"Temperature [K]\""
ylabel = " \"alpha(T) [1/K]\""
options = style+xlabel+ylabel
os.system("PLOT-plot.py linear-thermal-expansion "+options)
#-------------------------------------------------------------------------------