Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of the Fermi contact integrals over Cartesian
10 : !> Gaussian-type functions.
11 : !> \par Literature
12 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
13 : !> \par History
14 : !> \par Parameters
15 : !> - ax,ay,az : Angular momentum index numbers of orbital a.
16 : !> - bx,by,bz : Angular momentum index numbers of orbital b.
17 : !> - coset : Cartesian orbital set pointer.
18 : !> - dab : Distance between the atomic centers a and b.
19 : !> - l{a,b} : Angular momentum quantum number of shell a or b.
20 : !> - l{a,b}_max: Maximum angular momentum quantum number of shell a or b.
21 : !> - l{a,b}_min: Minimum angular momentum quantum number of shell a or b.
22 : !> - rab : Distance vector between the atomic centers a and b.
23 : !> - rpgf{a,b} : Radius of the primitive Gaussian-type function a or b.
24 : !> - sab : Shell set of overlap integrals.
25 : !> - zet{a,b} : Exponents of the Gaussian-type functions a or b.
26 : !> - zetp : Reciprocal of the sum of the exponents of orbital a and b.
27 : !> \author Matthias Krack (08.10.1999)
28 : ! **************************************************************************************************
29 : MODULE ai_fermi_contact
30 :
31 : USE kinds, ONLY: dp
32 : USE mathconstants, ONLY: pi
33 : USE orbital_pointers, ONLY: coset,&
34 : ncoset
35 : #include "../base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 :
39 : PRIVATE
40 :
41 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_fermi_contact'
42 :
43 : ! *** Public subroutines ***
44 : PUBLIC :: fermi_contact
45 :
46 : CONTAINS
47 :
48 : ! **************************************************************************************************
49 : !> \brief Purpose: Calculation of the two-center Fermi contact integrals
50 : !> 4/3*pi*[a|delta(r-c)|b] over Cartesian Gaussian-type functions.
51 : !> \param la_max ...
52 : !> \param la_min ...
53 : !> \param npgfa ...
54 : !> \param rpgfa ...
55 : !> \param zeta ...
56 : !> \param lb_max ...
57 : !> \param lb_min ...
58 : !> \param npgfb ...
59 : !> \param rpgfb ...
60 : !> \param zetb ...
61 : !> \param rac ...
62 : !> \param rbc ...
63 : !> \param dab ...
64 : !> \param fcab ...
65 : !> \param ldfc ...
66 : !> \date 27.02.2009
67 : !> \author VW
68 : !> \version 1.0
69 : ! **************************************************************************************************
70 9222 : SUBROUTINE fermi_contact(la_max, la_min, npgfa, rpgfa, zeta, &
71 18444 : lb_max, lb_min, npgfb, rpgfb, zetb, &
72 9222 : rac, rbc, dab, fcab, ldfc)
73 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa
74 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
75 : INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
76 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
77 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
78 : REAL(KIND=dp) :: dab
79 : INTEGER, INTENT(IN) :: ldfc
80 : REAL(KIND=dp), DIMENSION(ldfc, *), INTENT(INOUT) :: fcab
81 :
82 : INTEGER :: ax, ay, az, bx, by, bz, coa, cob, i, &
83 : ipgf, j, jpgf, la, lb, ma, mb, na, nb
84 : REAL(KIND=dp) :: dac2, dbc2, f0, fax, fay, faz, fbx, fby, &
85 : fbz
86 :
87 : ! *** Calculate some prefactors ***
88 :
89 9222 : dac2 = rac(1)**2 + rac(2)**2 + rac(3)**2
90 9222 : dbc2 = rbc(1)**2 + rbc(2)**2 + rbc(3)**2
91 :
92 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
93 :
94 9222 : na = 0
95 :
96 21050 : DO ipgf = 1, npgfa
97 :
98 11828 : nb = 0
99 :
100 27026 : DO jpgf = 1, npgfb
101 :
102 : ! *** Screening ***
103 :
104 15198 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
105 6888 : DO j = nb + 1, nb + ncoset(lb_max)
106 12801 : DO i = na + 1, na + ncoset(la_max)
107 9999 : fcab(i, j) = 0.0_dp
108 : END DO
109 : END DO
110 : nb = nb + ncoset(lb_max)
111 : CYCLE
112 : END IF
113 :
114 : ! *** Calculate some prefactors ***
115 :
116 12396 : f0 = 4.0_dp/3.0_dp*pi*EXP(-zeta(ipgf)*dac2 - zetb(jpgf)*dbc2)
117 :
118 : ! *** Calculate the primitive Fermi contact integrals ***
119 :
120 27668 : DO lb = lb_min, lb_max
121 45827 : DO bx = 0, lb
122 18159 : fbx = 1.0_dp
123 18159 : IF (bx > 0) fbx = (rbc(1))**bx
124 54477 : DO by = 0, lb - bx
125 21046 : bz = lb - bx - by
126 21046 : cob = coset(bx, by, bz)
127 21046 : mb = nb + cob
128 21046 : fby = 1.0_dp
129 21046 : IF (by > 0) fby = (rbc(2))**by
130 21046 : fbz = 1.0_dp
131 21046 : IF (bz > 0) fbz = (rbc(3))**bz
132 65949 : DO la = la_min, la_max
133 80245 : DO ax = 0, la
134 32455 : fax = fbx
135 32455 : IF (ax > 0) fax = fbx*(rac(1))**ax
136 97365 : DO ay = 0, la - ax
137 38166 : az = la - ax - ay
138 38166 : coa = coset(ax, ay, az)
139 38166 : ma = na + coa
140 38166 : fay = fby
141 38166 : IF (ay > 0) fay = fby*(rac(2))**ay
142 38166 : faz = fbz
143 38166 : IF (az > 0) faz = fbz*(rac(3))**az
144 :
145 70621 : fcab(ma, mb) = f0*fax*fay*faz
146 :
147 : END DO
148 : END DO
149 : END DO !la
150 :
151 : END DO
152 : END DO
153 : END DO !lb
154 :
155 24224 : nb = nb + ncoset(lb_max)
156 :
157 : END DO
158 :
159 21050 : na = na + ncoset(la_max)
160 :
161 : END DO
162 :
163 9222 : END SUBROUTINE fermi_contact
164 :
165 : END MODULE ai_fermi_contact
|
