|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00013 MODULE qs_kpp1_env_methods 00014 USE cell_types, ONLY: cell_type 00015 USE cp_dbcsr_interface, ONLY: cp_dbcsr_add,& 00016 cp_dbcsr_copy,& 00017 cp_dbcsr_init,& 00018 cp_dbcsr_scale,& 00019 cp_dbcsr_set 00020 USE cp_dbcsr_operations, ONLY: cp_dbcsr_allocate_matrix_set,& 00021 cp_dbcsr_deallocate_matrix_set 00022 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type 00023 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 00024 cp_print_key_should_output,& 00025 cp_print_key_unit_nr 00026 USE f77_blas 00027 USE hartree_local_methods, ONLY: Vh_1c_gg_integrals 00028 USE input_constants, ONLY: do_method_gapw,& 00029 do_method_gapw_xc,& 00030 tddfpt_excitations,& 00031 tddfpt_singlet,& 00032 tddfpt_triplet 00033 USE input_section_types, ONLY: section_get_ival,& 00034 section_get_rval,& 00035 section_vals_get,& 00036 section_vals_get_subs_vals,& 00037 section_vals_type,& 00038 section_vals_val_get 00039 USE kahan_sum, ONLY: accurate_sum 00040 USE kinds, ONLY: dp 00041 USE pw_env_types, ONLY: pw_env_get,& 00042 pw_env_type 00043 USE pw_methods, ONLY: pw_axpy,& 00044 pw_copy,& 00045 pw_integrate_function,& 00046 pw_transfer,& 00047 pw_zero 00048 USE pw_poisson_methods, ONLY: pw_poisson_solve 00049 USE pw_poisson_types, ONLY: pw_poisson_type 00050 USE pw_pool_types, ONLY: pw_pool_create_pw,& 00051 pw_pool_give_back_pw,& 00052 pw_pool_p_type,& 00053 pw_pool_type 00054 USE pw_types, ONLY: COMPLEXDATA1D,& 00055 REALDATA3D,& 00056 REALSPACE,& 00057 RECIPROCALSPACE,& 00058 pw_create,& 00059 pw_p_type,& 00060 pw_release,& 00061 pw_retain 00062 USE qs_energy_types, ONLY: allocate_qs_energy,& 00063 deallocate_qs_energy,& 00064 qs_energy_type 00065 USE qs_environment_types, ONLY: get_qs_env,& 00066 qs_environment_type 00067 USE qs_gapw_densities, ONLY: prepare_gapw_den 00068 USE qs_integrate_potential, ONLY: integrate_v_rspace 00069 USE qs_kpp1_env_types, ONLY: qs_kpp1_env_type 00070 USE qs_ks_atom, ONLY: update_ks_atom 00071 USE qs_ks_methods, ONLY: qs_ks_build_kohn_sham_matrix 00072 USE qs_p_env_types, ONLY: qs_p_env_type 00073 USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace 00074 USE qs_rho_types, ONLY: qs_rho_get,& 00075 qs_rho_type 00076 USE qs_vxc_atom, ONLY: calculate_xc_2nd_deriv_atom 00077 USE timings, ONLY: timeset,& 00078 timestop 00079 USE xc, ONLY: xc_calc_2nd_deriv,& 00080 xc_prep_2nd_deriv 00081 USE xc_derivative_set_types, ONLY: xc_dset_release 00082 USE xc_derivatives, ONLY: xc_functionals_get_needs 00083 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 00084 USE xc_rho_set_types, ONLY: xc_rho_set_create,& 00085 xc_rho_set_release,& 00086 xc_rho_set_type,& 00087 xc_rho_set_update 00088 #include "cp_common_uses.h" 00089 00090 IMPLICIT NONE 00091 00092 PRIVATE 00093 00094 LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE. 00095 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_kpp1_env_methods' 00096 INTEGER, PRIVATE, SAVE :: last_kpp1_id_nr=0 00097 00098 PUBLIC :: kpp1_create, & 00099 kpp1_calc_k_p_p1, & 00100 kpp1_calc_k_p_p1_fdiff, & 00101 kpp1_did_change 00102 00103 CONTAINS 00104 00105 ! ***************************************************************************** 00115 SUBROUTINE kpp1_create(kpp1_env,qs_env, error) 00116 TYPE(qs_kpp1_env_type), POINTER :: kpp1_env 00117 TYPE(qs_environment_type), POINTER :: qs_env 00118 TYPE(cp_error_type), INTENT(inout) :: error 00119 00120 CHARACTER(len=*), PARAMETER :: routineN = 'kpp1_create', 00121 routineP = moduleN//':'//routineN 00122 00123 INTEGER :: stat 00124 LOGICAL :: failure 00125 00126 failure=.FALSE. 00127 00128 ALLOCATE(kpp1_env,stat=stat) 00129 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00130 IF (.NOT. failure) THEN 00131 NULLIFY(kpp1_env%v_rspace, kpp1_env%v_ao, kpp1_env%drho_r,& 00132 kpp1_env%rho_set, & 00133 kpp1_env%deriv_set, kpp1_env%spin_pot, kpp1_env%grad_pot,& 00134 kpp1_env%ndiag_term) 00135 kpp1_env%ref_count=1 00136 last_kpp1_id_nr=last_kpp1_id_nr+1 00137 kpp1_env%id_nr=last_kpp1_id_nr 00138 kpp1_env%iter=0 00139 kpp1_env%print_count=0 00140 END IF 00141 END SUBROUTINE kpp1_create 00142 00143 ! ***************************************************************************** 00155 SUBROUTINE kpp1_calc_k_p_p1(kpp1_env, p_env, qs_env, k_p_p1, rho, rho1, rho1_xc, error) 00156 00157 TYPE(qs_kpp1_env_type), POINTER :: kpp1_env 00158 TYPE(qs_p_env_type), POINTER :: p_env 00159 TYPE(qs_environment_type), POINTER :: qs_env 00160 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00161 POINTER :: k_p_p1 00162 TYPE(qs_rho_type), POINTER :: rho, rho1 00163 TYPE(qs_rho_type), OPTIONAL, POINTER :: rho1_xc 00164 TYPE(cp_error_type), INTENT(inout) :: error 00165 00166 CHARACTER(len=*), PARAMETER :: routineN = 'kpp1_calc_k_p_p1', 00167 routineP = moduleN//':'//routineN 00168 00169 INTEGER :: excitations, handle, ispin, 00170 nspins, output_unit, 00171 res_etype, stat 00172 INTEGER, DIMENSION(2, 3) :: bo 00173 LOGICAL :: explicit, failure, gapw, 00174 gapw_xc, ionode, lsd, 00175 lsd_singlets 00176 REAL(KIND=dp) :: energy_hartree, 00177 energy_hartree_1c, fac 00178 TYPE(cell_type), POINTER :: cell 00179 TYPE(cp_logger_type), POINTER :: logger 00180 TYPE(pw_env_type), POINTER :: pw_env 00181 TYPE(pw_p_type) :: rho1_tot_gspace, 00182 v_hartree_gspace, 00183 v_hartree_rspace 00184 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho1_g_pw, rho1_r, rho1_r_pw, 00185 tau_pw, v_rspace_new, v_xc 00186 TYPE(pw_poisson_type), POINTER :: poisson_env 00187 TYPE(pw_pool_p_type), DIMENSION(:), 00188 POINTER :: pw_pools 00189 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00190 TYPE(section_vals_type), POINTER :: input, scf_section, 00191 xc_fun_section, xc_section 00192 TYPE(xc_rho_cflags_type) :: needs 00193 TYPE(xc_rho_set_type), POINTER :: rho1_set 00194 00195 CALL timeset(routineN,handle) 00196 failure=.FALSE. 00197 00198 NULLIFY(auxbas_pw_pool, pw_pools, pw_env, v_rspace_new, & 00199 rho1_r, rho1_g_pw, tau_pw, cell, v_xc, rho1_set,& 00200 poisson_env, input, scf_section) 00201 logger => cp_error_get_logger(error) 00202 ionode = logger%para_env%mepos==logger%para_env%source 00203 00204 CPPrecondition(ASSOCIATED(kpp1_env),cp_failure_level,routineP,error,failure) 00205 CPPrecondition(ASSOCIATED(k_p_p1),cp_failure_level,routineP,error,failure) 00206 CPPrecondition(ASSOCIATED(rho),cp_failure_level,routineP,error,failure) 00207 CPPrecondition(ASSOCIATED(rho1),cp_failure_level,routineP,error,failure) 00208 00209 IF (.NOT.failure) THEN 00210 CPPrecondition(kpp1_env%ref_count>0,cp_failure_level,routineP,error,failure) 00211 CALL kpp1_check_i_alloc(kpp1_env,qs_env=qs_env,error=error) 00212 CALL get_qs_env(qs_env=qs_env,& 00213 pw_env=pw_env,& 00214 cell=cell,& 00215 input=input,& 00216 error=error) 00217 00218 gapw=(section_get_ival(input,"DFT%QS%METHOD",error=error)==do_method_gapw) 00219 gapw_xc=(section_get_ival(input,"DFT%QS%METHOD",error=error)==do_method_gapw_xc) 00220 IF(gapw_xc) THEN 00221 CPPrecondition(ASSOCIATED(rho1_xc),cp_failure_level,routineP,error,failure) 00222 END IF 00223 00224 nspins = SIZE(k_p_p1) 00225 lsd = (nspins==2) 00226 00227 xc_section => section_vals_get_subs_vals(input,"DFT%XC",error=error) 00228 scf_section => section_vals_get_subs_vals(input,"DFT%SCF",error=error) 00229 CALL section_vals_val_get(input,"DFT%EXCITATIONS",& 00230 i_val=excitations,error=error) 00231 IF (excitations==tddfpt_excitations) THEN 00232 xc_section => section_vals_get_subs_vals(input,"DFT%TDDFPT%XC",error=error) 00233 !FM this check should already had happened and section made explicit, give an error? 00234 CALL section_vals_get(xc_section,explicit=explicit,error=error) 00235 IF (.NOT.explicit) THEN 00236 xc_section => section_vals_get_subs_vals(input,"DFT%XC",error=error) 00237 END IF 00238 END IF 00239 00240 CALL section_vals_val_get(input,"DFT%TDDFPT%LSD_SINGLETS",& 00241 l_val=lsd_singlets,error=error) 00242 CALL section_vals_val_get(input,"DFT%TDDFPT%RES_ETYPE",& 00243 i_val=res_etype,error=error) 00244 END IF 00245 00246 IF (.NOT.failure) THEN 00247 kpp1_env%iter=kpp1_env%iter+1 00248 END IF 00249 00250 ! gets the tmp grids 00251 IF (.NOT. failure) THEN 00252 CPPrecondition(ASSOCIATED(pw_env),cp_failure_level,routineP,error,failure) 00253 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool,& 00254 pw_pools=pw_pools, poisson_env=poisson_env,error=error) 00255 ALLOCATE(v_rspace_new(nspins), stat=stat) 00256 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00257 END IF 00258 IF (.NOT.failure) THEN 00259 CALL pw_pool_create_pw(auxbas_pw_pool,v_hartree_gspace%pw,& 00260 use_data = COMPLEXDATA1D,& 00261 in_space = RECIPROCALSPACE, error=error) 00262 CALL pw_pool_create_pw(auxbas_pw_pool,v_hartree_rspace%pw,& 00263 use_data = REALDATA3D,& 00264 in_space = REALSPACE, error=error) 00265 END IF 00266 00267 IF (gapw .OR. gapw_xc) & 00268 CALL prepare_gapw_den(qs_env,p_env%local_rho_set, do_rho0=(.NOT.gapw_xc), error=error) 00269 00270 ! *** calculate the hartree potential on the total density *** 00271 IF (.NOT. failure) THEN 00272 00273 CALL pw_pool_create_pw(auxbas_pw_pool, rho1_tot_gspace%pw,& 00274 use_data = COMPLEXDATA1D,& 00275 in_space = RECIPROCALSPACE, error=error) 00276 00277 CALL pw_copy(rho1%rho_g(1)%pw,rho1_tot_gspace%pw,error=error) 00278 DO ispin=2,nspins 00279 CALL pw_axpy(rho1%rho_g(ispin)%pw, rho1_tot_gspace%pw, error=error) 00280 END DO 00281 IF (gapw) & 00282 CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rho0_s_gs%pw, rho1_tot_gspace%pw,& 00283 error=error) 00284 00285 IF (cp_print_key_should_output(logger%iter_info,scf_section,"PRINT%TOTAL_DENSITIES",& 00286 error=error)/=0) THEN 00287 output_unit = cp_print_key_unit_nr(logger,scf_section,"PRINT%TOTAL_DENSITIES",& 00288 extension=".scfLog",error=error) 00289 CALL print_densities(kpp1_env, rho1, rho1_tot_gspace, output_unit, error=error) 00290 CALL cp_print_key_finished_output(output_unit,logger,scf_section,& 00291 "PRINT%TOTAL_DENSITIES", error=error) 00292 END IF 00293 00294 IF (.NOT.(nspins==1 .AND. excitations==tddfpt_excitations .AND. & 00295 res_etype /= tddfpt_singlet )) THEN 00296 CALL pw_poisson_solve(poisson_env,rho1_tot_gspace%pw, & 00297 energy_hartree, & 00298 v_hartree_gspace%pw,error=error) 00299 CALL pw_transfer(v_hartree_gspace%pw,v_hartree_rspace%pw,error=error) 00300 END IF 00301 00302 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho1_tot_gspace%pw,& 00303 error=error) 00304 00305 ! *** calculate the xc potential *** 00306 IF(gapw_xc) THEN 00307 CALL qs_rho_get(rho1_xc, rho_r=rho1_r, error=error) 00308 ELSE 00309 CALL qs_rho_get(rho1, rho_r=rho1_r, error=error) 00310 END IF 00311 00312 IF (nspins == 1 .AND. excitations==tddfpt_excitations .AND. & 00313 (lsd_singlets .OR. res_etype == tddfpt_triplet)) THEN 00314 00315 lsd = .TRUE. 00316 ALLOCATE(rho1_r_pw(2)) 00317 DO ispin=1, 2 00318 NULLIFY(rho1_r_pw(ispin)%pw) 00319 CALL pw_create(rho1_r_pw(ispin)%pw, rho1_r(1)%pw%pw_grid, & 00320 rho1_r(1)%pw%in_use, rho1_r(1)%pw%in_space,error=error) 00321 CALL pw_transfer(rho1_r(1)%pw, rho1_r_pw(ispin)%pw,error=error) 00322 END DO 00323 00324 ELSE 00325 00326 ALLOCATE(rho1_r_pw(nspins)) 00327 DO ispin=1, nspins 00328 rho1_r_pw(ispin)%pw => rho1_r(ispin)%pw 00329 CALL pw_retain(rho1_r_pw(ispin)%pw,error=error) 00330 END DO 00331 00332 END IF 00333 00334 NULLIFY(tau_pw) 00335 00336 !------! 00337 ! rho1 ! 00338 !------! 00339 bo = rho1_r(1)%pw%pw_grid%bounds_local 00340 ! create the place where to store the argument for the functionals 00341 CALL xc_rho_set_create(rho1_set, bo, & 00342 rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF",error=error), & 00343 drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF",error=error), & 00344 tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF",error=error), & 00345 error=error) 00346 00347 xc_fun_section => section_vals_get_subs_vals(xc_section,"XC_FUNCTIONAL",& 00348 error=error) 00349 needs=xc_functionals_get_needs(xc_fun_section,lsd,.TRUE.,error) 00350 00351 ! calculate the arguments needed by the functionals 00352 CALL xc_rho_set_update(rho1_set, rho1_r_pw, rho1_g_pw, tau_pw, needs,& 00353 section_get_ival(xc_section,"XC_GRID%XC_DERIV",error=error),& 00354 section_get_ival(xc_section,"XC_GRID%XC_SMOOTH_RHO",error=error),& 00355 cell, auxbas_pw_pool, error) 00356 00357 ALLOCATE(v_xc(nspins),stat=stat) 00358 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00359 DO ispin=1, nspins 00360 NULLIFY(v_xc(ispin)%pw) 00361 CALL pw_pool_create_pw(auxbas_pw_pool,v_xc(ispin)%pw,& 00362 use_data = REALDATA3D,& 00363 in_space = REALSPACE, error=error) 00364 CALL pw_zero(v_xc(ispin)%pw,error=error) 00365 END DO 00366 00367 fac=0._dp 00368 IF (nspins==1.AND.excitations==tddfpt_excitations) THEN 00369 IF (lsd_singlets) fac=1.0_dp 00370 IF (res_etype == tddfpt_triplet) fac=-1.0_dp 00371 END IF 00372 00373 CALL xc_calc_2nd_deriv(v_xc, kpp1_env%deriv_set, kpp1_env%rho_set, & 00374 rho1_set, auxbas_pw_pool,xc_section=xc_section,& 00375 tddfpt_fac=fac, error=error) 00376 00377 DO ispin=1,nspins 00378 v_rspace_new(ispin)%pw => v_xc(ispin)%pw 00379 END DO 00380 DEALLOCATE(v_xc,stat=stat) 00381 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00382 00383 IF (gapw) CALL calculate_xc_2nd_deriv_atom(p_env, qs_env, xc_section, error=error) 00384 00385 CALL xc_rho_set_release(rho1_set,error=error) 00386 00387 DO ispin=1,SIZE(rho1_r_pw) 00388 CALL pw_release(rho1_r_pw(ispin)%pw,error=error) 00389 END DO 00390 DEALLOCATE(rho1_r_pw, stat=stat) 00391 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00392 00393 !-------------------------------! 00394 ! Add both hartree and xc terms ! 00395 !-------------------------------! 00396 DO ispin=1,nspins 00397 00398 IF(gapw_xc) THEN 00399 ! XC and Hartree are integrated separatedly 00400 ! XC uses the sofft basis set only 00401 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d * & 00402 v_rspace_new(ispin)%pw%pw_grid%dvol 00403 00404 IF (excitations==tddfpt_excitations .AND. nspins==1) THEN 00405 00406 IF (.NOT.(lsd_singlets .OR. & 00407 res_etype == tddfpt_triplet)) THEN 00408 00409 v_rspace_new(1)%pw%cr3d = 2.0_dp * v_rspace_new(1)%pw%cr3d 00410 00411 END IF 00412 ! remove kpp1_env%v_ao and work directly on k_p_p1 ? 00413 CALL cp_dbcsr_set(kpp1_env%v_ao(ispin)%matrix,0.0_dp,error=error) 00414 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin),& 00415 p=rho%rho_ao(ispin),& 00416 h=kpp1_env%v_ao(ispin),& 00417 qs_env=qs_env,& 00418 calculate_forces=.FALSE.,gapw=gapw_xc,error=error) 00419 00420 ! add hartree only for SINGLETS 00421 IF (res_etype == tddfpt_singlet) THEN 00422 v_hartree_rspace%pw%cr3d = v_hartree_rspace%pw%cr3d * & 00423 v_hartree_rspace%pw%pw_grid%dvol 00424 v_rspace_new(1)%pw%cr3d = 2.0_dp * v_hartree_rspace%pw%cr3d 00425 00426 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin),& 00427 p=rho%rho_ao(ispin),& 00428 h=kpp1_env%v_ao(ispin),& 00429 qs_env=qs_env,& 00430 calculate_forces=.FALSE.,gapw=gapw,error=error) 00431 END IF 00432 ELSE 00433 ! remove kpp1_env%v_ao and work directly on k_p_p1 ? 00434 CALL cp_dbcsr_set(kpp1_env%v_ao(ispin)%matrix,0.0_dp,error=error) 00435 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin),& 00436 p=rho%rho_ao(ispin),& 00437 h=kpp1_env%v_ao(ispin),& 00438 qs_env=qs_env,& 00439 calculate_forces=.FALSE.,gapw=gapw_xc,error=error) 00440 00441 IF (ispin == 1) THEN 00442 v_hartree_rspace%pw%cr3d = v_hartree_rspace%pw%cr3d * & 00443 v_hartree_rspace%pw%pw_grid%dvol 00444 END IF 00445 v_rspace_new(ispin)%pw%cr3d = v_hartree_rspace%pw%cr3d 00446 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin),& 00447 p=rho%rho_ao(ispin),& 00448 h=kpp1_env%v_ao(ispin),& 00449 qs_env=qs_env,& 00450 calculate_forces=.FALSE.,gapw=gapw,error=error) 00451 END IF 00452 00453 ELSE 00454 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d * & 00455 v_rspace_new(ispin)%pw%pw_grid%dvol 00456 00457 IF (excitations==tddfpt_excitations .AND. nspins==1) THEN 00458 00459 IF (.NOT.(lsd_singlets .OR. & 00460 res_etype == tddfpt_triplet)) THEN 00461 00462 v_rspace_new(1)%pw%cr3d = 2.0_dp * v_rspace_new(1)%pw%cr3d 00463 00464 END IF 00465 00466 ! add hartree only for SINGLETS 00467 IF (res_etype == tddfpt_singlet) THEN 00468 v_hartree_rspace%pw%cr3d = v_hartree_rspace%pw%cr3d * & 00469 v_hartree_rspace%pw%pw_grid%dvol 00470 v_rspace_new(1)%pw%cr3d = v_rspace_new(1)%pw%cr3d + & 00471 2.0_dp * v_hartree_rspace%pw%cr3d 00472 END IF 00473 ELSE 00474 IF (ispin == 1) THEN 00475 v_hartree_rspace%pw%cr3d = v_hartree_rspace%pw%cr3d * & 00476 v_hartree_rspace%pw%pw_grid%dvol 00477 END IF 00478 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d + & 00479 v_hartree_rspace%pw%cr3d 00480 END IF 00481 00482 ! remove kpp1_env%v_ao and work directly on k_p_p1 ? 00483 CALL cp_dbcsr_set(kpp1_env%v_ao(ispin)%matrix,0.0_dp,error=error) 00484 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin),& 00485 p=rho%rho_ao(ispin),& 00486 h=kpp1_env%v_ao(ispin),& 00487 qs_env=qs_env,& 00488 calculate_forces=.FALSE.,gapw=gapw,error=error) 00489 00490 END IF 00491 00492 CALL cp_dbcsr_copy(k_p_p1(ispin)%matrix,kpp1_env%v_ao(ispin)%matrix,error=error) 00493 END DO 00494 00495 IF (gapw) THEN 00496 IF (.NOT. (excitations==tddfpt_excitations .AND. (nspins == 1 .AND. & 00497 res_etype == tddfpt_triplet))) THEN 00498 CALL Vh_1c_gg_integrals(qs_env,energy_hartree_1c, .TRUE., p_env=p_env,error=error) 00499 CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, & 00500 .FALSE., .TRUE., p_env=p_env, error=error) 00501 END IF 00502 ! *** Add single atom contributions to the KS matrix *** 00503 CALL update_ks_atom(qs_env,k_p_p1,rho%rho_ao,.FALSE.,.TRUE.,p_env,error) 00504 END IF 00505 00506 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_hartree_gspace%pw,& 00507 error=error) 00508 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_hartree_rspace%pw,& 00509 error=error) 00510 DO ispin=1,nspins 00511 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_rspace_new(ispin)%pw,& 00512 error=error) 00513 END DO 00514 DEALLOCATE(v_rspace_new, stat=stat) 00515 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00516 END IF 00517 00518 CALL timestop(handle) 00519 END SUBROUTINE kpp1_calc_k_p_p1 00520 00521 ! ***************************************************************************** 00541 SUBROUTINE kpp1_calc_k_p_p1_fdiff(kpp1_env,qs_env,k_p_p1,rho,rho1,& 00542 diff, error) 00543 TYPE(qs_kpp1_env_type), POINTER :: kpp1_env 00544 TYPE(qs_environment_type), POINTER :: qs_env 00545 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00546 POINTER :: k_p_p1 00547 TYPE(qs_rho_type), POINTER :: rho, rho1 00548 REAL(KIND=dp), INTENT(in), OPTIONAL :: diff 00549 TYPE(cp_error_type), INTENT(inout) :: error 00550 00551 CHARACTER(len=*), PARAMETER :: routineN = 'kpp1_calc_k_p_p1_fdiff', 00552 routineP = moduleN//':'//routineN 00553 00554 INTEGER :: ispin, nspins 00555 REAL(KIND=dp) :: my_diff 00556 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00557 POINTER :: ks_2, matrix_s, rho_ao 00558 TYPE(qs_energy_type), POINTER :: qs_energy 00559 00560 NULLIFY(ks_2,matrix_s,qs_energy) 00561 nspins = SIZE(k_p_p1) 00562 my_diff=1.0e-6_dp 00563 IF (PRESENT(diff)) my_diff=diff 00564 CALL allocate_qs_energy(qs_energy) 00565 00566 CALL qs_rho_get(rho, rho_ao=rho_ao,error=error) 00567 CALL get_qs_env(qs_env, & 00568 matrix_s=matrix_s,error=error) 00569 00570 ! rho = rho0+h/2*rho1 00571 my_diff=my_diff/2.0_dp 00572 DO ispin=1,SIZE(k_p_p1) 00573 CALL cp_dbcsr_add(rho%rho_ao(ispin)%matrix,rho1%rho_ao(ispin)%matrix,& 00574 alpha_scalar=1.0_dp,beta_scalar=my_diff,error=error) 00575 rho%rho_r(ispin)%pw%cr3d=rho%rho_r(ispin)%pw%cr3d+& 00576 my_diff*rho1%rho_r(ispin)%pw%cr3d 00577 rho%rho_g(ispin)%pw%cc=rho%rho_g(ispin)%pw%cc+& 00578 my_diff*rho1%rho_g(ispin)%pw%cc 00579 END DO 00580 00581 CALL qs_ks_build_kohn_sham_matrix(ks_env=qs_env%ks_env,& 00582 qs_env=qs_env,ks_matrix=k_p_p1,rho=rho,energy=qs_energy,& 00583 calculate_forces=.FALSE.,& 00584 just_energy=.FALSE.,error=error) 00585 00586 CALL cp_dbcsr_allocate_matrix_set(ks_2,nspins,error=error) 00587 DO ispin=1,nspins 00588 ALLOCATE(ks_2(ispin)%matrix) 00589 CALL cp_dbcsr_init(ks_2(ispin)%matrix,error=error) 00590 CALL cp_dbcsr_copy(ks_2(ispin)%matrix,matrix_s(1)%matrix,& 00591 name="tmp_ks2-"//ADJUSTL(cp_to_string(ispin)),error=error) 00592 END DO 00593 00594 ! rho = rho0-h/2*rho1 00595 my_diff=-2.0_dp*my_diff 00596 DO ispin=1,nspins 00597 CALL cp_dbcsr_add(rho%rho_ao(ispin)%matrix,rho1%rho_ao(ispin)%matrix,& 00598 alpha_scalar=1.0_dp,beta_scalar=my_diff,error=error) 00599 rho%rho_r(ispin)%pw%cr3d=rho%rho_r(ispin)%pw%cr3d+& 00600 my_diff*rho1%rho_r(ispin)%pw%cr3d 00601 rho%rho_g(ispin)%pw%cc=rho%rho_g(ispin)%pw%cc+& 00602 my_diff*rho1%rho_g(ispin)%pw%cc 00603 END DO 00604 00605 CALL qs_ks_build_kohn_sham_matrix(ks_env=qs_env%ks_env,& 00606 qs_env=qs_env,ks_matrix=ks_2,rho=rho,energy=qs_energy,& 00607 calculate_forces=.FALSE.,& 00608 just_energy=.FALSE.,error=error) 00609 00610 ! rho = rho0 00611 my_diff=-0.5_dp*my_diff 00612 DO ispin=1,nspins 00613 CALL cp_dbcsr_add(rho%rho_ao(ispin)%matrix,rho1%rho_ao(ispin)%matrix,& 00614 alpha_scalar=1.0_dp,beta_scalar=my_diff,error=error) 00615 rho%rho_r(ispin)%pw%cr3d=rho%rho_r(ispin)%pw%cr3d+& 00616 my_diff*rho1%rho_r(ispin)%pw%cr3d 00617 rho%rho_g(ispin)%pw%cc=rho%rho_g(ispin)%pw%cc+& 00618 my_diff*rho1%rho_g(ispin)%pw%cc 00619 END DO 00620 00621 ! k_p_p1=(H(rho0+h/2 rho1)-H(rho0-h/2 rho1))/h 00622 DO ispin=1,nspins 00623 CALL cp_dbcsr_add(k_p_p1(ispin)%matrix,ks_2(ispin)%matrix,& 00624 alpha_scalar=1.0_dp,beta_scalar=-1.0_dp,error=error) 00625 CALL cp_dbcsr_scale(k_p_p1(ispin)%matrix,alpha_scalar=0.5_dp/my_diff,error=error) 00626 END DO 00627 00628 CALL cp_dbcsr_deallocate_matrix_set(ks_2,error=error) 00629 CALL deallocate_qs_energy(qs_energy) 00630 END SUBROUTINE kpp1_calc_k_p_p1_fdiff 00631 00632 ! ***************************************************************************** 00642 SUBROUTINE kpp1_check_i_alloc(kpp1_env, qs_env, error) 00643 00644 TYPE(qs_kpp1_env_type), POINTER :: kpp1_env 00645 TYPE(qs_environment_type), POINTER :: qs_env 00646 TYPE(cp_error_type), INTENT(inout) :: error 00647 00648 CHARACTER(len=*), PARAMETER :: routineN = 'kpp1_check_i_alloc', 00649 routineP = moduleN//':'//routineN 00650 00651 INTEGER :: excitations, ispin, nspins, 00652 res_etype, stat 00653 LOGICAL :: explicit, failure, 00654 lsd_singlets 00655 TYPE(cell_type), POINTER :: cell 00656 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00657 POINTER :: matrix_s 00658 TYPE(pw_env_type), POINTER :: pw_env 00659 TYPE(pw_p_type), DIMENSION(:), POINTER :: my_rho_r, rho_r 00660 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00661 TYPE(qs_rho_type), POINTER :: rho 00662 TYPE(section_vals_type), POINTER :: input, xc_section 00663 00664 ! ------------------------------------------------------------------ 00665 00666 failure=.FALSE. 00667 NULLIFY(pw_env,auxbas_pw_pool,matrix_s,rho,rho_r,input) 00668 00669 CPPrecondition(ASSOCIATED(kpp1_env),cp_failure_level,routineP,error,failure) 00670 CPPrecondition(kpp1_env%ref_count>0,cp_failure_level,routineP,error,failure) 00671 00672 IF (.NOT. failure) THEN 00673 00674 CALL get_qs_env(qs_env, pw_env=pw_env,& 00675 matrix_s=matrix_s, cell=cell, input=input, error=error, rho=rho) 00676 00677 CALL qs_rho_get(rho, rho_r=rho_r, error=error) 00678 nspins=SIZE(rho_r) 00679 00680 CALL pw_env_get(pw_env, auxbas_pw_pool = auxbas_pw_pool, error=error) 00681 00682 IF (.NOT.ASSOCIATED(kpp1_env%v_rspace)) THEN 00683 ALLOCATE(kpp1_env%v_rspace(nspins),stat=stat) 00684 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00685 IF (.NOT.failure) THEN 00686 DO ispin=1,nspins 00687 CALL pw_pool_create_pw(auxbas_pw_pool, & 00688 kpp1_env%v_rspace(ispin)%pw,& 00689 use_data=REALDATA3D, in_space=REALSPACE,error=error) 00690 END DO 00691 END IF 00692 END IF 00693 00694 IF (.NOT.ASSOCIATED(kpp1_env%v_ao)) THEN 00695 CALL cp_dbcsr_allocate_matrix_set(kpp1_env%v_ao,nspins,error=error) 00696 DO ispin=1,nspins 00697 ALLOCATE(kpp1_env%v_ao(ispin)%matrix) 00698 CALL cp_dbcsr_init(kpp1_env%v_ao(ispin)%matrix,error=error) 00699 CALL cp_dbcsr_copy(kpp1_env%v_ao(ispin)%matrix,matrix_s(1)%matrix,& 00700 name="kpp1%v_ao-"//ADJUSTL(cp_to_string(ispin)),error=error) 00701 END DO 00702 END IF 00703 00704 IF (.not.ASSOCIATED(kpp1_env%deriv_set)) THEN 00705 00706 CALL section_vals_val_get(input,"DFT%EXCITATIONS",& 00707 i_val=excitations,error=error) 00708 CALL section_vals_val_get(input,"DFT%TDDFPT%LSD_SINGLETS",& 00709 l_val=lsd_singlets,error=error) 00710 CALL section_vals_val_get(input,"DFT%TDDFPT%RES_ETYPE",& 00711 i_val=res_etype,error=error) 00712 IF (nspins==1.AND.(excitations==tddfpt_excitations.and.& 00713 (lsd_singlets .OR. res_etype == tddfpt_triplet))) THEN 00714 ALLOCATE(my_rho_r(2),stat=stat) 00715 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00716 DO ispin=1,2 00717 CALL pw_pool_create_pw(auxbas_pw_pool,my_rho_r(ispin)%pw, & 00718 use_data=rho_r(1)%pw%in_use, in_space=rho_r(1)%pw%in_space,& 00719 error=error) 00720 my_rho_r(ispin)%pw%cr3d = 0.5_dp * rho_r(1)%pw%cr3d 00721 END DO 00722 ELSE 00723 ALLOCATE(my_rho_r(SIZE(rho_r)),stat=stat) 00724 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00725 DO ispin=1,SIZE(rho_r) 00726 my_rho_r(ispin)%pw => rho_r(ispin)%pw 00727 CALL pw_retain(my_rho_r(ispin)%pw,error=error) 00728 END DO 00729 END IF 00730 00731 xc_section => section_vals_get_subs_vals(input,"DFT%XC",error=error) 00732 CALL section_vals_val_get(input,"DFT%EXCITATIONS",& 00733 i_val=excitations,error=error) 00734 IF (excitations==tddfpt_excitations) THEN 00735 xc_section => section_vals_get_subs_vals(input,"DFT%TDDFPT%XC",error=error) 00736 !FM this check should already had happened and section made explicit, give an error? 00737 CALL section_vals_get(xc_section,explicit=explicit,error=error) 00738 IF (.NOT.explicit) THEN 00739 xc_section => section_vals_get_subs_vals(input,"DFT%XC",error=error) 00740 END IF 00741 END IF 00742 00743 CALL xc_prep_2nd_deriv(kpp1_env%deriv_set, kpp1_env%rho_set, & 00744 my_rho_r, auxbas_pw_pool, & 00745 xc_section=xc_section, cell=cell, error=error) 00746 00747 DO ispin=1,SIZE(my_rho_r) 00748 CALL pw_release(my_rho_r(ispin)%pw,error=error) 00749 END DO 00750 DEALLOCATE(my_rho_r,stat=stat) 00751 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00752 END IF 00753 END IF 00754 00755 END SUBROUTINE kpp1_check_i_alloc 00756 00757 ! ***************************************************************************** 00770 SUBROUTINE kpp1_did_change(kpp1_env, s_struct_changed, grid_changed,& 00771 psi0_changed, error) 00772 TYPE(qs_kpp1_env_type), POINTER :: kpp1_env 00773 LOGICAL, INTENT(in), OPTIONAL :: s_struct_changed, 00774 grid_changed, psi0_changed 00775 TYPE(cp_error_type), INTENT(inout) :: error 00776 00777 CHARACTER(len=*), PARAMETER :: routineN = 'kpp1_did_change', 00778 routineP = moduleN//':'//routineN 00779 00780 INTEGER :: i, stat 00781 LOGICAL :: failure, my_psi0_changed, 00782 my_s_struct_changed 00783 00784 failure=.FALSE. 00785 my_s_struct_changed=.FALSE. 00786 my_psi0_changed=.FALSE. 00787 00788 IF (PRESENT(s_struct_changed)) my_s_struct_changed=s_struct_changed 00789 IF (PRESENT(psi0_changed)) my_psi0_changed=psi0_changed 00790 00791 CPPrecondition(ASSOCIATED(kpp1_env),cp_failure_level,routineP,error,failure) 00792 IF (.NOT. failure) THEN 00793 CPPrecondition(kpp1_env%ref_count>0,cp_failure_level,routineP,error,failure) 00794 END IF 00795 IF (.not.failure) THEN 00796 IF (my_s_struct_changed) THEN 00797 IF (ASSOCIATED(kpp1_env%v_ao)) THEN 00798 CALL cp_dbcsr_deallocate_matrix_set(kpp1_env%v_ao,error=error) 00799 END IF 00800 END IF 00801 IF (my_s_struct_changed.or.my_psi0_changed) THEN 00802 IF (ASSOCIATED(kpp1_env%drho_r)) THEN 00803 DEALLOCATE(kpp1_env%drho_r, stat=stat) 00804 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00805 END IF 00806 IF (ASSOCIATED(kpp1_env%deriv_set)) THEN 00807 CALL xc_dset_release(kpp1_env%deriv_set, error=error) 00808 NULLIFY(kpp1_env%deriv_set) 00809 END IF 00810 IF (ASSOCIATED(kpp1_env%spin_pot)) THEN 00811 DEALLOCATE(kpp1_env%spin_pot, stat=stat) 00812 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00813 END IF 00814 IF (ASSOCIATED(kpp1_env%grad_pot)) THEN 00815 DEALLOCATE(kpp1_env%grad_pot, stat=stat) 00816 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00817 END IF 00818 IF (ASSOCIATED(kpp1_env%ndiag_term)) THEN 00819 DEALLOCATE(kpp1_env%ndiag_term, stat=stat) 00820 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00821 END IF 00822 CALL xc_rho_set_release(kpp1_env%rho_set,error=error) ! it would be better to pass a pw pool 00823 END IF 00824 IF (PRESENT(grid_changed)) THEN 00825 IF (grid_changed) THEN 00826 IF (ASSOCIATED(kpp1_env%v_rspace)) THEN 00827 DO i=1,SIZE(kpp1_env%v_rspace) 00828 CALL pw_release(kpp1_env%v_rspace(i)%pw,error=error) 00829 END DO 00830 DEALLOCATE(kpp1_env%v_rspace,stat=stat) 00831 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00832 END IF 00833 END IF 00834 END IF 00835 END IF 00836 END SUBROUTINE kpp1_did_change 00837 00838 ! ***************************************************************************** 00839 SUBROUTINE print_densities(kpp1_env, rho1, rho1_tot_gspace, output_unit, error) 00840 00841 TYPE(qs_kpp1_env_type), POINTER :: kpp1_env 00842 TYPE(qs_rho_type), POINTER :: rho1 00843 TYPE(pw_p_type), INTENT(IN) :: rho1_tot_gspace 00844 INTEGER :: output_unit 00845 TYPE(cp_error_type), INTENT(INOUT) :: error 00846 00847 CHARACTER(len=*), PARAMETER :: routineN = 'print_densities', 00848 routineP = moduleN//':'//routineN 00849 00850 REAL(KIND=dp) :: total_rho_gspace 00851 00852 total_rho_gspace = pw_integrate_function(rho1_tot_gspace%pw,isign=-1,error=error) 00853 IF (output_unit>0) THEN 00854 WRITE (UNIT=output_unit,FMT="(T3,A,T60,F20.10)")& 00855 "KPP1 total charge density (r-space):",& 00856 accurate_sum(rho1%tot_rho_r),& 00857 "KPP1 total charge density (g-space):",& 00858 total_rho_gspace 00859 END IF 00860 00861 END SUBROUTINE print_densities 00862 00863 END MODULE qs_kpp1_env_methods
1.7.3