MCNP5入门（一）

评论

MCNP5入门搬运学习，搬砖永无止境

写在前面

因为有童鞋要用MCNP做毕业设计，且求助于我（orz我何德何能！）

我对于MCNP也仅仅是停留在“勉强会用”这一水平，故从核技术论坛中搬运学习MCNP5的相关知识，并结合自己的理解，完成了此篇笔记。

AAAAND，希望大家能给自己的电脑装一个好用的文本编辑器（推荐vs code！yyds！）
vs code的好处就是它可以装很多的扩展，用来提高文件的可读性（代码高亮之类）。至于怎么安装扩展，可以百度vs code汉化教程，因为汉化就是给软件安装简体中文扩展。

注：此笔记中提及的软件等请自行寻找，恕不提供。

MCNP概述

MCNP文件结构

MCNP的输入文件名不可以多于七个字符，所以不建议对输入文件增加扩展名；MCNP的输出文件的扩展名为.o；MCNP内部计算过程的记录文件的扩展名为.r。*.r 文件一般用于断点计算，在基于上一次的计算结果继续模拟计算时，MCNP将不会读取输入文件，而只读取.r文件的信息。
除了上面三种常用的文件类型，MCNP还有两种不长用的输出文件。这两种文件只有用户通过输入文件要求MCNP输出这些文件时才会得到。其中一个为.w文件，这个文件是我们想通过MCNP输出surface source 的时候才会产生；另外一个为.p文件，这个文件时我们想通过MCNP追踪粒子的输运过程的时候才会得到。

以上几个文件是常用的用户文件，下面将介绍MCNP程序内部使用的若干文件。

  1. MCNP5.exe 这是MCNP的应用程序，此程序不能双击打开，只能在DOS（cmd）下运行。
  2. vise.exe 此程序为绘图程序，可以让user观察我们在输入文件中所建立的计算模型（相对于抽象的珊元曲面卡，能可视化检查真的很不错！）。
  3. MCNPDate 该目录存放的是MCNP的截面库。
  4. xs52 该文件说明了各个核素的属性，包括质量，截面库内各种反应类型对应的截面存放的位置等。
  5. X11.dll dll文件，链接动态数据库，重要性不言而喻。

为了更好的在dos下运行mcnp，我们可以编写一些bat文件。关于bat语法可以关注我之前的博客，或者csdn等寻找学习。

注，使用bat文件时应确保环境变量已经被添加，环境变量中path中的路径决定能不能运行.cmd

MCNP使用方法

上文我们说到 xs52 文件很重要，但如果你的MCNP为绿色免安装版本的话，还需要对该文件进行简单的修改才能让MCNP找到截面库。
我们用记事本打开xs52，可以看到第一行为datapath=X:\LANL\MCNPDATA\。我们需要正确的修改该路径，例如我将我的 LANL文件夹放在了D盘根目录，此时这里的路径就是：D:\LANL\MCNPDATA\
经过上述修改之后，我们就能使用MCNP程序来进行计算了。MCNP进行计算时会占据单个cpu的资源，大家的电脑都是多核的，不用担心MCNP会占用太多资源。若想提高计算效率可以使用mpi并行计算，我尝试过一次但并没有成功orz，大家感兴趣可以自行了解尝试。

MCNP 输入文件书写方法

输入文件规范

MCNP的输入文件在结构上基本上可以分为三块，分别为cell块、surface块和其他内容。块内除了注释外，不得有空行，而两个块之间只能有一个空行。通常情况下，每一行中只能有一个主导性的助记符，这个助记符可以有很多参数和附属的次要助记符。MCNP官方手册中称这样的助记符为卡(card)。每一列不得超过80个字符，若书写不下时，可以在行尾写入 &符号，并从下行继续书写。

输入文件的第一行不论写了什么都会被MCNP忽略掉，这一行的作用是方便用户给这个输入文件做一些备注。

输入文件中有两种注释方法，助记符分别为c (c后有一个空格)和$符号。在第一行写c就是把这一行都注释掉；在某一行的某列上写$就是注释该行该列后的内容。但是即使是注释， 也不能超过每一列80个字符（包括空格）的要求。

