Jens Gustedt

Futex based locks for C11's generic atomics *

Keywords: lock primitives, atomics, C11, futex, Linux primitives de verrouillage, opérations atomiques, C11, futex, Linux

We present a new algorithm and implementation of a lock primitive that is based on Linux' native lock interface, the futex system call. It allows us to assemble compiler support for atomic data structures that can not be handled through specific hardware instructions. Such a tool is needed for C11's atomics interface because here an Atomic qualification can be attached to almost any data type. Our lock data structure for that purpose meets very specific criteria concerning its field of operation and its performance. By that we are able to outperform gcc's libatomic library by around 60%.

INTRODUCTION

Only very recently (with C11, see JTC1/SC22/WG14) the C language has integrated threads and atomic operations into the core of the language. Support for these features is still partial: where the main open source compilers gcc and Clang now offer atomics, most Linux platforms still use glibc as their C library which does not implement C11 threads. Only platforms that are based on musl as C library, e.g Alpine, are feature complete.

The implementation of the C11 atomic interface typically sits in the middle between the implementation of the core language by the C compiler and the implementation of the C library. It needs compiler support for the individual atomic operations and a generic atomic lock interface in the C library when no low-level atomic instruction is available.

Since Linux' open source C libraries do not implement the generic interface, the compilers currently provide a library stub that implements the necessary lock by means of a combination of atomics and POSIX pthread_mutex_t.

By consensus in the community an API interface has emerged that both, gcc and clang, implement.

These library stubs can rely on all their knowledge of the architecture specific atomics for the non-contended part of the implementation, finally that is the context in which they were created. But since they are generic, in case of contention these library stubs rely on generic OS interfaces to put threads to sleep. From a point of view of the C language the natural interface for that would be C11 threads, but since these are still missing in glibc, they fall back to POSIX threads for which there are reliable and efficient implementations.

This situation is less than ideal, a language interface should map to platform specific properties and draw the best performance from any given hardware. In this work we present a specific algorithm for the generic lock that relies on a specific Linux utility, the futex system call.

A futex combines one atomic integer and OS scheduling. In our approach, we use one unsigned to implement the lock and a waiter count at the same time. The resulting data type of minimal size (32 bit on all architectures) and the algorithm can take advantage by minimizing the number of CPU to memory transfers. In most cases one such transfer is sufficient, where other algorithms have to update a lock and a waiter counter separately.

To our knowledge pursuing this approach to a complete solution is new. Previously urban myth had it that such approaches would risk deadlocks if confronted to heavy load, because repeated failures of calls to futex wait could lead to unbounded loops. We are able to prove that such unbounded loops will not happen for our algorithm. Also, our measurements have shown that such an approach can be very effective: failed system calls to futex wait are much less costly than commonly thought.

Our algorithm and its implementation is part of a larger open source project to provide the necessary interfaces (header files) and library support for C11's <stdatomic.h>.

It is available at http://stdatomic.gforge.inria.fr/. The benchmarks that are presented in this paper come with the reference implementation of Modular C, see http://cmod. gforge.inria.fr/. The code is functional to be used with gcc and clang, even for older version without full support for atomic operations. In a later stage, we intent to integrate the whole project into the musl C library.

TOOLS FOR DATA CONSISTENCY AND RACES

Data races are the most difficult challenge for parallel programming. They often lead to hard to trace erratic errors. The main problems are: atomicity:. Writes to memory may or may be not split by the hardware into several chunks that are written separately. A reader may occasionally see inconsistent values, for example a high and a low word of a wide integer that originate from different writes. divisibility:. A read-modify-write operations such as a simple i++ may be split into several CPU instructions, e.g a read, an increment and a write. If a first write of another threads falls sometime between the read and the write, the new write will erase the value with one that is outdated. memory ordering:. The results of read and writes to different objects in memory might be perceived by different threads in different order. Reordering of instructions may originate from optimization (compiler) or occur at run time (processor).

A data type that is guaranteed to avoid the first problem is called atomic. Before the C11 revision C only had one data type, sig atomic t, that had this guaranty to be read and written in one chunk. It can be used to communicate state between a user program and a signal handler.

For the second problem, divisibility, C had no standard tool. No other operation than read or write of sig atomic t was guaranteed to be indivisible. The indivisible operations that most computing hardware offer could not be accessed through language features. Usually they had to be programmed through extensions such as inline assembler or special compiler builtins.

Before C11, C also had no thread concept, so the memory ordering problem could not even be formalized within the vocabulary of the C standard. Obviously, it also could not provide tools to deal with it.

C atomics and its library interfaces

With modern multi-processor and multi-core hardware, parallel programming is an imperative for many if not most user applications that are used on a larger scale. Therefore it was crucial for C to provide the concepts and tools that are necessary to deal with it. So C11 introduces a lot of vocabulary and two optional features: threads through the <threads.h> interface, and atomics trough the <stdatomic.h> interface. Evidently here we are more interested in the latter, but it is important to note that both features need each other to unfold all of their potential.

C11 introduced a new qualifier, Atomic. A such qualified object guarantees that any read or write access to it is atomic in the sense we have defined above. This qualification also guarantees that between different threads all standard operations (defined through operators such as += or functional such as atomic exchange) are perceived as indivisible. Note well that this guarantee is only given between threads and in perception: in reality an operation can well be divided into several processor instructions and the perception guarantee does not extend to visibility between the main program and signal handlers. An operation that extends perception of indivisibility to signal handlers is called lock-free in the jargon of the C standard. Below we will see where this choice of words originates.

C11 also introduces different concepts of memory order for atomic operations. The whole of that specification is much too complex to unfold, here. In the following we will assume sequential consistency (memory order seq cst) for all atomic operations. This forces all atomic operations to appear totally ordered between all threads.

We will use the following atomic operations:

Inria atomic store: store a new value atomic exchange: store a new value and return previous atomic fetch add: add value and return previous atomic fetch sub: subtract value and return previous atomic compare exchange weak: compare to desired, then exchange, may fail

Atomic instructions on modern hardware

Almost since the beginning of modern computing, parallelism was implemented in hardware and the consistency problems that we introduced above became apparent.

An early overview of the problem had been given by [START_REF] Robert | What Are Race Conditions? Some Issues and Formalizations[END_REF]. Modern hardware (which almost always is inherently parallel) deals with this by providing special instructions, usually referred to as atomic instructions. It is important to note that these are not the same as the atomic operations on the level of the C language.

word size:. Usually atomic instructions are limited to word sized data types. Available on most platforms are instructions for 8, 16, 32 and 64 bit data types. Some also extend to 128 bit. primitives:. The instructions that are implemented may or may not directly correspond to atomic operations. E.g some CPU may have a proper instruction for the increment operator ++, e.g x86 64, on others, e.g arm, such an operation will have to be composed from primitives. boundedness:. Atomic instructions may give a guarantee to succeed within a bounded time (usually some memory cycles) or just return success or failure. For the latter, this may result in C level operations that have an unbounded response time. state:. Atomic instructions may operate on an internal state of the platform. E.g arm CPU work with a feature called monitors that memorize state of previous atomic access to memory.

Due to all these differences, programming with atomic instructions directly on assembly level is a mess, and in general it is very tedious to provide portability between different architectures. Such code has to ensure -the correct composition of atomic primitives to obtain sensible semantics of the operations, -the correct alignment of all atomic object such that they do not cross cache line boundaries, -the correct ordering of instructions, e.g it has to ensure that neighboring store instructions are not reordered by the CPU, -that the unboundedness of some operation may not result in application deadlocks, -and that the OS correctly restores processor state when the execution context switches from one thread to another or to a signal handler.

Luckily, C11 now ensures that only C compiler and C library Implementers have to consider all the glorious details of a specific architecture. One problem remains though, and this is what this paper is all about. Because of the limited word size for atomic instructions, the implemented compiler operations cannot just resort to a composition of atomic primitives on the atomic object itself. If an object is large, say 128 bit wide, or has a size that is not a power of 2, they must rely on external or internal locks that protect a critical section, CS. That is they need an auxiliary object that protects the data object by means of some lock primitives and by memorizing a state of the application.

Typically such locks can be made invisible between different threads, but remain visible between a thread and its signal handler. So the access to an object that is qualified with Atomic but that needs a lock for operation may be divisible with respect to a signal handler. This property is what coined C's terminology of lock-free that we already mentioned above.

Fast user space mutexes

In a singular toolbox called Fast User space muTEXes, futex for short, the Linux kernel combines two levels of operations for the implementation of lock primitives, see [START_REF] Hutton | Fuss, Futexes and Furwocks: Fast Userlevel Locking in Linux[END_REF][START_REF] Hart | A futex overview and update[END_REF]]:

(1) User space atomic integers with lock-free operations are used to regulate access to the lock as long as it is not congested. (2) Wait and wake-up system calls resolve conflicts when the lock is under congestion by multiple threads or processes. They relate to such integers by address (user space or kernel space addresses) and are guaranteed to be perceived as indivisible by the caller.

In the beginning, when futex were first introduced they needed non-standard features: assembly extensions for the atomic instructions, and a system call interface into the Linux kernel. Fortunately with the atomics interface of C11 we now have a standardized tool for the first. For the second, in the following we will assume that we have two library calls futex wait and futex wake at our disposal. This is only a very restricted set of the available futex operations that is simple to master and avoids most of the complexity of the interface as it was described by [START_REF] Drepper | Futexes are tricky[END_REF]]. The only difficulty that remains is that a call to futex wait is only successful if the kernel can guarantee that the value of the atomic object still is as expected. With these interfaces a simple but inefficient lock structure smpl could look as follows: Here the second parameter to futex wait guarantees that the thread will only be set to sleep if the value of the atomic object *lck still is prev. As a consequence the smpl lock function will iterate until the atomic exchange succeeds in modifying the value from a previous value of 0 to the value of 1.

