Bug fix for D4 code (definition of atomic type) (cp2k#3685)

hfp · Sep 20, 2024 · 7f5066a · 7f5066a
1 parent 46309a9
commit 7f5066a
Show file tree

Hide file tree

Showing 3 changed files with 22 additions and 8 deletions.
diff --git a/src/qs_dispersion_d4.F b/src/qs_dispersion_d4.F
@@ -88,7 +88,7 @@ SUBROUTINE calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calcula
 
       INTEGER                                            :: atoma, handle, i, iatom, ikind, mref, &
                                                             natom
-      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of
+      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, atomtype, kind_of
       INTEGER, DIMENSION(3)                              :: periodic
       LOGICAL                                            :: grad, use_virial
       LOGICAL, DIMENSION(3)                              :: lperiod
@@ -99,6 +99,7 @@ SUBROUTINE calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calcula
       TYPE(mp_para_env_type), POINTER                    :: para_env
       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
+      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
       TYPE(virial_type), POINTER                         :: virial
 
       CLASS(damping_param), ALLOCATABLE                  :: param
@@ -124,9 +125,12 @@ SUBROUTINE calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calcula
 
       !get information about particles
       natom = SIZE(particle_set)
-      ALLOCATE (xyz(3, natom))
+      ALLOCATE (xyz(3, natom), atomtype(natom))
+      CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
       DO iatom = 1, natom
          xyz(:, iatom) = particle_set(iatom)%r(:)
+         ikind = kind_of(iatom)
+         CALL get_qs_kind(qs_kind_set(ikind), zatom=atomtype(iatom))
       END DO
 
       !get information about cell / lattice
@@ -136,7 +140,7 @@ SUBROUTINE calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calcula
       lperiod(3) = periodic(3) == 1
 
       !prepare for the call to the dispersion function
-      CALL new(mol, kind_of, xyz, lattice=cell%hmat, periodic=lperiod)
+      CALL new(mol, atomtype, xyz, lattice=cell%hmat, periodic=lperiod)
       CALL new_d4_model(disp, mol)
       CALL get_rational_damping(dispersion_env%ref_functional, param, s9=dispersion_env%s9)
 
@@ -284,7 +288,7 @@ SUBROUTINE calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calcula
 
       END IF
 
-      DEALLOCATE (xyz)
+      DEALLOCATE (xyz, atomtype)
 
       CALL timestop(handle)
 

diff --git a/src/qs_kind_types.F b/src/qs_kind_types.F
@@ -373,6 +373,7 @@ END SUBROUTINE deallocate_qs_kind_set
 !> \param dftb_parameter ...
 !> \param xtb_parameter ...
 !> \param dftb3_param ...
+!> \param zatom ...
 !> \param zeff ...
 !> \param elec_conf ...
 !> \param mao ...
@@ -435,7 +436,7 @@ SUBROUTINE get_qs_kind(qs_kind, &
                           basis_set, basis_type, ncgf, nsgf, &
                           all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, &
                           se_parameter, dftb_parameter, xtb_parameter, &
-                          dftb3_param, zeff, elec_conf, mao, lmax_dftb, &
+                          dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, &
                           alpha_core_charge, ccore_charge, core_charge, core_charge_radius, &
                           paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, &
                           covalent_radius, vdw_radius, &
@@ -461,7 +462,9 @@ SUBROUTINE get_qs_kind(qs_kind, &
       TYPE(semi_empirical_type), OPTIONAL, POINTER       :: se_parameter
       TYPE(qs_dftb_atom_type), OPTIONAL, POINTER         :: dftb_parameter
       TYPE(xtb_atom_type), OPTIONAL, POINTER             :: xtb_parameter
-      REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: dftb3_param, zeff
+      REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: dftb3_param
+      INTEGER, INTENT(OUT), OPTIONAL                     :: zatom
+      REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: zeff
       INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: elec_conf
       INTEGER, INTENT(OUT), OPTIONAL                     :: mao, lmax_dftb
       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: alpha_core_charge, ccore_charge, &
@@ -506,6 +509,7 @@ SUBROUTINE get_qs_kind(qs_kind, &
 
       CHARACTER(LEN=default_string_length)               :: my_basis_type
       INTEGER                                            :: l
+      LOGICAL                                            :: found
       TYPE(gto_basis_set_type), POINTER                  :: tmp_basis_set
 
       ! Retrieve basis set from the kind container
@@ -621,6 +625,12 @@ SUBROUTINE get_qs_kind(qs_kind, &
          END IF
       END IF
 
+      IF (PRESENT(zatom)) THEN
+         ! Retrieve information on element
+         CALL get_ptable_info(qs_kind%element_symbol, ielement=zatom, found=found)
+         CPASSERT(found)
+      END IF
+
       IF (PRESENT(zeff)) THEN
          IF (ASSOCIATED(qs_kind%all_potential)) THEN
             CALL get_potential(potential=qs_kind%all_potential, zeff=zeff)

diff --git a/tests/QS/regtest-dft-vdw-corr-4/TEST_FILES b/tests/QS/regtest-dft-vdw-corr-4/TEST_FILES
@@ -3,7 +3,7 @@
 # e.g. 0 means do not compare anything, running is enough
 #      1 compares the last total energy in the file
 #      for details see cp2k/tools/do_regtest
-pbe_dftd4.inp                                          33    1.0E-14              -0.00141644869634
+pbe_dftd4.inp                                          33    1.0E-14              -0.00283102230260
 pbe_dftd4_force.inp                                    72    1.0E-07                     0.00007217 
-pbe_dftd4_stress.inp                                   31    1.0E-07             -5.16289312880E-03
+pbe_dftd4_stress.inp                                   31    1.0E-07             -2.14003785359E-02
 #EOF