1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6MODULE optimize_dmfet_potential 7 USE cp_blacs_env, ONLY: cp_blacs_env_create,& 8 cp_blacs_env_release,& 9 cp_blacs_env_type 10 USE cp_control_types, ONLY: dft_control_type 11 USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr 12 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 13 cp_fm_struct_release,& 14 cp_fm_struct_type 15 USE cp_fm_types, ONLY: cp_fm_copy_general,& 16 cp_fm_create,& 17 cp_fm_maxabsval,& 18 cp_fm_release,& 19 cp_fm_set_all,& 20 cp_fm_type 21 USE cp_gemm_interface, ONLY: cp_gemm 22 USE cp_para_types, ONLY: cp_para_env_type 23 USE dbcsr_api, ONLY: dbcsr_create,& 24 dbcsr_init_p,& 25 dbcsr_multiply,& 26 dbcsr_p_type,& 27 dbcsr_release,& 28 dbcsr_trace,& 29 dbcsr_type,& 30 dbcsr_type_no_symmetry 31 USE embed_types, ONLY: opt_dmfet_pot_type 32 USE input_section_types, ONLY: section_vals_type,& 33 section_vals_val_get 34 USE kinds, ONLY: dp 35 USE qs_environment_types, ONLY: get_qs_env,& 36 qs_environment_type 37 USE qs_mo_types, ONLY: get_mo_set,& 38 mo_set_p_type 39#include "./base/base_uses.f90" 40 41 IMPLICIT NONE 42 43 PRIVATE 44 45 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'optimize_dmfet_potential' 46 47 PUBLIC :: build_full_dm, prepare_dmfet_opt, release_dmfet_opt, subsys_spin, check_dmfet 48 49CONTAINS 50 51! ************************************************************************************************** 52!> \brief ... 53!> \param opt_dmfet ... 54!> \param opt_dmfet_section ... 55! ************************************************************************************************** 56 SUBROUTINE read_opt_dmfet_section(opt_dmfet, opt_dmfet_section) 57 TYPE(opt_dmfet_pot_type) :: opt_dmfet 58 TYPE(section_vals_type), POINTER :: opt_dmfet_section 59 60 ! Read keywords 61 62 CALL section_vals_val_get(opt_dmfet_section, "N_ITER", & 63 i_val=opt_dmfet%n_iter) 64 65 CALL section_vals_val_get(opt_dmfet_section, "TRUST_RAD", & 66 r_val=opt_dmfet%trust_rad) 67 68 CALL section_vals_val_get(opt_dmfet_section, "DM_CONV_MAX", & 69 r_val=opt_dmfet%conv_max) 70 71 CALL section_vals_val_get(opt_dmfet_section, "DM_CONV_INT", & 72 r_val=opt_dmfet%conv_int) 73 74 CALL section_vals_val_get(opt_dmfet_section, "BETA_DM_CONV_MAX", & 75 r_val=opt_dmfet%conv_max_beta) 76 77 CALL section_vals_val_get(opt_dmfet_section, "BETA_DM_CONV_INT", & 78 r_val=opt_dmfet%conv_int_beta) 79 80 CALL section_vals_val_get(opt_dmfet_section, "READ_DMFET_POT", & 81 l_val=opt_dmfet%read_dmfet_pot) 82 83 END SUBROUTINE read_opt_dmfet_section 84 85! ************************************************************************************************** 86!> \brief ... 87!> \param qs_env ... 88!> \return ... 89! ************************************************************************************************** 90 FUNCTION subsys_spin(qs_env) RESULT(subsys_open_shell) 91 TYPE(qs_environment_type), POINTER :: qs_env 92 LOGICAL :: subsys_open_shell 93 94 TYPE(dft_control_type), POINTER :: dft_control 95 96 NULLIFY (dft_control) 97 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control) 98 subsys_open_shell = .FALSE. 99 IF (dft_control%nspins == 2) subsys_open_shell = .TRUE. 100 101 END FUNCTION subsys_spin 102 103! ************************************************************************************************** 104!> \brief ... 105!> \param qs_env ... 106!> \param opt_dmfet ... 107!> \param opt_dmfet_section ... 108! ************************************************************************************************** 109 SUBROUTINE prepare_dmfet_opt(qs_env, opt_dmfet, opt_dmfet_section) 110 TYPE(qs_environment_type), POINTER :: qs_env 111 TYPE(opt_dmfet_pot_type) :: opt_dmfet 112 TYPE(section_vals_type), POINTER :: opt_dmfet_section 113 114 INTEGER :: diff_size, nao 115 TYPE(cp_blacs_env_type), POINTER :: blacs_env 116 TYPE(cp_fm_struct_type), POINTER :: fm_struct 117 TYPE(cp_fm_type), POINTER :: mo_coeff 118 TYPE(cp_para_env_type), POINTER :: para_env 119 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 120 121 ! Read the input 122 CALL read_opt_dmfet_section(opt_dmfet, opt_dmfet_section) 123 124 ! Get the orbital coefficients 125 CALL get_qs_env(qs_env=qs_env, mos=mos, para_env=para_env) 126 CALL get_mo_set(mo_set=mos(1)%mo_set, mo_coeff=mo_coeff, nao=nao) 127 128 ! Make cp_fm matrices 129 NULLIFY (blacs_env) 130 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env) 131 132 NULLIFY (opt_dmfet%dmfet_pot, opt_dmfet%dm_1, opt_dmfet%dm_2, opt_dmfet%dm_total, opt_dmfet%dm_diff) 133 NULLIFY (fm_struct) 134 135 CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, & 136 nrow_global=nao, ncol_global=nao) 137 CALL cp_fm_create(opt_dmfet%dmfet_pot, fm_struct, name="dmfet_pot") 138 CALL cp_fm_create(opt_dmfet%dm_subsys, fm_struct, name="dm_subsys") 139 CALL cp_fm_create(opt_dmfet%dm_total, fm_struct, name="dm_total") 140 CALL cp_fm_create(opt_dmfet%dm_diff, fm_struct, name="dm_diff") 141 142 CALL cp_fm_set_all(opt_dmfet%dmfet_pot, 0.0_dp) 143 CALL cp_fm_set_all(opt_dmfet%dm_subsys, 0.0_dp) 144 CALL cp_fm_set_all(opt_dmfet%dm_total, 0.0_dp) 145 CALL cp_fm_set_all(opt_dmfet%dm_diff, 0.0_dp) 146 147 ! Beta spin counterparts 148 IF (opt_dmfet%open_shell_embed) THEN 149 CALL cp_fm_create(opt_dmfet%dmfet_pot_beta, fm_struct, name="dmfet_pot_beta") 150 CALL cp_fm_create(opt_dmfet%dm_subsys_beta, fm_struct, name="dm_subsys_beta") 151 CALL cp_fm_create(opt_dmfet%dm_total_beta, fm_struct, name="dm_total_beta") 152 CALL cp_fm_create(opt_dmfet%dm_diff_beta, fm_struct, name="dm_diff_beta") 153 154 CALL cp_fm_set_all(opt_dmfet%dmfet_pot_beta, 0.0_dp) 155 CALL cp_fm_set_all(opt_dmfet%dm_subsys_beta, 0.0_dp) 156 CALL cp_fm_set_all(opt_dmfet%dm_total_beta, 0.0_dp) 157 CALL cp_fm_set_all(opt_dmfet%dm_diff_beta, 0.0_dp) 158 ENDIF 159 160 ! Release structure and context 161 CALL cp_fm_struct_release(fm_struct) 162 CALL cp_blacs_env_release(blacs_env) 163 164 ! Array to store functional values 165 ALLOCATE (opt_dmfet%w_func(opt_dmfet%n_iter)) 166 opt_dmfet%w_func = 0.0_dp 167 168 ! Allocate max_diff and int_diff 169 diff_size = 1 170 IF (opt_dmfet%open_shell_embed) diff_size = 2 171 ALLOCATE (opt_dmfet%max_diff(diff_size)) 172 ALLOCATE (opt_dmfet%int_diff(diff_size)) 173 174 END SUBROUTINE prepare_dmfet_opt 175 176! ************************************************************************************************** 177!> \brief ... 178!> \param opt_dmfet ... 179! ************************************************************************************************** 180 SUBROUTINE release_dmfet_opt(opt_dmfet) 181 TYPE(opt_dmfet_pot_type) :: opt_dmfet 182 183 CALL cp_fm_release(opt_dmfet%dmfet_pot) 184 CALL cp_fm_release(opt_dmfet%dm_subsys) 185 CALL cp_fm_release(opt_dmfet%dm_total) 186 CALL cp_fm_release(opt_dmfet%dm_diff) 187 188 IF (opt_dmfet%open_shell_embed) THEN 189 CALL cp_fm_release(opt_dmfet%dmfet_pot_beta) 190 CALL cp_fm_release(opt_dmfet%dm_subsys_beta) 191 CALL cp_fm_release(opt_dmfet%dm_total_beta) 192 CALL cp_fm_release(opt_dmfet%dm_diff_beta) 193 ENDIF 194 195 DEALLOCATE (opt_dmfet%w_func) 196 DEALLOCATE (opt_dmfet%max_diff) 197 DEALLOCATE (opt_dmfet%int_diff) 198 DEALLOCATE (opt_dmfet%all_nspins) 199 200 END SUBROUTINE release_dmfet_opt 201 202! ************************************************************************************************** 203!> \brief Builds density matrices from MO coefficients in full matrix format 204!> \param qs_env ... 205!> \param dm ... 206!> \param open_shell Subsystem is open shell 207!> \param open_shell_embed Embedding is open shell 208!> \param dm_beta ... 209! ************************************************************************************************** 210 SUBROUTINE build_full_dm(qs_env, dm, open_shell, open_shell_embed, dm_beta) 211 TYPE(qs_environment_type), POINTER :: qs_env 212 TYPE(cp_fm_type), POINTER :: dm 213 LOGICAL :: open_shell, open_shell_embed 214 TYPE(cp_fm_type), POINTER :: dm_beta 215 216 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_full_dm', routineP = moduleN//':'//routineN 217 218 INTEGER :: handle, homo, nao 219 REAL(KIND=dp) :: coeff 220 TYPE(cp_fm_type), POINTER :: mo_coeff 221 TYPE(cp_para_env_type), POINTER :: para_env 222 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 223 224 CALL timeset(routineN, handle) 225 226 coeff = 2.0_dp 227 IF (open_shell_embed) coeff = 1.0_dp 228 229 ! Get the orbital coefficients 230 CALL get_qs_env(qs_env=qs_env, mos=mos) 231 CALL get_mo_set(mo_set=mos(1)%mo_set, mo_coeff=mo_coeff, nao=nao, homo=homo) 232 233 ! Build the density matrix 234 CALL cp_gemm(transa="N", transb="T", m=nao, n=nao, & 235 k=homo, alpha=coeff, & 236 matrix_a=mo_coeff, matrix_b=mo_coeff, & 237 beta=0.0_dp, matrix_c=dm) 238 239 ! Open shell 240 IF (open_shell) THEN 241 CALL get_mo_set(mo_set=mos(2)%mo_set, mo_coeff=mo_coeff, nao=nao, homo=homo) 242 243 ! Build the density matrix 244 CALL cp_gemm(transa="N", transb="T", m=nao, n=nao, & 245 k=homo, alpha=coeff, & 246 matrix_a=mo_coeff, matrix_b=mo_coeff, & 247 beta=0.0_dp, matrix_c=dm_beta) 248 ENDIF 249 250 ! If embedding is open shell, but subsystem is not, copy alpha-spin DM intou beta-spin DM 251 IF (open_shell_embed .AND. (.NOT. open_shell)) THEN 252 CALL get_qs_env(qs_env=qs_env, para_env=para_env) 253 CALL cp_fm_copy_general(dm, dm_beta, para_env) 254 255 ENDIF 256 257 CALL timestop(handle) 258 259 END SUBROUTINE build_full_dm 260 261! ************************************************************************************************** 262!> \brief ... 263!> \param opt_dmfet ... 264!> \param qs_env ... 265! ************************************************************************************************** 266 SUBROUTINE check_dmfet(opt_dmfet, qs_env) 267 TYPE(opt_dmfet_pot_type) :: opt_dmfet 268 TYPE(qs_environment_type), POINTER :: qs_env 269 270 REAL(KIND=dp) :: max_diff, max_diff_beta, trace 271 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 272 TYPE(dbcsr_type), POINTER :: diff_dbcsr, dm_s 273 274 CALL get_qs_env(qs_env, matrix_s=matrix_s) 275 276 NULLIFY (diff_dbcsr) 277 CALL dbcsr_init_p(diff_dbcsr) 278 CALL dbcsr_create(matrix=diff_dbcsr, & 279 template=matrix_s(1)%matrix, & 280 matrix_type=dbcsr_type_no_symmetry) 281 282 NULLIFY (dm_s) 283 CALL dbcsr_init_p(dm_s) 284 CALL dbcsr_create(matrix=dm_s, & 285 template=matrix_s(1)%matrix, & 286 matrix_type=dbcsr_type_no_symmetry) 287 288 CALL copy_fm_to_dbcsr(opt_dmfet%dm_diff, diff_dbcsr, keep_sparsity=.FALSE.) 289 290 CALL dbcsr_multiply("N", "N", 1.0_dp, diff_dbcsr, matrix_s(1)%matrix, & 291 0.0_dp, dm_s) 292 293 CALL dbcsr_trace(dm_s, trace) 294 295 IF (opt_dmfet%open_shell_embed) THEN 296 CALL copy_fm_to_dbcsr(opt_dmfet%dm_diff_beta, diff_dbcsr, keep_sparsity=.FALSE.) 297 298 CALL dbcsr_multiply("N", "N", 1.0_dp, diff_dbcsr, matrix_s(1)%matrix, & 299 0.0_dp, dm_s) 300 301 CALL dbcsr_trace(dm_s, trace) 302 303 ENDIF 304 305 ! Release dbcsr objects 306 CALL dbcsr_release(diff_dbcsr) 307 CALL dbcsr_release(dm_s) 308 309 ! Find the absolute maximum value of the DM difference 310 CALL cp_fm_maxabsval(opt_dmfet%dm_diff, max_diff) 311 IF (opt_dmfet%open_shell_embed) THEN 312 CALL cp_fm_maxabsval(opt_dmfet%dm_diff_beta, max_diff_beta) 313 ENDIF 314 315 END SUBROUTINE check_dmfet 316 317END MODULE optimize_dmfet_potential 318