regtests.py 43.3 KB
Newer Older
1
2
3
4
5
6
7
8
import os
import unittest
import logging
import numpy as np
from cp2kparser import CP2KParser
from nomadcore.unit_conversion.unit_conversion import convert_unit


9
10
def get_result(folder, metaname=None):
    """Get the results from the calculation in the given folder. By default goes through different
11
12
13

    Args:
        folder: The folder relative to the directory of this script where the
14
            parsed calculation resides.
15
16
        metaname(str): Optional quantity to return. If not specified, returns
            the full dictionary of results.
17
18
    """
    dirname = os.path.dirname(__file__)
19
    filename = os.path.join("cp2k_{}".format(VERSION), dirname, folder, "unittest.out")
20
21
    parser = CP2KParser(None, debug=True, log_level=logging.CRITICAL)
    results = parser.parse(filename)
22

23
24
    if metaname is None:
        return results
25
    else:
26
        return results[metaname]
27
28


29
class TestErrors(unittest.TestCase):
30
    """Test error situations which may occur during the parsing.
31
32
    """
    def test_no_file(self):
33
34
        """File is not no present.
        """
35
        with self.assertRaises(ValueError):
36
            get_result("errors/no_file", "XC_functional")
37
38

    def test_invalid_file(self):
39
40
41
42
        """Main file is invalid.
        """
        with self.assertRaises(RuntimeError):
            get_result("errors/invalid_file", "XC_functional")
43
44

    def test_invalid_run_type(self):
45
46
47
48
49
        """Unrecognized run type.
        """
        with self.assertRaises(KeyError):
            get_result("errors/invalid_run_type", "XC_functional")

50
51
52
53
54
55
    def test_short_main_file(self):
        """Test what happens if main file is very small.
        """
        with self.assertRaises(RuntimeError):
            get_result("errors/short_main_file", "XC_functional")

56

57
58
59
class TestUnknownInput(unittest.TestCase):
    """Tests for cases where unknown information is encountered in the parsing.
    """
60
    def test_unknown_version(self):
61
62
        """Test how a new version is handled.
        """
63
64
        get_result("errors/unknown_version", "XC_functional")

65
    def test_unknown_input_keyword(self):
66
67
        """Test how an unknown input keyword is handled.
        """
68
69
70
        get_result("errors/unknown_input_keyword", "XC_functional")

    def test_unknown_input_section(self):
71
72
        """Test unknown input file section.
        """
73
74
75
        get_result("errors/unknown_input_section", "XC_functional")

    def test_unknown_input_section_parameter(self):
76
77
        """
        """
78
79
        get_result("errors/unknown_input_section_parameter", "XC_functional")

80

81
class TestXCFunctional(unittest.TestCase):
82
83
    """Tests that the XC functionals can be properly parsed.
    """
84
85
    def test_pade(self):
        xc = get_result("XC_functional/pade", "XC_functional")
86
        self.assertEqual(xc, "1*LDA_XC_TETER93")
87
88
89

    def test_lda(self):
        xc = get_result("XC_functional/lda", "XC_functional")
90
        self.assertEqual(xc, "1*LDA_XC_TETER93")
91
92
93

    def test_blyp(self):
        xc = get_result("XC_functional/blyp", "XC_functional")
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
        self.assertEqual(xc, "1*GGA_C_LYP+1*GGA_X_B88")

    def test_b3lyp(self):
        xc = get_result("XC_functional/b3lyp", "XC_functional")
        self.assertEqual(xc, "1*HYB_GGA_XC_B3LYP")

    def test_olyp(self):
        xc = get_result("XC_functional/olyp", "XC_functional")
        self.assertEqual(xc, "1*GGA_C_LYP+1*GGA_X_OPTX")

    def test_hcth120(self):
        xc = get_result("XC_functional/hcth120", "XC_functional")
        self.assertEqual(xc, "1*GGA_XC_HCTH_120")

    def test_pbe0(self):
        xc = get_result("XC_functional/pbe0", "XC_functional")
        self.assertEqual(xc, "1*HYB_GGA_XC_PBEH")

    def test_pbe(self):
        xc = get_result("XC_functional/pbe", "XC_functional")
        self.assertEqual(xc, "1*GGA_C_PBE+1*GGA_X_PBE")
115
116


117
118
119
120
121
122
123
124
125
126
127
128
129
class TestSCFConvergence(unittest.TestCase):
    """Tests whether the convergence status and number of SCF step can be
    parsed correctly.
    """
    def test_converged(self):
        result = get_result("convergence/converged", "single_configuration_calculation_converged")
        self.assertTrue(result)

    def test_non_converged(self):
        result = get_result("convergence/non_converged", "single_configuration_calculation_converged")
        self.assertFalse(result)


