Understanding: w_layers.n() > w_max_stacks_allowed

asbridger

Hi all,

Very much a TOPAS novice here, so I apologise if this is a common problem or I have just made an obvious mistake but I've spent a decent amount of time trying to solve this and can't seem to find a suitable solution.

The issue I am having is with the stacking fault functionality for more complex stacking dependencies. For my particular use case, this might not even be the correct way to approach things but I have made a more simple example of the problem and I was hoping someone could direct me as to what is going wrong and if there is any fix (potentially a keyword for a tolerance etc that I am missing?)

My example is just based on the cu_tutorial_01.inp

As far as I can tell my input file should be producing equivalent super cells to the tutorial (but maybe I am mistaken and this is the root of my issues). The change I have made is that instead of calculating transitions on a layer by layer basis (A to A, A to B etc.), I am calculating in groups of 3 layers (AAA to AAA, AAA to ABA, etc.). The overall probabilities for these transitions, I think, are unchanged, just individual transition probabilities are now much smaller.

Eg. to get AAA stacking on top of an already made AAA: in the tutorial this is 3 A to A transitions each with a probability of pa, in my example this is one transition of AAA to AAA with a probability of papapa. I think the problem seems to be related to the fact that in the example if pa = 0.9 there are only two transitions from each layer with probabilities 0.9 and 0.1, where as for my test there are now eight transitions with the smallest probability having a value of 0.001.

I have tried using:

pr_str { B layers_tol } 1

but this doesn't fix the issue, the refinement does not run and I still get:

Seq 87 : Error in generate_stack_sequences : w_layers.n() > w_max_stacks_allowed 35073 35072

I can set "pr_str { B layers_tol } 1.5" to get the code to run, but the fit to the data is no longer any good so I am assuming this is no longer a productive change for values above 1?

If anyone has any insight to what I am doing wrong/ what I need to change to get this to work, I would be incredibly grateful.

My full test input file is attached

cu-tutorial-triplestack.txt

9kB

and is also pasted below.

Thanks in advance,

Andy

