Energy Minimization

DEFINITIONS

s12=Sqrt[s1 s2]; e12=Sqrt[e1 e2];
v11=Compile[{{r,_Real}},4e1((s1/r)^6-1)(s1/r)^6];
dv11=Compile[{{r,_Real}},-24e1(2(s1/r)^6-1)(s1/r)^7];
v22=Compile[{{r,_Real}},4e2((s2/r)^6-1)(s2/r)^6];
dv22=Compile[{{r,_Real}},-24e2(2(s2/r)^6-1)(s2/r)^7];
v12=Compile[{{r,_Real}},4e12((s12/r)^6-1)(s12/r)^6];
dv12=Compile[{{r,_Real}},-24e12(2(s12/r)^6-1)(s12/r)^7];

Remove[e1,s1]

InputForm[Simplify[D[4((1/r)^6-1)(1/r)^6,r]]]

InputForm[Simplify[D[4((1/r)^8-1)(1/r)^8,r]]]

Remove[eforce];
Module[{na,nl,d,dxy,v11,dv11,v22,dv22,v12,dv12,s12,e12,ei,fij,e,f},
eforce[q_,l_,{s1_,e1_,s2_,e2_},neigh_:Automatic,pbc_:None]:=(
na=Length[q];
s12=Sqrt[s1 s2]; e12=Sqrt[e1 e2];
v11=Compile[{{r,_Real}},4e1((s1/r)^8-1)(s1/r)^8];
dv11=Compile[{{r,_Real}},-32e1(2(s1/r)^8-1)(s1/r)^9];
v22=Compile[{{r,_Real}},4e2((s2/r)^8-1)(s2/r)^8];
dv22=Compile[{{r,_Real}},-32e2(2(s2/r)^8-1)(s2/r)^9];
v12=Compile[{{r,_Real}},4e12((s12/r)^8-1)(s12/r)^8];
dv12=Compile[{{r,_Real}},-32e12(2(s12/r)^8-1)(s12/r)^9];
If[neigh===Automatic,nl=DeleteCases[Table[If[i==j,{},j],{i,1,na},{j,1,na}],{}];,
nl=neigh];
e=0; f=ConstantArray[{0,0},na];
Do[
If[i<nl[[i,j]],
dxy=q[[i]]-q[[nl[[i,j]]]]; 
If[pbc=!=None,
While[dxy[[1]]>pbc[[1]]/2,dxy[[1]]-=pbc[[1]]];
While[dxy[[1]]<-pbc[[1]]/2,dxy[[1]]+=pbc[[1]]];
While[dxy[[2]]>pbc[[2]]/2,dxy[[2]]-=pbc[[2]]];
While[dxy[[2]]<-pbc[[2]]/2,dxy[[2]]+=pbc[[2]]];
];
d=Sqrt[dxy.dxy];
If[l[[i]]==1 && l[[j]]==1,
ei=v11[d]; fij=-dxy/d*dv11[d];,
If[l[[i]]==2 && l[[j]]==2,
ei=v22[d]; fij=-dxy/d*dv22[d];,
ei=v12[d]; fij=-dxy/d*dv12[d];]
];
e+=ei;
f[[i]]+=fij;
f[[nl[[i,j]]]]-=fij;
],{i,1,na},{j,1,Length[nl[[i]]]}];
{e,f}
);
]

Tests

