Enable different BD values for BMHFTD potential

cp2k · Oct 19, 2021 · d062abb · d062abb
1 parent f0c271f
commit d062abb
Show file tree

Hide file tree

Showing 5 changed files with 63 additions and 42 deletions.
diff --git a/src/force_fields_input.F b/src/force_fields_input.F
@@ -532,7 +532,7 @@ SUBROUTINE set_IPBV_ff(at1, at2, ipbv)
          ipbv%a(15) = -79705131323.98_dp
       ELSE
          CPABORT("IPBV only for WATER")
-      ENDIF
+      END IF
    END SUBROUTINE set_IPBV_ff
 
 ! **************************************************************************************************
@@ -562,7 +562,7 @@ SUBROUTINE set_BMHFT_ff(at1, at2, ft)
          ft%d = cp_unit_to_cp2k(145.427_dp, "eV*angstrom^8")
       ELSE
          CPABORT("BMHFT only for NaCl")
-      ENDIF
+      END IF
 
    END SUBROUTINE set_BMHFT_ff
 
@@ -911,8 +911,12 @@ SUBROUTINE read_bmhftd_section(nonbonded, section, start)
          DIMENSION(:), POINTER                           :: atm_names
       INTEGER                                            :: i, isec, n_items, n_rep
       REAL(KIND=dp)                                      :: rcut
+      REAL(KIND=dp), DIMENSION(:), POINTER               :: bd_vals
+
+      NULLIFY (bd_vals)
 
       CALL section_vals_get(section, n_repetition=n_items)
+
       DO isec = 1, n_items
          CALL cite_reference(Tosi1964a)
          CALL cite_reference(Tosi1964b)
@@ -933,9 +937,21 @@ SUBROUTINE read_bmhftd_section(nonbonded, section, start)
                                       r_val=nonbonded%pot(start + isec)%pot%set(1)%ftd%c)
             CALL section_vals_val_get(section, "D", i_rep_section=isec, &
                                       r_val=nonbonded%pot(start + isec)%pot%set(1)%ftd%d)
-            CALL section_vals_val_get(section, "BD", i_rep_section=isec, &
-                                      r_val=nonbonded%pot(start + isec)%pot%set(1)%ftd%bd)
-
+            CALL section_vals_val_get(section, "BD", i_rep_section=isec, r_vals=bd_vals)
+            IF (ASSOCIATED(bd_vals)) THEN
+               SELECT CASE (SIZE(bd_vals))
+               CASE (0)
+                  CPABORT("No values specified for parameter BD in section &BMHFTD")
+               CASE (1)
+                  nonbonded%pot(start + isec)%pot%set(1)%ftd%bd(1:2) = bd_vals(1)
+               CASE (2)
+                  nonbonded%pot(start + isec)%pot%set(1)%ftd%bd(1:2) = bd_vals(1:2)
+               CASE DEFAULT
+                  CPABORT("Too many values specified for parameter BD in section &BMHFTD")
+               END SELECT
+            ELSE
+               CPABORT("Parameter BD in section &BMHFTD was not specified")
+            END IF
          ELSE
             CALL section_vals_val_get(section, "MAP_ATOMS", i_rep_section=isec, c_vals=atm_names)
             map_atoms = atm_names