130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
class TestForceFiles(unittest.TestCase):
    """Tests that different force files that can be output, can actually be
    found and parsed.
    """
    def test_single_point(self):

        result = get_result("force_file/single_point", "atom_forces")
        expected_result = convert_unit(
            np.array([
                [0.00000000, 0.00000000, 0.00000000],
                [0.00000000, 0.00000001, 0.00000001],
                [0.00000001, 0.00000001, 0.00000000],
                [0.00000001, 0.00000000, 0.00000001],
                [-0.00000001, -0.00000001, -0.00000001],
                [-0.00000001, -0.00000001, -0.00000001],
                [-0.00000001, -0.00000001, -0.00000001],
                [-0.00000001, -0.00000001, -0.00000001],
            ]),
            "forceAu"
        )
        self.assertTrue(np.array_equal(result, expected_result))


153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
class TestSelfInteractionCorrectionMethod(unittest.TestCase):
    """Tests that the self-interaction correction can be properly parsed.
    """
    def test_no(self):
        sic = get_result("sic/no", "self_interaction_correction_method")
        self.assertEqual(sic, "")

    def test_ad(self):
        sic = get_result("sic/ad", "self_interaction_correction_method")
        self.assertEqual(sic, "SIC_AD")

    def test_explicit_orbitals(self):
        sic = get_result("sic/explicit_orbitals", "self_interaction_correction_method")
        self.assertEqual(sic, "SIC_EXPLICIT_ORBITALS")

    def test_mauri_spz(self):
        sic = get_result("sic/mauri_spz", "self_interaction_correction_method")
        self.assertEqual(sic, "SIC_MAURI_SPZ")

    def test_mauri_us(self):
        sic = get_result("sic/mauri_us", "self_interaction_correction_method")
        self.assertEqual(sic, "SIC_MAURI_US")


177
178
179
180
181
class TestStressTensorMethods(unittest.TestCase):
    """Tests that the stress tensor can be properly parsed for different
    calculation methods.
    """
    def test_none(self):
182
        get_result("stress_tensor/none", "section_stress_tensor")
183
184

    def test_analytical(self):
185
        results = get_result("stress_tensor/analytical")
186
187
188
189
190
        method = results["stress_tensor_method"]
        results["stress_tensor"]
        self.assertEqual(method, "Analytical")

    def test_numerical(self):
191
        results = get_result("stress_tensor/numerical")
192
193
194
195
196
        method = results["stress_tensor_method"]
        results["stress_tensor"]
        self.assertEqual(method, "Numerical")

    def test_diagonal_analytical(self):
197
        results = get_result("stress_tensor/diagonal_analytical")
198
199
200
201
202
        method = results["stress_tensor_method"]
        results["stress_tensor"]
        self.assertEqual(method, "Diagonal analytical")

    def test_diagonal_numerical(self):
203
        results = get_result("stress_tensor/diagonal_numerical")
204
205
206
207
208
        method = results["stress_tensor_method"]
        results["stress_tensor"]
        self.assertEqual(method, "Diagonal numerical")


209
210
211
class TestConfigurationPeriodicDimensions(unittest.TestCase):
    """Tests that the self-interaction correction can be properly parsed.
    """
212
    folder = "configuration_periodic_dimensions/"
213
214

    def test_default(self):
215
        result = get_result(self.folder+"default", "configuration_periodic_dimensions")
216
217
218
        self.assertTrue(np.array_equal(result, np.array((True, True, True))))

    def test_none(self):
219
        result = get_result(self.folder+"none", "configuration_periodic_dimensions")
220
221
222
        self.assertTrue(np.array_equal(result, np.array((False, False, False))))

    def test_x(self):
223
        result = get_result(self.folder+"x", "configuration_periodic_dimensions")
224
225
226
        self.assertTrue(np.array_equal(result, np.array((True, False, False))))

    def test_y(self):
227
        result = get_result(self.folder+"y", "configuration_periodic_dimensions")
228
229
230
        self.assertTrue(np.array_equal(result, np.array((False, True, False))))

    def test_z(self):
231
        result = get_result(self.folder+"z", "configuration_periodic_dimensions")
232
233
234
        self.assertTrue(np.array_equal(result, np.array((False, False, True))))

    def test_xy(self):
235
        result = get_result(self.folder+"xy", "configuration_periodic_dimensions")
236
237
238
        self.assertTrue(np.array_equal(result, np.array((True, True, False))))

    def test_xyz(self):
239
        result = get_result(self.folder+"xyz", "configuration_periodic_dimensions")
240
241
242
        self.assertTrue(np.array_equal(result, np.array((True, True, True))))

    def test_xz(self):
243
        result = get_result(self.folder+"xz", "configuration_periodic_dimensions")
244
245
246
        self.assertTrue(np.array_equal(result, np.array((True, False, True))))

    def test_yz(self):
247
        result = get_result(self.folder+"yz", "configuration_periodic_dimensions")
248
249
250
        self.assertTrue(np.array_equal(result, np.array((False, True, True))))


251
252
253
class TestEnergyForce(unittest.TestCase):
    """Tests for a CP2K calculation with RUN_TYPE ENERGY_FORCE.
    """
254
255
    @classmethod
    def setUpClass(cls):
256
        cls.results = get_result("energy_force")
257
        # cls.results.print_summary()
258

259
    def test_energy_total_scf_iteration(self):
260
        result = self.results["energy_total_scf_iteration"]
