Editing Some scripts for file format conversion

QE use the xsf file format produced by pwi2xsf.sh and pwo2xsf.sh. They can be opened only by Xcrysden. If you need to open your structure with other tools you can convert it to a cif file using the BABEL converter (http://openbabel.org/wiki/Main_Page) the only problem is that they do it wrong... the converter is unable to write the atom elements of the xsf files to the cif file.

Here is a small script that correct the errors of the wrong cif file obtained with BABEL.

<pre>
#!/usr/bin/python2
#This script use the wrong output obtained from Babel when it tries to convert a xsf file to a cif file and convert it to a correct cif file 
#tha can be opened by Material studio
import os
import sys
import math

#function to count number of lines in file
def file_len(fname):
    with open(fname,"r") as f:
        for i, l in enumerate(f):
            pass
    return i + 1
    f.close

## program start here ##
#EDIT ONLY HERE (below)
name_prefix="ZnFeO2_spinel_001.out2"
#EDIT ONLY HERE (above)

filename=name_prefix+".xsf" #correct xsf 
filename1=name_prefix+".xsf.cif1" #wrong cif from Babel
output=name_prefix+".xsf.cif"   #right cif at the output 


#Step 1 read from xsf the Element names
nmax=file_len(filename)
f0=open(filename,"r")
i=0 #used as line counter
T=0 #used as logical switch
k=0 #used for PRIMVEC
A=[] #store name
while i<nmax:
    line=f0.readline().split()
    #print line
    if 'PRIMVEC' in line:
	i=i+1
	T=1
    elif 'PRIMCOORD' in line:
	i=i+1
	T=2
    elif T==1:
	if k==0:
	    a=float(line[0])
	    b=a
	    k=k+1
	elif k==1:
	    k=k+1
	elif k==2:
	    c=float(line[2])
	    k=k+1
	    T=0
	i=i+1   
    elif T==2:
	numatoms=int(line[0])
	T=3
	i=i+1
    elif T==3:
	A.append([line[0]])
	i=i+1
    else:
	i=i+1
f0.close

#Step 2 read the wrong cif file and write the right cif file
f1=open(filename1,"r")
out=open(output,"w")
i=0 #used as line counter
while i<1:
    line=f1.readline()
    print >> out, line,
    if '_atom_site_occupancy' in line:
        i=1 #coordinates start here, exit while loop
for i in range(numatoms):
    line=f1.readline().split()
    print >> out, "    "+str(A[i][0])+str(i), A[i][0], line[2], line[3], line[4], line[5]
f1.close
out.close
</pre>

Lazy script to convert cif file to QE input (must be updated)
<pre>
#!/usr/bin/python2
#This script prepare the input file for quantum espresso ATOMIC_POSITIONS from a .cif file for a geometry optimisation
#Not for cell optimisation. Replace relax by vc-relax for that and add the CELL section
import os
import sys
import math

#function to count number of lines in file
def file_len(fname):
    with open(fname,"r") as f:
        for i, l in enumerate(f):
            pass
    return i + 1
    f.close

## program start here ##
name_prefix="ZnFeO2_spinel_super222_001"

filename="fake.list" #changing the whole code is more bothersome than rewriting the whole code...
filename1=name_prefix+".cif"
#xyzname=name_prefix+".xyz"   #not needed 
output=name_prefix+".in"

#xyz=open(xyzname,"r")
#numatoms=int(xyz.readline())
#xyz.close
#numatoms=128  #Is now extracted from xyz file
#count number of atoms in cif file
f1=open(filename1,"r")
nmax1=file_len(filename1)
i=0
j=0
shift=1 #by default input list start at 1
for line in f1:
    if '_atom_site_occupancy' in line:
        i=1 #coordinates start here
    elif i == 1:
        j=j+1
    if 'Fe0' in line:
        print 'Fe0 found' 
        shift=0 #sometimes the list start at 0. I do not know why
    #print i,j,numatoms
f1.close
numatoms=j #taken from the cif file, xyz file not needed
print numatoms,shift
nmax=file_len(filename)
#print nmax

#create the list
uplist=[]
downlist=[]
freezelist=[]
T=0 #used for list type
i=0 #used for line index
f=open(filename,"r") #should be read-only

while i<nmax:
    line=f.readline().split()
    #print T
    #print line
    if 'up' in line:
        i=i+1
        T=1
    elif 'down' in line:
        i=i+1
        T=2
    elif 'freeze' in line:
        i=i+1
        T=3
    elif T==1:
        b=[int(l) for l in line]
        uplist = uplist + b
        i=i+1
    elif T==2:
        b=[int(l) for l in line]
        downlist = downlist + b
        i=i+1
    elif T==3:
        b=[int(l) for l in line]
        freezelist = freezelist + b
        i=i+1
    else:
        i=i+1
    #print i

f.close()

#Check if number is in both up & down list
for i in range(len(uplist)):
    if uplist[i] in downlist:
        print "BUG in list"


#print uplist
#print downlist
#print freezelist
f1=open(filename1,"r")
nmax1=file_len(filename1)
i=0
while i<1:
    line=f1.readline()
    if '_cell_length_a' in line:
        a=float(line.split()[1])
    if '_cell_length_b' in line:
        b=float(line.split()[1])
    if '_cell_length_c' in line:
        c=float(line.split()[1])
    if '_cell_angle_alpha' in line:
        cosBC=math.cos(math.radians(float(line.split()[1])))
        if -1.e-15 <cosBC < 1.e-15:
            cosBC=0
    if '_cell_angle_beta' in line:
        cosAC=math.cos(math.radians(float(line.split()[1])))
        if -1.e-15 < cosAC < 1.e-15:
            cosAC=0
    if '_cell_angle_gamma' in line:
        cosAB=math.cos(math.radians(float(line.split()[1])))
        if -1.e-15 < cosAB < 1.e-15:
            cosAB=0
    if '_atom_site_occupancy' in line:
        i=1 #coordinates start here, exit while loop
out=open(output,"w")
print >> out, " &CONTROL"
print >> out, "                 calculation = 'relax' ,"
print >> out, "                restart_mode = 'from_scratch' ,"
print >> out, "                  wf_collect = .true. ,"
print >> out, "                      prefix =","'"+name_prefix+"'", ","
print >> out, "                   verbosity = 'high' ,"
print >> out, "               etot_conv_thr = 1.0D-8 ,"
print >> out, "               forc_conv_thr = 1.0D-4 ,"
print >> out, "                     tstress = .true. ,"
print >> out, "                     tprnfor = .true. ,"
print >> out, "                    dipfield = .true. ,"
print >> out, "                       nstep = 400 ,"
print >> out, "                 max_seconds = 171000 , !47h30"
print >> out, " /"
print >> out, "&SYSTEM"
print >> out, "                      ibrav = 14,"
print >> out, "                          A =",a, ","
print >> out, "                          B =",b, ","
print >> out, "                          C =",c, ","
print >> out, "                      cosAB =",cosAB, ","
print >> out, "                      cosAC =",cosAC, ","
print >> out, "                      cosBC =",cosBC, ","
print >> out, "                        nat =",numatoms,","
print >> out, "                       ntyp = 3,"
print >> out, "                    ecutwfc = 40.0 , !J. Chem. Phys. 138, 194709"  #34
print >> out, "                    ecutrho = 320.0 ,!J. Chem. Phys. 138, 194709" #323
print >> out, "                  input_dft = 'PBE' ,"
print >> out, "                occupations = 'smearing' ,"
print >> out, "                    degauss = 0.02 ,"
print >> out, "                   smearing = 'gaussian' ,"
print >> out, "                      nspin = 2 ,"
print >> out, "    starting_magnetization(1) = 0.0,"
print >> out, "    starting_magnetization(2) = -0.0,"
print >> out, "    starting_magnetization(3) = 0.0,"
print >> out, "                  lda_plus_u = .true. ,"
print >> out, "             lda_plus_u_kind = 0 , !type 1 not implemented, ignores J"
print >> out, "                Hubbard_U(1) = 0.0,"
print >> out, "                Hubbard_U(2) = 4.2," 
print >> out, "               Hubbard_J0(1) = 0.0, !Is ignored!"
print >> out, "               Hubbard_J0(2) = 0.0, !Is ignored!"
print >> out, "                   !vdw_corr = Grimme-D2, ! DFT-D2 dispersion correction"
print >> out, "/"
print >> out, "&ELECTRONS"
print >> out, "           electron_maxstep = 800,"
print >> out, "              conv_thr_init = 1e-4 ,"
print >> out, "                   conv_thr = 1e-9 ,"
print >> out, "                startingpot = 'atomic' ,"
print >> out, "                startingwfc = 'random' ,"
print >> out, "               adaptive_thr = .true. ,"
print >> out, "                mixing_beta = 0.514,"
print >> out, "            diagonalization = 'david' ,"
print >> out, "/"
print >> out, "&IONS"
print >> out, "               ion_dynamics = 'bfgs' ,"
print >> out, "           trust_radius_min =  1.D-5 ,"
print >> out, "/"
print >> out, "ATOMIC_SPECIES"
print >> out, "  Zn   65.39000  Zn_pbe_v1.uspp.F.UPF"
print >> out, "  Fe   55.84500  Fe.pbe-sp-van_ak.UPF" 
print >> out, "   O   15.99990  O.pbe-van_ak.UPF"
print >> out, "ATOMIC_POSITIONS crystal"
for i in range(numatoms):
    line=f1.readline().split()
    if (i+shift) in uplist and (i+shift) not in freezelist: 
        print >> out, line[1]+"1", line[2], line[3], line[4] 
    elif (i+shift) in uplist and (i+shift) in freezelist: 
        print >> out, line[1]+"1", line[2], line[3], line[4], 0 ,0, 0
    elif (i+shift) in downlist and (i+shift) not in freezelist: 
        print >> out, line[1]+"2", line[2], line[3], line[4], 1 ,1, 1
    elif (i+shift) in downlist and (i+shift) in freezelist: 
        print >> out, line[1]+"2", line[2], line[3], line[4], 0 ,0, 0
    else:
        if (i+shift) in freezelist: 
            print >> out, line[1], line[2], line[3], line[4], 0 ,0, 0
        else:
            print >> out, line[1], line[2], line[3], line[4] 
print >> out, "K_POINTS Gamma"
#closing file
f1.close
f.close
out.close
</pre>