Anthony Scemama

Antonio Monari

Celestino Angeli

Stefano Borini

Stefano Evangelisti

Elda Rossi

Common Format for Quantum Chemistry Interoperability: Q5Cost format and library

A particular data format for interchange, specifically designed for Quantum Chemistry programs, was proposed within a COST activity. The new features of the recent version 1.0 are presented, together with the Q5Cost library that allows to use the data format within a user program. The problem of a general description of the wave function is presented as well as the adopted strategy for including it into the data format. Several utilities have been delivered for the use of the developers and the users of the library, as well as various converters for well known Quantum Chemistry codes.

Introduction

The present article describes the design, setup and upgrading of a data format and a FORTRAN library, Q5Cost, for the management of data produced by a generic abinitio Quantum Chemistry (QC) code. The data structure has been defined with the aim of making code interoperability easier in the scientific area of Quantum Chemistry. This activity has been carried out within the COST network, supporting cooperation among researchers in Europe: the former D23 "Metachem" action, now concluded, and the D37 "GridChem" action, started in July 2006 [1].

The first problem to be faced when integrating different QC codes in a common workflow is the lack of a standard and the different data formats adopted by each code. Our suggestion was to adopt such a standard, specifically designed for interchange, and to connect each code in the set through a converter. Of course not to invent "yet another format", we strongly tried to design a format as general and flexible as possible and to coordinate ourselves with other similar initiatives in the quantum chemistry context.

The first idea was to adopt two different technologies for describing large binary data and small coded information. HDF5 (Hierarchical Data Format) [START_REF]HDF5 a general purpose library and file format for storing scientific data[END_REF] was the technology of choice for the first type of data, and XML (eXtensible Markup Language) [START_REF]HDF5 a general purpose library and file format for storing scientific data[END_REF] for the second one.

Nevertheless, during the project development, important considerations lead us to the decision to collect all the data structures into a unique data format and, in particular, to design a single user interface for managing all the data. A first definition of the data model and a description of the library has been reported elsewhere [START_REF] Borini | [END_REF]4]. In this paper we want to report the state of the art of the library, now at the version 1.0, in addition to a critical examination of the improvements and modifications carried on it.

The next section will discuss the improvements carried on in the data format, in terms of "Geometry", "Basis set" and "Wave function". The same issues are then considered in the following section from the user interface point of view. The strategy adopted for storing the wave function is then discussed in detail, since this is a point that poses difficult problems due to the size and wide variety of possible wave functions. A number of utilities have been set up to facilitate the use of the Q5Cost library and are described in the next session. The utilities are addressed to the library developers but also to the programmers that are using the library in their programs to manage the files produced by a QC calculation. Wrappers are the last topic presented in this paper; they are specific translation programs that use the Q5Cost library for converting the proprietary data format into the Q5Cost data format and vice versa.

The data format

The structure of the first versions of the Q5Cost data format has been fully described elsewhere [START_REF] Borini | [END_REF]4], as well as its connections with the inherent HDF5 structure [START_REF]HDF5 a general purpose library and file format for storing scientific data[END_REF]. Here we just want to briefly recall the main features and present in details the modifications introduced with the present version.

Q5Cost is a common interchange format for data coming from ab initio quantum chemistry calculations. In this context ab initio refers to methods that are based on a wave function rather than an electronic density. In particular, the Density Functional Theory (DFT) methods at the moment are not supported by our format. This choice is justified by the fact that the majority of the codes involved in this project are based on Slater-Determinant expansions of the wave function.

The q5Cost format is flexible and extensible, so we can design it in an incremental way. The data format consists of a collection of chemical objects related within a hierarchical structure in a logical containment relationship as reported in Figure 1. Originally the Q5Cost data format was designed to handle only the large binary data coming from Quantum Chemistry, together with the information (metadata) needed for describing them. The data structures described by the format were basically integrals, on the Atomic and Molecular orbitals, the accompanying metadata entities like electrons number, symmetry indices, nuclear energy, molecular orbital labels, but also coefficients needed to define symmetry adapted orbitals or to transform the Atomic orbitals into Molecular orbitals.

Other small pieces of information, mainly ASCII coded in conventional QC data formats, are needed to define a QC system. They can be classified into three classes: Geometry, Basis set, spatial Symmetry of the system under investigation. In the original version of the data format those data were represented using an XML based language in a separated file [START_REF] Borini | [END_REF]. Starting from this version, we decided to merge them into the Q5Cost data format, in order to ensure a unique coherent structure grouping all the specifications on the chemical system. This is a new feature of the format that now includes the whole information about a molecular system.

Data represented into the Q5Cost format belong to two different categories:

1. Data: the large binary quantities used in quantum chemistry for representing integrals, properties and wave functions. They can be stored using matrices with an arbitrary number of indices (rank-n arrays), scale aggressively with the system size, and are normally accessed with a "chunked" approach (i.e., using welldefined blocks of data). 2. Metadata: simple and small pieces of data that describes and better defines the previous data. They represent well-known chemical entities like nuclear energy, molecular orbital labels, and molecular symmetry and can be stored as scalars, vectors or matrices. For example, the nuclear repulsion energy is a floating point scalar, molecular orbitals coefficients are an (N,M) floating point matrix, the associated orbital energies are a floating point vector, the molecular orbital labels are a vector of strings, and so on. Among the "data" entities we can include two-electron integrals, or atomic orbital overlap, but also other more application-specific objects, like the four particle density matrix. In the modern approaches based on localized orbitals, these matrices are of sparse nature; this encourages the storage of only nonzero elements, each one associated to n indices in the case of a rank-n array. This representation of the data, although not particularly efficient in terms of space occupation, is well-known by the interested parties, easy to debug, and already integrated in the current code-base both for memory representation and file storage of data.

