24小时热门版块排行榜    

查看: 2570  |  回复: 18
当前主题已经存档。

miaoskypop

铜虫 (小有名气)

[交流] 【分享】guassian使用经验

guassian使用经验
经过一段时间对GAUSSIAN的使用,发现不仅有助于加深对量子化学
的了解,而且可以巩固光谱学方面的理论知识。至少在研究光谱学理论
方面,GAUSSIAN是确实一个相当不错的东西。在使用中,也遇到很多问
题,而且在GAUSSIAN的帮助文件中没有涉及。这些问题有的是靠自己解
决的,有的是和别人探讨的,还有的是从网上或者文献中查到的。
    受本人研究方向所限,文中涉及的大多数都是处理分子激发态方面
的问题。如果在分子光谱学(包括红外、拉曼、可见紫外)方面的知识
有欠缺的话,建议先看一看哈里斯和伯特卢西写的《对称性与光谱学》
(胡玉才,戴寰译,高等教育出版社,1988),这是一本浅显易懂的速
成书籍。
    本人水平有限,对有些问题的理解可能有误,欢迎批评指正。
§1-1.内存与速度问题
    G98W默认的内存是128兆,我以前一直用64兆,计算机是P233,所
以奇慢无比。我后来发现,在计算之前把所有的杀毒程序都关闭(我的
计算机上装了三个),可以节省时间20%。后来又加了一条64兆内存,
时间就减少到原来的一半。如果把计算机换成P300,用128兆内存,只
需要从前五分之一的时间。而且大的内存支持更大的基组。所以计算机
配置可马虎不得。
§1-2.GAUSSIAN输入的一个注意事项
    输入的分子核坐标(包括键长和键角)必须是实型,不能是整型,
否则会报错。例如,键长2.0埃,可以写成r=2.或者r=2.0,不能写成
r=2。
§1-3.GAUSSIAN的缺陷
    对于含有重原子的分子来说,光谱项的自旋、轨道角动量已经失去
了意义,这个时候的光谱项应该用自旋-轨道耦合总角动量的各个分量
Ω表示。例如,一氯化铊分子的3Π电子态,在实验上观察到的是两个
分得很开的Ω分量:A0+态和B1态。但是GAUSSIAN仅能处理含有从H原子
到Cl原子的分子(因为不是重原子,所以没有必要),只能算是80年代
初期的水平。目前处理这类问题,除了自己编程,一般是使用MOLFDIR
程序。

[ Last edited by xuefei06 on 2009-2-10 at 18:27 ]
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

miaoskypop

铜虫 (小有名气)

(二)关于基组的设置
§2-1.基组
    基组有两种,一种是全电子基组,另一种是价电子基组。价电子基
组对内层的电子用包含了相对论效应的赝势(缩写PP,也称为模型势,
有效核势,缩写ECP)进行近似。以铟原子为例,可以输入下面的命令
来观察铟原子3-21G基组的形式:
# 3-21G GFPrint

In

0,2
In
输出结果为:
3-21G (6D, 7F)
S 3 1.00                                      N=1层,s轨道
.1221454700E+05 .6124760000E-01
.1848913600E+04 .3676754000E+00
.4063683300E+03 .6901359000E+00
SP 3 1.00                                     N=2层,s和p轨道
.5504422550E+03 -.1127094330E+00 .1523702990E+00
.1197743540E+03 .8344349530E-01 .6096507530E+00
.3866926950E+02 .9696880360E+00 .3970249500E+00
SP 3 1.00                                     N=3层,s和p轨道
.4702931320E+02 -.2758954250E+00 -.1408484750E+00
.2249642350E+02 .5977348150E-01 .5290866940E+00
.6697116970E+01 .1082147540E+01 .6620681110E+00
D 3 1.00                                      N=3层,d轨道
.1021735600E+03 .1205559000E+00
.2839463200E+02 .4884976000E+00
.8924804500E+01 .5850190000E+00
SP 3 1.00                                     N=4层,s和p轨道
.6572360380E+01 .4284830560E+00 .1091304620E-01
.2502157560E+01 -.4633643610E+00 .5036758870E+00
.9420245940E+00 -.8219679320E+00 .5581808650E+00
D 3 1.00                                      N=4层,d轨道
.4535363700E+01 .2508574000E+00
.1537148100E+01 .5693113000E+00
.4994922600E+00 .3840635000E+00
            以上都是内层电子轨道,每个轨道用三个高斯型函数拟合,
            也就是3-21G中的3。
