Einführung in Mathematica XI - 2

© Jürgen Schnack, 2012/04/12

Probleme mit gewöhnlichen Differentialgleichungen

Klassische Molekulardynamik - Argoncluster - Dynamik

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

Definitionen

Teilchenzahl, Potentialparameter, Plotoptionen

Teilchenzahl=6;
M=39.948`;
ϵ=10.3`;
σ=3.405`;
Rand=8.`;
PlotRand=Rand+4.`;

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: zwei 3-Cluster, SP1 bei {-Rand,0,0}, SP2 bei {+Rand,0,1}

(Subscript[X, 0])=Join[{-2.2066197828933727-Rand,-0.00335502942502179,5.669467807319734*^-7,
                   1.1066246747682547-Rand,-1.9086912132066345,0.0003251038313966591,
                   1.0999951081251174-Rand,1.913286299929468,-0.0003258801273337981},
                {-2.2066197828933727+Rand, 0.00335502942502179,5.669467807319734*^-7+1.0,
                   1.1066246747682547+Rand,1.9086912132066345,0.0003251038313966591+1.0,
                   1.0999951081251174+Rand,-1.913286299929468,-0.0003258801273337981+1.0}]

(Subscript[P, c])=60.0;
(Subscript[P, 0])={(Subscript[P, c]),0.0,0.0,(Subscript[P, c]),0.0,0.0,(Subscript[P, c]),0.0,0.0,
  -(Subscript[P, c]),0.0,0.0,-(Subscript[P, c]),0.0,0.0,-(Subscript[P, c]),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,200},MaxSteps->Infinity]

Bewegung darstellen

PlotOrte=Table[ParametricPlot3D[Table[X[3 (i-1)+j][t]/.Solution[[1]],{j,1,3}],{t,0.1`,190.`},DisplayFunction->Identity,PlotRange->{{-Rand,Rand},{-Rand,Rand},{-Rand,Rand}},
BaseStyle->{FontFamily->"Helvetica",FontSize->18}],{i,1,Teilchenzahl}];
FirstFrame=Show[PlotOrte,DisplayFunction->$DisplayFunction]

ColFrames=Table[FirstFrame,{i,1,31}];
For[i=1,i<31,i++,
ti=i*5.0;
ColFrames[[i]]=Show[Graphics3D[Join[
Table[{PointSize[0.05],Hue[0.0],Point[{X[3*(j-1)+1][ti]/.Solution[[1]],
X[3*(j-1)+2][ti]/.Solution[[1]],X[3*(j-1)+3][ti]/.Solution[[1]]}]},{j,1,Teilchenzahl/2}],
Table[{PointSize[0.05],Hue[0.6],Point[{X[3*(j-1)+1][ti]/.Solution[[1]],
X[3*(j-1)+2][ti]/.Solution[[1]],X[3*(j-1)+3][ti]/.Solution[[1]]}]},{j,Teilchenzahl/2+1,Teilchenzahl}],
{PointSize[0.00001],Hue[0.1],Point[{-PlotRand,-PlotRand,-PlotRand}]},
{PointSize[0.00001],Hue[0.1],Point[{PlotRand,PlotRand,PlotRand}]}],Axes->True]];
Print[ColFrames[[i]]]]

Export graphics

Export["col-060.gif",ColFrames,"gif",Loop->True]