[求助] 代码错误Incompatible ranks 0 and 2 in assignment

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