代码错误Incompatible ranks 0 and 2 in assignment

Kieran

大家好，
我写了段FORTRAN代码，编译的时候，总是出错，我没弄清楚哪里出了问题。所以我把我代码中有问题的部分粘贴到这里了，麻烦大家帮我看下，是哪里写得不对。谢谢啦。

我的有问题的代码部分如下。

[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

编译后，编译器总是提示下面的错误信息。
CALCULATION.f90:298:0:
ta4 = REAL(MATMUL(MATMUL(ta2,vel(1+(i-1)*num_wan:i*num_wan,:)),ta1))
1. Error: Incompatible ranks 0 and 2 in assignment at (1)

这个错误信息好像是说数组的维度不对等，但我仔细查看了下，也没看出为什么会不对等。能否麻烦大家帮我看下，我的代码哪里写得有问题呢？谢谢大家啦。

Kieran

fcode 发表于 2019-5-12 21:37
是 1*1 的二维数组。就好像 real :: a(1,1) 一样
你可以用
ta4 = sum(REAL(MATMUL(MATMUL(ta2,vel(1+(i-1)* ...

好的，我现在懂了。就是说MATMUL这个内建函数会把结果以矩阵的形式存储起来，即便最终得到的是一个数，也仍然会以1*1的矩阵存储起来。

非常感谢你的帮助。

fcode

是 1*1 的二维数组。就好像 real :: a(1,1) 一样
你可以用
ta4 = sum(REAL(MATMUL(MATMUL(ta2,vel(1+(i-1)*num_wan:i*num_wan,:)),ta1)))
得到单变量。
这个 sum 是假装的。对一个1*1的二维数组求和，就能得到它里面唯一的那个元素。

Kieran

fcode 发表于 2019-5-12 20:18
matmul 的结果是二维数组，而 ta4 是单变量。无法存储matmul的结果

谢谢你的回复。

不过我没理解，为什么会是二维数组。我的代码里vel是一个46*46的矩阵，ta1是一个46*1的列向量，ta2是一个1*46的矩阵。我首先是进行了ta2和vel的矩阵相乘（MATMUL(ta2,vel)）得到的应该是一个1*46的矩阵。再然后我又把这个结果和ta1矩阵相乘（MATMUL(MATMUL(ta2,vel),ta1)）。这样一来，我最终应该得到是一个数才对呀。

为什么最后是一个二维数组呢？还望不吝赐教。

fcode

matmul 的结果是二维数组，而 ta4 是单变量。无法存储matmul的结果