Informations

Projet 				: 	EPFL  Notebooks de démonstrations
Auteur		 		: 	Amina Matt
Date de création 		:	Novembre 2018
But					: 	Illustration du potentiel de Lennard-Jones
Historique des révisions	:	L’option “Index” dans la fonction ‘’Nearest’’ est disponible uniquement à partir de la version 11.1

Comment utiliser le Notebook?

Mac Os:  Pour évaluer toutes les cellules du notebook utiliser ⌘+A et  ⇧+↵ . Pour fermer la fenêtre de visualisation utiliser ⌘+W.

Windows:  Pour évaluer toutes les cellules du notebook utiliser Ctrl+A et  ⇧+↵ . Pour fermer la fenêtre de visualisation utiliser Alt+F4.

Pour voir le code utilisé double-cliquer sur le symbole  ou correspondant à la section souhaitée. Pour la refermer double cliquer sur   correspondant.

Pour en savoir plus sur la méthode de simulation utilisée cliquer sur les liens.

Le potentiel de Lennard-Jones : Atomes et structures

Théorie

Pour plus de détails, voir slides du Cours 3.1 Structure des matériaux

Lennard - Jones

The Lennard-Jones potential describes the interaction between atom, i.e. the energy depending on the distance. 
It is written as E(r)= ϵ ((a/(r^12))-(b/(r^6)))
it can be normalized and written as ((E(r))/(Subscript[E, min])) =((((Subscript[r, min])^12)/(r^12))-((2((Subscript[r, min])^6))/(r^6)))
using ρ =((Subscript[r, min])/r)  E(ρ) =((1/(ρ^12))-(2/(ρ^6)))=(1/(ρ^12))(1-2(ρ^6))

the forces then F(ρ)=((dE(ρ))/(dρ))= -12 ρ(((1-(ρ^6))/(ρ^14)))

Dynamique Moléculaire

Molecular dynamics (MD) is a method that allows computation of a system evolution, a system depending on time. 
Let’s explain what are the steps of a MD algorithm.

STEP 0. The system has an initial state, i.e. positions and initial velocities for each component.
STEP 1. The forces between the particles and the energies are calculated with these positions and velocities (this is where you need LJ potential)
STEP 2. The system is evolved in time, i.e. new positions and velocities, using these forces and energies (several formula exist, time step)
STEP 3. You can repeat 1-2 until the system (almost) does not change anymore (this is the tolerance criteria)

Critical input parameters

*Potential

*Updating algorithm

*Time step

*Tolerance

Code

Initialisation

Initialisation globale

Le potentiel de Lennard-Jones

lennardJones=a/r^12 - b/r^6;
force=-D[lennardJones,r]; (*La force est la derivée (en 1D) du potentiel*)

Le potentiel de Lennard-Jones normalisé