nx=24;ny=2;
Module[{},
BuildLattice[zz_,epitaxy_:False]:=(
e1=e2=1;s1=1;s2=s1*zz; z0=1.1116;z0=1.08865;
a1=s1 z0;a2=s2  z0;
b1=1.{1,0};b2={0.5,Sqrt[3]/2.};
 wx=a1*nx;
q1=Flatten[Table[a1(i b1+j b2 ),{i,1,nx},{j,-2,0}],1];
If[epitaxy,
q2=Flatten[Table[a1(i b1+j b2),{i,1,nx},{j,1,ny}],1];,
q2=Flatten[Table[a2(i b1+j b2),{i,1,nx*a1/a2},{j,1,ny}],1];
];
q2=(#1+{RandomReal[],RandomReal[]}*0.02)&/@q2;
myq=Join[q1,q2];
mynl=Join[Table[Length[q1]+j,{i,1,Length[q1]},{j,1,Length[q2]}],
Table[i,{j,1,Length[q2]},{i,1,Length[q1]+Length[q2]}]];
myl=Join[ConstantArray[1,Length[q1]],ConstantArray[2,Length[q2]]];
Do[While[myq[[i,1]]>wx/2,myq[[i,1]]-=wx];While[myq[[i,1]]<-wx/2,myq[[i,1]]+=wx],{i,1,Length[myq]}]
);
]

nx=12;ny=8

BuildLattice[1]

{mye,myf}=eforce[myq,myl,{s1,e1,s2,e2}];

myf=myq*0;

dx=0.1;myf=myq*0;
Do[
of=myf;
{mye,myf}=eforce[myq,myl,{s1,e1,s2,e2}];
If[Flatten[of].Flatten[myf]<0,dx*=0.8];
myq+=myf/Norm[myf]*dx;,{i,1,150}]

Dynamic[{dx,mye/Length[q2],Total[Flatten[myf]^2]}]

Dynamic[Graphics[{Arrowheads[0.01],Table[{If[myl[[i]]==1,Red,Blue],Disk[myq[[i]],0.25],Black,Arrow[{myq[[i]],myq[[i]]+myf[[i]]*5}]} ,{i,1,Length[myl]}]},PlotRange->All,ImageSize->Large]]

Graphics[{Point/@myq,Disk[myq[[110]],0.2],Disk[myq[[111]],0.2]}]

Norm[myq[[110]]-myq[[111]]]

2.^(1/6)

D[v[r],r]

Plot[{v[r],v[r] *10/(9+Exp[r])},{r,1,3},PlotRange->All]

lq=Flatten[Table[Random[]*0.0001+i*{1/2,Sqrt[3]/2}+j*{1,0},{i,0,10},{j,0,20}],1];

lq=1.Flatten[Table[Random[]*0.0001+i*{1/2,Sqrt[3]/2}+j*{1,0},{i,0,10},{j,0,20}],1];
lq=mkrot[lq,Pi*0.1];
lq=DeleteCases[({0,-1}+#1)&/@lq,{_,y_}/;y<0];
lq=Join[lq,({1,-1}+{1,-1}*#1)&/@lq];

lq={{0,0},{1,1},{1,-1}}

neigh=mknl[lq,4];
 dq=0.0001;
ffix=If[Length[#1]<3,1,1]&/@neigh;

ef=mkrigid[lq,eforce[lq,{v,v1},neigh]];
Do[
oq=lq;oe=ef[[1]];
lq+=(ffix*ef[[2]])*dq;
ef=mkrigid[lq,eforce[lq,{v,v1},neigh]];
If[ef[[1]]>oe, ef[[1]]=oe; lq=oq; dq*=0.5,dq*=1.2];
,{k,1,100}];

ef=mkrigid[lq,eforce[lq,{v,v1},neigh]];
g=h=ef[[2]];
Do[
oq=lq;oh=h;og=g;of=ef;
lq+=dq*h;
ef=mkrigid[lq,eforce[lq,{v,v1},neigh]];
If[ef[[1]]>of[[1]],ef=of;lq=oq; dq*=0.5,
dq*=1.2;
g=ef[[2]]; b=Flatten[g].(Flatten[g]-Flatten[og])/Flatten[og].Flatten[og];h=g+Max[0,b] h;
];
,{k,1,50}];

ef=eforce[lq,{v,v1},neigh];
g=h=ef[[2]];
Do[
oq=lq;oh=h;og=g;of=ef;
lq+=dq*h;
ef=eforce[lq,{v,v1},neigh];
If[ef[[1]]>of[[1]],ef=of;lq=oq; dq*=0.5; h=g,
dq*=1.2;
g=ef[[2]]; b=Flatten[g].(Flatten[g]-Flatten[og])/Flatten[og].Flatten[og];h=g+Max[0,b] h;
];
,{k,1,50}];

dq=0.001

lq+=(ef[[2]])*dq;

ef=eforce[lq,{v,v1},neigh];

Dynamic[{ef[[1]],Flatten[ef[[2]]].Flatten[ef[[2]]],dq,b}]
Dynamic[Show[Graphics[Table[{If[ffix[[i]]==0,Red,Black],Red,Circle[lq[[i]],0.5],Black,Arrowheads[0.03],Arrow[{lq[[i]],lq[[i]]+10*dq*ef[[2,i]]}]},{i,1,Length[lq]}]],Frame->True,AspectRatio->Automatic]]

ef=mkrigid[lq,eforce[lq,{v,v1},neigh]];

Graphics[Table[{If[ffix[[i]]==0,Red,Black],Point[lq[[i]]],Arrow[{lq[[i]],lq[[i]]+0.5ef[[2,i]]}]},{i,1,Length[lq]}],Frame->True]

Interactive

nx=24;ny=2;
Module[{},
BuildLattice[zz_,epitaxy_:False,e2e1_:1]:=(
e1=1; e2=e2e1*e1; s1=1;s2=s1*zz; z0=1.1116;
a1=s1 z0;a2=s2  z0;
b1=1.{1,0};b2={0.5,Sqrt[3]/2.};
 wx=a1*nx;
q1=Flatten[Table[a1(i b1+j b2 ),{i,1,nx},{j,-2,0}],1];
If[epitaxy,
q2=Flatten[Table[a1(i b1+j b2),{i,1,nx},{j,1,ny}],1];,
q2=Flatten[Table[a2(i b1+j b2),{i,1,nx*a1/a2},{j,1,ny}],1];
];
q2=(#1+{RandomReal[],RandomReal[]}*0.02)&/@q2;
myq=Join[q1,q2];
mynl=Join[Table[Length[q1]+j,{i,1,Length[q1]},{j,1,Length[q2]}],
Table[i,{j,1,Length[q2]},{i,1,Length[q1]+Length[q2]}]];
myl=Join[ConstantArray[1,Length[q1]],ConstantArray[2,Length[q2]]];
Do[While[myq[[i,1]]>wx/2,myq[[i,1]]-=wx];While[myq[[i,1]]<-wx/2,myq[[i,1]]+=wx],{i,1,Length[myq]}]
);
] 
go=False;

Manipulate[
Dynamic[fig=Graphics[{Do[While[myq[[i,1]]+wx dvx>wx/2,myq[[i,1]]-=wx];While[myq[[i,1]]+wx dvx<-wx/2,myq[[i,1]]+=wx],{i,1,Length[myq]}];
Arrowheads[0.01],Table[{If[myl[[i]]==1,Red,Blue],EdgeForm[Black],Disk[myq[[i]]+{wx dvx,0},0.25]},{i,1,Length[myl]}],
Black,Dashed,Line[{{-wx/2,-2.5},{-wx/2,4.5}}],Line[{{wx/2,-2.5},{wx/2,4.5}}],
If[go,Text[Style["Residual Force: "<>ToString[SetAccuracy[Sqrt[Total[Flatten[myf]^2]],3]],FontSize->12,FontFamily->"Helvetica"],Offset[{10,5},{-wx/2,4.5}],{-1,-1}],{}],
Text[Style["Lattice mismatch δ: "<>ToString[SetAccuracy[zz,3]],Bold,FontSize->12,FontFamily->"Helvetica"],Offset[{0,0},{0,4.5}],{0,-1}],
Text[Style["Energy/atom: "<>ToString[SetAccuracy[mye/Length[q2],3]],Bold,FontSize->12,FontFamily->"Helvetica"],Offset[{-10,5},{wx/2,4.5}],{1,-1}]
},PlotRange->All,ImageSize->10*72]]
,
{{nl,1,"N. layers"},{1,2,3,4}},
{{ebea,1.0,"Relative bond energy (Subscript[ϵ, β])/(Subscript[ϵ, α])"},0.2,5},
Grid[{{
Control[{{zz,SetAccuracy[24/25-1.,4],"Lattice mismatch"},SetAccuracy[1./{1,25/24,26/24,28/24,36/24,23/24,22/24,20/24,18/24}-1,4]}],
Control[{{epi,False,"Epitaxial"},{True,False}}]
}}]
,
Grid[{{Button["Build layer",
ny=nl;BuildLattice[zz+1,epi,ebea];
{mye,myf}=eforce[myq,myl,{s1,e1,s2,e2},mynl,{wx,100}];
,ImageSize->80,Method->"Queued"],
Button["Relax",dx=0.1;go=True;
While[go && Sqrt[Total[Flatten[myf]^2]]>0.01,
of=myf;
{mye,myf}=eforce[myq,myl,{s1,e1,s2,e2},mynl,{wx,100}];
If[Flatten[of].Flatten[myf]<0,dx*=0.7,dx*=1.05];
Do[myf[[i]]*=0,{i,1,Length[q1]}];
myq+=myf/Norm[Flatten[myf]]*dx;
];go=False;Do[While[myq[[i,1]]+wx dvx>wx/2,myq[[i,1]]-=wx];While[myq[[i,1]]+wx dvx<-wx/2,myq[[i,1]]+=wx],{i,1,Length[myq]}];
,ImageSize->80,Method->"Queued"],
Button["Stop",go=False,ImageSize->80,Method->"Preemptive"]}}],
Control[{{dvx,0,"Visualization shift"},-0.5,0.5}]
,SaveDefinitions->True
]