看起来这是代码改动带来的“口径变化”，不是画图出错。

你现在的图（中心 Γ 附近出现很小的一圈高值，半径随填充变小）正是论文 Fig.2 所描述的“低填充时、以 Γ 为中心的铁磁样自旋涨落主导、在一个小圆环上取最大”的形态。你之所以从原先的图变成现在这样，主要有两点代码层面的变化：

1. 你不再对完整的 4×4 非相互作用易感矩阵 $\chi_0(q)$ 求本征值，而是先取了它的自旋通道 2×2 子块再做特征分解。
   在 RPA5.for 里你把

• （旧）对 chi_0(4×4) 求本征值并取 maxval(w0) 的步骤注释掉了，
• （新）改成构造
  $\texttt{mag0}=\begin{pmatrix}\chi_0(1,1)&\chi_0(1,4)\\ \chi_0(4,1)&\chi_0(4,4)\end{pmatrix}$
  然后对这个 2×2 矩阵求本征值、取最大值 w0_max 再画图。
  这等价于投到“自旋（intra-sublattice）通道”，因此图形会从之前“混合了电荷/自旋信息的 4×4 结果”变成“在 q≈0 处强峰的小圆环”，这与论文结论一致（低填充时半径随填充缩小）。

2. 你在写文件时对结果取了负号：
   sus0.dat 是 write(..., -w0_max) 输出的。
   如果你的画图脚本按“$\lambda_{\max}(\chi^{(0)}_s)$”当作正值去画，这个负号会让色标范围和明暗反过来（但形状不变）。想与论文配色直观一致，可去掉这个负号，或在画图时再乘以 -1。

另外两点小提示（不影响形状，但会影响你看到的“圈的大小/亮度”）：

• 圈的半径随填充（化学势 μ）改变而变：填充越低，圈越小；你代码里当前设置 u=0（非相互作用）和某个 μ 值，对应的就是论文里“低填充”的情况。
• 你把 $(q_x,q_y)$ 映射到了六角晶格坐标 $(q_x - q_y/2,\ \sqrt{3}q_y/2)$，所以图边界是个正六边形——这是对蜂窝晶格第一布里渊区的正确展示，不是异常。

想恢复到“改动前”的那种 4×4 总易感本征值的图，只需把对 mag0(2×2) 的特征分解改回对 chi_0(4×4) 的特征分解（并去掉输出时的负号）；若只想看单个分量（例如 $\chi^{(0)}_{11}$），则直接把 -dreal(chi_0(1,1))（注意号）写出来画即可。

总之：你现在的结果其实更“物理”——它专门展示了自旋通道的最大本征值，因此在 Γ 附近出现小圆环峰，这是对的。如果你需要，我可以把 4×4 版本与 2×2 自旋通道版本并排画出来，方便直观比较差异。

好的，我把你上传的 RPA5.for 和你截图里的 MATLAB 脚本逐行看了一遍。你现在的图和文献 Fig.2 不一致，核心有四个原因（前两个是决定性的）：

---

结论先说

1. Lindhard 泡泡少了整体负号
   在 RPA5.for 里构造 $\chi^{(0)}$ 的地方使用了

$$\chi^{(0)}\sim \sum_k \frac{f(E_k)-f(E_{k+q})}{E_k-E_{k+q}+i0^+}$$

但标准定义前面有一个 负号：$\chi^{(0)}(q,0)=-\sum_k\frac{f(E_k)-f(E_{k+q})}{E_k-E_{k+q}}$。
你现在没有这个负号，导致非相互作用自旋易感性整体变负，所以 MATLAB 图的色条全是负值（见你的 colorbar：-0.022 到 -0.001），和论文正值配色相反。

2. 你输出的是 2×2 “自旋通道子块”的本征值，而 Fig.2 是 4×4 $\chi^{(0)}$ 的最大本征值
   你目前把 4×4 的 $\chi^{(0)}$ 取成
   mag0 = [[chi_0(1,1), chi_0(1,4)]; [chi_0(4,1), chi_0(4,4)]]
   （一个 2×2 子块）再做特征分解写进 sus0.dat。Fig.2 的说明是“noninteracting susceptibility matrix 的最大本征值”，这里指的是 4×4 $\chi^{(0)}$ 的最大本征值。两者形状相似（Γ 附近小圆环），但数值和配色刻度会不一样。