SP 2 1.00                                     N=5层,s和p轨道
.1001221380E+01 -.4364171950E+00 -.2316333510E-01
.1659704190E+00 .1189893480E+01 -.9903308880E+00
            这是价电子轨道,每个轨道用二个高斯型函数拟合,也就
            是3-21G中的2。
SP 1 1.00                                     N=6层,s和p轨道
.5433974090E-01 .1000000000E+01 .1000000000E+01
            这是空轨道,每个轨道用一个高斯型函数拟合,也就是
            3-21G中的1。
可见3-21G是全电子基组。再看看铟原子LanL2DZ基组的结果:
LANL2DZ (5D, 7F)
S 2 1.00
.4915000000E+00 -.4241868100E+01
.3404000000E+00 .4478482600E+01
S 1 1.00
.7740000000E-01 .1000000000E+01
P 2 1.00
.9755000000E+00 -.1226473000E+00
.1550000000E+00 .1041757100E+01
P 1 1.00
.4740000000E-01 .1000000000E+01
这些都是价电子轨道。在ECP上加上了双zeta(DZ)基组。所以LanL2DZ
是价电子基组。
    在GAUSSIAN中,CEP也是很重要的价电子基组,并且98版中对元素
的适用范围扩大到了86号元素Rn。需要注意的是,超过了Ar元素后,每
一种元素的CEP-4G,CEP-31G,CEP-121G是等价的(和加了极化*的CEP
还是有区别的)。SDD适用于除了Fr和Ra以外的所有元素,但是结果的
准确性并不太好。另外,在下面的两个网址,有几乎整个周期表元素的
赝势下载,不过我想不出那些短寿命的超铀元素有什么用:
http://www.mapr.ucl.ac.be/ABINIT
http://www.theochem.uni-stuttgart.de
    需要注意的是有些方法不需要输入基组,像各种半经验方法,G1,
G2等等。
§2-2.基组的选择
    《中国科学》B30(5),2000,419页的文章给出了一种选择基组
的方案,认为基函数的选择应该根据体系中不同原子的性质及实际的化
学环境:
    1。描述一般体系时可根据该原子在元素周期表的位置从左到右依
次增加基函数用量。
    2。对核负电原子,使用更多的基函数以及适当的极化函数或弥散
函数。对核正电原子,基函数用量可适当减少。
    3。对正常饱和共价键原子不需要增加极化或弥散函数,对氢键、
弱相互作用体系、官能团、零价或低价态金属原子等敏感体系,要增加
极化或弥散函数。
    这样就可以用适中的基组和可承受的计算量的到相当可靠的计算结
果。方案适用于HF,MPn,DFT等计算中。
    在GAUSSIAN中有关的关键字是ExtraBasis和Gen。二者都是在route
section指定基组的基础上,再加上附加的基函数。不同的是Gen可以自
己定义基组。实际上可并不那么简单。下面是一个使用ExtraBasis的例
子:
# HF/3-21G ExtraBasis
AgCl
molecule specification
Cl 0
6-311G(d,p)
****
    这里对分子中所有的原子都使用3-21G,但是对Cl原子使用更大的
基组6-311G(d,p)。
    如果计算过程中报错,一般是由于下面两点原因造成的:
1。并不是所有的基组都可以用(哪些基组适用于哪些原子需要去参考
帮助文件),即使是分别适用于不同原子的各个基组,放在一起也不
一定就行,可能有个匹配的问题。所以需要多试。
2。基组设置的顺序问题。还是看上面的例子,我认为route section
一行中设置的基组应该适用于分子中所有的原子(或者应该是比较小
的基组)——虽然在计算中并不对所有的原子都使用这个基组,而在
分子坐标下面定义的基组只要对该原子适用就可以了(一般来说是大
一些的基组)。比如,如果把上面的输入变成:
# HF/6-311G(d,p) ExtraBasis
AgCl
molecule specification
Ag 0
3-21G
****
虽然都是对Ag原子用3-21G基组,对Cl原子用6-311G(d,p)基组,但是可
能会报错。当然实际情况并不总是如此,Ag和Cl的基组哪个在前哪个在
后可能并没有太大关系,但是有时候注意基组顺序确实可以解决问题。
2楼2009-02-10 13:28:51
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

