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

1"""Functions to read from control file and from turbomole standard output"""

3import os

4import re

5import subprocess

6import warnings

8import numpy as np

10from ase import Atom, Atoms

11from ase.calculators.calculator import ReadError

12from ase.units import Bohr, Ha

15def execute_command(args):

16 """execute commands like sdg, eiger"""

17 proc = subprocess.Popen(args, stdout=subprocess.PIPE, encoding='ASCII')

18 stdout, stderr = proc.communicate()

19 return stdout

22def read_data_group(data_group):

23 """read a turbomole data group from control file"""

24 return execute_command(['sdg', data_group]).strip()

27def parse_data_group(dgr, dg_name):

28 """parse a data group"""

29 if len(dgr) == 0:

30 return None

31 dg_key = '$' + dg_name

32 if not dgr.startswith(dg_key):

33 raise ValueError(f'data group does not start with {dg_key}')

34 ndg = dgr.replace(dg_key, '').strip()

35 ndg = re.sub(r'=\s+', '=', re.sub(r'\s+=', '=', ndg))

36 if all(c not in ndg for c in ('\n', ' ', '=')):

37 return ndg

38 lsep = '\n' if '\n' in dgr else ' '

39 result = {}

40 lines = ndg.split(lsep)

41 for line in lines:

42 if len(line) == 0:

43 continue

44 ksep = '=' if '=' in line else None

45 fields = line.strip().split(ksep)

46 if len(fields) == 2:

47 result[fields[0]] = fields[1]

48 elif len(fields) == 1:

49 result[fields[0]] = True

50 return result

53def read_output(regex, path):

54 """collects all matching strings from the output"""

55 hitlist = []

56 checkfiles = []

57 for filename in os.listdir(path):

58 if filename.startswith('job.') or filename.endswith('.out'):

59 checkfiles.append(filename)

60 for filename in checkfiles:

61 with open(filename) as f:

62 lines = f.readlines()

63 for line in lines:

64 match = re.search(regex, line)

65 if match:

66 hitlist.append(match.group(1))

67 return hitlist

70def read_version(path):

71 """read the version from the tm output if stored in a file"""

72 versions = read_output(r'TURBOMOLE\s+V(\d+\.\d+)\s+', path)

73 if len(set(versions)) > 1:

74 warnings.warn('different turbomole versions detected')

75 version = list(set(versions))

76 elif len(versions) == 0:

77 warnings.warn('no turbomole version detected')

78 version = None

79 else:

80 version = versions[0]

81 return version

84def read_datetime(path):

85 """read the datetime of the most recent calculation

86 from the tm output if stored in a file

87 """

88 datetimes = read_output(

89 r'(\d{4}-[01]\d-[0-3]\d([T\s][0-2]\d:[0-5]'

90 r'\d:[0-5]\d\.\d+)?([+-][0-2]\d:[0-5]\d|Z)?)', path)

91 if len(datetimes) == 0:

92 warnings.warn('no turbomole datetime detected')

93 datetime = None

94 else:

95 # take the most recent time stamp

96 datetime = sorted(datetimes, reverse=True)[0]

97 return datetime

100def read_runtime(path):

101 """read the total runtime of calculations"""

102 hits = read_output(r'total wall-time\s+:\s+(\d+.\d+)\s+seconds', path)

103 if len(hits) == 0:

104 warnings.warn('no turbomole runtimes detected')

105 runtime = None

106 else:

107 runtime = np.sum([float(a) for a in hits])

108 return runtime

109

110

111def read_hostname(path):

112 """read the hostname of the computer on which the calc has run"""

113 hostnames = read_output(r'hostname is\s+(.+)', path)

114 if len(set(hostnames)) > 1:

115 warnings.warn('runs on different hosts detected')

116 hostname = list(set(hostnames))

117 else:

118 hostname = hostnames[0]

119 return hostname

120

121

122def read_convergence(restart, parameters):

123 """perform convergence checks"""

124 if restart:

125 if bool(len(read_data_group('restart'))):

126 return False

127 if bool(len(read_data_group('actual'))):

128 return False

129 if not bool(len(read_data_group('energy'))):

130 return False

131 if (os.path.exists('job.start') and

132 os.path.exists('GEO_OPT_FAILED')):

133 return False

134 return True

135

136 if parameters['task'] in ['optimize', 'geometry optimization']:

137 if os.path.exists('GEO_OPT_CONVERGED'):

138 return True

139 elif os.path.exists('GEO_OPT_FAILED'):

140 # check whether a failed scf convergence is the reason

141 checkfiles = []

142 for filename in os.listdir('.'):

143 if filename.startswith('job.'):

144 checkfiles.append(filename)

145 for filename in checkfiles:

146 for line in open(filename):

147 if 'SCF FAILED TO CONVERGE' in line:

148 # scf did not converge in some jobex iteration

149 if filename == 'job.last':

150 raise RuntimeError('scf failed to converge')

151 else:

152 warnings.warn('scf failed to converge')

153 warnings.warn('geometry optimization failed to converge')

154 return False

155 else:

156 raise RuntimeError('error during geometry optimization')

157 else:

158 if os.path.isfile('dscf_problem'):

159 raise RuntimeError('scf failed to converge')

160 else:

161 return True

162

163

164def read_run_parameters(results):

165 """read parameters set by define and not in self.parameters"""

166

167 if 'calculation parameters' not in results.keys():

168 results['calculation parameters'] = {}

169 parameters = results['calculation parameters']

170 dg = read_data_group('symmetry')

171 parameters['point group'] = str(dg.split()[1])

172 parameters['uhf'] = '$uhf' in read_data_group('uhf')

173 # Gaussian function type

174 gt = read_data_group('pople')

175 if gt == '':

176 parameters['gaussian type'] = 'spherical harmonic'

177 else:

178 gt = gt.split()[1]

179 if gt == 'AO':

180 parameters['gaussian type'] = 'spherical harmonic'

181 elif gt == 'CAO':

182 parameters['gaussian type'] = 'cartesian'

183 else:

184 parameters['gaussian type'] = None

185

186 nvibro = read_data_group('nvibro')

187 if nvibro:

188 parameters['nuclear degrees of freedom'] = int(nvibro.split()[1])

189

190

191def read_energy(results, post_HF):

192 """Read energy from Turbomole energy file."""

193 try:

194 with open('energy') as enf:

195 text = enf.read().lower()

196 except OSError:

197 raise ReadError('failed to read energy file')

198 if text == '':

199 raise ReadError('empty energy file')

200

201 lines = iter(text.split('\n'))

202

203 for line in lines:

204 if line.startswith('$end'):

205 break

206 elif line.startswith('$'):

207 pass

208 else:

209 energy_tmp = float(line.split()[1])

210 if post_HF:

211 energy_tmp += float(line.split()[4])

212 # update energy units

213 e_total = energy_tmp * Ha

214 results['total energy'] = e_total

215

216

217def read_occupation_numbers(results):

218 """read occupation numbers with module 'eiger' """

219 if 'molecular orbitals' not in results.keys():

220 return

221 mos = results['molecular orbitals']

222 lines = execute_command(['eiger', '--all', '--pview']).split('\n')

223 for line in lines:

224 regex = (

225 r'^\s+(\d+)\.\s([\sab])\s*(\d+)\s?(\w+)'

226 r'\s+(\d*\.*\d*)\s+([-+]?\d+\.\d*)'

227 )

228 match = re.search(regex, line)

229 if match:

230 orb_index = int(match.group(3))

231 if match.group(2) == 'a':

232 spin = 'alpha'

233 elif match.group(2) == 'b':

234 spin = 'beta'

235 else:

236 spin = None

237 ar_index = next(

238 index for (index, molecular_orbital) in enumerate(mos)

239 if (molecular_orbital['index'] == orb_index and

240 molecular_orbital['spin'] == spin)

241 )

242 mos[ar_index]['index by energy'] = int(match.group(1))

243 irrep = str(match.group(4))

244 mos[ar_index]['irreducible representation'] = irrep

245 if match.group(5) != '':

246 mos[ar_index]['occupancy'] = float(match.group(5))

247 else:

248 mos[ar_index]['occupancy'] = float(0)

249

250

251def read_mos(results):

252 """read the molecular orbital coefficients and orbital energies

253 from files mos, alpha and beta"""

254

255 results['molecular orbitals'] = []

256 mos = results['molecular orbitals']

257 keywords = ['scfmo', 'uhfmo_alpha', 'uhfmo_beta']

