现代分子理论与计算化学导论作业

《现代分子理论与计算化学导论》——课程大作业班级：xxxxxxx 姓名：小签牛学号：xxxxxxxxxx题目：在T*=1.5条件下，分别用分子模拟方法和微扰理论方法计算ρ*=0.02和0.85的体系的压力，并比较两种方法计算的结果。

Ⅰ.当T*=1.5、ρ*=0.02时的情况①由Monte Carlo模拟获得体系的内能、径向分布函数和压力，流体参数及模拟条件见contrifile文件；此时的contrifile文件为：---------------ENTER THE FOLLOWING IN LENNARD-JONES UNITS--------------------0.02 # Enter The Density1.5 # Enter The Temperature8.0 # Enter The Potential Cutoff Distance108 # Enter The Intial Molecular Number---------------ENTER THE SIMULATION STEP CONTROLLING PARAMETES--------------- 200000 # Enter Number Of Cycles400 # Enter Number Of Steps Between Output Lines400 # Enter Number Of Steps Between Data Saves400 # Enter Interval For Update Of Max. Displ..False. # Whether Read config. From Old Simulation Runconfig.dat # Enter The Configuration File Name---------------ENTER THE RADIAL DISTRIBUTION FUNCTION PARAMETES-------------- .True. # Whether Calculate The Radial Distribution Function0.01 # Enter The Radial Distribution Distance100000 # Enter Number Of Cycles Of Start Calculating The Radial Distributiongr0.02.dat # Enter The Radial Distribution File Name（运行程序见附件1）所得“result.dat”文件中的结果为：A VERAGES<V/N> = -0.149649<P> = 0.028542FLUCTUATIONSFLUCTUATION IN <V/N> = 0.029433 FLUCTUATION IN <P> = 0.002152END OF SIMULATIONMC 模拟径向分布函数图：利用MC 程序运行产生的数据文件“gr0.02.dat ”，使用origin8绘图软件作图：图1 MC 模拟法所得径向分布函数② 微扰理论法径向分布函数图：利用Lennard_Jones 程序运行产生的数据文件“gr0.02.dat ”，使用origin8绘图软件作图：123450.00.51.01.52.0g (r )rT *=1.5 ρ*=0.02图2 微扰理论法所得径向分布函数由径向分布函数计算体系的Helmholtz 自由能：根据*6ρπη=计算出η，进而由)1(3-422ηηη-=hs A 可计算出hs A 。

