Molecular dynamics tests of the Smoluchowski–Collins–Kimball model for fluorescence quenching of spherical molecules

Abstract

We test the Smoluchowski–Collins–Kimball (SCK) model of fluorescence quenching reaction in liquids. Our attention is focused on the description of diffusion controlled processes and we use the simplest microscopic model of binary de-excitation which occurs instantaneously. Molecular dynamics simulations have been performed for a wide range of densities and for various potentials of interparticle interactions. The large number of particles used (typically N = 681 472) allowed us to obtain quantitative results. The simulations show that at very short times the SCK model completely fails, especially if the distribution function of the reagents is not included in the model. Even if the short time data are excluded the diffusion constant obtained by fitting the simulation data with the SCK model may still be burdened with a 25% error when the distribution function is not taken into account. The error can be reduced to 10% level if this function is included in the model.

---

1. Introduction

Typical fluorescence quenching reaction for the A molecule consists of three processes: initial excitation A → A*, spontaneous de-excitation A* → A and de-excitation in the presence of a quencher, B, A* + B → A + B. In this work our attention is focused only on the third reaction, which is diffusion controlled. There are many theoretical descriptions of diffusion controlled reactions.^1–9 The most popular one is the Smoluchowski–Collins–Kimball (SCK)^2,3 model based on the Smoluchowski equation.^1,2 It gives a simple analytical formula for the reaction rate, which is widely used to interpret experimental data of fluorescence quenching.^10–14 The model returns the parameters of the reaction: intrinsic rate constant (k₀), the reaction radius (a) and the sum of the diffusion coefficients of the excited molecule and the quencher (D_S). Many experiments^15–18 show inconsistencies between the model and reality. The recent works of Krystkowiak and Maciejewski^17,18 report significant differences between the parameters obtained by fitting the experimental results of fluorescence quenching with the SCK model and the values obtained by different methods. On the other hand, the report on consistency between the model and experimental results can be also found in the literature.¹⁹

The purpose of our work is to test the SCK model applying it to fit simulation data of fluorescence quenching. We simulate reaction between A* and B using molecular dynamics (MD) and our assumptions match as closely as possible those of the SCK model: the microscopic model for reaction between A* and B molecules is the same as in the SCK, the reactants are spherical and the concentration of quencher is very low. Because most of the assumptions of the theory are satisfied the comparison of theory with simulation results gives us information on the scale and the origin of the inconsistencies coming from the application of the Smoluchowski equation. The simulations also allow us to estimate the errors caused by the neglect of spatial correlations between A* and B (g_A*B(r) = 1) typical in analysis of experimental data. Our simulations have been performed for a wide range of the parameters (the liquid density, D_S, a) and different types of intermolecular interactions so we believe that random coincidences are absent. The fact that the conditions of our computer experiments are much closer to the assumptions of the SCK model than it is for a real experiment allows us to treat the scale of deviations between the simulation results and the model as a minimum level of errors expected in applications of SCK to real experiments. Focusing our attention on the errors within the SCK model we have not concerned more realistic descriptions of the reaction cross section. Nevertheless, some of our results are also valid for the step function nonradiative lifetime model^2,20 (SFNL) which assumes more realistic reaction cross section than the SCK model.

A large scale of our simulations (typically, the total number of molecules N = 681 472) has allowed us to make the quantitative test. Simple, qualitative simulations on both the Smoluchowski^1,2 and the SCK model were performed over 10 years ago for a hard sphere liquid.^21,22 The main conclusion was²² that the SCK model gave surprisingly good results, much better than the Smoluchowski one. Zhou and Szabo²² showed also that the agreement between the model and the simulations improved significantly when the shape of the two-particle radial distribution function was taken into account. This conclusion was only qualitative one and no quantitative analysis of differences between the model and the simulations have been performed. Ten years ago the scale of simulations was not sufficient for such an analysis. As a result, the scale of errors resulting from the standard assumption that the liquid is structureless as well as from the simplifications of the SCK model itself is still unknown.

The paper is organized as follows. Section 2 defines the model, presents the computational details and methods of describing the results. Section 3 gives the results of simulations and discussions. It starts with the analysis (Section 3.1) of errors in results presented in Section 3.2. Section 3.2, describes an “idealized” situation when the time of measurement is long so we can divide the obtained intensity curve and discard a part of data inconsistent with the remaining ones. As a result, we are able to extract the errors related to the Smoluchowski equation for times when the relative motion of reactants is dominated by diffusion. We also show the errors resulting from the standard simplification that the liquid is structureless. The “idealized” situation allows us to generalize some of results to the case when the reaction is described by the SFNL model. Section 3.3 is concerned with the problem of the early phase of fluorescence quenching. We discuss here the influence of non-diffusive motion which dominates at short times. We also consider situation frequent in experiments on fluorescence quenching: the time of measurement is short and data cannot be divided into parts corresponding to short time processes and purely diffusive motion (long times). An influence of such situation on the parameters obtained from the fitting is presented. The summary and final remarks are presented in Section 4.

2. Model, computational aspects and interpretation of results

Following the SCK model, we consider spherical molecules in three dimensional bulk liquid and a fluorescence quenching reaction in the form:

A* + B → A + B + hν (2.1)

where A* is an excited A molecule and B is a quencher. We assume that both the excitation and the de-excitation (by the emission of photon with energy hν) do not influence either the motion nor the interaction of A with its surroundings. We also assume that the de-excitation process is instantaneous and, following the previous simulations on diffusion controlled reactions,^21,22 the reaction always occurs if the distance between A* and B becomes equal to the reaction radius, a. Therefore, the probability of the reaction, P, can be written as:If the concentrations [A*](t) and [B] are sufficiently small we can define the time-dependent rate constant of the reaction (2.1), k(t), by:

According to the SCK model, the probability of finding the A* particle around the B particle at the interparticle distance r and at time t, p(r,t), satisfies the Smoluchowski equation^1,2 for spherical molecules:

with the initial condition: p(r,0) = g_AB(r), where g_AB(r) is the equilibrium pair distribution function which, according to our assumptions, is identical to g_A*B(r). In the SCK model D_S is independent of r and equal to the sum of the diffusion coefficients of bulk liquid. The model also assumes the radiation boundary condition of Collins and Kimball:^2,3where k₀ is a constant. eqns. (2.4) and (2.5) give the following asymptotic formula for the time-dependent rate constant of the SCK model:^2,23 where:

In the case g_AB(r) ≡ 1, eqn. (2.6) is the exact solution of the model for a whole range of t.

The k₀ constant is determined, following Zhou and Szabo,²² from the kinetic theory of gases:^24,25

where T is temperature, k_B is the Boltzmann constant, and m_A, m_B are masses of the particles. The above formula is very general subject to the condition that the liquid is in (or sufficiently close to) the mechanical equilibrium state.

2.1. Simulations

The computer simulations were performed using standard molecular dynamics (MD) constant volume and energy (NVE) method.²⁶ The cubic box and the periodic boundary conditions were applied.²⁶ Most of the simulations were performed for the total number of particles N = 681 472. Several simulations for smaller N (110 592 and 343 000) were also carried out. During a single simulation run the particles interacted via the potential of the same kind. Four various kinds of spherical interparticle potentials were considered. Most of simulations were performed for the two following potentials:

Lennard-Jones²⁶ type (LJ):