258 spin = [None, 'alpha', 'beta']

259 converged = None

260

261 for index, keyword in enumerate(keywords):

262 flen = None

263 mo = {}

264 orbitals_coefficients_line = []

265 mo_string = read_data_group(keyword)

266 if mo_string == '':

267 continue

268 mo_string += '\n$end'

269 lines = mo_string.split('\n')

270 for line in lines:

271 if re.match(r'^\s*#', line):

272 continue

273 if 'eigenvalue' in line:

274 if len(orbitals_coefficients_line) != 0:

275 mo['eigenvector'] = orbitals_coefficients_line

276 mos.append(mo)

277 mo = {}

278 orbitals_coefficients_line = []

279 regex = (r'^\s*(\d+)\s+(\S+)\s+'

280 r'eigenvalue=([\+\-\d\.\w]+)\s')

281 match = re.search(regex, line)

282 mo['index'] = int(match.group(1))

283 mo['irreducible representation'] = str(match.group(2))

284 eig = float(re.sub('[dD]', 'E', match.group(3))) * Ha

285 mo['eigenvalue'] = eig

286 mo['spin'] = spin[index]

287 mo['degeneracy'] = 1

288 continue

289 if keyword in line:

290 # e.g. format(4d20.14)

291 regex = r'format$\d+[a-zA-Z](\d+)\.\d+$'

292 match = re.search(regex, line)

293 if match:

294 flen = int(match.group(1))

295 if ('scfdump' in line or 'expanded' in line or

296 'scfconv' not in line):

297 converged = False

298 continue

299 if '$end' in line:

300 if len(orbitals_coefficients_line) != 0:

301 mo['eigenvector'] = orbitals_coefficients_line

302 mos.append(mo)

303 break

304 sfields = [line[i:i + flen]

305 for i in range(0, len(line), flen)]

306 ffields = [float(f.replace('D', 'E').replace('d', 'E'))

307 for f in sfields]

308 orbitals_coefficients_line += ffields

309 return converged

310

311

312def read_basis_set(results):

313 """read the basis set"""

314 results['basis set'] = []

315 results['basis set formatted'] = {}

316 bsf = read_data_group('basis')

317 results['basis set formatted']['turbomole'] = bsf

318 lines = bsf.split('\n')

319 basis_set = {}

320 functions = []

321 function = {}

322 primitives = []

323 read_tag = False

324 read_data = False

325 for line in lines:

326 if len(line.strip()) == 0:

327 continue

328 if '$basis' in line:

329 continue

330 if '$end' in line:

331 break

332 if re.match(r'^\s*#', line):

333 continue

334 if re.match(r'^\s*\*', line):

335 if read_tag:

336 read_tag = False

337 read_data = True

338 else:

339 if read_data:

340 # end primitives

341 function['primitive functions'] = primitives

342 function['number of primitives'] = len(primitives)

343 primitives = []

344 functions.append(function)

345 function = {}

346 # end contracted

347 basis_set['functions'] = functions

348 functions = []

349 results['basis set'].append(basis_set)

350 basis_set = {}

351 read_data = False

352 read_tag = True

353 continue

354 if read_tag:

355 match = re.search(r'^\s*(\w+)\s+(.+)', line)

356 if match:

357 basis_set['element'] = match.group(1)

358 basis_set['nickname'] = match.group(2)

359 else:

360 raise RuntimeError('error reading basis set')

361 else:

362 match = re.search(r'^\s+(\d+)\s+(\w+)', line)

363 if match:

364 if len(primitives) > 0:

365 # end primitives

366 function['primitive functions'] = primitives

367 function['number of primitives'] = len(primitives)

368 primitives = []

369 functions.append(function)

370 function = {}

371 # begin contracted

372 function['shell type'] = str(match.group(2))

373 continue

374 regex = (

375 r'^\s*([-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)'

376 r'\s+([-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)'

377 )

378 match = re.search(regex, line)

379 if match:

380 exponent = float(match.group(1))

381 coefficient = float(match.group(3))

382 primitives.append(

383 {'exponent': exponent, 'coefficient': coefficient}

384 )

385

386

387def read_ecps(results):

388 """read the effective core potentials"""

389 ecpf = read_data_group('ecp')

390 if not bool(len(ecpf)):