miaoskypop

铜虫 (小有名气)

(三)用GAUSSIAN处理激发态
§3-1.输入文件中的电子态多重度
    输入文件中的多重度,指的是具有这种多重度的最低能量的电子态,
不一定就是基态,不过一般输入文件都用基态。激发态的多重度是在输
入文件的route section行控制的,如:singlets,triplets。虽然在
GAUSSIAN使用帮助中没有明说,但是从给出的例子中可以体会出来。另
外,在HyperChem的帮助文件中有类似的说明,虽然是针对HyperChem的,
但我想对GAUSSIAN一样适用。
§3-2.GAUSSIAN中研究激发态的方法
    GAUSSIAN可以用CIS或者ZINDO方法计算激发态,CIS能够研究的分
子种类与使用的基组有关,GAUSSIAN中的CIS实际是对分子激发态的零
级处理。有些从头计算程序包含计算激发态的一级组态相关(FOCI)和
二级组态相关(SOCI),我不清楚是不是和CIS相对应。但是,从结果
来看,FOCI和SOCI并没有比CIS改善多少。ZINDO方法适用范围窄,只能
用于周期表前48种元素构成的分子,这里不谈。98版还新增加了含时
(TD)方法处理激发态(包括含时Hatree-Fock和含时密度泛函)。
    我们一般计算的激发态称为Low-lying Excited States,也就是价
态。更高的电子态是Rydberg态。从头算法和其他半经验方法对空轨道,
特别是较高空轨道的计算结果不好,因而得到的Rydberg态的结果缺乏
明确的物理意义。
§3-3.含时方法(TD)输入方式
    帮助文件对TDDFT这么重要的方法说得真是太简单了,让人去参考
CIS的输入方法,连个例子都没有。CIS的输入方法是:
# CIS(Triplets,NStates=10)/3-21G Test
那就用B3LYP泛函试试看吧:
# TDB3LYP(Triplets,NStates=10)/3-21G Test
但是不行。其实正确的输入方法是:
# B3LYP/3-21G TD(Triplets,NStates=10)Test
这个问题浪费了我很长时间。为了表示它的重要,单独分为一节。
§3-4.CIS,TDHF,TDDFT方法激发能比较
    一般来说,CIS和TDDFT方法结果的准确度是比较差的,TDDFT的结
果是最好的。而各种TDDFT相比,B3LYP,B3P86,MPW1PW91是最好的(
依具体的分子和使用的基组而不同)。可以参考文献:
J. Chem. Phys., 111(7), 1999, 2889
Chem. Phys. Lett., 297(1-2), 1998, 60
J. Chem. Phys., 109(23), 1998, 10180
    这些文章都是计算C,H,O组成的有机物在基态平衡位置的垂直激
发能。我对含有重元素In的分子激发态进行计算,发现如果选好基组,
有些TDDFT结果和实验值的差距甚至小于500个波数,而CIS的结果则要
差5000多个波数。所以现在使用TDDFT方法的比较多。另外还有人用
TDDFT方法计算DNA这样的大分子。
§3-5.GAUSSIAN处理激发态的BUG
    1。用94版进行CIS计算,50-50选项不支持LanL2DZ和3-21G*基组。
98版则不会出现这个问题。这个问题可能还跟计算的分子有关。
    2。同时计算单重激发态和三重激发态(也就是使用50-50这个选项
),有时候会出现致命错误:某些电子态在某些位置的能量有不规则起
伏,而分开计算(使用Singlets或者Triplets选项)则不会出现这个问
题。我在计算中曾经遇见过两次。下面这个例子是从网上看到的,我想
这位德国的大哥也遇见了同样的问题。
    他对同一个任务执行了三种计算:
1) #Bp86/6-31G* td=(nstates=5,50-50) scf=direct density
2) #Bp86/6-31G* td=(triplets,nstates=1) scf=direct density
   --link1--
   #Bp86/6-31G* td=(nstates=10) scf=direct guess=read geom=check
   --link1--
   #Bp86/6-31G* td=(triplets,nstates=10) scf=direct guess=read
   geom=check