261
        expected_result = convert_unit(np.array(-32.2320848878), "hartree")
262
        self.assertTrue(np.array_equal(result[0], expected_result))
263

264
265
266
267
268
269
270
271
    def test_program_name(self):
        result = self.results["program_name"]
        self.assertEqual(result, "CP2K")

    def test_program_compilation_host(self):
        result = self.results["program_compilation_host"]
        self.assertEqual(result, "lenovo700")

272
273
274
275
    def test_scf_max_iteration(self):
        result = self.results["scf_max_iteration"]
        self.assertEqual(result, 300)

276
277
278
279
    def test_basis_set(self):
        section_basis_set = self.results["section_basis_set"][0]

        # Basis name
280
        name = section_basis_set["basis_set_name"]
281
282
283
        self.assertEqual(name, "DZVP-GTH-PADE_PW_150.0")

        # Basis kind
284
        kind = section_basis_set["basis_set_kind"]
285
286
287
        self.assertEqual(kind, "wavefunction")

        # Cell dependent basis mapping
288
        cell_basis_mapping = section_basis_set["mapping_section_basis_set_cell_dependent"]
289
290
291
        self.assertEqual(cell_basis_mapping, 0)

        # # Atom centered basis mapping
292
        atom_basis_mapping = section_basis_set["mapping_section_basis_set_atom_centered"]
293
294
        self.assertTrue(np.array_equal(atom_basis_mapping, np.array(8*[0])))

295
296
    def test_scf_threshold_energy_change(self):
        result = self.results["scf_threshold_energy_change"]
297
        self.assertEqual(result, convert_unit(1.00E-07, "hartree"))
298
299
300
301
302

    def test_number_of_spin_channels(self):
        result = self.results["number_of_spin_channels"]
        self.assertEqual(result, 1)

303
304
305
306
    def test_electronic_structure_method(self):
        result = self.results["electronic_structure_method"]
        self.assertEqual(result, "DFT")

307
308
309
310
311
312
313
314
315
316
    def test_energy_change_scf_iteration(self):
        energy_change = self.results["energy_change_scf_iteration"]
        expected_result = convert_unit(np.array(-3.22E+01), "hartree")
        self.assertTrue(np.array_equal(energy_change[0], expected_result))

    def test_energy_XC_scf_iteration(self):
        result = self.results["energy_XC_scf_iteration"]
        expected_result = convert_unit(np.array(-9.4555961214), "hartree")
        self.assertTrue(np.array_equal(result[0], expected_result))

317
    def test_energy_total(self):
318
        result = self.results["energy_total"]
319
        expected_result = convert_unit(np.array(-31.297885372811074), "hartree")
320
        self.assertTrue(np.array_equal(result, expected_result))
321

322
323
    def test_electronic_kinetic_energy(self):
        result = self.results["electronic_kinetic_energy"]
324
        expected_result = convert_unit(np.array(13.31525592466419), "hartree")
325
326
        self.assertTrue(np.array_equal(result, expected_result))

327
    def test_atom_forces(self):
328
        result = self.results["atom_forces"]
329
330
331
332
333
334
335
336
337
338
339
        expected_result = convert_unit(
            np.array([
                [0.00000000, 0.00000000, 0.00000000],
                [0.00000000, 0.00000001, 0.00000001],
                [0.00000001, 0.00000001, 0.00000000],
                [0.00000001, 0.00000000, 0.00000001],
                [-0.00000001, -0.00000001, -0.00000001],
                [-0.00000001, -0.00000001, -0.00000001],
                [-0.00000001, -0.00000001, -0.00000001],
                [-0.00000001, -0.00000001, -0.00000001],
            ]),
340
            "forceAu"
341
        )
342
        self.assertTrue(np.array_equal(result, expected_result))
343

344
    def test_atom_label(self):
345
        atom_labels = self.results["atom_labels"]
346
        expected_labels = np.array(8*["Si"])
Lauri Himanen's avatar
Lauri Himanen committed
347
348
349
        self.assertTrue(np.array_equal(atom_labels, expected_labels))

    def test_simulation_cell(self):
Lauri Himanen's avatar
Lauri Himanen committed
350
        cell = self.results["simulation_cell"]
Lauri Himanen's avatar
Lauri Himanen committed
351
352
353
354
355
356
357
358
        n_vectors = cell.shape[0]
        n_dim = cell.shape[1]
        self.assertEqual(n_vectors, 3)
        self.assertEqual(n_dim, 3)
        expected_cell = convert_unit(np.array([[5.431, 0, 0], [0, 5.431, 0], [0, 0, 5.431]]), "angstrom")
        self.assertTrue(np.array_equal(cell, expected_cell))

    def test_number_of_atoms(self):
Lauri Himanen's avatar
Lauri Himanen committed
359
        n_atoms = self.results["number_of_atoms"]
Lauri Himanen's avatar
Lauri Himanen committed
360
361
362
        self.assertEqual(n_atoms, 8)

    def test_atom_position(self):
363
        atom_position = self.results["atom_positions"]
