Chinaunix首页 | 论坛 | 博客
  • 博客访问: 39170
  • 博文数量: 10
  • 博客积分: 0
  • 博客等级: 民兵
  • 技术积分: 95
  • 用 户 组: 普通用户
  • 注册时间: 2015-10-18 15:21
文章分类

全部博文(10)

文章存档

2015年(10)

我的朋友

分类: C/C++

2015-10-27 00:16:30

  这是书中给出的第一个长案例,在Appendix B  Programs for the Examples中有源代码。使用了很多公用的程序包,这给我们分析带来不便。毕竟,我们不会再去使用他的那些程序包了。因此,要想重复这项工作,得做好切割。
  这个案例给我的启发相当多。最大的冲击,首先就是其标准化。程序中的注释相当多,相当翔实。这么大的一个程序,其主程序却没几行。分成很多的子程序和函数,每个子程序和函数都实现一部分相对独立的功能。对于我这种菜鸟而言,看到这么漂亮的程序,真是一种震撼!
  物理问题复述:
  势函数:Lennard-Jones势函数
  

  其最小值为,此时势阱深度为V0
  在量子力学中,该问题自然是用一维Schrodinger方程给出:
  ,m为约化质量。
  在这里使用量子化规则求解。本工作的重点之一便是求出各个能级En的数值
  其中,能量为:
  
  而且能量满足:
  分子的来回振荡,实现着动能和势能的相互转化。分子间距在rin和rout间振荡,当处在rin时,分子间距压缩到最小,动能为0;当处在rout时,分子间距拉伸到最大,动能为0。当处在时,分子动能最大。
  求解可知:,注意:动量是坐标r的函数!
  令
  ,则根据量子化条件有
  
  
  根据无参化处理
  
  均为无单位的参量。对于氢分子,我们取=21.7。
  求解过程应该是先求拐点值。拐点值是由能量En值给出的。当然,为了满足En小于零(束缚态),则要求拐点xin>1成立。其次,能量值En则是由量子化条件给出的。也就是先给出一个尝试解En,根据这个尝试解会求出相应的拐点,然后便可以求出积分值,进而与量子化条件进行比较,由此修改尝试解En。如此,周而复始。
  

