regtests.py 39.7 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
class TestUnknownInput(unittest.TestCase):
    """Tests for cases where unknown information is encountered in the parsing.
    """
54
    def test_unknown_version(self):
55
56
        """Test how a new version is handled.
        """
57
58
        get_result("errors/unknown_version", "XC_functional")

59
    def test_unknown_input_keyword(self):
60
61
        """Test how an unknown input keyword is handled.
        """
62
63
64
        get_result("errors/unknown_input_keyword", "XC_functional")

    def test_unknown_input_section(self):
65
66
        """Test unknown input file section.
        """
67
68
69
        get_result("errors/unknown_input_section", "XC_functional")

    def test_unknown_input_section_parameter(self):
70
71
        """
        """
72
73
        get_result("errors/unknown_input_section_parameter", "XC_functional")

74

75
class TestXCFunctional(unittest.TestCase):
76
77
    """Tests that the XC functionals can be properly parsed.
    """
78
79
    def test_pade(self):
        xc = get_result("XC_functional/pade", "XC_functional")
80
        self.assertEqual(xc, "1*LDA_XC_TETER93")
81
82
83

    def test_lda(self):
        xc = get_result("XC_functional/lda", "XC_functional")
84
        self.assertEqual(xc, "1*LDA_XC_TETER93")
85
86
87

    def test_blyp(self):
        xc = get_result("XC_functional/blyp", "XC_functional")
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
        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")
109
110


111
112
113
114
115
116
117
118
119
120
121
122
123
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)


124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
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))


147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
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")


171
172
173
174
175
class TestStressTensorMethods(unittest.TestCase):
    """Tests that the stress tensor can be properly parsed for different
    calculation methods.
    """
    def test_none(self):
176
        get_result("stress_tensor/none", "section_stress_tensor")
177
178

    def test_analytical(self):
179
        results = get_result("stress_tensor/analytical")
180
181
182
183
184
        method = results["stress_tensor_method"]
        results["stress_tensor"]
        self.assertEqual(method, "Analytical")

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

    def test_diagonal_analytical(self):
191
        results = get_result("stress_tensor/diagonal_analytical")
192
193
194
195
196
        method = results["stress_tensor_method"]
        results["stress_tensor"]
        self.assertEqual(method, "Diagonal analytical")

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


203
204
205
class TestConfigurationPeriodicDimensions(unittest.TestCase):
    """Tests that the self-interaction correction can be properly parsed.
    """
206
    folder = "configuration_periodic_dimensions/"
207
208

    def test_default(self):
209
        result = get_result(self.folder+"default", "configuration_periodic_dimensions")
210
211
212
        self.assertTrue(np.array_equal(result, np.array((True, True, True))))

    def test_none(self):
213
        result = get_result(self.folder+"none", "configuration_periodic_dimensions")
214
215
216
        self.assertTrue(np.array_equal(result, np.array((False, False, False))))

    def test_x(self):
217
        result = get_result(self.folder+"x", "configuration_periodic_dimensions")
218
219
220
        self.assertTrue(np.array_equal(result, np.array((True, False, False))))

    def test_y(self):
221
        result = get_result(self.folder+"y", "configuration_periodic_dimensions")
222
223
224
        self.assertTrue(np.array_equal(result, np.array((False, True, False))))

    def test_z(self):
225
        result = get_result(self.folder+"z", "configuration_periodic_dimensions")
226
227
228
        self.assertTrue(np.array_equal(result, np.array((False, False, True))))

    def test_xy(self):
229
        result = get_result(self.folder+"xy", "configuration_periodic_dimensions")
230
231
232
        self.assertTrue(np.array_equal(result, np.array((True, True, False))))

    def test_xyz(self):
233
        result = get_result(self.folder+"xyz", "configuration_periodic_dimensions")
234
235
236
        self.assertTrue(np.array_equal(result, np.array((True, True, True))))

    def test_xz(self):
