|
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 ! ***************************************************************************** 00034 MODULE qs_collocate_density 00035 USE atomic_kind_types, ONLY: atomic_kind_type,& 00036 get_atomic_kind,& 00037 get_atomic_kind_set 00038 USE basis_set_types, ONLY: get_gto_basis_set,& 00039 gto_basis_set_type 00040 USE cell_types, ONLY: cell_type,& 00041 pbc 00042 USE cp_control_types, ONLY: dft_control_type 00043 USE cp_dbcsr_interface, ONLY: cp_dbcsr_copy,& 00044 cp_dbcsr_get_block_p,& 00045 cp_dbcsr_init 00046 USE cp_dbcsr_operations, ONLY: cp_dbcsr_deallocate_matrix 00047 USE cp_dbcsr_types, ONLY: cp_dbcsr_type 00048 USE cp_fm_types, ONLY: cp_fm_get_element,& 00049 cp_fm_get_info,& 00050 cp_fm_type 00051 USE cp_para_types, ONLY: cp_para_env_type 00052 USE cube_utils, ONLY: compute_cube_center,& 00053 cube_info_type,& 00054 return_cube,& 00055 return_cube_nonortho 00056 USE d3_poly, ONLY: poly_cp2k2d3 00057 USE external_potential_types, ONLY: get_potential,& 00058 gth_potential_type 00059 USE gauss_colloc, ONLY: collocGauss 00060 USE gaussian_gridlevels, ONLY: gaussian_gridlevel,& 00061 gridlevel_info_type 00062 USE input_constants, ONLY: pw_interp,& 00063 spline3_pbc_interp,& 00064 use_aux_basis_set,& 00065 use_aux_fit_basis_set,& 00066 use_orb_basis_set,& 00067 use_ri_aux_basis_set 00068 USE input_section_types, ONLY: section_vals_get_subs_vals,& 00069 section_vals_type,& 00070 section_vals_val_get 00071 USE kinds, ONLY: dp,& 00072 dp_size,& 00073 int_8 00074 USE lgrid_types, ONLY: lgrid_allocate_grid,& 00075 lgrid_type 00076 USE mathconstants, ONLY: dfac,& 00077 fac,& 00078 pi,& 00079 twopi 00080 USE memory_utilities, ONLY: reallocate 00081 USE orbital_pointers, ONLY: coset,& 00082 indco,& 00083 nco,& 00084 ncoset,& 00085 nso,& 00086 nsoset 00087 USE orbital_transformation_matrices, ONLY: orbtramat 00088 USE particle_types, ONLY: particle_type 00089 USE pw_env_types, ONLY: pw_env_get,& 00090 pw_env_type 00091 USE pw_methods, ONLY: pw_axpy,& 00092 pw_copy,& 00093 pw_integrate_function,& 00094 pw_transfer,& 00095 pw_zero 00096 USE pw_pool_types, ONLY: pw_pool_create_pw,& 00097 pw_pool_give_back_pw,& 00098 pw_pool_p_type,& 00099 pw_pool_type,& 00100 pw_pools_create_pws,& 00101 pw_pools_give_back_pws 00102 USE pw_spline_utils, ONLY: pw_prolongate_s3 00103 USE pw_types, ONLY: COMPLEXDATA1D,& 00104 REALDATA3D,& 00105 REALSPACE,& 00106 RECIPROCALSPACE,& 00107 pw_p_type,& 00108 pw_type 00109 USE qs_environment_types, ONLY: get_qs_env,& 00110 qs_environment_type 00111 USE qs_interactions, ONLY: exp_radius_very_extended 00112 USE qs_modify_pab_block, ONLY: & 00113 FUNC_AB, FUNC_ADBmDAB, FUNC_ARB, FUNC_ARDBmDARB, FUNC_DABpADB, & 00114 FUNC_DADB, FUNC_DER, FUNC_DX, FUNC_DXDX, FUNC_DXDY, FUNC_DY, & 00115 FUNC_DYDY, FUNC_DYDZ, FUNC_DZ, FUNC_DZDX, FUNC_DZDZ, & 00116 prepare_adb_m_dab, prepare_arb, prepare_ardb_m_darb, & 00117 prepare_dab_p_adb, prepare_dadb, prepare_diadib, prepare_diiadiib, & 00118 prepare_dijadijb 00119 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type 00120 USE qs_rho_types, ONLY: qs_rho_type 00121 USE realspace_grid_types, ONLY: & 00122 realspace_grid_desc_p_type, realspace_grid_desc_type, & 00123 realspace_grid_p_type, realspace_grid_type, rs2pw, rs_grid_release, & 00124 rs_grid_retain, rs_grid_zero, rs_pw_transfer 00125 USE scptb_types, ONLY: get_scptb_parameter,& 00126 scp_vector_type,& 00127 scptb_parameter_type 00128 USE task_list_methods, ONLY: int2pair,& 00129 rs_distribute_matrix 00130 USE task_list_types, ONLY: task_list_type 00131 USE termination, ONLY: stop_memory,& 00132 stop_program 00133 USE timings, ONLY: timeset,& 00134 timestop 00135 #include "cp_common_uses.h" 00136 00137 IMPLICIT NONE 00138 00139 PRIVATE 00140 00141 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_collocate_density' 00142 ! *** Public subroutines *** 00143 00144 PUBLIC :: calculate_ppl_grid,& 00145 calculate_rho_core,& 00146 calculate_rho_single_gaussian,& 00147 calculate_rho_metal,& 00148 calculate_rho_resp_single,& 00149 calculate_rho_resp_all,& 00150 calculate_rho_elec,& 00151 calculate_drho_elec,& 00152 calculate_scp_charge,& 00153 calculate_wavefunction,& 00154 collocate_pgf_product_gspace,& 00155 collocate_pgf_product_rspace,& 00156 collocate_atomic_charge_density,& 00157 density_rs2pw,& 00158 calculate_rho_nlcc 00159 00160 INTEGER :: debug_count=0 00161 00162 CONTAINS 00163 00164 ! ***************************************************************************** 00167 SUBROUTINE calculate_scp_charge(scp,qs_env,scpv,error) 00168 00169 TYPE(pw_p_type), INTENT(INOUT) :: scp 00170 TYPE(qs_environment_type), POINTER :: qs_env 00171 TYPE(scp_vector_type), POINTER :: scpv 00172 TYPE(cp_error_type), INTENT(INOUT) :: error 00173 00174 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_scp_charge', 00175 routineP = moduleN//':'//routineN 00176 00177 INTEGER :: atom_a, handle, i, iatom, 00178 ierr, ii, ikind, ithread, j, 00179 jj, l, lmaxscp, natom, ni, 00180 nj, npme, nthread 00181 INTEGER(KIND=int_8) :: subpatch_pattern 00182 INTEGER, DIMENSION(:), POINTER :: atom_list, cores 00183 LOGICAL :: defined, failure 00184 REAL(KIND=dp) :: alpha, eps_rho_rspace, norm, 00185 pp 00186 REAL(KIND=dp), DIMENSION(3) :: ra 00187 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab 00188 TYPE(atomic_kind_type), DIMENSION(:), 00189 POINTER :: atomic_kind_set 00190 TYPE(atomic_kind_type), POINTER :: atomic_kind 00191 TYPE(cell_type), POINTER :: cell 00192 TYPE(cp_para_env_type), POINTER :: para_env 00193 TYPE(cube_info_type) :: cube_info 00194 TYPE(dft_control_type), POINTER :: dft_control 00195 TYPE(particle_type), DIMENSION(:), 00196 POINTER :: particle_set 00197 TYPE(pw_env_type), POINTER :: pw_env 00198 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00199 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 00200 TYPE(realspace_grid_type), POINTER :: rs_rho 00201 TYPE(scptb_parameter_type), POINTER :: scptb_kind 00202 00203 CALL timeset(routineN,handle) 00204 00205 failure = .FALSE. 00206 00207 NULLIFY(atomic_kind,cell,dft_control,pab,atomic_kind_set,particle_set,& 00208 atom_list,pw_env,rs_rho,auxbas_rs_desc,auxbas_pw_pool,cores) 00209 00210 CALL get_qs_env(qs_env=qs_env,& 00211 atomic_kind_set=atomic_kind_set,& 00212 dft_control=dft_control,& 00213 cell=cell,& 00214 particle_set=particle_set,& 00215 para_env=para_env,pw_env=pw_env,error=error) 00216 CALL pw_env_get(pw_env,auxbas_rs_desc=auxbas_rs_desc,auxbas_rs_grid=rs_rho,& 00217 auxbas_pw_pool=auxbas_pw_pool,error=error) 00218 cube_info=pw_env%cube_info(1) 00219 CALL rs_grid_retain(rs_rho, error=error) 00220 CALL rs_grid_zero(rs_rho) 00221 00222 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 00223 00224 DO ikind=1,SIZE(atomic_kind_set) 00225 00226 atomic_kind => atomic_kind_set(ikind) 00227 00228 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00229 natom=natom,atom_list=atom_list,& 00230 scptb_parameter=scptb_kind) 00231 CALL get_scptb_parameter(scptb_kind,defined=defined,lmaxscp=lmaxscp,ag=alpha) 00232 00233 IF (.NOT.defined) CYCLE 00234 00235 ni = ncoset(lmaxscp) 00236 ALLOCATE(pab(ni,1),STAT=ierr) 00237 CPPrecondition(ierr==0,cp_failure_level,routineP,error,failure) 00238 00239 nthread = 1 00240 ithread=0 00241 00242 ALLOCATE(cores(natom),STAT=ierr) 00243 CPPrecondition(ierr==0,cp_failure_level,routineP,error,failure) 00244 cores = 0 00245 npme = 0 00246 00247 DO iatom = 1, natom 00248 atom_a = atom_list(iatom) 00249 ra(:) = pbc(particle_set(atom_a)%r,cell) 00250 IF(rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN 00251 ! replicated realspace grid, split the atoms up between procs 00252 IF (MODULO(iatom,rs_rho%desc%group_size) == rs_rho % desc % my_pos ) THEN 00253 npme = npme + 1 00254 cores (npme) = iatom 00255 ENDIF 00256 ELSE 00257 npme = npme + 1 00258 cores (npme) = iatom 00259 ENDIF 00260 END DO 00261 00262 DO nj=1,npme 00263 00264 iatom = cores(nj) 00265 atom_a = atom_list(iatom) 00266 ra(:) = pbc(particle_set(atom_a)%r,cell) 00267 subpatch_pattern=0 00268 pab = 0._dp 00269 DO l=0,lmaxscp 00270 pp = (2._dp*l+3._dp)/2._dp 00271 norm = 2._dp**(l+2)/SQRT(pi)/dfac(2*l+1) 00272 norm = SQRT(0.25_dp*dfac(2*l+1)/pi) * norm*alpha**pp 00273 DO jj=1,nco(l) 00274 j = ncoset(l-1) + jj 00275 DO ii=1,nso(l) 00276 i = nsoset(l-1) + ii 00277 pab(j,1) = pab(j,1) + scpv%vector(ikind)%vmat(i,iatom)*& 00278 orbtramat(l)%s2c(ii,jj)*norm 00279 END DO 00280 END DO 00281 END DO 00282 CALL collocate_pgf_product_rspace(lmaxscp,alpha,0,0,0.0_dp,0,ra,& 00283 (/0.0_dp,0.0_dp,0.0_dp/),0.0_dp,1.0_dp,pab,0,0,rs_rho,& 00284 cell,cube_info,eps_rho_rspace,ga_gb_function=FUNC_AB,& 00285 ithread=ithread,use_subpatch=.TRUE.,subpatch_pattern=subpatch_pattern,error=error) 00286 00287 END DO 00288 00289 DEALLOCATE ( pab, cores, STAT=ierr ) 00290 CPPrecondition(ierr==0,cp_failure_level,routineP,error,failure) 00291 00292 END DO 00293 00294 CALL rs_pw_transfer(rs_rho,scp%pw,rs2pw,error=error) 00295 CALL rs_grid_release(rs_rho, error=error) 00296 00297 CALL timestop(handle) 00298 00299 END SUBROUTINE calculate_scp_charge 00300 ! ***************************************************************************** 00303 SUBROUTINE calculate_rho_nlcc(rho_nlcc,qs_env,error) 00304 00305 TYPE(pw_p_type), INTENT(INOUT) :: rho_nlcc 00306 TYPE(qs_environment_type), POINTER :: qs_env 00307 TYPE(cp_error_type), INTENT(INOUT) :: error 00308 00309 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_nlcc', 00310 routineP = moduleN//':'//routineN 00311 00312 INTEGER :: atom_a, handle, iatom, ierr, 00313 iexp_nlcc, ikind, ithread, j, 00314 n, natom, nc, nexp_nlcc, ni, 00315 npme, nthread 00316 INTEGER(KIND=int_8) :: subpatch_pattern 00317 INTEGER, DIMENSION(:), POINTER :: atom_list, cores, nct_nlcc 00318 LOGICAL :: failure, nlcc 00319 REAL(KIND=dp) :: alpha, eps_rho_rspace 00320 REAL(KIND=dp), DIMENSION(3) :: ra 00321 REAL(KIND=dp), DIMENSION(:), POINTER :: alpha_nlcc 00322 REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_nlcc, pab 00323 TYPE(atomic_kind_type), DIMENSION(:), 00324 POINTER :: atomic_kind_set 00325 TYPE(atomic_kind_type), POINTER :: atomic_kind 00326 TYPE(cell_type), POINTER :: cell 00327 TYPE(cp_para_env_type), POINTER :: para_env 00328 TYPE(cube_info_type) :: cube_info 00329 TYPE(dft_control_type), POINTER :: dft_control 00330 TYPE(gth_potential_type), POINTER :: gth_potential 00331 TYPE(particle_type), DIMENSION(:), 00332 POINTER :: particle_set 00333 TYPE(pw_env_type), POINTER :: pw_env 00334 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00335 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 00336 TYPE(realspace_grid_type), POINTER :: rs_rho 00337 00338 CALL timeset(routineN,handle) 00339 00340 failure = .FALSE. 00341 00342 NULLIFY(atomic_kind,cell,dft_control,pab,atomic_kind_set,particle_set,& 00343 atom_list,pw_env,rs_rho,auxbas_rs_desc,auxbas_pw_pool,cores) 00344 00345 CALL get_qs_env(qs_env=qs_env,& 00346 atomic_kind_set=atomic_kind_set,& 00347 cell=cell,& 00348 dft_control=dft_control,& 00349 particle_set=particle_set,& 00350 para_env=para_env,pw_env=pw_env,error=error) 00351 CALL pw_env_get(pw_env,auxbas_rs_desc=auxbas_rs_desc,auxbas_rs_grid=rs_rho,& 00352 auxbas_pw_pool=auxbas_pw_pool,error=error) 00353 cube_info=pw_env%cube_info(1) 00354 ! be careful in parallel nsmax is choosen with multigrid in mind! 00355 CALL rs_grid_retain(rs_rho, error=error) 00356 CALL rs_grid_zero(rs_rho) 00357 00358 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 00359 00360 DO ikind=1,SIZE(atomic_kind_set) 00361 00362 atomic_kind => atomic_kind_set(ikind) 00363 00364 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00365 natom=natom,& 00366 atom_list=atom_list,& 00367 gth_potential=gth_potential) 00368 00369 IF (.NOT.ASSOCIATED(gth_potential)) CYCLE 00370 CALL get_potential(potential=gth_potential,nlcc_present=nlcc,nexp_nlcc=nexp_nlcc,& 00371 alpha_nlcc=alpha_nlcc,nct_nlcc=nct_nlcc,cval_nlcc=cval_nlcc) 00372 00373 IF ( .NOT. nlcc ) CYCLE 00374 00375 DO iexp_nlcc=1,nexp_nlcc 00376 00377 alpha=alpha_nlcc(iexp_nlcc) 00378 nc=nct_nlcc(iexp_nlcc) 00379 00380 ni = ncoset(2*nc-2) 00381 ALLOCATE(pab(ni,1),STAT=ierr) 00382 CPPrecondition(ierr==0,cp_failure_level,routineP,error,failure) 00383 pab = 0._dp 00384 00385 nthread = 1 00386 ithread=0 00387 00388 CALL reallocate ( cores, 1, natom ) 00389 npme = 0 00390 cores = 0 00391 00392 ! prepare core function 00393 DO j=1,nc 00394 SELECT CASE (j) 00395 CASE (1) 00396 pab(1,1) = cval_nlcc(1,iexp_nlcc) 00397 CASE (2) 00398 n = coset(2,0,0) 00399 pab(n,1) = cval_nlcc(2,iexp_nlcc)/alpha**2 00400 n = coset(0,2,0) 00401 pab(n,1) = cval_nlcc(2,iexp_nlcc)/alpha**2 00402 n = coset(0,0,2) 00403 pab(n,1) = cval_nlcc(2,iexp_nlcc)/alpha**2 00404 CASE (3) 00405 n = coset(4,0,0) 00406 pab(n,1) = cval_nlcc(3,iexp_nlcc)/alpha**4 00407 n = coset(0,4,0) 00408 pab(n,1) = cval_nlcc(3,iexp_nlcc)/alpha**4 00409 n = coset(0,0,4) 00410 pab(n,1) = cval_nlcc(3,iexp_nlcc)/alpha**4 00411 n = coset(2,2,0) 00412 pab(n,1) = 2._dp*cval_nlcc(3,iexp_nlcc)/alpha**4 00413 n = coset(2,0,2) 00414 pab(n,1) = 2._dp*cval_nlcc(3,iexp_nlcc)/alpha**4 00415 n = coset(0,2,2) 00416 pab(n,1) = 2._dp*cval_nlcc(3,iexp_nlcc)/alpha**4 00417 CASE (4) 00418 n = coset(6,0,0) 00419 pab(n,1) = cval_nlcc(4,iexp_nlcc)/alpha**6 00420 n = coset(0,6,0) 00421 pab(n,1) = cval_nlcc(4,iexp_nlcc)/alpha**6 00422 n = coset(0,0,6) 00423 pab(n,1) = cval_nlcc(4,iexp_nlcc)/alpha**6 00424 n = coset(4,2,0) 00425 pab(n,1) = 3._dp*cval_nlcc(4,iexp_nlcc)/alpha**6 00426 n = coset(4,0,2) 00427 pab(n,1) = 3._dp*cval_nlcc(4,iexp_nlcc)/alpha**6 00428 n = coset(2,4,0) 00429 pab(n,1) = 3._dp*cval_nlcc(4,iexp_nlcc)/alpha**6 00430 n = coset(2,0,4) 00431 pab(n,1) = 3._dp*cval_nlcc(4,iexp_nlcc)/alpha**6 00432 n = coset(0,4,2) 00433 pab(n,1) = 3._dp*cval_nlcc(4,iexp_nlcc)/alpha**6 00434 n = coset(0,2,4) 00435 pab(n,1) = 3._dp*cval_nlcc(4,iexp_nlcc)/alpha**6 00436 n = coset(2,2,2) 00437 pab(n,1) = 6._dp*cval_nlcc(4,iexp_nlcc)/alpha**6 00438 CASE DEFAULT 00439 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00440 END SELECT 00441 END DO 00442 IF (dft_control%nspins==2) pab=pab*0.5_dp 00443 00444 DO iatom = 1, natom 00445 atom_a = atom_list(iatom) 00446 ra(:) = pbc(particle_set(atom_a)%r,cell) 00447 IF(rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN 00448 ! replicated realspace grid, split the atoms up between procs 00449 IF (MODULO(iatom,rs_rho%desc%group_size) == rs_rho % desc % my_pos ) THEN 00450 npme = npme + 1 00451 cores (npme) = iatom 00452 ENDIF 00453 ELSE 00454 npme = npme + 1 00455 cores (npme) = iatom 00456 ENDIF 00457 END DO 00458 00459 DO j=1,npme 00460 00461 iatom = cores(j) 00462 atom_a = atom_list(iatom) 00463 ra(:) = pbc(particle_set(atom_a)%r,cell) 00464 subpatch_pattern=0 00465 ni = 2*nc-2 00466 CALL collocate_pgf_product_rspace(ni,1/(2*alpha**2),0,0,0.0_dp,0,ra,& 00467 (/0.0_dp,0.0_dp,0.0_dp/),0.0_dp,1.0_dp,pab,0,0,rs_rho,& 00468 cell,cube_info,eps_rho_rspace,ga_gb_function=FUNC_AB,& 00469 ithread=ithread,use_subpatch=.TRUE.,subpatch_pattern=subpatch_pattern,error=error) 00470 00471 END DO 00472 00473 DEALLOCATE ( pab, STAT=ierr ) 00474 CPPrecondition(ierr==0,cp_failure_level,routineP,error,failure) 00475 00476 END DO 00477 00478 END DO 00479 00480 IF (ASSOCIATED(cores)) THEN 00481 DEALLOCATE (cores,STAT=ierr) 00482 CPPrecondition(ierr==0,cp_failure_level,routineP,error,failure) 00483 END IF 00484 00485 CALL rs_pw_transfer(rs_rho,rho_nlcc%pw,rs2pw,error=error) 00486 CALL rs_grid_release(rs_rho, error=error) 00487 00488 CALL timestop(handle) 00489 00490 END SUBROUTINE calculate_rho_nlcc 00491 00492 ! ***************************************************************************** 00495 SUBROUTINE calculate_ppl_grid(vppl,qs_env,error) 00496 00497 TYPE(pw_p_type), INTENT(INOUT) :: vppl 00498 TYPE(qs_environment_type), POINTER :: qs_env 00499 TYPE(cp_error_type), INTENT(INOUT) :: error 00500 00501 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ppl_grid', 00502 routineP = moduleN//':'//routineN 00503 00504 INTEGER :: atom_a, handle, iatom, ierr, 00505 ikind, ithread, j, lppl, n, 00506 natom, ni, npme, nthread 00507 INTEGER(KIND=int_8) :: subpatch_pattern 00508 INTEGER, DIMENSION(:), POINTER :: atom_list, cores 00509 LOGICAL :: failure 00510 REAL(KIND=dp) :: alpha, eps_rho_rspace 00511 REAL(KIND=dp), DIMENSION(3) :: ra 00512 REAL(KIND=dp), DIMENSION(:), POINTER :: cexp_ppl 00513 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab 00514 TYPE(atomic_kind_type), DIMENSION(:), 00515 POINTER :: atomic_kind_set 00516 TYPE(atomic_kind_type), POINTER :: atomic_kind 00517 TYPE(cell_type), POINTER :: cell 00518 TYPE(cp_para_env_type), POINTER :: para_env 00519 TYPE(cube_info_type) :: cube_info 00520 TYPE(dft_control_type), POINTER :: dft_control 00521 TYPE(gth_potential_type), POINTER :: gth_potential 00522 TYPE(particle_type), DIMENSION(:), 00523 POINTER :: particle_set 00524 TYPE(pw_env_type), POINTER :: pw_env 00525 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00526 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 00527 TYPE(realspace_grid_type), POINTER :: rs_rho 00528 00529 CALL timeset(routineN,handle) 00530 00531 failure = .FALSE. 00532 00533 NULLIFY(atomic_kind,cell,dft_control,pab,atomic_kind_set,particle_set,& 00534 atom_list,pw_env,rs_rho,auxbas_rs_desc,auxbas_pw_pool,cores) 00535 00536 CALL get_qs_env(qs_env=qs_env,& 00537 atomic_kind_set=atomic_kind_set,& 00538 cell=cell,& 00539 dft_control=dft_control,& 00540 particle_set=particle_set,& 00541 para_env=para_env,pw_env=pw_env,error=error) 00542 CALL pw_env_get(pw_env,auxbas_rs_desc=auxbas_rs_desc,auxbas_rs_grid=rs_rho,& 00543 auxbas_pw_pool=auxbas_pw_pool,error=error) 00544 cube_info=pw_env%cube_info(1) 00545 ! be careful in parallel nsmax is choosen with multigrid in mind! 00546 CALL rs_grid_retain(rs_rho, error=error) 00547 CALL rs_grid_zero(rs_rho) 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,alpha_ppl=alpha,nexp_ppl=lppl,cexp_ppl=cexp_ppl) 00562 00563 IF ( lppl <= 0 ) CYCLE 00564 00565 ni = ncoset(2*lppl-2) 00566 ALLOCATE(pab(ni,1),STAT=ierr) 00567 CPPrecondition(ierr==0,cp_failure_level,routineP,error,failure) 00568 pab = 0._dp 00569 00570 nthread = 1 00571 ithread=0 00572 00573 CALL reallocate ( cores, 1, natom ) 00574 npme = 0 00575 cores = 0 00576 00577 ! prepare core function 00578 DO j=1,lppl 00579 SELECT CASE (j) 00580 CASE (1) 00581 pab(1,1) = cexp_ppl(1) 00582 CASE (2) 00583 n = coset(2,0,0) 00584 pab(n,1) = cexp_ppl(2) 00585 n = coset(0,2,0) 00586 pab(n,1) = cexp_ppl(2) 00587 n = coset(0,0,2) 00588 pab(n,1) = cexp_ppl(2) 00589 CASE (3) 00590 n = coset(4,0,0) 00591 pab(n,1) = cexp_ppl(3) 00592 n = coset(0,4,0) 00593 pab(n,1) = cexp_ppl(3) 00594 n = coset(0,0,4) 00595 pab(n,1) = cexp_ppl(3) 00596 n = coset(2,2,0) 00597 pab(n,1) = 2._dp*cexp_ppl(3) 00598 n = coset(2,0,2) 00599 pab(n,1) = 2._dp*cexp_ppl(3) 00600 n = coset(0,2,2) 00601 pab(n,1) = 2._dp*cexp_ppl(3) 00602 CASE (4) 00603 n = coset(6,0,0) 00604 pab(n,1) = cexp_ppl(4) 00605 n = coset(0,6,0) 00606 pab(n,1) = cexp_ppl(4) 00607 n = coset(0,0,6) 00608 pab(n,1) = cexp_ppl(4) 00609 n = coset(4,2,0) 00610 pab(n,1) = 3._dp*cexp_ppl(4) 00611 n = coset(4,0,2) 00612 pab(n,1) = 3._dp*cexp_ppl(4) 00613 n = coset(2,4,0) 00614 pab(n,1) = 3._dp*cexp_ppl(4) 00615 n = coset(2,0,4) 00616 pab(n,1) = 3._dp*cexp_ppl(4) 00617 n = coset(0,4,2) 00618 pab(n,1) = 3._dp*cexp_ppl(4) 00619 n = coset(0,2,4) 00620 pab(n,1) = 3._dp*cexp_ppl(4) 00621 n = coset(2,2,2) 00622 pab(n,1) = 6._dp*cexp_ppl(4) 00623 CASE DEFAULT 00624 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00625 END SELECT 00626 END DO 00627 00628 DO iatom = 1, natom 00629 atom_a = atom_list(iatom) 00630 ra(:) = pbc(particle_set(atom_a)%r,cell) 00631 IF(rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN 00632 ! replicated realspace grid, split the atoms up between procs 00633 IF (MODULO(iatom,rs_rho%desc%group_size) == rs_rho % desc % my_pos ) THEN 00634 npme = npme + 1 00635 cores (npme) = iatom 00636 ENDIF 00637 ELSE 00638 npme = npme + 1 00639 cores (npme) = iatom 00640 ENDIF 00641 END DO 00642 00643 IF(npme.GT.0)THEN 00644 DO j=1,npme 00645 00646 iatom = cores(j) 00647 atom_a = atom_list(iatom) 00648 ra(:) = pbc(particle_set(atom_a)%r,cell) 00649 subpatch_pattern=0 00650 ni = 2*lppl-2 00651 CALL collocate_pgf_product_rspace(ni,alpha,0,0,0.0_dp,0,ra,& 00652 (/0.0_dp,0.0_dp,0.0_dp/),0.0_dp,1.0_dp,pab,0,0,rs_rho,& 00653 cell,cube_info,eps_rho_rspace,ga_gb_function=FUNC_AB,& 00654 ithread=ithread,use_subpatch=.TRUE.,subpatch_pattern=subpatch_pattern,error=error) 00655 00656 END DO 00657 ENDIF 00658 00659 DEALLOCATE ( pab, STAT=ierr ) 00660 CPPrecondition(ierr==0,cp_failure_level,routineP,error,failure) 00661 00662 END DO 00663 00664 IF (ASSOCIATED(cores)) THEN 00665 DEALLOCATE (cores,STAT=ierr) 00666 CPPrecondition(ierr==0,cp_failure_level,routineP,error,failure) 00667 END IF 00668 00669 CALL rs_pw_transfer(rs_rho,vppl%pw,rs2pw,error=error) 00670 CALL rs_grid_release(rs_rho, error=error) 00671 00672 CALL timestop(handle) 00673 00674 END SUBROUTINE calculate_ppl_grid 00675 00676 00677 ! ***************************************************************************** 00690 SUBROUTINE collocate_atomic_charge_density(total_rho, qs_env, error) 00691 00692 REAL(KIND=dp), INTENT(OUT) :: total_rho 00693 TYPE(qs_environment_type), POINTER :: qs_env 00694 TYPE(cp_error_type), INTENT(INOUT) :: error 00695 00696 CHARACTER(len=*), PARAMETER :: 00697 routineN = 'collocate_atomic_charge_density', 00698 routineP = moduleN//':'//routineN 00699 00700 INTEGER :: handle, i, iatom, ierr, igrid_level, ikind, ipgf, iset, 00701 ithread, maxco, na1, natom, ncoa, nkind, nseta, sgfa, unit_nr 00702 INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, 00703 npgfa, nsgfa 00704 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa 00705 LOGICAL :: failure, map_consistent 00706 REAL(KIND=dp) :: eps_rho_rspace 00707 REAL(KIND=dp), DIMENSION(3) :: ra 00708 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab, sphi_a, work, zeta 00709 TYPE(atomic_kind_type), DIMENSION(:), 00710 POINTER :: atomic_kind_set 00711 TYPE(atomic_kind_type), POINTER :: atomic_kind 00712 TYPE(cell_type), POINTER :: cell 00713 TYPE(cp_logger_type), POINTER :: logger 00714 TYPE(cube_info_type), DIMENSION(:), 00715 POINTER :: cube_info 00716 TYPE(dft_control_type), POINTER :: dft_control 00717 TYPE(gridlevel_info_type), POINTER :: gridlevel_info 00718 TYPE(gto_basis_set_type), POINTER :: aux_basis_set 00719 TYPE(particle_type), DIMENSION(:), 00720 POINTER :: particle_set 00721 TYPE(pw_env_type), POINTER :: pw_env 00722 TYPE(pw_p_type), DIMENSION(:), POINTER :: mgrid_gspace, mgrid_rspace 00723 TYPE(pw_pool_p_type), DIMENSION(:), 00724 POINTER :: pw_pools 00725 TYPE(qs_rho_type), POINTER :: rho_struct 00726 TYPE(realspace_grid_desc_p_type), 00727 DIMENSION(:), POINTER :: rs_descs 00728 TYPE(realspace_grid_p_type), 00729 DIMENSION(:), POINTER :: rs_rho 00730 00731 NULLIFY(aux_basis_set, atomic_kind_set, atomic_kind, npgfa, cell, particle_set, & 00732 sphi_a, rs_rho, pw_env, cube_info, dft_control, rho_struct, atom_list, & 00733 first_sgfa, pab, work) 00734 NULLIFY(gridlevel_info, rs_descs, pw_pools, mgrid_gspace, mgrid_rspace) 00735 00736 CALL timeset(routineN,handle) 00737 00738 logger => cp_error_get_logger(error) 00739 00740 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, cell=cell, & 00741 particle_set=particle_set, pw_env=pw_env, rho=rho_struct, & 00742 dft_control=dft_control, error=error) 00743 00744 cube_info => pw_env%cube_info 00745 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 00746 map_consistent = dft_control%qs_control%map_consistent 00747 ithread = 0 00748 gridlevel_info=>pw_env%gridlevel_info 00749 00750 ! *** set up the pw multi-grids *** ! 00751 CPPrecondition(ASSOCIATED(pw_env), cp_failure_level, routineP, error, failure) 00752 CALL pw_env_get(pw_env=pw_env, rs_descs=rs_descs, rs_grids=rs_rho, & 00753 pw_pools=pw_pools, error=error) 00754 00755 CALL pw_pools_create_pws(pw_pools, mgrid_rspace, & 00756 use_data=REALDATA3D, in_space=REALSPACE, & 00757 error=error) 00758 00759 CALL pw_pools_create_pws(pw_pools, mgrid_gspace, & 00760 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE, & 00761 error=error) 00762 00763 ! *** set up the rs multi-grids *** ! 00764 DO igrid_level = 1,gridlevel_info%ngrid_levels 00765 CALL rs_grid_retain(rs_rho(igrid_level)%rs_grid, error=error) 00766 CALL rs_grid_zero(rs_rho(igrid_level)%rs_grid) 00767 END DO 00768 00769 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, maxco=maxco) 00770 ALLOCATE (pab(maxco,1),STAT=ierr) 00771 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00772 "pab",maxco*dp_size) 00773 ALLOCATE (work(maxco,1),STAT=ierr) 00774 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00775 "work",maxco*dp_size) 00776 00777 nkind = SIZE(atomic_kind_set) 00778 00779 DO ikind = 1,nkind 00780 atomic_kind => atomic_kind_set(ikind) 00781 00782 CALL get_atomic_kind(atomic_kind=atomic_kind, aux_basis_set=aux_basis_set, natom=natom, & 00783 atom_list=atom_list) 00784 00785 CALL get_gto_basis_set(gto_basis_set=aux_basis_set, lmax=la_max, lmin=la_min, zet=zeta, & 00786 nset=nseta, npgf=npgfa, sphi=sphi_a, first_sgf=first_sgfa, nsgf_set=nsgfa) 00787 00788 DO iatom = 1,natom 00789 ! ra(:) = pbc(particle_set(iatom)%r, cell) 00790 ra(:) = pbc(particle_set(atom_list(iatom))%r, cell) 00791 00792 DO iset = 1,nseta 00793 00794 sgfa = first_sgfa(1,iset) 00795 ncoa = npgfa(iset)*ncoset(la_max(iset)) 00796 00797 DO i = 1,nsgfa(iset) 00798 work(i,1) = 1.0_dp 00799 END DO 00800 CALL dgemm("N","N",ncoa,1,nsgfa(iset),1.0_dp, sphi_a(1,sgfa),SIZE(sphi_a,1), & 00801 work(1,1),SIZE(work,1),0.0_dp,pab(1,1),SIZE(pab,1)) 00802 00803 DO ipgf = 1,npgfa(iset) 00804 00805 na1 = (ipgf-1)*ncoset(la_max(iset)) 00806 igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf,iset)) 00807 00808 CALL collocate_pgf_product_rspace(la_max=la_max(iset), zeta=zeta(ipgf,iset), & 00809 la_min=la_min(iset),& 00810 lb_max=0, zetb=0.0_dp, lb_min=0,& 00811 ra=ra,rab=(/0.0_dp,0.0_dp,0.0_dp/),rab2=0.0_dp,& 00812 scale=1.0_dp, pab=pab, o1=na1, o2=0,& 00813 rsgrid=rs_rho(igrid_level)%rs_grid,cell=cell,& 00814 cube_info=cube_info(igrid_level),& 00815 eps_rho_rspace=eps_rho_rspace,& 00816 ga_gb_function=FUNC_AB, ithread=ithread, & 00817 map_consistent=map_consistent,use_subpatch=.FALSE.,error=error) 00818 END DO 00819 END DO 00820 END DO 00821 END DO 00822 00823 DEALLOCATE (pab,STAT=ierr) 00824 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"pab") 00825 00826 DEALLOCATE (work,STAT=ierr) 00827 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"work") 00828 00829 IF (gridlevel_info%ngrid_levels==1) THEN 00830 CALL rs_pw_transfer(rs=rs_rho(1)%rs_grid, pw=qs_env%rho%rho_r(1)%pw, dir=rs2pw,error=error) 00831 CALL rs_grid_release(rs_rho(1)%rs_grid, error=error) 00832 CALL pw_transfer(qs_env%rho%rho_r(1)%pw, qs_env%rho%rho_g(1)%pw,error=error) 00833 IF (qs_env%rho%rho_r(1)%pw%pw_grid%spherical) THEN 00834 CALL pw_transfer(qs_env%rho%rho_g(1)%pw, qs_env%rho%rho_r(1)%pw,error=error) 00835 END IF 00836 ELSE 00837 DO igrid_level = 1,gridlevel_info%ngrid_levels 00838 CALL rs_pw_transfer(rs=rs_rho(igrid_level)%rs_grid, & 00839 pw=mgrid_rspace(igrid_level)%pw, dir=rs2pw, error=error) 00840 CALL rs_grid_release(rs_rho(igrid_level)%rs_grid, error=error) 00841 END DO 00842 00843 CALL pw_zero(qs_env%rho%rho_g(1)%pw,error=error) 00844 DO igrid_level=1, gridlevel_info%ngrid_levels 00845 CALL pw_transfer(mgrid_rspace(igrid_level)%pw, & 00846 mgrid_gspace(igrid_level)%pw,error=error) 00847 CALL pw_axpy(mgrid_gspace(igrid_level)%pw, qs_env%rho%rho_g(1)%pw,& 00848 error=error) 00849 END DO 00850 CALL pw_transfer(qs_env%rho%rho_g(1)%pw, qs_env%rho%rho_r(1)%pw,error=error) 00851 END IF 00852 00853 total_rho = pw_integrate_function(qs_env%rho%rho_r(1)%pw,isign=-1,& 00854 error=error) 00855 qs_env%rho%tot_rho_r = total_rho 00856 unit_nr=cp_logger_get_default_io_unit(logger) 00857 IF (unit_nr>0) THEN 00858 WRITE (unit_nr,*) "Total Rho =", total_rho 00859 END IF 00860 00861 ! *** give back the multi-grids *** ! 00862 CALL pw_pools_give_back_pws(pw_pools, mgrid_gspace, error=error) 00863 CALL pw_pools_give_back_pws(pw_pools, mgrid_rspace, error=error) 00864 00865 qs_env%rho%rho_r_valid = .TRUE. 00866 qs_env%rho%rho_g_valid = .TRUE. 00867 00868 CALL timestop(handle) 00869 00870 END SUBROUTINE collocate_atomic_charge_density 00871 00872 ! ***************************************************************************** 00875 SUBROUTINE calculate_rho_core(rho_core,total_rho,qs_env,only_nopaw,error) 00876 00877 TYPE(pw_p_type), INTENT(INOUT) :: rho_core 00878 REAL(KIND=dp), INTENT(OUT) :: total_rho 00879 TYPE(qs_environment_type), POINTER :: qs_env 00880 LOGICAL, INTENT(IN), OPTIONAL :: only_nopaw 00881 TYPE(cp_error_type), INTENT(INOUT) :: error 00882 00883 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_core', 00884 routineP = moduleN//':'//routineN 00885 00886 INTEGER :: atom_a, handle, iatom, ierr, 00887 ikind, ithread, j, natom, 00888 npme, nthread 00889 INTEGER(KIND=int_8) :: subpatch_pattern 00890 INTEGER, DIMENSION(:), POINTER :: atom_list, cores 00891 LOGICAL :: my_only_nopaw, paw_atom 00892 REAL(KIND=dp) :: alpha, eps_rho_rspace 00893 REAL(KIND=dp), DIMENSION(3) :: ra 00894 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab 00895 TYPE(atomic_kind_type), DIMENSION(:), 00896 POINTER :: atomic_kind_set 00897 TYPE(atomic_kind_type), POINTER :: atomic_kind 00898 TYPE(cell_type), POINTER :: cell 00899 TYPE(cp_para_env_type), POINTER :: para_env 00900 TYPE(cube_info_type) :: cube_info 00901 TYPE(dft_control_type), POINTER :: dft_control 00902 TYPE(particle_type), DIMENSION(:), 00903 POINTER :: particle_set 00904 TYPE(pw_env_type), POINTER :: pw_env 00905 TYPE(pw_p_type) :: rhoc_r 00906 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00907 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 00908 TYPE(realspace_grid_type), POINTER :: rs_rho 00909 00910 CALL timeset(routineN,handle) 00911 NULLIFY(atomic_kind,cell,dft_control,pab,atomic_kind_set,particle_set,& 00912 atom_list,pw_env,rs_rho,auxbas_rs_desc,auxbas_pw_pool,cores) 00913 ALLOCATE (pab(1,1),STAT=ierr) 00914 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00915 "pab",dp_size) 00916 00917 my_only_nopaw = .FALSE. 00918 IF(PRESENT(only_nopaw)) my_only_nopaw = only_nopaw 00919 00920 CALL get_qs_env(qs_env=qs_env,& 00921 atomic_kind_set=atomic_kind_set,& 00922 cell=cell,& 00923 dft_control=dft_control,& 00924 particle_set=particle_set,& 00925 para_env=para_env,pw_env=pw_env,error=error) 00926 CALL pw_env_get(pw_env,auxbas_rs_desc=auxbas_rs_desc,auxbas_rs_grid=rs_rho,& 00927 auxbas_pw_pool=auxbas_pw_pool,error=error) 00928 cube_info=pw_env%cube_info(1) 00929 ! be careful in parallel nsmax is choosen with multigrid in mind! 00930 CALL rs_grid_retain(rs_rho, error=error) 00931 CALL rs_grid_zero(rs_rho) 00932 00933 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 00934 00935 DO ikind=1,SIZE(atomic_kind_set) 00936 00937 atomic_kind => atomic_kind_set(ikind) 00938 00939 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00940 natom=natom,& 00941 paw_atom=paw_atom,& 00942 atom_list=atom_list,& 00943 alpha_core_charge=alpha,& 00944 ccore_charge=pab(1,1)) 00945 00946 IF(my_only_nopaw .AND. paw_atom ) CYCLE 00947 IF (alpha == 0.0_dp .OR. pab(1,1)== 0.0_dp) CYCLE 00948 00949 nthread = 1 00950 ithread=0 00951 00952 CALL reallocate ( cores, 1, natom ) 00953 npme = 0 00954 cores = 0 00955 00956 DO iatom = 1, natom 00957 atom_a = atom_list(iatom) 00958 ra(:) = pbc(particle_set(atom_a)%r,cell) 00959 IF(rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN 00960 ! replicated realspace grid, split the atoms up between procs 00961 IF (MODULO(iatom,rs_rho%desc%group_size) == rs_rho % desc % my_pos ) THEN 00962 npme = npme + 1 00963 cores (npme) = iatom 00964 ENDIF 00965 ELSE 00966 npme = npme + 1 00967 cores (npme) = iatom 00968 ENDIF 00969 END DO 00970 00971 IF(npme.GT.0)THEN 00972 DO j=1,npme 00973 00974 iatom = cores(j) 00975 atom_a = atom_list(iatom) 00976 ra(:) = pbc(particle_set(atom_a)%r,cell) 00977 subpatch_pattern=0 00978 CALL collocate_pgf_product_rspace(0,alpha,0,0,0.0_dp,0,ra,& 00979 (/0.0_dp,0.0_dp,0.0_dp/),0.0_dp,-1.0_dp,pab,0,0,rs_rho,& 00980 cell,cube_info,eps_rho_rspace,ga_gb_function=FUNC_AB,& 00981 ithread=ithread,use_subpatch=.TRUE.,subpatch_pattern=subpatch_pattern,error=error) 00982 00983 END DO 00984 ENDIF 00985 00986 END DO 00987 00988 IF (ASSOCIATED(cores)) THEN 00989 DEALLOCATE (cores,STAT=ierr) 00990 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"cores") 00991 END IF 00992 DEALLOCATE (pab,STAT=ierr) 00993 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"pab") 00994 00995 CALL pw_pool_create_pw(auxbas_pw_pool, rhoc_r%pw, & 00996 use_data=REALDATA3D,in_space=REALSPACE, error=error) 00997 00998 CALL rs_pw_transfer(rs_rho,rhoc_r%pw,rs2pw,error=error) 00999 CALL rs_grid_release(rs_rho, error=error) 01000 01001 total_rho = pw_integrate_function(rhoc_r%pw,isign=-1,error=error) 01002 01003 CALL pw_transfer(rhoc_r%pw,rho_core%pw,error=error) 01004 01005 CALL pw_pool_give_back_pw(auxbas_pw_pool, rhoc_r%pw, error=error) 01006 01007 CALL timestop(handle) 01008 01009 END SUBROUTINE calculate_rho_core 01010 ! ***************************************************************************** 01016 SUBROUTINE calculate_rho_single_gaussian(rho_gb,qs_env,iatom_in,error) 01017 01018 TYPE(pw_p_type), INTENT(INOUT) :: rho_gb 01019 TYPE(qs_environment_type), POINTER :: qs_env 01020 INTEGER, INTENT(IN) :: iatom_in 01021 TYPE(cp_error_type), INTENT(INOUT) :: error 01022 01023 CHARACTER(len=*), PARAMETER :: 01024 routineN = 'calculate_rho_single_gaussian', 01025 routineP = moduleN//':'//routineN 01026 01027 INTEGER :: atom_a, handle, iatom, ierr, 01028 npme 01029 INTEGER(KIND=int_8) :: subpatch_pattern 01030 REAL(KIND=dp) :: eps_rho_rspace 01031 REAL(KIND=dp), DIMENSION(3) :: ra 01032 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab 01033 TYPE(cell_type), POINTER :: cell 01034 TYPE(cp_para_env_type), POINTER :: para_env 01035 TYPE(dft_control_type), POINTER :: dft_control 01036 TYPE(pw_env_type), POINTER :: pw_env 01037 TYPE(pw_p_type) :: rhoc_r 01038 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 01039 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 01040 TYPE(realspace_grid_type), POINTER :: rs_rho 01041 01042 CALL timeset(routineN,handle) 01043 NULLIFY(cell,dft_control,pab,pw_env,rs_rho,auxbas_rs_desc,auxbas_pw_pool) 01044 ALLOCATE (pab(1,1),STAT=ierr) 01045 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01046 "pab",dp_size) 01047 01048 CALL get_qs_env(qs_env=qs_env,& 01049 cell=cell,& 01050 dft_control=dft_control,& 01051 para_env=para_env,pw_env=pw_env,error=error) 01052 CALL pw_env_get(pw_env,auxbas_rs_desc=auxbas_rs_desc,auxbas_rs_grid=rs_rho,& 01053 auxbas_pw_pool=auxbas_pw_pool,error=error) 01054 CALL rs_grid_retain(rs_rho, error=error) 01055 CALL rs_grid_zero(rs_rho) 01056 01057 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 01058 pab(1,1)=1.0_dp 01059 iatom=iatom_in 01060 01061 npme = 0 01062 01063 IF(rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN 01064 IF (MODULO(iatom,rs_rho%desc%group_size) == rs_rho % desc % my_pos ) THEN 01065 npme = npme + 1 01066 ENDIF 01067 ELSE 01068 npme = npme + 1 01069 ENDIF 01070 01071 IF(npme.GT.0)THEN 01072 atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom) 01073 ra(:) = pbc(qs_env%qmmm_env_qm%image_charge_pot%particles_all(atom_a)%r,cell) 01074 subpatch_pattern=0 01075 CALL collocate_pgf_product_rspace(0,qs_env%qmmm_env_qm%image_charge_pot%eta,& 01076 0,0,0.0_dp,0,ra,(/0.0_dp,0.0_dp,0.0_dp/),0.0_dp,1.0_dp,pab,0,0,rs_rho,& 01077 cell,pw_env%cube_info(1),eps_rho_rspace,ga_gb_function=FUNC_AB,& 01078 use_subpatch=.TRUE.,subpatch_pattern=subpatch_pattern,error=error) 01079 ENDIF 01080 01081 DEALLOCATE (pab,STAT=ierr) 01082 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"pab") 01083 01084 CALL pw_pool_create_pw(auxbas_pw_pool, rhoc_r%pw, & 01085 use_data=REALDATA3D,in_space=REALSPACE, error=error) 01086 01087 CALL rs_pw_transfer(rs_rho,rhoc_r%pw,rs2pw,error=error) 01088 CALL rs_grid_release(rs_rho, error=error) 01089 01090 CALL pw_transfer(rhoc_r%pw,rho_gb%pw,error=error) 01091 01092 CALL pw_pool_give_back_pw(auxbas_pw_pool, rhoc_r%pw, error=error) 01093 01094 CALL timestop(handle) 01095 01096 END SUBROUTINE calculate_rho_single_gaussian 01097 01098 ! ***************************************************************************** 01104 SUBROUTINE calculate_rho_metal(rho_metal,coeff,total_rho_metal,qs_env,error) 01105 01106 TYPE(pw_p_type), INTENT(INOUT) :: rho_metal 01107 REAL(KIND=dp), DIMENSION(:), POINTER :: coeff 01108 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: total_rho_metal 01109 TYPE(qs_environment_type), POINTER :: qs_env 01110 TYPE(cp_error_type), INTENT(INOUT) :: error 01111 01112 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_metal', 01113 routineP = moduleN//':'//routineN 01114 01115 INTEGER :: atom_a, handle, iatom, ierr, 01116 j, natom, npme 01117 INTEGER(KIND=int_8) :: subpatch_pattern 01118 INTEGER, DIMENSION(:), POINTER :: cores 01119 REAL(KIND=dp) :: eps_rho_rspace 01120 REAL(KIND=dp), DIMENSION(3) :: ra 01121 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab 01122 TYPE(cell_type), POINTER :: cell 01123 TYPE(cp_para_env_type), POINTER :: para_env 01124 TYPE(dft_control_type), POINTER :: dft_control 01125 TYPE(pw_env_type), POINTER :: pw_env 01126 TYPE(pw_p_type) :: rhoc_r 01127 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 01128 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 01129 TYPE(realspace_grid_type), POINTER :: rs_rho 01130 01131 CALL timeset(routineN,handle) 01132 01133 NULLIFY(cell,dft_control,pab,pw_env,rs_rho,auxbas_rs_desc,auxbas_pw_pool,& 01134 para_env,cores) 01135 01136 ALLOCATE (pab(1,1),STAT=ierr) 01137 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01138 "pab",dp_size) 01139 01140 CALL get_qs_env(qs_env=qs_env,& 01141 cell=cell,& 01142 dft_control=dft_control,& 01143 para_env=para_env,pw_env=pw_env,error=error) 01144 CALL pw_env_get(pw_env,auxbas_rs_desc=auxbas_rs_desc,auxbas_rs_grid=rs_rho,& 01145 auxbas_pw_pool=auxbas_pw_pool,error=error) 01146 CALL rs_grid_retain(rs_rho, error=error) 01147 CALL rs_grid_zero(rs_rho) 01148 01149 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 01150 pab(1,1)=1.0_dp 01151 01152 natom=SIZE(qs_env%qmmm_env_qm%image_charge_pot%image_mm_list) 01153 01154 CALL reallocate ( cores, 1, natom ) 01155 npme = 0 01156 cores = 0 01157 01158 DO iatom=1, natom 01159 IF(rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN 01160 IF (MODULO(iatom,rs_rho%desc%group_size) == rs_rho % desc % my_pos ) THEN 01161 npme = npme + 1 01162 cores (npme) = iatom 01163 ENDIF 01164 ELSE 01165 npme = npme + 1 01166 cores (npme) = iatom 01167 ENDIF 01168 ENDDO 01169 01170 IF(npme.GT.0)THEN 01171 DO j=1,npme 01172 iatom = cores(j) 01173 atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom) 01174 ra(:) = pbc(qs_env%qmmm_env_qm%image_charge_pot%particles_all(atom_a)%r,cell) 01175 subpatch_pattern=0 01176 CALL collocate_pgf_product_rspace(0,qs_env%qmmm_env_qm%image_charge_pot%eta,& 01177 0,0,0.0_dp,0,ra,(/0.0_dp,0.0_dp,0.0_dp/),0.0_dp,coeff(iatom),pab,0,0,rs_rho,& 01178 cell,pw_env%cube_info(1),eps_rho_rspace,ga_gb_function=FUNC_AB,& 01179 use_subpatch=.TRUE.,subpatch_pattern=subpatch_pattern,error=error) 01180 ENDDO 01181 ENDIF 01182 01183 DEALLOCATE (pab,cores,stat=ierr) 01184 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"pab") 01185 01186 CALL pw_pool_create_pw(auxbas_pw_pool, rhoc_r%pw, & 01187 use_data=REALDATA3D,in_space=REALSPACE, error=error) 01188 01189 CALL rs_pw_transfer(rs_rho,rhoc_r%pw,rs2pw,error=error) 01190 CALL rs_grid_release(rs_rho, error=error) 01191 01192 IF(PRESENT(total_rho_metal)) & 01193 !minus sign: account for the fact that rho_metal has opposite sign 01194 total_rho_metal = pw_integrate_function(rhoc_r%pw,isign=-1,error=error) 01195 01196 CALL pw_transfer(rhoc_r%pw,rho_metal%pw,error=error) 01197 CALL pw_pool_give_back_pw(auxbas_pw_pool, rhoc_r%pw, error=error) 01198 01199 CALL timestop(handle) 01200 01201 END SUBROUTINE calculate_rho_metal 01202 ! ***************************************************************************** 01208 SUBROUTINE calculate_rho_resp_single(rho_gb,qs_env,eta,iatom_in,error) 01209 01210 TYPE(pw_p_type), INTENT(INOUT) :: rho_gb 01211 TYPE(qs_environment_type), POINTER :: qs_env 01212 REAL(KIND=dp), INTENT(IN) :: eta 01213 INTEGER, INTENT(IN) :: iatom_in 01214 TYPE(cp_error_type), INTENT(INOUT) :: error 01215 01216 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_resp_single', 01217 routineP = moduleN//':'//routineN 01218 01219 INTEGER :: handle, iatom, ierr, npme 01220 INTEGER(KIND=int_8) :: subpatch_pattern 01221 REAL(KIND=dp) :: eps_rho_rspace 01222 REAL(KIND=dp), DIMENSION(3) :: ra 01223 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab 01224 TYPE(cell_type), POINTER :: cell 01225 TYPE(cp_para_env_type), POINTER :: para_env 01226 TYPE(dft_control_type), POINTER :: dft_control 01227 TYPE(particle_type), DIMENSION(:), 01228 POINTER :: particle_set 01229 TYPE(pw_env_type), POINTER :: pw_env 01230 TYPE(pw_p_type) :: rhoc_r 01231 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 01232 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 01233 TYPE(realspace_grid_type), POINTER :: rs_rho 01234 01235 CALL timeset(routineN,handle) 01236 NULLIFY(cell,dft_control,pab,pw_env,rs_rho,auxbas_rs_desc,auxbas_pw_pool,& 01237 particle_set) 01238 ALLOCATE (pab(1,1),STAT=ierr) 01239 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01240 "pab",dp_size) 01241 01242 CALL get_qs_env(qs_env=qs_env,& 01243 cell=cell,& 01244 dft_control=dft_control,& 01245 particle_set=particle_set,& 01246 para_env=para_env,pw_env=pw_env,error=error) 01247 CALL pw_env_get(pw_env,auxbas_rs_desc=auxbas_rs_desc,auxbas_rs_grid=rs_rho,& 01248 auxbas_pw_pool=auxbas_pw_pool,error=error) 01249 CALL rs_grid_retain(rs_rho, error=error) 01250 CALL rs_grid_zero(rs_rho) 01251 01252 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 01253 pab(1,1)=1.0_dp 01254 iatom=iatom_in 01255 01256 npme = 0 01257 01258 IF(rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN 01259 IF (MODULO(iatom,rs_rho%desc%group_size) == rs_rho % desc % my_pos ) THEN 01260 npme = npme + 1 01261 ENDIF 01262 ELSE 01263 npme = npme + 1 01264 ENDIF 01265 01266 IF(npme.GT.0)THEN 01267 ra(:) = pbc(particle_set(iatom)%r,cell) 01268 subpatch_pattern=0 01269 CALL collocate_pgf_product_rspace(0,eta,0,0,0.0_dp,0,ra,& 01270 (/0.0_dp,0.0_dp,0.0_dp/),0.0_dp,1.0_dp,pab,0,0,rs_rho,& 01271 cell,pw_env%cube_info(1),eps_rho_rspace,ga_gb_function=FUNC_AB,& 01272 use_subpatch=.TRUE.,subpatch_pattern=subpatch_pattern,error=error) 01273 ENDIF 01274 01275 DEALLOCATE (pab,STAT=ierr) 01276 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"pab") 01277 01278 CALL pw_pool_create_pw(auxbas_pw_pool, rhoc_r%pw, & 01279 use_data=REALDATA3D,in_space=REALSPACE, error=error) 01280 01281 CALL rs_pw_transfer(rs_rho,rhoc_r%pw,rs2pw,error=error) 01282 CALL rs_grid_release(rs_rho, error=error) 01283 01284 CALL pw_transfer(rhoc_r%pw,rho_gb%pw,error=error) 01285 01286 CALL pw_pool_give_back_pw(auxbas_pw_pool, rhoc_r%pw, error=error) 01287 01288 CALL timestop(handle) 01289 01290 END SUBROUTINE calculate_rho_resp_single 01291 ! ***************************************************************************** 01297 SUBROUTINE calculate_rho_resp_all(rho_resp,coeff,natom,eta,qs_env,error) 01298 01299 TYPE(pw_p_type), INTENT(INOUT) :: rho_resp 01300 REAL(KIND=dp), DIMENSION(:), POINTER :: coeff 01301 INTEGER, INTENT(IN) :: natom 01302 REAL(KIND=dp), INTENT(IN) :: eta 01303 TYPE(qs_environment_type), POINTER :: qs_env 01304 TYPE(cp_error_type), INTENT(INOUT) :: error 01305 01306 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_resp_all', 01307 routineP = moduleN//':'//routineN 01308 01309 INTEGER :: handle, iatom, ierr, j, npme 01310 INTEGER(KIND=int_8) :: subpatch_pattern 01311 INTEGER, DIMENSION(:), POINTER :: cores 01312 REAL(KIND=dp) :: eps_rho_rspace 01313 REAL(KIND=dp), DIMENSION(3) :: ra 01314 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab 01315 TYPE(cell_type), POINTER :: cell 01316 TYPE(cp_para_env_type), POINTER :: para_env 01317 TYPE(dft_control_type), POINTER :: dft_control 01318 TYPE(particle_type), DIMENSION(:), 01319 POINTER :: particle_set 01320 TYPE(pw_env_type), POINTER :: pw_env 01321 TYPE(pw_p_type) :: rhoc_r 01322 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 01323 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 01324 TYPE(realspace_grid_type), POINTER :: rs_rho 01325 01326 CALL timeset(routineN,handle) 01327 01328 NULLIFY(cell,dft_control,pab,pw_env,rs_rho,auxbas_rs_desc,auxbas_pw_pool,& 01329 para_env,particle_set) 01330 01331 ALLOCATE (pab(1,1),STAT=ierr) 01332 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01333 "pab",dp_size) 01334 01335 CALL get_qs_env(qs_env=qs_env,& 01336 cell=cell,& 01337 dft_control=dft_control,& 01338 particle_set=particle_set,& 01339 para_env=para_env,pw_env=pw_env,error=error) 01340 CALL pw_env_get(pw_env,auxbas_rs_desc=auxbas_rs_desc,auxbas_rs_grid=rs_rho,& 01341 auxbas_pw_pool=auxbas_pw_pool,error=error) 01342 CALL rs_grid_retain(rs_rho, error=error) 01343 CALL rs_grid_zero(rs_rho) 01344 01345 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 01346 pab(1,1)=1.0_dp 01347 01348 CALL reallocate ( cores, 1, natom ) 01349 npme = 0 01350 cores = 0 01351 01352 DO iatom=1, natom 01353 IF(rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN 01354 IF (MODULO(iatom,rs_rho%desc%group_size) == rs_rho % desc % my_pos ) THEN 01355 npme = npme + 1 01356 cores (npme) = iatom 01357 ENDIF 01358 ELSE 01359 npme = npme + 1 01360 cores (npme) = iatom 01361 ENDIF 01362 ENDDO 01363 01364 IF(npme.GT.0)THEN 01365 DO j=1,npme 01366 iatom = cores(j) 01367 ! atom_a = particle_set(iatom)%atom_index 01368 !atom_a = iatom 01369 !ra(:) = pbc(particle_set(atom_a)%r,cell) 01370 ra(:) = pbc(particle_set(iatom)%r,cell) 01371 subpatch_pattern=0 01372 CALL collocate_pgf_product_rspace(0,eta,& 01373 0,0,0.0_dp,0,ra,(/0.0_dp,0.0_dp,0.0_dp/),0.0_dp,coeff(iatom),pab,0,0,rs_rho,& 01374 cell,pw_env%cube_info(1),eps_rho_rspace,ga_gb_function=FUNC_AB,& 01375 use_subpatch=.TRUE.,subpatch_pattern=subpatch_pattern,error=error) 01376 ENDDO 01377 ENDIF 01378 01379 DEALLOCATE (pab,cores,stat=ierr) 01380 IF (ierr /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"pab") 01381 01382 CALL pw_pool_create_pw(auxbas_pw_pool, rhoc_r%pw, & 01383 use_data=REALDATA3D,in_space=REALSPACE, error=error) 01384 01385 CALL rs_pw_transfer(rs_rho,rhoc_r%pw,rs2pw,error=error) 01386 CALL rs_grid_release(rs_rho, error=error) 01387 01388 CALL pw_transfer(rhoc_r%pw,rho_resp%pw,error=error) 01389 CALL pw_pool_give_back_pw(auxbas_pw_pool, rhoc_r%pw, error=error) 01390 01391 CALL timestop(handle) 01392 01393 END SUBROUTINE calculate_rho_resp_all 01394 01395 ! ***************************************************************************** 01404 SUBROUTINE calculate_rho_elec(matrix_p, rho,rho_gspace, total_rho,& 01405 qs_env, soft_valid, compute_tau, compute_grad, basis_set_id, der_type, idir, task_list_external, error) 01406 01407 TYPE(cp_dbcsr_type), POINTER :: matrix_p 01408 TYPE(pw_p_type), INTENT(INOUT) :: rho, rho_gspace 01409 REAL(KIND=dp), INTENT(OUT) :: total_rho 01410 TYPE(qs_environment_type), POINTER :: qs_env 01411 LOGICAL, INTENT(IN), OPTIONAL :: soft_valid, compute_tau, 01412 compute_grad 01413 INTEGER, INTENT(IN), OPTIONAL :: basis_set_id, der_type, idir 01414 TYPE(task_list_type), OPTIONAL, POINTER :: task_list_external 01415 TYPE(cp_error_type), INTENT(inout) :: error 01416 01417 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_elec', 01418 routineP = moduleN//':'//routineN 01419 INTEGER, PARAMETER :: add_tasks = 1000, 01420 max_tasks = 2000 01421 REAL(kind=dp), PARAMETER :: mult_tasks = 2.0_dp 01422 01423 INTEGER :: bcol, brow, ga_gb_function, handle, i, iatom, iatom_old, 01424 igrid_level, igrid_level_dummy, ikind, ikind_old, ipair, ipgf, iset, 01425 iset_old, itask, ithread, jatom, jatom_old, jkind, jkind_old, jpgf, 01426 jset, jset_old, lb, lb_l, lbr, lbw, maxco, maxpgf, maxset, maxsgf, 01427 maxsgf_set, my_basis_set_id, my_der_type, my_idir, n, na1, na2, natoms, 01428 nb1, nb2, nblock, ncoa, ncob, nkind, nn, nr, nrlevel, nseta, nsetb, 01429 ntasks, nthread, nxy, nz, nzsize, segment, sgfa, sgfb, stat, ub 01430 INTEGER(kind=int_8), DIMENSION(:), 01431 POINTER :: atom_pair_recv, atom_pair_send 01432 INTEGER(kind=int_8), DIMENSION(:, :), 01433 POINTER :: tasks 01434 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, 01435 lb_min, npgfa, npgfb, nsgfa, 01436 nsgfb 01437 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 01438 LOGICAL :: atom_pair_changed, distributed_rs_grids, failure, found, 01439 map_consistent, my_compute_grad, my_compute_tau, my_soft, use_subpatch 01440 REAL(KIND=dp) :: eps_rho_rspace, rab2, scale, 01441 zetp 01442 REAL(KIND=dp), DIMENSION(3) :: ra, rab, rab_inv, rb 01443 REAL(KIND=dp), DIMENSION(:, :), POINTER :: dist_ab, p_block, pab, 01444 sphi_a, sphi_b, work, zeta, 01445 zetb 01446 REAL(KIND=dp), DIMENSION(:, :, :), 01447 POINTER :: pabt, workt 01448 TYPE(atomic_kind_type), DIMENSION(:), 01449 POINTER :: atomic_kind_set 01450 TYPE(atomic_kind_type), POINTER :: atomic_kind 01451 TYPE(cell_type), POINTER :: cell 01452 TYPE(cp_dbcsr_type), POINTER :: deltap 01453 TYPE(cp_para_env_type), POINTER :: para_env 01454 TYPE(cube_info_type), DIMENSION(:), 01455 POINTER :: cube_info 01456 TYPE(dft_control_type), POINTER :: dft_control 01457 TYPE(gridlevel_info_type), POINTER :: gridlevel_info 01458 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 01459 TYPE(lgrid_type), POINTER :: lgrid 01460 TYPE(particle_type), DIMENSION(:), 01461 POINTER :: particle_set 01462 TYPE(pw_env_type), POINTER :: pw_env 01463 TYPE(realspace_grid_desc_p_type), 01464 DIMENSION(:), POINTER :: rs_descs 01465 TYPE(realspace_grid_p_type), 01466 DIMENSION(:), POINTER :: rs_rho 01467 TYPE(section_vals_type), POINTER :: input, interp_section 01468 TYPE(task_list_type), POINTER :: task_list, task_list_soft 01469 01470 !$ INTEGER :: omp_get_thread_num, omp_get_max_threads 01471 01472 failure=.FALSE. 01473 NULLIFY(atomic_kind, cell, dft_control, orb_basis_set, deltap, & 01474 atomic_kind_set, particle_set, rs_rho, pw_env, rs_descs, & 01475 para_env, dist_ab, la_max, la_min, & 01476 lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, p_block, & 01477 sphi_a, sphi_b, zeta, zetb, first_sgfa, first_sgfb, tasks, pabt, & 01478 workt,lgrid) 01479 01480 debug_count=debug_count+1 01481 01482 ! by default, the full density is calculated 01483 my_soft=.FALSE. 01484 IF (PRESENT(soft_valid)) my_soft = soft_valid 01485 01486 ! by default, do not compute the kinetic energy density (tau) 01487 ! if compute_tau, all grids referening to rho are actually tau 01488 IF (PRESENT(compute_tau)) THEN 01489 my_compute_tau = compute_tau 01490 ELSE 01491 my_compute_tau = .FALSE. 01492 ENDIF 01493 01494 IF (PRESENT(compute_grad)) THEN 01495 my_compute_grad = compute_grad 01496 ELSE 01497 my_compute_grad = .FALSE. 01498 ENDIF 01499 01500 01501 CALL timeset(routineN,handle) 01502 my_der_type = 0 01503 my_idir = 0 01504 IF(PRESENT(der_type)) THEN 01505 my_der_type = der_type 01506 ga_gb_function = FUNC_DER(my_der_type) 01507 ELSE IF (my_compute_tau) THEN 01508 ga_gb_function = FUNC_DADB 01509 ELSE IF (my_compute_grad) THEN 01510 ga_gb_function = FUNC_DABpADB 01511 ELSE 01512 ga_gb_function = FUNC_AB 01513 ENDIF 01514 IF(PRESENT(idir)) my_idir = idir 01515 01516 IF( PRESENT(basis_set_id) ) THEN 01517 my_basis_set_id = basis_set_id 01518 ELSE 01519 my_basis_set_id = use_orb_basis_set 01520 END IF 01521 01522 01523 SELECT CASE (my_basis_set_id) 01524 CASE (use_orb_basis_set) 01525 CALL get_qs_env(qs_env=qs_env,& 01526 atomic_kind_set=atomic_kind_set,& 01527 cell=cell,& 01528 dft_control=dft_control,& 01529 particle_set=particle_set,& 01530 para_env=para_env,& 01531 task_list=task_list,& 01532 task_list_soft=task_list_soft,& 01533 input=input,& 01534 pw_env=pw_env,error=error) 01535 CASE (use_aux_fit_basis_set) 01536 CALL get_qs_env(qs_env=qs_env,& 01537 atomic_kind_set=atomic_kind_set,& 01538 cell=cell,& 01539 dft_control=dft_control,& 01540 particle_set=particle_set,& 01541 para_env=para_env,& 01542 task_list_aux_fit=task_list,& 01543 task_list_soft=task_list_soft,& 01544 input=input,& 01545 pw_env=pw_env,error=error) 01546 END SELECT 01547 01548 IF (PRESENT(task_list_external)) task_list=>task_list_external 01549 01550 ! *** assign from pw_env 01551 gridlevel_info=>pw_env%gridlevel_info 01552 cube_info=>pw_env%cube_info 01553 01554 ! *** Allocate work storage *** 01555 nthread = 1 01556 !$ nthread = omp_get_max_threads() 01557 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 01558 maxco=maxco,& 01559 maxsgf=maxsgf,& 01560 maxsgf_set=maxsgf_set,& 01561 basis_set_id=my_basis_set_id) 01562 nkind = SIZE(atomic_kind_set) 01563 CALL reallocate(pabt,1,maxco,1,maxco,0,nthread-1) 01564 CALL reallocate(workt,1,maxco,1,maxsgf_set,0,nthread-1) 01565 01566 ! find maximum numbers 01567 nkind = SIZE(atomic_kind_set) 01568 natoms = SIZE( particle_set ) 01569 maxset=0 01570 maxpgf=0 01571 DO ikind=1,nkind 01572 atomic_kind => atomic_kind_set(ikind) 01573 SELECT CASE( my_basis_set_id) 01574 CASE (use_orb_basis_set) 01575 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01576 softb = my_soft, & 01577 orb_basis_set=orb_basis_set) 01578 CASE (use_aux_fit_basis_set) 01579 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01580 softb = my_soft, & 01581 aux_fit_basis_set=orb_basis_set,& 01582 basis_set_id=my_basis_set_id) 01583 END SELECT 01584 01585 IF (.NOT.ASSOCIATED(orb_basis_set)) CYCLE 01586 01587 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 01588 npgf=npgfa, nset=nseta ) 01589 01590 maxset=MAX(nseta,maxset) 01591 maxpgf=MAX(MAXVAL(npgfa),maxpgf) 01592 END DO 01593 01594 ! get the task lists 01595 IF (my_soft) task_list=>task_list_soft 01596 CPPrecondition(ASSOCIATED(task_list),cp_failure_level,routineP,error,failure) 01597 tasks =>task_list%tasks 01598 dist_ab=>task_list%dist_ab 01599 atom_pair_send=>task_list%atom_pair_send 01600 atom_pair_recv=>task_list%atom_pair_recv 01601 ntasks=task_list%ntasks 01602 01603 ! *** set up the pw multi-grids 01604 CPPrecondition(ASSOCIATED(pw_env),cp_failure_level,routineP,error,failure) 01605 CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho, lgrid=lgrid, error=error) 01606 ! *** set up the rs multi-grids 01607 CPPrecondition(ASSOCIATED(rs_rho),cp_failure_level,routineP,error,failure) 01608 CPPrecondition(ASSOCIATED(rs_descs),cp_failure_level,routineP,error,failure) 01609 CPPrecondition(ASSOCIATED(lgrid),cp_failure_level,routineP,error,failure) 01610 distributed_rs_grids=.FALSE. 01611 01612 DO igrid_level=1,gridlevel_info%ngrid_levels 01613 CALL rs_grid_retain(rs_rho(igrid_level)%rs_grid, error=error) 01614 CALL rs_grid_zero(rs_rho(igrid_level)%rs_grid) 01615 IF ( rs_rho(igrid_level)%rs_grid%desc%distributed ) THEN 01616 distributed_rs_grids=.TRUE. 01617 ENDIF 01618 END DO 01619 01620 IF ( nthread > 1 ) THEN 01621 IF (.NOT. ASSOCIATED(lgrid%r)) THEN 01622 CALL lgrid_allocate_grid(lgrid, nthread, error) 01623 ENDIF 01624 END IF 01625 01626 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 01627 map_consistent = dft_control%qs_control%map_consistent 01628 ! *** Initialize working density matrix *** 01629 ! distributed rs grids require a matrix that will be changed 01630 ! whereas this is not the case for replicated grids 01631 IF (distributed_rs_grids) THEN 01632 NULLIFY(deltap) 01633 !CALL replicate_matrix(matrix_p,deltap,target_name="DeltaP",error=error) 01634 ALLOCATE(deltap) 01635 CALL cp_dbcsr_init(deltap, error=error) 01636 CALL cp_dbcsr_copy(deltap,matrix_p,name="DeltaP",error=error) 01637 ! this matrix has no strict sparsity pattern in parallel 01638 !deltap%sparsity_id=-1 01639 ELSE 01640 deltap=>matrix_p 01641 ENDIF 01642 01643 ! distribute the matrix 01644 IF (distributed_rs_grids) THEN 01645 CALL rs_distribute_matrix (rs_descs, deltap, atom_pair_send, atom_pair_recv, natoms, scatter=.TRUE., error=error) 01646 ENDIF 01647 01648 ! map all tasks on the grids 01649 01650 !$omp parallel default(none), & 01651 !$omp shared(ntasks,tasks,natoms,maxset,maxpgf,particle_set,pabt,workt), & 01652 !$omp shared(my_basis_set_id,my_soft,deltap,maxco,dist_ab,ncoset,nthread), & 01653 !$omp shared(cell,cube_info,eps_rho_rspace,ga_gb_function, my_idir,map_consistent), & 01654 !$omp shared(rs_rho,lgrid,gridlevel_info,task_list), & 01655 !$omp private(igrid_level,iatom,jatom,iset,jset,ipgf,jpgf,ikind,jkind,pab,work), & 01656 !$omp private(iatom_old,jatom_old,iset_old,jset_old,ikind_old,jkind_old), & 01657 !$omp private(brow,bcol,atomic_kind,orb_basis_set,first_sgfa,la_max,la_min), & 01658 !$omp private(npgfa,nseta,nsgfa,sphi_a,zeta,first_sgfb,lb_max,lb_min,npgfb), & 01659 !$omp private(nsetb,nsgfb,sphi_b,zetb,p_block,found), & 01660 !$omp private(atom_pair_changed,ncoa,sgfa,ncob,sgfb,rab,rab2,ra,rb,zetp), & 01661 !$omp private(na1,na2,nb1,nb2,scale,use_subpatch,rab_inv,ithread,error,lb,ub,n), & 01662 !$omp private(lb_l,nn,segment,itask,i,nz,nxy,nzsize,nrlevel,nblock,lbw,lbr,nr,igrid_level_dummy) 01663 01664 ithread = 0 01665 !$ ithread = omp_get_thread_num() 01666 pab => pabt(:,:,ithread) 01667 work => workt(:,:,ithread) 01668 01669 iatom_old = -1 ; jatom_old = -1 ; iset_old = -1 ; jset_old = -1 01670 ikind_old = -1 ; jkind_old = -1 01671 01672 ! Loop over each gridlevel first, then loop and load balance over atom pairs 01673 ! We only need a single lgrid, which we sum back onto the rsgrid after each 01674 ! grid level is completed 01675 loop_gridlevels: DO igrid_level=1,gridlevel_info%ngrid_levels 01676 01677 ! Only zero the region of the lgrid required for this grid level 01678 IF (nthread > 1) THEN 01679 lb = ithread*lgrid%ldim + 1 01680 ub = ithread*lgrid%ldim + rs_rho(igrid_level)%rs_grid%ngpts_local 01681 lgrid%r(lb:ub) = 0._dp 01682 END IF 01683 !$omp do schedule(dynamic, MAX(1,task_list%npairs(igrid_level)/(nthread*50))) 01684 loop_pairs: DO ipair = 1, task_list%npairs(igrid_level) 01685 loop_tasks: DO itask = task_list%taskstart(ipair,igrid_level), task_list%taskstop(ipair,igrid_level) 01686 !decode the atom pair and basis info (igrid_level_dummy equals do loop variable by construction). 01687 CALL int2pair(tasks(3,itask),igrid_level_dummy,iatom,jatom,iset,jset,ipgf,jpgf,natoms,maxset,maxpgf) 01688 ikind = particle_set(iatom)%atomic_kind%kind_number 01689 jkind = particle_set(jatom)%atomic_kind%kind_number 01690 01691 IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old) THEN 01692 01693 IF (iatom .NE. iatom_old) ra(:) = pbc(particle_set(iatom)%r,cell) 01694 01695 IF (iatom <= jatom) THEN 01696 brow = iatom 01697 bcol = jatom 01698 ELSE 01699 brow = jatom 01700 bcol = iatom 01701 END IF 01702 01703 IF(ikind .NE. ikind_old ) THEN 01704 SELECT CASE (my_basis_set_id) 01705 CASE (use_orb_basis_set) 01706 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& 01707 softb = my_soft, & 01708 orb_basis_set=orb_basis_set) 01709 CASE (use_aux_fit_basis_set) 01710 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& 01711 softb = my_soft, & 01712 aux_fit_basis_set=orb_basis_set,& 01713 basis_set_id=my_basis_set_id) 01714 END SELECT 01715 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 01716 first_sgf=first_sgfa,& 01717 lmax=la_max,& 01718 lmin=la_min,& 01719 npgf=npgfa,& 01720 nset=nseta,& 01721 nsgf_set=nsgfa,& 01722 sphi=sphi_a,& 01723 zet=zeta) 01724 ENDIF 01725 01726 IF (jkind .NE. jkind_old ) THEN 01727 SELECT CASE (my_basis_set_id) 01728 CASE (use_orb_basis_set) 01729 CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind,& 01730 softb = my_soft, & 01731 orb_basis_set=orb_basis_set) 01732 CASE (use_aux_fit_basis_set) 01733 CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind,& 01734 softb = my_soft, & 01735 aux_fit_basis_set=orb_basis_set,& 01736 basis_set_id=my_basis_set_id) 01737 END SELECT 01738 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 01739 first_sgf=first_sgfb,& 01740 lmax=lb_max,& 01741 lmin=lb_min,& 01742 npgf=npgfb,& 01743 nset=nsetb,& 01744 nsgf_set=nsgfb,& 01745 sphi=sphi_b,& 01746 zet=zetb) 01747 ENDIF 01748 01749 CALL cp_dbcsr_get_block_p(matrix=deltap,& 01750 row=brow,col=bcol,BLOCK=p_block,found=found) 01751 IF (.NOT.ASSOCIATED(p_block)) & 01752 CALL stop_program(routineN,moduleN,__LINE__,& 01753 "p_block not associated in deltap") 01754 01755 iatom_old = iatom 01756 jatom_old = jatom 01757 ikind_old = ikind 01758 jkind_old = jkind 01759 atom_pair_changed = .TRUE. 01760 01761 ELSE 01762 01763 atom_pair_changed = .FALSE. 01764 01765 ENDIF 01766 01767 IF (atom_pair_changed .OR. iset_old .NE. iset .OR. jset_old .NE. jset) THEN 01768 ncoa = npgfa(iset)*ncoset(la_max(iset)) 01769 sgfa = first_sgfa(1,iset) 01770 ncob = npgfb(jset)*ncoset(lb_max(jset)) 01771 sgfb = first_sgfb(1,jset) 01772 01773 IF (iatom <= jatom) THEN 01774 CALL dgemm("N","N",ncoa,nsgfb(jset),nsgfa(iset),& 01775 1.0_dp,sphi_a(1,sgfa),SIZE(sphi_a,1),& 01776 p_block(sgfa,sgfb),SIZE(p_block,1),& 01777 0.0_dp,work(1,1),maxco) 01778 CALL dgemm("N","T",ncoa,ncob,nsgfb(jset),& 01779 1.0_dp,work(1,1),maxco,& 01780 sphi_b(1,sgfb),SIZE(sphi_b,1),& 01781 0.0_dp,pab(1,1),maxco) 01782 ELSE 01783 CALL dgemm("N","N",ncob,nsgfa(iset),nsgfb(jset),& 01784 1.0_dp,sphi_b(1,sgfb),SIZE(sphi_b,1),& 01785 p_block(sgfb,sgfa),SIZE(p_block,1),& 01786 0.0_dp,work(1,1),maxco) 01787 CALL dgemm("N","T",ncob,ncoa,nsgfa(iset),& 01788 1.0_dp,work(1,1),maxco,& 01789 sphi_a(1,sgfa),SIZE(sphi_a,1),& 01790 0.0_dp,pab(1,1),maxco) 01791 END IF 01792 01793 iset_old = iset 01794 jset_old = jset 01795 ENDIF 01796 01797 rab(:) = dist_ab (:,itask) 01798 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 01799 rb(:) = ra(:) + rab(:) 01800 zetp = zeta(ipgf,iset) + zetb(jpgf,jset) 01801 01802 na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1 01803 na2 = ipgf*ncoset(la_max(iset)) 01804 nb1 = (jpgf - 1)*ncoset(lb_max(jset)) + 1 01805 nb2 = jpgf*ncoset(lb_max(jset)) 01806 01807 ! takes the density matrix symmetry in account, i.e. off-diagonal blocks need to be mapped 'twice' 01808 IF (iatom == jatom) THEN 01809 scale = 1.0_dp 01810 ELSE 01811 scale = 2.0_dp 01812 END IF 01813 01814 ! check whether we need to use fawzi's generalised collocation scheme 01815 IF(rs_rho(igrid_level)%rs_grid%desc%distributed)THEN 01816 !tasks(4,:) is 0 for replicated, 1 for distributed 2 for exceptional distributed tasks 01817 IF (tasks(4,itask) .EQ. 2 ) THEN 01818 use_subpatch = .TRUE. 01819 ELSE 01820 use_subpatch = .FALSE. 01821 ENDIF 01822 ELSE 01823 use_subpatch = .FALSE. 01824 ENDIF 01825 01826 IF (nthread > 1) THEN 01827 IF (iatom <= jatom) THEN 01828 CALL collocate_pgf_product_rspace(& 01829 la_max(iset),zeta(ipgf,iset),la_min(iset),& 01830 lb_max(jset),zetb(jpgf,jset),lb_min(jset),& 01831 ra,rab,rab2,scale,pab,na1-1,nb1-1,& 01832 rs_rho(igrid_level)%rs_grid,cell,cube_info(igrid_level),& 01833 eps_rho_rspace,& 01834 ga_gb_function=ga_gb_function, & 01835 idir=my_idir,& 01836 lgrid=lgrid,ithread=ithread, & 01837 map_consistent=map_consistent,use_subpatch=use_subpatch,& 01838 subpatch_pattern=tasks(6,itask),error=error) 01839 ELSE 01840 rab_inv=-rab 01841 CALL collocate_pgf_product_rspace(& 01842 lb_max(jset),zetb(jpgf,jset),lb_min(jset),& 01843 la_max(iset),zeta(ipgf,iset),la_min(iset),& 01844 rb,rab_inv,rab2,scale,pab,nb1-1,na1-1,& 01845 rs_rho(igrid_level)%rs_grid,cell,cube_info(igrid_level),& 01846 eps_rho_rspace,& 01847 ga_gb_function=ga_gb_function, & 01848 idir=my_idir,& 01849 lgrid=lgrid,ithread=ithread, & 01850 map_consistent=map_consistent,use_subpatch=use_subpatch,& 01851 subpatch_pattern=tasks(6,itask),error=error) 01852 END IF 01853 ELSE 01854 IF (iatom <= jatom) THEN 01855 CALL collocate_pgf_product_rspace(& 01856 la_max(iset),zeta(ipgf,iset),la_min(iset),& 01857 lb_max(jset),zetb(jpgf,jset),lb_min(jset),& 01858 ra,rab,rab2,scale,pab,na1-1,nb1-1,& 01859 rs_rho(igrid_level)%rs_grid,cell,cube_info(igrid_level),& 01860 eps_rho_rspace,& 01861 ga_gb_function=ga_gb_function, & 01862 idir=my_idir,& 01863 map_consistent=map_consistent,use_subpatch=use_subpatch,& 01864 subpatch_pattern=tasks(6,itask),error=error) 01865 ELSE 01866 rab_inv=-rab 01867 CALL collocate_pgf_product_rspace(& 01868 lb_max(jset),zetb(jpgf,jset),lb_min(jset),& 01869 la_max(iset),zeta(ipgf,iset),la_min(iset),& 01870 rb,rab_inv,rab2,scale,pab,nb1-1,na1-1,& 01871 rs_rho(igrid_level)%rs_grid,cell,cube_info(igrid_level),& 01872 eps_rho_rspace,& 01873 ga_gb_function=ga_gb_function, & 01874 idir=my_idir,& 01875 map_consistent=map_consistent,use_subpatch=use_subpatch,& 01876 subpatch_pattern=tasks(6,itask),error=error) 01877 END IF 01878 END IF 01879 END DO loop_tasks 01880 END DO loop_pairs 01881 !$omp end do 01882 01883 ! Now sum the thread-local grids back into the rs_grid (in parallel, each thread writes to a section of the rs_grid at a time) 01884 IF (nthread > 1) THEN 01885 nz = (1 + rs_rho(igrid_level)%rs_grid%ub_local(3) & 01886 - rs_rho(igrid_level)%rs_grid%lb_local(3)) 01887 nxy = (1 + rs_rho(igrid_level)%rs_grid%ub_local(1) & 01888 - rs_rho(igrid_level)%rs_grid%lb_local(1)) & 01889 *(1 + rs_rho(igrid_level)%rs_grid%ub_local(2) & 01890 - rs_rho(igrid_level)%rs_grid%lb_local(2)) 01891 01892 ! Work out the number of tree levels, and start the reduction 01893 01894 nrlevel = CEILING(LOG10(1.0_dp*nthread)/LOG10(2.0_dp)) 01895 01896 DO nr = 1, nrlevel 01897 nblock = 2**nr 01898 nzsize = MIN(nblock, nthread - (ithread/nblock)*nblock) 01899 itask = MODULO(ithread, nblock) 01900 01901 lb = (nz*itask)/nzsize 01902 ub = (nz*(itask+1))/nzsize - 1 01903 n = (ithread/nblock)*nblock 01904 lbw = 1 + n*lgrid%ldim + nxy*lb 01905 01906 n = n + nblock/2 01907 IF (n < nthread) THEN 01908 lbr = 1 + n*lgrid%ldim + nxy*lb 01909 01910 CALL daxpy(nxy*(1+ub-lb), 1.0_dp, lgrid%r(lbr), 1, lgrid%r(lbw), 1) 01911 END IF 01912 !$omp barrier 01913 END DO 01914 01915 ! Copy final result from first local grid to rs_grid 01916 lb = (nz*ithread)/nthread 01917 ub = (nz*(ithread+1))/nthread - 1 01918 nzsize = 1 + (ub - lb) 01919 01920 lb = lb + rs_rho(igrid_level)%rs_grid%lb_local(3) 01921 ub = ub + rs_rho(igrid_level)%rs_grid%lb_local(3) 01922 lbr = 1 + nxy*(nz*ithread/nthread) 01923 01924 CALL daxpy(nxy*nzsize, 1.0_dp, lgrid%r(lbr),1, & 01925 rs_rho(igrid_level)%rs_grid%r(:,:,lb:ub), 1) 01926 !$omp barrier 01927 END IF 01928 END DO loop_gridlevels 01929 !$omp end parallel 01930 01931 ! *** Release work storage *** 01932 01933 IF (distributed_rs_grids) THEN 01934 CALL cp_dbcsr_deallocate_matrix ( deltap ,error=error) 01935 ENDIF 01936 01937 DEALLOCATE (pabt,STAT=stat) 01938 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"pabt") 01939 DEALLOCATE (workt,STAT=stat) 01940 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"workt") 01941 01942 interp_section => section_vals_get_subs_vals(input,"DFT%MGRID%INTERPOLATOR",error=error) 01943 CALL density_rs2pw(pw_env,rs_rho,rho,rho_gspace,interp_section=interp_section,error=error) 01944 01945 total_rho = pw_integrate_function(rho%pw,isign=-1,error=error) 01946 CALL timestop(handle) 01947 END SUBROUTINE calculate_rho_elec 01948 01949 ! ***************************************************************************** 01954 SUBROUTINE calculate_drho_elec(matrix_p,drho,drho_gspace,qs_env,& 01955 soft_valid,basis_set_id,error) 01956 01957 TYPE(cp_dbcsr_type), POINTER :: matrix_p 01958 TYPE(pw_p_type), DIMENSION(:), 01959 INTENT(INOUT) :: drho, drho_gspace 01960 TYPE(qs_environment_type), POINTER :: qs_env 01961 LOGICAL, INTENT(IN), OPTIONAL :: soft_valid 01962 INTEGER, INTENT(IN), OPTIONAL :: basis_set_id 01963 TYPE(cp_error_type), INTENT(inout) :: error 01964 01965 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_drho_elec', 01966 routineP = moduleN//':'//routineN 01967 INTEGER, PARAMETER :: add_tasks = 1000, 01968 max_tasks = 2000 01969 REAL(kind=dp), PARAMETER :: mult_tasks = 2.0_dp 01970 01971 INTEGER :: bcol, brow, handle, i, iatom, iatom_old, idir, igrid_level, 01972 ikind, ikind_old, ipgf, iset, iset_old, itask, ithread, jatom, 01973 jatom_old, jkind, jkind_old, jpgf, jset, jset_old, maxco, maxpgf, 01974 maxset, maxsgf, maxsgf_set, my_basis_set_id, na1, na2, natoms, nb1, 01975 nb2, ncoa, ncob, nkind, nseta, nsetb, ntasks, nthread, sgfa, sgfb, stat 01976 INTEGER(kind=int_8), DIMENSION(:), 01977 POINTER :: atom_pair_recv, atom_pair_send 01978 INTEGER(kind=int_8), DIMENSION(:, :), 01979 POINTER :: tasks 01980 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, 01981 lb_min, npgfa, npgfb, nsgfa, 01982 nsgfb 01983 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 01984 LOGICAL :: atom_pair_changed, distributed_rs_grids, failure, found, 01985 map_consistent, my_soft, use_subpatch 01986 REAL(KIND=dp) :: eps_rho_rspace, rab2, scale, 01987 zetp 01988 REAL(KIND=dp), DIMENSION(3) :: ra, rab, rab_inv, rb 01989 REAL(KIND=dp), DIMENSION(:, :), POINTER :: dist_ab, p_block, pab, 01990 sphi_a, sphi_b, work, zeta, 01991 zetb 01992 REAL(KIND=dp), DIMENSION(:, :, :), 01993 POINTER :: pabt, workt 01994 TYPE(atomic_kind_type), DIMENSION(:), 01995 POINTER :: atomic_kind_set 01996 TYPE(atomic_kind_type), POINTER :: atomic_kind 01997 TYPE(cell_type), POINTER :: cell 01998 TYPE(cp_dbcsr_type), POINTER :: deltap 01999 TYPE(cp_para_env_type), POINTER :: para_env 02000 TYPE(cube_info_type), DIMENSION(:), 02001 POINTER :: cube_info 02002 TYPE(dft_control_type), POINTER :: dft_control 02003 TYPE(gridlevel_info_type), POINTER :: gridlevel_info 02004 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 02005 TYPE(neighbor_list_set_p_type), 02006 DIMENSION(:), POINTER :: sab_orb 02007 TYPE(particle_type), DIMENSION(:), 02008 POINTER :: particle_set 02009 TYPE(pw_env_type), POINTER :: pw_env 02010 TYPE(realspace_grid_desc_p_type), 02011 DIMENSION(:), POINTER :: rs_descs 02012 TYPE(realspace_grid_p_type), 02013 DIMENSION(:), POINTER :: rs_rho 02014 TYPE(section_vals_type), POINTER :: input, interp_section 02015 TYPE(task_list_type), POINTER :: task_list, task_list_soft 02016 02017 failure=.FALSE. 02018 NULLIFY(atomic_kind, cell, dft_control, orb_basis_set, deltap, & 02019 atomic_kind_set, sab_orb, particle_set, rs_rho, pw_env, rs_descs, & 02020 para_env, dist_ab, la_max, la_min, & 02021 lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, p_block, & 02022 sphi_a, sphi_b, zeta, zetb, first_sgfa, first_sgfb, tasks, pabt, & 02023 workt) 02024 02025 debug_count=debug_count+1 02026 02027 ! by default, the full density is calculated 02028 my_soft=.FALSE. 02029 IF (PRESENT(soft_valid)) my_soft = soft_valid 02030 02031 CALL timeset(routineN,handle) 02032 ! 02033 IF( PRESENT(basis_set_id) ) THEN 02034 my_basis_set_id = basis_set_id 02035 ELSE 02036 my_basis_set_id = use_orb_basis_set 02037 END IF 02038 02039 SELECT CASE (my_basis_set_id) 02040 CASE (use_orb_basis_set) 02041 CALL get_qs_env(qs_env=qs_env,& 02042 atomic_kind_set=atomic_kind_set,& 02043 cell=cell,& 02044 dft_control=dft_control,& 02045 particle_set=particle_set,& 02046 sab_orb=sab_orb,& 02047 para_env=para_env,& 02048 task_list=task_list,& 02049 task_list_soft=task_list_soft,& 02050 input=input,& 02051 pw_env=pw_env,error=error) 02052 CASE (use_aux_fit_basis_set) 02053 CALL get_qs_env(qs_env=qs_env,& 02054 atomic_kind_set=atomic_kind_set,& 02055 cell=cell,& 02056 dft_control=dft_control,& 02057 particle_set=particle_set,& 02058 sab_aux_fit=sab_orb,& 02059 para_env=para_env,& 02060 task_list_aux_fit=task_list,& 02061 task_list_soft=task_list_soft,& 02062 input=input,& 02063 pw_env=pw_env,error=error) 02064 END SELECT 02065 02066 ! *** assign from pw_env 02067 gridlevel_info=>pw_env%gridlevel_info 02068 cube_info=>pw_env%cube_info 02069 02070 ! *** Allocate work storage *** 02071 nthread = 1 02072 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 02073 maxco=maxco,& 02074 maxsgf=maxsgf,& 02075 maxsgf_set=maxsgf_set,& 02076 basis_set_id=my_basis_set_id) 02077 nkind = SIZE(atomic_kind_set) 02078 CALL reallocate(pabt,1,maxco,1,maxco,0,nthread-1) 02079 CALL reallocate(workt,1,maxco,1,maxsgf_set,0,nthread-1) 02080 02081 ! find maximum numbers 02082 nkind = SIZE(atomic_kind_set) 02083 natoms = SIZE( particle_set ) 02084 maxset=0 02085 maxpgf=0 02086 DO ikind=1,nkind 02087 atomic_kind => atomic_kind_set(ikind) 02088 SELECT CASE( my_basis_set_id) 02089 CASE (use_orb_basis_set) 02090 CALL get_atomic_kind(atomic_kind=atomic_kind,& 02091 softb = my_soft, & 02092 orb_basis_set=orb_basis_set) 02093 CASE (use_aux_fit_basis_set) 02094 CALL get_atomic_kind(atomic_kind=atomic_kind,& 02095 softb = my_soft, & 02096 aux_fit_basis_set=orb_basis_set,& 02097 basis_set_id=my_basis_set_id) 02098 END SELECT 02099 02100 IF (.NOT.ASSOCIATED(orb_basis_set)) CYCLE 02101 02102 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 02103 npgf=npgfa, nset=nseta ) 02104 02105 maxset=MAX(nseta,maxset) 02106 maxpgf=MAX(MAXVAL(npgfa),maxpgf) 02107 END DO 02108 02109 ! get the task lists 02110 IF (my_soft) task_list=>task_list_soft 02111 CPPrecondition(ASSOCIATED(task_list),cp_failure_level,routineP,error,failure) 02112 tasks =>task_list%tasks 02113 dist_ab=>task_list%dist_ab 02114 atom_pair_send=>task_list%atom_pair_send 02115 atom_pair_recv=>task_list%atom_pair_recv 02116 ntasks=task_list%ntasks 02117 02118 ! *** set up the rs multi-grids 02119 CPPrecondition(ASSOCIATED(pw_env),cp_failure_level,routineP,error,failure) 02120 CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho, error=error) 02121 DO igrid_level=1,gridlevel_info%ngrid_levels 02122 CALL rs_grid_retain(rs_rho(igrid_level)%rs_grid, error=error) 02123 distributed_rs_grids=rs_rho(igrid_level)%rs_grid%desc%distributed 02124 END DO 02125 02126 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 02127 map_consistent = dft_control%qs_control%map_consistent 02128 ! *** Initialize working density matrix *** 02129 ! distributed rs grids require a matrix that will be changed 02130 ! whereas this is not the case for replicated grids 02131 IF (distributed_rs_grids) THEN 02132 NULLIFY(deltap) 02133 ALLOCATE(deltap) 02134 CALL cp_dbcsr_init(deltap, error=error) 02135 CALL cp_dbcsr_copy(deltap,matrix_p,name="DeltaP",error=error) 02136 !CALL replicate_matrix(matrix_p,deltap,target_name="DeltaP",error=error) 02137 ! this matrix has no strict sparsity pattern in parallel 02138 !deltap%sparsity_id=-1 02139 ELSE 02140 deltap=>matrix_p 02141 ENDIF 02142 02143 ! distribute the matrix 02144 IF (distributed_rs_grids) THEN 02145 CALL rs_distribute_matrix (rs_descs, deltap, atom_pair_send, atom_pair_recv, & 02146 natoms, scatter=.TRUE., error=error) 02147 ENDIF 02148 02149 ! map all tasks on the grids 02150 02151 ithread = 0 02152 pab => pabt(:,:,ithread) 02153 work => workt(:,:,ithread) 02154 02155 loop_xyz: DO idir=1,3 02156 02157 DO igrid_level=1,gridlevel_info%ngrid_levels 02158 CALL rs_grid_zero(rs_rho(igrid_level)%rs_grid) 02159 END DO 02160 02161 iatom_old = -1 ; jatom_old = -1 ; iset_old = -1 ; jset_old = -1 02162 ikind_old = -1 ; jkind_old = -1 02163 loop_tasks: DO itask=1,ntasks 02164 02165 !decode the atom pair and basis info 02166 CALL int2pair(tasks(3,itask),igrid_level,iatom,jatom,iset,jset,ipgf,jpgf,natoms,maxset,maxpgf) 02167 ikind = particle_set(iatom)%atomic_kind%kind_number 02168 jkind = particle_set(jatom)%atomic_kind%kind_number 02169 02170 IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old) THEN 02171 02172 IF (iatom .NE. iatom_old) ra(:) = pbc(particle_set(iatom)%r,cell) 02173 02174 IF (iatom <= jatom) THEN 02175 brow = iatom 02176 bcol = jatom 02177 ELSE 02178 brow = jatom 02179 bcol = iatom 02180 END IF 02181 02182 IF(ikind .NE. ikind_old ) THEN 02183 SELECT CASE (my_basis_set_id) 02184 CASE (use_orb_basis_set) 02185 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& 02186 softb = my_soft, & 02187 orb_basis_set=orb_basis_set) 02188 CASE (use_aux_fit_basis_set) 02189 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& 02190 softb = my_soft, & 02191 aux_fit_basis_set=orb_basis_set,& 02192 basis_set_id=my_basis_set_id) 02193 END SELECT 02194 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 02195 first_sgf=first_sgfa,& 02196 lmax=la_max,& 02197 lmin=la_min,& 02198 npgf=npgfa,& 02199 nset=nseta,& 02200 nsgf_set=nsgfa,& 02201 sphi=sphi_a,& 02202 zet=zeta) 02203 ENDIF 02204 02205 IF (jkind .NE. jkind_old ) THEN 02206 SELECT CASE (my_basis_set_id) 02207 CASE (use_orb_basis_set) 02208 CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind,& 02209 softb = my_soft, & 02210 orb_basis_set=orb_basis_set) 02211 CASE (use_aux_fit_basis_set) 02212 CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind,& 02213 softb = my_soft, & 02214 aux_fit_basis_set=orb_basis_set,& 02215 basis_set_id=my_basis_set_id) 02216 END SELECT 02217 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 02218 first_sgf=first_sgfb,& 02219 lmax=lb_max,& 02220 lmin=lb_min,& 02221 npgf=npgfb,& 02222 nset=nsetb,& 02223 nsgf_set=nsgfb,& 02224 sphi=sphi_b,& 02225 zet=zetb) 02226 ENDIF 02227 02228 CALL cp_dbcsr_get_block_p(matrix=deltap,& 02229 row=brow,col=bcol,BLOCK=p_block,found=found) 02230 IF (.NOT.ASSOCIATED(p_block)) & 02231 CALL stop_program(routineN,moduleN,__LINE__,& 02232 "p_block not associated in deltap") 02233 02234 iatom_old = iatom 02235 jatom_old = jatom 02236 ikind_old = ikind 02237 jkind_old = jkind 02238 atom_pair_changed = .TRUE. 02239 02240 ELSE 02241 02242 atom_pair_changed = .FALSE. 02243 02244 ENDIF 02245 02246 IF (atom_pair_changed .OR. iset_old .NE. iset .OR. jset_old .NE. jset) THEN 02247 02248 ncoa = npgfa(iset)*ncoset(la_max(iset)) 02249 sgfa = first_sgfa(1,iset) 02250 ncob = npgfb(jset)*ncoset(lb_max(jset)) 02251 sgfb = first_sgfb(1,jset) 02252 02253 IF (iatom <= jatom) THEN 02254 CALL dgemm("N","N",ncoa,nsgfb(jset),nsgfa(iset),& 02255 1.0_dp,sphi_a(1,sgfa),SIZE(sphi_a,1),& 02256 p_block(sgfa,sgfb),SIZE(p_block,1),& 02257 0.0_dp,work(1,1),maxco) 02258 CALL dgemm("N","T",ncoa,ncob,nsgfb(jset),& 02259 1.0_dp,work(1,1),maxco,& 02260 sphi_b(1,sgfb),SIZE(sphi_b,1),& 02261 0.0_dp,pab(1,1),maxco) 02262 ELSE 02263 CALL dgemm("N","N",ncob,nsgfa(iset),nsgfb(jset),& 02264 1.0_dp,sphi_b(1,sgfb),SIZE(sphi_b,1),& 02265 p_block(sgfb,sgfa),SIZE(p_block,1),& 02266 0.0_dp,work(1,1),maxco) 02267 CALL dgemm("N","T",ncob,ncoa,nsgfa(iset),& 02268 1.0_dp,work(1,1),maxco,& 02269 sphi_a(1,sgfa),SIZE(sphi_a,1),& 02270 0.0_dp,pab(1,1),maxco) 02271 END IF 02272 02273 iset_old = iset 02274 jset_old = jset 02275 02276 ENDIF 02277 02278 rab(:) = dist_ab (:,itask) 02279 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 02280 rb(:) = ra(:) + rab(:) 02281 zetp = zeta(ipgf,iset) + zetb(jpgf,jset) 02282 02283 na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1 02284 na2 = ipgf*ncoset(la_max(iset)) 02285 nb1 = (jpgf - 1)*ncoset(lb_max(jset)) + 1 02286 nb2 = jpgf*ncoset(lb_max(jset)) 02287 02288 ! takes the density matrix symmetry in account, i.e. off-diagonal blocks need to be mapped 'twice' 02289 IF (iatom == jatom) THEN 02290 scale = 1.0_dp 02291 ELSE 02292 scale = 2.0_dp 02293 END IF 02294 02295 ! check whether we need to use fawzi's generalised collocation scheme 02296 IF(rs_rho(igrid_level)%rs_grid%desc%distributed)THEN 02297 !tasks(4,:) is 0 for replicated, 1 for distributed 2 for exceptional distributed tasks 02298 IF (tasks(4,itask) .EQ. 2 ) THEN 02299 use_subpatch = .TRUE. 02300 ELSE 02301 use_subpatch = .FALSE. 02302 ENDIF 02303 ELSE 02304 use_subpatch = .FALSE. 02305 ENDIF 02306 IF (iatom <= jatom) THEN 02307 CALL collocate_pgf_product_rspace(& 02308 la_max(iset),zeta(ipgf,iset),la_min(iset),& 02309 lb_max(jset),zetb(jpgf,jset),lb_min(jset),& 02310 ra,rab,rab2,scale,pab,na1-1,nb1-1,& 02311 rs_rho(igrid_level)%rs_grid,cell,cube_info(igrid_level),& 02312 eps_rho_rspace,& 02313 ga_gb_function=FUNC_DABpADB, & 02314 idir=idir, & 02315 map_consistent=map_consistent,use_subpatch=use_subpatch,subpatch_pattern=tasks(6,itask),error=error) 02316 ELSE 02317 rab_inv=-rab 02318 CALL collocate_pgf_product_rspace(& 02319 lb_max(jset),zetb(jpgf,jset),lb_min(jset),& 02320 la_max(iset),zeta(ipgf,iset),la_min(iset),& 02321 rb,rab_inv,rab2,scale,pab,nb1-1,na1-1,& 02322 rs_rho(igrid_level)%rs_grid,cell,cube_info(igrid_level),& 02323 eps_rho_rspace,& 02324 ga_gb_function=FUNC_DABpADB, & 02325 idir=idir, & 02326 map_consistent=map_consistent,use_subpatch=use_subpatch,subpatch_pattern=tasks(6,itask),error=error) 02327 END IF 02328 02329 END DO loop_tasks 02330 02331 interp_section => section_vals_get_subs_vals(input,"DFT%MGRID%INTERPOLATOR",error=error) 02332 CALL density_rs2pw_basic(pw_env,rs_rho,drho(idir),drho_gspace(idir),& 02333 interp_section=interp_section,error=error) 02334 02335 END DO loop_xyz 02336 02337 ! *** Release work storage *** 02338 IF (ASSOCIATED(rs_rho)) THEN 02339 DO i=1, SIZE(rs_rho) 02340 CALL rs_grid_release(rs_rho(i)%rs_grid, error=error) 02341 END DO 02342 END IF 02343 IF (distributed_rs_grids) THEN 02344 CALL cp_dbcsr_deallocate_matrix ( deltap ,error=error) 02345 ENDIF 02346 02347 DEALLOCATE (pabt,STAT=stat) 02348 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"pabt") 02349 DEALLOCATE (workt,STAT=stat) 02350 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"workt") 02351 02352 CALL timestop(handle) 02353 02354 END SUBROUTINE calculate_drho_elec 02355 02356 ! ***************************************************************************** 02369 SUBROUTINE calculate_wavefunction(mo_vectors,ivector,rho,rho_gspace, & 02370 atomic_kind_set,cell,dft_control,particle_set, & 02371 pw_env,basis_set_id, external_vector, error) 02372 02373 TYPE(cp_fm_type), POINTER :: mo_vectors 02374 INTEGER :: ivector 02375 TYPE(pw_p_type), INTENT(INOUT) :: rho, rho_gspace 02376 TYPE(atomic_kind_type), DIMENSION(:), 02377 POINTER :: atomic_kind_set 02378 TYPE(cell_type), POINTER :: cell 02379 TYPE(dft_control_type), POINTER :: dft_control 02380 TYPE(particle_type), DIMENSION(:), 02381 POINTER :: particle_set 02382 TYPE(pw_env_type), POINTER :: pw_env 02383 INTEGER, INTENT(IN), OPTIONAL :: basis_set_id 02384 REAL(KIND=dp), DIMENSION(:), OPTIONAL :: external_vector 02385 TYPE(cp_error_type), INTENT(inout) :: error 02386 02387 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_wavefunction', 02388 routineP = moduleN//':'//routineN 02389 02390 INTEGER :: dir, group, group_size, handle, i, iatom, igrid_level, ipgf, 02391 iset, maxco, maxsgf_set, my_basis_set_id, my_pos, na1, na2, nao, natom, 02392 ncoa, nseta, offset, sgfa, stat 02393 INTEGER, ALLOCATABLE, DIMENSION(:) :: where_is_the_point 02394 INTEGER, DIMENSION(3) :: lb, location, tp, ub 02395 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa 02396 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa 02397 LOGICAL :: failure, local, map_it_here 02398 REAL(KIND=dp) :: dab, eps_rho_rspace, rab2, 02399 scale, zetp 02400 REAL(KIND=dp), DIMENSION(3) :: ra, rab 02401 REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvector 02402 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab, sphi_a, work, zeta 02403 TYPE(cube_info_type), DIMENSION(:), 02404 POINTER :: cube_info 02405 TYPE(gridlevel_info_type), POINTER :: gridlevel_info 02406 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 02407 TYPE(pw_p_type), DIMENSION(:), POINTER :: mgrid_gspace, mgrid_rspace 02408 TYPE(pw_pool_p_type), DIMENSION(:), 02409 POINTER :: pw_pools 02410 TYPE(realspace_grid_desc_p_type), 02411 DIMENSION(:), POINTER :: rs_descs 02412 TYPE(REALSPACE_GRID_P_TYPE), 02413 DIMENSION(:), POINTER :: rs_rho 02414 02415 IF( PRESENT( basis_set_id ) ) THEN 02416 my_basis_set_id = basis_set_id 02417 ELSE 02418 my_basis_set_id = use_orb_basis_set 02419 END IF 02420 02421 failure=.FALSE. 02422 02423 CALL timeset(routineN,handle) 02424 02425 NULLIFY(eigenvector, orb_basis_set,& 02426 pab,work,la_max, la_min,& 02427 npgfa, nsgfa, & 02428 sphi_a, zeta, first_sgfa,& 02429 rs_rho,rs_descs,pw_pools,mgrid_rspace,mgrid_gspace) 02430 02431 IF (PRESENT(external_vector)) THEN 02432 nao=SIZE(external_vector) 02433 ALLOCATE (eigenvector(nao),STAT=stat) 02434 IF (stat.NE.0) CALL stop_memory(routineN,moduleN,__LINE__,"eigenvector",dp_size*nao) 02435 eigenvector=external_vector 02436 ELSE 02437 CALL cp_fm_get_info(matrix=mo_vectors,nrow_global=nao, error=error) 02438 ALLOCATE (eigenvector(nao),STAT=stat) 02439 IF (stat.NE.0) CALL stop_memory(routineN,moduleN,__LINE__,"eigenvector",dp_size*nao) 02440 DO i=1,nao 02441 CALL cp_fm_get_element(mo_vectors,i,ivector,eigenvector(i),local) 02442 ENDDO 02443 ENDIF 02444 02445 ! *** set up the pw multi-grids 02446 CPPrecondition(ASSOCIATED(pw_env),cp_failure_level,routineP,error,failure) 02447 CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho, pw_pools=pw_pools, & 02448 cube_info=cube_info, gridlevel_info=gridlevel_info, error=error) 02449 02450 CALL pw_pools_create_pws(pw_pools,mgrid_gspace,& 02451 use_data = COMPLEXDATA1D,& 02452 in_space = RECIPROCALSPACE, error=error) 02453 CALL pw_pools_create_pws(pw_pools,mgrid_rspace,& 02454 use_data = REALDATA3D,& 02455 in_space = REALSPACE, error=error) 02456 02457 ! *** set up rs multi-grids 02458 DO igrid_level=1,gridlevel_info%ngrid_levels 02459 CALL rs_grid_retain(rs_rho(igrid_level)%rs_grid, error=error) 02460 CALL rs_grid_zero(rs_rho(igrid_level)%rs_grid) 02461 END DO 02462 02463 02464 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 02465 ! *** Allocate work storage *** 02466 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 02467 maxco=maxco,& 02468 natom=natom,& 02469 maxsgf_set=maxsgf_set,& 02470 basis_set_id=my_basis_set_id) 02471 02472 ALLOCATE (pab(maxco,1),STAT=stat) 02473 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 02474 "pab",maxco*dp_size) 02475 ALLOCATE (work(maxco,1),STAT=stat) 02476 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 02477 "work",maxco*dp_size) 02478 02479 offset=0 02480 group=mgrid_rspace(1)%pw%pw_grid%para%group 02481 my_pos=mgrid_rspace(1)%pw%pw_grid%para%my_pos 02482 group_size=mgrid_rspace(1)%pw%pw_grid%para%group_size 02483 ALLOCATE(where_is_the_point(0:group_size-1)) 02484 02485 DO iatom=1,natom 02486 SELECT CASE (my_basis_set_id) 02487 CASE (use_orb_basis_set) 02488 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& 02489 orb_basis_set=orb_basis_set) 02490 CASE (use_aux_fit_basis_set) 02491 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& 02492 aux_fit_basis_set=orb_basis_set) 02493 CASE (use_aux_basis_set) 02494 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& 02495 aux_basis_set=orb_basis_set) 02496 CASE (use_ri_aux_basis_set) 02497 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& 02498 ri_aux_basis_set=orb_basis_set) 02499 END SELECT 02500 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 02501 first_sgf=first_sgfa,& 02502 lmax=la_max,& 02503 lmin=la_min,& 02504 npgf=npgfa,& 02505 nset=nseta,& 02506 nsgf_set=nsgfa,& 02507 sphi=sphi_a,& 02508 zet=zeta) 02509 ra(:) = pbc(particle_set(iatom)%r,cell) 02510 rab(:) = 0.0_dp 02511 rab2 = 0.0_dp 02512 dab = 0.0_dp 02513 02514 DO iset=1,nseta 02515 02516 ncoa = npgfa(iset)*ncoset(la_max(iset)) 02517 sgfa = first_sgfa(1,iset) 02518 02519 DO i=1,nsgfa(iset) 02520 work(i,1)=eigenvector(offset+i) 02521 ENDDO 02522 02523 CALL dgemm("N","N",ncoa,1,nsgfa(iset),& 02524 1.0_dp,sphi_a(1,sgfa),SIZE(sphi_a,1),& 02525 work(1,1),SIZE(work,1),& 02526 0.0_dp,pab(1,1),SIZE(pab,1)) 02527 02528 DO ipgf=1,npgfa(iset) 02529 02530 na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1 02531 na2 = ipgf*ncoset(la_max(iset)) 02532 02533 scale = 1.0_dp 02534 zetp = zeta(ipgf,iset) 02535 igrid_level = gaussian_gridlevel(gridlevel_info,zetp) 02536 02537 map_it_here=.FALSE. 02538 02539 IF (.NOT. ALL (rs_rho(igrid_level)%rs_grid%desc%perd == 1)) THEN 02540 DO dir = 1,3 02541 ! bounds of local grid (i.e. removing the 'wings'), if periodic 02542 tp(dir) = FLOOR(DOT_PRODUCT(cell%h_inv(dir,:),ra)*rs_rho(igrid_level)%rs_grid%desc%npts(dir)) 02543 tp(dir) = MODULO ( tp(dir), rs_rho(igrid_level)%rs_grid%desc%npts(dir) ) 02544 IF (rs_rho(igrid_level)%rs_grid%desc%perd(dir) .NE. 1) THEN 02545 lb(dir) = rs_rho(igrid_level)%rs_grid%lb_local ( dir ) + rs_rho(igrid_level)%rs_grid%desc%border 02546 ub(dir) = rs_rho(igrid_level)%rs_grid%ub_local ( dir ) - rs_rho(igrid_level)%rs_grid%desc%border 02547 ELSE 02548 lb(dir) = rs_rho(igrid_level)%rs_grid%lb_local ( dir ) 02549 ub(dir) = rs_rho(igrid_level)%rs_grid%ub_local ( dir ) 02550 ENDIF 02551 ! distributed grid, only map if it is local to the grid 02552 location(dir)=tp(dir)+rs_rho(igrid_level)%rs_grid%desc%lb(dir) 02553 ENDDO 02554 IF (lb(1)<=location(1) .AND. location(1)<=ub(1) .AND. & 02555 lb(2)<=location(2) .AND. location(2)<=ub(2) .AND. & 02556 lb(3)<=location(3) .AND. location(3)<=ub(3)) THEN 02557 map_it_here=.TRUE. 02558 ENDIF 02559 ELSE 02560 ! not distributed, just a round-robin distribution over the full set of CPUs 02561 IF (MODULO(offset,group_size)==my_pos) map_it_here=.TRUE. 02562 ENDIF 02563 02564 IF (map_it_here) CALL collocate_pgf_product_rspace(& 02565 la_max(iset),zeta(ipgf,iset),la_min(iset),& 02566 0,0.0_dp,0,& 02567 ra,rab,rab2,scale,pab,na1-1,0,& 02568 rs_rho(igrid_level)%rs_grid,cell,cube_info(igrid_level),& 02569 eps_rho_rspace,map_consistent=.TRUE.,ga_gb_function=FUNC_AB,error=error) 02570 02571 END DO 02572 02573 offset=offset+nsgfa(iset) 02574 02575 END DO 02576 02577 END DO 02578 02579 DO igrid_level=1,gridlevel_info%ngrid_levels 02580 CALL rs_pw_transfer(rs_rho(igrid_level)%rs_grid,& 02581 mgrid_rspace(igrid_level)%pw,rs2pw,error=error) 02582 CALL rs_grid_release(rs_rho(igrid_level)%rs_grid, error=error) 02583 ENDDO 02584 02585 CALL pw_zero(rho_gspace%pw,error=error) 02586 DO igrid_level=1,gridlevel_info%ngrid_levels 02587 CALL pw_transfer(mgrid_rspace(igrid_level)%pw,& 02588 mgrid_gspace(igrid_level)%pw,error=error) 02589 CALL pw_axpy(mgrid_gspace(igrid_level)%pw,rho_gspace%pw,& 02590 error=error) 02591 END DO 02592 02593 CALL pw_transfer(rho_gspace%pw,rho%pw,error=error) 02594 02595 ! Release work storage 02596 DEALLOCATE (eigenvector,STAT=stat) 02597 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"eigenvector") 02598 02599 DEALLOCATE (pab,STAT=stat) 02600 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"pab") 02601 02602 DEALLOCATE (work,STAT=stat) 02603 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"work") 02604 02605 ! give back the pw multi-grids 02606 CALL pw_pools_give_back_pws(pw_pools,mgrid_gspace,error=error) 02607 CALL pw_pools_give_back_pws(pw_pools,mgrid_rspace,error=error) 02608 02609 CALL timestop(handle) 02610 02611 END SUBROUTINE calculate_wavefunction 02612 02613 ! ***************************************************************************** 02616 SUBROUTINE collocate_pgf_product_rspace(la_max,zeta,la_min,& 02617 lb_max,zetb,lb_min,& 02618 ra,rab,rab2,scale,pab,o1,o2,& 02619 rsgrid,cell,cube_info,& 02620 eps_rho_rspace,ga_gb_function,& 02621 lgrid,ithread,& 02622 map_consistent,& 02623 collocate_rho0,& 02624 rpgf0_s,idir,ir,rsgauge,rsbuf,& 02625 use_subpatch,subpatch_pattern,error) 02626 02627 INTEGER, INTENT(IN) :: la_max 02628 REAL(KIND=dp), INTENT(IN) :: zeta 02629 INTEGER, INTENT(IN) :: la_min, lb_max 02630 REAL(KIND=dp), INTENT(IN) :: zetb 02631 INTEGER, INTENT(IN) :: lb_min 02632 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rab 02633 REAL(KIND=dp), INTENT(IN) :: rab2, scale 02634 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab 02635 INTEGER, INTENT(IN) :: o1, o2 02636 TYPE(realspace_grid_type), POINTER :: rsgrid 02637 TYPE(cell_type), POINTER :: cell 02638 TYPE(cube_info_type), INTENT(IN) :: cube_info 02639 REAL(KIND=dp), INTENT(IN) :: eps_rho_rspace 02640 INTEGER, INTENT(IN) :: ga_gb_function 02641 TYPE(lgrid_type), OPTIONAL :: lgrid 02642 INTEGER, INTENT(IN), OPTIONAL :: ithread 02643 LOGICAL, INTENT(IN), OPTIONAL :: map_consistent, collocate_rho0 02644 REAL(dp), INTENT(IN), OPTIONAL :: rpgf0_s 02645 INTEGER, INTENT(IN), OPTIONAL :: idir, ir 02646 TYPE(realspace_grid_type), POINTER, OPTIONAL :: rsgauge,rsbuf 02647 LOGICAL, OPTIONAL :: use_subpatch 02648 INTEGER(KIND=int_8), OPTIONAL, INTENT(IN):: subpatch_pattern 02649 TYPE(cp_error_type), INTENT(INOUT) :: error 02650 02651 CHARACTER(len=*), PARAMETER :: routineN = 'collocate_pgf_product_rspace', 02652 routineP = moduleN//':'//routineN 02653 02654 INTEGER :: cmax, gridbounds(2,3), i, ico, icoef, ider1, ider2, ig, ithread_l, 02655 jco, k, l, la_max_local, la_min_local, lb_max_local, lb_min_local, 02656 length, lx, lx_max, lxa, lxb, lxy, lxy_max, lxyz, lxyz_max, lya, lyb, 02657 lza, lzb, o1_local, o2_local, offset, start 02658 INTEGER, DIMENSION(3) :: cubecenter, lb_cube, ng, 02659 ub_cube 02660 INTEGER, DIMENSION(:), POINTER :: ly_max, lz_max, sphere_bounds 02661 LOGICAL :: failure, my_collocate_rho0, 02662 my_map_consistent,subpatch_collocate 02663 REAL(KIND=dp) :: a, b, binomial_k_lxa, binomial_l_lxb, cutoff, f, pg, 02664 prefactor, radius, rpg, ya, yap, yb, ybp, za, zap, zb, zbp, zetp 02665 REAL(KIND=dp), DIMENSION(3) :: dr, rap, rb, rbp, roffset, rp 02666 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab_local 02667 REAL(KIND=dp), DIMENSION(:, :, :), 02668 POINTER :: grid 02669 02670 INTEGER :: lxp,lyp,lzp,lp,lxpm,lypm,iaxis 02671 INTEGER, ALLOCATABLE, DIMENSION(:,:) :: map 02672 REAL(kind=dp) :: p_ele,ax,ay,az 02673 REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:,:,:) :: alpha 02674 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coef_xyz 02675 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coef_xyt 02676 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coef_xtt 02677 02678 REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:,:) :: pol_z 02679 REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:,:) :: pol_y 02680 REAL(kind=dp), ALLOCATABLE, DIMENSION(:,:) :: pol_x 02681 REAL(KIND=dp) :: t_exp_1,t_exp_2,t_exp_min_1,t_exp_min_2,t_exp_plus_1,t_exp_plus_2 02682 02683 REAL(kind=dp), POINTER, DIMENSION(:,:,:) :: grid_tmp,gauge 02684 02685 failure = .FALSE. 02686 02687 subpatch_collocate = .FALSE. 02688 02689 IF(PRESENT(use_subpatch)) THEN 02690 IF(use_subpatch)THEN 02691 subpatch_collocate = .TRUE. 02692 CPPrecondition(PRESENT(subpatch_pattern),cp_failure_level,routineP,error,failure) 02693 ENDIF 02694 ENDIF 02695 02696 IF (PRESENT(ithread)) THEN 02697 ithread_l=ithread 02698 ELSE 02699 ithread_l=0 02700 ENDIF 02701 02702 ! use identical radii for integrate and collocate ? 02703 IF (PRESENT(map_consistent)) THEN 02704 my_map_consistent=map_consistent 02705 ELSE 02706 my_map_consistent=.FALSE. 02707 ENDIF 02708 02709 IF (PRESENT(collocate_rho0).AND.PRESENT(rpgf0_s)) THEN 02710 my_collocate_rho0=collocate_rho0 02711 ELSE 02712 my_collocate_rho0=.FALSE. 02713 END IF 02714 02715 zetp = zeta + zetb 02716 f = zetb/zetp 02717 rap(:) = f*rab(:) 02718 rbp(:) = rap(:) - rab(:) 02719 rp(:) = ra(:) + rap(:) 02720 rb(:) = ra(:)+rab(:) 02721 02722 IF (my_map_consistent) THEN 02723 cutoff = 1.0_dp 02724 prefactor = EXP(-zeta*f*rab2) 02725 radius=exp_radius_very_extended(la_min,la_max,lb_min,lb_max,ra=ra,rb=rb,rp=rp,& 02726 zetp=zetp,eps=eps_rho_rspace,& 02727 prefactor=prefactor,cutoff=cutoff) 02728 prefactor = scale*EXP(-zeta*f*rab2) 02729 ELSE IF (my_collocate_rho0) THEN 02730 cutoff = 0.0_dp 02731 prefactor = 1.0_dp 02732 radius = rpgf0_s 02733 ELSE 02734 cutoff = 0.0_dp 02735 prefactor = scale*EXP(-zeta*f*rab2) 02736 radius=exp_radius_very_extended(la_min,la_max,lb_min,lb_max,pab,o1,o2,ra,rb,rp,& 02737 zetp,eps_rho_rspace,prefactor,cutoff) 02738 ENDIF 02739 02740 IF (radius .EQ. 0.0_dp ) THEN 02741 RETURN 02742 END IF 02743 02744 ! it's a choice to compute lX_min/max, pab here, 02745 ! this way we get the same radius as we use for the corresponding density 02746 SELECT CASE (ga_gb_function) 02747 CASE(FUNC_DADB) 02748 la_max_local=la_max+1 02749 la_min_local=MAX(la_min-1,0) 02750 lb_max_local=lb_max+1 02751 lb_min_local=MAX(lb_min-1,0) 02752 ! create a new pab_local so that mapping pab_local with pgf_a pgf_b 02753 ! is equivalent to mapping pab with 0.5 * (nabla pgf_a) . (nabla pgf_b) 02754 ! (ddx pgf_a ) (ddx pgf_b) = (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x})*(lbx pgf_{b-1x} - 2*zetb*pgf_{b+1x}) 02755 ! cleaner would possibly be to touch pzyx directly (avoiding the following allocate) 02756 ALLOCATE(pab_local(ncoset(la_max_local),ncoset(lb_max_local))) 02757 pab_local = 0.0_dp 02758 DO lxa=0,la_max 02759 DO lxb=0,lb_max 02760 DO lya=0,la_max-lxa 02761 DO lyb=0,lb_max-lxb 02762 DO lza=MAX(la_min-lxa-lya,0),la_max-lxa-lya 02763 DO lzb=MAX(lb_min-lxb-lyb,0),lb_max-lxb-lyb 02764 02765 ! this element of pab results in 12 elements of pab_local 02766 CALL prepare_dadb(pab_local,pab,lxa,lya,lza,lxb,lyb,lzb,o1,o2,zeta,zetb) 02767 02768 ENDDO 02769 ENDDO 02770 ENDDO 02771 ENDDO 02772 ENDDO 02773 ENDDO 02774 o1_local=0 02775 o2_local=0 02776 pab_local=pab_local * 0.5_dp 02777 CASE(FUNC_ADBmDAB) 02778 CPPrecondition(PRESENT(idir),cp_failure_level,routineP,error,failure) 02779 la_max_local=la_max+1 02780 la_min_local=MAX(la_min-1,0) 02781 lb_max_local=lb_max+1 02782 lb_min_local=MAX(lb_min-1,0) 02783 ! create a new pab_local so that mapping pab_local with pgf_a pgf_b 02784 ! is equivalent to mapping pab with 02785 ! pgf_a (nabla_{idir} pgf_b) - (nabla_{idir} pgf_a) pgf_b 02786 ! ( pgf_a ) (ddx pgf_b) - (ddx pgf_a)( pgf_b ) = 02787 ! pgf_a *(lbx pgf_{b-1x} - 2*zetb*pgf_{b+1x}) - 02788 ! (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x}) pgf_b 02789 02790 ALLOCATE(pab_local(ncoset(la_max_local),ncoset(lb_max_local))) 02791 pab_local = 0.0_dp 02792 DO lxa=0,la_max 02793 DO lxb=0,lb_max 02794 DO lya=0,la_max-lxa 02795 DO lyb=0,lb_max-lxb 02796 DO lza=MAX(la_min-lxa-lya,0),la_max-lxa-lya 02797 DO lzb=MAX(lb_min-lxb-lyb,0),lb_max-lxb-lyb 02798 ! this element of pab results in 4 elements of pab_local 02799 CALL prepare_adb_m_dab(pab_local,pab,idir,& 02800 lxa,lya,lza,lxb,lyb,lzb,o1,o2,zeta,zetb) 02801 END DO 02802 END DO 02803 END DO 02804 END DO 02805 END DO 02806 END DO 02807 o1_local=0 02808 o2_local=0 02809 CASE(FUNC_DABpADB) 02810 CPPrecondition(PRESENT(idir),cp_failure_level,routineP,error,failure) 02811 la_max_local=la_max+1 02812 la_min_local=MAX(la_min-1,0) 02813 lb_max_local=lb_max+1 02814 lb_min_local=MAX(lb_min-1,0) 02815 ! create a new pab_local so that mapping pab_local with pgf_a pgf_b 02816 ! is equivalent to mapping pab with 02817 ! pgf_a (nabla_{idir} pgf_b) + (nabla_{idir} pgf_a) pgf_b 02818 ! ( pgf_a ) (ddx pgf_b) + (ddx pgf_a)( pgf_b ) = 02819 ! pgf_a *(lbx pgf_{b-1x} + 2*zetb*pgf_{b+1x}) + 02820 ! (lax pgf_{a-1x} + 2*zeta*pgf_{a+1x}) pgf_b 02821 02822 ALLOCATE(pab_local(ncoset(la_max_local),ncoset(lb_max_local))) 02823 pab_local = 0.0_dp 02824 DO lxa=0,la_max 02825 DO lxb=0,lb_max 02826 DO lya=0,la_max-lxa 02827 DO lyb=0,lb_max-lxb 02828 DO lza=MAX(la_min-lxa-lya,0),la_max-lxa-lya 02829 DO lzb=MAX(lb_min-lxb-lyb,0),lb_max-lxb-lyb 02830 ! this element of pab results in 4 elements of pab_local 02831 CALL prepare_dab_p_adb(pab_local,pab,idir,& 02832 lxa,lya,lza,lxb,lyb,lzb,o1,o2,zeta,zetb) 02833 END DO 02834 END DO 02835 END DO 02836 END DO 02837 END DO 02838 END DO 02839 o1_local=0 02840 o2_local=0 02841 CASE(FUNC_DX,FUNC_DY,FUNC_DZ) 02842 ider1 = ga_gb_function - 500 02843 la_max_local=la_max+1 02844 la_min_local=MAX(la_min-1,0) 02845 lb_max_local=lb_max+1 02846 lb_min_local=MAX(lb_min-1,0) 02847 ! create a new pab_local so that mapping pab_local with pgf_a pgf_b 02848 ! is equivalent to mapping pab with 02849 ! d_{ider1} pgf_a d_{ider1} pgf_b 02850 ! dx pgf_a dx pgf_b = 02851 ! (lax pgf_{a-1x})*(lbx pgf_{b-1x}) - 2*zetb*lax*pgf_{a-1x}*pgf_{b+1x} - 02852 ! lbx pgf_{b-1x}*2*zeta*pgf_{a+1x}+ 4*zeta*zetab*pgf_{a+1x}pgf_{b+1x} 02853 02854 ALLOCATE(pab_local(ncoset(la_max_local),ncoset(lb_max_local))) 02855 pab_local = 0.0_dp 02856 DO lxa=0,la_max 02857 DO lxb=0,lb_max 02858 DO lya=0,la_max-lxa 02859 DO lyb=0,lb_max-lxb 02860 DO lza=MAX(la_min-lxa-lya,0),la_max-lxa-lya 02861 DO lzb=MAX(lb_min-lxb-lyb,0),lb_max-lxb-lyb 02862 ! this element of pab results in 4 elements of pab_local 02863 CALL prepare_dIadIb(pab_local,pab,ider1,& 02864 lxa,lya,lza,lxb,lyb,lzb,o1,o2,zeta,zetb) 02865 END DO 02866 END DO 02867 END DO 02868 END DO 02869 END DO 02870 END DO 02871 o1_local=0 02872 o2_local=0 02873 CASE(FUNC_DXDY,FUNC_DYDZ,FUNC_DZDX) 02874 ider1 = ga_gb_function - 600 02875 ider2 = ga_gb_function - 600 + 1 02876 IF(ider2>3) ider2 = ider1-2 02877 la_max_local=la_max+2 02878 la_min_local=MAX(la_min-2,0) 02879 lb_max_local=lb_max+2 02880 lb_min_local=MAX(lb_min-2,0) 02881 ! create a new pab_local so that mapping pab_local with pgf_a pgf_b 02882 ! is equivalent to mapping pab with 02883 ! d_{ider1} pgf_a d_{ider1} pgf_b 02884 ALLOCATE(pab_local(ncoset(la_max_local),ncoset(lb_max_local))) 02885 pab_local = 0.0_dp 02886 DO lxa=0,la_max 02887 DO lxb=0,lb_max 02888 DO lya=0,la_max-lxa 02889 DO lyb=0,lb_max-lxb 02890 DO lza=MAX(la_min-lxa-lya,0),la_max-lxa-lya 02891 DO lzb=MAX(lb_min-lxb-lyb,0),lb_max-lxb-lyb 02892 ! this element of pab results in 16 elements of pab_local 02893 CALL prepare_dijadijb (pab_local,pab,ider1,ider2,& 02894 lxa,lya,lza,lxb,lyb,lzb,o1,o2,zeta,zetb) 02895 END DO 02896 END DO 02897 END DO 02898 END DO 02899 END DO 02900 END DO 02901 o1_local=0 02902 o2_local=0 02903 CASE(FUNC_DXDX,FUNC_DYDY,FUNC_DZDZ) 02904 ider1 = ga_gb_function - 603 02905 la_max_local=la_max+2 02906 la_min_local=MAX(la_min-2,0) 02907 lb_max_local=lb_max+2 02908 lb_min_local=MAX(lb_min-2,0) 02909 ! create a new pab_local so that mapping pab_local with pgf_a pgf_b 02910 ! is equivalent to mapping pab with 02911 ! dd_{ider1} pgf_a dd_{ider1} pgf_b 02912 02913 ALLOCATE(pab_local(ncoset(la_max_local),ncoset(lb_max_local))) 02914 pab_local = 0.0_dp 02915 DO lxa=0,la_max 02916 DO lxb=0,lb_max 02917 DO lya=0,la_max-lxa 02918 DO lyb=0,lb_max-lxb 02919 DO lza=MAX(la_min-lxa-lya,0),la_max-lxa-lya 02920 DO lzb=MAX(lb_min-lxb-lyb,0),lb_max-lxb-lyb 02921 ! this element of pab results in 9 elements of pab_local 02922 CALL prepare_diiadiib (pab_local,pab,ider1,& 02923 lxa,lya,lza,lxb,lyb,lzb,o1,o2,zeta,zetb) 02924 END DO 02925 END DO 02926 END DO 02927 END DO 02928 END DO 02929 END DO 02930 o1_local=0 02931 o2_local=0 02932 CASE(FUNC_ARDBmDARB) 02933 CPPrecondition(PRESENT(idir),cp_failure_level,routineP,error,failure) 02934 CPPrecondition(PRESENT(ir),cp_failure_level,routineP,error,failure) 02935 la_max_local=la_max+1 02936 la_min_local=MAX(la_min-1,0) 02937 lb_max_local=lb_max+2 02938 lb_min_local=MAX(lb_min-1,0) 02939 ! create a new pab_local so that mapping pab_local with pgf_a pgf_b 02940 ! is equivalent to mapping pab with 02941 ! pgf_a (r-Rb)_{ir} (nabla_{idir} pgf_b) - (nabla_{idir} pgf_a) (r-Rb)_{ir} pgf_b 02942 ! ( pgf_a )(r-Rb)_{ir} (ddx pgf_b) - (ddx pgf_a) (r-Rb)_{ir} ( pgf_b ) = 02943 ! pgf_a *(lbx pgf_{b-1x+1ir} - 2*zetb*pgf_{b+1x+1ir}) - 02944 ! (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x}) pgf_{b+1ir} 02945 02946 ALLOCATE(pab_local(ncoset(la_max_local),ncoset(lb_max_local))) 02947 pab_local = 0.0_dp 02948 DO lxa=0,la_max 02949 DO lxb=0,lb_max 02950 DO lya=0,la_max-lxa 02951 DO lyb=0,lb_max-lxb 02952 DO lza=MAX(la_min-lxa-lya,0),la_max-lxa-lya 02953 DO lzb=MAX(lb_min-lxb-lyb,0),lb_max-lxb-lyb 02954 02955 ! this element of pab results in 4 elements of pab_local 02956 CALL prepare_ardb_m_darb(pab_local,pab,idir,ir,& 02957 lxa,lya,lza,lxb,lyb,lzb,o1,o2,zeta,zetb) 02958 END DO 02959 END DO 02960 END DO 02961 END DO 02962 END DO 02963 END DO 02964 o1_local=0 02965 o2_local=0 02966 CASE(FUNC_ARB) 02967 CPPrecondition(PRESENT(ir),cp_failure_level,routineP,error,failure) 02968 la_max_local=la_max 02969 la_min_local=la_min 02970 lb_max_local=lb_max+1 02971 lb_min_local=lb_min 02972 ! create a new pab_local so that mapping pab_local with pgf_a pgf_b 02973 ! is equivalent to mapping pab with 02974 ! pgf_a (r-Rb)_{ir} pgf_b = pgf_a * pgf_{b+1ir} 02975 ALLOCATE(pab_local(ncoset(la_max_local),ncoset(lb_max_local))) 02976 pab_local = 0.0_dp 02977 DO lxa=0,la_max 02978 DO lxb=0,lb_max 02979 DO lya=0,la_max-lxa 02980 DO lyb=0,lb_max-lxb 02981 DO lza=MAX(la_min-lxa-lya,0),la_max-lxa-lya 02982 DO lzb=MAX(lb_min-lxb-lyb,0),lb_max-lxb-lyb 02983 ! this element of pab results in 4 elements of pab_local 02984 CALL prepare_arb(pab_local,pab,ir,lxa,lya,lza,lxb,lyb,lzb,o1,o2) 02985 END DO 02986 END DO 02987 END DO 02988 END DO 02989 END DO 02990 END DO 02991 o1_local=0 02992 o2_local=0 02993 CASE(FUNC_AB) 02994 la_max_local=la_max 02995 la_min_local=la_min 02996 lb_max_local=lb_max 02997 lb_min_local=lb_min 02998 pab_local => pab 02999 o1_local=o1 03000 o2_local=o2 03001 CASE DEFAULT 03002 CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) 03003 END SELECT 03004 03005 ng(:) = rsgrid%desc%npts(:) 03006 grid => rsgrid%r(:,:,:) 03007 gridbounds(1,1)=LBOUND(GRID,1) 03008 gridbounds(2,1)=UBOUND(GRID,1) 03009 gridbounds(1,2)=LBOUND(GRID,2) 03010 gridbounds(2,2)=UBOUND(GRID,2) 03011 gridbounds(1,3)=LBOUND(GRID,3) 03012 gridbounds(2,3)=UBOUND(GRID,3) 03013 03014 03015 IF(PRESENT(ir).AND.PRESENT(rsgauge).AND.PRESENT(rsbuf)) THEN 03016 grid => rsbuf%r(:,:,:) 03017 grid_tmp => rsgrid%r(:,:,:) 03018 gauge => rsgauge%r(:,:,:) 03019 ENDIF 03020 03021 ! *** initialise the coefficient matrix, we transform the sum 03022 ! 03023 ! sum_{lxa,lya,lza,lxb,lyb,lzb} P_{lxa,lya,lza,lxb,lyb,lzb} (x-a_x)**lxa (y-a_y)**lya (z-a_z)**lza (x-b_x)**lxb (y-a_y)**lya (z-a_z)**lza 03024 ! 03025 ! into 03026 ! 03027 ! sum_{lxp,lyp,lzp} P_{lxp,lyp,lzp} (x-p_x)**lxp (y-p_y)**lyp (z-p_z)**lzp 03028 ! 03029 ! where p is center of the product gaussian, and lp = la_max + lb_max 03030 ! (current implementation is l**7) 03031 ! 03032 lp=la_max_local+lb_max_local 03033 ALLOCATE(coef_xyz(((lp+1)*(lp+2)*(lp+3))/6)) 03034 ALLOCATE(coef_xyt(((lp+1)*(lp+2))/2)) 03035 ALLOCATE(coef_xtt(0:lp)) 03036 ALLOCATE(alpha(0:lp,0:la_max_local,0:lb_max_local,3)) 03037 03038 ! 03039 ! compute polynomial expansion coefs -> (x-a)**lxa (x-b)**lxb -> sum_{ls} alpha(ls,lxa,lxb,1)*(x-p)**ls 03040 ! 03041 ! *** make the alpha matrix *** 03042 alpha(:,:,:,:)=0.0_dp 03043 DO iaxis=1,3 03044 DO lxa=0,la_max_local 03045 DO lxb=0,lb_max_local 03046 binomial_k_lxa=1.0_dp 03047 a=1.0_dp 03048 DO k=0,lxa 03049 binomial_l_lxb=1.0_dp 03050 b=1.0_dp 03051 DO l=0,lxb 03052 alpha(lxa-l+lxb-k,lxa,lxb,iaxis)=alpha(lxa-l+lxb-k,lxa,lxb,iaxis)+ & 03053 binomial_k_lxa*binomial_l_lxb*a*b 03054 binomial_l_lxb=binomial_l_lxb*REAL(lxb-l,dp)/REAL(l+1,dp) 03055 b=b*(rp(iaxis)-(ra(iaxis)+rab(iaxis))) 03056 ENDDO 03057 binomial_k_lxa=binomial_k_lxa*REAL(lxa-k,dp)/REAL(k+1,dp) 03058 a=a*(-ra(iaxis)+rp(iaxis)) 03059 ENDDO 03060 ENDDO 03061 ENDDO 03062 ENDDO 03063 03064 ! 03065 ! compute P_{lxp,lyp,lzp} given P_{lxa,lya,lza,lxb,lyb,lzb} and alpha(ls,lxa,lxb,1) 03066 ! use a three step procedure 03067 ! we don't store zeros, so counting is done using lxyz,lxy in order to have contiguous memory access in collocate_fast.F 03068 ! 03069 lxyz=0 03070 DO lzp=0,lp 03071 DO lyp=0,lp-lzp 03072 DO lxp=0,lp-lzp-lyp 03073 lxyz=lxyz+1 03074 coef_xyz(lxyz)=0.0_dp 03075 ENDDO 03076 ENDDO 03077 ENDDO 03078 DO lzb=0,lb_max_local 03079 DO lza=0,la_max_local 03080 lxy=0 03081 DO lyp=0,lp-lza-lzb 03082 DO lxp=0,lp-lza-lzb-lyp 03083 lxy=lxy+1 03084 coef_xyt(lxy)=0.0_dp 03085 ENDDO 03086 lxy=lxy+lza+lzb 03087 ENDDO 03088 DO lyb=0,lb_max_local-lzb 03089 DO lya=0,la_max_local-lza 03090 lxpm=(lb_max_local-lzb-lyb)+(la_max_local-lza-lya) 03091 coef_xtt(0:lxpm)=0.0_dp 03092 DO lxb=MAX(lb_min_local-lzb-lyb,0),lb_max_local-lzb-lyb 03093 DO lxa=MAX(la_min_local-lza-lya,0),la_max_local-lza-lya 03094 ico=coset(lxa,lya,lza) 03095 jco=coset(lxb,lyb,lzb) 03096 p_ele=prefactor*pab_local(o1_local+ico,o2_local+jco) 03097 DO lxp=0,lxa+lxb 03098 coef_xtt(lxp)=coef_xtt(lxp)+p_ele*alpha(lxp,lxa,lxb,1) 03099 ENDDO 03100 ENDDO 03101 ENDDO 03102 lxy=0 03103 DO lyp=0,lya+lyb 03104 DO lxp=0,lp-lza-lzb-lya-lyb 03105 lxy=lxy+1 03106 coef_xyt(lxy)=coef_xyt(lxy)+alpha(lyp,lya,lyb,2)*coef_xtt(lxp) 03107 ENDDO 03108 lxy=lxy+lza+lzb+lya+lyb-lyp 03109 ENDDO 03110 ENDDO 03111 ENDDO 03112 lxyz=0 03113 DO lzp=0,lza+lzb 03114 lxy=0 03115 DO lyp=0,lp-lza-lzb 03116 DO lxp=0,lp-lza-lzb-lyp 03117 lxy=lxy+1 ; lxyz=lxyz+1 03118 coef_xyz(lxyz)=coef_xyz(lxyz)+alpha(lzp,lza,lzb,3)*coef_xyt(lxy) 03119 ENDDO 03120 lxy=lxy+lza+lzb ; lxyz=lxyz+lza+lzb-lzp 03121 ENDDO 03122 DO lyp=lp-lza-lzb+1,lp-lzp 03123 DO lxp=0,lp-lyp-lzp 03124 lxyz=lxyz+1 03125 ENDDO 03126 ENDDO 03127 ENDDO 03128 ENDDO 03129 ENDDO 03130 03131 IF (subpatch_collocate) THEN 03132 CALL collocate_general_subpatch() 03133 ELSE 03134 IF (rsgrid%desc%orthorhombic ) THEN 03135 CALL collocate_ortho() 03136 ! CALL collocate_general() 03137 ELSE 03138 CALL collocate_general_wings() 03139 !CALL collocate_general_opt() 03140 END IF 03141 END IF 03142 03143 IF (ga_gb_function /= FUNC_AB) THEN 03144 DEALLOCATE(pab_local) 03145 ENDIF 03146 ! deallocation needed to pass around a pgi bug.. 03147 DEALLOCATE(alpha) 03148 DEALLOCATE(coef_xtt) 03149 DEALLOCATE(coef_xyt) 03150 DEALLOCATE(coef_xyz) 03151 03152 CONTAINS 03153 03154 ! 03155 ! this treats efficiently the orthogonal case 03156 ! 03157 ! ***************************************************************************** 03158 SUBROUTINE collocate_ortho() 03159 03160 03161 ! *** properties of the grid *** 03162 03163 ! notice we're in the ortho case 03164 dr(1) = rsgrid%desc%dh(1,1) 03165 dr(2) = rsgrid%desc%dh(2,2) 03166 dr(3) = rsgrid%desc%dh(3,3) 03167 03168 ! *** get the sub grid properties for the given radius *** 03169 CALL return_cube(cube_info,radius,lb_cube,ub_cube,sphere_bounds) 03170 cmax=MAXVAL(ub_cube) 03171 03172 ! *** position of the gaussian product 03173 ! 03174 ! this is the actual definition of the position on the grid 03175 ! i.e. a point rp(:) gets here grid coordinates 03176 ! MODULO(rp(:)/dr(:),ng(:))+1 03177 ! hence (0.0,0.0,0.0) in real space is rsgrid%lb on the rsgrid ((1,1,1) on grid) 03178 ! 03179 03180 03181 03182 ALLOCATE(map(-cmax:cmax,3)) 03183 CALL compute_cube_center(cubecenter,rsgrid%desc,zeta,zetb,ra,rab) 03184 roffset(:) = rp(:) - REAL(cubecenter(:),dp)*dr(:) 03185 ! *** a mapping so that the ig corresponds to the right grid point 03186 DO i=1,3 03187 IF ( rsgrid % desc % perd ( i ) == 1 ) THEN 03188 start=lb_cube(i) 03189 DO 03190 offset=MODULO(cubecenter(i)+start,ng(i))+1-start 03191 length=MIN(ub_cube(i),ng(i)-offset)-start 03192 DO ig=start,start+length 03193 map(ig,i) = ig+offset 03194 END DO 03195 IF (start+length.GE.ub_cube(i)) EXIT 03196 start=start+length+1 03197 END DO 03198 ELSE 03199 ! this takes partial grid + border regions into account 03200 offset=MODULO(cubecenter(i)+lb_cube(i)+rsgrid%desc%lb(i)-rsgrid%lb_local(i),ng(i))+1-lb_cube(i) 03201 ! check for out of bounds 03202 IF (ub_cube(i)+offset>UBOUND(grid,i).OR.lb_cube(i)+offset<LBOUND(grid,i)) THEN 03203 CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) 03204 ENDIF 03205 DO ig=lb_cube(i),ub_cube(i) 03206 map(ig,i) = ig+offset 03207 END DO 03208 END IF 03209 ENDDO 03210 ALLOCATE(pol_z(1:2,0:lp,-cmax:0)) 03211 ALLOCATE(pol_y(1:2,0:lp,-cmax:0)) 03212 ALLOCATE(pol_x(0:lp,-cmax:cmax)) 03213 03214 IF(PRESENT(ir).AND.PRESENT(rsgauge)) CALL collocate_ortho_set_to_0() 03215 03216 #include "prep.f90" 03217 03218 IF ( PRESENT ( lgrid ) ) THEN 03219 ig = lgrid%ldim * ithread_l + 1 03220 #include "call_collocate_omp.f90" 03221 ELSE 03222 #include "call_collocate.f90" 03223 END IF 03224 03225 IF(PRESENT(ir).AND.PRESENT(rsgauge)) CALL collocate_gauge_ortho() 03226 03227 ! deallocation needed to pass around a pgi bug.. 03228 DEALLOCATE(pol_z) 03229 DEALLOCATE(pol_y) 03230 DEALLOCATE(pol_x) 03231 DEALLOCATE(map) 03232 03233 END SUBROUTINE collocate_ortho 03234 03235 SUBROUTINE collocate_gauge_ortho() 03236 INTEGER :: i, igmax, igmin, j, j2, jg, 03237 jg2, jgmin, k, k2, kg, kg2, 03238 kgmin, sci 03239 REAL(KIND=dp) :: point(3,4), res(4), x, y, y2, 03240 z, z2 03241 03242 ! notice we're in the ortho case 03243 03244 dr(1) = rsgrid%desc%dh(1,1) 03245 dr(2) = rsgrid%desc%dh(2,2) 03246 dr(3) = rsgrid%desc%dh(3,3) 03247 ! 03248 sci=1 03249 kgmin=sphere_bounds(sci) 03250 sci=sci+1 03251 DO kg=kgmin,0 03252 kg2=1-kg 03253 k=map(kg,3) 03254 k2=map(kg2,3) 03255 jgmin=sphere_bounds(sci) 03256 sci=sci+1 03257 z = (REAL( kg,dp)+REAL(cubecenter(3),dp))*dr(3) 03258 z2 = (REAL(kg2,dp)+REAL(cubecenter(3),dp))*dr(3) 03259 DO jg=jgmin,0 03260 jg2=1-jg 03261 j=map(jg,2) 03262 j2=map(jg2,2) 03263 igmin=sphere_bounds(sci) 03264 sci=sci+1 03265 igmax=1-igmin 03266 y = (REAL( jg,dp)+REAL(cubecenter(2),dp))*dr(2) 03267 y2 = (REAL(jg2,dp)+REAL(cubecenter(2),dp))*dr(2) 03268 DO ig=igmin,igmax 03269 i=map(ig,1) 03270 x = (REAL(ig,dp)+REAL(cubecenter(1),dp))*dr(1) 03271 point(1,1) = x;point(2,1) = y ;point(3,1) = z 03272 point(1,2) = x;point(2,2) = y2;point(3,2) = z 03273 point(1,3) = x;point(2,3) = y ;point(3,3) = z2 03274 point(1,4) = x;point(2,4) = y2;point(3,4) = z2 03275 ! 03276 res(1) = ( point(ir,1) - rb(ir) ) - gauge(i, j, k) 03277 res(2) = ( point(ir,2) - rb(ir) ) - gauge(i,j2, k) 03278 res(3) = ( point(ir,3) - rb(ir) ) - gauge(i, j,k2) 03279 res(4) = ( point(ir,4) - rb(ir) ) - gauge(i,j2,k2) 03280 ! 03281 grid_tmp(i, j, k) = grid_tmp(i, j, k) + grid(i, j, k) * res(1) 03282 grid_tmp(i,j2, k) = grid_tmp(i,j2, k) + grid(i,j2, k) * res(2) 03283 grid_tmp(i, j,k2) = grid_tmp(i, j,k2) + grid(i, j,k2) * res(3) 03284 grid_tmp(i,j2,k2) = grid_tmp(i,j2,k2) + grid(i,j2,k2) * res(4) 03285 ENDDO 03286 ENDDO 03287 ENDDO 03288 END SUBROUTINE collocate_gauge_ortho 03289 03290 SUBROUTINE collocate_ortho_set_to_0() 03291 INTEGER :: i, igmax, igmin, j, j2, jg, 03292 jg2, jgmin, k, k2, kg, kg2, 03293 kgmin, sci 03294 03295 ! 03296 03297 dr(1) = rsgrid%desc%dh(1,1) 03298 dr(2) = rsgrid%desc%dh(2,2) 03299 dr(3) = rsgrid%desc%dh(3,3) 03300 ! 03301 sci=1 03302 kgmin=sphere_bounds(sci) 03303 sci=sci+1 03304 DO kg=kgmin,0 03305 kg2=1-kg 03306 k=map(kg,3) 03307 k2=map(kg2,3) 03308 jgmin=sphere_bounds(sci) 03309 sci=sci+1 03310 DO jg=jgmin,0 03311 jg2=1-jg 03312 j=map(jg,2) 03313 j2=map(jg2,2) 03314 igmin=sphere_bounds(sci) 03315 sci=sci+1 03316 igmax=1-igmin 03317 DO ig=igmin,igmax 03318 i=map(ig,1) 03319 grid(i, j, k) = 0.0_dp 03320 grid(i,j2, k) = 0.0_dp 03321 grid(i, j,k2) = 0.0_dp 03322 grid(i,j2,k2) = 0.0_dp 03323 ENDDO 03324 ENDDO 03325 ENDDO 03326 END SUBROUTINE collocate_ortho_set_to_0 03327 03328 ! 03329 ! this is a general 'optimized' routine to do the collocation 03330 ! 03331 ! ***************************************************************************** 03332 SUBROUTINE collocate_general_opt() 03333 03334 INTEGER :: i, i_index, il, ilx, ily, ilz, index_max(3), index_min(3), 03335 ismax, ismin, j, j_index, jl, jlx, jly, jlz, k, k_index, kl, klx, kly, 03336 klz, lpx, lpy, lpz, lx, ly, lz, offset(3) 03337 INTEGER, ALLOCATABLE, DIMENSION(:) :: grid_map 03338 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: coef_map 03339 REAL(KIND=dp) :: a, b, c, d, di, dip, dj, djp, dk, dkp, exp0i, exp1i, 03340 exp2i, gp(3), hmatgrid(3,3), pointj(3), pointk(3), res, v(3) 03341 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coef_ijk 03342 REAL(KIND=dp), ALLOCATABLE, 03343 DIMENSION(:, :, :) :: hmatgridp 03344 03345 ! 03346 ! transform P_{lxp,lyp,lzp} into a P_{lip,ljp,lkp} such that 03347 ! sum_{lxp,lyp,lzp} P_{lxp,lyp,lzp} (x-x_p)**lxp (y-y_p)**lyp (z-z_p)**lzp = 03348 ! sum_{lip,ljp,lkp} P_{lip,ljp,lkp} (i-i_p)**lip (j-j_p)**ljp (k-k_p)**lkp 03349 ! 03350 03351 ALLOCATE(coef_ijk(((lp+1)*(lp+2)*(lp+3))/6)) 03352 03353 ! aux mapping array to simplify life 03354 ALLOCATE(coef_map(0:lp,0:lp,0:lp)) 03355 coef_map=HUGE(coef_map) 03356 lxyz=0 03357 DO lzp=0,lp 03358 DO lyp=0,lp-lzp 03359 DO lxp=0,lp-lzp-lyp 03360 lxyz=lxyz+1 03361 coef_ijk(lxyz)=0.0_dp 03362 coef_map(lxp,lyp,lzp)=lxyz 03363 ENDDO 03364 ENDDO 03365 ENDDO 03366 03367 ! cell hmat in grid points 03368 hmatgrid=rsgrid%desc%dh 03369 03370 ! center in grid coords 03371 gp=MATMUL(rsgrid%desc%dh_inv,rp) 03372 cubecenter(:) = FLOOR(gp) 03373 03374 ! transform using multinomials 03375 ALLOCATE(hmatgridp(3,3,0:lp)) 03376 hmatgridp(:,:,0)=1.0_dp 03377 DO k=1,lp 03378 hmatgridp(:,:,k)=hmatgridp(:,:,k-1)*hmatgrid(:,:) 03379 ENDDO 03380 03381 lpx=lp 03382 DO klx=0,lpx 03383 DO jlx=0,lpx-klx 03384 DO ilx=0,lpx-klx-jlx 03385 lx=ilx+jlx+klx 03386 lpy=lp-lx 03387 DO kly=0,lpy 03388 DO jly=0,lpy-kly 03389 DO ily=0,lpy-kly-jly 03390 ly=ily+jly+kly 03391 lpz=lp-lx-ly 03392 DO klz=0,lpz 03393 DO jlz=0,lpz-klz 03394 DO ilz=0,lpz-klz-jlz 03395 lz=ilz+jlz+klz 03396 03397 il=ilx+ily+ilz 03398 jl=jlx+jly+jlz 03399 kl=klx+kly+klz 03400 coef_ijk(coef_map(il,jl,kl))=coef_ijk(coef_map(il,jl,kl))+ coef_xyz(coef_map(lx,ly,lz))* & 03401 hmatgridp(1,1,ilx) * hmatgridp(1,2,jlx) * hmatgridp(1,3,klx) * & 03402 hmatgridp(2,1,ily) * hmatgridp(2,2,jly) * hmatgridp(2,3,kly) * & 03403 hmatgridp(3,1,ilz) * hmatgridp(3,2,jlz) * hmatgridp(3,3,klz) * & 03404 fac(lx)*fac(ly)*fac(lz)/ & 03405 (fac(ilx)*fac(ily)*fac(ilz)*fac(jlx)*fac(jly)*fac(jlz)*fac(klx)*fac(kly)*fac(klz)) 03406 ENDDO 03407 ENDDO 03408 ENDDO 03409 ENDDO 03410 ENDDO 03411 ENDDO 03412 ENDDO 03413 ENDDO 03414 ENDDO 03415 03416 CALL return_cube_nonortho(cube_info,radius,index_min,index_max,rp) 03417 03418 offset(:)=MODULO(index_min(:)+rsgrid%desc%lb(:)-rsgrid%lb_local(:),ng(:))+1 03419 03420 ALLOCATE(grid_map(index_min(1):index_max(1))) 03421 DO i=index_min(1),index_max(1) 03422 grid_map(i)=MODULO(i,ng(1))+1 03423 IF (rsgrid % desc % perd ( 1 )==1) THEN 03424 grid_map(i)=MODULO(i,ng(1))+1 03425 ELSE 03426 grid_map(i)=i-index_min(1)+offset(1) 03427 ENDIF 03428 ENDDO 03429 03430 ! go over the grid, but cycle if the point is not within the radius 03431 DO k=index_min(3),index_max(3) 03432 dk=k-gp(3) 03433 pointk=hmatgrid(:,3)*dk 03434 03435 IF (rsgrid % desc % perd ( 3 )==1) THEN 03436 k_index=MODULO(k,ng(3))+1 03437 ELSE 03438 k_index=k-index_min(3)+offset(3) 03439 ENDIF 03440 03441 coef_xyt=0.0_dp 03442 lxyz = 0 03443 dkp=1.0_dp 03444 DO kl=0,lp 03445 lxy=0 03446 DO jl=0,lp-kl 03447 DO il=0,lp-kl-jl 03448 lxyz=lxyz+1 ; lxy=lxy+1 03449 coef_xyt(lxy)=coef_xyt(lxy)+coef_ijk(lxyz)*dkp 03450 ENDDO 03451 lxy=lxy+kl 03452 ENDDO 03453 dkp=dkp*dk 03454 ENDDO 03455 03456 DO j=index_min(2),index_max(2) 03457 dj=j-gp(2) 03458 pointj=pointk+hmatgrid(:,2)*dj 03459 IF (rsgrid % desc % perd ( 2 )==1) THEN 03460 j_index=MODULO(j,ng(2))+1 03461 ELSE 03462 j_index=j-index_min(2)+offset(2) 03463 ENDIF 03464 03465 coef_xtt=0.0_dp 03466 lxy=0 03467 djp=1.0_dp 03468 DO jl=0,lp 03469 DO il=0,lp-jl 03470 lxy=lxy+1 03471 coef_xtt(il)=coef_xtt(il)+coef_xyt(lxy)*djp 03472 ENDDO 03473 djp=djp*dj 03474 ENDDO 03475 03476 ! find bounds for the inner loop 03477 ! based on a quadratic equation in i 03478 ! a*i**2+b*i+c=radius**2 03479 v=pointj-gp(1)*hmatgrid(:,1) 03480 a=DOT_PRODUCT(hmatgrid(:,1),hmatgrid(:,1)) 03481 b=2*DOT_PRODUCT(v,hmatgrid(:,1)) 03482 c=DOT_PRODUCT(v,v) 03483 d=b*b-4*a*(c-radius**2) 03484 03485 IF (d<0) THEN 03486 CYCLE 03487 ELSE 03488 d=SQRT(d) 03489 ismin=CEILING((-b-d)/(2*a)) 03490 ismax=FLOOR((-b+d)/(2*a)) 03491 ENDIF 03492 ! prepare for computing -zetp*rsq 03493 a=-zetp*a 03494 b=-zetp*b 03495 c=-zetp*c 03496 i=ismin-1 03497 03498 ! the recursion relation might have to be done 03499 ! from the center of the gaussian (in both directions) 03500 ! instead as the current implementation from an edge 03501 exp2i=EXP((a*i+b)*i+c) 03502 exp1i=EXP(2*a*i+a+b) 03503 exp0i=EXP(2*a) 03504 03505 DO i=ismin,ismax 03506 di=i-gp(1) 03507 03508 ! polynomial terms 03509 res=0.0_dp 03510 dip=1.0_dp 03511 DO il=0,lp 03512 res=res+coef_xtt(il)*dip 03513 dip=dip*di 03514 ENDDO 03515 03516 ! the exponential recursion 03517 exp2i=exp2i*exp1i 03518 exp1i=exp1i*exp0i 03519 res=res*exp2i 03520 03521 i_index=grid_map(i) 03522 IF ( PRESENT ( lgrid ) ) THEN 03523 ig = lgrid%ldim * ithread_l + (k_index-1) * ng(2) * ng(1) + (j_index-1) * ng(1) + (i_index-1) + 1 03524 lgrid%r(ig)=lgrid%r(ig)+res 03525 ELSE 03526 grid(i_index,j_index,k_index)=grid(i_index,j_index,k_index)+res 03527 ENDIF 03528 ENDDO 03529 ENDDO 03530 ENDDO 03531 !t2=nanotime_ia32() 03532 !write(6,*) t2-t1 03533 ! deallocation needed to pass around a pgi bug.. 03534 DEALLOCATE(coef_ijk) 03535 DEALLOCATE(coef_map) 03536 DEALLOCATE(hmatgridp) 03537 DEALLOCATE(grid_map) 03538 03539 END SUBROUTINE collocate_general_opt 03540 03541 SUBROUTINE collocate_general_subpatch 03542 INTEGER, DIMENSION(2, 3) :: local_b 03543 INTEGER, DIMENSION(3) :: local_s, periodic 03544 REAL(dp), 03545 DIMENSION((lp+1)*(lp+2)*(lp+3)/6) :: poly_d3 03546 03547 periodic=1 ! cell%perd 03548 CALL poly_cp2k2d3(coef_xyz,lp,poly_d3,error) 03549 local_b(1,:)=rsgrid%lb_real-rsgrid%desc%lb 03550 local_b(2,:)=rsgrid%ub_real-rsgrid%desc%lb 03551 local_s=rsgrid%lb_real-rsgrid%lb_local 03552 IF (BTEST(subpatch_pattern,0)) local_b(1,1)=local_b(1,1)-rsgrid%desc%border 03553 IF (BTEST(subpatch_pattern,1)) local_b(2,1)=local_b(2,1)+rsgrid%desc%border 03554 IF (BTEST(subpatch_pattern,2)) local_b(1,2)=local_b(1,2)-rsgrid%desc%border 03555 IF (BTEST(subpatch_pattern,3)) local_b(2,2)=local_b(2,2)+rsgrid%desc%border 03556 IF (BTEST(subpatch_pattern,4)) local_b(1,3)=local_b(1,3)-rsgrid%desc%border 03557 IF (BTEST(subpatch_pattern,5)) local_b(2,3)=local_b(2,3)+rsgrid%desc%border 03558 IF (BTEST(subpatch_pattern,0)) local_s(1)=local_s(1)-rsgrid%desc%border 03559 IF (BTEST(subpatch_pattern,2)) local_s(2)=local_s(2)-rsgrid%desc%border 03560 IF (BTEST(subpatch_pattern,4)) local_s(3)=local_s(3)-rsgrid%desc%border 03561 IF ( PRESENT ( lgrid ) ) THEN 03562 CALL collocGauss(h=cell%hmat,h_inv=cell%h_inv,& 03563 grid=grid,poly=poly_d3,alphai=zetp,posi=rp,max_r2=radius*radius,& 03564 periodic=periodic,gdim=ng,local_bounds=local_b,local_shift=local_s,& 03565 lgrid=lgrid,error=error) 03566 ELSE 03567 CALL collocGauss(h=cell%hmat,h_inv=cell%h_inv,& 03568 grid=grid,poly=poly_d3,alphai=zetp,posi=rp,max_r2=radius*radius,& 03569 periodic=periodic,gdim=ng,local_bounds=local_b,local_shift=local_s,& 03570 error=error) 03571 END IF 03572 ! defaults: local_shift=(/0,0,0/),poly_shift=(/0.0_dp,0.0_dp,0.0_dp/),scale=1.0_dp, 03573 03574 END SUBROUTINE 03575 03576 SUBROUTINE collocate_general_wings 03577 INTEGER, DIMENSION(2, 3) :: local_b 03578 INTEGER, DIMENSION(3) :: periodic 03579 REAL(dp), 03580 DIMENSION((lp+1)*(lp+2)*(lp+3)/6) :: poly_d3 03581 REAL(dp), DIMENSION(3) :: local_shift, rShifted 03582 03583 periodic=1 ! cell%perd 03584 CALL poly_cp2k2d3(coef_xyz,lp,poly_d3,error) 03585 local_b(1,:)=0 03586 local_b(2,:)=MIN(rsgrid%desc%npts-1,rsgrid%ub_local-rsgrid%lb_local) 03587 local_shift=REAL(rsgrid%desc%lb-rsgrid%lb_local,dp)/REAL(rsgrid%desc%npts,dp) 03588 rShifted(1)=rp(1)+cell%hmat(1,1)*local_shift(1)& 03589 +cell%hmat(1,2)*local_shift(2)& 03590 +cell%hmat(1,3)*local_shift(3) 03591 rShifted(2)=rp(2)+cell%hmat(2,1)*local_shift(1)& 03592 +cell%hmat(2,2)*local_shift(2)& 03593 +cell%hmat(2,3)*local_shift(3) 03594 rShifted(3)=rp(3)+cell%hmat(3,1)*local_shift(1)& 03595 +cell%hmat(3,2)*local_shift(2)& 03596 +cell%hmat(3,3)*local_shift(3) 03597 IF ( PRESENT ( lgrid ) ) THEN 03598 CALL collocGauss(h=cell%hmat,h_inv=cell%h_inv,& 03599 grid=grid,poly=poly_d3,alphai=zetp,posi=rShifted,max_r2=radius*radius,& 03600 periodic=periodic,gdim=ng,local_bounds=local_b,& 03601 lgrid=lgrid,error=error) 03602 ELSE 03603 CALL collocGauss(h=cell%hmat,h_inv=cell%h_inv,& 03604 grid=grid,poly=poly_d3,alphai=zetp,posi=rShifted,max_r2=radius*radius,& 03605 periodic=periodic,gdim=ng,local_bounds=local_b,& 03606 error=error) 03607 END IF 03608 ! defaults: local_shift=(/0,0,0/),poly_shift=(/0.0_dp,0.0_dp,0.0_dp/),scale=1.0_dp, 03609 03610 END SUBROUTINE 03611 03612 ! 03613 ! this is a general 'reference' routine to do the collocation 03614 ! 03615 ! ***************************************************************************** 03616 SUBROUTINE collocate_general() 03617 03618 INTEGER :: i, index_max(3), 03619 index_min(3), ipoint(3), j, k 03620 REAL(KIND=dp) :: point(3) 03621 03622 ! still hard-wired (see MODULO) 03623 03624 CPPostcondition(ALL(rsgrid % desc % perd==1),cp_failure_level,routineP,error,failure) 03625 03626 CALL return_cube_nonortho(cube_info,radius,index_min,index_max,rp) 03627 03628 ! go over the grid, but cycle if the point is not within the radius 03629 DO k=index_min(3),index_max(3) 03630 DO j=index_min(2),index_max(2) 03631 DO i=index_min(1),index_max(1) 03632 ! point in real space 03633 point=MATMUL(cell%hmat,REAL((/i,j,k/),KIND=dp)/ng) 03634 ! skip if outside of the sphere 03635 IF (SUM((point-rp)**2)>radius**2) CYCLE 03636 ! point on the grid (including pbc) 03637 ipoint=MODULO((/i,j,k/),ng)+1 03638 ! add to grid 03639 IF ( PRESENT ( lgrid ) ) THEN 03640 ig = lgrid%ldim * ithread_l + ipoint(3) * ng(2) * ng(1) + ipoint(2) * ng(1) + ipoint(1) + 1 03641 lgrid%r(ig)=lgrid%r(ig)+primitive_value(point) 03642 ELSE 03643 grid(ipoint(1),ipoint(2),ipoint(3))=grid(ipoint(1),ipoint(2),ipoint(3))+primitive_value(point) 03644 ENDIF 03645 ENDDO 03646 ENDDO 03647 ENDDO 03648 03649 END SUBROUTINE collocate_general 03650 03651 ! ***************************************************************************** 03652 FUNCTION primitive_value(point) RESULT(res) 03653 REAL(KIND=dp) :: point(3), res 03654 03655 REAL(KIND=dp) :: dra(3), drap(3), drb(3), 03656 drbp(3), myexp 03657 03658 res=0.0_dp 03659 myexp=EXP(-zetp*SUM((point-rp)**2))*prefactor 03660 dra=point-ra 03661 drb=point-rb 03662 drap(1)=1.0_dp 03663 DO lxa=0,la_max_local 03664 drbp(1)=1.0_dp 03665 DO lxb=0,lb_max_local 03666 drap(2)=1.0_dp 03667 DO lya=0,la_max_local-lxa 03668 drbp(2)=1.0_dp 03669 DO lyb=0,lb_max_local-lxb 03670 drap(3)=1.0_dp 03671 DO lza=1,MAX(la_min_local-lxa-lya,0) 03672 drap(3)=drap(3)*dra(3) 03673 ENDDO 03674 DO lza=MAX(la_min_local-lxa-lya,0),la_max_local-lxa-lya 03675 drbp(3)=1.0_dp 03676 DO lzb=1,MAX(lb_min_local-lxb-lyb,0) 03677 drbp(3)=drbp(3)*drb(3) 03678 ENDDO 03679 DO lzb=MAX(lb_min_local-lxb-lyb,0),lb_max_local-lxb-lyb 03680 ico=coset(lxa,lya,lza) 03681 jco=coset(lxb,lyb,lzb) 03682 res=res+pab_local(ico+o1_local,jco+o2_local)*myexp*PRODUCT(drap)*PRODUCT(drbp) 03683 drbp(3)=drbp(3)*drb(3) 03684 ENDDO 03685 drap(3)=drap(3)*dra(3) 03686 ENDDO 03687 drbp(2)=drbp(2)*drb(2) 03688 ENDDO 03689 drap(2)=drap(2)*dra(2) 03690 ENDDO 03691 drbp(1)=drbp(1)*drb(1) 03692 ENDDO 03693 drap(1)=drap(1)*dra(1) 03694 ENDDO 03695 03696 END FUNCTION primitive_value 03697 03698 END SUBROUTINE collocate_pgf_product_rspace 03699 03700 ! ***************************************************************************** 03703 SUBROUTINE collocate_pgf_product_gspace(la_max,zeta,la_min,& 03704 lb_max,zetb,lb_min,& 03705 ra,rab,rab2,scale,pab,na,nb,& 03706 eps_rho_gspace,gsq_max,pw) 03707 03708 ! NOTE: this routine is much slower than collocate_pgf_product_rspace 03709 03710 INTEGER, INTENT(IN) :: la_max 03711 REAL(dp), INTENT(IN) :: zeta 03712 INTEGER, INTENT(IN) :: la_min, lb_max 03713 REAL(dp), INTENT(IN) :: zetb 03714 INTEGER, INTENT(IN) :: lb_min 03715 REAL(dp), DIMENSION(3), INTENT(IN) :: ra, rab 03716 REAL(dp), INTENT(IN) :: rab2, scale 03717 REAL(dp), DIMENSION(:, :), POINTER :: pab 03718 INTEGER, INTENT(IN) :: na, nb 03719 REAL(dp), INTENT(IN) :: eps_rho_gspace, gsq_max 03720 TYPE(pw_type), POINTER :: pw 03721 03722 CHARACTER(LEN=*), PARAMETER :: routineN = 'collocate_pgf_product_gspace', 03723 routineP = moduleN//':'//routineN 03724 03725 COMPLEX(dp), DIMENSION(3) :: phasefactor 03726 COMPLEX(dp), DIMENSION(:), POINTER :: rag, rbg 03727 COMPLEX(dp), DIMENSION(:, :, :, :), 03728 POINTER :: cubeaxis 03729 INTEGER :: ax, ay, az, bx, by, bz, handle, i, ico, ig, ig2, jco, jg, kg, 03730 la, lb, lb_cube_min, lb_grid, stat, ub_cube_max, ub_grid 03731 INTEGER, DIMENSION(3) :: lb_cube, ub_cube 03732 REAL(dp) :: f, fa, fb, pij, prefactor, 03733 rzetp, twozetp, zetp 03734 REAL(dp), DIMENSION(3) :: dg, expfactor, fap, fbp, rap, 03735 rbp, rp 03736 REAL(dp), DIMENSION(:), POINTER :: g 03737 03738 CALL timeset(routineN,handle) 03739 03740 dg(:) = twopi/(pw%pw_grid%npts(:)*pw%pw_grid%dr(:)) 03741 03742 zetp = zeta + zetb 03743 rzetp = 1.0_dp/zetp 03744 f = zetb*rzetp 03745 rap(:) = f*rab(:) 03746 rbp(:) = rap(:) - rab(:) 03747 rp(:) = ra(:) + rap(:) 03748 twozetp = 2.0_dp*zetp 03749 fap(:) = twozetp*rap(:) 03750 fbp(:) = twozetp*rbp(:) 03751 03752 prefactor = scale*SQRT((pi*rzetp)**3)*EXP(-zeta*f*rab2) 03753 phasefactor(:) = EXP(CMPLX(0.0_dp,-rp(:)*dg(:),KIND=dp)) 03754 expfactor(:) = EXP(-0.25*rzetp*dg(:)*dg(:)) 03755 03756 lb_cube(:) = pw%pw_grid%bounds(1,:) 03757 ub_cube(:) = pw%pw_grid%bounds(2,:) 03758 03759 lb_cube_min = MINVAL(lb_cube(:)) 03760 ub_cube_max = MAXVAL(ub_cube(:)) 03761 03762 NULLIFY (cubeaxis,g,rag,rbg) 03763 03764 CALL reallocate(cubeaxis,lb_cube_min,ub_cube_max,1,3,0,la_max,0,lb_max) 03765 CALL reallocate(g,lb_cube_min,ub_cube_max) 03766 CALL reallocate(rag,lb_cube_min,ub_cube_max) 03767 CALL reallocate(rbg,lb_cube_min,ub_cube_max) 03768 03769 lb_grid = LBOUND(pw%cc,1) 03770 ub_grid = UBOUND(pw%cc,1) 03771 03772 DO i=1,3 03773 03774 DO ig=lb_cube(i),ub_cube(i) 03775 ig2 = ig*ig 03776 cubeaxis(ig,i,0,0) = expfactor(i)**ig2*phasefactor(i)**ig 03777 END DO 03778 03779 IF (la_max > 0) THEN 03780 DO ig=lb_cube(i),ub_cube(i) 03781 g(ig) = REAL(ig,dp)*dg(i) 03782 rag(ig) = CMPLX(fap(i),-g(ig),KIND=dp) 03783 cubeaxis(ig,i,1,0) = rag(ig)*cubeaxis(ig,i,0,0) 03784 END DO 03785 DO la=2,la_max 03786 fa = REAL(la-1,dp)*twozetp 03787 DO ig=lb_cube(i),ub_cube(i) 03788 cubeaxis(ig,i,la,0) = rag(ig)*cubeaxis(ig,i,la-1,0) +& 03789 fa*cubeaxis(ig,i,la-2,0) 03790 END DO 03791 END DO 03792 IF (lb_max > 0) THEN 03793 fa = twozetp 03794 DO ig=lb_cube(i),ub_cube(i) 03795 rbg(ig) = CMPLX(fbp(i),-g(ig),KIND=dp) 03796 cubeaxis(ig,i,0,1) = rbg(ig)*cubeaxis(ig,i,0,0) 03797 cubeaxis(ig,i,1,1) = rbg(ig)*cubeaxis(ig,i,1,0) +& 03798 fa*cubeaxis(ig,i,0,0) 03799 END DO 03800 DO lb=2,lb_max 03801 fb = REAL(lb-1,dp)*twozetp 03802 DO ig=lb_cube(i),ub_cube(i) 03803 cubeaxis(ig,i,0,lb) = rbg(ig)*cubeaxis(ig,i,0,lb-1) +& 03804 fb*cubeaxis(ig,i,0,lb-2) 03805 cubeaxis(ig,i,1,lb) = rbg(ig)*cubeaxis(ig,i,1,lb-1) +& 03806 fb*cubeaxis(ig,i,1,lb-2) +& 03807 fa*cubeaxis(ig,i,0,lb-1) 03808 END DO 03809 END DO 03810 DO la=2,la_max 03811 fa = REAL(la,dp)*twozetp 03812 DO ig=lb_cube(i),ub_cube(i) 03813 cubeaxis(ig,i,la,1) = rbg(ig)*cubeaxis(ig,i,la,0) +& 03814 fa*cubeaxis(ig,i,la-1,0) 03815 END DO 03816 DO lb=2,lb_max 03817 fb = REAL(lb-1,dp)*twozetp 03818 DO ig=lb_cube(i),ub_cube(i) 03819 cubeaxis(ig,i,la,lb) = rbg(ig)*cubeaxis(ig,i,la,lb-1) +& 03820 fb*cubeaxis(ig,i,la,lb-2) +& 03821 fa*cubeaxis(ig,i,la-1,lb-1) 03822 END DO 03823 END DO 03824 END DO 03825 END IF 03826 ELSE 03827 IF (lb_max > 0) THEN 03828 DO ig=lb_cube(i),ub_cube(i) 03829 g(ig) = REAL(ig,dp)*dg(i) 03830 rbg(ig) = CMPLX(fbp(i),-g(ig),KIND=dp) 03831 cubeaxis(ig,i,0,1) = rbg(ig)*cubeaxis(ig,i,0,0) 03832 END DO 03833 DO lb=2,lb_max 03834 fb = REAL(lb-1,dp)*twozetp 03835 DO ig=lb_cube(i),ub_cube(i) 03836 cubeaxis(ig,i,0,lb) = rbg(ig)*cubeaxis(ig,i,0,lb-1) +& 03837 fb*cubeaxis(ig,i,0,lb-2) 03838 END DO 03839 END DO 03840 END IF 03841 END IF 03842 03843 END DO 03844 03845 DO la=0,la_max 03846 DO lb=0,lb_max 03847 IF (la + lb == 0) CYCLE 03848 fa = (1.0_dp/twozetp)**(la + lb) 03849 DO i=1,3 03850 DO ig=lb_cube(i),ub_cube(i) 03851 cubeaxis(ig,i,la,lb) = fa*cubeaxis(ig,i,la,lb) 03852 END DO 03853 END DO 03854 END DO 03855 END DO 03856 03857 ! Add the current primitive Gaussian function product to grid 03858 03859 DO ico=ncoset(la_min-1)+1,ncoset(la_max) 03860 03861 ax = indco(1,ico) 03862 ay = indco(2,ico) 03863 az = indco(3,ico) 03864 03865 DO jco=ncoset(lb_min-1)+1,ncoset(lb_max) 03866 03867 pij = prefactor*pab(na+ico,nb+jco) 03868 03869 IF (ABS(pij) < eps_rho_gspace) CYCLE 03870 03871 bx = indco(1,jco) 03872 by = indco(2,jco) 03873 bz = indco(3,jco) 03874 03875 DO i=lb_grid,ub_grid 03876 IF (pw%pw_grid%gsq(i) > gsq_max) CYCLE 03877 ig = pw%pw_grid%g_hat(1,i) 03878 jg = pw%pw_grid%g_hat(2,i) 03879 kg = pw%pw_grid%g_hat(3,i) 03880 pw%cc(i) = pw%cc(i) + pij*cubeaxis(ig,1,ax,bx)*& 03881 cubeaxis(jg,2,ay,by)*& 03882 cubeaxis(kg,3,az,bz) 03883 END DO 03884 03885 END DO 03886 03887 END DO 03888 03889 DEALLOCATE (cubeaxis,STAT=stat) 03890 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"cubeaxis") 03891 DEALLOCATE (g,STAT=stat) 03892 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"g") 03893 DEALLOCATE (rag,STAT=stat) 03894 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"rag") 03895 DEALLOCATE (rbg,STAT=stat) 03896 IF (stat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"rbg") 03897 03898 CALL timestop(handle) 03899 03900 END SUBROUTINE collocate_pgf_product_gspace 03901 03902 ! ***************************************************************************** 03910 SUBROUTINE density_rs2pw(pw_env,rs_rho,rho,rho_gspace,interp_section,error) 03911 03912 TYPE(pw_env_type), POINTER :: pw_env 03913 TYPE(realspace_grid_p_type), 03914 DIMENSION(:), POINTER :: rs_rho 03915 TYPE(pw_p_type), INTENT(INOUT) :: rho, rho_gspace 03916 TYPE(section_vals_type), OPTIONAL, 03917 POINTER :: interp_section 03918 TYPE(cp_error_type), INTENT(inout) :: error 03919 03920 CHARACTER(LEN=*), PARAMETER :: routineN = 'density_rs2pw', 03921 routineP = moduleN//':'//routineN 03922 03923 INTEGER :: handle, igrid_level, 03924 interp_kind 03925 LOGICAL :: failure 03926 TYPE(gridlevel_info_type), POINTER :: gridlevel_info 03927 TYPE(pw_p_type), DIMENSION(:), POINTER :: mgrid_gspace, mgrid_rspace 03928 TYPE(pw_pool_p_type), DIMENSION(:), 03929 POINTER :: pw_pools 03930 TYPE(realspace_grid_desc_p_type), 03931 DIMENSION(:), POINTER :: rs_descs 03932 03933 CALL timeset(routineN,handle) 03934 failure = .FALSE. 03935 NULLIFY(gridlevel_info, mgrid_gspace, mgrid_rspace, rs_descs, pw_pools) 03936 CPPrecondition(ASSOCIATED(pw_env),cp_failure_level,routineP,error,failure) 03937 CALL pw_env_get(pw_env, rs_descs=rs_descs, pw_pools=pw_pools, error=error) 03938 03939 gridlevel_info=>pw_env%gridlevel_info 03940 IF(PRESENT(interp_section)) THEN 03941 CALL section_vals_val_get(interp_section,"KIND",i_val=interp_kind,error=error) 03942 ELSE 03943 interp_kind = pw_interp 03944 END IF 03945 03946 CALL pw_pools_create_pws(pw_pools,mgrid_rspace,& 03947 use_data = REALDATA3D,& 03948 in_space = REALSPACE, error=error) 03949 03950 CALL pw_pools_create_pws(pw_pools,mgrid_gspace,& 03951 use_data = COMPLEXDATA1D,& 03952 in_space = RECIPROCALSPACE, error=error) 03953 03954 IF (gridlevel_info%ngrid_levels==1) THEN 03955 CALL rs_pw_transfer(rs_rho(1)%rs_grid,rho%pw,rs2pw,error=error) 03956 CALL rs_grid_release(rs_rho(1)%rs_grid, error=error) 03957 CALL pw_transfer(rho%pw,rho_gspace%pw,error=error) 03958 IF (rho%pw%pw_grid%spherical) THEN ! rho_gspace = rho 03959 CALL pw_transfer(rho_gspace%pw,rho%pw,error=error) 03960 ENDIF 03961 ELSE 03962 DO igrid_level=1,gridlevel_info%ngrid_levels 03963 CALL rs_pw_transfer(rs_rho(igrid_level)%rs_grid,& 03964 mgrid_rspace(igrid_level)%pw,rs2pw,error=error) 03965 CALL rs_grid_release(rs_rho(igrid_level)%rs_grid, error=error) 03966 ENDDO 03967 03968 ! we want both rho and rho_gspace, the latter for Hartree and co-workers. 03969 SELECT CASE(interp_kind) 03970 CASE(pw_interp) 03971 CALL pw_zero(rho_gspace%pw,error=error) 03972 DO igrid_level=1,gridlevel_info%ngrid_levels 03973 CALL pw_transfer(mgrid_rspace(igrid_level)%pw,& 03974 mgrid_gspace(igrid_level)%pw,error=error) 03975 CALL pw_axpy(mgrid_gspace(igrid_level)%pw,rho_gspace%pw,& 03976 error=error) 03977 END DO 03978 CALL pw_transfer(rho_gspace%pw,rho%pw,error=error) 03979 CASE(spline3_pbc_interp) 03980 DO igrid_level=gridlevel_info%ngrid_levels,2,-1 03981 CALL pw_prolongate_s3(mgrid_rspace(igrid_level)%pw,& 03982 mgrid_rspace(igrid_level-1)%pw,pw_pools(igrid_level)%pool,& 03983 interp_section,error=error) 03984 END DO 03985 CALL pw_copy(mgrid_rspace(1)%pw,rho%pw,error=error) 03986 CALL pw_transfer(rho%pw,rho_gspace%pw,error=error) 03987 CASE default 03988 CALL cp_unimplemented_error(routineN,"interpolator "//& 03989 cp_to_string(interp_kind),error=error) 03990 END SELECT 03991 END IF 03992 03993 ! *** give back the pw multi-grids 03994 CALL pw_pools_give_back_pws(pw_pools,mgrid_gspace,error=error) 03995 CALL pw_pools_give_back_pws(pw_pools,mgrid_rspace,error=error) 03996 CALL timestop(handle) 03997 03998 END SUBROUTINE density_rs2pw 03999 04000 ! ***************************************************************************** 04007 SUBROUTINE density_rs2pw_basic(pw_env,rs_rho,rho,rho_gspace,interp_section,error) 04008 04009 TYPE(pw_env_type), POINTER :: pw_env 04010 TYPE(realspace_grid_p_type), 04011 DIMENSION(:), POINTER :: rs_rho 04012 TYPE(pw_p_type), INTENT(INOUT) :: rho, rho_gspace 04013 TYPE(section_vals_type), OPTIONAL, 04014 POINTER :: interp_section 04015 TYPE(cp_error_type), INTENT(inout) :: error 04016 04017 CHARACTER(LEN=*), PARAMETER :: routineN = 'density_rs2pw_basic', 04018 routineP = moduleN//':'//routineN 04019 04020 INTEGER :: handle, igrid_level, 04021 interp_kind 04022 LOGICAL :: failure 04023 TYPE(gridlevel_info_type), POINTER :: gridlevel_info 04024 TYPE(pw_p_type), DIMENSION(:), POINTER :: mgrid_gspace, mgrid_rspace 04025 TYPE(pw_pool_p_type), DIMENSION(:), 04026 POINTER :: pw_pools 04027 TYPE(realspace_grid_desc_p_type), 04028 DIMENSION(:), POINTER :: rs_descs 04029 04030 CALL timeset(routineN,handle) 04031 failure = .FALSE. 04032 NULLIFY(gridlevel_info, mgrid_gspace, mgrid_rspace, rs_descs, pw_pools) 04033 CPPrecondition(ASSOCIATED(pw_env),cp_failure_level,routineP,error,failure) 04034 CALL pw_env_get(pw_env, rs_descs=rs_descs, pw_pools=pw_pools, error=error) 04035 04036 04037 gridlevel_info=>pw_env%gridlevel_info 04038 IF(PRESENT(interp_section)) THEN 04039 CALL section_vals_val_get(interp_section,"KIND",i_val=interp_kind,error=error) 04040 ELSE 04041 interp_kind = pw_interp 04042 END IF 04043 04044 CALL pw_pools_create_pws(pw_pools,mgrid_rspace,& 04045 use_data = REALDATA3D,& 04046 in_space = REALSPACE, error=error) 04047 04048 CALL pw_pools_create_pws(pw_pools,mgrid_gspace,& 04049 use_data = COMPLEXDATA1D,& 04050 in_space = RECIPROCALSPACE, error=error) 04051 04052 IF (gridlevel_info%ngrid_levels==1) THEN 04053 CALL rs_pw_transfer(rs_rho(1)%rs_grid,rho%pw,rs2pw,error=error) 04054 CALL pw_transfer(rho%pw,rho_gspace%pw,error=error) 04055 ELSE 04056 DO igrid_level=1,gridlevel_info%ngrid_levels 04057 CALL rs_pw_transfer(rs_rho(igrid_level)%rs_grid,& 04058 mgrid_rspace(igrid_level)%pw,rs2pw,error=error) 04059 ENDDO 04060 04061 ! we want both rho and rho_gspace, the latter for Hartree and co-workers. 04062 SELECT CASE(interp_kind) 04063 CASE(pw_interp) 04064 DO igrid_level=1,gridlevel_info%ngrid_levels 04065 CALL pw_transfer(mgrid_rspace(igrid_level)%pw,& 04066 mgrid_gspace(igrid_level)%pw,error=error) 04067 IF (igrid_level/=1) THEN 04068 CALL pw_axpy(mgrid_gspace(igrid_level)%pw,mgrid_gspace(1)%pw,& 04069 error=error) 04070 END IF 04071 END DO 04072 CALL pw_transfer(mgrid_gspace(1)%pw,rho%pw,error=error) 04073 CALL pw_transfer(mgrid_rspace(1)%pw,rho_gspace%pw,error=error) 04074 CASE(spline3_pbc_interp) 04075 DO igrid_level=gridlevel_info%ngrid_levels,2,-1 04076 CALL pw_prolongate_s3(mgrid_rspace(igrid_level)%pw,& 04077 mgrid_rspace(igrid_level-1)%pw,pw_pools(igrid_level)%pool,& 04078 interp_section,error=error) 04079 END DO 04080 CALL pw_copy(mgrid_rspace(1)%pw,rho%pw,error=error) 04081 CALL pw_transfer(rho%pw,rho_gspace%pw,error=error) 04082 CASE default 04083 CALL cp_unimplemented_error(routineN,"interpolator "//& 04084 cp_to_string(interp_kind),error=error) 04085 END SELECT 04086 END IF 04087 04088 ! *** give back the pw multi-grids 04089 CALL pw_pools_give_back_pws(pw_pools,mgrid_gspace,error=error) 04090 CALL pw_pools_give_back_pws(pw_pools,mgrid_rspace,error=error) 04091 CALL timestop(handle) 04092 04093 END SUBROUTINE density_rs2pw_basic 04094 04095 END MODULE qs_collocate_density
1.7.3