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