短时 AIMD

from pathlib import Path
import re,csv,json,sys
import numpy as np
r=Path(__file__).resolve().parent
BOHR=0.529177210903;RYEV=13.605693122994;RYTIME_FS=0.04837768653
scenarios=[('nvt-dt20-cg',20.,100),('nve-dt20-nosym',20.,50),('nve-dt10-nosym',10.,100)]
if '--nve-only' in sys.argv: scenarios=scenarios[1:]
all_summary={};all_rows={};all_xyz={}
for name,dt,nstep in scenarios:
    d=r/name;txt=(d/'al.md.out').read_text();inp=(d/'al.md.in').read_text()
    assert 'JOB DONE.' in txt and 'The maximum number of steps has been reached.' in txt
    assert (d/'al.md.err').stat().st_size==0
    assert not re.search(r'convergence NOT achieved|Error in routine',txt)
    starts=list(re.finditer(r'Entering Dynamics:\s+iteration\s*=\s*(\d+)',txt));assert len(starts)==nstep
    cellblock=inp.split('CELL_PARAMETERS angstrom')[1].split('K_POINTS')[0]
    cell=np.array([list(map(float,l.split())) for l in cellblock.strip().splitlines()])
    init=np.array([list(map(float,l.split()[1:4])) for l in inp.split('ATOMIC_POSITIONS crystal')[1].split('CELL_PARAMETERS')[0].strip().splitlines()])
    assert init.shape==(8,3) and cell.shape==(3,3)
    nconv=len(re.findall(r'convergence has been achieved in\s+(\d+) iterations',txt));assert nconv==nstep
    rows=[];positions=[init@cell]
    for i,m in enumerate(starts):
        step=int(m.group(1));assert step==i+1
        pre=txt[(starts[i-1].end() if i else 0):m.start()]
        end=starts[i+1].start() if i+1<len(starts) else len(txt)
        b=txt[m.start():end]
        potential=float(re.findall(r'!\s+total energy\s*=\s*([-\d.]+)\s+Ry',pre)[-1])
        kinetic=float(re.search(r'kinetic energy \(Ekin\)\s*=\s*([-\d.]+)',b).group(1))
        total=float(re.search(r'Ekin \+ Etot \(const\)\s*=\s*([-\d.]+)',b).group(1))
        temperature=float(re.search(r'temperature\s*=\s*([-\d.]+)\s*K',b).group(1))
        assert abs(total-(potential+kinetic))<2e-8
        assert 'eigenvalues not converged' not in re.split(r'iteration\s+#\s*\d+',pre)[-1]
        early_warnings=pre.count('eigenvalues not converged')
        conv=int(re.findall(r'convergence has been achieved in\s+(\d+) iterations',pre)[-1])
        residual=float(re.findall(r'estimated scf accuracy\s*<\s*([.\dEe+\-]+)\s*Ry',pre)[-1])
        assert residual<=1.0e-10
        posblock=b.split('ATOMIC_POSITIONS (crystal)')[1].strip().splitlines()[:8]
        frac=np.array([list(map(float,l.split()[1:4])) for l in posblock]);assert frac.shape==(8,3)
        positions.append(frac@cell)
        # Position-Verlet in QE 7.5 advances coordinates before output_tau.
        # Etot and the centred finite-difference velocity belong to the preceding position.
        rows.append(dict(step=step,energy_sample_time_fs=(step-1)*dt*RYTIME_FS,position_time_fs=step*dt*RYTIME_FS,temperature_K=temperature,potential_Ry=potential,kinetic_Ry=kinetic,total_Ry=total,scf_iterations=conv,early_diagonalization_warning_count=early_warnings,scf_estimated_accuracy_Ry=residual))
    t=np.array([v['energy_sample_time_fs'] for v in rows]);en=np.array([v['total_Ry'] for v in rows]);temps=np.array([v['temperature_K'] for v in rows]);de=(en-en[0])*RYEV*1000/8
    for v,x in zip(rows,de):v['total_change_meV_atom']=float(x)
    with (d/'thermo.csv').open('w') as f:
        w=csv.DictWriter(f,fieldnames=rows[0].keys());w.writeheader();w.writerows(rows)
    xyz=np.array(positions);np.savez_compressed(d/'trajectory.npz',positions_A=xyz,cell_A=cell,time_fs=np.arange(nstep+1)*dt*RYTIME_FS)
    lattice=' '.join(f'{v:.12f}' for v in cell.ravel())
    with (d/'trajectory.xyz').open('w') as f:
        for j,pos in enumerate(xyz):
            f.write('8\nLattice="'+lattice+'" Properties=species:S:1:pos:R:3 pbc="T T T" time_fs='+str(j*dt*RYTIME_FS)+'\n')
            f.writelines('Al '+' '.join(f'{x:.10f}' for x in p)+'\n' for p in pos)
    wall=re.findall(r'PWSCF\s*:\s*(.*?)\s+CPU\s+(.*?)\s+WALL',txt)[-1][1]
    summary=dict(nsteps=nstep,natoms=8,dt_Ry_au=dt,dt_fs=dt*RYTIME_FS,coordinate_end_time_fs=nstep*dt*RYTIME_FS,energy_end_time_fs=float(t[-1]),all_scf_converged=True,final_scf_iteration_diagonalization_warnings=0,early_scf_diagonalization_warnings=sum(v['early_diagonalization_warning_count'] for v in rows),scf_cycles=nconv,scf_iterations_min=min(x['scf_iterations'] for x in rows),scf_iterations_max=max(x['scf_iterations'] for x in rows),temperature_first_K=float(temps[0]),temperature_mean_K=float(np.mean(temps)),temperature_min_K=float(np.min(temps)),temperature_max_K=float(np.max(temps)),temperature_last_K=float(temps[-1]),energy_change_final_meV_atom=float(de[-1]),energy_peak_to_peak_meV_atom=float(np.ptp(de)),energy_linear_slope_meV_atom_ps=float(np.polyfit(t,de,1)[0]*1000),rms_displacement_final_A=float(np.sqrt(np.mean(np.sum((xyz[-1]-xyz[0])**2,axis=1)))),max_atom_displacement_A=float(np.max(np.linalg.norm(xyz[-1]-xyz[0],axis=1))),wall_time=wall,interpretation='short fixed-volume small-cell trajectory; not thermal-stability or equilibrium-property evidence')
    (d/'summary.json').write_text(json.dumps(summary,indent=2)+'\n');all_summary[name]=summary;all_rows[name]=rows;all_xyz[name]=xyz
