Skip to content

Commit

Permalink
Fist: Add GAL21 forcefield
Browse files Browse the repository at this point in the history
  • Loading branch information
@pclabaut
pclabaut committed Jul 9, 2021
1 parent 7a55edd commit c923b2c
Show file tree
Hide file tree
Showing 15 changed files with 3,615 additions and 63 deletions.
14 changes: 13 additions & 1 deletion src/common/bibliography.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ MODULE bibliography
Stoychev2016, Futera2017, Bailey2006, Papior2017, Lehtola2018, &
Brieuc2016, Barca2018, Scheiber2018, Huang2011, Heaton_Burgess2007, &
Schuett2018, Holmberg2018, Togo2018, Staub2019, Grimme2013, Grimme2016, &
Grimme2017, Kondov2007, Clabaut2020, &
Grimme2017, Kondov2007, Clabaut2020, Clabaut2021, &
Ren2011, Ren2013, Cohen2000, Rogers2002, Filippetti2000, &
Limpanuparb2011, Martin2003, Yin2017, Goerigk2017, &
Wilhelm2016a, Wilhelm2016b, Wilhelm2017, Wilhelm2018, Lass2018, cp2kqs2020, &
Expand Down Expand Up @@ -4118,6 +4118,18 @@ SUBROUTINE add_all_references()
"ER"), &
DOI="10.1021/acs.jctc.0c00091")