3. MATLAB 只画了中心 ±1 的小窗口
   R=1，而你的数据实际范围（你命令行里打印出来了）是
   q_x ∈ [-4.65, 4.68]， q_y ∈ [-2.67, 2.72]（单位弧度）。
   你用 xg,yg = linspace(-R,R,N) 把网格限制在 [-1,1]，因此只看到了 Γ 附近的一小块（虽然 Fig.2 重点也在中心，但想复现“整个第一布里渊区”的视觉，需要把 R 设成 BZ 的半径或把 q 归一化）。

4. 填充（化学势 μ）要对上 n=0.1
   论文 Fig.2 用的是 filling $n=0.1$。你现在 mu=-1.84d0，实际电子数（你代码里有 write(*,*) "doping=", ...）要核对是否就是 0.1。μ 不对时，那个“亮环”的半径会和论文不同。

---

具体定位（我在你的 RPA5.for 里看到的）

• 构造 $\chi^{(0)}$ 的循环（这里应该有 整体负号）：

  fortran
  chi(l1,l2,l3,l4,qx,qy) = chi(l1,l2,l3,l4,qx,qy)

& + unit(kx,ky,l4,a) * dconjg(unit(kx,ky,l1,a)) & * unit(mode(kx+qx,l),mode(ky+qy,l),l2,b) & * dconjg(unit(mode(kx+qx,l),mode(ky+qy,l),l3,b)) & * (nf(kx,ky,a)-nf(mode(kx+qx,l),mode(ky+qy,l),b)) & / (ek(kx,ky,a)-ek(mode(kx+qx,l),mode(ky+qy,l),b) + ii*eps) / n

text
这里应改成在等号右侧前面加一个负号（见下一节补丁）。

- 你输出 `sus0.dat` 的地方是写的 `w0_max`，没有人为再取负，所以数据本身就是负的，来源就是上面的符号。