还有三个常用的助记符，分别为r，i，j。r代表的是repeat，例如1 3r 就等于写了 1 1 1 1；i代表线性插入点，例如 1 2i 4 就等于写了 1 2 3 4；j代表的是jump，表示该位置使用默认值，例如1 j 2就等于写了 1 默认值 2，多用于描述重要性部分。

算例：7MeV轫致辐射X射线角分布、能谱

在书写MCNP输入文件前，我们首先想好要计算的模型，包括坐标系以及各个物体在这个坐标系内的相对位置。我的计算模型如下图所示。

计算模型中感兴区是半径为100cm的球内，源是能量为7MeV，位置在(-2,0,0)的单向电子束；X射线靶是半径（Y或者Z方向）1cm，高度(X方向)1cm的金；其他位置均为真空。记录面设置在圈定感兴区的球面上。

surface card

为了描述出上述模型，我们首先写surface模块，这里需要定义的surface有：

  1. 源所在面
  2. Au的三个表面（两个底面和一个侧面）
  3. 圈定感兴趣的球面
  4. 划分感兴趣区域表面的记录面

源所在面书写方法：

1
2
1    PX   -2  
C 第一个1是指这个面的序号，在cell模块会用到。PX为助记符，表示垂直于X轴的面。-2是该面在X轴上的轴距。

同理可得Au的三个面分别为：

1
2
3
4
2    PX   0
3    PX   1
4    CX   1   
C CX指的是轴心是X轴的圆柱面，1代表半径

圈定感兴趣的区域的面：

1
2
100  SO   100   
C 第一个100还是面的序号，定义面序号时不需要连续的定义，但要保证定义顺序递增。so代表的是圆心在原点的球面，后面的100就是指半径

除了上述的面，为了将记录面分割成若干个小面，我们需要额外定义一些分割用辅助面，此处我定义若干圆锥体来分割，以第一个为例：

1
2
30   KX   0 0.234567901 1   
C KX指的是轴心是X轴的圆锥面，0代表顶点在原点，0.234567901代表圆锥顶角的tan值平方。这样的圆锥有两叶，分别在Z>0和Z<0，所以最后的1代表的是我们所需要的Z>0的部分

cell card

有了这些surface之后我们就能来书写cell部分了。
栅元卡有四列数据构成：

  1. 栅元号，栅元号是定义栅元的标识符，相当于给栅元编个号，取个名字。栅元号取值区间[1,99999]。
  2. 材料号，材料号指明了该栅元是由几号材料组成的，同时需要在数据卡的材料卡里对各种材料进行详细说明。0代表真空，不用在数据卡中说明。
  3. 材料密度，指出了材料的密度，加负号是质量密度(g/cm3)，不加负号为质量密度(1024Atoms/cm3)
  4. 构成曲面，指出了栅元是由那些曲面包裹而成。除了最外部的空间，各个栅元都应该是由相应曲面包裹而成的封闭几何体。

假如1号曲面方程为f(x,y,z)=0，则f(x,y,z)<0的部分由“-1”表示，f(x,y,z)>0的部分直接由1表示。

描述栅元时还要用到布尔算符，布尔算符有三个：空格、冒号、井号。空格表示“交”，冒号表示“并”，井号表示“补”。比如100号栅元为“-1”和“-2”相交的空间。

回到我们的算例：

  1. 首先定义Au块(cell)，它由三个面围成，它们是序号分别为2，3，4。
    1. Au块由2号面的向右方向(或X轴的正方向)，3号面的向左方向（或X轴的负方向）以及4号面的内部围成。所以这个Au块（cell）的书写方法是：1 1 -19.3 2 -3 -4。
    2. 第一个1代表的是这个cell的序号；
    3. 第二个1代表的是这个cell所对应的材料种类为1，关于1号材料的定义将在后面介绍；
    4. 第三个-19.3中，负号代表的是实际密度，若没有负号则代表的是原子密度。所以-19.3代表的是Au块（cell）的密度是19.3g/cm3；
    5. 之后的2 -3 -4 描述的是Au块是由2号面的正方向，3号面的负方向，4号面的内部围成的。
  2. 源不是一个物体，所以不需要在cell块中定义。而在感兴趣内，除了Au块（cell），其他位置均为真空。
    1. 尽管是真空，我们仍需要在输入文件内把真空定义为一个物体。
    2. 定义的方法如下：2 0 -100 #1。
    3. 其中第一个2代表cell的序号；
    4. 第二个0代表定义的块内什么东西都没有；
    5. 后面的-100 #1限定了这个cell的几何空间，它表示这个块是在100号面的内部，（#助记符代表扣除的意思）除了1号cell的所有空间。
  3. 我们同样需要定义感兴区外的空间，定义的方法如下：100 0 100。
    1. 第一个100是这个cell的序号；
    2. 第二个0代表真空；
    3. 第三个100代表第100号面的外部。

