[C++] Schrodinger equation solver and nuclear density

[BRN]

Hi everybody!
I need to compute a C++ program for solve Schrodinger equation and calculate nuclear density.
My nucleus is made up of only neutrons immersed in a potential of a harmonic oscillator.
Schrodinger equation is:

$$[-\frac{\hbar^2}{2m}\triangledown^2+V_{HO}(r)]\psi=E\psi$$

with $$V_{HO}(r)=\frac{1}{2}m\omega^2r^2$$

At frist time I should use the HO eigenvalues formula

$$E_{nl}=\hbar\omega(2n+l+\frac{3}{2})$$

then, I should use a eigenvalues generator by shooting method.
Finally, my density equation is

$$n(r)=\frac{1}{4\pi r^2}\sum_\alpha(2j_\alpha+1)\psi_{\alpha}^2(r)$$

with $\alpha$ quantum number $n_{\alpha}, l{\alpha}, j_{\alpha}$

I have some problem to do it.
1) To solve the density, I need the solver to give me back a series of eigenfunctions for each value of $n$ and $l$. Right? But if I have the eigenvalues generated by the shooting method, I no longer have an equation that depends on $n$ and $l$, so what can I do?

2) The shooting method:

• I should provide two test eigenvalues, E1 and E2 and calculate the relative eigenfunctions.
• If they are in disagreement in $r_{max}$, then I define $E_t = (E1-E2) / 2$, calculate the relative eigenfunction and study the sign in $r_{max}$.
• If it is concordant, it replaces the eigenvalue between E1 and E2.
• One proceeds for further bisections until: $|\psi (r_{max}|<10^{-6}$

I have problem to implement this...

My code:
All parameters and constant are in "parameters.h"

//parameters.h

#pragma once
//#ifndef parameters_h
//#define parameters_h
#include<cmath>

namespace Parameters
{

const int NN=10; // Neutrons number
const int NP=10; // Protons number
const double mp= 1.6726219e-27 ; // Proton mass
const double mn= 1.6749273e-27 ; // Neutron mass
const double Rn=1e-9 ;

// Spin-Orbit potential parameters
//const double k0= ;
//const double r0= ;

// HO potential parameters
const double m=mn ;
const double f = 100;
//double r= ;

// Eigenvalues generator parameters
double eigenvalue1=0.93250 ;
double eigenvalue2=1.36256 ;
double error=10e8; //maybe you mean 1e-08?

// Other parameters
const unsigned int energyLevel = 1;
const int angularMomentum = 1;
const double PI = 4*atan(1);
const double hbar = 1.0545718*1e-34;
const double x_min=0; //Integration
const double x_max=3*Rn; //Interval

}

//#endif

Schrodinger solver class are "schroddy.h"

#pragma once

class InitialPot;

//hierarchy for eigenvalue objects

class Eigenvalues
{
public:
virtual double eigenvalue() const = 0;
virtual ~Eigenvalues();

virtual Eigenvalues* clone() const = 0;

private:

};class HarmonicEigenvalues: public Eigenvalues
{
public:
HarmonicEigenvalues(double omega, unsigned int n, int l);
double eigenvalue() const override;
//setEnergyLevel(unsigned int n);
//setAngularMomentum(int l);

Eigenvalues* clone() const override;

private:
HarmonicEigenvalues();
const double m_omega;

//if you want to be able to change quantum numbers within same object, implement set methods and remove const keyword
const unsigned int m_n;
const int m_l;
};class GenericEigenvalues: public Eigenvalues
{
public:
GenericEigenvalues (double eigval1, double eigval2, double eigenvalt, double err);
double eigenvalue() const override;

Eigenvalues* clone() const override;

private:
GenericEigenvalues();
const double m_eigval1;
const double m_eigval2;
const double m_eigenvalt;
const double m_err;
};//Shroedinger class + solverclass Schroddy
{
public:
Schroddy(const InitialPot& pot, const Eigenvalues& eigenval); //if you have more parameters to pass in, change this
double solveShroddyByRK(double x0, double x1, double psi0, double psiPrime0, unsigned long NSteps) const; //psiPrime0 is boundary condition on first derivative of eigenfunction

~Schroddy();

private:
const InitialPot* const m_pot;
const Eigenvalues* const m_eigenval;

};

and "schroddy.cpp"

//Class for Schrodinger equation solver
# include <math.h>
# include <stdlib.h>
# include <stdio.h>
# include "initpot.h"
# include "parameters.h"
# include "schroddy.h"
# include <cmath>
# include <vector>Eigenvalues::~Eigenvalues() {}

/*=======================================================================
* HO Eigenvalues generator
*=====================================================================*/

HarmonicEigenvalues::HarmonicEigenvalues(double omega, unsigned int n, int l): m_omega(omega), m_n(n), m_l(l) {}

double HarmonicEigenvalues::eigenvalue() const
{
return Parameters::hbar*m_omega*(2*m_n+m_l+(3/2));
}

Eigenvalues* HarmonicEigenvalues::clone() const
{
return new HarmonicEigenvalues(*this);
}/*========================================================================
* Eigenvalues generator by shooting method
*======================================================================*/

GenericEigenvalues::GenericEigenvalues (double eigval1, double eigval2, double eigenvalt, double err): m_eigval1(eigval1), m_eigval2(eigval2), m_eigenvalt(eigenvalt), m_err(err){}

double GenericEigenvalues::eigenvalue() const //this must return a double..
{
//m_eigval1=Parameters::eigenvalue1; //cannot assign to a const variable; don't you pass in these values to the constructor? why do it twice?
//m_eigval2=Parameters::eigenvalue2;
//m_err=Parameters::error;return 0.0;

}

Eigenvalues* GenericEigenvalues::clone() const
{
return new GenericEigenvalues(*this);
}

/*=========================================================================
* Schrodinger equation solver by Runge-Kutta method
*=======================================================================*/

Schroddy::Schroddy(const InitialPot& pot, const Eigenvalues& eigenval): m_pot(pot.clone()), m_eigenval(eigenval.clone()) {}

Schroddy::~Schroddy()
{
delete m_pot;
delete m_eigenval;
}

double Schroddy::solveShroddyByRK(double x0, double x1, double psi0, double psiPrime0, unsigned long NSteps) const
{
//implement runge-kutta here..
const double h = (x1 - x0)/NSteps;
const double factor = 2*Parameters::mn/(Parameters::hbar*Parameters::hbar);
const double eigenvalue = m_eigenval -> eigenvalue();

double runningX = x0, runningPsi = psi0, runningPsiPrime = psiPrime0;
std::vector<double> psiArray;
psiArray.push_back(psi0);

for (unsigned long i = 0; i < NSteps; ++i)
{

//compute decoupled RK factors
const double k1 = runningPsiPrime;
const double l1 = (m_pot -> potential(runningX) - eigenvalue)*runningPsi;
const double k2 = runningPsiPrime + h/2*l1;
const double l2 = (m_pot -> potential(runningX + h/2) - eigenvalue)*(runningPsi + h/2*k1);
const double k3 = runningPsiPrime + h/2*l2;
const double l3 = (m_pot -> potential(runningX + h/2) - eigenvalue)*(runningPsi + h/2*k2);
const double k4 = runningPsiPrime + h*l3;
const double l4 = (m_pot -> potential(runningX + h) - eigenvalue)*(runningPsi + h*k3);

//advance running variables and store intermediate results
runningPsi += h/6*factor*(k1 + 2*k2 + 2*k3 + k4);
runningPsiPrime =+ h/6*factor*(l1 + 2*l2 + 2*l3 + l4);
runningX += h;
psiArray.push_back(runningPsi);

}
//work out normalizatiion constant
double psiSquared = 0, normalPsi=0;
for (std::vector<double>::iterator it = psiArray.begin(); it != psiArray.end(); ++it)
psiSquared += (*it)*(*it); //derefence iterator at array's element and square (brackets are needed!)

const double scalar = psiSquared*h;

//return normalized eigenfunction (on interval [x0, x1]) and job done!
return normalPsi=runningPsi/sqrt(scalar);
}

Potential class are "initpot.h"

#pragma once

class InitialPot //pure virtual base class
{
public:
virtual double potential(double x) const = 0;
virtual ~InitialPot();

virtual InitialPot* clone() const = 0; //virtual constructor

private:

};

class HOPot: public InitialPot //derived class
{
public:
//=====================================================================
// HO potential
//=====================================================================
HOPot(double m, double omega);
double potential(double x) const override;

//double getOmega() const;
InitialPot* clone() const override;

private:
HOPot(); //same as before
const double m_m, m_omega;

and "initpot.cpp"

//Class for initial potentials

#include <cmath>
#include "initpot.h"

InitialPot::~InitialPot() {}HOPot::HOPot(double m, double omega): m_m(m), m_omega(omega) {}

double HOPot::potential(double x) const
{
double hopot=(1/2)*m_m*(m_omega*m_omega)*(x*x);

return hopot;
}

InitialPot* HOPot::clone() const
{
return new HOPot(*this); //return a derived class object through a base class pointer
}

Density class "density.h"

#pragma once

class Densities
{
public:
virtual double density() const = 0;
virtual ~Densities();

virtual Densities* clone() const = 0;

private:

};class Theoreticaldensity: public Densities
{
public:
Theoreticaldensity(double eigenfunc, double x);
double density() const override;

Densities* clone() const override;

private:
Theoreticaldensity();
const double m_x;
double m_eigenfunc;

and "density.cpp"

# include "initpot.h"
# include "parameters.h"
# include "schroddy.h"
# include "density.h"
# include <cmath>/*============================================
* Theoretical density
*==========================================*/
Densities::~Densities() {}

Theoreticaldensity::Theoreticaldensity(double eigenfunc, double x): m_eigenfunc(eigenfunc), m_x(x) {}

double Theoreticaldensity::density() const
{
//m_eigenfunc& schroddy::normalPsi;
return (1/(4*Parameters::PI*(m_x*m_x)))*Parameters::NN*(m_eigenfunc*m_eigenfunc);
}

Densities* Theoreticaldensity::clone() const
{
return new Theoreticaldensity(*this);
}

I hope I posted in the right session, I have not found a specific ...
Could someone help me? thanks!

 

[KataKoniK]

Hi there,

I am a fellow scientist and I would be happy to help you with your C++ program for solving the Schrodinger equation and calculating nuclear density.

To address your first question, yes, in order to solve the density equation, you will need to have a series of eigenfunctions for each value of n and l. In your code, you have already defined a class for generating eigenvalues using the shooting method, but you will also need to implement a class for generating eigenfunctions. This can be done by using the eigenvalues generated by the shooting method and solving the Schrodinger equation for each eigenvalue.

For your second question, I would suggest breaking down the shooting method into smaller steps and implementing them one by one. For example, you can start by defining the test eigenvalues and calculating the relative eigenfunctions. Then, you can work on checking the agreement of the eigenfunctions and replacing the eigenvalue between the two test eigenvalues. Finally, you can add the bisection loop and the condition for stopping the loop. If you are still having trouble implementing the shooting method, I would suggest reaching out to a colleague or a professor for assistance.

I hope this helps and good luck with your program!