Coverage for /builds/ase/ase/ase/io/acemolecule.py: 98.61%

1# fmt: off 

2 

3import numpy as np 

4 

5import ase.units 

6from ase.atoms import Atoms 

7from ase.calculators.singlepoint import SinglePointCalculator 

8from ase.data import chemical_symbols 

9from ase.io import read 

10 

11 

12def parse_geometry(filename): 

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

14 Parameters 

15 ========== 

16 filename: ACE-Molecule log file. 

17 

18 Returns 

19 ======= 

20 Dictionary of parsed geometry data. 

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

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

23 ''' 

24 with open(filename) as fd: 

25 lines = fd.readlines() 

26 start_line = 0 

27 end_line = 0 

28 for i, line in enumerate(lines): 

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

30 start_line = i 

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

32 end_line = i 

33 if start_line != end_line: 

34 break 

35 atoms = [] 

36 positions = [] 

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

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

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

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

41 positions.append(xyz) 

42 

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

44 

45 

46def read_acemolecule_out(filename): 

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

48 

49 Parameters 

50 ========== 

51 filename: ACE-Molecule log file. 

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

53 

54 Returns 

55 ======= 

56 - quantity = 'excitation-energy': 

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

58 without IndexError. 

59 - quantity = 'energy': 

60 returns energy as float value. 

61 - quantity = 'forces': 

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

63 - quantity = 'atoms': 

64 returns ASE atoms object. 

65 ''' 

66 data = parse_geometry(filename) 

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

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

69 atoms = Atoms(atom_symbol, positions=positions) 

70 energy = None 

71 forces = None 

72 excitation_energy = None 

73# results = {} 

74# if len(results)<1: 

75 with open(filename) as fd: 

76 lines = fd.readlines() 

77 

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

79 line = lines[i].split() 

80 if len(line) > 2: 

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

82 energy = float(line[3]) 

83 break 

84 # energy must be modified, hartree to eV 

85 energy *= ase.units.Hartree 

86 

87 forces = [] 

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

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

90 endline_num = i 

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

92 startline_num = i + 2 

93 for j in range(startline_num, endline_num): 

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

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

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

97 break 

98 if len(forces) <= 0: 

99 forces = None 

100 

101 # Set calculator to 

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

103 atoms.calc = calc 

104 

105 results = {} 

106 results['energy'] = energy 

107 results['atoms'] = atoms 

108 results['forces'] = forces 

109 results['excitation-energy'] = excitation_energy 

110 return results 

111 

112 

113def read_acemolecule_input(filename): 

114 '''Reads a ACE-Molecule input file 

115 Parameters 

116 ========== 

117 filename: ACE-Molecule input file name 

118 

119 Returns 

120 ======= 

121 ASE atoms object containing geometry only. 

122 ''' 

123 with open(filename) as fd: 

124 for line in fd: 

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

126 geometryfile = line.split()[1] 

127 break 

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

129 return atoms