2号cell和100号cell均为真空，区分哪一个是感兴区的方法是利用另外一个助记符imp，它代表的是importance，也就是粒子在该cell内的重要性。
这个助记符可以在每个cell的定义的末端写上，也可以在第三块中写入。
在cell末端写入的方法如：

1
2
3
4
2 0 -100 #1 imp:e,p=1
C 它代表了2号cell中电子(e)及光子(p)的重要性为1，即它是感兴区内的物体
100 0 100 imp:e,p=0
C 它代表了100号cell的电子及光子的重要性为0，即它在感兴区外。

感兴区与非感兴区并不需要一定是位置划分分明的两块区域，我们可以在感兴区内设定一个非感兴区，只要粒子被输运到该区域就会被杀死（类似于黑洞）。

数据卡

最后我们需要书写第三个模块，在这个模块内将指明模拟计算中所涉及的反应类型，记录方法，运算时间限定等。

  1. 模式卡：mode e p

    • 这个卡说明了计算中将涉及两种粒子，分别为电子及光子。（源粒子是电子，它轰击Au块会产生轫致辐射X射线）。

  2. 材料卡：m1 79197 -1.0

    • 其中m1代表第一号材料，同理会有m2、m3等材料;
    • 79197的格式是ZZAAA，即79是原子序数(Au)，197是质量数;
    • -1.0中负号代表的是原子个数的比重，正号则代表的是质量比重。
    • 材料卡的各种核素的比重的标识中并不要求归一化，MCNP会自动给用户归一化。也就是说 m1 79197 -2.0 79198 -2.0与m1 79197 -0.5 79198 -0.5 的写法对于MCNP而言是等效的，都说明1号材料中197Au与198Au各占一半，就像SRIM。

  3. 源定义：sdef par=3 sur=1 pos=-2 0 0 vec=1 0 0 dir=1 erg=7

    • 其中sdef是源定义的助记符代表source definition；
    • par代表源粒子的种类，1代表中子，2代表光子，3代表电子；
    • sur代表源所在的面，上面定义中sur=1说明源在序号为1的面上；
    • pos代表源所在的位置，按照上面的定义，源的位置为（-2,0,0）。请谨记若sur 和pos 同时被用上，请保证pos的点是在sur的面上，不然可能导致计算结果有误；
    • vec代表参考向量，向量的起始点是原点，而终止点是(1,0,0)。这说明了我们设定的参考向量就是X轴的正方向；
    • dir代表粒子出射方向与参考向量的夹角的余弦，上面定义中为1，说明了源粒子沿着X轴正方向出射；
    • erg代表源粒子的能量，按照上面的定义，源粒子的能量为7MeV。

    关于源的定义还有许多其他的助记符以及各种分布情况，请参考MCNP的官方使用手册（或者等我填坑(o^^o)）。

  4. 粒子重要性表示：imp:e 1 1 0

    • imp:e代表了这张卡描述的是各个cell中电子的重要性
    • 1 1 0代表了在cell块中写的各个cell（以书写顺序）的重要性分别为1、1、0。这说明了cell块中第三个描述的cell（cell的序号可能不是三，只是代表它是第三个被定义的cell。）是电子的黑洞，只要电子被输运到这个cell就会被杀死，停止输运。同理我们还可以定义imp:p 1 1 0。
    • 在前面已经描述过，这里的imp卡也可以在cell块中写入，不再赘述。

  5. 能量截断卡：cut:e j 1

    • 其中cut代表能量截断卡的助记符；
    • e代表电子；
    • j代表第一个参数选用默认值；
    • 1代表截断能量为1MeV。
    • 这个卡的作用是当输运过程中，电子的能量小于1MeV，就杀死这个电子。
    • 这样做的目的是，若我们只关心高能（>1MeV）的X射线，而能量小于1MeV的电子已经不会再产生能量大于1MeV的X射线了，所以我们可以停止这些低能电子输运来节省计算机机时而且不会影响计算结果。
    • 同理我们还可以写cut:p j 1

  6. 记录方法描述卡1：F1:p 100

    • 其中F1代表的是第一种记录模式，也就是通过某个面的粒子的个数（相对于一个源粒子）；
    • p代表记录的是光子；
    • 100代表记录面是序号为100的面。
    • 故这个记录卡的作用是记录所有通过100号面的光子。

  7. 记录方法描述2：F2:p 100

    • 其中F2代表的是第二种记录模式，也就是通过某个面的粒子的注量，它等于通过面的粒子个数（相对于一个源粒子）除以该面的面积。

  8. 记录面分割：Fs2 -30 -31 -32 -33 -34 -35 -36 -37 -38 -39 40 41 42 43 44

    • 其中Fs2代表的是分割第二号记录卡的记录面（也就是分割记录第100号面）；
    • 后面的-30 等代表分割方法，输出时MCNP会给出:
      1. 100号面上，从30号面的内部（圆锥体右叶的内部）出射的粒子
      2. 100号面上，从30号面的外部且40号面的内部出射的粒子
      3. 100号面上，从30号面的外部且31号面的外部且32号面的内部出射的粒子
      4. 100号面上，从30号面的外部且31号面的外部且32号面的外部且33号面的内部出射的粒子。
      5. 以此类推

  9. 能量分割：我们可以对任意一个记录卡分别设置能量分割卡或所有记录卡同意设置相同的能量分割卡。

    • 这里我们将对2号记录卡设置一个能量分割卡：E2 1 19i 7
    • 其中E2代表这是对2号记录卡的能量分割卡，后面的1 19i 7代表的是分割方法，也就是在1MeV至7MeV之间插入19个点（若把E2写成E0则代表对所有记录卡使用相同的能量分割方法。）。

  10. 计算时间设定卡：ctme 10

    • 这个卡说明了要求MCNP计算10分钟。
    • 可以想象，计算的时间越长，模拟的源粒子的个数就会越多，模拟结果就会越接近它的期望值，统计涨落就会越小。所以计算时间需要用户根据自己希望得到的统计可信度，设定计算时间。
    • 除了ctme还有利用模拟源粒子的个数进行时间限制的方法如NPS 10000，它说明了需要MCNP模拟10000个源粒子的输运过程（也可以写成NPS 1E4）。

