From 6d5d130f8038bb10555409dbcd6dfd727a6e202d Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Mon, 1 Jun 2026 13:43:28 +0800 Subject: [PATCH 1/4] Fix: change kF_r type from complex to double in cal_CD_potential kF_r is always real-valued (derived from std::abs), so using std::complex wastes half the memory. Co-Authored-By: Claude Opus 4.8 --- source/source_pw/module_ofdft/evolve_ofdft.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/source_pw/module_ofdft/evolve_ofdft.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp index 4157dc6516a..51e2fe4076f 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.cpp +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -183,7 +183,7 @@ void Evolve_OFDFT::cal_CD_potential(std::vector>& psi_, std::vector> rCurrent_x(nrxx); std::vector> rCurrent_y(nrxx); std::vector> rCurrent_z(nrxx); - std::vector> kF_r(nrxx); + std::vector kF_r(nrxx); std::vector> rCDPotential(nrxx); for (int ir = 0; ir < nrxx; ++ir) From 178eb2f5b59bee070498fc5292f318ddab03cb8a Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Mon, 1 Jun 2026 13:43:59 +0800 Subject: [PATCH 2/4] Fix: guard against division by zero in cal_CD_potential When kF_r is zero (wavefunction is zero at a grid point), the denominator 2*kF^2 becomes zero, producing inf which std::isnan() does not catch. Also, the previous approach zeroed the entire rpot when NaN was detected, discarding Hartree/XC/TF contributions. Now skip the CD term when kF_r is too small, avoiding both issues. Co-Authored-By: Claude Opus 4.8 --- source/source_pw/module_ofdft/evolve_ofdft.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/source/source_pw/module_ofdft/evolve_ofdft.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp index 51e2fe4076f..7f276781ddb 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.cpp +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -226,11 +226,10 @@ void Evolve_OFDFT::cal_CD_potential(std::vector>& psi_, for (int ir = 0; ir < nrxx; ++ir) { - rpot(is, ir) -= mCD_para*2.0*std::real(rCDPotential[ir])*std::pow(ModuleBase::PI,3) - / (2.0*std::pow(std::real(kF_r[ir]),2)); - if (std::isnan(rpot(is, ir))) + if (kF_r[ir] > 1e-12) { - rpot(is, ir)=0.0; + rpot(is, ir) -= mCD_para*2.0*std::real(rCDPotential[ir])*std::pow(ModuleBase::PI,3) + / (2.0*kF_r[ir]*kF_r[ir]); } } } From 9097e5602d2246e84becc4ee005c80116d2908c6 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Fri, 19 Jun 2026 11:46:44 +0800 Subject: [PATCH 3/4] Add TDOFDFT example and kinetic energy post-processing tool - Add examples/35_tdofdft/ with Al_proton ESP calculation (01_gs: ground-state OFDFT, 02_td: TDOFDT proton dynamics) - Add tools/02_postprocessing/tdofdft-tools/get_kinetic.py to extract proton kinetic energy curve from MD_dump Co-Authored-By: Claude --- .../35_tdofdft/01_Al_proton_ESP/01_gs/INPUT | 29 ++ .../35_tdofdft/01_Al_proton_ESP/01_gs/KPT | 4 + .../35_tdofdft/01_Al_proton_ESP/01_gs/STRU | 273 +++++++++++++++++ .../35_tdofdft/01_Al_proton_ESP/01_gs/run.sh | 18 ++ .../35_tdofdft/01_Al_proton_ESP/02_td/INPUT | 43 +++ .../35_tdofdft/01_Al_proton_ESP/02_td/KPT | 4 + .../35_tdofdft/01_Al_proton_ESP/02_td/STRU | 279 ++++++++++++++++++ .../35_tdofdft/01_Al_proton_ESP/02_td/run.sh | 18 ++ examples/35_tdofdft/README | 90 ++++++ .../tdofdft-tools/get_kinetic.py | 79 +++++ 10 files changed, 837 insertions(+) create mode 100644 examples/35_tdofdft/01_Al_proton_ESP/01_gs/INPUT create mode 100644 examples/35_tdofdft/01_Al_proton_ESP/01_gs/KPT create mode 100644 examples/35_tdofdft/01_Al_proton_ESP/01_gs/STRU create mode 100644 examples/35_tdofdft/01_Al_proton_ESP/01_gs/run.sh create mode 100644 examples/35_tdofdft/01_Al_proton_ESP/02_td/INPUT create mode 100644 examples/35_tdofdft/01_Al_proton_ESP/02_td/KPT create mode 100644 examples/35_tdofdft/01_Al_proton_ESP/02_td/STRU create mode 100644 examples/35_tdofdft/01_Al_proton_ESP/02_td/run.sh create mode 100644 examples/35_tdofdft/README create mode 100644 tools/02_postprocessing/tdofdft-tools/get_kinetic.py diff --git a/examples/35_tdofdft/01_Al_proton_ESP/01_gs/INPUT b/examples/35_tdofdft/01_Al_proton_ESP/01_gs/INPUT new file mode 100644 index 00000000000..542846c5edb --- /dev/null +++ b/examples/35_tdofdft/01_Al_proton_ESP/01_gs/INPUT @@ -0,0 +1,29 @@ +INPUT_PARAMETERS +#Parameters (1.General) +suffix autotest +calculation scf +esolver_type ofdft + +symmetry 1 +pseudo_dir ../../../tests/PP_ORB +pseudo_rcut 16 +nspin 1 +cal_force 1 +test_force 1 +cal_stress 1 +test_stress 1 + +#Parameters (2.Iteration) +ecutwfc 44 +scf_nmax 50 + +#OFDFT +of_kinetic tf+ +of_method tn +of_conv energy +of_tole 2e-6 + +#Parameters (3.Basis) +basis_type pw + +out_chg 1 diff --git a/examples/35_tdofdft/01_Al_proton_ESP/01_gs/KPT b/examples/35_tdofdft/01_Al_proton_ESP/01_gs/KPT new file mode 100644 index 00000000000..c289c0158aa --- /dev/null +++ b/examples/35_tdofdft/01_Al_proton_ESP/01_gs/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/examples/35_tdofdft/01_Al_proton_ESP/01_gs/STRU b/examples/35_tdofdft/01_Al_proton_ESP/01_gs/STRU new file mode 100644 index 00000000000..fc3340f9f57 --- /dev/null +++ b/examples/35_tdofdft/01_Al_proton_ESP/01_gs/STRU @@ -0,0 +1,273 @@ +ATOMIC_SPECIES +Al 26.981539 al.lda.lps blps + +LATTICE_CONSTANT +30.613575 // add lattice constant + +LATTICE_VECTORS +1.000000000000 0.000000000000 0.000000000000 +0.000000000000 1.000000000000 0.000000000000 +0.000000000000 0.000000000000 1.000000000000 + +ATOMIC_POSITIONS +Direct + +Al +0 +256 +0 0 0 0 0 0 +0 0.125 0.125 0 0 0 +0.125 0 0.125 0 0 0 +0.125 0.125 0 0 0 0 +0 0 0.25 0 0 0 +0 0.125 0.375 0 0 0 +0.125 0 0.375 0 0 0 +0.125 0.125 0.25 0 0 0 +0 0 0.5 0 0 0 +0 0.125 0.625 0 0 0 +0.125 0 0.625 0 0 0 +0.125 0.125 0.5 0 0 0 +0 0 0.75 0 0 0 +0 0.125 0.875 0 0 0 +0.125 0 0.875 0 0 0 +0.125 0.125 0.75 0 0 0 +0 0.25 0 0 0 0 +0 0.375 0.125 0 0 0 +0.125 0.25 0.125 0 0 0 +0.125 0.375 0 0 0 0 +0 0.25 0.25 0 0 0 +0 0.375 0.375 0 0 0 +0.125 0.25 0.375 0 0 0 +0.125 0.375 0.25 0 0 0 +0 0.25 0.5 0 0 0 +0 0.375 0.625 0 0 0 +0.125 0.25 0.625 0 0 0 +0.125 0.375 0.5 0 0 0 +0 0.25 0.75 0 0 0 +0 0.375 0.875 0 0 0 +0.125 0.25 0.875 0 0 0 +0.125 0.375 0.75 0 0 0 +0.25 0 0 0 0 0 +0.25 0.125 0.125 0 0 0 +0.375 0 0.125 0 0 0 +0.375 0.125 0 0 0 0 +0.25 0 0.25 0 0 0 +0.25 0.125 0.375 0 0 0 +0.375 0 0.375 0 0 0 +0.375 0.125 0.25 0 0 0 +0.25 0 0.5 0 0 0 +0.25 0.125 0.625 0 0 0 +0.375 0 0.625 0 0 0 +0.375 0.125 0.5 0 0 0 +0.25 0 0.75 0 0 0 +0.25 0.125 0.875 0 0 0 +0.375 0 0.875 0 0 0 +0.375 0.125 0.75 0 0 0 +0.25 0.25 0 0 0 0 +0.25 0.375 0.125 0 0 0 +0.375 0.25 0.125 0 0 0 +0.375 0.375 0 0 0 0 +0.25 0.25 0.25 0 0 0 +0.25 0.375 0.375 0 0 0 +0.375 0.25 0.375 0 0 0 +0.375 0.375 0.25 0 0 0 +0.25 0.25 0.5 0 0 0 +0.25 0.375 0.625 0 0 0 +0.375 0.25 0.625 0 0 0 +0.375 0.375 0.5 0 0 0 +0.25 0.25 0.75 0 0 0 +0.25 0.375 0.875 0 0 0 +0.375 0.25 0.875 0 0 0 +0.375 0.375 0.75 0 0 0 +0.5 0 0 0 0 0 +0.5 0.125 0.125 0 0 0 +0.625 0 0.125 0 0 0 +0.625 0.125 0 0 0 0 +0.5 0 0.25 0 0 0 +0.5 0.125 0.375 0 0 0 +0.625 0 0.375 0 0 0 +0.625 0.125 0.25 0 0 0 +0.5 0 0.5 0 0 0 +0.5 0.125 0.625 0 0 0 +0.625 0 0.625 0 0 0 +0.625 0.125 0.5 0 0 0 +0.5 0 0.75 0 0 0 +0.5 0.125 0.875 0 0 0 +0.625 0 0.875 0 0 0 +0.625 0.125 0.75 0 0 0 +0.5 0.25 0 0 0 0 +0.5 0.375 0.125 0 0 0 +0.625 0.25 0.125 0 0 0 +0.625 0.375 0 0 0 0 +0.5 0.25 0.25 0 0 0 +0.5 0.375 0.375 0 0 0 +0.625 0.25 0.375 0 0 0 +0.625 0.375 0.25 0 0 0 +0.5 0.25 0.5 0 0 0 +0.5 0.375 0.625 0 0 0 +0.625 0.25 0.625 0 0 0 +0.625 0.375 0.5 0 0 0 +0.5 0.25 0.75 0 0 0 +0.5 0.375 0.875 0 0 0 +0.625 0.25 0.875 0 0 0 +0.625 0.375 0.75 0 0 0 +0.75 0 0 0 0 0 +0.75 0.125 0.125 0 0 0 +0.875 0 0.125 0 0 0 +0.875 0.125 0 0 0 0 +0.75 0 0.25 0 0 0 +0.75 0.125 0.375 0 0 0 +0.875 0 0.375 0 0 0 +0.875 0.125 0.25 0 0 0 +0.75 0 0.5 0 0 0 +0.75 0.125 0.625 0 0 0 +0.875 0 0.625 0 0 0 +0.875 0.125 0.5 0 0 0 +0.75 0 0.75 0 0 0 +0.75 0.125 0.875 0 0 0 +0.875 0 0.875 0 0 0 +0.875 0.125 0.75 0 0 0 +0.75 0.25 0 0 0 0 +0.75 0.375 0.125 0 0 0 +0.875 0.25 0.125 0 0 0 +0.875 0.375 0 0 0 0 +0.75 0.25 0.25 0 0 0 +0.75 0.375 0.375 0 0 0 +0.875 0.25 0.375 0 0 0 +0.875 0.375 0.25 0 0 0 +0.75 0.25 0.5 0 0 0 +0.75 0.375 0.625 0 0 0 +0.875 0.25 0.625 0 0 0 +0.875 0.375 0.5 0 0 0 +0.75 0.25 0.75 0 0 0 +0.75 0.375 0.875 0 0 0 +0.875 0.25 0.875 0 0 0 +0.875 0.375 0.75 0 0 0 +0 0.5 0 0 0 0 +0 0.625 0.125 0 0 0 +0.125 0.5 0.125 0 0 0 +0.125 0.625 0 0 0 0 +0 0.5 0.25 0 0 0 +0 0.625 0.375 0 0 0 +0.125 0.5 0.375 0 0 0 +0.125 0.625 0.25 0 0 0 +0 0.5 0.5 0 0 0 +0 0.625 0.625 0 0 0 +0.125 0.5 0.625 0 0 0 +0.125 0.625 0.5 0 0 0 +0 0.5 0.75 0 0 0 +0 0.625 0.875 0 0 0 +0.125 0.5 0.875 0 0 0 +0.125 0.625 0.75 0 0 0 +0 0.75 0 0 0 0 +0 0.875 0.125 0 0 0 +0.125 0.75 0.125 0 0 0 +0.125 0.875 0 0 0 0 +0 0.75 0.25 0 0 0 +0 0.875 0.375 0 0 0 +0.125 0.75 0.375 0 0 0 +0.125 0.875 0.25 0 0 0 +0 0.75 0.5 0 0 0 +0 0.875 0.625 0 0 0 +0.125 0.75 0.625 0 0 0 +0.125 0.875 0.5 0 0 0 +0 0.75 0.75 0 0 0 +0 0.875 0.875 0 0 0 +0.125 0.75 0.875 0 0 0 +0.125 0.875 0.75 0 0 0 +0.25 0.5 0 0 0 0 +0.25 0.625 0.125 0 0 0 +0.375 0.5 0.125 0 0 0 +0.375 0.625 0 0 0 0 +0.25 0.5 0.25 0 0 0 +0.25 0.625 0.375 0 0 0 +0.375 0.5 0.375 0 0 0 +0.375 0.625 0.25 0 0 0 +0.25 0.5 0.5 0 0 0 +0.25 0.625 0.625 0 0 0 +0.375 0.5 0.625 0 0 0 +0.375 0.625 0.5 0 0 0 +0.25 0.5 0.75 0 0 0 +0.25 0.625 0.875 0 0 0 +0.375 0.5 0.875 0 0 0 +0.375 0.625 0.75 0 0 0 +0.25 0.75 0 0 0 0 +0.25 0.875 0.125 0 0 0 +0.375 0.75 0.125 0 0 0 +0.375 0.875 0 0 0 0 +0.25 0.75 0.25 0 0 0 +0.25 0.875 0.375 0 0 0 +0.375 0.75 0.375 0 0 0 +0.375 0.875 0.25 0 0 0 +0.25 0.75 0.5 0 0 0 +0.25 0.875 0.625 0 0 0 +0.375 0.75 0.625 0 0 0 +0.375 0.875 0.5 0 0 0 +0.25 0.75 0.75 0 0 0 +0.25 0.875 0.875 0 0 0 +0.375 0.75 0.875 0 0 0 +0.375 0.875 0.75 0 0 0 +0.5 0.5 0 0 0 0 +0.5 0.625 0.125 0 0 0 +0.625 0.5 0.125 0 0 0 +0.625 0.625 0 0 0 0 +0.5 0.5 0.25 0 0 0 +0.5 0.625 0.375 0 0 0 +0.625 0.5 0.375 0 0 0 +0.625 0.625 0.25 0 0 0 +0.5 0.5 0.5 0 0 0 +0.5 0.625 0.625 0 0 0 +0.625 0.5 0.625 0 0 0 +0.625 0.625 0.5 0 0 0 +0.5 0.5 0.75 0 0 0 +0.5 0.625 0.875 0 0 0 +0.625 0.5 0.875 0 0 0 +0.625 0.625 0.75 0 0 0 +0.5 0.75 0 0 0 0 +0.5 0.875 0.125 0 0 0 +0.625 0.75 0.125 0 0 0 +0.625 0.875 0 0 0 0 +0.5 0.75 0.25 0 0 0 +0.5 0.875 0.375 0 0 0 +0.625 0.75 0.375 0 0 0 +0.625 0.875 0.25 0 0 0 +0.5 0.75 0.5 0 0 0 +0.5 0.875 0.625 0 0 0 +0.625 0.75 0.625 0 0 0 +0.625 0.875 0.5 0 0 0 +0.5 0.75 0.75 0 0 0 +0.5 0.875 0.875 0 0 0 +0.625 0.75 0.875 0 0 0 +0.625 0.875 0.75 0 0 0 +0.75 0.5 0 0 0 0 +0.75 0.625 0.125 0 0 0 +0.875 0.5 0.125 0 0 0 +0.875 0.625 0 0 0 0 +0.75 0.5 0.25 0 0 0 +0.75 0.625 0.375 0 0 0 +0.875 0.5 0.375 0 0 0 +0.875 0.625 0.25 0 0 0 +0.75 0.5 0.5 0 0 0 +0.75 0.625 0.625 0 0 0 +0.875 0.5 0.625 0 0 0 +0.875 0.625 0.5 0 0 0 +0.75 0.5 0.75 0 0 0 +0.75 0.625 0.875 0 0 0 +0.875 0.5 0.875 0 0 0 +0.875 0.625 0.75 0 0 0 +0.75 0.75 0 0 0 0 +0.75 0.875 0.125 0 0 0 +0.875 0.75 0.125 0 0 0 +0.875 0.875 0 0 0 0 +0.75 0.75 0.25 0 0 0 +0.75 0.875 0.375 0 0 0 +0.875 0.75 0.375 0 0 0 +0.875 0.875 0.25 0 0 0 +0.75 0.75 0.5 0 0 0 +0.75 0.875 0.625 0 0 0 +0.875 0.75 0.625 0 0 0 +0.875 0.875 0.5 0 0 0 +0.75 0.75 0.75 0 0 0 +0.75 0.875 0.875 0 0 0 +0.875 0.75 0.875 0 0 0 +0.875 0.875 0.75 0 0 0 diff --git a/examples/35_tdofdft/01_Al_proton_ESP/01_gs/run.sh b/examples/35_tdofdft/01_Al_proton_ESP/01_gs/run.sh new file mode 100644 index 00000000000..7623ec13913 --- /dev/null +++ b/examples/35_tdofdft/01_Al_proton_ESP/01_gs/run.sh @@ -0,0 +1,18 @@ +#!/bin/bash + +ABACUS_PATH=$(awk -F "=" '$1=="ABACUS_PATH"{print $2}' ../../../../SETENV) +ABACUS_NPROCS=$(awk -F "=" '$1=="ABACUS_NPROCS"{print $2}' ../../../../SETENV) +ABACUS_THREADS=$(awk -F "=" '$1=="ABACUS_THREADS"{print $2}' ../../../../SETENV) + +OMP_NUM_THREADS=${ABACUS_THREADS} mpirun -np ${ABACUS_NPROCS} ${ABACUS_PATH} | tee relax.output + +if [[ ! -f scf.output ]] || + [[ ! -f OUT.ABACUS/running_scf.log ]] || + [[ ! ( "$(tail -1 OUT.ABACUS/running_scf.log)" == " Total Time :"* ) ]] +then + echo "job failed!" + exit 1 +else + echo "job succeeded!" + exit 0 +fi diff --git a/examples/35_tdofdft/01_Al_proton_ESP/02_td/INPUT b/examples/35_tdofdft/01_Al_proton_ESP/02_td/INPUT new file mode 100644 index 00000000000..1a33f282846 --- /dev/null +++ b/examples/35_tdofdft/01_Al_proton_ESP/02_td/INPUT @@ -0,0 +1,43 @@ +INPUT_PARAMETERS +#Parameters (1.General) +suffix autotest +calculation md +esolver_type tdofdft + +symmetry 1 +pseudo_dir ../../../tests/PP_ORB +pseudo_rcut 16 +nspin 1 +cal_force 1 +test_force 1 +cal_stress 1 +test_stress 1 + +#Parameters (2.Iteration) +ecutwfc 44 +scf_nmax 50 + +#OFDFT +of_kinetic tf+ +of_method tn +of_conv energy +of_tole 2e-6 +of_cd 1 +of_mcd_alpha 2.0 + +#Parameters (3.Basis) +basis_type pw + +init_vel 1 +init_chg file +read_file_dir ./ + +md_restart 0 +md_type nve +md_nstep 4500 +md_dt 0.0005 +md_dumpfreq 1 + +# Total charge: 256 Al * 3 valence electrons = 768 +# (proton contributes +1 charge, no valence electrons) +nelec 768 diff --git a/examples/35_tdofdft/01_Al_proton_ESP/02_td/KPT b/examples/35_tdofdft/01_Al_proton_ESP/02_td/KPT new file mode 100644 index 00000000000..c289c0158aa --- /dev/null +++ b/examples/35_tdofdft/01_Al_proton_ESP/02_td/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/examples/35_tdofdft/01_Al_proton_ESP/02_td/STRU b/examples/35_tdofdft/01_Al_proton_ESP/02_td/STRU new file mode 100644 index 00000000000..ed7e153c324 --- /dev/null +++ b/examples/35_tdofdft/01_Al_proton_ESP/02_td/STRU @@ -0,0 +1,279 @@ +ATOMIC_SPECIES +Al 26.981539 al.lda.lps blps +H 1.008 H.pz-vbc.UPF 1/r + +LATTICE_CONSTANT +30.613575 // add lattice constant + +LATTICE_VECTORS +1.000000000000 0.000000000000 0.000000000000 +0.000000000000 1.000000000000 0.000000000000 +0.000000000000 0.000000000000 1.000000000000 + +ATOMIC_POSITIONS +Direct + +Al +0 +256 +0 0 0 0 0 0 +0 0.125 0.125 0 0 0 +0.125 0 0.125 0 0 0 +0.125 0.125 0 0 0 0 +0 0 0.25 0 0 0 +0 0.125 0.375 0 0 0 +0.125 0 0.375 0 0 0 +0.125 0.125 0.25 0 0 0 +0 0 0.5 0 0 0 +0 0.125 0.625 0 0 0 +0.125 0 0.625 0 0 0 +0.125 0.125 0.5 0 0 0 +0 0 0.75 0 0 0 +0 0.125 0.875 0 0 0 +0.125 0 0.875 0 0 0 +0.125 0.125 0.75 0 0 0 +0 0.25 0 0 0 0 +0 0.375 0.125 0 0 0 +0.125 0.25 0.125 0 0 0 +0.125 0.375 0 0 0 0 +0 0.25 0.25 0 0 0 +0 0.375 0.375 0 0 0 +0.125 0.25 0.375 0 0 0 +0.125 0.375 0.25 0 0 0 +0 0.25 0.5 0 0 0 +0 0.375 0.625 0 0 0 +0.125 0.25 0.625 0 0 0 +0.125 0.375 0.5 0 0 0 +0 0.25 0.75 0 0 0 +0 0.375 0.875 0 0 0 +0.125 0.25 0.875 0 0 0 +0.125 0.375 0.75 0 0 0 +0.25 0 0 0 0 0 +0.25 0.125 0.125 0 0 0 +0.375 0 0.125 0 0 0 +0.375 0.125 0 0 0 0 +0.25 0 0.25 0 0 0 +0.25 0.125 0.375 0 0 0 +0.375 0 0.375 0 0 0 +0.375 0.125 0.25 0 0 0 +0.25 0 0.5 0 0 0 +0.25 0.125 0.625 0 0 0 +0.375 0 0.625 0 0 0 +0.375 0.125 0.5 0 0 0 +0.25 0 0.75 0 0 0 +0.25 0.125 0.875 0 0 0 +0.375 0 0.875 0 0 0 +0.375 0.125 0.75 0 0 0 +0.25 0.25 0 0 0 0 +0.25 0.375 0.125 0 0 0 +0.375 0.25 0.125 0 0 0 +0.375 0.375 0 0 0 0 +0.25 0.25 0.25 0 0 0 +0.25 0.375 0.375 0 0 0 +0.375 0.25 0.375 0 0 0 +0.375 0.375 0.25 0 0 0 +0.25 0.25 0.5 0 0 0 +0.25 0.375 0.625 0 0 0 +0.375 0.25 0.625 0 0 0 +0.375 0.375 0.5 0 0 0 +0.25 0.25 0.75 0 0 0 +0.25 0.375 0.875 0 0 0 +0.375 0.25 0.875 0 0 0 +0.375 0.375 0.75 0 0 0 +0.5 0 0 0 0 0 +0.5 0.125 0.125 0 0 0 +0.625 0 0.125 0 0 0 +0.625 0.125 0 0 0 0 +0.5 0 0.25 0 0 0 +0.5 0.125 0.375 0 0 0 +0.625 0 0.375 0 0 0 +0.625 0.125 0.25 0 0 0 +0.5 0 0.5 0 0 0 +0.5 0.125 0.625 0 0 0 +0.625 0 0.625 0 0 0 +0.625 0.125 0.5 0 0 0 +0.5 0 0.75 0 0 0 +0.5 0.125 0.875 0 0 0 +0.625 0 0.875 0 0 0 +0.625 0.125 0.75 0 0 0 +0.5 0.25 0 0 0 0 +0.5 0.375 0.125 0 0 0 +0.625 0.25 0.125 0 0 0 +0.625 0.375 0 0 0 0 +0.5 0.25 0.25 0 0 0 +0.5 0.375 0.375 0 0 0 +0.625 0.25 0.375 0 0 0 +0.625 0.375 0.25 0 0 0 +0.5 0.25 0.5 0 0 0 +0.5 0.375 0.625 0 0 0 +0.625 0.25 0.625 0 0 0 +0.625 0.375 0.5 0 0 0 +0.5 0.25 0.75 0 0 0 +0.5 0.375 0.875 0 0 0 +0.625 0.25 0.875 0 0 0 +0.625 0.375 0.75 0 0 0 +0.75 0 0 0 0 0 +0.75 0.125 0.125 0 0 0 +0.875 0 0.125 0 0 0 +0.875 0.125 0 0 0 0 +0.75 0 0.25 0 0 0 +0.75 0.125 0.375 0 0 0 +0.875 0 0.375 0 0 0 +0.875 0.125 0.25 0 0 0 +0.75 0 0.5 0 0 0 +0.75 0.125 0.625 0 0 0 +0.875 0 0.625 0 0 0 +0.875 0.125 0.5 0 0 0 +0.75 0 0.75 0 0 0 +0.75 0.125 0.875 0 0 0 +0.875 0 0.875 0 0 0 +0.875 0.125 0.75 0 0 0 +0.75 0.25 0 0 0 0 +0.75 0.375 0.125 0 0 0 +0.875 0.25 0.125 0 0 0 +0.875 0.375 0 0 0 0 +0.75 0.25 0.25 0 0 0 +0.75 0.375 0.375 0 0 0 +0.875 0.25 0.375 0 0 0 +0.875 0.375 0.25 0 0 0 +0.75 0.25 0.5 0 0 0 +0.75 0.375 0.625 0 0 0 +0.875 0.25 0.625 0 0 0 +0.875 0.375 0.5 0 0 0 +0.75 0.25 0.75 0 0 0 +0.75 0.375 0.875 0 0 0 +0.875 0.25 0.875 0 0 0 +0.875 0.375 0.75 0 0 0 +0 0.5 0 0 0 0 +0 0.625 0.125 0 0 0 +0.125 0.5 0.125 0 0 0 +0.125 0.625 0 0 0 0 +0 0.5 0.25 0 0 0 +0 0.625 0.375 0 0 0 +0.125 0.5 0.375 0 0 0 +0.125 0.625 0.25 0 0 0 +0 0.5 0.5 0 0 0 +0 0.625 0.625 0 0 0 +0.125 0.5 0.625 0 0 0 +0.125 0.625 0.5 0 0 0 +0 0.5 0.75 0 0 0 +0 0.625 0.875 0 0 0 +0.125 0.5 0.875 0 0 0 +0.125 0.625 0.75 0 0 0 +0 0.75 0 0 0 0 +0 0.875 0.125 0 0 0 +0.125 0.75 0.125 0 0 0 +0.125 0.875 0 0 0 0 +0 0.75 0.25 0 0 0 +0 0.875 0.375 0 0 0 +0.125 0.75 0.375 0 0 0 +0.125 0.875 0.25 0 0 0 +0 0.75 0.5 0 0 0 +0 0.875 0.625 0 0 0 +0.125 0.75 0.625 0 0 0 +0.125 0.875 0.5 0 0 0 +0 0.75 0.75 0 0 0 +0 0.875 0.875 0 0 0 +0.125 0.75 0.875 0 0 0 +0.125 0.875 0.75 0 0 0 +0.25 0.5 0 0 0 0 +0.25 0.625 0.125 0 0 0 +0.375 0.5 0.125 0 0 0 +0.375 0.625 0 0 0 0 +0.25 0.5 0.25 0 0 0 +0.25 0.625 0.375 0 0 0 +0.375 0.5 0.375 0 0 0 +0.375 0.625 0.25 0 0 0 +0.25 0.5 0.5 0 0 0 +0.25 0.625 0.625 0 0 0 +0.375 0.5 0.625 0 0 0 +0.375 0.625 0.5 0 0 0 +0.25 0.5 0.75 0 0 0 +0.25 0.625 0.875 0 0 0 +0.375 0.5 0.875 0 0 0 +0.375 0.625 0.75 0 0 0 +0.25 0.75 0 0 0 0 +0.25 0.875 0.125 0 0 0 +0.375 0.75 0.125 0 0 0 +0.375 0.875 0 0 0 0 +0.25 0.75 0.25 0 0 0 +0.25 0.875 0.375 0 0 0 +0.375 0.75 0.375 0 0 0 +0.375 0.875 0.25 0 0 0 +0.25 0.75 0.5 0 0 0 +0.25 0.875 0.625 0 0 0 +0.375 0.75 0.625 0 0 0 +0.375 0.875 0.5 0 0 0 +0.25 0.75 0.75 0 0 0 +0.25 0.875 0.875 0 0 0 +0.375 0.75 0.875 0 0 0 +0.375 0.875 0.75 0 0 0 +0.5 0.5 0 0 0 0 +0.5 0.625 0.125 0 0 0 +0.625 0.5 0.125 0 0 0 +0.625 0.625 0 0 0 0 +0.5 0.5 0.25 0 0 0 +0.5 0.625 0.375 0 0 0 +0.625 0.5 0.375 0 0 0 +0.625 0.625 0.25 0 0 0 +0.5 0.5 0.5 0 0 0 +0.5 0.625 0.625 0 0 0 +0.625 0.5 0.625 0 0 0 +0.625 0.625 0.5 0 0 0 +0.5 0.5 0.75 0 0 0 +0.5 0.625 0.875 0 0 0 +0.625 0.5 0.875 0 0 0 +0.625 0.625 0.75 0 0 0 +0.5 0.75 0 0 0 0 +0.5 0.875 0.125 0 0 0 +0.625 0.75 0.125 0 0 0 +0.625 0.875 0 0 0 0 +0.5 0.75 0.25 0 0 0 +0.5 0.875 0.375 0 0 0 +0.625 0.75 0.375 0 0 0 +0.625 0.875 0.25 0 0 0 +0.5 0.75 0.5 0 0 0 +0.5 0.875 0.625 0 0 0 +0.625 0.75 0.625 0 0 0 +0.625 0.875 0.5 0 0 0 +0.5 0.75 0.75 0 0 0 +0.5 0.875 0.875 0 0 0 +0.625 0.75 0.875 0 0 0 +0.625 0.875 0.75 0 0 0 +0.75 0.5 0 0 0 0 +0.75 0.625 0.125 0 0 0 +0.875 0.5 0.125 0 0 0 +0.875 0.625 0 0 0 0 +0.75 0.5 0.25 0 0 0 +0.75 0.625 0.375 0 0 0 +0.875 0.5 0.375 0 0 0 +0.875 0.625 0.25 0 0 0 +0.75 0.5 0.5 0 0 0 +0.75 0.625 0.625 0 0 0 +0.875 0.5 0.625 0 0 0 +0.875 0.625 0.5 0 0 0 +0.75 0.5 0.75 0 0 0 +0.75 0.625 0.875 0 0 0 +0.875 0.5 0.875 0 0 0 +0.875 0.625 0.75 0 0 0 +0.75 0.75 0 0 0 0 +0.75 0.875 0.125 0 0 0 +0.875 0.75 0.125 0 0 0 +0.875 0.875 0 0 0 0 +0.75 0.75 0.25 0 0 0 +0.75 0.875 0.375 0 0 0 +0.875 0.75 0.375 0 0 0 +0.875 0.875 0.25 0 0 0 +0.75 0.75 0.5 0 0 0 +0.75 0.875 0.625 0 0 0 +0.875 0.75 0.625 0 0 0 +0.875 0.875 0.5 0 0 0 +0.75 0.75 0.75 0 0 0 +0.75 0.875 0.875 0 0 0 +0.875 0.75 0.875 0 0 0 +0.875 0.875 0.75 0 0 0 + +H +0 +1 +0.1875 0.1875 0.50 m 0 0 1 v 0 0 1.3 diff --git a/examples/35_tdofdft/01_Al_proton_ESP/02_td/run.sh b/examples/35_tdofdft/01_Al_proton_ESP/02_td/run.sh new file mode 100644 index 00000000000..7623ec13913 --- /dev/null +++ b/examples/35_tdofdft/01_Al_proton_ESP/02_td/run.sh @@ -0,0 +1,18 @@ +#!/bin/bash + +ABACUS_PATH=$(awk -F "=" '$1=="ABACUS_PATH"{print $2}' ../../../../SETENV) +ABACUS_NPROCS=$(awk -F "=" '$1=="ABACUS_NPROCS"{print $2}' ../../../../SETENV) +ABACUS_THREADS=$(awk -F "=" '$1=="ABACUS_THREADS"{print $2}' ../../../../SETENV) + +OMP_NUM_THREADS=${ABACUS_THREADS} mpirun -np ${ABACUS_NPROCS} ${ABACUS_PATH} | tee relax.output + +if [[ ! -f scf.output ]] || + [[ ! -f OUT.ABACUS/running_scf.log ]] || + [[ ! ( "$(tail -1 OUT.ABACUS/running_scf.log)" == " Total Time :"* ) ]] +then + echo "job failed!" + exit 1 +else + echo "job succeeded!" + exit 0 +fi diff --git a/examples/35_tdofdft/README b/examples/35_tdofdft/README new file mode 100644 index 00000000000..a2cb3675cda --- /dev/null +++ b/examples/35_tdofdft/README @@ -0,0 +1,90 @@ +# TDOFDFT Examples + +This directory contains examples for performing Time-Dependent Orbital-Free Density +Functional Theory (TDOFDFT) calculations in ABACUS. + +## What is TDOFDFT? + +TDOFDFT extends OFDFT to the time domain by directly propagating the time-dependent +Pauli-like equations. It is particularly suited for studying electronic stopping power (ESP), +warm dense matter (WDM), and other high-energy-density physics scenarios where the +computational cost of KS-DFT becomes prohibitive. + +### Key Features: +- Linear-scaling with system size (no orbitals, only electron density) +- Natural description of free-electron-like systems (e.g., metals, plasmas) +- Supports Thomas-Fermi (TF), von Weizsäcker (VW), and other kinetic energy functionals +- Current-density-dependent (CD) kinetic potential for improved dynamic response +- Mermin/CD correction (mcd) for finite-temperature effects + +## Example Included + +### 1. 01_Al_proton_ESP +- **System**: Al FCC 4×4×4 supercell (256 atoms) + a moving proton +- **Purpose**: Calculate electronic stopping power by tracking the proton's kinetic energy + loss over time via TDOFDFT +- **Workflow**: + 1. Run `01_gs/` to get the ground-state charge density (`chg.cube`) + 2. Copy `chg.cube` from `01_gs/OUT.autotest/` to `02_td/` + 3. Run `02_td/` for the time-dependent proton dynamics + 4. Use `get_kinetic.py` from `tools/02_postprocessing/tdofdft-tools/` to extract the + proton kinetic energy curve from `OUT.autotest/MD_dump` + +## Input Files + +Each subdirectory contains: +- `INPUT`: Calculation parameters +- `KPT`: k-point sampling +- `STRU`: Atomic structure and pseudopotentials +- `run.sh`: Script to run the calculation + +## Key INPUT Parameters + +- `esolver_type ofdft`: Ground-state OFDFT +- `esolver_type tdofdft`: Time-dependent OFDFT +- `of_kinetic tf+`: Use TF + VW kinetic functional (coefficients both = 1) +- `of_cd 1`: Enable current-density-dependent kinetic potential +- `of_mcd_alpha 2.0`: MCD correction factor (see Li et al., CPB 2026) +- `nelec 768`: For Al 4×4×4 (256 atoms × 3 valence electrons); no electron added for proton +- `init_chg file`: Read initial charge density from `chg.cube` + +## How to Run + +1. First, ground state: + ```bash + cd 01_Al_proton_ESP/01_gs + bash run.sh + ``` + +2. Copy charge density: + ```bash + cp 01_gs/OUT.autotest/chg.cube 02_td/ + ``` + +3. Then, time-dependent: + ```bash + cd 01_Al_proton_ESP/02_td + bash run.sh + ``` + +4. Extract kinetic energy curve: + ```bash + python3 ../../../../tools/02_postprocessing/tdofdft-tools/get_kinetic.py + ``` + +## References + +- White, A. J., et al. "Time-dependent orbital-free density functional theory for + electronic stopping power: Comparison to the Mermin-Kohn-Sham theory at high + temperatures." Physical Review B 98.14 (2018): 144302. +- Li, Yuanbo, et al. "Real-Time Time-Dependent Orbital-Free Density Functional + Theory with Revised Dynamic Kinetic Energy Potential: A Study of Electronic + Stopping Power." Chinese Physics B (2026). + +## Notes + +- The ground-state calculation uses a local pseudopotential (blps) for Al +- The proton (H) uses `1/r` local potential only — it carries no electrons +- Proton initial velocity is set to 1.3 a.u. along z in STRU +- A larger `md_nstep` (~4500+) is needed for sufficient stopping dynamics +- `pseudo_rcut 16` is set for the large Al supercell diff --git a/tools/02_postprocessing/tdofdft-tools/get_kinetic.py b/tools/02_postprocessing/tdofdft-tools/get_kinetic.py new file mode 100644 index 00000000000..23fa5a265da --- /dev/null +++ b/tools/02_postprocessing/tdofdft-tools/get_kinetic.py @@ -0,0 +1,79 @@ +#!/usr/bin/env python3 +""" +Extract proton kinetic energy from TDOFDFT MD_dump output. + +Usage: + python3 get_kinetic.py [md_dump_file] [output_file] + +Default: + md_dump_file = OUT.autotest/MD_dump + output_file = kinetic.txt + +Output format: + time(fs) kinetic_energy(eV) +""" + +import sys +import os + +# ========== Parameters ========== +md_dump_file = "OUT.autotest/MD_dump" +output_file = "kinetic.txt" +md_dt = 0.0005 # fs, must match md_dt in INPUT +m_H = 1.008 # amu, proton mass +# ================================ + +if len(sys.argv) >= 2: + md_dump_file = sys.argv[1] +if len(sys.argv) >= 3: + output_file = sys.argv[2] + +# Conversion factor: E_k(eV) = 0.5 * m(amu) * v(Ang/fs)^2 * factor +amu_to_kg = 1.66053906660e-27 +ang_fs_to_m_s = 1e5 # 1 Ang/fs = 10^5 m/s +eV_to_J = 1.602176634e-19 +factor = amu_to_kg * ang_fs_to_m_s**2 / eV_to_J # ≈ 103.6427 + +if not os.path.exists(md_dump_file): + print(f"Error: {md_dump_file} not found!") + sys.exit(1) + +with open(md_dump_file, "r") as f: + lines = f.readlines() + +kinetic_data = [] +i = 0 +while i < len(lines): + line = lines[i].strip() + if line.startswith("MDSTEP:"): + step = int(line.split(":")[1].strip()) + time_fs = step * md_dt + + i += 1 + # Skip lattice + virial + header lines to reach atom data + while i < len(lines) and "INDEX" not in lines[i] and "MDSTEP:" not in lines[i]: + i += 1 + i += 1 # skip INDEX header + + # Atom data: INDEX LABEL x y z Fx Fy Fz Vx Vy Vz + # Find the last H atom (proton) + last_vx = last_vy = last_vz = 0.0 + while i < len(lines) and not lines[i].strip().startswith("MDSTEP:"): + parts = lines[i].split() + if len(parts) >= 11 and parts[1] == "H": + last_vx = float(parts[8]) + last_vy = float(parts[9]) + last_vz = float(parts[10]) + i += 1 + + v_sq = last_vx**2 + last_vy**2 + last_vz**2 + ek_eV = 0.5 * m_H * v_sq * factor + kinetic_data.append((time_fs, ek_eV)) + continue + i += 1 + +with open(output_file, "w") as f: + for t, ek in kinetic_data: + f.write(f"{t:.4f}\t{ek:.10f}\n") + +print(f"Done! {len(kinetic_data)} steps written to {output_file}") From 322fe98227ccc38032e5370c6d5e3238ee6ad5d0 Mon Sep 17 00:00:00 2001 From: lyb9812 Date: Fri, 19 Jun 2026 11:55:04 +0800 Subject: [PATCH 4/4] Update TDOFDFT README: refine description and terminology Co-Authored-By: Claude --- examples/35_tdofdft/README | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/examples/35_tdofdft/README b/examples/35_tdofdft/README index a2cb3675cda..1a0ea93bb53 100644 --- a/examples/35_tdofdft/README +++ b/examples/35_tdofdft/README @@ -5,17 +5,13 @@ Functional Theory (TDOFDFT) calculations in ABACUS. ## What is TDOFDFT? -TDOFDFT extends OFDFT to the time domain by directly propagating the time-dependent -Pauli-like equations. It is particularly suited for studying electronic stopping power (ESP), -warm dense matter (WDM), and other high-energy-density physics scenarios where the +TDOFDFT extends OFDFT to the time domain by directly propagating the Madelung wavefunction. It is particularly suited for studying electronic stopping power (ESP) of warm dense matter (WDM), and other high-energy-density physics scenarios where the computational cost of KS-DFT becomes prohibitive. ### Key Features: - Linear-scaling with system size (no orbitals, only electron density) -- Natural description of free-electron-like systems (e.g., metals, plasmas) - Supports Thomas-Fermi (TF), von Weizsäcker (VW), and other kinetic energy functionals -- Current-density-dependent (CD) kinetic potential for improved dynamic response -- Mermin/CD correction (mcd) for finite-temperature effects +- Current-dependent (CD) kinetic potential for improved dynamic response ## Example Included