サンプル C コンピュータープログラム(H2)

トップページ
分子結合のページに戻る。

下のソースプログラムをそのままテキストエディタ(メモ帳など)にコピー and ペースト すれば、簡単にコンパイルと実行できる。
(この プログラムは単純なC言語なので、このテキストエディタを "ファイル名.c" とセーブしてコンパイルしてほしい。)
ここでは 1 MM = 1 × 10-14 meter という新しい単位を使用している。

このサンプルプログラムでは、実行すると、最初に H2 の核間距離 ( MM ) を入力する。
その後、電子1の y座標 "r" ( MM )、x 座標 "a" ( MM ) を 入力する。
これら入力値から、このプログラムは H2 の結合エネルギー (eV)、原子核 1 に作用する トータルの力 (= nucforce )、電子1に作用する力の x, y 成分 (= elefx, elefy ) を出力する。
これらの力は 水素原子の基底状態の電子と原子核間の力との比で表されている。

電子に作用する力と その速度から 1軌道に含まれるドブロイ波を計算する。
( ここでは ビリアル定理 E = -T = 1/2 V を使用している。)
"another x" は 電子1が 力の方向へ移動した先の もう1つの x 座標である。
それから x が another x の座標のときの 上記の計算を行う。
"avebinding" は x 座標が "a" と "another x" のときの 平均の結合エネルギーであり、"avewave" は x 座標が "a" と "another x" のときの平均のドブロイ波である。
"after-elefx" は 力の方向に移動した後の 電子1に作用する力の x 成分である。
最初の x 座標は +100 まで自動的に計算ごとに増加していく。


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

void calc(double ca);
void main(void);
  
  double ke,r,nuc, rra,rrb,rrc;
  double poten,ppot, kinetic,velo,binding,elefx,elefy,eleforce;
  double radius, debroglie, wave, nucforce, anotherx;
  double me=9.1093826e-31;
  double pai=3.141592653589793; 
  double epsi=8.85418781787346e-12;
  double h=6.62606896e-34; 
  double ele=1.60217653e-19;  
  double boh=5292.0*5292.0*1.0e-28;  /* boh = ( Bohr radius )^2 */ 
     

void calc(double ca )         /* argument ca = a  */
  {
    rra=sqrt(ca*ca+r*r);       /* rra =distance between electron 1 and n1 */
    rrb=sqrt(r*r+(nuc-ca)*(nuc-ca));  /* rrb=  between electron 1 and n2 */
    rrc=sqrt(4*r*r+(nuc-2*ca)*(nuc-2*ca));  /* between two electrons */
 
    poten=ke*(-2.0/rra-2.0/rrb+1.0/nuc+1.0/rrc);   /* poten=potential energy (J) */
    ppot=poten*6.241509e18;       /* ppot = potential energy (eV) */ 
                                  
    kinetic=-0.5*poten;      /* kinetic = total kinetic energy of two electrons (J) */
    velo=sqrt(kinetic/me);  /* velo = electron's velocity (m/s) */ 
    binding=-ppot*0.5-13.606*2;       /* binding energy (eV) of H2 */

                            /* foeces acting on electron 1 ( x, -y directions ) */
   elefx= -ca/(rra*rra*rra) + (nuc-ca)/(rrb*rrb*rrb) - (nuc-2*ca)/(rrc*rrc*rrc);
   elefy= r/(rra*rra*rra) + r/(rrb*rrb*rrb) - (2*r)/(rrc*rrc*rrc);
   eleforce = ke*sqrt(elefx*elefx + elefy*elefy);  /* total force */

   radius=(me*velo*velo)/eleforce;  /* rotation radius from centrifugal force */
   debroglie = h/(me*velo);       /* de Broglie wavelength of electron */
   wave=(2*pai*radius)/debroglie;  /* wave's number in one orbit */

   nucforce= ca/(rra*rra*rra)-1.0/(nuc*nuc)+(nuc-ca)/(rrb*rrb*rrb);
   nucforce=nucforce * boh;            /* total force acting on nucleus 1 */

   anotherx=ca+(2.0*r*elefx)/elefy;    /* another x coordinate in direction of force  */

  }


void main(void) 
 {
   int i;
   double a,aa,nnuc,rr, wwave, bbind;
     
                                      /* input  nuc, x, y coordinate */
   
    printf(" Internuclear distance (MM) of H2 molecule ? ");  
    scanf("%lf",&nnuc);
 
    printf(" r (MM) ? = y coordinate of electron 1 ");  
    scanf("%lf",&rr);

    printf(" a (MM) ? = x coordinate of electron 1 ");  
    scanf("%lf", &aa);
    printf("                        \n");

    r = rr * 1.0e-14; nuc = nnuc * 1.0e-14;  /*  change MM to meter */ 
    ke = (ele*ele)/(4.0*pai*epsi);
 
   for (i=1; i < 10 ;i++) {      /* repeat until a=initial a+100 */

    a=aa*1.0e-14;
    calc(a);                    /*  to void calc()    */

   elefx=elefx * boh; elefy = elefy * boh; 
   a=anotherx; anotherx=anotherx*1.0e14;     /* a = another x  */
   printf("a:%.1f ", aa);
   printf(" binding: %.4f ", binding); 
   printf(" elefx: %.3f ", elefx);
   printf(" elefy: %.3f ", elefy);
   printf("nucforce: %.3f \n", nucforce);
   printf(" another x : %.2f ", anotherx);
  

    wwave = wave; bbind = binding;
    calc(a);                         /*  repeat calcu, when a = anotherx   */
                                    /*  average waves and binding energy  */
    wave=(wave+wwave)/2.0; binding=(binding+bbind)/2.0;

    elefx=elefx * boh;
    printf("  after-elefx: %.3f ", elefx);   /* elefx at another x */
    printf(" avebinding: %.4f ", binding);
    printf(" avewave: %.4f \n", wave);     
    printf("                        \n");

   aa=aa+10;}                             
     
   }