Chapter1, Example1:半经典量子化条件求分子振动问题-一根柳条-ChinaUnix博客

想得美

首页　| 　博文目录　| 　关于我

一根柳条

博客访问： 40901
博文数量： 10
博客积分： 0
博客等级：民兵
技术积分： 95
用户组：普通用户
注册时间： 2015-10-18 15:21

文章分类

全部博文（10）

Gnuplot（0）
Computational&nb（6）
Fortran（4）
未分配的博文（0）

文章存档

2015年（10）

我的朋友

windhawk

相关博文

Chapter1, Example1:半经典量子化条件求分子振动问题

分类： C/C++

2015-10-27 00:16:30

　　这是书中给出的第一个长案例，在Appendix B Programs for the Examples中有源代码。使用了很多公用的程序包，这给我们分析带来不便。毕竟，我们不会再去使用他的那些程序包了。因此，要想重复这项工作，得做好切割。
　　这个案例给我的启发相当多。最大的冲击，首先就是其标准化。程序中的注释相当多，相当翔实。这么大的一个程序，其主程序却没几行。分成很多的子程序和函数，每个子程序和函数都实现一部分相对独立的功能。对于我这种菜鸟而言，看到这么漂亮的程序，真是一种震撼！
　　物理问题复述：
　　势函数：Lennard-Jones势函数
　　 This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

　　其最小值为

,此时势阱深度为V₀。
　　在量子力学中，该问题自然是用一维Schrodinger方程给出：
　　 This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

，m为约化质量。
　　在这里使用量子化规则求解。本工作的重点之一便是求出各个能级En的数值。
　　其中，能量为：
　　 This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

，
　　而且能量满足： This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

　　分子的来回振荡，实现着动能和势能的相互转化。分子间距在r_in和r_out间振荡，当处在r_in时，分子间距压缩到最小，动能为0；当处在r_out时，分子间距拉伸到最大，动能为0。当处在 This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

时，分子动能最大。
　　求解可知： This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

，注意：动量是坐标r的函数！
　　令
　　 This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

，则根据量子化条件有
　　 This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

　　根据无参化处理
　　

，

均为无单位的参量。对于氢分子，我们取 This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

＝21.7。
　　求解过程应该是先求拐点值。拐点值是由能量En值给出的。当然，为了满足En小于零（束缚态），则要求拐点x_in>1成立。其次，能量值En则是由量子化条件给出的。也就是先给出一个尝试解En，根据这个尝试解会求出相应的拐点，然后便可以求出积分值，进而与量子化条件进行比较，由此修改尝试解En。如此，周而复始。
　　

点击(此处)折叠或打开