where α = 1 and ε, σ are standard parameters of the potential; and the potential with an exponential repulsive term (R6e):

The c_i and d_i constants from eqns. (2.9) and (2.10) were adjusted so the derivative of u(r) (force) was differentiable at R_S.

The cut-off distances, R_C, were low. For the L-J type R_C was always equal to 1.6 (R_S = 1.25). For R6e three values of R_C were applied: 1.65 (R_S = 1.40), 1.95 (R_S = 1.65), 2.15 (R_S = 1.85). Such low values of R_C made the time of computation sufficiently short. Any qualitative influence of the cut-off distance on the investigated problem was not observed.

Two simulation runs for non-attracting potentials were also performed: soft spheres (SS, eqn. (2.9) for α = 1, R_C = 1.15, R_S = 1.0), and hard soft spheres (HSS–(2.9) for α = 5, R_C = 1.03, R_S = 0.985). The Newton equations of motion were solved applying the velocity Verlet algorithm^26,27 with the time step Δt = 0.01σ(m/ε)^1/2 except of HSS for which the time step was four times shorter. The temperature was evaluated directly from the average of the total kinetic energy²⁶ of simulated systems.

For all the investigated systems the reagents (A, A*, and B whose total number was N_R) were mechanically identical (identical ε, σ, and mass m) and differed only by a chemical identification parameter. The reaction (2.1) were realized by simple re-labeling A* to A so the simulated system was always in the mechanical equilibrium state. At the beginning of each evolution (after equilibration) some of the reagents were randomly labeled as the B particles. In order to keep our simulations consistent with the SCK model the considered molar fraction of B was never larger than 0.0025 (see Table 1). The remaining particles were labeled as the A* particles but if the condition in eqn. (2.2) was fulfilled they were instantaneously relabeled as A (the reaction had occurred before the start of the simulation run). The mechanical identity of A and A* allows us to regard the single simulation as a set of N_A* (N_X is the number of X-type particles) independent evolution runs for one excited particle (A*) in a mixture of the B and A particles (A* + A does not lead to de-excitation). The names of the reagents (A, A*, B) were stored in an array (of the length of N_R). Mechanical identity of A and B enabled us to lessen the statistical error by applying procedure similar to that proposed by Gorecki.^28,29 During a single simulation run we considered simultaneously several sets of the reagents (n_set is the number of the sets), each set described by a different array. Therefore, a single simulation run was approximately equivalent to N_Rn_set evolutions for single A* particle.

Table 1 List of all simulations presented in the work. The columns: Run, the number of the simulation run; Type, the symbol of the interparticle potential. The characters followed by R6e for # 7–10 show the value of ε_AC; R_C, cut-off radius; ρ, number density; k_BT, reduced temperature; x_B, molar fraction; N, total number of particles; n_set, number of sets of the B particles taken for the run; t_T, total time of the simulation; D₀, the diffusion coefficient evaluated from eqns. (2.12) and (2.13) ΔD(1.0, 1.13)–deviation of D_fit (≡ D_fit/D₀ − 1) for a₀/σ_AB = 1.0, 1.13 respectively. The index ‘g ≡ 1’ means g_AB(r) ≡ 1, the index ‘SIM’, g_AB taken from the simulation (by eqns. (2.14)). The mark — in the last column means that the value has not been “measured”

Run	Type	R C	ρ	k B T	10^3xB	N	n set	t T	D 0	ΔD(1.0, 1.13)g≡1	ΔD(1.0, 1.13)SIM
a The mean value from the eight independent simulation runs.
#1	LJ	1.60	0.5787	1.00	1.00	681 472	4	464	0.1971	0.165, 0.371	0.084, −0.117
#2	LJ	1.60	0.7830	1.00	1.00	681 472	4	576	0.0794	0.058, −0.009	−0.085, −0.131
#3	HSS	1.03	0.9358	1.00	1.00	343 000	6	576	0.0236	−0.052, −0.122	−0.086, −0.074
#4	SS	1.15	0.9358	1.00	1.50	343 000	6	576	0.0401	−0.002, −0.116	−0.072, −0.088
#5a	LJ	1.60	0.9358	1.00	2.50	681 472	4	656	0.02782	−0.020, −0.101	−0.068, −0.061
#5b	LJ	1.60	0.9358	1.00	1.25	681 472	4	2184	0.02797	−0.032, −0.107	−0.076, −0.069
#5c	LJ	1.60	0.9358	1.00	0.625	681 472	4	738	0.02781	−0.022, −0.101	−0.068, −0.061
#5d–k	LJ	1.60	0.9358	1.00	0.10	110 592	8	1000	0.02763^a	−0.009^a, −0.089^a	—
#6	R6e	2.15	0.9358	1.00	1.00	681 472	4	1080	0.01330	−0.065, −0.106	−0.038, −0.002
#7a	R6e1.0	1.95	0.8594	1.00	0.50	681 472	6	3360	0.00287	−0.198, −0.233	−0.040, 0.030
#7b	R6e1.0	1.95	0.8594	1.00	0.50	681 472	4	2864	0.00287	−0.213, −0.247	−0.056, 0.017
#8a	R6e1.1	1.95	0.8594	1.25	1.50	681 472	2	2048	0.00616	−0.215, −0.266	−0.072, −0.025
#8b	R6e1.1	1.95	0.8594	1.25	0.50	343 000	6	2736	0.00611	−0.201, −0.243	−0.056, −0.001
#9	R6e1.0	1.65	0.8594	1.50	0.50	681 472	6	1312	0.00942	−0.168, −0.235	−0.106, −0.074
#10a	R6e1.5	1.65	0.8594	1.50	0.70	681 472	6	2320	0.00557	−0.218, −0.246	−0.022, 0.026
#10b	R6e1.5	1.65	0.8594	1.50	0.70	343 000	6	3160	0.00559	−0.213, −0.239	−0.019, 0.032

---

From the mechanical point of view the reagent-only system is a simple liquid. To make our computer experiments more realistic we considered also the mixtures composed of the heavy components (10% molar fraction) and the light solvent, C, (90%). The parameters of the heavy particles were: σ_A = σ_B = 1.35σ_C and m_A = m_B = (σ_A/σ_C)³m_C. The σ_AC parameter (identical to σ_BC) of the A–C (or B–C) interaction satisfied the Lorentz–Berthelot rule:³⁰

σ_AC = (σ_AA + σ_CC)/2 (2.11)

The energy parameter ε_XY, was equal to 1.0 for all the A–A, B–B, A–B (the simple liquids as well as the mixtures) and C–C interactions. For a few mixtures the A–C and B–C interactions were strengthened by using the values of energy parameters, ε_AC and ε_BC, higher than 1.0. The numerical values presented below are expressed in reduced units (i.e. for: σ = ε = m = 1.0) of A (or B) for the simple liquids and C for the mixtures. If the other units are used, it is explicitly written in the text.

The reaction condition, eqn. (2.2), was checked at each time step. A comparison with additional simulations for Δt = 0.005σ(m/ε)^1/2 showed that the error resulting from a finite value of Δt was much lower than the errors discussed in Section 3.1. During each simulation run a set of 21 fixed values of a, from σ_AB to over 3σ_AB, was considered and evolutions corresponding to 21 different reaction radii were obtained. The mechanical identity between A and A* enabled us to apply an efficient procedure analogous with the simultaneous simulations for different sets of reagents. In the following we denote the reaction radius for a given simulation as a₀ to emphasis that this is a constant value and a is a parameter of SCK theory and it may be used as an argument of minimization procedure.