A great part of these data objects can be described by a "generic property" object, provided that we define the matrix rank, the involved operator(s) and basis functions. Nevertheless, since some of these "properties" are well-known chemical entities, and chemists are used to referring to them by name, we provide a specific library access to most of them (overlap, one-electron integrals, two-electron integrals...), in addition to an interface to the "generic property" for handling other properties not explicitly provided by the library. This should ensure both ease of use and general adaptability of the library to either alternative or future theoretical developments.

All the chemical objects described in the Q5Cost data format are related within a hierarchical structure, and logical containment relations can be defined for them. A first (root) container, named System, represents the molecular system as defined by its structural data (chemical composition and spatial geometry). All the general metadata, that cannot be integrated into a specific domain, can be associated to this container. A system can contain several "Domains". The role of the Domain is to group together entities whose indices conceptually refer to the same kind of functions. In the present version three Domains have been defined: 1. Atomic Orbital (AO): refers to the data defined on the AO basis, overlap, oneelectron integrals, two-electron integrals and the generic property, i.e. any other property that can be described on the AO basis (dipole moment integrals, for example). This domain contains also the definition of the Basis Set. 2. Molecular Orbital (MO): refers to the data defined on the MO basis, one-and twoelectron integrals and the generic property. This domain contains also the transformation matrix needed to define the MO on the AO basis. 3. Wave Function (WF): refers to the definition of the wave function. This will be described if full details in the next section.

4.

One-and two-electron integrals and the generic property are present in both AO and MO domain but they are different quantities. The first set is defined on the AO basis and the second is obtained from the first one thought the transformation matrix, i.e. the matrix containing the molecular orbitals coefficients.

In the data format are present, with respect to the previous version, two additional groups of information about Geometry and Basis set.

Geometry is specified by giving the Cartesian coordinates, the atomic number and label of each atom unique by symmetry. Only a single geometry can be stored, at present, in a given Q5Cost file and this group is contained in the System container. We are considering the possibility to store multiple geometries in a single file, in order to facilitate the treatment of problems such as the study of potential energy surfaces or direct dynamics.

The basis set information are specific of the AO domain, so in a given Q5Cost file several AO definitions are possible. The AO domain (like all the other domains) can be splitted into several sub-domains, each for every different AO choice. The different AO integrals can arise from different basis set or from different orbital types (Atomic Orbitals vs. Symmetry Adapted Linear Combination of Atomic Orbitals)

. Concerning the wave function, complex issues arise because no clear agreement is reached on the many different wave functions that could be stored. A general agreement has been found that the lowest level representation of a wave function is a multi-reference configuration interaction, where each determinant is described together with a proper coefficient. This allows a complete description of the wave function, and should therefore be generally present among the stored information. We do however realize the following points: 1. A multireference description lacks any additional information about the true nature of the wave function (examples like HF, CAS, RAS, CAS+S etc..). Description of higher level information is important to prevent loss. Q5Cost must provide access to proper storing of diversified information to fully qualify the described entity.

A design assumption about the concept of Domain is that Properties living into a

given Domain have indexes referring to a particular concept belonging to that Domain. As an example, the molecular orbitals one-electron integrals refer to molecular orbitals, and they belong to the MO Domain..

The Library

The new interface for the user provides access for handling geometry and basis set information. Furthermore, storage of the wave function is discussed.

The "Geom" routines allow the client to store and retrieve the molecular geometry. The geometry is considered top level information of the chemical system under investigation, therefore is associated with the System. At present only one System, and consequently only one geometry, per Q5Cost file is allowed. The restriction to a single System per file will be possibly lifted in future releases of the library. This will require, however, a slightly modified programming interface, and considerations about inter-System characteristics, currently outside the scope of our efforts. The actual data handled by these routines are the Cartesian coordinates, the atomic numbers and the labels for all the atoms, unique by symmetry, contained in the system. The only metadata needed by these routines is the number of atoms. If in future releases other coordinate types are going to be supported in addition to the Cartesian, specific metainformation will be added.

Information about the Basis Set is associated to the AO domain, since multiple basis sets can exist in a given Q5Cost file. Its presence is bound to the presence of a Geometry description. The routines allow the user to store and retrieve the exponents and contraction coefficients of each contracted function defined for each angular momentum on each atom of the system. The actual data handled by the basis set routines are the exponents and the coefficients defining the contraction. The required metadata include the type of angular functions, angular momentum, and the number of functions for each contraction. A possible link to the EMSL data bank [START_REF]Extensible Computational Chemistry Environment Basis Set Database[END_REF] is advisable in the future. Basis set names like "6-31G*" could also be an option, however, these names do not represent a standard nor in terms of naming, neither in terms of actual contents among different quantum chemistry programs. As a consequence, the plain use of basis sets' simple name should be discouraged, because it makes difficult to check for the consistency between two stored bases, and hence preventing easy interoperability.

The routines for storing and retrieving the wave function are still a matter of debate in the working group, as better described in the following section. A tentative set of routines are available anyway for a Configuration Interaction type Wave Function (CI-WF), in terms of a linear combination of Slater determinants made of Molecular orbitals. The data to be considered are the coefficients (real numbers) and the Slater determinants that are part of the combination. The metadata needed depends on the type of WF to be described and includes the energy and the core energy associated to the WF, the Molecular orbitals to be used for the determinants and their classification.

The Wave function

The storage of the Wave Function (WF) is a particularly relevant point, and poses different problems because of its size and the wide variety of possible WF's. The wave function is usually a very large object, and its storage often represents the largest amounts of data needs in a quantum-chemistry calculation. Several options are possible, each one typically suited for a particular type of WF calculation.

