-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsd_core.f90
1624 lines (1623 loc) · 67.5 KB
/
sd_core.f90
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
!
! 14 March 2009, SP - Initial version, from fission of superdyson.f90
! 09 September 2009, SP - Try to exploit frozen cores and identical MO spaces
! 21 May 2010, SP - Incorporate new GAMESS import routines, incl. 1-e
! integrals. This should allow calculation of completely
! general matrix elements, with no restriction on basis
! sets or molecular geometries.
! 12 July 2010, SP - Add an option not to calculate transition dipoles
! with the new-style input
! 18 June 2011, SP - Add option to disable S^2 and normalization checks
! 05 July 2013, SP - Added support for even newer-style namelist input
! and access to more one-electron matrix elements
!
! Core routines for dealing with 1-electron operators between
! general CI wavefunctions.
!
! WARNING: The format of the input determinant list is perhaps misleading.
! The signs of the amplitudes for the determinants are expected
! for the alpha/beta ordering of the molecular orbitals (ie, the
! determinants should have all orbitals associated with spin alpha
! at the beginning, followed by all spin-beta orbitals). This is
! the sign choice used by GAMESS CI routines and internally by my
! routines here. The commonly used "textbook" sign choice (with
! orbitals appearing in the fixed canonical order, and alpha/beta
! spins intemixed) can be obtained by multiplying amplitudes by
! (-1)**parity, where parity is the number of permutations leading
! from the alpha/beta to the canonical ordering.
!
module sd_core
use accuracy
use timer
use import_gamess
use lapack
use math
use block_diag
use block_matrices
use block_determinant
use sort_tools
implicit none
public ! The sd_core module is not supposed to lead an
! independent existence - it is always a part of
! a specific integral package. Hence, everything
! must be public.
private sd_v2
!
character(len=clen), save :: rcsid_sd_core = "$Id: sd_core.f90,v 1.31 2022/08/02 14:12:59 ps Exp $"
!
integer, parameter :: lk = kind(.true.)
integer, parameter :: lchar = 1024
!
! Local data for dyson
!
integer(ik) :: verbose = 1 ! Debug level
integer(ik), parameter :: ioscr = 10 ! Local use in subroutines
real(rk), parameter :: eps_s = 1e-4_rk ! Tolerance for self-overlap matrix elements
real(rk) :: eps_cdet = 1e-5_rk ! Cut-off for determinant weights product
real(rk) :: eps_integral = 1e-10_rk
! Tolerance for blocking overlap and property matrices
!
integer(ik) :: naos_bra = -1 ! Total number of AOs in the "reference"
integer(ik) :: naos_ket = -1 ! Total number of AOs in the "ion"
! -1 means that the value will be obtained from
! the MOs input block
integer(ik) :: nmobra = 0 ! Total number of orbitals for the
integer(ik) :: nmoket = 0 ! reference and the ion
integer(ik) :: ndetbra = 0 ! Total number of determinants for the
integer(ik) :: ndetket = 0 ! reference and the ion
real(rk) :: sign_bra = 1.0_rk ! Overall scale factor to apply to the bra wavefunction
real(rk) :: sign_ket = 1.0_rk ! Overall scale factor to apply to the ket wavefunction
!
integer(ik) :: nelbra ! Number of electrons in the reference
integer(ik) :: nelket ! Number of electrons in the ion
logical(lk) :: have_overlaps ! Overlap integrals are available; MO sets may differ
logical(lk) :: have_operator = .true.
! One-electron integrals are available, calculate
! corresponding matrix elements
logical(lk) :: mode_overlap ! Same number of electrons in both wfs
logical(lk) :: same_basis ! Is the basis set identical for the ref and ion?
logical(lk) :: same_mos ! MOs used to reprent the two WFs are identical
integer(ik) :: frozen_orbitals ! Number of identical doubly-occupied orbitals in
! all determinants.
logical(lk) :: dont_check_s2 = .false.
! True if normalization and S^2 checks are disabled
character(len=20) :: task = 'braket' ! Task to perform; currently one of:
! 'braket' or '1-rdm'
character(len=20) :: braket_operator &
= 'dipole' ! Operator to apply for task='bracket'. Can be one of:
! 'none' -
! 'kinetic' - kinetic energy (-0.5 \Nabla)
! 'dipole' - coordinate expectation
! 'velocity' - velocity \Delta (ie momentum expectation
! sans -I*hbar)
! 'r-d/dr' - r_i (d/d r_j) operator. The i index runs fast
! Operators below require coordinates of the special centre:
! '1/r' - Coulomb nuclear attraction
! 'r' - Linear attraction
! 'r**2' - Harmonic attraction
! Operators below require both coordinates of the centre, and
! Gaussian exponent:
! 'gauss' - Gaussian attraction
! 'gauss r' - Gaussian-damped linear attraction
! See import_gamess.f90 for the complete list of
! operators potentially available in the integrals
! package.
! Regardless of the choice of operator, overlap
! matrix element will be evaluated as well.
real(rk) :: op_centre(3) = 0. ! Coordinates of the special centre (if any)
real(rk) :: op_gauss = 1. ! Gaussian exponent (if any)
integer(ik) :: op_components = 3 ! Number of components in the operator; must
! match the choice of operator in braket_operator above
logical :: use_symmetry = .true.
! Set to true to try to use structure of the overlap
! matrix for eliminating unnecessary work. There is
! a slight overhead associated with this, so there
! may be cases where symmetry use is not beneficial.
logical :: detail_timers = .false.
! Record timing of low-level routines, potentially with timing
! overhead exceeding the actual work done
!
! Composition of the determinants is encoded by:
! 2 = alpha+beta
! +1 = alpha
! -1 = beta
! 0 = unoccupied
!
integer(ik), allocatable :: occbra(:,:) ! Occupation numbers (0/1) for reference and
integer(ik), allocatable :: occket(:,:) ! ion. The last index is determinant.
real(rk), allocatable :: cdetref(:) ! Weights of the determinants for the reference
real(rk), allocatable :: cdetion(:) ! and the ion
!
! For some of our routines (analyze_spinorbit_integral_blocks, which is in a critical path)
! it is advantageous to record sorted lists of orbitals separately. We also need to record
! the sorting sequence, so that we can flip back and forth.
!
! The spinbra, spinket, ord_bra, and ord_ket arrays are only needed if symmetry is on.
!
integer(ik), allocatable :: scnt_bra (:,:) ! Counts of spin-orbitals for the bra determinants
! First index: 1=alpha, 2=beta
! Second index: the determinant
integer(ik), allocatable :: scnt_ket (:,:) ! ditto for the ket
integer(ik), allocatable :: spinbra (:,:,:) ! First index: up to nelbra spin-orbitals, sorted, grouped by spin
! Second index: 1=alpha, 2=beta
! Third index: determinant
integer(ik), allocatable :: spinket (:,:,:) ! ditto, for the ket
integer(ik), allocatable :: ord_bra (:,:,:) ! Original spin-MO index for each sorted spin-MO in spinbra
integer(ik), allocatable :: ord_ket (:,:,:) ! ditto, ket
!
real(rk), allocatable :: cbra (:,:) ! Reference spatial orbitals
real(rk), allocatable :: cket (:,:) ! Ion spatial orbitals
real(rk), allocatable :: sao_ref(:,:) ! Reference overlap matrix in AO basis
real(rk), allocatable :: sao_ion(:,:) ! Ion overlap matrix in AO basis
real(rk), allocatable :: sao (:,:) ! Cross-system overlap matrix in AO basis
real(rk), allocatable :: dipao(:,:,:) ! Dipole matrix in AO basis, X/Y/Z along last ind.
!
character(len=lchar) :: comment = ' ' ! Comment line
character(len=14), allocatable &
:: aolabels(:) ! Atomic orbital labels, for human readability
!
! File names
!
character(len=lchar) :: file_sao = ' ' ! File name for the overlaps and (optionally) dipoles
character(len=lchar) :: file_cbra = ' ' ! File name for the reference MOs
character(len=lchar) :: file_cket = ' ' ! File name for the ion MOs
character(len=lchar) :: file_detbra = ' ' ! File name for the reference CI coefficients
character(len=lchar) :: file_detket = ' ' ! File name for the ion CI coeffficients
!
! The parent MO index is on the left; the ion MO index is on the right.
!
real(rk), allocatable :: sx (:,:) ! 1-particle overlaps between MOs of reference and target
real(rk), allocatable :: dipx (:,:,:) ! Dipole matrix elements between the MOs. Last index is X/Y/Z
! Could be some other operator, too (see braket_operator),
! but since the logic does not change, who cares?
type(bm_structure) :: sx_ms ! Block structure of the overlap matrix
type(bm_structure) :: dipx_ms ! Block structure of the property matrices
!
! Stuff for the new-style input:
!
logical(lk) :: own_integrals = .true.
! We'll calculate our own integrals using GAMESS ckpt.
type(gam_structure) :: g_bra ! Structure and basis description from GAMESS checkpoint
type(gam_structure) :: g_ket ! ditto, the ion
!
namelist /sd_v2/ verbose, &
task, braket_operator, op_centre, op_gauss, &
comment, &
nmobra, nmoket, ndetbra, ndetket, eps_cdet, &
dont_check_s2, &
file_cbra, file_cket, file_detbra, file_detket, &
eps_integral, use_symmetry, detail_timers, &
sign_bra, sign_ket
contains
subroutine read_input_new
integer(ik) :: info
!
call math_init
call TimerStart('read_input_new')
!
read (input,nml=sd_v2,iostat=info)
if (info/=0) then
write (out,"(/'**** ENCOUNTERED ERROR ',i4,' READING INPUT NAMELIST ****'/)") info
end if
write (out,"(/'==== Options ====')")
write (out,nml=sd_v2)
write (out,"()")
!
if (nmobra <=0 .or. nmoket<=0 .or. ndetbra<=0 .or. ndetket<=0 ) then
write (out,"('ERROR: All of nmobra, nmoket, ndetbra, and ndetket must be given.')")
stop 'sd_core%read_input_new - missing values with no defaults'
end if
!
call TimerStop('read_input_new')
end subroutine read_input_new
!
subroutine read_core_input
call math_init
call TimerStart('read_core_input')
!
! LINE 0: comment, to be added to the output
!
read(input,"(a)") comment
!
! LINE 1: naos nmobra nmoket ndetbra ndetket cutoff
! If naos is set to -1, this will trigger the new-style
! input format. naos of -2 will be internally converted
! to -1 together with have_operator = .false.
!
read(input,*) naos_bra, nmobra, nmoket, ndetbra, ndetket, eps_cdet
!
naos_ket = naos_bra
own_integrals = .true.
dont_check_s2 = .true.
select case (naos_bra)
case ( -1) ; have_overlaps = .true. ; have_operator = .true. ; dont_check_s2 = .false.
case ( -2) ; have_overlaps = .true. ; have_operator = .false. ; dont_check_s2 = .false.
case (-11) ; have_overlaps = .true. ; have_operator = .true. ;
case (-12) ; have_overlaps = .true. ; have_operator = .false. ;
case default
own_integrals = .false.
end select
!
! LINE 2: name of the file (in GAMESS NPRINT=3 format) containing
! the overlap integrals. Not used in the new format.
!
if (.not.own_integrals) then
read(input,*) file_sao
end if
!
! LINE 3: name of the file (in GAMESS PUNCH format) containing reference MOs
! In the old format, everything up to and including the $VEC line
! and everyting starting at the $END line should be removed.
! In the new format, the file should contain a C1-symmetric geometry
! with an inline basis specification (see import_gamess.f90).
!
read(input,*) file_cbra
!
! LINE 4: name of the file (in GAMESS PUNCH format) containing ion MOs
!
read(input,*) file_cket
!
! LINE 5: name of the file containing determinants for the parent
! Composition of the determinants for the parent wavefunction.
! Each determinant is defined by:
! weight (real)
! orbital occupations (nmobra values; 2, 0, +1, or -1)
!
read(input,*) file_detbra
!
! LINE 6: Name of the file, containing ion wavefunction. This is not
! produced directly by GAMESS; See det_conv.awk for an example
! of how to post-process GAMESS output for some interesting
! cases.
!
read(input,*) file_detket
call TimerStop('read_core_input')
end subroutine read_core_input
!
subroutine math_init
real(rk) :: dummy
!
if (log10(huge(1._rk))>=4930._rk) then
! Quad precision; don't forget "factorial_slack"
dummy = MathFactorial(1750_ik-5_ik)
dummy = MathLogFactorial(20000_ik)
else if (log10(huge(1._rk))>=308._rk) then
! Double precision
dummy = MathFactorial(170_ik-5_ik)
dummy = MathLogFactorial(10000_ik)
else
! Single precision
dummy = MathFactorial(34_ik-5_ik)
dummy = MathLogFactorial(500_ik)
end if
dummy = MathDoubleFactorial(80)
end subroutine math_init
!
subroutine initialize_core_data
integer(ik) :: detref, detion, status
integer(ik) :: ic, jc
!
call TimerStart('initialize_core_data')
!
! Echo core input parameters
!
if (verbose>=0) then
write(out,"(a)") trim(comment)
write (out,"(' task = ',a)") trim(task)
write (out,"(' operator = ',a)") trim(braket_operator)
write (out,"(' naos = ',i6)") naos_bra
write (out,"(' nmobra = ',i6)") nmobra
write (out,"(' nmoket = ',i6)") nmoket
write (out,"(' ndetbra = ',i6)") ndetbra
write (out,"(' ndetket = ',i6)") ndetket
write (out,"(' eps_cdet = ',g12.6)") eps_cdet
write (out,"(' dont_check_s2 = ',l8)") dont_check_s2
write (out,"(' own_integrals = ',l8)") own_integrals
!
if (.not.own_integrals) then
write (out,"(' AO file = ',a)") trim(file_sao)
end if
write (out,"(' Bra MO file = ',a)") trim(file_cbra)
write (out,"(' Ket MO file = ',a)") trim(file_cket)
write (out,"(' Bra determinants = ',a)") trim(file_detbra)
write (out,"(' Ket determinants = ',a)") trim(file_detket)
end if
!
! If we are using the new format, we'll have to parse checkpoint
! files before doing anything else.
!
if (own_integrals) then
call gamess_load_orbitals(file=trim(file_cbra),structure=g_bra)
call gamess_load_orbitals(file=trim(file_cket),structure=g_ket)
naos_bra = g_bra%nbasis
naos_ket = g_ket%nbasis
write (out,"('Number of AOs in the reference: ',i6)") naos_bra
write (out,"(' Number of AOs in the ion: ',i6)") naos_ket
end if
!
! How many components does our operator have?
!
select case (braket_operator)
case default
op_components = 1
case ('none')
op_components = 0
have_operator = .false.
case ('dipole','velocity')
op_components = 3
case ('r-d/dr')
op_components = 9
end select
!
! Allocate storage for orbital quantities now.
!
allocate (occbra(nmobra,ndetbra),cbra(naos_bra,nmobra),cdetref(ndetbra))
allocate (occket(nmoket,ndetket),cket(naos_ket,nmoket),cdetion(ndetket))
allocate (sao(naos_bra,naos_ket),dipao(naos_bra,naos_ket,op_components))
allocate (sao_ref(naos_bra,naos_bra),sao_ion(naos_ket,naos_ket))
allocate (sx(nmobra,nmoket),dipx(nmobra,nmoket,op_components))
allocate (aolabels(naos_bra))
!
! 1e integrals in the AO basis
!
if (own_integrals) then
have_overlaps = .true.
call gamess_1e_integrals('AO OVERLAP', sao, g_bra,g_ket)
if (have_operator) then
select case (braket_operator)
case default
write (out,"('initialize_core_data: operator ',a,' is not recognized')") trim(braket_operator)
stop 'sd_core%initialize_core_data - bad operator'
case ('kinetic')
call gamess_1e_integrals('AO KINETIC', dipao(:,:,1),g_bra,g_ket)
case ('dipole')
call gamess_1e_integrals('AO DIPOLE X', dipao(:,:,1),g_bra,g_ket)
call gamess_1e_integrals('AO DIPOLE Y', dipao(:,:,2),g_bra,g_ket)
call gamess_1e_integrals('AO DIPOLE Z', dipao(:,:,3),g_bra,g_ket)
case ('velocity')
call gamess_1e_integrals('AO D/DX', dipao(:,:,1),g_bra,g_ket)
call gamess_1e_integrals('AO D/DY', dipao(:,:,2),g_bra,g_ket)
call gamess_1e_integrals('AO D/DZ', dipao(:,:,3),g_bra,g_ket)
case ('r-d/dr')
fill_rddr: do jc=1,3
do ic=1,3
call gamess_1e_integrals('AO R-D/DR',dipao(:,:,ic+3*(jc-1)),g_bra,g_ket,op_index=(/ic,jc/))
end do
end do fill_rddr
case ('1/r')
call gamess_1e_integrals('AO 3C 1/R', dipao(:,:,1),g_bra,g_ket,op_xyz=op_centre)
case ('r')
call gamess_1e_integrals('AO 3C R', dipao(:,:,1),g_bra,g_ket,op_xyz=op_centre)
case ('r**2')
call gamess_1e_integrals('AO 3C R**2', dipao(:,:,1),g_bra,g_ket,op_xyz=op_centre)
case ('gauss')
call gamess_1e_integrals('AO 3C GAUSS', dipao(:,:,1),g_bra,g_ket,op_xyz=op_centre,op_param=(/op_gauss/))
case ('gauss r')
call gamess_1e_integrals('AO 3C GAUSS R',dipao(:,:,1),g_bra,g_ket,op_xyz=op_centre,op_param=(/op_gauss/))
end select
end if
call gamess_1e_integrals('AO OVERLAP', sao_ref, g_bra,g_bra)
call gamess_1e_integrals('AO OVERLAP', sao_ion, g_ket,g_ket)
aolabels(:) = g_bra%bas_labels
else ! .not.own_integrals; can only handle overlaps
have_overlaps = .true.
call read_AO_integrals(file_sao,' OVERLAP MATRIX',sao,aolabels,found=have_overlaps)
sao_ref = sao
sao_ion = sao
!
! Dipole integrals may or may not be present on input, depending on the
! type of caluclation used to generate the output. If no dipoles are
! provided, we won't be able to calculate transition matrix elements
! or exchange corrections.
!
braket_operator = 'dipole'
have_operator = .true.
call read_AO_integrals(file_sao,' X DIPOLE INTEGRALS',dipao(:,:,1),found=have_operator)
call read_AO_integrals(file_sao,' Y DIPOLE INTEGRALS',dipao(:,:,2),found=have_operator)
call read_AO_integrals(file_sao,' Z DIPOLE INTEGRALS',dipao(:,:,3),found=have_operator)
end if ! own_integrals
!
! MO coeffcients. If using the new format, these have been loaded before.
! Otherwise, read them now.
!
if (own_integrals) then
if (nmobra>g_bra%nvectors) then
write (out,"('Need ',i6,' MOs in the reference, but only ',i6,' were found in ',a)") &
nmobra, g_bra%nvectors, trim(file_cbra)
stop 'sd_core%read_core_input - not enough ref MOs'
end if
if (nmoket>g_ket%nvectors) then
write (out,"('Need ',i6,' MOs in the ion, but only ',i6,' were found in ',a)") &
nmoket, g_ket%nvectors, trim(file_cket)
stop 'sd_core%read_core_input - not enough ion MOs'
end if
cbra = g_bra%vectors(:,1:nmobra)
cket = g_ket%vectors(:,1:nmoket)
else ! .not.own_integrals
call read_orbitals(file_cbra,nmobra,cbra)
call read_orbitals(file_cket,nmoket,cket)
end if
!
! If orbitals are identical, I'd like to know this.
!
call compare_orbitals(same_mos)
if (.not.same_mos .and. .not.have_overlaps) then
write (out,"('Reference and ion wavefunctions do not share the MO space,"// &
" and no overlap integrals are available')")
stop 'insufficient data'
end if
!
! Ready to read in the reterminants
!
open(ioscr,file=file_detbra,action='read',status='old',recl=lchar,iostat=status)
if (status/=0) then
write (out,"('Cannot open reference (bra) wavefunction ',a)") trim(file_detbra)
stop 'file_detbra'
end if
read_ref_dets: do detref=1,ndetbra
read(ioscr,*) cdetref(detref), occbra(:,detref)
end do read_ref_dets
close (ioscr)
open(ioscr,file=file_detket,action='read',status='old',recl=lchar,iostat=status)
if (status/=0) then
write (out,"('Cannot open ion (ket) wavefunction ',a)") trim(file_detket)
stop 'file_detket'
end if
read_ion_dets: do detion=1,ndetket
read(ioscr,*) cdetion(detion), occket(:,detion)
end do read_ion_dets
close (ioscr)
call TimerStop('initialize_core_data')
!
! We can now finally know the number of electrons ...
!
nelbra = sum(abs(occbra(:,1)))
nelket = sum(abs(occket(:,1)))
!
! Precompute sorted spin-orbital lists for symmetry handling
!
if (use_symmetry) then
call TimerStart('initialize_symmetry')
allocate (scnt_bra(2,ndetbra),spinbra(nelbra,2,ndetbra),ord_bra(nelbra,2,ndetbra))
allocate (scnt_ket(2,ndetket),spinket(nelket,2,ndetket),ord_ket(nelket,2,ndetket))
call symmetry_sort_spinorbitals(nelbra,ndetbra,occbra,scnt_bra,spinbra,ord_bra)
call symmetry_sort_spinorbitals(nelket,ndetket,occket,scnt_ket,spinket,ord_ket)
call TimerStop('initialize_symmetry')
end if
end subroutine initialize_core_data
!
subroutine check_input_consistency
integer(ik) :: detref, detion
!
call TimerStart('check_input_consistency')
!
! Check that determinants have matching number of electrons, and
! that ion has one less electron than the reference
!
check_ref_dets: do detref=1,ndetbra
if ( any(occbra(:,detref)/=0 .and. occbra(:,detref)/=2 .and. abs(occbra(:,detref))/=1) ) then
write (out,"('Reference determinant ',i6,' has invalid occupations')") detref
stop 'bad ref det'
end if
if (sum(abs(occbra(:,detref)))/=nelbra) then
write (out,"('Reference determinant ',i6,' has invalid electron count')") detref
stop 'bad ref nel'
end if
end do check_ref_dets
!
check_ion_dets: do detion=1,ndetket
if ( any(occket(:,detion)/=0 .and. occket(:,detion)/=2 .and. abs(occket(:,detion))/=1) ) then
write (out,"('Reference determinant ',i6,' has invalid occupations')") detion
stop 'bad ref det'
end if
if (sum(abs(occket(:,detion)))/=nelket) then
write (out,"('Reference determinant ',i6,' has invalid electron count')") detion
stop 'bad ref nel'
end if
end do check_ion_dets
!
write (out,"('Number of electrons in reference and ion: ',2i6)") nelbra, nelket
if (nelket ==nelbra) then
write (out,"('Overlap mode')")
mode_overlap = .true.
else if (nelket+1==nelbra) then
write (out,"('Dyson mode')")
mode_overlap = .false.
else
stop 'bad electron counts'
end if
!
! Check normalization of the reference and ions' MOs - this should
! guard against any stupid errors in reading the data.
!
call check_orthonormality('parent',nmobra,cbra,sao_ref)
call check_orthonormality('ion ',nmoket,cket,sao_ion)
!
! Check total wavefunction norm and s^2 expectation value
!
if (dont_check_s2) then
write (out,"(/'WARNING: CHECK FOR S^2 AND NORMALIZATION IS OFF'/)")
else
call check_s2('parent',cdetref,occbra)
call check_s2('ion ',cdetion,occket)
end if
!
call TimerStop('check_input_consistency')
end subroutine check_input_consistency
!
! Find the extent of the frozen core, taking into account the MO overlap
! integrals and the determinant populations.
!
subroutine locate_frozen_core
integer(ik) :: ref_docc, ion_docc ! Number of doubly occupied MOs
integer(ik) :: ndocc, imo
real(rk) :: err
!
call TimerStart('locate_frozen_core')
ref_docc = count_doubly_occupied(occbra)
ion_docc = count_doubly_occupied(occket)
ndocc = min(ref_docc,ion_docc)
scan_overlap: do imo=1,ndocc
err = abs(sx(imo,imo)-1.0_rk)
if (err>spacing(1e3_rk)) then
ndocc = imo - 1
exit scan_overlap
end if
end do scan_overlap
frozen_orbitals = ndocc
write (out,"(' Doubly occupied reference orbitals: ',i6)") ref_docc
write (out,"(' Doubly occupied ion orbitals: ',i6)") ion_docc
write (out,"(' Frozen core orbitals: ',i6)") frozen_orbitals
call TimerStop('locate_frozen_core')
end subroutine locate_frozen_core
!
! Count the number of doubly-occupied MOs in all determinants in the list
!
function count_doubly_occupied(dets) result(ndocc)
integer(ik), intent(in) :: dets(:,:)
integer(ik) :: ndocc
!
integer(ik) :: nmos, ndets, idet, imo
!
nmos = size(dets,dim=1)
ndets = size(dets,dim=2)
ndocc = nmos
scan_dets: do idet=1, ndets
scan_mos: do imo=1, nmos
if (dets(imo,idet)/=2) then
ndocc = min(ndocc,imo-1)
exit scan_mos
end if
end do scan_mos
end do scan_dets
end function count_doubly_occupied
!
! Read one-electron integrals from GAMESS NPRINT=3 output
!
subroutine read_AO_integrals(name,section,xao,aolabels,found)
character(len=*), intent(in) :: name ! File containing integrals
character(len=*), intent(in) :: section ! Section header
real(rk), intent(out) :: xao(:,:) ! Matrix containing the same
character(len=*), intent(out), optional :: aolabels(:) ! AO labels, for human-readable output
logical(lk), intent(inout), optional :: found ! Set to .FALSE. if section was not found
!
character(len=lchar) :: buf
character(len=lchar) :: errbuf
integer(ik) :: status, file_line
logical(lk) :: go
integer(ik) :: row, column, base_column
!
if (verbose>=1) write (out,"('Reading integral section ',a,' from ',a)") trim(section), trim(name)
file_line = 0
error_block: do
errbuf = 'Error opening '//trim(name)//' for reading'
open(ioscr,file=name,action='read',status='old',recl=lchar,iostat=status)
if (status/=0) exit error_block
!
! Skip until the header in front of the overlap integrals matrix
!
errbuf = 'Error searching for "'//trim(section)//'" in '//trim(name)
go = .true.
skip_header: do while(go)
buf = ' '
file_line = file_line + 1
read(ioscr,"(a)",iostat=status) buf
if (status/=0) exit error_block
go = (buf/=section)
end do skip_header
!
! Read the integrals. The whole thing will end with a line
! containing a word 'INTEGRALS'
!
errbuf = 'Error processing '//trim(section)//' from '//trim(name)
go = .true.
read_integrals: do while(go)
buf = ' '
file_line = file_line + 1
read(ioscr,"(a)",iostat=status) buf
if (status/=0) exit error_block
if (buf==' ') cycle read_integrals
!
if (index(buf,'INTEGRALS')/=0) then
go = .false.
cycle read_integrals
end if
!
if (buf(1:10)==' ') then
!
! This is the header - we'll read the first integer, which
! is our base column for the next few lines.
!
read(buf,*,iostat=status) base_column
if (status/=0) exit error_block
if (base_column<=0 .or. base_column>naos_bra) then ! Old-style; naos_bra == naos_ket
errbuf = 'Invalid base column value'
exit error_block
end if
else
!
! Integral lines. The first integer is the orbital index; then
! we have labels (to be ignored) and the overlap integrals.
!
read(buf,"(i5)",iostat=status) row
if (status/=0) exit error_block
if (row<=0 .or. row>naos_bra) then ! Old-style; naos_bra == naos_ket
errbuf = 'Invalid row value'
exit error_block
end if
if (base_column==1 .and. present(aolabels)) then
read(buf,"(5x,a11)",iostat=status) aolabels(row)
end if
read(buf,"(15x,5f11.6)",iostat=status) &
(xao(row,column),column=base_column,min(row,base_column+4))
if (status/=0) exit error_block
!
! Add the transpose
!
xao(base_column:min(row,base_column+4),row) = xao(row,base_column:min(row,base_column+4))
end if
end do read_integrals
!
! Close the file, and we are done
!
errbuf = 'Error closing '//trim(name)
close(ioscr,iostat=status)
if (status/=0) exit error_block
return
end do error_block
write (out,"(' Line ',i10,': ',a)") file_line, trim(errbuf)
if (present(found)) then
found = .false.
else
stop 'read_AO_integrals'
end if
end subroutine read_AO_integrals
!
! Load a set of molecular orbitals.
!
subroutine read_orbitals(name,nmos,c)
character(len=*), intent(in) :: name ! Name of the file
integer(ik), intent(in) :: nmos ! Number of orbitals to read
real(rk), intent(out) :: c(:,:) ! Space for the orbitals
!
integer(ik) :: status, file_line
character(len=lchar) :: errbuf
integer(ik) :: imo, minao, maxao
integer(ik) :: chk_mo, chk_line, line
!
if (verbose>=1) write (out,"('Reading orbitals from ',a)") trim(name)
file_line = 0
error_block: do
errbuf = 'Error opening '//trim(name)//' for reading'
open(ioscr,file=name,action='read',status='old',recl=lchar,iostat=status)
if (status/=0) exit error_block
!
errbuf = 'Error reading orbitals from '//trim(name)
mos_loop: do imo=1,nmos
!
! The orbitals are split over a number of lines, in a rather
! ugly format - so we'll have to jump through a few hoops.
!
minao = 1
line = 1
line_loop: do while(minao<=naos_bra) ! Note that with the old input, naos_bra and naos_ket are identical
maxao = min(minao+4,naos_bra)
file_line = file_line + 1
read(ioscr,"(i2,i3,5e15.0)",iostat=status) chk_mo, chk_line, c(minao:maxao,imo)
if (status/=0) exit error_block
if (mod(chk_mo,100)/=mod(imo,100) .or. chk_line/=line) then
write (errbuf,"('In file ',a,' chk_mo = ',i5,' imo = ',i5,' chk_line = ',i5,' line = ',i5)") &
trim(name), chk_mo, imo, chk_line, line
exit error_block
end if
minao = maxao + 1
line = line + 1
end do line_loop
end do mos_loop
!
errbuf = 'Error closing '//trim(name)
close(ioscr,iostat=status)
return
end do error_block
write (out,"(' Line ',i10,': ',a)") file_line, trim(errbuf)
stop 'read_orbitals'
end subroutine read_orbitals
!
! Compare orbital sets for being identical - this permits many short-cuts
!
subroutine compare_orbitals(same)
logical(lk), intent(out) :: same ! Orbital sets are identical
!
integer(ik) :: nmo, iat, ish, is, p1, p2, p1i, p2i
real(rk) :: diff, maxc
!
same = .false.
same_basis = .false.
check_orbitals: do
!
! If using the new input, we have to make sure that geometries and the
! number of AOs are the same.
!
if (own_integrals) then
write (out,"(' Reference atoms = ',i8)") g_bra%natoms
write (out,"(' Ion atoms = ',i8)") g_ket%natoms
if (g_bra%natoms /= g_ket%natoms) exit check_orbitals
!
scan_atoms: do iat=1,g_bra%natoms
diff = maxval(abs(g_bra%atoms(iat)%xyz-g_ket%atoms(iat)%xyz))
if (diff>spacing(100._rk)) then
write (out,"('Reference atom ',i4,' is at ',3f16.8)") iat, g_bra%atoms(iat)%xyz
write (out,"(' Ion atom ',i4,' is at ',3f16.8)") iat, g_ket%atoms(iat)%xyz
exit check_orbitals
end if
if (g_bra%atoms(iat)%nshell/=g_ket%atoms(iat)%nshell) then
write (out,"('Reference atom ',i4,' has ',i4,' basis shells')") iat, g_bra%atoms(iat)%nshell
write (out,"(' Ion atom ',i4,' has ',i4,' basis shells')") iat, g_ket%atoms(iat)%nshell
exit check_orbitals
end if
scan_shells: do is=1,g_bra%atoms(iat)%nshell
if ( g_bra%atoms(iat)%sh_l(is) /= g_ket%atoms(iat)%sh_l(is) ) then
write (out,"('Reference atom ',i4,' shell ',i3,' has angular momentum ',i2)") iat, is, g_bra%atoms(iat)%sh_l(is)
write (out,"(' Ion atom ',i4,' shell ',i3,' has angular momentum ',i2)") iat, is, g_ket%atoms(iat)%sh_l(is)
exit check_orbitals
end if
p1 = g_bra%atoms(iat)%sh_p(is)
p1i = g_ket%atoms(iat)%sh_p(is)
p2 = g_bra%atoms(iat)%sh_p(is+1) - 1
p2i = g_ket%atoms(iat)%sh_p(is+1) - 1
if ( p1/=p1i .or. p2/=p2i ) then
write (out,"('Reference atom ',i4,' shell ',i4,' spans primitives ',i4,' - ',i4)") iat, is, p1, p2
write (out,"(' Ion atom ',i4,' shell ',i4,' spans primitives ',i4,' - ',i4)") iat, is, p1i, p2i
exit check_orbitals
end if
maxc = max(maxval(abs(g_bra%atoms(iat)%p_zet(p1:p2))),maxval(abs(g_ket%atoms(iat)%p_zet(p1:p2))))
diff = maxval(abs(g_bra%atoms(iat)%p_zet(p1:p2)-g_ket%atoms(iat)%p_zet(p1:p2)))
if (diff>100._rk*maxc) then
write (out,"('Reference and ion atoms ',i4,' shell ',i4,' have different primitive exponents')") iat, ish
exit check_orbitals
end if
maxc = max(maxval(abs(g_bra%atoms(iat)%p_c(p1:p2))),maxval(abs(g_ket%atoms(iat)%p_c(p1:p2))))
diff = maxval(abs(g_bra%atoms(iat)%p_c(p1:p2)-g_ket%atoms(iat)%p_c(p1:p2)))
if (diff>100._rk*maxc) then
write (out,"('Reference and ion atoms ',i4,' shell ',i4,' have different contraction weights')") iat, ish
exit check_orbitals
end if
end do scan_shells
end do scan_atoms
end if
!
same_basis = .true.
nmo = min(nmoket,nmobra)
maxc = max(maxval(abs(cbra)),maxval(abs(cket)))
diff = maxval(abs(cbra(:,:nmo)-cket(:,:nmo)))
same = (diff <= 100._rk * spacing(maxc))
write (out,"('Comparing ',i6,' lowest MOs out of ',i6,'/',i6)") nmo, nmobra, nmoket
write (out,"(' Max absolute coefficient = ',g12.6)") maxc
write (out,"('Max coefficient difference = ',g12.6)") diff
exit check_orbitals
end do check_orbitals
write (out,"(' The two sets match = ',l8)") same
end subroutine compare_orbitals
!
! Calculate and check overlaps between MOs - we must get a unit matrix
! AO overlaps are in sao(:,:).
!
subroutine check_orthonormality(title,nmos,mos,sao)
character(len=*), intent(in) :: title ! For error messages
integer(ik), intent(in) :: nmos ! Number of MOs
real(rk), intent(in) :: mos(:,:) ! MO coefficients
real(rk), intent(in) :: sao(:,:) ! Overlap matrix (possibly ion or ref-specific)
!
integer(ik) :: lmo, rmo ! MO indices
real(rk) :: ov ! Calculated MO overlap
real(rk) :: refval ! Expected MO overlap
!
if (.not.have_overlaps) then
write (out,"('Atomic overlap integrals not present; unable to check '" // &
",a,' orbitals for orthonormality')") trim(title)
return
end if
!
left_scan: do lmo=1,nmos
right_scan: do rmo=1,nmos
ov = dot_product(mos(:,lmo),matmul(sao,mos(:,rmo)))
!
refval = 0._rk
if (lmo==rmo) refval = 1._rk
!
if (abs(ov-refval)>eps_s) then
write (out,"(5x,'warning: ',a,' <',i5,'|',i5,'> = ',g16.9,' (expected ',g16.9,')')") &
trim(title), lmo, rmo, ov, refval
end if
end do right_scan
end do left_scan
!
end subroutine check_orthonormality
!
! Check norm and S^2
!
subroutine check_s2(title,c,d)
character(len=*), intent(in) :: title ! Name of the w.f.
real(rk), intent(in) :: c( :) ! Amplitudes
integer(ik), intent(in) :: d(:,:) ! Determinants
!
real(rk) :: ov, s2, s2exact
integer(ik) :: is
!
call calculate_s2(c,d,ov,s2)
!
write (out,"(/1x,' ',a,' norm is: ',f12.8 )") trim(title), ov
write (out,"( 1x,a,' S^2 expectation is: ',f12.8, '(',f12.8,')'/)") trim(title), s2, s2/ov
!
if (abs(ov-1._rk)>0.03_rk) then
write (out,"('sd_core%check_s2 - unnormalized wavefunction')")
stop 'sd_core%check_s2 - unnormalized wavefunction'
end if
!
scan_spins: do is=0,10
s2exact = 0.25_rk * is * (is+2)
if (abs(s2exact-s2/ov)<0.03_rk) return
end do scan_spins
stop 'sd_core%check_s2 - wavefunction is not an S^2 eigenfunction'
end subroutine check_s2
!
! Give user some idea as to how similar the MOs are between parent
! and the ion.
!
subroutine report_mo_correlations(smo)
real(rk), intent(in) :: smo(:,:) ! Cross-system overlap integrals, MO basis
!
integer(ik) :: refmo, simmo
!
write (out,"(/t5,'Reference - ion MO correlations')")
write (out,"(2x,a7,2x,a7,2x,a12)") 'Ref. MO', 'Ion MO', 'Overlap', &
'-------', '------', '-------'
mo_ref: do refmo=1,nmobra
simmo = maxloc(abs(smo(refmo,:)),dim=1)
write (out,"(2x,i7,2x,i7,2x,f12.5)") refmo, simmo, smo(refmo,simmo)
end do mo_ref
write (out,"()")
end subroutine report_mo_correlations
!
! Construct a one-particle operator matrix elements in cross-state
! molecular orbital basis. Reference state is on the left; target
! is on the right; all alpha-spin appears before all beta-spin
! spin-orbitals.
!
subroutine mat_ao_to_mo(title,xmo,xao)
character(len=*), intent(in) :: title ! Title; identification only
real(rk), intent(out) :: xmo(:,:) ! Matrix elements in the MO basis
real(rk), intent(in) :: xao(:,:) ! Matrix elements in the AO basis
!
integer(ik) :: refmo, ionmo ! Orbitals of the reference and ion
real(rk) :: val
!
call TimerStart('mat_ao_to_mo')
!
! This is a very inefficient way of making the transformation.
! Much better perfomance can be obtained with MATMUL and a buffer or DGEMM
!
if (verbose>=2) write (out,"(/'Cross-state ',a,' elements')") trim(title)
mo_ref: do refmo=1,nmobra
mo_ion: do ionmo=1,nmoket
val = dot_product(cbra(:,refmo),matmul(xao,cket(:,ionmo)))
if (verbose>=2) then
write (out,"(5x,a,' <',i5,'|',i5,'> = ',f25.12)") trim(title), refmo, ionmo, val
end if
xmo(refmo,ionmo) = val
end do mo_ion
end do mo_ref
if (verbose>=2) write (out,"()")
!
call TimerStop('mat_ao_to_mo')
end subroutine mat_ao_to_mo
!
! Transform 1-electron quantity back from MO to the AO basis.
! Should reverse the action mat_ao_to_mo
!
subroutine mat_mo_to_ao(title,xmo,xao)
character(len=*), intent(in) :: title ! Title; identification only
real(rk), intent(in) :: xmo(:,:) ! Matrix elements in the MO basis
real(rk), intent(out) :: xao(:,:) ! Matrix elements in the AO basis
!
integer(ik) :: refao, ionao ! Atomic orbitals of the reference and the ion
real(rk) :: val
!
call TimerStart('mat_mo_to_ao')
!
! This is a very inefficient way of making the transformation.
! Much better perfomance can be obtained with MATMUL and a buffer or DGEMM
!
if (verbose>=2) write (out,"(/'Cross-state ',a,' elements')") trim(title)
ao_ref: do refao=1,naos_bra
ao_ion: do ionao=1,naos_ket
val = dot_product(cbra(refao,:),matmul(xmo,cket(ionao,:)))
if (verbose>=2) then
write (out,"(5x,a,' <',i5,'|',i5,'> = ',f25.12)") trim(title), refao, ionao, val
end if
xao(refao,ionao) = val
end do ao_ion
end do ao_ref
if (verbose>=2) write (out,"()")
!
call TimerStop('mat_mo_to_ao')
end subroutine mat_mo_to_ao
!
! Initialize unit matrix
!
subroutine mat_unit(m)
real(rk), intent(out) :: m(:,:)
integer(ik) :: i
m = 0._rk
do i=1,min(size(m,dim=1),size(m,dim=2))
m(i,i) = 1._rk
end do
end subroutine mat_unit
!
! Calculate the distance between two determinants
!
function determinant_distance(ref,ion) result(n)
integer(ik), intent(in) :: ref(:), ion(:)
integer(ik) :: n
!
integer(ik) :: len_r, len_i, len_c, idet
integer(ik) :: nr, ni
!
len_r = size(ref)
len_i = size(ion)
len_c = min(len_r,len_i)
n = 0
sum_common: do idet=1,len_c
nr = ref(idet)
ni = ion(idet)
if (nr==ni) cycle sum_common
!
if (abs(nr)/=1) then
n = n + abs(nr - abs(ni))
else if (abs(ni)/=1) then
n = n + abs(ni - abs(nr))
else