（凝聚态物理专业论文）gpu运算在stm模拟中的应用.pdf

中国科技大学硕士学位论文 摘要 i 摘摘 要要 近年来，图形处理器（graphics processing unit,gpu）的发展相当迅速，gpu 的运算能力及存储带宽均已远远超过目前主流cpu。 将gpu作为cpu的协处理器 完成大规模数据密集型的计算任务，相对于集群和超级计算机的实现，具有很高 性能和性价比。 现在，stm在表面科学研究中应用非常广泛，其中对stm图像的理论模拟十 分关键。但是做大量图像模拟时，运算量会比较大，耗时较长。因此，本文基于 cuda（computer unifie device architecture，统一计算设备架构）平台，利用 gpu来加速stm图像模拟。 同时，我们还对金属酞菁（mpc）分子的负微分电阻效应（ndr效应）做了 一些研究，基于我们开发的stm模拟方法，模拟了金属酞菁分子的sts谱，进一 步研究了ndr效应的形成机制及应用。对负微分电阻效应的研究有一定的促进 作用。 第一章，我们首先从gpu的硬件特征入手介绍gpu的基本架构，以及gpu的 硬件模型的发展。然后我们讲述了gpu用于通用计算从基本的图像处理到 gpgpu通用计算再到cuda的发展历程。最后我们详述了cuda的设计原理， cuda软件的架构，cuda的线程模型，存储结构，编译方法等。 第二章，我们首先介绍stm的基本工作原理，接着我们对stm模拟的几种基 本方法做了概述。首先，我们讲述了修正的bardeen方法（mba方法），以及这 种方法的一个特例，也就是应用比较广泛的tesoff-hammann方法，最后介绍了 基于bardeen方法发展的一个新的stm模拟方法，这种方法也是本文gpu模拟所 采用的方法。 第三章，我们首先通过两个测试体系验证了程序的正确性。接着我们介绍了 程序的基本算法，并分析了程序中部分重要的代码。然后对程序各部分的加速进 行测试，得到不同部分在gpu中的加速情况。最后，对整体加速，我们比较了两 种算法。一种应用fftw傅里叶变换库，另一种应用cuda的傅里叶变换。我们 通过不断加大程序的运算量来考察gpu的总体加速效果。 第四章，我们基于现在的程序对负微分电阻效应进行一些研究。我们首先阐 述了负微分电阻效应的几种机理， 重点介绍了以前胡振芃等人的理论工作。 然后， 我们基于此基础模拟的几种酞菁分子做了理论计算，解释了模拟的结果，并分析 各个ndr效应的产生机理以及将来可能的应用。接着我们计算了自旋极化的电 子遂穿情况，更加细致的说明了ndr效应的产生机理。 中国科技大学硕士学位论文 摘要 ii 关键字关键字: 图形处理器 cuda stm模拟 修正的bardeen方法 酞菁分子 负微分 电阻效应 中国科技大学硕士学位论文 abstract iii abstract in rencent years, gpu (graphics processing unit) has developed very rapidly. gpu has a large advantage to cpu in computation capacity and bandwidth. using gpu, together with cpu, in large-scale data-intensive tasks has better performance and is more cost-effective than clusters and supercomputers. stm has been widely used in surface science research, with the help of stm image simulation. however, intensive stm simulation can become very time consuming. in this thesis, we focuse on accelerating stm simulation with gpu using cuda (computer unifie device architecture), which is a gpu computation platform. we also study the ndr effect in mpc moleculars. we simulate the sts of mpc with our stm simulation method. mechanism and application of the ndr effect are discussed. in the first chapter, we introduce the gpu hardware firstly, including the gpu basic architecture and its recent development. then, we describe how gpu is used as a general computation tool from gpgpu to cuda. finally, we introduced in details about cuda, including the basic design principles, software architecture, thread mode, memory structure, and debug principles. in the second chapter, we first talk about the working principles of stm, and then introduce several basic stm simulation methods briefly. they includes mba (modified bardeen approximation) method and the tesoff-hammann method which can be considered as a special case of the mba method. in the end of the chapter, we introduce a new stm simulation method which developed from bardeen method. in the third chapter, we first verify the correctness of our program through two testing systems. then, we introduce the algorithm of our program, and analyze some important lines of the codes. finally, we test the speed-up ratio of each part, and also the total acceleration of our program. in this part, we compare two implementations, one of which use fftw, and the other use cufft. the final chapter is on the ndr effect. we first introduce several mechanisms of ndr, focusing on the zhenpeng hus work in 2007. and then, we made some simulations based on hus work. we explain the simulation results, and forecast the applications of ndr. we also consider the polarization situation. 中国科技大学硕士学位论文 abstract iv key words: gpu(graphics processing unit) cuda(computer unifie device architecture) stm simulation mba(modified bardeen approximation) mpc molecular ndr effect 中国科学技术大学学位论文原创性声明 本人声明所呈交的学位论文,是本人在导师指导下进行研究工作所取得的成 果。 除已特别加以标注和致谢的地方外，论文中不包含任何他人已经发表或撰写 过的研究成果。 与我一同工作的同志对本研究所做的贡献均已在论文中作了明确 的说明。 作者签名：_ 签字日期：_ 中国科学技术大学学位论文授权使用声明 作为申请学位的条件之一， 学位论文著作权拥有者授权中国科学技术大学拥 有学位论文的部分使用权，即：学校有权按有关规定向国家有关部门或机构送交 论文的复印件和电子版，允许论文被查阅和借阅，可以将学位论文编入中国学 位论文全文数据库等有关数据库进行检索，可以采用影印、缩印或扫描等复制 手段保存、 汇编学位论文。 本人提交的电子文档的内容和纸质论文的内容相一致。 保密的学位论文在解密后也遵守此规定。 公开 保密（_年） 作者签名：_ 导师签名：_ 签字日期：_ 签字日期：_ 中国科技大学硕士学位论文 第一章 绪论 1 第一章第一章 绪论绪论 1.1 基于图形接口的通用 gpu 运算 gpu（graphics processing unit，图形处理器）主要任务是加速图形图像处理 速度，是显卡中非常重要的组成部分。在20世纪80年代之前gpu还没有出现，计 算机图形处理的任务主要由cpu负责完成，后来ge（geometry engine）1图形 处理器的出现，使得图形处理任务从cpu独立出来转向gpu。 随着gpu硬件上的发展，gpu已经不止局限于图形处理。随着其性能和可编 程架构的完善，以及存储带宽的加大，将gpu应用于通用计算，使gpu作为cpu 的协处理器进行科学运算已经有了长足的发展。 由于gpu在硬件结构上与cpu的 差别，主流gpu 都采用了统一架构单元，拥有强大的可编程流处理器阵列以及 较大的存储器带宽，这些条件使得gpu更加适合于浮点运算。如图1.1和图1.2 3 是gpu与cpu在浮点运算能力和存储带宽的对比。由图1.1可以看出，gpu的单 精度浮点运算能力远远高于同时期主流cpu，运算效率相差10倍左右。同时双精 度浮点计算能力也相当于同时期cpu的两倍以上。图1.2显示了gpu的存储器带 宽相当于同时期cpu的5倍左右。 将gpu从传统的图形处理器的概念中解放出来，使之应用于数值计算，能够 大大提升运算效率，这也成为当今计算机领域的一个研究热点。 图图1.1 gpu与cpu的计单精度浮点计算能力发展趋势(摘自cuda program guide 3.1) 中国科技大学硕士学位论文 第一章 绪论 2 图图1.2 gpu 与cpu 的存储器带宽发展趋势(摘自cuda program guide 3.1) 注：注：cpu存储器带宽指的是cpu与内存（通过dimm 插槽与主板相连并具有可扩展性）的 数据传输带宽，gpu存储器带宽指的是gpu与显存（固化在显卡pcb 板上的gddr 存 储器颗粒）的数据传输带宽。 1.2 stm 模拟 1982年， ibm公司在瑞士zurich实验室的科学家binnig 456和rohrer 等人根 据量子力学中电子隧穿原理，设计制作了第一台扫描隧道显微镜（stm）。这种 新型显微术拥有很高的分辨率，如图1.3所示，stm的横向分辨率达到0.1 nm，而 纵向分辨率更达到0.01nm， 这样的分辨率足以使我们观察原子和分子的表面特征， 这也使得人类第一次真实的认识了原子和分子。 这是以往所有检测手段不能匹敌 的。 stm不仅能直接反映样品实空间的信息， 而且能够反映固体表面电荷密度的 分部情况。stm超高的分辨率对于我们分析研究表面缺陷、原子分子的吸附、固 体电子态等有极大的推动作用。此外，由于stm对环境的要求并不严格，使用范 围比较广泛，这样可以进行大气、液体等特殊环境下的研究。stm从分辨率到适 用条件相对于其他显微镜都有绝对的优势， 使之成为纳米研究中不可缺少的工具。 应用stm做出的科研成果也层出不穷。 中国科技大学硕士学位论文 第一章 绪论 3 图图1.3 探测手段的应用范围的比较(摘自赵爱迪博士论文) stm发明以后，针对stm做的一系列理论研究从来没有停止过，从简单的 tessof-hamann近似，到bardeen方法，以及green函数法等，随着计算机技术的 发展， 应用这些方法对体系做的stm模拟也不断的发展。 但是总体来讲要得到一 幅stm的模拟图像， 需要处理针尖和样品费米面附近的所有波函数， 如果对于比 较大的体系，计算量是非常大的，也比较耗时。2005年paz等人7，发展了一种 stm模拟方法，该方法通过快速傅里叶变换，能够大大提高stm模拟的效率。 胡振芃8后来又对该方法做了一定的改进，使之更具有物理意义。本文就是基于 胡的方法，对该方法做了gpu加速处理，大大提高了模拟速度。 1.3 ndr 效应 负微分电阻效应（ndr效应），即半导体器件中电流随电压增加而减小，而 呈现电阻增大的现象。这种效应基于量子力学中的电子隧穿原理，隧道二极管具 有超短的响应时间 (10-1310-16 s)， 在高速电路中有着广泛的应用。 在研究ndr 效应中， 用stm的点模式测量样品中某点的i-v谱(sts谱)是一个重要的研究手段。 对于ndr效应的成因有各种不同的解释，陈岚等9利用sts谱的模拟提出了新的 ndr效应的产生机理。本文继续他们的工作，用我们自编的程序模拟sts谱，进 一步阐述该机理。 中国科技大学硕士学位论文 第二章 gpu 通用计算基本理论 5 第二章第二章 gpu 通用计算基本理论通用计算基本理论 2.1 引言 随着图形处理器硬件的发展，gpu通用计算(gpgpu)越来越受到科学计算领 域专家学者的关注。图形处理器(graphics processing unit，gpu)最先由nvidia 公司于1999年推出。第一代图形处理器geforce256，首次实现了三维几何变换 和光照计算(transformlighting，t 其中 为波矢，电子运动与经典情况相同。电子能够自 由的通过势垒发生电子转移。 当电子能量 ，也就是处于经典禁区时，方程的解为： ； 其中 。电子出现在势垒x点处的几率为 ，可见电子有遂穿的几率。 现在我们把上面的理论应用到如图 3.4 所示的理想一维金属真空金属 隧道结模型，这个模型对应 stm 中针尖真空样品隧道结，可以得到其中隧 穿现象的定性描述。当所加偏压 v 远小于功函数（设针尖与样品功函数相同） 中国科技大学硕士学位论文 第三章 stm 模拟基础 25 时，只有处于 与 之间本征能量为 .的样品电子态为 的电子能够 隧穿进入针尖，其出现在 x = d 处的几率可表示为： (3.2) 其中 为样品电子态 在表面处的数值，衰减因子为： (3.3) 将 与 之间的所有的表面电子态相加，便得到遂穿电流的值，即： (3.4) 根据样品的局域电子态密度ldos的定义： (3.5) 如果v足够小以至于 范围内电子态密度没有明显变化，式2.5的求和可 方便地改写为 处的样品ldos： (3.6) 图图3.4 一维金属（样品）真空金属（针尖）隧穿模型（摘自胡振芃博士论文） 对于常用的金属表面功函数约5 左右，代入式3.6就可以得到隧道电流随针 尖样品距离指数变化的衰减因子 。也就是说，距离每增大1，隧道 电流就减小 倍。可见遂穿电流对针尖和样品的距离非常敏感。 考虑到费米面电子态密度： (3.7) 隧道电流可表达为： (3.8) 中国科技大学硕士学位论文 第三章 stm 模拟基础 26 以上的一维近似处理方法可以用来定性地说明stm 中隧道电流i 与针尖 样品间距d 和样品费米面附近电子态密度的依赖关系。 3.4 微扰理论近似和 bardeen 公式 简单的一维隧穿模式只能对stm的工作原理做定性的解释，但是往往实验中 都是涉及的三维复杂体系。简单的一位模型显然不能满足要求。据此chen 24等 基于bardeen 23发展的微扰理论发展了修正的bardeen方法（modified bardeen approach, mba）。其基本原理为：考虑针尖和样品相互有一定距离，只存在量 子隧穿效应， 没有直接的相互作用，可以把针尖和样品看做近似独立的两个弱耦 合子体系。通过微扰理论来处理这两者之间的耦合，如图3.5所示。这种处理方 法利用了含时微扰理论，最早是由bardeen提出，因此被称为bardeen 方法 （bardeen approach），stm发明以后，chen等将该方法应用于stm体系，并对 其做了一定的修正，使其更加适用于stm 体系，也就是修正的bardeen方法 （mba）。 图图3.5 bardeen 微扰理论示意图（源自chen 1993） 注注：(a) 样品针尖遂穿的耦合系统。(b) 和 (c) 将耦合系统拆分成两个子系统，通过子系统 波函数的重叠，用femi黄金规则计算隧道电流。 使用mba处理stm体系时，先分别求解两个子系统的schrdinger 方程，得 两个子系统中电子的两组定态解 和 （s表示样品，t表示针尖）。然后将 两组定态波函数代入遂穿矩阵元mst： (3.9) 中国科技大学硕士学位论文 第三章 stm 模拟基础 27 其中 为针尖和样品之间任意一个分隔面。从 态到 态的遂穿几率可以由费米 黄金规则（fermi golden rule）得出： (3.10) 当偏压为v 时，对于所有可能发生遂穿的状态求和便可计算总的隧道电流： (3.11) 其中 为费米分布函数， 和 是针尖和样品的 态密度 （dos） 。 考虑 小于测量所要求的能量分辨率， 隧道电流便可简化为： (3.12) 进一步，我们假设矩阵元 在无明显变化，则： (3.13) 这就是著名的bardeen公式。上式表明遂穿电流是由两个子体系态密度的卷积所 决定的，并且样品和针尖电子结构对隧道电流的贡献是对称的。这一公式已经被 证明对大多数 stm 情况都是适用的，被广泛用于隧道电流的实验解释。 3.5 tersoff-hamann (t-h）近似 bardeen方法虽然能够解释大部分的实验结果， 但是该方法计算相对复杂， 需 要计算针尖和样品费米面附近的所有波函数，下面我们介绍的tersoff-hamamm 方法是一种相对简化的stm模拟方法。即假设针尖是具有曲率半径为r的金属， 将针尖波函数简化为s波来处理。这种方法由tersoff 25,26 和hamann提出，因此 成为tersoff-hamamm近似。 由于针尖波函数统一用s波来表示，低偏压时，（3.13）式可以简化为： (3.14) 即隧道电流正比于针尖曲率中心 处费米能级附近的样品的局域态密度。由 上式可以认为恒流模式的stm图像是对样品表面在费米能级附近ldos的描述。 虽然这种模型对针尖模型做了简化，但是利用t-h近似的方法，可以方便的定性 解释stm实验中获得的图像细节， 分析样品电子结构， 近二十年来取得了巨大成 功，例如表面吸附分子、异质结界面、半导体表面等的stm图像都能应用这个近 似方法得到较好的解释。 实际上tersoff-hamann近似只是前面所讲的mba方法的 一个特例。因此，当实验现象不能由t-h近似解释，而必须考虑针尖的实际结构 时，就必须回归到mba方法。 中国科技大学硕士学位论文 第三章 stm 模拟基础 28 图图3.6 s波针尖模型（摘自李斌 2002） 3.6 本文所用 stm 模拟原理21 由上一节可知，虽然tersoff-hamann近似简洁直观，成功解释了很多问题， 但当需要考虑实际体系针尖效应时，tersoff-hamann 近似基本失效，需要使用 chen24发展的mba 方法。由3.4节不难发现，mba方法的精髓在于求解隧穿矩 阵元，将针尖的影响考虑到模拟过程中。隧穿矩阵元m 的表达式为： （3.15） 其中， 为样品波函数， 为针尖波函数，表示局域在真空区域内针尖与样 品之间任意的一个分隔面。由于分隔面的选取的任意性，根据分隔面选取方法 不同就产生了不同的算法。李斌 27 博士将选为一个垂直于针尖取向的水平平 面，利用解析的高斯函数计算隧穿矩阵元m，解释了一系列实验现象。 paz 28等人则采用了一些近似，求解隧穿矩阵元。基本思想是：隧穿过程可 以看作一个电子从样品的等电荷密度面s 通过真空传播到针尖的等电荷密度面 t 的过程。考虑针尖和样品距离较远，两者之间的真空区域可以近似看作一个 方势垒（图3.7 a）。两个曲面的电荷密度都等于 （图3.7 b）。 中国科技大学硕士学位论文 第三章 stm 模拟基础 29 图图3.7 波函数传播示意图（取自paz 2005） 通过这样的近似，样品在 处的波函数 ( )可通过 处的波函数 ( )经真 空格林函数 传播获得： (3.16) 其中， 满足 , ； , 其中 是广义梯度， 可使 只在等电荷 密度面上不为零。 遂穿矩阵元m可表示为： (3.17) 其中， 和 是 和 的傅里叶变换， 是 的傅里 叶变换， 是针尖所在的位置。 这种方法主要的特点是可以通过快速傅立叶变换的方法求出针尖相对与样 品在空间各个位置时的隧穿矩阵元，大大提高了计算效率。但是，这个方法中有 两个不严格的因素。第一，他们将样品波函数直接通过真空格林函数传播，全部 按照s 波形式传播，消弭了样品不同对称性轨道的性质。并且，只考虑等电荷密 度面上的样品波函数的遂穿，丢失了很多其他波函数遂穿的信息。第二， 这 中国科技大学硕士学位论文 第三章 stm 模拟基础 30 一项物理意义和数学意义都不明确。因此，胡振芃21借鉴他们使用的傅立叶变 换思想，重新发展了一种算法。 根据式(3.15)，由于分隔面的选取在针尖和样品间是任意的，为了突出针 尖对样品的影响胡振芃选取一个将针尖包络的曲面进行计算。 这样做的好处是可 以将面积分变为体积分，应用快速傅立叶变换减少计算量。在针尖部分，选取的 曲面也同样选择一个等电荷密度面 。 波函数的获得是直接通过dmol计算软件 的计算获得的。考虑针尖和样品间距（510 ），计算中数值基组的截断半径 选为5.5,这样便能够完全描述实验中的状况。 该算法可以表示成下面的形式： (3.18) 其中， 具有明确的数学意义和物理意义,是由 确定的针尖等电荷密度面上的单位法向量。 得到遂穿矩阵元后，我们可以根据式（3.11）求出隧道电流。为了方便计算我们 可以将（3.11）式变型为： (3.19) 其中 函数类似于(3.11)式中 和 的作用，因为dos相当于将分立的分 子轨道按lorentz展开求和得到的结果，而 函数也是考虑样品和针尖定态（也即 分子轨道）能量本征值的扩展，我们的模拟计算中我们主要应用(3.19)式，我们 也用lorentz函数代替 函数的作用。 中国科技大学硕士学位论文 第四章 cuda实现stm模拟及优化结果 31 第四章第四章 cuda 实现实现 stm 模拟及优化结果模拟及优化结果 我们计算的测试体系采用的针尖分别是一个w原子和h原子，样品分别选择 为金属酞菁分子（copc）（图4.1）以及苯分子（c6h6）。 图图 4.1 模拟所用 copc 分子 计算时选取lda的vwn泛函描述体系的交换相关能，采用极化扩展的双数 值精度基组（dnp）作为电子基组。考虑原子的相对论效应和并减小计算量，在 计算中还采用了软件自带的dspp （density functional semi-core pseudopotentials） 赝势。由于都是单分子，因此我们计算中只考虑一个k点点的贡献，数值积 分网格选用medium。电子自洽计算的能量收敛标准为110-5 hartree。结构优化 时能量的收敛标准为510-5 hartree，能量梯度的收敛标准为510-3 hartree/, 原 子位移的收敛标准为510-3 。 4.1 程序正确性验证 在应用模拟程序之前我们先用dmol计算软件，算出针尖和样品的波函数. 中国科技大学硕士学位论文 第四章 cuda实现stm模拟及优化结果 32 图4.2是用h原子的homo轨道做针尖对c6h6的lumo轨道的模拟结果，图 4.2a是用胡振芃师兄程序的模拟结果，图4.2b是用自己的cuda代码模拟的结果。 可以看出结果完全吻合。 为了让模拟结果更直观，我们用 w 原子做针尖，copc 分子做样品，模拟偏 压选择 5v, 针尖和样品距离为 3 ，如图 4.3 是 copc 分子占据态的模拟结果。 由图 4.2 和图 4.3 可以看出在我们的程序移植是正确的，能够模拟得到正确 的结果。 图图 4.2 源程序与 cuda 程序测试结果 注：注： (a)源程序测试结果 (b) cuda 测试结果 图图 4.3 copc 分子 92 个占据态模拟结果，针尖为 w 原子费米面以上的 4 个空态 (a) (b) 中国科技大学硕士学位论文 第四章 cuda实现stm模拟及优化结果 33 4.2 算法流程及部分代码解析 4.2.1 算法流程 本文首先根据胡振芃师兄发展的stm模拟方法，在原来的基础上加入了电流 的计算，将程序移植到gpu上进行加速测试。具体的算法流程如图4.4所示。 我们首先读出样品波函数，并对样品波函数做傅里叶变换，然后读出针尖波 函数，根据上一章中 ，我们接着计算针尖点和密度 的梯度，继而求c。同样根据上一章中a,b的定义我们求出a,b，并对a，b做傅 里叶变换。这样根据式（3.18）我们得到了遂穿矩阵元m的所有要素，接下来可 以求空间各点的遂穿矩阵元m，并根据式（3.19）求各点的遂穿电流。 我们通过测试，在程序中数据拷贝的时间真正在用于gpu运算时间的10倍， 因此我们尽量将所有的运算都放到gpu中进行运算， 这样做的好处是通过一次拷 贝数据就能在gpu中进行几乎所有的运算， 不会出现在gpu和cpu之间频繁传输 数据的现象。我们分别测量了在gpu中运算的时间，并和c程序进行比较，表3.1 所示。对于mulab,mulas是对傅里叶变换后做单位统一处理，分别对数组元 素乘以一个空间单位体积。 divat是在做傅里叶反变换后做的归一化处理，对 每个元素除以整个变换空间的体积。 中国科技大学硕士学位论文 第四章 cuda实现stm模拟及优化结果 34 读取样品波函数读取样品波函数 对样品波函数做傅里叶对样品波函数做傅里叶 变换变换（ffts） 读取针尖波函数并求电读取针尖波函数并求电 荷密度荷密度 求针尖求针尖 电荷密度的梯电荷密度的梯 度度(grandientden) 求求c（genc） 计算针尖波函数的梯计算针尖波函数的梯 度度（grandientt） 求求a（gena） 求求b（genb） 求遂穿矩阵元求遂穿矩阵元m(mat) 对每个针尖样品波函对每个针尖样品波函 数的电流数的电流i（current） 得到总电流得到总电流i 对对a b 做傅里叶变换做傅里叶变换 fftab 图图4.4 stm模拟算法流程图 注：注：其中蓝框部分移植到了gpu中 中国科技大学硕士学位论文 第四章 cuda实现stm模拟及优化结果 35 4.2.2 代码解析 对应上面的运算流程我们对主要的代码做简要的解析： / operate the samples wave function for(o = 0; o ns; o+) / 读取样品波函数，此部分不能在gpu中操作 fp2 = fopen(fileso, r) ; for(n = 1; n = 5; n+) fgets(pstr,500, fp2); for( i = 0; i = n ; i+) for(j = 0; j = m; j +) for(k = 0; k = l; k +) fscanf(fp2, %f n , if (i!= n wfk + m*j + m*n*i.y = 0.0f; fclose(fp2); / 将波函数用cudamemcpy函数传到device端，对读取的样品波函数应用cufft 做傅里叶变换即流程图中cuffts段。 int size = sizeof(complex)* l* m* n; cudamemcpy(d_wf, wf, size, cudamemcpyhosttodevice); cuffthandle plan; cufftplan3d( cufftexecc2c(plan, (cufftcomplex *)d_wf, (cufftcomplex *)d_as, cufft_forward); cufftdestroy(plan); for(p=0; p nt; p+) / 读取针尖波函数并求电荷密度 fp2 = fopen(filetp, r) ; for(n = 1; n = 5; n+) fgets(pstr,255, fp2); 中国科技大学硕士学位论文 第四章 cuda实现stm模拟及优化结果 36 for (i = 0; i = n; i+ ) for (j = 0; j = m; j+) for(k = 0; k = l; k+) fscanf(fp2, %fn, if(i!=n wfk + m * j + m * n * i.y = 0.0f; denk + m * j + m * n * i.x = t * t; denk + m * j + m * n * i.y = 0.0f; fclose(fp2) ; / 将针尖电荷密度传到device端，并求梯度，即流程图中grandientden，下面用 到的函数都在主函数外独立定义。 cudamemcpy(d_den, den, size, cudamemcpyhosttodevice); gradient(d_den, a, b, c, l, m, n, d_gx, d_gy, d_gz); / 产生c函数，即genc genc(d_cx, d_cy, d_cz, d_gx, d_gy, d_gz, d_den, den_0 ); / 将针尖波函数传入device端，并求梯度，即grandientt cudamemcpy(d_wf, wf, size, cudamemcpyhosttodevice); gradient(d_wf, a, b, c, l, m, n, d_gx, d_gy, d_gz); / 用得到c函数以及波函数梯度得到a和b gena(d_cx, d_cy, d_cz, l, m, n, d_gx, d_gy, d_gz, d_at); genb(d_cx, d_cy, d_cz, l, m, n, d_btx, d_bty, d_btz, d_wf); / 由cufft对a b进行傅里叶变换, 即fftab cufftplan3d( cufftexecc2c(plan, (cufftcomplex *)d_btx, (cufftcomplex *)d_cx, cufft_forward); cufftdestroy(plan); cufftplan3d( cufftexecc2c(plan, (cufftcomplex *)d_at, (cufftcomplex *)d_wf, cufft_forward); 中国科技大学硕士学位论文 第四章 cuda实现stm模拟及优化结果 37 cufftdestroy(plan); cufftplan3d( cufftexecc2c(plan, (cufftcomplex *)d_bty, (cufftcomplex *)d_cy, cufft_forward); cufftdestroy(plan); cufftplan3d( cufftexecc2c(plan, (cufftcomplex *)d_btz, (cufftcomplex *)d_cz, cufft_forward); cufftdestroy(plan); / 由公式（2.18）求遂穿矩阵，即流程图中mat段 computematrix(d_at, d_btx, d_bty, d_btz, d_as, d_wf, hx, hy, hz, l, m, n); cufftplan3d( cufftexecc2c(plan, (cufftcomplex *)d_wf, (cufftcomplex *)d_at, cufft_inverse); cufftdestroy(plan); / 由公式（2.19）计算遂穿电流，即流程图中current段 computecurrent(d_current, d_at, es, et, v); cudamemcpy(current, d_current, size2, cudamemcpydevicetohost); / 各个 kernel 函数 / 梯度函数 _global_ void gradient(complex *d_f, float a, float b, float c, int l, int m, int n, complex *d_gx, complex *d_gy, complex *d_gz) int n1, n2, n3, i, j, k, m1, m2, m3 ; float hx, hy, hz ; hx = a/l ; hy = b/m ; hz = c/n ; int id = griddim.x * blockdim.x * blockidx.y + blockdim.x * blockidx.x + threadidx.x; n1 = blockidx.y - 1; n2 = blockidx.x - 1; n3 = threadidx.x - 1; if(n1 = -1) n1 = n - 1; 中国科技大学硕士学位论文 第四章 cuda实现stm模拟及优化结果 38 if(n2 = -1) n2 = m - 1; if(n3 = -1) n3 = l - 1; m1 = blockidx.y + 1; m2 = blockidx.x + 1; m3 = threadidx.x + 1; if(m1 = n) m1 = 0; if(m2 = m) m2 = 0; if(m3 = l) m3 = 0; d_gxid.x = (d_fm3 + griddim.x * blockidx.x + griddim.x * blockdim.x * blockidx.y.x - d_fn3 + blockdim.x * blockidx.x + griddim.x * blockdim.x * blockidx.y.x) / (2 * hx); d_gyid.x = (d_fthreadidx.x + blockdim.x * m2 + griddim.x * blockdim.x * blockidx.y.x - d_fthreadidx.x + blockdim.x * n2 + griddim.x * blockdim.x * blockidx.y.x) / (2 * hy); d_gzid.x = (d_f griddim.x * blockdim.x * m1 + blockdim.x * blockidx.x + threadidx.x.x - d_f griddim.x * blockdim.x * n1 + blockdim.x * blockidx.x + threadidx.x.x) / (2 * hz)