Coverage for /builds/ase/ase/ase/calculators/turbomole/reader.py : 9.73%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""Functions to read from control file and from turbomole standard output"""
3import os
4import re
5import warnings
6import subprocess
7import numpy as np
8from ase import Atom, Atoms
9from ase.units import Ha, Bohr
10from ase.calculators.calculator import ReadError
13def execute_command(args):
14 """execute commands like sdg, eiger"""
15 proc = subprocess.Popen(args, stdout=subprocess.PIPE, encoding='ASCII')
16 stdout, stderr = proc.communicate()
17 return stdout
20def read_data_group(data_group):
21 """read a turbomole data group from control file"""
22 return execute_command(['sdg', data_group]).strip()
25def parse_data_group(dgr, dg_name):
26 """parse a data group"""
27 if len(dgr) == 0:
28 return None
29 dg_key = '$' + dg_name
30 if not dgr.startswith(dg_key):
31 raise ValueError(f'data group does not start with {dg_key}')
32 ndg = dgr.replace(dg_key, '').strip()
33 ndg = re.sub(r'=\s+', '=', re.sub(r'\s+=', '=', ndg))
34 if all(c not in ndg for c in ('\n', ' ', '=')):
35 return ndg
36 lsep = '\n' if '\n' in dgr else ' '
37 result = {}
38 lines = ndg.split(lsep)
39 for line in lines:
40 if len(line) == 0:
41 continue
42 ksep = '=' if '=' in line else None
43 fields = line.strip().split(ksep)
44 if len(fields) == 2:
45 result[fields[0]] = fields[1]
46 elif len(fields) == 1:
47 result[fields[0]] = True
48 return result
51def read_output(regex, path):
52 """collects all matching strings from the output"""
53 hitlist = []
54 checkfiles = []
55 for filename in os.listdir(path):
56 if filename.startswith('job.') or filename.endswith('.out'):
57 checkfiles.append(filename)
58 for filename in checkfiles:
59 with open(filename, 'rt') as f:
60 lines = f.readlines()
61 for line in lines:
62 match = re.search(regex, line)
63 if match:
64 hitlist.append(match.group(1))
65 return hitlist
68def read_version(path):
69 """read the version from the tm output if stored in a file"""
70 versions = read_output(r'TURBOMOLE\s+V(\d+\.\d+)\s+', path)
71 if len(set(versions)) > 1:
72 warnings.warn('different turbomole versions detected')
73 version = list(set(versions))
74 elif len(versions) == 0:
75 warnings.warn('no turbomole version detected')
76 version = None
77 else:
78 version = versions[0]
79 return version
82def read_datetime(path):
83 """read the datetime of the most recent calculation
84 from the tm output if stored in a file
85 """
86 datetimes = read_output(
87 r'(\d{4}-[01]\d-[0-3]\d([T\s][0-2]\d:[0-5]'
88 r'\d:[0-5]\d\.\d+)?([+-][0-2]\d:[0-5]\d|Z)?)', path)
89 if len(datetimes) == 0:
90 warnings.warn('no turbomole datetime detected')
91 datetime = None
92 else:
93 # take the most recent time stamp
94 datetime = sorted(datetimes, reverse=True)[0]
95 return datetime
98def read_runtime(path):
99 """read the total runtime of calculations"""
100 hits = read_output(r'total wall-time\s+:\s+(\d+.\d+)\s+seconds', path)
101 if len(hits) == 0:
102 warnings.warn('no turbomole runtimes detected')
103 runtime = None
104 else:
105 runtime = np.sum([float(a) for a in hits])
106 return runtime
109def read_hostname(path):
110 """read the hostname of the computer on which the calc has run"""
111 hostnames = read_output(r'hostname is\s+(.+)', path)
112 if len(set(hostnames)) > 1:
113 warnings.warn('runs on different hosts detected')
114 hostname = list(set(hostnames))
115 else:
116 hostname = hostnames[0]
117 return hostname
120def read_convergence(restart, parameters):
121 """perform convergence checks"""
122 if restart:
123 if bool(len(read_data_group('restart'))):
124 return False
125 if bool(len(read_data_group('actual'))):
126 return False
127 if not bool(len(read_data_group('energy'))):
128 return False
129 if (os.path.exists('job.start') and
130 os.path.exists('GEO_OPT_FAILED')):
131 return False
132 return True
134 if parameters['task'] in ['optimize', 'geometry optimization']:
135 if os.path.exists('GEO_OPT_CONVERGED'):
136 return True
137 elif os.path.exists('GEO_OPT_FAILED'):
138 # check whether a failed scf convergence is the reason
139 checkfiles = []
140 for filename in os.listdir('.'):
141 if filename.startswith('job.'):
142 checkfiles.append(filename)
143 for filename in checkfiles:
144 for line in open(filename):
145 if 'SCF FAILED TO CONVERGE' in line:
146 # scf did not converge in some jobex iteration
147 if filename == 'job.last':
148 raise RuntimeError('scf failed to converge')
149 else:
150 warnings.warn('scf failed to converge')
151 warnings.warn('geometry optimization failed to converge')
152 return False
153 else:
154 raise RuntimeError('error during geometry optimization')
155 else:
156 if os.path.isfile('dscf_problem'):
157 raise RuntimeError('scf failed to converge')
158 else:
159 return True
162def read_run_parameters(results):
163 """read parameters set by define and not in self.parameters"""
165 if 'calculation parameters' not in results.keys():
166 results['calculation parameters'] = {}
167 parameters = results['calculation parameters']
168 dg = read_data_group('symmetry')
169 parameters['point group'] = str(dg.split()[1])
170 parameters['uhf'] = '$uhf' in read_data_group('uhf')
171 # Gaussian function type
172 gt = read_data_group('pople')
173 if gt == '':
174 parameters['gaussian type'] = 'spherical harmonic'
175 else:
176 gt = gt.split()[1]
177 if gt == 'AO':
178 parameters['gaussian type'] = 'spherical harmonic'
179 elif gt == 'CAO':
180 parameters['gaussian type'] = 'cartesian'
181 else:
182 parameters['gaussian type'] = None
184 nvibro = read_data_group('nvibro')
185 if nvibro:
186 parameters['nuclear degrees of freedom'] = int(nvibro.split()[1])
189def read_energy(results, post_HF):
190 """Read energy from Turbomole energy file."""
191 try:
192 with open('energy', 'r') as enf:
193 text = enf.read().lower()
194 except IOError:
195 raise ReadError('failed to read energy file')
196 if text == '':
197 raise ReadError('empty energy file')
199 lines = iter(text.split('\n'))
201 for line in lines:
202 if line.startswith('$end'):
203 break
204 elif line.startswith('$'):
205 pass
206 else:
207 energy_tmp = float(line.split()[1])
208 if post_HF:
209 energy_tmp += float(line.split()[4])
210 # update energy units
211 e_total = energy_tmp * Ha
212 results['total energy'] = e_total
215def read_occupation_numbers(results):
216 """read occupation numbers with module 'eiger' """
217 if 'molecular orbitals' not in results.keys():
218 return
219 mos = results['molecular orbitals']
220 lines = execute_command(['eiger', '--all', '--pview']).split('\n')
221 for line in lines:
222 regex = (
223 r'^\s+(\d+)\.\s([\sab])\s*(\d+)\s?(\w+)'
224 r'\s+(\d*\.*\d*)\s+([-+]?\d+\.\d*)'
225 )
226 match = re.search(regex, line)
227 if match:
228 orb_index = int(match.group(3))
229 if match.group(2) == 'a':
230 spin = 'alpha'
231 elif match.group(2) == 'b':
232 spin = 'beta'
233 else:
234 spin = None
235 ar_index = next(
236 index for (index, molecular_orbital) in enumerate(mos)
237 if (molecular_orbital['index'] == orb_index and
238 molecular_orbital['spin'] == spin)
239 )
240 mos[ar_index]['index by energy'] = int(match.group(1))
241 irrep = str(match.group(4))
242 mos[ar_index]['irreducible representation'] = irrep
243 if match.group(5) != '':
244 mos[ar_index]['occupancy'] = float(match.group(5))
245 else:
246 mos[ar_index]['occupancy'] = float(0)
249def read_mos(results):
250 """read the molecular orbital coefficients and orbital energies
251 from files mos, alpha and beta"""
253 results['molecular orbitals'] = []
254 mos = results['molecular orbitals']
255 keywords = ['scfmo', 'uhfmo_alpha', 'uhfmo_beta']
256 spin = [None, 'alpha', 'beta']
257 converged = None
259 for index, keyword in enumerate(keywords):
260 flen = None
261 mo = {}
262 orbitals_coefficients_line = []
263 mo_string = read_data_group(keyword)
264 if mo_string == '':
265 continue
266 mo_string += '\n$end'
267 lines = mo_string.split('\n')
268 for line in lines:
269 if re.match(r'^\s*#', line):
270 continue
271 if 'eigenvalue' in line:
272 if len(orbitals_coefficients_line) != 0:
273 mo['eigenvector'] = orbitals_coefficients_line
274 mos.append(mo)
275 mo = {}
276 orbitals_coefficients_line = []
277 regex = (r'^\s*(\d+)\s+(\S+)\s+'
278 r'eigenvalue=([\+\-\d\.\w]+)\s')
279 match = re.search(regex, line)
280 mo['index'] = int(match.group(1))
281 mo['irreducible representation'] = str(match.group(2))
282 eig = float(re.sub('[dD]', 'E', match.group(3))) * Ha
283 mo['eigenvalue'] = eig
284 mo['spin'] = spin[index]
285 mo['degeneracy'] = 1
286 continue
287 if keyword in line:
288 # e.g. format(4d20.14)
289 regex = r'format\(\d+[a-zA-Z](\d+)\.\d+\)'
290 match = re.search(regex, line)
291 if match:
292 flen = int(match.group(1))
293 if ('scfdump' in line or 'expanded' in line or
294 'scfconv' not in line):
295 converged = False
296 continue
297 if '$end' in line:
298 if len(orbitals_coefficients_line) != 0:
299 mo['eigenvector'] = orbitals_coefficients_line
300 mos.append(mo)
301 break
302 sfields = [line[i:i + flen]
303 for i in range(0, len(line), flen)]
304 ffields = [float(f.replace('D', 'E').replace('d', 'E'))
305 for f in sfields]
306 orbitals_coefficients_line += ffields
307 return converged
310def read_basis_set(results):
311 """read the basis set"""
312 results['basis set'] = []
313 results['basis set formatted'] = {}
314 bsf = read_data_group('basis')
315 results['basis set formatted']['turbomole'] = bsf
316 lines = bsf.split('\n')
317 basis_set = {}
318 functions = []
319 function = {}
320 primitives = []
321 read_tag = False
322 read_data = False
323 for line in lines:
324 if len(line.strip()) == 0:
325 continue
326 if '$basis' in line:
327 continue
328 if '$end' in line:
329 break
330 if re.match(r'^\s*#', line):
331 continue
332 if re.match(r'^\s*\*', line):
333 if read_tag:
334 read_tag = False
335 read_data = True
336 else:
337 if read_data:
338 # end primitives
339 function['primitive functions'] = primitives
340 function['number of primitives'] = len(primitives)
341 primitives = []
342 functions.append(function)
343 function = {}
344 # end contracted
345 basis_set['functions'] = functions
346 functions = []
347 results['basis set'].append(basis_set)
348 basis_set = {}
349 read_data = False
350 read_tag = True
351 continue
352 if read_tag:
353 match = re.search(r'^\s*(\w+)\s+(.+)', line)
354 if match:
355 basis_set['element'] = match.group(1)
356 basis_set['nickname'] = match.group(2)
357 else:
358 raise RuntimeError('error reading basis set')
359 else:
360 match = re.search(r'^\s+(\d+)\s+(\w+)', line)
361 if match:
362 if len(primitives) > 0:
363 # end primitives
364 function['primitive functions'] = primitives
365 function['number of primitives'] = len(primitives)
366 primitives = []
367 functions.append(function)
368 function = {}
369 # begin contracted
370 function['shell type'] = str(match.group(2))
371 continue
372 regex = (
373 r'^\s*([-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)'
374 r'\s+([-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)'
375 )
376 match = re.search(regex, line)
377 if match:
378 exponent = float(match.group(1))
379 coefficient = float(match.group(3))
380 primitives.append(
381 {'exponent': exponent, 'coefficient': coefficient}
382 )
385def read_ecps(results):
386 """read the effective core potentials"""
387 ecpf = read_data_group('ecp')
388 if not bool(len(ecpf)):
389 results['ecps'] = None
390 results['ecps formatted'] = None
391 return
392 results['ecps'] = []
393 results['ecps formatted'] = {}
394 results['ecps formatted']['turbomole'] = ecpf
395 lines = ecpf.split('\n')
396 ecp = {}
397 groups = []
398 group = {}
399 terms = []
400 read_tag = False
401 read_data = False
402 for line in lines:
403 if len(line.strip()) == 0:
404 continue
405 if '$ecp' in line:
406 continue
407 if '$end' in line:
408 break
409 if re.match(r'^\s*#', line):
410 continue
411 if re.match(r'^\s*\*', line):
412 if read_tag:
413 read_tag = False
414 read_data = True
415 else:
416 if read_data:
417 # end terms
418 group['terms'] = terms
419 group['number of terms'] = len(terms)
420 terms = []
421 groups.append(group)
422 group = {}
423 # end group
424 ecp['groups'] = groups
425 groups = []
426 results['ecps'].append(ecp)
427 ecp = {}
428 read_data = False
429 read_tag = True
430 continue
431 if read_tag:
432 match = re.search(r'^\s*(\w+)\s+(.+)', line)
433 if match:
434 ecp['element'] = match.group(1)
435 ecp['nickname'] = match.group(2)
436 else:
437 raise RuntimeError('error reading ecp')
438 else:
439 regex = r'ncore\s*=\s*(\d+)\s+lmax\s*=\s*(\d+)'
440 match = re.search(regex, line)
441 if match:
442 ecp['number of core electrons'] = int(match.group(1))
443 ecp['maximum angular momentum number'] = \
444 int(match.group(2))
445 continue
446 match = re.search(r'^(\w(\-\w)?)', line)
447 if match:
448 if len(terms) > 0:
449 # end terms
450 group['terms'] = terms
451 group['number of terms'] = len(terms)
452 terms = []
453 groups.append(group)
454 group = {}
455 # begin group
456 group['title'] = str(match.group(1))
457 continue
458 regex = (r'^\s*([-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)\s+'
459 r'(\d)\s+([-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)')
460 match = re.search(regex, line)
461 if match:
462 terms.append(
463 {
464 'coefficient': float(match.group(1)),
465 'power of r': float(match.group(3)),
466 'exponent': float(match.group(4))
467 }
468 )
471def read_forces(results, natoms):
472 """Read forces from Turbomole gradient file."""
473 dg = read_data_group('grad')
474 if len(dg) == 0:
475 return
476 file = open('gradient', 'r')
477 lines = file.readlines()
478 file.close()
480 forces = np.array([[0, 0, 0]])
482 nline = len(lines)
483 iline = -1
485 for i in range(nline):
486 if 'cycle' in lines[i]:
487 iline = i
489 if iline < 0:
490 raise RuntimeError('Please check TURBOMOLE gradients')
492 # next line
493 iline += natoms + 1
494 # $end line
495 nline -= 1
496 # read gradients
497 for i in range(iline, nline):
498 line = lines[i].replace('D', 'E')
499 tmp = np.array([[float(f) for f in line.split()[0:3]]])
500 forces = np.concatenate((forces, tmp))
501 # Note the '-' sign for turbomole, to get forces
502 forces = -np.delete(forces, np.s_[0:1], axis=0) * Ha / Bohr
503 results['energy gradient'] = (-forces).tolist()
504 return forces
507def read_gradient(results):
508 """read all information in file 'gradient'"""
509 grad_string = read_data_group('grad')
510 if len(grad_string) == 0:
511 return
512# try to reuse ase:
513# structures = read('gradient', index=':')
514 lines = grad_string.split('\n')
515 history = []
516 image = {}
517 gradient = []
518 atoms = Atoms()
519 (cycle, energy, norm) = (None, None, None)
520 for line in lines:
521 # cycle lines
522 regex = (
523 r'^\s*cycle =\s*(\d+)\s+'
524 r'SCF energy =\s*([-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)\s+'
525 r'\|dE\/dxyz\| =\s*([-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)'
526 )
527 match = re.search(regex, line)
528 if match:
529 if len(atoms):
530 image['optimization cycle'] = cycle
531 image['total energy'] = energy
532 image['gradient norm'] = norm
533 image['energy gradient'] = gradient
534 history.append(image)
535 image = {}
536 atoms = Atoms()
537 gradient = []
538 cycle = int(match.group(1))
539 energy = float(match.group(2)) * Ha
540 norm = float(match.group(4)) * Ha / Bohr
541 continue
542 # coordinate lines
543 regex = (
544 r'^\s*([-+]?[0-9]*\.?[0-9]+([eEdD][-+]?[0-9]+)?)'
545 r'\s+([-+]?[0-9]*\.?[0-9]+([eEdD][-+]?[0-9]+)?)'
546 r'\s+([-+]?[0-9]*\.?[0-9]+([eEdD][-+]?[0-9]+)?)'
547 r'\s+(\w+)'
548 )
549 match = re.search(regex, line)
550 if match:
551 x = float(match.group(1)) * Bohr
552 y = float(match.group(3)) * Bohr
553 z = float(match.group(5)) * Bohr
554 symbol = str(match.group(7)).capitalize()
556 if symbol == 'Q':
557 symbol = 'X'
558 atoms += Atom(symbol, (x, y, z))
560 continue
561 # gradient lines
562 regex = (
563 r'^\s*([-+]?[0-9]*\.?[0-9]+([eEdD][-+]?[0-9]+)?)'
564 r'\s+([-+]?[0-9]*\.?[0-9]+([eEdD][-+]?[0-9]+)?)'
565 r'\s+([-+]?[0-9]*\.?[0-9]+([eEdD][-+]?[0-9]+)?)'
566 )
567 match = re.search(regex, line)
568 if match:
569 gradx = float(match.group(1).replace('D', 'E')) * Ha / Bohr
570 grady = float(match.group(3).replace('D', 'E')) * Ha / Bohr
571 gradz = float(match.group(5).replace('D', 'E')) * Ha / Bohr
572 gradient.append([gradx, grady, gradz])
574 image['optimization cycle'] = cycle
575 image['total energy'] = energy
576 image['gradient norm'] = norm
577 image['energy gradient'] = gradient
578 history.append(image)
579 results['geometry optimization history'] = history
582def read_hessian(results, noproj=False):
583 """Read in the hessian matrix"""
584 results['hessian matrix'] = {}
585 results['hessian matrix']['array'] = []
586 results['hessian matrix']['units'] = '?'
587 results['hessian matrix']['projected'] = True
588 results['hessian matrix']['mass weighted'] = True
589 dg = read_data_group('nvibro')
590 if len(dg) == 0:
591 return
592 nvibro = int(dg.split()[1])
593 results['hessian matrix']['dimension'] = nvibro
594 row = []
595 key = 'hessian'
596 if noproj:
597 key = 'npr' + key
598 results['hessian matrix']['projected'] = False
599 lines = read_data_group(key).split('\n')
600 for line in lines:
601 if key in line:
602 continue
603 fields = line.split()
604 row.extend(fields[2:len(fields)])
605 if len(row) == nvibro:
606 # check whether it is mass-weighted
607 float_row = [float(element) for element in row]
608 results['hessian matrix']['array'].append(float_row)
609 row = []
612def read_normal_modes(results, noproj=False):
613 """Read in vibrational normal modes"""
614 results['normal modes'] = {}
615 results['normal modes']['array'] = []
616 results['normal modes']['projected'] = True
617 results['normal modes']['mass weighted'] = True
618 results['normal modes']['units'] = '?'
619 dg = read_data_group('nvibro')
620 if len(dg) == 0:
621 return
622 nvibro = int(dg.split()[1])
623 results['normal modes']['dimension'] = nvibro
624 row = []
625 key = 'vibrational normal modes'
626 if noproj:
627 key = 'npr' + key
628 results['normal modes']['projected'] = False
629 lines = read_data_group(key).split('\n')
630 for line in lines:
631 if key in line:
632 continue
633 if '$end' in line:
634 break
635 fields = line.split()
636 row.extend(fields[2:len(fields)])
637 if len(row) == nvibro:
638 # check whether it is mass-weighted
639 float_row = [float(element) for element in row]
640 results['normal modes']['array'].append(float_row)
641 row = []
644def read_vibrational_reduced_masses(results):
645 """Read vibrational reduced masses"""
646 results['vibrational reduced masses'] = []
647 dg = read_data_group('vibrational reduced masses')
648 if len(dg) == 0:
649 return
650 lines = dg.split('\n')
651 for line in lines:
652 if '$vibrational' in line:
653 continue
654 if '$end' in line:
655 break
656 fields = [float(element) for element in line.split()]
657 results['vibrational reduced masses'].extend(fields)
660def read_vibrational_spectrum(results, noproj=False):
661 """Read the vibrational spectrum"""
662 results['vibrational spectrum'] = []
663 key = 'vibrational spectrum'
664 if noproj:
665 key = 'npr' + key
666 lines = read_data_group(key).split('\n')
667 for line in lines:
668 dictionary = {}
669 regex = (
670 r'^\s+(\d+)\s+(\S*)\s+([-+]?\d+\.\d*)'
671 r'\s+(\d+\.\d*)\s+(\S+)\s+(\S+)'
672 )
673 match = re.search(regex, line)
674 if match:
675 dictionary['mode number'] = int(match.group(1))
676 dictionary['irreducible representation'] = str(match.group(2))
677 dictionary['frequency'] = {
678 'units': 'cm^-1',
679 'value': float(match.group(3))
680 }
681 dictionary['infrared intensity'] = {
682 'units': 'km/mol',
683 'value': float(match.group(4))
684 }
686 if match.group(5) == 'YES':
687 dictionary['infrared active'] = True
688 elif match.group(5) == 'NO':
689 dictionary['infrared active'] = False
690 else:
691 dictionary['infrared active'] = None
693 if match.group(6) == 'YES':
694 dictionary['Raman active'] = True
695 elif match.group(6) == 'NO':
696 dictionary['Raman active'] = False
697 else:
698 dictionary['Raman active'] = None
700 results['vibrational spectrum'].append(dictionary)
703def read_ssquare(results):
704 """Read the expectation value of S^2 operator"""
705 s2_string = read_data_group('ssquare from dscf')
706 if s2_string == '':
707 return
708 string = s2_string.split('\n')[1]
709 ssquare = float(re.search(r'^\s*(\d+\.*\d*)', string).group(1))
710 results['ssquare from scf calculation'] = ssquare
713def read_dipole_moment(results):
714 """Read the dipole moment"""
715 dip_string = read_data_group('dipole')
716 if dip_string == '':
717 return
718 lines = dip_string.split('\n')
719 for line in lines:
720 regex = (
721 r'^\s+x\s+([-+]?\d+\.\d*)\s+y\s+([-+]?\d+\.\d*)'
722 r'\s+z\s+([-+]?\d+\.\d*)\s+a\.u\.'
723 )
724 match = re.search(regex, line)
725 if match:
726 dip_vec = [float(match.group(c)) for c in range(1, 4)]
727 regex = r'^\s+\| dipole \| =\s+(\d+\.*\d*)\s+debye'
728 match = re.search(regex, line)
729 if match:
730 dip_abs_val = float(match.group(1))
731 results['electric dipole moment'] = {}
732 results['electric dipole moment']['vector'] = {
733 'array': dip_vec,
734 'units': 'a.u.'
735 }
736 results['electric dipole moment']['absolute value'] = {
737 'value': dip_abs_val,
738 'units': 'Debye'
739 }
742def read_charges(filename, natoms):
743 """read partial charges on atoms from an ESP fit"""
744 charges = None
745 if os.path.exists(filename):
746 with open(filename, 'r') as infile:
747 lines = infile.readlines()
748 oklines = None
749 for n, line in enumerate(lines):
750 if 'atom radius/au charge' in line:
751 oklines = lines[n + 1:n + natoms + 1]
752 if oklines is not None:
753 qm_charges = [float(line.split()[3]) for line in oklines]
754 charges = np.array(qm_charges)
755 return charges
758def read_point_charges():
759 """read point charges from previous calculation"""
760 pcs = read_data_group('point_charges')
761 lines = pcs.split('\n')[1:]
762 (charges, positions) = ([], [])
763 for line in lines:
764 columns = [float(col) for col in line.strip().split()]
765 positions.append([col * Bohr for col in columns[0:3]])
766 charges.append(columns[3])
767 return charges, positions