237
        result = get_result(self.folder+"xz", "configuration_periodic_dimensions")
238
239
240
        self.assertTrue(np.array_equal(result, np.array((True, False, True))))

    def test_yz(self):
241
        result = get_result(self.folder+"yz", "configuration_periodic_dimensions")
242
243
244
        self.assertTrue(np.array_equal(result, np.array((False, True, True))))


245
246
247
class TestEnergyForce(unittest.TestCase):
    """Tests for a CP2K calculation with RUN_TYPE ENERGY_FORCE.
    """
248
249
    @classmethod
    def setUpClass(cls):
250
        cls.results = get_result("energy_force")
251
        # cls.results.print_summary()
252

253
    def test_energy_total_scf_iteration(self):
254
        result = self.results["energy_total_scf_iteration"]
255
        expected_result = convert_unit(np.array(-32.2320848878), "hartree")
256
        self.assertTrue(np.array_equal(result[0], expected_result))
257

258
259
260
261
262
263
264
265
    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")

266
267
268
269
    def test_scf_max_iteration(self):
        result = self.results["scf_max_iteration"]
        self.assertEqual(result, 300)

270
271
272
273
    def test_basis_set(self):
        section_basis_set = self.results["section_basis_set"][0]

        # Basis name
274
        name = section_basis_set["basis_set_name"]
275
276
277
        self.assertEqual(name, "DZVP-GTH-PADE_PW_150.0")

        # Basis kind
278
        kind = section_basis_set["basis_set_kind"]
279
280
281
        self.assertEqual(kind, "wavefunction")

        # Cell dependent basis mapping
282
        cell_basis_mapping = section_basis_set["mapping_section_basis_set_cell_dependent"]
283
284
285
        self.assertEqual(cell_basis_mapping, 0)

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

289
290
    def test_scf_threshold_energy_change(self):
        result = self.results["scf_threshold_energy_change"]
291
        self.assertEqual(result, convert_unit(1.00E-07, "hartree"))
292
293
294
295
296

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

297
298
299
300
    def test_electronic_structure_method(self):
        result = self.results["electronic_structure_method"]
        self.assertEqual(result, "DFT")

301
302
303
304
305
306
307
308
309
310
    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))

311
    def test_energy_total(self):
312
        result = self.results["energy_total"]
313
        expected_result = convert_unit(np.array(-31.297885372811074), "hartree")
314
        self.assertTrue(np.array_equal(result, expected_result))
315

316
317
    def test_electronic_kinetic_energy(self):
        result = self.results["electronic_kinetic_energy"]
318
        expected_result = convert_unit(np.array(13.31525592466419), "hartree")
319
320
        self.assertTrue(np.array_equal(result, expected_result))

321
    def test_atom_forces(self):
322
        result = self.results["atom_forces"]
323
324
325
326
327
328
329
330
331
332
333
        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],
            ]),
334
            "forceAu"
335
        )
336
        self.assertTrue(np.array_equal(result, expected_result))
337

338
    def test_atom_label(self):
339
        atom_labels = self.results["atom_labels"]
340
        expected_labels = np.array(8*["Si"])
Lauri Himanen's avatar
Lauri Himanen committed
341
342
343
        self.assertTrue(np.array_equal(atom_labels, expected_labels))

    def test_simulation_cell(self):
Lauri Himanen's avatar
Lauri Himanen committed
344
        cell = self.results["simulation_cell"]
Lauri Himanen's avatar
Lauri Himanen committed
345
346
347
348
349
350
351
352
        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
353
        n_atoms = self.results["number_of_atoms"]
Lauri Himanen's avatar
Lauri Himanen committed
354
355
356
        self.assertEqual(n_atoms, 8)

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

361
    def test_x_cp2k_filenames(self):
362
        input_filename = self.results["x_cp2k_input_filename"]
363
364
365
        expected_input = "si_bulk8.inp"
        self.assertTrue(input_filename, expected_input)

