Regiospecific Addition of La@C82 with Alkylidene, Alkyl Radical and Carbanion: Prediction of Regioisomers by Using "Paired Interacting Orbitals" (PIO) Analysis

Akinobu SHIGA, Midori O. ISHITSUKA and Takeshi AKASAKA

Return

1 Introduction

Endohedral metallofullerenes (EMFs) are especially expected to have a wide range of potential applications in nanomaterials science. Chemical modification of EMFs is very important for their applications. The reactivities of EMFs are quite different to those of empty fullerenes. Several monoadditions to EMFs take place with remarkable regioselectivity and afford few or even only one regioisomer. Nagase et al, reported the structures and molecular properties, charges and p-orbital axis vector (POAV) values of EMFs by using DFT calculation [1]. Prediction of the regioselectivity of EMFs is a very complicated problem, therefore, it is necessary for the prediction to apply reactivity indexes from many different points of view. Orbital interaction is one of the powerful candidates. Fujimoto et al. [2] proposed the paired interacting orbitals (PIO) calculation that is a method for unequivocally determining the orbitals which should play dominant roles in chemical interactions between two systems. Applying this method to the (PH3)2Pd/C60 system and the C60/Cu(111) surface system, they reported that the interaction of the former is mainly expressed by only two PIOs and that of the latter is mainly expressed by three PIOs [3, 4]. Since the molecular symmetry of La@C82 is C2v, La@C82 possesses 24 non-equivalent carbon atoms. Regioselective reaction takes place according to the chemical natures of reagents and the reactivities of these non equivalent carbons. Here, we applied the PIO analysis for addition of adamantylidene(Ad), methyl anion(CH3-), diethyl bromomalonate anion(DBEM-) and methyl radical(.CH3) to La@C82 to get a reactivity index of the regioselective addition.

2 Models and Calculation Methods

2. 1 Models

The models employed here are La@C82 and four reagents: Ad, CH3-, DBEM- and .CH3. We present the structure of La@C82 in supporting information of reference[1]. We optimized the structures of these reagents by RHF method (SBKJC ECP basis sets) using PC GAMESS developed by Dr. Alex A. Granovsky of the Laboratory of Chemical Cybernetics at Moscow State University [5].
We analyzed an early stage of the addition of each reagent to all the twenty four non-equivalent carbons of the La@C82.
We determined plausible reaction paths of methylidene and CH3- addition on ethylene and then, selected a model state of the early stage of the addition by considering the change of the PIO analysis of each model state on the reaction path.
We generated a model of the early stage of the addition of each reagent to the La@C82 by placing the active site carbon atom of each reagent at the point, the distance of which is 2.500A above the carbon atom of La@C82. Twenty four non-equivalent carbon atoms of La@C82, having C2v symmetry, are shown in Figure 1. Examples of the 3D structures of the model of the early stage of each addition, [C], are shown in Figure 2.

2. 2 PIO calculations

We applied paired interacting orbitals (PIO) calculation, proposed by Fujimoto et al. [2], to the combined system [C]. Here, [C] is the model of the early stage of the addition, while [A] is the La@C82 portion and [B] is the reagent portion. The geometries of [A] and [B] are the same as those in the complex [C] ([A-B]@[C]). The extended Huckel MOs of [A], [B] and [C] are calculated. The extended Huckel parameters are given in Appendix. PIOs are obtained by applying the procedure that was proposed by Fujimoto et al. [2]. It is summarized as follows:
1) we expand the MOs of a complex in terms of the MOs of two fragment species, to determine the expansion coefficients ci, f, cm+j, f and dk, f, dn+l, f in Eq.1

2) we construct an interaction matrix P which represents the interaction between the MOs of the fragment [A] and the MOs of the fragment [B]

in which

here q is an electron number,
3) we obtain transformation matrix UA (for A) and UB (for B) by

here g is an eigen value,
4) and finally we obtain the PIOs by Eq. 5 and Eq. 6

