24 USE dbcsr_api,
ONLY: &
25 dbcsr_add, dbcsr_binary_read, dbcsr_checksum, dbcsr_complete_redistribute, &
26 dbcsr_copy_into_existing, dbcsr_create, dbcsr_desymmetrize, dbcsr_distribution_get, &
27 dbcsr_distribution_new, dbcsr_distribution_type, dbcsr_dot, dbcsr_filter, &
28 dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
29 dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
30 dbcsr_release, dbcsr_reserve_diag_blocks, dbcsr_scale, dbcsr_set, dbcsr_type
67 #include "./base/base_uses.f90"
73 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'pao_methods'
92 TYPE(pao_env_type),
POINTER :: pao
93 TYPE(qs_environment_type),
POINTER :: qs_env
95 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_init_kinds'
97 INTEGER :: handle, i, ikind, pao_basis_size
98 TYPE(gto_basis_set_type),
POINTER :: basis_set
99 TYPE(pao_descriptor_type),
DIMENSION(:),
POINTER :: pao_descriptors
101 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
103 CALL timeset(routinen, handle)
104 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
106 DO ikind = 1,
SIZE(qs_kind_set)
108 basis_set=basis_set, &
109 pao_basis_size=pao_basis_size, &
111 pao_descriptors=pao_descriptors)
113 IF (pao_basis_size < 1)
THEN
115 CALL set_qs_kind(qs_kind_set(ikind), pao_basis_size=basis_set%nsgf)
122 DO i = 1,
SIZE(pao_descriptors)
123 pao_descriptors(i)%beta_radius =
exp_radius(0, pao_descriptors(i)%beta, pao%eps_pgf, 1.0_dp)
124 pao_descriptors(i)%screening_radius =
exp_radius(0, pao_descriptors(i)%screening, pao%eps_pgf, 1.0_dp)
127 CALL timestop(handle)
135 TYPE(pao_env_type),
POINTER :: pao
137 INTEGER :: iatom, natoms
138 INTEGER,
DIMENSION(:),
POINTER :: pao_basis, param_cols, param_rows, &
141 CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=pri_basis, col_blk_size=pao_basis)
142 cpassert(
SIZE(pao_basis) ==
SIZE(pri_basis))
143 natoms =
SIZE(pao_basis)
145 CALL dbcsr_get_info(pao%matrix_X, row_blk_size=param_rows, col_blk_size=param_cols)
146 cpassert(
SIZE(param_rows) == natoms .AND.
SIZE(param_cols) == natoms)
148 IF (pao%iw_atoms > 0)
THEN
150 WRITE (pao%iw_atoms,
"(A,I7,T20,A,I3,T45,A,I3,T65,A,I3)") &
151 " PAO| atom: ", iatom, &
152 " prim_basis: ", pri_basis(iatom), &
153 " pao_basis: ", pao_basis(iatom), &
154 " pao_params: ", (param_cols(iatom)*param_rows(iatom))
165 TYPE(pao_env_type),
POINTER :: pao
166 TYPE(qs_environment_type),
POINTER :: qs_env
168 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_build_orthogonalizer'
170 INTEGER :: acol, arow, handle, i, iatom, j, k, n
173 REAL(
dp),
DIMENSION(:),
POINTER :: evals
174 REAL(
dp),
DIMENSION(:, :),
POINTER :: a, block_n, block_n_inv, block_s
175 TYPE(dbcsr_iterator_type) :: iter
176 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s
178 CALL timeset(routinen, handle)
182 CALL dbcsr_create(pao%matrix_N, template=matrix_s(1)%matrix, name=
"PAO matrix_N")
183 CALL dbcsr_reserve_diag_blocks(pao%matrix_N)
185 CALL dbcsr_create(pao%matrix_N_inv, template=matrix_s(1)%matrix, name=
"PAO matrix_N_inv")
186 CALL dbcsr_reserve_diag_blocks(pao%matrix_N_inv)
190 CALL dbcsr_iterator_start(iter, pao%matrix_N)
191 DO WHILE (dbcsr_iterator_blocks_left(iter))
192 CALL dbcsr_iterator_next_block(iter, arow, acol, block_n)
193 iatom = arow; cpassert(arow == acol)
195 CALL dbcsr_get_block_p(matrix=pao%matrix_N_inv, row=iatom, col=iatom, block=block_n_inv, found=found)
196 cpassert(
ASSOCIATED(block_n_inv))
198 CALL dbcsr_get_block_p(matrix=matrix_s(1)%matrix, row=iatom, col=iatom, block=block_s, found=found)
199 cpassert(
ASSOCIATED(block_s))
201 n =
SIZE(block_s, 1); cpassert(
SIZE(block_s, 1) ==
SIZE(block_s, 2))
202 ALLOCATE (a(n, n), evals(n))
212 w = 1.0_dp/sqrt(evals(k))
216 block_n(i, j) = block_n(i, j) + w*a(i, k)*a(j, k)
217 block_n_inv(i, j) = block_n_inv(i, j) + v*a(i, k)*a(j, k)
221 DEALLOCATE (a, evals)
223 CALL dbcsr_iterator_stop(iter)
227 CALL dbcsr_create(pao%matrix_N_diag, &
228 name=
"PAO matrix_N_diag", &
229 dist=pao%diag_distribution, &
230 template=matrix_s(1)%matrix)
231 CALL dbcsr_reserve_diag_blocks(pao%matrix_N_diag)
232 CALL dbcsr_complete_redistribute(pao%matrix_N, pao%matrix_N_diag)
233 CALL dbcsr_create(pao%matrix_N_inv_diag, &
234 name=
"PAO matrix_N_inv_diag", &
235 dist=pao%diag_distribution, &
236 template=matrix_s(1)%matrix)
237 CALL dbcsr_reserve_diag_blocks(pao%matrix_N_inv_diag)
238 CALL dbcsr_complete_redistribute(pao%matrix_N_inv, pao%matrix_N_inv_diag)
240 CALL timestop(handle)
249 TYPE(pao_env_type),
POINTER :: pao
250 TYPE(qs_environment_type),
POINTER :: qs_env
252 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_build_selector'
254 INTEGER :: acol, arow, handle, i, iatom, ikind, m, &
256 INTEGER,
DIMENSION(:),
POINTER :: blk_sizes_aux, blk_sizes_pri
257 REAL(
dp),
DIMENSION(:, :),
POINTER :: block_y
258 TYPE(dbcsr_iterator_type) :: iter
259 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s
260 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
261 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
263 CALL timeset(routinen, handle)
268 qs_kind_set=qs_kind_set, &
269 particle_set=particle_set)
271 CALL dbcsr_get_info(matrix_s(1)%matrix, col_blk_size=blk_sizes_pri)
273 ALLOCATE (blk_sizes_aux(natoms))
275 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
276 CALL get_qs_kind(qs_kind_set(ikind), pao_basis_size=m)
278 IF (blk_sizes_pri(iatom) < m) &
279 cpabort(
"PAO basis size exceeds primary basis size.")
280 blk_sizes_aux(iatom) = m
283 CALL dbcsr_create(pao%matrix_Y, &
284 template=matrix_s(1)%matrix, &
286 row_blk_size=blk_sizes_pri, &
287 col_blk_size=blk_sizes_aux, &
289 DEALLOCATE (blk_sizes_aux)
291 CALL dbcsr_reserve_diag_blocks(pao%matrix_Y)
295 CALL dbcsr_iterator_start(iter, pao%matrix_Y)
296 DO WHILE (dbcsr_iterator_blocks_left(iter))
297 CALL dbcsr_iterator_next_block(iter, arow, acol, block_y)
301 block_y(i, i) = 1.0_dp
304 CALL dbcsr_iterator_stop(iter)
307 CALL timestop(handle)
316 TYPE(pao_env_type),
POINTER :: pao
317 TYPE(qs_environment_type),
POINTER :: qs_env
319 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_build_diag_distribution'
321 INTEGER :: handle, iatom, natoms, pgrid_cols, &
323 INTEGER,
DIMENSION(:),
POINTER :: diag_col_dist, diag_row_dist
324 TYPE(dbcsr_distribution_type) :: main_dist
325 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s
327 CALL timeset(routinen, handle)
329 CALL get_qs_env(qs_env, natom=natoms, matrix_s=matrix_s)
332 CALL dbcsr_get_info(matrix=matrix_s(1)%matrix, distribution=main_dist)
333 CALL dbcsr_distribution_get(main_dist, nprows=pgrid_rows, npcols=pgrid_cols)
336 ALLOCATE (diag_row_dist(natoms), diag_col_dist(natoms))
338 diag_row_dist(iatom) = mod(iatom - 1, pgrid_rows)
339 diag_col_dist(iatom) = mod((iatom - 1)/pgrid_rows, pgrid_cols)
343 CALL dbcsr_distribution_new(pao%diag_distribution, template=main_dist, &
344 row_dist=diag_row_dist, col_dist=diag_col_dist)
346 DEALLOCATE (diag_row_dist, diag_col_dist)
348 CALL timestop(handle)
357 TYPE(pao_env_type),
POINTER :: pao
358 TYPE(qs_environment_type),
POINTER :: qs_env
360 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_build_matrix_X'
362 INTEGER :: handle, iatom, ikind, natoms
363 INTEGER,
DIMENSION(:),
POINTER :: col_blk_size, row_blk_size
364 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
366 CALL timeset(routinen, handle)
370 particle_set=particle_set)
373 ALLOCATE (row_blk_size(natoms), col_blk_size(natoms))
376 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
381 CALL dbcsr_create(pao%matrix_X, &
382 name=
"PAO matrix_X", &
383 dist=pao%diag_distribution, &
385 row_blk_size=row_blk_size, &
386 col_blk_size=col_blk_size)
387 DEALLOCATE (row_blk_size, col_blk_size)
389 CALL dbcsr_reserve_diag_blocks(pao%matrix_X)
390 CALL dbcsr_set(pao%matrix_X, 0.0_dp)
392 CALL timestop(handle)
401 TYPE(pao_env_type),
POINTER :: pao
402 TYPE(qs_environment_type),
POINTER :: qs_env
404 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
405 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s
406 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
410 atomic_kind_set=atomic_kind_set, &
411 qs_kind_set=qs_kind_set)
414 CALL dbcsr_create(pao%matrix_H0, &
415 name=
"PAO matrix_H0", &
416 dist=pao%diag_distribution, &
417 template=matrix_s(1)%matrix)
418 CALL dbcsr_reserve_diag_blocks(pao%matrix_H0)
438 TYPE(pao_env_type),
POINTER :: pao
439 TYPE(ls_scf_env_type) :: ls_scf_env
440 REAL(kind=
dp),
INTENT(IN) :: new_energy
441 LOGICAL,
INTENT(OUT) :: is_converged
443 REAL(kind=
dp) :: energy_diff, loop_eps, now, time_diff
446 energy_diff = new_energy - pao%energy_prev
447 pao%energy_prev = new_energy
449 time_diff = now - pao%step_start_time
450 pao%step_start_time = now
453 loop_eps = pao%norm_G/ls_scf_env%nelectron_total
454 is_converged = loop_eps < pao%eps_pao
456 IF (pao%istep > 1)
THEN
457 IF (pao%iw > 0)
WRITE (pao%iw, *)
"PAO| energy improvement:", energy_diff
461 IF (pao%iw > 0)
WRITE (pao%iw,
'(A,I6,11X,F20.9,1X,E10.3,1X,E10.3,1X,F9.3)') &
466 pao%linesearch%step_size, &
479 TYPE(pao_env_type),
POINTER :: pao
480 TYPE(qs_environment_type),
POINTER :: qs_env
481 TYPE(ls_scf_env_type),
TARGET :: ls_scf_env
482 REAL(kind=
dp),
INTENT(OUT) :: energy
484 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_calc_energy'
486 INTEGER :: handle, ispin
487 REAL(kind=
dp) :: penalty, trace_ph
489 CALL timeset(routinen, handle)
492 CALL pao_calc_ab(pao, qs_env, ls_scf_env, gradient=.false., penalty=penalty)
495 CALL pao_rebuild_s(qs_env, ls_scf_env)
498 CALL pao_dm_trs4(qs_env, ls_scf_env)
502 DO ispin = 1, ls_scf_env%nspins
503 CALL dbcsr_dot(ls_scf_env%matrix_p(ispin), ls_scf_env%matrix_ks(ispin), trace_ph)
504 energy = energy + trace_ph
508 energy = energy + penalty
512 WRITE (pao%iw, *)
"PAO| energy:", energy,
"penalty:", penalty
514 CALL timestop(handle)
522 TYPE(ls_scf_env_type) :: ls_scf_env
524 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_check_trace_PS'
526 INTEGER :: handle, ispin
527 REAL(kind=
dp) :: tmp, trace_ps
528 TYPE(dbcsr_type) :: matrix_s_desym
530 CALL timeset(routinen, handle)
531 CALL dbcsr_create(matrix_s_desym, template=ls_scf_env%matrix_s, matrix_type=
"N")
532 CALL dbcsr_desymmetrize(ls_scf_env%matrix_s, matrix_s_desym)
535 DO ispin = 1, ls_scf_env%nspins
536 CALL dbcsr_dot(ls_scf_env%matrix_p(ispin), matrix_s_desym, tmp)
537 trace_ps = trace_ps + tmp
540 CALL dbcsr_release(matrix_s_desym)
542 IF (abs(ls_scf_env%nelectron_total - trace_ps) > 0.5) &
543 cpabort(
"Number of electrons wrong. Trace(PS) ="//cp_to_string(trace_ps))
545 CALL timestop(handle)
553 SUBROUTINE pao_read_preopt_dm(pao, qs_env)
554 TYPE(pao_env_type),
POINTER :: pao
555 TYPE(qs_environment_type),
POINTER :: qs_env
557 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_read_preopt_dm'
559 INTEGER :: handle, ispin
560 REAL(kind=
dp) :: cs_pos
561 TYPE(dbcsr_distribution_type) :: dist
562 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s, rho_ao
563 TYPE(dbcsr_type) :: matrix_tmp
564 TYPE(dft_control_type),
POINTER :: dft_control
565 TYPE(qs_energy_type),
POINTER :: energy
566 TYPE(qs_rho_type),
POINTER :: rho
568 CALL timeset(routinen, handle)
571 dft_control=dft_control, &
578 IF (dft_control%nspins /= 1) cpabort(
"open shell not yet implemented")
580 CALL dbcsr_get_info(matrix_s(1)%matrix, distribution=dist)
582 DO ispin = 1, dft_control%nspins
583 CALL dbcsr_binary_read(pao%preopt_dm_file, matrix_new=matrix_tmp, distribution=dist)
584 cs_pos = dbcsr_checksum(matrix_tmp, pos=.true.)
585 IF (pao%iw > 0)
WRITE (pao%iw, *)
"PAO| Read restart DM "// &
586 trim(pao%preopt_dm_file)//
" with checksum: ", cs_pos
587 CALL dbcsr_copy_into_existing(rho_ao(ispin)%matrix, matrix_tmp)
588 CALL dbcsr_release(matrix_tmp)
595 just_energy=.false., print_active=.true.)
596 IF (pao%iw > 0)
WRITE (pao%iw, *)
"PAO| Quickstep energy from restart density:", energy%total
598 CALL timestop(handle)
600 END SUBROUTINE pao_read_preopt_dm
607 SUBROUTINE pao_rebuild_s(qs_env, ls_scf_env)
608 TYPE(qs_environment_type),
POINTER :: qs_env
609 TYPE(ls_scf_env_type),
TARGET :: ls_scf_env
611 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_rebuild_S'
614 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s
616 CALL timeset(routinen, handle)
618 CALL dbcsr_release(ls_scf_env%matrix_s_inv)
619 CALL dbcsr_release(ls_scf_env%matrix_s_sqrt)
620 CALL dbcsr_release(ls_scf_env%matrix_s_sqrt_inv)
625 CALL timestop(handle)
626 END SUBROUTINE pao_rebuild_s
633 SUBROUTINE pao_dm_trs4(qs_env, ls_scf_env)
634 TYPE(qs_environment_type),
POINTER :: qs_env
635 TYPE(ls_scf_env_type),
TARGET :: ls_scf_env
637 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_dm_trs4'
639 CHARACTER(LEN=default_path_length) :: project_name
640 INTEGER :: handle, ispin, nelectron_spin_real, nspin
642 REAL(kind=
dp) :: homo_spin, lumo_spin, mu_spin
643 TYPE(cp_logger_type),
POINTER :: logger
644 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_ks
646 CALL timeset(routinen, handle)
648 project_name = logger%iter_info%project_name
649 nspin = ls_scf_env%nspins
653 CALL matrix_qs_to_ls(ls_scf_env%matrix_ks(ispin), matrix_ks(ispin)%matrix, &
654 ls_scf_env%ls_mstruct, covariant=.true.)
656 nelectron_spin_real = ls_scf_env%nelectron_spin(ispin)
657 IF (ls_scf_env%nspins == 1) nelectron_spin_real = nelectron_spin_real/2
659 ls_scf_env%matrix_s_sqrt_inv, &
660 nelectron_spin_real, ls_scf_env%eps_filter, homo_spin, lumo_spin, mu_spin, &
661 dynamic_threshold=.false., converged=converged, &
662 max_iter_lanczos=ls_scf_env%max_iter_lanczos, &
663 eps_lanczos=ls_scf_env%eps_lanczos)
664 IF (.NOT. converged) cpabort(
"TRS4 did not converge")
667 IF (nspin == 1)
CALL dbcsr_scale(ls_scf_env%matrix_p(1), 2.0_dp)
669 CALL timestop(handle)
670 END SUBROUTINE pao_dm_trs4
679 TYPE(pao_env_type),
POINTER :: pao
680 TYPE(qs_environment_type),
POINTER :: qs_env
681 TYPE(ls_scf_env_type),
TARGET :: ls_scf_env
683 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_check_grad'
685 INTEGER :: handle, i, iatom, j, natoms
686 INTEGER,
DIMENSION(:),
POINTER :: blk_sizes_col, blk_sizes_row
688 REAL(
dp) :: delta, delta_max, eps, gij_num
689 REAL(
dp),
DIMENSION(:, :),
POINTER :: block_g, block_x
690 TYPE(ls_mstruct_type),
POINTER :: ls_mstruct
691 TYPE(mp_para_env_type),
POINTER :: para_env
693 IF (pao%check_grad_tol < 0.0_dp)
RETURN
695 CALL timeset(routinen, handle)
697 ls_mstruct => ls_scf_env%ls_mstruct
699 CALL get_qs_env(qs_env, para_env=para_env, natom=natoms)
701 eps = pao%num_grad_eps
704 CALL dbcsr_get_info(pao%matrix_X, col_blk_size=blk_sizes_col, row_blk_size=blk_sizes_row)
708 IF (pao%iw > 0)
WRITE (pao%iw, *)
'PAO| checking gradient of atom ', iatom
709 CALL dbcsr_get_block_p(matrix=pao%matrix_X, row=iatom, col=iatom, block=block_x, found=found)
711 IF (
ASSOCIATED(block_x))
THEN
712 CALL dbcsr_get_block_p(matrix=pao%matrix_G, row=iatom, col=iatom, block=block_g, found=found)
713 cpassert(
ASSOCIATED(block_g))
716 DO i = 1, blk_sizes_row(iatom)
717 DO j = 1, blk_sizes_col(iatom)
718 SELECT CASE (pao%num_grad_order)
720 gij_num = -eval_point(block_x, i, j, -eps, pao, ls_scf_env, qs_env)
721 gij_num = gij_num + eval_point(block_x, i, j, +eps, pao, ls_scf_env, qs_env)
722 gij_num = gij_num/(2.0_dp*eps)
725 gij_num = eval_point(block_x, i, j, -2_dp*eps, pao, ls_scf_env, qs_env)
726 gij_num = gij_num - 8_dp*eval_point(block_x, i, j, -1_dp*eps, pao, ls_scf_env, qs_env)
727 gij_num = gij_num + 8_dp*eval_point(block_x, i, j, +1_dp*eps, pao, ls_scf_env, qs_env)
728 gij_num = gij_num - eval_point(block_x, i, j, +2_dp*eps, pao, ls_scf_env, qs_env)
729 gij_num = gij_num/(12.0_dp*eps)
732 gij_num = -1_dp*eval_point(block_x, i, j, -3_dp*eps, pao, ls_scf_env, qs_env)
733 gij_num = gij_num + 9_dp*eval_point(block_x, i, j, -2_dp*eps, pao, ls_scf_env, qs_env)
734 gij_num = gij_num - 45_dp*eval_point(block_x, i, j, -1_dp*eps, pao, ls_scf_env, qs_env)
735 gij_num = gij_num + 45_dp*eval_point(block_x, i, j, +1_dp*eps, pao, ls_scf_env, qs_env)
736 gij_num = gij_num - 9_dp*eval_point(block_x, i, j, +2_dp*eps, pao, ls_scf_env, qs_env)
737 gij_num = gij_num + 1_dp*eval_point(block_x, i, j, +3_dp*eps, pao, ls_scf_env, qs_env)
738 gij_num = gij_num/(60.0_dp*eps)
741 cpabort(
"Unsupported numerical derivative order: "//cp_to_string(pao%num_grad_order))
744 IF (
ASSOCIATED(block_x))
THEN
745 delta = abs(gij_num - block_g(i, j))
746 delta_max = max(delta_max, delta)
753 CALL para_env%max(delta_max)
754 IF (pao%iw > 0)
WRITE (pao%iw, *)
'PAO| checked gradient, max delta:', delta_max
755 IF (delta_max > pao%check_grad_tol)
CALL cp_abort(__location__, &
756 "Analytic and numeric gradients differ too much:"//cp_to_string(delta_max))
758 CALL timestop(handle)
772 FUNCTION eval_point(block_X, i, j, eps, pao, ls_scf_env, qs_env)
RESULT(energy)
773 REAL(
dp),
DIMENSION(:, :),
POINTER :: block_x
774 INTEGER,
INTENT(IN) :: i, j
775 REAL(
dp),
INTENT(IN) :: eps
776 TYPE(pao_env_type),
POINTER :: pao
777 TYPE(ls_scf_env_type),
TARGET :: ls_scf_env
778 TYPE(qs_environment_type),
POINTER :: qs_env
783 IF (
ASSOCIATED(block_x))
THEN
784 old_xij = block_x(i, j)
785 block_x(i, j) = block_x(i, j) + eps
792 IF (
ASSOCIATED(block_x))
THEN
793 block_x(i, j) = old_xij
796 END FUNCTION eval_point
804 TYPE(qs_environment_type),
POINTER :: qs_env
805 TYPE(ls_scf_env_type),
TARGET :: ls_scf_env
807 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_store_P'
809 INTEGER :: handle, ispin, istore
810 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s
811 TYPE(dft_control_type),
POINTER :: dft_control
812 TYPE(ls_mstruct_type),
POINTER :: ls_mstruct
813 TYPE(pao_env_type),
POINTER :: pao
815 IF (ls_scf_env%scf_history%nstore == 0)
RETURN
816 CALL timeset(routinen, handle)
817 ls_mstruct => ls_scf_env%ls_mstruct
818 pao => ls_scf_env%pao_env
819 CALL get_qs_env(qs_env, dft_control=dft_control, matrix_s=matrix_s)
821 ls_scf_env%scf_history%istore = ls_scf_env%scf_history%istore + 1
822 istore = mod(ls_scf_env%scf_history%istore - 1, ls_scf_env%scf_history%nstore) + 1
823 IF (pao%iw > 0)
WRITE (pao%iw, *)
"PAO| Storing density matrix for ASPC guess in slot:", istore
826 IF (ls_scf_env%scf_history%istore <= ls_scf_env%scf_history%nstore)
THEN
827 DO ispin = 1, dft_control%nspins
828 CALL dbcsr_create(ls_scf_env%scf_history%matrix(ispin, istore), template=matrix_s(1)%matrix)
835 DO ispin = 1, dft_control%nspins
837 CALL matrix_ls_to_qs(ls_scf_env%scf_history%matrix(ispin, istore), ls_scf_env%matrix_p(ispin), &
838 ls_scf_env%ls_mstruct, covariant=.false., keep_sparsity=.false.)
841 CALL timestop(handle)
851 TYPE(pao_env_type),
POINTER :: pao
852 TYPE(qs_environment_type),
POINTER :: qs_env
853 TYPE(ls_scf_env_type),
TARGET :: ls_scf_env
855 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_guess_initial_P'
859 CALL timeset(routinen, handle)
861 IF (ls_scf_env%scf_history%istore > 0)
THEN
862 CALL pao_aspc_guess_p(pao, qs_env, ls_scf_env)
863 pao%need_initial_scf = .true.
865 IF (len_trim(pao%preopt_dm_file) > 0)
THEN
866 CALL pao_read_preopt_dm(pao, qs_env)
867 pao%need_initial_scf = .false.
868 pao%preopt_dm_file =
""
871 IF (pao%iw > 0)
WRITE (pao%iw,
'(A,F20.9)') &
872 " PAO| Energy from initial atomic guess:", ls_scf_env%energy_init
873 pao%need_initial_scf = .true.
877 CALL timestop(handle)
887 SUBROUTINE pao_aspc_guess_p(pao, qs_env, ls_scf_env)
888 TYPE(pao_env_type),
POINTER :: pao
889 TYPE(qs_environment_type),
POINTER :: qs_env
890 TYPE(ls_scf_env_type),
TARGET :: ls_scf_env
892 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_aspc_guess_P'
894 INTEGER :: handle, iaspc, ispin, istore, naspc
896 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s
897 TYPE(dbcsr_type) :: matrix_p
898 TYPE(dft_control_type),
POINTER :: dft_control
899 TYPE(ls_mstruct_type),
POINTER :: ls_mstruct
901 CALL timeset(routinen, handle)
902 ls_mstruct => ls_scf_env%ls_mstruct
903 cpassert(ls_scf_env%scf_history%istore > 0)
906 CALL get_qs_env(qs_env, dft_control=dft_control, matrix_s=matrix_s)
908 IF (pao%iw > 0)
WRITE (pao%iw, *)
"PAO| Calculating initial guess with ASPC"
910 CALL dbcsr_create(matrix_p, template=matrix_s(1)%matrix)
912 naspc = min(ls_scf_env%scf_history%istore, ls_scf_env%scf_history%nstore)
913 DO ispin = 1, dft_control%nspins
915 CALL dbcsr_set(matrix_p, 0.0_dp)
917 alpha = (-1.0_dp)**(iaspc + 1)*real(iaspc, kind=
dp)* &
919 istore = mod(ls_scf_env%scf_history%istore - iaspc, ls_scf_env%scf_history%nstore) + 1
920 CALL dbcsr_add(matrix_p, ls_scf_env%scf_history%matrix(ispin, istore), 1.0_dp, alpha)
924 CALL matrix_qs_to_ls(ls_scf_env%matrix_p(ispin), matrix_p, ls_scf_env%ls_mstruct, covariant=.false.)
927 CALL dbcsr_release(matrix_p)
930 DO ispin = 1, dft_control%nspins
931 IF (dft_control%nspins == 1)
CALL dbcsr_scale(ls_scf_env%matrix_p(ispin), 0.5_dp)
933 CALL dbcsr_filter(ls_scf_env%matrix_p(ispin), ls_scf_env%eps_filter**(2.0_dp/3.0_dp))
936 CALL purify_mcweeny(ls_scf_env%matrix_p(ispin:ispin), ls_scf_env%matrix_s, ls_scf_env%eps_filter, 10)
937 IF (dft_control%nspins == 1)
CALL dbcsr_scale(ls_scf_env%matrix_p(ispin), 2.0_dp)
943 CALL ls_scf_dm_to_ks(qs_env, ls_scf_env, ls_scf_env%energy_init, iscf=0)
945 CALL timestop(handle)
946 END SUBROUTINE pao_aspc_guess_p
954 TYPE(qs_environment_type),
POINTER :: qs_env
955 TYPE(ls_scf_env_type),
TARGET :: ls_scf_env
957 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_add_forces'
959 INTEGER :: handle, iatom, natoms
960 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: forces
961 TYPE(mp_para_env_type),
POINTER :: para_env
962 TYPE(pao_env_type),
POINTER :: pao
963 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
965 CALL timeset(routinen, handle)
966 pao => ls_scf_env%pao_env
968 IF (pao%iw > 0)
WRITE (pao%iw, *)
"PAO| Adding forces."
970 IF (pao%max_pao /= 0)
THEN
971 IF (pao%penalty_strength /= 0.0_dp) &
972 cpabort(
"PAO forces require PENALTY_STRENGTH or MAX_PAO set to zero")
973 IF (pao%linpot_regu_strength /= 0.0_dp) &
974 cpabort(
"PAO forces require LINPOT_REGULARIZATION_STRENGTH or MAX_PAO set to zero")
975 IF (pao%regularization /= 0.0_dp) &
976 cpabort(
"PAO forces require REGULARIZATION or MAX_PAO set to zero")
981 particle_set=particle_set, &
984 ALLOCATE (forces(natoms, 3))
985 CALL pao_calc_ab(pao, qs_env, ls_scf_env, gradient=.true., forces=forces)
987 IF (
SIZE(pao%ml_training_set) > 0) &
990 CALL para_env%sum(forces)
992 particle_set(iatom)%f = particle_set(iatom)%f + forces(iatom, :)
997 CALL timestop(handle)
All kind of helpful little routines.
real(kind=dp) function, public exp_radius(l, alpha, threshold, prefactor, epsabs, epsrel, rlow)
The radius of a primitive Gaussian function for a given threshold is calculated. g(r) = prefactor*r**...
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 kuhne2007
integer, save, public kolafa2004
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
lower level routines for linear scaling SCF
subroutine, public ls_scf_init_matrix_s(matrix_s, ls_scf_env)
initialize S matrix related properties (sqrt, inverse...) Might be factored-out since this seems comm...
subroutine, public density_matrix_trs4(matrix_p, matrix_ks, matrix_s_sqrt_inv, nelectron, threshold, e_homo, e_lumo, e_mu, dynamic_threshold, matrix_ks_deviation, max_iter_lanczos, eps_lanczos, converged)
compute the density matrix using a trace-resetting algorithm
Routines for a linear scaling quickstep SCF run based on the density matrix, with a focus on the inte...
subroutine, public matrix_ls_to_qs(matrix_qs, matrix_ls, ls_mstruct, covariant, keep_sparsity)
second link to QS, copy a LS matrix to QS matrix used to isolate QS style matrices from LS style will...
subroutine, public ls_scf_dm_to_ks(qs_env, ls_scf_env, energy_new, iscf)
use the density matrix in ls_scf_env to compute the new energy and KS matrix
subroutine, public matrix_qs_to_ls(matrix_ls, matrix_qs, ls_mstruct, covariant)
first link to QS, copy a QS matrix to LS matrix used to isolate QS style matrices from LS style will ...
subroutine, public ls_scf_qs_atomic_guess(qs_env, energy)
get an atomic initial guess
Types needed for a linear scaling quickstep SCF run based on the density matrix.
Routines useful for iterative matrix calculations.
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_path_length
Machine interface based on Fortran 2003 and POSIX.
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Collection of simple mathematical functions and subroutines.
elemental real(kind=dp) function, public binomial(n, k)
The binomial coefficient n over k for 0 <= k <= n is calculated, otherwise zero is returned.
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Interface to the message passing library MPI.
Methods used by pao_main.F.
subroutine, public pao_build_orthogonalizer(pao, qs_env)
Constructs matrix_N and its inverse.
subroutine, public pao_build_core_hamiltonian(pao, qs_env)
Creates the matrix_H0 which contains the core hamiltonian.
subroutine, public pao_guess_initial_p(pao, qs_env, ls_scf_env)
Provide an initial guess for the density matrix.
subroutine, public pao_test_convergence(pao, ls_scf_env, new_energy, is_converged)
Test whether the PAO optimization has reached convergence.
subroutine, public pao_add_forces(qs_env, ls_scf_env)
Calculate the forces contributed by PAO.
subroutine, public pao_check_trace_ps(ls_scf_env)
Ensure that the number of electrons is correct.
subroutine, public pao_init_kinds(pao, qs_env)
Initialize qs kinds.
subroutine, public pao_build_matrix_x(pao, qs_env)
Creates the matrix_X.
subroutine, public pao_check_grad(pao, qs_env, ls_scf_env)
Debugging routine for checking the analytic gradient.
subroutine, public pao_print_atom_info(pao)
Prints a one line summary for each atom.
subroutine, public pao_store_p(qs_env, ls_scf_env)
Stores density matrix as initial guess for next SCF optimization.
subroutine, public pao_build_selector(pao, qs_env)
Build rectangular matrix to converert between primary and PAO basis.
subroutine, public pao_calc_energy(pao, qs_env, ls_scf_env, energy)
Calculate the pao energy.
subroutine, public pao_build_diag_distribution(pao, qs_env)
Creates new DBCSR distribution which spreads diagonal blocks evenly across ranks.
Main module for PAO Machine Learning.
subroutine, public pao_ml_forces(pao, qs_env, matrix_G, forces)
Calculate forces contributed by machine learning.
Front-End for any PAO parametrization.
subroutine, public pao_calc_ab(pao, qs_env, ls_scf_env, gradient, penalty, forces)
Takes current matrix_X and calculates the matrices A and B.
subroutine, public pao_param_count(pao, qs_env, ikind, nparams)
Returns the number of parameters for given atomic kind.
Factory routines for potentials used e.g. by pao_param_exp and pao_ml.
Types used by the PAO machinery.
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.
Routines to somehow generate an initial guess.
subroutine, public calculate_atomic_fock_matrix(matrix_f, atomic_kind_set, qs_kind_set, ounit)
returns a block diagonal fock matrix.
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 set_qs_kind(qs_kind, paw_atom, ghost, floating, hard_radius, hard0_radius, covalent_radius, vdw_radius, lmax_rho0, zeff, no_optimize, dispersion, u_minus_j, reltmat, dftb_parameter, xtb_parameter, elec_conf, pao_basis_size)
Set the components of an atomic kind data set.
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
subroutine, public qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, print_active)
updates the Kohn Sham matrix of the given qs_env (facility method)
subroutine, public qs_ks_did_change(ks_env, s_mstruct_changed, rho_changed, potential_changed, full_reset)
tells that some of the things relevant to the ks calculation did change. has to be called when change...
methods of the rho structure (defined in qs_rho_types)
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
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...