Guess on parameter values for orig. code

cp2k · Mar 19, 2019 · 2c6ccf5 · 2c6ccf5
1 parent 70c07f1
commit 2c6ccf5
Show file tree

Hide file tree

Showing 5 changed files with 19 additions and 10 deletions.
diff --git a/src/qs_environment.F b/src/qs_environment.F
@@ -666,9 +666,6 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
             CALL init_potential(qs_kind%all_potential, itype="BARE", &
                                 zeff=qs_kind%xtb_parameter%zeff, zeff_correction=zeff_correction)
             qs_kind%xtb_parameter%zeff = qs_kind%xtb_parameter%zeff-zeff_correction
-!qxtb
-! set econf ?
-!qxtb
          END DO
          ! check for Ewald
          IF (xtb_control%do_ewald) THEN

diff --git a/src/xtb_coulomb.F b/src/xtb_coulomb.F
@@ -463,8 +463,8 @@ SUBROUTINE build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
          CALL get_xtb_atom_param(xtb_kind, xgamma=xgamma(ikind), zeff=zeffk(ikind))
       END DO
       ! Diagonal 3rd order correction (DFTB3)
-      CALL build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeffk, &
-                                calculate_forces, just_energy)
+!deb  CALL build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeffk, &
+!deb                            calculate_forces, just_energy)
       DEALLOCATE (zeffk, xgamma)
 
       DEALLOCATE (gmcharge, gchrg, atom_of_kind, kind_of)

diff --git a/src/xtb_matrices.F b/src/xtb_matrices.F
@@ -385,15 +385,16 @@ SUBROUTINE build_xtb_matrices(qs_env, para_env, calculate_forces)
          END IF
 
          ! Calculate Pi = Pia * Pib (Eq. 11)
+         ! factor 0.01 from original Grimme code, not in paper?
          rcovab = rcova+rcovb
          rrab = SQRT(dr/rcovab)
          DO i = 1, nsa
-            pia(i) = 1._dp+kpolya(i)*rrab
-            dpia(i) = 0.5_dp*kpolya(i)/rrab
+            pia(i) = 1._dp+0.01_dp*kpolya(i)*rrab
+            dpia(i) = 0.5_dp*0.01_dp*kpolya(i)/rrab
          END DO
          DO i = 1, nsb
-            pib(i) = 1._dp+kpolyb(i)*rrab
-            dpib(i) = 0.5_dp*kpolyb(i)/rrab
+            pib(i) = 1._dp+0.01_dp*kpolyb(i)*rrab
+            dpib(i) = 0.5_dp*0.01_dp*kpolyb(i)/rrab
          END DO
 
          IF (iatom == jatom .AND. dr < 0.001_dp) THEN
@@ -416,7 +417,9 @@ SUBROUTINE build_xtb_matrices(qs_env, para_env, calculate_forces)
             ! get KAB
             kab = kabset(ikind, jkind)
             ! get Fen = (1+ken*deltaEN^2)
-            fen = 1.0_dp+ken*(ena-enb)**2
+            ! original code reads: (1-0.01*ken*den^2)
+            ! fen = 1.0_dp+ken*(ena-enb)**2
+            fen = 1.0_dp - 0.01_dp*ken*(ena-enb)**2
             DO j = 1, natorb_b
                lb = laob(j)
                nb = naob(j)

diff --git a/src/xtb_parameters.F b/src/xtb_parameters.F
@@ -151,6 +151,9 @@ SUBROUTINE xtb_parameters_init(param, element_symbol, parameter_file_path, param
             found = .TRUE.
             CALL parser_get_object(parser, param%eta)
             CALL parser_get_object(parser, param%xgamma)
+!xtb
+            param%xgamma = param%xgamma/evolt
+!xtb
             CALL parser_get_object(parser, param%alpha)
             CALL parser_get_object(parser, param%zneff)
             DO i = 1, 5

diff --git a/tests/xTB/regtest-1/ch2o.inp b/tests/xTB/regtest-1/ch2o.inp
@@ -20,6 +20,12 @@
       SCF_GUESS MOPAC
       MAX_SCF  20
     &END SCF
+    &POISSON
+      &EWALD
+       EWALD_TYPE SPME
+       GMAX 25
+      &END EWALD
+    &END POISSON
   &END DFT
   &SUBSYS
     &CELL