The N×M (N<M) orbital interactions in the complex C can thus be reduced to the interactions of N PIOs, N indicating the smaller of the numbers of MOs of the two fragments, A and B.
The eigen values mean the contribution of each PIO to the interaction. The larger the eigen value, the larger is the contribution of the PIO to the interaction.
Although PIO analysis was proposed originally for ab initio calculations, we already reported that this approach was also useful in analyzing the results of extended Huckel MO calculations and was suitable for understanding reaction procedures in terms of orbital interactions, especially in catalytic systems large in size [6]. We also reported that relative catalytic activities deduced from PIO analysis based on extended Huckel MO calculations coincided with the order of activation energy obtained from ab initio MO calculations [7]. Extended Huckel parameters are given in Appendix. All calculations were carried out on LUMMOXTM system [8].

3 Elucidation of a plausible reaction path of methylidene addition on ethylene

3. 1 Search of a methylidene approaching route to ethylene

As shown in Figure 3 and Table 1, we put a methylidene carbon atom above the molecular plane of ethylene and searched several points above the ethylene by changing the distance r between mC and Ca of ethylene and the angle q : (mC,Ca,Cb).


Figure 1. Twenty four non equivalent carbons in La@C82


Figure 2. The 3D structures of the model of the early stage of each addition, [C], a): Ad addition on carbon(7), b): DBEM- addition on carbon(23), c): .CH3 addition on carbon(10)


Figure 3. Schematic illustration of methylidene approaching position to the ethylene

Table 1. The methylidene approaching position
point no(p1)(p2)(p3)(p4)(p5)(p6)(p7)(p8)(CH2)3
r(A)3.0002.7002.6002.5002.5002.3002.2002.1001.513
q(°)909090907070707060

The DE(kcal/mol) energies of each point model are shown Table 2. Figure 4 is the energy diagram of them.

Table 2. Energy of each point No.(pn)
DE (kcal/mol) = E(pn) - Estarting CH2 and C2H4
start(p1)triplet(p2)triplet(p3)triplet(p4)triplet(p5)triplet(p5)singlet(p6)triplet(p6)singlet(p7)singlet(p8)singlet(CH2)3
0.0+4.59+7.40+9.30+12.9+19.6+24.9+29.3+22.7+18.4+14.3-5.99


Figure 4. Energy diagrams of each model

According to the approaching of the methylidene to the ethylene, the DE of triplet state is increasing. The DE of triplet state and that of singlet state is crossing over between the point (p5) and the point (p6). After that, the DE of singlet state is decreasing toward the energy of final product, cyclopropane.

3. 2 PIO analysis of each point model on the methylidene and Ad approaching route to the ethylene

As already described in 2.2, we divide each point model[C] into two fragments, [A] and [B], here, [A] is the ethylene and [B] is the methylidene, and PIO calculation is executed.
Eigen values
Eigen values of each model are shown in Table 3. Eigen values of PIOs express the contribution of the PIOs to the interaction between [A] and [B]. The larger the eigen values, the larger is the contribution of the PIO to the interaction.

Table 3. The eigen values of PIOs of each point model
ModelPIO-1PIO-2PIO-3PIO-4PIO-5...
(p1)0.04280.00310.00060.00030.0000...
(p2)0.17990.05930.00520.00100.0000...
(p3)0.58380.06060.00940.00150.0000...
(p4)0.63380.05950.01240.00230.0000...
(p5)0.66260.16030.02290.00340.0000...
(p6)0.60150.16660.04500.01390.0003...
(p7)0.60470.23390.05370.01630.0004...
(p8)0.60020.30680.06150.02000.0004...
Ad addition
(p4)0.32400.04200.00840.00220.0007...
(p5)0.61100.30470.02990.00390.0011...

Table 3 tells us that the main interaction between methylidene and ethylene is expressed in the PIO-1. Contour maps of PIO-1 of each point are shown in Figure 5. It is clearly observed out-of-phase interaction between methylidene and Cb of ethylene in (p1) and in (p2), however, PIO-1 turns in-phase interaction between methylidene LUMO and p-orbital of ethylene in (p3) and after that the in-phase interactions are gradually increased in (p4) to (p8).
Although contour maps of PIO-1 of (p4) of Ad still shows out-of-phase interaction between Ad and Cb of ethylene, that of PIO-1 turns in-phase one in the case of (p5) of Ad.


Figure 5. Contour maps of PIO-1 of model of each point

4 Elucidation of a plausible reaction path of carbanion addition on ethylene

4. 1 Search of carbanion approaching route to ethylene

As described in 3.1, we put a carbon atom of CH3- above the molecular plane of ethylene and searched several points in Table 4.

