ヒュッケル分子軌道法の計算プログラム - PHP版
化合物名 : " . $title . "");
print("原子数 : ". $na ."");
print("環状 : ". $cyclo ."");
// 永年行列
// Zero Clear
for($i=0 ; $i < $na ; $i++) {
for($j=0 ; $j < $na; $j++) {
$h[$i][$j] = 0;
$v[$i][$j] = 0;
}
}
// 結合している原子の番号の組 (i,j)
for($i=0 ; $i < $na-1; $i++) {
$h[$i][$i+1] = 1.0 ;
$h[$i+1][$i] = 1.0 ;
}
if(($na > 2) && ($cyclo == 1 )) {
$h[0][$na-1] = 1.0 ;
$h[$na-1][0] = 1.0 ;
}
// ヤコビ法による行列の対角化
// var $sm, $hj, $hk, $ta, $tb ;
// var $sn, $cs, $hl, $hjl, $vj, $vk ;
define("TH", 0.00001); // 収束判定のしきい値
define("ITMAX",100); // 収束判定のしきい値
$it = 0 ;
$d = 1.0 ;
// 単位行列
for($i=0; $i < $na; $i++) {
$v[$i][$i] = 1.0 ;
}
// 収束判定
while (($d > TH) && ($it < ITMAX)) {
$d = abs($h[1][0]);
for($i=1; $i < $na ; $i++) {
for($j=0; $j < $i; $j++) {
if(abs($h[$i][$j]) > $d) $d = abs($h[$i][$j]);
}
}
if($d > TH) {
$sm = 0.1*$d ;
$it++ ;
for($j=1; $j < $na; $j++) {
for($k=0; $k < $j; $k++) {
if(abs($h[$j][$k]) > $sm) {
$hj = $h[$j][$j] ;
$hk = $h[$k][$k] ;
$tb = ($hj - $hk)/(2.0*$h[$k][$j]);
if($tb <= 0)
$ta = $tb - sqrt(1.0+$tb*$tb);
else
$ta = $tb + sqrt(1.0+$tb*$tb);
$cs = 1.0/sqrt(1.0+$ta*$ta);
$sn = $cs*$ta ;
for($l=0; $l < $na ; $l++) {
if($l == $k) {
$h[$k][$k] = $hk*$cs*$cs+$h[$j][$k]*2.0*$sn*$cs+$hj*$sn*$sn ;
$h[$j][$j] = $hk+$hj-$h[$k][$k] ;
}
elseif($l != $j) {
$hkl = $h[$k][$l] ;
$hjl = $h[$j][$l] ;
$h[$k][$l] = $hkl*$cs+$hjl*$sn ;
$h[$j][$l] = $hjl*$cs-$hkl*$sn ;
$h[$l][$k] = $h[$k][$l] ;
$h[$l][$j] = $h[$j][$l];
}
}
$h[$j][$k] = $h[$j][$k]*($cs*$cs-$sn*$sn)+($hj-$hk)*$sn*$cs ;
$h[$k][$j] = $h[$j][$k] ;
for($i=0; $i < $na ; $i++) {
$vj = $v[$i][$j] ;
$vk = $v[$i][$k] ;
$v[$i][$j] = $vj*$cs - $vk*$sn ;
$v[$i][$k] = $vj*$sn + $vk*$cs ;
}
}
}
}
}
}
// 固有値と固有ベクトルを固有値の昇順に並べかえる
// var $imax, $emax, $temp;
for($i=0; $i < $na; $i++) {
$emax = $h[$i][$i] ;
$imax = $i ;
for($j=$i; $j < $na ; $j++) {
if($h[$j][$j] > $emax) {
$imax = $j ;
$emax = $h[$j][$j] ;
}
}
if($i < $imax) {
$temp = $h[$i][$i];
$h[$i][$i]= $h[$imax][$imax] ;
$h[$imax][$imax] = $temp ;
for($j=0; $j < $na; $j++) {
$temp = $v[$j][$i];
$v[$j][$i]= $v[$j][$imax] ;
$v[$j][$imax] = $temp ;
}
}
}
// 固有値と固有ベクトルの印刷
define("COLUMN", 6); // 一行に印刷する項目数
$first = 0 ;
print(" -- Eigen values and Eigen vector --
");
while ($first < $na) {
$last = $first + COLUMN - 1 ;
print("
");
if($last > $na-1 ) $last = $na-1 ;
print(" ");
for($i=$first; $i <= $last; $i++){
printf("%10.4f",$h[$i][$i]);
}
print("
");
for($k=0; $k < $na ; $k++) {
printf("%2d", $k+1 );
for($j=$first; $j <= $last; $j++) {
printf("%10.4f",$v[$k][$j] );
}
print("
");
}
print("
");
$first = $last + 1 ;
}
print("
");
?>