1import numpy as np 

2import ase.units 

3from ase.atoms import Atoms 

4from ase.calculators.singlepoint import SinglePointCalculator 

5from ase.io import read 

6from ase.data import chemical_symbols 

7 

8 

9def parse_geometry(filename): 

10 '''Read atoms geometry from ACE-Molecule log file and put it to self.data. 

11 Parameters 

12 ========== 

13 filename: ACE-Molecule log file. 

14 

15 Returns 

16 ======= 

17 Dictionary of parsed geometry data. 

18 retval["Atomic_numbers"]: list of atomic numbers 

19 retval["Positions"]: list of [x, y, z] coordinates for each atoms. 

20 ''' 

21 with open(filename, 'r') as fd: 

22 lines = fd.readlines() 

23 start_line = 0 

24 end_line = 0 

25 for i, line in enumerate(lines): 

26 if line == '==================== Atoms =====================\n': 

27 start_line = i 

28 if start_line != 0 and len(line.split('=')) > 3: 

29 end_line = i 

30 if not start_line == end_line: 

31 break 

32 atoms = [] 

33 positions = [] 

34 for i in range(start_line + 1, end_line): 

35 atomic_number = lines[i].split()[0] 

36 atoms.append(str(chemical_symbols[int(atomic_number)])) 

37 xyz = [float(n) for n in lines[i].split()[1:4]] 

38 positions.append(xyz) 

39 

40 return {"Atomic_numbers": atoms, "Positions": positions} 

41 

42 

43def read_acemolecule_out(filename): 

44 '''Interface to ACEMoleculeReader, return values for corresponding quantity 

45 

46 Parameters 

47 ========== 

48 filename: ACE-Molecule log file. 

49 quantity: One of atoms, energy, forces, excitation-energy. 

50 

51 Returns 

52 ======= 

53 - quantity = 'excitation-energy': 

54 returns None. This is placeholder function to run TDDFT calculations 

55 without IndexError. 

56 - quantity = 'energy': 

57 returns energy as float value. 

58 - quantity = 'forces': 

59 returns force of each atoms as numpy array of shape (natoms, 3). 

60 - quantity = 'atoms': 

61 returns ASE atoms object. 

62 ''' 

63 data = parse_geometry(filename) 

64 atom_symbol = np.array(data["Atomic_numbers"]) 

65 positions = np.array(data["Positions"]) 

66 atoms = Atoms(atom_symbol, positions=positions) 

67 energy = None 

68 forces = None 

69 excitation_energy = None 

70# results = {} 

71# if len(results)<1: 

72 with open(filename, 'r') as fd: 

73 lines = fd.readlines() 

74 

75 for i in range(len(lines) - 1, 1, -1): 

76 line = lines[i].split() 

77 if len(line) > 2: 

78 if line[0] == 'Total' and line[1] == 'energy': 

79 energy = float(line[3]) 

80 break 

81 # energy must be modified, hartree to eV 

82 energy *= ase.units.Hartree 

83 

84 forces = [] 

85 for i in range(len(lines) - 1, 1, -1): 

86 if '!============================' in lines[i]: 

87 endline_num = i 

88 if '! Force:: List of total force in atomic unit' in lines[i]: 

89 startline_num = i + 2 

90 for j in range(startline_num, endline_num): 

91 forces.append(lines[j].split()[3:6]) 

92 convert = ase.units.Hartree / ase.units.Bohr 

93 forces = np.array(forces, dtype=float) * convert 

94 break 

95 if not len(forces) > 0: 

96 forces = None 

97 

98 # Set calculator to 

99 calc = SinglePointCalculator(atoms, energy=energy, forces=forces) 

100 atoms.calc = calc 

101 

102 results = {} 

103 results['energy'] = energy 

104 results['atoms'] = atoms 

105 results['forces'] = forces 

106 results['excitation-energy'] = excitation_energy 

107 return results 

108 

109 

110def read_acemolecule_input(filename): 

111 '''Reads a ACE-Molecule input file 

112 Parameters 

113 ========== 

114 filename: ACE-Molecule input file name 

115 

116 Returns 

117 ======= 

118 ASE atoms object containing geometry only. 

119 ''' 

120 with open(filename, 'r') as fd: 

121 for line in fd: 

122 if len(line.split('GeometryFilename')) > 1: 

123 geometryfile = line.split()[1] 

124 break 

125 atoms = read(geometryfile, format='xyz') 

126 return atoms