Code generation with Yacas

[JOEL LONTSI]

Hi Grzegorz,

I discover Yacas because I am currently interested in Automatic code generation for two-electron repulsions integrals in Quantum Chemistry. I came across your paper "Automatic Code Generation for Quantum Chemistry Applications". I would like to try the simplest of the examples described therein namely, 'The SVWN exchange-correlation functional'. I have installed yacas and do understand the Yacas implementation of this functional. My question now is: "Which yacas function or command should be use to generate C or Fortran code from the SVWN function described in the paper. I read the documentation from yacas.org but I could not find any such function or command.

Cheers,

Joel

[Grzegorz Mazur]

Hi Joel,

We're beginning a new academic year just after we moved the department to the new campus, so things are a bit hectic now. I'll answer very shortly now to give you a way to start doing anything at all and I'll follow up in next few days with some more advanced stuff.

Let's start with defining the Slater exchange functional:

In>  XSXSlater(rhoa,rhob) := -(3. / 4.) * (2. * (3. / Pi))^(1./3.) * (rhoa^(4./3.) + rhob^(4./3.))

Out> True

Get the expression in yacas:

In> XSlater(rhoa,rhob)

Out> -((rhoa^(4./3.)+rhob^(4./3.))*3.*(6./Pi)^(1./3.))/4.

Get the same expression in C:

In> CForm(XSlater(rhoa,rhob))

Out> " - ( ( pow(rhoa, 4. / 3.) + pow(rhob, 4. / 3.))  * 3. * pow(6. / Pi, 1. / 3.))  / 4."

Calculate the exchange potential (with respect to alpha density):

In> D(rhoa)XSlater(rhoa,rhob)

Out> -(3.*(6./Pi)^(1./3.)*4.*rhoa^(1./3.))/12.

And translate it to C:

In> CForm(D(rhoa)XSlater(rhoa,rhob))

Out> " - ( 3. * pow(6. / Pi, 1. / 3.) * 4. * pow(rhoa, 1. / 3.))  / 12."

This should give you the basic idea about how the automatic generation works. Next time I'll try to cover code templates, common subexpression elimination and explicit unfolding of recursively defined formulas (the latter is what you may need for generating code for integrals).

Cheers,

Grzesiek