|
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 ! ***************************************************************************** 00009 MODULE atom_kind_orbitals 00010 USE atom_electronic_structure, ONLY: calculate_atom 00011 USE atom_fit, ONLY: atom_fit_density 00012 USE atom_operators, ONLY: atom_int_release,& 00013 atom_int_setup,& 00014 atom_ppint_release,& 00015 atom_ppint_setup,& 00016 atom_relint_release,& 00017 atom_relint_setup 00018 USE atom_types, ONLY: & 00019 CGTO_BASIS, Clementi_geobas, GTO_BASIS, atom_basis_type, & 00020 atom_integrals, atom_orbitals, atom_potential_type, atom_type, & 00021 create_atom_orbs, create_atom_type, release_atom_basis, & 00022 release_atom_potential, release_atom_type, set_atom 00023 USE atom_utils, ONLY: get_maxl_occ,& 00024 get_maxn_occ 00025 USE atomic_kind_types, ONLY: atomic_kind_type,& 00026 get_atomic_kind 00027 USE basis_set_types, ONLY: get_gto_basis_set,& 00028 gto_basis_set_type 00029 USE external_potential_types, ONLY: all_potential_type,& 00030 get_potential,& 00031 gth_potential_type 00032 USE input_constants, ONLY: do_analytic,& 00033 do_dkh2_atom,& 00034 do_gapw_log,& 00035 do_nonrel_atom,& 00036 do_numeric,& 00037 do_rks_atom,& 00038 gth_pseudo,& 00039 no_pseudo 00040 USE input_section_types, ONLY: section_vals_type 00041 USE kinds, ONLY: dp 00042 USE periodic_table, ONLY: ptable 00043 USE physcon, ONLY: bohr 00044 USE qs_grid_atom, ONLY: allocate_grid_atom,& 00045 create_grid_atom,& 00046 grid_atom_type 00047 #include "cp_common_uses.h" 00048 00049 IMPLICIT NONE 00050 00051 PRIVATE 00052 00053 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_kind_orbitals' 00054 00055 PUBLIC :: calculate_atomic_orbitals, calculate_atomic_density 00056 00057 ! ***************************************************************************** 00058 00059 CONTAINS 00060 00061 ! ***************************************************************************** 00062 SUBROUTINE calculate_atomic_orbitals (atomic_kind,iunit,pmat,ispin,confine,& 00063 xc_section,error) 00064 TYPE(atomic_kind_type), POINTER :: atomic_kind 00065 INTEGER, INTENT(IN), OPTIONAL :: iunit 00066 REAL(KIND=dp), DIMENSION(:, :), 00067 OPTIONAL, POINTER :: pmat 00068 INTEGER, INTENT(IN) :: ispin 00069 LOGICAL, INTENT(IN), OPTIONAL :: confine 00070 TYPE(section_vals_type), OPTIONAL, 00071 POINTER :: xc_section 00072 TYPE(cp_error_type), INTENT(INOUT) :: error 00073 00074 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_atomic_orbitals', 00075 routineP = moduleN//':'//routineN 00076 00077 INTEGER :: i, ierr, ii, ipgf, j, k, k1, k2, l, ll, lm, m, mb, mo, n, ne, 00078 nexp_lpot, nexp_lsd, nexp_nlcc, ngp, nj, nn, nr, ns, nset, nsgf, 00079 quadtype, z 00080 INTEGER, DIMENSION(0:3) :: econfx 00081 INTEGER, DIMENSION(0:3, 10) :: ncalc, ncore, nelem 00082 INTEGER, DIMENSION(0:3, 100) :: set_index, shell_index 00083 INTEGER, DIMENSION(:), POINTER :: econf, lmax, lmin, nct_lpot, 00084 nct_lsd, nct_nlcc, npgf, 00085 nppnl, nshell, ppeconf 00086 INTEGER, DIMENSION(:, :), POINTER :: addel, first_sgf, laddel, 00087 last_sgf, ls, naddel 00088 LOGICAL :: bs_occupation, 00089 failure = .FALSE., 00090 lpot_present, lsd_present, 00091 nlcc_present 00092 REAL(KIND=dp) :: ac, al, ear, rk, scal, zeff 00093 REAL(KIND=dp), DIMENSION(:), POINTER :: alpha_lpot, alpha_lsd, 00094 alpha_nlcc, ap, ce 00095 REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_lpot, cval_lsd, 00096 cval_nlcc, zet 00097 REAL(KIND=dp), DIMENSION(:, :, :), 00098 POINTER :: gcc, hp 00099 TYPE(all_potential_type), POINTER :: all_potential 00100 TYPE(atom_basis_type), POINTER :: basis 00101 TYPE(atom_integrals), POINTER :: integrals 00102 TYPE(atom_orbitals), POINTER :: orbitals 00103 TYPE(atom_potential_type), POINTER :: potential 00104 TYPE(atom_type), POINTER :: atom 00105 TYPE(grid_atom_type), POINTER :: grid 00106 TYPE(gth_potential_type), POINTER :: gth_potential 00107 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 00108 00109 NULLIFY(atom) 00110 CALL create_atom_type(atom,error) 00111 00112 IF(PRESENT(xc_section)) THEN 00113 atom%xc_section => xc_section 00114 ELSE 00115 NULLIFY(atom%xc_section) 00116 END IF 00117 00118 CALL get_atomic_kind(atomic_kind,z=z,zeff=zeff) 00119 NULLIFY(all_potential,gth_potential,orb_basis_set) 00120 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00121 all_potential=all_potential,& 00122 gth_potential=gth_potential,& 00123 orb_basis_set=orb_basis_set,& 00124 bs_occupation=bs_occupation,& 00125 addel=addel,laddel=laddel,naddel=naddel) 00126 CPPostcondition(ASSOCIATED(orb_basis_set), cp_failure_level, routineP, error, failure) 00127 00128 atom%z = z 00129 CALL set_atom(atom,& 00130 pp_calc=ASSOCIATED(gth_potential),& 00131 method_type=do_rks_atom,& 00132 relativistic=do_nonrel_atom,& 00133 coulomb_integral_type=do_numeric,& 00134 exchange_integral_type=do_numeric,& 00135 error=error) 00136 00137 ALLOCATE (potential,basis,integrals,STAT=ierr) 00138 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00139 00140 IF ( ASSOCIATED(gth_potential) ) THEN 00141 potential%ppot_type=gth_pseudo 00142 IF ( PRESENT(confine) ) THEN 00143 potential%confinement=confine 00144 ELSE 00145 potential%confinement=.TRUE. 00146 END IF 00147 potential%acon=0.1_dp 00148 potential%rcon=2._dp*ptable(z)%vdw_radius*bohr 00149 potential%ncon=2 00150 00151 CALL get_potential(gth_potential,& 00152 zeff=zeff,& 00153 elec_conf=ppeconf,& 00154 alpha_core_charge=ac,& 00155 nexp_ppl=ne,& 00156 cexp_ppl=ce,& 00157 lppnl=lm,& 00158 nprj_ppnl=nppnl,& 00159 alpha_ppnl=ap,& 00160 hprj_ppnl=hp) 00161 00162 potential%gth_pot%zion = zeff 00163 potential%gth_pot%rc = SQRT(0.5_dp/ac) 00164 potential%gth_pot%ncl = ne 00165 potential%gth_pot%cl(:) = 0._dp 00166 IF (ac > 0._dp) THEN 00167 DO i=1,ne 00168 potential%gth_pot%cl(i) = ce(i)/(2._dp*ac)**(i-1) 00169 END DO 00170 END IF 00171 !extended type 00172 potential%gth_pot%lpotextended = .FALSE. 00173 potential%gth_pot%lsdpot = .FALSE. 00174 potential%gth_pot%nlcc = .FALSE. 00175 potential%gth_pot%nexp_lpot = 0 00176 potential%gth_pot%nexp_lsd = 0 00177 potential%gth_pot%nexp_nlcc = 0 00178 CALL get_potential(gth_potential,& 00179 lpot_present=lpot_present,& 00180 lsd_present=lsd_present,& 00181 nlcc_present=nlcc_present) 00182 IF (lpot_present) THEN 00183 CALL get_potential(gth_potential,& 00184 nexp_lpot=nexp_lpot,& 00185 alpha_lpot=alpha_lpot,& 00186 nct_lpot=nct_lpot,& 00187 cval_lpot=cval_lpot) 00188 potential%gth_pot%lpotextended = .TRUE. 00189 potential%gth_pot%nexp_lpot = nexp_lpot 00190 potential%gth_pot%alpha_lpot(1:nexp_lpot) = SQRT(0.5_dp/alpha_lpot(1:nexp_lpot)) 00191 potential%gth_pot%nct_lpot(1:nexp_lpot) = nct_lpot(1:nexp_lpot) 00192 DO j=1,nexp_lpot 00193 ac = alpha_lpot(j) 00194 DO i=1,4 00195 potential%gth_pot%cval_lpot(i,j) = cval_lpot(i,j)/(2._dp*ac)**(i-1) 00196 END DO 00197 END DO 00198 END IF 00199 IF (lsd_present) THEN 00200 CALL get_potential(gth_potential,& 00201 nexp_lsd=nexp_lsd,& 00202 alpha_lsd=alpha_lsd,& 00203 nct_lsd=nct_lsd,& 00204 cval_lsd=cval_lsd) 00205 potential%gth_pot%lsdpot = .TRUE. 00206 potential%gth_pot%nexp_lsd = nexp_lsd 00207 potential%gth_pot%alpha_lsd(1:nexp_lsd) = SQRT(0.5_dp/alpha_lsd(1:nexp_lsd)) 00208 potential%gth_pot%nct_lsd(1:nexp_lsd) = nct_lsd(1:nexp_lsd) 00209 DO j=1,nexp_lpot 00210 ac = alpha_lsd(j) 00211 DO i=1,4 00212 potential%gth_pot%cval_lsd(i,j) = cval_lsd(i,j)/(2._dp*ac)**(i-1) 00213 END DO 00214 END DO 00215 END IF 00216 00217 ! nonlocal part 00218 potential%gth_pot%nl(:) = 0 00219 potential%gth_pot%rcnl(:) = 0._dp 00220 potential%gth_pot%hnl(:,:,:)= 0._dp 00221 DO l=0,lm 00222 n = nppnl(l) 00223 potential%gth_pot%nl(l) = n 00224 potential%gth_pot%rcnl(l) = SQRT(0.5_dp/ap(l)) 00225 potential%gth_pot%hnl(1:n,1:n,l)= hp(1:n,1:n,l) 00226 END DO 00227 00228 IF (nlcc_present) THEN 00229 CALL get_potential(gth_potential,& 00230 nexp_nlcc=nexp_nlcc,& 00231 alpha_nlcc=alpha_nlcc,& 00232 nct_nlcc=nct_nlcc,& 00233 cval_nlcc=cval_nlcc) 00234 potential%gth_pot%nlcc = .TRUE. 00235 potential%gth_pot%nexp_nlcc = nexp_nlcc 00236 potential%gth_pot%alpha_nlcc(1:nexp_nlcc) = alpha_nlcc(1:nexp_nlcc) 00237 potential%gth_pot%nct_nlcc(1:nexp_nlcc) = nct_nlcc(1:nexp_nlcc) 00238 potential%gth_pot%cval_nlcc(1:4,1:nexp_nlcc) = cval_nlcc(1:4,1:nexp_nlcc) 00239 END IF 00240 00241 CALL set_atom(atom,zcore=NINT(zeff),potential=potential,error=error) 00242 ELSE 00243 potential%ppot_type=no_pseudo 00244 IF ( PRESENT(confine) ) THEN 00245 potential%confinement=confine 00246 ELSE 00247 potential%confinement=.FALSE. 00248 END IF 00249 potential%acon=0.1_dp 00250 potential%rcon=2._dp*ptable(z)%vdw_radius*bohr 00251 potential%ncon=2 00252 CALL set_atom(atom,zcore=z,potential=potential,error=error) 00253 END IF 00254 00255 CALL get_gto_basis_set(orb_basis_set,& 00256 nset=nset,nshell=nshell,npgf=npgf,lmin=lmin,lmax=lmax,l=ls,nsgf=nsgf,zet=zet,gcc=gcc,& 00257 first_sgf=first_sgf,last_sgf=last_sgf) 00258 00259 NULLIFY(grid) 00260 ngp = 400 00261 quadtype = do_gapw_log 00262 CALL allocate_grid_atom(grid,error) 00263 CALL create_grid_atom(grid,ngp,1,1,0,quadtype) 00264 grid%nr = ngp 00265 basis%grid => grid 00266 00267 NULLIFY(basis%am,basis%cm,basis%as,basis%ns,basis%bf,basis%dbf) 00268 basis%basis_type = CGTO_BASIS 00269 basis%eps_eig = 1.e-12_dp 00270 00271 ! fill in the basis data structures 00272 set_index = 0 00273 shell_index = 0 00274 basis%nprim = 0 00275 basis%nbas = 0 00276 DO i=1,nset 00277 DO j=lmin(i),MIN(lmax(i),3) 00278 basis%nprim(j) = basis%nprim(j) + npgf(i) 00279 END DO 00280 DO j=1,nshell(i) 00281 l = ls(j,i) 00282 IF ( l <= 3 ) THEN 00283 basis%nbas(l) = basis%nbas(l)+1 00284 k = basis%nbas(l) 00285 CPPostcondition(k<=100, cp_failure_level, routineP, error, failure) 00286 set_index(l,k) = i 00287 shell_index(l,k) = j 00288 END IF 00289 END DO 00290 END DO 00291 00292 nj = MAXVAL(basis%nprim) 00293 ns = MAXVAL(basis%nbas) 00294 ALLOCATE (basis%am(nj,0:3),STAT=ierr) 00295 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00296 basis%am = 0._dp 00297 ALLOCATE (basis%cm(nj,ns,0:3),STAT=ierr) 00298 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00299 basis%cm = 0._dp 00300 DO j=0,3 00301 nj = 0 00302 ns = 0 00303 DO i=1,nset 00304 IF (j >= lmin(i) .AND. j <= lmax(i)) THEN 00305 DO ipgf = 1,npgf(i) 00306 basis%am(nj+ipgf,j) = zet(ipgf,i) 00307 END DO 00308 DO ii=1,nshell(i) 00309 IF ( ls(ii,i) == j ) THEN 00310 ns = ns + 1 00311 DO ipgf=1,npgf(i) 00312 basis%cm(nj+ipgf,ns,j) = gcc(ipgf,ii,i) 00313 END DO 00314 END IF 00315 END DO 00316 nj = nj + npgf(i) 00317 END IF 00318 END DO 00319 END DO 00320 00321 ! initialize basis function on a radial grid 00322 nr = basis%grid%nr 00323 m = MAXVAL(basis%nbas) 00324 ALLOCATE (basis%bf(nr,m,0:3),STAT=ierr) 00325 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00326 ALLOCATE (basis%dbf(nr,m,0:3),STAT=ierr) 00327 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00328 00329 basis%bf = 0._dp 00330 basis%dbf = 0._dp 00331 DO l=0,3 00332 DO i=1,basis%nprim(l) 00333 al = basis%am(i,l) 00334 DO k=1,nr 00335 rk = basis%grid%rad(k) 00336 ear = EXP(-al*basis%grid%rad(k)**2) 00337 DO j=1,basis%nbas(l) 00338 basis%bf(k,j,l) = basis%bf(k,j,l) + rk**l * ear*basis%cm(i,j,l) 00339 basis%dbf(k,j,l) = basis%dbf(k,j,l) & 00340 + ( REAL(l,dp)*rk**(l-1) - 2._dp*al*rk**(l+1) ) * ear*basis%cm(i,j,l) 00341 END DO 00342 END DO 00343 END DO 00344 END DO 00345 00346 CALL set_atom(atom,basis=basis,error=error) 00347 00348 ! optimization defaults 00349 atom%optimization%damping = 0.2_dp 00350 atom%optimization%eps_scf = 1.e-6_dp 00351 atom%optimization%eps_diis = 100._dp 00352 atom%optimization%max_iter = 50 00353 atom%optimization%n_diis = 5 00354 00355 ! electronic state 00356 nelem = 0 00357 ncore = 0 00358 ncalc = 0 00359 NULLIFY(econf) 00360 IF ( ASSOCIATED(gth_potential) ) THEN 00361 CALL get_potential(gth_potential,elec_conf=econf) 00362 econfx = 0 00363 econfx(0:SIZE(econf)-1) = econf 00364 IF ( SUM(econf) >= 0 ) THEN 00365 DO l=0,3 00366 ll = 2*(2*l+1) 00367 nn = ptable(z)%e_conv(l)-econfx(l) 00368 ii = 0 00369 DO 00370 ii = ii + 1 00371 IF(nn <= ll) THEN 00372 ncore(l,ii) = nn 00373 EXIT 00374 ELSE 00375 ncore(l,ii) = ll 00376 nn = nn - ll 00377 END IF 00378 END DO 00379 END DO 00380 DO l=0,3 00381 ll = 2*(2*l+1) 00382 nn = ptable(z)%e_conv(l) 00383 ii = 0 00384 DO 00385 ii = ii + 1 00386 IF(nn <= ll) THEN 00387 nelem(l,ii) = nn 00388 EXIT 00389 ELSE 00390 nelem(l,ii) = ll 00391 nn = nn - ll 00392 END IF 00393 END DO 00394 END DO 00395 ncalc = nelem - ncore 00396 ELSE 00397 ncore = 0 00398 ncalc = 0 00399 DO l=0,3 00400 ll = 2*(2*l+1) 00401 nn = ABS(econfx(l)) 00402 ii = 0 00403 DO 00404 ii = ii + 1 00405 IF(nn <= ll) THEN 00406 ncalc(l,ii) = -nn 00407 EXIT 00408 ELSE 00409 ncalc(l,ii) = -ll 00410 nn = nn - ll 00411 END IF 00412 END DO 00413 END DO 00414 END IF 00415 ELSE 00416 DO l=0,3 00417 ll = 2*(2*l+1) 00418 nn = ptable(z)%e_conv(l) 00419 ii = 0 00420 DO 00421 ii = ii + 1 00422 IF(nn <= ll) THEN 00423 nelem(l,ii) = nn 00424 EXIT 00425 ELSE 00426 nelem(l,ii) = ll 00427 nn = nn - ll 00428 END IF 00429 END DO 00430 END DO 00431 ncalc = nelem - ncore 00432 END IF 00433 00434 IF(bs_occupation) THEN 00435 ! readjust the occupation number of the atomic orbitals 00436 ! according to the changes required from input in order to bias the initial guess 00437 DO i=1,SIZE(addel,1) 00438 ne=addel(i,ispin) 00439 l=laddel(i,ispin) 00440 nn=naddel(i,ispin)-l 00441 IF(ne/=0) THEN 00442 IF(nn==0) THEN 00443 DO ii=SIZE(nelem,2),1,-1 00444 IF(ncalc(l,ii)>0) THEN 00445 IF((ncalc(l,ii)+ne) < 2*(2*l+1)+1) THEN 00446 ncalc(l,ii) = ncalc(l,ii)+ne 00447 nn = ii 00448 ELSE 00449 ncalc(l,ii+1) = ncalc(l,ii+1)+ne 00450 nn = ii + 1 00451 END IF 00452 EXIT 00453 ELSE IF (ii==1) THEN 00454 ncalc(l,ii)=ncalc(l,ii)+ne 00455 nn = ii 00456 END IF 00457 END DO 00458 ELSE 00459 ncalc(l,nn) = ncalc(l,nn) + ne 00460 END IF 00461 IF( ncalc(l,nn)<0 ) THEN 00462 ncalc(l,nn) = 0 00463 END IF 00464 END IF 00465 END DO 00466 END IF 00467 00468 IF ( atomic_kind%ghost ) THEN 00469 nelem = 0 00470 ncore = 0 00471 ncalc = 0 00472 END IF 00473 00474 ALLOCATE (atom%state,STAT=ierr) 00475 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00476 00477 atom%state%core = 0._dp 00478 atom%state%core(0:3,1:7) = REAL(ncore(0:3,1:7),dp) 00479 atom%state%occ = 0._dp 00480 atom%state%occ(0:3,1:7) = REAL(ncalc(0:3,1:7),dp) 00481 atom%state%occupation = 0._dp 00482 atom%state%multiplicity = -1 00483 DO l=0,3 00484 k = 0 00485 DO i=1,7 00486 IF ( ncalc(l,i) > 0 ) THEN 00487 k = k + 1 00488 atom%state%occupation(l,k) = REAL(ncalc(l,i),dp) 00489 END IF 00490 END DO 00491 END DO 00492 00493 00494 atom%state%maxl_occ = get_maxl_occ(atom%state%occupation) 00495 atom%state%maxn_occ = get_maxn_occ(atom%state%occupation) 00496 atom%state%maxl_calc = atom%state%maxl_occ 00497 atom%state%maxn_calc = atom%state%maxn_occ 00498 00499 ! calculate integrals 00500 ! general integrals 00501 CALL atom_int_setup(integrals,basis,potential=atom%potential,& 00502 eri_coulomb=(atom%coulomb_integral_type==do_analytic),& 00503 eri_exchange=(atom%exchange_integral_type==do_analytic),error=error) 00504 ! potential 00505 CALL atom_ppint_setup(integrals,basis,potential=atom%potential,error=error) 00506 ! relativistic correction terms 00507 NULLIFY(integrals%tzora,integrals%hdkh) 00508 CALL atom_relint_setup(integrals,basis,atom%relativistic,zcore=REAL(atom%zcore,dp),error=error) 00509 CALL set_atom(atom,integrals=integrals,error=error) 00510 00511 NULLIFY(orbitals) 00512 mo = MAXVAL(atom%state%maxn_calc) 00513 mb = MAXVAL(atom%basis%nbas) 00514 CALL create_atom_orbs(orbitals,mb,mo,error) 00515 CALL set_atom(atom,orbitals=orbitals,error=error) 00516 00517 IF(PRESENT(iunit)) THEN 00518 CALL calculate_atom(atom,iunit,error=error) 00519 ELSE 00520 CALL calculate_atom(atom,-1,error=error) 00521 END IF 00522 00523 IF (PRESENT(pmat)) THEN 00524 ! recover density matrix in CP2K/GPW order and normalization 00525 IF(ASSOCIATED(pmat)) THEN 00526 DEALLOCATE (pmat,STAT=ierr) 00527 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00528 END IF 00529 ALLOCATE (pmat(nsgf,nsgf),STAT=ierr) 00530 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00531 pmat = 0._dp 00532 DO l=0,3 00533 ll = 2*l 00534 DO k1=1,atom%basis%nbas(l) 00535 DO k2=1,atom%basis%nbas(l) 00536 scal=SQRT(atom%integrals%ovlp(k1,k1,l)*atom%integrals%ovlp(k2,k2,l))/REAL(2*l+1,KIND=dp) 00537 i=first_sgf(shell_index(l,k1),set_index(l,k1)) 00538 j=first_sgf(shell_index(l,k2),set_index(l,k2)) 00539 DO m=0,ll 00540 pmat(i+m,j+m) = atom%orbitals%pmat(k1,k2,l)*scal 00541 END DO 00542 END DO 00543 END DO 00544 ENDDO 00545 END IF 00546 00547 ! clean up 00548 CALL atom_int_release(integrals,error) 00549 CALL atom_ppint_release(integrals,error) 00550 CALL atom_relint_release(integrals,error) 00551 CALL release_atom_basis(basis,error) 00552 CALL release_atom_potential(potential,error) 00553 CALL release_atom_type(atom,error) 00554 00555 DEALLOCATE (potential,basis,integrals,STAT=ierr) 00556 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00557 00558 END SUBROUTINE calculate_atomic_orbitals 00559 00560 ! ***************************************************************************** 00561 00562 SUBROUTINE calculate_atomic_density (density,atomic_kind,ngto,iunit,allelectron,confine,error) 00563 REAL(KIND=dp), DIMENSION(:, :), 00564 INTENT(INOUT) :: density 00565 TYPE(atomic_kind_type), POINTER :: atomic_kind 00566 INTEGER, INTENT(IN) :: ngto 00567 INTEGER, INTENT(IN), OPTIONAL :: iunit 00568 LOGICAL, INTENT(IN), OPTIONAL :: allelectron, confine 00569 TYPE(cp_error_type), INTENT(INOUT) :: error 00570 00571 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_atomic_density', 00572 routineP = moduleN//':'//routineN 00573 INTEGER, PARAMETER :: num_gto = 40 00574 00575 INTEGER :: i, ierr, ii, j, k, l, ll, lm, m, mb, mo, n, ne, nexp_lpot, 00576 nexp_lsd, nexp_nlcc, ngp, nn, nr, quadtype, relativistic, z 00577 INTEGER, DIMENSION(0:3) :: econfx, starti 00578 INTEGER, DIMENSION(0:3, 10) :: ncalc, ncore, nelem 00579 INTEGER, DIMENSION(:), POINTER :: econf, nct_lpot, nct_lsd, 00580 nct_nlcc, nppnl, ppeconf 00581 LOGICAL :: failure = .FALSE., 00582 lpot_present, lsd_present, 00583 nlcc_present 00584 REAL(KIND=dp) :: ac, al, aval, cc, cval, ear, 00585 rk, xx, zeff 00586 REAL(KIND=dp), DIMENSION(:), POINTER :: alpha_lpot, alpha_lsd, 00587 alpha_nlcc, ap, ce 00588 REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_lpot, cval_lsd, cval_nlcc 00589 REAL(KIND=dp), DIMENSION(:, :, :), 00590 POINTER :: hp 00591 REAL(KIND=dp), DIMENSION(num_gto+2) :: results 00592 TYPE(all_potential_type), POINTER :: all_potential 00593 TYPE(atom_basis_type), POINTER :: basis 00594 TYPE(atom_integrals), POINTER :: integrals 00595 TYPE(atom_orbitals), POINTER :: orbitals 00596 TYPE(atom_potential_type), POINTER :: potential 00597 TYPE(atom_type), POINTER :: atom 00598 TYPE(grid_atom_type), POINTER :: grid 00599 TYPE(gth_potential_type), POINTER :: gth_potential 00600 00601 NULLIFY(atom) 00602 CALL create_atom_type(atom,error) 00603 00604 CALL get_atomic_kind(atomic_kind,z=z,zeff=zeff) 00605 NULLIFY(all_potential, gth_potential) 00606 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00607 all_potential=all_potential,& 00608 gth_potential=gth_potential) 00609 00610 IF(PRESENT(allelectron)) THEN 00611 IF(allelectron) THEN 00612 NULLIFY(gth_potential) 00613 zeff = z 00614 END IF 00615 END IF 00616 00617 CPPrecondition(ngto<=num_gto, cp_failure_level, routineP, error, failure) 00618 00619 IF ( ASSOCIATED(gth_potential) ) THEN 00620 ! PP calculation are non-relativistic 00621 relativistic = do_nonrel_atom 00622 ELSE 00623 ! AE calculations use DKH2 00624 relativistic = do_dkh2_atom 00625 END IF 00626 00627 atom%z = z 00628 CALL set_atom(atom,& 00629 pp_calc=ASSOCIATED(gth_potential),& 00630 method_type=do_rks_atom,& 00631 relativistic=relativistic,& 00632 coulomb_integral_type=do_numeric,& 00633 exchange_integral_type=do_numeric,& 00634 error=error) 00635 00636 ALLOCATE (potential,basis,integrals,STAT=ierr) 00637 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00638 00639 IF ( ASSOCIATED(gth_potential) ) THEN 00640 potential%ppot_type=gth_pseudo 00641 IF ( PRESENT(confine) ) THEN 00642 potential%confinement=confine 00643 ELSE 00644 potential%confinement=.TRUE. 00645 END IF 00646 potential%acon=0.1_dp 00647 potential%rcon=2._dp*ptable(z)%vdw_radius*bohr 00648 potential%ncon=2 00649 00650 CALL get_potential(gth_potential,& 00651 zeff=zeff,& 00652 elec_conf=ppeconf,& 00653 alpha_core_charge=ac,& 00654 nexp_ppl=ne,& 00655 cexp_ppl=ce,& 00656 lppnl=lm,& 00657 nprj_ppnl=nppnl,& 00658 alpha_ppnl=ap,& 00659 hprj_ppnl=hp) 00660 00661 potential%gth_pot%zion = zeff 00662 potential%gth_pot%rc = SQRT(0.5_dp/ac) 00663 potential%gth_pot%ncl = ne 00664 potential%gth_pot%cl(:) = 0._dp 00665 IF (ac > 0._dp) THEN 00666 DO i=1,ne 00667 potential%gth_pot%cl(i) = ce(i)/(2._dp*ac)**(i-1) 00668 END DO 00669 END IF 00670 !extended type 00671 potential%gth_pot%lpotextended = .FALSE. 00672 potential%gth_pot%lsdpot = .FALSE. 00673 potential%gth_pot%nlcc = .FALSE. 00674 potential%gth_pot%nexp_lpot = 0 00675 potential%gth_pot%nexp_lsd = 0 00676 potential%gth_pot%nexp_nlcc = 0 00677 CALL get_potential(gth_potential,& 00678 lpot_present=lpot_present,& 00679 lsd_present=lsd_present,& 00680 nlcc_present=nlcc_present) 00681 IF (lpot_present) THEN 00682 CALL get_potential(gth_potential,& 00683 nexp_lpot=nexp_lpot,& 00684 alpha_lpot=alpha_lpot,& 00685 nct_lpot=nct_lpot,& 00686 cval_lpot=cval_lpot) 00687 potential%gth_pot%lpotextended = .TRUE. 00688 potential%gth_pot%nexp_lpot = nexp_lpot 00689 potential%gth_pot%alpha_lpot(1:nexp_lpot) = SQRT(0.5_dp/alpha_lpot(1:nexp_lpot)) 00690 potential%gth_pot%nct_lpot(1:nexp_lpot) = nct_lpot(1:nexp_lpot) 00691 DO j=1,nexp_lpot 00692 ac = alpha_lpot(j) 00693 DO i=1,4 00694 potential%gth_pot%cval_lpot(i,j) = cval_lpot(i,j)/(2._dp*ac)**(i-1) 00695 END DO 00696 END DO 00697 END IF 00698 IF (lsd_present) THEN 00699 CALL get_potential(gth_potential,& 00700 nexp_lsd=nexp_lsd,& 00701 alpha_lsd=alpha_lsd,& 00702 nct_lsd=nct_lsd,& 00703 cval_lsd=cval_lsd) 00704 potential%gth_pot%lsdpot = .TRUE. 00705 potential%gth_pot%nexp_lsd = nexp_lsd 00706 potential%gth_pot%alpha_lsd(1:nexp_lsd) = SQRT(0.5_dp/alpha_lsd(1:nexp_lsd)) 00707 potential%gth_pot%nct_lsd(1:nexp_lsd) = nct_lsd(1:nexp_lsd) 00708 DO j=1,nexp_lpot 00709 ac = alpha_lsd(j) 00710 DO i=1,4 00711 potential%gth_pot%cval_lsd(i,j) = cval_lsd(i,j)/(2._dp*ac)**(i-1) 00712 END DO 00713 END DO 00714 END IF 00715 00716 ! nonlocal part 00717 potential%gth_pot%nl(:) = 0 00718 potential%gth_pot%rcnl(:) = 0._dp 00719 potential%gth_pot%hnl(:,:,:)= 0._dp 00720 DO l=0,lm 00721 n = nppnl(l) 00722 potential%gth_pot%nl(l) = n 00723 potential%gth_pot%rcnl(l) = SQRT(0.5_dp/ap(l)) 00724 potential%gth_pot%hnl(1:n,1:n,l)= hp(1:n,1:n,l) 00725 END DO 00726 00727 IF (nlcc_present) THEN 00728 CALL get_potential(gth_potential,& 00729 nexp_nlcc=nexp_nlcc,& 00730 alpha_nlcc=alpha_nlcc,& 00731 nct_nlcc=nct_nlcc,& 00732 cval_nlcc=cval_nlcc) 00733 potential%gth_pot%nlcc = .TRUE. 00734 potential%gth_pot%nexp_nlcc = nexp_nlcc 00735 potential%gth_pot%alpha_nlcc(1:nexp_nlcc) = alpha_nlcc(1:nexp_nlcc) 00736 potential%gth_pot%nct_nlcc(1:nexp_nlcc) = nct_nlcc(1:nexp_nlcc) 00737 potential%gth_pot%cval_nlcc(1:4,1:nexp_nlcc) = cval_nlcc(1:4,1:nexp_nlcc) 00738 END IF 00739 00740 CALL set_atom(atom,zcore=NINT(zeff),potential=potential,error=error) 00741 ELSE 00742 potential%ppot_type=no_pseudo 00743 IF ( PRESENT(confine) ) THEN 00744 potential%confinement=confine 00745 ELSE 00746 potential%confinement=.FALSE. 00747 END IF 00748 potential%acon=0.1_dp 00749 potential%rcon=2._dp*ptable(z)%vdw_radius*bohr 00750 potential%ncon=2 00751 CALL set_atom(atom,zcore=z,potential=potential,error=error) 00752 END IF 00753 00754 ! atomic grid 00755 NULLIFY(grid) 00756 ngp = 400 00757 quadtype = do_gapw_log 00758 CALL allocate_grid_atom(grid,error) 00759 CALL create_grid_atom(grid,ngp,1,1,0,quadtype) 00760 grid%nr = ngp 00761 basis%grid => grid 00762 00763 NULLIFY(basis%am,basis%cm,basis%as,basis%ns,basis%bf,basis%dbf) 00764 00765 ! fill in the basis data structures 00766 basis%eps_eig = 1.e-12_dp 00767 basis%basis_type = GTO_BASIS 00768 CALL Clementi_geobas(z,cval,aval,basis%nbas,starti,error) 00769 basis%nprim = basis%nbas 00770 m = MAXVAL(basis%nbas) 00771 ALLOCATE (basis%am(m,0:3),STAT=ierr) 00772 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00773 basis%am = 0._dp 00774 DO l=0,3 00775 DO i=1,basis%nbas(l) 00776 ll = i - 1 + starti(l) 00777 basis%am(i,l) = aval * cval**(ll) 00778 END DO 00779 END DO 00780 00781 basis%geometrical = .TRUE. 00782 basis%aval = aval 00783 basis%cval = cval 00784 basis%start = starti 00785 00786 ! initialize basis function on a radial grid 00787 nr = basis%grid%nr 00788 m = MAXVAL(basis%nbas) 00789 ALLOCATE (basis%bf(nr,m,0:3),STAT=ierr) 00790 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00791 ALLOCATE (basis%dbf(nr,m,0:3),STAT=ierr) 00792 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00793 basis%bf = 0._dp 00794 basis%dbf = 0._dp 00795 DO l=0,3 00796 DO i=1,basis%nbas(l) 00797 al = basis%am(i,l) 00798 DO k=1,nr 00799 rk = basis%grid%rad(k) 00800 ear = EXP(-al*basis%grid%rad(k)**2) 00801 basis%bf(k,i,l) = rk**l * ear 00802 basis%dbf(k,i,l) = ( REAL(l,dp)*rk**(l-1) - 2._dp*al*rk**(l+1) ) * ear 00803 END DO 00804 END DO 00805 END DO 00806 00807 CALL set_atom(atom,basis=basis,error=error) 00808 00809 ! optimization defaults 00810 atom%optimization%damping = 0.2_dp 00811 atom%optimization%eps_scf = 1.e-6_dp 00812 atom%optimization%eps_diis = 100._dp 00813 atom%optimization%max_iter = 50 00814 atom%optimization%n_diis = 5 00815 00816 ! electronic state 00817 nelem = 0 00818 ncore = 0 00819 ncalc = 0 00820 NULLIFY(econf) 00821 IF ( ASSOCIATED(gth_potential) ) THEN 00822 CALL get_potential(gth_potential,elec_conf=econf) 00823 econfx = 0 00824 econfx(0:SIZE(econf)-1) = econf 00825 IF ( SUM(econf) >= 0 ) THEN 00826 DO l=0,3 00827 ll = 2*(2*l+1) 00828 nn = ptable(z)%e_conv(l)-econfx(l) 00829 ii = 0 00830 DO 00831 ii = ii + 1 00832 IF(nn <= ll) THEN 00833 ncore(l,ii) = nn 00834 EXIT 00835 ELSE 00836 ncore(l,ii) = ll 00837 nn = nn - ll 00838 END IF 00839 END DO 00840 END DO 00841 DO l=0,3 00842 ll = 2*(2*l+1) 00843 nn = ptable(z)%e_conv(l) 00844 ii = 0 00845 DO 00846 ii = ii + 1 00847 IF(nn <= ll) THEN 00848 nelem(l,ii) = nn 00849 EXIT 00850 ELSE 00851 nelem(l,ii) = ll 00852 nn = nn - ll 00853 END IF 00854 END DO 00855 END DO 00856 ncalc = nelem - ncore 00857 ELSE 00858 ncore = 0 00859 ncalc = 0 00860 DO l=0,3 00861 ll = 2*(2*l+1) 00862 nn = ABS(econfx(l)) 00863 ii = 0 00864 DO 00865 ii = ii + 1 00866 IF(nn <= ll) THEN 00867 ncalc(l,ii) = -nn 00868 EXIT 00869 ELSE 00870 ncalc(l,ii) = -ll 00871 nn = nn - ll 00872 END IF 00873 END DO 00874 END DO 00875 END IF 00876 ELSE 00877 DO l=0,3 00878 ll = 2*(2*l+1) 00879 nn = ptable(z)%e_conv(l) 00880 ii = 0 00881 DO 00882 ii = ii + 1 00883 IF(nn <= ll) THEN 00884 nelem(l,ii) = nn 00885 EXIT 00886 ELSE 00887 nelem(l,ii) = ll 00888 nn = nn - ll 00889 END IF 00890 END DO 00891 END DO 00892 ncalc = nelem - ncore 00893 END IF 00894 00895 IF ( atomic_kind%ghost ) THEN 00896 nelem = 0 00897 ncore = 0 00898 ncalc = 0 00899 END IF 00900 00901 ALLOCATE (atom%state,STAT=ierr) 00902 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00903 00904 atom%state%core = 0._dp 00905 atom%state%core(0:3,1:7) = REAL(ncore(0:3,1:7),dp) 00906 atom%state%occ = 0._dp 00907 atom%state%occ(0:3,1:7) = REAL(ncalc(0:3,1:7),dp) 00908 atom%state%occupation = 0._dp 00909 atom%state%multiplicity = -1 00910 DO l=0,3 00911 k = 0 00912 DO i=1,7 00913 IF ( ncalc(l,i) > 0 ) THEN 00914 k = k + 1 00915 atom%state%occupation(l,k) = REAL(ncalc(l,i),dp) 00916 END IF 00917 END DO 00918 END DO 00919 00920 atom%state%maxl_occ = get_maxl_occ(atom%state%occupation) 00921 atom%state%maxn_occ = get_maxn_occ(atom%state%occupation) 00922 atom%state%maxl_calc = atom%state%maxl_occ 00923 atom%state%maxn_calc = atom%state%maxn_occ 00924 00925 ! calculate integrals 00926 ! general integrals 00927 CALL atom_int_setup(integrals,basis,potential=atom%potential,& 00928 eri_coulomb=(atom%coulomb_integral_type==do_analytic),& 00929 eri_exchange=(atom%exchange_integral_type==do_analytic),error=error) 00930 ! potential 00931 CALL atom_ppint_setup(integrals,basis,potential=atom%potential,error=error) 00932 ! relativistic correction terms 00933 NULLIFY(integrals%tzora,integrals%hdkh) 00934 CALL atom_relint_setup(integrals,basis,atom%relativistic,zcore=REAL(atom%zcore,dp),error=error) 00935 CALL set_atom(atom,integrals=integrals,error=error) 00936 00937 NULLIFY(orbitals) 00938 mo = MAXVAL(atom%state%maxn_calc) 00939 mb = MAXVAL(atom%basis%nbas) 00940 CALL create_atom_orbs(orbitals,mb,mo,error) 00941 CALL set_atom(atom,orbitals=orbitals,error=error) 00942 00943 IF(PRESENT(iunit)) THEN 00944 CALL calculate_atom(atom,iunit,error=error) 00945 CALL atom_fit_density (atom,ngto,0,iunit,results=results,error=error) 00946 ELSE 00947 CALL calculate_atom(atom,-1,error=error) 00948 CALL atom_fit_density (atom,ngto,0,-1,results=results,error=error) 00949 END IF 00950 00951 xx = results(1) 00952 cc = results(2) 00953 DO i=1,ngto 00954 density(i,1) = xx * cc**i 00955 density(i,2) = results(2+i) 00956 END DO 00957 00958 ! clean up 00959 CALL atom_int_release(integrals,error) 00960 CALL atom_ppint_release(integrals,error) 00961 CALL atom_relint_release(integrals,error) 00962 CALL release_atom_basis(basis,error) 00963 CALL release_atom_potential(potential,error) 00964 CALL release_atom_type(atom,error) 00965 00966 DEALLOCATE (potential,basis,integrals,STAT=ierr) 00967 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00968 00969 END SUBROUTINE calculate_atomic_density 00970 00971 ! ***************************************************************************** 00972 00973 END MODULE atom_kind_orbitals
1.7.3