在这个文件中只使用了两个记录方法，分别为F1和F2。
MCNP中一共有8中记录方法，我们可以参考使用手册，选用适合计算模型的记录卡。

此外一种记录卡可以有多个记录面，只要保证记录卡的个位数是8种记录卡的一种，例如F11和F21都是第一种记录类型，即通过某个记录面的粒子个数。

另外请谨记，在使用F1和F2的记录卡时，记录面一定要是构成某个cell的一个面，不然MCNP不会给出记录结果，也就是说如果用户随意定义了一个面(如上面定义的某个分割面)，而且这个面没有被使用在定义某个cell上，那么我们不能期望MCNP可以记录通过这个面的粒子信息。
使用分割面时，其实还是记录通过那个记录面的粒子个数，分割面只是把那个记录面分割为若个块而已。

总结-输入文件

总结上文，得输入文件为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
X-ray Ener=7MeV
c cell card
1    1  -19.3   2 -3 -4   $ Au      
2    0 -100 #1 
100   0 100

c surface card
c radiation source
1    px     -2            $ 放射源所在面
2    px     0      
3    px     1      
4    cx     1      
30   kx  0  0.234567901   1
31   kx  0  0.5625        1
32   kx  0  1.040816327   1
33   kx  0  1.777777778   1
34   kx  0  3             1
35   kx  0  5.25          1
36   kx  0  10.11111111   1
37   kx  0  24            1
38   kx  0  99            1
39   kx  0  1e33          1
40   kx  0  99           -1
41   kx  0  24           -1
42   kx  0  10.11111111  -1
43   kx  0  5.25         -1
44   kx  0  3            -1
45   kx  0  1.777777778  -1
46   kx  0  1.040816327  -1
47   kx  0  0.5625       -1
48   kx  0  0.234567901  -1
49   kx  0  0            -1
100  so    100

