Example for Feynman-Hibbs potential

Pretty Set regtest values Adjust regtests Update regtest Update regtest tolerance
cp2k · Jul 9, 2021 · 50617a1 · 50617a1
1 parent f51128f
commit 50617a1
Show file tree

Hide file tree

Showing 7 changed files with 108 additions and 20 deletions.
diff --git a/src/particle_methods.F b/src/particle_methods.F
@@ -123,7 +123,7 @@ SUBROUTINE get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf
          DO iparticle = 1, nparticle
             CALL get_atomic_kind(particle_set(iparticle)%atomic_kind, kind_number=ikind)
             IF (PRESENT(basis)) THEN
-               IF(ASSOCIATED(basis(ikind)%gto_basis_set)) THEN
+               IF (ASSOCIATED(basis(ikind)%gto_basis_set)) THEN
                   CALL get_gto_basis_set(gto_basis_set=basis(ikind)%gto_basis_set, nsgf=ns)
                ELSE
                   ns = 0

diff --git a/src/qs_environment.F b/src/qs_environment.F
@@ -711,7 +711,7 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
                CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, "ORB")
             ELSE
                CALL set_qs_kind(qs_kind, ghost=.TRUE.)
-               IF(ASSOCIATED(qs_kind%all_potential)) THEN
+               IF (ASSOCIATED(qs_kind%all_potential)) THEN
                   DEALLOCATE (qs_kind%all_potential%elec_conf)
                   DEALLOCATE (qs_kind%all_potential)
                END IF
@@ -960,7 +960,7 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
          ! cutoff radius
          DO ikind = 1, nkind
             qs_kind => qs_kind_set(ikind)
-            IF(qs_kind%xtb_parameter%defined) THEN
+            IF (qs_kind%xtb_parameter%defined) THEN
                CALL get_qs_kind(qs_kind, basis_set=tmp_basis_set)
                CALL get_gto_basis_set(tmp_basis_set, kind_radius=qs_kind%xtb_parameter%rcut)
             ELSE

diff --git a/src/xtb_matrices.F b/src/xtb_matrices.F
@@ -1391,10 +1391,10 @@ SUBROUTINE nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_
                                ' Error '//TRIM(this_error)//' in computing numerical derivatives larger then'// &
                                TRIM(def_error)//' .')
                END IF
- 
+
                atom_i = atom_of_kind(iatom)
                atom_j = atom_of_kind(jatom)
- 
+
                fij(1:3) = dedf*rij(1:3)/pot%set(kk)%gp%values(1)
 
                force(ikind)%repulsive(:, atom_i) = force(ikind)%repulsive(:, atom_i) - fij(:)
@@ -1420,7 +1420,7 @@ SUBROUTINE nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_
       ENDIF
 
       CALL timestop(handle)
- 
+
    END SUBROUTINE nonbonded_correction
 ! **************************************************************************************************
 !> \brief ...

diff --git a/src/xtb_parameters.F b/src/xtb_parameters.F
@@ -508,7 +508,7 @@ SUBROUTINE init_xtb_basis(param, gto_basis_set, ngauss)
             NULLIFY (sto_basis_set)
             CALL allocate_sto_basis_set(sto_basis_set)
             nshell = param%nshell
-   
+
             ALLOCATE (symbol(1:nshell))
             symbol = ""
             DO i = 1, nshell
@@ -540,7 +540,7 @@ SUBROUTINE init_xtb_basis(param, gto_basis_set, ngauss)
             CALL deallocate_sto_basis_set(sto_basis_set)
             DEALLOCATE (symbol, nq, lq, zet)
          END IF
-   
+
       ELSE
          CPABORT("The pointer param is not associated")
       END IF

diff --git a/tests/xTB/regtest-1/TEST_FILES b/tests/xTB/regtest-1/TEST_FILES
@@ -17,4 +17,7 @@ h2o-atprop0.inp                            1    1.0E-12            -187.45499307
 si_geo.inp                                 1    1.0E-12             -14.55678348978414
 si_kp.inp                                  1    1.0E-12             -14.75431646565682
 h2o_dimer.inp                              1    1.0E-12             -11.54506384130837
+ghost.inp                                  1    1.0E-12              -1.03458109975585
+ghost_md.inp                               1    1.0E-08              -2.072098719     
+fhpot.inp                                  1    1.0E-07              -2.071697129     
 #EOF