coarse=np.array([v['total_change_meV_atom'] for v in all_rows['nve-dt20-nosym']]);fine=np.array([v['total_change_meV_atom'] for v in all_rows['nve-dt10-nosym'][::2]])
assert coarse.shape==fine.shape
matched=dict(matched_energy_times=50,energy_time_end_fs=all_rows['nve-dt20-nosym'][-1]['energy_sample_time_fs'],max_energy_change_difference_meV_atom=float(np.max(np.abs(coarse-fine))),final_position_same_time_fs=all_summary['nve-dt20-nosym']['coordinate_end_time_fs'],rms_position_difference_at_same_end_A=float(np.sqrt(np.mean(np.sum((all_xyz['nve-dt20-nosym'][-1]-all_xyz['nve-dt10-nosym'][-1])**2,axis=1)))))
all_summary['nve_matched_time_comparison']=matched
all_summary['time_axis_note']='QE7.5 position Verlet prints advanced coordinates at n*dt; Etot and centred velocity in that block refer to preceding geometry at (n-1)*dt. CSV keeps both clocks.'
(r/'summary.json').write_text(json.dumps(all_summary,indent=2)+'\n')
print('case                steps  dt_fs    coord_end_fs  energy_range_meV_atom  mean_T_K')
for name,s in all_summary.items():
    if not isinstance(s,dict) or 'nsteps' not in s:continue
    print(f"{name:20s} {s['nsteps']:3d}  {s['dt_fs']:.6f}  {s['coordinate_end_time_fs']:9.6f}  {s['energy_peak_to_peak_meV_atom']:13.6f}       {s['temperature_mean_K']:.3f}")
