/*
Microcanonical Molecular Dynamics
Narejeno po F77 programu
microcanonical molecular dynamics
1.0 Dieter W. Heermann 1985
Vhodni parametri so:
npart - število delcev (mora biti večkratnik števila 4)
side - dolžina stranice kocke v sigma enotah
den - reducirana gostota
tref - reducirana temperatura
rcoff - odrez (cutoff) potenciala v sigma enotah
h - osnovni časovni korak
irep - popravi hitrost vsak irep-ti časovni korak
istop - kje prenehaj popravljati hitrost
timemx - število integracijskih korakov
*/
function molecular_dynamics ($npart, $den, $side, $tref, $rcoff, $h, $irep, $istop=0, $timemx, $koliko, $md) {
// -------------------------------------------
// najava spremenljivk
// -------------------------------------------
// sile
$fx = array();
$fy = array();
$fz = array();
// položaji
$x = array();
$y = array();
$z = array();
// hitrosti
$vh = array();
// za izpis podatkov
$podatki = array();
// preveri, da ni preveč korakov (max 2500)
if ($timemx > 2500) $timemx = 2500;
// popravi število delcev, da je deljivo s 4
$npart = $npart - ($npart % 4);
if ($npart < 4) $npart = 4;
if($md) {
echo "-- Microcanonical Molecular Dynamics --\n\n";
echo "Število delcev je: ".$npart."\n";
echo "Velikost stranice kocke je: ".$side."\n";
echo "Odrez LJ potenciala je: ".$rcoff."\n";
echo "Število korakov je: ".$timemx."\n";
echo "Osnovni časovni korak je: ".$h."\n";
echo "Popravim hitrosti vsak ".$irep."-ti korak\n";
echo "Preneham popravljat po ".$istop." korakih\n";
}
else {
echo "-- Brownian Dynamics --\n\n";
echo "Število delcev je: ".$npart."\n";
echo "Velikost stranice kocke je: ".$side."\n";
echo "Odrez LJ potenciala je: ".$rcoff."\n";
echo "Število korakov je: ".$timemx."\n";
echo "Osnovni časovni korak je: ".$h."\n";
echo "Popravim hitrosti vsak ".$irep."-ti korak\n";
}
// -------------------------------------------
// nastavi parametre
// -------------------------------------------
$a = $side / 4.0;
$sideh = $side * 0.5;
$hsq = $h*$h;
$hsq2 = $hsq * 0.5; // zakaj ne $hsq/2.0 ?
$rcoffs = $rcoff * $rcoff;
$tscale = 16.0 / (1.0 * $npart - 1.0);
$vaver = 1.13 * sqrt($tref / 24.0);
$n3 = 3 * $npart;
$iof1 = $npart; // odmik za for zanke
$iof2 = 2 * $npart; // odmik za for zanke
// -------------------------------------------
// nastavi začetno konfiguracijo
// -------------------------------------------
// postavi atome v kocki
$ijk = 0;
for ($lg = 0; $lg <= 3; $lg++) { //za $lg 0 do 3
for ($i = 0; $i <= 3; $i++) {
for ($j = 0; $j <= 3; $j++) {
for ($k = 0; $k <= 3; $k++) {
$ijk++;
if ($lg < 2) {
$x[$ijk] = $i * $a + $lg * $a * 0.5; // $lg je tukaj 0, 1
$y[$ijk] = $j * $a + $lg * $a * 0.5; // $lg je tukaj 0, 1
$z[$ijk] = $k * $a;
}
else {
$x[$ijk] = $i * $a + (3 - $lg) * $a * 0.5; // $lg izraz je tukaj 1, 0
$y[$ijk] = $j * $a + ($lg - 2) * $a * 0.5; // $lg izraz je tukaj 0, 1
$z[$ijk] = $k * $a + $a * 0.5;
}
} //$k
} //$j
} //$i
} //$lg
// določi hitrosti, porazdeljene normalno
$vh = mxwell($n3, $h, $tref);
// počisti sile
for ($i = 1; $i <= $npart; $i++) {
$fx[$i] = 0.0; // popravi polje $f na 0 (mogoče bi to naredil direkt)
$fy[$i] = 0.0;
$fz[$i] = 0.0;
}
// -------------------------------------------
// Začni dejansko molekularno dinamiko
//
// Enačbe gibanja so integrirane z uporabo "sumarne oblike".
// Vsak irep-ti korak so hitrosti popravljene, da dajo določeno reducirano temperaturo (tref)
// -------------------------------------------
echo "\n\n- Začnem računat MD -\n";
for ($clock = 1; $clock <= $timemx; $clock++) {
// napreduj položaj za en osnovni časovni korak
for ($i = 1; $i <= $npart; $i++) {
$x[$i] = $x[$i] + $vh[$i] + $fx[$i];
$y[$i] = $y[$i] + $vh[$i+$iof2] + $fy[$i];
$z[$i] = $z[$i] + $vh[$i+$iof3] + $fz[$i];
}
// dodaj perodične mejne pogoje (da delci "krožijo" v škatli)
for ($i = 1; $i <= $npart; $i++) {
if ($x[$i] < 0) $x[$i] = $x[$i] + $side;
if ($x[$i] > $side) $x[$i] = $x[$i] - $side;
if ($y[$i] < 0) $y[$i] = $y[$i] + $side;
if ($y[$i] > $side) $y[$i] = $y[$i] - $side;
if ($z[$i] < 0) $z[$i] = $z[$i] + $side;
if ($z[$i] > $side) $z[$i] = $z[$i] - $side;
}
// izračunaj delne hitrosti
for ($i = 1; $i <= $npart; $i++) {
$vh[$i] = $vh[$i] + $fx[$i];
$vh[$i+$iof1] = $vh[$i+$iof1] + $fy[$i];
$vh[$i+$iof2] = $vh[$i+$iof2] + $fz[$i];
}
// ta del izračuna sile na delcih ...
$vir = 0.0;
$epot = 0.0;
for ($i = 1; $i <= $npart; $i++) {
$fx[$i] = 0.0; // popravi polje $f na 0 (mogoče bi to naredil direkt)
$fy[$i] = 0.0;
$fz[$i] = 0.0;
}
for ($i = 1; $i <= $npart; $i++) {
$xi = $x[$i];
$yi = $y[$i];
$zi = $z[$i];
for ($j = $i+1; $j < $npart; $j++) {
$xx = $xi - $x[$j];
$yy = $yi - $y[$j];
$zz = $zi - $z[$j];
if ($xx < $sideh*-1) $xx = $xx + $side;
if ($xx > $sideh) $xx = $xx - $side;
if ($yy < $sideh*-1) $yy = $yy + $side;
if ($yy > $sideh) $yy = $yy - $side;
if ($zz < $sideh*-1) $zz = $zz + $side;
if ($zz > $sideh) $zz = $zz - $side;
$rd = $xx * $xx + $yy * $yy + $zz * $zz;
if ($rd > $rcoffs) break 2; // pojdi ven iz obeh for loopov
$epot = $epot + pow($rd, (-6.0)) - pow($rd, (-3.0));
$r148 = pow($rd, (-7.0)) - 0.5 * pow($rd, (-4.0));
$vir = $vir - $rd * $r148;
$kx = $xx * $r148;
$fx[$i] = $fx[$i] + $kx;
$fx[$j] = $fx[$j] - $kx;
$ky = $yy * $r148;
$fy[$i] = $fy[$i] + $ky;
$fy[$j] = $fy[$j] - $ky;
$kz = $zz * $r148;
$fz[$i] = $fz[$i] + $kz;
$fz[$j] = $fz[$j] - $kz;
} //$i
} //$j
for ($i = 1; $i <= $npart; $i++) {
$fx[$i] = $fx[$i] * $hsq2;
$fy[$i] = $fy[$i] * $hsq2;
$fz[$i] = $fz[$i] * $hsq2;
}
// ... konec računanja sil na delcih
// izračunaj hitrosti
for ($i = 1; $i <= $npart; $i++) {
$vh[$i] = $vh[$i] + $fx[$i];
$vh[$i+$iof1] = $vh[$i+$iof1] + $fy[$i];
$vh[$i+$iof2] = $vh[$i+$iof2] + $fz[$i];
}
// izračunaj kinetično energijo
$ekin = 0.0;
for ($i = 1; $i <= $n3; $i++) {
$ekin = $ekin + $vh[$i] * $vh[$i];
}
$ekin = $ekin / $hsq;
// izračunaj povprečno hitrost
$vel = 0.0;
$count = 0.0;
for ($i = 1; $i <= $npart; $i++) {
$vx = $vh[$i] * $vh[$i];
$vy = $vh[$i+$iof1] * $vh[$i+$iof1];
$vz = $vh[$i+$iof2] * $vh[$i+$iof2];
$sq = sqrt($vx + $vy + $vz);
$sqt = $sq / $h;
if ($sqt > $vaver) $count = $count + 1;
$vel = $vel + $sq;
}// $i
$vel = $vel / $h;
// izberi tip simulacije
if ($md) {
if ($clock < $istop) {
// normaliziraj hitrosti da dobiš nakazano referenčno temperaturo
if (($clock % $irep) == 0) {
echo "\n Korak ".$clock."\n";
echo "Popravljam hitrosti ...\n";
$ts = $tscale * $ekin;
echo "Temperatura pred spremembo je ".$ts."\n";
$sc = $tref / $ts;
$sc = sqrt($sc);
echo "Faktor spremembe je ".$sc."\n";
for ($i = 1; $i <= $n3; $i++) {
$vh[$i] = $vh[$i] * $sc;
}
$ekin = $tref / $tscale; // to je zmeraj enako
}
} // if clock
}
// simulacija je BD
else {
// normaliziraj hitrosti da dobiš nakazano referenčno temperaturo
if (($clock % $irep) == 0) {
echo " Korak ".$clock."\n";
echo "Zamenjam hitrosti ...\n";
$vh = mxwell($n3, $h, $tref);
$ekin = $tref / $tscale;
}
}
// izračunaj razne količine
if (($clock % $koliko) == 0) {
$ek = 24.0 * $ekin;
$epot = 4.0 * $epot;
$etot = $ek + $epot;
$temp = $tscale * $ekin;
$pres = $den * 16.0 * ($ekin - $vir) / $npart;
$vel = $vel / $npart;
//$rp = ($count / $npart) * 100.0;
/*
echo "korak ".$clock."\n";
echo "kin. energija: ".$ek."\n";
echo "pot. energija: ".$epot."\n";
echo "tot. energjia: ".$etot."\n";
echo "temperatura: ".$temp."\n";
echo "pritisk: ".$pres."\n";
echo "povp. hitrost: ".$vel."\n";
//echo "rp: ".$rp."\n\n";
*/
$ek_out = $ek_out. "[".$clock.", ".$ek."], ";
$epot_out = $epot_out. "[".$clock.", ".$epot."], ";
$etot_out = $etot_out. "[".$clock.", ".$etot."], ";
$temp_out = $temp_out. "[".$clock.", ".$temp."], ";
$pres_out = $pres_out. "[".$clock.", ".$pres."], ";
$vel_out = $vel_out. "[".$clock.", ".$vel."], ";
}
} //$clock
$podatki["ek"] = $ek_out;
$podatki["epot"] = $epot_out;
$podatki["etot"] = $etot_out;
$podatki["temp"] = $temp_out;
$podatki["pres"] = $pres_out;
$podatki["vel"] = $vel_out;
return $podatki;
} //molecular_dynamics()