Three-component reaction involving isoquinoline and dimethyl acethylenedicarboxylate in the presence of indole: Theoretical and experimental investigations of the reaction mechanism

Abstract

The reaction kinetics among isoquinoline, dimethyl acetylenedicarboxylate, and indole (as NH-acid) were investigated using ultraviolet (UV) spectrophotometry. The reaction rate equation was obtained, the dependence of the reaction rate on different reactants was determined, and the overall rate constant (k_ov) was calculated. By studying the effects of solvent, temperature, and concentration on the reaction rate, some useful information was obtained. A logical mechanism consistent with the experimental observations was proposed. Also, comprehensive theoretical studies were performed to evaluate the potential energy surfaces of all structures that participated in the reaction mechanism. Finally, the proposed mechanism was confirmed by the obtained results and the probable and logical reaction paths and also a correct product configuration were suggested based on the theoretical results.

Keywords

NH-acid ultraviolet spectrophotometry isoquinoline kinetic parameters acetylenic ester theoretical study

Introduction

Multicomponent reactions (MCRs) have been recognized as strong tools in modern drug discovery procedure since these types of reactions make it easy and rapid to access the structural variation and complexity of compounds through single-step conversions.^1,2 Such reactions are rapid, well-diversified, efficient, atom-economic, and environment-friendly.^3,4 Due to the advantages mentioned above, developing new MCRs with environmentally benign protocols has been recognized as one of the most important topics in synthetic chemistry.^5–7 MCRs display many of the most desired features for an organic synthesis,^8–10 and they have been intensively studied in recent years.^11–14 Heterocyclic rings are the fundamental components in the skeleton of more than half of the biologically active compounds produced by nature.¹⁵ As an example, aza-aromatic compounds that are activated by acyl chlorides or dimethyl acetylenedicarboxylate (DMAD) are important intermediates for the synthesis of a variety of biologically active nitrogen-containing alkaloids.^16,17 This is why great efforts have been made to discover and optimize methods to construct new heterocycles.

Isoquinoline skeleton can be found in a large number of naturally occurring^18–21 and biologically synthesized active heterocyclic compounds.^22,23 In particular, 1,2-dihydroisoquinoline derivatives act as delivery systems that transport drugs through the highly impermeable blood−brain barrier.^24,25 These compounds also exhibit sedative,²⁶ antidepressants,²⁷ antitumor, and antimicrobial activities.^28,29 With wide kinetic studies been done in the synthesis of phosphorus ylides,^30–34 the structural effect of the reactant (different NH-acids and dialkyl acetylenedicarboxylates) on the reaction rate has been well investigated. Changes in the reactant structure that participates in the reaction, eventuates changes in the overall reaction rate. Evaluating these changes could reveal many facts about the confirmation of the reaction mechanism. A number of our projects, including the theoretical and experimental studies of chemical reaction kinetics, have been completed.^35–41 Now, by the synthesis of the novel compounds via a pseudo-three-component reaction, in the presence of isoquinoline 1, DMAD 2, and different NH-acids 3,⁴² an opportunity has been achieved to investigate the effect of changing the structure of the first component participating in the reaction (i.e. isoquinoline instead of triphenylphosphine) along with the effect of changing the structure of the NH-acids on the kinetic parameters and reaction mechanism. Using available information and based on well-established chemistry of trivalent phosphorus and nitrogen nucleophiles,^43–46 it is reasonable that the product results from initially adding the isoquinoline to DMAD (2) and subsequent protonation of the 1:1 adduct by the NH-acid 3 to form the product (Figure 1). To gain further insight into the reaction mechanism, a kinetic study, including determining the preferred kinetic path, investigation of intermediate structures and transition states in the reaction path, determination of the kinetic and thermodynamic stabilities of products, recognition of rate-determining step, calculation of the reaction rate and finally, identifying and confirming the reaction mechanism are the issues that have been first explored for this class of reactions using ultraviolet (UV) spectrophotometry. Also, another goal of this project was to investigate the effect of substituted groups on the structure of NH-acid. Therefore, to investigate the effect of the substituted methyl group as an electron donor group, the results of this study will be compared with the results obtained in the presence of NH-acid 3-methyl indole³⁶ and will be examined.

Figure 1.

The reaction among isoquinoline 1, dimethyl acetylenedicarboxylate 2, and indole 3 for the synthesis of N-amino esters.

