收敛测试

先看一份完成的 SCF 输出的开头。这里能确认版本、4 个 MPI 进程、读到的输入文件，以及实际使用的晶格、原子数和截断。

[preston@preston-System-Product-Name si-pbe]$ head -n 82 scf/scf.out

     Program PWSCF v.7.5 starts on 22Sep2026 at 21:32:23

     This program is part of the open-source Quantum ESPRESSO suite
     for quantum simulation of materials; please cite
         "P. Giannozzi et al., J. Phys.:Condens. Matter 21 395502 (2009);
         "P. Giannozzi et al., J. Phys.:Condens. Matter 29 465901 (2017);
         "P. Giannozzi et al., J. Chem. Phys. 152 154105 (2020);
          URL http://www.quantum-espresso.org",
     in publications or presentations arising from this work. More details at
     http://www.quantum-espresso.org/quote

     Parallel version (MPI), running on     4 processors

     MPI processes distributed on     1 nodes
     8639 MiB available memory on the printing compute node when the environment starts

     Reading input from scf.in

     Current dimensions of program PWSCF are:
     Max number of different atomic species (ntypx) = 10
     Max number of k-points (npk) =  40000
     Max angular momentum in pseudopotentials (lmaxx) =  4

     R & G space division:  proc/nbgrp/npool/nimage =       4
     Subspace diagonalization in iterative solution of the eigenvalue problem:
     a serial algorithm will be used


     Parallelization info
     --------------------
     sticks:   dense  smooth     PW     G-vecs:    dense   smooth      PW
     Min         571     214     63                18093     4156     682
     Max         572     217     64                18095     4157     686
     Sum        2287     859    253                72377    16625    2733

     Using Slab Decomposition



     bravais-lattice index     =            2
     lattice parameter (alat)  =      10.2000  a.u.
     unit-cell volume          =     265.3020 (a.u.)^3
     number of atoms/cell      =            2
     number of atomic types    =            1
     number of electrons       =         8.00
     number of Kohn-Sham states=            4
     kinetic-energy cutoff     =      60.0000  Ry
     charge density cutoff     =     640.0000  Ry
     scf convergence threshold =      1.0E-10
     mixing beta               =       0.7000
     number of iterations used =            8  plain     mixing
     Exchange-correlation= PBE
                           (   1   4   3   4   0   0   0)

     celldm(1)=  10.200000  celldm(2)=   0.000000  celldm(3)=   0.000000
     celldm(4)=   0.000000  celldm(5)=   0.000000  celldm(6)=   0.000000

     crystal axes: (cart. coord. in units of alat)
               a(1) = (  -0.500000   0.000000   0.500000 )
               a(2) = (   0.000000   0.500000   0.500000 )
               a(3) = (  -0.500000   0.500000   0.000000 )

     reciprocal axes: (cart. coord. in units 2 pi/alat)
               b(1) = ( -1.000000 -1.000000  1.000000 )
               b(2) = (  1.000000  1.000000  1.000000 )
               b(3) = ( -1.000000  1.000000 -1.000000 )


     PseudoPot. # 1 for Si read from file:
     ../pseudo/Si.pbe-n-rrkjus_psl.1.0.0.UPF
     MD5 check sum: fa25574f73a70a4139f2adfbefec430c
     Pseudo is Ultrasoft + core correction, Zval =  4.0
     Generated using "atomic" code by A. Dal Corso  v.6.3
     Using radial grid of 1141 points,  6 beta functions with:
                l(1) =   0
                l(2) =   0
                l(3) =   1
                l(4) =   1
                l(5) =   2
                l(6) =   2
     Q(r) pseudized with 0 coefficients
[preston@preston-System-Product-Name si-pbe]$

继续向下，电子迭代从初始电荷开始。estimated scf accuracy 会随迭代下降；中间没有感叹号的能量是当前迭代值，不能和已经完成的另一个计算直接放进收敛表。