The second parameter of futex wake corresponds to the maximal number of threads that are to be woken up. So here, the thread that holds the lock restores the object *lck to the value 0 and wakes up one possible waiter.

Both functions as described above are simplistic and not very efficient. The first, smpl lock, is inefficient because each failed attempt to acquire the lock will result in a call into the OS kernel, even if the lock would be available almost instantly. The second, smpl unlock, tries to wake up another thread without any knowledge if there even is such a thread that is waiting for it.

To avoid these two shortcomings, system libraries that implement locks (such as e.g glibc and musl) usually combine two strategies: Inria -A first spinning phase attempts the atomic operation several times. Thereby an application with a very short CS can mostly avoid sending threads into sleep. -They use at least two Atomic objects, one for the lock itself and a second one that counts the waiters. By checking if the counter is 0, this allows to avoid useless calls to futex wake.

Even though these additions enlarge the lock data structure and add one atomic operation to the unlock function these strategies have proven to be much more efficient then our simplistic versions, above.

A NEW GENERIC LOCK ALGORITHM USING FUTEX SYSTEM CALLS

To construct and adapted lock data structure for our situation, we want to have the following properties:

-The size of the data structure should be minimal: it should use just one 32 bit word as is needed by the futex calls. -When there is no contention, the number of atomic operations should be minimal.

That is one such operation for each, lock and unlock, should suffice in that case. -The procedure should be efficient, that is it should not unnecessarily waste resources. In particular, threads that have no chance to acquire the lock should be put into an OS sleep state. -The procedure should be optimized for the use case, namely a CS that is very short, just the time to a copy operation of a small object. -If the number of threads is bounded, the procedure should be deadlock free.

Typical OS lock structures also have the problem of fairness, that is that they should ensure that no single thread has an waiting time for the lock that is much longer than for others. For a lock that is used for atomics, this is only a secondary concern if we can guarantee liveness. As we will see below in the benchmark section, congestion here only is noticeable for an extremely high load, for which we may not expect any reasonable application responsiveness.

The algorithm

For our strategy, see Listing 1, we use a single unsigned value that at the same time holds the lock bit (HO bit) and a 31 bit counter.1 That counter is not viewed as a counter of the threads that are in a kernel wait, but counts the number of threads inside the critical section. An update of the counter part is done once when a thread enters the CS. Compared to the number of times the counter is accessed under congestion such events are relatively rare. Thereby we save memory bandwidth for the update, and we also avoid too much interaction between the different threads that compete for the lock.

(1) A thread is on the fast path for the lock when the overall value is 0. The lock can be acquired with one atomic operation which sets the counter and the lock bit simultaneously, if successful. (2) Otherwise, we increment the lock value atomically and enter an acquisition loop.

(a) First, we spin E times (E determined below) to set the HO bit as well, and thus acquire the lock. (b) If that fails, we suppose that the lock is under congestion and we go into a futex wait. Going into the futex wait may fail if the value changes, but since additional threads only change the counter when they arrive, this will not happen too often and the thread goes to sleep, eventually. Unlocking is a very simple operation. The locker has contributed ftx set(1u) to the value, and just has to decrement the value atomically by that amount. The return value of the operation reveals if other threads still are in the CS, and a futex wake call can be placed accordingly.

void ftx unlock ( f t x * lck ) { unsigned prev = atomic fetch sub ( lck , f t x s e t (1 u ) ) ; i f ( prev != f t x s e t (1 u ) ) futex wake( lck , 1) ; }

Analysis

It is relatively easy to see that this new strategy provides a functional lock primitive using just a 32 bit data structure and one atomic operation for fast ftx lock and ftx unlock. It remains to show that it cannot deadlock. The worst case scenario for our use of our lock primitive is that the thread that holds the lock, say T 0 , is unscheduled while inside the CS. Suppose further that there are N other threads that are ready to be scheduled, and that once they are scheduled they start to compete for the lock.

Different quantities are interesting for an analysis of the runtime behavior of the algorithm. We can control one of them, namely the time t mono that a scheduled thread spends spinning before trying to switch to futex wait. Three others are platform dependent: t fail . is the maximum of two system specific times: the time a thread T 1 may either spend in a failed attempt to futex wait or that the system needs to put T 1 to sleep and start another thread T 2 . Inria P . is the number of processor cores, which is viewed to be equal to the maximum number of threads that are scheduled simultaneously. t para . is the time that P threads need for a spinning phase that they perform in parallel.

A value t para close to t mono indicates a perfect parallelism, a value of P ⋅ t mono means that there is none at all. Usually it will be greater than t mono , e.g because of memory contention or contention on other shared resources (execution pipelines, caches). We derive some other quantities from the above: P . given as P ⋅tmono tpara is the parallelism of the platform. E. given as tmono tpara = P P is the efficiency of the platform.

For example, on a modern hyperthreaded machine with 4 cores in total, P is typically between 2.5 and 3, E is between 0.625 and 0.75.

Remark 3.1. On a platform where P is close to one, the spinning phase of the algorithm should entirely be skipped. This is simply because there no other thread can make progress while a thread is spinning. Thus spinning would just waste resources and the state of the application would not progress. So from now on we may assume that P ≥ 1 + for some reasonable value of > 0.

LEMMA 3.2. that no other threads are unscheduled, after at most t para + (P -1) ⋅ t fail seconds a first thread successfully calls futex wait.

PROOF. For the first term, observe that after t para time, at least one thread has finished the spinning phase, and attempts futex wait.

While no thread is unscheduled at most P scheduled threads can enter the CS. There are at most P -1 atomic increments that change the futex value. Thus the first thread that enters the CS will need at most t para time for spinning and then futex wait may fail at most P -1 times in a row. This already shows that, provided no other descheduling takes place, our algorithm is deadlock-free. Now, once a thread successfully goes into futex wait a new thread T P can be scheduled, compete for the lock and change the futex value. It may disturb all other threads that are trying to go into futex wait, forcing them to restart their attempt.

Remark 3.3. Provided that no threads are unscheduled otherwise, that there are always P threads inside the CS and that at least one of them has finished spinning, after a time of t fail another threads succeeds his call to futex wait.

That is, under these circumstances we have a stable regime where each t f ail seconds a thread enters futex wait.

To be able to ensure that there is always at least one thread that has finished spinning, we observe that if t para ≤ t fail (or equivalently t mono ≤ E ⋅ t fail ) a newly scheduled thread has finished spinning when the next thread successfully goes into futex wait.

LEMMA 3.4. Provided that no threads are unscheduled otherwise, that there are always P threads inside the CS and that t para ≤ t fail , threads succeed calls to futex wait at a rate of 1 t fail per second.

Or, roughly if P ≪ N the time for all threads to calm down and successfully call futex wait is N ⋅ t fail . THEOREM 3.5. Let be T 0 a thread out of N ≫ P that is unscheduled when holding the lock. Provided that none of the threads is unscheduled by other means and that t para ≤ t fail , after a time of N ⋅ t fail the application makes progress.

PROOF. This progress can be of two forms. Either there is another thread than T 0 that does not enter the the CS and thus progresses the application, or T 0 will be rescheduled and finishes its CS.

The time t mono has not only an influence for this worst case, but is also responsible for the response time in the non-congested situation. The longer we spin, the higher the probability to acquire the lock without going into futex wait. So the best compromise would be to choose t mono = E ⋅ t fail . Practically a factor of 0.9 always guarantees liveness of the application and shows good performance on average.

BENCHMARKS

The framework

We have run a long series of benchmarks to validate the approach. The code for the benchmark is integrated in p11 with comes with Modular C, see Cmod. For compilation of that benchmark we also need a C11 compliant library, that has C11 threads, and a C11 compiler that also has gcc extension. We used musl for the first and gcc and clang for the latter.

The implementation of our algorithm is a bit more sophisticated than what may appear above. In particular it takes care of reducing the number of atomic operations to a minimum and to use memory ordering for the locks that is adapted to the case.

The test program

The test in p11 is called p11#test#lifo. It is based on a stack implementation (Last In First Out) that uses an atomic pair of pointers for the head to avoid the ABA problem, [IBM 1983;[START_REF] Maged | ABA Prevention Using Single-Word Instructions[END_REF]. For the benchmarks, the size of the atomic data structure has been chosen in such a way that the generic atomic functions based on locks are chosen.

The idea of this benchmark is to have a application that runs under full load, stresses the platform with a lot of allocations and deallocations and in the middle of that does a lot of locking and unlocking. It works as follows:

It creates or deletes a random number of list elements for the lifo inside a loop. All test runs last 10s and were repeated at least 10 times. The measure that is reported is the number of list elements that have been handled per second on average.

The parameters of the runs are the number of threads that run in parallel, values ranging from 1 up to 256. Different lock primitives can be chosen at compile time to protect the head of the LIFO: futex. the futex based algorithm described here mutex. based on a standard mutex musl. musl's low-level lock/unlock functions spin. a spin lock using atomic exchange native. the compilers "native" generic lock, also a mutex The compiler is gcc version 5.2. Here, the native C library is musl and so the atomic library of gcc are compatible. Therefore the benchmarks include "native".

4.3.2. A x86 64 machine with 2x2 hyperthreaded cores. This is a i7-4600U CPU at 2.10GHz and with 8 GiB of RAM. The OS is Debian Linux, with glibc as native library. The processor has atomic instructions for word sizes up to 128 bit. The compiler is gcc version 4.9. Because musl is not the native C library, the atomic library of gcc is not compatible. Therefore the benchmark "native" is missing.

Performance comparison

Fig. 1(a) shows the results on the arm platform. We see that all lock implementations allow for an acceleration of the application when a small number of threads is used. But what is also clear that the "native" lock performs worst for the case that is the most interesting: the range where each thread has its own CPU core at its disposal. Even the "mutex" lock performs better.

