|
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 MODULE qs_rho0_types 00007 00008 USE cp_units, ONLY: cp_unit_from_cp2k 00009 USE erf_fn, ONLY: erf 00010 USE f77_blas 00011 USE kinds, ONLY: dp,& 00012 dp_size,& 00013 int_size 00014 USE mathconstants, ONLY: fourpi,& 00015 pi,& 00016 rootpi 00017 USE memory_utilities, ONLY: reallocate 00018 USE pw_types, ONLY: pw_p_type,& 00019 pw_release 00020 USE qs_grid_atom, ONLY: grid_atom_type 00021 USE qs_rho_atom_types, ONLY: rho_atom_coeff 00022 USE termination, ONLY: stop_memory,& 00023 stop_program 00024 USE whittaker, ONLY: whittaker_c0a,& 00025 whittaker_ci 00026 #include "cp_common_uses.h" 00027 00028 IMPLICIT NONE 00029 00030 PRIVATE 00031 00032 ! *** Global parameters (only in this module) 00033 00034 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho0_types' 00035 00036 ! *** Define multipole type *** 00037 00038 ! ***************************************************************************** 00039 TYPE mpole_rho_atom 00040 REAL(dp), DIMENSION(:), POINTER :: Qlm_h, 00041 Qlm_s, 00042 Qlm_tot, 00043 Qlm_car 00044 REAL(dp) :: Qlm_z 00045 END TYPE mpole_rho_atom 00046 00047 ! ***************************************************************************** 00048 TYPE mpole_gau_overlap 00049 REAL(dp), DIMENSION(:,:,:), POINTER :: Qlm_gg 00050 REAL(dp), DIMENSION(:,:), POINTER :: g0_h, vg0_h 00051 REAL(dp) :: rpgf0_h, rpgf0_s 00052 END TYPE mpole_gau_overlap 00053 00054 ! ***************************************************************************** 00055 TYPE rho0_mpole_type 00056 TYPE(mpole_rho_atom), DIMENSION(:), POINTER :: mp_rho 00057 TYPE(mpole_gau_overlap), DIMENSION(:), 00058 POINTER :: mp_gau 00059 REAL(dp) :: zet0_h, 00060 total_rho0_h 00061 REAL(dp) :: max_rpgf0_s 00062 REAL(dp), DIMENSION(:), POINTER :: norm_g0l_h 00063 INTEGER, DIMENSION(:), POINTER :: lmax0_kind 00064 INTEGER :: lmax_0,igrid_zet0_s 00065 TYPE(pw_p_type), POINTER :: rho0_s_rs, 00066 rho0_s_gs 00067 END TYPE rho0_mpole_type 00068 00069 ! ***************************************************************************** 00070 TYPE rho0_atom_type 00071 TYPE(rho_atom_coeff), POINTER :: rho0_rad_h, 00072 vrho0_rad_h 00073 END TYPE rho0_atom_type 00074 00075 ! Public Types 00076 00077 PUBLIC :: mpole_rho_atom, mpole_gau_overlap, & 00078 rho0_atom_type, rho0_mpole_type 00079 00080 ! Public Subroutine 00081 00082 PUBLIC :: allocate_multipoles, allocate_rho0_mpole, & 00083 allocate_rho0_atom, allocate_rho0_atom_rad, & 00084 deallocate_rho0_atom, deallocate_rho0_mpole, & 00085 calculate_g0, get_rho0_mpole, initialize_mpole_rho, & 00086 write_rho0_info 00087 00088 CONTAINS 00089 00090 ! ***************************************************************************** 00091 SUBROUTINE allocate_multipoles(mp_rho,natom,mp_gau,nkind) 00092 00093 TYPE(mpole_rho_atom), DIMENSION(:), 00094 POINTER :: mp_rho 00095 INTEGER, INTENT(IN) :: natom 00096 TYPE(mpole_gau_overlap), DIMENSION(:), 00097 POINTER :: mp_gau 00098 INTEGER, INTENT(IN) :: nkind 00099 00100 CHARACTER(len=*), PARAMETER :: routineN = 'allocate_multipoles', 00101 routineP = moduleN//':'//routineN 00102 00103 INTEGER :: iat, ikind, istat 00104 00105 IF(ASSOCIATED(mp_rho)) THEN 00106 CALL deallocate_mpole_rho(mp_rho) 00107 END IF 00108 00109 ALLOCATE (mp_rho(natom),STAT=istat) 00110 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00111 "mp_rho",natom*int_size) 00112 00113 DO iat = 1,natom 00114 NULLIFY(mp_rho(iat)%Qlm_h) 00115 NULLIFY(mp_rho(iat)%Qlm_s) 00116 NULLIFY(mp_rho(iat)%Qlm_tot) 00117 NULLIFY(mp_rho(iat)%Qlm_car) 00118 END DO 00119 00120 IF(ASSOCIATED(mp_gau)) THEN 00121 CALL deallocate_mpole_gau(mp_gau) 00122 END IF 00123 00124 ALLOCATE (mp_gau(nkind),STAT=istat) 00125 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00126 "mp_gau",nkind*int_size) 00127 00128 DO ikind = 1,nkind 00129 NULLIFY(mp_gau(ikind)%Qlm_gg) 00130 NULLIFY(mp_gau(ikind)%g0_h) 00131 NULLIFY(mp_gau(ikind)%vg0_h) 00132 mp_gau(ikind)%rpgf0_h = 0.0_dp 00133 mp_gau(ikind)%rpgf0_s = 0.0_dp 00134 END DO 00135 00136 END SUBROUTINE allocate_multipoles 00137 00138 ! ***************************************************************************** 00139 SUBROUTINE allocate_rho0_atom(rho0_set,natom) 00140 00141 TYPE(rho0_atom_type), DIMENSION(:), 00142 POINTER :: rho0_set 00143 INTEGER :: natom 00144 00145 CHARACTER(len=*), PARAMETER :: routineN = 'allocate_rho0_atom', 00146 routineP = moduleN//':'//routineN 00147 00148 INTEGER :: iat, istat 00149 00150 IF(ASSOCIATED(rho0_set)) THEN 00151 CALL deallocate_rho0_atom(rho0_set) 00152 END IF 00153 00154 ALLOCATE (rho0_set(natom),STAT=istat) 00155 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00156 "rho0_atom_set",natom*int_size) 00157 00158 DO iat = 1,natom 00159 NULLIFY(rho0_set(iat)%rho0_rad_h) 00160 NULLIFY(rho0_set(iat)%vrho0_rad_h) 00161 ENDDO 00162 00163 END SUBROUTINE allocate_rho0_atom 00164 00165 ! ***************************************************************************** 00166 SUBROUTINE allocate_rho0_atom_rad(rho0_atom,nr,nchannels) 00167 00168 TYPE(rho0_atom_type) :: rho0_atom 00169 INTEGER :: nr, nchannels 00170 00171 CHARACTER(len=*), PARAMETER :: routineN = 'allocate_rho0_atom_rad', 00172 routineP = moduleN//':'//routineN 00173 00174 INTEGER :: istat 00175 00176 ALLOCATE(rho0_atom%rho0_rad_h,STAT=istat) 00177 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00178 "rho0_rad_h",int_size) 00179 00180 NULLIFY(rho0_atom%rho0_rad_h%r_coef) 00181 ALLOCATE (rho0_atom%rho0_rad_h%r_coef(1:nr,1:nchannels),STAT=istat) 00182 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00183 "rho0_rad_h",int_size) 00184 rho0_atom%rho0_rad_h%r_coef = 0.0_dp 00185 00186 ALLOCATE(rho0_atom%vrho0_rad_h,STAT=istat) 00187 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00188 "vrho0_rad_h",int_size) 00189 00190 NULLIFY(rho0_atom%vrho0_rad_h%r_coef) 00191 ALLOCATE (rho0_atom%vrho0_rad_h%r_coef(1:nr,1:nchannels),STAT=istat) 00192 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00193 "vrho0_rad_h",int_size) 00194 rho0_atom%vrho0_rad_h%r_coef = 0.0_dp 00195 00196 END SUBROUTINE allocate_rho0_atom_rad 00197 00198 ! ***************************************************************************** 00199 SUBROUTINE allocate_rho0_mpole(rho0,error) 00200 00201 TYPE(rho0_mpole_type), POINTER :: rho0 00202 TYPE(cp_error_type), INTENT(inout) :: error 00203 00204 CHARACTER(len=*), PARAMETER :: routineN = 'allocate_rho0_mpole', 00205 routineP = moduleN//':'//routineN 00206 00207 INTEGER :: istat 00208 00209 IF(ASSOCIATED(rho0)) THEN 00210 CALL deallocate_rho0_mpole(rho0,error=error) 00211 END IF 00212 00213 ALLOCATE (rho0,STAT=istat) 00214 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00215 "rho0_mpole",int_size) 00216 00217 NULLIFY(rho0%mp_rho) 00218 NULLIFY(rho0%mp_gau) 00219 NULLIFY(rho0%norm_g0l_h) 00220 NULLIFY(rho0%lmax0_kind) 00221 NULLIFY(rho0%rho0_s_rs) 00222 NULLIFY(rho0%rho0_s_gs) 00223 00224 END SUBROUTINE allocate_rho0_mpole 00225 00226 ! ***************************************************************************** 00227 SUBROUTINE calculate_g0(rho0_mpole,grid_atom,ik) 00228 00229 TYPE(rho0_mpole_type), POINTER :: rho0_mpole 00230 TYPE(grid_atom_type), POINTER :: grid_atom 00231 INTEGER :: ik 00232 00233 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_g0', 00234 routineP = moduleN//':'//routineN 00235 00236 INTEGER :: ir, istat, l, lmax, nr 00237 REAL(dp) :: c1, prefactor, root_z_h, z_h 00238 REAL(dp), ALLOCATABLE, DIMENSION(:) :: erf_z_h, gexp, gh_tmp, int1, 00239 int2 00240 00241 nr = grid_atom%nr 00242 lmax = rho0_mpole%lmax0_kind(ik) 00243 z_h = rho0_mpole%zet0_h 00244 root_z_h = SQRT(z_h) 00245 00246 ! Allocate g0 00247 CALL reallocate(rho0_mpole%mp_gau(ik)%g0_h,1,nr,0,lmax) 00248 CALL reallocate(rho0_mpole%mp_gau(ik)%vg0_h,1,nr,0,lmax) 00249 00250 ALLOCATE(gexp(nr),gh_tmp(nr),erf_z_h(nr),int1(nr),int2(nr),STAT=istat) 00251 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00252 "gexp,gh_tmp,erf_z_h,in1,int2",6*nr*dp_size) 00253 00254 gh_tmp(1:nr) = EXP(-z_h*grid_atom%rad2(1:nr)) 00255 00256 DO ir = 1,nr 00257 erf_z_h(ir) = erf(grid_atom%rad(ir)*root_z_h) 00258 END DO 00259 00260 DO ir =1,nr 00261 IF(gh_tmp(ir) < 1.0E-30_dp) gh_tmp(ir) = 0.0_dp 00262 END DO 00263 00264 gexp(1:nr) = gh_tmp(1:nr) 00265 rho0_mpole%mp_gau(ik)%g0_h(1:nr,0) = gh_tmp(1:nr)* & 00266 rho0_mpole%norm_g0l_h(0) 00267 CALL whittaker_c0a(int1,grid_atom%rad,gh_tmp,erf_z_h,z_h,0,0,nr) 00268 CALL whittaker_ci(int2,grid_atom%rad,gh_tmp,z_h,0,nr) 00269 00270 prefactor = fourpi*rho0_mpole%norm_g0l_h(0) 00271 00272 c1 = SQRT(pi*pi*pi/(z_h*z_h*z_h))*rho0_mpole%norm_g0l_h(0) 00273 00274 DO ir = 1,nr 00275 rho0_mpole%mp_gau(ik)%vg0_h(ir,0) = c1*erf_z_h(ir)*grid_atom%oorad2l(ir,1) 00276 END DO 00277 00278 DO l = 1,lmax 00279 gh_tmp(1:nr) = gh_tmp(1:nr)*grid_atom%rad(1:nr) 00280 rho0_mpole%mp_gau(ik)%g0_h(1:nr,l) = gh_tmp(1:nr)* & 00281 rho0_mpole%norm_g0l_h(l) 00282 00283 prefactor = fourpi/(2.0_dp*l+1.0_dp)*rho0_mpole%norm_g0l_h(l) 00284 CALL whittaker_c0a(int1,grid_atom%rad,gexp,erf_z_h,z_h,l,l,nr) 00285 DO ir = 1,nr 00286 rho0_mpole%mp_gau(ik)%vg0_h(ir,l) = prefactor*(int1(ir)+& 00287 int2(ir)*grid_atom%rad2l(ir,l)) 00288 END DO 00289 00290 END DO ! l 00291 00292 DEALLOCATE (gexp,erf_z_h,gh_tmp,int1,int2,STAT=istat) 00293 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00294 "gexp,gh_tmp,erf_z_h,int1,int2") 00295 END SUBROUTINE calculate_g0 00296 00297 ! ***************************************************************************** 00298 SUBROUTINE deallocate_mpole_gau(mp_gau) 00299 00300 TYPE(mpole_gau_overlap), DIMENSION(:), 00301 POINTER :: mp_gau 00302 00303 CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_mpole_gau', 00304 routineP = moduleN//':'//routineN 00305 00306 INTEGER :: ikind, istat, nkind 00307 00308 istat = 0 00309 nkind = SIZE(mp_gau) 00310 00311 DO ikind = 1,nkind 00312 00313 IF(ASSOCIATED(mp_gau(ikind)%Qlm_gg)) THEN 00314 DEALLOCATE(mp_gau(ikind)%Qlm_gg,STAT=istat) 00315 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00316 "Qlm_gg") 00317 END IF 00318 00319 DEALLOCATE(mp_gau(ikind)%g0_h, STAT=istat) 00320 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00321 "mp_gau(ikind)%g0_h") 00322 00323 DEALLOCATE(mp_gau(ikind)%vg0_h, STAT=istat) 00324 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00325 "mp_gau(ikind)%vg0_h") 00326 END DO 00327 00328 DEALLOCATE(mp_gau, STAT=istat) 00329 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00330 "mp_gau") 00331 00332 END SUBROUTINE deallocate_mpole_gau 00333 00334 ! ***************************************************************************** 00335 SUBROUTINE deallocate_mpole_rho(mp_rho) 00336 00337 TYPE(mpole_rho_atom), DIMENSION(:), 00338 POINTER :: mp_rho 00339 00340 CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_mpole_rho', 00341 routineP = moduleN//':'//routineN 00342 00343 INTEGER :: iat, istat, natom 00344 00345 natom = SIZE(mp_rho) 00346 00347 DO iat = 1,natom 00348 DEALLOCATE(mp_rho(iat)%Qlm_h,STAT=istat) 00349 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00350 "Qlm_h") 00351 DEALLOCATE(mp_rho(iat)%Qlm_s,STAT=istat) 00352 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00353 "Qlm_s") 00354 DEALLOCATE(mp_rho(iat)%Qlm_tot,STAT=istat) 00355 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00356 "Qlm_tot") 00357 DEALLOCATE(mp_rho(iat)%Qlm_car,STAT=istat) 00358 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00359 "Qlm_car") 00360 END DO 00361 00362 DEALLOCATE(mp_rho,STAT=istat) 00363 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00364 "mp_rho") 00365 00366 END SUBROUTINE deallocate_mpole_rho 00367 00368 ! ***************************************************************************** 00369 SUBROUTINE deallocate_rho0_atom(rho0_atom_set) 00370 00371 TYPE(rho0_atom_type), DIMENSION(:), 00372 POINTER :: rho0_atom_set 00373 00374 CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_rho0_atom', 00375 routineP = moduleN//':'//routineN 00376 00377 INTEGER :: iat, istat, natom 00378 00379 IF (ASSOCIATED(rho0_atom_set)) THEN 00380 00381 natom = SIZE(rho0_atom_set) 00382 00383 DO iat = 1,natom 00384 IF(ASSOCIATED(rho0_atom_set(iat)%rho0_rad_h)) THEN 00385 DEALLOCATE (rho0_atom_set(iat)%rho0_rad_h%r_coef) 00386 DEALLOCATE (rho0_atom_set(iat)%rho0_rad_h) 00387 ENDIF 00388 IF(ASSOCIATED(rho0_atom_set(iat)%vrho0_rad_h)) THEN 00389 DEALLOCATE (rho0_atom_set(iat)%vrho0_rad_h%r_coef) 00390 DEALLOCATE (rho0_atom_set(iat)%vrho0_rad_h) 00391 ENDIF 00392 ENDDO 00393 00394 DEALLOCATE (rho0_atom_set,STAT=istat) 00395 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00396 "rho0_atom_set") 00397 ELSE 00398 CALL stop_program(routineN,moduleN,__LINE__,& 00399 "The pointer rho0_atom_set is not associated and "//& 00400 "cannot be deallocated") 00401 END IF 00402 00403 END SUBROUTINE deallocate_rho0_atom 00404 ! ***************************************************************************** 00405 SUBROUTINE deallocate_rho0_mpole(rho0,error) 00406 00407 TYPE(rho0_mpole_type), POINTER :: rho0 00408 TYPE(cp_error_type), INTENT(inout) :: error 00409 00410 CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_rho0_mpole', 00411 routineP = moduleN//':'//routineN 00412 00413 INTEGER :: istat 00414 00415 IF (ASSOCIATED(rho0)) THEN 00416 00417 IF (ASSOCIATED(rho0%mp_gau)) CALL deallocate_mpole_gau(rho0%mp_gau) 00418 00419 IF (ASSOCIATED(rho0%mp_rho)) CALL deallocate_mpole_rho(rho0%mp_rho) 00420 00421 IF (ASSOCIATED(rho0%lmax0_kind)) THEN 00422 DEALLOCATE(rho0%lmax0_kind,STAT=istat) 00423 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00424 "rho0%lmax0_kind") 00425 END IF 00426 00427 IF (ASSOCIATED(rho0%norm_g0l_h)) THEN 00428 DEALLOCATE(rho0%norm_g0l_h, STAT=istat) 00429 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00430 "rho0%norm_g0l_h") 00431 END IF 00432 00433 IF (ASSOCIATED(rho0%rho0_s_rs)) THEN 00434 CALL pw_release(rho0%rho0_s_rs%pw,error=error) 00435 DEALLOCATE(rho0%rho0_s_rs, STAT=istat) 00436 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00437 "rho0%rho0_s_rs") 00438 ENDIF 00439 00440 IF (ASSOCIATED(rho0%rho0_s_gs)) THEN 00441 CALL pw_release(rho0%rho0_s_gs%pw,error=error) 00442 DEALLOCATE(rho0%rho0_s_gs, STAT=istat) 00443 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00444 "rho0%rho0_s_gs") 00445 00446 ENDIF 00447 DEALLOCATE (rho0,STAT=istat) 00448 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00449 "rho0_mpole") 00450 ELSE 00451 CALL stop_program(routineN,moduleN,__LINE__,& 00452 "The pointer rho0 is not associated and "//& 00453 "cannot be deallocated") 00454 END IF 00455 00456 END SUBROUTINE deallocate_rho0_mpole 00457 00458 ! ***************************************************************************** 00459 SUBROUTINE get_rho0_mpole(rho0_mpole, g0_h, vg0_h, iat, ikind, lmax_0, l0_ikind,& 00460 mp_gau_ikind, mp_rho, norm_g0l_h, & 00461 Qlm_gg, Qlm_car, Qlm_tot, & 00462 zet0_h, igrid_zet0_s, rpgf0_h,rpgf0_s, & 00463 max_rpgf0_s, rho0_s_rs, rho0_s_gs) 00464 00465 TYPE(rho0_mpole_type), POINTER :: rho0_mpole 00466 REAL(dp), DIMENSION(:, :), OPTIONAL, 00467 POINTER :: g0_h, vg0_h 00468 INTEGER, INTENT(IN), OPTIONAL :: iat, ikind 00469 INTEGER, INTENT(OUT), OPTIONAL :: lmax_0, l0_ikind 00470 TYPE(mpole_gau_overlap), OPTIONAL, 00471 POINTER :: mp_gau_ikind 00472 TYPE(mpole_rho_atom), DIMENSION(:), 00473 OPTIONAL, POINTER :: mp_rho 00474 REAL(dp), DIMENSION(:), OPTIONAL, 00475 POINTER :: norm_g0l_h 00476 REAL(dp), DIMENSION(:, :, :), OPTIONAL, 00477 POINTER :: Qlm_gg 00478 REAL(dp), DIMENSION(:), OPTIONAL, 00479 POINTER :: Qlm_car, Qlm_tot 00480 REAL(dp), INTENT(OUT), OPTIONAL :: zet0_h 00481 INTEGER, INTENT(OUT), OPTIONAL :: igrid_zet0_s 00482 REAL(dp), INTENT(OUT), OPTIONAL :: rpgf0_h, rpgf0_s, max_rpgf0_s 00483 TYPE(pw_p_type), OPTIONAL, POINTER :: rho0_s_rs, rho0_s_gs 00484 00485 CHARACTER(len=*), PARAMETER :: routineN = 'get_rho0_mpole', 00486 routineP = moduleN//':'//routineN 00487 00488 IF (ASSOCIATED(rho0_mpole)) THEN 00489 00490 IF(PRESENT(lmax_0)) lmax_0 = rho0_mpole%lmax_0 00491 IF(PRESENT(mp_rho)) mp_rho => rho0_mpole%mp_rho 00492 IF(PRESENT(norm_g0l_h)) norm_g0l_h => rho0_mpole%norm_g0l_h 00493 IF(PRESENT(zet0_h)) zet0_h = rho0_mpole%zet0_h 00494 IF(PRESENT(igrid_zet0_s)) igrid_zet0_s = rho0_mpole%igrid_zet0_s 00495 IF(PRESENT(max_rpgf0_s)) max_rpgf0_s = rho0_mpole%max_rpgf0_s 00496 IF(PRESENT(rho0_s_rs)) rho0_s_rs => rho0_mpole%rho0_s_rs 00497 IF(PRESENT(rho0_s_gs)) rho0_s_gs => rho0_mpole%rho0_s_gs 00498 00499 IF(PRESENT(ikind)) THEN 00500 IF(PRESENT(l0_ikind)) l0_ikind = rho0_mpole%lmax0_kind(ikind) 00501 IF(PRESENT(mp_gau_ikind)) mp_gau_ikind => rho0_mpole%mp_gau(ikind) 00502 IF(PRESENT(g0_h)) g0_h => rho0_mpole%mp_gau(ikind)%g0_h 00503 IF(PRESENT(vg0_h)) vg0_h => rho0_mpole%mp_gau(ikind)%vg0_h 00504 IF(PRESENT(Qlm_gg)) Qlm_gg => rho0_mpole%mp_gau(ikind)%Qlm_gg 00505 IF(PRESENT(rpgf0_h)) rpgf0_h = rho0_mpole%mp_gau(ikind)%rpgf0_h 00506 IF(PRESENT(rpgf0_s)) rpgf0_s = rho0_mpole%mp_gau(ikind)%rpgf0_s 00507 END IF 00508 IF(PRESENT(iat)) THEN 00509 IF(PRESENT(Qlm_car)) Qlm_car => rho0_mpole%mp_rho(iat)%Qlm_car 00510 IF(PRESENT(Qlm_tot)) Qlm_tot => rho0_mpole%mp_rho(iat)%Qlm_tot 00511 END IF 00512 00513 ELSE 00514 CALL stop_program(routineN,moduleN,__LINE__,& 00515 "The pointer rho0_mpole is not associated") 00516 END IF 00517 00518 END SUBROUTINE get_rho0_mpole 00519 00520 ! ***************************************************************************** 00521 SUBROUTINE initialize_mpole_rho(mp_rho,nchan_s,nchan_c,zeff,tddft) 00522 00523 TYPE(mpole_rho_atom) :: mp_rho 00524 INTEGER, INTENT(IN) :: nchan_s, nchan_c 00525 REAL(KIND=dp), INTENT(IN) :: zeff 00526 LOGICAL, OPTIONAL :: tddft 00527 00528 CHARACTER(len=*), PARAMETER :: routineN = 'initialize_mpole_rho', 00529 routineP = moduleN//':'//routineN 00530 00531 LOGICAL :: my_tddft 00532 00533 my_tddft = .FALSE. 00534 IF (PRESENT(tddft)) my_tddft = tddft 00535 00536 CALL reallocate(mp_rho%Qlm_h,1,nchan_s) 00537 CALL reallocate(mp_rho%Qlm_s,1,nchan_s) 00538 CALL reallocate(mp_rho%Qlm_tot,1,nchan_s) 00539 CALL reallocate(mp_rho%Qlm_car,1,nchan_c) 00540 00541 mp_rho%Qlm_h = 0.0_dp 00542 mp_rho%Qlm_s = 0.0_dp 00543 mp_rho%Qlm_tot = 0.0_dp 00544 mp_rho%Qlm_car = 0.0_dp 00545 IF (.NOT.my_tddft) THEN 00546 mp_rho%Qlm_z = -2.0_dp*rootpi*Zeff 00547 ELSE 00548 mp_rho%Qlm_z = 0.0_dp 00549 END IF 00550 00551 END SUBROUTINE initialize_mpole_rho 00552 00553 ! ***************************************************************************** 00554 SUBROUTINE write_rho0_info(rho0_mpole,unit_str,output_unit,error) 00555 00556 TYPE(rho0_mpole_type), POINTER :: rho0_mpole 00557 CHARACTER(LEN=*), INTENT(IN) :: unit_str 00558 INTEGER, INTENT(in) :: output_unit 00559 TYPE(cp_error_type), INTENT(inout) :: error 00560 00561 CHARACTER(len=*), PARAMETER :: routineN = 'write_rho0_info', 00562 routineP = moduleN//':'//routineN 00563 00564 INTEGER :: ikind, l, nkind 00565 REAL(dp) :: conv 00566 00567 IF (ASSOCIATED(rho0_mpole)) THEN 00568 conv = cp_unit_from_cp2k(1.0_dp,TRIM(unit_str),error=error) 00569 00570 WRITE (UNIT=output_unit,FMT="(/,T5,A,/)") & 00571 "*** Compensation density charges data set ***" 00572 WRITE (UNIT=output_unit,FMT="(T2,A,T35,f16.10)")& 00573 "- Rho0 exponent :",rho0_mpole%zet0_h 00574 WRITE (UNIT=output_unit,FMT="(T2,A,T35,I5)")& 00575 "- Global max l :",rho0_mpole%lmax_0 00576 00577 WRITE (UNIT=output_unit,FMT="(T2,A)")& 00578 "- Normalization constants for g0" 00579 DO l = 0,rho0_mpole%lmax_0 00580 WRITE (UNIT=output_unit,FMT="(T20,A,T31,I2,T38,A,f15.5)")& 00581 "ang. mom.= ", l, " hard= ", rho0_mpole%norm_g0l_h(l) 00582 END DO 00583 00584 nkind = SIZE(rho0_mpole%lmax0_kind,1) 00585 DO ikind = 1, nkind 00586 WRITE (UNIT=output_unit,FMT="(/,T2,A,T55,I2)")& 00587 "- rho0 max L and radii in "//TRIM(unit_str)//& 00588 " for the atom kind :", ikind 00589 00590 WRITE (UNIT=output_unit,FMT="(T2,T20,A,T55,I5)")& 00591 "=> l max :",rho0_mpole%lmax0_kind(ikind) 00592 00593 WRITE (UNIT=output_unit,FMT="(T2,T20,A,T55,f20.10)")& 00594 "=> max radius of g0: ",& 00595 rho0_mpole%mp_gau(ikind)%rpgf0_h*conv 00596 END DO ! ikind 00597 00598 ELSE 00599 WRITE(UNIT=output_unit,FMT="(/,T5,A,/)") & 00600 ' WARNING: I cannot print rho0, it is not associated' 00601 END IF 00602 00603 END SUBROUTINE write_rho0_info 00604 END MODULE qs_rho0_types
1.7.3