Table 4. The CH3- approaching position
point no(p9)(p10)(p11)(p12)(p13)(p14)(p15)
r(A)3.0003.0002.7002.7002.5002.5001.900
q(°)901099010990109109

The DE(kcal/mol) energies of each point model, shown in Table 5, tells us that as the Ca in final product is sp3 carbon, it is easier to reach the final product via the approaching route of CH3- with q:109°,(p10), (p12), (p14), than via the route with q:90°,(p9), (p11), (p13). It is suggested that the TS of CH3- addition locates near the point (p14).

Table 5. Energy of each point No.(pn)
DE (kcal/mol) = E(pn) - Estarting CH2 and C2H4
start(p9)(p10)(p11)(p12)(p13)(p14)(p15)
0.0+10.0+1.59+13.8+0.84+16.6+0.48-11.1

4. 2 PIO analysis of several point models on CH3- and DBEM- approaching route to the ethylene

We selected four point models, (p11), (p12), (p13) and (p14). PIO calculation is executed as described in 2.2, here, [A] is the ethylene and [B] is the carbanions.
Eigen values
Eigen values, shown in Table 6, tell us that the main interaction between carbanion and ethylene is mainly expressed in the PIO-1.

Table 6. The eigen values of PIOs of each point model of carbanion approaching route
ModelPIO-1PIO-2PIO-3PIO-4PIO-5...
(p11)0.16640.01370.00130.00110.0002...
(p12)0.21000.01190.00110.00100.0001...
(p13)0.32600.02860.00290.00250.0003...
(p14)0.36500.02310.00260.00230.0002...
DBEM- addition
(p14)0.18370.03070.01000.00170.0002...

Contour maps of PIO-1 of each point are shown in Figure 6.


Figure 6. Contour maps of PIO-1 of model of each point

We can observe in-phase interaction between carbanion HOMO and p*-orbital of ethylene in each model.
Approaching near to the TS, the reagents incline toward C-C double bond to form cycopropanes in the case of alkylidene addition, on the other, the reagents incline in the opposite direction of C-C bond to make sp3 Ca in the case of carbanion addition. Since every Ca of La@C82 possesses three double bonds, a perpendicular position up to Ca is neutral for all three bond directions. From above results and consideration, we determined a model of the early stage of the addition of each reagent to the La@C82 by placing the active site carbon atom of each reagent at the point, the distance of which is 2.500A above the carbon atom of La@C82.

5 PIO analysis of each addition reaction

5. 1 PIO analysis of the intermediate models of the addition of Ad to La@C82

PIO calculation is executed as described in 2.2, here, [A] is the La@C82 and [B] is the Ad.
Eigen values
Eigen values of each model are summarized in Table 7.

Table 7. The eigen values of PIOs of each intermediate model of Ad addition
PIO-1PIO-2PIO-3...PIO-54
Modela spinb spina spinb spina spinb spin...a spinb spin
(2)0.0890.1480.0430.0740.0050.009...0.0000.000
(3)0.1660.1810.0440.0310.0050.008...0.0000.000
(4)0.1320.1100.0610.0360.0070.008...0.0000.000
(5)0.0830.0840.0350.0500.0050.007...0.0000.000
(6)0.1520.1100.0330.0570.0050.007...0.0000.000
(7)0.1280.1160.0280.0680.0040.008...0.0000.000
(8)0.1040.1460.0620.0260.0060.010...0.0000.000
(9)0.1790.1350.0280.0810.0050.008...0.0000.000
(10)0.1810.1300.0360.0770.0060.008...0.0000.000
(11)0.0770.0800.0330.0400.0050.006...0.0000.000
(12)0.0860.0890.0420.0530.0060.007...0.0000.000
(13)0.1500.0980.0420.0380.0060.005...0.0000.000
(14)0.2020.1820.0460.0270.0050.009...0.0000.000
(15)0.1500.0840.0350.0300.0040.004...0.0000.000
(16)0.0770.0770.0260.0250.0050.005...0.0000.000
(17)0.1030.1340.0540.0570.0070.008...0.0000.000
(18)0.1620.1100.0620.0600.0080.008...0.0000.000
(19)0.0740.0740.0350.0350.0050.005...0.0000.000
(20)0.0870.0920.0410.0520.0050.006...0.0000.000
(21)0.1380.1070.0360.0540.0060.007...0.0000.000
(22)0.1790.1150.0480.0610.0070.007...0.0000.000
(23)0.1630.1140.0660.0590.0080.007...0.0000.000
(24)0.0810.0810.0250.0250.0050.005...0.0000.000
(25)0.0840.0840.0300.0320.0040.004...0.0000.000

