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 compute mulliken charges
10 : !> we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
11 : !> \author Joost VandeVondele March 2003
12 : ! **************************************************************************************************
13 : MODULE mulliken
14 : USE atomic_charges, ONLY: print_atomic_charges
15 : USE cp_control_types, ONLY: mulliken_restraint_type
16 : USE cp_dbcsr_api, ONLY: &
17 : dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
18 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_type
19 : USE kinds, ONLY: dp
20 : USE message_passing, ONLY: mp_para_env_type
21 : USE particle_types, ONLY: particle_type
22 : USE qs_kind_types, ONLY: qs_kind_type
23 : #include "./base/base_uses.f90"
24 :
25 : IMPLICIT NONE
26 :
27 : PRIVATE
28 :
29 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mulliken'
30 :
31 : ! *** Public subroutines ***
32 :
33 : PUBLIC :: mulliken_charges, ao_charges, mulliken_restraint, compute_bond_order
34 : PUBLIC :: compute_charges, atom_trace
35 :
36 : INTERFACE mulliken_charges
37 : MODULE PROCEDURE mulliken_charges_a, mulliken_charges_b, mulliken_charges_c, &
38 : mulliken_charges_s, &
39 : mulliken_charges_akp, mulliken_charges_bkp, mulliken_charges_ckp
40 : END INTERFACE
41 :
42 : INTERFACE ao_charges
43 : MODULE PROCEDURE ao_charges_1, ao_charges_2, ao_charges_kp, ao_charges_kp_2
44 : END INTERFACE
45 :
46 : CONTAINS
47 :
48 : ! **************************************************************************************************
49 : !> \brief computes the energy and density matrix derivate of a constraint on the
50 : !> mulliken charges
51 : !>
52 : !> optional outputs:
53 : !> computes energy (added)
54 : !> contribution to KS matrix (added)
55 : !> contribution to W matrix (added)
56 : !> \param mulliken_restraint_control additional parameters needed to control the restraint
57 : !> \param para_env para_env of the matrices
58 : !> \param s_matrix ,p_matrix : containing the respective quantities
59 : !> \param p_matrix ...
60 : !> \param energy ...
61 : !> \param order_p ...
62 : !> \param ks_matrix ...
63 : !> \param w_matrix ...
64 : !> \par History
65 : !> 06.2004 created [Joost VandeVondele]
66 : !> \note
67 : !> contribution to the KS matrix is derivative wrt P
68 : !> contribution to the W matrix is derivate wrt S (sign?)
69 : !> needed for orbital and ionic forces respectively
70 : ! **************************************************************************************************
71 36 : SUBROUTINE mulliken_restraint(mulliken_restraint_control, para_env, &
72 : s_matrix, p_matrix, energy, order_p, ks_matrix, w_matrix)
73 : TYPE(mulliken_restraint_type), INTENT(IN) :: mulliken_restraint_control
74 : TYPE(mp_para_env_type), POINTER :: para_env
75 : TYPE(dbcsr_type), POINTER :: s_matrix
76 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
77 : REAL(KIND=dp), OPTIONAL :: energy, order_p
78 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
79 : POINTER :: ks_matrix, w_matrix
80 :
81 : INTEGER :: blk, iblock_col, iblock_row, ispin, &
82 : nblock, nspin
83 : LOGICAL :: found
84 : REAL(kind=dp) :: mult, my_energy, my_order_p
85 36 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges, charges_deriv, ks_block, &
86 36 : p_block, s_block, w_block
87 : TYPE(dbcsr_iterator_type) :: iter
88 :
89 : ! here we get the numbers for charges
90 :
91 36 : nspin = SIZE(p_matrix)
92 36 : CALL dbcsr_get_info(s_matrix, nblkrows_total=nblock)
93 :
94 144 : ALLOCATE (charges(nblock, nspin))
95 108 : ALLOCATE (charges_deriv(nblock, nspin))
96 36 : CALL compute_charges(p_matrix, s_matrix, charges, para_env)
97 : !
98 : ! this can be used to check the correct implementation of the derivative
99 : ! CALL rf_deriv_check(mulliken_restraint_control,charges)
100 : !
101 : CALL restraint_functional(mulliken_restraint_control, &
102 36 : charges, charges_deriv, my_energy, my_order_p)
103 :
104 36 : IF (PRESENT(order_p)) THEN
105 30 : order_p = my_order_p
106 : END IF
107 36 : IF (PRESENT(energy)) THEN
108 30 : energy = my_energy
109 : END IF
110 :
111 36 : IF (PRESENT(ks_matrix)) THEN
112 :
113 54 : DO ispin = 1, nspin
114 36 : CALL dbcsr_iterator_start(iter, s_matrix)
115 414 : DO WHILE (dbcsr_iterator_blocks_left(iter))
116 378 : CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, s_block, blk)
117 : CALL dbcsr_get_block_p(matrix=ks_matrix(ispin)%matrix, &
118 378 : row=iblock_row, col=iblock_col, BLOCK=ks_block, found=found)
119 :
120 378 : IF (.NOT. (ASSOCIATED(s_block) .AND. ASSOCIATED(ks_block))) THEN
121 0 : CPABORT("Unexpected s / ks structure")
122 : END IF
123 : mult = 0.5_dp*charges_deriv(iblock_row, ispin) + &
124 378 : 0.5_dp*charges_deriv(iblock_col, ispin)
125 :
126 51030 : ks_block = ks_block + mult*s_block
127 :
128 : END DO
129 90 : CALL dbcsr_iterator_stop(iter)
130 : END DO
131 :
132 : END IF
133 :
134 36 : IF (PRESENT(w_matrix)) THEN
135 :
136 18 : DO ispin = 1, nspin
137 12 : CALL dbcsr_iterator_start(iter, p_matrix(ispin)%matrix)
138 138 : DO WHILE (dbcsr_iterator_blocks_left(iter))
139 126 : CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, p_block, blk)
140 : CALL dbcsr_get_block_p(matrix=w_matrix(ispin)%matrix, &
141 126 : row=iblock_row, col=iblock_col, BLOCK=w_block, found=found)
142 :
143 : ! we can cycle if a block is not present
144 126 : IF (.NOT. (ASSOCIATED(w_block) .AND. ASSOCIATED(p_block))) CYCLE
145 :
146 : ! minus sign relates to convention for W
147 : mult = -0.5_dp*charges_deriv(iblock_row, ispin) &
148 126 : - 0.5_dp*charges_deriv(iblock_col, ispin)
149 :
150 17010 : w_block = w_block + mult*p_block
151 :
152 : END DO
153 30 : CALL dbcsr_iterator_stop(iter)
154 : END DO
155 :
156 : END IF
157 :
158 36 : DEALLOCATE (charges)
159 36 : DEALLOCATE (charges_deriv)
160 :
161 72 : END SUBROUTINE mulliken_restraint
162 :
163 : ! **************************************************************************************************
164 : !> \brief computes energy and derivatives given a set of charges
165 : !> this implementation uses the spin density on a number of atoms
166 : !> as a penalty function
167 : !> \param mulliken_restraint_control ...
168 : !> \param charges (nblock,nspin)
169 : !> \param charges_deriv derivate wrt the corresponding charge entry
170 : !> \param energy ...
171 : !> \param order_p ...
172 : !> \par History
173 : !> 06.2004 created [Joost VandeVondele]
174 : !> 02.2005 added more general form [Joost VandeVondele]
175 : !> \note
176 : !> should be easy to adapt for other specialized cases
177 : ! **************************************************************************************************
178 36 : SUBROUTINE restraint_functional(mulliken_restraint_control, charges, &
179 : charges_deriv, energy, order_p)
180 : TYPE(mulliken_restraint_type), INTENT(IN) :: mulliken_restraint_control
181 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges, charges_deriv
182 : REAL(KIND=dp), INTENT(OUT) :: energy, order_p
183 :
184 : INTEGER :: I
185 : REAL(KIND=dp) :: dum
186 :
187 540 : charges_deriv = 0.0_dp
188 36 : order_p = 0.0_dp
189 :
190 72 : DO I = 1, mulliken_restraint_control%natoms
191 : order_p = order_p + charges(mulliken_restraint_control%atoms(I), 1) &
192 72 : - charges(mulliken_restraint_control%atoms(I), 2) ! spin density on the relevant atoms
193 : END DO
194 : ! energy
195 36 : energy = mulliken_restraint_control%strength*(order_p - mulliken_restraint_control%target)**2
196 : ! derivative
197 36 : dum = 2*mulliken_restraint_control%strength*(order_p - mulliken_restraint_control%target)
198 72 : DO I = 1, mulliken_restraint_control%natoms
199 36 : charges_deriv(mulliken_restraint_control%atoms(I), 1) = dum
200 72 : charges_deriv(mulliken_restraint_control%atoms(I), 2) = -dum
201 : END DO
202 36 : END SUBROUTINE restraint_functional
203 :
204 : ! **************************************************************************************************
205 : !> \brief compute the mulliken charges
206 : !> \param p_matrix , s_matrix, para_env ...
207 : !> \param s_matrix ...
208 : !> \param charges previously allocated with the right size (natom,nspin)
209 : !> \param para_env ...
210 : !> \par History
211 : !> 06.2004 created [Joost VandeVondele]
212 : !> \note
213 : !> charges are computed per spin in the LSD case
214 : ! **************************************************************************************************
215 126370 : SUBROUTINE compute_charges(p_matrix, s_matrix, charges, para_env)
216 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
217 : TYPE(dbcsr_type), POINTER :: s_matrix
218 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges
219 : TYPE(mp_para_env_type), POINTER :: para_env
220 :
221 : INTEGER :: blk, iblock_col, iblock_row, ispin, nspin
222 : LOGICAL :: found
223 : REAL(kind=dp) :: mult
224 126370 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, s_block
225 : TYPE(dbcsr_iterator_type) :: iter
226 :
227 126370 : nspin = SIZE(p_matrix)
228 :
229 1314356 : charges = 0.0_dp
230 262364 : DO ispin = 1, nspin
231 135994 : CALL dbcsr_iterator_start(iter, s_matrix)
232 5270161 : DO WHILE (dbcsr_iterator_blocks_left(iter))
233 5134167 : NULLIFY (s_block, p_block)
234 5134167 : CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, s_block, blk)
235 : CALL dbcsr_get_block_p(matrix=p_matrix(ispin)%matrix, &
236 5134167 : row=iblock_row, col=iblock_col, BLOCK=p_block, found=found)
237 :
238 : ! we can cycle if a block is not present
239 5134167 : IF (.NOT. found) CYCLE
240 5134167 : IF (.NOT. (ASSOCIATED(s_block) .AND. ASSOCIATED(p_block))) CYCLE
241 :
242 5134167 : IF (iblock_row .EQ. iblock_col) THEN
243 : mult = 0.5_dp ! avoid double counting of diagonal blocks
244 : ELSE
245 4608019 : mult = 1.0_dp
246 : END IF
247 : charges(iblock_row, ispin) = charges(iblock_row, ispin) + &
248 80766525 : mult*SUM(p_block*s_block)
249 : charges(iblock_col, ispin) = charges(iblock_col, ispin) + &
250 80902519 : mult*SUM(p_block*s_block)
251 :
252 : END DO
253 398358 : CALL dbcsr_iterator_stop(iter)
254 : END DO
255 2502342 : CALL para_env%sum(charges)
256 :
257 126370 : END SUBROUTINE compute_charges
258 :
259 : ! **************************************************************************************************
260 : !> \brief compute the mulliken charges for single matrices
261 : !> \param p_matrix , s_matrix, para_env ...
262 : !> \param s_matrix ...
263 : !> \param charges previously allocated with the right size (natom,nspin)
264 : !> \param para_env ...
265 : !> \par History
266 : !> 06.2004 created [Joost VandeVondele]
267 : ! **************************************************************************************************
268 30 : SUBROUTINE compute_charges_single(p_matrix, s_matrix, charges, para_env)
269 : TYPE(dbcsr_type) :: p_matrix, s_matrix
270 : REAL(KIND=dp), DIMENSION(:) :: charges
271 : TYPE(mp_para_env_type), POINTER :: para_env
272 :
273 : INTEGER :: blk, iblock_col, iblock_row
274 : LOGICAL :: found
275 : REAL(kind=dp) :: mult
276 30 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, s_block
277 : TYPE(dbcsr_iterator_type) :: iter
278 :
279 124 : charges = 0.0_dp
280 30 : CALL dbcsr_iterator_start(iter, s_matrix)
281 128 : DO WHILE (dbcsr_iterator_blocks_left(iter))
282 98 : NULLIFY (s_block, p_block)
283 98 : CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, s_block, blk)
284 : CALL dbcsr_get_block_p(matrix=p_matrix, &
285 98 : row=iblock_row, col=iblock_col, BLOCK=p_block, found=found)
286 :
287 : ! we can cycle if a block is not present
288 98 : IF (.NOT. found) CYCLE
289 98 : IF (.NOT. (ASSOCIATED(s_block) .AND. ASSOCIATED(p_block))) CYCLE
290 :
291 98 : IF (iblock_row .EQ. iblock_col) THEN
292 : mult = 0.5_dp ! avoid double counting of diagonal blocks
293 : ELSE
294 51 : mult = 1.0_dp
295 : END IF
296 5690 : charges(iblock_row) = charges(iblock_row) + mult*SUM(p_block*s_block)
297 5720 : charges(iblock_col) = charges(iblock_col) + mult*SUM(p_block*s_block)
298 : END DO
299 30 : CALL dbcsr_iterator_stop(iter)
300 :
301 218 : CALL para_env%sum(charges)
302 :
303 30 : END SUBROUTINE compute_charges_single
304 :
305 : ! **************************************************************************************************
306 : !> \brief compute the mulliken charge derivatives
307 : !> \param p_matrix , s_matrix, para_env ...
308 : !> \param s_matrix ...
309 : !> \param charges ...
310 : !> \param dcharges previously allocated with the right size (natom,3)
311 : !> \param para_env ...
312 : !> \par History
313 : !> 01.2012 created [JHU]
314 : ! **************************************************************************************************
315 0 : SUBROUTINE compute_dcharges(p_matrix, s_matrix, charges, dcharges, para_env)
316 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix, s_matrix
317 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges, dcharges
318 : TYPE(mp_para_env_type), POINTER :: para_env
319 :
320 : INTEGER :: blk, iblock_col, iblock_row, ider, &
321 : ispin, nspin
322 : LOGICAL :: found
323 : REAL(kind=dp) :: mult
324 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ds_block, p_block, s_block
325 : TYPE(dbcsr_iterator_type) :: iter
326 :
327 0 : nspin = SIZE(p_matrix)
328 :
329 0 : charges = 0.0_dp
330 0 : dcharges = 0.0_dp
331 0 : DO ispin = 1, nspin
332 0 : CALL dbcsr_iterator_start(iter, s_matrix(1)%matrix)
333 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
334 0 : NULLIFY (s_block, p_block)
335 0 : CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, s_block, blk)
336 : CALL dbcsr_get_block_p(matrix=p_matrix(ispin)%matrix, &
337 0 : row=iblock_row, col=iblock_col, BLOCK=p_block, found=found)
338 :
339 : ! we can cycle if a block is not present
340 0 : IF (.NOT. found) CYCLE
341 0 : IF (.NOT. (ASSOCIATED(s_block) .AND. ASSOCIATED(p_block))) CYCLE
342 :
343 0 : IF (iblock_row .EQ. iblock_col) THEN
344 : mult = 0.5_dp ! avoid double counting of diagonal blocks
345 : ELSE
346 0 : mult = 1.0_dp
347 : END IF
348 0 : charges(iblock_row, ispin) = charges(iblock_row, ispin) + mult*SUM(p_block*s_block)
349 0 : charges(iblock_col, ispin) = charges(iblock_col, ispin) + mult*SUM(p_block*s_block)
350 0 : DO ider = 1, 3
351 : CALL dbcsr_get_block_p(matrix=s_matrix(ider + 1)%matrix, &
352 0 : row=iblock_row, col=iblock_col, BLOCK=ds_block, found=found)
353 0 : dcharges(iblock_row, ider) = dcharges(iblock_row, ider) + mult*SUM(p_block*ds_block)
354 0 : dcharges(iblock_col, ider) = dcharges(iblock_col, ider) + mult*SUM(p_block*ds_block)
355 : END DO
356 :
357 : END DO
358 0 : CALL dbcsr_iterator_stop(iter)
359 : END DO
360 0 : CALL para_env%sum(charges)
361 0 : CALL para_env%sum(dcharges)
362 :
363 0 : END SUBROUTINE compute_dcharges
364 :
365 : ! **************************************************************************************************
366 : !> \brief print the mulliken charges to scr on ionode
367 : !> \param p_matrix , s_matrix, para_env ...
368 : !> \param s_matrix ...
369 : !> \param para_env ...
370 : !> \param particle_set (needed for Z)
371 : !> \param qs_kind_set ...
372 : !> \param scr unit for output
373 : !> \param title ...
374 : !> \par History
375 : !> 06.2004 adapted to remove explicit matrix multiply [Joost VandeVondele]
376 : ! **************************************************************************************************
377 0 : SUBROUTINE mulliken_charges_a(p_matrix, s_matrix, para_env, particle_set, &
378 : qs_kind_set, scr, title)
379 :
380 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
381 : TYPE(dbcsr_type), POINTER :: s_matrix
382 : TYPE(mp_para_env_type), POINTER :: para_env
383 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
384 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
385 : INTEGER :: scr
386 : CHARACTER(LEN=*) :: title
387 :
388 : CHARACTER(len=*), PARAMETER :: routineN = 'mulliken_charges_a'
389 :
390 : INTEGER :: handle, nblock, nspin
391 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges
392 :
393 0 : CALL timeset(routineN, handle)
394 :
395 0 : CPASSERT(ASSOCIATED(p_matrix))
396 0 : CPASSERT(ASSOCIATED(s_matrix))
397 : ! here we get the numbers for charges
398 0 : nspin = SIZE(p_matrix)
399 0 : CALL dbcsr_get_info(s_matrix, nblkrows_total=nblock)
400 0 : ALLOCATE (charges(nblock, nspin))
401 :
402 0 : CALL compute_charges(p_matrix, s_matrix, charges, para_env)
403 :
404 0 : CALL print_atomic_charges(particle_set, qs_kind_set, scr, title, electronic_charges=charges)
405 :
406 0 : DEALLOCATE (charges)
407 :
408 0 : CALL timestop(handle)
409 :
410 0 : END SUBROUTINE mulliken_charges_a
411 :
412 : ! **************************************************************************************************
413 : !> \brief ...
414 : !> \param p_matrix ...
415 : !> \param s_matrix ...
416 : !> \param para_env ...
417 : !> \param mcharge ...
418 : ! **************************************************************************************************
419 52 : SUBROUTINE mulliken_charges_b(p_matrix, s_matrix, para_env, mcharge)
420 :
421 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
422 : TYPE(dbcsr_type), POINTER :: s_matrix
423 : TYPE(mp_para_env_type), POINTER :: para_env
424 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: mcharge
425 :
426 : CHARACTER(len=*), PARAMETER :: routineN = 'mulliken_charges_b'
427 :
428 : INTEGER :: handle
429 :
430 52 : CALL timeset(routineN, handle)
431 :
432 52 : IF (ASSOCIATED(p_matrix) .AND. ASSOCIATED(s_matrix)) THEN
433 52 : CALL compute_charges(p_matrix, s_matrix, mcharge, para_env)
434 : END IF
435 :
436 52 : CALL timestop(handle)
437 :
438 52 : END SUBROUTINE mulliken_charges_b
439 :
440 : ! **************************************************************************************************
441 : !> \brief ...
442 : !> \param p_matrix ...
443 : !> \param s_matrix ...
444 : !> \param para_env ...
445 : !> \param mcharge ...
446 : !> \param dmcharge ...
447 : ! **************************************************************************************************
448 0 : SUBROUTINE mulliken_charges_c(p_matrix, s_matrix, para_env, mcharge, dmcharge)
449 :
450 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix, s_matrix
451 : TYPE(mp_para_env_type), POINTER :: para_env
452 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: mcharge, dmcharge
453 :
454 : CHARACTER(len=*), PARAMETER :: routineN = 'mulliken_charges_c'
455 :
456 : INTEGER :: handle
457 :
458 0 : CALL timeset(routineN, handle)
459 :
460 0 : IF (ASSOCIATED(p_matrix) .AND. ASSOCIATED(s_matrix)) THEN
461 0 : CALL compute_dcharges(p_matrix, s_matrix, mcharge, dmcharge, para_env)
462 : END IF
463 :
464 0 : CALL timestop(handle)
465 :
466 0 : END SUBROUTINE mulliken_charges_c
467 :
468 : ! **************************************************************************************************
469 : !> \brief ...
470 : !> \param p_matrix ...
471 : !> \param s_matrix ...
472 : !> \param para_env ...
473 : !> \param mcharge ...
474 : ! **************************************************************************************************
475 30 : SUBROUTINE mulliken_charges_s(p_matrix, s_matrix, para_env, mcharge)
476 :
477 : TYPE(dbcsr_type) :: p_matrix, s_matrix
478 : TYPE(mp_para_env_type), POINTER :: para_env
479 : REAL(KIND=dp), DIMENSION(:) :: mcharge
480 :
481 : CHARACTER(len=*), PARAMETER :: routineN = 'mulliken_charges_s'
482 :
483 : INTEGER :: handle
484 :
485 30 : CALL timeset(routineN, handle)
486 :
487 30 : CALL compute_charges_single(p_matrix, s_matrix, mcharge, para_env)
488 :
489 30 : CALL timestop(handle)
490 :
491 30 : END SUBROUTINE mulliken_charges_s
492 :
493 : ! **************************************************************************************************
494 : !> \brief print the mulliken charges to scr on ionode
495 : !> \param p_matrix_kp ...
496 : !> \param s_matrix_kp ...
497 : !> \param para_env ...
498 : !> \param particle_set (needed for Z)
499 : !> \param qs_kind_set ...
500 : !> \param scr unit for output
501 : !> \param title ...
502 : !> \par History
503 : !> 06.2004 adapted to remove explicit matrix multiply [Joost VandeVondele]
504 : ! **************************************************************************************************
505 0 : SUBROUTINE mulliken_charges_akp(p_matrix_kp, s_matrix_kp, para_env, particle_set, &
506 : qs_kind_set, scr, title)
507 :
508 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: p_matrix_kp, s_matrix_kp
509 : TYPE(mp_para_env_type), POINTER :: para_env
510 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
511 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
512 : INTEGER :: scr
513 : CHARACTER(LEN=*) :: title
514 :
515 : CHARACTER(len=*), PARAMETER :: routineN = 'mulliken_charges_akp'
516 :
517 : INTEGER :: handle, ic, nblock, nspin
518 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges, charges_im
519 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
520 : TYPE(dbcsr_type), POINTER :: s_matrix
521 :
522 0 : CALL timeset(routineN, handle)
523 :
524 0 : CPASSERT(ASSOCIATED(p_matrix_kp))
525 0 : CPASSERT(ASSOCIATED(s_matrix_kp))
526 :
527 0 : nspin = SIZE(p_matrix_kp, 1)
528 0 : CALL dbcsr_get_info(s_matrix_kp(1, 1)%matrix, nblkrows_total=nblock)
529 0 : ALLOCATE (charges(nblock, nspin), charges_im(nblock, nspin))
530 0 : charges = 0.0_dp
531 :
532 0 : DO ic = 1, SIZE(s_matrix_kp, 2)
533 0 : NULLIFY (p_matrix, s_matrix)
534 0 : p_matrix => p_matrix_kp(:, ic)
535 0 : s_matrix => s_matrix_kp(1, ic)%matrix
536 0 : charges_im = 0.0_dp
537 0 : CALL compute_charges(p_matrix, s_matrix, charges_im, para_env)
538 0 : charges(:, :) = charges(:, :) + charges_im(:, :)
539 : END DO
540 :
541 0 : CALL print_atomic_charges(particle_set, qs_kind_set, scr, title, electronic_charges=charges)
542 :
543 0 : DEALLOCATE (charges, charges_im)
544 :
545 0 : CALL timestop(handle)
546 :
547 0 : END SUBROUTINE mulliken_charges_akp
548 :
549 : ! **************************************************************************************************
550 : !> \brief ...
551 : !> \param p_matrix_kp ...
552 : !> \param s_matrix_kp ...
553 : !> \param para_env ...
554 : !> \param mcharge ...
555 : ! **************************************************************************************************
556 17568 : SUBROUTINE mulliken_charges_bkp(p_matrix_kp, s_matrix_kp, para_env, mcharge)
557 :
558 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: p_matrix_kp, s_matrix_kp
559 : TYPE(mp_para_env_type), POINTER :: para_env
560 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: mcharge
561 :
562 : CHARACTER(len=*), PARAMETER :: routineN = 'mulliken_charges_bkp'
563 :
564 : INTEGER :: handle, ic, natom, nspin
565 17568 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: mcharge_im
566 17568 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
567 : TYPE(dbcsr_type), POINTER :: s_matrix
568 :
569 17568 : CALL timeset(routineN, handle)
570 :
571 17568 : IF (ASSOCIATED(p_matrix_kp) .AND. ASSOCIATED(s_matrix_kp)) THEN
572 :
573 250508 : mcharge = 0.0_dp
574 17568 : natom = SIZE(mcharge, 1)
575 17568 : nspin = SIZE(mcharge, 2)
576 70272 : ALLOCATE (mcharge_im(natom, nspin))
577 :
578 143850 : DO ic = 1, SIZE(s_matrix_kp, 2)
579 126282 : NULLIFY (p_matrix, s_matrix)
580 126282 : p_matrix => p_matrix_kp(:, ic)
581 126282 : s_matrix => s_matrix_kp(1, ic)%matrix
582 143850 : IF (ASSOCIATED(p_matrix) .AND. ASSOCIATED(s_matrix)) THEN
583 126282 : CALL compute_charges(p_matrix, s_matrix, mcharge_im, para_env)
584 2627112 : mcharge(:, :) = mcharge(:, :) + mcharge_im(:, :)
585 : END IF
586 : END DO
587 :
588 17568 : DEALLOCATE (mcharge_im)
589 :
590 : END IF
591 :
592 17568 : CALL timestop(handle)
593 :
594 17568 : END SUBROUTINE mulliken_charges_bkp
595 :
596 : ! **************************************************************************************************
597 : !> \brief ...
598 : !> \param p_matrix_kp ...
599 : !> \param s_matrix_kp ...
600 : !> \param para_env ...
601 : !> \param mcharge ...
602 : !> \param dmcharge ...
603 : ! **************************************************************************************************
604 0 : SUBROUTINE mulliken_charges_ckp(p_matrix_kp, s_matrix_kp, para_env, &
605 : mcharge, dmcharge)
606 :
607 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: p_matrix_kp, s_matrix_kp
608 : TYPE(mp_para_env_type), POINTER :: para_env
609 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: mcharge, dmcharge
610 :
611 : CHARACTER(len=*), PARAMETER :: routineN = 'mulliken_charges_ckp'
612 :
613 : INTEGER :: handle, ic, natom, nder, nspin
614 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dmcharge_im, mcharge_im
615 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix, s_matrix
616 :
617 0 : CALL timeset(routineN, handle)
618 :
619 0 : IF (ASSOCIATED(p_matrix_kp) .AND. ASSOCIATED(s_matrix_kp)) THEN
620 :
621 0 : mcharge = 0.0_dp
622 0 : dmcharge = 0.0_dp
623 0 : natom = SIZE(mcharge, 1)
624 0 : nspin = SIZE(mcharge, 2)
625 0 : nder = SIZE(dmcharge, 2)
626 0 : ALLOCATE (mcharge_im(natom, nspin), dmcharge_im(natom, nder))
627 :
628 0 : DO ic = 1, SIZE(s_matrix_kp, 2)
629 : NULLIFY (p_matrix, s_matrix)
630 0 : p_matrix => p_matrix_kp(:, ic)
631 0 : s_matrix => s_matrix_kp(:, ic)
632 0 : IF (ASSOCIATED(p_matrix) .AND. ASSOCIATED(s_matrix)) THEN
633 0 : CALL compute_dcharges(p_matrix, s_matrix, mcharge_im, dmcharge_im, para_env)
634 0 : mcharge(:, :) = mcharge(:, :) + mcharge_im(:, :)
635 0 : dmcharge(:, :) = dmcharge(:, :) + dmcharge_im(:, :)
636 : END IF
637 : END DO
638 :
639 0 : DEALLOCATE (mcharge_im, dmcharge_im)
640 :
641 : END IF
642 :
643 0 : CALL timestop(handle)
644 :
645 0 : END SUBROUTINE mulliken_charges_ckp
646 :
647 : ! **************************************************************************************************
648 : !> \brief compute Mayer bond orders for a single spin channel
649 : !> for complete result sum up over all spins and multiply by Nspin
650 : !> \param psmat ...
651 : !> \param spmat ...
652 : !> \param bond_order ...
653 : !> \par History
654 : !> 12.2016 created [JGH]
655 : ! **************************************************************************************************
656 0 : SUBROUTINE compute_bond_order(psmat, spmat, bond_order)
657 : TYPE(dbcsr_type) :: psmat, spmat
658 : REAL(KIND=dp), DIMENSION(:, :) :: bond_order
659 :
660 : INTEGER :: blk, iat, jat
661 : LOGICAL :: found
662 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ps, sp
663 : TYPE(dbcsr_iterator_type) :: iter
664 :
665 0 : CALL dbcsr_iterator_start(iter, psmat)
666 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
667 0 : NULLIFY (ps, sp)
668 0 : CALL dbcsr_iterator_next_block(iter, iat, jat, ps, blk)
669 : CALL dbcsr_get_block_p(matrix=spmat, &
670 0 : row=iat, col=jat, BLOCK=sp, found=found)
671 0 : IF (.NOT. found) CYCLE
672 0 : IF (.NOT. (ASSOCIATED(sp) .AND. ASSOCIATED(ps))) CYCLE
673 :
674 0 : bond_order(iat, jat) = bond_order(iat, jat) + SUM(ps*sp)
675 : END DO
676 0 : CALL dbcsr_iterator_stop(iter)
677 :
678 0 : END SUBROUTINE compute_bond_order
679 :
680 : ! **************************************************************************************************
681 : !> \brief compute the mulliken charges for a single atom (AO resolved)
682 : !> \param p_matrix , s_matrix, para_env ...
683 : !> \param s_matrix ...
684 : !> \param charges ...
685 : !> \param iatom ...
686 : !> \param para_env ...
687 : !> \par History
688 : !> 06.2004 created [Joost VandeVondele]
689 : !> 10.2018 adapted [JGH]
690 : !> \note
691 : ! **************************************************************************************************
692 0 : SUBROUTINE ao_charges_1(p_matrix, s_matrix, charges, iatom, para_env)
693 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
694 : TYPE(dbcsr_type), POINTER :: s_matrix
695 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: charges
696 : INTEGER, INTENT(IN) :: iatom
697 : TYPE(mp_para_env_type), POINTER :: para_env
698 :
699 : CHARACTER(len=*), PARAMETER :: routineN = 'ao_charges_1'
700 :
701 : INTEGER :: blk, handle, i, iblock_col, iblock_row, &
702 : ispin, j, nspin
703 : LOGICAL :: found
704 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, s_block
705 : TYPE(dbcsr_iterator_type) :: iter
706 :
707 0 : CALL timeset(routineN, handle)
708 :
709 0 : nspin = SIZE(p_matrix)
710 0 : charges = 0.0_dp
711 0 : DO ispin = 1, nspin
712 0 : CALL dbcsr_iterator_start(iter, s_matrix)
713 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
714 0 : NULLIFY (s_block, p_block)
715 0 : CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, s_block, blk)
716 : CALL dbcsr_get_block_p(matrix=p_matrix(ispin)%matrix, &
717 0 : row=iblock_row, col=iblock_col, BLOCK=p_block, found=found)
718 :
719 : ! we can cycle if a block is not present
720 0 : IF (.NOT. found) CYCLE
721 0 : IF (.NOT. (ASSOCIATED(s_block) .AND. ASSOCIATED(p_block))) CYCLE
722 :
723 0 : IF (iblock_row == iatom) THEN
724 0 : DO j = 1, SIZE(p_block, 2)
725 0 : DO i = 1, SIZE(p_block, 1)
726 0 : charges(i) = charges(i) + p_block(i, j)*s_block(i, j)
727 : END DO
728 : END DO
729 0 : ELSEIF (iblock_col == iatom) THEN
730 0 : DO j = 1, SIZE(p_block, 2)
731 0 : DO i = 1, SIZE(p_block, 1)
732 0 : charges(j) = charges(j) + p_block(i, j)*s_block(i, j)
733 : END DO
734 : END DO
735 : END IF
736 : END DO
737 0 : CALL dbcsr_iterator_stop(iter)
738 : END DO
739 0 : CALL para_env%sum(charges)
740 :
741 0 : CALL timestop(handle)
742 :
743 0 : END SUBROUTINE ao_charges_1
744 :
745 : ! **************************************************************************************************
746 : !> \brief compute the mulliken charges (AO resolved)
747 : !> \param p_matrix , s_matrix, para_env ...
748 : !> \param s_matrix ...
749 : !> \param charges ...
750 : !> \param para_env ...
751 : !> \par History
752 : !> 06.2004 created [Joost VandeVondele]
753 : !> 10.2018 adapted [JGH]
754 : !> \note
755 : ! **************************************************************************************************
756 327800 : SUBROUTINE ao_charges_2(p_matrix, s_matrix, charges, para_env)
757 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
758 : TYPE(dbcsr_type), POINTER :: s_matrix
759 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
760 : TYPE(mp_para_env_type), POINTER :: para_env
761 :
762 : CHARACTER(len=*), PARAMETER :: routineN = 'ao_charges_2'
763 :
764 : INTEGER :: blk, handle, i, iblock_col, iblock_row, &
765 : ispin, j, nspin
766 : LOGICAL :: found
767 327800 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, s_block
768 : TYPE(dbcsr_iterator_type) :: iter
769 :
770 327800 : CALL timeset(routineN, handle)
771 :
772 327800 : nspin = SIZE(p_matrix)
773 14057130 : charges = 0.0_dp
774 712202 : DO ispin = 1, nspin
775 384402 : CALL dbcsr_iterator_start(iter, s_matrix)
776 7960453 : DO WHILE (dbcsr_iterator_blocks_left(iter))
777 7576051 : NULLIFY (s_block, p_block)
778 7576051 : CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, s_block, blk)
779 : CALL dbcsr_get_block_p(matrix=p_matrix(ispin)%matrix, &
780 7576051 : row=iblock_row, col=iblock_col, BLOCK=p_block, found=found)
781 :
782 : ! we can cycle if a block is not present
783 7576051 : IF (.NOT. found) CYCLE
784 7576051 : IF (.NOT. (ASSOCIATED(s_block) .AND. ASSOCIATED(p_block))) CYCLE
785 :
786 36457538 : DO j = 1, SIZE(p_block, 2)
787 195094039 : DO i = 1, SIZE(p_block, 1)
788 187517988 : charges(i, iblock_row) = charges(i, iblock_row) + p_block(i, j)*s_block(i, j)
789 : END DO
790 : END DO
791 7960453 : IF (iblock_col /= iblock_row) THEN
792 29957177 : DO j = 1, SIZE(p_block, 2)
793 153913599 : DO i = 1, SIZE(p_block, 1)
794 147487066 : charges(j, iblock_col) = charges(j, iblock_col) + p_block(i, j)*s_block(i, j)
795 : END DO
796 : END DO
797 : END IF
798 :
799 : END DO
800 1096604 : CALL dbcsr_iterator_stop(iter)
801 : END DO
802 27786460 : CALL para_env%sum(charges)
803 :
804 327800 : CALL timestop(handle)
805 :
806 327800 : END SUBROUTINE ao_charges_2
807 :
808 : ! **************************************************************************************************
809 : !> \brief ...
810 : !> \param p_matrix_kp ...
811 : !> \param s_matrix_kp ...
812 : !> \param charges ...
813 : !> \param iatom ...
814 : !> \param para_env ...
815 : ! **************************************************************************************************
816 0 : SUBROUTINE ao_charges_kp(p_matrix_kp, s_matrix_kp, charges, iatom, para_env)
817 :
818 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: p_matrix_kp, s_matrix_kp
819 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: charges
820 : INTEGER, INTENT(IN) :: iatom
821 : TYPE(mp_para_env_type), POINTER :: para_env
822 :
823 : CHARACTER(len=*), PARAMETER :: routineN = 'ao_charges_kp'
824 :
825 : INTEGER :: handle, ic, na
826 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: charge_im
827 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
828 : TYPE(dbcsr_type), POINTER :: s_matrix
829 :
830 0 : CALL timeset(routineN, handle)
831 :
832 0 : IF (ASSOCIATED(p_matrix_kp) .AND. ASSOCIATED(s_matrix_kp)) THEN
833 :
834 0 : charges = 0.0_dp
835 0 : na = SIZE(charges, 1)
836 0 : ALLOCATE (charge_im(na))
837 :
838 0 : DO ic = 1, SIZE(s_matrix_kp, 2)
839 0 : NULLIFY (p_matrix, s_matrix)
840 0 : p_matrix => p_matrix_kp(:, ic)
841 0 : s_matrix => s_matrix_kp(1, ic)%matrix
842 0 : IF (ASSOCIATED(p_matrix) .AND. ASSOCIATED(s_matrix)) THEN
843 0 : CALL ao_charges_1(p_matrix, s_matrix, charge_im, iatom, para_env)
844 0 : charges(:) = charges(:) + charge_im(:)
845 : END IF
846 : END DO
847 :
848 0 : DEALLOCATE (charge_im)
849 :
850 : END IF
851 :
852 0 : CALL timestop(handle)
853 :
854 0 : END SUBROUTINE ao_charges_kp
855 :
856 : ! **************************************************************************************************
857 : !> \brief ...
858 : !> \param p_matrix_kp ...
859 : !> \param s_matrix_kp ...
860 : !> \param charges ...
861 : !> \param para_env ...
862 : ! **************************************************************************************************
863 2856 : SUBROUTINE ao_charges_kp_2(p_matrix_kp, s_matrix_kp, charges, para_env)
864 :
865 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: p_matrix_kp, s_matrix_kp
866 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
867 : TYPE(mp_para_env_type), POINTER :: para_env
868 :
869 : CHARACTER(len=*), PARAMETER :: routineN = 'ao_charges_kp_2'
870 :
871 : INTEGER :: handle, ic, na, nb
872 2856 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: charge_im
873 2856 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
874 : TYPE(dbcsr_type), POINTER :: s_matrix
875 :
876 2856 : CALL timeset(routineN, handle)
877 :
878 2856 : IF (ASSOCIATED(p_matrix_kp) .AND. ASSOCIATED(s_matrix_kp)) THEN
879 :
880 120606 : charges = 0.0_dp
881 2856 : na = SIZE(charges, 1)
882 2856 : nb = SIZE(charges, 2)
883 11424 : ALLOCATE (charge_im(na, nb))
884 :
885 312072 : DO ic = 1, SIZE(s_matrix_kp, 2)
886 309216 : NULLIFY (p_matrix, s_matrix)
887 309216 : p_matrix => p_matrix_kp(:, ic)
888 309216 : s_matrix => s_matrix_kp(1, ic)%matrix
889 312072 : IF (ASSOCIATED(p_matrix) .AND. ASSOCIATED(s_matrix)) THEN
890 309216 : CALL ao_charges_2(p_matrix, s_matrix, charge_im, para_env)
891 13064726 : charges(:, :) = charges(:, :) + charge_im(:, :)
892 : END IF
893 : END DO
894 :
895 2856 : DEALLOCATE (charge_im)
896 :
897 : END IF
898 :
899 2856 : CALL timestop(handle)
900 :
901 5712 : END SUBROUTINE ao_charges_kp_2
902 :
903 : ! **************************************************************************************************
904 : !> \brief Compute partial trace of product of two matrices
905 : !> \param amat ...
906 : !> \param bmat ...
907 : !> \param factor ...
908 : !> \param atrace ...
909 : !> \par History
910 : !> 06.2004 created [Joost VandeVondele]
911 : !> \note
912 : !> charges are computed per spin in the LSD case
913 : ! **************************************************************************************************
914 736 : SUBROUTINE atom_trace(amat, bmat, factor, atrace)
915 : TYPE(dbcsr_type), POINTER :: amat, bmat
916 : REAL(kind=dp), INTENT(IN) :: factor
917 : REAL(KIND=dp), DIMENSION(:), POINTER :: atrace
918 :
919 : INTEGER :: blk, iblock_col, iblock_row, nblock
920 : LOGICAL :: found
921 : REAL(kind=dp) :: btr, mult
922 368 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a_block, b_block
923 : TYPE(dbcsr_iterator_type) :: iter
924 :
925 368 : CALL dbcsr_get_info(bmat, nblkrows_total=nblock)
926 368 : CPASSERT(nblock == SIZE(atrace))
927 :
928 368 : CALL dbcsr_iterator_start(iter, bmat)
929 279476 : DO WHILE (dbcsr_iterator_blocks_left(iter))
930 279108 : CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, b_block, blk)
931 : CALL dbcsr_get_block_p(matrix=amat, &
932 279108 : row=iblock_row, col=iblock_col, BLOCK=a_block, found=found)
933 :
934 : ! we can cycle if a block is not present
935 279108 : IF (.NOT. (ASSOCIATED(b_block) .AND. ASSOCIATED(a_block))) CYCLE
936 :
937 279108 : IF (iblock_row .EQ. iblock_col) THEN
938 : mult = 0.5_dp ! avoid double counting of diagonal blocks
939 : ELSE
940 271596 : mult = 1.0_dp
941 : END IF
942 2182926 : btr = factor*mult*SUM(a_block*b_block)
943 279108 : atrace(iblock_row) = atrace(iblock_row) + btr
944 279108 : atrace(iblock_col) = atrace(iblock_col) + btr
945 :
946 : END DO
947 368 : CALL dbcsr_iterator_stop(iter)
948 :
949 368 : END SUBROUTINE atom_trace
950 :
951 : ! **************************************************************************************************
952 :
953 : END MODULE mulliken
|
