Fortran Coder

查看: 206|回复: 4

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

[复制链接]

6

帖子

2

主题

0

精华

入门

F 币
45 元
贡献
23 点
发表于 2019-5-12 17:19:39 | 显示全部楼层 |阅读模式
大家好,
我写了段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)

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

回复

使用道具 举报

1418

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
698 元
贡献
530 点

新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2019-5-12 20:18:53 | 显示全部楼层
matmul 的结果是二维数组,而 ta4 是单变量。无法存储matmul的结果

6

帖子

2

主题

0

精华

入门

F 币
45 元
贡献
23 点
 楼主| 发表于 2019-5-12 20:54:28 | 显示全部楼层
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))。这样一来,我最终应该得到是一个数才对呀。

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

1418

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
698 元
贡献
530 点

新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2019-5-12 21:37:59 | 显示全部楼层
是 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的二维数组求和,就能得到它里面唯一的那个元素。

6

帖子

2

主题

0

精华

入门

F 币
45 元
贡献
23 点
 楼主| 发表于 2019-5-13 12:58:48 | 显示全部楼层
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的矩阵存储起来。

非常感谢你的帮助。
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2019-8-21 12:34

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表