!ccccccccccccccccccccccccccccccccccccccccccccccccc!
module init
!初始化参数，并给出欢迎界面
!ccccccccccccccccccccccccccccccccccccccccccccccccc!
! Global variables
　implicit none
public ! 凡是没特别声明的参数都公开
　real,parameter :: potmin=2.0**(1.0/6.0) ! 阱底处的x值
　real,parameter :: gammavalue=21.7
　real,parameter :: etol=0.0005 ! 寻找能量时的最小误差
　real,parameter :: xtol=0.0005 ! 寻找拐点时的最小误差
　integer,parameter :: npts=100 ! 积分过程分段数目
　integer :: nlevel ! 束缚态总数，该数值由gamma值决定。
　integer,parameter :: Maxlevel=100 !限定束缚态最大数目
　real,parameter :: pi=4.0*atan(1.0)
　common nlevel
! Local variables
end module
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
!cccccccccccccccccccccccccccccccccccccccccccccccc!
program example1
! example1.f90
! main function
!cccccccccccccccccccccccccccccccccccccccccccccccc!
　use init
　call archon
stop
end program example1
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
!ccccccccccccccccccccccccccccccccccccccccccccccccc!
subroutine archon
!此处调用search子程序，得到第n个束缚态
!而束缚态的数目是在这个子程序中确定的，
!search子程序只是按照输入参量（从本子程序中输入的）计算相应的束缚态能量En
!
!ccccccccccccccccccccccccccccccccccccccccccccccccc!
! Global variables
　use init
! Input/Output variables
! Local variables
　real :: Eo
　real :: s
　real :: E1
　real :: x1,x2
　real :: F1
　integer :: Ilevel
　real :: energy(0:Maxlevel)=0.0
　real :: xin(0:Maxlevel)=0.0
　real :: xout(0:Maxlevel)=0.0
!ccccccccccccccccccccccccccccccccccccccccccccccccc!
! calculate total number of levels
　E0=-etol
　call actions(E0,x1,x2,s)
　nlevel=Int(s/pi-0.5)+1
　if((nlevel-1)>Maxlevel)then
　　write(*,*)"gamma value is to large!"
　　return
　end if
!
　E1=-1.0
　F1=-pi/2.0
　do Ilevel=0,nlevel-1,1
　　call search(Ilevel,E1,F1,x1,x2)
　　energy(Ilevel)=E1
　　xin(Ilevel)=x1
　　xout(Ilevel)=x2
　　F1=F1-pi
　end do
　write(*,*)"the number of bound states:",nlevel
　do Ilevel=0,nlevel-1,1
　　write(*,*)"energy(",Ilevel,"):",energy(Ilevel)
　end do
return
end subroutine
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
!ccccccccccccccccccccccccccccccccccccccccccccccccc!
subroutine search(N,E1,F1,x1,x2)
!寻找第n个束缚态
!本子程序还要调用actions子程序，本程序给出一个尝试能量E，然后actions便计算出s和x1,x2
!本程序再根据actions返回的计算值s,x1,x2，判断该尝试能量是否合理。
!如果不合理，再尝试新的能量。如此，周而复始。
!ccccccccccccccccccccccccccccccccccccccccccccccccc!
! Global variables
　use init
! Input/Output variables
　integer :: N ! Current level (Input)
　real :: E1,E2 ! trial energies (I/O)
　real :: F1,F2 ! f=action/2-pi*(n+1/2)(I/O)
　real :: s ! action (Output),用于存储actions计算出的s数值
　real :: x1,x2 ! turning points (Output)
! Local variables
　real :: de
!cccccccccccccccccccccccccccccccccccccccccccccccccc!
　E2=E1+abs(E1)/4.0 !使用割线法求根，此处便是寻找的下一个点，二者连接，求第三点
　de=2.0*etol
　do while(abs(de)>etol)
　　call actions(E2,x1,x2,s)
　　F2=s-(N+0.5)*pi
　　if(F1/=F2)then
　　　de=-F2*(E2-E1)/(F2-F1)
　　else
　　　de=0.0
　　end if
　　E1=E2
　　F1=F2
　　E2=E1+de
　　if(E2>0)E2=-etol
　end do
return
end subroutine
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
!ccccccccccccccccccccccccccccccccccccccccccccccccc!
subroutine actions(E,x1,x2,s)
!该子程序是求拐点的，原理便是pot(x)-E=0.
!根据输入的能量E值来寻根，确定拐点。
! 此外，该程序还要求出作用量s
!cccccccccccccccccccccccccccccccccccccccccccccccccc!
! Global variables:变量的分类这一点做的很漂亮，注释的很翔实。
　use init
! Input/Output variables:
　real :: E ! energy (Input)
　real :: s !action (Output)
　real :: x1,x2 ! turning points (Output) 漂亮！
! Local variables:
　real :: dx ! increment in turning search
　real :: H ! quadrature step size, H=(x2-x1)/npts
　real :: x,total
　integer :: Ix
　real,external :: pot
!cccccccccccccccccccccccccccccccccccccccccccccccccc!
　x1=potmin 　!从最小值点往左找拐点
　dx=0.1
　do while(dx>xtol)
　　x1=x1-dx
　　if(pot(x1)>E)then
　　　x1=x1+dx
　　　dx=dx/2.0
　　end if
　end do
!
　x2=potmin
　dx=0.1
　do while(dx>xtol)
　　x2=x2+dx
　　if(pot(x2)>E)then
　　　x2=x2-dx
　　　dx=dx/2.0
　　end if
　end do
!
! Simpson's rule from x1+H to x2-H
　H=(x2-x1)/real(npts)
　total=0.0
　x=x1+H
　do Ix=1,npts/2,1
　　total=total+H/3.0*(sqrt(E-pot(x-H))+4.0*sqrt(E-pot(x))+sqrt(E-pot(x+H)))
　　x=x+2*H
　end do
!
　s=total*gammavalue
return
end
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
!cccccccccccccccccccccccccccccccccccccccccccccccccc!
real function pot(x)
! the Lennard-Jones potential
!cccccccccccccccccccccccccccccccccccccccccccccccccc!
! Passed variables
　real :: x
　pot=4*(x**(-12)-x**(-6))
return
end
!ccccccccccccccccccccccccccccccccccccccccccccccccccc!

　　第一次捣鼓这么复杂的代码，竟然就运行通了。