We also see that musl's internal lock structure shows a drastic performance loss when it comes to congestion. This is due to a switch of the spinning strategy: as soon as congestion is detected, spinning is abandoned and threads directly attempt futex wait. This is meant to ensure fairness of lock acquisition, but as we can see for our use case it has a dramatic impact on the application throughput.

Fig. 1(b) shows the relative performance of the same experiments, where the "mutex" implementation is taken as a base. We see that our new implementation is about 60% better than the "native" version, or 40% than a direct implementation with mutex. It combines the good performance of a spinlock for the less congested range with a good policy for strong congestion.

To finish let us consider the x86 64 platform, Fig. 2. Although it has more compute power than the other, the atomics of the hardware are much less performing. This is due to the fact that here an atomic instruction usually enforces a complete synchronization at a cost of about 50 CPU cycles. Basically, the CPU is blocked for this number of cycles. Compared to that, in a monitor based approach as on the arm architecture part of these cycles can be used for other computations. So on the x86 64 platform any atomic operation incurs a strong latency penalty. Thereby, our application is not even able to accelerate for 2, 3 or 4 threads as it was the case on arm. In the contrary it even decelerates. Nevertheless the relative performance difference between the different lock implementations look very similar.

CONCLUSION

We have presented a new locking algorithm that combines consequent use of C11 atomics with Linux' futex system calls. We have proven that it is deadlock free. is deadlock free, and that it shows better performance than other lock implementations. This is not surprising, an implementation that is tuned for the purpose (very short CS) and that may avoid stacked calls into the C library should always perform better than a generic one. Surprising to us was the wide performance gap between the implementations.

By pursuing this research we learned to mistrust some of the urban legends that turn around atomics, futexes and lock structures in general. At least when we stick to the basics (futex wait and futex wake) and if we have a decent interface for atomics, programming them is not as difficult as the legends suggest. Also using a system call is not so much worse that spinning around an atomic access. The performance factor between the two is only about 10, and so spinlocks in the order of 10 should be sufficient in many cases.

This support library is now available as open source at http://stdatomic.gforge.inria. fr. We hope to integrate it into the C library that we used for most of our experiments, musl.

  Listing 1. A lock algorithm using a single Atomic unsigned. typedef Atomic(unsigned) f t x ; # define f t x s e t ( VAL ) (0 x80000000u | ( VAL ) ) # define f t x l k d ( VAL ) (0 x80000000u & ( VAL ) ) void f t x l o c k ( f t x * lck ) { unsigned cur = 0; i f (! atomic compare exchange weak( lck , & cur , f t x s e t (1) ) ) { cur = f t x f e t c h a d d ( lck , 1) + 1; for (;;) { while (! f t x l k d ( cur ) ) { i f (atomic compare exchange weak( lck , & cur , f t x s e t ( cur ) ) ) return ; for (unsigned i = 0; i < E && f t x l k d ( cur ) ; i ++) cur = atomic load ( lck ) ; } while ( f t x l k d ( cur ) ) { futex wait ( lck , cur ) ; cur = atomic load ( lck ) ; } } } }

4. 3 .

 3 Fig. 1. benchmarks on arm

  Fig. 2. benchmarks on x86 64

RR n°8818

On Linux, unsigned is always 32 bit wide.RR n°8818


Bastien Mussard
email: bastien.mussard@upmc.fr

János G Ángyán
email: janos.angyan@univ-lorraine.fr

Prof Péter Surján

Bastien Mussard

Local Random Phase Approximation with Projected Oscillator Orbitals

Keywords: RPA, oscillator orbitals, London dispersion energy, dispersion coefficient, local correlation method 1 Introduction

An approximation to the many-body London dispersion energy in molecular systems is expressed as a functional of the occupied orbitals only. The method is based on the local-RPA theory. The occupied orbitals are localized molecular orbitals and the virtual space is described by projected oscillator orbitals, i.e. functions obtained by multiplying occupied localized orbitals with solid spherical harmonic polynomials having their origin at the orbital centroids. Since we are interested in the long-range part of the correlation energy, responsible for dispersion forces, the electron repulsion is approximated by its multipolar expansion. This procedure leads to a fully non-empirical longrange correlation energy expression. Molecular dispersion coefficients calculated from determinant wave functions obtained by a range-separated hybrid method reproduce experimental values with less than 15% error.

approximation (RPA) is situated on the highest, fifth rung of Jacob's ladder, which leads from the simplest Hartree level towards the "heaven" corresponding to the exact solution of the Schrödinger equation. When stepping upwards on Jacob's ladder, one uses more and more ingredients of the Kohn-Sham single determinant. Starting from the lowest rungs, the density, its gradient, the full set of occupied orbitals are successively necessary for the construction of the functional. At the highest rung DFA is usually based on many-body methods, which require the knowledge of the complete set of occupied and virtual orbitals. Such methods have the drawback that the size of the virtual orbital space can be very large even in a moderately sized atomic orbital basis. In the case of plane wave calculations the virtual space can become even prohibitively large. One solution to keep the size of matrices in reasonable limits makes recourse to an auxiliary basis set to expand the occupied-virtual product functions. Such approaches are known in quantum chemistry as resolution of identity [START_REF] Eichkorn | [END_REF] or density-fitting [3,4] methods. Similar advantages can be achieved by Cholesky decomposition [5] techniques. In plane wave calculations the plane wave basis itself can be used to expand the product states [6]. Further gain can be achieved by projection methods, which avoid any explicit reference to virtual orbitals. Such a technique, named projective dielectric eigenpotential (PDEP) method, has been successfully applied to construct the dielectric function in plane wave RPA calculations [7][8][9]. Recently, Rocca succeeded to reduce further the size of the problem by a prescreened set of virtuals [10].

The purpose of the present work is to demonstrate that an approximate variant of the RPA (and also of the MP2) correlation energy can be expressed by using quantities that are computable from occupied orbitals alone. In this context, one should mention the beautiful result, which has been obtained by Surján, who has shown that the MP2 correlation energy, exempt of any further approximation in a given basis set, can be reformulated as a functional of the Hartree-Fock density matrix, i.e. using exclusively the occupied orbitals [11]. Our main interest is not to reproduce the full correlation energy with high numerical precision, but we focus our attention to a well-defined part of it, namely on the longrange dynamical correlation energy, which is usually taken responsible for the London dispersion forces.

It is now well-documented that most of the conventional density functional calculations in the Kohn-Sham framework, unless special corrections are added to the total energy, are unable to grasp the physics of these longrange forces. Various fairly successful dispersion correction schemes are known, but although most of them were claimed to posses an essentially ab initio character, they use, without exception, some "external" data, like atomic polarizabilities, atomic radii, etc. . . [12][13][14][15][16][17][18]. Even if these quantities are usually not taken from experimental sources, and rather originate from ab initio computations, their presence deprives these theories of their self-contained character. Hence, although we do not pretend that the rather drastic approximations which are going to be implemented in the following allow us to achieve results of a quality comparable to the precision attained by carefully fine-tuned methodologies, we argue that the ingredients of the present approach originate from a controlled series of approximations and do not use any "external" inputs. Moreover, in contrast to the relatively costly and sophisticated methods, like RPA, we do not need virtual orbitals, i.e. we stay on the fourth rung of Jacob's ladder.

Our approach relies on the use of localized occupied molecular orbitals (LMOs) that can be obtained relatively easily by a unitary transformation in the subspace of occupied orbitals, according to either an external or an internal localization criterion [START_REF] Surján | Theoretical Models of Chemical Bonding[END_REF]. The localization of virtual orbitals is much more difficult, since the usual localization criteria for the occupied orbitals lead often to divergent results. It is to be noted that recently a significant progress has been reported [START_REF] Høyvik | [END_REF] for the efficient localization of virtual orbitals. However, we follow another strategy here and build excited determinants using localized functions which are able to span the essential part of the virtual space. One of the most popular local correlation approaches in this spirit consists in using the atomic orbital (AO) basis functions to represent the virtual space. They are made orthogonal to the occupied space by projection, leading to the projected AOs (PAOs) techniques [21][22][23]. The locality of these functions is guaranteed by construction, even if it may be somewhat deteriorated by the projection procedure.

In the present work we are going to revisit and explore a quite old idea of Foster and Boys from the early sixties [24][25][START_REF] Boys | Quantum Theory of Atoms, Molecules, and the Solid State, A Tribute to[END_REF]. The main concern of these authors was to construct a set of virtual orbitals directly from the set of occupied LMOs by multiplying them with solid spherical harmonic functions centered on the barycenter of the LMO. The orthogonality of these new functions can be ensured by a projection procedure. Boys and Foster called these new functions, obtained after multiplication, oscillator orbitals (OOs), and after removing the components of the OOs in the space of the occupied orbitals they may be called projected oscillator orbitals (POOs). Very few articles in the literature mention Boys' oscillator orbitals [START_REF] Alagona | [END_REF], probably because it had no particular numerical advantage in high-precision configuration interaction calculations and its practical implementation raised a number of complications which could be avoided by more straightforward algorithms, like the use of the full set of virtual molecular orbitals (VMOs). We have found only a single, very recent article, which referred to the notion of oscillator orbitals [28] as a useful concept, but not as a practical computational tool. To the best of our knowledge, the mathematical implications of using oscillator orbitals to define the virtual space has never been rigorously studied. Such an analysis is beyond the scope of the present study: it is going to be the subject of a forthcoming publication.

The POOs are non-orthogonal among each other, which is at the origin one of the complications mentioned above. This problem can be handled just like in the case of the PAOs, therefore a theory of electron correlation based on POOs can follow a similar reasoning as local correlation methods using PAOs [START_REF] Knowles | Modern Methods and Algorithms of Quantum Chemistry[END_REF]. In particular, in this paper, the RPA method will be reformulated for a virtual solace constructed from POOs.

