|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00012 MODULE external_potential_types 00013 USE bibliography, ONLY: Goedecker1996,& 00014 Hartwigsen1998,& 00015 Krack2000,& 00016 Krack2005,& 00017 cite_reference 00018 USE cp_linked_list_val, ONLY: cp_sll_val_next,& 00019 cp_sll_val_type 00020 USE cp_para_types, ONLY: cp_para_env_type 00021 USE cp_parser_methods, ONLY: parser_get_next_line,& 00022 parser_get_object,& 00023 parser_search_string,& 00024 parser_test_next_token 00025 USE cp_parser_types, ONLY: cp_parser_type,& 00026 parser_create,& 00027 parser_release 00028 USE f77_blas 00029 USE input_section_types, ONLY: section_vals_get,& 00030 section_vals_list_get,& 00031 section_vals_type,& 00032 section_vals_val_get,& 00033 section_vals_val_set 00034 USE input_val_types, ONLY: val_get,& 00035 val_type 00036 USE kinds, ONLY: default_path_length,& 00037 default_string_length,& 00038 dp 00039 USE mathconstants, ONLY: dfac,& 00040 fac,& 00041 pi,& 00042 rootpi 00043 USE mathlib, ONLY: symmetrize_matrix 00044 USE memory_utilities, ONLY: reallocate 00045 USE orbital_pointers, ONLY: co,& 00046 coset,& 00047 init_orbital_pointers,& 00048 nco,& 00049 ncoset,& 00050 nso 00051 USE orbital_transformation_matrices, ONLY: orbtramat 00052 USE periodic_table, ONLY: ptable 00053 USE string_utilities, ONLY: remove_word,& 00054 uppercase 00055 USE termination, ONLY: stop_program 00056 USE timings, ONLY: timeset,& 00057 timestop 00058 #include "cp_common_uses.h" 00059 00060 IMPLICIT NONE 00061 00062 PRIVATE 00063 00064 ! Global parameters 00065 00066 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'external_potential_types' 00067 00068 ! Define the all-electron potential type 00069 00070 ! Literature: M. Krack and M. Parrinello, 00071 ! Phys. Chem. Chem. Phys. 2, 2105 (2000) 00072 00073 ! ***************************************************************************** 00074 TYPE all_potential_type 00075 !MK PRIVATE 00076 CHARACTER(LEN=default_string_length) :: name 00077 CHARACTER(LEN=default_string_length), 00078 DIMENSION(2) :: description 00079 REAL(KIND = dp) :: alpha_core_charge, 00080 ccore_charge, 00081 core_charge_radius, 00082 zeff, zeff_correction 00083 INTEGER :: z 00084 INTEGER, DIMENSION(:), POINTER :: elec_conf 00085 END TYPE all_potential_type 00086 00087 ! Define the effective charge & inducible dipole potential type (for Fist) 00088 ! ***************************************************************************** 00089 TYPE fist_potential_type 00090 PRIVATE 00091 CHARACTER(LEN=default_string_length) :: name 00092 CHARACTER(LEN=default_string_length), 00093 DIMENSION(1) :: description 00094 REAL(KIND = dp) :: apol, cpol, mm_radius, qeff, 00095 qmmm_corr_radius, qmmm_radius 00096 00097 END TYPE fist_potential_type 00098 00099 ! Define the GTH potential type 00100 00101 ! Literature: - S. Goedecker, M. Teter and J. Hutter, 00102 ! Phys. Rev. B 54, 1703 (1996) 00103 ! - C. Hartwigsen, S. Goedecker and J. Hutter, 00104 ! Phys. Rev. B 58, 3641 (1998) 00105 ! - M. Krack, 00106 ! Theor. Chem. Acc. 114, 145 (2005) 00107 00108 ! ***************************************************************************** 00109 TYPE gth_potential_type 00110 !PRIVATE 00111 CHARACTER(LEN=default_string_length) :: name 00112 CHARACTER(LEN=default_string_length), 00113 DIMENSION(4) :: description 00114 REAL(KIND = dp) :: alpha_core_charge, 00115 alpha_ppl,ccore_charge, 00116 cerf_ppl,zeff, 00117 core_charge_radius, 00118 ppl_radius,ppnl_radius, 00119 zeff_correction 00120 INTEGER :: lppnl,lprj_ppnl_max, 00121 nexp_ppl,nppnl, 00122 nprj_ppnl_max,z 00123 REAL(KIND = dp), DIMENSION(:), POINTER :: alpha_ppnl,cexp_ppl 00124 INTEGER, DIMENSION(:), POINTER :: elec_conf 00125 ! nonlocal projectors 00126 INTEGER, DIMENSION(:), POINTER :: nprj_ppnl 00127 REAL(KIND = dp), DIMENSION(:,:), POINTER :: cprj,cprj_ppnl,vprj_ppnl 00128 REAL(KIND = dp), DIMENSION(:,:,:), POINTER :: hprj_ppnl 00129 ! type extensions 00130 ! NLCC 00131 LOGICAL :: nlcc 00132 INTEGER :: nexp_nlcc 00133 REAL(KIND = dp), DIMENSION(:), POINTER :: alpha_nlcc 00134 INTEGER, DIMENSION(:), POINTER :: nct_nlcc 00135 REAL(KIND = dp), DIMENSION(:,:), POINTER :: cval_nlcc 00136 ! LSD potential 00137 LOGICAL :: lsdpot 00138 INTEGER :: nexp_lsd 00139 REAL(KIND = dp), DIMENSION(:), POINTER :: alpha_lsd 00140 INTEGER, DIMENSION(:), POINTER :: nct_lsd 00141 REAL(KIND = dp), DIMENSION(:,:), POINTER :: cval_lsd 00142 ! extended local potential 00143 LOGICAL :: lpotextended 00144 INTEGER :: nexp_lpot 00145 REAL(KIND = dp), DIMENSION(:), POINTER :: alpha_lpot 00146 INTEGER, DIMENSION(:), POINTER :: nct_lpot 00147 REAL(KIND = dp), DIMENSION(:,:), POINTER :: cval_lpot 00148 00149 END TYPE gth_potential_type 00150 00151 TYPE all_potential_p_type 00152 TYPE(all_potential_type), POINTER :: all_potential 00153 END TYPE all_potential_p_type 00154 00155 TYPE gth_potential_p_type 00156 TYPE(gth_potential_type), POINTER :: gth_potential 00157 END TYPE gth_potential_p_type 00158 00159 ! ***************************************************************************** 00160 00161 ! Public subroutines 00162 PUBLIC :: allocate_potential,& 00163 deallocate_potential,& 00164 get_potential,& 00165 init_potential,& 00166 read_potential,& 00167 set_potential,& 00168 set_default_all_potential,& 00169 write_potential 00170 00171 ! Public data types 00172 00173 PUBLIC :: all_potential_type,& 00174 fist_potential_type,& 00175 gth_potential_type 00176 PUBLIC :: all_potential_p_type,& 00177 gth_potential_p_type 00178 00179 INTERFACE allocate_potential 00180 MODULE PROCEDURE allocate_all_potential,& 00181 allocate_fist_potential,& 00182 allocate_gth_potential 00183 END INTERFACE 00184 00185 INTERFACE deallocate_potential 00186 MODULE PROCEDURE deallocate_all_potential,& 00187 deallocate_fist_potential,& 00188 deallocate_gth_potential 00189 END INTERFACE 00190 00191 INTERFACE get_potential 00192 MODULE PROCEDURE get_all_potential,& 00193 get_fist_potential,& 00194 get_gth_potential 00195 END INTERFACE 00196 00197 INTERFACE init_potential 00198 MODULE PROCEDURE init_all_potential,& 00199 init_gth_potential 00200 END INTERFACE 00201 00202 INTERFACE read_potential 00203 MODULE PROCEDURE read_all_potential,& 00204 read_gth_potential 00205 END INTERFACE 00206 00207 INTERFACE set_potential 00208 MODULE PROCEDURE set_all_potential,& 00209 set_fist_potential,& 00210 set_gth_potential 00211 END INTERFACE 00212 00213 INTERFACE write_potential 00214 MODULE PROCEDURE write_all_potential,& 00215 write_gth_potential 00216 END INTERFACE 00217 00218 CONTAINS 00219 00220 ! ***************************************************************************** 00226 SUBROUTINE allocate_all_potential(potential,error) 00227 TYPE(all_potential_type), POINTER :: potential 00228 TYPE(cp_error_type), INTENT(inout) :: error 00229 00230 CHARACTER(len=*), PARAMETER :: routineN = 'allocate_all_potential', 00231 routineP = moduleN//':'//routineN 00232 00233 INTEGER :: stat 00234 LOGICAL :: failure 00235 00236 failure = .FALSE. 00237 IF (ASSOCIATED(potential)) CALL deallocate_potential(potential,error) 00238 00239 ALLOCATE (potential,STAT=stat) 00240 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00241 00242 NULLIFY (potential%elec_conf) 00243 00244 potential%description(1) = "All-electron potential" 00245 potential%description(2) = "Krack, Parrinello, PCCP 2, 2105 (2000)" 00246 00247 END SUBROUTINE allocate_all_potential 00248 00249 ! ***************************************************************************** 00254 SUBROUTINE allocate_fist_potential(potential,error) 00255 TYPE(fist_potential_type), POINTER :: potential 00256 TYPE(cp_error_type), INTENT(inout) :: error 00257 00258 CHARACTER(len=*), PARAMETER :: routineN = 'allocate_fist_potential', 00259 routineP = moduleN//':'//routineN 00260 00261 INTEGER :: stat 00262 LOGICAL :: failure 00263 00264 failure = .FALSE. 00265 IF (ASSOCIATED(potential)) CALL deallocate_potential(potential,error) 00266 00267 ALLOCATE (potential,STAT=stat) 00268 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00269 00270 potential%apol = 0.0_dp 00271 potential%cpol = 0.0_dp 00272 potential%mm_radius = 0.0_dp 00273 potential%qeff = 0.0_dp 00274 potential%qmmm_radius = 0.0_dp 00275 potential%qmmm_corr_radius = 0.0_dp 00276 00277 potential%description(1) = "Effective charge and inducible dipole potential" 00278 00279 END SUBROUTINE allocate_fist_potential 00280 00281 ! ***************************************************************************** 00287 SUBROUTINE allocate_gth_potential(potential,error) 00288 TYPE(gth_potential_type), POINTER :: potential 00289 TYPE(cp_error_type), INTENT(inout) :: error 00290 00291 CHARACTER(len=*), PARAMETER :: routineN = 'allocate_gth_potential', 00292 routineP = moduleN//':'//routineN 00293 00294 INTEGER :: stat 00295 LOGICAL :: failure 00296 00297 failure = .FALSE. 00298 IF (ASSOCIATED(potential)) CALL deallocate_potential(potential,error) 00299 00300 ALLOCATE (potential,STAT=stat) 00301 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00302 00303 NULLIFY (potential%alpha_ppnl) 00304 NULLIFY (potential%cexp_ppl) 00305 NULLIFY (potential%elec_conf) 00306 NULLIFY (potential%nprj_ppnl) 00307 NULLIFY (potential%cprj) 00308 NULLIFY (potential%cprj_ppnl) 00309 NULLIFY (potential%vprj_ppnl) 00310 NULLIFY (potential%hprj_ppnl) 00311 00312 NULLIFY (potential%alpha_lpot) 00313 NULLIFY (potential%nct_lpot) 00314 NULLIFY (potential%cval_lpot) 00315 NULLIFY (potential%alpha_lsd) 00316 NULLIFY (potential%nct_lsd) 00317 NULLIFY (potential%cval_lsd) 00318 NULLIFY (potential%alpha_nlcc) 00319 NULLIFY (potential%nct_nlcc) 00320 NULLIFY (potential%cval_nlcc) 00321 00322 potential%description(1) = "Goedecker-Teter-Hutter pseudopotential" 00323 potential%description(2) = "Goedecker et al., PRB 54, 1703 (1996)" 00324 potential%description(3) = "Hartwigsen et al., PRB 58, 3641 (1998)" 00325 potential%description(4) = "Krack, TCA 114, 145 (2005)" 00326 00327 END SUBROUTINE allocate_gth_potential 00328 00329 ! ***************************************************************************** 00335 SUBROUTINE deallocate_all_potential(potential,error) 00336 TYPE(all_potential_type), POINTER :: potential 00337 TYPE(cp_error_type), INTENT(inout) :: error 00338 00339 CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_all_potential', 00340 routineP = moduleN//':'//routineN 00341 00342 INTEGER :: stat 00343 LOGICAL :: failure 00344 00345 failure = .FALSE. 00346 IF (ASSOCIATED(potential)) THEN 00347 DEALLOCATE (potential%elec_conf,STAT=stat) 00348 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00349 DEALLOCATE (potential,STAT=stat) 00350 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00351 ELSE 00352 CALL stop_program(routineN,moduleN,__LINE__,& 00353 "The pointer potential is not associated.") 00354 END IF 00355 00356 END SUBROUTINE deallocate_all_potential 00357 00358 ! ***************************************************************************** 00363 SUBROUTINE deallocate_fist_potential(potential,error) 00364 TYPE(fist_potential_type), POINTER :: potential 00365 TYPE(cp_error_type), INTENT(inout) :: error 00366 00367 CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_fist_potential', 00368 routineP = moduleN//':'//routineN 00369 00370 INTEGER :: stat 00371 LOGICAL :: failure 00372 00373 failure = .FALSE. 00374 IF (ASSOCIATED(potential)) THEN 00375 ! Nothing exciting here yet. 00376 DEALLOCATE (potential,STAT=stat) 00377 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00378 ELSE 00379 CALL stop_program(routineN,moduleN,__LINE__,& 00380 "The pointer potential is not associated.") 00381 END IF 00382 00383 END SUBROUTINE deallocate_fist_potential 00384 00385 ! ***************************************************************************** 00391 SUBROUTINE deallocate_gth_potential(potential,error) 00392 TYPE(gth_potential_type), POINTER :: potential 00393 TYPE(cp_error_type), INTENT(inout) :: error 00394 00395 CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_gth_potential', 00396 routineP = moduleN//':'//routineN 00397 00398 INTEGER :: stat 00399 LOGICAL :: failure 00400 00401 failure = .FALSE. 00402 00403 IF (ASSOCIATED(potential)) THEN 00404 00405 DEALLOCATE (potential%elec_conf,STAT=stat) 00406 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00407 ! *** Deallocate the parameters of the local part *** 00408 00409 IF (ASSOCIATED(potential%cexp_ppl)) THEN 00410 DEALLOCATE (potential%cexp_ppl,STAT=stat) 00411 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00412 END IF 00413 00414 ! *** Deallocate the parameters of the non-local part *** 00415 IF (ASSOCIATED(potential%alpha_ppnl)) THEN 00416 DEALLOCATE (potential%alpha_ppnl,STAT=stat) 00417 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00418 DEALLOCATE (potential%cprj,STAT=stat) 00419 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00420 DEALLOCATE (potential%cprj_ppnl,STAT=stat) 00421 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00422 DEALLOCATE (potential%hprj_ppnl,STAT=stat) 00423 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00424 DEALLOCATE (potential%nprj_ppnl,STAT=stat) 00425 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00426 DEALLOCATE (potential%vprj_ppnl,STAT=stat) 00427 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00428 END IF 00429 00430 IF (ASSOCIATED(potential%alpha_lpot)) THEN 00431 DEALLOCATE (potential%alpha_lpot,STAT=stat) 00432 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00433 DEALLOCATE (potential%nct_lpot,STAT=stat) 00434 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00435 DEALLOCATE (potential%cval_lpot,STAT=stat) 00436 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00437 END IF 00438 00439 IF (ASSOCIATED(potential%alpha_lsd)) THEN 00440 DEALLOCATE (potential%alpha_lsd,STAT=stat) 00441 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00442 DEALLOCATE (potential%nct_lsd,STAT=stat) 00443 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00444 DEALLOCATE (potential%cval_lsd,STAT=stat) 00445 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00446 END IF 00447 00448 IF (ASSOCIATED(potential%alpha_nlcc)) THEN 00449 DEALLOCATE (potential%alpha_nlcc,STAT=stat) 00450 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00451 DEALLOCATE (potential%nct_nlcc,STAT=stat) 00452 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00453 DEALLOCATE (potential%cval_nlcc,STAT=stat) 00454 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00455 END IF 00456 00457 DEALLOCATE (potential,STAT=stat) 00458 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00459 ELSE 00460 00461 CALL stop_program(routineN,moduleN,__LINE__,& 00462 "The pointer potential is not associated.") 00463 END IF 00464 00465 END SUBROUTINE deallocate_gth_potential 00466 00467 ! ***************************************************************************** 00473 SUBROUTINE get_all_potential(potential,name,alpha_core_charge,& 00474 ccore_charge,core_charge_radius,z,zeff,& 00475 zeff_correction, elec_conf) 00476 TYPE(all_potential_type), POINTER :: potential 00477 CHARACTER(LEN=default_string_length), 00478 INTENT(OUT), OPTIONAL :: name 00479 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: alpha_core_charge, 00480 ccore_charge, 00481 core_charge_radius 00482 INTEGER, INTENT(OUT), OPTIONAL :: z 00483 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: zeff, zeff_correction 00484 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: elec_conf 00485 00486 CHARACTER(len=*), PARAMETER :: routineN = 'get_all_potential', 00487 routineP = moduleN//':'//routineN 00488 00489 IF (ASSOCIATED(potential)) THEN 00490 00491 IF (PRESENT(name)) name = potential%name 00492 IF (PRESENT(alpha_core_charge))& 00493 alpha_core_charge = potential%alpha_core_charge 00494 IF (PRESENT(ccore_charge)) ccore_charge = potential%ccore_charge 00495 IF (PRESENT(core_charge_radius))& 00496 core_charge_radius = potential%core_charge_radius 00497 IF (PRESENT(z)) z = potential%z 00498 IF (PRESENT(zeff)) zeff = potential%zeff 00499 IF (PRESENT(zeff_correction)) zeff_correction = potential%zeff_correction 00500 IF (PRESENT(elec_conf)) elec_conf => potential%elec_conf 00501 00502 ELSE 00503 00504 CALL stop_program(routineN,moduleN,__LINE__,& 00505 "The pointer potential is not associated.") 00506 00507 END IF 00508 00509 END SUBROUTINE get_all_potential 00510 00511 ! ***************************************************************************** 00517 SUBROUTINE get_fist_potential(potential, name, apol, cpol, mm_radius, qeff, & 00518 qmmm_corr_radius, qmmm_radius) 00519 TYPE(fist_potential_type), POINTER :: potential 00520 CHARACTER(LEN=default_string_length), 00521 INTENT(OUT), OPTIONAL :: name 00522 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: apol, cpol, mm_radius, qeff, 00523 qmmm_corr_radius, qmmm_radius 00524 00525 CHARACTER(len=*), PARAMETER :: routineN = 'get_fist_potential', 00526 routineP = moduleN//':'//routineN 00527 00528 IF (ASSOCIATED(potential)) THEN 00529 00530 IF (PRESENT(name)) name = potential%name 00531 IF (PRESENT(apol)) apol = potential%apol 00532 IF (PRESENT(cpol)) cpol = potential%cpol 00533 IF (PRESENT(mm_radius)) mm_radius = potential%mm_radius 00534 IF (PRESENT(qeff)) qeff = potential%qeff 00535 IF (PRESENT(qmmm_corr_radius)) qmmm_corr_radius = potential%qmmm_corr_radius 00536 IF (PRESENT(qmmm_radius)) qmmm_radius = potential%qmmm_radius 00537 00538 ELSE 00539 00540 CALL stop_program(routineN,moduleN,__LINE__,& 00541 "The pointer potential is not associated.") 00542 00543 END IF 00544 00545 END SUBROUTINE get_fist_potential 00546 00547 ! ***************************************************************************** 00553 SUBROUTINE get_gth_potential(potential,name,alpha_core_charge,& 00554 alpha_ppl,ccore_charge,cerf_ppl,& 00555 core_charge_radius,ppl_radius,ppnl_radius,& 00556 lppnl,lprj_ppnl_max,nexp_ppl,nppnl,& 00557 nprj_ppnl_max,z,zeff,zeff_correction,& 00558 ppl_present,ppnl_present,& 00559 alpha_ppnl,cexp_ppl,elec_conf,nprj_ppnl,cprj,& 00560 cprj_ppnl,vprj_ppnl,hprj_ppnl,& 00561 lpot_present,nexp_lpot,alpha_lpot,nct_lpot,cval_lpot,& 00562 lsd_present,nexp_lsd,alpha_lsd,nct_lsd,cval_lsd,& 00563 nlcc_present,nexp_nlcc,alpha_nlcc,nct_nlcc,cval_nlcc) 00564 00565 TYPE(gth_potential_type), POINTER :: potential 00566 CHARACTER(LEN=default_string_length), 00567 INTENT(OUT), OPTIONAL :: name 00568 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: alpha_core_charge, alpha_ppl, 00569 ccore_charge, cerf_ppl, core_charge_radius, ppl_radius, ppnl_radius 00570 INTEGER, INTENT(OUT), OPTIONAL :: lppnl, lprj_ppnl_max, 00571 nexp_ppl, nppnl, 00572 nprj_ppnl_max, z 00573 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: zeff, zeff_correction 00574 LOGICAL, INTENT(OUT), OPTIONAL :: ppl_present, ppnl_present 00575 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00576 POINTER :: alpha_ppnl, cexp_ppl 00577 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: elec_conf, nprj_ppnl 00578 REAL(KIND=dp), DIMENSION(:, :), 00579 OPTIONAL, POINTER :: cprj, cprj_ppnl, vprj_ppnl 00580 REAL(KIND=dp), DIMENSION(:, :, :), 00581 OPTIONAL, POINTER :: hprj_ppnl 00582 LOGICAL, INTENT(OUT), OPTIONAL :: lpot_present 00583 INTEGER, INTENT(OUT), OPTIONAL :: nexp_lpot 00584 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00585 POINTER :: alpha_lpot 00586 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: nct_lpot 00587 REAL(KIND=dp), DIMENSION(:, :), 00588 OPTIONAL, POINTER :: cval_lpot 00589 LOGICAL, INTENT(OUT), OPTIONAL :: lsd_present 00590 INTEGER, INTENT(OUT), OPTIONAL :: nexp_lsd 00591 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00592 POINTER :: alpha_lsd 00593 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: nct_lsd 00594 REAL(KIND=dp), DIMENSION(:, :), 00595 OPTIONAL, POINTER :: cval_lsd 00596 LOGICAL, INTENT(OUT), OPTIONAL :: nlcc_present 00597 INTEGER, INTENT(OUT), OPTIONAL :: nexp_nlcc 00598 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00599 POINTER :: alpha_nlcc 00600 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: nct_nlcc 00601 REAL(KIND=dp), DIMENSION(:, :), 00602 OPTIONAL, POINTER :: cval_nlcc 00603 00604 CHARACTER(len=*), PARAMETER :: routineN = 'get_gth_potential', 00605 routineP = moduleN//':'//routineN 00606 00607 IF (ASSOCIATED(potential)) THEN 00608 00609 IF (PRESENT(name)) name = potential%name 00610 IF (PRESENT(alpha_core_charge))& 00611 alpha_core_charge = potential%alpha_core_charge 00612 IF (PRESENT(alpha_ppl)) alpha_ppl = potential%alpha_ppl 00613 IF (PRESENT(ccore_charge)) ccore_charge = potential%ccore_charge 00614 IF (PRESENT(cerf_ppl)) cerf_ppl = potential%cerf_ppl 00615 IF (PRESENT(core_charge_radius))& 00616 core_charge_radius = potential%core_charge_radius 00617 IF (PRESENT(ppl_radius)) ppl_radius = potential%ppl_radius 00618 IF (PRESENT(ppnl_radius)) ppnl_radius = potential%ppnl_radius 00619 IF (PRESENT(lppnl)) lppnl = potential%lppnl 00620 IF (PRESENT(lprj_ppnl_max)) lprj_ppnl_max = potential%lprj_ppnl_max 00621 IF (PRESENT(nexp_ppl)) nexp_ppl = potential%nexp_ppl 00622 IF (PRESENT(nppnl)) nppnl = potential%nppnl 00623 IF (PRESENT(nprj_ppnl_max)) nprj_ppnl_max = potential%nprj_ppnl_max 00624 IF (PRESENT(z)) z = potential%z 00625 IF (PRESENT(zeff)) zeff = potential%zeff 00626 IF (PRESENT(zeff_correction)) zeff_correction = potential%zeff_correction 00627 IF (PRESENT(ppl_present)) ppl_present = (potential%nexp_ppl > 0) 00628 IF (PRESENT(ppnl_present)) ppnl_present = (potential%nppnl > 0) 00629 IF (PRESENT(alpha_ppnl)) alpha_ppnl => potential%alpha_ppnl 00630 IF (PRESENT(cexp_ppl)) cexp_ppl => potential%cexp_ppl 00631 IF (PRESENT(elec_conf)) elec_conf => potential%elec_conf 00632 IF (PRESENT(nprj_ppnl)) nprj_ppnl => potential%nprj_ppnl 00633 IF (PRESENT(cprj)) cprj => potential%cprj 00634 IF (PRESENT(cprj_ppnl)) cprj_ppnl => potential%cprj_ppnl 00635 IF (PRESENT(vprj_ppnl)) vprj_ppnl => potential%vprj_ppnl 00636 IF (PRESENT(hprj_ppnl)) hprj_ppnl => potential%hprj_ppnl 00637 00638 IF (PRESENT(lpot_present)) lpot_present = potential%lpotextended 00639 IF (PRESENT(nexp_lpot)) nexp_lpot = potential%nexp_lpot 00640 IF (PRESENT(alpha_lpot)) alpha_lpot => potential%alpha_lpot 00641 IF (PRESENT(nct_lpot)) nct_lpot => potential%nct_lpot 00642 IF (PRESENT(cval_lpot)) cval_lpot => potential%cval_lpot 00643 00644 IF (PRESENT(lsd_present)) lsd_present = potential%lsdpot 00645 IF (PRESENT(nexp_lsd)) nexp_lsd = potential%nexp_lsd 00646 IF (PRESENT(alpha_lsd)) alpha_lsd => potential%alpha_lsd 00647 IF (PRESENT(nct_lsd)) nct_lsd => potential%nct_lsd 00648 IF (PRESENT(cval_lsd)) cval_lsd => potential%cval_lsd 00649 00650 IF (PRESENT(nlcc_present)) nlcc_present = potential%nlcc 00651 IF (PRESENT(nexp_nlcc)) nexp_nlcc = potential%nexp_nlcc 00652 IF (PRESENT(alpha_nlcc)) alpha_nlcc => potential%alpha_nlcc 00653 IF (PRESENT(nct_nlcc)) nct_nlcc => potential%nct_nlcc 00654 IF (PRESENT(cval_nlcc)) cval_nlcc => potential%cval_nlcc 00655 00656 ELSE 00657 00658 CALL stop_program(routineN,moduleN,__LINE__,& 00659 "The pointer potential is not associated.") 00660 00661 END IF 00662 00663 END SUBROUTINE get_gth_potential 00664 00665 ! ***************************************************************************** 00674 SUBROUTINE init_cprj_ppnl(potential,error) 00675 TYPE(gth_potential_type), POINTER :: potential 00676 TYPE(cp_error_type), INTENT(inout) :: error 00677 00678 CHARACTER(len=*), PARAMETER :: routineN = 'init_cprj_ppnl', 00679 routineP = moduleN//':'//routineN 00680 00681 INTEGER :: cpx, cpy, cpz, cx, cy, cz, 00682 ico, iprj, iprj_ppnl, l, lp, 00683 lprj_ppnl, nprj, px, py, pz 00684 REAL(KIND=dp) :: alpha_ppnl, cp 00685 00686 nprj = 0 00687 00688 DO l=0,potential%lppnl 00689 alpha_ppnl = potential%alpha_ppnl(l) 00690 DO iprj_ppnl=1,potential%nprj_ppnl(l) 00691 lp = iprj_ppnl - 1 00692 lprj_ppnl = l + 2*lp 00693 cp = SQRT(2.0_dp**(2.0_dp*REAL(lprj_ppnl,dp) + 3.5_dp)* 00694 alpha_ppnl**(REAL(lprj_ppnl,dp) + 1.5_dp)/ 00695 (rootpi*dfac(2*lprj_ppnl + 1))) 00696 potential%cprj_ppnl(iprj_ppnl,l) = cp 00697 DO cx=0,l 00698 DO cy=0,l-cx 00699 cz = l - cx - cy 00700 iprj = nprj + co(cx,cy,cz) 00701 DO px=0,lp 00702 DO py=0,lp-px 00703 pz = lp - px - py 00704 cpx = cx + 2*px 00705 cpy = cy + 2*py 00706 cpz = cz + 2*pz 00707 ico = coset(cpx,cpy,cpz) 00708 potential%cprj(ico,iprj) = cp*fac(lp)/(fac(px)*fac(py)*fac(pz)) 00709 END DO 00710 END DO 00711 END DO 00712 END DO 00713 nprj = nprj + nco(l) 00714 END DO 00715 END DO 00716 00717 END SUBROUTINE init_cprj_ppnl 00718 00719 ! ***************************************************************************** 00725 SUBROUTINE init_gth_potential(potential,error) 00726 TYPE(gth_potential_type), POINTER :: potential 00727 TYPE(cp_error_type), INTENT(inout) :: error 00728 00729 CHARACTER(len=*), PARAMETER :: routineN = 'init_gth_potential', 00730 routineP = moduleN//':'//routineN 00731 00732 INTEGER :: handle 00733 00734 IF (.NOT.ASSOCIATED(potential)) RETURN 00735 00736 CALL timeset(routineN,handle) 00737 00738 IF (potential%nppnl > 0) THEN 00739 00740 ! *** Initialise the projector coefficients of the *** 00741 ! *** non-local part of the GTH pseudopotential and *** 00742 ! *** the transformation matrices "pgf" -> "prj_ppnl" *** 00743 00744 CALL init_cprj_ppnl(potential,error) 00745 00746 ! *** Initialise the h(i,j) projector coefficients of *** 00747 ! *** the non-local part of the GTH pseudopotential *** 00748 00749 CALL init_vprj_ppnl(potential,error) 00750 00751 END IF 00752 00753 CALL timestop(handle) 00754 00755 END SUBROUTINE init_gth_potential 00756 00757 ! ***************************************************************************** 00764 SUBROUTINE init_vprj_ppnl(potential,error) 00765 TYPE(gth_potential_type), POINTER :: potential 00766 TYPE(cp_error_type), INTENT(inout) :: error 00767 00768 CHARACTER(len=*), PARAMETER :: routineN = 'init_vprj_ppnl', 00769 routineP = moduleN//':'//routineN 00770 00771 INTEGER :: i, ico, iprj, iprj_ppnl, iso, 00772 j, jco, jprj, jprj_ppnl, l, 00773 nprj 00774 00775 nprj = 0 00776 00777 DO l=0,potential%lppnl 00778 DO iprj_ppnl=1,potential%nprj_ppnl(l) 00779 iprj = nprj + (iprj_ppnl - 1)*nco(l) 00780 DO jprj_ppnl=1,potential%nprj_ppnl(l) 00781 jprj = nprj + (jprj_ppnl - 1)*nco(l) 00782 DO ico=1,nco(l) 00783 i = iprj + ico 00784 DO jco=1,nco(l) 00785 j = jprj + jco 00786 DO iso=1,nso(l) 00787 potential%vprj_ppnl(i,j) = potential%vprj_ppnl(i,j) +& 00788 orbtramat(l)%slm(iso,ico)*& 00789 potential%hprj_ppnl(iprj_ppnl,& 00790 jprj_ppnl,l)*& 00791 orbtramat(l)%slm(iso,jco) 00792 END DO 00793 END DO 00794 END DO 00795 END DO 00796 END DO 00797 nprj = nprj + potential%nprj_ppnl(l)*nco(l) 00798 END DO 00799 00800 END SUBROUTINE init_vprj_ppnl 00801 00802 ! ***************************************************************************** 00803 SUBROUTINE init_all_potential(potential,itype,zeff,zeff_correction,error) 00804 00805 TYPE(all_potential_type), POINTER :: potential 00806 CHARACTER(LEN=*), OPTIONAL :: itype 00807 REAL(KIND=dp), OPTIONAL :: zeff, zeff_correction 00808 TYPE(cp_error_type), INTENT(inout) :: error 00809 00810 CHARACTER(len=*), PARAMETER :: routineN = 'init_all_potential', 00811 routineP = moduleN//':'//routineN 00812 00813 INTEGER :: dz, handle 00814 00815 IF (.NOT.ASSOCIATED(potential)) RETURN 00816 00817 CALL timeset(routineN,handle) 00818 00819 IF ( PRESENT (zeff) ) potential%zeff = zeff 00820 IF ( PRESENT (zeff_correction) ) potential%zeff_correction = zeff_correction 00821 dz = potential%z - INT(potential%zeff - potential%zeff_correction) 00822 SELECT CASE (dz) 00823 CASE DEFAULT 00824 CASE (2) 00825 potential%elec_conf(0) = potential%elec_conf(0) - 2 00826 CASE (10) 00827 potential%elec_conf(0) = potential%elec_conf(0) - 4 00828 potential%elec_conf(1) = potential%elec_conf(1) - 6 00829 CASE (18) 00830 potential%elec_conf(0) = potential%elec_conf(0) - 6 00831 potential%elec_conf(1) = potential%elec_conf(1) - 12 00832 CASE (28) 00833 potential%elec_conf(0) = potential%elec_conf(0) - 6 00834 potential%elec_conf(1) = potential%elec_conf(1) - 12 00835 potential%elec_conf(2) = potential%elec_conf(2) - 10 00836 CASE (30) 00837 potential%elec_conf(0) = potential%elec_conf(0) - 8 00838 potential%elec_conf(1) = potential%elec_conf(1) - 12 00839 potential%elec_conf(2) = potential%elec_conf(2) - 10 00840 CASE (36) 00841 potential%elec_conf(0) = potential%elec_conf(0) - 8 00842 potential%elec_conf(1) = potential%elec_conf(1) - 18 00843 potential%elec_conf(2) = potential%elec_conf(2) - 10 00844 CASE (46) 00845 potential%elec_conf(0) = potential%elec_conf(0) - 8 00846 potential%elec_conf(1) = potential%elec_conf(1) - 18 00847 potential%elec_conf(2) = potential%elec_conf(2) - 20 00848 CASE (48) 00849 potential%elec_conf(0) = potential%elec_conf(0) - 10 00850 potential%elec_conf(1) = potential%elec_conf(1) - 18 00851 potential%elec_conf(2) = potential%elec_conf(2) - 20 00852 CASE (54) 00853 potential%elec_conf(0) = potential%elec_conf(0) - 10 00854 potential%elec_conf(1) = potential%elec_conf(1) - 24 00855 potential%elec_conf(2) = potential%elec_conf(2) - 20 00856 CASE (68) 00857 potential%elec_conf(0) = potential%elec_conf(0) - 10 00858 potential%elec_conf(1) = potential%elec_conf(1) - 24 00859 potential%elec_conf(2) = potential%elec_conf(2) - 20 00860 potential%elec_conf(3) = potential%elec_conf(3) - 14 00861 CASE (78) 00862 potential%elec_conf(0) = potential%elec_conf(0) - 10 00863 potential%elec_conf(1) = potential%elec_conf(1) - 24 00864 potential%elec_conf(2) = potential%elec_conf(2) - 30 00865 potential%elec_conf(3) = potential%elec_conf(3) - 14 00866 CASE (80) 00867 potential%elec_conf(0) = potential%elec_conf(0) - 12 00868 potential%elec_conf(1) = potential%elec_conf(1) - 24 00869 potential%elec_conf(2) = potential%elec_conf(2) - 30 00870 potential%elec_conf(3) = potential%elec_conf(3) - 14 00871 CASE (86) 00872 potential%elec_conf(0) = potential%elec_conf(0) - 12 00873 potential%elec_conf(1) = potential%elec_conf(1) - 30 00874 potential%elec_conf(2) = potential%elec_conf(2) - 30 00875 potential%elec_conf(3) = potential%elec_conf(3) - 14 00876 CASE (100) 00877 potential%elec_conf(0) = potential%elec_conf(0) - 12 00878 potential%elec_conf(1) = potential%elec_conf(1) - 30 00879 potential%elec_conf(2) = potential%elec_conf(2) - 30 00880 potential%elec_conf(3) = potential%elec_conf(3) - 28 00881 END SELECT 00882 00883 IF ( PRESENT (itype) ) THEN 00884 IF ( itype == "BARE" ) THEN 00885 potential%description(1) = "Bare Coulomb Potential" 00886 IF ( dz > 0 ) THEN 00887 potential%description(2) = "Valence charge only" 00888 ELSE 00889 potential%description(2) = "Full atomic charge" 00890 END IF 00891 ENDIF 00892 END IF 00893 00894 CALL timestop(handle) 00895 00896 END SUBROUTINE init_all_potential 00897 00898 ! ***************************************************************************** 00904 SUBROUTINE read_all_potential(element_symbol,potential_name,potential,zeff_correction,& 00905 para_env, dft_section, potential_section, error) 00906 00907 CHARACTER(LEN=*), INTENT(IN) :: element_symbol, potential_name 00908 TYPE(all_potential_type), POINTER :: potential 00909 REAL(KIND=dp) :: zeff_correction 00910 TYPE(cp_para_env_type), POINTER :: para_env 00911 TYPE(section_vals_type), POINTER :: dft_section, potential_section 00912 TYPE(cp_error_type), INTENT(inout) :: error 00913 00914 CHARACTER(len=*), PARAMETER :: routineN = 'read_all_potential', 00915 routineP = moduleN//':'//routineN 00916 00917 CHARACTER(LEN=240) :: line 00918 CHARACTER(LEN=242) :: line2 00919 CHARACTER(len=5*default_string_length) :: line_att 00920 CHARACTER(len=default_path_length) :: potential_file_name 00921 CHARACTER(LEN=LEN(element_symbol)) :: symbol 00922 CHARACTER(LEN=LEN(element_symbol)+2) :: symbol2 00923 CHARACTER(LEN=LEN(potential_name)) :: apname 00924 CHARACTER(LEN=LEN(potential_name)+2) :: apname2 00925 INTEGER :: irep, l, stat, strlen1, 00926 strlen2 00927 INTEGER, DIMENSION(:), POINTER :: elec_conf 00928 LOGICAL :: failure, found, is_ok, match, 00929 read_from_input 00930 REAL(KIND=dp) :: alpha, r 00931 TYPE(cp_parser_type), POINTER :: parser 00932 TYPE(cp_sll_val_type), POINTER :: list 00933 TYPE(val_type), POINTER :: val 00934 00935 line2 = "" 00936 symbol2 = "" 00937 apname2 = "" 00938 NULLIFY(parser) 00939 failure = .FALSE. 00940 CALL cite_reference(Krack2000) 00941 00942 potential%name = potential_name 00943 read_from_input = .FALSE. 00944 CALL section_vals_get(potential_section,explicit=read_from_input, error=error) 00945 IF (.NOT.read_from_input) THEN 00946 CALL section_vals_val_get(dft_section,"POTENTIAL_FILE_NAME",c_val=potential_file_name,& 00947 error=error) 00948 CALL parser_create(parser,potential_file_name,para_env=para_env,error=error) 00949 END IF 00950 00951 ! *** Search for the requested potential in the potential file *** 00952 ! *** until the potential is found or the end of file is reached *** 00953 00954 apname = potential_name 00955 symbol = element_symbol 00956 irep = 0 00957 search_loop: DO 00958 IF (read_from_input) THEN 00959 NULLIFY(list, val) 00960 found = .TRUE. 00961 CALL section_vals_list_get(potential_section,"_DEFAULT_KEYWORD_",list=list,error=error) 00962 ELSE 00963 CALL parser_search_string(parser,TRIM(apname),.TRUE.,found,line,error=error) 00964 END IF 00965 IF (found) THEN 00966 CALL uppercase(symbol) 00967 CALL uppercase(apname) 00968 00969 IF (read_from_input) THEN 00970 match = .TRUE. 00971 ELSE 00972 ! Check both the element symbol and the atomic potential name 00973 match = .FALSE. 00974 CALL uppercase(line) 00975 line2 = " "//line//" " 00976 symbol2 = " "//TRIM(symbol)//" " 00977 apname2 = " "//TRIM(apname)//" " 00978 strlen1 = LEN_TRIM(symbol2) + 1 00979 strlen2 = LEN_TRIM(apname2) + 1 00980 00981 IF ( (INDEX(line2,symbol2(:strlen1)) > 0).AND.& 00982 (INDEX(line2,apname2(:strlen2)) > 0)) match = .TRUE. 00983 END IF 00984 IF (match) THEN 00985 ! Read the electronic configuration 00986 NULLIFY (elec_conf) 00987 l = 0 00988 CALL reallocate(elec_conf,0,l) 00989 IF (read_from_input) THEN 00990 is_ok=cp_sll_val_next(list,val,error=error) 00991 IF (.NOT.is_ok) CALL stop_program(routineN,moduleN, __LINE__,& 00992 "Error reading the Potential from input file!!") 00993 CALL val_get(val,c_val=line_att,error=error) 00994 READ(line_att,*)elec_conf(l) 00995 CALL remove_word(line_att) 00996 DO WHILE (LEN_TRIM(line_att) /= 0) 00997 l = l + 1 00998 CALL reallocate(elec_conf,0,l) 00999 READ(line_att,*)elec_conf(l) 01000 CALL remove_word(line_att) 01001 END DO 01002 ELSE 01003 CALL parser_get_object(parser,elec_conf(l),newline=.TRUE.,error=error) 01004 DO WHILE (parser_test_next_token(parser,error=error) == "INT") 01005 l = l + 1 01006 CALL reallocate(elec_conf,0,l) 01007 CALL parser_get_object(parser,elec_conf(l),error=error) 01008 END DO 01009 irep = irep + 1 01010 WRITE(line_att,'(100(1X,I0))')elec_conf 01011 CALL section_vals_val_set(potential_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,& 01012 c_val=TRIM(line_att), error=error) 01013 END IF 01014 01015 CALL reallocate(potential%elec_conf,0,l) 01016 potential%elec_conf(:) = elec_conf(:) 01017 01018 potential%zeff_correction = zeff_correction 01019 potential%zeff = REAL(SUM(elec_conf),dp)+zeff_correction 01020 01021 DEALLOCATE (elec_conf,STAT=stat) 01022 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01023 01024 ! Read r(loc) to define the exponent of the core charge 01025 ! distribution and calculate the corresponding coefficient 01026 01027 IF (read_from_input) THEN 01028 is_ok=cp_sll_val_next(list,val,error=error) 01029 IF (.NOT.is_ok) CALL stop_program(routineN,moduleN,__LINE__,& 01030 "Error reading the Potential from input file!!") 01031 CALL val_get(val,c_val=line_att,error=error) 01032 READ(line_att,*)r 01033 ELSE 01034 CALL parser_get_object(parser,r,newline=.TRUE.,error=error) 01035 irep = irep + 1 01036 WRITE(line_att,'(E24.16)')r 01037 CALL section_vals_val_set(potential_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,& 01038 c_val=TRIM(line_att), error=error) 01039 END IF 01040 alpha = 1.0_dp/(2.0_dp*r**2) 01041 01042 potential%alpha_core_charge = alpha 01043 potential%ccore_charge = potential%zeff*SQRT((alpha/pi)**3) 01044 01045 EXIT search_loop 01046 END IF 01047 ELSE 01048 ! Stop program, if the end of file is reached 01049 CALL stop_program(routineN,moduleN,__LINE__,& 01050 "The requested atomic potential <"//& 01051 TRIM(potential_name)//& 01052 "> was not found in the potential file <"//& 01053 TRIM(potential_file_name)//">") 01054 END IF 01055 END DO search_loop 01056 01057 IF (.NOT.read_from_input) THEN 01058 ! Dump the potential info the in potential section 01059 IF (match) THEN 01060 irep = irep + 1 01061 WRITE(line_att,'(A)')" # Potential name: "//apname2(:strlen2)//" for symbol: "//symbol2(:strlen1) 01062 CALL section_vals_val_set(potential_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,& 01063 c_val=TRIM(line_att), error=error) 01064 irep = irep + 1 01065 WRITE(line_att,'(A)')" # Potential read from the potential filename: "//TRIM(potential_file_name) 01066 CALL section_vals_val_set(potential_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,& 01067 c_val=TRIM(line_att), error=error) 01068 END IF 01069 CALL parser_release(parser,error=error) 01070 END IF 01071 END SUBROUTINE read_all_potential 01072 01073 ! ***************************************************************************** 01084 SUBROUTINE read_gth_potential(element_symbol,potential_name,potential,zeff_correction,& 01085 para_env, dft_section, potential_section, error) 01086 01087 CHARACTER(LEN=*), INTENT(IN) :: element_symbol, potential_name 01088 TYPE(gth_potential_type), POINTER :: potential 01089 REAL(KIND=dp), INTENT(IN) :: zeff_correction 01090 TYPE(cp_para_env_type), POINTER :: para_env 01091 TYPE(section_vals_type), POINTER :: dft_section, potential_section 01092 TYPE(cp_error_type), INTENT(inout) :: error 01093 01094 CHARACTER(len=*), PARAMETER :: routineN = 'read_gth_potential', 01095 routineP = moduleN//':'//routineN 01096 01097 CHARACTER(LEN=240) :: line 01098 CHARACTER(LEN=242) :: line2 01099 CHARACTER(len=5*default_string_length) :: line_att 01100 CHARACTER(len=default_path_length) :: potential_file_name 01101 CHARACTER(LEN=LEN(element_symbol)) :: symbol 01102 CHARACTER(LEN=LEN(element_symbol)+2) :: symbol2 01103 CHARACTER(LEN=LEN(potential_name)) :: apname 01104 CHARACTER(LEN=LEN(potential_name)+2) :: apname2 01105 INTEGER :: i, ic, ipot, irep, istr, j, l, lppnl, lprj_ppnl_max, maxlppl, 01106 n, nppnl, nprj_ppnl, nprj_ppnl_max, stat, strlen1, strlen2 01107 INTEGER, DIMENSION(:), POINTER :: elec_conf 01108 LOGICAL :: failure, found, is_ok, match, 01109 read_from_input 01110 REAL(KIND=dp) :: alpha, ci, r, rc2 01111 REAL(KIND=dp), DIMENSION(:), POINTER :: tmp_vals 01112 REAL(KIND=dp), DIMENSION(:, :, :), 01113 POINTER :: hprj_ppnl 01114 TYPE(cp_parser_type), POINTER :: parser 01115 TYPE(cp_sll_val_type), POINTER :: list 01116 TYPE(val_type), POINTER :: val 01117 01118 line2 = "" 01119 symbol2 = "" 01120 apname2 = "" 01121 NULLIFY(parser,tmp_vals) 01122 failure = .FALSE. 01123 CALL cite_reference(Goedecker1996) 01124 CALL cite_reference(Hartwigsen1998) 01125 CALL cite_reference(Krack2005) 01126 01127 potential%name = potential_name 01128 read_from_input = .FALSE. 01129 CALL section_vals_get(potential_section,explicit=read_from_input, error=error) 01130 IF (.NOT.read_from_input) THEN 01131 CALL section_vals_val_get(dft_section,"POTENTIAL_FILE_NAME",c_val=potential_file_name,& 01132 error=error) 01133 CALL parser_create(parser,potential_file_name,para_env=para_env,error=error) 01134 END IF 01135 01136 !initialize extended form 01137 potential%lpotextended = .FALSE. 01138 potential%nexp_lpot = 0 01139 potential%lsdpot = .FALSE. 01140 potential%nexp_lsd = 0 01141 potential%nlcc = .FALSE. 01142 potential%nexp_nlcc = 0 01143 01144 ! *** Search for the requested potential in the potential file *** 01145 ! *** until the potential is found or the end of file is reached *** 01146 01147 apname = potential_name 01148 symbol = element_symbol 01149 irep = 0 01150 search_loop: DO 01151 IF (read_from_input) THEN 01152 NULLIFY(list, val) 01153 found = .TRUE. 01154 CALL section_vals_list_get(potential_section,"_DEFAULT_KEYWORD_",list=list,error=error) 01155 ELSE 01156 CALL parser_search_string(parser,TRIM(apname),.TRUE.,found,line,error=error) 01157 END IF 01158 IF (found) THEN 01159 CALL uppercase(symbol) 01160 CALL uppercase(apname) 01161 IF (read_from_input) THEN 01162 match = .TRUE. 01163 ELSE 01164 ! *** Check both the element symbol and the atomic potential name *** 01165 match = .FALSE. 01166 CALL uppercase(line) 01167 line2 = " "//line//" " 01168 symbol2 = " "//TRIM(symbol)//" " 01169 apname2 = " "//TRIM(apname)//" " 01170 strlen1 = LEN_TRIM(symbol2) + 1 01171 strlen2 = LEN_TRIM(apname2) + 1 01172 01173 IF ( (INDEX(line2,symbol2(:strlen1)) > 0).AND.& 01174 (INDEX(line2,apname2(:strlen2)) > 0) ) match = .TRUE. 01175 END IF 01176 IF (match) THEN 01177 ! *** Read the electronic configuration *** 01178 NULLIFY (elec_conf) 01179 l = 0 01180 CALL reallocate(elec_conf,0,l) 01181 IF (read_from_input) THEN 01182 is_ok=cp_sll_val_next(list,val,error=error) 01183 IF (.NOT.is_ok) CALL stop_program(routineN,moduleN,__LINE__,& 01184 "Error reading the Potential from input file!!") 01185 CALL val_get(val,c_val=line_att,error=error) 01186 READ(line_att,*)elec_conf(l) 01187 CALL remove_word(line_att) 01188 DO WHILE (LEN_TRIM(line_att) /= 0) 01189 l = l + 1 01190 CALL reallocate(elec_conf,0,l) 01191 READ(line_att,*)elec_conf(l) 01192 CALL remove_word(line_att) 01193 END DO 01194 ELSE 01195 CALL parser_get_object(parser,elec_conf(l),newline=.TRUE.,error=error) 01196 DO WHILE (parser_test_next_token(parser,error=error) == "INT") 01197 l = l + 1 01198 CALL reallocate(elec_conf,0,l) 01199 CALL parser_get_object(parser,elec_conf(l),error=error) 01200 END DO 01201 irep = irep + 1 01202 WRITE(line_att,'(100(1X,I0))')elec_conf 01203 CALL section_vals_val_set(potential_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,& 01204 c_val=TRIM(line_att), error=error) 01205 END IF 01206 01207 CALL reallocate(potential%elec_conf,0,l) 01208 potential%elec_conf(:) = elec_conf(:) 01209 01210 potential%zeff_correction = zeff_correction 01211 potential%zeff = REAL(SUM(elec_conf),dp)+zeff_correction 01212 01213 DEALLOCATE (elec_conf,STAT=stat) 01214 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01215 01216 ! *** Read r(loc) to define the exponent of the core charge *** 01217 ! *** distribution and calculate the corresponding coefficient *** 01218 01219 IF (read_from_input) THEN 01220 is_ok=cp_sll_val_next(list,val,error=error) 01221 IF (.NOT.is_ok) CALL stop_program(routineN,moduleN,__LINE__,& 01222 "Error reading the Potential from input file!!") 01223 CALL val_get(val,c_val=line_att,error=error) 01224 READ(line_att,*)r 01225 CALL remove_word(line_att) 01226 ELSE 01227 line_att = "" 01228 CALL parser_get_object(parser,r,newline=.TRUE.,error=error) 01229 istr=LEN_TRIM(line_att)+1 01230 WRITE(line_att(istr:),'(E24.16)')r 01231 END IF 01232 alpha = 1.0_dp/(2.0_dp*r**2) 01233 01234 potential%alpha_core_charge = alpha 01235 potential%ccore_charge = potential%zeff*SQRT((alpha/pi)**3) 01236 01237 potential%alpha_ppl = alpha 01238 potential%cerf_ppl = potential%zeff*SQRT((alpha/pi)**3) 01239 01240 ! *** Read the parameters for the local part *** 01241 ! *** of the GTH pseudopotential (ppl) *** 01242 01243 IF (read_from_input) THEN 01244 READ(line_att,*)n 01245 CALL remove_word(line_att) 01246 ELSE 01247 CALL parser_get_object(parser,n,error=error) 01248 istr=LEN_TRIM(line_att)+1 01249 WRITE(line_att(istr:),'(1X,I0)')n 01250 END IF 01251 potential%nexp_ppl = n 01252 CALL reallocate(potential%cexp_ppl,1,n) 01253 01254 DO i=1,n 01255 IF (read_from_input) THEN 01256 READ(line_att,*) ci 01257 CALL remove_word(line_att) 01258 ELSE 01259 CALL parser_get_object(parser,ci,error=error) 01260 istr=LEN_TRIM(line_att)+1 01261 WRITE(line_att(istr:),'(E24.16)') ci 01262 END IF 01263 rc2 = (2.0_dp*potential%alpha_ppl) 01264 potential%cexp_ppl(i) = rc2**(i - 1)*ci 01265 END DO 01266 01267 IF (.NOT.read_from_input) THEN 01268 irep = irep + 1 01269 CALL section_vals_val_set(potential_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,& 01270 c_val=TRIM(line_att), error=error) 01271 line_att = "" 01272 ELSE 01273 IF (LEN_TRIM(line_att)/= 0) THEN 01274 CALL stop_program(routineN,moduleN,__LINE__,& 01275 "Error reading the Potential from input file!!") 01276 END IF 01277 END IF 01278 maxlppl = 2*(n - 1) 01279 01280 IF (maxlppl > -1) CALL init_orbital_pointers(maxlppl) 01281 01282 ! Read extended form of GTH pseudopotential 01283 ! local potential, NLCC, LSD potential 01284 IF (read_from_input) THEN 01285 DO 01286 is_ok=cp_sll_val_next(list,val,error=error) 01287 CPPostcondition(is_ok, cp_failure_level, routineP, error, failure) 01288 CALL val_get(val,c_val=line_att,error=error) 01289 IF(INDEX(line_att,"LPOT") /= 0) THEN 01290 potential%lpotextended = .TRUE. 01291 CALL remove_word(line_att) 01292 READ(line_att,*) potential%nexp_lpot 01293 n=potential%nexp_lpot 01294 maxlppl = 2*(n - 1) 01295 IF (maxlppl > -1) CALL init_orbital_pointers(maxlppl) 01296 NULLIFY(potential%alpha_lpot,potential%nct_lpot,potential%cval_lpot) 01297 CALL reallocate(potential%alpha_lpot,1,n) 01298 CALL reallocate(potential%nct_lpot,1,n) 01299 CALL reallocate(potential%cval_lpot,1,4,1,n) 01300 DO ipot=1,potential%nexp_lpot 01301 is_ok=cp_sll_val_next(list,val,error=error) 01302 CPPostcondition(is_ok, cp_failure_level, routineP, error, failure) 01303 CALL val_get(val,c_val=line_att,error=error) 01304 READ(line_att,*) r 01305 potential%alpha_lpot(ipot) = 0.5_dp/(r*r) 01306 CALL remove_word(line_att) 01307 READ(line_att,*) potential%nct_lpot(ipot) 01308 CALL remove_word(line_att) 01309 DO ic=1,potential%nct_lpot(ipot) 01310 READ(line_att,*) ci 01311 rc2 = (2._dp*potential%alpha_lpot(ipot))**(ic-1) 01312 potential%cval_lpot(ic,ipot) = ci*rc2 01313 CALL remove_word(line_att) 01314 END DO 01315 END DO 01316 ELSEIF(INDEX(line_att,"NLCC") /= 0) THEN 01317 potential%nlcc = .TRUE. 01318 CALL remove_word(line_att) 01319 READ(line_att,*) potential%nexp_nlcc 01320 n=potential%nexp_nlcc 01321 NULLIFY(potential%alpha_nlcc,potential%nct_nlcc,potential%cval_nlcc) 01322 CALL reallocate(potential%alpha_nlcc,1,n) 01323 CALL reallocate(potential%nct_nlcc,1,n) 01324 CALL reallocate(potential%cval_nlcc,1,4,1,n) 01325 DO ipot=1,potential%nexp_nlcc 01326 is_ok=cp_sll_val_next(list,val,error=error) 01327 CPPostcondition(is_ok, cp_failure_level, routineP, error, failure) 01328 CALL val_get(val,c_val=line_att,error=error) 01329 READ(line_att,*) potential%alpha_nlcc(ipot) 01330 CALL remove_word(line_att) 01331 READ(line_att,*) potential%nct_nlcc(ipot) 01332 CALL remove_word(line_att) 01333 DO ic=1,potential%nct_nlcc(ipot) 01334 READ(line_att,*) potential%cval_nlcc(ic,ipot) 01335 !make it compatible with bigdft style 01336 potential%cval_nlcc(ic,ipot)=potential%cval_nlcc(ic,ipot)/(4.0_dp*pi) 01337 CALL remove_word(line_att) 01338 END DO 01339 END DO 01340 ELSEIF(INDEX(line_att,"LSD") /= 0) THEN 01341 potential%lsdpot = .TRUE. 01342 CALL remove_word(line_att) 01343 READ(line_att,*) potential%nexp_lsd 01344 n=potential%nexp_lsd 01345 NULLIFY(potential%alpha_lsd,potential%nct_lsd,potential%cval_lsd) 01346 CALL reallocate(potential%alpha_lsd,1,n) 01347 CALL reallocate(potential%nct_lsd,1,n) 01348 CALL reallocate(potential%cval_lsd,1,4,1,n) 01349 DO ipot=1,potential%nexp_lsd 01350 is_ok=cp_sll_val_next(list,val,error=error) 01351 CPPostcondition(is_ok, cp_failure_level, routineP, error, failure) 01352 CALL val_get(val,c_val=line_att,error=error) 01353 READ(line_att,*) r 01354 potential%alpha_lsd(ipot) = 0.5_dp/(r*r) 01355 CALL remove_word(line_att) 01356 READ(line_att,*) potential%nct_lsd(ipot) 01357 CALL remove_word(line_att) 01358 DO ic=1,potential%nct_lsd(ipot) 01359 READ(line_att,*) ci 01360 rc2 = (2._dp*potential%alpha_lsd(ipot))**(ic-1) 01361 potential%cval_lsd(ic,ipot) = ci*rc2 01362 CALL remove_word(line_att) 01363 END DO 01364 END DO 01365 ELSE 01366 EXIT 01367 END IF 01368 END DO 01369 ELSE 01370 DO 01371 CALL parser_get_next_line(parser,1,error=error) 01372 IF(parser_test_next_token(parser,error=error) == "INT") THEN 01373 EXIT 01374 ELSEIF(parser_test_next_token(parser,error=error) == "STR") THEN 01375 CALL parser_get_object(parser,line,error=error) 01376 IF(INDEX(LINE,"LPOT") /= 0) THEN 01377 ! local potential 01378 potential%lpotextended = .TRUE. 01379 CALL parser_get_object(parser,potential%nexp_lpot,error=error) 01380 n = potential%nexp_lpot 01381 NULLIFY(potential%alpha_lpot,potential%nct_lpot,potential%cval_lpot) 01382 CALL reallocate(potential%alpha_lpot,1,n) 01383 CALL reallocate(potential%nct_lpot,1,n) 01384 CALL reallocate(potential%cval_lpot,1,4,1,n) 01385 ! add to input section 01386 WRITE(line_att,'(A,1X,I0)') "LPOT",n 01387 irep = irep + 1 01388 CALL section_vals_val_set(potential_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,& 01389 c_val=TRIM(line_att), error=error) 01390 DO ipot=1,potential%nexp_lpot 01391 CALL parser_get_object(parser,r,newline=.TRUE.,error=error) 01392 potential%alpha_lpot(ipot) = 0.5_dp/(r*r) 01393 CALL parser_get_object(parser,potential%nct_lpot(ipot),error=error) 01394 CALL reallocate(tmp_vals,1,potential%nct_lpot(ipot)) 01395 DO ic=1,potential%nct_lpot(ipot) 01396 CALL parser_get_object(parser,ci,error=error) 01397 tmp_vals(ic)=ci 01398 rc2 = (2._dp*potential%alpha_lpot(ipot))**(ic-1) 01399 potential%cval_lpot(ic,ipot) = ci*rc2 01400 END DO 01401 ! add to input section 01402 irep = irep + 1 01403 WRITE(line_att,'(E24.16,1X,I0,100(1X,E24.16))') r,potential%nct_lpot(ipot), & 01404 tmp_vals(1:potential%nct_lpot(ipot)) 01405 CALL section_vals_val_set(potential_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,& 01406 c_val=TRIM(line_att), error=error) 01407 END DO 01408 ELSEIF(INDEX(LINE,"NLCC") /= 0) THEN 01409 ! NLCC 01410 potential%nlcc = .TRUE. 01411 CALL parser_get_object(parser,potential%nexp_nlcc,error=error) 01412 n = potential%nexp_nlcc 01413 NULLIFY(potential%alpha_nlcc,potential%nct_nlcc,potential%cval_nlcc) 01414 CALL reallocate(potential%alpha_nlcc,1,n) 01415 CALL reallocate(potential%nct_nlcc,1,n) 01416 CALL reallocate(potential%cval_nlcc,1,4,1,n) 01417 ! add to input section 01418 WRITE(line_att,'(A,1X,I0)') "NLCC",n 01419 irep = irep + 1 01420 CALL section_vals_val_set(potential_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,& 01421 c_val=TRIM(line_att), error=error) 01422 DO ipot=1,potential%nexp_nlcc 01423 CALL parser_get_object(parser,potential%alpha_nlcc(ipot),newline=.TRUE.,error=error) 01424 CALL parser_get_object(parser,potential%nct_nlcc(ipot),error=error) 01425 CALL reallocate(tmp_vals,1,potential%nct_nlcc(ipot)) 01426 DO ic=1,potential%nct_nlcc(ipot) 01427 CALL parser_get_object(parser,potential%cval_nlcc(ic,ipot),error=error) 01428 tmp_vals(ic)=potential%cval_nlcc(ic,ipot) 01429 !make it compatible with bigdft style 01430 potential%cval_nlcc(ic,ipot)=potential%cval_nlcc(ic,ipot)/(4.0_dp*pi) 01431 END DO 01432 ! add to input section 01433 irep = irep + 1 01434 WRITE(line_att,'(E24.16,1X,I0,100(1X,E24.16))') potential%alpha_nlcc(ipot), potential%nct_nlcc(ipot), & 01435 tmp_vals(1:potential%nct_nlcc(ipot)) 01436 CALL section_vals_val_set(potential_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,& 01437 c_val=TRIM(line_att), error=error) 01438 END DO 01439 ELSEIF(INDEX(LINE,"LSD") /= 0) THEN 01440 ! LSD potential 01441 potential%lsdpot = .TRUE. 01442 CALL parser_get_object(parser,potential%nexp_lsd,error=error) 01443 n = potential%nexp_lsd 01444 NULLIFY(potential%alpha_lsd,potential%nct_lsd,potential%cval_lsd) 01445 CALL reallocate(potential%alpha_lsd,1,n) 01446 CALL reallocate(potential%nct_lsd,1,n) 01447 CALL reallocate(potential%cval_lsd,1,4,1,n) 01448 ! add to input section 01449 WRITE(line_att,'(A,1X,I0)') "LSD",n 01450 irep = irep + 1 01451 CALL section_vals_val_set(potential_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,& 01452 c_val=TRIM(line_att), error=error) 01453 DO ipot=1,potential%nexp_lsd 01454 CALL parser_get_object(parser,r,newline=.TRUE.,error=error) 01455 potential%alpha_lsd(ipot) = 0.5_dp/(r*r) 01456 CALL parser_get_object(parser,potential%nct_lsd(ipot),error=error) 01457 CALL reallocate(tmp_vals,1,potential%nct_lsd(ipot)) 01458 DO ic=1,potential%nct_lsd(ipot) 01459 CALL parser_get_object(parser,ci,error=error) 01460 tmp_vals(ic)=ci 01461 rc2 = (2._dp*potential%alpha_lsd(ipot))**(ic-1) 01462 potential%cval_lsd(ic,ipot) = ci*rc2 01463 END DO 01464 ! add to input section 01465 irep = irep + 1 01466 WRITE(line_att,'(E24.16,1X,I0,100(1X,E24.16))') r,potential%nct_lsd(ipot), & 01467 tmp_vals(1:potential%nct_lsd(ipot)) 01468 CALL section_vals_val_set(potential_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,& 01469 c_val=TRIM(line_att), error=error) 01470 END DO 01471 ELSE 01472 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 01473 END IF 01474 ELSE 01475 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 01476 END IF 01477 END DO 01478 END IF 01479 ! *** Read the parameters for the non-local *** 01480 ! *** part of the GTH pseudopotential (ppnl) *** 01481 IF (read_from_input) THEN 01482 READ(line_att,*)n 01483 CALL remove_word(line_att) 01484 ELSE 01485 CALL parser_get_object(parser,n,error=error) 01486 irep = irep + 1 01487 WRITE(line_att,'(1X,I0)')n 01488 CALL section_vals_val_set(potential_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,& 01489 c_val=TRIM(line_att), error=error) 01490 END IF 01491 potential%lppnl = n - 1 01492 potential%nppnl = 0 01493 01494 potential%lprj_ppnl_max = n - 1 01495 potential%nprj_ppnl_max = 0 01496 01497 IF (n > 0) THEN 01498 01499 lppnl = potential%lppnl 01500 nppnl = potential%nppnl 01501 01502 CALL init_orbital_pointers(lppnl) 01503 01504 NULLIFY (hprj_ppnl) 01505 01506 ! Load the parameter for n non-local projectors 01507 01508 CALL reallocate(potential%alpha_ppnl,0,lppnl) 01509 CALL reallocate(potential%nprj_ppnl,0,lppnl) 01510 01511 lprj_ppnl_max = -1 01512 nprj_ppnl_max = 0 01513 01514 DO l=0,lppnl 01515 IF (read_from_input) THEN 01516 is_ok=cp_sll_val_next(list,val,error=error) 01517 IF (.NOT.is_ok) CALL stop_program(routineN,moduleN,__LINE__,& 01518 "Error reading the Potential from input file!!") 01519 CALL val_get(val,c_val=line_att,error=error) 01520 READ(line_att,*)r 01521 CALL remove_word(line_att) 01522 READ(line_att,*)nprj_ppnl 01523 CALL remove_word(line_att) 01524 ELSE 01525 line_att ="" 01526 CALL parser_get_object(parser,r,newline=.TRUE.,error=error) 01527 CALL parser_get_object(parser,nprj_ppnl,error=error) 01528 istr=LEN_TRIM(line_att)+1 01529 WRITE(line_att(istr:),'(E24.16,1X,I0)')r, nprj_ppnl 01530 END IF 01531 IF (r==0.0_dp.AND.nprj_ppnl/=0) THEN 01532 CALL stop_program(routineN,moduleN,__LINE__,& 01533 "An error was detected in the atomic potential <"//& 01534 TRIM(potential_name)//& 01535 "> potential file <"//& 01536 TRIM(potential_file_name)//">",para_env) 01537 END IF 01538 potential%alpha_ppnl(l) = 0.0_dp 01539 IF (r/=0.0_dp.AND.n/=0) potential%alpha_ppnl(l) = 1.0_dp/(2.0_dp*r**2) 01540 potential%nprj_ppnl(l) = nprj_ppnl 01541 nppnl = nppnl + nprj_ppnl*nco(l) 01542 IF (nprj_ppnl > nprj_ppnl_max) THEN 01543 nprj_ppnl_max = nprj_ppnl 01544 CALL reallocate(hprj_ppnl,1,nprj_ppnl_max,& 01545 1,nprj_ppnl_max,& 01546 0,lppnl) 01547 END IF 01548 DO i=1,nprj_ppnl 01549 IF (i == 1) THEN 01550 IF (read_from_input) THEN 01551 READ(line_att,*)hprj_ppnl(i,i,l) 01552 CALL remove_word(line_att) 01553 ELSE 01554 CALL parser_get_object(parser,hprj_ppnl(i,i,l),error=error) 01555 istr=LEN_TRIM(line_att)+1 01556 WRITE(line_att(istr:),'(E24.16)')hprj_ppnl(i,i,l) 01557 END IF 01558 ELSE 01559 IF (read_from_input) THEN 01560 IF (LEN_TRIM(line_att)/=0) CALL stop_program(routineN,moduleN,__LINE__,& 01561 "Error reading the Potential from input file!!") 01562 is_ok=cp_sll_val_next(list,val,error=error) 01563 IF (.NOT.is_ok) CALL stop_program(routineN,moduleN,__LINE__,& 01564 "Error reading the Potential from input file!!") 01565 CALL val_get(val,c_val=line_att,error=error) 01566 READ(line_att,*)hprj_ppnl(i,i,l) 01567 CALL remove_word(line_att) 01568 ELSE 01569 irep = irep + 1 01570 CALL section_vals_val_set(potential_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,& 01571 c_val=TRIM(line_att), error=error) 01572 line_att ="" 01573 CALL parser_get_object(parser,hprj_ppnl(i,i,l),newline=.TRUE.,error=error) 01574 istr=LEN_TRIM(line_att)+1 01575 WRITE(line_att(istr:),'(E24.16)')hprj_ppnl(i,i,l) 01576 END IF 01577 END IF 01578 DO j=i+1,nprj_ppnl 01579 IF (read_from_input) THEN 01580 READ(line_att,*)hprj_ppnl(i,j,l) 01581 CALL remove_word(line_att) 01582 ELSE 01583 CALL parser_get_object(parser,hprj_ppnl(i,j,l),error=error) 01584 istr=LEN_TRIM(line_att)+1 01585 WRITE(line_att(istr:),'(E24.16)')hprj_ppnl(i,j,l) 01586 END IF 01587 END DO 01588 END DO 01589 IF (.NOT.read_from_input) THEN 01590 irep = irep + 1 01591 CALL section_vals_val_set(potential_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,& 01592 c_val=TRIM(line_att), error=error) 01593 line_att ="" 01594 ELSE 01595 IF (LEN_TRIM(line_att) /= 0) THEN 01596 CALL stop_program(routineN,moduleN,__LINE__,& 01597 "Error reading the Potential from input file!!") 01598 END IF 01599 END IF 01600 IF (nprj_ppnl > 1) THEN 01601 CALL symmetrize_matrix(hprj_ppnl(:,:,l),"upper_to_lower") 01602 END IF 01603 lprj_ppnl_max = MAX(lprj_ppnl_max,l + 2*(nprj_ppnl - 1)) 01604 END DO 01605 01606 potential%nppnl = nppnl 01607 CALL init_orbital_pointers(lprj_ppnl_max) 01608 01609 potential%lprj_ppnl_max = lprj_ppnl_max 01610 potential%nprj_ppnl_max = nprj_ppnl_max 01611 CALL reallocate(potential%hprj_ppnl,1,nprj_ppnl_max,& 01612 1,nprj_ppnl_max,& 01613 0,lppnl) 01614 potential%hprj_ppnl(:,:,:) = hprj_ppnl(:,:,:) 01615 01616 CALL reallocate(potential%cprj,1,ncoset(lprj_ppnl_max),1,nppnl) 01617 CALL reallocate(potential%cprj_ppnl,1,nprj_ppnl_max,0,lppnl) 01618 CALL reallocate(potential%vprj_ppnl,1,nppnl,1,nppnl) 01619 01620 DEALLOCATE (hprj_ppnl,STAT=stat) 01621 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01622 END IF 01623 EXIT search_loop 01624 END IF 01625 ELSE 01626 ! Stop program, if the end of file is reached 01627 CALL stop_program(routineN,moduleN,__LINE__,& 01628 "The requested atomic potential <"//& 01629 TRIM(potential_name)//& 01630 "> was not found in the potential file <"//& 01631 TRIM(potential_file_name)//">",para_env) 01632 END IF 01633 END DO search_loop 01634 01635 IF (.NOT.read_from_input) THEN 01636 ! Dump the potential info the in potential section 01637 IF (match) THEN 01638 irep = irep + 1 01639 WRITE(line_att,'(A)')" # Potential name: "//apname2(:strlen2)//" for symbol: "//symbol2(:strlen1) 01640 CALL section_vals_val_set(potential_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,& 01641 c_val=TRIM(line_att), error=error) 01642 irep = irep + 1 01643 WRITE(line_att,'(A)')" # Potential read from the potential filename: "//TRIM(potential_file_name) 01644 CALL section_vals_val_set(potential_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,& 01645 c_val=TRIM(line_att), error=error) 01646 END IF 01647 CALL parser_release(parser,error=error) 01648 END IF 01649 01650 IF (ASSOCIATED(tmp_vals)) DEALLOCATE(tmp_vals) 01651 01652 END SUBROUTINE read_gth_potential 01653 01654 ! ***************************************************************************** 01655 SUBROUTINE set_default_all_potential(potential,z,zeff_correction,error) 01656 01657 TYPE(all_potential_type), POINTER :: potential 01658 INTEGER, INTENT(IN) :: z 01659 REAL(KIND=dp) :: zeff_correction 01660 TYPE(cp_error_type), INTENT(inout) :: error 01661 01662 CHARACTER(len=*), PARAMETER :: routineN = 'set_default_all_potential', 01663 routineP = moduleN//':'//routineN 01664 01665 CHARACTER(LEN=default_string_length) :: name 01666 INTEGER, DIMENSION(:), POINTER :: elec_conf 01667 REAL(KIND=dp) :: alpha, alpha_core_charge, 01668 ccore_charge, 01669 core_charge_radius, r, zeff 01670 01671 ALLOCATE ( elec_conf(0:3) ) 01672 elec_conf(0:3) = ptable(z)%e_conv(0:3) 01673 zeff = REAL(SUM ( elec_conf ),dp)+zeff_correction 01674 name = ptable(z)%name 01675 01676 r = ptable(z)%covalent_radius * 0.5_dp 01677 r = MAX ( r, 0.2_dp ) 01678 r = MIN ( r, 1.0_dp ) 01679 alpha = 1.0_dp/(2.0_dp*r**2) 01680 01681 core_charge_radius = r 01682 alpha_core_charge = alpha 01683 ccore_charge = zeff*SQRT((alpha/pi)**3) 01684 01685 CALL set_all_potential(potential,& 01686 name=name,& 01687 alpha_core_charge=alpha_core_charge,& 01688 ccore_charge=ccore_charge,& 01689 core_charge_radius=core_charge_radius,& 01690 z=z,& 01691 zeff=zeff,& 01692 zeff_correction=zeff_correction,& 01693 elec_conf=elec_conf) 01694 01695 DEALLOCATE ( elec_conf ) 01696 01697 END SUBROUTINE set_default_all_potential 01698 01699 ! ***************************************************************************** 01705 SUBROUTINE set_all_potential(potential,name,alpha_core_charge,& 01706 ccore_charge,core_charge_radius,z,zeff,& 01707 zeff_correction,elec_conf) 01708 01709 TYPE(all_potential_type), POINTER :: potential 01710 CHARACTER(LEN=default_string_length), 01711 INTENT(IN), OPTIONAL :: name 01712 REAL(KIND=dp), INTENT(IN), OPTIONAL :: alpha_core_charge, 01713 ccore_charge, 01714 core_charge_radius 01715 INTEGER, INTENT(IN), OPTIONAL :: z 01716 REAL(KIND=dp), INTENT(IN), OPTIONAL :: zeff, zeff_correction 01717 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: elec_conf 01718 01719 CHARACTER(len=*), PARAMETER :: routineN = 'set_all_potential', 01720 routineP = moduleN//':'//routineN 01721 01722 IF (ASSOCIATED(potential)) THEN 01723 01724 IF (PRESENT(name)) potential%name = name 01725 IF (PRESENT(alpha_core_charge))& 01726 potential%alpha_core_charge = alpha_core_charge 01727 IF (PRESENT(ccore_charge)) potential%ccore_charge = ccore_charge 01728 IF (PRESENT(core_charge_radius))& 01729 potential%core_charge_radius = core_charge_radius 01730 IF (PRESENT(z)) potential%z = z 01731 IF (PRESENT(zeff)) potential%zeff = zeff 01732 IF (PRESENT(zeff_correction)) potential%zeff_correction = zeff_correction 01733 IF (PRESENT(elec_conf)) THEN 01734 IF ( .NOT. ASSOCIATED( potential%elec_conf ) ) THEN 01735 CALL reallocate(potential%elec_conf,0,SIZE(elec_conf)-1) 01736 END IF 01737 potential%elec_conf(:) = elec_conf(:) 01738 END IF 01739 01740 ELSE 01741 01742 CALL stop_program(routineN,moduleN,__LINE__,& 01743 "The pointer potential is not associated") 01744 01745 END IF 01746 01747 END SUBROUTINE set_all_potential 01748 01749 ! ***************************************************************************** 01755 SUBROUTINE set_fist_potential(potential, apol, cpol, qeff, mm_radius, & 01756 qmmm_corr_radius, qmmm_radius) 01757 TYPE(fist_potential_type), POINTER :: potential 01758 REAL(KIND=dp), INTENT(IN), OPTIONAL :: apol, cpol, qeff, mm_radius, 01759 qmmm_corr_radius, qmmm_radius 01760 01761 CHARACTER(len=*), PARAMETER :: routineN = 'set_fist_potential', 01762 routineP = moduleN//':'//routineN 01763 01764 IF (ASSOCIATED(potential)) THEN 01765 01766 IF (PRESENT(apol)) potential%apol = apol 01767 IF (PRESENT(cpol)) potential%cpol = cpol 01768 IF (PRESENT(mm_radius)) potential%mm_radius = mm_radius 01769 IF (PRESENT(qeff)) potential%qeff = qeff 01770 IF (PRESENT(qmmm_corr_radius)) potential%qmmm_corr_radius = qmmm_corr_radius 01771 IF (PRESENT(qmmm_radius)) potential%qmmm_radius = qmmm_radius 01772 01773 ELSE 01774 01775 CALL stop_program(routineN,moduleN,__LINE__,& 01776 "The pointer potential is not associated") 01777 01778 END IF 01779 01780 END SUBROUTINE set_fist_potential 01781 01782 ! ***************************************************************************** 01788 SUBROUTINE set_gth_potential(potential,name,alpha_core_charge,alpha_ppl,& 01789 ccore_charge,cerf_ppl,core_charge_radius,& 01790 ppl_radius,ppnl_radius,lppnl,lprj_ppnl_max,& 01791 nexp_ppl,nppnl,nprj_ppnl_max,z,zeff,zeff_correction,& 01792 alpha_ppnl,cexp_ppl,elec_conf,nprj_ppnl,cprj,cprj_ppnl,& 01793 vprj_ppnl,hprj_ppnl) 01794 TYPE(gth_potential_type), POINTER :: potential 01795 CHARACTER(LEN=default_string_length), 01796 INTENT(IN), OPTIONAL :: name 01797 REAL(KIND=dp), INTENT(IN), OPTIONAL :: alpha_core_charge, alpha_ppl, 01798 ccore_charge, cerf_ppl, core_charge_radius, ppl_radius, ppnl_radius 01799 INTEGER, INTENT(IN), OPTIONAL :: lppnl, lprj_ppnl_max, 01800 nexp_ppl, nppnl, 01801 nprj_ppnl_max, z 01802 REAL(KIND=dp), INTENT(IN), OPTIONAL :: zeff, zeff_correction 01803 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 01804 POINTER :: alpha_ppnl, cexp_ppl 01805 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: elec_conf, nprj_ppnl 01806 REAL(KIND=dp), DIMENSION(:, :), 01807 OPTIONAL, POINTER :: cprj, cprj_ppnl, vprj_ppnl 01808 REAL(KIND=dp), DIMENSION(:, :, :), 01809 OPTIONAL, POINTER :: hprj_ppnl 01810 01811 CHARACTER(len=*), PARAMETER :: routineN = 'set_gth_potential', 01812 routineP = moduleN//':'//routineN 01813 01814 IF (ASSOCIATED(potential)) THEN 01815 01816 IF (PRESENT(name)) potential%name = name 01817 IF (PRESENT(alpha_core_charge))& 01818 potential%alpha_core_charge = alpha_core_charge 01819 IF (PRESENT(alpha_ppl)) potential%alpha_ppl = alpha_ppl 01820 IF (PRESENT(ccore_charge)) potential%ccore_charge = ccore_charge 01821 IF (PRESENT(cerf_ppl)) potential%cerf_ppl = cerf_ppl 01822 IF (PRESENT(core_charge_radius))& 01823 potential%core_charge_radius = core_charge_radius 01824 IF (PRESENT(ppl_radius)) potential%ppl_radius = ppl_radius 01825 IF (PRESENT(ppnl_radius)) potential%ppnl_radius = ppnl_radius 01826 IF (PRESENT(lppnl)) potential%lppnl = lppnl 01827 IF (PRESENT(lprj_ppnl_max)) potential%lprj_ppnl_max = lprj_ppnl_max 01828 IF (PRESENT(nexp_ppl)) potential%nexp_ppl = nexp_ppl 01829 IF (PRESENT(nppnl)) potential%nppnl = nppnl 01830 IF (PRESENT(nprj_ppnl_max)) potential%nprj_ppnl_max = nprj_ppnl_max 01831 IF (PRESENT(z)) potential%z = z 01832 IF (PRESENT(zeff)) potential%zeff = zeff 01833 IF (PRESENT(zeff_correction)) potential%zeff_correction = zeff_correction 01834 IF (PRESENT(alpha_ppnl)) potential%alpha_ppnl => alpha_ppnl 01835 IF (PRESENT(cexp_ppl)) potential%cexp_ppl => cexp_ppl 01836 IF (PRESENT(elec_conf)) potential%elec_conf => elec_conf 01837 IF (PRESENT(nprj_ppnl)) potential%nprj_ppnl => nprj_ppnl 01838 IF (PRESENT(cprj)) potential%cprj => cprj 01839 IF (PRESENT(cprj_ppnl)) potential%cprj_ppnl => cprj_ppnl 01840 IF (PRESENT(vprj_ppnl)) potential%vprj_ppnl => vprj_ppnl 01841 IF (PRESENT(hprj_ppnl)) potential%hprj_ppnl => hprj_ppnl 01842 01843 ELSE 01844 01845 CALL stop_program(routineN,moduleN,__LINE__,& 01846 "The pointer potential is not associated") 01847 01848 END IF 01849 01850 END SUBROUTINE set_gth_potential 01851 01852 ! ***************************************************************************** 01853 SUBROUTINE write_all_potential(potential,output_unit,error) 01854 01855 ! Write an atomic all-electron potential data set to the output unit. 01856 01857 ! - Creation (09.02.2002,MK) 01858 01859 TYPE(all_potential_type), POINTER :: potential 01860 INTEGER, INTENT(in) :: output_unit 01861 TYPE(cp_error_type), INTENT(inout) :: error 01862 01863 CHARACTER(LEN=20) :: string 01864 01865 IF (output_unit > 0.AND.ASSOCIATED(potential)) THEN 01866 WRITE (UNIT=output_unit,FMT="(/,T6,A,T41,A40,/)")& 01867 "Potential information for",ADJUSTR(TRIM(potential%name)) 01868 WRITE (UNIT=output_unit,FMT="(T8,A,T41,A40)")& 01869 "Description: ",TRIM(potential%description(1)),& 01870 " ",TRIM(potential%description(2)) 01871 WRITE (UNIT=output_unit,FMT="(/,T8,A,T69,F12.6)")& 01872 "Gaussian exponent of the core charge distribution: ",& 01873 potential%alpha_core_charge 01874 WRITE (UNIT=string,FMT="(5I4)") potential%elec_conf 01875 WRITE (UNIT=output_unit,FMT="(T8,A,T61,A20)")& 01876 "Electronic configuration (s p d ...):",& 01877 ADJUSTR(TRIM(string)) 01878 END IF 01879 01880 END SUBROUTINE write_all_potential 01881 01882 ! ***************************************************************************** 01883 SUBROUTINE write_gth_potential(potential,output_unit,error) 01884 01885 ! Write an atomic GTH potential data set to the output unit. 01886 ! - Creation (09.02.2002,MK) 01887 01888 TYPE(gth_potential_type), POINTER :: potential 01889 INTEGER, INTENT(in) :: output_unit 01890 TYPE(cp_error_type), INTENT(inout) :: error 01891 01892 CHARACTER(LEN=20) :: string 01893 INTEGER :: i, j, l 01894 REAL(KIND=dp) :: r 01895 01896 IF (output_unit>0.AND.ASSOCIATED(potential)) THEN 01897 WRITE (UNIT=output_unit,FMT="(/,T6,A,T41,A40,/)")& 01898 "Potential information for",ADJUSTR(TRIM(potential%name)) 01899 WRITE (UNIT=output_unit,FMT="(T8,A,T41,A40)")& 01900 "Description: ",ADJUSTR(TRIM(potential%description(1))),& 01901 " ",ADJUSTR(TRIM(potential%description(2))),& 01902 " ",ADJUSTR(TRIM(potential%description(3))),& 01903 " ",ADJUSTR(TRIM(potential%description(4))) 01904 WRITE (UNIT=output_unit,FMT="(/,T8,A,T69,F12.6)")& 01905 "Gaussian exponent of the core charge distribution: ",& 01906 potential%alpha_core_charge 01907 WRITE (UNIT=string,FMT="(5I4)") potential%elec_conf 01908 WRITE (UNIT=output_unit,FMT="(T8,A,T61,A20)")& 01909 "Electronic configuration (s p d ...):",& 01910 ADJUSTR(TRIM(string)) 01911 01912 r = 1.0_dp/SQRT(2.0_dp*potential%alpha_ppl) 01913 01914 WRITE (UNIT=output_unit,FMT="(/,T8,A,/,/,T27,A,/,T21,5F12.6)")& 01915 "Parameters of the local part of the GTH pseudopotential:",& 01916 "rloc C1 C2 C3 C4",& 01917 r,(potential%cexp_ppl(i)*r**(2*(i-1)),i=1,potential%nexp_ppl) 01918 01919 IF (potential%lppnl > -1) THEN 01920 WRITE (UNIT=output_unit,FMT="(/,T8,A,/,/,T20,A,/)")& 01921 "Parameters of the non-local part of the GTH pseudopotential:",& 01922 "l r(l) h(i,j,l)" 01923 DO l=0,potential%lppnl 01924 r = SQRT(0.5_dp/potential%alpha_ppnl(l)) 01925 WRITE (UNIT=output_unit,FMT="(T19,I2,5F12.6)")& 01926 l,r,(potential%hprj_ppnl(1,j,l),j=1,potential%nprj_ppnl(l)) 01927 DO i=2,potential%nprj_ppnl(l) 01928 WRITE (UNIT=output_unit,FMT="(T33,4F12.6)")& 01929 (potential%hprj_ppnl(i,j,l),j=1,potential%nprj_ppnl(l)) 01930 END DO 01931 END DO 01932 END IF 01933 END IF 01934 01935 END SUBROUTINE write_gth_potential 01936 01937 ! ***************************************************************************** 01938 01939 END MODULE external_potential_types
1.7.3