391 results['ecps'] = None

392 results['ecps formatted'] = None

393 return

394 results['ecps'] = []

395 results['ecps formatted'] = {}

396 results['ecps formatted']['turbomole'] = ecpf

397 lines = ecpf.split('\n')

398 ecp = {}

399 groups = []

400 group = {}

401 terms = []

402 read_tag = False

403 read_data = False

404 for line in lines:

405 if len(line.strip()) == 0:

406 continue

407 if '$ecp' in line:

408 continue

409 if '$end' in line:

410 break

411 if re.match(r'^\s*#', line):

412 continue

413 if re.match(r'^\s*\*', line):

414 if read_tag:

415 read_tag = False

416 read_data = True

417 else:

418 if read_data:

419 # end terms

420 group['terms'] = terms

421 group['number of terms'] = len(terms)

422 terms = []

423 groups.append(group)

424 group = {}

425 # end group

426 ecp['groups'] = groups

427 groups = []

428 results['ecps'].append(ecp)

429 ecp = {}

430 read_data = False

431 read_tag = True

432 continue

433 if read_tag:

434 match = re.search(r'^\s*(\w+)\s+(.+)', line)

435 if match:

436 ecp['element'] = match.group(1)

437 ecp['nickname'] = match.group(2)

438 else:

439 raise RuntimeError('error reading ecp')

440 else:

441 regex = r'ncore\s*=\s*(\d+)\s+lmax\s*=\s*(\d+)'

442 match = re.search(regex, line)

443 if match:

444 ecp['number of core electrons'] = int(match.group(1))

445 ecp['maximum angular momentum number'] = \

446 int(match.group(2))

447 continue

448 match = re.search(r'^(\w(\-\w)?)', line)

449 if match:

450 if len(terms) > 0:

451 # end terms

452 group['terms'] = terms

453 group['number of terms'] = len(terms)

454 terms = []

455 groups.append(group)

456 group = {}

457 # begin group

458 group['title'] = str(match.group(1))

459 continue

460 regex = (r'^\s*([-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)\s+'

461 r'(\d)\s+([-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)')

462 match = re.search(regex, line)

463 if match:

464 terms.append(

465 {

466 'coefficient': float(match.group(1)),

467 'power of r': float(match.group(3)),

468 'exponent': float(match.group(4))

469 }

470 )

471

472

473def read_forces(results, natoms):

474 """Read forces from Turbomole gradient file."""

475 dg = read_data_group('grad')

476 if len(dg) == 0:

477 return

478 file = open('gradient')

479 lines = file.readlines()

480 file.close()

481

482 forces = np.array([[0, 0, 0]])

483

484 nline = len(lines)

485 iline = -1

486

487 for i in range(nline):

488 if 'cycle' in lines[i]:

489 iline = i

490

491 if iline < 0:

492 raise RuntimeError('Please check TURBOMOLE gradients')

493

494 # next line

495 iline += natoms + 1

496 # $end line

497 nline -= 1

498 # read gradients

499 for i in range(iline, nline):

500 line = lines[i].replace('D', 'E')

501 tmp = np.array([[float(f) for f in line.split()[0:3]]])

502 forces = np.concatenate((forces, tmp))

503 # Note the '-' sign for turbomole, to get forces

504 forces = -np.delete(forces, np.s_[0:1], axis=0) * Ha / Bohr

505 results['energy gradient'] = (-forces).tolist()

506 return forces

507

508

509def read_gradient(results):

510 """read all information in file 'gradient'"""

511 grad_string = read_data_group('grad')

512 if len(grad_string) == 0:

513 return

514# try to reuse ase:

515# structures = read('gradient', index=':')

516 lines = grad_string.split('\n')

517 history = []

518 image = {}

519 gradient = []

520 atoms = Atoms()

521 (cycle, energy, norm) = (None, None, None)

522 for line in lines:

523 # cycle lines

524 regex = (

525 r'^\s*cycle =\s*(\d+)\s+'

526 r'SCF energy =\s*([-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)\s+'

527 r'\|dE\/dxyz\| =\s*([-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)'

528 )

529 match = re.search(regex, line)

530 if match:

531 if len(atoms):

532 image['optimization cycle'] = cycle

533 image['total energy'] = energy

534 image['gradient norm'] = norm

