/ase/calculators/siesta/mbpt_lcao_utils.py

https://gitlab.com/ansobolev/ase · Python · 286 lines · 164 code · 65 blank · 57 comment · 36 complexity · 3d2af8f21b30969cc76285f9f4a72bff MD5 · raw file

  1. from __future__ import division
  2. import numpy as np
  3. import re
  4. def read_file(fname):
  5. """
  6. read the file fname and return a list of the lines.
  7. """
  8. f = open(fname, 'r')
  9. LINE = list()
  10. for line in f:
  11. LINE.append(line)
  12. return LINE
  13. def delete_blanc(L):
  14. """
  15. delete the blank space from a string
  16. """
  17. compt = 0
  18. while L[compt] == ' ' or L[compt] == '\t':
  19. compt = compt + 1
  20. L = L[compt:len(L)]
  21. return L
  22. def read_number(L):
  23. compt = 0
  24. while L[compt] != ' ' and compt < (len(L) - 1):
  25. compt = compt + 1
  26. nb1 = float(L[0:compt + 1])
  27. L = L[compt:len(L)]
  28. return nb1, L
  29. def recover_data_string(fname, string):
  30. L = read_file(fname)
  31. for i in L:
  32. print(i[0:len(string)], string)
  33. if i[0:len(string)] == string:
  34. v = str2float(i)
  35. return v
  36. def save_line(Lnb, Cnb, LINE):
  37. number = LINE[Lnb]
  38. nombre = list()
  39. for i in range(Cnb):
  40. number = delete_blanc(number)
  41. nb1, number = read_number(number)
  42. nombre.append(nb1)
  43. return nombre
  44. def dim_y(Cnb, L):
  45. if len(L) == Cnb * 10 + Cnb * 6 + 1:
  46. nb_col = Cnb
  47. else:
  48. nb_col = Cnb + 1
  49. return nb_col
  50. def read_color_file(fname):
  51. L = read_file(fname)
  52. atom = list()
  53. for i in range(len(L)):
  54. nb = np.array(str2float(L[i]))
  55. species = recover_specie(L[i])
  56. atom.append([species, nb])
  57. return atom
  58. def readSiestaFA(fname):
  59. L = read_file(fname)
  60. Forces = []
  61. for i in range(1, len(L)):
  62. Forces.append([i, np.array(str2float(L[i])[1:4])])
  63. return Forces
  64. def readBasis_spec(fname, nb_specie):
  65. """
  66. Example Basis_specs from siesta output
  67. <basis_specs>
  68. ===============================================================================
  69. H Z= 1 Mass= 1.0100 Charge= 0.17977+309
  70. Lmxo=0 Lmxkb= 2 BasisType=split Semic=F
  71. L=0 Nsemic=0 Cnfigmx=1
  72. n=1 nzeta=2 polorb=1
  73. splnorm: 0.15000
  74. vcte: 0.0000
  75. rinn: 0.0000
  76. rcs: 0.0000 0.0000
  77. lambdas: 1.0000 1.0000
  78. -------------------------------------------------------------------------------
  79. L=0 Nkbl=1 erefs: 0.17977+309
  80. L=1 Nkbl=1 erefs: 0.17977+309
  81. L=2 Nkbl=1 erefs: 0.17977+309
  82. ===============================================================================
  83. </basis_specs>
  84. """
  85. L = read_file(fname)
  86. specie_charac = {}
  87. line = 0
  88. len_basis = len('<basis_specs>')
  89. i = 0
  90. while i < nb_specie:
  91. if L[line][0:len_basis] == '<basis_specs>':
  92. i = i + 1
  93. info = str2float(L[line + 2])
  94. if L[line + 2][1] == ' ':
  95. specie_charac[
  96. L[line + 2][0]] = {'Z': info[0],
  97. 'Mass': info[1], 'Charge': info[2]}
  98. else:
  99. specie_charac[
  100. L[line + 2][0:2]] = {'Z': info[0],
  101. 'Mass': info[1], 'Charge': info[2]}
  102. while L[line][0:len('</basis_specs>')] != '</basis_specs>':
  103. line = line + 1
  104. else:
  105. line = line + 1
  106. return specie_charac
  107. def str2float(string):
  108. numeric_const_pattern = r"""
  109. [-+]? # optional sign
  110. (?:
  111. (?: \d* \. \d+ ) # .1 .12 .123 etc 9.1 etc 98.1 etc
  112. |
  113. (?: \d+ \.? ) # 1. 12. 123. etc 1 12 123 etc
  114. )
  115. # followed by optional exponent part if desired
  116. (?: [Ee] [+-]? \d+ ) ?
  117. """
  118. rx = re.compile(numeric_const_pattern, re.VERBOSE)
  119. nb = rx.findall(string)
  120. for i in enumerate(nb):
  121. nb[i[0]] = float(i[1])
  122. return np.array(nb)
  123. def str2int(string):
  124. numeric_const_pattern = r"""
  125. [-+]? # optional sign
  126. (?:
  127. (?: \d* \. \d+ ) # .1 .12 .123 etc 9.1 etc 98.1 etc
  128. |
  129. (?: \d+ \.? ) # 1. 12. 123. etc 1 12 123 etc
  130. )
  131. # followed by optional exponent part if desired
  132. (?: [Ee] [+-]? \d+ ) ?
  133. """
  134. rx = re.compile(numeric_const_pattern, re.VERBOSE)
  135. nb = rx.findall(string)
  136. for i in enumerate(nb):
  137. nb[i[0]] = int(i[1])
  138. return np.array(nb)
  139. def recover_specie(string):
  140. """
  141. Select specie in a string of caractere from
  142. a .xyz file
  143. Input parameters:
  144. string (str): the string to analyse
  145. Output parameter:
  146. string_p (str): the specie
  147. """
  148. specie = list()
  149. comp = 0
  150. letter = string[0]
  151. if letter == ' ':
  152. while letter == ' ' or comp >= len(string):
  153. letter = string[comp]
  154. comp = comp + 1
  155. while letter != ' ' or comp >= len(string):
  156. letter = string[comp]
  157. specie.append(letter)
  158. comp = comp + 1
  159. else:
  160. while letter != ' ' or comp >= len(string):
  161. letter = string[comp]
  162. specie.append(letter)
  163. comp = comp + 1
  164. specie.remove(' ')
  165. string_p = ''
  166. for i in specie:
  167. string_p = string_p + i
  168. return string_p
  169. def delete_number_string(string):
  170. # not working for exponential expression
  171. nb_list = ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '.', '-']
  172. L = list()
  173. for i in string:
  174. L.append(i)
  175. L_P = list()
  176. for i in enumerate(L):
  177. inside = False
  178. for j in nb_list:
  179. if i[1] == j:
  180. inside = True
  181. if not inside:
  182. L_P.append(i[1])
  183. string_p = ''
  184. for i in L_P:
  185. if i != ' ':
  186. string_p = string_p + i
  187. return string_p
  188. def pol2cross_sec(p, omg):
  189. """
  190. Convert the polarizability in au to cross section in nm**2
  191. INPUT PARAMETERS:
  192. -----------------
  193. p (np array): polarizability from mbpt_lcao calc
  194. omg (np.array): frequency range in eV
  195. OUTPUT_PARAMETERS:
  196. ------------------
  197. sigma (np array): cross section in nm**2
  198. """
  199. c = 137 # speed of the light in au
  200. omg = omg * 0.036749309 # to convert from eV to Hartree
  201. sigma = 4 * np.pi * omg * p / (c) # bohr**2
  202. sigma = sigma * (0.052917725)**2 # nm**2
  203. return sigma
  204. def interpolate(x, y, nb_pts):
  205. """
  206. perform a 1D spline interpolation.
  207. INPUT PARAMETERS
  208. ----------------
  209. x (1D np array) : the original abscisse
  210. y (1D np array) : the original data
  211. nb_pts (integer): the number of points for the interpolation
  212. OUTPUT PARAMETERS
  213. -----------------
  214. xnew (1D np array) : the spline abscisse
  215. ynew (1D np array) : the spline approximations
  216. """
  217. from scipy import interpolate
  218. tck = interpolate.splrep(x, y, s=0)
  219. xnew = np.linspace(x[0], x[x.shape[0] - 1], nb_pts)
  220. ynew = interpolate.splev(xnew, tck, der=0)
  221. return xnew, ynew