|
CP2K 2.5 (Revision 12981)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00015 MODULE xc 00016 USE cell_types, ONLY: cell_type 00017 USE cp_array_r_utils, ONLY: cp_3d_r_p_type 00018 USE cp_linked_list_xc_deriv, ONLY: cp_sll_xc_deriv_next,& 00019 cp_sll_xc_deriv_type 00020 USE input_constants, ONLY: & 00021 xc_debug_new_routine, xc_deriv_nn10_smooth, xc_deriv_nn50_smooth, & 00022 xc_deriv_pw, xc_deriv_spline2, xc_deriv_spline2_smooth, & 00023 xc_deriv_spline3, xc_deriv_spline3_smooth, xc_new_f_routine, & 00024 xc_rho_nn10, xc_rho_nn50, xc_rho_no_smooth, xc_rho_spline2_smooth, & 00025 xc_rho_spline3_smooth, xc_test_lsd_f_routine 00026 USE input_section_types, ONLY: section_get_ival,& 00027 section_get_rval,& 00028 section_vals_get_subs_vals,& 00029 section_vals_type,& 00030 section_vals_val_get 00031 USE kahan_sum, ONLY: accurate_sum 00032 USE kinds, ONLY: default_path_length,& 00033 dp 00034 USE message_passing, ONLY: mp_sum 00035 USE pw_grid_types, ONLY: PW_MODE_DISTRIBUTED,& 00036 pw_grid_type 00037 USE pw_methods, ONLY: pw_axpy,& 00038 pw_copy,& 00039 pw_derive,& 00040 pw_transfer,& 00041 pw_zero 00042 USE pw_pool_types, ONLY: pw_pool_create_pw,& 00043 pw_pool_give_back_cr3d,& 00044 pw_pool_give_back_pw,& 00045 pw_pool_type 00046 USE pw_spline_utils, ONLY: & 00047 nn10_coeffs, nn10_deriv_coeffs, nn50_coeffs, nn50_deriv_coeffs, & 00048 pw_nn_deriv_r, pw_nn_smear_r, pw_spline2_deriv_g, & 00049 pw_spline2_interpolate_values_g, pw_spline3_deriv_g, & 00050 pw_spline3_interpolate_values_g, pw_spline_scale_deriv, & 00051 spline2_coeffs, spline2_deriv_coeffs, spline3_coeffs, & 00052 spline3_deriv_coeffs 00053 USE pw_types, ONLY: COMPLEXDATA1D,& 00054 REALDATA3D,& 00055 REALSPACE,& 00056 RECIPROCALSPACE,& 00057 pw_create,& 00058 pw_p_type,& 00059 pw_release,& 00060 pw_type 00061 USE timings, ONLY: timeset,& 00062 timestop 00063 USE virial_types, ONLY: virial_type 00064 USE xc_derivative_desc, ONLY: MAX_DERIVATIVE_DESC_LENGTH,& 00065 MAX_LABEL_LENGTH 00066 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 00067 xc_dset_create,& 00068 xc_dset_get_derivative,& 00069 xc_dset_release 00070 USE xc_derivative_types, ONLY: xc_derivative_get,& 00071 xc_derivative_type 00072 USE xc_derivatives, ONLY: xc_functionals_eval,& 00073 xc_functionals_get_needs 00074 USE xc_rho_set_types, ONLY: xc_rho_set_create,& 00075 xc_rho_set_get,& 00076 xc_rho_set_release,& 00077 xc_rho_set_type,& 00078 xc_rho_set_update 00079 #include "cp_common_uses.h" 00080 00081 IMPLICIT NONE 00082 PRIVATE 00083 PUBLIC :: xc_vxc_pw_create1, xc_vxc_pw_create, & 00084 xc_rho_set_and_dset_create, xc_exc_calc,& 00085 xc_calc_2nd_deriv, xc_prep_2nd_deriv, divide_by_norm_drho 00086 00087 LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE. 00088 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc' 00089 00090 CONTAINS 00091 00092 ! ***************************************************************************** 00118 SUBROUTINE xc_vxc_pw_create1(vxc_rho,vxc_tau,rho_r,rho_g,tau,exc,xc_section,& 00119 cell,pw_pool,error,virial) 00120 TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau, rho_r, 00121 rho_g, tau 00122 REAL(KIND=dp), INTENT(out) :: exc 00123 TYPE(section_vals_type), POINTER :: xc_section 00124 TYPE(cell_type), POINTER :: cell 00125 TYPE(pw_pool_type), POINTER :: pw_pool 00126 TYPE(cp_error_type), INTENT(inout) :: error 00127 TYPE(virial_type), OPTIONAL, POINTER :: virial 00128 00129 CHARACTER(len=*), PARAMETER :: routineN = 'xc_vxc_pw_create1', 00130 routineP = moduleN//':'//routineN 00131 00132 INTEGER :: f_routine 00133 LOGICAL :: failure 00134 00135 failure=.FALSE. 00136 00137 CPPrecondition(ASSOCIATED(rho_r),cp_failure_level,routineP,error,failure) 00138 CPPrecondition(ASSOCIATED(xc_section),cp_failure_level,routineP,error,failure) 00139 CPPrecondition(.NOT.ASSOCIATED(vxc_rho),cp_failure_level,routineP,error,failure) 00140 CPPrecondition(.NOT.ASSOCIATED(vxc_tau),cp_failure_level,routineP,error,failure) 00141 IF (.NOT.failure) THEN 00142 00143 CALL section_vals_val_get(xc_section,"FUNCTIONAL_ROUTINE",& 00144 i_val=f_routine,error=error) 00145 SELECT CASE(f_routine) 00146 CASE(xc_new_f_routine) 00147 CALL xc_vxc_pw_create(vxc_rho=vxc_rho, vxc_tau=vxc_tau,tau=tau,& 00148 rho_r=rho_r, rho_g=rho_g, exc=exc, xc_section=xc_section,& 00149 cell=cell, pw_pool=pw_pool,error=error,virial=virial) 00150 CASE(xc_debug_new_routine) 00151 CALL xc_vxc_pw_create_debug(vxc_rho=vxc_rho, vxc_tau=vxc_tau,tau=tau,& 00152 rho_r=rho_r, rho_g=rho_g, exc=exc, xc_section=xc_section,& 00153 cell=cell, pw_pool=pw_pool, error=error) 00154 CASE(xc_test_lsd_f_routine) 00155 CALL xc_vxc_pw_create_test_lsd(vxc_rho=vxc_rho, vxc_tau=vxc_tau,& 00156 tau=tau, rho_r=rho_r, rho_g=rho_g, exc=exc, & 00157 xc_section=xc_section, cell=cell, pw_pool=pw_pool,& 00158 error=error) 00159 CASE default 00160 END SELECT 00161 END IF 00162 00163 END SUBROUTINE xc_vxc_pw_create1 00164 00165 ! ***************************************************************************** 00188 SUBROUTINE xc_vxc_pw_create_test_lsd(vxc_rho,vxc_tau,rho_r,rho_g,tau,& 00189 exc,xc_section, cell,pw_pool, error) 00190 TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau, rho_r, 00191 rho_g, tau 00192 REAL(KIND=dp), INTENT(out) :: exc 00193 TYPE(section_vals_type), POINTER :: xc_section 00194 TYPE(cell_type), POINTER :: cell 00195 TYPE(pw_pool_type), POINTER :: pw_pool 00196 TYPE(cp_error_type), INTENT(inout) :: error 00197 00198 CHARACTER(len=*), PARAMETER :: routineN = 'xc_vxc_pw_create_test_lsd', 00199 routineP = moduleN//':'//routineN 00200 00201 CHARACTER(len=default_path_length) :: filename 00202 CHARACTER(len=MAX_LABEL_LENGTH), 00203 DIMENSION(:), POINTER :: split_desc 00204 INTEGER :: i, ii, ispin, j, k, stat 00205 INTEGER, DIMENSION(2, 3) :: bo 00206 LOGICAL :: failure 00207 REAL(kind=dp) :: diff, exc2, maxdiff, tmp 00208 REAL(kind=dp), DIMENSION(:, :, :), 00209 POINTER :: pot, pot2, pot3 00210 TYPE(cp_sll_xc_deriv_type), POINTER :: deriv_iter 00211 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho2_g, rho2_r, tau2, 00212 vxc_rho2, vxc_tau2 00213 TYPE(xc_derivative_set_type), POINTER :: dSet1, dSet2 00214 TYPE(xc_derivative_type), POINTER :: deriv, deriv2, deriv3 00215 TYPE(xc_rho_set_type), POINTER :: rho_set1, rho_set2 00216 00217 failure=.FALSE. 00218 NULLIFY(vxc_rho2,vxc_tau2,tau2,dSet1,dSet2,rho_set1,rho_set2,split_desc,pot,pot3,pot3,& 00219 deriv,deriv2,deriv3,rho2_g) 00220 00221 IF (.NOT. failure) THEN 00222 bo = rho_r(1)%pw%pw_grid%bounds_local 00223 00224 ALLOCATE(rho2_r(2), stat=stat) 00225 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00226 DO ispin=1,2 00227 NULLIFY(rho2_r(ispin)%pw) 00228 CALL pw_pool_create_pw(pw_pool,rho2_r(ispin)%pw,in_space=REALSPACE,& 00229 use_data=REALDATA3D, error=error) 00230 END DO 00231 DO k=bo(1,3),bo(2,3) 00232 DO j=bo(1,2),bo(2,2) 00233 DO i=bo(1,1),bo(2,1) 00234 tmp=rho_r(1)%pw%cr3d(i,j,k)*0.5 00235 rho2_r(1)%pw%cr3d(i,j,k)=tmp 00236 rho2_r(2)%pw%cr3d(i,j,k)=tmp 00237 END DO 00238 END DO 00239 END DO 00240 00241 IF (ASSOCIATED(tau)) THEN 00242 ALLOCATE(tau2(2),stat=stat) 00243 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00244 DO ispin=1,2 00245 NULLIFY(tau2(ispin)%pw) 00246 CALL pw_pool_create_pw(pw_pool,tau2(ispin)%pw,in_space=REALSPACE,& 00247 use_data=REALDATA3D, error=error) 00248 END DO 00249 00250 DO k=bo(1,3),bo(2,3) 00251 DO j=bo(1,2),bo(2,2) 00252 DO i=bo(1,1),bo(2,1) 00253 tmp=tau(1)%pw%cr3d(i,j,k)*0.5 00254 tau2(1)%pw%cr3d(i,j,k)=tmp 00255 tau2(2)%pw%cr3d(i,j,k)=tmp 00256 END DO 00257 END DO 00258 END DO 00259 END IF 00260 00261 PRINT *, "about to calculate xc (lda)" 00262 CALL xc_rho_set_and_dset_create(rho_r=rho_r, rho_g=rho_g,& 00263 tau=tau,xc_section=xc_section,& 00264 cell=cell, pw_pool=pw_pool,rho_set=rho_set1,& 00265 deriv_set=dSet1, deriv_order=1,& 00266 needs_basic_components=.FALSE.,error=error) 00267 CALL xc_vxc_pw_create(rho_r=rho_r, rho_g=rho_g,tau=tau,& 00268 vxc_rho=vxc_rho,vxc_tau=vxc_tau, exc=exc, xc_section=xc_section,& 00269 cell=cell, pw_pool=pw_pool, & 00270 error=error) 00271 PRINT *, "did calculate xc (lda)" 00272 PRINT *, "about to calculate xc (lsd)" 00273 CALL xc_rho_set_and_dset_create(rho_set=rho_set2,deriv_set=dSet2,& 00274 rho_r=rho2_r, rho_g=rho2_g,tau=tau2, xc_section=xc_section,& 00275 cell=cell, pw_pool=pw_pool, deriv_order=1,& 00276 needs_basic_components=.FALSE.,error=error) 00277 CALL xc_vxc_pw_create(rho_r=rho2_r, rho_g=rho2_g,tau=tau2,& 00278 vxc_rho=vxc_rho2,vxc_tau=vxc_tau2,exc=exc2, xc_section=xc_section,& 00279 cell=cell, pw_pool=pw_pool, error=error) 00280 PRINT *, "did calculate xc (new)" 00281 PRINT *, "at (0,0,0) rho_r=",rho_r(1)%pw%cr3d(0,0,0),& 00282 "rho2_r(1)=",rho2_r(1)%pw%cr3d(0,0,0),& 00283 "rho2_r(2)=",rho2_r(2)%pw%cr3d(0,0,0),& 00284 "rho_r_sm=",rho_set1%rho(0,0,0), "rhoa2_r_sm=",rho_set2%rhoa(0,0,0),& 00285 "rhob2_r_sm=",rho_set2%rhob(0,0,0) 00286 OPEN(unit=120,file="rho.bindata",status="unknown",access='sequential',& 00287 form="unformatted",action="write") 00288 pot => rho_set1%rho 00289 WRITE(unit=120) pot(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(2,3)) 00290 CLOSE(unit=120) 00291 OPEN(unit=120,file="rhoa.bindata",status="unknown",access='sequential',& 00292 form="unformatted",action="write") 00293 pot => rho_set2%rhoa 00294 WRITE(unit=120) pot(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(2,3)) 00295 CLOSE(unit=120) 00296 OPEN(unit=120,file="rhob.bindata",status="unknown",access='sequential',& 00297 form="unformatted",action="write") 00298 pot => rho_set2%rhob 00299 WRITE(unit=120) pot(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(2,3)) 00300 CLOSE(unit=120) 00301 OPEN(unit=120,file="ndrho.bindata",status="unknown",access='sequential',& 00302 form="unformatted",action="write") 00303 pot => rho_set1%norm_drho 00304 WRITE(unit=120) pot(:,:,bo(2,3)) 00305 CLOSE(unit=120) 00306 OPEN(unit=120,file="ndrhoa.bindata",status="unknown",access='sequential',& 00307 form="unformatted",action="write") 00308 pot => rho_set2%norm_drhoa 00309 WRITE(unit=120) pot(:,:,bo(2,3)) 00310 CLOSE(unit=120) 00311 OPEN(unit=120,file="ndrhob.bindata",status="unknown",access='sequential',& 00312 form="unformatted",action="write") 00313 pot => rho_set2%norm_drhob 00314 WRITE(unit=120) pot(:,:,bo(2,3)) 00315 CLOSE(unit=120) 00316 IF (rho_set1%has%tau) THEN 00317 OPEN(unit=120,file="tau.bindata",status="unknown",access='sequential',& 00318 form="unformatted",action="write") 00319 pot => rho_set1%tau 00320 WRITE(unit=120) pot(:,:,bo(2,3)) 00321 CLOSE(unit=120) 00322 END IF 00323 IF (rho_set2%has%tau_spin) THEN 00324 OPEN(unit=120,file="tau_a.bindata",status="unknown",access='sequential',& 00325 form="unformatted",action="write") 00326 pot => rho_set2%tau_a 00327 WRITE(unit=120) pot(:,:,bo(2,3)) 00328 CLOSE(unit=120) 00329 OPEN(unit=120,file="tau_v.bindata",status="unknown",access='sequential',& 00330 form="unformatted",action="write") 00331 pot => rho_set2%tau_b 00332 WRITE(unit=120) pot(:,:,bo(2,3)) 00333 CLOSE(unit=120) 00334 END IF 00335 OPEN(unit=120,file="vxc.bindata",status="unknown",access='sequential',& 00336 form="unformatted",action="write") 00337 pot => vxc_rho(1)%pw%cr3d 00338 WRITE(unit=120) pot(:,:,bo(2,3)) 00339 CLOSE(unit=120) 00340 OPEN(unit=120,file="vxc2.bindata",status="unknown",access='sequential',& 00341 form="unformatted",action="write") 00342 pot => vxc_rho2(1)%pw%cr3d 00343 WRITE(unit=120) pot(:,:,bo(2,3)) 00344 CLOSE(unit=120) 00345 IF (ASSOCIATED(vxc_tau)) THEN 00346 OPEN(unit=120,file="vxc_tau.bindata",status="unknown",access='sequential',& 00347 form="unformatted",action="write") 00348 pot => vxc_tau(1)%pw%cr3d 00349 WRITE(unit=120) pot(:,:,bo(2,3)) 00350 CLOSE(unit=120) 00351 END IF 00352 IF (ASSOCIATED(vxc_tau2)) THEN 00353 OPEN(unit=120,file="vxc_tau2_a.bindata",status="unknown",access='sequential',& 00354 form="unformatted",action="write") 00355 pot => vxc_tau2(1)%pw%cr3d 00356 WRITE(unit=120) pot(:,:,bo(2,3)) 00357 CLOSE(unit=120) 00358 OPEN(unit=120,file="vxc_tau2_b.bindata",status="unknown",access='sequential',& 00359 form="unformatted",action="write") 00360 pot => vxc_tau2(2)%pw%cr3d 00361 WRITE(unit=120) pot(:,:,bo(2,3)) 00362 CLOSE(unit=120) 00363 END IF 00364 00365 PRINT *,"calc diff on vxc" 00366 maxDiff=0.0_dp 00367 DO ispin=1,1 00368 ii=0 00369 DO k=bo(1,3),bo(2,3) 00370 DO j=bo(1,2),bo(2,2) 00371 DO i=bo(1,1),bo(2,1) 00372 ii=ii+1 00373 diff=ABS(vxc_rho(ispin)%pw%cr3d(i,j,k)-& 00374 vxc_rho2(ispin)%pw%cr3d(i,j,k)) 00375 IF (ii==1) THEN 00376 PRINT *,"vxc",ispin,"=",vxc_rho(ispin)%pw%cr3d(i,j,k),"vs",vxc_rho2(ispin)%pw%cr3d(i,j,k),"diff=",diff 00377 END IF 00378 IF (maxDiff<diff)THEN 00379 maxDiff=diff 00380 PRINT *, "diff=",diff," at ",i,",",j,",",k,& 00381 " spin=",ispin,"rho=",rho_set1%rho(i,j,k),& 00382 " ndrho=",rho_set1%norm_drho(i,j,k) 00383 END IF 00384 END DO 00385 END DO 00386 END DO 00387 END DO 00388 PRINT *,"diff exc=",ABS(exc-exc2),"diff vxc=",maxdiff 00389 ! CPPostcondition(maxdiff<5.e-11,cp_failure_level,routineP,error,failure) 00390 ! CPPostcondition(ABS(exc-exc2)<1.e-14,cp_failure_level,routineP,error,failure) 00391 00392 IF (ASSOCIATED(vxc_tau)) THEN 00393 PRINT *,"calc diff on vxc_tau" 00394 maxDiff=0.0_dp 00395 DO ispin=1,1 00396 ii=0 00397 DO k=bo(1,3),bo(2,3) 00398 DO j=bo(1,2),bo(2,2) 00399 DO i=bo(1,1),bo(2,1) 00400 ii=ii+1 00401 diff=ABS(vxc_tau(ispin)%pw%cr3d(i,j,k)-& 00402 vxc_tau2(ispin)%pw%cr3d(i,j,k)) 00403 IF (ii==1) THEN 00404 PRINT *,"vxc_tau",ispin,"=",vxc_tau(ispin)%pw%cr3d(i,j,k),"vs",vxc_tau2(ispin)%pw%cr3d(i,j,k),"diff=",diff 00405 END IF 00406 IF (maxDiff<diff)THEN 00407 maxDiff=diff 00408 PRINT *, "diff=",diff," at ",i,",",j,",",k,& 00409 " spin=",ispin,"rho=",rho_set1%rho(i,j,k),& 00410 " ndrho=",rho_set1%norm_drho(i,j,k) 00411 END IF 00412 END DO 00413 END DO 00414 END DO 00415 END DO 00416 PRINT *,"diff exc=",ABS(exc-exc2),"diff vxc_tau=",maxdiff 00417 END IF 00418 deriv_iter => dSet1%derivs 00419 DO WHILE (cp_sll_xc_deriv_next(deriv_iter,el_att=deriv,error=error)) 00420 CALL xc_derivative_get(deriv,& 00421 split_desc=split_desc,deriv_data=pot,& 00422 error=error) 00423 SELECT CASE (SIZE(split_desc)) 00424 CASE(0) 00425 filename="e_0.bindata" 00426 deriv2 => xc_dset_get_derivative(dSet2, "", error=error) 00427 CASE(1) 00428 filename="e_"//TRIM(split_desc(1))//".bindata" 00429 IF (split_desc(1)=="rho") THEN 00430 deriv2 => xc_dset_get_derivative(dSet2, "(rhoa)", error=error) 00431 ELSEIF (split_desc(1)=="tau") THEN 00432 deriv2 => xc_dset_get_derivative(dSet2,"(tau_a)",error=error) 00433 ELSEIF (split_desc(1)=="norm_drho") THEN 00434 deriv2 => xc_dset_get_derivative(dSet2, "(norm_drhoa)", error=error) 00435 deriv3 => xc_dset_get_derivative(dSet2, "(norm_drho)", error=error) 00436 IF (ASSOCIATED(deriv3)) THEN 00437 IF (ASSOCIATED(deriv2)) THEN 00438 CALL xc_derivative_get(deriv2,& 00439 deriv_data=pot2,& 00440 error=error) 00441 CALL xc_derivative_get(deriv3,& 00442 deriv_data=pot3,& 00443 error=error) 00444 pot2=pot2+pot3 00445 ELSE 00446 deriv2 => deriv3 00447 END IF 00448 NULLIFY(deriv3,pot2,pot3) 00449 END IF 00450 ELSE 00451 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00452 END IF 00453 CASE default 00454 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00455 END SELECT 00456 CALL xc_derivative_get(deriv2,& 00457 deriv_data=pot2,& 00458 error=error) 00459 PRINT *, "checking ",filename 00460 maxDiff=0.0_dp 00461 DO k=bo(1,3),bo(2,3) 00462 DO j=bo(1,2),bo(2,2) 00463 DO i=bo(1,1),bo(2,1) 00464 diff=ABS(pot(i,j,k)-pot2(i,j,k)) 00465 IF (maxDiff<diff) THEN 00466 maxDiff=diff 00467 PRINT *, "ediff(",i,j,k,")=",maxDiff,& 00468 "rho=",rho_set1%rho(i,j,k),& 00469 "ndrho=",rho_set1%norm_drho(i,j,k) 00470 END IF 00471 END DO 00472 END DO 00473 END DO 00474 PRINT *,"maxdiff ",filename,"=",maxDiff 00475 OPEN (unit=120,file=TRIM(filename),status="unknown",& 00476 access='sequential',& 00477 form="unformatted") 00478 WRITE (unit=120) pot(:,:,bo(2,3)) 00479 CLOSE (unit=120) 00480 END DO 00481 deriv_iter => dSet2%derivs 00482 DO WHILE (cp_sll_xc_deriv_next(deriv_iter,el_att=deriv,error=error)) 00483 CALL xc_derivative_get(deriv,& 00484 split_desc=split_desc,deriv_data=pot,& 00485 error=error) 00486 SELECT CASE (SIZE(split_desc)) 00487 CASE(0) 00488 filename="e_0-2.bindata" 00489 CASE(1) 00490 filename="e_"//TRIM(split_desc(1))//"-2.bindata" 00491 CASE default 00492 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00493 END SELECT 00494 OPEN (unit=120,file=TRIM(filename),status="unknown",& 00495 access='sequential',& 00496 form="unformatted") 00497 WRITE (unit=120) pot(:,:,bo(2,3)) 00498 CLOSE (unit=120) 00499 END DO 00500 CALL xc_rho_set_release(rho_set1,error=error) 00501 CALL xc_rho_set_release(rho_set2,error=error) 00502 CALL xc_dset_release(dSet2,error=error) 00503 CALL xc_dset_release(dSet1, error=error) 00504 DO ispin=1,2 00505 CALL pw_pool_give_back_pw(pw_pool,rho2_r(ispin)%pw,& 00506 error=error) 00507 CALL pw_pool_give_back_pw(pw_pool,vxc_rho2(ispin)%pw,& 00508 error=error) 00509 IF (ASSOCIATED(vxc_tau2)) THEN 00510 CALL pw_pool_give_back_pw(pw_pool,vxc_tau2(ispin)%pw,& 00511 error=error) 00512 END IF 00513 END DO 00514 DEALLOCATE(vxc_rho2,rho2_r,rho2_g, stat=stat) 00515 CPPostcondition(stat==0,cp_warning_level,routineP,error,failure) 00516 IF (ASSOCIATED(vxc_tau2)) THEN 00517 DEALLOCATE(vxc_tau2,stat=stat) 00518 CPPostcondition(stat==0,cp_warning_level,routineP,error,failure) 00519 END IF 00520 00521 END IF 00522 END SUBROUTINE xc_vxc_pw_create_test_lsd 00523 00524 ! ***************************************************************************** 00547 SUBROUTINE xc_vxc_pw_create_debug(vxc_rho,vxc_tau,rho_r,rho_g,tau,exc,& 00548 xc_section, cell,pw_pool,error) 00549 TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau, rho_r, 00550 rho_g, tau 00551 REAL(KIND=dp), INTENT(out) :: exc 00552 TYPE(section_vals_type), POINTER :: xc_section 00553 TYPE(cell_type), POINTER :: cell 00554 TYPE(pw_pool_type), POINTER :: pw_pool 00555 TYPE(cp_error_type), INTENT(inout) :: error 00556 00557 CHARACTER(len=*), PARAMETER :: routineN = 'xc_vxc_pw_create_debug', 00558 routineP = moduleN//':'//routineN 00559 00560 CHARACTER(len=default_path_length) :: filename 00561 CHARACTER(len=MAX_LABEL_LENGTH), 00562 DIMENSION(:), POINTER :: split_desc 00563 INTEGER :: i, ispin, j, k 00564 INTEGER, DIMENSION(2, 3) :: bo 00565 LOGICAL :: failure 00566 REAL(kind=dp), DIMENSION(:, :, :), 00567 POINTER :: pot 00568 TYPE(cp_logger_type), POINTER :: logger 00569 TYPE(cp_sll_xc_deriv_type), POINTER :: deriv_iter 00570 TYPE(xc_derivative_set_type), POINTER :: dSet1 00571 TYPE(xc_derivative_type), POINTER :: deriv 00572 TYPE(xc_rho_set_type), POINTER :: rho_set1 00573 00574 failure=.FALSE. 00575 NULLIFY(dSet1,rho_set1,split_desc,pot,& 00576 deriv) 00577 logger => cp_error_get_logger(error) 00578 00579 IF (.NOT. failure) THEN 00580 bo = rho_r(1)%pw%pw_grid%bounds_local 00581 00582 CALL xc_rho_set_and_dset_create(rho_r=rho_r, rho_g=rho_g,& 00583 tau=tau,xc_section=xc_section,& 00584 cell=cell, pw_pool=pw_pool,rho_set=rho_set1,& 00585 deriv_set=dSet1, deriv_order=1,& 00586 needs_basic_components=.FALSE.,error=error) 00587 00588 ! outputs 0,:,: plane 00589 IF (bo(1,1)<=0.AND.0<=bo(2,1)) THEN 00590 IF (rho_set1%has%rho_spin) THEN 00591 OPEN(unit=120,file="rhoa.bindata",status="unknown",access='sequential',& 00592 form="unformatted",action="write") 00593 pot => rho_set1%rhoa 00594 WRITE(unit=120) pot(0,:,:) 00595 CLOSE(unit=120) 00596 OPEN(unit=120,file="rhob.bindata",status="unknown",access='sequential',& 00597 form="unformatted",action="write") 00598 pot => rho_set1%rhob 00599 WRITE(unit=120) pot(0,:,:) 00600 CLOSE(unit=120) 00601 END IF 00602 IF (rho_set1%has%norm_drho) THEN 00603 OPEN(unit=120,file="ndrho.bindata",status="unknown",access='sequential',& 00604 form="unformatted",action="write") 00605 pot => rho_set1%norm_drho 00606 WRITE(unit=120) pot(0,:,:) 00607 CLOSE(unit=120) 00608 END IF 00609 IF (rho_set1%has%norm_drho_spin) THEN 00610 OPEN(unit=120,file="ndrhoa.bindata",status="unknown",access='sequential',& 00611 form="unformatted",action="write") 00612 pot => rho_set1%norm_drhoa 00613 WRITE(unit=120) pot(0,:,:) 00614 CLOSE(unit=120) 00615 OPEN(unit=120,file="ndrhob.bindata",status="unknown",access='sequential',& 00616 form="unformatted",action="write") 00617 pot => rho_set1%norm_drhob 00618 WRITE(unit=120) pot(0,:,:) 00619 CLOSE(unit=120) 00620 END IF 00621 IF (rho_set1%has%rho) THEN 00622 OPEN(unit=120,file="rho.bindata",status="unknown",access='sequential',& 00623 form="unformatted",action="write") 00624 pot => rho_set1%rho 00625 WRITE(unit=120) pot(0,:,:) 00626 CLOSE(unit=120) 00627 END IF 00628 IF (rho_set1%has%tau) THEN 00629 OPEN(unit=120,file="tau.bindata",status="unknown",access='sequential',& 00630 form="unformatted",action="write") 00631 pot => rho_set1%tau 00632 WRITE(unit=120) pot(0,:,:) 00633 CLOSE(unit=120) 00634 END IF 00635 IF (rho_set1%has%tau_spin) THEN 00636 OPEN(unit=120,file="tau_a.bindata",status="unknown",access='sequential',& 00637 form="unformatted",action="write") 00638 pot => rho_set1%tau_a 00639 WRITE(unit=120) pot(0,:,:) 00640 CLOSE(unit=120) 00641 OPEN(unit=120,file="tau_b.bindata",status="unknown",access='sequential',& 00642 form="unformatted",action="write") 00643 pot => rho_set1%tau_b 00644 WRITE(unit=120) pot(0,:,:) 00645 CLOSE(unit=120) 00646 END IF 00647 00648 deriv_iter => dSet1%derivs 00649 DO WHILE (cp_sll_xc_deriv_next(deriv_iter,el_att=deriv,error=error)) 00650 CALL xc_derivative_get(deriv,& 00651 split_desc=split_desc,deriv_data=pot,& 00652 error=error) 00653 SELECT CASE (SIZE(split_desc)) 00654 CASE(0) 00655 filename="e_0.bindata" 00656 CASE(1) 00657 filename="e_"//TRIM(split_desc(1))//".bindata" 00658 CASE default 00659 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00660 END SELECT 00661 OPEN (unit=120,file=TRIM(filename),status="unknown",& 00662 access='sequential',& 00663 form="unformatted") 00664 WRITE (unit=120) pot(0,:,:) 00665 CLOSE (unit=120) 00666 END DO 00667 END IF 00668 00669 CALL xc_vxc_pw_create(vxc_rho=vxc_rho,vxc_tau=vxc_tau,& 00670 rho_r=rho_r, rho_g=rho_g,tau=tau,& 00671 exc=exc, xc_section=xc_section,& 00672 cell=cell, pw_pool=pw_pool,& 00673 error=error) 00674 00675 ! outputs 0,:,: plane 00676 IF (bo(1,1)<=0.AND.0<=bo(2,1)) THEN 00677 IF (ASSOCIATED(vxc_rho)) THEN 00678 DO ispin=1,SIZE(vxc_rho) 00679 WRITE (filename,"('vxc-',i1,'.bindata')") ispin 00680 OPEN(unit=120,file=filename,status="unknown",access='sequential',& 00681 form="unformatted",action="write") 00682 pot => vxc_rho(ispin)%pw%cr3d 00683 WRITE(unit=120) pot(0,:,:) 00684 CLOSE(unit=120) 00685 00686 pot => vxc_rho(ispin)%pw%cr3d 00687 DO k=bo(1,3),bo(2,3) 00688 DO j=bo(1,2),bo(2,2) 00689 DO i=bo(1,1),bo(2,1) 00690 IF (ABS(pot(i,j,k))>10.0_dp) THEN 00691 WRITE(cp_logger_get_default_unit_nr(logger),& 00692 "(' vxc_',i1,'(',i6,',',i6,',',i6,')=',e11.4)",& 00693 advance="no") ispin,i,j,k,pot(i,j,k) 00694 IF (rho_set1%has%rho_spin) THEN 00695 WRITE(cp_logger_get_default_unit_nr(logger),& 00696 "(' rho=(',e11.4,',',e11.4,')')",advance="no")& 00697 rho_set1%rhoa(i,j,k), rho_set1%rhob(i,j,k) 00698 ELSE IF (rho_set1%has%rho) THEN 00699 WRITE(cp_logger_get_default_unit_nr(logger),& 00700 "(' rho=',e11.4)",advance="no") rho_set1%rho(i,j,k) 00701 END IF 00702 IF (rho_set1%has%norm_drho_spin) THEN 00703 WRITE(cp_logger_get_default_unit_nr(logger),& 00704 "(' ndrho=(',e11.4,',',e11.4,')')",advance="no")& 00705 rho_set1%norm_drhoa(i,j,k), & 00706 rho_set1%norm_drhob(i,j,k) 00707 ELSE IF (rho_set1%has%norm_drho) THEN 00708 WRITE(cp_logger_get_default_unit_nr(logger),& 00709 "(' ndrho=',e11.4)",advance="no") rho_set1%norm_drho(i,j,k) 00710 END IF 00711 IF (rho_set1%has%tau_spin) THEN 00712 WRITE(cp_logger_get_default_unit_nr(logger),& 00713 "(' tau=(',e11.4,',',e11.4,')')",advance="no")& 00714 rho_set1%tau_a(i,j,k), & 00715 rho_set1%tau_b(i,j,k) 00716 ELSE IF (rho_set1%has%tau) THEN 00717 WRITE(cp_logger_get_default_unit_nr(logger),& 00718 "(' tau=',e11.4)",advance="no") rho_set1%tau(i,j,k) 00719 END IF 00720 00721 deriv_iter => dSet1%derivs 00722 DO WHILE (cp_sll_xc_deriv_next(deriv_iter,el_att=deriv,error=error)) 00723 CALL xc_derivative_get(deriv,& 00724 split_desc=split_desc,deriv_data=pot,& 00725 error=error) 00726 SELECT CASE (SIZE(split_desc)) 00727 CASE(0) 00728 filename=" e_0" 00729 CASE(1) 00730 filename=" e_"//TRIM(split_desc(1)) 00731 CASE default 00732 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00733 END SELECT 00734 WRITE (unit=cp_logger_get_default_unit_nr(logger),& 00735 fmt="(a,'=',e11.4)",advance="no") & 00736 TRIM(filename),pot(i,j,k) 00737 END DO 00738 00739 WRITE(cp_logger_get_default_unit_nr(logger),& 00740 "()") 00741 END IF 00742 END DO 00743 END DO 00744 END DO 00745 END DO 00746 END IF 00747 IF (ASSOCIATED(vxc_tau)) THEN 00748 DO ispin=1,SIZE(vxc_tau) 00749 WRITE (filename,"('vxc_tau_',i1,'.bindata')") ispin 00750 OPEN(unit=120,file=filename,status="unknown",access='sequential',& 00751 form="unformatted",action="write") 00752 pot => vxc_tau(ispin)%pw%cr3d 00753 WRITE(unit=120) pot(0,:,:) 00754 CLOSE(unit=120) 00755 00756 pot => vxc_tau(ispin)%pw%cr3d 00757 DO k=bo(1,3),bo(2,3) 00758 DO j=bo(1,2),bo(2,2) 00759 DO i=bo(1,1),bo(2,1) 00760 IF (ABS(pot(i,j,k))>10.0_dp) THEN 00761 WRITE(cp_logger_get_default_unit_nr(logger),& 00762 "('vxc_tau_',i1,'(',i6,',',i6,',',i6,')=',e11.4)",& 00763 advance="no") ispin,i,j,k,pot(i,j,k) 00764 IF (rho_set1%has%rho_spin) THEN 00765 WRITE(cp_logger_get_default_unit_nr(logger),& 00766 "(' rho=(',e11.4,',',e11.4,')')",advance="no")& 00767 rho_set1%rhoa(i,j,k), rho_set1%rhob(i,j,k) 00768 ELSE IF (rho_set1%has%rho) THEN 00769 WRITE(cp_logger_get_default_unit_nr(logger),& 00770 "(' rho=',e11.4)",advance="no") rho_set1%rho(i,j,k) 00771 END IF 00772 IF (rho_set1%has%norm_drho_spin) THEN 00773 WRITE(cp_logger_get_default_unit_nr(logger),& 00774 "(' ndrho=(',e11.4,',',e11.4,')')",advance="no")& 00775 rho_set1%norm_drhoa(i,j,k), & 00776 rho_set1%norm_drhob(i,j,k) 00777 ELSE IF (rho_set1%has%norm_drho) THEN 00778 WRITE(cp_logger_get_default_unit_nr(logger),& 00779 "(' ndrho=',e11.4)",advance="no") rho_set1%norm_drho(i,j,k) 00780 END IF 00781 IF (rho_set1%has%tau_spin) THEN 00782 WRITE(cp_logger_get_default_unit_nr(logger),& 00783 "(' tau=(',e11.4,',',e11.4,')')",advance="no")& 00784 rho_set1%tau_a(i,j,k), & 00785 rho_set1%tau_b(i,j,k) 00786 ELSE IF (rho_set1%has%tau) THEN 00787 WRITE(cp_logger_get_default_unit_nr(logger),& 00788 "(' tau=',e11.4)",advance="no") rho_set1%tau(i,j,k) 00789 END IF 00790 00791 deriv_iter => dSet1%derivs 00792 DO WHILE (cp_sll_xc_deriv_next(deriv_iter,el_att=deriv,error=error)) 00793 CALL xc_derivative_get(deriv,& 00794 split_desc=split_desc,deriv_data=pot,& 00795 error=error) 00796 SELECT CASE (SIZE(split_desc)) 00797 CASE(0) 00798 filename=" e_0" 00799 CASE(1) 00800 filename=" e_"//TRIM(split_desc(1)) 00801 CASE default 00802 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00803 END SELECT 00804 WRITE (unit=cp_logger_get_default_unit_nr(logger),& 00805 fmt="(a,'=',e11.4)",advance="no") & 00806 TRIM(filename),pot(i,j,k) 00807 END DO 00808 00809 WRITE(cp_logger_get_default_unit_nr(logger),"()") 00810 END IF 00811 END DO 00812 END DO 00813 END DO 00814 END DO 00815 END IF 00816 END IF 00817 00818 CALL xc_dset_release(dSet1, error=error) 00819 CALL xc_rho_set_release(rho_set1,error=error) 00820 END IF 00821 END SUBROUTINE xc_vxc_pw_create_debug 00822 00823 ! ***************************************************************************** 00852 SUBROUTINE xc_rho_set_and_dset_create(rho_set,deriv_set,deriv_order,& 00853 rho_r,rho_g,tau,xc_section,cell,pw_pool,& 00854 needs_basic_components,error) 00855 00856 TYPE(xc_rho_set_type), POINTER :: rho_set 00857 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00858 INTEGER, INTENT(in) :: deriv_order 00859 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r, rho_g, tau 00860 TYPE(section_vals_type), POINTER :: xc_section 00861 TYPE(cell_type), POINTER :: cell 00862 TYPE(pw_pool_type), POINTER :: pw_pool 00863 LOGICAL, INTENT(in) :: needs_basic_components 00864 TYPE(cp_error_type), INTENT(inout) :: error 00865 00866 CHARACTER(len=*), PARAMETER :: routineN = 'xc_rho_set_and_dset_create', 00867 routineP = moduleN//':'//routineN 00868 00869 INTEGER :: handle, nspins 00870 LOGICAL :: failure, lsd 00871 TYPE(section_vals_type), POINTER :: xc_fun_sections 00872 00873 CALL timeset(routineN,handle) 00874 failure=.FALSE. 00875 00876 CPPrecondition(.NOT.ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 00877 CPPrecondition(.NOT.ASSOCIATED(deriv_set),cp_failure_level,routineP,error,failure) 00878 CPPrecondition(ASSOCIATED(pw_pool),cp_failure_level,routineP,error,failure) 00879 00880 IF (.NOT. failure) THEN 00881 nspins=SIZE(rho_r) 00882 lsd=(nspins/=1) 00883 END IF 00884 00885 IF (.NOT.failure) THEN 00886 00887 xc_fun_sections => section_vals_get_subs_vals(xc_section,"XC_FUNCTIONAL",& 00888 error=error) 00889 00890 CALL xc_dset_create(deriv_set, pw_pool, error=error) 00891 00892 CALL xc_rho_set_create(rho_set,& 00893 rho_r(1)%pw%pw_grid%bounds_local,& 00894 rho_cutoff=section_get_rval(xc_section,"density_cutoff",error),& 00895 drho_cutoff=section_get_rval(xc_section,"gradient_cutoff",error),& 00896 tau_cutoff=section_get_rval(xc_section,"tau_cutoff",error),& 00897 error=error) 00898 00899 CALL xc_rho_set_update(rho_set, rho_r, rho_g, tau, & 00900 xc_functionals_get_needs(xc_fun_sections,lsd,needs_basic_components,error),& 00901 section_get_ival(xc_section,"XC_GRID%XC_DERIV",error),& 00902 section_get_ival(xc_section,"XC_GRID%XC_SMOOTH_RHO",error),& 00903 cell,pw_pool, error=error) 00904 00905 CALL xc_functionals_eval(xc_fun_sections, & 00906 lsd=lsd,& 00907 rho_set=rho_set, & 00908 deriv_set=deriv_set,& 00909 deriv_order=deriv_order, & 00910 error=error) 00911 00912 END IF 00913 00914 CALL timestop(handle) 00915 00916 END SUBROUTINE xc_rho_set_and_dset_create 00917 00918 ! ***************************************************************************** 00933 SUBROUTINE smooth_cutoff(pot,rho,rhoa,rhob,rho_cutoff,& 00934 rho_smooth_cutoff_range, e_0, e_0_scale_factor,error) 00935 REAL(kind=dp), DIMENSION(:, :, :), 00936 POINTER :: pot, rho, rhoa, rhob 00937 REAL(kind=dp), INTENT(in) :: rho_cutoff, 00938 rho_smooth_cutoff_range 00939 REAL(kind=dp), DIMENSION(:, :, :), 00940 OPTIONAL, POINTER :: e_0 00941 REAL(kind=dp), INTENT(in), OPTIONAL :: e_0_scale_factor 00942 TYPE(cp_error_type), INTENT(inout) :: error 00943 00944 CHARACTER(len=*), PARAMETER :: routineN = 'smooth_cutoff', 00945 routineP = moduleN//':'//routineN 00946 00947 INTEGER :: i, j, k 00948 INTEGER, DIMENSION(2, 3) :: bo 00949 LOGICAL :: failure 00950 REAL(kind=dp) :: my_e_0_scale_factor, my_rho, my_rho_n, my_rho_n2, 00951 rho_smooth_cutoff, rho_smooth_cutoff_2, rho_smooth_cutoff_range_2 00952 00953 failure=.FALSE. 00954 CPPrecondition(ASSOCIATED(pot),cp_failure_level,routineP,error,failure) 00955 bo(1,:)=LBOUND(pot) 00956 bo(2,:)=UBOUND(pot) 00957 my_e_0_scale_factor=1.0_dp 00958 IF (PRESENT(e_0_scale_factor)) my_e_0_scale_factor=e_0_scale_factor 00959 IF (.NOT. failure) THEN 00960 rho_smooth_cutoff=rho_cutoff*rho_smooth_cutoff_range 00961 rho_smooth_cutoff_2=(rho_cutoff+rho_smooth_cutoff)/2 00962 rho_smooth_cutoff_range_2=rho_smooth_cutoff_2-rho_cutoff 00963 00964 IF (rho_smooth_cutoff_range>0.0_dp) THEN 00965 IF (PRESENT(e_0)) THEN 00966 CPPrecondition(ASSOCIATED(e_0),cp_failure_level,routineP,error,failure) 00967 IF (ASSOCIATED(rho)) THEN 00968 !$omp parallel do default(none) shared(bo,e_0,pot,rho,& 00969 !$omp rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2,& 00970 !$omp rho_smooth_cutoff_range_2,my_e_0_scale_factor)& 00971 !$omp private(k,j,i,my_rho,my_rho_n,my_rho_n2) 00972 DO k = bo(1,3), bo(2,3) 00973 DO j = bo(1,2), bo(2,2) 00974 DO i = bo(1,1), bo(2,1) 00975 my_rho=rho(i,j,k) 00976 IF (my_rho<rho_smooth_cutoff) THEN 00977 IF (my_rho<rho_cutoff) THEN 00978 pot(i,j,k)=0.0_dp 00979 ELSEIF (my_rho<rho_smooth_cutoff_2) THEN 00980 my_rho_n=(my_rho-rho_cutoff)/rho_smooth_cutoff_range_2 00981 my_rho_n2=my_rho_n*my_rho_n 00982 pot(i,j,k)=pot(i,j,k)*& 00983 my_rho_n2*(my_rho_n-0.5_dp*my_rho_n2)+& 00984 my_e_0_scale_factor*e_0(i,j,k)*& 00985 my_rho_n2*(3.0_dp-2.0_dp*my_rho_n)& 00986 /rho_smooth_cutoff_range_2 00987 ELSE 00988 my_rho_n=2.0_dp-(my_rho-rho_cutoff)/rho_smooth_cutoff_range_2 00989 my_rho_n2=my_rho_n*my_rho_n 00990 pot(i,j,k)=pot(i,j,k)*& 00991 (1.0_dp-my_rho_n2*(my_rho_n-0.5_dp*my_rho_n2))& 00992 +my_e_0_scale_factor*e_0(i,j,k)*& 00993 my_rho_n2*(3.0_dp-2.0_dp*my_rho_n)& 00994 /rho_smooth_cutoff_range_2 00995 END IF 00996 END IF 00997 END DO 00998 END DO 00999 END DO 01000 ELSE 01001 !$omp parallel do default(none) shared(bo,pot,e_0,rhoa,rhob,& 01002 !$omp rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2,& 01003 !$omp rho_smooth_cutoff_range_2,my_e_0_scale_factor)& 01004 !$omp private(k,j,i,my_rho,my_rho_n,my_rho_n2) 01005 DO k = bo(1,3), bo(2,3) 01006 DO j = bo(1,2), bo(2,2) 01007 DO i = bo(1,1), bo(2,1) 01008 my_rho=rhoa(i,j,k)+rhob(i,j,k) 01009 IF (my_rho<rho_smooth_cutoff) THEN 01010 IF (my_rho<rho_cutoff) THEN 01011 pot(i,j,k)=0.0_dp 01012 ELSEIF (my_rho<rho_smooth_cutoff_2) THEN 01013 my_rho_n=(my_rho-rho_cutoff)/rho_smooth_cutoff_range_2 01014 my_rho_n2=my_rho_n*my_rho_n 01015 pot(i,j,k)=pot(i,j,k)*& 01016 my_rho_n2*(my_rho_n-0.5_dp*my_rho_n2)+& 01017 my_e_0_scale_factor*e_0(i,j,k)*& 01018 my_rho_n2*(3.0_dp-2.0_dp*my_rho_n)& 01019 /rho_smooth_cutoff_range_2 01020 ELSE 01021 my_rho_n=2.0_dp-(my_rho-rho_cutoff)/rho_smooth_cutoff_range_2 01022 my_rho_n2=my_rho_n*my_rho_n 01023 pot(i,j,k)=pot(i,j,k)*& 01024 (1.0_dp-my_rho_n2*(my_rho_n-0.5_dp*my_rho_n2))& 01025 +my_e_0_scale_factor*e_0(i,j,k)*& 01026 my_rho_n2*(3.0_dp-2.0_dp*my_rho_n)& 01027 /rho_smooth_cutoff_range_2 01028 END IF 01029 END IF 01030 END DO 01031 END DO 01032 END DO 01033 END IF 01034 ELSE 01035 IF (ASSOCIATED(rho)) THEN 01036 !$omp parallel do default(none) shared(bo,pot,& 01037 !$omp rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2,& 01038 !$omp rho_smooth_cutoff_range_2,rho)& 01039 !$omp private(k,j,i,my_rho,my_rho_n,my_rho_n2) 01040 DO k = bo(1,3), bo(2,3) 01041 DO j = bo(1,2), bo(2,2) 01042 DO i = bo(1,1), bo(2,1) 01043 my_rho=rho(i,j,k) 01044 IF (my_rho<rho_smooth_cutoff) THEN 01045 IF (my_rho<rho_cutoff) THEN 01046 pot(i,j,k)=0.0_dp 01047 ELSEIF (my_rho<rho_smooth_cutoff_2) THEN 01048 my_rho_n=(my_rho-rho_cutoff)/rho_smooth_cutoff_range_2 01049 my_rho_n2=my_rho_n*my_rho_n 01050 pot(i,j,k)=pot(i,j,k)*& 01051 my_rho_n2*(my_rho_n-0.5_dp*my_rho_n2) 01052 ELSE 01053 my_rho_n=2.0_dp-(my_rho-rho_cutoff)/rho_smooth_cutoff_range_2 01054 my_rho_n2=my_rho_n*my_rho_n 01055 pot(i,j,k)=pot(i,j,k)*& 01056 (1.0_dp-my_rho_n2*(my_rho_n-0.5_dp*my_rho_n2)) 01057 END IF 01058 END IF 01059 END DO 01060 END DO 01061 END DO 01062 ELSE 01063 CPPrecondition(ASSOCIATED(rhoa),cp_failure_level,routineP,error,failure) 01064 CPPrecondition(ASSOCIATED(rhob),cp_failure_level,routineP,error,failure) 01065 !$omp parallel do default(none) shared(bo,pot,& 01066 !$omp rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2,& 01067 !$omp rho_smooth_cutoff_range_2,rhoa,rhob)& 01068 !$omp private(k,j,i,my_rho,my_rho_n,my_rho_n2) 01069 DO k = bo(1,3), bo(2,3) 01070 DO j = bo(1,2), bo(2,2) 01071 DO i = bo(1,1), bo(2,1) 01072 my_rho=rhoa(i,j,k)+rhob(i,j,k) 01073 IF (my_rho<rho_smooth_cutoff) THEN 01074 IF (my_rho<rho_cutoff) THEN 01075 pot(i,j,k)=0.0_dp 01076 ELSEIF (my_rho<rho_smooth_cutoff_2) THEN 01077 my_rho_n=(my_rho-rho_cutoff)/rho_smooth_cutoff_range_2 01078 my_rho_n2=my_rho_n*my_rho_n 01079 pot(i,j,k)=pot(i,j,k)*& 01080 my_rho_n2*(my_rho_n-0.5_dp*my_rho_n2) 01081 ELSE 01082 my_rho_n=2.0_dp-(my_rho-rho_cutoff)/rho_smooth_cutoff_range_2 01083 my_rho_n2=my_rho_n*my_rho_n 01084 pot(i,j,k)=pot(i,j,k)*& 01085 (1.0_dp-my_rho_n2*(my_rho_n-0.5_dp*my_rho_n2)) 01086 END IF 01087 END IF 01088 END DO 01089 END DO 01090 END DO 01091 END IF 01092 END IF 01093 END IF 01094 END IF 01095 END SUBROUTINE smooth_cutoff 01096 01097 ! ***************************************************************************** 01128 SUBROUTINE xc_vxc_pw_create(vxc_rho,vxc_tau,exc,rho_r,rho_g,tau,xc_section,& 01129 cell,pw_pool,error,virial) 01130 TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau 01131 REAL(KIND=dp), INTENT(out) :: exc 01132 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r, rho_g, tau 01133 TYPE(section_vals_type), POINTER :: xc_section 01134 TYPE(cell_type), POINTER :: cell 01135 TYPE(pw_pool_type), POINTER :: pw_pool 01136 TYPE(cp_error_type), INTENT(inout) :: error 01137 TYPE(virial_type), OPTIONAL, POINTER :: virial 01138 01139 CHARACTER(len=*), PARAMETER :: routineN = 'xc_vxc_pw_create', 01140 routineP = moduleN//':'//routineN 01141 CHARACTER(len=30), DIMENSION(2), PARAMETER :: 01142 norm_drho_spin_name = (/ "(norm_drhoa)", "(norm_drhob)" /) 01143 01144 CHARACTER& 01145 (len=MAX_DERIVATIVE_DESC_LENGTH) :: desc 01146 INTEGER :: handle, i, idir, ispin, j, jdir, k, n_deriv, npoints, nspins, 01147 order, stat, xc_deriv_method_id, xc_rho_smooth_id 01148 INTEGER, DIMENSION(2, 3) :: bo 01149 INTEGER, DIMENSION(3, 3) :: nd, nd_laplace 01150 LOGICAL :: dealloc_pw_to_deriv, failure, 01151 has_laplace, has_tau, lsd, 01152 use_virial, zero_result 01153 REAL(KIND=dp) :: density_smooth_cut_range, 01154 drho_cutoff, my_rho, ndr, 01155 rho_cutoff 01156 REAL(kind=dp), DIMENSION(:, :, :), 01157 POINTER :: deriv_data, norm_drho, 01158 norm_drho_spin, rho, rhoa, 01159 rhob, tmp_cr3d 01160 TYPE(cp_3d_r_p_type), DIMENSION(:), 01161 POINTER :: drho, drho_spin, drhoa, drhob 01162 TYPE(cp_sll_xc_deriv_type), POINTER :: pos 01163 TYPE(pw_grid_type), POINTER :: pw_grid 01164 TYPE(pw_p_type), DIMENSION(2) :: vxc_to_deriv 01165 TYPE(pw_p_type), DIMENSION(3) :: pw_to_deriv, pw_to_deriv_rho 01166 TYPE(pw_type), POINTER :: tmp_g, tmp_r, virial_pw, vxc_g 01167 TYPE(xc_derivative_set_type), POINTER :: deriv_set 01168 TYPE(xc_derivative_type), POINTER :: deriv_att 01169 TYPE(xc_rho_set_type), POINTER :: rho_set 01170 01171 CALL timeset(routineN,handle) 01172 failure=.FALSE. 01173 NULLIFY(tmp_g, tmp_r, vxc_g, norm_drho_spin, norm_drho, drho_spin, drhoa, & 01174 drhob, pos, deriv_set, rho_set, virial_pw) 01175 nd = RESHAPE ((/1,0,0,0,1,0,0,0,1/),(/3,3/)) 01176 DO idir=1,3 01177 NULLIFY(pw_to_deriv(idir)%pw, pw_to_deriv_rho(idir)%pw) 01178 END DO 01179 DO i=1,2 01180 NULLIFY(vxc_to_deriv(i)%pw) 01181 END DO 01182 01183 pw_grid => rho_r(1)%pw%pw_grid 01184 01185 CPPrecondition(ASSOCIATED(xc_section),cp_failure_level,routineP,error,failure) 01186 CPPrecondition(ASSOCIATED(pw_pool),cp_failure_level,routineP,error,failure) 01187 CPPrecondition(.NOT.ASSOCIATED(vxc_rho),cp_failure_level,routineP,error,failure) 01188 CPPrecondition(.NOT.ASSOCIATED(vxc_tau),cp_failure_level,routineP,error,failure) 01189 IF (.NOT. failure) THEN 01190 nspins=SIZE(rho_r) 01191 CPPrecondition(ASSOCIATED(rho_r(SIZE(rho_r))%pw),cp_failure_level,routineP,error,failure) 01192 lsd=(nspins/=1) 01193 IF (lsd) THEN 01194 CPPrecondition(nspins==2,cp_failure_level,routineP,error,failure) 01195 END IF 01196 END IF 01197 01198 IF (PRESENT(virial)) THEN 01199 use_virial = virial%pv_calculate.AND.(.NOT.virial%pv_numer) 01200 ELSE 01201 use_virial = .FALSE. 01202 ENDIF 01203 01204 IF (.NOT.failure) THEN 01205 bo = rho_r(1)%pw%pw_grid%bounds_local 01206 npoints=(bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1) 01207 01208 ! calculate the potential derivatives 01209 CALL xc_rho_set_and_dset_create(rho_set=rho_set,deriv_set=deriv_set,& 01210 deriv_order=1, rho_r=rho_r, rho_g=rho_g, tau=tau,& 01211 xc_section=xc_section,& 01212 cell=cell, pw_pool=pw_pool,& 01213 needs_basic_components=.TRUE.,& 01214 error=error) 01215 01216 END IF 01217 01218 IF (.NOT.failure) THEN 01219 CALL section_vals_val_get(xc_section,"XC_GRID%XC_DERIV",& 01220 i_val=xc_deriv_method_id,error=error) 01221 CALL section_vals_val_get(xc_section,"XC_GRID%XC_SMOOTH_RHO",& 01222 i_val=xc_rho_smooth_id,error=error) 01223 CALL section_vals_val_get(xc_section,"DENSITY_SMOOTH_CUTOFF_RANGE",& 01224 r_val=density_smooth_cut_range,error=error) 01225 01226 CALL xc_rho_set_get(rho_set,rho_cutoff=rho_cutoff,& 01227 drho_cutoff=drho_cutoff, error=error) 01228 01229 has_tau=.FALSE. 01230 has_laplace = .FALSE. 01231 ! check for unknown derivatives 01232 n_deriv=0 01233 pos => deriv_set%derivs 01234 DO WHILE (cp_sll_xc_deriv_next(pos,el_att=deriv_att,error=error)) 01235 CALL xc_derivative_get(deriv_att,order=order,& 01236 desc=desc,& 01237 error=error) 01238 IF (order==1) THEN 01239 IF (lsd) THEN 01240 SELECT CASE(desc) 01241 CASE("(rho)","(rhoa)","(rhob)","(norm_drho)","(norm_drhoa)",& 01242 "(norm_drhob)") 01243 n_deriv=n_deriv+1 01244 CASE("(tau)","(tau_a)","(tau_b)") 01245 has_tau=.TRUE. 01246 n_deriv=n_deriv+1 01247 CASE("(laplace_rhoa)","(laplace_rhob)") 01248 has_laplace = .TRUE. 01249 n_deriv = n_deriv + 1 01250 CASE default 01251 !FM if you are looking at this error probably you are missing the 01252 !FM cross term (drhoa_drhob), I never got round to implement it, 01253 !FM in the functionals that I have implemented I use norm_drho 01254 !FM instead, either do the same or implement it, it shouldn't be 01255 !FM difficult (I am just lazy ;) 01256 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 01257 "unknown functional derivative (LSD): '"//& 01258 TRIM(desc)//"' in "//& 01259 CPSourceFileRef,& 01260 error,failure) 01261 END SELECT 01262 ELSE 01263 SELECT CASE(desc) 01264 CASE("(rho)","(norm_drho)") 01265 n_deriv=n_deriv+1 01266 CASE("(tau)") 01267 has_tau=.TRUE. 01268 n_deriv=n_deriv+1 01269 CASE("(laplace_rho)") 01270 has_laplace = .TRUE. 01271 n_deriv = n_deriv + 1 01272 CASE default 01273 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 01274 "unknown functional derivative (LDA): '"//& 01275 TRIM(desc)//"' in "//& 01276 CPSourceFileRef,& 01277 error,failure) 01278 END SELECT 01279 END IF 01280 END IF 01281 END DO 01282 END IF 01283 01284 IF (.NOT.failure) THEN 01285 ALLOCATE(vxc_rho(nspins),stat=stat) 01286 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01287 DO ispin=1,nspins 01288 NULLIFY(vxc_rho(ispin)%pw) 01289 END DO 01290 01291 CALL xc_rho_set_get(rho_set,rho=rho,rhoa=rhoa,rhob=rhob,& 01292 can_return_null=.TRUE., error=error) 01293 01294 ! recover the vxc arrays 01295 IF (lsd) THEN 01296 deriv_att => xc_dset_get_derivative(deriv_set, "(rhoa)",error=error) 01297 IF (ASSOCIATED(deriv_att)) THEN 01298 CALL pw_create(vxc_rho(1)%pw,pw_grid=pw_grid,& 01299 cr3d_ptr=deriv_att%deriv_data,& 01300 use_data=REALDATA3D,in_space=REALSPACE,& 01301 error=error) 01302 NULLIFY(deriv_att%deriv_data) 01303 ELSE 01304 CALL pw_pool_create_pw(pw_pool,vxc_rho(1)%pw,& 01305 use_data=REALDATA3D,in_space=REALSPACE,error=error) 01306 CALL pw_zero(vxc_rho(1)%pw, error=error) 01307 END IF 01308 01309 deriv_att => xc_dset_get_derivative(deriv_set, "(rhob)",error=error) 01310 IF (ASSOCIATED(deriv_att)) THEN 01311 CALL pw_create(vxc_rho(2)%pw,pw_grid=pw_grid,& 01312 cr3d_ptr=deriv_att%deriv_data,& 01313 use_data=REALDATA3D,in_space=REALSPACE,& 01314 error=error) 01315 NULLIFY(deriv_att%deriv_data) 01316 ELSE 01317 CALL pw_pool_create_pw(pw_pool,vxc_rho(2)%pw,& 01318 use_data=REALDATA3D,in_space=REALSPACE,error=error) 01319 CALL pw_zero(vxc_rho(2)%pw, error=error) 01320 END IF 01321 01322 ELSE 01323 deriv_att => xc_dset_get_derivative(deriv_set, "(rho)",error=error) 01324 IF (ASSOCIATED(deriv_att)) THEN 01325 CALL pw_create(vxc_rho(1)%pw,pw_grid=pw_grid,& 01326 cr3d_ptr=deriv_att%deriv_data,& 01327 use_data=REALDATA3D,in_space=REALSPACE,& 01328 error=error) 01329 NULLIFY(deriv_att%deriv_data) 01330 ELSE 01331 CALL pw_pool_create_pw(pw_pool,vxc_rho(1)%pw,& 01332 use_data=REALDATA3D,in_space=REALSPACE,error=error) 01333 CALL pw_zero(vxc_rho(1)%pw, error=error) 01334 END IF 01335 END IF 01336 01337 deriv_att => xc_dset_get_derivative(deriv_set, "(rho)",error=error) 01338 IF (ASSOCIATED(deriv_att)) THEN 01339 IF (lsd) THEN ! otherwise already taken care in vxc recovery 01340 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,& 01341 error=error) 01342 !$omp parallel do default(none) shared(bo,vxc_rho,deriv_data)& 01343 !$omp private(k,j,i) 01344 DO k = bo(1,3), bo(2,3) 01345 DO j = bo(1,2), bo(2,2) 01346 DO i = bo(1,1), bo(2,1) 01347 vxc_rho(1)%pw%cr3d(i,j,k) = vxc_rho(1)%pw%cr3d(i,j,k)+deriv_data(i,j,k) 01348 vxc_rho(2)%pw%cr3d(i,j,k) = vxc_rho(2)%pw%cr3d(i,j,k)+deriv_data(i,j,k) 01349 END DO 01350 END DO 01351 END DO 01352 CALL pw_pool_give_back_cr3d(pw_pool,deriv_att%deriv_data,& 01353 error=error) 01354 NULLIFY(deriv_att%deriv_data) 01355 END IF 01356 END IF 01357 01358 ! rhoa,rhob already taken care of in vxc recovery 01359 01360 deriv_att => xc_dset_get_derivative(deriv_set, "(norm_drho)",error=error) 01361 IF (ASSOCIATED(deriv_att)) THEN 01362 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,& 01363 error=error) 01364 01365 CALL xc_rho_set_get(rho_set,norm_drho=norm_drho,& 01366 drho=drho,drhoa=drhoa,drhob=drhob,rho_cutoff=rho_cutoff,& 01367 drho_cutoff=drho_cutoff,& 01368 can_return_null=.TRUE., error=error) 01369 01370 IF (ASSOCIATED(norm_drho)) THEN 01371 !$omp parallel do default(none) shared(bo,deriv_data,norm_drho,drho_cutoff)& 01372 !$omp private(k,j,i) 01373 DO k = bo(1,3), bo(2,3) 01374 DO j = bo(1,2), bo(2,2) 01375 DO i = bo(1,1), bo(2,1) 01376 deriv_data(i,j,k) = -deriv_data(i,j,k)/MAX(norm_drho(i,j,k),drho_cutoff) 01377 END DO 01378 END DO 01379 END DO 01380 ELSE IF (ASSOCIATED(drho)) THEN 01381 !$omp parallel do default(none) shared(bo,deriv_data,drho,drho_cutoff)& 01382 !$omp private(k,j,i,ndr) 01383 DO k = bo(1,3), bo(2,3) 01384 DO j = bo(1,2), bo(2,2) 01385 DO i = bo(1,1), bo(2,1) 01386 ndr = SQRT(drho(1)%array(i,j,k)**2 +& 01387 drho(2)%array(i,j,k)**2 +& 01388 drho(3)%array(i,j,k)**2) 01389 deriv_data(i,j,k) = -deriv_data(i,j,k)/MAX(ndr,drho_cutoff) 01390 END DO 01391 END DO 01392 END DO 01393 ELSE 01394 CPPrecondition(ASSOCIATED(drhoa),cp_failure_level,routineP,error,failure) 01395 CPPrecondition(ASSOCIATED(drhob),cp_failure_level,routineP,error,failure) 01396 !$omp parallel do default(none) shared(bo,deriv_data,drhoa,drhob,drho_cutoff)& 01397 !$omp private(k,j,i,ndr) 01398 DO k = bo(1,3), bo(2,3) 01399 DO j = bo(1,2), bo(2,2) 01400 DO i = bo(1,1), bo(2,1) 01401 ndr = SQRT((drhoa(1)%array(i,j,k) + drhob(1)%array(i,j,k))**2 +& 01402 (drhoa(2)%array(i,j,k) + drhob(2)%array(i,j,k))**2 +& 01403 (drhoa(3)%array(i,j,k) + drhob(3)%array(i,j,k))**2) 01404 deriv_data(i,j,k) = -deriv_data(i,j,k)/MAX(ndr,drho_cutoff) 01405 END DO 01406 END DO 01407 END DO 01408 END IF 01409 01410 IF (ASSOCIATED(drho).AND.ASSOCIATED(drho(1)%array)) THEN 01411 CPPrecondition(ASSOCIATED(deriv_data),cp_failure_level,routineP,error,failure) 01412 IF (use_virial) THEN 01413 CALL pw_pool_create_pw(pw_pool,virial_pw,& 01414 use_data=REALDATA3D,& 01415 in_space=REALSPACE,& 01416 error=error) 01417 CALL pw_zero(virial_pw, error=error) 01418 DO idir=1,3 01419 DO k=bo(1,3),bo(2,3) 01420 DO j=bo(1,2),bo(2,2) 01421 DO i=bo(1,1),bo(2,1) 01422 virial_pw%cr3d(i,j,k) = drho(idir)%array(i,j,k)*deriv_data(i,j,k) 01423 END DO 01424 END DO 01425 END DO 01426 DO jdir=1,idir 01427 virial%pv_xc(idir,jdir) = pw_grid%dvol*& 01428 accurate_sum(virial_pw%cr3d(:,:,:)*& 01429 drho(jdir)%array(:,:,:)) 01430 virial%pv_xc(jdir,idir) = virial%pv_xc(idir,jdir) 01431 END DO 01432 END DO 01433 CALL pw_pool_give_back_pw(pw_pool,virial_pw,error=error) 01434 END IF ! use_virial 01435 DO idir=1,3 01436 CALL pw_create(pw_to_deriv_rho(idir)%pw,pw_grid=pw_grid,& 01437 cr3d_ptr=drho(idir)%array,& 01438 use_data=REALDATA3D,in_space=REALSPACE,& 01439 error=error) 01440 CPPrecondition(ASSOCIATED(drho(idir)%array),cp_failure_level,routineP,error,failure) 01441 !$omp parallel do default(none) shared(bo,deriv_data,drho,idir,use_virial,virial_pw)& 01442 !$omp private(k,j,i) 01443 DO k=bo(1,3),bo(2,3) 01444 DO j=bo(1,2),bo(2,2) 01445 DO i=bo(1,1),bo(2,1) 01446 drho(idir)%array(i,j,k) = drho(idir)%array(i,j,k)*deriv_data(i,j,k) 01447 END DO 01448 END DO 01449 END DO 01450 NULLIFY (drho(idir)%array) 01451 END DO 01452 ELSE 01453 IF (use_virial) THEN 01454 CALL pw_pool_create_pw(pw_pool,virial_pw,& 01455 use_data=REALDATA3D,& 01456 in_space=REALSPACE,& 01457 error=error) 01458 CALL pw_zero(virial_pw, error=error) 01459 DO idir=1,3 01460 DO k=bo(1,3),bo(2,3) 01461 DO j=bo(1,2),bo(2,2) 01462 DO i=bo(1,1),bo(2,1) 01463 my_rho = drhoa(idir)%array(i,j,k) + drhob(idir)%array(i,j,k) 01464 virial_pw%cr3d(i,j,k) = my_rho*deriv_data(i,j,k) 01465 END DO 01466 END DO 01467 END DO 01468 DO jdir=1,idir 01469 virial%pv_xc(idir,jdir) = pw_grid%dvol*accurate_sum(virial_pw%cr3d(:,:,:)*& 01470 (drhoa(jdir)%array(:,:,:) + drhob(jdir)%array(:,:,:))) 01471 virial%pv_xc(jdir,idir) = virial%pv_xc(idir,jdir) 01472 END DO 01473 END DO 01474 CALL pw_pool_give_back_pw(pw_pool,virial_pw,error=error) 01475 END IF ! use_virial 01476 DO idir=1,3 01477 CALL pw_pool_create_pw(pw_pool,pw_to_deriv_rho(idir)%pw,& 01478 use_data=REALDATA3D,in_space=REALSPACE,& 01479 error=error) 01480 !$omp parallel do default(none) shared(bo,deriv_data,& 01481 !$omp pw_to_deriv_rho,drhoa,drhob,idir)& 01482 !$omp private(k,j,i,my_rho) 01483 DO k=bo(1,3),bo(2,3) 01484 DO j=bo(1,2),bo(2,2) 01485 DO i=bo(1,1),bo(2,1) 01486 my_rho = drhoa(idir)%array(i,j,k) + drhob(idir)%array(i,j,k) 01487 pw_to_deriv_rho(idir)%pw%cr3d(i,j,k) = my_rho*deriv_data(i,j,k) 01488 END DO 01489 END DO 01490 END DO 01491 END DO 01492 END IF 01493 01494 CALL pw_pool_give_back_cr3d(pw_pool, deriv_att%deriv_data, error=error) 01495 01496 END IF 01497 01498 DO ispin=1,nspins 01499 01500 IF (lsd) THEN 01501 IF (ispin==1) THEN 01502 CALL xc_rho_set_get(rho_set,drhoa=drho_spin,& 01503 can_return_null=.TRUE.,error=error) 01504 CALL xc_rho_set_get(rho_set,norm_drhoa=norm_drho_spin,& 01505 can_return_null=.TRUE., error=error) 01506 ELSE 01507 CALL xc_rho_set_get(rho_set,drhob=drho_spin,& 01508 can_return_null=.TRUE.,error=error) 01509 CALL xc_rho_set_get(rho_set,norm_drhob=norm_drho_spin,& 01510 can_return_null=.TRUE., error=error) 01511 END IF 01512 01513 deriv_att => xc_dset_get_derivative(deriv_set, norm_drho_spin_name(ispin),error=error) 01514 IF (ASSOCIATED(deriv_att)) THEN 01515 CPPrecondition(lsd,cp_failure_level,routineP,error,failure) 01516 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,& 01517 error=error) 01518 CALL pw_create(tmp_r,pw_grid,& 01519 use_data=REALDATA3D,in_space=REALSPACE,& 01520 cr3d_ptr=deriv_data, error=error) 01521 01522 IF (ASSOCIATED(norm_drho_spin)) THEN 01523 !$omp parallel do default(none) shared(bo,deriv_data,norm_drho_spin,drho_cutoff)& 01524 !$omp private(k,j,i) 01525 DO k = bo(1,3), bo(2,3) 01526 DO j = bo(1,2), bo(2,2) 01527 DO i = bo(1,1), bo(2,1) 01528 deriv_data(i,j,k) = -deriv_data(i,j,k)/& 01529 MAX(norm_drho_spin(i,j,k),drho_cutoff) 01530 END DO 01531 END DO 01532 END DO 01533 ELSE 01534 !$omp parallel do default(none) shared(bo,deriv_data,drho_spin,drho_cutoff)& 01535 !$omp private(k,j,i, ndr) 01536 DO k = bo(1,3), bo(2,3) 01537 DO j = bo(1,2), bo(2,2) 01538 DO i = bo(1,1), bo(2,1) 01539 ndr = SQRT(drho_spin(1)%array(i,j,k)**2 +& 01540 drho_spin(2)%array(i,j,k)**2 +& 01541 drho_spin(3)%array(i,j,k)**2) 01542 deriv_data(i,j,k) = -deriv_data(i,j,k)/MAX(ndr,drho_cutoff) 01543 END DO 01544 END DO 01545 END DO 01546 END IF 01547 01548 IF (use_virial) THEN 01549 CALL pw_pool_create_pw(pw_pool,virial_pw,& 01550 use_data=REALDATA3D,& 01551 in_space=REALSPACE,& 01552 error=error) 01553 CALL pw_zero(virial_pw, error=error) 01554 DO idir=1,3 01555 DO k=bo(1,3),bo(2,3) 01556 DO j=bo(1,2),bo(2,2) 01557 DO i=bo(1,1),bo(2,1) 01558 virial_pw%cr3d(i,j,k) = drho_spin(idir)%array(i,j,k)*deriv_data(i,j,k) 01559 END DO 01560 END DO 01561 END DO 01562 DO jdir=1,idir 01563 virial%pv_xc(idir,jdir) = virial%pv_xc(idir,jdir) + pw_grid%dvol*& 01564 accurate_sum(virial_pw%cr3d(:,:,:)*& 01565 drho_spin(jdir)%array(:,:,:)) 01566 virial%pv_xc(jdir,idir) = virial%pv_xc(idir,jdir) 01567 END DO 01568 END DO 01569 CALL pw_pool_give_back_pw(pw_pool,virial_pw,error=error) 01570 END IF ! use_virial 01571 01572 vxc_to_deriv(ispin)%pw => tmp_r 01573 NULLIFY(tmp_r,deriv_att%deriv_data) 01574 01575 DO idir=1,3 01576 CPPrecondition(ASSOCIATED(drho_spin(idir)%array),cp_failure_level,routineP,error,failure) 01577 CPPrecondition(ASSOCIATED(vxc_to_deriv(ispin)%pw%cr3d),cp_failure_level,routineP,error,failure) 01578 !$omp parallel do default(none) shared(bo,deriv_data,drho_spin,& 01579 !$omp ispin,idir,vxc_to_deriv,drho_cutoff)& 01580 !$omp private(k,j,i) 01581 DO k=bo(1,3),bo(2,3) 01582 DO j=bo(1,2),bo(2,2) 01583 DO i=bo(1,1),bo(2,1) 01584 drho_spin(idir)%array(i,j,k)= vxc_to_deriv(ispin)%pw%cr3d(i,j,k)*& 01585 drho_spin(idir)%array(i,j,k) 01586 END DO 01587 END DO 01588 END DO 01589 01590 CALL pw_create(pw_to_deriv(idir)%pw,pw_grid=pw_grid,& 01591 cr3d_ptr=drho_spin(idir)%array,& 01592 use_data=REALDATA3D,in_space=REALSPACE,& 01593 error=error) 01594 NULLIFY(drho_spin(idir)%array) 01595 END DO 01596 01597 dealloc_pw_to_deriv=.TRUE. 01598 CALL pw_pool_give_back_pw(pw_pool,vxc_to_deriv(ispin)%pw,error=error) 01599 END IF ! deriv_att 01600 01601 END IF ! LSD 01602 01603 IF (ASSOCIATED(pw_to_deriv_rho(1)%pw)) THEN 01604 IF (.NOT.ASSOCIATED(pw_to_deriv(1)%pw)) THEN 01605 DO idir=1,3 01606 pw_to_deriv(idir)%pw => pw_to_deriv_rho(idir)%pw 01607 END DO 01608 dealloc_pw_to_deriv=(.NOT.lsd).OR.(ispin==2) 01609 IF (dealloc_pw_to_deriv) THEN 01610 DO idir=1,3 01611 NULLIFY(pw_to_deriv_rho(idir)%pw) 01612 END DO 01613 END IF 01614 ELSE 01615 DO idir=1,3 01616 CALL pw_axpy(pw_to_deriv_rho(idir)%pw,pw_to_deriv(idir)%pw, error=error) 01617 IF (ispin==2) THEN 01618 CALL pw_pool_give_back_pw(pw_pool,pw_to_deriv_rho(idir)%pw,& 01619 error=error) 01620 END IF 01621 END DO 01622 END IF 01623 END IF 01624 01625 IF (ASSOCIATED(pw_to_deriv(1)%pw)) THEN 01626 ! partial integration 01627 IF (xc_deriv_method_id/=xc_deriv_pw) THEN 01628 CALL pw_spline_scale_deriv(pw_to_deriv, cell=cell,& 01629 transpose=.TRUE.,& 01630 error=error) 01631 END IF 01632 01633 IF (xc_deriv_method_id==xc_deriv_pw.OR.& 01634 xc_deriv_method_id==xc_deriv_spline3.OR.& 01635 xc_deriv_method_id==xc_deriv_spline2) THEN 01636 01637 IF (.NOT.ASSOCIATED(vxc_g)) THEN 01638 CALL pw_pool_create_pw(pw_pool,vxc_g,& 01639 use_data=COMPLEXDATA1D,in_space=RECIPROCALSPACE,& 01640 error=error) 01641 zero_result=.TRUE. 01642 ELSE 01643 zero_result=.FALSE. 01644 END IF 01645 01646 DO idir = 1,3 01647 IF (zero_result .AND. idir==1) THEN 01648 tmp_g => vxc_g 01649 ELSE 01650 CALL pw_pool_create_pw(pw_pool,tmp_g,& 01651 use_data=COMPLEXDATA1D,in_space=RECIPROCALSPACE,& 01652 error=error) 01653 END IF 01654 01655 CALL pw_transfer ( pw_to_deriv(idir)%pw, tmp_g , error=error) 01656 01657 SELECT CASE(xc_deriv_method_id) 01658 CASE (xc_deriv_pw) 01659 CALL pw_derive ( tmp_g, nd(:,idir) , error=error) 01660 CASE (xc_deriv_spline2) 01661 CALL pw_spline2_interpolate_values_g(tmp_g,error=error) 01662 CALL pw_spline2_deriv_g ( tmp_g, idir=idir, error=error ) 01663 CASE (xc_deriv_spline3) 01664 CALL pw_spline3_interpolate_values_g(tmp_g,error=error) 01665 CALL pw_spline3_deriv_g ( tmp_g, idir=idir, error=error ) 01666 CASE default 01667 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 01668 END SELECT 01669 01670 IF (zero_result .AND. idir==1) THEN 01671 NULLIFY(tmp_g) 01672 ELSE 01673 CALL pw_axpy ( tmp_g, vxc_g , error=error) 01674 CALL pw_pool_give_back_pw(pw_pool,tmp_g,error=error) 01675 END IF 01676 IF(dealloc_pw_to_deriv) THEN 01677 CALL pw_pool_give_back_pw(pw_pool,pw_to_deriv(idir)%pw,error=error) 01678 END IF 01679 END DO 01680 ! transfer vxc in real space 01681 CALL pw_pool_create_pw(pw_pool, tmp_r,& 01682 use_data=REALDATA3D, in_space=REALSPACE,& 01683 error=error) 01684 CALL pw_transfer ( vxc_g, tmp_r , error=error) 01685 CALL pw_axpy ( tmp_r, vxc_rho(ispin)%pw , error=error) 01686 CALL pw_pool_give_back_pw(pw_pool, tmp_r, error=error) 01687 CALL pw_pool_give_back_pw(pw_pool,vxc_g,error=error) 01688 ELSE 01689 tmp_r => vxc_rho(ispin)%pw 01690 DO idir=1,3 01691 SELECT CASE(xc_deriv_method_id) 01692 CASE (xc_deriv_spline2_smooth) 01693 CALL pw_nn_deriv_r ( pw_in=pw_to_deriv(idir)%pw,& 01694 pw_out=tmp_r,coeffs=spline2_deriv_coeffs,& 01695 idir=idir, error=error ) 01696 CASE (xc_deriv_spline3_smooth) 01697 CALL pw_nn_deriv_r ( pw_in=pw_to_deriv(idir)%pw,& 01698 pw_out=tmp_r,coeffs=spline3_deriv_coeffs,& 01699 idir=idir, error=error ) 01700 CASE (xc_deriv_nn10_smooth) 01701 CALL pw_nn_deriv_r ( pw_in=pw_to_deriv(idir)%pw,& 01702 pw_out=tmp_r,coeffs=nn10_deriv_coeffs,& 01703 idir=idir, error=error ) 01704 CASE (xc_deriv_nn50_smooth) 01705 CALL pw_nn_deriv_r ( pw_in=pw_to_deriv(idir)%pw,& 01706 pw_out=tmp_r,coeffs=nn50_deriv_coeffs,& 01707 idir=idir, error=error ) 01708 CASE default 01709 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 01710 END SELECT 01711 IF (dealloc_pw_to_deriv) THEN 01712 CALL pw_pool_give_back_pw(pw_pool,pw_to_deriv(idir)%pw,error=error) 01713 END IF 01714 END DO 01715 NULLIFY(tmp_r) 01716 END IF 01717 END IF 01718 01719 ! ** Add laplace part to vxc_rho 01720 IF( has_laplace ) THEN 01721 nd_laplace = RESHAPE((/2,0,0,0,2,0,0,0,2/),(/3,3/)) 01722 IF(lsd) THEN 01723 IF( ispin == 1) THEN 01724 deriv_att => xc_dset_get_derivative(deriv_set, "(laplace_rhoa)", error=error) 01725 CPPrecondition(ASSOCIATED(deriv_att),cp_failure_level,routineP,error,failure) 01726 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,& 01727 error=error) 01728 ELSE 01729 deriv_att => xc_dset_get_derivative(deriv_set, "(laplace_rhob)", error=error) 01730 CPPrecondition(ASSOCIATED(deriv_att),cp_failure_level,routineP,error,failure) 01731 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,& 01732 error=error) 01733 END IF 01734 01735 ELSE 01736 deriv_att => xc_dset_get_derivative(deriv_set, "(laplace_rho)", error=error) 01737 CPPrecondition(ASSOCIATED(deriv_att),cp_failure_level,routineP,error,failure) 01738 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,& 01739 error=error) 01740 END IF 01741 DO idir=1,3 01742 CALL pw_pool_create_pw(pw_pool,pw_to_deriv(idir)%pw,& 01743 use_data=REALDATA3D,in_space=REALSPACE,& 01744 error=error) 01745 CALL pw_zero(pw_to_deriv(idir)%pw, error=error) 01746 !$omp parallel do default(none) shared(bo,deriv_data,& 01747 !$omp pw_to_deriv,idir)& 01748 !$omp private(k,j,i) 01749 DO k = bo(1,3), bo(2,3) 01750 DO j = bo(1,2), bo(2,2) 01751 DO i = bo(1,1), bo(2,1) 01752 pw_to_deriv(idir)%pw%cr3d(i,j,k)=deriv_data(i,j,k) 01753 END DO 01754 END DO 01755 END DO 01756 CALL pw_pool_create_pw(pw_pool,tmp_g,& 01757 use_data=COMPLEXDATA1D,in_space=RECIPROCALSPACE,& 01758 error=error) 01759 CALL pw_zero(tmp_g, error=error) 01760 CALL pw_transfer ( pw_to_deriv(idir)%pw, tmp_g , error=error) 01761 01762 SELECT CASE(xc_deriv_method_id) 01763 CASE (xc_deriv_pw) 01764 CALL pw_derive ( tmp_g, nd_laplace(:,idir) , error=error) 01765 CASE default 01766 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 01767 END SELECT 01768 ! Add this to the potential 01769 CALL pw_pool_create_pw(pw_pool, tmp_r,& 01770 use_data=REALDATA3D, in_space=REALSPACE,& 01771 error=error) 01772 CALL pw_zero(tmp_r, error=error) 01773 CALL pw_transfer ( tmp_g, tmp_r , error=error) 01774 01775 CALL pw_axpy ( tmp_r, vxc_rho(ispin)%pw , error=error) 01776 CALL pw_pool_give_back_pw(pw_pool, tmp_r, error=error) 01777 CALL pw_pool_give_back_pw(pw_pool,pw_to_deriv(idir)%pw,error=error) 01778 CALL pw_pool_give_back_pw(pw_pool,tmp_g,error=error) 01779 END DO 01780 END IF 01781 01782 01783 IF (pw_grid%spherical) THEN 01784 ! filter vxc 01785 CALL pw_pool_create_pw(pw_pool,vxc_g,& 01786 use_data=COMPLEXDATA1D,in_space=RECIPROCALSPACE,& 01787 error=error) 01788 CALL pw_transfer ( vxc_rho(ispin)%pw, vxc_g , error=error) 01789 CALL pw_transfer ( vxc_g,vxc_rho(ispin)%pw , error=error) 01790 CALL pw_pool_give_back_pw(pw_pool,vxc_g,error=error) 01791 END IF 01792 CALL smooth_cutoff(pot=vxc_rho(ispin)%pw%cr3d,rho=rho,rhoa=rhoa,rhob=rhob,& 01793 rho_cutoff=rho_cutoff*density_smooth_cut_range,& 01794 rho_smooth_cutoff_range=density_smooth_cut_range,& 01795 error=error) 01796 01797 ! final smoothing if rho was smoothed 01798 IF (xc_rho_smooth_id/=xc_rho_no_smooth) THEN 01799 CALL pw_pool_create_pw(pw_pool, tmp_r,& 01800 use_data=REALDATA3D, in_space=REALSPACE,& 01801 error=error) 01802 CALL pw_zero(tmp_r, error=error) 01803 SELECT CASE(xc_rho_smooth_id) 01804 CASE (xc_rho_spline2_smooth) 01805 CALL pw_nn_smear_r(pw_in=vxc_rho(ispin)%pw,pw_out=tmp_r,& 01806 coeffs=spline2_coeffs,error=error) 01807 CASE (xc_rho_spline3_smooth) 01808 CALL pw_nn_smear_r(pw_in=vxc_rho(ispin)%pw,pw_out=tmp_r,& 01809 coeffs=spline3_coeffs,error=error) 01810 CASE (xc_rho_nn10) 01811 CALL pw_nn_smear_r(pw_in=vxc_rho(ispin)%pw,pw_out=tmp_r,& 01812 coeffs=nn10_coeffs,error=error) 01813 CASE (xc_rho_nn50) 01814 CALL pw_nn_smear_r(pw_in=vxc_rho(ispin)%pw,pw_out=tmp_r,& 01815 coeffs=nn50_coeffs,error=error) 01816 CASE default 01817 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 01818 END SELECT 01819 deriv_data => vxc_rho(ispin)%pw%cr3d 01820 vxc_rho(ispin)%pw%cr3d => tmp_r%cr3d 01821 tmp_r%cr3d => deriv_data 01822 CALL pw_pool_give_back_pw(pw_pool,tmp_r,error=error) 01823 END IF 01824 END DO 01825 01826 DO idir=1,3 01827 CPPostcondition(.NOT.ASSOCIATED(pw_to_deriv(idir)%pw),cp_warning_level,routineP,error,failure) 01828 CPPostcondition(.NOT.ASSOCIATED(pw_to_deriv_rho(idir)%pw),cp_warning_level,routineP,error,failure) 01829 END DO 01830 01831 ! 0-deriv -> value of exc 01832 ! this has to be kept consistent with xc_exc_calc 01833 IF (n_deriv>0) THEN 01834 deriv_att => xc_dset_get_derivative(deriv_set, "", error=error) 01835 CPPrecondition(ASSOCIATED(deriv_att),cp_failure_level,routineP,error,failure) 01836 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,& 01837 error=error) 01838 01839 CALL pw_create(tmp_r,pw_grid,& 01840 use_data=REALDATA3D,in_space=REALSPACE,& 01841 cr3d_ptr=deriv_data, error=error) 01842 NULLIFY(tmp_r%cr3d) 01843 CALL pw_release(tmp_r,error=error) 01844 01845 CALL smooth_cutoff(pot=deriv_data,rho=rho,rhoa=rhoa,rhob=rhob,& 01846 rho_cutoff=rho_cutoff,& 01847 rho_smooth_cutoff_range=density_smooth_cut_range,& 01848 error=error) 01849 01850 exc = accurate_sum ( deriv_data )*pw_grid%dvol 01851 IF ( pw_grid%para%mode == PW_MODE_DISTRIBUTED ) THEN 01852 CALL mp_sum ( exc, pw_grid%para%group ) 01853 END IF 01854 ELSE 01855 exc=0.0_dp 01856 END IF 01857 01858 CALL xc_rho_set_release(rho_set, pw_pool=pw_pool, error=error) 01859 01860 ! tau part 01861 IF (has_tau) THEN 01862 ALLOCATE(vxc_tau(nspins), stat=stat) 01863 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01864 DO ispin=1,nspins 01865 NULLIFY(vxc_tau(ispin)%pw) 01866 END DO 01867 IF (lsd) THEN 01868 deriv_att => xc_dset_get_derivative(deriv_set, "(tau_a)", error=error) 01869 IF (ASSOCIATED(deriv_att)) THEN 01870 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,& 01871 error=error) 01872 01873 CALL pw_create(vxc_tau(1)%pw,pw_grid,& 01874 use_data=REALDATA3D,in_space=REALSPACE,& 01875 cr3d_ptr=deriv_data, error=error) 01876 NULLIFY(deriv_att%deriv_data) 01877 END IF 01878 deriv_att => xc_dset_get_derivative(deriv_set, "(tau_b)", error=error) 01879 IF (ASSOCIATED(deriv_att)) THEN 01880 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,& 01881 error=error) 01882 01883 CALL pw_create(vxc_tau(2)%pw,pw_grid,& 01884 use_data=REALDATA3D,in_space=REALSPACE,& 01885 cr3d_ptr=deriv_data, error=error) 01886 NULLIFY(deriv_att%deriv_data) 01887 END IF 01888 deriv_att => xc_dset_get_derivative(deriv_set, "(tau)", error=error) 01889 IF (ASSOCIATED(deriv_att)) THEN 01890 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,& 01891 error=error) 01892 IF (ASSOCIATED(vxc_tau(1)%pw)) THEN 01893 DO ispin=1,2 01894 CPPrecondition(ASSOCIATED(vxc_tau(ispin)%pw),cp_failure_level,routineP,error,failure) 01895 tmp_cr3d => vxc_tau(ispin)%pw%cr3d 01896 CALL daxpy(npoints,1.0_dp,deriv_data,1,tmp_cr3d,1) 01897 END DO 01898 ELSE 01899 CALL pw_create(vxc_tau(1)%pw,pw_grid,& 01900 use_data=REALDATA3D,in_space=REALSPACE,& 01901 cr3d_ptr=deriv_data, error=error) 01902 NULLIFY(deriv_att%deriv_data) 01903 CALL pw_pool_create_pw(pw_pool, vxc_tau(2)%pw,& 01904 use_data=REALDATA3D, in_space=REALSPACE,& 01905 error=error) 01906 CALL pw_copy(vxc_tau(1)%pw,vxc_tau(2)%pw, error=error) 01907 END IF 01908 END IF 01909 ELSE 01910 deriv_att => xc_dset_get_derivative(deriv_set, "(tau)", error=error) 01911 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,& 01912 error=error) 01913 CALL pw_create(vxc_tau(1)%pw,pw_grid,& 01914 use_data=REALDATA3D,in_space=REALSPACE,& 01915 cr3d_ptr=deriv_data, error=error) 01916 NULLIFY(deriv_att%deriv_data) 01917 END IF 01918 DO ispin=1,nspins 01919 CPPostcondition(ASSOCIATED(vxc_tau(ispin)%pw),cp_failure_level,routineP,error,failure) 01920 END DO 01921 END IF 01922 CALL xc_dset_release(deriv_set, error=error) 01923 END IF 01924 01925 CALL timestop(handle) 01926 01927 END SUBROUTINE xc_vxc_pw_create 01928 01929 ! ***************************************************************************** 01940 FUNCTION xc_exc_calc(rho_r,rho_g,tau,xc_section, cell,pw_pool,& 01941 error) RESULT(exc) 01942 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r, rho_g, tau 01943 TYPE(section_vals_type), POINTER :: xc_section 01944 TYPE(cell_type), POINTER :: cell 01945 TYPE(pw_pool_type), POINTER :: pw_pool 01946 TYPE(cp_error_type), INTENT(inout) :: error 01947 REAL(kind=dp) :: exc 01948 01949 CHARACTER(len=*), PARAMETER :: routineN = 'xc_exc_calc', 01950 routineP = moduleN//':'//routineN 01951 01952 INTEGER :: handle 01953 LOGICAL :: failure 01954 REAL(dp) :: density_smooth_cut_range, 01955 rho_cutoff 01956 REAL(dp), DIMENSION(:, :, :), POINTER :: e_0 01957 TYPE(xc_derivative_set_type), POINTER :: deriv_set 01958 TYPE(xc_derivative_type), POINTER :: deriv 01959 TYPE(xc_rho_set_type), POINTER :: rho_set 01960 01961 CALL timeset(routineN,handle) 01962 NULLIFY(rho_set, deriv_set, deriv,e_0) 01963 failure=.FALSE. 01964 exc=0.0_dp 01965 01966 ! this has to be consistent with what is done in xc_vxc_create 01967 CALL xc_rho_set_and_dset_create(rho_set=rho_set,& 01968 deriv_set=deriv_set,deriv_order=0,& 01969 rho_r=rho_r,rho_g=rho_g,tau=tau,xc_section=xc_section,& 01970 cell=cell,pw_pool=pw_pool,& 01971 needs_basic_components=.FALSE.,error=error) 01972 deriv => xc_dset_get_derivative(deriv_set,"",error=error) 01973 01974 IF (ASSOCIATED(deriv)) THEN 01975 CALL xc_derivative_get(deriv, deriv_data=e_0, error=error) 01976 01977 CALL section_vals_val_get(xc_section,"DENSITY_CUTOFF",& 01978 r_val=rho_cutoff,error=error) 01979 CALL section_vals_val_get(xc_section,"DENSITY_SMOOTH_CUTOFF_RANGE",& 01980 r_val=density_smooth_cut_range,error=error) 01981 CALL smooth_cutoff(pot=e_0,rho=rho_set%rho,& 01982 rhoa=rho_set%rhoa,rhob=rho_set%rhob,& 01983 rho_cutoff=rho_cutoff,& 01984 rho_smooth_cutoff_range=density_smooth_cut_range,& 01985 error=error) 01986 01987 exc = accurate_sum ( e_0 )*rho_r(1)%pw%pw_grid%dvol 01988 IF ( rho_r(1)%pw%pw_grid%para%mode == PW_MODE_DISTRIBUTED ) THEN 01989 CALL mp_sum ( exc, rho_r(1)%pw%pw_grid%para%group ) 01990 END IF 01991 01992 CALL xc_rho_set_release(rho_set, pw_pool=pw_pool, error=error) 01993 CALL xc_dset_release(deriv_set, error=error) 01994 END IF 01995 CALL timestop(handle) 01996 END FUNCTION xc_exc_calc 01997 01998 ! ***************************************************************************** 02032 SUBROUTINE xc_calc_2nd_deriv(v_xc, deriv_set, rho_set, rho1_set, & 02033 pw_pool, xc_section, gapw, vxg, tddfpt_fac, error) 02034 02035 TYPE(pw_p_type), DIMENSION(:), POINTER :: v_xc 02036 TYPE(xc_derivative_set_type), POINTER :: deriv_set 02037 TYPE(xc_rho_set_type), POINTER :: rho_set, rho1_set 02038 TYPE(pw_pool_type), POINTER :: pw_pool 02039 TYPE(section_vals_type), POINTER :: xc_section 02040 LOGICAL, INTENT(IN), OPTIONAL :: gapw 02041 REAL(kind=dp), DIMENSION(:, :, :, :), 02042 OPTIONAL, POINTER :: vxg 02043 REAL(kind=dp), INTENT(in), OPTIONAL :: tddfpt_fac 02044 TYPE(cp_error_type), INTENT(inout) :: error 02045 02046 CHARACTER(len=*), PARAMETER :: routineN = 'xc_calc_2nd_deriv', 02047 routineP = moduleN//':'//routineN 02048 02049 CHARACTER& 02050 (len=MAX_DERIVATIVE_DESC_LENGTH) :: desc 02051 CHARACTER(len=MAX_LABEL_LENGTH), 02052 DIMENSION(:), POINTER :: split_desc 02053 INTEGER :: handle, i, ia, idir, ir, 02054 ispin, j, k, nspins, order, 02055 stat, xc_deriv_method_id 02056 INTEGER, DIMENSION(2, 3) :: bo 02057 LOGICAL :: failure, gradient_f, lsd, 02058 my_gapw, unknown_deriv 02059 REAL(KIND=dp) :: dr1dr, fac, gradient_cut 02060 REAL(kind=dp), DIMENSION(:, :, :), 02061 POINTER :: deriv_data, e_drhoa, e_drhob, 02062 e_norm_drho, rho1, rho1a, 02063 rho1b 02064 TYPE(cp_3d_r_p_type), DIMENSION(:), 02065 POINTER :: drho, drho1, drho1a, drho1b, 02066 drhoa, drhob 02067 TYPE(cp_sll_xc_deriv_type), POINTER :: pos 02068 TYPE(pw_p_type) :: v_drho 02069 TYPE(pw_p_type), DIMENSION(:), POINTER :: tmp_a, tmp_b, tmp_r, tmp_r2 02070 TYPE(xc_derivative_type), POINTER :: deriv_att 02071 02072 ! PARAMETER 02073 02074 CALL timeset(routineN,handle) 02075 02076 failure=.FALSE. 02077 NULLIFY(tmp_r, tmp_r2, tmp_a, tmp_b, e_drhoa, e_drhob, e_norm_drho, & 02078 split_desc) 02079 02080 my_gapw = .FALSE. 02081 IF (PRESENT(gapw)) my_gapw = gapw 02082 02083 CPPrecondition(ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 02084 CPPrecondition(ASSOCIATED(rho1_set),cp_failure_level,routineP,error,failure) 02085 CPPrecondition(ASSOCIATED(v_xc),cp_failure_level,routineP,error,failure) 02086 CPPrecondition(ASSOCIATED(xc_section),cp_failure_level,routineP,error,failure) 02087 IF (my_gapw) THEN 02088 CPPrecondition(PRESENT(vxg),cp_failure_level,routineP,error,failure) 02089 END IF 02090 02091 IF (.NOT. failure) THEN 02092 CALL section_vals_val_get(xc_section,"XC_GRID%XC_DERIV",& 02093 i_val=xc_deriv_method_id,error=error) 02094 CALL xc_rho_set_get(rho_set,drho_cutoff=gradient_cut,error=error) 02095 nspins = SIZE(v_xc) 02096 lsd = ASSOCIATED(rho_set%rhoa) 02097 fac = 0.0_dp 02098 IF (PRESENT(tddfpt_fac)) fac=tddfpt_fac 02099 02100 ALLOCATE(tmp_r(nspins), tmp_r2(nspins), tmp_a(nspins), tmp_b(nspins),stat=stat) 02101 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02102 DO ispin=1, nspins 02103 NULLIFY(tmp_r(ispin)%pw, tmp_r2(ispin)%pw, tmp_a(ispin)%pw, tmp_b(nspins)%pw) 02104 END DO 02105 02106 END IF 02107 02108 bo = rho_set%local_bounds 02109 02110 gradient_f=.FALSE. 02111 ! check for unknown derivatives 02112 pos => deriv_set%derivs 02113 DO WHILE (cp_sll_xc_deriv_next(pos,el_att=deriv_att,error=error)) 02114 unknown_deriv=.FALSE. 02115 CALL xc_derivative_get(deriv_att,order=order,& 02116 desc=desc, split_desc=split_desc, error=error) 02117 SELECT CASE(order) 02118 CASE(1) 02119 IF (lsd) THEN 02120 SELECT CASE(split_desc(1)) 02121 CASE("rho","rhoa","rhob") 02122 CASE("norm_drho","norm_drhoa","norm_drhob") 02123 gradient_f=.TRUE. 02124 CASE default 02125 unknown_deriv=.TRUE. 02126 END SELECT 02127 ELSE 02128 SELECT CASE(split_desc(1)) 02129 CASE("rho") 02130 CASE("norm_drho") 02131 gradient_f=.TRUE. 02132 CASE default 02133 unknown_deriv=.TRUE. 02134 END SELECT 02135 END IF 02136 CASE(2) 02137 IF (lsd) THEN 02138 SELECT CASE(split_desc(1)) 02139 CASE("rhoa","rhob") 02140 SELECT CASE(split_desc(2)) 02141 CASE("rhoa","rhob") 02142 CASE("norm_drhoa", "norm_drhob","norm_drho") 02143 gradient_f=.TRUE. 02144 CASE default 02145 unknown_deriv=.TRUE. 02146 END SELECT 02147 CASE("norm_drho","norm_drhoa", "norm_drhob") 02148 gradient_f=.TRUE. 02149 SELECT CASE(split_desc(2)) 02150 CASE("rhoa","rhob","norm_drhoa", "norm_drhob","norm_drho") 02151 CASE default 02152 unknown_deriv=.TRUE. 02153 END SELECT 02154 CASE default 02155 unknown_deriv=.TRUE. 02156 END SELECT 02157 ELSE 02158 SELECT CASE(split_desc(1)) 02159 CASE("rho") 02160 SELECT CASE(split_desc(2)) 02161 CASE("rho") 02162 CASE("norm_drho") 02163 gradient_f=.TRUE. 02164 CASE default 02165 unknown_deriv=.TRUE. 02166 END SELECT 02167 CASE("norm_drho") 02168 gradient_f=.TRUE. 02169 SELECT CASE(split_desc(2)) 02170 CASE("rho","norm_drho") 02171 CASE default 02172 unknown_deriv=.TRUE. 02173 END SELECT 02174 CASE default 02175 unknown_deriv=.TRUE. 02176 END SELECT 02177 END IF 02178 END SELECT 02179 CALL cp_assert(.NOT.unknown_deriv,cp_failure_level,& 02180 cp_assertion_failed,routineP,& 02181 "unknown functional derivative (LSD="//TRIM(cp_to_string(lsd))//& 02182 "): '"//TRIM(desc)//"' in "//& 02183 CPSourceFileRef,& 02184 error,failure) 02185 END DO 02186 02187 IF (lsd) THEN 02188 02189 !-------------------! 02190 ! UNrestricted case ! 02191 !-------------------! 02192 02193 CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, error=error) 02194 02195 IF (gradient_f) THEN 02196 CALL xc_rho_set_get(rho_set, drhoa=drhoa, drhob=drhob, error=error) 02197 CALL xc_rho_set_get(rho1_set,drhoa=drho1a, drhob=drho1b, error=error) 02198 DO ispin=1, nspins 02199 IF (ASSOCIATED(pw_pool)) THEN 02200 CALL pw_pool_create_pw(pw_pool,tmp_a(ispin)%pw,& 02201 use_data = REALDATA3D,& 02202 in_space = REALSPACE, error=error) 02203 CALL pw_zero(tmp_a(ispin)%pw, error=error) 02204 CALL pw_pool_create_pw(pw_pool,tmp_b(ispin)%pw,& 02205 use_data = REALDATA3D,& 02206 in_space = REALSPACE, error=error) 02207 CALL pw_zero(tmp_b(ispin)%pw, error=error) 02208 CALL pw_pool_create_pw(pw_pool,tmp_r(ispin)%pw,& 02209 use_data = REALDATA3D,& 02210 in_space = REALSPACE, error=error) 02211 CALL pw_zero(tmp_r(ispin)%pw, error=error) 02212 ELSE 02213 ALLOCATE(tmp_a(ispin)%pw, tmp_b(ispin)%pw, tmp_r(ispin)%pw, stat=stat) 02214 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02215 ALLOCATE(tmp_a(ispin)%pw%cr3d(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3)), & 02216 tmp_b(ispin)%pw%cr3d(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3)), & 02217 tmp_r(ispin)%pw%cr3d(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3)), & 02218 stat=stat) 02219 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02220 tmp_a(ispin)%pw%cr3d = 0.0_dp 02221 tmp_b(ispin)%pw%cr3d = 0.0_dp 02222 END IF 02223 END DO 02224 END IF 02225 02226 deriv_att=> xc_dset_get_derivative(deriv_set, "(rhoa)(rhoa)", & 02227 error=error) 02228 IF (ASSOCIATED(deriv_att)) THEN 02229 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02230 error=error) 02231 !$omp parallel do private(k,j,i) 02232 DO k = bo(1,3), bo(2,3) 02233 DO j = bo(1,2), bo(2,2) 02234 DO i = bo(1,1), bo(2,1) 02235 v_xc(1)%pw%cr3d(i,j,k) = v_xc(1)%pw%cr3d(i,j,k) + & 02236 deriv_data(i,j,k) * rho1a(i,j,k) 02237 END DO 02238 END DO 02239 END DO 02240 END IF 02241 02242 IF (nspins /= 1) THEN 02243 deriv_att=> xc_dset_get_derivative(deriv_set, "(rhob)(rhob)", & 02244 error=error) 02245 IF (ASSOCIATED(deriv_att)) THEN 02246 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02247 error=error) 02248 !$omp parallel do private(k,j,i) 02249 DO k = bo(1,3), bo(2,3) 02250 DO j = bo(1,2), bo(2,2) 02251 DO i = bo(1,1), bo(2,1) 02252 v_xc(2)%pw%cr3d(i,j,k) = v_xc(2)%pw%cr3d(i,j,k) + & 02253 deriv_data(i,j,k) * rho1b(i,j,k) 02254 END DO 02255 END DO 02256 END DO 02257 END IF 02258 END IF 02259 02260 deriv_att=> xc_dset_get_derivative(deriv_set, "(rhoa)(rhob)", & 02261 error=error) 02262 IF (ASSOCIATED(deriv_att)) THEN 02263 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02264 error=error) 02265 !$omp parallel do private(k,j,i) 02266 DO k = bo(1,3), bo(2,3) 02267 DO j = bo(1,2), bo(2,2) 02268 DO i = bo(1,1), bo(2,1) 02269 IF (nspins /= 1) THEN 02270 v_xc(1)%pw%cr3d(i,j,k) = v_xc(1)%pw%cr3d(i,j,k) + & 02271 deriv_data(i,j,k) * rho1b(i,j,k) 02272 v_xc(2)%pw%cr3d(i,j,k) = v_xc(2)%pw%cr3d(i,j,k) + & 02273 deriv_data(i,j,k) * rho1a(i,j,k) 02274 ELSE 02275 v_xc(1)%pw%cr3d(i,j,k) = v_xc(1)%pw%cr3d(i,j,k) + & 02276 fac * deriv_data(i,j,k) * rho1b(i,j,k) 02277 END IF 02278 END DO 02279 END DO 02280 END DO 02281 END IF 02282 02283 deriv_att=> xc_dset_get_derivative(deriv_set, "(rhoa)(norm_drhoa)", & 02284 error=error) 02285 IF (ASSOCIATED(deriv_att)) THEN 02286 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02287 error=error) 02288 !$omp parallel do private(k,j,i,dr1dr,idir) 02289 DO k = bo(1,3), bo(2,3) 02290 DO j = bo(1,2), bo(2,2) 02291 DO i = bo(1,1), bo(2,1) 02292 dr1dr=0._dp 02293 DO idir=1,3 02294 dr1dr=dr1dr+drhoa(idir)%array(i,j,k)*drho1a(idir)%array(i,j,k) 02295 END DO 02296 v_xc(1)%pw%cr3d(i,j,k) = v_xc(1)%pw%cr3d(i,j,k) + & 02297 deriv_data(i,j,k) * dr1dr 02298 tmp_a(1)%pw%cr3d(i,j,k) = tmp_a(1)%pw%cr3d(i,j,k) - & 02299 deriv_data(i,j,k) * rho1a(i,j,k) 02300 END DO 02301 END DO 02302 END DO 02303 END IF 02304 02305 IF (nspins /= 1) THEN 02306 deriv_att=> xc_dset_get_derivative(deriv_set, "(rhob)(norm_drhob)", error=error) 02307 IF (ASSOCIATED(deriv_att)) THEN 02308 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02309 error=error) 02310 !$omp parallel do private(k,j,i,dr1dr,idir) 02311 DO k = bo(1,3), bo(2,3) 02312 DO j = bo(1,2), bo(2,2) 02313 DO i = bo(1,1), bo(2,1) 02314 dr1dr=0._dp 02315 DO idir=1,3 02316 dr1dr=dr1dr+drhob(idir)%array(i,j,k)*drho1b(idir)%array(i,j,k) 02317 END DO 02318 v_xc(2)%pw%cr3d(i,j,k) = v_xc(2)%pw%cr3d(i,j,k) + & 02319 deriv_data(i,j,k) * dr1dr 02320 tmp_b(2)%pw%cr3d(i,j,k) = tmp_b(2)%pw%cr3d(i,j,k) - & 02321 deriv_data(i,j,k) * rho1b(i,j,k) 02322 END DO 02323 END DO 02324 END DO 02325 END IF 02326 END IF 02327 02328 deriv_att=> xc_dset_get_derivative(deriv_set, "(rhoa)(norm_drhob)", & 02329 error=error) 02330 IF (ASSOCIATED(deriv_att)) THEN 02331 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02332 error=error) 02333 !$omp parallel do private(k,j,i,dr1dr,idir) 02334 DO k = bo(1,3), bo(2,3) 02335 DO j = bo(1,2), bo(2,2) 02336 DO i = bo(1,1), bo(2,1) 02337 dr1dr=0._dp 02338 DO idir=1,3 02339 dr1dr=dr1dr+drhob(idir)%array(i,j,k)*drho1b(idir)%array(i,j,k) 02340 END DO 02341 IF (nspins /= 1) THEN 02342 v_xc(1)%pw%cr3d(i,j,k) = v_xc(1)%pw%cr3d(i,j,k) + & 02343 deriv_data(i,j,k) * dr1dr 02344 tmp_b(2)%pw%cr3d(i,j,k) = tmp_b(2)%pw%cr3d(i,j,k) - & 02345 deriv_data(i,j,k) * rho1a(i,j,k) 02346 ELSE 02347 v_xc(1)%pw%cr3d(i,j,k) = v_xc(1)%pw%cr3d(i,j,k) + & 02348 fac * deriv_data(i,j,k) * dr1dr 02349 END IF 02350 END DO 02351 END DO 02352 END DO 02353 END IF 02354 02355 deriv_att=> xc_dset_get_derivative(deriv_set, "(rhob)(norm_drhoa)", & 02356 error=error) 02357 IF (ASSOCIATED(deriv_att)) THEN 02358 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02359 error=error) 02360 !$omp parallel do private(k,j,i,dr1dr,idir) 02361 DO k = bo(1,3), bo(2,3) 02362 DO j = bo(1,2), bo(2,2) 02363 DO i = bo(1,1), bo(2,1) 02364 IF (nspins /= 1) THEN 02365 dr1dr=0._dp 02366 DO idir=1,3 02367 dr1dr=dr1dr+drhoa(idir)%array(i,j,k)*drho1a(idir)%array(i,j,k) 02368 END DO 02369 v_xc(2)%pw%cr3d(i,j,k) = v_xc(2)%pw%cr3d(i,j,k) + & 02370 deriv_data(i,j,k) * dr1dr 02371 tmp_a(1)%pw%cr3d(i,j,k) = tmp_a(1)%pw%cr3d(i,j,k) - & 02372 deriv_data(i,j,k) * rho1b(i,j,k) 02373 ELSE 02374 tmp_a(1)%pw%cr3d(i,j,k) = tmp_a(1)%pw%cr3d(i,j,k) - & 02375 fac * deriv_data(i,j,k) * rho1b(i,j,k) 02376 END IF 02377 END DO 02378 END DO 02379 END DO 02380 END IF 02381 02382 deriv_att=> xc_dset_get_derivative(deriv_set, "(rhoa)(norm_drho)", & 02383 error=error) 02384 IF (ASSOCIATED(deriv_att)) THEN 02385 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02386 error=error) 02387 !$omp parallel do private(k,j,i,dr1dr,idir) 02388 DO k = bo(1,3), bo(2,3) 02389 DO j = bo(1,2), bo(2,2) 02390 DO i = bo(1,1), bo(2,1) 02391 IF (nspins /= 1) THEN 02392 dr1dr=0._dp 02393 DO idir=1,3 02394 dr1dr=dr1dr+drhob(idir)%array(i,j,k)*drho1a(idir)%array(i,j,k) 02395 END DO 02396 v_xc(1)%pw%cr3d(i,j,k) = v_xc(1)%pw%cr3d(i,j,k) + & 02397 deriv_data(i,j,k) * dr1dr 02398 dr1dr=0._dp 02399 DO idir=1,3 02400 dr1dr=dr1dr+drhoa(idir)%array(i,j,k)*drho1b(idir)%array(i,j,k) 02401 END DO 02402 v_xc(1)%pw%cr3d(i,j,k) = v_xc(1)%pw%cr3d(i,j,k) + & 02403 deriv_data(i,j,k) * dr1dr 02404 tmp_b(1)%pw%cr3d(i,j,k) = tmp_b(1)%pw%cr3d(i,j,k) - & 02405 deriv_data(i,j,k) * rho1a(i,j,k) 02406 tmp_a(2)%pw%cr3d(i,j,k) = tmp_a(2)%pw%cr3d(i,j,k) - & 02407 deriv_data(i,j,k) * rho1a(i,j,k) 02408 ELSE 02409 dr1dr=0._dp 02410 DO idir=1,3 02411 dr1dr = dr1dr + drhob(idir)%array(i,j,k)*drho1a(idir)%array(i,j,k) + & 02412 fac * drhoa(idir)%array(i,j,k)*drho1b(idir)%array(i,j,k) 02413 END DO 02414 v_xc(1)%pw%cr3d(i,j,k) = v_xc(1)%pw%cr3d(i,j,k) + & 02415 deriv_data(i,j,k) * dr1dr 02416 tmp_b(1)%pw%cr3d(i,j,k) = tmp_b(1)%pw%cr3d(i,j,k) - & 02417 deriv_data(i,j,k) * rho1a(i,j,k) 02418 END IF 02419 END DO 02420 END DO 02421 END DO 02422 END IF 02423 02424 deriv_att=> xc_dset_get_derivative(deriv_set, "(norm_drhoa)(norm_drho)", & 02425 error=error) 02426 IF (ASSOCIATED(deriv_att)) THEN 02427 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02428 error=error) 02429 !$omp parallel do private(k,j,i,dr1dr,idir) 02430 DO k = bo(1,3), bo(2,3) 02431 DO j = bo(1,2), bo(2,2) 02432 DO i = bo(1,1), bo(2,1) 02433 IF (nspins /= 1) THEN 02434 dr1dr=0._dp 02435 DO idir=1,3 02436 dr1dr=dr1dr+drhob(idir)%array(i,j,k)*drho1a(idir)%array(i,j,k) 02437 END DO 02438 tmp_a(1)%pw%cr3d(i,j,k) = tmp_a(1)%pw%cr3d(i,j,k) - & 02439 deriv_data(i,j,k) * dr1dr 02440 dr1dr=0._dp 02441 DO idir=1,3 02442 dr1dr=dr1dr+drhoa(idir)%array(i,j,k)*drho1b(idir)%array(i,j,k) 02443 END DO 02444 tmp_a(1)%pw%cr3d(i,j,k) = tmp_a(1)%pw%cr3d(i,j,k) - & 02445 deriv_data(i,j,k) * dr1dr 02446 dr1dr=0._dp 02447 DO idir=1,3 02448 dr1dr=dr1dr+drhoa(idir)%array(i,j,k)*drho1a(idir)%array(i,j,k) 02449 END DO 02450 tmp_b(1)%pw%cr3d(i,j,k) = tmp_b(1)%pw%cr3d(i,j,k) - & 02451 deriv_data(i,j,k) * dr1dr 02452 dr1dr=0._dp 02453 DO idir=1,3 02454 dr1dr=dr1dr+drhoa(idir)%array(i,j,k)*drho1a(idir)%array(i,j,k) 02455 END DO 02456 tmp_a(2)%pw%cr3d(i,j,k) = tmp_a(2)%pw%cr3d(i,j,k) - & 02457 deriv_data(i,j,k) * dr1dr 02458 ELSE 02459 dr1dr=0._dp 02460 DO idir=1,3 02461 dr1dr=dr1dr+drhob(idir)%array(i,j,k)*drho1a(idir)%array(i,j,k) + & 02462 fac * drhoa(idir)%array(i,j,k)*drho1b(idir)%array(i,j,k) 02463 END DO 02464 tmp_a(1)%pw%cr3d(i,j,k) = tmp_a(1)%pw%cr3d(i,j,k) - & 02465 deriv_data(i,j,k) * dr1dr 02466 dr1dr=0._dp 02467 DO idir=1,3 02468 dr1dr=dr1dr+drhoa(idir)%array(i,j,k)*drho1a(idir)%array(i,j,k) 02469 END DO 02470 tmp_b(1)%pw%cr3d(i,j,k) = tmp_b(1)%pw%cr3d(i,j,k) - & 02471 deriv_data(i,j,k) * dr1dr 02472 END IF 02473 END DO 02474 END DO 02475 END DO 02476 END IF 02477 02478 deriv_att=> xc_dset_get_derivative(deriv_set, "(rhob)(norm_drho)", & 02479 error=error) 02480 IF (ASSOCIATED(deriv_att)) THEN 02481 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02482 error=error) 02483 !$omp parallel do private(k,j,i,dr1dr,idir) 02484 DO k = bo(1,3), bo(2,3) 02485 DO j = bo(1,2), bo(2,2) 02486 DO i = bo(1,1), bo(2,1) 02487 IF (nspins /= 1) THEN 02488 dr1dr=0._dp 02489 DO idir=1,3 02490 dr1dr=dr1dr+drhoa(idir)%array(i,j,k)*drho1b(idir)%array(i,j,k) 02491 END DO 02492 v_xc(2)%pw%cr3d(i,j,k) = v_xc(2)%pw%cr3d(i,j,k) + & 02493 deriv_data(i,j,k) * dr1dr 02494 dr1dr=0._dp 02495 DO idir=1,3 02496 dr1dr=dr1dr+drhob(idir)%array(i,j,k)*drho1a(idir)%array(i,j,k) 02497 END DO 02498 v_xc(2)%pw%cr3d(i,j,k) = v_xc(2)%pw%cr3d(i,j,k) + & 02499 deriv_data(i,j,k) * dr1dr 02500 tmp_a(2)%pw%cr3d(i,j,k) = tmp_a(2)%pw%cr3d(i,j,k) - & 02501 deriv_data(i,j,k) * rho1b(i,j,k) 02502 tmp_b(1)%pw%cr3d(i,j,k) = tmp_b(1)%pw%cr3d(i,j,k) - & 02503 deriv_data(i,j,k) * rho1b(i,j,k) 02504 ELSE 02505 tmp_b(1)%pw%cr3d(i,j,k) = tmp_b(1)%pw%cr3d(i,j,k) - & 02506 fac * deriv_data(i,j,k) * rho1b(i,j,k) 02507 END IF 02508 END DO 02509 END DO 02510 END DO 02511 END IF 02512 02513 deriv_att=> xc_dset_get_derivative(deriv_set, "(norm_drhob)(norm_drho)", & 02514 error=error) 02515 IF (ASSOCIATED(deriv_att)) THEN 02516 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02517 error=error) 02518 !$omp parallel do private(k,j,i,dr1dr,idir) 02519 DO k = bo(1,3), bo(2,3) 02520 DO j = bo(1,2), bo(2,2) 02521 DO i = bo(1,1), bo(2,1) 02522 IF (nspins /= 1) THEN 02523 dr1dr=0._dp 02524 DO idir=1,3 02525 dr1dr=dr1dr+drhoa(idir)%array(i,j,k)*drho1b(idir)%array(i,j,k) 02526 END DO 02527 tmp_b(2)%pw%cr3d(i,j,k) = tmp_b(2)%pw%cr3d(i,j,k) - & 02528 deriv_data(i,j,k) * dr1dr 02529 dr1dr=0._dp 02530 DO idir=1,3 02531 dr1dr=dr1dr+drhob(idir)%array(i,j,k)*drho1a(idir)%array(i,j,k) 02532 END DO 02533 tmp_b(2)%pw%cr3d(i,j,k) = tmp_b(2)%pw%cr3d(i,j,k) - & 02534 deriv_data(i,j,k) * dr1dr 02535 dr1dr=0._dp 02536 DO idir=1,3 02537 dr1dr=dr1dr+drhob(idir)%array(i,j,k)*drho1b(idir)%array(i,j,k) 02538 END DO 02539 tmp_a(2)%pw%cr3d(i,j,k) = tmp_a(2)%pw%cr3d(i,j,k) - & 02540 deriv_data(i,j,k) * dr1dr 02541 dr1dr=0._dp 02542 DO idir=1,3 02543 dr1dr=dr1dr+drhob(idir)%array(i,j,k)*drho1b(idir)%array(i,j,k) 02544 END DO 02545 tmp_b(1)%pw%cr3d(i,j,k) = tmp_b(1)%pw%cr3d(i,j,k) - & 02546 deriv_data(i,j,k) * dr1dr 02547 ELSE 02548 dr1dr=0._dp 02549 DO idir=1,3 02550 dr1dr=dr1dr+drhob(idir)%array(i,j,k)*drho1b(idir)%array(i,j,k) 02551 END DO 02552 tmp_b(1)%pw%cr3d(i,j,k) = tmp_b(1)%pw%cr3d(i,j,k) - & 02553 fac*deriv_data(i,j,k) * dr1dr 02554 END IF 02555 END DO 02556 END DO 02557 END DO 02558 END IF 02559 02560 deriv_att=> xc_dset_get_derivative(deriv_set, "(norm_drhoa)(norm_drhoa)", & 02561 error=error) 02562 IF (ASSOCIATED(deriv_att)) THEN 02563 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02564 error=error) 02565 !$omp parallel do private(k,j,i,dr1dr,idir) 02566 DO k = bo(1,3), bo(2,3) 02567 DO j = bo(1,2), bo(2,2) 02568 DO i = bo(1,1), bo(2,1) 02569 dr1dr=0._dp 02570 DO idir=1,3 02571 dr1dr=dr1dr+drhoa(idir)%array(i,j,k)*drho1a(idir)%array(i,j,k) 02572 END DO 02573 tmp_a(1)%pw%cr3d(i,j,k) = tmp_a(1)%pw%cr3d(i,j,k) - & 02574 deriv_data(i,j,k) * dr1dr 02575 END DO 02576 END DO 02577 END DO 02578 END IF 02579 02580 IF (nspins /= 1) THEN 02581 deriv_att=> xc_dset_get_derivative(deriv_set, "(norm_drhob)(norm_drhob)", error=error) 02582 IF (ASSOCIATED(deriv_att)) THEN 02583 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02584 error=error) 02585 !$omp parallel do private(k,j,i,dr1dr,idir) 02586 DO k = bo(1,3), bo(2,3) 02587 DO j = bo(1,2), bo(2,2) 02588 DO i = bo(1,1), bo(2,1) 02589 dr1dr=0._dp 02590 DO idir=1,3 02591 dr1dr=dr1dr+drhob(idir)%array(i,j,k)*drho1b(idir)%array(i,j,k) 02592 END DO 02593 tmp_b(2)%pw%cr3d(i,j,k) = tmp_b(2)%pw%cr3d(i,j,k) - & 02594 deriv_data(i,j,k) * dr1dr 02595 END DO 02596 END DO 02597 END DO 02598 END IF 02599 END IF 02600 02601 deriv_att=> xc_dset_get_derivative(deriv_set, "(norm_drhoa)(norm_drhob)", & 02602 error=error) 02603 IF (ASSOCIATED(deriv_att)) THEN 02604 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02605 error=error) 02606 !$omp parallel do private(k,j,i,dr1dr,idir) 02607 DO k = bo(1,3), bo(2,3) 02608 DO j = bo(1,2), bo(2,2) 02609 DO i = bo(1,1), bo(2,1) 02610 dr1dr=0._dp 02611 DO idir=1,3 02612 dr1dr=dr1dr+drhob(idir)%array(i,j,k)*drho1b(idir)%array(i,j,k) 02613 END DO 02614 IF (nspins /= 1) THEN 02615 tmp_a(1)%pw%cr3d(i,j,k) = tmp_a(1)%pw%cr3d(i,j,k) - & 02616 deriv_data(i,j,k) * dr1dr 02617 dr1dr=0._dp 02618 DO idir=1,3 02619 dr1dr=dr1dr+drhoa(idir)%array(i,j,k)*drho1a(idir)%array(i,j,k) 02620 END DO 02621 tmp_b(2)%pw%cr3d(i,j,k) = tmp_b(2)%pw%cr3d(i,j,k) - & 02622 deriv_data(i,j,k) * dr1dr 02623 ELSE 02624 tmp_a(1)%pw%cr3d(i,j,k) = tmp_a(1)%pw%cr3d(i,j,k) - & 02625 fac * deriv_data(i,j,k) * dr1dr 02626 END IF 02627 END DO 02628 END DO 02629 END DO 02630 END IF 02631 02632 deriv_att=> xc_dset_get_derivative(deriv_set, "(norm_drhoa)", & 02633 error=error) 02634 IF (ASSOCIATED(deriv_att)) THEN 02635 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02636 error=error) 02637 CALL xc_derivative_get(deriv_att, deriv_data=e_drhoa, & 02638 error=error) 02639 !$omp parallel do private(k,j,i,dr1dr,idir) 02640 DO k = bo(1,3), bo(2,3) 02641 DO j = bo(1,2), bo(2,2) 02642 DO i = bo(1,1), bo(2,1) 02643 dr1dr=0._dp 02644 DO idir=1,3 02645 dr1dr=dr1dr+drhoa(idir)%array(i,j,k)*drho1a(idir)%array(i,j,k) 02646 END DO 02647 tmp_a(1)%pw%cr3d(i,j,k) = tmp_a(1)%pw%cr3d(i,j,k) + & 02648 deriv_data(i,j,k) * dr1dr 02649 END DO 02650 END DO 02651 END DO 02652 END IF 02653 02654 IF (nspins /= 1) THEN 02655 deriv_att=> xc_dset_get_derivative(deriv_set, "(norm_drhob)", & 02656 error=error) 02657 IF (ASSOCIATED(deriv_att)) THEN 02658 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02659 error=error) 02660 CALL xc_derivative_get(deriv_att, deriv_data=e_drhob, & 02661 error=error) 02662 !$omp parallel do private(k,j,i,dr1dr,idir) 02663 DO k = bo(1,3), bo(2,3) 02664 DO j = bo(1,2), bo(2,2) 02665 DO i = bo(1,1), bo(2,1) 02666 dr1dr=0._dp 02667 DO idir=1,3 02668 dr1dr=dr1dr+drhob(idir)%array(i,j,k)*drho1b(idir)%array(i,j,k) 02669 END DO 02670 tmp_b(2)%pw%cr3d(i,j,k) = tmp_b(2)%pw%cr3d(i,j,k) + & 02671 deriv_data(i,j,k) * dr1dr 02672 END DO 02673 END DO 02674 END DO 02675 END IF 02676 END IF 02677 02678 deriv_att=> xc_dset_get_derivative(deriv_set, "(norm_drho)(norm_drho)", & 02679 error=error) 02680 IF (ASSOCIATED(deriv_att)) THEN 02681 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02682 error=error) 02683 !$omp parallel do private(k,j,i,dr1dr,idir) 02684 DO k = bo(1,3), bo(2,3) 02685 DO j = bo(1,2), bo(2,2) 02686 DO i = bo(1,1), bo(2,1) 02687 IF (nspins /= 1) THEN 02688 dr1dr=0._dp 02689 DO idir=1,3 02690 dr1dr=dr1dr+drhob(idir)%array(i,j,k)*drho1a(idir)%array(i,j,k) 02691 END DO 02692 tmp_b(1)%pw%cr3d(i,j,k) = tmp_b(1)%pw%cr3d(i,j,k) - & 02693 deriv_data(i,j,k) * dr1dr 02694 tmp_a(2)%pw%cr3d(i,j,k) = tmp_a(2)%pw%cr3d(i,j,k) - & 02695 deriv_data(i,j,k) * dr1dr 02696 dr1dr=0._dp 02697 DO idir=1,3 02698 dr1dr=dr1dr+drhoa(idir)%array(i,j,k)*drho1b(idir)%array(i,j,k) 02699 END DO 02700 tmp_b(1)%pw%cr3d(i,j,k) = tmp_b(1)%pw%cr3d(i,j,k) - & 02701 deriv_data(i,j,k) * dr1dr 02702 tmp_a(2)%pw%cr3d(i,j,k) = tmp_a(2)%pw%cr3d(i,j,k) - & 02703 deriv_data(i,j,k) * dr1dr 02704 ELSE 02705 dr1dr=0._dp 02706 DO idir=1,3 02707 dr1dr=dr1dr+drhob(idir)%array(i,j,k)*drho1a(idir)%array(i,j,k) + & 02708 fac * drhoa(idir)%array(i,j,k)*drho1b(idir)%array(i,j,k) 02709 END DO 02710 tmp_b(1)%pw%cr3d(i,j,k) = tmp_b(1)%pw%cr3d(i,j,k) - & 02711 deriv_data(i,j,k) * dr1dr 02712 END IF 02713 END DO 02714 END DO 02715 END DO 02716 END IF 02717 02718 deriv_att=> xc_dset_get_derivative(deriv_set, "(norm_drho)", & 02719 error=error) 02720 IF (ASSOCIATED(deriv_att)) THEN 02721 CALL xc_derivative_get(deriv_att, deriv_data=e_norm_drho, & 02722 error=error) 02723 END IF 02724 02725 IF (gradient_f) THEN 02726 02727 IF (my_gapw) THEN 02728 !$omp parallel do private(ia,idir,ispin,ir) 02729 DO ir = bo(1,2), bo(2,2) 02730 DO ia = bo(1,1), bo(2,1) 02731 DO idir = 1,3 02732 DO ispin=1, nspins 02733 vxg(idir,ia,ir,ispin) = & 02734 tmp_a(ispin)%pw%cr3d(ia,ir,1) * drhoa(idir)%array(ia,ir,1) + & 02735 tmp_b(ispin)%pw%cr3d(ia,ir,1) * drhob(idir)%array(ia,ir,1) 02736 END DO 02737 IF (ASSOCIATED(e_drhoa)) THEN 02738 vxg(idir,ia,ir,1) = vxg(idir,ia,ir,1) - & 02739 e_drhoa(ia,ir,1) * drho1a(idir)%array(ia,ir,1) 02740 END IF 02741 IF (nspins /= 1 .AND. ASSOCIATED(e_drhob)) THEN 02742 vxg(idir,ia,ir,2) = vxg(idir,ia,ir,2) - & 02743 e_drhob(ia,ir,1) * drho1b(idir)%array(ia,ir,1) 02744 END IF 02745 IF (ASSOCIATED(e_norm_drho)) THEN 02746 IF (nspins /= 1) THEN 02747 vxg(idir,ia,ir,1) = vxg(idir,ia,ir,1) - & 02748 e_norm_drho(ia,ir,1) * drho1b(idir)%array(ia,ir,1) 02749 vxg(idir,ia,ir,2) = vxg(idir,ia,ir,2) - & 02750 e_norm_drho(ia,ir,1) * drho1a(idir)%array(ia,ir,1) 02751 ELSE 02752 vxg(idir,ia,ir,1) = vxg(idir,ia,ir,1) - & 02753 fac * e_norm_drho(ia,ir,1) * drho1b(idir)%array(ia,ir,1) 02754 END IF 02755 END IF 02756 END DO 02757 END DO 02758 END DO 02759 ELSE 02760 02761 ! partial integration 02762 DO idir=1, 3 02763 02764 DO ispin=1, nspins 02765 !$omp parallel do private(k,j,i,dr1dr) 02766 DO k = bo(1,3), bo(2,3) 02767 DO j = bo(1,2), bo(2,2) 02768 DO i = bo(1,1), bo(2,1) 02769 tmp_r(ispin)%pw%cr3d(i,j,k) = & 02770 tmp_a(ispin)%pw%cr3d(i,j,k) * drhoa(idir)%array(i,j,k) + & 02771 tmp_b(ispin)%pw%cr3d(i,j,k) * drhob(idir)%array(i,j,k) 02772 END DO 02773 END DO 02774 END DO 02775 IF (ASSOCIATED(e_drhoa)) THEN 02776 !$omp parallel do private(k,j,i) 02777 DO k = bo(1,3), bo(2,3) 02778 DO j = bo(1,2), bo(2,2) 02779 DO i = bo(1,1), bo(2,1) 02780 tmp_r(1)%pw%cr3d(i,j,k) = tmp_r(1)%pw%cr3d(i,j,k) - & 02781 e_drhoa(i,j,k) * drho1a(idir)%array(i,j,k) 02782 END DO 02783 END DO 02784 END DO 02785 END IF 02786 IF (nspins /= 1 .AND. ASSOCIATED(e_drhob)) THEN 02787 !$omp parallel do private(k,j,i) 02788 DO k = bo(1,3), bo(2,3) 02789 DO j = bo(1,2), bo(2,2) 02790 DO i = bo(1,1), bo(2,1) 02791 tmp_r(2)%pw%cr3d(i,j,k) = tmp_r(2)%pw%cr3d(i,j,k) - & 02792 e_drhob(i,j,k) * drho1b(idir)%array(i,j,k) 02793 END DO 02794 END DO 02795 END DO 02796 END IF 02797 IF (ASSOCIATED(e_norm_drho)) THEN 02798 !$omp parallel do private(k,j,i) 02799 DO k = bo(1,3), bo(2,3) 02800 DO j = bo(1,2), bo(2,2) 02801 DO i = bo(1,1), bo(2,1) 02802 IF (nspins /= 1) THEN 02803 tmp_r(1)%pw%cr3d(i,j,k) = tmp_r(1)%pw%cr3d(i,j,k) - & 02804 e_norm_drho(i,j,k) * drho1b(idir)%array(i,j,k) 02805 tmp_r(2)%pw%cr3d(i,j,k) = tmp_r(2)%pw%cr3d(i,j,k) - & 02806 e_norm_drho(i,j,k) * drho1a(idir)%array(i,j,k) 02807 ELSE 02808 tmp_r(1)%pw%cr3d(i,j,k) = tmp_r(1)%pw%cr3d(i,j,k) - & 02809 fac * e_norm_drho(i,j,k) * drho1b(idir)%array(i,j,k) 02810 END IF 02811 END DO 02812 END DO 02813 END DO 02814 END IF 02815 END DO 02816 02817 DO ispin=1, nspins 02818 SELECT CASE(xc_deriv_method_id) 02819 CASE (xc_deriv_spline2_smooth) 02820 CALL pw_nn_deriv_r ( pw_in=tmp_r(ispin)%pw,& 02821 pw_out=v_xc(ispin)%pw,coeffs=spline2_deriv_coeffs,& 02822 idir=idir, error=error ) 02823 CASE (xc_deriv_spline3_smooth) 02824 CALL pw_nn_deriv_r ( pw_in=tmp_r(ispin)%pw,& 02825 pw_out=v_xc(ispin)%pw,coeffs=spline3_deriv_coeffs,& 02826 idir=idir, error=error ) 02827 CASE (xc_deriv_nn10_smooth) 02828 CALL pw_nn_deriv_r ( pw_in=tmp_r(ispin)%pw,& 02829 pw_out=v_xc(ispin)%pw,coeffs=nn10_deriv_coeffs,& 02830 idir=idir, error=error ) 02831 CASE (xc_deriv_nn50_smooth) 02832 CALL pw_nn_deriv_r ( pw_in=tmp_r(ispin)%pw,& 02833 pw_out=v_xc(ispin)%pw,coeffs=nn50_deriv_coeffs,& 02834 idir=idir, error=error ) 02835 CASE default 02836 CALL cp_unimplemented_error(fromWhere=routineP, & 02837 message="XC_DERIV method not implemented for GPW-LSD!",& 02838 error=error, error_level=cp_failure_level) 02839 END SELECT 02840 END DO ! ispin 02841 02842 END DO ! idir 02843 02844 END IF 02845 02846 DO ispin=1, nspins 02847 IF (ASSOCIATED(pw_pool)) THEN 02848 CALL pw_pool_give_back_pw(pw_pool,tmp_a(ispin)%pw,error=error) 02849 CALL pw_pool_give_back_pw(pw_pool,tmp_b(ispin)%pw,error=error) 02850 CALL pw_pool_give_back_pw(pw_pool,tmp_r(ispin)%pw,error=error) 02851 ELSE 02852 DEALLOCATE(tmp_a(ispin)%pw%cr3d, tmp_b(ispin)%pw%cr3d, tmp_r(ispin)%pw%cr3d, stat=stat) 02853 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02854 DEALLOCATE(tmp_a(ispin)%pw, tmp_b(ispin)%pw, tmp_r(ispin)%pw, stat=stat) 02855 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02856 END IF 02857 END DO 02858 02859 END IF ! gradient_f 02860 02861 ELSE 02862 02863 !-----------------! 02864 ! restricted case ! 02865 !-----------------! 02866 02867 CALL xc_rho_set_get(rho1_set,rho=rho1, error=error) 02868 02869 IF (gradient_f) THEN 02870 CALL xc_rho_set_get(rho_set,drho=drho, error=error) 02871 CALL xc_rho_set_get(rho1_set,drho=drho1, error=error) 02872 IF (ASSOCIATED(pw_pool)) THEN 02873 CALL pw_pool_create_pw(pw_pool,v_drho%pw,& 02874 use_data = REALDATA3D,& 02875 in_space = REALSPACE, error=error) 02876 CALL pw_zero(v_drho%pw, error=error) 02877 ELSE 02878 ALLOCATE(v_drho%pw, stat=stat) 02879 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02880 ALLOCATE(v_drho%pw%cr3d(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3)), & 02881 stat=stat) 02882 v_drho%pw%cr3d = 0.0_dp 02883 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02884 END IF 02885 END IF 02886 02887 deriv_att=> xc_dset_get_derivative(deriv_set, "(rho)(rho)", & 02888 error=error) 02889 IF (ASSOCIATED(deriv_att)) THEN 02890 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02891 error=error) 02892 !$omp parallel do private(k,j,i) 02893 DO k = bo(1,3), bo(2,3) 02894 DO j = bo(1,2), bo(2,2) 02895 DO i = bo(1,1), bo(2,1) 02896 v_xc(1)%pw%cr3d(i,j,k)=& 02897 deriv_data(i,j,k)*rho1(i,j,k) 02898 END DO 02899 END DO 02900 END DO 02901 END IF 02902 02903 deriv_att=> xc_dset_get_derivative(deriv_set, "(rho)(norm_drho)", & 02904 error=error) 02905 IF (ASSOCIATED(deriv_att)) THEN 02906 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02907 error=error) 02908 !$omp parallel do private(k,j,i,idir,dr1dr) 02909 DO k = bo(1,3), bo(2,3) 02910 DO j = bo(1,2), bo(2,2) 02911 DO i = bo(1,1), bo(2,1) 02912 dr1dr=0._dp 02913 DO idir=1,3 02914 dr1dr=dr1dr+drho(idir)%array(i,j,k)*drho1(idir)%array(i,j,k) 02915 END DO 02916 v_xc(1)%pw%cr3d(i,j,k)=v_xc(1)%pw%cr3d(i,j,k)+& 02917 deriv_data(i,j,k)*dr1dr 02918 v_drho%pw%cr3d(i,j,k) = -1._dp * deriv_data(i,j,k)*rho1(i,j,k) 02919 END DO 02920 END DO 02921 END DO 02922 END IF 02923 02924 deriv_att=> xc_dset_get_derivative(deriv_set, "(norm_drho)(norm_drho)", & 02925 error=error) 02926 IF (ASSOCIATED(deriv_att)) THEN 02927 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02928 error=error) 02929 !$omp parallel do private(k,j,i,idir,dr1dr) 02930 DO k = bo(1,3), bo(2,3) 02931 DO j = bo(1,2), bo(2,2) 02932 DO i = bo(1,1), bo(2,1) 02933 dr1dr=0._dp 02934 DO idir=1,3 02935 dr1dr=dr1dr+drho(idir)%array(i,j,k)*drho1(idir)%array(i,j,k) 02936 END DO 02937 v_drho%pw%cr3d(i,j,k) = v_drho%pw%cr3d(i,j,k) - deriv_data(i,j,k)*dr1dr 02938 END DO 02939 END DO 02940 END DO 02941 END IF 02942 02943 deriv_att=> xc_dset_get_derivative(deriv_set, "(norm_drho)", & 02944 error=error) 02945 IF (ASSOCIATED(deriv_att)) THEN 02946 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data, & 02947 error=error) 02948 CALL xc_derivative_get(deriv_att, deriv_data=e_norm_drho, & 02949 error=error) 02950 !$omp parallel do private(k,j,i,idir,dr1dr) 02951 DO k = bo(1,3), bo(2,3) 02952 DO j = bo(1,2), bo(2,2) 02953 DO i = bo(1,1), bo(2,1) 02954 dr1dr=0._dp 02955 DO idir=1,3 02956 dr1dr=dr1dr+drho(idir)%array(i,j,k)*drho1(idir)%array(i,j,k) 02957 END DO 02958 IF (rho_set%norm_drho(i,j,k) > gradient_cut) THEN 02959 dr1dr = dr1dr / (rho_set%norm_drho(i,j,k))**2 02960 v_drho%pw%cr3d(i,j,k) = v_drho%pw%cr3d(i,j,k) + deriv_data(i,j,k)*dr1dr 02961 END IF 02962 END DO 02963 END DO 02964 END DO 02965 END IF 02966 02967 IF (gradient_f) THEN 02968 02969 IF (my_gapw) THEN 02970 02971 DO idir=1,3 02972 !$omp parallel do private(ia,ir) 02973 DO ia = bo(1,1), bo(2,1) 02974 DO ir = bo(1,2), bo(2,2) 02975 vxg(idir,ia,ir,1) = drho(idir)%array(ia,ir,1)*v_drho%pw%cr3d(ia,ir,1) 02976 IF (ASSOCIATED(e_norm_drho)) THEN 02977 vxg(idir,ia,ir,1) = vxg(idir,ia,ir,1) - drho1(idir)%array(ia,ir,1)*e_norm_drho(ia,ir,1) 02978 END IF 02979 END DO 02980 END DO 02981 END DO 02982 02983 ELSE 02984 ! partial integration 02985 02986 ! this does not work with non orthorombic cells 02987 ! (you will have to use a vector of pw with 3 components) 02988 IF (ASSOCIATED(pw_pool)) THEN 02989 CALL pw_pool_create_pw(pw_pool,tmp_r(1)%pw,& 02990 use_data = REALDATA3D,& 02991 in_space = REALSPACE, error=error) 02992 ELSE 02993 ALLOCATE(tmp_r(1)%pw, stat=stat) 02994 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02995 ALLOCATE(tmp_r(1)%pw%cr3d(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3)), & 02996 stat=stat) 02997 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02998 END IF 02999 03000 DO idir=1,3 03001 !$omp parallel do private(k,j,i) 03002 DO k = bo(1,3), bo(2,3) 03003 DO j = bo(1,2), bo(2,2) 03004 DO i = bo(1,1), bo(2,1) 03005 tmp_r(1)%pw%cr3d(i,j,k) = drho(idir)%array(i,j,k)*v_drho%pw%cr3d(i,j,k)-& 03006 drho1(idir)%array(i,j,k)*deriv_data(i,j,k) 03007 END DO 03008 END DO 03009 END DO 03010 03011 SELECT CASE(xc_deriv_method_id) 03012 CASE (xc_deriv_spline2_smooth) 03013 CALL pw_nn_deriv_r ( pw_in=tmp_r(1)%pw,& 03014 pw_out=v_xc(1)%pw,coeffs=spline2_deriv_coeffs,& 03015 idir=idir, error=error ) 03016 CASE (xc_deriv_spline3_smooth) 03017 CALL pw_nn_deriv_r ( pw_in=tmp_r(1)%pw,& 03018 pw_out=v_xc(1)%pw,coeffs=spline3_deriv_coeffs,& 03019 idir=idir, error=error ) 03020 CASE (xc_deriv_nn10_smooth) 03021 CALL pw_nn_deriv_r ( pw_in=tmp_r(1)%pw,& 03022 pw_out=v_xc(1)%pw,coeffs=nn10_deriv_coeffs,& 03023 idir=idir, error=error ) 03024 CASE (xc_deriv_nn50_smooth) 03025 CALL pw_nn_deriv_r ( pw_in=tmp_r(1)%pw,& 03026 pw_out=v_xc(1)%pw,coeffs=nn50_deriv_coeffs,& 03027 idir=idir, error=error ) 03028 CASE default 03029 CALL cp_unimplemented_error(fromWhere=routineP, & 03030 message="XC_DERIV method not implemented for GPW!",& 03031 error=error, error_level=cp_failure_level) 03032 END SELECT 03033 03034 END DO 03035 03036 IF (ASSOCIATED(pw_pool)) THEN 03037 CALL pw_pool_give_back_pw(pw_pool,tmp_r(1)%pw,error=error) 03038 ELSE 03039 DEALLOCATE(tmp_r(1)%pw%cr3d, stat=stat) 03040 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 03041 DEALLOCATE(tmp_r(1)%pw, stat=stat) 03042 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 03043 END IF 03044 END IF 03045 IF (ASSOCIATED(pw_pool)) THEN 03046 CALL pw_pool_give_back_pw(pw_pool,v_drho%pw,error=error) 03047 ELSE 03048 DEALLOCATE(v_drho%pw%cr3d, stat=stat) 03049 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 03050 DEALLOCATE(v_drho%pw, stat=stat) 03051 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 03052 END IF 03053 03054 END IF 03055 03056 END IF 03057 03058 DEALLOCATE(tmp_r, tmp_r2, tmp_a, tmp_b, stat=stat) 03059 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 03060 03061 CALL timestop(handle) 03062 END SUBROUTINE xc_calc_2nd_deriv 03063 03064 ! ***************************************************************************** 03083 SUBROUTINE xc_prep_2nd_deriv(deriv_set, & 03084 rho_set, rho_r, pw_pool, xc_section, cell, error) 03085 03086 TYPE(xc_derivative_set_type), POINTER :: deriv_set 03087 TYPE(xc_rho_set_type), POINTER :: rho_set 03088 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r 03089 TYPE(pw_pool_type), POINTER :: pw_pool 03090 TYPE(section_vals_type), POINTER :: xc_section 03091 TYPE(cell_type), POINTER :: cell 03092 TYPE(cp_error_type), INTENT(inout) :: error 03093 03094 CHARACTER(len=*), PARAMETER :: routineN = 'xc_prep_2nd_deriv', 03095 routineP = moduleN//':'//routineN 03096 03097 INTEGER :: handle, ispin, nspins, stat 03098 INTEGER, DIMENSION(2, 3) :: bo 03099 LOGICAL :: failure, lsd 03100 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g, rho_r_pw, tau 03101 03102 CALL timeset(routineN,handle) 03103 03104 failure=.FALSE. 03105 03106 CPPrecondition(.NOT.ASSOCIATED(rho_set),cp_failure_level,routineP,error,failure) 03107 CPPrecondition(ASSOCIATED(xc_section),cp_failure_level,routineP,error,failure) 03108 CPPrecondition(ASSOCIATED(pw_pool),cp_failure_level,routineP,error,failure) 03109 03110 IF (.NOT. failure) THEN 03111 nspins = SIZE(rho_r) 03112 lsd = (nspins /= 1) 03113 END IF 03114 03115 ALLOCATE(rho_r_pw(nspins), stat=stat) 03116 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 03117 DO ispin=1, nspins 03118 rho_r_pw(ispin)%pw => rho_r(ispin)%pw 03119 END DO 03120 03121 NULLIFY(rho_g, tau) 03122 CALL xc_rho_set_and_dset_create(rho_set,deriv_set,2,& 03123 rho_r_pw,rho_g,tau,xc_section,cell,pw_pool,& 03124 needs_basic_components=.TRUE.,error=error) 03125 03126 DEALLOCATE(rho_r_pw, stat=stat) 03127 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 03128 03129 bo = rho_r(1)%pw%pw_grid%bounds_local 03130 03131 CALL divide_by_norm_drho(deriv_set, rho_set, lsd, error) 03132 03133 CALL timestop(handle) 03134 END SUBROUTINE xc_prep_2nd_deriv 03135 03136 ! ***************************************************************************** 03137 SUBROUTINE divide_by_norm_drho(deriv_set, rho_set, lsd, error) 03138 03139 TYPE(xc_derivative_set_type), POINTER :: deriv_set 03140 TYPE(xc_rho_set_type), POINTER :: rho_set 03141 LOGICAL, INTENT(IN) :: lsd 03142 TYPE(cp_error_type), INTENT(inout) :: error 03143 03144 CHARACTER(len=*), PARAMETER :: routineN = 'divide_by_norm_drho', 03145 routineP = moduleN//':'//routineN 03146 03147 CHARACTER& 03148 (len=MAX_DERIVATIVE_DESC_LENGTH) :: desc 03149 CHARACTER(len=MAX_LABEL_LENGTH), 03150 DIMENSION(:), POINTER :: split_desc 03151 INTEGER :: i, idesc, j, k, order 03152 INTEGER, DIMENSION(2, 3) :: bo 03153 LOGICAL :: failure 03154 REAL(KIND=dp) :: drho_cutoff 03155 TYPE(cp_sll_xc_deriv_type), POINTER :: pos 03156 TYPE(xc_derivative_type), POINTER :: deriv_att 03157 03158 ! check for unknown derivatives and divide by norm_drho where necessary 03159 03160 failure = .FALSE. 03161 03162 bo = rho_set%local_bounds 03163 CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff, error=error) 03164 03165 pos => deriv_set%derivs 03166 DO WHILE (cp_sll_xc_deriv_next(pos,el_att=deriv_att,error=error)) 03167 CALL xc_derivative_get(deriv_att,order=order,& 03168 desc=desc, split_desc=split_desc,& 03169 error=error) 03170 IF (order==1 .OR. order==2) THEN 03171 DO idesc=1,SIZE(split_desc) 03172 SELECT CASE(split_desc(idesc)) 03173 CASE("norm_drho") 03174 !$omp parallel do private(i,j,k) 03175 DO k = bo(1,3), bo(2,3) 03176 DO j = bo(1,2), bo(2,2) 03177 DO i = bo(1,1), bo(2,1) 03178 deriv_att%deriv_data(i,j,k) = deriv_att%deriv_data(i,j,k) / & 03179 MAX(rho_set%norm_drho(i,j,k), drho_cutoff) 03180 END DO 03181 END DO 03182 END DO 03183 CASE("norm_drhoa") 03184 !$omp parallel do private(i,j,k) 03185 DO k = bo(1,3), bo(2,3) 03186 DO j = bo(1,2), bo(2,2) 03187 DO i = bo(1,1), bo(2,1) 03188 deriv_att%deriv_data(i,j,k) = deriv_att%deriv_data(i,j,k) / & 03189 MAX(rho_set%norm_drhoa(i,j,k), drho_cutoff) 03190 END DO 03191 END DO 03192 END DO 03193 CASE("norm_drhob") 03194 !$omp parallel do private(i,j,k) 03195 DO k = bo(1,3), bo(2,3) 03196 DO j = bo(1,2), bo(2,2) 03197 DO i = bo(1,1), bo(2,1) 03198 deriv_att%deriv_data(i,j,k) = deriv_att%deriv_data(i,j,k) / & 03199 MAX(rho_set%norm_drhob(i,j,k), drho_cutoff) 03200 END DO 03201 END DO 03202 END DO 03203 CASE("rho") 03204 CALL cp_assert(.NOT.lsd,cp_failure_level,cp_assertion_failed,routineP,& 03205 "rho not handled in lsd: '"//& 03206 TRIM(desc)//"' in "//& 03207 CPSourceFileRef,& 03208 error,failure) 03209 CASE("rhoa","rhob") 03210 CASE default 03211 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 03212 "unhandled derivative: '"//& 03213 TRIM(split_desc(idesc))//"' in '"//& 03214 TRIM(desc)//"' in "//& 03215 CPSourceFileRef,& 03216 error,failure) 03217 END SELECT 03218 END DO 03219 END IF 03220 END DO 03221 03222 END SUBROUTINE divide_by_norm_drho 03223 03224 ! ***************************************************************************** 03225 SUBROUTINE pw_smooth(pw_in,pw_out) 03226 TYPE(pw_type), POINTER :: pw_in, pw_out 03227 03228 INTEGER :: bo(2,3), i, il, ir, j, jl, 03229 jr, k, kl, kr, method, n(3), 03230 nc(3), p, q, r 03231 REAL(KIND=dp) :: alpha, beta, dist, dr(3), 03232 radius, sigma, sum 03233 REAL(KIND=dp), ALLOCATABLE, 03234 DIMENSION(:, :, :) :: Kernel 03235 03236 n(1:3) = pw_in%pw_grid%npts_local (1:3) 03237 dr(:) = pw_in%pw_grid%dr(:) 03238 bo = pw_in%pw_grid%bounds_local 03239 03240 method = 1 ! hard coded right now, like everything in here 03241 03242 SELECT CASE(method) 03243 CASE(1) ! just some averaging over neighbors, very fast 03244 alpha=1.0_dp 03245 beta =0.1_dp 03246 sum = alpha + 6*beta 03247 alpha = alpha/sum 03248 beta = beta/sum 03249 DO k = bo(1,3), bo(2,3) 03250 DO j = bo(1,2), bo(2,2) 03251 DO i = bo(1,1), bo(2,1) 03252 ir = MODULO(( i + 1 ) - bo(1,1),n(1))+bo(1,1) 03253 il = MODULO(( i - 1 ) - bo(1,1),n(1))+bo(1,1) 03254 jr = MODULO(( j + 1 ) - bo(1,2),n(2))+bo(1,2) 03255 jl = MODULO(( j - 1 ) - bo(1,2),n(2))+bo(1,2) 03256 kr = MODULO(( k + 1 ) - bo(1,3),n(3))+bo(1,3) 03257 kl = MODULO(( k - 1 ) - bo(1,3),n(3))+bo(1,3) 03258 pw_out%cr3d(i,j,k) = alpha*pw_in%cr3d(i,j,k)+beta*( & 03259 pw_in%cr3d(il,j,k)+pw_in%cr3d(ir,j,k)+ & 03260 pw_in%cr3d(i,jl,k)+pw_in%cr3d(i,jr,k)+ & 03261 pw_in%cr3d(i,j,kl)+pw_in%cr3d(i,j,kr)) 03262 END DO 03263 END DO 03264 END DO 03265 CASE(2) ! allowing for a more advanced functional form and wider mesh for averaging 03266 ! gets *very* slow rapidly. A g-space smoother would be possible 03267 ! however, this will most likely not be positive definite 03268 radius=0.5_dp 03269 sigma =0.1_dp 03270 nc(:)=CEILING(radius/dr(:)) 03271 ALLOCATE(Kernel(-nc(1):nc(1),-nc(2):nc(2),-nc(3):nc(3))) 03272 sum = 0.0_dp 03273 DO r=-nc(3),nc(3) 03274 DO q=-nc(2),nc(2) 03275 DO p=-nc(1),nc(1) 03276 dist=SQRT((r*dr(3))**2+(q*dr(2))**2+(p*dr(1))**2) 03277 Kernel(p,q,r)=EXP(-(dist/sigma)**2) 03278 sum = sum + Kernel(p,q,r) 03279 ENDDO 03280 ENDDO 03281 ENDDO 03282 ! normalize to 1 exactly. 03283 DO r=-nc(3),nc(3) 03284 DO q=-nc(2),nc(2) 03285 DO p=-nc(1),nc(1) 03286 Kernel(p,q,r)=Kernel(p,q,r)/sum 03287 ENDDO 03288 ENDDO 03289 ENDDO 03290 pw_out%cr3d(:,:,:) = 0.0_dp 03291 DO r=-nc(3),nc(3) 03292 DO q=-nc(2),nc(2) 03293 DO k = bo(1,3), bo(2,3) 03294 kr = MODULO(( k + r )- bo(1,3),n(3))+bo(1,3) 03295 DO j = bo(1,2), bo(2,2) 03296 jr = MODULO(( j + q )- bo(1,2),n(2))+bo(1,2) 03297 DO i = bo(1,1), bo(2,1) 03298 DO p=-nc(1),nc(1) 03299 ir = MODULO(( i + p )- bo(1,1),n(1))+bo(1,1) 03300 pw_out%cr3d(i,j,k) = pw_out%cr3d(i,j,k) + & 03301 Kernel(p,q,r)*pw_in%cr3d(ir,jr,kr) 03302 ENDDO 03303 END DO 03304 END DO 03305 END DO 03306 ENDDO 03307 ENDDO 03308 03309 DEALLOCATE(Kernel) 03310 END SELECT 03311 03312 END SUBROUTINE pw_smooth 03313 03314 END MODULE xc 03315
1.7.3