r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import numpy as np

3from ase import Atoms

4from ase.units import Bohr, Ry

9 try:

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

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

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

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')

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])

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()

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]])

179 else:

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

181 return cell