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

49 

50 

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 

66 

67 

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 

80 

81 

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 

96 

97 

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 

107 

108 

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 

118 

119 

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 

133 

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 

160 

161 

162def read_run_parameters(results): 

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

164 

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 

183 

184 nvibro = read_data_group('nvibro') 

185 if nvibro: 

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

187 

188 

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

198 

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

200 

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 

213 

214 

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) 

247 

248 

249def read_mos(results): 

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

251 from files mos, alpha and beta""" 

252 

253 results['molecular orbitals'] = [] 

254 mos = results['molecular orbitals'] 

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

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

257 converged = None 

258 

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 

308 

309 

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 ) 

383 

384 

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 ) 

469 

470 

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

479 

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

481 

482 nline = len(lines) 

483 iline = -1 

484 

485 for i in range(nline): 

486 if 'cycle' in lines[i]: 

487 iline = i 

488 

489 if iline < 0: 

490 raise RuntimeError('Please check TURBOMOLE gradients') 

491 

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 

505 

506 

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

555 

556 if symbol == 'Q': 

557 symbol = 'X' 

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

559 

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

573 

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 

580 

581 

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

610 

611 

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

642 

643 

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) 

658 

659 

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 } 

685 

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 

692 

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 

699 

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

701 

702 

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 

711 

712 

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 } 

740 

741 

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 

756 

757 

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