Error in PDF refinement beqs for equal form factors

sebastianm

Hi everyone,

I have a specific problem when using a user defined macro (found in the Topas Academic wiki) for displacement parameters in PDF refinements. Here is a minimal example:

xdd PS3SC_SDD360_60s-00103.xy

pdf_data

convolute_Qmax_Sinc( !Qmax, 19.3)
dQ_damp (,Qdamp)

prm delt1 1.66356` min 0.01
prm !delt2 0.0
prm !Qbroad 0.02816 min 0.001
prm !Qdamp 0.02377 min 0.01

str
space_group "182"
Hexagonal( @ 5.781930`, @ 12.624767`)
spherical_damping(r, 413.73481`)

scale @ 0.546491048`
site Mn1 x =1/3; y =2/3; z 0.25 occ Mn+2 1.2 Beqpdfgui(bmn, 1.17378`)
site Nb1 x 0 y 0 z 0 occ Nb+5 1 Beqpdfgui(bnb, 0.31098`)
site Nb2 x =1/3; y =2/3; z 0.99768 occ Nb 1 Beqpdfgui(bnb, 0.31098`)
site S1 x 0.33422 y 0.00043 z 0.37600 occ S 1 Beqpdfgui(bs, 0.63027`)

macro Beqpdfgui(c, v) ' ##### Defines the r-dependent ADPs as a function of mean squared displacements, linear #####
{ ' ##### and quadratic sharpening due to correlated motion, and angular dependence of Q resolution #####
#m_argu c
If_Prm_Eqn_Rpt(c, v, min 1e-6 max 5 val_on_continue = Val Rand(.1, 2); )
beq = CeV(c, v) Sqrt(1- Min((delt1/X + delt2/X^2 - Qbroad^2 X^2),1));
}

macro dQ_damp(dQ,dQv) ' ##### Multiplies PDF by a Gaussian envelope function to simulate finite Q resolution effects #####
{
#m_argu dQ 'Gaussian damping factor dur to instrumental resolution of S(Q) data
If_Prm_Eqn_Rpt(dQ, dQv, min 0.001 max 0.2, del 0.0001)
scale_phase_X = Exp(-0.5 X^2 (CeV(dQ,dQv))^2);
}

There are two distinct sites for Nb. The input file works as posted above. However when using 'occ Nb 1' on both Nb-sites (instead of Nb/Nb+5) the exception is thrown.

Does anybody have an idea how to solve this?

All best
Sebastian

rowlesmr

What is the exact text of the exception? Is it an actual exception in a popup box, or is it an internal error in topas that prints to the output screen?

Can you post your data too, so we can run the input file?

Matthew

johnsoevans

Sebastian,

Two comments.

1. I can't see any reason why having Nb and Nb5+ on the different sties would matter. I tried your file in topas v7 and it seemed to run ok against a dummy dataset with Nb5+ on both sites. However I also got an exception in v6.

2. Where did your beq macro come from? It's not one of the ones in Phil Chater's pdf.inc at: https://github.com/pachater/topas/blob/master/pdf.inc. Can you get your file to run with one of the standard beq macros?

If it ends up being a v6 vs v7 issue you'll probably have to contact Alan for an answer.

John

sebastianm

Hi,

you can find the data file attached. The error is a popup box saying 'Exception caught'.

@ John: Thanks for trying the file on both Topas versions. In this simple case the work around with Nb/Nb5+ is fine. However, having five sites with the same elements but giving each a different scattering factor is kinda tedious.
I got the macro from a workshop tutorial, it seems to be a simplified version of the macro beq_PDFfit2 by Phil Chater. That macro throws the same error when having two identical scattering factors.

I just found out one more thing: The problem can be solved when using the same scattering factor (Nb on both sites) but a different occupancy (1.0 and 0.99). That makes it easier to work with, still I wonder what causes the problem. I will ask Alan about that.

Sebastian

rowlesmr

There is something going on.

I get no issues in v7, but there are in v6 - That's an Alan thing.

In v6, if you remove all the BEQ macros, and have both sites as "Nb 1", it runs. (the model is wrong, but it runs).

I don't know enough about PDF to comment further.

.

Also, why does Mn occ == 1.2?

alancoelho

Hi Sebastian

I removed the convolute_Qmax_Sinc and spherical_damping macros (as they don't exist in the INP file) and then ran your INP file in TOPAS-Academic Version 6 and didn't get an error; ie.

TOPAS-64 Version 6 (c) 1992-2016 Alan A. Coelho
Maximum number of threads 8
Constant x-axis step size detected for xy/xye file
Internally changing to constant x-axis step size for efficiency
Loading xyz's for 182 from file C:\c\t5-ta6\ta-web2\sg\182.sg
Number of independent parameters : 7
Space group 182 is not centrosymmetric
Number of equivalent positions 12
G(r), approximate number of interactions: 1774560
0 Time 0.07 Rwp 22.790 0.000 MC 0.00 0
Constant x-axis step size detected for xy/xye file
Internally changing to constant x-axis step size for efficiency
1 Time 0.14 Rwp 21.532 -1.259 MC 0.16 1
2 Time 0.17 Rwp 21.374 -0.157 MC 0.05 1
3 Time 0.23 Rwp 21.374 -0.001 MC 29.84 3
--- 0.238 seconds ---
File C:\tmp\sm\aac.out updated
with parameters from last iteration

Send me the two macros and I will run the INP file again with then included.
cheers
alan

rowlesmr

Hi Alan

Those two macros are in pdf.inc

Did you try with both Nb sites with "occ Nb 1"?

Matthew

alancoelho

Thanks Matthew

OK, there's a problem is that the beq equation is being evaluated at X=0. This is a bug that I think it has been fixed in V7.

Anyway, the work around is as follows:

macro X_ { (X + 1e-5) }

macro Beqpdfgui(c, v)
{                    
	#m_argu c
	If_Prm_Eqn_Rpt(c, v, min 1e-6 max 5 val_on_continue = Val Rand(.1, 2); )
	beq = CeV(c, v) Sqrt(1- Min((delt1/X_ + delt2/X_^2 - Qbroad^2 X_^2),1));
}

sebastianm

Dear Alan,

thanks a lot for having a look and fixing the problem.

brycemullens

Hi all,

I'm not too sure if this thread is still active, but I have a similar issue when fitting XPDF data that’s a bit more odd.

I'm trying to do a PDF fit using Diamond data. There are two possible space groups I'm trying to fit: I41/a and I-4. The I-4 seems to have a similar problem mentioned in this thread, with the same error message.

	#ifdef I_4_double_scheelite
	str
		scale =scale_I_4;
		a =a_lat_I_4; 
		b=Get(a);  
		c =c_lat_I_4;
		al 90.0     
		be 90.0     
		ga 90.0     
		space_group "I -4"

		site Bi1	x = 0;		y = 1/2;		z = 1/4;		occ Bi Bi_1_occ 0.5			beq_A  
		site Bi2	x = 1/2;		y = 1/2;		z = 0;		occ Bi = 1 - Bi_1_occ;: 0.5		beq_A  
		site Na1	x = 0;		y = 1/2;		z = 1/4;		occ Na = 1 - Bi_1_occ;: 0.5	beq_A  
		site Na2	x = 1/2;		y = 1/2;		z = 0;		occ Na = Bi_1_occ;: 0.5		beq_A  
		site W1	x = 0;		y = 0;		z = 0;		occ Mo 1.0					beq_B 
		site W2	x = 1/2;		y = 0;		z = 1/4;		occ Mo 1.0					beq_B 
		site O1	x @ 0.3472	y @ 0.2318     z @ 0.1646		occ O  1.0					beq_O 
		site O2	x @ 0.8456	y @ 0.2461     z @ 0.0846		occ O  1.0					beq_O 

	macro fit_type {_tetra_I_4_double}
	#endif

I’m using the exact same .inp file, and I noticed it runs differently depending on which version of TOPAS I use. It also only seems to be a problem with XPDF data.
Topas V6, academic, XPDF, I-4: works
Topas V6, commercial, XPDF, I-4: does not work
Topas V6, commercial, NPDF, I-4: works

As already mentioned in this thread, commenting out 'convolute_Qmax_Sinc' makes it work.

Any ideas what could be going on here? Thanks!

rowlesmr

Hi Bryce

There could be a division by zero error. The version of the convolute_Qmax_Sinc macro that comes with TOPAS7 is at the bottom of this post; can you check the versions of TOPAS that you have to see if their macros are the same?

If they are the same, then try the same fix as previously given, and define your own macro and use that:

macro & X_ { X + 0.00001 }

macro bryce_convolute_Qmax_Sinc(Qmax,Qmaxv)
{
	#m_argu Qmax
	If_Prm_Eqn_Rpt(Qmax, Qmaxv, min 1 max 100, del 0.01)
	local #m_unique conv_max = (5 CeV(Qmax,Qmaxv) - Mod(5 CeV(Qmax,Qmaxv),1))/5 2 Pi / CeV(Qmax,Qmaxv);
	pdf_convolute = If(Abs(X_) > Yobs_dx_at(Xo),(Sin(CeV(Qmax,Qmaxv) X_)/(X_)),CeV(Qmax,Qmaxv));
		min_X = Min(-Xo,-conv_max) ;
		max_X = Max( Xo, conv_max) ;
}

Maybe.

Give it a go, and see what happens. If you are able to share the entire inp file and data, I can have a go on v7 when I get home.

Matthew

---------
Just for the record:

macro convolute_Qmax_Sinc(Qmax,Qmaxv)
{
	'Convolution to include Sinc function due to termeination of S(Q) at finite Qmax during Fourier transform
	'Use in combination with start_X 1 to avoid very low-r data 
	'Ref: J. S. Chung and M. F. Thorpe, Phys. Rev. B, 55(3), 1997, 1545
	#m_argu Qmax
	If_Prm_Eqn_Rpt(Qmax, Qmaxv, min 1 max 100, del 0.01)
	local #m_unique conv_max = (5 CeV(Qmax,Qmaxv) - Mod(5 CeV(Qmax,Qmaxv),1))/5 2 Pi / CeV(Qmax,Qmaxv);
	pdf_convolute = If(Abs(X) > Yobs_dx_at(Xo),(Sin(CeV(Qmax,Qmaxv) X)/(X)),CeV(Qmax,Qmaxv));
		min_X = Min(-Xo,-conv_max) ;
		max_X = Max( Xo, conv_max) ;
}

brycemullens

Hi Matthew,

Thanks for the reply. Unfortunately, I didn't have much luck trying to define my own macro.

I'd be interested to see if it runs on V7. I've attached my .inp file, as well as the datafile. If you wouldn't mind giving trying to run it, that would be great!

Cheers,

Bryce

brycemullens

...and here's the datafile itself.

rowlesmr

Hi Bryce

It all appears to work on mine (6 & 7). I also realised that my Topas 6 is academic, so I can't test that pathway.

You missed a bit in your macros: `(Sin(CeV(Qmax,Qmaxv) X)/(X_))` you needed to make the X in the denominator the X_ variable to avoid the (possible) divide by zero error.

You could also try trplaceing all the X variables in the pdf_convolute equation to be X_s.

Matthew