366
        bs_filename = self.results["x_cp2k_basis_set_filename"]
367
368
369
        expected_bs = "../BASIS_SET"
        self.assertEqual(bs_filename, expected_bs)

370
        geminal_filename = self.results["x_cp2k_geminal_filename"]
371
372
373
        expected_geminal = "BASIS_GEMINAL"
        self.assertEqual(geminal_filename, expected_geminal)

374
        potential_filename = self.results["x_cp2k_potential_filename"]
375
376
377
        expected_potential = "../GTH_POTENTIALS"
        self.assertEqual(potential_filename, expected_potential)

378
        mm_potential_filename = self.results["x_cp2k_mm_potential_filename"]
379
380
381
        expected_mm_potential = "MM_POTENTIAL"
        self.assertEqual(mm_potential_filename, expected_mm_potential)

382
        coordinate_filename = self.results["x_cp2k_coordinate_filename"]
383
384
        expected_coordinate = "__STD_INPUT__"
        self.assertEqual(coordinate_filename, expected_coordinate)
385

386
    def test_target_multiplicity(self):
387
        multiplicity = self.results["spin_target_multiplicity"]
388
389
390
391
392
393
        self.assertEqual(multiplicity, 1)

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

394
395
    def test_section_basis_set_atom_centered(self):
        basis = self.results["section_basis_set_atom_centered"][0]
396
397
        name = basis["basis_set_atom_centered_short_name"]
        number = basis["basis_set_atom_number"]
398
399
400
        self.assertEquals(name, "DZVP-GTH-PADE")
        self.assertEquals(number, 14)

Lauri Himanen's avatar
Lauri Himanen committed
401
402
    def test_section_basis_set_cell_dependent(self):
        basis = self.results["section_basis_set_cell_dependent"][0]
403
        cutoff = basis["basis_set_planewave_cutoff"]
Lauri Himanen's avatar
Lauri Himanen committed
404
405
        self.assertEquals(cutoff, convert_unit(300.0, "hartree"))

406
407
    def test_section_method_atom_kind(self):
        kind = self.results["section_method_atom_kind"][0]
408
409
        self.assertEqual(kind["method_atom_kind_atom_number"], 14)
        self.assertEqual(kind["method_atom_kind_label"], "Si")
410
411
412

    def test_section_method_basis_set(self):
        kind = self.results["section_method_basis_set"][0]
413
414
        self.assertEqual(kind["method_basis_set_kind"], "wavefunction")
        self.assertEqual(kind["number_of_basis_sets_atom_centered"], 1)
415
        self.assertTrue(np.array_equal(kind["mapping_section_method_basis_set_atom_centered"], np.array([[0, 0]])))
416

417
418
419
420
421
    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):
422
        result = self.results["number_of_scf_iterations"]
423
424
425
426
427
428
429
        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):
430
        result = self.results["single_configuration_calculation_to_system_ref"]
431
432
        self.assertEqual(result, 0)

433
    def test_stress_tensor(self):
434
        result = self.results["stress_tensor"]
435
436
        expected_result = convert_unit(
            np.array([
437
438
439
                [7.77640934, -0.00000098, -0.00000099],
                [-0.00000098, 7.77640935, -0.00000101],
                [-0.00000099, -0.00000101, 7.77640935],
440
441
442
443
444
            ]),
            "GPa"
        )
        self.assertTrue(np.array_equal(result, expected_result))

445
446
    def test_stress_tensor_eigenvalues(self):
        result = self.results["x_cp2k_stress_tensor_eigenvalues"]
447
        expected_result = convert_unit(np.array([7.77640735, 7.77641033, 7.77641036]), "GPa")
