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

1# fmt: off

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

5import os

6import re

7import subprocess

8import warnings

10import numpy as np

12from ase import Atom, Atoms

13from ase.calculators.calculator import ReadError

14from ase.units import Bohr, Ha

17def execute_command(args):

18 """execute commands like sdg, eiger"""

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

20 stdout, _stderr = proc.communicate()

21 return stdout

24def read_data_group(data_group):

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

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

29def parse_data_group(dgr, dg_name):

30 """parse a data group"""

31 if len(dgr) == 0:

32 return None

33 dg_key = '$' + dg_name

34 if not dgr.startswith(dg_key):

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

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

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

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

39 return ndg

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

41 result = {}

42 lines = ndg.split(lsep)

43 for line in lines:

44 if len(line) == 0:

45 continue

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

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

48 if len(fields) == 2:

49 result[fields[0]] = fields[1]

50 elif len(fields) == 1:

51 result[fields[0]] = True

52 return result

55def read_output(regex, path):

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

57 hitlist = []

58 checkfiles = []

59 for filename in os.listdir(path):

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

61 checkfiles.append(filename)

62 for filename in checkfiles:

63 with open(filename) as f:

64 lines = f.readlines()

65 for line in lines:

66 match = re.search(regex, line)

67 if match:

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

69 return hitlist

72def read_version(path):

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

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

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

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

77 version = list(set(versions))

78 elif len(versions) == 0:

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

80 version = None

81 else:

82 version = versions[0]

83 return version

86def read_datetime(path):

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

88 from the tm output if stored in a file

89 """

90 datetimes = read_output(

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

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

93 if len(datetimes) == 0:

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

95 datetime = None

96 else:

97 # take the most recent time stamp

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

99 return datetime

100

101

102def read_runtime(path):

103 """read the total runtime of calculations"""

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

105 if len(hits) == 0:

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

107 runtime = None

108 else:

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

110 return runtime

111

112

113def read_hostname(path):

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

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

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

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

118 hostname = list(set(hostnames))

119 else:

120 hostname = hostnames[0]

121 return hostname

122

123

124def read_convergence(restart, parameters):

125 """perform convergence checks"""

126 if restart:

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

128 return False

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

130 return False

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

132 return False

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

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

135 return False

136 return True

137

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

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

140 return True

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

142 # check whether a failed scf convergence is the reason

143 checkfiles = []

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

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

146 checkfiles.append(filename)

147 for filename in checkfiles:

148 for line in open(filename):

149 if 'SCF FAILED TO CONVERGE' in line:

150 # scf did not converge in some jobex iteration

151 if filename == 'job.last':

152 raise RuntimeError('scf failed to converge')

153 else:

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

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

156 return False

157 else:

158 raise RuntimeError('error during geometry optimization')

159 else:

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

161 raise RuntimeError('scf failed to converge')

162 else:

163 return True

164

165

166def read_run_parameters(results):

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

168

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

170 results['calculation parameters'] = {}

171 parameters = results['calculation parameters']

172 dg = read_data_group('symmetry')

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

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

175 # Gaussian function type

176 gt = read_data_group('pople')

177 if gt == '':

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

179 else:

180 gt = gt.split()[1]

181 if gt == 'AO':

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

183 elif gt == 'CAO':

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

185 else:

186 parameters['gaussian type'] = None

187

188 nvibro = read_data_group('nvibro')

189 if nvibro:

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

191

192

193def read_energy(results, post_HF):

194 """Read energy from Turbomole energy file."""

195 try:

196 with open('energy') as enf:

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

198 except OSError:

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

200 if text == '':

201 raise ReadError('empty energy file')

202

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

204

205 for line in lines:

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

207 break

208 elif line.startswith('$'):

209 pass

210 else:

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

212 if post_HF:

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

214 # update energy units

215 e_total = energy_tmp * Ha

216 results['total energy'] = e_total

217

218

219def read_occupation_numbers(results):

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

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

222 return

223 mos = results['molecular orbitals']

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

225 for line in lines:

226 regex = (

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

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

229 )

230 match = re.search(regex, line)

231 if match:

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

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

234 spin = 'alpha'

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

236 spin = 'beta'

237 else:

238 spin = None

239 ar_index = next(

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

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

242 molecular_orbital['spin'] == spin)

243 )

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

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

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

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

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

249 else:

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

251

252

253def read_mos(results):

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

255 from files mos, alpha and beta"""

