Coverage for /builds/ase/ase/ase/io/wien2k.py: 67.83%

1# fmt: off 

2 

3import numpy as np 

4 

5from ase import Atoms 

6from ase.units import Bohr, Ry 

7from ase.utils import reader, writer 

8 

9 

10def read_scf(filename): 

11 try: 

12 with open(filename + '.scf') as fd: 

13 pip = fd.readlines() 

14 ene = [] 

15 for line in pip: 

16 if line[0:4] == ':ENE': 

17 ene.append(float(line[43:59]) * Ry) 

18 return ene 

19 except Exception: # XXX Which kind of exception exactly? 

20 return None 

21 

22 

23@reader 

24def read_struct(fd, ase=True): 

25 pip = fd.readlines() 

26 lattice = pip[1][0:3] 

27 nat = int(pip[1][27:30]) 

28 cell = np.zeros(6) 

29 for i in range(6): 

30 cell[i] = float(pip[3][0 + i * 10:10 + i * 10]) 

31 cell[0:3] = cell[0:3] * Bohr 

32 if lattice == 'P ': 

33 lattice = 'P' 

34 elif lattice == 'H ': 

35 lattice = 'P' 

36 cell[3:6] = [90.0, 90.0, 120.0] 

37 elif lattice == 'R ': 

38 lattice = 'R' 

39 elif lattice == 'F ': 

40 lattice = 'F' 

41 elif lattice == 'B ': 

42 lattice = 'I' 

43 elif lattice == 'CXY': 

44 lattice = 'C' 

45 elif lattice == 'CXZ': 

46 lattice = 'B' 

47 elif lattice == 'CYZ': 

48 lattice = 'A' 

49 else: 

50 raise RuntimeError('TEST needed') 

51 pos = np.array([]) 

52 atomtype = [] 

53 rmt = [] 

54 neq = np.zeros(nat) 

55 iline = 4 

56 indif = 0 

57 for iat in range(nat): 

58 indifini = indif 

59 if len(pos) == 0: 

60 pos = np.array([[float(pip[iline][12:22]), 

61 float(pip[iline][25:35]), 

62 float(pip[iline][38:48])]]) 

63 else: 

64 pos = np.append(pos, np.array([[float(pip[iline][12:22]), 

65 float(pip[iline][25:35]), 

66 float(pip[iline][38:48])]]), 

67 axis=0) 

68 indif += 1 

69 iline += 1 

70 neq[iat] = int(pip[iline][15:17]) 

71 iline += 1 

72 for _ in range(1, int(neq[iat])): 

73 pos = np.append(pos, np.array([[float(pip[iline][12:22]), 

74 float(pip[iline][25:35]), 

75 float(pip[iline][38:48])]]), 

76 axis=0) 

77 indif += 1 

78 iline += 1 

79 for _ in range(indif - indifini): 

80 atomtype.append(pip[iline][0:2].replace(' ', '')) 

81 rmt.append(float(pip[iline][43:48])) 

82 iline += 4 

83 if ase: 

84 cell2 = coorsys(cell) 

85 atoms = Atoms(atomtype, pos, pbc=True) 

86 atoms.set_cell(cell2, scale_atoms=True) 

87 cell2 = np.dot(c2p(lattice), cell2) 

88 if lattice == 'R': 

89 atoms.set_cell(cell2, scale_atoms=True) 

90 else: 

91 atoms.set_cell(cell2) 

92 return atoms 

93 else: 

94 return cell, lattice, pos, atomtype, rmt 

95 

96 

97@writer 

98def write_struct(fd, atoms2=None, rmt=None, lattice='P', zza=None): 

99 atoms = atoms2.copy() 

100 atoms.wrap() 

101 fd.write('ASE generated\n') 

102 nat = len(atoms) 

103 if rmt is None: 

104 rmt = [2.0] * nat 

105 fd.write(lattice + 

106 ' LATTICE,NONEQUIV.ATOMS:%3i\nMODE OF CALC=RELA\n' % nat) 

107 cell = atoms.get_cell() 

108 metT = np.dot(cell, np.transpose(cell)) 

109 cell2 = cellconst(metT) 

110 cell2[0:3] = cell2[0:3] / Bohr 

