\documentclass[]{report}
\usepackage{fontspec}
\usepackage{listings}
\usepackage{amsmath}
\usepackage{xcolor}
\usepackage{hyperref}
\usepackage{graphicx}
\hypersetup{
    bookmarks=true,         % show bookmarks bar?
    unicode=false,          % non-Latin characters in Acrobat’s bookmarks
    pdftoolbar=true,        % show Acrobat’s toolbar?
    pdfmenubar=true,        % show Acrobat’s menu?
    pdffitwindow=false,     % window fit to page when opened
    pdfstartview={FitH},    % fits the width of the page to the window
    pdftitle={My title},    % title
    pdfauthor={Author},     % author
    pdfsubject={Subject},   % subject of the document
    pdfcreator={Creator},   % creator of the document
    pdfproducer={Producer}, % producer of the document
    pdfkeywords={keyword1} {key2} {key3}, % list of keywords
    pdfnewwindow=true,      % links in new window
    colorlinks=true,       % false: boxed links; true: colored links
    linkcolor=blue,          % color of internal links (change box color with linkbordercolor)
    citecolor=green,        % color of links to bibliography
    filecolor=magenta,      % color of file links
    urlcolor=cyan           % color of external links
}
\definecolor{listingsstringcolor}{rgb}{0,0.5,0}
\definecolor{listingskeywordcolor}{rgb}{0.5,0.5,0.0}
\definecolor{listingsbasiccolor}{rgb}{0.0,0.0,0.0}
% \definecolor{listingskeywordcolor}{rgb}{0.0,0.0,0.7}
\definecolor{listingsnumbercolor}{rgb}{0.0,0.0,1.0}
\definecolor{listingscommentcolor}{rgb}{0.4,0.4,0.4}
\definecolor{listingsbackgroundcolor}{rgb}{0.975,0.975,0.975}
\definecolor{listingsrulecolor}{rgb}{0.86,0.86,0.86}
\definecolor{listingsidentifiercolor}{rgb}{0.5,0.0,0.5}
% \definecolor{listingsvariablecolor}{rgb}{0.5,0.0,0.5}

\renewcommand{\thesection}{\arabic{section}}

\lstset {
  language=c++,
  numbers=none,
  breaklines=true,
  tabsize=2,
  backgroundcolor=\color{listingsbackgroundcolor},
  breakatwhitespace=false,         % sets if automatic breaks should only happen at whitespace
  breaklines=true,                 % sets automatic line breaking
  %   numbers=left,                    % where to put the line-numbers; possible values are (none, 
  %%left, right)
  numbersep=5pt,                   % how far the line-numbers are from the code
  frame=lrtb,                    % adds a frame around the code
  framexleftmargin=7pt,
  framexrightmargin=7pt,
  framextopmargin=7pt,
  framexbottommargin=7pt,
  xleftmargin=18pt,
  xrightmargin=18pt,
  rulecolor=\color{listingsrulecolor},
  tabsize=2,                       % sets default tabsize to 2 spaces
  literate={å}{{\aa}}1 {æ}{{\ae}}1 {ø}{{\oslash}}1,
  showstringspaces=false,
  captionpos=b,
  basicstyle=\color{listingsbasiccolor}\footnotesize\sffamily,
  keywordstyle=\color{listingskeywordcolor}\footnotesize\sffamily,
  stringstyle=\color{listingsstringcolor}\footnotesize\sffamily,
  commentstyle=\color{listingscommentcolor}\footnotesize\sffamily,
  numberstyle=\color{listingsnumbercolor}\footnotesize\sffamily,
  identifierstyle=\color{listingsidentifiercolor}\footnotesize\sffamily,
}

\newfontfamily\listingsfont[Scale=0.85]{Droid Sans Mono}
\lstset {
  basicstyle=\color{listingsbasiccolor}\footnotesize\listingsfont,
  keywordstyle=\color{listingskeywordcolor}\footnotesize\listingsfont,
  stringstyle=\color{listingsstringcolor}\footnotesize\listingsfont,
  commentstyle=\color{listingscommentcolor}\footnotesize\listingsfont,
  numberstyle=\color{listingsnumbercolor}\footnotesize\listingsfont,
  identifierstyle=\color{listingsidentifiercolor}\footnotesize\listingsfont,
}

%opening
\title{}
\author{}

