From 9ec4c7d1d7cadce45c5566cc8db765df7b7152bb Mon Sep 17 00:00:00 2001 From: Juerg Hutter Date: Thu, 10 Dec 2020 15:23:09 +0100 Subject: [PATCH] Update GRB Basis: allow for 0 polarization functions (#1214) * xTB pairpotential (A. Hehn) * Adjust regtest * Adjust regtests * Allow GRB Basis without polarization functions --- src/atom_grb.F | 73 +++++++++++++++------------- tests/ATOM/regtest-pseudo/C_grb2.inp | 47 ++++++++++++++++++ tests/ATOM/regtest-pseudo/TEST_FILES | 1 + 3 files changed, 88 insertions(+), 33 deletions(-) create mode 100644 tests/ATOM/regtest-pseudo/C_grb2.inp diff --git a/src/atom_grb.F b/src/atom_grb.F index 388225630f..615bf95ccf 100644 --- a/src/atom_grb.F +++ b/src/atom_grb.F @@ -411,16 +411,6 @@ CONTAINS CALL deallocate_orbital_pointers CALL deallocate_spherical_harmonics - ! generate polarization sets - IF (iw > 0) THEN - WRITE (iw, '(/,A)') " Polarization basis set " - END IF - maxl = atom%state%maxl_occ - CALL section_vals_val_get(grb_section, "NUM_GTO_POLARIZATION", i_val=num_gto) - CPASSERT(num_gto > 0) - num_pol = num_gto - ALLOCATE (pbasis(num_gto, num_gto, 0:7), alp(num_gto)) - pbasis = 0.0_dp ! get density maximum ALLOCATE (rho(basis%grid%nr)) CALL calculate_atom(atom, iw=0, noguess=.TRUE.) @@ -428,25 +418,37 @@ CONTAINS n = SUM(MAXLOC(rho(:))) rmax = basis%grid%rad(n) DEALLOCATE (rho) - ! optimize exponents - lval = maxl + 1 - zval = SQRT(REAL(2*lval + 2, dp))*REAL(lval + 1, dp)/(2._dp*rmax) - aval = atom%basis%am(1, 0) - cval = 2.5_dp - rconf = atom%potential%scon - CALL atom_fit_pol(zval, rconf, lval, aval, cval, num_gto, iw, powell_section) - ! calculate contractions - DO i = 1, num_gto - alp(i) = aval*cval**(i - 1) - END DO - ALLOCATE (rho(num_gto)) - DO l = maxl + 1, MIN(maxl + num_gto, 7) - zval = SQRT(REAL(2*l + 2, dp))*REAL(l + 1, dp)/(2._dp*rmax) - CALL hydrogenic(zval, rconf, l, alp, num_gto, rho, pbasis(:, :, l)) - IF (iw > 0) WRITE (iw, '(T5,A,i5,T66,A,F10.4)') & - " Polarization basis set contraction for lval=", l, "zval=", zval - END DO - DEALLOCATE (rho) + + ! generate polarization sets + maxl = atom%state%maxl_occ + CALL section_vals_val_get(grb_section, "NUM_GTO_POLARIZATION", i_val=num_gto) + IF (num_gto > 0) THEN + IF (iw > 0) THEN + WRITE (iw, '(/,A)') " Polarization basis set " + END IF + num_pol = num_gto + ALLOCATE (pbasis(num_gto, num_gto, 0:7), alp(num_gto)) + pbasis = 0.0_dp + ! optimize exponents + lval = maxl + 1 + zval = SQRT(REAL(2*lval + 2, dp))*REAL(lval + 1, dp)/(2._dp*rmax) + aval = atom%basis%am(1, 0) + cval = 2.5_dp + rconf = atom%potential%scon + CALL atom_fit_pol(zval, rconf, lval, aval, cval, num_gto, iw, powell_section) + ! calculate contractions + DO i = 1, num_gto + alp(i) = aval*cval**(i - 1) + END DO + ALLOCATE (rho(num_gto)) + DO l = maxl + 1, MIN(maxl + num_gto, 7) + zval = SQRT(REAL(2*l + 2, dp))*REAL(l + 1, dp)/(2._dp*rmax) + CALL hydrogenic(zval, rconf, l, alp, num_gto, rho, pbasis(:, :, l)) + IF (iw > 0) WRITE (iw, '(T5,A,i5,T66,A,F10.4)') & + " Polarization basis set contraction for lval=", l, "zval=", zval + END DO + DEALLOCATE (rho) + END IF ! generate valence expansion sets maxl = atom%state%maxl_occ @@ -674,14 +676,19 @@ CONTAINS CALL grb_print_basis(header=basline, nprim=num_pol, nbas=nbas, al=alp, gcc=pbasis, iunit=iunit) END DO ! extension set - basline(1) = "" - WRITE (basline(1), "(T2,A,T5,A)") ADJUSTL(ptable(atom_ref%z)%symbol), TRIM(ADJUSTL(basname))//"-EXT" - CALL grb_print_basis(header=basline, nprim=next_prim(0), nbas=next_bas, al=ale, gcc=ebasis, iunit=iunit) + IF (SUM(next_bas) > 0) THEN + basline(1) = "" + WRITE (basline(1), "(T2,A,T5,A)") ADJUSTL(ptable(atom_ref%z)%symbol), TRIM(ADJUSTL(basname))//"-EXT" + CALL grb_print_basis(header=basline, nprim=next_prim(0), nbas=next_bas, al=ale, gcc=ebasis, iunit=iunit) + END IF ! CALL close_file(unit_number=iunit) ! clean up - DEALLOCATE (wfn, rbasis, qbasis, ebasis, pbasis, ale, alp) + IF (ALLOCATED(pbasis)) DEALLOCATE (pbasis) + IF (ALLOCATED(alp)) DEALLOCATE (alp) + IF (ALLOCATED(ebasis)) DEALLOCATE (ebasis) + DEALLOCATE (wfn, rbasis, qbasis, ale) DO ider = 0, 10 IF (ASSOCIATED(vbasis(ider)%basis)) THEN diff --git a/tests/ATOM/regtest-pseudo/C_grb2.inp b/tests/ATOM/regtest-pseudo/C_grb2.inp new file mode 100644 index 0000000000..d7bcc89f98 --- /dev/null +++ b/tests/ATOM/regtest-pseudo/C_grb2.inp @@ -0,0 +1,47 @@ +&GLOBAL + PROGRAM_NAME ATOM +&END GLOBAL +&ATOM + ELEMENT C + + RUN_TYPE ENERGY + + ELECTRON_CONFIGURATION CORE 2s2 2p2 + CORE [He] + MAX_ANGULAR_MOMENTUM 3 + + &METHOD + METHOD_TYPE KOHN-SHAM + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END METHOD + + &PRINT + &GEOMETRICAL_RESPONSE_BASIS + NUM_GTO_CORE 6 + NUM_GTO_EXTENDED 4 + EXTENSION_BASIS 0 0 + CONFINEMENT 8.0 + NAME_BODY GRB-q4 + &END + &END + + &POTENTIAL + PSEUDO_TYPE GTH + >H_POTENTIAL + 2 2 0 0 + 0.33842037096913 2 -8.80041712690571 1.33249221738633 + 1 + 0.30268766321861 1 9.60946227253083 + &END + CONFINEMENT_TYPE BARRIER + CONFINEMENT 200. 4.0 12.0 + &END POTENTIAL + &POWELL + ACCURACY 1.e-14 + STEP_SIZE 0.20 + &END + +&END ATOM diff --git a/tests/ATOM/regtest-pseudo/TEST_FILES b/tests/ATOM/regtest-pseudo/TEST_FILES index 8aa32a9492..d6482e2dae 100644 --- a/tests/ATOM/regtest-pseudo/TEST_FILES +++ b/tests/ATOM/regtest-pseudo/TEST_FILES @@ -20,5 +20,6 @@ C_basis1.inp 35 1.e-12 C_basis2.inp 35 1.e-12 -37.728835366125 C_basis3.inp 35 1.e-12 -13.923782735900 C_grb.inp 35 1.e-12 -5.360945665362 +C_grb2.inp 35 1.e-12 -5.360945665362 Li_NLCC.inp 61 1.e-12 0.5527955347 #EOF