forked from Jellby/Pegamoid
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpegamoid.py
executable file
·7413 lines (6955 loc) · 283 KB
/
pegamoid.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
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
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
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
671
672
673
674
675
676
677
678
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
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
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
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
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
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
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
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from __future__ import division, absolute_import, print_function
from builtins import bytes, dict, int, range, super
__name__ = 'Pegamoid'
__author__ = u'Ignacio Fdez. Galván'
__copyright__ = u'Copyright © 2018–2020'
__license__ = 'GPL v3.0'
__version__ = '2.6.1'
import sys
try:
from qtpy.QtCore import Qt, QObject, QThread, QEvent, QSettings
from qtpy.QtWidgets import *
from qtpy.QtGui import QPixmap, QIcon, QKeySequence, QColor, QPalette
import qtpy
v = qtpy.PYQT_VERSION
if (v is None):
v = qtpy.PYSIDE_VERSION
QtVersion = '{0} {1} (Qt {2})'.format(qtpy.API_NAME, v, qtpy.QT_VERSION)
if (qtpy.API in qtpy.PYSIDE2_API):
print('Unfortunately, VTK does not support PySide2 yet')
sys.exit(1)
except:
try:
from PyQt5.QtCore import Qt, QObject, QThread, QEvent, QSettings, PYQT_VERSION_STR, QT_VERSION_STR
from PyQt5.QtWidgets import *
from PyQt5.QtGui import QPixmap, QIcon, QKeySequence, QColor, QPalette
QtVersion = 'PyQt5 {0} (Qt {1})'.format(PYQT_VERSION_STR, QT_VERSION_STR)
except ImportError:
try:
from PyQt4.QtCore import Qt, QObject, QThread, QEvent, QSettings, PYQT_VERSION_STR, QT_VERSION_STR
from PyQt4.QtGui import *
QtVersion = 'PyQt4 {0} (Qt {1})'.format(PYQT_VERSION_STR, QT_VERSION_STR)
except ImportError:
print('{0} needs python Qt bindings.'.format(__name__))
print('Please install at least one of: PyQt4, PyQt5, PySide')
sys.exit(1)
import os
import vtk.qt
try:
if os.environ.get('PEGAMOID_NO_QGL', None):
raise ImportError
vtk.qt.QVTKRWIBase = 'QGLWidget'
from vtk.qt.QVTKRenderWindowInteractor import QVTKRenderWindowInteractor
QtVersion += ' with QtOpenGL'
except ImportError:
vtk.qt.QVTKRWIBase = 'QWidget'
from vtk.qt.QVTKRenderWindowInteractor import QVTKRenderWindowInteractor
import vtk
from vtk.util import numpy_support
import h5py
import numpy as np
from fractions import Fraction
import os.path
import codecs
import re
import struct
import traceback
import time
from copy import deepcopy
from socket import gethostname
from datetime import datetime
from tempfile import mkdtemp, TemporaryFile
from shutil import rmtree
from functools import partial
from collections import OrderedDict
try:
from itertools import zip_longest
except ImportError:
from itertools import izip_longest as zip_longest
try:
from ttfquery._scriptregistry import registry
except ImportError:
pass
icondata = codecs.decode(b'''
iVBORw0KGgoAAAANSUhEUgAAADAAAAAwCAYAAABXAvmHAAAABGdBTUEAALGPC/xhBQAAAAFzUkdC
AK7OHOkAAAAgY0hSTQAAeiYAAICEAAD6AAAAgOgAAHUwAADqYAAAOpgAABdwnLpRPAAAAAZiS0dE
AP8A/wD/oL2nkwAAAAlwSFlzAAAASAAAAEgARslrPgAAC31JREFUaN7tmWusXNdVx/9r7X0ec+bM
nTv36Vdsx4kNseu0xE1CadqGggCpEhKilaBSJb4UCZSvVAjxiCpRUD/QSAUVgYQKAhpRJPhQ5QsK
tCqRiJs4IbUdP5N7k2vfx9w775kzZ5+91+LDdRGIYDtc+gD5L62zpTl7pP2bNXvttdYG7ume7un/
tOjdfuGL+uuYwyzODl7B0A0pqABQCARBA8pQwrPHMPQwxBA+AGEqYCZAFRIA9UriFH4k8F0hvyVw
1wOGLzgdfbMUignq9H8HYOajCTQgy05HJ+uPxJocM+N4ny1qec1lplbFFHsDEwCVoEEqVOpQaemn
KIdTBYHEEkhBAAhCjApGHKxOkMhY0zCURtXzc35blv22NKXUrwG4ufmHozsC2DtNkIlCA97jd+Rz
1Va4whldAIdVCZN21XD9elQvMspcjLiKOQl1kBBIJFKVnGFZYYwlBjODDUAWqnFQSb34vJSyOZHJ
wtAP9036xf3l2/54tS5vUoSbd+OBOwJwSpBSj8CgASBRRV2D1qE6YmaXUKIZZZRSyjFib8kGAyME
Us6MKgS3Fs8EMgAihcYePnHq4qlOY1a2Xiuq2LkwkCp09XByv8Ghp2ew9vTg9uu73cv4PoP4sAFZ
OswpeUqo5BiOY/IUUYjYSoRILeyu0e5oYDQiqwaAAYNA+h0DSAEo6a6nCBQIFFjZg1HBoALh0OaX
xub2q7sLAIoI238xIYqwnxNynJCjhCpKyHPMYsiqhVFDBgZGDQyYDAwZEBiMGAwLBkGh2N2WqrsP
FQCiQADgFfBQVERUQbGUPmhrqnv8C5EB7BzHZGmBUnIUw3FCnhMEGxmxsMpkwGBl4t0RrAklEsFK
0FI9GB7KgJJCVHdjkShUBSIKEVEJguBVEVS00qBNEOoIGO1tDxjANCiGQcYxKk4oUIzACUnEkUYU
g8FKYCWQGjZqiPXq9OpyZrIHm9xcyihzMcUrlVbXAkIQiKqqBgQJGiRoCAE+BA1evQT1EK1QU685
BJt78wAT1MISIyZLgSwFikk4YokpRowIBgYEwDBrvzO059sXPtw4lP1obvOGpYgN2RqDagI577R6
FtCuh1evXj28VPBSoQqVVCJOVZ2qlhqrQ6q85yiku2cFgcC7AZAiqDEGESKKKEJEFkyMahzo7D+8
+pPmITmzFM1fz6mxU6esiilORHVfBX+qxPRTTt2fJZSU34Go1KlTpz54lRIkpZJM1Uihhva6iVUA
9QhQeAAEBpEhMDNZsmxh2cJSVq/pua9cPDrpTT8UN6NLM3bm9ZSSaxb2moG5Ysj8KxO/WGq5sFlu
nlGFeng4ODjsIgQn0KmyFmqkUJZCRQrdowcCIBNxGlBAkILAxGBDxhhYNrBkyFJEEXobg/eaEzSZ
WWisNanZjSgaKdR79UagLqiPCpm2up3+D8+kg5fjlkGlFUqdogwlpBAjEzVholEYqcpYpmDaowec
wrfFaaVdqTSGIALBMhljYYwlYyyskb5asbqYHop6rXS2TJFWBtYTKAikcijdUIejbtnpDN+YzA7W
R3UwKEC41JJ84U0YaRRGEoe+pmEgLvR1FAayNwCZAstP5UEdVmSiiXrEUIoYFBkykYWNLKxNksQY
w6WxnOYmz2IkCQFRgBgPb6da2h3d4fZqJy3OV5ERKwFiAjx78TZMNAp9SUJXMr8jue9I162H/p4B
qs2A8TkHmeiFMFArhdY1aA1AyuDUkEks2ThJE9tqzr5RbvjWxE8WlHQWoDwg1AudZl3tZBuDzZnB
S8X9YYW2Fk41Sxdc5NRF3vlYxpqFvuZVW2bcjTDrt+T60i/Xi9DXvQEAgO8IfE9e8zvS9TsyH8ba
9MHXBVInUJ1hMplKfubJ06uyRpP1tc3DQwwWnJYtr75ZatnsSne2e3VwpHgtLLf2Nc+Fuo8mOknH
fpxW/ZD7jjSrdphzN8Ji+aavuZvhhc0vjlBcrPYOEHqKT157+KZvh+fLN/1ytREWp+OyNdHJjFPX
DPCzQaR58PSyObb/yNmdV4f7N4qN/QWK+RJlq8Bktj/uL47OuVM1Si+d/tUH1gfj4cwgDOvTnput
NsKCWwvLbiUcnF72R8sV/3q54s+69XAriuNOZ+0d0umR4tKzbfi+vkXAhzim/ZSjlEYIHDMbGMtE
VoLG9/3QgcHGt9tZx3UOZwfSqYK4i262dm7zweo86g//wvF/wgGxvdDPx8Nx0234pWpNDpar4ej0
sj9eXPaJ35Cnuc7XZXB3BY25m0kyVtQeiobuZnhLCv0oBIsSiZTpVKvEcaCQBPgkcEiXDy/2tq52
5oc8aDWsT7e32nOdb7l9Rx87eK72fuu23fbMcDBqlevVsnsrHCrf9Eenl6sT00u+VW3I78tIn9Py
7hZ/1wDqAb8tOPR7zbXe30+vhZ48Hvp6zJehNtVpPLHjZGKLbEpFw9eqPM8SLL+y+ZEn7PSJR5C/
xzfo5vixqNMZdxZG7fG+6Yo/XF719xcX/fHi29WJ4nXvq3X5XOjLV8H/nrZ+F2piBvIPJCiv+8PR
Ev+iXTJPxPu4aZeNRAtc2hnjbd34eF2aJzbdqZOPpC6PsujlN/qDlxJcrcZiQkfialsivyVabYVu
1ZazoStflUIvfU+K+v+U69Vo3szwEa7Tfk5plhJkiMjMAY88tN98at9RM/UTra29HVZfvRm+UDnt
6lTHMtFeGOuWTHQDigHwbn7zd1lS3jZVKnRbirD9Hz984oMJfKWHGgYPxBP9SDnScRjIHw1fq750
8oTFxZXw/W2r3ElzvzODZ0rCZyEfa0n8Y82ji752cv7Pi/7ojc/87DP4KfrxH1yA4383DxmpjRb5
g3bJPNk63No4yAeWGtLIUyR/WUd2fopSn1n4yg8GwCe3fhrgGjpyg85vXjANzQ4Z4p+hiB7nGk0j
ExUNbqSLZvFgi+dsSslzBHpuS9qbYxkp3yryCRCBSolSrvpVzamGVw+u/c8BTq8toZCCDtiDScZZ
M0HcMmRmGTxDQEOBXKF1haZBJRGEKIjUVOUYMc8ZsDdixZKlmCNOOE0yruUJkoxA7QLTG2MZlRWq
SjTsFvRQB1BJoMIQjwxs38BsM5kbBKzlnHfaph2ef+w8Nq5svDPAxzYeRUyxaYetMwR+NOPseErJ
gRhx08DUGBQRseFbvYfdYg0gIjYwZGFhyKoBE+82VIhodwZUabczoXTLdtsruwUfAUQKZYWwQhkg
YlAAMBXVrZf+9MLXX/iN175gGlSEob5zFJrneRjwPsflUxOd7CcQq4KVlCKKQk55yLmBGqVkYZXA
nkAlgPEtmwAoADgCVUQkBGIDNgyTMLhmyNQNTMPA5EycESghENMuAIIGLlGaQid2ooX16nPnXHPn
ci9Chr+y82Y1DP07A8xyEwwzX8LF4zAZe/Woc10WeUnmeS7UqDYxZNoEWgPobQA3CbRNwIDABYOr
mGLfpKacit+rP5L8BLarK3jR/TO1ZZsqVAbQiMA1AjcMeMGADzGZoww6RuD9RNQAYASBSy2jAYbJ
xnCTi6LsRkcptymjXPlvzoEaZbAw9ZzySa51t2SX/ZJZGqWUvk2gixX8lVLL9YEMxznn+gf/+CeY
fPzuzoxbJgCqW57aAbDyx73PvlRqSU2aySxFBw3MSUP8PoP4gZjiRo7ctmZa2636N//mBjYHt93E
v7nzaUQUf3ggg48HDb2E4wsKvFJqufrM/F+Xn9h6En+7/I3vWq//851fw/uSh7HuN9KY4mOG+AMM
fj+De1D63Z9rfGIULy2iast/BfitzqcxVYc6ZY/3pLd/qtMXezrYTJHIl5ee+55fXHy5/3k4OM4p
vy8ie8TA/IsC7ucbv/LOqYRC8bI7h0fjM69MtPgWAHl26fnv283LLzU/AwDy252nVhWyeive3buS
uqd7uqd7uqf/X/o3WcSqZKuF0icAAAAASUVORK5CYII=
''', 'base64')
boxicondata = codecs.decode(b'''
iVBORw0KGgoAAAANSUhEUgAAAEAAAAApCAQAAAAvm+fHAAAABGdBTUEAALGPC/xhBQAAAAFzUkdC
AK7OHOkAAAAgY0hSTQAAeiYAAICEAAD6AAAAgOgAAHUwAADqYAAAOpgAABdwnLpRPAAAAAJiS0dE
AP+Hj8y/AAAACXBIWXMAAC38AAAt/AGuw+yYAAAGBklEQVRYw7XYe2zW1RkH8A+XchEITO5Mpitg
AkPlIsMtUJ2XlEF0QwkXCTjcFi6C2cAlyhwiGl0QMTqUWtkUzDBKQCbCRAdxdIEIyE1BjDC0JTIu
HQVboLT07A9ef/29Lb++Bd3z/nd+53nO857v+Z7n+R4uxrJ08mOTLLJdqUIvG+VKjf2fraE2fmC4
WZbZpViFEP1O22uJe/XS/NtetoGWuhvid161xeHYshVKBcFRxaoEwTn/sdZMg12uwTdfupF+pnhR
gYPORMueU2K3N80xXL6gwkR9/cZKhVFyJ33oOXd+U1gGOxjb5DJ7rfGU8QboIAtt/EuwxxWgie7G
yLdLWcrjrP1eN9l1Wlza1r8QWz444DED0hAe6pRgXo1z0sEt5njfsRQsVY5ab7abtdfwYhLItl+w
z4YYwoetM8tN2mmgoUWC4wZd0LulPqZY5t/ORju4U55RumlSvwQmqVJlqub6mWa5zyOES+3wvKkK
BavrPPVZrnKXBbb5Kjq4X3jT/fprVffyLbwrOKB7FCrbSC/YkTr5wTlVgjf00TIjmG3leNh7DjuX
guW/CjwhV6ckWHKcELxY43MD7dzkEesciRAutsHjbtMxI8LN9fYrf/WZ8uj22OMvxrm6NizPCMrk
JoRqaqIKwfGInqd8bJGxemiakdxdDDXXlsi30kF56ZO+Z6+gQJvETV0gOGG02z1rq5NRqCJvmW6g
1omeLWS71QxvR2AGwYn0aRNUCmbUwZB9gnVaoYHvGOQh7zgUIXzcJnMN00WjlEczXd1oijwFitKu
tfNwrE7H6m1BkZ51MiSYmjbWTE8TLLE3Cn/GXkv80j3me9d+palzE1Qptc9a84yzUlBmVDzUDYoF
r0TZZ2JITYSHmWuT49HtcS76v2cU2WChyXJ01Qyd7RJs1D4e5I+C0+5I/P85SgT5dZ761gaabrXT
gkqHbbbYDLm6aZFWpsaqEDwYd+3iI8EH2iYGny8oM6Qel9n1jgrW6K31BdNtaoXgS9fEB+9WIZiZ
GLRrBobE7feCs8Ykfu/viGCprHhOywWHXJvolIkh1dbWZsEunRNnPCood1ftnF6L51TjNluVgSHV
9jOnBU8mfu9gm+BDHeKDswXlRiQ6nWfI4kSGVFtjiwXFbkicMUK5YHbtnLal55RmT2ZgSLX1clCw
KrFaZnlNcES/2jk9mhi0c0aGVNsDgkq/SPx+rUOC5fHakWWp4Ij+iU6ZGFJt59u1T3RNnDFTUOHu
+NA1vqyZU1pL3tvfBUddX48EhigTzK+DIR/UZsiDggpja9SubnLNsMRmh1UKTltjhoFa19F4N5Qv
KDE4ccYdtRnS3kbBTp1TtSvHZAttqFG7vq53JTZ5yjDfvSAfejggWJvYDTfySm2GjFYmWGm8p62t
Vbv2W+tp97jX4rR696lXL6CHpgmqTEz8/z0V1WbIakFQXqt25cVqV7yj2Rird4e84yGDUnqolfWC
fb6fmMAMQaUJ6YMnYv1JqVWmu0V2jdqVXu9+aLq3FKmM9NBWz7rdaCcFf0r0bKPgQgzJczAKVW6r
uYbGOhqJ3WEPY73kY6eifStJCbZmCT65F2ZIE1cb5892R6HKfWapX+vtsoy6uaPbPO6fMT10xDqP
+Il2NXbia4bkJIXqJNcTCtL00HsedqO2GRVvS328IaiKTlKpHRYaKTsqb90zMCRlrfR3vxW+iPTQ
V7ZZYISrEislXGa1oNB9Ftge9b0VPrfCNP00N1WVKpPqJ8+ayDZKnp1RqLMOWOY+fRP00CDHBS9p
mBIxs/wjpoeKbbBPsE/2xb2MtHez2dY7GiF8zPvmuFWHGq3WPEGZn6Z1EAM85kCa0n7+0p4vWrjO
JK/bH1O8H8k3RveUtLrCnqhdy9LRAOPNs8an0ZtBEBxMUNT1tMauNNxzMT1UodDf/FZfk1QI8t1p
jpV2K6lxrRXIM0W/erQz9Xi+uNxgM62N6aFix1InP/6OdNgWSzwgV3ctv42Xo5o9Yi8TLPGJ07FN
rnDMTsv8wc/10ubiXkYuFZaRXlao1HaLTPQjneqkay37H9gL/mliu8SBAAAAAElFTkSuQmCC
''', 'base64')
# For CIELab <-> RGB conversion: http://colorizer.org/
background_color = {
'F': (0.449, 0.448, 0.646), # CIELab(50,12,-27)
'I': (0.711, 0.704, 0.918), # CIELab(75,12,-27) Hue ~ 240
'1': (0.605, 0.759, 0.682), # CIELab(75,-17,5.5) Hue ~ 150
'2': (0.569, 0.775, 0.569), # CIELab(75,-27.5,21) Hue ~ 120
'3': (0.654, 0.757, 0.548), # CIELab(75,-18.5,24) Hue ~ 90
'S': (0.883, 0.673, 0.670), # CIELab(75,19.5,8) Hue ~ 0
'D': (0.609, 0.417, 0.416), # CIELab(50,19.5,8)
'?': (0.466, 0.466, 0.466) # CIELab(50,0,0)
}
angstrom = 1.88972612462 # = 1 / 0.529177210904
# Atomic radii are only available up to Z=109 in VTK, apparently
maxZ = 109
# Symmetry multiplication table
symmult = np.array([
[0, 1, 2, 3, 4, 5, 6, 7],
[1, 0, 3, 2, 5, 4, 7, 6],
[2, 3, 0, 1, 6, 7, 4, 5],
[3, 2, 1, 0, 7, 6, 5, 4],
[4, 5, 6, 7, 0, 1, 2, 3],
[5, 4, 7, 6, 1, 0, 3, 2],
[6, 7, 4, 5, 2, 3, 0, 1],
[7, 6, 5, 4, 3, 2, 1, 0],
])
#===============================================================================
# Class for orbitals defined in term of basis functions, which can be computed
# at arbitrary points in space.
# The basis functions and orbitals can be read from HDF5 and Molden formats.
# Orbital coefficients can be read from InpOrb format if an HDF5 file has been
# read before.
class Orbitals(object):
def __init__(self, orbfile, ftype):
self.inporb = None
self.file = orbfile
self.type = ftype
self.eps = np.finfo(np.float).eps
self.wf = 'SCF'
if (self.type == 'hdf5'):
self.inporb = 'gen'
self.h5file = self.file
self.read_h5_basis()
self.read_h5_MO()
self.get_orbitals('State', 0)
elif (self.type == 'molden'):
self.read_molden_basis()
self.read_molden_MO()
self.normalize_rad()
# Read basis set from an HDF5 file
def read_h5_basis(self):
with h5py.File(self.file, 'r') as f:
sym = f.attrs['NSYM']
self.N_bas = f.attrs['NBAS']
self.irrep = [i.decode('ascii').strip() for i in f.attrs['IRREP_LABELS']]
# First read the centers and their properties
if (sym > 1):
labels = f['DESYM_CENTER_LABELS'][:]
try:
charges = f['DESYM_CENTER_ATNUMS'][:]
except KeyError:
charges = f['DESYM_CENTER_CHARGES'][:]
coords = f['DESYM_CENTER_COORDINATES'][:]
self.mat = np.reshape(f['DESYM_MATRIX'][:], (sum(self.N_bas), sum(self.N_bas))).T
else:
labels = f['CENTER_LABELS'][:]
try:
charges = f['CENTER_ATNUMS'][:]
except KeyError:
charges = f['CENTER_CHARGES'][:]
coords = f['CENTER_COORDINATES'][:]
charges = charges.astype(np.int).clip(0, maxZ)
self.centers = [{'name':str(l.decode('ascii')).strip(), 'Z':q, 'xyz':x} for l,q,x in zip(labels, charges, coords)]
self.geomcenter = (np.amin(coords, axis=0) + np.amax(coords, axis=0))/2
# Then read the primitives and assign them to the centers
prims = f['PRIMITIVES'][:] # (exponent, coefficient)
prids = f['PRIMITIVE_IDS'][:] # (center, l, shell)
# The basis_id contains negative l if the shell is Cartesian
if (sym > 1):
basis_function_ids = 'DESYM_BASIS_FUNCTION_IDS'
else:
basis_function_ids = 'BASIS_FUNCTION_IDS'
bf_id = np.rec.fromrecords(np.insert(f[basis_function_ids][:], 4, -1, axis=1), names='c, s, l, m, tl') # (center, shell, l, m, true-l)
bf_cart = set([(b['c'], b['l'], b['s']) for b in bf_id if (b['l'] < 0)])
for i,c in enumerate(self.centers):
c['bf_ids'] = np.where(bf_id['c'] == i+1)[0].tolist()
# Add contaminants, which are found as lower l basis functions after higher l ones
# The "tl" field means the "l" from which exponents and coefficients are to be taken, or "true l"
ii = [sum(self.N_bas[:i]) for i in range(len(self.N_bas))]
if (sym > 1):
sbf_id = np.rec.fromrecords(np.insert(f['BASIS_FUNCTION_IDS'][:], 4, -1, axis=1), names='c, s, l, m, tl')
else:
sbf_id = bf_id
for i,nb in zip(ii, self.N_bas):
prev = {'c': -1, 'l': -1, 's': -1, 'tl': 0, 'm': None}
for b in sbf_id[i:i+nb]:
if (b['c']==prev['c'] and abs(b['l']) < abs(prev['tl'])):
b['tl'] = prev['tl']
else:
b['tl'] = b['l']
prev = {n: b[n] for n in b.dtype.names}
if (sym > 1):
for i,b in enumerate(self.mat.T):
bb = np.array(sbf_id[i])
nz = np.nonzero(b)[0]
assert (np.all(bf_id[nz][['l','s','m']] == bb[['l','s','m']]))
for j in nz:
bf_id[j]['tl'] = bb['tl']
# Workaround for bug in HDF5 files where p-type contaminants did all have m=0
p_shells, p0_counts = np.unique([np.array(b)[['c','s','tl']].copy() for b in bf_id if ((b['l']==1) and (b['m']==0))], return_counts=True)
if (np.any(p0_counts > 1)):
if (sym > 1):
# can't fix it with symmetry
error = 'Bad m for p contaminants. The file could have been created by a buggy or unsupported OpenMolcas version'
raise Exception(error)
else:
m = -1
for i in np.where(bf_id['l']==1)[0]:
bi = list(p_shells).index(bf_id[i][['c','s','tl']])
if (p0_counts[bi] > 1):
bf_id[i]['m'] = m
m += 1
if (m > 1):
m = -1
# Count the number of m per basis to make sure it matches with the expected type
counts = {}
for b in bf_id:
key = (b['c'], b['l'], b['s'], b['tl'])
counts[key] = counts.get(key, 0)+1
for f,n in counts.items():
l = f[1]
if (((l >= 0) and (n != 2*l+1)) or ((l < 0) and (n != (-l+1)*(-l+2)/2))):
error = 'Inconsistent basis function IDs. The file could have been created by a buggy or unsupported OpenMolcas version'
raise Exception(error)
# Maximum angular momentum in the whole basis set,
maxl = max([p[1] for p in prids])
for i,c in enumerate(self.centers):
c['basis'] = []
c['cart'] = {}
for l in range(maxl+1):
ll = []
# number of shells for this l and center
maxshell = max([0] + [p[2] for p in prids if ((p[0] == i+1) and (p[1] == l))])
for s in range(maxshell):
# find out if this is a Cartesian shell (if the l is negative)
# note that Cartesian shells never have (nor are) contaminants,
# and since contaminants come after regular shells,
# it should be safe to save just l and s
if ((i+1, -l, s+1) in bf_cart):
c['cart'][(l, s)] = True
# get exponents and coefficients
ll.append([0, [pp.tolist() for p,pp in zip(prids, prims) if ((p[0] == i+1) and (p[1] == l) and (p[2] == s+1))]])
c['basis'].append(ll)
# Add contaminant shells, that is, additional shells for lower l, with exponents and coefficients
# from a higher l, and with some power of r**2
for l in range(maxl-1):
# find basis functions for this center and l, where l != tl
cont = [(b['l'],b['tl'],b['s']) for b in bf_id[np.logical_and(bf_id['c']==i+1, bf_id['l']==l)] if (b['l'] != b['tl'])]
# get a sorted unique set
cont = sorted(set(cont))
# copy the exponents and coefficients from the higher l and set the power of r**2
for j in cont:
new = deepcopy(c['basis'][j[1]][j[2]-1])
new[0] = (j[1]-j[0])//2
c['basis'][l].append(new)
# At this point each center[i]['basis'] is a list of maxl items, one for each value of l,
# each item is a list of shells,
# each item is [power of r**2, primitives],
# each "primitives" is a list of [exponent, coefficient]
# Now get the indices for sorting all the basis functions (2l+1 or (l+1)(l+2)/2 for each shell)
# by center, l, m, "true l", shell
# To get the correct sorting for Cartesian shells, invert l
for b in bf_id:
if (b['l'] < 0):
b['l'] *= -1
self.bf_sort = np.argsort(bf_id, order=('c', 'l', 'm', 'tl', 's'))
# And sph_c can be computed
self.set_sph_c(maxl)
# center of atoms with basis
nb = [isEmpty(c['basis']) for c in self.centers]
if (any(nb) and not all(nb)):
xyz = [c['xyz'] for i,c in enumerate(self.centers) if not nb[i]]
self.geomcenter = (np.amin(xyz, axis=0) + np.amax(xyz, axis=0))/2
# Reading the basis set invalidates the orbitals, if any
self.base_MO = None
self.MO = None
self.MO_a = None
self.MO_b = None
self.current_orbs = None
# Read molecular orbitals from an HDF5 file
def read_h5_MO(self):
with h5py.File(self.file, 'r') as f:
# Read the orbital properties
if ('MO_ENERGIES' in f):
mo_en = f['MO_ENERGIES'][:]
mo_oc = f['MO_OCCUPATIONS'][:]
mo_cf = f['MO_VECTORS'][:]
if ('MO_TYPEINDICES' in f):
mo_ti = f['MO_TYPEINDICES'][:]
else:
mo_ti = [b'?' for i in mo_oc]
else:
mo_en = []
mo_oc = []
mo_cf = []
mo_ti = []
if ('MO_ALPHA_ENERGIES' in f):
mo_en_a = f['MO_ALPHA_ENERGIES'][:]
mo_oc_a = f['MO_ALPHA_OCCUPATIONS'][:]
mo_cf_a = f['MO_ALPHA_VECTORS'][:]
if ('MO_ALPHA_TYPEINDICES' in f):
mo_ti_a = f['MO_ALPHA_TYPEINDICES'][:]
else:
mo_ti_a = [b'?' for i in mo_oc_a]
else:
mo_en_a = []
mo_oc_a = []
mo_cf_a = []
mo_ti_a = []
if ('MO_BETA_ENERGIES' in f):
mo_en_b = f['MO_BETA_ENERGIES'][:]
mo_oc_b = f['MO_BETA_OCCUPATIONS'][:]
mo_cf_b = f['MO_BETA_VECTORS'][:]
if ('MO_BETA_TYPEINDICES' in f):
mo_ti_b = f['MO_BETA_TYPEINDICES'][:]
else:
mo_ti_b = [b'?' for i in mo_oc_b]
else:
mo_en_b = []
mo_oc_b = []
mo_cf_b = []
mo_ti_b = []
mo_ti = [str(i.decode('ascii')) for i in mo_ti]
mo_ti_a = [str(i.decode('ascii')) for i in mo_ti_a]
mo_ti_b = [str(i.decode('ascii')) for i in mo_ti_b]
self.base_MO = {}
self.base_MO[0] = [{'ene':e, 'occup':o, 'type':t} for e,o,t in zip(mo_en, mo_oc, mo_ti)]
self.base_MO['a'] = [{'ene':e, 'occup':o, 'type':t} for e,o,t in zip(mo_en_a, mo_oc_a, mo_ti_a)]
self.base_MO['b'] = [{'ene':e, 'occup':o, 'type':t} for e,o,t in zip(mo_en_b, mo_oc_b, mo_ti_b)]
# Read the coefficients
ii = [sum(self.N_bas[:i]) for i in range(len(self.N_bas))]
j = 0
for i,b,s in zip(ii, self.N_bas, self.irrep):
for orb,orb_a,orb_b in zip_longest(self.base_MO[0][i:i+b], self.base_MO['a'][i:i+b], self.base_MO['b'][i:i+b]):
if (orb):
orb['sym'] = s
orb['coeff'] = np.zeros(sum(self.N_bas))
orb['coeff'][i:i+b] = mo_cf[j:j+b]
if (orb_a):
orb_a['sym'] = s
orb_a['coeff'] = np.zeros(sum(self.N_bas))
orb_a['coeff'][i:i+b] = mo_cf_a[j:j+b]
if (orb_b):
orb_b['sym'] = s
orb_b['coeff'] = np.zeros(sum(self.N_bas))
orb_b['coeff'][i:i+b] = mo_cf_b[j:j+b]
j += b
# Desymmetrize the MOs
if (len(self.N_bas) > 1):
for orbs in self.base_MO.values():
for orb in orbs:
orb['coeff'] = np.dot(self.mat, orb['coeff'])
self.roots = [(0, 'Average')]
self.sdm = None
self.tdm = None
self.tsdm = None
mod = None
if ('MOLCAS_MODULE' in f.attrs):
mod = f.attrs['MOLCAS_MODULE'].decode('ascii')
if (mod == 'CASPT2'):
self.wf = 'PT2'
self.roots[0] = (0, 'Reference')
if ('DENSITY_MATRIX' in f):
rootids = f.attrs['STATE_ROOTID'][:]
# For MS-CASPT2, the densities are SS, but the energies are MS,
# so take the energies from the effective Hamiltonian matrix instead
if ('H_EFF' in f):
H_eff = f['H_EFF']
self.roots.extend([(i+1, '{0}: {1:.6f}'.format(r, e)) for i,(r,e) in enumerate(zip(rootids, np.diag(H_eff)))])
self.msroots = [(0, 'Reference')]
self.msroots.extend([(i+1, '{0}: {1:.6f}'.format(i+1, e)) for i,e in enumerate(f['STATE_PT2_ENERGIES'])])
else:
self.roots.extend([(i+1, '{0}: {1:.6f}'.format(r, e)) for i,(r,e) in enumerate(zip(rootids, f['STATE_PT2_ENERGIES']))])
elif ('SFS_TRANSITION_DENSITIES' in f):
# This is a RASSI-like calculation, with TDMs
self.wf = 'SI'
self.roots = [(0, 'Average')]
mult = f.attrs['STATE_SPINMULT']
sym = f.attrs['STATE_IRREPS']
ene = f['SFS_ENERGIES'][:]
sup = [u'⁰', u'¹', u'²', u'³', u'⁴', u'⁵', u'⁶', u'⁷', u'⁸', u'⁹']
rootlist = []
for i,(e,m,s) in enumerate(zip(ene, mult, sym)):
try:
l = list(self.irrep[s-1])
except IndexError:
l = ['?']
l[0] = l[0].upper()
l = ''.join([sup[int(d)] for d in str(m)] + l)
rootlist.append((i+1, u'{0}: ({1}) {2:.6f}'.format(i+1, l, e)))
self.ene_idx = np.argsort(ene)
for i in self.ene_idx:
self.roots.append(rootlist[i])
self.tdm = True
if ('SFS_TRANSITION_SPIN_DENSITIES' in f):
self.sdm = True
if (any(mult != 1)):
self.tsdm = True
# find out which transitions are actually nonzero (stored)
self.have_tdm = np.zeros((len(ene), len(ene)), dtype=bool)
tdm = f['SFS_TRANSITION_DENSITIES']
if ('SFS_TRANSITION_SPIN_DENSITIES' in f):
tsdm = f['SFS_TRANSITION_SPIN_DENSITIES']
else:
tsdm = np.zeros_like(self.have_tdm)
for j in range(len(ene)):
for i in range(j):
if ((not np.allclose(tdm[i,j], 0)) or (not np.allclose(tsdm[i,j], 0))):
self.have_tdm[i,j] = True
self.have_tdm[j,i] = True
if (not np.any(self.have_tdm)):
self.tdm = False
else:
if ('DENSITY_MATRIX' in f):
rootids = [i+1 for i in range(f.attrs['NROOTS'])]
self.roots.extend([(i+1, '{0}: {1:.6f}'.format(r, e)) for i,(r,e) in enumerate(zip(rootids, f['ROOT_ENERGIES']))])
if ('SPINDENSITY_MATRIX' in f):
sdm = f['SPINDENSITY_MATRIX']
if (not np.allclose(sdm, 0)):
self.sdm = True
if ('TRANSITION_DENSITY_MATRIX' in f):
tdm = f['TRANSITION_DENSITY_MATRIX']
if (not np.allclose(tdm, 0)):
self.tdm = True
if ('TRANSITION_SPIN_DENSITY_MATRIX' in f):
tsdm = f['TRANSITION_SPIN_DENSITY_MATRIX']
if (not np.allclose(tsdm, 0)):
self.tsdm = True
self.tdm = True
# find out which transitions are actually nonzero (stored)
n = int(np.sqrt(1+8*tdm.shape[0])+1)//2
self.have_tdm = np.zeros((n, n), dtype=bool)
if ('TRANSITION_SPIN_DENSITY_MATRIX' not in f):
tsdm = np.zeros(n*(n-1)//2)
for i in range(n):
for j in range(i):
m = j*(j-1)//2+i
if ((not np.allclose(tdm[m], 0)) or (not np.allclose(tsdm[m], 0))):
self.have_tdm[i,j] = True
self.have_tdm[j,i] = True
if (not np.any(self.have_tdm)):
self.tdm = False
if ('WFA' in f):
self.wfa_orbs = []
for i in f['WFA'].keys():
match = re.match('DESYM_(.*)_VECTORS', i)
if (match):
self.wfa_orbs.append(match.group(1))
# Read the optional notes
if ('Pegamoid_notes' in f):
self.notes = [str(i.decode('ascii')) for i in f['Pegamoid_notes'][:]]
# Obtain a different type of MO orbitals (only for HDF5 files)
def get_orbitals(self, density, root):
if (self.file != self.h5file):
return
if ((density, root) == self.current_orbs):
return
# in a PT2 calculation, density matrices include all orbitals (not only active)
fix_frozen = False
if (self.wf == 'PT2'):
if (density == 'State'):
fix_frozen = True
tp_act = set([o['type'] for o in self.base_MO[0]])
tp_act -= set(['F', 'D'])
elif (self.wf == 'SI'):
tp_act = ['?']
else:
tp_act = ['1', '2', '3']
sym = 0
transpose = False
sdm = False
# Depending on the type of density selected, read the matrix
# and choose the appropriate method to get the orbitals
if (density == 'State'):
if (root == 0):
if (self.wf == 'SI'):
algo = 'eigsort'
n = len(self.roots)-1
with h5py.File(self.file, 'r') as f:
dm = f['SFS_TRANSITION_DENSITIES'][0,0,:]
if (self.sdm):
sdm = f['SFS_TRANSITION_SPIN_DENSITIES'][0,0,:]
for i in range(1, n):
dm += f['SFS_TRANSITION_DENSITIES'][i,i,:]
if (self.sdm):
sdm += f['SFS_TRANSITION_SPIN_DENSITIES'][i,i,:]
dm /= n
if (self.sdm):
sdm /= n
else:
if (self.sdm):
algo = 'eig'
dm = np.diag([o['occup'] for o in self.base_MO[0] if (o['type'] in tp_act)])
with h5py.File(self.file, 'r') as f:
sdm = np.mean(f['SPINDENSITY_MATRIX'], axis=0)
else:
algo = 'non'
self.MO = self.base_MO[0]
self.MO_a = self.base_MO['a']
self.MO_b = self.base_MO['b']
else:
algo = 'eig'
with h5py.File(self.file, 'r') as f:
if (self.wf == 'SI'):
algo += 'sort'
dm = f['SFS_TRANSITION_DENSITIES'][root-1,root-1,:]
if (self.sdm):
sdm = f['SFS_TRANSITION_SPIN_DENSITIES'][root-1,root-1,:]
else:
dm = f['DENSITY_MATRIX'][root-1,:]
if (self.sdm):
sdm = f['SPINDENSITY_MATRIX'][root-1,:]
elif (density == 'Spin'):
algo = 'eig'
if (root == 0):
if (self.wf == 'SI'):
algo += 'sort'
n = len(self.roots)-1
with h5py.File(self.file, 'r') as f:
dm = f['SFS_TRANSITION_SPIN_DENSITIES'][0,0,:]
for i in range(1, n):
dm += f['SFS_TRANSITION_SPIN_DENSITIES'][i,i,:]
dm /= n
else:
with h5py.File(self.file, 'r') as f:
dm = np.mean(f['SPINDENSITY_MATRIX'], axis=0)
else:
with h5py.File(self.file, 'r') as f:
if (self.wf == 'SI'):
algo += 'sort'
dm = f['SFS_TRANSITION_SPIN_DENSITIES'][root-1,root-1,:]
else:
dm = f['SPINDENSITY_MATRIX'][root-1,:]
elif (density.split()[0] in ['Difference', 'Transition']):
r1 = root[0] - 1
r2 = root[1] - 1
if (density == 'Difference'):
algo = 'eig'
with h5py.File(self.file, 'r') as f:
if (self.wf == 'SI'):
algo += 'sort'
dm = f['SFS_TRANSITION_DENSITIES'][r2,r2,:] - f['SFS_TRANSITION_DENSITIES'][r1,r1,:]
if (self.sdm):
sdm = f['SFS_TRANSITION_SPIN_DENSITIES'][r2,r2,:] - f['SFS_TRANSITION_SPIN_DENSITIES'][r1,r1,:]
else:
dm = f['DENSITY_MATRIX'][r2,:] - f['DENSITY_MATRIX'][r1,:]
if (self.sdm):
sdm = f['SPINDENSITY_MATRIX'][r2,:] - f['SPINDENSITY_MATRIX'][r1,:]
elif (density.startswith('Transition')):
if (r1 > r2):
r1, r2 = (r2, r1)
transpose = True
algo = 'svd'
fact = 0
if ('(alpha)' in density):
fact = 1
elif ('(beta)' in density):
fact = -1
with h5py.File(self.file, 'r') as f:
if (self.wf == 'SI'):
# find the symmetry of the transition
sym = f.attrs['STATE_IRREPS']
sym = symmult[sym[r1]-1, sym[r2]-1]
dm = f['SFS_TRANSITION_DENSITIES'][r1,r2,:]
if (fact != 0):
dm = (dm + fact * f['SFS_TRANSITION_SPIN_DENSITIES'][r1,r2,:]) / np.sqrt(2)
else:
n = r2*(r2-1)//2+r1
dm = f['TRANSITION_DENSITY_MATRIX'][n,:]
if (fact != 0):
dm = (dm + fact * f['TRANSITION_SPIN_DENSITY_MATRIX'][n,:]) / np.sqrt(2)
elif (density == 'WFA'):
algo = 'non'
label = self.wfa_orbs[root]
new_MO = []
with h5py.File(self.h5file, 'r') as f:
occ = f['WFA/DESYM_{0}_OCCUPATIONS'.format(label)][:]
norb = len(occ)
vec = np.reshape(f['WFA/DESYM_{0}_VECTORS'.format(label)], (norb, -1))
for i,o in enumerate(occ):
new_MO.append({'coeff': vec[i,:], 'occup': o, 'type': '?', 'ene': 0.0, 'sym': 'z'})
# if these are NTOs, split in hole and particle
ntos = ('_NTO' in label) and (norb%2 == 0)
if (ntos):
for i in range(norb//2):
if (not np.isclose(occ[i], -occ[-i-1], atol=self.eps)):
ntos = False
break
if (ntos):
self.MO = []
self.MO_a = sorted(new_MO[int(norb/2):], key=lambda i: i['occup'], reverse=True)
self.MO_b = sorted(new_MO[:int(norb/2)], key=lambda i: i['occup'])
for o in self.MO_b:
o['occup'] *= -1
else:
self.MO = new_MO
self.MO_a = []
self.MO_b = []
if (sdm is not False):
if (np.allclose(sdm, 0)):
sdm = False
# In RASSI, DMs are stored in (symmetrized) AO basis
if ((self.wf == 'SI') and ('non' not in algo)):
with h5py.File(self.file, 'r') as f:
S = f['AO_OVERLAP_MATRIX'][:]
tot = sum(self.N_bas)
full_S = np.zeros((tot, tot))
full_D = np.zeros((tot, tot))
if (sdm is not False):
full_sD = np.zeros((tot, tot))
i = 0
k = 0
for s1,n1 in enumerate(self.N_bas):
j1 = np.sum(self.N_bas[:s1])
full_S[j1:j1+n1,j1:j1+n1] = np.reshape(S[i:i+n1*n1], (n1, n1))
i += n1*n1
s2 = np.flatnonzero(symmult[s1,:] == sym)[0]
n2 = self.N_bas[s2]
j2 = np.sum(self.N_bas[:s2])
full_D[j2:j2+n2,j1:j1+n1] = np.reshape(dm[k:k+n1*n2], (n2, n1))
if (sdm is not False):
full_sD[j2:j2+n2,j1:j1+n1] = np.reshape(sdm[k:k+n1*n2], (n2, n1))
k += n1*n2
S = full_S
dm = full_D
if (sdm is not False):
sdm = full_sD
s, V = np.linalg.eigh(S)
sqrtS = np.dot(V * np.sqrt(s), V.T)
# convert to Löwdin-orthogonalized symmetrized AOs
dm = np.dot(sqrtS, np.dot(dm, sqrtS))
if (sdm is not False):
sdm = np.dot(sqrtS, np.dot(sdm, sqrtS))
X = np.dot(V / np.sqrt(s), V.T)
AO = True
else:
AO = False
# Now obtain the MOs from the density matrix
mix_threshold = 1.0-1e-4
# as eigenvectors
if ('eig' in algo):
dmlist = [dm, None, None]
MOlist = ['MO', 'MO_a', 'MO_b']
if (sdm is not False):
dmlist[1] = 0.5 * (dm + sdm) # alpha density
dmlist[2] = 0.5 * (dm - sdm) # beta density
for dm_,MO_ in zip(dmlist, MOlist):
if (dm_ is None):
setattr(self, MO_, [])
continue
# In AO (RASSI), there is no "base" set of MOs, use Löwdin symmetrized AOs
if (AO):
new_MO = []
for s,nbas in enumerate(self.N_bas):
new_MO.extend([{'sym': self.irrep[s], 'type': '?', 'ene': 0.0} for i in range(nbas)])
for i,o in enumerate(new_MO):
o['coeff'] = X[:,i]
act = new_MO
else:
new_MO = deepcopy(self.base_MO[0])
act = [o for o in new_MO if (o['type'] in tp_act)]
if (density in ['Difference', 'Spin']):
for o in new_MO:
o['hide'] = True
o['occup'] = 0.0
o['ene'] = 0.0
# If density matrix is one-dimensional (e.g. from CASPT2), it is stored by symmetry blocks,
# reconstruct the full matrix
if (len(dm_.shape) == 1):
full_dm = np.zeros((len(act), len(act)))
nMO = [(sum(self.N_bas[:i]), sum(self.N_bas[:i+1])) for i in range(len(self.N_bas))]
j = 0
k = 0
for i,nbas in zip(nMO, self.N_bas):
n = len([o for o in new_MO[i[0]:i[1]] if (o['type'] in tp_act)])
j1 = int(n*(n+1)/2)
sym_dm = np.zeros((n, n))
sym_dm[np.tril_indices(n, 0)] = dm_[j:j+j1]
sym_dm = sym_dm + np.tril(sym_dm, -1).T
full_dm[k:k+n,k:k+n] = sym_dm
j += j1
k += n
dm_ = full_dm
for s in set([o['sym'] for o in act]):
symidx = [i for i,o in enumerate(act) if (o['sym'] == s)]
symdm = dm_[np.ix_(symidx,symidx)]
symact = [act[i] for i in symidx]
occ, vec = np.linalg.eigh(symdm)
# Match similar orbitals
if ('sort' in algo):
idx = occ.argsort()[::-1]
else:
freevec = [1 for i in symidx]
idx = [-1 for i in symidx]
for i in range(len(symidx)):
idx[i] = np.argmax(abs(vec[i,:]*freevec))
if (vec[i,idx[i]] < 0.0):
vec[:,idx[i]] *= -1
freevec[idx[i]] = 0
occ = occ[idx]
vec = vec[:,idx]
# Check contributions for each orbital, if there are significant
# contributions from several types, mark this type as '?'
mix = []
for i in range(vec.shape[1]):
types = []
for a in tp_act:
if (sum([j**2 for n,j in enumerate(vec[:,i]) if (symact[n]['type']==a)]) > mix_threshold):
types.append(a)
mix.append(types[0] if (len(types) == 1) else '?')
new_coeff = np.dot(vec.T, [o['coeff'] for o in symact])
# If in AO, desymmetrize the resulting MOs
if (AO):
if (len(self.N_bas) > 1):
new_coeff = np.dot(new_coeff, self.mat.T)
occ[np.abs(occ) < 10*self.eps] = 0.0
for o,n,c,t in zip(symact, occ, new_coeff, mix):
o['hide'] = False
o['occup'] = n
o['coeff'] = c
o['type'] = t
o['ene'] = 0.0
if (AO and (density in ['Difference', 'Spin'])):
for o in symact:
#if (np.abs(o['occup']) < self.eps):
if (np.abs(o['occup']) < 1e-6):
o['hide'] = True
if (fix_frozen):
for o in new_MO:
if (o['type'] == 'F'):
o['occup'] = 2.0
setattr(self, MO_, new_MO)
# as singular vectors
elif ('svd' in algo):
if (AO):
new_MO_l = []
new_MO_r = []
for s,nbas in enumerate(self.N_bas):
new_MO_l.extend([{'sym': self.irrep[s], 'type': '?', 'ene': 0.0} for i in range(nbas)])
new_MO_r.extend([{'sym': self.irrep[s], 'type': '?', 'ene': 0.0} for i in range(nbas)])
for i,(o_l,o_r) in enumerate(zip(new_MO_l, new_MO_r)):
o_l['coeff'] = X[:,i]
o_r['coeff'] = X[:,i]
act_l = new_MO_l
act_r = new_MO_r
else:
new_MO_l = deepcopy(self.base_MO[0])
new_MO_r = deepcopy(self.base_MO[0])
act_l = [o for o in new_MO_l if (o['type'] in tp_act)]
act_r = [o for o in new_MO_r if (o['type'] in tp_act)]
for o in new_MO_l + new_MO_r:
o['hide'] = True
o['occup'] = 0.0
o['ene'] = 0.0
rl_sort = np.array([-1 for i in act_l])
for s in set([o['sym'] for o in act_l]):
s1 = self.irrep.index(s)
s2 = np.flatnonzero(symmult[s1,:] == sym)[0]
symidx1 = [i for i,o in enumerate(act_l) if (o['sym'] == s)]
symidx2 = [i for i,o in enumerate(act_r) if (o['sym'] == self.irrep[s2])]
symdm = dm[np.ix_(symidx1,symidx2)]
symact_l = [act_l[i] for i in symidx1]
symact_r = [act_r[i] for i in symidx2]
if (not np.allclose(symdm, 0)):
vec_l, occ, vec_r = np.linalg.svd(symdm)
vec_r = vec_r.T
mix_l = []
mix_r = []
for i in range(vec_l.shape[1]):
types = []
for a in tp_act:
if (sum([j**2 for n,j in enumerate(vec_l[:,i]) if (symact_l[n]['type']==a)]) > mix_threshold):
types.append(a)
mix_l.append(types[0] if (len(types) == 1) else '?')
for i in range(vec_r.shape[1]):
types = []
for a in tp_act:
if (sum([j**2 for n,j in enumerate(vec_r[:,i]) if (symact_r[n]['type']==a)]) > mix_threshold):
types.append(a)
mix_r.append(types[0] if (len(types) == 1) else '?')
new_coeff_l = np.dot(vec_l.T, [o['coeff'] for o in symact_l])
new_coeff_r = np.dot(vec_r.T, [o['coeff'] for o in symact_r])
for i,j in zip(symidx1, symidx2):
rl_sort[j] = i
else:
continue
if (AO):
if (len(self.N_bas) > 1):
new_coeff_l = np.dot(new_coeff_l, self.mat.T)
new_coeff_r = np.dot(new_coeff_r, self.mat.T)
occ[np.abs(occ) < 10*self.eps] = 0.0
for o_l,o_r,n,cl,cr,tl,tr in zip(symact_l, symact_r, occ, new_coeff_l, new_coeff_r, mix_l, mix_r):
o_l['hide'] = False
o_l['occup'] = n**2/2
o_l['coeff'] = cl
o_l['type'] = tl
o_r['hide'] = False
o_r['occup'] = n**2/2
o_r['coeff'] = cr
o_r['type'] = tr
if (AO):
for o in symact_l+symact_r:
#if (np.abs(o['occup']) < self.eps):
if (np.abs(o['occup']) < 1e-6):
o['hide'] = True
# reorder the right NTOs to match them with the left, since they may come from different symmetries
resort = deepcopy(new_MO_r)
for i in np.flatnonzero(rl_sort >=0):
resort[i] = new_MO_r[rl_sort[i]]
new_MO_r = resort
self.MO = []
if (transpose):
self.MO_a = new_MO_l
self.MO_b = new_MO_r
else:
self.MO_a = new_MO_r
self.MO_b = new_MO_l
self.current_orbs = (density, root)
# Read basis set from a Molden file
def read_molden_basis(self):
with open(self.file, 'r') as f:
# Molden supports up to g functions, and by default all are Cartesian
ang_labels = 'spdfg'
cart = [False, True, True, True, True]
# Specify the order of the Cartesian components according to the convention (see ang)
order = []
order.append([0])
order.append([-1, 0, 1])
order.append([-2, 1, 3, -1, 0, 2])
order.append([-3, 3, 6, 0, -2, -1, 2, 5, 4, 1])
order.append([-4, 6, 10, -3, -2, 2, 7, 5, 9, -1, 1, 8, 0, 3, 4])
maxl = 0
done = True
if (re.search(r'\[MOLDEN FORMAT\]', f.readline(), re.IGNORECASE)):
done = False
num = None
line = ' '
while ((not done) and (line != '')):
line = f.readline()
# Read the geometry
if re.search(r'\[N_ATOMS\]', line, re.IGNORECASE):
num = int(f.readline())
elif re.search(r'\[ATOMS\]', line, re.IGNORECASE):
unit = 1
if (re.search(r'Angs', line, re.IGNORECASE)):
unit = angstrom
self.centers = []
if (num is None):
num = 0
while True:
save = f.tell()
try:
l, _, q, x, y, z = f.readline().split()
q = min(maxZ, max(0, int(q)))
self.centers.append({'name':l, 'Z':q, 'xyz':np.array([fortran_float(x), fortran_float(y), fortran_float(z)])*unit})
num += 1
except:
f.seek(save)
break
else:
for i in range(num):
l, _, q, x, y, z = f.readline().split()
q = min(maxZ, max(0, int(q)))
self.centers.append({'name':l, 'Z':q, 'xyz':np.array([fortran_float(x), fortran_float(y), fortran_float(z)])*unit})
self.geomcenter = (np.amin([c['xyz'] for c in self.centers], axis=0) + np.amax([c['xyz'] for c in self.centers], axis=0))/2
# Read tags for spherical shells
elif re.search(r'\[5D\]', line, re.IGNORECASE):
cart[2] = False
cart[3] = False
elif re.search(r'\[5D7F\]', line, re.IGNORECASE):
cart[2] = False
cart[3] = False
elif re.search(r'\[5D10F\]', line, re.IGNORECASE):
cart[2] = False
cart[3] = True
elif re.search(r'\[7F\]', line, re.IGNORECASE):
cart[3] = False
elif re.search(r'\[9G\]', line, re.IGNORECASE):
cart[4] = False
# Make sure we read the basis after the Cartesian types are known
# (we still assume [MO] will be after all this)
elif re.search(r'\[GTO\]', line, re.IGNORECASE):
save_GTO = f.tell()
# Read basis functions: a series of blank-separated blocks
# starting with center number and followed by all the shells,
# each is the angular momentum letter and number of primitives,
# plus this number of exponents and coefficients.
elif re.search(r'\[MO\]', line, re.IGNORECASE):
save_MO = f.tell()
f.seek(save_GTO)
bf_id = []
while (True):
save = f.tell()
# First find out if this is another center, or the basis set
# specification has finished
try:
n = int(f.readline().split()[0])
except:
f.seek(save)
break
try:
self.centers[n-1]['cart'] = []
except IndexError:
error = 'Basis functions found for unspecified centers. Broken molden file?'
raise Exception(error)
basis = {}
# Read the shells for this center
while (True):
try: