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

@@ -411,42 +411,44 @@ SUBROUTINE atom_grb_construction(atom_info, atom_section, iw)
       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.)
       CALL atom_density(rho(:), atom%orbitals%pmat, atom%basis, maxl, typ="RHO")
       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 @@ SUBROUTINE atom_grb_construction(atom_info, atom_section, iw)
          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

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
+    &GTH_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

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