diff --git a/tests/xTB/regtest-1/fhpot.inp b/tests/xTB/regtest-1/fhpot.inp
@@ -0,0 +1,92 @@
+&FORCE_EVAL
+  &DFT
+    &QS
+      METHOD xTB
+      &xTB
+         DO_EWALD T
+         DO_NONBONDED T
+         &NONBONDED
+            # Feynman-Hibbs-Potential
+            # E : LJ epsilon value [Hartree]
+            # S : LJ sigma value [Bohr]
+            # K : hbar**2 / 24*k (mu: reduced mass, k: Boltzmann constant)
+            # T : Temperature [Kelvin]
+            # mu: reduced mass [atomic units]
+            &GENPOT
+               ATOMS X X  
+               FUNCTION 24*E*( K/(mu*T*S*S)*( 22*(S/r)**14 - 5*(S/r)**8) + (K/(mu*T*S*S))**2*( 2002*(S/r)**16 - 140*(S/r)**10) )
+               PARAMETERS E S K T mu
+               VALUES  8.2337E-5 4.830 7.217825 20.0 0.5
+               VARIABLES r
+            &END GENPOT
+        &END NONBONDED
+      &END XTB
+    &END QS
+    &SCF
+      SCF_GUESS MOPAC
+      MAX_SCF  100
+      EPS_SCF 1.e-8
+    &END SCF
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC 10.0 10.0 10.0
+    &END CELL
+    &COORD
+     UNIT bohr
+     H     0.000000    0.000000   20.000000 H2
+     H     1.360000    0.000000   20.000000 H2
+     X     0.680000    0.000000   20.000000 H2
+     H    -4.000000    0.000000   20.000000 H2
+     H    -4.000000    1.360000   20.000000 H2
+     X    -4.000000    0.680000   20.000000 H2
+    &END COORD
+    &COLVAR
+      &DISTANCE
+        ATOMS 1 2
+      &END DISTANCE
+    &END COLVAR
+  &END SUBSYS
+&END FORCE_EVAL
+&GLOBAL
+  PROJECT H2
+  RUN_TYPE MD
+  PRINT_LEVEL MEDIUM
+&END GLOBAL
+&MOTION
+  &MD
+    ENSEMBLE NVE
+    STEPS 10
+    TIMESTEP 1.0
+    TEMPERATURE 100
+  &END MD
+  &CONSTRAINT
+    &COLLECTIVE OFF
+      MOLNAME H2
+      COLVAR 1
+      INTERMOLECULAR FALSE
+      TARGET [BOHR] 1.40
+    &END COLLECTIVE
+    &VIRTUAL_SITE
+      MOLNAME H2
+      INTERMOLECULAR FALSE
+      ATOMS 3 2 1 2
+      PARAMETERS 0.250000000 0.250000000
+    &END VIRTUAL_SITE
+  &END CONSTRAINT
+  &PRINT
+    &FORCES ON
+      &EACH
+        MD 2
+      &END EACH
+      FORMAT XYZ
+    &END FORCES
+    &VELOCITIES ON
+      &EACH
+        MD 2
+      &END EACH
+      FORMAT XYZ
+    &END VELOCITIES
+  &END PRINT
+&END MOTION
+
diff --git a/tests/xTB/regtest-1/ghost_md.inp b/tests/xTB/regtest-1/ghost_md.inp
@@ -7,17 +7,10 @@
          DO_NONBONDED T
          &NONBONDED
             &GENPOT
-               ATOMS X X  
-               FUNCTION Aparam*exp(-Bparam*r)-Cparam/r**8
-               PARAMETERS Aparam Bparam Cparam
-               VALUES  1.0 1.441 6703.192
-               VARIABLES r
-            &END GENPOT
-            &GENPOT
-               ATOMS X X  
-               FUNCTION Aparam*exp(-Bparam*r)-Cparam/r**8
-               PARAMETERS Aparam Bparam Cparam
-               VALUES  1.0 1.441 6703.192
+               ATOMS X X
+               FUNCTION 4*E*( (S/r)**8 - (S/r)**6 )
+               PARAMETERS E S
+               VALUES  8.2337E-5 4.830
                VARIABLES r
             &END GENPOT
         &END NONBONDED
@@ -63,7 +56,7 @@
     TEMPERATURE 100
   &END MD
   &CONSTRAINT
-    &COLLECTIVE OFF
+    &COLLECTIVE 
       MOLNAME H2
       COLVAR 1
       INTERMOLECULAR FALSE