'input file for Cu data simulated in diffax with cu.dat 'this is the standard file distributed with diffax 'simulation from 2 to 170 degrees with step 0.01 and full adaptive quadrature options 'xye file generated from column 3 of cu.spc using I^0.5 for error iters 100 seed xdd cu.xye LP_Factor(0) bkg 1 start_X 10 finish_X 149.6 rebin_with_dx_of 0.01 lam ymin_on_ymax 0.0001 la 1 lo 1.5418 lh 1e-5 Zero_Error(@,-0.00733)
weighting 1
str
for strs {
space_group P1
a 2.5560 b 2.55560 c = Get(generated_c); ga 120
prm !pa 0.9 'probability of having a single AA or BB type transition
prm !h 6.26
prm s 4875.84407` min 1e-15 scale = s 1e-6 / (Nv Nstr);
generate_stack_sequences {
number_of_sequences Nstr 100
number_of_stacks_per_sequence Nv 133

    'pr_str { B layers_tol } 1.5 'can set this to 1.5 to get it to run
    
    Transition(AAA, AAA, h) 
      'all AA or BB transitions
      to AAA = pa*pa*pa;   a_add = 0; b_add = 0;              
      'one AB/BA type transition
      to AAB = pa*pa*(1-pa); a_add = 0; b_add = 0; 
      to ABB = pa*(1-pa)*pa; a_add = 0; b_add = 0; 
      to BBB = (1-pa)*pa*pa;   a_add = 0; b_add = 0;   
      'two AB/BA type transitions
      to ABA = pa*(1-pa)*(1-pa); a_add = 0; b_add = 0; 
      to BAA = (1-pa)*(1-pa)*pa; a_add = 0; b_add = 0; 
      to BBA = (1-pa)*pa*(1-pa); a_add = 0; b_add = 0; 
      'all AB/BA type transitions
      to BAB = (1-pa)*(1-pa)*(1-pa); a_add = 0; b_add = 0; 
      
    Transition(ABA, ABA, h) 
      'all AA or BB transitions
      to AAA = pa*pa*pa;   a_add = 0; b_add = 0;              
      'one AB/BA type transition
      to AAB = pa*pa*(1-pa); a_add = 0; b_add = 0; 
      to ABB = pa*(1-pa)*pa; a_add = 0; b_add = 0; 
      to BBB = (1-pa)*pa*pa;   a_add = 0; b_add = 0;   
      'two AB/BA type transitions
      to ABA = pa*(1-pa)*(1-pa); a_add = 0; b_add = 0; 
      to BAA = (1-pa)*(1-pa)*pa; a_add = 0; b_add = 0; 
      to BBA = (1-pa)*pa*(1-pa); a_add = 0; b_add = 0; 
      'all AB/BA type transitions
      to BAB = (1-pa)*(1-pa)*(1-pa); a_add = 0; b_add = 0;         

    Transition(BAA, BAA, h) 
      'all AA or BB transitions
      to AAA = pa*pa*pa;   a_add = 0; b_add = 0;              
      'one AB/BA type transition
      to AAB = pa*pa*(1-pa); a_add = 0; b_add = 0; 
      to ABB = pa*(1-pa)*pa; a_add = 0; b_add = 0; 
      to BBB = (1-pa)*pa*pa;   a_add = 0; b_add = 0;   
      'two AB/BA type transitions
      to ABA = pa*(1-pa)*(1-pa); a_add = 0; b_add = 0; 
      to BAA = (1-pa)*(1-pa)*pa; a_add = 0; b_add = 0; 
      to BBA = (1-pa)*pa*(1-pa); a_add = 0; b_add = 0; 
      'all AB/BA type transitions
      to BAB = (1-pa)*(1-pa)*(1-pa); a_add = 0; b_add = 0;    
  
    Transition(BBA, BBA, h) 
      'all AA or BB transitions
      to AAA = pa*pa*pa;   a_add = 0; b_add = 0;              
      'one AB/BA type transition
      to AAB = pa*pa*(1-pa); a_add = 0; b_add = 0; 
      to ABB = pa*(1-pa)*pa; a_add = 0; b_add = 0; 
      to BBB = (1-pa)*pa*pa;   a_add = 0; b_add = 0;   
      'two AB/BA type transitions
      to ABA = pa*(1-pa)*(1-pa); a_add = 0; b_add = 0; 
      to BAA = (1-pa)*(1-pa)*pa; a_add = 0; b_add = 0; 
      to BBA = (1-pa)*pa*(1-pa); a_add = 0; b_add = 0; 
      'all AB/BA type transitions
      to BAB = (1-pa)*(1-pa)*(1-pa); a_add = 0; b_add = 0;    
      
      
      
    Transition(AAB, AAB, h) 
      'all AA or BB transitions
      to BBB = pa*pa*pa;   a_add = 0; b_add = 0;              
      'one AB/BA type transition
      to BBA = pa*pa*(1-pa); a_add = 0; b_add = 0; 
      to BAA = pa*(1-pa)*pa; a_add = 0; b_add = 0; 
      to AAA = (1-pa)*pa*pa;   a_add = 0; b_add = 0;   
      'two AB/BA type transitions
      to BAB = pa*(1-pa)*(1-pa); a_add = 0; b_add = 0; 
      to ABB = (1-pa)*(1-pa)*pa; a_add = 0; b_add = 0; 
      to AAB = (1-pa)*pa*(1-pa); a_add = 0; b_add = 0; 
      'all AB/BA type transitions
      to ABA = (1-pa)*(1-pa)*(1-pa); a_add = 0; b_add = 0; 
      
      
    Transition(ABB, ABB, h) 
      'all AA or BB transitions
      to BBB = pa*pa*pa;   a_add = 0; b_add = 0;              
      'one AB/BA type transition
      to BBA = pa*pa*(1-pa); a_add = 0; b_add = 0; 
      to BAA = pa*(1-pa)*pa; a_add = 0; b_add = 0; 
      to AAA = (1-pa)*pa*pa;   a_add = 0; b_add = 0;   
      'two AB/BA type transitions
      to BAB = pa*(1-pa)*(1-pa); a_add = 0; b_add = 0; 
      to ABB = (1-pa)*(1-pa)*pa; a_add = 0; b_add = 0; 
      to AAB = (1-pa)*pa*(1-pa); a_add = 0; b_add = 0; 
      'all AB/BA type transitions
      to ABA = (1-pa)*(1-pa)*(1-pa); a_add = 0; b_add = 0; 
      
    Transition(BAB, BAB, h) 
      'all AA or BB transitions
      to BBB = pa*pa*pa;   a_add = 0; b_add = 0;              
      'one AB/BA type transition
      to BBA = pa*pa*(1-pa); a_add = 0; b_add = 0; 
      to BAA = pa*(1-pa)*pa; a_add = 0; b_add = 0; 
      to AAA = (1-pa)*pa*pa;   a_add = 0; b_add = 0;   
      'two AB/BA type transitions
      to BAB = pa*(1-pa)*(1-pa); a_add = 0; b_add = 0; 
      to ABB = (1-pa)*(1-pa)*pa; a_add = 0; b_add = 0; 
      to AAB = (1-pa)*pa*(1-pa); a_add = 0; b_add = 0; 
      'all AB/BA type transitions
      to ABA = (1-pa)*(1-pa)*(1-pa); a_add = 0; b_add = 0; 
      
    Transition(BBB, BBB, h) 
      'all AA or BB transitions
      to BBB = pa*pa*pa;   a_add = 0; b_add = 0;              
      'one AB/BA type transition
      to BBA = pa*pa*(1-pa); a_add = 0; b_add = 0; 
      to BAA = pa*(1-pa)*pa; a_add = 0; b_add = 0; 
      to AAA = (1-pa)*pa*pa;   a_add = 0; b_add = 0;   
      'two AB/BA type transitions
      to BAB = pa*(1-pa)*(1-pa); a_add = 0; b_add = 0; 
      to ABB = (1-pa)*(1-pa)*pa; a_add = 0; b_add = 0; 
      to AAB = (1-pa)*pa*(1-pa); a_add = 0; b_add = 0; 
      'all AB/BA type transitions
      to ABA = (1-pa)*(1-pa)*(1-pa); a_add = 0; b_add = 0; 
      
 
      
  'site Cu1 x =1/3; y =2/3; z = 0.0/Nv; occ Cu 1 beq @ -0.49524` layer A 
  'site Cu1 x =2/3; y =1/3; z = 0.0/Nv; occ Cu 1 beq @ -0.49524` layer B 
  
  site Cu1 x =1/3; y =2/3; z = 0.0/Nv;   occ Cu 1 beq @ -0.49524`  layer AAA
  site Cu2 x =1/3; y =2/3; z = (1/3)/Nv; occ Cu 1 beq @ -0.49524`  layer AAA
  site Cu3 x =1/3; y =2/3; z = (2/3)/Nv; occ Cu 1 beq @ -0.49524`  layer AAA

  site Cu1 x =2/3; y =1/3; z = 0.0/Nv;   occ Cu 1 beq @ -0.49524`  layer BAA
  site Cu2 x =1/3; y =2/3; z = (1/3)/Nv; occ Cu 1 beq @ -0.49524`  layer BAA
  site Cu3 x =1/3; y =2/3; z = (2/3)/Nv; occ Cu 1 beq @ -0.49524`  layer BAA      

  site Cu1 x =1/3; y =2/3; z = 0.0/Nv;   occ Cu 1 beq @ -0.49524`  layer ABA
  site Cu2 x =2/3; y =1/3; z = (1/3)/Nv; occ Cu 1 beq @ -0.49524`  layer ABA
  site Cu3 x =1/3; y =2/3; z = (2/3)/Nv; occ Cu 1 beq @ -0.49524`  layer ABA      

  site Cu1 x =1/3; y =2/3; z = 0.0/Nv;   occ Cu 1 beq @ -0.49524`  layer AAB
  site Cu2 x =1/3; y =2/3; z = (1/3)/Nv; occ Cu 1 beq @ -0.49524`  layer AAB
  site Cu3 x =2/3; y =1/3; z = (2/3)/Nv; occ Cu 1 beq @ -0.49524`  layer AAB     
  
  site Cu1 x =1/3; y =2/3; z = 0.0/Nv;   occ Cu 1 beq @ -0.49524`  layer ABB
  site Cu2 x =2/3; y =1/3; z = (1/3)/Nv; occ Cu 1 beq @ -0.49524`  layer ABB
  site Cu3 x =2/3; y =1/3; z = (2/3)/Nv; occ Cu 1 beq @ -0.49524`  layer ABB      

  site Cu1 x =2/3; y =1/3; z = 0.0/Nv;   occ Cu 1 beq @ -0.49524`  layer BAB
  site Cu2 x =1/3; y =2/3; z = (1/3)/Nv; occ Cu 1 beq @ -0.49524`  layer BAB
  site Cu3 x =2/3; y =1/3; z = (2/3)/Nv; occ Cu 1 beq @ -0.49524`  layer BAB      

  site Cu1 x =2/3; y =1/3; z = 0.0/Nv;   occ Cu 1 beq @ -0.49524`  layer BBA
  site Cu2 x =2/3; y =1/3; z = (1/3)/Nv; occ Cu 1 beq @ -0.49524`  layer BBA
  site Cu3 x =1/3; y =2/3; z = (2/3)/Nv; occ Cu 1 beq @ -0.49524`  layer BBA 
  
  site Cu1 x =2/3; y =1/3; z = 0.0/Nv;   occ Cu 1 beq @ -0.49524`  layer BBB
  site Cu2 x =2/3; y =1/3; z = (1/3)/Nv; occ Cu 1 beq @ -0.49524`  layer BBB
  site Cu3 x =2/3; y =1/3; z = (2/3)/Nv; occ Cu 1 beq @ -0.49524`  layer BBB       
  
	peak_buffer_based_on = Xo; peak_buffer_based_on_tol 0.1
	TCHZ_Peak_Type( , 0.06032, ,-0.01272, ,-0.01003,,0, , 0.12249, , 0.01501) 'Instrument function refined with Nstr=100, Nv=5000
	SF_smooth(@, 1, 1)      
	}
	
	
	} 
 
	'{{{ output data 
	xdd_out temp.xyd append load out_record out_fmt out_eqn 
{ 
    " %11.5f  " = X; 
    " %11.5f  " = Yobs; 
    " %11.5f  " = Ycalc; 
    " %11.5f\n  " = Yobs-Ycalc; 
} 
'}}}

