scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Phonon in the noncollinear and spin-orbit case. Changes in PW. 

angle1, angle2 and starting_magnetization are saved in the punch file.
The transformation of angle1 and angle2 to radiants is done in input.f90.
Clean_up of sum_band.


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@3765 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
dalcorso 2007-02-08 12:47:41 +00:00
parent f4b24be5e8
commit cfef6c3be5
9 changed files with 388 additions and 135 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
Modules/xml_io_base.f90
View File

@ -36,7 +36,8 @@ MODULE xml_io_base
save_print_counter, read_print_counter, set_kpoints_vars, & save_print_counter, read_print_counter, set_kpoints_vars, &
write_header, & write_header, &
write_cell, write_ions, write_symmetry, write_planewaves, & write_cell, write_ions, write_symmetry, write_planewaves, &
write_efield, write_spin, write_xc, write_occ, write_bz, & write_efield, write_spin, write_init_mag, write_xc, &
write_occ, write_bz, &
write_phonon, write_rho_xml, write_wfc, write_eig, & write_phonon, write_rho_xml, write_wfc, write_eig, &
read_wfc, read_rho_xml read_wfc, read_rho_xml
! !
@ -921,6 +922,34 @@ MODULE xml_io_base
END SUBROUTINE write_spin END SUBROUTINE write_spin
! !
!------------------------------------------------------------------------ !------------------------------------------------------------------------
SUBROUTINE write_init_mag(starting_magnetization, angle1, angle2, ntyp )
!------------------------------------------------------------------------
USE constants, ONLY : pi
IMPLICIT NONE
INTEGER, INTENT(IN):: ntyp
REAL(DP), INTENT(IN) :: starting_magnetization(ntyp), &
angle1(ntyp), angle2(ntyp)
INTEGER :: ityp
!
CALL iotk_write_begin( iunpun, "STARTING_MAG" )
CALL iotk_write_dat( iunpun, "NTYP", ntyp)
DO ityp=1,ntyp
CALL iotk_write_dat( iunpun, "STARTING_MAGNETIZATION", &
starting_magnetization(ityp) )
CALL iotk_write_dat( iunpun, "ANGLE1", &
angle1(ityp)*180.d0/pi )
CALL iotk_write_dat( iunpun, "ANGLE2", &
angle2(ityp)*180.d0/pi )
END DO
CALL iotk_write_end( iunpun, "STARTING_MAG" )
!
RETURN
!
END SUBROUTINE write_init_mag
!
!------------------------------------------------------------------------
SUBROUTINE write_xc( dft, nsp, lda_plus_u, & SUBROUTINE write_xc( dft, nsp, lda_plus_u, &
Hubbard_lmax, Hubbard_l, Hubbard_U, Hubbard_alpha ) Hubbard_lmax, Hubbard_l, Hubbard_U, Hubbard_alpha )
!------------------------------------------------------------------------ !------------------------------------------------------------------------

2
PW/Makefile
View File

@ -219,6 +219,8 @@ symtns.o \
symvect.o \ symvect.o \
symz.o \ symz.o \
tabd.o \ tabd.o \
transform_becsum_so.o \
transform_becsum_nc.o \
trntns.o \ trntns.o \
trnvecc.o \ trnvecc.o \
trnvect.o \ trnvect.o \

18
PW/input.f90
View File

