[Fortran] 纯文本查看 复制代码
SUBROUTINE SPIN_ORBIT_TORQUE(efermi,rea_cons,tb,evc,evl,num_wan,monk,kp,sot_total)
!Parameters used in the main code
DOUBLE PRECISION :: efermi !Fermi energy level
DOUBLE PRECISION :: rea_cons(3,3) !Lattice constant of unit cell in real space
DOUBLE PRECISION :: rec_cons(3,3) !Vectors of unit cell in the reciprocal space
DOUBLE COMPLEX, ALLOCATABLE :: tb(:,:) !Array to store the extracted tight binding Hamiltonian matrix
DOUBLE COMPLEX, ALLOCATABLE :: evc(:,:) !Array to store the Eigen vector matrix
DOUBLE PRECISION, ALLOCATABLE :: evl(:,:) !Array to store the Eigen value
INTEGER :: num_wan !Number of Wannier function
INTEGER :: monk(2) !Monkhorst k point mesh
DOUBLE PRECISION, ALLOCATABLE :: kp(:,:) !Array to store the coodinates of k point
DOUBLE PRECISION :: sot_total !Total Spin Orbit Torque of the whole system
INTEGER :: i, j, k, l
DOUBLE PRECISION :: boltzm !Boltzmann constant
DOUBLE PRECISION :: temper !Temperature in Fermi Dirac Distribution
DOUBLE PRECISION :: h_bar !Planck constant over pi
DOUBLE PRECISION :: ele !Electron quantity
DOUBLE PRECISION :: bro !Broadening value from defects
DOUBLE PRECISION :: vol !Volume or Area of the unit cell
DOUBLE PRECISION, ALLOCATABLE :: exc(:) !Exchange potential
DOUBLE PRECISION :: lec(3) !Length of the three boundaries in the triangle within the
!in-plane lattice vector in real space
DOUBLE PRECISION :: dc !Transit value
DOUBLE COMPLEX, ALLOCATABLE :: vel(:,:) !Velocity operator of Hamiltonian matrix
DOUBLE COMPLEX :: pc(2) !Complex number 1.0d0 and i used in this subroutine
DOUBLE PRECISION :: sx, sy, sz !Expectation value of spin operator
DOUBLE COMPLEX, ALLOCATABLE :: ta1(:,:), ta2(:,:) !Transit vector for spin density calculation
DOUBLE COMPLEX, ALLOCATABLE :: ta3(:,:) !Transit vector for spin density calculation
DOUBLE PRECISION :: ta4 !Transit value for the expectation value of velocity operator
DOUBLE PRECISION, ALLOCATABLE :: intr(:,:) !Spin density from intra band part
DOUBLE PRECISION, ALLOCATABLE :: fd(:) !Fermi Dirac distribution for each eigen state
DOUBLE PRECISION, ALLOCATABLE :: inte(:,:) !Spin density from inter band part
DOUBLE PRECISION, ALLOCATABLE :: sd_tot(:,:) !Total Spin density
DOUBLE PRECISION, ALLOCATABLE :: sot(:) !Spin Orbit Torque on each k point
!Complex constant in this subroutine
pc(1) = (1.0d0, 0.0d0)
pc(2) = (0.0d0, 1.0d0)
!
!Parameterizing the variables
boltzm = 8.6173303d-5 !Boltzmann constant (unit: eV/K)
temper = 0.0d0 !Temperature in Fermi Dirac Distribution (unit: K)
h_bar = 6.582119514d-16 !Planck constant over pi (unit: eV*s)
bro = 1.5d-2 !Broadening from defects
ele = 1.6021766208d-19 !Elementary charge (unit: Coulombs)
!
!Computing Volume or Area of the unit cell based on Heron's Formula
lec(1) = DSQRT((rea_cons(1,1) - 0.0d0) ** 2 &
+ (rea_cons(1,2) - 0.0d0) ** 2 + (rea_cons(1,3) - 0.0d0) ** 2)
lec(2) = DSQRT((rea_cons(2,1) - 0.0d0) ** 2 &
+ (rea_cons(2,2) - 0.0d0) ** 2 + (rea_cons(2,3) - 0.0d0) ** 2)
lec(3) = DSQRT((rea_cons(1,1) - rea_cons(2,1)) ** 2 &
+ (rea_cons(1,2) - rea_cons(2,2)) ** 2 + (rea_cons(1,3) - rea_cons(2,3)) ** 2)
dc = (lec(1) + lec(2) + lec(3)) / 2
vol = DSQRT(dc * (dc - lec(1)) * (dc - lec(2)) * (dc - lec(3)))
!
!Allocating the array to store the exchange potential on each k point
ALLOCATE (exc(monk(1)))
!
!Allocating the transit arrays for the spin density calculation
ALLOCATE (ta1(num_wan,1))
ALLOCATE (ta2(1,num_wan))
ALLOCATE (ta3(1,num_wan))
!
!Calculating the total Exchange potential on each k point
exc = 0.0d0
DO i = 1, monk(1), 1
DO j = 1, num_wan / 2, 1
exc(i) = exc(i) + evl(i,j) - evl(i,j+num_wan/2)
END DO
END DO
!
!Allocating the array to store the velocity operator of Hamiltonian matrix
ALLOCATE (vel(monk(1)*num_wan,num_wan))
!
!Computing velocity operator of Hamiltonian matrix
DO i = 1, monk(1), 1
IF (i .EQ. 1) THEN
dc = kp(2,4) * rec_cons(1,1) + kp(2,5) * rec_cons(2,1) + kp(2,6) * rec_cons(3,1)&
- (kp(monk(1)-1,4) * (-1.0d0) * rec_cons(1,1) + kp(monk(1)-1,5) * rec_cons(2,1)&
+ kp(monk(1)-1,6) * rec_cons(3,1))
vel(1:num_wan,:) = evc(1+num_wan:2*num_wan,:) - evc(1+(monk(1)-2)*num_wan:(monk(1)-1)*num_wan,:)&
/ dc
ELSE IF (i .EQ. monk(1)) THEN
dc = kp(2,4) * 2.0d0 * rec_cons(1,1) + kp(2,5) * rec_cons(2,1) + kp(2,6) * rec_cons(3,1)&
- (kp(monk(1)-1,4) * rec_cons(1,1) + kp(monk(1)-1,5) * rec_cons(2,1)&
+ kp(monk(1)-1,6) *rec_cons(3,1))
vel(1+(monk(1)-1)*num_wan:monk(1)*num_wan,:) = evc(1+num_wan:2*num_wan,:) - evc(1+(monk(1)-2)&
*num_wan:(monk(1)-1)*num_wan,:) / dc
ELSE
dc = kp(i+1,1) - kp(i-1,1)
vel(1+(i-1)*num_wan:i*num_wan,:) = evc(1+i*num_wan:(i+1)*num_wan,:)&
- evc(1+(i-2)*num_wan:(i-1)*num_wan,:) / dc
END IF
END DO
!
!Allocating the array to store the total spin density
ALLOCATE (sd_tot(monk(1),3))
!
!Allocating the arrays to store the spin density from intra band part
ALLOCATE (intr(monk(1)*num_wan,3))
!
!Calculating the 'Intra-band' part of the spin density
intr = 0.0d0
sd_tot = 0.0d0
DO i = 1, monk(1), 1
DO j = 1, num_wan, 1
DO k = 1, num_wan / 2, 1
sx = DCONJG(evc(k+(i-1)*num_wan+num_wan/2,j)) * evc(k+(i-1)*num_wan,j)&
+ DCONJG(evc(k+(i-1)*num_wan,j)) * evc(k+(i-1)*num_wan+num_wan/2,j)
sy = REAL(pc(2) * DCONJG(evc(k+(i-1)*num_wan+num_wan/2,j)) * evc(k+(i-1)*num_wan,j)&
- pc(2) * DCONJG(evc(k+(i-1)*num_wan,j)) * evc(k+(i-1)*num_wan+num_wan/2,j))
sz = DCONJG(evc(k+(i-1)*num_wan,j)) * evc(k+(i-1)*num_wan,j)&
- DCONJG(evc(k+(i-1)*num_wan+num_wan/2,j)) * evc(k+(i-1)*num_wan+num_wan/2,j)
END DO
ta1(:,1) = evc(1+(i-1)*num_wan:i*num_wan,j)
ta2 = DCONJG(TRANSPOSE(ta1))
ta4 = REAL(MATMUL(MATMUL(ta2,vel(1+(i-1)*num_wan:i*num_wan,:)),ta1)) ！就是这里出了问题，但我没想清楚为什么不对。
intr(j+(i-1)*num_wan,1) = intr(j+(i-1)*num_wan,1) + sx * ta4&
* 1.0d0 / ((evl(i,j) - efermi) **2 + bro **2)
intr(j+(i-1)*num_wan,2) = intr(j+(i-1)*num_wan,2) + sy * ta4&
* 1.0d0 / ((evl(i,j) - efermi) **2 + bro **2)
intr(j+(i-1)*num_wan,3) = intr(j+(i-1)*num_wan,3) + sz * ta4&
* 1.0d0 / ((evl(i,j) - efermi) **2 + bro **2)
END DO
intr(1+(i-1)*num_wan:i*num_wan,:) = intr(1+(i-1)*num_wan:i*num_wan,:) * ele&
* h_bar / 2.0d0 / bro / vol
DO j = 1, num_wan, 1
sd_tot(i,1) = sd_tot(i,1) + intr(j+(i-1)*num_wan,1)
sd_tot(i,2) = sd_tot(i,2) + intr(j+(i-1)*num_wan,2)
sd_tot(i,3) = sd_tot(i,3) + intr(j+(i-1)*num_wan,3)
END DO
END DO