@@ -1651,7 +1667,7 @@ SUBROUTINE read_apol_section(apol_atm, apol, damping_list, section, &
 
       IF (n_damp > 0) THEN
          ALLOCATE (damping_list(1:n_damp))
-      ENDIF
+      END IF
 
 ! *** Reads DIPOLE sections *****
       start_damp = 0
@@ -1989,11 +2005,11 @@ SUBROUTINE read_torsions_section(torsion_kind, torsion_a, torsion_b, torsion_c,
             IF (torsion_phi0(start + isec) .NE. 0.0_dp) THEN
                CALL cp_warn(__LOCATION__, "PHI0 parameter was non-zero "// &
                             "for an OPLS-type TORSION.  It will be ignored.")
-            ENDIF
+            END IF
             IF (MODULO(torsion_m(start + isec), 2) .EQ. 0) THEN
                ! For even M, negate the cosine using a Pi phase factor
                torsion_phi0(start + isec) = pi
-            ENDIF
+            END IF
             ! the K parameter appears as K/2 in the OPLS parameterisation
             torsion_k(start + isec) = torsion_k(start + isec)*0.5_dp
          END IF

diff --git a/src/input_cp2k_mm.F b/src/input_cp2k_mm.F
@@ -1842,37 +1842,39 @@ SUBROUTINE create_BMHFTD_section(section)
       CALL keyword_release(keyword)
 
       CALL keyword_create(keyword, __LOCATION__, name="A", &
-                          description="Defines the A parameter of the dispersion-damped Fumi-Tosi Potential", &
+                          description="Defines the A parameter of the dispersion-damped Fumi-Tosi potential", &
                           usage="A {real}", type_of_var=real_t, &
                           n_var=1, unit_str="hartree")
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
       CALL keyword_create(keyword, __LOCATION__, name="B", &
-                          description="Defines the B parameter of the dispersion-damped Fumi-Tosi Potential", &
+                          description="Defines the B parameter of the dispersion-damped Fumi-Tosi potential", &
                           usage="B {real}", type_of_var=real_t, &
                           n_var=1, unit_str="angstrom^-1")
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
       CALL keyword_create(keyword, __LOCATION__, name="C", &
-                          description="Defines the C parameter of the dispersion-damped Fumi-Tosi Potential", &
+                          description="Defines the C parameter of the dispersion-damped Fumi-Tosi potential", &
                           usage="C {real}", type_of_var=real_t, &
                           n_var=1, unit_str="hartree*angstrom^6")
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
       CALL keyword_create(keyword, __LOCATION__, name="D", &
-                          description="Defines the D parameter of the dispersion-damped Fumi-Tosi Potential", &
+                          description="Defines the D parameter of the dispersion-damped Fumi-Tosi potential", &
                           usage="D {real}", type_of_var=real_t, &
                           n_var=1, unit_str="hartree*angstrom^8")
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
       CALL keyword_create(keyword, __LOCATION__, name="BD", &
-                          description="Defines the BD parameter of the dispersion-damped Fumi-Tosi Potential", &
-                          usage="D {real}", type_of_var=real_t, &
-                          n_var=1, unit_str="angstrom^-1")
+                          description="Defines the BD parameters of the dispersion-damped Fumi-Tosi potential. "// &
+                          "One or two parameter values are expected. If only one value is provided, then this "// &
+                          "value will be used both for the 6th and the 8th order term.", &
+                          usage="BD {real} {real}", type_of_var=real_t, &
+                          n_var=-1, unit_str="angstrom^-1")
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 

diff --git a/src/pair_potential.F b/src/pair_potential.F
@@ -179,7 +179,7 @@ SUBROUTINE spline_nonbond_control(spline_env, potparm, atomic_kind_set, &
                      ip = 0
                   END IF
                END IF
-            ENDIF
+            END IF
             ! Setup of Exclusion Types
             pot%no_pp = .TRUE.
             pot%no_mb = .TRUE.
@@ -314,7 +314,7 @@ SUBROUTINE get_spline_cutoff(hicut, locut, found_locut, pot, do_zbl, &
             x = x - dx2
          END DO
          x = x + dx2
-      ENDDO
+      END DO
       locut = locut_found
 
    END SUBROUTINE get_spline_cutoff
@@ -380,7 +380,7 @@ SUBROUTINE generate_spline_low(spl_p, npoints, locut, hicut, eps_spline, &
          fixed_spline_points = .FALSE.
          npoints = 20
          IF (.NOT. found_locut) npoints = 2
-      ENDIF
+      END IF
       spline_data => spl_p(1)%spline_data
       DO WHILE (.TRUE.)
          CALL init_splinexy(spline_data, npoints + 1)
@@ -609,7 +609,7 @@ SUBROUTINE get_nonbond_storage(spline_env, potparm, atomic_kind_set, do_zbl, &
          CASE (ft_type)
             nvar = 4 + nvar
          CASE (ftd_type)
-            nvar = 5 + nvar
+            nvar = 6 + nvar
          CASE (ip_type)
             nvar = 3 + nvar
          CASE (b4_type)
@@ -694,7 +694,8 @@ SUBROUTINE get_nonbond_storage(spline_env, potparm, atomic_kind_set, do_zbl, &
                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%ftd%B
                         pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%ftd%C
                         pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%ftd%D
-                        pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%ftd%BD
+                        pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%ftd%BD(1)
+                        pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%ftd%BD(2)
                      CASE (ip_type)
                         pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%ipbv%rcore
                         pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%ipbv%m
@@ -1034,7 +1035,7 @@ SUBROUTINE set_potparm_index(potparm, my_index, pot_target, ntype, tmpij_out, &
          pot%spl_f%rcutsq_f = 1.0_dp
          pot%spl_f%rscale = 1.0_dp
          pot%spl_f%fscale = 1.0_dp
-      ENDDO
+      END DO
 
       IF (ANY(potential_single_allocation == pot_target)) THEN
          DO i = 1, nvalues
@@ -1062,7 +1063,7 @@ SUBROUTINE set_potparm_index(potparm, my_index, pot_target, ntype, tmpij_out, &
                   pot%spl_f%fscale(1) = l_epsilon/m_epsilon
                END IF
             END IF
-         ENDDO
+         END DO
       END IF
 
       DO i = 1, nvalues

diff --git a/src/pair_potential_types.F b/src/pair_potential_types.F
@@ -136,7 +136,7 @@ MODULE pair_potential_types
       REAL(KIND=dp) :: B
       REAL(KIND=dp) :: C
       REAL(KIND=dp) :: D
-      REAL(KIND=dp) :: BD
+      REAL(KIND=dp), DIMENSION(2) :: BD
    END TYPE ftd_pot_type
 
 ! **************************************************************************************************
@@ -432,7 +432,7 @@ SUBROUTINE compare_pot(pot1, pot2, compare)
                 (pot1%set(i)%ftd%B == pot2%set(i)%ftd%B) .AND. &
                 (pot1%set(i)%ftd%C == pot2%set(i)%ftd%C) .AND. &
                 (pot1%set(i)%ftd%D == pot2%set(i)%ftd%D) .AND. &
-                (pot1%set(i)%ftd%BD == pot2%set(i)%ftd%BD)) mycompare = .TRUE.
+                (ALL(pot1%set(i)%ftd%BD(:) == pot2%set(i)%ftd%BD(:)))) mycompare = .TRUE.
          CASE (ip_type)
             IF ((SUM(ABS(pot1%set(i)%ipbv%a - pot2%set(i)%ipbv%a)) == 0.0_dp) .AND. &
                 (pot1%set(i)%ipbv%rcore == pot2%set(i)%ipbv%rcore) .AND. &

diff --git a/src/pair_potential_util.F b/src/pair_potential_util.F
@@ -16,6 +16,7 @@ MODULE pair_potential_util
    USE fparser,                         ONLY: EvalErrType,&
                                               evalf
    USE kinds,                           ONLY: dp
+   USE mathconstants,                   ONLY: ifac
    USE pair_potential_types,            ONLY: &
         b4_type, bm_type, ea_type, ft_type, ftd_type, gp_type, gw_type, ip_type, lj_charmm_type, &
         lj_type, not_initialized, pair_potential_single_type, wl_type
@@ -48,8 +49,8 @@ FUNCTION ener_pot(pot, r, energy_cutoff) RESULT(value)
       REAL(KIND=dp)                                      :: value
 
       INTEGER                                            :: i, index, j, n
-      REAL(KIND=dp)                                      :: damp, dampexp, dampsum, f6, f8, &
-                                                            factorial, lvalue, pp, qq, scale, xf
+      REAL(KIND=dp)                                      :: bd6, bd8, dampsum, f6, f8, lvalue, pp, &
+                                                            qq, scale, xf
 
       value = 0.0_dp
       DO j = 1, SIZE(pot%type)
@@ -74,7 +75,7 @@ FUNCTION ener_pot(pot, r, energy_cutoff) RESULT(value)
             ELSE
                ! use a linear potential
                lvalue = pot%set(j)%ipbv%m*r + pot%set(j)%ipbv%b
-            ENDIF
+            END IF
             lvalue = lvalue
          ELSE IF (pot%type(j) == wl_type) THEN
             lvalue = pot%set(j)%willis%a*EXP(-pot%set(j)%willis%b*r) - pot%set(j)%willis%c/r**6
@@ -85,31 +86,32 @@ FUNCTION ener_pot(pot, r, energy_cutoff) RESULT(value)
          ELSE IF (pot%type(j) == ft_type) THEN
             lvalue = pot%set(j)%ft%a*EXP(-pot%set(j)%ft%b*r) - pot%set(j)%ft%c/r**6 - pot%set(j)%ft%d/r**8
          ELSE IF (pot%type(j) == ftd_type) THEN
+            ! Compute 6th order dispersion correction term
+            bd6 = pot%set(j)%ftd%bd(1)
             dampsum = 1.0_dp
             xf = 1.0_dp
-            factorial = 1.0_dp
-            damp = pot%set(j)%ftd%bd
-            dampexp = dexp(-damp*r)
             DO i = 1, 6
-               xf = xf*damp*r
-               factorial = factorial*REAL(i, KIND=dp)
-               dampsum = dampsum + xf/factorial
-            ENDDO
-            f6 = 1.0_dp - dampexp*dampsum
-            DO i = 7, 8
-               xf = xf*damp*r
-               factorial = factorial*REAL(i, KIND=dp)
-               dampsum = dampsum + xf/factorial
-            ENDDO
-            f8 = 1.0_dp - dampexp*dampsum
+               xf = xf*bd6*r
+               dampsum = dampsum + xf*ifac(i)
+            END DO
+            f6 = 1.0_dp - EXP(-bd6*r)*dampsum
+            ! Compute 8th order dispersion correction term
+            bd8 = pot%set(j)%ftd%bd(2)
+            dampsum = 1.0_dp
+            xf = 1.0_dp
+            DO i = 1, 8
+               xf = xf*bd8*r
+               dampsum = dampsum + xf*ifac(i)
+            END DO
+            f8 = 1.0_dp - EXP(-bd8*r)*dampsum
             lvalue = pot%set(j)%ftd%a*EXP(-pot%set(j)%ftd%b*r) - f6*pot%set(j)%ftd%c/r**6 - f8*pot%set(j)%ftd%d/r**8
          ELSE IF (pot%type(j) == ea_type) THEN
             index = INT(r/pot%set(j)%eam%drar) + 1
             IF (index > pot%set(j)%eam%npoints) THEN
                index = pot%set(j)%eam%npoints
             ELSEIF (index < 1) THEN
                index = 1
-            ENDIF
+            END IF
             qq = r - pot%set(j)%eam%rval(index)
             pp = pot%set(j)%eam%phi(index) + &
                  qq*pot%set(j)%eam%phip(index)