|
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 ! ***************************************************************************** 00028 MODULE qs_integrate_potential 00029 USE atomic_kind_types, ONLY: atomic_kind_type,& 00030 get_atomic_kind,& 00031 get_atomic_kind_set 00032 USE atprop_types, ONLY: atprop_array_init 00033 USE basis_set_types, ONLY: get_gto_basis_set,& 00034 gto_basis_set_type 00035 USE cell_types, ONLY: cell_type,& 00036 pbc 00037 USE cp_array_r_utils, ONLY: cp_2d_r_p_type 00038 USE cp_control_types, ONLY: dft_control_type 00039 USE cp_dbcsr_interface, ONLY: & 00040 cp_dbcsr_col_block_sizes, cp_dbcsr_copy, cp_dbcsr_create, & 00041 cp_dbcsr_distribution, cp_dbcsr_finalize, cp_dbcsr_get_block_p, & 00042 cp_dbcsr_get_data_size, cp_dbcsr_get_matrix_type, & 00043 cp_dbcsr_get_num_blocks, cp_dbcsr_init, cp_dbcsr_row_block_sizes, & 00044 cp_dbcsr_work_create 00045 USE cp_dbcsr_operations, ONLY: cp_dbcsr_add_block_node,& 00046 cp_dbcsr_deallocate_matrix 00047 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type,& 00048 cp_dbcsr_type 00049 USE cp_para_types, ONLY: cp_para_env_type 00050 USE cube_utils, ONLY: cube_info_type 00051 USE dbcsr_dist_operations 00052 USE dbcsr_methods, ONLY: dbcsr_distribution_has_threads 00053 USE dbcsr_types, ONLY: dbcsr_distribution_obj 00054 USE external_potential_types, ONLY: get_potential,& 00055 gth_potential_type 00056 USE gaussian_gridlevels, ONLY: gridlevel_info_type 00057 USE input_constants, ONLY: pw_interp,& 00058 spline3_pbc_interp,& 00059 use_aux_fit_basis_set,& 00060 use_orb_basis_set 00061 USE input_section_types, ONLY: section_vals_get_subs_vals,& 00062 section_vals_type,& 00063 section_vals_val_get 00064 USE kinds, ONLY: dp,& 00065 int_8 00066 USE mathconstants, ONLY: dfac,& 00067 pi 00068 USE memory_utilities, ONLY: reallocate 00069 USE orbital_pointers, ONLY: coset,& 00070 nco,& 00071 ncoset,& 00072 nso,& 00073 nsoset 00074 USE orbital_transformation_matrices, ONLY: orbtramat 00075 USE particle_types, ONLY: particle_type 00076 USE pw_env_types, ONLY: pw_env_get,& 00077 pw_env_type 00078 USE pw_methods, ONLY: pw_copy,& 00079 pw_transfer,& 00080 pw_zero 00081 USE pw_pool_types, ONLY: pw_pool_p_type,& 00082 pw_pools_create_pws,& 00083 pw_pools_give_back_pws 00084 USE pw_spline_utils, ONLY: pw_restrict_s3 00085 USE pw_types, ONLY: COMPLEXDATA1D,& 00086 REALDATA3D,& 00087 REALSPACE,& 00088 RECIPROCALSPACE,& 00089 pw_p_type 00090 USE qs_environment_types, ONLY: get_qs_env,& 00091 qs_environment_type 00092 USE qs_force_types, ONLY: qs_force_type 00093 USE qs_integrate_potential_low, ONLY: integrate_pgf_product_rspace 00094 USE realspace_grid_types, ONLY: pw2rs,& 00095 realspace_grid_desc_p_type,& 00096 realspace_grid_desc_type,& 00097 realspace_grid_p_type,& 00098 realspace_grid_type,& 00099 rs_grid_release,& 00100 rs_grid_retain,& 00101 rs_pw_transfer 00102 USE scptb_types, ONLY: get_scptb_parameter,& 00103 scp_vector_type,& 00104 scptb_parameter_type 00105 USE task_list_methods, ONLY: int2pair,& 00106 rs_distribute_matrix 00107 USE task_list_types, ONLY: task_list_type 00108 USE termination, ONLY: stop_program 00109 USE timings, ONLY: timeset,& 00110 timestop 00111 USE virial_types, ONLY: virial_type 00112 #include "cp_common_uses.h" 00113 00114 IMPLICIT NONE 00115 00116 PRIVATE 00117 00118 INTEGER :: debug_count=0 00119 00120 LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.FALSE. 00121 00122 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_integrate_potential' 00123 00124 ! *** Public subroutines *** 00125 00126 PUBLIC :: integrate_v_rspace,& 00127 integrate_v_core_rspace,& 00128 integrate_ppl_rspace,& 00129 integrate_scp_rspace,& 00130 integrate_pgf_product_rspace,& 00131 potential_pw2rs, integrate_rho_nlcc 00132 00133 CONTAINS 00134 00135 ! ***************************************************************************** 00138 SUBROUTINE integrate_scp_rspace(scp_pot,qs_env,scpv,calculate_forces,error) 00139 00140 TYPE(pw_p_type), INTENT(INOUT) :: scp_pot 00141 TYPE(qs_environment_type), POINTER :: qs_env 00142 TYPE(scp_vector_type), POINTER :: scpv 00143 LOGICAL, INTENT(IN) :: calculate_forces 00144 TYPE(cp_error_type), INTENT(inout) :: error 00145 00146 CHARACTER(len=*), PARAMETER :: routineN = 'integrate_scp_rspace', 00147 routineP = moduleN//':'//routineN 00148 00149 INTEGER :: atom_a, handle, i, iatom, 00150 ierr, ii, ikind, j, jj, l, 00151 lmaxscp, natom_of_kind, ni, 00152 nj, npme, stat 00153 INTEGER, DIMENSION(:), POINTER :: atom_list, cores 00154 LOGICAL :: defined, failure, use_virial 00155 REAL(KIND=dp) :: alpha, dvol, eps_rho_rspace, 00156 norm, pp 00157 REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra 00158 REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b 00159 REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab 00160 TYPE(atomic_kind_type), DIMENSION(:), 00161 POINTER :: atomic_kind_set 00162 TYPE(atomic_kind_type), POINTER :: atomic_kind 00163 TYPE(cell_type), POINTER :: cell 00164 TYPE(cp_para_env_type), POINTER :: para_env 00165 TYPE(dft_control_type), POINTER :: dft_control 00166 TYPE(particle_type), DIMENSION(:), 00167 POINTER :: particle_set 00168 TYPE(pw_env_type), POINTER :: pw_env 00169 TYPE(qs_force_type), DIMENSION(:), 00170 POINTER :: force 00171 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 00172 TYPE(realspace_grid_type), POINTER :: rs_v 00173 TYPE(scptb_parameter_type), POINTER :: scptb_kind 00174 TYPE(virial_type), POINTER :: virial 00175 00176 CALL timeset(routineN,handle) 00177 00178 failure=.FALSE. 00179 NULLIFY(pw_env,auxbas_rs_desc,cores) 00180 00181 CALL get_qs_env(qs_env=qs_env,pw_env=pw_env,error=error) 00182 CALL pw_env_get(pw_env=pw_env,auxbas_rs_desc=auxbas_rs_desc, & 00183 auxbas_rs_grid=rs_v,error=error) 00184 CALL rs_grid_retain(rs_v,error=error) 00185 00186 CALL rs_pw_transfer(rs_v,scp_pot%pw,pw2rs,error=error) 00187 00188 CALL get_qs_env(qs_env=qs_env,& 00189 atomic_kind_set=atomic_kind_set,& 00190 cell=cell,& 00191 dft_control=dft_control,& 00192 particle_set=particle_set,& 00193 para_env=para_env,pw_env=pw_env,& 00194 force=force,virial=virial,error=error) 00195 00196 use_virial = virial%pv_availability.AND.(.NOT.virial%pv_numer).AND.calculate_forces 00197 00198 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 00199 dvol = scp_pot%pw%pw_grid%dvol 00200 00201 DO ikind=1,SIZE(atomic_kind_set) 00202 00203 atomic_kind => atomic_kind_set(ikind) 00204 00205 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00206 natom=natom_of_kind,atom_list=atom_list,& 00207 scptb_parameter=scptb_kind) 00208 CALL get_scptb_parameter(scptb_kind,defined=defined,lmaxscp=lmaxscp,ag=alpha) 00209 IF (.NOT.defined) CYCLE 00210 00211 ni = ncoset(lmaxscp) 00212 ALLOCATE(hab(ni,1),pab(ni,1),STAT=stat) 00213 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00214 pab = 0._dp 00215 hab = 0._dp 00216 00217 ALLOCATE(cores(natom_of_kind),STAT=ierr) 00218 CPPrecondition(ierr==0,cp_failure_level,routineP,error,failure) 00219 cores = 0 00220 npme = 0 00221 00222 DO iatom = 1, natom_of_kind 00223 atom_a = atom_list(iatom) 00224 ra(:) = pbc(particle_set(atom_a)%r,cell) 00225 IF(rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN 00226 ! replicated realspace grid, split the atoms up between procs 00227 IF (MODULO(iatom,rs_v%desc%group_size) == rs_v % desc % my_pos ) THEN 00228 npme = npme + 1 00229 cores (npme) = iatom 00230 ENDIF 00231 ELSE 00232 npme = npme + 1 00233 cores (npme) = iatom 00234 ENDIF 00235 END DO 00236 00237 DO nj=1,npme 00238 00239 iatom = cores(nj) 00240 atom_a = atom_list(iatom) 00241 ra(:) = pbc(particle_set(atom_a)%r,cell) 00242 hab(:,1) = 0.0_dp 00243 IF (calculate_forces) THEN 00244 force_a(:) = 0.0_dp 00245 force_b(:) = 0.0_dp 00246 IF (use_virial) THEN 00247 my_virial_a = 0.0_dp 00248 my_virial_b = 0.0_dp 00249 END IF 00250 END IF 00251 00252 CALL integrate_pgf_product_rspace(lmaxscp,alpha,0,& 00253 0,0.0_dp,0,ra,(/0.0_dp,0.0_dp,0.0_dp/),0.0_dp,& 00254 rs_v,cell,pw_env%cube_info(1),hab,pab=pab,o1=0,o2=0,& 00255 eps_gvg_rspace=eps_rho_rspace,& 00256 calculate_forces=calculate_forces,force_a=force_a,& 00257 force_b=force_b,use_virial=use_virial,my_virial_a=my_virial_a,& 00258 my_virial_b=my_virial_b,use_subpatch=.TRUE.,subpatch_pattern=0_int_8,error=error) 00259 00260 DO l=0,lmaxscp 00261 pp = (2._dp*l+3._dp)/2._dp 00262 norm = 2._dp**(l+2)/SQRT(pi)/dfac(2*l+1) * dvol 00263 norm = SQRT(0.25_dp*dfac(2*l+1)/pi) * norm*alpha**pp 00264 DO jj=1,nco(l) 00265 j = ncoset(l-1) + jj 00266 DO ii=1,nso(l) 00267 i = nsoset(l-1) + ii 00268 scpv%vector(ikind)%vmat(i,iatom) = scpv%vector(ikind)%vmat(i,iatom) + & 00269 hab(j,1)*orbtramat(l)%c2s(ii,jj)*norm 00270 END DO 00271 END DO 00272 END DO 00273 00274 IF (calculate_forces) THEN 00275 force(ikind)%gth_ppl(:,iatom) = force(ikind)%gth_ppl(:,iatom) + force_a(:)*dvol 00276 IF (use_virial) THEN 00277 virial%pv_virial = virial%pv_virial + my_virial_a*dvol 00278 END IF 00279 END IF 00280 END DO 00281 00282 DEALLOCATE(cores,hab,pab,STAT=stat) 00283 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00284 00285 END DO 00286 00287 CALL rs_grid_release(rs_v, error=error) 00288 00289 CALL timestop(handle) 00290 00291 END SUBROUTINE integrate_scp_rspace 00292 ! ***************************************************************************** 00295 SUBROUTINE integrate_ppl_rspace(rho_rspace,qs_env,error) 00296 TYPE(pw_p_type), INTENT(INOUT) :: rho_rspace 00297 TYPE(qs_environment_type), POINTER :: qs_env 00298 TYPE(cp_error_type), INTENT(inout) :: error 00299 00300 CHARACTER(len=*), PARAMETER :: routineN = 'integrate_ppl_rspace', 00301 routineP = moduleN//':'//routineN 00302 00303 INTEGER :: atom_a, handle, iatom, ikind, 00304 j, lppl, n, natom_of_kind, 00305 ni, npme, stat 00306 INTEGER, DIMENSION(:), POINTER :: atom_list, cores 00307 LOGICAL :: failure, use_virial 00308 REAL(KIND=dp) :: alpha, eps_rho_rspace 00309 REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra 00310 REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b 00311 REAL(KIND=dp), DIMENSION(:), POINTER :: cexp_ppl 00312 REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab 00313 TYPE(atomic_kind_type), DIMENSION(:), 00314 POINTER :: atomic_kind_set 00315 TYPE(atomic_kind_type), POINTER :: atomic_kind 00316 TYPE(cell_type), POINTER :: cell 00317 TYPE(cp_para_env_type), POINTER :: para_env 00318 TYPE(dft_control_type), POINTER :: dft_control 00319 TYPE(gth_potential_type), POINTER :: gth_potential 00320 TYPE(particle_type), DIMENSION(:), 00321 POINTER :: particle_set 00322 TYPE(pw_env_type), POINTER :: pw_env 00323 TYPE(qs_force_type), DIMENSION(:), 00324 POINTER :: force 00325 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 00326 TYPE(realspace_grid_type), POINTER :: rs_v 00327 TYPE(virial_type), POINTER :: virial 00328 00329 CALL timeset(routineN,handle) 00330 00331 failure=.FALSE. 00332 NULLIFY(pw_env,auxbas_rs_desc,cores) 00333 00334 CALL get_qs_env(qs_env=qs_env,pw_env=pw_env,error=error) 00335 CALL pw_env_get(pw_env=pw_env,auxbas_rs_desc=auxbas_rs_desc, & 00336 auxbas_rs_grid=rs_v,error=error) 00337 CALL rs_grid_retain(rs_v,error=error) 00338 00339 CALL rs_pw_transfer(rs_v,rho_rspace%pw,pw2rs,error=error) 00340 00341 CALL get_qs_env(qs_env=qs_env,& 00342 atomic_kind_set=atomic_kind_set,& 00343 cell=cell,& 00344 dft_control=dft_control,& 00345 particle_set=particle_set,& 00346 para_env=para_env,pw_env=pw_env,& 00347 force=force,virial=virial,error=error) 00348 00349 use_virial = virial%pv_availability.AND.(.NOT.virial%pv_numer) 00350 00351 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 00352 00353 DO ikind=1,SIZE(atomic_kind_set) 00354 00355 atomic_kind => atomic_kind_set(ikind) 00356 00357 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00358 natom=natom_of_kind,& 00359 atom_list=atom_list,& 00360 gth_potential=gth_potential) 00361 00362 IF (.NOT.ASSOCIATED(gth_potential)) CYCLE 00363 CALL get_potential(potential=gth_potential,alpha_ppl=alpha,nexp_ppl=lppl,cexp_ppl=cexp_ppl) 00364 00365 IF ( lppl <= 0 ) CYCLE 00366 00367 ni = ncoset(2*lppl-2) 00368 ALLOCATE(hab(ni,1),pab(ni,1),STAT=stat) 00369 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00370 pab = 0._dp 00371 00372 CALL reallocate ( cores, 1, natom_of_kind ) 00373 npme = 0 00374 cores = 0 00375 00376 ! prepare core function 00377 DO j=1,lppl 00378 SELECT CASE (j) 00379 CASE (1) 00380 pab(1,1) = cexp_ppl(1) 00381 CASE (2) 00382 n = coset(2,0,0) 00383 pab(n,1) = cexp_ppl(2) 00384 n = coset(0,2,0) 00385 pab(n,1) = cexp_ppl(2) 00386 n = coset(0,0,2) 00387 pab(n,1) = cexp_ppl(2) 00388 CASE (3) 00389 n = coset(4,0,0) 00390 pab(n,1) = cexp_ppl(3) 00391 n = coset(0,4,0) 00392 pab(n,1) = cexp_ppl(3) 00393 n = coset(0,0,4) 00394 pab(n,1) = cexp_ppl(3) 00395 n = coset(2,2,0) 00396 pab(n,1) = 2._dp*cexp_ppl(3) 00397 n = coset(2,0,2) 00398 pab(n,1) = 2._dp*cexp_ppl(3) 00399 n = coset(0,2,2) 00400 pab(n,1) = 2._dp*cexp_ppl(3) 00401 CASE (4) 00402 n = coset(6,0,0) 00403 pab(n,1) = cexp_ppl(4) 00404 n = coset(0,6,0) 00405 pab(n,1) = cexp_ppl(4) 00406 n = coset(0,0,6) 00407 pab(n,1) = cexp_ppl(4) 00408 n = coset(4,2,0) 00409 pab(n,1) = 3._dp*cexp_ppl(4) 00410 n = coset(4,0,2) 00411 pab(n,1) = 3._dp*cexp_ppl(4) 00412 n = coset(2,4,0) 00413 pab(n,1) = 3._dp*cexp_ppl(4) 00414 n = coset(2,0,4) 00415 pab(n,1) = 3._dp*cexp_ppl(4) 00416 n = coset(0,4,2) 00417 pab(n,1) = 3._dp*cexp_ppl(4) 00418 n = coset(0,2,4) 00419 pab(n,1) = 3._dp*cexp_ppl(4) 00420 n = coset(2,2,2) 00421 pab(n,1) = 6._dp*cexp_ppl(4) 00422 CASE DEFAULT 00423 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00424 END SELECT 00425 END DO 00426 00427 DO iatom = 1, natom_of_kind 00428 atom_a = atom_list(iatom) 00429 ra(:) = pbc(particle_set(atom_a)%r,cell) 00430 IF(rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN 00431 ! replicated realspace grid, split the atoms up between procs 00432 IF (MODULO(iatom,rs_v%desc%group_size) == rs_v % desc % my_pos ) THEN 00433 npme = npme + 1 00434 cores (npme) = iatom 00435 ENDIF 00436 ELSE 00437 npme = npme + 1 00438 cores (npme) = iatom 00439 ENDIF 00440 END DO 00441 00442 DO j=1,npme 00443 00444 iatom = cores(j) 00445 atom_a = atom_list(iatom) 00446 ra(:) = pbc(particle_set(atom_a)%r,cell) 00447 hab(:,1) = 0.0_dp 00448 force_a(:) = 0.0_dp 00449 force_b(:) = 0.0_dp 00450 IF (use_virial) THEN 00451 my_virial_a = 0.0_dp 00452 my_virial_b = 0.0_dp 00453 END IF 00454 ni = 2*lppl-2 00455 00456 CALL integrate_pgf_product_rspace(ni,alpha,0,& 00457 0,0.0_dp,0,ra,(/0.0_dp,0.0_dp,0.0_dp/),0.0_dp,& 00458 rs_v,cell,pw_env%cube_info(1),hab,pab=pab,o1=0,o2=0,& 00459 eps_gvg_rspace=eps_rho_rspace,& 00460 calculate_forces=.TRUE.,force_a=force_a,& 00461 force_b=force_b,use_virial=use_virial,my_virial_a=my_virial_a,& 00462 my_virial_b=my_virial_b,use_subpatch=.TRUE.,subpatch_pattern=0_int_8,error=error) 00463 00464 force(ikind)%gth_ppl(:,iatom) =& 00465 force(ikind)%gth_ppl(:,iatom) + force_a(:)*rho_rspace%pw%pw_grid%dvol 00466 00467 IF (use_virial) THEN 00468 virial%pv_virial = virial%pv_virial + my_virial_a*rho_rspace%pw%pw_grid%dvol 00469 CALL cp_unimplemented_error(fromWhere=routineP, & 00470 message="Virial not debuged for CORE_PPL", & 00471 error=error, error_level=cp_failure_level) 00472 END IF 00473 END DO 00474 00475 DEALLOCATE(hab,pab,STAT=stat) 00476 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00477 00478 END DO 00479 00480 CALL rs_grid_release(rs_v, error=error) 00481 00482 DEALLOCATE(cores,STAT=stat) 00483 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00484 00485 CALL timestop(handle) 00486 00487 END SUBROUTINE integrate_ppl_rspace 00488 00489 ! ***************************************************************************** 00492 SUBROUTINE integrate_rho_nlcc(rho_rspace,qs_env,error) 00493 TYPE(pw_p_type), INTENT(INOUT) :: rho_rspace 00494 TYPE(qs_environment_type), POINTER :: qs_env 00495 TYPE(cp_error_type), INTENT(inout) :: error 00496 00497 CHARACTER(len=*), PARAMETER :: routineN = 'integrate_rho_nlcc', 00498 routineP = moduleN//':'//routineN 00499 00500 INTEGER :: atom_a, handle, iatom, 00501 iexp_nlcc, ikind, ithread, j, 00502 n, natom, nc, nexp_nlcc, ni, 00503 npme, nthread, stat 00504 INTEGER, DIMENSION(:), POINTER :: atom_list, cores, nct_nlcc 00505 LOGICAL :: failure, nlcc, use_virial 00506 REAL(KIND=dp) :: alpha, eps_rho_rspace 00507 REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra 00508 REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b 00509 REAL(KIND=dp), DIMENSION(:), POINTER :: alpha_nlcc 00510 REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_nlcc, hab, pab 00511 TYPE(atomic_kind_type), DIMENSION(:), 00512 POINTER :: atomic_kind_set 00513 TYPE(atomic_kind_type), POINTER :: atomic_kind 00514 TYPE(cell_type), POINTER :: cell 00515 TYPE(cp_para_env_type), POINTER :: para_env 00516 TYPE(dft_control_type), POINTER :: dft_control 00517 TYPE(gth_potential_type), POINTER :: gth_potential 00518 TYPE(particle_type), DIMENSION(:), 00519 POINTER :: particle_set 00520 TYPE(pw_env_type), POINTER :: pw_env 00521 TYPE(qs_force_type), DIMENSION(:), 00522 POINTER :: force 00523 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 00524 TYPE(realspace_grid_type), POINTER :: rs_v 00525 TYPE(virial_type), POINTER :: virial 00526 00527 CALL timeset(routineN,handle) 00528 00529 failure=.FALSE. 00530 NULLIFY(pw_env,auxbas_rs_desc,cores) 00531 00532 CALL get_qs_env(qs_env=qs_env,pw_env=pw_env,error=error) 00533 CALL pw_env_get(pw_env=pw_env,auxbas_rs_desc=auxbas_rs_desc, & 00534 auxbas_rs_grid=rs_v,error=error) 00535 CALL rs_grid_retain(rs_v,error=error) 00536 00537 CALL rs_pw_transfer(rs_v,rho_rspace%pw,pw2rs,error=error) 00538 00539 CALL get_qs_env(qs_env=qs_env,& 00540 atomic_kind_set=atomic_kind_set,& 00541 cell=cell,& 00542 dft_control=dft_control,& 00543 particle_set=particle_set,& 00544 para_env=para_env,pw_env=pw_env,& 00545 force=force,virial=virial,error=error) 00546 00547 use_virial = virial%pv_availability.AND.(.NOT.virial%pv_numer) 00548 00549 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 00550 00551 DO ikind=1,SIZE(atomic_kind_set) 00552 00553 atomic_kind => atomic_kind_set(ikind) 00554 00555 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00556 natom=natom,& 00557 atom_list=atom_list,& 00558 gth_potential=gth_potential) 00559 00560 IF (.NOT.ASSOCIATED(gth_potential)) CYCLE 00561 CALL get_potential(potential=gth_potential,nlcc_present=nlcc,nexp_nlcc=nexp_nlcc,& 00562 alpha_nlcc=alpha_nlcc,nct_nlcc=nct_nlcc,cval_nlcc=cval_nlcc) 00563 00564 IF ( .NOT. nlcc ) CYCLE 00565 00566 DO iexp_nlcc=1,nexp_nlcc 00567 00568 alpha=alpha_nlcc(iexp_nlcc) 00569 nc=nct_nlcc(iexp_nlcc) 00570 00571 ni = ncoset(2*nc-2) 00572 00573 nthread = 1 00574 ithread=0 00575 00576 ALLOCATE(hab(ni,1),pab(ni,1),STAT=stat) 00577 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00578 pab = 0._dp 00579 00580 CALL reallocate ( cores, 1, natom ) 00581 npme = 0 00582 cores = 0 00583 00584 ! prepare core function 00585 DO j=1,nc 00586 SELECT CASE (j) 00587 CASE (1) 00588 pab(1,1) = cval_nlcc(1,iexp_nlcc) 00589 CASE (2) 00590 n = coset(2,0,0) 00591 pab(n,1) = cval_nlcc(2,iexp_nlcc)/alpha**2 00592 n = coset(0,2,0) 00593 pab(n,1) = cval_nlcc(2,iexp_nlcc)/alpha**2 00594 n = coset(0,0,2) 00595 pab(n,1) = cval_nlcc(2,iexp_nlcc)/alpha**2 00596 CASE (3) 00597 n = coset(4,0,0) 00598 pab(n,1) = cval_nlcc(3,iexp_nlcc)/alpha**4 00599 n = coset(0,4,0) 00600 pab(n,1) = cval_nlcc(3,iexp_nlcc)/alpha**4 00601 n = coset(0,0,4) 00602 pab(n,1) = cval_nlcc(3,iexp_nlcc)/alpha**4 00603 n = coset(2,2,0) 00604 pab(n,1) = 2._dp*cval_nlcc(3,iexp_nlcc)/alpha**4 00605 n = coset(2,0,2) 00606 pab(n,1) = 2._dp*cval_nlcc(3,iexp_nlcc)/alpha**4 00607 n = coset(0,2,2) 00608 pab(n,1) = 2._dp*cval_nlcc(3,iexp_nlcc)/alpha**4 00609 CASE (4) 00610 n = coset(6,0,0) 00611 pab(n,1) = cval_nlcc(4,iexp_nlcc)/alpha**6 00612 n = coset(0,6,0) 00613 pab(n,1) = cval_nlcc(4,iexp_nlcc)/alpha**6 00614 n = coset(0,0,6) 00615 pab(n,1) = cval_nlcc(4,iexp_nlcc)/alpha**6 00616 n = coset(4,2,0) 00617 pab(n,1) = 3._dp*cval_nlcc(4,iexp_nlcc)/alpha**6 00618 n = coset(4,0,2) 00619 pab(n,1) = 3._dp*cval_nlcc(4,iexp_nlcc)/alpha**6 00620 n = coset(2,4,0) 00621 pab(n,1) = 3._dp*cval_nlcc(4,iexp_nlcc)/alpha**6 00622 n = coset(2,0,4) 00623 pab(n,1) = 3._dp*cval_nlcc(4,iexp_nlcc)/alpha**6 00624 n = coset(0,4,2) 00625 pab(n,1) = 3._dp*cval_nlcc(4,iexp_nlcc)/alpha**6 00626 n = coset(0,2,4) 00627 pab(n,1) = 3._dp*cval_nlcc(4,iexp_nlcc)/alpha**6 00628 n = coset(2,2,2) 00629 pab(n,1) = 6._dp*cval_nlcc(4,iexp_nlcc)/alpha**6 00630 CASE DEFAULT 00631 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00632 END SELECT 00633 END DO 00634 IF(dft_control%nspins==2)pab=pab*0.5_dp 00635 00636 DO iatom = 1, natom 00637 atom_a = atom_list(iatom) 00638 ra(:) = pbc(particle_set(atom_a)%r,cell) 00639 IF(rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN 00640 ! replicated realspace grid, split the atoms up between procs 00641 IF (MODULO(iatom,rs_v%desc%group_size) == rs_v % desc % my_pos ) THEN 00642 npme = npme + 1 00643 cores (npme) = iatom 00644 ENDIF 00645 ELSE 00646 npme = npme + 1 00647 cores (npme) = iatom 00648 ENDIF 00649 END DO 00650 00651 DO j=1,npme 00652 00653 iatom = cores(j) 00654 atom_a = atom_list(iatom) 00655 ra(:) = pbc(particle_set(atom_a)%r,cell) 00656 hab(:,1) = 0.0_dp 00657 force_a(:) = 0.0_dp 00658 force_b(:) = 0.0_dp 00659 IF (use_virial) THEN 00660 my_virial_a = 0.0_dp 00661 my_virial_b = 0.0_dp 00662 END IF 00663 ni = 2*nc-2 00664 00665 CALL integrate_pgf_product_rspace(ni,1/(2*alpha**2),0,& 00666 0,0.0_dp,0,ra,(/0.0_dp,0.0_dp,0.0_dp/),0.0_dp,& 00667 rs_v,cell,pw_env%cube_info(1),hab,pab=pab,o1=0,o2=0,& 00668 eps_gvg_rspace=eps_rho_rspace,& 00669 calculate_forces=.TRUE.,force_a=force_a,& 00670 force_b=force_b,use_virial=use_virial,my_virial_a=my_virial_a,& 00671 my_virial_b=my_virial_b,use_subpatch=.TRUE.,subpatch_pattern=0_int_8,error=error) 00672 00673 force(ikind)%gth_nlcc(:,iatom) =& 00674 force(ikind)%gth_nlcc(:,iatom) + force_a(:)*rho_rspace%pw%pw_grid%dvol 00675 00676 IF (use_virial) THEN 00677 virial%pv_virial = virial%pv_virial + my_virial_a*rho_rspace%pw%pw_grid%dvol 00678 END IF 00679 END DO 00680 00681 DEALLOCATE(hab,pab,STAT=stat) 00682 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00683 00684 END DO 00685 00686 END DO 00687 00688 CALL rs_grid_release(rs_v, error=error) 00689 00690 DEALLOCATE(cores,STAT=stat) 00691 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00692 00693 CALL timestop(handle) 00694 00695 END SUBROUTINE integrate_rho_nlcc 00696 00697 ! ***************************************************************************** 00701 SUBROUTINE integrate_v_core_rspace(v_rspace,qs_env,error) 00702 TYPE(pw_p_type), INTENT(INOUT) :: v_rspace 00703 TYPE(qs_environment_type), POINTER :: qs_env 00704 TYPE(cp_error_type), INTENT(inout) :: error 00705 00706 CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_core_rspace', 00707 routineP = moduleN//':'//routineN 00708 00709 INTEGER :: atom_a, handle, iatom, ikind, 00710 j, natom, natom_of_kind, 00711 npme, stat 00712 INTEGER, DIMENSION(:), POINTER :: atom_list, cores 00713 LOGICAL :: atenergy, failure, paw_atom, 00714 skip_fcore, use_virial 00715 REAL(KIND=dp) :: alpha_core_charge, 00716 ccore_charge, eps_rho_rspace 00717 REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra 00718 REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b 00719 REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab 00720 TYPE(atomic_kind_type), DIMENSION(:), 00721 POINTER :: atomic_kind_set 00722 TYPE(atomic_kind_type), POINTER :: atomic_kind 00723 TYPE(cell_type), POINTER :: cell 00724 TYPE(cp_para_env_type), POINTER :: para_env 00725 TYPE(dft_control_type), POINTER :: dft_control 00726 TYPE(particle_type), DIMENSION(:), 00727 POINTER :: particle_set 00728 TYPE(pw_env_type), POINTER :: pw_env 00729 TYPE(qs_force_type), DIMENSION(:), 00730 POINTER :: force 00731 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 00732 TYPE(realspace_grid_type), POINTER :: rs_v 00733 TYPE(virial_type), POINTER :: virial 00734 00735 CALL timeset(routineN,handle) 00736 00737 !If gapw, check for gpw kinds 00738 skip_fcore = .FALSE. 00739 IF(qs_env%dft_control%qs_control%gapw) THEN 00740 IF(.NOT. qs_env%dft_control%qs_control%gapw_control%nopaw_as_gpw) skip_fcore = .TRUE. 00741 END IF 00742 ! atomic energy contributions 00743 atenergy = .FALSE. 00744 IF(ASSOCIATED(qs_env%atprop)) THEN 00745 atenergy=qs_env%atprop%energy 00746 END IF 00747 00748 failure = .FALSE. 00749 IF(.NOT. skip_fcore) THEN 00750 NULLIFY(pw_env,auxbas_rs_desc,hab,pab,cores) 00751 ALLOCATE (hab(1,1),STAT=stat) 00752 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00753 ALLOCATE (pab(1,1),STAT=stat) 00754 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00755 00756 CALL get_qs_env(qs_env=qs_env,pw_env=pw_env,error=error) 00757 CALL pw_env_get(pw_env=pw_env,auxbas_rs_desc=auxbas_rs_desc, & 00758 auxbas_rs_grid=rs_v,error=error) 00759 CALL rs_grid_retain(rs_v,error=error) 00760 00761 CALL rs_pw_transfer(rs_v,v_rspace%pw,pw2rs,error=error) 00762 00763 CALL get_qs_env(qs_env=qs_env,& 00764 atomic_kind_set=atomic_kind_set,& 00765 cell=cell,& 00766 dft_control=dft_control,& 00767 particle_set=particle_set,& 00768 para_env=para_env,pw_env=pw_env,& 00769 force=force,virial=virial,error=error) 00770 00771 IF(atenergy) THEN 00772 natom = SIZE(particle_set) 00773 CALL atprop_array_init(qs_env%atprop%ateb,natom,error) 00774 END IF 00775 00776 use_virial = virial%pv_availability.AND.(.NOT.virial%pv_numer) 00777 00778 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 00779 00780 DO ikind=1,SIZE(atomic_kind_set) 00781 00782 atomic_kind => atomic_kind_set(ikind) 00783 00784 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00785 natom=natom_of_kind,& 00786 paw_atom=paw_atom, & 00787 atom_list=atom_list,& 00788 alpha_core_charge=alpha_core_charge,& 00789 ccore_charge=ccore_charge) 00790 00791 IF(paw_atom) THEN 00792 force(ikind)%rho_core(:,:) = 0.0_dp 00793 CYCLE 00794 END IF 00795 pab(1,1) = -ccore_charge 00796 00797 IF (alpha_core_charge == 0.0_dp .OR. pab(1,1)== 0.0_dp) CYCLE 00798 00799 CALL reallocate ( cores, 1, natom_of_kind ) 00800 npme = 0 00801 cores = 0 00802 00803 DO iatom = 1, natom_of_kind 00804 atom_a = atom_list(iatom) 00805 ra(:) = pbc(particle_set(atom_a)%r,cell) 00806 IF(rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN 00807 ! replicated realspace grid, split the atoms up between procs 00808 IF (MODULO(iatom,rs_v%desc%group_size) == rs_v % desc % my_pos ) THEN 00809 npme = npme + 1 00810 cores (npme) = iatom 00811 ENDIF 00812 ELSE 00813 npme = npme + 1 00814 cores (npme) = iatom 00815 ENDIF 00816 END DO 00817 00818 DO j=1,npme 00819 00820 iatom = cores(j) 00821 atom_a = atom_list(iatom) 00822 ra(:) = pbc(particle_set(atom_a)%r,cell) 00823 hab(1,1) = 0.0_dp 00824 force_a(:) = 0.0_dp 00825 force_b(:) = 0.0_dp 00826 IF (use_virial) THEN 00827 my_virial_a = 0.0_dp 00828 my_virial_b = 0.0_dp 00829 END IF 00830 00831 CALL integrate_pgf_product_rspace(0,alpha_core_charge,0,& 00832 0,0.0_dp,0,ra,(/0.0_dp,0.0_dp,0.0_dp/),0.0_dp,& 00833 rs_v,cell,pw_env%cube_info(1),hab,pab=pab,o1=0,o2=0,& 00834 eps_gvg_rspace=eps_rho_rspace,& 00835 calculate_forces=.TRUE.,force_a=force_a,& 00836 force_b=force_b,use_virial=use_virial,my_virial_a=my_virial_a,& 00837 my_virial_b=my_virial_b,use_subpatch=.TRUE.,subpatch_pattern=0_int_8,error=error) 00838 00839 force(ikind)%rho_core(:,iatom) =& 00840 force(ikind)%rho_core(:,iatom) + force_a(:) 00841 00842 IF (use_virial) THEN 00843 virial%pv_virial = virial%pv_virial + my_virial_a 00844 END IF 00845 IF (atenergy) THEN 00846 qs_env%atprop%ateb(iatom) = qs_env%atprop%ateb(iatom) + 0.5_dp*hab(1,1)*pab(1,1) 00847 END IF 00848 00849 END DO 00850 00851 END DO 00852 00853 CALL rs_grid_release(rs_v, error=error) 00854 00855 DEALLOCATE (hab,pab,cores,STAT=stat) 00856 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00857 00858 END IF 00859 00860 CALL timestop(handle) 00861 00862 END SUBROUTINE integrate_v_core_rspace 00863 00864 ! ***************************************************************************** 00880 SUBROUTINE integrate_v_rspace(v_rspace, p, h,qs_env,calculate_forces,compute_tau,gapw,& 00881 basis_set_id, pw_env_external, task_list_external, error) 00882 00883 TYPE(pw_p_type) :: v_rspace 00884 TYPE(cp_dbcsr_p_type), INTENT(IN), 00885 OPTIONAL :: p 00886 TYPE(cp_dbcsr_p_type), INTENT(INOUT) :: h 00887 TYPE(qs_environment_type), POINTER :: qs_env 00888 LOGICAL, INTENT(IN) :: calculate_forces 00889 LOGICAL, INTENT(IN), OPTIONAL :: compute_tau, gapw 00890 INTEGER, INTENT(IN), OPTIONAL :: basis_set_id 00891 TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external 00892 TYPE(task_list_type), OPTIONAL, POINTER :: task_list_external 00893 TYPE(cp_error_type), INTENT(inout) :: error 00894 00895 CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_rspace', 00896 routineP = moduleN//':'//routineN 00897 INTEGER, PARAMETER :: add_tasks = 1000, 00898 max_tasks = 3000 00899 REAL(kind=dp), PARAMETER :: mult_tasks = 2.0_dp 00900 00901 INTEGER :: atom_a, atom_b, bcol, brow, handle, i, iatom, idir, 00902 igrid_level, ikind, ikind_old, ilevel, ipair, ipgf, ipgf_new, iset, 00903 iset_new, iset_old, itask, ithread, jatom, jkind, jkind_old, jpgf, 00904 jpgf_new, jset, jset_new, jset_old, maxco, maxpgf, maxset, maxsgf_set, 00905 my_basis_set_id, na1, na2, natom, nb1, nb2, ncoa, ncob, nkind, nseta, 00906 nsetb, nthread, offs_dv, sgfa, sgfb, stat 00907 INTEGER(KIND=int_8), DIMENSION(:), 00908 POINTER :: atom_pair_recv, atom_pair_send 00909 INTEGER(kind=int_8), DIMENSION(:, :), 00910 POINTER :: tasks 00911 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind 00912 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: block_touched 00913 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, 00914 lb_min, npgfa, npgfb, nsgfa, 00915 nsgfb 00916 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 00917 LOGICAL :: atom_pair_changed, atom_pair_done, distributed_grids, failure, 00918 found, h_duplicated, had_thread_dist, map_consistent, my_compute_tau, 00919 my_gapw, new_set_pair_coming, p_duplicated, pab_required, scatter, 00920 use_subpatch, use_virial 00921 REAL(KIND=dp) :: dab, eps_gvg_rspace, rab2, 00922 zetp 00923 REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra, rab, 00924 rab_inv, rb 00925 REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b 00926 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 00927 REAL(KIND=dp), DIMENSION(:, :), POINTER :: dist_ab, h_block, hab, 00928 p_block, pab, rpgfa, rpgfb, 00929 sphi_a, sphi_b, work, zeta, 00930 zetb 00931 REAL(KIND=dp), DIMENSION(:, :, :), 00932 POINTER :: habt, hadb, hdab, pabt, workt 00933 REAL(kind=dp), DIMENSION(:, :, :, :), 00934 POINTER :: hadbt, hdabt 00935 TYPE(atomic_kind_type), DIMENSION(:), 00936 POINTER :: atomic_kind_set 00937 TYPE(atomic_kind_type), POINTER :: atomic_kind 00938 TYPE(cell_type), POINTER :: cell 00939 TYPE(cp_2d_r_p_type), DIMENSION(3) :: dv_block 00940 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00941 POINTER :: ddv 00942 TYPE(cp_dbcsr_type), POINTER :: deltap, dh 00943 TYPE(cp_para_env_type), POINTER :: para_env 00944 TYPE(cube_info_type), DIMENSION(:), 00945 POINTER :: cube_info 00946 TYPE(dbcsr_distribution_obj), POINTER :: dist 00947 TYPE(dft_control_type), POINTER :: dft_control 00948 TYPE(gridlevel_info_type), POINTER :: gridlevel_info 00949 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 00950 TYPE(particle_type), DIMENSION(:), 00951 POINTER :: particle_set 00952 TYPE(pw_env_type), POINTER :: pw_env 00953 TYPE(qs_force_type), DIMENSION(:), 00954 POINTER :: force 00955 TYPE(realspace_grid_desc_p_type), 00956 DIMENSION(:), POINTER :: rs_descs 00957 TYPE(realspace_grid_p_type), 00958 DIMENSION(:), POINTER :: rs_v 00959 TYPE(section_vals_type), POINTER :: input, interp_section 00960 TYPE(task_list_type), POINTER :: task_list, task_list_soft 00961 TYPE(virial_type), POINTER :: virial 00962 00963 !$ INTEGER :: omp_get_max_threads, omp_get_thread_num 00964 00965 failure=.FALSE. 00966 NULLIFY(pw_env, rs_descs, tasks, dist_ab) 00967 00968 debug_count=debug_count+1 00969 00970 offs_dv=0 00971 my_compute_tau = .FALSE. 00972 my_gapw = .FALSE. 00973 IF (PRESENT(compute_tau)) my_compute_tau = compute_tau 00974 IF (PRESENT(gapw)) my_gapw = gapw 00975 IF (PRESENT(basis_set_id)) THEN 00976 my_basis_set_id = basis_set_id 00977 ELSE 00978 my_basis_set_id = use_orb_basis_set 00979 END IF 00980 00981 SELECT CASE (my_basis_set_id) 00982 CASE (use_orb_basis_set) 00983 CALL get_qs_env(qs_env=qs_env,& 00984 atomic_kind_set=atomic_kind_set,& 00985 cell=cell,& 00986 dft_control=dft_control,& 00987 particle_set=particle_set,& 00988 para_env=para_env,& 00989 input=input,& 00990 task_list=task_list,& 00991 task_list_soft=task_list_soft,& 00992 force=force,pw_env=pw_env,& 00993 virial=virial,error=error) 00994 CASE (use_aux_fit_basis_set) 00995 CALL get_qs_env(qs_env=qs_env,& 00996 atomic_kind_set=atomic_kind_set,& 00997 cell=cell,& 00998 dft_control=dft_control,& 00999 particle_set=particle_set,& 01000 para_env=para_env,& 01001 input=input,& 01002 task_list_aux_fit=task_list,& 01003 task_list_soft=task_list_soft,& 01004 force=force,pw_env=pw_env,& 01005 virial=virial,error=error) 01006 END SELECT 01007 01008 ! maybe don't even pass qs_env ??? 01009 IF (PRESENT(pw_env_external)) pw_env=>pw_env_external 01010 IF (PRESENT(task_list_external)) task_list=>task_list_external 01011 01012 IF (my_compute_tau) THEN 01013 CALL timeset(routineN,handle) 01014 ELSE 01015 CALL timeset(routineN,handle) 01016 END IF 01017 01018 ! get the task lists 01019 IF (my_gapw) task_list=>task_list_soft 01020 CPPrecondition(ASSOCIATED(task_list),cp_failure_level,routineP,error,failure) 01021 tasks =>task_list%tasks 01022 dist_ab=>task_list%dist_ab 01023 atom_pair_send=>task_list%atom_pair_send 01024 atom_pair_recv=>task_list%atom_pair_recv 01025 01026 CPPrecondition(ASSOCIATED(pw_env),cp_failure_level,routineP,error,failure) 01027 CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_v, error=error) 01028 DO i=1,SIZE(rs_v) 01029 CALL rs_grid_retain(rs_v(i)%rs_grid,error=error) 01030 END DO 01031 01032 ! *** assign from pw_env 01033 gridlevel_info=>pw_env%gridlevel_info 01034 cube_info=>pw_env%cube_info 01035 01036 interp_section => section_vals_get_subs_vals(input,"DFT%MGRID%INTERPOLATOR",& 01037 error=error) 01038 CALL potential_pw2rs(rs_v,v_rspace,pw_env,interp_section,error) 01039 01040 ! *** having the potential on the rs_multigrids, just integrate ... 01041 nkind = SIZE(atomic_kind_set) 01042 natom = SIZE(particle_set) 01043 use_virial = virial%pv_availability.AND.(.NOT.virial%pv_numer) 01044 01045 IF (calculate_forces) THEN 01046 ALLOCATE (atom_of_kind(natom),STAT=stat) 01047 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01048 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 01049 atom_of_kind=atom_of_kind) 01050 END IF 01051 01052 map_consistent=dft_control%qs_control%map_consistent 01053 IF (map_consistent) THEN 01054 eps_gvg_rspace = dft_control%qs_control%eps_rho_rspace ! needs to be consistent with rho_rspace 01055 ELSE 01056 eps_gvg_rspace = dft_control%qs_control%eps_gvg_rspace 01057 ENDIF 01058 01059 pab_required = PRESENT(p) .AND. (calculate_forces .OR. .NOT. map_consistent) 01060 01061 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 01062 maxco=maxco,& 01063 maxsgf_set=maxsgf_set,& 01064 basis_set_id=my_basis_set_id) 01065 01066 distributed_grids = .FALSE. 01067 DO igrid_level = 1, gridlevel_info%ngrid_levels 01068 IF ( rs_v(igrid_level)%rs_grid%desc%distributed ) THEN 01069 distributed_grids = .TRUE. 01070 ENDIF 01071 ENDDO 01072 01073 h_duplicated = .FALSE. 01074 dh => h%matrix 01075 IF (distributed_grids) THEN 01076 NULLIFY ( dh ) 01077 h_duplicated = .TRUE. 01078 ALLOCATE(dh) 01079 CALL cp_dbcsr_init(dh, error=error) 01080 CALL cp_dbcsr_create(dh, 'LocalH', & 01081 cp_dbcsr_distribution (h%matrix),& 01082 cp_dbcsr_get_matrix_type (h%matrix), cp_dbcsr_row_block_sizes(h%matrix),& 01083 cp_dbcsr_col_block_sizes(h%matrix), cp_dbcsr_get_num_blocks(h%matrix),& 01084 cp_dbcsr_get_data_size(h%matrix),& 01085 error=error) 01086 01087 END IF 01088 01089 p_duplicated = .FALSE. 01090 IF ( pab_required ) THEN 01091 deltap => p%matrix 01092 IF (distributed_grids) THEN 01093 p_duplicated = .TRUE. 01094 NULLIFY ( deltap ) 01095 ALLOCATE(deltap) 01096 CALL cp_dbcsr_init(deltap, error=error) 01097 CALL cp_dbcsr_copy(deltap,p%matrix,name="LocalP",error=error) 01098 END IF 01099 END IF 01100 01101 nthread = 1 01102 !$ nthread = omp_get_max_threads() 01103 01104 ! *** Allocate work storage *** 01105 NULLIFY ( pabt, habt, workt ) 01106 CALL reallocate(habt,1,maxco,1,maxco,0,nthread) 01107 CALL reallocate(workt,1,maxco,1,maxsgf_set,0,nthread) 01108 IF (pab_required) THEN 01109 CALL reallocate(pabt,1,maxco,1,maxco,0,nthread) 01110 ELSE 01111 IF (calculate_forces) THEN 01112 CALL stop_program(routineN,moduleN,__LINE__,& 01113 "Need p for forces") 01114 END IF 01115 END IF 01116 01117 NULLIFY(hdabt,hadbt,hdab,hadb) 01118 01119 ! get maximum numbers 01120 natom = SIZE( particle_set ) 01121 maxset=0 01122 maxpgf=0 01123 DO ikind=1,nkind 01124 atomic_kind => atomic_kind_set(ikind) 01125 SELECT CASE (my_basis_set_id) 01126 CASE (use_orb_basis_set) 01127 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01128 softb = my_gapw, & 01129 orb_basis_set=orb_basis_set) 01130 CASE (use_aux_fit_basis_set) 01131 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01132 softb = my_gapw, & 01133 aux_fit_basis_set=orb_basis_set,& 01134 basis_set_id=my_basis_set_id) 01135 END SELECT 01136 01137 IF (.NOT.ASSOCIATED(orb_basis_set)) CYCLE 01138 01139 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 01140 npgf=npgfa, nset=nseta ) 01141 01142 maxset=MAX(nseta,maxset) 01143 maxpgf=MAX(MAXVAL(npgfa),maxpgf) 01144 END DO 01145 01146 IF (distributed_grids .AND. pab_required) THEN 01147 CALL rs_distribute_matrix (rs_descs, deltap, atom_pair_send, atom_pair_recv, natom, scatter=.TRUE., error=error) 01148 ENDIF 01149 01150 IF (debug_this_module) & 01151 ALLOCATE(block_touched(natom,natom)) 01152 01153 !$omp parallel default(none), & 01154 !$omp shared(workt,habt,hdabt,hadbt,pabt,tasks,particle_set,natom,maxset), & 01155 !$omp shared(maxpgf,my_basis_set_id,my_gapw,dh,ddv,deltap,use_virial), & 01156 !$omp shared(pab_required,calculate_forces,ncoset,rs_v,cube_info,my_compute_tau), & 01157 !$omp shared(map_consistent,eps_gvg_rspace,force,virial,cell,atom_of_kind,dist_ab), & 01158 !$omp shared(gridlevel_info,task_list,failure,block_touched,nthread), & 01159 !$omp private(ithread,work,hab,hdab,hadb,pab,iset_old,jset_old), & 01160 !$omp private(ikind_old,jkind_old,iatom,jatom,iset,jset,ikind,jkind,ilevel,ipgf,jpgf), & 01161 !$omp private(brow,bcol,orb_basis_set,first_sgfa,la_max,la_min,npgfa,nseta,nsgfa), & 01162 !$omp private(rpgfa,set_radius_a,sphi_a,zeta,first_sgfb,lb_max,lb_min,npgfb), & 01163 !$omp private(nsetb,nsgfb,rpgfb,set_radius_b,sphi_b,zetb,found,error,atom_a,atom_b), & 01164 !$omp private(force_a,force_b,my_virial_a,my_virial_b,atom_pair_changed,h_block), & 01165 !$omp private(dv_block,p_block,ncoa,sgfa,ncob,sgfb,rab,rab2,ra,rb,zetp,dab,igrid_level), & 01166 !$omp private(na1,na2,nb1,nb2,use_subpatch,rab_inv,new_set_pair_coming,atom_pair_done), & 01167 !$omp private(iset_new,jset_new,ipgf_new,jpgf_new,idir,dist), & 01168 !$omp private(had_thread_dist,itask) 01169 01170 ithread = 0 01171 !$ ithread = omp_get_thread_num() 01172 work => workt(:,:,ithread) 01173 hab => habt(:,:,ithread) 01174 IF (pab_required) THEN 01175 pab => pabt(:,:,ithread) 01176 END IF 01177 01178 iset_old = -1 ; jset_old = -1 01179 ikind_old = -1 ; jkind_old = -1 01180 01181 01182 ! Here we loop over gridlevels first, finalising the matrix after each grid level is 01183 ! completed. On each grid level, we loop over atom pairs, which will only access 01184 ! a single block of each matrix, so with OpenMP, each matrix block is only touched 01185 ! by a single thread for each grid level 01186 loop_gridlevels: DO igrid_level = 1, gridlevel_info%ngrid_levels 01187 01188 CALL cp_dbcsr_work_create(dh,work_mutable=.TRUE.,n=nthread,error=error) 01189 !$ dist => dh%matrix%m%dist 01190 !$ CALL cp_assert (dbcsr_distribution_has_threads(dist), cp_fatal_level,& 01191 !$ cp_internal_error, routineN, "No thread distribution defined.",& 01192 !$ error=error) 01193 !$omp barrier 01194 01195 IF (debug_this_module) THEN 01196 !$omp single 01197 block_touched = -1 01198 !$omp end single 01199 !$omp flush 01200 END IF 01201 01202 !$omp do schedule (dynamic, MAX(1,task_list%npairs(igrid_level)/(nthread*50))) 01203 loop_pairs: DO ipair = 1, task_list%npairs(igrid_level) 01204 loop_tasks: DO itask = task_list%taskstart(ipair,igrid_level), task_list%taskstop(ipair,igrid_level) 01205 01206 CALL int2pair(tasks(3,itask),ilevel,iatom,jatom,iset,jset,ipgf,jpgf,natom,maxset,maxpgf) 01207 01208 ! At the start of a block of tasks, get atom data (and kind data, if needed) 01209 IF (itask .EQ. task_list%taskstart(ipair,igrid_level) ) THEN 01210 01211 ikind = particle_set(iatom)%atomic_kind%kind_number 01212 jkind = particle_set(jatom)%atomic_kind%kind_number 01213 01214 ra(:) = pbc(particle_set(iatom)%r,cell) 01215 01216 IF (iatom <= jatom) THEN 01217 brow = iatom 01218 bcol = jatom 01219 ELSE 01220 brow = jatom 01221 bcol = iatom 01222 END IF 01223 01224 IF (ikind .NE. ikind_old ) THEN 01225 SELECT CASE (my_basis_set_id) 01226 CASE (use_orb_basis_set) 01227 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& 01228 softb = my_gapw, & 01229 orb_basis_set=orb_basis_set) 01230 CASE (use_aux_fit_basis_set) 01231 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& 01232 softb = my_gapw, & 01233 aux_fit_basis_set=orb_basis_set,& 01234 basis_set_id = my_basis_set_id) 01235 END SELECT 01236 01237 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 01238 first_sgf=first_sgfa,& 01239 lmax=la_max,& 01240 lmin=la_min,& 01241 npgf=npgfa,& 01242 nset=nseta,& 01243 nsgf_set=nsgfa,& 01244 pgf_radius=rpgfa,& 01245 set_radius=set_radius_a,& 01246 sphi=sphi_a,& 01247 zet=zeta) 01248 ENDIF 01249 01250 IF (jkind .NE. jkind_old ) THEN 01251 SELECT CASE (my_basis_set_id) 01252 CASE (use_orb_basis_set) 01253 CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind,& 01254 softb = my_gapw, & 01255 orb_basis_set=orb_basis_set) 01256 CASE (use_aux_fit_basis_set) 01257 CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind,& 01258 softb = my_gapw, & 01259 aux_fit_basis_set=orb_basis_set,& 01260 basis_set_id=my_basis_set_id) 01261 END SELECT 01262 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 01263 first_sgf=first_sgfb,& 01264 lmax=lb_max,& 01265 lmin=lb_min,& 01266 npgf=npgfb,& 01267 nset=nsetb,& 01268 nsgf_set=nsgfb,& 01269 pgf_radius=rpgfb,& 01270 set_radius=set_radius_b,& 01271 sphi=sphi_b,& 01272 zet=zetb) 01273 01274 ENDIF 01275 01276 IF (debug_this_module) THEN 01277 !$omp critical (block_touched_critical) 01278 IF ((block_touched(brow,bcol).NE.ithread) .AND. (block_touched(brow,bcol).NE. -1) ) THEN 01279 CALL stop_program(routineN,moduleN,__LINE__,& 01280 "Block has been modified by another thread") 01281 END IF 01282 block_touched(brow,bcol) = ithread 01283 !$omp end critical (block_touched_critical) 01284 END IF 01285 01286 NULLIFY(h_block) 01287 CALL cp_dbcsr_get_block_p(dh,brow,bcol,h_block,found) 01288 IF (.NOT.ASSOCIATED(h_block)) THEN 01289 CALL cp_dbcsr_add_block_node ( dh, brow, bcol, h_block ,error=error) 01290 END IF 01291 01292 IF (pab_required) THEN 01293 CALL cp_dbcsr_get_block_p(matrix=deltap,& 01294 row=brow,col=bcol,BLOCK=p_block,found=found) 01295 01296 IF (.NOT.ASSOCIATED(p_block)) THEN 01297 CALL stop_program(routineN,moduleN,__LINE__,& 01298 "p_block not associated in deltap") 01299 END IF 01300 END IF 01301 01302 IF (calculate_forces) THEN 01303 atom_a = atom_of_kind(iatom) 01304 atom_b = atom_of_kind(jatom) 01305 force_a(:) = 0.0_dp 01306 force_b(:) = 0.0_dp 01307 ENDIF 01308 IF (use_virial) THEN 01309 my_virial_a = 0.0_dp 01310 my_virial_b = 0.0_dp 01311 ENDIF 01312 01313 ikind_old = ikind 01314 jkind_old = jkind 01315 01316 atom_pair_changed = .TRUE. 01317 01318 ELSE 01319 01320 atom_pair_changed = .FALSE. 01321 01322 ENDIF 01323 01324 IF (atom_pair_changed .OR. iset_old .NE. iset .OR. jset_old .NE. jset) THEN 01325 01326 ncoa = npgfa(iset)*ncoset(la_max(iset)) 01327 sgfa = first_sgfa(1,iset) 01328 ncob = npgfb(jset)*ncoset(lb_max(jset)) 01329 sgfb = first_sgfb(1,jset) 01330 IF (pab_required) THEN 01331 IF (iatom <= jatom) THEN 01332 CALL dgemm("N","N",ncoa,nsgfb(jset),nsgfa(iset),& 01333 1.0_dp,sphi_a(1,sgfa),SIZE(sphi_a,1),& 01334 p_block(sgfa,sgfb),SIZE(p_block,1),& 01335 0.0_dp,work(1,1),SIZE(work,1)) 01336 CALL dgemm("N","T",ncoa,ncob,nsgfb(jset),& 01337 1.0_dp,work(1,1),SIZE(work,1),& 01338 sphi_b(1,sgfb),SIZE(sphi_b,1),& 01339 0.0_dp,pab(1,1),SIZE(pab,1)) 01340 ELSE 01341 CALL dgemm("N","N",ncob,nsgfa(iset),nsgfb(jset),& 01342 1.0_dp,sphi_b(1,sgfb),SIZE(sphi_b,1),& 01343 p_block(sgfb,sgfa),SIZE(p_block,1),& 01344 0.0_dp,work(1,1),SIZE(work,1)) 01345 CALL dgemm("N","T",ncob,ncoa,nsgfa(iset),& 01346 1.0_dp,work(1,1),SIZE(work,1),& 01347 sphi_a(1,sgfa),SIZE(sphi_a,1),& 01348 0.0_dp,pab(1,1),SIZE(pab,1)) 01349 END IF 01350 END IF 01351 01352 IF (iatom<=jatom) THEN 01353 hab(1:ncoa,1:ncob) = 0._dp 01354 ELSE 01355 hab(1:ncob,1:ncoa) = 0._dp 01356 ENDIF 01357 01358 iset_old = iset 01359 jset_old = jset 01360 01361 ENDIF 01362 01363 rab(1) = dist_ab (1,itask) 01364 rab(2) = dist_ab (2,itask) 01365 rab(3) = dist_ab (3,itask) 01366 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 01367 rb(1) = ra(1) + rab(1) 01368 rb(2) = ra(2) + rab(2) 01369 rb(3) = ra(3) + rab(3) 01370 zetp = zeta(ipgf,iset) + zetb(jpgf,jset) 01371 dab=SQRT(rab2) 01372 01373 na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1 01374 na2 = ipgf*ncoset(la_max(iset)) 01375 nb1 = (jpgf - 1)*ncoset(lb_max(jset)) + 1 01376 nb2 = jpgf*ncoset(lb_max(jset)) 01377 01378 ! check whether we need to use fawzi's generalised collocation scheme 01379 IF(rs_v(igrid_level)%rs_grid%desc%distributed)THEN 01380 !tasks(4,:) is 0 for replicated, 1 for distributed 2 for exceptional distributed tasks 01381 IF (tasks(4,itask) .EQ. 2 ) THEN 01382 use_subpatch = .TRUE. 01383 ELSE 01384 use_subpatch = .FALSE. 01385 ENDIF 01386 ELSE 01387 use_subpatch = .FALSE. 01388 ENDIF 01389 01390 IF (pab_required) THEN 01391 IF (iatom <= jatom) THEN 01392 CALL integrate_pgf_product_rspace(& 01393 la_max(iset),zeta(ipgf,iset),la_min(iset),& 01394 lb_max(jset),zetb(jpgf,jset),lb_min(jset),& 01395 ra,rab,rab2,rs_v(igrid_level)%rs_grid,cell,& 01396 cube_info(igrid_level),& 01397 hab,pab=pab,o1=na1-1,o2=nb1-1, & 01398 eps_gvg_rspace=eps_gvg_rspace,& 01399 calculate_forces=calculate_forces,& 01400 force_a=force_a,force_b=force_b,& 01401 compute_tau=my_compute_tau,map_consistent=map_consistent,& 01402 use_virial=use_virial,my_virial_a=my_virial_a,& 01403 my_virial_b=my_virial_b,use_subpatch=use_subpatch,subpatch_pattern=tasks(6,itask),error=error) 01404 ELSE 01405 rab_inv=-rab 01406 CALL integrate_pgf_product_rspace(& 01407 lb_max(jset),zetb(jpgf,jset),lb_min(jset),& 01408 la_max(iset),zeta(ipgf,iset),la_min(iset),& 01409 rb,rab_inv,rab2,rs_v(igrid_level)%rs_grid,cell,& 01410 cube_info(igrid_level),& 01411 hab,pab=pab,o1=nb1-1,o2=na1-1, & 01412 eps_gvg_rspace=eps_gvg_rspace,& 01413 calculate_forces=calculate_forces,& 01414 force_a=force_b,force_b=force_a,& 01415 compute_tau=my_compute_tau,map_consistent=map_consistent,& 01416 use_virial=use_virial,my_virial_a=my_virial_b,& 01417 my_virial_b=my_virial_a,use_subpatch=use_subpatch,subpatch_pattern=tasks(6,itask),error=error) 01418 END IF 01419 ELSE 01420 IF (iatom <= jatom) THEN 01421 CALL integrate_pgf_product_rspace(& 01422 la_max(iset),zeta(ipgf,iset),la_min(iset),& 01423 lb_max(jset),zetb(jpgf,jset),lb_min(jset),& 01424 ra,rab,rab2,rs_v(igrid_level)%rs_grid,cell,& 01425 cube_info(igrid_level),& 01426 hab,o1=na1-1,o2=nb1-1,& 01427 eps_gvg_rspace=eps_gvg_rspace,& 01428 calculate_forces=calculate_forces,& 01429 force_a=force_a,force_b=force_b,& 01430 compute_tau=my_compute_tau,& 01431 map_consistent=map_consistent,use_subpatch=use_subpatch,subpatch_pattern=tasks(6,itask),error=error) 01432 ELSE 01433 rab_inv=-rab 01434 CALL integrate_pgf_product_rspace(& 01435 lb_max(jset),zetb(jpgf,jset),lb_min(jset),& 01436 la_max(iset),zeta(ipgf,iset),la_min(iset),& 01437 rb,rab_inv,rab2,rs_v(igrid_level)%rs_grid,cell,& 01438 cube_info(igrid_level),& 01439 hab,o1=nb1-1,o2=na1-1,& 01440 eps_gvg_rspace=eps_gvg_rspace,& 01441 calculate_forces=calculate_forces,& 01442 force_a=force_b,force_b=force_a, & 01443 compute_tau=my_compute_tau,& 01444 map_consistent=map_consistent,use_subpatch=use_subpatch,subpatch_pattern=tasks(6,itask),error=error) 01445 END IF 01446 END IF 01447 01448 new_set_pair_coming=.FALSE. 01449 atom_pair_done = .FALSE. 01450 IF (itask < task_list%taskstop(ipair,igrid_level)) THEN 01451 CALL int2pair(tasks(3,itask+1),ilevel,iatom,jatom,iset_new,jset_new,ipgf_new,jpgf_new,natom,maxset,maxpgf) 01452 IF (iset_new .NE. iset .OR. jset_new .NE. jset) THEN 01453 new_set_pair_coming=.TRUE. 01454 ENDIF 01455 ELSE 01456 ! do not forget the last block 01457 new_set_pair_coming=.TRUE. 01458 atom_pair_done = .TRUE. 01459 ENDIF 01460 01461 ! contract the block into h if we're done with the current set pair 01462 IF (new_set_pair_coming) THEN 01463 IF (iatom <= jatom) THEN 01464 CALL dgemm("N","N",ncoa,nsgfb(jset),ncob,& 01465 1.0_dp,hab(1,1),SIZE(hab,1),& 01466 sphi_b(1,sgfb),SIZE(sphi_b,1),& 01467 0.0_dp,work(1,1),SIZE(work,1)) 01468 CALL dgemm("T","N",nsgfa(iset),nsgfb(jset),ncoa,& 01469 1.0_dp,sphi_a(1,sgfa),SIZE(sphi_a,1),& 01470 work(1,1),SIZE(work,1),& 01471 1.0_dp,h_block(sgfa,sgfb),SIZE(h_block,1)) 01472 ELSE 01473 CALL dgemm("N","N",ncob,nsgfa(iset),ncoa,& 01474 1.0_dp,hab(1,1),SIZE(hab,1),& 01475 sphi_a(1,sgfa),SIZE(sphi_a,1),& 01476 0.0_dp,work(1,1),SIZE(work,1)) 01477 CALL dgemm("T","N",nsgfb(jset),nsgfa(iset),ncob,& 01478 1.0_dp,sphi_b(1,sgfb),SIZE(sphi_b,1),& 01479 work(1,1),SIZE(work,1),& 01480 1.0_dp,h_block(sgfb,sgfa),SIZE(h_block,1)) 01481 END IF 01482 END IF 01483 01484 IF (atom_pair_done) THEN 01485 !$omp critical(force_critical) 01486 IF (calculate_forces) THEN 01487 force(ikind)%rho_elec(:,atom_a) =& 01488 force(ikind)%rho_elec(:,atom_a) + 2.0_dp*force_a(:) 01489 IF (iatom /= jatom ) THEN 01490 force(jkind)%rho_elec(:,atom_b) =& 01491 force(jkind)%rho_elec(:,atom_b) + 2.0_dp*force_b(:) 01492 END IF 01493 ENDIF 01494 IF (use_virial) THEN 01495 IF (use_virial .AND. calculate_forces) THEN 01496 virial%pv_virial = virial%pv_virial + 2.0_dp*my_virial_a 01497 IF (iatom /= jatom) THEN 01498 virial%pv_virial = virial%pv_virial + 2.0_dp*my_virial_b 01499 END IF 01500 END IF 01501 END IF 01502 !$omp end critical (force_critical) 01503 ENDIF 01504 END DO loop_tasks 01505 END DO loop_pairs 01506 !$omp end do 01507 01508 CALL cp_dbcsr_finalize(dh, error=error) 01509 01510 END DO loop_gridlevels 01511 01512 !$omp end parallel 01513 01514 IF (debug_this_module) & 01515 DEALLOCATE(block_touched) 01516 01517 IF ( h_duplicated ) THEN 01518 ! Reconstruct H matrix if using distributed RS grids 01519 ! note send and recv direction reversed WRT collocate 01520 scatter = .FALSE. 01521 CALL rs_distribute_matrix (rs_descs, dh, atom_pair_recv, atom_pair_send,& 01522 natom, scatter, error, h%matrix) 01523 CALL cp_dbcsr_deallocate_matrix ( dh ,error=error) 01524 01525 ELSE 01526 NULLIFY ( dh, ddv ) 01527 END IF 01528 01529 IF ( p_duplicated ) THEN 01530 CALL cp_dbcsr_deallocate_matrix ( deltap ,error=error) 01531 ELSE 01532 NULLIFY ( deltap ) 01533 END IF 01534 01535 ! *** Release work storage *** 01536 01537 DEALLOCATE (habt,workt,STAT=stat) 01538 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01539 01540 IF ( pab_required ) THEN 01541 DEALLOCATE (pabt,STAT=stat) 01542 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01543 END IF 01544 01545 IF (ASSOCIATED(rs_v)) THEN 01546 DO i=1,SIZE(rs_v) 01547 CALL rs_grid_release(rs_v(i)%rs_grid, error=error) 01548 END DO 01549 END IF 01550 01551 IF (calculate_forces) THEN 01552 DEALLOCATE (atom_of_kind,STAT=stat) 01553 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01554 END IF 01555 01556 CALL timestop(handle) 01557 01558 END SUBROUTINE integrate_v_rspace 01559 01560 ! ***************************************************************************** 01572 SUBROUTINE potential_pw2rs(rs_v,v_rspace,pw_env,interp_section,error) 01573 01574 TYPE(realspace_grid_p_type), 01575 DIMENSION(:), POINTER :: rs_v 01576 TYPE(pw_p_type), INTENT(IN) :: v_rspace 01577 TYPE(pw_env_type), POINTER :: pw_env 01578 TYPE(section_vals_type), POINTER :: interp_section 01579 TYPE(cp_error_type), INTENT(inout) :: error 01580 01581 CHARACTER(len=*), PARAMETER :: routineN = 'potential_pw2rs', 01582 routineP = moduleN//':'//routineN 01583 01584 INTEGER :: auxbas_grid, handle, 01585 igrid_level, interp_kind 01586 REAL(KIND=dp) :: scale 01587 TYPE(gridlevel_info_type), POINTER :: gridlevel_info 01588 TYPE(pw_p_type), DIMENSION(:), POINTER :: mgrid_gspace, mgrid_rspace 01589 TYPE(pw_pool_p_type), DIMENSION(:), 01590 POINTER :: pw_pools 01591 01592 CALL timeset(routineN,handle) 01593 01594 ! *** set up of the potential on the multigrids 01595 CALL pw_env_get(pw_env, pw_pools=pw_pools, gridlevel_info=gridlevel_info, & 01596 auxbas_grid = auxbas_grid, error=error) 01597 01598 CALL pw_pools_create_pws(pw_pools,mgrid_rspace,& 01599 use_data = REALDATA3D,& 01600 in_space = REALSPACE, error=error) 01601 01602 ! use either realspace or fft techniques to get the potential on the rs multigrids 01603 CALL section_vals_val_get(interp_section,"KIND",i_val=interp_kind,error=error) 01604 SELECT CASE(interp_kind) 01605 CASE (pw_interp) 01606 CALL pw_pools_create_pws(pw_pools,mgrid_gspace,& 01607 use_data = COMPLEXDATA1D,& 01608 in_space = RECIPROCALSPACE, error=error) 01609 CALL pw_transfer(v_rspace%pw,mgrid_gspace(auxbas_grid)%pw,error=error) 01610 DO igrid_level=1,gridlevel_info%ngrid_levels 01611 IF ( igrid_level /= auxbas_grid ) THEN 01612 CALL pw_copy(mgrid_gspace(auxbas_grid)%pw,mgrid_gspace(igrid_level)%pw,& 01613 error=error) 01614 CALL pw_transfer(mgrid_gspace(igrid_level)%pw,mgrid_rspace(igrid_level)%pw,& 01615 error=error) 01616 ELSE 01617 IF (mgrid_gspace(auxbas_grid)%pw%pw_grid%spherical) THEN 01618 CALL pw_transfer(mgrid_gspace(auxbas_grid)%pw,mgrid_rspace(auxbas_grid)%pw,& 01619 error=error) 01620 ELSE ! fft forward + backward should be identical 01621 CALL pw_copy(v_rspace%pw,mgrid_rspace(auxbas_grid)%pw,error=error) 01622 ENDIF 01623 ENDIF 01624 ! *** Multiply by the grid volume element ratio *** 01625 IF ( igrid_level /= auxbas_grid ) THEN 01626 scale = mgrid_rspace(igrid_level)%pw%pw_grid%dvol/& 01627 mgrid_rspace(auxbas_grid)%pw%pw_grid%dvol 01628 mgrid_rspace(igrid_level)%pw%cr3d = & 01629 scale*mgrid_rspace(igrid_level)%pw%cr3d 01630 END IF 01631 END DO 01632 CALL pw_pools_give_back_pws(pw_pools,mgrid_gspace,error=error) 01633 CASE(spline3_pbc_interp) 01634 CALL pw_copy(v_rspace%pw,mgrid_rspace(1)%pw,error=error) 01635 DO igrid_level=1,gridlevel_info%ngrid_levels-1 01636 CALL pw_zero(mgrid_rspace(igrid_level+1)%pw,error=error) 01637 CALL pw_restrict_s3(mgrid_rspace(igrid_level)%pw,& 01638 mgrid_rspace(igrid_level+1)%pw,pw_pools(igrid_level+1)%pool,& 01639 interp_section,error=error) 01640 ! *** Multiply by the grid volume element ratio 01641 mgrid_rspace(igrid_level+1) % pw % cr3d = & 01642 mgrid_rspace(igrid_level+1) % pw % cr3d * 8._dp 01643 END DO 01644 CASE default 01645 CALL cp_unimplemented_error(routineN,"interpolation not supported "//& 01646 cp_to_string(interp_kind),error=error) 01647 END SELECT 01648 01649 DO igrid_level=1,gridlevel_info%ngrid_levels 01650 CALL rs_pw_transfer(rs_v(igrid_level)%rs_grid,& 01651 mgrid_rspace(igrid_level)%pw,pw2rs,error=error) 01652 ENDDO 01653 ! *** give back the pw multi-grids 01654 CALL pw_pools_give_back_pws(pw_pools,mgrid_rspace,error=error) 01655 01656 CALL timestop(handle) 01657 01658 END SUBROUTINE potential_pw2rs 01659 01660 END MODULE qs_integrate_potential
1.7.3