print(f"All {sum(n for _,_,n in scenarios)} SCF cycles converged; no final-iteration diagonalization warning; {len(scenarios)} native JOB DONE endings.")
print('Early SCF diagonalization warnings:', {n: all_summary[n]['early_scf_diagonalization_warnings'] for n,_,_ in scenarios})
print('NVE position difference at the common final time:',matched['rms_position_difference_at_same_end_A'],'angstrom RMS')
print('NVE total-energy conservation tested only over ~48 fs; no thermodynamic convergence claim.')

# Use energy_sample_time_fs, not position_time_fs.
if (!exists("data_root")) data_root = "."
if (!exists("out_root")) out_root = "."
coarse = data_root."/aimd/nve-dt20-nosym/thermo.csv"
fine = data_root."/aimd/nve-dt10-nosym/thermo.csv"
set encoding utf8
set datafile separator ","
set datafile columnheaders
stats fine using ($1==1 ? $5 : 1/0) nooutput
p0 = STATS_mean
stats fine using ($1==1 ? $6 : 1/0) nooutput
k0 = STATS_mean
conv = 13.605693122994*1000/8
set border 3
set tics out nomirror
set key top left
do for [ext in "svg png pdf"] {
    if (ext eq "svg") { set terminal svg size 1000,400 enhanced font "DejaVu Sans,11" }
    if (ext eq "png") { set terminal pngcairo size 1000,400 enhanced font "DejaVu Sans,11" }
    if (ext eq "pdf") { set terminal pdfcairo size 10,4 enhanced font "DejaVu Sans,11" }
    set output out_root."/al-nve-gnuplot.".ext
    set size 1,1
    set origin 0,0
    unset title
    set multiplot layout 1,2 margins 0.11,0.98,0.18,0.82 spacing 0.14,0.05 title "8-atom Al | existing short NVE records"
    set title "(a) 50 matched energy times"
    set xlabel "Energy sample time (fs)"
    set ylabel "Change of total energy (meV/atom)"
    set xrange [0:47.4101328]
    set yrange [-0.002:0.022]
    plot coarse using 2:11 with linespoints pt 7 ps 0.3 lw 1 lc rgb "#D55E00" title "dt = 0.967554 fs", fine every 2 using 2:11 with linespoints pt 5 ps 0.3 lw 1 lc rgb "#CC79A7" title "dt = 0.483777 fs"
    set title "(b) Energy exchange, fine step"
    set ylabel "Change from first sample (meV/atom)"
    set xrange [0:47.8939097]
    set yrange [*:*]
    plot fine using 2:(($5-p0)*conv) with lines lw 1 lc rgb "#009E73" title "Potential", fine using 2:(($6-k0)*conv) with lines lw 1 lc rgb "#D55E00" title "Kinetic"
    unset multiplot
    unset output
}
print "Wrote al-nve-gnuplot.svg/.png/.pdf; panel (a) uses 50 shared energy times"

"""Inspect existing Al trajectory frames and match their native energy clock."""
from pathlib import Path
import csv, itertools, json, sys
import numpy as np

