Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of spin orbit integrals over Cartesian Gaussian-type functions
10 : !> \par Literature
11 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
12 : !> \par History
13 : !> none
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 : !> - dac : Distance between the atomic centers a and c.
20 : !> - dbc : Distance between the atomic centers b and c.
21 : !> - l{a,b,c} : Angular momentum quantum number of shell a, b or c.
22 : !> - l{a,b,c}_max: Maximum angular momentum quantum number of shell a, b or c.
23 : !> - l{a,b,c}_min: Minimum angular momentum quantum number of shell a, b or c.
24 : !> - ncoset : Number of orbitals in a Cartesian orbital set.
25 : !> - npgf{a,b} : Degree of contraction of shell a or b.
26 : !> - rab : Distance vector between the atomic centers a and b.
27 : !> - rab2 : Square of the distance between the atomic centers a and b.
28 : !> - rac : Distance vector between the atomic centers a and c.
29 : !> - rac2 : Square of the distance between the atomic centers a and c.
30 : !> - rbc : Distance vector between the atomic centers b and c.
31 : !> - rbc2 : Square of the distance between the atomic centers b and c.
32 : !> - rpgf{a,b,c} : Radius of the primitive Gaussian-type function a, b or c.
33 : !> - zet{a,b,c} : Exponents of the Gaussian-type functions a, b or c.
34 : !> - zetp : Reciprocal of the sum of the exponents of orbital a and b.
35 : !> \author VW
36 : ! **************************************************************************************************
37 : MODULE ai_spin_orbit
38 : USE ai_cartesian_operators, ONLY: apply_cartesian_operator_pair,&
39 : cartesian_ladder_map
40 : USE ai_os_rr, ONLY: os_rr_coul
41 : USE kinds, ONLY: dp
42 : USE mathconstants, ONLY: pi
43 : USE orbital_pointers, ONLY: ncoset
44 : #include "../base/base_uses.f90"
45 :
46 : IMPLICIT NONE
47 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_spin_orbit'
48 : PRIVATE
49 :
50 : ! *** Public subroutines ***
51 :
52 : PUBLIC :: pso
53 :
54 : CONTAINS
55 :
56 : ! **************************************************************************************************
57 : !> \brief Calculation of the primitive paramagnetic spin orbit integrals over
58 : !> Cartesian Gaussian-type functions.
59 : !> \param la_max ...
60 : !> \param la_min ...
61 : !> \param npgfa ...
62 : !> \param rpgfa ...
63 : !> \param zeta ...
64 : !> \param lb_max ...
65 : !> \param lb_min ...
66 : !> \param npgfb ...
67 : !> \param rpgfb ...
68 : !> \param zetb ...
69 : !> \param rac ...
70 : !> \param rbc ...
71 : !> \param rab ...
72 : !> \param vab ...
73 : !> \param ldrr1 ...
74 : !> \param ldrr2 ...
75 : !> \param rr ...
76 : !> \date 02.03.2009
77 : !> \author VW
78 : !> \version 1.0
79 : ! **************************************************************************************************
80 81476 : SUBROUTINE pso(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, &
81 40738 : rac, rbc, rab, vab, ldrr1, ldrr2, rr)
82 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa
83 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
84 : INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
85 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
86 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc, rab
87 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: vab
88 : INTEGER, INTENT(IN) :: ldrr1, ldrr2
89 : REAL(dp), DIMENSION(0:ldrr1-1, ldrr2, *), &
90 : INTENT(INOUT) :: rr
91 :
92 : INTEGER, DIMENSION(2, 6), PARAMETER :: &
93 : components = RESHAPE([2, 3, 3, 2, 3, 1, 1, 3, 1, 2, 2, 1], [2, 6])
94 :
95 : INTEGER :: axis, first_a, first_b, i, ipgf, j, &
96 : jpgf, na, nb, ncol, nrow
97 40738 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: map_a, map_b
98 : REAL(dp) :: dab, f0, rab2, xhi, zet, zetab
99 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: work
100 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: derivatives
101 : REAL(dp), DIMENSION(3) :: rap, rbp, rcp
102 :
103 40738 : first_a = ncoset(la_min - 1) + 1
104 40738 : first_b = ncoset(lb_min - 1) + 1
105 40738 : nrow = ncoset(la_max) - first_a + 1
106 40738 : ncol = ncoset(lb_max) - first_b + 1
107 285166 : ALLOCATE (map_a(4, nrow, 3), map_b(4, ncol, 3))
108 325904 : ALLOCATE (work(nrow, ncoset(lb_max + 1)), derivatives(nrow, ncol, 6))
109 162952 : DO axis = 1, 3
110 122214 : CALL cartesian_ladder_map(la_max, la_min, axis, map_a(:, :, axis), source_lmin=0)
111 162952 : CALL cartesian_ladder_map(lb_max, lb_min, axis, map_b(:, :, axis), source_lmin=0)
112 : END DO
113 :
114 : ! *** Calculate the distance of the centers a and c ***
115 :
116 40738 : rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
117 40738 : dab = SQRT(rab2)
118 :
119 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
120 :
121 40738 : na = 0
122 :
123 92858 : DO ipgf = 1, npgfa
124 :
125 52120 : nb = 0
126 :
127 118782 : DO jpgf = 1, npgfb
128 :
129 : ! *** Screening ***
130 :
131 66662 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
132 29344 : DO j = nb + 1, nb + ncoset(lb_max)
133 52413 : DO i = na + 1, na + ncoset(la_max)
134 23069 : vab(i, j, 1) = 0.0_dp
135 23069 : vab(i, j, 2) = 0.0_dp
136 40267 : vab(i, j, 3) = 0.0_dp
137 : END DO
138 : END DO
139 12146 : nb = nb + ncoset(lb_max)
140 12146 : CYCLE
141 : END IF
142 :
143 : ! *** Calculate some prefactors ***
144 :
145 54516 : zetab = zeta(ipgf)*zetb(jpgf)
146 54516 : zet = zeta(ipgf) + zetb(jpgf)
147 54516 : xhi = zetab/zet
148 218064 : rap = zetb(jpgf)*rab/zet
149 218064 : rbp = -zeta(ipgf)*rab/zet
150 218064 : rcp = -(zeta(ipgf)*rac + zetb(jpgf)*rbc)/zet
151 :
152 54516 : f0 = 2.0_dp*SQRT(zet/pi)*(pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
153 :
154 : ! *** Calculate the recurrence relation ***
155 :
156 54516 : CALL os_rr_coul(rap, la_max + 1, rbp, lb_max + 1, rcp, zet, ldrr1, ldrr2, rr)
157 :
158 : ! pso = (T_yz-T_zy, T_zx-T_xz, T_xy-T_yx), T_ij = D_Ai D_Bj V.
159 : CALL apply_cartesian_operator_pair(rr(0, 1:ncoset(la_max + 1), 1:ncoset(lb_max + 1)), &
160 : derivatives, map_a, map_b, components, &
161 163548 : [.TRUE., .TRUE.], [zeta(ipgf), zetb(jpgf)], [.FALSE., .TRUE.], work)
162 218064 : DO axis = 1, 3
163 : vab(na + first_a:na + ncoset(la_max), nb + first_b:nb + ncoset(lb_max), axis) = &
164 992436 : f0*(derivatives(:, :, 2*axis - 1) - derivatives(:, :, 2*axis))
165 : END DO
166 :
167 106636 : nb = nb + ncoset(lb_max)
168 :
169 : END DO
170 :
171 92858 : na = na + ncoset(la_max)
172 :
173 : END DO
174 :
175 40738 : END SUBROUTINE pso
176 :
177 : END MODULE ai_spin_orbit
|
