Einführung in Mathematica X-I

© Jürgen Schnack, 2023/04/13

Probleme mit gewöhnlichen Differentialgleichungen

Klassische Molekulardynamik - Argoncluster - Grundzustände

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

Definitionen

Teilchenzahl, Masse, Potentialparameter des Lennard-Jones-Potentials

Teilchenzahl=13;
M=39.948`;
ϵ=10.3`;
σ=3.405`;

V[r_]=4*ϵ*(((σ^12)/(r^12))-((σ^6)/(r^6)));
Plot[V[r],{r,0,10},PlotRange->{-12,12},
PlotStyle->{{Thickness[0.01],Red}},
AxesLabel->{"x","y"},
BaseStyle->{FontFamily->"Helvetica",FontSize->18}]

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 
und vordefinierter Orientierung (Warum?!)

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 +
(X[2])^2 +(X[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 (aus kleinerem Cluster + neues Atom)

(*(Subscript[X, 0])=Flatten[Append[4.0*PolyhedronData["Cuboctahedron","VertexCoordinates"],{0.0,0.0,0.0}]]*)
(Subscript[X, 0])={-3.64549,-0.402536,0.000351715,-1.51974,-1.19507,-3.12045,-1.51923,-1.19524,3.12174,-1.93696,2.4603,-1.93587,-1.93687,2.46117,1.93667,-1.26344,-3.44246,0.000100732,1.27089,3.4824,-0.00027573,1.93704,-2.46226,-1.93647,1.9383,-2.46287,1.93647,1.52548,1.20183,-3.14507,1.52742,1.20277,3.14627,3.68548,0.422335,-0.00054033,-0.0628789,-0.0389676,-0.0029622}

(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]];   (* Schrittweite beim Gradientenverfahren *)
(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];
(* Print[H[P,X]//.pxrule];Print[(Subscript[X, n])]; *)
η=Min[0.1/Max[Abs[XGrad//.pxrule]],0.01];(* das kann man diskutieren *)
(*η=0.01;*)
(*Print[η];*)
(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[Chop[(Subscript[X, n])]];

Test Schwerpunkt

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

Aha, Schwerpunkt liegt jetzt wie gewollt im Ursprung.

Plot

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

ListPointPlot3D[Orte,PlotStyle->Thickness[1.`],PlotRange->{{-10,10},{-10,10},{-10,10}},
BaseStyle->{FontFamily->"Helvetica",FontSize->18}]

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},PlotPoints->{15,15,15},Axes->True,
BaseStyle->{FontFamily->"Helvetica",FontSize->18}]

... und Schluss.