DDSCAT7.1编译及使用

DDSCAT是DDA（分立偶极子近似）方法的数值仿真程序，由普林斯顿大学的Bruce T. Draine和加利福尼亚大学的Piotr J. Flatau开发，十多年来广受计算物理学工作者的欢迎。目前的版本是7.3，项目地址：

https://code.google.com/p/ddscat/

我这里只有2010年的7.1版，两者有稍微的差别，比如7.3增加了一些新的可选算法，修改了配置文件的一些格式。

在他们项目主页上可以下载到源代码和Windows的相关版本。如果你不想折腾，就直接用编译好的程序吧。

编译

DDSCAT是用Fortran写成的，编译时应该使用Fortran编译器。我在Windows下使用的是MingW，不过默认情况下是不包含gFortran的，应该在MingW Shell里输入：

mingws-get install gfortran

速度正常的话马上就能安装好。

同样在这个shell里cd到下载解压后的ddscat/src里，手动修改Makefile文件。

Makefile里加#号的行都是注释，不会被读取执行。

搜索并定位到 OPENMP    = -fopenmp这一行，这是通知编译器在编译时连接openmp的，它是一个多线程的库，目前的电脑一般都是多核的，所以都应当加上这个选项。

也可以开启Intel的MPI支持，不过我没有安装这个库，没有使用。

定位到：

# 1.  gfortran compiler

#     sp + no MKL + no OpenMP + no MPI

下面这一段就是设置gFortran编译参数的，其中precisition是设置浮点去处是双精度还是单精度，双精度用dp，单精度用sp。我使用了默认的sp。

在OPENMP = 后面加上-fopenmp，表示使用openmp。

在LFLAGS = 这一行的行尾加上 -lgomp -lpthread，表示链接时连接mp和pthread库，这是多线程程序必需的。

好了，保存这个文件，在Mingw Shell里输入make all，如果没有错误，就说明编译成功了。如果有错误，看看是不是你的电脑缺少某些组件，比如gFortran没有安装等。

编译成功后会多出四个可执行文件：ddscat.exe，calltarget.exe，ddpol.exe，ddfield.exe。其中ddscat.exe是主程序。

下面看看ddscat/examples_exp文件夹，有一个是diel，包含了几种物质的介电常数。另外几个大写字母的文件夹是例子。你可以将ddscat.exe拷贝到相应的文件夹里，在命令行里输入ddscat，就可以正常运行，生成一些数据文件。

我们先不管例子里怎么做的，先写一个例子来说明下如果编写配置文件和处理数据。

新建一个文件夹，把ddscat.exe那四个可执行文件都拷贝到里面，然后将那些例子里任何一个ddscat.par文件拷贝过来，用文本编辑器打开，推荐ultraedit，最好不用要记事本。

[plain] view plain copy

