//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*= // Program radschro //=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*= #include "cip.h" #include "nxgraph.h" // If you want to use gnuplot to plot the radial wave function // uncomment below line and re-direct the output of this program // into a proper file. // //#define USE_GNUPLOT //---- definition of fundermental constants //---- physical settings const double Z = 1.0; //---- experimental settings const double DX = 0.0625; // size of a cell in x-space const int R_NUM = 256; // total number of grids const int R_FAR = R_NUM-1; // index of farthest grid const int R_O = 160; // index of Bohr radius //---- graphics settings const int WIN_WIDTH = 512; const int WIN_HEIGHT = 256; const int BASE_Y = 160; const double MAG_Y = 196.0; const double OUTPUT_R_MAX = 64.0; // max radius output [a.u] //---- table of Radius class RTable { double r[R_NUM]; public: RTable( void ){ for( int i=0 ; i<=R_FAR ; i++ ){ r[i] = (1.0/Z)*exp(DX*(i-R_O)); } } inline double operator [] ( int i ){ return r[i]; } } R; //---- declaration of class Potential class Potential { double U[R_NUM]; // rV, potential muliplied by radius public: Potential( void ); double operator [] ( int i ){ return U[i]; } }; //---- sets potential U=rV as hydrogenic potential Potential::Potential( void ) { for( int i=0 ; i<=R_FAR ; i++ ){ U[i] = -Z; // Coulomb potential of nuclear } } //---- declaration of class Orbital class Orbital { int n, l; // quantum numbers double E; // energy of this orbital double G[R_NUM], F[R_NUM]; // wave function public: void Assign( int _n, int _l ); void Adjust( void ); void Draw( void ); void Output( void ); private: inline double CalcDG( int i, double P[], double G, double F ); inline double CalcDF( int i, double P[], double G, double F ); int PreInteg( double P[] ); int HammingA( double P[], int i_fit ); int HammingB( double P[], int i_fit ); void Normalize( double a, int i_fit ); }; //---- definition of class Orbital //-- sets quantum numbers void Orbital::Assign( int _n, int _l ) { n = _n; l = _l; E = -0.5*Z*Z/(n*n)-0.01; // a guess of the true energy } //-- adjusts energy so that an eigen wave function can stand void Orbital::Adjust( void ) { static double P[R_NUM]; double GA, GB, LogA, LogB, dE; int true_node = n-l-1, node, i_fit; do{ i_fit = PreInteg( P ); // prepares integration node = HammingA( P, i_fit ); // integrates from origin GA = G[i_fit]; // stores G at joint grid LogA = F[i_fit]/(R[i_fit]*G[i_fit]); // stores logarithmic derivative of G at joint grid node += HammingB( P, i_fit ); // integrates from far GB = G[i_fit]; // stores G at joint grid LogB = F[i_fit]/(R[i_fit]*G[i_fit]); // stores logarithmic derivative of G at joint grid Normalize( GA/GB, i_fit ); // connects and normalize dE = -0.5*sqr(G[i_fit])*(LogB-LogA); // estimates energy modification if( node != true_node ){ // if node is violated dE = double(true_node-node)/(n*n*n); // rough modification }else if( fabs(dE) > 1.0e-2 ){ // if still large error dE *= 0.5; // smaller modification } E+=dE; // modifies energy if( E>0.0 ){ // it can't be binding state fprintf( stderr, "Energy can not converge.\n" ); exit(1); } }while( fabs(dE)>1e-8 ); // repeats until the modification is sufficiently small fprintf( stderr, "Orbital energy of (%d %d) has converged to %18.16f [a.u]\n", n, l, E ); } //---- prepares for integration int Orbital::PreInteg( double P[] ) { extern Potential U; int i; // sets parameters for differential equation for( i=0 ; i<=R_FAR ; i++ ){ P[i] = 2*R[i]*U[i] - 2*sqr(R[i])*E + l*(l+1); } // sets G,F at near origin for( i=0 ; i<4 ; i++ ){ G[i] = pow(R[i],l+1); F[i] = (l+1)*G[i]; } // sets G,F at far away double a = -sqrt(-2*E); for( i=R_FAR ; i>R_FAR-4 ; i-- ){ G[i] = exp(a*R[i]); F[i] = (a*R[i])*G[i]; } // searches the cross grid and reports it as the joint grid for( i=R_FAR ; P[i]>0.0 && i>=0 ; i-- ); return i; } //---- represents derivative of G inline double Orbital::CalcDG( int i, double P[], double G, double F ){ return F; } //---- represents derivative of F inline double Orbital::CalcDF( int i, double P[], double G, double F ){ return P[i]*G + F; } //---- integrates from origin to the fitting grid int Orbital::HammingA( double P[], int i_fit ) { int i, node=0; double Gp, Fp, Gm, Fm, Gc, Fc, Gpc=0.0, Fpc=0.0, dGm, dFm; static double dG[R_NUM], dF[R_NUM]; // prepares derivatives at grids near origin for( i=0 ; i<4 ; i++ ){ dG[i] = CalcDG( i, P, G[i], F[i] ); dF[i] = CalcDF( i, P, G[i], F[i] ); } // integrates from the origin to the fitting grid for( ; i<=i_fit ; i++ ){ // calculates predictors Gp = G[i-4] + (4.0*DX/3.0)*( 2*dG[i-1] - dG[i-2] + 2*dG[i-3] ); Fp = F[i-4] + (4.0*DX/3.0)*( 2*dF[i-1] - dF[i-2] + 2*dF[i-3] ); // calculates modifiers Gm = Gp - (112.0/121.0)*( Gpc ); Fm = Fp - (112.0/121.0)*( Fpc ); // calculates derivative at modifiers dGm = CalcDG( i, P, Gm, Fm ); dFm = CalcDF( i, P, Gm, Fm ); // calculates correctors Gc = (9.0/8.0)*(G[i-1]) - (1.0/8.0)*(G[i-3]) + (3.0*DX/8.0)*(dGm + 2*dG[i-1] - dG[i-2]); Fc = (9.0/8.0)*(F[i-1]) - (1.0/8.0)*(F[i-3]) + (3.0*DX/8.0)*(dFm + 2*dF[i-1] - dF[i-2]); // calculates difference between predictors and correctors Gpc = Gp - Gc; Fpc = Fp - Fc; // calculates the next step G[i] = Gc + (9.0/121.0)*Gpc; F[i] = Fc + (9.0/121.0)*Fpc; // calculates derivative at the next step dG[i] = CalcDG( i, P, G[i], F[i] ); dF[i] = CalcDF( i, P, G[i], F[i] ); if( G[i]*G[i-1] < 0 ) node++; } return node; } //---- integrates from far to the fitting grid int Orbital::HammingB( double P[], int i_fit ) { int i, node=0; double Gp, Fp, Gm, Fm, Gc, Fc, Gpc=0.0, Fpc=0.0, dGm, dFm; static double dG[R_NUM], dF[R_NUM]; // prepares derivatives at far grids for( i=R_FAR ; i>R_FAR-4 ; i-- ){ dG[i] = CalcDG( i, P, G[i], F[i] ); dF[i] = CalcDF( i, P, G[i], F[i] ); } // integrates from far to the fitting grid for( ; i>=i_fit ; i-- ){ // calculates predictors Gp = G[i+4] - (4.0*DX/3.0)*( 2*dG[i+1] - dG[i+2] + 2*dG[i+3] ); Fp = F[i+4] - (4.0*DX/3.0)*( 2*dF[i+1] - dF[i+2] + 2*dF[i+3] ); // calculates modifiers Gm = Gp - (112.0/121.0)*( Gpc ); Fm = Fp - (112.0/121.0)*( Fpc ); // calculates derivative at modifiers dGm = CalcDG( i, P, Gm, Fm ); dFm = CalcDF( i, P, Gm, Fm ); // calculates correctors Gc = (9.0/8.0)*(G[i+1]) - (1.0/8.0)*(G[i+3]) - (3.0*DX/8.0)*(dGm + 2*dG[i+1] - dG[i+2]); Fc = (9.0/8.0)*(F[i+1]) - (1.0/8.0)*(F[i+3]) - (3.0*DX/8.0)*(dFm + 2*dF[i+1] - dF[i+2]); // calculates difference between predictors and correctors Gpc = Gp - Gc; Fpc = Fp - Fc; // calculates the next step G[i] = Gc + (9.0/121.0)*Gpc; F[i] = Fc + (9.0/121.0)*Fpc; // calculates derivative at the next step dG[i] = CalcDG( i, P, G[i], F[i] ); dF[i] = CalcDF( i, P, G[i], F[i] ); if( G[i]*G[i+1] < 0 ) node++; } return node; } //---- normalizes wave function void Orbital::Normalize( double a, int i_fit ) { int i; double sum, norm; for( i=i_fit ; i<=R_FAR ; i++ ){ G[i] *= a; F[i] *= a; } for( i=0, sum=0.0 ; i<=R_FAR ; i++ ){ sum += sqr(G[i])*R[i]*DX; } for( i=0, norm=sqrt(1.0/sum) ; i<=R_FAR ; i++ ){ G[i] *= norm; F[i] *= norm; } } //---- draws wave function void Orbital::Draw( void ) { const double MAG_X = (double)WIN_WIDTH/OUTPUT_R_MAX; NXDrawMoveto(); for( int i=0 ; i<=R_FAR ; i++ ){ NXDrawLineto( int(R[i]*MAG_X), BASE_Y - int(G[i]*MAG_Y) ); } } //---- outputs wave function void Orbital::Output( void ) { for( int i=0 ; i<=R_FAR ; i++ ){ printf("%f %f\n", R[i], G[i] ); if( R[i]>OUTPUT_R_MAX ) break; } printf("\n"); } //---- definition of fundermental static Potential U; static Orbital O; //---- main function int main( void ) { NXOpenWindow( "Radial Wave Function of Hydrogen Atom", WIN_WIDTH, WIN_HEIGHT ); NXClearWindow(); for( int n=1 ; n<=3 ; n++ ){ for( int l=0 ; l