下のソースプログラムをそのままテキストエディタ（メモ帳など）にコピー and ペースト すれば、簡単にコンパイルと実行できる。
(このプログラム は 単純なＣ言語 なので、このテキストエディタを "filename.c" とセーブしてコンパイルしてほしい。)
ここでは 1 MM = 1 × 10-14 meter という新しい単位を使用している。
このサンプルプログラムでは、実行すると、最初に H2 の核間距離 ( MM ) を入力する。
その後、電子１の ｙ座標 "r" ( MM )、x 座標 "a" ( MM ) を 入力する。
これら入力値から、このプログラムは H2 の結合エネルギー (eV)、原子核 1 に作用する トータルの力 (= nucforce )、電子１に作用する力の x, y 成分 (= elefx, elefy ) を出力する。
これらの力は 水素原子の基底状態の電子と原子核間の力との比で表されている。
電子に作用する力と その速度から １軌道に含まれるドブロイ波を計算する。
( ここでは ビリアル定理 E = -T = 1/2 V を使用している。)
"another x" は 電子１が 力の方向へ移動した先の もう１つの x 座標である。
そして 今度は "a" の代わりに "another x" を x 座標として上記の計算を繰り返す。
"avebinding" は x が "a" と "another x" のときのを 平均した結合エネルギーである。
"avewave" は a と another 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);
printf(" wave: %.4f ", wave);
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;
printf(" avebinding: %.4f ", binding);
printf(" avewave: %.4f \n", wave);
printf(" \n");
aa=aa+10;}
}