HP bugfix

Doc/release-notes

@@ -1,6 +1,10 @@
 Fixed in development version:
   * Symmetrization glitches on many processors (reported by Yunqi Shao)
   * Crash for noncolinear PAW case on many processors (reported by I. Timrov)
+  * The HP code for the magnetic noncollinear case had a bug
+    when symmetry was used (some symmetry arrays were not re-initialized
+    when perturbing Hubbard atoms) and the magnetization direction flip was 
+    missing in the Hubbard term (for the second Sternheimer equation) (fixed by L. Binci)
 
 New in development version:
   * Interface between Wannier90 and DFT+U to use Wannier functions as

HP/Doc/INPUT_HP.def

@@ -322,5 +322,15 @@ input_description -distribution {Quantum ESPRESSO} -package PWscf -program hp.x
             }
         }
 
+        var no_metq0 -type LOGICAL {
+            default { .false. }
+            info {
+                If .true. the metallic response term at q=0 is ignored
+                (i.e. the last term in Eq. (22) in PRB 103, 045141 (2021)).
+                This is useful for magnetic insulators to avoid the divergence
+                of the calculation.
+            }
+        }
+
 }  
 }

HP/examples/example02/README

@@ -1,5 +1,6 @@
-In order to compute U for magnetic insulators one needs to follow
+There are two ways how to compute U for magnetic insulators.
-special procedure. Since we do not know a priori the magnetization
+
+OPTION 1: Since we do not know a priori the magnetization
 of the system, we perform a calculation with a finite 
 starting_magnetization (in the range from -1 to +1) which requires 
 to treat the system of interest as a metal (i.e. to use a smearing). 
@@ -52,3 +53,8 @@ DFPT calculation. The solution to this problem is the following:
      constrained magnetization). 
 
 3. Perform the linear-response calculation of U using the HP code.
+
+OPTION 2: Perform step 1 above. Then jump directly to step 3 by setting
+          no_metq0 = .true. in the input file for the HP calculation. 
+          This will skip the metallic response term at q=0 (i.e. the last 
+          term in Eq. (22) in PRB 103, 045141 (2021)).

HP/examples/example02/run_example2

