加载中…
个人资料
微博
  • 博客等级：
  • 博客积分：
  • 博客访问：
  • 关注人气：
  • 获赠金笔：0支
  • 赠出金笔：0支
  • 荣誉徽章：
正文 字体大小：大 中 小

高斯展宽，洛伦慈展宽，Voigt展宽的fortran code

(2015-04-13 07:05:44)
标签：

杂谈

这里谈的展宽，是对量子化学正则分析计算出来的分子的线性光谱（红外，拉曼等）展开成带状峰。这是一种经验处理，只是为了视觉上更接近于实验得到的光谱。既然实验得到的是带状光谱，那为什么量子化学正则分析计算得到是线性光谱呢？原因很简单，量子化学计算的气相分子在0K下的振动光谱，是在势能面的最低点根据二次简谐近似得到的振动频率。这个频率-吸收强度数据对，体现到谱图上就是一根线。实际的光谱总是在一定的温度下的凝聚相（固体或液体）下的振动光谱，它是极大量分子（～Na, 10^23）的平均结果，每个分子都具有大小不同的动能，处于不断运动、碰撞中，分布在各种激发态能级中。这些都会或多或少地改变其指定模式的振动频率，使得它们或大或小，按照一定的概率分布在特征频率的两侧，形成峰带。展宽是多种效应的总和效果。


按照量子电动力学，这些展宽效应分为两类：非均展宽，和均匀展宽。这个均和非均，并非是指展宽的峰形是否对称，而是指效应本身是否对所有的分子造成的后果是否一样。比如，杂质效应(dopant effect)，是指分子受到其中混杂的少量异种分子的作用后能级的变化。这种效应必然随着与杂质分子的距离而变化，因而属于非均展宽。然而，如果是分子在体系压力下造成的能级变化，那么对于所有的分子来说，振动频率受到影响的机会是均等的，因而属于均匀展宽。另外一个均匀展宽的例子是多普勒效应，即高速向你运动的物体，测得的频率会升高，而远离你的物体，频率会降低。非均展宽的概率分布符合gaussian分布，均匀展宽的概率分布则是Lorentzian型的。如果多种不同效应之间互相独立，则综合效果符合Voigt峰形。实际上各效应之间的耦合是不可避免的。但是Voigt展宽仍然是最符合实际的。


下面提供了3种展宽的fortran code。如果你要使用它，需要注意的是，它对输入数据的格式有特殊的要求，即频率-强度数据必须按照频率值的均匀间隔给出。比如，你高斯计算结果得到两组频率-强度值： 200， 12； 230，5；那么你应该这样准备输入文件：

...

190 0

200 12

210 0

220 0

230 5

240 0

...

如果你的频率值是 198.5呢？你必须减小频率间隔到你需要的精度，然后，原始频率值取整到最近的值。


为什么会这样呢？因为我写这个code并非为了处理量子化学计算结果，而是对分子动力学模拟得到的红外光谱结果进行平滑(Smoothing)。Smoothing同样是一种展宽。因此下列code只处理x值有规律的输入文件。但是很容易改写。


三个code仅仅是计算公式不一样，差别只是几行而已。不过，Voigt展宽后的结果最好。


