CRYSTALLITE

A program to grow a digital 3D crystal segment, with growth based on spiral defects. Based in part on the FORTRAN progam FACE and the Mathematica program BLOCK by Dan Walker and Jeanne Paquette, and in part on the Mathematica program Simple Volumetric Modeling ("SMV") (http://library.wolfram.com/examples/3DHead/)
Version 1a    Begun 1 August 2001  
Version 1b    Begun 24 August 2001: Length of defect and boundary at which defect terminates are fixed at time 0.
Version 1c    Begun 24 August 2001: laymax made expandable; graphics cleaned up.
Synopsis
Crystal and solution are represented by an array xtl of cellular automata. xtl begins as a strip along the defect, and grows as required. The defect length and the boundary against which the defect terminates are fixed at time 0.
The state of each automaton depends on whether it is solution or crystal. The solution state is represented by the code 0 (red Hue). The crystal state is represented by a 2 digit code which depends on the Hue[h] value we want to assign to a particular ion (CO3, Ca , or Mn) and growth direction (base case or growth towards the east, south, or southeast --- N, W, WN, SWN, and WNE sites). See Mathematic Book, p. 486 and p. 994 (colour plate) for details on Hue[h].
Growth is based on a screw defect, which generates an Archimedes-spiral growth mound. The approach taken here is to ignore the defect within the 3D array itself. Each layer of the array represents a "spiral layer" of the crystal. The defect is simply stored as coordinates. This data is utilized in assigning probability of transition (see below) and for visualization.
Transitions happen at the crystal/solution interface. This requires maintaining a 2D array surface of the current interface. We only implement crystallization transitions here. 
For crystallization transitions, a candidate solution automaton at the interface is chosen at random. First, the candidate is assessed for eligibility:  all nearest-neighbour sites in the underlying layer must be occupied (i.e., no transitions are permitted if they would create a "cliff" more than one growth unit high). The candidate is then assigned: (1) a ion type, based on position (CO3 or cation) and, if cation, Ca or Mn, based on desired Ca:Mn ratio (e.g., 10:1); and (2) a probability of attachment. The probability of attachment is based on: (1) ion type, especially Ca vs Mn; (2) whether or not the candidate lies on the defect line; and (3) the configuration of the 4 nearest neighbours (15 possibilites, based on combinations of N, E, S, and W). The transition occurs if the probability of attachment is greater than a randomly-chosen value.
[For dissolution transitions, the process would be reversed. A candidate crystal automaton at the interface would be chosen at random, and assigned a probability of detachment, based on (1) state of the candidate; (2) whether or not the candidate lies on the defect line; and (3) the configuration of the 4 nearest neighbours. The transition would occur if the probability of detachment is greater than a randomly-chosen value.]
[Transition probabilities could also be made time-dependent, e.g., to represent changing saturation in the solution.]
Note that a single Ca growth unit is about 10^-10 m (1 Angstrom), so we would have to get up to 1000x1000 growth units (which is feasible) just to show an area 10^-7 m (0.1 microns) on a side.

\( (*SET\ VARIABLES\ FOR\ THIS\ RUN*) 

 (*Set\ the\ working\ directory*) 
\(SetDirectory["D://CRYSTALLITE/30aug01"];\)\  (*this\ determines\ where\ output\ files\ will\ be\ saved*) 

 (*For\ this\ first\ pass, \ 
      the\ defect\ is\ 6  
        E*) 

 (*This\ section\ sets\ up\ the\ defect\ as\ 6  E*) \)
defectlength = 500 ;
defectdowncol = 2;
defectupcol = defectdowncol - 1 ;
defectrowmax = defectlength;
\(\(defectrowmin\  = 
        1;\)

 (*This\ controls\ the\ number\ of\ iterations*) \)
maxhits = 20000000; 
hitcounterctrl := 
        maxhits/10;

  
mnca = 
        0.1 ;

 (*This is the list of attachment probabilities for different sites and ions .  
          C is for carbonate,  K is for calcium,  
      M is for manganese .  
          N indicates that the site north of the candidate site is occupied,  
      EN indicates that sites to the east and north of the candidate site are occupied .  Ordering convention :  
          clockwise from first occupied site.*) 
 (*THESE HAVE BEEN ALTERED 30  AUG,  
      TO TRY TO GET A POLYGONIZED MOUND*) 
NC = 0.3;
EC = 0.01;
SC = EC;
WC = NC;
NEC = 0.015;
ESC = 0.015;
SWC = NEC;
WNC = 0.9; 
NSC = 0.015;
EWC = NSC;
NESC = 0.02;
ESWC = 0.02;
SWNC = 0.999;
WNEC = 0.999;
NESWC = 0.05;
NK = 0.3; 
EK = 0.3;
SK = EK;
WK = NK;
NEK = 0.9;
ESK = 0.9;
SWK = NEK;
WNK = 0.9;
NSK = 0.9;
EWK = NSK;
NESK = 0.99;
ESWK = 0.99;
SWNK = 0.99;
WNEK = 0.99;
NESWK = 0.9999;
NM = 0.1; 
EM = 0.1;
SM = EM;
WM = NM;
NEM = 0.9;
ESM = 0.9;
SWM = NEM;
WNM = 0.9; 
NSM = 0.9;
EWM = NSM;
NESM = 0.99;
ESWM = 0.99;
SWNM = 0.99;
WNEM = 0.99;
NESWM = 
        0.999;

  
solncolour = 
        50 ;  
CO3colour = 
        15;   
Cacolour = 65;  
Mncolour = 85 ;  
secolour = 
        2 ; (*should be 1 or 2*) 

 (*INITIALIZE ARRAYS AND COUNTERS*) 

laymax = 2;
rowmax = defectlength;
colmax = 
        2;

  
xtl = 
        Table [solncolour, {layer, laymax},  {row, rowmax}, {col, 
            colmax}];
  
surface = 
        Table[0, {row, rowmax}, {col, 
            colmax}];

  

\(\(Do[
 (*Let\ us\ know\ how\ many\ iterations\ have\ been\ done, \ 
          and\ save\ what' 
            s\ been\ done*) 
If[
Mod\ [i, 
              hitcounterctrl] == 
            1, 
Print["Total no. of hits: ", 
            maxhits, ". Hits completed: ", i - 1, ". Layers: ", 
            laymax - 1, ". Rows: ", rowmax, ". Cols: ", 
            colmax]; 
Put[
            xtl, "xtl_out.m"]; 
Put[
            surface, "surface_out.m"]; 
Put[
            colmax, "colmax_out.m"]; 
Put[
            rowmax, "rowmax_out.m"]; 
Put[
            laymax, "laymax_out.m"]; 
Put[
            defectlength, "defectlength_out.m"]; 
Put[
            defectrowmax, "defectrowmax_out.m"]; 
Put[
            defectrowmin, "defectrowmin_out.m"]; 
Put[
            defectdowncol, "defectdowncol_out.m"]; 
Put[
            defectupcol, "defectupcol_out.m"]; 
Put[
            i - 1, "hits_out.m"];
]; 

 (*This\ section\ picks\ a\ random\ candidate\ site*) 
colcand = 
          Random[Integer, {1, colmax}]; 
rowcand = 
          Random[Integer, {1, rowmax}]; 
basecand = \ 
          surface[\([rowcand, colcand]\)]; 
laycand\  = 
          basecand + 
            1; 

 (*This\ section\ figures\ out\ the\ neighbours\ of\ the\ candidate\ site*) 
 (*\(...\) \(neighbour\)\(\ \)\(to\)\(\ \)\(the\)\(\ \)\(north\)*) 
If[rowcand < rowmax, 
          Ncand = surface[\([rowcand + 1, colcand]\)] - basecand, 
          Ncand = 0]; 
 (*\(...\) \(neighbour\)\(\ \)\(to\)\(\ \)\(the\)\(\ \)\(south\)*) 
If[rowcand\  > 1, 
          Scand = surface[\([rowcand - 1, colcand]\)] - basecand, 
          Scand = 0]; 
 (*\(...\) \(neighbour\)\(\ \)\(to\)\(\ \)\(the\)\(\ \)\(east\)*) 
Which[
\((colcand < 
                  colmax\  && 
                colcand ≠ defectupcol)\) || \((colcand == 
                  defectupcol && rowcand > defectrowmax)\), 
          Ecand = surface[\([rowcand, colcand + 1]\)] - basecand, 
          colcand == defectupcol && rowcand <= defectrowmax, 
          Ecand = surface[\([rowcand, colcand + 1]\)] - basecand - 
              1, 
colcand == colmax, 
          Ecand = 0
]; 
 (*\(\(...\) \(neighbour\)\(\ \)\(to\)\(\ \)\(the\)\(\ \)\(west\)\) ... *) 
Which[
\((colcand > 1\  && 
                colcand ≠ defectdowncol)\) || \((colcand == 
                  defectdowncol && rowcand > defectrowmax)\)\ , 
          Wcand = surface[\([rowcand, colcand - 1]\)] - 
              basecand, 
colcand == defectdowncol && 
            rowcand <= defectrowmax, 
          Wcand = surface[\([rowcand, colcand - 1]\)] - basecand + 
              1, 
colcand == 1, 
          Wcand = 0
]; 

 (*This\ drops\ out\ of\ the\ loop\ if\ there\ is\ nothing\ to\ stick\ to*) 
If[
          Ncand == 0 && Ecand == 0 && Scand == 0 && Wcand == 0, 
          Continue[]]; 

 (*This\ drops\ out\ of\ the\ loop\ if\ there\ are\ any\  - 1\ candidates, \ i . e . , \ 
          if\ sticking\ would\ create\ a\ cliff\ 2\ units\ high*) 
If[Ncand < 0 || Ecand < 0 || Scand < 0 || Wcand < 0, 
          Continue[]]; 

 (*Set\ the\ iontype*) 
If[OddQ[colcand + rowcand + laycand], 
          iontype = CO3colour, 
          If[Random[] > mnca, iontype = Cacolour, 
            iontype = 
              Mncolour]]; 

 (*Reset\ probability\ of\ attachment\ to\ zero*) \ 
pa = 
          0.0; 

 (*Determine\ probability\ of\ attachment, \ 
          depending\ on\ ion\ \ type\ and\ \(\(neighbours\)\(.\)\)*) 
If[iontype == CO3colour, 
          Which[\((Ncand == 1 && Ecand == 0 && Scand == 0 && 
                Wcand == 0)\), {pa = NC, 
              iontype = iontype + secolour}, \((Ecand == 1 && 
                Scand == 0 && Wcand == 0 && Ncand == 0)\), 
            pa = EC, \((Scand == 1 && Wcand == 0 && Ncand == 0 && 
                Ecand == 0)\), 
            pa = SC, \((Wcand == 1 && Ncand == 0 && Ecand == 0 && 
                Scand == 0)\), {pa = WC, 
              iontype = iontype + secolour}, \((Ncand == 1 && 
                Ecand == 1 && Scand == 0 && Ecand == 0)\), 
            pa = NEC, \((Ecand == 1 && Scand == 1 && Wcand == 0 && 
                Ncand == 0)\), 
            pa = ESC, \((Scand == 1 && Wcand == 1 && Ncand == 0 && 
                Ecand == 0)\), 
            pa = SWC, \((Wcand == 1 && Ncand == 1 && Scand == 0 && 
                Ecand == 0)\), {pa = WNC, 
              iontype = iontype + secolour}, \((Ncand == 1 && 
                Scand == 1 && Ecand == 0 && Wcand == 0)\), 
            pa = NSC, \((Ecand == 1 && Wcand == 1 && Ncand == 0 && 
                Scand == 0)\), 
            pa = EWC, \((Ncand == 1 && Ecand\  == 1 && Scand == 1 && 
                Wcand == 0)\), 
            pa = NESC, \((Ecand == 0 && Scand\  == 0 && Wcand == 3 && 
                Ncand == 0)\), 
            pa = ESWC, \((Scand == 1 && Wcand\  == 1 && 
                Ncand == 1 && Ecand == 0)\), {pa = SWNC, 
              iontype = iontype + secolour}, \((Wcand == 1 && Ncand\  == 1 && 
                Ecand == 1 && Scand == 0)\), {pa = WNEC, 
              iontype = iontype + secolour}, \((Ncand == 1 && Ecand\  == 1 && 
                Scand == 1 && Wcand == 1)\), 
            pa = NESWC]]; 
If[iontype == Cacolour, 
          Which[\((Ncand == 1 && Ecand == 0 && Scand == 0 && 
                Wcand == 0)\), {pa = NK, 
              iontype = iontype + secolour}, \((Ecand == 1 && 
                Scand == 0 && Wcand == 0 && Ncand == 0)\), 
            pa = EK, \((Scand == 1 && Wcand == 0 && Ncand == 0 && 
                Ecand == 0)\), 
            pa = SK, \((Wcand == 1 && Ncand == 0 && Ecand == 0 && 
                Scand == 0)\), {pa = WK, 
              iontype = iontype + secolour}, \((Ncand == 1 && 
                Ecand == 1 && Scand == 0 && Ecand == 0)\), 
            pa = NEK, \((Ecand == 1 && Scand == 1 && Wcand == 0 && 
                Ncand == 0)\), 
            pa = ESK, \((Scand == 1 && Wcand == 1 && Ncand == 0 && 
                Ecand == 0)\), 
            pa = SWK, \((Wcand == 1 && Ncand == 1 && Scand == 0 && 
                Ecand == 0)\), {pa = WNK, 
              iontype = iontype + secolour}, \((Ncand == 1 && 
                Scand == 1 && Ecand == 0 && Wcand == 0)\), 
            pa = NSK, \((Ecand == 1 && Wcand == 1 && Ncand == 0 && 
                Scand == 0)\), 
            pa = EWK, \((Ncand == 1 && Ecand\  == 1 && Scand == 1 && 
                Wcand == 0)\), 
            pa = NESK, \((Ecand == 1 && Scand\  == 1 && Wcand == 1 && 
                Ncand == 0)\), 
            pa = ESWK, \((Scand == 1 && Wcand == 1 && Ncand == 1 && 
                Ecand == 0)\), {pa = SWNK, 
              iontype = iontype + secolour}, \((Wcand == 1 && Ncand\  == 1 && 
                Ecand == 1 && Scand == 0)\), {pa = WNEK, 
              iontype = iontype + secolour}, \((Ncand == 0 && Ecand\  == 0 && 
                Scand == 0 && Wcand == 4)\), 
            pa = NESWK]]; 
If[iontype == Mncolour, 
          Which[\((Ncand == 1 && Ecand == 0 && Scand == 0 && 
                Wcand == 0)\), {pa = NM, 
              iontype = iontype + secolour}, \((Ecand == 1 && 
                Scand == 0 && Wcand == 0 && Ncand == 0)\), 
            pa = EM, \((Scand == 1 && Wcand == 0 && Ncand == 0 && 
                Ecand == 0)\), 
            pa = SM, \((Wcand == 1 && Ncand == 0 && Ecand == 0 && 
                Scand == 0)\), {pa = WM, 
              iontype = iontype + secolour}, \((Ncand == 1 && 
                Ecand == 1 && Scand == 0 && Ecand == 0)\), 
            pa = NEM, \((Ecand == 1 && Scand == 1 && Wcand == 0 && 
                Ncand == 0)\), 
            pa = ESM, \((Scand == 1 && Wcand == 1 && Ncand == 0 && 
                Ecand == 0)\), 
            pa = SWM, \((Wcand == 1 && Ncand == 1 && Scand == 0 && 
                Ecand == 0)\), {pa = WNM, 
              iontype = iontype + secolour}, \((Ncand == 1 && 
                Scand == 1 && Ecand == 0 && Wcand == 0)\), 
            pa = NSM, \((Ecand == 1 && Wcand == 1 && Ncand == 0 && 
                Scand == 0)\), 
            pa = EWM, \((Ncand == 1 && Ecand\  == 1 && Scand == 1 && 
                Wcand == 0)\), 
            pa = NESM, \((Ecand == 1 && Scand\  == 1 && Wcand == 1 && 
                Ncand == 0)\), 
            pa = ESWM, \((Scand == 1 && Wcand\  == 1 && Ncand == 1 && 
                Ecand == 0)\), {pa = SWNM, 
              iontype = iontype + secolour}, \((Wcand == 1 && Ncand\  == 1 && 
                Ecand == 1 && Scand == 0)\), {pa = WNEM, 
              iontype = iontype + secolour}, \((Ncand == 1 && Ecand\  == 1 && 
                Scand == 1 && Wcand == 1)\), 
            pa = NESWM]]; 

 (*Determine\ whether\ the\ ion\ does\ attach, \ 
          and\ update\ arrays\ if\ it\ does*) 
If[
          Random[] ≤ 
            pa, 
\((xtl[\([laycand, rowcand, colcand]\)] = 
              iontype; \ surface[\([rowcand, colcand]\)] = laycand)\), 
          Continue[]
]; 

 (*This\ section\ gives\ an\ animation\ frame\ for\ every\ added\ ion, \ 
          as\ a\ bug\ check*) 
 (*\(If[
              maxhits ≤ 1000, 
p1 = 
                ListDensityPlot[Mod[xtl[\([1, All, All]\)], 10]/10 +  .5, 
                  colourplan, 
                  Epilog -> {Line[{{defectupcol, 
                            defectrowmin - 1}, {defectupcol, 
                            defectrowmax}}]}]; 
p2 = 
                ListDensityPlot[Mod[xtl[\([2, All, All]\)], 10]/10 +  .5, 
                  colourplan, 
                  Epilog -> {Line[{{defectupcol, 
                            defectrowmin - 1}, {defectupcol, 
                            defectrowmax}}]}]; 
p3 = 
                ListDensityPlot[Mod[xtl[\([3, All, All]\)], 10]/10 +  .5, 
                  colourplan, 
                  Epilog -> {Line[{{defectupcol, 
                            defectrowmin - 1}, {defectupcol, 
                            defectrowmax}}]}]; 
Show[
                GraphicsArray[{p1, p2, p3}, 
                  ImageSize -> 72*colmax]]; 
Clear[p1, 
                p2, p3];
];\)*) 

 (*Now\ add\ new\ row\ or\ column\ if\ needed*) 
If\ [
          colcand == colmax, 
colmax = 
            colmax + 1; 
xtl = 
            PadRight[xtl, {laymax, rowmax, colmax}, 
              solncolour]; 
surface = 
            PadRight[surface, {rowmax, colmax}, 
              0];
]; 
If\ [
          rowcand == rowmax, 
rowmax = 
            rowmax + 1; 
xtl = 
            PadRight[xtl, {laymax, rowmax, colmax}, 
              solncolour]; 
surface = 
            PadRight[surface, {rowmax, colmax}, 
              0];
]; 
If\ [
          colcand == 1, 
colmax = 
            colmax + 1; 
defectdowncol = 
            defectdowncol + 1; 
defectupcol = 
            defectupcol + 1; 
xtl = 
            PadLeft[xtl, {laymax, rowmax, colmax}, 
              solncolour]; 
surface = 
            PadLeft[surface, {rowmax, colmax}, 
              0];
]; 
 (*The\ following\ is\ commented\ out\ because\ this\ boundary\ is\ fixed . \ 
              This\ is\ a\ temporary\ fix, \ 
          and\ applies\ only\ to\ cases\ in\ which\ the\ defect\ terminates\ at\ the\ south\ boundary*) 
 (*\(If\ [
              rowcand == 1, 
rowmax = 
                rowmax + 1; 
defectrowmax = 
                defectrowmax + 1; 
xtl = 
                PadLeft[xtl, {laymax, rowmax, colmax}, 
                  solncolour]; 
surface = 
                PadLeft[surface, {rowmax, colmax}, 
                  0];
];\)*) 
If\ [
          laycand == laymax, 
laymax = 
            laymax + 1; 
xtl = 
            PadRight[xtl, {laymax, rowmax, colmax}, 
              solncolour];
];
 (*This\ is\ the\ end\ of\ the\ iteration*) 
, {i, 
          maxhits}];\)

 (*PUT\ THE\ IMPORTANT\ STUFF\ AWAY*) \)
Put[xtl, "xtl_out.m"];
Put[surface, "surface_out.m"];
Put[colmax, "colmax_out.m"];
Put[rowmax, "rowmax_out.m"];
Put[laymax, "laymax_out.m"];
Put[
        defectlength, "defectlength_out.m"];
Put[
        defectrowmax, "defectrowmax_out.m"];
Put[
        defectrowmin, "defectrowmin_out.m"];
Put[
        defectdowncol, "defectdowncol_out.m"];
Put[
        defectupcol, "defectupcol_out.m"];
Put[maxhits, "hits_out.m"];