[preston@preston-System-Product-Name si-pbe]$ grep -A28 'Self-consistent Calculation' scf/scf.out
     Self-consistent Calculation

     iteration #  1     ecut=    60.00 Ry     beta= 0.70
     Davidson diagonalization with overlap
     ethr =  1.00E-02,  avg # of iterations =  2.0

     Threshold (ethr) on eigenvalues was too large:
     Diagonalizing with lowered threshold

     Davidson diagonalization with overlap
     ethr =  6.40E-04,  avg # of iterations =  1.5

     total cpu time spent up to now is        1.2 secs

     total energy              =     -22.83581101 Ry
     estimated scf accuracy    <       0.05540955 Ry

     iteration #  2     ecut=    60.00 Ry     beta= 0.70
     Davidson diagonalization with overlap
     ethr =  6.93E-04,  avg # of iterations =  1.0

     total cpu time spent up to now is        1.5 secs

     total energy              =     -22.83770169 Ry
     estimated scf accuracy    <       0.00288592 Ry

     iteration #  3     ecut=    60.00 Ry     beta= 0.70
     Davidson diagonalization with overlap
     ethr =  3.61E-05,  avg # of iterations =  2.8
[preston@preston-System-Product-Name si-pbe]$

提取带感叹号的最终总能量时，同时核对相邻的 convergence has been achieved 信息和末轮误差。这里还会给出分项能量和力；一份完整 OUT 不只有一列最终数字。

[preston@preston-System-Product-Name si-pbe]$ grep -A20 '!    total energy' scf/scf.out
!    total energy              =     -22.83859230 Ry
     estimated scf accuracy    <          4.3E-11 Ry

     The total energy is the sum of the following terms:
     one-electron contribution =       5.30781076 Ry
     hartree contribution      =       1.08523150 Ry
     xc contribution           =     -12.33187599 Ry
     ewald contribution        =     -16.89975857 Ry

     convergence has been achieved in   9 iterations

     Forces acting on atoms (cartesian axes, Ry/au):

     atom    1 type  1   force =    -0.00000000   -0.00000000   -0.00000000
     atom    2 type  1   force =     0.00000000    0.00000000    0.00000000

     Total force =     0.000000     Total SCF correction =     0.000000


     Computing stress (Cartesian axis) and pressure

[preston@preston-System-Product-Name si-pbe]$

末尾用于确认程序走到结束；PWSCF 行给出 CPU 和 WALL 时间，JOB DONE. 与前面的电子收敛信息一起读取。最早几次 scf.err 保留 X11 授权提示；对应 OUT 的电子迭代已收敛并正常结束。

[preston@preston-System-Product-Name si-pbe]$ tail -n 12 scf/scf.out
     interpolate  :      0.12s CPU      0.13s WALL (      10 calls)

     Parallel routines

     PWSCF        :      9.22s CPU     13.73s WALL


   This run was terminated on:  21:32:37  22Sep2026

=------------------------------------------------------------------------------=
   JOB DONE.
=------------------------------------------------------------------------------=
[preston@preston-System-Product-Name si-pbe]$

# Exact existing CSV; no smoothing, fitting, interpolation or new DFT.
set datafile separator ','
set encoding utf8
set border linewidth 1
set tics out nomirror
set style line 1 lc rgb '#205493' lw 1.6 pt 7 ps 1.0
set style line 2 lc rgb '#ae3b2d' lw 1.5 dt 2
set style line 3 lc rgb '#237747' lw 1.6 pt 6 ps 1.6
unset key
set logscale y
set yrange [0.001:120]
set ylabel '|E-E(ref)| (meV/atom)'
parameters='ecutwfc ecutrho kmesh'
labels='ecutwfc ecutrho k-mesh'
do for [fmt in 'svg png pdf'] {
 if (fmt eq 'svg') { set terminal svg size 1140,430 enhanced font 'Arial,13' }
 if (fmt eq 'png') { set terminal pngcairo size 1140,430 enhanced font 'DejaVu Sans,13' }
 if (fmt eq 'pdf') { set terminal pdfcairo size 11.4in,4.3in enhanced font 'DejaVu Sans,13' }
 set output 'convergence.'.fmt
 set multiplot layout 1,3 margins 0.08,0.985,0.16,0.82 spacing 0.095,0.04 title 'Si: independent scans; highest measured point in each scan is its reference'
 do for [i=1:3] {
  par=word(parameters,i)
  set xlabel (i<3 ? word(labels,i).' (Ry)' : 'Uniform mesh edge N (N x N x N)')
  set title sprintf('(%s) %s',word('a b c',i),word(labels,i))
  if (i==1) { set xrange [40:80]; set xtics 40,10,80; set ylabel "|E-E(ref)| (meV/atom)" }
  if (i==2) { set xrange [320:640]; set xtics 320,160,640; unset ylabel }
  if (i==3) { set xrange [4:14]; set xtics 4,2,14; unset ylabel }
  set key top right font ',10'
  plot 'convergence.csv' using (strcol(1) eq par && $6>0 ? $3 : 1/0):6 with linespoints ls 1 title 'Sampled residual', \
       1 with lines ls 2 title '1 meV/atom', \
       'convergence.csv' using (strcol(1) eq par && strcol(8) eq 'True' && $6>0 ? $3 : 1/0):6 with points ls 3 title 'Selected'
 }
 unset multiplot
 unset output
}

