GPAW状态方程与Δ曲线比较

from pathlib import Path
import os, json, hashlib, numpy as np
from ase.build import bulk
from ase.eos import EquationOfState, birchmurnaghan
from ase.utils.deltacodesdft import delta
from scipy.integrate import quad
from gpaw import GPAW, PW, FermiDirac
R = Path(__file__).resolve().parent / 'results'
R.mkdir()
p = json.loads((R.parent / 'inputs.json').read_text())
checks = []
rows = []

def ck(n, b, **kw):
    checks.append(dict(name=n, passed=bool(b), **kw))
for f, h in p['source_SHA'].items():
    assert hashlib.sha256(Path(f).read_bytes()).hexdigest() == h
for f, h in p['old_sources_SHA'].items():
    assert hashlib.sha256(Path(f).read_bytes()).hexdigest() == h
old = []
for row in p['old_curve']:
    old.append(row)
    ck('old_converged:' + str(row['scale']), row['passed'])
for scale in p['scales']:
    D = R / f'v{scale:.2f}'
    D.mkdir()
    a = bulk('Si', 'diamond', a=5.43)
    a.set_cell(a.cell * scale ** (1 / 3), scale_atoms=True)
    c = GPAW(legacy_gpaw=True, mode=PW(300), xc='PBE', kpts=(4, 4, 4), nbands=12, occupations=FermiDirac(0.01), convergence=dict(energy=1e-06, density=1e-05, eigenstates=1e-08), maxiter=200, txt=str(D / 'native.txt'))
    a.calc = c
    e = float(a.get_potential_energy())
    c.write(str(D / 'state.gpw'), mode='all')
    force = a.get_forces()
    rho = c.get_all_electron_density(gridrefinement=1)
    ne = float(rho.sum() * a.get_volume() / rho.size)
    occs = np.array([c.get_occupation_numbers(kpt=k) for k in range(c.wfs.kd.nibzkpts)])
    ck(str(scale) + ':AE28', abs(ne - 28) < 0.03, value=ne)
    ck(str(scale) + ':valence8', abs(occs.sum() - 8) < 1e-07)
    ck(str(scale) + ':symmetric_force', abs(force).max() < 0.02)
    ck(str(scale) + ':true_volume_not_lattice_scale', abs(a.get_volume() / p['reference_cell_volume_A3'] - scale) < 1e-12)
    ck(str(scale) + ':finite', np.isfinite([e, ne, *force.ravel()]).all() and np.isfinite(rho).all())
    np.savez_compressed(D / 'fields.npz', AE_total_density_e_A3=rho, cell_A=a.cell.array, positions_A=a.positions, forces_eV_A=force, eigenvalues_eV=np.array([c.get_eigenvalues(kpt=k) for k in range(c.wfs.kd.nibzkpts)]), weighted_occupations=occs, IBZ_weights=c.get_k_point_weights())
    row = dict(scale=scale, volume_A3=a.get_volume(), energy_eV=e, atoms=2, AE_electrons=ne, scf_iterations=c.get_number_of_iterations())
    rows.append(row)
    (R / 'PARTIAL-RESULTS.json').write_text(json.dumps(rows, indent=2) + '\n')
    print(scale, e, a.get_volume(), flush=True)
fits = {}
for label, data in [('PW250-k2-reused', old), ('PW300-k4', rows)]:
    v = np.array([r['volume_A3'] / 2 for r in data])
    e = np.array([r['energy_eV'] / 2 for r in data])
    fit = EquationOfState(v, e, eos='birchmurnaghan')
    v0, e0, b0 = fit.fit()
    e0, b0, bp, v0 = map(float, fit.eos_parameters)
    pred = birchmurnaghan(v, *fit.eos_parameters)
    rmse = float(np.sqrt(np.mean((pred - e) ** 2)))
    ck(label + ':minimum_in_sampled_volume_range', min(v) < v0 < max(v))
    ck(label + ':positive_bulk_modulus', b0 > 0)
    ck(label + ':fit_RMSE', rmse < 0.005, value=rmse)
    fits[label] = dict(v0_A3_atom=v0, e0_eV_atom=e0, B0_eV_A3=b0, Bprime=bp, fit_RMSE_eV_atom=rmse, volumes_A3_atom=v.tolist(), energies_eV_atom=e.tolist())