535 image['energy gradient'] = gradient

536 history.append(image)

537 image = {}

538 atoms = Atoms()

539 gradient = []

540 cycle = int(match.group(1))

541 energy = float(match.group(2)) * Ha

542 norm = float(match.group(4)) * Ha / Bohr

543 continue

544 # coordinate lines

545 regex = (

546 r'^\s*([-+]?[0-9]*\.?[0-9]+([eEdD][-+]?[0-9]+)?)'

547 r'\s+([-+]?[0-9]*\.?[0-9]+([eEdD][-+]?[0-9]+)?)'

548 r'\s+([-+]?[0-9]*\.?[0-9]+([eEdD][-+]?[0-9]+)?)'

549 r'\s+(\w+)'

550 )

551 match = re.search(regex, line)

552 if match:

553 x = float(match.group(1)) * Bohr

554 y = float(match.group(3)) * Bohr

555 z = float(match.group(5)) * Bohr

556 symbol = str(match.group(7)).capitalize()

557

558 if symbol == 'Q':

559 symbol = 'X'

560 atoms += Atom(symbol, (x, y, z))

561

562 continue

563 # gradient lines

564 regex = (

565 r'^\s*([-+]?[0-9]*\.?[0-9]+([eEdD][-+]?[0-9]+)?)'

566 r'\s+([-+]?[0-9]*\.?[0-9]+([eEdD][-+]?[0-9]+)?)'

567 r'\s+([-+]?[0-9]*\.?[0-9]+([eEdD][-+]?[0-9]+)?)'

568 )

569 match = re.search(regex, line)

570 if match:

571 gradx = float(match.group(1).replace('D', 'E')) * Ha / Bohr

572 grady = float(match.group(3).replace('D', 'E')) * Ha / Bohr

573 gradz = float(match.group(5).replace('D', 'E')) * Ha / Bohr

574 gradient.append([gradx, grady, gradz])

575

576 image['optimization cycle'] = cycle

577 image['total energy'] = energy

578 image['gradient norm'] = norm

579 image['energy gradient'] = gradient

580 history.append(image)

581 results['geometry optimization history'] = history

582

583

584def read_hessian(results, noproj=False):

585 """Read in the hessian matrix"""

586 results['hessian matrix'] = {}

587 results['hessian matrix']['array'] = []

588 results['hessian matrix']['units'] = '?'

589 results['hessian matrix']['projected'] = True

590 results['hessian matrix']['mass weighted'] = True

591 dg = read_data_group('nvibro')

592 if len(dg) == 0:

593 return

594 nvibro = int(dg.split()[1])

595 results['hessian matrix']['dimension'] = nvibro

596 row = []

597 key = 'hessian'

598 if noproj:

599 key = 'npr' + key

600 results['hessian matrix']['projected'] = False

601 lines = read_data_group(key).split('\n')

602 for line in lines:

603 if key in line:

604 continue

605 fields = line.split()

606 row.extend(fields[2:len(fields)])

607 if len(row) == nvibro:

608 # check whether it is mass-weighted

609 float_row = [float(element) for element in row]

610 results['hessian matrix']['array'].append(float_row)

611 row = []

612

613

614def read_normal_modes(results, noproj=False):

615 """Read in vibrational normal modes"""

616 results['normal modes'] = {}

617 results['normal modes']['array'] = []

618 results['normal modes']['projected'] = True

619 results['normal modes']['mass weighted'] = True

620 results['normal modes']['units'] = '?'

621 dg = read_data_group('nvibro')

622 if len(dg) == 0:

623 return

624 nvibro = int(dg.split()[1])

625 results['normal modes']['dimension'] = nvibro

626 row = []

627 key = 'vibrational normal modes'

628 if noproj:

629 key = 'npr' + key

630 results['normal modes']['projected'] = False

631 lines = read_data_group(key).split('\n')

632 for line in lines:

633 if key in line:

634 continue

635 if '$end' in line:

636 break

637 fields = line.split()

638 row.extend(fields[2:len(fields)])

639 if len(row) == nvibro:

640 # check whether it is mass-weighted

641 float_row = [float(element) for element in row]

642 results['normal modes']['array'].append(float_row)

643 row = []

644

645

646def read_vibrational_reduced_masses(results):

647 """Read vibrational reduced masses"""

