- Email: [email protected]
Non-orthogonal basis functions in the solution of diffusional problems
89
Chemical Physics 106 (1986) 89-101 forth-HoIland, Amsterdam
NON-OR~OGONAL BASIS FUNCTIONS ~IF~SIONAL PROBLEMS
IN THE SOLUTION
OF
Giorgio MORO Dipariimento
di Chimica Fisica, Universiib di Padova, 35131 Padua, Italy
Received 19 December 1985
The diffusion equation in presence of strong potential gradients can be conveniently solved by means of a suitable set of non-orthogonal basis functions that mimic the asymptotic behaviour of the eigensolutions of the time evolution operator. This method is discussed in the context of the calculation of spectroscopic observables by means of the Lanczos algorithm, that is by representing the spectral density as a continued fraction. The modified Lanczos algorithm, which takes into account the non-orthogonality of the basis functions, is analysed in detail. In order to test the performance of the new procedure, the planar rotor subject to a cosine potential is considered, and the set of non-orthogonal functions having the correct asymptotic behaviour is derived in both cases of single and double minimum. The results of the numerical calculations clearly show the advantages of the proposed method when dealing with systems characterized by hindered motions.
Nature provides a variety of situations where the molecular motion is strongly hindered by a local field. As examples one can consider conformational motions [l-4], the rotational and the translational motions of solutes in anisotropic liquids like nematics or smectics [5,6], or the reorientational motion of molecules in solids [7]. Often an adequate model for hindered motions is supplied by diffusion-type equations having a non-uniform stationary solution defined in terms of a potential function which takes into account the effects of the local field. Closely related are the studies of unimolecular kinetics in condensed phases where diffusion equations or FokkerPlanck equations are used with potential functions having narrow wells which identify the reactants and the products [8-121. In all these cases, it is necessary to work out the solutions of the diffusion equation with strong potentials, pertinent to the experimental observables. The available methods fall mainly in two categories: the numericrti imputation of the observables, or the analytical calculation of their asymptotic behaviour in the limit of infinite gradi-
ents in the potential function. Since the landmark contribution of Kramers [S], much work has been devoted to the latter aspect of the problem with a variety of methods [13-151, such as the calculation of the mean exit time (16-181. In contrast, comparatively less attention has been given to numerical methods (for some recent works, see refs. [19,20]). However, the two approaches are complementary. In fact the asymptotic analysis is essential in clarifying the peculiar features of the diffusion equation, like its relation with the kinetic description and with the random walk equation [21,22], due to the presence of steep potentials. On the other hand, the numerical solution is required if accurate theoretical estimates of experimentally observable quantities are necessary. A variety of numerical methods are available for solving diffusion-Ike equations. As an example, the brownian molecular dynamics is often used in the studies [1,23,24] of the conformational motion in polymers. In this case, one generates a representative ensemble of stochastic trajectories from the corresponding Langevin equation. Other methods are based on the calculation of the experimental observables from the matrix representation of the time evolution operator. In particu-
0301-0104/86/$03.50 0 Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
90
G. More / Non-orthogonal
ham Junctions in diffusional problems
lar, the Lanczos algorithm [25-271 has proven to be a valuable tool when dealing with spectroscopic or scattering-type measurements, where the spectral densities represent the observables of the physical system [28,29]. The advantages of the method have already been discussed together with its relation with the continued fraction representation encountered in statistical physics [30]. The standard implementation of the Lanczos algorithm uses an orthonormal set of basis functions that normally are the eigenfunctions of the diffusion operator in the absence of the potential. However, as the potential increases in strength, a larger basis set and a progressively heavier computation effort are needed for reproducing the eigensolutions of the operator. In order to illustrate the last consideration, let us examine the equilibrium distribution function which is the stationary eigensolution of the diffusion equation. In the presence of strong potentials, the equilibrium distribution is a sharply peaked function, since it is non-negligible only near the minima of the potential. If this eigensolution is expanded on a set of functions having a smooth dependence on the coordinates (say the circular functions for the planar rotor problem or the Wigner rotation matrices for the rotational diffusion of a molecule), a large number of terms is necessary to ensure convergence. With increasing potential strengths, the equilibrium distribution tends to resemble a delta function (or a sum of delta functions), and so the basis set must be infinitely expanded to represent accurately the stationary solution. The growth of the basis size with the strength of the potential constitutes the main handicap of the numerical solution of diffusional models in the study of hindered motions. It would be possible to overcome this difficulty, if one can use basis functions that mimic the known asymptotic solutions of the diffusion operator. In such a case, it is reasonable to expect that few basis functions are required for an accurate calculation of the spectroscopic observables. In a short communication [31], a set of functions having this feature, the so-called modified Wigner functions set, has been proposed for the specific problem of the reorientational motion of rigid molecules in nematic liquid
crystals. The asymptotic behaviour of the solutions of the diffusion equation was reproduced by the equilibrium distribution taken as common factor of the basis functions. Since these functions are not mutually orthogonal, an extension of the Lanczos algorithm to a non-orthonormal basis set was reported in the same work. A related approach has been considered by Blackmore and Shizgal [19]. They have solved the diffusional problem for a particular bistable potential on a polynomial basis, by taking the equilibrium probability as weight function. Their treatment, however, is based on the analytical orthogonalization of the polynomial basis. The orthogonalization procedure is not easily applicable to problems characterized by potential functions which are not expressed as simple polynomials, but rather by trigonometric functions. In general, one should use basis functions defined in terms of the equilibrium distribution, but avoid the unnecessary analytical orthogonalization by adopting the numerical procedure naturally provided by the modified Lanczos algorithm [29,31]. In some recent studies [21,22] concerning the transition processes amongst potential wells, another type of non-orthogonal functions, called localized functions, was considered. These functions represent the minimal basis set required to generate the random walk equation from a complete description of the diffusion motion, in the limit of infinite potential barriers. With finite potential barriers, the localized functions allow one to obtain a very good approximation of the low-frequency branch of the spectrum of the diffusion operator. On the other hand, the use of the localized functions for a “numerically exact” (i.e. within an arbitrary accuracy) calculation of the spectral densities was not considered in refs. [21,22]. The use of basis functions which mimic the asymptotic behaviour of the solutions of the diffusion equations constitutes a promising method in treating molecular hindered motions, since it should minimize the size of the matrices to be handled. One would like to have orthogonal bases in order to apply standard numerical techniques. But this in general is unlikely to occur, unless particular forms of the potential are chosen [19].
G. Moro / Non-orthogonal basis funciions in diffusional problems
Therefore, the calculation of the observables from non-orthogonal basis sets appears a necessary task. Actually, one can resort to the Gram-Schmidt orthogonalization procedure applied to the ensemble of functions having the correct asymptotic behaviour, but an alternative and preferable route is represented by a suitable modification of the Lanczos algorithm [31], because of its intrinsic efficiency in computing spectral densities [28,29]. We recall here that the Lanczos algorithm itself is an orthonormalization method applied to the moments of an operator acting on a given starting vector [32]. In this article, the use of non-orthogonal functions in the planar rotor problem will be analysed in detail. The choice of this particular problem has been dictated, apart from its relation to chain dynamics [2-41, by the simple form of the diffusion equation, in the belief that the advantages and the disadvantages of the method can thus manifest clearly. Of course the ultimate and future goal is the application of the method to complicated problems, like multidimensional diffusion or Fokker-Planck equations, where one cannot afford the standard implementation of the Lanczos algorithm in presence of a strong potential because of the extent of the basis set. In section 2, the modified Lanczos algorithm will be presented in the context of the calculation of observables given by spectral densities. The application of the method to single minimum potentials will be considered in section 3. It will be shown that a choice of basis functions closely related to those already presented in the context of rotational motion in nematics [31], is adequate to efficiently solve the problem. These functions, in fact, are suitable to represent small amplitude motions around the minimum of the potential. The diffusion equation with a double minimum potential will be treated in section 4. In this case, two independent physical processes are present: the small-angle displacements about the potential minima and the transitions between potential wells. It will be shown that when dealing with spectral densities influenced by the transition kinetics, we need to insert among the basis elements functions which are similar to the localized functions previously considered [21,22]. The per-
91
formances of the non-orthogonal basis functions will be discussed and compared with the results of the standard orthogonal representation of the diffusion operator.
2. The Lanczos
algorithm with non-orthogonal
basis functions
The spectral density enters as a fundamental quantity in the interpretation of spectroscopic [33,34] or scattering-type [35] measurements. In particular its frequency dependence determines the observable effects of the molecular motion. Therefore the calculation of spectral densities is an important tool in providing an “interface” between the experiments and the statistical physics which supplies the possible models of molecular motion. The Lanczos algorithm [25-271 represents a convenient method for calculating spectral densities since their frequency dependence is generated from the continued fraction representation [28-301, in analogy with the projection method proposed by Mori [36]. Moreover the intrinsic efficiency of the Lanczos algorithm with respect to the computer resources (i.e. computation time and size of the memory) is stressed when the spectral lineshapes are calculated instead of the eigenvalues [28,29]. The calculation of spectral densities by means of the Lanczos algorithm with an orthonormal basis has been already analyzed in the past [28]. Here we shall summarize its application to a physical problem described by the classical stochastic variable x = x(t). For a given function F(x) representing the observable probed in a spectroscopic experiment, the spectral density J(w), that is the Fourier-Laplace transform of the correlation function F( x (t ))*F( x(0)) [33,34], is written as J(w)=(FI(io+r)-‘JFP),
(I) where P(x) represents the equilibrium distribution function and the scalar product is defined as an integration on the variable x. In the relation (l), r denotes the time evolution operator for the non-equilibrium probability p ( x, t )
G. More / Non-orthogonal basis functions in diffusionalproblems
92
a+,
g/at=
-Q+,
(2)
t).
We consider only markovian and stationary processes, so the equilibrium distribution function is the unique stationary solution of r, rP=Q.
(3)
In the calculation of spectral densities, venient to use the symmetrized time operator defined in the relation
it is conevolution
i; = p-wrpw
(4)
There are different ways to implement the Lanczos procedure, depending on the symmetry of i; [29,37]. In the following we shall specialize the method to self-adjoint operators, like the diffusion operator with reflecting or periodic boundary conditions. For the spectral density given in eq. (l), an orthonormal set of functions c$~(x) is generated by means of the Lanczos method, starting from function @r(x) defined as +1 =
FP"2(
I’) -1’2,
1F
where the average 1F I* is calculated with the equilibrium distribution function P(x). The following three-term recursive relation is written [29,30]:
where CY,is the diagonal matrix element of f in the I#+ representation and &+i is the real coefficient derived from the normalization conditions applied to r&+ 1. The operator p has a tridiagonal matrix representation on the +k basis, therefore the spectral density can be written in terms of the following continued fraction [38]
iw+cu,
-
I P2
I2
io+cu,-
: -
I I% I 2 iw + (Ye-
(7) In the usual
implementation
R,,=
of the Lanczos
represet of
(8)
If’lfd.
(f,
By indicating with xk the column matrix defined by the expansion coefficients of C#Qin the f/ basis, (9)
G!i = C(4,fi i the recursive P n+lX,+l=
relation P
-
(6) becomes
%l)%
-
P,LI.
(10)
This is the standard Lanczos relation which can be easily implemented in a computer program for calculating the ensemble of coefficients C.X~ and flj [26-281. The frequency dependence of the spectral density is then derived from the continued fraction of eq. (7). In a recent communication [31], a modification of the Lanczos algorithm that permits the use of a non-orthogonal basis set, was proposed. We shall now analyse in more detail the implementation of this new method, with particular emphasis on the truncation of the, in principle, infinite set of basis functions. The truncation problem was not considered in the previous communication. Let f/ denote the non-orthogonal basis functions, and S the corresponding normalization matrix, s,j'=
(11)
lh’>.
It is implicitly assumed that the functions f/ are linearly independent and that they constitute a complete set for representing the eigensolutions of f. The functions r’m can be formally decomposed as
(12)
Ffm= Cn/r,mfn. m
1 J(4/lF12=
algorithm [28,29], one considers the matrix sentation R of r on a given orthonormal basis functions 6,
It is easily shown that eq. (10) continues to hold, with the column matrices x, defined according to eq. (9) if the matrix R is replaced by the matrix M, that is P n+1x,+1=
CM
It is convenient
(13)
-%1)%-&x,-I-
to introduce
the auxiliary
column
G. Moro / Non-orthogonal basis functions in diffusional problems
matrices
yk,
9N=
Yk=sq,
(14
such that the normalization condition on +k and the diagonal element (Ye can be calculated from the relations Y,$k
= 1,
CQ = yiMx,.
05) 06)
The sequence of operations of the standard Lanczos algorithm [26-281 can easily be adapted to the previous relations. In each iteration, one starts with the column matrices x~_~, xk and yk and proceeds as follows: (1) calculate simultaneously CQ and &+ ,xk + , , that is the right-hand side of eq. (13); (2) calculate from eq. (14) the column matrix &+l Yk+l; (3) calculate, by taking into account eq. (15), P k + , and normalize xk + 1 and yk + , . Note that the first step can be implemented like in the standard Lanczos algorithm by superimposing the column matrix Pk+l~k+l to x~_~. Moreover, it is clear from the sequence of operations, that only one array is necessary for storing column matrices y, and yk+i. Therefore, the storage is increased by one column matrix (yk) and one square matrix (S) with respect to the standard Lanczos algorithm. As regards the time-consuming operations, the number of multiplication of square matrices by column matrices is doubled, as effect of computing yk+, at each iteration. Of course the global efficiency of the algorithm, in both the standard and the modified form, is strongly dependent on the sparsity of the square matrices [26-281. The set of functions fj should be infinite if it must satisfy the condition of completeness with respect to the ensemble of the solutions of the diffusion operator. Of course, in computer applications, only finite matrices can be handled, in correspondence of a subspace cN spanned by functions fi for j = 1,2,. . . , N. Therefore the operator f in eq. (6) should be substituted with its projected form
E Ifi
)(SN1)j,'(f/'
93
I
9
(18)
j,j'=l
where S, is the first N X N block of the normalization matrix. Also the function FP*/*, which determines C#Qaccording to eq. (5), should be projected onto eN, but usually it is convenient to include the function FP”* in the finite basis set. In eq. (13) which implements the Lanczos procedure, the infinite matrix M is then replaced by the N x N matrix M N given implicitly by the equation
(19) One easily recognizes that the elements of M, are in general different from the corresponding elements of M, therefore M, cannot be generated simply by truncating M. By taking into account eq. (18) for the projection operator, MN can be explicitly written as M, = S,‘R,,
(20)
where R, is the first N X N block of the matrix defined in eq. (8). The use of eq. (19) for calculating the elements of M, is convenient when p’n can be analytically written as a linear combination of basis functions (see the next section for an example). Otherwise M,,, must be derived from eq. (20). Jones and co-workers in their calculations of the electron charge density in solids [39,40], have implemented the Lanczos algorithm according to the second procedure. In particular they avoided the timeconsuming inversion of S, by computing, at each iteration of eq. (13), the array S,‘R,x,, from the solution of the linear equation system with coefficients given by S,.
3. The single minimum planar rotor The time evolution operator r for the diffusional motion of the variable y describing the rotation of a dumbbell in a plane, is written as r=
-(a/ay)PD(a/ay)P-‘,
(21)
G. Mom / Non-orthogonal
94
basis functions
where D is the diffusion coefficient and P(y) is the equilibrium distribution function that can be defined in terms of the potential V(y) acting on the dumbbell, P(y)
= exp( - V/k,T)/L2”exp(
- V/k,T)
This type of equation is encountered in the theory of internal motion of molecular conformers [3]. We do not consider here the angular dependence of D; however, once a model for D(y) has been chosen, the following analysis can be easily adapted to the specific case. In this section we shall analyse the case of a symmetric potential (23)
having a single minimum centered at y = 0 and with V(0) = 0. The asymptotic expansion of both static (mean values) and dynamic observables (correlation functions or spectral densities) are obtained by introducing a fictitious parameter h (taken as unity at the end) as factor of the potential I/ and considering the limit for X -+ 00. The value of the observables in the limit of highly anisotropic potential is then derived from the leading terms in the expansion in powers of X-‘12. The symmetrized diffusion operator, in the approximation of zeroth order, becomes F = zJ(x2-
I -
a2/ax2),
(24)
where y = { XV”(0)/(2k,T)}
-“2x,
u = DXV”(0)/(2k,T),
(254 (25b)
and V”(0) indicating the second derivative of the potential calculated in the location of the minimum. The corresponding asymptotic eigenfunctions q,* are written in terms of the hermitean functions H,(x), \k,”= ( AV”(0)/2nk,T)“4 \k: = ( j! 2’) -I’*H,(x)
exp( -x2/2), q;(x).
problems
potential minimum. The asymptotic behaviour of the spectral density given in eq. (l), is then derived from the following representation of FP”‘: (27)
dy. (22)
V(Y) = v(-Y)
in diffusionul
(26a) (26b)
Such approximate solutions of the diffusion operator represent the small-angle motions around the
where coefficients cj are calculated from the Taylor expansion of F(y) around y = 0. Since the coefficient c, is of the order of h-‘J-“/2, the functions !$A contribute to the spectral density with decreasing weights given by powers of l/X. In other words, by considering the approximated diffusion operator of eq. (24) in the subspace e”, spanned by functions $, !I$, . . . , !Pi-‘, that is in the ensemble of polynomials of (N - 1)th degree in the variable x multiplied by @, one takes into account in the spectral density the first N contributions with respect to the powers of the asymptotic parameters 1 /A. Although the functions q; reproduce well the behaviour of the solutions of the diffusion operator around y = 0 with steep potentials, they are not useful in numerical calculation of spectral densities. In fact they cannot be used for generating a matrix representation of the complete diffusion operator of eq. (21), because the periodic boundary condition is not ensured. In ref. [31], a set of functions that mimic the asymptotic solutions, still preserving boundary conditions, has been proposed for the rotational diffusion of an axially symmetric molecule in nematic liquid crystals. The same method can be applied to the simpler diffusion problem considered here, in this way enlightening its rationale. The basic idea is to take advantage of the knowledge of the stationary solution of r, that is the subspace CL previP(Y) ‘I* 3 for approximating ously defined, with a subspace zN made by functions which satisfy periodic boundary conditions. in all generThe subspace cN can be constructed ality from the linear combination of functions f, = g, P’/~, fi = g, P112,. . . , fN = gNP112, where the g,(y) are a set of functions having a well defined Taylor expansion around y = 0. Of course subspaced ek and eN are not strictly comparable, because they refer to functions having different boundary conditions. Nevertheless, they are equivalent in the asymptotic limit. In fact, by making
G. More / Non-orthogonal
basis functions in diffurional problems
the substitution eq. (25a) and expanding functions up to order (N - l), the fi as powers of X-l/’ ensemble of functions pN_I(x)!Pt(x), where pN_ 1( x) indicates a generic polynomial of (N 1)th degree in the parameter x, is recovered. Therefore, the leading contributions with respect to the asymptotic parameter l/h are taken into account, and the ensemble of functions fJ previously defined represents a convenient basis set in the calculation of the spectral densities. We emphasize that a precise correspondence between each function fi and the asymptotic eigenfunction @,” is not necessary, since our objective is to provide an optimal subspace to the numerical solution of the diffusion operator by means of the Lanczos algorithm. Apart from the condition of linear independence, we have left much freedom in choosing the functions g,. The eigenfunctions of the isotropic diffusion operator, that is the sine and cosine functions, are a natural choice in order to represent conveniently the behaviour of the operator also in the limit of weak potentials. The resulting functions fJ are non-orthogonal. This is the price that must be paid in order to mimic the asymptotic behaviour of the solutions of the diffusion operator. On the other hand, the modified Lanczos algorithm discussed in the previous section, can handle such type of basis functions. In the remainder of this section, we shall discuss in detail the application of the method to the calculation of the spectral density of the function F(y) = cos y - cos y, with a cosine-type potential
composed
as follows,
=_i’Df-
p&
(A/2)(1
- cos y).
The following set of non-orthonormal tions has been used: f, = (cos( jy)
- cos( jy
)) P”2.
(jDA/4)(f,+i
(30)
(31) It is known that the Lanczos algorithm can reproduce correctly the spectral density with a number of iterations less than the dimension N of the matrix [28,29]. Also in the present calculations we have found the same behaviour, but in the follow-
Log
E
10 N
0 0
(28)
0 OO
basis func-
0 -2--
0
n
0
(29)
Evans, in a recent analysis of the isomerization of butane [4], has considered a basis set related to the functions of eq. (29). However, he calculates the matrix representation of the time evolution operator with the corresponding orthonormal set generated by means of the Gram-Schmidt orthonormalization procedure. Since the function l;fi can be analytically de-
-&-i),
the elements of M, are conveniently computed by means of eq. (19). In particular, the projection should be taken into account exoperator 9, plicitly only for j = N, where the solution of a linear system of equations having (fN+, I&) as constant terms, is required according to eq. (18). For characterizing the efficiency of the method, a criterion is needed for the error EN of the approximated spectral density JN (w), calculated with a N-dimensional basis set, with respect to the exact J(w). Like in ref. [28], we use the normalized difference in the area of the real part of the spectral density,
0 --
V/k,T=
95
-4--
0
n
0
.
5
10 N
Fig. 1. Relative error in the area’ EN (logarithmic scale) as function of the size N of the basis set, for the spectral density relative to cos y - cos y with single minimum potential eq. (28) and A = 20. (0) Orthogonal basis functions; (A) nonorthogonal basis functions.
96
G. More / Non-orthogonal basis functrons in diffusional problems
ing discussion we disregard this aspect by letting n = N. In fig. 1 the error E, is plotted as function of N for both the non-orthogonal and orthogonal basis set. In the second case the standard Lanczos algorithm has been used with a basis set of normalized cos(jy) functions. The value of A used (A = 20) corresponds to a very steep potential. A large number of orthogonal basis functions is necessary for obtaining a reasonable accuracy in the spectral density, because many cosine functions are needed to represent FP’/‘, as a consequence of the factor P ‘I2 characterized by a sharp peak centered at y = 0. A comparable accuracy in the results is achieved with one or two non-orthogonal basis functions only. The growth of the size of the orthogonal basis set with increasing A, is manifest from fig. 2 where the least N which ensures an error EN G 10.e3, is plotted against A. On the contrary, the number of non-orthonormal basis functions decreases with A, starting from A = 9. In fact one can easily show that function f, alone predicts the leading contribution to the spectral density in the limit of an infinite A. The relative advantages of the two methods can be understood from the trends shown in fig. 2, by taking into account that the modified Lanczos algorithm discussed in the previous section, is
always less efficient than the standard Lanczos algorithm in treating matrices of the same dimension, as a consequence of algebraic operations involving the normalization matrix S,. Therefore the use of non-orthogonal functions becomes convenient only when it allows a reduction of the basis set. It is clear from the behaviour displayed in fig. 2, that this condition should be always reached with a sufficiently large value of A. Of course stronger potentials lead to even greater savings in computation time. On the other hand a general criterion for the relative convenience of the non-orthogonal basis set with respect to the, normally used, orthogonal ones cannot be given, for it depends on the sparsity of matrices M, and S, appropriate to a particular diffusional problem. Finally, we mention that also in the case of odd functions of y, trends like those shown in figs. 1 and 2 are obtained with non-orthonormal functions f, = sin( jy) PI/‘. 4. The double minimum planar rotator In this method of applied to specifically potential V/k,T=
8
0
0
n
n
0 N 0 4
0 0
n n
n
n LJ 0 :;;;;; 0
-10
1 A
20
Fig. 2. Least-size N of tie basis set that assures a 10m3 accuracy in the calculation of the special density relative to cost -COS Y, as function of the parameter A of the potential given in eq. (28). (0) Orthogonal basis functions; (A) nonorthogonal basis functions.
section we intend to show how the non-orthogonal basis functions can be multiminima potentials, by considering the planar rotor with the two-minima
(A/2)(1
- cos 2y) = A sin2y,
(32)
with a positive A. Since the diffusion operator is invariant under symmetry operations of the point group CzV, the calculation of the spectral densities can be factorized in the four functional subspaces pertinent to the irreducible representations. The asymptotic analysis of the small amplitude motions can be applied independently to each potential minimum and, in strict analogy with the treatment illustrated in section 3, suitable sets of nonorthogonal basis function are generated. These functions can be factorized in the four subspaces as follows: A, :
f, = [ cos(2jy)
A,:
f, = sin(2 jy) P”2,
- cos(2jy)]
P1’2,
G. More / Non-orihogonal basis functions in diffusional problems
B, :
f, = cos[(2 j - l)y]
PI”,
B,:
f, = sin[(2 j - l)y]
PI”,
where the index j runs over positive integer numbers. We do not report here specific results for functions belonging to representations A,, A, and B,, since the same trends as discussed in the previous section are recovered and thereby identical conclusions can be drawn. A completely different behaviour is found instead with functions of B, symmetry. In fig. 3 the error EN on the spectral density relative to function F = cos y is displayed for some values of the size N of the basis set in both cases of orthogonal (circles) and nonorthogonal (triangles) functions. The ensemble of functions 7~~~‘~ cos(2j - l)y, with j = l-N, was taken as orthogonal basis set. A large value of A (A = 12) was used in order to report the behaviour for highly anisotropic systems. A striking difference with respect to fig. 1 is manifest: only a little improvement on the calculation is achieved when passing from an orthogonal to a non-orthogonal basis set with the same dimension. It has been also found that the size of both types of basis sets must be increased with larger values of A, in order to calculate spectral densities with the same accuracy. A similar behaviour is reported in the Log
0--
N
A
A
O
A
0
A
0 A A
-2 --
work of Blackmore and Shizgal about the numerical solution of the diffusion equation of a bistable potential (compare, in particular, tables II and III of ref. [19]). In fact, as mentioned in section 1, an analogy exists between our non-orthogonal functions fi and their orthogonal polynomials constructed with the equilibrium distribution as weight function. The inefficiency of the non-orthogonal functions of B, symmetry can be explained by taking into account the peculiar features of the bistable system. There is, in fact, an eigensolution ?P, of the diffusion operator having the eigenvalue (Y, which vanishes exponentially like an Arrhenius factor exp( - E/k,T), with an increasing energy barrier E between the minima of the potential [13,14]. The solution q, is representative of the slow process of transition between the potential wells. The decay behaviour of the correlation function associated to cos y is dominated by the root (Y, which, instead, does not contribute for symmetry reasons to the correlation functions relative to the other representations A,, A, and B, [21]. In fact *I belongs to the B, representation because it is the antisymmetric counterpart of the stationary solution ‘k, = P”‘. It can be more conveniently written as *, = g(u>P”‘,
(33)
where, in the limit of high-energy barriers, g(y) has a step-like behaviour, that is g - 1 (g = - 1) for A in a wide interval centered at y = 0 (y = n) and it changes abruptly for y = +~/2. In analogy with the bistable system defined in an infinite domain [13], the asymptotic analysis leads to a g(y) with an error function shape,
E I0
97
0
g(y)r=
-erf[A1/2(y-7r/2)],
(34)
A q A -4
0
-rJ 0
I
I
A
I
5
10 N
Fig. 3. function relative A = 12. functions
Relative error in the area EN (logarithmic scale) as of the size N of the basis set for the spectral density to cos y. with the bistable potential eq. (32) and (0) Orthogonal basis functions; (A) non-orthogonal j,; (0) basis set with f, and functions h,.
for 0 6 y i 71. It is evident that the function f, = *i only in a (cos y ) P”2 mimics eigensolution small interval around y = 0 (y = IT) where cos y = 1 (cos y = -l), but it is unable to predict the error function behaviour near potential maxima. Remarkable improvements are not expected from the other functions 6.. The inefficiency of the non-orthogonal functions A in the calculation of spectral densities relative to the function cos y should be ascribed
to their inadequacy in representing *i. Therefore, one needs to include, in the basis set, functions which approximate the eigensolution !P,, that is functions describing the transition between the potential wells. In some recent works 121,221, a class of such functions, called localized functions, has been introduced. The localized functions are constructed by matching in the locations of the potential minima, the functions g,(y)
= j%(+‘da, (35) ?J which are locally exact solutions of the differential equation PgP ‘P = 0.
(36)
One easily recognizes that the functions g, of eq. (35) behave like the error function of eq. (34) in correspondence with the maxima of the potential. Moreover, the ensemble of localized functions have been identified in refs. [21,22], with the least set of functions which generates the random walk equation amongst the stable states (i.e. states localized at the potential wells) from the diffusion equation. In the present work we prefer, instead, to work with a generalized form of localized functions which have continuous derivatives. They are obtained from the approximation of function g( y > g(y) = ri2G(cr)
P(cu)-‘dru,
(37)
with function G(a) obeying to the following conditions: (i) G(a) must belong to symmet~ subspace %* (ii) G(a) must behave like a non-vanishing constant near the potential maxima. The former ensures that g(y) belongs to the B, irreducible representation while the latter leads, in the asymptotic limit, to the correct error function behaviour of eq. (34). All the functions sin(2j 1)a satisfy these conditions, and one can insert the following functions hi, hi(y)=P’/‘III/‘sin(2jY
l>a P(a)-‘da,
among the elements of the basis set. h, constitutes the best trial function
(38)
for @,,
because the eigensolution cos y with the lowest non-vanishing eigenvalue is recovered for A = 0. Since both limits of small and large values of A are reproduced, one expects that also in the middle range of potential strengths h, approximates well the behaviour of !P,. As a test, we have calculated the approximate eigenvalue (Y~from the expectation value of i; with function h,, ai = (hi I f’l W/‘(h, 14) = D((sin2y}~(y)-‘)/(h,
I h,),
(39)
and compared it with the true eigenvalue obtained from numerical diagonalization of I=. In fig. 4 the results are displayed as a function of A, together with the leading term of the asymptotic expansion of ai, CY,= (4DA/a)
e-‘.
(40)
The numerical results for some selected values of A are reported in table 1. It can be clearly seen that relation (39) yields a good estimate, within an error of I%, of the eigenvalue (Y, in the whole range of potential strengths. h, can now be used, together with the functions 4, for calculating the spectral density by means of the modified Lanczos algorithm. Of course f, must always be included in the basis set, since it represents the starting vector in the calculation of the spectral density relative to function cos y. Since I
I
I
0
I
b.
1
2
3
Fig. 4. Lowest positive root IX, of the diffusion operator with the bistable potential given in eq. (32), as function of the parameter A (diffusion coefficient D =I). (0) Exact values: approximated values calculated from eq. (39); (- - -) c-) asymptotic relation eq. (40).
G. Moro / Non-orthogonal basisfunctionsin dijfusional problems Table 1 Lowest positive root of the diffusion operator with the bistable potential of eq. (32), as function of A (diffusion coefficient D =l) A
Exact CX,
Approximate
0 1 2 3 4 5 10 20
l.OtM 0.5876 0.3237 0.1676 0.8198 x 0.3823 x 0.5463 x 0.5112~
1.000 0.5890 0.3257 0.1689 0.8264x 10-l 0.3850 x10-’ 0.5473 x10-3 0.5113 x 10-r
10-l 10-t lo- 3 lo-’
(Y, ‘) (0.0) (0.4684) (0.3446) (0.1902) (0.9328x 10-l) (0.4290 x 10 - ’ ) (0.5781 x 10-3) (0.5249x10-‘)
a) (Y, calculated by numerical integration of eq. (39). asymptotic values obtained from eq. (40) are reported tween parentheses.
The be-
one cannot extend, in a simple way, eq. (30) to the function h,, the matrix M, is calculated from eq. (20). The efficiency of this type of non-orthogonal basis functions appear clear from the results shown in fig. 3. The two-dimensional basis set constituted by functions h, and fi, gives a spectral density as accurate as that obtained with seven functions fi alone. The addition of other functions hj to the basis set improves considerably the accuracy of the calculated spectral density. In fig. 3 only the result for the basis set formed by ft, h, and h, is reported, because wider sets obtained by adding more hj functions lead to a value of EN below the reported scale. We mention also that the efficiency of the basis sets made with functions h,, becomes more pronounced when the parameter A is increased. Therefore, all considerations of the previous section about the convenience of the modified Lanczos algorithm with respect to its standard form, can be extended to the bistable diffusion equation. The results reported in table 1 prove the validity of our choice of basis functions. In fact, it shows that a single function is sufficient to calculate with reasonable accuracy the lowest positive root. We remind that this eigenvalue determines the long-time behaviour of the bistable system and it is identified with the observed rate in model system for the chemical kinetics [8,9,21]. One could compare this result with that obtained by Black-
99
more and Shizgal [19] in their study of a related bistable system in one dimension. From tables II and III of ref. [19], it appears that 10 and 50 orthogonal polynomials are necessary to compute the same root with a comparable accuracy, for barrier heights of 2.5 and 25 k,T units respectively. One can conclude that a reduction of orders of magnitude in the size of the basis set can be achieved. by properly taking into account the asymptotic behaviour of the diffusion equation.
5. Conclusion The results presented in this work show that, in regimes of strongly anisotropic potentials characterized by large barriers, a proper choice of the non-orthogonal basis functions drastically reduces the subspace in which the diffusion operator must be solved numerically. In truth, this conclusion is supported so far by applications to simple diffusional problems, like the planar rotator considered here or the rotational motion in nematic liquid crystals analyzed in ref. [3]. On the other hand, our objective was to develop the methodology, leaving to future works its application to problems characterized by large-scale computations. However, we expect that also in complex systems, like those described by multidimensional diffusion equations, it should be possible to construct from the asymptotic analysis of the time evolution operator [18], a proper set of functions that mimic the fundamental processes determining the spectroscopic observables. The examples considered in this work tell us that the non-orthogonal functions must be chosen careffilly if the optimal basis set is required. In particular, the use of the equilibrium distribution probability as weight function, is sufficient only when dealing with the small amplitude motions around the minima of the potential. If a given spectral density is influenced by the transition process amongst potential wells, then one should include in the basis set specific functions which generate the random walk equation from the more general diffusion equation [21,22]. The proposed modification of the Lanczos algorithm allows the calculation of the spectral den-
100
G. More / Non-orthogonal
ham Junctrons
sities by means of non-orthogonal functions. One should realize that the computational procedure is more complicated than the standard Lanczos algorithm. Therefore the use of non-orthogonal functions becomes convenient only when it leads to a reduction on the size of the matrices to be handled. The numerical results previously discussed show that such a situation is always reached by increasing the energy barrier of the potential. Finally, we would like to suggest two fields of future applications of the non-orthogonal functions method with the modified Lanczos algorithm: (i) the Fokker-Planck equation in the low friction regime where very large basis sets are needed (see refs. [41,42] for an analysis of the asymptotic behaviour of the Fokker-Planck operator); (ii) the multidimensional diffusion equation pertinent to chain dynamics [l], where it is immediately evident the analogy between the planar rotor problem and the motion around each internal degree of freedom of the molecule. The chain dynamics is a manifest example supporting the need for efficient methods for representing diffusion operators. Let us consider the internal motion of the butane molecule under the action of the Ryckaert-Bellemans potential [22]. We have found that a minimal basis of 20 elements is necessary when the standard circular functions are used. Therefore, we expect that matrices of the order of 20 N should be solved for alkyl chains having N internal degrees of freedom. In this way, problems with more than three degrees of freedom could not be handled even with large computers. It follows that only with an accurate choice of the basis functions, which minimizes the representation for a single degree of freedom, it will be possible to tackle the complete analysis of the dynamics of large molecules. In conclusion, it is conceivable that, by properly taking into account the asymptotic solutions of the diffusion operator, very few basis functions might be strictly necessary to compute the dynamical observables of a physical system.
n diffusiond
problem
Acknowledgement This work has been supported by the Italian Ministry of Public Education and in part by the National Research Council through its Centro Studi sugli Stati Molecolari Radicalici ed Eccitati. The author is indebted to Professor P.L. Nordio for stimulating discussions and for advice on the revision of the manuscript.
References PI M. Fixman, J. Chem. Phys. 69 (1978) 1527, 1538. 121J. Skolnick and E. Helfand, J. Chem. Phys. 72 (1980) 5489. [31 D.C. Knauss and G.T. Evans, J. Chem. Phys. 72 (1980) 1499. 141 G.T. Evans, J. Chem. Phys. 78 (1983) 4963. 151 G. Moro, U. Segre and P.L. Nordio, in: Nuclear magnetic resonance of liquid crystals, ed. J.W. Emsley (Reidel, Dordrecht, 1985). [61 G. Moro and P.L. Nordio, J. Phys. Chem. 89 (1985) 997. rotations in molecular crystals [71 W. Press, Single-particle (Springer, Berlin, 1981). 181 H.A. Kramers, Physica 7 (1940) 284. 191 D. Chandler, J. Chem. Phys. 68 (1978) 2959. HOI J.L. Skinner and P.G. Wolynes, J. Chem. Phys. 69 (1978) 2143. [JLI B. Carmeli and A. Nitzan, J. Chem. Phys. 79 (1983) 393; 80 (1984) 3596. 1121 R.F. Grote and J.T. Hynes, J. Chem. Phys. 77 (1982) 3736. 1131 R.S. Larson and M.D. Kostin, J. Chem. Phys. 69 (1978) 4821. 1141 H. Tomita, A. Ito and H. Kidachi, Progr. Theoret. Phys. 56 (1976) 786. I151 B. Caroli, C. Caroli and B. Roulet, J. Stat. Phys. 21 (1979) 415. [161 K. Lindeberg and V. Seshadri, J. Chem. Phys. 71 (1979) 4075. [171 J.M. Deutch, J. Chem. Phys. 73 (1980) 4700. of stochastic differen[181 Z. Schuss, Theory and applications tial equations (Wiley, New York, 1980). 1191 R. Blackmore and B. Shizgal, Phys. Rev. A31 (1985) 1855. [201 W. Nadler and K. Schulten, J. Chem. Phys. 82 (1985) 151. 1211 G. Moro and P.L. Nordio, Mol. Phys. 56 (1985) 255. I221 G. Moro and P.L. Nordio, Mol. Phys. 57 (1986) 947. Chem. 1231 R.M. Levy, M. Karplus and J.A. McCammon, Phys. Letters 65 (1979) 4. 1241 T.A. Weber and E. Helfand, J. Phys. Chem. 87 (1983) 2881. [25] C. Lanczos, J. Res. Nat]. Bur. Std. 45 (1950) 255; 49 (1952) 33.
G. More / No-o~~~ogonai basis functions in ~iff~ionaI problems
[26] B.N. Parlett, The symmetric eigenvalue problem (Prentice-Hall, Englewood Cliffs, 1980). (271 J. Cullum and R. Willoughby, Lanczos algorithm for large symmetric eigenvalue computations, Vol. 1. Theory (Birkhauser, Bascl, 1985). (281 G. Moro and 3.H. Freed, J. Chem. Phys. 74 (f981) 3757. [29] G. Moro and J.H. Freed, The Lanczos algorithm in molecular dynamics: calculation of spectral densities, Proceedings of the Workshop in Large Eigenvalue Problems, Oberlech, Austria (1985), eds. J. Cullum and R.A. Willoughby. 1301 G. Moro and J.H. Freed, J. Chem. Phys. 75 (1981) 3157. 1311 G. Moro and P.L. Nordio. Chem. Phys. Letters 96 (1983) 192. [32] Yu.V. Vorobiev, Method of moments in applied mathematics (Gordon and Breach, New York, 1965).
101
[33] R. Kubo, J. Phys. Sot. Japan 12 (1957) 570. [34] B.J. Beme, in: Physical chemistry, Vol. 8B, eds. H. Eyring, D. Henderson and W. Jost (Academic Press, New York, 1971) p. 540. [35] L. van Hove, Phys. Rev. 95 (1954) 249. [36] H. Mori, Progr. Theoret. Phys. 34 (1965) 399. [37] W.A. Wassam, 3. Chem. Phys. 82 (1985) 3371, 3386. [38] H.S. Wall, Analytic theory of continued fractions (Van Nostrand, Princeton, 1948). [39] R. Jones and M.W. Lewis, Phil. Mag. 49B (1984) 95. [40] R. Jones, in: The recursion method and its application, eds. D.G. Pettifor and D.L. Weaire (Springer, Berlin, 1985) p. 132. [41] H.D. Vollmer and H. Risken, Physica 1lOA (1982) 106. [42] G. Moro, Chem. Phys. 96 (1985) 287.
Chemical Physics 106 (1986) 89-101 forth-HoIland, Amsterdam
NON-OR~OGONAL BASIS FUNCTIONS ~IF~SIONAL PROBLEMS
IN THE SOLUTION
OF
Giorgio MORO Dipariimento
di Chimica Fisica, Universiib di Padova, 35131 Padua, Italy
Received 19 December 1985
The diffusion equation in presence of strong potential gradients can be conveniently solved by means of a suitable set of non-orthogonal basis functions that mimic the asymptotic behaviour of the eigensolutions of the time evolution operator. This method is discussed in the context of the calculation of spectroscopic observables by means of the Lanczos algorithm, that is by representing the spectral density as a continued fraction. The modified Lanczos algorithm, which takes into account the non-orthogonality of the basis functions, is analysed in detail. In order to test the performance of the new procedure, the planar rotor subject to a cosine potential is considered, and the set of non-orthogonal functions having the correct asymptotic behaviour is derived in both cases of single and double minimum. The results of the numerical calculations clearly show the advantages of the proposed method when dealing with systems characterized by hindered motions.
Nature provides a variety of situations where the molecular motion is strongly hindered by a local field. As examples one can consider conformational motions [l-4], the rotational and the translational motions of solutes in anisotropic liquids like nematics or smectics [5,6], or the reorientational motion of molecules in solids [7]. Often an adequate model for hindered motions is supplied by diffusion-type equations having a non-uniform stationary solution defined in terms of a potential function which takes into account the effects of the local field. Closely related are the studies of unimolecular kinetics in condensed phases where diffusion equations or FokkerPlanck equations are used with potential functions having narrow wells which identify the reactants and the products [8-121. In all these cases, it is necessary to work out the solutions of the diffusion equation with strong potentials, pertinent to the experimental observables. The available methods fall mainly in two categories: the numericrti imputation of the observables, or the analytical calculation of their asymptotic behaviour in the limit of infinite gradi-
ents in the potential function. Since the landmark contribution of Kramers [S], much work has been devoted to the latter aspect of the problem with a variety of methods [13-151, such as the calculation of the mean exit time (16-181. In contrast, comparatively less attention has been given to numerical methods (for some recent works, see refs. [19,20]). However, the two approaches are complementary. In fact the asymptotic analysis is essential in clarifying the peculiar features of the diffusion equation, like its relation with the kinetic description and with the random walk equation [21,22], due to the presence of steep potentials. On the other hand, the numerical solution is required if accurate theoretical estimates of experimentally observable quantities are necessary. A variety of numerical methods are available for solving diffusion-Ike equations. As an example, the brownian molecular dynamics is often used in the studies [1,23,24] of the conformational motion in polymers. In this case, one generates a representative ensemble of stochastic trajectories from the corresponding Langevin equation. Other methods are based on the calculation of the experimental observables from the matrix representation of the time evolution operator. In particu-
0301-0104/86/$03.50 0 Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
90
G. More / Non-orthogonal
ham Junctions in diffusional problems
lar, the Lanczos algorithm [25-271 has proven to be a valuable tool when dealing with spectroscopic or scattering-type measurements, where the spectral densities represent the observables of the physical system [28,29]. The advantages of the method have already been discussed together with its relation with the continued fraction representation encountered in statistical physics [30]. The standard implementation of the Lanczos algorithm uses an orthonormal set of basis functions that normally are the eigenfunctions of the diffusion operator in the absence of the potential. However, as the potential increases in strength, a larger basis set and a progressively heavier computation effort are needed for reproducing the eigensolutions of the operator. In order to illustrate the last consideration, let us examine the equilibrium distribution function which is the stationary eigensolution of the diffusion equation. In the presence of strong potentials, the equilibrium distribution is a sharply peaked function, since it is non-negligible only near the minima of the potential. If this eigensolution is expanded on a set of functions having a smooth dependence on the coordinates (say the circular functions for the planar rotor problem or the Wigner rotation matrices for the rotational diffusion of a molecule), a large number of terms is necessary to ensure convergence. With increasing potential strengths, the equilibrium distribution tends to resemble a delta function (or a sum of delta functions), and so the basis set must be infinitely expanded to represent accurately the stationary solution. The growth of the basis size with the strength of the potential constitutes the main handicap of the numerical solution of diffusional models in the study of hindered motions. It would be possible to overcome this difficulty, if one can use basis functions that mimic the known asymptotic solutions of the diffusion operator. In such a case, it is reasonable to expect that few basis functions are required for an accurate calculation of the spectroscopic observables. In a short communication [31], a set of functions having this feature, the so-called modified Wigner functions set, has been proposed for the specific problem of the reorientational motion of rigid molecules in nematic liquid
crystals. The asymptotic behaviour of the solutions of the diffusion equation was reproduced by the equilibrium distribution taken as common factor of the basis functions. Since these functions are not mutually orthogonal, an extension of the Lanczos algorithm to a non-orthonormal basis set was reported in the same work. A related approach has been considered by Blackmore and Shizgal [19]. They have solved the diffusional problem for a particular bistable potential on a polynomial basis, by taking the equilibrium probability as weight function. Their treatment, however, is based on the analytical orthogonalization of the polynomial basis. The orthogonalization procedure is not easily applicable to problems characterized by potential functions which are not expressed as simple polynomials, but rather by trigonometric functions. In general, one should use basis functions defined in terms of the equilibrium distribution, but avoid the unnecessary analytical orthogonalization by adopting the numerical procedure naturally provided by the modified Lanczos algorithm [29,31]. In some recent studies [21,22] concerning the transition processes amongst potential wells, another type of non-orthogonal functions, called localized functions, was considered. These functions represent the minimal basis set required to generate the random walk equation from a complete description of the diffusion motion, in the limit of infinite potential barriers. With finite potential barriers, the localized functions allow one to obtain a very good approximation of the low-frequency branch of the spectrum of the diffusion operator. On the other hand, the use of the localized functions for a “numerically exact” (i.e. within an arbitrary accuracy) calculation of the spectral densities was not considered in refs. [21,22]. The use of basis functions which mimic the asymptotic behaviour of the solutions of the diffusion equations constitutes a promising method in treating molecular hindered motions, since it should minimize the size of the matrices to be handled. One would like to have orthogonal bases in order to apply standard numerical techniques. But this in general is unlikely to occur, unless particular forms of the potential are chosen [19].
G. Moro / Non-orthogonal basis funciions in diffusional problems
Therefore, the calculation of the observables from non-orthogonal basis sets appears a necessary task. Actually, one can resort to the Gram-Schmidt orthogonalization procedure applied to the ensemble of functions having the correct asymptotic behaviour, but an alternative and preferable route is represented by a suitable modification of the Lanczos algorithm [31], because of its intrinsic efficiency in computing spectral densities [28,29]. We recall here that the Lanczos algorithm itself is an orthonormalization method applied to the moments of an operator acting on a given starting vector [32]. In this article, the use of non-orthogonal functions in the planar rotor problem will be analysed in detail. The choice of this particular problem has been dictated, apart from its relation to chain dynamics [2-41, by the simple form of the diffusion equation, in the belief that the advantages and the disadvantages of the method can thus manifest clearly. Of course the ultimate and future goal is the application of the method to complicated problems, like multidimensional diffusion or Fokker-Planck equations, where one cannot afford the standard implementation of the Lanczos algorithm in presence of a strong potential because of the extent of the basis set. In section 2, the modified Lanczos algorithm will be presented in the context of the calculation of observables given by spectral densities. The application of the method to single minimum potentials will be considered in section 3. It will be shown that a choice of basis functions closely related to those already presented in the context of rotational motion in nematics [31], is adequate to efficiently solve the problem. These functions, in fact, are suitable to represent small amplitude motions around the minimum of the potential. The diffusion equation with a double minimum potential will be treated in section 4. In this case, two independent physical processes are present: the small-angle displacements about the potential minima and the transitions between potential wells. It will be shown that when dealing with spectral densities influenced by the transition kinetics, we need to insert among the basis elements functions which are similar to the localized functions previously considered [21,22]. The per-
91
formances of the non-orthogonal basis functions will be discussed and compared with the results of the standard orthogonal representation of the diffusion operator.
2. The Lanczos
algorithm with non-orthogonal
basis functions
The spectral density enters as a fundamental quantity in the interpretation of spectroscopic [33,34] or scattering-type [35] measurements. In particular its frequency dependence determines the observable effects of the molecular motion. Therefore the calculation of spectral densities is an important tool in providing an “interface” between the experiments and the statistical physics which supplies the possible models of molecular motion. The Lanczos algorithm [25-271 represents a convenient method for calculating spectral densities since their frequency dependence is generated from the continued fraction representation [28-301, in analogy with the projection method proposed by Mori [36]. Moreover the intrinsic efficiency of the Lanczos algorithm with respect to the computer resources (i.e. computation time and size of the memory) is stressed when the spectral lineshapes are calculated instead of the eigenvalues [28,29]. The calculation of spectral densities by means of the Lanczos algorithm with an orthonormal basis has been already analyzed in the past [28]. Here we shall summarize its application to a physical problem described by the classical stochastic variable x = x(t). For a given function F(x) representing the observable probed in a spectroscopic experiment, the spectral density J(w), that is the Fourier-Laplace transform of the correlation function F( x (t ))*F( x(0)) [33,34], is written as J(w)=(FI(io+r)-‘JFP),
(I) where P(x) represents the equilibrium distribution function and the scalar product is defined as an integration on the variable x. In the relation (l), r denotes the time evolution operator for the non-equilibrium probability p ( x, t )
G. More / Non-orthogonal basis functions in diffusionalproblems
92
a+,
g/at=
-Q+,
(2)
t).
We consider only markovian and stationary processes, so the equilibrium distribution function is the unique stationary solution of r, rP=Q.
(3)
In the calculation of spectral densities, venient to use the symmetrized time operator defined in the relation
it is conevolution
i; = p-wrpw
(4)
There are different ways to implement the Lanczos procedure, depending on the symmetry of i; [29,37]. In the following we shall specialize the method to self-adjoint operators, like the diffusion operator with reflecting or periodic boundary conditions. For the spectral density given in eq. (l), an orthonormal set of functions c$~(x) is generated by means of the Lanczos method, starting from function @r(x) defined as +1 =
FP"2(
I’) -1’2,
1F
where the average 1F I* is calculated with the equilibrium distribution function P(x). The following three-term recursive relation is written [29,30]:
where CY,is the diagonal matrix element of f in the I#+ representation and &+i is the real coefficient derived from the normalization conditions applied to r&+ 1. The operator p has a tridiagonal matrix representation on the +k basis, therefore the spectral density can be written in terms of the following continued fraction [38]
iw+cu,
-
I P2
I2
io+cu,-
: -
I I% I 2 iw + (Ye-
(7) In the usual
implementation
R,,=
of the Lanczos
represet of
(8)
If’lfd.
(f,
By indicating with xk the column matrix defined by the expansion coefficients of C#Qin the f/ basis, (9)
G!i = C(4,fi i the recursive P n+lX,+l=
relation P
-
(6) becomes
%l)%
-
P,LI.
(10)
This is the standard Lanczos relation which can be easily implemented in a computer program for calculating the ensemble of coefficients C.X~ and flj [26-281. The frequency dependence of the spectral density is then derived from the continued fraction of eq. (7). In a recent communication [31], a modification of the Lanczos algorithm that permits the use of a non-orthogonal basis set, was proposed. We shall now analyse in more detail the implementation of this new method, with particular emphasis on the truncation of the, in principle, infinite set of basis functions. The truncation problem was not considered in the previous communication. Let f/ denote the non-orthogonal basis functions, and S the corresponding normalization matrix, s,j'=
(11)
lh’>.
It is implicitly assumed that the functions f/ are linearly independent and that they constitute a complete set for representing the eigensolutions of f. The functions r’m can be formally decomposed as
(12)
Ffm= Cn/r,mfn. m
1 J(4/lF12=
algorithm [28,29], one considers the matrix sentation R of r on a given orthonormal basis functions 6,
It is easily shown that eq. (10) continues to hold, with the column matrices x, defined according to eq. (9) if the matrix R is replaced by the matrix M, that is P n+1x,+1=
CM
It is convenient
(13)
-%1)%-&x,-I-
to introduce
the auxiliary
column
G. Moro / Non-orthogonal basis functions in diffusional problems
matrices
yk,
9N=
Yk=sq,
(14
such that the normalization condition on +k and the diagonal element (Ye can be calculated from the relations Y,$k
= 1,
CQ = yiMx,.
05) 06)
The sequence of operations of the standard Lanczos algorithm [26-281 can easily be adapted to the previous relations. In each iteration, one starts with the column matrices x~_~, xk and yk and proceeds as follows: (1) calculate simultaneously CQ and &+ ,xk + , , that is the right-hand side of eq. (13); (2) calculate from eq. (14) the column matrix &+l Yk+l; (3) calculate, by taking into account eq. (15), P k + , and normalize xk + 1 and yk + , . Note that the first step can be implemented like in the standard Lanczos algorithm by superimposing the column matrix Pk+l~k+l to x~_~. Moreover, it is clear from the sequence of operations, that only one array is necessary for storing column matrices y, and yk+i. Therefore, the storage is increased by one column matrix (yk) and one square matrix (S) with respect to the standard Lanczos algorithm. As regards the time-consuming operations, the number of multiplication of square matrices by column matrices is doubled, as effect of computing yk+, at each iteration. Of course the global efficiency of the algorithm, in both the standard and the modified form, is strongly dependent on the sparsity of the square matrices [26-281. The set of functions fj should be infinite if it must satisfy the condition of completeness with respect to the ensemble of the solutions of the diffusion operator. Of course, in computer applications, only finite matrices can be handled, in correspondence of a subspace cN spanned by functions fi for j = 1,2,. . . , N. Therefore the operator f in eq. (6) should be substituted with its projected form
E Ifi
)(SN1)j,'(f/'
93
I
9
(18)
j,j'=l
where S, is the first N X N block of the normalization matrix. Also the function FP*/*, which determines C#Qaccording to eq. (5), should be projected onto eN, but usually it is convenient to include the function FP”* in the finite basis set. In eq. (13) which implements the Lanczos procedure, the infinite matrix M is then replaced by the N x N matrix M N given implicitly by the equation
(19) One easily recognizes that the elements of M, are in general different from the corresponding elements of M, therefore M, cannot be generated simply by truncating M. By taking into account eq. (18) for the projection operator, MN can be explicitly written as M, = S,‘R,,
(20)
where R, is the first N X N block of the matrix defined in eq. (8). The use of eq. (19) for calculating the elements of M, is convenient when p’n can be analytically written as a linear combination of basis functions (see the next section for an example). Otherwise M,,, must be derived from eq. (20). Jones and co-workers in their calculations of the electron charge density in solids [39,40], have implemented the Lanczos algorithm according to the second procedure. In particular they avoided the timeconsuming inversion of S, by computing, at each iteration of eq. (13), the array S,‘R,x,, from the solution of the linear equation system with coefficients given by S,.
3. The single minimum planar rotor The time evolution operator r for the diffusional motion of the variable y describing the rotation of a dumbbell in a plane, is written as r=
-(a/ay)PD(a/ay)P-‘,
(21)
G. Mom / Non-orthogonal
94
basis functions
where D is the diffusion coefficient and P(y) is the equilibrium distribution function that can be defined in terms of the potential V(y) acting on the dumbbell, P(y)
= exp( - V/k,T)/L2”exp(
- V/k,T)
This type of equation is encountered in the theory of internal motion of molecular conformers [3]. We do not consider here the angular dependence of D; however, once a model for D(y) has been chosen, the following analysis can be easily adapted to the specific case. In this section we shall analyse the case of a symmetric potential (23)
having a single minimum centered at y = 0 and with V(0) = 0. The asymptotic expansion of both static (mean values) and dynamic observables (correlation functions or spectral densities) are obtained by introducing a fictitious parameter h (taken as unity at the end) as factor of the potential I/ and considering the limit for X -+ 00. The value of the observables in the limit of highly anisotropic potential is then derived from the leading terms in the expansion in powers of X-‘12. The symmetrized diffusion operator, in the approximation of zeroth order, becomes F = zJ(x2-
I -
a2/ax2),
(24)
where y = { XV”(0)/(2k,T)}
-“2x,
u = DXV”(0)/(2k,T),
(254 (25b)
and V”(0) indicating the second derivative of the potential calculated in the location of the minimum. The corresponding asymptotic eigenfunctions q,* are written in terms of the hermitean functions H,(x), \k,”= ( AV”(0)/2nk,T)“4 \k: = ( j! 2’) -I’*H,(x)
exp( -x2/2), q;(x).
problems
potential minimum. The asymptotic behaviour of the spectral density given in eq. (l), is then derived from the following representation of FP”‘: (27)
dy. (22)
V(Y) = v(-Y)
in diffusionul
(26a) (26b)
Such approximate solutions of the diffusion operator represent the small-angle motions around the
where coefficients cj are calculated from the Taylor expansion of F(y) around y = 0. Since the coefficient c, is of the order of h-‘J-“/2, the functions !$A contribute to the spectral density with decreasing weights given by powers of l/X. In other words, by considering the approximated diffusion operator of eq. (24) in the subspace e”, spanned by functions $, !I$, . . . , !Pi-‘, that is in the ensemble of polynomials of (N - 1)th degree in the variable x multiplied by @, one takes into account in the spectral density the first N contributions with respect to the powers of the asymptotic parameters 1 /A. Although the functions q; reproduce well the behaviour of the solutions of the diffusion operator around y = 0 with steep potentials, they are not useful in numerical calculation of spectral densities. In fact they cannot be used for generating a matrix representation of the complete diffusion operator of eq. (21), because the periodic boundary condition is not ensured. In ref. [31], a set of functions that mimic the asymptotic solutions, still preserving boundary conditions, has been proposed for the rotational diffusion of an axially symmetric molecule in nematic liquid crystals. The same method can be applied to the simpler diffusion problem considered here, in this way enlightening its rationale. The basic idea is to take advantage of the knowledge of the stationary solution of r, that is the subspace CL previP(Y) ‘I* 3 for approximating ously defined, with a subspace zN made by functions which satisfy periodic boundary conditions. in all generThe subspace cN can be constructed ality from the linear combination of functions f, = g, P’/~, fi = g, P112,. . . , fN = gNP112, where the g,(y) are a set of functions having a well defined Taylor expansion around y = 0. Of course subspaced ek and eN are not strictly comparable, because they refer to functions having different boundary conditions. Nevertheless, they are equivalent in the asymptotic limit. In fact, by making
G. More / Non-orthogonal
basis functions in diffurional problems
the substitution eq. (25a) and expanding functions up to order (N - l), the fi as powers of X-l/’ ensemble of functions pN_I(x)!Pt(x), where pN_ 1( x) indicates a generic polynomial of (N 1)th degree in the parameter x, is recovered. Therefore, the leading contributions with respect to the asymptotic parameter l/h are taken into account, and the ensemble of functions fJ previously defined represents a convenient basis set in the calculation of the spectral densities. We emphasize that a precise correspondence between each function fi and the asymptotic eigenfunction @,” is not necessary, since our objective is to provide an optimal subspace to the numerical solution of the diffusion operator by means of the Lanczos algorithm. Apart from the condition of linear independence, we have left much freedom in choosing the functions g,. The eigenfunctions of the isotropic diffusion operator, that is the sine and cosine functions, are a natural choice in order to represent conveniently the behaviour of the operator also in the limit of weak potentials. The resulting functions fJ are non-orthogonal. This is the price that must be paid in order to mimic the asymptotic behaviour of the solutions of the diffusion operator. On the other hand, the modified Lanczos algorithm discussed in the previous section, can handle such type of basis functions. In the remainder of this section, we shall discuss in detail the application of the method to the calculation of the spectral density of the function F(y) = cos y - cos y, with a cosine-type potential
composed
as follows,
=_i’Df-
p&
(A/2)(1
- cos y).
The following set of non-orthonormal tions has been used: f, = (cos( jy)
- cos( jy
)) P”2.
(jDA/4)(f,+i
(30)
(31) It is known that the Lanczos algorithm can reproduce correctly the spectral density with a number of iterations less than the dimension N of the matrix [28,29]. Also in the present calculations we have found the same behaviour, but in the follow-
Log
E
10 N
0 0
(28)
0 OO
basis func-
0 -2--
0
n
0
(29)
Evans, in a recent analysis of the isomerization of butane [4], has considered a basis set related to the functions of eq. (29). However, he calculates the matrix representation of the time evolution operator with the corresponding orthonormal set generated by means of the Gram-Schmidt orthonormalization procedure. Since the function l;fi can be analytically de-
-&-i),
the elements of M, are conveniently computed by means of eq. (19). In particular, the projection should be taken into account exoperator 9, plicitly only for j = N, where the solution of a linear system of equations having (fN+, I&) as constant terms, is required according to eq. (18). For characterizing the efficiency of the method, a criterion is needed for the error EN of the approximated spectral density JN (w), calculated with a N-dimensional basis set, with respect to the exact J(w). Like in ref. [28], we use the normalized difference in the area of the real part of the spectral density,
0 --
V/k,T=
95
-4--
0
n
0
.
5
10 N
Fig. 1. Relative error in the area’ EN (logarithmic scale) as function of the size N of the basis set, for the spectral density relative to cos y - cos y with single minimum potential eq. (28) and A = 20. (0) Orthogonal basis functions; (A) nonorthogonal basis functions.
96
G. More / Non-orthogonal basis functrons in diffusional problems
ing discussion we disregard this aspect by letting n = N. In fig. 1 the error E, is plotted as function of N for both the non-orthogonal and orthogonal basis set. In the second case the standard Lanczos algorithm has been used with a basis set of normalized cos(jy) functions. The value of A used (A = 20) corresponds to a very steep potential. A large number of orthogonal basis functions is necessary for obtaining a reasonable accuracy in the spectral density, because many cosine functions are needed to represent FP’/‘, as a consequence of the factor P ‘I2 characterized by a sharp peak centered at y = 0. A comparable accuracy in the results is achieved with one or two non-orthogonal basis functions only. The growth of the size of the orthogonal basis set with increasing A, is manifest from fig. 2 where the least N which ensures an error EN G 10.e3, is plotted against A. On the contrary, the number of non-orthonormal basis functions decreases with A, starting from A = 9. In fact one can easily show that function f, alone predicts the leading contribution to the spectral density in the limit of an infinite A. The relative advantages of the two methods can be understood from the trends shown in fig. 2, by taking into account that the modified Lanczos algorithm discussed in the previous section, is
always less efficient than the standard Lanczos algorithm in treating matrices of the same dimension, as a consequence of algebraic operations involving the normalization matrix S,. Therefore the use of non-orthogonal functions becomes convenient only when it allows a reduction of the basis set. It is clear from the behaviour displayed in fig. 2, that this condition should be always reached with a sufficiently large value of A. Of course stronger potentials lead to even greater savings in computation time. On the other hand a general criterion for the relative convenience of the non-orthogonal basis set with respect to the, normally used, orthogonal ones cannot be given, for it depends on the sparsity of matrices M, and S, appropriate to a particular diffusional problem. Finally, we mention that also in the case of odd functions of y, trends like those shown in figs. 1 and 2 are obtained with non-orthonormal functions f, = sin( jy) PI/‘. 4. The double minimum planar rotator In this method of applied to specifically potential V/k,T=
8
0
0
n
n
0 N 0 4
0 0
n n
n
n LJ 0 :;;;;; 0
-10
1 A
20
Fig. 2. Least-size N of tie basis set that assures a 10m3 accuracy in the calculation of the special density relative to cost -COS Y, as function of the parameter A of the potential given in eq. (28). (0) Orthogonal basis functions; (A) nonorthogonal basis functions.
section we intend to show how the non-orthogonal basis functions can be multiminima potentials, by considering the planar rotor with the two-minima
(A/2)(1
- cos 2y) = A sin2y,
(32)
with a positive A. Since the diffusion operator is invariant under symmetry operations of the point group CzV, the calculation of the spectral densities can be factorized in the four functional subspaces pertinent to the irreducible representations. The asymptotic analysis of the small amplitude motions can be applied independently to each potential minimum and, in strict analogy with the treatment illustrated in section 3, suitable sets of nonorthogonal basis function are generated. These functions can be factorized in the four subspaces as follows: A, :
f, = [ cos(2jy)
A,:
f, = sin(2 jy) P”2,
- cos(2jy)]
P1’2,
G. More / Non-orihogonal basis functions in diffusional problems
B, :
f, = cos[(2 j - l)y]
PI”,
B,:
f, = sin[(2 j - l)y]
PI”,
where the index j runs over positive integer numbers. We do not report here specific results for functions belonging to representations A,, A, and B,, since the same trends as discussed in the previous section are recovered and thereby identical conclusions can be drawn. A completely different behaviour is found instead with functions of B, symmetry. In fig. 3 the error EN on the spectral density relative to function F = cos y is displayed for some values of the size N of the basis set in both cases of orthogonal (circles) and nonorthogonal (triangles) functions. The ensemble of functions 7~~~‘~ cos(2j - l)y, with j = l-N, was taken as orthogonal basis set. A large value of A (A = 12) was used in order to report the behaviour for highly anisotropic systems. A striking difference with respect to fig. 1 is manifest: only a little improvement on the calculation is achieved when passing from an orthogonal to a non-orthogonal basis set with the same dimension. It has been also found that the size of both types of basis sets must be increased with larger values of A, in order to calculate spectral densities with the same accuracy. A similar behaviour is reported in the Log
0--
N
A
A
O
A
0
A
0 A A
-2 --
work of Blackmore and Shizgal about the numerical solution of the diffusion equation of a bistable potential (compare, in particular, tables II and III of ref. [19]). In fact, as mentioned in section 1, an analogy exists between our non-orthogonal functions fi and their orthogonal polynomials constructed with the equilibrium distribution as weight function. The inefficiency of the non-orthogonal functions of B, symmetry can be explained by taking into account the peculiar features of the bistable system. There is, in fact, an eigensolution ?P, of the diffusion operator having the eigenvalue (Y, which vanishes exponentially like an Arrhenius factor exp( - E/k,T), with an increasing energy barrier E between the minima of the potential [13,14]. The solution q, is representative of the slow process of transition between the potential wells. The decay behaviour of the correlation function associated to cos y is dominated by the root (Y, which, instead, does not contribute for symmetry reasons to the correlation functions relative to the other representations A,, A, and B, [21]. In fact *I belongs to the B, representation because it is the antisymmetric counterpart of the stationary solution ‘k, = P”‘. It can be more conveniently written as *, = g(u>P”‘,
(33)
where, in the limit of high-energy barriers, g(y) has a step-like behaviour, that is g - 1 (g = - 1) for A in a wide interval centered at y = 0 (y = n) and it changes abruptly for y = +~/2. In analogy with the bistable system defined in an infinite domain [13], the asymptotic analysis leads to a g(y) with an error function shape,
E I0
97
0
g(y)r=
-erf[A1/2(y-7r/2)],
(34)
A q A -4
0
-rJ 0
I
I
A
I
5
10 N
Fig. 3. function relative A = 12. functions
Relative error in the area EN (logarithmic scale) as of the size N of the basis set for the spectral density to cos y. with the bistable potential eq. (32) and (0) Orthogonal basis functions; (A) non-orthogonal j,; (0) basis set with f, and functions h,.
for 0 6 y i 71. It is evident that the function f, = *i only in a (cos y ) P”2 mimics eigensolution small interval around y = 0 (y = IT) where cos y = 1 (cos y = -l), but it is unable to predict the error function behaviour near potential maxima. Remarkable improvements are not expected from the other functions 6.. The inefficiency of the non-orthogonal functions A in the calculation of spectral densities relative to the function cos y should be ascribed
to their inadequacy in representing *i. Therefore, one needs to include, in the basis set, functions which approximate the eigensolution !P,, that is functions describing the transition between the potential wells. In some recent works 121,221, a class of such functions, called localized functions, has been introduced. The localized functions are constructed by matching in the locations of the potential minima, the functions g,(y)
= j%(+‘da, (35) ?J which are locally exact solutions of the differential equation PgP ‘P = 0.
(36)
One easily recognizes that the functions g, of eq. (35) behave like the error function of eq. (34) in correspondence with the maxima of the potential. Moreover, the ensemble of localized functions have been identified in refs. [21,22], with the least set of functions which generates the random walk equation amongst the stable states (i.e. states localized at the potential wells) from the diffusion equation. In the present work we prefer, instead, to work with a generalized form of localized functions which have continuous derivatives. They are obtained from the approximation of function g( y > g(y) = ri2G(cr)
P(cu)-‘dru,
(37)
with function G(a) obeying to the following conditions: (i) G(a) must belong to symmet~ subspace %* (ii) G(a) must behave like a non-vanishing constant near the potential maxima. The former ensures that g(y) belongs to the B, irreducible representation while the latter leads, in the asymptotic limit, to the correct error function behaviour of eq. (34). All the functions sin(2j 1)a satisfy these conditions, and one can insert the following functions hi, hi(y)=P’/‘III/‘sin(2jY
l>a P(a)-‘da,
among the elements of the basis set. h, constitutes the best trial function
(38)
for @,,
because the eigensolution cos y with the lowest non-vanishing eigenvalue is recovered for A = 0. Since both limits of small and large values of A are reproduced, one expects that also in the middle range of potential strengths h, approximates well the behaviour of !P,. As a test, we have calculated the approximate eigenvalue (Y~from the expectation value of i; with function h,, ai = (hi I f’l W/‘(h, 14) = D((sin2y}~(y)-‘)/(h,
I h,),
(39)
and compared it with the true eigenvalue obtained from numerical diagonalization of I=. In fig. 4 the results are displayed as a function of A, together with the leading term of the asymptotic expansion of ai, CY,= (4DA/a)
e-‘.
(40)
The numerical results for some selected values of A are reported in table 1. It can be clearly seen that relation (39) yields a good estimate, within an error of I%, of the eigenvalue (Y, in the whole range of potential strengths. h, can now be used, together with the functions 4, for calculating the spectral density by means of the modified Lanczos algorithm. Of course f, must always be included in the basis set, since it represents the starting vector in the calculation of the spectral density relative to function cos y. Since I
I
I
0
I
b.
1
2
3
Fig. 4. Lowest positive root IX, of the diffusion operator with the bistable potential given in eq. (32), as function of the parameter A (diffusion coefficient D =I). (0) Exact values: approximated values calculated from eq. (39); (- - -) c-) asymptotic relation eq. (40).
G. Moro / Non-orthogonal basisfunctionsin dijfusional problems Table 1 Lowest positive root of the diffusion operator with the bistable potential of eq. (32), as function of A (diffusion coefficient D =l) A
Exact CX,
Approximate
0 1 2 3 4 5 10 20
l.OtM 0.5876 0.3237 0.1676 0.8198 x 0.3823 x 0.5463 x 0.5112~
1.000 0.5890 0.3257 0.1689 0.8264x 10-l 0.3850 x10-’ 0.5473 x10-3 0.5113 x 10-r
10-l 10-t lo- 3 lo-’
(Y, ‘) (0.0) (0.4684) (0.3446) (0.1902) (0.9328x 10-l) (0.4290 x 10 - ’ ) (0.5781 x 10-3) (0.5249x10-‘)
a) (Y, calculated by numerical integration of eq. (39). asymptotic values obtained from eq. (40) are reported tween parentheses.
The be-
one cannot extend, in a simple way, eq. (30) to the function h,, the matrix M, is calculated from eq. (20). The efficiency of this type of non-orthogonal basis functions appear clear from the results shown in fig. 3. The two-dimensional basis set constituted by functions h, and fi, gives a spectral density as accurate as that obtained with seven functions fi alone. The addition of other functions hj to the basis set improves considerably the accuracy of the calculated spectral density. In fig. 3 only the result for the basis set formed by ft, h, and h, is reported, because wider sets obtained by adding more hj functions lead to a value of EN below the reported scale. We mention also that the efficiency of the basis sets made with functions h,, becomes more pronounced when the parameter A is increased. Therefore, all considerations of the previous section about the convenience of the modified Lanczos algorithm with respect to its standard form, can be extended to the bistable diffusion equation. The results reported in table 1 prove the validity of our choice of basis functions. In fact, it shows that a single function is sufficient to calculate with reasonable accuracy the lowest positive root. We remind that this eigenvalue determines the long-time behaviour of the bistable system and it is identified with the observed rate in model system for the chemical kinetics [8,9,21]. One could compare this result with that obtained by Black-
99
more and Shizgal [19] in their study of a related bistable system in one dimension. From tables II and III of ref. [19], it appears that 10 and 50 orthogonal polynomials are necessary to compute the same root with a comparable accuracy, for barrier heights of 2.5 and 25 k,T units respectively. One can conclude that a reduction of orders of magnitude in the size of the basis set can be achieved. by properly taking into account the asymptotic behaviour of the diffusion equation.
5. Conclusion The results presented in this work show that, in regimes of strongly anisotropic potentials characterized by large barriers, a proper choice of the non-orthogonal basis functions drastically reduces the subspace in which the diffusion operator must be solved numerically. In truth, this conclusion is supported so far by applications to simple diffusional problems, like the planar rotator considered here or the rotational motion in nematic liquid crystals analyzed in ref. [3]. On the other hand, our objective was to develop the methodology, leaving to future works its application to problems characterized by large-scale computations. However, we expect that also in complex systems, like those described by multidimensional diffusion equations, it should be possible to construct from the asymptotic analysis of the time evolution operator [18], a proper set of functions that mimic the fundamental processes determining the spectroscopic observables. The examples considered in this work tell us that the non-orthogonal functions must be chosen careffilly if the optimal basis set is required. In particular, the use of the equilibrium distribution probability as weight function, is sufficient only when dealing with the small amplitude motions around the minima of the potential. If a given spectral density is influenced by the transition process amongst potential wells, then one should include in the basis set specific functions which generate the random walk equation from the more general diffusion equation [21,22]. The proposed modification of the Lanczos algorithm allows the calculation of the spectral den-
100
G. More / Non-orthogonal
ham Junctrons
sities by means of non-orthogonal functions. One should realize that the computational procedure is more complicated than the standard Lanczos algorithm. Therefore the use of non-orthogonal functions becomes convenient only when it leads to a reduction on the size of the matrices to be handled. The numerical results previously discussed show that such a situation is always reached by increasing the energy barrier of the potential. Finally, we would like to suggest two fields of future applications of the non-orthogonal functions method with the modified Lanczos algorithm: (i) the Fokker-Planck equation in the low friction regime where very large basis sets are needed (see refs. [41,42] for an analysis of the asymptotic behaviour of the Fokker-Planck operator); (ii) the multidimensional diffusion equation pertinent to chain dynamics [l], where it is immediately evident the analogy between the planar rotor problem and the motion around each internal degree of freedom of the molecule. The chain dynamics is a manifest example supporting the need for efficient methods for representing diffusion operators. Let us consider the internal motion of the butane molecule under the action of the Ryckaert-Bellemans potential [22]. We have found that a minimal basis of 20 elements is necessary when the standard circular functions are used. Therefore, we expect that matrices of the order of 20 N should be solved for alkyl chains having N internal degrees of freedom. In this way, problems with more than three degrees of freedom could not be handled even with large computers. It follows that only with an accurate choice of the basis functions, which minimizes the representation for a single degree of freedom, it will be possible to tackle the complete analysis of the dynamics of large molecules. In conclusion, it is conceivable that, by properly taking into account the asymptotic solutions of the diffusion operator, very few basis functions might be strictly necessary to compute the dynamical observables of a physical system.
n diffusiond
problem
Acknowledgement This work has been supported by the Italian Ministry of Public Education and in part by the National Research Council through its Centro Studi sugli Stati Molecolari Radicalici ed Eccitati. The author is indebted to Professor P.L. Nordio for stimulating discussions and for advice on the revision of the manuscript.
References PI M. Fixman, J. Chem. Phys. 69 (1978) 1527, 1538. 121J. Skolnick and E. Helfand, J. Chem. Phys. 72 (1980) 5489. [31 D.C. Knauss and G.T. Evans, J. Chem. Phys. 72 (1980) 1499. 141 G.T. Evans, J. Chem. Phys. 78 (1983) 4963. 151 G. Moro, U. Segre and P.L. Nordio, in: Nuclear magnetic resonance of liquid crystals, ed. J.W. Emsley (Reidel, Dordrecht, 1985). [61 G. Moro and P.L. Nordio, J. Phys. Chem. 89 (1985) 997. rotations in molecular crystals [71 W. Press, Single-particle (Springer, Berlin, 1981). 181 H.A. Kramers, Physica 7 (1940) 284. 191 D. Chandler, J. Chem. Phys. 68 (1978) 2959. HOI J.L. Skinner and P.G. Wolynes, J. Chem. Phys. 69 (1978) 2143. [JLI B. Carmeli and A. Nitzan, J. Chem. Phys. 79 (1983) 393; 80 (1984) 3596. 1121 R.F. Grote and J.T. Hynes, J. Chem. Phys. 77 (1982) 3736. 1131 R.S. Larson and M.D. Kostin, J. Chem. Phys. 69 (1978) 4821. 1141 H. Tomita, A. Ito and H. Kidachi, Progr. Theoret. Phys. 56 (1976) 786. I151 B. Caroli, C. Caroli and B. Roulet, J. Stat. Phys. 21 (1979) 415. [161 K. Lindeberg and V. Seshadri, J. Chem. Phys. 71 (1979) 4075. [171 J.M. Deutch, J. Chem. Phys. 73 (1980) 4700. of stochastic differen[181 Z. Schuss, Theory and applications tial equations (Wiley, New York, 1980). 1191 R. Blackmore and B. Shizgal, Phys. Rev. A31 (1985) 1855. [201 W. Nadler and K. Schulten, J. Chem. Phys. 82 (1985) 151. 1211 G. Moro and P.L. Nordio, Mol. Phys. 56 (1985) 255. I221 G. Moro and P.L. Nordio, Mol. Phys. 57 (1986) 947. Chem. 1231 R.M. Levy, M. Karplus and J.A. McCammon, Phys. Letters 65 (1979) 4. 1241 T.A. Weber and E. Helfand, J. Phys. Chem. 87 (1983) 2881. [25] C. Lanczos, J. Res. Nat]. Bur. Std. 45 (1950) 255; 49 (1952) 33.
G. More / No-o~~~ogonai basis functions in ~iff~ionaI problems
[26] B.N. Parlett, The symmetric eigenvalue problem (Prentice-Hall, Englewood Cliffs, 1980). (271 J. Cullum and R. Willoughby, Lanczos algorithm for large symmetric eigenvalue computations, Vol. 1. Theory (Birkhauser, Bascl, 1985). (281 G. Moro and 3.H. Freed, J. Chem. Phys. 74 (f981) 3757. [29] G. Moro and J.H. Freed, The Lanczos algorithm in molecular dynamics: calculation of spectral densities, Proceedings of the Workshop in Large Eigenvalue Problems, Oberlech, Austria (1985), eds. J. Cullum and R.A. Willoughby. 1301 G. Moro and J.H. Freed, J. Chem. Phys. 75 (1981) 3157. 1311 G. Moro and P.L. Nordio. Chem. Phys. Letters 96 (1983) 192. [32] Yu.V. Vorobiev, Method of moments in applied mathematics (Gordon and Breach, New York, 1965).
101
[33] R. Kubo, J. Phys. Sot. Japan 12 (1957) 570. [34] B.J. Beme, in: Physical chemistry, Vol. 8B, eds. H. Eyring, D. Henderson and W. Jost (Academic Press, New York, 1971) p. 540. [35] L. van Hove, Phys. Rev. 95 (1954) 249. [36] H. Mori, Progr. Theoret. Phys. 34 (1965) 399. [37] W.A. Wassam, 3. Chem. Phys. 82 (1985) 3371, 3386. [38] H.S. Wall, Analytic theory of continued fractions (Van Nostrand, Princeton, 1948). [39] R. Jones and M.W. Lewis, Phil. Mag. 49B (1984) 95. [40] R. Jones, in: The recursion method and its application, eds. D.G. Pettifor and D.L. Weaire (Springer, Berlin, 1985) p. 132. [41] H.D. Vollmer and H. Risken, Physica 1lOA (1982) 106. [42] G. Moro, Chem. Phys. 96 (1985) 287.