a, b = fits.values()
args = [a['v0_A3_atom'], a['B0_eV_A3'], a['Bprime'], b['v0_A3_atom'], b['B0_eV_A3'], b['Bprime']]
native = float(delta(*args))
vref = (args[0] + args[3]) / 2
lo = 0.94 * vref
hi = 1.06 * vref
v = lo + (np.arange(100) + 0.5) * (hi - lo) / 100
ea = birchmurnaghan(v, 0, args[1], args[2], args[0])
eb = birchmurnaghan(v, 0, args[4], args[5], args[3])
independent = float(np.sqrt(np.mean((ea - eb) ** 2)))
ck('native_ASE100point_RMS_independent', abs(native - independent) < 1e-12)
ck('native_symmetric_exchange', abs(native - float(delta(*args[3:], *args[:3]))) < 1e-12)
ck('native_self_delta_zero', abs(float(delta(*args[:3], *args[:3]))) < 1e-14)

def diff(v):
    return float(birchmurnaghan(v, 0, args[1], args[2], args[0]) - birchmurnaghan(v, 0, args[4], args[5], args[3]))
mean = quad(diff, lo, hi, epsabs=1e-13)[0] / (hi - lo)
centered = np.sqrt(quad(lambda x: (diff(x) - mean) ** 2, lo, hi, epsabs=1e-14)[0] / (hi - lo))
dense = np.linspace(lo, hi, 10001)
de = birchmurnaghan(dense, 0, args[1], args[2], args[0]) - birchmurnaghan(dense, 0, args[4], args[5], args[3])
dm = np.trapezoid(de, dense) / (hi - lo)
dc = np.sqrt(np.trapezoid((de - dm) ** 2, dense) / (hi - lo))
shift = 4.7
shiftmean = mean + shift
shifted = np.sqrt(quad(lambda x: (diff(x) + shift - shiftmean) ** 2, lo, hi, epsabs=1e-14)[0] / (hi - lo))
ck('explicit_mean_centered_continuum_dense_quadrature', abs(centered - dc) < 1e-07)
ck('explicit_mean_centered_energy_offset_invariance', abs(centered - shifted) < 1e-12)
for label, data in [('coarse', old), ('fine', rows)]:
    sampled = [r['volume_A3'] / 2 for r in data]
    ck(label + ':whole_integral_window_inside_actual_samples', min(sampled) <= lo and max(sampled) >= hi, sampled_range=[min(sampled), max(sampled)], integral_range=[lo, hi])
for f, h in p['old_sources_SHA'].items():
    ck('original_EOS_source_immutable:' + f, hashlib.sha256(Path(f).read_bytes()).hexdigest() == h)
for f, h in p['source_SHA'].items():
    ck('library_immutable:' + f, hashlib.sha256(Path(f).read_bytes()).hexdigest() == h)
out = dict(passed=all((x['passed'] for x in checks)), checks=len(checks), items=checks, old_curve=old, new_curve=rows, fits=fits, summary=dict(native_ASE_delta_eV_atom=native, native_ASE_delta_meV_atom=1000 * native, explicit_mean_centered_continuum_eV_atom=float(centered), explicit_mean_centered_continuum_meV_atom=float(centered * 1000), integration_range_A3_atom=[lo, hi]), whole_group_complete=False, ASE_native_minima_aligned_100point_RMS_and_mean_centered_continuum_distinguished=True, no_all_element_DeltaCodesDFT_or_convergence_claim=True)
(R / 'VALIDATION.json').write_text(json.dumps(out, indent=2) + '\n')
print('418', out['passed'], out['checks'], out['summary'], flush=True)
assert out['passed']