root = Path(sys.argv[1]) if len(sys.argv) > 1 else Path(".")
out = Path(sys.argv[2]) if len(sys.argv) > 2 else Path("structure-observations")
out.mkdir(parents=True, exist_ok=True)
cases = ("nve-dt20-nosym", "nve-dt10-nosym")
summary = {}
end_frames = {}
for case in cases:
    folder = root / "aimd" / case
    thermo = list(csv.DictReader((folder / "thermo.csv").open()))
    with np.load(folder / "trajectory.npz") as trajectory:
        positions = trajectory["positions_A"]
        cell = trajectory["cell_A"]
        time = trajectory["time_fs"]
    assert positions.shape == (len(thermo) + 1, 8, 3)
    assert cell.shape == (3, 3) and np.isfinite(positions).all()
    assert all(abs(float(row["energy_sample_time_fs"]) - time[j]) < 1e-8 for j, row in enumerate(thermo))
    # Use explicit periodic images; component-wise rounding is unsuitable for this oblique cell.
    candidates = []
    for i, j, shift in itertools.product(range(8), range(8), itertools.product((-1, 0, 1), repeat=3)):
        if i == j and shift == (0, 0, 0):
            continue
        vector = positions[0, j] - positions[0, i] + np.array(shift) @ cell
        candidates.append((float(np.linalg.norm(vector)), i, j, shift))
    nearest = min(item[0] for item in candidates)
    bonds = []
    for distance, i, j, shift in candidates:
        reverse = (j, i, tuple(-n for n in shift))
        if abs(distance - nearest) < 1e-7 and (i, j, shift) < reverse:
            bonds.append((i, j, np.array(shift) @ cell))
    assert len(bonds) == 48
    displacement = positions - positions[0]
    rms = np.sqrt(np.mean(np.sum(displacement**2, axis=2), axis=1))
    lengths = np.array([[np.linalg.norm(frame[j] - frame[i] + translation)
                         for i, j, translation in bonds] for frame in positions])
    rows = []
    for frame, t in enumerate(time):
        rows.append(dict(time_fs=float(t),rms_A=float(rms[frame]),
                         bond_min_A=float(lengths[frame].min()),bond_max_A=float(lengths[frame].max()),
                         bond_mean_A=float(lengths[frame].mean()),
                         temperature_K=float(thermo[frame]["temperature_K"]) if frame < len(thermo) else ""))
    with (out / (case + "-structure.csv")).open("w") as file:
        writer = csv.DictWriter(file, fieldnames=rows[0].keys())
        writer.writeheader()
        writer.writerows(rows)
    cold = min(range(len(thermo)), key=lambda j: float(thermo[j]["temperature_K"]))
    e0 = thermo[0]
    row = thermo[cold]
    factor = 13.605693122994 * 1000 / 8
    summary[case] = dict(nframes=len(time),initial_neighbor_pairs=48,initial_neighbor_A=nearest,
                         rms_end_A=float(rms[-1]),rms_peak_A=float(rms.max()),
                         rms_peak_time_fs=float(time[rms.argmax()]),
                         end_bond_min_A=float(lengths[-1].min()),end_bond_max_A=float(lengths[-1].max()),
                         end_bond_mean_A=float(lengths[-1].mean()),
                         cold_energy_time_fs=float(row["energy_sample_time_fs"]),
                         cold_temperature_K=float(row["temperature_K"]),cold_coordinate_frame=cold,
                         cold_rms_A=float(rms[cold]),cold_bond_min_A=float(lengths[cold].min()),
                         cold_bond_max_A=float(lengths[cold].max()),
                         cold_potential_change_meV_atom=(float(row["potential_Ry"])-float(e0["potential_Ry"]))*factor,
                         cold_kinetic_change_meV_atom=(float(row["kinetic_Ry"])-float(e0["kinetic_Ry"]))*factor,
                         cold_total_change_meV_atom=(float(row["total_Ry"])-float(e0["total_Ry"]))*factor)
    end_frames[case] = positions[-1]
    print(f"{case}: frames={len(time)} initial_pairs=48 initial_nn={nearest:.8f} A")
    print(f"  RMS_end={rms[-1]:.8f} A RMS_peak={rms.max():.8f} A at {time[rms.argmax()]:.8f} fs")
    print(f"  end_initial_neighbor_range={lengths[-1].min():.8f}..{lengths[-1].max():.8f} A")
difference = end_frames[cases[0]] - end_frames[cases[1]]
summary["matched_end_rms_difference_A"] = float(np.sqrt(np.mean(np.sum(difference**2, axis=1))))
summary["definition"] = "48 unique periodic-image pairs in the initial first neighbor shell, tracked without changing pair identity; range is not a confidence interval or a coordination-number test."
summary["time_axis"] = "Coordinate frame j at j*dt matches thermo energy row j at j*dt; the final coordinate frame has no matching printed energy/temperature sample."
(out / "structure-summary.json").write_text(json.dumps(summary, indent=2) + "\n")
print(f"matched_end_RMS_difference={summary['matched_end_rms_difference_A']:.10e} A")

