14 USE dbcsr_api,
ONLY: &
15 dbcsr_complete_redistribute, dbcsr_create, dbcsr_distribution_type, dbcsr_get_block_p, &
16 dbcsr_get_info, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
17 dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
18 dbcsr_release, dbcsr_reserve_diag_blocks, dbcsr_type
31 #include "./base/base_uses.f90"
37 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'pao_param_equi'
49 TYPE(pao_env_type),
POINTER :: pao
51 IF (pao%precondition) &
52 cpabort(
"PAO preconditioning not supported for selected parametrization.")
72 TYPE(qs_environment_type),
POINTER :: qs_env
73 INTEGER,
INTENT(IN) :: ikind
74 INTEGER,
INTENT(OUT) :: nparams
76 INTEGER :: pao_basis_size, pri_basis_size
77 TYPE(gto_basis_set_type),
POINTER :: basis_set
78 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
80 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
82 basis_set=basis_set, &
83 pao_basis_size=pao_basis_size)
84 pri_basis_size = basis_set%nsgf
86 nparams = pao_basis_size*pri_basis_size
96 TYPE(pao_env_type),
POINTER :: pao
97 TYPE(qs_environment_type),
POINTER :: qs_env
99 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_param_initguess_equi'
101 INTEGER :: acol, arow, handle, i, iatom, m, n
102 INTEGER,
DIMENSION(:),
POINTER :: blk_sizes_pao, blk_sizes_pri
104 REAL(
dp),
DIMENSION(:),
POINTER :: h_evals
105 REAL(
dp),
DIMENSION(:, :),
POINTER :: a, block_h0, block_n, block_n_inv, &
106 block_x, h, h_evecs, v0
107 TYPE(dbcsr_iterator_type) :: iter
109 CALL timeset(routinen, handle)
111 CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=blk_sizes_pri, col_blk_size=blk_sizes_pao)
116 CALL dbcsr_iterator_start(iter, pao%matrix_X)
117 DO WHILE (dbcsr_iterator_blocks_left(iter))
118 CALL dbcsr_iterator_next_block(iter, arow, acol, block_x)
119 iatom = arow; cpassert(arow == acol)
121 CALL dbcsr_get_block_p(matrix=pao%matrix_H0, row=iatom, col=iatom, block=block_h0, found=found)
122 CALL dbcsr_get_block_p(matrix=pao%matrix_N_diag, row=iatom, col=iatom, block=block_n, found=found)
123 CALL dbcsr_get_block_p(matrix=pao%matrix_N_inv_diag, row=iatom, col=iatom, block=block_n_inv, found=found)
124 cpassert(
ASSOCIATED(block_h0) .AND.
ASSOCIATED(block_n) .AND.
ASSOCIATED(block_n_inv))
126 n = blk_sizes_pri(iatom)
127 m = blk_sizes_pao(iatom)
134 h = matmul(matmul(block_n, block_h0 + v0), block_n)
137 ALLOCATE (h_evecs(n, n), h_evals(n))
143 a = matmul(block_n_inv, h_evecs(:, 1:m))
147 a(:, i) = a(:, i)/norm2(a(:, i))
150 block_x = reshape(a, (/n*m, 1/))
151 DEALLOCATE (h, v0, a, h_evecs, h_evals)
154 CALL dbcsr_iterator_stop(iter)
157 CALL timestop(handle)
170 TYPE(pao_env_type),
POINTER :: pao
171 TYPE(qs_environment_type),
POINTER :: qs_env
172 TYPE(ls_scf_env_type),
TARGET :: ls_scf_env
173 LOGICAL,
INTENT(IN) :: gradient
174 REAL(
dp),
INTENT(INOUT),
OPTIONAL :: penalty
176 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_calc_AB_equi'
178 INTEGER :: acol, arow, group_handle, handle, i, &
182 REAL(
dp),
DIMENSION(:),
POINTER :: anna_evals
183 REAL(
dp),
DIMENSION(:, :),
POINTER :: anna, anna_evecs, anna_inv, block_a, &
184 block_b, block_g, block_ma, block_mb, &
185 block_n, block_x, d, g, m1, m2, m3, &
187 TYPE(dbcsr_distribution_type) :: main_dist
188 TYPE(dbcsr_iterator_type) :: iter
189 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s
190 TYPE(dbcsr_type) :: matrix_g_nondiag, matrix_ma, matrix_mb, &
192 TYPE(ls_mstruct_type),
POINTER :: ls_mstruct
193 TYPE(mp_comm_type) :: group
195 CALL timeset(routinen, handle)
196 ls_mstruct => ls_scf_env%ls_mstruct
204 CALL dbcsr_get_info(matrix=matrix_s(1)%matrix, distribution=main_dist)
205 CALL dbcsr_create(matrix_x_nondiag, &
206 name=
"PAO matrix_X_nondiag", &
208 template=pao%matrix_X)
209 CALL dbcsr_reserve_diag_blocks(matrix_x_nondiag)
210 CALL dbcsr_complete_redistribute(pao%matrix_X, matrix_x_nondiag)
214 CALL dbcsr_create(matrix_g_nondiag, &
215 name=
"PAO matrix_G_nondiag", &
217 template=pao%matrix_G)
218 CALL dbcsr_reserve_diag_blocks(matrix_g_nondiag)
226 CALL dbcsr_iterator_start(iter, matrix_x_nondiag)
227 DO WHILE (dbcsr_iterator_blocks_left(iter))
228 CALL dbcsr_iterator_next_block(iter, arow, acol, block_x)
229 iatom = arow; cpassert(arow == acol)
230 CALL dbcsr_get_block_p(matrix=ls_mstruct%matrix_A, row=iatom, col=iatom, block=block_a, found=found)
231 cpassert(
ASSOCIATED(block_a))
232 CALL dbcsr_get_block_p(matrix=ls_mstruct%matrix_B, row=iatom, col=iatom, block=block_b, found=found)
233 cpassert(
ASSOCIATED(block_b))
234 CALL dbcsr_get_block_p(matrix=pao%matrix_N, row=iatom, col=iatom, block=block_n, found=found)
235 cpassert(
ASSOCIATED(block_n))
239 block_a = reshape(block_x, (/n, m/))
242 IF (
PRESENT(penalty))
THEN
244 w = 1.0_dp - sum(block_a(:, i)**2)
245 penalty = penalty + pao%penalty_strength*w**2
249 ALLOCATE (nn(n, n), anna(m, m))
250 nn = matmul(block_n, block_n)
251 anna = matmul(matmul(transpose(block_a), nn), block_a)
254 ALLOCATE (anna_evecs(m, m), anna_evals(m))
255 anna_evecs(:, :) = anna
257 IF (minval(abs(anna_evals)) < 1e-10_dp) cpabort(
"PAO basis singualar.")
260 ALLOCATE (anna_inv(m, m))
261 anna_inv(:, :) = 0.0_dp
263 w = 1.0_dp/anna_evals(k)
266 anna_inv(i, j) = anna_inv(i, j) + w*anna_evecs(i, k)*anna_evecs(j, k)
272 block_b = matmul(matmul(nn, block_a), anna_inv)
276 CALL dbcsr_get_block_p(matrix=matrix_g_nondiag, row=iatom, col=iatom, block=block_g, found=found)
277 cpassert(
ASSOCIATED(block_g))
278 CALL dbcsr_get_block_p(matrix=matrix_ma, row=iatom, col=iatom, block=block_ma, found=found)
279 CALL dbcsr_get_block_p(matrix=matrix_mb, row=iatom, col=iatom, block=block_mb, found=found)
285 IF (
PRESENT(penalty))
THEN
287 w = 1.0_dp - sum(block_a(:, i)**2)
288 g(:, i) = -4.0_dp*pao%penalty_strength*w*block_a(:, i)
292 IF (
ASSOCIATED(block_ma))
THEN
296 IF (
ASSOCIATED(block_mb))
THEN
297 g = g + matmul(matmul(nn, block_mb), anna_inv)
300 ALLOCATE (d(m, m), m1(m, m), m2(m, m), m3(m, m), m4(m, m), m5(m, m))
304 denom = anna_evals(i) - anna_evals(j)
306 d(i, i) = -1.0_dp/anna_evals(i)**2
307 ELSE IF (abs(denom) > 1e-10_dp)
THEN
308 d(i, j) = (1.0_dp/anna_evals(i) - 1.0_dp/anna_evals(j))/denom
315 m1 = matmul(matmul(transpose(block_a), nn), block_mb)
316 m2 = matmul(matmul(transpose(anna_evecs), m1), anna_evecs)
318 m4 = matmul(matmul(anna_evecs, m3), transpose(anna_evecs))
319 m5 = 0.5_dp*(m4 + transpose(m4))
320 g = g + 2.0_dp*matmul(matmul(nn, block_a), m5)
322 DEALLOCATE (d, m1, m2, m3, m4, m5)
325 block_g = reshape(g, (/n*m, 1/))
329 DEALLOCATE (nn, anna, anna_evecs, anna_evals, anna_inv)
331 CALL dbcsr_iterator_stop(iter)
335 IF (
PRESENT(penalty))
THEN
336 CALL dbcsr_get_info(pao%matrix_X, group=group_handle)
337 CALL group%set_handle(group_handle)
338 CALL group%sum(penalty)
341 CALL dbcsr_release(matrix_x_nondiag)
344 CALL dbcsr_complete_redistribute(matrix_g_nondiag, pao%matrix_G)
345 CALL dbcsr_release(matrix_g_nondiag)
346 CALL dbcsr_release(matrix_ma)
347 CALL dbcsr_release(matrix_mb)
350 CALL timestop(handle)
Types needed for a linear scaling quickstep SCF run based on the density matrix.
Defines the basic variable types.
integer, parameter, public dp
Collection of simple mathematical functions and subroutines.
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.
Equivariant parametrization.
subroutine, public pao_param_count_equi(qs_env, ikind, nparams)
Returns the number of parameters for given atomic kind.
subroutine, public pao_param_finalize_equi()
Finalize equivariant parametrization.
subroutine, public pao_param_initguess_equi(pao, qs_env)
Fills matrix_X with an initial guess.
subroutine, public pao_param_init_equi(pao)
Initialize equivariant parametrization.
subroutine, public pao_calc_ab_equi(pao, qs_env, ls_scf_env, gradient, penalty)
Takes current matrix_X and calculates the matrices A and B.
Common routines for PAO parametrizations.
subroutine, public pao_calc_grad_lnv_wrt_ab(qs_env, ls_scf_env, matrix_Ma, matrix_Mb)
Helper routine, calculates partial derivative dE/dA and dE/dB. As energy functional serves the defini...
Factory routines for potentials used e.g. by pao_param_exp and pao_ml.
subroutine, public pao_guess_initial_potential(qs_env, iatom, block_V)
Makes an educated guess for the initial potential based on positions of neighboring atoms.
Types used by the PAO machinery.
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.