Einführung in Mathematica XI - 2

© Jürgen Schnack, 2015/11/25

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=26;
M=39.948`;
ϵ=10.3`;
σ=3.405`;
Rand=8.`;
Stossparameter=0.0;
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,Stossparameter}

Ort0={{-3.645491367896175`,-0.4025355059440745`,0.00035171472177316894`},{-1.5197412540850084`,-1.1950704380347066`,-3.1204494059839227`},{-1.5192289177526144`,-1.1952419605765854`,3.121738245661295`},{-1.9369628205990985`,2.4602961122749374`,-1.9358697043990056`},{-1.9368733422316688`,2.461165032453709`,1.936674524541644`},{-1.2634413608573518`,-3.442463085419254`,0.0001007317665096586`},{1.2708934526771403`,3.4824001614755664`,-0.0002757303556482697`},{1.9370365366408806`,-2.4622641134334913`,-1.9364683131311224`},{1.938304285920814`,-2.462872795907553`,1.9364667455181817`},{1.5254783034436967`,1.201828495033298`,-3.1450692426903046`},{1.527421327794014`,1.2027686165909806`,3.146266587118005`},{3.685484087443547`,0.42233530440944306`,-0.0005403296620990516`},{-0.0628789304981757`,-0.038967563852657414`,-0.0029622016378999727`},{-3.645491367896175`,-0.4025355059440745`,0.00035171472177316894`},{-1.5197412540850084`,-1.1950704380347066`,-3.1204494059839227`},{-1.5192289177526144`,-1.1952419605765854`,3.121738245661295`},{-1.9369628205990985`,2.4602961122749374`,-1.9358697043990056`},{-1.9368733422316688`,2.461165032453709`,1.936674524541644`},{-1.2634413608573518`,-3.442463085419254`,0.0001007317665096586`},{1.2708934526771403`,3.4824001614755664`,-0.0002757303556482697`},{1.9370365366408806`,-2.4622641134334913`,-1.9364683131311224`},{1.938304285920814`,-2.462872795907553`,1.9364667455181817`},{1.5254783034436967`,1.201828495033298`,-3.1450692426903046`},{1.527421327794014`,1.2027686165909806`,3.146266587118005`},{3.685484087443547`,0.42233530440944306`,-0.0005403296620990516`},{-0.0628789304981757`,-0.038967563852657414`,-0.0029622016378999727`}};

For[k=1,k<Teilchenzahl/2+1,k++,Ort0[[k,1]]=Ort0[[k,1]]-Rand];
For[k=1,k<Teilchenzahl/2+1,k++,Ort0[[k,3]]=Ort0[[k,3]]+Stossparameter];
For[k=Teilchenzahl/2+1,k<Teilchenzahl+1,k++,Ort0[[k,1]]=Ort0[[k,1]]+Rand];

(Subscript[X, 0])=Flatten[Ort0]

(Subscript[P, c])=60.0;
Impulse0=Table[{(Subscript[P, c]),0.0,0.0},{k,1,Teilchenzahl}];
For[k=Teilchenzahl/2+1,k<Teilchenzahl+1,k++,Impulse0[[k,1]]=-1.0*Impulse0[[k,1]]]
(Subscript[P, 0])=Flatten[Impulse0]

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,300},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,251}];
For[i=1,i<251,i++,
ti=i*1.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]]];

ListAnimate[ColFrames]