Usually a WF is called of Configuration Interaction (CI) type if it is expressed as a linear combination of Slater determinants (thereafter "determinants"). Therefore, a CI WF is expressed in terms of a list of Slater determinants and the corresponding coefficients of the linear expansion. We notice that this is, by no means, the only possibility. In fact, many other different WFs are currently used in Quantum Chemistry, depending on the level of theory used to study the system. For instance: Coupled Cluster (CC) Contracted CI (C-CI), Density-Matrix Renormalization Group (DMRG), and many others. At the moment, we decided to restrict our work on CI WF because the interest of all the partners of the project is on this type of WF. We plan in the future to extend our definition to more general WF's. In particular, CC or C-CI wave functions could be easily stored by simple generalization of the procedure used for a CI WF; some works is in progress on this point. Things are more complex for DMRG, and we do not plan, for the moment, to deal with this kind of formalism.

Even with this restriction, the different ways a WF can be defined are rather numerous. 1. SCF (or HF): If a single determinant is used to describe the WF, this is called of Self Consistent Field (SCF), or, equivalently, Hartree-Fock (HF) type. One should notice that the Molecular Orbitals are usually variationally optimized in order to define an SCF WF. 2. FCI: If all the Slater determinants are considered, the WF is called of Full Configuration Interaction (FCI) type. 3. CAS: In the Complete Active Space Self-Consistent Field (CAS-SCF) formalism, the orbitals are classified into three disjoint classes whose union gives the total set of orbitals:

-Occupied Orbitals (O), if they are doubly occupied in all CAS determinants;

-Virtual Orbitals (V), if they are empty in all CAS determinants:

-Active Orbitals (A), if their occupation varies in the set of the CAS determinants, i.e., each A orbital can be either singly or doubly occupied or empty in the CAS determinants. 4. It is important to note that in the special case when all the orbitals are in the Active class (O and V classes empty) the CAS-type WF degenerate into a FCI WF. In other words, the FCI type WF is a special case of the more general CAS-type WF when the O and V classes are empty. On the other hand, when the A class is empty and all orbitals are either in the Occupied or Virtual classes, the CAS-type WF become an SCF one. In conclusion we can assume the CAS-type WF to be a more general type of CI wave function with FCI and SCF wave functions as special cases. 5. TCI: in the Truncated CI (TCI) formalism, the determinants are obtained by exciting a given number of electrons from a set of "Reference" determinants. In the most common case, the reference determinants have a CAS structure, and the excitations are restricted to Singles (S) or Singles+Doubles (SD). One obtains then the CAS-S and CAS-SD CI spaces, respectively. 6. Selected-CI: in this case the set of determinants can be completely arbitrary. This is the case, for instance, if the determinants are derived from a TCI or a FCI expansion by selecting those determinants whose coefficients have an absolute value larger than a given threshold. In the first four cases, it is possible to define an ordering of the Slater determinants, and the WF is completely defined by specifying -The WF nature (SCF, FCI, CAS-CI, etc), and the corresponding information that uniquely specify it (level of truncation in TCI, nature of the CAS, etc) -The ordering used for the determinants; -The corresponding CI coefficients.

In the other case, described at point 5., when this is not possible, for example because the determinants involved in the CI expansion are completely arbitrary or their generation schema is too complex, it is necessary to explicitly define all the determinants. In order to describe the CI WF from a general point of view, the simplest way is to give the determinants list without relying on a possible order. The price to pay, of course, is that this piece of information can be very large and space demanding but also redundant in cases like WF's of type 1-4.

Generally speaking, the WF storage requires two data structures: 1. the list of the determinants used to define the WF; 2. the corresponding list of the coefficients of the CI expansion. The second structure is simply a vector of real numbers, while the first one is a vector of more complicated objects, since there are several not trivial ways to code a Slater determinant.

The first and common way to define a Slater determinant on the basis of Molecular Orbitals is to indicate, for each electron in the system, the spin orbital that hosts it, as shown in Figure 2. A similar way, that sometimes can be more effective, is to use vectors running on the spinorbitals and reporting the occupation (1 if the orbital is occupied by an electron, 0 if empty), as shown in Figure 3. When working with a WF expanded on determinants obtained by an excitation process on CAS determinants, a hole/particle formalism is more effective to describe the occupation of the orbitals in the Occupied and Virtual classes. In this formalism it is advisable to define a reference occupation pattern in the Occupied and Virtual classes (the trivial choice is occupation 2 and 0 for the two classes, respectively), and simply define the occupation in these classes of any other determinant. The differences can be "holes" (an orbital that contained an electron in the reference is now empty) or "particles" (an orbital that was empty in the reference is now full).

 

Num. of electrons

Quite often this formalism is used to take care of the orbitals that are in the Occupied or Virtual classes. In the O class, the number of holes n h in spinorbitals is specified, followed by the list of the corresponding holes:

n h , (h(i),i=1, n h )
The same is done for the spinorbitals in the V class where the number of particles np is specified, followed by the list of particles:

n p , (p(i),i=1,n p )

Utilities

A number of utilities have been set up to facilitate the use of the Q5Cost library and data format. The utilities are addressed to the library developers or to the library installers; some of the utilities are for the use of the final users, the programmers that  

Num. of MO

Orbital occupation 0 = empty 1 = occupied need to include calls to the routines in their programs and manage the data files produced by a QC calculation.

The user interface to the data format is in Fortran 90, because of a great flexibility of this language and its large use for chemical programming. Nevertheless, other interfaces are available for those programmers which prefer Fortran 77 or other languages. A Python script reads the library code and builds automatically an interface for calling Q5Cost functions in Fortran 77 codes, which exposes the library functionalities to Fortran77 and C programmers. The same script builds also a C++ interface for the use of the library in a C++ environment which is very common in visualization software. In turn, starting from the C++ interface and using SWIG [START_REF]SWIG is an interface compiler that connects programs written in C and C++ with scripting languages such as Perl, Python, Ruby, and Tcl[END_REF], another interface is built for the Python programming language, a well suited environment for writing wrappers.