# Existing Al coordinate frames; the pair interval is a geometric range, not uncertainty.
if (!exists("data_root")) data_root = "structure-observations"
if (!exists("out_root")) out_root = data_root
coarse = data_root."/nve-dt20-nosym-structure.csv"
fine = data_root."/nve-dt10-nosym-structure.csv"
set encoding utf8
set datafile separator ","
set datafile columnheaders
stats fine using ($1==0 ? $5 : 1/0) nooutput
r0 = STATS_mean
set border 3
set tics out nomirror
do for [ext in "svg png pdf"] {
    if (ext eq "svg") { set terminal svg size 1100,430 enhanced font "DejaVu Sans,11" }
    if (ext eq "png") { set terminal pngcairo size 1100,430 enhanced font "DejaVu Sans,11" }
    if (ext eq "pdf") { set terminal pdfcairo size 11,4.3 enhanced font "DejaVu Sans,11" }
    set output out_root."/al-nve-structure.".ext
    set size 1,1
    set origin 0,0
    unset title
    set multiplot layout 1,2 margins 0.09,0.98,0.18,0.83 spacing 0.14,0.05 title "8-atom Al | recorded positions, fixed cell"
    set title "(a) Motion from the shared initial frame"
    set xlabel "Coordinate time (fs)"
    set ylabel "RMS displacement (Å)"
    set xrange [0:48.37768653]
    set yrange [0:0.11]
    set key top left
    plot coarse using 1:2 with lines lw 1.5 lc rgb "#D55E00" title "dt = 0.967554 fs", fine using 1:2 with lines dt 2 lw 1.5 lc rgb "#CC79A7" title "dt = 0.483777 fs"
    set title "(b) Initial neighbor pairs, fine step"
    set ylabel "Al-Al pair length (Å)"
    set yrange [*:*]
    set key top left
    set style fill transparent solid 0.2 noborder
    plot fine using 1:3:4 with filledcurves lc rgb "#009E73" title "min-max of 48 pairs", fine using 1:5 with lines lw 1.5 lc rgb "#009E73" title "mean", r0 with lines dt 3 lc rgb "#777777" title "initial distance"
    unset multiplot
    unset output
}
print "Wrote al-nve-structure.svg/.png/.pdf from all recorded coordinate frames"

"""Plot the verified QE Al AIMD records from the Al bundle root."""

from atlas_plot_style import install as install_atlas_style
install_atlas_style()
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
r=Path(__file__).resolve().parent;d=r/'aimd';out=r/'figures';out.mkdir(exist_ok=True)
plt.rcParams.update({'font.size':10,'axes.spines.top':False,'axes.spines.right':False,'savefig.dpi':220})
colors=['#009e73','#d55e00','#cc79a7']
def data(name):return np.genfromtxt(d/name/'thermo.csv',delimiter=',',names=True)
def save(fig,name):
    fig.tight_layout();fig.savefig(out/(name+'.png'),bbox_inches='tight');fig.savefig(out/(name+'.pdf'),bbox_inches='tight');plt.close(fig)
nvt=data('nvt-dt20-cg');nve20=data('nve-dt20-nosym');nve10=data('nve-dt10-nosym')
fig,ax=plt.subplots(2,1,figsize=(7.4,6),sharex=True)
ax[0].plot(nvt['energy_sample_time_fs'],nvt['temperature_K'],color=colors[0],label='SVR target 300 K; dt = 0.968 fs')
ax[0].axhline(300,color='#777',ls='--',lw=1,label='Target')
ax[0].set(ylabel='Instantaneous temperature (K)');ax[0].legend(frameon=False,fontsize=9)
for arr,c,label in [(nve20,colors[1],'NVE dt = 0.968 fs'),(nve10,colors[2],'NVE dt = 0.484 fs')]:
    ax[1].plot(arr['energy_sample_time_fs'],arr['temperature_K'],color=c,label=label)
