写这个系列第五篇的时候又有不少人问核壳结构怎么建模。之前几篇讲的主要是球、圆柱这类单材质目标用DDSCAT内置的形状生成器就能搞定但一到核壳结构就得自己动手写target文件。离散偶极近似DDA本身不关心目标在几何上长什么样子它只认一个个离散后的偶极子位置所以核壳球、核壳圆柱这类多材料复合结构核心工作就是把目标几何翻译成DDSCAT能读的target.out。这篇把多核壳球和多核壳圆柱的Matlab建模代码完整放出来参数怎么改、格式怎么对应、有什么坑一次说清楚。1. 核壳结构为什么要自己写模型1.1 DDSCAT内置形状与FROM_FILE模式用DDSCAT的人应该都知道程序自带了几种基础形状各向同性球、椭球、圆柱、立方体、四面体等。你在ddscat.par里把目标类型选成对应的形状程序内部会自动把偶极子填满整个几何体。这种方式的优点是配置简单几行参数就能跑缺点是形状固定、材料只能有一种。哪怕你只是想要一个普通的核壳球——也就是内核一种折射率、外壳另一种折射率——内置选项也完全帮不上忙因为程序在填充偶极子的时候只会用一种材料。这时候就要走DDSCAT支持的FROM_FILE模式。简单说就是我们用Matlab或者其他工具把目标区域里的所有偶极子位置和材料类型算好写成一个文本文件让DDSCAT直接读取。程序不再负责“生成形状”只负责“读取形状”然后把电磁场求解跑完。所有复杂几何包括核壳、多核、多层、形状不规则、混合材料本质上都可以用这种方式描述。1.2 target.out文件格式DDSCAT的自定义目标文件格式不同小版本可能有细微差别但核心结构是固定的分四个部分。这里以最常见的格式说明。第一部分是第一行包含四个整数NA、NX、NY、NZ。NA是目标包含的偶极子总数NX、NY、NZ分别是包住目标的最小长方体在x、y、z三个方向上的格点数。第二部分是接下来三行每行三个浮点数分别对应三个基矢A1、A2、A3。对立方晶格来说三个基矢就是沿坐标轴的三个正交矢量数值等于偶极子间距d。比如d1nm那么三行就是“1 0 0”“0 1 0”“0 0 1”。第三部分是NA行数据每行四个整数前三个是格点索引ix、iy、iz第四列是材料编号。第k个偶极子的实际物理坐标就是r A1 * ix A2 * iy A3 * iz这就是为什么索引不需要从某个特定数字开始它只是一个格点序号真正的空间位置由基矢和序号共同决定。在立方晶格下更直观第ix个点在x方向离原点ix*d的位置iy和iz同理。第四列材料编号在最简单的情况下全部写1表示整个目标单一材料。如果要做核壳结构就要让内核位置的材料编号为1外壳位置的材料编号为2这样DDSCAT就知道不同偶极子用不同的介电函数去算极化率。1.3 多材料区分与ddscat.par的关系生成target.out只是第一步DDSCAT如何知道编号1、编号2各自对应什么折射率这就要在ddscat.par里做两件事。第一把介电函数数量设成2不同版本关键字不同通常是NCOMP或者类似参数也就是告诉程序“这个目标不是单一材料”。第二准备两个折射率文件或者介电函数表分别对应材料1和材料2在par文件里按序号填进去。有的版本直接在par文件里写折射率随波长的表格有的版本指向外部数据文件操作逻辑都一样第1个文件对应材料编号1第2个对应材料编号2。我见过不少人卡在这个环节target文件明明写好了跑出来的结果却是完全均匀的球。排查下来基本都是par文件里只保留了一个介电函数表程序只认材料编号1外壳的编号2全部当成真空或者默认材料处理。所以写完target文件之后第一件事就是检查par文件里的介电函数数量是不是大于等于target文件里用到的最大材料编号。2. 多核壳球模型设计Matlab实现2.1 几何设计与判据推导先明确什么叫多核壳球。我这里的定义是一个球形外壳内部包含多个小球形内核外壳和内核对不同波长的折射率不同甚至内核之间材料也可以不同。这类结构在催化剂载体、表面增强拉曼散射基底、多色发光纳米颗粒里都很常见。建模逻辑其实非常简单核心就是一个三维距离判断。对于空间中任意一个格点先判断它到球心的距离是不是大于外壳半径大于就丢弃这个点不属于目标然后判断它是否落在任意一个内核球里面判断方法是计算该点到所有内核球心的距离只要有一个距离小于内核半径就标记为材料1剩下的点都在外壳内标记为材料2。判断顺序很关键一定要先判断外壳再判断内核因为内核区域在物理上属于外壳内部的一个子区域不能提前排掉。内核个数和位置的设置是这个模型最灵活的部分。内核位置可以用N行3列的数组表示每一行是一个内核球心的坐标。要注意的是所有内核区域不能重叠并且每个内核必须完整地落在球形外壳内部。如果内核贴到外壳边界甚至凸出去物理上就会出现一个“破了洞”的壳层仿真结果会变得很奇怪而且不容易发现。2.2 完整Matlab代码这里给出一个可以直接运行的Matlab脚本。参数部分放在最前面方便按自己需求修改。代码在MATLAB R2019b之后的大部分版本上都可以直接跑不依赖任何工具箱纯基础语法。% % 多核壳球模型生成 DDSCAT target.out % 输出文件: target_multi_core_shell.out % clear; close all; clc; %% 1. 模型参数 d 1.0; % 偶极子间距晶格常数 R_shell 24; % 球形外壳半径单位为d的倍数 R_core 8; % 内核球半径 N_core 4; % 内核个数 % 内核球心位置每行一个 [x y z] core_centers [ 10, 6, 4; ... -10, 4, -6; ... 0, -10, 8; ... 5, -5, -12 ]; %% 2. 网格范围设置 % NX, NY, NZ 是包住目标的最小长方体在三个方向的格点数 NX round(2*R_shell) 5; NY round(2*R_shell) 5; NZ round(2*R_shell) 5; % 几何中心对应的索引偏移 ox (NX-1)/2; oy (NY-1)/2; oz (NZ-1)/2; % 预分配存储最大不超过包围盒体积 max_dip NX * NY * NZ; target zeros(max_dip, 4); cnt 0; %% 3. 扫描所有格点 for ix 1:NX x (ix-1-ox)*d; for iy 1:NY y (iy-1-oy)*d; for iz 1:NZ z (iz-1-oz)*d; % 3.1 外壳球判断 r_sq x^2 y^2 z^2; if r_sq R_shell^2 continue; % 在目标外部跳过 end % 3.2 内核判断 in_core false; for nc 1:N_core dx x - core_centers(nc,1); dy y - core_centers(nc,2); dz z - core_centers(nc,3); if dx^2 dy^2 dz^2 R_core^2 in_core true; break; end end cnt cnt 1; if in_core target(cnt, :) [ix-1, iy-1, iz-1, 1]; % 材料1内核 else target(cnt, :) [ix-1, iy-1, iz-1, 2]; % 材料2外壳 end end end end target(cnt1:end, :) []; fprintf(偶极子总数: %d\n, cnt); %% 4. 输出 target.out fname target_multi_core_shell.out; fid fopen(fname, w); fprintf(fid, %d %d %d %d\n, cnt, NX, NY, NZ); fprintf(fid, %e %e %e\n, d, 0, 0); fprintf(fid, %e %e %e\n, 0, d, 0); fprintf(fid, %e %e %e\n, 0, 0, d); for k 1:cnt fprintf(fid, %d %d %d %d\n, ... target(k,1), target(k,2), target(k,3), target(k,4)); end fclose(fid); fprintf(目标文件已生成: %s\n, fname); %% 5. 可视化检查 figure; scatter3(target(:,1), target(:,2), target(:,3), 6, target(:,4), filled); axis equal; grid on; colormap([0.8 0.2 0.2; 0.2 0.2 0.8]); colorbar; xlabel(ix); ylabel(iy); zlabel(iz); title(多核壳球 target 可视化);2.3 代码逻辑解读这段代码里藏了几个细节值得展开说一下。索引偏移这块上面代码里x(ix-1-ox)*d其中ox(NX-1)/2。这个写法是把NX个格点排列在0到NX-1的索引区间然后通过偏移ox让物理坐标的零点落在包围盒中心。这么做的好处是球心刚好和坐标原点重合判断距离的时候可以直接用x²y²z²不用额外减球心坐标。如果你想让球心放在别的位置直接改ox、oy、oz这三个偏移量就行。内核判断用了提前break一旦确认当前点在任意一个内核内部就立即跳出循环不再继续判断后面的内核。内核数量少的时候影响不大但如果内核数量多、又在大循环里跑这个break能省下不少时间。我试过内核数量超过二十个、网格规模接近百万量级的情况提前break和不加break的耗时差距非常明显基本能快两三倍。预分配用的是包围盒体积NXNYNZ因为实际偶极子数一定小于等于这个数。跑完之后用target(cnt1:end, :) []把多余的行删掉。这只是为了在Matlab里避免循环内动态增长数组带来的卡顿不影响最终结果。可视化部分是很多人容易忽略的。生成target文件之后我强烈建议先scatter3画一下用材料编号着色。这样一眼就能看出内核是不是完整包在外壳里面、有没有位置偏移、边界是不是干净。我通常先看可视化的结果确认没问题再拿去做DDSCAT仿真。省一次错误仿真就省几十分钟到几小时。3. 多核壳圆柱模型设计Matlab实现3.1 圆柱壳与球核的几何判据多核壳圆柱稍微复杂一点因为圆柱的形状判断比球复杂。我这里的模型是一个沿z轴方向的圆柱形外壳内部放置多个球形内核。圆柱半径为R_shell高度为H_shellz坐标从-H_shell/2到H_shell/2。圆柱判据是两个条件同时满足径向距离rho sqrt(x² y²)不超过R_shell轴向距离|z|不超过H_shell/2。只有两个条件都满足这个格点才属于外壳区域。然后内核判断和球形情况一样只要落在任意一个球形核内就标记为内核材料。需要特别强调的是内核球体的位置必须是圆柱外壳的一个子集。因为圆柱的约束条件和球不同一个内核可能在三个方向上都满足圆柱的径向约束但在轴向方向超出H_shell/2范围就会裸露出一部分在壳体外面。设计内核位置的时候要留出足够余量最简单的做法是让内核球心距离圆柱轴线的距离加上内核半径小于圆柱半径同时内核球心的z坐标绝对值加上内核半径小于H_shell/2。如果实际模型里形状偏长比如圆柱半径只有几个d、高度却有好几十个d那么网格重心偏移也要相应调整。代码里ox和oy由半径决定oz由高度决定别把三条轴的包围盒都设成一样的值不然整个网格区域里大量格点都在目标外部白算一圈循环。3.2 核心代码改动多核壳圆柱的Matlab代码整体结构和球版本差不多关键是几何判据不同。这里给出完整代码。% % 多核壳圆柱模型生成 DDSCAT target.out % 圆柱轴沿z方向 % 输出文件: target_multi_core_cylinder.out % clear; close all; clc; %% 1. 模型参数 d 1.0; % 偶极子间距 R_shell 20; % 圆柱外壳半径 H_shell 60; % 圆柱外壳高度沿z方向 R_core 6; % 内核球半径 N_core 4; % 内核个数 % 内核球心位置需保证所有内核都完整包在圆柱内 core_centers [ 8, 0, 10; ... -8, 0, -10; ... 0, 8, 20; ... 0, -8, -20 ]; %% 2. 网格范围设置 NX round(2*R_shell) 5; NY round(2*R_shell) 5; NZ round(H_shell) 5; ox (NX-1)/2; oy (NY-1)/2; oz (NZ-1)/2; max_dip NX * NY * NZ; target zeros(max_dip, 4); cnt 0; %% 3. 扫描所有格点 for ix 1:NX x (ix-1-ox)*d; for iy 1:NY y (iy-1-oy)*d; for iz 1:NZ z (iz-1-oz)*d; % 3.1 圆柱外壳判断 rho2 x^2 y^2; if rho2 R_shell^2 || abs(z) H_shell/2 continue; end % 3.2 内核判断 in_core false; for nc 1:N_core dx x - core_centers(nc,1); dy y - core_centers(nc,2); dz z - core_centers(nc,3); if dx^2 dy^2 dz^2 R_core^2 in_core true; break; end end cnt cnt 1; if in_core target(cnt, :) [ix-1, iy-1, iz-1, 1]; else target(cnt, :) [ix-1, iy-1, iz-1, 2]; end end end end target(cnt1:end, :) []; fprintf(偶极子总数: %d\n, cnt); %% 4. 输出 target.out fname target_multi_core_cylinder.out; fid fopen(fname, w); fprintf(fid, %d %d %d %d\n, cnt, NX, NY, NZ); fprintf(fid, %e %e %e\n, d, 0, 0); fprintf(fid, %e %e %e\n, 0, d, 0); fprintf(fid, %e %e %e\n, 0, 0, d); for k 1:cnt fprintf(fid, %d %d %d %d\n, ... target(k,1), target(k,2), target(k,3), target(k,4)); end fclose(fid); fprintf(目标文件已生成: %s\n, fname); %% 5. 可视化检查 figure; scatter3(target(:,1), target(:,2), target(:,3), 6, target(:,4), filled); axis equal; grid on; colormap([0.8 0.2 0.2; 0.2 0.2 0.8]); colorbar; xlabel(ix); ylabel(iy); zlabel(iz); title(多核壳圆柱 target 可视化);3.3 扩展到圆柱形核如果你实际的模型不是球形核而是多个圆柱形核放在同一个圆柱外壳里那判断逻辑只需要把球判据换成圆柱判据就行。比如第n个核是一个半径rcore_n、高度hcore_n、轴方向沿z的小圆柱那么判据就是(x-cx_n)² (y-cy_n)² ≤ rcore_n² 且 |z-cz_n| ≤ hcore_n/2两个条件同时满足才标记为内核材料。如果是沿x轴或者y轴取向的圆柱核就交换对应的坐标公式本质思路完全一样。在代码里实现的时候把3.2那一整块循环体替换成圆柱判断即可。这种多圆柱核结构其实在某些光纤传感、纳米柱阵列模型里会出现不过注意方向对齐很重要。如果几个圆柱核的轴方向不一样比如一个横着一个竖着那么判断条件不能共用同一组公式得给每个核单独写判断逻辑。这种复杂情况我建议单独封一个判断函数比如isInsideCore(x,y,z,coreParams)每个核传不同的参数类型进去代码会干净很多。4. 生成target后如何调用DDSCAT4.1 偶极子数量与有效半径估算很多新手第一次生成target文件后会看到一大串偶极子数量直接懵了这么多偶极子要多少内存能不能跑起来这里给一个简单有效的估算方法。偶极子总数NA近似等于目标实际体积除以偶极子元胞体积d³。球壳的有效体积要按外壳体积减去内核体积来算。比如R_shell24、R_core8、四个内核的球形目标外壳体积近似为4/3π(24³ - 4×8³)算下来大约5.7万偶极子数也就在这个量级附近。有效半径由ae (3V/4π)^(1/3)给出这个参数在ddscat.par里要用到建议先算好。网格规模、偶极子数和内存占用的大致关系我整理成了表格。这里的“内存”是个人电脑上跑DDSCAT的经验值考虑了极化率张量、电场迭代等主要数组的开销。外壳半径/格点数偶极子数估算典型内存占用是否适合个人机R10约4200约1-2 MB轻松R20约33500约10-20 MB轻松R30约113000约40-70 MB轻松R50约524000约200-300 MB可以跑R80约2145000约800 MB-1.2 GB比较吃力实际内存还受迭代算法和输出选项影响但用这个表做预判足够。我的经验是偶极子数超过两百万之后迭代收敛会明显变慢最好先在低分辨率下跑通流程再逐步加密网格。4.2 ddscat.par关键配置target文件生成之后仿真阶段的配置主要集中在ddscat.par里。不同DDSCAT版本的关键字名称有变化但逻辑一致这里说几个核心点。第一目标类型选FROM_FILE把生成的目标文件名填进去。有的版本需要在par文件中写文件路径有的版本是通过命令行参数指定按你手头版本的manual操作。第二介电函数数量改成2对应target文件里的材料编号1和2。漏掉这一步是核壳仿真翻车的第一大原因一定检查一下par文件里实际加载的介电函数表是不是真的有两个。第三有效半径ae填入估算值。这个值作为归一化尺度影响后续消光效率、吸收效率等结果的定义算错会导致所有输出截面数值整体偏移。建议用Matlab脚本直接算出ae再填避免手算出错。第四入射波长要和d的单位匹配。如果d1nm波长就写多少纳米如果d1μm波长就写多少微米。单位不一致的话仿真结果看起来像是跑出来了实际物理意义完全不对。4.3 网格规模与仿真时间权衡网格加密带来的偶极子数量增长是三次方关系。d从1变成0.5三个方向各翻倍偶极子数量变成8倍仿真时间往往不止8倍因为迭代收敛所需步数也可能增加。所以确定d的时候要先想清楚目标结构的特征尺寸。一个实用的建议是先跑一个d比较粗的target比如核半径和壳半径只有两三个格点确认从建模到仿真整个链路没问题。然后逐步加密。这样做的好处是如果前面配置有错比如par文件介电函数没设对在小规模下几分钟就能发现而不是等两个小时的仿真结束后才发现结果不合理。我在调新结构的时候基本都会这么做很少上来就直接跑细网格。5. 常见问题与排查技巧5.1 偶极子数量爆炸或内存不足偶极子数量太多导致内存溢出是最常见的问题。如果你发现目标文件体积很大、偶极子数远超预期先别急着加内存检查三件事一是包围盒NX、NY、NZ是不是开得过大导致大量空格点被扫进来二是d的值是不是定得太小结合目标实际尺寸重新确认一遍三是目标填充阶段是否误把外部区域也写进了target文件通过可视化检查一眼就能看出来。5.2 核壳边界锯齿与最小偶极子间距要求DDA要求偶极子间距d满足d ≤ λ/(2π|m|)左右这是精度层面的约束。但如果d相对核壳厚度来说太大边界就会非常粗糙内核和外壳的交界呈现明显阶梯状。这种锯齿边界在远场散射计算里影响不大但在近场增强、吸收效率峰位这类对边界敏感的计算里会带来不可忽略的误差。解决方法是减小d。比如外壳半径30nm、内核半径10nm的结构d1nm时外壳层只有20个格点厚度d0.5nm时就有40个格点边界精度显著提升。代价是偶极子数变成8倍需要结合计算机性能和精度需求做取舍。如果发现d1时外壳厚度只有5个格点左右我建议一定要加密否则“壳层”在离散后会变得非常模糊。5.3 材料区分失效很多人跑完发现核壳结构的散射谱和裸球几乎一样第一反应是target文件错了但问题往往出在par文件。核心检查点就是par文件的介电函数数量。如果你在target文件里写了材料编号1和2但par文件里只配置了一个介电函数表DDSCAT要么报错要么把所有偶极子都当材料1处理外壳就“不存在”了。另外不同DDSCAT版本对composition index的定义有细微差别。老版本有的需要把材料编号写在某一列固定位置新版本可能支持更多列扩展信息。写代码之前翻一下手动确认你用的版本ENCODING格式。我见过有用户在7.3版本里用老格式的target文件偶极子坐标完全错乱仿真结果自然不可能对。5.4 可视化验证与调试技巧生成target文件之后我一直习惯用scatter3做可视化再叠加一个透明网格或者等值面来检查。Matlab里直接用scatter3画几万个点并不慢而且能很直观地看到目标形状和材料分布。调试内核位置时可以先把内核球心坐标单独打印出来用scatter3画在外壳内部肉眼检查是否越界。如果内核比较多也可以用矩阵运算一次性判断所有内核是否都满足边界条件内核球心到外壳球心的距离加上内核半径必须小于外壳半径这个条件一旦不满足就在命令行打印警告。在脚本里加这种自动检查能少踩很多坑。还有一个实用小技巧每次修改参数后生成的target文件最好带上参数后缀比如target_R25_r8_4core.out。这样后面跑仿真的时候一看文件名就知道用的什么几何。我刚开始做这个系列时吃过亏同一目录下覆盖了多个target文件最后跑的哪个自己都分不清浪费时间重跑了好几组。我对核壳结构建模最大的感受是几何生成本身并不难难的是把几何数据、材料配置、仿真参数三者的对应关系理清楚。target文件里第四列写的是“材料编号”par文件里配置的是“材料物理参数”两者通过编号关联任何一个环节错位整个仿真结果就没有意义。如果你现在正卡在核壳模型上按这篇的思路把几何判据、代码、par配置三个环节挨个过一遍基本上能把问题定位出来。实际项目中如果遇到更特殊的核壳复合结构比如多层壳、随机多核、嵌入非球形核也无非是把这套判据逻辑扩展一下核心思路不变。