Creating a template file

A general template file for a Runge-Kutta function subroutine for the n body problem is detailed in this section. This involves evaluating the first derivatives of the Hamiltonian. The subroutine below is written so that derivatives with respect to components of position are evaluated first, followed by derivatives of components of momentum.

In[4]:= !!rksub.mf

        subroutine evalf(f,q,p)
        implicit double precision (a-h,o-z)
        dimension q(<* 

(***** Input degrees of freedom (dimension/2) *****)

                      d = 3;

        (***** End of required input *****)

(* Ensure that arrays are protected form N during translation. *)

SetOptions[FortranAssign,AssignToArray->{q,p}];

(* Automatically evaluate variables and n-body Hamiltonian *)

qvars = Array[q,d];    pvars = Array[p,d];
H = pvars.pvars/2 + 1/Sqrt[qvars.qvars];
d *>),p(<* d *>),f(<* 2 d *>)
<* (***** Calculate derivatives of the Hamiltonian. *****)

FortranAssign[
    f, 
    Flatten[{Outer[D,{H},pvars],- Outer[D,{H},qvars]}]
] *>
        return
        end

Appropriate array dimensions are specified by Mathematica. This means that precise control can be exercised over the generated FORTRAN subroutine, by allocating only as much stack space as the problem dictates. This is useful for large scale problems such as molecular dynamics simulations.

An attraction of this approach is that it may be generalised to account for any functional form of Hamiltonian. In order to do so, the Hamiltonian must be specified by the user in the template file, along with the number of degrees of freedom. Mathematica can then be used to Splice the results into a file containing FORTRAN source code for the problem.

In[5]:= Splice["rksub.mf",FormatType->OutputForm];

In[6]:= !!rksub.f

        subroutine evalf(f,q,p)
        implicit double precision (a-h,o-z)
        dimension q(3),p(3),f(6)
        f(1)=p(1)
        f(2)=p(2)
        f(3)=p(3)
        f(4)=q(1)/sqrt((q(1)**2+q(2)**2+q(3)**2)**3)
        f(5)=q(2)/sqrt((q(1)**2+q(2)**2+q(3)**2)**3)
        f(6)=q(3)/sqrt((q(1)**2+q(2)**2+q(3)**2)**3)
        return
        end