There are 165 occupied MOs in which more than forty p type MOs are included, SOMO, and 171 unoccupied MOs in which more than twenty p* type MOs are included, in the extended Huckel MOs of La@C82 and there are 27 occupied MOs and 27 unoccupied MOs in the same MOs of Ad. The 337 × 54 orbital interactions are compactly represented in 54 PIOs. Moreover, the eigen values tell us that the interaction between the La@C82 and the Ad is mainly expressed in only one PIO, PIO-1.
The typical PIO-1, for example, model(3), (6), (7), (9), (10), (13), (14), (15), (18), (21), (22), and (23) are composed from the SOMO and occupied orbitals near SOMO of the La@C82 and the LUMO of the Ad.
Contour maps of PIO-1
Contour maps of PIO-1 of model(6), (7), (9) and (14) are shown in Figure 7. It can be seen in each PIO-1 that occupied orbitals of the La@C82 are localized on the carbon atom of the reaction center and interacted on the LUMO of the Ad.


Figure 7. Contour maps of PIO-1 (a-spin part) of a): model(6), b): model(7), c): model(9), and d): model(14)

5. 2 PIO analysis of the models of the addition of CH3- to La@C82

Eigen values
Eigen values of each model of the CH3- addition are summarized in Table 8.
Here, the 337 × 7 orbital interactions are compactly represented in 7 PIOs and the eigen values tell us that the interaction between the La@C82 and the CH3- is mainly expressed in only one PIO, PIO-1. The typical PIO-1, for example, model(3), (6), (7), (9), (10), (13), (14), (15), (18), (21), (22), and (23) are composed from the SOMO and, or unoccupied orbitals near SOMO of the La@C82 and the HOMO of the CH3-.

Table 8. The eigen values of PIOs of each model of the early stage of CH3- addition
PIO-1PIO-2PIO-3...PIO-7
Modela spinb spina spinb spina spinb spin...a spinb spin
(2)0.0660.0720.0050.0060.0010.001...0.0000.000
(3)0.0630.1400.0050.0090.0020.002...0.0000.000
(4)0.0740.1000.0060.0080.0010.001...0.0000.000
(5)0.0710.0860.0060.0070.0020.002...0.0000.000
(6)0.0670.1360.0050.0080.0010.001...0.0000.000
(7)0.0750.1450.0060.0080.0020.002...0.0000.000
(8)0.0720.1150.0060.0080.0020.002...0.0000.000
(9)0.0690.1570.0060.0080.0010.001...0.0000.000
(10)0.0710.1530.0060.0080.0010.001...0.0000.000
(11)0.0690.0710.0060.0060.0010.001...0.0000.000
(12)0.0800.0960.0070.0070.0010.001...0.0000.000
(13)0.0880.1180.0070.0080.0010.001...0.0000.000
(14)0.0730.1500.0060.0080.0010.009...0.0000.000
(15)0.0790.1090.0070.0080.0010.001...0.0000.000
(16)0.0750.0750.0060.0060.0010.001...0.0000.000
(17)0.0920.1040.0070.0080.0010.001...0.0000.000
(18)0.0840.1310.0070.0080.0010.001...0.0000.000
(19)0.0720.0720.0060.0060.0010.001...0.0000.000
(20)0.0990.1020.0070.0080.0010.001...0.0000.000
(21)0.0730.1240.0060.0080.0010.001...0.0000.000
(22)0.0750.1400.0060.0080.0010.001...0.0000.000
(23)0.1050.1380.0070.0080.0010.001...0.0000.000
(24)0.0780.0780.0060.0060.0010.001...0.0000.000
(25)0.0880.0890.0070.0070.0010.001...0.0000.000

Contour maps of PIO-1
Contour maps of PIO-1 of model(10), (14) and (23)are shown in Figure 8.
It can be seen in each PIO-1 that unoccupied orbitals of the La@C82 are localized on the carbon atom of the reaction center and interacted on the HOMO of the CH3-.


Figure 8. Contour maps of PIO-1 (a-spin part) of a): model(10), b): model(14) and c): model(23)