It is worthwhile to mention that the projected oscillator orbitals bear some similarities to the trial perturbed wave function in the variation-perturbational technique associated with the names of Kirkwood [START_REF] Kirkwood | [END_REF], Pople and Schofield [31] (KPS), to calculate multipole molecular polarizabilities. The closely related Karplus-Kolker [32,33] (KK) method and its variants [34,35] use a similar Ansatz for the perturbed orbitals. In these latter methods, which were formulated originally as simplified perturbed Hartree-Fock theories, the first order perturbed wave function is a determinant with first order orbitals ψ (1) i which are taken in the following product form [36]

ψ (1) i = g i ψ i -k ψ k | g i ψ i ψ k ,
where g i are linear combination of some analytically defined functions, like polynomials. As we shall see, the principal difference of this Ansatz and the POOs is that in the former case one multiplies the occupied canonical orbitals with the function g i , while the oscillator orbitals are constructed from localized orbitals.

As mentioned previously, our main focus is the modeling of London dispersion forces. It has been demonstrated in our earlier works [37][38][39] that the essential physical ingredients of London dispersion forces are contained in the range-separated hybrid RPA method, where the shortrange correlation effects are described within a DFA and the long-range exchange and correlation are handled at the long-range Hartree-Fock and long-range RPA levels, respectively. Among numerous possible formulations of the RPA [39,40], we have chosen to adopt the variant based on the ring-diagram approximation to the coupled cluster doubles theory. The relevant amplitude equations will be rewritten with the help of POOs leading to simplified working equations, which do not refer to virtual orbitals explicitly. For the sake of comparison we are going to study the case where the POOs are expanded in the virtual orbital space. The long-range electron repulsion integrals, appearing in the range-separated correlation energy expression, can be reasonably approximated by a truncated multipole expansion. It means that in addition to the well-known improved convergence properties of the correlation energy with respect to the size of the basis set, one is able to control the convergence through the selection of the multipolar nature of the excitations, leading to a possibility of further computational gain.

The exploitation of localized orbitals for dispersion energy calculations has already been proposed since the early works on local correlation methods [41][42][43][44][45]. In classical and semiclassical models most often the atoms are selected as force centers; only a few works exploit the advantages related to the use of two-center localized orbitals and lone pairs. A notable exception is the recent work of Silvestrelli and coworkers [46][47][48][49][50], who adapted the Tkatchenko-Scheffler model [16] for maximally localized Wannier functions, which are essentially Boys' localized orbitals for solids. It is worthwhile to mention that one of the very first use of the bond polarizabilities as interacting units for the decryption of London dispersion forces has been suggested as early as in 1969 by Claverie and Rein [51]; see also [START_REF] Claverie | Intermolecular Interactions: From Diatomics to Biopolymers[END_REF].

Our approach, at least in its simplest form, is situated somewhere between classical models and fully quantum local correlation methods and can be considered (practically in all its forms) as a coarse-grained nonlocal dispersion functional formulated exclusively on the basis of ground state densities and occupied orbitals. It will be shown how the various matrix elements can be expressed from occupied orbital quantities only. As a numerical illustration, molecular C 6 dispersion coefficients will be calculated from localized orbital contributions and compared to experimental reference data. The paper will be closed by a discussion of possible future developments.

Theory

Projected oscillator orbitals

The a posteriori localization of the subspace of the occupied orbitals is a relatively standard procedure, which can be achieved following a large variety of localization criteria (for a succinct overview, see Ref. [START_REF] Alcoba | [END_REF]). In the context of correlation energy calculations, i.e. in various "local correlation approaches", the most widely used localization methods are based either on the criterion of Foster and Boys [24] or that proposed by Pipek and Mezey [54]. For reasons which become clearer below, in the present work we will use the Foster-Boys localization criterion, which can be expressed in various equivalent forms [START_REF] Boys | Quantum Theory of Atoms, Molecules, and the Solid State, A Tribute to[END_REF]. In its the most suggestive formulation, the Foster-Boys' localization procedure consists in the maximization of the squared distance between the centroids of the orbitals: max

         occ i< j | φ i | r |φ i -φ j | r |φ j | 2          . ( 1
)
Another form of the localization criterion, which is strictly equivalent to the previous one, corresponds to the minimization of the sum of quadratic orbital spreads min

       occ i φ i | r2 |φ i -| φ i | r |φ i | 2        . (2)
As it has been demonstrated by Resta [55], the previous minimization implies that the sum of the spherically averaged squared off-diagonal matrix elements of the position operator is minimal too. This last property of the Boys' localized orbitals is going to be useful in the development of the present model. Any set of localized orbitals obtained by a unitary transformation from a set of occupied orbitals spans the same invariant subspace as the generalized Kohn-Sham operator, f µ , and satisfies the equation:

f µ φ µ i = occ j ε µ i j φ µ j , (3)
where f µ φ i = t + vne + vH + vlr,µ x,HF + vsr,µ xc,DFA is the rangeseparated hybrid operator. In this expression t is the kinetic energy, vne is the nuclear attraction, vH is the full-range Hartree, vlr,µ

x,HF is the long-range Hartree-Fock (non-local) exchange, and vsr,µ xc,DFA is the short-range exchange-correlation potential operator. The range-separation parameter, µ, defined below, cuts the electron repulsion terms to short-and long-range components. For µ = 0 one recovers the fullrange density functional approximation (DFA), while for µ → ∞ one obtains the full-range Hartree-Fock theory.

As outlined in the introduction, inspired by the original idea suggested first by Foster and Boys [24] and refined later by Boys [START_REF] Boys | Quantum Theory of Atoms, Molecules, and the Solid State, A Tribute to[END_REF], we propose here to construct localized virtual orbitals by multiplying the localized occupied orbital φ i (r) by solid spherical harmonics having their origin at the barycenter of the localized occupied orbital. According to Boys [START_REF] Boys | Quantum Theory of Atoms, Molecules, and the Solid State, A Tribute to[END_REF], this definition can be made independent of the orientation of the coordinate system by choosing local coordinate axes which are parallel to the principal axes of the tensor of the moment of inertia of the charge distribution φ * i (r)φ i (r) (see Appendix A).

In the following we are going to elaborate the theory for the simplest case, when these oscillator orbitals are generated by first order solid spherical harmonic polynomials, i.e. the i-th LMO is multiplied by (r α -D i α ), where rα is the α = x, y, z component of the position operator and D i α is a component of the position vector pointing to the centroid of the i-th LMO, defined as

D i α = φ i | rα |φ i .
It is possible to generate oscillator orbitals by higher order spherical harmonics too, which is left for forthcoming work. We denote the POO by | φi α , where the index i α refers to the fact that the OO has been generated from the i-th LMO by using the rα function. For the sake of the simplicity of the formulae, the POOs will be expressed in the laboratory frame; the expressions for the dipolar POOs in a local frame are shown in the Appendix A. The POO reads in the laboratory frame as

| φi α = Î - occ m |φ m φ m | rα -D i α |φ i = Q rα |φ i , (4)
with Q = ( Îocc m |φ m φ m | ), the projector onto the virtual space.

In Eq. 4 the "pure" (OO) component of the POO is rα |φ i , while the orthogonalization tails stem from the term occ m i |φ m φ m | rα |φ i . In this sense the "locality" of the OO is somewhat deteriorated, since we have contributions from each of the other LMOs. At this point the Boys localization criterion will be at our advantage, since it ensures that the sum of the off-diagonal elements of the x, y and z operators taken between the occupied orbitals, and appearing in the orthogonalization tails, be minimized [55]. In this sense, the Boys-localization scheme seems to be naturally adapted for the construction of dipolar oscillator orbitals.

To illustrate the concept of dipolar oscillator orbitals, two examples are taken from the oxygen lone pair and the C-H bonding orbitals of the formaldehyde molecule, described in the above-defined local frame. Figs 1 and2 show the localized orbitals and the three projected dipolar oscillator orbitals having a nodal surface intersecting the orbital centroid and oriented in the three Cartesian coordinate directions of the local coordinate system, x, y and z. It is quite clear that the node coincides with the region of the highest electron density of the orbital and in this sense ensures an optimal description of the correlation. Higher order polynomials generate virtual orbitals with further nodes.

From now on, we will use the simplified notations |φ i ≡ |i and | φi α ≡ |i α , in other words the subscript α on the orbital index indicates that it is an oscillator orbital. We designate the occupied (canonical or localized) molecular orbitals as i, j, k, . . . and the canonical virtual molecular orbitals (VMOs) as a, b, c, . . . . The oscillator orbitals are non-orthogonal among each other; their overlap integral can be evaluated using the idempotency of the projectors:

S i α , j β = i α | j β = i| rα rβ | j - occ m i| rα |m m| rβ | j . (5)
Note that the overlap matrix S is of size N POO × N POO , where N POO is the number of projected oscillator orbitals. The POOs can be expanded in terms of a set of orthonormalized virtual orbitals (e.g. the set of canonical virtuals). Although later we eliminate explicit reference to the set of virtual orbitals of the Fock/Kohn-Sham operator (i.e. everything will be written only in terms of the occupied orbitals or equivalently in terms of the corresponding density matrix), with the help of the resolution of identity, we give the explicit form of the coefficient matrix linking the POOs with the virtuals:

|i α = all p |p p| Q rα |i = virt a |a a| rα |i = virt a |a V a i α . (6)
The matrix V is constructed simply from the elements of the occupied/virtual block of the position operator, V a i α = a| rα |i . The overlap matrix of the expanded POOs can be written in terms of the coefficient matrix V:

