Commit 16edbc1a authored by dalcorso's avatar dalcorso

Phonon in the noncollinear and spin-orbit case. New routines.


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@3767 c92efa57-630b-4861-b058-cf58834340f0
parent 69c1fd13
......@@ -13,6 +13,7 @@ adddvscf.o \
addnlcc.o \
addnlcc_zstar_eu_us.o \
addusdbec.o \
addusdbec_nc.o \
addusddens.o \
addusddense.o \
addusdynmat.o \
......@@ -32,6 +33,7 @@ compute_alphasum.o \
compute_becalp.o \
compute_becsum.o \
compute_drhous.o \
compute_drhous_nc.o \
compute_dvloc.o \
compute_nldyn.o \
compute_weight.o \
......@@ -68,7 +70,9 @@ generate_dynamical_matrix_c.o \
gmressolve_all.o \
h_psiq.o \
incdrhoscf.o \
incdrhoscf_nc.o \
incdrhous.o \
incdrhous_nc.o \
io_pattern.o \
localdos.o \
newdq.o \
......@@ -85,6 +89,7 @@ print_clock_ph.o \
psidspsi.o \
psymdvscf.o \
psyme.o \
psym_dmag.o \
punch_plot_e.o \
punch_plot_ph.o \
q_points.o \
......@@ -93,6 +98,7 @@ random_matrix.o \
rotate_and_add_dyn.o \
set_asr_c.o \
set_drhoc.o \
set_int12_nc.o \
set_irr.o \
set_irr_mode.o \
set_irr_nosym.o \
......@@ -107,11 +113,18 @@ star_q.o \
stop_ph.o \
sym_and_write_zue.o \
sym_def.o \
sym_dmag.o \
symdvscf.o \
symdyn_munu.o \
symdynph_gq.o \
syme.o \
symm.o \
transform_int_so.o \
transform_int_nc.o \
transform_alphasum_nc.o \
transform_alphasum_so.o \
transform_dbecsum_so.o \
transform_dbecsum_nc.o \
tra_write_matrix.o \
trntnsc.o \
xk_wk_collect.o \
......
!
! Copyright (C) 2007 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 addusdbec_nc (ik, wgt, psi, dbecsum_nc)
!----------------------------------------------------------------------
!
! This routine adds to the dbecsum the term which correspond to this
! k point. After the accumulation the additional part of the charge
! is computed in addusddens.
!
#include "f_defs.h"
USE ions_base, ONLY : nat, ityp, ntyp => nsp
use pwcom
USE noncollin_module, ONLY : noncolin, npol
USE kinds, only : DP
USE uspp_param, only: nh, tvanp, nhm
use phcom
implicit none
!
! the dummy variables
!
complex(DP) :: dbecsum_nc (nhm,nhm,nat,nspin), psi(npwx*npol,nbnd)
! inp/out: the sum kv of bec *
! input : contains delta psi
integer :: ik
! input: the k point
real(DP) :: wgt
! input: the weight of this k point
!
! here the local variables
!
integer :: na, nt, ih, jh, ibnd, ikk, ikb, jkb, startb, &
lastb, ijkb0, is1, is2, ijs
! counter on atoms
! counter on atomic type
! counter on solid beta functions
! counter on solid beta functions
! counter on the bands
! the real k point
! counter on solid becp
! counter on solid becp
! composite index for dbecsum
! divide among processors the sum
! auxiliary variable for counting
real(DP) :: w
! the weight of the k point
complex(DP), allocatable :: dbecq_nc(:,:,:)
! the change of becq
if (.not.okvan) return
call start_clock ('addusdbec_nc')
allocate (dbecq_nc( nkb,npol, nbndx))
w = 0.5d0 * wgt * omega
if (lgamma) then
ikk = ik
else
ikk = 2 * ik - 1
endif
!
! First compute the product of psi and vkb
!
call ccalbec_nc (nkb, npwx, npwq, npol, nbnd, dbecq_nc, vkb, psi)
!
! And then we add the product to becsum
!
! Band parallelization: each processor takes care of its slice of bands
!
call divide (nbnd_occ (ikk), startb, lastb)
!
ijkb0 = 0
do nt = 1, ntyp
if (tvanp (nt) ) then
do na = 1, nat
if (ityp (na) .eq.nt) then
do ih = 1, nh (nt)
ikb = ijkb0 + ih
do jh = 1, nh (nt)
jkb = ijkb0 + jh
DO ibnd = startb, lastb
ijs=0
DO is1=1,npol
DO is2=1,npol
ijs=ijs+1
dbecsum_nc(ih,jh,na,ijs)=dbecsum_nc(ih,jh,na,ijs)+&
w*CONJG(becp1_nc(ikb,is1,ibnd,ik)) &
*dbecq_nc(jkb,is2,ibnd)
ENDDO
ENDDO
ENDDO
enddo
enddo
ijkb0 = ijkb0 + nh (nt)
endif
enddo
else
do na = 1, nat
if (ityp (na) .eq.nt) ijkb0 = ijkb0 + nh (nt)
enddo
endif
enddo
!
deallocate (dbecq_nc)
call stop_clock ('addusdbec_nc')
return
end subroutine addusdbec_nc
!
! Copyright (C) 2001 PWSCF 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 .
!
#include "f_defs.h"
!
!-----------------------------------------------------------------------
subroutine compute_drhous_nc (drhous, dbecsum, wgg, becq, alpq)
!-----------------------------------------------------------------------
!
! This routine computes the part of the change of the charge density
! which is due to the orthogonalization constraint on wavefunctions
!
!
!
USE ions_base, ONLY : nat
use pwcom
USE noncollin_module, ONLY : noncolin, npol
USE wavefunctions_module, ONLY: evc
USE io_files, ONLY: iunigk
USE kinds, only : DP
USE uspp_param, ONLY: nhm
use phcom
implicit none
!
! the dummy variables
!
complex(DP) :: dbecsum (nhm, nhm, nat, nspin, 3 * nat), &
drhous (nrxx, nspin, 3 * nat), becq (nkb, npol, nbnd, nksq), &
alpq (nkb, npol, nbnd, 3, nksq)
!output:the derivative of becsum
! output: add the orthogonality term
! input: the becp with psi_{k+q}
! input: the alphap with psi_{k+q}
real(DP) :: wgg (nbnd, nbnd, nksq)
! input: the weights
integer :: ipert, mode, ik, ikq, ikk, is, ig, nu_i, ibnd, ios
! counter on the pertubations
! counter on the modes
! counter on k points
! the point k+q
! record for wfcs at k point
! counter on spin
! counter on g vectors
! counter on modes
! counter on the bands
! integer variable for I/O control
real(DP) :: weight
! the weight of the k point
complex(DP), allocatable :: evcr (:,:,:)
! the wavefunctions in real space
if (.not.okvan) return
call start_clock ('com_drhous')
allocate (evcr( nrxxs, npol, nbnd))
!
drhous(:,:,:) = (0.d0, 0.d0)
dbecsum = (0.d0, 0.d0)
if (nksq.gt.1) rewind (unit = iunigk)
do ik = 1, nksq
if (nksq.gt.1) then
read (iunigk, err = 110, iostat = ios) npw, igk
110 call errore ('compute_drhous', 'reading igk', abs (ios) )
endif
if (lgamma) then
ikk = ik
ikq = ik
npwq = npw
else
ikk = 2 * ik - 1
ikq = ikk + 1
endif
weight = wk (ikk)
if (lsda) current_spin = isk (ikk)
if (.not.lgamma.and.nksq.gt.1) then
read (iunigk, err = 210, iostat = ios) npwq, igkq
210 call errore ('compute_drhous', 'reading igkq', abs (ios) )
endif
!
! For each k point we construct the beta functions
!
call init_us_2 (npwq, igkq, xk (1, ikq), vkb)
!
! Read the wavefunctions at k and transform to real space
!
call davcio (evc, lrwfc, iuwfc, ikk, - 1)
evcr = (0.d0, 0.d0)
do ibnd = 1, nbnd
do ig = 1, npw
evcr (nls (igk (ig) ), 1, ibnd) = evc (ig, ibnd)
evcr (nls (igk (ig) ), 2, ibnd) = evc (ig+npwx, ibnd)
enddo
call cft3s (evcr (1, 1, ibnd), nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s,+2)
call cft3s (evcr (1, 2, ibnd),nr1s,nr2s,nr3s,nrx1s,nrx2s,nrx3s,+2)
enddo
!
! Read the wavefunctions at k+q
!
if (.not.lgamma.and.nksq.gt.1) call davcio (evq, lrwfc, iuwfc, ikq, -1)
!
! And compute the contribution of this k point to the change of
! the charge density
!
do nu_i = 1, 3 * nat
call incdrhous_nc (drhous (1, 1, nu_i), weight, ik, &
dbecsum (1, 1, 1, 1, nu_i), evcr, wgg, becq, alpq, nu_i)
enddo
enddo
deallocate(evcr)
call stop_clock ('com_drhous')
return
end subroutine compute_drhous_nc
!
! Copyright (C) 2001 PWSCF 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 incdrhoscf_nc (drhoscf, weight, ik, dbecsum)
!-----------------------------------------------------------------------
!
! This routine computes the change of the charge density due to the
! perturbation. It is called at the end of the computation of the
! change of the wavefunction for a given k point.
!
!
#include "f_defs.h"
USE ions_base, ONLY : nat
use pwcom
USE noncollin_module, ONLY : noncolin, npol
USE wavefunctions_module, ONLY: evc
USE kinds, only : DP
USE uspp_param, ONLY: nhm
use phcom
implicit none
integer :: ik
! input: the k point
real(DP) :: weight
! input: the weight of the k point
complex(DP) :: drhoscf (nrxx,nspin), dbecsum (nhm,nhm,nat,nspin)
! output: the change of the charge densit
! inp/out: the accumulated dbec
!
! here the local variable
!
real(DP) :: wgt
! the effective weight of the k point
complex(DP), allocatable :: psi (:,:), dpsic (:,:)
! the wavefunctions in real space
! the change of wavefunctions in real space
integer :: ibnd, jbnd, ikk, ir, ig
! counters
call start_clock ('incdrhoscf')
allocate (dpsic( nrxxs, npol))
allocate (psi ( nrxxs, npol))
wgt = 2.d0 * weight / omega
if (lgamma) then
ikk = ik
else
ikk = 2 * ik - 1
endif
!
! dpsi contains the perturbed wavefunctions of this k point
! evc contains the unperturbed wavefunctions of this k point
!
do ibnd = 1, nbnd_occ (ikk)
psi = (0.d0, 0.d0)
do ig = 1, npw
psi (nls (igk (ig) ), 1) = evc (ig, ibnd)
psi (nls (igk (ig) ), 2) = evc (ig+npwx, ibnd)
enddo
call cft3s (psi, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, + 2)
call cft3s (psi(1,2), nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, + 2)
dpsic = (0.d0, 0.d0)
do ig = 1, npwq
dpsic (nls (igkq (ig)), 1 ) = dpsi (ig, ibnd)
dpsic (nls (igkq (ig)), 2 ) = dpsi (ig+npwx, ibnd)
enddo
call cft3s (dpsic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, + 2)
call cft3s (dpsic(1,2), nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, + 2)
do ir = 1, nrxxs
drhoscf(ir,1)=drhoscf(ir,1)+wgt*(CONJG(psi(ir,1))*dpsic(ir,1) + &
CONJG(psi(ir,2))*dpsic(ir,2) )
enddo
IF (domag) THEN
do ir = 1, nrxxs
drhoscf(ir,2)=drhoscf (ir,2) + wgt *(CONJG(psi(ir,1))*dpsic(ir,2)+ &
CONJG(psi(ir,2))*dpsic(ir,1) )
drhoscf(ir,3)=drhoscf (ir,3) + wgt *(CONJG(psi(ir,1))*dpsic(ir,2)- &
CONJG(psi(ir,2))*dpsic(ir,1) ) * (0.d0,-1.d0)
drhoscf(ir,4)=drhoscf (ir,4) + wgt *(CONJG(psi(ir,1))*dpsic(ir,1)- &
CONJG(psi(ir,2))*dpsic(ir,2) )
enddo
END IF
enddo
call addusdbec_nc (ik, wgt, dpsi, dbecsum)
deallocate (psi)
deallocate (dpsic)
call stop_clock ('incdrhoscf')
return
end subroutine incdrhoscf_nc
!
! Copyright (C) 2001 PWSCF 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 incdrhous_nc (drhoscf, weight, ik, dbecsum, evcr, wgg, becq, &
alpq, mode)
!-----------------------------------------------------------------------
!
! This routine computes the change of the charge density
! and of the magnetization due
! to the displacement of the augmentation charge. Only the
! smooth part is computed here.
!
#include "f_defs.h"
USE ions_base, ONLY : ntyp => nsp, nat, ityp
use pwcom
USE noncollin_module, ONLY : npol
USE kinds, only : DP
USE uspp_param, ONLY: nhm, nh
use phcom
implicit none
integer :: ik, mode
! input: the k point
! input: the mode which is computed
! input: the quantity to compute (1 charge, 2-4 magnetization)
real(DP) :: weight, wgg (nbnd, nbnd, nksq)
! input: the weight of the k point
! input: the weights
complex(DP) :: evcr (nrxxs, npol, nbnd), drhoscf (nrxx,nspin), &
dbecsum(nhm, nhm, nat, nspin), becq (nkb, npol, nbnd, nksq), &
alpq (nkb, npol, nbnd, 3, nksq)
! input: the wavefunctions at k in real
! output: the change of the charge densi
! inp/out: the accumulated dbec
! input: the becp with psi_{k+q}
! input: the alphap with psi_{k+
!
! here the local variable
!
real(DP) :: wgt
! the effective weight of the k point
complex(DP), allocatable :: ps1 (:,:), dpsir (:,:)
! auxiliary space
! the change of wavefunctions in real sp
integer :: ibnd, jbnd, nt, na, mu, ih, jh, ikb, jkb, ijkb0, &
startb, lastb, ipol, ikk, ir, ig, ijs, is1, is2
! counters
call start_clock ('incdrhous')
allocate (dpsir(nrxxs,npol))
allocate (ps1 (nbnd, nbnd))
call divide (nbnd, startb, lastb)
ps1 (:,:) = (0.d0, 0.d0)
if (lgamma) then
ikk = ik
else
ikk = 2 * ik - 1
endif
!
! Here we prepare the two terms
!
ijkb0 = 0
do nt = 1, ntyp
do na = 1, nat
if (ityp (na) == nt) then
mu = 3 * (na - 1)
if (abs(u(mu+1,mode)) + abs(u(mu+2,mode)) &
+ abs(u(mu+3,mode)) > 1.0d-12) then
do ih = 1, nh (nt)
ikb = ijkb0 + ih
do jh = 1, nh (nt)
jkb = ijkb0 + jh
do ibnd = 1, nbnd
do jbnd = startb, lastb
do ipol = 1, 3
mu = 3 * (na - 1) + ipol
IF (lspinorb) THEN
ijs=0
DO is1=1,npol
DO is2=1,npol
ijs=ijs+1
ps1(ibnd,jbnd)=ps1(ibnd,jbnd)- &
qq_so(ih,jh,ijs,nt) * &
(alphap_nc(jkb,is2,ibnd,ipol,ik)*&
CONJG(becq(ikb,is1,jbnd,ik)) + &
becp1_nc(jkb,is2,ibnd,ik) * &
CONJG(alpq(ikb,is1,jbnd,ipol,ik)) )* &
wgg (ibnd, jbnd, ik) * u (mu, mode)
END DO
END DO
ELSE
ps1(ibnd,jbnd)=ps1(ibnd,jbnd)-qq(ih,jh,nt)*&
(alphap_nc(ikb,1,ibnd,ipol,ik) *&
CONJG(becq(jkb,1,jbnd,ik)) + &
becp1_nc(ikb,1,ibnd,ik) * &
CONJG(alpq(jkb,1,jbnd,ipol,ik)) + &
alphap_nc(ikb,2,ibnd,ipol,ik) *&
CONJG(becq(jkb,2,jbnd,ik)) + &
becp1_nc(ikb,2,ibnd,ik) * &
CONJG(alpq(jkb,2,jbnd,ipol,ik)) )* &
wgg (ibnd, jbnd, ik) * u (mu, mode)
END IF
enddo
enddo
enddo
enddo
enddo
endif
ijkb0 = ijkb0 + nh (nt)
endif
enddo
enddo
#ifdef __PARA
call reduce (2 * nbnd * nbnd, ps1)
#endif
dpsi (:,:) = (0.d0, 0.d0)
wgt = 2.d0 * weight / omega
do ibnd = 1, nbnd_occ (ikk)
do jbnd = 1, nbnd
call ZAXPY (npwx*npol,ps1(ibnd,jbnd),evq(1,jbnd),1,dpsi(1,ibnd), 1)
enddo
dpsir = (0.d0, 0.d0)
do ig = 1, npwq
dpsir(nls(igkq(ig)),1) = dpsi (ig, ibnd)
dpsir(nls(igkq(ig)),2) = dpsi (ig+npwx, ibnd)
enddo
call cft3s (dpsir, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, + 2)
call cft3s (dpsir(1,2), nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, + 2)
do ir = 1, nrxxs
drhoscf(ir,1)=drhoscf(ir,1)+wgt* &
(dpsir(ir,1)*CONJG(evcr(ir,1,ibnd))+ &
dpsir(ir,2)*CONJG(evcr(ir,2,ibnd)) )
IF (domag) THEN
drhoscf(ir,2)=drhoscf(ir,2)+ &
wgt*(dpsir(ir,1)*CONJG(evcr(ir,2,ibnd))+ &
dpsir(ir,2)*CONJG(evcr(ir,1,ibnd)))
drhoscf(ir,3)=drhoscf(ir,3)+ &
wgt*(dpsir(ir,1)*CONJG(evcr(ir,2,ibnd)) - &
dpsir(ir,2)*CONJG(evcr(ir,1,ibnd)) ) *(0.d0,-1.d0)
drhoscf(ir,4)=drhoscf(ir,4)+wgt* &
(dpsir(ir,1)*CONJG(evcr(ir,1,ibnd)) - &
dpsir(ir,2)*CONJG(evcr(ir,2,ibnd)) )
END IF
enddo
enddo
call addusdbec_nc (ik, wgt, dpsi, dbecsum)
deallocate (ps1)
deallocate (dpsir)
call stop_clock ('incdrhous')
return
end subroutine incdrhous_nc
!
! Copyright (C) 2001 PWSCF 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 .
!
#include "f_defs.h"
!
!-----------------------------------------------------------------------
SUBROUTINE psym_dmag (nper, irr, dvtosym)
!-----------------------------------------------------------------------
!
! ... p-symmetrize the charge density.
!
USE pwcom
USE kinds, ONLY : DP
USE phcom
USE mp_global, ONLY : me_pool
USE pfft, ONLY : npp, ncplane
!
IMPLICIT NONE
!
INTEGER :: nper, irr
! the number of perturbations
! the representation under consideration
COMPLEX(DP) :: dvtosym (nrxx, nspin, nper)
! the potential to symmetrize
!-local variable
!
#if defined (__PARA)
!
INTEGER :: i, is, iper, npp0
COMPLEX(DP), ALLOCATABLE :: ddvtosym (:,:,:)
! the potential to symm
IF (nsymq.EQ.1.AND. (.NOT.minus_q) ) RETURN
CALL start_clock ('psym_dmag')
ALLOCATE (ddvtosym ( nrx1 * nrx2 * nrx3, nspin, nper))
npp0 = 1
DO i = 1, me_pool
npp0 = npp0 + npp (i) * ncplane
ENDDO
DO iper = 1, nper
DO is = 1, nspin
CALL cgather_sym (dvtosym (1, is, iper), ddvtosym (1, is, iper) )
ENDDO
ENDDO
CALL sym_dmag (nper, irr, ddvtosym)
DO iper = 1, nper
DO is = 1, nspin
CALL ZCOPY (npp (me_pool+1) * ncplane, ddvtosym (npp0, is, iper), &
1, dvtosym (1, is, iper), 1)
ENDDO
ENDDO
DEALLOCATE (ddvtosym)
CALL stop_clock ('psym_dmag')
#endif
RETURN
END SUBROUTINE psym_dmag
!
! Copyright (C) 2007 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_int12_nc(iflag)
!----------------------------------------------------------------------------
!
! This is a driver to call the routines that rotate and multiply
! by the Pauli matrices the integrals.
!
USE ions_base, ONLY : nat, ntyp => nsp, ityp
USE spin_orb, ONLY : lspinorb, so
USE uspp_param, only: tvanp
USE phus, ONLY : int1, int2, int1_nc, int2_so
IMPLICIT NONE
INTEGER :: iflag
INTEGER :: np, na
int1_nc=(0.d0,0.d0)
IF (lspinorb) int2_so=(0.d0,0.d0)
DO np = 1, ntyp
IF ( tvanp(np) ) THEN
DO na = 1, nat
IF (ityp(na)==np) THEN
IF (so(np)) THEN
CALL transform_int1_so(int1,na,iflag)
CALL transform_int2_so(int2,na,iflag)
ELSE
CALL transform_int1_nc(int1,na,iflag)
IF (lspinorb) CALL transform_int2_nc(int2,na,iflag)
END IF
END IF
END DO
END IF
END DO
END SUBROUTINE set_int12_nc
!----------------------------------------------------------------------------
SUBROUTINE set_int3_nc(npe)
!----------------------------------------------------------------------------
USE ions_base, ONLY : nat, ntyp => nsp, ityp
USE spin_orb, ONLY : so
USE uspp_param, only: tvanp
USE phus, ONLY : int3, int3_nc
IMPLICIT NONE
INTEGER :: npe
INTEGER :: np, na
int3_nc=(0.d0,0.d0)
DO np = 1, ntyp
IF