Materials and methods

DMAD, isoquinoline, and indole, purchased from Fluka (Buchs, Switzerland), were used without further purifications. All the other pure solvents including dichloromethane and acetone were purchased from Merck (Darmstadt, Germany). A Cary UV/Vis spectrophotometer model Bio-300 (Australia) was also applied throughout this work. All geometrical structures were optimized at B3LYP/6-311++G (2d, 2p) level of theory using the Gamess suite software package.⁴⁷ The corresponding frequencies of the structures were estimated at the same level of theory to check the stationary points without imaginary frequencies and the transition states with only one imaginary frequency. Also, the intrinsic reaction coordinate (IRC) approach^48,49 was followed to ensure that the given transition state is related to the corresponding reactants and products. The solvent effect on complex stability was examined using the B3LYP/6-311++G (2d, 2p) method by applying the polarizable continuum model (PCM).

Results and discussion

Experimental studies

The first step in investigating the reaction kinetics is to determine the appropriate wavelength. So, for this purpose, the UV spectra of 3-mM dichloromethane solution of compounds 1−3 were recorded separately over the wavelength range of 190–400 nm (Figure 2). In the second step, the UV spectrum was recorded for the reaction among a higher concentration of 1 (30 mM) with 2 and 3, every 25 min (Figure 3). A comparison of the results obtained in Figures 1 and 2 shows that the spectrum recorded in the range of 330–380 nm is part of the product 4 spectrum where the reactant lacks any absorption spectrum. It should be noted that a specific concentration of reactants must be used, which provides a normal absorption of the product (0.8–1.2) in the wavelength range of 330–380 nm. Selecting 370 nm in the wavelength range of 330–380 nm, which has a good absorption value (0.9), would be a suitable choice for kinetic measurements.

Figure 2.

The plot of absorbance versus wavelength of the solution 3-mM isoquinoline 1 (red line), dimethyl acetylene-dicarboxylate 2 (green line), and indole 3 (blue line) in dichloromethane.

Figure 3.

Expanded section of UV spectra (an increasing absorption) over the wavelength range 280–400 nm for the reaction among 1 (30 mM), 2 (3 mM), and 3 (3 mM) in dichloromethane.

The reaction kinetics were then followed by plotting UV absorbance versus time for the reaction among 1 (30 mM), 2 (3 mM), and 3 (3 mM) at 15.0°C. The time dependency of this absorbance is represented in Figure 4. The absorbance at reaction completing time (A_∞) can be obtained from this plot at t = 950 min. Under this condition, the rate equation can be expressed as

Rate = k_{obs} {[2]}^{β} {[3]}^{γ} and k_{obs} = k_{ov} {[1]}^{α}

(1)

Figure 4.

The experimental absorbance changes (dotted line) versus time accompanied by the second-order fit curve (solid line) for the reaction among 1 (30 mM), 2 (3 mM), and 3 (3 mM) at 15.0°C and 370 nm in dichloromethane.

Using nonlinear regression analysis,⁵⁰ the data of original experimental absorbance versus time provided a pseudo-first-order curve (solid line) that exactly coincides with the experimental curve (dotted line), as shown in Figure 4. The pseudo-first-order rate constant (k_obs) was also calculated at 15.0°C. The observed pseudo-first-order rate constant (k_obs) indicates that the order of reaction concerning the sum of compounds 2 and 3 is equal to one (β + γ = 1).

Continuing the experiments, the reaction among 1 (3 mM), 2 (30 mM), and 3 (3 mM) at 15.0°C was conducted, in which the rate equation can be expressed as

Rate = k'_{obs} {[1]}^{α} {[3]}^{γ} and k'_{obs} = k_{ov} {[2]}^{β}

(2)