#!/usr/bin/env python3
"""Read the supplied Si QE SCF files and tabulate total-energy changes.

Usage: python3 analyze_si_convergence.py si-pbe --outdir reproduced
Python standard library only. SHA256SUMS.raw is read beside this script.
"""
from __future__ import annotations
import argparse
import csv
import hashlib
import json
import math
import re
from pathlib import Path

RY_EV = 13.6056931229905
AXES = {
    'ecutwfc': ['cutoff40', 'cutoff50', 'cutoff60', 'cutoff70', 'cutoff80'],
    'ecutrho': ['rho320', 'rho480', 'scf'],
    'kmesh': ['k4', 'k6', 'k8', 'k10', 'k12', 'k14'],
}
LOGISTICS = {'prefix', 'outdir', 'pseudo_dir', 'wfcdir'}


def sha256(path):
    return hashlib.sha256(path.read_bytes()).hexdigest()


def verify_sources(package):
    records = []
    for line in (package / 'SHA256SUMS.raw').read_text().splitlines():
        digest, name = line.split('  ', 1)
        path = package / name
        if not path.is_file() or sha256(path) != digest:
            raise ValueError(f'Raw-file checksum mismatch: {name}')
        records.append({'file': name, 'sha256': digest})
    return records


def read_input(path):
    lines = [re.sub(r'\s+', ' ', line.split('!', 1)[0].strip()).lower()
             for line in path.read_text().splitlines() if line.split('!', 1)[0].strip()]
    params = {}
    for line in lines:
        match = re.fullmatch(r'([a-z][a-z0-9_]*)\s*=\s*(.*?)\s*,?', line)
        if match:
            params[match[1]] = match[2].rstrip(',').strip().strip("'\"")
    index = lines.index('k_points automatic')
    grid = tuple(int(value) for value in lines[index + 1].split())
    if len(grid) != 6 or len(set(grid[:3])) != 1:
        raise ValueError(f'{path}: expected an n x n x n automatic grid')
    return params, grid, lines, index


