File : lattice-ops.adb
with Messages, Lin.Ops, Lin.IO, Volumes, Lattice.IO;
use Messages, Lin.Ops, Lin.IO, Volumes, Lattice.IO;
package body Lattice.Ops is
procedure InitPatch(N1,N2,N3,LP: in Positive) is
B: BallPtr;
begin
if BPatch /= null then FreePatch; end if;
SetDim(N1,N2,N3);
LayerPeriod := LP;
if Verbose>0 then
if LP=2 then Show("initializing patch for HCP lattice"); end if;
if LP=3 then Show("initializing patch for FCC lattice"); end if;
end if;
if BPatch = null then
BPatch := new BallPatch(-Dim1 .. Dim1,-Dim2 .. Dim2,-Dim3 .. Dim3);
end if;
for I1 in -Dim1 .. Dim1 loop
for I2 in -Dim2 .. Dim2 loop
for I3 in -Dim3 .. Dim3 loop
if FitsType((I1,I2,I3)) then
B := new Ball;
B.ID := 0;
B.IP := (I1,I2,I3);
BPatch(I1,I2,I3) := B;
end if;
end loop;
end loop;
end loop;
for I1 in -Dim1 .. Dim1 loop
for I2 in -Dim2 .. Dim2 loop
for I3 in -Dim3 .. Dim3 loop
B := BallAt((I1,I2,I3));
if B /= null then
GetNbors(B.IP,B.Nb);
end if;
end loop;
end loop;
end loop;
Symmetries := Symm;
end InitPatch;
procedure FreePatch is
begin
if BPatch=null then return; end if;
for I1 in -Dim1 .. Dim1 loop
for I2 in -Dim2 .. Dim2 loop
for I3 in -Dim3 .. Dim3 loop
if BPatch(I1,I2,I3) /= null then Free(BPatch(I1,I2,I3)); end if;
end loop;
end loop;
end loop;
Free(Symmetries);
end FreePatch;
function BallAt(L: IPoint) return BallPtr is
begin
if not InPatch(L) then return null; end if;
return BPatch(L(1),L(2),L(3));
end BallAt;
procedure GetNbors(Q: in IPoint; N: out Nbors) is
F: constant IPoint := ( 1, 1,0);
G: constant IPoint := ( 2,-1,0);
H: constant IPoint := (-1, 2,0); -- H=F-G
P: IPoint;
begin
N(1) := BallAt(Q-F);
N(2) := BallAt(Q-G);
N(3) := BallAt(Q+H);
N(4) := BallAt(Q+F);
N(5) := BallAt(Q+G);
N(6) := BallAt(Q-H);
P := Q+Pt(1,0,-1);
if BallAt(P)=null then
P := Q+Pt(0,1,-1);
N(9) := BallAt(P-H);
else
N(9) := BallAt(P-G);
end if;
N(7) := BallAt(P);
N(8) := BallAt(P-F);
P := Q+Pt(1,0,1);
if BallAt(P)=null then
P := Q+Pt(0,1,1);
N(12) := BallAt(P-H);
else
N(12) := BallAt(P-G);
end if;
N(10) := BallAt(P);
N(11) := BallAt(P-F);
end GetNbors;
procedure AddMovable(L: in IPoint; B: in out Blob; OK: out Boolean) is
C: constant BallPtr := BallAt(L);
N,Nfix: Integer;
begin
CountNbors(C,N,Nfix);
OK := N=12;
if OK and then C.ID=0 then
B.M := B.M+Nfix; -- number of constraints
B.L := B.L+1;
B.V(B.L) := C;
C.ID := B.L;
end if;
end AddMovable;
procedure SetMovable(B: in PolyPtr; BM: in out Blob; OK: out Boolean) is
begin
ClearMovable(BM);
for N in 1 .. B.L loop
AddMovable(UnPack(B.V(N)),BM,OK);
exit when not OK;
end loop;
end SetMovable;
procedure ClearMovable(B: in out Blob) is
begin
for N in B.V'Range loop
if B.V(N) /= null then
B.V(N).ID := 0;
end if;
end loop;
B.M := 0;
B.L := 0;
end ClearMovable;
procedure ClearAllMovable is
B: BallPtr;
begin
for I1 in -Dim1 .. Dim1 loop
for I2 in -Dim2 .. Dim2 loop
for I3 in -Dim3 .. Dim3 loop
B := BallAt((I1,I2,I3));
if B /= null then B.Id := 0; end if;
end loop;
end loop;
end loop;
end ClearAllMovable;
procedure CheckMovable(M: in Positive; B: in PolyPtr) is
begin
if B.V'Last>99 then
Show("B.V'Last =",B.V'Last);
Show("Error: Poly too large");
end if;
for N in 1 .. B.L loop
if not Check(B.V(N)) then
Show("B.L =",B.L);
Show("Error: ball not in patch");
end if;
end loop;
if B.L<M then
Show("Error: unfilled Blob"); return;
end if;
for N in 1 .. M loop
if B.V(N)<Fac3 then
Show("Error: ball should be movable");
end if;
end loop;
for N in M+1 .. B.L loop
if B.V(N) >= Fac3 then
Show("Error: ball should be fixed");
end if;
end loop;
end CheckMovable;
procedure CheckMovable(M: in Positive; P: in PolysPtr; PLast: in Natural := 999) is
begin
for K in 1 .. Integer'Min(P'Last,PLast) loop
CheckMovable(M,P(K));
end loop;
end CheckMovable;
procedure AddConstraint(B0,B1: in BallPtr; I: in Integer; A: in out RMatPtr) is
DP: constant IPoint := B1.IP-B0.IP;
begin
A(I,0) := Half*(DP*DP);
A(I,3*B0.ID-2) := Pt(1,0,0)*DP;
A(I,3*B0.ID-1) := Pt(0,1,0)*DP;
A(I,3*B0.ID) := Pt(0,0,1)*DP;
if B1.ID=0 then return; end if;
A(I,3*B1.ID-2) := -A(I,3*B0.ID-2);
A(I,3*B1.ID-1) := -A(I,3*B0.ID-1);
A(I,3*B1.ID) := -A(I,3*B0.ID);
end AddConstraint;
function Constraints(Mov: Blob) return RMatPtr is
I: Integer := 0;
B0,B1: BallPtr;
A: RMatPtr := new RMat(1 .. Mov.M,0 .. 3*Mov.L);
begin
SetZero(A.all);
for M in 1 .. Mov.L loop
B0 := Mov.V(M);
for N in 1 .. 12 loop
B1 := B0.Nb(N);
if B0.ID>B1.ID then
I := I+1;
AddConstraint(B0,B1,I,A);
end if;
end loop;
end loop;
return A;
end Constraints;
procedure Cleanup(B: in PolyPtr) is
L: Integer := 1;
begin
if B.L<2 then return; end if;
SortInt(B.V,1,B.L);
for K in 2 .. B.L loop
if B.V(L) > B.V(K) then
L := L+1;
if K>L then B.V(L) := B.V(K); end if;
end if;
end loop;
B.L := L;
end Cleanup;
procedure Cleanup(B: in out PolysPtr; Size: in out Natural; AddUp: Boolean := True) is
--- assuming components are already cleaned up
Kmax: constant Integer := Size;
C: PolyPtr;
begin
if Size<2 then return; end if;
SortPolyPtr(B.all,1,Size);
Size := 1;
for K in 2 .. Kmax loop
if B(Size) > B(K) then
Size := Size+1;
if K>Size then
C := B(K);
B(K) := B(Size);
B(Size) := C;
end if;
elsif AddUp then
B(Size).M := B(Size).M+B(K).M;
end if;
end loop;
end Cleanup;
function FindConnAll(NM: Positive) return PolysPtr is
PDim: constant Positive := 2**16;
Size: Integer := 0;
P: PolysPtr := new Polys(1 .. PDim);
Mov: Balls(1 .. NM);
procedure AddP is
begin
Size := Size+1;
for N in 1 .. NM loop
Pack(Mov(N).IP,1,P(Size).V(N));
end loop;
P(Size).L := NM;
SortInt(P(Size).V,1,NM);
if Size=P'Last then Cleanup(P,Size,False); end if;
end AddP;
procedure AddMovable(K: in Integer) is
begin
for M in 1 .. K-1 loop
for N in 1 .. 12 loop
Mov(K) := Mov(M).Nb(N);
if Mov(K).ID=0 then
if K=NM then AddP;
else
Mov(K).ID := K;
AddMovable(K+1);
Mov(K).ID := 0;
end if;
end if;
end loop;
end loop;
end AddMovable;
begin
ClearAllMovable;
for N in 1 .. PDim loop
P(N) := new Poly(NM);
end loop;
Mov(1) := BallAt((0,0,0));
Mov(1).ID := 1;
AddMovable(2);
Cleanup(P,Size,False);
Trim(P,Size);
ClearAllMovable;
return P;
end FindConnAll;
function Addable(K: Int; B: PolyPtr) return Boolean is
Eps: constant Rep := 0.001;
begin
return B.L=0 or else abs(RDist(K,B)-One)<Eps;
end Addable;
pragma Inline (Addable);
function FindSubAll(NM: Positive; B: PolyPtr) return PolysPtr is
PDim: constant Positive := 2**12;
Size: Integer := 0;
P: PolysPtr := new Polys(1 .. PDim);
A: PolyPtr := new Poly(NM);
procedure AddP is
begin
Size := Size+1;
P(Size) := new Poly(NM);
P(Size).all := A.all;
SortInt(P(Size).V,1,NM);
end AddP;
procedure AddToBall(K: in Integer) is
begin
if K=NM then AddP; return; end if;
for I in 1 .. B.L loop
if Addable(B.V(I),A) then
A.L := A.L+1;
A.V(A.L) := B.V(I);
AddToBall(K+1);
A.L := A.L-1;
end if;
end loop;
end AddToBall;
begin
AddToBall(0);
Free(A);
Cleanup(P,Size,False);
Trim(P,Size);
return P;
end FindSubAll;
procedure AddNbors(B: in out PolyPtr) is
L: constant Integer := B.L;
M: Integer := 0;
A: BallPtr;
begin
for K in 1 .. L loop
if B.V(K) >= Fac3 then M := M+1; end if;
end loop;
AddSpace(12*M,B);
M := L;
for K in 1 .. L loop
if B.V(K) >= Fac3 then
A := BallAt(Unpack(B.V(K)));
for N in 1 .. 12 loop
M := M+1;
Pack(A.Nb(N).IP,A.Nb(N).ID,B.V(M));
end loop;
end if;
end loop;
B.L := M;
end AddNbors;
procedure AddNbors(P: in PolysPtr) is
begin
for K in P'Range loop
AddNbors(P(K));
end loop;
end AddNbors;
procedure RemoveNbors(B: in PolyPtr) is
begin
B.L := NumMovable(B);
end RemoveNbors;
procedure RemoveNbors(P: in PolysPtr) is
begin
for K in P'Range loop
RemoveNbors(P(K));
end loop;
end RemoveNbors;
function ToInt(I: Integer) return Int is
begin
if I >= 0 then
return Int(I);
else
return Int(I-Integer'First)-(Int(-1-Integer'First)+1);
end if;
end ToInt;
procedure ScrewHCP(L: in IPoint; B: in Poly; C: out Poly) is
begin
if (L(3) mod 2)=0 then Translate(L,B,C); return; end if;
CopyParam(B,C);
declare
N: Natural;
P: IPoint;
begin
for I in 1 .. B.L loop
UnPack(B.V(I),P,N);
Pack((P(2)+L(2),P(1)+L(1),P(3)+L(3)),N,C.V(I));
end loop;
end;
end ScrewHCP;
procedure NormalizeHCP(B: in out PolyPtr) is
begin
ScrewHCP(-UnPack(B.V(1)),B.all,B.all);
SortInt(B.V,1,B.L);
declare
A: constant Poly(B.V'Last) := B.all;
C: PolyPtr := new Poly(B.V'Last);
begin
for I in 1 .. A.L loop
exit when A.V(I)<Fac3;
if A.V(I) /= IZero then
ScrewHCP(-UnPack(A.V(I)),A,C.all);
SortInt(C.V,1,A.L);
if C>B then B.all := C.all; end if;
end if;
end loop;
Free(C);
end;
end NormalizeHCP;
procedure Normalize(B: in out PolyPtr) is
begin
SortInt(B.V,1,B.L);
if LayerPeriod=2 then NormalizeHCP(B); return; end if;
declare
T: constant Int := IZero-(B.V(1) and Mask3);
begin
for N in 1 .. B.L loop
B.V(N) := B.V(N)+T;
end loop;
end;
end Normalize;
procedure Normalize(P: in PolysPtr) is
begin
for K in P'Range loop
Normalize(P(K));
end loop;
end Normalize;
function Group(Gen: IMapVec) return IMapVecPtr is
Count: Integer := -1;
G: IMapVec(0 .. 1023);
procedure NextG(A: in IMap) is
begin
for N in 0 .. Count loop
if G(N)=A then return; end if;
end loop;
Count := Count+1;
G(Count) := A;
for N in 1 .. Gen'Last loop
NextG(A*Gen(N));
end loop;
end NextG;
begin
SetOne(G(0));
NextG(G(0));
if Verbose>0 then Show("symmetry group has order" & Count'Img & "+1"); end if;
declare
Grp: IMapVecPtr := new IMapVec(0 .. Count);
begin
for N in 0 .. Count loop
Grp(N) := G(N);
end loop;
return Grp;
end;
end Group;
function SymmHCP return IMapVecPtr is
Gen: IMapVec(1 .. 3); -- generators of group
begin
SetZero(Gen(1)); -- Pz
SetRow(1,Pt( 1, 0, 0),Gen(1));
SetRow(2,Pt( 0, 1, 0),Gen(1));
SetRow(3,Pt( 0, 0,-1),Gen(1));
SetZero(Gen(2)); -- Ry
SetRow(1,Pt( 0,-1, 0),Gen(2));
SetRow(2,Pt(-1, 0, 0),Gen(2));
SetRow(3,Pt( 0, 0,-1),Gen(2));
SetZero(Gen(3));
SetRow(1,Pt( 0, 1, 0),Gen(3));
SetRow(2,Pt(-1,-1, 0),Gen(3));
SetRow(3,Pt( 0, 0, 1),Gen(3));
return Group(Gen);
end SymmHCP;
function SymmFCC return IMapVecPtr is
Gen: IMapVec(1 .. 3); -- generators of group
begin
SetZero(Gen(1)); -- Rx
SetRow(1,Pt( 0, 1, 0),Gen(1));
SetRow(2,Pt( 1, 0, 0),Gen(1));
SetRow(3,Pt( 0, 0,-1),Gen(1));
SetZero(Gen(2)); -- Pz*Rz
SetRow(1,Pt( 1, 1, 0),Gen(2));
SetRow(2,Pt(-1, 0, 0),Gen(2));
SetRow(3,Pt( 0, 0,-1),Gen(2));
SetZero(Gen(3)); -- Pz*Ra
SetRow(1,Pt( 2, 1, 4),Gen(3));
SetRow(2,Pt( 1, 2,-4),Gen(3));
SetRow(3,Pt( 1,-1,-1),Gen(3)); Divide(Gen(3),3);
return Group(Gen);
end SymmFCC;
function Symm return IMapVecPtr is
begin
if LayerPeriod=2 then return SymmHCP; end if;
if LayerPeriod=3 then return SymmFCC; end if;
Message("Lattice.Ops.Symm",Not_Implemented);
return null;
end Symm;
procedure Transform(A: in IMap; B: in Poly; C: out Poly) is
M: Integer;
P: IPoint;
begin
CopyParam(B,C);
for K in 1 .. B.L loop
UnPack(B.V(K),P,M);
Pack(A*P,M,C.V(K));
end loop;
SortInt(C.V,1,C.L);
end Transform;
procedure Symmetrize(B: in out PolyPtr) is
--- assuming B sorted
A: constant Poly(B.V'Last) := B.all;
C: PolyPtr := new Poly(B.V'Last);
begin
for I in 1 .. Symmetries'Last loop
Transform(Symmetries(I),A,C.all);
Normalize(C);
if C>B then B.all := C.all; end if;
end loop;
Free(C);
end Symmetrize;
procedure Symmetrize(P: in out PolysPtr) is
--- assuming P is normalized
Size: Integer := P'Last;
begin
for K in 1 .. Size loop
Symmetrize(P(K));
end loop;
Cleanup(P,Size);
Trim(P,Size);
end Symmetrize;
function FindConnClass(NM: Positive) return PolysPtr is
P: PolysPtr := FindConnAll(NM);
Size: Integer := P'Last;
begin
if Verbose>0 then Show("starting with:" & Size'Img & " configurations"); end if;
AddNbors(P);
Normalize(P);
Cleanup(P,Size,False);
if Verbose>0 then Show("after translation:" & Size'Img & " configurations"); end if;
Trim(P,Size);
Symmetrize(P);
Size := P'Last;
if Verbose>0 then Show("after symmetrizing:" & Size'Img & " configurations"); end if;
return P;
end FindConnClass;
function FindSubClass(NM: Positive; B: PolyPtr) return PolysPtr is
P: PolysPtr := FindSubAll(NM,B);
Size: Integer := P'Last;
begin
if Verbose>1 then Show("starting with:" & Size'Img & " configurations"); end if;
AddNbors(P);
Normalize(P);
Cleanup(P,Size);
if Verbose>1 then Show("after translation:" & Size'Img & " configurations"); end if;
Trim(P,Size);
Symmetrize(P);
Size := P'Last;
if Verbose>1 then Show("after symmetrizing:" & Size'Img & " configurations"); end if;
return P;
end FindSubClass;
procedure AddEntropies(N: in Positive; Pn,P: in PolysPtr) is
R: Rep;
B: PolysPtr;
begin
if N=1 then
R := BasicVol;
for K in P'Range loop
P(K).R := P(K).R-Rep(P(K).L)*R;
end loop;
return;
end if;
for K in P'Range loop
B := FindSubClass(N,P(K));
RemoveNbors(B);
for I in B'Range loop
GetR(B(I),Pn,R);
P(K).R := P(K).R-Rep(B(I).M)*R;
end loop;
FreeAll(B);
end loop;
end AddEntropies;
function Entropies(NM: in Positive) return PExpPtr is
E: PExpPtr;
begin
E := StoredVols(LayerPeriod,NM);
for N in 1 .. NM loop
for K in N+1 .. NM loop
AddEntropies(N,E(N),E(K));
end loop;
end loop;
return E;
end Entropies;
procedure EntropySums(NM: in Positive) is
E: PExpPtr := Entropies(NM);
R,Rn: Rep := Zero;
begin
if LayerPeriod=2 then Show("HCP lattice:"); end if;
if LayerPeriod=3 then Show("FCC lattice:"); end if;
for N in 1 .. NM loop
Rn := Zero;
for K in E(N)'Range loop
Rn := Rn+Rep(E(N)(K).M)*E(N)(K).R;
end loop;
R := R+Rn;
Show("size" & N'Img & ": correction =" & Rn'Img,False);
Show(" entropy =" & R'Img);
end loop;
end EntropySums;
end Lattice.Ops;