1import numpy as np 

2 

3from ase import Atoms 

4from ase.units import Bohr, Ry 

5from ase.utils import reader, writer 

6 

7 

8def read_scf(filename): 

9 try: 

10 with open(filename + '.scf', 'r') as fd: 

11 pip = fd.readlines() 

12 ene = [] 

13 for line in pip: 

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

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

16 return ene 

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

18 return None 

19 

20 

21@reader 

22def read_struct(fd, ase=True): 

23 pip = fd.readlines() 

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

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

26 cell = np.zeros(6) 

27 for i in range(6): 

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

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

30 if lattice == 'P ': 

31 lattice = 'P' 

32 elif lattice == 'H ': 

33 lattice = 'P' 

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

35 elif lattice == 'R ': 

36 lattice = 'R' 

37 elif lattice == 'F ': 

38 lattice = 'F' 

39 elif lattice == 'B ': 

40 lattice = 'I' 

41 elif lattice == 'CXY': 

42 lattice = 'C' 

43 elif lattice == 'CXZ': 

44 lattice = 'B' 

45 elif lattice == 'CYZ': 

46 lattice = 'A' 

47 else: 

48 raise RuntimeError('TEST needed') 

49 pos = np.array([]) 

50 atomtype = [] 

51 rmt = [] 

52 neq = np.zeros(nat) 

53 iline = 4 

54 indif = 0 

55 for iat in range(nat): 

56 indifini = indif 

57 if len(pos) == 0: 

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

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

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

61 else: 

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

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

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

65 axis=0) 

66 indif += 1 

67 iline += 1 

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

69 iline += 1 

70 for ieq in range(1, int(neq[iat])): 

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

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

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

74 axis=0) 

75 indif += 1 

76 iline += 1 

77 for i in range(indif - indifini): 

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

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

80 iline += 4 

81 if ase: 

82 cell2 = coorsys(cell) 

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

84 atoms.set_cell(cell2, scale_atoms=True) 

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

86 if lattice == 'R': 

87 atoms.set_cell(cell2, scale_atoms=True) 

88 else: 

89 atoms.set_cell(cell2) 

90 return atoms 

91 else: 

92 return cell, lattice, pos, atomtype, rmt 

93 

94 

95@writer 

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

97 atoms = atoms2.copy() 

98 atoms.wrap() 

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

100 nat = len(atoms) 

101 if rmt is None: 

102 rmt = [2.0] * nat 

103 fd.write(lattice + 

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

105 cell = atoms.get_cell() 

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

107 cell2 = cellconst(metT) 

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

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

110 if zza is None: 

111 zza = atoms.get_atomic_numbers() 

112 for ii in range(nat): 

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

114 pos = atoms.get_scaled_positions()[ii] 

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

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

117 zz = zza[ii] 

118 if zz > 71: 

119 ro = 0.000005 

120 elif zz > 36: 

121 ro = 0.00001 

122 elif zz > 18: 

123 ro = 0.00005 

124 else: 

125 ro = 0.0001 

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

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

128 fd.write('LOCAL ROT MATRIX: %9.7f %9.7f %9.7f\n' % (1.0, 0.0, 0.0)) 

129 fd.write(' %9.7f %9.7f %9.7f\n' % (0.0, 1.0, 0.0)) 

130 fd.write(' %9.7f %9.7f %9.7f\n' % (0.0, 0.0, 1.0)) 

131 fd.write(' 0\n') 

132 

133 

134def cellconst(metT): 

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

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

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

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

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

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

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

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

143 

144 

145def coorsys(latconst): 

146 a = latconst[0] 

147 b = latconst[1] 

148 c = latconst[2] 

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

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

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

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

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

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

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

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

157 ]).transpose() 

158 

159 

160def c2p(lattice): 

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

162 if lattice == 'P': 

163 cell = np.eye(3) 

164 elif lattice == 'F': 

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

166 elif lattice == 'I': 

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

168 elif lattice == 'C': 

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

170 elif lattice == 'B': 

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

172 elif lattice == 'A': 

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

174 elif lattice == 'R': 

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

176 [-1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0], 

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

178 

179 else: 

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

181 return cell