In the previous releases, we remarked that the users had some difficulties to link all the necessary libraries (HDF5, HDF5-fortran, libz, libm, Q5Cost) in their codes, because the linking process needs a specific order of the libraries in the command line. To simplify the linking of the codes, we chose to put all the needed libraries into one big static library, containing also the Fortran77 and the C++ interfaces. Now, only one file needs to be added in the linking process.

The user documentation is automatically extracted from the source code by a Python tool that can generate by choice a standard HTML file or a text file. This tool is designed in a modular way in order to make it possible to include new formats when needed.

As far as the installation of the library on personal systems is concerned, an automatic configuration script is now available to prepare the library makefile. It takes into account the specificity of the computing platform, like the compilers, the location of the libraries, and verifies if the needed underlying HDF5 layer is present. This configuration script is based on autoconf , an extensible package of macros that produces shell scripts to automatically configure software source code packages and guarantees the portability on different architectures; autoconf has been developed within the GNU project [START_REF]Gnu autoconf[END_REF] Two useful utilities for the final user are available to inspect the q5 file without writing an ad-hoc program. As HDF5 is a binary format, one can argue that the 'small data' such as the molecular geometry, the basis set or the molecular orbitals can't be read easily by humans. On the contrary, since it is based on HDF5, it contains its own description and it is quite easy, using standard HDF5 instruments, to extract a complete description in a human readable form.

-q5dump is a simple tool that reports the content of a Q5Cost file, listing the data and the metadata contained in it.

-q5edit is a more complete and attractive tool. Using this Python interactive program, the user can navigate through the file, listing the content and modifying the data using his favorite text editor, defined by an environment variable. The Python interface to Q5Cost is used to access the file, and the "ncurses" library [START_REF]Gnu ncurses is a programming library allowing the programmer to write text user interfaces in a terminal-independent manner[END_REF] is used for the interaction of the user with the program.

Wrappers

The Q5Cost library, as a tool for helping data exchange among different chemistry programs, has been used to build specific translation programs [START_REF] Borini | [END_REF] to/from the proprietary data format of each of the involved programs. These translation programs, also called "wrappers", are small programs that, using both the program specific I/O library and the Q5Cost library, take care of converting the proprietary data format into the Q5Cost data format and vice versa.

Several QC codes have been wrapped up to now, both general purpose and more specific "home-made" ones usually written and maintained by the partners of the project. In the first class we recall for instance Dalton [START_REF]DALTON a molecular electronic structure program[END_REF], Molcas [START_REF] Andersson | Molcas version 6.2[END_REF] and GAMESS US [START_REF] Schmidt | General Atomic and Molecular Electronic Structure System[END_REF][START_REF] Gordon | Theory and Applications of Computational Chemistry: the first forty years[END_REF], in the second class the Bologna FCI [START_REF] Bendazzoli | [END_REF]14] code and the Toulouse chain [15,16] among others. In many cases the wrappers have been written in collaboration with research groups outside the authors of this paper. Here we are simply giving a brief summary, while a more exhaustive and detailed description will be the subject of a forthcoming future paper.

Wrappers are, anyway, normally written and maintained by the author of the chemistry programs. They are external programs that take care to translate the proprietary data to the Q5Cost format (output wrapper) and the common format to the proprietary format (input wrapper). In other cases, the authors liked better to include the Q5Cost routines in the code itself, so to allow the program to read and write "natively" the Q5Cost format. This strategy was used, for example, for the FCI code and for the GAMESS US package: it is more simple and direct of the previous one but it requires a greater involvement of the authors of the codes and it could be difficult in the case of large commercial packages.

The programs that can now use the Q5Cost data format are reported in the following list. In some cases the implementation is still experimental, and has not been released yet in the official versions.

-Dalton [START_REF]DALTON a molecular electronic structure program[END_REF] -GAMESS US [START_REF] Schmidt | General Atomic and Molecular Electronic Structure System[END_REF][START_REF] Gordon | Theory and Applications of Computational Chemistry: the first forty years[END_REF] -Molcas [START_REF] Andersson | Molcas version 6.2[END_REF] -Columbus/ACES II [START_REF] Columbus | an ab initio electronic structure program[END_REF][START_REF] Aces | [END_REF] -Molekel [19] -Toulouse Chain [15,16] -Ferrara 4-indices transformation [START_REF] Cimiraglia | Angeli private communication[END_REF] -Paris 4-indices transformation [START_REF] Rehinardt | [END_REF] -Bologna FCI code [START_REF] Bendazzoli | [END_REF]14] -QMC=Chem For the Bologna FCI code [START_REF] Bendazzoli | [END_REF]14], the Q5Cost format is now the native data format. Thanks to the integration of FCI with the other codes it was possible to face new scientific problems: the study of the basis set superposition error on the Ne dimer [START_REF] Monari | [END_REF], the metal insulator transition [23] and the behavior of small alkali metal clusters [24]. The large systems treated involved a FCI space larger than 1 billion determinants and were good benchmark for the Q5Cost data format and library, which proved to be effective, usable and flexible. The performance outcome, both in terms of disk occupation and access time to I/O, was acceptable too, confirming the previous tests [4].

The Toulouse chain [15,16] has been only partially integrated. The complete support is however strongly required because it will open the way to the integration of an effective 4-indices transformation, like the one developed in Paris [START_REF] Rehinardt | [END_REF]. The performance of the integral transformation is of crucial importance when local approaches are involved. GAMESS US [START_REF] Schmidt | General Atomic and Molecular Electronic Structure System[END_REF][START_REF] Gordon | Theory and Applications of Computational Chemistry: the first forty years[END_REF] is being integrated in collaboration with the authors, by including the Q5Cost routines directly in the program I/O interface. The Dalton [START_REF]DALTON a molecular electronic structure program[END_REF] package was originally integrated by using an external wrapper, at the moment we are working on the inclusion of the Q5Cost routines in the package itself.

