Coverage for /builds/kinetik161/ase/ase/io/prismatic.py: 98.25%

57 statements  

« prev     ^ index     » next       coverage.py v7.2.7, created at 2023-12-10 11:04 +0000

1"""Module to read and write atoms in xtl file format for the prismatic and 

2computem software. 

3 

4See https://prism-em.com/docs-inputs for an example of this format and the 

5documentation of prismatic. 

6 

7See https://sourceforge.net/projects/computem/ for the source code of the 

8computem software. 

9""" 

10 

11import numpy as np 

12 

13from ase.atoms import Atoms, symbols2numbers 

14from ase.utils import reader 

15 

16from .utils import verify_cell_for_export, verify_dictionary 

17 

18 

19@reader 

20def read_prismatic(fd): 

21 r"""Import prismatic and computem xyz input file as an Atoms object. 

22 

23 Reads cell, atom positions, occupancies and Debye Waller factor. 

24 The occupancy values and the Debye Waller factors are obtained using the 

25 `get_array` method and the `occupancies` and `debye_waller_factors` keys, 

26 respectively. The root means square (RMS) values from the 

27 prismatic/computem xyz file are converted to Debye-Waller factors (B) in Ų 

28 by: 

29 

30 .. math:: 

31 

32 B = RMS^2 * 8\pi^2 

33 

34 """ 

35 # Read comment: 

36 fd.readline() 

37 

38 # Read unit cell parameters: 

39 cellpar = [float(i) for i in fd.readline().split()] 

40 

41 # Read all data at once 

42 # Use genfromtxt instead of loadtxt to skip last line 

43 read_data = np.genfromtxt(fname=fd, skip_footer=1) 

44 # Convert from RMS to Debye-Waller factor 

45 RMS = read_data[:, 5]**2 * 8 * np.pi**2 

46 

47 atoms = Atoms(symbols=read_data[:, 0], 

48 positions=read_data[:, 1:4], 

49 cell=cellpar, 

50 ) 

51 atoms.set_array('occupancies', read_data[:, 4]) 

52 atoms.set_array('debye_waller_factors', RMS) 

53 

54 return atoms 

55 

56 

57class XYZPrismaticWriter: 

58 """ See the docstring of the `write_prismatic` function. 

59 """ 

60 

61 def __init__(self, atoms, debye_waller_factors=None, comments=None): 

62 verify_cell_for_export(atoms.get_cell()) 

63 

64 self.atoms = atoms.copy() 

65 self.atom_types = set(atoms.symbols) 

66 self.comments = comments 

67 

68 self.occupancies = self._get_occupancies() 

69 debye_waller_factors = self._get_debye_waller_factors( 

70 debye_waller_factors) 

71 self.RMS = np.sqrt(debye_waller_factors / (8 * np.pi**2)) 

72 

73 def _get_occupancies(self): 

74 if 'occupancies' in self.atoms.arrays: 

75 occupancies = self.atoms.get_array('occupancies', copy=False) 

76 else: 

77 occupancies = np.ones_like(self.atoms.numbers) 

78 

79 return occupancies 

80 

81 def _get_debye_waller_factors(self, DW): 

82 if np.isscalar(DW): 

83 if len(self.atom_types) > 1: 

84 raise ValueError('This cell contains more then one type of ' 

85 'atoms and the Debye-Waller factor needs to ' 

86 'be provided for each atom using a ' 

87 'dictionary.') 

88 DW = np.ones_like(self.atoms.numbers) * DW 

89 elif isinstance(DW, dict): 

90 verify_dictionary(self.atoms, DW, 'DW') 

91 # Get the arrays of DW from mapping the DW defined by symbol 

92 DW = {symbols2numbers(k)[0]: v for k, v in DW.items()} 

93 DW = np.vectorize(DW.get)(self.atoms.numbers) 

94 else: 

95 for name in ['DW', 'debye_waller_factors']: 

96 if name in self.atoms.arrays: 

97 DW = self.atoms.get_array(name) 

98 

99 if DW is None: 

100 raise ValueError('Missing Debye-Waller factors. It can be ' 

101 'provided as a dictionary with symbols as key or ' 

102 'can be set for each atom by using the ' 

103 '`set_array("debye_waller_factors", values)` of ' 

104 'the `Atoms` object.') 

105 

106 return DW 

107 

108 def _get_file_header(self): 

109 # 1st line: comment line 

110 if self.comments is None: 

111 s = "{} atoms with chemical formula: {}.".format( 

112 len(self.atoms), 

113 self.atoms.get_chemical_formula()) 

114 else: 

115 s = self.comments 

116 

117 s = s.strip() 

118 s += " generated by the ase library.\n" 

119 # 2nd line: lattice parameter 

120 s += "{} {} {}".format( 

121 *self.atoms.cell.cellpar()[:3]) 

122 

123 return s 

124 

125 def write_to_file(self, f): 

126 data_array = np.vstack((self.atoms.numbers, 

127 self.atoms.positions.T, 

128 self.occupancies, 

129 self.RMS) 

130 ).T 

131 

132 np.savetxt(fname=f, 

133 X=data_array, 

134 fmt='%.6g', 

135 header=self._get_file_header(), 

136 newline='\n', 

137 footer='-1', 

138 comments='' 

139 ) 

140 

141 

142def write_prismatic(fd, *args, **kwargs): 

143 r"""Write xyz input file for the prismatic and computem software. The cell 

144 needs to be orthorhombric. If the cell contains the `occupancies` and 

145 `debye_waller_factors` arrays (see the `set_array` method to set them), 

146 these array will be written to the file. 

147 If the occupancies is not specified, the default value will be set to 0. 

148 

149 Parameters: 

150 

151 atoms: Atoms object 

152 

153 debye_waller_factors: float or dictionary of float or None (optional). 

154 If float, set this value to all 

155 atoms. If dictionary, each atom needs to specified with the symbol as 

156 key and the corresponding Debye-Waller factor as value. 

157 If None, the `debye_waller_factors` array of the Atoms object needs to 

158 be set (see the `set_array` method to set them), otherwise raise a 

159 ValueError. Since the prismatic/computem software use root means square 

160 (RMS) displacements, the Debye-Waller factors (B) needs to be provided 

161 in Ų and these values are converted to RMS displacement by: 

162 

163 .. math:: 

164 

165 RMS = \sqrt{\frac{B}{8\pi^2}} 

166 

167 Default is None. 

168 

169 comment: str (optional) 

170 Comments to be written in the first line of the file. If not 

171 provided, write the total number of atoms and the chemical formula. 

172 

173 """ 

174 

175 writer = XYZPrismaticWriter(*args, **kwargs) 

176 writer.write_to_file(fd)