1# -*- coding: utf-8 -*-
2#
3# Copyright (c) 2018, the cclib development team
4#
5# This file is part of cclib (http://cclib.github.io) and is distributed under
6# the terms of the BSD 3-Clause License.
7
8"""Test logfiles with vibration output in cclib"""
9
10import os
11import unittest
12
13from skip import skipForParser, skipForLogfile
14
15
16__filedir__ = os.path.realpath(os.path.dirname(__file__))
17
18
19class GenericIRTest(unittest.TestCase):
20    """Generic vibrational frequency unittest"""
21
22    # Unit tests should normally give this value for the largest IR intensity.
23    max_IR_intensity = 100
24
25    # Unit tests may give these values for the largest force constant and reduced mass, respectively.
26    max_force_constant = 10.0
27    max_reduced_mass = 6.9
28
29    # reference zero-point correction from Gaussian 16 dvb_ir.out
30    zpve = 0.1771
31
32    def setUp(self):
33        """Initialize the number of vibrational frequencies on a per molecule basis"""
34        self.numvib = 3*len(self.data.atomnos) - 6
35
36    def testbasics(self):
37        """Are basic attributes correct?"""
38        self.assertEqual(self.data.natom, 20)
39
40    def testvibdisps(self):
41        """Are the dimensions of vibdisps consistent with numvib x N x 3"""
42        self.assertEqual(len(self.data.vibfreqs), self.numvib)
43        self.assertEqual(self.data.vibdisps.shape,
44                         (self.numvib, len(self.data.atomnos), 3))
45
46    def testlengths(self):
47        """Are the lengths of vibfreqs and vibirs (and if present, vibsyms, vibfconnsts and vibrmasses) correct?"""
48        self.assertEqual(len(self.data.vibfreqs), self.numvib)
49        if hasattr(self.data, 'vibirs'):
50            self.assertEqual(len(self.data.vibirs), self.numvib)
51        if hasattr(self.data, 'vibsyms'):
52            self.assertEqual(len(self.data.vibsyms), self.numvib)
53        if hasattr(self.data, 'vibfconsts'):
54            self.assertEqual(len(self.data.vibfconsts), self.numvib)
55        if hasattr(self.data, 'vibrmasses'):
56            self.assertEqual(len(self.data.vibrmasses), self.numvib)
57
58    def testfreqval(self):
59        """Is the highest freq value 3630 +/- 200 wavenumber?"""
60        self.assertAlmostEqual(max(self.data.vibfreqs), 3630, delta=200)
61
62    @skipForParser('Psi4', 'Psi cannot print IR intensities')
63    def testirintens(self):
64        """Is the maximum IR intensity 100 +/- 10 km/mol?"""
65        self.assertAlmostEqual(max(self.data.vibirs), self.max_IR_intensity, delta=10)
66
67    @skipForParser('ADF', 'ADF cannot print force constants')
68    @skipForParser('DALTON', 'DALTON cannot print force constants')
69    @skipForParser('GAMESS', 'GAMESS-US cannot print force constants')
70    @skipForParser('GAMESSUK', 'GAMESS-UK cannot print force constants')
71    @skipForParser('Molcas', 'Molcas cannot print force constants')
72    @skipForParser('Molpro', 'Molpro cannot print force constants')
73    @skipForParser('ORCA', 'ORCA cannot print force constants')
74    @skipForParser('Turbomole', 'Turbomole cannot print force constants')
75    @skipForLogfile('Jaguar/Jaguar4.2', 'Data file does not contain force constants')
76    @skipForLogfile('Psi4/Psi4-1.0', 'Data file contains vibrational info with cartesian coordinates')
77    def testvibfconsts(self):
78        """Is the maximum force constant 10. +/- 0.1 mDyn/angstrom?"""
79        self.assertAlmostEqual(max(self.data.vibfconsts), self.max_force_constant, delta=0.1)
80
81    @skipForParser('ADF', 'ADF cannot print reduced masses')
82    @skipForParser('DALTON', 'DALTON cannot print reduced masses')
83    @skipForParser('GAMESSUK', 'GAMESSUK cannot print reduced masses')
84    @skipForParser('Molpro', 'Molpro cannot print reduced masses')
85    @skipForParser('ORCA', 'ORCA cannot print reduced masses')
86    @skipForLogfile('GAMESS/PCGAMESS', 'Data file does not contain reduced masses')
87    @skipForLogfile('Psi4/Psi4-1.0', 'Data file does not contain reduced masses')
88    def testvibrmasses(self):
89        """Is the maximum reduced mass 6.9 +/- 0.1 daltons?"""
90        self.assertAlmostEqual(max(self.data.vibrmasses), self.max_reduced_mass, delta=0.1)
91
92    @skipForParser('Psi3', 'not implemented yet')
93    def testzeropointcorrection(self):
94        """Is the zero-point correction correct?"""
95        self.assertAlmostEqual(self.data.zpve, self.zpve, delta=1.0e-3)
96
97
98class ADFIRTest(GenericIRTest):
99    """Customized vibrational frequency unittest"""
100
101    # ???
102    def testzeropointcorrection(self):
103        """Is the zero-point correction correct?"""
104        self.assertAlmostEqual(self.data.zpve, self.zpve, delta=1.0e-2)
105
106
107class FireflyIRTest(GenericIRTest):
108    """Customized vibrational frequency unittest"""
109
110    max_IR_intensity = 135
111    # ???
112    zpve = 0.1935
113
114
115class GaussianIRTest(GenericIRTest):
116    """Customized vibrational frequency unittest"""
117
118    def testvibsyms(self):
119        """Is the length of vibsyms correct?"""
120        self.assertEqual(len(self.data.vibsyms), self.numvib)
121
122    def testzeropointcorrection(self):
123        # reference zero-point correction from dvb_ir.out
124        zpve = 0.1771
125        """Is the zero-point correction correct?"""
126        self.assertAlmostEqual(self.data.zpve, zpve, delta=0.001)
127
128    entropy_places = 6
129    enthalpy_places = 3
130    freeenergy_places = 3
131
132    def testtemperature(self):
133        """Is the temperature 298.15 K?"""
134        self.assertAlmostEqual(298.15, self.data.temperature)
135
136    def testpressure(self):
137        """Is the pressure 1 atm?"""
138        self.assertAlmostEqual(1, self.data.pressure)
139
140    def testentropy(self):
141         """Is the entropy reasonable"""
142         self.assertAlmostEqual(0.0001462623335480945, self.data.entropy, self.entropy_places)
143
144    def testenthalpy(self):
145         """Is the enthalpy reasonable"""
146         self.assertAlmostEqual(-382.12130688525264, self.data.enthalpy, self.enthalpy_places)
147
148    def testfreeenergy(self):
149         """Is the freeenergy reasonable"""
150         self.assertAlmostEqual(-382.164915, self.data.freeenergy, self.freeenergy_places)
151
152    def testfreeenergyconsistency(self):
153        """Does G = H - TS hold"""
154        self.assertAlmostEqual(self.data.enthalpy - self.data.temperature * self.data.entropy, self.data.freeenergy, self.freeenergy_places)
155
156class JaguarIRTest(GenericIRTest):
157    """Customized vibrational frequency unittest"""
158
159    # Jagur outputs vibrational info with cartesian coordinates
160    max_force_constant = 3.7
161    max_reduced_mass = 2.3
162
163    def testvibsyms(self):
164        """Is the length of vibsyms correct?"""
165        self.assertEqual(len(self.data.vibsyms), self.numvib)
166
167
168class MolcasIRTest(GenericIRTest):
169    """Customized vibrational frequency unittest"""
170
171    max_IR_intensity = 65
172    zpve = 0.1783
173
174    entropy_places = 6
175    enthalpy_places = 3
176    freeenergy_places = 3
177
178    def testtemperature(self):
179        """Is the temperature 298.15 K?"""
180        self.assertAlmostEqual(298.15, self.data.temperature)
181
182    def testpressure(self):
183        """Is the pressure 1 atm?"""
184        self.assertAlmostEqual(1, self.data.pressure)
185
186    def testentropy(self):
187         """Is the entropy reasonable"""
188         self.assertAlmostEqual(0.00013403320476271246, self.data.entropy, self.entropy_places)
189
190    def testenthalpy(self):
191         """Is the enthalpy reasonable"""
192         self.assertAlmostEqual(-382.11385, self.data.enthalpy, self.enthalpy_places)
193
194    def testfreeenergy(self):
195         """Is the freeenergy reasonable"""
196         self.assertAlmostEqual(-382.153812, self.data.freeenergy, self.freeenergy_places)
197
198    def testfreeenergyconsistency(self):
199        """Does G = H - TS hold"""
200        self.assertAlmostEqual(self.data.enthalpy - self.data.temperature * self.data.entropy, self.data.freeenergy, self.freeenergy_places)
201
202
203class OrcaIRTest(GenericIRTest):
204    """Customized vibrational frequency unittest"""
205
206    # ORCA has a bug in the intensities for version < 4.0
207    max_IR_intensity = 215
208    zpve = 0.1921
209
210    enthalpy_places = 3
211    entropy_places = 6
212    freeenergy_places = 3
213
214    def testtemperature(self):
215        """Is the temperature 298.15 K?"""
216        self.assertAlmostEqual(298.15, self.data.temperature)
217
218    def testpressure(self):
219        """Is the pressure 1 atm?"""
220        self.assertAlmostEqual(1, self.data.pressure)
221
222    def testenthalpy(self):
223         """Is the enthalpy reasonable"""
224         self.assertAlmostEqual(-381.85224835, self.data.enthalpy, self.enthalpy_places)
225
226    def testentropy(self):
227         """Is the entropy reasonable"""
228         self.assertAlmostEqual(0.00012080325339594164, self.data.entropy, self.entropy_places)
229
230    def testfreeenergy(self):
231         """Is the freeenergy reasonable"""
232         self.assertAlmostEqual(-381.88826585, self.data.freeenergy, self.freeenergy_places)
233
234    def testfreeenergyconsistency(self):
235        """Does G = H - TS hold"""
236        self.assertAlmostEqual(self.data.enthalpy - self.data.temperature * self.data.entropy, self.data.freeenergy, self.freeenergy_places)
237
238
239class QChemIRTest(GenericIRTest):
240    """Customized vibrational frequency unittest"""
241
242    enthalpy_places = 3
243    entropy_places = 6
244    freeenergy_places = 3
245
246    def testtemperature(self):
247        """Is the temperature 298.15 K?"""
248        self.assertEqual(298.15, self.data.temperature)
249
250    def testpressure(self):
251        """Is the pressure 1 atm?"""
252        self.assertAlmostEqual(1, self.data.pressure)
253
254    def testenthalpy(self):
255         """Is the enthalpy reasonable"""
256         self.assertAlmostEqual(0.1871270552135131, self.data.enthalpy, self.enthalpy_places)
257
258    def testentropy(self):
259         """Is the entropy reasonable"""
260         self.assertAlmostEqual(0.00014667348271900577, self.data.entropy, self.entropy_places)
261
262    def testfreeenergy(self):
263         """Is the freeenergy reasonable"""
264         self.assertAlmostEqual(0.14339635634084155, self.data.freeenergy, self.freeenergy_places)
265
266    # Molecular mass of DVB in mD.
267    molecularmass = 130078.25
268
269    def testatommasses(self):
270        """Do the atom masses sum up to the molecular mass (130078.25+-0.1mD)?"""
271        mm = 1000*sum(self.data.atommasses)
272        self.assertAlmostEqual(mm, 130078.25, delta=0.1, msg = "Molecule mass: %f not 130078 +- 0.1mD" % mm)
273
274    def testhessian(self):
275        """Do the frequencies from the Hessian match the printed frequencies?"""
276
277    def testfreeenergyconsistency(self):
278        """Does G = H - TS hold"""
279        self.assertAlmostEqual(self.data.enthalpy - self.data.temperature * self.data.entropy, self.data.freeenergy, self.freeenergy_places)
280
281
282class GamessIRTest(GenericIRTest):
283    """Customized vibrational frequency unittest"""
284    # Molecular mass of DVB in mD.
285    molecularmass = 130078.25
286    enthalpy_places = 3
287    entropy_places = 6
288    freeenergy_places = 3
289
290    def testatommasses(self):
291        """Do the atom masses sum up to the molecular mass (130078.25+-0.1mD)?"""
292        mm = 1000*sum(self.data.atommasses)
293        self.assertAlmostEqual(mm, 130078.25, delta=0.1, msg = "Molecule mass: %f not 130078 +- 0.1mD" % mm)
294
295    def testtemperature(self):
296        """Is the temperature 298.15 K?"""
297        self.assertAlmostEqual(298.15, self.data.temperature)
298
299    def testpressure(self):
300        """Is the pressure 1 atm?"""
301        self.assertAlmostEqual(1, self.data.pressure)
302
303    def testenthalpy(self):
304         """Is the enthalpy reasonable"""
305         self.assertAlmostEqual(-381.86372805188300, self.data.enthalpy, self.enthalpy_places)
306
307    def testentropy(self):
308         """Is the entropy reasonable"""
309         self.assertAlmostEqual(0.00014875961938, self.data.entropy, self.entropy_places)
310
311    def testfreeenergy(self):
312         """Is the freeenergy reasonable"""
313         self.assertAlmostEqual(-381.90808120060200, self.data.freeenergy, self.freeenergy_places)
314
315    def testfreeenergyconsistency(self):
316        """Does G = H - TS hold"""
317        self.assertAlmostEqual(self.data.enthalpy - self.data.temperature * self.data.entropy, self.data.freeenergy, self.freeenergy_places)
318
319
320class Psi4IRTest(GenericIRTest):
321    """Customized vibrational frequency unittest"""
322
323    # RHF is used for Psi4 IR test data instead of B3LYP
324    max_force_constant = 9.37
325    zpve = 0.1917
326
327
328class TurbomoleIRTest(GenericIRTest):
329    """Customized vibrational frequency unittest"""
330
331    # ???
332    zpve = 0.1725
333
334
335class GenericIRimgTest(unittest.TestCase):
336    """Generic imaginary vibrational frequency unittest"""
337
338    def setUp(self):
339        """Initialize the number of vibrational frequencies on a per molecule basis"""
340        self.numvib = 3*len(self.data.atomnos) - 6
341
342    def testvibdisps(self):
343        """Are the dimensions of vibdisps consistent with numvib x N x 3"""
344        self.assertEqual(self.data.vibdisps.shape,
345                         (self.numvib, len(self.data.atomnos), 3))
346
347    def testlengths(self):
348        """Are the lengths of vibfreqs and vibirs correct?"""
349        self.assertEqual(len(self.data.vibfreqs), self.numvib)
350        self.assertEqual(len(self.data.vibirs), self.numvib)
351
352    def testfreqval(self):
353        """Is the lowest freq value negative?"""
354        self.assertTrue(self.data.vibfreqs[0] < 0)
355
356##    def testmaxvibdisps(self):
357##        """What is the maximum value of displacement for a H vs a C?"""
358##        Cvibdisps = compress(self.data.atomnos==6, self.data.vibdisps, 1)
359##        Hvibdisps = compress(self.data.atomnos==1, self.data.vibdisps, 1)
360##        self.assertEqual(max(abs(Cvibdisps).flat), 1.0)
361
362
363class GenericRamanTest(unittest.TestCase):
364    """Generic Raman unittest"""
365
366    # This value is in amu.
367    max_raman_intensity = 575
368
369    def setUp(self):
370        """Initialize the number of vibrational frequencies on a per molecule basis"""
371        self.numvib = 3*len(self.data.atomnos) - 6
372
373    def testlengths(self):
374        """Is the length of vibramans correct?"""
375        self.assertEqual(len(self.data.vibramans), self.numvib)
376
377    # The tolerance for this number has been increased, since ORCA
378    # failed to make it inside +/-5, but it would be nice in the future
379    # to determine is it's not too much work whether this is due to
380    # algorithmic differences, or to differences in the input basis set
381    # or coordinates. The first would be OK, but in the second case the
382    # unit test jobs should be made more comparable. With cclib, we first
383    # of all want to succeed in parsing, but would also like to remain
384    # as comparable between programs as possible (for these tests).
385    # Note also that this value is adjusted for Gaussian and DALTON - why?
386    def testramanintens(self):
387        """Is the maximum Raman intensity correct?"""
388        self.assertAlmostEqual(max(self.data.vibramans), self.max_raman_intensity, delta=8)
389
390        # We used to test this, but it seems to vary wildly between
391        # programs... perhaps we could use it if we knew why...
392        #self.assertInside(self.data.vibramans[1], 2.6872, 0.0001)
393
394    def testvibdisps(self):
395        """Is the length and value of vibdisps correct?"""
396        assert hasattr(self.data, "vibdisps")
397        assert len(self.data.vibdisps) == 54
398
399
400class DALTONRamanTest(GenericRamanTest):
401    """Customized Raman unittest"""
402
403    max_raman_intensity = 745
404
405
406class GaussianRamanTest(GenericRamanTest):
407    """Customized Raman unittest"""
408
409    max_raman_intensity = 1066
410
411
412class OrcaRamanTest(GenericRamanTest):
413    """Customized Raman unittest"""
414
415    max_raman_intensity = 1048
416
417
418class QChemRamanTest(GenericRamanTest):
419    """Customized Raman unittest"""
420
421    max_raman_intensity = 588
422
423if __name__=="__main__":
424
425    import sys
426    sys.path.insert(1, os.path.join(__filedir__, ".."))
427
428    from test_data import DataSuite
429    suite = DataSuite(['vib'])
430    suite.testall()
431