2013年8月27日火曜日

電子の軌道計算

最近、大学時代に勉強したけどわからなかったことをインターネットで勉強してます。
まずこれ、シュレーディンガー方程式。

http://ja.wikipedia.org/wiki/%E3%82%B7%E3%83%A5%E3%83%A
C%E3%83%BC%E3%83%87%E3%82%A3%E3%83%B3%E3%82%AC%E3%83%BC%E6%96%B9%E7%A8%8B%E5%BC%8F
 
大学時代に、ハミルトニアンとか難解な数式がたくさんあってさっぱり意味がわからなかんですが、最近はネット上にいろいろな資料が乗っていて、一週間くらい毎日朝から深夜までひたすら物理を勉強してたら、ついにわかった。


方程式の解を求めるのは至難の業なのですが、誰かが説いた解があれば、それをもとにいろいろな計算ができるようです。

水素原子の周りをまわる電子の軌道の波動関数を計算してみました。



 

わかってしまえば意外と簡単じゃん。
こんなふうに簡単に電子の軌道を計算して図にすることができます。
次はタンパク質とか水分子の軌道計算してみようかなぁ。


-----------------------------------

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

double PI = 3.141592653589793238462;
double Factorial(int n);
double Factorial2(int n, int m);
const int n_max = 4;
const int kizami = 200;
const double L = 80.0;
double Laguerre( int k, int n, double x);           //ラゲール陪関数
double Legendre( int m, int l, double x);           //ルジャンドル陪関数

double bai=20000;


void* cbmp(int w,int h){
unsigned char* ret;
int sz;
sz=w*h*4+54;
ret=malloc(sz);
memset(ret,0,sz);

ret[0]='B';
ret[1]='M';

ret[2]=sz &255;
ret[3]=(sz>>8) &255;
ret[4]=(sz>>16) &255;
ret[5]=(sz>>24) &255;

ret[10]=54;

ret[14]=40;

ret[18]=w &255;
ret[19]=(w>>8) &255;

ret[22]=(-h) &255;
ret[23]=((-h)>>8) &255;
ret[24]=255;
ret[25]=255;

ret[26]=1;
ret[28]=32;



return (void*)ret;

}

void fbmp(void* vp)
{
if(vp==NULL)return;
free(vp);
}

void wbmp(void* vp,char* fn)
{
int sz;
unsigned char* p;
FILE* fp;

if(vp==NULL || fn==NULL)return;

p=vp;
p+=2;

sz=(*p);
sz|=(*(p+1))<<8;
sz|=(*(p+2))<<16;
sz|=(*(p+3))<<24;

fp=fopen(fn,"wb");
if(fp==NULL)return;
fwrite(vp,1,sz,fp);
fclose(fp);

}

void pbmp(void*vp, int x,int y,int b,int g,int r)
{
int off;
unsigned char* p;
int w;


p=vp;

w=p[18] | p[19]<<8;

off=y*w*4+x*4+54;
if(r<0)r=0;
if(r>255)r=255;
if(g<0)g=0;
if(g>255)g=255;
if(b<0)b=0;
if(b>255)b=255;

p[off]=b;
p[off+1]=g;
p[off+2]=r;
p[off+3]=255;

}

void abmp(void*vp, int x,int y,int b,int g,int r)
{
int off;
unsigned char* p;
int w;


p=vp;

w=p[18] | p[19]<<8;

off=y*w*4+x*4+54;
if(r<0)r=0;
if(r>255)r=255;
if(g<0)g=0;
if(g>255)g=255;
if(b<0)b=0;
if(b>255)b=255;

b+=p[off];
g+=p[off+1];
r+=p[off+2];

if(r<0)r=0;
if(r>255)r=255;
if(g<0)g=0;
if(g>255)g=255;
if(b<0)b=0;
if(b>255)b=255;

p[off]=b;
p[off+1]=g;
p[off+2]=r;

}