mode e p
m1   79197 -1.0           $ Au        
sdef par=3 sur=1 pos=-2 0 0 vec=1 0 0 dir=1 erg=7   
imp:e 1 1 0
imp:p 1 1 0
cut:e j 1
cut:p j 1
F1:p 100
F2:p 100
Fs2 -30 -31 -32 -33 -34 -35 -36 -37 -38 -39 40 41 42 43 44 &
 45 46 47 48  
E2 1 19i 7
ctme 10

注意，该输入文件共有两个空行，用来分割 cell card & surface card & data card。MCNP对于输入文件格式要求很严格，应注意。

关于vise.exe

我们可以利用vised.exe观察这个已经写好的输入文件的模型，如果书写有误，我们需要修改输入文件。
有错误时会导致vised.exe无法打开文件并被强制关闭。
出现这种情况时，请看输入文件的书写规范，如空行是否多了或少了等。

若打开后，观察到明确的红线，代表两个物体分享了一部分空间，这在MCNP中是不允许的。

在书写复杂的模型时，vised.exe可以帮助我们确定写出来的模型与我们想像的模型是一致的。我们应反复地观察各个方向的视图，从而确保输入文件的正确性，不然计算出来的结果就很有可能是错误的。

Vised.exe是windows下的应用程序，我们双击它并用它打开我们已经书写好的输入文件，就可以观察写好的模型了，操作简单不做赘述。

输出文件解读

利用MCNP 计算上面的输入文件可以得到同名的、拓展名为o的文件。使用文本编辑器打开它，搜索关键字：tally 找到如下部分：

1
2
3
4
5
6
1tally   1        nps =  2958959
           tally type 1    number of particles crossing a surface.                        
           tally for  photons  
 
 surface  100                                                                                                          
                 2.90042E-01 0.0010

这是输出文件中关于记录信息的结果部分。

  • nps 代表了一共运行了2958959个粒子的输运过程

  • tally type 1 代表记录类型是1，也就是number of particles crossing a surface

  • tally for photons 代表记录的是光子

  • surface 100 代表记录的是第100号面

  • 记录结果是：2.90042E-01 0.0010

    • 它代表有1个7MeV的电子轰击1cm的Au时，会有0.29个能量大于1MeV（我们在cut中限定了能量）光子透过第100号记录面，而这个数据的相对误差为0.001。

第二号记录结果如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
1tally   2        nps =  2958959
           tally type 2    particle flux averaged over a surface.       units   1/cm**2   
           tally for  photons  
           areas   
                surface:      100                                                                                  
               segment
                 1       6.28319E+03
                 2       6.28319E+03
                 3       6.28319E+03
                 4       6.28319E+03
                 5       6.28319E+03
                 6       6.28319E+03
                 7       6.28319E+03
                 8       6.28319E+03
                 9       6.28319E+03
                10       6.28319E+03
                11       6.28319E+03
                12       6.28319E+03
                13       6.28319E+03
                14       6.28319E+03
                15       6.28319E+03
                16       6.28319E+03
                17       6.28319E+03
                18       6.28319E+03
                19       6.28319E+03
                20       6.28319E+03
                ......
  • nps代表一共运行了2958959个源粒子的输运过程。
  • Tally type 1 代表记录类型为2，也就是particle flux averaged over a surface，请注意它的单位是1/cm2。
  • Tally for photons代表的是记录的是光子。
  • Surface 100 代表记录面是第100号面。
  • Segment 代表的是利用分割面分割后，100号面的各个子面的面积。

    由于我们使用的等立体角分割，所以各个子面的面积均相同