448
449
450
451
452
        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([
453
454
455
            [0.57490332, -0.79965737, -0.17330395],
            [0.57753686, 0.54662171, -0.60634634],
            [0.57960102, 0.24850110, 0.77608624],
456
457
458
459
460
        ])
        self.assertTrue(np.array_equal(result, expected_result))

    def test_stress_tensor_determinant(self):
        result = self.results["x_cp2k_stress_tensor_determinant"]
461
        expected_result = convert_unit(4.70259243E+02, "GPa^3")
462
463
464
465
        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"]
466
        expected_result = convert_unit(7.77640934E+00, "GPa")
467
468
        self.assertTrue(np.array_equal(result, expected_result))

469
470
    def test_program_basis_set_type(self):
        result = self.results["program_basis_set_type"]
471
        self.assertEqual(result, "gaussians")
472
473
474


class TestPreprocessor(unittest.TestCase):
475
476
477
    """Tests that the parser can read input files with preprocessor
    declarations, such as variables.
    """
478
    def test_include(self):
479
        result = get_result("input_preprocessing/include", "x_cp2k_input_GLOBAL.PRINT_LEVEL")
480
481
482
        self.assertEqual(result, "LOW")

    def test_variable(self):
483
        result = get_result("input_preprocessing/variable", "x_cp2k_input_GLOBAL.PROJECT_NAME")
484
485
        self.assertEqual(result, "variable_test")

486
    def test_variable_multiple(self):
487
        result = get_result("input_preprocessing/variable_multiple", "x_cp2k_input_FORCE_EVAL.DFT.MGRID.CUTOFF")
488
        self.assertEqual(result, "50")
489

Lauri Himanen's avatar
Lauri Himanen committed
490
    def test_comments(self):
491
        result = get_result("input_preprocessing/comments", "x_cp2k_input_FORCE_EVAL.DFT.MGRID.CUTOFF")
492
        self.assertEqual(result, "120")
Lauri Himanen's avatar
Lauri Himanen committed
493

494
    def test_tabseparator(self):
495
        result = get_result("input_preprocessing/tabseparator", "x_cp2k_input_FORCE_EVAL.DFT.MGRID.CUTOFF")
496
        self.assertEqual(result, "120")
497

498

499
class TestGeoOpt(unittest.TestCase):
500
501
    """Tests that geometry optimizations are correctly parsed.
    """
502
503
    @classmethod
    def setUpClass(cls):
504
        cls.results = get_result("geo_opt/cg")
505
506
507
508
509

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

510
    def test_number_of_frames_in_sequence(self):
511
        result = self.results["number_of_frames_in_sequence"]
512
513
        self.assertEqual(result, 7)

514
515
516
517
518
519
    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"]
520
        expected_result = np.array([0, 1, 2, 3, 4, 5, 6])
521
522
        self.assertTrue(np.array_equal(result, expected_result))

523
524
525
526
527
528
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
    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,
557
                -17.20977820662248,
558
559
560
561
562
            ]),
            "hartree"
        )
        self.assertTrue(np.array_equal(result, expected_result))

563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
    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"
        )
582
583
        result_start = result[0, :, :]
        result_end = result[-1, :, :]
584
585
        self.assertTrue(np.array_equal(result_start, expected_start))
        self.assertTrue(np.array_equal(result_end, expected_end))
586
587


588
class TestGeoOptTrajFormats(unittest.TestCase):
589
590
    """Different trajectory formats in geometry optimization.
    """
591
    def test_xyz(self):
592

593
        result = get_result("geo_opt/geometry_formats/xyz", "atom_positions")
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
        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"
        )
610
611
        result_start = result[0, :, :]
        result_end = result[-1, :, :]
612
613
614
615
        self.assertTrue(np.array_equal(result_start, expected_start))
        self.assertTrue(np.array_equal(result_end, expected_end))

    def test_pdb(self):
616
        result = get_result("geo_opt/geometry_formats/pdb", "atom_positions")
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
        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"
        )
633
634
        result_start = result[0, :, :]
        result_end = result[-1, :, :]