johnsoevans

I've only got very general advice (I haven't had time to look through your INP file in detail), so it's not going to be a very helpful response:

Stacking fault problems are hard and it's really easy to make a mistake with your probabilities and end up with an unexpected sequence.
As you've done, it's useful to check your probabilities on end-members to make sure you have things correct. e.g. a stacking fault sequence for AAAA...., BBBB.... and ABABAB.... should look just like the simple crystallographic models.
If you think there's a fault in the TOPAS calculation you could try using diffax to do the same calculations. Every thing I've checked has been equivalent.
Try one of the save_sequences commands to output your stacking sequence to a file. Is it what you expect?
Sometimes viewing stacking sequences by reading them into excel and auto-colouring the cells can be a quick way to spot problems. It's easier to spot a changing colour sequence than a series of ABBAAABBB.
It won't answer your problem, but reading the diffax manual can be helpful (it's available online if you google). It discusses various examples and how to set probabilities up.

asbridger

Thanks for the help, the general issue I outlined was a TOPAS related one, but it turns out there were also so other errors in my INP file that your advice was very useful for solving.

All sorted now but I figure I should post an explanation of the issue here, in case anyone else runs into the same problem and stumbles across this:

From what I understand (please someone correct me if this is wrong) the problem is basically coming from an internal check in TOPAS when it generates the stacks.

It generates a stack of the requested length, Nv, using a transition matrix generated from the individual transition probabilities.

The software then checks that the number of each transition generated is consistent with the expected number from the probabilities, within a tolerance q (as explained in Figure 1 of https://journals.iucr.org/j/issues/2016/05/00/po5075/index.html). If this check is not met a new layer is added to the stack and the check is repeated. Once this check is met the last Nv layers are used.

This q can be modified with the line: pr_str { B layers_tol } 0.5, setting 0.5 to whatever q value required.

This is great for individual transitions but if you make composite layers which account for some of these transitions within them, sometimes the check is never met.

I figure this is due to the relatively small sampling relative to some of the transition probabilities.

The analogy that makes sense to me is that if you simulate 16 coin tosses one at a time you will get some sequence HHTTHTHTHHTTTHTT...
This will have an expected 4 HH transitions, 4 HT transitions, 4 TT transitions and 4 TT transitions and that will be approximately simulated in your sequence (here it is 2HH, 5HT, 5TH {wrapping around}, 4TT), so check passes

If instead you try and simulate 16 coin tosses but in batches of 4 you can get an entirely equivalent sequence built up from HHTT-HTHT-HHTT-THTT. The difference is that there are many more possible transitions and they all have the same expected value (probability of any transition is 0.5⁴ so 1/16, expected probability of each transition when 4 are calculated is therefore 0.25) but the generated stack matrix has 4 of these transitions with a 1 and everything else 0, this won't be within tolerance so the check will fail.

As both are equivalent sequences they should give the same diffraction pattern so should be equally valid fits to the data.

In the extreme where many, many sets of these 4 batch simulations are produced, the expected probabilities will again match the generated one and pass the checks, but the size of numbers in combinatorics will probably push the number of required samples outside of what is computationally feasible.

If you are also trying to do something similar to me, Alan was kind enough to update my version of TOPAS so this check can be bypassed, but that is not currently in the standard version. It may be sufficient to just set q to be arbitrarily large: pr_str { B layers_tol } 1e9 , in order to also bypass the check, but I'm not sure this will work in all situations.

It does work in the case of my example inp file which I'll attach in case it is of any use

cu-tutorial-triplestack.txt

10kB