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