- 你现在用的是自旋通道 2×2 子块：
```fortran
mag0(1,1)=chi_0(1,1)
mag0(1,2)=chi_0(1,4)
mag0(2,1)=chi_0(4,1)
mag0(2,2)=chi_0(4,4)
! 对 mag0 做特征分解 -> w2(:) -> w0_max = max(w2)
write(37,*) qx_hex, qy_hex, w0_max

如果要 一模一样复现 Fig.2，应对 4×4 的 chi_0 做特征分解并写最大本征值。

---

可直接粘贴的修复补丁（Fortran）

目标：① 把 $\chi^{(0)}$ 的符号修正；② 输出 4×4 的最大本征值
（注：只为复现 Fig.2 的“非相互作用”图；后续做 RPA 时请保持同一套符号约定）

1. 给 Lindhard 泡泡加上负号（把 + 改为 -）：

fortran
! 原来
chi(...) = chi(...) 
& + unit(...) * dconjg(unit(...)) * unit(...) * dconjg(unit(...))
& * (nf(...) - nf'(...)) / (ek(...) - ek'(...) + ii*eps) / n

! 修改后（注意就是把第一行的 + 换成 -）
chi(...) = chi(...) 
& - unit(...) * dconjg(unit(...)) * unit(...) * dconjg(unit(...))
& * (nf(...) - nf'(...)) / (ek(...) - ek'(...) + ii*eps) / n

2. 用 4×4 的 chi_0 求最大本征值并写出（替换你现在的 2×2 块）：

fortran
! ---- 在写 sus0.dat 之前，加入：对 4x4 chi_0 做特征分解 ----
chi0_real = dble(chi_0)
chi0_imag = dimag(chi_0)
call ch(4,4, chi0_real, chi0_imag, w0, matz, zr0, zi0, fv10, fv20, fm10, ierr0)

w0_max = max( max(w0(1), w0(2)), max(w0(3), w0(4)) )

! 坐标换算保持不变
px = qx1 * 2.d0 * pi / l
py = qy1 * 2.d0 * pi / l

! 非相互作用应该是正的，不要再取负号
write(37,*) px - py/2.d0,  py*dsqrt(3.d0)/2.d0,  w0_max

如果你后面还有 RPA 的 chi_uc = - chi_uc 之类的步骤，那是相互作用通道里的约定，不影响这个“非相互作用”输出；但一旦你把 $\chi^{(0)}$ 的符号修正了，整套 RPA 公式要统一同一套记号（建议等复现 Fig.2 后再核对 RPA 的号）。

3. 确保填充 n=0.1
   你代码里有

fortran
write(*,*) "doping=", count*2.d0/nn

用一个二分或牛顿法把 mu 调到 doping ≈ 0.1。最简单做法：在 mu 外套一层循环，每次根据 doping-0.1 调整 mu，直到 |误差| < 1e-4。这样圆环半径才和 Fig.2 一致。

---

MATLAB 脚本需要的小改动

1. 先把数值变成正值（用上面 Fortran 修正后就不必在 MATLAB 里再乘 -1）。如果你暂时不想改 Fortran，也可以临时在读完数据后加一行：

matlab
val = -val;   % 临时翻转符号，仅用于画图

2. 坐标归一化到“hex 单位”再取 R=1（推荐）：
   Fig.2 画的是整个第一布里渊区并且把 K 点半径归一化到 1。对蜂窝晶格，$|\Gamma \to K| = 4\pi/3$。因此：

matlab
RK = 4*pi/3;              % Γ->K 的半径
qx = qx / RK; 
qy = qy / RK;
R  = 1;                   % 现在 R=1 就是第一 BZ 边界

然后你现有的六边形 mask / 边界绘制函数就完全对齐 Fig.2 了，坐标轴 -1~1 也有物理意义。

3. 色轴正向、去掉异常值剪裁
   修正符号后，设置：

matlab
valid_data = Vg(~isnan(Vg));
caxis([0, 1.05*max(valid_data)]);
colormap(parula);   % 视觉更接近论文配色（不强制）

你脚本里“>10×median 的异常值剪裁”可能会把亮环抹掉，建议先注释掉这段。

---

验收要点（你可以照这个 checklist 自测）

• 运行 Fortran：终端打印的 doping 接近 0.1。
• 打开 sus0.dat：第三列（\lambda_\max) 为正，数量级 ~10^-2 到 10^-1（与你参数相关，不必和论文一模一样）。
• MATLAB 画图：
  • 颜色条从 0 开始到正的最大值；
  • 整个 第一布里渊区（正六边形边界）可见；
  • Γ 点周围有一条小圆环的高值带，半径随填充变小（如果你改 mu 往低填充走，亮环会往 Γ 收缩），这与 Fig.2 一致。

---

如果你只想“快速得到像 Fig.2 的图”（最小改动）

• 临时法（不改物理定义）：
  只在写文件前 w0_max = -w0_max，并在 MATLAB 里把 R 设成 4*pi/3 的归一化版本（或者先把 qx,qy 除以 4*pi/3 再令 R=1）。这会马上得到“正值 + 全 BZ + 小圆环”的图。但长期建议按上面“正式修复”把 $\chi^{(0)}$ 的负号放在正确的位置，并用 4×4 本征值来复现 Fig.2。

---

需要的话，我可以把两种版本（4×4 与 2×2 子块）的图在同一配色和坐标下并排画出来，帮你确认差别。

下面给你一份可以直接粘贴到 RPA5.for 的修改版片段，用来：

• ① 纠正 Lindhard 泡泡的整体负号；
• ② 把输出从你现在用的 2×2 自旋通道子块，改成 完整 4×4 $\chi^{(0)}$ 的最大本征值；
• ③ 写出坐标的六角投影（与你的 MATLAB 脚本一致），且不再取负号。

说明：下面的调用使用你文件里已有的 EISPACK 例程 ch（对复 Hermitian 矩阵求特征值）。如果你用的是 LAPACK，也给了等价的 zheev 写法（可二选一）。

---

1) 在累加 $\chi^{(0)}$ 的地方加“负号”

把你原来对 chi(...) 的累加这一行（注意你文件中的变量名可能是 chi, chi_0 等，只要是“k 和 k+q 循环里更新 $\chi^{(0)}$”的那一行）：

fortran
! 原来的（少了一个整体负号）
chi(l1,l2,l3,l4,qx,qy) = chi(l1,l2,l3,l4,qx,qy)
& + unit(kx,ky,l4,a) * dconjg(unit(kx,ky,l1,a))
& * unit(mode(kx+qx,l),mode(ky+qy,l),l2,b)
& * dconjg(unit(mode(kx+qx,l),mode(ky+qy,l),l3,b))
& * (nf(kx,ky,a)-nf(mode(kx+qx,l),mode(ky+qy,l),b))
& / (ek(kx,ky,a)-ek(mode(kx+qx,l),mode(ky+qy,l),b) + ii*eps) / n

改成（仅把第 2 行最前面的 + 改为 -）：

fortran
! 修改后（标准 Lindhard：前面有一个负号）
chi(l1,l2,l3,l4,qx,qy) = chi(l1,l2,l3,l4,qx,qy)
& - unit(kx,ky,l4,a) * dconjg(unit(kx,ky,l1,a))
& * unit(mode(kx+qx,l),mode(ky+qy,l),l2,b)
& * dconjg(unit(mode(kx+qx,l),mode(ky+qy,l),l3,b))
& * (nf(kx,ky,a)-nf(mode(kx+qx,l),mode(ky+qy,l),b))
& / (ek(kx,ky,a)-ek(mode(kx+qx,l),mode(ky+qy,l),b) + ii*eps) / n

---

2) 用 4×4 的 $\chi^{(0)}$ 求最大本征值并写出

把你现在那段“构造 mag0(2×2) → 求特征值 → write(37,*) ...”整块替换为下面这段。
该代码假设你当前 q 点循环里已经有 chi_0(4,4)、格点到六角坐标的 px,py，并正准备 write(37,*)。

先在变量声明区加上这些声明（若已有相同名字可忽略或合并）：

fortran
      real*8  chi0r(4,4), chi0i(4,4), w0(4), w0_max
      real*8  zr(4,4), zi(4,4), fv1(4), fv2(4), fm1(4)
      integer matz, ierr, i, j

然后在计算/写文件处放入（EISPACK ch 版本）：

fortran
C === 把 4x4 的 chi_0 拆成实部/虚部（Hermitian；允许有很小的数值不对称）
      do i=1,4
        do j=1,4
          chi0r(i,j) = dreal(chi_0(i,j))
          chi0i(i,j) = dimag(chi_0(i,j))
        end do
      end do

C === 求 4x4 的本征值（不需要特征向量）
      matz = 0
      call ch(4,4, chi0r, chi0i, w0, matz, zr, zi, fv1, fv2, fm1, ierr)
      if (ierr .ne. 0) then
         write(*,*) 'WARN: eigen-decomp failed, ierr=', ierr
      end if

C === 取最大本征值（非相互作用应为正，切勿再取负号）
      w0_max = w0(1)
      do i=2,4
        if (w0(i) .gt. w0_max) w0_max = w0(i)
      end do

C === 把 (px,py) 投到六角坐标系并写出 (qx_hex, qy_hex, lambda_max)
C     这里 px,py 是你之前的 2*pi/L 物理单位；若已是 hex 单位则直接用
      qx_hex = px - py/2.d0
      qy_hex = py*dsqrt(3.d0)/2.d0

      write(37,'(3(1x,es22.14))') qx_hex, qy_hex, w0_max

如果你更喜欢 LAPACK（zheev）而不是 ch，把上面这段换成下面这种（把 use lapack 换成你工程里包含方式；注意需要链接 LAPACK）：

fortran
C === LAPACK：对复 Hermitian 矩阵求本征值
      complex*16 chi0z(4,4), work(1)
      integer    lwork, info

      do i=1,4
        do j=1,4
          chi0z(i,j) = chi_0(i,j)
        end do
      end do

C --- workspace query
      lwork = -1
      call zheev('N','U',4, chi0z, 4, w0, work, lwork, info)
      lwork = max(16, int(dble(work)))
      deallocate(work)
      allocate(work(lwork))
C --- real call
      call zheev('N','U',4, chi0z, 4, w0, work, lwork, info)
      if (info .ne. 0) write(*,*) 'WARN: zheev failed, info=', info

      w0_max = w0(4)     ! zheev 返回升序，最大值是 w0(4)

      qx_hex = px - py/2.d0
      qy_hex = py*dsqrt(3.d0)/2.d0
      write(37,'(3(1x,es22.14))') qx_hex, qy_hex, w0_max

☑️ 记住：不要在写出时再加负号。你的 MATLAB 脚本也就无需再 val=-val。

---

3) 两个小提示（不改也能用，但建议）

• 把 $\chi^{(0)}$ 做一次 Hermitian 对称化（可减少数值噪声）：

  fortran
  do i=1,4
    do j=1,4
      chi0r(i,j) = 0.5d0*( dreal(chi_0(i,j)) + dreal(chi_0(j,i)) )
      chi0i(i,j) = 0.5d0*( dimag(chi_0(i,j)) - dimag(chi_0(j,i)) )
    end do
  end do

  这样 ch/zheev 都更稳定。

• 确保填充 $n=0.1$：如果你要和论文 Fig.2 半径一致，用一个简单的外层循环调 mu 直到 abs(doping-0.1d0) < 1d-4。

---

最后检查（快速自测）

1. 重新跑 RPA5.for：终端打印的 doping 接近 0.1；
2. 打开 sus0.dat：第三列全是正数；
3. 用你现有 MATLAB 脚本画：Γ 点出现小圆环高值带，色条为正、六角边界内图形与 Fig.2 一致。

需要我把你现在的 2×2 版本和新的 4×4 版本在同一坐标/色轴下并排画一个对比图的话，告诉我你想用的 μ（或 n），我这边直接给出对比绘图脚本即可。

你这些错误都是固定格式 Fortran（.for）排版+重复声明导致的，并不是算法问题。

发生了什么

• .for 默认是 fixed-form。在 fixed-form 里：

  • 第 1–5 列是语句标签栏，只能是空格或数字；
  • 第 6 列是续行标志（非空即把本行当作上一行的续行）；
  • 第 7 列以后才是代码。

• 你把 real(8) 直接从第 1 或第 5 列开始写，编译器把字母 r/i 认为出现在“标签栏/续行栏”，于是报：

  Illegal character in statement label field [r],[i]
  

  并且由于第 6 列非空，本行被当作上一行的续行，把

  text
  ..., w0_max
  real(8) zr(4,4), ...

  拼成了 w0_maxal（real 的 al 接在了前一行末尾），于是后面一串语法/重复维度错误全被连锁触发。

• 同时你又重复声明了已经在本程序里出现过的名字（例如 w0、ierr、i、j），因此有：

  text
  The dimensions of this array have been defined more than once [W0]
  This name has already been assigned a data type [IERR/I/J]

---

两种修法（二选一）

方案 A：继续用 fixed-form（.for）

1. 把新加的声明行起始至少缩进 6 个空格（从第 7 列开始），并且不要在第 6 列放任何字符（否则会被当续行）。
2. 不要重复声明程序里已经声明过的变量（如 w0、ierr、i、j）。若不确定，给新变量换个不冲突的名字。

示例（可直接替换你的那两行声明）：

fortran
      real*8  chi0r(4,4), chi0i(4,4), eval0(4), eval0_max
      real*8  zr(4,4), zi(4,4), fv1(4), fv2(4), fm1(4)
      integer matz
C     (如程序中已有 IERR/I/J 的声明，这里不要再声明它们)

后续求特征值的调用相应改名（把原来的 w0 换成 eval0，w0_max 换成 eval0_max）：

fortran
      matz = 0
      call ch(4,4, chi0r, chi0i, eval0, matz, zr, zi, fv1, fv2, fm1, ierr)
      eval0_max = eval0(1)
      do ii = 2,4
         if (eval0(ii) .gt. eval0_max) eval0_max = eval0(ii)
      end do
      write(37,'(3(1x,es22.14))') qx_hex, qy_hex, eval0_max

注意：行首必须从第 7 列开始，第 6 列保持空格，这样就不会再被当续行。

方案 B：改用 free-form（推荐）

把文件后缀改成 .f90 或用编译选项让编译器按 free-form 解析：

• gfortran -ffree-form RPA.f90
• ifort -free RPA.f90

然后把声明写成普通格式（不需要固定列对齐，也可用 & 续行）：

fortran
real(8)  :: chi0r(4,4), chi0i(4,4), eval0(4), eval0_max
real(8)  :: zr(4,4), zi(4,4), fv1(4), fv2(4), fm1(4)
integer  :: matz
! integer :: ierr, i, j   ! 若前面已有声明则不要重复

---

顺手避免的两个坑

• 不要再次声明已存在的数组/变量：如果程序里已声明了 w0(4)，就改名用 eval0(4) 或把原 w0 直接复用，不要再 real*8 w0(4) 一次。
• 续行（仅 fixed-form）：如果一行太长需要续行，把下一行第 6 列放任意非空字符（常用 &），并从第 7 列继续写；否则保持第 6 列为空表示新的一行。

---

一段“固定格式”安全示例（可直接抄用）

fortran
C --- declarations (fixed-form; start at col 7)
      real*8  chi0r(4,4), chi0i(4,4), eval0(4), eval0_max
      real*8  zr(4,4), zi(4,4), fv1(4), fv2(4), fm1(4)
      integer matz

C --- build chi0r/chi0i ...
      do i=1,4
         do j=1,4
            chi0r(i,j) = dreal(chi_0(i,j))
            chi0i(i,j) = dimag(chi_0(i,j))
         end do
      end do

C --- eigenvalues (EISPACK ch)
      matz = 0
      call ch(4,4, chi0r, chi0i, eval0, matz, zr, zi, fv1, fv2, fm1, ierr)

      eval0_max = eval0(1)
      do ii=2,4
         if (eval0(ii) .gt. eval0_max) eval0_max = eval0(ii)
      end do

      qx_hex = px - py/2.d0
      qy_hex = py*dsqrt(3.d0)/2.d0
      write(37,'(3(1x,es22.14))') qx_hex, qy_hex, eval0_max

按以上修改后，再编译即可消除这些错误。如果还有报“重复维度/已声明”的错误，把报错的那个名字在整个同一程序/子程序里搜一下，删掉重复的 real*8/integer 声明或换个新名字即可。

这些报错都是重复声明导致的（同一个作用域里把同名变量/数组又声明了一次）。
报错里点名的有：w0、w0_max、zr、zi、fv1、fv2、fm1、matz——说明这些名字在你的程序/子程序里之前已经声明过（可能是你先做过 2×2 本征值时就声明了），你现在又声明了一遍，所以编译器直接报：

• The dimensions of this array have been defined more than once.
• This name has already been assigned a data type.

另外，如果你的文件是 .for（fixed-form），记得从第 7 列开始写代码；第 6 列不要放字符，否则会被当成上一行的续行，容易连带出更多语法错。

---

修改方法（直接可用的补丁，避免与旧变量重名）

目标：计算 4×4 的 $\chi^{(0)}$ 最大本征值并写出；全部使用全新变量名，不和你现有的 w0、zr 等冲突。以下为 fixed-form 写法（行首 6 个空格）。

1) 变量声明（放在声明区，且第 7 列开始）

fortran
      real*8  chi0r4(4,4), chi0i4(4,4), eval4(4), eval4_max
      real*8  zrr4(4,4),  zir4(4,4),  wv1(4),  wv2(4),  wm1(4)
      integer imatz

若你的程序里已经有 ierr、i、j 的声明，就不要在这里再声明它们；下面代码直接复用你原来的 i,j,ierr。

2) 构造实/虚部并对称化（可放在每个 q 点的循环里）

fortran
C --- 把 4x4 chi_0 分解为实部/虚部，并做一次 Hermitian 对称化（更稳）
      do i = 1,4
         do j = 1,4
            chi0r4(i,j) = 0.5d0 * ( dreal(chi_0(i,j)) + dreal(chi_0(j,i)) )
            chi0i4(i,j) = 0.5d0 * ( dimag(chi_0(i,j)) - dimag(chi_0(j,i)) )
         end do
      end do

3) 求本征值（EISPACK ch 例程）

fortran
      imatz = 0
      call ch(4,4, chi0r4, chi0i4, eval4, imatz, zrr4, zir4, wv1, wv2, wm1, ierr)

      eval4_max = eval4(1)
      do i = 2,4
         if (eval4(i) .gt. eval4_max) eval4_max = eval4(i)
      end do

4) 写出到 sus0.dat（保持你的六角坐标映射，不要再取负号）

fortran
C --- 六角坐标（与你原来一致）
      qx_hex = px - py/2.d0
      qy_hex = py * dsqrt(3.d0) / 2.d0

      write(37,'(3(1x,es22.14))') qx_hex, qy_hex, eval4_max

上面代码完全避免使用 w0、w0_max、zr、zi、fv1、fv2、fm1、matz 这些已存在的名字，因此不会再触发重复声明错误。

---

如果你更喜欢 free-form（推荐）

把文件改名为 .f90 或用编译选项（gfortran: -ffree-form，ifort: -free），然后用现代声明（不受列位限制）：

fortran
real(8)  :: chi0r4(4,4), chi0i4(4,4), eval4(4), eval4_max
real(8)  :: zrr4(4,4),  zir4(4,4),  wv1(4),  wv2(4),  wm1(4)
integer  :: imatz

其余计算部分与上文相同。

---

额外提醒（避免后续再踩坑）

1. 不要重复声明：同一程序/子程序里，一个名字只能声明一次。想再做一套运算时，换新名字（如 eval4 而不是 w0）。
2. fixed-form 缩进：声明/代码从第 7 列开始；第 6 列非空表示续行。
3. 物理符号：前面我们讨论过 Lindhard 泡泡需要整体负号。别忘了在累加 $\chi^{(0)}$ 的地方把 + 改为 -（只改那一行的号即可），这样写出的本征值应为正值，与你 MATLAB 的配色一致。

照此改完再编译，以上这些“重复定义/已分配类型”的错误就会消失了。若仍有报错，把报错名字在文件里全局搜一下，确认只保留一次声明即可。