The data of original experimental absorbance versus time were used to provide a pseudo-first order fit curve (solid line) at 370 nm, which exactly coincides with the experimental curve (dotted line). In this case, the observed pseudo-first order rate constant $(k'_{obs})$ indicates that the order of reaction concerning the sum of compounds 1 and 3 (α + γ = 1) is equal to one. Given (β + γ = 1) and (α + γ = 1), it is obvious that γ = 0 and the order of the reaction with respect to indole 3 is equal to zero. The above results show that the overall reaction order is equal to two (α + β + γ = 2).

To further confirm the above results, all kinetics were studied in the presence of isoquinoline 1 at different concentrations (ranging from 3.6 × 10⁻² to 5.4 × 10⁻² M) with a constant concentration of 2 and 3 (3 mM). The straight line obtained from the plot of observed pseudo-first-order rate constant (k_obs) versus [1] (Figure 5) proved that the reaction is of first-order concerning compound 1 (α = 1). The overall rate constant (k_ov) was calculated based on the slope of the line in Figure 5 and was reported in Table 1. The same results obtained in the presence of DMAD 2 at different concentrations (ranging from 3.6 × 10⁻² to 5.4 × 10⁻² M) indicating the reaction is of first-order to the compound 2 (β = 1). As a result, since (β + γ = 1) and (α + γ = 1) (as determined previously), it is obvious that γ = 0 and order of reaction to indole 3 must be equal to zero. All the above results indicate that the overall reaction order is equal to two (α + β + γ = 2).

Figure 5.

Observed rate constant $(k_{obs})$ versus concentration of reactant 1 at 370 nm and 15.0°C in dichloromethane. [1] = 3 × 10⁻²–5.4 × 10⁻² M and [2] = [3] = 3 × 10⁻³ M.

Table 1.

Values of the overall second-order rate constant (k₂) for the reactions of 1–3 in the presence of solvents such as dichloromethane and acetone, respectively, at all temperatures, investigated.

Solvent	k₂/M⁻¹ min⁻¹
	15.0°C	20.0°C	25.0°C	30.0°C
Dichloromethane	5.5 ± 0.1	8.4 ± 0.2	12.2 ± 0.4	17.8 ± 0.4
Acetone	6.4 ± 0.1	9.0 ± 0.3	13.3 ± 0.2	19.1 ± 0.3

To determine the effects of temperature and solvent environment on the reaction rate, a set of experiments was conducted at different temperatures and solvent polarities with other conditions being the same as those of the previous experiments. For this purpose, acetone, with the dielectric constant of 20.7 was replaced by dichloromethane with a dielectric constant of 8.9. Based on the results, the rate of reaction was accelerated in an environment with a higher dielectric constant (acetone). Also, at all the investigated temperatures, the results indicate an increase in reaction rates at higher temperatures (Table 1). In the studied range of temperatures, the reaction activation energy was calculated from the slope of the pseudo-second-order rate constant (log k_ov) versus reciprocal temperature (1/T).

Rate equation

Based on the above results, a simplified proposed reaction mechanism is depicted in Figure 6. Unlike the experimental results, it may be assumed that the second step is the rate-determining step for the proposed mechanism. In this case, the rate equation can be expressed as follows

Rate = \frac{k_{2} k_{3} [1] [2] [3]}{k_{- 2}}

(3)

Figure 6.

Proposed mechanism for the reaction among 1, 2, and 3 for the generation of N-amino esters 4.

The final equation indicates that the overall order of the reaction is three which is not compatible with experimental overall order of reaction (= 2). Besides, according to this equation, the order of reaction to indole 3 is one, whereas it was shown to be equal to zero. For this reason, it appeared that the second step is fast.

If we assume that the third step (rate constant k₄) is the rate-determining step for the proposed mechanism, in this case, there are two ionic species to consider in the rate-determining step (I2 and N⁻) (Figure 6), have full positive and negative charges and form very powerful ion–dipole bonds to the acetone as a higher dielectric constant solvent. However, the transition state for the reaction between two ions carries a dispersed charge and bonding of solvent (acetone) to this dispersed charge would be much weaker than the reactants. The solvent, thus stabilizes the ionic species, more than it would the transition state. Therefore, the activation energy could be raised and slow down the reaction. However, in practice, acetone speeds up the reaction and for this reason, the third step could not be the rate-determining step. Furthermore, the rate law can be expressed as

Rate = k_{2} k_{3} k_{4} [1] [2] [3]

(4)

The final equation also indicates that the overall order of the reaction is three, and the order of reaction for indole 3 is one, which is not compatible with experimental results.

If the first step (rate constant k₂) was the rate-determining step, in this case, bonding of solvent to the dispersed charge in the transition state is much stronger than the uncharged reactants. The solvent thus stabilizes the transition state more than it does the reactants and speeds up the reaction. Our experimental results show that the solvent with higher dielectric constant exerts a powerful effect on the rate of reaction (see Table 1). Results indicate that the first step (k₂) of the proposed mechanism could be the rate-determining step. However, a good kinetic description of the experimental result using a mechanistic scheme based upon the steady-state approximation is frequently taken as an evidence of its validity. By application of this, the rate of the formation of product 4 from the reaction mechanism (Figure 6) is given by

\frac{d [4]}{dt} = \frac{d [product]}{dt} = rate = k_{4} [I_{2}] [N^{-}]

(5)

By replacing the value of [I1] and [I2]

Rate = k_{2} [1] [2]

(6)

This equation is compatible with the results obtained by UV spectrophotometry. For equation (6) that shows the overall reaction rate, the activation parameters involving ΔG^#, ΔS^#, and ΔH^# could be now calculated for the first step (rate-determining step), as an elementary reaction, based on Eyring equation. The results are reported in Table 2.

Table 2.

Activation parameters for the reactions 1–3.

Solvent	Ea (kJ mol⁻¹)	ΔH^# (kJ mol⁻¹)	ΔS^# (J mol⁻¹ K⁻¹)	ΔG^# (kJ mol⁻¹)
Dichloromethane	56.6 ± 2.3	54.2 ± 2.2	−76.5 ± 2.8	76.8 ± 2.7
Acetone	53.3 ± 1.8	50.8 ± 1.8	−86.9 ± 2.7	76.5 ± 2.6

The activation parameters (ΔG^#, ΔS^#, and ΔH^#) of the first step (rate-determining step), as an elementary reaction, can be calculated based on the Eyring equation. The results are reported in Table 2.

Theoretical study

The overall order of reaction and the order of reaction to different reactants, solvent, and temperature effect and also the detection of rate-determining step, and finally, the reaction mechanism as well are determined in the experimental section. Now, to gain more insight into the reaction mechanism, we intend to consider the mechanism of this reaction using theoretical calculations.

In this section, we will evaluate the potential energy levels, different kinetic paths, and energy levels of different products. Also, we intend to determine the preferred thermodynamic product, the kinetic preferred path, and to confirm the rate-determining step, as outlined in the experimental section. All structures were optimized by density functional theory at B3LYP/6-311++G (2d, 2p) level of theory. Two different pathways (a and b) were predicted for this reaction. Potential energy profiles for the pathways a and b are presented in Figure 7. Also, the optimized geometries of all structures are presented in Figure 8. As can be seen, the reaction path b includes four steps. The first step initiates with a nucleophilic attack of atom N17 in the structure of isoquinoline (R2) to atom C8 of DMAD (R1) to form I1, via passing through Ts1. The second step includes an intramolecular proton transfer from atom C32 to atom C9, via Ts2b, to form I2b. Then, I3b structure can be formed by the hydrogen bonding interaction of I2b and R3. The third step includes an intermolecular proton transfer from atom N34 of indole (R3) to atom C32 in the structure of Ts3b, to form I4b. In the fourth step, once the indole ion attacks atoms C18 or C21, stable products, P1 and P2 are formed, respectively. The fourth step proceeds without any activation barrier. This result seems to be due to the high tendency of indole ion for bonding to C18 or C21. The two pathways, a and b, also lead to the formation of identical products. The path a includes three steps with the first step being the same as that of the path b and the second step including a direct intermolecular proton transfer from atom N34 in indole (R3) to atom C9 within the structure of Ts2a, to form I3a. The very high tendency of the indole ion towards forming N34-C18 or N34-C21 bonds caused the third step of the reaction to proceed without producing any activation barrier, as same as the path b.

Figure 7.

The potential energy profile of the reaction in paths a and b (data of dichloromethane are in parentheses and data of acetone are in italic form).

Figure 8.

Optimized geometries of all structures in pathways a and b of the reaction accompanied by geometrical parameters.

Investigation of the potential energy surfaces in the pathways a and b shows that the first step of the two pathways a and b with an energy barrier of 70.68 kJ/mol was recognized as the rate-determining step. Because this step of the reaction is a common step for two kinetic paths a and b, no kinetic preference was observed in these pathways. Also, because both kinetic paths lead to the same products (P1 or P2), no thermodynamic preference was observed. Despite these, the path a, because of less exhausted energy, is more preferable. In addition to the experimental results which imply the absence of the product P2, theoretical studies also show that the product P2 is considerably more unstable than P1, by the value of 103.32 kJ/mol, which can be a factor in the absence of this product.⁴² Also, the geometrical structures in Figure 8 show that C18 is more accessible for nucleophilic attack in the last step of the reaction.

To check the effect of solvent on the potential energy surfaces, dichloromethane was subjected to condensed phase calculations carried out using the PCM. As shown in Figure 6, the energy barrier of the first and the third step of the pathway b has decreased, but that of the second step considerably increased, both in dichloromethane and acetone. In this condition, the second step of the reaction is determined to be the rate-determining step which is in contrast to what was obtained in the experimental section and the gas phase. Both energy barriers of the pathway a decreased and results show that the first step of the reaction is determined to be the rate-determining step. Results indicate that pathway a has a considerable kinetic preference rather than that of b. Also, the only path is compatible with the increasing of the reaction rate with increasing of dielectric constant of the solvent in experimental results. Similar results were obtained from theoretical investigations. The product P2 is more unstable than the P1, both in dichloromethane and acetone, by the value of 89.96 and 90.11 kJ/mol, respectively and it suggests that products P2 cannot exist. Although the calculations in the gas phase do not show any superiority between the two paths, but in the solution phase, the path a is preferable because of the matching with experimental results and kinetic superiority to path b. Theoretical results illustrate how the first step of the reaction is the rate-determining step and show the effect of solvents with higher dielectric constant in reducing the energy barrier of this reaction.

Conclusion

Mechanistic investigation of the reaction among isoquinoline 1, DMAD 2 in the presence of NH-acid, such as indole 3 was investigated both experimentally and theoretically. The results are summarized below.

The kinetics of the reaction were studied by measuring the increase in the product absorption rate. The overall reaction order follows the pseudo-second-order kinetics with the first-order dependent on the concentration of isoquinoline 1 and DMAD 2. The reaction activation parameters were calculated by measuring the rate constant at different temperatures. Recognizing the first step of the reaction (k₁) as the rate-determining step, a mechanism was proposed for the reaction based upon experimental and theoretical results. For all reactions, the rates were increased by increasing the dielectric constant of the solvent. The effect of temperature on the reaction rate is much more than that of the solvent dielectric. The computational results in solution do not show any thermodynamic priority between the two studied kinetic pathways, but the pathway a, because of the more kinetic preference, is more preferable than the pathway b. Comparing the results obtained in this project with the results obtained in our previous project in the presence of isoquinoline 1, DMAD 2, and 3-methyl indole 3³⁶ shows that as expected, the reaction in the presence of the methyl electron donor group has more activation energy. The electron donor groups reduce the acidic strength of NH-acid. Therefore, it has slowed down the overall reaction rate.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors sincerely thank Payame Noor University for providing financial support for this work.

References

Dömling

Chem Rev 2006; 106: 17–89.

Lieby-Muller

Constantieux

Rodriguez

J Am Chem Soc 2005; 127: 17176–17177.

Pando

Stark

Dankert

, et al. J Am Chem Soc 2011; 133: 7692–7695.

Shaabani

Rezayan

Keshipour

, et al. Org Lett 2009; 11: 3342–3345.

Ganem

Acc Chem Res 2009; 42: 463–472.

Candeias

Veiros

Carlos

, et al. Eur J Org Chem 2009; 12: 1859–1863.

Ning

Jiang

Zhang

, et al. Green Chem 2010; 12: 357–1361.

Wender

Haddy

Wright

DL.

Chem Ind 1997; 765: 767–769.

Maghsoodlou

Habibi-Khorassani

Rofouei

, et al. Arkivoc 2006; xii: 145–151.

10.

Dastoorani

Khalilzadeh

Khaleghi

, et al. New J Chem 2019; 43: 6615–6621

11.

Ugi

Dömling

Werner

J Heterocyclic Chem 2000; 37: 647–658.

12.

Yavari

Anary-Abbasinejad

Monatshefte

AA.

Fur Chemie 2002; 133: 1331–1336.

13.

Adib

Yavari

Mollahosseini

Tetrahedron Lett 2004; 45: 1803–1805.

14.

Nair

Sreekanth

Abhilash

, et al. Synthesis 2003; 12: 1895–1902.

15.

Katritzky

Rees

CW.

Comprehensive heterocyclic. New York: Pergamon Press, 1984.

16.

Yavari

Hossaini

Sabbaghan

, et al. Tetrahedron 2007; 63: 9423–9428.

17.

Nair

Devipriya

Eringathodi

Tetrahedron Lett 2007; 48: 3667.

18.

Díaz

Miguel

Lavilla

J Org Chem 2004; 69: 3550–3553.

19.

Adib

Mollahosseini

Yavari

, et al. Synthesis 2004; 06: 861–864.

20.

Yavari

Hossaini

Sabbaghan

Tetrahedron Lett 2006; 47: 6037–6040.

21.

Gilchrist

TL.

Heterocyclic chemistry. New York: Wiley, 1985.

22.

Michael

JP.

Nat Prod Rep 2002; 19: 742–760.

23.

Scott

Williams

RM.

Chem Rev 2002; 102: 1669–1730.

24.

Sheha

El-Koussi

Farag

Arch Pharm 2003; 336: 47–52.

25.

Mahmoud

Aboul-Fadl

Shehq

, et al. Arch Pharm 2003; 336: 573–584.

26.

Lukevics

Segal

Zablotskaya

, et al. Molecules 1997; 2: 180–185.

27.

Maryanoff

McComsey

Gardocki

, et al. J Med Chem 1987; 30:1433–1454.

28.

Tietze

Rackehmann

Miller

Chem Eur J 2004; 10: 2722–2731.

29.

Knölker

Agarwal

Tetrahedron Lett 2005; 46: 1173–1175.

30.

Habibi-Khorassani

Maghsoodlou

Ebrahimi

, et al. J Solution Chem 2007; 36: 1117–1127.

31.

Habibi-Khorassani

Ebrahimi

Maghsoodlou

, et al. Curr Org Chem 2011; 15: 942–952.

32.

Habibi Khorassani

Maghsoodlou

Ebrahimi

, et al. Prog React Kinet Mech 2005; 30: 127–143.

33.

Habibi-Khorassani

Maghsoodlou

Zakarianezhad

, et al. Heteroat Chem 2008; 19: 723–732.

34.

Zakarianezhad

Habibi-Khorassani

Ebrahimi

, et al. Heteroat Chem 2010; 21: 462–474.

35.

Zakarianezhad

Makiabadi

Habibi-Khorassani

, et al. Res Chem Intermed 2018; 44: 2653–2665.

36.

Zakarianezhad

Masoodi

Shool

, et al. Int J Chem Kinet 2016; 48: 770–778.

37.

Zakarianezhad

Shool

Iran Chem Commun 2017; 5: 301–307.

38.

Zakarianezhad

Mohammadi

Phys Chem Res 2016; 4: 221–230.

39.

Zakarianezhad

Masoodi

Shool

Phosp Sulf Silicon 2016; 19: 1063–1069.

40.

Zakarianezhad

Mohammaddadi

J Sulf Chem 2015; 36: 422–433.

41.

Mohtat

Siadati

Khalilzadeh

MA.

Prog React Kinet Mech 2019, 44: 213–221.

42.

Nassiri

Maghsoodlou

Heydari

, et al. Mol Divers 2008; 12: 111–117.

43.

Cadogan

JIG

. Organophosphorus reagents in organic synthesis. New York: Academic Press, 1979.

44.

Engel

Synthesis of carbon-phosphors bonds. Boca Raton, FL: CRC Press, 1988.

45.

Hudson

. The chemistry of organophosphorus compounds. In: Hartley

(ed.) Primary, secondary and tertiary phosphates and heterocyclic organophosphorus (III) compounds. New York: Wiley, 1990.

46.

Corbridge

DEC

. Phosphorus: An outline of chemistry, biochemistry, and uses. Amsterdam: Elsevier, 1995.

47.

Frisch

Trucks

Schlegel

, et al. GAUSSIAN 09. Pittsburgh, PA: Gaussian Inc., 2009.

48.

Schmidt

Baldridge

Boatz

Elbert

Gordon

Jensen

Koseki

Matsunaga

Nguyen

Windus

Dupuis

Montgomery

. J Comput Chem 1993; 14: 1347–1363.

49.

Gonzalez

Schlegel

HB.

J Phys Chem 1990; 94: 5523–5527.

50.

Schwartz

Gelb

RI.

Anal Chem 1978; 50: 1592–1594.