Another activity in place regards Molekel [19], a free access molecular visualization program being developed by the Swiss National Supercomputing Centre (CSCS). The inclusion of a visualization program is valuable since many of the programs in the set do not have any graphical interface. For including Molekel in the workflow the strategy was to write a converted to the OpenBabel format [START_REF] Openbabel | The Blue Obelisk --Interoperability in Chemical Informatics[END_REF], a format that allows an immediate connection not only to Molekel but with several other environments that are already integrated. The Molcas [START_REF] Andersson | Molcas version 6.2[END_REF] wrapper is in good shape: at the moment, everything can be translated in Q5Cost format but the wave function expansion which will be the next improvement. A wrapper was also written for the QMC=Chem program, which is a Quantum Monte Carlo program developed at LCPQ, Toulouse.

Conclusion

The design and development of the new 1.0 version of the Q5Cost data format and library has been presented. The format and the library have been designed to facilitate the exchange of information between different Quantum Chemistry codes, thus strongly enhancing interoperability between codes. The ultimate goal of the reported activity is the construction of a grid-based distributed meta-laboratory, where multiple QC codes can be used together in a workflow finalized to solve a chemical problem.

The new features of the version 1.0 of the library have been presented: Geometry, Basis set and Wave Function storage, in addition to new tools in support of developers and end users, such as the binding to programming languages different than Fortran 90. A brief analysis of the wrappers and the codes which have been interfaced up to now has been presented too. The data format and the library appear as valuable tools to make the inter code communication possible and easy; some preliminary applications proved them to be effective, usable and also efficient, with respect to other Input/Output technologies.

The inclusion of the new features, in particular the Wave Function storage, is expected to allow a more widespread diffusion of Q5Cost, with the possibility to face new original scientific problems.

Glossary

Fig. 1 .

 1 Fig. 1. Schema of the Q5Cost data formatThe green boxes are containers (groups in HDF5 terminology), and the orange boxes are data (data sets), while the metadata (attributes describing the data) are listed in the white labels. The light green boxes are the new entities included in the 1.0 version of the data format.Originally the Q5Cost data format was designed to handle only the large binary data coming from Quantum Chemistry, together with the information (metadata) needed for describing them. The data structures described by the format were basically integrals, on the Atomic and Molecular orbitals, the accompanying metadata entities like electrons number, symmetry indices, nuclear energy, molecular orbital labels, but also coefficients needed to define symmetry adapted orbitals or to transform the Atomic orbitals into Molecular orbitals.Other small pieces of information, mainly ASCII coded in conventional QC data formats, are needed to define a QC system. They can be classified into three classes: Geometry, Basis set, spatial Symmetry of the system under investigation. In the original version of the data format those data were represented using an XML based language in a separated file[START_REF] Borini | [END_REF]. Starting from this version, we decided to merge them into the Q5Cost data format, in order to ensure a unique coherent structure grouping all the specifications on the chemical system. This is a new feature of the format that now includes the whole information about a molecular system.Data represented into the Q5Cost format belong to two different categories:

Fig. 2 .

 2 Fig. 2. Data structure for defining a Slater determinant based on the electrons in the system

Fig. 3 .

 3 Fig. 3. Data structure for defining a Slater determinant based on the MO in the system

-

  CAS = Complete Active Space -RAS = Restricted Active Space -CAS-SCF = Complete Active Space Self-Consistent Field -SCF = Self-Consistent Field -HF = Hartree-Fock -WF = Wave Function -CI = Configuration Interaction -SR = Single Reference -MR = Multi Reference -FCI = Full Configuration Interaction -O = Occupied Orbitals -V = Virtual Orbitals -A = Active Orbitals -I = Inactive Orbitals -CC = Coupled Cluster -C-CI = Contracted Configuration Interaction -DMRG = Density-Matrix Renormalization Group

Acknowledgments

Support from COST in Chemistry action D37 is gratefully acknowledged. This work has been also financed by the Universities of Bologna and Ferrara and by the Italian Ministry of Research and University MIUR under the project PRIN 2006 "Molecular Quantum Mechanics: Computational Methods and Analysis of Novel Phenomena".

The French CNR is gratefully acknowledged too for partial funding.


Julien Francken
email: julien.francken@arcelormittal.com

Didier Maquin

José Ragot

Bertrand Bèle
email: bertrand.bele@arcelormittal.com

Simultaneous data reconciliation and parameter estimation. Application to a basic oxygen furnace

In the steel industry, the determination of the control system set-points of batch processes is a common problem. It consists in adjusting the set-points in order to reach the given product specifications thanks to a process model. Small changes in operating conditions may impact final product quality. This is particularly true for the Basic Oxygen Furnace (BOF) where the information collected during a specific batch serves to adjust the set-points of the next batch. For being able to control that type of process, measurements must be made coherent and it may be convenient to use data reconciliation procedure. The proposed paper describes a method allowing simultaneous data reconciliation and model parameter estimation. Parameter estimation results can either be used to update the process model or to detect abnormal parameter variations due, e.g. to fouling, corrosion, degradation of parts of the process.

I. INTRODUCTION

A. Context

In the steel industry, the determination of the control system set-points of batch processes is a common problem. It consists in adjusting the set-points in order to reach the given product specifications thanks to a process model. Batch processes are characterized by prescribed processing of raw materials into refined products. The objective is to produce products of the desired quality. However, small changes in operating conditions may impact final product quality. Moreover, batch processes with the same trajectory generally exhibit some degree of batch-to-batch variation. This is particularly true for the Basic Oxygen Furnace (BOF) where the information collected during a specific batch serves to adjust the set-points of the next batch. Usually measurements collected from an industrial process contain random errors, systematic biases or gross errors due to the quality of the instrumentation chain and its environment. Based on first principles (mass and energy balances, chemical equilibrium relationships, etc), modelfiltering approach, commonly called data reconciliation may be used to provide coherent data. However, the main hypothesis underlying these methods is that the model describing the system is perfectly known which is not the case in real world. The proposed paper describes a method allowing simultaneous data reconciliation and model parameter estimation.