1. ' ========== Parameter file for v7.1 ==================='   
2. '**** Preliminaries ****'  
3. 'NOTORQ' = CMTORQ*6 (DOTORQ, NOTORQ) -- either do or skip torque calculations  
4. 'PBCGS2' = CMDSOL*6 (PBCGS2, PBCGST, PETRKP) -- select solution method  
5. 'GPFAFT' = CMETHD*6 (GPFAFT, FFTWJ, CONVEX)  
6. 'GKDLDR' = CALPHA*6 (LATTDR, SCLDR)  
7. 'NOTBIN' = CBINFLAG (ALLBIN, ORIBIN, NOTBIN)  
8. '**** Initial Memory Allocation ****'  
9. 100 100 100 = dimensioning allowance for target generation  
10. '**** Target Geometry and Composition ****'  
11. 'CYLNDRCAP' = CSHAPE*9 shape directive  
12. 40 10  = shape parameters 1 - 2  
13. 1         = NCOMP = number of dielectric materials  
14. './Au_evap' = file with refractive index 1  
15. '**** Error Tolerance ****'  
16. 1.00e-5 = TOL = MAX ALLOWED (NORM OF |G>=AC|E>-ACA|X>)/(NORM OF AC|E>)  
17. '**** Integration limiter for PBC calculations ****'  
18. 1.00e-3 = GAMMA (1e-2 is normal, 3e-3 for greater accuracy)  
19. '**** Angular resolution for calculation of <cos>, etc. ****'  
20. 0.5 = ETASCA (number of angles is proportional to [(3+x)/ETASCA]^2 )  
21. '**** Wavelengths (micron) ****'  
22. 0.4 0.9 50 'LIN' = wavelengths (first,last,how many,how=LIN,INV,LOG)  
23. '**** Effective Radii (micron) **** '  
24. 0.001 0.001 1 'LIN' = aeff (first,last,how many,how=LIN,INV,LOG)  
25. '**** Define Incident Polarizations ****'  
26. (0.,0) (1.,0.) (0.,0.) = Polarization state e01 (k along x axis)  
27. 2 = IORTH  (=1 to do only pol. state e01; =2 to also do orth. pol. state)  
28. '**** Specify which output files to write ****'  
29. 1 = IWRKSC (=0 to suppress, =1 to write ".sca" file for each target orient.  
30. 1 = IWRPOL (=0 to suppress, =1 to write ".pol" file for each (BETA,THETA)  
31. '**** Specify Target Rotations ****'  
32. 0.    0.   1  = BETAMI, BETAMX, NBETA  (beta=rotation around a1)  
33. 90.    90.   1  = THETMI, THETMX, NTHETA (theta=angle between a1 and k)  
34. 0.    0.   1  = PHIMIN, PHIMAX, NPHI (phi=rotation angle of a1 around k)  
35. '**** Specify first IWAV, IRAD, IORI (normally 0 0 0) ****'  
36. 0   0   0    = first IWAV, first IRAD, first IORI (0 0 0 to begin fresh)  
37. '**** Select Elements of S_ij Matrix to Print ****'  
38. 6   = NSMELTS = number of elements of S_ij to print (not more than 9)  
39. 11 12 21 22 31 41   = indices ij of elements to print  
40. '**** Specify Scattered Directions ****'  
41. 'TFRAME' = CMDFRM (LFRAME, TFRAME for Lab Frame or Target Frame)  
42. 1 = NPLANES = number of scattering cones  
43. 0.  0. 180. 0.05  = OrderM thetamin thetamax dtheta for scattering cone 1  

这个文件分好几段，每一段前面都有*******围起来的注释。我不知道源代码里是如何读取参数的，不过最好不要改动那些你用不到的东西。第一部分是Preliminaries，是一些计算中的选项，如调用哪个算法进行迭代和FFT，是否输出二进制文件、一般没有必要改动。第二部分是Initial Memory Allocation，顾名思义是初始时的内存分配量，文档里说好像是多分配一些好，计算时会重新动态分配，因此也最好维持这个值不要改动，反正几十M内存谁也不缺。第三部分是Target Geometry and Composition，它决定了模型的参数，因此我们设置的模型应该体现在这里。'CYLNDRCAP' = CSHAPE*9 shape directive CSHAPE是使用的图形，有几十种可以选，具体类型可以在文档里查到。DDSCAT里的模型分为两大类——finite target和periodic target，分别是单个的和阵列结构，两者的配置文件写法有较大差别，我只看了单独的。单独的CSHAPE种类有二十多种，这里的CYLNDRCAP表示一种棒状结构：圆柱体两端有半球覆盖。40 10 = shape parameters 1 - 2 表示使用图形的参数，不同的图形参数都不相同，应该在文档里查阅后再填写。如上面的CYLNDRCAP里有两个参数：圆柱体的长度和直径。这里的40和10表示有多少个单位长度，一个单位长度表示分割成偶极子的尺寸。1 = NCOMP = number of dielectric materials表示有多少种介质材料，对棒来说只有一种，其他的如三个椭球相接的可以有好几种，那么这里的数值就比1多，比如是2的话需要定义每种材料的介电常数。'./Au_evap' = file with refractive index 1表示第1种材料的介电常数文件在本目录下的Au_evap文件里。这个文件是从diel里拷贝到。可以打开看看它的格式，应该很容易就能看懂。如果有好几种材料，要在下面加上第2种、第3种……的对应文件。第四部分是Error Tolerance，设置误差精度，保持默认即可。第五部分是Integration limiter for PBC calculations，周期性结构限定，对于单个的物体使用不上。第六部分是Angular resolution for calculation of , etc.，是在计算余弦平均值采用的精度，值越小越精确，一般保持默认即可，如果根本使用不上可以设置一个很大值，如10。第七部分是Wavelengths (micron)，设置计算时的波长范围。前三个值，分别是波长最小值、最大值、波长取值的数量。如果只有一个值，则前两个值相同，最后一个值取1。单位是微米。最后一个值表示最小值到最大值的分割方式，有LIN（线性）、INV（倒数）和LOG（对数）三种。第八部分是Effective Radii (micron)，它将确定分割物体的格点的大小。这里需要填写的数值是“有效半径”，即“一个格点体积对应的球体的半径”。格式同波长里一样，一般只取一种格点，所以前两个数值相同，第三个是1，最后一个不要改。第九部分是Define Incident Polarizations，入射光偏振。在DDSCAT中，有两种坐标系，一种是实验室坐标系，其X轴正向是入射方向。一种是物体的固有坐标系，一般其X轴是物体最长的方向，比如棒的纵轴方向。在建立模型前，最好先在纸上画出大致的方位。这里的三个括号分别是实验室坐标系下x，y和z的方向，默认偏振方向是y轴。下一行IORTH为1表示只计算一个入射偏振值，为2时表示再计算其垂直方向的值。第十部分是Specify which output files to write，输出文件的类型，保持默认吧。第十一部分是Specify Target Rotations，是三个角度值，分别相当于刚体力学中物体固有坐标系相对实验室坐标系的旋转角度，在第28页有详细叙述。如果你的物体不是对称的，最好让它的最长轴不与实验室坐标系的x轴重合，这样只要改变上面入射光的偏振就可以调节不同偏振的入射情况了。三个值里也都是最小值、最大值和间隔数量的顺序。第十二部分是Specify first IWAV, IRAD, IORI (normally 0 0 0)三个计算的初始值，保证默认。第十三部分是Select Elements of S_ij Matrix to Print，输出S矩阵的数量，最大为9，保持默认为6，第二行也别改，一般这就是最常用的了。第十四部分是Specify Scattered Directions，散射方向的选择，保持默认即为实验室坐标系，180度范围。修改完成后保存，在命令行里运行ddscat，即会计算出相应的值，同时输出几个文件。这个配置文件在我的电脑上跑了十分钟。今天先写这么多。