-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathimport_gamess.f90
2817 lines (2815 loc) · 136 KB
/
import_gamess.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
!
! Import simple GAMESS checkpoint files.
!
! Restrictions:
! Only C1 symmetry (no symmetry) is allowed.
! Basis sets must be in-line.
! Only simple, single-SCF files are supported.
! Only one $VEC section is supported.
!
! In most cases, these restrictions mean that the checkpoint file
! will have to be edited by hand to be loadable by this import filter.
!
! The module is derived from my (much more capable) gamess-dx.c filter.
!
! Added Feb 2010 (Michael Spanner): Support for loading natural orbital occupation numbers
! and rdm_sv values from the GAMESS checkpoint files.
! Added March 2010 (Michael Spanner): Support for rotating the gamess data relative to multigrid
!
! Added Aug 2011 (SP): Support for 3-centre 1e integral evaluation. The 3-centre integrals
! implementation is essentiall Obara-Saika, with the extensions to
! the supported operator classes based on Ahlrichs PCCP 8, 3072 (2006).
! Note that OS and A use slightly different definitions of the primitive
! integrals, which adds a tiny bit of complexity.
!
! WARNING: Non-unit rotation matrix is only implemented for
! evaluation of orbitals on grid. Obara-Saika integrals
! or ECP integrals will prodice incorrect results.
!
! WARNING: R/(R**2+A**2) and A/(R**2+A**2) integrals have limited range and accuracy.
! Please make sure you understand the limits when you use those!
!
! Added May 2012 (SP): Support for 2e integral evaluation, following Ahlrichs' PCCP.
! Support for kinetic-energy integrals, following Obara & Saika.
!
! WARNING: 2e integral implementation is extremely inefficient.
!
! Added Dec 21, 2012 (SP): Support for matrix elements of the planewave
!
! Feb 06, 2013 (SP): Changed batching of the integrals, to improve 2e efficiency
! for heavily-contracted basis sets (think correlation-consistent
! basis sets).
!
! Jul 17, 2013 (SP): Added support for r_i (d/d r_j) integrals.
!
! Apr 23, 2014 (SP): Added support for (most) integrals in quadruple precision.
!
! In order to remain compatible with historical usage in multigrid, we
! maintain a default structure context; all calls which do not specify
! (one of the) context(s) default to that one.
!
module import_gamess
use accuracy
use constants
use math
use lapack
use timer
use gamess_internal
use os_integral_operators
!$ use OMP_LIB
implicit none
private
! public os_basic_integral ! Debugging!
public gamess_load_orbitals
public gamess_evaluate_functions
public gamess_oversample_functions
public gamess_load_natocc
public gamess_load_rdmsv
public gamess_report_nuclei
public gamess_destroy
public gamess_1e_integrals
public gamess_print_1e_integrals
public gam_atom, gam_structure
public gamess_oversample
public gamess_2e_info
public gamess_2e_integrals
public rcsid_import_gamess
!
character(len=clen), save :: rcsid_import_gamess = "$Id: import_gamess.f90,v 1.15 2022/08/02 08:51:57 ps Exp $"
!
! Conditional compilation is needed for the interfaces:
! the _quad versions differ only in their argument kinds; if we
! compile with real and quad kinds being the same, the interface
! part will fail!
!
interface gamess_print_1e_integrals
module procedure gamess_print_1e_integrals_real
module procedure gamess_print_1e_integrals_complex
!*qd module procedure gamess_print_1e_integrals_quad
end interface gamess_print_1e_integrals
!
interface gamess_1e_integrals
module procedure gamess_1e_integrals_real
module procedure gamess_1e_integrals_complex
!*qd module procedure gamess_1e_integrals_quad
end interface gamess_1e_integrals
!
interface gamess_2e_integrals
module procedure gamess_2e_integrals_real
!*qd module procedure gamess_2e_integrals_quad
end interface gamess_2e_integrals
!
! Interfaces below are for the internal routines; they exist to make it possible to
! share code between different-precision versions of the routines.
!
interface os_1e_matrix
module procedure os_1e_matrix_real
module procedure os_1e_matrix_complex
!*qd module procedure os_1e_matrix_quad
end interface os_1e_matrix
!
interface os_1e_contraction
module procedure os_1e_contraction_real
module procedure os_1e_contraction_complex
!*qd module procedure os_1e_contraction_quad
end interface os_1e_contraction
!
interface os_1e_overlap
module procedure os_1e_overlap_real
!*qd module procedure os_1e_overlap_quad
end interface os_1e_overlap
!
interface os_1e_dipole
module procedure os_1e_dipole_real
!*qd module procedure os_1e_dipole_quad
end interface os_1e_dipole
!
interface os_1e_grad
module procedure os_1e_grad_real
!*qd module procedure os_1e_grad_quad
end interface os_1e_grad
!
interface os_1e_rddr
module procedure os_1e_rddr_real
!*qd module procedure os_1e_rddr_quad
end interface os_1e_rddr
!
interface os_1e_kinetic
module procedure os_1e_kinetic_real
!*qd module procedure os_1e_kinetic_quad
end interface os_1e_kinetic
!
interface os_1e_3c
module procedure os_1e_3c_real
!*qd module procedure os_1e_3c_quad
end interface os_1e_3c
!
interface os_2e_batch
module procedure os_2e_batch_real
!*qd module procedure os_2e_batch_quad
end interface os_2e_batch
!
interface os_2e_primitives
module procedure os_2e_primitives_real
!*qd module procedure os_2e_primitives_quad
end interface os_2e_primitives
!
interface os_basic_integral
module procedure os_basic_integral_real
!*qd module procedure os_basic_integral_quad
end interface os_basic_integral
!
interface os_boys_table
module procedure os_boys_table_real
!*qd module procedure os_boys_table_quad
end interface os_boys_table
!
interface os_common_primitives
module procedure os_common_primitives_real
!*qd module procedure os_common_primitives_quad
end interface os_common_primitives
!
interface gamess_evaluate_functions
module procedure gamess_evaluate_functions_full
module procedure gamess_evaluate_functions_simple
end interface gamess_evaluate_functions
!
! Internal data structures and magic constants
!
integer(ik), parameter :: verbose = 0 ! 0 is too verbose for reading files repeatedly
!
logical, save :: oversample = .true. ! Activate oversampling
!
type(gam_structure), save, target :: gam_def ! Default GAMESS structure description
!
contains
!
! External interfaces
!
subroutine gamess_oversample(sample)
logical, intent(in) :: sample
!
oversample = sample
end subroutine gamess_oversample
!
! Report positions and charges of the nuclei loaded alongside the MOs
!
subroutine gamess_report_nuclei(natom,xyzq,structure,true_charge)
integer(ik), intent(out) :: natom ! Number of atoms in the molecule loaded by the
! last call to gamess_load_orbitals, or 0 if none
real(rk), intent(out), optional :: xyzq(:,:) ! Coordinates and charges of the atoms
type(gam_structure), intent(in), &
target, optional :: structure ! Context to use
logical, intent(in), optional :: true_charge ! Report true nuclear charge, including the ECP charge
!
type(gam_structure), pointer :: gam
integer(ik) :: na
logical :: add_ecp
!
if (present(structure)) then
gam => structure
else
gam => gam_def
end if
natom = gam%natoms
!
add_ecp = .false.
if (present(true_charge)) add_ecp = true_charge
!
if (present(xyzq)) then
if ( size(xyzq,dim=1)<4 .or. size(xyzq,dim=2)<gam%natoms ) then
stop 'import_gamess%gamess_report_nuclei - buffer too small'
end if
xyzq(1,1:gam%natoms) = real(gam%atoms(1:gam%natoms)%xyz(1),kind=kind(xyzq)) / abohr
xyzq(2,1:gam%natoms) = real(gam%atoms(1:gam%natoms)%xyz(2),kind=kind(xyzq)) / abohr
xyzq(3,1:gam%natoms) = real(gam%atoms(1:gam%natoms)%xyz(3),kind=kind(xyzq)) / abohr
xyzq(4,1:gam%natoms) = real(gam%atoms(1:gam%natoms)%znuc,kind=kind(xyzq))
if (add_ecp) then
xyzq(4,1:gam%natoms) = xyzq(4,1:gam%natoms) + nint(gam%atoms(1:gam%natoms)%ecp_zcore)
end if
!
! Rotate coordinates
!
rotate_atoms: do na=1,gam%natoms
xyzq(1:3,na) = real(matmul(gam%rotmat,xyzq(1:3,na)),kind=kind(xyzq))
end do rotate_atoms
end if
end subroutine gamess_report_nuclei
!
! Destroy data context
!
subroutine gamess_destroy(structure)
type(gam_structure), intent(inout), target, optional :: structure ! Context to use
!
type(gam_structure), pointer :: gam
integer(ik) :: alloc
!
if (present(structure)) then
gam => structure
else
gam => gam_def
end if
if (allocated(gam%atoms)) then
deallocate(gam%atoms,stat=alloc)
if (alloc/=0) then
write (out,"('gamess_import%gamess_destroy: Error ',i6,' deallocating atoms table')") alloc
stop 'gamess_import%gamess_destroy: Error deallocating atoms table'
end if
end if
if (allocated(gam%bas_labels)) then
deallocate(gam%bas_labels,stat=alloc)
if (alloc/=0) then
write (out,"('gamess_import%gamess_destroy: Error ',i6,' deallocating basis labels table')") alloc
stop 'gamess_import%gamess_destroy: Error deallocating basis labels table'
end if
end if
if (allocated(gam%vectors)) then
deallocate(gam%vectors,stat=alloc)
if (alloc/=0) then
write (out,"('gamess_import%gamess_destroy: Error ',i6,' deallocating vectors table')") alloc
stop 'gamess_import%gamess_destroy: Error deallocating vectors table'
end if
end if
if (allocated(gam%batches)) then
deallocate(gam%batches,stat=alloc)
if (alloc/=0) then
write (out,"('gamess_import%gamess_destroy: Error ',i6,' deallocating batches table')") alloc
stop 'gamess_import%gamess_destroy: Error deallocating batches table'
end if
end if
gam%natoms = 0
gam%nbasis = 0
gam%nvectors = 0
gam%nbatches = 0
end subroutine gamess_destroy
!
! Load structure file and evaluate the MOs. This subroutine is a little
! complicated, since it can perform three separate functions:
!
! 1. Parse GAMESS checkpoint file, and remember the result. In this case,
! it should be called with the argument "file" (and possibly "structure")
! 2. Evaluate orbitals on grid. In this case, the arguments "mos", "dst",
! "coord", and "grid" (and possibly "structure") should all be present,
! but "file" should not be.
! 3. Perform (1), then (2), then deallocate the orbitals. In this case, all
! arguments (except possibly "structure") should be present.
!
! Historically, only the #3 semantics was used in multigrid.
!
! See also gamess_evaluate_functions() and gamess_oversample_functions
! for less integrated interfaces
!
subroutine gamess_load_orbitals(file,mos,dst,coord,grid,structure,rot,dx,max_mos)
character(len=*), intent(in), optional :: file ! Name of the GAMESS checkpoint file. WARNING: Presence of the
! file= argument will affect semantics of the rot= argument
! below
integer(ik), intent(in), optional :: mos(:) ! Indices of the orbitals to be computed
integer(ik), intent(in), optional :: dst(:) ! Last indices of the grid array, receiving
! the MOs.
real(rk), intent(in), optional :: coord(:,:,:,:) ! Coordinates of the points, at which MOs
! must be evaluated. First index is X,Y,Z.
complex(rk), intent(out), optional :: grid(:,:,:,:) ! Data field for the MOs. First three
! indices match last three indices of coord().
! The MO mos(i) should end up in grid(:,:,:,dst(i))
type(gam_structure), intent(inout), target, optional &
:: structure ! Context to use
real(rk), intent(in), optional :: rot(:,:) ! Rotation matrix to apply to the structure. It is
! perfectly OK to use different rotation matrix
! in different calls using the same context.
! If the rotation matrix is not given for a call where
! file= is specified, it will be reset to unit matrix.
! If the rotation matrix is not given for a pure-evaluation
! call (file= not specified), rotation matrix from the
! previous call will be used.
real(rk), intent(in), optional :: dx(:) ! Grid spacing, used for oversampling. Supplying all-negative
! input, or omitting this argument, suppresses oversampling.
integer(ik), intent(in), optional :: max_mos ! Maximum number of MOs of a single spin. Total number of MOs
! may be twice as much; the second set would corresponds to
! the opposite spin. Orbital counter in the input file for the
! second set will restart from 1.
!
type(gam_structure), pointer :: gam
integer(ik) :: ios
integer(ik) :: max_mos_ ! Either max_mos, or 0 if max_mos is not present
logical :: have_geo, have_vec, have_ecp
logical :: load_mos, evaluate_mos
logical :: eof
real(xrk) :: det_rot
!
if (present(structure)) then
gam => structure
else
gam => gam_def
end if
max_mos_ = 0
if (present(max_mos)) max_mos_ = max_mos
!
! Decide on the call semantics
!
load_mos = present(file)
evaluate_mos = present(mos)
!
! Deal with the rotation matrix
!
if (present(rot)) then
gam%rotmat = real(rot,kind=kind(gam%rotmat))
gam%rotmat_rk = real(gam%rotmat,kind=kind(gam%rotmat_rk))
else if (load_mos) then
gam%rotmat = 0.0_xrk
gam%rotmat(1,1) = 1.0_xrk
gam%rotmat(2,2) = 1.0_xrk
gam%rotmat(3,3) = 1.0_xrk
gam%rotmat_rk = real(gam%rotmat,kind=kind(gam%rotmat_rk))
end if
!
det_rot = lapack_determinant(gam%rotmat)
if (abs(abs(det_rot)-1)>spacing(100._rk)) then
write (out,"('Rotation matrix has non-unit determinant: ',g20.13)") det_rot
stop 'import_gamess%gamess_load_orbitals: Bad rotation matrix'
end if
!
! Grid spacing; if not supplied upon orbital load, disable oversampling.
!
if (present(dx)) then
gam%dx = dx
else if (load_mos) then
gam%dx = -1._rk
end if
!
! Make sure present arguments are consistent with the semantics
!
if (evaluate_mos) then
if (.not.present(dst) .or. .not.present(coord) .or. .not.present(grid)) then
write (out,"('import_gamess%gamess_load_orbitals: Arguments needed for orbital evaluation are missing')")
stop 'import_gamess%gamess_load_orbitals: Arguments needed for orbital evaluation are missing'
end if
!
! A little bit of sanity checking on dimensions
!
if ( (size(mos)/=size(dst)) .or. &
(any(dst<1)) .or. &
(any(dst>size(grid,dim=4))) .or. &
(size(coord,dim=1)/=3) .or. &
(size(coord,dim=2)/=size(grid,dim=1)) .or. &
(size(coord,dim=3)/=size(grid,dim=2)) .or. &
(size(coord,dim=4)/=size(grid,dim=3)) ) then
stop 'import_gamess%gamess_load_orbitals - Argument array sizes are not compatible'
end if
end if
!
! If no implicit load is performed, orbitals must already be present
!
if (evaluate_mos .and. .not.load_mos) then
if (.not.allocated(gam%atoms) .or. .not.allocated(gam%vectors)) then
write (out,"('import_gamess%gamess_load_orbitals: basis set and/or orbitals are not present')")
stop 'import_gamess%gamess_load_orbitals: basis set and/or orbitals are not present'
end if
end if
!
! After this point, we should not have major surprises accessing out data.
!
!
! Open GAMESS .DAT file, and parse the first $DATA and $VEC sections we see.
!
if (load_mos) then
call gamess_destroy(gam)
call TimerStart('GAMESS import - parse')
open (gam_file,file=file,action='read',position='rewind',status='old',iostat=ios)
if (ios/=0) then
write (out,"('import_gamess%gamess_load_orbitals - error ',i8,' opening ',a)") ios, trim(file)
stop 'import_gamess%gamess_load_orbitals - nofile'
end if
have_geo = .false.
have_vec = .false.
have_ecp = .false.
scan_lines: do while (.not.have_geo .or. .not.have_vec .or. .not.have_ecp)
call gam_readline(eof)
if (eof) exit scan_lines
!
if (gam_line_buf==' $DATA ') then
if (have_geo) stop 'import_gamess%gamess_load_orbitals - multiple $DATA ?!'
call parse_geometry_data(gam)
have_geo = .true.
end if
if (gam_line_buf(1:5)==' $ECP') then
if (have_ecp) stop 'import_gamess%gamess_load_orbitals - multiple $ECP ?!'
if (.not.have_geo) stop 'import_gamess%gamess_load_orbitals - not $DATA before $ECP ?!'
call parse_ecp_data(gam)
have_ecp = .true.
end if
if (gam_line_buf(1:5)==' $VEC') then
if (have_vec) stop 'import_gamess%gamess_load_orbitals - multiple $VEC ?!'
if (.not.have_geo) stop 'import_gamess%gamess_load_orbitals - no $DATA before $VEC ?!'
call parse_orbital_data(gam,max_mos_)
have_vec = .true.
end if
end do scan_lines
close(gam_file,iostat=ios)
call TimerStop('GAMESS import - parse')
end if ! load_mos
!
! MO vectors are loaded; proceed to evaluate MOs on grid.
!
if (evaluate_mos) then
call TimerStart('GAMESS import - orbitals')
call build_orbitals(gam,mos,dst,coord,grid)
call TimerStop('GAMESS import - orbitals')
end if ! evaluate_mos
!
! Special logic for the legacy semantics - MO matrix may be quite large,
! and it won't be needed after this point.
!
if (load_mos .and. evaluate_mos) then
deallocate (gam%vectors)
end if
end subroutine gamess_load_orbitals
!
! Evaluate basis functions and their gradients on grid.
! If oversampling is necessary, see gamess_oversample_functions()
!
subroutine gamess_evaluate_functions_full(xyz,basval,structure,rot)
real(rk), intent(in) :: xyz(:) ! Coordinates of the point, at which AOs must be evaluated.
real(rk), intent(out) :: basval(:,:) ! Basis functions and gradients at a grid point. The first
! index is (function,d/dx,d/dy,d/dz) or just (function)
type(gam_structure), intent(inout), target, optional &
:: structure ! Context to use
real(rk), intent(in) :: rot(:,:) ! Rotation matrix to apply to the structure. It is
! perfectly OK to use different rotation matrix
! in different calls using the same context.
! If the rotation matrix is not given, rotation matrix
! from the previous call will be used.
!
type(gam_structure), pointer :: gam
real(xrk) :: det_rot
!
call TimerStart('GAMESS Evaluate functions')
!
if (present(structure)) then
gam => structure
else
gam => gam_def
end if
!
! Deal with the rotation matrix
!
gam%rotmat = real(rot,kind=kind(gam%rotmat))
det_rot = lapack_determinant(gam%rotmat)
if (abs(det_rot-1)>spacing(100._rk)) then
write (out,"('Rotation matrix has non-unit determinant: ',g20.13)") det_rot
stop 'import_gamess%gamess_evaluate_functions: Bad rotation matrix'
end if
gam%rotmat_rk = real(gam%rotmat,kind=kind(gam%rotmat_rk))
!
if (size(xyz)/=3 .or. all(size(basval,dim=1)/=(/1,4/)) .or. size(basval,dim=2)/=gam%nbasis) then
stop 'import_gamess%gamess_evaluate_functions - Argument array sizes are not valid'
end if
!
if (.not.allocated(gam%atoms)) then
write (out,"('import_gamess%gamess_evaluate_functions: basis set is not present')")
stop 'import_gamess%gamess_evaluate_functions: basis set is not present'
end if
!
! After this point, we should not have major surprises accessing out data.
!
if (size(basval,dim=1)==1) then
call evaluate_basis_functions(gam,xyz,basval(1,:))
else
call evaluate_basis_functions_and_gradients(gam,xyz,basval)
end if
call TimerStop('GAMESS Evaluate functions')
!
end subroutine gamess_evaluate_functions_full
!
subroutine gamess_evaluate_functions_simple(xyz,basval,structure)
real(rk), intent(in) :: xyz(:) ! Coordinates of the point, at which AOs must be evaluated.
real(rk), intent(out) :: basval(:,:) ! Basis functions and gradients at a grid point. The first
! index is (function,d/dx,d/dy,d/dz) or just (function)
type(gam_structure), intent(in), target, optional &
:: structure ! Context to use
!
type(gam_structure), pointer :: gam
!
! call TimerStart('GAMESS Evaluate functions (simple)')
!
if (present(structure)) then
gam => structure
else
gam => gam_def
end if
!
if (size(xyz)/=3 .or. all(size(basval,dim=1)/=(/1,4/)) .or. size(basval,dim=2)/=gam%nbasis) then
stop 'import_gamess%gamess_evaluate_functions_simple - Argument array sizes are not valid'
end if
!
if (.not.allocated(gam%atoms)) then
write (out,"('import_gamess%gamess_evaluate_functions_simple: basis set is not present')")
stop 'import_gamess%gamess_evaluate_functions_simple: basis set is not present'
end if
!
! After this point, we should not have major surprises accessing out data.
!
if (size(basval,dim=1)==1) then
call evaluate_basis_functions(gam,xyz,basval(1,:))
else
call evaluate_basis_functions_and_gradients(gam,xyz,basval)
end if
! call TimerStop('GAMESS Evaluate functions (simple)')
!
end subroutine gamess_evaluate_functions_simple
!
! Evaluate basis functuions, overssampling them if necessary.
! Note that the oversampling procedure is NOT a full equivalent of the build_orbitals(),
! because it operates on the AOs and the MOs.
!
subroutine gamess_oversample_functions(gam,xyz,basval)
type(gam_structure), intent(in) :: gam ! Context to use
real(rk), intent(in) :: xyz(:) ! Cartesian coordinates of the point, at which AOs
! must be evaluated.
real(rk), intent(out) :: basval(:) ! Data field for the AOs.
!
real(rk) :: r_cut ! Max distance from the nearest nucleus to trigger oversampling
real(rk) :: rot_xyz(3) ! Rotated coordinates in the molecular frame
real(rk) :: r12 ! Distance to the nearest nucleus
integer(ik) :: num_special_max ! We don't need this value; it will be ignored
! Quantities below are used only if oversampling turned out to be necessary
integer(ik) :: ord(3) ! Required oversampling order along each Cartesian direction
real(rk) :: cpt(3)
real(rk), allocatable :: ptx(:), wgx(:) ! Integration abscissas and weights, X direction
real(rk), allocatable :: pty(:), wgy(:) ! Integration abscissas and weights, Y direction
real(rk), allocatable :: ptz(:), wgz(:) ! Integration abscissas and weights, Z direction
real(rk), allocatable :: pv(:) ! MO value at sub-sample grid point
real(rk), allocatable :: spv(:), s2pv(:) ! Accumulated average AO value and its square
real(rk) :: wgt, swgt
integer(ik) :: alloc
integer(ik) :: ix, iy, iz, iat
!
if (size(xyz)/=3 .or. size(basval)/=gam%nbasis) then
write (out,"('gamess_oversample_functions: Bad dimensions: size(xyz) = ',i0,' size(basval) = ',i0,' nbasis = ',i0)") &
size(xyz), size(basval), gam%nbasis
stop 'import_gamess%gamess_oversample_functions - dimensions do not match'
end if
!
! We don't want to oversample everywhere - it could get rather expensive!
! If we are sufficiently far from a nucleus, let's just return the value at the centre.
!
call prepare_for_oversampling(gam,num_special_max,r_cut)
rot_xyz = matmul(transpose(real(gam%rotmat,kind=kind(xyz))),xyz)
r12 = huge(1._rk)
get_r12: do iat=1,gam%natoms
r12 = min(r12,sum((rot_xyz(:)-gam%atoms(iat)%xyz_rk(:)/abohr)**2))
end do get_r12
r12 = sqrt(r12)
!
if (r12>r_cut) then
call evaluate_basis_functions(gam,xyz,basval)
else
call get_oversampling_order(ord,gam,xyz)
allocate (ptx(ord(1)),wgx(ord(1)), pty(ord(2)),wgy(ord(2)), ptz(ord(3)),wgz(ord(3)), &
pv(gam%nbasis), spv(gam%nbasis), s2pv(gam%nbasis),stat=alloc)
if (alloc/=0) then
write (out,"('gamess_oversample_functions: Allocation failed. code = ',i8,' ord = ',3i8,' nbasis = ',i8)") &
alloc, ord, gam%nbasis
stop 'import_gamess%gamess_oversample_functions - no memory'
end if
call MathGetQuadrature('Legendre',order=ord(1),x=ptx,w=wgx)
call MathGetQuadrature('Legendre',order=ord(2),x=pty,w=wgy)
call MathGetQuadrature('Legendre',order=ord(3),x=ptz,w=wgz)
!
swgt = 0 ; spv = 0 ; s2pv = 0
samples_z: do iz=1,ord(3)
cpt(3) = xyz(3) + 0.5_rk*gam%dx(3) * ptz(iz)
samples_y: do iy=1,ord(2)
cpt(2) = xyz(2) + 0.5_rk*gam%dx(2) * pty(iy)
samples_x: do ix=1,ord(1)
cpt(1) = xyz(1) + 0.5_rk*gam%dx(1) * ptx(ix)
wgt = wgx(ix) * wgy(iy) * wgz(iz)
call evaluate_basis_functions(gam,cpt,pv)
spv = spv + wgt * pv
s2pv = s2pv + wgt * pv**2
swgt = swgt + wgt
end do samples_x
end do samples_y
end do samples_z
!
spv = spv / swgt
s2pv = sqrt(s2pv / swgt)
basval = sign(s2pv,spv)
!
deallocate (ptx,wgx,pty,wgy,ptz,wgz,pv,spv,s2pv)
end if
end subroutine gamess_oversample_functions
!
! Load the natural orbital occupation numbers
!
subroutine gamess_load_natocc(file,natural_occ,natural_count)
character(len=*), intent(in) :: file ! Name of the GAMESS checkpoint file
real(rk), intent(out) :: natural_occ(:) ! Natural orbital occupation numbers
integer(ik), intent(out) :: natural_count ! Number of natural orbitals in the density matrix
!
integer(ik) :: ios
logical :: have_natocc
!
! Open GAMESS .DAT file, and parse the first $OCCNO section we see
!
call TimerStart('GAMESS import natocc')
open (gam_file,file=file,action='read',position='rewind',status='old',iostat=ios)
if (ios/=0) then
write (out,"('import_gamess%gamess_load_orbitals - error ',i8,' opening ',a)") ios, trim(file)
stop 'import_gamess%gamess_load_natocc - nofile'
end if
have_natocc = .false.
scan_lines: do while (.not.have_natocc)
call gam_readline
if (gam_line_buf==' $OCCNO ') then
call parse_natocc_data(natural_occ,natural_count)
have_natocc = .true.
end if
end do scan_lines
close(gam_file,iostat=ios)
call TimerStop('GAMESS import natocc')
end subroutine gamess_load_natocc
!
! Load the transition density decomposition singular values
!
subroutine gamess_load_rdmsv(file,rdm_sv,rdm_count)
character(len=*), intent(in) :: file ! Name of the GAMESS checkpoint file
real(rk), intent(out) :: rdm_sv(:) ! Density matrix singular values
integer(ik), intent(out) :: rdm_count ! Number of rdm orbitals
!
integer(ik) :: ios
logical :: have_rdm
!
! Open GAMESS .DAT file, and parse the first $RDMPCE section we see
!
call TimerStart('GAMESS import rdmsv')
open (gam_file,file=file,action='read',position='rewind',status='old',iostat=ios)
if (ios/=0) then
write (out,"('import_gamess%gamess_load_orbitals - error ',i8,' opening ',a)") ios, trim(file)
stop 'import_gamess%gamess_load_rdm_sv - nofile'
end if
have_rdm = .false.
scan_lines: do while (.not.have_rdm)
call gam_readline
if (gam_line_buf==' $RDMPCE') then
call parse_rdmsv_data(rdm_sv,rdm_count)
have_rdm = .true.
end if
end do scan_lines
close(gam_file,iostat=ios)
call TimerStop('GAMESS import rdmsv')
end subroutine gamess_load_rdmsv
!
! Evaluate 1-electron integrals. Not much is implemented at the moment ...
! Numerical efficiency is NOT an issue - in every conceivable usage case,
! the overall cost will be dominated by something else.
!
! Note that complex operators are handled in gamess_1e_integrals_complex below
!
subroutine gamess_1e_integrals_real(what,v,bra,ket,op_xyz,op_param,op_index)
character(len=*), intent(in) :: what ! What to evaluate
real(rk), intent(out) :: v(:,:) ! Buffer for the results
type(gam_structure), intent(in), target, optional :: bra
type(gam_structure), intent(in), target, optional :: ket
real(rk), intent(in), optional :: op_xyz(:) ! Position of the operator (3c integrals)
real(rk), intent(in), optional :: op_param(:) ! Additional parameters for the operator
! Required size and allowed values depend on
! the operator; see below.
integer(ik), intent(in), optional :: op_index(:) ! Additional (integer) parameters for the operator.
!
include 'import_gamess_1e_integrals_common.f90'
end subroutine gamess_1e_integrals_real
!
subroutine gamess_1e_integrals_quad(what,v,bra,ket,op_xyz,op_param,op_index)
character(len=*), intent(in) :: what ! What to evaluate
real(xrk), intent(out) :: v(:,:) ! Buffer for the results
type(gam_structure), intent(in), target, optional :: bra
type(gam_structure), intent(in), target, optional :: ket
real(xrk), intent(in), optional :: op_xyz(:) ! Position of the operator (3c integrals)
real(xrk), intent(in), optional :: op_param(:) ! Additional parameters for the operator
! Required size and allowed values depend on
! the operator; see below.
integer(ik), intent(in), optional :: op_index(:) ! Additional (integer) parameters for the operator.
!
include 'import_gamess_1e_integrals_common.f90'
end subroutine gamess_1e_integrals_quad
!
! 1-electron integrals for complex operators.
!
subroutine gamess_1e_integrals_complex(what,v,bra,ket,op_xyz,op_param,op_index)
character(len=*), intent(in) :: what ! What to evaluate
complex(rk), intent(out) :: v(:,:) ! Buffer for the results
type(gam_structure), intent(in), target, optional :: bra
type(gam_structure), intent(in), target, optional :: ket
real(rk), intent(in), optional :: op_xyz(:) ! Position of the operator (3c integrals)
real(rk), intent(in), optional :: op_param(:) ! Additional parameters for the operator
! Required size and allowed values depend on
! the operator; see below.
integer(ik), intent(in), optional :: op_index(:) ! Additional (integer) parameters for the operator.
!
character(len=40) :: lwhat ! Local copy of "what"
type(gam_structure), pointer :: l, r
type(gam_operator_data) :: op_data ! Operator data
!
call TimerStart('GAMESS '//trim(what))
l => gam_def ; r => gam_def
if (present(bra)) l => bra
if (present(ket)) r => ket
!
! What is it we want to evaluate?
!
lwhat = what
select case (lwhat(1:3))
case ('AO ')
if ( (size(v,dim=1)<l%nbasis) .or. (size(v,dim=2)<r%nbasis) ) then
write (out,"('import_gamess%gamess_1e_integrals_complex: Output buffer is not large enough')")
stop 'import_gamess%gamess_1e_integrals_complex: Output buffer is not large enough'
end if
end select
!
! 3-centre integrals need position of the operator; make sure it was supplied
! Operator-specific parameters will be handled later.
!
if (lwhat(4:5)=='3C') then
if (.not.present(op_xyz)) then
write (out,"('import_gamess%gamess_1e_integrals_complex: operator location missing for ',a)") trim(lwhat)
stop 'import_gamess%gamess_1e_integrals_complex: Missing operator position'
end if
op_data%op_xyz = real(op_xyz,kind=kind(op_data%op_xyz))
end if
!
op_data%op_name = lwhat(4:)
!
select case (lwhat)
case default
write (out,"('import_gamess%gamess_1e_integrals_complex: ',a,' are not implemented')") trim(lwhat)
stop 'import_gamess%gamess_1e_integrals_complex - integral not implemented'
!
! First the integrals where no additional parmaters are needed:
!
! 'AO PLANEWAVE' - Complex planewave operator [Exp(I K R)]; needs K vector
! 'AO R PLANEWAVE' - Component of R (X, Y, or Z), times the planewave.
! Neeeds K vector and the X/Y/Z index.
!
case ('AO PLANEWAVE')
if (.not.present(op_param)) then
write (out,"('import_gamess%gamess_1e_integrals_complex: exponent missing for ',a)") trim(lwhat)
stop 'import_gamess%gamess_1e_integrals_complex: Missing operator parameter'
end if
if (size(op_param)/=3) then
write (out,"('import_gamess%gamess_1e_integrals_complex: wrong number of parameters for ',a)") trim(lwhat)
stop 'import_gamess%gamess_1e_integrals_complex: Wrong parameter count'
end if
op_data%kvec = real(op_param(1:3),kind=kind(op_data%kvec))
call os_1e_matrix(op_data,l,r,v)
case ('AO R PLANEWAVE')
if (.not.present(op_param)) then
write (out,"('import_gamess%gamess_1e_integrals_complex: exponent missing for ',a)") trim(lwhat)
stop 'import_gamess%gamess_1e_integrals_complex: Missing operator parameter'
end if
if (size(op_param)/=3) then
write (out,"('import_gamess%gamess_1e_integrals_complex: wrong number of parameters for ',a)") trim(lwhat)
stop 'import_gamess%gamess_1e_integrals_complex: Wrong parameter count'
end if
op_data%kvec = real(op_param(1:3),kind=kind(op_data%kvec))
if (.not.present(op_index)) then
write (out,"('import_gamess%gamess_1e_integrals_complex: index missing for ',a)") trim(lwhat)
stop 'import_gamess%gamess_1e_integrals_complex: Missing operator parameter'
end if
if (size(op_index)/=1) then
write (out,"('import_gamess%gamess_1e_integrals_complex: wrong number of parameters for ',a)") trim(lwhat)
stop 'import_gamess%gamess_1e_integrals_complex: Wrong parameter count'
end if
op_data%op_i = op_index(1)
call os_1e_matrix(op_data,l,r,v)
end select
call TimerStop('GAMESS '//trim(what))
end subroutine gamess_1e_integrals_complex
!
! Print 1-electron integrals in a somewhat reasonable form
!
subroutine gamess_print_1e_integrals_real(v,bra,ket,symmetry,heading)
real(rk), intent(in) :: v(:,:) ! Integrals to print
type(gam_structure), intent(in), target, optional :: bra
type(gam_structure), intent(in), target, optional :: ket
character(len=*), intent(in), optional :: symmetry ! Expected symmetry of matrix elements, either
! 'HERMITIAN' (default), 'ANTI-HERMITIAN', or 'NONE'
character(len=*), intent(in), optional :: heading ! Headings to include, either 'BOTH', 'LEFT', or 'NONE'
! (the latter options are useful for MO printing)
!
include 'import_gamess_print_1e_integrals_common.f90'
end subroutine gamess_print_1e_integrals_real
!
! Print 1-electron integrals in a somewhat reasonable form (quad precision)
!
subroutine gamess_print_1e_integrals_quad(v,bra,ket,symmetry,heading)
real(xrk), intent(in) :: v(:,:) ! Integrals to print
type(gam_structure), intent(in), target, optional :: bra
type(gam_structure), intent(in), target, optional :: ket
character(len=*), intent(in), optional :: symmetry ! Expected symmetry of matrix elements, either
! 'HERMITIAN' (default), 'ANTI-HERMITIAN', or 'NONE'
character(len=*), intent(in), optional :: heading ! Headings to include, either 'BOTH', 'LEFT', or 'NONE'
! (the latter options are useful for MO printing)
!
include 'import_gamess_print_1e_integrals_common.f90'
end subroutine gamess_print_1e_integrals_quad
!
subroutine gamess_print_1e_integrals_complex(v,bra,ket,symmetry,heading)
complex(rk), intent(in) :: v(:,:) ! Integrals to print
type(gam_structure), intent(in), target, optional :: bra
type(gam_structure), intent(in), target, optional :: ket
character(len=*), intent(in), optional :: symmetry ! Expected symmetry of matrix elements, either
! 'HERMITIAN' (default), 'ANTI-HERMITIAN', 'SYMMETRIC', or 'NONE'
character(len=*), intent(in), optional :: heading ! Headings to include, either 'BOTH', 'LEFT', or 'NONE'
! (the latter options are useful for MO printing)
!
type(gam_structure), pointer :: l, r
integer(ik) :: ir, ic, ic_p1, ic_pn
integer(ik), parameter :: cpp = 5 ! Columns per page
logical :: usegfmt
character(len=20) :: sym, head
character(len=20) :: lab
!
call TimerStart('GAMESS print 1e')
l => gam_def ; r => gam_def
if (present(bra)) l => bra
if (present(ket)) r => ket
!
! Symmetry check
!
sym = 'HERMITIAN'
if (present(symmetry)) sym = symmetry
if (size(v,dim=1)==size(v,dim=2)) then
write (out,"('Expected integral symmetry: ',a)") trim(sym)
select case (sym)
case default
write (out,"('Symmetry code ',a,' is not recognized')") trim(sym)
stop 'import_gamess%gamess_print_1e_integrals_real - bad argument'
case ('NONE')
case ('HERMITIAN')
write (out,"('Maximum deviation from symmetry = ',g12.5)") maxval(abs(v-conjg(transpose(v))))
write (out,"(' Maximum deviation at index = ',2i6 )") maxloc(abs(v-conjg(transpose(v))))
case ('ANTI-HERMITIAN')
write (out,"('Maximum deviation from symmetry = ',g12.5)") maxval(abs(v+conjg(transpose(v))))
write (out,"(' Maximum deviation at index = ',2i6 )") maxloc(abs(v+conjg(transpose(v))))
case ('SYMMETRIC')
write (out,"('Maximum deviation from symmetry = ',g12.5)") maxval(abs(v-transpose(v)))
write (out,"(' Maximum deviation at index = ',2i6 )") maxloc(abs(v-transpose(v)))
end select
write (out,"(' Largest matrix element = ',g12.5)") maxval(abs(v))
write (out,"(' Largest matrix element is at index = ',2i6)") maxloc(abs(v))
end if
!
head = 'BOTH'
if (present(heading)) head = heading
usegfmt = any(abs(v)>=1e5) .or. all(abs(v)<=1e-3)
print_pages: do ic_p1=1,r%nbasis,cpp
ic_pn = min(r%nbasis,ic_p1+cpp-1)
write (out,"(t19,5(2x,i12,1x,12x))") (ic, ic=ic_p1,ic_pn)
if (head=='BOTH') write (out,"(t25,5(2x,a13,1x,11x))") (r%bas_labels(ic), ic=ic_p1,ic_pn)
print_rows: do ir=1,l%nbasis
lab = ' '
if (head/='NONE') lab = l%bas_labels(ir)
if (usegfmt) then
write (out,"(1x,i5,1x,a13,1x,5(2x,g12.5,1x,g12.5))") ir, lab, v(ir,ic_p1:ic_pn)
else
write (out,"(1x,i5,1x,a13,1x,5(2x,f12.5,1x,g12.5))") ir, lab, v(ir,ic_p1:ic_pn)
endif
end do print_rows
write (out,"()")
end do print_pages
!
call TimerStop('GAMESS print 1e')
end subroutine gamess_print_1e_integrals_complex
!
! Return key information on 2e integrals and integral batches
!
subroutine gamess_2e_info(what,gam,op_iparam)
character(len=*), intent(in) :: what ! What to do, see below
type(gam_structure), intent(inout), target, optional :: gam ! Context to work on
integer(ik), intent(inout), optional :: op_iparam(:) ! Additional parameters
!
type(gam_structure), pointer :: g
integer(ik) :: batch
!
call TimerStart('GAMESS 2e '//trim(what))
g => gam_def
if (present(gam)) g => gam
!
select case (what)
case default
write (out,"('import_gamess%gamess_2e_info: ',a,' are not implemented')") trim(what)
stop 'import_gamess%gamess_2e_info - operation not implemented'
case ('GET INFO')
call need_iparam(3_ik)
op_iparam(1) = g%nbatches ! Output: number of batches
op_iparam(2) = g%max_batch ! Output: number of orbitals in the largest batch
op_iparam(3) = g%nbasis ! Output: total number of atomic orbitals
case ('GET BATCH INFO')
call need_iparam(4_ik)
batch = op_iparam(1) ! Input: desired batch
if (batch<1 .or. batch>g%nbatches) then
write (out,"('import_gamess%gamess_2e_info: batch ',i10,' is outside of the valid range')") batch
stop 'import_gamess%gamess_2e_info - bad batch'
end if
op_iparam(2) = g%batches(4,batch) ! Output: size of this batch
op_iparam(3) = g%batches(3,batch) ! Output: first orbital of this batch
op_iparam(4) = g%batches(1,batch) ! Output: atom this batch belongs to
end select
!
call TimerStop('GAMESS 2e '//trim(what))
!
contains
!
subroutine need_iparam(count)
integer(ik), intent(in) :: count ! Required number of entries in op_iparam
!
if (.not.present(op_iparam)) stop 'import_gamess%gamess_2e_info - required parameter op_iparam missing'
if (size(op_iparam)<count) stop 'import_gamess%gamess_2e_info - op_iparam is too small'
end subroutine need_iparam
end subroutine gamess_2e_info
!
! Evaluate a batch of 2e integrals. Due to the number of 2e integrals in
! a typical calculation, it is impractical to evaluate them all in-memory,
! as we do for the 1e integrals above.
!
subroutine gamess_2e_integrals_real(what,v,batch,a,b,c,d,op_param,accuracy)
character(len=*), intent(in) :: what ! What to evaluate
real(rk), intent(out) :: v(:,:,:,:) ! Buffer for the results
! Integrals are in the charge-cloud order: that is,
! first two indices correspond to the first variable;
! the last two indices correspond to the second variable
integer(ik), intent(in) :: batch(:) ! Batch indices
type(gam_structure), intent(in), target, optional :: a, b, c, d ! Contexts to work on. In principle, we can handle
! the situation with every index being on a separate
! molecule - not that this is hard for an AO integral
real(rk), intent(in), optional :: op_param(:) ! Additional parameters for the operator
real(rk), intent(in), optional :: accuracy ! Desired accuracy of the integrals; contributions
! smaller than this can be neglected.
! The default is to be as accurate as possible.
! Negative accuracy will restore the default behaviour
!
include 'import_gamess_gamess_2e_integrals_common.f90'
end subroutine gamess_2e_integrals_real
!
subroutine gamess_2e_integrals_quad(what,v,batch,a,b,c,d,op_param,accuracy)
character(len=*), intent(in) :: what ! What to evaluate
real(xrk), intent(out) :: v(:,:,:,:) ! Buffer for the results
! Integrals are in the charge-cloud order: that is,
! first two indices correspond to the first variable;
! the last two indices correspond to the second variable
integer(ik), intent(in) :: batch(:) ! Batch indices