(git:374b731)
Loading...
Searching...
No Matches
mulliken.F
Go to the documentation of this file.
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! **************************************************************************************************
13MODULE mulliken
14 USE atomic_charges, ONLY: print_atomic_charges
15 USE cp_control_types, ONLY: mulliken_restraint_type
16 USE 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
46CONTAINS
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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: charges, charges_deriv, ks_block, &
86 p_block, s_block, w_block
87 TYPE(dbcsr_iterator_type) :: iter
88
89! here we get the numbers for charges
90
91 nspin = SIZE(p_matrix)
92 CALL dbcsr_get_info(s_matrix, nblkrows_total=nblock)
93
94 ALLOCATE (charges(nblock, nspin))
95 ALLOCATE (charges_deriv(nblock, nspin))
96 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 charges, charges_deriv, my_energy, my_order_p)
103
104 IF (PRESENT(order_p)) THEN
105 order_p = my_order_p
106 END IF
107 IF (PRESENT(energy)) THEN
108 energy = my_energy
109 END IF
110
111 IF (PRESENT(ks_matrix)) THEN
112
113 DO ispin = 1, nspin
114 CALL dbcsr_iterator_start(iter, s_matrix)
115 DO WHILE (dbcsr_iterator_blocks_left(iter))
116 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 row=iblock_row, col=iblock_col, block=ks_block, found=found)
119
120 IF (.NOT. (ASSOCIATED(s_block) .AND. ASSOCIATED(ks_block))) THEN
121 cpabort("Unexpected s / ks structure")
122 END IF
123 mult = 0.5_dp*charges_deriv(iblock_row, ispin) + &
124 0.5_dp*charges_deriv(iblock_col, ispin)
125
126 ks_block = ks_block + mult*s_block
127
128 END DO
129 CALL dbcsr_iterator_stop(iter)
130 END DO
131
132 END IF
133
134 IF (PRESENT(w_matrix)) THEN
135
136 DO ispin = 1, nspin
137 CALL dbcsr_iterator_start(iter, p_matrix(ispin)%matrix)
138 DO WHILE (dbcsr_iterator_blocks_left(iter))
139 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 row=iblock_row, col=iblock_col, block=w_block, found=found)
142
143 ! we can cycle if a block is not present
144 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 - 0.5_dp*charges_deriv(iblock_col, ispin)
149
150 w_block = w_block + mult*p_block
151
152 END DO
153 CALL dbcsr_iterator_stop(iter)
154 END DO
155
156 END IF
157
158 DEALLOCATE (charges)
159 DEALLOCATE (charges_deriv)
160
161 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 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 charges_deriv = 0.0_dp
188 order_p = 0.0_dp
189
190 DO i = 1, mulliken_restraint_control%natoms
191 order_p = order_p + charges(mulliken_restraint_control%atoms(i), 1) &
192 - charges(mulliken_restraint_control%atoms(i), 2) ! spin density on the relevant atoms
193 END DO
194 ! energy
195 energy = mulliken_restraint_control%strength*(order_p - mulliken_restraint_control%target)**2
196 ! derivative
197 dum = 2*mulliken_restraint_control%strength*(order_p - mulliken_restraint_control%target)
198 DO i = 1, mulliken_restraint_control%natoms
199 charges_deriv(mulliken_restraint_control%atoms(i), 1) = dum
200 charges_deriv(mulliken_restraint_control%atoms(i), 2) = -dum
201 END DO
202 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: p_block, s_block
225 TYPE(dbcsr_iterator_type) :: iter
226
227 nspin = SIZE(p_matrix)
228
229 charges = 0.0_dp
230 DO ispin = 1, nspin
231 CALL dbcsr_iterator_start(iter, s_matrix)
232 DO WHILE (dbcsr_iterator_blocks_left(iter))
233 NULLIFY (s_block, p_block)
234 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 row=iblock_row, col=iblock_col, block=p_block, found=found)
237
238 ! we can cycle if a block is not present
239 IF (.NOT. found) cycle
240 IF (.NOT. (ASSOCIATED(s_block) .AND. ASSOCIATED(p_block))) cycle
241
242 IF (iblock_row .EQ. iblock_col) THEN
243 mult = 0.5_dp ! avoid double counting of diagonal blocks
244 ELSE
245 mult = 1.0_dp
246 END IF
247 charges(iblock_row, ispin) = charges(iblock_row, ispin) + &
248 mult*sum(p_block*s_block)
249 charges(iblock_col, ispin) = charges(iblock_col, ispin) + &
250 mult*sum(p_block*s_block)
251
252 END DO
253 CALL dbcsr_iterator_stop(iter)
254 END DO
255 CALL para_env%sum(charges)
256
257 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: p_block, s_block
277 TYPE(dbcsr_iterator_type) :: iter
278
279 charges = 0.0_dp
280 CALL dbcsr_iterator_start(iter, s_matrix)
281 DO WHILE (dbcsr_iterator_blocks_left(iter))
282 NULLIFY (s_block, p_block)
283 CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, s_block, blk)
284 CALL dbcsr_get_block_p(matrix=p_matrix, &
285 row=iblock_row, col=iblock_col, block=p_block, found=found)
286
287 ! we can cycle if a block is not present
288 IF (.NOT. found) cycle
289 IF (.NOT. (ASSOCIATED(s_block) .AND. ASSOCIATED(p_block))) cycle
290
291 IF (iblock_row .EQ. iblock_col) THEN
292 mult = 0.5_dp ! avoid double counting of diagonal blocks
293 ELSE
294 mult = 1.0_dp
295 END IF
296 charges(iblock_row) = charges(iblock_row) + mult*sum(p_block*s_block)
297 charges(iblock_col) = charges(iblock_col) + mult*sum(p_block*s_block)
298 END DO
299 CALL dbcsr_iterator_stop(iter)
300
301 CALL para_env%sum(charges)
302
303 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: ds_block, p_block, s_block
325 TYPE(dbcsr_iterator_type) :: iter
326
327 nspin = SIZE(p_matrix)
328
329 charges = 0.0_dp
330 dcharges = 0.0_dp
331 DO ispin = 1, nspin
332 CALL dbcsr_iterator_start(iter, s_matrix(1)%matrix)
333 DO WHILE (dbcsr_iterator_blocks_left(iter))
334 NULLIFY (s_block, p_block)
335 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 row=iblock_row, col=iblock_col, block=p_block, found=found)
338
339 ! we can cycle if a block is not present
340 IF (.NOT. found) cycle
341 IF (.NOT. (ASSOCIATED(s_block) .AND. ASSOCIATED(p_block))) cycle
342
343 IF (iblock_row .EQ. iblock_col) THEN
344 mult = 0.5_dp ! avoid double counting of diagonal blocks
345 ELSE
346 mult = 1.0_dp
347 END IF
348 charges(iblock_row, ispin) = charges(iblock_row, ispin) + mult*sum(p_block*s_block)
349 charges(iblock_col, ispin) = charges(iblock_col, ispin) + mult*sum(p_block*s_block)
350 DO ider = 1, 3
351 CALL dbcsr_get_block_p(matrix=s_matrix(ider + 1)%matrix, &
352 row=iblock_row, col=iblock_col, block=ds_block, found=found)
353 dcharges(iblock_row, ider) = dcharges(iblock_row, ider) + mult*sum(p_block*ds_block)
354 dcharges(iblock_col, ider) = dcharges(iblock_col, ider) + mult*sum(p_block*ds_block)
355 END DO
356
357 END DO
358 CALL dbcsr_iterator_stop(iter)
359 END DO
360 CALL para_env%sum(charges)
361 CALL para_env%sum(dcharges)
362
363 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 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 REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges
392
393 CALL timeset(routinen, handle)
394
395 cpassert(ASSOCIATED(p_matrix))
396 cpassert(ASSOCIATED(s_matrix))
397 ! here we get the numbers for charges
398 nspin = SIZE(p_matrix)
399 CALL dbcsr_get_info(s_matrix, nblkrows_total=nblock)
400 ALLOCATE (charges(nblock, nspin))
401
402 CALL compute_charges(p_matrix, s_matrix, charges, para_env)
403
404 CALL print_atomic_charges(particle_set, qs_kind_set, scr, title, electronic_charges=charges)
405
406 DEALLOCATE (charges)
407
408 CALL timestop(handle)
409
410 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 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 CALL timeset(routinen, handle)
431
432 IF (ASSOCIATED(p_matrix) .AND. ASSOCIATED(s_matrix)) THEN
433 CALL compute_charges(p_matrix, s_matrix, mcharge, para_env)
434 END IF
435
436 CALL timestop(handle)
437
438 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 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 CALL timeset(routinen, handle)
459
460 IF (ASSOCIATED(p_matrix) .AND. ASSOCIATED(s_matrix)) THEN
461 CALL compute_dcharges(p_matrix, s_matrix, mcharge, dmcharge, para_env)
462 END IF
463
464 CALL timestop(handle)
465
466 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 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 CALL timeset(routinen, handle)
486
487 CALL compute_charges_single(p_matrix, s_matrix, mcharge, para_env)
488
489 CALL timestop(handle)
490
491 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 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 REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges, charges_im
519 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
520 TYPE(dbcsr_type), POINTER :: s_matrix
521
522 CALL timeset(routinen, handle)
523
524 cpassert(ASSOCIATED(p_matrix_kp))
525 cpassert(ASSOCIATED(s_matrix_kp))
526
527 nspin = SIZE(p_matrix_kp, 1)
528 CALL dbcsr_get_info(s_matrix_kp(1, 1)%matrix, nblkrows_total=nblock)
529 ALLOCATE (charges(nblock, nspin), charges_im(nblock, nspin))
530 charges = 0.0_dp
531
532 DO ic = 1, SIZE(s_matrix_kp, 2)
533 NULLIFY (p_matrix, s_matrix)
534 p_matrix => p_matrix_kp(:, ic)
535 s_matrix => s_matrix_kp(1, ic)%matrix
536 charges_im = 0.0_dp
537 CALL compute_charges(p_matrix, s_matrix, charges_im, para_env)
538 charges(:, :) = charges(:, :) + charges_im(:, :)
539 END DO
540
541 CALL print_atomic_charges(particle_set, qs_kind_set, scr, title, electronic_charges=charges)
542
543 DEALLOCATE (charges, charges_im)
544
545 CALL timestop(handle)
546
547 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 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 REAL(KIND=dp), DIMENSION(:, :), POINTER :: mcharge_im
566 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
567 TYPE(dbcsr_type), POINTER :: s_matrix
568
569 CALL timeset(routinen, handle)
570
571 IF (ASSOCIATED(p_matrix_kp) .AND. ASSOCIATED(s_matrix_kp)) THEN
572
573 mcharge = 0.0_dp
574 natom = SIZE(mcharge, 1)
575 nspin = SIZE(mcharge, 2)
576 ALLOCATE (mcharge_im(natom, nspin))
577
578 DO ic = 1, SIZE(s_matrix_kp, 2)
579 NULLIFY (p_matrix, s_matrix)
580 p_matrix => p_matrix_kp(:, ic)
581 s_matrix => s_matrix_kp(1, ic)%matrix
582 IF (ASSOCIATED(p_matrix) .AND. ASSOCIATED(s_matrix)) THEN
583 CALL compute_charges(p_matrix, s_matrix, mcharge_im, para_env)
584 mcharge(:, :) = mcharge(:, :) + mcharge_im(:, :)
585 END IF
586 END DO
587
588 DEALLOCATE (mcharge_im)
589
590 END IF
591
592 CALL timestop(handle)
593
594 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 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 REAL(KIND=dp), DIMENSION(:, :), POINTER :: dmcharge_im, mcharge_im
615 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix, s_matrix
616
617 CALL timeset(routinen, handle)
618
619 IF (ASSOCIATED(p_matrix_kp) .AND. ASSOCIATED(s_matrix_kp)) THEN
620
621 mcharge = 0.0_dp
622 dmcharge = 0.0_dp
623 natom = SIZE(mcharge, 1)
624 nspin = SIZE(mcharge, 2)
625 nder = SIZE(dmcharge, 2)
626 ALLOCATE (mcharge_im(natom, nspin), dmcharge_im(natom, nder))
627
628 DO ic = 1, SIZE(s_matrix_kp, 2)
629 NULLIFY (p_matrix, s_matrix)
630 p_matrix => p_matrix_kp(:, ic)
631 s_matrix => s_matrix_kp(:, ic)
632 IF (ASSOCIATED(p_matrix) .AND. ASSOCIATED(s_matrix)) THEN
633 CALL compute_dcharges(p_matrix, s_matrix, mcharge_im, dmcharge_im, para_env)
634 mcharge(:, :) = mcharge(:, :) + mcharge_im(:, :)
635 dmcharge(:, :) = dmcharge(:, :) + dmcharge_im(:, :)
636 END IF
637 END DO
638
639 DEALLOCATE (mcharge_im, dmcharge_im)
640
641 END IF
642
643 CALL timestop(handle)
644
645 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: ps, sp
663 TYPE(dbcsr_iterator_type) :: iter
664
665 CALL dbcsr_iterator_start(iter, psmat)
666 DO WHILE (dbcsr_iterator_blocks_left(iter))
667 NULLIFY (ps, sp)
668 CALL dbcsr_iterator_next_block(iter, iat, jat, ps, blk)
669 CALL dbcsr_get_block_p(matrix=spmat, &
670 row=iat, col=jat, block=sp, found=found)
671 IF (.NOT. found) cycle
672 IF (.NOT. (ASSOCIATED(sp) .AND. ASSOCIATED(ps))) cycle
673
674 bond_order(iat, jat) = bond_order(iat, jat) + sum(ps*sp)
675 END DO
676 CALL dbcsr_iterator_stop(iter)
677
678 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 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 REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, s_block
705 TYPE(dbcsr_iterator_type) :: iter
706
707 CALL timeset(routinen, handle)
708
709 nspin = SIZE(p_matrix)
710 charges = 0.0_dp
711 DO ispin = 1, nspin
712 CALL dbcsr_iterator_start(iter, s_matrix)
713 DO WHILE (dbcsr_iterator_blocks_left(iter))
714 NULLIFY (s_block, p_block)
715 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 row=iblock_row, col=iblock_col, block=p_block, found=found)
718
719 ! we can cycle if a block is not present
720 IF (.NOT. found) cycle
721 IF (.NOT. (ASSOCIATED(s_block) .AND. ASSOCIATED(p_block))) cycle
722
723 IF (iblock_row == iatom) THEN
724 DO j = 1, SIZE(p_block, 2)
725 DO i = 1, SIZE(p_block, 1)
726 charges(i) = charges(i) + p_block(i, j)*s_block(i, j)
727 END DO
728 END DO
729 ELSEIF (iblock_col == iatom) THEN
730 DO j = 1, SIZE(p_block, 2)
731 DO i = 1, SIZE(p_block, 1)
732 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 CALL dbcsr_iterator_stop(iter)
738 END DO
739 CALL para_env%sum(charges)
740
741 CALL timestop(handle)
742
743 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: p_block, s_block
768 TYPE(dbcsr_iterator_type) :: iter
769
770 CALL timeset(routinen, handle)
771
772 nspin = SIZE(p_matrix)
773 charges = 0.0_dp
774 DO ispin = 1, nspin
775 CALL dbcsr_iterator_start(iter, s_matrix)
776 DO WHILE (dbcsr_iterator_blocks_left(iter))
777 NULLIFY (s_block, p_block)
778 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 row=iblock_row, col=iblock_col, block=p_block, found=found)
781
782 ! we can cycle if a block is not present
783 IF (.NOT. found) cycle
784 IF (.NOT. (ASSOCIATED(s_block) .AND. ASSOCIATED(p_block))) cycle
785
786 DO j = 1, SIZE(p_block, 2)
787 DO i = 1, SIZE(p_block, 1)
788 charges(i, iblock_row) = charges(i, iblock_row) + p_block(i, j)*s_block(i, j)
789 END DO
790 END DO
791 IF (iblock_col /= iblock_row) THEN
792 DO j = 1, SIZE(p_block, 2)
793 DO i = 1, SIZE(p_block, 1)
794 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 CALL dbcsr_iterator_stop(iter)
801 END DO
802 CALL para_env%sum(charges)
803
804 CALL timestop(handle)
805
806 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 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 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: charge_im
827 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
828 TYPE(dbcsr_type), POINTER :: s_matrix
829
830 CALL timeset(routinen, handle)
831
832 IF (ASSOCIATED(p_matrix_kp) .AND. ASSOCIATED(s_matrix_kp)) THEN
833
834 charges = 0.0_dp
835 na = SIZE(charges, 1)
836 ALLOCATE (charge_im(na))
837
838 DO ic = 1, SIZE(s_matrix_kp, 2)
839 NULLIFY (p_matrix, s_matrix)
840 p_matrix => p_matrix_kp(:, ic)
841 s_matrix => s_matrix_kp(1, ic)%matrix
842 IF (ASSOCIATED(p_matrix) .AND. ASSOCIATED(s_matrix)) THEN
843 CALL ao_charges_1(p_matrix, s_matrix, charge_im, iatom, para_env)
844 charges(:) = charges(:) + charge_im(:)
845 END IF
846 END DO
847
848 DEALLOCATE (charge_im)
849
850 END IF
851
852 CALL timestop(handle)
853
854 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: charge_im
873 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
874 TYPE(dbcsr_type), POINTER :: s_matrix
875
876 CALL timeset(routinen, handle)
877
878 IF (ASSOCIATED(p_matrix_kp) .AND. ASSOCIATED(s_matrix_kp)) THEN
879
880 charges = 0.0_dp
881 na = SIZE(charges, 1)
882 nb = SIZE(charges, 2)
883 ALLOCATE (charge_im(na, nb))
884
885 DO ic = 1, SIZE(s_matrix_kp, 2)
886 NULLIFY (p_matrix, s_matrix)
887 p_matrix => p_matrix_kp(:, ic)
888 s_matrix => s_matrix_kp(1, ic)%matrix
889 IF (ASSOCIATED(p_matrix) .AND. ASSOCIATED(s_matrix)) THEN
890 CALL ao_charges_2(p_matrix, s_matrix, charge_im, para_env)
891 charges(:, :) = charges(:, :) + charge_im(:, :)
892 END IF
893 END DO
894
895 DEALLOCATE (charge_im)
896
897 END IF
898
899 CALL timestop(handle)
900
901 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: a_block, b_block
923 TYPE(dbcsr_iterator_type) :: iter
924
925 CALL dbcsr_get_info(bmat, nblkrows_total=nblock)
926 cpassert(nblock == SIZE(atrace))
927
928 CALL dbcsr_iterator_start(iter, bmat)
929 DO WHILE (dbcsr_iterator_blocks_left(iter))
930 CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, b_block, blk)
931 CALL dbcsr_get_block_p(matrix=amat, &
932 row=iblock_row, col=iblock_col, block=a_block, found=found)
933
934 ! we can cycle if a block is not present
935 IF (.NOT. (ASSOCIATED(b_block) .AND. ASSOCIATED(a_block))) cycle
936
937 IF (iblock_row .EQ. iblock_col) THEN
938 mult = 0.5_dp ! avoid double counting of diagonal blocks
939 ELSE
940 mult = 1.0_dp
941 END IF
942 btr = factor*mult*sum(a_block*b_block)
943 atrace(iblock_row) = atrace(iblock_row) + btr
944 atrace(iblock_col) = atrace(iblock_col) + btr
945
946 END DO
947 CALL dbcsr_iterator_stop(iter)
948
949 END SUBROUTINE atom_trace
950
951! **************************************************************************************************
952
953END MODULE mulliken
mulliken::ao_charges
Definition mulliken.F:42
mulliken::mulliken_charges
Definition mulliken.F:36
atomic_charges
simple routine to print charges for all atomic charge methods (currently mulliken,...
Definition atomic_charges.F:14
atomic_charges::print_atomic_charges
subroutine, public print_atomic_charges(particle_set, qs_kind_set, scr, title, electronic_charges, atomic_charges)
generates a unified output format for atomic charges
Definition atomic_charges.F:48
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::sp
integer, parameter, public sp
Definition kinds.F:33
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
mulliken
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
Definition mulliken.F:13
mulliken::mulliken_restraint
subroutine, public mulliken_restraint(mulliken_restraint_control, para_env, s_matrix, p_matrix, energy, order_p, ks_matrix, w_matrix)
computes the energy and density matrix derivate of a constraint on the mulliken charges
Definition mulliken.F:73
mulliken::compute_bond_order
subroutine, public compute_bond_order(psmat, spmat, bond_order)
compute Mayer bond orders for a single spin channel for complete result sum up over all spins and mul...
Definition mulliken.F:657
mulliken::compute_charges
subroutine, public compute_charges(p_matrix, s_matrix, charges, para_env)
compute the mulliken charges
Definition mulliken.F:216
mulliken::atom_trace
subroutine, public atom_trace(amat, bmat, factor, atrace)
Compute partial trace of product of two matrices.
Definition mulliken.F:915
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
cp_control_types::mulliken_restraint_type
Definition cp_control_types.F:248
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171