pp.x: fix gaussian cube file parsing.

aiida_quantumespresso/parsers/pp.py

@@ -217,7 +217,7 @@ def detect_important_message(logs, line):
  split_line = line.split('=')
  if 'negative/imaginary' in line: # QE6.1
  output_dict['negative_core_charge'] = float(split_line[-1].split()[0])
- output_dict['imaginary_core_charge'] = float(split_line.split()[-1])
+ output_dict['imaginary_core_charge'] = float(split_line[-1].split()[-1])
  else: # QE6.4
  output_dict['negative_core_charge'] = float(split_line[1])
  if 'Min, Max, imaginary charge:' in line:
@@ -335,10 +335,15 @@ def parse_gaussian(self, data_file_str):
  """
  lines = data_file_str.splitlines()
 
+ title = lines[0]
+ comment = lines[1]
+
  atoms_line = lines[2].split()
- atoms = int(atoms_line[0]) # The number of atoms listed in the file
- header = lines[:6 + atoms] # Header of the file: comments, the voxel, and the number of atoms and datapoints
- data_lines = lines[6 + atoms:] # The actual data: atoms and volumetric data
+ natoms = int(atoms_line[0]) # The number of atoms listed in the file
+ origin = np.array(atoms_line[1:], dtype=float)
+
+ header = lines[:6 + natoms] # Header of the file: comments, the voxel, and the number of atoms and datapoints
+ data_lines = lines[6 + natoms:] # The actual data: atoms and volumetric data
 
  # Parse the declared dimensions of the volumetric data
  x_line = header[3].split()
@@ -352,34 +357,31 @@ def parse_gaussian(self, data_file_str):
  voxel_array = np.array([[x_line[1], x_line[2], x_line[3]], [y_line[1], y_line[2], y_line[3]],
  [z_line[1], z_line[2], z_line[3]]],
  dtype=np.float64)
+ atomic_numbers = np.empty(natoms, int)
+ coordinates = np.empty((natoms, 3))
+ for i in range(natoms):
+ line = header[6 + i].split()
+ atomic_numbers[i] = int(line[0])
+ coordinates[i] = [float(s) for s in line[2:]]
 
  # Get the volumetric data
- data_array = np.zeros((xdim, ydim, zdim))
-
- # The data is organised into columns with a fixed number of values
- # Gather up all the data points into one list for easier unpacking
- datalist = []
+ data_array = np.empty(xdim * ydim * zdim, dtype=float)
+ cursor = 0
  for line in data_lines:
- for i in range(0, len(line), 13):
- data_point = line[i:i + 13].strip()
- if data_point != '':
- datalist.append(float(data_point))
-
- # Unpack the list and repack as a 3D array
- # Note the unusual indexing: cube files run over the z index first, then y and x.
- # E.g. The first volumetric data point is x,y,z = (0,0,0) and the second is (0,0,1)
- for i in range(0, xdim):
- for j in range(0, ydim):
- for k in range(0, zdim):
- data_array[i, j, k] = (datalist[(i * ydim * zdim) + (j * zdim) + k])
+ ls = line.split()
+ data_array[cursor:cursor + len(ls)] = ls
+ cursor += len(ls)
+ data_array = data_array.reshape((xdim, ydim, zdim))
 
  coordinates_units = 'bohr'
  data_units = self.units_dict[self.output_parameters['plot_num']]
 
  arraydata = orm.ArrayData()
  arraydata.set_array('voxel', voxel_array)
  arraydata.set_array('data', data_array)
- arraydata.set_array('coordinates_units', np.array(coordinates_units))
  arraydata.set_array('data_units', np.array(data_units))
+ arraydata.set_array('coordinates_units', np.array(coordinates_units))
+ arraydata.set_array('coordinates', coordinates)
+ arraydata.set_array('atomic_numbers', atomic_numbers)
 
  return arraydata

tests/parsers/fixtures/pp/default_3d_multiple/aiida.filplot_K001_B001aiida.fileout

@@ -5,4 +5,4 @@
  1 0.141231 0.000000 0.141231
  1 0.141231 0.141231 0.000000
  14 14.000000 10.168616 10.168616 10.168616
- 0.16567E-01 0.14498E-01 0.96968E-02
+ 0.16567E-01

tests/parsers/fixtures/pp/default_3d_multiple/aiida.filplot_K001_B002aiida.fileout

@@ -5,4 +5,4 @@
  1 0.141231 0.000000 0.141231
  1 0.141231 0.141231 0.000000
  14 14.000000 2.542154 2.542154 2.542154
- 0.22404E-01 0.19330E-01 0.12458E-01
+ 0.22404E-01

tests/parsers/test_pp.py

@@ -271,7 +271,7 @@ def test_pp_default_3d(
  assert calcfunction.is_finished_ok, calcfunction.exit_message
  assert 'output_parameters' in results
  assert 'output_data' in results
- assert len(results['output_data'].get_arraynames()) == 4
+ assert len(results['output_data'].get_arraynames()) == 6
  data_array = results['output_data'].get_array('data').flatten()
  voxel_array = results['output_data'].get_array('voxel').flatten()
  data_units_array = results['output_data'].get_array('data_units')
@@ -313,7 +313,7 @@ def test_pp_default_3d_keep_plot_file(
  assert calcfunction.is_finished_ok, calcfunction.exit_message
  assert 'output_parameters' in results
  assert 'output_data' in results
- assert len(results['output_data'].get_arraynames()) == 4
+ assert len(results['output_data'].get_arraynames()) == 6
 
 
 def test_pp_default_3d_multiple(aiida_profile, generate_calc_job_node, generate_parser, generate_inputs_3d):
@@ -334,7 +334,7 @@ def test_pp_default_3d_multiple(aiida_profile, generate_calc_job_node, generate_
  for key in ['K001_B001', 'K001_B002']:
  assert key in results['output_data_multiple']
  node = results['output_data_multiple'][key]
- assert len(node.get_arraynames()) == 4
+ assert len(node.get_arraynames()) == 6
 
 
 def test_pp_default_3d_failed_missing(