According to our assumptions D_S = 2D, where D is the diffusion coefficient of the A (or B) particle of a bulk liquid. The value of the coefficient, denoted by D₀ (indexed by “0” because we treat the value as a constant parameter), was “measured” during simulations by using the Einstein formula:²⁶

where

The limit in (2.12) was obtained by analyzing the time evolution of 〈Δr²(t)〉. This was not difficult because due to high N and the lengths of the simulation runs the fluctuations of 〈Δr²(t)〉/t were small. As an example the plot of 〈Δr²(t)〉/6t as a function of time (the simulation run labelled #10b in Table 1) is shown in Fig. 1. The presented simulation run has been performed for the lowest numbers of reagents (here N_R = 0.1N = 34 300) so the fluctuations are much higher than for a typical run discussed in the work.

Fig. 1 Asymptotic behavior of 〈Δr²(t)〉/6t, eqn. (2.13), as a function of time, t. The data from the simulation run marked by #10b in Table 1. The horizontal line corresponds to the value adopted as D₀ (=0.00559).

The radial distribution function of the reagents was obtained directly from simulation data using the standard formula:²³

where n_R(r₁,r₂) is the total number of pair of the reagents for which the interparticle distance r satisfies r₁ < r < r₂and 〈 〉_t means an average over time. Δr was always equal to 0.0025σ_AB. Typical g_AB(r) function is presented in Fig. 2.

Fig. 2 The partial radial distribution function of the reagents A* and B (g_AB(r), eqn. (2.14)) as a function of the interparticle distance, r. The data from # 8a in Table 1.

2.2. Interpretation of results

In this subsection we present the methods of describing and interpreting the results of simulations. Experimental data of fluorescence quenching are interpreted by analyzing time evolution of the survival probability. For example, the survival probabilities of A* as functions of time obtained from the simulations (solid line) and predicted by the SCK model (dashed line) for D_S = 2D₀ and g_AB(r) ≡ 1 are presented in Fig. 3 (a₀ = σ_AB) and Fig. 4 (a₀ = 3.27σ_AB). The probabilities predicted by the SCK model have been obtained by numerical integration of formula (2.3) with k(t) given by (2.6) (for g_AB(r) ≡ 1, (2.6) is the exact solution of (2.4)). It can be noticed that the inconsistency between the model and the simulations is quite important not only in the vicinity of t = 0, but also for large t which is especially evident in Fig. 3. The differences between the solid and the dashed lines presented in Figs 3 and 4 indicate that there is a significant difference between k_SCK(t) (from (2.6)) and the time dependent rate constant, k(t), characterizing the reaction simulated by MD method. Detailed discussion on the difference is presented in Sections 3.2 and 3.3.

Fig. 3 The survival probability (N_A*(t)/N_A*(0)) as a function of t for a₀ = σ_AB. The solid line–the data from # 8a in Table 1; the dashed line, the SCK model eqn. (2.6) for g_AB(r) ≡ 1, D_S = 2D₀, and a = a₀.

Fig. 4 The survival probability (N_A*(t)/N_A*(0)) as a function of t for a₀ = 3.27σ_AB. Notation as in Fig. 3.

In this work the simulation data are analyzed by calculating the time dependent rate constant, k(t), evaluated from eqn. (2.3) by omitting the higher order terms in δt:

where: ρ is the numerical density (N/V), and x_B is the number fraction of B particles (N_B/N). As an example, Fig. 5 compares k(t) as a function of t obtained from the simulation (empty circles) with k_SCK(t) (dashed line) for the case presented in Fig. 3.

Fig. 5 The time dependent rate constant, k(t), as a function of t for the case shown in Fig. 3. The circles: the results of simulation (evaluated from eqn. (2.15)) the dashed line: the SCK model for g_AB(r) ≡ 1 and D_S = 2D₀; the solid line: g_AB(r) ≡ 1 and D_S taken as double the “idealized” value of D_fit.

In order to describe the differences between the SCK model and the simulations quantitatively, we fitted the values of k(t) obtained from simulation runs using (2.15) to (2.6) by the least square method. The weights of “experimental” points were proportional to N_A*(t). The minimization procedure was applied with respect to one of the parameters (D (≡ D_S/2) at given a or a at given D), or with respect to both a and D. In further considerations, D_fit and a_fit denote the values corresponding to the best fit. The approximation in eqn. (2.15) may lead to considerably large errors if random fluctuations in N_A* are high. But in our case eqn. (2.15) worked very well. We checked that for a wide range of δt (in eqn. (2.15)) its value did not have any significant influence on results of the minimizations.

2.2.1 “Idealized” values. Fig. 5 shows a large inconsistency between the simulation data and the SCK model for short times. This inconsistency is mainly due to non-diffusive processes that are not considered by the SCK model. The problem is described and analyzed in Section 3.3. For very large t the process of A* approaching B can be treated as purely diffusive one. In order to estimate the errors in k_SCK(t) related to diffusive motion only, we introduce the “idealized” values of D_fit and a_fit. These values (see Table 1 and Figs. 6–11) are defined as the limits attained for t_D → ∞, where t_D (discarded time) is the upper limit of the time interval {0, t_D} that was not taken into account in the minimization. The results and discussion on the “idealized” values is given in Section 3.2. As an example, the solid line in Fig. 5 presents the k_SCK(t) curve for a = a₀ (=σ_AB) and D equal to the “idealized” value of D_fit. It is clearly seen that this curve fits the asymptotic data much better than the SCK formula for D₀ (the dashed line). The “idealized” values are defined as the limits for t_D → ∞ and we use such notation being aware that the total time of simulation, t_T, however very long, is finite. The limit has been determined by analyzing the dependence of D_fit on t_D. The problem is that the interval {t_D, t_T} should be sufficiently long to obtain reasonable results of the minimization. Thus, the value of t_D must be significantly lower than t_T. The problem has been overcome by increasing t_T to very large values. The influence of t_T on the determination of the limit is also discussed in Section 3.3.

Fig. 6 The “idealized” values. D_fit/D₀ − 1 at fixed a = a₀ as a function of a₀ for various x_B for #5. Symbols: triangle down–#5a (x_B = 0.0025), square–#5b (x_B = 0.00125), triangle up–#5c (x_B = 0.000625), black circles connected by a line–5d–k (x_B = 0.0001). The error bars have been estimated from the mean square deviation.

For a typical experiment on fluorescence quenching all the measurements from t = 0 to t = t_T are taken into account in the minimization (i.e.t_D = 0). The “idealized” values have been obtained from the minimization over the time interval {t_D → ∞, t_T} and, as a result, they are not burdened with an error due to non-diffusive processes that dominate k(t) for short times. Therefore, the errors of the “idealized” values should be treated as the minimum values of errors of parameters obtained from minimization. This error is not small. In the case presented in Fig. 5 the “idealized” value of D_fit is of about 22% lower than D₀ (see also Section 3.2). In the case of g_AB(r) ≠ 1 eqn. (2.6) is valid only for sufficiently large t. For low D and a close to σ_AB, the convergence of eqn. (2.6) to the exact solution of eqn. (2.4) is poor. But the procedure of determination of the “idealized” values guarantees that the difference between k(t) evaluated from (2.6) and its real value is negligible. This has been confirmed by comparing k(t) obtained from the numerical solution of eqns. (2.4) and (2.5) with k_SCK(t) at a₀ = σ_AB (the worst convergence) for all the simulation points. The highest relative deviations in k_SCK(t) found was only 0.0075 for t = t_T/3 and 0.0035 for t = t_T.

The “idealized” values calculated within the SCK model may be used for the SFNL model, which assumes a more realistic form of reaction than eqn. (2.2). Szabo has shown^2,20 that, for g_AB(r) ≡ 1, the formula for k(t) of SFNL model transforms into that of k_SCK(t) if t/τ → ∞ where τ is the nonradiative lifetime. As a consequence, we generalize some our results concerning the “idealized” values obtained from the SCK model to the SFNL one (Section 3.2).

3. Results of simulations and discussions

The parameters and some results of the simulations discussed in the work are presented in Table 1. The first column (Run) presents the number of the simulation run and the second column (Type) the symbol of the potential used. For the mixtures there are additional characters following the symbol (here R6e), which give the value of ε_AC (e.g. R6e1.1 means: the mixture, the potential (2.10), and ε_AC = ε_BC = 1.1). In the following text we refer to a particular simulation run by its number in Table 1. The last two columns of Table 1 present pairs of values of ΔD (≡ D_fit/D₀ − 1) for a₀/σ_AB = 1.0 and 1.13 (separated by comas) where D_fit are the “idealized” values. The second last column gives the values for g_AB(r) ≡ 1, and the last one: for g_AB(r) calculated from the simulation. All of the evolutions presented in Table 1 have been performed independently except of #7a and #7b. The two results were obtained on the basis of the same evolution run for different moments of initialization of reaction (2.1). The time difference between #7a and #7b is 496 time units.

Further discussion of the simulations results is mainly focused on the minimization with respect to D. The value of D_fit/D₀ − 1 will be considered as a measure of the quality of the SCK approximation. Such an approach is motivated by the fact that the dependence of D on the distance between reactants is also interesting from theoretical point of view (see refs. 31,32). The second reason is given in Section 3.2: for the case of g_AB(r) ≡ 1 (which is a standard assumption for fitting experimental data with eqn. (2.6)) D_fit and a_fit are strongly correlated and, as a result, deviations of D_fit from D₀ give us also information about the deviations of a_fit. For a typical fluorescence quenching the reaction radius, a, is close to σ_AB. Thus, our attention is mainly focused on a < 1.5σ_AB (called further the main region of a).

3.1. Errors of the “idealized” values

The SCK theory assumes that the interactions between the B particles are negligible which is satisfied if the value of x_B is very low. In order to estimate the possible influence of x_B as well as of N (which may always occur in MD simulations) on results eleven evolutions have been performed for the same ρ and k_BT. Eight of them, presented in Table 1 as a single run (#5d–k), were performed for the same x_B and N and differed by the initial conditions only.

First, let us discuss if the value of D₀ obtained in the simulations is reliable. This is an important question because we treat D₀ as a reference value so its error directly burdens the “idealized” values. We have noticed a weak dependence of D₀ on N, but the effect is on the level of statistical errors and its influence can be neglected. One can see this dependence in Table 1 by considering the difference between D₀ for #5a, b, c and D₀ for #5d–k. The statistical error in D₀ for sets #5d–k (sigma of the Gauss distribution) has been estimated from the mean square deviation as 5 × 10⁻⁵. This value is too small to explain the difference between #5a–c and #5d–k by a statistical deviation only. Nevertheless, the difference in D₀ is small and, taking into account that the possible influence of N on D decreases with an increase of N, the effect for N = 681 472 becomes negligible.

The value 5 × 10⁻⁵ seems to underestimate the error of D₀ of #5d–k if compared with the dispersion of the results of #5a, b, c. The difference between D₀ of #5b and #5a,c does not come from different t_T. The values of D₀ of #5b obtained by taking the limit in the Einstein formula for t = 656 (cf. #5a) and t = 738 (cf. #5c) are very close to that for t = t_T (=2184). The dispersion of D₀ for #5a, b, c shows that the random (statistical) deviation of D₀ is more important than the hypothetical dependence on N. Examination of the results for #5a, b, c, #8a, b, and #10a, b enables us to estimate the error in D₀ as not larger than 1%.

The “idealized” values of D_fit obtained by the minimization at the fixed a = a₀ and g_AB(r) ≡ 1 for #5a, b, c and #5d–k are presented in Fig. 6. The range of error bars in the Figure has been evaluated from the mean square deviation of the results of #5d–k. Fig. 6 shows that the dependence of D_fit on x_B, however noticeable, is very weak. For the main region (a < 1.5σ) the maximal difference between the results for a very low x_B (the curve–#5d–k) and for low x_B (the symbols–#5a, b, c) amounts to 0.023 (#5b and #5d–k at a₀ = σ) and, typically it is of about 0.01 - 0.015. Out of the main region, the differences increase with a₀ but they do not exceed 0.06. The A and B particles of the mixtures are about 2.5 times larger than in the case of the simple liquids. But, the simulations for the mixtures have been usually performed for much lower x_B. Therefore, the deviations due to a finite value of x_B should be similar to that presented in Fig. 6. This agrees with the differences in ΔD between #8a (x_B = 0.0015) and #8b (x_B = 0.0005) presented in Table 1.

The random errors (error bars) from Fig. 6 are very small (form 0.004 to 0.010 for a₀ = 1.0, 3.125, respectively) in spite of the total number of the B particles considered in #5d–k is 704 only; it is a few times less than for any of the remaining simulations. Therefore, we can assume that the random (statistical) errors for the remaining simulations are comparable to these from Fig. 6. In order to confirm it, the statistical errors (corresponding to sigma in the Gauss distribution) have been additionally estimated from the mean square deviation of the results of 6 sets of #10a. For a low a the obtained values are only slightly higher than that in Fig. 6. The highest statistical error for a₀ < 1.5σ_AB amounts to 0.0055 for g_AB(r) ≡ 1 and 0.0065 for g_AB(r) taken from the simulation. We should note that the latter value contains also the error resulting from taking the limit for t_D → ∞ because for each set of the reagents the limit has been determined independently. The errors sharply increases with a and for a₀ = 3.27σ_AB (the highest value of a for the mixtures) the estimated statistical error amounts to 0.027. The statistical errors estimated for low a₀ agree with the differences in D_fit between #10a and #10b presented in Table 1. But, the differences between #7a and #7b are nearly three times higher (the differences between #8a and #8b are also high but here x_B is not the same) which shows that the adopted statistical errors may be slightly underestimated.

In the paragraph above we estimated the errors of “idealized” values. Although the estimation is very rough, the values of errors are much smaller than typical values of D_fit/D₀-1 resulting form our simulations. Therefore, we can treat the results shown in Sec. 3.2 as a quantitative test of the SCK model. Section 3.3 presents the values of D_fit and a_fit that have been obtained for a finite t_D without considering the limit. The errors of these values are not analyzed here but we believe that they are of the similar order as the errors of the “idealized” values. Therefore, taking into account that, the deviations between the SCK model and the simulations for the case presented in Sec. 3.3 are significantly larger than that analyzed in Sec. 3.2 we can treat the values of D_fit and a_fit presented in Sec. 3.3 also as quantitative ones but, their errors may be significantly larger then the errors of the “idealized” values.

3.2. The errors of the SCK model for long time phase of fluorescence quenching

All values of D_fit and a_fit presented in this subsection are the “idealized” values obtained by taking the limit for t_D → ∞. The results of the minimization at fixed a = a₀ are presented as D_fitversusa₀ in Figs. 7–9. We show the results for: Fig. 7, the simple liquids for various densities at a fixed interparticle potential; Fig. 8, the simple liquids for different interparticle potentials at a fixed density; Fig. 9–the mixtures for various temperatures, solute-solvent interactions, and R_C. The filled symbols always show the results for g_AB(r) ≡ 1. The empty symbols connected by a line show the results obtained by including g_AB(r) (evaluated from eqn. (2.14)) in eqn. (2.6). First, general conclusion from Figs. 7–9 agrees with the results of Zhou and Szabo.²² An inclusion of real g_AB(r) significantly improves an agreement between the simulation data and the SCK model.

Fig. 7 The “idealized” values. D_fit/D₀–1 at fixed a = a₀ as a function of a₀ for the L-J potential and various densities. Symbols: triangle up–#1 (ρ = 0.5787), square–#2 (ρ = 0.7830), triangle down–#5b (ρ = 0.9358). Filled symbols–g_AB(r) ≡ 1.0, empty symbols connected by a solid line–g_AB(r) taken from the simulations.

Fig. 8 The “idealized” values. D_fit/D₀–1 at fixed a = a₀ as a function of a₀ for various potentials at the same state point (ρ = 0.9358, k_BT = 1.0). Symbols: triangle up–#3 (HSS), square–#4 (SS), triangle down–#5b (LJ), circle–#6 (R6e). Remaining notation as in Fig. 7.

Fig. 9 The “idealized” values. D_fit/D₀–1 at fixed a = a₀ as a function of a₀ for mixtures at various temperatures (k_BT), cut-off distances (R_C) and solute-solvent interactions (ε_AC). Symbols: triangle up–#7, square–#8b, triangle down–#9, circle–#10. Remaining notation as in Fig. 7.

For g_AB(r) ≡ 1 the main conclusion is that for all the cases except of the largest diffusion (triangle up in Fig. 7, the case is discussed separately) the relative deviations of D_fit from D₀ (D_fit/D₀ − 1) can be described by similar curves with a minimum at a/σ_AB ≈ 1.4 for the simple liquids and a/σ_AB ≈ 1.2 for the mixtures. The minimum values are about:–0.21 and–0.25 respectively. Both the shape of D_fit/D₀ − 1 curves and the minimum values are similar in spite of significant differences in densities and interparticle interactions. Therefore, we believe that the presented dependence of D_fit/D₀ on a₀ is generally held for spherical molecules if D₀ is small. Some differences between the simple liquids and the mixtures occur in D_fit/D₀ − 1 at a₀ ≈ σ_AB. For the simple liquids (Fig. 8) the deviation at a₀ = σ_AB is much closer to 0 than the value at the minimum. For the mixtures (Fig. 9), the deviations at a₀ = σ_AB are strongly negative and comparable with the value at the minimum (see also Table 1). According to our considerations from Sec. 2.2.1, we can also predict that, if g_AB(r) ≡ 1 then the minimum values of deviations of D_fit from D₀ for the SFNL model are very close to the values from Figs. 7–9.

Fig. 7 shows that, the deviations in D_fit/D₀-1 for g_AB(r) ≡ 1 at high D₀ (#1 and, to some extent, #2) are qualitatively different form the discussed above. They are characterized not only by a minimum, but also by a maximum. For #1 the maximum is very sharp (+0.55 at a₀ = 1.063σ_AB). The effect is additionally strengthened because ∂k_SCK/∂D_S → 0 for D_S → ∞. The sharp maximum is a consequence of g_AB(r) ≡ 1 because for g_AB(r) taken from the simulation D_fit/D₀ − 1 at a₀ = 1.063σ_AB decreases to −.075 only.

The deep minimum which appears in D_fit/D₀ curves obtained for g_AB(r) ≡ 1 (black symbols) comes from the presence of the first coordination shell in the structure of real liquids. The shell is manifested by a sharp maximum in g_AB(r) at r slightly above σ_AB (e.g.Fig. 2). The screening by the first coordination shell decreases the probability of reaction between A* with B. Due to the finite size of particles the screening is effective at a range of distances larger than the width of the first peak of g_AB(r) (it can be seen as the first minimum of g_AB(r)). As a result, D_fit calculated with the assumption g_AB(r) ≡ 1 for values of a₀ in the range where screening is effective are significantly reduced if compared with the real ones. At very short distances the probability of collision with a screening particle decreases and the contribution to the transport process coming from ballistic motions becomes significant. As a result, the rate of approaches between A* and B increases. This explains why, for simple liquids, the values of D_fit/D₀ are significantly higher for a₀ = σ_AB than at the minimum (see Figs. 7, 8 and ΔD(1.0,1.13)_g≡1 from Table 1). For the mixtures (Fig. 9) the increase in D_fit/D₀ for a₀ approaching σ_AB is much less pronounced. We think that it is related to the presence of solvent. The C molecules, which are smaller than A (and B) ones, can efficiently screen B molecules at distances very close to σ_AB (which is impossible for the simple liquids) and so reduce the increase of D_fit at small a. It is confirmed by the shape of g_BC(r) function which has a sharp maximum for r ≈ 0.9σ_AB. The maximum in the local probability density at small intermolecular distances is also present if molecules are non-spherical. Therefore, the decrease of D_fit if a is close to the distance of the maximum is expected in analysis of experimental results. Such tendency (D_fit lower than the value of diffusion coefficient) has been recently observed.¹⁷ However, the values of D_fit analyzed in this experiment were not the “idealized” ones and, as it is shown in Section 3.3, they may be burdened with large errors due to non-diffusive processes neglected by the SCK model.

Figs. 7–9 show that, the result of including the real g_AB(r) in (2.6) depends on a₀. In the main region (a₀ < 1.5σ_AB), the inclusion significantly improves the agreement. The effect is especially strong for the potential R6e (the most realistic potential) for which D_fit/D₀ − 1 for real g_AB(r) is typically a few times lower than that for g_AB(r) ≡ 1 (compare ΔD from Table 1). Figs. 7–9 show also that, opposite to the case of g_AB(r) ≡ 1, the curves obtained for real g_AB(r) are similar only out of the main region of a. In the main region, the shapes of the curves visibly depend on the interparticle potential (Fig. 8) as well as on ρ (compare D_fit/D₀ for L-J potential in Figs. 7 and 8). All simulations for the mixtures were performed for the same density and for the fixed potential type, eqn. (2.10). As a result, the curves have similar shapes but the extreme values are different. Changes of ε_AC from 1.0 to 1.5 (#9 and #10 in Table 1), as well as the change of k_BT from 1.0 to 1.5 or R_C from 1.65 to 1.95, do not modify significantly the functional dependence of D_fit/D₀ on a. The values of D_fit/D₀ increase with the increase of liquid bonding (decrease of D₀). In spite of significant differences between parameters of intermolecular interactions we have not observed qualitative changes in the behavior of D_fit for the mixtures. We think it is related to a high temperature of the system. Simulations at distinctly lower temperatures were beyond our computational facilities. In such case t_T must be very large to attain equilibrium and perform credible measurements at a very low diffusion.

For a typical liquid system, the assumption g_AB(r) ≡ 1 is a very rough simplification so the model improves significantly when the real g_AB(r) is taken into account. The Smoluchowski approach treats reactants as large molecules moving in a Brownian medium and intermolecular interactions are incorporated only viag_AB(r). This is also an approximation because it treats the interactions in a statistically averaged way. Nevertheless, it describes the presence of local correlations more realistically than the assumption that they are absent (g_AB(r) ≡ 1). The observed deviations in D_fit/D₀ as well as the differences between these deviations observed for various systems are the consequence of higher order correlations between molecules. In order to describe them one may consider distance dependent diffusion, D_S(r) in eqn. (2.6). Some theoretical works on the distance dependent diffusion, D(r), were presented^31,32 but the problem has not been solved yet. The results presented in Figs. 6–9 suggest that, especially for strongly bonded liquids with a low diffusion, the inclusion of a correct form of g_AB function is more important than taking into account D_S is a function of r. Out of the main region the D_fit curves from Figs. 6–9 qualitatively agree with the form of D(r) proposed by Northrup and Hynes.³¹ However, Fig. 7 shows that the deviations of D_fit from D₀ depend on ρ, which was not considered in the mentioned paper. A more detailed analysis of the D(r) function on the basis of the presented results is impossible. D_fit(a₀) obtained here from the fit of k_SCK(t) for constant D_S may be treated only as a very rough estimation of D_S(a₀) because if D_S is a function of r, the formula for k(t) is slightly different²³ than eqn. (2.6) taken to the minimization.

All the a_fit/a₀ curves obtained from the minimization fixed D = D₀ and g_AB(r) ≡ 1 are qualitatively similar to the curves from Figs. 6–9. As an example, the “idealized” values of a_fit/a₀-1 and D_fit/D₀ − 1 as a function of a₀ for g_AB(r) ≡ 1 are presented in Fig. 10. The both curves are obtained from the minimization keeping the other parameter constant. The correlation between D_fit and a_fit is evident. The curves differ only by amplitudes, which for a_fit/a₀-1 typically amounts to 70–80% of that of D_fit/D₀ − 1.

Fig. 10 The “idealized” values. D_fit/D₀ − 1 at fixed a = a₀ (triangles down) and a_fit/a₀ − 1 at fixed D = D₀ (triangles up) for g_AB(r) ≡ 1 as a function of a₀ for #6.

The minimization with respect to both D and a have been performed only for g_AB(r) ≡ 1. As example, the “idealized” values of D_fit/D₀ and a_fit/a₀ as functions of a₀ for #5b and #9 are presented in Fig. 11. The functional dependencies are typical for the simulated liquids. A comparison of Fig. 11 with Figs. 8 and 9 shows that, for a low a₀, the minimization with respect to both D and a gives usually lower deviations than the fit for a single parameter. This is the result of the correlations between D_fit and a_fit (see Fig. 10). The correlation can be also noticed analyzing the form of eqn. (2.6). For example, for low D_S and t → ∞, k_SCK(t) ≈ 4πD_Sa and the both parameters are fully correlated.

Fig. 11 The “idealized” values. Simultaneous fitting. D_fit/D₀ − 1 (filled symbols) and a_fit/a₀ − 1 (empty symbols) for g_AB(r) ≡ 1 as a function of a₀ for two simulations: square, #5b, triangle down, #9.

3.3. The early phase of fluorescence quenching process

In order to compare the results presented in this subsection with a real experiment we assume the following values of the parameters of the A (≡ B) particles:

m_A = 200 amu, σ_AB = 8 × 10⁻⁸ cm, ε_AB = 300k_B (3.1)

The mass and σ-parameter presented above have been adopted to be similar to the parameters of the molecules form the experiments of Krystkowiak and Maciejewski.³³ In the following we present the dimensional values for the parameters (3.1) in square brackets.

Figs. 12 and 13 show k(t)/k(0) as a function of time for the early phase of fluorescence quenching for a₀ = σ_AB. The values of k(0) = k_SCK(0) have been evaluated from eqns. (2.5) and (2.8) taking g_AB(σ_AB) from the simulations. k(t)/k(0) obtained from the simulations (circles) is compared with the theoretical results obtained from the numerical solution of eqns. (2.4) and (2.5) with the real g_AB(r) taken from simulations (the solid line) and the formula (2.6) for g_AB(r) ≡ 1 (the dashed line). The theoretical curves have been calculated for D_S = 2D₀ and a = σ_AB. According to Table 1 the values of ΔD for #6 (Fig. 12) and #7a (Fig. 13) are small except of #7a for g_AB(r) ≡ 1. As a result, except of the case mentioned, the influence of the errors coming from the fact that D_S is a function of r (see Section 3.2) can be neglected. The differences between the solid curve and the circles come from the fact that the SCK model neglects non-diffusive motion. For the case shown on Fig. 12 (larger D₀) the curve obtained by the numerical solution matches the “experimental” points in a very wide range of times and the model fails just for t < t_B (in this case ≈ (m_A/ε_A)^1/2σ_AB [7 ps]). This break down is also seen in Fig. 13 at t = t_B close to that for #6 from Fig. 12. The break down can be explained considering the ballistic motions. At the very early stage of the reaction some A* and B molecules can come closer than a without colliding with another molecule (i.e.via ballistic motion). Such way of transport is much faster than diffusive one. As a result, for very short times, k(t) obtained from the simulations is significantly higher than that resulting from the model (which takes into account diffusion only). The effect lasts for a very short time because the time of the mean free flight for a dense liquid is very low (a few times shorter than (m_A/ε_A)^1/2σ_AB, which agrees with the value of t_B). On Fig. 13 the differences between the solid curve and the circles for t < t_B are larger than in Fig. 12 (the contribution from diffusive process decreases with D₀) and the inconsistency between the model and simulations spreads for the time of an order of magnitude longer than t_B. This inconsistency probably results from the higher order correlations between the reagents (for #7 the correlations are much stronger than for #6). However the effect is weaker than for t < t_B.

Fig. 12 k(t)/k(0) as a function of time for the early phase of fluorescence quenching for simulation run #6 at a₀ = σ_AB (circles) and the SCK model for D_S = 2D₀ and a = σ_AB. The solid line: the SCK model for g_AB(r) taken from simulations (from the numerical solution of eqns. (2.4) and (2.5)). The dashed line: the SCK model for g_AB(r) ≡ 1.

Fig. 13 k(t)/k(0) as a function of time for the early phase of fluorescence quenching of simulation run #7a at a₀ = σ_AB and the SCK model for D_S = 2D₀ and a = σ_AB. Notation as in Fig. 12.

The inconsistency between the solid curves and the dashed ones in Figs. 12 and 13 is not smaller than the one coming from the neglect of non-diffusive processes (the solid line and the simulation data) and, for both #6 and #7a, the effect is observed for times an order of magnitude larger than t_B. The inconsistency is an obvious consequence of the assumption that g_AB(r) ≡ 1. From the physical point of view the most important conclusion is that using the simplification g_AB(r) ≡ 1 it is not possible to describe the short time dependence of the “experimental” k(t) even in a qualitative way. This means that such simplification neglects some basic processes important for a description of the early stage of fluorescence quenching.

In the following we discuss the influence of the inconsistencies discussed above on the results of minimization. The inconsistency can be eliminated taking the limit t_D → ∞ (the “idealized” values) but in a real experiment this is frequently impossible. Our attention is focused especially on the case of t_D = 0 which means that all the obtained experimental data are treated as a whole, which is typical for analysis of experimental results. All the results of the minimization presented below have been obtained for g_AB(r) ≡ 1. In a real experiment, the influence of initial stage of the fluorescence quenching on results of the minimization strongly depends on the parameters describing the system and the experimental setup, potentials of intermolecular interactions, the de-excitation process studied, and other experimental conditions, which are beyond the considered system. Therefore, we have decided to give two examples only (Figs. 14, 15) and present conclusions of a general character.

Fig. 14 D _fit/D₀ − 1 at fixed a = a₀ = σ_AB as a function of t_D for various t_T for the mixture #8a. D₀ = 0.00616 [6.36 × 10⁻⁶ cm² s⁻¹]. Horizontal line–the limit of D_fit(t_D) attained for t_T ≥ 800 [2700 ps] (the critical value). Lines with symbols: circles–t_T = 25 [84 ps], triangles down–t_T = 100 [340 ps], squares–t_T = 400 [1350 ps].

Fig. 15 Simultaneous fitting. D_fit/D₀ − 1 (empty symbols) and a_fit/a₀–1 (filled symbols) as functions of t_D for various t_T for #5b. D₀ = 0.02797 [2.50 × 10⁻⁵ cm² s⁻¹]. Lines and symbols: triangles down, t_T = 100 [716 ps]; squares, t_T = 420 [3000 ps]; triangles up, t_T = 2184 [15 600 ps].

The dependence of D_fit/D₀ − 1 at fixed a for #8a on t_D for three different total times of measurement, t_T, are presented in Fig. 14. The deviations of D_fit from D₀ at t_D = 0 are extremely large. The effect is only due to inconsistency of the model with experiment, i.e. from the fact that the ballistic motion is neglected. The simulated system is always in the mechanical equilibrium state and reaction (2.1) does not disturb it anyway (see Section 2.1). The values of D_fit/D₀ − 1 strongly decrease with an increase of t_T, which is a consequence of the increased relative contribution coming from the diffusive processes. Analogically, a comparison between results of different simulation runs show that the deviation at t_D = 0 decreases with an increase of D₀ (compare also Figs. 12 and 13). We have also found that D_fit as a function of t_D strongly depends on a. For example, for the simulated mixtures, D_fit/D₀-1 for t_D = 0 changes its value from a large and positive at a₀/σ_AB = 1.0 (e.g.Fig. 14) to a large and negative at a₀/σ_AB = 1.27. In a real experiment a is never determined exactly (the molecules are not spherical ones) so, we are unable to estimate even a tendency and say if the effect overestimates or underestimates the resulting values of D_fit. The curve marked with triangles down in Fig. 14 has been obtained at the conditions similar to that of the experiment of Krystkowiak and Maciejewski.^17,18 The corresponding deviations at t_D = 0 amount to about 50% which is over two times higher than the deviation for the “idealized” values from Sec. 3.2.

The convergence in Fig. 14 seems to be very fast but as one can see the curves converge to the values different than the limit for t_D → ∞ and if t_T is small the difference may be quite large (see, especially the curve for t_T = 100). The effect occurs if t_T is lower than some critical value (as an example the value is given in the description of Figs. 14). A precise analysis shows that actually, for small t_T the value of D_fit systematically changes with t_D and the limit does not exist. But the rate of change is sometimes very low and in a real experiment the fact that the limit has not been attained yet can be missed especially because the analyzed curve loses its stability with t_D approaching t_T.

Fig. 15 presents the results of the minimization for both D and a for #5b. The shapes of the curves and the time scale are similar to that in Fig. 14 but the amplitudes are much higher in spite of D₀ being about four times larger than in #8a. The most surprising result is that for the highest t_T (=2184 [15 600 ps]) the deviations at t_D = 0 attain nearly 100%. It means that even if the data are collected for a very long time (for #5b, t_T = 2184 → (D₀t_T)/σ² ≈ 61) the contribution from the very early stage of quenching process is still important and may significantly influence the fit. This fact questions the credibility of results obtained when real experimental data are analyzed as a whole (i.e.t_D = 0). In the units resulting from (3.1) both D₀ (=[2.5 × 10⁻⁵ cm² s⁻¹]) and t_T (the minimal one = 100 [720 ps]) are significantly larger than those describing the system studied by Krystkowiak and Maciejewski.^17,18,33 For lower t_T and D₀ (as in their experiment) the deviations should be higher. Therefore, large errors discussed in their work¹⁸ may have the same origin as the inconsistencies between the SCK model and our simulations.

The main technical problem of fitting both parameters at once is not the convergence itself (as in the case on Fig. 14) but a lack of stability which is a consequence of the correlation between D and a (see also Fig. 10). If t_T is not sufficiently larger than t_D then the fitted values fluctuate with opposite phases and the amplitude of the fluctuations increases with increasing t_D (i.e.t_D approaching t_T). This effect appears for much lower t_D than the instability of the minimization with respect to one of the parameters.

4. Summary and final conclusions

In this work we have presented the results of molecular dynamics tests of the SCK model of diffusion controlled de-excitation process for spherical molecules. The simulations have been performed for various potentials describing interparticle interactions. We have found that the shape of potential has a little influence on the results thus we believe that the conclusions presented below can be generalized to other molecules of spherical or close to spherical shape. Our MD simulations were performed assuming that the reaction is instantaneous which agrees with the assumptions of the SCK model (the assumption (2.2)). As a result, we were able extract the errors coming only from the description based on the Smoluchowski equation. The agreement of the reaction radius a_fit and the diffusion coefficient D_fit obtained from the time dependent rate constant, k(t), with their values from simulations was considered as a qualitative measure of the model accuracy. In most cases studied here the inconsistencies appeared to be large and some times (Sec.3.3) the errors in D_fit and a_fit were of the same order of magnitude as the reference values. This shows that typical applications of the SCK model of fluorescence quenching process (i.e. with the assumption g_AB(r) ≡ 1) should be treated only as a rough simplification, even if the assumption on instantaneous reaction is justified.

According to results of the work a successful application of the SCK model to interpret experimental data (to extract a_fit and D_fit) mainly depends on possibility of discarding the data coming from the early stage of fluorescence quenching. If this condition is fulfilled (here, the “idealized” situation–Sec. 3.2) then the most important conclusions for the main region of a (σ_AB < a < 1.5σ_AB) are the following:

(i) If g_AB(r) ≡ 1 and diffusion is small then the values of D_fit/D₀ − 1 for a fixed a are almost always negative and can attain even 25%. This conclusion concerns also the SFNL model because it is asymptotically equivalent with the SCK one. For the case of a large diffusion (here, D₀ = 0.197σ(ε/m)^1/2) the maximum deviation exceeds +50%.

(ii) The deviations significantly decrease if we include the information on g_AB(r) into the model. However, the errors of about 10% may still occur. We expect similar improvement for the SNFL model. The inclusion of g_AB(r) is much more physical than the assumption g_AB(r) ≡ 1.

(iii) The relative deviations of a_fit at fixed D and g_AB(r) ≡ 1 are strongly correlated to that of D_fit at fixed a; they have the same sign and their amplitude is decreased by about 20–30%. As a consequence, the results of the minimization with respect of both D and a are usually closer to the real values than the results obtained from single parameter minimization at fixed a or D.

We expect that the long time behavoiur of k(t) does not depend on a particular form of the microscopic reaction cross section and its asymptotic form should be similar to an SCK model with “effective” k₀ and a. Therefore our conclusions for errors in the “idealized” values (in particular the negative deviation in D_fit if g_AB(r) ≡ 1) should apply to any reaction.

The results presented in Sec. 3.3 show that the SCK model fails completely for a very short time. If the distribution function of the reagents is incorporated into the model than in the considered examples (Figs. 12 and 13) significant differences in the k(t) curves are observed for times shorter than t_B ≈ (m_A/ε_A)^1/2σ_AB[7 ps]. If g_AB(r) ≡ 1 then the qualitative inconsistency between the model and simulations spreads out for times of an order of magnitude longer than t_B. Therefore, if one treats the obtained experimental data as a whole (i.e.t_D = 0) then the assumption g_AB(r) ≡ 1 does lead to large errors in the parameters obtained from the minimization. The errors are especially high if D is small and the time of measurement is short. The example form Sec. 3.3 shows that for t_T = 100 [340 ps] and D₀ = 0.00616 [6.36 × 10⁻⁶ cm² s⁻¹] the deviations for t_D = 0 attain 50%. If the minimization is performed with respect both D and a, without discarding results for short times, then very high deviations occur even if D₀ is not small and t_T is large (100% of error for (D₀t_T)/σ² ≈ 61).

In this work we estimate the scale of errors when the SCK model is applied to describe a reaction which agrees as much as possible with the assumptions of the model. The spherical shapes of molecules as well as the assumption on instantaneous de-excitation (2.2) are only rough simplifications. In real experiments on fluorescence quenching the molecules are never spherical ones and the de-excitation process itself lasts for a period of time. Therefore, the real errors, especially of the “idealized” values, should be higher than the estimated here, because many potential sources of errors are not taken into account. On the other hand, due to non-instantaneous de-excitation, the contribution coming from non-diffusive processes spreads out over time and, as a result, the dependence of the parameters on t_D may be less sharp than that for the SCK model (cf.Figs. 14, 15). Therefore, the influence of t_D on the parameters obtained for the real systems may be weaker than that presented in the work. This fact should not influence the two general conclusions important from practical point of view. First–the SCK model fails if t is very close to 0 but the resulting errors can be decreased (however not fully eliminated) if the shape of g_AB(r) is included into the model. Second–assuming g_AB(r) ≡ 1, the obtained results will be always only a rough estimation even if we fulfill all the conditions of the “idealized” situation.

The results of our work allows us to conclude that the shape of the intensity curve for very short times is the most important problem that influences the interpretation of experimental results on fluorescence quenching when the time of experiment is short. The problem can be solved in two possible ways. The first approach is to describe the quenching by a more general theory which includes the ballistic motion. The second approach is based on the fact that by analyzing the experimental data as a function of t_D (in order to eliminate the dominant part of the errors coming from the neglect of non-diffusive processes) and including g_AB we can decrease the errors to quite reasonable level. This approach requires sophisticated data analysis methods. We must also take into account that in a real system the g_AB function is not spherical. To conclude, the second approach is less efficient (e.g. the errors are still non-negligible) and less elegant than the first one but it seems to be easier for practical implementation.

Acknowledgements

The authors are grateful to Professor A. Maciejewski and Dr J. Kubicki for introduction to the experimental problems related to fluorescence quenching and discussions.

References

1. M. V. Smoluchowski, Z. Phys. Chem., 1917, 92, 129.
2. A. Szabo, J. Phys. Chem., 1989, 93, 6929.
3. F. C. Collins and G. E. Kimball, J. Colloid Sci., 1949, 4, 425.
4. R. M. Noyes, J. Chem. Phys., 1954, 22, 1349.
5. H. V. Beijeren, W. Dong and L. Bocquet, J. Chem. Phys., 2001, 114, 6265.
6. B. U. Felderhof and R. B. Jones, J. Chem. Phys., 1995, 103, 10 201.
7. B. U. Felderhof and R. B. Jones, J. Chem. Phys., 1997, 106, 967.
8. I. V. Gopich, A. A. Kipriyanov and A. B. Doktorov, J. Chem. Phys., 1999, 110, 10 888.
9. J. Sung and S. Lee, J. Chem. Phys., 1999, 111, 10 159.
10. T. L. Nemzek and W. R. Ware, J. Chem. Phys., 1975, 62, 477.
11. J. R. Lakowicz, M. L. Johnson, I. Gryczynski, N. Joshi and G. Laczko, J. Phys. Chem., 1987, 91, 3277.
12. R. Das and N. Periasamy, Chem. Phys., 1989, 136, 361.
13. A. D. Scully, S. Hirayama, D. Hachisu and T. Tominaga, J. Phys. Chem., 1992, 96, 7333.
14. A. D. Scully, S. Hirayama, K. Fukishima and T. Tominaga, J. Phys. Chem., 1993, 97, 10 524.
15. D. D. Eads, B. G. Dismer and G. R. Fleming, J. Chem. Phys., 1990, 93, 1136.
16. M. Sikorski, E. Krystkowiak and R. P. Steer, J. Photochem. Photobiol. A, 1998, 117, 1.
17. E. Krystkowiak and A. Maciejewski, J. Chem. Phys., 2002, 117, 2246.
18. E. Krystkowiak and A. Maciejewski, J. Chem. Phys., 2002, 117, 5802.
19. B. Cohen, D. Huppert and N. Agmon, J. Phys. Chem. A, 2001, 105, 7165.
20. H-X. Zhou and A. Szabo, J. Chem. Phys., 1990, 92, 3874.
21. W. Dong, F. Baros and J. C. Andre, J. Chem. Phys., 1989, 91, 4643.
22. H-X. Zhou and A. Szabo, J. Chem. Phys., 1991, 95, 5948.
23. J. B. Pedresen and P. Sibani, J. Chem. Phys., 1981, 75, 5368.
24. R. M. Noyes, Prog. React. Kinet., 1961, 1, 129.
25. D. A. McQuarrie, Statistical Mechanics, Harper and Row, New York, 1976.
26. M. P. Allen and D. J. Tildesley, Computer Simulations of Liquids, Oxford University, Oxford, 1987.
27. L. Verlet, Phys. Rev., 1967, 159, 98.
28. J. Gorecki, Mol. Phys. Reports, 1995, 10, 48.
29. B. Nowakowski, J. Gorecki and A. L. Kawczynski, J. Phys. Chem., 1998, 102, 7250.
30. J. S. Rowlinson and F. L. Swinton, Liquids and Liquid Mixtures, Butterworth, London, 1982.
31. S. H. Northrup and J. T. Hynes, J. Chem. Phys., 1979, 71, 871.
32. G. T. Evans and B. Kumar, J. Chem. Phys., 1988, 90, 1804.
33. A. Maciejewski, personal communication.

---