|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00012 MODULE qs_linres_epr_nablavks 00013 USE atomic_kind_types, ONLY: atomic_kind_type,& 00014 get_atomic_kind 00015 USE cell_types, ONLY: cell_type 00016 USE cp_control_types, ONLY: dft_control_type 00017 USE cp_output_handling, ONLY: cp_p_file,& 00018 cp_print_key_finished_output,& 00019 cp_print_key_should_output,& 00020 cp_print_key_unit_nr 00021 USE cp_para_types, ONLY: cp_para_env_type 00022 USE cp_subsys_types, ONLY: cp_subsys_get,& 00023 cp_subsys_type 00024 USE erf_fn, ONLY: erf,& 00025 erfc 00026 USE external_potential_types, ONLY: all_potential_type,& 00027 get_potential,& 00028 gth_potential_type 00029 USE f77_blas 00030 USE hartree_local_methods, ONLY: calculate_Vh_1center 00031 USE input_section_types, ONLY: section_get_ivals,& 00032 section_vals_get_subs_vals,& 00033 section_vals_type,& 00034 section_vals_val_get 00035 USE kinds, ONLY: dp 00036 USE mathconstants, ONLY: rootpi,& 00037 twopi 00038 USE particle_list_types, ONLY: particle_list_type 00039 USE particle_types, ONLY: particle_type 00040 USE pw_env_types, ONLY: pw_env_get,& 00041 pw_env_type 00042 USE pw_methods, ONLY: pw_axpy,& 00043 pw_copy,& 00044 pw_derive,& 00045 pw_transfer,& 00046 pw_zero 00047 USE pw_poisson_methods, ONLY: pw_poisson_solve 00048 USE pw_poisson_types, ONLY: pw_poisson_type 00049 USE pw_pool_types, ONLY: pw_pool_create_pw,& 00050 pw_pool_give_back_pw,& 00051 pw_pool_type 00052 USE pw_types, ONLY: COMPLEXDATA1D,& 00053 REALDATA3D,& 00054 REALSPACE,& 00055 RECIPROCALSPACE,& 00056 pw_p_type,& 00057 pw_type 00058 USE qs_environment_types, ONLY: get_qs_env,& 00059 qs_environment_type 00060 USE qs_gapw_densities, ONLY: prepare_gapw_den 00061 USE qs_grid_atom, ONLY: grid_atom_type 00062 USE qs_harmonics_atom, ONLY: harmonics_atom_type 00063 USE qs_ks_methods, ONLY: calc_rho_tot_gspace 00064 USE qs_linres_types, ONLY: epr_env_type,& 00065 get_epr_env,& 00066 nablavks_atom_type 00067 ! R0 00068 USE qs_local_rho_types, ONLY: rhoz_type 00069 USE qs_rho0_types, ONLY: rho0_atom_type 00070 USE qs_rho_atom_methods, ONLY: calculate_rho_atom_coeff 00071 USE qs_rho_atom_types, ONLY: get_rho_atom,& 00072 rho_atom_coeff,& 00073 rho_atom_type 00074 ! R0 00075 USE qs_rho_types, ONLY: qs_rho_p_type,& 00076 qs_rho_type 00077 USE qs_vxc, ONLY: qs_vxc_create 00078 USE qs_vxc_atom, ONLY: calculate_vxc_atom 00079 USE realspace_grid_cube, ONLY: pw_to_cube 00080 USE termination, ONLY: stop_memory 00081 USE util, ONLY: get_limit 00082 #include "cp_common_uses.h" 00083 00084 IMPLICIT NONE 00085 00086 PRIVATE 00087 PUBLIC :: epr_nablavks 00088 00089 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_epr_nablavks' 00090 00091 CONTAINS 00092 00093 ! ***************************************************************************** 00099 SUBROUTINE epr_nablavks(epr_env,qs_env,error) 00100 00101 TYPE(epr_env_type) :: epr_env 00102 TYPE(qs_environment_type), POINTER :: qs_env 00103 TYPE(cp_error_type), INTENT(INOUT), 00104 OPTIONAL :: error 00105 00106 CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_nablavks', 00107 routineP = moduleN//':'//routineN 00108 00109 CHARACTER(LEN=80) :: ext, filename 00110 COMPLEX(KIND=dp) :: gtemp 00111 INTEGER :: bo_atom(2), ia, iat, iatom, idir, iexp, ig, ikind, ir, iso, 00112 ispin, istat, ix, iy, iz, natom, nexp_ppl, nkind, nspins, output_unit, 00113 unit_nr 00114 INTEGER, DIMENSION(2, 3) :: bo 00115 INTEGER, DIMENSION(:), POINTER :: atom_list 00116 LOGICAL :: failure, gapw, gapw_xc, 00117 gth_gspace, ionode, 00118 make_soft, paw_atom 00119 REAL(KIND=dp) :: alpha, alpha_core, arg, charge, ehartree, exc, exp_rap, 00120 gapw_max_alpha, hard_value, soft_value, sqrt_alpha, sqrt_rap 00121 REAL(KIND=dp), ALLOCATABLE, 00122 DIMENSION(:, :) :: vh1_rad_h, vh1_rad_s 00123 REAL(KIND=dp), DIMENSION(3) :: rap, ratom, roffset, rpoint 00124 REAL(KIND=dp), DIMENSION(:), POINTER :: cexp_ppl, rho_rad_z 00125 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rho_rad_0 00126 TYPE(all_potential_type), POINTER :: all_potential 00127 TYPE(atomic_kind_type), DIMENSION(:), 00128 POINTER :: atomic_kind_set 00129 TYPE(atomic_kind_type), POINTER :: atomic_kind 00130 TYPE(cell_type), POINTER :: cell 00131 TYPE(cp_logger_type), POINTER :: logger 00132 TYPE(cp_para_env_type), POINTER :: para_env 00133 TYPE(cp_subsys_type), POINTER :: subsys 00134 TYPE(dft_control_type), POINTER :: dft_control 00135 TYPE(grid_atom_type), POINTER :: grid_atom 00136 TYPE(gth_potential_type), POINTER :: gth_potential 00137 TYPE(harmonics_atom_type), POINTER :: harmonics 00138 TYPE(nablavks_atom_type), DIMENSION(:), 00139 POINTER :: nablavks_atom_set 00140 TYPE(particle_list_type), POINTER :: particles 00141 TYPE(particle_type), DIMENSION(:), 00142 POINTER :: particle_set 00143 TYPE(pw_env_type), POINTER :: pw_env 00144 TYPE(pw_p_type), DIMENSION(:), POINTER :: v_rspace_new, v_tau_rspace 00145 TYPE(pw_p_type), POINTER :: rho_tot_gspace, v_coulomb_gspace, 00146 v_coulomb_gtemp, v_coulomb_rtemp, v_hartree_gspace, v_hartree_gtemp, 00147 v_hartree_rtemp, v_xc_gtemp, v_xc_rtemp, wf_r 00148 TYPE(pw_poisson_type), POINTER :: poisson_env 00149 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00150 TYPE(pw_type), POINTER :: pwx, pwy, pwz 00151 TYPE(qs_rho_p_type), DIMENSION(:, :), 00152 POINTER :: nablavks_set 00153 TYPE(qs_rho_type), POINTER :: rho 00154 TYPE(rho0_atom_type), DIMENSION(:), 00155 POINTER :: rho0_atom_set 00156 TYPE(rho_atom_coeff), DIMENSION(:), 00157 POINTER :: rho_rad_h, rho_rad_s 00158 TYPE(rho_atom_coeff), DIMENSION(:, :), 00159 POINTER :: nablavks_vec_rad_h, 00160 nablavks_vec_rad_s 00161 TYPE(rho_atom_type), DIMENSION(:), 00162 POINTER :: rho_atom_set 00163 TYPE(rho_atom_type), POINTER :: rho_atom 00164 TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set 00165 TYPE(section_vals_type), POINTER :: g_section, input, lr_section, 00166 xc_section 00167 00168 ! R0 00169 ! R0 00170 ! R0 00171 ! R0 00172 ! R0 00173 ! R0 00174 00175 NULLIFY(auxbas_pw_pool) 00176 NULLIFY(cell) 00177 NULLIFY(dft_control) 00178 NULLIFY(g_section) 00179 NULLIFY(logger) 00180 NULLIFY(lr_section) 00181 NULLIFY(nablavks_set) 00182 NULLIFY(nablavks_atom_set) ! R0 00183 NULLIFY(nablavks_vec_rad_h) ! R0 00184 NULLIFY(nablavks_vec_rad_s) ! R0 00185 NULLIFY(para_env) 00186 NULLIFY(particle_set) 00187 NULLIFY(particles) 00188 NULLIFY(poisson_env) 00189 NULLIFY(pw_env) 00190 NULLIFY(pwx) ! R0 00191 NULLIFY(pwy) ! R0 00192 NULLIFY(pwz) ! R0 00193 NULLIFY(rho) 00194 NULLIFY(rho0_atom_set) 00195 NULLIFY(rho_atom_set) 00196 NULLIFY(rhoz_set) 00197 NULLIFY(subsys) 00198 NULLIFY(v_rspace_new) 00199 NULLIFY(v_tau_rspace) 00200 NULLIFY(xc_section) 00201 NULLIFY(input) 00202 00203 logger => cp_error_get_logger(error) 00204 lr_section => section_vals_get_subs_vals(qs_env%input,"PROPERTIES%LINRES",error=error) 00205 ionode = logger%para_env%mepos==logger%para_env%source 00206 00207 output_unit = cp_print_key_unit_nr(logger,lr_section,"PRINT%PROGRAM_RUN_INFO",& 00208 extension=".linresLog",error=error) 00209 00210 failure = .FALSE. 00211 00212 ! ------------------------------------- 00213 ! Read settings 00214 ! ------------------------------------- 00215 00216 g_section => section_vals_get_subs_vals(lr_section,& 00217 "EPR%PRINT%G_TENSOR",error=error) 00218 00219 CALL section_vals_val_get(g_section,"gapw_max_alpha",r_val=gapw_max_alpha, error=error) 00220 00221 gth_gspace = .TRUE. 00222 00223 ! ------------------------------------- 00224 ! Get nablavks arrays 00225 ! ------------------------------------- 00226 00227 CALL get_epr_env(epr_env,nablavks_set=nablavks_set,& ! R0 00228 nablavks_atom_set=nablavks_atom_set,& ! R0 00229 error=error) ! R0 00230 00231 DO ispin = 1,SIZE(nablavks_set,2) 00232 DO idir = 1,SIZE(nablavks_set,1) 00233 CALL pw_zero(nablavks_set(idir,ispin)%rho%rho_r(1)%pw, error=error) 00234 END DO 00235 END DO 00236 00237 pwx => nablavks_set(1,1)%rho%rho_r(1)%pw 00238 pwy => nablavks_set(2,1)%rho%rho_r(1)%pw 00239 pwz => nablavks_set(3,1)%rho%rho_r(1)%pw 00240 00241 roffset = -REAL(MODULO(pwx%pw_grid%npts,2),dp)*pwx%pw_grid%dr/2.0_dp 00242 00243 ! ------------------------------------- 00244 ! Get grids / atom info 00245 ! ------------------------------------- 00246 00247 CALL get_qs_env(qs_env=qs_env,& 00248 atomic_kind_set=atomic_kind_set,& 00249 input=input,& 00250 cell=cell,& 00251 dft_control=dft_control,& 00252 para_env=para_env,& 00253 particle_set=particle_set,& 00254 pw_env=pw_env,& 00255 rho=rho,& 00256 rho_atom_set=rho_atom_set,& 00257 rho0_atom_set=rho0_atom_set,& 00258 rhoz_set=rhoz_set,& 00259 subsys=subsys,& 00260 error=error) 00261 00262 CALL pw_env_get(pw_env,auxbas_pw_pool=auxbas_pw_pool,& 00263 poisson_env=poisson_env,& 00264 error=error) 00265 00266 CALL cp_subsys_get(subsys,particles=particles,error=error) 00267 00268 gapw = dft_control%qs_control%gapw 00269 gapw_xc = dft_control%qs_control%gapw_xc 00270 nkind = SIZE(atomic_kind_set) 00271 nspins = dft_control%nspins 00272 00273 ! ------------------------------------- 00274 ! Add Hartree potential 00275 ! ------------------------------------- 00276 00277 ALLOCATE(v_hartree_gspace,STAT=istat) 00278 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00279 "v_hartree_gspace",0) 00280 ALLOCATE(v_hartree_gtemp,STAT=istat) 00281 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00282 "v_hartree_gtemp",0) 00283 ALLOCATE(v_hartree_rtemp,STAT=istat) 00284 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00285 "v_hartree_rtemp",0) 00286 ALLOCATE(rho_tot_gspace,STAT=istat) 00287 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00288 "rho_tot_gspace",0) 00289 00290 CALL pw_pool_create_pw(auxbas_pw_pool,v_hartree_gspace%pw, & 00291 use_data=COMPLEXDATA1D,in_space=RECIPROCALSPACE,& 00292 error=error) 00293 CALL pw_pool_create_pw(auxbas_pw_pool,v_hartree_gtemp%pw, & 00294 use_data=COMPLEXDATA1D,in_space=RECIPROCALSPACE,& 00295 error=error) 00296 CALL pw_pool_create_pw(auxbas_pw_pool,v_hartree_rtemp%pw,& 00297 use_data=REALDATA3D,in_space=REALSPACE,& 00298 error=error) 00299 CALL pw_pool_create_pw(auxbas_pw_pool,rho_tot_gspace%pw,& 00300 use_data=COMPLEXDATA1D,in_space=RECIPROCALSPACE,& 00301 error=error) 00302 00303 IF (gapw) THEN 00304 ! need to rebuild the coeff ! 00305 CALL calculate_rho_atom_coeff(qs_env,rho%rho_ao,error=error) 00306 CALL prepare_gapw_den(qs_env,do_rho0=.TRUE.,error=error) 00307 END IF 00308 00309 CALL pw_zero(rho_tot_gspace%pw, error=error) 00310 00311 CALL calc_rho_tot_gspace(rho_tot_gspace,qs_env,rho,& 00312 skip_nuclear_density=.NOT. gapw,error=error) 00313 00314 CALL pw_poisson_solve(poisson_env,rho_tot_gspace%pw,ehartree,& 00315 v_hartree_gspace%pw,error=error) 00316 00317 ! ------------------------------------- 00318 ! Atomic grids part 00319 ! ------------------------------------- 00320 00321 IF (gapw) THEN 00322 00323 DO ikind = 1,nkind ! loop over atom types 00324 00325 NULLIFY(atomic_kind) 00326 NULLIFY(atom_list) 00327 NULLIFY(grid_atom) 00328 NULLIFY(harmonics) 00329 NULLIFY(rho_rad_z) 00330 00331 atomic_kind => atomic_kind_set(ikind) 00332 rho_rad_z => rhoz_set(ikind)%r_coef 00333 00334 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00335 atom_list=atom_list,& 00336 grid_atom=grid_atom,& 00337 harmonics=harmonics,& 00338 natom=natom,& 00339 paw_atom=paw_atom,& 00340 zeff=charge,& 00341 alpha_core_charge=alpha_core) 00342 00343 IF (paw_atom) THEN 00344 00345 ALLOCATE(vh1_rad_h(grid_atom%nr,harmonics%max_iso_not0),STAT=istat) 00346 ALLOCATE(vh1_rad_s(grid_atom%nr,harmonics%max_iso_not0),STAT=istat) 00347 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00348 00349 ! DO iat = 1, natom ! natom = # atoms for ikind 00350 ! 00351 ! iatom = atom_list(iat) 00352 ! ratom = particle_set(iatom)%r 00353 ! 00354 ! DO ig = v_hartree_gspace%pw%pw_grid%first_gne0,v_hartree_gspace%pw%pw_grid%ngpts_cut_local 00355 ! 00356 ! gtemp = fourpi * charge / cell%deth * & 00357 ! EXP ( - v_hartree_gspace%pw%pw_grid%gsq(ig) / (4.0_dp * alpha_core) ) & 00358 ! / v_hartree_gspace%pw%pw_grid%gsq(ig) 00359 ! 00360 ! arg = DOT_PRODUCT(v_hartree_gspace%pw%pw_grid%g(:,ig),ratom) 00361 ! 00362 ! gtemp = gtemp * CMPLX(COS(arg),-SIN(arg),KIND=dp) 00363 ! 00364 ! v_hartree_gspace%pw%cc(ig) = v_hartree_gspace%pw%cc(ig) + gtemp 00365 ! END DO 00366 ! IF ( v_hartree_gspace%pw%pw_grid%have_g0 ) v_hartree_gspace%pw%cc(1) = 0.0_dp 00367 ! 00368 ! END DO 00369 00370 bo_atom = get_limit(natom, para_env%num_pe, para_env%mepos) 00371 00372 DO iat = bo_atom(1),bo_atom(2) ! natomkind = # atoms for ikind 00373 00374 iatom = atom_list(iat) 00375 ratom = particle_set(iatom)%r 00376 00377 nablavks_vec_rad_h => nablavks_atom_set(iatom)%nablavks_vec_rad_h 00378 nablavks_vec_rad_s => nablavks_atom_set(iatom)%nablavks_vec_rad_s 00379 00380 DO ispin = 1,nspins 00381 DO idir = 1,3 00382 nablavks_vec_rad_h(idir,ispin)%r_coef(:,:) = 0.0_dp 00383 nablavks_vec_rad_s(idir,ispin)%r_coef(:,:) = 0.0_dp 00384 END DO ! idir 00385 END DO ! ispin 00386 00387 rho_atom => rho_atom_set(iatom) 00388 NULLIFY(rho_rad_h,rho_rad_s,rho_rad_0) 00389 CALL get_rho_atom(rho_atom=rho_atom,rho_rad_h=rho_rad_h,& 00390 rho_rad_s=rho_rad_s) 00391 rho_rad_0 => rho0_atom_set(iatom)%rho0_rad_h%r_coef 00392 vh1_rad_h = 0.0_dp 00393 vh1_rad_s = 0.0_dp 00394 00395 CALL calculate_Vh_1center(vh1_rad_h,vh1_rad_s,rho_rad_h,rho_rad_s,rho_rad_0,rho_rad_z,grid_atom) 00396 00397 DO ir = 2,grid_atom%nr 00398 00399 IF (grid_atom%rad(ir) >= atomic_kind%hard_radius) CYCLE 00400 00401 DO ia = 1,grid_atom%ng_sphere 00402 00403 ! hard part 00404 00405 DO idir= 1,3 00406 hard_value = 0.0_dp 00407 DO iso = 1,harmonics%max_iso_not0 00408 hard_value = hard_value + & 00409 vh1_rad_h(ir,iso) * harmonics%dslm_dxyz(idir,ia,iso) + & 00410 harmonics%slm(ia,iso) * & 00411 ( vh1_rad_h(ir - 1,iso) - vh1_rad_h(ir,iso) ) / & 00412 ( grid_atom%rad(ir - 1) - grid_atom%rad(ir) ) * & 00413 ( harmonics%a(idir,ia) ) 00414 END DO 00415 nablavks_vec_rad_h(idir,1)%r_coef(ir,ia) = hard_value 00416 END DO 00417 00418 ! soft part 00419 00420 DO idir= 1,3 00421 soft_value = 0.0_dp 00422 DO iso = 1,harmonics%max_iso_not0 00423 soft_value = soft_value + & 00424 vh1_rad_s(ir,iso) * harmonics%dslm_dxyz(idir,ia,iso) + & 00425 harmonics%slm(ia,iso) * & 00426 ( vh1_rad_s(ir - 1,iso) - vh1_rad_s(ir,iso) ) / & 00427 ( grid_atom%rad(ir - 1) - grid_atom%rad(ir) ) * & 00428 ( harmonics%a(idir,ia) ) 00429 END DO 00430 nablavks_vec_rad_s(idir,1)%r_coef(ir,ia) = soft_value 00431 END DO 00432 00433 END DO ! ia 00434 00435 END DO ! ir 00436 00437 DO idir = 1,3 00438 nablavks_vec_rad_h(idir,2)%r_coef(:,:) = nablavks_vec_rad_h(idir,1)%r_coef(:,:) 00439 nablavks_vec_rad_s(idir,2)%r_coef(:,:) = nablavks_vec_rad_s(idir,1)%r_coef(:,:) 00440 END DO 00441 00442 END DO ! iat 00443 00444 DEALLOCATE(vh1_rad_h,STAT=istat) 00445 DEALLOCATE(vh1_rad_s,STAT=istat) 00446 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00447 00448 END IF ! paw_atom 00449 00450 END DO ! ikind 00451 00452 END IF ! gapw 00453 00454 CALL pw_copy(v_hartree_gspace%pw, v_hartree_gtemp%pw, error=error) 00455 CALL pw_derive(v_hartree_gtemp%pw, (/1,0,0/) , error=error) 00456 CALL pw_transfer(v_hartree_gtemp%pw, v_hartree_rtemp%pw, error=error) 00457 CALL pw_copy(v_hartree_rtemp%pw, pwx, error=error) 00458 00459 CALL pw_copy(v_hartree_gspace%pw, v_hartree_gtemp%pw, error=error) 00460 CALL pw_derive(v_hartree_gtemp%pw, (/0,1,0/) , error=error) 00461 CALL pw_transfer(v_hartree_gtemp%pw, v_hartree_rtemp%pw, error=error) 00462 CALL pw_copy(v_hartree_rtemp%pw, pwy, error=error) 00463 00464 CALL pw_copy(v_hartree_gspace%pw, v_hartree_gtemp%pw, error=error) 00465 CALL pw_derive(v_hartree_gtemp%pw, (/0,0,1/) , error=error) 00466 CALL pw_transfer(v_hartree_gtemp%pw, v_hartree_rtemp%pw, error=error) 00467 CALL pw_copy(v_hartree_rtemp%pw, pwz, error=error) 00468 00469 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_hartree_gspace%pw,error=error) 00470 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_hartree_gtemp%pw,error=error) 00471 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_hartree_rtemp%pw,error=error) 00472 CALL pw_pool_give_back_pw(auxbas_pw_pool,rho_tot_gspace%pw,error=error) 00473 00474 DEALLOCATE(v_hartree_gspace,STAT=istat) 00475 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00476 "v_hartree_gspace") 00477 DEALLOCATE(v_hartree_gtemp,STAT=istat) 00478 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00479 "v_hartree_gtemp") 00480 DEALLOCATE(v_hartree_rtemp,STAT=istat) 00481 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00482 "v_hartree_rtemp") 00483 DEALLOCATE(rho_tot_gspace,STAT=istat) 00484 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00485 "rho_tot_gspace") 00486 00487 ! ------------------------------------- 00488 ! Add Coulomb potential 00489 ! ------------------------------------- 00490 00491 DO ikind = 1,nkind ! loop over atom types 00492 00493 NULLIFY(atom_list) 00494 NULLIFY(atomic_kind) 00495 NULLIFY(grid_atom) 00496 NULLIFY(harmonics) 00497 00498 atomic_kind => atomic_kind_set(ikind) 00499 00500 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00501 all_potential=all_potential,& 00502 atom_list=atom_list,& 00503 grid_atom=grid_atom,& 00504 gth_potential=gth_potential,& 00505 harmonics=harmonics,& 00506 natom=natom,& 00507 paw_atom=paw_atom) 00508 00509 IF (ASSOCIATED(gth_potential)) THEN 00510 00511 NULLIFY(cexp_ppl) 00512 00513 CALL get_potential(potential=gth_potential,& 00514 zeff=charge,& 00515 alpha_ppl=alpha,& 00516 nexp_ppl=nexp_ppl,& 00517 cexp_ppl=cexp_ppl) 00518 00519 sqrt_alpha = SQRT(alpha) 00520 00521 IF (gapw .AND. paw_atom .AND. alpha > gapw_max_alpha) THEN 00522 make_soft = .TRUE. 00523 ELSE 00524 make_soft = .FALSE. 00525 END IF 00526 00527 ! ------------------------------------- 00528 ! PW grid part 00529 ! ------------------------------------- 00530 00531 IF (gth_gspace) THEN 00532 00533 ALLOCATE(v_coulomb_gspace,STAT=istat) 00534 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00535 "v_coulomb_gspace",0) 00536 ALLOCATE(v_coulomb_gtemp,STAT=istat) 00537 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00538 "v_coulomb_gtemp",0) 00539 ALLOCATE(v_coulomb_rtemp,STAT=istat) 00540 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00541 "v_coulomb_rtemp",0) 00542 00543 CALL pw_pool_create_pw(auxbas_pw_pool,v_coulomb_gspace%pw, & 00544 use_data=COMPLEXDATA1D,in_space=RECIPROCALSPACE,& 00545 error=error) 00546 CALL pw_pool_create_pw(auxbas_pw_pool,v_coulomb_gtemp%pw, & 00547 use_data=COMPLEXDATA1D,in_space=RECIPROCALSPACE,& 00548 error=error) 00549 CALL pw_pool_create_pw(auxbas_pw_pool,v_coulomb_rtemp%pw,& 00550 use_data=REALDATA3D,in_space=REALSPACE,& 00551 error=error) 00552 00553 CALL pw_zero(v_coulomb_gspace%pw, error=error) 00554 00555 DO iat = 1, natom ! natom = # atoms for ikind 00556 00557 iatom = atom_list(iat) 00558 ratom = particle_set(iatom)%r 00559 00560 DO ig = v_coulomb_gspace%pw%pw_grid%first_gne0,v_coulomb_gspace%pw%pw_grid%ngpts_cut_local 00561 gtemp = 0.0_dp 00562 ! gtemp = - fourpi * charge / cell%deth * & 00563 ! EXP ( - v_coulomb_gspace%pw%pw_grid%gsq(ig) / (4.0_dp * alpha) ) & 00564 ! / v_coulomb_gspace%pw%pw_grid%gsq(ig) 00565 00566 IF (.NOT. make_soft) THEN 00567 00568 SELECT CASE (nexp_ppl) 00569 CASE(1) 00570 gtemp = gtemp + & 00571 (twopi)**(1.5_dp) / ( cell%deth * (2.0_dp * alpha)**(1.5_dp) ) * & 00572 EXP ( - v_coulomb_gspace%pw%pw_grid%gsq(ig) / (4.0_dp * alpha) ) * ( & 00573 ! C1 00574 + cexp_ppl(1) & 00575 ) 00576 CASE(2) 00577 gtemp = gtemp + & 00578 (twopi)**(1.5_dp) / ( cell%deth * (2.0_dp * alpha)**(1.5_dp) ) * & 00579 EXP ( - v_coulomb_gspace%pw%pw_grid%gsq(ig) / (4.0_dp * alpha) ) * ( & 00580 ! C1 00581 + cexp_ppl(1) & 00582 ! C2 00583 + cexp_ppl(2) / (2.0_dp * alpha) * & 00584 ( 3.0_dp - v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp * alpha) ) & 00585 ) 00586 CASE(3) 00587 gtemp = gtemp + & 00588 (twopi)**(1.5_dp) / ( cell%deth * (2.0_dp * alpha)**(1.5_dp) ) * & 00589 EXP ( - v_coulomb_gspace%pw%pw_grid%gsq(ig) / (4.0_dp * alpha) ) * ( & 00590 ! C1 00591 + cexp_ppl(1) & 00592 ! C2 00593 + cexp_ppl(2) / (2.0_dp * alpha) * & 00594 ( 3.0_dp - v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp * alpha) ) & 00595 ! C3 00596 + cexp_ppl(3) / (2.0_dp * alpha)**2 * & 00597 ( 15.0_dp - 10.0_dp * v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp * alpha) & 00598 + ( v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp * alpha) )**2 ) & 00599 ) 00600 CASE(4) 00601 gtemp = gtemp + & 00602 (twopi)**(1.5_dp) / ( cell%deth * (2.0_dp * alpha)**(1.5_dp) ) * & 00603 EXP ( - v_coulomb_gspace%pw%pw_grid%gsq(ig) / (4.0_dp * alpha) ) * ( & 00604 ! C1 00605 + cexp_ppl(1) & 00606 ! C2 00607 + cexp_ppl(2) / (2.0_dp * alpha) * & 00608 ( 3.0_dp - v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp * alpha) ) & 00609 ! C3 00610 + cexp_ppl(3) / (2.0_dp * alpha)**2 * & 00611 ( 15.0_dp - 10.0_dp * v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp * alpha) & 00612 + ( v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp * alpha) )**2 ) & 00613 ! C4 00614 + cexp_ppl(4) / (2.0_dp * alpha)**3 * & 00615 ( 105.0_dp - 105.0_dp * v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp * alpha) & 00616 + 21.0_dp * ( v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp * alpha) )**2 & 00617 - ( v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp * alpha) )**3 ) & 00618 ) 00619 END SELECT 00620 00621 END IF 00622 00623 arg = DOT_PRODUCT(v_coulomb_gspace%pw%pw_grid%g(:,ig),ratom) 00624 00625 gtemp = gtemp * CMPLX(COS(arg),-SIN(arg),KIND=dp) 00626 v_coulomb_gspace%pw%cc(ig) = v_coulomb_gspace%pw%cc(ig) + gtemp 00627 END DO 00628 IF ( v_coulomb_gspace%pw%pw_grid%have_g0 ) v_coulomb_gspace%pw%cc(1) = 0.0_dp 00629 00630 END DO 00631 00632 CALL pw_copy(v_coulomb_gspace%pw, v_coulomb_gtemp%pw,error=error) 00633 CALL pw_derive(v_coulomb_gtemp%pw, (/1,0,0/),error=error ) 00634 CALL pw_transfer(v_coulomb_gtemp%pw, v_coulomb_rtemp%pw,error=error) 00635 CALL pw_axpy(v_coulomb_rtemp%pw, pwx,error=error) 00636 00637 CALL pw_copy(v_coulomb_gspace%pw, v_coulomb_gtemp%pw,error=error) 00638 CALL pw_derive(v_coulomb_gtemp%pw, (/0,1,0/),error=error ) 00639 CALL pw_transfer(v_coulomb_gtemp%pw, v_coulomb_rtemp%pw,error=error) 00640 CALL pw_axpy(v_coulomb_rtemp%pw, pwy,error=error) 00641 00642 CALL pw_copy(v_coulomb_gspace%pw, v_coulomb_gtemp%pw,error=error) 00643 CALL pw_derive(v_coulomb_gtemp%pw, (/0,0,1/),error=error ) 00644 CALL pw_transfer(v_coulomb_gtemp%pw, v_coulomb_rtemp%pw,error=error) 00645 CALL pw_axpy(v_coulomb_rtemp%pw, pwz,error=error) 00646 00647 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_coulomb_gspace%pw,error=error) 00648 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_coulomb_gtemp%pw,error=error) 00649 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_coulomb_rtemp%pw,error=error) 00650 00651 DEALLOCATE(v_coulomb_gspace,STAT=istat) 00652 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00653 "v_coulomb_gspace") 00654 DEALLOCATE(v_coulomb_gtemp,STAT=istat) 00655 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00656 "v_coulomb_gtemp") 00657 DEALLOCATE(v_coulomb_rtemp,STAT=istat) 00658 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00659 "v_coulomb_rtemp") 00660 ELSE 00661 00662 ! Attic of the atomic parallellisation 00663 ! 00664 ! bo(2) 00665 ! bo = get_limit(natom, para_env%num_pe, para_env%mepos) 00666 ! DO iat = bo(1),bo(2) ! natom = # atoms for ikind 00667 ! DO ix = lbound(pwx%cr3d,1), ubound(pwx%cr3d,1) 00668 ! DO iy = lbound(pwx%cr3d,2), ubound(pwx%cr3d,2) 00669 ! DO iz = lbound(pwx%cr3d,3), ubound(pwx%cr3d,3) 00670 00671 bo = pwx%pw_grid%bounds_local 00672 00673 DO iat = 1, natom ! natom = # atoms for ikind 00674 00675 iatom = atom_list(iat) 00676 ratom = particle_set(iatom)%r 00677 00678 DO ix = bo(1,1),bo(2,1) 00679 DO iy = bo(1,2),bo(2,2) 00680 DO iz = bo(1,3),bo(2,3) 00681 rpoint = (/REAL(ix,dp)*pwx%pw_grid%dr(1), 00682 REAL(iy,dp)*pwx%pw_grid%dr(2), 00683 REAL(iz,dp)*pwx%pw_grid%dr(3)/) 00684 rpoint = rpoint + roffset 00685 rap = rpoint - ratom 00686 rap(1)=MODULO(rap(1),cell%hmat(1,1))-cell%hmat(1,1)/2._dp 00687 rap(2)=MODULO(rap(2),cell%hmat(2,2))-cell%hmat(2,2)/2._dp 00688 rap(3)=MODULO(rap(3),cell%hmat(3,3))-cell%hmat(3,3)/2._dp 00689 sqrt_rap = SQRT(DOT_PRODUCT(rap,rap)) 00690 exp_rap = EXP( - alpha * sqrt_rap**2 ) 00691 IF (sqrt_rap < 1.e-10_dp ) sqrt_rap = 1.e-10_dp 00692 ! d_x 00693 00694 pwx%cr3d(ix,iy,iz) = pwx%cr3d(ix,iy,iz) + charge * ( & 00695 - 2.0_dp * sqrt_alpha * EXP( - sqrt_rap**2 * sqrt_alpha**2 ) * rap(1) & 00696 / ( rootpi * sqrt_rap**2 ) & 00697 + erf( sqrt_rap * sqrt_alpha ) * rap(1) & 00698 / sqrt_rap**3 ) 00699 00700 ! d_y 00701 00702 pwy%cr3d(ix,iy,iz) = pwy%cr3d(ix,iy,iz) + charge * ( & 00703 - 2.0_dp * sqrt_alpha * EXP( - sqrt_rap**2 * sqrt_alpha**2 ) * rap(2) & 00704 / ( rootpi * sqrt_rap**2 ) & 00705 + erf( sqrt_rap * sqrt_alpha ) * rap(2) & 00706 / sqrt_rap**3 ) 00707 00708 ! d_z 00709 00710 pwz%cr3d(ix,iy,iz) = pwz%cr3d(ix,iy,iz) + charge * ( & 00711 - 2.0_dp * sqrt_alpha * EXP( - sqrt_rap**2 * sqrt_alpha**2 ) * rap(3) & 00712 / ( rootpi * sqrt_rap**2 ) & 00713 + erf( sqrt_rap * sqrt_alpha ) * rap(3) & 00714 / sqrt_rap**3 ) 00715 00716 IF (make_soft) CYCLE 00717 00718 ! d_x 00719 00720 DO iexp = 1,nexp_ppl 00721 pwx%cr3d(ix,iy,iz) = pwx%cr3d(ix,iy,iz) + ( & 00722 - 2.0_dp * alpha * rap(1) * exp_rap * & 00723 cexp_ppl(iexp) * ( sqrt_rap**2 )**(iexp - 1) ) 00724 IF (iexp > 1) THEN 00725 pwx%cr3d(ix,iy,iz) = pwx%cr3d(ix,iy,iz) + ( & 00726 2.0_dp * exp_rap * cexp_ppl(iexp) * & 00727 ( sqrt_rap**2 )**(iexp - 2) * REAL(iexp - 1,dp) * rap(1) ) 00728 END IF 00729 END DO 00730 00731 ! d_y 00732 00733 DO iexp = 1,nexp_ppl 00734 pwy%cr3d(ix,iy,iz) = pwy%cr3d(ix,iy,iz) + ( & 00735 - 2.0_dp * alpha * rap(2) * exp_rap * & 00736 cexp_ppl(iexp) * ( sqrt_rap**2 )**(iexp - 1) ) 00737 IF (iexp > 1) THEN 00738 pwy%cr3d(ix,iy,iz) = pwy%cr3d(ix,iy,iz) + ( & 00739 2.0_dp * exp_rap * cexp_ppl(iexp) * & 00740 ( sqrt_rap**2 )**(iexp - 2) * REAL(iexp - 1,dp) * rap(2) ) 00741 END IF 00742 END DO 00743 00744 ! d_z 00745 00746 DO iexp = 1,nexp_ppl 00747 pwz%cr3d(ix,iy,iz) = pwz%cr3d(ix,iy,iz) + ( & 00748 - 2.0_dp * alpha * rap(3) * exp_rap * & 00749 cexp_ppl(iexp) * ( sqrt_rap**2 )**(iexp - 1) ) 00750 IF (iexp > 1) THEN 00751 pwz%cr3d(ix,iy,iz) = pwz%cr3d(ix,iy,iz) + ( & 00752 2.0_dp * exp_rap * cexp_ppl(iexp) * & 00753 ( sqrt_rap**2 )**(iexp - 2) * REAL(iexp - 1,dp) * rap(3) ) 00754 END IF 00755 END DO 00756 00757 END DO ! iz 00758 END DO ! iy 00759 END DO ! ix 00760 END DO ! iat 00761 END IF ! gth_gspace 00762 00763 ! ------------------------------------- 00764 ! Atomic grids part 00765 ! ------------------------------------- 00766 00767 IF (gapw .AND. paw_atom) THEN 00768 00769 bo_atom = get_limit(natom, para_env%num_pe, para_env%mepos) 00770 00771 DO iat = bo_atom(1),bo_atom(2) ! natom = # atoms for ikind 00772 00773 iatom = atom_list(iat) 00774 00775 nablavks_vec_rad_h => nablavks_atom_set(iatom)%nablavks_vec_rad_h 00776 nablavks_vec_rad_s => nablavks_atom_set(iatom)%nablavks_vec_rad_s 00777 00778 DO ir = 1,grid_atom%nr 00779 00780 IF (grid_atom%rad(ir) >= atomic_kind%hard_radius) CYCLE 00781 00782 exp_rap = EXP( - alpha * grid_atom%rad(ir)**2 ) 00783 00784 DO ia = 1,grid_atom%ng_sphere 00785 00786 DO idir = 1,3 00787 hard_value = 0.0_dp 00788 hard_value = charge * ( & 00789 - 2.0_dp * sqrt_alpha * EXP( - grid_atom%rad(ir)**2 * sqrt_alpha**2 ) & 00790 * grid_atom%rad(ir)*harmonics%a(idir,ia) & 00791 / ( rootpi * grid_atom%rad(ir)**2 ) & 00792 + erf( grid_atom%rad(ir) * sqrt_alpha ) & 00793 * grid_atom%rad(ir)*harmonics%a(idir,ia) & 00794 / grid_atom%rad(ir)**3 ) 00795 soft_value = hard_value 00796 DO iexp = 1,nexp_ppl 00797 hard_value = hard_value + ( & 00798 - 2.0_dp * alpha * grid_atom%rad(ir)*harmonics%a(idir,ia) & 00799 * exp_rap * cexp_ppl(iexp) * ( grid_atom%rad(ir)**2 )**(iexp - 1) ) 00800 IF (iexp > 1) THEN 00801 hard_value = hard_value + ( & 00802 2.0_dp * exp_rap * cexp_ppl(iexp) & 00803 * ( grid_atom%rad(ir)**2 )**(iexp - 2) * REAL(iexp - 1,dp) 00804 * grid_atom%rad(ir)*harmonics%a(idir,ia) ) 00805 END IF 00806 END DO 00807 nablavks_vec_rad_h(idir,1)%r_coef(ir,ia) = & 00808 nablavks_vec_rad_h(idir,1)%r_coef(ir,ia) + hard_value 00809 IF (make_soft) THEN 00810 nablavks_vec_rad_s(idir,1)%r_coef(ir,ia) = & 00811 nablavks_vec_rad_s(idir,1)%r_coef(ir,ia) + soft_value 00812 ELSE 00813 nablavks_vec_rad_s(idir,1)%r_coef(ir,ia) = & 00814 nablavks_vec_rad_s(idir,1)%r_coef(ir,ia) + hard_value 00815 END IF 00816 END DO 00817 00818 END DO ! ia 00819 END DO ! ir 00820 00821 DO ispin = 2,nspins 00822 DO idir = 1,3 00823 nablavks_vec_rad_h(idir,ispin)%r_coef(:,:) = nablavks_vec_rad_h(idir,1)%r_coef(:,:) 00824 nablavks_vec_rad_s(idir,ispin)%r_coef(:,:) = nablavks_vec_rad_s(idir,1)%r_coef(:,:) 00825 END DO 00826 END DO 00827 00828 END DO 00829 00830 END IF 00831 00832 ELSE IF (ASSOCIATED(all_potential)) THEN 00833 00834 CALL get_potential(potential=all_potential,& 00835 alpha_core_charge=alpha,& 00836 zeff=charge) 00837 00838 sqrt_alpha = SQRT(alpha) 00839 00840 ! ------------------------------------- 00841 ! Atomic grids part 00842 ! ------------------------------------- 00843 00844 bo_atom = get_limit(natom, para_env%num_pe, para_env%mepos) 00845 00846 DO iat = bo_atom(1),bo_atom(2) ! natom = # atoms for ikind 00847 00848 iatom = atom_list(iat) 00849 00850 nablavks_vec_rad_h => nablavks_atom_set(iatom)%nablavks_vec_rad_h 00851 00852 DO ir = 1,grid_atom%nr 00853 00854 IF (grid_atom%rad(ir) >= atomic_kind%hard_radius) CYCLE 00855 00856 DO ia = 1,grid_atom%ng_sphere 00857 00858 DO idir = 1,3 00859 hard_value = 0.0_dp 00860 hard_value = charge * ( & 00861 2.0_dp * sqrt_alpha * EXP( - grid_atom%rad(ir)**2 * sqrt_alpha**2 ) & 00862 * grid_atom%rad(ir)*harmonics%a(idir,ia) & 00863 / ( rootpi * grid_atom%rad(ir)**2 ) & 00864 + erfc( grid_atom%rad(ir) * sqrt_alpha ) & 00865 * grid_atom%rad(ir)*harmonics%a(idir,ia) & 00866 / grid_atom%rad(ir)**3 ) 00867 nablavks_vec_rad_h(idir,1)%r_coef(ir,ia) = & 00868 nablavks_vec_rad_h(idir,1)%r_coef(ir,ia) + hard_value 00869 END DO 00870 00871 END DO ! ia 00872 END DO ! ir 00873 00874 DO ispin = 2,nspins 00875 DO idir = 1,3 00876 nablavks_vec_rad_h(idir,ispin)%r_coef(:,:) = nablavks_vec_rad_h(idir,1)%r_coef(:,:) 00877 END DO 00878 END DO 00879 00880 END DO 00881 00882 ELSE 00883 CYCLE 00884 END IF 00885 00886 END DO 00887 00888 DO idir = 1,3 00889 CALL pw_copy(nablavks_set(idir,1)%rho%rho_r(1)%pw,nablavks_set(idir,2)%rho%rho_r(1)%pw,& 00890 error=error) 00891 END DO 00892 00893 ! ------------------------------------- 00894 ! Add V_xc potential 00895 ! ------------------------------------- 00896 00897 ALLOCATE(v_xc_gtemp,STAT=istat) 00898 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00899 "v_xc_gtemp",0) 00900 ALLOCATE(v_xc_rtemp,STAT=istat) 00901 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00902 "v_xc_rtemp",0) 00903 00904 CALL pw_pool_create_pw(auxbas_pw_pool,v_xc_gtemp%pw, & 00905 use_data=COMPLEXDATA1D,in_space=RECIPROCALSPACE,& 00906 error=error) 00907 CALL pw_pool_create_pw(auxbas_pw_pool,v_xc_rtemp%pw,& 00908 use_data=REALDATA3D,in_space=REALSPACE,& 00909 error=error) 00910 00911 xc_section => section_vals_get_subs_vals(input, "PROPERTIES%LINRES%EPR%PRINT%G_TENSOR%XC", error=error) 00912 00913 IF (gapw_xc) THEN 00914 CALL qs_vxc_create(qs_env=qs_env, rho_struct=qs_env%rho_xc, xc_section=xc_section,& 00915 vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE., & 00916 error=error) 00917 ELSE 00918 CALL qs_vxc_create(qs_env=qs_env, rho_struct=qs_env%rho, xc_section=xc_section,& 00919 vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE., & 00920 error=error) 00921 END IF 00922 00923 IF (ASSOCIATED(v_rspace_new)) THEN 00924 00925 DO ispin = 1,nspins 00926 00927 CALL pw_transfer(v_rspace_new(ispin)%pw, v_xc_gtemp%pw, error=error) 00928 CALL pw_derive(v_xc_gtemp%pw, (/1,0,0/) , error=error) 00929 CALL pw_transfer(v_xc_gtemp%pw, v_xc_rtemp%pw, error=error) 00930 CALL pw_axpy(v_xc_rtemp%pw, nablavks_set(1,ispin)%rho%rho_r(1)%pw, error=error) 00931 00932 CALL pw_transfer(v_rspace_new(ispin)%pw, v_xc_gtemp%pw, error=error) 00933 CALL pw_derive(v_xc_gtemp%pw, (/0,1,0/) , error=error) 00934 CALL pw_transfer(v_xc_gtemp%pw, v_xc_rtemp%pw, error=error) 00935 CALL pw_axpy(v_xc_rtemp%pw, nablavks_set(2,ispin)%rho%rho_r(1)%pw, error=error) 00936 00937 CALL pw_transfer(v_rspace_new(ispin)%pw, v_xc_gtemp%pw, error=error) 00938 CALL pw_derive(v_xc_gtemp%pw, (/0,0,1/) , error=error) 00939 CALL pw_transfer(v_xc_gtemp%pw, v_xc_rtemp%pw, error=error) 00940 CALL pw_axpy(v_xc_rtemp%pw, nablavks_set(3,ispin)%rho%rho_r(1)%pw, error=error) 00941 00942 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_rspace_new(ispin)%pw,& 00943 error=error) 00944 00945 END DO 00946 00947 DEALLOCATE(v_rspace_new,stat=istat) 00948 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00949 "v_rspace_new") 00950 END IF 00951 00952 IF (ASSOCIATED(v_tau_rspace)) THEN 00953 00954 DO ispin = 1,nspins 00955 00956 CALL pw_transfer(v_tau_rspace(ispin)%pw, v_xc_gtemp%pw, error=error) 00957 CALL pw_derive(v_xc_gtemp%pw, (/1,0,0/) , error=error) 00958 CALL pw_transfer(v_xc_gtemp%pw, v_xc_rtemp%pw, error=error) 00959 CALL pw_axpy(v_xc_rtemp%pw, nablavks_set(1,ispin)%rho%rho_r(1)%pw, error=error) 00960 00961 CALL pw_transfer(v_tau_rspace(ispin)%pw, v_xc_gtemp%pw, error=error) 00962 CALL pw_derive(v_xc_gtemp%pw, (/0,1,0/) , error=error) 00963 CALL pw_transfer(v_xc_gtemp%pw, v_xc_rtemp%pw, error=error) 00964 CALL pw_axpy(v_xc_rtemp%pw, nablavks_set(2,ispin)%rho%rho_r(1)%pw, error=error) 00965 00966 CALL pw_transfer(v_tau_rspace(ispin)%pw, v_xc_gtemp%pw, error=error) 00967 CALL pw_derive(v_xc_gtemp%pw, (/0,0,1/) , error=error) 00968 CALL pw_transfer(v_xc_gtemp%pw, v_xc_rtemp%pw, error=error) 00969 CALL pw_axpy(v_xc_rtemp%pw, nablavks_set(3,ispin)%rho%rho_r(1)%pw, error=error) 00970 00971 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_tau_rspace(ispin)%pw,& 00972 error=error) 00973 00974 END DO 00975 00976 DEALLOCATE(v_tau_rspace,stat=istat) 00977 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00978 "v_tau_rspace") 00979 END IF 00980 00981 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_xc_gtemp%pw,error=error) 00982 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_xc_rtemp%pw,error=error) 00983 00984 DEALLOCATE(v_xc_gtemp,STAT=istat) 00985 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00986 "v_xc_gtemp") 00987 DEALLOCATE(v_xc_rtemp,STAT=istat) 00988 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00989 "v_xc_rtemp") 00990 00991 IF (gapw .OR. gapw_xc) THEN 00992 CALL calculate_vxc_atom(qs_env=qs_env,energy_only=.FALSE.,& 00993 gradient_atom_set=nablavks_atom_set,& 00994 error=error) 00995 END IF 00996 00997 ! ------------------------------------- 00998 ! Write Nabla V_KS (local) to cubes 00999 ! ------------------------------------- 01000 01001 IF (BTEST(cp_print_key_should_output(logger%iter_info,lr_section,& 01002 "EPR%PRINT%NABLAVKS_CUBES",error=error),cp_p_file)) THEN 01003 ALLOCATE(wf_r,STAT=istat) 01004 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01005 "wf_r",0) 01006 CALL pw_pool_create_pw(auxbas_pw_pool,wf_r%pw,& 01007 use_data = REALDATA3D,& 01008 in_space = REALSPACE, error=error) 01009 DO idir = 1,3 01010 CALL pw_zero(wf_r%pw, error=error) 01011 CALL pw_copy(nablavks_set(idir,1)%rho%rho_r(1)%pw,wf_r%pw, error=error) ! RA 01012 filename="nablavks" 01013 WRITE(ext,'(a2,I1,a5)') "_d",idir,".cube" 01014 unit_nr=cp_print_key_unit_nr(logger,lr_section,"EPR%PRINT%NABLAVKS_CUBES",& 01015 extension=TRIM(ext),middle_name=TRIM(filename),& 01016 log_filename=.FALSE.,file_position="REWIND",error=error) 01017 CALL pw_to_cube(wf_r%pw,unit_nr,"NABLA V_KS ",& 01018 particles=particles,& 01019 stride=section_get_ivals(lr_section,& 01020 "EPR%PRINT%NABLAVKS_CUBES%STRIDE",error=error),& 01021 error=error) 01022 CALL cp_print_key_finished_output(unit_nr,logger,lr_section,& 01023 "EPR%PRINT%NABLAVKS_CUBES",error=error) 01024 END DO 01025 CALL pw_pool_give_back_pw(auxbas_pw_pool,wf_r%pw,error=error) 01026 DEALLOCATE(wf_r,STAT=istat) 01027 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01028 "wf_r") 01029 END IF 01030 01031 CALL cp_print_key_finished_output(output_unit,logger,lr_section,& 01032 "PRINT%PROGRAM_RUN_INFO",error=error) 01033 01034 END SUBROUTINE epr_nablavks 01035 01036 END MODULE qs_linres_epr_nablavks 01037
1.7.3