Ü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, 3*Teilchenzahl];
Array[X, 3*Teilchenzahl];

Hamilton-Funktion in P und X mit Schwerpunkt bei Null(!)

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

Definiere Ableitungsfelder

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

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

Grundzustand

Startkonfiguration

(Subscript[X, 0]) = {-1.0, 0.0, 0.0, 1.0, 0.0, 0.0}

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

Definiere Ersetzungsregeln: zum konkreten Ausrechnen von H oder seinen Gradienten
müssen den Variablen P und X konkrete Zahlen zugewiesen werden

prule = 
        Table[P[i] -> (Subscript[P, 0])[[i]], {i, 1, 
            3*Teilchenzahl}];
xrule = 
        Table[X[i] -> (Subscript[X, 0])[[i]], {i, 1, 
            3*Teilchenzahl}];
pxrule = Union[prule, xrule];

H[P, X] //. pxrule

PGrad //. pxrule

XG = XGrad//.pxrule

Kraft = 
      Table[ (∑_(i = 1)^Teilchenzahl)(XG[[3*(i - 1) + 
                  j]]),   {j, 1, 3}]

Primitives Gradientenverfahren

η = 1/Max[Abs[XGrad //. pxrule]];     
(Subscript[P, n]) = 
        Table[(Subscript[P, 0])[[i]], {i, 1, 
            3*Teilchenzahl}];
(Subscript[P, n + 1]) = 
        Table[0, {i, 1, 3*Teilchenzahl}];
(Subscript[X, n]) = 
        Table[(Subscript[X, 0])[[i]], {i, 1, 
            3*Teilchenzahl}];
(Subscript[X, n + 1]) = 
        Table[0, {i, 1, 
            3*Teilchenzahl}];

For[k = 1, 
        k < 1001, k++, 
prule = 
          Table[P[i] -> (Subscript[P, n])[[i]], {i, 1, 
              3*Teilchenzahl}]; 
xrule = 
          Table[X[i] -> (Subscript[X, n])[[i]], {i, 1, 
              3*Teilchenzahl}]; 
pxrule = 
          Union[prule, xrule]; 
  
η = 
          Min[0.1/Max[Abs[XGrad //. pxrule]], 
            0.01]; 
  
(Subscript[X, n + 1]) = 
          (Subscript[X, n]) - η*XGrad //. pxrule; 
(Subscript[X, n]) = 
          (Subscript[X, n + 1])
];
prule = 
        Table[P[i] -> (Subscript[P, n])[[i]], {i, 1, 
            3*Teilchenzahl}];
xrule = 
        Table[X[i] -> (Subscript[X, n])[[i]], {i, 1, 
            3*Teilchenzahl}];
pxrule = Union[prule, xrule];
Print[H[P, X] //. pxrule] ; Print[
        (Subscript[X, n])] ;

Test Schwerpunkt

Schwerpunkt = 
      Table[ (1/Teilchenzahl)    (∑_(i = 1)^Teilchenzahl)((Subscript[X, n])[[(i - 1)*3 + j]]) , {j, 1, 
          3}]

Plot

Orte = 
      Table[(Subscript[X, n])[[(i - 1)*3 + j]], {i, 1, Teilchenzahl}, {j, 1, 
          3}]

ScatterPlot3D[Orte, PlotStyle -> Thickness[1.0], 
        PlotRange -> {{-10, 10}, {-10, 10}, {-10, 
              10}}];

Smoozy[x_, y_, z_] = 
      10*(∑_(i = 1)^Teilchenzahl)((E^(-((((x - (Subscript[X, n])[[(i - 1)*3 + 1]]))^2) +    (((y - (Subscript[X, n])[[(i - 1)*3 + 2]]))^2) +   (((z - (Subscript[X, n])[[(i - 1)*3 + 3]]))^2)  ))))

ContourPlot3D[
        Smoozy[x, y, z], {x, -10, 10}, {y, -10, 10}, {z, -10, 
          10}, 
Contours -> {0.5}, 
        Maxecursions -> 10, PlotPoints -> {15, 15, 15}, 
        Axes -> True];