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

1import numpy as np 

2 

3import ase.units 

4from ase.atoms import Atoms 

5from ase.calculators.singlepoint import SinglePointCalculator 

6from ase.data import chemical_symbols 

7from ase.io import read 

8 

9 

10def parse_geometry(filename): 

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

12 Parameters 

13 ========== 

14 filename: ACE-Molecule log file. 

15 

16 Returns 

17 ======= 

18 Dictionary of parsed geometry data. 

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

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

21 ''' 

22 with open(filename) as fd: 

23 lines = fd.readlines() 

24 start_line = 0 

25 end_line = 0 

26 for i, line in enumerate(lines): 

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

28 start_line = i 

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

30 end_line = i 

31 if not start_line == end_line: 

32 break 

33 atoms = [] 

34 positions = [] 

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

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

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

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

39 positions.append(xyz) 

40 

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

42 

43 

44def read_acemolecule_out(filename): 

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

46 

47 Parameters 

48 ========== 

49 filename: ACE-Molecule log file. 

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

51 

52 Returns 

53 ======= 

54 - quantity = 'excitation-energy': 

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

56 without IndexError. 

57 - quantity = 'energy': 

58 returns energy as float value. 

59 - quantity = 'forces': 

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

61 - quantity = 'atoms': 

62 returns ASE atoms object. 

63 ''' 

64 data = parse_geometry(filename) 

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

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

67 atoms = Atoms(atom_symbol, positions=positions) 

68 energy = None 

69 forces = None 

70 excitation_energy = None 

71# results = {} 

72# if len(results)<1: 

73 with open(filename) as fd: 

74 lines = fd.readlines() 

75 

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

77 line = lines[i].split() 

78 if len(line) > 2: 

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

80 energy = float(line[3]) 

81 break 

82 # energy must be modified, hartree to eV 

83 energy *= ase.units.Hartree 

84 

85 forces = [] 

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

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

88 endline_num = i 

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

90 startline_num = i + 2 

91 for j in range(startline_num, endline_num): 

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

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

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

95 break 

96 if not len(forces) > 0: 

97 forces = None 

98 

99 # Set calculator to 

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

101 atoms.calc = calc 

102 

103 results = {} 

104 results['energy'] = energy 

105 results['atoms'] = atoms 

106 results['forces'] = forces 

107 results['excitation-energy'] = excitation_energy 

108 return results 

109 

110 

111def read_acemolecule_input(filename): 

112 '''Reads a ACE-Molecule input file 

113 Parameters 

114 ========== 

115 filename: ACE-Molecule input file name 

116 

117 Returns 

118 ======= 

119 ASE atoms object containing geometry only. 

120 ''' 

121 with open(filename) as fd: 

122 for line in fd: 

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

124 geometryfile = line.split()[1] 

125 break 

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

127 return atoms