S i α , j β = i α | j β = virt ab V † i α a a | b V b j β = (V † V) i α j β . (7)
Higher order oscillator orbitals, not used in the present work, can be generated in an analogous manner, using higher order solid spherical harmonic functions.

Ring CCD-RPA equations with POOs

In the ring CCD (ring coupled cluster double excitations) formulation [39], the general RPA correlation energy (direct-RPA or RPA-exchange) is a sum of paircontributions attributed to a pair of occupied (localized) orbitals: where the amplitudes T i j satisfy the Riccati equations, which can be written in terms of orthogonalized occupied i, j, k, . . . and virtual a, b, c, . . . orbitals as:

E RPA c = 1 2 occ i j tr B i j T i j , (8)
R i j = B i j + (( + A)T) i j + (T(A + )) i j + (TBT) i j = 0, (9)
with the matrix elements:

i j ab = δ i j f ab -f i j δ ab A i j ab = K i j ab -J i j ab = a j | ib -a j | bi (10)
B i j ab = K i j ab -K i j ab = i j | ab -i j | ba ,
where f is the fock matrix of the fock operator f and with the two-electron integrals of spin-orbitals written with the physicists' notation:

i j | ab = φ * i (r)φ * j (r )w(r, r )φ a (r)φ b (r ) drdr , (11)
where w(r, r ) = |r -r| -1 is the Coulomb electron repulsion interaction.

Note that the size of the matrices in Eq. 8 and 9 is, for each pair [i j], N virt × N virt and that in this notation the matrix multiplications are understood as, for example: (T(A + )) i j ab = mc T im ac A m j cb + T im ac m j cb . Using the transformation rule between the amplitudes in the VMOs and in the POO basis, T i j = VT i j POO V † (see Appendix B), the Riccati equations can be recast in the POO basis as (again, see Appendix B):

R i j POO = B i j POO + ( POO + A im POO ) T m j POO S POO + S POO T im POO ( POO + A m j POO ) + S POO T im POO B mn POO T n j POO S POO = 0. ( 12
)
This corresponds to a local formulation of the ring CCD amplitudes equations, and the dimension of the matrices, emphasized by the subscript POO, is merely N POO × N POO . As explained in Appendix C, this type of Riccati equations can be solved iteratively in a pseudo-canonical basis.

Local excitation approximation

Since the excitations are limited to "pair-domains", the effective dimension of the equations for a pair is actually roughly independent from the size of the system, just like in any local correlation procedure. As the simplest approximation, one can take only the three excitations to the dipolar POOs generated by a selected LMO, i.e. for each pair of LMOs [i j] we have the local excitations i → i α and j → j β , leading to a 3 × 3 problem to solve and iterate on. Within this approximation, all matrices involved in the derivations can be fully characterized by two occupied LMO indices and two cartesian components, α and β (the subscript POO is omitted from now on):

(A) i j i α j β ≡ (A) i j αβ and: (S) i α j β ≡ (S) i j αβ (B) i j i α j β ≡ (B) i j αβ (f) i α j β ≡ (f) i j αβ ( 13
) (T) i j i α j β ≡ (T) i j αβ
With this in mind, and with the additional approximation which consists in neglecting the overlap between POOs coming from different LMOs, i.e. (S) i j αβ ≈ δ i j (S) ii αβ , the direct RPA Riccati equations of Eq. 12 become: R i j = B i j + f ii T i j S j jf ii S ii T i j S j j + A im T m j S j j + S ii T i j f j j -S ii T i j S j j f j j + S ii T im A m j + S ii T im B mn T n j S j j = 0.

In the above equation, we have written explicitly the fock matrix contributions and used implicit summation conventions over m and n. A detailed derivation of Eq. 14 from Eq. 12 is shown in Appendix D. These Riccati equations can be solved by a transformation to the pseudo-canonical basis, as described in Appendix C.

Multipole approximation for the long-range two-electron integrals in the POO basis

In the context of range-separation, and in the spirit of constructing an approximate theory which takes advantage of the localized character of the occupied molecular orbitals, we are going to proceed via a multipole expansion of the long-range two-electron integrals.

The matrices A i j POO and B i j POO will be reinterpreted in terms of long-range two-electron integrals, i.e. w(r, r ) will be replaced by w lr (r, r ) = erf µ|r -r| |r -r| -1 in Eq. 11. They read respectively as (see Eq. 10 and the transformation described in Appendix B):

K i j m α n β -J i j m α n β = m α j | in β lr -m α j | n β i lr ( 15
)
K i j m α n β -K i j m α n β = i j | m α n β lr -i j | n β m α lr . (16)
Note that these integrals could be calculated by using the POO to VMO transformation of Eq. 6. However, such an expression is not in harmony with our goal of getting rid of virtual orbitals, since it requires the full set of integrals transformed in occupied and canonical VMOs with an additional two-index transformation. We could formally eliminate virtual molecular orbitals by applying the resolution of identity, but in this case we would be faced with new type of two-electron integrals, in addition to the usual ones generated by the Coulomb interaction |rr | -1 , namely integrals generated by rα |rr | -1 , r β |rr | -1 and rα r β |rr | -1 . Therefore we are going to proceed by a multipole expansion technique.

The expansion centre for the multipole expansion will be chosen at the centroid of the LMOs, i.e. in this example at D i and D j . Using the second order long-range interaction tensor L i j (D i j ), with D i j = D i -D j (see Appendix E), we have

K i j m α n β = m α j | in β lr ≈ γδ m α | rγ |i L i j γδ j| rδ |n β + higher multipole terms. (17)
A truly remarkable formal result emerging from the framework of oscillator orbitals is that the rγ matrix element between the POO m α and the LMO i that appears in the previous equation is nothing else but the overlap between the POOs m α and i γ (cf. Eq. 5):

m α | rγ |i = m| rα rγ |i - occ n m| rα |n n| rγ |i = S m α i γ , (18)
so that the matrix element K i j m α n β simply reads:

K i j m α n β = S m α i γ L i j γδ S j δ n β . (19)
After applying the local excitation approximation, this bielectronic integral becomes even simpler, according to the following expression:

K i j = S ii L i j S j j . ( 20
)
In the case of direct RPA, only the K i j two-electron integrals are needed. For the more general exchange RPA (RPAx) case, most of the electron repulsion integrals, m α j | n β i lr , can be neglected in the multipole approximation, since they correspond to the interaction of overlap charge densities formed by localized orbitals in different domains. Nevertheless, integrals of the type m α j | m β j lr should be kept: they describe the interaction of the overlap charge densities of the j-th LMO and the m α POO, which is a typical long-range Coulomb interaction. Similar considerations hold for the K i j matrix elements, which correspond to an exchange integral involving overlap charge densities of orbitals belonging to different domains and can be neglected at this point (a few integrals will however survive). The RPAx variant of the model will be considered in more details in forthcoming works.

Spherical average approximation

The previously discussed 3 × 3 matrices can be easily replaced by scalar quantities, if one considers a spherical average of the POO overlap and fock matrices:

S ii αβ ≈ 1 3 s i δ αβ with s i = α S ii αα , (21)
and:

f ii αβ ≈ 1 3 f i δ αβ with f i = α f ii αα . (22)
In this diagonal approximation, and in the case of direct RPA where A = B = K, the Riccati equations of Eq. 14 supposing implicit summations on m and n become: R i j = s i s j L i j + f i s jf ii s i s j T i j + 1 3 s i s m s j L im T m j + s i f js i s j f j j T i j + 1 3 s i s m s j T im L m j + 1 3 2 s i s m s n s j T im L mn T n j = 0.

This set of equations can be solved directly, i.e. without proceeding by the pseudo-canonical transformation described in Appendix C for the more general case. The only quantities needed are the spherically averaged s i and f i associated to localized orbitals and the long-range dipole-dipole tensors. The update formula to get the n-th approximation to the amplitude matrix element is

T i j (n) αβ = s i s j L i j αβ + ∆R i j αβ (T (n-1) ) ∆ i s j + s i ∆ j , (24)
with:

∆ i = f ii s i -f i , (25)
and:

∆R i j (T) = 1 3 s i s m s j L im T m j + 1 3 s i s m s j T im L m j + 1 3 2 s i s m s n s j T im L mn T n j . ( 26
)
Pursuing with the local excitation and the spherical average approximations, the long-range correlation energy is given by the following spin-adapted expression:

E RPA,lr c = 4 9 occ i j s i s j tr L i j T i j . (27)

Bond-bond C 6 coefficients

Using the first order amplitudes, i.e. the amplitudes obtained in the first iteration step during the solution of Eq. 24:

T i j (1) = s i s j ∆ i s j + s i ∆ j L i j , (28)
the second order long-range correlation energy becomes:

E (2),lr c = 4 9 occ i j s i s j s i s j ∆ i s j + s i ∆ j tr L i j L i j . ( 29
)
This expression describes the correlation energy as a pairwise additive quantity made up from bond-bond contributions. The summation over the components of the longrange interaction tensors gives (see Appendix E):

tr L i j L i j = 6 D i j6 F µ damp (D i j ), (30)
and allows us to cast the long-range correlation energy in a familiar form, as:

E (2),lr c = occ i j C i j 6 D i j 6 F µ damp (D i j ), (31)
where C i j 6 is the dispersion coefficient between the i and j LMOs:

C i j 6 = 8 3
s i s j s i s j ∆ i s j + s i ∆ j . (32)
Note that the above dispersion coefficient corresponds to a single-term approximation to the bare (non-interacting) spherically averaged dipolar dynamic polarizability associated with the localized orbital i,

