| 6 | 1/1 | 返回列表 |
| 查看: 1924 | 回復(fù): 5 | ||
那年的小艾金蟲(chóng) (初入文壇)
|
[求助]
F軌道電子云分布 已有1人參與
|
| F軌道的電子云分布應(yīng)該怎么畫(huà),用一些簡(jiǎn)單的軟件可以做么?我想畫(huà)出非共線情況下的F軌道的電子云分布,有沒(méi)有大俠會(huì)的,望指教一二。還有想問(wèn)下網(wǎng)絡(luò)上的一些軌道電子云彩圖中對(duì)稱結(jié)構(gòu)中用的兩種顏色是不是區(qū)別自旋向上向下的? |
木蟲(chóng) (著名寫(xiě)手)
|
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ù)值寫(xiě)入cube文件,然后就可以用vesta之類的軟件打開(kāi)看了。 網(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ù)。不是自旋。 |
木蟲(chóng) (著名寫(xiě)手)
|
這是之前寫(xiě)的一個(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; } } |
木蟲(chóng) (著名寫(xiě)手)
金蟲(chóng) (初入文壇)
| 6 | 1/1 | 返回列表 |
| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 300求調(diào)劑,材料科學(xué)英一數(shù)二 +9 | leaflight 2026-03-24 | 9/450 |
|
|---|---|---|---|---|
|
[考研] 327求調(diào)劑 +6 | 汲亦昊 2026-03-29 | 6/300 |
|
|
[考研] 343求調(diào)劑 +6 | 愛(ài)羈絆 2026-03-29 | 6/300 |
|
|
[考研] 2026年華南師范大學(xué)歡迎化學(xué),化工,生物,生醫(yī)工等專業(yè)優(yōu)秀學(xué)子加入! +3 | llss0711 2026-03-28 | 6/300 |
|
|
[考研] 數(shù)一英一271專碩(085401)求調(diào)劑,可跨 +7 | 前行必有光 2026-03-28 | 8/400 |
|
|
[考研] 316求調(diào)劑 +7 | 江辭666 2026-03-26 | 7/350 |
|
|
[考研] 求調(diào)劑 +7 | 爭(zhēng)取九點(diǎn)睡 2026-03-28 | 8/400 |
|
|
[考研]
|
y7czhao 2026-03-26 | 10/500 |
|
|
[考研] 286求調(diào)劑 +4 | 丟掉懶惰 2026-03-27 | 7/350 |
|
|
[考研] 275求調(diào)劑 +10 | jjjjjjjjjjl 2026-03-27 | 10/500 |
|
|
[考博] 26申博 +3 | 加油沖。 2026-03-26 | 3/150 |
|
|
[考研] 08開(kāi)頭275求調(diào)劑 +4 | 拉誰(shuí)不重要 2026-03-26 | 4/200 |
|
|
[考研] 求調(diào)劑,一志愿 南京航空航天大學(xué)大學(xué) ,080500材料科學(xué)與工程學(xué)碩 +4 | @taotao 2026-03-26 | 5/250 |
|
|
[考研] 336材料求調(diào)劑 +7 | 陳瀅瑩 2026-03-26 | 9/450 |
|
|
[考研] 341求調(diào)劑 +7 | 青檸檬1 2026-03-26 | 7/350 |
|
|
[考研] 281求調(diào)劑 +3 | 亞克西good 2026-03-26 | 5/250 |
|
|
[考研] 機(jī)械學(xué)碩310分,數(shù)一英一,一志愿211本科雙非找調(diào)劑信息 +3 | @357 2026-03-25 | 3/150 |
|
|
[考研] 309求調(diào)劑 +4 | gajsj 2026-03-25 | 5/250 |
|
|
[考研] 347求調(diào)劑 +4 | L when 2026-03-25 | 4/200 |
|
|
[考研] 材料/農(nóng)業(yè)專業(yè),07/08開(kāi)頭均可,過(guò)線就行 +3 | 呵唔哦豁 2026-03-23 | 4/200 |
|