111 fd.write(('%10.6f' * 6) % tuple(cell2) + '\n') 

112 if zza is None: 

113 zza = atoms.get_atomic_numbers() 

114 for ii in range(nat): 

115 fd.write('ATOM %3i: ' % (ii + 1)) 

116 pos = atoms.get_scaled_positions()[ii] 

117 fd.write('X=%10.8f Y=%10.8f Z=%10.8f\n' % tuple(pos)) 

118 fd.write(' MULT= 1 ISPLIT= 1\n') 

119 zz = zza[ii] 

120 if zz > 71: 

121 ro = 0.000005 

122 elif zz > 36: 

123 ro = 0.00001 

124 elif zz > 18: 

125 ro = 0.00005 

126 else: 

127 ro = 0.0001 

128 fd.write('%-10s NPT=%5i R0=%9.8f RMT=%10.4f Z:%10.5f\n' % 

129 (atoms.get_chemical_symbols()[ii], 781, ro, rmt[ii], zz)) 

130 fd.write(f'LOCAL ROT MATRIX: {1.0:9.7f} {0.0:9.7f} {0.0:9.7f}\n') 

131 fd.write(f' {0.0:9.7f} {1.0:9.7f} {0.0:9.7f}\n') 

132 fd.write(f' {0.0:9.7f} {0.0:9.7f} {1.0:9.7f}\n') 

133 fd.write(' 0\n') 

134 

135 

136def cellconst(metT): 

137 """ metT=np.dot(cell,cell.T) """ 

138 aa = np.sqrt(metT[0, 0]) 

139 bb = np.sqrt(metT[1, 1]) 

140 cc = np.sqrt(metT[2, 2]) 

141 gamma = np.arccos(metT[0, 1] / (aa * bb)) / np.pi * 180.0 

142 beta = np.arccos(metT[0, 2] / (aa * cc)) / np.pi * 180.0 

143 alpha = np.arccos(metT[1, 2] / (bb * cc)) / np.pi * 180.0 

144 return np.array([aa, bb, cc, alpha, beta, gamma]) 

145 

146 

147def coorsys(latconst): 

148 a = latconst[0] 

149 b = latconst[1] 

150 c = latconst[2] 

151 cal = np.cos(latconst[3] * np.pi / 180.0) 

152 cbe = np.cos(latconst[4] * np.pi / 180.0) 

153 cga = np.cos(latconst[5] * np.pi / 180.0) 

154 sga = np.sin(latconst[5] * np.pi / 180.0) 

155 return np.array([[a, b * cga, c * cbe], 

156 [0, b * sga, c * (cal - cbe * cga) / sga], 

157 [0, 0, c * np.sqrt(1 - cal**2 - cbe**2 - cga**2 + 

158 2 * cal * cbe * cga) / sga] 

159 ]).transpose() 

160 

161 

162def c2p(lattice): 

163 """ apply as eg. cell2 = np.dot(c2p('F'), cell)""" 

164 if lattice == 'P': 

165 cell = np.eye(3) 

166 elif lattice == 'F': 

167 cell = np.array([[0.0, 0.5, 0.5], [0.5, 0.0, 0.5], [0.5, 0.5, 0.0]]) 

168 elif lattice == 'I': 

169 cell = np.array([[-0.5, 0.5, 0.5], [0.5, -0.5, 0.5], [0.5, 0.5, -0.5]]) 

170 elif lattice == 'C': 

171 cell = np.array([[0.5, 0.5, 0.0], [0.5, -0.5, 0.0], [0.0, 0.0, -1.0]]) 

172 elif lattice == 'B': 

173 cell = np.array([[0.5, 0.0, 0.5], [0.0, 1.0, 0.0], [0.5, 0.0, -0.5]]) 

174 elif lattice == 'A': 

175 cell = np.array([[-1.0, 0.0, 0.0], [0.0, -0.5, 0.5], [0.0, 0.5, 0.5]]) 

176 elif lattice == 'R': 

177 cell = np.array([[2.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0], 

178 [-1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0], 

179 [-1.0 / 3.0, -2.0 / 3.0, 1.0 / 3.0]]) 

180 

181 else: 

182 raise ValueError('lattice is ' + lattice + '!') 

183 return cell