Lauri Himanen's avatar
Lauri Himanen committed
364
365
366
        expected_position = convert_unit(np.array([4.073023, 4.073023, 1.357674]), "angstrom")
        self.assertTrue(np.array_equal(atom_position[-1, :], expected_position))

367
    def test_x_cp2k_filenames(self):
368
        input_filename = self.results["x_cp2k_input_filename"]
369
370
371
        expected_input = "si_bulk8.inp"
        self.assertTrue(input_filename, expected_input)

372
        bs_filename = self.results["x_cp2k_basis_set_filename"]
373
374
375
        expected_bs = "../BASIS_SET"
        self.assertEqual(bs_filename, expected_bs)

376
        geminal_filename = self.results["x_cp2k_geminal_filename"]
377
378
379
        expected_geminal = "BASIS_GEMINAL"
        self.assertEqual(geminal_filename, expected_geminal)

380
        potential_filename = self.results["x_cp2k_potential_filename"]
381
382
383
        expected_potential = "../GTH_POTENTIALS"
        self.assertEqual(potential_filename, expected_potential)

384
        mm_potential_filename = self.results["x_cp2k_mm_potential_filename"]
385
386
387
        expected_mm_potential = "MM_POTENTIAL"
        self.assertEqual(mm_potential_filename, expected_mm_potential)

388
        coordinate_filename = self.results["x_cp2k_coordinate_filename"]
389
390
        expected_coordinate = "__STD_INPUT__"
        self.assertEqual(coordinate_filename, expected_coordinate)
391

392
    def test_target_multiplicity(self):
393
        multiplicity = self.results["spin_target_multiplicity"]
394
395
396
397
398
399
        self.assertEqual(multiplicity, 1)

    def test_total_charge(self):
        charge = self.results["total_charge"]
        self.assertEqual(charge, 0)

400
401
    def test_section_basis_set_atom_centered(self):
        basis = self.results["section_basis_set_atom_centered"][0]
402
403
        name = basis["basis_set_atom_centered_short_name"]
        number = basis["basis_set_atom_number"]
404
405
406
        self.assertEquals(name, "DZVP-GTH-PADE")
        self.assertEquals(number, 14)

Lauri Himanen's avatar
Lauri Himanen committed
407
408
    def test_section_basis_set_cell_dependent(self):
        basis = self.results["section_basis_set_cell_dependent"][0]
409
        cutoff = basis["basis_set_planewave_cutoff"]
Lauri Himanen's avatar
Lauri Himanen committed
410
411
        self.assertEquals(cutoff, convert_unit(300.0, "hartree"))

412
413
    def test_section_method_atom_kind(self):
        kind = self.results["section_method_atom_kind"][0]
414
415
        self.assertEqual(kind["method_atom_kind_atom_number"], 14)
        self.assertEqual(kind["method_atom_kind_label"], "Si")
416
417
418

    def test_section_method_basis_set(self):
        kind = self.results["section_method_basis_set"][0]
419
420
        self.assertEqual(kind["method_basis_set_kind"], "wavefunction")
        self.assertEqual(kind["number_of_basis_sets_atom_centered"], 1)
421
        self.assertTrue(np.array_equal(kind["mapping_section_method_basis_set_atom_centered"], np.array([[0, 0]])))
422

423
424
425
426
427
    def test_single_configuration_calculation_converged(self):
        result = self.results["single_configuration_calculation_converged"]
        self.assertTrue(result)

    def test_scf_dft_number_of_iterations(self):
428
        result = self.results["number_of_scf_iterations"]
429
430
431
432
433
434
435
        self.assertEqual(result, 10)

    def test_single_configuration_to_calculation_method_ref(self):
        result = self.results["single_configuration_to_calculation_method_ref"]
        self.assertEqual(result, 0)

    def test_single_configuration_calculation_to_system_description_ref(self):
436
        result = self.results["single_configuration_calculation_to_system_ref"]
437
438
        self.assertEqual(result, 0)

439
    def test_stress_tensor(self):
440
        result = self.results["stress_tensor"]
441
442
        expected_result = convert_unit(
            np.array([
443
444
445
                [7.77640934, -0.00000098, -0.00000099],
                [-0.00000098, 7.77640935, -0.00000101],
                [-0.00000099, -0.00000101, 7.77640935],
446
447
448
449
450
            ]),
            "GPa"
        )
        self.assertTrue(np.array_equal(result, expected_result))

451
452
    def test_stress_tensor_eigenvalues(self):
        result = self.results["x_cp2k_stress_tensor_eigenvalues"]
453
        expected_result = convert_unit(np.array([7.77640735, 7.77641033, 7.77641036]), "GPa")
454
455
456
457
458
        self.assertTrue(np.array_equal(result, expected_result))

    def test_stress_tensor_eigenvectors(self):
        result = self.results["x_cp2k_stress_tensor_eigenvectors"]
        expected_result = np.array([
459
460
461
            [0.57490332, -0.79965737, -0.17330395],
            [0.57753686, 0.54662171, -0.60634634],
            [0.57960102, 0.24850110, 0.77608624],
462
463
464
465
466
        ])
        self.assertTrue(np.array_equal(result, expected_result))

    def test_stress_tensor_determinant(self):
        result = self.results["x_cp2k_stress_tensor_determinant"]
