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 Smooth band transition between population and coherence relaxation times.
10 : ! **************************************************************************************************
11 : MODULE nlo_T1_T2_band_transition
12 : USE ieee_arithmetic, ONLY: ieee_is_finite
13 : USE kinds, ONLY: dp
14 : USE mathconstants, ONLY: gaussi,&
15 : pi
16 : USE nlo_derivatives, ONLY: derivatives_index,&
17 : nlo_derivative_count,&
18 : nlo_derivative_indices,&
19 : nlo_n_derivatives
20 : #include "../base/base_uses.f90"
21 :
22 : IMPLICIT NONE
23 :
24 : PRIVATE
25 :
26 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nlo_T1_T2_band_transition'
27 :
28 : PUBLIC :: nlo_band_transition_type, nlo_band_transition_prepare, &
29 : nlo_band_transition_correction, &
30 : nlo_band_transition_rate
31 :
32 : TYPE :: nlo_band_transition_type
33 : REAL(dp), ALLOCATABLE :: energies(:), pair_rate(:, :)
34 : REAL(dp) :: width = 0.0_dp, rate_difference = 0.0_dp
35 : COMPLEX(dp), ALLOCATABLE :: integration_points(:), integration_weights(:)
36 : END TYPE nlo_band_transition_type
37 :
38 : CONTAINS
39 :
40 : ! **************************************************************************************************
41 : !> \brief Prepare Γ_H(X)=γ₂X+(γ₁−γ₂)exp[−(ad_H/w)²]X, ad_H(X)=[H,X].
42 : !> \param energies Band energies in Hartree
43 : !> \param width Energy scale of the smooth transition in Hartree
44 : !> \param rate1 Population damping energy in Hartree
45 : !> \param rate2 Coherence damping energy in Hartree
46 : !> \param relaxation Prepared band transition state
47 : !> \note Central rates use EXP exactly. Positive grades use Cauchy's formula
48 : !> on two horizontal contours and Gauss-Hermite integration. The contour
49 : !> shifts and factors are generated from the Hermite Jacobi matrix; no
50 : !> fitted or tabulated coefficients enter the relaxation model.
51 : !> Taylor derivatives follow analytically from the resolvent equation.
52 : !> Contour points share one workspace; lower grades are reconstructed
53 : !> once per requested grade from the already solved physical density.
54 : ! **************************************************************************************************
55 5 : SUBROUTINE nlo_band_transition_prepare(energies, width, rate1, rate2, relaxation)
56 : REAL(dp), INTENT(IN) :: energies(:), width, rate1, rate2
57 : TYPE(nlo_band_transition_type), INTENT(INOUT) :: relaxation
58 :
59 : INTEGER :: iband, jband, nband
60 : REAL(dp) :: difference, gaussian
61 :
62 : IF (.NOT. ALL(ieee_is_finite(energies)) .OR. .NOT. ieee_is_finite(width) .OR. &
63 : .NOT. ieee_is_finite(rate1) .OR. .NOT. ieee_is_finite(rate2) .OR. &
64 135 : width <= 0.0_dp .OR. rate1 <= 0.0_dp .OR. rate2 < rate1) THEN
65 0 : CPABORT('Invalid NLO band-transition relaxation parameters.')
66 : END IF
67 5 : IF (ALLOCATED(relaxation%pair_rate)) DEALLOCATE (relaxation%pair_rate)
68 5 : IF (ALLOCATED(relaxation%energies)) DEALLOCATE (relaxation%energies)
69 15 : ALLOCATE (relaxation%energies(SIZE(energies)))
70 135 : relaxation%energies(:) = energies
71 5 : relaxation%width = width
72 5 : relaxation%rate_difference = rate1 - rate2
73 5 : CALL nlo_gaussian_quadrature(relaxation%integration_points, relaxation%integration_weights)
74 5 : nband = SIZE(energies)
75 5 : CPASSERT(nband >= 2)
76 20 : ALLOCATE (relaxation%pair_rate(nband, nband))
77 135 : DO jband = 1, nband
78 3515 : DO iband = 1, nband
79 3380 : difference = energies(iband) - energies(jband)
80 3380 : gaussian = 0.0_dp
81 3380 : IF (width >= ABS(difference)/26.0_dp) gaussian = EXP(-(difference/width)**2)
82 3510 : relaxation%pair_rate(iband, jband) = rate1 + (rate2 - rate1)*(1.0_dp - gaussian)
83 : END DO
84 : END DO
85 5 : END SUBROUTINE nlo_band_transition_prepare
86 :
87 : ! **************************************************************************************************
88 : !> \brief Return the central relaxation rate for one band pair.
89 : !> \param relaxation Prepared band transition state
90 : !> \param iband First band index
91 : !> \param jband Second band index
92 : !> \return Central damping rate of the selected band pair
93 : ! **************************************************************************************************
94 47320 : FUNCTION nlo_band_transition_rate(relaxation, iband, jband) RESULT(rate)
95 : TYPE(nlo_band_transition_type), INTENT(IN) :: relaxation
96 : INTEGER, INTENT(IN) :: iband, jband
97 : REAL(dp) :: rate
98 :
99 47320 : CPASSERT(ALLOCATED(relaxation%pair_rate))
100 47320 : CPASSERT(iband >= 1)
101 47320 : CPASSERT(iband <= SIZE(relaxation%pair_rate, 1))
102 47320 : CPASSERT(jband >= 1)
103 47320 : CPASSERT(jband <= SIZE(relaxation%pair_rate, 2))
104 47320 : rate = relaxation%pair_rate(iband, jband)
105 47320 : END FUNCTION nlo_band_transition_rate
106 :
107 : ! **************************************************************************************************
108 : !> \brief Sum relaxation corrections for one complete Cartesian Taylor grade.
109 : !> \param relaxation Prepared band transition state
110 : !> \param hamiltonian Hamiltonian Taylor coefficients in the band basis
111 : !> \param density Physical Taylor coefficients; only lower grades are read
112 : !> \param grade Requested total Cartesian derivative order, between one and four
113 : !> \param correction Corrections of this grade in nlo_derivative_indices order
114 : !> \note For Y_l=(ad_H−w z_l)⁻¹ρ and R_l,α=−Σ_{0<β≤α}[H_β,Y_l,α−β],
115 : !> Y_l,α=(ρ_α+R_l,α)/(Δε−w z_l). Only grades below the requested grade
116 : !> are needed in R_l,α. Reconstruct them for one contour point, accumulate
117 : !> all corrections of the requested grade, then reuse the workspace.
118 : !> No band-sized array carries a contour-point dimension. At grade four,
119 : !> 20 lower-grade matrices replace 35 matrices for each of 96 points.
120 : ! **************************************************************************************************
121 10 : SUBROUTINE nlo_band_transition_correction(relaxation, hamiltonian, density, grade, correction)
122 : TYPE(nlo_band_transition_type), INTENT(IN) :: relaxation
123 : COMPLEX(dp), INTENT(IN) :: hamiltonian(:, :, :), density(:, :, :)
124 : INTEGER, INTENT(IN) :: grade
125 : COMPLEX(dp), INTENT(OUT) :: correction(:, :, :)
126 :
127 : COMPLEX(dp) :: weight
128 10 : COMPLEX(dp), ALLOCATABLE :: inverse(:, :), lower_density(:, :, :), &
129 10 : matrix_term(:, :), remainder(:, :)
130 : INTEGER :: alpha(3, nlo_n_derivatives), ia, iband, ibeta, igamma, ipoint, jband, nband, &
131 : ncoeff, nlower, partner(nlo_n_derivatives, nlo_n_derivatives)
132 :
133 10 : CPASSERT(grade >= 1 .AND. grade <= 4)
134 10 : CPASSERT(ALLOCATED(relaxation%pair_rate))
135 10 : nband = SIZE(relaxation%pair_rate, 1)
136 10 : nlower = nlo_derivative_count(grade - 1)
137 10 : ncoeff = nlo_derivative_count(grade)
138 10 : CPASSERT(SIZE(hamiltonian, 1) == nband)
139 10 : CPASSERT(SIZE(hamiltonian, 2) == nband)
140 10 : CPASSERT(SIZE(hamiltonian, 3) >= ncoeff)
141 10 : CPASSERT(SIZE(density, 1) == nband)
142 10 : CPASSERT(SIZE(density, 2) == nband)
143 10 : CPASSERT(SIZE(density, 3) >= nlower)
144 10 : CPASSERT(SIZE(correction, 1) == nband)
145 10 : CPASSERT(SIZE(correction, 2) == nband)
146 10 : CPASSERT(SIZE(correction, 3) == ncoeff - nlower)
147 : ALLOCATE (inverse(nband, nband), lower_density(nband, nband, nlower), &
148 120 : remainder(nband, nband), matrix_term(nband, nband))
149 10 : CALL nlo_derivative_indices(alpha)
150 10 : partner = 0
151 121 : DO ia = 2, ncoeff
152 1255 : DO ibeta = 2, ia
153 4647 : partner(ibeta, ia) = derivatives_index(alpha, alpha(:, ia) - alpha(:, ibeta))
154 : END DO
155 : END DO
156 45705 : correction = CMPLX(0.0_dp, 0.0_dp, dp)
157 970 : DO ipoint = 1, SIZE(relaxation%integration_points)
158 25920 : DO jband = 1, nband
159 674880 : DO iband = 1, nband
160 : inverse(iband, jband) = 1.0_dp/ &
161 : (relaxation%energies(iband) - relaxation%energies(jband) - &
162 673920 : relaxation%width*relaxation%integration_points(ipoint))
163 : END DO
164 : END DO
165 960 : weight = relaxation%width*relaxation%integration_weights(ipoint)
166 12586 : DO ia = 1, ncoeff
167 11616 : remainder = CMPLX(0.0_dp, 0.0_dp, dp)
168 120480 : DO ibeta = 2, ia
169 108864 : igamma = partner(ibeta, ia)
170 108864 : IF (igamma == 0) CYCLE
171 647515296 : matrix_term(:, :) = MATMUL(hamiltonian(:, :, ibeta), lower_density(:, :, igamma))
172 24903072 : remainder(:, :) = remainder - matrix_term
173 647515296 : matrix_term(:, :) = MATMUL(lower_density(:, :, igamma), hamiltonian(:, :, ibeta))
174 24988128 : remainder(:, :) = remainder + matrix_term
175 : END DO
176 8166048 : remainder(:, :) = remainder*inverse
177 12576 : IF (ia <= nlower) THEN
178 3779328 : lower_density(:, :, ia) = density(:, :, ia)*inverse + remainder
179 : ELSE
180 4386720 : correction(:, :, ia - nlower) = correction(:, :, ia - nlower) + weight*remainder
181 : END IF
182 : END DO
183 : END DO
184 45705 : correction = CMPLX(relaxation%rate_difference, 0.0_dp, dp)*correction
185 10 : END SUBROUTINE nlo_band_transition_correction
186 :
187 : ! **************************************************************************************************
188 : !> \brief Generate integration points and weights for the Gaussian matrix function.
189 : !> \param points Complex contour integration points
190 : !> \param weights Weights for the complex contour integration
191 : !> \note f(x)=exp(-x²) equals the difference of Cauchy integrals on the
192 : !> lines Im(z)=±h. Along each line, f(x±ih)=exp(h²)exp(-x²)exp(∓2ixh).
193 : !> Eigenvalues of the Hermite Jacobi matrix are quadrature nodes; the
194 : !> squared first-row eigenvector entries determine quadrature weights.
195 : !> The numerical rule uses 48 nodes per line at height h=2.5. These
196 : !> choices resolve f and its first four derivatives to about 2e-14
197 : !> scaled absolute error in checks over −10≤x≤10; they are not model
198 : !> parameters. Points and weights belong to the calling transition state.
199 : ! **************************************************************************************************
200 5 : SUBROUTINE nlo_gaussian_quadrature(points, weights)
201 : COMPLEX(dp), ALLOCATABLE, INTENT(OUT) :: points(:), weights(:)
202 :
203 : INTEGER, PARAMETER :: n_hermite = 48
204 : REAL(dp), PARAMETER :: height = 2.5_dp
205 :
206 : COMPLEX(dp) :: factor
207 : INTEGER :: info, ipoint
208 : REAL(dp) :: jacobi(n_hermite, n_hermite), &
209 : nodes(n_hermite), quadrature_weight, &
210 : work(3*n_hermite), x
211 :
212 5 : ALLOCATE (points(2*n_hermite), weights(2*n_hermite))
213 5 : jacobi = 0.0_dp
214 240 : DO ipoint = 1, n_hermite - 1
215 235 : jacobi(ipoint, ipoint + 1) = SQRT(REAL(ipoint, dp)/2.0_dp)
216 240 : jacobi(ipoint + 1, ipoint) = jacobi(ipoint, ipoint + 1)
217 : END DO
218 5 : CALL dsyev('V', 'U', n_hermite, jacobi, n_hermite, nodes, work, SIZE(work), info)
219 5 : IF (info /= 0) CPABORT('NLO Gaussian contour quadrature failed.')
220 245 : DO ipoint = 1, n_hermite
221 240 : x = nodes(ipoint)
222 240 : quadrature_weight = SQRT(pi)*jacobi(1, ipoint)**2
223 240 : factor = quadrature_weight*EXP(height**2)/(2.0_dp*pi*gaussi)
224 240 : points(ipoint) = CMPLX(x, -height, dp)
225 240 : points(ipoint + n_hermite) = CMPLX(x, height, dp)
226 240 : weights(ipoint) = -factor*EXP(2.0_dp*gaussi*x*height)
227 245 : weights(ipoint + n_hermite) = factor*EXP(-2.0_dp*gaussi*x*height)
228 : END DO
229 5 : END SUBROUTINE nlo_gaussian_quadrature
230 :
231 0 : END MODULE nlo_T1_T2_band_transition
|
