+- Python Forum (https://python-forum.io)
+-- Forum: Python Coding (https://python-forum.io/forum-7.html)
+--- Forum: General Coding Help (https://python-forum.io/forum-8.html)
+--- Thread: modify script. (/thread-22598.html)
Pages:
1
2
modify script. - MuhammadTalha - Nov-19-2019
Hi i want to modify python script, so it can not look for a specific file.
In my case i want it to not look for .dvscf1.
Regards
Muhammad Talha
RE: modify script. - Gribouillis - Nov-19-2019
Without the script, it is impossible to say how to modify it.
RE: modify script. - MuhammadTalha - Nov-19-2019
previously i am useing "python pp.py"
RE: modify script. - buran - Nov-19-2019
how you expect to get help on script we haven't seen? post your code in python tags, explain what the goal is.
RE: modify script. - DeaD_EyE - Nov-19-2019
My car doesn't work.
My car is blue.
What is the problem?
RE: modify script. - ChislaineWijdeven - Nov-20-2019
Hi! We really want to help you here, however, you need to show what you're working on (post the code!).
RE: modify script. - MuhammadTalha - Nov-21-2019
hi,
sorry to all, ChislaineWijdeven,DeaD_EyE,Gribouillis, as i am new in this field, so got confused
here is the python script.
#!/usr/bin/python
#
# Post-processing script from of PH data in format used by EPW
# 14/07/2015 - Creation of the script - Samuel Ponce
# 14/03/2018 - Automatically reads the number of q-points - Michael Waters
# 14/03/2018 - Detect if SOC is included in the calculation - Samuel Ponce
# 13/11/2018 - Write dyn files in xml format for SOC case - Shunhong Zhang (USTC)
#
import numpy as np
import os
from xml.dom import minidom
# Convert the dyn files to the xml form, for SOC case - Shunhong Zhang (USTC)
def dyn2xml(prefix):
ndyn=int(os.popen('head -2 {0}.dyn0|tail -1'.format(prefix)).read())
for idyn in range(1,ndyn+1):
print '{0}.dyn{1} to {0}.dyn_q{1}.xml'.format(prefix,idyn)
dynmat=dyn(prefix,idyn)
dynmat._write_xml()
def get_geom_info():
if os.path.isfile('ph.out')==False:
print 'cannot extract geometry info from ph.out'
return 1
else:
volm=float(os.popen('grep -a volume ph.out 2>/dev/null|tail -1').readline().split()[-2])
get_at=os.popen('grep -a -A 3 "crystal axes" ph.out 2>/dev/null|tail -3').readlines()
at=np.array([[float(item) for item in line.split()[3:6]] for line in get_at])
get_bg=os.popen('grep -a -A 3 "reciprocal axes" ph.out 2>/dev/null|tail -3').readlines()
bg=np.array([[float(item) for item in line.split()[3:6]] for line in get_bg])
return volm,at,bg
class dyn(object):
def __init__(self,prefix,idyn):
self._prefix=prefix
self._idyn=idyn
fil='{0}.dyn{1}'.format(prefix,idyn)
f=open(fil)
self._comment=f.readline()
f.readline()
line=f.readline().split()
self._ntype=int(line[0])
self._natom=int(line[1])
self._ibrav=int(line[2])
self._nspin=1
self._cell_dim=np.array([float(ii) for ii in line[3:]])
self._volm=0
self._at=np.zeros((3,3),float)
self._bg=np.zeros((3,3),float)
try: self._volm,self._at,self._bg = get_geom_info()
except: print 'warning: lattice info not found'
self._species=[];
self._mass=[]
for i in range(self._ntype):
line=f.readline().split()
self._species.append(line[1].strip("'"))
self._mass.append(float(line[-1])/911.4442) # normalize to atomic mass
self._atom_type=np.zeros(self._natom,int)
self._pos=np.zeros((self._natom,3),float)
for i in range(self._natom):
line=f.readline().split()
self._atom_type[i]=int(line[1])
for j in range(3): self._pos[i,j]=float(line[j+2])
self._nqpt=int(os.popen('grep -c "Dynamical Matrix" {0}'.format(fil)).read().split()[0])
self._qpt=[]
self._dynmat=np.zeros((self._nqpt,self._natom,self._natom,3,3,2),float)
f.readline()
for iqpt in range(self._nqpt):
f.readline();
f.readline()
line=f.readline().split()
self._qpt.append(np.array([float(item) for item in line[3:6]]))
f.readline()
for i in range(self._natom):
for j in range(self._natom):
f.readline()
data=np.fromfile(f,sep=' ',count=18,dtype=float).reshape(3,3,2)
self._dynmat[iqpt,i,j]=data
self._qpt=np.array(self._qpt)
for i in range(5): f.readline()
self._freq=np.zeros((self._natom*3,2),float)
self._disp=np.zeros((self._natom*3,self._natom,3,2),float)
for i in range(self._natom*3):
line=f.readline().split()
self._freq[i,0]=float(line[4])
self._freq[i,1]=float(line[7])
for j in range(self._natom):
line=f.readline().split()[1:-1]
data=np.array([float(item) for item in line]).reshape(3,2)
self._disp[i,j]=data
def _write_xml(self):
doc=minidom.Document()
root = doc.createElement('Root')
doc.appendChild(root)
geom_info=doc.createElement('GEOMETRY_INFO')
tags=('NUMBER_OF_TYPES','NUMBER_OF_ATOMS','BRAVAIS_LATTICE_INDEX','SPIN_COMPONENTS')
numbers=(self._ntype,self._natom,self._ibrav,self._nspin)
for i,(tag,num) in enumerate(zip(tags,numbers)):
inode=doc.createElement(tag)
inode.setAttribute('type','integer')
inode.setAttribute('size','1')
inode.text=num
inode.appendChild(doc.createTextNode(str(num)))
geom_info.appendChild(inode)
cell_dim=doc.createElement('CELL_DIMENSIONS')
cell_dim.setAttribute('type','real')
cell_dim.setAttribute('size','6')
for i in range(6):
cell_dim.appendChild(doc.createTextNode('{0:16.10f}'.format(self._cell_dim[i])))
geom_info.appendChild(cell_dim)
tags=['AT','BG']
for tag,lat in zip(tags,(self._at,self._bg)):
inode=doc.createElement(tag)
inode.setAttribute('type','real')
inode.setAttribute('size','9')
inode.setAttribute('columns','3')
for i in range(3):
text=' '.join(['{0:16.10f}'.format(item) for item in lat[i]])
inode.appendChild(doc.createTextNode(text))
geom_info.appendChild(inode)
volm=doc.createElement('UNIT_CELL_VOLUME_AU')
volm.setAttribute('type','real')
volm.setAttribute('size','1')
volm.appendChild(doc.createTextNode('{0:16.10f}'.format(self._volm)))
geom_info.appendChild(volm)
for itype in range(self._ntype):
nt=doc.createElement('TYPE_NAME.{0}'.format(itype+1))
nt.setAttribute('type','character')
nt.setAttribute('size','1')
nt.setAttribute('len','3')
nt.appendChild(doc.createTextNode('{0}'.format(self._species[itype])))
na=doc.createElement('MASS.{0}'.format(itype+1))
na.setAttribute('type','real')
na.setAttribute('size','1')
na.appendChild(doc.createTextNode('{0:16.10f}'.format(self._mass[itype])))
geom_info.appendChild(nt)
geom_info.appendChild(na)
for iat in range(self._natom):
at=doc.createElement('ATOM.{0}'.format(iat+1))
at.setAttribute('SPECIES','{0}'.format(self._species[self._atom_type[iat]-1]))
at.setAttribute('INDEX',str(iat+1))
pos=' '.join(['{0:16.10f}'.format(item) for item in self._pos[iat]])
at.setAttribute('TAU',pos)
geom_info.appendChild(at)
nqpt=doc.createElement('NUMBER_OF_Q')
nqpt.setAttribute('type','integer')
nqpt.setAttribute('size','1')
nqpt.appendChild(doc.createTextNode(str(self._nqpt)))
geom_info.appendChild(nqpt)
root.appendChild(geom_info)
for iqpt in range(self._nqpt):
dynmat=doc.createElement('DYNAMICAL_MAT_.{0}'.format(iqpt+1))
qpt=doc.createElement('Q_POINT')
qpt.setAttribute('type','real')
qpt.setAttribute('size','3')
qpt.setAttribute('columns','3')
tnode=doc.createTextNode(' '.join(['{0:16.10f}'.format(item) for item in self._qpt[iqpt]]))
qpt.appendChild(tnode)
dynmat.appendChild(qpt)
for iat in range(self._natom):
for jat in range(self._natom):
ph=doc.createElement('PHI.{0}.{1}'.format(iat+1,jat+1))
ph.setAttribute('type','complex')
ph.setAttribute('size','9')
ph.setAttribute('columns','3')
for i in range(3):
for j in range(3):
text='{0:16.10f} {1:16.10f}'.format(self._dynmat[iqpt,iat,jat,i,j,0],self._dynmat[iqpt,iat,jat,i,j,1])
ph.appendChild(doc.createTextNode(text))
dynmat.appendChild(ph)
root.appendChild(dynmat)
mode=doc.createElement('FREQUENCIES_THZ_CMM1')
for iomega in range(self._natom*3):
inode=doc.createElement('OMEGA.{0}'.format(iomega+1))
inode.setAttribute('type','real')
inode.setAttribute('size','2')
inode.setAttribute('columns','2')
inode.appendChild(doc.createTextNode('{0:16.10f} {1:16.10f}'.format(self._freq[iomega,0],self._freq[iomega,1])))
idisp=doc.createElement('DISPLACEMENT.{0}'.format(iomega+1))
idisp.setAttribute('tpye','complex')
idisp.setAttribute('size','3')
for iat in range(self._natom):
for j in range(3):
tnode=doc.createTextNode('{0:16.10f} {1:16.10f}'.format(self._disp[iomega,iat,j,0],self._disp[iomega,iat,j,1]))
idisp.appendChild(tnode)
mode.appendChild(inode)
mode.appendChild(idisp)
root.appendChild(mode)
fp = open('{0}.dyn_q{1}.xml'.format(self._prefix,self._idyn), 'w')
doc.writexml(fp, addindent=' ', newl='\n')
# Return the number of q-points in the IBZ
def get_nqpt(prefix):
fname = '_ph0/' +prefix+'.phsave/control_ph.xml'
fid = open(fname,'r')
lines = fid.readlines() # these files are relatively small so reading the whole thing shouldn't be an issue
fid.close()
line_number_of_nqpt = 0
while 'NUMBER_OF_Q_POINTS' not in lines[line_number_of_nqpt]: # increment to line of interest
line_number_of_nqpt +=1
line_number_of_nqpt +=1 # its on the next line after that text
nqpt = int(lines[line_number_of_nqpt])
return nqpt
# Check if the calculation include SOC
def hasSOC(prefix):
fname = prefix+'.save/data-file-schema.xml'
xmldoc = minidom.parse(fname)
item = xmldoc.getElementsByTagName('spinorbit')[0]
lSOC = item.childNodes[0].data
return lSOC
# Check if the calculation was done in sequential
def isSEQ(prefix):
fname = '_ph0/'+str(prefix)+'.dvscf'
if (os.path.isfile(fname)):
lseq = True
else:
lseq = False
return lseq
# Enter the number of irr. q-points
user_input = raw_input('Enter the prefix used for PH calculations (e.g. diam)\n')
prefix = str(user_input)
# Test if SOC
SOC = hasSOC(prefix)
# If SOC detected, but dyn is not in XML and we want to convert it
if SOC=='true':
user_input = raw_input('Calculation with SOC detected. Do you want to convert dyn in XML format [y/n]?\n')
if str(user_input) == 'y':
dyn2xml(prefix)
os.system('mv {0}.dyn*.xml save'.format(prefix))
# If no SOC detected, do you want to convert into XML format
if SOC=='false':
user_input = raw_input('Calculation without SOC detected. Do you want to convert to xml anyway [y/n]?\n')
if str(user_input) == 'y':
SOC = 'true'
dyn2xml(prefix)
os.system('mv {0}.dyn*.xml save'.format(prefix))
# Test if seq. or parallel run
SEQ = isSEQ(prefix)
if True: # this gets the nqpt from the outputfiles
nqpt = get_nqpt(prefix)
else:
# Enter the number of irr. q-points
user_input = raw_input('Enter the number of irreducible q-points\n')
nqpt = user_input
try:
nqpt = int(user_input)
except ValueError:
raise Exception('The value you enter is not an integer!')
os.system('mkdir save 2>/dev/null')
for iqpt in np.arange(1,nqpt+1):
label = str(iqpt)
# Case calculation in seq.
if SEQ:
# Case with SOC
if SOC == 'true':
os.system('cp '+prefix+'.dyn0 '+prefix+'.dyn0.xml')
os.system('cp '+prefix+'.dyn'+str(iqpt)+'.xml save/'+prefix+'.dyn_q'+label+'.xml')
if (iqpt == 1):
os.system('cp _ph0/'+prefix+'.dvscf* save/'+prefix+'.dvscf_q'+label)
os.system('cp -r _ph0/'+prefix+'.phsave save/')
os.system('cp '+prefix+'.fc.xml save/ifc.q2r.xml')
else:
os.system('cp _ph0/'+prefix+'.q_'+str(iqpt)+'/'+prefix+'.dvscf* save/'+prefix+'.dvscf_q'+label)
os.system('rm _ph0/'+prefix+'.q_'+str(iqpt)+'/*wfc*' )
# Case without SOC
if SOC == 'false':
os.system('cp '+prefix+'.dyn'+str(iqpt)+' save/'+prefix+'.dyn_q'+label)
if (iqpt == 1):
os.system('cp _ph0/'+prefix+'.dvscf save/'+prefix+'.dvscf_q'+label)
os.system('cp -r _ph0/'+prefix+'.phsave save/')
os.system('cp '+prefix+'.fc save/ifc.q2r')
else:
os.system('cp _ph0/'+prefix+'.q_'+str(iqpt)+'/'+prefix+'.dvscf save/'+prefix+'.dvscf_q'+label)
os.system('rm _ph0/'+prefix+'.q_'+str(iqpt)+'/*wfc*' )
else:
# Case with SOC
if SOC == 'true':
os.system('cp '+prefix+'.dyn0 '+prefix+'.dyn0.xml')
os.system('cp '+prefix+'.dyn'+str(iqpt)+'.xml save/'+prefix+'.dyn_q'+label+'.xml')
if (iqpt == 1):
os.system('cp _ph0/'+prefix+'.dvscf1 save/'+prefix+'.dvscf_q'+label)
os.system('cp -r _ph0/'+prefix+'.phsave save/')
os.system('cp '+prefix+'.fc.xml save/ifc.q2r.xml')
else:
os.system('cp _ph0/'+prefix+'.q_'+str(iqpt)+'/'+prefix+'.dvscf1 save/'+prefix+'.dvscf_q'+label)
os.system('rm _ph0/'+prefix+'.q_'+str(iqpt)+'/*wfc*' )
# Case without SOC
if SOC == 'false':
os.system('cp '+prefix+'.dyn'+str(iqpt)+' save/'+prefix+'.dyn_q'+label)
if (iqpt == 1):
os.system('cp _ph0/'+prefix+'.dvscf1 save/'+prefix+'.dvscf_q'+label)
os.system('cp -r _ph0/'+prefix+'.phsave save/')
os.system('cp '+prefix+'.fc save/ifc.q2r')
else:
os.system('cp _ph0/'+prefix+'.q_'+str(iqpt)+'/'+prefix+'.dvscf1 save/'+prefix+'.dvscf_q'+label)
os.system('rm _ph0/'+prefix+'.q_'+str(iqpt)+'/*wfc*' )
RE: modify script. - DeaD_EyE - Nov-21-2019
Holy.... i've never seen so often the use of
os.system.This whole script could be refactored to an easier version, without calling
the whole time foreign tools on the system.
Also pathlib allows a better abstraction of paths and has useful methods.
Using the tools
tail, head, grep inside a Python program, is not very good.In this case the program relies on the os and the existence of this tools (yes, they are standard on all Linux Systems).
If you want to look for another files, just change the
.dvscf1 to something else and hope that it works.Otherwise you have to dive deep into the code and first you should correct the sins from the other guys who made it.
RE: modify script. - MuhammadTalha - Nov-21-2019
Hi DeaD_EyE,
so nice of u. i have tried changing .dvscf1 before asking here but its giving same error.
you are right i have to contact with the person that wrote this codding.
so nice of u for your kindness :)
RE: modify script. - Gribouillis - Nov-21-2019
MuhammadTalha Wrote:i have tried changing .dvscf1 before asking here but its giving same error.Please post the whole error message sent by python, it would help a lot identifying the problem.