(git:0de0cc2)
mao_basis.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 Calculate MAO's and analyze wavefunctions
10 !> \par History
11 !> 03.2016 created [JGH]
12 !> 12.2016 split into four modules [JGH]
13 !> \author JGH
14 ! **************************************************************************************************
15 MODULE mao_basis
17  USE basis_set_types, ONLY: gto_basis_set_p_type,&
18  gto_basis_set_type
19  USE cp_control_types, ONLY: dft_control_type
22  USE dbcsr_api, ONLY: dbcsr_create,&
23  dbcsr_distribution_type,&
24  dbcsr_p_type,&
25  dbcsr_reserve_diag_blocks,&
26  dbcsr_type_no_symmetry
27  USE kinds, ONLY: dp
28  USE mao_methods, ONLY: mao_build_q
29  USE mao_optimizer, ONLY: mao_optimize
31  USE particle_types, ONLY: particle_type
32  USE qs_environment_types, ONLY: get_qs_env,&
33  qs_environment_type
34  USE qs_kind_types, ONLY: get_qs_kind,&
35  qs_kind_type
36  USE qs_ks_types, ONLY: get_ks_env,&
37  qs_ks_env_type
38  USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
42  USE qs_rho_types, ONLY: qs_rho_get,&
43  qs_rho_type
44 #include "./base/base_uses.f90"
45 
46  IMPLICIT NONE
47  PRIVATE
48 
49  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mao_basis'
50 
51  PUBLIC :: mao_generate_basis
52 
53 ! **************************************************************************************************
54 
55 CONTAINS
56 
57 ! **************************************************************************************************
58 !> \brief ...
59 !> \param qs_env ...
60 !> \param mao_coef ...
61 !> \param ref_basis_set ...
62 !> \param pmat_external ...
63 !> \param smat_external ...
64 !> \param molecular ...
65 !> \param max_iter ...
66 !> \param eps_grad ...
67 !> \param nmao_external ...
68 !> \param eps1_mao ...
69 !> \param iolevel ...
70 !> \param unit_nr ...
71 ! **************************************************************************************************
72  SUBROUTINE mao_generate_basis(qs_env, mao_coef, ref_basis_set, pmat_external, smat_external, &
73  molecular, max_iter, eps_grad, nmao_external, &
74  eps1_mao, iolevel, unit_nr)
75  TYPE(qs_environment_type), POINTER :: qs_env
76  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef
77  CHARACTER(len=*), OPTIONAL :: ref_basis_set
78  TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
79  POINTER :: pmat_external, smat_external
80  LOGICAL, INTENT(IN), OPTIONAL :: molecular
81  INTEGER, INTENT(IN), OPTIONAL :: max_iter
82  REAL(kind=dp), INTENT(IN), OPTIONAL :: eps_grad
83  INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: nmao_external
84  REAL(kind=dp), INTENT(IN), OPTIONAL :: eps1_mao
85  INTEGER, INTENT(IN), OPTIONAL :: iolevel, unit_nr
86 
87  CHARACTER(len=*), PARAMETER :: routinen = 'mao_generate_basis'
88 
89  CHARACTER(len=10) :: mao_basis_set
90  INTEGER :: handle, iab, iatom, ikind, iolev, ispin, &
91  iw, mao_max_iter, natom, nbas, &
92  nimages, nkind, nmao, nspin
93  INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
94  LOGICAL :: do_nmao_external, molecule
95  REAL(kind=dp) :: electra(2), eps1, eps_filter, eps_fun, &
96  mao_eps_grad
97  TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
98  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_q, matrix_smm, matrix_smo
99  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
100  TYPE(dft_control_type), POINTER :: dft_control
101  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: mao_basis_set_list, orb_basis_set_list
102  TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
103  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
104  POINTER :: smm_list, smo_list
105  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
106  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
107  TYPE(qs_kind_type), POINTER :: qs_kind
108  TYPE(qs_ks_env_type), POINTER :: ks_env
109  TYPE(qs_rho_type), POINTER :: rho
110 
111  CALL timeset(routinen, handle)
112 
113  ! k-points?
114  CALL get_qs_env(qs_env, dft_control=dft_control)
115  nimages = dft_control%nimages
116  cpassert(nimages == 1)
117 
118  ! output
119  IF (PRESENT(unit_nr)) THEN
120  iw = unit_nr
121  ELSE
122  iw = -1
123  END IF
124  IF (PRESENT(iolevel)) THEN
125  iolev = iolevel
126  ELSE
127  iolev = 1
128  END IF
129  IF (iolevel == 0) THEN
130  iw = -1
131  END IF
132 
133  ! molecules
134  IF (PRESENT(molecular)) THEN
135  molecule = molecular
136  ELSE
137  molecule = .false.
138  END IF
139 
140  ! iterations
141  IF (PRESENT(max_iter)) THEN
142  mao_max_iter = max_iter
143  ELSE
144  mao_max_iter = 0
145  END IF
146 
147  ! threshold
148  IF (PRESENT(eps_grad)) THEN
149  mao_eps_grad = eps_grad
150  ELSE
151  mao_eps_grad = 0.00001_dp
152  END IF
153  eps_fun = 10._dp*mao_eps_grad
154 
155  do_nmao_external = .false.
156  ! external number of MAOs per atom
157  IF (PRESENT(nmao_external)) THEN
158  do_nmao_external = .true.
159  END IF
160 
161  ! mao_threshold
162  IF (PRESENT(eps1_mao)) THEN
163  eps1 = eps1_mao
164  ELSE
165  eps1 = 1000._dp
166  END IF
167 
168  IF (iw > 0) THEN
169  WRITE (iw, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
170  WRITE (unit=iw, fmt="(T37,A)") "MAO BASIS"
171  WRITE (iw, '(T2,A)') '!-----------------------------------------------------------------------------!'
172  END IF
173 
174  ! Reference basis set
175  IF (PRESENT(ref_basis_set)) THEN
176  mao_basis_set = ref_basis_set
177  ELSE
178  mao_basis_set = "ORB"
179  END IF
180 
181  CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
182  nkind = SIZE(qs_kind_set)
183  ALLOCATE (mao_basis_set_list(nkind), orb_basis_set_list(nkind))
184  DO ikind = 1, nkind
185  qs_kind => qs_kind_set(ikind)
186  NULLIFY (mao_basis_set_list(ikind)%gto_basis_set)
187  NULLIFY (orb_basis_set_list(ikind)%gto_basis_set)
188  NULLIFY (basis_set_a, basis_set_b)
189  CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type="ORB")
190  IF (ASSOCIATED(basis_set_a)) orb_basis_set_list(ikind)%gto_basis_set => basis_set_a
191  CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_b, basis_type=mao_basis_set)
192  IF (ASSOCIATED(basis_set_b)) mao_basis_set_list(ikind)%gto_basis_set => basis_set_b
193  END DO
194  IF (iw > 0) THEN
195  DO ikind = 1, nkind
196  IF (.NOT. ASSOCIATED(mao_basis_set_list(ikind)%gto_basis_set)) THEN
197  WRITE (unit=iw, fmt="(T2,A,I4)") &
198  "WARNING: No MAO basis set associated with Kind ", ikind
199  ELSE
200  IF (do_nmao_external) THEN
201  nmao = nmao_external(ikind)
202  ELSE
203  CALL get_qs_kind(qs_kind_set(ikind), mao=nmao)
204  END IF
205  nbas = mao_basis_set_list(ikind)%gto_basis_set%nsgf
206  WRITE (unit=iw, fmt="(T2,A,I4,T30,A,I2,T56,A,I10)") &
207  "MAO basis set Kind ", ikind, " MinNum of MAO:", nmao, " Number of BSF:", nbas
208  END IF
209  END DO
210  END IF
211 
212  ! neighbor lists
213  NULLIFY (smm_list, smo_list)
214  CALL setup_neighbor_list(smm_list, mao_basis_set_list, molecular=molecule, qs_env=qs_env)
215  CALL setup_neighbor_list(smo_list, mao_basis_set_list, orb_basis_set_list, &
216  molecular=molecule, qs_env=qs_env)
217 
218  ! overlap matrices
219  NULLIFY (matrix_smm, matrix_smo)
220  CALL get_qs_env(qs_env, ks_env=ks_env)
221  CALL build_overlap_matrix_simple(ks_env, matrix_smm, &
222  mao_basis_set_list, mao_basis_set_list, smm_list)
223  CALL build_overlap_matrix_simple(ks_env, matrix_smo, &
224  mao_basis_set_list, orb_basis_set_list, smo_list)
225 
226  ! get reference density matrix and overlap matrix
227  IF (PRESENT(pmat_external)) THEN
228  matrix_p => pmat_external
229  ELSE
230  CALL get_qs_env(qs_env, rho=rho)
231  CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
232  END IF
233  IF (PRESENT(smat_external)) THEN
234  matrix_s => smat_external
235  ELSE
236  CALL get_qs_env(qs_env, matrix_s_kp=matrix_s)
237  END IF
238 
239  nspin = SIZE(matrix_p, 1)
240  eps_filter = 0.0_dp
241  ! Q matrix
242  CALL mao_build_q(matrix_q, matrix_p, matrix_s, matrix_smm, matrix_smo, smm_list, electra, eps_filter)
243 
244  ! MAO matrices
245  CALL get_qs_env(qs_env=qs_env, natom=natom)
246  CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, dbcsr_dist=dbcsr_dist)
247  NULLIFY (mao_coef)
248  CALL dbcsr_allocate_matrix_set(mao_coef, nspin)
249  ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
250  CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
251  basis=mao_basis_set_list)
252  IF (do_nmao_external) THEN
253  DO iatom = 1, natom
254  CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
255  kind_number=ikind)
256  col_blk_sizes(iatom) = nmao_external(ikind)
257  END DO
258  ELSE
259  CALL get_particle_set(particle_set, qs_kind_set, nmao=col_blk_sizes)
260  END IF
261 
262  ! check if MAOs have been specified
263  DO iab = 1, natom
264  IF (col_blk_sizes(iab) < 0) &
265  cpabort("Minimum number of MAOs has to be specified in KIND section for all elements")
266  END DO
267  DO ispin = 1, nspin
268  ! coefficients
269  ALLOCATE (mao_coef(ispin)%matrix)
270  CALL dbcsr_create(matrix=mao_coef(ispin)%matrix, &
271  name="MAO_COEF", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
272  row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes, nze=0)
273  CALL dbcsr_reserve_diag_blocks(matrix=mao_coef(ispin)%matrix)
274  END DO
275  DEALLOCATE (row_blk_sizes, col_blk_sizes)
276 
277  ! optimize MAOs
278  CALL mao_optimize(mao_coef, matrix_q, matrix_smm, electra, mao_max_iter, mao_eps_grad, eps1, &
279  iolev, iw)
280 
281  ! Deallocate the neighbor list structure
282  CALL release_neighbor_list_sets(smm_list)
283  CALL release_neighbor_list_sets(smo_list)
284 
285  DEALLOCATE (mao_basis_set_list, orb_basis_set_list)
286 
287  IF (ASSOCIATED(matrix_smm)) CALL dbcsr_deallocate_matrix_set(matrix_smm)
288  IF (ASSOCIATED(matrix_smo)) CALL dbcsr_deallocate_matrix_set(matrix_smo)
289  IF (ASSOCIATED(matrix_q)) CALL dbcsr_deallocate_matrix_set(matrix_q)
290 
291  CALL timestop(handle)
292 
293  END SUBROUTINE mao_generate_basis
294 
295 END MODULE mao_basis
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.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Calculate MAO's and analyze wavefunctions.
Definition: mao_basis.F:15
subroutine, public mao_generate_basis(qs_env, mao_coef, ref_basis_set, pmat_external, smat_external, molecular, max_iter, eps_grad, nmao_external, eps1_mao, iolevel, unit_nr)
...
Definition: mao_basis.F:75
Calculate MAO's and analyze wavefunctions.
Definition: mao_methods.F:15
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.
Definition: mao_methods.F:717
Calculate MAO's and analyze wavefunctions.
Definition: mao_optimizer.F:15
subroutine, public mao_optimize(mao_coef, matrix_q, matrix_smm, electra, max_iter, eps_grad, eps1, iolevel, iw)
...
Definition: mao_optimizer.F:55
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.
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.
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)
...
Definition: qs_ks_types.F:330
Define the neighbor list data types and the corresponding functionality.
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
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.
Definition: qs_overlap.F:19
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.
Definition: qs_overlap.F:558
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
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...
Definition: qs_rho_types.F:229