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

2 

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 

11 

12 

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 

18 

19 

20def read_data_group(data_group): 

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

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

23 

24 

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 

47 

48 

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 

64 

65 

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 

78 

79 

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 

94 

95 

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 

105 

106 

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 

116 

117 

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 

131 

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 

158 

159 

160def read_run_parameters(results): 

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

162 

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 

181 

182 nvibro = read_data_group('nvibro') 

183 if nvibro: 

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

185 

186 

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') 

196 

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

198 

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 

211 

212 

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) 

245 

246 

247def read_mos(results): 

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

249 from files mos, alpha and beta""" 

250 

251 results['molecular orbitals'] = [] 

252 mos = results['molecular orbitals'] 

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

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

255 converged = None 

256 

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 

306 

307 

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 ) 

381 

382 

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 ) 

467 

468 

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() 

477 

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

479 

480 nline = len(lines) 

481 iline = -1 

482 

483 for i in range(nline): 

484 if 'cycle' in lines[i]: 

485 iline = i 

486 

487 if iline < 0: 

488 raise RuntimeError('Please check TURBOMOLE gradients') 

489 

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 

503 

504 

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() 

553 

554 if symbol == 'Q': 

555 symbol = 'X' 

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

557 

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]) 

571 

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 

578 

579 

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 = [] 

608 

609 

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 = [] 

640 

641 

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) 

656 

657 

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 } 

683 

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 

690 

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 

697 

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

699 

700 

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 

709 

710 

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 } 

738 

739 

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 

754 

755 

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