ax[1].set(xlabel='Time of the sampled energy/velocity (fs)',ylabel='Instantaneous temperature (K)');ax[1].legend(frameon=False,fontsize=9)
fig.suptitle('8-atom periodic Al | short trajectories, not a thermal-stability test',fontsize=11);save(fig,'aimd-temperature')
fig,ax=plt.subplots(1,2,figsize=(10,4))
for arr,c,label in [(nve20,colors[1],'dt = 0.968 fs'),(nve10,colors[2],'dt = 0.484 fs')]:
    ax[0].plot(arr['energy_sample_time_fs'],arr['total_change_meV_atom'],color=c,label=label)
ax[0].set(xlabel='Time (fs)',ylabel='Change of total energy (meV/atom)');ax[0].legend(frameon=False)
conv=13.605693122994*1000/8
for field,c,label in [('kinetic_Ry',colors[1],'Kinetic'),('potential_Ry',colors[0],'Potential')]:
    ax[1].plot(nve10['energy_sample_time_fs'],(nve10[field]-nve10[field][0])*conv,color=c,label=label)
ax[1].set(xlabel='Time (fs)',ylabel='Energy change (meV/atom)');ax[1].legend(frameon=False)
fig.suptitle('NVE step-size check over ~48 fs | identical initial positions and velocities',fontsize=11);save(fig,'aimd-energy')
fig,ax=plt.subplots(figsize=(7.2,3.7))
for name,c,label in [('nvt-dt20-cg',colors[0],'SVR'),('nve-dt20-nosym',colors[1],'NVE dt = 0.968 fs'),('nve-dt10-nosym',colors[2],'NVE dt = 0.484 fs')]:
    a=np.load(d/name/'trajectory.npz');p=a['positions_A'];rms=np.sqrt(np.mean(np.sum((p-p[0])**2,axis=2),axis=1))
    ax.plot(a['time_fs'],rms,color=c,label=label)
ax.set(xlabel='Coordinate time (fs)',ylabel='RMS displacement from initial positions (Å)',title='Short-time displacement; no diffusion or long-term stability inference')
ax.legend(frameon=False);save(fig,'aimd-displacement')
print('Wrote aimd-temperature, aimd-energy, aimd-displacement as PNG and PDF')

"""Export existing Al NVE frames with their full periodic cells."""
from pathlib import Path
import sys
import numpy as np
from ase.io import read, write

root = Path(sys.argv[1]) if len(sys.argv) > 1 else Path(".")
out = Path(sys.argv[2]) if len(sys.argv) > 2 else Path("nve-frames")
out.mkdir(parents=True, exist_ok=True)
coarse = root / "aimd/nve-dt20-nosym/trajectory.xyz"
fine = root / "aimd/nve-dt10-nosym/trajectory.xyz"
initial = read(coarse, index=0, format="extxyz")
initial_fine = read(fine, index=0, format="extxyz")
assert np.allclose(initial.positions, initial_fine.positions, atol=1e-10, rtol=0)
assert np.allclose(initial.cell, initial_fine.cell, atol=1e-10, rtol=0)
frames = [("initial.POSCAR", initial),
          ("nve-dt20-end.POSCAR", read(coarse, index=-1, format="extxyz")),
          ("nve-dt10-end.POSCAR", read(fine, index=-1, format="extxyz"))]
for name, atoms in frames:
    assert len(atoms) == 8 and atoms.get_chemical_symbols() == ["Al"] * 8
    assert atoms.pbc.all() and np.allclose(atoms.cell, initial.cell, atol=1e-10, rtol=0)
    if name != "initial.POSCAR":
        assert abs(float(atoms.info["time_fs"]) - 48.37768653) < 1e-8
    path = out / name
    write(path, atoms, format="vasp", direct=True, vasp5=True, sort=False)
    restored = read(path, format="vasp")
    assert np.allclose(restored.cell, atoms.cell, atol=1e-10, rtol=0)
    assert np.allclose(restored.positions, atoms.positions, atol=1e-10, rtol=0)
    print(name, "coordinate_time_fs =", atoms.info["time_fs"])