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 and return values for corresponding quantity 

45 Parameters 

46 ========== 

47 filename: ACE-Molecule log file. 

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

49 

50 Returns 

51 ======= 

52 - quantity = 'excitation-energy': 

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

54 without IndexError. 

55 - quantity = 'energy': 

56 returns energy as float value. 

57 - quantity = 'forces': 

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

59 - quantity = 'atoms': 

60 returns ASE atoms object. 

61 ''' 

62 data = parse_geometry(filename) 

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

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

65 atoms = Atoms(atom_symbol, positions=positions) 

66 energy = None 

67 forces = None 

68 excitation_energy = None 

69# results = {} 

70# if len(results)<1: 

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

72 lines = fd.readlines() 

73 

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

75 line = lines[i].split() 

76 if len(line) > 2: 

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

78 energy = float(line[3]) 

79 break 

80 # energy must be modified, hartree to eV 

81 energy *= ase.units.Hartree 

82 

83 forces = [] 

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

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

86 endline_num = i 

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

88 startline_num = i + 2 

89 for j in range(startline_num, endline_num): 

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

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

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

93 break 

94 if not len(forces) > 0: 

95 forces = None 

96 

97 # Set calculator to 

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

99 atoms.calc = calc 

100 

101 results = {} 

102 results['energy'] = energy 

103 results['atoms'] = atoms 

104 results['forces'] = forces 

105 results['excitation-energy'] = excitation_energy 

106 return results 

107 

108 

109def read_acemolecule_input(filename): 

110 '''Reads a ACE-Molecule input file 

111 Parameters 

112 ========== 

113 filename: ACE-Molecule input file name 

114 

115 Returns 

116 ======= 

117 ASE atoms object containing geometry only. 

118 ''' 

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

120 for line in fd: 

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

122 geometryfile = line.split()[1] 

123 break 

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

125 return atoms