运行结果为：
　　the number of bound states:           5
　　energy( 0 ): 　-0.772370696
　　energy( 1 ):　 -0.422584146
　　 energy( 2 ):　 -0.195745081
　　energy( 3 ):　 -6.75332770E-02
　　energy( 4 ): 　-1.26941651E-02
　　为了多输出一些结果，略微修改一下。输出信息是：
　　阱口处：x1:   1.00058699     x2:   4.47167873
　　总的束缚态数目：the number of bound states: 5
　　energy( 0 ): -0.772370696
　　xin( 0 ):   1.05214953
　　xout( 0 ):   1.24980581
＝＝＝＝＝＝＝＝＝＝＝＝＝＝
　　energy(1 ): -0.422584146
　　xin( 1 ):   1.02168071
　　 xout( 1 ):   1.42324340
＝＝＝＝＝＝＝＝＝＝＝＝＝＝
　　 energy( 2 ): -0.195745081
　　xin( 2 ):   1.00918067
　　xout( 2 ):   1.63886845
＝＝＝＝＝＝＝＝＝＝＝＝＝＝
　　energy( 3 ): -6.75332770E-02
　　xin(3 ):   1.00293064
　　 xout( 3 ):   1.96855593
＝＝＝＝＝＝＝＝＝＝＝＝＝＝
　　energy( 4 ): -1.26941651E-02
　　xin( 4 ):   1.00058699
　　xout(4 ):   2.60683656
　　至少从内拐点和外拐点的信息上看还是合理的，能量越高，内拐点位置越接近１.0，但永远不会达到１.0。外拐点，则随着能级的升高而变大。　
　　此外，书中还提到了Morse势函数。这种势函数是分子力的一种。
　　 This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

　　当要求在位置零点r=0时，势函数为排斥势，即势函数大于零。则需要满足关系式： This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

　　在对函数进行无参化处理之后， This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

。而且容易看出，该势函数形状与Lennard-Jones势函数差不多，只是在位置零点r=0时为有限值（且大于零），而且极小值点处洼的更厉害！
　　其他参数为：
　　 This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

，　

　　当取

时，出现15个束缚态，然后调整 This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

的数值，来看得到的结果。第一组结果是 This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

得到的，第二组是

得到的。

只列出一组数据吧，别的意思差不多，粘贴起来挺麻烦。
the number of bound states:          15
energy(0): -0.937868953                   energy(0): -0.937868953
xin(0): 0.471125007                          xin( 0 ): 0.571124911
xout(0): 0.978937566                       xout( 0 ):   1.07893741
==============================================
energy( 1 ): -0.819242954
xin( 1 ): 0.339093745
xout( 1 ):   1.24612498
energy(2 ): -0.708644629
xin(2 ): 0.261750013
xout( 2 ):   1.46878135
energy(3 ): -0.606048405
xin(3 ): 0.206281260
xout(3 ):   1.68050015
energy(4 ): -0.511439860
xin(4 ): 0.163312510
xout(4 ):   1.89300025
energy(5 ): -0.424902737
xin(5 ): 0.128937513
xout(5 ):   2.11253142
energy(6 ): -0.346358627
xin(6 ): 0.100812525
xout(6 ):   2.34534383
energy(7 ): -0.275833398
xin(7 ):   7.73750171E-02
xout(7 ):   2.59612465
energy(8 ): -0.213712454
xin(8 ):   5.86250164E-02
xout(8 ):   2.87034345
energy(9 ): -0.158712506
xin(9 ):   4.22187671E-02
xout(9 ):   3.18440557
energy(10 ): -0.112374529
xin(10 ):   2.97187679E-02
xout(10 ):   3.54221773
energy(11 ): -7.39284754E-02
xin(11 ):   1.95625164E-02
xout(11 ):   3.97112346
energy(12 ): -4.35847901E-02
xin(12 ):   1.09687662E-02
xout(12 ):   4.50784159
energy(13 ): -2.13147216E-02
xin(13 ):   5.50001580E-03
xout(13 ):   5.22893524
energy(14 ): -6.77488185E-03
xin(14 ):   1.59376580E-03
xout(14 ):   6.37893438
　　从中可以看出一个很有意思的现象。改变 This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

值可以改变束缚态数目，但是改变 This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.

仅是把内拐点和外拐点平移了而已，对能量没有任何改变！这在分析之前是没有想到的。

阅读(3664) | 评论(0) | 转发(0) |

上一篇： Chapter 1 Basic Mathematical Operations（3）寻根

下一篇：关于Fortran中子程序与函数不同点的一些认识

给主人留下些什么吧！~~

感谢所有关心和支持过ChinaUnix的朋友们

16024965号-6