(git:58e3e09)
pao_param_equi.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 Equivariant parametrization
10 !> \author Ole Schuett
11 ! **************************************************************************************************
13  USE basis_set_types, ONLY: gto_basis_set_type
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
19  USE dm_ls_scf_types, ONLY: ls_mstruct_type,&
20  ls_scf_env_type
21  USE kinds, ONLY: dp
22  USE mathlib, ONLY: diamat_all
23  USE message_passing, ONLY: mp_comm_type
26  USE pao_types, ONLY: pao_env_type
27  USE qs_environment_types, ONLY: get_qs_env,&
28  qs_environment_type
29  USE qs_kind_types, ONLY: get_qs_kind,&
30  qs_kind_type
31 #include "./base/base_uses.f90"
32 
33  IMPLICIT NONE
34 
35  PRIVATE
36 
37  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_param_equi'
38 
41 
42 CONTAINS
43 
44 ! **************************************************************************************************
45 !> \brief Initialize equivariant parametrization
46 !> \param pao ...
47 ! **************************************************************************************************
48  SUBROUTINE pao_param_init_equi(pao)
49  TYPE(pao_env_type), POINTER :: pao
50 
51  IF (pao%precondition) &
52  cpabort("PAO preconditioning not supported for selected parametrization.")
53 
54  END SUBROUTINE pao_param_init_equi
55 
56 ! **************************************************************************************************
57 !> \brief Finalize equivariant parametrization
58 ! **************************************************************************************************
60 
61  ! Nothing to do.
62 
63  END SUBROUTINE pao_param_finalize_equi
64 
65 ! **************************************************************************************************
66 !> \brief Returns the number of parameters for given atomic kind
67 !> \param qs_env ...
68 !> \param ikind ...
69 !> \param nparams ...
70 ! **************************************************************************************************
71  SUBROUTINE pao_param_count_equi(qs_env, ikind, nparams)
72  TYPE(qs_environment_type), POINTER :: qs_env
73  INTEGER, INTENT(IN) :: ikind
74  INTEGER, INTENT(OUT) :: nparams
75 
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
79 
80  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
81  CALL get_qs_kind(qs_kind_set(ikind), &
82  basis_set=basis_set, &
83  pao_basis_size=pao_basis_size)
84  pri_basis_size = basis_set%nsgf
85 
86  nparams = pao_basis_size*pri_basis_size
87 
88  END SUBROUTINE pao_param_count_equi
89 
90 ! **************************************************************************************************
91 !> \brief Fills matrix_X with an initial guess
92 !> \param pao ...
93 !> \param qs_env ...
94 ! **************************************************************************************************
95  SUBROUTINE pao_param_initguess_equi(pao, qs_env)
96  TYPE(pao_env_type), POINTER :: pao
97  TYPE(qs_environment_type), POINTER :: qs_env
98 
99  CHARACTER(len=*), PARAMETER :: routinen = 'pao_param_initguess_equi'
100 
101  INTEGER :: acol, arow, handle, i, iatom, m, n
102  INTEGER, DIMENSION(:), POINTER :: blk_sizes_pao, blk_sizes_pri
103  LOGICAL :: found
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
108 
109  CALL timeset(routinen, handle)
110 
111  CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=blk_sizes_pri, col_blk_size=blk_sizes_pao)
112 
113 !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,qs_env,blk_sizes_pri,blk_sizes_pao) &
114 !$OMP PRIVATE(iter,arow,acol,iatom,n,m,i,found) &
115 !$OMP PRIVATE(block_X,block_H0,block_N,block_N_inv,A,H,H_evecs,H_evals,V0)
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)
120 
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))
125 
126  n = blk_sizes_pri(iatom) ! size of primary basis
127  m = blk_sizes_pao(iatom) ! size of pao basis
128 
129  ALLOCATE (v0(n, n))
130  CALL pao_guess_initial_potential(qs_env, iatom, v0)
131 
132  ! construct H
133  ALLOCATE (h(n, n))
134  h = matmul(matmul(block_n, block_h0 + v0), block_n) ! transform into orthonormal basis
135 
136  ! diagonalize H
137  ALLOCATE (h_evecs(n, n), h_evals(n))
138  h_evecs = h
139  CALL diamat_all(h_evecs, h_evals)
140 
141  ! use first m eigenvectors as initial guess
142  ALLOCATE (a(n, m))
143  a = matmul(block_n_inv, h_evecs(:, 1:m))
144 
145  ! normalize vectors
146  DO i = 1, m
147  a(:, i) = a(:, i)/norm2(a(:, i))
148  END DO
149 
150  block_x = reshape(a, (/n*m, 1/))
151  DEALLOCATE (h, v0, a, h_evecs, h_evals)
152 
153  END DO
154  CALL dbcsr_iterator_stop(iter)
155 !$OMP END PARALLEL
156 
157  CALL timestop(handle)
158 
159  END SUBROUTINE pao_param_initguess_equi
160 
161 ! **************************************************************************************************
162 !> \brief Takes current matrix_X and calculates the matrices A and B.
163 !> \param pao ...
164 !> \param qs_env ...
165 !> \param ls_scf_env ...
166 !> \param gradient ...
167 !> \param penalty ...
168 ! **************************************************************************************************
169  SUBROUTINE pao_calc_ab_equi(pao, qs_env, ls_scf_env, gradient, penalty)
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
175 
176  CHARACTER(len=*), PARAMETER :: routinen = 'pao_calc_AB_equi'
177 
178  INTEGER :: acol, arow, group_handle, handle, i, &
179  iatom, j, k, m, n
180  LOGICAL :: found
181  REAL(dp) :: denom, w
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, &
186  m4, m5, nn
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, &
191  matrix_x_nondiag
192  TYPE(ls_mstruct_type), POINTER :: ls_mstruct
193  TYPE(mp_comm_type) :: group
194 
195  CALL timeset(routinen, handle)
196  ls_mstruct => ls_scf_env%ls_mstruct
197 
198  IF (gradient) THEN
199  CALL pao_calc_grad_lnv_wrt_ab(qs_env, ls_scf_env, matrix_ma, matrix_mb)
200  END IF
201 
202  ! Redistribute matrix_X from diag_distribution to distribution of matrix_s.
203  CALL get_qs_env(qs_env, matrix_s=matrix_s)
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", &
207  dist=main_dist, &
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)
211 
212  ! Compuation of matrix_G uses distr. of matrix_s, afterwards we redistribute to diag_distribution.
213  IF (gradient) THEN
214  CALL dbcsr_create(matrix_g_nondiag, &
215  name="PAO matrix_G_nondiag", &
216  dist=main_dist, &
217  template=pao%matrix_G)
218  CALL dbcsr_reserve_diag_blocks(matrix_g_nondiag)
219  END IF
220 
221 !$OMP PARALLEL DEFAULT(NONE) &
222 !$OMP SHARED(pao,ls_mstruct,matrix_X_nondiag,matrix_G_nondiag,matrix_Ma,matrix_Mb,gradient,penalty) &
223 !$OMP PRIVATE(iter,arow,acol,iatom,found,n,m,w,i,j,k,denom) &
224 !$OMP PRIVATE(NN,ANNA,ANNA_evals,ANNA_evecs,ANNA_inv,D,G,M1,M2,M3,M4,M5) &
225 !$OMP PRIVATE(block_X,block_A,block_B,block_N,block_Ma, block_Mb, block_G)
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))
236 
237  n = SIZE(block_a, 1) ! size of primary basis
238  m = SIZE(block_a, 2) ! size of pao basis
239  block_a = reshape(block_x, (/n, m/))
240 
241  ! restrain pao basis vectors to unit norm
242  IF (PRESENT(penalty)) THEN
243  DO i = 1, m
244  w = 1.0_dp - sum(block_a(:, i)**2)
245  penalty = penalty + pao%penalty_strength*w**2
246  END DO
247  END IF
248 
249  ALLOCATE (nn(n, n), anna(m, m))
250  nn = matmul(block_n, block_n) ! it's actually S^{-1}
251  anna = matmul(matmul(transpose(block_a), nn), block_a)
252 
253  ! diagonalize ANNA
254  ALLOCATE (anna_evecs(m, m), anna_evals(m))
255  anna_evecs(:, :) = anna
256  CALL diamat_all(anna_evecs, anna_evals)
257  IF (minval(abs(anna_evals)) < 1e-10_dp) cpabort("PAO basis singualar.")
258 
259  ! build ANNA_inv
260  ALLOCATE (anna_inv(m, m))
261  anna_inv(:, :) = 0.0_dp
262  DO k = 1, m
263  w = 1.0_dp/anna_evals(k)
264  DO i = 1, m
265  DO j = 1, m
266  anna_inv(i, j) = anna_inv(i, j) + w*anna_evecs(i, k)*anna_evecs(j, k)
267  END DO
268  END DO
269  END DO
270 
271  !B = 1/S * A * 1/(A^T 1/S A)
272  block_b = matmul(matmul(nn, block_a), anna_inv)
273 
274  ! TURNING POINT (if calc grad) ------------------------------------------
275  IF (gradient) THEN
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)
280  ! don't check ASSOCIATED(block_M), it might have been filtered out.
281 
282  ALLOCATE (g(n, m))
283  g(:, :) = 0.0_dp
284 
285  IF (PRESENT(penalty)) THEN
286  DO i = 1, m
287  w = 1.0_dp - sum(block_a(:, i)**2)
288  g(:, i) = -4.0_dp*pao%penalty_strength*w*block_a(:, i)
289  END DO
290  END IF
291 
292  IF (ASSOCIATED(block_ma)) THEN
293  g = g + block_ma
294  END IF
295 
296  IF (ASSOCIATED(block_mb)) THEN
297  g = g + matmul(matmul(nn, block_mb), anna_inv)
298 
299  ! calculate derivatives dAA_inv/ dAA
300  ALLOCATE (d(m, m), m1(m, m), m2(m, m), m3(m, m), m4(m, m), m5(m, m))
301 
302  DO i = 1, m
303  DO j = 1, m
304  denom = anna_evals(i) - anna_evals(j)
305  IF (i == j) THEN
306  d(i, i) = -1.0_dp/anna_evals(i)**2 ! diagonal elements
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
309  ELSE
310  d(i, j) = -1.0_dp ! limit according to L'Hospital's rule
311  END IF
312  END DO
313  END DO
314 
315  m1 = matmul(matmul(transpose(block_a), nn), block_mb)
316  m2 = matmul(matmul(transpose(anna_evecs), m1), anna_evecs)
317  m3 = m2*d ! Hadamard product
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)
321 
322  DEALLOCATE (d, m1, m2, m3, m4, m5)
323  END IF
324 
325  block_g = reshape(g, (/n*m, 1/))
326  DEALLOCATE (g)
327  END IF
328 
329  DEALLOCATE (nn, anna, anna_evecs, anna_evals, anna_inv)
330  END DO
331  CALL dbcsr_iterator_stop(iter)
332 !$OMP END PARALLEL
333 
334  ! sum penalty energies across ranks
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)
339  END IF
340 
341  CALL dbcsr_release(matrix_x_nondiag)
342 
343  IF (gradient) THEN
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)
348  END IF
349 
350  CALL timestop(handle)
351 
352  END SUBROUTINE pao_calc_ab_equi
353 
354 END MODULE pao_param_equi
Types needed for a linear scaling quickstep SCF run based on the density matrix.
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
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...
Definition: mathlib.F:376
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.
Definition: pao_types.F:12
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.
Definition: qs_kind_types.F:23
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.