计算材料学计算材料学 晶体结构计算、分子动力学晶体结构计算、分子动力学主讲人：黄远主讲人：黄远天津大学材料科学与工程学院天津大学材料科学与工程学院2021/6/161主要内容主要内容晶体学中的对称操作元素：晶体学中的对称操作元素：旋转轴、倒反中心、镜面、反轴、映轴、螺旋旋转轴、倒反中心、镜面、反轴、映轴、螺旋轴和滑移面轴和滑移面晶体学点群晶体学点群，晶系和点阵型式晶系和点阵型式空间群及其应用：空间群及其应用： 空间群符号，等效点系，分数坐标，不对称单空间群符号，等效点系，分数坐标，不对称单位位预备知知识：对称性称性-从点从点阵到空到空间群群2021/6/162一、晶体学中的对称操作元素一、晶体学中的对称操作元素分子和晶体都是对称图像。对称图像是一个能经过分子和晶体都是对称图像。对称图像是一个能经过不改变不改变其中任何两点间距离的操作其中任何两点间距离的操作后后复原复原的图像，这样的操作称的图像，这样的操作称为为对称操作对称操作。在操作中保持空间中至少一个点不动的对称操作称为在操作中保持空间中至少一个点不动的对称操作称为点对点对称操作称操作（点式操作点式操作），如简单旋转和镜像转动，如简单旋转和镜像转动( (反映和倒反映和倒反反););使空间中所有点都运动的对称操作称为使空间中所有点都运动的对称操作称为非点式操作非点式操作，如平移，螺旋转动和滑移反映。如平移，螺旋转动和滑移反映。 2021/6/163对称操作对称操作： 一个物体运动或变换，使得变换后的物体与变换一个物体运动或变换，使得变换后的物体与变换前不可区分（复原，重合）。前不可区分（复原，重合）。对称元素对称元素：在对称操作中保持不变的几何图型：点、轴或面。：在对称操作中保持不变的几何图型：点、轴或面。2021/6/1641 1、点式操作、点式操作（1 1）全同操作）全同操作(1) 全同操作全同操作(Identity)，符号表示，符号表示为为1(E)，对应对应于物体不于物体不动动的的对对称操作，相称操作，相对应对应的的变换变换矩矩阵为单阵为单位矩位矩阵阵。矩矩阵表示表示注意：注意：符号表示符号表示为国国际符号也称符号也称为赫赫尔曼曼- -毛古因毛古因Hermann-Hermann-MauguinMauguin符号，括号内符号，括号内为熊夫利斯熊夫利斯Schnflies Schnflies 符号。符号。2021/6/165（2 2）旋转轴）旋转轴 旋转轴旋转轴(旋转轴旋转轴)：绕某轴反时针旋转绕某轴反时针旋转 =360/n度，度，n称称为旋转轴的次数为旋转轴的次数(或重数或重数)，符号为符号为n (Cn)。其变换矩阵为：。其变换矩阵为：2021/6/166变换矩阵的获得变换矩阵的获得2021/6/167晶体中的旋转轴限制晶体中的旋转轴限制1 1、平移对称性对旋转轴的次数、平移对称性对旋转轴的次数n n有很大的限制，证明在晶体学有很大的限制，证明在晶体学中只能出中只能出n=1,2,3,4,6n=1,2,3,4,6的旋转轴。的旋转轴。2 2、写出沿三个坐标轴、写出沿三个坐标轴X X，Y Y和和Z Z的的4 4次旋转轴的表示矩阵次旋转轴的表示矩阵。2021/6/168（3 3）倒反中心）倒反中心( (Inversion Inversion centercenter) ) 倒反中心：倒反中心：也称为反演中心或对称中心也称为反演中心或对称中心( (Center of Center of symmetrysymmetry) )，它的操作是通过一个点的倒反，它的操作是通过一个点的倒反( (反演反演) )，使空间，使空间点的每一个位置由坐标为点的每一个位置由坐标为(x(x、y y， z) z)变换到变换到(- x(- x， - y - y， - z)- z)。符号为。符号为1(i)1(i)，变换矩阵为，变换矩阵为2021/6/169（4 4）反映面）反映面-镜面镜面 反映面，反映面，也称镜面也称镜面，反映操作是从空间某一点向反映面反映操作是从空间某一点向反映面引垂线，并延长该垂线到反映面的另一侧，在延长线上取引垂线，并延长该垂线到反映面的另一侧，在延长线上取一点，使其到反映面的距离等于原来点到反映面的距离。一点，使其到反映面的距离等于原来点到反映面的距离。符号为符号为m ( )。 为了表示反映面的方向，可以在其符号后面标以该面的为了表示反映面的方向，可以在其符号后面标以该面的法线。如法线为法线。如法线为010的反映面，可记为的反映面，可记为m 010。m010(x、y，z)=(x，-y，z)2021/6/1610镜面类型和矩阵表示镜面类型和矩阵表示关于对称平面（或镜面）关于对称平面（或镜面）的反的反映映，可以平行于，可以平行于(vertical(vertical，v v)或或 垂直于垂直于(horizontal(horizontal， h h)主轴。主轴。 在二个在二个C C2 2轴之间角平分线的一个垂直平面叫作双面镜面，轴之间角平分线的一个垂直平面叫作双面镜面，d d （dihedralplanedihedralplane）。）。通通过yz面的反映面的反映。2021/6/1611（5 5）旋转倒反轴）旋转倒反轴- -反轴反轴旋转倒反轴，简称反轴旋转倒反轴，简称反轴 (Axis of inversion ，Rotoinversion axis)，其对称操作是先进行旋转操作其对称操作是先进行旋转操作(n)后立刻再进行倒反操后立刻再进行倒反操作作，这样的复合操作称为记为这样的复合操作称为记为组合成合成这种复合操作的每一个操作本身不一定是种复合操作的每一个操作本身不一定是对称称操作。其矩操作。其矩阵表示表示为：2021/6/1612（6）旋转反映轴）旋转反映轴-映轴映轴旋转反映旋转反映 S Sn n，包括绕对称轴的逆时针旋转，包括绕对称轴的逆时针旋转360/n360/n，接着，接着作垂直反射。作垂直反射。符号为符号为(Sn)，设对称轴沿，设对称轴沿001方向，其矩阵表示为：方向，其矩阵表示为：2021/6/1613 （7 7）反轴和映轴间的对应关系）反轴和映轴间的对应关系旋转倒反轴和旋转反映轴之间存在简单的一一对应关系，旋旋转倒反轴和旋转反映轴之间存在简单的一一对应关系，旋转角度为转角度为 的反轴和旋转角为的反轴和旋转角为( (-p-p) )的映轴是等价的对称轴。的映轴是等价的对称轴。 这一关系也很容易从他们的表示矩阵看出。所以这一关系也很容易从他们的表示矩阵看出。所以1 1次，次， 2 2次，次， 3 3次，次， 4 4次和次和6 6次反轴分别等价于次反轴分别等价于2 2次，次， 1 1次，次， 6 6次，次， 4 4次和次和3 3次映轴。次映轴。2021/6/16142 2、非点式对称操作、非点式对称操作非点式对称操作：是由点式操作与平移操作复合后形成的新非点式对称操作：是由点式操作与平移操作复合后形成的新的对称操作的对称操作- -、平移和旋转复合形成能导出螺旋旋转，平移和旋转复合形成能导出螺旋旋转，、平移和反映复合能导出滑移反映。平移和反映复合能导出滑移反映。2021/6/1615（1 1）螺旋轴）螺旋轴螺旋轴：螺旋轴：先绕轴进行逆时针方向先绕轴进行逆时针方向360/n度的旋转，接着作平度的旋转，接着作平行于该轴的平移行于该轴的平移，平移量为平移量为(p/n) t，这里这里t是平行于转轴方向是平行于转轴方向的最短的晶格平移矢量的最短的晶格平移矢量， n称为螺旋轴的次数，称为螺旋轴的次数， ( (n可以取可以取值值2,3,4,6)，而而p只取小于只取小于n的整数。所以可以有以下的整数。所以可以有以下11种螺旋种螺旋轴：轴：21，31，32，41，42，43，61，62，63，64，65。2021/6/1616螺旋轴螺旋轴 2 21 1，3 31 1 ，3 32 2 ，6 63 32021/6/1617螺旋轴螺旋轴4 41 1，4 42 2 ，4 43 341和和43彼此对映。彼此对映。当其中之一是左手当其中之一是左手螺旋时，另一个为螺旋时，另一个为右手螺旋。右手螺旋。2021/6/1618螺旋轴螺旋轴6 61 1，6 62 2，6 63 3，6 64 42021/6/1619（2 2）滑移面）滑移面滑移反映面，滑移反映面，(滑移面滑移面)简称滑移面简称滑移面，其对称操作是沿滑移面其对称操作是沿滑移面进行镜面反映操作，然后接着进行与平行于滑移面的一个方进行镜面反映操作，然后接着进行与平行于滑移面的一个方向的平移。向的平移。沿滑移面的平移量等于该方向点阵平移周期的一半。沿滑移面的平移量等于该方向点阵平移周期的一半。2021/6/1620滑移反射滑移反射不对称单位先经滑移面反映，然后沿平行与滑移面的方向平移。不对称单位先经滑移面反映，然后沿平行与滑移面的方向平移。 2021/6/1621滑移面分类滑移面分类轴向滑移面：沿晶轴轴向滑移面：沿晶轴(a、b， c)方向滑移方向滑移; ;对角滑移面：沿晶胞面对角线或体对角线方向滑移，平移分对角滑移面：沿晶胞面对角线或体对角线方向滑移，平移分量为对角线一半量为对角线一半; ;金刚石滑移面：沿晶胞面对角线或体对角线方向滑移，平移金刚石滑移面：沿晶胞面对角线或体对角线方向滑移，平移分量对角线分量对角线1/4的对角滑移面。的对角滑移面。只有在体心或面心点阵中出现只有在体心或面心点阵中出现。2021/6/1622镜面和滑移面镜面和滑移面镜面或滑移面的符号。面或滑移面的符号。 (在左在左边： 沿沿镜面的面的边缘看。看。 在右在右边是沿垂直于是沿垂直于镜面的方向面的方向观看。看。 箭箭头表示平移方向。表示平移方向。a a， b b， c c是平行于是平行于单胞胞边的滑移。的滑移。 n n是是对角滑移，在两个角滑移，在两个方向都滑移方向都滑移单胞胞长度度的一半。的一半。 d d是是类似似n n的的对角滑移，角滑移，但但这里在每个方向移里在每个方向移动单胞胞边长的的1/41/4。2021/6/16231、群的定义、群的定义 假设假设G G是由一些元素组成的集合，即是由一些元素组成的集合，即G= G= ， g g，。 在在G G中定义了一种规则中定义了一种规则( (操作、运算，群的乘法操作、运算，群的乘法) )。 如果如果G G对这对这种合成规则满足以下四个条件：种合成规则满足以下四个条件： a)a)封闭性：封闭性： G G中任意两个元素的乘积仍然属于中任意两个元素的乘积仍然属于G G。 b) b)结合律：结合律：二、晶体学中的群论二、晶体学中的群论2021/6/1624c)单位元素。单位元素。集合集合G中存在一个单位元素中存在一个单位元素e，对任意元素，对任意元素，有有d)可逆性。可逆性。对任意元素对任意元素，存在逆元素，存在逆元素，使，使则称集合则称集合G为一个群。为一个群。2021/6/16252 2、晶体学点群、晶体学点群定义定义晶体中满足群的性质定义的点对称操作的集合称作晶体学点群。晶体中满足群的性质定义的点对称操作的集合称作晶体学点群。点对称操作的共同特征是进行操作后物体中至少有一个点是不动点对称操作的共同特征是进行操作后物体中至少有一个点是不动的。的。如果把点对称操作元素按所有可能组合起来，则一共可以得出如果把点对称操作元素按所有可能组合起来，则一共可以得出32种不同的组合方式，称为种不同的组合方式，称为32个晶体学点群。个晶体学点群。2021/6/1626（2 2）点群的）点群的SchnfliesSchnflies符号符号C Cn n: : 具有一个具有一个n n次旋转轴的点群。次旋转轴的点群。C Cnhnh: : 具有一个具有一个n n次旋转轴和一个垂直于该轴的镜面的点群。次旋转轴和一个垂直于该轴的镜面的点群。C Cnvnv: : 具有一个具有一个n n次旋转轴和次旋转轴和n n个通过该轴的镜面的点群。个通过该轴的镜面的点群。D Dn n: : 具有一个具有一个n n次旋转主轴和次旋转主轴和n n个垂直该轴的二次轴的点群。个垂直该轴的二次轴的点群。S Sn n: :具有一个具有一个n n次反轴的点群。次反轴的点群。T:T:具有具有4 4个个3 3次轴和次轴和4 4个个2 2次轴的正四面体点群。次轴的正四面体点群。O:O:具有具有3 3个个4 4次轴次轴,4,4个个3 3次轴和次轴和6 6个个2 2次轴的八面体点群。次轴的八面体点群。2021/6/1627（3 3）3232种点群的表示符号种点群的表示符号C1,C2,C3,C4,C6;1,2,3,4,6-旋旋转轴(C=cyclic)；C1h=Cs,C2h,C3h,C4h,C6h;m,2/m,3/m,4/m,6/m-旋旋转轴加上垂直于加上垂直于该轴的的对称平面；称平面；C2v,C3v,C4v,C6v;mm2,3m,4mm,6mm-旋旋转轴加通加通过该轴的的镜面；面；S2=Ci,S4,S6=C3d;-旋旋转反演反演轴；2021/6/1628D D2 2,D,D3 3,D,D4 4,D,D6 6; ; 222,32,422,622222,32,422,622 - -旋旋转轴(n)(n)加加n n个垂直于个垂直于该轴的二次的二次轴； D D2h2h,D,D3h3h,D,D4h4h,D,D6h6h; ;mmm,3/mm,4/mm,6/mmmmmm,3/mm,4/mm,6/mmm - -旋旋转轴(n)(n)加加n n个垂直于个垂直于该轴的二次的二次轴和和镜面；面；D D2d2d,D,D3d3d; ; -42m,-3m -42m,-3m - -D D群附加群附加对角角竖直平面：直平面： T, TT, Th h, O, T, O, Td d, O, Oh h; ; 23,m3,432,-43m,m3m23,m3,432,-43m,m3m - -立方体群立方体群(T=tetrahedral(T=tetrahedral， O=octahedral) O=octahedral)2021/6/1629晶体点群的Schnflies和国际符号2021/6/1630（4 4）各点群的对称元素和所在方向）各点群的对称元素和所在方向晶系晶系第一位第一位第二位第二位第三位第三位点群点群可能对称可能对称元素元素方向方向可能对称可能对称元素元素方向方向可能对称可能对称元素元素方向方向三斜三斜1, 1任意任意无无无无1， 1单斜单斜2,m,2/mY无无无无2,m,2/m正交正交2，mX2，mY2，mZ222,mm2,mmm四方四方4, 4，4/m Z无，无，2，mX无，无，2，m底对底对角线角线4, 4，4/m，422，4mm, 42m,4/mmm三方三方3, 3Z无，无，2，mX无无3, 3,32,3m, 3m六方六方6, 6,6/mZ无，无，2，mX无，无，2，m底对底对角线角线6, 6,6/m,622,6mm, 62m,6/mmm立方立方2,m,4, 4X3, 3体对体对角线角线无，无，2，m面对面对角线角线23,m3,432, 43m,m3m2021/6/1631（5 5）空间点阵型式）空间点阵型式-布拉伐点阵布拉伐点阵空间点阵按点群对称性和带心的模式一共可以产生空间点阵按点群对称性和带心的模式一共可以产生1414种型式，种型式，称为称为1414种布拉伐点阵或布拉伐点阵。布拉伐点阵表示出所属空种布拉伐点阵或布拉伐点阵。布拉伐点阵表示出所属空间群的平移子群。间群的平移子群。 Bravais Bravais点阵点阵 描述点阵的纯平移对称。描述点阵的纯平移对称。 实质上通过指定实质上通过指定BravaisBravais点阵，指定了单胞点阵，指定了单胞( (晶系晶系) )的形状和带心的型式。的形状和带心的型式。 2021/6/16321414种空间点阵型式示意图种空间点阵型式示意图(14(14个个BravaisBravais点阵点阵) )2021/6/16331414种种BravaisBravais点阵的点阵参数点阵的点阵参数2021/6/16343 3、空间群、空间群(Space Group)(Space Group)空间群是点阵、平移群空间群是点阵、平移群( (滑移面和螺旋轴滑移面和螺旋轴) )和点群的组合。和点群的组合。 230230个空间群是由个空间群是由1414个个BravaisBravais点阵与点阵与3232个晶体点群系统组合个晶体点群系统组合而成。而成。参见： http:/asdp.bio.bnl.gov/asda/Libraries/sgtable.html2021/6/1635（1）从点群到空间群）从点群到空间群7个晶系个晶系旋旋转，反射，反演，反射，反演平移平移螺旋螺旋轴，滑移面，滑移面3232个点群个点群14种种Bravais格子格子230个空个空间群群（按照晶胞的特征按照晶胞的特征对称元素分称元素分类）2021/6/1636（2）空间群分布）空间群分布三斜晶系：三斜晶系：2个；单斜晶系：个；单斜晶系：13个个正交晶系：正交晶系：59个；个；三方晶系：三方晶系：25四方晶系：四方晶系：68个；六方晶系：个；六方晶系：27个个立方晶系：立方晶系：36个。个。有对称中心有对称中心90个，无对称中心个，无对称中心140个。个。73个个symmorphic (点式点式) ，157个个non-symmorphic。2021/6/1637（3 3）了解）了解Herman-MauguinHerman-Mauguin空间群符号空间群符号空间群是经常用简略空间群是经常用简略Herman-Mauguin符号符号(即即Pnma、I4/mmm等等)来指定。来指定。在简略符号中包含能产生所有其余在简略符号中包含能产生所有其余对称元素所必需的最少对称元素。对称元素所必需的最少对称元素。从简略从简略H-M符号，我们可以确定晶系、符号，我们可以确定晶系、Bravais点阵、点点阵、点群和某些对称元素的存在和取向（反之亦然群和某些对称元素的存在和取向（反之亦然）。2021/6/1638（4 4）空间群符号）空间群符号L LS S1 1S S2 2S S3 3运用以下规则，可以从对称元素获得运用以下规则，可以从对称元素获得H-M空间群符号。空间群符号。 .第一字母第一字母(L)是点阵描述符号，指明点阵带心类型：是点阵描述符号，指明点阵带心类型：P， I， F， C， A， B。 .其于三个符号其于三个符号(S1S2S3)表示在特定方向（对每种晶系分别规定）表示在特定方向（对每种晶系分别规定）上的对称元素。上的对称元素。 .如果没有二义性可能，常用符号的省略形式如果没有二义性可能，常用符号的省略形式 (如如Pm，而不用，而不用写成写成P1m1)。 * 由于不同的晶轴选择和标记，同一个空间群可能有几种不同由于不同的晶轴选择和标记，同一个空间群可能有几种不同的符号。如的符号。如P21/c,如滑移面选为在如滑移面选为在a方向，符号为方向，符号为P21/a;如滑移如滑移面选为对角滑移，符号为面选为对角滑移，符号为P21/n。2021/6/1639（5 5）从空间群符号辨认晶系和点群）从空间群符号辨认晶系和点群1）立方）立方第第2个个对称符号：对称符号： 3 或或 3 (如：如： Ia3, Pm3m, Fd3m) 2）四方）四方第第1个个对称符号：对称符号： 4, 4 , 41, 42 或或 43 (如：如： P41212, I4/m, P4/mcc) 3）六方）六方第第1个个对称符号：对称符号： 6, 6 , 61, 62, 63, 64 或或 65 (如：如： P6mm, P63/mcm) 4）三方）三方第第1个个对称符号：对称符号： 3, 3 ，31 或或 32 (如：如： P31m, R3, R3c, P312) 5）正交）正交点阵符号后的全部点阵符号后的全部三个符号三个符号是镜面，滑移面，是镜面，滑移面，2次旋转次旋转轴或轴或2次螺旋次螺旋(即即Pnma, Cmc21, Pnc2）6）单斜）单斜点阵符号后有唯一的镜面、滑移面、点阵符号后有唯一的镜面、滑移面、2次旋转或者螺旋次旋转或者螺旋轴，或者轴轴，或者轴/平面符号平面符号(即即Cc、P2、P21/n)7）三斜）三斜点阵符号后是点阵符号后是1或或(- 1)2021/6/1640从空间群符号从空间群符号确定确定点群点群点群可以从简略点群可以从简略H-M符号通过下列变换得出：符号通过下列变换得出：1）把所有滑移面全部转换成镜面；）把所有滑移面全部转换成镜面；2）把所有螺旋轴全部转换成旋转轴。）把所有螺旋轴全部转换成旋转轴。例如例如：空间群空间群=Pnma点群点群=mmm空间群空间群=I 4c2点群点群= 4m2空间群空间群=P42/n点群点群=4/m2021/6/1641晶系对称方向第一第二第三三斜无单斜b 010正交a 100b 010c 001四方c 001a 100/010a+b 110六方c 001a 100/0102a+b 120三方 (R)a+b+c 111a-b 11 0立方a 100/010/001a+b+c 111a+b 1102021/6/1642（6 6）X-X-射线结晶学国际表射线结晶学国际表 ( () )提供的信息的是：提供的信息的是：1 1）空间群的国际符号为）空间群的国际符号为2 2）SchoenfliesSchoenflies符号符号3 3）晶系）晶系4 4）晶类）晶类5）一般等效点图：）一般等效点图：单胞的投影，包含所有等效点位置。单胞的投影，包含所有等效点位置。“+”表示表示z0，“-“表示表示ztmax结束结束YN2021/6/161021.2NEV系综分子动力学过程系综分子动力学过程有确定的粒子数有确定的粒子数N、体积、体积V和能量和能量E的系综的系综-NEV系综系综1、相互作用势及作用在粒子上的力、相互作用势及作用在粒子上的力2021/6/16103分子动力学单位分子动力学单位(CGS)在采用在采用Lennard-Jones势后，所有的量应该采取标势后，所有的量应该采取标度后的形式进行度后的形式进行长度的标度单位为长度的标度单位为 能量的标度单位为能量的标度单位为 标度后的标度后的Lennard-Jones势能为势能为2、2021/6/16104时间的标度单位为时间的标度单位为于是有于是有温度的标度单位为温度的标度单位为于是有于是有分子密度为分子密度为压强为压强为2021/6/161053、运动方程及求解、运动方程及求解NEV系综中粒子遵循的牛顿运动方程组如下：系综中粒子遵循的牛顿运动方程组如下：采用差分方法推导后的数值求解形式如下：采用差分方法推导后的数值求解形式如下：Verlet算法形式的改进算法形式的改进2021/6/161064、势能作用的截断距离、势能作用的截断距离基本元胞中的一个粒子只同基本元胞中的另外基本元胞中的一个粒子只同基本元胞中的另外N-1个个粒子中每个粒子或其最邻近影像发生相互作用，实际粒子中每个粒子或其最邻近影像发生相互作用，实际上这是根据条件上这是根据条件来截断位势（意味着大于来截断位势（意味着大于L/2的粒子的相互作用小得可的粒子的相互作用小得可以忽略）以忽略）2021/6/161075、边界条件、边界条件（1）自由边界）自由边界-大型的自由分子模拟大型的自由分子模拟（2）固定边界）固定边界-常用于点缺陷等性质的研究常用于点缺陷等性质的研究（3）柔性边界）柔性边界-用于模拟缺陷延伸的情况，比如位错运动用于模拟缺陷延伸的情况，比如位错运动（4）周期性边界）周期性边界-使用广泛，消除表面效应使用广泛，消除表面效应处理周期性边界条件，如果有一个粒子处理周期性边界条件，如果有一个粒子穿过基本元胞的一个表面，那么这个粒子就穿过基本元胞的一个表面，那么这个粒子就穿过对面的墙重新进入元胞，速度不变。穿过对面的墙重新进入元胞，速度不变。A（x）=A（x+nL），），n=（n1，n2，n3）2021/6/161086、标度因子（、标度因子（ScalingFactor）对于对于NEV系综，其有确定的能量，但不可能给出精确的初始系综，其有确定的能量，但不可能给出精确的初始条件，这时需要先给出一个合理的初始条件，然后在模拟的条件，这时需要先给出一个合理的初始条件，然后在模拟的过程中逐渐调节系统能量达到给定值。具体的步骤如下：过程中逐渐调节系统能量达到给定值。具体的步骤如下：（1）首先，将运动方程组求解出若干步的结果）首先，将运动方程组求解出若干步的结果（2）然后，计算出动能和位能）然后，计算出动能和位能（3）假如总能量不等于给定恒定值，则通过调整速度来实现）假如总能量不等于给定恒定值，则通过调整速度来实现总能量不变的目的总能量不变的目的（4）具体的做法是将速度乘以一个标度因子）具体的做法是将速度乘以一个标度因子 ，即：，即：2021/6/161097、一些力学量的计算、一些力学量的计算2021/6/161101.3分子力场分子力场量子化学计算分子结构和原子、分子间相互作用比较准确，量子化学计算分子结构和原子、分子间相互作用比较准确，但是很慢；而采用分子力场计算就会很快，因为分子力场并不但是很慢；而采用分子力场计算就会很快，因为分子力场并不计算电子相互作用，它是对分子结构的一种简化模型，所以计计算电子相互作用，它是对分子结构的一种简化模型，所以计算很快。在这个模型中，它把组成分子的原子看成是由弹簧连算很快。在这个模型中，它把组成分子的原子看成是由弹簧连接起来的球，然后用简单的数学函数来描述球与球之间的相互接起来的球，然后用简单的数学函数来描述球与球之间的相互作用。作用。力场用简单的数学函数描述原子间作用，称为分子力场，又叫力场用简单的数学函数描述原子间作用，称为分子力场，又叫分子力学力场。采用分子力场的分子模拟称为经典分子模拟。分子力学力场。采用分子力场的分子模拟称为经典分子模拟。一个力场通常包括三个部分：原子类型，势函数，和力场一个力场通常包括三个部分：原子类型，势函数，和力场参数。参数。2021/6/16111例如：氢分子，看做有弹簧链接的两个球的话，可以用胡克例如：氢分子，看做有弹簧链接的两个球的话，可以用胡克定律描述两个氢原子间的能量：定律描述两个氢原子间的能量：E=k*(b-b0)2其中，其中，b表示两氢原子间距离，表示两氢原子间距离，b0表示平衡时原子间距，表示平衡时原子间距，k为为键能系数，键能系数，b0和和K称为力场参数。称为力场参数。更复杂一点可以用四次方表达：更复杂一点可以用四次方表达：E=K1*(b-b0)2+K2*(b-b0)3+K3*(b-b0)4，更多的参数可以获得对成键分子的更精确的描述。这是描更多的参数可以获得对成键分子的更精确的描述。这是描述成键作用，不成键的原子间的相互作用则采用述成键作用，不成键的原子间的相互作用则采用Legendre-Jones函数，或者函数，或者Bukingham函数描述。函数描述。2021/6/16112分子力场有时被称为势函数。以下是一般分子力场势函数包括分子力场有时被称为势函数。以下是一般分子力场势函数包括的几个部分：的几个部分：（1）描述分子内成键作用的项）描述分子内成键作用的项a、键伸缩能：构成分子的各个化学键在键轴方向上的伸缩、键伸缩能：构成分子的各个化学键在键轴方向上的伸缩运动所引起的能量变化运动所引起的能量变化b、键角弯曲能：键角变化引起的分子能量变化、键角弯曲能：键角变化引起的分子能量变化c、二面角扭曲能：单键旋转引起分子骨架扭曲所产生的能、二面角扭曲能：单键旋转引起分子骨架扭曲所产生的能量变化量变化d、交叉能量项：上述作用之间耦合引起的能量变化、交叉能量项：上述作用之间耦合引起的能量变化（2）描述分子间作用的项）描述分子间作用的项非键相互作用：包括范德华力、静电相互作用等与能量有非键相互作用：包括范德华力、静电相互作用等与能量有关的非键相互作用关的非键相互作用2021/6/16113COMPASS力场的模型体系力场的模型体系COMPASS力场是第一个把以往分别处理的有机分子体系的力场力场是第一个把以往分别处理的有机分子体系的力场与无机分子体系的力场统一的分子力场。与无机分子体系的力场统一的分子力场。COMPASS力场能够模拟小分子与高分子，一些金属离子、金属力场能够模拟小分子与高分子，一些金属离子、金属氧化物与金属。在处理有机与无机体系时，采用分类别处理的方式，氧化物与金属。在处理有机与无机体系时，采用分类别处理的方式，不同的体系采用不同的模型，即使对于两类体系的混合，仍然能够不同的体系采用不同的模型，即使对于两类体系的混合，仍然能够采用合理的模型描述。采用合理的模型描述。在在COMPASS力场中，有下列几种不同模型体系，描述不同的分力场中，有下列几种不同模型体系，描述不同的分子体系。子体系。2021/6/16114(1)共价键模型共价键模型对于所有的有机分子和无机共价键分子体系，包括小分子和对于所有的有机分子和无机共价键分子体系，包括小分子和高分子，高分子，COMPASS采用采用CFF分子力场的基本模式。分子力场的基本模式。函数形式由两部分组成，键合项和非键合项：函数形式由两部分组成，键合项和非键合项：a、键合项包括对角和非对角的交叉偶合项式、键合项包括对角和非对角的交叉偶合项式(1)，分别为，分别为键键伸缩能伸缩能Eb，键角弯曲能键角弯曲能Eq，键扭转能键扭转能Ef，键角面外弯曲能键角面外弯曲能E以及以及它们之间的它们之间的相互偶合能相互偶合能Ecross。b、非键合项包括、非键合项包括范德华能范德华能Eij，见式见式(2)和和库仑能库仑能Eelec，见见式式(3)，分别用来计算距离两个或两个以上的原子对或不同分子，分别用来计算距离两个或两个以上的原子对或不同分子链上的原子对的相互作用力。链上的原子对的相互作用力。2021/6/16115 （1）2021/6/16116(2) 范德华相互作用采用范德华相互作用采用Lennard-Jones-9-6Lennard-Jones-9-6函数，对于不同函数，对于不同种原子对的之间的参数种原子对的之间的参数r ri,ji,j， i,ji,j可以在同种原子对参数可以在同种原子对参数 、r r0 0的基础上，采用六次方平均的方法计算不同种原子对的参的基础上，采用六次方平均的方法计算不同种原子对的参数：数：（2）2021/6/16117（3）静电相互作用通过原子剩余电荷来描述。它采用键增量静电相互作用通过原子剩余电荷来描述。它采用键增量d dijij表示原表示原子子j j对原子对原子i i的剩余电荷。在键增量的剩余电荷。在键增量ijij中，中，i i表示电荷受体，表示电荷受体，j j表示表示电荷给体，例如电荷给体，例如d dCHCH=-0.053,=-0.053,表示碳氢键中氢原子对碳原子的剩余电表示碳氢键中氢原子对碳原子的剩余电荷为荷为-0.053-0.053，电荷键增量，电荷键增量d dijij是通过从头计算静电势能得到的。原子是通过从头计算静电势能得到的。原子i i的剩余电荷是所有与其成键的电荷键增量的剩余电荷是所有与其成键的电荷键增量d dijij的总和。的总和。 2021/6/16118（2）离子模型）离子模型COMPASS采用了一个简单的离子模型，它包括：采用了一个简单的离子模型，它包括：静电项式静电项式(2)和范德华项式和范德华项式(3)。在此模型中，每一个原子作为一个非键合的粒子在此模型中，每一个原子作为一个非键合的粒子，既离子之间，既离子之间没有特殊的拓扑连接方式与特殊的距离。没有特殊的拓扑连接方式与特殊的距离。而粒子的堆砌结构是倚赖与不同离子之间的的强的静电引力和而粒子的堆砌结构是倚赖与不同离子之间的的强的静电引力和范德华斥力决定的。范德华斥力决定的。在该模型中，不采用键端电荷差确定原子电荷在该模型中，不采用键端电荷差确定原子电荷的方法。的方法。2021/6/16119（3）准离子模型）准离子模型所谓所谓“准离子模型准离子模型”，即处理的原子体系介乎于纯离子与，即处理的原子体系介乎于纯离子与共价键体系之间。共价键体系之间。这里仍然采用了在离子模型中使用的方式，把每个原子描这里仍然采用了在离子模型中使用的方式，把每个原子描述成非键合的粒子。述成非键合的粒子。该模型的能量表达是由两部分组成：直接接触原子间的能该模型的能量表达是由两部分组成：直接接触原子间的能量和不接触原子间的能量：量和不接触原子间的能量：对于不接触原子间的能量，仍用式（对于不接触原子间的能量，仍用式（2）与式（）与式（3）。对于）。对于直接接触原子间的能量，采用摩斯色散函数与静电相互作用直接接触原子间的能量，采用摩斯色散函数与静电相互作用项的加和，见式（项的加和，见式（4）。）。2021/6/16120fs是一个连续转换函数，将短程的摩斯函数转变为长程的色是一个连续转换函数，将短程的摩斯函数转变为长程的色散函数散函数（4）2021/6/16121（4）金属模型金属模型用分子力场的方法准确地描述金属与金属合金是一个很用分子力场的方法准确地描述金属与金属合金是一个很难的课题。难的课题。COMPASS仍使用比较粗糙的模型来描述纯金属仍使用比较粗糙的模型来描述纯金属的本体。所有参数来自对体系的晶格能与结构的能量优化的的本体。所有参数来自对体系的晶格能与结构的能量优化的结果。结果。金属原子仍然处理成非键合的粒子，其相互作用通过金属原子仍然处理成非键合的粒子，其相互作用通过LJ9-6函数计算。其参数的贡献既包括接触原子间的，也包括函数计算。其参数的贡献既包括接触原子间的，也包括远程原子间的。远程原子间的。也就是说，这两种相互作用采用同一参数、也就是说，这两种相互作用采用同一参数、同一函数计算，这种近似方法的优劣需要实践的考察。同一函数计算，这种近似方法的优劣需要实践的考察。2021/6/16122C=CMOLECULARDYNAMICS-分子动力学分子动力学CCMICROCANONICALENSEMBLE-微正则系统微正则系统CCAPPLICATIONTOARGON.THELENNARD-JONESPOTENTIALISTRUNCATEDCATRCOFFANDNOTSMOOTHLYCONTINUEDTOZERO.INITIALLYTHECNPARTPARTICLESAREPLACEDONANFCCLATTICE.THEVELOCITYSCAREDRANFROMABOLTZMANDISTRIBUTIONWITHTEMPERATURETREF.C=分子动力学部分分子动力学部分1.3分子动力学程序讲解分子动力学程序讲解该部分主要为程序的介绍该部分主要为程序的介绍及说（申）明及说（申）明2021/6/16123CCINPUTPARAMETERSAREASFOLLOWSCCNPARTNUMBEROFPARTICLES(MUSTBEAMULTIPLEOF4)粒子数粒子数CSIDESIDELENGTHOFTHECUBICALBOXINSIGMAUINTS-立方盒子边长立方盒子边长CTREFREDUCEDTEMPERATURE-约化温度约化温度CRCOFFCUTOFFOFTHEPOTENTIALINSIGMAUINTS-势能截断距离势能截断距离CHBASICTIMESTEP-基本时间步基本时间步CIREPVELOCITIESSCALINGEVERYIREPTHTIMESTEP-速度标度因子速度标度因子CISTOPSTOPOFSCALINGOFTHEVELOCITIESATISTOPCTIMEMXNUMBEROFINTEGRATIONSTEPS-停止标度标志停止标度标志CISEEDSEEDFORTHERANDOMNUMBERGENERATOR-随机数产生种子随机数产生种子CCVERSION1.0ASOFAUGUST1985CCDIETERW.HEERMANN该部分主要为程序中各参数该部分主要为程序中各参数含义的介绍。含义的介绍。距离大于距离大于rc的粒子的的粒子的相互作用力可以忽略相互作用力可以忽略2021/6/16124CREALX(768),VH(768),F(768)CCREALFY(256),FZ(256),Y(256),Z(256)CREALH,HSQ,HSQ2,TREFREALKX,KY,KZCINTEGERCLOCK,TIMEMXCEQUIVALENCE(FY,F(257),(FZ,F(513),(Y,X(257),(Z,X(513)CC该部分主要为对变量进行该部分主要为对变量进行数据类型定义数据类型定义2021/6/16125C-CCDEFINITIONOFTHESIMULATIONPARAMETERSCC-CNPART=256;DEN=0.83134SIDE=6.75284TREF=0.722RCOFF=2.5H=0.064IREP=50ISTOP=500TIMEMX=3000CISEED=4771该部分主要为对变量进行该部分主要为对变量进行赋值赋值2021/6/16126CC-CCCSETTHEOTHERPARAMETERSCCC-CWRITE(*,*)MOLECULARDYNAMICSSIMULATIONPROGRAMWRITE(*,*)WRITE(*,*)NUMBEROFPARTICLEIS,NPARTWRITE(*,*)SIDELENGTHOFTHEBOXIS,SIDEWRITE(*,*)CUTOFFIS,RCOFFWRITE(*,*)REDUCEDTEMPERATUREIS,TREFWRITE(*,*)BASICTIMESTEPIS,HWRITE(*,*)CNPART粒子数，粒子数，SIDE-盒子边长，盒子边长，RCOFF-势能作用截止距离，势能作用截止距离，TREF-约化温度，约化温度，H-基本时间步基本时间步2021/6/16127OPEN(1,FILE=A.TXT,STATUS=UNKNOWN)A=SIDE/4.0SIDEH=SIDE*0.5HSQ=H*HHSQ2=HSQ*0.5NPARTM=NPART-1RCOFFS=RCOFF*RCOFFTSCALE=16.0/(1.0*NPART-1.0)VAVER=1.13*SQRT(TREF/24.0)CN3=3*NPARTIOF1=NPARTIOF2=2*NPARTCCCALLRANSET(ISEED)C为系综达到预定平均温度为系综达到预定平均温度要进行标度并监视平均速度要进行标度并监视平均速度TSCALE-标度因子的一部分标度因子的一部分2021/6/16128CC-CCTHISPARTOFTHEPROGRAMPREPARESTHEINITIALCONFIGURATIONCC-C说明程序这部分是进行说明程序这部分是进行256个粒子初始位置的设定个粒子初始位置的设定2021/6/16129IJK=0DO10LG=0,1DO10I=0,3DO10J=0,3DO10K=0,3IJK=IJK+1X(IJK)=I*A+LG*A*0.5Y(IJK)=J*A+LG*A*0.5Z(IJK)=K*A10CONTINUEDO15LG=1,2DO15I=0,3DO15J=0,3DO15K=0,3IJK=IJK+1X(IJK)=I*A+(2-LG)*A*0.5Y(IJK)=J*A+(LG-1)*A*0.5Z(IJK)=K*A+A*0.515CONTINUECSETUPFCCLATTICEFORTHEATOMSINSIDETHEBOXC=C256个粒子初始位置的设定个粒子初始位置的设定是分布在面心立方上，是分布在面心立方上，A=SIDE/4.02021/6/16130CCASSIGNVELOCITIESDISTRIBUTEDNORMALLYC=CCALLMXWELL(VH,N3,H,TREF)给给256个粒子设定初始的速度，个粒子设定初始的速度，速度分布应符合麦克斯韦分布，速度分布应符合麦克斯韦分布，这个设定初始速度的子程序为这个设定初始速度的子程序为MXWELL(VH,N3,H,TREF)。DO50I=1,N3F(I)=0.050CONTINUEN3=3*NPART2021/6/16131C-CCSTARTOFTHEACTUALMOLECULARDYNAMICSPROGRAM.CPRENTICCTHEEQUATIONSOFMOTIONAREINTEGRATEDUSINGTHESUMMEDFORMC(G.DAHLQUISTANDA.BJOERK,NUMERICALMETHODS,EHALLC(1974).EVERYIREPTHSTEPTHEVELOCITIESARERESCALEDSOASCTOGIVETHESPECIFIEDTEMPERATURETREFCCVERSION1.0AUGUST1985CDIETERW.HEERMANNCC-C程序说明，表示以下将程序说明，表示以下将开始一个真正的分子开始一个真正的分子动力学程序动力学程序并告诉程序版本并告诉程序版本2021/6/16132DO200CLOCK=1,TIMEMXCC=ADVANCEPOSITIONSONEBASICTIMESTEP=CDO210I=1,N3X(I)=X(I)+VH(I)+F(I)210CONTINUEX(I)VH(I)F(I)2021/6/16133CC=APPLYPERIODICBOUNDARYCONDITIONS=CDO215I=1,N3IF(X(I).LT.0)X(I)=X(I)+SIDEIF(X(I).GT.SIDE)X(I)=X(I)-SIDE215CONTINUEC处理周期性边界条件处理周期性边界条件如果有一个粒子穿过基本元胞的如果有一个粒子穿过基本元胞的一个表面，那么这个粒子就穿过一个表面，那么这个粒子就穿过对面的墙重新进入元胞，速度不变对面的墙重新进入元胞，速度不变2021/6/16134C=COMPUTETHEPARTIALVELOCITIES=CDO220I=1,N3VH(I)=VH(I)+F(I)220CONTINUE计算该积分计算该积分步的速度步的速度C=CTHISPARTCOMPUTESTHEFORCEONTHEPARTICLESC=CCCVIR=0.0EPOT=0.0DO225I=1,N3F(I)=0.0225CONTINUE计算作用在粒子上的力计算作用在粒子上的力EPOT势能势能初始作用在粒子上的力为初始作用在粒子上的力为02021/6/16135CDO270I=1,NPARTX1=X(I)Y1=Y(I)Z1=Z(I)DO270J=I+1,NPARTXX=X1-X(J)YY=Y1-Y(J)ZZ=Z1-Z(J)CIF(XX.LT.-SIDEH)XX=XX+SIDEIF(XX.GT.SIDEH)XX=XX-SIDECIF(YY.LT.-SIDEH)YY=YY+SIDEIF(YY.GT.SIDEH)YY=YY-SIDECIF(ZZ.LT.-SIDEH)ZZ=ZZ+SIDEIF(ZZ.GT.SIDEH)ZZ=ZZ-SIDECRD=XX*XX+YY*YY+ZZ*ZZIF(RD.GT.RCOFFS)GOTO270C计算粒子之间的距离，计算粒子之间的距离，并判断其是否超过了并判断其是否超过了位势截断（止）距离位势截断（止）距离2021/6/16136EPOT=EPOT+RD*(-6.0)-RD*(-3.0)R148=RD*(-7.0)-0.5*RD*(-4.0)VIR=VIR-RD*R148KX=XX*R148F(I)=F(I)+KXF(J)=F(J)-KXKY=YY*R148FY(I)=FY(I)+KYFY(J)=FY(J)-KYKZ=ZZ*R148FZ(I)=FZ(I)+KZFZ(J)=FZ(J)-KZ270CONTINUE2021/6/16137DO275I=1,N3F(I)=F(I)*HSQ2275CONTINUECC=CCENDOFTHEFORCECALCULATIONCC力乘上力乘上（1/2）h2C=COMPUTETHEVELOCITES=DO300I=1,N3VH(I)=VH(I)+F(I)300CONTINUECC=COMPUTETHEKINETICENERGY=CEKIN=0.0DO305I=1,N3EKIN=EKIN+VH(I)*VH(I)305CONTINUEEKIN=EKIN/HSQHSQ=h22021/6/16138C=COMPUTEAVERAGEVELOCITY=CVEL=0.0COUNT=0.0DO306I=1,NPARTVX=VH(I)*VH(I)VY=VH(I+IOF1)*VH(I+IOF1)VZ=VH(I+IOF2)*VH(I+IOF2)SQ=SQRT(VX+VY+VZ)SQT=SQ/HIF(SQT.GT.VAVER)COUNT=COUNT+1VEL=VEL+SQ306CONTINUEVEL=VEL/HN3=3*NPARTIOF1=NPARTIOF2=2*NPART计算平均速度，并进行监视计算平均速度，并进行监视VAVER=1.13*SQRT(TREF/24.0)2021/6/16139IF(CLOCK.LT.ISTOP)THENC=NORMALIZETHEVELOCITIESTOOBTAINTHE=C=SPECIFIEDREFERENCETEMPERATURE=IF(CLOCK/IREP)*IREP.EQ.CLOCK)THENWRITE(*,*)VELOCITYADJUSTMENTTS=TSCALE*EKINWRITE(*,*)TEMPERATUREBEFORESCALINGIS,TSSC=TREF/TSSC=SQRT(SC)WRITE(*,*)SCALEFACTORIS,SCDO310I=1,N3VH(I)=VH(I)*SC310CONTINUEEKIN=TREF/TSCALEENDIFENDIFCTSCALE=16.0/(1.0*NPART-1.0)2021/6/16140C=CCOMPUTEVARIOUSQUANTITIESC=CIF(CLOCK/2)*2.EQ.CLOCK)THENEK=24.0*EKINEPOT=4.0*EPOTETOT=EK+EPOTTEMP=TSCALE*EKINPRES=DEN*16.0*(EKIN-VIR)/NPARTVEL=VEL/NPARTRP=(COUNT/256.0)*100.0CWRITE(*,6000)CLOCK,EK,EPOT,ETOT,TEMP,PRES,VEL,RPWRITE(1,6000)CLOCK,EK,EPOT,ETOTENDIFC200CONTINUE6000FORMAT(1X,I8,F15.6,F15.6,F15.6)CSTOPEND计算别的量，包括总动能、计算别的量，包括总动能、总势能、温度、压强等总势能、温度、压强等2021/6/16141CC=CCMXWELLRETURNSUPONCALLINGAMAXWELLDISTRIBUTEDDEVIATESCFORTHESPECIFIEDTEMPERATURETREF.ALLDEVIATESARESCALEDBYCAFACTORCCCALLINGPARAMETERSAREASFOLLOWSCCVHVECTOROFLENGTHNPART(MUSTBEAMULTIPLEOF2)CN3VECTORLENGTHCHSCALINGFACTORWITHWHICHALLELEMENTSOFVHAREMULTIPLIEDCTREFTEMPERATURECCVERSION1.0ASOFAUGUST1985CDIETERW.HEERMANNCC=C2021/6/16142SUBROUTINEMXWELL(VH,N3,H,TREF)REALVH(N3)CREALU1,U2,V1,V2,S,RISEED=-1CNPART=N3/3IOF1=NPARTIOF2=2*NPARTTSCALE=16.0/(1.0*NPART-1.0)CDO10I=1,N3,21U1=RNDF(ISEED)U2=RNDF(ISEED)CV1=2.0*U1-1.0V2=2.0*U2-1.0S=V1*V1+V2*V2CIF(S.GE.1.0)GOTO1R=-2.0*LOG(S)/SVH(I)=V1*SQRT(R)VH(I+1)=V2*SQRT(R)10CONTINUE产生正态随机数的一个方法是极坐标产生正态随机数的一个方法是极坐标方法。优点是同时产生两个独立的正方法。优点是同时产生两个独立的正态分布随机变量，而实际上并不增加态分布随机变量，而实际上并不增加计算时间。步骤如下：计算时间。步骤如下：（1）产生两个在）产生两个在0，1区间上均匀区间上均匀分布的独立随机变量分布的独立随机变量U1，U2；（2）令）令V1=2U1-1V2=2U2-1；（3）计算）计算S=V12+V22；（4）若）若S 1，则回到第，则回到第1步；步；（5）否则）否则R=-2lnS/SX1=V1*SQRT(-2lnS)/SX2=V2*SQRT(-2lnS)/S随机数产生随机数产生子程序子程序2021/6/16143EKIN=0.0SP=0.0DO20I=1,NPARTSP=SP+VH(I)20CONTINUESP=SP/NPARTDO21I=1,NPARTVH(I)=VH(I)-SPEKIN=EKIN+VH(I)*VH(I)21CONTINUEWRITE(*,*)TOTALLINEARMOMENTUMINXDIRECTIONIS,SP正态随机分布向麦克斯正态随机分布向麦克斯韦分布转变（韦分布转变（x方向方向）0vf(v)麦克斯韦分布麦克斯韦分布x(x)正态分布正态分布2021/6/16144SP=0.0DO22I=IOF1+1,IOF2SP=SP+VH(I)22CONTINUESP=SP/NPARTDO23I=IOF1+1,IOF2VH(I)=VH(I)-SPEKIN=EKIN+VH(I)*VH(I)23CONTINUEWRITE(*,*)TOTALLINEARMOMENTUMINYDIRECTIONIS,SP正态随机分布向麦克正态随机分布向麦克斯韦分布转变（斯韦分布转变（y 方方向向）如果进行以下变换：如果进行以下变换：2021/6/16145SP=0.0DO24I=IOF2+1,N3SP=SP+VH(I)24CONTINUESP=SP/NPARTDO25I=IOF2+1,N3VH(I)=VH(I)-SPEKIN=EKIN+VH(I)*VH(I)25CONTINUEWRITE(*,*)TOTALLINEARMOMENTUMINZDIRECTIONIS,SPWRITE(*,*)正态随机分布向麦克斯正态随机分布向麦克斯韦分布转变（韦分布转变（z方向方向）则可以发现，只要将正态则可以发现，只要将正态分布变量分布变量x（y、z）减去减去x（y、z）的数学期望的数学期望 就得就得到了麦克斯分布，而数学到了麦克斯分布，而数学期望就是均值。期望就是均值。2021/6/16146WRITE(*,*)VELOCITYADJUSTMENTTS=TSCALE*EKINWRITE(*,*)TEMPERATUREBEFORESCALINGIS,TSSC=TREF/TSSC=SQRT(SC)WRITE(*,*)SCALEFACTORIS,SCSC=SC*HDO30I=1,N3VH(I)=VH(I)*SC30CONTINUEEND初始的速度也需要初始的速度也需要标度，并确定标度标度，并确定标度因子使其等于总能因子使其等于总能量，但这并不改变量，但这并不改变分布规律分布规律2021/6/16147FUNCTIONRNDF(IDUM)DIMENSIONR(97)CPARAMETER(M1=259200,IA1=7141,IC1=54773,RM1=1./M1)CPARAMETER(M2=134456,IA2=8121,IC2=28411,RM2=1./M2)PARAMETER(M1=259200,IA1=7141,IC1=54773,RM1=3.85802E-6)PARAMETER(M2=134456,IA2=8121,IC2=28411,RM2=7.43738E-6)PARAMETER(M3=243000,IA3=4561,IC3=51349)DATAIFF/0/IF(IDUM.LT.0.OR.IFF.EQ.0)THENIFF=1IX1=MOD(IC1-IDUM,M1)IX1=MOD(IA1*IX1+IC1,M1)IX2=MOD(IX1,M2)IX1=MOD(IA1*IX1+IC1,M1)IX3=MOD(IX1,M3)DO11J=1,97IX1=MOD(IA1*IX1+IC1,M1)IX2=MOD(IA2*IX2+IC2,M2)R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM111CONTINUEIDUM=1ENDIF此处为随机数产生程序此处为随机数产生程序（函数）（函数）2021/6/16148IX1=MOD(IA1*IX1+IC1,M1)IX2=MOD(IA2*IX2+IC2,M2)IX3=MOD(IA3*IX3+IC3,M3)J=1+(97*IX3)/M3IF(J.GT.97.OR.J.LT.1)PAUSERNDF=R(J)R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1RETURNEND此处为随机数产生程序此处为随机数产生程序（函数）（函数）2021/6/161492021/6/16150 结束语结束语若有不当之处，请指正，谢谢！若有不当之处，请指正，谢谢！