From 090f3054e4aa67ea3c38287712412c35e208d9e0 Mon Sep 17 00:00:00 2001
From: giannozz <giannozz@c92efa57-630b-4861-b058-cf58834340f0>
Date: Mon, 29 Oct 2012 06:24:45 +0000
Subject: [PATCH] Converter from CIF format to QE input data, and from CASTEP
 data. Both are experimental and little-tested. Courtesy of Carlo Nervi

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@9590 c92efa57-630b-4861-b058-cf58834340f0
---
 PW/tools/castep2qe.sh | 136 ++++++++++++++++++++++++++++++
 PW/tools/cif2qe.sh    | 192 ++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 328 insertions(+)
 create mode 100755 PW/tools/castep2qe.sh
 create mode 100755 PW/tools/cif2qe.sh

diff --git a/PW/tools/castep2qe.sh b/PW/tools/castep2qe.sh
new file mode 100755
index 000000000..0b0a9676c
--- /dev/null
+++ b/PW/tools/castep2qe.sh
@@ -0,0 +1,136 @@
+!
+! Copyright (C) 2012 Carlo Nervi
+! This file is distributed under the terms of the
+! GNU General Public License. See the file `License'
+! in the root directory of the present distribution,
+! or http://www.gnu.org/copyleft/gpl.txt .
+!
+#!/bin/bash
+
+if [ $# != 1 ]; then
+ printf "Usage: castep2qe.sh File (.cell extension not required!)\n Requires File.cell\n"
+ exit
+fi
+
+if [ ! -f $1.cell ]; then
+ echo "Error. Cannot find file $1.cell"
+ exit
+fi
+
+#
+# 1 Bohr = 0.529177249 Angstrom
+#
+
+awk -v FILE="$1" '
+
+BEGIN {bohr = 0.52917720859;}
+
+/%BLOCK LATTICE_CART/ {
+  for (i=0; i<3; i++) {
+    getline;
+    cell_px[i] = $1;
+    cell_py[i] = $2;
+    cell_pz[i] = $3;
+  }
+ }
+
+/%BLOCK POSITIONS_FRAC/ {
+  natom=0;
+  getline;
+  do {
+    frac[natom++] = $0;
+    getline;
+  } while ($0 !~ /%ENDBLOCK POSITIONS_FRAC/)
+}
+
+/%BLOCK KPOINTS_LIST/ {
+  nks=0;
+  getline;
+  do {
+    kpout[nks++] = sprintf ("   %19.15f   %19.15f   %19.15f     %19.15f", $1, $2, $3, $4);
+    getline;
+  } while ($0 !~ /%ENDBLOCK KPOINTS_LIST/)
+}
+
+/%BLOCK SPECIES_MASS/ {
+  nspec=0;
+  getline;
+  do {
+    species_symb[nspec] = $1;
+    species_mass[nspec] = $2;
+    nspec++;
+    getline;
+  } while ($0 !~ /%ENDBLOCK SPECIES_MASS/)
+}
+
+/%BLOCK SPECIES_POT/ {
+  npot=0;
+  getline;
+  do {
+    pot_symb[npot] = $1;
+    pot_ps[npot] = $2;
+    npot++;
+    getline;
+  } while ($0 !~ /%ENDBLOCK SPECIES_POT/)
+}
+
+END {
+print "&CONTROL"
+print "                       title = \x027" FILE "\x027"
+print "                 calculation = \x027" "relax\x027"
+print "                restart_mode = \x027" "from_scratch\x027"
+print "                      outdir = \x27" "./1\x027"
+print "                  pseudo_dir = \x027" "../PP/atompaw\x027"
+print "                      prefix = \x027" FILE "\x027"
+print "                     disk_io = \x027" "none\x027"
+print "                   verbosity = \x027" "default\x027"
+print "               etot_conv_thr = 0.00001"
+print "               forc_conv_thr = 0.0001"
+print "                       nstep = 680"
+print "                     tstress = .true."
+print "                     tprnfor = .true."
+print " /"
+print " &SYSTEM"
+print "                       ibrav = 0"
+print "                         nat = " natom
+print "                        ntyp = " nspec
+print "                     ecutwfc = 60"
+print "                     ecutrho = 600"
+print "!                      london = .true."
+print " /"
+print " &ELECTRONS"
+print "            electron_maxstep = 200"
+print "                    conv_thr = 1.0D-7"
+print "              diago_thr_init = 1e-4"
+print "                 startingpot = \x027" "atomic\x027"
+print "                 startingwfc = \x027" "atomic\x027"
+print "                 mixing_mode = \x027" "plain\x027"
+print "                 mixing_beta = 0.5"
+print "                 mixing_ndim = 8"
+print "             diagonalization = \x027" "david\x027"
+print " /"
+print "&IONS"
+print "                ion_dynamics = \x027" "bfgs\x027"
+print " /"
+print "\n"
+print "ATOMIC_SPECIES"
+if (npot != nspec) {
+  print "Error; npot [" npot "] != nspec [" nspec "]"
+  exit;
+}
+for (i=0; i< npot; i++) printf "  %3s  %14.10f  %s\n", species_symb[i], species_mass[i], pot_ps[i];
+
+print "\nATOMIC_POSITIONS crystal"
+for (i=0; i< natom; i++) print frac[i];
+
+print "\nK_POINTS crystal"
+print nks;
+for (i=0; i<nks; i++) print kpout[i];
+
+print "\nCELL_PARAMETERS"
+for (i=0; i<3; i++) printf "  %19.15f   %19.15f   %19.15f\n", cell_px[i]/bohr, cell_py[i]/bohr, cell_pz[i]/bohr;
+
+print "\n\n\n"
+
+}
+' $1.cell
diff --git a/PW/tools/cif2qe.sh b/PW/tools/cif2qe.sh
new file mode 100755
index 000000000..4984da303
--- /dev/null
+++ b/PW/tools/cif2qe.sh
@@ -0,0 +1,192 @@
+!
+! Copyright (C) 2012 Carlo Nervi
+! This file is distributed under the terms of the
+! GNU General Public License. See the file `License'
+! in the root directory of the present distribution,
+! or http://www.gnu.org/copyleft/gpl.txt .
+!
+#!/bin/bash
+
+#
+# Convertitore da file cif a Quantum Espresso (C) Carlo Nervi 2012
+#
+#
+
+if [ $# != 1 ]; then
+ printf "Usage: cif2qe.sh File (.cif extension not required!)\n Requires File.cif\n"
+ exit
+fi
+
+if [ ! -f $1.cif ]; then
+ echo "Error. Cannot find file $1.cif"
+ exit
+fi
+
+awk -v FILE="$1" '
+
+BEGIN {
+ bohr = 0.52917720859
+ nfield=split("H He Li Be B C N O F Ne Na Mg Al Si P S Cl Ar K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe Cs Ba La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn Fr Ra Ac Th Pa U Np Pu Am Cm Bk Cf Es Fm Md No Lr Rf Db Sg Bh Hs Mt", AtomSymb, " ")
+ split("1.0079 4.0026 6.941 9.0122 10.811 12.0107 14.0067 15.9994 18.9984 20.1797 22.9897 24.305 26.9815 28.0855 30.9738 32.065 35.453 39.948 39.0983 40.078 44.9559 47.867 50.9415 51.9961 54.938 55.845 58.9332 58.6934 63.546 65.39 69.723 72.64 74.9216 78.96 79.904 83.8 85.4678 87.62 88.9059 91.224 92.9064 95.94 98 101.07 102.906 106.42 107.868 112.411 114.818 118.71 121.76 127.6 126.904 131.293 132.905 137.327 138.905 140.116 140.908 144.24 145 150.36 151.964 157.25 158.925 162.5 164.93 167.259 168.934 173.04 174.967 178.49 180.948 183.84 186.207 190.23 192.217 195.078 196.966 200.59 204.383 207.2 208.98 209 210 222 223 226 227 232.038 231.036 238.029 237 244 243 247 247 251 252 257 258 259 262 261 262 266 264 277 268", AtomMass, " ")
+ for (i=1; i<=nfield; i++) Atoms[AtomSymb[i]] = AtomMass[i]
+}
+
+function eval(cmd,expression) {
+ expression = "awk \047BEGIN{printf \"%19.15f\", "cmd" }\047"
+ expression|getline l
+ return l     
+}
+
+function norma(pat, repl, str) {
+ gsub(pat, repl, str)
+ gsub(/--/, "+", str)
+ return eval(str)
+}
+
+# Salta il commento
+/^\#/ { next
+}
+
+/loop_/ { loop_switch=1; next }
+
+
+/^_/ {
+ jvar=0
+ if (loop_switch==1) {
+  Var[ivar]=$1
+  ivar++
+  Num_Var=ivar
+# Debug:
+#  print "D: " $1 "=" ivar
+  if ($1 ~ /^_symmetry/) nsymm=0
+  if ($1 ~ /^_atom/) natom=0
+#  if ($1 == "_symmetry_equiv_pos_as_xyz") I_Symm=ivar
+ } else {
+   ivar=0
+   loop_switch=0
+   Var2[$1]=$2
+ }
+}
+
+/^[^_]/ {
+ ivar=0
+ if (loop_switch==1 || loop_switch==2) {
+  loop_switch=2
+  for (i=0; i<Num_Var; i++) {
+    Array[Var[i]][jvar]=$(i+1)
+# Debug:
+#    printf "Array[" Var[i] "][" jvar "]=  [" $(i+1) "]   "
+  }
+# Debug:
+#  printf "\n"
+  jvar++
+  if (Var[0] ~ /^_atom/) natom++
+  if (Var[0] ~ /^_symmetry/) nsymm++
+ } else {
+  jvar=0
+  loop_switch=0
+ }
+}
+
+END {
+ a=strtonum(Var2["_cell_length_a"])
+ b=strtonum(Var2["_cell_length_b"])
+ c=strtonum(Var2["_cell_length_c"])
+ alpha=strtonum(Var2["_cell_angle_alpha"])/180.0*3.14159265359
+ beta=strtonum(Var2["_cell_angle_beta"])/180.0*3.14159265359
+ gamma=strtonum(Var2["_cell_angle_gamma"])/180.0*3.14159265359
+
+ ntyp=0
+ for (i=0; i<natom; i++) {
+  if (Type_Atom[Array["_atom_site_type_symbol"][i]] != 1) {
+    Type_Atom[Array["_atom_site_type_symbol"][i]]=1
+    AtomTyp[ntyp]=Array["_atom_site_type_symbol"][i]
+# Debug:
+#    print "D: [" Array["_atom_site_type_symbol"][i] "] [" Atoms[Array["_atom_site_type_symbol"][i]] "]  AtomTyp[" ntyp "]=[" AtomTyp[ntyp] "]"
+    ntyp++
+  }
+ }
+ print "&CONTROL"
+ print "                       title = \x027" FILE "\x027"
+# print "!   _symmetry_space_group_name_H-M = " Var2["_symmetry_space_group"]
+# print "!   _symmetry_Int_Tables_number = " Var2["_symmetry_Int_Tables_number"]
+# print "!   _symmetry_cell_setting = " Var2["_symmetry_cell_setting"]
+ print "                 calculation = \x027" "relax\x027"
+ print "                restart_mode = \x027" "from_scratch\x027"
+ print "                      outdir = \x027" "./1\x027"
+ print "                  pseudo_dir = \x027" "../PP/atompaw\x027"
+ print "                      prefix = \x027" FILE "\x027"
+ print "                     disk_io = \x027" "none\x027"
+ print "                   verbosity = \x027" "default\x027"
+ print "               etot_conv_thr = 0.00001"
+ print "               forc_conv_thr = 0.0001"
+ print "                       nstep = 680"
+ print "                     tstress = .true."
+ print "                     tprnfor = .true."
+ print " /"
+ print " &SYSTEM"
+ print "                       ibrav = 0"
+ print "                         nat = " natom*nsymm
+ print "                        ntyp = " ntyp
+ print "                     ecutwfc = 60"
+ print "                     ecutrho = 600"
+ print "!                      london = .true."
+ print " /"
+ print " &ELECTRONS"
+ print "            electron_maxstep = 200"
+ print "                    conv_thr = 1.0D-7"
+ print "              diago_thr_init = 1e-4"
+ print "                 startingpot = \x027" "atomic\x027"
+ print "                 startingwfc = \x027" "atomic\x027"
+ print "                 mixing_mode = \x027" "plain\x027"
+ print "                 mixing_beta = 0.5"
+ print "                 mixing_ndim = 8"
+ print "             diagonalization = \x027" "david\x027"
+ print " /"
+ print "&IONS"
+ print "                ion_dynamics = \x027" "bfgs\x027"
+ print " /"
+ print "\n"
+ print "ATOMIC_SPECIES"
+ for (i=0; i<ntyp; i++) printf "  %3s  %14.10f  %s.pbe-van_ak.UPF\n", AtomTyp[i], Atoms[AtomTyp[i]], AtomTyp[i];
+
+ print "\nATOMIC_POSITIONS crystal"
+ for (j=0; j<nsymm; j++) {
+  for (i=0; i<natom; i++) {
+    if (split(Array["_symmetry_equiv_pos_as_xyz"][j], Tmp, ",") != 3) {
+      print "Error in _symmetry_equiv_pos_as_xyz. Number of fields != 3: [1]=" Tmp[1] " [2]=" Tmp[2] " [3]=" Tmp[3]
+      print "D: " Array["_symmetry_equiv_pos_as_xyz"][j] "  " Tmp[1] "  " Tmp[2] "  " Tmp[3]
+      exit
+    }
+    p_X=norma("x", Array["_atom_site_fract_x"][i], Tmp[1])
+    p_Y=norma("y", Array["_atom_site_fract_y"][i], Tmp[2])
+    p_Z=norma("z", Array["_atom_site_fract_z"][i], Tmp[3])
+    
+    printf "%2s   %19.15f   %19.15f   %19.15f\n", Array["_atom_site_type_symbol"][i], p_X, p_Y, p_Z
+  }
+ }
+
+ print "\nK_POINTS crystal"
+ print nks;
+ for (i=0; i<nks; i++) print kpout[i];
+
+ cell_px[0]=a
+ cell_py[0]=0.0
+ cell_pz[0]=0.0
+ cell_px[1]=b*cos(gamma)
+ cell_py[1]=b*sin(gamma)
+ cell_pz[1]=0.0
+ cell_px[2]=c*cos(beta)
+ cell_py[2]=c*(cos(alpha)-cos(beta)*cos(gamma))/sin(gamma)
+ cell_pz[2]=c*sqrt(1.0 - cos(alpha)^2 - cos(beta)^2 - cos(gamma)^2 + 2*cos(alpha)*cos(beta)*cos(gamma))/sin(gamma)
+
+ print "\nCELL_PARAMETERS"
+# Debug:
+# print "a=" a "  b=" b "  c=" c "   alfa=" alpha "   beta=" beta "   gamma=" gamma
+# for (i=0; i<3; i++) printf "  %19.15f   %19.15f   %19.15f\n", cell_px[i], cell_py[i], cell_pz[i];
+ for (i=0; i<3; i++) printf "  %19.15f   %19.15f   %19.15f\n", cell_px[i]/bohr, cell_py[i]/bohr, cell_pz[i]/bohr;
+
+ print "\n\n"
+
+} ' $1.cif
+