点击(此处)折叠或打开

  1. !ccccccccccccccccccccccccccccccccccccccccccccccccc!
  2. module init
  3. !初始化参数,并给出欢迎界面
  4. !ccccccccccccccccccccccccccccccccccccccccccccccccc!
  5. ! Global variables
  6.  implicit none
  7. public ! 凡是没特别声明的参数都公开
  8.  real,parameter :: potmin=2.0**(1.0/6.0) ! 阱底处的x值
  9.  real,parameter :: gammavalue=21.7
  10.  real,parameter :: etol=0.0005 ! 寻找能量时的最小误差
  11.  real,parameter :: xtol=0.0005 ! 寻找拐点时的最小误差
  12.  integer,parameter :: npts=100 ! 积分过程分段数目
  13.  integer :: nlevel ! 束缚态总数,该数值由gamma值决定。
  14.  integer,parameter :: Maxlevel=100 !限定束缚态最大数目
  15.  real,parameter :: pi=4.0*atan(1.0)
  16.  common nlevel
  17. ! Local variables

  18. end module
  19. !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
  20. !cccccccccccccccccccccccccccccccccccccccccccccccc!
  21. program example1
  22. ! example1.f90
  23. ! main function
  24. !cccccccccccccccccccccccccccccccccccccccccccccccc!
  25.  use init
  26.  call archon

  27. stop

  28. end program example1
  29. !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
  30. !ccccccccccccccccccccccccccccccccccccccccccccccccc!
  31. subroutine archon
  32. !此处调用search子程序,得到第n个束缚态
  33. !而束缚态的数目是在这个子程序中确定的,
  34. !search子程序只是按照输入参量(从本子程序中输入的)计算相应的束缚态能量En
  35. !
  36. !ccccccccccccccccccccccccccccccccccccccccccccccccc!
  37. ! Global variables
  38.  use init
  39. ! Input/Output variables

  40. ! Local variables
  41.  real :: Eo
  42.  real :: s
  43.  real :: E1
  44.  real :: x1,x2
  45.  real :: F1
  46.  integer :: Ilevel
  47.  real :: energy(0:Maxlevel)=0.0
  48.  real :: xin(0:Maxlevel)=0.0
  49.  real :: xout(0:Maxlevel)=0.0
  50. !ccccccccccccccccccccccccccccccccccccccccccccccccc!
  51. ! calculate total number of levels
  52.  E0=-etol
  53.  call actions(E0,x1,x2,s)
  54.  nlevel=Int(s/pi-0.5)+1
  55.  if((nlevel-1)>Maxlevel)then
  56.   write(*,*)"gamma value is to large!"
  57.   return
  58.  end if
  59. !
  60.  E1=-1.0
  61.  F1=-pi/2.0

  62.  do Ilevel=0,nlevel-1,1
  63.   call search(Ilevel,E1,F1,x1,x2)
  64.   energy(Ilevel)=E1
  65.   xin(Ilevel)=x1
  66.   xout(Ilevel)=x2
  67.   F1=F1-pi
  68.  end do
  69.  write(*,*)"the number of bound states:",nlevel
  70.  do Ilevel=0,nlevel-1,1
  71.   write(*,*)"energy(",Ilevel,"):",energy(Ilevel)
  72.  end do
  73. return
  74. end subroutine
  75. !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
  76. !ccccccccccccccccccccccccccccccccccccccccccccccccc!
  77. subroutine search(N,E1,F1,x1,x2)
  78. !寻找第n个束缚态
  79. !本子程序还要调用actions子程序,本程序给出一个尝试能量E,然后actions便计算出s和x1,x2
  80. !本程序再根据actions返回的计算值s,x1,x2,判断该尝试能量是否合理。
  81. !如果不合理,再尝试新的能量。如此,周而复始。
  82. !ccccccccccccccccccccccccccccccccccccccccccccccccc!
  83. ! Global variables
  84.  use init
  85. ! Input/Output variables
  86.  integer :: N ! Current level (Input)
  87.  real :: E1,E2 ! trial energies (I/O)
  88.  real :: F1,F2 ! f=action/2-pi*(n+1/2)(I/O)
  89.  real :: s ! action (Output),用于存储actions计算出的s数值
  90.  real :: x1,x2 ! turning points (Output)
  91. ! Local variables
  92.  real :: de
  93. !cccccccccccccccccccccccccccccccccccccccccccccccccc!
  94.  E2=E1+abs(E1)/4.0 !使用割线法求根,此处便是寻找的下一个点,二者连接,求第三点
  95.  de=2.0*etol

  96.  do while(abs(de)>etol)
  97.   call actions(E2,x1,x2,s)
  98.   F2=s-(N+0.5)*pi
  99.   if(F1/=F2)then
  100.    de=-F2*(E2-E1)/(F2-F1)
  101.   else
  102.    de=0.0
  103.   end if

  104.   E1=E2
  105.   F1=F2
  106.   E2=E1+de
  107.   if(E2>0)E2=-etol

  108.  end do

  109. return
  110. end subroutine
  111. !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!

  112. !ccccccccccccccccccccccccccccccccccccccccccccccccc!
  113. subroutine actions(E,x1,x2,s)
  114. !该子程序是求拐点的,原理便是pot(x)-E=0.
  115. !根据输入的能量E值来寻根,确定拐点。
  116. ! 此外,该程序还要求出作用量s
  117. !cccccccccccccccccccccccccccccccccccccccccccccccccc!
  118. ! Global variables:变量的分类这一点做的很漂亮,注释的很翔实。
  119.  use init
  120. ! Input/Output variables:
  121.  real :: E ! energy (Input)
  122.  real :: s !action (Output)
  123.  real :: x1,x2 ! turning points (Output) 漂亮!
  124. ! Local variables:
  125.  real :: dx ! increment in turning search
  126.  real :: H ! quadrature step size, H=(x2-x1)/npts
  127.  real :: x,total
  128.  integer :: Ix
  129.  real,external :: pot
  130. !cccccccccccccccccccccccccccccccccccccccccccccccccc!
  131.  x1=potmin  !从最小值点往左找拐点
  132.  dx=0.1

  133.  do while(dx>xtol)
  134.   x1=x1-dx
  135.   if(pot(x1)>E)then
  136.    x1=x1+dx
  137.    dx=dx/2.0
  138.   end if
  139.  end do
  140. !
  141.  x2=potmin
  142.  dx=0.1

  143.  do while(dx>xtol)
  144.   x2=x2+dx
  145.   if(pot(x2)>E)then
  146.    x2=x2-dx
  147.    dx=dx/2.0
  148.   end if
  149.  end do
  150. !
  151. ! Simpson's rule from x1+H to x2-H
  152.  H=(x2-x1)/real(npts)
  153.  total=0.0
  154.  x=x1+H
  155.  do Ix=1,npts/2,1
  156.   total=total+H/3.0*(sqrt(E-pot(x-H))+4.0*sqrt(E-pot(x))+sqrt(E-pot(x+H)))
  157.   x=x+2*H
  158.  end do
  159. !
  160.  s=total*gammavalue
  161. return

  162. end
  163. !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
  164. !cccccccccccccccccccccccccccccccccccccccccccccccccc!
  165. real function pot(x)
  166. ! the Lennard-Jones potential
  167. !cccccccccccccccccccccccccccccccccccccccccccccccccc!
  168. ! Passed variables
  169.  real :: x

  170.  pot=4*(x**(-12)-x**(-6))

  171. return

  172. end
  173. !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  
  至少从内拐点和外拐点的信息上看还是合理的,能量越高,内拐点位置越接近1.0,但永远不会达到1.0。外拐点,则随着能级的升高而变大。  
  此外,书中还提到了Morse势函数。这种势函数是分子力的一种。
  
  当要求在位置零点r=0时,势函数为排斥势,即势函数大于零。则需要满足关系式:
  在对函数进行无参化处理之后,。而且容易看出,该势函数形状与Lennard-Jones势函数差不多,只是在位置零点r=0时为有限值(且大于零),而且极小值点处洼的更厉害!
  其他参数为:
  , , 
  当取时,出现15个束缚态,然后调整的数值,来看得到的结果。第一组结果是得到的,第二组是得到的。只列出一组数据吧,别的意思差不多,粘贴起来挺麻烦。
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    
  从中可以看出一个很有意思的现象。改变值可以改变束缚态数目,但是改变仅是把内拐点和外拐点平移了而已,对能量没有任何改变!这在分析之前是没有想到的。
阅读(3457) | 评论(0) | 转发(0) |
给主人留下些什么吧!~~