@ -282,6 +282,18 @@ SUBROUTINE iosys()
tfixed_occ = .FALSE. tfixed_occ = .FALSE.
ltetra = .FALSE. ltetra = .FALSE.
! !
IF (noncolin) THEN
DO nt = 1, ntyp
!
angle1(nt) = pi * angle1(nt) / 180.D0
angle2(nt) = pi * angle2(nt) / 180.D0
!
END DO
ELSE
angle1=0.d0
angle2=0.d0
ENDIF
SELECT CASE( TRIM( occupations ) ) SELECT CASE( TRIM( occupations ) )
CASE( 'fixed' ) CASE( 'fixed' )
! !
@ -425,8 +437,8 @@ SUBROUTINE iosys()
! !
DO nt = 1, ntyp DO nt = 1, ntyp
! !
theta = pi * angle1(nt) / 180.D0 theta = angle1(nt)
phi = pi * angle2(nt) / 180.D0 phi = angle2(nt)
! !
mcons(1,nt) = starting_magnetization(nt) * SIN( theta ) * COS( phi ) mcons(1,nt) = starting_magnetization(nt) * SIN( theta ) * COS( phi )
mcons(2,nt) = starting_magnetization(nt) * SIN( theta ) * SIN( phi ) mcons(2,nt) = starting_magnetization(nt) * SIN( theta ) * SIN( phi )
@ -444,7 +456,7 @@ SUBROUTINE iosys()
! !
DO nt = 1, ntyp DO nt = 1, ntyp
! !
theta = pi * angle1(nt) / 180.D0 theta = angle1(nt)
! !
mcons(3,nt) = cos(theta) mcons(3,nt) = cos(theta)
! !

75
PW/pw_restart.f90
View File