def read_run(directory, source_names, package):
    inp = directory / 'scf.in'
    out = directory / 'scf.out'
    err = directory / 'scf.err'
    for path in (inp, out, err):
        if path.exists() and path.relative_to(package).as_posix() not in source_names:
            raise ValueError(f'Raw file missing from SHA256SUMS.raw: {path.name}')
    params, grid, lines, index = read_input(inp)
    text = out.read_text(errors='replace') if out.is_file() else ''
    energies = re.findall(r'^\s*!\s*total energy\s*=\s*([-+\d.eEdD]+)\s+Ry', text, re.M)
    completed = 'JOB DONE.' in text and len(re.findall('convergence has been achieved', text)) == 1
    reasons = []
    if not out.is_file():
        reasons.append('missing scf.out')
    if len(energies) != 1 or not completed:
        reasons.append('one converged SCF energy and normal end required')
    energy = float(energies[0].replace('D', 'E').replace('d', 'e')) if len(energies) == 1 else None
    if energy is not None and not math.isfinite(energy):
        reasons.append('nonfinite energy')
    version = re.search(r'Program PWSCF v\.([^\s]+)', text)
    if text:
        echoes = {'nat': r'number of atoms/cell\s*=\s*(\d+)',
                  'ecutwfc': r'kinetic-energy cutoff\s*=\s*([-+\d.]+)',
                  'ecutrho': r'charge density cutoff\s*=\s*([-+\d.]+)'}
        for key, pattern in echoes.items():
            match = re.search(pattern, text)
            if not match or not math.isclose(float(match[1]), float(params[key]), abs_tol=1e-6):
                reasons.append(f'{key} input/output mismatch')
        if not version or version[1] != '7.5':
            reasons.append('expected QE 7.5')
        if not re.search(r'Exchange-correlation\s*=\s*PBE', text):
            reasons.append('expected PBE output')
    if params.get('calculation') != 'scf' or int(params['nat']) != 2:
        reasons.append('expected a two-atom SCF input')
    stderr = [line.strip() for line in err.read_text(errors='replace').splitlines() if line.strip()] if err.is_file() else []
    errclass = ('empty' if not stderr else 'X11 authorization notice' if all(
        line == 'Authorization required, but no authorization protocol specified' for line in stderr
    ) else 'inspect stderr')
    if errclass == 'inspect stderr':
        reasons.append('unclassified stderr')
    return {'directory': directory.name, 'params': params, 'grid': grid,
            'lines': lines, 'grid_line': index + 1, 'energy': energy,
            'complete': not reasons, 'reasons': '; '.join(reasons),
            'stderr': errclass, 'input_sha256': sha256(inp),
            'output_sha256': sha256(out) if out.is_file() else ''}


def signature(run, axis):
    kept = []
    for index, line in enumerate(run['lines']):
        assignment = re.match(r'([a-z][a-z0-9_]*)\s*=', line)
        if assignment and (assignment[1] in LOGISTICS or assignment[1] == axis):
            continue
        if axis == 'kmesh' and index == run['grid_line']:
            kept.append(' '.join(str(value) for value in run['grid'][3:]))
        else:
            kept.append(line)
    return tuple(kept)


def write_csv(path, rows):
    with path.open('w', newline='', encoding='utf-8') as handle:
        writer = csv.DictWriter(handle, fieldnames=list(rows[0]))
        writer.writeheader()
        writer.writerows(rows)


