|
|
【答案】应助回帖
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
lijunlhc: 金币+20, ★★★很有帮助 2014-07-01 08:57:25
只有类氢原子才有1s 2s 2p这样的解。实际体系中由于势场不再是中心势，角动量不再是守恒量，不再存在s,p这样的态。
这是一个输出类氢原子波函数的C程序源码：
/******************************************************************************
This program generates wavefunction for hydrogen-like atoms using analytical
solution to Shrodinger's equation.
This program requires no command-line parameters, but an input file called
orb.inp. This input file must be in the form:
---------------------------------------------
H4f2.real.cube // output filename
1 // which part to output, 1 for real, -1 for imaginary
1 // Zval of the nuclea
4 3 2 // n l m
20.0 20.0 20.0 // size of the box, which will be multiplied by 2
0.20 0.20 0.20 // resolution
----------------------------------------------
******************************************************************************/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
double factor( int );
double ReNormFactor( int, int, int );
double KummerFunc( int, int, double );
double RadialFunc( int, int, int, double );
double RealSphereHarm( int, int, double, double );
double ImagSphereHarm( int, int, double, double );
int checkinp( int, int, int );
double factor( int n )
{
if ( n == 0 )
return 1.0;
else{
double y;
int i;
y = 1.0;
for( i = 1; i < n + 1; i ++ )
y = y * (double)(i);
return y;
}
}
double ReNormFactor( int Z, int n, int l )
{
double A, B, C;
A = 2.0 / ( (double)(n*n) * factor(2*l+1));
B = sqrt( factor(n+l) / factor(n-l-1) );
C = sqrt( (double)(Z*Z*Z) );
return A * B * C;
}
double KummerFunc( int alpha, int gamma, double ita )
{
double sum, ck;
int k;
sum = 1.0;
ck = 1.0;
for( k = 1; k < -alpha+1; k ++ )
{
ck = ck * (double)(alpha+k-1) / (double)(gamma+k-1) * (1.0/(double)(k));
sum = sum+ ck * pow( ita, (double)(k) );
}
return sum;
}
double RadialFunc( int Z, int n, int l, double r )
{
double ita, NormFactor, A, B, C;
ita = 2.0 * (double)(Z) / (double)(n) * r;
NormFactor = ReNormFactor( Z, n, l );
A = exp( -0.5 * ita );
B = pow( ita, (double)(l) );
C = KummerFunc( -n+l+1, 2*l+2, ita );
return NormFactor * A * B * C;
}
double RealSphereHarm( int l, int m, double theta, double phi )
{
if (l == 0)
return 2.820947917738782e-1;
else if (l == 1)
switch(m)
{
case -1: return 3.454941494713355e-1 * cos(m*phi) * sin(theta); break;
case 0: return 4.886025119029199e-1 * cos(m*phi) * cos(theta); break;
case 1: return -3.454941494713355e-1 * cos(m*phi) * sin(theta); break;
default: return 0.0; break;
}
else if (l == 2)
switch(m)
{
case -2: return 3.862742020231896e-1 * cos(m*phi) * sin(theta) * sin(theta); break;
case -1: return 7.725484040463792e-1 * cos(m*phi) * sin(theta) * cos(theta); break;
case 0: return 3.153915652525200e-1 * cos(m*phi) * ( 3.0*cos(theta)*cos(theta) - 1.0 ); break;
case 1: return -7.725484040463792e-1 * cos(m*phi) * sin(theta) * cos(theta); break;
case 2: return 3.862742020231896e-1 * cos(m*phi) * sin(theta) * sin(theta); break;
default: return 0.0; break;
}
else if (l == 3)
switch(m)
{
case -3: return 4.172238236327841e-1 * cos(m*phi) * sin(theta) * sin(theta) * sin(theta); break;
case -2: return 1.021985476433282 * cos(m*phi) * sin(theta) * sin(theta) * cos(theta); break;
case -1: return 3.231801841141507e-1 * cos(m*phi) * sin(theta) * ( 5.0*cos(theta)*cos(theta) - 1.0 ); break;
case 0: return 3.731763325901154e-1 * cos(m*phi) * cos(theta) * ( 5.0*cos(theta)*cos(theta) - 3.0 ); break;
case 1: return -3.231801841141507e-1 * cos(m*phi) * sin(theta) * ( 5.0*cos(theta)*cos(theta) - 1.0 ); break;
case 2: return 1.021985476433282 * cos(m*phi) * sin(theta) * sin(theta) * cos(theta); break;
case 3: return -4.172238236327841e-1 * cos(m*phi) * sin(theta) * sin(theta) * sin(theta); break;
default: return 0.0; break;
}
else
return 0.0;
}
double ImagSphereHarm( int l, int m, double theta, double phi )
{
if (l == 0)
return 0.0;
else if (l == 1)
switch(m)
{
case -1: return 3.454941494713355e-1 * sin(m*phi) * sin(theta); break;
case 0: return 4.886025119029199e-1 * sin(m*phi) * cos(theta); break;
case 1: return -3.454941494713355e-1 * sin(m*phi) * sin(theta); break;
default: return 0.0; break;
}
else if (l == 2)
switch(m)
{
case -2: return 3.862742020231896e-1 * sin(m*phi) * sin(theta) * sin(theta); break;
case -1: return 7.725484040463792e-1 * sin(m*phi) * sin(theta) * cos(theta); break;
case 0: return 3.153915652525200e-1 * sin(m*phi) * ( 3.0*cos(theta)*cos(theta) - 1.0 ); break;
case 1: return -7.725484040463792e-1 * sin(m*phi) * sin(theta) * cos(theta); break;
case 2: return 3.862742020231896e-1 * sin(m*phi) * sin(theta) * sin(theta); break;
default: return 0.0; break;
}
else if (l == 3)
switch(m)
{
case -3: return 4.172238236327841e-1 * sin(m*phi) * sin(theta) * sin(theta) * sin(theta); break;
case -2: return 1.021985476433282 * sin(m*phi) * sin(theta) * sin(theta) * cos(theta); break;
case -1: return 3.231801841141507e-1 * sin(m*phi) * sin(theta) * ( 5.0*cos(theta)*cos(theta) - 1.0 ); break;
case 0: return 3.731763325901154e-1 * sin(m*phi) * cos(theta) * ( 5.0*cos(theta)*cos(theta) - 3.0 ); break;
case 1: return -3.231801841141507e-1 * sin(m*phi) * sin(theta) * ( 5.0*cos(theta)*cos(theta) - 1.0 ); break;
case 2: return 1.021985476433282 * sin(m*phi) * sin(theta) * sin(theta) * cos(theta); break;
case 3: return -4.172238236327841e-1 * sin(m*phi) * sin(theta) * sin(theta) * sin(theta); break;
default: return 0.0; break;
}
else
return 0.0;
}
int checkinp( int l, int m, int flavor )
{
if ( l < 0 )
return -1;
else if ( l > 3 )
return -2;
else if ( m < 0 && -m > l )
return -3;
else if ( m >= 0 && m > l )
return -3;
else if ( flavor != -1 && flavor != 1 )
return -4;
else
return 0;
}
int main( void )
{
// these variables are related to orb.inp
char cubename[32];
int flavor;
int Z;
int n, l, m;
double xmax, ymax, zmax;
double hx, hy, hz;
// these variables are related to mesh-info
int NGX, NGY, NGZ, NX, NY, NZ;
int i, j, k;
double orgx, orgy, orgz;
double x, y, z, r, theta, phi, gridval;
// files used through the program
FILE *cube, *inp;
// program status
int status;
// global constant
const double PI = 3.14159265358979;
// read in 'orb.inp'
inp = fopen("orb.inp","r" ;
fscanf( inp, "%s\n", &cubename );
fscanf( inp, "%d\n", &flavor );
fscanf( inp, "%d\n", &Z );
fscanf( inp, "%d%d%d\n", &n, &l, &m );
fscanf( inp, "%lf%lf%lf\n", &xmax, &ymax, &zmax );
fscanf( inp, "%lf%lf%lf\n", &hx, &hy, &hz );
fclose( inp );
// check input parameters
status = checkinp( l, m, flavor );
if ( status != 0 )
{
printf( "\nProgram exited with error code:" );
printf( "%4d\n", status );
printf( "Check your input" );
printf("\n";
return status;
}
else
{
// calculate NGX, NGY, NGZ
NGX = (int)(xmax/hx);
NGY = (int)(ymax/hy);
NGZ = (int)(zmax/hz);
// origin of grid
orgx = -hx * (NGX - 1) + hx * 0.5;
orgy = -hy * (NGY - 1) + hy * 0.5;
orgz = -hz * (NGZ - 1) + hz * 0.5;
// size of grid
NX = 2 * NGX;
NY = 2 * NGY;
NZ = 2 * NGZ;
// open cube file
cube = fopen(cubename,"w";
// write first 2 lines to cube
if ( flavor == 1)
fprintf( cube, "Real part of atomic orbitals of H-like atoms\n";
else
fprintf( cube, "Imag part of atomic orbitals of H-like atoms\n";
fprintf( cube, "Generated by programm orb.x\n";
// write natom and origin of data
fprintf( cube, "%5d%12.6f%12.6f%12.6f\n", 1, orgx, orgy, orgz );
// write mesh info
fprintf( cube, "%5d%12.6f%12.6f%12.6f\n", NX, hx, 0.0, 0.0 );
fprintf( cube, "%5d%12.6f%12.6f%12.6f\n", NY, 0.0, hy, 0.0 );
fprintf( cube, "%5d%12.6f%12.6f%12.6f\n", NZ, 0.0, 0.0, hz );
// write atoms info
fprintf( cube, "%5d%12.6f%12.6f%12.6f%12.6f\n", Z, 1.0*Z, 0.0, 0.0, 0.0 );
// entering main loop
if ( flavor == 1 )
{
for ( i = 1; i <= NX; i ++ )
for ( j = 1; j <= NY; j ++ )
for( k = 1; k <= NZ; k ++ )
{
x = orgx + hx * ( i - 1 );
y = orgy + hy * ( j - 1 );
z = orgz + hz * ( k - 1 );
r = sqrt( x*x + y*y + z*z );
theta = acos( z/r );
if ( y >= 0.0 )
phi = acos( x / sqrt(x*x+y*y) );
else
phi = 2.0 * PI - acos( x / sqrt(x*x+y*y) );
gridval = RadialFunc( Z, n, l, r ) * RealSphereHarm( l, m, theta, phi );
fprintf( cube, "%13.5e", gridval );
if ( k == NZ || k%6==0 )
fprintf( cube, "\n" );
}
}
else
{
for ( i = 1; i <= NX; i ++ )
for ( j = 1; j <= NY; j ++ )
for( k = 1; k <= NZ; k ++ )
{
x = orgx + hx * ( i - 1 );
y = orgy + hy * ( j - 1 );
z = orgz + hz * ( k - 1 );
r = sqrt( x*x + y*y + z*z );
theta = acos( z/r );
if ( y >= 0.0 )
phi = acos( x / sqrt(x*x+y*y) );
else
phi = 2.0 * PI - acos( x / sqrt(x*x+y*y) );
gridval = RadialFunc( Z, n, l, r ) * ImagSphereHarm( l, m, theta, phi );
fprintf( cube, "%13.5e", gridval );
if ( k == NZ || k%6==0 )
fprintf( cube, "\n" );
}
}
fclose(cube);
return 0;
}
}
=================================================== |
|
回复此楼