5. 3 Overlap population(OP) of PIO-1 of Ad, CH3-, DBEM- and .CH3 addition

The OP of the PIO-1 of each addition are summarized in Table 9 together with the values of charge and spin density. The larger the OP, the larger is the reactivity of the reaction.

Table 9. Overlap Population of PIO-1 of each addition
Model Xa):AdCH3-DBEM.CH3chargeb)spin densityb)
(2)_X-0.0129-0.0554-0.045
(3)_X0.0298-0.0203-0.099
(4)_X0.0007-0.0292-0.033
(5)_X-0.0090-0.0344-0.092
(6)_X0.0429-0.0138-0.071
(7)_X0.0439-0.00580.0344-0.1360.002
(8)_X0.0166-0.01180.0120-0.1700.016
(9)_X0.0433-0.00300.0499-0.0220.026
(10)_X0.0414+0.00050.0493-0.0200.032
(11)_X-0.0113-0.0465-0.061
(12)_X-0.0110-0.0235-0.047
(13)_X0.0415-0.0073-0.021
(14)_X0.0626+0.00030.04460.0000.023
(15)_X0.0324-0.0217-0.012
(16)_X-0.0326-0.0308-0.036
(17)_X-0.0109-0.0115-0.006
(18)_X0.0262-0.00460.0210-0.0060.024
(19)_X-0.0114-0.0348-0.026
(20)_X-0.0111-0.00360.002
(21)_X0.0234-0.0169-0.027
(22)_X0.0336-0.0024-0.02970.03350.0040.030
(23)_X0.0292+0.0323-0.02200.03900.0060.025
(24)_X-0.0278-0.0278-0.025
(25)_X-0.0010-0.0231-0.0226-0.0060.000
a) Model(m) means carbon number in Figure 1 and X means Ad, CH3-, DBEM, or CH3..b) Ref. [1]

We can estimate a regioselectivity of each addition from the OP; 1) order of the Ad addition: (14) > (7) > (9) > (6) > (13) > (10) ... , 2) order of the CH3-, DBEM- addition: (23) > (10) > (14) > (22) > (9) > ... , 3) order of the .CH3 addition: (9) (10) > (14) > (23) ... .
Ad addition will be expected to take place on the carbon atom having minus charge. The charges of the above carbon are (14):0.000, (7):-0.136, (9):-0.022, (6):-0.071, (13):-0.021, and (10):-0.021. The experimental result of the Ad addition takes place on the carbon(7). The carbon having the largest minus charge is (8): -0.170, however its OP is very small, 0.0166. A straight forward relationship between the OP and the charge is not seen. This is discussed in the next section.
CH3- and DBEM- addition will be expected to take place on the carbon atom having plus charge. The charges of the carbon (23), (10), (14), (22) and (9) are 0.006, -0.020, 0.000, 0.004 and -0.022, respectively,. The experimental result of the DBEM- addition takes place on the carbon(23). Prediction from OP and that from charge coincides with each other, however the charges of (10) and (14) are more negative than that of (22). The relationship between the OP and the charge is also not clear. This is also discussed in the next section.
.CH3 addition will be expected to take place on the carbon atom having large spin density. The spin densities of the carbon (9), (10), (14) and (23) are 0.026, 0.032, 0.023 and 0.025, respectively,. A fairly parallel relationship between the OP and the spin density is observed.

6 Discussion

In the case of Ad addition, the OP of model(14): 0.0626 is markedly larger than that of model(7): 0.0439, however the Ad addition takes place on the carbon(7). The OP represents interaction strength of each PIO, however, it does not give the information of contribution of each atom in the PIO. We can see the contribution of each atom in the PIO by examining the atomic compositions of the PIO. Atomic compositions of the PIO-1 of each addition are summarized in Table 10. In the case of Ad addition, the contribution of the addition center carbon of model(7) is 30.1%, the largest one, whereas those of model(14) and model(9) are only 19.7% and 16.4%, respectively. The addition center is more localized on the carbon(7) than on the carbon(14) and carbon(9). Moreover, the contribution of the addition center carbon of model(8) is only 9.9%. This is the reason why though having the carbon(8) with the largest minus charge, the Ad addition to the carbon(8) is unfavorable.
In the case of CH3- addition, the atomic composition of model(23) is 12.1%, the largest, wheras that of model(22) is 8.9%. In the case of .CH3 addition, the atomic compositions of model(9), model(10), model(14) and model(23) are almost the same as each other. Taking account of the atomic composition of PIO, we can more properly estimate the regioselective addition of Ad, CH3-, DBEM- and .CH3* by using the OP of PIO-1 of them.

