test_interface/
test_interface.rs

1#![allow(clippy::excessive_precision)]
2
3use approx::assert_abs_diff_eq;
4use dftd4::prelude::*;
5
6fn test_rational_damping_noargs() {
7    // Check constructor of damping parameters for insufficient arguments
8    assert!(DFTD4RationalDampingParamBuilder::default().build().is_err());
9
10    let builder = DFTD4RationalDampingParamBuilder::default().a1(0.4).a2(5.0);
11    assert!(builder.build().is_err_and(|e| match &e {
12        DFTD4Error::BuilderError(f) => f.field_name() == "s8",
13        _ => false,
14    }));
15
16    let builder = DFTD4RationalDampingParamBuilder::default().s8(1.0).a2(5.0);
17    assert!(builder.build().is_err_and(|e| match &e {
18        DFTD4Error::BuilderError(f) => f.field_name() == "a1",
19        _ => false,
20    }));
21
22    let builder = DFTD4RationalDampingParamBuilder::default().s8(1.0).a1(0.4);
23    assert!(builder.build().is_err_and(|e| match &e {
24        DFTD4Error::BuilderError(f) => f.field_name() == "a2",
25        _ => false,
26    }));
27}
28
29fn test_structure() {
30    // Check if the molecular structure data is working as expected
31    let numbers = vec![6, 7, 6, 7, 6, 6, 6, 8, 7, 6, 8, 7, 6, 6, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1];
32    #[rustfmt::skip]
33    let positions = vec![
34         2.02799738646442,  0.09231312124713, -0.14310895950963,
35         4.75011007621000,  0.02373496014051, -0.14324124033844,
36         6.33434307654413,  2.07098865582721, -0.14235306905930,
37         8.72860718071825,  1.38002919517619, -0.14265542523943,
38         8.65318821103610, -1.19324866489847, -0.14231527453678,
39         6.23857175648671, -2.08353643730276, -0.14218299370797,
40         5.63266886875962, -4.69950321056008, -0.13940509630299,
41         3.44931709749015, -5.48092386085491, -0.14318454855466,
42         7.77508917214346, -6.24427872938674, -0.13107140408805,
43        10.30229550927022, -5.39739796609292, -0.13672168520430,
44        12.07410272485492, -6.91573621641911, -0.13666499342053,
45        10.70038521493902, -2.79078533715849, -0.14148379504141,
46        13.24597858727017, -1.76969072232377, -0.14218299370797,
47         7.40891694074004, -8.95905928176407, -0.11636933482904,
48         1.38702118184179,  2.05575746325296, -0.14178615122154,
49         1.34622199478497, -0.86356704498496,  1.55590600570783,
50         1.34624089204623, -0.86133716815647, -1.84340893849267,
51         5.65596919189118,  4.00172183859480, -0.14131371969009,
52        14.67430918222276, -3.26230980007732, -0.14344911021228,
53        13.50897177220290, -0.60815166181684,  1.54898960808727,
54        13.50780014200488, -0.60614855212345, -1.83214617078268,
55         5.41408424778406, -9.49239668625902, -0.11022772492007,
56         8.31919801555568, -9.74947502841788,  1.56539243085954,
57         8.31511620712388, -9.76854236502758, -1.79108242206824,
58    ];
59
60    // Constructor should raise an error for nuclear fusion input
61    assert!(DFTD4Structure::new_f(&numbers, &vec![0.0; 72], None, None, None)
62        .is_err_and(|e| e.to_string().contains("Too close interatomic distances found")));
63
64    // The Rust struct should protect from garbage input like this
65    assert!(DFTD4Structure::new_f(&[1, 1, 1], &positions, None, None, None)
66        .is_err_and(|e| e.to_string().contains("Invalid dimension for positions")));
67
68    // Also check for sane coordinate input
69    assert!(DFTD4Structure::new_f(&numbers, &[0.0; 7], None, None, None)
70        .is_err_and(|e| e.to_string().contains("Invalid dimension for positions")));
71
72    // Construct real molecule
73    let mut mol = DFTD4Structure::new(&numbers, &positions, None, None, None);
74
75    // Try to update a structure with mismatched coordinates
76    assert!(mol
77        .update_f(&[0.0; 7], None)
78        .is_err_and(|e| e.to_string().contains("Invalid dimension for positions")));
79
80    // Try to add a mismatched lattice
81    assert!(mol
82        .update_f(&positions, Some(&[0.0; 7]))
83        .is_err_and(|e| e.to_string().contains("Invalid dimension for lattice")));
84
85    // Try to update a structure with nuclear fusion coordinates
86    assert!(mol
87        .update_f(&[0.0; 72], None)
88        .is_err_and(|e| e.to_string().contains("Too close interatomic distances found"),));
89}
90
91fn test_blypd4() {
92    // Use BLYP-D4 for a mindless molecule
93    let thr = 1.0e-7;
94
95    let numbers = vec![1, 1, 6, 5, 1, 15, 8, 17, 13, 15, 5, 1, 9, 15, 1, 15];
96    #[rustfmt::skip]
97    let positions = vec![
98         2.79274810283778,  3.82998228828316, -2.79287054959216,
99        -1.43447454186833,  0.43418729987882,  5.53854345129809,
100        -3.26268343665218, -2.50644032426151, -1.56631149351046,
101         2.14548759959147, -0.88798018953965, -2.24592534506187,
102        -4.30233097423181, -3.93631518670031, -0.48930754109119,
103         0.06107643564880, -3.82467931731366, -2.22333344469482,
104         0.41168550401858,  0.58105573172764,  5.56854609916143,
105         4.41363836635653,  3.92515871809283,  2.57961724984000,
106         1.33707758998700,  1.40194471661647,  1.97530004949523,
107         3.08342709834868,  1.72520024666801, -4.42666116106828,
108        -3.02346932078505,  0.04438199934191, -0.27636197425010,
109         1.11508390868455, -0.97617412809198,  6.25462847718180,
110         0.61938955433011,  2.17903547389232, -6.21279842416963,
111        -2.67491681346835,  3.00175899761859,  1.05038813614845,
112        -4.13181080289514, -2.34226739863660, -3.44356159392859,
113         2.85007173009739, -2.64884892757600,  0.71010806424206,
114    ];
115
116    let model = DFTD4Model::new(&numbers, &positions, None, None, None);
117
118    let param = DFTD4Param::load_rational_damping("blyp", false);
119    let res = model.get_dispersion(&param, false);
120
121    assert_abs_diff_eq!(res.energy, -0.06991716314879085, epsilon = thr);
122
123    let res = model.get_dispersion(&param, true);
124
125    assert_abs_diff_eq!(res.energy, -0.06991716314879085, epsilon = thr);
126}
127
128#[cfg(feature = "d4s")]
129fn test_tpssd4s() {
130    // Use TPSS-D4S for a mindless molecule
131    let thr = 1.0e-7;
132
133    let numbers = vec![1, 1, 6, 5, 1, 15, 8, 17, 13, 15, 5, 1, 9, 15, 1, 15];
134    #[rustfmt::skip]
135    let positions = vec![
136         2.79274810283778,  3.82998228828316, -2.79287054959216,
137        -1.43447454186833,  0.43418729987882,  5.53854345129809,
138        -3.26268343665218, -2.50644032426151, -1.56631149351046,
139         2.14548759959147, -0.88798018953965, -2.24592534506187,
140        -4.30233097423181, -3.93631518670031, -0.48930754109119,
141         0.06107643564880, -3.82467931731366, -2.22333344469482,
142         0.41168550401858,  0.58105573172764,  5.56854609916143,
143         4.41363836635653,  3.92515871809283,  2.57961724984000,
144         1.33707758998700,  1.40194471661647,  1.97530004949523,
145         3.08342709834868,  1.72520024666801, -4.42666116106828,
146        -3.02346932078505,  0.04438199934191, -0.27636197425010,
147         1.11508390868455, -0.97617412809198,  6.25462847718180,
148         0.61938955433011,  2.17903547389232, -6.21279842416963,
149        -2.67491681346835,  3.00175899761859,  1.05038813614845,
150        -4.13181080289514, -2.34226739863660, -3.44356159392859,
151         2.85007173009739, -2.64884892757600,  0.71010806424206,
152    ];
153
154    let model = DFTD4Model::new_d4s(&numbers, &positions, None, None, None);
155
156    let param = DFTD4Param::load_rational_damping("tpss", false);
157    let res = model.get_dispersion(&param, false);
158
159    assert_abs_diff_eq!(res.energy, -0.046233140236052253, epsilon = thr);
160
161    let res = model.get_dispersion(&param, true);
162
163    assert_abs_diff_eq!(res.energy, -0.046233140236052253, epsilon = thr);
164}
165
166fn test_pbed4() {
167    // Use PBE-D4 for a mindless molecule
168    let thr = 1.0e-7;
169
170    let numbers = vec![1, 9, 15, 13, 1, 1, 13, 5, 3, 15, 8, 1, 1, 5, 16, 1];
171    #[rustfmt::skip]
172    let positions = vec![
173        -2.14132037405479, -1.34402701877044, -2.32492500904728,
174         4.46671289205392, -2.04800110524830,  0.44422406067087,
175        -4.92212517643478, -1.73734240529793,  0.96890323821450,
176        -1.30966093045696, -0.52977363497805,  3.44453452239668,
177        -4.34208759006189, -4.30470270977329,  0.39887431726215,
178         0.61788392767516,  2.62484136683297, -3.28228926932647,
179         4.23562873444840, -1.68839322682951, -3.53824299552792,
180         2.23130060612446,  1.93579813100155, -1.80384647554323,
181        -2.32285463652832,  2.90603947535842, -1.39684847191937,
182         2.34557941578250,  2.86074312333371,  1.82827238641666,
183        -3.66431367659153, -0.42910188232667, -1.81957402856634,
184        -0.34927881505446, -1.75988134003940,  5.98017466326572,
185         0.29500802281217, -2.00226104143537,  0.53023447931897,
186         2.10449364205058, -0.56741404446633,  0.30975625014335,
187        -1.59355304432499,  3.69176153150419,  2.87878226787916,
188         4.34858700256050,  2.39171478113440, -2.61802993563738,
189    ];
190
191    let model = DFTD4Model::new(&numbers, &positions, None, None, None);
192
193    let param = DFTD4Param::load_rational_damping("pbe", false);
194    let res = model.get_dispersion(&param, false);
195
196    assert_abs_diff_eq!(res.energy, -0.028415184156428127, epsilon = thr);
197
198    let res = model.get_dispersion(&param, true);
199
200    assert_abs_diff_eq!(res.energy, -0.028415184156428127, epsilon = thr);
201}
202
203#[cfg(feature = "d4s")]
204fn test_r2scan3c() {
205    // Use r2SCAN-3c for a mindless molecule
206    let thr = 1.0e-8;
207
208    let numbers = vec![1, 9, 15, 13, 1, 1, 13, 5, 3, 15, 8, 1, 1, 5, 16, 1];
209    #[rustfmt::skip]
210    let positions = vec![
211        -2.14132037405479, -1.34402701877044, -2.32492500904728,
212         4.46671289205392, -2.04800110524830,  0.44422406067087,
213        -4.92212517643478, -1.73734240529793,  0.96890323821450,
214        -1.30966093045696, -0.52977363497805,  3.44453452239668,
215        -4.34208759006189, -4.30470270977329,  0.39887431726215,
216         0.61788392767516,  2.62484136683297, -3.28228926932647,
217         4.23562873444840, -1.68839322682951, -3.53824299552792,
218         2.23130060612446,  1.93579813100155, -1.80384647554323,
219        -2.32285463652832,  2.90603947535842, -1.39684847191937,
220         2.34557941578250,  2.86074312333371,  1.82827238641666,
221        -3.66431367659153, -0.42910188232667, -1.81957402856634,
222        -0.34927881505446, -1.75988134003940,  5.98017466326572,
223         0.29500802281217, -2.00226104143537,  0.53023447931897,
224         2.10449364205058, -0.56741404446633,  0.30975625014335,
225        -1.59355304432499,  3.69176153150419,  2.87878226787916,
226         4.34858700256050,  2.39171478113440, -2.61802993563738,
227    ];
228
229    let model = DFTD4Model::custom_d4s(&numbers, &positions, 2.0, 1.0, None, None, None);
230
231    let param = DFTD4RationalDampingParamBuilder::default()
232        .s8(0.0)
233        .a1(0.42)
234        .a2(5.65)
235        .s9(2.0)
236        .build()
237        .unwrap()
238        .new_param();
239    let res = model.get_dispersion(&param, false);
240
241    assert_abs_diff_eq!(res.energy, -0.008016697276824889, epsilon = thr);
242
243    let res = model.get_dispersion(&param, true);
244
245    assert_abs_diff_eq!(res.energy, -0.008016697276824889, epsilon = thr);
246}
247
248fn test_pair_resolved() {
249    // Calculate pairwise resolved dispersion energy for a molecule
250    let thr = 1.0e-8;
251
252    let numbers = vec![16, 16, 16, 16, 16, 16, 16, 16];
253    #[rustfmt::skip]
254    let positions = vec![
255        -4.15128787379191,  1.71951973863958, -0.93066267097296,
256        -4.15128787379191, -1.71951973863958,  0.93066267097296,
257        -1.71951973863958, -4.15128787379191, -0.93066267097296,
258         1.71951973863958, -4.15128787379191,  0.93066267097296,
259         4.15128787379191, -1.71951973863958, -0.93066267097296,
260         4.15128787379191,  1.71951973863958,  0.93066267097296,
261         1.71951973863958,  4.15128787379191, -0.93066267097296,
262        -1.71951973863958,  4.15128787379191,  0.93066267097296,
263    ];
264    #[rustfmt::skip]
265    let pair_disp2 = vec![
266        -0.00000000e-0, -5.80599854e-4, -4.74689854e-4, -2.11149449e-4, -1.63163128e-4, -2.11149449e-4, -4.74689854e-4, -5.80599854e-4,
267        -5.80599854e-4, -0.00000000e-0, -5.80599854e-4, -4.74689854e-4, -2.11149449e-4, -1.63163128e-4, -2.11149449e-4, -4.74689854e-4,
268        -4.74689854e-4, -5.80599854e-4, -0.00000000e-0, -5.80599854e-4, -4.74689854e-4, -2.11149449e-4, -1.63163128e-4, -2.11149449e-4,
269        -2.11149449e-4, -4.74689854e-4, -5.80599854e-4, -0.00000000e-0, -5.80599854e-4, -4.74689854e-4, -2.11149449e-4, -1.63163128e-4,
270        -1.63163128e-4, -2.11149449e-4, -4.74689854e-4, -5.80599854e-4, -0.00000000e-0, -5.80599854e-4, -4.74689854e-4, -2.11149449e-4,
271        -2.11149449e-4, -1.63163128e-4, -2.11149449e-4, -4.74689854e-4, -5.80599854e-4, -0.00000000e-0, -5.80599854e-4, -4.74689854e-4,
272        -4.74689854e-4, -2.11149449e-4, -1.63163128e-4, -2.11149449e-4, -4.74689854e-4, -5.80599854e-4, -0.00000000e-0, -5.80599854e-4,
273        -5.80599854e-4, -4.74689854e-4, -2.11149449e-4, -1.63163128e-4, -2.11149449e-4, -4.74689854e-4, -5.80599854e-4, -0.00000000e-0,
274    ];
275    #[rustfmt::skip]
276    let pair_disp3 = vec![
277        0.00000000e-0, 3.39353850e-7, 8.74462839e-7, 1.17634100e-6, 9.86937310e-7, 1.17634100e-6, 8.74462839e-7, 3.39353850e-7,
278        3.39353850e-7, 0.00000000e-0, 3.39353850e-7, 8.74462839e-7, 1.17634100e-6, 9.86937310e-7, 1.17634100e-6, 8.74462839e-7,
279        8.74462839e-7, 3.39353850e-7, 0.00000000e-0, 3.39353850e-7, 8.74462839e-7, 1.17634100e-6, 9.86937310e-7, 1.17634100e-6,
280        1.17634100e-6, 8.74462839e-7, 3.39353850e-7, 0.00000000e-0, 3.39353850e-7, 8.74462839e-7, 1.17634100e-6, 9.86937310e-7,
281        9.86937310e-7, 1.17634100e-6, 8.74462839e-7, 3.39353850e-7, 0.00000000e-0, 3.39353850e-7, 8.74462839e-7, 1.17634100e-6,
282        1.17634100e-6, 9.86937310e-7, 1.17634100e-6, 8.74462839e-7, 3.39353850e-7, 0.00000000e-0, 3.39353850e-7, 8.74462839e-7,
283        8.74462839e-7, 1.17634100e-6, 9.86937310e-7, 1.17634100e-6, 8.74462839e-7, 3.39353850e-7, 0.00000000e-0, 3.39353850e-7,
284        3.39353850e-7, 8.74462839e-7, 1.17634100e-6, 9.86937310e-7, 1.17634100e-6, 8.74462839e-7, 3.39353850e-7, 0.00000000e-0,
285    ];
286
287    let model = DFTD4Model::new(&numbers, &positions, None, None, None);
288
289    let param = DFTD4RationalDampingParamBuilder::default()
290        .s8(1.20065498)
291        .a1(0.40085597)
292        .a2(5.02928789)
293        .build()
294        .unwrap()
295        .new_param();
296    let res = model.get_pairwise_dispersion(&param);
297
298    for (a, b) in res.pair_energy2.iter().zip(pair_disp2.iter()) {
299        assert_abs_diff_eq!(a, b, epsilon = thr);
300    }
301
302    for (a, b) in res.pair_energy3.iter().zip(pair_disp3.iter()) {
303        assert_abs_diff_eq!(a, b, epsilon = thr);
304    }
305}
306
307fn test_properties() {
308    // Calculate dispersion related properties for a molecule
309    let thr = 1.0e-7;
310
311    let numbers = vec![7, 7, 7, 7, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1];
312    #[rustfmt::skip]
313    let positions = vec![
314        1.97420621099560, 1.97415497783241, 1.97424596974304,
315        6.82182427659395, 2.87346383480995, 7.72099517560089,
316        7.72104957181201, 6.82177051521773, 2.87336561318016,
317        2.87343220660781, 7.72108897828386, 6.82187093171878,
318        3.51863272100286, 2.63865333484548, 1.00652979981286,
319        2.63877594964754, 1.00647313885594, 3.51882748086447,
320        1.00639728563189, 3.51850454450845, 2.63869202592387,
321        8.36624975982697, 2.20896711017229, 8.68870955681018,
322        7.48639684558259, 3.84114715917956, 6.17640982573725,
323        5.85401675167715, 1.32911569888797, 7.05654606696031,
324        7.05646299938990, 5.85409590282274, 1.32879923864813,
325        8.68882633853582, 8.36611541129785, 2.20894120662207,
326        3.84121223226912, 6.17673669892998, 7.48629723649480,
327        1.32897854262127, 7.05658604099926, 5.85414031368096,
328        2.20884896069885, 8.68875820985799, 8.36643568423387,
329        6.17659142004652, 7.48627051643848, 3.84109594690835,
330    ];
331    #[rustfmt::skip]
332    let lattice = vec![
333        9.69523775911749, 0.0, 0.0,
334        0.0, 9.69523775911749, 0.0,
335        0.0, 0.0, 9.69523775911749,
336    ];
337    #[rustfmt::skip]
338    let cn = [
339        2.5775156287150218,
340        2.5775155620078536,
341        2.5775157938667150,
342        2.5775157704485387,
343        0.8591731475439074,
344        0.8591680526841657,
345        0.8591744284869478,
346        0.8591732359038715,
347        0.8591678667769283,
348        0.8591744593270527,
349        0.8591684383407867,
350        0.8591754863625011,
351        0.8591751690053771,
352        0.8591719636497469,
353        0.8591686377934131,
354        0.8591718691634261,
355    ];
356    #[rustfmt::skip]
357    let charges = [
358        -0.86974543285813199,
359        -0.86974501326130316,
360        -0.86974658178316333,
361        -0.86974643466661739,
362         0.28992010784471012,
363         0.28989847135893759,
364         0.28992738637932841,
365         0.28992042841052362,
366         0.28989746859382759,
367         0.28992753735979965,
368         0.28990032127437559,
369         0.28993196014690176,
370         0.28993044356403513,
371         0.28991421277784613,
372         0.28990119160907674,
373         0.28991393324985348,
374    ];
375    #[rustfmt::skip]
376    let alpha = [
377        9.9853045768095097,
378        9.9853025181287016,
379        9.9853102162086920,
380        9.9853094944097140,
381        1.3513315505023076,
382        1.3513939255952379,
383        1.3513106323935196,
384        1.3513306244831862,
385        1.3513968091133417,
386        1.3513101978580895,
387        1.3513885997538553,
388        1.3512974505020905,
389        1.3513018164322617,
390        1.3513485145877058,
391        1.3513860914395939,
392        1.3513493246534483,
393    ];
394
395    let model = DFTD4Model::new(&numbers, &positions, None, Some(&lattice), None);
396
397    let res = model.get_properties();
398
399    for (a, b) in res.cn.iter().zip(cn.iter()) {
400        assert_abs_diff_eq!(a, b, epsilon = thr);
401    }
402
403    for (a, b) in res.charges.iter().zip(charges.iter()) {
404        assert_abs_diff_eq!(a, b, epsilon = thr);
405    }
406
407    for (a, b) in res.alpha.iter().zip(alpha.iter()) {
408        assert_abs_diff_eq!(a, b, epsilon = thr);
409    }
410}
411
412#[test]
413fn test() {
414    test_rational_damping_noargs();
415    test_structure();
416    test_blypd4();
417    #[cfg(feature = "d4s")]
418    test_tpssd4s();
419    test_pbed4();
420    #[cfg(feature = "d4s")]
421    test_r2scan3c();
422    test_pair_resolved();
423    test_properties();
424}
425
426fn main() {
427    test_rational_damping_noargs();
428    test_structure();
429    test_blypd4();
430    #[cfg(feature = "d4s")]
431    test_tpssd4s();
432    test_pbed4();
433    #[cfg(feature = "d4s")]
434    test_r2scan3c();
435    test_pair_resolved();
436    test_properties();
437}