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

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