Coverage for gpaw/scf.py: 99%
159 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1import time
3import numpy as np
5from gpaw import KohnShamConvergenceError
6from gpaw.convergence_criteria import check_convergence
7from gpaw.forces import calculate_forces
8from gpaw.directmin.scf_helper import do_if_converged, check_eigensolver_state
11class SCFLoop:
12 """Self-consistent field loop."""
13 def __init__(self, criteria, maxiter=100, niter_fixdensity=None):
14 self.criteria = criteria
15 self.maxiter = maxiter
16 self.niter_fixdensity = niter_fixdensity
17 self.niter = None
18 self.reset()
19 self.converged = False
20 self.eigensolver_name = None
21 self.fix_fermi_level = False
23 def __str__(self):
24 s = 'Convergence criteria:\n'
25 for criterion in self.criteria.values():
26 if criterion.description is not None:
27 s += ' ' + criterion.description + '\n'
28 s += f' Maximum number of scf [iter]ations: {self.maxiter}'
29 s += ("\n (Square brackets indicate name in SCF output, whereas a 'c'"
30 " in\n the SCF output indicates the quantity has converged.)\n")
31 return s
33 def write(self, writer):
34 writer.write(converged=self.converged)
36 def read(self, reader):
37 if reader.version < 4:
38 self.converged = reader.scf.converged
39 else:
40 self.converged = True
42 def reset(self):
43 for criterion in self.criteria.values():
44 criterion.reset()
45 self.converged = False
46 self.eigensolver_name = None
48 def irun(self, wfs, ham, dens, log, callback):
50 self.eigensolver_name = getattr(wfs.eigensolver, "name", None)
51 check_eigensolver_state(self.eigensolver_name, wfs, ham, dens, log)
52 self.niter = 1
53 converged = False
55 while self.niter <= self.maxiter:
56 restart = self.iterate_eigensolver(wfs, ham, dens, log)
58 ctx = self.check_convergence(
59 dens, ham, wfs, log, callback)
60 yield ctx
62 if restart:
63 log('MOM has detected variational collapse, '
64 'occupied orbitals have changed')
66 converged = (self.converged and
67 self.niter >= self.niter_fixdensity)
68 if converged:
69 do_if_converged(
70 self.eigensolver_name, wfs, ham, dens, log
71 )
72 break
74 self.update_ham_and_dens(wfs, ham, dens)
75 self.niter += 1
77 # Don't fix the density in the next step.
78 self.niter_fixdensity = 0
80 if not converged:
81 self.not_converged(dens, ham, wfs, log)
83 def log(self, log, converged_items, entries, context):
84 """Output from each iteration."""
85 write_iteration(self.criteria, converged_items, entries, context, log)
87 def check_convergence(self, dens, ham, wfs, log, callback):
89 context = SCFEvent(dens=dens, ham=ham, wfs=wfs, niter=self.niter,
90 log=log)
92 # Converged?
93 # entries: for log file, per criteria
94 # converged_items: True/False, per criteria
95 self.converged, converged_items, entries = check_convergence(
96 self.criteria, context)
98 callback(self.niter)
99 self.log(log, converged_items, entries, context)
100 return context
102 def not_converged(self, dens, ham, wfs, log):
104 context = SCFEvent(dens=dens, ham=ham, wfs=wfs, niter=self.niter,
105 log=log)
106 eigerr = self.criteria['eigenstates'].get_error(context)
107 if not np.isfinite(eigerr):
108 msg = 'Not enough bands for ' + wfs.eigensolver.nbands_converge
109 log(msg)
110 log.fd.flush()
111 raise KohnShamConvergenceError(msg)
112 log(oops, flush=True)
113 raise KohnShamConvergenceError(
114 'Did not converge! See text output for help.')
116 def iterate_eigensolver(self, wfs, ham, dens, log):
118 restart = False
119 if self.eigensolver_name == 'etdm-lcao':
120 wfs.eigensolver.iterate(ham, wfs, dens)
121 restart = wfs.eigensolver.check_mom(wfs, dens)
122 e_entropy = 0.0
123 kin_en_using_band = False
124 elif self.eigensolver_name == 'etdm-fdpw':
125 if not wfs.eigensolver.initialized:
126 wfs.eigensolver.initialize_dm_helper(wfs, ham, dens, log)
127 wfs.eigensolver.iterate(ham, wfs, dens, log)
128 restart = wfs.eigensolver.check_restart(wfs)
129 e_entropy = 0.0
130 kin_en_using_band = False
131 else:
132 wfs.eigensolver.iterate(ham, wfs)
133 e_entropy = wfs.calculate_occupation_numbers(self.fix_fermi_level)
134 kin_en_using_band = True
136 if hasattr(wfs.eigensolver, 'e_sic'):
137 e_sic = wfs.eigensolver.e_sic
138 else:
139 e_sic = 0.0
141 ham.get_energy(
142 e_entropy, wfs, kin_en_using_band=kin_en_using_band, e_sic=e_sic)
144 return restart
146 def update_ham_and_dens(self, wfs, ham, dens):
148 to_update = self.niter > self.niter_fixdensity and not dens.fixed
149 if self.eigensolver_name == 'etdm-lcao' \
150 or self.eigensolver_name == 'etdm-fdpw' \
151 or not to_update:
152 ham.npoisson = 0
153 else:
154 dens.update(wfs)
155 ham.update(dens)
158def write_iteration(criteria, converged_items, entries, ctx, log):
159 custom = (set(criteria) -
160 {'energy', 'eigenstates', 'density'})
162 eigensolver_name = getattr(ctx.wfs.eigensolver, "name", None)
163 print_iloop = False
164 if eigensolver_name == 'etdm-fdpw':
165 if ctx.wfs.eigensolver.iloop is not None or \
166 ctx.wfs.eigensolver.outer_iloop is not None:
167 print_iloop = True
169 if ctx.niter == 1:
170 header1 = (' {:<4s} {:>8s} {:>12s} '
171 .format('iter', 'time', 'total'))
172 header2 = (' {:>4s} {:>8s} {:>12s} '
173 .format('', '', 'energy'))
174 header1 += 'log10-change:'
175 header2 += ' eigst dens '
176 for name in custom:
177 criterion = criteria[name]
178 header1 += ' ' * 7
179 header2 += f'{criterion.tablename:>5s} '
180 if ctx.wfs.nspins == 2:
181 header1 += '{:>8s} '.format('magmom')
182 header2 += '{:>8s} '.format('')
183 if print_iloop:
184 header1 += '{:>12s} '.format('iter')
185 header2 += '{:>12s} '.format('inner loop')
186 log(header1.rstrip())
187 log(header2.rstrip())
189 def format_conv(fmt: str, name: str) -> str:
190 """Add "c" to number and color it green if converged."""
191 txt = fmt.format(entries.get(name, ''))
192 if converged_items.get(name):
193 return log.green + txt + log.reset + 'c '
194 return txt + ' '
196 # Iterations and time.
197 now = time.localtime()
198 line = ('iter:{:4d} {:02d}:{:02d}:{:02d} '
199 .format(ctx.niter, *now[3:6]))
201 # Energy.
202 line += format_conv('{:>12s}', 'energy')
204 # Eigenstates.
205 line += format_conv('{:>6s}', 'eigenstates')
207 # Density.
208 line += format_conv('{:>5s}', 'density')
210 # Custom criteria (optional).
211 for name in custom:
212 line += format_conv('{:>5s}', name)
214 # Magnetic moment (optional).
215 if ctx.wfs.nspins == 2 or not ctx.wfs.collinear:
216 totmom_v, _ = ctx.dens.calculate_magnetic_moments()
217 if ctx.wfs.collinear:
218 line += f' {totmom_v[2]:+.4f}'
219 else:
220 line += ' {:+.1f},{:+.1f},{:+.1f}'.format(*totmom_v)
222 # Inner loop etdm
223 if print_iloop:
224 iloop_counter = (ctx.wfs.eigensolver.eg_count_iloop +
225 ctx.wfs.eigensolver.eg_count_outer_iloop)
226 line += ('{:12d}'.format(iloop_counter))
228 log(line.rstrip())
229 log.fd.flush()
232class SCFEvent:
233 """Object to pass the state of the SCF cycle to a convergence-checking
234 function."""
236 def __init__(self, dens, ham, wfs, niter, log):
237 self.dens = dens
238 self.ham = ham
239 self.wfs = wfs
240 self.niter = niter
241 self.log = log
243 def calculate_forces(self):
244 with self.wfs.timer('Forces'):
245 F_av = calculate_forces(self.wfs, self.dens, self.ham)
246 return F_av
249oops = """
250Did not converge!
252Here are some tips:
2541) Make sure the geometry and spin-state is physically sound.
2552) Use less aggressive density mixing.
2563) Solve the eigenvalue problem more accurately at each scf-step.
2574) Use a smoother distribution function for the occupation numbers.
2585) Try adding more empty states.
2596) Use enough k-points.
2607) Don't let your structure optimization algorithm take too large steps.
2618) Solve the Poisson equation more accurately.
2629) Better initial guess for the wave functions.
264See details here:
266 https://gpaw.readthedocs.io/documentation/convergence.html
268"""