From a4435247fe41e18cd58e86ce6552df10c9e07c52 Mon Sep 17 00:00:00 2001
From: giannozz
Date: Fri, 3 May 2019 16:15:01 +0200
Subject: [PATCH 1/4] Variable "npw" made explicitly local in a few place.
Reminder: global variable "npw" must not be use,d I want to get rid of it
---
HP/src/hp_efermi_shift.f90 | 1 -
PP/src/ppacf.f90 | 3 ++-
PW/src/stres_mgga.f90 | 4 ++--
3 files changed, 4 insertions(+), 4 deletions(-)
diff --git a/HP/src/hp_efermi_shift.f90 b/HP/src/hp_efermi_shift.f90
index 0a1b52bcb..da74e6d63 100644
--- a/HP/src/hp_efermi_shift.f90
+++ b/HP/src/hp_efermi_shift.f90
@@ -31,7 +31,6 @@ SUBROUTINE hp_ef_shift (drhoscf, ldos, ldoss, dos_ef, dbecsum, becsum1)
USE gvect, ONLY : gg
USE buffers, ONLY : get_buffer, save_buffer
USE lsda_mod, ONLY : nspin
- USE wvfct, ONLY : npw, npwx, et
USE klist, ONLY : degauss, ngauss, ngk
USE noncollin_module, ONLY : noncolin, npol, nspin_mag, nspin_lsda
USE mp_bands, ONLY : intra_bgrp_comm
diff --git a/PP/src/ppacf.f90 b/PP/src/ppacf.f90
index a55d6649f..1e4141c1b 100644
--- a/PP/src/ppacf.f90
+++ b/PP/src/ppacf.f90
@@ -51,7 +51,7 @@ PROGRAM do_ppacf
USE funct, ONLY : xc,xc_spin,gcxc,gcx_spin,gcc_spin,dft_is_nonlocc,nlc
USE funct, ONLY : get_iexch, get_icorr, get_igcx, get_igcc
USE funct, ONLY : set_exx_fraction,set_auxiliary_flags,enforce_input_dft
- USE wvfct, ONLY : npw, npwx
+ USE wvfct, ONLY : npwx
USE environment, ONLY : environment_start, environment_end
USE kernel_table, ONLY : Nqs, vdw_table_name, kernel_file_name
USE vdW_DF, ONLY : get_potential, vdW_energy
@@ -110,6 +110,7 @@ PROGRAM do_ppacf
REAL(DP), ALLOCATABLE :: tot_grad_rho(:,:),grad_rho(:,:,:)
REAL(DP), ALLOCATABLE :: tot_rho(:)
+ INTEGER :: npw ! number of plane waves
INTEGER :: ik ! counter on k points
INTEGER, ALLOCATABLE :: igk_buf(:)
REAL(dp), ALLOCATABLE :: gk(:) ! work space
diff --git a/PW/src/stres_mgga.f90 b/PW/src/stres_mgga.f90
index 3334bba27..4a964c834 100644
--- a/PW/src/stres_mgga.f90
+++ b/PW/src/stres_mgga.f90
@@ -23,7 +23,7 @@ SUBROUTINE stres_mgga( sigmaxc )
USE klist, ONLY : nks, xk, ngk
USE buffers, ONLY : get_buffer
USE io_files, ONLY : iunwfc, nwordwfc
- USE wvfct, ONLY : nbnd, npwx, npw, wg
+ USE wvfct, ONLY : nbnd, npwx, wg
USE lsda_mod, ONLY : lsda, nspin, current_spin, isk
USE fft_interfaces, ONLY : fwfft, invfft
USE fft_base, ONLY : dfftp, dffts
@@ -37,7 +37,7 @@ SUBROUTINE stres_mgga( sigmaxc )
!
! Internal variables
!
- INTEGER :: ix, iy, ir, ipol, iss, incr, ibnd, ik
+ INTEGER :: ix, iy, ir, ipol, iss, incr, ibnd, ik, npw
INTEGER :: ipol2xy(3,3)
REAL(DP), PARAMETER :: epsr = 1.0d-6, epsg = 1.0d-10, e2 = 2.d0
COMPLEX(DP), ALLOCATABLE :: gradwfc (:,:), crosstaus(:,:,:)
From 32fa25445dc68a98f565793ee20bea9d614c0129 Mon Sep 17 00:00:00 2001
From: Oleksandr Motornyi
Date: Tue, 12 Mar 2019 16:34:27 +0000
Subject: [PATCH 2/4] spin-orbit coupling with uspp in turbo_eels code
done by: Oleksandr Motornyi, Andrea Dal Corso, Nathalie Vast.
---
Doc/release-notes | 3 +
LR_Modules/lr_sm1_psi.f90 | 977 ++++++-----------------
LR_Modules/lrcom.f90 | 1 +
TDDFPT/src/Makefile | 4 +
TDDFPT/src/compute_intq.f90 | 89 +++
TDDFPT/src/dveqpsi_us_only.f90 | 110 +++
TDDFPT/src/lr_addus_dvpsi.f90 | 185 ++---
TDDFPT/src/lr_alloc_init.f90 | 7 +-
TDDFPT/src/lr_apply_liouvillian.f90 | 4 +-
TDDFPT/src/lr_apply_liouvillian_eels.f90 | 2 +-
TDDFPT/src/lr_dealloc.f90 | 3 +-
TDDFPT/src/lr_dvpsi_e.f90 | 4 +-
TDDFPT/src/lr_dvpsi_eels.f90 | 5 +-
TDDFPT/src/lr_readin.f90 | 2 -
TDDFPT/src/lr_restart.f90 | 2 +
TDDFPT/src/lr_solve_e.f90 | 19 +-
TDDFPT/src/lr_variables.f90 | 7 +-
TDDFPT/src/make.depend | 21 +
TDDFPT/src/set_intq_nc.f90 | 35 +
TDDFPT/src/transform_intq_nc.f90 | 39 +
TDDFPT/src/transform_intq_so.f90 | 58 ++
21 files changed, 738 insertions(+), 839 deletions(-)
create mode 100644 TDDFPT/src/compute_intq.f90
create mode 100644 TDDFPT/src/dveqpsi_us_only.f90
create mode 100644 TDDFPT/src/set_intq_nc.f90
create mode 100644 TDDFPT/src/transform_intq_nc.f90
create mode 100644 TDDFPT/src/transform_intq_so.f90
diff --git a/Doc/release-notes b/Doc/release-notes
index 2406f1e56..50c92d9e3 100644
--- a/Doc/release-notes
+++ b/Doc/release-notes
@@ -1,3 +1,6 @@
+New in development branch:
+ * turbo_eels code of TDDFPT module now works with ultrasoft pseudopotentials and spin-orbit coupling together (Oleksandr Motornyi, Andrea Dal Corso, Nathalie Vast). lr_sm1_psi.f90 of LR_Modules is rewritten and simplified.
+
Problems fixed in development branch :
* Bug in spin-polarized meta-GGA (noticed by Shoaib Muhammad,
Sungkyunkwan U.)
diff --git a/LR_Modules/lr_sm1_psi.f90 b/LR_Modules/lr_sm1_psi.f90
index 0a8923b77..d4edaead9 100644
--- a/LR_Modules/lr_sm1_psi.f90
+++ b/LR_Modules/lr_sm1_psi.f90
@@ -6,20 +6,14 @@
! or http://www.gnu.org/copyleft/gpl.txt .
!
!----------------------------------------------------------------------------
-SUBROUTINE lr_sm1_psi (recalculate, ik, lda, n, m, psi, spsi)
+SUBROUTINE lr_sm1_psi_tpw (ik, lda, n, m, psi, spsi)
!----------------------------------------------------------------------------
!
! This subroutine applies the S^{-1} matrix to m wavefunctions psi
! and puts the results in spsi.
! See Eq.(13) in B. Walker and R. Gebauer, JCP 127, 164106 (2007).
- ! Requires the products of psi with all beta functions
- ! in array becp(nkb,m) (calculated in h_psi or by ccalbec).
!
- ! INPUT: recalculate Decides if the overlap of beta
- ! functions is recalculated or not.
- ! This is needed e.g. if ions are moved
- ! and the overlap changes accordingly.
- ! ik k point under consideration
+ ! INPUT: ik k point under consideration
! lda leading dimension of arrays psi, spsi
! n true dimension of psi, spsi
! m number of states psi
@@ -30,38 +24,36 @@ SUBROUTINE lr_sm1_psi (recalculate, ik, lda, n, m, psi, spsi)
! Original routine written by R. Gebauer
! Modified by Osman Baris Malcioglu (2009)
! Modified by Iurii Timrov (2013)
+ ! Simplified and generalized to the relativistic case by Andrea Dal Corso (2018)
!
USE kinds, ONLY : DP
USE control_flags, ONLY : gamma_only
+ USE klist, ONLY : xk, igk_k, ngk
+ USE qpoint, ONLY : ikks, ikqs, nksq
USE uspp, ONLY : okvan, vkb, nkb, qq_nt
USE uspp_param, ONLY : nh, upf
- USE ions_base, ONLY : ityp,nat,ntyp=>nsp
+ USE ions_base, ONLY : ityp, nat, ntyp=>nsp
+ USE becmod, ONLY : bec_type, becp, calbec
USE mp, ONLY : mp_sum
- USE mp_bands, ONLY : intra_bgrp_comm
- USE noncollin_module, ONLY : noncolin, npol
- USE matrix_inversion
+ USE mp_global, ONLY : intra_bgrp_comm
+ USE noncollin_module, ONLY : noncolin, npol, nspin_mag
!
IMPLICIT NONE
- LOGICAL, INTENT(in) :: recalculate
- INTEGER, INTENT(in) :: lda, n, m, ik
+ INTEGER, INTENT(in) :: ik, lda,n,m
COMPLEX(DP), INTENT(in) :: psi(lda*npol,m)
COMPLEX(DP), INTENT(out) :: spsi(lda*npol,m)
!
- LOGICAL :: recalc
- !
- CALL start_clock( 'lr_sm1_psi' )
- !
- recalc = recalculate
+ CALL start_clock( 'lr_sm1_psi_tpw' )
!
IF ( gamma_only ) THEN
CALL sm1_psi_gamma()
ELSEIF (noncolin) THEN
- CALL errore( 'lr_sm1_psi', 'Noncollinear case is not implemented', 1 )
+ CALL sm1_psi_nc()
ELSE
CALL sm1_psi_k()
ENDIF
!
- CALL stop_clock( 'lr_sm1_psi' )
+ CALL stop_clock( 'lr_sm1_psi_tpw' )
!
RETURN
!
@@ -77,8 +69,9 @@ CONTAINS
! outside of this routine.
!
USE becmod, ONLY : bec_type,becp,calbec
- USE realus, ONLY : real_space, invfft_orbital_gamma, initialisation_level, &
- fwfft_orbital_gamma, calbec_rs_gamma, add_vuspsir_gamma, &
+ USE realus, ONLY : real_space, invfft_orbital_gamma, &
+ initialisation_level, fwfft_orbital_gamma, &
+ calbec_rs_gamma, add_vuspsir_gamma, &
v_loc_psir, s_psir_gamma
USE lrus, ONLY : bbg
!
@@ -86,10 +79,9 @@ CONTAINS
!
! ... local variables
!
- INTEGER :: ikb, jkb, ih, jh, na, nt, ijkb0, ibnd, ii
+ INTEGER :: ibnd
! counters
REAL(DP), ALLOCATABLE :: ps(:,:)
- LOGICAL, SAVE :: first_entry = .true.
!
! Initialize spsi : spsi = psi
!
@@ -97,115 +89,7 @@ CONTAINS
!
IF ( nkb == 0 .OR. .NOT. okvan ) RETURN
!
- ! If this is the first entry, we calculate and save the coefficients B from Eq.(15)
- ! B. Walker and R. Gebauer, J. Chem. Phys. 127, 164106 (2007).
- ! If this is not the first entry, we do not recalculate the coefficients B but
- ! use the ones which were already calculated and saved (if recalc=.false.).
- !
- IF (first_entry) THEN
- !
- IF (allocated(bbg)) DEALLOCATE(bbg)
- first_entry = .false.
- recalc = .true.
- !
- ENDIF
- !
- IF (.not.allocated(bbg)) recalc = .true.
- IF (recalc .and. allocated(bbg)) DEALLOCATE(bbg)
- !
- IF (recalc) THEN
- !
- ALLOCATE(bbg(nkb,nkb))
- bbg = 0.d0
- !
- ! The beta-functions vkb must have beed calculated elsewhere.
- ! So there is no need to call the routine init_us_2.
- !
- ! Calculate the coefficients B_ij defined by Eq.(15).
- ! B_ij = , where beta(i) = vkb(i).
- !
- CALL calbec (n, vkb, vkb, bbg(:,:), nkb)
- !
- ALLOCATE(ps(nkb,nkb))
- !
- ps(:,:) = 0.d0
- !
- ! Calculate the product of q_nm and B_ij of Eq.(16).
- !
- ijkb0 = 0
- DO nt=1,ntyp
- IF (upf(nt)%tvanp) THEN
- DO na=1,nat
- IF(ityp(na)==nt) THEN
- DO ii=1,nkb
- DO jh=1,nh(nt)
- jkb=ijkb0 + jh
- DO ih=1,nh(nt)
- ikb = ijkb0 + ih
- ps(ikb,ii) = ps(ikb,ii) + &
- & qq_nt(ih,jh,nt) * bbg(jkb,ii)
- ENDDO
- ENDDO
- ENDDO
- ijkb0 = ijkb0+nh(nt)
- ENDIF
- ENDDO
- ELSE
- DO na = 1, nat
- IF ( ityp(na) == nt ) ijkb0 = ijkb0 + nh(nt)
- ENDDO
- ENDIF
- ENDDO
- !
- ! Add identity to q_nm * B_ij [see Eq.(16)].
- ! ps = (1 + q*B)
- !
- DO ii=1,nkb
- ps(ii,ii) = ps(ii,ii) + 1.d0
- ENDDO
- !
- ! Invert matrix: (1 + q*B)^{-1}
- !
- CALL invmat( nkb, ps )
- !
- ! Use the array bbg as a workspace in order to save memory
- !
- bbg(:,:) = 0.d0
- !
- ! Finally, let us calculate lambda_nm = -(1+q*B)^{-1} * q
- !
- ijkb0 = 0
- DO nt=1,ntyp
- IF (upf(nt)%tvanp) THEN
- DO na=1,nat
- IF(ityp(na)==nt) THEN
- DO ii=1,nkb
- DO jh=1,nh(nt)
- jkb=ijkb0 + jh
- DO ih=1,nh(nt)
- ikb = ijkb0 + ih
- bbg(ii,jkb) = bbg(ii,jkb) &
- & - ps(ii,ikb) * qq_nt(ih,jh,nt)
- ENDDO
- ENDDO
- ENDDO
- ijkb0 = ijkb0+nh(nt)
- ENDIF
- ENDDO
- ELSE
- DO na = 1, nat
- IF ( ityp(na) == nt ) ijkb0 = ijkb0 + nh(nt)
- ENDDO
- ENDIF
- ENDDO
- DEALLOCATE(ps)
- ENDIF
- !
- ! Compute the product of the beta-functions vkb with the functions psi,
- ! and put the result in becp.
- ! becp(ikb,jbnd) = \sum_G vkb^*(ikb,G) psi(G,jbnd) =
- !
- IF (real_space) THEN
+ IF (real_space) THEN
!
DO ibnd=1,m,2
CALL invfft_orbital_gamma(psi,ibnd,m)
@@ -213,7 +97,7 @@ CONTAINS
ENDDO
!
ELSE
- CALL calbec(n,vkb,psi,becp,m)
+ CALL calbec(n,vkb,psi,becp,m)
ENDIF
!
! Use the array ps as a workspace
@@ -221,7 +105,7 @@ CONTAINS
ps(:,:) = 0.D0
!
! Now let us apply the operator S^{-1}, given by Eq.(13), to the functions psi.
- ! Let's do this in 2 steps.
+ ! Let's DO this in 2 steps.
!
! Step 1 : ps = lambda *
! Here, lambda = bbg, and = becp%r
@@ -237,255 +121,8 @@ CONTAINS
RETURN
!
END SUBROUTINE sm1_psi_gamma
-
+
SUBROUTINE sm1_psi_k()
- !-----------------------------------------------------------------------
- !
- ! k-points version
- ! Note: the array bbk must be deallocated somewhere
- ! outside of this routine.
- !
- USE becmod, ONLY : bec_type,becp,calbec
- USE klist, ONLY : nks, xk, ngk, igk_k
- USE lrus, ONLY : bbk
- !
- IMPLICIT NONE
- !
- ! ... local variables
- !
- INTEGER :: ikb, jkb, ih, jh, na, nt, ijkb0, ibnd, ii, ik1
- ! counters
- COMPLEX(DP), ALLOCATABLE :: ps(:,:)
- !
- ! Initialize spsi : spsi = psi
- !
- CALL ZCOPY( lda*npol*m, psi, 1, spsi, 1 )
- !
- IF ( nkb == 0 .OR. .NOT. okvan ) RETURN
- !
- ! If this is the first entry, we calculate and save the coefficients B from Eq.(15)
- ! B. Walker and R. Gebauer, J. Chem. Phys. 127, 164106 (2007).
- ! If this is not the first entry, we do not recalculate the coefficients B
- ! (they are already allocated) but use the ones which were already calculated
- ! and saved (if recalc=.false.).
- !
- IF (.NOT.ALLOCATED(bbk)) recalc = .true.
- IF (recalc .AND. ALLOCATED(bbk)) DEALLOCATE(bbk)
- !
- ! If recalc=.true. we (re)calculate the coefficients lambda_nm defined by Eq.(16).
- !
- IF (recalc) THEN
- !
- ALLOCATE(bbk(nkb,nkb,nks))
- bbk = (0.d0,0.d0)
- !
- ALLOCATE(ps(nkb,nkb))
- !
- DO ik1 = 1, nks
- !
- ! Calculate beta-functions vkb for a given k point.
- !
- CALL init_us_2(ngk(ik1),igk_k(:,ik1),xk(1,ik1),vkb)
- !
- ! Calculate the coefficients B_ij defined by Eq.(15).
- ! B_ij = , where beta(i) = vkb(i).
- !
- CALL zgemm('C','N',nkb,nkb,ngk(ik1),(1.d0,0.d0),vkb, &
- & lda,vkb,lda,(0.d0,0.d0),bbk(1,1,ik1),nkb)
- !
-#if defined(__MPI)
- CALL mp_sum(bbk(:,:,ik1), intra_bgrp_comm)
-#endif
- !
- ps(:,:) = (0.d0,0.d0)
- !
- ! Calculate the product of q_nm and B_ij of Eq.(16).
- !
- ijkb0 = 0
- DO nt=1,ntyp
- IF (upf(nt)%tvanp) THEN
- DO na=1,nat
- IF(ityp(na)==nt) THEN
- DO ii=1,nkb
- DO jh=1,nh(nt)
- jkb=ijkb0 + jh
- DO ih=1,nh(nt)
- ikb = ijkb0 + ih
- ps(ikb,ii) = ps(ikb,ii) + &
- & bbk(jkb,ii, ik1) * qq_nt(ih,jh,nt)
- ENDDO
- ENDDO
- ENDDO
- ijkb0 = ijkb0+nh(nt)
- ENDIF
- ENDDO
- ELSE
- DO na = 1, nat
- IF ( ityp(na) == nt ) ijkb0 = ijkb0 + nh(nt)
- ENDDO
- ENDIF
- ENDDO
- !
- ! Add an identity to q_nm * B_ij [see Eq.(16)].
- ! ps = (1 + q*B)
- !
- DO ii = 1, nkb
- ps(ii,ii) = ps(ii,ii) + (1.d0,0.d0)
- ENDDO
- !
- ! Invert matrix: (1 + q*B)^{-1}
- !
- CALL invmat( nkb, ps )
- !
- ! Use the array bbk as a work space in order to save memory.
- !
- bbk(:,:,ik1) = (0.d0,0.d0)
- !
- ! Finally, let us calculate lambda_nm = -(1+q*B)^{-1} * q
- !
- ijkb0 = 0
- DO nt=1,ntyp
- IF (upf(nt)%tvanp) THEN
- DO na=1,nat
- IF(ityp(na)==nt) THEN
- DO ii=1,nkb
- DO jh=1,nh(nt)
- jkb=ijkb0 + jh
- DO ih=1,nh(nt)
- ikb = ijkb0 + ih
- bbk(ii,jkb,ik1) = bbk(ii,jkb,ik1) &
- & - ps(ii,ikb) * qq_nt(ih,jh,nt)
- ENDDO
- ENDDO
- ENDDO
- ijkb0 = ijkb0+nh(nt)
- ENDIF
- ENDDO
- ELSE
- DO na = 1, nat
- IF ( ityp(na) == nt ) ijkb0 = ijkb0 + nh(nt)
- ENDDO
- ENDIF
- ENDDO
- !
- ENDDO ! loop on k points
- DEALLOCATE(ps)
- !
- ENDIF
- !
- IF (n.NE.ngk(ik)) CALL errore( 'sm1_psi_k', &
- & 'Mismatch in the number of plane waves', 1 )
- !
- ! Calculate beta-functions vkb for a given k point 'ik'.
- !
- CALL init_us_2(n,igk_k(:,ik),xk(1,ik),vkb)
- !
- ! Compute the product of the beta-functions vkb with the functions psi
- ! at point k, and put the result in becp%k.
- ! becp%k(ikb,jbnd) = \sum_G vkb^*(ikb,G) psi(G,jbnd) =
- !
- CALL calbec(n,vkb,psi,becp,m)
- !
- ! Use ps as a work space.
- !
- ALLOCATE(ps(nkb,m))
- ps(:,:) = (0.d0,0.d0)
- !
- ! Apply the operator S^{-1}, given by Eq.(13), to the functions psi.
- ! Let's do this in 2 steps.
- !
- ! Step 1 : calculate the product
- ! ps = lambda *
- !
- DO ibnd=1,m
- DO jkb=1,nkb
- DO ii=1,nkb
- ps(jkb,ibnd) = ps(jkb,ibnd) + bbk(jkb,ii,ik) * becp%k(ii,ibnd)
- ENDDO
- ENDDO
- ENDDO
- !
- ! Step 2 : |spsi> = S^{-1} * |psi> = |psi> + ps * |beta>
- !
- CALL ZGEMM( 'N', 'N', n, m, nkb, (1.D0, 0.D0), vkb, &
- & lda, ps, nkb, (1.D0, 0.D0), spsi, lda )
- !
- DEALLOCATE(ps)
- !
- RETURN
- !
- END SUBROUTINE sm1_psi_k
-
-END SUBROUTINE lr_sm1_psi
-
-
-!----------------------------------------------------------------------------
-SUBROUTINE lr_sm1_psiq (recalculate, ik, lda, n, m, psi, spsi)
- !----------------------------------------------------------------------------
- !
- ! This subroutine applies the S^{-1} matrix to m wavefunctions psi
- ! and puts the results in spsi.
- ! See Eq.(13) in B. Walker and R. Gebauer, JCP 127, 164106 (2007).
- ! Requires the products of psi with all beta functions
- ! in array becp(nkb,m) (calculated in h_psi or by ccalbec)
- !
- ! INPUT: recalculate Decides if the overlap of beta
- ! functions is recalculated or not.
- ! This is needed e.g. if ions are moved
- ! and the overlap changes accordingly.
- ! ik k point under consideration
- ! lda leading dimension of arrays psi, spsi
- ! n true dimension of psi, spsi
- ! m number of states psi
- ! psi the wavefunction to which the S^{-1}
- ! matrix is applied
- ! OUTPUT: spsi = S^{-1}*psi
- !
- ! Written by Iurii Timrov (2016) on the basis of lr_sm1_psi.
- ! Note: The difference of this routine from lr_sm1_psi is that
- ! it can be used when there is a perturbation with the finite
- ! (transferred) momentum q (e.g. by PHonon and turboEELS codes).
- !
- USE kinds, ONLY : DP
- USE control_flags, ONLY : gamma_only
- USE klist, ONLY : xk, igk_k, ngk
- USE qpoint, ONLY : ikks, ikqs, nksq
- USE uspp, ONLY : okvan, vkb, nkb, qq_nt
- USE uspp_param, ONLY : nh, upf
- USE ions_base, ONLY : ityp,nat,ntyp=>nsp
- USE becmod, ONLY : bec_type, becp, calbec
- USE mp, ONLY : mp_sum
- USE mp_bands, ONLY : intra_bgrp_comm
- USE noncollin_module, ONLY : noncolin, npol, nspin_mag
- USE matrix_inversion
- !
- IMPLICIT NONE
- LOGICAL, INTENT(in) :: recalculate
- INTEGER, INTENT(in) :: lda, n, m, ik
- COMPLEX(DP), INTENT(in) :: psi(lda*npol,m)
- COMPLEX(DP), INTENT(out) :: spsi(lda*npol,m)
- !
- LOGICAL :: recalc
- !
- CALL start_clock( 'lr_sm1_psiq' )
- !
- recalc = recalculate
- !
- IF ( gamma_only ) THEN
- CALL errore( 'lr_sm1_psiq', 'gamma_only is not supported', 1 )
- ELSEIF (noncolin) THEN
- CALL sm1_psiq_nc()
- ELSE
- CALL sm1_psiq_k()
- ENDIF
- !
- CALL stop_clock( 'lr_sm1_psiq' )
- !
- RETURN
- !
-CONTAINS
- !
- SUBROUTINE sm1_psiq_k()
!-----------------------------------------------------------------------
!
! k-points version
@@ -498,7 +135,7 @@ CONTAINS
!
! ... local variables
!
- INTEGER :: ikb, jkb, ih, jh, na, nt, ijkb0, ibnd, ii
+ INTEGER :: na, nt, ibnd, ii, jkb
INTEGER :: ik1, & ! dummy index for k points
ikk, & ! index of the point k
ikq, & ! index of the point k+q
@@ -511,131 +148,6 @@ CONTAINS
!
IF ( nkb == 0 .OR. .NOT. okvan ) RETURN
!
- ! If this is the first entry, we calculate and save the coefficients B from Eq.(15)
- ! B. Walker and R. Gebauer, J. Chem. Phys. 127, 164106 (2007).
- ! If this is not the first entry, we do not recalculate the coefficients B
- ! (they are already allocated) but use the ones which were already
- ! calculated and saved (if recalc=.false.).
- !
- IF (.NOT.ALLOCATED(bbk)) recalc = .true.
- IF (recalc .AND. ALLOCATED(bbk)) DEALLOCATE(bbk)
- !
- ! If recalc=.true. we (re)calculate the coefficients lambda_nm defined by Eq.(16).
- !
- IF ( recalc ) THEN
- !
- ALLOCATE(bbk(nkb,nkb,nksq))
- bbk = (0.0d0,0.0d0)
- !
- ALLOCATE(ps(nkb,nkb))
- !
- DO ik1 = 1, nksq
- !
- ikk = ikks(ik1)
- ikq = ikqs(ik1)
- npwq = ngk(ikq)
- !
- ! Calculate beta-functions vkb for a given k+q point.
- !
- CALL init_us_2 (npwq, igk_k(1,ikq), xk(1,ikq), vkb)
- !
- ! Calculate the coefficients B_ij defined by Eq.(15).
- ! B_ij = , where beta(i) = vkb(i).
- !
- CALL zgemm('C','N',nkb,nkb,npwq,(1.d0,0.d0),vkb, &
- & lda,vkb,lda,(0.d0,0.d0),bbk(1,1,ik1),nkb)
- !
-#if defined(__MPI)
- CALL mp_sum(bbk(:,:,ik1), intra_bgrp_comm)
-#endif
- !
- ps(:,:) = (0.d0,0.d0)
- !
- ! Calculate the product of q_nm and B_ij of Eq.(16).
- !
- ijkb0 = 0
- do nt=1,ntyp
- if (upf(nt)%tvanp) then
- do na=1,nat
- if(ityp(na).eq.nt) then
- do ii=1,nkb
- do jh=1,nh(nt)
- !
- jkb=ijkb0 + jh
- !
- do ih=1,nh(nt)
- !
- ikb = ijkb0 + ih
- !
- ps(ikb,ii) = ps(ikb,ii) + bbk(jkb,ii,ik1)*qq_nt(ih,jh,nt)
- !
- enddo
- enddo
- enddo
- ijkb0 = ijkb0+nh(nt)
- endif
- enddo
- else
- DO na = 1, nat
- IF ( ityp(na) == nt ) ijkb0 = ijkb0 + nh(nt)
- END DO
- endif
- enddo
- !
- ! Add an identity to q_nm * B_ij [see Eq.(16)].
- ! ps = (1 + q*B)
- !
- DO ii=1,nkb
- ps(ii,ii) = ps(ii,ii) + (1.d0,0.d0)
- ENDDO
- !
- ! Invert matrix: (1 + q*B)^{-1}
- !
- CALL invmat( nkb, ps )
- !
- ! Use the array bbk as a work space in order to save the memory.
- !
- bbk(:,:,ik1) = (0.d0,0.d0)
- !
- ! Finally, let us calculate lambda_nm = -(1+q*B)^{-1} * q
- ! Let us use the array bbk and put there the values of the lambda
- ! coefficients.
- !
- ijkb0 = 0
- do nt=1,ntyp
- if (upf(nt)%tvanp) then
- do na=1,nat
- if(ityp(na).eq.nt) then
- do ii=1,nkb
- do jh=1,nh(nt)
- !
- jkb=ijkb0 + jh
- !
- do ih=1,nh(nt)
- !
- ikb = ijkb0 + ih
- !
- bbk(ii,jkb,ik1) = bbk(ii,jkb,ik1) - ps(ii,ikb) * qq_nt(ih,jh,nt)
- !
- enddo
- enddo
- enddo
- ijkb0 = ijkb0+nh(nt)
- endif
- enddo
- else
- DO na = 1, nat
- IF ( ityp(na) == nt ) ijkb0 = ijkb0 + nh(nt)
- END DO
- endif
- enddo
- !
- ENDDO ! loop on k points
- !
- DEALLOCATE(ps)
- !
- ENDIF
- !
! Now set up the indices ikk and ikq such that they
! correspond to the points k and k+q using the index ik,
! which was passed to this routine as an input.
@@ -662,52 +174,44 @@ CONTAINS
ps(:,:) = (0.d0,0.d0)
!
! Apply the operator S^{-1}, given by Eq.(13), to the functions psi.
- ! Let's do this in 2 steps.
+ ! Let's DO this in 2 steps.
!
! Step 1 : calculate the product
! ps = lambda *
!
- DO ibnd=1,m
- DO jkb=1,nkb
- DO ii=1,nkb
- ps(jkb,ibnd) = ps(jkb,ibnd) + bbk(jkb,ii,ik) * becp%k(ii,ibnd)
- ENDDO
- ENDDO
- ENDDO
+ CALL ZGEMM( 'N', 'N', nkb, m, nkb, (1.D0, 0.D0), bbk(1,1,ik), &
+ nkb, becp%k, nkb, (1.D0, 0.D0), ps, nkb )
!
! Step 2 : |spsi> = S^{-1} * |psi> = |psi> + ps * |beta>
!
- CALL zgemm( 'N', 'N', n, m, nkb, (1.D0, 0.D0), vkb, &
- & lda, ps, nkb, (1.D0, 0.D0), spsi, lda )
+ CALL ZGEMM( 'N', 'N', n, m, nkb, (1.D0, 0.D0), vkb, &
+ lda, ps, nkb, (1.D0, 0.D0), spsi, lda )
!
DEALLOCATE(ps)
!
RETURN
!
-END SUBROUTINE sm1_psiq_k
+END SUBROUTINE sm1_psi_k
-SUBROUTINE sm1_psiq_nc()
+SUBROUTINE sm1_psi_nc()
!-----------------------------------------------------------------------
!
! Noncollinear case
- ! Note: 1) the implementation of this routine is not finished...
- ! 2) the array bbnc must be deallocated somewhere
- ! outside of this routine.
!
USE uspp, ONLY : qq_so
+ USE lrus, ONLY : bbnc_sm1
USE spin_orb, ONLY : lspinorb
- USE lrus, ONLY : bbnc
!
IMPLICIT NONE
!
! ... local variables
!
- INTEGER :: ikb, jkb, ih, jh, na, nt, ijkb0, ibnd, ii, ipol
+ INTEGER :: ibnd, ii, jkb, ipol, jpol, iis, jkbs
INTEGER :: ik1, & ! dummy index for k points
ikk, & ! index of the point k
ikq, & ! index of the point k+q
npwq ! number of the plane-waves at point k+q
- COMPLEX(DP), ALLOCATABLE :: ps(:,:,:)
+ COMPLEX(DP), ALLOCATABLE :: ps(:,:)
!
! Initialize spsi : spsi = psi
!
@@ -715,175 +219,6 @@ SUBROUTINE sm1_psiq_nc()
!
IF ( nkb == 0 .OR. .NOT. okvan ) RETURN
!
- CALL errore( 'sm1_psiq_nc', 'USPP + noncolin is not implemented', 1 )
- !
- ! If this is the first entry, we calculate and save the coefficients B from Eq.(15)
- ! B. Walker and R. Gebauer, J. Chem. Phys. 127, 164106 (2007).
- ! If this is not the first entry, we do not recalculate the coefficients B
- ! (they are already allocated) but use the ones which were already calculated
- ! and saved (if recalc=.false.).
- !
- IF (.NOT.ALLOCATED(bbnc)) recalc = .true.
- IF (recalc .AND. ALLOCATED(bbnc)) DEALLOCATE(bbnc)
- !
- ! If recalc=.true. we (re)calculate the coefficients lambda_nm defined by Eq.(16).
- !
- IF ( recalc ) THEN
- !
- ALLOCATE(bbnc(nkb,nkb,nspin_mag,nksq))
- bbnc = (0.0d0,0.0d0)
- !
- ALLOCATE(ps(nkb,nkb,nspin_mag))
- !
- DO ik1 = 1, nksq
- !
- ikk = ikks(ik1)
- ikq = ikqs(ik1)
- npwq = ngk(ikq)
- !
- ! Calculate beta-functions vkb for a given k+q point.
- !
- CALL init_us_2 (npwq, igk_k(1,ikq), xk(1,ikq), vkb)
- !
- ! Calculate the coefficients B_ij defined by Eq.(15).
- ! B_ij = , where beta(i) = vkb(i).
- !
- CALL ZGEMM('C','N',nkb,nkb,npwq,(1.d0,0.d0),vkb, &
- & lda,vkb,lda,(0.d0,0.d0),bbnc(1,1,1,ik1),nkb)
- !
- IF (lspinorb) THEN
- bbnc(:,:,2,ik1) = bbnc(:,:,1,ik1)
- bbnc(:,:,3,ik1) = bbnc(:,:,1,ik1)
- bbnc(:,:,4,ik1) = bbnc(:,:,1,ik1)
- ENDIF
- !
-#if defined(__MPI)
- CALL mp_sum(bbnc(:,:,:,ik1), intra_bgrp_comm)
-#endif
- !
- ps(:,:,:) = (0.d0,0.d0)
- !
- ! Calculate the product of q_nm and B_ij of Eq.(16).
- !
- ijkb0 = 0
- do nt=1,ntyp
- if (upf(nt)%tvanp) then
- do na=1,nat
- if(ityp(na).eq.nt) then
- do ii=1,nkb
- do jh=1,nh(nt)
- !
- jkb=ijkb0 + jh
- !
- do ih=1,nh(nt)
- !
- ikb = ijkb0 + ih
- !
- if (lspinorb) then
- do ipol=1,4
- ps(ikb,ii,ipol) = ps(ikb,ii,ipol) + &
- & bbnc(jkb,ii,ipol,ik1) * qq_so(ih,jh,ipol,nt)
- enddo
- else
- ps(ikb,ii,1) = ps(ikb,ii,1) + &
- & bbnc(jkb,ii,1,ik1) * qq_nt(ih,jh,nt)
- endif
- !
- enddo
- enddo
- enddo
- ijkb0 = ijkb0+nh(nt)
- endif
- enddo
- else
- DO na = 1, nat
- IF ( ityp(na) == nt ) ijkb0 = ijkb0 + nh(nt)
- END DO
- endif
- enddo
- !
- ! Add an identity to q_nm * B_ij [see Eq.(16)].
- ! ps = (1 + q*B)
- !
- DO ii=1,nkb
- IF (lspinorb) THEN
- DO ipol=1,4
- ps(ii,ii,ipol) = ps(ii,ii,ipol) + (1.d0,0.d0)
- ENDDO
- ELSE
- ps(ii,ii,1) = ps(ii,ii,1) + (1.d0,0.d0)
- ENDIF
- ENDDO
- !
- ! Invert matrix: (1 + q*B)^{-1}
- ! WARNING: How to do the invertion of the matrix in the spin-orbit case,
- ! when there are 4 components?
- !
- IF (lspinorb) THEN
- DO ipol=1,4
- CALL invmat( nkb, ps(:,:,ipol) )
- ENDDO
- ELSE
- CALL invmat ( nkb, ps(:,:,1) )
- ENDIF
- !
- ! Finally, let us calculate lambda_nm = -(1+q*B)^{-1} * q
- !
- ! Use the array bbnc as a workspace and put there
- ! the values of the lambda coefficients.
- !
- bbnc(:,:,:,ik1) = (0.d0,0.d0)
- !
- ijkb0 = 0
- do nt=1,ntyp
- if (upf(nt)%tvanp) then
- do na=1,nat
- if(ityp(na).eq.nt) then
- do ii=1,nkb
- do jh=1,nh(nt)
- jkb=ijkb0 + jh
- do ih=1,nh(nt)
- ikb = ijkb0 + ih
- if (lspinorb) then
- bbnc(ii,jkb,1,ik1) = bbnc(ii,jkb,1,ik1) &
- & - ps(ii,ikb,1) * qq_so(ih,jh,1,nt) &
- & - ps(ii,ikb,2) * qq_so(ih,jh,3,nt)
- !
- bbnc(ii,jkb,2,ik1) = bbnc(ii,jkb,2,ik1) &
- & - ps(ii,ikb,1) * qq_so(ih,jh,2,nt) &
- & - ps(ii,ikb,2) * qq_so(ih,jh,4,nt)
- !
- bbnc(ii,jkb,3,ik1) = bbnc(ii,jkb,3,ik1) &
- & - ps(ii,ikb,3) * qq_so(ih,jh,1,nt) &
- & - ps(ii,ikb,4) * qq_so(ih,jh,3,nt)
- !
- bbnc(ii,jkb,4,ik1) = bbnc(ii,jkb,4,ik1) &
- & - ps(ii,ikb,3) * qq_so(ih,jh,2,nt) &
- & - ps(ii,ikb,4) * qq_so(ih,jh,4,nt)
- else
- bbnc(ii,jkb,1,ik1) = bbnc(ii,jkb,1,ik1) &
- & - ps(ii,ikb,1)*qq_nt(ih,jh,nt)
- endif
- !
- enddo
- enddo
- enddo
- ijkb0 = ijkb0+nh(nt)
- endif
- enddo
- else
- DO na = 1, nat
- IF ( ityp(na) == nt ) ijkb0 = ijkb0 + nh(nt)
- END DO
- endif
- enddo
- !
- enddo ! loop on k points
- !
- deallocate(ps)
- !
- endif
- !
! Now set up the indices ikk and ikq such that they
! correspond to the points k and k+q using the index ik,
! which was passed to this routine as an input.
@@ -891,7 +226,7 @@ SUBROUTINE sm1_psiq_nc()
ikk = ikks(ik)
ikq = ikqs(ik)
!
- IF (n.NE.ngk(ikq)) CALL errore( 'sm1_psiq_nc', &
+ IF (n/=ngk(ikq)) CALL errore( 'sm1_psiq_nc', &
& 'Mismatch in the number of plane waves', 1 )
!
! Calculate beta-functions vkb for a given k+q point.
@@ -906,37 +241,17 @@ SUBROUTINE sm1_psiq_nc()
!
! Use ps as a work space.
!
- ALLOCATE(ps(nkb,npol,m))
- ps(:,:,:) = (0.d0,0.d0)
+ ALLOCATE(ps(nkb*npol,m))
+ ps(:,:) = (0.d0,0.d0)
!
! Apply the operator S^{-1}, given by Eq.(13), to the functions psi.
- ! Let's do this in 2 steps.
+ ! Let's DO this in 2 steps.
!
! Step 1 : calculate the product
! ps = lambda *
!
- DO ibnd=1,m
- DO jkb=1,nkb
- DO ii=1,nkb
- IF (lspinorb) THEN
- !
- ps(jkb,1,ibnd) = ps(jkb,1,ibnd) &
- & + bbnc(jkb,ii,1,ik) * becp%nc(ii,1,ibnd) &
- & + bbnc(jkb,ii,2,ik) * becp%nc(ii,2,ibnd)
- !
- ps(jkb,2,ibnd) = ps(jkb,2,ibnd) &
- & + bbnc(jkb,ii,3,ik) * becp%nc(ii,1,ibnd) &
- & + bbnc(jkb,ii,4,ik) * becp%nc(ii,2,ibnd)
- !
- ELSE
- DO ipol=1,npol
- ps(jkb,ipol,ibnd) = ps(jkb,ipol,ibnd) &
- & + bbnc(jkb,ii,1,ik) * becp%nc(ii,ipol,ibnd)
- ENDDO
- ENDIF
- ENDDO
- ENDDO
- ENDDO
+ CALL ZGEMM( 'N', 'N', nkb*npol, m, nkb*npol, (1.D0, 0.D0), bbnc_sm1(1,1,ik), &
+ nkb*npol, becp%nc, nkb*npol, (1.D0, 0.D0), ps, nkb*npol )
!
! Step 2 : |spsi> = S^{-1} * |psi> = |psi> + ps * |beta>
!
@@ -947,6 +262,218 @@ SUBROUTINE sm1_psiq_nc()
!
RETURN
!
-END SUBROUTINE sm1_psiq_nc
+END SUBROUTINE sm1_psi_nc
-END SUBROUTINE lr_sm1_psiq
+END SUBROUTINE lr_sm1_psi_tpw
+
+SUBROUTINE lr_sm1_initialize()
+!
+! This routine initializes the coefficients of S^-1 and saves them in
+! bbg, bbk or bbnc
+!
+USE kinds, ONLY : DP
+USE control_flags, ONLY : gamma_only
+USE lrus, ONLY : bbg,bbk, bbnc_sm1
+USE klist, ONLY : xk, ngk, igk_k
+USE qpoint, ONLY : nksq, ikks, ikqs
+USE wvfct, ONLY : npwx
+USE becmod, ONLY : calbec
+USE uspp, ONLY : vkb, nkb, qq_nt, qq_so
+USE uspp_param, ONLY : nh, upf
+USE ions_base, ONLY : ityp,nat,ntyp=>nsp
+USE mp, ONLY : mp_sum
+USE mp_global, ONLY : intra_bgrp_comm
+USE noncollin_module, ONLY : noncolin, npol
+USE matrix_inversion, ONLY : invmat
+
+IMPLICIT NONE
+!
+INTEGER :: ik1, ikk, ikq, npw, npwq, na, nt, ikb, jkb, ijkb0, ih, jh, ii, &
+ ipol, jpol, kpol, ijs, iis, ikbs, jkbs, iks, kjs
+REAL(DP), ALLOCATABLE :: psr(:,:)
+COMPLEX(DP), ALLOCATABLE :: ps(:,:), bbnc_nc(:,:)
+
+IF (gamma_only) THEN
+ bbg = 0.0d0
+ ALLOCATE(psr(nkb,nkb))
+ psr(:,:) = (0.d0,0.d0)
+ELSE
+ IF (noncolin) THEN
+ ALLOCATE(bbnc_nc(nkb,nkb))
+ bbnc_nc = (0.0d0,0.0d0)
+ ELSE
+ bbk = (0.0d0,0.0d0)
+ ENDIF
+ ALLOCATE(ps(nkb*npol,nkb*npol))
+ ps(:,:) = (0.d0,0.d0)
+ENDIF
+
+DO ik1 = 1, nksq
+ !
+ ikk = ikks(ik1)
+ ikq = ikqs(ik1)
+ npw = ngk(ikk)
+ npwq = ngk(ikq)
+ !
+ ! Calculate beta-functions vkb for a given k+q point.
+ !
+ CALL init_us_2 (npwq, igk_k(1,ikq), xk(1,ikq), vkb)
+ !
+ ! Calculate the coefficients B_ij defined by Eq.(15).
+ ! B_ij = , where beta(i) = vkb(i).
+ !
+ IF (gamma_only) THEN
+ CALL calbec (npw, vkb, vkb, bbg, nkb)
+ ELSEIF (noncolin) THEN
+ CALL zgemm('C','N',nkb,nkb,npwq,(1.d0,0.d0),vkb, &
+ npwx,vkb,npwx,(0.d0,0.d0),bbnc_nc,nkb)
+ CALL mp_sum(bbnc_nc, intra_bgrp_comm)
+ ELSE
+ CALL zgemm('C','N',nkb,nkb,npwq,(1.d0,0.d0),vkb, &
+ npwx,vkb,npwx,(0.d0,0.d0),bbk(1,1,ik1),nkb)
+ CALL mp_sum(bbk(:,:,ik1), intra_bgrp_comm)
+ ENDIF
+ !
+ !
+ ! Calculate the product of q_nm and B_ij of Eq.(16).
+ !
+ ijkb0 = 0
+ DO nt=1,ntyp
+ IF (upf(nt)%tvanp) THEN
+ DO na=1,nat
+ IF (ityp(na)==nt) THEN
+ DO ii=1,nkb
+ DO jh=1,nh(nt)
+ !
+ jkb=ijkb0 + jh
+ !
+ DO ih=1,nh(nt)
+ !
+ ikb = ijkb0 + ih
+ !
+ IF (gamma_only) THEN
+ psr(ikb,ii) = psr(ikb,ii) + bbg(jkb,ii) &
+ * qq_nt(ih,jh,nt)
+ ELSEIF(noncolin) THEN
+ ijs=0
+ DO ipol=1, npol
+ ikbs = ikb + nkb * ( ipol - 1 )
+ DO jpol=1, npol
+ iis = ii + nkb * ( jpol - 1 )
+ ijs = ijs + 1
+ ps(ikbs,iis) = ps(ikbs,iis) + &
+ bbnc_nc(jkb,ii)*qq_so(ih,jh,ijs,nt)
+ ENDDO
+ ENDDO
+ ELSE
+ ps(ikb,ii) = ps(ikb,ii) + bbk(jkb,ii,ik1) &
+ * qq_nt(ih,jh,nt)
+ ENDIF
+ !
+ ENDDO
+ ENDDO
+ ENDDO
+ ijkb0 = ijkb0+nh(nt)
+ ENDIF
+ ENDDO
+ ELSE
+ DO na = 1, nat
+ IF ( ityp(na) == nt ) ijkb0 = ijkb0 + nh(nt)
+ ENDDO
+ ENDIF
+ ENDDO
+ !
+ ! Add an identity to q_nm * B_ij [see Eq.(16)].
+ ! ps = (1 + q*B)
+ !
+ IF (gamma_only) THEN
+ DO ii=1,nkb
+ psr(ii,ii) = psr(ii,ii) + 1.d0
+ ENDDO
+ ELSE
+ DO ii=1,nkb*npol
+ ps(ii,ii) = ps(ii,ii) + (1.d0,0.d0)
+ ENDDO
+ ENDIF
+ !
+ ! Invert matrix: (1 + q*B)^{-1}
+ !
+ IF (gamma_only) THEN
+ CALL invmat( nkb, psr )
+ ELSE
+ CALL invmat( nkb*npol, ps )
+ ENDIF
+ !
+ ! Use the array bbk as a work space in order to save the memory.
+ !
+ IF (gamma_only) THEN
+ bbg(:,:)=0.0_DP
+ ELSEIF (noncolin) THEN
+ bbnc_sm1(:,:,ik1) = (0.d0,0.d0)
+ ELSE
+ bbk(:,:,ik1) = (0.d0,0.d0)
+ ENDIF
+ !
+ ! Finally, let us calculate lambda_nm = -(1+q*B)^{-1} * q
+ ! Let us use the array bbk and put there the values of the lambda
+ ! coefficients.
+ !
+ ijkb0 = 0
+ DO nt=1,ntyp
+ IF (upf(nt)%tvanp) THEN
+ DO na=1,nat
+ IF (ityp(na)==nt) THEN
+ DO ii=1,nkb*npol
+ DO jh=1,nh(nt)
+ !
+ jkb=ijkb0 + jh
+ !
+ DO ih=1,nh(nt)
+ !
+ ikb = ijkb0 + ih
+ !
+ IF (gamma_only) THEN
+ bbg(ii,jkb) = bbg(ii,jkb) &
+ - psr(ii,ikb) * qq_nt(ih,jh,nt)
+
+ ELSEIF (noncolin) THEN
+ kjs = 0
+ DO kpol=1,npol
+ ikbs = ikb + nkb * (kpol-1)
+ DO jpol=1,npol
+ jkbs=jkb + nkb * (jpol-1)
+ kjs=kjs+1
+ bbnc_sm1(ii,jkbs,ik1) = &
+ bbnc_sm1(ii,jkbs,ik1) - &
+ ps(ii,ikbs)*qq_so(ih,jh,kjs,nt)
+ ENDDO
+ ENDDO
+ ELSE
+ bbk(ii,jkb,ik1) = bbk(ii,jkb,ik1) - &
+ ps(ii,ikb) * qq_nt(ih,jh,nt)
+ ENDIF
+ !
+ ENDDO
+ ENDDO
+ ENDDO
+ ijkb0 = ijkb0+nh(nt)
+ ENDIF
+ ENDDO
+ ELSE
+ DO na = 1, nat
+ IF ( ityp(na) == nt ) ijkb0 = ijkb0 + nh(nt)
+ ENDDO
+ ENDIF
+ ENDDO
+ !
+ENDDO ! loop on k points
+!
+IF (gamma_only) THEN
+ DEALLOCATE(psr)
+ELSE
+ IF (noncolin) DEALLOCATE(bbnc_nc)
+ DEALLOCATE(ps)
+ENDIF
+
+RETURN
+END SUBROUTINE lr_sm1_initialize
diff --git a/LR_Modules/lrcom.f90 b/LR_Modules/lrcom.f90
index 816a12c22..2be1c6499 100644
--- a/LR_Modules/lrcom.f90
+++ b/LR_Modules/lrcom.f90
@@ -154,6 +154,7 @@ MODULE lrus
COMPLEX (DP), ALLOCATABLE :: bbk(:,:,:) ! nkb, nkb, nks)
! for k points
COMPLEX (DP), ALLOCATABLE :: bbnc(:,:,:,:) ! nkb, nkb, nspin_mag, nks)
+ COMPLEX (DP), ALLOCATABLE :: bbnc_sm1(:,:,:) ! nkb, nkb, nspin_mag, nks)
! for the noncollinear case
! bbg = < beta^N_i | beta^P_j >
! bbg/bbk/bbnc are the scalar products of beta functions
diff --git a/TDDFPT/src/Makefile b/TDDFPT/src/Makefile
index 828103a5c..9c7f4df7b 100644
--- a/TDDFPT/src/Makefile
+++ b/TDDFPT/src/Makefile
@@ -29,6 +29,10 @@ lr_normalise.o \
lr_lanczos.o \
lr_apply_liouvillian.o \
lr_dv_setup.o \
+compute_intq.o \
+set_intq_nc.o \
+transform_intq_nc.o \
+transform_intq_so.o \
lr_solve_e.o \
lr_dvpsi_e.o \
stop_lr.o \
diff --git a/TDDFPT/src/compute_intq.f90 b/TDDFPT/src/compute_intq.f90
new file mode 100644
index 000000000..efee9937e
--- /dev/null
+++ b/TDDFPT/src/compute_intq.f90
@@ -0,0 +1,89 @@
+!
+! Copyright (C) 2001-2008 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,
+! or http://www.gnu.org/copyleft/gpl.txt .
+!
+!
+!----------------------------------------------------------------------
+SUBROUTINE compute_intq
+ !----------------------------------------------------------------------
+ !
+ ! This routine computes the contribution of the selfconsistent
+ ! change of the potential to the known part of the linear
+ ! system and adds it to dvpsi.
+ !
+ USE kinds, ONLY : DP
+ USE ions_base, ONLY : nat, ityp, ntyp => nsp
+ USE noncollin_module, ONLY : noncolin
+ USE cell_base, ONLY : omega
+ USE uspp, ONLY : okvan
+ USE uspp_param, ONLY : upf, lmaxq, nh, nhm
+
+ USE lr_variables, ONLY : intq
+ USE qpoint, ONLY : xq, eigqts
+
+ IMPLICIT NONE
+
+ INTEGER :: na, ig, nt, ir, ih, jh
+ ! countera
+
+ REAL(DP), ALLOCATABLE :: ylmk0 (:,:)
+ ! the modulus of q+G
+ ! the values of q+G
+ ! the spherical harmonics
+
+
+ ! work space
+ COMPLEX(DP) :: qgm(1), aux1
+ REAL(DP) :: qmod(1), zero(3,1), qg(3,1)
+
+ IF (.NOT.okvan) RETURN
+ CALL start_clock ('compute_intq')
+
+ intq (:,:,:) = (0.D0, 0.0D0)
+ ALLOCATE (ylmk0(1 , lmaxq * lmaxq))
+ !
+ ! first compute the spherical harmonics
+ !
+ zero=0.0_DP
+ CALL setqmod (1, xq, zero, qmod, qg)
+ CALL ylmr2 (lmaxq * lmaxq, 1, qg, qmod, ylmk0)
+ qmod(1) = SQRT (qmod(1) )
+
+ DO nt = 1, ntyp
+ IF (upf(nt)%tvanp ) THEN
+ DO ih = 1, nh (nt)
+ DO jh = ih, nh (nt)
+ CALL qvan2 (1, ih, jh, nt, qmod, qgm, ylmk0)
+ DO na = 1, nat
+ IF (ityp (na) == nt) THEN
+ aux1 = qgm(1) * eigqts(na)
+ intq(ih,jh,na) = omega * CONJG(aux1)
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
+ DO na = 1, nat
+ IF (ityp(na) == nt) THEN
+ !
+ ! We use the symmetry properties of the ps factor
+ !
+ DO ih = 1, nh (nt)
+ DO jh = ih, nh (nt)
+ intq(jh,ih,na) = intq(ih,jh,na)
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDDO
+ ENDIF
+ ENDDO
+
+ IF (noncolin) CALL set_intq_nc()
+
+ DEALLOCATE (ylmk0)
+
+ CALL stop_clock ('compute_intq')
+ RETURN
+END SUBROUTINE compute_intq
diff --git a/TDDFPT/src/dveqpsi_us_only.f90 b/TDDFPT/src/dveqpsi_us_only.f90
new file mode 100644
index 000000000..5704a2cea
--- /dev/null
+++ b/TDDFPT/src/dveqpsi_us_only.f90
@@ -0,0 +1,110 @@
+!
+! Copyright (C) 2001-2019 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,
+! or http://www.gnu.org/copyleft/gpl.txt .
+!
+!----------------------------------------------------------------------
+SUBROUTINE dveqpsi_us_only (npwq, ik)
+ !----------------------------------------------------------------------
+ !
+ ! This routine computes the contribution of the Fourier transform
+ ! of the augmentation function at the given q, and adds it to
+ ! dvpsi.
+ !
+ USE kinds, ONLY : DP
+ USE uspp_param, ONLY : upf, nh
+ USE uspp, ONLY : vkb, okvan
+! modules from pwcom
+ USE lsda_mod, ONLY : lsda, current_spin, isk
+ USE ions_base, ONLY : ntyp => nsp, nat, ityp
+ USE wvfct, ONLY : nbnd, npwx
+ USE noncollin_module, ONLY : noncolin, npol
+! modules from phcom
+ USE qpoint, ONLY : ikks
+ USE optical, ONLY : intq, intq_nc
+ USE lrus, ONLY : becp1
+ USE eqv, ONLY : dvpsi
+ IMPLICIT NONE
+ !
+ ! The dummy variables
+ !
+ INTEGER :: ik, npwq
+ ! input: the k point
+ ! input: number of plane waves at the q point
+
+ ! And the local variables
+ !
+ INTEGER :: na, nt, ibnd, ih, jh, ijkb0, ikk, ikb, jkb, is, js, ijs
+ ! counter on atoms
+ ! counter on atomic types
+ ! counter on bands
+ ! counter on beta functions
+ ! counter on beta functions
+ ! auxiliary variable for indexing
+ ! counter on the k points
+ ! counter on vkb
+ ! counter on vkb
+ COMPLEX(DP) :: sum0, sum_nc(npol)
+ ! auxiliary variable
+
+ IF (.NOT.okvan) RETURN
+ CALL start_clock ('dveqpsi_us_only')
+ ikk = ikks(ik)
+ IF (lsda) current_spin = isk (ikk)
+ ijkb0 = 0
+ DO nt = 1, ntyp
+ IF (upf(nt)%tvanp ) THEN
+ DO na = 1, nat
+ IF (ityp (na)==nt) THEN
+ !
+ ! we multiply the integral for the becp term and the beta_n
+ !
+ DO ibnd = 1, nbnd
+ DO ih = 1, nh (nt)
+ ikb = ijkb0 + ih
+ IF (noncolin) THEN
+ sum_nc = (0.d0, 0.d0)
+ ELSE
+ sum0 = (0.d0, 0.d0)
+ END IF
+ DO jh = 1, nh (nt)
+ jkb = ijkb0 + jh
+ IF (noncolin) THEN
+ ijs=0
+ DO is=1,npol
+ DO js=1,npol
+ ijs=ijs+1
+ sum_nc(is)=sum_nc(is)+ &
+ intq_nc(ih,jh,na,ijs)* &
+ becp1(ik)%nc(jkb, js, ibnd)
+ ENDDO
+ ENDDO
+ ELSE
+ sum0 = sum0 + intq (ih, jh, na)* &
+ becp1(ik)%k(jkb, ibnd)
+ ENDIF
+ ENDDO
+ IF (noncolin) THEN
+ CALL zaxpy(npwq,sum_nc(1),vkb(1,ikb),1,dvpsi(1,ibnd),1)
+ CALL zaxpy(npwq,sum_nc(2),vkb(1,ikb),1, &
+ dvpsi(1+npwx,ibnd),1)
+ ELSE
+ CALL zaxpy(npwq,sum0,vkb(1,ikb),1,dvpsi(1,ibnd),1)
+ ENDIF
+ ENDDO
+ ENDDO
+ ijkb0 = ijkb0 + nh (nt)
+ ENDIF
+ ENDDO
+ ELSE
+ DO na = 1, nat
+ IF (ityp (na)==nt) ijkb0 = ijkb0 + nh (nt)
+ ENDDO
+ ENDIF
+ ENDDO
+
+ CALL stop_clock ('dveqpsi_us_only')
+ RETURN
+END SUBROUTINE dveqpsi_us_only
diff --git a/TDDFPT/src/lr_addus_dvpsi.f90 b/TDDFPT/src/lr_addus_dvpsi.f90
index de40222ad..f2aaee9d8 100644
--- a/TDDFPT/src/lr_addus_dvpsi.f90
+++ b/TDDFPT/src/lr_addus_dvpsi.f90
@@ -1,12 +1,13 @@
!
-! Copyright (C) 2001-2016 Quantum ESPRESSO group
+! Copyright (C) 2001-2019 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,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!-----------------------------------------------------------------------------
-SUBROUTINE lr_addus_dvpsi ( ik, lda, n, m, psi, dpsi )
+SUBROUTINE lr_addus_dvpsi (npwq, ik,psi,dvpsi)
+ !----------------------------------------------------------------------
!---------------------------------------------------------------------------
!
! ... Calculate the ultrasoft term of the perturbation exp(iq*r),
@@ -20,118 +21,102 @@ SUBROUTINE lr_addus_dvpsi ( ik, lda, n, m, psi, dpsi )
! ... m number of bands of psi
!
! Written by Iurii Timrov (2015)
+ ! Generalized to the relativistic case by Andrea Dal Corso (2018)
!
- USE kinds, ONLY : DP
- USE ions_base, ONLY : nat, ityp, ntyp => nsp
- USE uspp, ONLY : okvan
- USE uspp_param, ONLY : upf, lmaxq, nh
- USE paw_variables, ONLY : okpaw
- USE noncollin_module, ONLY : npol, noncolin
- USE uspp, ONLY : vkb, nkb, indv_ijkb0
- USE cell_base, ONLY : omega
- USE lrus, ONLY : becp1
- USE qpoint, ONLY : xq, eigqts
- !
+ USE kinds, ONLY : DP
+ USE uspp_param, ONLY : upf, nh
+ USE uspp, ONLY : vkb, okvan
+ USE lsda_mod, ONLY : lsda, current_spin, isk
+ USE ions_base, ONLY : ntyp => nsp, nat, ityp
+ USE wvfct, ONLY : nbnd, npwx
+ USE noncollin_module, ONLY : noncolin, npol
+ USE qpoint, ONLY : ikks
+ USE lr_variables, ONLY : intq, intq_nc
+ USE lrus, ONLY : becp1
IMPLICIT NONE
!
- INTEGER, INTENT(in) :: ik, lda, n, m
- COMPLEX(DP), INTENT(in) :: psi(lda*npol,m)
+ ! The dummy variables
+ !
+ COMPLEX(DP), INTENT(in) :: psi(npwx*npol,nbnd)
! input: wavefunction u_n,k
- COMPLEX(DP), INTENT(out) :: dpsi(lda*npol,m)
+ COMPLEX(DP), INTENT(out) :: dvpsi(npwx*npol,nbnd)
! output: sum of the input wavefunction and the USPP term
+ INTEGER :: ik, npwq
+ ! input: the k point
+ ! input: number of plane waves at the q point
+
+ ! And the local variables
!
- ! the local variables
- !
- INTEGER :: na, nt, ih, jh
+ INTEGER :: na, nt, ibnd, ih, jh, ijkb0, ikk, ikb, jkb, is, js, ijs
! counter on atoms
- ! counter on atomic type
+ ! counter on atomic types
+ ! counter on bands
! counter on beta functions
! counter on beta functions
- !
- REAL(DP) :: qmod ! the modulus of q
- REAL(DP), ALLOCATABLE :: ylmk0(:) ! the spherical harmonics
- COMPLEX(DP) :: qgm(1) ! FFT of Q(r)
- COMPLEX(DP), ALLOCATABLE :: ps(:,:), qqc(:,:),qqc_d(:,:) ! auxiliary arrays
- !
- dpsi = psi
- !
+ ! auxiliary variable for indexing
+ ! counter on the k points
+ ! counter on vkb
+ ! counter on vkb
+ COMPLEX(DP) :: sum0, sum_nc(npol)
+ ! auxiliary variable
+
IF (.NOT.okvan) RETURN
- !
CALL start_clock ('lr_addus_dvpsi')
- !
- IF (noncolin) CALL errore( 'lr_addus_dvpsi', 'Noncollinear case is not supported', 1 )
- IF (okpaw) CALL errore( 'lr_addus_dvpsi', 'PAW is not supported', 1 )
- !
- ALLOCATE (ylmk0(lmaxq*lmaxq))
- ALLOCATE (ps(nkb,m))
- ps(:,:) = (0.d0, 0.d0)
- !
- qmod = xq(1)**2 + xq(2)**2 + xq(3)**2
- !
- ! Calculate sphrecial harmonics
- !
- CALL ylmr2 (lmaxq*lmaxq, 1, xq, qmod, ylmk0)
- !
- qmod = sqrt(qmod)
- !
+ dvpsi=psi
+ ikk = ikks(ik)
+ IF (lsda) current_spin = isk (ikk)
+ ijkb0 = 0
DO nt = 1, ntyp
- IF (upf(nt)%tvanp) THEN
- !
- ! Calculate the Fourier transform of the Q functions.
- !
- ALLOCATE (qqc(nh(nt),nh(nt)))
- ALLOCATE (qqc_d(nh(nt),nh(nt)))
- qqc(:,:) = (0.d0, 0.d0)
- !
- DO ih = 1, nh(nt)
- DO jh = ih, nh(nt)
- !
- CALL qvan2 (1, ih, jh, nt, qmod, qgm(1), ylmk0)
- qqc(ih,jh) = omega * qgm(1)
- qqc(jh,ih) = qqc(ih,jh)
- !
+ IF (upf(nt)%tvanp ) THEN
+ DO na = 1, nat
+ IF (ityp (na)==nt) THEN
+ !
+ ! we multiply the integral for the becp term and the beta_n
+ !
+ DO ibnd = 1, nbnd
+ DO ih = 1, nh (nt)
+ ikb = ijkb0 + ih
+ IF (noncolin) THEN
+ sum_nc = (0.d0, 0.d0)
+ ELSE
+ sum0 = (0.d0, 0.d0)
+ END IF
+ DO jh = 1, nh (nt)
+ jkb = ijkb0 + jh
+ IF (noncolin) THEN
+ ijs=0
+ DO is=1,npol
+ DO js=1,npol
+ ijs=ijs+1
+ sum_nc(is)=sum_nc(is)+ &
+ intq_nc(ih,jh,na,ijs)* &
+ becp1(ik)%nc(jkb, js, ibnd)
+ ENDDO
+ ENDDO
+ ELSE
+ sum0 = sum0 + intq (ih, jh, na)* &
+ becp1(ik)%k(jkb, ibnd)
+ ENDIF
+ ENDDO
+ IF (noncolin) THEN
+ CALL zaxpy(npwq,sum_nc(1),vkb(1,ikb),1,dvpsi(1,ibnd),1)
+ CALL zaxpy(npwq,sum_nc(2),vkb(1,ikb),1, &
+ dvpsi(1+npwx,ibnd),1)
+ ELSE
+ CALL zaxpy(npwq,sum0,vkb(1,ikb),1,dvpsi(1,ibnd),1)
+ ENDIF
+ ENDDO
+ ENDDO
+ ijkb0 = ijkb0 + nh (nt)
+ ENDIF
ENDDO
- ENDDO
- !
- ! Calculate a product of Q^*(q) with
- !
- DO na = 1, nat
- IF (ityp (na).eq.nt) THEN
- !
- qqc_d = qqc * eigqts(na)
- !
- CALL ZGEMM('C','N', nh(nt), m, nh(nt), (1.0_dp,0.0_dp), &
- & qqc_d, nh(nt), becp1(ik)%k(indv_ijkb0(na)+1,1), nkb, &
- & (0.0_dp,0.0_dp), ps(indv_ijkb0(na)+1,1), nkb )
- !
- ENDIF
- ENDDO
- !
- DEALLOCATE (qqc)
- DEALLOCATE (qqc_d)
- !
- ENDIF
+ ELSE
+ DO na = 1, nat
+ IF (ityp (na)==nt) ijkb0 = ijkb0 + nh (nt)
+ ENDDO
+ ENDIF
ENDDO
- !
- ! Sum up the normal and ultrasoft terms.
- !
- IF ( m == 1 ) THEN
- !
- CALL ZGEMV( 'N', n, nkb, ( 1.D0, 0.D0 ), vkb, lda, &
- & ps, 1, ( 1.D0, 0.D0 ), dpsi, 1 )
- !
- ELSE
- !
- CALL ZGEMM( 'N', 'N', n, m, nkb, ( 1.D0, 0.D0 ), vkb, lda, &
- & ps, nkb, ( 1.D0, 0.D0 ), dpsi, lda )
- !
- ENDIF
- !
- DEALLOCATE (ylmk0)
- DEALLOCATE (ps)
- !
+
CALL stop_clock ('lr_addus_dvpsi')
- !
RETURN
- !
END SUBROUTINE lr_addus_dvpsi
diff --git a/TDDFPT/src/lr_alloc_init.f90 b/TDDFPT/src/lr_alloc_init.f90
index 67795b195..b8359ee1e 100644
--- a/TDDFPT/src/lr_alloc_init.f90
+++ b/TDDFPT/src/lr_alloc_init.f90
@@ -15,6 +15,7 @@ SUBROUTINE lr_alloc_init()
USE kinds, ONLY : dp
USE ions_base, ONLY : nat
USE uspp, ONLY : nkb, okvan
+ USE lrus, ONLY : bbk, bbnc_sm1
USE uspp_param, ONLY : nhm
USE fft_base, ONLY : dfftp, dffts
USE klist, ONLY : nks
@@ -153,12 +154,13 @@ SUBROUTINE lr_alloc_init()
ENDIF
!
! Optical case: allocate the response charge-density
- !
+ !
IF (.NOT.eels) THEN
!
IF (gamma_only) THEN
ALLOCATE(rho_1(dfftp%nnr,nspin_mag))
rho_1(:,:)=0.0d0
+ ALLOCATE(bbg(nkb,nkb))
ELSE
ALLOCATE(rho_1c(dfftp%nnr,nspin_mag))
rho_1c(:,:)=(0.0d0,0.0d0)
@@ -181,6 +183,9 @@ SUBROUTINE lr_alloc_init()
IF (noncolin) THEN
ALLOCATE (int3_nc(nhm,nhm,nat,nspin,1))
int3_nc = (0.0d0, 0.0d0)
+ ALLOCATE(bbnc_sm1(nkb*npol, nkb*npol, nksq))
+ ELSE
+ ALLOCATE(bbk(nkb, nkb, nksq))
ENDIF
!
ENDIF
diff --git a/TDDFPT/src/lr_apply_liouvillian.f90 b/TDDFPT/src/lr_apply_liouvillian.f90
index 8f6433b8e..5417a9783 100644
--- a/TDDFPT/src/lr_apply_liouvillian.f90
+++ b/TDDFPT/src/lr_apply_liouvillian.f90
@@ -296,8 +296,10 @@ SUBROUTINE lr_apply_liouvillian( evc1, evc1_new, interaction )
!
DO ik = 1, nks
!
- CALL lr_sm1_psi (.FALSE., ik, npwx, ngk(ik), nbnd, &
+ CALL lr_sm1_psi_tpw ( ik, npwx, ngk(ik), nbnd, &
& sevc1_new(1,1,ik), evc1_new(1,1,ik))
+! CALL lr_sm1_psi (.FALSE., ik, npwx, ngk(ik), nbnd, &
+! & sevc1_new(1,1,ik), evc1_new(1,1,ik))
!
ENDDO
!
diff --git a/TDDFPT/src/lr_apply_liouvillian_eels.f90 b/TDDFPT/src/lr_apply_liouvillian_eels.f90
index 7281a3fda..dd61ae5a0 100644
--- a/TDDFPT/src/lr_apply_liouvillian_eels.f90
+++ b/TDDFPT/src/lr_apply_liouvillian_eels.f90
@@ -291,7 +291,7 @@ SUBROUTINE lr_apply_liouvillian_eels ( evc1, evc1_new, interaction )
! evc1_new = S^{-1} * sevc1_new
! If not ultrasoft: evc1_new = sevc1_new
!
- CALL lr_sm1_psiq (.FALSE., ik, npwx, npwq, nbnd_occ(ikk), &
+ CALL lr_sm1_psi_tpw ( ik, npwx, npwq, nbnd_occ(ikk), &
& sevc1_new(1,1,ik), evc1_new(1,1,ik))
!
ENDDO ! loop on ik
diff --git a/TDDFPT/src/lr_dealloc.f90 b/TDDFPT/src/lr_dealloc.f90
index d478cc57e..dd6172284 100644
--- a/TDDFPT/src/lr_dealloc.f90
+++ b/TDDFPT/src/lr_dealloc.f90
@@ -25,7 +25,7 @@ SUBROUTINE lr_dealloc()
USE lr_exx_kernel, ONLY : lr_exx_dealloc
USE becmod, ONLY : bec_type, becp, deallocate_bec_type
USE lrus, ONLY : int3, int3_nc, becp1, &
- & bbg, bbk, bbnc
+ & bbg, bbk, bbnc_sm1, bbnc
USE qpoint, ONLY : ikks, ikqs, eigqts
USE eqv, ONLY : dmuxc, evq, dpsi, dvpsi
USE control_lr, ONLY : nbnd_occ
@@ -53,6 +53,7 @@ SUBROUTINE lr_dealloc()
IF (allocated(bbg)) DEALLOCATE(bbg)
IF (allocated(bbk)) DEALLOCATE(bbk)
IF (allocated(bbnc)) DEALLOCATE(bbnc)
+ IF (allocated(bbnc_sm1)) DEALLOCATE(bbnc_sm1)
!
IF (project) THEN
DEALLOCATE(F)
diff --git a/TDDFPT/src/lr_dvpsi_e.f90 b/TDDFPT/src/lr_dvpsi_e.f90
index 8a115782d..22f4406f0 100644
--- a/TDDFPT/src/lr_dvpsi_e.f90
+++ b/TDDFPT/src/lr_dvpsi_e.f90
@@ -147,7 +147,9 @@ SUBROUTINE lr_dvpsi_e(ik,ipol,dvpsi)
!
IF (okvan) THEN
ALLOCATE (spsi ( npwx*npol, nbnd))
- CALL lr_sm1_psi (.TRUE.,ik,npwx,ngk(ik),nbnd,dvpsi,spsi)
+ CALL lr_sm1_initialize()
+ CALL lr_sm1_psi_tpw (ik,npwx,ngk(ik),nbnd,dvpsi,spsi)
+! CALL lr_sm1_psi(.TRUE.,ik,npwx,ngk(ik),nbnd,dvpsi,spsi)
dvpsi(:,:) = spsi(:,:)
DEALLOCATE(spsi)
ENDIF
diff --git a/TDDFPT/src/lr_dvpsi_eels.f90 b/TDDFPT/src/lr_dvpsi_eels.f90
index 43050bb91..d27b3aea2 100644
--- a/TDDFPT/src/lr_dvpsi_eels.f90
+++ b/TDDFPT/src/lr_dvpsi_eels.f90
@@ -118,7 +118,7 @@ SUBROUTINE lr_dvpsi_eels (ik, dvpsi1, dvpsi2)
! Calculate beta-functions vkb at point k+q
CALL init_us_2(npwq, igk_k(1,ikq), xk(1,ikq), vkb)
!
- CALL lr_addus_dvpsi (ik, npwx, npwq, nbnd_occ(ikk), dvpsi1, dpsi)
+ CALL lr_addus_dvpsi ( npwq,ik, dvpsi1, dpsi)
!
dvpsi1 = dpsi
dpsi(:,:) = (0.d0, 0.d0)
@@ -139,7 +139,8 @@ SUBROUTINE lr_dvpsi_eels (ik, dvpsi1, dvpsi2)
!
dpsi(:,:) = (0.0d0, 0.0d0)
!
- CALL lr_sm1_psiq (.TRUE., ik, npwx, npwq, nbnd_occ(ikk), dvpsi1, dpsi)
+ CALL lr_sm1_initialize()
+ CALL lr_sm1_psi_tpw ( ik, npwx, npwq, nbnd_occ(ikk), dvpsi1, dpsi)
!
dvpsi1 = dpsi
!
diff --git a/TDDFPT/src/lr_readin.f90 b/TDDFPT/src/lr_readin.f90
index c07b82ca7..9fcd9c4f9 100644
--- a/TDDFPT/src/lr_readin.f90
+++ b/TDDFPT/src/lr_readin.f90
@@ -48,7 +48,6 @@ SUBROUTINE lr_readin
USE esm, ONLY : do_comp_esm
USE qpoint, ONLY : xq
USE io_rho_xml, ONLY : write_scf
- USE noncollin_module, ONLY : noncolin
USE mp_bands, ONLY : ntask_groups
USE constants, ONLY : eps4
USE control_lr, ONLY : lrpa
@@ -534,7 +533,6 @@ CONTAINS
!
IF (eels) THEN
!
- IF (okvan .AND. noncolin) CALL errore( 'lr_readin', 'Ultrasoft PP + noncolin is not fully implemented', 1 )
IF (gamma_only) CALL errore( 'lr_readin', 'gamma_only is not supported', 1 )
!
! Tamm-Dancoff approximation is not recommended to be used with EELS, and
diff --git a/TDDFPT/src/lr_restart.f90 b/TDDFPT/src/lr_restart.f90
index 83e1be32a..ea6ca6711 100644
--- a/TDDFPT/src/lr_restart.f90
+++ b/TDDFPT/src/lr_restart.f90
@@ -19,6 +19,7 @@ SUBROUTINE lr_restart(iter_restart,rflag)
USE control_flags, ONLY : gamma_only
USE klist, ONLY : nks, xk, ngk, igk_k
USE io_files, ONLY : tmp_dir, prefix, diropn, wfc_dir
+ USE uspp, ONLY : okvan
USE lr_variables, ONLY : itermax, evc1, evc1_old, &
restart, nwordrestart, iunrestart,project,nbnd_total,F, &
bgz_suffix, beta_store, gamma_store, zeta_store, norm0, &
@@ -74,6 +75,7 @@ SUBROUTINE lr_restart(iter_restart,rflag)
ENDIF
!
! Reading Lanczos coefficients
+ IF (okvan) CALL lr_sm1_initialize()
!
IF (eels) THEN
filename = trim(prefix) // trim(bgz_suffix) // trim("dat")
diff --git a/TDDFPT/src/lr_solve_e.f90 b/TDDFPT/src/lr_solve_e.f90
index cfdf04606..68c2b3c0e 100644
--- a/TDDFPT/src/lr_solve_e.f90
+++ b/TDDFPT/src/lr_solve_e.f90
@@ -27,9 +27,9 @@ SUBROUTINE lr_solve_e
USE klist, ONLY : nks, xk, ngk, igk_k, degauss
USE lr_variables, ONLY : nwordd0psi, iund0psi,LR_polarization, test_case_no, &
& n_ipol, evc0, d0psi, d0psi2, evc1, lr_verbosity, &
- & d0psi_rs, eels, lr_exx
- USE lsda_mod, ONLY : lsda, isk, current_spin
- USE uspp, ONLY : vkb
+ & d0psi_rs, eels, lr_exx,intq, intq_nc
+ USE lsda_mod, ONLY : lsda, isk, current_spin,nspin
+ USE uspp, ONLY : vkb, okvan
USE wvfct, ONLY : nbnd, npwx, et, current_k
USE control_flags, ONLY : gamma_only
USE wavefunctions, ONLY : evc
@@ -37,7 +37,10 @@ SUBROUTINE lr_solve_e
USE mp, ONLY : mp_max, mp_min, mp_barrier
USE control_lr, ONLY : alpha_pv
USE qpoint, ONLY : nksq
- USE noncollin_module, ONLY : npol
+ USE noncollin_module, ONLY : npol,noncolin
+ USE uspp_param, ONLY : nhm
+ USE ions_base, ONLY : nat
+
!
IMPLICIT NONE
INTEGER :: ibnd, ik, is, ip
@@ -55,6 +58,14 @@ SUBROUTINE lr_solve_e
IF (eels) THEN
!
! EELS case
+
+ IF (okvan) THEN
+ ALLOCATE (intq (nhm, nhm, nat) )
+ IF (noncolin) THEN
+ ALLOCATE(intq_nc( nhm, nhm, nat, nspin))
+ ENDIF
+ CALL compute_intq()
+ ENDIF
!
DO ik = 1, nksq
CALL lr_dvpsi_eels(ik, d0psi(:,:,ik,1), d0psi2(:,:,ik,1))
diff --git a/TDDFPT/src/lr_variables.f90 b/TDDFPT/src/lr_variables.f90
index 7e8f0d199..e078cd40d 100644
--- a/TDDFPT/src/lr_variables.f90
+++ b/TDDFPT/src/lr_variables.f90
@@ -75,7 +75,12 @@ MODULE lr_variables
COMPLEX(kind=dp), ALLOCATABLE :: R(:,:,:) ! The oscillator strength from valence state (first index)
! to conduction state (second index), for each polarization
! direction (third index).
- !
+
+
+
+ COMPLEX (DP), ALLOCATABLE :: &
+ intq(:,:,:), &! nhm, nhm, nat), integral of e^iqr Q
+ intq_nc(:,:,:,:) ! nhm, nhm, nat, nspin), integral of
! Lanczos Matrix
!
!
diff --git a/TDDFPT/src/make.depend b/TDDFPT/src/make.depend
index f38f16c86..848d52961 100644
--- a/TDDFPT/src/make.depend
+++ b/TDDFPT/src/make.depend
@@ -569,3 +569,24 @@ stop_lr.o : ../../Modules/noncol.o
stop_lr.o : ../../PW/src/buffers.o
stop_lr.o : ../../PW/src/pwcom.o
stop_lr.o : lr_variables.o
+compute_intq.o : ../../LR_Modules/lrcom.o
+compute_intq.o : ../../Modules/cell_base.o
+compute_intq.o : ../../Modules/ions_base.o
+compute_intq.o : ../../Modules/kind.o
+compute_intq.o : ../../Modules/noncol.o
+compute_intq.o : ../../Modules/uspp.o
+compute_intq.o : lr_variables.o
+set_intq_nc.o : ../../Modules/ions_base.o
+set_intq_nc.o : ../../Modules/uspp.o
+set intq_nc.o : lr_variables.o
+transform_intq_nc.o : ../../Modules/ions_base.o
+transform_intq_nc.o : ../../Modules/kind.o
+transform_intq_nc.o : ../../Modules/uspp.o
+transform_intq_nc.o : lr_variables.o
+transform_intq_so.o : ../../Modules/ions_base.o
+transform_intq_so.o : ../../Modules/kind.o
+transform_intq_so.o : ../../Modules/noncol.o
+transform_intq_so.o : ../../Modules/uspp.o
+transform_intq_so.o : ../../PW/src/pwcom.o
+transform_intq_so.o : lr_variables.o
+
diff --git a/TDDFPT/src/set_intq_nc.f90 b/TDDFPT/src/set_intq_nc.f90
new file mode 100644
index 000000000..d2037e42d
--- /dev/null
+++ b/TDDFPT/src/set_intq_nc.f90
@@ -0,0 +1,35 @@
+!
+! Copyright (C) 2001-2019 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,
+! or http://www.gnu.org/copyleft/gpl.txt .
+
+!----------------------------------------------------------------------------
+SUBROUTINE set_intq_nc()
+!----------------------------------------------------------------------------
+USE ions_base, ONLY : nat, ntyp => nsp, ityp
+USE uspp_param, ONLY : upf
+USE lr_variables, ONLY : intq, intq_nc
+IMPLICIT NONE
+INTEGER :: npe
+INTEGER :: np, na
+
+intq_nc=(0.d0,0.d0)
+DO np = 1, ntyp
+ IF ( upf(np)%tvanp ) THEN
+ DO na = 1, nat
+ IF (ityp(na)==np) THEN
+ IF (upf(np)%has_so) THEN
+ CALL transform_intq_so(intq,na)
+ ELSE
+ CALL transform_intq_nc(intq,na)
+ END IF
+ END IF
+ END DO
+ END IF
+END DO
+
+RETURN
+END SUBROUTINE set_intq_nc
+!
diff --git a/TDDFPT/src/transform_intq_nc.f90 b/TDDFPT/src/transform_intq_nc.f90
new file mode 100644
index 000000000..684663df1
--- /dev/null
+++ b/TDDFPT/src/transform_intq_nc.f90
@@ -0,0 +1,39 @@
+!
+! Copyright (C) 2001-2019 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,
+! or http://www.gnu.org/copyleft/gpl.txt .
+!
+
+!----------------------------------------------------------------------------
+SUBROUTINE transform_intq_nc(intq,na)
+!----------------------------------------------------------------------------
+!
+! This routine multiply intq by the identity and the Pauli
+! matrices and saves it in intq_nc.
+!
+USE kinds, ONLY : DP
+USE ions_base, ONLY : nat, ityp
+USE uspp_param, ONLY : nh, nhm
+USE lr_variables, ONLY : intq_nc
+!
+IMPLICIT NONE
+
+INTEGER :: na
+COMPLEX(DP) :: intq(nhm,nhm,nat)
+!
+! ... local variables
+!
+INTEGER :: ih, jh, np
+
+np=ityp(na)
+DO ih = 1, nh(np)
+ DO jh = 1, nh(np)
+ intq_nc(ih,jh,na,1)=intq(ih,jh,na)
+ intq_nc(ih,jh,na,4)=intq(ih,jh,na)
+ END DO
+END DO
+
+RETURN
+END SUBROUTINE transform_intq_nc
diff --git a/TDDFPT/src/transform_intq_so.f90 b/TDDFPT/src/transform_intq_so.f90
new file mode 100644
index 000000000..f2e85615b
--- /dev/null
+++ b/TDDFPT/src/transform_intq_so.f90
@@ -0,0 +1,58 @@
+!
+! Copyright (C) 2001-2019 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,
+! or http://www.gnu.org/copyleft/gpl.txt .
+!
+!----------------------------------------------------------------------------
+SUBROUTINE transform_intq_so(intq,na)
+!----------------------------------------------------------------------------
+!
+! This routine multiply intq by the identity and the Pauli
+! matrices, rotate it as appropriate for the spin-orbit case
+! and saves it in intq_nc.
+!
+USE kinds, ONLY : DP
+USE ions_base, ONLY : nat, ityp
+USE uspp_param, ONLY : nh, nhm
+USE noncollin_module, ONLY : npol, nspin_mag
+USE spin_orb, ONLY : fcoef, domag
+USE lr_variables, ONLY : intq_nc
+!
+IMPLICIT NONE
+
+COMPLEX(DP) :: intq(nhm,nhm,nat)
+INTEGER :: na
+!
+! ... local variables
+!
+INTEGER :: ih, jh, lh, kh, np, npert, is1, is2, ijs
+LOGICAL :: same_lj
+
+np=ityp(na)
+DO ih = 1, nh(np)
+ DO kh = 1, nh(np)
+ IF (same_lj(kh,ih,np)) THEN
+ DO jh = 1, nh(np)
+ DO lh= 1, nh(np)
+ IF (same_lj(lh,jh,np)) THEN
+ ijs=0
+ DO is1=1,npol
+ DO is2=1,npol
+ ijs=ijs+1
+ intq_nc(ih,jh,na,ijs)= &
+ intq_nc(ih,jh,na,ijs) + intq (kh,lh,na)* &
+ (fcoef(ih,kh,is1,1,np)*fcoef(lh,jh,1,is2,np) + &
+ fcoef(ih,kh,is1,2,np)*fcoef(lh,jh,2,is2,np) )
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDIF
+ ENDDO
+ENDDO
+ !
+RETURN
+END SUBROUTINE transform_intq_so
From e9fc7aebefdc2a0d0db0319e72bcfbc16c9b8544 Mon Sep 17 00:00:00 2001
From: Pietro Delugas
Date: Fri, 3 May 2019 16:50:26 +0200
Subject: [PATCH 3/4] updated dependencies and fixed few things
* in TDDFPT/src/lr_alloc_init.f90 use bbg allocatable from lrus
* in TDDFPT/src/dveqpsi_us_only.f90 commented the use of inexistent
optical module
---
LR_Modules/make.depend | 2 +-
TDDFPT/src/dveqpsi_us_only.f90 | 2 +-
TDDFPT/src/lr_alloc_init.f90 | 2 +-
TDDFPT/src/make.depend | 32 ++++++++++++++++++--------------
4 files changed, 21 insertions(+), 17 deletions(-)
diff --git a/LR_Modules/make.depend b/LR_Modules/make.depend
index 2d776b8af..088f05044 100644
--- a/LR_Modules/make.depend
+++ b/LR_Modules/make.depend
@@ -261,7 +261,7 @@ lr_sm1_psi.o : ../Modules/control_flags.o
lr_sm1_psi.o : ../Modules/invmat.o
lr_sm1_psi.o : ../Modules/ions_base.o
lr_sm1_psi.o : ../Modules/kind.o
-lr_sm1_psi.o : ../Modules/mp_bands.o
+lr_sm1_psi.o : ../Modules/mp_global.o
lr_sm1_psi.o : ../Modules/noncol.o
lr_sm1_psi.o : ../Modules/uspp.o
lr_sm1_psi.o : ../PW/src/pwcom.o
diff --git a/TDDFPT/src/dveqpsi_us_only.f90 b/TDDFPT/src/dveqpsi_us_only.f90
index 5704a2cea..e19ceee95 100644
--- a/TDDFPT/src/dveqpsi_us_only.f90
+++ b/TDDFPT/src/dveqpsi_us_only.f90
@@ -23,7 +23,7 @@ SUBROUTINE dveqpsi_us_only (npwq, ik)
USE noncollin_module, ONLY : noncolin, npol
! modules from phcom
USE qpoint, ONLY : ikks
- USE optical, ONLY : intq, intq_nc
+! USE optical, ONLY : intq, intq_nc
USE lrus, ONLY : becp1
USE eqv, ONLY : dvpsi
IMPLICIT NONE
diff --git a/TDDFPT/src/lr_alloc_init.f90 b/TDDFPT/src/lr_alloc_init.f90
index b8359ee1e..b33064ce5 100644
--- a/TDDFPT/src/lr_alloc_init.f90
+++ b/TDDFPT/src/lr_alloc_init.f90
@@ -15,7 +15,7 @@ SUBROUTINE lr_alloc_init()
USE kinds, ONLY : dp
USE ions_base, ONLY : nat
USE uspp, ONLY : nkb, okvan
- USE lrus, ONLY : bbk, bbnc_sm1
+ USE lrus, ONLY : bbk, bbg, bbnc_sm1
USE uspp_param, ONLY : nhm
USE fft_base, ONLY : dfftp, dffts
USE klist, ONLY : nks
diff --git a/TDDFPT/src/make.depend b/TDDFPT/src/make.depend
index 848d52961..77c1e3488 100644
--- a/TDDFPT/src/make.depend
+++ b/TDDFPT/src/make.depend
@@ -8,13 +8,26 @@ bcast_lr_input.o : ../../UtilXlib/mp.o
bcast_lr_input.o : lr_charg_resp.o
bcast_lr_input.o : lr_dav_variables.o
bcast_lr_input.o : lr_variables.o
+compute_intq.o : ../../LR_Modules/lrcom.o
+compute_intq.o : ../../Modules/cell_base.o
+compute_intq.o : ../../Modules/ions_base.o
+compute_intq.o : ../../Modules/kind.o
+compute_intq.o : ../../Modules/noncol.o
+compute_intq.o : ../../Modules/uspp.o
+compute_intq.o : lr_variables.o
+dveqpsi_us_only.o : ../../LR_Modules/lrcom.o
+dveqpsi_us_only.o : ../../Modules/ions_base.o
+dveqpsi_us_only.o : ../../Modules/kind.o
+dveqpsi_us_only.o : ../../Modules/noncol.o
+dveqpsi_us_only.o : ../../Modules/uspp.o
+dveqpsi_us_only.o : ../../PW/src/pwcom.o
lr_addus_dvpsi.o : ../../LR_Modules/lrcom.o
-lr_addus_dvpsi.o : ../../Modules/cell_base.o
lr_addus_dvpsi.o : ../../Modules/ions_base.o
lr_addus_dvpsi.o : ../../Modules/kind.o
lr_addus_dvpsi.o : ../../Modules/noncol.o
-lr_addus_dvpsi.o : ../../Modules/paw_variables.o
lr_addus_dvpsi.o : ../../Modules/uspp.o
+lr_addus_dvpsi.o : ../../PW/src/pwcom.o
+lr_addus_dvpsi.o : lr_variables.o
lr_alloc_init.o : ../../LR_Modules/lrcom.o
lr_alloc_init.o : ../../Modules/becmod.o
lr_alloc_init.o : ../../Modules/control_flags.o
@@ -407,7 +420,6 @@ lr_readin.o : ../../Modules/io_global.o
lr_readin.o : ../../Modules/kind.o
lr_readin.o : ../../Modules/mp_bands.o
lr_readin.o : ../../Modules/mp_global.o
-lr_readin.o : ../../Modules/noncol.o
lr_readin.o : ../../Modules/paw_variables.o
lr_readin.o : ../../Modules/recvec.o
lr_readin.o : ../../Modules/uspp.o
@@ -558,6 +570,9 @@ sd0psi.o : ../../Modules/uspp.o
sd0psi.o : ../../PW/src/pwcom.o
sd0psi.o : ../../PW/src/realus.o
sd0psi.o : lr_variables.o
+set_intq_nc.o : ../../Modules/ions_base.o
+set_intq_nc.o : ../../Modules/uspp.o
+set_intq_nc.o : lr_variables.o
stop_lr.o : ../../Modules/cell_base.o
stop_lr.o : ../../Modules/environment.o
stop_lr.o : ../../Modules/io_files.o
@@ -569,16 +584,6 @@ stop_lr.o : ../../Modules/noncol.o
stop_lr.o : ../../PW/src/buffers.o
stop_lr.o : ../../PW/src/pwcom.o
stop_lr.o : lr_variables.o
-compute_intq.o : ../../LR_Modules/lrcom.o
-compute_intq.o : ../../Modules/cell_base.o
-compute_intq.o : ../../Modules/ions_base.o
-compute_intq.o : ../../Modules/kind.o
-compute_intq.o : ../../Modules/noncol.o
-compute_intq.o : ../../Modules/uspp.o
-compute_intq.o : lr_variables.o
-set_intq_nc.o : ../../Modules/ions_base.o
-set_intq_nc.o : ../../Modules/uspp.o
-set intq_nc.o : lr_variables.o
transform_intq_nc.o : ../../Modules/ions_base.o
transform_intq_nc.o : ../../Modules/kind.o
transform_intq_nc.o : ../../Modules/uspp.o
@@ -589,4 +594,3 @@ transform_intq_so.o : ../../Modules/noncol.o
transform_intq_so.o : ../../Modules/uspp.o
transform_intq_so.o : ../../PW/src/pwcom.o
transform_intq_so.o : lr_variables.o
-
From 0b0d934815b26ce0f35490903a7828a856db9744 Mon Sep 17 00:00:00 2001
From: Nikolas Garofil
Date: Mon, 6 May 2019 17:56:21 +0200
Subject: [PATCH 4/4] fix fetching Wannier90
---
install/plugins_list | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/install/plugins_list b/install/plugins_list
index 714528ce2..19e03358e 100644
--- a/install/plugins_list
+++ b/install/plugins_list
@@ -28,7 +28,7 @@ YAMBO_URL=http://www.yambo-code.org/files/yambo-stable-latest.tar.gz
#
# Package maintainer:
W90=wannier90-3.0.0
-W90_URL=https://github.com/wannier-developers/wannier90/archive/v3.0.0.tar.gz
+W90_URL=https://codeload.github.com/wannier-developers/wannier90/tar.gz/v3.0.0
#
# Package maintainer: Layla Martin-Samos
SAX=sax-2.0.3