3) #Bp86/6-31G* td=(nstates=10) scf=direct
   #Bp86/6-31G* td=(triplets,nstates=10) scf=direct density
    第一种计算产生10个态(5S和5T),后两种产生20个。关键字
Density对结果没有影响。结果是:
                 1)        2)        3)
Triplet-B1U    3.1077    3.1077    3.1077
Triplet-B3G    3.9568    3.9568    3.9568
Triplet-AG     4.1612    4.1612    4.1612
Triplet-B2U    4.1769    4.1769    4.1769
Triplet-B3G    4.2101    4.2101    4.2101
Singlet-B3G    4.3609    4.3609    4.3609
Singlet-B1U    4.3788   >4.4864   >4.4864
Singlet-B2U    4.4864   >4.6409   >4.6409
Triplet-AG     4.5752   >4.7144   >4.7144
Singlet-B3G    4.6409   >  --     >5.0010
    通过比较发现,对前六个态,三种处理的结果相同,但是后面的态,
使用了50-50的方法就和另外两种方法的结果不一样了。
    我个人认为,尽量不要使用50-50这个选项。推荐单重态、三重态
分开计算,虽谈多花一倍时间,但是结果保险。
§3-6.画激发态的势能曲线
    最笨的方法莫过于算完一个点后,更改输入文件的坐标,计算下一
个点了。这样费时费力。在优化基态几何构型的时候,可以使用Scan关
键字进行扫描,实际上Scan对激发态一样适用。这样开上一晚上计算机,
到早上所有激发态的势能曲线数据就全有了。例如计算从R=0.6埃到1.2
埃氢分子的10个单重激发态:
# RCIS(NStates=10)/6-311G Scan
H2 Excited States
0,1
H1
H2 H1 r
r 0.6 60 0.01
    需要注意三点:
    1。结果给的是垂直激发能,必须换算成相对于基态平衡位置的能
量。
    2。扫面范围不要离基态平衡位置太远,一般从0.75Re到2Re就可以
了,太远了会报错。
    3。某些理论和某些基组的搭配可能不支持Scan,这时候就只好手
工控制了。
    最后是数据处理,这一步是最枯燥的了。数据处理需要用到Excel、
Oringin、Matlab。这些软件用不着钻研很深,知道几个最基本的命令
就足够用了。
§3-7.高对称性的分子使用CIS方法
    默认的CIS计算求出最低的三个激发态。以它为例,GAUSSIAN在初
