| 6 | 1/1 | 返回列表 |
| 查看: 1920 | 回復(fù): 5 | ||
那年的小艾金蟲 (初入文壇)
|
[求助]
F軌道電子云分布 已有1人參與
|
| F軌道的電子云分布應(yīng)該怎么畫,用一些簡(jiǎn)單的軟件可以做么?我想畫出非共線情況下的F軌道的電子云分布,有沒(méi)有大俠會(huì)的,望指教一二。還有想問(wèn)下網(wǎng)絡(luò)上的一些軌道電子云彩圖中對(duì)稱結(jié)構(gòu)中用的兩種顏色是不是區(qū)別自旋向上向下的? |
木蟲 (著名寫手)
|
F軌道實(shí)際上是角動(dòng)量模方算符L^2對(duì)應(yīng)的l=3的本征態(tài)波函數(shù)。在中心勢(shì)場(chǎng)(比如氫原子)中,L^2與哈密頓量對(duì)易,而在固體或分子中,L^2與哈密頓量不一定對(duì)易。一般程序輸出的是哈密頓量與某個(gè)算符的本征態(tài)(對(duì)固體為平移算符),對(duì)這些本征態(tài)無(wú)法斷定l的取值。 如果是單原子,可以根據(jù)氫原子薛定諤方程解析解求出在一系列格點(diǎn)上的波函數(shù)值寫入cube文件,然后就可以用vesta之類的軟件打開看了。 網(wǎng)上的程序給出的電子云圖有兩種:一種直接輸出|n,l,m>,另一種輸出|n,l,m>的線性組合。第一種波函數(shù)值是復(fù)數(shù),作圖時(shí)給出的是實(shí)部或虛部;第二種波函數(shù)值是實(shí)數(shù),常說(shuō)的px,py,pz就是這種。這兩種情況下顏色一般用來(lái)區(qū)分波函數(shù)值的正負(fù)。不是自旋。 |
木蟲 (著名寫手)
|
這是之前寫的一個(gè)輸出cube格式的程序源代碼: /****************************************************************************** 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; } } |
木蟲 (著名寫手)
金蟲 (初入文壇)
| 6 | 1/1 | 返回列表 |
| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 0856求調(diào)劑 +13 | zhn03 2026-03-25 | 14/700 |
|
|---|---|---|---|---|
|
[考研] 求調(diào)劑 +7 | 爭(zhēng)取九點(diǎn)睡 2026-03-28 | 8/400 |
|
|
[考研]
|
y7czhao 2026-03-26 | 10/500 |
|
|
[考研] 317求調(diào)劑 +6 | 十閑wx 2026-03-24 | 6/300 |
|
|
[考研] 295求調(diào)劑 +5 | 1428151015 2026-03-27 | 6/300 |
|
|
[考研] 265求調(diào)劑 +8 | 小木蟲085600 2026-03-27 | 8/400 |
|
|
[考研] 070300化學(xué)求調(diào)劑 +4 | 起個(gè)名咋這么難 2026-03-27 | 4/200 |
|
|
[考研] 266分求材料化工冶金礦業(yè)等專業(yè)的調(diào)劑 +4 | 哇呼哼呼哼 2026-03-26 | 4/200 |
|
|
[考研] 274求調(diào)劑 +17 | 顧九笙要謙虛 2026-03-24 | 23/1150 |
|
|
[考研] 化學(xué)308分求調(diào)劑 +8 | 你好明天你好 2026-03-23 | 9/450 |
|
|
[考研] 314求調(diào)劑 +3 | 溪云珂 2026-03-26 | 3/150 |
|
|
[考研] 317求調(diào)劑 +7 | 蛋黃咸肉粽 2026-03-26 | 7/350 |
|
|
[考研] 343求調(diào)劑 +4 | 贈(zèng)我一本書 2026-03-23 | 4/200 |
|
|
[考研] 化學(xué)工程085602 305分求調(diào)劑 +17 | RichLi_ 2026-03-25 | 17/850 |
|
|
[考研] 一志愿南航 335分 | 0856材料化工 | GPA 4.07 | 有科研經(jīng)歷 +6 | cccchenso 2026-03-23 | 6/300 |
|
|
[考研] 284求調(diào)劑 +15 | Zhao anqi 2026-03-22 | 15/750 |
|
|
[考研] 277分求調(diào)劑,跨調(diào)材料 +3 | 考研調(diào)劑lxh 2026-03-24 | 3/150 |
|
|
[考研] 一志愿國(guó)科過(guò)程所081700,274求調(diào)劑 +3 | 三水研0水立方 2026-03-23 | 3/150 |
|
|
[考研] 一志愿重慶大學(xué)085700資源與環(huán)境,總分308求調(diào)劑 +7 | 墨墨漠 2026-03-23 | 8/400 |
|
|
[考研] 求老師收我 +3 | zzh16938784 2026-03-23 | 3/150 |
|