\begin{document}

\chapter*{Week 11-12}

In the coming weeks, we are focusing on implementing a Hartree-Fock solver for Gaussian Type 
Orbitals (GTOs) in FYS4411. This will allow you to calculate the energies (and other observables) 
for much larger systems, such as complex molecules, or simply the molecular orbitals of systems 
such as H$_{2}$O:
%
\begin{center}
\includegraphics[width=0.6\textwidth]{mo_h2o.pdf}
\end{center}
%
This document is written to give you some helpful hints about how you could move forward, but the 
final sections are not yet complete, so expect some updates as we progress.

\section{Create class HartreeFockSolver}
Hopefully most of you have a working Hartree-Fock solver for atoms with hydrogenic basis functions.
Now we need to generalize this solver to larger systems. The first step is to create a class
HartreeFockSolver that may contain your currect code or act as a skeleton for the code to come.

\section{Create class Primitive}
To get you started we will give you a complete example of how class Primitive can be written.
The header file is:
\begin{lstlisting}
#ifndef PRIMITIVE_H
#define PRIMITIVE_H

#include <armadillo>

using namespace arma;

class Primitive
{
public:
    explicit Primitive(double weight, int xExponent,
                       int yExponent, int zExponent, double exponent,
                       vec nucleusPosition);

    double exponent() const;
    int zExponent() const;
    int yExponent() const;
    int xExponent() const;
    double weight() const;
    const vec& nucleusPosition() const;

private:
    double m_weight;
    int m_xExponent;
    int m_yExponent;
    int m_zExponent;
    double m_exponent;
    vec m_nucleusPosition;
};

#endif // PRIMITIVE_H
\end{lstlisting}
and the .cpp file is:
\begin{lstlisting}
#include "primitive.h"

Primitive::Primitive(double weight,
                     int xExponent, int yExponent, int zExponent, 
                     double exponent, vec nucleusPosition) :
    m_weight(weight),
    m_xExponent(xExponent),
    m_yExponent(yExponent),
    m_zExponent(zExponent),
    m_exponent(exponent),
    m_nucleusPosition(nucleusPosition)
{
		// If you haven't seen the above way of doing initilization, check out
    // http://www.informit.com/guides/content.aspx?g=cplusplus&seqNum=172
}

double Primitive::exponent() const
{
    return m_exponent;
}

int Primitive::zExponent() const
{
    return m_zExponent;
}

int Primitive::yExponent() const
{
    return m_yExponent;
}

int Primitive::xExponent() const
{
    return m_xExponent;
}

double Primitive::weight() const
{
    return m_weight;
}

const vec& Primitive::nucleusPosition() const
{
    return m_nucleusPosition;
}
\end{lstlisting}

\section{Create class Integrator}
Create the skeleton of class Integrator, which will calculate all the integrals. Each function in
the integrator should receive primitive objects as parameters.

\section{Write function overlap}
The simplest function in class Integrator is the overlap integral between two primitives:
\begin{equation}
S_{ab} = \langle G_a|G_b\rangle,
\end{equation}
where
\begin{equation}
\begin{split}
G_a & = G_{ikm}(a, \mathbf r_A) = x_A^i y_A^k z_A^m \exp(-a r_A^2) \\
G_b & = G_{jln}(b, \mathbf r_B) = x_B^j y_B^l z_B^n \exp(-b r_B^2).
\end{split}
\end{equation}
To do this integral you have to calculate the coefficients $E^{ij}_0$, $E^{kl}_0$ and $E^{mn}_0$
by recursion:
\begin{equation}
S_{ab} = E^{ij}_0 E^{kl}_0 E^{mn}_0\left(\frac{\pi}{p}\right)^{3/2},
\end{equation}
where
\begin{equation}
E^{i=0,j=0}_0 = \exp\left(-\frac{ab}{a+b}X_{AB}^2\right),
\end{equation}
and equivalently for the other two directions.

\subsection{Calculate simple overlap}
In the simplest case when both primitives are s-type Gaussians, that is, when all exponents
are equal to zero ($i=j=k=l=m=n=0$), you don't need to do any recursion. You can use this to
start doing simple tests with different primitives.
%TODO Give values to check against

\subsection{Write function setupHermiteCoefficients ($E^{ij}_t$-factors)}
To the general integrals you have to calculate all the Hermite coefficients ($E^{ij}_t$-factors),
see for example Helgaker's slides page 16.