635
636
637
638
        self.assertTrue(np.array_equal(result_start, expected_start))
        self.assertTrue(np.array_equal(result_end, expected_end))

    def test_dcd(self):
639
        result = get_result("geo_opt/geometry_formats/dcd", "atom_positions")
640
641
642
643
        frames = result.shape[0]
        self.assertEqual(frames, 7)


644
class TestGeoOptOptimizers(unittest.TestCase):
645
646
    """Different optimization methods in gemoetry optimization.
    """
647
648
649
650
651
652
653
654
    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")

655
656

class TestGeoOptTrajectory(unittest.TestCase):
657
658
659
    """Tests that the print settings for geometry optimization are handled
    correctly.
    """
660
661
662
663
    def test_each_and_add_last(self):
        """Test that the EACH and ADD_LAST settings affect the parsing
        correctly.
        """
664
        results = get_result("geo_opt/each")
665
666
667
668
669

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

        i_conf = 0
670
        for calc in single_conf:
671
            system_index = calc["single_configuration_calculation_to_system_ref"]
672
673
674
            system = systems[system_index]

            if i_conf == 0 or i_conf == 2 or i_conf == 4:
675
                with self.assertRaises(KeyError):
676
                    pos = system["atom_positions"]
677
            else:
678
                pos = system["atom_positions"]
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
                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
723
724

class TestMD(unittest.TestCase):
725
726
    """Molecular dynamics tests.
    """
Lauri Himanen's avatar
Lauri Himanen committed
727
728
    @classmethod
    def setUpClass(cls):
729
        cls.results = get_result("md/nve")
Lauri Himanen's avatar
Lauri Himanen committed
730
731
732
733
734
735
736
737
738
        cls.temp = convert_unit(
            np.array([
                300.000000000,
                275.075405378,
                235.091633019,
                202.752506973,
                192.266488819,
                201.629598676,
                218.299664775,
739
                230.324748557,
Lauri Himanen's avatar
Lauri Himanen committed
740
741
742
743
744
745
746
747
748
749
750
751
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
                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
794

795
796
797
798
799
    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))

800
801
802
803
804
805
806
807
808
809
810
811
812
    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
813
    def test_ensemble_type(self):
814
        result = self.results["ensemble_type"]
Lauri Himanen's avatar
Lauri Himanen committed
815
816
        self.assertEqual(result, "NVE")

817
818
819
820
    def test_sampling_method(self):
        result = self.results["sampling_method"]
        self.assertEqual(result, "molecular_dynamics")

Lauri Himanen's avatar
Lauri Himanen committed
821
822
823
824
    def test_number_of_frames_in_sequence(self):
        result = self.results["number_of_frames_in_sequence"]
        self.assertEqual(result, 11)

825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
    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"
        )
849
850
        self.assertTrue(np.array_equal(result[0, :], expected_start))
        self.assertTrue(np.array_equal(result[-1, :], expected_end))
851

852
853
854
855
856
857
858
859
860
861
862
    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],
            ]),
863
            "bohr*(hbar/hartree)^-1"
864
865
866
867
868
869
870
871
872
873
        )
        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],
            ]),
874
            "bohr*(hbar/hartree)^-1"
875
876
        )

877
878
        self.assertTrue(np.array_equal(result[0, :], expected_start))
        self.assertTrue(np.array_equal(result[-1, :], expected_end))
879

880
881
    def test_frame_sequence_potential_energy(self):
        result = self.results["frame_sequence_potential_energy"]
Lauri Himanen's avatar
Lauri Himanen committed
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
        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"]
898
899
        expected_result = convert_unit(
            np.array([
Lauri Himanen's avatar
Lauri Himanen committed
900
901
902
903
904
905
906
907
908
909
910
                0.000000,
                0.500000,
                1.000000,
                1.500000,
                2.000000,
                2.500000,
                3.000000,
                3.500000,
                4.000000,
                4.500000,
                5.000000,
911
            ]),
Lauri Himanen's avatar
Lauri Himanen committed
912
            "fs"
913
914
915
        )
        self.assertTrue(np.array_equal(result, expected_result))

