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 Does all kind of post scf calculations for semi-empirical
10 : !> \par History
11 : !> Started printing preliminary stuff for MO_CUBES and MO requires some
12 : !> more work to complete all other functionalities
13 : !> - Revise MO information printout (10.05.2021, MK)
14 : !> \author Teodoro Laino (07.2008)
15 : ! **************************************************************************************************
16 : MODULE qs_scf_post_se
17 :
18 : USE ai_moments, ONLY: moment
19 : USE atomic_kind_types, ONLY: atomic_kind_type,&
20 : get_atomic_kind
21 : USE basis_set_types, ONLY: gto_basis_set_type
22 : USE cell_types, ONLY: cell_type,&
23 : pbc
24 : USE cp_control_types, ONLY: dft_control_type
25 : USE cp_dbcsr_api, ONLY: dbcsr_get_block_p,&
26 : dbcsr_p_type
27 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
28 : USE cp_log_handling, ONLY: cp_get_default_logger,&
29 : cp_logger_get_default_io_unit,&
30 : cp_logger_type
31 : USE cp_output_handling, ONLY: cp_p_file,&
32 : cp_print_key_finished_output,&
33 : cp_print_key_should_output,&
34 : cp_print_key_unit_nr
35 : USE cp_result_methods, ONLY: cp_results_erase,&
36 : put_results
37 : USE cp_result_types, ONLY: cp_result_type
38 : USE eeq_method, ONLY: eeq_print
39 : USE input_section_types, ONLY: section_get_ival,&
40 : section_vals_get,&
41 : section_vals_get_subs_vals,&
42 : section_vals_type,&
43 : section_vals_val_get
44 : USE kinds, ONLY: default_string_length,&
45 : dp
46 : USE mathconstants, ONLY: twopi,&
47 : z_one,&
48 : z_zero
49 : USE message_passing, ONLY: mp_para_env_type
50 : USE moments_utils, ONLY: get_reference_point
51 : USE orbital_pointers, ONLY: coset,&
52 : ncoset
53 : USE particle_list_types, ONLY: particle_list_type
54 : USE particle_types, ONLY: particle_type
55 : USE physcon, ONLY: debye
56 : USE qs_environment_types, ONLY: get_qs_env,&
57 : qs_environment_type
58 : USE qs_kind_types, ONLY: get_qs_kind,&
59 : qs_kind_type
60 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
61 : USE qs_ks_types, ONLY: qs_ks_did_change
62 : USE qs_rho_types, ONLY: qs_rho_get,&
63 : qs_rho_type
64 : USE qs_scf_output, ONLY: qs_scf_write_mos
65 : USE qs_scf_types, ONLY: qs_scf_env_type
66 : USE qs_subsys_types, ONLY: qs_subsys_get,&
67 : qs_subsys_type
68 : USE semi_empirical_types, ONLY: get_se_param,&
69 : semi_empirical_type
70 : #include "./base/base_uses.f90"
71 :
72 : IMPLICIT NONE
73 : PRIVATE
74 :
75 : ! Global parameters
76 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_post_se'
77 : PUBLIC :: scf_post_calculation_se
78 :
79 : CONTAINS
80 :
81 : ! **************************************************************************************************
82 : !> \brief collects possible post - scf calculations and prints info / computes properties.
83 : !> specific for Semi-empirical calculations
84 : !> \param qs_env the qs_env in which the qs_env lives
85 : !> \par History
86 : !> 07.2008 created [tlaino] - Split from qs_scf_post (general)
87 : !> \author tlaino
88 : !> \note
89 : !> this function changes mo_eigenvectors and mo_eigenvalues, depending on the print keys.
90 : !> In particular, MO_CUBES causes the MOs to be rotated to make them eigenstates of the KS
91 : !> matrix, and mo_eigenvalues is updated accordingly. This can, for unconverged wavefunctions,
92 : !> change afterwards slightly the forces (hence small numerical differences between MD
93 : !> with and without the debug print level). Ideally this should not happen...
94 : ! **************************************************************************************************
95 7360 : SUBROUTINE scf_post_calculation_se(qs_env)
96 :
97 : TYPE(qs_environment_type), POINTER :: qs_env
98 :
99 : CHARACTER(len=*), PARAMETER :: routineN = 'scf_post_calculation_se'
100 :
101 : INTEGER :: handle, output_unit
102 : LOGICAL :: explicit, my_localized_wfn
103 : TYPE(cp_logger_type), POINTER :: logger
104 : TYPE(mp_para_env_type), POINTER :: para_env
105 : TYPE(particle_list_type), POINTER :: particles
106 : TYPE(qs_rho_type), POINTER :: rho
107 : TYPE(qs_subsys_type), POINTER :: subsys
108 : TYPE(section_vals_type), POINTER :: input, print_key, wfn_mix_section
109 :
110 3680 : CALL timeset(routineN, handle)
111 :
112 : ! Writes the data that is already available in qs_env
113 3680 : CALL write_available_results(qs_env)
114 :
115 3680 : my_localized_wfn = .FALSE.
116 3680 : NULLIFY (rho, subsys, particles, input, print_key, para_env)
117 :
118 3680 : logger => cp_get_default_logger()
119 3680 : output_unit = cp_logger_get_default_io_unit(logger)
120 :
121 3680 : CPASSERT(ASSOCIATED(qs_env))
122 : ! Here we start with data that needs a postprocessing...
123 : CALL get_qs_env(qs_env, &
124 : rho=rho, &
125 : input=input, &
126 : subsys=subsys, &
127 3680 : para_env=para_env)
128 3680 : CALL qs_subsys_get(subsys, particles=particles)
129 :
130 : ! Compute Atomic Charges
131 3680 : CALL qs_scf_post_charges(input, logger, qs_env, rho, para_env)
132 :
133 : ! Moments of charge distribution
134 3680 : CALL qs_scf_post_moments(input, logger, qs_env)
135 :
136 : ! MO_CUBES
137 : print_key => section_vals_get_subs_vals(section_vals=input, &
138 3680 : subsection_name="DFT%PRINT%MO_CUBES")
139 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
140 4 : CPWARN("Printing of MO cube files not implemented for Semi-Empirical method.")
141 : END IF
142 :
143 : ! STM
144 : print_key => section_vals_get_subs_vals(section_vals=input, &
145 3680 : subsection_name="DFT%PRINT%STM")
146 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
147 4 : CPWARN("STM not implemented for Semi-Empirical method.")
148 : END IF
149 :
150 : ! DFT+U
151 : print_key => section_vals_get_subs_vals(section_vals=input, &
152 3680 : subsection_name="DFT%PRINT%PLUS_U")
153 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
154 4 : CPWARN("DFT+U not available for Semi-Empirical method.")
155 : END IF
156 :
157 : ! Kinetic Energy
158 : print_key => section_vals_get_subs_vals(section_vals=input, &
159 3680 : subsection_name="DFT%PRINT%KINETIC_ENERGY")
160 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
161 4 : CPWARN("Kinetic energy not available for Semi-Empirical method.")
162 : END IF
163 :
164 : ! Wavefunction mixing
165 3680 : wfn_mix_section => section_vals_get_subs_vals(input, "DFT%PRINT%WFN_MIX")
166 3680 : CALL section_vals_get(wfn_mix_section, explicit=explicit)
167 3680 : IF (explicit .AND. .NOT. qs_env%run_rtp) THEN
168 0 : CPWARN("Wavefunction mixing not implemented for Semi-Empirical method.")
169 : END IF
170 :
171 : ! Print coherent X-ray diffraction spectrum
172 : print_key => section_vals_get_subs_vals(section_vals=input, &
173 3680 : subsection_name="DFT%PRINT%XRAY_DIFFRACTION_SPECTRUM")
174 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
175 4 : CPWARN("XRAY_DIFFRACTION_SPECTRUM not implemented for Semi-Empirical calculations!")
176 : END IF
177 :
178 : ! Calculation of Electric Field Gradients
179 : print_key => section_vals_get_subs_vals(section_vals=input, &
180 3680 : subsection_name="DFT%PRINT%ELECTRIC_FIELD_GRADIENT")
181 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
182 4 : CPWARN("ELECTRIC_FIELD_GRADIENT not implemented for Semi-Empirical calculations!")
183 : END IF
184 :
185 : ! Calculation of EPR Hyperfine Coupling Tensors
186 : print_key => section_vals_get_subs_vals(section_vals=input, &
187 3680 : subsection_name="DFT%PRINT%HYPERFINE_COUPLING_TENSOR")
188 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), &
189 : cp_p_file)) THEN
190 4 : CPWARN("HYPERFINE_COUPLING_TENSOR not implemented for Semi-Empirical calculations!")
191 : END IF
192 :
193 3680 : CALL timestop(handle)
194 :
195 3680 : END SUBROUTINE scf_post_calculation_se
196 :
197 : ! **************************************************************************************************
198 : !> \brief Computes and prints electric dipole moments
199 : !> We use the approximation for NDDO from
200 : !> Pople and Beveridge, Approximate Molecular Orbital Theory,
201 : !> Mc Graw Hill 1970
202 : !> mu = \sum_A [ Q_A * R_a + Tr(P_A*D_A) ]
203 : !> \param input ...
204 : !> \param logger ...
205 : !> \param qs_env the qs_env in which the qs_env lives
206 : ! **************************************************************************************************
207 3680 : SUBROUTINE qs_scf_post_moments(input, logger, qs_env)
208 : TYPE(section_vals_type), POINTER :: input
209 : TYPE(cp_logger_type), POINTER :: logger
210 : TYPE(qs_environment_type), POINTER :: qs_env
211 :
212 : CHARACTER(LEN=default_string_length) :: description, dipole_type
213 : COMPLEX(KIND=dp) :: dzeta, zeta
214 : COMPLEX(KIND=dp), DIMENSION(3) :: dggamma, dzphase, ggamma, zphase
215 : INTEGER :: i, iat, iatom, ikind, ix, j, nat, natom, &
216 : natorb, nkind, nspin, reference, &
217 : unit_nr
218 : LOGICAL :: do_berry, found
219 : REAL(KIND=dp) :: charge_tot, ci(3), dci(3), dipole(3), dipole_deriv(3), drcc(3), dria(3), &
220 : dtheta, gvec(3), q, rcc(3), ria(3), tcharge(2), theta, tmp(3), via(3), zeff
221 3680 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ncharge
222 3680 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mom
223 3680 : REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point
224 3680 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock
225 3680 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
226 : TYPE(cell_type), POINTER :: cell
227 : TYPE(cp_result_type), POINTER :: results
228 3680 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p
229 : TYPE(dft_control_type), POINTER :: dft_control
230 : TYPE(gto_basis_set_type), POINTER :: basis_set
231 3680 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
232 3680 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
233 : TYPE(qs_rho_type), POINTER :: rho
234 : TYPE(section_vals_type), POINTER :: print_key
235 : TYPE(semi_empirical_type), POINTER :: se_kind
236 :
237 3680 : NULLIFY (results)
238 7360 : print_key => section_vals_get_subs_vals(input, "DFT%PRINT%MOMENTS")
239 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
240 : ! Dipole Moments
241 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MOMENTS", &
242 26 : extension=".data", middle_name="se_dipole", log_filename=.FALSE.)
243 :
244 : ! Reference point
245 26 : reference = section_get_ival(print_key, keyword_name="REFERENCE")
246 26 : NULLIFY (ref_point)
247 26 : description = '[DIPOLE]'
248 26 : CALL section_vals_val_get(print_key, "REF_POINT", r_vals=ref_point)
249 26 : CALL section_vals_val_get(print_key, "PERIODIC", l_val=do_berry)
250 : CALL get_reference_point(rcc, drcc, qs_env=qs_env, reference=reference, &
251 26 : ref_point=ref_point)
252 : !
253 26 : NULLIFY (particle_set)
254 : CALL get_qs_env(qs_env=qs_env, &
255 : rho=rho, &
256 : cell=cell, &
257 : atomic_kind_set=atomic_kind_set, &
258 : natom=natom, &
259 : qs_kind_set=qs_kind_set, &
260 : particle_set=particle_set, &
261 : results=results, &
262 26 : dft_control=dft_control)
263 :
264 26 : CALL qs_rho_get(rho, rho_ao=matrix_p)
265 26 : nspin = SIZE(matrix_p)
266 26 : nkind = SIZE(atomic_kind_set)
267 : ! net charges
268 78 : ALLOCATE (ncharge(natom))
269 238 : ncharge = 0.0_dp
270 82 : DO ikind = 1, nkind
271 56 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat)
272 56 : CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind)
273 56 : CALL get_se_param(se_kind, zeff=zeff, natorb=natorb)
274 350 : DO iatom = 1, nat
275 212 : iat = atomic_kind_set(ikind)%atom_list(iatom)
276 212 : tcharge = 0.0_dp
277 424 : DO i = 1, nspin
278 : CALL dbcsr_get_block_p(matrix=matrix_p(i)%matrix, row=iat, col=iat, &
279 212 : block=pblock, found=found)
280 424 : IF (found) THEN
281 455 : DO j = 1, natorb
282 455 : tcharge(i) = tcharge(i) + pblock(j, j)
283 : END DO
284 : END IF
285 : END DO
286 692 : ncharge(iat) = zeff - SUM(tcharge)
287 : END DO
288 : END DO
289 : ! Contributions from net atomic charges
290 : ! Dipole deriv will be the derivative of the Dipole(dM/dt=\sum e_j v_j)
291 26 : dipole_deriv = 0.0_dp
292 26 : dipole = 0.0_dp
293 26 : IF (do_berry) THEN
294 4 : dipole_type = "periodic (Berry phase)"
295 16 : rcc = pbc(rcc, cell)
296 4 : charge_tot = 0._dp
297 150 : charge_tot = SUM(ncharge)
298 64 : ria = twopi*MATMUL(cell%h_inv, rcc)
299 16 : zphase = CMPLX(COS(ria), SIN(ria), dp)**charge_tot
300 :
301 64 : dria = twopi*MATMUL(cell%h_inv, drcc)
302 16 : dzphase = charge_tot*CMPLX(-SIN(ria), COS(ria), dp)**(charge_tot - 1.0_dp)*dria
303 :
304 16 : ggamma = z_one
305 4 : dggamma = z_zero
306 16 : DO ikind = 1, SIZE(atomic_kind_set)
307 12 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat)
308 162 : DO i = 1, nat
309 146 : iat = atomic_kind_set(ikind)%atom_list(i)
310 584 : ria = particle_set(iat)%r(:)
311 584 : ria = pbc(ria, cell)
312 584 : via = particle_set(iat)%v(:)
313 146 : q = ncharge(iat)
314 596 : DO j = 1, 3
315 1752 : gvec = twopi*cell%h_inv(j, :)
316 1752 : theta = SUM(ria(:)*gvec(:))
317 1752 : dtheta = SUM(via(:)*gvec(:))
318 438 : zeta = CMPLX(COS(theta), SIN(theta), KIND=dp)**(-q)
319 438 : dzeta = -q*CMPLX(-SIN(theta), COS(theta), KIND=dp)**(-q - 1.0_dp)*dtheta
320 438 : dggamma(j) = dggamma(j)*zeta + ggamma(j)*dzeta
321 584 : ggamma(j) = ggamma(j)*zeta
322 : END DO
323 : END DO
324 : END DO
325 16 : dggamma = dggamma*zphase + ggamma*dzphase
326 16 : ggamma = ggamma*zphase
327 16 : IF (ALL(REAL(ggamma, KIND=dp) /= 0.0_dp)) THEN
328 16 : tmp = AIMAG(ggamma)/REAL(ggamma, KIND=dp)
329 16 : ci = ATAN(tmp)
330 : dci = (1.0_dp/(1.0_dp + tmp**2))* &
331 : (AIMAG(dggamma)*REAL(ggamma, KIND=dp) - AIMAG(ggamma)* &
332 16 : REAL(dggamma, KIND=dp))/(REAL(ggamma, KIND=dp))**2
333 64 : dipole = MATMUL(cell%hmat, ci)/twopi
334 64 : dipole_deriv = MATMUL(cell%hmat, dci)/twopi
335 : END IF
336 : ELSE
337 22 : dipole_type = "non-periodic"
338 88 : DO i = 1, natom
339 : ! no pbc(particle_set(i)%r(:),cell) so that the total dipole is the sum of the molecular dipoles
340 264 : ria = particle_set(i)%r(:)
341 66 : q = ncharge(i)
342 264 : dipole = dipole - q*(ria - rcc)
343 286 : dipole_deriv(:) = dipole_deriv(:) - q*(particle_set(i)%v(:) - drcc)
344 : END DO
345 : END IF
346 : ! Contributions from atomic polarization
347 : ! No contribution to dipole derivatives
348 82 : DO ikind = 1, nkind
349 56 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat)
350 56 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set)
351 56 : CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind)
352 56 : CALL get_se_param(se_kind, natorb=natorb)
353 280 : ALLOCATE (mom(natorb, natorb, 3))
354 2180 : mom = 0.0_dp
355 56 : CALL atomic_moments(mom, basis_set)
356 268 : DO iatom = 1, nat
357 212 : iat = atomic_kind_set(ikind)%atom_list(iatom)
358 480 : DO i = 1, nspin
359 : CALL dbcsr_get_block_p(matrix=matrix_p(i)%matrix, row=iat, col=iat, &
360 212 : block=pblock, found=found)
361 424 : IF (found) THEN
362 139 : CPASSERT(natorb == SIZE(pblock, 1))
363 139 : ix = coset(1, 0, 0) - 1
364 1479 : dipole(1) = dipole(1) + SUM(pblock*mom(:, :, ix))
365 139 : ix = coset(0, 1, 0) - 1
366 1479 : dipole(2) = dipole(2) + SUM(pblock*mom(:, :, ix))
367 139 : ix = coset(0, 0, 1) - 1
368 1479 : dipole(3) = dipole(3) + SUM(pblock*mom(:, :, ix))
369 : END IF
370 : END DO
371 : END DO
372 194 : DEALLOCATE (mom)
373 : END DO
374 26 : CALL cp_results_erase(results=results, description=description)
375 : CALL put_results(results=results, description=description, &
376 26 : values=dipole(1:3))
377 26 : IF (unit_nr > 0) THEN
378 : WRITE (unit_nr, '(/,T2,A,T31,A50)') &
379 24 : 'SE_DIPOLE| Dipole type', ADJUSTR(TRIM(dipole_type))
380 : WRITE (unit_nr, '(T2,A,T30,3(1X,F16.8))') &
381 24 : 'SE_DIPOLE| Moment [a.u.]', dipole(1:3)
382 : WRITE (unit_nr, '(T2,A,T30,3(1X,F16.8))') &
383 96 : 'SE_DIPOLE| Moment [Debye]', dipole(1:3)*debye
384 : WRITE (unit_nr, '(T2,A,T30,3(1X,F16.8))') &
385 24 : 'SE_DIPOLE| Derivative [a.u.]', dipole_deriv(1:3)
386 : END IF
387 52 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
388 : END IF
389 :
390 7360 : END SUBROUTINE qs_scf_post_moments
391 :
392 : ! **************************************************************************************************
393 : !> \brief Computes the dipole integrals for an atom (a|x|b), a,b on atom A
394 : !> \param mom ...
395 : !> \param basis_set ...
396 : ! **************************************************************************************************
397 56 : SUBROUTINE atomic_moments(mom, basis_set)
398 : REAL(KIND=dp), DIMENSION(:, :, :) :: mom
399 : TYPE(gto_basis_set_type), POINTER :: basis_set
400 :
401 : INTEGER :: i, iset, jset, ncoa, ncob, nm, nset, &
402 : sgfa, sgfb
403 56 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgf, nsgf
404 56 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf
405 56 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
406 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
407 : REAL(KIND=dp), DIMENSION(3) :: rac, rbc
408 56 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf, sphi, zet
409 :
410 56 : rac = 0.0_dp
411 56 : rbc = 0.0_dp
412 :
413 56 : first_sgf => basis_set%first_sgf
414 56 : la_max => basis_set%lmax
415 56 : la_min => basis_set%lmin
416 56 : npgf => basis_set%npgf
417 56 : nset = basis_set%nset
418 56 : nsgf => basis_set%nsgf_set
419 56 : rpgf => basis_set%pgf_radius
420 56 : sphi => basis_set%sphi
421 56 : zet => basis_set%zet
422 :
423 56 : nm = 0
424 142 : DO iset = 1, nset
425 86 : ncoa = npgf(iset)*ncoset(la_max(iset))
426 142 : nm = MAX(nm, ncoa)
427 : END DO
428 448 : ALLOCATE (mab(nm, nm, 4), work(nm, nm))
429 :
430 142 : DO iset = 1, nset
431 86 : ncoa = npgf(iset)*ncoset(la_max(iset))
432 86 : sgfa = first_sgf(1, iset)
433 288 : DO jset = 1, nset
434 146 : ncob = npgf(jset)*ncoset(la_max(jset))
435 146 : sgfb = first_sgf(1, jset)
436 : !*** Calculate the primitive integrals ***
437 : CALL moment(la_max(iset), npgf(iset), zet(:, iset), rpgf(:, iset), la_min(iset), &
438 146 : la_max(jset), npgf(jset), zet(:, jset), rpgf(:, jset), 1, rac, rbc, mab)
439 : !*** Contraction step ***
440 670 : DO i = 1, 3
441 : CALL dgemm("N", "N", ncoa, nsgf(jset), ncob, 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
442 438 : sphi(1, sgfb), SIZE(sphi, 1), 0.0_dp, work(1, 1), SIZE(work, 1))
443 : CALL dgemm("T", "N", nsgf(iset), nsgf(jset), ncoa, 1.0_dp, sphi(1, sgfa), SIZE(sphi, 1), &
444 584 : work(1, 1), SIZE(work, 1), 1.0_dp, mom(sgfa, sgfb, i), SIZE(mom, 1))
445 : END DO
446 : END DO
447 : END DO
448 56 : DEALLOCATE (mab, work)
449 :
450 56 : END SUBROUTINE atomic_moments
451 : ! **************************************************************************************************
452 : !> \brief Computes and Prints Atomic Charges with several methods
453 : !> \param input ...
454 : !> \param logger ...
455 : !> \param qs_env the qs_env in which the qs_env lives
456 : !> \param rho ...
457 : !> \param para_env ...
458 : ! **************************************************************************************************
459 3680 : SUBROUTINE qs_scf_post_charges(input, logger, qs_env, rho, para_env)
460 : TYPE(section_vals_type), POINTER :: input
461 : TYPE(cp_logger_type), POINTER :: logger
462 : TYPE(qs_environment_type), POINTER :: qs_env
463 : TYPE(qs_rho_type), POINTER :: rho
464 : TYPE(mp_para_env_type), POINTER :: para_env
465 :
466 : CHARACTER(LEN=2) :: aname
467 : INTEGER :: i, iat, iatom, ikind, j, nat, natom, &
468 : natorb, nkind, nspin, unit_nr
469 : LOGICAL :: found
470 : REAL(KIND=dp) :: npe, zeff
471 3680 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mcharge
472 3680 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: charges
473 3680 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock
474 3680 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
475 3680 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p
476 : TYPE(dft_control_type), POINTER :: dft_control
477 3680 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
478 3680 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
479 : TYPE(section_vals_type), POINTER :: print_key
480 : TYPE(semi_empirical_type), POINTER :: se_kind
481 :
482 3680 : NULLIFY (particle_set)
483 : CALL get_qs_env(qs_env=qs_env, &
484 : atomic_kind_set=atomic_kind_set, &
485 : natom=natom, &
486 : qs_kind_set=qs_kind_set, &
487 : particle_set=particle_set, &
488 3680 : dft_control=dft_control)
489 :
490 : ! Compute the mulliken charges
491 3680 : print_key => section_vals_get_subs_vals(input, "DFT%PRINT%MULLIKEN")
492 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
493 2866 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MULLIKEN", extension=".mulliken", log_filename=.FALSE.)
494 2866 : CALL qs_rho_get(rho, rho_ao=matrix_p)
495 2866 : nspin = SIZE(matrix_p)
496 2866 : npe = REAL(para_env%num_pe, KIND=dp)
497 17196 : ALLOCATE (charges(natom, nspin), mcharge(natom))
498 17342 : charges = 0.0_dp
499 14134 : mcharge = 0.0_dp
500 : ! calculate atomic charges
501 2866 : nkind = SIZE(atomic_kind_set)
502 8810 : DO ikind = 1, nkind
503 5944 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat)
504 5944 : CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind)
505 5944 : CALL get_se_param(se_kind, zeff=zeff, natorb=natorb)
506 26022 : DO iatom = 1, nat
507 11268 : iat = atomic_kind_set(ikind)%atom_list(iatom)
508 22794 : DO i = 1, nspin
509 : CALL dbcsr_get_block_p(matrix=matrix_p(i)%matrix, row=iat, col=iat, &
510 11526 : block=pblock, found=found)
511 22794 : IF (found) THEN
512 20982 : DO j = 1, natorb
513 20982 : charges(iat, i) = charges(iat, i) + pblock(j, j)
514 : END DO
515 : END IF
516 : END DO
517 28738 : mcharge(iat) = zeff/npe - SUM(charges(iat, 1:nspin))
518 : END DO
519 : END DO
520 : !
521 2866 : CALL para_env%sum(charges)
522 2866 : CALL para_env%sum(mcharge)
523 : !
524 2866 : IF (unit_nr > 0) THEN
525 1444 : WRITE (UNIT=unit_nr, FMT="(/,/,T2,A)") "POPULATION ANALYSIS"
526 1444 : IF (nspin == 1) THEN
527 : WRITE (UNIT=unit_nr, FMT="(/,T2,A,T70,A)") &
528 1402 : " # Atom Element Kind Atomic population", " Net charge"
529 4312 : DO ikind = 1, nkind
530 2910 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat)
531 2910 : CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind, element_symbol=aname)
532 12760 : DO iatom = 1, nat
533 5538 : iat = atomic_kind_set(ikind)%atom_list(iatom)
534 : WRITE (UNIT=unit_nr, &
535 : FMT="(T2,I7,6X,A2,3X,I6,T39,F12.6,T69,F12.6)") &
536 8448 : iat, aname, ikind, charges(iat, 1), mcharge(iat)
537 : END DO
538 : END DO
539 : WRITE (UNIT=unit_nr, &
540 : FMT="(T2,A,T39,F12.6,T69,F12.6,/)") &
541 12478 : "# Total charge", SUM(charges(:, 1)), SUM(mcharge(:))
542 : ELSE
543 : WRITE (UNIT=unit_nr, FMT="(/,T2,A)") &
544 42 : "# Atom Element Kind Atomic population (alpha,beta) Net charge Spin moment"
545 126 : DO ikind = 1, nkind
546 84 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=nat)
547 84 : CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind, element_symbol=aname)
548 339 : DO iatom = 1, nat
549 129 : iat = atomic_kind_set(ikind)%atom_list(iatom)
550 : WRITE (UNIT=unit_nr, &
551 : FMT="(T2,I6,5X,A2,2X,I6,T29,4(1X,F12.6))") &
552 213 : iat, aname, ikind, charges(iat, 1:2), mcharge(iat), charges(iat, 1) - charges(iat, 2)
553 : END DO
554 : END DO
555 : WRITE (UNIT=unit_nr, &
556 : FMT="(T2,A,T29,4(1X,F12.6),/)") &
557 429 : "# Total charge and spin", SUM(charges(:, 1)), SUM(charges(:, 2)), SUM(mcharge(:))
558 : END IF
559 : END IF
560 :
561 2866 : CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MULLIKEN")
562 :
563 2866 : DEALLOCATE (charges, mcharge)
564 : END IF
565 :
566 : ! EEQ Charges
567 3680 : print_key => section_vals_get_subs_vals(input, "DFT%PRINT%EEQ_CHARGES")
568 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
569 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%EEQ_CHARGES", &
570 4 : extension=".eeq", log_filename=.FALSE.)
571 4 : CALL eeq_print(qs_env, unit_nr, print_level=1, ext=.FALSE.)
572 4 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
573 : END IF
574 :
575 : ! Compute the Lowdin charges
576 3680 : print_key => section_vals_get_subs_vals(input, "DFT%PRINT%LOWDIN")
577 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
578 4 : CPWARN("Lowdin charges not available for semi-empirical calculations!")
579 : END IF
580 :
581 : ! Hirshfeld charges
582 3680 : print_key => section_vals_get_subs_vals(input, "DFT%PRINT%HIRSHFELD")
583 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
584 2866 : CPWARN("Hirshfeld charges not available for semi-empirical calculations!")
585 : END IF
586 :
587 : ! MAO
588 3680 : print_key => section_vals_get_subs_vals(input, "DFT%PRINT%MAO_ANALYSIS")
589 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
590 4 : CPWARN("MAO analysis not available for semi-empirical calculations!")
591 : END IF
592 :
593 : ! ED
594 3680 : print_key => section_vals_get_subs_vals(input, "DFT%PRINT%ENERGY_DECOMPOSITION_ANALYSIS")
595 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
596 4 : CPWARN("ED analysis not available for semi-empirical calculations!")
597 : END IF
598 :
599 7360 : END SUBROUTINE qs_scf_post_charges
600 :
601 : ! **************************************************************************************************
602 : !> \brief Write QS results always available (if switched on through the print_keys)
603 : !> \param qs_env the qs_env in which the qs_env lives
604 : ! **************************************************************************************************
605 7360 : SUBROUTINE write_available_results(qs_env)
606 : TYPE(qs_environment_type), POINTER :: qs_env
607 :
608 : CHARACTER(len=*), PARAMETER :: routineN = 'write_available_results'
609 :
610 : INTEGER :: after, handle, ispin, iw, output_unit
611 : LOGICAL :: omit_headers
612 : TYPE(cp_logger_type), POINTER :: logger
613 3680 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv, rho_ao
614 : TYPE(dft_control_type), POINTER :: dft_control
615 : TYPE(mp_para_env_type), POINTER :: para_env
616 : TYPE(particle_list_type), POINTER :: particles
617 3680 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
618 : TYPE(qs_rho_type), POINTER :: rho
619 : TYPE(qs_scf_env_type), POINTER :: scf_env
620 : TYPE(qs_subsys_type), POINTER :: subsys
621 : TYPE(section_vals_type), POINTER :: dft_section, input
622 :
623 3680 : CALL timeset(routineN, handle)
624 3680 : NULLIFY (dft_control, particle_set, rho, ks_rmpv, dft_section, input, &
625 3680 : particles, subsys, para_env, rho_ao)
626 3680 : logger => cp_get_default_logger()
627 3680 : output_unit = cp_logger_get_default_io_unit(logger)
628 :
629 3680 : CPASSERT(ASSOCIATED(qs_env))
630 : CALL get_qs_env(qs_env, &
631 : dft_control=dft_control, &
632 : particle_set=particle_set, &
633 : rho=rho, &
634 : matrix_ks=ks_rmpv, &
635 : input=input, &
636 : subsys=subsys, &
637 : scf_env=scf_env, &
638 3680 : para_env=para_env)
639 3680 : CALL qs_subsys_get(subsys, particles=particles)
640 3680 : CALL qs_rho_get(rho, rho_ao=rho_ao)
641 :
642 : ! Print MO information if requested
643 3680 : CALL qs_scf_write_mos(qs_env, scf_env, final_mos=.TRUE.)
644 :
645 : ! Aat the end of SCF printout the projected DOS for each atomic kind
646 3680 : dft_section => section_vals_get_subs_vals(input, "DFT")
647 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%PDOS") &
648 : , cp_p_file)) THEN
649 4 : CPWARN("PDOS not implemented for Semi-Empirical calculations!")
650 : END IF
651 :
652 : ! Print the total density (electronic + core charge)
653 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
654 : "DFT%PRINT%TOT_DENSITY_CUBE"), cp_p_file)) THEN
655 4 : CPWARN("TOT_DENSITY_CUBE not implemented for Semi-Empirical calculations!")
656 : END IF
657 :
658 : ! Write cube file with electron density
659 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
660 : "DFT%PRINT%E_DENSITY_CUBE"), cp_p_file)) THEN
661 4 : CPWARN("E_DENSITY_CUBE not implemented for Semi-Empirical calculations!")
662 : END IF ! print key
663 :
664 : ! Write cube file with EFIELD
665 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
666 : "DFT%PRINT%EFIELD_CUBE"), cp_p_file)) THEN
667 4 : CPWARN("EFIELD_CUBE not implemented for Semi-Empirical calculations!")
668 : END IF ! print key
669 :
670 : ! Write cube file with ELF
671 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
672 : "DFT%PRINT%ELF_CUBE"), cp_p_file)) THEN
673 4 : CPWARN("ELF function not implemented for Semi-Empirical calculations!")
674 : END IF ! print key
675 :
676 : ! Print the hartree potential
677 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
678 : "DFT%PRINT%V_HARTREE_CUBE"), cp_p_file)) THEN
679 4 : CPWARN("V_HARTREE_CUBE not implemented for Semi-Empirical calculations!")
680 : END IF
681 :
682 : ! Print the XC potential
683 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
684 : "DFT%PRINT%V_XC_CUBE"), cp_p_file)) THEN
685 4 : CPWARN("V_XC_CUBE not available for Semi-Empirical calculations!")
686 : END IF
687 :
688 : ! Write the density matrix
689 3680 : CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
690 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
691 : "DFT%PRINT%AO_MATRICES/DENSITY"), cp_p_file)) THEN
692 : iw = cp_print_key_unit_nr(logger, input, "DFT%PRINT%AO_MATRICES/DENSITY", &
693 780 : extension=".Log")
694 780 : CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
695 780 : after = MIN(MAX(after, 1), 16)
696 1566 : DO ispin = 1, dft_control%nspins
697 : CALL cp_dbcsr_write_sparse_matrix(rho_ao(ispin)%matrix, 4, after, qs_env, &
698 1566 : para_env, output_unit=iw, omit_headers=omit_headers)
699 : END DO
700 : CALL cp_print_key_finished_output(iw, logger, input, &
701 780 : "DFT%PRINT%AO_MATRICES/DENSITY")
702 : END IF
703 :
704 : ! The Kohn-Sham matrix itself
705 3680 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
706 : "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX"), cp_p_file)) THEN
707 632 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.FALSE.)
708 632 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
709 : iw = cp_print_key_unit_nr(logger, input, "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX", &
710 632 : extension=".Log")
711 632 : CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
712 632 : after = MIN(MAX(after, 1), 16)
713 : CALL cp_dbcsr_write_sparse_matrix(ks_rmpv(1)%matrix, 4, after, qs_env, &
714 632 : para_env, output_unit=iw, omit_headers=omit_headers)
715 : CALL cp_print_key_finished_output(iw, logger, input, &
716 632 : "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX")
717 : END IF
718 :
719 3680 : CALL timestop(handle)
720 :
721 3680 : END SUBROUTINE write_available_results
722 :
723 : END MODULE qs_scf_post_se
|
