Skip to content

Commit

Permalink
Add an implementation of the GAL19 forcefield and suitable regtests
Browse files Browse the repository at this point in the history
  • Loading branch information
@pclabaut @oschuett
pclabaut authored and oschuett committed Jun 23, 2020
1 parent b021919 commit b3d5a84
Show file tree
Hide file tree
Showing 15 changed files with 3,417 additions and 46 deletions.
20 changes: 19 additions & 1 deletion src/common/bibliography.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ MODULE bibliography
Stoychev2016, Futera2017, Bailey2006, Papior2017, Lehtola2018, &
Brieuc2016, Barca2018, Scheiber2018, Huang2011, Heaton_Burgess2007, &
Schuett2018, Holmberg2018, Togo2018, Staub2019, Grimme2013, Grimme2016, &
Grimme2017, Kondov2007, &
Grimme2017, Kondov2007, Clabaut2020, &
Ren2011, Ren2013, Cohen2000, Rogers2002, Filippetti2000, Paziani2006, &
Toulouse2004, Limpanuparb2011, Martin2003, Yin2017, Goerigk2017, &
Wilhelm2016a, Wilhelm2016b, Wilhelm2017, Wilhelm2018, Lass2018, cp2kqs2020
Expand Down Expand Up @@ -4047,6 +4047,24 @@ SUBROUTINE add_all_references()
"ER"), &
DOI="10.1021/acs.jctc.8b00957")

CALL add_reference(key=Clabaut2020, ISI_record=s2a( &
"AU Clabaut, P", &
" Fleurat-Lessard, P", &
" Michel, C", &
" Steinmann, SN", &
"AF Clabaut, P", &
" Fleurat-Lessard, P", &
" Michel, C", &
" Steinmann, SN", &
"TI Ten Facets, One Force Field: The GAL19 Force Field for Water-Noble Metal Interfaces", &
"SO Journal of Chemical Theory and Computation", &
"SN None", &
"PY 2020", &
"VL None", &
"DI 10.1021/acs.jctc.0c00091", &
"ER"), &
DOI="10.1021/acs.jctc.0c00091")

