Übung 6

Klassische Molekulardynamik - Argoncluster

Energien in meV,
Längen in A,
Geschwindigkeiten in A/ps,
Massen in u

Definitionen

Teilchenzahl

Teilchenzahl = 2;
M = 39.948;
ϵ = 10.3;
σ = 3.405;
<< Graphics`Graphics3D`
<< Graphics`ContourPlot3D`

Hamilton-Funktion und verallgemeinerte Gradienten in P und X

P und X werden als Array definiert, damit weiß Mathematica, was das für Objekte sind

Array[P[t], 3*Teilchenzahl];
Array[X[t], 3*Teilchenzahl];

Hamilton-Funktion in P und X

H[P_, X_] = 
        (0.1036446/(2 M)) (∑_(i = 1)^(3*Teilchenzahl))(P[
                      i][t]*P[i][t])       + 
          4*ϵ*
(∑_(i = 1)^Teilchenzahl)((∑_(j = i + 1)^Teilchenzahl)(((σ^12)/((((((X[3*i - 2][t] - X[3*j - 2][t]))^2) +   (((X[3*i - 1][t] - X[3*j - 1][t]))^2) + (((X[3*i][t] - X[3*j][t]))^2)))^6)) - 
((σ^6)/((((((X[3*i - 2][t] - X[3*j - 2][t]))^2) +   (((X[3*i - 1][t] - X[3*j - 1][t]))^2) + (((X[3*i][t] - X[3*j][t]))^2)))^3))));

Definiere Ableitungsfelder

PGrad = 
      Table[D[H[P, X], P[i][t]], {i, 1, 3*Teilchenzahl}]

XGrad = 
        Table[D[H[P, X], X[i][t]], {i, 1, 3*Teilchenzahl}];

Dynamik

Differentialgleichungen

EOMP=Table[(P[i])'[t]==-XGrad[[i]],{i, 1, 3*Teilchenzahl}];

EOMX=Table[(X[i])'[t]==PGrad[[i]],{i, 1, 3*Teilchenzahl}];

Anfangsbedingungen

(Subscript[X, 0]) = {-5.0, -2.0, 0.0, 5.0, 2.0, 
          0.0};
(Subscript[P, 0]) = {40.0, 0.0, 0.0, -40.0, 0.0, 
          0.0};

ABP = 
        Table[P[i][0] == (Subscript[P, 0])[[i]], {i, 1, 
            3*Teilchenzahl}];
ABX = 
        Table[X[i][0] == (Subscript[X, 0])[[i]], {i, 1, 
            3*Teilchenzahl}];

Alle Gleichungen

Equations = Join[EOMP, EOMX, ABP, ABX];

Alle Impulse und Koordinaten

AlleP = 
        Table[P[i], {i, 1, 3*Teilchenzahl}];
AlleX = 
        Table[X[i], {i, 1, 
            3*Teilchenzahl}];

AlleKoordinaten = Join[AlleP, AlleX]

DGL lösen

Solution = NDSolve[Equations, AlleKoordinaten, {t, 0, 1000}]

Bewegung darstellen

P2 = 
        ParametricPlot3D[{{X[1][t] /. Solution[[1]], X[2][t] /. 
                Solution[[1]], X[3][t] /. 
                Solution[[1]]}, 
{X[4][t] /. 
                Solution[[1]], X[5][t] /. 
                Solution[[1]], X[6][t] /. Solution[[1]]}}, {t, 
            0.1, 80.0}, 
PlotRange -> {{-5, 
                5}, {-5, 5}, {-5, 5}}];

PlotOrte = 
        Table[ParametricPlot3D[
            Table[X[3*(i - 1) + j][t] /. Solution[[1]], {j, 1, 
                3}], {t, 0.1, 
              80.0}, 
DisplayFunction -> Identity, 
            PlotRange -> {{-5, 5}, {-5, 5}, {-5, 5}}], {i, 
            1, Teilchenzahl}];
Show[PlotOrte, 
        DisplayFunction -> $DisplayFunction];