Coverage for gpaw/borncharges.py: 94%

1import numpy as np 

2from gpaw.mpi import world 

3from gpaw.berryphase import polarization_phase, ionic_phase 

4from ase.parallel import paropen, parprint 

5from ase.io.jsonio import write_json, read_json 

6from pathlib import Path 

7 

8 

9def born_charges_wf(atoms, calc, delta=0.01, cleanup=False, 

10 ionic_only=False, out='born_charges.json'): 

11 

12 # generate displacement dictionary 

13 disps_av = _all_disp(atoms, delta) 

14 

15 # carry out polarization phase calculation 

16 # for each displacement 

17 phases_c = {} 

18 for dlabel in disps_av: 

19 ia, iv, sign, delta = disps_av[dlabel] 

20 atoms_d = displace_atom(atoms, ia, iv, sign, delta) 

21 check_distance_to_non_pbc_boundary(atoms_d) 

22 

23 if not ionic_only: 

24 

25 # proper polarization phase calculation 

26 gpw_wfs = Path(dlabel + '.gpw') 

27 berryname = Path(dlabel + '_berry-phases.json') 

28 if not berryname.is_file(): 

29 if not gpw_wfs.is_file(): 

30 

31 # run calculations 

32 atoms_d.calc = calc 

33 assert is_symmetry_off(atoms_d.calc), 'Set symmetry off' 

34 atoms_d.get_potential_energy() 

35 

36 # write wavefunctions 

37 atoms_d.calc.write(gpw_wfs, 'all') 

38 

39 # dict with entries phase_c, electronic_phase_c 

40 # atomic_phase_c, dipole_moment_c 

41 phase_c = polarization_phase(gpw_wfs, comm=world) 

42 

43 # only master rank should write 

44 with paropen(berryname, 'w') as fd: 

45 write_json(fd, phase_c) 

46 

47 else: 

48 # all ranks can read 

49 with open(berryname, 'r') as fd: 

50 phase_c = read_json(fd) 

51 

52 if cleanup: 

53 if berryname.is_file(): 

54 # remove gpw file 

55 if world.rank == 0: 

56 gpw_wfs.unlink() 

57 else: 

58 # only atomic contribution considered 

59 # for unexpensive testing only 

60 phase_c = ionic_phase(atoms_d) 

61 

62 phases_c[dlabel] = phase_c['phase_c'] 

63 

64 results = born_charges(atoms, disps_av, phases_c, check=(not ionic_only)) 

65 with paropen(out, 'w') as fd: 

66 write_json(fd, results) 

67 

68 return results 

69 

70 

71def is_symmetry_off(calc): 

72 params = calc.parameters 

73 if calc.old: 

74 if 'symmetry' in params: 

75 return params['symmetry'] == 'off' 

76 else: 

77 return False 

78 else: 

79 # new: 

80 return (not params.symmetry.point_group and 

81 not params.symmetry.time_reversal) 

82 

83 

84def born_charges(atoms, disps_av, phases_c, check=True): 

85 

86 natoms = len(atoms) 

87 cell_cv = atoms.get_cell() 

88 vol = abs(np.linalg.det(cell_cv)) 

89 sym_a = atoms.get_chemical_symbols() 

90 

91 ndisp = len(disps_av) 

92 parprint('Not using symmetry: ndisp:', ndisp) 

93 

94 # obtain phi(dr) map 

95 phi_ascv = np.zeros((natoms, 2, 3, 3), float) 

96 for dlabel in disps_av: 

97 ia, iv, sign, delta = disps_av[dlabel] 

98 isign = [None, 1, 0][sign] 

99 phi_ascv[ia, isign, :, iv] = phases_c[dlabel] 

100 

101 # calculate dphi / dr 

102 # exploit +- displacement 

103 dphi_acv = phi_ascv[:, 1] - phi_ascv[:, 0] 

104 # mod 2 pi 

105 mod_acv = np.round(dphi_acv / (2 * np.pi)) * 2 * np.pi 

106 dphi_acv -= mod_acv 

107 # transform to cartesian 

108 dphi_avv = np.array([np.dot(dphi_cv.T, cell_cv).T for dphi_cv in dphi_acv]) 

109 dphi_dr_avv = dphi_avv / (2.0 * delta) 

110 

111 # calculate polarization change and born charges 

112 dP_dr_avv = dphi_dr_avv / (2 * np.pi * vol) 

113 Z_avv = dP_dr_avv * vol 

114 

115 if check: 

116 # check acoustic sum rule: sum_a Z_aij = 0 for all i,j 

117 asr_vv = np.sum(Z_avv, axis=0) 

118 asr_dev = np.abs(asr_vv).max() / natoms 

119 assert asr_dev < 1e-1, f'Acoustic sum rule violated: {asr_vv}' 

120 

121 # correct to match acoustic sum rule 

122 Z_avv -= asr_vv[None, :, :] / natoms 

123 

124 results = {'Z_avv': Z_avv, 'sym_a': sym_a} 

125 

126 return results 

127 

128 

129def _cartesian_label(ia, iv, sign): 

130 """Generate name from (ia, iv, sign). 

131 ia ... atomic_index 

132 iv ... cartesian_index 

133 sign ... +- 

134 """ 

135 

136 sym_v = 'xyz'[iv] 

137 sym_s = ' +-'[sign] 

138 return f'{ia}{sym_v}{sym_s}' 

139 

140 

141def _all_avs(atoms): 

142 """Generate ia, iv, sign for all displacements.""" 

143 for ia in range(len(atoms)): 

144 for iv in range(3): 

145 for sign in [-1, 1]: 

146 yield (ia, iv, sign) 

147 

148 

149def _all_disp(atoms, delta): 

150 all_disp = {} 

151 for dd, avs in enumerate(_all_avs(atoms)): 

152 dd = int(dd) 

153 lavs = _cartesian_label(*avs) 

154 label = f'disp_{dd:03d}_' + lavs 

155 all_disp[label] = (*avs, delta) 

156 return all_disp 

157 

158 

159def displace_atom(atoms, ia, iv, sign, delta): 

160 new_atoms = atoms.copy() 

161 pos_av = new_atoms.get_positions() 

162 pos_av[ia, iv] += sign * delta 

163 new_atoms.set_positions(pos_av) 

164 return new_atoms 

165 

166 

167def check_distance_to_non_pbc_boundary(atoms, eps=1): 

168 dist_a = distance_to_non_pbc_boundary(atoms) 

169 if dist_a is not None and np.any(dist_a < eps): 

170 raise AtomsTooCloseToBoundary( 

171 'The atoms are too close to a non-pbc boundary ' 

172 'which creates problems when using a dipole correction. ' 

173 f'Please center the atoms in the unit-cell. Distances: {dist_a}.' 

174 ) 

175 

176 

177def distance_to_non_pbc_boundary(atoms): 

178 pbc_c = atoms.get_pbc() 

179 if pbc_c.all(): 

180 return None 

181 cell_cv = atoms.get_cell() 

182 pos_ac = atoms.get_scaled_positions() 

183 pos_ac -= np.round(pos_ac) 

184 posnonpbc_av = np.dot(pos_ac[:, ~pbc_c], cell_cv[~pbc_c]) 

185 dist_to_cell_edge_a = np.linalg.norm(posnonpbc_av, axis=1) 

186 return dist_to_cell_edge_a 

187 

188 

189class AtomsTooCloseToBoundary(Exception): 

190 pass