|
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 ! ***************************************************************************** 00012 MODULE xc_rho_set_types 00013 USE cell_types, ONLY: cell_type 00014 USE cp_array_r_utils, ONLY: cp_3d_r_p_type 00015 USE f77_blas 00016 USE input_constants, ONLY: & 00017 xc_deriv_collocate, xc_deriv_nn10_smooth, xc_deriv_nn50_smooth, & 00018 xc_deriv_pw, xc_deriv_spline2, xc_deriv_spline2_smooth, & 00019 xc_deriv_spline3, xc_deriv_spline3_smooth, xc_rho_nn10, xc_rho_nn50, & 00020 xc_rho_no_smooth, xc_rho_spline2_smooth, xc_rho_spline3_smooth 00021 USE kinds, ONLY: dp 00022 USE pw_methods, ONLY: pw_copy,& 00023 pw_derive,& 00024 pw_transfer,& 00025 pw_zero 00026 USE pw_pool_types, ONLY: pw_pool_create_cr3d,& 00027 pw_pool_create_pw,& 00028 pw_pool_give_back_cr3d,& 00029 pw_pool_give_back_pw,& 00030 pw_pool_type 00031 USE pw_spline_utils, ONLY: & 00032 nn10_coeffs, nn10_deriv_coeffs, nn50_coeffs, nn50_deriv_coeffs, & 00033 pw_nn_deriv_r, pw_nn_smear_r, pw_spline2_deriv_g, & 00034 pw_spline2_interpolate_values_g, pw_spline3_deriv_g, & 00035 pw_spline3_interpolate_values_g, pw_spline_scale_deriv, & 00036 spline2_coeffs, spline2_deriv_coeffs, spline3_coeffs, & 00037 spline3_deriv_coeffs 00038 USE pw_types, ONLY: COMPLEXDATA1D,& 00039 REALDATA3D,& 00040 REALSPACE,& 00041 RECIPROCALSPACE,& 00042 pw_p_type,& 00043 pw_type 00044 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_equal,& 00045 xc_rho_cflags_setall,& 00046 xc_rho_cflags_type 00047 #include "cp_common_uses.h" 00048 00049 IMPLICIT NONE 00050 PRIVATE 00051 LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.FALSE. 00052 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_rho_set_types' 00053 INTEGER, SAVE :: last_rho_set_id=0 00054 00055 PUBLIC :: xc_rho_set_type 00056 PUBLIC :: xc_rho_set_create, xc_rho_set_retain, xc_rho_set_release,& 00057 xc_rho_set_update, xc_rho_set_get 00058 00059 ! ***************************************************************************** 00094 TYPE xc_rho_set_type 00095 INTEGER :: ref_count, id_nr 00096 INTEGER, DIMENSION(2,3) :: local_bounds 00097 REAL(kind=dp) :: rho_cutoff, drho_cutoff, tau_cutoff 00098 TYPE(xc_rho_cflags_type) :: owns, has 00099 ! for spin restricted systems 00100 REAL(KIND = dp), DIMENSION(:,:,:), POINTER :: rho 00101 TYPE(cp_3d_r_p_type), DIMENSION(3) :: drho 00102 REAL(KIND = dp), DIMENSION(:,:,:), POINTER :: norm_drho 00103 REAL(KIND = dp), DIMENSION(:,:,:), POINTER :: rho_1_3 00104 REAL(kind = dp), DIMENSION(:,:,:), POINTER :: tau 00105 ! for UNrestricted systems 00106 REAL(KIND = dp), DIMENSION(:,:,:), POINTER :: rhoa, rhob 00107 TYPE(cp_3d_r_p_type), DIMENSION(3) :: drhoa, drhob 00108 REAL(KIND = dp), DIMENSION(:,:,:), POINTER :: norm_drhoa, norm_drhob 00109 REAL(KIND = dp), DIMENSION(:,:,:), POINTER :: drhoa_drhob 00110 REAL(kind = dp), DIMENSION(:,:,:), POINTER :: rhoa_1_3, rhob_1_3 00111 REAL(kind = dp), DIMENSION(:,:,:), POINTER :: tau_a, tau_b 00112 REAL(kind = dp), DIMENSION(:,:,:), POINTER :: laplace_rho, laplace_rhoa, laplace_rhob 00113 END TYPE xc_rho_set_type 00114 00115 CONTAINS 00116 00117 ! ***************************************************************************** 00123 SUBROUTINE xc_rho_set_create(rho_set,local_bounds,rho_cutoff,drho_cutoff,& 00124 tau_cutoff,error) 00125 TYPE(xc_rho_set_type), POINTER :: rho_set 00126 INTEGER, DIMENSION(2, 3), INTENT(in) :: local_bounds 00127 REAL(kind=dp), INTENT(in), OPTIONAL :: rho_cutoff, drho_cutoff, 00128 tau_cutoff 00129 TYPE(cp_error_type), INTENT(inout) :: error 00130 00131 CHARACTER(len=*), PARAMETER :: routineN = 'xc_rho_set_create', 00132 routineP = moduleN//':'//routineN 00133 00134 INTEGER :: i, stat 00135 LOGICAL :: failure 00136 00137 failure=.FALSE. 00138 00139 CPPrecondition(.NOT.ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00140 ALLOCATE(rho_set, stat=stat) 00141 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00142 IF (.NOT. failure) THEN 00143 rho_set%ref_count=1 00144 last_rho_set_id=last_rho_set_id+1 00145 rho_set%id_nr=last_rho_set_id 00146 rho_set%rho_cutoff=EPSILON(0.0_dp) 00147 IF (PRESENT(rho_cutoff)) rho_set%rho_cutoff=rho_cutoff 00148 rho_set%drho_cutoff=EPSILON(0.0_dp) 00149 IF (PRESENT(drho_cutoff)) rho_set%drho_cutoff=drho_cutoff 00150 rho_set%tau_cutoff=EPSILON(0.0_dp) 00151 IF (PRESENT(tau_cutoff)) rho_set%tau_cutoff=tau_cutoff 00152 rho_set%local_bounds=local_bounds 00153 CALL xc_rho_cflags_setall(rho_set%owns,.TRUE.,error=error) 00154 CALL xc_rho_cflags_setall(rho_set%has,.FALSE.,error=error) 00155 NULLIFY(rho_set%rho) 00156 DO i=1,3 00157 NULLIFY(rho_set%drho(i)%array) 00158 END DO 00159 NULLIFY(rho_set%rho_1_3) 00160 NULLIFY(rho_set%norm_drho,rho_set%rhoa,rho_set%rhob) 00161 DO i=1,3 00162 NULLIFY(rho_set%drhoa(i)%array, rho_set%drhob(i)%array) 00163 END DO 00164 NULLIFY(rho_set%norm_drhoa, rho_set%norm_drhob, & 00165 rho_set%drhoa_drhob,rho_set%rhoa_1_3,rho_set%rhob_1_3,& 00166 rho_set%tau,rho_set%tau_a,rho_set%tau_b, rho_set%laplace_rho, rho_set%laplace_rhoa,& 00167 rho_set%laplace_rhob) 00168 END IF 00169 END SUBROUTINE xc_rho_set_create 00170 00171 ! ***************************************************************************** 00177 SUBROUTINE xc_rho_set_retain(rho_set, error) 00178 TYPE(xc_rho_set_type), POINTER :: rho_set 00179 TYPE(cp_error_type), INTENT(inout) :: error 00180 00181 CHARACTER(len=*), PARAMETER :: routineN = 'xc_rho_set_retain', 00182 routineP = moduleN//':'//routineN 00183 00184 LOGICAL :: failure 00185 00186 failure=.FALSE. 00187 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00188 IF (.NOT. failure) THEN 00189 CPPreconditionNoFail(rho_set%ref_count>0,cp_failure_level,routineP,error) 00190 rho_set%ref_count=rho_set%ref_count+1 00191 END IF 00192 END SUBROUTINE xc_rho_set_retain 00193 00194 ! ***************************************************************************** 00201 SUBROUTINE xc_rho_set_release(rho_set, pw_pool, error) 00202 TYPE(xc_rho_set_type), POINTER :: rho_set 00203 TYPE(pw_pool_type), OPTIONAL, POINTER :: pw_pool 00204 TYPE(cp_error_type), INTENT(inout) :: error 00205 00206 CHARACTER(len=*), PARAMETER :: routineN = 'xc_rho_set_release', 00207 routineP = moduleN//':'//routineN 00208 00209 INTEGER :: i, stat 00210 LOGICAL :: failure 00211 00212 failure=.FALSE. 00213 00214 IF (ASSOCIATED(rho_set)) THEN 00215 CPPreconditionNoFail(rho_set%ref_count>0,cp_failure_level,routineP,error) 00216 rho_set%ref_count=rho_set%ref_count-1 00217 IF (rho_set%ref_count==0) THEN 00218 IF (PRESENT(pw_pool)) THEN 00219 IF (ASSOCIATED(pw_pool)) THEN 00220 rho_set%ref_count=1 00221 CALL xc_rho_set_clean(rho_set,pw_pool,error) 00222 rho_set%ref_count=0 00223 ELSE 00224 CPPrecondition(.FALSE.,cp_warning_level,routineP,error,failure) 00225 END IF 00226 END IF 00227 00228 rho_set%local_bounds(1,:)=-HUGE(0) ! we want to crash... 00229 rho_set%local_bounds(1,:)=HUGE(0) 00230 IF (rho_set%owns%rho .AND. ASSOCIATED(rho_set%rho)) THEN 00231 DEALLOCATE(rho_set%rho, stat=stat) 00232 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00233 END IF 00234 IF (rho_set%owns%rho_spin) THEN 00235 IF (ASSOCIATED(rho_set%rhoa)) THEN 00236 DEALLOCATE(rho_set%rhoa, stat=stat) 00237 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00238 END IF 00239 IF (ASSOCIATED(rho_set%rhob)) THEN 00240 DEALLOCATE(rho_set%rhob, stat=stat) 00241 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00242 END IF 00243 END IF 00244 IF (rho_set%owns%rho_1_3.AND.ASSOCIATED(rho_set%rho_1_3)) THEN 00245 DEALLOCATE(rho_set%rho_1_3, stat=stat) 00246 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00247 END IF 00248 IF (rho_set%owns%rho_spin) THEN 00249 IF (ASSOCIATED(rho_set%rhoa_1_3)) THEN 00250 DEALLOCATE(rho_set%rhoa_1_3, stat=stat) 00251 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00252 END IF 00253 IF (ASSOCIATED(rho_set%rhob_1_3)) THEN 00254 DEALLOCATE(rho_set%rhob_1_3, stat=stat) 00255 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00256 END IF 00257 END IF 00258 IF (rho_set%owns%drho) THEN 00259 DO i=1,3 00260 IF (ASSOCIATED(rho_set%drho(i)%array)) THEN 00261 DEALLOCATE(rho_set%drho(i)%array, stat=stat) 00262 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00263 END IF 00264 END DO 00265 END IF 00266 IF (rho_set%owns%drho_spin) THEN 00267 DO i=1,3 00268 IF (ASSOCIATED(rho_set%drhoa(i)%array)) THEN 00269 DEALLOCATE(rho_set%drhoa(i)%array, stat=stat) 00270 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00271 END IF 00272 IF (ASSOCIATED(rho_set%drhob(i)%array)) THEN 00273 DEALLOCATE(rho_set%drhob(i)%array, stat=stat) 00274 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00275 END IF 00276 END DO 00277 END IF 00278 IF (rho_set%owns%laplace_rho.AND.ASSOCIATED(rho_set%laplace_rho)) THEN 00279 DEALLOCATE(rho_set%laplace_rho, stat=stat) 00280 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00281 END IF 00282 00283 IF (rho_set%owns%norm_drho.AND.ASSOCIATED(rho_set%norm_drho)) THEN 00284 DEALLOCATE(rho_set%norm_drho, stat=stat) 00285 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00286 END IF 00287 IF (rho_set%owns%laplace_rho_spin) THEN 00288 IF (ASSOCIATED(rho_set%laplace_rhoa)) THEN 00289 DEALLOCATE(rho_set%laplace_rhoa, stat=stat) 00290 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00291 END IF 00292 IF (ASSOCIATED(rho_set%laplace_rhob)) THEN 00293 DEALLOCATE(rho_set%laplace_rhob, stat=stat) 00294 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00295 END IF 00296 END IF 00297 00298 IF (rho_set%owns%norm_drho_spin) THEN 00299 IF (ASSOCIATED(rho_set%norm_drhoa)) THEN 00300 DEALLOCATE(rho_set%norm_drhoa, stat=stat) 00301 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00302 END IF 00303 IF (ASSOCIATED(rho_set%norm_drhob)) THEN 00304 DEALLOCATE(rho_set%norm_drhob, stat=stat) 00305 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00306 END IF 00307 END IF 00308 IF (rho_set%owns%drhoa_drhob.AND.ASSOCIATED(rho_set%drhoa_drhob)) THEN 00309 DEALLOCATE(rho_set%drhoa_drhob, stat=stat) 00310 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00311 END IF 00312 IF (rho_set%owns%tau.AND.ASSOCIATED(rho_set%tau)) THEN 00313 DEALLOCATE(rho_set%tau, stat=stat) 00314 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00315 END IF 00316 IF (rho_set%owns%tau_spin) THEN 00317 IF (ASSOCIATED(rho_set%tau_a)) THEN 00318 DEALLOCATE(rho_set%tau_a, stat=stat) 00319 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00320 END IF 00321 IF (ASSOCIATED(rho_set%tau_b)) THEN 00322 DEALLOCATE(rho_set%tau_b, stat=stat) 00323 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00324 END IF 00325 END IF 00326 DEALLOCATE(rho_set, stat=stat) 00327 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00328 END IF 00329 END IF 00330 NULLIFY(rho_set) 00331 END SUBROUTINE xc_rho_set_release 00332 00333 ! ***************************************************************************** 00344 SUBROUTINE xc_rho_set_get(rho_set, can_return_null, rho, drho, norm_drho,& 00345 rhoa, rhob, norm_drhoa, norm_drhob, drhoa_drhob,rho_1_3,rhoa_1_3,& 00346 rhob_1_3,laplace_rho,laplace_rhoa,laplace_rhob,drhoa,drhob,rho_cutoff,& 00347 drho_cutoff,tau_cutoff,tau,tau_a,tau_b,local_bounds,error) 00348 TYPE(xc_rho_set_type), POINTER :: rho_set 00349 LOGICAL, INTENT(in), OPTIONAL :: can_return_null 00350 REAL(KIND=dp), DIMENSION(:, :, :), 00351 OPTIONAL, POINTER :: rho 00352 TYPE(cp_3d_r_p_type), DIMENSION(:), 00353 OPTIONAL, POINTER :: drho 00354 REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, POINTER :: norm_drho, rhoa, 00355 rhob, norm_drhoa, norm_drhob, drhoa_drhob, rho_1_3, rhoa_1_3, rhob_1_3, 00356 laplace_rho, laplace_rhoa, laplace_rhob 00357 TYPE(cp_3d_r_p_type), DIMENSION(:), 00358 OPTIONAL, POINTER :: drhoa, drhob 00359 REAL(kind=dp), INTENT(out), OPTIONAL :: rho_cutoff, drho_cutoff, 00360 tau_cutoff 00361 REAL(KIND=dp), DIMENSION(:, :, :), 00362 OPTIONAL, POINTER :: tau, tau_a, tau_b 00363 INTEGER, DIMENSION(:, :), OPTIONAL, 00364 POINTER :: local_bounds 00365 TYPE(cp_error_type), INTENT(inout) :: error 00366 00367 CHARACTER(len=*), PARAMETER :: routineN = 'xc_rho_set_get', 00368 routineP = moduleN//':'//routineN 00369 00370 INTEGER :: i 00371 LOGICAL :: failure, my_can_return_null 00372 00373 failure=.FALSE. 00374 my_can_return_null=.FALSE. 00375 IF (PRESENT(can_return_null)) my_can_return_null=can_return_null 00376 00377 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00378 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00379 IF (.NOT. failure) THEN 00380 IF (PRESENT(rho)) THEN 00381 rho => rho_set%rho 00382 CPPrecondition(my_can_return_null.OR.ASSOCIATED(rho),cp_failure_level,routineP,error,failure) 00383 END IF 00384 IF (PRESENT(drho)) THEN 00385 drho => rho_set%drho 00386 IF (.NOT.my_can_return_null) THEN 00387 DO i=1,3 00388 CPPrecondition(ASSOCIATED(rho_set%drho(i)%array),cp_failure_level,routineP,error,failure) 00389 END DO 00390 END IF 00391 END IF 00392 IF (PRESENT(norm_drho)) THEN 00393 norm_drho => rho_set%norm_drho 00394 CPPrecondition(my_can_return_null.OR.ASSOCIATED(norm_drho),cp_failure_level,routineP,error,failure) 00395 END IF 00396 IF (PRESENT(laplace_rho)) THEN 00397 laplace_rho => rho_set%laplace_rho 00398 CPPrecondition(my_can_return_null.OR.ASSOCIATED(laplace_rho),cp_failure_level,routineP,error,failure) 00399 END IF 00400 IF (PRESENT(rhoa)) THEN 00401 rhoa => rho_set%rhoa 00402 CPPrecondition(my_can_return_null.OR.ASSOCIATED(rhoa),cp_failure_level,routineP,error,failure) 00403 END IF 00404 IF (PRESENT(rhob)) THEN 00405 rhob => rho_set%rhob 00406 CPPrecondition(my_can_return_null.OR.ASSOCIATED(rhob),cp_failure_level,routineP,error,failure) 00407 END IF 00408 IF (PRESENT(drhoa)) THEN 00409 drhoa => rho_set%drhoa 00410 IF (.NOT.my_can_return_null) THEN 00411 DO i=1,3 00412 CPAssert(ASSOCIATED(rho_set%drhoa(i)%array),cp_failure_level,routineP,error,failure) 00413 END DO 00414 END IF 00415 END IF 00416 IF (PRESENT(drhob)) THEN 00417 drhob => rho_set%drhob 00418 IF (.NOT.my_can_return_null) THEN 00419 DO i=1,3 00420 CPPrecondition(ASSOCIATED(rho_set%drhob(i)%array),cp_failure_level,routineP,error,failure) 00421 END DO 00422 END IF 00423 END IF 00424 IF (PRESENT(laplace_rhoa)) THEN 00425 laplace_rhoa => rho_set%laplace_rhoa 00426 CPPrecondition(my_can_return_null.OR.ASSOCIATED(laplace_rhoa),cp_failure_level,routineP,error,failure) 00427 END IF 00428 IF (PRESENT(laplace_rhob)) THEN 00429 laplace_rhob => rho_set%laplace_rhob 00430 CPPrecondition(my_can_return_null.OR.ASSOCIATED(laplace_rhob),cp_failure_level,routineP,error,failure) 00431 END IF 00432 IF (PRESENT(norm_drhoa)) THEN 00433 norm_drhoa => rho_set%norm_drhoa 00434 CPPrecondition(my_can_return_null.OR.ASSOCIATED(norm_drhoa),cp_failure_level,routineP,error,failure) 00435 END IF 00436 IF (PRESENT(norm_drhob)) THEN 00437 norm_drhob => rho_set%norm_drhob 00438 CPPrecondition(my_can_return_null.OR.ASSOCIATED(norm_drhob),cp_failure_level,routineP,error,failure) 00439 END IF 00440 IF (PRESENT(drhoa_drhob)) THEN 00441 drhoa_drhob => rho_set%drhoa_drhob 00442 CPPrecondition(my_can_return_null.OR.ASSOCIATED(drhoa_drhob),cp_failure_level,routineP,error,failure) 00443 END IF 00444 IF (PRESENT(rho_1_3)) THEN 00445 rho_1_3 => rho_set%rho_1_3 00446 CPPrecondition(my_can_return_null.OR.ASSOCIATED(rho_1_3),cp_failure_level,routineP,error,failure) 00447 END IF 00448 IF (PRESENT(rhoa_1_3)) THEN 00449 rhoa_1_3 => rho_set%rhoa_1_3 00450 CPPrecondition(my_can_return_null.OR.ASSOCIATED(rhoa_1_3),cp_failure_level,routineP,error,failure) 00451 END IF 00452 IF (PRESENT(rhob_1_3)) THEN 00453 rhob_1_3 => rho_set%rhob_1_3 00454 CPPrecondition(my_can_return_null.OR.ASSOCIATED(rhob_1_3),cp_failure_level,routineP,error,failure) 00455 END IF 00456 IF (PRESENT(tau)) THEN 00457 tau => rho_set%tau 00458 CPPrecondition(my_can_return_null.OR.ASSOCIATED(tau),cp_failure_level,routineP,error,failure) 00459 END IF 00460 IF (PRESENT(tau_a)) THEN 00461 tau_a => rho_set%tau_a 00462 CPPrecondition(my_can_return_null.OR.ASSOCIATED(tau_a),cp_failure_level,routineP,error,failure) 00463 END IF 00464 IF (PRESENT(tau_b)) THEN 00465 tau_b => rho_set%tau_b 00466 CPPrecondition(my_can_return_null.OR.ASSOCIATED(tau_b),cp_failure_level,routineP,error,failure) 00467 END IF 00468 IF (PRESENT(rho_cutoff)) rho_cutoff=rho_set%rho_cutoff 00469 IF (PRESENT(drho_cutoff)) drho_cutoff=rho_set%drho_cutoff 00470 IF (PRESENT(tau_cutoff)) tau_cutoff=rho_set%tau_cutoff 00471 IF (PRESENT(local_bounds)) local_bounds => rho_set%local_bounds 00472 END IF 00473 END SUBROUTINE xc_rho_set_get 00474 00475 ! ***************************************************************************** 00484 SUBROUTINE xc_rho_set_clean(rho_set,pw_pool,error) 00485 TYPE(xc_rho_set_type), POINTER :: rho_set 00486 TYPE(pw_pool_type), POINTER :: pw_pool 00487 TYPE(cp_error_type), INTENT(inout) :: error 00488 00489 CHARACTER(len=*), PARAMETER :: routineN = 'xc_rho_set_clean', 00490 routineP = moduleN//':'//routineN 00491 00492 INTEGER :: idir 00493 LOGICAL :: failure 00494 00495 failure=.FALSE. 00496 00497 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00498 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00499 00500 IF (.NOT. failure) THEN 00501 00502 IF (rho_set%owns%rho) THEN 00503 CALL pw_pool_give_back_cr3d(pw_pool,rho_set%rho, & 00504 accept_non_compatible=.TRUE., error=error) 00505 ELSE 00506 NULLIFY(rho_set%rho) 00507 END IF 00508 IF (rho_set%owns%rho_1_3) THEN 00509 CALL pw_pool_give_back_cr3d(pw_pool,rho_set%rho_1_3,& 00510 accept_non_compatible=.TRUE.,error=error) 00511 ELSE 00512 NULLIFY(rho_set%rho_1_3) 00513 END IF 00514 IF (rho_set%owns%drho) THEN 00515 DO idir=1,3 00516 CALL pw_pool_give_back_cr3d(pw_pool,rho_set%drho(idir)%array,& 00517 accept_non_compatible=.TRUE.,error=error) 00518 END DO 00519 ELSE 00520 DO idir=1,3 00521 NULLIFY(rho_set%drho(idir)%array) 00522 END DO 00523 END IF 00524 IF (rho_set%owns%norm_drho) THEN 00525 CALL pw_pool_give_back_cr3d(pw_pool,rho_set%norm_drho,& 00526 accept_non_compatible=.TRUE.,error=error) 00527 ELSE 00528 NULLIFY(rho_set%norm_drho) 00529 END IF 00530 IF (rho_set%owns%laplace_rho) THEN 00531 CALL pw_pool_give_back_cr3d(pw_pool,rho_set%laplace_rho,& 00532 accept_non_compatible=.TRUE.,error=error) 00533 ELSE 00534 NULLIFY(rho_set%laplace_rho) 00535 END IF 00536 IF (rho_set%owns%tau) THEN 00537 CALL pw_pool_give_back_cr3d(pw_pool, rho_set%tau,& 00538 accept_non_compatible=.TRUE., error=error) 00539 ELSE 00540 NULLIFY(rho_set%tau) 00541 END IF 00542 IF (rho_set%owns%rho_spin) THEN 00543 CALL pw_pool_give_back_cr3d(pw_pool,rho_set%rhoa,& 00544 accept_non_compatible=.TRUE., error=error) 00545 CALL pw_pool_give_back_cr3d(pw_pool,rho_set%rhob,& 00546 accept_non_compatible=.TRUE., error=error) 00547 ELSE 00548 NULLIFY(rho_set%rhoa,rho_set%rhob) 00549 END IF 00550 IF (rho_set%owns%rho_spin_1_3) THEN 00551 CALL pw_pool_give_back_cr3d(pw_pool,rho_set%rhoa_1_3,& 00552 accept_non_compatible=.TRUE.,error=error) 00553 CALL pw_pool_give_back_cr3d(pw_pool,rho_set%rhob_1_3,& 00554 accept_non_compatible=.TRUE.,error=error) 00555 ELSE 00556 NULLIFY(rho_set%rhoa_1_3,rho_set%rhob_1_3) 00557 END IF 00558 IF (rho_set%owns%drho_spin) THEN 00559 DO idir=1,3 00560 CALL pw_pool_give_back_cr3d(pw_pool,rho_set%drhoa(idir)%array,& 00561 accept_non_compatible=.TRUE.,error=error) 00562 CALL pw_pool_give_back_cr3d(pw_pool,rho_set%drhob(idir)%array,& 00563 accept_non_compatible=.TRUE.,error=error) 00564 END DO 00565 ELSE 00566 DO idir=1,3 00567 NULLIFY(rho_set%drhoa(idir)%array,rho_set%drhob(idir)%array) 00568 END DO 00569 END IF 00570 IF (rho_set%owns%laplace_rho_spin) THEN 00571 CALL pw_pool_give_back_cr3d(pw_pool,rho_set%laplace_rhoa,& 00572 accept_non_compatible=.TRUE.,error=error) 00573 CALL pw_pool_give_back_cr3d(pw_pool,rho_set%laplace_rhob,& 00574 accept_non_compatible=.TRUE.,error=error) 00575 ELSE 00576 NULLIFY(rho_set%laplace_rhoa, rho_set%laplace_rhob) 00577 END IF 00578 IF (rho_set%owns%norm_drho_spin) THEN 00579 CALL pw_pool_give_back_cr3d(pw_pool,rho_set%norm_drhoa,& 00580 accept_non_compatible=.TRUE.,error=error) 00581 CALL pw_pool_give_back_cr3d(pw_pool,rho_set%norm_drhob,& 00582 accept_non_compatible=.TRUE.,error=error) 00583 ELSE 00584 NULLIFY(rho_set%norm_drhoa, rho_set%norm_drhob) 00585 END IF 00586 IF (rho_set%owns%drhoa_drhob) THEN 00587 CALL pw_pool_give_back_cr3d(pw_pool, rho_set%drhoa_drhob,& 00588 accept_non_compatible=.TRUE., error=error) 00589 ELSE 00590 NULLIFY(rho_set%drhoa_drhob) 00591 END IF 00592 IF (rho_set%owns%tau_spin) THEN 00593 CALL pw_pool_give_back_cr3d(pw_pool, rho_set%tau_a,& 00594 accept_non_compatible=.TRUE., error=error) 00595 CALL pw_pool_give_back_cr3d(pw_pool, rho_set%tau_b,& 00596 accept_non_compatible=.TRUE., error=error) 00597 ELSE 00598 NULLIFY(rho_set%tau_a,rho_set%tau_b) 00599 END IF 00600 00601 CALL xc_rho_cflags_setall(rho_set%has,.FALSE.,error=error) 00602 CALL xc_rho_cflags_setall(rho_set%owns,.FALSE.,error=error) 00603 00604 END IF 00605 END SUBROUTINE xc_rho_set_clean 00606 00607 ! ***************************************************************************** 00621 SUBROUTINE xc_rho_set_update(rho_set, rho_r, rho_g, tau, needs,& 00622 xc_deriv_method_id,xc_rho_smooth_id,cell,pw_pool,error) 00623 TYPE(xc_rho_set_type), POINTER :: rho_set 00624 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r, rho_g, tau 00625 TYPE(xc_rho_cflags_type), INTENT(in) :: needs 00626 INTEGER, INTENT(IN) :: xc_deriv_method_id, 00627 xc_rho_smooth_id 00628 TYPE(cell_type), POINTER :: cell 00629 TYPE(pw_pool_type), POINTER :: pw_pool 00630 TYPE(cp_error_type), INTENT(inout) :: error 00631 00632 CHARACTER(len=*), PARAMETER :: routineN = 'xc_rho_set_update', 00633 routineP = moduleN//':'//routineN 00634 REAL(KIND=dp), PARAMETER :: f13 = (1.0_dp/3.0_dp) 00635 00636 INTEGER :: i, idir, ispin, j, k, nspins 00637 INTEGER, DIMENSION(3, 3) :: nd, nd_laplace 00638 LOGICAL :: failure, gradient_f, 00639 my_rho_g_local, 00640 my_rho_r_local, needs_rho_g 00641 REAL(kind=dp) :: rho_cutoff 00642 TYPE(pw_p_type), DIMENSION(2) :: my_rho_r 00643 TYPE(pw_p_type), DIMENSION(3) :: drho_r_att 00644 TYPE(pw_p_type), DIMENSION(3, 2) :: drho_r, laplace_rho_r 00645 TYPE(pw_type), POINTER :: my_rho_g, tmp_g 00646 00647 failure=.FALSE. 00648 DO ispin=1,2 00649 NULLIFY(my_rho_r(ispin)%pw) 00650 DO idir=1,3 00651 NULLIFY(drho_r(idir,ispin)%pw) 00652 END DO 00653 END DO 00654 DO idir=1,3 00655 NULLIFY(drho_r_att(idir)%pw) 00656 END DO 00657 NULLIFY(tmp_g,my_rho_g) 00658 nd = RESHAPE ((/1,0,0,0,1,0,0,0,1/),(/3,3/)) 00659 nd_laplace = RESHAPE((/2,0,0,0,2,0,0,0,2/),(/3,3/)) 00660 00661 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00662 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00663 CALL cp_assert(ALL(rho_set%local_bounds==pw_pool%pw_grid%bounds_local),& 00664 cp_failure_level,cp_assertion_failed,routineP,& 00665 "pw_pool cr3d have different size than expected",error,failure) 00666 nspins = SIZE(rho_r) 00667 rho_set%local_bounds=rho_r(1)%pw%pw_grid%bounds_local 00668 rho_cutoff=0.5*rho_set%rho_cutoff 00669 00670 my_rho_g_local=.FALSE. 00671 ! some checks 00672 SELECT CASE(nspins) 00673 CASE(1) 00674 CPPrecondition(SIZE(rho_r)==1,cp_failure_level,routineP,error,failure) 00675 CPPrecondition(ASSOCIATED(rho_r(1)%pw),cp_failure_level,routineP,error,failure) 00676 CPPrecondition(rho_r(1)%pw%in_use==REALDATA3D,cp_failure_level,routineP,error,failure) 00677 CPPrecondition(.NOT.needs%rho_spin,cp_failure_level,routineP,error,failure) 00678 CPPrecondition(.NOT.needs%drho_spin,cp_failure_level,routineP,error,failure) 00679 CPPrecondition(.NOT.needs%norm_drho_spin,cp_failure_level,routineP,error,failure) 00680 CPPrecondition(.NOT.needs%drhoa_drhob,cp_failure_level,routineP,error,failure) 00681 CPPrecondition(.NOT.needs%rho_spin_1_3,cp_failure_level,routineP,error,failure) 00682 CASE(2) 00683 CPPrecondition(SIZE(rho_r)==2,cp_failure_level,routineP,error,failure) 00684 CPPrecondition(ASSOCIATED(rho_r(1)%pw),cp_failure_level,routineP,error,failure) 00685 CPPrecondition(ASSOCIATED(rho_r(2)%pw),cp_failure_level,routineP,error,failure) 00686 CPPrecondition(rho_r(1)%pw%in_use==REALDATA3D,cp_failure_level,routineP,error,failure) 00687 CPPrecondition(rho_r(2)%pw%in_use==REALDATA3D,cp_failure_level,routineP,error,failure) 00688 CASE default 00689 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00690 END SELECT 00691 00692 CALL xc_rho_set_clean(rho_set,pw_pool=pw_pool,error=error) 00693 00694 gradient_f=(needs%drho_spin.OR.needs%norm_drho_spin.OR.& 00695 needs%drhoa_drhob.OR.needs%drho.OR.needs%norm_drho.OR.needs%laplace_rho & 00696 .OR.needs%laplace_rho_spin) 00697 needs_rho_g=(xc_deriv_method_id==xc_deriv_spline3.OR.& 00698 xc_deriv_method_id==xc_deriv_spline2.OR.& 00699 xc_deriv_method_id==xc_deriv_pw) 00700 DO ispin=1,nspins 00701 ! introduce a smoothing kernel on the density 00702 IF (xc_rho_smooth_id==xc_rho_no_smooth) THEN 00703 my_rho_r_local=.FALSE. 00704 my_rho_r(ispin)%pw => rho_r(ispin)%pw 00705 IF (needs_rho_g) THEN 00706 IF (ASSOCIATED(rho_g)) THEN 00707 my_rho_g_local=.FALSE. 00708 my_rho_g => rho_g(ispin)%pw 00709 END IF 00710 END IF 00711 00712 my_rho_r_local=.TRUE. 00713 CALL pw_pool_create_pw(pw_pool, my_rho_r(ispin)%pw,& 00714 use_data=REALDATA3D, in_space=REALSPACE, & 00715 error=error) 00716 CALL pw_copy(rho_r(ispin)%pw,my_rho_r(ispin)%pw,error=error) 00717 ELSE 00718 my_rho_r_local=.TRUE. 00719 CALL pw_pool_create_pw(pw_pool, my_rho_r(ispin)%pw,& 00720 use_data=REALDATA3D, in_space=REALSPACE, & 00721 error=error) 00722 00723 SELECT CASE(xc_rho_smooth_id) 00724 CASE (xc_rho_no_smooth) 00725 CALL pw_copy(rho_r(ispin)%pw,my_rho_r(ispin)%pw,error=error) 00726 CASE (xc_rho_spline2_smooth) 00727 CALL pw_zero(my_rho_r(ispin)%pw,error=error) 00728 CALL pw_nn_smear_r(pw_in=rho_r(ispin)%pw,& 00729 pw_out=my_rho_r(ispin)%pw,& 00730 coeffs=spline2_coeffs, error=error) 00731 CASE (xc_rho_spline3_smooth) 00732 CALL pw_zero(my_rho_r(ispin)%pw,error=error) 00733 CALL pw_nn_smear_r(pw_in=rho_r(ispin)%pw,& 00734 pw_out=my_rho_r(ispin)%pw,& 00735 coeffs=spline3_coeffs, error=error) 00736 CASE (xc_rho_nn10) 00737 CALL pw_zero(my_rho_r(ispin)%pw,error=error) 00738 CALL pw_nn_smear_r(pw_in=rho_r(ispin)%pw,& 00739 pw_out=my_rho_r(ispin)%pw,& 00740 coeffs=nn10_coeffs, error=error) 00741 CASE (xc_rho_nn50) 00742 CALL pw_zero(my_rho_r(ispin)%pw,error=error) 00743 CALL pw_nn_smear_r(pw_in=rho_r(ispin)%pw,& 00744 pw_out=my_rho_r(ispin)%pw,& 00745 coeffs=nn50_coeffs, error=error) 00746 CASE default 00747 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00748 END SELECT 00749 END IF 00750 00751 IF (gradient_f) THEN ! calculate the grad of rho 00752 ! normally when you need the gradient you need the whole gradient 00753 ! (for the partial integration) 00754 ! deriv rho 00755 DO idir=1,3 00756 NULLIFY(drho_r(idir,ispin)%pw) 00757 CALL pw_pool_create_pw(pw_pool,drho_r(idir,ispin)%pw, & 00758 use_data=REALDATA3D, in_space=REALSPACE, & 00759 error=error) 00760 END DO 00761 IF (needs_rho_g) THEN 00762 IF (.NOT.ASSOCIATED(my_rho_g)) THEN 00763 my_rho_g_local=.TRUE. 00764 CALL pw_pool_create_pw(pw_pool, my_rho_g,& 00765 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE, & 00766 error=error) 00767 CALL pw_transfer(my_rho_r(ispin)%pw,my_rho_g,error=error) 00768 END IF 00769 CALL pw_pool_create_pw(pw_pool, tmp_g,& 00770 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE, & 00771 error=error) 00772 SELECT CASE(xc_deriv_method_id) 00773 CASE (xc_deriv_pw) 00774 DO idir=1,3 00775 CALL pw_copy ( my_rho_g, tmp_g ,error=error) 00776 CALL pw_derive ( tmp_g, nd(:,idir) ,error=error) 00777 CALL pw_transfer ( tmp_g, drho_r(idir,ispin)%pw ,error=error) 00778 END DO 00779 IF(needs%laplace_rho.OR.needs%laplace_rho_spin) THEN 00780 DO idir=1,3 00781 NULLIFY(laplace_rho_r(idir,ispin)%pw) 00782 CALL pw_pool_create_pw(pw_pool,laplace_rho_r(idir,ispin)%pw, & 00783 use_data=REALDATA3D, in_space=REALSPACE, & 00784 error=error) 00785 CALL pw_copy ( my_rho_g, tmp_g ,error=error) 00786 CALL pw_derive ( tmp_g, nd_laplace(:,idir) ,error=error) 00787 CALL pw_transfer ( tmp_g, laplace_rho_r(idir,ispin)%pw ,error=error) 00788 END DO 00789 END IF 00790 CASE (xc_deriv_spline2) 00791 IF (.NOT.my_rho_g_local) THEN 00792 CALL pw_pool_create_pw(pw_pool, my_rho_g,& 00793 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE, & 00794 error=error) 00795 my_rho_g_local=.TRUE. 00796 CALL pw_copy(rho_g(ispin)%pw, my_rho_g,error=error) 00797 END IF 00798 CALL pw_spline2_interpolate_values_g(my_rho_g,error=error) 00799 DO idir=1,3 00800 CALL pw_copy ( my_rho_g, tmp_g ,error=error) 00801 CALL pw_spline2_deriv_g ( tmp_g, idir=idir, error=error ) 00802 CALL pw_transfer ( tmp_g, drho_r(idir,ispin)%pw ,error=error) 00803 END DO 00804 CASE (xc_deriv_spline3) 00805 IF (.NOT.my_rho_g_local) THEN 00806 CALL pw_pool_create_pw(pw_pool, my_rho_g,& 00807 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE, & 00808 error=error) 00809 CALL pw_copy(rho_g(ispin)%pw, my_rho_g,error=error) 00810 my_rho_g_local=.TRUE. 00811 END IF 00812 CALL pw_spline3_interpolate_values_g(my_rho_g,error=error) 00813 DO idir=1,3 00814 CALL pw_copy ( my_rho_g, tmp_g ,error=error) 00815 CALL pw_spline3_deriv_g ( tmp_g, idir=idir, error=error ) 00816 CALL pw_transfer ( tmp_g, drho_r(idir,ispin)%pw ,error=error) 00817 END DO 00818 CASE (xc_deriv_collocate) 00819 DO idir=1,3 00820 CALL pw_copy ( my_rho_g, tmp_g ,error=error) 00821 CALL pw_derive ( tmp_g, nd(:,idir) ,error=error) 00822 CALL pw_transfer ( tmp_g, drho_r(idir,ispin)%pw ,error=error) 00823 END DO 00824 CALL cp_unimplemented_error(fromWhere=routineP, & 00825 message="Drho collocation not implemented", & 00826 error=error, error_level=cp_failure_level) 00827 CASE default 00828 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00829 END SELECT 00830 CALL pw_pool_give_back_pw(pw_pool, tmp_g ,error=error) 00831 IF (my_rho_g_local) THEN 00832 my_rho_g_local=.FALSE. 00833 CALL pw_pool_give_back_pw(pw_pool, my_rho_g ,error=error) 00834 END IF 00835 ELSE 00836 SELECT CASE(xc_deriv_method_id) 00837 CASE (xc_deriv_spline2_smooth) 00838 DO idir=1,3 00839 CALL pw_zero(drho_r(idir,ispin)%pw,error=error) 00840 CALL pw_nn_deriv_r(pw_in=my_rho_r(ispin)%pw,& 00841 pw_out=drho_r(idir,ispin)%pw,& 00842 coeffs=spline2_deriv_coeffs, idir=idir, error=error) 00843 END DO 00844 CASE (xc_deriv_spline3_smooth) 00845 DO idir=1,3 00846 CALL pw_zero(drho_r(idir,ispin)%pw,error=error) 00847 CALL pw_nn_deriv_r(pw_in=my_rho_r(ispin)%pw,& 00848 pw_out=drho_r(idir,ispin)%pw,& 00849 coeffs=spline3_deriv_coeffs, idir=idir, error=error) 00850 END DO 00851 CASE (xc_deriv_nn10_smooth) 00852 DO idir=1,3 00853 CALL pw_zero(drho_r(idir,ispin)%pw,error=error) 00854 CALL pw_nn_deriv_r(pw_in=my_rho_r(ispin)%pw,& 00855 pw_out=drho_r(idir,ispin)%pw,& 00856 coeffs=nn10_deriv_coeffs, idir=idir, error=error) 00857 END DO 00858 CASE (xc_deriv_nn50_smooth) 00859 DO idir=1,3 00860 CALL pw_zero(drho_r(idir,ispin)%pw,error=error) 00861 CALL pw_nn_deriv_r(pw_in=my_rho_r(ispin)%pw,& 00862 pw_out=drho_r(idir,ispin)%pw,& 00863 coeffs=nn50_deriv_coeffs, idir=idir, error=error) 00864 END DO 00865 CASE (xc_deriv_collocate) 00866 DO idir=1,3 00867 CALL pw_copy ( my_rho_g, tmp_g ,error=error) 00868 CALL pw_derive ( tmp_g, nd(:,idir) ,error=error) 00869 CALL pw_transfer ( tmp_g, drho_r(idir,ispin)%pw ,error=error) 00870 END DO 00871 CALL cp_unimplemented_error(fromWhere=routineP, & 00872 message="Drho collocation not implemented", & 00873 error=error, error_level=cp_failure_level) 00874 CASE default 00875 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00876 END SELECT 00877 END IF 00878 00879 IF (xc_deriv_method_id/=xc_deriv_pw) THEN 00880 DO idir=1,3 00881 drho_r_att(idir)%pw => drho_r(idir,ispin)%pw 00882 END DO 00883 CALL pw_spline_scale_deriv(drho_r_att, cell=cell,& 00884 error=error) 00885 END IF 00886 00887 END IF 00888 00889 END DO 00890 00891 SELECT CASE(nspins) 00892 CASE(1) 00893 IF (needs%rho_1_3) THEN 00894 CALL pw_pool_create_cr3d(pw_pool, rho_set%rho_1_3, & 00895 error=error) 00896 !$omp parallel do default(none) private(i,j,k) shared(rho_set,my_rho_r) 00897 DO k=rho_set%local_bounds(1,3),rho_set%local_bounds(2,3) 00898 DO j=rho_set%local_bounds(1,2),rho_set%local_bounds(2,2) 00899 DO i=rho_set%local_bounds(1,1),rho_set%local_bounds(2,1) 00900 rho_set%rho_1_3(i,j,k)=MAX(my_rho_r(1)%pw%cr3d(i,j,k),0.0_dp)**f13 00901 END DO 00902 END DO 00903 END DO 00904 rho_set%owns%rho_1_3=.TRUE. 00905 rho_set%has%rho_1_3=.TRUE. 00906 END IF 00907 IF (needs%rho) THEN 00908 rho_set%rho => my_rho_r(1)%pw%cr3d 00909 IF (my_rho_r_local) NULLIFY(my_rho_r(1)%pw%cr3d) 00910 rho_set%owns%rho=my_rho_r_local 00911 rho_set%has%rho=.TRUE. 00912 END IF 00913 IF (needs%norm_drho) THEN 00914 CALL pw_pool_create_cr3d(pw_pool, rho_set%norm_drho, & 00915 error=error) 00916 !$omp parallel do default(none) private(i,j,k) shared(rho_set,drho_r) 00917 DO k=rho_set%local_bounds(1,3),rho_set%local_bounds(2,3) 00918 DO j=rho_set%local_bounds(1,2),rho_set%local_bounds(2,2) 00919 DO i=rho_set%local_bounds(1,1),rho_set%local_bounds(2,1) 00920 rho_set%norm_drho(i,j,k)=SQRT(& 00921 drho_r(1,1)%pw%cr3d(i,j,k)**2+& 00922 drho_r(2,1)%pw%cr3d(i,j,k)**2+& 00923 drho_r(3,1)%pw%cr3d(i,j,k)**2) 00924 END DO 00925 END DO 00926 END DO 00927 rho_set%owns%norm_drho=.TRUE. 00928 rho_set%has%norm_drho=.TRUE. 00929 END IF 00930 IF (needs%laplace_rho) THEN 00931 CALL pw_pool_create_cr3d(pw_pool, rho_set%laplace_rho, & 00932 error=error) 00933 !$omp parallel do default(none) private(i,j,k) shared(rho_set,laplace_rho_r) 00934 DO k=rho_set%local_bounds(1,3),rho_set%local_bounds(2,3) 00935 DO j=rho_set%local_bounds(1,2),rho_set%local_bounds(2,2) 00936 DO i=rho_set%local_bounds(1,1),rho_set%local_bounds(2,1) 00937 rho_set%laplace_rho(i,j,k)=& 00938 laplace_rho_r(1,1)%pw%cr3d(i,j,k)+& 00939 laplace_rho_r(2,1)%pw%cr3d(i,j,k)+& 00940 laplace_rho_r(3,1)%pw%cr3d(i,j,k) 00941 END DO 00942 END DO 00943 END DO 00944 rho_set%owns%laplace_rho=.TRUE. 00945 rho_set%has%laplace_rho=.TRUE. 00946 END IF 00947 00948 IF (needs%drho) THEN 00949 DO idir=1,3 00950 rho_set%drho(idir)%array => drho_r(idir,1)%pw%cr3d 00951 NULLIFY(drho_r(idir,1)%pw%cr3d) 00952 END DO 00953 rho_set%owns%drho=.TRUE. 00954 rho_set%has%drho=.TRUE. 00955 END IF 00956 CASE(2) 00957 IF (needs%rho) THEN 00958 ! this should basically never be the case unless you use LDA functionals 00959 ! with LSD 00960 00961 CALL pw_pool_create_cr3d(pw_pool,rho_set%rho,error=error) 00962 !assume that the bounds are the same? 00963 !$omp parallel do default(none) private(i,j,k) shared(rho_set,my_rho_r) 00964 DO k=rho_set%local_bounds(1,3),rho_set%local_bounds(2,3) 00965 DO j=rho_set%local_bounds(1,2),rho_set%local_bounds(2,2) 00966 DO i=rho_set%local_bounds(1,1),rho_set%local_bounds(2,1) 00967 rho_set%rho(i,j,k)=my_rho_r(1)%pw%cr3d(i,j,k)+& 00968 my_rho_r(2)%pw%cr3d(i,j,k) 00969 END DO 00970 END DO 00971 END DO 00972 rho_set%owns%rho=.TRUE. 00973 rho_set%has%rho=.TRUE. 00974 END IF 00975 IF (needs%rho_1_3) THEN 00976 CALL pw_pool_create_cr3d(pw_pool,rho_set%rho_1_3,error=error) 00977 !assume that the bounds are the same? 00978 !$omp parallel do default(none) private(i,j,k) shared(rho_set,my_rho_r) 00979 DO k=rho_set%local_bounds(1,3),rho_set%local_bounds(2,3) 00980 DO j=rho_set%local_bounds(1,2),rho_set%local_bounds(2,2) 00981 DO i=rho_set%local_bounds(1,1),rho_set%local_bounds(2,1) 00982 rho_set%rho_1_3(i,j,k)=MAX(my_rho_r(1)%pw%cr3d(i,j,k)+& 00983 my_rho_r(2)%pw%cr3d(i,j,k),0.0_dp)**f13 00984 END DO 00985 END DO 00986 END DO 00987 rho_set%owns%rho_1_3=.TRUE. 00988 rho_set%has%rho_1_3=.TRUE. 00989 END IF 00990 IF (needs%rho_spin_1_3) THEN 00991 CALL pw_pool_create_cr3d(pw_pool,rho_set%rhoa_1_3,error=error) 00992 !assume that the bounds are the same? 00993 !$omp parallel do default(none) private(i,j,k) shared(rho_set,my_rho_r) 00994 DO k=rho_set%local_bounds(1,3),rho_set%local_bounds(2,3) 00995 DO j=rho_set%local_bounds(1,2),rho_set%local_bounds(2,2) 00996 DO i=rho_set%local_bounds(1,1),rho_set%local_bounds(2,1) 00997 rho_set%rhoa_1_3(i,j,k)=MAX(my_rho_r(1)%pw%cr3d(i,j,k),0.0_dp)**f13 00998 END DO 00999 END DO 01000 END DO 01001 CALL pw_pool_create_cr3d(pw_pool,rho_set%rhob_1_3,error=error) 01002 !assume that the bounds are the same? 01003 !$omp parallel do default(none) private(i,j,k) shared(rho_set,my_rho_r) 01004 DO k=rho_set%local_bounds(1,3),rho_set%local_bounds(2,3) 01005 DO j=rho_set%local_bounds(1,2),rho_set%local_bounds(2,2) 01006 DO i=rho_set%local_bounds(1,1),rho_set%local_bounds(2,1) 01007 rho_set%rhob_1_3(i,j,k)=MAX(my_rho_r(2)%pw%cr3d(i,j,k),0.0_dp)**f13 01008 END DO 01009 END DO 01010 END DO 01011 rho_set%owns%rho_spin_1_3=.TRUE. 01012 rho_set%has%rho_spin_1_3=.TRUE. 01013 END IF 01014 IF (needs%rho_spin) THEN 01015 01016 rho_set%rhoa => my_rho_r(1)%pw%cr3d 01017 IF (my_rho_r_local) NULLIFY(my_rho_r(1)%pw%cr3d) 01018 01019 rho_set%rhob => my_rho_r(2)%pw%cr3d 01020 IF (my_rho_r_local) NULLIFY(my_rho_r(2)%pw%cr3d) 01021 01022 rho_set%owns%rho_spin=my_rho_r_local 01023 rho_set%has%rho_spin=.TRUE. 01024 END IF 01025 IF (needs%norm_drho) THEN 01026 01027 CALL pw_pool_create_cr3d(pw_pool, rho_set%norm_drho, & 01028 error=error) 01029 !$omp parallel do default(none) private(i,j,k) shared(rho_set,drho_r) 01030 DO k=rho_set%local_bounds(1,3),rho_set%local_bounds(2,3) 01031 DO j=rho_set%local_bounds(1,2),rho_set%local_bounds(2,2) 01032 DO i=rho_set%local_bounds(1,1),rho_set%local_bounds(2,1) 01033 rho_set%norm_drho(i,j,k)=SQRT(& 01034 (drho_r(1,1)%pw%cr3d(i,j,k)+drho_r(1,2)%pw%cr3d(i,j,k))**2+& 01035 (drho_r(2,1)%pw%cr3d(i,j,k)+drho_r(2,2)%pw%cr3d(i,j,k))**2+& 01036 (drho_r(3,1)%pw%cr3d(i,j,k)+drho_r(3,2)%pw%cr3d(i,j,k))**2) 01037 END DO 01038 END DO 01039 END DO 01040 01041 rho_set%owns%norm_drho=.TRUE. 01042 rho_set%has%norm_drho=.TRUE. 01043 END IF 01044 IF (needs%norm_drho_spin) THEN 01045 01046 CALL pw_pool_create_cr3d(pw_pool, rho_set%norm_drhoa, & 01047 error=error) 01048 !$omp parallel do default(none) private(i,j,k) shared(rho_set,drho_r) 01049 DO k=rho_set%local_bounds(1,3),rho_set%local_bounds(2,3) 01050 DO j=rho_set%local_bounds(1,2),rho_set%local_bounds(2,2) 01051 DO i=rho_set%local_bounds(1,1),rho_set%local_bounds(2,1) 01052 rho_set%norm_drhoa(i,j,k)=SQRT(& 01053 drho_r(1,1)%pw%cr3d(i,j,k)**2+& 01054 drho_r(2,1)%pw%cr3d(i,j,k)**2+& 01055 drho_r(3,1)%pw%cr3d(i,j,k)**2) 01056 END DO 01057 END DO 01058 END DO 01059 01060 CALL pw_pool_create_cr3d(pw_pool, rho_set%norm_drhob, & 01061 error=error) 01062 rho_set%owns%norm_drho_spin=.TRUE. 01063 !$omp parallel do default(none) private(i,j,k) shared(rho_set,drho_r) 01064 DO k=rho_set%local_bounds(1,3),rho_set%local_bounds(2,3) 01065 DO j=rho_set%local_bounds(1,2),rho_set%local_bounds(2,2) 01066 DO i=rho_set%local_bounds(1,1),rho_set%local_bounds(2,1) 01067 rho_set%norm_drhob(i,j,k)=SQRT(& 01068 drho_r(1,2)%pw%cr3d(i,j,k)**2+& 01069 drho_r(2,2)%pw%cr3d(i,j,k)**2+& 01070 drho_r(3,2)%pw%cr3d(i,j,k)**2) 01071 END DO 01072 END DO 01073 END DO 01074 01075 rho_set%owns%norm_drho_spin=.TRUE. 01076 rho_set%has%norm_drho_spin=.TRUE. 01077 END IF 01078 IF (needs%laplace_rho_spin) THEN 01079 CALL pw_pool_create_cr3d(pw_pool, rho_set%laplace_rhoa, & 01080 error=error) 01081 !$omp parallel do default(none) private(i,j,k) shared(rho_set,laplace_rho_r) 01082 DO k=rho_set%local_bounds(1,3),rho_set%local_bounds(2,3) 01083 DO j=rho_set%local_bounds(1,2),rho_set%local_bounds(2,2) 01084 DO i=rho_set%local_bounds(1,1),rho_set%local_bounds(2,1) 01085 rho_set%laplace_rhoa(i,j,k)=& 01086 laplace_rho_r(1,1)%pw%cr3d(i,j,k)+& 01087 laplace_rho_r(2,1)%pw%cr3d(i,j,k)+& 01088 laplace_rho_r(3,1)%pw%cr3d(i,j,k) 01089 END DO 01090 END DO 01091 END DO 01092 01093 CALL pw_pool_create_cr3d(pw_pool, rho_set%laplace_rhob, & 01094 error=error) 01095 rho_set%owns%laplace_rho_spin=.TRUE. 01096 !$omp parallel do default(none) private(i,j,k) shared(rho_set,laplace_rho_r) 01097 DO k=rho_set%local_bounds(1,3),rho_set%local_bounds(2,3) 01098 DO j=rho_set%local_bounds(1,2),rho_set%local_bounds(2,2) 01099 DO i=rho_set%local_bounds(1,1),rho_set%local_bounds(2,1) 01100 rho_set%laplace_rhob(i,j,k)=& 01101 laplace_rho_r(1,2)%pw%cr3d(i,j,k)+& 01102 laplace_rho_r(2,2)%pw%cr3d(i,j,k)+& 01103 laplace_rho_r(3,2)%pw%cr3d(i,j,k) 01104 END DO 01105 END DO 01106 END DO 01107 01108 rho_set%owns%laplace_rho_spin=.TRUE. 01109 rho_set%has%laplace_rho_spin=.TRUE. 01110 END IF 01111 IF (needs%drhoa_drhob) THEN 01112 CALL pw_pool_create_cr3d(pw_pool, rho_set%drhoa_drhob, & 01113 error=error) 01114 !$omp parallel do default(none) private(i,j,k) shared(rho_set,drho_r) 01115 DO k=rho_set%local_bounds(1,3),rho_set%local_bounds(2,3) 01116 DO j=rho_set%local_bounds(1,2),rho_set%local_bounds(2,2) 01117 DO i=rho_set%local_bounds(1,1),rho_set%local_bounds(2,1) 01118 rho_set%drhoa_drhob(i,j,k)=& 01119 drho_r(1,1)%pw%cr3d(i,j,k)*drho_r(1,2)%pw%cr3d(i,j,k)+& 01120 drho_r(2,1)%pw%cr3d(i,j,k)*drho_r(2,2)%pw%cr3d(i,j,k)+& 01121 drho_r(3,1)%pw%cr3d(i,j,k)*drho_r(3,2)%pw%cr3d(i,j,k) 01122 END DO 01123 END DO 01124 END DO 01125 rho_set%owns%drhoa_drhob=.TRUE. 01126 rho_set%has%drhoa_drhob=.TRUE. 01127 END IF 01128 IF (needs%drho) THEN 01129 ! this should basically never be the case unless you use LDA functionals 01130 ! with LSD 01131 DO idir=1,3 01132 CALL pw_pool_create_cr3d(pw_pool,rho_set%drho(idir)%array,& 01133 error=error) 01134 !assume that the bounds are the same? 01135 !$omp parallel do default(none) private(i,j,k) shared(rho_set,drho_r,idir) 01136 DO k=rho_set%local_bounds(1,3),rho_set%local_bounds(2,3) 01137 DO j=rho_set%local_bounds(1,2),rho_set%local_bounds(2,2) 01138 DO i=rho_set%local_bounds(1,1),rho_set%local_bounds(2,1) 01139 rho_set%drho(idir)%array(i,j,k)=& 01140 drho_r(idir,1)%pw%cr3d(i,j,k)+& 01141 drho_r(idir,2)%pw%cr3d(i,j,k) 01142 END DO 01143 END DO 01144 END DO 01145 END DO 01146 rho_set%owns%drho=.TRUE. 01147 rho_set%has%drho=.TRUE. 01148 END IF 01149 IF (needs%drho_spin) THEN 01150 DO idir=1,3 01151 rho_set%drhoa(idir)%array => drho_r(idir,1)%pw%cr3d 01152 NULLIFY(drho_r(idir,1)%pw%cr3d) 01153 rho_set%drhob(idir)%array => drho_r(idir,2)%pw%cr3d 01154 NULLIFY(drho_r(idir,2)%pw%cr3d) 01155 END DO 01156 rho_set%owns%drho_spin=.TRUE. 01157 rho_set%has%drho_spin=.TRUE. 01158 END IF 01159 END SELECT 01160 ! post cleanup 01161 DO ispin=1,nspins 01162 DO idir=1,3 01163 IF(needs%laplace_rho.OR.needs%laplace_rho_spin) THEN 01164 CALL pw_pool_give_back_pw(pw_pool, laplace_rho_r(idir,ispin)%pw, & 01165 accept_non_compatible=.TRUE., error=error) 01166 END IF 01167 CALL pw_pool_give_back_pw(pw_pool, drho_r(idir,ispin)%pw, & 01168 accept_non_compatible=.TRUE., error=error) 01169 END DO 01170 END DO 01171 IF (my_rho_r_local) THEN 01172 DO ispin=1,nspins 01173 CALL pw_pool_give_back_pw(pw_pool, my_rho_r(ispin)%pw,& 01174 accept_non_compatible=.TRUE., error=error) 01175 END DO 01176 END IF 01177 01178 ! tau part 01179 IF (needs%tau.OR.needs%tau_spin) THEN 01180 CPPrecondition(ASSOCIATED(tau),cp_failure_level,routineP,error,failure) 01181 DO ispin=1,nspins 01182 CPPrecondition(ASSOCIATED(tau(ispin)%pw),cp_failure_level,routineP,error,failure) 01183 CPPrecondition(ASSOCIATED(tau(ispin)%pw%cr3d),cp_failure_level,routineP,error,failure) 01184 END DO 01185 END IF 01186 IF (needs%tau) THEN 01187 IF (nspins==2) THEN 01188 CALL pw_pool_create_cr3d(pw_pool,rho_set%tau,& 01189 error=error) 01190 !$omp parallel do default(none) private(i,j,k) shared(rho_set,tau) 01191 DO k=rho_set%local_bounds(1,3),rho_set%local_bounds(2,3) 01192 DO j=rho_set%local_bounds(1,2),rho_set%local_bounds(2,2) 01193 DO i=rho_set%local_bounds(1,1),rho_set%local_bounds(2,1) 01194 01195 rho_set%tau(i,j,k)=& 01196 tau(1)%pw%cr3d(i,j,k)+& 01197 tau(2)%pw%cr3d(i,j,k) 01198 END DO 01199 END DO 01200 END DO 01201 rho_set%owns%tau=.TRUE. 01202 01203 ELSE 01204 rho_set%tau => tau(1)%pw%cr3d 01205 rho_set%owns%tau=.FALSE. 01206 END IF 01207 rho_set%has%tau=.TRUE. 01208 END IF 01209 IF (needs%tau_spin) THEN 01210 CPPrecondition(nspins==2,cp_failure_level,routineP,error,failure) 01211 rho_set%tau_a => tau(1)%pw%cr3d 01212 rho_set%tau_b => tau(2)%pw%cr3d 01213 rho_set%owns%tau_spin=.FALSE. 01214 rho_set%has%tau_spin=.TRUE. 01215 END IF 01216 01217 CPPostcondition(xc_rho_cflags_equal(rho_set%has,needs,error=error),cp_failure_level,routineP,error,failure) 01218 01219 END SUBROUTINE xc_rho_set_update 01220 01221 END MODULE xc_rho_set_types
1.7.3