@ -40,6 +40,7 @@ MODULE pw_restart
lpw_read = .FALSE., & lpw_read = .FALSE., &
lions_read = .FALSE., & lions_read = .FALSE., &
lspin_read = .FALSE., & lspin_read = .FALSE., &
lstarting_mag_read = .FALSE., &
lxc_read = .FALSE., & lxc_read = .FALSE., &
locc_read = .FALSE., & locc_read = .FALSE., &
lbz_read = .FALSE., & lbz_read = .FALSE., &
@ -85,7 +86,8 @@ MODULE pw_restart
USE spin_orb, ONLY : lspinorb, domag USE spin_orb, ONLY : lspinorb, domag
USE symme, ONLY : nsym, invsym, s, ftau, irt, t_rev USE symme, ONLY : nsym, invsym, s, ftau, irt, t_rev
USE char, ONLY : sname USE char, ONLY : sname
USE lsda_mod, ONLY : nspin, isk, lsda USE lsda_mod, ONLY : nspin, isk, lsda, starting_magnetization
USE noncollin_module, ONLY : angle1, angle2
USE ions_base, ONLY : amass USE ions_base, ONLY : amass
USE funct, ONLY : get_dft_name USE funct, ONLY : get_dft_name
USE scf, ONLY : rho USE scf, ONLY : rho
@ -332,6 +334,8 @@ MODULE pw_restart
! !
CALL write_spin( lsda, noncolin, npol, lspinorb, domag ) CALL write_spin( lsda, noncolin, npol, lspinorb, domag )
! !
CALL write_init_mag(starting_magnetization, angle1, angle2, nsp )
!
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! ... EXCHANGE_CORRELATION ! ... EXCHANGE_CORRELATION
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
@ -807,7 +811,7 @@ MODULE pw_restart
INTEGER, INTENT(OUT) :: ierr INTEGER, INTENT(OUT) :: ierr
! !
CHARACTER(LEN=256) :: dirname CHARACTER(LEN=256) :: dirname
LOGICAL :: lexist, lcell, lpw, lions, lspin, & LOGICAL :: lexist, lcell, lpw, lions, lspin, linit_mag, &
lxc, locc, lbz, lbs, lwfc, & lxc, locc, lbz, lbs, lwfc, &
lsymm, lph, lrho, lefield lsymm, lph, lrho, lefield
! !
@ -827,6 +831,7 @@ MODULE pw_restart
lpw = .FALSE. lpw = .FALSE.
lions = .FALSE. lions = .FALSE.
lspin = .FALSE. lspin = .FALSE.
linit_mag = .FALSE.
lxc = .FALSE. lxc = .FALSE.
locc = .FALSE. locc = .FALSE.
lbz = .FALSE. lbz = .FALSE.
@ -868,6 +873,7 @@ MODULE pw_restart
lpw = .TRUE. lpw = .TRUE.
lions = .TRUE. lions = .TRUE.
lspin = .TRUE. lspin = .TRUE.
linit_mag = .TRUE.
lxc = .TRUE. lxc = .TRUE.
locc = .TRUE. locc = .TRUE.
lbz = .TRUE. lbz = .TRUE.
@ -882,6 +888,7 @@ MODULE pw_restart
lpw = .TRUE. lpw = .TRUE.
lions = .TRUE. lions = .TRUE.
lspin = .TRUE. lspin = .TRUE.
linit_mag = .TRUE.
lxc = .TRUE. lxc = .TRUE.
locc = .TRUE. locc = .TRUE.
lbz = .TRUE. lbz = .TRUE.
@ -898,6 +905,7 @@ MODULE pw_restart
lpw_read = .FALSE. lpw_read = .FALSE.
lions_read = .FALSE. lions_read = .FALSE.
lspin_read = .FALSE. lspin_read = .FALSE.
lstarting_mag_read = .FALSE.
lxc_read = .FALSE. lxc_read = .FALSE.
locc_read = .FALSE. locc_read = .FALSE.
lbz_read = .FALSE. lbz_read = .FALSE.
@ -937,6 +945,12 @@ MODULE pw_restart
IF ( ierr > 0 ) RETURN IF ( ierr > 0 ) RETURN
! !
END IF END IF
IF (linit_mag) THEN
!
CALL read_init_mag( dirname, ierr )
IF ( ierr > 0 ) RETURN
!
ENDIF
IF ( lxc ) THEN IF ( lxc ) THEN
! !
CALL read_xc( dirname, ierr ) CALL read_xc( dirname, ierr )
@ -1772,6 +1786,63 @@ MODULE pw_restart
! !
END SUBROUTINE read_spin END SUBROUTINE read_spin
! !
!--------------------------------------------------------------------------
SUBROUTINE read_init_mag( dirname, ierr )
!------------------------------------------------------------------------
!
USE constants, ONLY : pi
USE lsda_mod, ONLY : starting_magnetization
USE noncollin_module, ONLY : angle1, angle2
!
IMPLICIT NONE
!
CHARACTER(LEN=*), INTENT(IN) :: dirname
INTEGER, INTENT(OUT) :: ierr
!
LOGICAL :: found
INTEGER :: ityp, ntyp
!
IF ( lstarting_mag_read ) RETURN
!
IF ( ionode ) &
CALL iotk_open_read( iunpun, FILE = TRIM( dirname ) // '/' // &
& TRIM( xmlpun ), BINARY = .FALSE., IERR = ierr )
!
CALL mp_bcast( ierr, ionode_id, intra_image_comm )
!
IF ( ierr > 0 ) RETURN
!
IF ( ionode ) THEN
!
CALL iotk_scan_begin( iunpun, "STARTING_MAG" )
!
CALL iotk_scan_dat( iunpun, "NTYP", ntyp )
!
DO ityp=1,ntyp
CALL iotk_scan_dat( iunpun, "STARTING_MAGNETIZATION", &
starting_magnetization(ityp) )
CALL iotk_scan_dat( iunpun, "ANGLE1", angle1(ityp) )
CALL iotk_scan_dat( iunpun, "ANGLE2", angle2(ityp) )
angle1(ityp)=angle1(ityp)*pi/180.d0
angle2(ityp)=angle2(ityp)*pi/180.d0
END DO
CALL iotk_scan_end( iunpun, "STARTING_MAG" )
!
CALL iotk_close_read( iunpun )
!
END IF
!
CALL mp_bcast( starting_magnetization, ionode_id, intra_image_comm )
CALL mp_bcast( angle1, ionode_id, intra_image_comm )
CALL mp_bcast( angle2, ionode_id, intra_image_comm )
!
lstarting_mag_read = .TRUE.
!
RETURN
!
END SUBROUTINE read_init_mag
!
!------------------------------------------------------------------------ !------------------------------------------------------------------------
SUBROUTINE read_xc( dirname, ierr ) SUBROUTINE read_xc( dirname, ierr )
!------------------------------------------------------------------------ !------------------------------------------------------------------------