\subsection{Calculate general overlap integrals}
Once you have the function setupHermiteCoefficients, you are ready to calculate the general overlap
integrals.

We have selected three values that you may test against using 
\href{http://www.uio.no/studier/emner/matnat/fys/FYS4411/v14/guides/unit-tests/index.html}{UnitTest++}:
\begin{lstlisting}
// PrimitiveA
a = 0.2;
weight = 1;
i = k = m = 0;
A = {1.2, 2.3, 3.4};

// PrimitiveB
b = 0.3;
weight = 1;
j = l = n = 0;
B = {-1.3, 1.4, -2.4};

Primitive primitiveA(weight,i,k,m,a,A);
Primitive primitiveB(weight,j,l,n,b,B);

CHECK_CLOSE(1.191723635809e-01, integrator.overlapIntegral(primitiveA, primitiveB), 1e-5);
\end{lstlisting}
\begin{lstlisting}
// PrimitiveA:
a = 0.2;
weight = 1;
i = m = 0;
k = 1;
A = {1.2, 2.3, 3.4};

// PrimitiveB:
b = 0.3;
weight = 1;
j = 0;
l = n = 1;
B = {-1.3, 1.4, -2.4};

Primitive primitiveA(weight,i,k,m,a,A);
Primitive primitiveB(weight,j,l,n,b,B);

CHECK_CLOSE(2.227321941537e-01, integrator.overlapIntegral(primitiveA, primitiveB), 1e-5);
\end{lstlisting}
\begin{lstlisting}
// PrimitiveA:
a = 0.2;
weight = 1;
i = m = 0;
k = 2;
A = {1.2, 2.3, 3.4};

// PrimitiveB:
b = 0.3;
weight = 1;
j = l = 1;
n = 0;
B = {-1.3, 1.4, -2.4};

Primitive primitiveA(weight,i,k,m,a,A);
Primitive primitiveB(weight,j,l,n,b,B);

CHECK_CLOSE(-7.329386373895e-02, integrator.overlapIntegral(primitiveA, primitiveB), 1e-5); 
\end{lstlisting}

\section{Write function kinetic}
The kinetic integrals are given as simple products of overlap integrals:
\begin{equation}
T_{ab} = -\frac{1}{2}(T_{ij}S_{kl}S_{mn} + S_{ij}T_{kl}S_{mn} + S_{ij}S_{kl}T_{mn}),
\end{equation}
where
\begin{equation}
T_{ij} = 4b^2S_{i,j+2} - 2b(2j+1)S_{ij} + j(j-1)S_{i,j-2}.
\end{equation}
Here are some tests that you may run to check that your kinetic integrator works. 
\begin{lstlisting}
// PrimitiveA
a = 0.2;
weight = 1;
i = k = m = 0;
A = {1.2, 2.3, 3.4};

// PrimitiveB
b = 0.3;
weight = 1;
j = l = n = 0;
B = {-1.3, 1.4, -2.4};

Primitive primitiveA(weight,i,k,m,a,A);
Primitive primitiveB(weight,j,l,n,b,B);

CHECK_CLOSE(-9.678702680582e-02, integrator.kineticIntegral(primitiveA, primitiveB), 1e-5); 
\end{lstlisting}
\begin{lstlisting}
// PrimitiveA:
a = 0.2;
weight = 1;
i = m = 0;
k = 1;
A = {1.2, 2.3, 3.4};

// PrimitiveB:
b = 0.3;
weight = 1;
j = 0;
l = n = 1;
B = {-1.3, 1.4, -2.4};

Primitive primitiveA(weight,i,k,m,a,A);
Primitive primitiveB(weight,j,l,n,b,B);

CHECK_CLOSE(-8.688217105502e-02, integrator.kineticIntegral(primitiveA, primitiveB), 1e-5); 
\end{lstlisting}
\begin{lstlisting}
// PrimitiveA:
a = 0.2;
weight = 1;
i = m = 0;
k = 2;
A = {1.2, 2.3, 3.4};

// PrimitiveB:
b = 0.3;
weight = 1;
j = l = 1;
n = 0;
B = {-1.3, 1.4, -2.4};

Primitive primitiveA(weight,i,k,m,a,A);
Primitive primitiveB(weight,j,l,n,b,B);

CHECK_CLOSE(-1.598401092187e-02, integrator.kineticIntegral(primitiveA, primitiveB), 1e-5); 
\end{lstlisting}
To obtain more tests, you may calculate some new integrals numerically using for instance Python, 
Matlab or your own C++ code.

\chapter*{Week 13-14}
\section{Create class BoysFunction}
For your convenience we have created a class BoysFunction that you can used. It is defined as:
\begin{lstlisting}
#ifndef BOYSFUNCTION_H
#define BOYSFUNCTION_H

#include <armadillo>

using namespace arma;


class BoysFunction
{
public:
    BoysFunction(int angMomMax);
    void setx(double x);
    double returnValue(int n);
private:
    double tabulated(int n, double x);
    double asymptotic(int n, double x);
    double factorial2(int n);

    mat m_Ftabulated;
    vec m_F;
    int m_nMax;
};

#endif // BOYSFUNCTION_H
\end{lstlisting}
with the following implementation:
\begin{lstlisting}
#include "boysfunction.h"

BoysFunction::BoysFunction(int angMomMax)
{
    m_nMax = 4*angMomMax;
    m_F = zeros<vec>(4*angMomMax+1);
    m_Ftabulated.load("boys_tabulated.dat");
}

double BoysFunction::returnValue(int n)
{
    return m_F(n);
}

void BoysFunction::setx(double x)
{
    if (x <= 50){
        m_F(m_nMax) = tabulated(m_nMax, x);
    } else {
        m_F(m_nMax) = asymptotic(m_nMax, x);
    }

    double ex = exp(-x);

    for(int n = m_nMax; n > 0; n--){
        m_F(n-1) = (2*x*m_F(n) + ex)/(2*n - 1);
    }
}

double BoysFunction::tabulated(int n, double x)
{
    int nxVals = m_Ftabulated.n_rows;
    double dx = 50.0/(nxVals - 1); // x-spacing in tabulated values
    int xIndex = int ((x + 0.5*dx)/dx);
    double xt = xIndex*dx;     // tabulated x-value
    double Dx = x - xt;        // difference between actual and tabulated x-value

    double value = 0;
    double factorial = 1;

    for(int k = 0; k < 7; k++){
        if(k != 0){
            factorial *= k;
        }
        value += m_Ftabulated(xIndex, n+k)*pow(-Dx,k)/factorial;
    }

    return value;
}

double BoysFunction::asymptotic(int n, double x)
{
    return factorial2(n)*sqrt(M_PI/(pow(x,2*n+1)))/(pow(2,n+1));
}

double BoysFunction::factorial2(int n)
{
    double value = 1;
    double i = 1;

    while(i < 2*n - 1){
        i += 2;
        value *= i;
    }
    return value;
}
\end{lstlisting}
This class uses a file named boys\_tabulated.dat which you can find \href{http://www.uio.no/studier/emner/matnat/fys/FYS4411/v14/guides/gto/boys_tabulated.dat}{here}.
Use this class in the following way:
\begin{lstlisting}
 BoysFunction boys(2); // Boys function for system with l <= 2
                       // This will allow calculations up to F_8(x)
 boys.setx(1.4); // Will evaluate the Boys function for x = 1.4
 cout << "F_0 (1.4) = " << boys.returnValue(0) << endl;
 cout << "F_3 (1.4) = " << boys.returnValue(3) << endl;
\end{lstlisting}


\section{Write function setupHermiteIntegrals ($R_{tuv}$-factors)}
This function should be implemented in class Integrator.
Make a function which calulates the Hermite integrals $R_{tuv}$ using the recursioin relations
from the slides. These integrals may be stored in an armadillo object:
\begin{lstlisting}
field<mat> R;
R.set_size(4);
for(int n = 0; n < 4; n++){
    R(n) = zeros(4, 4, 4);
} // R is now 4x4x4x4
\end{lstlisting}


\section{Write function nucleiElectronIntegral}
This function should be implemented in class Integrator.
Make a function that recieves two primitives and the nucleus position and returns the value of the integral.
Alternatively (as is shown in the test below) you may have functions such as setPrimitiveA etc. in the integrator class.
The result must be multiplied by the charge and coefficients of the primitives the function. Here are some tests you can use:
\begin{lstlisting}
 TEST(GTOnuclearAttractionIntegral)
{
    /*
     * test case:   electron-nucleus integral
     * max angular momentum:    2
     *
     * source:
     *      numerical integration
     *      python scripts
     * */

    Integrator integrator;

    rowvec posA = {1.2, 2.3, 3.4};
    rowvec posB = {-1.3, 1.4, -2.4};
    rowvec posC = {2.3, 0.9, 3.2};


    PrimitiveGTO primitiveA(0.2, 1.0);
    PrimitiveGTO primitiveB(0.3, 1.0);


    integrator.setNucleusPosition(posC);
    
    
    primitiveA.setPowers({0,0,0}); primitiveB.setPowers({0,0,0});
integrator.setPrimitiveA(primitiveA);integrator.setPrimitiveB(primitiveB);
    CHECK_CLOSE(2.788948987251e-02, integrator.nuclearAttractionIntegral(), 1e-5);

    primitiveA.setPowers({0,0,0}); primitiveB.setPowers({0,0,1});
integrator.setPrimitiveA(primitiveA);integrator.setPrimitiveB(primitiveB);
    CHECK_CLOSE(6.971203468743e-02, integrator.nuclearAttractionIntegral(), 1e-5);

    primitiveA.setPowers({0,0,0}); primitiveB.setPowers({0,0,2});
integrator.setPrimitiveA(primitiveA);integrator.setPrimitiveB(primitiveB);
    CHECK_CLOSE(2.024071525839e-01, integrator.nuclearAttractionIntegral(), 1e-5);

    primitiveA.setPowers({0,0,0}); primitiveB.setPowers({0,1,0});
integrator.setPrimitiveA(primitiveA);integrator.setPrimitiveB(primitiveB);
    CHECK_CLOSE(8.727033700014e-03, integrator.nuclearAttractionIntegral(), 1e-5);

    primitiveA.setPowers({0,0,0}); primitiveB.setPowers({0,1,1});
integrator.setPrimitiveA(primitiveA);integrator.setPrimitiveB(primitiveB);
    CHECK_CLOSE(2.134361291529e-02, integrator.nuclearAttractionIntegral(), 1e-5);

    primitiveA.setPowers({0,0,0}); primitiveB.setPowers({0,2,0});
integrator.setPrimitiveA(primitiveA);integrator.setPrimitiveB(primitiveB);
    CHECK_CLOSE(2.921666495443e-02, integrator.nuclearAttractionIntegral(), 1e-5);

    primitiveA.setPowers({2,0,0}); primitiveB.setPowers({0,2,0});
integrator.setPrimitiveA(primitiveA);integrator.setPrimitiveB(primitiveB);
    CHECK_CLOSE(8.419480232293e-02, integrator.nuclearAttractionIntegral(), 1e-5);

    primitiveA.setPowers({2,0,0}); primitiveB.setPowers({1,0,0});
integrator.setPrimitiveA(primitiveA);integrator.setPrimitiveB(primitiveB);
    CHECK_CLOSE(1.481688684288e-02, integrator.nuclearAttractionIntegral(), 1e-5);

    primitiveA.setPowers({2,0,0}); primitiveB.setPowers({1,0,1});
integrator.setPrimitiveA(primitiveA);integrator.setPrimitiveB(primitiveB);
    CHECK_CLOSE(3.878852644576e-02, integrator.nuclearAttractionIntegral(), 1e-5);

    primitiveA.setPowers({2,0,0}); primitiveB.setPowers({1,1,0});
integrator.setPrimitiveA(primitiveA);integrator.setPrimitiveB(primitiveB);
    CHECK_CLOSE(4.176920693786e-03, integrator.nuclearAttractionIntegral(), 1e-5);

    primitiveA.setPowers({2,0,0}); primitiveB.setPowers({2,0,0});
integrator.setPrimitiveA(primitiveA);integrator.setPrimitiveB(primitiveB);
    CHECK_CLOSE(6.422210627967e-02, integrator.nuclearAttractionIntegral(), 1e-5);

} 
\end{lstlisting}
This code will not compile as written, but should give you some values to check against.

\section{Write function electronElectronIntegral}
This function should be implemented in class Integrator. Make a function that receives four primitives
and returns the value of the integral.


\section{Create class ElectronSystem}

\section{Write function nucleiInteractionEnergy}


\chapter{Week 15-16}


\end{document}