467
        expected_result = convert_unit(4.70259243E+02, "GPa^3")
468
469
470
471
        self.assertTrue(np.array_equal(result, expected_result))

    def test_stress_tensor_one_third_of_trace(self):
        result = self.results["x_cp2k_stress_tensor_one_third_of_trace"]
472
        expected_result = convert_unit(7.77640934E+00, "GPa")
473
474
        self.assertTrue(np.array_equal(result, expected_result))

475
476
    def test_program_basis_set_type(self):
        result = self.results["program_basis_set_type"]
477
        self.assertEqual(result, "gaussians")
478
479
480


class TestPreprocessor(unittest.TestCase):
481
482
483
    """Tests that the parser can read input files with preprocessor
    declarations, such as variables.
    """
484
    def test_include(self):
485
        result = get_result("input_preprocessing/include", "x_cp2k_input_GLOBAL.PRINT_LEVEL")
486
487
488
        self.assertEqual(result, "LOW")

    def test_variable(self):
489
        result = get_result("input_preprocessing/variable", "x_cp2k_input_GLOBAL.PROJECT_NAME")
490
491
        self.assertEqual(result, "variable_test")

492
    def test_variable_multiple(self):
493
        result = get_result("input_preprocessing/variable_multiple", "x_cp2k_input_FORCE_EVAL.DFT.MGRID.CUTOFF")
494
        self.assertEqual(result, "50")
495

Lauri Himanen's avatar
Lauri Himanen committed
496
    def test_comments(self):
497
        result = get_result("input_preprocessing/comments", "x_cp2k_input_FORCE_EVAL.DFT.MGRID.CUTOFF")
498
        self.assertEqual(result, "120")
Lauri Himanen's avatar
Lauri Himanen committed
499

500
    def test_tabseparator(self):
501
        result = get_result("input_preprocessing/tabseparator", "x_cp2k_input_FORCE_EVAL.DFT.MGRID.CUTOFF")
502
        self.assertEqual(result, "120")
503

504

505
class TestGeoOpt(unittest.TestCase):
506
507
    """Tests that geometry optimizations are correctly parsed.
    """
508
509
    @classmethod
    def setUpClass(cls):
510
        cls.results = get_result("geo_opt/cg")
511
512
513
514
515

    def test_geometry_optimization_converged(self):
        result = self.results["geometry_optimization_converged"]
        self.assertTrue(result)

516
    def test_number_of_frames_in_sequence(self):
517
        result = self.results["number_of_frames_in_sequence"]
518
519
        self.assertEqual(result, 7)

520
521
522
523
524
525
    def test_frame_sequence_to_sampling_ref(self):
        result = self.results["frame_sequence_to_sampling_ref"]
        self.assertEqual(result, 0)

    def test_frame_sequence_local_frames_ref(self):
        result = self.results["frame_sequence_local_frames_ref"]
526
        expected_result = np.array([0, 1, 2, 3, 4, 5, 6])
527
528
        self.assertTrue(np.array_equal(result, expected_result))

529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
    def test_sampling_method(self):
        result = self.results["sampling_method"]
        self.assertEqual(result, "geometry_optimization")

    def test_geometry_optimization_method(self):
        result = self.results["geometry_optimization_method"]
        self.assertEqual(result, "conjugate_gradient")

    def test_geometry_optimization_geometry_change(self):
        result = self.results["geometry_optimization_geometry_change"]
        expected_result = convert_unit(
            0.0010000000,
            "bohr"
        )
        self.assertEqual(result, expected_result)

    def test_geometry_optimization_threshold_force(self):
        result = self.results["geometry_optimization_threshold_force"]
        expected_result = convert_unit(
            0.0010000000,
            "bohr^-1*hartree"
        )
        self.assertEqual(result, expected_result)

    def test_frame_sequence_potential_energy(self):
        result = self.results["frame_sequence_potential_energy"]
        expected_result = convert_unit(
            np.array([
                -17.1534159246,
                -17.1941015290,
                -17.2092321965,
                -17.2097667733,
                -17.2097743028,
                -17.2097743229,
563
                -17.20977820662248,
564
565
566
567
568
            ]),
            "hartree"
        )
        self.assertTrue(np.array_equal(result, expected_result))

569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
    def test_atom_positions(self):
        result = self.results["atom_positions"]
        expected_start = convert_unit(
            np.array([
                [12.2353220000, 1.3766420000, 10.8698800000],
                [12.4175775999, 2.2362362573, 11.2616216864],
                [11.9271436933, 1.5723516602, 10.0115134757],
            ]),
            "angstrom"
        )

        expected_end = convert_unit(
            np.array([
                [12.2353220000, 1.3766420000, 10.8698800000],
                [12.4958164689, 2.2307248873, 11.3354322515],
                [11.9975558616, 1.5748085240, 10.0062792262],
            ]),
            "angstrom"
        )
588
589
        result_start = result[0, :, :]
        result_end = result[-1, :, :]