648 results['vibrational reduced masses'] = []

649 dg = read_data_group('vibrational reduced masses')

650 if len(dg) == 0:

651 return

652 lines = dg.split('\n')

653 for line in lines:

654 if '$vibrational' in line:

655 continue

656 if '$end' in line:

657 break

658 fields = [float(element) for element in line.split()]

659 results['vibrational reduced masses'].extend(fields)

660

661

662def read_vibrational_spectrum(results, noproj=False):

663 """Read the vibrational spectrum"""

664 results['vibrational spectrum'] = []

665 key = 'vibrational spectrum'

666 if noproj:

667 key = 'npr' + key

668 lines = read_data_group(key).split('\n')

669 for line in lines:

670 dictionary = {}

671 regex = (

672 r'^\s+(\d+)\s+(\S*)\s+([-+]?\d+\.\d*)'

673 r'\s+(\d+\.\d*)\s+(\S+)\s+(\S+)'

674 )

675 match = re.search(regex, line)

676 if match:

677 dictionary['mode number'] = int(match.group(1))

678 dictionary['irreducible representation'] = str(match.group(2))

679 dictionary['frequency'] = {

680 'units': 'cm^-1',

681 'value': float(match.group(3))

682 }

683 dictionary['infrared intensity'] = {

684 'units': 'km/mol',

685 'value': float(match.group(4))

686 }

687

688 if match.group(5) == 'YES':

689 dictionary['infrared active'] = True

690 elif match.group(5) == 'NO':

691 dictionary['infrared active'] = False

692 else:

693 dictionary['infrared active'] = None

694

695 if match.group(6) == 'YES':

696 dictionary['Raman active'] = True

697 elif match.group(6) == 'NO':

698 dictionary['Raman active'] = False

699 else:

700 dictionary['Raman active'] = None

701

702 results['vibrational spectrum'].append(dictionary)

703

704

705def read_ssquare(results):

706 """Read the expectation value of S^2 operator"""

707 s2_string = read_data_group('ssquare from dscf')

708 if s2_string == '':

709 return

710 string = s2_string.split('\n')[1]

711 ssquare = float(re.search(r'^\s*(\d+\.*\d*)', string).group(1))

712 results['ssquare from scf calculation'] = ssquare

713

714

715def read_dipole_moment(results):

716 """Read the dipole moment"""

717 dip_string = read_data_group('dipole')

718 if dip_string == '':

719 return

720 lines = dip_string.split('\n')

721 for line in lines:

722 regex = (

723 r'^\s+x\s+([-+]?\d+\.\d*)\s+y\s+([-+]?\d+\.\d*)'

724 r'\s+z\s+([-+]?\d+\.\d*)\s+a\.u\.'

725 )

726 match = re.search(regex, line)

727 if match:

728 dip_vec = [float(match.group(c)) for c in range(1, 4)]

729 regex = r'^\s+\| dipole \| =\s+(\d+\.*\d*)\s+debye'

730 match = re.search(regex, line)

731 if match:

732 dip_abs_val = float(match.group(1))

733 results['electric dipole moment'] = {}

734 results['electric dipole moment']['vector'] = {

735 'array': dip_vec,

736 'units': 'a.u.'

737 }

738 results['electric dipole moment']['absolute value'] = {

739 'value': dip_abs_val,

740 'units': 'Debye'

741 }

742

743

744def read_charges(filename, natoms):

745 """read partial charges on atoms from an ESP fit"""

746 charges = None

747 if os.path.exists(filename):

748 with open(filename) as infile:

749 lines = infile.readlines()

750 oklines = None

751 for n, line in enumerate(lines):

752 if 'atom radius/au charge' in line:

753 oklines = lines[n + 1:n + natoms + 1]

754 if oklines is not None:

755 qm_charges = [float(line.split()[3]) for line in oklines]

756 charges = np.array(qm_charges)

757 return charges

758

759

760def read_point_charges():

761 """read point charges from previous calculation"""

762 pcs = read_data_group('point_charges')

763 lines = pcs.split('\n')[1:]

764 (charges, positions) = ([], [])

765 for line in lines:

766 columns = [float(col) for col in line.strip().split()]

767 positions.append([col * Bohr for col in columns[0:3]])

768 charges.append(columns[3])

769 return charges, positions