-
Notifications
You must be signed in to change notification settings - Fork 0
/
EuphonicDisp2SqwTest.m
106 lines (90 loc) · 4.75 KB
/
EuphonicDisp2SqwTest.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
classdef EuphonicDisp2SqwTest < matlab.mock.TestCase
properties(TestParameter)
% pars, intensity_scale, expected_output_filename
scale_pars = {...
{150., 150., 'expected_cut1d_disp2sqw_eval.sqw'}, ...
{[120., 0.9], 120., 'expected_cut1d_disp2sqw_eval_fscale0p9.sqw'}, ...
{{'frequency_scale', 0.9, 'intensity_scale', 2.5e3}, 2.5e3, 'expected_cut1d_disp2sqw_eval_fscale0p9.sqw'}};
fitting_scale_pars = {...
{200., 200, 'expected_cut1d_disp2sqw_tobyfit.sqw'}, ...
{[1e3, 0.9], 1e3, 'expected_cut1d_disp2sqw_tobyfit_fscale0p9.sqw'}};
end
methods(TestClassSetup)
function setLinuxLibFlags(testCase)
if ~ispc
py.sys.setdlopenflags(int32(10)); % Needed on Linux due to linker issues
end
end
end
methods(Test)
function testQuartzCoherentCrystalDisp2sqweval(testCase, scale_pars)
[pars, iscale, expected_out_fname] = scale_pars{:};
disp('Running testQuartzCoherentCrystalDisp2sqweval...');
disp('Reading sqw...');
ws = read_sqw('quartz/cut1d.sqw')
% Set up simulation
disp('Setting up CoherentCrystal...');
effective_fwhm = 1;
fc = euphonic.ForceConstants.from_castep('quartz/quartz.castep_bin');
euobj = euphonic.CoherentCrystal( ...
fc, 'debye_waller_grid', [6 6 6], 'temperature', 100, ...
'negative_e', true, 'chunk', 1000, 'use_c', true, ...
'dipole_parameter', 0.75);
% Run simulation
disp('Running disp2sqw_eval');
wsim = disp2sqw_eval(ws, @euobj.horace_disp, pars, effective_fwhm);
disp('Reading expected sqw...');
expected_wsim = read_sqw(['quartz', filesep, expected_out_fname]);
disp('Testing result...');
import matlab.unittest.constraints.IsEqualTo
import matlab.unittest.constraints.AbsoluteTolerance
import matlab.unittest.constraints.RelativeTolerance
bounds = AbsoluteTolerance(1e-7*iscale*mean(expected_wsim.data.s)) | RelativeTolerance(1e-5);
testCase.verifyThat(wsim.data.s, ...
IsEqualTo(iscale*expected_wsim.data.s, 'within', bounds));
end
function testQuartzCoherentCrystalDisp2sqwTobyfit(testCase, fitting_scale_pars)
[pars, iscale, expected_out_fname] = fitting_scale_pars{:};
FIXED_SEED = 101;
[rng_state, old_rng_state] = seed_rng(FIXED_SEED);
clean_up = onCleanup(@() rng(old_rng_state));
fprintf('RNG seed: %i\n', rng_state.Seed)
disp('Running testQuartzCoherentCrystalDisp2sqwTobyfit...');
disp('Reading sqw...');
ws = read_sqw('quartz/cut1d.sqw');
% Set up simulation
disp('Setting up CoherentCrystal...');
intrinsic_fwhm = 0.1;
fc = euphonic.ForceConstants.from_castep('quartz/quartz.castep_bin');
euobj = euphonic.CoherentCrystal( ...
fc, 'debye_waller_grid', [6 6 6], 'temperature', 100, ...
'negative_e', true, 'chunk', 10000, 'use_c', true, ...
'dipole_parameter', 0.75);
% Resets the RNG - PR920 somehow reduced the number of rand
% calls before TobyFit was called putting the RNG out of sync
rand(1,20);
% Run simulation with resolution convolution
disp('Running disp2sqw_eval with resolution convolution');
is_crystal = true; xgeom = [0,0,1]; ygeom = [0,1,0]; shape = 'cuboid'; shape_pars = [0.01,0.05,0.01];
sample = IX_sample(is_crystal, xgeom, ygeom, shape, shape_pars);
sample.alatt = [4.9, 4.9, 5.4];
sample.angdeg = [90, 90, 120];
ws = set_sample(ws, sample);
ei = 40; freq = 400; chopper = 'g';
ws = set_instrument(ws, merlin_instrument(ei, freq, chopper));
kk = tobyfit(ws);
kk = kk.set_fun(@disp2sqw, {@euobj.horace_disp, pars, [intrinsic_fwhm]});
kk = kk.set_mc_points(3);
wsim = kk.simulate('fore');
disp('Reading expected sqw...');
expected_wsim = read_sqw(['quartz', filesep, expected_out_fname]);
disp('Testing result...');
import matlab.unittest.constraints.IsEqualTo
import matlab.unittest.constraints.AbsoluteTolerance
import matlab.unittest.constraints.RelativeTolerance
bounds = AbsoluteTolerance(1e-7*iscale*mean(expected_wsim.data.s)) | RelativeTolerance(2e-4);
testCase.verifyThat(wsim.data.s, ...
IsEqualTo(iscale*expected_wsim.data.s, 'within', bounds));
end
end
end