def scf_history(path):
    history = []
    for match in re.finditer(r'iteration #\s*(\d+)(.*?)(?=iteration #|\Z)', path.read_text(), re.S):
        energy = re.search(r'(?:!\s*)?total energy\s*=\s*([-+\d.]+)\s+Ry', match[2])
        accuracy = re.search(r'estimated scf accuracy\s*<\s*([-+\d.Ee]+)\s+Ry', match[2])
        if not energy or not accuracy:
            raise ValueError(f'incomplete SCF iteration {match[1]}')
        history.append({'iteration': int(match[1]), 'energy_Ry_per_cell': float(energy[1]),
                        'estimated_scf_accuracy_Ry': float(accuracy[1])})
    return history


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('raw', type=Path)
    parser.add_argument('--outdir', type=Path, default=Path('results'))
    parser.add_argument('--tolerance', type=float, default=1.0, help='meV/atom')
    args = parser.parse_args()
    if not math.isfinite(args.tolerance) or args.tolerance <= 0:
        raise ValueError('tolerance must be finite and positive')
    raw = args.raw.resolve()
    package = Path(__file__).resolve().parent
    sources = verify_sources(package)
    source_names = {row['file'] for row in sources}
    runs = [read_run(path.parent, source_names, package) for path in sorted(raw.glob('*/scf.in'))]
    by_name = {run['directory']: run for run in runs}
    base = by_name['scf']
    for name in ('cutoff60', 'k8'):
        other = by_name[name]
        if signature(other, '') != signature(base, '') or not other['complete'] or not base['complete']:
            raise ValueError(f'{name}: baseline copy differs')
        if abs(other['energy'] - base['energy']) > 1e-8:
            raise ValueError(f'{name}: baseline energies differ beyond output precision')
    rows, selections = [], []
    for axis, names in AXES.items():
        members = [by_name[name] for name in names]
        if not all(run['complete'] for run in members):
            raise ValueError(f'{axis}: a required scan point is incomplete')
        if len({signature(run, axis) for run in members}) != 1:
            raise ValueError(f'{axis}: fixed input settings differ')
        reference = members[-1]['energy']
        series = []
        previous = None
        for run in members:
            setting = run['grid'][0] if axis == 'kmesh' else float(run['params'][axis])
            adjacent = None if previous is None else abs(run['energy'] - previous) * RY_EV * 1000 / 2
            series.append({'parameter': axis, 'directory': run['directory'], 'setting': setting,
                'natoms': 2, 'energy_Ry_per_cell': run['energy'],
                'delta_meV_per_atom': abs(run['energy'] - reference) * RY_EV * 1000 / 2,
                'adjacent_meV_per_atom': adjacent, 'selected': False,
                'input_sha256': run['input_sha256'], 'output_sha256': run['output_sha256']})
            previous = run['energy']
        chosen = next(i for i, row in enumerate(series) if row['delta_meV_per_atom'] <= args.tolerance
                      and all(later['adjacent_meV_per_atom'] <= args.tolerance for later in series[i + 1:]))
        series[chosen]['selected'] = True
        selections.append({'parameter': axis, 'selected': series[chosen]['setting'],
            'reference': series[-1]['setting'], 'delta_meV_per_atom': series[chosen]['delta_meV_per_atom']})
        rows.extend(series)
    inventory = [{key:run[key] for key in ('directory','complete','reasons','stderr','input_sha256','output_sha256')} for run in runs]
    args.outdir.mkdir(parents=True, exist_ok=True)
    write_csv(args.outdir / 'convergence.csv', rows)
    write_csv(args.outdir / 'run-inventory.csv', inventory)
    write_csv(args.outdir / 'scf-history.csv', scf_history(raw/'scf/scf.out'))
    summary = {'QE_version': '7.5', 'Ry_to_eV': RY_EV, 'natoms': 2,
        'tolerance_meV_per_atom': args.tolerance, 'candidate_runs': len(runs),
        'complete_runs': sum(run['complete'] for run in runs), 'table_rows': len(rows),
        'reference': 'highest sampled setting in each axis', 'selections': selections,
        'independent_scans': True, 'joint_selected_settings_calculated': False}
    (args.outdir/'summary.json').write_text(json.dumps(summary,ensure_ascii=False,indent=2)+'\n')
    (args.outdir/'source-files.json').write_text(json.dumps(sources,indent=2)+'\n')
    report = ['# Si 固定晶胞总能量扫描', '',
        f'QE 7.5；两个原子/原胞；1 Ry = {RY_EV} eV。', '',
        '| 参数 | 设置 | 总能量 (Ry/原胞) | 与最高点差值 (meV/atom) | 相邻变化 (meV/atom) | 选择 |',
        '| --- | ---: | ---: | ---: | ---: | --- |']
    for row in rows:
        adjacent = '—' if row['adjacent_meV_per_atom'] is None else f"{row['adjacent_meV_per_atom']:.6f}"
        report.append(f"| {row['parameter']} | {row['setting']:g} | {row['energy_Ry_per_cell']:.8f} | {row['delta_meV_per_atom']:.6f} | {adjacent} | {'✓' if row['selected'] else ''} |")
    report.extend(['', f'按 {args.tolerance:g} meV/atom 比较线，选择同时满足参照差和后续相邻差的最低采样点。',
        '参照是每组最高已测设置。三组分别固定其余参数；所选最低设置尚未组合为同一次计算。',
        'conv_thr 控制单次电子自洽；力、应力及后续性质按各自目标量继续比较。', ''])
    (args.outdir/'energy-report.md').write_text('\n'.join(report),encoding='utf-8')
    print(f"raw_checksums={len(sources)} complete_runs={summary['complete_runs']} excluded_runs={len(runs)-summary['complete_runs']} table_rows={len(rows)}")
    for item in selections:
        setting = f"{int(item['selected'])}x{int(item['selected'])}x{int(item['selected'])}" if item['parameter']=='kmesh' else f"{item['selected']:g} Ry"
        print(f"{item['parameter']}: {setting}; difference to reference = {item['delta_meV_per_atom']:.6f} meV/atom")
    print(f'Results: {args.outdir}')


if __name__ == '__main__':
    main()