Table 10. Atomic compositions (%)a) of PIO-1 of each addition
Model X b):AdCH3-.CH3
(6)_X27.8
(7)_X30.1
(8)_X9.9
(9)_X16.47.914.6
(10)_X16.28.715.9
(13)_X26.6
(14)_X19.78.514.8
(22)_X8.9
(23)_X12.115.3
a) an average value of a-spin and b-spinb) Model(m) means carbon number in Figure 1 and X means Ad, CH3-, DBEM, or CH3..

7 Conclusion

Regioselective addition of Ad, CH3- ,DBEM- and .CH3 to La@C82 was studied by using PIO analysis proposed by Fujimoto et al. [2].
PIO analysis of the above addition reactions reveals that the interaction between the reagents and the La@C82 is mainly expressed in only one PIO, PIO-1.
The OP of the POI-1 can be a proper reactivity index of the regioselective addition by taking account of the atomic composition of the PIO.
The orders of regioselectivity of each addition is as follows; 1) Ad addition: (7) > (6) > ... , 2) CH3-, DBEM- addition: (23) > (22) > ... , 3) .CH3 addition: (9) (10) ... . The order of Ad and DBEM- addition are coincident with the experimental results.
PIO analysis based on extended Huckel MOs is a powerful tool for prediction of the reactivities of the reaction systems large in size, such as ones including EMFs.

8 Appendix

Table 11. Extended Huckel Parameters
Orbital ; ( Hii(eV) , z1,)
H1s;( -13.60, 1.300),
C2s;( -21.40, 1.625), C2p;( -11.40, 1.625),
O2s;( -32.30, 2.275), O2p;( -14.80, 2.275),
Br4s;( -24.00, 2.638), Br4p;( -12.30, 2.257),
La6s;( -5.577, 1.533), La6p;( -3.480, 1.198), La5d;( -7.88, 2.274)

References

[ 1] L. Feng, T. Nakahodo, T. Wakahara, T. Tsuchiya, Y. Maeda, T. Akasaka, T. Kato, E. Horn, K. Yoza, N. Mizorogi, S. Nagase, J. Am. Chem. Soc., 127, 17136 (2005).
 [ 2] (a) H. Fujimoto, T. Yamasaki, H. Mizutani, N. Koga, J. Am. Chem. Soc., 107, 6157 (1985).
(b) H.Fujimoto, Acc. Chem. Res., 20, 448 (1987).
 [ 3] H. Fujimoto, Y. Nakano, K. Fukui, J. Mol. Struct., 300, 425 (1993).
 [ 4] A. Ogawa, M. Tacibana, M. Kondo, K. Yoshizawa, H. Fujimoto, R. Hoffmann, J. Phys. Chem. B, 107, 12672 (2003).
 [ 5] http://classic.chem.msu.su/gran/gamess/index.old.html
[ 6] (a) A. Shiga, Y. Kurusu, J. Mol. Catal., A 241, 199 (2005).
(b) T. Motoki, A. Shiga, J. Computational Chem., 25, 106 (2004).
 [ 7] (a) A. Shiga, J. Kojima, T. Sasaki, Y. Kikuzono, J. Organometal. Chem., 345, 275 (1988).
(b) A. Shiga, H. Kawamura, T. Ebara, T. Sasaki, Y. Kikuzono, J. Organometal. Chem., 366, 95 (1989).
(c) A. Shiga, J. Macromol.Sci., A34, 1867 (1997).
(d) A. Shiga, J. Molecular Catal., A 146, 325 (1999).
(e) J. Handzlik, A. Shiga, J. Kondziolka, J. Mol. Catal., A 284, 8 (2008).
 [ 8] T. Motoki, A. Shiga, J. Computational Chem., 25, 106 (2004).
LUMMOXTM, sold by Ryouka System Inc. (Tokyo), correspondent K. Chiba e-mail:chibao@rsi.co.jp
http://www.rsi.co.jp/kagaku/cs/pio/index.html

Return