附件1. gaussian-broadening.f90

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
!---------------------------
PROGRAM GaussianBroadening
!---------------------------
!
! Purpose : Perform a Gaussian broadening on a set of impulse signal data as input.
!           Sigma will be calculated from the hight of the impulses signal.
!           You must input the range and the step of the output data in the proper order.
!           No ouput file to be written, since it is designed to output to standard screen 
!           so as to work with gnuplot, as usage shown below.
!           if you want to get an output file with final data, uncommnent related 2 lines.
!           you would get 'g-broadened.dat'
!
! UsageEx : plot '< gbroad.x input.dat width Xmin Xmax Xstep < ' u 1:2 ...
!           in which, width     : FWHM, full width at half maximum, e.g., 3.0
!                     Xmin, Xmax: the range of X in output, e.g., 500 3000
!                     Xstep     : the increment of X in the output, e.g., 2
!
! Note    : all unit of X are the same with input, e.g., cm-1
!
! Author  : Xijun Wang, 2012.12.26
!
                                    
 implicit none
                                    
 integer, parameter     :: dp = kind(1.0d0)
 character(len=20)      :: arg, input
 character(len=100)     :: temp
 real(dp), allocatable  :: X0(:), Y0(:) ! 0 - input, no 0 - output
 integer,parameter      :: inputfile = 10, outputfile = 20
 real(dp) ::  X, Y, Xmin, Xmax, Ysum, Yavg, sigma,  pi ! sigma2 is sigma square
 integer  :: step, i, j, nline, stat  ! nline: number of lines
 logical  :: alive
                                    
 CALL getarg(1, arg)
 input = trim(arg)
                                    
 inquire(file=input, exist=alive)
 if( .not. alive) then
   write(*,*) input, "does not exist! "
   stop
 end if
                                    
 call getarg(2, arg)
 read(arg, *) Xmin
 call getarg(3, arg)
 read(arg, *) Xmax
 call getarg(4, arg)
 read(arg, *) step
 call getarg(5, arg)
 read(arg, *) sigma
                                    
 ! open and count number of lines in input file
 open(unit=inputfile, file=input, access="sequential", status="old")
                                    
 nline = 0
 do
     read(unit=inputfile, FMT=*, END=100) temp
     nline = nline + 1
 end do
 100 continue
                                    
 rewind(inputfile)
                                    
 ! allocate memory for arrays X0, Y0
 allocate(X0(1:nline), Y0(1:nline))
                                    
 ! read in data from input file
 do i = 1, nline
     read(unit=inputfile,FMT=*,iostat=stat) X0(i), Y0(i)
 end do
                                    
 ! Uncomment the following line if you want to output to file
 ! open(unit=outputfile,file="g-broadended.dat", status='replace', action='write')
                                    
 pi = 2.0 * acos(0.0_dp)
                                    
 X = Xmin
 do while(X .le. Xmax)
     Y = 0.0
     do i = 1, nline
       if( abs(X - X0(i)) .le. 3 * sigma ) then
         Y = Y + Y0(i)/(sigma*sqrt(2*pi)) * exp(-1.0*(X - X0(i))**2.0/(2.0*sigma*sigma) )
       end if
     end do
   ! uncomment the following line if you want to output to file
   ! write(unit=outputfile,fmt="(F6.2,1X,F8.6)") X, Y
     write(*,fmt="(F9.1,1X,F15.8)") X, Y
     X = X + step
 end do
                                    
 ! release memory 
 deallocate(X0, Y0)
                                    
 stop
END PROGRAM GaussianBroadening


附件2. l-broad.f90: Lorentzian broadening

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
!---------------------------
PROGRAM LorentzianBroadening
!---------------------------
!
! Purpose : Perform a Lorentzian broadening on a set of impulse signal data as input.
!           Sigma will be calculated from the hight of the impulses signal.
!           You must input the range and the step of the output data in the proper order.
!           No ouput file to be written, since it is designed to output to standard screen 
!           so as to work with gnuplot, as usage shown below.
!           if you want to get an output file with final data, uncommnent related 2 lines.
!           you would get 'g-broadened.dat'
!
! UsageEx : plot '< gbroad.x input.dat width Xmin Xmax Xstep < ' u 1:2 ...
!           in which, width     : FWHM, full width at half maximum, e.g., 3.0
!                     Xmin, Xmax: the range of X in output, e.g., 500 3000
!                     Xstep     : the increment of X in the output, e.g., 2
!
! Note    : all unit of X are the same with input, e.g., cm-1
!
! Author  : Xijun wang, 2012.12.26
!
                                   
 implicit none
                                   
 integer, parameter     :: dp = kind(1.0d0)
 character(len=20)      :: arg, input
 character(len=100)     :: temp
 real(dp), allocatable  :: X0(:), Y0(:) ! 0 - input, no 0 - output
 integer,parameter      :: inputfile = 10, outputfile = 20
 real(dp) ::  X, Y, Xmin, Xmax, Ysum, Yavg, sigma,  pi ! sigma2 is sigma square
 integer  :: step, i, j, nline, stat  ! nline: number of lines
 logical  :: alive
                                   
 CALL getarg(1, arg)
 input = trim(arg)
                                   
 inquire(file=input, exist=alive)
 if( .not. alive) then
   write(*,*) input, "does not exist! "
   stop
 end if
                                   
 call getarg(2, arg)
 read(arg, *) Xmin
 call getarg(3, arg)
 read(arg, *) Xmax
 call getarg(4, arg)
 read(arg, *) step
 call getarg(5, arg)
 read(arg, *) sigma
                                   
 ! open and count number of lines in input file
 open(unit=inputfile, file=input, access="sequential", status="old")
                                   
 nline = 0
 do
     read(unit=inputfile, FMT=*, END=100) temp
     nline = nline + 1
 end do
 100 continue
                                   
 rewind(inputfile)
                                   
 ! allocate memory for arrays X0, Y0
 allocate(X0(1:nline), Y0(1:nline))

喜欢

0

阅读┊ 收藏 ┊ 喜欢▼ ┊打印┊举报/Report
  

新浪BLOG意见反馈留言板　欢迎批评指正

新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 产品答疑

Copyright © 1996 - 2022 SINA Corporation, All Rights Reserved

新浪公司 版权所有