B. The Basic Oxygen Furnace

Basic oxygen furnace is a method of primary steelmaking in which carbon-rich molten iron is made into refined steel. The vast majority of steel manufactured in the world is produced using the basic oxygen furnace. Modern furnaces will take a charge of iron of up to 350 tons and convert it into steel in less than 40 minutes. By blowing oxygen through molten pig iron, the carbon content of the alloy is lowered and changes the material into low-carbon steel. The basic oxygen furnace actions are scheduled as follows: molten iron from a blast furnace is poured into a large refractory-lined container. Then the container is filled with the required ingredients which quantities are computed thanks to a setup charge balance model. As the required thermal energy is produced during the process, it is relevant to maintain a proper charge balance. A water-cooled lance is so lowered down into the vessel to blow pure oxygen onto the steel and iron, igniting the carbon dissolved in the steel and burning it to form carbon oxide, causing the temperature to rise to about 1700 • C. This melts the scrap, lowers the carbon content of the molten iron and helps remove unwanted chemical elements. Other fluid components are added to form slag which absorbs impurities of the steelmaking process and emulsifies to make easier the refining process. At the end of the blowing cycle, the temperature is measured and samples are taken to give a chemical analysis. The main functions of the BOF are to decarburize and remove phosphorus from the hot metal, and to optimize the steel temperature so that any further treatments prior to casting can be performed with minimal reheating or cooling of the steel. These actions are essentially performed by blowing oxygen with high pressure on hot metal and adding scraps and iron ore. The exothermic oxidation reactions that occur during BOF process generate a lot of heat energy, more than is necessary to attain the target steel temperature. This extra heat is used to melt scrap and/or iron ore additions.

For the BOF, the setup problem consists to determine the quantity of iron ore to add and the oxygen volume to blow to reaching the target of Carbon rate in steel and the temperature defined for each heat (batch) by the given product specification.

This task is done by a static charge computation based on a model formed by comprehensive heat and mass balances. These balances between input and output products describe the system. For example, a simplified iron mass balance equation could be as follows:

[%F e] Scp Q Scp + [%F e] HM Q HM -[%F e] Stl Q Stl = 0 (1)

C. Main objectives and proposed method

Data reconciliation and parameter estimation are fundamental components to real time optimization of industrial processes. Great efforts have been made to develop models for various processes. However, as mentioned by Rolandi and Romagnoli [START_REF] Rolandi | Simultaneous dynamic validation/identification of mechanistic process models and reconciliation of industrial process data[END_REF], physical and chemical process phenomena are complex and difficult to model conceptually and mathematically. Indeed, thermodynamic and transport properties and reaction rates are difficult to characterize experimentally and are subject to parametric uncertainties. Moreover, even when measurements are sufficiently numerous, raw processinstrumentation data are also corrupted by different errors. Indeed, data sets taken at a certain operation point are not model-consistent. Online optimization provides a mean to run a process near its optimum operating condition by providing real-time computed optimal set-points to the control system [START_REF] Faber | An optimization framework for parameter estimation of largscale systems[END_REF]. However, the implemented algorithms need to be fed by coherent data. The more classical way to deal with that problem consists to reconcile process data based on the minimization of measurement errors subject to satisfying the model constraints and next estimate the model parameters from these reconciled values. An alternative to this twostep approach consists to develop simultaneous strategies for data reconciliation and parameter estimation [START_REF] Rod | Iterative estimation of model parameters when measurements of all variables are subject to error[END_REF], [START_REF] Tjoa | Simultaneous strategies for data reconciliation and gross error detection of nonlinear systems[END_REF], [START_REF] Joe | A robust strategy for joint data reconciliation and parameter estimation[END_REF]. The proposed paper is dedicated to the presentation of such method.

On a mathematical point of view, data reconciliation is generally based on the assumption that the measurement errors have Gaussian probability density function (pdf) with zero mean. As some model parameters are uncertain, they are also modeled as random Gaussian variables.

Simultaneous state variable and parameter estimates are obtained using a maximum likelihood estimation approach applied on a sliding time window of observation data. To make less sensitive the parameter estimation to measurement uncertainties, the parameters are considered as constant value on a given time window.

Parameter estimation results can either be used to update process model or to detect abnormal variation due, e.g. to fouling, corrosion, degradation of parts of the process.

II. SIMULTANEOUS DATA RECONCILIATION AND PARAMETER ESTIMATION

The proposed paper describes a method allowing simultaneous data reconciliation and model parameter estimation in order to provide model-consistent data. In real process, model parameters are subjected to smooth evolution/degradation, that's why simultaneous state variable and parameter estimates will be extracted from a sliding time window of observation data. For each observation data window, parameters are assumed to be constant whereas state variables evolve around different operating points. This computation artifice will make less sensitive the parameter estimates to measurement uncertainties.

The method of simultaneous data reconciliation and parameter estimation on a given window of length N is firstly presented and the corresponding algorithm is outlined. Next, the practical implementation using a sliding window is described A. The process model Consider the system described by:

• a set of nonlinear equations in relation to a vector of state variables x * i and a parameter vector with true value a

* F (x * i , a * ) = 0, i ∈ [1, N ] (2)
x * i ∈ IR v , a * ∈ IR p , F : IR v×p → IR n

• observation equations

x i = x * i + ε xi , i ∈ [1, N ] (3)
• a priori knowledge of the parameters

a = a * + ε a (4)
The x i measurements of the true values x * i are available for the N realizations. For the parameter, the a priori knowledge is formulated as a kind of "observation equation". The parameter is defined on the basis of a nominal value and a "pseudonoise". This form allows to express in the same manner, the state variable and parameters probability density functions, given as follows:

