Kernphysik: Spaltung (nur zur Veranschaulichung)

© Jürgen Schnack, 2019/10/14

Klassische Molekulardynamik - Argoncluster - Dynamik

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

Definitionen

Teilchenzahl, Potentialparameter, Plotoptionen

Teilchenzahl=48;
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:

OrtA={-3.3953477844147337`,0.07219186654122567`,1.1443747097698573`,-2.1832565376120403`,-2.0766286266854466`,-4.400974058664173`,-0.19583736110049624`,-1.872677535394477`,1.7046726055572645`,-1.9993856529072038`,1.9729152633359246`,-4.453904430501098`,-0.13974667432058782`,1.8698419502981813`,1.7301511166652925`,-1.1765099058303183`,-5.504563782855491`,0.901432778295945`,-0.7156450977674557`,5.4824922261125755`,0.6531641389926675`,0.6892372375048098`,-3.2977884757039306`,-1.5562940219891077`,3.7933834282313557`,-1.7401005545236115`,2.563971428436241`,0.9859081552797554`,3.0075151195582572`,-1.6886034812205093`,3.4041595510480556`,2.01561274826802`,2.7170909420779203`,6.8747639224952755`,0.5072505544474243`,2.6470738272982253`,-0.9096160435651298`,-0.07727912013785714`,-1.5425572002618901`,1.1566920484144045`,-0.16412413773288947`,4.714843002092007`,-3.0144906898284436`,-2.915585305575937`,-0.9243024637298645`,-2.7790275921407313`,2.993708735421117`,-0.9928701683712257`,2.534770142545313`,-0.10543765685376816`,-0.4081636608581103`,-2.4634918548301226`,-8.633616828904627`,-0.7298256365846633`,-1.9563820160648695`,8.58906497035623`,-1.0391341018877156`,1.0039863504089486`,-4.084471601521381`,-5.120508675506795`,1.3639261301664591`,3.643757978433851`,-5.268203038067103`,1.9194188288898646`,-4.926121710571862`,2.7594237758326847`,1.354589164881452`,5.037609245659121`,3.6586442234194956`,2.929873637889217`,-3.2465103847539045`,5.964525733033533`,2.6838634969747646`,2.893752718288783`,6.396969561191676`,3.797730700033279`,-2.20674876726963`,-3.370200940300861`,4.002690614835908`,1.5514548048412848`,-3.4865730690149195`,-4.640576121510384`,0.08667699739058811`,-2.433248496155547`,-0.5745378734848432`,-2.032970205053393`,-7.7868129052065775`,-0.4491930087778154`,1.6697267174968702`,-7.8408173040581755`,4.845508231363351`,-0.030370699170077998`,5.726948624627625`,3.9896997458824957`,-4.246537554960656`,-0.23180519708411282`,3.0705116499436973`,5.148646640884568`,0.40050620143527793`,0.9943976891044248`,-0.19905439838176978`,-4.713863189224971`,-4.936723882458472`,-5.849290709662667`,0.10763900015703462`,-4.518446593244514`,6.012841201653011`,0.1121962360647093`,-6.390449332314956`,-1.9652360173344203`,0.32605890844677254`,-6.209804869367994`,2.330575149162433`,0.30314868614211243`,6.093306939538736`,-1.1967925724556527`,-0.5911971984929649`,5.515499412917931`,2.4569119818992404`,-0.27496908498013495`,-2.5518923797825774`,-0.30139335641462156`,4.620953605290757`,-3.679244831236425`,-3.3352978580559633`,2.706537865121849`,-3.363718397221114`,3.4714901903217545`,2.6945722178977465`,-1.8062372838035627`,-5.688688226099672`,-3.0094423887674893`,-1.2838884020460826`,5.551748832169752`,-3.19976780171883`,-0.7126873602908802`,-3.545773138460912`,5.001210207469414`,-0.9577795324277445`,2.8757029302960593`,5.483580168475712`,-2*Rand,0.0,0.0};

Ort0=Table[OrtA[[(i-1)*3+j]],{i,1,Teilchenzahl},{j,1,3}]

(* 
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])=120.0;
Impulse0=Table[{0.0,0.0,0.0},{k,1,Teilchenzahl}];
For[k=1,k<Teilchenzahl,k++,If[Ort0[[k,3]]>0,Impulse0[[k,3]]=(Subscript[P, c]),Impulse0[[k,3]]=-(Subscript[P, c])]];
Impulse0[[Teilchenzahl,1]]=1.0*(Subscript[P, c]);
(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,201}];
For[i=1,i<201,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-1}],
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,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]