Ehrenfest.md and qs_moments (#5019)

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

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

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

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)