CALL add_reference(key=Clabaut2021, ISI_record=s2a( &
"AU Clabaut, P", &
"AF Clabaut, P", &
"TI Ph.D. Solvation and adsorptions at the solid/water interface : Developments and applications", &
"SO ", &
"SN None", &
"PY 2021", &
"VL None", &
"DI None", &
"ER"), &
DOI="None")

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 @@ -40,7 +40,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: gal_type,&
USE pair_potential_types, ONLY: gal21_type,&
gal_type,&
pair_potential_pp_type,&
siepmann_type,&
tersoff_type
Expand Down Expand Up @@ -224,6 +225,9 @@ SUBROUTINE list_control(atomic_kind_set, particle_set, local_particles, &
IF (ANY(potparm%pot(ikind, jkind)%pot%type == gal_type)) THEN
full_nl(ikind, jkind) = .TRUE.
END IF
IF (ANY(potparm%pot(ikind, jkind)%pot%type == gal21_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 @@ -18,7 +18,8 @@ MODULE fist_nonbond_env_types
fist_neighbor_type
USE kinds, ONLY: default_string_length,&
dp
USE pair_potential_types, ONLY: gal_type,&
USE pair_potential_types, ONLY: gal21_type,&
gal_type,&
pair_potential_pp_release,&
pair_potential_pp_type,&
siepmann_type,&
Expand Down Expand Up @@ -409,6 +410,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 == gal21_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
15 changes: 5 additions & 10 deletions src/fist_nonbond_force.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,14 +38,9 @@ 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: gal_type,&
nosh_nosh,&
nosh_sh,&
pair_potential_pp_type,&
pair_potential_single_type,&
sh_sh,&
siepmann_type,&
tersoff_type
USE pair_potential_types, ONLY: &
gal21_type, gal_type, nosh_nosh, nosh_sh, pair_potential_pp_type, &
pair_potential_single_type, sh_sh, siepmann_type, tersoff_type
USE particle_types, ONLY: particle_type
USE shell_potential_types, ONLY: get_shell,&
shell_kind_type
Expand Down Expand Up @@ -234,7 +229,7 @@ SUBROUTINE force_nonbond(fist_nonbond_env, ewald_env, particle_set, cell, &
fac_ei = fac_kind
fac_vdw = fac_kind
full_nl = ANY(pot%type == tersoff_type) .OR. ANY(pot%type == siepmann_type) &
.OR. ANY(pot%type == gal_type)
.OR. ANY(pot%type == gal_type) .OR. ANY(pot%type == gal21_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 @@ -728,7 +723,7 @@ 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) &
.OR. ANY(pot%type == gal_type)
.OR. ANY(pot%type == gal_type) .OR. ANY(pot%type == gal21_type)
IF ((.NOT. full_nl) .AND. (atom_a == atom_b)) THEN
fac_ei = fac_ei*0.5_dp
END IF
Expand Down
145 changes: 138 additions & 7 deletions src/force_fields_input.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
! **************************************************************************************************
MODULE force_fields_input
USE bibliography, ONLY: Clabaut2020,&
Clabaut2021,&
Siepmann1995,&
Tersoff1988,&
Tosi1964a,&
Expand Down Expand Up @@ -63,9 +64,10 @@ MODULE force_fields_input
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, 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
ft_type, ftd_type, gal21_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,&
shell_p_type
USE string_utilities, ONLY: uppercase
Expand Down Expand Up @@ -98,8 +100,8 @@ SUBROUTINE read_force_field_section1(ff_section, mm_section, ff_type, para_env)
TYPE(force_field_type), INTENT(INOUT) :: ff_type
TYPE(cp_para_env_type), POINTER :: para_env

INTEGER :: nb4, nbends, nbm, nbmhft, nbmhftd, nbonds, nchg, neam, ngal, ngd, ngp, nimpr, &
nipbv, nlj, nopbend, nquip, nshell, nsiepmann, ntersoff, ntors, ntot, nubs, nwl
INTEGER :: nb4, nbends, nbm, nbmhft, nbmhftd, nbonds, nchg, neam, ngal, ngal21, 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 @@ -252,17 +254,25 @@ SUBROUTINE read_force_field_section1(ff_section, mm_section, ff_type, para_env)
CALL read_gal_section(inp_info%nonbonded, tmp_section2, ntot, tmp_section2)
END IF
tmp_section2 => section_vals_get_subs_vals(tmp_section, "GAL21")
CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=ngal21)
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 + ngal21, gal21=.TRUE.)
CALL read_gal21_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 + ngal
ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff + ngal + ngal21
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 + ngal + nsiepmann
ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff + ngal + ngal21 + 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 @@ -1340,6 +1350,127 @@ SUBROUTINE read_gal_section(nonbonded, section, start, gal_section)
END DO
END SUBROUTINE read_gal_section
! **************************************************************************************************
!> \brief Reads the gal21 section
!> \param nonbonded ...
!> \param section ...
!> \param start ...
!> \param gal21_section ...
!> \author Clabaut Paul
! **************************************************************************************************
SUBROUTINE read_gal21_section(nonbonded, section, start, gal21_section)
TYPE(pair_potential_p_type), POINTER :: nonbonded
TYPE(section_vals_type), POINTER :: section
INTEGER, INTENT(IN) :: start
TYPE(section_vals_type), POINTER :: gal21_section
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(Clabaut2021)
CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
nonbonded%pot(start + isec)%pot%type = gal21_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)%gal21%met1 = atm_names(1)
nonbonded%pot(start + isec)%pot%set(1)%gal21%met2 = atm_names(2)
CALL section_vals_val_get(gal21_section, "epsilon", i_rep_section=isec, r_vals=rvalues)
nonbonded%pot(start + isec)%pot%set(1)%gal21%epsilon1 = rvalues(1)
nonbonded%pot(start + isec)%pot%set(1)%gal21%epsilon2 = rvalues(2)
nonbonded%pot(start + isec)%pot%set(1)%gal21%epsilon3 = rvalues(3)
CALL section_vals_val_get(gal21_section, "bxy", i_rep_section=isec, r_vals=rvalues)
nonbonded%pot(start + isec)%pot%set(1)%gal21%bxy1 = rvalues(1)
nonbonded%pot(start + isec)%pot%set(1)%gal21%bxy2 = rvalues(2)
CALL section_vals_val_get(gal21_section, "bz", i_rep_section=isec, r_vals=rvalues)
nonbonded%pot(start + isec)%pot%set(1)%gal21%bz1 = rvalues(1)
nonbonded%pot(start + isec)%pot%set(1)%gal21%bz2 = rvalues(2)
CALL section_vals_val_get(gal21_section, "r", i_rep_section=isec, r_vals=rvalues)
nonbonded%pot(start + isec)%pot%set(1)%gal21%r1 = rvalues(1)
nonbonded%pot(start + isec)%pot%set(1)%gal21%r2 = rvalues(2)
CALL section_vals_val_get(gal21_section, "a1", i_rep_section=isec, r_vals=rvalues)
nonbonded%pot(start + isec)%pot%set(1)%gal21%a11 = rvalues(1)
nonbonded%pot(start + isec)%pot%set(1)%gal21%a12 = rvalues(2)
nonbonded%pot(start + isec)%pot%set(1)%gal21%a13 = rvalues(3)
CALL section_vals_val_get(gal21_section, "a2", i_rep_section=isec, r_vals=rvalues)
nonbonded%pot(start + isec)%pot%set(1)%gal21%a21 = rvalues(1)
nonbonded%pot(start + isec)%pot%set(1)%gal21%a22 = rvalues(2)
nonbonded%pot(start + isec)%pot%set(1)%gal21%a23 = rvalues(3)
CALL section_vals_val_get(gal21_section, "a3", i_rep_section=isec, r_vals=rvalues)
nonbonded%pot(start + isec)%pot%set(1)%gal21%a31 = rvalues(1)
nonbonded%pot(start + isec)%pot%set(1)%gal21%a32 = rvalues(2)
nonbonded%pot(start + isec)%pot%set(1)%gal21%a33 = rvalues(3)
CALL section_vals_val_get(gal21_section, "a4", i_rep_section=isec, r_vals=rvalues)
nonbonded%pot(start + isec)%pot%set(1)%gal21%a41 = rvalues(1)
nonbonded%pot(start + isec)%pot%set(1)%gal21%a42 = rvalues(2)
nonbonded%pot(start + isec)%pot%set(1)%gal21%a43 = rvalues(3)
CALL section_vals_val_get(gal21_section, "A", i_rep_section=isec, r_vals=rvalues)
nonbonded%pot(start + isec)%pot%set(1)%gal21%AO1 = rvalues(1)
nonbonded%pot(start + isec)%pot%set(1)%gal21%AO2 = rvalues(2)
CALL section_vals_val_get(gal21_section, "B", i_rep_section=isec, r_vals=rvalues)
nonbonded%pot(start + isec)%pot%set(1)%gal21%BO1 = rvalues(1)
nonbonded%pot(start + isec)%pot%set(1)%gal21%BO2 = rvalues(2)
CALL section_vals_val_get(gal21_section, "C", i_rep_section=isec, &
r_val=nonbonded%pot(start + isec)%pot%set(1)%gal21%c)
CALL section_vals_val_get(gal21_section, "AH", i_rep_section=isec, r_vals=rvalues)
nonbonded%pot(start + isec)%pot%set(1)%gal21%AH1 = rvalues(1)
nonbonded%pot(start + isec)%pot%set(1)%gal21%AH2 = rvalues(2)
CALL section_vals_val_get(gal21_section, "BH", i_rep_section=isec, r_vals=rvalues)
nonbonded%pot(start + isec)%pot%set(1)%gal21%BH1 = rvalues(1)
nonbonded%pot(start + isec)%pot%set(1)%gal21%BH2 = rvalues(2)
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)%gal21%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)%gal21%gcn(iatom) = rval
END DO
CALL section_vals_val_get(gal21_section, "Fit_express", i_rep_section=isec, &
l_val=nonbonded%pot(start + isec)%pot%set(1)%gal21%express)
! ! In case it is defined override the standard specification of RCUT
CALL section_vals_val_get(gal21_section, "RCUT", i_rep_section=isec, n_rep_val=n_rep)
IF (n_rep == 1) THEN
CALL section_vals_val_get(gal21_section, "RCUT", i_rep_section=isec, r_val=rcut)
nonbonded%pot(start + isec)%pot%rcutsq = rcut**2
nonbonded%pot(start + isec)%pot%set(1)%gal21%rcutsq = rcut**2
END IF
END DO
END SUBROUTINE read_gal21_section
! **************************************************************************************************
!> \brief Reads the siepmann section
!> \param nonbonded ...
Expand Down

0 comments on commit c923b2c

Please sign in to comment.