|
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 00014 MODULE xc_xlda_hole_t_c_lr 00015 USE f77_blas 00016 USE input_section_types, ONLY: section_vals_type,& 00017 section_vals_val_get 00018 USE kinds, ONLY: dp 00019 USE mathconstants, ONLY: pi 00020 USE mathlib, ONLY: expint 00021 USE timings, ONLY: timeset,& 00022 timestop 00023 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 00024 xc_dset_get_derivative 00025 USE xc_derivative_types, ONLY: xc_derivative_get,& 00026 xc_derivative_type 00027 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 00028 USE xc_rho_set_types, ONLY: xc_rho_set_get,& 00029 xc_rho_set_type 00030 #include "cp_common_uses.h" 00031 00032 IMPLICIT NONE 00033 00034 PRIVATE 00035 00036 ! *** Global parameters *** 00037 00038 PUBLIC :: xlda_hole_t_c_lr_lda_eval, xlda_hole_t_c_lr_lda_info,& 00039 xlda_hole_t_c_lr_lsd_eval, xlda_hole_t_c_lr_lsd_info,& 00040 xlda_hole_t_c_lr_lda_calc_0 00041 00042 REAL(KIND=dp), PARAMETER :: A = 1.0161144_dp, 00043 B = -0.37170836_dp, 00044 C = -0.077215461_dp, 00045 D = 0.57786348_dp, 00046 E = -0.051955731_dp 00047 00048 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xlda_hole_t_c_lr' 00049 00050 CONTAINS 00051 00052 ! ***************************************************************************** 00066 SUBROUTINE xlda_hole_t_c_lr_lda_info ( xlda_params, reference, shortform, needs, max_deriv, error) 00067 TYPE(section_vals_type), POINTER :: xlda_params 00068 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 00069 TYPE(xc_rho_cflags_type), 00070 INTENT(inout), OPTIONAL :: needs 00071 INTEGER, INTENT(out), OPTIONAL :: max_deriv 00072 TYPE(cp_error_type), INTENT(inout) :: error 00073 00074 IF ( PRESENT ( reference ) ) THEN 00075 reference = "{LDA version}" 00076 END IF 00077 IF ( PRESENT ( shortform ) ) THEN 00078 shortform = "{LDA}" 00079 END IF 00080 IF ( PRESENT(needs) ) THEN 00081 needs%rho=.TRUE. 00082 END IF 00083 IF (PRESENT(max_deriv)) max_deriv=1 00084 00085 END SUBROUTINE xlda_hole_t_c_lr_lda_info 00086 00087 ! ***************************************************************************** 00101 SUBROUTINE xlda_hole_t_c_lr_lsd_info ( xlda_params, reference, shortform, needs, max_deriv, error) 00102 TYPE(section_vals_type), POINTER :: xlda_params 00103 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 00104 TYPE(xc_rho_cflags_type), 00105 INTENT(inout), OPTIONAL :: needs 00106 INTEGER, INTENT(out), OPTIONAL :: max_deriv 00107 TYPE(cp_error_type), INTENT(inout) :: error 00108 00109 IF ( PRESENT ( reference ) ) THEN 00110 reference = "{LSD version}" 00111 END IF 00112 IF ( PRESENT ( shortform ) ) THEN 00113 shortform = "{LSD}" 00114 END IF 00115 IF ( PRESENT(needs) ) THEN 00116 needs%rho_spin=.TRUE. 00117 END IF 00118 IF (PRESENT(max_deriv)) max_deriv=1 00119 00120 END SUBROUTINE xlda_hole_t_c_lr_lsd_info 00121 00122 ! ***************************************************************************** 00137 SUBROUTINE xlda_hole_t_c_lr_lda_eval ( rho_set, deriv_set, order, params, error ) 00138 00139 TYPE(xc_rho_set_type), POINTER :: rho_set 00140 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00141 INTEGER, INTENT(IN) :: order 00142 TYPE(section_vals_type), POINTER :: params 00143 TYPE(cp_error_type), INTENT(inout) :: error 00144 00145 CHARACTER(len=*), PARAMETER :: routineN = 'xlda_hole_t_c_lr_lda_eval', 00146 routineP = moduleN//':'//routineN 00147 00148 INTEGER :: handle, npoints, stat 00149 INTEGER, DIMENSION(:, :), POINTER :: bo 00150 LOGICAL :: failure 00151 REAL(kind=dp) :: epsilon_rho, R, sx 00152 REAL(kind=dp), DIMENSION(:, :, :), 00153 POINTER :: dummy, e_0, e_rho, rho 00154 TYPE(xc_derivative_type), POINTER :: deriv 00155 00156 CALL timeset(routineN,handle) 00157 failure=.FALSE. 00158 00159 NULLIFY(bo) 00160 00161 CALL section_vals_val_get(params,"SCALE_X",r_val=sx,error=error) 00162 CALL section_vals_val_get(params,"CUTOFF_RADIUS",r_val=R,error=error) 00163 00164 00165 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00166 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00167 CPPrecondition(ASSOCIATED(deriv_set),cp_failure_level,routineP,error,failure) 00168 CPPrecondition(deriv_set%ref_count>0,cp_failure_level,routineP,error,failure) 00169 00170 IF (.NOT.failure) THEN 00171 00172 CALL xc_rho_set_get(rho_set,rho=rho,& 00173 local_bounds=bo,rho_cutoff=epsilon_rho,& 00174 error=error) 00175 npoints=(bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1) 00176 00177 ! meaningful default for the arrays we don't need: let us make compiler 00178 ! and debugger happy... 00179 IF (cp_debug) THEN 00180 ALLOCATE(dummy(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3)),stat=stat) 00181 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00182 ELSE 00183 dummy=> rho 00184 END IF 00185 00186 e_0 => dummy 00187 e_rho => dummy 00188 00189 IF (order>=0) THEN 00190 deriv => xc_dset_get_derivative(deriv_set,"",& 00191 allocate_deriv=.TRUE., error=error) 00192 CALL xc_derivative_get(deriv,deriv_data=e_0,error=error) 00193 END IF 00194 IF (order>=1.OR.order==-1) THEN 00195 deriv => xc_dset_get_derivative(deriv_set,"(rho)",& 00196 allocate_deriv=.TRUE.,error=error) 00197 CALL xc_derivative_get(deriv,deriv_data=e_rho,error=error) 00198 END IF 00199 IF (order>1.OR.order<-1) THEN 00200 CALL cp_unimplemented_error(fromWhere=routineP, & 00201 message="derivatives bigger than 1 not implemented", & 00202 error=error, error_level=cp_failure_level) 00203 END IF 00204 00205 IF ( R == 0.0_dp ) THEN 00206 CALL cp_unimplemented_error(fromWhere=routineP, & 00207 message="Cutoff_Radius 0.0 not implemented", & 00208 error=error, error_level=cp_failure_level) 00209 END IF 00210 CALL xlda_hole_t_c_lr_lda_calc(npoints,order,rho=rho,& 00211 e_0=e_0,e_rho=e_rho,& 00212 epsilon_rho=epsilon_rho,& 00213 sx=sx,R=R,error=error) 00214 00215 IF (cp_debug) THEN 00216 DEALLOCATE(dummy,stat=stat) 00217 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00218 ELSE 00219 NULLIFY(dummy) 00220 END IF 00221 END IF 00222 CALL timestop(handle) 00223 00224 END SUBROUTINE xlda_hole_t_c_lr_lda_eval 00225 00226 ! ***************************************************************************** 00232 SUBROUTINE xlda_hole_t_c_lr_lda_calc(npoints, order, rho, e_0, e_rho,& 00233 epsilon_rho,sx, R,& 00234 error) 00235 00236 INTEGER, INTENT(in) :: npoints, order 00237 REAL(kind=dp), DIMENSION(1:npoints), 00238 INTENT(inout) :: rho, e_0, e_rho 00239 REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, R 00240 TYPE(cp_error_type), INTENT(inout) :: error 00241 00242 CHARACTER(len=*), PARAMETER :: routineN = 'xlda_hole_t_c_lr_lda_calc', 00243 routineP = moduleN//':'//routineN 00244 00245 INTEGER :: ip 00246 REAL(dp) :: my_rho 00247 00248 !$omp parallel do default(none) & 00249 !$omp shared(npoints, rho, epsilon_rho, order, e_0, e_rho) & 00250 !$omp shared(sx, r, error) & 00251 !$omp private(ip, my_rho) 00252 00253 DO ip =1,npoints 00254 my_rho = MAX(rho(ip),0.0_dp) 00255 IF(my_rho > epsilon_rho) THEN 00256 CALL xlda_hole_t_c_lr_lda_calc_0(order, my_rho, e_0(ip), e_rho(ip),& 00257 sx, R, error) 00258 END IF 00259 END DO 00260 00261 !$omp end parallel do 00262 00263 END SUBROUTINE xlda_hole_t_c_lr_lda_calc 00264 00265 ! ***************************************************************************** 00271 SUBROUTINE xlda_hole_t_c_lr_lda_calc_0(order, rho, e_0, e_rho,& 00272 sx, R,& 00273 error) 00274 INTEGER, INTENT(IN) :: order 00275 REAL(KIND=dp), INTENT(IN) :: rho 00276 REAL(kind=dp), INTENT(INOUT) :: e_0, e_rho 00277 REAL(KIND=dp), INTENT(IN) :: sx, R 00278 TYPE(cp_error_type), INTENT(inout) :: error 00279 00280 CHARACTER(len=*), PARAMETER :: routineN = 'xlda_hole_t_c_lr_lda_calc_0', 00281 routineP = moduleN//':'//routineN 00282 00283 REAL(KIND=dp) :: t1, t12, t14, t15, t19, t2, t22, t23, t24, t25, t3, t32, 00284 t33, t36, t4, t41, t46, t5, t6, t62, t64, t67, t68, t7, t82, t86, t9, 00285 t91, t95 00286 00287 IF( order >= 0 ) THEN 00288 t1 = rho ** 2 00289 t2 = t1 * pi 00290 t3 = 3 ** (0.1e1_dp / 0.3e1_dp) 00291 t4 = pi ** 2 00292 t5 = t4 * rho 00293 t6 = t5 ** (0.1e1_dp / 0.3e1_dp) 00294 t7 = t6 ** 2 00295 t9 = t3 / t7 00296 t12 = LOG(R * t3 * t6) 00297 t14 = R ** 2 00298 t15 = t14 ** 2 00299 t19 = 0.1e1_dp / D 00300 t22 = t3 ** 2 00301 t23 = t22 * t7 00302 t24 = D * t14 * t23 00303 t25 = EXP(-t24) 00304 t32 = 9 + 4 * A * t14 * t23 00305 t33 = LOG(t32) 00306 t36 = D ** 2 00307 t41 = expint(1, t24) 00308 t46 = 0.1e1_dp / t36 00309 t62 = LOG(0.2e1_dp) 00310 t64 = LOG(A) 00311 t67 = A * t12 + 0.3e1_dp / 0.2e1_dp * E * t15 * t3 * t6 * t5 * t19 * t25 & 00312 - A * t33 / 0.2e1_dp + E / t36 / D * t25 + A * t41 / 0.2e1_dp + E * t14 & 00313 * t22 * t7 * t46 * t25 + B * t19 * t25 / 0.2e1_dp + C * t46 * t25 / 0.2e1_dp & 00314 + C * t14 * t22 * t7 * t19 * t25 / 0.2e1_dp + A * t62 + A * t64 & 00315 / 0.2e1_dp 00316 t68 = t9 * t67 00317 e_0 = e_0 + (0.2e1_dp / 0.3e1_dp * t2 * t68) * sx 00318 END IF 00319 IF( order >=1 .OR. order ==-1) THEN 00320 t82 = A / rho 00321 t86 = t4 ** 2 00322 t91 = A ** 2 00323 t95 = 0.1e1_dp / t6 * t4 00324 e_rho = e_rho + (0.4e1_dp / 0.3e1_dp * rho * pi * t68 - 0.4e1_dp / 0.9e1_dp * t1 * t4 * pi & 00325 * t3 / t7 / t5 * t67 + 0.2e1_dp / 0.3e1_dp * t2 * t9 * (t82 / 0.3e1_dp - & 00326 0.3e1_dp * E * t15 * t14 * t86 * rho * t25 - 0.4e1_dp / 0.3e1_dp * t91 * t14 & 00327 * t22 * t95 / t32 - t82 * t25 / 0.3e1_dp - B * t14 * t22 * t95 * t25 & 00328 / 0.3e1_dp - C * t15 * t3 * t6 * t4 * t25) ) * sx 00329 END IF 00330 00331 00332 END SUBROUTINE xlda_hole_t_c_lr_lda_calc_0 00333 00334 ! ***************************************************************************** 00350 SUBROUTINE xlda_hole_t_c_lr_lsd_eval ( rho_set, deriv_set, order, params, error ) 00351 00352 TYPE(xc_rho_set_type), POINTER :: rho_set 00353 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00354 INTEGER, INTENT(IN) :: order 00355 TYPE(section_vals_type), POINTER :: params 00356 TYPE(cp_error_type), INTENT(inout) :: error 00357 00358 CHARACTER(len=*), PARAMETER :: routineN = 'xlda_hole_t_c_lr_lsd_eval', 00359 routineP = moduleN//':'//routineN 00360 00361 INTEGER :: handle, npoints, stat 00362 INTEGER, DIMENSION(:, :), POINTER :: bo 00363 LOGICAL :: failure 00364 REAL(kind=dp) :: epsilon_rho, R, sx 00365 REAL(kind=dp), DIMENSION(:, :, :), 00366 POINTER :: dummy, e_0, e_rhoa, e_rhob, 00367 rhoa, rhob 00368 TYPE(xc_derivative_type), POINTER :: deriv 00369 00370 CALL timeset(routineN,handle) 00371 failure=.FALSE. 00372 00373 NULLIFY(bo) 00374 00375 CALL section_vals_val_get(params,"SCALE_X",r_val=sx,error=error) 00376 CALL section_vals_val_get(params,"CUTOFF_RADIUS",r_val=R,error=error) 00377 00378 00379 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00380 CPPrecondition(rho_set%ref_count>0,cp_failure_level,routineP,error,failure) 00381 CPPrecondition(ASSOCIATED(deriv_set),cp_failure_level,routineP,error,failure) 00382 CPPrecondition(deriv_set%ref_count>0,cp_failure_level,routineP,error,failure) 00383 00384 IF (.NOT.failure) THEN 00385 00386 CALL xc_rho_set_get(rho_set,rhoa=rhoa, rhob=rhob,& 00387 local_bounds=bo,rho_cutoff=epsilon_rho,& 00388 error=error) 00389 npoints=(bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1) 00390 00391 ! meaningful default for the arrays we don't need: let us make compiler 00392 ! and debugger happy... 00393 IF (cp_debug) THEN 00394 ALLOCATE(dummy(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3)),stat=stat) 00395 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00396 ELSE 00397 dummy=> rhoa 00398 END IF 00399 00400 e_0 => dummy 00401 e_rhoa => dummy 00402 e_rhob => dummy 00403 00404 IF (order>=0) THEN 00405 deriv => xc_dset_get_derivative(deriv_set,"",& 00406 allocate_deriv=.TRUE., error=error) 00407 CALL xc_derivative_get(deriv,deriv_data=e_0,error=error) 00408 END IF 00409 IF (order>=1.OR.order==-1) THEN 00410 deriv => xc_dset_get_derivative(deriv_set,"(rhoa)",& 00411 allocate_deriv=.TRUE.,error=error) 00412 CALL xc_derivative_get(deriv,deriv_data=e_rhoa,error=error) 00413 deriv => xc_dset_get_derivative(deriv_set,"(rhob)",& 00414 allocate_deriv=.TRUE.,error=error) 00415 CALL xc_derivative_get(deriv,deriv_data=e_rhob,error=error) 00416 END IF 00417 IF (order>1.OR.order<-1) THEN 00418 CALL cp_unimplemented_error(fromWhere=routineP, & 00419 message="derivatives bigger than 2 not implemented", & 00420 error=error, error_level=cp_failure_level) 00421 END IF 00422 IF ( R == 0.0_dp ) THEN 00423 CALL cp_unimplemented_error(fromWhere=routineP, & 00424 message="Cutoff_Radius 0.0 not implemented", & 00425 error=error, error_level=cp_failure_level) 00426 END IF 00427 00428 !$omp parallel default(none) & 00429 !$omp shared(npoints, order, rhoa, e_0, e_rhoa, epsilon_rho) & 00430 !$omp shared(sx, r, error, rhob, e_rhob) 00431 00432 CALL xlda_hole_t_c_lr_lsd_calc(npoints,order,rho=rhoa,& 00433 e_0=e_0,e_rho=e_rhoa,& 00434 epsilon_rho=epsilon_rho,& 00435 sx=sx,R=R,error=error) 00436 00437 CALL xlda_hole_t_c_lr_lsd_calc(npoints,order,rho=rhob,& 00438 e_0=e_0,e_rho=e_rhob,& 00439 epsilon_rho=epsilon_rho,& 00440 sx=sx,R=R,error=error) 00441 !$omp end parallel 00442 00443 IF (cp_debug) THEN 00444 DEALLOCATE(dummy,stat=stat) 00445 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00446 ELSE 00447 NULLIFY(dummy) 00448 END IF 00449 END IF 00450 CALL timestop(handle) 00451 00452 END SUBROUTINE xlda_hole_t_c_lr_lsd_eval 00453 00454 ! ***************************************************************************** 00460 SUBROUTINE xlda_hole_t_c_lr_lsd_calc(npoints, order, rho, e_0, e_rho,& 00461 epsilon_rho,sx, R,& 00462 error) 00463 00464 INTEGER, INTENT(in) :: npoints, order 00465 REAL(kind=dp), DIMENSION(1:npoints), 00466 INTENT(inout) :: rho, e_0, e_rho 00467 REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, R 00468 TYPE(cp_error_type), INTENT(inout) :: error 00469 00470 CHARACTER(len=*), PARAMETER :: routineN = 'xlda_hole_t_c_lr_lsd_calc', 00471 routineP = moduleN//':'//routineN 00472 00473 INTEGER :: ip 00474 REAL(dp) :: e_tmp, my_rho 00475 00476 !$omp do 00477 00478 DO ip =1,npoints 00479 my_rho = 2.0_dp*MAX(rho(ip),0.0_dp) 00480 IF(my_rho > epsilon_rho) THEN 00481 e_tmp = 0.0_dp 00482 CALL xlda_hole_t_c_lr_lda_calc_0(order, my_rho, e_tmp , e_rho(ip),& 00483 sx, R, error) 00484 e_0(ip) = e_0(ip) + 0.5_dp * e_tmp 00485 END IF 00486 END DO 00487 00488 !$omp end do 00489 00490 END SUBROUTINE xlda_hole_t_c_lr_lsd_calc 00491 END MODULE xc_xlda_hole_t_c_lr 00492
1.7.3