*2*)()(2T dr r r u r g A LJ LJ ⎪⎭⎫ ⎝⎛=⎰∞σπρ可转换为微分迭代求和获得。

最后由AA A LJhsex+=即可获得Helmholtz 自由能由Helmholtz 自由能计算体系的压缩因子，进而得到压力，所涉及公式如下：1*ex *+∂∂=ρρA Z ；T Z P ***=ρ（运行程序见附件2）所得 “result.dat”文件中的结果为：Helmholtz Energy <A> : -0.024717 Compressibility Factor <Z> : 0.975845 Pressure <P> : 0.029275123456780.00.51.01.52.0g (r )rT *=1.5 ρ*=0.02③ 结果分析对于r 比较小的情况，g(r)主要表征的是原子的堆积状况及各个键之间的距离。

对于长程的性质，由于对于给定的距离找到原子的几率基本上相同，所以g(r)随着r 的增大而变得平缓，最后趋向于恒定值。

通常g(r) 归一化的条件为r 趋向于无穷大 时，g(r)趋向于1。

分析图1和图2，可以看出：由Monte Carlo 模拟与用微扰理论法所得的径向分布函数走势基本一致。

而且随着r 的增大，g(r)趋向于恒定值1.通过MC 模拟的方法和微扰理论的方法都计算出了体系的压力，分别为0.028542和0.029275，通过MC 模拟法所得压力与微扰理论法所得压力的相对误差为2.50%，可见两种方法的偏差比较小。

从而，在同等温度下，当密度相对较小时，两种方法所得结果接近。

偏差分析： 用MC 模拟时，根据所设定的温度、密度和每个分子在的不同位置时的相互作用力计算出多个压力P ，当计算完成后输出的P 是一个加权平均值。

微扰理论方法是通过简单问题的精确解来求得复杂问题的近似解。

对于Helmholtz 状态函数的计算，P *的计算涉及到*2*/)()(2T dr r r u r g A LJ LJ ⎪⎭⎫ ⎝⎛=⎰∞σπρ、L J hs ex A A A +=、1**+∂∂=ρρex A Z 和***T Z P ρ=, A LJ 计算过程中采用微分迭代，dr 选取length/n=15/1024=0.015，以及在求解*ρ∂∂exA 也是用的(A2-A1) /(rou2-rou1)即(A2-A1)/（2step ）近似求解。

对于本身有误差的引用很有可能进一步的扩大了误差，这在一定程度上影响了积分的准确性；诸多方面造成了偏差的存在。

Ⅱ.当T*=1.5、ρ*=0.85时的情况① 由Monte Carlo 模拟获得体系的内能、径向分布函数和压力，流体参数及模拟条件见contrifile 文件；此时的contrifile 文件为：---------------ENTER THE FOLLOWING IN LENNARD-JONES UNITS--------------------0.85 # Enter The Density 1.5# Enter The Temperature8.0 # Enter The Potential Cutoff Distance108 # Enter The Intial Molecular Number---------------ENTER THE SIMULATION STEP CONTROLLING PARAMETES--------------- 200000 # Enter Number Of Cycles400 # Enter Number Of Steps Between Output Lines400 # Enter Number Of Steps Between Data Saves400 # Enter Interval For Update Of Max. Displ..False. # Whether Read config. From Old Simulation Runconfig.dat # Enter The Configuration File Name---------------ENTER THE RADIAL DISTRIBUTION FUNCTION PARAMETES-------------- .True. # Whether Calculate The Radial Distribution Function0.01 # Enter The Radial Distribution Distance100000 # Enter Number Of Cycles Of Start Calculating The Radial Distribution gr0.85.dat # Enter The Radial Distribution File Name（运行程序见附件1）所得“result.dat”文件中的结果为：A VERAGES<V/N> = -5.096817<P> = 4.958224FLUCTUATIONSFLUCTUATION IN <V/N> = 0.139289FLUCTUATION IN <P> = 0.660584END OF SIMULATIONMC模拟径向分布函数图：利用MC程序运行产生的数据文件“gr0.85.dat”，使用origin8绘图软件作图：2.5T*=1.5 ρ*=0.852.01.5)r(g1.00.50.00.5 1.0 1.5 2.0 2.5r图3 MC模拟法所得径向分布函数②微扰理论法径向分布函数图：利用Lennard_Jones程序运行产生的数据文件“gr0.85.dat”，使用origin8绘图软件作图：2.5T*=1.5 ρ*=0.852.01.5)r(g1.00.50.001234567r图4 微扰理论法所得径向分布函数由径向分布函数计算体系的Helmholtz 自由能；根据*6ρπη=计算出η，进而由)1(3-422ηηη-=hs A 可计算出hs A 。

*2*)()(2T dr r r u r g A LJ LJ ⎪⎭⎫ ⎝⎛=⎰∞σπρ可转换为迭代求和获得。

最后由AA A LJhsex +=即可获得Helmholtz 自由能。

由Helmholtz 自由能计算体系的压缩因子，进而得到压力，所涉及的公式如下：1*ex *+∂∂=ρρA Z ；T Z P ***=ρ（运行程序见附件3）所得 “result.dat”文件中的结果为：Helmholtz Energy <A> : 1.655380 Compressibility Factor <Z> : 6.903019 Pressure <P> : 8.801349③ 结果分析对于r 比较小的情况，g(r)主要表征的是原子的堆积状况及各个键之间的距离。