α i 0 (iω) ≈ 4 3 ω i ω 2 i + ω 2 s i , (33)
where ω i = ∆ i /s i is an effective energy denominator and the quantity s i stands for the second cumulant moment (spread) of the localized orbital. Note that this possibility to approximate α i 0 (iω) only as a function of objects like f i and s i is a direct consequence of the remarkable feature that the second moment between an LMO and a POO corresponds to an overlap between two POOs (see Eq.

It is easy to verify that with the help of the Casimir-Polder formula

C i j 6 = 3 π ∞ 0 dω α i 0 (iω) α j 0 (iω), ( 34
)
that it is indeed the C i j 6 dispersion coefficient which is recovered from the polarizabilities defined in Eq. 33. This simple model for the dynamic polarizability associated to an LMO, deduced from first principles, can be the starting point of alternative dispersion energy expressions, e.g. based on the modeling of the dielectric matrix of the system or using the plasmonic energy expression.

Preliminary results: molecular C 6 coefficients

In order to have a broad idea about the appropriateness of this simple dispersion energy correction, we have calculated the molecular C 6 coefficients for a series of homodimers as the sum of the atom-atom dispersion coefficients given by Eq. 32.

The matrix elements between POOs, S ii αβ and f ii αβ , which are needed to calculate the scalars s i and f i , can be obtained directly by manipulating the matrix representation of the operators. Such a procedure leads to what we will call the "matrix algebra" expressions (denoted by [M]), of the following form:

s i [M] = α i| rα Q rα |i = virt a | i| r |a | 2 , (35)
and (see Appendix F):

f i [M] = virt ab α i| rα |a f ab b| rα |i . (36)
Alternatively, these matrix elements can be obtained through the application of commutator relationships, therefore this latter option will be referred to as the "operator algebra" approach (denoted by [O]). They take the form:

s i [O] = α i| rα Q rα |i = i| r2 |i - occ m | i| r |m | 2 , (37)
and (again, see Appendix F):

f i [O] = 3 2 + 1 2 occ m f im m| r2 |i + i| r2 |m f mi - occ mn α i| rα |m f mn n| rα |i . (38)
Four different Fock/Kohn-Sham operators have been applied to obtain the orbitals, which are subsequently localized by the standard Foster-Boys procedure. In addition to the local/semi-local functionals LDA and PBE, the rangeseparated hybrid RSHLDA [37,56] with a range-separation parameter of µ = 0.5 a.u. as well as the standard restricted Hartree-Fock (RSH) method were used. The notations LDA[M] and LDA[O] refer to the procedure applied to obtain the matrix elements: either by the matrix algebra [M] or by the operator algebra [O] method. All calculations were performed with the aug-cc-pVTZ basis set, using the MOL-PRO quantum chemical program package [57]. The matrix elements were obtained by the MATROP facility of MOL-PRO [57]; the C 6 coefficients were calculated by Mathematica.

H 2 H F H 2 O N 2 C O N H 3 C H 4 H C l C O 2 H 2 C O N 2 O C 2 H 2 H B r H 2 S C H 3 O H S O 2 C 2 H 4 C H 3 N H 2 S iH 4 C 2 H 6 C l 2 C H 3 C H O C O S C H 3 O C H 3 C 3 H 6 C S 2 C C l 4
Error in the molecular C 6 coefficients (%)

LDA[M] PBE[M] RHF[M] RSHLDA[M]
Fig. 3 Percentage errors in calculated molecular C 6 coefficients obtained with different methods using the matrix algebra approach.

The results and their statistical analysis are collected in Table 1 for a set of small molecules taken from the database compiled by Tkatchenko and Scheffler [16], as used in [58].

The experimental dispersion coefficients have been determined from dipole oscillator strength distributions (DOSD) [59][60][61][62][63][64][65][66][67][68][69][70][71]. The percentage errors of some of the methods (LDA[M], PBE[M], RHF[M] and RSHLDA[M]) are shown in Figure 3.

Besides Eq. 32, the C 6 coefficients were also calculated from an iterative procedure, where the amplitude matrices T i j were updated according to the first two lines of Eq. 23. Such a procedure corresponds roughly to a local MP2 iteration, and the methods are labelled LDA2, PBE2, RHF2 and RSHLDA2. These results are summarized in Table 2, where a detailed statistical analysis is presented for all the computational results.

The dispersion coefficients obtained from LDA and PBE orbitals are strongly overestimated. It is not really surprising in view of the rather diffuse nature of DFA orbitals and their tendency to underestimate the occupied/virtual gap. Due to the fact that LDA and PBE lead to local potentials, the operator and matrix algebra results are very similar: their difference is smaller than the supposed experimental uncertainty of the DOSD dispersion coefficients (few percent). It is quite clear from Figure 3 that the errors for the molecules containing second-row elements are considerably higher than for most of the molecules composed of H, C, N, O and F atoms. Notable exceptions are the N 2 and N 2 O molecules.

The performance of the method is significantly better for orbitals obtained by nonlocal exchange: RHF and RSHLDA. In the latter case the long-range exchange is nonlocal, and only the short-range exchange is described by a short-range functional. The best performance has been achieved by the RSHLDA[M] method. In contrast to the pure DFA calculations, the matrix and operator algebra methods differ for the RHF and RSH methods significantly: the mean absolute error and the standard deviation of the error is increased by a factor of 2 or 3. This deterioration of the quality of the results reflects the fact that in the operator algebra approach the commutator of the position and non-local exchange operators is neglected (see Appendix F), while this effect is automatically taken into account in the matrix algebra calculations. In view of the simplicity of the model and of its fully ab initio character, the best ME%E of about 11% seems to be very promising and indicates that further work in this direction is worthwhile.

Conclusions, perspectives

It has been shown that, using projected dipolar oscillator orbitals to represent the virtual space in a localized orbital context, the equations involved in long-range ring coupled cluster doubles type RPA calculations can be formulated without explicit knowledge of the virtual orbital set. The POO virtuals have been constructed directly from the localized occupied orbitals. The fockian matrix elements and the electron repulsion integrals were evaluated using only the elements of the occupied block of the Fock/Kohn-Sham matrix in LMO basis and from simple multipole integrals between occupied LMOs. An interesting feature of the model is that, as far as one uses Boys localized orbitals and firstorder dipolar oscillator orbitals, practically all the emerging quantities can be expressed by the overlap integral between two oscillator orbitals, which happens to be the matrix element of the dipole moment fluctuation operator. Various levels of approximations have been considered for the long-range RPA energy leading finally to a pairwise additive dispersion energy formula with a non-emprirical expression for the bond-bond C 6 dispersion coefficient. At this simplest level a straightforward relationship has been unraveled between the ring coupled cluster and the dielectric matrix formulations of the long-range RPA correlation energy.

Our derivation, starting from the quantum chemical RPA correlation energy and pursuing a hierarchy of simplifying assumptions, has led us to a dispersion energy expression which is in a straightforward analogy with classical van der Waals energy formulae. Our procedure produces an explicit model, derived from the wave function of the system, for the dynamical polarizabilities associated with the building blocks, which are bonds, lone pairs and in general localized electron pairs. Thus one arrives to the quantum chemical analogs of the "quantum harmonic oscillators" (QHO) appearing in the semiclassical theory of dispersion forces, elaborated within a RPA framework by Tkatchenko, Ambrosetti, di Stasio and their coworkers [72,73].

A promising perspective of the projected oscillator orbital approach concerns the description of dispersion forces in plane wave calculations for solids. Much effort has been spent recently in finding a compact representation of the dielectric matrix in plane wave calculations, [10] but it still remains a bottleneck for a really fast non-empirical calculation of the long-range correlation energy. Projected oscillator orbitals may offer an opportunity to construct extremely compact representations of a part of the conduction band, which contributes the most to the local dipolar excitations at the origin of van der Waals forces. However, this approach is probably limited to relatively large gap semiconductors, where the construction of the maximally localized Wannier functions (MLWF), [74,75] which are the solid-state analogs of the Boys localized orbitals, is a convergent procedure. Such a methodology could become a fully non-empirical variant of Silvestrelli's extension of the Tkatchenko-Scheffler approach for MLWFs [46][47][48][49][50].

The main purpose of the present work has been to describe the principal features of the formalism, without extensive numerical applications. However, the first rudimentary numerical tests on the molecular C 6 coefficients have indicated that the results are quite plausible in spite of the simplifying approximations and it is reasonable to expect that more sophisticated variants of the method will improve the quality of the model. In view of the modest computational costs and the fully ab initio character of the projected oscillator orbital approach applied at various approximation levels of the RPA, which is able to describe dispersion forces even beyond a pairwise additive scheme, the full numerical implementation of the presently outlined methodology have certainly a great potential for the treatment of London dispersion forces in the context of density functional theory.

A Dipolar oscillator orbitals in local frame

Let R R R i be the rotation matrix which transforms an arbitrary vector v from the laboratory frame to the vector v loc in the local frame defined as the principal axes of the second moment tensor of the charge distribution associated with a given localized occupied orbital

i, R R R i • v = v loc .
The expression in the local frame of a POO i α constructed from the LMO i is then:

|i loc α =        Î - occ. m |m m|        R R R i • (r -D i ) α |i = R R R i • r α |i - occ. m R R R i • m| r |i α |m = β R i αβ r β |i - occ. m β R i αβ m| r β |i |m (39)

B Riccati equations in POO basis

The first order wave function Ψ (1) That is to say that:

Ψ (1) = occ i j virt ab T i j ab | . . . ab . . . | ≈ occ i j POO mαnβ T i j mαnβ | . . . m α n β . . . |. (40)
The canonical virtual orbitals and the POOs in question are related by (see Eq.

and leads to the transformation rule between the amplitudes in the VMO and in the POO basis:

T i j = V T i j POO V † . ( 43
)
Multiplication of the Riccati equations of Eq. 9 by V † and V from the left and from the right respectively, and expressing the amplitudes in POOs using Eq. 43 leads to:

V † R i j V = V † B i j V + V † ( + A) im V T m j POO V † V + V † V T im POO V † ( + A) m j V + V † V T im POO V † B mn V T n j POO V † V = 0, (44)
where implicit summation conventions are supposed on m and n. Recognizing the expression for the overlap matrix, S POO = V † V, we obtain:

R i j POO = B i j POO + ( POO + A POO ) im T m j POO S POO + S POO T im POO ( POO + A POO ) m j + S POO T im POO B mn POO T n j POO S POO = 0, (45)
which defines R i j POO , i j POO , A i j POO and B i j POO and which are the Riccati equations seen in Eq. 12.

C Solution of the Riccati equations in POO basis

To derive the iterative resolution of the Riccati equations seen in Eq. 12, we write explicitly the fock matrix contributions hidden in the matrix . The matrix elements in canonical virtual orbitals, i j ab , read:

i j ab = δ i j f ab -δ ab f i j , (46)
so that the matrix element in POOs is (we omit the "POO" indices):

i j mαnβ = V † mαa i j ab V bnβ = δ i j f mαnβ -S mαnβ f i j . ( 47
)
The terms in the Riccati equations containing the matrix then read (we use implicit summations over m and n):

im T m j S = f T i j S -f im S T m j S (48)
S T im m j = S T i j f -S T im S f m j .

We this in mind, the Riccati equations of Eq. 12 yield:

R i j = B i j + f -f ii S T i j S + S T i j f -S f j j - m i f im S T m j S - m j S T im S f m j + A im T m j S + S T im A m j + S T im B mn T n j S = 0. (50)
Remember that the matrices are of dimension N POO ×N POO . Due to the nonorthogonality of the POOs and the non diagonal structure of the fock matrix, the usual simple updating scheme for the solution of the Riccati equations should be modified in a similar fashion as in the local coupled cluster theory [START_REF] Knowles | Modern Methods and Algorithms of Quantum Chemistry[END_REF]. The fock matrix in the basis of the POOs will be diagonalized by the matrix X obtained from the solution of the generalized eigenvalue problem:

f X = S X ε ε ε. (51)
Note that the transformation X † f X does not brings us back to the canonical virtual orbitals. We can write the transformation by the orthogonal matrix X as X † a iα f iα jβ X jβ b = δ ab ε b , where a and b are pseudocanonical virtual orbitals that diagonalize the fock matrix expressed in POOs. The Riccati equations of Eq. 50 are transformed separately for each pair [i j] in the basis of the pseudo-canonical virtual orbitals that diagonalize f POO :

X † R i j X = X † B i j X + X † f -f ii X † S T i j S X + X † S T i j f X -S X f j j - m i f im X † S T m j S X - m j X † S T im S X f m j + X † A im T m j S X + X † S T im A m j X + X † S T im B mn T n j S X = 0, (52)
which can be simplified by the application of the generalized eigenvalue equation Eq. 51 and the use of the relationships I = S X X † = X X † S:

R i j = B i j + ε ε ε -f ii I T i j + T i j ε ε ε -f j j I - m i f im T m j - m j T im f m j + A im T m j + T im A m j + T im B mn T n j = 0, (53)
with the notations:

R i j = X † R i j X A i j = X † A i j X B i j = X † B i j X T i j = X † S T i j S X.
The new Riccati equations of Eq. 53 can be solved by the iteration formula:

T i j (n) ab = - B i j ab + ∆R i j ab (T (n-1) ) ε a -f ii + ε b -f j j , (54)
where

∆R i j (T) is ∆R i j (T) = - m i f im T m j - m j T im f m j + A im T m j + T im A m j + T im B mn T n j . (55)
As presented here, the update of the "non-diagonal" part of the residue is done in the pseudo-canonical basis. After convergence, we could transform the amplitudes back to the original POO basis according to:

T i j = X † S -1 T i j S X -1 = S -1 (X † ) -1 T i j X -1 S -1 = X X † (X † ) -1 T i j X -1 X X † = X T i j X † . (56)
However, this back-transformation is not necessary since the correlation energy can be obtained directly in the pseudo-canonical basis, as:

occ i j tr B i j T i j = occ i j tr X † B i j X X † S T i j S X = occ i j tr B i j T i j S X X † = occ i j tr B i j T i j . (57)

D Riccati equations in the local excitation approximation

The local excitation approximation imposes that in the matrices R, , A and B, the excitations remain on the same localized orbitals. In this approximation the Riccati equations of Eq. 12, with explicit virtual indexes, read: R i j iα jβ = B i j iα jβ + ( + A) im iαmγ T m j mγ pδ S pδ jβ + S iα pγ T im pγmδ ( + A) m j mδ jβ + S iα pγ T im pγmδ B mn mδnτ T n j nτqζ S qζ jβ = 0.

In this context, the terms containing the matrix are (with explicit POO indexes):

im iαmγ T m j mγ pδ S pδ jβ = f iαiγ T i j iγ pδ S pδ jβf im S iαmγ T m j mγ pδ S pδ jβ

S iα pγ T im pγmδ m j mδ jβ = S iα pγ T i j pγ jδ f jδ jβ -S iα pγ T im pγmδ S mδ jβ f m j .

Inserting this in Eq. 58 and using the shorthand notation R i j iα jβ ≡ R i j αβ , B i j iα jβ ≡ B i j αβ , f iα jβ ≡ f i j αβ and S iα jβ ≡ S i j αβ , (note the matrix elements T i j iγ pδ cannot yet be translated to the shorthand notation) one obtains:

R i j αβ = B i j αβ + f ii αγ T i j iγ pδ S p j δβ -f im S im αγ T m j mγ pδ S p j δβ + A im αγ T m j mγ pδ S p j δβ + S ip αγ T i j pγ jδ f j j δβ -S ip αγ T im pγmδ S m j δβ f m j + S ip αγ T im pγmδ A m j δβ + S ip αγ T im pγmδ B mn δτ T n j nτqζ S q j ζβ = 0. ( 61
)
It is then a further approximation to tell that the POOs coming from different LMOs have a negligible overlap, i.e. that S i j αβ = δ i j S ii αβ . The Riccati equations become: R i j αβ = B i j αβ + f ii αγ T i j iγ jδ S j j δβf ii S ii αγ T i j iγ jδ S j j δβ + A im αγ T m j mγ jδ S j j δβ + S ii αγ T i j iγ jδ f j j δβ -S ii αγ T i j iγ jδ S j j δβ f j j + S ii αγ T im iγmδ A m j δβ + S ii αγ T im iγmδ B mn δτ T n j nτ jζ S j j ζβ = 0,

which, in turn, allows us to use the shorthand notation T iα jβ ≡ T i j αβ to arrive at: R i j = B i j + f ii T i j S j jf ii S ii T i j S j j + A im T m j S j j + S ii T i j f j j -S ii T i j S j j f j j + S ii T im A m j + S ii T im B mn T n j S j j = 0,

E Screened dipole interaction tensor

Any interaction L(r) can be expanded in multipole series using a double Taylor expansion around appropriately selected centers, here D i and D j , such that, with r = r ir j = (r i -D i ) + D i j -(r j -D j ) where

D i j = D i -D j : L(r) = L i j (D i j ) + α (r i α -D i α ) L i j α (D i j ) + α (r j α -D j α ) L i j α (D i j ) + αβ (r i α -D i α )(r j β -D j β )L i j αβ (D i j ) + . . . , (64)
where the definitions of L i j α (D i j ), L i j αβ (D i j ) are obvious. For example, in the case of the long-range interaction, L(r) will be defined according to the RSH theory as:

L(r) = erf(µr) r , (65)
with r = |r|. The multipolar expansion of the long-range interaction leads to the following first and second order interaction tensors:

L i j α (D i j ) = - D i j α D i j 3 1 - 2 √ π D i j µ e -µ 2 D i j 2 -erf(µD i j ) ( 66
) L i j αβ (D i j ) = 3 D i j α D i j β D i j 5 erf(µD i j ) - 2 3 √ π D i j µ e -µ 2 D i j 2 3 + 2D i j 2 µ 2 - δ αβ D i j 2 D i j 5 erf(µD i j ) - 2 √ π D i j µ e -µ 2 D i j 2 . ( 67
)
Remembering that the full-range Coulomb interaction tensor reads T i j αβ (D i j ) = 3 D i j α D i j β -δ αβ D i j 2 D i j -5 , the long-range interaction tensor can be written in an alternate form which clearly shows the damped dipole-dipole interaction contribution: L i j αβ (D i j ) = T i j αβ (D i j ) erf(µ D i j ) -

2 3 √ π D i j µ e -µ 2 D i j 2 3 + 2D i j 2 µ 2 -δ αβ e -µ 2 D i j 2 4 µ 3 3 √ π . ( 68
)
The trace of the tensor product (used for the spherically averaged C 6 ) then reads:

αβ L i j αβ L i j αβ = 6
D i j 6 4e -2D i j 2 µ 2 D i j µ D i j µ 3 + 4D i j 2 µ 2 + 2D i j 4 µ 4 3π -3 + 2D i j 2 µ 2 erf(D i j µ)

3 √ π + erf(D i j µ) 2 = 6 D i j 6 F µ damp (D i j ). ( 69
)

F Fock matrix element in POO basis

The occupied-occupied block of the fock matrix, f i j , is known. The POOs are orthogonal to the occupied subspace of the original basis set, they satisfy the local Brillouin theorem, i.e. the occupied-virtual block is zero. As a result, in the local excitation approximation, we need only to deal with the fock matrix elements f ii αβ :

f ii αβ = i α | f |i β = i| rα Q f Q rβ |i , (70)
from which we directly derive the quantity f i [M] of Eq. 36:

f i [M] = α f ii αα = virt ab i| rα |a f ab b| rα |i . (71)
Since we would like to express everything in occupied orbitals, we expand the projector Q and use that the occupied-virtual block of the fock matrix is zero to obtain the following expression:

f ii αβ = i| rα f rβ |i -
Note that this holds provided that the fockian contains only local potential terms, which commute with the coordinate operator: see later for the more general case. In this special case, the double commutator can be written as

The diagonal matrix element of the triple operator product is then:

i| rα f rβ |i = 1 2 δ αβ + 1 2 i| rα rβ f |i + i| f rβ rα |i . (76)
Since the localized orbitals satisfy local Brillouin theorem, we finally obtain for the matrix elements of the fock operator with multiplicative potential (typically Kohn-Sham operator with local or semi-local functionals) between two oscillator orbitals:

f ii αβ = 1 2 δ αβ + 1
From this, we obtain the quantity f i [O] of Eq. 38:

f i [O] = α f ii αα = 3 2 + 1
In the more general case, i.e. when the fockian contains a nonlocal exchange operator, like in hybrid DFT and in Hartree-Fock calculations, the relation seen Eq. 73 does not hold any more and the commutator of the position operator with the fockian contains an exchange contribution [76,77], which gives rise to an additional term:

) 80
where Prr is the permutation operator that changes the coordinates r appearing after K to r, and we recall that w(r, r ) is the two-electron interaction. Hence, the diagonal blocks of the POO fockian in the general case can be written as:

i α | f |i β = 1 2 δ αβ
In the present work, the exchange contribution, which is present only in the case of Hartree-Fock or hybrid density functional fockians and which would give rise to non-conventional two-electron integrals, is not treated explicitly. Possible approximate solutions for this problem will be discussed in forthcoming works. Although we do not use the elements of the off-diagonal (i j) blocks of the POO fock operator, for the sake of completeness we give its expression: (82)

To derive this, instead of the double commutator of Eq. 73, one needs to consider the product of the commutator with a coordinate operator

rα [r β , f ] = rα [r β , T ] -rα [r β , K], ( 83
)
where T is the kinetic energy operator, and [r β , T ] = ∇β .

Fig. 2

 2 Fig. 1 Projected Oscillator Orbitals ((b),(c),(d)) generated by projection of the products of first order solid spherical harmonics polynomials with the O lone pair orbital of H 2 C--O seen in (a). The harmonics are aligned with the local frame axes, i.e. the principal axes of the tensor of the moment of inertia of the charge distribution of the oxygen lone pair orbital. The green dot indicates the position of the LMO centroid.

  6): |m α = virt. a |a V amα and: |n β = virt. b |b V bnβ . (41) This allows us to write: mαnβ | . . . m α n β . . . nβ b | . . . ab . . . |,

  occ mn i| rα |m f mn n| rβ |i .(72)In order to transform the triple operator product, rα f rβ , let us consider the following double commutator:rα , [r β , f ] = -δ αβ .

  rα , [r β , f ] = rα rβ frα f rβrβ f rα + f rβ rα = -δ αβ ,(74)which allows us to express the two triple products: rα f rβ + rβ f rα = δ αβ + rα rβ f + f rβ rα .

  |m f mi + f im m| rα rβ |i -occ mn i| rα |m f mn n| rβ |i .

  r2 |m f mi + f im m| r2 |i -occ mn α i| rα |m f mn n| rα |i .

  i| rα , [r β , K] |i = occ m im| rαr α w(r, r ) rβr β |mi ,(79)where the nonlocal exchange operator is defined asK = occ m dr φ † m (r ) w(r, r ) Prr φ m (r ), (

  rαr α w(r, r ) rβr β |mi |m f mi + f im m| rα rβ |i -occ mn i| rα |m f mn n| rβ |i .

  i α | f | j β = -i| rα ∇β | j +occ m im| rα w(r, r ) rβr β |m j + occ m i| rα rβ |m f mi occ mn i| rα |m f mn n| rβ |i .

Table 1

 1 Molecular C 6 coefficients from the dipolar oscillator orbital method, using LDA, PBE, RHF and sr-LDA/lr-RHF determinants in augcc-pVTZ basis set and Boys localized orbitals. The matrix elements are calculated with a matrix algebra [M] and operator algebra [O] approach, respectively. All results are in atomic units. Values highlighted as bold indicate the best agreement with experiment..

	Molecule	Ref.	LDA[M]	LDA[O]	PBE[M]	PBE[O]	RHF[M]	RHF[O] RSHLDA[M] RSHLDA[O]
	H 2	12.1	18.4	18.4	17.2	17.2	9.6	16.2	10.6	16.9
	HF	19.0	20.7	21.5	20.7	21.6	12.4	17.2	14.9	19.3
	H 2 O	45.3	55.6	56.5	55.1	56.1	31.9	45.4	36.4	48.9
	N 2	73.3	118.3	118.8	117.6	118.0	75.2	109.0	81.1	110.9
	CO	81.4	120.4	120.9	119.1	119.7	70.9	104.1	79.1	109.6
	NH 3	89.0	118.8	119.9	116.7	117.8	66.3	98.6	72.6	102.0
	CH 4	129.7	192.4	193.0	185.9	186.5	105.3	163.5	114.5	168.2
	HCl	130.4	208.9	211.8	205.1	208.5	121.5	186.3	126.7	184.8
	CO 2	158.7	234.9	237.0	233.2	235.4	130.9	184.0	153.3	202.8
	H 2 CO	165.2	205.7	207.2	202.8	204.3	115.2	168.3	129.4	178.4
	N 2 O	184.9	317.1	319.3	315.3	317.5	178.9	252.3	206.0	274.7
	C 2 H 2	204.1	343.5	345.2	340.4	342.2	206.4	316.1	209.7	306.5
	HBr	216.6	303.3	325.9	302.5	325.4	188.3	293.5	188.2	282.7
	H 2 S	216.8	392.8	397.7	382.8	388.2	213.0	339.2	220.5	335.5
	CH 3 OH	222.0	303.8	305.3	297.0	298.5	166.8	246.2	187.2	260.2
	SO 2	294.0	542.6	554.9	539.8	552.5	284.9	416.7	329.5	456.2
	C 2 H 4	300.2	466.1	467.8	456.7	458.4	266.3	406.4	279.4	405.8
	CH 3 NH 2	303.8	440.6	442.3	429.8	431.4	240.2	360.3	263.2	373.3
	SiH 4	343.9	639.6	655.3	598.4	613.9	280.0	484.6	310.9	513.0
	C 2 H 6	381.9	579.2	580.9	560.9	562.5	313.5	480.9	341.2	495.5
	Cl 2	389.2	727.4	735.4	714.3	723.7	401.3	607.5	424.3	612.2
	CH 3 CHO	401.7	627.5	630.4	613.3	616.0	333.2	493.6	372.6	521.7
	COS	402.2	845.7	855.8	843.9	853.7	456.3	689.2	495.7	713.5
	CH 3 OCH 3	534.1	781.4	784.3	758.4	761.1	415.2	619.4	464.8	654.3
	C 3 H 6	662.1	1045.9	1049.4	1018.6	1021.8	571.2	868.2	609.8	881.2
	CS 2	871.1	2099.4	2119.5	2073.4	2094.5	1085.6	1658.0	1144.0	1686.2
	CCl 4	2024.1	3831.3	3861.0	3750.4	3784.0	2004.8	3007.0	2135.5	3051.9
	MAD%E		59.84	61.69	56.71	58.62	15.22	33.75	11.84	37.60
	STD%E		28.08	28.25	27.76	27.98	9.88	21.15	7.24	20.86
	CSSD%E		67.14	68.92	64.11	65.97	18.39	40.37	14.07	43.62
	150
	100
	50
	0
	-50

Table 2

 2 Detailed statistical analysis of the methods to obtain POO C 6 coefficients.

		LDA[M]	PBE[M]	RHF[M]	RSHLDA[M]
	MA%E	59.8	56.7	15.2	11.8
	STD%E	28.1	27.8	9.9	7.2
	CSSD%E	67.1	64.1	18.4	14.1
	MED%E	55.3	52.1	17.1	11.4
	MAX%E	141	138	24.6	31.3
	MIN%E	9.1	8.8	-34.7	-21.7
		LDA2[M] PBE2[M] RHF2[M] RSHLDA2[M]
	MA%E	64.3	61	17.1	17.1
	STD%E	35.7	35.3	12.3	12.3
	CSSD%E	74.7	71.5	21.3	21.3
	MED%E	53.5	50.6	15.3	15.3
	MAX%E	142.1	148.7	54.3	54.3
	MIN%E	-4.7	-5.1	-36.6	-36.6
		LDA[O]	PBE[O]	RHF[O]	RSHLDA[O]
	MA%E	61.7	58.6	33.7	37.6
	STD%E	28.3	28	21.2	20.9
	CSSD%E	68.9	66	40.4	43.6
	MED%E	55.8	52.7	34	34.6
	MAX%E	143.3	140.4	90.3	93.6
	MIN%E	13.4	13.4	-9.7	1.7
		LDA2[O]	PBE2[O]	RHF2[O]	RSHLDA2[O]
	MA%E	65.9	62.7	51	51
	STD%E	36.3	36.9	47.7	47.7
	CSSD%E	76.3	73.8	70.5	70.5
	MED%E	54.7	51.1	34	34
	MAX%E	143.3	152.9	214.8	214.8
	MIN%E	-1.1	-1.2	-10.6	-10.6

  can be written in terms of Slater determinants | . . . ab . . . | formed with LMOs and canonical virtual orbitals a and b on the one hand, and on the other hand in terms of Slater determinants | . . . m α n β . . . | formed with LMOs and POOs m α and n β .

Acknowledgements J.G.A. thanks Prof. Péter Surján (Budapest), to whom this article is dedicated, the fruitful discussions during an early stage of this work.