始的猜测中选择电子态数目两倍的向量(也就是6个)进行迭代,直到
有三个收敛。对于大多数分子体系来说,这种默认的向量数目足够了。
但是对于高对称性的分子,则需要特殊的方案,否则不会得到比较高的
收敛态。方法是增加NStates的数量,以增加初始向量的数量。推荐
NStates值为最大阿贝尔点群操作的数量,也就是输出文件在Standard
orientation之前的Symmetry部分(在输入文件中使用#就可以得到,#T
不行)。这是一个计算苯分子的例子:
Stoichiometry   C6H6
Framework group    D6H[3C2'(HC,CH)]
Deg. of freedom   2
Full point group   D2H   NOp     8
Largest concise Abelian subgroup   D2  NOp   4
           Standard orientation
3楼2009-02-10 13:29:49
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

miaoskypop

铜虫 (小有名气)

(四)求分子的光谱常数
§4-1.求离解能
    首先必须注意,原子化能和离解能不是同一个概念。帮助文件和
Exploring Chemistry with Electronic Structure Methods一书中的
解释并不确切。计算原子化能和离解能是不同的两个过程。以三原子分
子GaCl2为例,计算离解能的过程是:
GaCl2(X2A1)---->Cl(2P)+GaCl(X1Σ+)
过程中的能量变化称为离解能。很明显,离解能还可以细分为第一步离
解能、第二步离解能......,上面的是计算第一步离解能,计算第二步
离解能的过程是:
GaCl(X1Σ+)---->Cl(2P)+Ga(2P)
而计算原子化能的过程是:
GaCl2(X2A1)---->Cl(2P)+2Ga(2P)
可见还是有区别的。不过对于双原子分子来说,这两个概念是相同的。
    第二个需要注意的是,离解能一般用符号D0或者De表示(前者不包
括零点能,后者包括),但是在分子光谱中,还有一组D0和De。为避免
混淆,这里需要作特别说明。
    第一组D0和De,也就是离解能,数值比较大,数量级一般在0.1到
1eV。这个参数多出现在真空-紫外、振动光谱中。有一个公式:
              ωe^2
    ωeχe = -------       (^是指数符号)
               4De
其中的De就是离解能。
    第二组D0和De称为离心畸变常数,数值很小,作为转动光谱中的微
扰项。如:
                         h^3
    De = -----------------------------------
          128(π^6)(μ^3)(c^3)(re^6)(ωe^2)
其中的De就是离心畸变常数。
    GAUSSIAN中有个求PH2分子的原子化能的例子,可以作为求离解能
的参考。也就是文件e7-01。结果是
原子化能 = E(P)(Hartree)+2E(H)(Hartree)-E(PH2)(Hartree)
           -ZPE(Hartree)
         = (-341.25930)+(-0.50027)*2-(-342.50942)-(0.01322)
           (Hartree)             其中的零点能ZPE经过了矫正
         = 148.3 Kcal/mol
和实验值144.7 Kcal/mol的差距相当小了。

§4-2.求电离能
    电离能缩写为I.P.。在GAUSSIAN中有两种求I.P.的方法。
法一:见文件e7-03,这个例子求PH2分子的电离能。结果为:
         E               ZPE
PH2+   -342.14416      0.01347 (ZPE经过矫正)
PH2    -342.50942      0.01322 (这是上面求原子化能的结果)
I.P. = (E(PH2+)+ZPE(PH2+))-(E(PH2)+ZPE(PH2))
     = (-342.14416+0.01347)-(-342.50942+0.01322)(Hartree)
     = 0.36551(Hartree)
     = 9.95(eV)
法二:OVGF方法。帮助文件里有详细的使用说明。

§4-3.求激发态的光谱学参数Te,ωe,Re
有两种方法。
    1。CIS可以用结构优化的方法。Starl的文章里说得很详细了,这
里不再细讲。仅仅需要补充的是,如何强制保留设置后的波函:
    在使用Guess=Alter设置波函之后,程序总会回到设置前的波函。
这个时候需要使用SCF关键字中的Symm选项,强制程序保留设置以后的
波函。使用方式是在route section部分添加SCF=Symm。
    2。势能曲线扫描法。
    TDDFT方法不能计算激发态的波函,所以不能用来研究激发态的平
衡构型和分子属性。上面的结构优化方法用在TDDFT中就不适合了。但
是可以用势能曲线扫描的方法。如何扫描激发态的势能曲线见前面的§
3-6节。
    以双原子分子为例,有了势能曲线以后,就可以找到势能最低点Re
(为了精确,这一区域的点要取密一些),Re位置的能量相对于基态平
衡位置能量的差值就是电子态高度Te。势阱的深度就是离解能De。
    这些都简单,关键是求振动参数ωe。求ωe也有两种方法。一种是
用Morse势函数拟和。Morse势函数形式为:
    V(R) = Te+De{1-EXP[-α(R-Re)]}^2    (对基态Te是0)
Te、De、Re是已知的,拟和出了α后,用公式
    α = ωe ((2(π^2)c*μ)/(h*De))^0.5    (高斯单位制)
就可以得到ωe。
    这种方法太麻烦,我们不这样做。实际只要把扫描出来的Morse势
函数对R求二阶导数,在Re处的导数值f2和ωe的关系有:
          4(π^2)(ωe^2)c*μ
    f2 = --------------------
              h*(10^18)
其中c、h、μ是国际单位制,ωe的单位是波数cm^-1。另外在扫描的势
函数中,De的单位是波数cm^-1,R和Re的单位是埃。否则上面的f2需要
乘上个系数。
    求f2用Origin软件就可以做。注意要取Re位置的f2。有了f2,就可
以算出ωe,省去了计算De和拟合α的步骤,避免引入更大的误差。

§4-4.如何获得光谱项的对称性
    如果用GAUSSIAN计算比较轻的分子,一般都能给出正确的激发态光
谱项对称性,但是对含有重原子的分子,可能得不到对称性。这个时候
就需要根据轨道跃迁的的信息获得对称性。
    下面是用CIS方法LanL2DZ基组计算InCl分子的例子,输出文件给出
InCl的电子轨道信息:
Occupied  (SG) (SG) (PI) (PI) (SG)
Virtual   (PI) (PI) (SG) (PI) (PI) (SG) (PI) (PI) (SG) (SG)
          (SG)
这里结果只给出前4个单重激发态:
CIS wavefunction symmetry could not be determined.
Excited State   1:   Singlet-?Sym    4.9226 eV  251.87 nm
  f=0.3243
       5 ->  6          .43587
       5 ->  7          .54909
This state for optimization and/or second-order correction.
Copying the Cisingles density for this state as the 1-particle
RhoCI density.
CIS wavefunction symmetry could not be determined.
Excited State   2:   Singlet-?Sym    4.9226 eV  251.87 nm
  f=0.3243
       5 ->  6          .54909
       5 ->  7         -.43587
Excited State   3:   Singlet-SG     6.7609 eV  183.38 nm
  f=0.2830
3 ->  6         -.28038
       4 ->  7         -.28038
       5 ->  8          .57818
CIS wavefunction symmetry could not be determined.
Excited State   4:   Singlet-?Sym    6.8241 eV  181.68 nm
  f=0.0000
       3 ->  7          .49690
       4 ->  6         -.49690
    由于In原子是重原子,所以大多数电子态并没有给出对称性(第三
个激发态除外)。我们发现第一激发态是从第五个轨道跃迁到第六、七
个轨道,而第六、七两个轨道是简并的(PI、DELT轨道都是二重简并的
)。也就是说第一激发态的电子组态是:
InCl: (内层电子)(1σ)2 (2σ)2 (1π)4 (3σ)1 (2π)1
    这个组态耦合出来的单重态光谱项有一种(当然还有一个三重的):
1Π,所以就是1Π态。再往下看,第二激发态和第一激发态的电子组态、
能量完全相同,说明这是一个简并态,而实际上除了Σ态以外的态都是
二度简并的。所以第一激发态和第二激发态实际都是InCl的第一激发态
1Π态。
    第三个激发态给出对称性了。如果没给,方法和上面的一样,先找
到它的电子组态,它是由两个组态叠加而成的:
InCl: (内层电子)(1σ)2 (2σ)2 (1π)3 (3σ)2 (2π)1
InCl: (内层电子)(1σ)2 (2σ)2 (1π)4 (3σ)1 (4σ)1
    这个组态耦合出来的单重态光谱项有一种(当然还有一个三重的):
1Π,所以就是1Π态。再往下看,第二激发态和第一激发态的电子组态、
能量完全相同,说明这是一个简并态,而实际上除了Σ态以外的态都是
二度简并的。所以第一激发态和第二激发态实际都是InCl的第一激发态
1Π态。
    第一个组态可以耦合出的单重态光谱项有:1Δ,1Σ+,1Σ-。因
为这不是一个简并态,所以只能是1Σ+和1Σ-其中之一。这就需要进一
步知道+、-对称性。方法是看跃迁系数:
       3 ->  6         -.28038
       4 ->  7         -.28038
由于正负号相同,所以是1Σ+态。那么第四激发态必然是1Σ-了。
    另外还可以看第二个组态耦合出来的光谱项,只有1Σ+,因此第三
个激发态必然是1Σ+。两个结果一致。
4楼2009-02-10 13:30:24
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

sizj

金虫 (正式写手)


1

hao ding
5楼2009-02-14 09:28:47
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

erylingjet

铁杆木虫 (著名写手)

0.5

支持鼓励,呵呵

不错
6楼2009-02-14 16:47:41
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhangwc5020

银虫 (小有名气)

评价

感激五内
7楼2009-02-14 17:15:13
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lxl8506

铜虫 (正式写手)

辛苦了,非常感谢
8楼2009-02-14 19:09:14
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

diytkib891

好人啊!
9楼2009-02-28 00:13:14
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

a349857460

木虫 (著名写手)

10楼2009-02-28 12:28:44
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 miaoskypop 的主题更新
普通表情 高级回复 (可上传附件)
信息提示
请填处理意见