590
591
        self.assertTrue(np.array_equal(result_start, expected_start))
        self.assertTrue(np.array_equal(result_end, expected_end))
592
593


594
class TestGeoOptTrajFormats(unittest.TestCase):
595
596
    """Different trajectory formats in geometry optimization.
    """
597
598
    def test_default(self):
        result = get_result("geo_opt/geometry_formats/default", "atom_positions")
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
        expected_start = convert_unit(
            np.array([
                [12.2353220000, 1.3766420000, 10.8698800000],
                [12.4175624065, 2.2362390825, 11.2616392180],
                [11.9271777126, 1.5723402996, 10.0115089094],
            ]),
            "angstrom"
        )
        expected_end = convert_unit(
            np.array([
                [12.2353220000, 1.3766420000, 10.8698800000],
                [12.4957995882, 2.2307218433, 11.3354453867],
                [11.9975764125, 1.5747996320, 10.0062529540],
            ]),
            "angstrom"
        )
615
616
        result_start = result[0, :, :]
        result_end = result[-1, :, :]
617
618
619
        self.assertTrue(np.array_equal(result_start, expected_start))
        self.assertTrue(np.array_equal(result_end, expected_end))

620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
    # def test_xyz(self):

        # result = get_result("geo_opt/geometry_formats/xyz", "atom_positions")
        # expected_start = convert_unit(
            # np.array([
                # [12.2353220000, 1.3766420000, 10.8698800000],
                # [12.4175624065, 2.2362390825, 11.2616392180],
                # [11.9271777126, 1.5723402996, 10.0115089094],
            # ]),
            # "angstrom"
        # )
        # expected_end = convert_unit(
            # np.array([
                # [12.2353220000, 1.3766420000, 10.8698800000],
                # [12.4957995882, 2.2307218433, 11.3354453867],
                # [11.9975764125, 1.5747996320, 10.0062529540],
            # ]),
            # "angstrom"
        # )
        # result_start = result[0, :, :]
        # result_end = result[-1, :, :]
        # self.assertTrue(np.array_equal(result_start, expected_start))
        # self.assertTrue(np.array_equal(result_end, expected_end))

    # def test_pdb(self):
        # result = get_result("geo_opt/geometry_formats/pdb", "atom_positions")
        # expected_start = convert_unit(
            # np.array([
                # [12.235, 1.377, 10.870],
                # [12.418, 2.236, 11.262],
                # [11.927, 1.572, 10.012],
            # ]),
            # "angstrom"
        # )
        # expected_end = convert_unit(
            # np.array([
                # [12.235, 1.377, 10.870],
                # [12.496, 2.231, 11.335],
                # [11.998, 1.575, 10.006],
            # ]),
            # "angstrom"
        # )
        # result_start = result[0, :, :]
        # result_end = result[-1, :, :]
        # self.assertTrue(np.array_equal(result_start, expected_start))
        # self.assertTrue(np.array_equal(result_end, expected_end))

    # def test_dcd(self):
        # result = get_result("geo_opt/geometry_formats/dcd", "atom_positions")
        # frames = result.shape[0]
        # self.assertEqual(frames, 7)
671
672


673
class TestGeoOptOptimizers(unittest.TestCase):
674
675
    """Different optimization methods in gemoetry optimization.
    """
676
677
678
679
680
681
682
683
    def test_bfgs(self):
        result = get_result("geo_opt/bfgs", "geometry_optimization_method")
        self.assertEqual(result, "bfgs")

    def test_lbfgs(self):
        result = get_result("geo_opt/lbfgs", "geometry_optimization_method")
        self.assertEqual(result, "bfgs")

684
685

class TestGeoOptTrajectory(unittest.TestCase):
686
687
688
    """Tests that the print settings for geometry optimization are handled
    correctly.
    """
689
690
691
692
    def test_each_and_add_last(self):
        """Test that the EACH and ADD_LAST settings affect the parsing
        correctly.
        """
693
        results = get_result("geo_opt/each")
694
695
696
697
698

        single_conf = results["section_single_configuration_calculation"]
        systems = results["section_system"]

        i_conf = 0
699
        for calc in single_conf:
700
            system_index = calc["single_configuration_calculation_to_system_ref"]
701
702
703
            system = systems[system_index]

            if i_conf == 0 or i_conf == 2 or i_conf == 4:
704
                with self.assertRaises(KeyError):
705
                    pos = system["atom_positions"]
706
            else:
