/****************************************************************************** Title : numerov_harmonic_oscillator.c Author : Stewart Weiss (Adaptation of algorithm in Levine, Quantum Chemistry to C ) Created on : September 10, 2008 Description : User guided search for solution to Hamiltonian Purpose : To solve harmonic oscillator equation numerically Build with : gcc -lm Usage : numerov [-hv] [-x leftvalue ] [-s segsize ] [ -N intvls ] [ -e quantum ] [output file]\n" -h : help\n" -x leftvalue : set leftmost x-coordinate to leftvalue\n" -s segsize : set segment size to segsize\n" -N intvls : set number of intervals to intvls\n" -e quantum : set energy level to quantum (0,1,2,...)\n" output file : file to write array of psi values\n" ******************************************************************************* Notes: This is based on the recurrence relation 2*Psi[n] - Psi[n-1] + 5*G[n]*Psi[n]*s*s/6 + G[n-1]*Psi[n-1]*s*s/12 Psi[n+1] = ------------------------------------------------------------------- 1 - G[n+1]*s*s/12 where G[n] is defined as G(x[n]) = number_of_intervals*(2*V(x[n]) - 2*E) Since V(x[n]) = (1/2)k*x[n]*x[n], G[n] is estimated as x[n]*x[n] - 2*E ******************************************************************************/ #include #include #include #include #include #define FALSE 0 #define TRUE 1 static const char *help = " [-h] [-x leftvalue] [-s segsize] [-N intvls] [-e quantum] [output file]\n" " -h : help\n" " -x leftvalue : set leftmost x-coordinate to leftvalue\n" " -s segsize : set segment size to segsize\n" " -N intvls : set number of intervals to intvls\n" " -e quantum : set energy level to quantum (0,1,2,...)\n" " output file : file to write array of psi values\n" "\nNotes\n" " By default, output will be written to the terminal.\n" " Any values not supplied on the command line will be prompted for by" " the program.\n" " Make sure that leftvalue is NEGATIVE and that " " 2*leftvalue + segsize*intvls = 0\n" " so that the arrays are symmetric around the origin." ; /*****************************************************************************/ /* Utility Function Prototypes */ /*****************************************************************************/ /*****************************************************************************/ /* Functions for Interface */ /*****************************************************************************/ /* * Parses the command line and fills the file_flags. It is given the * addresses of the argc and argv parameters so that it can update them * if needed. */ void get_options( int argc, char **argv, double *xval, double *segSize, int *numSegs, int *energylevel, FILE **fileptr ); // Waits until one of y,Y,n, or N is entered. If it returns without one of // these, it returns EOF int get_valid_char(); // Displays a help message void usage( char * prog); // Write the values of a[] to the file stream fp void write_values (FILE* fp, double a[], int n); /*****************************************************************************/ /* Functions for Numerov */ /*****************************************************************************/ /* * Engine of algorithm -- given energy guess, segment size, number of * intervals, as well as the array of x coordinates (x[]), and the specific * function that encapsulates the potential energy function (g[]), this * computes approximate values of psi. * This computation is dimensionless: the values of all variables have no * dimension. */ void compute_psi ( double energy_guess, double segment_size, int n_intervals, double x[], double g[], double psi[] ); // Stores into a_sqr[] the squares of values in a[] void compute_squares( double a[], double a_sqr[], int n); // Returns the sum of the elements in array a[] double array_sum ( double a[], int n ); // Modifies a[] so that a[i] = a[i]*normal void normalize_psi ( double a[], double a_sqr[], int n, double segment_size ); // Counts the number of nodes in the wave function int count_nodes ( double a[], int num_intervals ); /*****************************************************************************/ /* Main */ /*****************************************************************************/ int main( int argc, char* argv[]) { double *x; // points at which to estimate Psi double *G_r; // G = 2V_r - 2E_r at each point double *Psi_r; // Psi_r at each point double *Psi_r_sqr; // Psi_r^2 at each point double starting_position = 0; // leftmost point to evaluate psi double segment_size = 0; // size of segment dx for integration double energy_guess; // guessed value of energy eigenvalue int number_of_intervals = 0; // number of intervals for integration int number_of_nodes; // number of nodes in wave function int energy_level = -1; // quantum number for energy, 0,1,2... double sum_of_probabilities; // sum of psi_r_sqr values - should be 1 double normalizing_const; // to force sum of prob's to be 1.0 FILE *fp_out; // file stream to write psi_r array char reply; fp_out = stdout; // default output stream is terminal // First, retrieve command line options and file argument, if any get_options(argc, argv, &starting_position, &segment_size, &number_of_intervals, &energy_level, &fp_out); // Now check which parameters need to be obtained interactively if ( 0 > energy_level ) { printf("Enter the energy level: "); scanf( "%d", &energy_level ); } if ( 0 == starting_position ) { starting_position = -1 * ceil(sqrt( 2*energy_level)); printf("Enter the value of left point x_0 (recommend %f): ", starting_position); scanf ( "%lf", &starting_position ); } if ( 0 == segment_size ) { printf("Enter the value of segment size: "); scanf( "%lf", &segment_size ); } if ( 0 == number_of_intervals ) { number_of_intervals = (int) 2*ceil((-starting_position)/segment_size); printf("Enter the number of intervals (recommend %d): ", number_of_intervals ); scanf( "%d", &number_of_intervals ); } // Let user know parameters we're working with printf("x[0] = %f, segsize = %f, numintervals = %d, energy = %d\n", starting_position, segment_size, number_of_intervals, energy_level); // Allocate memory for arrays x = (double*) calloc(number_of_intervals+1, sizeof(double)); G_r = (double*) calloc(number_of_intervals+1, sizeof(double)); Psi_r = (double*) calloc(number_of_intervals+1, sizeof(double)); Psi_r_sqr = (double*) calloc(number_of_intervals+1, sizeof(double)); // If allocating any of them failed, then quit if ( 0 == x || 0 == G_r || 0 == Psi_r || 0 == Psi_r_sqr ) { perror("calloc()"); exit(1); } // Initialize arrays x[0] = starting_position; x[1] = x[0] + segment_size; Psi_r[0] = 0; Psi_r[1] = 0.0001 * (energy_level %2 ? -1 : 1); energy_guess = energy_level + 0.5; while ( TRUE ) { // compute psi array compute_psi(energy_guess, segment_size, number_of_intervals, x, G_r, Psi_r ); // normalize so that psi^2 is a probability density normalize_psi(Psi_r, Psi_r_sqr, number_of_intervals+1, segment_size); // count how many nodes are in the wave function number_of_nodes = count_nodes( Psi_r, number_of_intervals ); // Give user a chance to modify the energy guess printf("Energy Guess = %16.9f; Nodes = %2d; Psi_r[%d] = %16.9f\n", energy_guess, number_of_nodes, number_of_intervals, Psi_r[number_of_intervals]); printf("Try again?[y/n]"); fflush(stdout); reply = get_valid_char(); if ( (reply == 'n' ) || ( reply == 'N' ) ) break; printf("New energy guess: "); scanf( "%12lf", &energy_guess ); } // write the array of psi values to the output stream fprintf(fp_out, "Array of psi values:\n"); write_values (fp_out, Psi_r, number_of_intervals+1 ); // clean up free(x); free(G_r); free(Psi_r); free(Psi_r_sqr); return 0; } /******************************************************************************/ void compute_psi ( double energy_guess, double segment_size, int n_intervals, double x[], double g[], double psi[] ) { int i; double segsize_sqr = segment_size*segment_size/12; g[0] = x[0]*x[0] - 2*energy_guess; g[1] = x[1]*x[1] - 2*energy_guess; for ( i = 1; i <= n_intervals-1; i++ ) { x[i+1] = x[i] + segment_size; g[i+1] = x[i+1]*x[i+1] - 2*energy_guess; psi[i+1] = ( -psi[i-1] + 2*psi[i] + 10*g[i] * psi[i] * segsize_sqr + g[i-1] * psi[i-1] * segsize_sqr ) / ( 1.0 - g[i+1] * segsize_sqr); } } /******************************************************************************/ void normalize_psi ( double a[], double a_sqr[], int n, double segment_size ) { int k; double sum_of_probabilities; double normalizing_const; compute_squares(a, a_sqr, n); sum_of_probabilities = array_sum(a_sqr, n) * segment_size; normalizing_const = 1.0 / sqrt(sum_of_probabilities); for ( k = 0; k < n; k++) { a[k] = normalizing_const*a[k]; } } /******************************************************************************/ void compute_squares( double a[], double a_sqr[], int n) { int k; for ( k = 0; k < n; k++) { a_sqr[k] = a[k]*a[k]; } } double array_sum ( double a[], int n ) { int k; double sum = 0.0; for ( k = 0; k < n; k++) { sum += a[k]; } return sum; } /******************************************************************************/ /****************************************************************************** A node is a point at which the wave function crosses the x-axis between two extremes. During the approximation algorithm, it is possible for the wave function to cross the axis in the classically prohibited region as it approaches its zeros, simply because the estimate is bad. This function counts only true crossings by counting only those crossings such that the first derivate changes at some point after the crossing, i.e., delta-psi changes. ******************************************************************************/ int count_nodes ( double a[], int num_intervals ) { int k; int count = 0; int crossing_found = 0; // 1 when a crossing is found int direction_of_crossing = 0; // +1 if upward, -1 if downward // start at a[1] since a[0] == 0 for ( k = 1; k < num_intervals; k++ ) { if ( crossing_found ) { if ( (a[k+1] - a[k]) * direction_of_crossing < 0 ) { // a reversal of direction has occurred so the // previous crossing was a real node count++; crossing_found = 0; // reset flag direction_of_crossing = 0; // reset flag } } if ( a[k]*a[k+1] < 0 ) { // a new crossing was just found crossing_found = 1; direction_of_crossing = (a[k] < a[k+1])? 1 : -1 ; } } return count; } /******************************************************************************/ void get_options( int argc, char **argv, double *xval, double *segSize, int *numSegs, int *energylevel, FILE **fileptr ) { int ch; char options[] = ":hx:s:e:N:"; /* turn off error messages by getopt() */ opterr = 0; /* process the command line options */ while (TRUE) { ch = getopt(argc, argv, options); if ( -1 == ch ) break; switch ( ch ) { case 'h' : usage((argv)[0]); exit(0); case 'x': *xval = strtod(optarg,'\0'); break; case 's': *segSize = strtod(optarg,'\0'); break; case 'N': *numSegs = strtol(optarg,'\0',0); break; case 'e': *energylevel = strtol(optarg,'\0',0); break; case '?' : default: fprintf (stderr, "?? bad option 0%o ??\n", ch); break; } } if ( optind < argc ) if ( NULL == (*fileptr = fopen(argv[optind], "w")) ) { perror(argv[optind]); exit(1); } } /******************************************************************************/ int get_valid_char() { int n; unsigned char c; while ( ( c = fgetc(stdin) ) && strchr("yYnN",c) == NULL) /* no -op */ ; return c; } /******************************************************************************/ void usage( char * prog) { printf("Usage: %s", prog); printf("%s\n", help); } /******************************************************************************/ void write_values (FILE* fp, double a[], int n) { int i; for ( i = 0; i < n; i++ ) fprintf(fp, "%12.8f\n", a[i]); } /******************************************************************************/