|
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 ! ***************************************************************************** 00014 MODULE xray_diffraction 00015 00016 USE atomic_kind_types, ONLY: atomic_kind_type,& 00017 get_atomic_kind 00018 USE basis_set_types, ONLY: get_gto_basis_set,& 00019 gto_basis_set_type 00020 USE bibliography, ONLY: Krack2002,& 00021 cite_reference 00022 USE cell_types, ONLY: cell_type,& 00023 pbc 00024 USE cp_control_types, ONLY: dft_control_type 00025 USE cp_para_types, ONLY: cp_para_env_type 00026 USE kinds, ONLY: dp,& 00027 int_8 00028 USE mathconstants, ONLY: twopi 00029 USE memory_utilities, ONLY: reallocate 00030 USE message_passing, ONLY: mp_gather,& 00031 mp_gatherv 00032 USE orbital_pointers, ONLY: nco,& 00033 ncoset,& 00034 nso,& 00035 nsoset 00036 USE orbital_transformation_matrices, ONLY: orbtramat 00037 USE particle_types, ONLY: particle_type 00038 USE paw_proj_set_types, ONLY: get_paw_proj_set,& 00039 paw_proj_set_type 00040 USE physcon, ONLY: angstrom 00041 USE pw_env_types, ONLY: pw_env_get,& 00042 pw_env_type 00043 USE pw_grids, ONLY: get_pw_grid_info 00044 USE pw_methods, ONLY: pw_axpy,& 00045 pw_integrate_function,& 00046 pw_scale,& 00047 pw_transfer,& 00048 pw_zero 00049 USE pw_pool_types, ONLY: pw_pool_create_pw,& 00050 pw_pool_give_back_pw,& 00051 pw_pool_type 00052 USE pw_types, ONLY: COMPLEXDATA1D,& 00053 RECIPROCALSPACE,& 00054 pw_p_type 00055 USE qs_collocate_density, ONLY: collocate_pgf_product_gspace 00056 USE qs_environment_types, ONLY: get_qs_env,& 00057 qs_environment_type 00058 USE qs_rho_atom_types, ONLY: get_rho_atom,& 00059 rho_atom_coeff,& 00060 rho_atom_type 00061 USE qs_rho_types, ONLY: qs_rho_type 00062 USE timings, ONLY: timeset,& 00063 timestop 00064 USE util, ONLY: sort 00065 #include "cp_common_uses.h" 00066 00067 IMPLICIT NONE 00068 00069 PRIVATE 00070 00071 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xray_diffraction' 00072 00073 PUBLIC :: calculate_rhotot_elec_gspace,& 00074 xray_diffraction_spectrum 00075 00076 CONTAINS 00077 00078 ! ***************************************************************************** 00084 SUBROUTINE xray_diffraction_spectrum(qs_env,unit_number,q_max,error) 00085 00086 TYPE(qs_environment_type), POINTER :: qs_env 00087 INTEGER, INTENT(IN) :: unit_number 00088 REAL(KIND=dp), INTENT(IN) :: q_max 00089 TYPE(cp_error_type), INTENT(INOUT) :: error 00090 00091 CHARACTER(LEN=*), PARAMETER :: routineN = 'xray_diffraction_spectrum', 00092 routineP = moduleN//':'//routineN 00093 INTEGER, PARAMETER :: nblock = 100 00094 00095 INTEGER :: group, handle, i, ig, 00096 ig_shell, ipe, ishell, istat, 00097 jg, ng, npe, nshell, 00098 nshell_gather, source 00099 INTEGER(KIND=int_8) :: ngpts 00100 INTEGER, DIMENSION(3) :: npts 00101 INTEGER, DIMENSION(:), POINTER :: aux_index, ng_shell, 00102 ng_shell_gather, nshell_pe, 00103 offset_pe 00104 LOGICAL :: failure 00105 REAL(KIND=dp) :: cutoff, f, f2, q, rho_hard, 00106 rho_soft, rho_total 00107 REAL(KIND=dp), DIMENSION(3) :: dg, dr 00108 REAL(KIND=dp), DIMENSION(:), POINTER :: f2sum, f2sum_gather, f4sum, 00109 f4sum_gather, fmax, fmax_gather, fmin, fmin_gather, fsum, fsum_gather, 00110 gsq, q_shell, q_shell_gather 00111 TYPE(atomic_kind_type), DIMENSION(:), 00112 POINTER :: atomic_kind_set 00113 TYPE(cp_para_env_type), POINTER :: para_env 00114 TYPE(dft_control_type), POINTER :: dft_control 00115 TYPE(particle_type), DIMENSION(:), 00116 POINTER :: particle_set 00117 TYPE(pw_env_type), POINTER :: pw_env 00118 TYPE(pw_p_type) :: rhotot_elec_gspace 00119 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00120 TYPE(qs_rho_type), POINTER :: rho 00121 TYPE(rho_atom_type), DIMENSION(:), 00122 POINTER :: rho_atom_set 00123 00124 failure = .FALSE. 00125 00126 CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure) 00127 00128 IF (.NOT.failure) THEN 00129 00130 CALL timeset(routineN,handle) 00131 00132 NULLIFY (atomic_kind_set) 00133 NULLIFY (aux_index) 00134 NULLIFY (auxbas_pw_pool) 00135 NULLIFY (dft_control) 00136 NULLIFY (f2sum) 00137 NULLIFY (f2sum_gather) 00138 NULLIFY (f4sum) 00139 NULLIFY (f4sum_gather) 00140 NULLIFY (fmax) 00141 NULLIFY (fmax_gather) 00142 NULLIFY (fmin) 00143 NULLIFY (fmin_gather) 00144 NULLIFY (fsum) 00145 NULLIFY (fsum_gather) 00146 NULLIFY (gsq) 00147 NULLIFY (ng_shell) 00148 NULLIFY (ng_shell_gather) 00149 NULLIFY (nshell_pe) 00150 NULLIFY (offset_pe) 00151 NULLIFY (para_env) 00152 NULLIFY (particle_set) 00153 NULLIFY (pw_env) 00154 NULLIFY (q_shell) 00155 NULLIFY (q_shell_gather) 00156 NULLIFY (rho) 00157 NULLIFY (rho_atom_set) 00158 00159 CALL cite_reference(Krack2002) 00160 00161 CALL get_qs_env(qs_env=qs_env,& 00162 atomic_kind_set=atomic_kind_set,& 00163 dft_control=dft_control,& 00164 para_env=para_env,& 00165 particle_set=particle_set,& 00166 pw_env=pw_env,& 00167 rho=rho,& 00168 rho_atom_set=rho_atom_set,& 00169 error=error) 00170 00171 CALL pw_env_get(pw_env=pw_env,& 00172 auxbas_pw_pool=auxbas_pw_pool,& 00173 error=error) 00174 00175 group = para_env%group 00176 npe = para_env%num_pe 00177 source = para_env%source 00178 00179 ! Plane waves grid to assemble the total electronic density 00180 00181 CALL pw_pool_create_pw(pool=auxbas_pw_pool,& 00182 pw=rhotot_elec_gspace%pw,& 00183 use_data=COMPLEXDATA1D,& 00184 in_space=RECIPROCALSPACE,& 00185 error=error) 00186 CALL pw_zero(rhotot_elec_gspace%pw,error=error) 00187 00188 CALL get_pw_grid_info(pw_grid=rhotot_elec_gspace%pw%pw_grid,& 00189 dr=dr,& 00190 npts=npts,& 00191 cutoff=cutoff,& 00192 ngpts=ngpts,& 00193 gsquare=gsq,& 00194 error=error) 00195 00196 dg(:) = twopi/(npts(:)*dr(:)) 00197 00198 ! Build the total electronic density in reciprocal space 00199 00200 CALL calculate_rhotot_elec_gspace(qs_env=qs_env,& 00201 auxbas_pw_pool=auxbas_pw_pool,& 00202 rhotot_elec_gspace=rhotot_elec_gspace,& 00203 q_max=q_max,& 00204 rho_hard=rho_hard,& 00205 rho_soft=rho_soft,& 00206 error=error) 00207 00208 rho_total = rho_hard + rho_soft 00209 00210 ! Calculate the coherent X-ray spectrum 00211 00212 ! Now we have to gather the data from all processes, since each 00213 ! process has only worked his sub-grid 00214 00215 ! Scan the g-vector shells 00216 00217 CALL reallocate(q_shell,1,nblock) 00218 CALL reallocate(ng_shell,1,nblock) 00219 00220 ng = SIZE(gsq) 00221 00222 jg = 1 00223 nshell = 1 00224 q_shell(1) = SQRT(gsq(1)) 00225 ng_shell(1) = 1 00226 00227 DO ig=2,ng 00228 CPPostcondition(gsq(ig)>=gsq(jg),cp_failure_level,routineP,error,failure) 00229 IF (ABS(gsq(ig) - gsq(jg)) > 1.0E-12_dp) THEN 00230 nshell = nshell + 1 00231 IF (nshell > SIZE(q_shell)) THEN 00232 CALL reallocate(q_shell,1,SIZE(q_shell)+nblock) 00233 CALL reallocate(ng_shell,1,SIZE(ng_shell)+nblock) 00234 END IF 00235 q = SQRT(gsq(ig)) 00236 IF (q > q_max) THEN 00237 nshell = nshell - 1 00238 EXIT 00239 END IF 00240 q_shell(nshell) = q 00241 ng_shell(nshell) = 1 00242 jg = ig 00243 ELSE 00244 ng_shell(nshell) = ng_shell(nshell) + 1 00245 END IF 00246 END DO 00247 00248 CALL reallocate(q_shell,1,nshell) 00249 CALL reallocate(ng_shell,1,nshell) 00250 CALL reallocate(fmin,1,nshell) 00251 CALL reallocate(fmax,1,nshell) 00252 CALL reallocate(fsum,1,nshell) 00253 CALL reallocate(f2sum,1,nshell) 00254 CALL reallocate(f4sum,1,nshell) 00255 00256 ig = 0 00257 DO ishell=1,nshell 00258 fmin(ishell) = HUGE(0.0_dp) 00259 fmax(ishell) = 0.0_dp 00260 fsum(ishell) = 0.0_dp 00261 f2sum(ishell) = 0.0_dp 00262 f4sum(ishell) = 0.0_dp 00263 DO ig_shell=1,ng_shell(ishell) 00264 f = ABS(rhotot_elec_gspace%pw%cc(ig+ig_shell)) 00265 fmin(ishell) = MIN(fmin(ishell),f) 00266 fmax(ishell) = MAX(fmax(ishell),f) 00267 fsum(ishell) = fsum(ishell) + f 00268 f2 = f*f 00269 f2sum(ishell) = f2sum(ishell) + f2 00270 f4sum(ishell) = f4sum(ishell) + f2*f2 00271 END DO 00272 ig = ig + ng_shell(ishell) 00273 END DO 00274 00275 CALL reallocate(nshell_pe,0,npe-1) 00276 CALL reallocate(offset_pe,0,npe-1) 00277 00278 ! Root (source) process gathers the number of shell of each process 00279 00280 CALL mp_gather(nshell,nshell_pe,source,group) 00281 00282 ! Only the root process which has to print the full spectrum has to 00283 ! allocate here the receive buffers with their real sizes 00284 00285 IF (unit_number > 0) THEN 00286 nshell_gather = SUM(nshell_pe) 00287 offset_pe(0) = 0 00288 DO ipe=1,npe-1 00289 offset_pe(ipe) = offset_pe(ipe-1) + nshell_pe(ipe-1) 00290 END DO 00291 ELSE 00292 nshell_gather = 1 ! dummy value for the non-root processes 00293 END IF 00294 00295 CALL reallocate(q_shell_gather,1,nshell_gather) 00296 CALL reallocate(ng_shell_gather,1,nshell_gather) 00297 CALL reallocate(fmin_gather,1,nshell_gather) 00298 CALL reallocate(fmax_gather,1,nshell_gather) 00299 CALL reallocate(fsum_gather,1,nshell_gather) 00300 CALL reallocate(f2sum_gather,1,nshell_gather) 00301 CALL reallocate(f4sum_gather,1,nshell_gather) 00302 00303 CALL mp_gatherv(q_shell,q_shell_gather,nshell_pe,offset_pe,source,group) 00304 CALL mp_gatherv(ng_shell,ng_shell_gather,nshell_pe,offset_pe,source,group) 00305 CALL mp_gatherv(fmax,fmax_gather,nshell_pe,offset_pe,source,group) 00306 CALL mp_gatherv(fmin,fmin_gather,nshell_pe,offset_pe,source,group) 00307 CALL mp_gatherv(fsum,fsum_gather,nshell_pe,offset_pe,source,group) 00308 CALL mp_gatherv(f2sum,f2sum_gather,nshell_pe,offset_pe,source,group) 00309 CALL mp_gatherv(f4sum,f4sum_gather,nshell_pe,offset_pe,source,group) 00310 00311 IF (ASSOCIATED(offset_pe)) THEN 00312 DEALLOCATE (offset_pe,STAT=istat) 00313 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00314 END IF 00315 00316 IF (ASSOCIATED(nshell_pe)) THEN 00317 DEALLOCATE (nshell_pe,STAT=istat) 00318 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00319 END IF 00320 00321 ! Print X-ray diffraction spectrum (I/O node only) 00322 00323 IF (unit_number > 0) THEN 00324 00325 CALL reallocate(aux_index,1,nshell_gather) 00326 00327 ! Sort the gathered shells 00328 00329 CALL sort(q_shell_gather,nshell_gather,aux_index) 00330 00331 ! Allocate final arrays of sufficient size, i.e. nshell_gather 00332 ! is always greater or equal the final nshell value 00333 00334 CALL reallocate(q_shell,1,nshell_gather) 00335 CALL reallocate(ng_shell,1,nshell_gather) 00336 CALL reallocate(fmin,1,nshell_gather) 00337 CALL reallocate(fmax,1,nshell_gather) 00338 CALL reallocate(fsum,1,nshell_gather) 00339 CALL reallocate(f2sum,1,nshell_gather) 00340 CALL reallocate(f4sum,1,nshell_gather) 00341 00342 jg = 1 00343 nshell = 1 00344 q_shell(1) = q_shell_gather(1) 00345 i = aux_index(1) 00346 ng_shell(1) = ng_shell_gather(i) 00347 fmin(1) = fmin_gather(i) 00348 fmax(1) = fmax_gather(i) 00349 fsum(1) = fsum_gather(i) 00350 f2sum(1) = f2sum_gather(i) 00351 f4sum(1) = f4sum_gather(i) 00352 00353 DO ig=2,nshell_gather 00354 i = aux_index(ig) 00355 IF (ABS(q_shell_gather(ig) - q_shell_gather(jg)) > 1.0E-12_dp) THEN 00356 nshell = nshell + 1 00357 q_shell(nshell) = q_shell_gather(ig) 00358 ng_shell(nshell) = ng_shell_gather(i) 00359 fmin(nshell) = fmin_gather(i) 00360 fmax(nshell) = fmax_gather(i) 00361 fsum(nshell) = fsum_gather(i) 00362 f2sum(nshell) = f2sum_gather(i) 00363 f4sum(nshell) = f4sum_gather(i) 00364 jg = ig 00365 ELSE 00366 ng_shell(nshell) = ng_shell(nshell) + ng_shell_gather(i) 00367 fmin(nshell) = MIN(fmin(nshell),fmin_gather(i)) 00368 fmax(nshell) = MAX(fmax(nshell),fmax_gather(i)) 00369 fsum(nshell) = fsum(nshell) + fsum_gather(i) 00370 f2sum(nshell) = f2sum(nshell) + f2sum_gather(i) 00371 f4sum(nshell) = f4sum(nshell) + f4sum_gather(i) 00372 END IF 00373 END DO 00374 00375 ! The auxiliary index array is no longer needed now 00376 00377 IF (ASSOCIATED(aux_index)) THEN 00378 DEALLOCATE (aux_index,STAT=istat) 00379 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00380 END IF 00381 00382 ! Allocate the final arrays for printing with their real size 00383 00384 CALL reallocate(q_shell,1,nshell) 00385 CALL reallocate(ng_shell,1,nshell) 00386 CALL reallocate(fmin,1,nshell) 00387 CALL reallocate(fmax,1,nshell) 00388 CALL reallocate(fsum,1,nshell) 00389 CALL reallocate(f2sum,1,nshell) 00390 CALL reallocate(f4sum,1,nshell) 00391 00392 ! Write the X-ray diffraction spectrum to the specified file 00393 00394 WRITE (UNIT=unit_number,FMT="(A)")& 00395 "#",& 00396 "# Coherent X-ray diffraction spectrum",& 00397 "#" 00398 WRITE (UNIT=unit_number,FMT="(A,1X,F20.10)")& 00399 "# Soft electronic charge (G-space) :",rho_soft,& 00400 "# Hard electronic charge (G-space) :",rho_hard,& 00401 "# Total electronic charge (G-space):",rho_total,& 00402 "# Density cutoff [Rydberg] :",2.0_dp*cutoff,& 00403 "# q(min) [1/Angstrom] :",q_shell(2)/angstrom,& 00404 "# q(max) [1/Angstrom] :",q_shell(nshell)/angstrom,& 00405 "# q(max) [1/Angstrom] (requested) :",q_max/angstrom 00406 WRITE (UNIT=unit_number,FMT="(A,2X,I8)")& 00407 "# Number of g-vectors (grid points):",ngpts,& 00408 "# Number of g-vector shells :",nshell 00409 WRITE (UNIT=unit_number,FMT="(A,3(1X,I6))")& 00410 "# Grid size (a,b,c) :",npts(1:3) 00411 WRITE (UNIT=unit_number,FMT="(A,3F7.3)")& 00412 "# dg [1/Angstrom] :",dg(1:3)/angstrom,& 00413 "# dr [Angstrom] :",dr(1:3)*angstrom 00414 WRITE (UNIT=unit_number,FMT="(A)")& 00415 "#",& 00416 "# shell points q [1/A] <|F(q)|^2> Min(|F(q)|) "//& 00417 " Max(|F(q)|) <|F(q)|>^2 <|F(q)|^4>" 00418 00419 DO ishell=1,nshell 00420 WRITE (UNIT=unit_number,FMT="(T2,I6,2X,I6,5(1X,F15.6),1X,ES15.6)")& 00421 ishell,& 00422 ng_shell(ishell),& 00423 q_shell(ishell)/angstrom,& 00424 f2sum(ishell)/REAL(ng_shell(ishell),KIND=dp), 00425 fmin(ishell), 00426 fmax(ishell), 00427 (fsum(ishell)/REAL(ng_shell(ishell),KIND=dp))**2, 00428 f4sum(ishell)/REAL(ng_shell(ishell),KIND=dp) 00429 END DO 00430 00431 END IF 00432 00433 ! Release work storage 00434 00435 IF (ASSOCIATED(fmin)) THEN 00436 DEALLOCATE (fmin,STAT=istat) 00437 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00438 END IF 00439 00440 IF (ASSOCIATED(fmax)) THEN 00441 DEALLOCATE (fmax,STAT=istat) 00442 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00443 END IF 00444 00445 IF (ASSOCIATED(fsum)) THEN 00446 DEALLOCATE (fsum,STAT=istat) 00447 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00448 END IF 00449 00450 IF (ASSOCIATED(f2sum)) THEN 00451 DEALLOCATE (f2sum,STAT=istat) 00452 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00453 END IF 00454 00455 IF (ASSOCIATED(f4sum)) THEN 00456 DEALLOCATE (f4sum,STAT=istat) 00457 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00458 END IF 00459 00460 IF (ASSOCIATED(ng_shell)) THEN 00461 DEALLOCATE (ng_shell,STAT=istat) 00462 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00463 END IF 00464 00465 IF (ASSOCIATED(q_shell)) THEN 00466 DEALLOCATE (q_shell,STAT=istat) 00467 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00468 END IF 00469 00470 IF (ASSOCIATED(fmin_gather)) THEN 00471 DEALLOCATE (fmin_gather,STAT=istat) 00472 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00473 END IF 00474 00475 IF (ASSOCIATED(fmax_gather)) THEN 00476 DEALLOCATE (fmax_gather,STAT=istat) 00477 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00478 END IF 00479 00480 IF (ASSOCIATED(fsum_gather)) THEN 00481 DEALLOCATE (fsum_gather,STAT=istat) 00482 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00483 END IF 00484 00485 IF (ASSOCIATED(f2sum_gather)) THEN 00486 DEALLOCATE (f2sum_gather,STAT=istat) 00487 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00488 END IF 00489 00490 IF (ASSOCIATED(f4sum_gather)) THEN 00491 DEALLOCATE (f4sum_gather,STAT=istat) 00492 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00493 END IF 00494 00495 IF (ASSOCIATED(ng_shell_gather)) THEN 00496 DEALLOCATE (ng_shell_gather,STAT=istat) 00497 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00498 END IF 00499 00500 IF (ASSOCIATED(q_shell_gather)) THEN 00501 DEALLOCATE (q_shell_gather,STAT=istat) 00502 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00503 END IF 00504 00505 CALL pw_pool_give_back_pw(auxbas_pw_pool,& 00506 rhotot_elec_gspace%pw,& 00507 error=error) 00508 00509 CALL timestop(handle) 00510 00511 END IF ! failure 00512 00513 END SUBROUTINE xray_diffraction_spectrum 00514 00515 ! ***************************************************************************** 00522 SUBROUTINE calculate_rhotot_elec_gspace(qs_env,auxbas_pw_pool,& 00523 rhotot_elec_gspace,q_max,rho_hard,& 00524 rho_soft,fsign,error) 00525 00526 TYPE(qs_environment_type), POINTER :: qs_env 00527 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00528 TYPE(pw_p_type), INTENT(INOUT) :: rhotot_elec_gspace 00529 REAL(KIND=dp), INTENT(IN) :: q_max 00530 REAL(KIND=dp), INTENT(OUT) :: rho_hard, rho_soft 00531 REAL(KIND=dp), INTENT(IN), OPTIONAL :: fsign 00532 TYPE(cp_error_type), INTENT(INOUT) :: error 00533 00534 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_rhotot_elec_gspace', 00535 routineP = moduleN//':'//routineN 00536 00537 INTEGER :: atom, handle, iatom, ico, ico1_pgf, ico1_set, ikind, ipgf, 00538 iset, iso, iso1_pgf, iso1_set, ison, ispin, istat, jco, jco1_pgf, 00539 jco1_set, jpgf, jset, jso, jso1_pgf, jso1_set, json, la, lb, maxco, 00540 maxso, na, natom, nb, ncoa, ncob, ncotot, nkind, nsatbas, nset, nsoa, 00541 nsob, nsotot, nspin 00542 INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf, 00543 o2nindex 00544 LOGICAL :: failure, orthorhombic, 00545 paw_atom 00546 REAL(KIND=dp) :: alpha, eps_rho_gspace, 00547 rho_total, scale, volume 00548 REAL(KIND=dp), DIMENSION(3) :: ra 00549 REAL(KIND=dp), DIMENSION(:, :), POINTER :: delta_cpc, pab, work, zet 00550 TYPE(atomic_kind_type), DIMENSION(:), 00551 POINTER :: atomic_kind_set 00552 TYPE(atomic_kind_type), POINTER :: atomic_kind 00553 TYPE(cell_type), POINTER :: cell 00554 TYPE(dft_control_type), POINTER :: dft_control 00555 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 00556 TYPE(particle_type), DIMENSION(:), 00557 POINTER :: particle_set 00558 TYPE(paw_proj_set_type), POINTER :: paw_proj 00559 TYPE(pw_p_type) :: rho_elec_gspace 00560 TYPE(qs_rho_type), POINTER :: rho 00561 TYPE(rho_atom_coeff), DIMENSION(:), 00562 POINTER :: cpc_h, cpc_s 00563 TYPE(rho_atom_type), DIMENSION(:), 00564 POINTER :: rho_atom_set 00565 TYPE(rho_atom_type), POINTER :: rho_atom 00566 00567 failure = .FALSE. 00568 00569 CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure) 00570 CPPrecondition(ASSOCIATED(auxbas_pw_pool),cp_failure_level,routineP,error,failure) 00571 00572 IF (.NOT.failure) THEN 00573 00574 CALL timeset(routineN,handle) 00575 00576 NULLIFY (atom_list) 00577 NULLIFY (atomic_kind) 00578 NULLIFY (atomic_kind_set) 00579 NULLIFY (cell) 00580 NULLIFY (cpc_h) 00581 NULLIFY (cpc_s) 00582 NULLIFY (delta_cpc) 00583 NULLIFY (dft_control) 00584 NULLIFY (lmax) 00585 NULLIFY (lmin) 00586 NULLIFY (npgf) 00587 NULLIFY (orb_basis_set) 00588 NULLIFY (pab) 00589 NULLIFY (particle_set) 00590 NULLIFY (rho) 00591 NULLIFY (rho_atom) 00592 NULLIFY (rho_atom_set) 00593 NULLIFY (work) 00594 NULLIFY (zet) 00595 00596 CALL get_qs_env(qs_env=qs_env,& 00597 atomic_kind_set=atomic_kind_set,& 00598 cell=cell,& 00599 dft_control=dft_control,& 00600 particle_set=particle_set,& 00601 rho=rho,& 00602 rho_atom_set=rho_atom_set,& 00603 error=error) 00604 00605 eps_rho_gspace = dft_control%qs_control%eps_rho_gspace 00606 nkind = SIZE(atomic_kind_set) 00607 nspin = dft_control%nspins 00608 00609 ! Load the soft contribution of the electronic density 00610 00611 CALL pw_pool_create_pw(pool=auxbas_pw_pool,& 00612 pw=rho_elec_gspace%pw,& 00613 use_data=COMPLEXDATA1D,& 00614 in_space=RECIPROCALSPACE,& 00615 error=error) 00616 00617 CALL pw_zero(rhotot_elec_gspace%pw,error=error) 00618 00619 DO ispin=1,nspin 00620 CALL pw_zero(rho_elec_gspace%pw,error=error) 00621 CALL pw_transfer(rho%rho_r(ispin)%pw,rho_elec_gspace%pw,debug=.FALSE.,error=error) 00622 IF (PRESENT(fsign).AND.(ispin == 2)) THEN 00623 alpha = fsign 00624 ELSE 00625 alpha = 1.0_dp 00626 END IF 00627 CALL pw_axpy(rho_elec_gspace%pw,rhotot_elec_gspace%pw,alpha=alpha,error=error) 00628 END DO 00629 00630 ! Release the auxiliary PW grid for the calculation of the soft 00631 ! contribution 00632 00633 CALL pw_pool_give_back_pw(pool=auxbas_pw_pool,& 00634 pw=rho_elec_gspace%pw,& 00635 error=error) 00636 00637 rho_soft = pw_integrate_function(rhotot_elec_gspace%pw,isign=-1,error=error) 00638 00639 CALL get_pw_grid_info(pw_grid=rhotot_elec_gspace%pw%pw_grid,& 00640 vol=volume,& 00641 orthorhombic=orthorhombic,& 00642 error=error) 00643 CPPostcondition(orthorhombic,cp_failure_level,routineP,error,failure) 00644 00645 CALL pw_scale(rhotot_elec_gspace%pw,volume,error=error) 00646 00647 ! Add the hard contribution of the electronic density 00648 00649 ! Each process has to loop over all PAW atoms, since the g-space grid 00650 ! is already distributed over all processes 00651 00652 DO ikind=1,nkind 00653 00654 atomic_kind => atomic_kind_set(ikind) 00655 00656 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00657 atom_list=atom_list,& 00658 natom=natom,& 00659 orb_basis_set=orb_basis_set,& 00660 paw_proj_set=paw_proj,& 00661 paw_atom=paw_atom) 00662 00663 IF (.NOT.paw_atom) CYCLE ! no PAW atom: nothing to do 00664 00665 CALL get_paw_proj_set(paw_proj_set=paw_proj,o2nindex=o2nindex,nsatbas=nsatbas) 00666 00667 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 00668 lmax=lmax,& 00669 lmin=lmin,& 00670 maxco=maxco,& 00671 maxso=maxso,& 00672 npgf=npgf,& 00673 nset=nset,& 00674 zet=zet) 00675 00676 ncotot = maxco*nset 00677 nsotot = maxso*nset 00678 CALL reallocate(delta_cpc,1,nsatbas,1,nsatbas) 00679 CALL reallocate(pab,1,ncotot,1,ncotot) 00680 CALL reallocate(work,1,maxso,1,maxco) 00681 00682 DO iatom=1,natom 00683 00684 atom = atom_list(iatom) 00685 rho_atom => rho_atom_set(atom) 00686 00687 CALL get_rho_atom(rho_atom=rho_atom,& 00688 cpc_h=cpc_h,& 00689 cpc_s=cpc_s) 00690 00691 ra(:) = pbc(particle_set(iatom)%r,cell) 00692 00693 delta_cpc = 0.0_dp 00694 00695 DO ispin=1,nspin 00696 IF (PRESENT(fsign).AND.(ispin == 2)) THEN 00697 alpha = fsign 00698 ELSE 00699 alpha = 1.0_dp 00700 END IF 00701 delta_cpc = delta_cpc + alpha*(cpc_h(ispin)%r_coef - cpc_s(ispin)%r_coef) 00702 END DO 00703 00704 scale = 1.0_dp 00705 00706 DO iset=1,nset 00707 ico1_set = (iset - 1)*maxco + 1 00708 iso1_set = (iset - 1)*maxso + 1 00709 ncoa = ncoset(lmax(iset)) 00710 nsoa = nsoset(lmax(iset)) 00711 DO jset=1,nset 00712 jco1_set = (jset - 1)*maxco + 1 00713 jso1_set = (jset - 1)*maxso + 1 00714 ncob = ncoset(lmax(jset)) 00715 nsob = nsoset(lmax(jset)) 00716 DO ipgf=1,npgf(iset) 00717 ico1_pgf = ico1_set + (ipgf - 1)*ncoa 00718 iso1_pgf = iso1_set + (ipgf - 1)*nsoa 00719 DO jpgf=1,npgf(jset) 00720 jco1_pgf = jco1_set + (jpgf - 1)*ncob 00721 jso1_pgf = jso1_set + (jpgf - 1)*nsob 00722 ico = ico1_pgf + ncoset(lmin(iset)-1) 00723 iso = iso1_pgf + nsoset(lmin(iset)-1) 00724 00725 ! Transformation spherical to Cartesian 00726 00727 DO la=lmin(iset),lmax(iset) 00728 jco = jco1_pgf + ncoset(lmin(jset)-1) 00729 jso = jso1_pgf + nsoset(lmin(jset)-1) 00730 DO lb=lmin(jset),lmax(jset) 00731 ison = o2nindex(iso) 00732 json = o2nindex(jso) 00733 CALL dgemm("N","N",nso(la),nco(lb),nso(lb),1.0_dp,& 00734 delta_cpc(ison,json),SIZE(delta_cpc,1),& 00735 orbtramat(lb)%slm,nso(lb),0.0_dp,work,& 00736 maxso) 00737 CALL dgemm("T","N",nco(la),nco(lb),nso(la),1.0_dp,& 00738 orbtramat(la)%slm,nso(la),work,maxso,& 00739 0.0_dp,pab(ico,jco),SIZE(pab,1)) 00740 jco = jco + nco(lb) 00741 jso = jso + nso(lb) 00742 END DO ! next lb 00743 ico = ico + nco(la) 00744 iso = iso + nso(la) 00745 END DO ! next la 00746 00747 ! Collocate current product of primitive Cartesian functions 00748 00749 na = ico1_pgf - 1 00750 nb = jco1_pgf - 1 00751 00752 CALL collocate_pgf_product_gspace(& 00753 la_max=lmax(iset),& 00754 zeta=zet(ipgf,iset),& 00755 la_min=lmin(iset),& 00756 lb_max=lmax(jset),& 00757 zetb=zet(jpgf,jset),& 00758 lb_min=lmin(jset),& 00759 ra=ra,& 00760 rab=(/0.0_dp,0.0_dp,0.0_dp/),& 00761 rab2=0.0_dp,& 00762 scale=scale,& 00763 pab=pab,& 00764 na=na,& 00765 nb=nb,& 00766 eps_rho_gspace=eps_rho_gspace,& 00767 gsq_max=q_max*q_max,& 00768 pw=rhotot_elec_gspace%pw) 00769 00770 END DO ! next primitive Gaussian function "jpgf" 00771 END DO ! next primitive Gaussian function "ipgf" 00772 END DO ! next shell set "jset" 00773 END DO ! next shell set "iset" 00774 END DO ! next atom "iatom" of atomic kind "ikind" 00775 END DO ! next atomic kind "ikind" 00776 00777 rho_total = pw_integrate_function(rhotot_elec_gspace%pw,isign=-1,error=error)/volume 00778 00779 rho_hard = rho_total - rho_soft 00780 00781 ! Release work storage 00782 00783 IF (ASSOCIATED(delta_cpc)) THEN 00784 DEALLOCATE (delta_cpc,STAT=istat) 00785 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00786 END IF 00787 00788 IF (ASSOCIATED(work)) THEN 00789 DEALLOCATE (work,STAT=istat) 00790 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00791 END IF 00792 00793 IF (ASSOCIATED(pab)) THEN 00794 DEALLOCATE (pab,STAT=istat) 00795 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00796 END IF 00797 00798 CALL timestop(handle) 00799 00800 END IF ! failure 00801 00802 END SUBROUTINE calculate_rhotot_elec_gspace 00803 00804 END MODULE xray_diffraction
1.7.3