3.5.1 A toy problem

We use a very simple potential for two particles, using VdW and Coulomb forces only.

We want to generate a program, that will just calculate the potential, its derivatives (i.e. the gradient) and its second derivatives (the Hessian).

We'll do this interactively in the following interactive exercise :

The window below is initialized with Maple commands to generate C code for the calculation of a van der Waals/Coulomb potential. You may edit the commands. Press "submit" to send them to the server. Results will be displayed in a new window. Plot results will be represented by a button in the result window. Click the button to open the plot window.

Cf. the Maple front-end for more information on how to use this Maple-interface.

The energy due to the VdW-forces is defined to be $\frac{a}{\rho^{12}}%% +\frac{b}{\rho^{6}},$ and we know, that it has a minimum of -1 at a distance $\rho=3.4$

This means we'll have two conditions: The derivative of the VdW-force with respect to $\rho$ has to be 0 at $\rho=3.4,$ because it is a minimum. And the VdW-force at $\rho=3.4$ has to be -1. See top part of the Maple activity :

##########################################
# Computation of the van der Waals force #
##########################################
vdw := a/rho^12 + b/rho^6;
# Minimum = -1 at 3.4
eqns := subs( rho = 34/10,{ diff( vdw, rho ) = 0, vdw = -1} );
solution:=solve( eqns, {a,b} );
vdw := subs(solution,vdw);
plot( vdw, rho=2..6, y=-1.5..2 );

Finally we have to find the coefficients $a$ and $b$, substitute them and get the exact expression for the VdW-force. Should you be in doubt, whether this result is correct, you can always plot it. (Compare Plot 1 in the result of the Maple activity above !)

We can convince ourselves that our result is correct: The minimum is indeed at $\rho=3.4$ and has value -1 and the shape of the function is typical for a van der Waals-force. So with this plot we get confidence, that our formula is correct.

Then we do the same for the Coulomb-forces, i.e the charges. For attraction, i.e. different charges, the constant is $k=100$, see the second block of the Maple-activity:

#####################################
# Computation of the Coulumb force  #
#####################################
Coul := Ch1*Ch2*k/rho;
Coulab := subs( Ch1=1, Ch2=-1, k=100, Coul );
Pot := vdw+Coulab;
plot( Pot, rho=2..6, y=-40..50 );

The total potential, i.e. the sum of the VdW and Coulomb potential is plotted in Plot 2 of the above results of our Maple-run!

The two atoms have different charges, so they attract each other, so their minimal energy position should be closer together compared with the VdW-forces alone. And indeed, it is closer together at approx. $\rho=3.$

The third block

###########################################################
# Computing the total potential, its gradient and hessian #
###########################################################
Pot := subs( rho = ((x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2 ) ^ (1/2), Pot );
Vars := [op(indets(Pot,name))];
res := poten=Pot:
for i to 6 do
    res := res, Grad[i] = diff(Pot,Vars[i]);
    for j from i to 6 do
        res := res, Hess[i,j] = 
                   diff(Pot,Vars[i],Vars[j]) od
    od:
res := [res]:  codegen[cost](res);

shows the replacement of coordinates and distances for rho. The distance $\rho,$ is not a simple distance, it is $\rho=\sqrt{\left( x_{2}-x_{1}\right) ^{2}+\left( y_{2}-y_{1}\right) ^{2}+\left( z_{2}%% -z_{1}\right) ^{2}}.$ So when we substitute this, the formulas become larger.

The variables in this problem are $x_{1},y_{1,}z_{1,}x_{2,}y_{2,}z_{2}$ and for doing minimization, we have to compute the derivatives of this function. So we compute the gradient, which is the derivative of the potential with regard to each variable, and we compute the Hessian, which is the second derivative.

Now when we take this function (Pot), which is fairly small, and compute the derivatives twice with respect to all variables, the result is very large. It takes 500 subtractions, more than 1100 multiplications, more than 100 divisions and 100 additions, all in all about 2000 operations to compute the expression res(ult) with all its derivatives. Remember: This is still our toy example!

The fourth block shows the optimization:

############################
# Optimization of the code #
############################
opt := codegen[optimize](res,tryhard):  codegen[cost](opt);
codegen[C]( [opt] );

And here is, where optimization really pays, because all of these expressions are highly repetitive. The same subexpression appears again and again and again. And so, if we optimize symbolically, then instead of 2000 operations we'll have only 150 operations. Naturally we have more assignments now, since we use more temporary variables.

All such programs typically look the same: The first part is assigning a lot of temporary variables. These temporary assignments are typically very simple and they are reused.

In this case we use two criteria for optimization:

• Any repeated sub-expression is assigned to a temporary variable and hence computed only once.

• Any temporary variable used only once is eliminated by substitution.

Every temporary variable will be used at least twice, and no expression that appears on the right hand side will appear twice, because if it appeared twice, we could have stored it in a temporary variable.

Summary:

1. We started from formulas (from the theory).

2. We constructed potentials.

3. We constructed a program, which was optimized in two steps.

What confidence do we have, that this program gives the original expression?

The beauty of working in a symbolic CAS is that you can get your optimized function and execute it symbolically. This is shown in the last block of the Maple-code:

#########################################
# Symbolic execution of the program and #
# equality testing with the original    #
#########################################
:= poten:
for i from nops([opt]) by -1 to 1 do
    t := subs(opt[i], t ) od:
t;
testeq( eval( subs( pow = ( (a,b)->a^b ), t)) = Pot );

To execute it symbolically we assign, sequentially each temporary variable its right-hand-side value. So we get in the end a potential function which should be the original one. We can shortcircuit this whole process of generating the program and we find out what is the result of the computation of the C-program.

The loop is a very compact way of representing all the assignments. Basically we use the assignments as substitutions, start with the original potential and substitute every temporary variable into it. At the end we get an expression which doesn't look exactly like the original one.

In particular, since it is a C-program, it contains the function pow, which means power, and it also has a different aspect. However, when we substiute pow for powering and evaluate everything and then test equivalence we get equality.