707
                pos = system["atom_positions"]
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
                if i_conf == 1:
                    expected_pos = convert_unit(
                        np.array([
                            [12.2353220000, 1.3766420000, 10.8698800000],
                            [12.4618486015, 2.2314871691, 11.3335607388],
                            [11.9990227122, 1.5776813026, 10.0384213366],
                        ]),
                        "angstrom"
                    )
                    self.assertTrue(np.array_equal(pos, expected_pos))
                if i_conf == 3:
                    expected_pos = convert_unit(
                        np.array([
                            [12.2353220000, 1.3766420000, 10.8698800000],
                            [12.4962705528, 2.2308411983, 11.3355758433],
                            [11.9975151486, 1.5746309898, 10.0054430868],
                        ]),
                        "angstrom"
                    )
                    self.assertTrue(np.array_equal(pos, expected_pos))
                if i_conf == 5:
                    expected_pos = convert_unit(
                        np.array([
                            [12.2353220000, 1.3766420000, 10.8698800000],
                            [12.4958168364, 2.2307249171, 11.3354322532],
                            [11.9975556812, 1.5748088251, 10.0062793864],
                        ]),
                        "angstrom"
                    )
                    self.assertTrue(np.array_equal(pos, expected_pos))

                if i_conf == 6:
                    expected_pos = convert_unit(
                        np.array([
                            [12.2353220000, 1.3766420000, 10.8698800000],
                            [12.4958164689, 2.2307248873, 11.3354322515],
                            [11.9975558616, 1.5748085240, 10.0062792262],
                        ]),
                        "angstrom"
                    )
                    self.assertTrue(np.array_equal(pos, expected_pos))

            i_conf += 1

Lauri Himanen's avatar
Lauri Himanen committed
752

753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
class TestRestart(unittest.TestCase):
    """Used to test that restarts from a previous run are parsed correctly.
    """
    @classmethod
    def setUpClass(cls):
        cls.results = get_result("restart")

    def test_program_name(self):
        result = self.results["program_name"]
        self.assertEqual(result, "CP2K")

    def test_frames(self):
        """Test that the number of frames equals the actual number of steps
        that are present in the output file.
        """
        scc = self.results["section_single_configuration_calculation"]
        n_scc = len(scc)
        n_frames = self.results["number_of_frames_in_sequence"]
        self.assertEqual(n_scc, n_frames)


class TestMDPrintSettings(unittest.TestCase):
    """Tests that MD print settings are respected.
    """
    @classmethod
    def setUpClass(cls):
        cls.results = get_result("md/print_settings")

    def test_traj_print(self):
        result = self.results["number_of_frames_in_sequence"]
        self.assertEqual(result, 11)

        sccs = self.results["section_single_configuration_calculation"]
        systems = self.results["section_system"]
        for i_scc, scc in enumerate(sccs):

            # Check that the configuration was read correctly from the XYZ file
            # by taking into account the print settings
            if i_scc % 3 == 0:

                # Check that positions are found
                ref_system = scc["single_configuration_calculation_to_system_ref"]
                system = systems[ref_system]
                system["atom_positions"]

        # Check that frames from ener file are read correctly with interval of 3
        frames = [0, 3, 6, 9]
        self.results["frame_sequence_temperature"]
        temp_frames = self.results["frame_sequence_temperature_frames"]
        pot_frames = self.results["frame_sequence_potential_energy_frames"]
        kin_frames = self.results["frame_sequence_kinetic_energy_frames"]
        con_frames = self.results["frame_sequence_conserved_quantity_frames"]
        self.assertTrue(np.array_equal(temp_frames, frames))
        self.assertTrue(np.array_equal(pot_frames, frames))
        self.assertTrue(np.array_equal(kin_frames, frames))
        self.assertTrue(np.array_equal(con_frames, frames))


Lauri Himanen's avatar
Lauri Himanen committed
811
class TestMD(unittest.TestCase):
812
813
    """Molecular dynamics tests.
    """
Lauri Himanen's avatar
Lauri Himanen committed
814
815
    @classmethod
    def setUpClass(cls):
816
        cls.results = get_result("md/nve")
Lauri Himanen's avatar
Lauri Himanen committed
817
818
819
820
821
822
823
824
825
        cls.temp = convert_unit(
            np.array([
                300.000000000,
                275.075405378,
                235.091633019,
                202.752506973,
                192.266488819,
                201.629598676,
                218.299664775,
826
                230.324748557,
Lauri Himanen's avatar
Lauri Himanen committed
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
                232.691881533,
                226.146979313,
                213.165337396,
            ]),
            "K"
        )
        cls.cons = convert_unit(
            np.array([
                -34.323271136,
                -34.323245645,
                -34.323206964,
                -34.323183380,
                -34.323187747,
                -34.323208962,
                -34.323227533,
                -34.323233583,
                -34.323230715,
                -34.323227013,
                -34.323224123,
            ]),
            "hartree"
        )
        cls.pot = convert_unit(
            np.array([
                -34.330396471,
                -34.329778993,
                -34.328790653,
                -34.327998978,
                -34.327754290,
                -34.327997890,
                -34.328412394,
                -34.328704052,
                -34.328757407,
                -34.328598255,
                -34.328287038,
            ]),
            "hartree"
        )
        cls.kin = convert_unit(
            np.array([
                0.007125335,
                0.006533348,
                0.005583688,
                0.004815598,
                0.004566544,
                0.004788928,
                0.005184860,
                0.005470470,
                0.005526692,
                0.005371243,
                0.005062914,
            ]),
            "hartree"
        )
Lauri Himanen's avatar
Lauri Himanen committed
881