CALL add_reference(key=Richters2018, ISI_record=s2a( &
"AU Richters, D", &
" Lass, M", &
Expand Down
6 changes: 5 additions & 1 deletion src/fist_neighbor_list_control.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,8 @@ MODULE fist_neighbor_list_control
section_vals_val_get
USE kinds, ONLY: dp
USE message_passing, ONLY: mp_max
USE pair_potential_types, ONLY: pair_potential_pp_type,&
USE pair_potential_types, ONLY: gal_type,&
pair_potential_pp_type,&
siepmann_type,&
tersoff_type
USE particle_types, ONLY: particle_type
Expand Down Expand Up @@ -215,6 +216,9 @@ SUBROUTINE list_control(atomic_kind_set, particle_set, local_particles, &
IF (ANY(potparm%pot(ikind, jkind)%pot%type == siepmann_type)) THEN
full_nl(ikind, jkind) = .TRUE.
END IF
IF (ANY(potparm%pot(ikind, jkind)%pot%type == gal_type)) THEN
full_nl(ikind, jkind) = .TRUE.
END IF
full_nl(jkind, ikind) = full_nl(ikind, jkind)
END DO
END DO
Expand Down
7 changes: 6 additions & 1 deletion src/fist_nonbond_env_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,8 @@ MODULE fist_nonbond_env_types
fist_neighbor_type
USE kinds, ONLY: default_string_length,&
dp
USE pair_potential_types, ONLY: pair_potential_pp_release,&
USE pair_potential_types, ONLY: gal_type,&
pair_potential_pp_release,&
pair_potential_pp_type,&
siepmann_type,&
tersoff_type
Expand Down Expand Up @@ -408,6 +409,10 @@ SUBROUTINE init_fist_nonbond_env(fist_nonbond_env, atomic_kind_set, &
fist_nonbond_env%ij_kind_full_fac(idim, jdim) = 0.5_dp
fist_nonbond_env%ij_kind_full_fac(idim, jdim) = 0.5_dp
END IF
IF (ANY(potparm%pot(idim, jdim)%pot%type == gal_type)) THEN
fist_nonbond_env%ij_kind_full_fac(idim, jdim) = 0.5_dp
fist_nonbond_env%ij_kind_full_fac(idim, jdim) = 0.5_dp
END IF
ELSE
! In case we don't use potparm for initialization let's account
! only for the real-space part of the Ewald sum.
Expand Down
9 changes: 6 additions & 3 deletions src/fist_nonbond_force.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,8 @@ MODULE fist_nonbond_force
USE mathlib, ONLY: matvec_3x3
USE message_passing, ONLY: mp_sum
USE pair_potential_coulomb, ONLY: potential_coulomb
USE pair_potential_types, ONLY: nosh_nosh,&
USE pair_potential_types, ONLY: gal_type,&
nosh_nosh,&
nosh_sh,&
pair_potential_pp_type,&
pair_potential_single_type,&
Expand Down Expand Up @@ -230,7 +231,8 @@ SUBROUTINE force_nonbond(fist_nonbond_env, ewald_env, particle_set, cell, &
! Determine the scaling factors
fac_ei = fac_kind
fac_vdw = fac_kind
full_nl = ANY(pot%type == tersoff_type) .OR. ANY(pot%type == siepmann_type)
full_nl = ANY(pot%type == tersoff_type) .OR. ANY(pot%type == siepmann_type) &
.OR. ANY(pot%type == gal_type)
IF ((.NOT. full_nl) .AND. (atom_a == atom_b)) THEN
fac_ei = 0.5_dp*fac_ei
fac_vdw = 0.5_dp*fac_vdw
Expand Down Expand Up @@ -724,7 +726,8 @@ SUBROUTINE bonded_correct_gaussian(fist_nonbond_env, atomic_kind_set, &

! Determine the scaling factors
fac_ei = ij_kind_full_fac(kind_a, kind_b)
full_nl = ANY(pot%type == tersoff_type) .OR. ANY(pot%type == siepmann_type)
full_nl = ANY(pot%type == tersoff_type) .OR. ANY(pot%type == siepmann_type) &
.OR. ANY(pot%type == gal_type)
IF ((.NOT. full_nl) .AND. (atom_a == atom_b)) THEN
fac_ei = fac_ei*0.5_dp
END IF
Expand Down
122 changes: 116 additions & 6 deletions src/force_fields_input.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,12 +16,15 @@
!> \author CJM
! **************************************************************************************************
MODULE force_fields_input
USE bibliography, ONLY: Siepmann1995,&
USE bibliography, ONLY: Clabaut2020,&
Siepmann1995,&
Tersoff1988,&
Tosi1964a,&
Tosi1964b,&
Yamada2000,&
cite_reference
USE cp_linked_list_input, ONLY: cp_sll_val_next,&
cp_sll_val_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type,&
cp_to_string
Expand All @@ -46,16 +49,19 @@ MODULE force_fields_input
USE force_fields_util, ONLY: get_generic_info
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_list_get,&
section_vals_type,&
section_vals_val_get
USE input_val_types, ONLY: val_get,&
val_type
USE kinds, ONLY: default_string_length,&
dp
USE mathconstants, ONLY: pi
USE mathlib, ONLY: invert_matrix
USE memory_utilities, ONLY: reallocate
USE pair_potential_types, ONLY: &
b4_type, bm_type, do_potential_single_allocation, ea_type, eam_pot_type, ft_pot_type, &
ft_type, ftd_type, gp_type, gw_type, ip_type, ipbv_pot_type, lj_charmm_type, &
ft_type, ftd_type, gal_type, gp_type, gw_type, ip_type, ipbv_pot_type, lj_charmm_type, &
no_potential_single_allocation, pair_potential_p_type, pair_potential_reallocate, &
potential_single_allocation, quip_type, siepmann_type, tersoff_type, wl_type
USE shell_potential_types, ONLY: shell_p_create,&
Expand Down Expand Up @@ -93,8 +99,8 @@ SUBROUTINE read_force_field_section1(ff_section, mm_section, ff_type, para_env)
CHARACTER(len=*), PARAMETER :: routineN = 'read_force_field_section1', &
routineP = moduleN//':'//routineN

INTEGER :: nb4, nbends, nbm, nbmhft, nbmhftd, nbonds, nchg, neam, ngd, ngp, nimpr, nipbv, &
nlj, nopbend, nquip, nshell, nsiepmann, ntersoff, ntors, ntot, nubs, nwl
INTEGER :: nb4, nbends, nbm, nbmhft, nbmhftd, nbonds, nchg, neam, ngal, ngd, ngp, nimpr, &
nipbv, nlj, nopbend, nquip, nshell, nsiepmann, ntersoff, ntors, ntot, nubs, nwl
LOGICAL :: explicit, unique_spline
REAL(KIND=dp) :: min_eps_spline_allowed
TYPE(input_info_type), POINTER :: inp_info
Expand Down Expand Up @@ -238,17 +244,26 @@ SUBROUTINE read_force_field_section1(ff_section, mm_section, ff_type, para_env)
CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + ntersoff, tersoff=.TRUE.)
CALL read_tersoff_section(inp_info%nonbonded, tmp_section2, ntot, tmp_section2)
END IF
tmp_section2 => section_vals_get_subs_vals(tmp_section, "GAL19")
CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=ngal)
ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff
IF (explicit) THEN
CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + ngal, gal=.TRUE.)
CALL read_gal_section(inp_info%nonbonded, tmp_section2, ntot, tmp_section2)
END IF
tmp_section2 => section_vals_get_subs_vals(tmp_section, "SIEPMANN")
CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nsiepmann)
ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff
ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff + ngal
IF (explicit) THEN
CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nsiepmann, siepmann=.TRUE.)
CALL read_siepmann_section(inp_info%nonbonded, tmp_section2, ntot, tmp_section2)
END IF
tmp_section2 => section_vals_get_subs_vals(tmp_section, "quip")
CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nquip)
ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff + nsiepmann
ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff + ngal + nsiepmann
IF (explicit) THEN
CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nquip, quip=.TRUE.)
CALL read_quip_section(inp_info%nonbonded, tmp_section2, ntot)
Expand Down Expand Up @@ -1275,6 +1290,101 @@ SUBROUTINE read_tersoff_section(nonbonded, section, start, tersoff_section)
END DO
END SUBROUTINE read_tersoff_section
! **************************************************************************************************
!> \brief Reads the gal19 section
!> \param nonbonded ...
!> \param section ...
!> \param start ...
!> \param gal_section ...
!> \author Clabaut Paul
! **************************************************************************************************
SUBROUTINE read_gal_section(nonbonded, section, start, gal_section)
TYPE(pair_potential_p_type), POINTER :: nonbonded
TYPE(section_vals_type), POINTER :: section
INTEGER, INTENT(IN) :: start
TYPE(section_vals_type), POINTER :: gal_section
CHARACTER(len=*), PARAMETER :: routineN = 'read_gal_section', &
routineP = moduleN//':'//routineN
CHARACTER(LEN=default_string_length), &
DIMENSION(:), POINTER :: atm_names
INTEGER :: iatom, isec, n_items, n_rep, nval
LOGICAL :: is_ok
REAL(KIND=dp) :: rcut, rval
REAL(KIND=dp), DIMENSION(:), POINTER :: rvalues
TYPE(cp_sll_val_type), POINTER :: list
TYPE(section_vals_type), POINTER :: subsection
TYPE(val_type), POINTER :: val
CALL section_vals_get(section, n_repetition=n_items)
DO isec = 1, n_items
CALL cite_reference(Clabaut2020)
CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
nonbonded%pot(start + isec)%pot%type = gal_type
nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
CALL section_vals_val_get(section, "METALS", i_rep_section=isec, c_vals=atm_names)
nonbonded%pot(start + isec)%pot%set(1)%gal%met1 = atm_names(1)
nonbonded%pot(start + isec)%pot%set(1)%gal%met2 = atm_names(2)
CALL section_vals_val_get(gal_section, "epsilon", i_rep_section=isec, &
r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%epsilon)
CALL section_vals_val_get(gal_section, "bxy", i_rep_section=isec, &
r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%bxy)
CALL section_vals_val_get(gal_section, "bz", i_rep_section=isec, &
r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%bz)
CALL section_vals_val_get(gal_section, "r", i_rep_section=isec, r_vals=rvalues)
nonbonded%pot(start + isec)%pot%set(1)%gal%r1 = rvalues(1)
nonbonded%pot(start + isec)%pot%set(1)%gal%r2 = rvalues(2)
CALL section_vals_val_get(gal_section, "a1", i_rep_section=isec, &
r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%a1)
CALL section_vals_val_get(gal_section, "a2", i_rep_section=isec, &
r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%a2)
CALL section_vals_val_get(gal_section, "a3", i_rep_section=isec, &
r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%a3)
CALL section_vals_val_get(gal_section, "a4", i_rep_section=isec, &
r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%a4)
CALL section_vals_val_get(gal_section, "A", i_rep_section=isec, &
r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%a)
CALL section_vals_val_get(gal_section, "B", i_rep_section=isec, &
r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%b)
CALL section_vals_val_get(gal_section, "C", i_rep_section=isec, &
r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%c)
NULLIFY (list)
subsection => section_vals_get_subs_vals(section, "GCN", i_rep_section=isec)
CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", n_rep_val=nval)
ALLOCATE (nonbonded%pot(start + isec)%pot%set(1)%gal%gcn(nval))
CALL section_vals_list_get(subsection, "_DEFAULT_KEYWORD_", list=list)
DO iatom = 1, nval
! we use only the first default_string_length characters of each line
is_ok = cp_sll_val_next(list, val)
CALL val_get(val, r_val=rval)
! assign values
nonbonded%pot(start + isec)%pot%set(1)%gal%gcn(iatom) = rval
END DO
!WRITE(*,*) "Youhou", SIZE(nonbonded%pot(start+isec)%pot%set(1)%gal%gcn)
CALL section_vals_val_get(gal_section, "Fit_express", i_rep_section=isec, &
l_val=nonbonded%pot(start + isec)%pot%set(1)%gal%express)
! ! In case it is defined override the standard specification of RCUT
CALL section_vals_val_get(gal_section, "RCUT", i_rep_section=isec, n_rep_val=n_rep)
IF (n_rep == 1) THEN
CALL section_vals_val_get(gal_section, "RCUT", i_rep_section=isec, r_val=rcut)
nonbonded%pot(start + isec)%pot%rcutsq = rcut**2
nonbonded%pot(start + isec)%pot%set(1)%gal%rcutsq = rcut**2
END IF
END DO
END SUBROUTINE read_gal_section
! **************************************************************************************************
!> \brief Reads the siepmann section
!> \param nonbonded ...
Expand Down

0 comments on commit b3d5a84

Please sign in to comment.