よーめー日記: 8月 2013

電子の軌道の計算方法を勉強していたら、原子間の共有結合の結合角度と結合距離さえわかればどんな分子でも構造を表現できることがわかりました。

さらに調べると、世の中にはProtain Data Bank(http://www.rcsb.org/pdb/home/home.do)というところにたんぱく質と遺伝子の情報が全部数値化されて公開されているようです。

表示は、PyMol(http://www.pymol.org/)という分子モデルビュワーがあってそれを使うときれいなCGが作れるようです。

たんぱく質は数十から数百のアミノ酸のできていて、そのアミノ酸同士の結合パターンは10種類。
n個のアミノ酸からできていると10のn乗回の計算をしないと、それぞれの部分がどの結合になるのかコンピュータで構造を求められないのもわかりました。だから解析装置で解析しないといけないのね。

ということで、僕の長年の夢である、ヒトのヘモグロビンのデータ(http://www.rcsb.org/pdb/explore.do?structureId=1LFQ)を取得して、たんぱく質の立体構造を出力してみました。

黄色の色のアンテナの網みたいな部分に酸素がくっつくみたい。

いろいろなウイルスとかのデータもあって結構面白い。
僕が大学生のときにこんな面白いものがあったらもっと勉強したのに。

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

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;

}

よーめー日記

2013年8月31日土曜日

透き通るような青

2013年8月28日水曜日

たんぱく質の構造

2013年8月27日火曜日

電子の軌道計算