28 USE dbcsr_api,
ONLY: &
29 dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_distribution_type, dbcsr_dot, &
30 dbcsr_get_block_diag, dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, &
31 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
32 dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_replicate_all, &
33 dbcsr_reserve_diag_blocks, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric
57 neighbor_list_iterator_p_type,&
59 neighbor_list_set_p_type,&
65 #include "./base/base_uses.f90"
71 REAL(KIND=
dp),
DIMENSION(:, :),
ALLOCATABLE :: mat
74 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'mao_wfn_analysis'
89 TYPE(qs_environment_type),
POINTER :: qs_env
90 TYPE(section_vals_type),
POINTER :: input_section
91 INTEGER,
INTENT(IN) :: unit_nr
93 CHARACTER(len=*),
PARAMETER :: routinen =
'mao_analysis'
95 CHARACTER(len=2) :: element_symbol, esa, esb, esc
96 INTEGER :: fall, handle, ia, iab, iabc, iatom, ib, ic, icol, ikind, irow, ispin, jatom, &
97 mao_basis, max_iter, me, na, nab, nabc, natom, nb, nc, nimages, nspin, ssize
98 INTEGER,
DIMENSION(:),
POINTER :: col_blk_sizes, mao_blk, mao_blk_sizes, &
99 orb_blk, row_blk_sizes
100 LOGICAL :: analyze_ua, explicit, fo,
for, fos, &
101 found, neglect_abc, print_basis
102 REAL(kind=
dp) :: deltaq, electra(2), eps_ab, eps_abc, eps_filter, eps_fun, eps_grad, epsx, &
103 senabc, senmax, threshold, total_charge, total_spin, ua_charge(2), zeff
104 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: occnuma, occnumabc, qab, qmatab, qmatac, &
105 qmatbc, raq, sab, selnabc, sinv, &
106 smatab, smatac, smatbc, uaq
107 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: occnumab, selnab
108 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block, cmao, diag, qblka, qblkb, qblkc, &
109 rblkl, rblku, sblk, sblka, sblkb, sblkc
110 TYPE(block_type),
ALLOCATABLE,
DIMENSION(:) :: rowblock
111 TYPE(cp_blacs_env_type),
POINTER :: blacs_env
112 TYPE(dbcsr_distribution_type),
POINTER :: dbcsr_dist
113 TYPE(dbcsr_iterator_type) :: dbcsr_iter
114 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: mao_coef, mao_dmat, mao_qmat, mao_smat, &
115 matrix_q, matrix_smm, matrix_smo
116 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks, matrix_p, matrix_s
117 TYPE(dbcsr_type) :: amat, axmat, cgmat, cholmat, crumat, &
118 qmat, qmat_diag, rumat, smat_diag, &
120 TYPE(dft_control_type),
POINTER :: dft_control
121 TYPE(gto_basis_set_p_type),
DIMENSION(:),
POINTER :: mao_basis_set_list, orb_basis_set_list
122 TYPE(kpoint_type),
POINTER :: kpoints
123 TYPE(mp_para_env_type),
POINTER :: para_env
124 TYPE(neighbor_list_iterator_p_type), &
125 DIMENSION(:),
POINTER :: nl_iterator
126 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
127 POINTER :: sab_all, sab_orb, smm_list, smo_list
128 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
129 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
130 TYPE(qs_ks_env_type),
POINTER :: ks_env
131 TYPE(qs_rho_type),
POINTER :: rho
136 IF (.NOT. explicit)
RETURN
138 CALL timeset(routinen, handle)
140 IF (unit_nr > 0)
THEN
141 WRITE (unit_nr,
'(/,T2,A)')
'!-----------------------------------------------------------------------------!'
142 WRITE (unit=unit_nr, fmt=
"(T36,A)")
"MAO ANALYSIS"
143 WRITE (unit=unit_nr, fmt=
"(T12,A)")
"Claus Ehrhardt and Reinhart Ahlrichs, TCA 68:231-245 (1985)"
144 WRITE (unit_nr,
'(T2,A)')
'!-----------------------------------------------------------------------------!'
162 CALL get_qs_env(qs_env, dft_control=dft_control)
163 nimages = dft_control%nimages
164 IF (nimages > 1)
THEN
165 IF (unit_nr > 0)
THEN
166 WRITE (unit=unit_nr, fmt=
"(T2,A)") &
167 "K-Points: MAO's determined and analyzed using Gamma-Point only."
172 NULLIFY (mao_basis_set_list, orb_basis_set_list)
174 unit_nr, print_basis)
177 NULLIFY (smm_list, smo_list)
182 NULLIFY (matrix_smm, matrix_smo)
185 mao_basis_set_list, mao_basis_set_list, smm_list)
187 mao_basis_set_list, orb_basis_set_list, smo_list)
190 CALL get_qs_env(qs_env, rho=rho, matrix_s_kp=matrix_s)
192 nspin =
SIZE(matrix_p, 1)
195 IF (nimages == 1)
THEN
196 CALL mao_build_q(matrix_q, matrix_p, matrix_s, matrix_smm, matrix_smo, smm_list, electra, eps_filter)
198 CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, kpoints=kpoints)
199 CALL mao_build_q(matrix_q, matrix_p, matrix_s, matrix_smm, matrix_smo, smm_list, electra, eps_filter, &
200 nimages=nimages, kpoints=kpoints, matrix_ks=matrix_ks, sab_orb=sab_orb)
208 IF (iatom <= jatom)
THEN
215 CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
216 row=irow, col=icol, block=block, found=found)
217 IF (.NOT. found) fall = fall + 1
221 CALL get_qs_env(qs_env=qs_env, para_env=para_env)
222 CALL para_env%sum(fall)
223 IF (unit_nr > 0 .AND. fall > 0)
THEN
224 WRITE (unit=unit_nr, fmt=
"(/,T2,A,/,T2,A,/)") &
225 "Warning: Extended MAO basis used with original basis filtered density matrix", &
226 "Warning: Possible errors can be controlled with EPS_PGF_ORB"
230 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom)
231 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, dbcsr_dist=dbcsr_dist)
234 ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
236 basis=mao_basis_set_list)
240 IF (col_blk_sizes(iab) < 0) &
241 cpabort(
"Number of MAOs has to be specified in KIND section for all elements")
245 ALLOCATE (mao_coef(ispin)%matrix)
246 CALL dbcsr_create(matrix=mao_coef(ispin)%matrix, &
247 name=
"MAO_COEF", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
248 row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes, nze=0)
249 CALL dbcsr_reserve_diag_blocks(matrix=mao_coef(ispin)%matrix)
251 DEALLOCATE (row_blk_sizes, col_blk_sizes)
255 CALL mao_optimize(mao_coef, matrix_q, matrix_smm, electra, max_iter, eps_grad, epsx, &
260 qs_kind_set, unit_nr, para_env)
263 NULLIFY (mao_dmat, mao_smat, mao_qmat)
267 CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
269 ALLOCATE (mao_dmat(ispin)%matrix)
270 CALL dbcsr_create(mao_dmat(ispin)%matrix, name=
"MAO density", dist=dbcsr_dist, &
271 matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
272 col_blk_size=col_blk_sizes, nze=0)
273 ALLOCATE (mao_smat(ispin)%matrix)
274 CALL dbcsr_create(mao_smat(ispin)%matrix, name=
"MAO overlap", dist=dbcsr_dist, &
275 matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
276 col_blk_size=col_blk_sizes, nze=0)
277 ALLOCATE (mao_qmat(ispin)%matrix)
278 CALL dbcsr_create(mao_qmat(ispin)%matrix, name=
"MAO covar density", dist=dbcsr_dist, &
279 matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
280 col_blk_size=col_blk_sizes, nze=0)
282 CALL dbcsr_create(amat, name=
"MAO overlap", template=mao_dmat(1)%matrix)
283 CALL dbcsr_create(tmat, name=
"MAO Overlap Inverse", template=amat)
284 CALL dbcsr_create(qmat, name=
"MAO covar density", template=amat)
285 CALL dbcsr_create(cgmat, name=
"TEMP matrix", template=mao_coef(1)%matrix)
286 CALL dbcsr_create(axmat, name=
"TEMP", template=amat, matrix_type=dbcsr_type_no_symmetry)
289 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, matrix_smm(1)%matrix, mao_coef(ispin)%matrix, &
291 CALL dbcsr_multiply(
"T",
"N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, amat)
294 CALL invert_hotelling(tmat, amat, threshold, norm_convergence=1.e-4_dp, silent=.true.)
295 CALL dbcsr_copy(mao_smat(ispin)%matrix, amat)
297 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, matrix_q(ispin)%matrix, mao_coef(ispin)%matrix, &
298 0.0_dp, cgmat, filter_eps=eps_filter)
299 CALL dbcsr_multiply(
"T",
"N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, &
300 0.0_dp, qmat, filter_eps=eps_filter)
301 CALL dbcsr_copy(mao_qmat(ispin)%matrix, qmat)
303 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, qmat, tmat, 0.0_dp, axmat, filter_eps=eps_filter)
304 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, tmat, axmat, 0.0_dp, mao_dmat(ispin)%matrix, &
305 filter_eps=eps_filter)
307 CALL dbcsr_release(amat)
308 CALL dbcsr_release(tmat)
309 CALL dbcsr_release(qmat)
310 CALL dbcsr_release(cgmat)
311 CALL dbcsr_release(axmat)
315 CALL dbcsr_dot(mao_dmat(ispin)%matrix, mao_smat(ispin)%matrix, ua_charge(ispin))
316 ua_charge(ispin) = electra(ispin) - ua_charge(ispin)
318 IF (unit_nr > 0)
THEN
321 WRITE (unit=unit_nr, fmt=
"(T2,A,T32,A,i2,T55,A,F12.8)") &
322 "Unassigned charge",
"Spin ", ispin,
"delta charge =", ua_charge(ispin)
330 ALLOCATE (occnuma(natom, nspin))
334 CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, &
335 row=iatom, col=iatom, block=block, found=found)
337 DO iab = 1,
SIZE(block, 1)
338 occnuma(iatom, ispin) = occnuma(iatom, ispin) + block(iab, iab)
343 CALL para_env%sum(occnuma)
346 ALLOCATE (occnumab(natom, natom, nspin))
349 CALL dbcsr_create(qmat_diag, name=
"MAO diagonal density", template=mao_dmat(1)%matrix)
350 CALL dbcsr_create(smat_diag, name=
"MAO diagonal overlap", template=mao_dmat(1)%matrix)
352 CALL dbcsr_get_block_diag(mao_qmat(ispin)%matrix, qmat_diag)
353 CALL dbcsr_replicate_all(qmat_diag)
354 CALL dbcsr_get_block_diag(mao_smat(ispin)%matrix, smat_diag)
355 CALL dbcsr_replicate_all(smat_diag)
357 DO ib = ia + 1, natom
359 CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, &
360 row=ia, col=ib, block=block, found=found)
362 CALL para_env%sum(iab)
364 IF (iab == 0 .AND. para_env%is_source())
THEN
367 occnumab(ia, ib, ispin) = occnuma(ia, ispin) + occnuma(ib, ispin)
368 occnumab(ib, ia, ispin) = occnuma(ia, ispin) + occnuma(ib, ispin)
374 ALLOCATE (sab(nab, nab), qab(nab, nab), sinv(nab, nab))
376 qab(1:na, na + 1:nab) = block(1:na, 1:nb)
377 qab(na + 1:nab, 1:na) = transpose(block(1:na, 1:nb))
378 CALL dbcsr_get_block_p(matrix=qmat_diag, row=ia, col=ia, block=diag, found=fo)
380 qab(1:na, 1:na) = diag(1:na, 1:na)
381 CALL dbcsr_get_block_p(matrix=qmat_diag, row=ib, col=ib, block=diag, found=fo)
383 qab(na + 1:nab, na + 1:nab) = diag(1:nb, 1:nb)
385 CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, &
386 row=ia, col=ib, block=block, found=fo)
388 sab(1:na, na + 1:nab) = block(1:na, 1:nb)
389 sab(na + 1:nab, 1:na) = transpose(block(1:na, 1:nb))
390 CALL dbcsr_get_block_p(matrix=smat_diag, row=ia, col=ia, block=diag, found=fo)
392 sab(1:na, 1:na) = diag(1:na, 1:na)
393 CALL dbcsr_get_block_p(matrix=smat_diag, row=ib, col=ib, block=diag, found=fo)
395 sab(na + 1:nab, na + 1:nab) = diag(1:nb, 1:nb)
397 sinv(1:nab, 1:nab) = sab(1:nab, 1:nab)
400 occnumab(ia, ib, ispin) = sum(qab*sinv)
401 occnumab(ib, ia, ispin) = occnumab(ia, ib, ispin)
403 DEALLOCATE (sab, qab, sinv)
407 CALL dbcsr_release(qmat_diag)
408 CALL dbcsr_release(smat_diag)
410 CALL para_env%sum(occnumab)
413 ALLOCATE (selnab(natom, natom, nspin))
417 DO ib = ia + 1, natom
418 selnab(ia, ib, ispin) = occnuma(ia, ispin) + occnuma(ib, ispin) - occnumab(ia, ib, ispin)
419 selnab(ib, ia, ispin) = selnab(ia, ib, ispin)
424 IF (.NOT. neglect_abc)
THEN
426 nabc = (natom*(natom - 1)*(natom - 2))/6
427 ALLOCATE (occnumabc(nabc, nspin))
430 CALL dbcsr_create(qmat_diag, name=
"MAO diagonal density", template=mao_dmat(1)%matrix)
431 CALL dbcsr_create(smat_diag, name=
"MAO diagonal overlap", template=mao_dmat(1)%matrix)
433 CALL dbcsr_get_block_diag(mao_qmat(ispin)%matrix, qmat_diag)
434 CALL dbcsr_replicate_all(qmat_diag)
435 CALL dbcsr_get_block_diag(mao_smat(ispin)%matrix, smat_diag)
436 CALL dbcsr_replicate_all(smat_diag)
439 CALL dbcsr_get_block_p(matrix=qmat_diag, row=ia, col=ia, block=qblka, found=fo)
441 CALL dbcsr_get_block_p(matrix=smat_diag, row=ia, col=ia, block=sblka, found=fo)
444 DO ib = ia + 1, natom
446 IF (selnab(ia, ib, ispin) < eps_abc)
THEN
447 iabc = iabc + (natom - ib)
450 CALL dbcsr_get_block_p(matrix=qmat_diag, row=ib, col=ib, block=qblkb, found=fo)
452 CALL dbcsr_get_block_p(matrix=smat_diag, row=ib, col=ib, block=sblkb, found=fo)
456 ALLOCATE (qmatab(na, nb), smatab(na, nb))
457 CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, row=ia, col=ib, &
458 block=block, found=found)
460 IF (found) qmatab(1:na, 1:nb) = block(1:na, 1:nb)
461 CALL para_env%sum(qmatab)
462 CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, row=ia, col=ib, &
463 block=block, found=found)
465 IF (found) smatab(1:na, 1:nb) = block(1:na, 1:nb)
466 CALL para_env%sum(smatab)
467 DO ic = ib + 1, natom
469 IF ((selnab(ia, ic, ispin) < eps_abc) .OR. (selnab(ib, ic, ispin) < eps_abc))
THEN
473 CALL dbcsr_get_block_p(matrix=qmat_diag, row=ic, col=ic, block=qblkc, found=fo)
475 CALL dbcsr_get_block_p(matrix=smat_diag, row=ic, col=ic, block=sblkc, found=fo)
478 ALLOCATE (qmatac(na, nc), smatac(na, nc))
479 CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, row=ia, col=ic, &
480 block=block, found=found)
482 IF (found) qmatac(1:na, 1:nc) = block(1:na, 1:nc)
483 CALL para_env%sum(qmatac)
484 CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, row=ia, col=ic, &
485 block=block, found=found)
487 IF (found) smatac(1:na, 1:nc) = block(1:na, 1:nc)
488 CALL para_env%sum(smatac)
489 ALLOCATE (qmatbc(nb, nc), smatbc(nb, nc))
490 CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, row=ib, col=ic, &
491 block=block, found=found)
493 IF (found) qmatbc(1:nb, 1:nc) = block(1:nb, 1:nc)
494 CALL para_env%sum(qmatbc)
495 CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, row=ib, col=ic, &
496 block=block, found=found)
498 IF (found) smatbc(1:nb, 1:nc) = block(1:nb, 1:nc)
499 CALL para_env%sum(smatbc)
502 ALLOCATE (sab(nabc, nabc), sinv(nabc, nabc), qab(nabc, nabc))
504 qab(1:na, 1:na) = qblka(1:na, 1:na)
505 qab(na + 1:nab, na + 1:nab) = qblkb(1:nb, 1:nb)
506 qab(nab + 1:nabc, nab + 1:nabc) = qblkc(1:nc, 1:nc)
507 qab(1:na, na + 1:nab) = qmatab(1:na, 1:nb)
508 qab(na + 1:nab, 1:na) = transpose(qmatab(1:na, 1:nb))
509 qab(1:na, nab + 1:nabc) = qmatac(1:na, 1:nc)
510 qab(nab + 1:nabc, 1:na) = transpose(qmatac(1:na, 1:nc))
511 qab(na + 1:nab, nab + 1:nabc) = qmatbc(1:nb, 1:nc)
512 qab(nab + 1:nabc, na + 1:nab) = transpose(qmatbc(1:nb, 1:nc))
514 sab(1:na, 1:na) = sblka(1:na, 1:na)
515 sab(na + 1:nab, na + 1:nab) = sblkb(1:nb, 1:nb)
516 sab(nab + 1:nabc, nab + 1:nabc) = sblkc(1:nc, 1:nc)
517 sab(1:na, na + 1:nab) = smatab(1:na, 1:nb)
518 sab(na + 1:nab, 1:na) = transpose(smatab(1:na, 1:nb))
519 sab(1:na, nab + 1:nabc) = smatac(1:na, 1:nc)
520 sab(nab + 1:nabc, 1:na) = transpose(smatac(1:na, 1:nc))
521 sab(na + 1:nab, nab + 1:nabc) = smatbc(1:nb, 1:nc)
522 sab(nab + 1:nabc, na + 1:nab) = transpose(smatbc(1:nb, 1:nc))
524 sinv(1:nabc, 1:nabc) = sab(1:nabc, 1:nabc)
528 me = mod(iabc, para_env%num_pe)
529 IF (me == para_env%mepos)
THEN
530 occnumabc(iabc, ispin) = sum(qab*sinv)
532 occnumabc(iabc, ispin) = 0.0_dp
535 DEALLOCATE (sab, sinv, qab)
536 DEALLOCATE (qmatac, smatac)
537 DEALLOCATE (qmatbc, smatbc)
539 DEALLOCATE (qmatab, smatab)
542 CALL dbcsr_release(qmat_diag)
543 CALL dbcsr_release(smat_diag)
545 CALL para_env%sum(occnumabc)
548 IF (.NOT. neglect_abc)
THEN
550 nabc = (natom*(natom - 1)*(natom - 2))/6
551 ALLOCATE (selnabc(nabc, nspin))
556 DO ib = ia + 1, natom
557 DO ic = ib + 1, natom
559 IF (occnumabc(iabc, ispin) >= 0.0_dp)
THEN
560 selnabc(iabc, ispin) = occnuma(ia, ispin) + occnuma(ib, ispin) + occnuma(ic, ispin) - &
561 occnumab(ia, ib, ispin) - occnumab(ia, ic, ispin) - occnumab(ib, ic, ispin) + &
562 occnumabc(iabc, ispin)
571 ALLOCATE (raq(natom, nspin))
575 raq(ia, ispin) = occnuma(ia, ispin)
577 raq(ia, ispin) = raq(ia, ispin) - 0.5_dp*selnab(ia, ib, ispin)
580 IF (.NOT. neglect_abc)
THEN
583 DO ib = ia + 1, natom
584 DO ic = ib + 1, natom
586 raq(ia, ispin) = raq(ia, ispin) + selnabc(iabc, ispin)/3._dp
587 raq(ib, ispin) = raq(ib, ispin) + selnabc(iabc, ispin)/3._dp
588 raq(ic, ispin) = raq(ic, ispin) + selnabc(iabc, ispin)/3._dp
597 deltaq = (electra(ispin) - sum(raq(1:natom, ispin))) - ua_charge(ispin)
598 IF (unit_nr > 0)
THEN
599 WRITE (unit=unit_nr, fmt=
"(T2,A,T32,A,i2,T55,A,F12.8)") &
600 "Cutoff error on charge",
"Spin ", ispin,
"error charge =", deltaq
605 ALLOCATE (uaq(natom, nspin))
608 CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env)
609 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb, sab_all=sab_all)
610 CALL dbcsr_get_info(mao_coef(1)%matrix, row_blk_size=mao_blk_sizes, &
611 col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
612 CALL dbcsr_get_info(matrix_s(1, 1)%matrix, row_blk_size=row_blk_sizes)
613 CALL dbcsr_create(amat, name=
"temp", template=matrix_s(1, 1)%matrix)
614 CALL dbcsr_create(tmat, name=
"temp", template=mao_coef(1)%matrix)
616 CALL dbcsr_get_block_diag(matrix_smm(1)%matrix, smat_diag)
617 CALL dbcsr_replicate_all(smat_diag)
619 ALLOCATE (orb_blk(natom), mao_blk(natom))
621 orb_blk = row_blk_sizes
622 mao_blk = row_blk_sizes
623 mao_blk(ia) = col_blk_sizes(ia)
624 CALL dbcsr_create(sumat, name=
"Smat", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
625 row_blk_size=mao_blk, col_blk_size=mao_blk, nze=0)
627 CALL dbcsr_create(cholmat, name=
"Cholesky matrix", dist=dbcsr_dist, &
628 matrix_type=dbcsr_type_no_symmetry, row_blk_size=mao_blk, col_blk_size=mao_blk, nze=0)
629 CALL dbcsr_create(rumat, name=
"Rmat", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
630 row_blk_size=orb_blk, col_blk_size=mao_blk, nze=0)
632 CALL dbcsr_create(crumat, name=
"Rmat*Umat", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
633 row_blk_size=orb_blk, col_blk_size=mao_blk, nze=0)
635 ALLOCATE (rowblock(natom))
637 na = mao_blk_sizes(ia)
638 nb = row_blk_sizes(ib)
639 ALLOCATE (rowblock(ib)%mat(na, nb))
640 rowblock(ib)%mat = 0.0_dp
641 CALL dbcsr_get_block_p(matrix=matrix_smo(1)%matrix, row=ia, col=ib, &
642 block=block, found=found)
643 IF (found) rowblock(ib)%mat(1:na, 1:nb) = block(1:na, 1:nb)
644 CALL para_env%sum(rowblock(ib)%mat)
648 CALL dbcsr_copy(tmat, mao_coef(ispin)%matrix)
649 CALL dbcsr_replicate_all(tmat)
650 CALL dbcsr_iterator_start(dbcsr_iter, matrix_s(1, 1)%matrix)
651 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
652 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, block)
653 CALL dbcsr_get_block_p(matrix=sumat, row=iatom, col=jatom, block=sblk, found=fos)
655 CALL dbcsr_get_block_p(matrix=rumat, row=iatom, col=jatom, block=rblku, found=
for)
657 CALL dbcsr_get_block_p(matrix=rumat, row=jatom, col=iatom, block=rblkl, found=
for)
659 CALL dbcsr_get_block_p(matrix=tmat, row=ia, col=ia, block=cmao, found=found)
661 IF (iatom /= ia .AND. jatom /= ia)
THEN
665 rblkl = transpose(block)
666 ELSE IF (iatom /= ia)
THEN
667 rblkl = transpose(block)
668 sblk = matmul(transpose(rowblock(iatom)%mat), cmao)
670 ELSE IF (jatom /= ia)
THEN
672 sblk = matmul(transpose(cmao), rowblock(jatom)%mat)
673 rblkl = transpose(sblk)
675 CALL dbcsr_get_block_p(matrix=smat_diag, row=ia, col=ia, block=block, found=found)
677 sblk = matmul(transpose(cmao), matmul(block, cmao))
678 rblku = matmul(transpose(rowblock(ia)%mat), cmao)
681 CALL dbcsr_iterator_stop(dbcsr_iter)
683 CALL dbcsr_desymmetrize(sumat, cholmat)
688 transa=
"N", para_env=para_env, blacs_env=blacs_env)
690 CALL dbcsr_multiply(
"N",
"T", 1.0_dp, crumat, crumat, 0.0_dp, amat, &
691 filter_eps=eps_filter)
693 CALL dbcsr_dot(matrix_p(ispin, 1)%matrix, amat, uaq(ia, ispin))
694 uaq(ia, ispin) = uaq(ia, ispin) - electra(ispin)
697 CALL dbcsr_release(sumat)
698 CALL dbcsr_release(cholmat)
699 CALL dbcsr_release(rumat)
700 CALL dbcsr_release(crumat)
703 DEALLOCATE (rowblock(ib)%mat)
705 DEALLOCATE (rowblock)
707 CALL dbcsr_release(smat_diag)
708 CALL dbcsr_release(amat)
709 CALL dbcsr_release(tmat)
710 DEALLOCATE (orb_blk, mao_blk)
713 raq(1:natom, 1:nspin) = raq(1:natom, 1:nspin) - uaq(1:natom, 1:nspin)
715 deltaq = electra(ispin) - sum(raq(1:natom, ispin))
716 IF (unit_nr > 0)
THEN
717 WRITE (unit=unit_nr, fmt=
"(T2,A,T32,A,i2,T55,A,F12.8)") &
718 "Charge/Atom redistributed",
"Spin ", ispin,
"delta charge =", &
719 (deltaq + ua_charge(ispin))/real(natom, kind=
dp)
724 IF (unit_nr > 0)
THEN
726 WRITE (unit_nr,
"(/,T2,A,T40,A,T75,A)")
"MAO atomic charges ",
"Atom",
"Charge"
728 WRITE (unit_nr,
"(/,T2,A,T40,A,T55,A,T70,A)")
"MAO atomic charges ",
"Atom",
"Charge",
"Spin Charge"
731 deltaq = electra(ispin) - sum(raq(1:natom, ispin))
732 raq(:, ispin) = raq(:, ispin) + deltaq/real(natom, kind=
dp)
734 total_charge = 0.0_dp
738 element_symbol=element_symbol, kind_number=ikind)
741 WRITE (unit_nr,
"(T30,I6,T42,A2,T69,F12.6)") iatom, element_symbol, zeff - raq(iatom, 1)
742 total_charge = total_charge + (zeff - raq(iatom, 1))
744 WRITE (unit_nr,
"(T30,I6,T42,A2,T48,F12.6,T69,F12.6)") iatom, element_symbol, &
745 zeff - raq(iatom, 1) - raq(iatom, 2), raq(iatom, 1) - raq(iatom, 2)
746 total_charge = total_charge + (zeff - raq(iatom, 1) - raq(iatom, 2))
747 total_spin = total_spin + (raq(iatom, 1) - raq(iatom, 2))
751 WRITE (unit_nr,
"(T2,A,T69,F12.6)")
"Total Charge", total_charge
753 WRITE (unit_nr,
"(T2,A,T49,F12.6,T69,F12.6)")
"Total Charge", total_charge, total_spin
759 IF (unit_nr > 0)
THEN
761 WRITE (unit_nr,
"(/,T2,A,T40,A,T75,A)")
"MAO hypervalent charges ",
"Atom",
"Charge"
763 WRITE (unit_nr,
"(/,T2,A,T40,A,T55,A,T70,A)")
"MAO hypervalent charges ",
"Atom", &
764 "Charge",
"Spin Charge"
766 total_charge = 0.0_dp
770 element_symbol=element_symbol)
772 WRITE (unit_nr,
"(T30,I6,T42,A2,T69,F12.6)") iatom, element_symbol, uaq(iatom, 1)
773 total_charge = total_charge + uaq(iatom, 1)
775 WRITE (unit_nr,
"(T30,I6,T42,A2,T48,F12.6,T69,F12.6)") iatom, element_symbol, &
776 uaq(iatom, 1) + uaq(iatom, 2), uaq(iatom, 1) - uaq(iatom, 2)
777 total_charge = total_charge + uaq(iatom, 1) + uaq(iatom, 2)
778 total_spin = total_spin + uaq(iatom, 1) - uaq(iatom, 2)
782 WRITE (unit_nr,
"(T2,A,T69,F12.6)")
"Total Charge", total_charge
784 WRITE (unit_nr,
"(T2,A,T49,F12.6,T69,F12.6)")
"Total Charge", total_charge, total_spin
790 IF (unit_nr > 0)
THEN
792 WRITE (unit_nr,
"(/,T2,A,T31,A,T40,A,T78,A)")
"Shared electron numbers ",
"Atom",
"Atom",
"SEN"
794 WRITE (unit_nr,
"(/,T2,A,T31,A,T40,A,T51,A,T63,A,T71,A)")
"Shared electron numbers ",
"Atom",
"Atom", &
795 "SEN(1)",
"SEN(2)",
"SEN(total)"
798 DO ib = ia + 1, natom
799 CALL get_atomic_kind(atomic_kind=particle_set(ia)%atomic_kind, element_symbol=esa)
800 CALL get_atomic_kind(atomic_kind=particle_set(ib)%atomic_kind, element_symbol=esb)
802 IF (selnab(ia, ib, 1) > eps_ab)
THEN
803 WRITE (unit_nr,
"(T26,I6,' ',A2,T35,I6,' ',A2,T69,F12.6)") ia, esa, ib, esb, selnab(ia, ib, 1)
806 IF ((selnab(ia, ib, 1) + selnab(ia, ib, 2)) > eps_ab)
THEN
807 WRITE (unit_nr,
"(T26,I6,' ',A2,T35,I6,' ',A2,T45,3F12.6)") ia, esa, ib, esb, &
808 selnab(ia, ib, 1), selnab(ia, ib, 2), (selnab(ia, ib, 1) + selnab(ia, ib, 2))
815 IF (.NOT. neglect_abc)
THEN
817 IF (unit_nr > 0)
THEN
818 WRITE (unit_nr,
"(/,T2,A,T40,A,T49,A,T58,A,T78,A)")
"Shared electron numbers ABC", &
819 "Atom",
"Atom",
"Atom",
"SEN"
823 DO ib = ia + 1, natom
824 DO ic = ib + 1, natom
826 senabc = sum(selnabc(iabc, :))
827 senmax = max(senmax, senabc)
828 IF (senabc > eps_abc)
THEN
829 CALL get_atomic_kind(atomic_kind=particle_set(ia)%atomic_kind, element_symbol=esa)
830 CALL get_atomic_kind(atomic_kind=particle_set(ib)%atomic_kind, element_symbol=esb)
831 CALL get_atomic_kind(atomic_kind=particle_set(ic)%atomic_kind, element_symbol=esc)
832 WRITE (unit_nr,
"(T35,I6,' ',A2,T44,I6,' ',A2,T53,I6,' ',A2,T69,F12.6)") &
833 ia, esa, ib, esb, ic, esc, senabc
838 WRITE (unit_nr,
"(T2,A,T69,F12.6)")
"Maximum SEN value calculated", senmax
842 IF (unit_nr > 0)
THEN
843 WRITE (unit_nr,
'(/,T2,A)') &
844 '!---------------------------END OF MAO ANALYSIS-------------------------------!'
848 DEALLOCATE (occnuma, occnumab, selnab, raq, uaq)
849 IF (.NOT. neglect_abc)
THEN
850 DEALLOCATE (occnumabc, selnabc)
857 DEALLOCATE (mao_basis_set_list, orb_basis_set_list)
868 CALL timestop(handle)
for(int lxp=0;lxp<=lp;lxp++)
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public ehrhardt1985
integer, save, public heinzmann1976
methods related to the blacs parallel environment
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Interface to (sca)lapack for the Cholesky based procedures.
subroutine, public cp_dbcsr_cholesky_decompose(matrix, n, para_env, blacs_env)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
subroutine, public cp_dbcsr_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa, para_env, blacs_env)
...
Routines that link DBCSR and CP2K concepts together.
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
DBCSR operations in CP2K.
Routines useful for iterative matrix calculations.
subroutine, public invert_hotelling(matrix_inverse, matrix, threshold, use_inv_as_guess, norm_convergence, filter_eps, accelerator_order, max_iter_lanczos, eps_lanczos, silent)
invert a symmetric positive definite matrix by Hotelling's method explicit symmetrization makes this ...
Defines the basic variable types.
integer, parameter, public dp
Types and basic routines needed for a kpoint calculation.
Calculate MAO's and analyze wavefunctions.
Calculate MAO's and analyze wavefunctions.
subroutine, public mao_reference_basis(qs_env, mao_basis, mao_basis_set_list, orb_basis_set_list, iunit, print_basis)
Define the MAO reference basis set.
subroutine, public mao_basis_analysis(mao_coef, matrix_smm, mao_basis_set_list, particle_set, qs_kind_set, unit_nr, para_env)
Analyze the MAO basis, projection on angular functions.
subroutine, public mao_build_q(matrix_q, matrix_p, matrix_s, matrix_smm, matrix_smo, smm_list, electra, eps_filter, nimages, kpoints, matrix_ks, sab_orb)
Calculte the Q=APA(T) matrix, A=(MAO,ORB) overlap.
Calculate MAO's and analyze wavefunctions.
subroutine, public mao_optimize(mao_coef, matrix_q, matrix_smm, electra, max_iter, eps_grad, eps1, iolevel, iw)
...
Calculate MAO's and analyze wavefunctions.
subroutine, public mao_analysis(qs_env, input_section, unit_nr)
...
Collection of simple mathematical functions and subroutines.
subroutine, public invmat_symm(a, cholesky_triangle)
returns inverse of real symmetric, positive definite matrix
Interface to the message passing library MPI.
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Define the data structure for the particle information.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Define the neighbor list data types and the corresponding functionality.
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Generate the atomic neighbor lists.
subroutine, public setup_neighbor_list(ab_list, basis_set_a, basis_set_b, qs_env, mic, symmetric, molecular, operator_type)
Build a neighborlist.
Calculation of overlap matrix, its derivatives and forces.
subroutine, public build_overlap_matrix_simple(ks_env, matrix_s, basis_set_list_a, basis_set_list_b, sab_nl)
Calculation of the overlap matrix over Cartesian Gaussian functions.
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...