Lauri Himanen's avatar
Lauri Himanen committed
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
    def test_frame_sequence_potential_energy_stats(self):
        result = self.results["frame_sequence_potential_energy_stats"]
        expected_result = np.array([self.pot.mean(), self.pot.std()])
        self.assertTrue(np.array_equal(result, expected_result))

    def test_frame_sequence_kinetic_energy_stats(self):
        result = self.results["frame_sequence_kinetic_energy_stats"]
        expected_result = np.array([self.kin.mean(), self.kin.std()])
        self.assertTrue(np.array_equal(result, expected_result))

    def test_frame_sequence_conserved_quantity_stats(self):
        result = self.results["frame_sequence_conserved_quantity_stats"]
        expected_result = np.array([self.cons.mean(), self.cons.std()])
        self.assertTrue(np.array_equal(result, expected_result))

    def test_frame_sequence_temperature_stats(self):
        result = self.results["frame_sequence_temperature_stats"]
        expected_result = np.array([self.temp.mean(), self.temp.std()])
        self.assertTrue(np.array_equal(result, expected_result))


class TestMDEnsembles(unittest.TestCase):
938
939
    """Different ensembles in MD.
    """
Lauri Himanen's avatar
Lauri Himanen committed
940
941
942
    @classmethod
    def setUpClass(cls):
        cls.pressure = convert_unit(
943
            np.array([
944
945
                -0.192828092559E+04,
                -0.145371071470E+04,
Lauri Himanen's avatar
Lauri Himanen committed
946
                -0.210098903760E+03,
947
948
                0.167260570313E+04,
                0.395562042841E+04,
Lauri Himanen's avatar
Lauri Himanen committed
949
                0.630374855942E+04,
950
951
                0.836906136786E+04,
                0.983216022830E+04,
Lauri Himanen's avatar
Lauri Himanen committed
952
                0.104711540465E+05,
953
                0.102444821550E+05,
Lauri Himanen's avatar
Lauri Himanen committed
954
                0.931695792434E+04,
955
            ]),
Lauri Himanen's avatar
Lauri Himanen committed
956
            "bar"
957
        )
Lauri Himanen's avatar
Lauri Himanen committed
958
959

    def test_nvt(self):
960
        results = get_result("md/nvt")
Lauri Himanen's avatar
Lauri Himanen committed
961
962
963
964
        ensemble = results["ensemble_type"]
        self.assertEqual(ensemble, "NVT")

    def test_npt(self):
965
        results = get_result("md/npt")
Lauri Himanen's avatar
Lauri Himanen committed
966
967
968
969
970
971
972
973
974
        ensemble = results["ensemble_type"]
        self.assertEqual(ensemble, "NPT")

        pressure = results["frame_sequence_pressure"]
        self.assertTrue(np.array_equal(pressure, self.pressure))

        pressure_stats = results["frame_sequence_pressure_stats"]
        expected_pressure_stats = np.array([self.pressure.mean(), self.pressure.std()])
        self.assertTrue(np.array_equal(pressure_stats, expected_pressure_stats))
Lauri Himanen's avatar
Lauri Himanen committed
975

976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
        simulation_cell = results["simulation_cell"]
        expected_cell_start = convert_unit(
            np.array(
                [[
                    6.0000000000,
                    0.0000000000,
                    0.0000000000,
                ], [
                    0.0000000000,
                    6.0000000000,
                    0.0000000000,
                ], [
                    0.0000000000,
                    0.0000000000,
                    6.0000000000,
                ]]),
            "angstrom"
        )
        expected_cell_end = convert_unit(
            np.array(
                [[
                    5.9960617905,
                    -0.0068118798,
                    -0.0102043036,
                ], [