Gleichgewichtsdynamik eines Coulomb-Gases im harmonischen Oszillator

© Jürgen Schnack, 2019/10/22

Klassische Molekulardynamik - Argoncluster - Dynamik

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

Definitionen

Teilchenzahl, Potentialparameter, Plotoptionen

Teilchenzahl=3;
M=39.948`;
ϵ=10.3`;
σ=3.405`;
Rand=10.`;
PlotRand=Rand;
TMax=5000;
Temperatur=10.0;
(Subscript[M, η])=3.1415925;

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];
Array[η[t],1];

Hamilton-Funktion in P und X

H[P_,X_]= (∑_(i=1)^(3*Teilchenzahl))(P[i][t]*P[i][t])      +1.0*
(∑_(i=1)^Teilchenzahl)((∑_(j=i+1)^Teilchenzahl)((1.0/((((((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)))^(1/2)))))+
1.0*(∑_(i=1)^(3*Teilchenzahl))(X[i][t]*X[i][t])

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]]-η[1][t]*P[i][t],{i,1,3*Teilchenzahl}];

EOMX=Table[(X[i])'[t]==PGrad[[i]],{i,1,3*Teilchenzahl}];

EOMT=Table[(η[k])'[t]==(1/(Subscript[M, η]))( (∑_(i=1)^(3*Teilchenzahl))(P[i][t]*P[i][t])-1.5*Teilchenzahl*Temperatur),{k,1,1}];

Anfangsbedingungen:

Ort0={{0.0,0.0,-2.0},{1.0,0.0,0.0},{0.0,1.0,0.0}};

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

(Subscript[P, c])=10.0*Sqrt[Temperatur];
Impulse0=(Subscript[P, c])*Table[{RandomReal[]-0.5,RandomReal[]-0.5,RandomReal[]-0.5},{k,1,Teilchenzahl}];
MeanP=Table[Sum[Impulse0[[k,i]],{k,1,Teilchenzahl}]/Teilchenzahl,{i,1,3}];
For[k=1,k<Teilchenzahl+1,k++,
Impulse0[[k,1]]=Impulse0[[k,1]]-MeanP[[1]];Impulse0[[k,2]]=Impulse0[[k,2]]-MeanP[[2]];Impulse0[[k,3]]=Impulse0[[k,3]]-MeanP[[3]]];
(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}];
ABT=Table[η[i][0]==0.0,{i,1,1}];

Alle Gleichungen

Equations=Join[EOMP,EOMX,EOMT,ABP,ABX,ABT];

Alle Impulse und Koordinaten

AlleP=Table[P[i],{i,1,3*Teilchenzahl}];
AlleX=Table[X[i],{i,1,3*Teilchenzahl}];
AlleT=Table[η[i],{i,1,1}];

AlleKoordinaten=Join[AlleP,AlleX,AlleT];

DGL lösen

Solution=NDSolve[Equations,AlleKoordinaten,{t,0,TMax},MaxSteps->Infinity];

Bewegung darstellen

Print[ParametricPlot3D[Evaluate[Table[X[3 (i-1)+j][t],{i,1,Teilchenzahl},{j,1,3}]/.Solution[[1]]],{t,0.1`,TMax},
DisplayFunction->Identity,PlotRange->{{-Rand,Rand},{-Rand,Rand},{-Rand,Rand}},
BaseStyle->{FontFamily->"Helvetica",FontSize->18}]];

PlotEta=Plot[η[1][t]/.Solution[[1]],{t,0.1`,TMax},BaseStyle->{FontFamily->"Helvetica",FontSize->18}]

Mittlere kinetische Energie (Was erwartet man?)

EkinMittel=Sum[(∑_(i=1)^(3*Teilchenzahl))(P[i][ti]*P[i][ti])/.Solution[[1]],{ti,1,TMax}]/TMax

EkinMittel-1.5*Teilchenzahl*Temperatur

Mittlere potentielle Energie im harmonischen Oszillator (Was erwartet man?)

EHOMittel=Sum[(∑_(i=1)^(3*Teilchenzahl))(X[i][ti]*X[i][ti])/.Solution[[1]],{ti,1,TMax}]/TMax