10
PW/setup.f90
View File

@ -141,12 +141,12 @@ SUBROUTINE setup()
! !
! ... transform angles to radiants ! ... transform angles to radiants
! !
DO nt = 1, ntyp ! DO nt = 1, ntyp
! !
angle1(nt) = pi * angle1(nt) / 180.D0 ! angle1(nt) = pi * angle1(nt) / 180.D0
angle2(nt) = pi * angle2(nt) / 180.D0 ! angle2(nt) = pi * angle2(nt) / 180.D0
! ! !
END DO ! END DO
! !
! ... Set the nomag variable to make a spin-orbit calculation with zero ! ... Set the nomag variable to make a spin-orbit calculation with zero
! ... magnetization ! ... magnetization

4
PW/sgama.f90
View File

@ -118,8 +118,8 @@ subroutine sgama (nrot, nat, s, sname, t_rev, at, bg, tau, ityp, nsym,&
call smallg_q (xq, modenum, at, bg, nrot, s, ftau, sym, minus_q) call smallg_q (xq, modenum, at, bg, nrot, s, ftau, sym, minus_q)
IF (noncolin.and.domag) THEN IF (noncolin.and.domag) THEN
minus_q=.false. minus_q=.false.
IF ( ABS(DOT_PRODUCT(xq,xq)) > 1.0D-07 ) CALL errore ('sgama', & ! IF ( ABS(DOT_PRODUCT(xq,xq)) > 1.0D-07 ) CALL errore ('sgama', &
'phonon not implemented with non collinear magnetism', 1) ! 'phonon not implemented with non collinear magnetism', 1)
! If somebody want to implement phononic calculation in non ! If somebody want to implement phononic calculation in non
! collinear magnetic case he has to pay attention to the fact ! collinear magnetic case he has to pay attention to the fact
! that in non collinear case the symmetry k -> -k is not ! that in non collinear case the symmetry k -> -k is not

161
PW/sum_band.f90
View File

@ -371,15 +371,19 @@ SUBROUTINE sum_band()
COMPLEX(DP), ALLOCATABLE :: becp(:,:), becp_nc(:,:,:) COMPLEX(DP), ALLOCATABLE :: becp(:,:), becp_nc(:,:,:)
! contains <beta|psi> ! contains <beta|psi>
! !
COMPLEX(DP), ALLOCATABLE :: be1(:,:), be2(:,:) COMPLEX(DP), ALLOCATABLE :: becsum_nc(:,:,:,:)
! !
INTEGER :: ipol, kh, kkb, is1, is2 INTEGER :: ipol, kh, kkb, is1, is2, js
! !
IF (okvan) THEN
IF (noncolin) THEN IF (noncolin) THEN
ALLOCATE(becsum_nc(nhm*(nhm+1)/2,nat,npol,npol))
becsum_nc=(0.d0, 0.d0)
ALLOCATE( becp_nc( nkb, npol, nbnd ) ) ALLOCATE( becp_nc( nkb, npol, nbnd ) )
IF (lspinorb) ALLOCATE(be1(nhm,2), be2(nhm,2))
ELSE ELSE
ALLOCATE( becp( nkb, nbnd ) ) ALLOCATE( becp( nkb, nbnd ) )
END IF
ENDIF ENDIF
! !
! ... here we sum for each k point the contribution ! ... here we sum for each k point the contribution
@ -504,30 +508,6 @@ SUBROUTINE sum_band()
! !
IF (ityp(na)==np) THEN IF (ityp(na)==np) THEN
! !
IF (so(np)) THEN
be1=(0.d0,0.d0)
be2=(0.d0,0.d0)
DO ih = 1, nh(np)
ikb = ijkb0 + ih
DO kh = 1, nh(np)
IF ((nhtol(kh,np)==nhtol(ih,np)).AND. &
(nhtoj(kh,np)==nhtoj(ih,np)).AND. &
(indv(kh,np)==indv(ih,np))) THEN
kkb=ijkb0 + kh
DO is1=1,2
DO is2=1,2
be1(ih,is1)=be1(ih,is1)+ &
fcoef(ih,kh,is1,is2,np)* &
becp_nc(kkb,is2,ibnd)
be2(ih,is1)=be2(ih,is1)+ &
fcoef(kh,ih,is2,is1,np)* &
CONJG(becp_nc(kkb,is2,ibnd))
END DO
END DO
END IF
END DO
END DO
END IF
ijh = 1 ijh = 1
! !
DO ih = 1, nh(np) DO ih = 1, nh(np)
@ -536,46 +516,24 @@ SUBROUTINE sum_band()
! !
IF (noncolin) THEN IF (noncolin) THEN
! !
IF (so(np)) THEN DO is=1,npol
becsum(ijh,na,1)=becsum(ijh,na,1)+ w1*& !
(be1(ih,1)*be2(ih,1)+ be1(ih,2)*be2(ih,2)) DO js=1,npol
IF (domag) THEN becsum_nc(ijh,na,is,js) = &
becsum(ijh,na,2)=becsum(ijh,na,2)+ w1*& becsum_nc(ijh,na,is,js)+w1 * &
(be1(ih,2)*be2(ih,1)+ be1(ih,1)*be2(ih,2)) CONJG(becp_nc(ikb,is,ibnd)) * &
becsum(ijh,na,3)=becsum(ijh,na,3)+ & becp_nc(ikb,js,ibnd)
w1*(0.d0,-1.d0)* & END DO
(be1(ih,2)*be2(ih,1)-be1(ih,1)*be2(ih,2)) !
becsum(ijh,na,4)=becsum(ijh,na,4)+ w1* & END DO
(be1(ih,1)*be2(ih,1)-be1(ih,2)*be2(ih,2)) !
ENDIF
ELSE
becsum(ijh,na,1) = becsum(ijh,na,1) &
+ w1*( CONJG(becp_nc(ikb,1,ibnd)) &
*becp_nc(ikb,1,ibnd) &
+ CONJG(becp_nc(ikb,2,ibnd)) &
*becp_nc(ikb,2,ibnd) )
IF (domag) THEN
becsum(ijh,na,2)=becsum(ijh,na,2) &
+ w1*(CONJG(becp_nc(ikb,2,ibnd)) &
*becp_nc(ikb,1,ibnd) &
+ CONJG(becp_nc(ikb,1,ibnd)) &
*becp_nc(ikb,2,ibnd) )
becsum(ijh,na,3)=becsum(ijh,na,3) &
+ w1*2.d0 &
*AIMAG(CONJG(becp_nc(ikb,1,ibnd))* &
becp_nc(ikb,2,ibnd) )
becsum(ijh,na,4) = becsum(ijh,na,4) &
+ w1*( CONJG(becp_nc(ikb,1,ibnd)) &
*becp_nc(ikb,1,ibnd) &
- CONJG(becp_nc(ikb,2,ibnd)) &
*becp_nc(ikb,2,ibnd) )
END IF
END IF
ELSE ELSE
!
becsum(ijh,na,current_spin) = & becsum(ijh,na,current_spin) = &
becsum(ijh,na,current_spin) + & becsum(ijh,na,current_spin) + &
w1 * DBLE( CONJG( becp(ikb,ibnd) ) * & w1 * DBLE( CONJG( becp(ikb,ibnd) ) * &
becp(ikb,ibnd) ) becp(ikb,ibnd) )
!
END IF END IF
! !
ijh = ijh + 1 ijh = ijh + 1
@ -585,53 +543,18 @@ SUBROUTINE sum_band()
jkb = ijkb0 + jh jkb = ijkb0 + jh
! !
IF (noncolin) THEN IF (noncolin) THEN
IF (so(np)) THEN !
becsum(ijh,na,1)=becsum(ijh,na,1)+ w1*( & DO is=1,npol
(be1(jh,1)*be2(ih,1)+be1(jh,2)*be2(ih,2))+ & !
(be1(ih,1)*be2(jh,1)+be1(ih,2)*be2(jh,2))) DO js=1,npol
IF (domag) THEN becsum_nc(ijh,na,is,js) = &
becsum(ijh,na,2)=becsum(ijh,na,2)+w1*( & becsum_nc(ijh,na,is,js) + w1 * &
(be1(jh,2)*be2(ih,1)+be1(jh,1)*be2(ih,2))+& CONJG(becp_nc(ikb,is,ibnd)) * &
(be1(ih,2)*be2(jh,1)+be1(ih,1)*be2(jh,2))) becp_nc(jkb,js,ibnd)
becsum(ijh,na,3)=becsum(ijh,na,3)+ & END DO
w1*(0.d0,-1.d0)*((be1(jh,2)*& !
be2(ih,1)-be1(jh,1)*be2(ih,2))+ & END DO
(be1(ih,2)*be2(jh,1)-be1(ih,1)*& !
be2(jh,2)) )
becsum(ijh,na,4)=becsum(ijh,na,4)+ &
w1*((be1(jh,1)*be2(ih,1)- &
be1(jh,2)*be2(ih,2))+ &
(be1(ih,1)*be2(jh,1)- &
be1(ih,2)*be2(jh,2)) )
END IF
ELSE
becsum(ijh,na,1)= becsum(ijh,na,1)+ &
w1*2.d0* &
DBLE(CONJG(becp_nc(ikb,1,ibnd))* &
becp_nc(jkb,1,ibnd) + &
CONJG(becp_nc(ikb,2,ibnd))* &
becp_nc(jkb,2,ibnd) )
IF (domag) THEN
becsum(ijh,na,2)=becsum(ijh,na,2)+ &
w1*2.d0* &
DBLE(CONJG(becp_nc(ikb,2,ibnd))* &
becp_nc(jkb,1,ibnd) + &
CONJG(becp_nc(ikb,1,ibnd))* &
becp_nc(jkb,2,ibnd) )
becsum(ijh,na,3)=becsum(ijh,na,3)+ &
w1*2.d0* &
AIMAG(CONJG(becp_nc(ikb,1,ibnd))* &
becp_nc(jkb,2,ibnd) + &
CONJG(becp_nc(ikb,1,ibnd))* &
becp_nc(jkb,2,ibnd) )
becsum(ijh,na,4)=becsum(ijh,na,4)+ &
w1*2.d0* &
DBLE(CONJG(becp_nc(ikb,1,ibnd))* &
becp_nc(jkb,1,ibnd) - &
CONJG(becp_nc(ikb,2,ibnd))* &
becp_nc(jkb,2,ibnd) )
END IF
END IF
ELSE ELSE
! !
becsum(ijh,na,current_spin) = & becsum(ijh,na,current_spin) = &
@ -669,13 +592,31 @@ SUBROUTINE sum_band()
CALL stop_clock( 'becsum' ) CALL stop_clock( 'becsum' )
! !
END DO k_loop END DO k_loop
IF (noncolin.and.okvan) THEN
DO np = 1, ntyp
IF ( tvanp(np) ) THEN
DO na = 1, nat
IF (ityp(na)==np) THEN
IF (so(np)) THEN
CALL transform_becsum_so(becsum_nc,becsum,na)
ELSE
CALL transform_becsum_nc(becsum_nc,becsum,na)
END IF
END IF
END DO
END IF
END DO
END IF
! !
IF (okvan) THEN
IF (noncolin) THEN IF (noncolin) THEN
DEALLOCATE( becsum_nc )
DEALLOCATE( becp_nc ) DEALLOCATE( becp_nc )
IF (lspinorb) DEALLOCATE(be1, be2)
ELSE ELSE
DEALLOCATE( becp ) DEALLOCATE( becp )
ENDIF ENDIF
END IF
! !
RETURN RETURN
! !

69
PW/transform_becsum_nc.f90 Normal file
View File

@ -0,0 +1,69 @@
!
! Copyright (C) 2006 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 .
!
#include "f_defs.h"
!
!----------------------------------------------------------------------------
SUBROUTINE transform_becsum_nc(becsum_nc,becsum,na)
!----------------------------------------------------------------------------
!
! This routine multiply becsum_nc by the identity and the Pauli
! matrices and saves it in becsum for the calculation of
! augmentation charge and magnetization.
!
USE kinds, ONLY : DP
USE ions_base, ONLY : nat, ntyp => nsp, ityp
USE uspp_param, ONLY : nh, nhm
USE lsda_mod, ONLY : nspin
USE noncollin_module, ONLY : npol
USE spin_orb, ONLY : domag
!
IMPLICIT NONE
COMPLEX(DP) :: becsum_nc(nhm*(nhm+1)/2,nat,npol,npol)
REAL(DP) :: becsum(nhm*(nhm+1)/2,nat,nspin)
INTEGER :: na
!
! ... local variables
!
INTEGER :: ih, jh, ijh, np
np=ityp(na)
ijh=1
DO ih = 1, nh(np)
becsum(ijh,na,1)= becsum(ijh,na,1)+ &
becsum_nc(ijh,na,1,1)+becsum_nc(ijh,na,2,2)
IF (domag) THEN
becsum(ijh,na,2)= becsum(ijh,na,2)+ &
becsum_nc(ijh,na,1,2)+becsum_nc(ijh,na,2,1)
becsum(ijh,na,3)= becsum(ijh,na,3)+(0.d0,-1.d0)* &
(becsum_nc(ijh,na,1,2)-becsum_nc(ijh,na,2,1))
becsum(ijh,na,4)= becsum(ijh,na,4)+ &
becsum_nc(ijh,na,1,1)-becsum_nc(ijh,na,2,2)
END IF
ijh=ijh+1
DO jh = ih+1, nh(np)
becsum(ijh,na,1)= becsum(ijh,na,1) + &
(becsum_nc(ijh,na,1,1)+becsum_nc(ijh,na,2,2)) &
+ CONJG(becsum_nc(ijh,na,1,1)+becsum_nc(ijh,na,2,2))
IF (domag) THEN
becsum(ijh,na,2)= becsum(ijh,na,2) + &
becsum_nc(ijh,na,1,2)+becsum_nc(ijh,na,2,1) &
+ CONJG(becsum_nc(ijh,na,2,1)+becsum_nc(ijh,na,1,2))
becsum(ijh,na,3)= becsum(ijh,na,3) +(0.d0,-1.d0)* &
(becsum_nc(ijh,na,1,2)-becsum_nc(ijh,na,2,1) &
+ CONJG(becsum_nc(ijh,na,2,1)-becsum_nc(ijh,na,1,2)) )
becsum(ijh,na,4)= becsum(ijh,na,4) + &
(becsum_nc(ijh,na,1,1)-becsum_nc(ijh,na,2,2)) &
+ CONJG(becsum_nc(ijh,na,1,1)-becsum_nc(ijh,na,2,2))
END IF
ijh=ijh+1
END DO
END DO
RETURN
END SUBROUTINE transform_becsum_nc

129
PW/transform_becsum_so.f90 Normal file
View File

@ -0,0 +1,129 @@
!
! Copyright (C) 2006 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 .
!
#include "f_defs.h"
!
!----------------------------------------------------------------------------
SUBROUTINE transform_becsum_so(becsum_nc,becsum,na)
!----------------------------------------------------------------------------
!
! This routine multiply becsum_nc by the identity and the Pauli
! matrices, rotate it as appropriate for the spin-orbit case
! and saves it in becsum for the calculation of
! augmentation charge and magnetization.
!
USE kinds, ONLY : DP
USE ions_base, ONLY : nat, ntyp => nsp, ityp
USE uspp_param, ONLY : nh, nhm
USE lsda_mod, ONLY : nspin
USE noncollin_module, ONLY : npol
USE spin_orb, ONLY : fcoef, domag
!
IMPLICIT NONE
COMPLEX(DP) :: becsum_nc(nhm*(nhm+1)/2,nat,npol,npol)
REAL(DP) :: becsum(nhm*(nhm+1)/2,nat,nspin)
INTEGER :: na
!
! ... local variables
!
INTEGER :: ih, jh, lh, kh, ijh, np, is1, is2
INTEGER, ALLOCATABLE :: ijh_save(:,:)
COMPLEX(DP) :: fac
INTEGER :: find_ijh, ijh_l
LOGICAL :: same_lj
ALLOCATE(ijh_save(nhm,nhm))
np=ityp(na)
DO ih=1,nh(np)
DO jh=1,nh(np)
ijh_save(ih,jh)=find_ijh(ih,jh,nh(np))
END DO
END DO
DO ih = 1, nh(np)
DO jh = 1, nh(np)
ijh=ijh_save(ih,jh)
DO kh = 1, nh(np)
IF (same_lj(kh,ih,np)) THEN
DO lh=1,nh(np)
IF (same_lj(lh,jh,np)) THEN
ijh_l=ijh_save(kh,lh)
DO is1=1,npol
DO is2=1,npol
IF (kh <= lh) THEN
fac=becsum_nc(ijh_l,na,is1,is2)
ELSE
fac=CONJG(becsum_nc(ijh_l,na,is2,is1))
ENDIF
becsum(ijh,na,1)=becsum(ijh,na,1) + fac * &
(fcoef(kh,ih,is1,1,np)*fcoef(jh,lh,1,is2,np) + &
fcoef(kh,ih,is1,2,np)*fcoef(jh,lh,2,is2,np) )
IF (domag) THEN
becsum(ijh,na,2)=becsum(ijh,na,2)+fac * &
(fcoef(kh,ih,is1,1,np)*fcoef(jh,lh,2,is2,np) +&
fcoef(kh,ih,is1,2,np)*fcoef(jh,lh,1,is2,np) )
becsum(ijh,na,3)=becsum(ijh,na,3)+fac*(0.d0,-1.d0)*&
(fcoef(kh,ih,is1,1,np)*fcoef(jh,lh,2,is2,np) - &
fcoef(kh,ih,is1,2,np)*fcoef(jh,lh,1,is2,np) )
becsum(ijh,na,4)=becsum(ijh,na,4) + fac * &
(fcoef(kh,ih,is1,1,np)*fcoef(jh,lh,1,is2,np) - &
fcoef(kh,ih,is1,2,np)*fcoef(jh,lh,2,is2,np) )
END IF
END DO
END DO
END IF
END DO
END IF
END DO
END DO
END DO
!
DEALLOCATE(ijh_save)
RETURN
END SUBROUTINE transform_becsum_so
FUNCTION same_lj(ih,jh,np)
USE uspp, ONLY : nhtol, nhtoj, indv
IMPLICIT NONE
LOGICAL :: same_lj
INTEGER :: ih, jh, np
same_lj = ((nhtol(ih,np)==nhtol(jh,np)).AND. &
(ABS(nhtoj(ih,np)-nhtoj(jh,np))<1.d8).AND. &
(indv(ih,np)==indv(jh,np)) )
RETURN
END FUNCTION same_lj
FUNCTION find_ijh(ih,jh,nh)
IMPLICIT NONE
INTEGER :: find_ijh, ih, jh, nh
INTEGER :: ih0, jh0, ijh, i, j
if (jh > ih) then
ih0=ih
jh0=jh
else
ih0=jh
jh0=ih
end if
ijh=0
do i=1, ih0-1
do j=i, nh
ijh=ijh+1
enddo
enddo
do j=ih0, jh0
ijh=ijh+1
enddo
find_ijh=ijh
end function find_ijh