256

257 results['molecular orbitals'] = []

258 mos = results['molecular orbitals']

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

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

261 converged = None

262

263 for index, keyword in enumerate(keywords):

264 flen = None

265 mo = {}

266 orbitals_coefficients_line = []

267 mo_string = read_data_group(keyword)

268 if mo_string == '':

269 continue

270 mo_string += '\n$end'

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

272 for line in lines:

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

274 continue

275 if 'eigenvalue' in line:

276 if len(orbitals_coefficients_line) != 0:

277 mo['eigenvector'] = orbitals_coefficients_line

278 mos.append(mo)

279 mo = {}

280 orbitals_coefficients_line = []

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

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

283 match = re.search(regex, line)

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

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

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

287 mo['eigenvalue'] = eig

288 mo['spin'] = spin[index]

289 mo['degeneracy'] = 1

290 continue

291 if keyword in line:

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

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

294 match = re.search(regex, line)

295 if match:

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

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

298 'scfconv' not in line):

299 converged = False

300 continue

301 if '$end' in line:

302 if len(orbitals_coefficients_line) != 0:

303 mo['eigenvector'] = orbitals_coefficients_line

304 mos.append(mo)

305 break

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

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

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

309 for f in sfields]

310 orbitals_coefficients_line += ffields

311 return converged

312

313

314def read_basis_set(results):

315 """read the basis set"""

316 results['basis set'] = []

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

318 bsf = read_data_group('basis')

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

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

321 basis_set = {}

322 functions = []

323 function = {}

324 primitives = []

325 read_tag = False

326 read_data = False

327 for line in lines:

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

329 continue

330 if '$basis' in line:

331 continue

332 if '$end' in line:

333 break

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

335 continue

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

337 if read_tag:

338 read_tag = False

339 read_data = True

340 else:

341 if read_data:

342 # end primitives

343 function['primitive functions'] = primitives

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

345 primitives = []

346 functions.append(function)

347 function = {}

348 # end contracted

349 basis_set['functions'] = functions

350 functions = []

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

352 basis_set = {}

353 read_data = False

354 read_tag = True

355 continue

356 if read_tag:

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

358 if match:

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

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

361 else:

362 raise RuntimeError('error reading basis set')

363 else:

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

365 if match:

366 if len(primitives) > 0:

367 # end primitives

368 function['primitive functions'] = primitives

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

370 primitives = []

371 functions.append(function)

372 function = {}

373 # begin contracted

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

375 continue

376 regex = (

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

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

379 )

380 match = re.search(regex, line)

381 if match:

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

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

384 primitives.append(

385 {'exponent': exponent, 'coefficient': coefficient}

386 )

387

388

389def read_ecps(results):

390 """read the effective core potentials"""

391 ecpf = read_data_group('ecp')

392 if not bool(len(ecpf)):

393 results['ecps'] = None

394 results['ecps formatted'] = None

395 return

396 results['ecps'] = []

397 results['ecps formatted'] = {}

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

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

400 ecp = {}

401 groups = []

402 group = {}

403 terms = []

404 read_tag = False

405 read_data = False

406 for line in lines:

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

408 continue

409 if '$ecp' in line:

410 continue

411 if '$end' in line:

412 break

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

414 continue

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

416 if read_tag:

417 read_tag = False

418 read_data = True

419 else:

420 if read_data:

421 # end terms

422 group['terms'] = terms

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

424 terms = []

425 groups.append(group)

426 group = {}

427 # end group

428 ecp['groups'] = groups

429 groups = []

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

431 ecp = {}

432 read_data = False

433 read_tag = True

434 continue

435 if read_tag:

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

437 if match:

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

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

440 else:

441 raise RuntimeError('error reading ecp')

442 else:

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

444 match = re.search(regex, line)

445 if match:

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

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

448 int(match.group(2))

449 continue

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

451 if match:

452 if len(terms) > 0:

453 # end terms

454 group['terms'] = terms

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

456 terms = []

457 groups.append(group)

458 group = {}

459 # begin group

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

461 continue

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

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

464 match = re.search(regex, line)

465 if match:

466 terms.append(

467 {

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

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

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

471 }

472 )

473

474

475def read_forces(results, natoms):

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

477 dg = read_data_group('grad')

478 if len(dg) == 0:

479 return None

480 with open('gradient') as file:

481 lines = file.readlines()

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