(*La résolution du système à deux équations permet de trouver l'expression de a et b en fonction de rmin et emin.*)

equationDistanceMinimum=(force/.r-> rMin)==0; (* /. est une règle, i.e. r prend la valeur rMin*)

equationEnergieMinimum=(lennardJones/.r-> rMin)==-eMin;

criticalPointSolutions=Solve[{equationDistanceMinimum,equationEnergieMinimum},{a,b}];

lennardJonesCritical=Expand[lennardJones/.criticalPointSolutions];
lennardJonesScaled=Expand[lennardJonesCritical/eMin];

Fonctions globale

Fonctions d’energie et de la force

energy[ρ_]=lennardJonesScaled/.r->ρ rMin (*la fonction énergie est l'equivalent normalisé de l'expression initiale du potentiel de Lennard-Jones. L'énergie est exprimée par rapport à l'énergie minimale et ρ la distance par rapport à la distance minimale.*)
forceScaled[ρ_]=D[energy[ρ],ρ];

Le potentiel de Lennard-Jones, ainsi que la force, peuvent être séparé en deux parties : attractive et répulsive.

energyAttractive[ρ_]=lennardJonesScaled[[1,1]]/.r->ρ rMin;
energyRepulsive[ρ_]=lennardJonesScaled[[1,2]]/.r->ρ rMin;
forceScaledAttractive[ρ_]=(forceScaled[ρ][[1,2]])/.r->ρ rMin;
forceScaledRepulsive[ρ_]=(forceScaled[ρ][[1,1]])/.r->ρ rMin;

Algorithme de dynamique moléculaire

velocityVerlet[{positions_, oldVelocities_,oldForces_, oldEnergies_}, dt_] :=Module[{deltaPositions,newpositions,newVelocities,newForces, newEnergies},
deltaPositions = (oldVelocities dt + 0.5 oldForces dt^2);
newpositions = (positions + deltaPositions);
{newForces,newEnergies}=particleForcesAndEnergies[newpositions];
newVelocities =  (oldVelocities + 0.5(oldForces + newForces)dt);
{ newpositions, newVelocities,newForces,newEnergies}
]

Démonstration 1

Fonctions

Graphiques

baseGraphicEnergy[ρValue_]:=Module[
{x=ρValue,y=energy[ρValue][[1]]},Plot[energy[ρ],{ρ,0.6,1.8},
(* apparence\[Implies] *)Frame->True,FrameStyle->Directive[Black,Thick],BaseStyle->16,ImageSize->500,PlotStyle->{Red,Thick},PlotRange->{{0.75,1.5},{-3,6}}, FrameLabel->{"Distance Interatomique normalisée, ρ = r/(Subscript[r, min])","Energie Potentielle normalisée, Φ = E/(Subscript[Ε, min])"},PlotLabel->"Potentiel de Lennard-Jones",Epilog->{{Gray,Dashed,Thick,Line[{{0.75,y},{x,y}}],Line[{{x,y},{x,-3}}]},
{Red,PointSize[0.025],Point[{1,-1}]},{Gray,PointSize[0.025],Point[{x,y}]},Inset["Minimum\n du \n potentiel \n ↓",{0.996,0.7}]}]]

baseGraphicAttraction=Plot[energyAttractive[ρ],{ρ,0.7,1.5},PlotStyle->{Blue,Thick}];
baseGraphicRepulsion=Plot[energyRepulsive[ρ],{ρ,0.7,1.5},PlotStyle->{Green,Thick}];

Force

baseGraphicForce[ρValue_]:=Module[
{x=ρValue,y=forceScaled[ρValue][[1]]},Plot[forceScaled[ρ],{ρ,0.8,1.5},Frame->True,FrameStyle->Directive[Black,Thick],BaseStyle->16,ImageSize->500,PlotStyle->{Red,Thick},PlotLabel->"Force exercée, F",PlotRange->{{0.75,1.5},{-30,10}},FrameLabel->{"Distance Interatomique normalisee, ρ = r/(Subscript[r, min])","Force, F"},Epilog->{{Gray,Dashed,Thick,Line[{{0.75,y},{x,y}}],Line[{{x,y},{x,-40}}]},{Red,PointSize[0.025],Point[{1,0}]},{Gray,PointSize[0.025],Point[{x,y}]},Inset["Force\n minimum \n ↓",{0.997,5.5}]}]]

baseGraphicAttractionForce=Plot[forceScaledAttractive[ρ],{ρ,0.7,1.5},PlotStyle->{Blue,Thick}];
baseGraphicRepulsionForce=Plot[forceScaledRepulsive[ρ],{ρ,0.7,1.5},PlotStyle->{Green,Thick}];

Illustration du potentiel entre 2 atomes

pairOfAtoms[ρValue_]:=Module[{c1={-ρValue/2,0},c2={ρValue/2,0}},link=Plot[TriangleWave[{-0.2,0.2},x],{x,-ρValue/2,ρValue/2},Axes-> False];Show[Graphics[{
{Transparent,Rectangle[{-2,-2},{2,2}]},
Disk[c1,0.15],Disk[c2,0.15],(*atomes*)
(*{Line[{c1,c2}]},(*Lien entre les atomes*)*)
{Line[{c1,c1-{0,0.5}}]},
{Line[{c2,c2-{0,0.5}}]},
{Dashed,Line[{{-0.5,0},{-0.5,-0.5}}]},
{Dashed,Line[{{0.5,0},{0.5,-0.5}}]},
{Dashed,{Arrowheads[{-0.04,.04}],Arrow[{{-0.5,-0.5},{+0.5,-0.5}}]}},(*Distance minimale*)
Style[Text["r\!\(\*SubscriptBox[\(\\ \), \(min\)]\)",{0,-0.7}],16],
{If[forceScaled[ρValue][[1]]≥ 0,Blue,Green],Thickness[0.009],
(*La longueur des fleches est proportionelle a la force exercee entre les atomes*)Arrow[{c1+{0,0.5},c1+{forceScaled[ρValue][[1]]/5,0.5}}],Arrow[{c2+{0,0.5},c2-{forceScaled[ρValue][[1]]/5,-0.5}}]}
},ImageSize->350,AspectRatio->1],link,Graphics[{Disk[c1,0.15],Disk[c2,0.15]}](*atomes*)]]

background=Graphics[Transparent,ImageSize->650];

Visualisation

visualisation1 =
Manipulate[
Grid[{{Show[
baseGraphicEnergy[ρ],If[flagAttraction,baseGraphicAttraction,background],If[flagRepulsion,baseGraphicRepulsion,background]],
Show[baseGraphicForce[ρ],If[flagAttractionForce,baseGraphicAttractionForce,background],If[flagRepulsionForce,baseGraphicRepulsionForce,background]]}}],Style["Initialisation",Bold],
Button["Energie minimum",ρ=1,ImageSize->{300,50},ImageMargins->{{80,5},{1,2}}],
Style["Paramètres",Bold],{{ρ,1,"Distance"},0.95,1.3},
Style["Illustration",Bold],
Style["Le potentiel décrit les interactions entre paire d'atomes.\nNB: La force exercée est positive en traction et négative en compression. \nLa force illustrée sur le graphique est la force exercée,\nles flèches bleu indiquent la force de réaction."(*legende modifiable*),16],
Dynamic[Show[pairOfAtoms[ρ],ImageSize-> 400]],
Delimiter,Style["Affichage",Bold],
{{flagAttraction,False,"Energie d'attraction"},{False,True}},{{flagRepulsion,False,"Energie de répulsion"},{False,True}},
{{flagAttractionForce,False,"Force attractive"},{False,True}},{{flagRepulsionForce,False,"Force répulsive"},{False,True}},
Delimiter,LabelStyle->Directive[FontSize-> 16],ControlPlacement->Left];

Démonstration 2: Approximations et oscillations

Fonctions

energyAndForce[ε0_,r0_,color_]:=Module[{},
Show[Plot[ ε0(((((r0/r)))^12) -2((((r0/r)))^6))(*Expression du potentiel de Lennard-Jones avec r0= distance minimale et ε0 l'energie minimale*),{r,0.5,150},
(*Apparance\[Implies]*)
Frame->True,FrameStyle->Directive[Black,Thick],BaseStyle->16,ImageSize->650,PlotStyle->{color,Thick},PlotLabel->"",PlotRange->{{r0/2,2*r0},{-2*ε0,2*ε0}},FrameLabel->{"Distance Interatomique r [Å]","Energy potentielle E [kcal/mol]"}]]]

Approximations

Le developpement en series entiere, i.e. expansion de Taylor, permet d' approximer la courbe avec une parabole ou un polynome cubique.

parabolaApproximation = Normal@Series[energy[ρ],{ρ,1,2}][[1]]
cubicApproximation =Normal@ Series[energy[ρ],{ρ,1,3}][[1]];

Oscillations for the approximation

fEqualsMA= mu D[ρ[t],{t,2}] == (12/((ρ[t])^13))-(12/((ρ[t])^7));

mu=1;charateristicPeriod = 2 Pi/Sqrt[72];

We will construct a numerical solution for a starting point near equilibrium and compute it for two of our estimated periods.  If we are doing everything right so far, our numerical solution should show al
most exactly two oscillations:

numSol =NDSolve[{fEqualsMA,ρ'[0]==0,ρ[0]==1.1},ρ,{t,0,10 charateristicPeriod}];

numericalSolutionFromVariableStartingPoint[initialValue_]:=  ρ/.NDSolve[{fEqualsMA,ρ'[0]==0,ρ[0]==initialValue},ρ,{t,0,40}][[1]]

Animation

visualisation2=Manipulate[
animationExample = numericalSolutionFromVariableStartingPoint[startingPoint];
tPeriod = t/.FindRoot[animationExample[t]==1.5,{t,6,3,9}]//Quiet;
image=
Show[
Plot[{energy[ρ],
If[flagParabola,parabolaApproximation,background],
If[flagCubic,cubicApproximation,background]},
{ρ,0.5,1.5},
Epilog-> {PointSize[0.03],Black, Point[{animationExample[t],(energy[ρρ]/.ρρ->animationExample[t])[[1]]}]}(*Oscillation animation*),
PlotRange->{{0.75,1.5},{-1.25,2}},
(*Apparence*)
Frame-> True,FrameLabel->{"Distance Interatomique normalisee, ρ = r/(Subscript[r, min])","Energy Potentielle normalisee, Φ = E/(Subscript[Ε, min])"},
ImageSize-> Large,  BaseStyle-> {FontSize-> 16}, PlotStyle->{{Thick,Red},{Dashed,Thick,Blue},{Dashed,Thick,Purple}}, PlotLabel->"Approximation du potentiel de LJ autour du minimum",PlotLegends->{Style["Lennard-Jones",16],Style["Approximation parabolique",16],Style["Approximation cubique",16]}]],
Style["Approximations",16,Bold],
(*Controles*)
{{flagParabola,False,"Approximation parabolique"},{False,True}},{{flagCubic,False,"Approximation cubique"},{False,True}},
(*Zoom*)
Style["Zoom",16],
Dynamic[
Plot[{energy[ρ],
If[flagParabola,parabolaApproximation,background],
If[flagCubic,cubicApproximation,background]},
{ρ,0.5,1.5}, PlotStyle->{{Thick,Red},{Dashed,Thick,Blue},{Dashed,Thick,Purple}},Frame-> True,
PlotRange->{{0.925,1.075},{-1.25,0}}]],
Delimiter,
(*Controles*)
Style["Oscillations",16,Bold],
{{startingPoint,0.98,"Position Initiale"},0.98,1.2},{t,0,tPeriod },
Delimiter,
(*Question*)
Style["Questions",16,Bold],
Style["Pouvez-vous deduire une consequence de l'assymetrie du potentiel de LJ? "],
ControlPlacement->Left,LabelStyle->Directive[FontSize-> 16]
];

Exercice 1 : Matériaux et énergies de liaisons

exercise1=
Manipulate[
(*If[ε1>=ε2,*)Show[energyAndForce[ε1,r01,Blue],energyAndForce[ε2,r02,Red],PlotRange->{{Min[r01,r02]/2,2*Max[r01,r02]},{-2*Max[ε1,ε2],2*Max[ε1,ε2]}}](*,Show[energyAndForce[ε2,r02,Red],energyAndForce[ε1,r01,Blue]]]*),(*La fonction energyAndForce trace la courbe de Lennard-Jones en fonction de parametres ε et r.*)
Dynamic@"Exercise 1",Style["Illustre different types de liaisons chimiques/ ou materiau sur le graphique"],Delimiter,
Style["Constantes: \n Le potentiel de LJ peut être écrit sous plusieurs formes:\n E = (a/(r^12))-(b/(r^6))  ou E = (Subscript[ε, min]) [(((Subscript[r, min])/r))^12 -2(((Subscript[r, min])/r))^6].\n
Avec les valeurs empiriques des constantes (Subscript[ε, min]) en [kJ (mol^(-1))] et (Subscript[r, min]) en [Å].\nLes constantes a, b et (Subscript[ε, min]), (Subscript[r, min]) sont reliées par: a=4(Subscript[ε, min])((Subscript[r, min])^12) et b=4(Subscript[ε, min])((Subscript[r, min])^6)."],
Style["Matériau 1",Bold, Blue],{ε1,1},{r01,1},Delimiter,Style["Matériau 2",Bold,Red],{ε2,1},{r02,1},Delimiter,
Style["Questions"],
Style["Choisis deux matériaux et compare leur courbe de potentiel\n 
Substance | (Subscript[ε, 0]) | (Subscript[r, 0])
Xe | 1.77 | 4.10  \n(Subscript[H, 2])O | 15.32 | 3.02  \n(Subscript[N, 2]) | 21.79 | 3.92    \nZn | 62.40 | 2.46\n(Subscript[CO, 2]) | 41.27 | 4.49    \n(Subscript[TiCl, 4]) | 95.71 | 6.72 \n
Quelle variable est liée a l'énergie de liaison?"],
LabelStyle->Directive[FontSize-> 16],ControlPlacement->Left];

(*CreateDocument[exercise1,WindowSize-> All]*)

(*Xe --> Van der Waals: <10 kJ/mol  
Hydrogen: 40 kJ/mol et 2.5 Å. 
Metallique: 25-200 kJ/mol 
Covalente:125-300 
Ionique:150-370*)

Exercice 2 : Force et potentiel

lennardJonesFormula=StringJoin["Potentiel de Lennard-Jones : \n \n     ϕ = (Subscript[e, min])[(((Subscript[r, min])/r))^12 - 2 (((Subscript[r, min])/r))^6] \n "];
myForce[myForceExpression_,ε0_,r00_]:=Module[{},
Show[Plot[myForceExpression/.{emin->ε0,rmin->r00},{r,1,5},Frame->True,PlotRange->{{1,5},{-2,2}},FrameStyle->Directive[Black,Thick],BaseStyle->16,ImageSize->650,PlotStyle->{Black,Thick},PlotLabel->"",FrameLabel->{"Distance Interatomique r [Å]","Force [10^-9 N]"}]]]

answer=Plot[-e0*1.6(-((12 (r0^12))/(r^13))+((12 (r0^6))/(r^7)))/.{e0->1,r0->3},{r,1,5},PlotRange->{{1,5},{-2,2}},PlotStyle->{Thick,Red,Dashed},PlotLegends-> Placed["\n F=-(Subscript[e, min])(-\!\(\*FractionBox[\(12\\ \*SuperscriptBox[SubscriptBox[\(r\), \(min\)], \(12\)]\), SuperscriptBox[\(r\), \(13\)]]\)+\!\(\*FractionBox[\(12\\ \*SuperscriptBox[SubscriptBox[\(r\), \(min\)], \(6\)]\), SuperscriptBox[\(r\), \(7\)]]\))",{Right,Top}]];

exercise2=
Manipulate[Module[{emin=1,rmin=3},
Show[
myForce[Force(*/.{e0->0.066,r0->3}*),emin,rmin],
If[flagAnswer,answer,background]]],
Style["Exercise 2",Bold],
Style["Calculez à la main l'expression de la force de réaction. Entrez la formule ainsi obtenue dans le champ suivant et comparez avec la réponse"],Delimiter,Dynamic@lennardJonesFormula,{Force},Dynamic@"NB: Gardez les mêmes noms de variable: r, rmin et emin.",Delimiter,{{flagAnswer,False,"Voir la réponse"},{False,True}},Delimiter,
Style["Commentaire:\nF>0 en traction et F<0 en compression.\nLa force de réaction est celle exprimée comme F=-(dE/dr).\nExemple: si la distance r entre les atomes est plus petite que (Subscript[r, min]),\nla force externe qui les rapproche est une force de compression\n(F<0) et la force de réaction qui tend a les repousser à leur\nposition d'équilibre est une force de traction (F>0)."],
LabelStyle->Directive[FontSize-> 16],ControlPlacement->Left];

(*CreateDocument[exercise2,WindowSize->All]*)

Exercice 3 : Énergie d’un système de particules

particleForcesAndEnergiesCompiled= Compile[{{positions,_Real,2}}, 
Block[{distanceMatSq,distanceMatSixth, differenceMat,ret},
differenceMat = Outer
[Subtract,N[positions],N[positions],1] ;
distanceMatSq =Map[#.#&,differenceMat,{2}] + IdentityMatrix[Length[positions]];
distanceMatSixth= distanceMatSq^3;
Flatten[{(* Here it is where LJ comes into place -> *)
(*Forces*)Total[-12 differenceMat ( 1-distanceMatSixth)/distanceMatSq^7],
(*Energies*)Total[(1-2distanceMatSixth)/(distanceMatSixth^2)]
}]
]
];

particleForcesAndEnergies[positions_]:= Module[{res},
res = particleForcesAndEnergiesCompiled[positions];
{Partition[Take[res,2 Length[positions]],2],
Take[res,- Length[positions]]}
]

Initial hexagonal lattice

hexagonalLatticeVectors =N@ {{1,0},{1/2,Sqrt[3]/2}}
pts =(#. hexagonalLatticeVectors  )&/@Tuples[Range[-45,45],2];
ptsGraphic =Graphics[Point[pts], Axes->True, ImageSize->Medium];

Pruning to a block (for subsequent pruning)

pts = Cases[pts,{x_,y_}/;(-20 ≤  x ≤  20 && -10 < y < 25)]; 
pts = Cases[pts,{x_,y_}/;(-5 ≤  x ≤  5 && -5 < y < 5)];
ptsGraphic =Graphics[Point[pts], Axes->True, ImageSize->Large];

ptsSquare=(#. {{1,0},{0,1}}  )&/@Tuples[Range[-4.5,4.5],2];

{forces,energies}=particleForcesAndEnergies[pts];

m=7;

circularParticle[ radius_]:= Select[Tuples[Range[-Round@radius,Round@radius],2],(Norm[#]≤   radius)&]

displayParticles[positions_, energies_,m_]:=Module[{},
Graphics[{
  Transpose[
   {
ColorData["Rainbow"][Sqrt@Clip[Rescale[#, {Min@energies, Max@energies}],{0,1}]] & /@ energies,
    Disk[#, 1/2] & /@ positions}
   ]},PlotRange->m,ImageSize->400]]

colorbar[energies_]:=Module[{allrange,allRangeColor,allcommand,allText,alldisk},
allrange=Table[i,{i,Min@energies, Max@energies,0.55}];
allRangeColor=ColorData["Rainbow"][Sqrt@Clip[Rescale[#, {Min@energies, Max@energies}],{0,1}]] & /@ allrange;
allcommand=Partition[Riffle[
   allRangeColor,  Disk[{0,#+.1}, 0.5] & /@allrange],2];
allText=Graphics[{White,Text[
Style[NumberForm[#,{2,1}],Medium]
,{0,#}]}]& /@ allrange;
alldisk=Thread[Graphics[allcommand]];
Show[alldisk,allText]
]

exercise3=Manipulate[


If[flagPoints,points=pts,points=ptsSquare];
If[flagPointsCircular,points=circularParticle[radius]];

voisin=Nearest[points,p][[1]];
n=Position[points,voisin][[1]];
{forces,energies}=particleForcesAndEnergies[points];
totalEnergy=Total@energies;
numberOfAtoms=Length@energies;
displayParticles[points,energies,7],Style["Exercice 3\n",16],
Dynamic[Style[StringJoin["Energie = ",ToString[NumberForm[energies[[n]][[1]],3]]],Bold,16]],
Dynamic@If[True,Style[StringJoin["Energie totale : ", ToString@totalEnergy,"   Nombres d'atomes : ",ToString@numberOfAtoms],16],""],
(*{{flagPointsTotalInfos,True,Style["Details",16]},{True,False}},*)
(*Style[" Utilise le bouton Détails  pour visualiser l'énergie totale du systeme et le nombre d'atomes.",14],*)
Delimiter,
Style["Utilise ton curseur pour afficher l'énergie de chaque atome."],
Style["Question",Italic,16],
Style["Comment l'énergie par atome varie en fonction de la position de l'atome?",16],
Delimiter,
Style["Utilise les options suivantes pour modifier la configuration du système"],
{{flagPoints,False,Style["Hexagonal",16]},{True,False}},
{{flagPointsCircular,False,Style["Particule",16]},{True,False}}
,{{radius,1,Style["Rayon",16]},1,5,1},

Style["Questions",16,Italic],
Style["Comment l'énergie par atome varie en fonction de la taille de la particule?",16],
Style["Comment l'énergie varie en fonction de la configuration?"],
{{p,{0,0}},Locator},
LabelStyle-> 16,ControlPlacement->Left];

(*CreateDocument[exercise3,WindowSize-> All]*)

Démonstration 3: Énergie potentielle et cinétique avec dynamique moléculaire

Module[{
 collision={{0,-6}},
initialVelocity={{0,10}},
dt=.01,
iter=0(*forcesInitiales,energiesInitiales,velocitiesInitiales,allParticlesInitiales*)},
pointsColl=(#. {{1,0},{0,1}}  )&/@Tuples[Range[-4.5,4.5],2];
pointsCollHexTmp=(#. (N@ {{1,0},{1/2,Sqrt[3]/2}})  )&/@Tuples[Range[-45,45],2];
pointsCollHex = Cases[pointsCollHexTmp,{x_,y_}/;(-3 ≤  x ≤  3 && -4< y < 2)];
(*Dynamic@If[hexagonalFlag,tousPoints=pointsCollHex,tousPoints=pointsColl];*)
allParticles=Catenate[{pointsCollHex,collision}];
allParticlesInitiales=allParticles;
{forces,energies}=particleForcesAndEnergies[allParticles];
{forcesInitiales,energiesInitiales}=particleForcesAndEnergies[allParticles];
velocitiesParticles=ConstantArray[{0,0},Length@pointsCollHex];
velocitiesInitiales=ConstantArray[{0,0},Length@pointsCollHex];
velocities=Catenate[{velocitiesParticles,initialVelocity}]];

Pour éviter que l’utilisateur doive executer la simulation en dynamique, l’animation est enregistree ici.

Clear@listFrames
listFrames={};
Do[
listFrames=AppendTo[listFrames,Grid[{{displayParticles[allParticles,energies,{{-10,10},{-10,7}}],colorbar[energiesInitiales]}}]];
{allParticles,velocities,forces,energies}=velocityVerlet[{allParticles,velocities,forces,energies},0.01];
,55]

visualisation3=Manipulate[
Dynamic@simulation,
(*Controls*)
Style["Exercise 4",16],
Style["Cette visualisation est une simulation de collision.\nLorsque la particule collisione le système,\nun réarrangement se produit.",16],
Delimiter,
(*{{dt,0.01,Intervalle},0.01,0.10},
{{iter,5,Iterations},1,10,1},*)
Button["Simulation",simulation=ListAnimate[listFrames,DefaultDuration->10,AnimationRepetitions->1,AppearanceElements->None],Method-> "Preemptive",ImageSize->{250,50}],
Button["Réinitialisation",simulation=ListAnimate[{Grid[{{displayParticles[allParticlesInitiales,energiesInitiales,{{-10,10},{-10,7}}],colorbar[energiesInitiales]}}]}],Method-> "Preemptive",ImageSize->{250,50}],
Style["Question",Italic,16],
Style["Quelle autre énergie entre en jeu dans cette simulation?",16]
,ControlPlacement->Left,LabelStyle->16];

(*CreateDocument[visualisation3,WindowSize->All]*)

Exercice 4 : Calcul du module de Young et de la contrainte, E et σ

Clear@r0;Clear@ε0;

(*Starting from the Lennard-Jones potential, the relationship can be calculated as the variation of the force around the equilibrium position, i.e. the seconde derivative of lennard-jones potential. It represent the curvature of the well.*)

lennardJonesWithUnits=ε0 (((((r0/r)))^12)-2 ((((r0/r)))^6));

forceWithUnits=D[lennardJonesWithUnits,r];

curvatureWithUnits=D[lennardJonesWithUnits,{r,2}];

modulusAroundMinimalDistance=curvatureWithUnits/.r->r0;

(*Starting from the parabolic approximation of the potential, the force is then a linear relationship between the displacement, i.e. r-(Subscript[r, 0 ])and the force. 
F = k * (r-(Subscript[r, 0])) *)

approxLennardJonesWithUnits=Normal@Series[ε0 (((((r0/r)))^12)-2 ((((r0/r)))^6)),{r,r0,2}];

approxForceWithUnits=D[approxLennardJonesWithUnits,r];

approxCurvatureWithUnits=D[approxLennardJonesWithUnits,{r,2}];

pairWithSpring[ρValue_,εValue_]:=Module[{c1={-ρValue/2,0},c2={ρValue/2,0}},
link=Plot[TriangleWave[{-0.2,0.2},x],{x,-ρValue/2,ρValue/2},Axes-> False,PlotStyle->{Blue,Thickness[εValue/(ρValue^2*60)]}];
Show[Graphics[{
{Transparent,Rectangle[{-2,-2},{2,2}]},
(*{Line[{c1,c2}]},(*Lien entre les atomes*)*)
{Line[{c1,c1-{0,0.5}}]},
{Line[{c2,c2-{0,0.5}}]},
{Dashed,{Arrowheads[{-0.04,.04}],Arrow[{c1-{0,0.5},c2-{0,0.5}}]}},(*Distance minimale*)
Style[Text["(Subscript[r, 0])",{0,-0.7}],16]
},ImageSize->350,AspectRatio->1],link,Graphics[{Disk[c1,0.25],Disk[c2,0.25]}](*atomes*)]]

exercise4=Manipulate[
maxForce=FindMaximum[forceWithUnits/.{ε0->εChosen,r0-> r0Chosen},r]//Quiet;
maxsigma=maxForce[[1]]*0.160/( r0Chosen^2 );
Show[
Plot[{
(*approxLennardJonesWithUnits/.{ε0->εChosen,r0-> r0Chosen},*)approxForceWithUnits/.{ε0->εChosen,r0-> r0Chosen},lennardJonesWithUnits/.{ε0->εChosen,r0-> r0Chosen},forceWithUnits/.{ε0->εChosen,r0-> r0Chosen}},
{r,0.5`,10},
Frame->True,FrameStyle->Directive[Black,Thick],BaseStyle->16,ImageSize->650,PlotStyle->{(*{Dashed,Blue},*){Thickness[0.007],Blue},{Dashed,Red},{Red}},PlotLabel->"",PlotLegends->{(*"Approximation parabolique du potentiel LJ",*)"Force du potential approximé","Potentiel de LJ","Force LJ"},PlotRange->{{1/2,2.5 },{-100,300 }},FrameLabel->{"Distance Interatomique r [Å]","Energie potentielle E [kJ/mol] \n Force F [kN/mol]"},Epilog-> Inset[Style[StringJoin["E=\!\(\*FractionBox[\(\(\\ \)\(72 \*SubscriptBox[\(ε\), \(0\)]\)\), SuperscriptBox[SubscriptBox[\(r\), \(0\)], \(3\)]]\)=",ToString[NumberForm[72 εChosen 0.160/ (r0Chosen^3),{5,2}]]," TPa\n (SuperStar[(Subscript[σ, max])])=((Subscript[F, max])/((Subscript[r, 0])^2))=",ToString[NumberForm[maxsigma,{5,2}]]," TPa"],24],
Scaled[{0.75,0.75}]] ],
Graphics[{Dashed,Thickness[0.005],GrayLevel[0.2],Line[{{0,maxForce[[1]]},{maxForce[[2,1,2]],maxForce[[1]]}}]}]],
(*Controls*)
Style["Exercise 4"],
Style["La valeur du module de Young et de la contrainte\nmaximale théorique dépendent de la dérivée du \npotentiel à la position d'équilibre."],
{{εChosen,50,"(Subscript[ε, 0])"},25,75},{{r0Chosen,1,"(Subscript[r, 0])"},0.7,1.3},
Delimiter,Dynamic[pairWithSpring[r0Chosen,εChosen/200]],ControlPlacement->Left,LabelStyle->16];

(*CreateDocument[exercise4,WindowSize->All]*)

Visualisation finale

CreateDocument[
Manipulate[
Which[
Select==1,visualisation1,
Select==2,exercise1,
Select==3,exercise2,
Select==4,exercise3,
Select==5,visualisation3,
Select==6,visualisation2,
Select==7,exercise4
],
{{Select,1},
{1-> "Potentiel de Lennard-Jones: Force et énergie",
2 -> "Energie de liaisons pour différents matériaux",
3 -> "Expression mathématique de la force",
4 -> "Energie totale et energie de surface",
5 -> "Energie potentielle et énergie cinétique",
6 -> "Approximation parabolique et cubique du potentiel LJ",
7-> "Calcul du module de Young et de la contrainte, E et σ"}}]
(*,WindowSize->Full*)
(*Commenter le (*,WindowSize->Full*) pour faire apparaître la fenêtre dans le même bureau*)]

Références