p xi = 1 (2π) v/2 |V | 1/2 exp - 1 2 (x * i -x i ) T V -1 (x * i -x i ) (5)
p a = 1 (2π) p/2 |W | 1/2 exp - 1 2 (a * -a) T W -1 (a * -a) (6)
with V and W the variance-covariance matrices of the state variables and parameters respectively. The measurement and parameter error distributions are then assumed independent random variables.

B. State and parameter estimation method

The problem issue concerns the estimation of the true value of the state variables and the model parameters based on the knowledge of the measurements (3) on the given horizon [1, N ], the a priori knowledge (4) of the parameter and the process model [START_REF] Arora | Redescending estimators for data reconciliation and parameter estimation[END_REF]. The estimation can be performed by the maximum likelihood principle. Taking into account that measurement errors are independent, the likelihood function V is the product of the probability density functions:

V = N i=1 p xi p a (7)
Estimations xi , â of the true values maximize the likelihood function V in relation to x * i and a * with respect to the constraint (2) satisfaction.

Let Φ be the Lagrange function relative to the log-likelihood function and the constraint that have to be satisfied:

Φ = ln V + N i=1 λ T i F (x * i , a * ) (8) with λ i ∈ IR n .
Genrally, this nonlinear optimisation problem must be solved using an iterative algorithm solution. In this case, instead of a global system linearization arised from Lagrange regularity conditions, it is advisable to previously linearize the constraint equations.

For all the realizations i ∈ [1, N ], at the iteration j, we assume that a solution pair {x i,j , âj } is available. From this solution pair, estimates can be improved based on a first order Taylor series expansion in the neighborhood of this solution pair. Assuming the convergence of the estimation process, the successive solution pairs constitute a convergent series toward the problem solution.

At the iteration j + 1, for a solution pair {x i,j+1 , âj+1 }, in the neighborhood of {x i,j , âj }, the constraint first order Taylor series is:

F (x i,j+1 , âj+1 ) =F (x i,j , âj ) + G x (x i,j , âj )(x i,j+1 -xi,j ) + G a (x i,j , âj )(â j+1 -âj ) (9)
with the constraint gradients as follows:

G x (x i,j , âj ) = ∂F (x * i , a * ) ∂x * T i x * i =xi,j ,a * =âj (10a) G a (x i,j , âj ) = ∂F (x * i , a * ) ∂a * T x * i =xi,j ,a * =âj (10b)
For more readability, these last expressions will be denoted F i , G ix , G ia omitting the dependancy with regard to the quantities at iteration index j. At iteration j + 1, estimation problem reduces to the search of the Φ j+1 Lagrange function extremum in relation to x * i,j+1 , a * j+1 , and λ i :

Φ j+1 = Φ 0 - 1 2 (a * j+1 -a) T W -1 (a * j+1 -a) - 1 2 N i=1 (x * i,j+1 -x i ) T V -1 (x * i,j+1 -x i ) (11)
+ N i=1 λ T i F i + G ix (x * i,j+1 -xi,j ) + G ia (a * j+1 -âj )
where Φ 0 is the constant component of Lagrange function.

This Lagrange function presents an extremum for x * i,j+1 = xi,j+1 and a * j+1 = âj+1 if:

∂Φ j+1 ∂x * i,j+1 = 0 ∂Φ j+1 ∂a * j+1 = 0 (12) ∂Φ j+1 ∂λ i = 0
Equation ( 12) can be transformed into the following system:

-V -1 (x i,j+1 -x i ) + G T ix λ i = 0 (13) -W -1 (â j+1 -a) + N i=1 G T ia λ i = 0 (14)
F i + G ix (x i,j+1 -xi,j ) + G ia (â j+1 -âj ) = 0 (15)
From ( 13), ( 14) and (15), one can deduce:

xi,j+1 = x i + V G T ix λ i (16) âj+1 = a + W N i=1 G T ia λ i (17) G ix xi,j+1 = G ix xi,j -F i -G ia (â j+1 -âj ) (18)
The Lagrange parameters can be deduced from ( 16) and (18):

λ i = (G ix V G T ix ) -1 (G ix (x i,j -x i ) -F i -G ia (â j+1 -âj )) (19) Substituting (19) in (17) gives: âj+1 = a + W N i=1 S i (G ix (x i,j -x i ) -F i -G ia (â j+1 -âj )) (20)
with S i defined by:

S i = G T ia (G ix V G T ix ) -1 (21)
Let us define R as follows:

R = I + W N i=1 G T ia (G ix V G T ix ) -1 G ia (22)
Assuming that the matrix R is regular,

âj+1 =R -1 W N i=1 S i (G ix (x i,j -x i ) -F i + G ia âj ) + R -1 a (23)
Finally the expression of the state estimation is obtained from ( 16) and (19):

xi,j+1 = x i + V G T ix (G ix V G T ix ) -1 G ix (x i,j -x i ) -F i -G ia (â j+1 -âj ) (24)
The estimations xi and â reliant to x i and a are obtained using successive iterations provided algorithm convergence (the analysis of this latter can be based on the norms of the gradient matrices G ix and G ia ).

C. Algorithm

For an observation data window of length N , an algorithm based on the method developped in the previous section could be formulated in the following way:

• Initialization j = 0, ∀i ∈ [1; N ], xi,j = x i and âj = a • Repeat
Compute F i , G ix and G ia using ( 9), (10a) and (10b). Compute âj+1 using (23) Compute xi,j+1 using (24) j = j + 1 until norms of G ix and G ia are below a given threshold.

• Update the parameter and variable estimates

â = âj+1 xi = xi,j+1 , ∀i ∈ [1; N ]
This algorithm gives coherent parameter and state estimates for a given observation data window of length N . Let us now examine the practical implementation of this algorithm.

D. Practical implementation -Estimation on a sliding window

