From 0498dc7f2dcc2bb2136fc3f83d0ee48762f9116c Mon Sep 17 00:00:00 2001 From: Shridhar Shanbhag <68288744+ShridharShanbhag@users.noreply.github.com> Date: Fri, 27 Mar 2026 14:46:43 +0100 Subject: [PATCH] Ehrenfest.md and qs_moments (#5019) --- docs/methods/sampling/ehrenfest.md | 29 +++++++++---------- src/qs_moments.F | 10 ++++--- .../regtest-moments-kpoints/TEST_FILES.toml | 2 +- tests/matchers.py | 2 +- 4 files changed, 21 insertions(+), 22 deletions(-) diff --git a/docs/methods/sampling/ehrenfest.md b/docs/methods/sampling/ehrenfest.md index aca7a8fc8f..ca0719d83d 100644 --- a/docs/methods/sampling/ehrenfest.md +++ b/docs/methods/sampling/ehrenfest.md @@ -33,9 +33,9 @@ radiation. The light is then treated classically using, either the length gauge systems, or the velocity gauge for periodic systems. A rather general derivation for Ehrenfest (non-adiabatic quantum-classical molecular) dynamics can -be obtained starting from the action of a system. For a situation in which the the electrons are -treated quantum mechanically while the nuclei are treated classically the total action can be -written as the sum of the two environments +be obtained starting from the action of a system. For a situation in which the electrons are treated +quantum mechanically, while the nuclei are treated classically the total action can be written as +the sum of the two environments $$ A = A_c + A_q \qquad A_c = \int_{t_0}^{t_f} \left[ \sum_A \frac{M_A}{2}\dot{\bf R}_A -U({\bf R},t)\right]dt @@ -51,10 +51,10 @@ $$ Evaluating these expressions in the framework of TDDFT, the equations of motion become $$ -M\ddot{\bf R} = -\frac{\partial}{\partial {\bf R}} U({\bf R},t) - \sum_j\left\langle \Psi^j \left| \frac{\partial}{\partial{\bf R}} V_{\text{int}}({\bf r},{\bf R})\right| \Psi^j \right\rangle +M\ddot{\bf R} = -\frac{\partial}{\partial {\bf R}} U({\bf R},t) - \sum_j\left\langle \psi^j \left| \frac{\partial}{\partial{\bf R}} V_{\text{int}}({\bf r},{\bf R})\right| \psi^j \right\rangle $$ -for the nuclear motion, wile for the electrons the time dependent SE as given above is used. These +for the nuclear motion, while for the electrons the time dependent SE as given above is used. These equations are valid for the exact wavefunctions. In the Kohn-Sham approach, the wavefunctions are replaced by a linear combination of basis functions. For plane waves the equations remain the same, since they do not depend on the nuclear coordinates and therefore are independent of the change of @@ -104,8 +104,8 @@ usual [WFN_RESTART_FILE_NAME](#CP2K_INPUT.FORCE_EVAL.DFT.WFN_RESTART_FILE_NAME) Three different propagators are available in CP2K, the enforced time-reversible symmetry propagator (ETRS), the exponential midpoint (EM) propagator, and the Crank-Nicholson propagator which can be -seen as a first order Pad'e approximation of the EM propagator. The ETRS approach starts with an -exponential approximation to the evolution operator $\hat{U}(t,0)=\exp{-it\hat{H}(0)}$, to then +seen as a first order Padé approximation of the EM propagator. The ETRS approach starts with an +exponential approximation to the evolution operator $\hat{U}(t,0)=\exp{[-it\hat{H}(0)]}$, to then compute the final time-reversible and unitary propagator self-consistently. In the real time propagation scheme (fixed ionic positions) the self consistent solution only involves the calculation of the new Kohn-Sham matrix for the propagated coefficients. For Ehrenfest MD, the @@ -123,10 +123,10 @@ given by `EPS_ITER`. In terms of computational cost, the most expensive part is the evaluation of the matrix exponential in the propagator. Four different methods among those listed in have been implemented in CP2K, i.e., -the Taylor expansion, he diagonalization, the Pad'e approximation, and the Arnoldi subspace +the Taylor expansion, he diagonalization, the Padé approximation, and the Arnoldi subspace iteration, to be selected via the `MAT_ESP` input key. The Arnoldi method often provides a superior performance. Comparing the theoretical scaling, the Arnoldi method is expected to be about 5 times -as fast as Pad'e or Taylor. However, the Pade approximation can be sometime the faster and more +as fast as Padé or Taylor. However, the Padé approximation can be sometime the faster and more stable choice than the Arnoldi method (e.g., large time step). Since propagation schemes require an iterative procedure, it is convenient to apply an extrapolation scheme in order to speed up convergence. For Ehrenfest dynamics, an extrapolation based on the product of the density and the @@ -140,7 +140,7 @@ keyword). ## The external electric field -One of the most relevant application domains for RT-TDDFT is the study of light?matter interactions, +One of the most relevant application domains for RT-TDDFT is the study of light-matter interactions, e.g., in the field of spectroscopy, excited state dynamics and radiation damage. To mimic these phenomena, at any time during the propagation, it is possible to apply a time dependent electric field ${\bf E}(t)$. The applied field is in general modulated by an envelope function, @@ -166,15 +166,12 @@ $$ In the time dependent KS equations, the vector potential ${\bf A}(t)$ appears in the kinetic energy term and, in the case where non-local pseudopotentials are used, the gauge field also transforms the -electron?ion interaction. To use this representation the `VELOCITY_GAUGE` input keyword has to be +electron-ion interaction. To use this representation the `VELOCITY_GAUGE` input keyword has to be activated. The total energy varies with time as the applied external field interacts with the system. To monitor the time evolution of the field as well as of the various terms contributing to the total electronic energy, the corresponding print key can be activated from the [REAL_TIME_PROPAGATION](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION) section. -with $\Delta C$ as the difference of coefficient matrices in two successive steps, and $\epsilon$ -given by `EPS_ITER`. - ## Resonant Excitation The input file provided below starts a RT--TDDFT simulation for an isolated carbon monoxide molecule @@ -203,7 +200,7 @@ while monitoring density differences is usually helpful to detect charge transfe By projecting the time-dependent molecular orbitals onto some reference states (e.g., the initial MOs) means computing the overlap between the propagated orbital -$\psi_i({\bf r}, t) = \sum_\alpha C_{i\alpha}(t)\phi({\bf r})$ and any reference orbital +$\psi_i({\bf r}, t) = \sum_\alpha C_{i\alpha}(t)\phi_\alpha({\bf r})$ and any reference orbital $\psi^{\text{ref}}_m$. For instance, considering as reference orbitals the static unoccupied ones, the quantity @@ -497,7 +494,7 @@ XAS_TDP will save one .wfn file per printed excited state, and the projection du by: $$ -n_{\omega}^i(t) = |<\omega | \psi_i(t)>|^2 = |\sum_{ab} \left( C_\omega^a \right)^* c_i^b(t) S_{ab} |^2 +n_{\omega}^i(t) = |\langle\omega | \psi_i(t)\rangle|^2 = |\sum_{ab} \left( C_\omega^a \right)^* c_i^b(t) S_{ab} |^2 $$ where $C_{\omega\alpha}$ is the $a^{\text{th}}$ atomic coefficient of the excited state found in the diff --git a/src/qs_moments.F b/src/qs_moments.F index f19cf48ba3..248bb69267 100644 --- a/src/qs_moments.F +++ b/src/qs_moments.F @@ -3641,6 +3641,7 @@ CONTAINS n_img_scf, n_img_all, nao, & num_pe, num_copy, mepos, n, m, mu, & ispin, nspin + INTEGER, DIMENSION(3) :: periodic INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_all INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_all REAL(KIND=dp), DIMENSION(3), OPTIONAL :: rcc @@ -3680,7 +3681,7 @@ CONTAINS mos=mos) CALL get_mo_set(mo_set=mos(1), nao=nao) - CALL get_cell(cell=cell, h=hmat) + CALL get_cell(cell=cell, h=hmat, periodic=periodic) nspin = SIZE(matrix_ks_kp, 1) nkp = SIZE(xkp, 2) @@ -3720,8 +3721,8 @@ CONTAINS IF (calc_bc) ALLOCATE (berry_c(nspin, num_copy, 3, nao), source=0.0_dp) !$OMP PARALLEL DEFAULT(NONE) PRIVATE(ikp, S_k, H_k, eigenvals, C_k, ispin, n, m, & -!$OMP dS_dk_i, dH_dk_i, D_k, dip, bc, C_dS_C, C_dH_C, CDC, tmp_max, phase) & -!$OMP SHARED(num_pe, mepos, dipole, berry_c, nao, nspin, & +!$OMP i_dir, dS_dk_i, dH_dk_i, D_k, dip, bc, C_dS_C, C_dH_C, CDC, tmp_max, phase) & +!$OMP SHARED(num_pe, mepos, dipole, berry_c, nao, nspin, periodic, & !$OMP nkp, xkp, S_rs, H_rs, D_rs, index_to_cell_all, hmat, calc_bc) ALLOCATE (dS_dk_i(nao, nao), C_dS_C(nao, nao), dH_dk_i(nao, nao), C_dH_C(nao, nao), source=z_zero) ALLOCATE (CDC(nao, nao), dip(3, nao, nao), S_k(nao, nao), H_k(nao, nao), source=z_zero) @@ -3757,6 +3758,7 @@ CONTAINS DO i_dir = 1, 3 ! d^x, d^y, d^z + IF (periodic(i_dir) == 0) CYCLE ! ∇ S(k) = Σ_R iR S^R e^(ikR), ∇ H(k) = Σ_R iR H^R e^(ikR) CALL rs_to_kp(S_rs(1, :, :, :), dS_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat) CALL rs_to_kp(H_rs(ispin, :, :, :), dH_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat) @@ -3783,7 +3785,7 @@ CONTAINS END DO END DO ! Compute the Berry curvature from the dipoles - ! Ω^γ_n = Σ_m 2*Im[d^α_nm (d^β_mn)*], where, α, β, γ belong to {x, y, z} + ! Ω^γ_n = Σ_m 2*Im[d^α_nm d^β_mn], where, α, β, γ belong to {x, y, z} IF (calc_bc) THEN bc = 0.0_dp DO i_dir = 1, 3 diff --git a/tests/QS/regtest-moments-kpoints/TEST_FILES.toml b/tests/QS/regtest-moments-kpoints/TEST_FILES.toml index df309a4b8b..38e19868bf 100644 --- a/tests/QS/regtest-moments-kpoints/TEST_FILES.toml +++ b/tests/QS/regtest-moments-kpoints/TEST_FILES.toml @@ -13,5 +13,5 @@ # and ensure matching of Berry Curvature "C2_pbe_moment_kpset.inp" = [{matcher="BC_near_K_point", tol=1.0E-01, ref=-7100}] # compute dipole for open shell system of monolayer CrSBr -"CrSBr_open_shell.inp" = [{matcher="Dipole_for_CrSBr", tol=1.0E-06, ref=-1.25}] +"CrSBr_open_shell.inp" = [{matcher="Dipole_for_CrSBr", tol=1.0E-06, ref=0.298E-01}] #EOF diff --git a/tests/matchers.py b/tests/matchers.py index 2fce794d6f..6647828e92 100644 --- a/tests/matchers.py +++ b/tests/matchers.py @@ -274,7 +274,7 @@ registry["M127"] = GenericMatcher(r"Checksum (Acoustic Sum Rule):", col=5) registry["Dipole_at_kp_1"] = GenericMatcher(r" 1 1 2", col=4) # Dipole moment calculated at a specific k-point (-0.375,-0.375, 0.00) -registry["Dipole_for_CrSBr"] = GenericMatcher(r" 1 31 32", col=8) +registry["Dipole_for_CrSBr"] = GenericMatcher(r" 1 31 32", col=4) # Berry curvature calculated from dipoles near K point in graphene BZ registry["BC_near_K_point"] = GenericMatcher(r" 1 4", col=5)