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

57 statements  

« prev     ^ index     » next       coverage.py v7.5.3, created at 2025-08-02 00:12 +0000

1# fmt: off 

2 

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

4computem software. 

5 

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

7documentation of prismatic. 

8 

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

10computem software. 

11""" 

12 

13import numpy as np 

14 

15from ase.atoms import Atoms, symbols2numbers 

16from ase.utils import reader 

17 

18from .utils import verify_cell_for_export, verify_dictionary 

19 

20 

21@reader 

22def read_prismatic(fd): 

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

24 

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

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

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

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

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

30 by: 

31 

32 .. math:: 

33 

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

35 

36 """ 

37 # Read comment: 

38 fd.readline() 

39 

40 # Read unit cell parameters: 

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

42 

43 # Read all data at once 

44 # Use genfromtxt instead of loadtxt to skip last line 

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

46 # Convert from RMS to Debye-Waller factor 

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

48 

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

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

51 cell=cellpar, 

52 ) 

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

54 atoms.set_array('debye_waller_factors', RMS) 

55 

56 return atoms 

57 

58 

59class XYZPrismaticWriter: 

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

61 """ 

62 

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

64 verify_cell_for_export(atoms.get_cell()) 

65 

66 self.atoms = atoms.copy() 

67 self.atom_types = set(atoms.symbols) 

68 self.comments = comments 

69 

70 self.occupancies = self._get_occupancies() 

71 debye_waller_factors = self._get_debye_waller_factors( 

72 debye_waller_factors) 

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

74 

75 def _get_occupancies(self): 

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

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

78 else: 

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

80 

81 return occupancies 

82 

83 def _get_debye_waller_factors(self, DW): 

84 if np.isscalar(DW): 

85 if len(self.atom_types) > 1: 

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

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

88 'be provided for each atom using a ' 

89 'dictionary.') 

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

91 elif isinstance(DW, dict): 

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

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

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

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

96 else: 

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

98 if name in self.atoms.arrays: 

99 DW = self.atoms.get_array(name) 

100 

101 if DW is None: 

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

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

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

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

106 'the `Atoms` object.') 

107 

108 return DW 

109 

110 def _get_file_header(self): 

111 # 1st line: comment line 

112 if self.comments is None: 

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

114 len(self.atoms), 

115 self.atoms.get_chemical_formula()) 

116 else: 

117 s = self.comments 

118 

119 s = s.strip() 

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

121 # 2nd line: lattice parameter 

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

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

124 

125 return s 

126 

127 def write_to_file(self, f): 

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

129 self.atoms.positions.T, 

130 self.occupancies, 

131 self.RMS) 

132 ).T 

133 

134 np.savetxt(fname=f, 

135 X=data_array, 

136 fmt='%.6g', 

137 header=self._get_file_header(), 

138 newline='\n', 

139 footer='-1', 

140 comments='' 

141 ) 

142 

143 

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

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

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

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

148 these array will be written to the file. 

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

150 

151 Parameters: 

152 

153 atoms: Atoms object 

154 

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

156 If float, set this value to all 

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

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

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

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

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

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

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

164 

165 .. math:: 

166 

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

168 

169 Default is None. 

170 

171 comment: str (optional) 

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

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

174 

175 """ 

176 

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

178 writer.write_to_file(fd)