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