In order to be able to monitor the time evolution of the model parameters, the proposed algorithm is implemented using a sliding window. Let us consider a first observation window of length N . Let us denote a (N ) the a priori nominal value of the parameter vector (given by the user) for this observation window. From the knowledge of X 1 , . . . , X N and a (N ) , the proposed algorithm provides the estimates X1 , . . . , XN and â(N) . The observation window is then slid from one observation. Based on the knowledge of the measurements X 2 , . . . , X N +1 and the previously estimated value of parameter vector â(N) which serves as the nominal value for this new observation data window, the algorithm is used to provide the estimates XN+1 and â(N+1) . This process, which is depicted on figure 1 is re-iterated based on the availability of new sets of measurements along the time.

As previously said, the determination of the control system set-points of batch processes in the steel industry is difficult to solve. Indeed, small changes in the operating conditions

Time 1 N 2 N+1 3 N+2 X 1 , . . . , X N , a (N ) → X1 , . . . , XN , â(N) X 2 , . . . , X N +1 , â(N) → XN+1 , â(N+1) X 3 , . . . , X N +2 , â(N+1) → XN+2 , â(N+2)
Fig. 1. Principle of estimation on a sliding window may impact final product quality and some degree of batch-tobatch variation exists. Therefore, model adaptation algorithm are frequently implemented for that kind of processes. Clearly, the performances of the existing model adaptation algorithm will be enhanced when fed by the coherent estimated values provided by the proposed method.

III. EXAMPLE: SIMULTANEOUS STATE AND PARAMETER

ESTIMATIONS OF A SIMPLIFIED BOF PROCESS The proposed simultaneous state and parameter estimation method has been applied to a simplified nonlinear process model of a basic oxygen furnace (BOF).

A. The process model

Due to the random nature of chemical reaction in the non homogeneous slag or to the wear-out refractory-lined container or other many causes, the process cannot accurately be described with equilibrium balances. Therefore some of the model parameters need to be estimated without forget to use model consistent datas.

As described in section I-A, the model is a set of mass and heat balances. Let the nonlinear equation system be a simplified model of basic oxygen furnace described as next:

0.5x * 1 + (-3 + x * 2 )x * 3 + (a * 1 -x * 4 )x * 5 = 0 3x * 1 + (0.25x * 2 x * 4 -x * 5 )x * 3 + 9 = 0 x * 1 -0.5x * 2 x * 3 + x * 4 + a * 2 x * 5 -1 = 0
The first equation represents the iron mass balance where x 1 , x 3 , x 5 are material quantities, x 2 and x 4 iron mass percentage and a 1 a parameter which evolves along the time. The second equation describes a heat balance and the third another chemical element mass balance.

We have at our disposal state variable measurements x i for 1000 realizations. Their accuracies (constant for the whole realizations) are given in the table I. For that example, it is desired to monitor the values of the two parameters a * 1 and a * 2 , the others being assumed constant. A priori knowledge on a 1 and a 2 parameters are gathered in table II.

B. Results

To validate the proposed method, a data base composed of measurement data and true values of evolving parameters have been created for 1000 realizations. The simultaneous 2 and 3 when the observation data window length is fixed to N = 20, parameters estimates are near to the true value. Figure 4 shows a zoom (realizations 50 to 500) on the estimation of parameter a 2 when using different lengths for the sliding time window. The filtering capacity of the proposed estimator is pointed out. The more the window length is bigger, the more the parameter estimation filter the noise measurement, but, in the same time, the more the estimation is delayed. Figures 5,[START_REF] Rolandi | Simultaneous dynamic validation/identification of mechanistic process models and reconciliation of industrial process data[END_REF][START_REF] Tjoa | Simultaneous strategies for data reconciliation and gross error detection of nonlinear systems[END_REF]8,9 show the measurements of the state variables and their estimations. As the measurements are unbiased, the estimations are closed to the measurements.

IV. CONCLUSION

In this paper, a general methodology for simultaneous data reconciliation and parameter estimation for nonlinear model has been proposed. Data reconciliation based on balances was performed to obtain model-consistent measurement data and simultaneously parameter estimates. The use of a sliding window on which parameters are considered as constant values allows the desensitization of their estimations with regard measurement uncertainties. It is interesting to note that the knowledge of parameter distribution errors is not a limiting factor, indeed it would be easy to modify the algorithm in order to run without any parameter knowledge. The next step would be to include some robustness properties with regard gross measurement errors using a so-called contaminated distribution in the maximum likelihood estimation [START_REF] Arora | Redescending estimators for data reconciliation and parameter estimation[END_REF], [START_REF] Alhaj-Dibo | Data reconciliation: a robust approach using contaminated distribution[END_REF]. Finally the proposed method must be evaluated on real process model with real data in terms of process model adjustment along the successive heat to observe the impact of the setpoints adjustment on the successive batch processes.

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. a 2 parameter estimation

Fig. 6 .

 6 Fig. 6. x 2 state estimation and measurement

Fig. 7 .

 7 Fig. 7. x 3 state estimation and measurement

Fig. 9 .

 9 Fig. 9. x 5 state estimation and measurement

  with • [%F e] Scp : the iron mass percentage in scraps • [%F e] HM : the iron mass percentage in hot metal • [%F e] Stl : the iron mass percentage in steel • Q Scp : the scrap quantity to convert • Q HM : the hot metal quantity to convert • Q Stl : the converted steel quantity

	For each heat, the system state variables have a different
	operating point.

  Fig.2. a 1 parameter estimation state and parameter estimation algorithm described above was performed based on the measurement data from the whole realizations, the knowledge of the model, and only parameter nominal values given in table II. Because data proceed from created data base, we can easily observe the performance of the parameter estimation by matching parameter estimates with the true values. Results are given in figure

		2.8
		2.6
		2.4
		2.2
	value	2
		1.8
		1.6
		1.4	True value
			Estimation with horizon N=20
		0	200	400	600	800	1000
				realization