@@ -0,0 +1,152 @@
+#!/bin/sh
+
+# run from directory where this script is
+cd `dirname $0`
+EXAMPLE_DIR=`pwd`
+
+# check whether ECHO has the -e option
+if test "`echo -e`" = "-e" ; then ECHO=echo ; else ECHO="echo -e" ; fi
+
+$ECHO
+$ECHO "$EXAMPLE_DIR : starting"
+$ECHO
+$ECHO "This example shows how to use pw.x and hp.x to calculate"
+$ECHO "the Hubbard U parameter for Ni in NiO."
+
+# set the needed environment variables
+. ../../../environment_variables
+
+# required executables and pseudopotentials
+BIN_LIST="pw.x hp.x"
+PSEUDO_LIST="Ni.pbesol-n-rrkjus_psl.0.1.UPF O.pbesol-n-rrkjus_psl.0.1.UPF"
+
+$ECHO
+$ECHO "  executables directory: $BIN_DIR"
+$ECHO "  pseudo directory:      $PSEUDO_DIR"
+$ECHO "  temporary directory:   $TMP_DIR"
+$ECHO
+$ECHO "  checking that needed directories and files exist...\c"
+
+# check for directories
+for DIR in "$BIN_DIR" "$PSEUDO_DIR" ; do
+    if test ! -d $DIR ; then
+        $ECHO
+        $ECHO "ERROR: $DIR not existent or not a directory"
+        $ECHO "Aborting"
+        exit 1
+    fi
+done
+for DIR in "$TMP_DIR" "$EXAMPLE_DIR/results" ; do
+    if test ! -d $DIR ; then
+        mkdir $DIR
+    fi
+done
+cd $EXAMPLE_DIR/results
+
+# check for executables
+for FILE in $BIN_LIST ; do
+    if test ! -x $BIN_DIR/$FILE ; then
+        $ECHO
+        $ECHO "ERROR: $BIN_DIR/$FILE not existent or not executable"
+        $ECHO "Aborting"
+        exit 1
+    fi
+done
+$ECHO " done"
+
+# check for pseudopotentials
+for FILE in $PSEUDO_LIST ; do
+    if test ! -r $PSEUDO_DIR/$FILE ; then
+       $ECHO
+       $ECHO "Downloading $FILE to $PSEUDO_DIR...\c"
+            $WGET $PSEUDO_DIR/$FILE $NETWORK_PSEUDO/$FILE 2> /dev/null
+    fi
+    if test $? != 0; then
+        $ECHO
+        $ECHO "ERROR: $PSEUDO_DIR/$FILE not existent or not readable"
+        $ECHO "Aborting"
+        exit 1
+    fi
+done
+$ECHO " done"
+
+# how to run executables
+PW_COMMAND="$PARA_PREFIX $BIN_DIR/pw.x $PARA_POSTFIX"
+HP_COMMAND="$PARA_PREFIX $BIN_DIR/hp.x $PARA_POSTFIX"
+$ECHO
+$ECHO "  running pw.x as:   $PW_COMMAND"
+$ECHO "  running hp.x as:   $HP_COMMAND"
+$ECHO
+
+# clean TMP_DIR
+$ECHO "  cleaning $TMP_DIR...\c"
+rm -rf $TMP_DIR/*
+$ECHO " done"
+
+PREFIX='NiO'
+
+# First self-consistent calculation
+cat > $PREFIX.scf.in << EOF
+ &control
+    calculation='scf'
+    restart_mode='from_scratch',
+    prefix='$PREFIX'
+    pseudo_dir = '$PSEUDO_DIR/'
+    outdir='$TMP_DIR/'
+    verbosity='high'
+ /
+ &system
+    ibrav = 0,
+    celldm(1) = 7.88,
+    nat = 4,
+    ntyp = 3,
+    ecutwfc = 50.0,
+    ecutrho = 400.0,
+    occupations = 'smearing',
+    smearing = 'gauss',
+    degauss = 0.001,
+    nspin = 2,
+    starting_magnetization(1) =  0.5,
+    starting_magnetization(2) = -0.5,
+ /
+ &electrons
+    conv_thr =  1.d-10
+    mixing_beta = 0.7
+ /
+ATOMIC_SPECIES
+ Ni1  58.693  Ni.pbesol-n-rrkjus_psl.0.1.UPF
+ Ni2  58.693  Ni.pbesol-n-rrkjus_psl.0.1.UPF
+ O    16.000  O.pbesol-n-rrkjus_psl.0.1.UPF
+ATOMIC_POSITIONS {alat}
+ Ni1  0.0000000   0.0000000   0.0000000
+ Ni2  0.5000000   0.5000000   0.0000000
+ O    0.5000000   0.0000000   0.0000000
+ O    1.0000000   0.5000000   0.0000000
+CELL_PARAMETERS {alat}
+ 1.00000000  0.50000000  0.50000000
+ 0.50000000  1.00000000  0.50000000
+ 0.50000000  0.50000000  1.00000000
+K_POINTS {automatic}
+ 2 2 2 0 0 0
+HUBBARD {ortho-atomic}
+U Ni1-3d 0.0001
+U Ni2-3d 0.0001
+EOF
+$ECHO "  Running the SCF calculation for $PREFIX..."
+$PW_COMMAND < $PREFIX.scf.in > $PREFIX.scf.out
+$ECHO " done"
+
+# Perform the linear-response calculation
+cat > $PREFIX.hp.in << EOF
+ &inputhp
+   prefix = '$PREFIX',
+   outdir = '$TMP_DIR/',
+   nq1 = 2, nq2 = 2, nq3 = 2,
+   conv_thr_chi = 1.0d-8,
+   iverbosity = 2
+   no_metq0 = .true.
+ /
+EOF
+$ECHO "  Running the linear-response calculation of Hubbard U..."
+$HP_COMMAND < $PREFIX.hp.in > $PREFIX.hp.out
+$ECHO " done"

HP/src/hp_allocate_q.f90

@@ -25,8 +25,8 @@ subroutine hp_allocate_q
   USE lrus,                 ONLY : becp1
   USE eqv,                  ONLY : dpsi, evq, dmuxc, dvpsi
   USE control_lr,           ONLY : lgamma
-  USE ldaU,                 ONLY : Hubbard_lmax, nwfcU
+  USE ldaU,                 ONLY : Hubbard_lmax, nwfcU, lda_plus_u_kind, max_num_neighbors
-  USE ldaU_lr,              ONLY : swfcatomk, swfcatomkpq
+  USE ldaU_lr,              ONLY : swfcatomk, swfcatomkpq, vh_u_save, vh_uv_save
   USE qpoint_aux,           ONLY : becpt
   USE lr_nc_mag,            ONLY : deeq_nc_save
   USE uspp_param,           ONLY : nhm 
@@ -70,6 +70,12 @@ subroutine hp_allocate_q
      ALLOCATE (swfcatomkpq(npwx*npol,nwfcU))
   ENDIF
   !
+  IF (lda_plus_u_kind == 0) THEN
+     ALLOCATE (vh_u_save(2*Hubbard_lmax+1, 2*Hubbard_lmax+1, nspin, nat, 2))
+  ELSEIF (lda_plus_u_kind == 2) THEN
+     ALLOCATE (vh_uv_save(2*Hubbard_lmax+1, 2*Hubbard_lmax+1, max_num_neighbors, nat, nspin, 2))
+  ENDIF
+  !
   RETURN
   !
 end subroutine hp_allocate_q

HP/src/hp_bcast_input.f90

@@ -1,5 +1,5 @@
 !
-! Copyright (C) 2001-2018 Quantum ESPRESSO group
+! Copyright (C) 2001-2025 Quantum ESPRESSO group
 ! This file is distributed under the terms of the
 ! GNU General Public License. See the file `License'
 ! in the root directory of the present distribution,
@@ -27,7 +27,8 @@ SUBROUTINE hp_bcast_input ( )
                                background, compute_hp, sum_pertq, perturb_only_atom,   &
                                determine_num_pert_only, skip_equivalence_q, niter_max, &
                                disable_type_analysis, docc_thr, num_neigh, lmin, rmax, &
-                               nmix, nq1, nq2, nq3, dist_thr, determine_q_mesh_only
+                               nmix, nq1, nq2, nq3, dist_thr, determine_q_mesh_only,   &
+                               no_metq0
   !
   IMPLICIT NONE
   !
@@ -45,6 +46,7 @@ SUBROUTINE hp_bcast_input ( )
   CALL mp_bcast (determine_num_pert_only, meta_ionode_id, world_comm)
   CALL mp_bcast (determine_q_mesh_only, meta_ionode_id, world_comm)
   CALL mp_bcast (disable_type_analysis, meta_ionode_id, world_comm)
+  CALL mp_bcast (no_metq0, meta_ionode_id, world_comm)
   !
   ! Integers
   !

HP/src/hp_check_type.f90

@@ -1,5 +1,5 @@
 !
-! Copyright (C) 2001-2021 Quantum ESPRESSO group
+! Copyright (C) 2001-2025 Quantum ESPRESSO group
 ! This file is distributed under the terms of the
 ! GNU General Public License. See the file `License'
 ! in the root directory of the present distribution,
@@ -23,12 +23,16 @@ subroutine hp_check_type(na)
   !   other atoms (this is needed in order to mimic an isolated
   !   perturbation in a supercell approach).
   ! 
+  USE kinds,              ONLY : DP
   USE ions_base,          ONLY : ityp, nat, ntyp => nsp, tau
   USE io_global,          ONLY : stdout
   USE symm_base,          ONLY : nsym, set_sym, ft, nofrac
-  USE noncollin_module,   ONLY : nspin_mag, m_loc
+  USE noncollin_module,   ONLY : nspin_mag, m_loc, ux, noncolin, &
+                                 domag, angle1, angle2
   USE fft_base,           ONLY : dfftp
   USE ldaU_hp,            ONLY : recalc_sym
+  USE xc_lib,             ONLY : xclib_dft_is
+  USE lsda_mod,           ONLY : starting_magnetization
   !
   IMPLICIT NONE
   !
@@ -87,6 +91,19 @@ subroutine hp_check_type(na)
      IF (.not.ALLOCATED(m_loc)) ALLOCATE(m_loc(3,nat))
      m_loc(:,:) = 0.0d0 
      !
+     IF (noncolin.and.domag) THEN
+        DO na_ = 1, nat
+           m_loc(1,na_) = starting_magnetization(ityp(na_)) * &
+                       SIN( angle1(ityp(na_)) ) * COS( angle2(ityp(na_)) )
+           m_loc(2,na_) = starting_magnetization(ityp(na_)) * &
+                       SIN( angle1(ityp(na_)) ) * SIN( angle2(ityp(na_)) )
+           m_loc(3,na_) = starting_magnetization(ityp(na_)) * &
+                       COS( angle1(ityp(na_)) )
+        END DO
+        ux=0.0_DP
+        if (xclib_dft_is('gradient')) call compute_ux(m_loc,ux,nat)
+     ENDIF
+     !
      CALL set_sym (nat, tau, ityp, nspin_mag, m_loc)
      !
      IF (nsym > nsym_old) THEN

HP/src/hp_dealloc_q.f90

@@ -22,7 +22,7 @@ SUBROUTINE hp_dealloc_q()
                                   & dvxc_s, vsgga, segni
   USE eqv,                 ONLY : dmuxc, dpsi, dvpsi, evq
   USE control_lr,          ONLY : lgamma, nbnd_occ
-  USE ldaU_lr,             ONLY : swfcatomk, swfcatomkpq
+  USE ldaU_lr,             ONLY : swfcatomk, swfcatomkpq, vh_u_save, vh_uv_save
   USE qpoint_aux,          ONLY : ikmks, ikmkmqs, becpt
   USE lr_nc_mag,           ONLY : deeq_nc_save
   !
@@ -62,6 +62,9 @@ SUBROUTINE hp_dealloc_q()
   endif
   if (allocated(deeq_nc_save)) deallocate(deeq_nc_save)
   !
+  if (allocated(vh_u_save))     deallocate (vh_u_save)
+  if (allocated(vh_uv_save))    deallocate (vh_uv_save)
+  !
   ! GGA-specific arrays
   !
   if (allocated(dvxc_rr))         deallocate (dvxc_rr)

HP/src/hp_readin.f90

@@ -1,5 +1,5 @@
 !
-! Copyright (C) 2001-2023 Quantum ESPRESSO group
+! Copyright (C) 2001-2025 Quantum ESPRESSO group
 ! This file is distributed under the terms of the
 ! GNU General Public License. See the file `License'
 ! in the root directory of the present distribution,
@@ -27,7 +27,7 @@ SUBROUTINE hp_readin()
                                skip_equivalence_q, tmp_dir_save, niter_max, dist_thr,  &
                                disable_type_analysis, docc_thr, num_neigh, lmin, rmax, &
                                nmix, nq1, nq2, nq3, alpha_mix, start_q, last_q, maxter, &
-                               determine_q_mesh_only
+                               determine_q_mesh_only, no_metq0
   USE paw_variables,    ONLY : okpaw
   !
   IMPLICIT NONE
@@ -44,7 +44,7 @@ SUBROUTINE hp_readin()
                          niter_max, alpha_mix, nmix, compute_hp, perturb_only_atom,   &
                          start_q, last_q, sum_pertq, ethr_nscf, num_neigh, lmin,      &
                          determine_num_pert_only, disable_type_analysis, docc_thr,    &
-                         determine_q_mesh_only
+                         determine_q_mesh_only, no_metq0
   !
   ! Note: meta_ionode is a single processor that reads the input
   !       Data read from input is subsequently broadcast to all processors
@@ -87,6 +87,7 @@ SUBROUTINE hp_readin()
   nmix               = 4     
   max_seconds        = 1.E+7_DP
   lrpa               = .FALSE.   ! Needed in dv_of_drho
+  no_metq0           = .FALSE.
   CALL get_environment_variable( 'ESPRESSO_TMPDIR', outdir )
   IF ( TRIM( outdir ) == ' ' ) outdir = './'
   !

HP/src/hp_setup_q.f90

@@ -1,5 +1,5 @@
 !
-! Copyright (C) 2001-2023 Quantum ESPRESSO group
+! Copyright (C) 2001-2025 Quantum ESPRESSO group
 ! This file is distributed under the terms of the
 ! GNU General Public License. See the file `License'
 ! in the root directory of the present distribution,
@@ -66,6 +66,7 @@ SUBROUTINE hp_setup_q()
   USE control_flags,    ONLY : noinv
   USE eqv,              ONLY : dmuxc
   USE qpoint,           ONLY : xq
+  USE xc_lib,           ONLY : xclib_dft_is
   USE control_lr,       ONLY : lgamma
   USE lr_symm_base,     ONLY : gi, gimq, irotmq, minus_q, invsymq, nsymq, rtau
   USE ldaU_hp,          ONLY : niter_max, alpha_mix, skip_equivalence_q
@@ -73,10 +74,13 @@ SUBROUTINE hp_setup_q()
   USE control_flags,    ONLY : modenum
   USE lr_nc_mag,        ONLY : deeq_nc_save
   USE dfunct,           ONLY : newd
+  USE ldaU_lr,          ONLY : vh_u_save, vh_uv_save
+  USE ldaU,             ONLY : lda_plus_u_kind, Hubbard_lmax, max_num_neighbors, nsg, v_nsg
   !
   IMPLICIT NONE
   INTEGER :: ir, isym, ik, it, na
   LOGICAL :: sym(48), magnetic_sym
+  COMPLEX(DP), ALLOCATABLE :: ns_nc(:,:,:,:,:), nsg_nc(:,:,:,:,:,:)
   !
   CALL start_clock ('hp_setup_q')
   !
@@ -104,6 +108,7 @@ SUBROUTINE hp_setup_q()
                      COS( angle1(ityp(na)) )
      END DO
      ux=0.0_DP
+     if (xclib_dft_is('gradient')) call compute_ux(m_loc,ux,nat)
      !
      ! Change the sign of the magnetic field in the screened US coefficients
      ! and save also the coefficients computed with -B_xc.
@@ -117,6 +122,33 @@ SUBROUTINE hp_setup_q()
        deeq_nc(:,:,:,:)=deeq_nc_save(:,:,:,:,1)
      ENDIF
   ENDIF
+  !
+  ! Calculate the unperturbed Hubbard potential
+  ! with the reversed Hubbard magnetization
+  ! and save it. Used for the unperturbed
+  ! Hamiltonian in the Sternheimer linear system.
+  !
+  IF (noncolin .and. domag) THEN
+    IF (lda_plus_u_kind == 0) THEN
+        vh_u_save = (0.d0, 0.d0)
+        ALLOCATE (ns_nc(2*Hubbard_lmax+1, 2*Hubbard_lmax+1, nspin, nat, 1))
+        ns_nc = (0.d0, 0.d0)
+        vh_u_save(:,:,:,:,1) = v%ns_nc(:,:,:,:)
+        ns_nc(:,:,:,:,1) = rho%ns_nc(:,:,:,:)
+        CALL revert_mag_u ( ns_nc(:,:,:,:,1) )
+        CALL calc_vh_u (ns_nc(:,:,:,:,1), vh_u_save(:,:,:,:,2))
+        DEALLOCATE(ns_nc)
+     ELSEIF(lda_plus_u_kind == 2) THEN
+        vh_uv_save = (0.d0, 0.d0)
+        ALLOCATE (nsg_nc(2*Hubbard_lmax+1, 2*Hubbard_lmax+1, max_num_neighbors, nat, nspin, 1))
+        nsg_nc = (0.d0, 0.d0)
+        vh_uv_save(:,:,:,:,:,1) = v_nsg(:,:,:,:,:)
+        nsg_nc(:,:,:,:,:,1) = nsg(:,:,:,:,:)
+        CALL revert_mag_uv ( nsg_nc(:,:,:,:,:,1) )
+        CALL calc_vh_uv (nsg_nc(:,:,:,:,:,1), vh_uv_save(:,:,:,:,:,2))
+        DEALLOCATE(nsg_nc)
+     ENDIF
+  ENDIF
   ! 
   ! 4) Compute the derivative of the XC potential (dmuxc)
   !

HP/src/hp_solve_linear_system.f90

@@ -1,5 +1,5 @@
 !
-! Copyright (C) 2001-2023 Quantum ESPRESSO group
+! Copyright (C) 2001-2025 Quantum ESPRESSO group
 ! This file is distributed under the terms of the
 ! GNU General Public License. See the file `License'
 ! in the root directory of the present distribution,
@@ -22,7 +22,7 @@ SUBROUTINE hp_solve_linear_system (na, iq)
   USE wavefunctions,        ONLY : evc
   USE klist,                ONLY : lgauss, ltetra, nelec, ngk
   USE gvecs,                ONLY : doublegrid
-  USE scf,                  ONLY : rho
+  USE scf,                  ONLY : rho, v, vrs
   USE fft_base,             ONLY : dfftp, dffts
   USE lsda_mod,             ONLY : lsda, current_spin, isk
   USE wvfct,                ONLY : nbnd, npwx
@@ -45,10 +45,10 @@ SUBROUTINE hp_solve_linear_system (na, iq)
   USE fft_helper_subroutines
   USE fft_interfaces,       ONLY : fft_interpolate
   USE lr_symm_base,         ONLY : irotmq, minus_q, nsymq, rtau
-  USE ldaU_lr,              ONLY : dnsscf
+  USE ldaU_lr,              ONLY : dnsscf, vh_u_save, vh_uv_save
   USE ldaU_hp,              ONLY : thresh_init, dns0, trace_dns_tot_old, &
                                    conv_thr_chi_best, iter_best, niter_max, nmix, &
-                                   alpha_mix, code, lrdvwfc, iudvwfc
+                                   alpha_mix, code, lrdvwfc, iudvwfc, no_metq0
   USE apply_dpot_mod,       ONLY : apply_dpot_allocate, apply_dpot_deallocate
   USE efermi_shift,         ONLY : ef_shift, def
   USE response_kernels,     ONLY : sternheimer_kernel
@@ -56,6 +56,7 @@ SUBROUTINE hp_solve_linear_system (na, iq)
   USE lsda_mod,             ONLY : nspin
   USE lr_nc_mag,            ONLY : lr_apply_time_reversal, deeq_nc_save, int3_nc_save
   USE lr_symm_base,         ONLY : lr_npert, upert, upert_mq
+  USE ldaU,                 ONLY : lda_plus_u_kind, nsg, v_nsg
   !
   IMPLICIT NONE
   !
@@ -199,6 +200,15 @@ SUBROUTINE hp_solve_linear_system (na, iq)
   !
   lmetq0 = (lgauss .OR. ltetra) .AND. lgamma
   !
+  ! If the user specified in the input that no_metq0=.true. this means
+  ! we want to force remove the metallic term at q=0 (e.g. for magnetic insulators
+  ! when the smearing is used). We remove the last term in Eq. (22) 
+  ! in PRB 103, 045141 (2021) (and also for the response charge density), 
+  ! otherwise this term can diverge for magnetic insulators
+  ! due to a division by DOS(E_Fermi) which can be vanishing.
+  !
+  IF (no_metq0) lmetq0 = .FALSE.
+  !
   IF (lmetq0) THEN
      ALLOCATE (ldos (dfftp%nnr, nspin_mag))
      ALLOCATE (ldoss(dffts%nnr, nspin_mag))
@@ -253,6 +263,16 @@ SUBROUTINE hp_solve_linear_system (na, iq)
      IF (noncolin) dbecsum_nc = (0.d0, 0.d0)
      !
      DO isolv = 1, nsolv
+        !
+        !  change the sign of the magnetic field if required
+        !
+        IF (isolv == 2) THEN
+           IF (lda_plus_u_kind == 0) THEN
+              v%ns_nc (:,:,:,:) = vh_u_save(:,:,:,:,2)
+           ELSEIF (lda_plus_u_kind == 2) THEN
+              v_nsg(:,:,:,:,:) = vh_uv_save(:,:,:,:,:,2)
+           ENDIF
+        ENDIF
         !
         ! set threshold for iterative solution of the linear system
         !
@@ -278,6 +298,16 @@ SUBROUTINE hp_solve_linear_system (na, iq)
            WRITE(stdout, '(6x, "sternheimer_kernel not converged. Try to increase thresh_init.")')
         ENDIF
         !
+        !  reset the original magnetic field if it was changed
+        !
+        IF (isolv == 2) THEN
+           IF (lda_plus_u_kind == 0) THEN
+               v%ns_nc (:,:,:,:) = vh_u_save(:,:,:,:,1)
+           ELSEIF(lda_plus_u_kind == 2) THEN
+               v_nsg (:,:,:,:,:) = vh_uv_save(:,:,:,:,:,1)
+           ENDIF
+        ENDIF
+        !
      ENDDO ! isolv
      !
      IF (nsolv==2) THEN

HP/src/hp_sym_dmag.f90

@@ -1,5 +1,5 @@
 !
-! Copyright (C) 2001-2023 Quantum ESPRESSO group
+! Copyright (C) 2001-2025 Quantum ESPRESSO group
 ! This file is distributed under the terms of the
 ! GNU General Public License. See the file `License'
 ! in the root directory of the present distribution,
@@ -251,6 +251,8 @@ subroutine hp_sym_dmag (dmagtosym)
                 !
                 if (t_rev(isym) == 1) then 
                    mag(:) = conjg(mag(:)) * phase2 (isym)
+                else
+                   mag(:) = mag(:) * phase2 (isym)
                 end if
                 !
                 dmagsym(i,j,k,1)=dmagsym(i,j,k,1)+mag(1)

HP/src/hpcom.f90

@@ -1,5 +1,5 @@
 !
-! Copyright (C) 2001-2022 Quantum ESPRESSO group
+! Copyright (C) 2001-2025 Quantum ESPRESSO group
 ! This file is distributed under the terms of the
 ! GNU General Public License. See the file `License'
 ! in the root directory of the present distribution,
@@ -29,6 +29,7 @@ MODULE ldaU_hp
              determine_q_mesh_only,   &     ! If .true. determine the q mesh for a given perturbed atom and exit
              determine_num_pert_only, &     ! If .true. determine only which atoms must be perterbed
              skip_equivalence_q,      &     ! If .true. the full frid of q points will be used
+             no_metq0,                &     ! If .true. the metallic response term at q=0 is ignored
              disable_type_analysis,   &     ! If .true. disable the algorithm which detects whether
                                             ! there are atoms of the same type but with different occupations
              skip_atom(500)                 ! If .true. no LR calculation will be performed 

LR_Modules/adddvhubscf.f90

@@ -1,5 +1,5 @@
 !                                         
-! Copyright (C) 2001-2018 Quantum ESPRESSO
+! Copyright (C) 2001-2015 Quantum ESPRESSO
 ! This file is distributed under the terms
 ! GNU General Public License. See the file
 ! in the root directory of the present dis
@@ -208,3 +208,272 @@ SUBROUTINE adddvhubscf (ipert, ik)
   RETURN
   ! 
 END SUBROUTINE adddvhubscf
+
+SUBROUTINE revert_mag_u (ns_nc)
+   !
+   ! This routine reverts the sign of the Hubbard magnetization,
+   ! to be used in the noncollinear magnetic case.
+   ! Written by L. Binci (2023)
+   !
+   USE kinds,                ONLY : DP
+   USE ions_base,            ONLY : nat, ityp
+   USE lsda_mod,             ONLY : nspin
+   USE ldaU,                 ONLY : lda_plus_u_kind, Hubbard_l, is_hubbard, &
+                                    Hubbard_lmax
+   USE noncollin_module,     ONLY : npol
+   !
+   IMPLICIT NONE
+   INTEGER :: na, nt, m1, m2, is, is1, is2, ldim
+   COMPLEX(DP) :: ns, ms(3)
+   COMPLEX(DP), INTENT(INOUT) :: ns_nc(2*Hubbard_lmax+1, 2*Hubbard_lmax+1, nspin, nat)
+   !
+   DO na = 1, nat
+      nt = ityp(na)
+      IF (is_hubbard(nt)) THEN
+          ldim = 2 * Hubbard_l(nt) + 1
+          DO m1 = 1, ldim
+              DO m2 = 1, ldim
+                   !
+                   ! charge
+                   ns = (ns_nc(m1,m2,1,na) + ns_nc(m1,m2,4,na))
+                   ! magnetization
+                   ms(1) =  (ns_nc(m1,m2,2,na)  + ns_nc(m1,m2,3,na))
+                   ms(2) =  -(0.d0,1.0d0)*(ns_nc(m1,m2,2,na)  - ns_nc(m1,m2,3,na))
+                   ms(3) =  (ns_nc(m1,m2,1,na)  - ns_nc(m1,m2,4,na))
+                   !
+                   ms(:) = -ms(:)
+                   !
+                   ns_nc(m1,m2,1,na) = 0.5d0*( ns + ms(3))
+                   ns_nc(m1,m2,2,na) = 0.5d0*( ms(1) + (0.d0,1.0d0)*ms(2) )
+                   ns_nc(m1,m2,3,na) = 0.5d0*( ms(1) - (0.d0,1.0d0)*ms(2) )
+                   ns_nc(m1,m2,4,na) = 0.5d0*( ns - ms(3) )
+              ENDDO
+          ENDDO
+      ENDIF
+   ENDDO
+   !
+   RETURN
+   !
+ END SUBROUTINE revert_mag_u
+
+ SUBROUTINE calc_vh_u( ns, v_hub)
+   !
+   !! Recomputes the noncollinear Hubbard potential.
+   !! Used for a (spin-)Hubbard matrix with reversed magnetization
+   !! Similar to v_hubbard_nc within PW/src/v_of_rho.f90
+   !
+   USE kinds,                ONLY : DP
+   USE ions_base,            ONLY : nat, ityp
+   USE ldaU,                 ONLY : Hubbard_lmax, Hubbard_l, Hubbard_U
+   USE lsda_mod,             ONLY : nspin
+   USE control_flags,        ONLY : iverbosity, dfpt_hub
+   USE io_global,            ONLY : stdout
+   !
+   IMPLICIT NONE
+   !
+   COMPLEX(DP), INTENT(IN) :: ns(2*Hubbard_lmax+1,2*Hubbard_lmax+1,nspin,nat)
+   !! Occupation matrix
+   COMPLEX(DP), INTENT(INOUT) :: v_hub(2*Hubbard_lmax+1,2*Hubbard_lmax+1,nspin,nat)
+   !! Hubbard potential
+   INTEGER  :: is, is1, na, nt, m1, m2
+   ! ... local variables
+   !
+   v_hub(:,:,:,:) = 0.d0
+   !
+   DO na = 1, nat
+      nt = ityp (na)
+      IF (Hubbard_U(nt) /= 0.d0) THEN
+         !
+         ! is=1 and is=4 are diagonal components
+         ! is=2 and is=3 are off-diagonal components
+         DO is = 1, nspin
+            IF (is == 2) THEN
+             is1 = 3
+            ELSEIF (is == 3) THEN
+             is1 = 2
+            ELSE
+             is1 = is ! is=1 or is=4
+            ENDIF
+            !
+            IF (is1 == is) THEN
+               !
+               ! Non spin-flip contribution (is=1 and is=4)
+               ! (diagonal [spin indexes] occupancy matrices)
+               DO m1 = 1, 2*Hubbard_l(nt) + 1
+                  v_hub(m1,m1,is,na) = v_hub(m1,m1,is,na) + 0.5D0*Hubbard_U(nt)
+                  DO m2 = 1, 2 * Hubbard_l(nt) + 1
+                     v_hub(m1,m2,is,na) = v_hub(m1,m2,is,na) - Hubbard_U(nt) * ns(m2,m1,is,na)
+                  ENDDO
+               ENDDO
+            ELSE
+               !
+               ! Spin-flip contribution (is=2 and is=3)
+               ! (NON-diagonal [spin indexes] occupancy matrices)
+               DO m1 = 1, 2*Hubbard_l(nt) + 1
+                  DO m2 = 1, 2 * Hubbard_l(nt) + 1
+                     v_hub(m1,m2,is,na) = v_hub(m1,m2,is,na) - Hubbard_U(nt)*ns(m2,m1,is1,na)
+                  ENDDO
+               ENDDO
+            ENDIF
+         ENDDO !is
+      ENDIF
+   ENDDO
+   !
+   RETURN
+   !
+ END SUBROUTINE calc_vh_u
+ !--------------------------------------------------------------------------
+
+ SUBROUTINE revert_mag_uv (nsg_nc)
+   !
+   ! This routine reverts the sign of the Hubbard magnetization,
+   ! to be used in the noncollinear magnetic case.
+   ! Written by L. Binci (2023)
+   !
+   USE kinds,                ONLY : DP
+   USE ions_base,            ONLY : nat, ityp
+   USE lsda_mod,             ONLY : nspin
+   USE ldaU,                 ONLY : lda_plus_u_kind, Hubbard_l, is_hubbard, &
+                                    Hubbard_lmax, ldmx_tot, max_num_neighbors,ldim_u,&
+                                    neighood, at_sc
+   USE noncollin_module,     ONLY : npol
+   USE control_lr,           ONLY : nbnd_occ
+   !
+   IMPLICIT NONE
+   INTEGER :: na, nt, m1, m2, is, is1, is2, ldim, ibnd, nt1, nt2, na1, na2, viz,&
+              ldim1, ldim2, eq_na2, isp1, isp2
+   INTEGER, EXTERNAL :: find_viz
+   COMPLEX(DP) :: ns, ms(3)
+   COMPLEX(DP), INTENT(INOUT) :: nsg_nc(ldmx_tot, ldmx_tot, max_num_neighbors, nat, nspin)
+   !
+   DO na1 = 1, nat
+      nt1 = ityp (na1)
+      IF ( ldim_u(nt1).GT.0 ) THEN
+         ldim1 = ldim_u(nt1)
+         DO viz = 1, neighood(na1)%num_neigh
+            na2 = neighood(na1)%neigh(viz)
+            eq_na2 = at_sc(na2)%at
+            nt2 = ityp (eq_na2)
+            ldim2 = ldim_u(nt2)
+            !IF (na1.GT.na2) THEN
+            !   DO m1 = 1, ldim1
+            !      DO m2 = 1, ldim2
+            !         DO is1 = 1, npol
+            !            DO is2 = 1, npol
+            !               isp1 = is2 + npol*(is1-1)
+            !               isp2 = is1 + npol*(is2-1)
+            !               nsg_nc(m2,m1,viz,na1,isp1) = &
+            !                  CONJG(nsg_nc(m1,m2,find_viz(na2,na1),na2,isp2))
+            !            ENDDO
+            !         ENDDO
+            !      ENDDO
+            !   ENDDO
+            !ELSE
+               DO m1 = 1, ldim1
+                  DO m2 = 1, ldim2
+                     ! charge
+                     ns = (nsg_nc(m2,m1,viz,na1,1) + nsg_nc(m2,m1,viz,na1,4))
+                     ! magnetization
+                     ms(1) = (nsg_nc(m2,m1,viz,na1,2) + nsg_nc(m2,m1,viz,na1,3))
+                     ms(2) = -(0.d0,1.0d0)*(nsg_nc(m2,m1,viz,na1,2) - nsg_nc(m2,m1,viz,na1,3))
+                     ms(3) = (nsg_nc(m2,m1,viz,na1,1) - nsg_nc(m2,m1,viz,na1,4))
+                     !
+                     ms(:) = -ms(:)
+                     !
+                     nsg_nc(m2,m1,viz,na1,1) = 0.5d0*( ns + ms(3))
+                     nsg_nc(m2,m1,viz,na1,2) = 0.5d0*( ms(1) + (0.d0,1.0d0)*ms(2) )
+                     nsg_nc(m2,m1,viz,na1,3) = 0.5d0*( ms(1) - (0.d0,1.0d0)*ms(2) )
+                     nsg_nc(m2,m1,viz,na1,4) = 0.5d0*( ns - ms(3) )
+                  ENDDO
+               ENDDO
+            !ENDIF
+         ENDDO
+      ENDIF
+   ENDDO
+   !
+   RETURN
+   !
+ END SUBROUTINE revert_mag_uv
+
+ SUBROUTINE calc_vh_uv( nsg, v_hub)
+   !
+   !! Recomputes the noncollinear Hubbard potential.
+   !! Used for a (spin-)Hubbard matrix with reversed magnetization
+   !! Similar to v_hubbard_nc within PW/src/v_of_rho.f90
+   !
+   USE kinds,             ONLY : DP
+   USE ions_base,         ONLY : nat, ityp
+   USE ldaU,              ONLY : Hubbard_l, Hubbard_alpha, Hubbard_J0, Hubbard_beta,   &
+                                 ldim_u, ldmx_tot, max_num_neighbors, at_sc, neighood, &
+                                 Hubbard_V, Hubbard_alpha_back, is_hubbard, is_hubbard_back
+   USE lsda_mod,          ONLY : nspin
+   USE control_flags,     ONLY : iverbosity, dfpt_hub
+   USE io_global,         ONLY : stdout
+   USE noncollin_module,  ONLY : npol
+   !
+   IMPLICIT NONE
+   !
+   COMPLEX(DP), INTENT(IN)  :: nsg  (ldmx_tot, ldmx_tot, max_num_neighbors, nat, nspin)
+   COMPLEX(DP), INTENT(OUT) :: v_hub(ldmx_tot, ldmx_tot, max_num_neighbors, nat, nspin)
+   !
+   ! Local variables
+   !
+   INTEGER :: is, is1, isop, na, na1, na2, nt, nt1, nt2, m1, m2, viz, equiv_na2
+   INTEGER, EXTERNAL :: find_viz
+   !
+   v_hub(:,:,:,:,:) = (0.d0, 0.d0)
+   !
+   DO na1 = 1, nat
+      nt1 = ityp(na1)
+      IF ( is_hubbard(nt1) ) THEN
+         DO is = 1, nspin
+            IF (is == 2) THEN
+               is1 = 3
+            ELSEIF (is == 3) THEN
+               is1 = 2
+            ELSE
+               is1 = is
+            ENDIF
+            DO viz = 1, neighood(na1)%num_neigh
+               na2 = neighood(na1)%neigh(viz)
+               equiv_na2 = at_sc(na2)%at
+               nt2 = ityp(equiv_na2)
+               !
+               IF (is_hubbard(nt2) .AND. &
+                   (Hubbard_V(na1,na2,1).NE.0.d0) ) THEN
+                   !
+                   ! Here no need to use is1: complex conjugation is enough
+                   ! For both standard and background states of a center atom
+                   DO m1 = 1, ldim_u(nt1)
+                      ! For both standard and background states of the neighbor atom
+                      DO m2 = 1, ldim_u(nt2)
+                         v_hub(m2,m1,viz,na1,is) = - CONJG(nsg(m2,m1,viz,na1,is)) * Hubbard_V(na1,na2,1)
+                      ENDDO
+                   ENDDO
+                   !
+                   IF ( na1.EQ.na2 .AND. is1.EQ.is) THEN
+                      !
+                      na = find_viz(na1,na1)
+                      ! This is the diagonal term (like in the DFT+U only case)
+                      DO m1 = 1, ldim_u(nt1)
+                         v_hub(m1,m1,na,na1,is) = v_hub(m1,m1,na,na1,is) &
+                                        + Hubbard_V(na1,na1,1) * 0.5d0
+                      ENDDO
+                      !
+                   ENDIF
+                   !
+               ENDIF
+               !
+            ENDDO ! viz
+            !
+         ENDDO ! is
+         !
+      ENDIF
+      !
+   ENDDO ! na1
+   !
+   !
+   RETURN
+   !
+ END SUBROUTINE calc_vh_uv
+ !--------------------------------------------------------------------------

LR_Modules/lrcom.f90

@@ -1,5 +1,5 @@
 !
-! Copyright (C) 2001-2019 Quantum ESPRESSO group
+! Copyright (C) 2001-2025 Quantum ESPRESSO group
 ! This file is distributed under the terms of the
 ! GNU General Public License. See the file `License'
 ! in the root directory of the present distribution,
@@ -261,6 +261,11 @@ MODULE ldaU_lr
   COMPLEX(DP), ALLOCATABLE :: dnsscf(:,:,:,:,:)
   !! SCF derivative of ns
   !
+  COMPLEX(DP), ALLOCATABLE :: vh_u_save(:,:,:,:,:)
+  COMPLEX(DP), ALLOCATABLE :: vh_uv_save(:,:,:,:,:,:)
+  ! to save the two unperturbed Hubbard U and V potentials;
+  ! one normal, and the other with the time-reversed m_hubb
+  !
   COMPLEX(DP), ALLOCATABLE, TARGET :: swfcatomk(:,:)
   !! S * atomic wfc at k
   COMPLEX(DP), POINTER :: swfcatomkpq(:,:)