882
883
884
885
886
    def test_number_of_atoms(self):
        result = self.results["number_of_atoms"]
        expected_result = np.array(11*[6])
        self.assertTrue(np.array_equal(result, expected_result))

887
888
889
890
891
892
893
894
895
896
897
898
899
    def test_simulation_cell(self):
        result = self.results["simulation_cell"]
        self.assertEqual(len(result), 11)
        expected_start = convert_unit(
            np.array([
                [9.853, 0, 0],
                [0, 9.853, 0],
                [0, 0, 9.853],
            ]),
            "angstrom"
        )
        self.assertTrue(np.array_equal(result[0], expected_start))

Lauri Himanen's avatar
Lauri Himanen committed
900
    def test_ensemble_type(self):
901
        result = self.results["ensemble_type"]
Lauri Himanen's avatar
Lauri Himanen committed
902
903
        self.assertEqual(result, "NVE")

904
905
906
907
    def test_sampling_method(self):
        result = self.results["sampling_method"]
        self.assertEqual(result, "molecular_dynamics")

Lauri Himanen's avatar
Lauri Himanen committed
908
909
910
911
    def test_number_of_frames_in_sequence(self):
        result = self.results["number_of_frames_in_sequence"]
        self.assertEqual(result, 11)

912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
    def test_atom_positions(self):
        result = self.results["atom_positions"]
        expected_start = convert_unit(
            np.array([
                [2.2803980000, 9.1465390000, 5.0886960000],
                [1.2517030000, 2.4062610000, 7.7699080000],
                [1.7620190000, 9.8204290000, 5.5284540000],
                [3.0959870000, 9.1070880000, 5.5881860000],
                [0.5541290000, 2.9826340000, 8.0820240000],
                [1.7712570000, 2.9547790000, 7.1821810000],
            ]),
            "angstrom"
        )
        expected_end = convert_unit(
            np.array([
                [2.2916014875, 9.1431763260, 5.0868100688],
                [1.2366834078, 2.4077552776, 7.7630044423],
                [1.6909790671, 9.8235337924, 5.5042564094],
                [3.1130341664, 9.0372111810, 5.6100739746],
                [0.5652070478, 3.0441761067, 8.1734257299],
                [1.8669280879, 2.9877213524, 7.2364955946],
            ]),
            "angstrom"
        )
936
937
        self.assertTrue(np.array_equal(result[0, :], expected_start))
        self.assertTrue(np.array_equal(result[-1, :], expected_end))
938

939
940
941
942
943
944
945
946
947
948
949
    def test_atom_velocities(self):
        result = self.results["atom_velocities"]
        expected_start = convert_unit(
            np.array([
                [0.0000299284, 0.0000082360, -0.0000216368],
                [-0.0001665963, 0.0001143863, -0.0000622640],
                [-0.0005732926, -0.0003112611, -0.0007149779],
                [0.0013083605, -0.0009262219, 0.0006258560],
                [0.0012002313, -0.0003701042, 0.0002810523],
                [0.0002340810, -0.0003388418, 0.0011398583],
            ]),
950
            "bohr*(hbar/hartree)^-1"
951
952
953
954
955
956
957
958
959
960
        )
        expected_end = convert_unit(
            np.array([
                [0.0001600263, -0.0000383308, 0.0000153662],
                [-0.0001269381, -0.0000005151, -0.0000726214],
                [0.0000177093, -0.0003257814, -0.0000257852],
                [-0.0015067045, -0.0001700489, -0.0003651605],
                [0.0000307926, 0.0006886719, 0.0008431321],
                [0.0007424681, 0.0003614127, 0.0005749089],
            ]),
961
            "bohr*(hbar/hartree)^-1"
962
963
        )

964
965
        self.assertTrue(np.array_equal(result[0, :], expected_start))
        self.assertTrue(np.array_equal(result[-1, :], expected_end))
966

967
968
    def test_frame_sequence_potential_energy(self):
        result = self.results["frame_sequence_potential_energy"]
Lauri Himanen's avatar
Lauri Himanen committed
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
        self.assertTrue(np.array_equal(result, self.pot))

    def test_frame_sequence_kinetic_energy(self):
        result = self.results["frame_sequence_kinetic_energy"]
        self.assertTrue(np.array_equal(result, self.kin))

    def test_frame_sequence_conserved_quantity(self):
        result = self.results["frame_sequence_conserved_quantity"]
        self.assertTrue(np.array_equal(result, self.cons))

    def test_frame_sequence_temperature(self):
        result = self.results["frame_sequence_temperature"]
        self.assertTrue(np.array_equal(result, self.temp))

    def test_frame_sequence_time(self):
        result = self.results["frame_sequence_time"]
985
986
        expected_result = convert_unit(
            np.array([
Lauri Himanen's avatar
Lauri Himanen committed
987
988
989
990
991
992
993
994
995
996
997
                0.000000,
                0.500000,
                1.000000,
                1.500000,
                2.000000,
                2.500000,
                3.000000,
                3.500000,
                4.000000,
                4.500000,
                5.000000,
998
            ]),
Lauri Himanen's avatar
Lauri Himanen committed
999
            "fs"
1000
        )