其他记录结果：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
surface  100                                                                                                          
 segment:         -30                                                                                                  
      energy
    1.0000E+00   0.00000E+00 0.0000
    1.3000E+00   1.58985E-06 0.0058
    1.6000E+00   1.25669E-06 0.0065
    1.9000E+00   1.01556E-06 0.0073
    2.2000E+00   8.05683E-07 0.0082
    2.5000E+00   6.50129E-07 0.0091
    2.8000E+00   5.29430E-07 0.0101
    3.1000E+00   4.41488E-07 0.0110
    3.4000E+00   3.71080E-07 0.0120
    3.7000E+00   3.09977E-07 0.0132
    4.0000E+00   2.59094E-07 0.0144
    4.3000E+00   2.11708E-07 0.0159
    4.6000E+00   1.80941E-07 0.0172
    4.9000E+00   1.59103E-07 0.0184
    5.2000E+00   1.31403E-07 0.0202
    5.5000E+00   1.04079E-07 0.0227
    5.8000E+00   8.97713E-08 0.0245
    6.1000E+00   6.86866E-08 0.0280
    6.4000E+00   4.88390E-08 0.0332
    6.7000E+00   2.82384E-08 0.0436
    7.0000E+00   1.05423E-08 0.0714
      total      8.26230E-06 0.0025
 
 surface  100                                                                                                          
 segment:          30      -31                                                                                         
      energy
    1.0000E+00   0.00000E+00 0.0000
    1.3000E+00   9.48975E-07 0.0075
    1.6000E+00   7.33985E-07 0.0086
    1.9000E+00   5.56701E-07 0.0098
    2.2000E+00   4.35679E-07 0.0111
    2.5000E+00   3.35311E-07 0.0127
    2.8000E+00   2.73079E-07 0.0140
    3.1000E+00   2.10847E-07 0.0160
    .....

其中：

1
2
surface  100                                                                                                          
        segment:         -30

代表统计的是从记录面100且从分割面30的内部（圆锥右叶的内部）通过的光子。
同理，

1
2
surface  100                                                                                                          
         segment:          30      -31

代表的是统计的是从记录面100且从分割面30的外部且从分割面31的内部通过的光子。

对于：

1
2
3
1.0000E+00   0.00000E+00 0.0000
1.3000E+00   1.58985E-06 0.0058
1.6000E+00   1.25669E-06 0.0065

最左边一列代表的是能量，第二列代表的是概率，第三列代表的是第二列的相对误差。

对于这三行统计结果的解读如下，在相应的记录面内，

  • 能量为0~1MeV的光子的注量为0，这是因为我们把能量小于1MeV的光子都杀死了（利用cut卡）。
  • 能量为1~1.3MeV的光子的注量为1.58985E-6，单位是1/cm2，且是相对于一个电子。
  • 能量为1.3~1.6MeV的光子的注量为1.25669E-06，单位是1/cm2，且是相对于一个电子。
  • 以此类推。

结束语

感觉没啥可写的了，剩下的都在手册里或者不难找到。有啥疑问我们评论区聊（评论托管服务器在国外，可能会偶尔渲染不出来大家可以试着多进几次）。

关于使用bat文件配置环境变量的问题

我们可以使用快捷键 WIN+S 搜索环境变量，打开编辑环境变量进入如下截面：

首先双击 Path，然后点击新建，将bat文件注意是bat文件所在的路径添加进去。该操作的意义是只有在Path内的路径才能使用cmd。
bat文件可以不放在MCNP根目录，只要确保bat文件与输入文件在同一目录下即可。

最后的最后，附上我写的bat文件(可以删除中间文件)，大家看看能不能用～

calculate：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
::设置输入文件和输出结果文件名
set inp=test.i
set outp=test.o

::删除旧文件
if exist  %outp%    del %outp%

::调用MCNP开始计算
call mcnp5.exe i=%inp% o=%outp%

::删除各类中间备份转储文件(应该是这么几个)
del runtp*
del dcinp
del ddinp
del dfill
del inp 
del out*
del *.sav

draw：

1
2
3
4
5
6
7
8
9
10
11
call vised.exe

del *.o
del runtp*

del dcinp
del ddinp
del dfill
del inp 
del out*
del *.sav

Paneer一个成熟稳重的沙雕网友

搬砖永无止境