コンピュータによる赤外スペクトルからの構造推定

田辺 和俊, 松本 高利, 伊藤 祥司, 上坂 博亨, 都築 誠二, 田村 禎夫, 佐伯 慎之助, 小野 修一郎

1 緒 言

2 帰納的方法−スペクトルデータベース

Table 1. Average ranks of correct answers

	Scoring equation	Average rank
A	SIh / SIs × 100	27.1
B	[SIh - SIn] / SIs × 100	10.4
C	[SIh - SIn'] / SIs × 100	6.3
D	[SIh - S|Is - Ih|] / SIs × 100	17.4
E	[SIh - SIn - S|Is - Ih|] / SIs × 100	7.5
F	[SIh - SIn' - S|Is - Ih|] / SIs × 100	5.6
G	[SIs - SIn] / SIs × 100	21.5
H	[SIs - SIn'] / SIs × 100	4.5
I	[SIs - S|Is - Ih|] / SIs × 100	14.0
J	[SIs - SIn - S|Is - Ih|] / SIs × 100	9.8
K	[SIs - SIn' - S|Is - Ih|] / SIs × 100	6.0
L	[NsSIsIh - (SIs)(SIh)] / [NsSIs2 - (SIs)2]1/2 / [NsSIh2 - (SIh)2]1/2 × 100	37.9
M	[NfSIsIh - (SIs)(SIf)] / [NfSIs2 - (SIs)2]1/2 / [NfSIf2 - (SIf)2]1/2 × 100	1.7

I: peak intensity (absorbance), N: number of peaks, s: sample (input), f: file, h: hit peaks among file peaks, n: non-hit peaks among file peaks, n': non-hit peaks stronger than the weakest sample peak.

また、最近では化学反応プロセスのモニタリングが関心を集めており [19]、そのためには反応途中の不安定化学種を同定する必要がある。しかし、スペクトルデータベースは元来、安定化合物について構築されることが多いので、不安定化学種のような未知化合物の関与する問題には無力である。
しかし、最大の問題点はデータベースの規模の問題である。データベース中の標準スペクトルと照合して正答が出るためにはその化合物のスペクトルがデータベースに集積されていることが不可欠である。データベースに集積されていない化合物はいかにしても正答が得られないからである。したがって、正答率を上げるためには可能な限りの化合物のスペクトルを集積してデータベースの規模をできるだけ大きくすることが必要となる。
ところが、これまでに化学者が論文で扱ってきた化合物の総数は 2,000 万種類以上にも達しており、その中で現実社会で流通している化合物だけでも数十万種類以上になるといわれている。それに対して、これまでに構築されているスペクトルデータベースは最大規模のものでも数万化合物程度であり、正答率の点からは十分な規模ではない。我々の SDBS も現在約 3 万化合物にとどまっており、大規模なスペクトルデータベースの構築はきわめて困難な問題である。

3 経験的手法

3. 1 相関表

3. 2 パターン認識

Table 2. Number of compounds which contain the functional group (Nc), number of compounds which contain it but are estimated not to contain (Ncn), number of compounds which do not contain it (Nn), number of compounds which do not contain it but are estimated to contain (Nnc), and error rate % (E)

Functional group	Nc	Ncn	Nn	Nnc	E
C=O	4570	342	5046	389	7.6
C=C	2135	967	7481	1364	24.2
C≡N	223	68	9393	150	2.3
C-O-C	1264	459	8352	752	12.6
Benzene ring	4572	698	5044	705	14.6
NO2	485	109	9131	212	3.3

しかし、この正答率は教師データに対する学習分類率であり、未知データに対する予測率はこの数値よりかなり低下すると思われる。なぜならば、赤外スペクトルパターンと構造との関係は本来、線形ではなく非線形であり、線形学習機械で取り扱うのは限界があると考えられるからである。

3. 3 ニューラルネットワーク

3. 3. 1 入力データ

Table 3. Main functional groups, sub-functional groups, structures, numbers of substances that contain those groups (N), and identification accuracies (%) for present (P) and absent (A) groups

Main functional group	Sub-functional group	Structure	N	P	A
AMIDE			1024	83.0	80.7
	AMIDE1	CO-NH2	136	91.2	85.3
	AMIDE2	CO-NHR	608	81.9	82.9
	AMIDE3	CO-NR2	214	72.9	75.7
	AMIDE4	beta-lactam	12	83.3	33.3
	AMIDE5	gamma-lactam	10	80.0	80.0
	AMIDE6	delta-lactam	2	100.0	100.0
	AMIDE7	R-O-CO-NHR	4	100.0	0.0
	AMIDE8	NHR-CO-NHR	132	89.4	72.7
AMINO			682	80.6	81.2
	AMINO1	R-NH2	232	84.5	86.2
	AMINO2	Ar-NH2	188	88.3	81.9
	AMINO3	R-NH-R	82	75.6	75.6
	AMINO4	R-NH-Ar	50	80.0	76.0
	AMINO5	Ar-NH-Ar	16	100.0	62.5
	AMINO6	NR3	152	61.8	63.2
BENZ			886	72.0	75.6
	BENZ1	mono-substituted benzene ring	262	87.0	87.8
	BENZ2	o-di-substituted benzene ring	256	70.3	82.8
	BENZ3	m-di-substituted benzene ring	72	86.1	66.7
	BENZ4	p-di-substituted benzene ring	188	78.7	87.2
	BENZ5	1,2,4-tri-substituted benzene ring	124	66.1	75.8
	BENZ6	1,2,3-tri-substituted benzene ring	76	76.3	68.4
	BENZ7	1,3,5-tri-substituted benzene ring	10	20.0	60.0
CDC			1304	63.8	66.6
	CDC1	CH2=CHR	186	79.6	86.0
	CDC2	CH2=CR2	80	77.5	75.0
	CDC3	CHR=CHR(cis)	462	66.2	64.5
	CDC4	CHR=CHR(trans)	84	83.3	73.8
	CDC5	CHR=CR2	528	60.6	65.9
	CDC6	CR2=CR2	148	44.6	74.3
	CDC7	C=C=C	2	100.0	100.0
	CDC8	C≡C	86	90.7	76.7
CHO			154	72.7	80.5
	CHO1	R-CHO	34	88.2	70.6
	CHO2	C=C-CHO	30	100.0	86.7
	CHO3	Ar-CHO	90	66.7	77.8
COOH			742	91.1	86.5
	COOH1	R-COOH	302	74.2	82.8
	COOH2	CX-COOH	24	91.7	83.3
	COOH3	C=C-COOH	62	90.3	71.0
	COOH4	Ar-COOH	174	88.5	81.6
	COOH5	amino acid	182	96.7	95.6
CTN			670	74.3	98.2
	CTN1	R-C≡N	292	63.0	95.9
	CTN2	C=C-C≡N	162	86.4	93.8
	CTN3	Ar-C≡N	220	63.6	96.4
ESTER			518	89.2	89.2
	ESTER1	R-COO-R	304	86.8	95.4
	ESTER2	C=C-COO-R	76	94.7	63.2
	ESTER3	R-COO-C=C	12	83.3	100.0
	ESTER4	C=C-COO-C=C	4	100.0	50.0
	ESTER5	Ar-COO-R	84	90.5	78.6
	ESTER6	R-COO-Ar	24	91.7	100.0
	ESTER7	Ar-COO-Ar	6	66.7	100.0
	ESTER8	C=C-COO-Ar	6	100.0	0.0
ETHER			486	68.3	77.8
	ETHER1	R-O-R	236	89.8	83.9
	ETHER2	R-O-C=C	14	85.7	85.7
	ETHER3	R-O-Ar	194	74.2	91.8
	ETHER4	C=C-O-C=C	28	64.3	78.6
	ETHER5	Ar-O-Ar	12	33.3	100.0
	ETHER6	C=C-O-Ar	12	83.3	66.7
HALO			906	57.2	68.0
	HALO1	C-F	96	52.1	77.1
	HALO2	C-Cl	548	55.5	65.0
	HALO3	C-Br	252	43.7	66.7
	HALO4	C-I	64	37.5	59.4
KETONE			632	71.2	76.9
	KETONE1	R-CO-R	154	72.7	92.2
	KETONE2	R-CO-C=C	136	64.7	75.0
	KETONE3	R-CO-CX	16	87.5	37.5
	KETONE4	R-CO-CO-R	12	83.3	66.7
	KETONE5	5-membered ring ketone	38	89.5	89.5
	KETONE6	R-CO-Ar	158	62.0	84.8
	KETONE7	Ar-CO-Ar	42	85.7	95.2
	KETONE8	p-quinone	78	79.5	79.5
NITRO			1028	90.7	94.4
	NITRO1	R-NO2	34	94.1	100.0
	NITRO2	C=C-NO2	36	44.4	77.8
	NITRO3	Ar-NO2	948	87.8	92.4
	NITRO4	>N-NO2	10	40.0	80.0
	NITRO5	O-NO2	2	100.0	0.0
OH			772	76.4	85.0
	OH1	CH2R-OH	248	83.1	87.9
	OH2	CHR2-OH	248	73.4	90.3
	OH3	CR3-OH	48	79.2	62.5
	OH4	Ar-OH	322	68.3	84.5
PHOS			482	75.5	80.5
	PHOS1	H-P(=O)(OR)2	8	75.0	50.0
	PHOS2	R-P(=O)(OR)2	226	82.3	85.0
	PHOS3	P(OR)3	4	100.0	100.0
	PHOS4	others	244	68.9	81.1
SULF			926	68.3	57.2
	SULF1	SO3H	86	81.4	76.7
	SULF2	SO3R	34	76.5	88.2
	SULF3	SO2-O-	6	33.3	100.0
	SULF4	SO3	136	83.8	94.1
	SULF5	SH	108	66.7	46.3
	SULF6	-S-	326	58.9	59.5
	SULF7	SO2-	102	76.5	88.2
	SULF8	-S-S-	44	45.5	45.5
	SULF9	>C=S	108	55.6	61.1
	SULF10	-SO-O-	6	66.7	66.7
	SULF11	-SO-	4	100.0	100.0
Overall			5944	74.5	80.3

R: alkyl group，Ar: aromatic group

赤外スペクトルの波数範囲は通常4,000〜400 cm^-1であり、SDBS の赤外スペクトルもこの波数範囲のデータを集積している。一方、スペクトル中のピーク幅 (半値幅) はピークによって様々であるが、ニトリル基のように鋭いピークでは 2〜3 cm^-1 程度のものもある。したがって、このような鋭いピークをもつ赤外スペクトルパターンを忠実に再現するためには 1 cm^-1 程度の波数間隔が必要となり、ニューラルネットワークに入力する波数点は 3,600 にもなる。しかし、本研究で用いたニューラルネットワークのソフトウエアはユニット数の上限が 255 であり、全波数範囲のスペクトルパターンを 1 cm^-1 の波数間隔で入力することは不可能である。
この問題が赤外スペクトルなどのスペクトルパターンをニューラルネットワークで解析する場合の最大の問題であり、これまでの論文でも入力点を減らすためのスペクトルの特徴抽出の前処理方法として高速フーリエ変換 (FFT)、wavelet 変換、Hadamard 変換など色々な方法が試みられている。しかし、我々はスペクトル解析の専門家が構造推定を行う際に果たして FFT や wavelet 変換に類する特徴抽出操作を行っているのかどうかという点に疑問を抱き、むしろ専門家は実測の赤外スペクトルそのものに対してパターン認識を行っているのではないかと考えた。
そこで、本研究では専門家のそのような思考過程にできるだけ近い構造推定システムの開発を目指して、ニューラルネットワークに入力するデータとしてはこのような特徴抽出の変換操作をせず、実測のスペクトルパターンを入力することにした。したがって、ニューラルネットワークに入力する波数点が 255 以下になるように波数範囲と波数間隔を設定する必要がある。むしろ、入力点が 255 一杯になるよう波数範囲を広くとるよりも狭い波数範囲の方が正答率が高くなる可能性もある。
すなわち、Colthupの相関表などに示されているように各官能基のピークが現れる波数範囲はほぼ一定であり、入力データをそれらの波数範囲に限定することも考えられる。なぜならば、スペクトル解析の専門家が構造推定を行う場合、全波数範囲のスペクトルパターンに着目しているのではなく、各官能基ごとに特定の波数範囲のパターンに注目しているのではないかと考えられるからである。したがって、ニューラルネットワークに入力する場合も必要最小限度の波数範囲のスペクトルデータを入力する方が有効である可能性もある。
そこで、15種類の各官能基について、SDBSに集積されている赤外スペクトルデータからその官能基を含む化合物と含まない化合物をそれぞれ数百以上抽出した。次に、それぞれについて平均スペクトル P、A を計算し、両者の差スペクトル D を計算した。この方法で得られた各官能基ごとのスペクトル P、A、D を Figure 1 に示す。

Figure 1. Averaged spectra for substances that contain (thin) and that do not contain (dotted) functional groups and difference spectra (thick).

この差スペクトル D は各官能基に特有のスペクトルパターンであり、Colthupの相関表をコンピュータを用いてより定量的に表現したものと考えることができる。その吸光度が正の波数領域はその官能基の存在を、また吸光度が負の領域はその官能基の不在を示す特徴的な波数範囲であると考えられる。したがって、この差スペクトルの吸光度の絶対値が高い波数範囲をニューラルネットワークに入力すべき波数範囲と考えることができる。しかし、上に記したように、波数範囲は広い方が必ずしも適切とは言えない可能性がある。
そこで、この差スペクトルのデータを用いて波数範囲の影響を調べるために、官能基 BENZ について波数範囲を 4 種類設定し、正答率への影響を調べた。その結果、波数範囲が最も狭い場合に正答率が最大になり、波数範囲はできるだけ狭い方が適切であることが分かった。そこで、この方針に基づき、Figure 1 の差スペクトルのデータを参考にして各官能基の波数範囲と波数間隔を決定した。

3. 3. 2 ニューラルネットワークの構造

Figure 2. Two types of neural network structures

このニューラルネットワークを用いて Table 3 に示した合計 90 種類の官能基を推定するために本研究では Figure 2 に示す 2 段構造のニューラルネットワークシステムを用いた。すなわち、初段のニューラルネットワークでは15種類の主官能基の判定を行い、初段のニューラルネットワークの出力層に結合している次段のニューラルネットワークでは副官能基の判定を行った。
赤外スペクトルを入力して15種類の官能基の有無を推定するニューラルネットワークを構築する場合、Figure 2 に示す (a)、(b) の2種類の構造のニューラルネットワークが考えられる。構造 (a) は、初段では1個のニューラルネットワークを構築し、出力層のユニット数を15にしてそれぞれに主官能基を割り当て、全てのデータを一括して用いて学習を行う方法である。構造 (b) は、初段の15種類の主官能基ごとに1個のニューラルネットワークを構築し、出力層のユニット数は1にして、それぞれの主官能基ごとに学習を行う方法である。
両方の方法でニューラルネットワークを構築、学習し、テストを行った。その他の条件は以下の検討における最適条件を用いたが、この点は以下の検討においても同様である。幾つかの官能基について検討した結果、(b) の場合に正答率が高いことが分かった。そこで、以下の検討においてはこの構造を採用した。

3. 3. 3 ニューラルネットワークの学習とテスト

3. 3. 4 ニューラルネットワークによるテスト結果

3. 4 経験的方法の問題点

Figure 3. Infrared spectrum of 4-chloro-3,5-dimethylphenol

このように、ある官能基が分子内に存在するにも拘わらず、その特性吸収帯が赤外スペクトルに現れないために、ニューラルネットワークによる官能基推定の正答率が平均として 77.2 % という低い数値になったと考えられる。このような現象を与える原因としては

1. 分子の対称性によりその官能基を含む振動が赤外不活性、または振動双極子能率の値が小さくなり、赤外スペクトルにほとんど現れない。
   
2. その官能基による特性吸収帯が現れる波数領域で他の非常に強いピークと偶然重なったために、その官能基の存在が推定できない。
   
3. 赤外スペクトルは一般にスペクトルの中で最も強いピークの透過率が5% 程度になるように試料セルの厚さや溶液濃度を調整して測定するので、ある官能基の特性吸収帯がある程度の強度で現れたとしても、他の領域で極めて強いピークが存在すると、そのピークの強度で全体を規格化するために、他のピークの強度は低くなり、最初の官能基が存在しないという推定になる。
   

4 理論的方法

4. 1 方 法

4. 2 結果と考察

Table 4. Uncertainty class (UC), observed frequencies (cm^-1), means and standard deviations of frequency ratios

		UC	Obsd Freq	Mean	Std Dev			UC	Obsd Freq	Mean	Std Dev
CH4	n4	C	1306	1.0476	0.0475	CH2Cl2	n4	B	282	0.9981	0.0604
	n2	A	1534	1.0312	0.0367		n3	B	717	0.9776	0.0649
	n1	A	2917	1.0425	0.0252		n9	B	758	0.9646	0.0930
	n3	A	3019	1.0456	0.0250		n7	B	898	1.0150	0.0460
CH3D	n6	C	1155	1.0451	0.0455		n5	B	1153	1.0320	0.0523
	n3	C	1300	1.0485	0.0473		n8	B	1268	1.0384	0.0521
	n5	C	1471	1.0346	0.0382		n2	C	1467	1.0067	0.0459
	n2	A	2200	1.0393	0.0259		n1	B	2999	1.0486	0.0324
	n1	E	2945	1.0442	0.0259		n6	B	3040	1.0626	0.0328
	n4	B	3017	1.0463	0.0260	CHDCl2	n6	B	283	0.9902	0.0599
CH2D2	n4	C	1033	1.0340	0.0433		n5	B	692	0.9815	0.0621
	n7	C	1090	1.0419	0.0470		n9	B	738	0.9665	0.0868
	n9	C	1234	1.0435	0.0470		n4	C	778	1.0147	0.0481
	n5	C	1333	1.0283	0.0366		n8	A	890	1.0215	0.0568
	n3	C	1436	1.0353	0.0397		n7	A	1223	1.0342	0.0523
	n2	C	2202	1.0169	0.0253		n3	C	1282	1.0286	0.0470
	n8	C	2234	1.0461	0.0262		n2	B	2249	1.0418	0.0319
	n1	C	2974	1.0439	0.0258		n1	B	3024	1.0565	0.0322
	n6	C	3013	1.0476	0.0260	CD2Cl2	n4	C	282	0.9893	0.0605
CHD3	n3	C	1003	1.0400	0.0469		n3	B	687	0.9807	0.0642
	n6	C	1036	1.0330	0.0429		n7	D	712	1.0106	0.0467
	n5	C	1291	1.0349	0.0400		n9	B	727	0.9700	0.0857
	n2	C	2142	1.0245	0.0255		n5	C	826	1.0243	0.0523
	n4	C	2263	1.0327	0.0259		n8	B	957	1.0321	0.0593
	n1	C	2993	1.0462	0.0259		n2	D	1052	1.0297	0.0478
CD4	n4	B	996	1.0390	0.0468		n1	B	2205	1.0344	0.0325
	n2	B	1092	1.0251	0.0366		n6	C	2304	1.0457	0.0321
	n1	B	2109	1.0201	0.0254	CHCl3	n6	B	261	0.9957	0.0608
	n3	A	2259	1.0345	0.0260		n3	C	363	1.0018	0.0605
CH3Cl	n3	A	732	0.9721	0.0701		n2	B	680	0.9730	0.0597
	n6	A	1017	1.0249	0.0459		n5	B	774	0.9588	0.0987
	n2	A	1355	1.0385	0.0479		n4	B	1220	1.0290	0.0563
	n5	A	1452	1.0344	0.0424		n1	B	3034	1.0570	0.0340
	n1	E	2937	1.0523	0.0299	CDCl3	n6	C	262	0.9888	0.0606
	n4	B	3039	1.0521	0.0301		n3	C	369	0.9805	0.0590
CD3Cl	n3	A	701	0.9736	0.0662		n2	B	659	0.9753	0.0604
	n6	A	768	1.0169	0.0458		n5	B	749	0.9626	0.0901
	n2	A	1029	1.0295	0.0494		n4	B	914	1.0177	0.0644
	n5	A	1060	1.0269	0.0420		n1	C	2266	1.0430	0.0334
	n1	A	2160	1.0238	0.0292	CCl4	n2	C	217	1.0030	0.0645
	n4	A	2283	1.0412	0.0298		n4	C	314	0.9997	0.0587
							n1	C	459	0.9753	0.0720
							n3	E	776	0.9717	0.1060

Uncertainty class (UC) taken from ref [46].

84個の振動に共通して、振動数の高い振動ほど振動数比も高いという傾向が認められる。分子軌道法による振動数の計算値を補正するscale factorとして種々の値がこれまでに提案されている [48 - 52]。従来のscale factorは振動数に関わらず一定値を用いるのに対し、最近、多数の分子についてのDFT法による振動数の計算結果から、振動数と1次式の関係のscale factorを用いる linear scaling 法が提案された [51, 52]。我々の結果はRHFとMP2を含む計算法全般についてこの方法と同様な傾向があることを示している。
振動形ごとに振動数比をまとめると、C-H 伸縮振動では1.0425-1.0626、C-D 伸縮振動では1.0169-1.0461、C-H およびC-D 変角振動では1.0067-1.0485、C-Cl 伸縮振動では0.9588-0.9815、C-Cl 変角振動では0.9805-1.0030となり、それぞれかなり狭い範囲に収まる。特に興味深いのは塩素原子が関与する振動の振動数比であり、振動数が高い C-Cl 伸縮振動よりも振動数が低い C-Cl 変角振動の振動数比が大きくなった。
この結果は振動数が高いほど振動数比も高い関係にあるとする上記のlinear scaling法とは逆の結果である。我々の結果は塩素化メタン分子13種の84個の振動についての結果であるのに対し、linear scaling法の論文でははるかに多種多様の分子のはるかに多数の振動を取り扱っており、単純な比較はできない。しかし、振動数が高い C-Cl 伸縮振動よりも振動数が低い C-Cl 変角振動の振動数比が大きいという我々の結果は、塩素を含む分子の振動の非調和性について有益な情報を提供していると思われる。
一方、標準偏差は C-Cl 伸縮振動が最大で、C-Cl変角振動や C-H または C-D 変角振動がそれに次ぎ、C-H または C-D 伸縮振動で最小、という傾向が見られる。このことは、種々の計算法のどれを用いても C-H または C-D 伸縮振動は振動数比がほぼ一定値になるのに対し、塩素原子が関与する振動の振動数比は計算法によってかなりばらつくことを意味する。
すなわち、塩素原子が関与する振動の振動数比は基底関数により大きく変動し、6-31G、6-311G、D95などの基底関数を使った場合は分極関数を含む基底関数を使った場合に比べかなり低い値になる。一方、分極関数を含む基底関数を使った場合には基底関数の影響はごくわずかであり、塩素原子が関与する振動の振動数の計算には分極関数を含む基底関数を使う必要があることが分かる。また、塩素原子が関与する振動の振動数は計算法によっても大きく変動し、RHFやMP2に比べてDFT法は低目の振動数を与える傾向が認められるが、この傾向は文献 [51, 52] においても指摘されている。
赤外強度についても、振動数の場合と同様、各振動ごとに165種類の計算法について求めた強度比の平均値と標準偏差を求めた結果を Table 5 に示す。赤外強度の計算値は振動数の場合よりも実測値との差が大きい。赤外不活性および観測されない振動を除く72個の振動についての強度比の平均値は 1.6693、標準偏差は2.2624となる。上記のように、調和振動の近似の下で計算される振動数あるいは赤外強度の計算値と、振動の非調和性が含まれている実測値との比の大半は分子振動の非調和性によるものと考えられる。

Table 5. Observed intensities (km/mol), their errors (km/mol), means and standard deviations of intensity ratios

		Obsd Int	Error	Mean	Std Dev			Obsd Int	Error	Mean	Std Dev
CH4	n4	34.00	2.00	1.5545	0.6922	CH2Cl2	n4	0.60	0.10	0.8921	0.3646
	n3	74.00	5.00	1.1845	0.3685		n3	8.00	0.40	1.7307	0.4624
CH3D	n6	15.70	2.00	0.7410	0.3285		n9	95.00	8.00	1.6580	0.2374
	n3	12.20	1.50	1.4853	0.6584		n7	1.20	0.10	1.5378	0.9025
	n5	4.10	1.00	0.7697	0.3773		n8	26.60	1.20	1.9008	0.5058
	n2	7.20	1.00	1.1487	0.4001		n2	0.60	0.10	0.2190	0.3512
	n1	7.20	2.00	1.2204	0.3830		n1	6.90	0.50	1.2068	0.5075
	n4	53.00	2.00	0.5452	0.1730		n6	0.00
CH2D2	n4	5.20	1.50	1.4577	0.6316	CHDCl2	n6	0.60	0.10	0.9143	0.3763
	n7	8.10	2.00	1.4943	0.6740		n5	5.70	0.40	2.0031	0.5601
	n9	11.10	2.00	1.4773	0.6523		n9	72.00	3.00	1.8017	0.3914
	n3	4.40	1.50	1.5018	0.6967		n4	2.80	1.00	0.8065	0.2446
	n2	5.60	1.50	0.8440	0.3020		n8	29.00	0.70	1.3647	0.3941
	n8	12.20	2.00	1.0631	0.3607		n7	17.00	0.40	2.0500	0.5569
	n1	13.80	2.00	1.1684	0.3687		n3	0.60	0.10	1.5514	0.6371
	n6	26.00	2.00	1.1002	0.3494		n2	2.40	0.10	1.7031	0.7487
CHD3	n3	7.00	0.50	1.4838	0.6657		n1	3.10	0.10	1.2846	0.6501
	n6	8.10	0.50	1.8310	0.7745	CD2Cl2	n4	0.60	0.10	0.9332	0.3768
	n5	7.10	0.50	1.7175	0.8150		n3	8.00	1.50	1.5128	0.3973
	n2	2.50	1.00	0.8239	0.3025		n7	0.00
	n4	22.70	1.50	1.1348	0.3858		n9	67.00	3.00	1.7922	0.4022
	n1	18.30	1.50	1.2269	0.3890		n8	50.00	2.00	1.6114	0.3760
CD4	n4	20.00	1.00	1.5942	0.7056		n2	0.20	0.10	3.6932	2.9995
	n3	35.00	1.50	1.0963	0.3736		n1	4.30	0.20	1.5944	0.5997
CH3Cl	n3	22.00	1.00	1.3736	0.4285		n6	0.00
	n6	4.10	0.50	1.8184	0.8486	CHCl3	n6	0.05	0.03	4.3290	6.1091
	n2	7.50	0.50	2.2169	0.8416		n3	0.50	0.10	0.7144	0.5820
	n5	12.00	1.00	1.0708	0.2770		n2	4.40	0.60	1.3646	0.5768
	n1	18.00	1.00	1.4064	0.2831		n5	226.00	1.10	1.5551	0.2475
	n4	9.70	1.00	1.5198	0.7351		n4	31.40	2.60	1.7132	0.4066
CD3Cl	n3	15.30	1.00	1.5084	0.4997		n1	0.32	0.10	2.2949	2.5719
	n6	1.10	0.50	1.6358	1.0940	CDCl3	n6	0.05	0.03	4.2205	5.9592
	n2	10.70	1.00	1.7684	0.6245		n3	0.30	0.20	1.2813	1.0256
	n5	6.90	1.00	0.9723	0.2640		n2	5.20	0.10	1.0554	0.4471
	n1	11.30	0.50	1.5961	0.3439		n5	157.00	7.00	1.7301	0.4563
	n4	5.10	0.50	1.5924	0.7755		n4	87.50	1.70	1.4365	0.3701
							n1	0.06	0.04	6.4807	12.0997
						CCl4	n4	0.20	0.05	6.7870	7.4114
							n3	322.00	60.00	1.6749	0.2921

Observed intensities and their errors taken from ref. [43].

したがって、上記の振動数比1.0205は振動数に対する非調和性の影響が僅か2.05 % にすぎないことを意味する。それに対して、強度比1.6693という結果は赤外強度に対する非調和性の影響が非常に大きいこと、およびその影響が赤外強度を減少させる方向にあることを意味している。振動数に対する非調和性の影響に関する研究は多いが、赤外強度に対する非調和性の影響に関する研究はないので、この結果を検証できる実験データはないが、興味深い結果である。
赤外強度の計算値と実測値の違いのもう1つの原因として考えられるのが電子相関の効果である。本研究の計算結果から計算法ごとに強度比の平均値を求めると、RHFは2.1788、B3LYPは1.4878、BLYPは1.7240、BPW91は1.5467、MP2は1.4600となり、電子相関が取り入れられていないRHFでは計算値が実測値の2倍以上の値になる。赤外強度をMP3以上の電子相関補正法で計算した論文は見あたらないので、DFT法やMP2における強度比がほぼ1.5という値が振動の非調和性によるものか、電子相関によるものかどうかは不明であるが、この点も興味ある結果である。
赤外強度の場合は振動数の場合と異なり、実測強度の測定誤差が大きな振動が多く、また強度比の標準偏差が大きい振動も多いので、強度に対する振動形の影響は認めにくい。しかし、注意すべきはCHCl₃のn₆、CDCl₃のn₆とn₁、CCl₄のn₄における標準偏差の異常な大きさである。このことはこれらの赤外強度が計算法によって大きく変動することを意味する。事実、標準偏差が最大のCDCl₃のn₁では強度比はMP2/D95での72.9からMP2/cc-pVTZでの0.000まで大きな差がある。この振動はC-D伸縮振動であり、塩素原子が関与する振動ではない。しかし、この分子では塩素原子が3個存在するため、C-D結合の伸縮振動における電子密度の変化量 (赤外強度に比例する) が基底関数、特に電子相関の効果によって大きく影響されるためであると考えられる。
分子軌道法による理論計算が赤外スペクトルをどの程度まで再現するかという観点からは、振動数や赤外強度の計算値の絶対値が実測値と大きくかけ離れていても問題ではなく、重要なことは計算値と実測値の比率が一定であることである。すなわち、その比率が一定であれば比率をscale factorとして計算値を補正すれば実測値が再現できることになる。
この観点から各計算法の精度を比較するために、基底関数と計算法の組み合わせについて、全振動 (赤外強度については不活性および観測できない振動は除く) についての振動数比と強度比の平均値と標準誤差を計算した。ただし、実測値には測定誤差が存在するので、この誤差に見合った重みを考慮して平均値と標準偏差を計算した。すなわち、振動数の場合は文献[46] に記載されている振動数の不確かさのクラスによって測定誤差 Dn を A: 1.0、B: 2.0、C: 4.5、D: 10.5 cm^-1 とし、これに基づく重み w_i を
w_i = (1/Dn)²
により計算した。
一方、赤外強度の場合は強度の実測値 I とその測定誤差 DI に対して重み w_i は
w_i = (I/DI) ²
により計算した。これらの重み w_i を用い、振動数比または強度比 x_i の平均 M と標準偏差 SD を
M = Sw_i x_i /Sw_i
SD = [Sw_i (M - x_i)² /Sw_i]^1/2
により計算した。この標準偏差を計算法の精度の目安となる計算誤差とみなした。種々の計算法の内で最適法を決めるために、この標準偏差を全振動について平均した。また、最適の計算法の選定に当たっては計算時間も評価基準の1つになるので、全分子の計算時間 (CPU時間) を平均した。このようにして求めた各計算法の標準偏差と計算時間を Table 6 に示す。この結果から計算法の影響について幾つかの傾向を見いだすことができる。

Table 6. Standard deviations of frequency (SDF) and intensity (SDI) ratios and CPU times (T, sec)

	RHF			B3LYP			BLYP
	SDF	SDI	T	SDF	SDI	T	SDF	SDI	T
6-31G	0.0543	4.9585	162.0	0.0620	0.8870	417.5	0.0758	0.8995	398.0
6-31G(d)	0.0215	2.3211	235.8	0.0261	0.9114	694.9	0.0426	1.8811	645.0
6-31G(d,p)	0.0201	2.2419	212.7	0.0246	0.8272	755.8	0.0416	1.7879	732.7
6-31G(3df,3pd)	0.0181	0.9604	2455.8	0.0214	0.4714	5241.3	0.0373	1.3340	4218.6
6-31+G	0.0541	5.2328	182.3	0.0622	1.0530	588.0	0.0763	0.9934	545.3
6-31+G(d)	0.0216	2.0680	275.3	0.0266	0.6031	962.4	0.0432	1.9392	894.4
6-31+G(d,p)	0.0197	1.9505	257.6	0.0251	0.5638	1037.2	0.0420	1.9195	985.6
6-31+G(3df,3pd)	0.0185	0.8047	3215.5	0.0218	0.7091	6508.6	0.0380	1.7851	5398.4
6-31++G	0.0537	4.9785	198.9	0.0622	0.9183	627.1	0.0767	0.9479	577.9
6-31++G(d)	0.0214	2.0093	280.5	0.0266	0.6094	1008.7	0.0432	1.9663	931.2
6-31++G(d,p)	0.0196	1.8694	262.9	0.0251	0.5754	1083.5	0.0421	2.0024	1024.4
6-31++G(3df,3pd)	0.0181	0.8122	3320.9	0.0217	0.6969	6762.8	0.0381	1.8894	5578.7
6-311G	0.0511	2.9923	226.7	0.0647	1.0823	649.6	0.0797	2.4455	591.9
6-311G(d)	0.0231	2.6583	324.7	0.0278	1.7864	1106.9	0.0441	4.4690	987.0
6-311G(d,p)	0.0207	1.8726	318.3	0.0259	1.5580	1152.7	0.0426	4.4133	1023.6
6-311G(3df,3pd)	0.0198	0.9954	3945.0	0.0225	1.3899	7014.0	0.0397	2.5214	5689.2
6-311+G	0.0480	4.2303	269.5	0.0613	0.7404	899.2	0.0772	1.2433	832.8
6-311+G(d)	0.0228	1.5954	411.2	0.0266	0.8633	1495.3	0.0428	2.8614	1319.8
6-311+G(d,p)	0.0204	1.4028	1058.6	0.0246	0.8431	1542.1	0.0416	3.0586	1369.6
6-311+G(3df,3pd)	0.0195	0.7765	5390.5	0.0224	1.0011	8661.6	0.0395	2.4934	6853.3
6-311++G	0.0476	4.3371	275.6	0.0611	0.6981	952.0	0.0760	1.1935	884.0
6-311++G(d)	0.0225	1.5250	418.5	0.0264	0.8512	1553.0	0.0423	2.8163	1374.4
6-311++G(d,p)	0.0203	1.3453	410.5	0.0245	0.8370	1614.2	0.0409	2.8303	1439.4
6-311++G(3df,3pd)	0.0195	0.7879	8851.9	0.0223	1.0253	8925.2	0.0389	2.1609	7084.0
D95	0.0556	7.1656	193.6	0.0590	5.5569	485.0	0.0680	6.2332	453.3
D95(d)	0.0216	8.3824	402.1	0.0225	4.2274	804.2	0.0371	4.4182	744.3
D95(3df,3pd)	0.0184	0.9399	3022.4	0.0211	0.4982	5870.1	0.0375	1.1316	4803.9
cc-pVDZ	0.0161	1.3587	468.6	0.0251	0.7192	1205.1	0.0381	1.8259	1148.9
cc-pVTZ	0.0187	1.0471	3016.8	0.0262	0.5896	6480.7	0.0423	1.9499	5688.1
cc-pVQZ	0.0188	0.8817	37667.9	0.0244	0.5736	40888.3	0.0409	1.8497	34250.1
AUG-cc-pVDZ	0.0159	0.9569	815.4	0.0266	0.8379	2568.7	0.0444	2.0494	2459.5
AUG-cc-pVTZ	0.0185	0.8538	15834.2	0.0270	0.8558	21269.5	0.0420	2.1714	17523.1
AUG-cc-pVQZ	0.0186	0.8231	187538.6	0.0243	0.7902	207582.3	0.0398	2.0411	141621.1

Table 6B. Standard deviations of frequency (SDF) and intensity (SDI) ratios and CPU times (T, sec)

	BPW91			MP2
	SDF	SDI	T	SDF	SDI	T
6-31G	0.0618	0.7662	422.2	0.0631	2.0596	401.6
6-31G(d)	0.0279	1.1719	674.5	0.0119	1.2526	674.7
6-31G(d,p)	0.0274	1.0858	760.2	0.0124	1.8623	1062.7
6-31G(3df,3pd)	0.0241	0.7282	4531.1	0.0140	0.4840	23526.0
6-31+G	0.0622	0.8435	582.5	0.0627	2.0183	671.7
6-31+G(d)	0.0288	1.0948	942.1	0.0121	0.9868	1401.1
6-31+G(d,p)	0.0282	1.0756	1025.7	0.0133	0.9977	1752.7
6-31+G(3df,3pd)	0.0245	0.9432	5554.1	0.0146	0.7267	35256.2
6-31++G	0.0621	0.7912	616.6	0.0627	1.8551	693.6
6-31++G(d)	0.0287	1.0713	984.8	0.0128	0.9131	1675.7
6-31++G(d,p)	0.0282	1.0599	1079.2	0.0144	0.9577	1811.5
6-31++G(3df,3pd)	0.0244	0.9405	5751.3	0.0143	0.7240	36377.7
6-311G	0.0640	2.0806	619.0	0.0584	1.5001	998.9
6-311G(d)	0.0288	3.1829	1020.4	0.0134	1.3524	2143.1
6-311G(d,p)	0.0282	2.9490	1090.0	0.0162	2.0261	2335.0
6-311G(3df,3pd)	0.0251	1.4643	5796.5	0.0120	1.4083	39071.5
6-311+G	0.0611	0.9072	866.7	0.0543	1.3863	1842.6
6-311+G(d)	0.0276	1.6028	1385.4	0.0126	0.5928	3467.3
6-311+G(d,p)	0.0271	1.5587	1436.4	0.0160	0.7798	20654.4
6-311+G(3df,3pd)	0.0250	1.1392	6919.5	0.0122	0.9132	38668.9
6-311++G	0.0608	0.8459	921.0	0.0538	1.3970	2480.4
6-311++G(d)	0.0275	1.5724	1439.6	0.0126	0.5751	3587.5
6-311++G(d,p)	0.0270	1.5485	1512.3	0.0162	0.7331	20804.4
6-311++G(3df,3pd)	0.0250	1.1446	7191.3	0.0121	0.9301	39308.2
D95	0.0603	5.7000	479.4	0.0607	8.6427	1023.4
D95(d)	0.0252	4.2460	781.0	0.0144	5.2768	15689.8
D95(3df,3pd)	0.0244	0.5802	4904.1	0.0150	0.6629	22594.7
cc-pVDZ	0.0289	1.1041	1187.8	0.0183	1.6280	24458.2
cc-pVTZ	0.0285	1.0033	5796.2	0.0145	0.5230	31159.1
cc-pVQZ	0.0273	0.9122	33913.6	0.0147	0.5230	226014.8
AUG-cc-pVDZ	0.0302	1.0043	2557.5	0.0195	0.5668	7349.4
AUG-cc-pVTZ	0.0292	1.1123	17806.9	0.0114	0.7771	76835.4
AUG-cc-pVQZ	0.0270	0.9855	138130.3	0.0143	0.3287	549134.6

まず第1は振動数の計算誤差に対する基底関数の影響である。たとえば、計算法RHFにおいて基底関数を6-31G、6-31G(d)、6-31G(d,p)、6-31G(3df,3pd) と変えると計算誤差は0.0543、0.0215、0.0201、0.0181 と大幅に減少する。これに対して、同じRHFにおいて基底関数を6-31G、6-31+G、6-31++G、6-311G、6-311+G、6-311++G と変えても計算誤差は0.0543、0.0541、0.0537、0.0511、0.0480、0.0476 と変化は小さい。これは上記のように、本研究の塩素原子を含む分子では分極関数の影響が大きいが、diffuse 関数の影響は小さいことを示している。
また、電子相関の効果としてRHFをMP2に変えると計算誤差は全般的にはほぼ半減しているが、分極関数を含まない6-31G、6-31+G、6-31++G、6-311G、6-311+G、6-311++G ではRHFとMP2とで計算誤差が減少せず、むしろMP2の方が若干ではあるが誤差が大きい。これは電子相関の補正の効果を引き出すためには分極関数を取り入れた基底関数を用いる必要があることを示しているが、MP2で誤差が増えた原因は不明である。
第2点は振動数の計算誤差に対する計算法の影響である。Table 6 の結果から各計算法における計算誤差の平均を算出すると、RHFでは0.0266、B3LYPでは0.0325、BLYPでは0.0482、BPW91では0.0344、MP2では0.0237となる。この結果から、まず本研究の場合には RHF が DFT 法より良い精度を与えることが分かる。多数の分子についての振動数の計算結果から、RHF が最も精度が悪く、MP2 がそれに次ぎ、DFT が最も良い精度であるとの報告がある [49]。彼らと我々では分子の数と種類が違うので結果の違いは当然であるが、上記のようにDFT法は塩素原子の関与する振動の振動数を低く与える傾向があり、これが誤差を大きくしている。
また、MP2 は RHF より計算時間がほぼ2倍程度かかるのに計算精度はそれほど向上しないことが分かる。さらに、DFT 法の中ではこれまで B3LYP が比較的よい結果を与えるとされているが、本研究の結果では BPW91 も B3LYP に匹敵する誤差になった。この結果は文献 [52] に示されている B3LYP とB3PW91 が最も誤差が小さいという結果と対応している。
第3点は計算精度と計算時間の関係であり、全般的に基底関数や計算法の精度を上げれば、計算時間は長くなるが振動数の計算誤差は減少するという傾向が認められる。しかし、大幅な計算時間がかかるcc-pVQZ、AUG-cc-pVTZ、AUG-cc-pVQZ などの基底関数を用いた場合に必ずしも最高の計算精度が得られるとは限らない点である。Table 6 の結果において、各計算法の中で最高の計算精度が得られた基底関数はRHFではAUG-cc-pVDZ、B3LYPではD95(3df,3pd)、BLYPではD95(d)、BPW91では6-31G(3df,3pd)、MP2ではAUG-cc-pVTZ であり、これらの基底関数の場合の計算時間は各計算法の中で最長ではない。したがって、計算精度が最高で、かつ計算時間も適切な最適な計算法が存在することが分かる。
第4点は赤外強度に関する全般的傾向として、上記の振動数に関する傾向ほど明確ではないが、ほぼ同様の傾向が認められる。この場合に特筆すべき点として、計算法 B3LYP、BPW91 では基底関数 6-31G(3pd,3df) との組み合わせにおいて赤外強度の計算誤差が最小となっている。振動数の再現においては DFT法、特にB3LYP がよい結果を与えることはこれまでの研究から判明しているが、赤外強度の再現に関しては種々の計算法を比較した研究はこれまで皆無である。本研究では塩素化メタン分子しか扱っていないので、一般の分子についてこの結論があてはまるかどうかは不明であるが、比較的計算時間の短い基底関数 6-31G(3pd,3df) と B3LYP、BPW91の組み合わせが赤外強度の計算法として最適であることは興味深く、他の多くの分子についての検討が待たれる。
Table 6の結果から最適の計算法を決定することができるが、計算の目的によって最適の計算法を選定する基準が異なる。赤外スペクトルを理論計算する目的はいろいろあるが、実測の赤外スペクトルに現れているピークがどの振動モードによるかを決める振動の帰属を行う場合には、振動数の計算値が実測値にできるだけ近く、計算値の絶対的な精度が高い計算法であることが要求される。このような場合は一般に少数の分子について計算を行うことが多いので計算時間はあまり問題にならない。
そこで、この観点から Table 6 の結果の中から標準偏差の最も小さい計算法を探すと、振動数に対しては AUG-cc-pVTZ の基底関数を用いて MP2 計算を行うと標準偏差 0.0114、すなわち相対誤差1.14 % で実測振動数が再現できることになる。この相対誤差 1.14 % という数値は C-H 伸縮振動が現れる 3000 cm^-1 領域において平均 34 cm^-1 の計算誤差を与えることを意味するので、実測の赤外ピークの帰属には絶対的な精度として満足できる計算法とはいえない。
また、赤外強度に対しては AUG-cc-pVQZ の基底関数を用いて MP2 計算を行うと標準偏差 0.3287、すなわち相対誤差 32.87 % で実測強度が再現できることになる。赤外強度の測定誤差は実測強度の通常10 % 程度なので、強度の計算誤差 32.87 % も絶対的な誤差としては必ずしも満足できる値ではない。したがって、このような実測ピークの帰属を行う場合には、これらの計算法を用いて振動数と赤外強度を計算し、本研究で行ったように振動数だけでなく赤外強度の計算値も考慮に入れて慎重に帰属を行うことが必要であると結論できる。
以上の場合は計算値の絶対的な精度の高さが要求されるので、現状の計算法でこのような問題を解決することは難題である。しかし、分子の環境の変化などに由来する実測スペクトルの変化を解析する場合には、振動数や赤外強度の相対的な変化が問題となるので、上記ほどの精度の高さは必要ない。たとえば、孤立分子と分子錯体について振動数を分子軌道法で計算し、マトリックス単離法で観測した赤外スペクトルとの比較から、分子錯体の存在およびその構造を推定した研究がある [53 - 55]。このような問題の場合には振動数の絶対値は問題でなく、相対的なシフトが計算で説明できればよいので、上記の計算法で十分解決可能と考えられる。
赤外スペクトルを理論計算する3番目の目的として、緒言で記した未知化学種の同定がある。この場合には無数の化学種の中から実測の赤外スペクトルに一致する計算スペクトルを示す化学種を探索することが必要となる。スペクトルの一致は振動数や赤外強度の数値ではなくスペクトルの図形として判断するので、計算精度はあまり問題にならない。本研究で検討した計算法の中で最も近似の低い計算法は基底関数 6-31G と計算法 RHF の組み合わせであるが、この場合の振動数の計算誤差 5.43 % はスペクトル図形の再現には十分な精度である。
一方、赤外強度の計算誤差 495.85 % は測定誤差と比較すると満足できる精度ではない。しかし、未知化学種の赤外スペクトルを測定する場合、スペクトル中の最大ピークが通常 5 % 程度の透過率になるようセルの厚さを調節して測定が行われるので、その絶対的な強度は問題にならない。したがって、このような目的では相対的な強度が再現できれば十分であり、この点で上記の計算誤差はこの目的には十分な計算精度である。むしろ、この場合は莫大な数の化学種について赤外スペクトルを計算する必要があるので、最大の選定基準は計算時間である。上記の計算法での平均計算時間は162.0秒であり、多数の化学種の赤外スペクトルを計算できる方法とはいえない。この目的には同程度の計算精度でかつ計算時間が圧倒的に短い計算法の開発が望まれる。

4. 3 理論的方法の問題点

5 考察―各種の方法の比較

6 結 言

一連の本研究を行うにあたり、SDBSの構築やニューラルネットワークシステムの評価などに貢献していただいた沢山の方々に感謝する。

参考文献