int main(){

  double r_min = 0.0;
  double r_max =  40.0;
  const int Z = 1;

  int l,m,n;




  for(n=1; n<=n_max; n++)
  {
    for(l=0; l<n; l++)
    {
      for(m=0; m<=l; m++)
      {
        void* fp;
        char str[200];
int ix,iy;

sprintf(str, "Psi%d_%d_%d.bmp", n, l, m);
fp=cbmp(kizami,kizami);

 for(ix=0; ix<kizami; ix++ )
        {
          for(iy=0; iy<kizami; iy++ )
          {
            double x = - L/2.0 + L * (double)(ix)/(double)(kizami);
            double y = - L/2.0 + L * (double)(iy)/(double)(kizami);
            double z =0;
            double r = sqrt(pow(x,2)+pow(y,2)+pow(z,2));
            double theta = acos(z/r);
            double phi;
double R,Y;
int cc;


            if(y>=0) phi = acos(x/r);
            else phi = 2.0*PI - acos(x/r);
            if(r==0){
              theta = 0;
              phi = 0;
            }
            R = 2.0/pow((double)(n),2) * sqrt(Factorial(n-l-1)/Factorial(n+l))
* exp(-r/(double)(n)) * pow( 2.0*r/(double)(n),l) 
*Laguerre(2*l+1, n-l-1, 2.0*r/(double)(n))  ;
            Y = pow(-1.0, m) * sqrt((double)(l)+1.0/2.0) 
* Factorial(l-m)/Factorial(l+m) * Legendre(m,l,cos(theta));
            if(m==0) {
              //fprintf(fp,"%f %f %f\n",x,y,R*Y);
cc=(int)(R*Y*bai);
if(cc>=0)
pbmp(fp,ix,iy,0,0,cc);
else
pbmp(fp,ix,iy,-cc,0,0);
            }else{
              double Yx = Y * cos((double)(m)*phi);
              double Yy = Y * sin((double)(m)*phi);
              //fprintf(fp,"%f %f %f %f\n",x,y,R * Yx,R * Yy);
cc=(int)(R*Yx*bai);
if(cc>=0)
pbmp(fp,ix,iy,0,0,cc);
else
pbmp(fp,ix,iy,-cc,0,0);
cc=(int)(R*Yx*bai);
if(cc>=0)
abmp(fp,ix,iy,0,0,cc);
else
abmp(fp,ix,iy,-cc,0,0);
            }
          }
        }
wbmp(fp,str);
        fbmp(fp);
      }
    }
  }
}
double Laguerre( int k, int n, double x)
{
  double sum=0;
  int m;
  for(m=0; m<=n; m++){
    sum += pow(-x,m)*Factorial(n+k)/Factorial(n-m)/Factorial(k+m)/Factorial(m);
  }
  return sum;
}
double Legendre( int m, int l, double x)
{
  int mm = abs(m);
  int ll;
  double r0, r1, r2;

  if( mm>l )  return 0;
    r0 = 0.0;
    r1 = pow(1.0-x*x, (double)(mm)/2.0) * Factorial2(2*mm-1, 2);
    if( mm==l && m>=0) return r1;
    if( mm==l && m<0)  return r1 * pow(-1.0,mm) * Factorial(l-mm)/Factorial(l+mm) ;
    for(ll = mm+1; ll<=l; ll++ ){
      r2 = ((2.0*ll-1.0)*x*r1 - (ll+mm-1.0)*r0)/(ll-mm);
      r0 = r1;
      r1 = r2;
    }
  if(m>=0) return r2;
  else return r2 * pow(-1.0,mm) * Factorial(l-mm)/Factorial(l+mm) ;
}
double Factorial(int n)
{
  double F = 1.0;
  int i;

  if( n<=0 ) return 1.0;
  for(i=n; i>=2; i--){
    F *= (double)(i);    
  }
  return F;
}
double Factorial2(int n, int m)
{
  double F = 1.0;
  int i;
  if( n<=0 ) return 1.0;
  for(i=n; i>=2; i=i-m){
    F *= (double)(i);    
  }
  return F;
}


0 件のコメント:

コメントを投稿