diff --git a/Doc/release-notes b/Doc/release-notes index bbe7cccea..dd8998810 100644 --- a/Doc/release-notes +++ b/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 diff --git a/HP/Doc/INPUT_HP.def b/HP/Doc/INPUT_HP.def index eebec0687..e3cdc4034 100644 --- a/HP/Doc/INPUT_HP.def +++ b/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. + } + } + } } diff --git a/HP/examples/example02/README b/HP/examples/example02/README index 4c3af8e7a..a17793044 100644 --- a/HP/examples/example02/README +++ b/HP/examples/example02/README @@ -1,5 +1,6 @@ -In order to compute U for magnetic insulators one needs to follow -special procedure. Since we do not know a priori the magnetization +There are two ways how to compute U for magnetic insulators. + +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)). diff --git a/HP/examples/example02/run_example b/HP/examples/example02/run_example1 similarity index 100% rename from HP/examples/example02/run_example rename to HP/examples/example02/run_example1 diff --git a/HP/examples/example02/run_example2 b/HP/examples/example02/run_example2 new file mode 100755 index 000000000..c3a21d881 --- /dev/null +++ b/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" diff --git a/HP/src/hp_allocate_q.f90 b/HP/src/hp_allocate_q.f90 index 038571159..0c7717176 100644 --- a/HP/src/hp_allocate_q.f90 +++ b/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_lr, ONLY : swfcatomk, swfcatomkpq + USE ldaU, ONLY : Hubbard_lmax, nwfcU, lda_plus_u_kind, max_num_neighbors + 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 diff --git a/HP/src/hp_bcast_input.f90 b/HP/src/hp_bcast_input.f90 index 722c3d9c0..e69b95dcb 100644 --- a/HP/src/hp_bcast_input.f90 +++ b/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 ! diff --git a/HP/src/hp_check_type.f90 b/HP/src/hp_check_type.f90 index 988a59cb1..139909e11 100644 --- a/HP/src/hp_check_type.f90 +++ b/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 diff --git a/HP/src/hp_dealloc_q.f90 b/HP/src/hp_dealloc_q.f90 index af06079cf..0d52dd237 100644 --- a/HP/src/hp_dealloc_q.f90 +++ b/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) diff --git a/HP/src/hp_readin.f90 b/HP/src/hp_readin.f90 index 47a508a5a..6e24f1177 100644 --- a/HP/src/hp_readin.f90 +++ b/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 = './' ! diff --git a/HP/src/hp_setup_q.f90 b/HP/src/hp_setup_q.f90 index 935aa9d53..1fff574d8 100644 --- a/HP/src/hp_setup_q.f90 +++ b/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) ! diff --git a/HP/src/hp_solve_linear_system.f90 b/HP/src/hp_solve_linear_system.f90 index fedaf9345..b66102683 100644 --- a/HP/src/hp_solve_linear_system.f90 +++ b/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 diff --git a/HP/src/hp_sym_dmag.f90 b/HP/src/hp_sym_dmag.f90 index 794799044..7b61d9fbd 100644 --- a/HP/src/hp_sym_dmag.f90 +++ b/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) diff --git a/HP/src/hpcom.f90 b/HP/src/hpcom.f90 index 1f9977a19..bd6daa5a9 100644 --- a/HP/src/hpcom.f90 +++ b/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 diff --git a/LR_Modules/adddvhubscf.f90 b/LR_Modules/adddvhubscf.f90 index 87e5dac40..4ffdf0a6e 100644 --- a/LR_Modules/adddvhubscf.f90 +++ b/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 + !-------------------------------------------------------------------------- diff --git a/LR_Modules/lrcom.f90 b/LR_Modules/lrcom.f90 index a3ed4075e..d92c9844c 100644 --- a/LR_Modules/lrcom.f90 +++ b/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(:,:)