From c923b2c4bf6affc9555ff972ba7c13d028cc6381 Mon Sep 17 00:00:00 2001 From: pclabaut <60351749+pclabaut@users.noreply.github.com> Date: Fri, 9 Jul 2021 14:51:05 +0200 Subject: [PATCH] Fist: Add GAL21 forcefield --- src/common/bibliography.F | 14 +- src/fist_neighbor_list_control.F | 6 +- src/fist_nonbond_env_types.F | 7 +- src/fist_nonbond_force.F | 15 +- src/force_fields_input.F | 145 +++- src/input_cp2k_mm.F | 157 +++- src/manybody_gal21.F | 937 ++++++++++++++++++++++++ src/manybody_potential.F | 186 ++++- src/pair_potential.F | 43 +- src/pair_potential_types.F | 272 ++++++- tests/Fist/regtest-gal21/TEST_FILES | 7 + tests/Fist/regtest-gal21/coord.inc | 870 ++++++++++++++++++++++ tests/Fist/regtest-gal21/gcn.inc | 870 ++++++++++++++++++++++ tests/Fist/regtest-gal21/test_gal21.inp | 148 ++++ tests/TEST_DIRS | 1 + 15 files changed, 3615 insertions(+), 63 deletions(-) create mode 100644 src/manybody_gal21.F create mode 100644 tests/Fist/regtest-gal21/TEST_FILES create mode 100644 tests/Fist/regtest-gal21/coord.inc create mode 100644 tests/Fist/regtest-gal21/gcn.inc create mode 100644 tests/Fist/regtest-gal21/test_gal21.inp diff --git a/src/common/bibliography.F b/src/common/bibliography.F index d9375e52e0..c336b8e0bc 100644 --- a/src/common/bibliography.F +++ b/src/common/bibliography.F @@ -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, & @@ -4118,6 +4118,18 @@ CONTAINS "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", & diff --git a/src/fist_neighbor_list_control.F b/src/fist_neighbor_list_control.F index 36351a1f04..078e669629 100644 --- a/src/fist_neighbor_list_control.F +++ b/src/fist_neighbor_list_control.F @@ -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 @@ -224,6 +225,9 @@ CONTAINS 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 diff --git a/src/fist_nonbond_env_types.F b/src/fist_nonbond_env_types.F index 0ff7c0c870..a1c4990121 100644 --- a/src/fist_nonbond_env_types.F +++ b/src/fist_nonbond_env_types.F @@ -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,& @@ -409,6 +410,10 @@ CONTAINS 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. diff --git a/src/fist_nonbond_force.F b/src/fist_nonbond_force.F index 1297cd6701..842039bf09 100644 --- a/src/fist_nonbond_force.F +++ b/src/fist_nonbond_force.F @@ -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 @@ -234,7 +229,7 @@ CONTAINS 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 @@ -728,7 +723,7 @@ CONTAINS ! 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 diff --git a/src/force_fields_input.F b/src/force_fields_input.F index 30cb4e3930..f93afc8b9f 100644 --- a/src/force_fields_input.F +++ b/src/force_fields_input.F @@ -19,6 +19,7 @@ ! ************************************************************************************************** MODULE force_fields_input USE bibliography, ONLY: Clabaut2020,& + Clabaut2021,& Siepmann1995,& Tersoff1988,& Tosi1964a,& @@ -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 @@ -98,8 +100,8 @@ CONTAINS 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 @@ -252,9 +254,17 @@ CONTAINS 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) @@ -262,7 +272,7 @@ CONTAINS 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) @@ -1340,6 +1350,127 @@ CONTAINS 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 ... diff --git a/src/input_cp2k_mm.F b/src/input_cp2k_mm.F index 125dfc2128..db4fd6ba3a 100644 --- a/src/input_cp2k_mm.F +++ b/src/input_cp2k_mm.F @@ -15,8 +15,8 @@ ! ************************************************************************************************** MODULE input_cp2k_mm USE bibliography, ONLY: & - Clabaut2020, Devynck2012, Dick1958, Foiles1986, Mitchell1993, QUIP_ref, Siepmann1995, & - Tersoff1988, Tosi1964a, Tosi1964b, Yamada2000 + Clabaut2020, Clabaut2021, Devynck2012, Dick1958, Foiles1986, Mitchell1993, QUIP_ref, & + Siepmann1995, Tersoff1988, Tosi1964a, Tosi1964b, Yamada2000 USE cp_output_handling, ONLY: cp_print_key_section_create,& debug_print_level,& high_print_level,& @@ -1210,6 +1210,10 @@ CONTAINS CALL section_add_subsection(section, subsection) CALL section_release(subsection) + CALL create_Gal21_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + END SUBROUTINE create_NONBONDED_section ! ************************************************************************************************** @@ -2515,6 +2519,151 @@ CONTAINS n_var=1, unit_str="angstrom") CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="Fit_express", & + description="Demands the particular output needed to a least square fit", & + usage="Fit_express TRUE", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + CALL create_GCN_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + + END SUBROUTINE create_Gal_section + +! ************************************************************************************************** +!> \brief This section specifies the input parameters for GAL21 +!> potential type +!> (??) +!> \param section ... +! ************************************************************************************************** + SUBROUTINE create_Gal21_section(section) + TYPE(section_type), POINTER :: section + + TYPE(keyword_type), POINTER :: keyword + TYPE(section_type), POINTER :: subsection + + CPASSERT(.NOT. ASSOCIATED(section)) + CALL section_create(section, __LOCATION__, name="GAL21", & + description="Implementation of the GAL21 forcefield, see associated paper", & + citations=(/Clabaut2021/), n_keywords=1, n_subsections=1, repeats=.TRUE.) + + NULLIFY (keyword, subsection) + + CALL keyword_create(keyword, __LOCATION__, name="ATOMS", & + description="Defines the atomic kind involved in the nonbond potential", & + usage="ATOMS {KIND1} {KIND2}", type_of_var=char_t, & + n_var=2) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="METALS", & + description="Defines the two atomic kinds to be considered as part of the metallic phase in the system", & + usage="METALS {KIND1} {KIND2} ..", type_of_var=char_t, & + n_var=2) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="epsilon", & + description="Defines the epsilon parameter of GAL21 potential", & + usage="epsilon {real} {real} {real}", type_of_var=real_t, & + n_var=3, unit_str="kcalmol") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="bxy", & + description="Defines the b perpendicular parameter of GAL21 potential", & + usage="bxy {real} {real}", type_of_var=real_t, & + n_var=2, unit_str="angstrom^-2") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="bz", & + description="Defines the b parallel parameter of GAL21 potential", & + usage="bz {real} {real}", type_of_var=real_t, & + n_var=2, unit_str="angstrom^-2") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="r", & + description="Defines the R_0 parameters of GAL21 potential for the two METALS."// & + "This is the only parameter that is shared between the two section of the forcefield in the case of two metals (alloy)."// & + "If one metal only is present, a second number should be given but won't be read", & + usage="r {real} {real}", type_of_var=real_t, n_var=2, unit_str="angstrom") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="a1", & + description="Defines the a1 parameter of GAL21 potential", & + usage="a1 {real} {real} {real}", type_of_var=real_t, & + n_var=3, unit_str="kcalmol") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="a2", & + description="Defines the a2 parameter of GAL21 potential", & + usage="a2 {real} {real} {real}", type_of_var=real_t, & + n_var=3, unit_str="kcalmol") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="a3", & + description="Defines the a3 parameter of GAL21 potential", & + usage="a3 {real} {real} {real}", type_of_var=real_t, & + n_var=3, unit_str="kcalmol") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="a4", & + description="Defines the a4 parameter of GAL21 potential", & + usage="a4 {real} {real} {real}", type_of_var=real_t, & + n_var=3, unit_str="kcalmol") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="A", & + description="Defines the A parameter of GAL21 potential", & + usage="A {real} {real}", type_of_var=real_t, & + n_var=2, unit_str="kcalmol") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="B", & + description="Defines the B parameter of GAL21 potential", & + usage="B {real} {real}", type_of_var=real_t, & + n_var=2, 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 GAL21 potential", & + usage="C {real}", type_of_var=real_t, & + n_var=1, unit_str="angstrom^6*kcalmol") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="AH", & + description="Defines the AH parameter of GAL21 potential", & + usage="AH {real} {real}", type_of_var=real_t, & + n_var=2, unit_str="kcalmol") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="BH", & + description="Defines the BH parameter of GAL21 potential", & + usage="BH {real} {real}", type_of_var=real_t, & + n_var=2, unit_str="angstrom^-1") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="RCUT", & + description="Defines the cutoff parameter of GAL21 potential", & + usage="RCUT {real}", type_of_var=real_t, & + default_r_val=cp_unit_to_cp2k(value=3.2_dp, & + unit_str="angstrom"), & + n_var=1, unit_str="angstrom") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) CALL keyword_create(keyword, __LOCATION__, name="Fit_express", & description="Demands the particular output needed to a least square fit", & @@ -2527,10 +2676,10 @@ CONTAINS CALL section_add_subsection(section, subsection) CALL section_release(subsection) - END SUBROUTINE create_Gal_section + END SUBROUTINE create_Gal21_section ! ************************************************************************************************** -!> \brief This section specifies the input parameters for the subsection GCN of GAL19 +!> \brief This section specifies the input parameters for the subsection GCN of GAL19 and GAL21 !> potential type !> (??) !> \param section ... diff --git a/src/manybody_gal21.F b/src/manybody_gal21.F new file mode 100644 index 0000000000..8866cc1907 --- /dev/null +++ b/src/manybody_gal21.F @@ -0,0 +1,937 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2021 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Implementation of the GAL21 potential +!> +!> \author Clabaut Paul +! ************************************************************************************************** +MODULE manybody_gal21 + + USE atomic_kind_types, ONLY: get_atomic_kind + USE cell_types, ONLY: cell_type,& + pbc + USE cp_log_handling, ONLY: cp_get_default_logger,& + cp_logger_type,& + cp_to_string + USE cp_output_handling, ONLY: cp_print_key_finished_output,& + cp_print_key_unit_nr + USE cp_para_types, ONLY: cp_para_env_type + USE fist_neighbor_list_types, ONLY: fist_neighbor_type,& + neighbor_kind_pairs_type + USE fist_nonbond_env_types, ONLY: pos_type + USE input_section_types, ONLY: section_vals_type + USE kinds, ONLY: dp + USE mathlib, ONLY: matvec_3x3 + USE message_passing, ONLY: mp_sum + USE pair_potential_types, ONLY: gal21_pot_type,& + gal21_type,& + pair_potential_pp_type,& + pair_potential_single_type + USE particle_types, ONLY: particle_type + USE util, ONLY: sort +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + PUBLIC :: setup_gal21_arrays, destroy_gal21_arrays, & + gal21_energy, gal21_forces, & + print_nr_ions_gal21 + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'manybody_gal21' + +CONTAINS + +! ************************************************************************************************** +!> \brief Main part of the energy evaluation of GAL2119 +!> \param pot_loc value of total potential energy +!> \param gal21 all parameters of GAL2119 +!> \param r_last_update_pbc position of every atoms on previous frame +!> \param iparticle first index of the atom of the evaluated pair +!> \param jparticle second index of the atom of the evaluated pair +!> \param cell dimension of the pbc cell +!> \param particle_set full list of atoms of the system +!> \param mm_section ... +!> \author Clabaut Paul - 2019 - ENS de Lyon +! ************************************************************************************************** + SUBROUTINE gal21_energy(pot_loc, gal21, r_last_update_pbc, iparticle, jparticle, & + cell, particle_set, mm_section) + + REAL(KIND=dp), INTENT(OUT) :: pot_loc + TYPE(gal21_pot_type), POINTER :: gal21 + TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc + INTEGER, INTENT(IN) :: iparticle, jparticle + TYPE(cell_type), POINTER :: cell + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(section_vals_type), POINTER :: mm_section + + CHARACTER(LEN=*), PARAMETER :: routineN = 'gal21_energy', routineP = moduleN//':'//routineN + + CHARACTER(LEN=2) :: element_symbol + INTEGER :: index_outfile + REAL(KIND=dp) :: anglepart, AO, BO, bxy, bz, cosalpha, & + drji2, eps, nvec(3), rji(3), sinalpha, & + sum_weight, Vang, Vgaussian, VH, VTT, & + weight + TYPE(cp_logger_type), POINTER :: logger + + pot_loc = 0.0_dp + CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, & + element_symbol=element_symbol) !Read the atom type of i + + IF (element_symbol == "O") THEN !To avoid counting two times each pair + + rji(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell) !Vector in pbc from j to i + + IF (.NOT. ALLOCATED(gal21%n_vectors)) THEN !First calling of the forcefield only + ALLOCATE (gal21%n_vectors(3, SIZE(particle_set))) + gal21%n_vectors(:, :) = 0.0_dp + END IF + + IF (gal21%express) THEN + logger => cp_get_default_logger() + index_outfile = cp_print_key_unit_nr(logger, mm_section, & + "PRINT%PROGRAM_RUN_INFO", extension=".mmLog") + IF (index_outfile > 0) WRITE (index_outfile, *) "GCN", gal21%gcn(jparticle) + CALL cp_print_key_finished_output(index_outfile, logger, mm_section, & + "PRINT%PROGRAM_RUN_INFO") + END IF + + !Build epsilon attraction and the parameters of the gaussian attraction as a function of gcn + eps = gal21%epsilon1*gal21%gcn(jparticle)*gal21%gcn(jparticle) + gal21%epsilon2*gal21%gcn(jparticle) + gal21%epsilon3 + bxy = gal21%bxy1 + gal21%bxy2*gal21%gcn(jparticle) + bz = gal21%bz1 + gal21%bz2*gal21%gcn(jparticle) + + !Angular dependance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + Vang = 0.0_dp + + !Calculation of the normal vector centered on the Me atom of the pair, only the first time that an interaction with the metal atom of the pair is evaluated + IF (gal21%n_vectors(1, jparticle) == 0.0_dp .AND. & + gal21%n_vectors(2, jparticle) == 0.0_dp .AND. & + gal21%n_vectors(3, jparticle) == 0.0_dp) THEN + gal21%n_vectors(:, jparticle) = normale(gal21, r_last_update_pbc, jparticle, & + particle_set, cell) + END IF + + nvec(:) = gal21%n_vectors(:, jparticle) !Else, retrive it, should not have moved sinc metal is supposed to be frozen + + !Calculation of the sum of the expontial weights of each Me surrounding the principal one + sum_weight = somme(gal21, r_last_update_pbc, iparticle, particle_set, cell) + + !Exponential damping weight for angular dependance + weight = EXP(-SQRT(DOT_PRODUCT(rji, rji))/gal21%r1) + + !Calculation of the truncated fourier series of the water-dipole/surface-normal angle + anglepart = 0.0_dp + VH = 0.0_dp + CALL angular(anglepart, VH, gal21, r_last_update_pbc, iparticle, jparticle, cell, particle_set, nvec, & + .TRUE., mm_section) + + !Build the complete angular potential while avoiding division by 0 + IF (weight /= 0) THEN + Vang = weight*weight*anglepart/sum_weight + IF (gal21%express) THEN + logger => cp_get_default_logger() + index_outfile = cp_print_key_unit_nr(logger, mm_section, & + "PRINT%PROGRAM_RUN_INFO", extension=".mmLog") + IF (index_outfile > 0) WRITE (index_outfile, *) "Fermi", weight*weight/sum_weight + CALL cp_print_key_finished_output(index_outfile, logger, mm_section, & + "PRINT%PROGRAM_RUN_INFO") + END IF + END IF + !END Angular%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + !Attractive Gaussian %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + Vgaussian = 0.0_dp + drji2 = DOT_PRODUCT(rji, rji) + !Alpha is the angle of the Me-O vector with the normale vector. Used for gaussian attaction + + cosalpha = DOT_PRODUCT(rji, nvec)/SQRT(drji2) + IF (cosalpha < -1.0_dp) cosalpha = -1.0_dp + IF (cosalpha > +1.0_dp) cosalpha = +1.0_dp + sinalpha = SIN(ACOS(cosalpha)) + + !Gaussian component of the energy + Vgaussian = -1.0_dp*eps*EXP(-bz*drji2*cosalpha*cosalpha & + - bxy*drji2*sinalpha*sinalpha) + !END Gaussian%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + AO = gal21%AO1 + gal21%AO2*gal21%gcn(jparticle) + BO = gal21%BO1 + gal21%BO2*gal21%gcn(jparticle) + + !Tang and toennies potential for physisorption + VTT = AO*EXP(-BO*SQRT(drji2)) - (1.0 - EXP(-BO*SQRT(drji2)) & + - BO*SQRT(drji2)*EXP(-BO*SQRT(drji2)) & + - (((BO*SQRT(drji2))**2)/2)*EXP(-BO*SQRT(drji2)) & + - (((BO*SQRT(drji2))**3)/6)*EXP(-BO*SQRT(drji2)) & + - (((BO*SQRT(drji2))**4)/24)*EXP(-BO*SQRT(drji2)) & + - (((BO*SQRT(drji2))**5)/120)*EXP(-BO*SQRT(drji2)) & + - (((BO*SQRT(drji2))**6)/720)*EXP(-BO*SQRT(drji2))) & + *gal21%c/(SQRT(drji2)**6) + + !For fit purpose only + IF (gal21%express) THEN + logger => cp_get_default_logger() + index_outfile = cp_print_key_unit_nr(logger, mm_section, & + "PRINT%PROGRAM_RUN_INFO", extension=".mmLog") + IF (index_outfile > 0) WRITE (index_outfile, *) "Gau", -1.0_dp*EXP(-bz*drji2*cosalpha*cosalpha & + - bxy*drji2*sinalpha*sinalpha) + IF (weight == 0 .AND. index_outfile > 0) WRITE (index_outfile, *) "Fermi 0" + IF (index_outfile > 0) WRITE (index_outfile, *) "expO", EXP(-BO*SQRT(drji2)) + IF (index_outfile > 0) WRITE (index_outfile, *) "cstpart", -(1.0 - EXP(-BO*SQRT(drji2)) & + - BO*SQRT(drji2)*EXP(-BO*SQRT(drji2)) & + - (((BO*SQRT(drji2))**2)/2)*EXP(-BO*SQRT(drji2)) & + - (((BO*SQRT(drji2))**3)/6)*EXP(-BO*SQRT(drji2)) & + - (((BO*SQRT(drji2))**4)/24)*EXP(-BO*SQRT(drji2)) & + - (((BO*SQRT(drji2))**5)/120)*EXP(-BO*SQRT(drji2)) & + - (((BO*SQRT(drji2))**6)/720)*EXP(-BO*SQRT(drji2))) & + *gal21%c/(SQRT(drji2)**6) + IF (index_outfile > 0) WRITE (index_outfile, *) "params_lin_eps", gal21%epsilon1, gal21%epsilon2, gal21%epsilon3 + IF (index_outfile > 0) WRITE (index_outfile, *) "params_lin_A0", AO + CALL cp_print_key_finished_output(index_outfile, logger, mm_section, & + "PRINT%PROGRAM_RUN_INFO") + END IF + !Compute the total energy + pot_loc = Vgaussian + Vang + VTT + VH + + END IF + + END SUBROUTINE gal21_energy + +! ************************************************************************************************** +! The idea is to build a vector normal to the local surface by using the symetry of the surface that +! make the opposite vectors compensate themself. The vector is therefore in the direction of the +! missing atoms of a large coordination sphere +! ************************************************************************************************** +!> \brief ... +!> \param gal21 ... +!> \param r_last_update_pbc ... +!> \param jparticle ... +!> \param particle_set ... +!> \param cell ... +!> \return ... +!> \retval normale ... +! ************************************************************************************************** + FUNCTION normale(gal21, r_last_update_pbc, jparticle, particle_set, cell) + TYPE(gal21_pot_type), POINTER :: gal21 + TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc + INTEGER, INTENT(IN) :: jparticle + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(cell_type), POINTER :: cell + REAL(KIND=dp) :: normale(3) + + CHARACTER(LEN=*), PARAMETER :: routineN = 'normale', routineP = moduleN//':'//routineN + + CHARACTER(LEN=2) :: element_symbol_k + INTEGER :: kparticle, natom + REAL(KIND=dp) :: drjk, rjk(3) + + natom = SIZE(particle_set) + normale(:) = 0.0_dp + + DO kparticle = 1, natom !Loop on every atom of the system + IF (kparticle == jparticle) CYCLE !Avoid the principal Me atom (j) in the counting + CALL get_atomic_kind(atomic_kind=particle_set(kparticle)%atomic_kind, & + element_symbol=element_symbol_k) + IF (element_symbol_k /= gal21%met1 .AND. element_symbol_k /= gal21%met2) CYCLE !Keep only metals + rjk(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(kparticle)%r(:), cell) + drjk = SQRT(DOT_PRODUCT(rjk, rjk)) + IF (drjk > gal21%rcutsq) CYCLE !Keep only those within square root of the force-field cutoff distance of the metallic atom of the evaluated pair + normale(:) = normale(:) - rjk(:)/(drjk*drjk*drjk*drjk*drjk) !Build the normal, vector by vector + END DO + + ! Normalisation of the vector + normale(:) = normale(:)/SQRT(DOT_PRODUCT(normale, normale)) + + END FUNCTION normale + +! ************************************************************************************************** +! Scan all the Me atoms that have been counted in the O-Me paires and sum their exponential weights +! ************************************************************************************************** +!> \brief ... +!> \param gal21 ... +!> \param r_last_update_pbc ... +!> \param iparticle ... +!> \param particle_set ... +!> \param cell ... +!> \return ... +!> \retval somme ... +! ************************************************************************************************** + FUNCTION somme(gal21, r_last_update_pbc, iparticle, particle_set, cell) + TYPE(gal21_pot_type), POINTER :: gal21 + TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc + INTEGER, INTENT(IN) :: iparticle + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(cell_type), POINTER :: cell + REAL(KIND=dp) :: somme + + CHARACTER(LEN=*), PARAMETER :: routineN = 'somme', routineP = moduleN//':'//routineN + + CHARACTER(LEN=2) :: element_symbol_k + INTEGER :: kparticle, natom + REAL(KIND=dp) :: rki(3) + + natom = SIZE(particle_set) + somme = 0.0_dp + + DO kparticle = 1, natom !Loop on every atom of the system + CALL get_atomic_kind(atomic_kind=particle_set(kparticle)%atomic_kind, & + element_symbol=element_symbol_k) + IF (element_symbol_k /= gal21%met1 .AND. element_symbol_k /= gal21%met2) CYCLE !Keep only metals + rki(:) = pbc(r_last_update_pbc(kparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell) + IF (SQRT(DOT_PRODUCT(rki, rki)) > gal21%rcutsq) CYCLE !Keep only those within cutoff distance of the oxygen atom of the evaluated pair (the omega ensemble) + IF (element_symbol_k == gal21%met1) somme = somme + EXP(-SQRT(DOT_PRODUCT(rki, rki))/gal21%r1) !Build the sum of the exponential weights + IF (element_symbol_k == gal21%met2) somme = somme + EXP(-SQRT(DOT_PRODUCT(rki, rki))/gal21%r2) !Build the sum of the exponential weights + END DO + + END FUNCTION somme + +! ************************************************************************************************** + +! ************************************************************************************************** +! Compute the angular dependance (on theta) of the forcefield +! ************************************************************************************************** +!> \brief ... +!> \param anglepart ... +!> \param VH ... +!> \param gal21 ... +!> \param r_last_update_pbc ... +!> \param iparticle ... +!> \param jparticle ... +!> \param cell ... +!> \param particle_set ... +!> \param nvec ... +!> \param energy ... +!> \param mm_section ... +!> \return ... +!> \retval angular ... +! ************************************************************************************************** + SUBROUTINE angular(anglepart, VH, gal21, r_last_update_pbc, iparticle, jparticle, cell, & + particle_set, nvec, energy, mm_section) + REAL(KIND=dp) :: anglepart, VH + TYPE(gal21_pot_type), POINTER :: gal21 + TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc + INTEGER, INTENT(IN) :: iparticle, jparticle + TYPE(cell_type), POINTER :: cell + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + REAL(KIND=dp), DIMENSION(3) :: nvec + LOGICAL :: energy + TYPE(section_vals_type), POINTER :: mm_section + + CHARACTER(LEN=*), PARAMETER :: routineN = 'angular', routineP = moduleN//':'//routineN + + CHARACTER(LEN=2) :: element_symbol + INTEGER :: count_h, iatom, index_h1, index_h2, & + index_outfile, natom + REAL(KIND=dp) :: a1, a2, a3, a4, BH, costheta, & + h_max_dist, rih(3), rih1(3), rih2(3), & + rix(3), rjh1(3), rjh2(3), theta + TYPE(cp_logger_type), POINTER :: logger + + count_h = 0 + index_h1 = 0 + index_h2 = 0 + h_max_dist = 2.1_dp ! 1.1 angstrom + natom = SIZE(particle_set) + + DO iatom = 1, natom !Loop on every atom of the system + CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, & + element_symbol=element_symbol) + IF (element_symbol /= "H") CYCLE !Kepp only hydrogen + rih(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(iatom)%r(:), cell) + IF (SQRT(DOT_PRODUCT(rih, rih)) >= h_max_dist) CYCLE !Keep only hydrogen that are bounded to the considered O + count_h = count_h + 1 + IF (count_h == 1) THEN + index_h1 = iatom + ELSEIF (count_h == 2) THEN + index_h2 = iatom + ENDIF + ENDDO + + ! Abort if the oxygen is not part of a water molecule (2 H) + IF (count_h /= 2) THEN + CALL cp_abort(__LOCATION__, & + " Error: Found "//cp_to_string(count_h)//" H atoms for O atom "//cp_to_string(iparticle)) + ENDIF + + a1 = gal21%a11 + gal21%a12*gal21%gcn(jparticle) + gal21%a13*gal21%gcn(jparticle)*gal21%gcn(jparticle) + a2 = gal21%a21 + gal21%a22*gal21%gcn(jparticle) + gal21%a23*gal21%gcn(jparticle)*gal21%gcn(jparticle) + a3 = gal21%a31 + gal21%a32*gal21%gcn(jparticle) + gal21%a33*gal21%gcn(jparticle)*gal21%gcn(jparticle) + a4 = gal21%a41 + gal21%a42*gal21%gcn(jparticle) + gal21%a43*gal21%gcn(jparticle)*gal21%gcn(jparticle) + + rih1(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell) + rih2(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell) + rix(:) = rih1(:) + rih2(:) ! build the dipole vector rix of the H2O molecule + costheta = DOT_PRODUCT(rix, nvec)/SQRT(DOT_PRODUCT(rix, rix)) + IF (costheta < -1.0_dp) costheta = -1.0_dp + IF (costheta > +1.0_dp) costheta = +1.0_dp + theta = ACOS(costheta) ! Theta is the angle between the normal to the surface and the dipole + anglepart = a1*costheta + a2*COS(2.0_dp*theta) + a3*COS(3.0_dp*theta) & + + a4*COS(4.0_dp*theta) ! build the fourier series + + BH = gal21%BH1 + gal21%gcn(jparticle)*gal21%BH2 + + rjh1(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell) + rjh2(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell) + VH = (gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)*(EXP(-BH*SQRT(DOT_PRODUCT(rjh1, rjh1))) + & + EXP(-BH*SQRT(DOT_PRODUCT(rjh2, rjh2)))) + + ! For fit purpose + IF (gal21%express .AND. energy) THEN + logger => cp_get_default_logger() + index_outfile = cp_print_key_unit_nr(logger, mm_section, & + "PRINT%PROGRAM_RUN_INFO", extension=".mmLog") + + IF (index_outfile > 0) WRITE (index_outfile, *) "Fourier", costheta, COS(2.0_dp*theta), COS(3.0_dp*theta), & + COS(4.0_dp*theta) !, theta + IF (index_outfile > 0) WRITE (index_outfile, *) "H_rep", EXP(-BH*SQRT(DOT_PRODUCT(rjh1, rjh1))) + & + EXP(-BH*SQRT(DOT_PRODUCT(rjh2, rjh2))) + !IF (index_outfile > 0) WRITE (index_outfile, *) "H_r6", -1/DOT_PRODUCT(rjh1,rjh1)**3 -1/DOT_PRODUCT(rjh2,rjh2)**3 + + CALL cp_print_key_finished_output(index_outfile, logger, mm_section, & + "PRINT%PROGRAM_RUN_INFO") + ENDIF + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief forces generated by the GAL2119 potential +!> \param gal21 all parameters of GAL2119 +!> \param r_last_update_pbc position of every atoms on previous frame +!> \param iparticle first index of the atom of the evaluated pair +!> \param jparticle second index of the atom of the evaluated pair +!> \param f_nonbond all the forces applying on the system +!> \param pv_nonbond ... +!> \param use_virial request of usage of virial (for barostat) +!> \param cell dimension of the pbc cell +!> \param particle_set full list of atoms of the system +!> \author Clabaut Paul - 2019 - ENS de Lyon +! ************************************************************************************************** + SUBROUTINE gal21_forces(gal21, r_last_update_pbc, iparticle, jparticle, f_nonbond, pv_nonbond, & + use_virial, cell, particle_set) + TYPE(gal21_pot_type), POINTER :: gal21 + TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc + INTEGER, INTENT(IN) :: iparticle, jparticle + REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond, pv_nonbond + LOGICAL, INTENT(IN) :: use_virial + TYPE(cell_type), POINTER :: cell + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + + CHARACTER(LEN=*), PARAMETER :: routineN = 'gal21_forces', routineP = moduleN//':'//routineN + + CHARACTER(LEN=2) :: element_symbol + REAL(KIND=dp) :: anglepart, AO, BO, bxy, bz, cosalpha, dGauss(3), drji, drjicosalpha(3), & + drjisinalpha(3), dTT(3), dweight(3), eps, nvec(3), prefactor, rji(3), rji_hat(3), & + sinalpha, sum_weight, Vgaussian, VH, weight + TYPE(section_vals_type), POINTER :: mm_section + + CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, & + element_symbol=element_symbol) + + IF (element_symbol == "O") THEN !To avoid counting two times each pair + + rji(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell) + drji = SQRT(DOT_PRODUCT(rji, rji)) + rji_hat(:) = rji(:)/drji ! hat = pure directional component of a given vector + + IF (.NOT. ALLOCATED(gal21%n_vectors)) THEN !First calling of the forcefield only + ALLOCATE (gal21%n_vectors(3, SIZE(particle_set))) + gal21%n_vectors(:, :) = 0.0_dp + END IF + + !Build epsilon attraction and the a parameters of the Fourier serie as quadratic fucntion of gcn + eps = gal21%epsilon1*gal21%gcn(jparticle)*gal21%gcn(jparticle) + gal21%epsilon2*gal21%gcn(jparticle) + gal21%epsilon3 + bxy = gal21%bxy1 + gal21%bxy2*gal21%gcn(jparticle) + bz = gal21%bz1 + gal21%bz2*gal21%gcn(jparticle) + + !Angular dependance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + !Calculation of the normal vector centered on the Me atom of the pair, only the first time that an interaction with the metal atom of the pair is evaluated + IF (gal21%n_vectors(1, jparticle) == 0.0_dp .AND. & + gal21%n_vectors(2, jparticle) == 0.0_dp .AND. & + gal21%n_vectors(3, jparticle) == 0.0_dp) THEN + gal21%n_vectors(:, jparticle) = normale(gal21, r_last_update_pbc, jparticle, & + particle_set, cell) + END IF + + nvec(:) = gal21%n_vectors(:, jparticle) !Else, retrive it, should not have moved sinc metal is supposed to be frozen + + !Calculation of the sum of the expontial weights of each Me surrounding the principal one + sum_weight = somme(gal21, r_last_update_pbc, iparticle, particle_set, cell) + + !Exponential damping weight for angular dependance + weight = EXP(-drji/gal21%r1) + dweight(:) = 1.0_dp/gal21%r1*weight*rji_hat(:) !Derivativ of it + + !Calculation of the truncated fourier series of the water-dipole/surface-normal angle + NULLIFY (mm_section) + anglepart = 0.0_dp + VH = 0.0_dp + CALL angular(anglepart, VH, gal21, r_last_update_pbc, iparticle, jparticle, cell, particle_set, nvec, & + .FALSE., mm_section) + + !Build the average of the exponential weight while avoiding division by 0 + IF (weight /= 0) THEN + ! Calculate the first component of the derivativ of the angular term + f_nonbond(1:3, iparticle) = f_nonbond(1:3, iparticle) + 2.0_dp*dweight(1:3)*weight* & + anglepart/sum_weight + + IF (use_virial) THEN + pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rji(1)*2.0_dp*dweight(1:3)*weight* & + anglepart/sum_weight + pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rji(2)*2.0_dp*dweight(1:3)*weight* & + anglepart/sum_weight + pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rji(3)*2.0_dp*dweight(1:3)*weight* & + anglepart/sum_weight + END IF + + ! Calculate the second component of the derivativ of the angular term + CALL somme_d(gal21, r_last_update_pbc, iparticle, jparticle, & + f_nonbond, pv_nonbond, use_virial, particle_set, cell, anglepart, sum_weight) + + prefactor = (-1.0_dp)*weight*weight/sum_weight ! Avoiding division by 0 + + ! Calculate the third component of the derivativ of the angular term + CALL angular_d(gal21, r_last_update_pbc, iparticle, jparticle, & + f_nonbond, pv_nonbond, use_virial, prefactor, cell, particle_set, nvec) + END IF + !END Angular%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + !Attractive Gaussian %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + !Alpha is the angle of the Me-O vector with the normale vector. Used for gaussian attaction + cosalpha = DOT_PRODUCT(rji, nvec)/drji + IF (cosalpha < -1.0_dp) cosalpha = -1.0_dp + IF (cosalpha > +1.0_dp) cosalpha = +1.0_dp + sinalpha = SIN(ACOS(cosalpha)) + + !Gaussian component of the energy + Vgaussian = -1.0_dp*eps*EXP(-bz*DOT_PRODUCT(rji, rji)*cosalpha*cosalpha & + - bxy*DOT_PRODUCT(rji, rji)*sinalpha*sinalpha) + + ! Calculation of partial derivativ of the gaussian components + drjicosalpha(:) = rji_hat(:)*cosalpha + nvec(:) - cosalpha*rji_hat(:) + drjisinalpha(:) = rji_hat(:)*sinalpha - (cosalpha/sinalpha)*(nvec(:) - cosalpha*rji_hat(:)) + dGauss(:) = (-1.0_dp*bz*2*drji*cosalpha*drjicosalpha - & + 1.0_dp*bxy*2*drji*sinalpha*drjisinalpha)*Vgaussian*(-1.0_dp) + + ! Force due to gaussian term + f_nonbond(1:3, iparticle) = f_nonbond(1:3, iparticle) + dGauss(1:3) + + IF (use_virial) THEN + pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rji(1)*dGauss(1:3) + pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rji(2)*dGauss(1:3) + pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rji(3)*dGauss(1:3) + END IF + !END Gaussian%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + AO = gal21%AO1 + gal21%AO2*gal21%gcn(jparticle) + BO = gal21%BO1 + gal21%BO2*gal21%gcn(jparticle) + + !Derivativ of the Tang and Toennies term + dTT(:) = (-(AO*BO + (BO**7)*gal21%c/720)*EXP(-BO*drji) + 6*(gal21%c/drji**7)* & + (1.0 - EXP(-BO*drji) & + - BO*drji*EXP(-BO*drji) & + - (((BO*drji)**2)/2)*EXP(-BO*drji) & + - (((BO*drji)**3)/6)*EXP(-BO*drji) & + - (((BO*drji)**4)/24)*EXP(-BO*drji) & + - (((BO*drji)**5)/120)*EXP(-BO*drji) & + - (((BO*drji)**6)/720)*EXP(-BO*drji)) & + )*rji_hat(:) + + ! Force of Tang & Toennies + f_nonbond(1:3, iparticle) = f_nonbond(1:3, iparticle) - dTT(1:3) + + IF (use_virial) THEN + pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) - rji(1)*dTT(1:3) + pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) - rji(2)*dTT(1:3) + pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) - rji(3)*dTT(1:3) + END IF + + END IF + + END SUBROUTINE gal21_forces +! ************************************************************************************************** +! Derivativ of the second component of angular dependance +! ************************************************************************************************** + +! ************************************************************************************************** +!> \brief ... +!> \param gal21 ... +!> \param r_last_update_pbc ... +!> \param iparticle ... +!> \param jparticle ... +!> \param f_nonbond ... +!> \param pv_nonbond ... +!> \param use_virial ... +!> \param particle_set ... +!> \param cell ... +!> \param anglepart ... +!> \param sum_weight ... +! ************************************************************************************************** + SUBROUTINE somme_d(gal21, r_last_update_pbc, iparticle, jparticle, & + f_nonbond, pv_nonbond, use_virial, particle_set, cell, anglepart, sum_weight) + TYPE(gal21_pot_type), POINTER :: gal21 + TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc + INTEGER, INTENT(IN) :: iparticle, jparticle + REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond, pv_nonbond + LOGICAL, INTENT(IN) :: use_virial + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(cell_type), POINTER :: cell + REAL(KIND=dp), INTENT(IN) :: anglepart, sum_weight + + CHARACTER(LEN=*), PARAMETER :: routineN = 'somme_d', routineP = moduleN//':'//routineN + + CHARACTER(LEN=2) :: element_symbol_k + INTEGER :: kparticle, natom + REAL(KIND=dp) :: drki, dwdr(3), rji(3), rki(3), & + rki_hat(3), weight_rji + + rji(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell) + weight_rji = EXP(-SQRT(DOT_PRODUCT(rji, rji))/gal21%r1) + + natom = SIZE(particle_set) + DO kparticle = 1, natom !Loop on every atom of the system + CALL get_atomic_kind(atomic_kind=particle_set(kparticle)%atomic_kind, & + element_symbol=element_symbol_k) + IF (element_symbol_k /= gal21%met1 .AND. element_symbol_k /= gal21%met2) CYCLE !Keep only metals + rki(:) = pbc(r_last_update_pbc(kparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell) + IF (SQRT(DOT_PRODUCT(rki, rki)) > gal21%rcutsq) CYCLE !Keep only those within cutoff distance of the oxygen atom of the evaluated pair (the omega ensemble) + drki = SQRT(DOT_PRODUCT(rki, rki)) + rki_hat(:) = rki(:)/drki + + IF (element_symbol_k == gal21%met1) dwdr(:) = (-1.0_dp)*(1.0_dp/gal21%r1)*EXP(-drki/gal21%r1)*rki_hat(:) !Build the sum of derivativs + IF (element_symbol_k == gal21%met2) dwdr(:) = (-1.0_dp)*(1.0_dp/gal21%r2)*EXP(-drki/gal21%r2)*rki_hat(:) + + f_nonbond(1:3, iparticle) = f_nonbond(1:3, iparticle) + dwdr(1:3)*weight_rji & + *weight_rji*anglepart/(sum_weight**2) + + IF (use_virial) THEN + pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rki(1)*dwdr(1:3)*weight_rji & + *weight_rji*anglepart/(sum_weight**2) + pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rki(2)*dwdr(1:3)*weight_rji & + *weight_rji*anglepart/(sum_weight**2) + pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rki(3)*dwdr(1:3)*weight_rji & + *weight_rji*anglepart/(sum_weight**2) + END IF + + END DO + + END SUBROUTINE somme_d + +! ************************************************************************************************** +! Derivativ of the third component of angular term +! ************************************************************************************************** +!> \brief ... +!> \param gal21 ... +!> \param r_last_update_pbc ... +!> \param iparticle ... +!> \param jparticle ... +!> \param f_nonbond ... +!> \param pv_nonbond ... +!> \param use_virial ... +!> \param prefactor ... +!> \param cell ... +!> \param particle_set ... +!> \param nvec ... +! ************************************************************************************************** + SUBROUTINE angular_d(gal21, r_last_update_pbc, iparticle, jparticle, f_nonbond, & + pv_nonbond, use_virial, prefactor, cell, particle_set, nvec) + TYPE(gal21_pot_type), POINTER :: gal21 + TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc + INTEGER, INTENT(IN) :: iparticle, jparticle + REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond, pv_nonbond + LOGICAL, INTENT(IN) :: use_virial + REAL(KIND=dp), INTENT(IN) :: prefactor + TYPE(cell_type), POINTER :: cell + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + REAL(KIND=dp), DIMENSION(3) :: nvec + + CHARACTER(LEN=*), PARAMETER :: routineN = 'angular_d', routineP = moduleN//':'//routineN + + CHARACTER(LEN=2) :: element_symbol + INTEGER :: count_h, iatom, index_h1, index_h2, natom + REAL(KIND=dp) :: a1, a2, a3, a4, BH, costheta, & + dsumdtheta, h_max_dist, theta + REAL(KIND=dp), DIMENSION(3) :: dangular, dcostheta, rih, rih1, rih2, & + rix, rix_hat, rjh1, rjh2, rji, rji_hat + + count_h = 0 + index_h1 = 0 + index_h2 = 0 + h_max_dist = 2.1_dp ! 1.1 angstrom + natom = SIZE(particle_set) + + DO iatom = 1, natom !Loop on every atom of the system + CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, & + element_symbol=element_symbol) + IF (element_symbol /= "H") CYCLE !Kepp only hydrogen + rih(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(iatom)%r(:), cell) + IF (SQRT(DOT_PRODUCT(rih, rih)) >= h_max_dist) CYCLE !Keep only hydrogen that are bounded to the considered O + count_h = count_h + 1 + IF (count_h == 1) THEN + index_h1 = iatom + ELSEIF (count_h == 2) THEN + index_h2 = iatom + ENDIF + ENDDO + + ! Abort if the oxygen is not part of a water molecule (2 H) + IF (count_h /= 2) THEN + CALL cp_abort(__LOCATION__, & + " Error: Found "//cp_to_string(count_h)//" H atoms for O atom "//cp_to_string(iparticle)) + ENDIF + + a1 = gal21%a11 + gal21%a12*gal21%gcn(jparticle) + gal21%a13*gal21%gcn(jparticle)*gal21%gcn(jparticle) + a2 = gal21%a21 + gal21%a22*gal21%gcn(jparticle) + gal21%a23*gal21%gcn(jparticle)*gal21%gcn(jparticle) + a3 = gal21%a31 + gal21%a32*gal21%gcn(jparticle) + gal21%a33*gal21%gcn(jparticle)*gal21%gcn(jparticle) + a4 = gal21%a41 + gal21%a42*gal21%gcn(jparticle) + gal21%a43*gal21%gcn(jparticle)*gal21%gcn(jparticle) + + rji(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell) + rji_hat(:) = rji(:)/SQRT(DOT_PRODUCT(rji, rji)) ! hat = pure directional component of a given vector + + !dipole vector rix of the H2O molecule + rih1(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell) + rih2(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell) + rix(:) = rih1(:) + rih2(:) ! build the dipole vector rix of the H2O molecule + rix_hat(:) = rix(:)/SQRT(DOT_PRODUCT(rix, rix)) ! hat = pure directional component of a given vector + costheta = DOT_PRODUCT(rix, nvec)/SQRT(DOT_PRODUCT(rix, rix)) ! Theta is the angle between the normal to the surface and the dipole + IF (costheta < -1.0_dp) costheta = -1.0_dp + IF (costheta > +1.0_dp) costheta = +1.0_dp + theta = ACOS(costheta) ! Theta is the angle between the normal to the surface and the dipole + + ! Calculation of partial derivativ of the angular components + dsumdtheta = -1.0_dp*a1*SIN(theta) - a2*2.0_dp*SIN(2.0_dp*theta) - & + a3*3.0_dp*SIN(3.0_dp*theta) - a4*4.0_dp*SIN(4.0_dp*theta) + dcostheta(:) = (1.0_dp/SQRT(DOT_PRODUCT(rix, rix)))*(nvec(:) - costheta*rix_hat(:)) + dangular(:) = prefactor*dsumdtheta*(-1.0_dp/SIN(theta))*dcostheta(:) + + !Force due to the third component of the derivativ of the angular term + f_nonbond(1:3, iparticle) = f_nonbond(1:3, iparticle) - dangular(1:3)*2.0_dp !(one per H) + f_nonbond(1:3, index_h1) = f_nonbond(1:3, index_h1) + dangular(1:3) + f_nonbond(1:3, index_h2) = f_nonbond(1:3, index_h2) + dangular(1:3) + + IF (use_virial) THEN + pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rix(1)*dangular(1:3) + pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rix(2)*dangular(1:3) + pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rix(3)*dangular(1:3) + END IF + + BH = gal21%BH1 + gal21%gcn(jparticle)*gal21%BH2 + + rjh1(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell) + f_nonbond(1:3, index_h1) = f_nonbond(1:3, index_h1) + (gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* & + BH*EXP(-BH*SQRT(DOT_PRODUCT(rjh1, rjh1)))*rjh1(:)/SQRT(DOT_PRODUCT(rjh1, rjh1)) + + IF (use_virial) THEN + pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rjh1(1)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* & + BH*EXP(-BH*SQRT(DOT_PRODUCT(rjh1, rjh1)))) & + *rjh1(:)/SQRT(DOT_PRODUCT(rjh1, rjh1)) + pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rjh1(2)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* & + BH*EXP(-BH*SQRT(DOT_PRODUCT(rjh1, rjh1)))) & + *rjh1(:)/SQRT(DOT_PRODUCT(rjh1, rjh1)) + pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rjh2(3)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* & + BH*EXP(-BH*SQRT(DOT_PRODUCT(rjh1, rjh1)))) & + *rjh1(:)/SQRT(DOT_PRODUCT(rjh1, rjh1)) + END IF + + rjh2(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell) + f_nonbond(1:3, index_h2) = f_nonbond(1:3, index_h2) + ((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* & + BH*EXP(-BH*SQRT(DOT_PRODUCT(rjh2, rjh2)))) & + *rjh2(:)/SQRT(DOT_PRODUCT(rjh2, rjh2)) + + IF (use_virial) THEN + pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rjh2(1)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* & + BH*EXP(-BH*SQRT(DOT_PRODUCT(rjh2, rjh2)))) & + *rjh2(:)/SQRT(DOT_PRODUCT(rjh2, rjh2)) + pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rjh2(2)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* & + BH*EXP(-BH*SQRT(DOT_PRODUCT(rjh2, rjh2)))) & + *rjh2(:)/SQRT(DOT_PRODUCT(rjh2, rjh2)) + pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rjh2(3)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* & + BH*EXP(-BH*SQRT(DOT_PRODUCT(rjh2, rjh2)))) & + *rjh2(:)/SQRT(DOT_PRODUCT(rjh2, rjh2)) + END IF + + END SUBROUTINE angular_d + +! ************************************************************************************************** +!> \brief ... +!> \param nonbonded ... +!> \param potparm ... +!> \param glob_loc_list ... +!> \param glob_cell_v ... +!> \param glob_loc_list_a ... +!> \param cell ... +!> \par History +! ************************************************************************************************** + SUBROUTINE setup_gal21_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, & + glob_loc_list_a, cell) + TYPE(fist_neighbor_type), POINTER :: nonbonded + TYPE(pair_potential_pp_type), POINTER :: potparm + INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list + REAL(KIND=dp), DIMENSION(:, :), POINTER :: glob_cell_v + INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a + TYPE(cell_type), POINTER :: cell + + CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_gal21_arrays', & + routineP = moduleN//':'//routineN + + INTEGER :: handle, i, iend, igrp, ikind, ilist, & + ipair, istart, jkind, nkinds, npairs, & + npairs_tot + INTEGER, DIMENSION(:), POINTER :: work_list, work_list2 + INTEGER, DIMENSION(:, :), POINTER :: list + REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi + REAL(KIND=dp), DIMENSION(:, :), POINTER :: rwork_list + TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair + TYPE(pair_potential_single_type), POINTER :: pot + + CPASSERT(.NOT. ASSOCIATED(glob_loc_list)) + CPASSERT(.NOT. ASSOCIATED(glob_loc_list_a)) + CPASSERT(.NOT. ASSOCIATED(glob_cell_v)) + CALL timeset(routineN, handle) + npairs_tot = 0 + nkinds = SIZE(potparm%pot, 1) + DO ilist = 1, nonbonded%nlists + neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) + npairs = neighbor_kind_pair%npairs + IF (npairs == 0) CYCLE + Kind_Group_Loop1: DO igrp = 1, neighbor_kind_pair%ngrp_kind + istart = neighbor_kind_pair%grp_kind_start(igrp) + iend = neighbor_kind_pair%grp_kind_end(igrp) + ikind = neighbor_kind_pair%ij_kind(1, igrp) + jkind = neighbor_kind_pair%ij_kind(2, igrp) + pot => potparm%pot(ikind, jkind)%pot + npairs = iend - istart + 1 + IF (pot%no_mb) CYCLE + DO i = 1, SIZE(pot%type) + IF (pot%type(i) == gal21_type) npairs_tot = npairs_tot + npairs + END DO + END DO Kind_Group_Loop1 + END DO + ALLOCATE (work_list(npairs_tot)) + ALLOCATE (work_list2(npairs_tot)) + ALLOCATE (glob_loc_list(2, npairs_tot)) + ALLOCATE (glob_cell_v(3, npairs_tot)) + ! Fill arrays with data + npairs_tot = 0 + DO ilist = 1, nonbonded%nlists + neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) + npairs = neighbor_kind_pair%npairs + IF (npairs == 0) CYCLE + Kind_Group_Loop2: DO igrp = 1, neighbor_kind_pair%ngrp_kind + istart = neighbor_kind_pair%grp_kind_start(igrp) + iend = neighbor_kind_pair%grp_kind_end(igrp) + ikind = neighbor_kind_pair%ij_kind(1, igrp) + jkind = neighbor_kind_pair%ij_kind(2, igrp) + list => neighbor_kind_pair%list + cvi = neighbor_kind_pair%cell_vector + pot => potparm%pot(ikind, jkind)%pot + npairs = iend - istart + 1 + IF (pot%no_mb) CYCLE + CALL matvec_3x3(cell_v, cell%hmat, cvi) + DO i = 1, SIZE(pot%type) + ! gal21 + IF (pot%type(i) == gal21_type) THEN + DO ipair = 1, npairs + glob_loc_list(:, npairs_tot + ipair) = list(:, istart - 1 + ipair) + glob_cell_v(1:3, npairs_tot + ipair) = cell_v(1:3) + END DO + npairs_tot = npairs_tot + npairs + END IF + END DO + END DO Kind_Group_Loop2 + END DO + ! Order the arrays w.r.t. the first index of glob_loc_list + CALL sort(glob_loc_list(1, :), npairs_tot, work_list) + DO ipair = 1, npairs_tot + work_list2(ipair) = glob_loc_list(2, work_list(ipair)) + END DO + glob_loc_list(2, :) = work_list2 + DEALLOCATE (work_list2) + ALLOCATE (rwork_list(3, npairs_tot)) + DO ipair = 1, npairs_tot + rwork_list(:, ipair) = glob_cell_v(:, work_list(ipair)) + END DO + glob_cell_v = rwork_list + DEALLOCATE (rwork_list) + DEALLOCATE (work_list) + ALLOCATE (glob_loc_list_a(npairs_tot)) + glob_loc_list_a = glob_loc_list(1, :) + CALL timestop(handle) + END SUBROUTINE setup_gal21_arrays + +! ************************************************************************************************** +!> \brief ... +!> \param glob_loc_list ... +!> \param glob_cell_v ... +!> \param glob_loc_list_a ... +! ************************************************************************************************** + SUBROUTINE destroy_gal21_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a) + INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list + REAL(KIND=dp), DIMENSION(:, :), POINTER :: glob_cell_v + INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a + + CHARACTER(LEN=*), PARAMETER :: routineN = 'destroy_gal21_arrays', & + routineP = moduleN//':'//routineN + + IF (ASSOCIATED(glob_loc_list)) THEN + DEALLOCATE (glob_loc_list) + END IF + IF (ASSOCIATED(glob_loc_list_a)) THEN + DEALLOCATE (glob_loc_list_a) + END IF + IF (ASSOCIATED(glob_cell_v)) THEN + DEALLOCATE (glob_cell_v) + END IF + + END SUBROUTINE destroy_gal21_arrays + +! ************************************************************************************************** +!> \brief prints the number of OH- ions or H3O+ ions near surface +!> \param nr_ions number of ions +!> \param mm_section ... +!> \param para_env ... +!> \param print_oh flag indicating if number OH- is printed +!> \param print_h3o flag indicating if number H3O+ is printed +!> \param print_o flag indicating if number O^(2-) is printed +! ************************************************************************************************** + SUBROUTINE print_nr_ions_gal21(nr_ions, mm_section, para_env, print_oh, & + print_h3o, print_o) + INTEGER, INTENT(INOUT) :: nr_ions + TYPE(section_vals_type), POINTER :: mm_section + TYPE(cp_para_env_type), OPTIONAL, POINTER :: para_env + LOGICAL, INTENT(IN) :: print_oh, print_h3o, print_o + + CHARACTER(LEN=*), PARAMETER :: routineN = 'print_nr_ions_gal21', & + routineP = moduleN//':'//routineN + + INTEGER :: iw + TYPE(cp_logger_type), POINTER :: logger + + NULLIFY (logger) + + CALL mp_sum(nr_ions, para_env%group) + logger => cp_get_default_logger() + + iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%PROGRAM_RUN_INFO", & + extension=".mmLog") + + IF (iw > 0 .AND. nr_ions > 0 .AND. print_oh) THEN + WRITE (iw, '(/,A,T71,I10,/)') " gal21: number of OH- ions at surface", nr_ions + ENDIF + IF (iw > 0 .AND. nr_ions > 0 .AND. print_h3o) THEN + WRITE (iw, '(/,A,T71,I10,/)') " gal21: number of H3O+ ions at surface", nr_ions + ENDIF + IF (iw > 0 .AND. nr_ions > 0 .AND. print_o) THEN + WRITE (iw, '(/,A,T71,I10,/)') " gal21: number of O^2- ions at surface", nr_ions + ENDIF + + CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%PROGRAM_RUN_INFO") + + END SUBROUTINE print_nr_ions_gal21 + +END MODULE manybody_gal21 diff --git a/src/manybody_potential.F b/src/manybody_potential.F index a9349acb79..6093d2841f 100644 --- a/src/manybody_potential.F +++ b/src/manybody_potential.F @@ -32,6 +32,10 @@ MODULE manybody_potential gal_energy,& gal_forces,& setup_gal_arrays + USE manybody_gal21, ONLY: destroy_gal21_arrays,& + gal21_energy,& + gal21_forces,& + setup_gal21_arrays USE manybody_quip, ONLY: quip_add_force_virial,& quip_energy_store_force_virial USE manybody_siepmann, ONLY: destroy_siepmann_arrays,& @@ -47,9 +51,9 @@ MODULE manybody_potential USE mathlib, ONLY: matvec_3x3 USE message_passing, ONLY: mp_sum USE pair_potential_types, ONLY: & - ea_type, eam_pot_type, gal_pot_type, gal_type, pair_potential_pp_type, & - pair_potential_single_type, quip_type, siepmann_pot_type, siepmann_type, tersoff_pot_type, & - tersoff_type + ea_type, eam_pot_type, gal21_pot_type, gal21_type, gal_pot_type, gal_type, & + pair_potential_pp_type, pair_potential_single_type, quip_type, siepmann_pot_type, & + siepmann_type, tersoff_pot_type, tersoff_type USE particle_types, ONLY: particle_type USE util, ONLY: sort #include "./base/base_uses.f90" @@ -96,8 +100,8 @@ CONTAINS nloc_size, npairs, nparticle, nparticle_local, nr_h3O, nr_o, nr_oh, nunique INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a, work_list INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list, list, sort_list - LOGICAL :: any_gal, any_quip, any_siepmann, & - any_tersoff + LOGICAL :: any_gal, any_gal21, any_quip, & + any_siepmann, any_tersoff REAL(KIND=dp) :: drij, embed, pot_loc, pot_quip, qr, & rab2_max, rij(3) REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi @@ -106,6 +110,7 @@ CONTAINS TYPE(eam_pot_type), POINTER :: eam TYPE(eam_type), DIMENSION(:), POINTER :: eam_data TYPE(fist_neighbor_type), POINTER :: nonbonded + TYPE(gal21_pot_type), POINTER :: gal21 TYPE(gal_pot_type), POINTER :: gal TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair TYPE(pair_potential_pp_type), POINTER :: potparm @@ -114,10 +119,11 @@ CONTAINS TYPE(siepmann_pot_type), POINTER :: siepmann TYPE(tersoff_pot_type), POINTER :: tersoff - NULLIFY (eam, siepmann, tersoff, gal) + NULLIFY (eam, siepmann, tersoff, gal, gal21) any_tersoff = .FALSE. any_siepmann = .FALSE. any_gal = .FALSE. + any_gal21 = .FALSE. any_quip = .FALSE. CALL timeset(routineN, handle) CALL fist_nonbond_env_get(fist_nonbond_env, r_last_update_pbc=r_last_update_pbc, & @@ -164,6 +170,7 @@ CONTAINS any_quip = any_quip .OR. ANY(pot%type == quip_type) any_siepmann = any_siepmann .OR. ANY(pot%type == siepmann_type) any_gal = any_gal .OR. ANY(pot%type == gal_type) + any_gal21 = any_gal21 .OR. ANY(pot%type == gal21_type) END DO END DO CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, natom_types=nkinds) @@ -421,6 +428,86 @@ CONTAINS CALL destroy_gal_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a) END IF + !GAL21 POTENTIAL + IF (any_gal21) THEN + NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a) + CALL setup_gal21_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell) + DO ilist = 1, nonbonded%nlists + neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) + npairs = neighbor_kind_pair%npairs + IF (npairs == 0) CYCLE + Kind_Group_Loop_5: DO igrp = 1, neighbor_kind_pair%ngrp_kind + istart = neighbor_kind_pair%grp_kind_start(igrp) + iend = neighbor_kind_pair%grp_kind_end(igrp) + ikind = neighbor_kind_pair%ij_kind(1, igrp) + jkind = neighbor_kind_pair%ij_kind(2, igrp) + list => neighbor_kind_pair%list + cvi = neighbor_kind_pair%cell_vector + pot => potparm%pot(ikind, jkind)%pot + DO i = 1, SIZE(pot%type) + IF (pot%type(i) /= gal21_type) CYCLE + rab2_max = pot%set(i)%gal21%rcutsq + CALL matvec_3x3(cell_v, cell%hmat, cvi) + pot => potparm%pot(ikind, jkind)%pot + gal21 => pot%set(i)%gal21 + npairs = iend - istart + 1 + IF (npairs /= 0) THEN + ALLOCATE (sort_list(2, npairs), work_list(npairs)) + sort_list = list(:, istart:iend) + ! Sort the list of neighbors, this increases the efficiency for single + ! potential contributions + CALL sort(sort_list(1, :), npairs, work_list) + DO ipair = 1, npairs + work_list(ipair) = sort_list(2, work_list(ipair)) + END DO + sort_list(2, :) = work_list + ! find number of unique elements of array index 1 + nunique = 1 + DO ipair = 1, npairs - 1 + IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1 + END DO + ipair = 1 + junique = sort_list(1, ipair) + ifirst = 1 + DO iunique = 1, nunique + atom_a = junique + IF (glob_loc_list_a(ifirst) > atom_a) CYCLE + DO mpair = ifirst, SIZE(glob_loc_list_a) + IF (glob_loc_list_a(mpair) == atom_a) EXIT + END DO + ifirst = mpair + DO mpair = ifirst, SIZE(glob_loc_list_a) + IF (glob_loc_list_a(mpair) /= atom_a) EXIT + END DO + ilast = mpair - 1 + nloc_size = 0 + IF (ifirst /= 0) nloc_size = ilast - ifirst + 1 + DO WHILE (ipair <= npairs) + IF (sort_list(1, ipair) /= junique) EXIT + atom_b = sort_list(2, ipair) + ! Energy terms + pot_loc = 0.0_dp + rij(:) = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v + drij = DOT_PRODUCT(rij, rij) + ipair = ipair + 1 + IF (drij > rab2_max) CYCLE + drij = SQRT(drij) + CALL gal21_energy(pot_loc, gal21, r_last_update_pbc, atom_a, atom_b, & + cell, particle_set, mm_section) + + pot_manybody = pot_manybody + pot_loc + END DO + ifirst = ilast + 1 + IF (ipair <= npairs) junique = sort_list(1, ipair) + END DO + DEALLOCATE (sort_list, work_list) + END IF + END DO + END DO Kind_Group_Loop_5 + END DO + CALL destroy_gal21_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a) + END IF + CALL timestop(handle) END SUBROUTINE energy_manybody @@ -453,8 +540,8 @@ CONTAINS INTEGER, ALLOCATABLE, DIMENSION(:, :) :: eam_kinds_index INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a, work_list INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list, list, sort_list - LOGICAL :: any_gal, any_quip, any_siepmann, & - any_tersoff + LOGICAL :: any_gal, any_gal21, any_quip, & + any_siepmann, any_tersoff REAL(KIND=dp) :: f_eam, fac, fr(3), ptens11, ptens12, ptens13, ptens21, ptens22, ptens23, & ptens31, ptens32, ptens33, rab(3), rab2, rab2_max, rtmp(3) REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi @@ -462,6 +549,7 @@ CONTAINS TYPE(eam_pot_type), POINTER :: eam_a, eam_b TYPE(eam_type), DIMENSION(:), POINTER :: eam_data TYPE(fist_neighbor_type), POINTER :: nonbonded + TYPE(gal21_pot_type), POINTER :: gal21 TYPE(gal_pot_type), POINTER :: gal TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair TYPE(pair_potential_pp_type), POINTER :: potparm @@ -474,8 +562,9 @@ CONTAINS any_quip = .FALSE. any_siepmann = .FALSE. any_gal = .FALSE. + any_gal21 = .FALSE. CALL timeset(routineN, handle) - NULLIFY (eam_a, eam_b, tersoff, siepmann, gal) + NULLIFY (eam_a, eam_b, tersoff, siepmann, gal, gal21) CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, potparm=potparm, & natom_types=nkinds, eam_data=eam_data, r_last_update_pbc=r_last_update_pbc) @@ -527,6 +616,7 @@ CONTAINS any_tersoff = any_tersoff .OR. ANY(pot%type == tersoff_type) any_siepmann = any_siepmann .OR. ANY(pot%type == siepmann_type) any_gal = any_gal .OR. ANY(pot%type == gal_type) + any_gal21 = any_gal21 .OR. ANY(pot%type == gal21_type) i = eam_kinds_index(ikind, jkind) IF (i == -1) CYCLE ! EAM @@ -816,6 +906,84 @@ CONTAINS CALL destroy_gal_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a) END IF + ! GAL21 potential.. + IF (any_gal21) THEN + NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a) + CALL setup_gal21_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell) + DO ilist = 1, nonbonded%nlists + neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) + npairs = neighbor_kind_pair%npairs + IF (npairs == 0) CYCLE + Kind_Group_Loop6: DO igrp = 1, neighbor_kind_pair%ngrp_kind + istart = neighbor_kind_pair%grp_kind_start(igrp) + iend = neighbor_kind_pair%grp_kind_end(igrp) + ikind = neighbor_kind_pair%ij_kind(1, igrp) + jkind = neighbor_kind_pair%ij_kind(2, igrp) + list => neighbor_kind_pair%list + cvi = neighbor_kind_pair%cell_vector + pot => potparm%pot(ikind, jkind)%pot + + IF (pot%no_mb) CYCLE + rab2_max = pot%rcutsq + CALL matvec_3x3(cell_v, cell%hmat, cvi) + DO i = 1, SIZE(pot%type) + ! GAL21 + IF (pot%type(i) == gal21_type) THEN + npairs = iend - istart + 1 + gal21 => pot%set(i)%gal21 + ALLOCATE (sort_list(2, npairs), work_list(npairs)) + sort_list = list(:, istart:iend) + ! Sort the list of neighbors, this increases the efficiency for single + ! potential contributions + CALL sort(sort_list(1, :), npairs, work_list) + DO ipair = 1, npairs + work_list(ipair) = sort_list(2, work_list(ipair)) + END DO + sort_list(2, :) = work_list + ! find number of unique elements of array index 1 + nunique = 1 + DO ipair = 1, npairs - 1 + IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1 + END DO + ipair = 1 + junique = sort_list(1, ipair) + ifirst = 1 + DO iunique = 1, nunique + atom_a = junique + IF (glob_loc_list_a(ifirst) > atom_a) CYCLE + DO mpair = ifirst, SIZE(glob_loc_list_a) + IF (glob_loc_list_a(mpair) == atom_a) EXIT + END DO + ifirst = mpair + DO mpair = ifirst, SIZE(glob_loc_list_a) + IF (glob_loc_list_a(mpair) /= atom_a) EXIT + END DO + ilast = mpair - 1 + nloc_size = 0 + IF (ifirst /= 0) nloc_size = ilast - ifirst + 1 + DO WHILE (ipair <= npairs) + IF (sort_list(1, ipair) /= junique) EXIT + atom_b = sort_list(2, ipair) + ! Derivative terms + rtmp = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v + ipair = ipair + 1 + IF (DOT_PRODUCT(rtmp, rtmp) <= gal21%rcutsq) THEN + CALL gal21_forces(gal21, r_last_update_pbc, & + atom_a, atom_b, f_nonbond, pv_nonbond, use_virial, & + cell, particle_set) + END IF + END DO + ifirst = ilast + 1 + IF (ipair <= npairs) junique = sort_list(1, ipair) + END DO + DEALLOCATE (sort_list, work_list) + END IF + END DO + END DO Kind_Group_Loop6 + END DO + CALL destroy_gal_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a) + END IF + IF (use_virial) THEN pv_nonbond(1, 1) = pv_nonbond(1, 1) + ptens11 pv_nonbond(1, 2) = pv_nonbond(1, 2) + ptens12 diff --git a/src/pair_potential.F b/src/pair_potential.F index 087776c430..71f2d86c84 100644 --- a/src/pair_potential.F +++ b/src/pair_potential.F @@ -29,10 +29,10 @@ MODULE pair_potential default_string_length,& dp USE pair_potential_types, ONLY: & - b4_type, bm_type, compare_pot, ea_type, ft_type, ftd_type, gal_type, gp_type, gw_type, & - ip_type, list_pot, lj_charmm_type, lj_type, multi_type, nn_type, pair_potential_pp_type, & - pair_potential_single_type, potential_single_allocation, quip_type, siepmann_type, & - tersoff_type, wl_type + b4_type, bm_type, compare_pot, ea_type, ft_type, ftd_type, gal21_type, gal_type, gp_type, & + gw_type, ip_type, list_pot, lj_charmm_type, lj_type, multi_type, nn_type, & + pair_potential_pp_type, pair_potential_single_type, potential_single_allocation, & + quip_type, siepmann_type, tersoff_type, wl_type USE pair_potential_util, ONLY: ener_pot,& ener_zbl,& zbl_matching_polinomial @@ -194,6 +194,8 @@ CONTAINS pot%no_mb = .FALSE. CASE (gal_type) pot%no_mb = .FALSE. + CASE (gal21_type) + pot%no_mb = .FALSE. CASE (nn_type) ! Do nothing.. CASE DEFAULT @@ -622,6 +624,8 @@ CONTAINS nvar = 5 + nvar CASE (gal_type) nvar = 12 + nvar + CASE (gal21_type) + nvar = 30 + nvar CASE (nn_type) nvar = nvar CASE DEFAULT @@ -745,6 +749,37 @@ CONTAINS pot_par(nk, 10) = potparm%pot(i, j)%pot%set(1)%gal%a pot_par(nk, 11) = potparm%pot(i, j)%pot%set(1)%gal%b pot_par(nk, 12) = potparm%pot(i, j)%pot%set(1)%gal%c + CASE (gal21_type) + pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%gal21%epsilon1 + pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%gal21%epsilon2 + pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%gal21%epsilon3 + pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%gal21%bxy1 + pot_par(nk, 5) = potparm%pot(i, j)%pot%set(1)%gal21%bxy2 + pot_par(nk, 6) = potparm%pot(i, j)%pot%set(1)%gal21%bz1 + pot_par(nk, 7) = potparm%pot(i, j)%pot%set(1)%gal21%bz2 + pot_par(nk, 8) = potparm%pot(i, j)%pot%set(1)%gal21%r1 + pot_par(nk, 9) = potparm%pot(i, j)%pot%set(1)%gal21%r2 + pot_par(nk, 10) = potparm%pot(i, j)%pot%set(1)%gal21%a11 + pot_par(nk, 11) = potparm%pot(i, j)%pot%set(1)%gal21%a12 + pot_par(nk, 12) = potparm%pot(i, j)%pot%set(1)%gal21%a13 + pot_par(nk, 13) = potparm%pot(i, j)%pot%set(1)%gal21%a21 + pot_par(nk, 14) = potparm%pot(i, j)%pot%set(1)%gal21%a22 + pot_par(nk, 15) = potparm%pot(i, j)%pot%set(1)%gal21%a23 + pot_par(nk, 16) = potparm%pot(i, j)%pot%set(1)%gal21%a31 + pot_par(nk, 17) = potparm%pot(i, j)%pot%set(1)%gal21%a32 + pot_par(nk, 18) = potparm%pot(i, j)%pot%set(1)%gal21%a33 + pot_par(nk, 19) = potparm%pot(i, j)%pot%set(1)%gal21%a41 + pot_par(nk, 20) = potparm%pot(i, j)%pot%set(1)%gal21%a42 + pot_par(nk, 21) = potparm%pot(i, j)%pot%set(1)%gal21%a43 + pot_par(nk, 22) = potparm%pot(i, j)%pot%set(1)%gal21%AO1 + pot_par(nk, 23) = potparm%pot(i, j)%pot%set(1)%gal21%AO2 + pot_par(nk, 24) = potparm%pot(i, j)%pot%set(1)%gal21%BO1 + pot_par(nk, 25) = potparm%pot(i, j)%pot%set(1)%gal21%BO2 + pot_par(nk, 26) = potparm%pot(i, j)%pot%set(1)%gal21%c + pot_par(nk, 27) = potparm%pot(i, j)%pot%set(1)%gal21%AH1 + pot_par(nk, 28) = potparm%pot(i, j)%pot%set(1)%gal21%AH2 + pot_par(nk, 29) = potparm%pot(i, j)%pot%set(1)%gal21%BH1 + pot_par(nk, 30) = potparm%pot(i, j)%pot%set(1)%gal21%BH2 CASE (nn_type) ! no checks CASE DEFAULT diff --git a/src/pair_potential_types.F b/src/pair_potential_types.F index 809d62100a..a36430f38a 100644 --- a/src/pair_potential_types.F +++ b/src/pair_potential_types.F @@ -48,9 +48,10 @@ MODULE pair_potential_types ftd_type = 12, & siepmann_type = 13, & gal_type = 14, & - quip_type = 15 + quip_type = 15, & + gal21_type = 16 - INTEGER, PUBLIC, PARAMETER, DIMENSION(16) :: list_pot = (/nn_type, & + INTEGER, PUBLIC, PARAMETER, DIMENSION(17) :: list_pot = (/nn_type, & lj_type, & lj_charmm_type, & ft_type, & @@ -65,7 +66,8 @@ MODULE pair_potential_types ftd_type, & siepmann_type, & gal_type, & - quip_type/) + quip_type, & + gal21_type/) ! Shell model INTEGER, PUBLIC, PARAMETER :: nosh_nosh = 0, & nosh_sh = 1, & @@ -99,7 +101,8 @@ MODULE pair_potential_types quip_pot_type, & tersoff_pot_type, & siepmann_pot_type, & - gal_pot_type + gal_pot_type, & + gal21_pot_type PUBLIC :: pair_potential_lj_create PUBLIC :: compare_pot @@ -256,6 +259,47 @@ MODULE pair_potential_types LOGICAL :: express END TYPE gal_pot_type +! ************************************************************************************************** + + TYPE gal21_pot_type + CHARACTER(LEN=2) :: met1 + CHARACTER(LEN=2) :: met2 + REAL(KIND=dp) :: epsilon1 + REAL(KIND=dp) :: epsilon2 + REAL(KIND=dp) :: epsilon3 + REAL(KIND=dp) :: bxy1 + REAL(KIND=dp) :: bxy2 + REAL(KIND=dp) :: bz1 + REAL(KIND=dp) :: bz2 + REAL(KIND=dp) :: r1 + REAL(KIND=dp) :: r2 + REAL(KIND=dp) :: a11 + REAL(KIND=dp) :: a12 + REAL(KIND=dp) :: a13 + REAL(KIND=dp) :: a21 + REAL(KIND=dp) :: a22 + REAL(KIND=dp) :: a23 + REAL(KIND=dp) :: a31 + REAL(KIND=dp) :: a32 + REAL(KIND=dp) :: a33 + REAL(KIND=dp) :: a41 + REAL(KIND=dp) :: a42 + REAL(KIND=dp) :: a43 + REAL(KIND=dp) :: AO1 + REAL(KIND=dp) :: AO2 + REAL(KIND=dp) :: BO1 + REAL(KIND=dp) :: BO2 + REAL(KIND=dp) :: c + REAL(KIND=dp) :: AH1 + REAL(KIND=dp) :: AH2 + REAL(KIND=dp) :: BH1 + REAL(KIND=dp) :: BH2 + REAL(KIND=dp), POINTER, DIMENSION(:) :: gcn + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: n_vectors + REAL(KIND=dp) :: rcutsq + LOGICAL :: express + END TYPE gal21_pot_type + ! ************************************************************************************************** TYPE pot_set_type @@ -273,6 +317,7 @@ MODULE pair_potential_types TYPE(tersoff_pot_type), POINTER :: tersoff TYPE(siepmann_pot_type), POINTER :: siepmann TYPE(gal_pot_type), POINTER :: gal + TYPE(gal21_pot_type), POINTER :: gal21 TYPE(ftd_pot_type), POINTER :: ftd END TYPE pot_set_type @@ -437,6 +482,40 @@ CONTAINS (pot1%set(i)%gal%express .EQV. & pot2%set(i)%gal%express) .AND. & (pot1%set(i)%gal%rcutsq == pot2%set(i)%gal%rcutsq)) mycompare = .TRUE. + CASE (gal21_type) + IF ((pot1%set(i)%gal21%epsilon1 == pot2%set(i)%gal21%epsilon1) .AND. & + (pot1%set(i)%gal21%epsilon2 == pot2%set(i)%gal21%epsilon2) .AND. & + (pot1%set(i)%gal21%epsilon3 == pot2%set(i)%gal21%epsilon3) .AND. & + (pot1%set(i)%gal21%bxy1 == pot2%set(i)%gal21%bxy1) .AND. & + (pot1%set(i)%gal21%bxy2 == pot2%set(i)%gal21%bxy1) .AND. & + (pot1%set(i)%gal21%bz1 == pot2%set(i)%gal21%bz1) .AND. & + (pot1%set(i)%gal21%bz2 == pot2%set(i)%gal21%bz2) .AND. & + (pot1%set(i)%gal21%r1 == pot2%set(i)%gal21%r1) .AND. & + (pot1%set(i)%gal21%r2 == pot2%set(i)%gal21%r2) .AND. & + (pot1%set(i)%gal21%a11 == pot2%set(i)%gal21%a11) .AND. & + (pot1%set(i)%gal21%a12 == pot2%set(i)%gal21%a12) .AND. & + (pot1%set(i)%gal21%a13 == pot2%set(i)%gal21%a13) .AND. & + (pot1%set(i)%gal21%a21 == pot2%set(i)%gal21%a21) .AND. & + (pot1%set(i)%gal21%a22 == pot2%set(i)%gal21%a22) .AND. & + (pot1%set(i)%gal21%a23 == pot2%set(i)%gal21%a23) .AND. & + (pot1%set(i)%gal21%a31 == pot2%set(i)%gal21%a31) .AND. & + (pot1%set(i)%gal21%a32 == pot2%set(i)%gal21%a32) .AND. & + (pot1%set(i)%gal21%a33 == pot2%set(i)%gal21%a33) .AND. & + (pot1%set(i)%gal21%a41 == pot2%set(i)%gal21%a41) .AND. & + (pot1%set(i)%gal21%a42 == pot2%set(i)%gal21%a42) .AND. & + (pot1%set(i)%gal21%a43 == pot2%set(i)%gal21%a43) .AND. & + (pot1%set(i)%gal21%AO1 == pot2%set(i)%gal21%AO1) .AND. & + (pot1%set(i)%gal21%AO2 == pot2%set(i)%gal21%AO2) .AND. & + (pot1%set(i)%gal21%BO1 == pot2%set(i)%gal21%BO1) .AND. & + (pot1%set(i)%gal21%BO2 == pot2%set(i)%gal21%BO2) .AND. & + (pot1%set(i)%gal21%c == pot2%set(i)%gal21%c) .AND. & + (pot1%set(i)%gal21%AH1 == pot2%set(i)%gal21%AH1) .AND. & + (pot1%set(i)%gal21%AH2 == pot2%set(i)%gal21%AH2) .AND. & + (pot1%set(i)%gal21%BH1 == pot2%set(i)%gal21%BH1) .AND. & + (pot1%set(i)%gal21%BH2 == pot2%set(i)%gal21%BH2) .AND. & + (pot1%set(i)%gal21%express .EQV. & + pot2%set(i)%gal21%express) .AND. & + (pot1%set(i)%gal21%rcutsq == pot2%set(i)%gal21%rcutsq)) mycompare = .TRUE. END SELECT mycompare = mycompare .AND. & @@ -484,6 +563,7 @@ CONTAINS potparm%set(i)%tersoff, & potparm%set(i)%siepmann, & potparm%set(i)%gal, & + potparm%set(i)%gal21, & potparm%set(i)%ftd) END DO CALL pair_potential_single_clean(potparm) @@ -529,6 +609,7 @@ CONTAINS CALL pair_potential_tersoff_clean(potparm%set(i)%tersoff) CALL pair_potential_siepmann_clean(potparm%set(i)%siepmann) CALL pair_potential_gal_clean(potparm%set(i)%gal) + CALL pair_potential_gal21_clean(potparm%set(i)%gal21) END DO END SUBROUTINE pair_potential_single_clean @@ -582,6 +663,7 @@ CONTAINS CALL pair_potential_tersoff_copy(potparm_source%set(i)%tersoff, potparm_dest%set(i)%tersoff) CALL pair_potential_siepmann_copy(potparm_source%set(i)%siepmann, potparm_dest%set(i)%siepmann) CALL pair_potential_gal_copy(potparm_source%set(i)%gal, potparm_dest%set(i)%gal) + CALL pair_potential_gal21_copy(potparm_source%set(i)%gal21, potparm_dest%set(i)%gal21) END DO END SUBROUTINE pair_potential_single_copy @@ -627,6 +709,7 @@ CONTAINS (ASSOCIATED(potparm_dest%set(1)%gp)) .OR. & (ASSOCIATED(potparm_dest%set(1)%tersoff)) .OR. & (ASSOCIATED(potparm_dest%set(1)%siepmann)) .OR. & + (ASSOCIATED(potparm_dest%set(1)%gal)) .OR. & (ASSOCIATED(potparm_dest%set(1)%gal)) IF (.NOT. check) THEN allocate_new = .TRUE. @@ -676,6 +759,7 @@ CONTAINS CALL pair_potential_tersoff_copy(potparm_tmp%set(i)%tersoff, potparm_dest%set(i)%tersoff) CALL pair_potential_siepmann_copy(potparm_tmp%set(i)%siepmann, potparm_dest%set(i)%siepmann) CALL pair_potential_gal_copy(potparm_tmp%set(i)%gal, potparm_dest%set(i)%gal) + CALL pair_potential_gal21_copy(potparm_tmp%set(i)%gal21, potparm_dest%set(i)%gal21) END DO CALL pair_potential_single_release(potparm_tmp) END IF @@ -708,6 +792,7 @@ CONTAINS CALL pair_potential_tersoff_copy(potparm_source%set(j)%tersoff, potparm_dest%set(i)%tersoff) CALL pair_potential_siepmann_copy(potparm_source%set(j)%siepmann, potparm_dest%set(i)%siepmann) CALL pair_potential_gal_copy(potparm_source%set(j)%gal, potparm_dest%set(i)%gal) + CALL pair_potential_gal21_copy(potparm_source%set(j)%gal21, potparm_dest%set(i)%gal21) END DO END SUBROUTINE pair_potential_single_add @@ -739,6 +824,7 @@ CONTAINS CALL pair_potential_tersoff_release(potparm%set(i)%tersoff) CALL pair_potential_siepmann_release(potparm%set(i)%siepmann) CALL pair_potential_gal_release(potparm%set(i)%gal) + CALL pair_potential_gal21_release(potparm%set(i)%gal21) END DO DEALLOCATE (potparm%type) DEALLOCATE (potparm%set) @@ -913,20 +999,22 @@ CONTAINS !> \param tersoff ... !> \param siepmann ... !> \param gal ... +!> \param gal21 ... !> \author Teodoro Laino [Teo] 11.2005 ! ************************************************************************************************** SUBROUTINE pair_potential_reallocate(p, lb1_new, ub1_new, lj, lj_charmm, williams, & goodwin, eam, quip, bmhft, bmhftd, ipbv, buck4r, buckmo, gp, tersoff, & - siepmann, gal) + siepmann, gal, gal21) TYPE(pair_potential_p_type), POINTER :: p INTEGER, INTENT(IN) :: lb1_new, ub1_new LOGICAL, INTENT(IN), OPTIONAL :: lj, lj_charmm, williams, goodwin, eam, & quip, bmhft, bmhftd, ipbv, buck4r, & - buckmo, gp, tersoff, siepmann, gal + buckmo, gp, tersoff, siepmann, gal, & + gal21 INTEGER :: i, ipot, lb1_old, std_dim, ub1_old - LOGICAL :: check, lbmhft, lbmhftd, lbuck4r, lbuckmo, leam, lgal, lgoodwin, lgp, lipbv, llj, & - llj_charmm, lquip, lsiepmann, ltersoff, lwilliams + LOGICAL :: check, lbmhft, lbmhftd, lbuck4r, lbuckmo, leam, lgal, lgal21, lgoodwin, lgp, & + lipbv, llj, llj_charmm, lquip, lsiepmann, ltersoff, lwilliams TYPE(pair_potential_p_type), POINTER :: work NULLIFY (work) @@ -946,96 +1034,103 @@ CONTAINS ltersoff = .FALSE.; IF (PRESENT(tersoff)) ltersoff = tersoff lsiepmann = .FALSE.; IF (PRESENT(siepmann)) lsiepmann = siepmann lgal = .FALSE.; IF (PRESENT(gal)) lgal = gal + lgal21 = .FALSE.; IF (PRESENT(gal21)) lgal21 = gal21 IF (llj) THEN ipot = lj_type check = .NOT. (llj_charmm .OR. lwilliams .OR. lgoodwin .OR. leam .OR. lquip .OR. lbmhft & .OR. lbmhftd .OR. lipbv .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff & - .OR. lsiepmann .OR. lgal) + .OR. lsiepmann .OR. lgal .OR. lgal21) CPASSERT(check) END IF IF (llj_charmm) THEN ipot = lj_charmm_type check = .NOT. (llj .OR. lwilliams .OR. lgoodwin .OR. leam .OR. lquip .OR. lbmhft .OR. lbmhftd .OR. lipbv & - .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal) + .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal .OR. lgal21) CPASSERT(check) END IF IF (lwilliams) THEN ipot = wl_type check = .NOT. (llj .OR. llj_charmm .OR. lgoodwin .OR. leam .OR. lquip .OR. lbmhft .OR. lbmhftd .OR. lipbv & - .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal) + .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal .OR. lgal21) CPASSERT(check) END IF IF (lgoodwin) THEN ipot = gw_type check = .NOT. (llj .OR. llj_charmm .OR. lwilliams .OR. leam .OR. lquip .OR. lbmhft .OR. lbmhftd .OR. lipbv & - .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal) + .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal .OR. lgal21) CPASSERT(check) END IF IF (leam) THEN ipot = ea_type check = .NOT. (llj .OR. llj_charmm .OR. lgoodwin .OR. lwilliams .OR. lquip .OR. lbmhft .OR. lbmhftd .OR. lipbv & - .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal) + .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal .OR. lgal21) CPASSERT(check) END IF IF (lquip) THEN ipot = quip_type check = .NOT. (llj .OR. llj_charmm .OR. lgoodwin .OR. lwilliams .OR. leam .OR. lbmhft .OR. lbmhftd .OR. lipbv & - .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal) + .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal .OR. lgal21) CPASSERT(check) END IF IF (lbmhft) THEN ipot = ft_type check = .NOT. (llj .OR. llj_charmm .OR. lgoodwin .OR. leam .OR. lquip .OR. lwilliams .OR. lbmhftd .OR. lipbv & - .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal) + .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal .OR. lgal21) CPASSERT(check) END IF IF (lbmhftd) THEN ipot = ftd_type check = .NOT. (llj .OR. llj_charmm .OR. lgoodwin .OR. leam .OR. lquip .OR. lwilliams .OR. lbmhft .OR. lipbv & - .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal) + .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal .OR. lgal21) CPASSERT(check) END IF IF (lipbv) THEN ipot = ip_type check = .NOT. (llj .OR. llj_charmm .OR. lgoodwin .OR. leam .OR. lquip .OR. lwilliams .OR. lbmhft .OR. lbmhftd & - .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal) + .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal .OR. lgal21) CPASSERT(check) END IF IF (lbuck4r) THEN ipot = b4_type check = .NOT. (llj .OR. llj_charmm .OR. lgoodwin .OR. leam .OR. lquip .OR. lwilliams .OR. lbmhft .OR. lbmhftd & - .OR. lipbv .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal) + .OR. lipbv .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal .OR. lgal21) CPASSERT(check) END IF IF (lbuckmo) THEN ipot = bm_type check = .NOT. (llj .OR. llj_charmm .OR. lgoodwin .OR. leam .OR. lquip .OR. lwilliams .OR. lbmhft .OR. lbmhftd & - .OR. lipbv .OR. lbuck4r .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal) + .OR. lipbv .OR. lbuck4r .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal .OR. lgal21) CPASSERT(check) END IF IF (ltersoff) THEN ipot = tersoff_type check = .NOT. (llj .OR. llj_charmm .OR. lgoodwin .OR. leam .OR. lquip .OR. lwilliams .OR. lbmhft .OR. lbmhftd & - .OR. lipbv .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. lsiepmann .OR. lgal) + .OR. lipbv .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. lsiepmann .OR. lgal .OR. lgal21) CPASSERT(check) END IF IF (lsiepmann) THEN ipot = siepmann_type check = .NOT. (llj .OR. llj_charmm .OR. lgoodwin .OR. leam .OR. lquip .OR. lwilliams .OR. lbmhft .OR. lbmhftd & - .OR. lipbv .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lgal) + .OR. lipbv .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lgal .OR. lgal21) CPASSERT(check) END IF IF (lgal) THEN ipot = gal_type check = .NOT. (llj .OR. llj_charmm .OR. lgoodwin .OR. leam .OR. lquip .OR. lwilliams .OR. lbmhft .OR. lbmhftd & - .OR. lipbv .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann) + .OR. lipbv .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal21) + CPASSERT(check) + END IF + IF (lgal21) THEN + ipot = gal21_type + check = .NOT. (llj .OR. llj_charmm .OR. lgoodwin .OR. leam .OR. lquip .OR. lwilliams .OR. lbmhft .OR. lbmhftd & + .OR. lipbv .OR. lbuck4r .OR. lbuckmo .OR. lgp .OR. ltersoff .OR. lsiepmann .OR. lgal) CPASSERT(check) END IF IF (lgp) THEN ipot = gp_type check = .NOT. (llj .OR. llj_charmm .OR. lgoodwin .OR. leam .OR. lwilliams .OR. lbmhft .OR. lbmhftd & - .OR. lipbv .OR. lbuck4r .OR. lbuckmo .OR. ltersoff .OR. lsiepmann .OR. lgal) + .OR. lipbv .OR. lbuck4r .OR. lbuckmo .OR. ltersoff .OR. lsiepmann .OR. lgal .OR. lgal21) CPASSERT(check) END IF @@ -1095,6 +1190,8 @@ CONTAINS CALL pair_potential_siepmann_create(p%pot(i)%pot%set(std_dim)%siepmann) CASE (gal_type) CALL pair_potential_gal_create(p%pot(i)%pot%set(std_dim)%gal) + CASE (gal21_type) + CALL pair_potential_gal21_create(p%pot(i)%pot%set(std_dim)%gal21) END SELECT NULLIFY (p%pot(i)%pot%spl_f) NULLIFY (p%pot(i)%pot%pair_spline_data) @@ -2013,7 +2110,7 @@ CONTAINS END SUBROUTINE pair_potential_siepmann_release ! ************************************************************************************************** -!> \brief Creates the GAL17 potential type +!> \brief Creates the GAL19 potential type !> (??) !> \param gal ... ! ************************************************************************************************** @@ -2059,7 +2156,7 @@ CONTAINS END SUBROUTINE pair_potential_gal_copy ! ************************************************************************************************** -!> \brief Creates the GAL17 potential type +!> \brief Creates the GAL19 potential type !> (??) !> \param gal ... ! ************************************************************************************************** @@ -2085,7 +2182,7 @@ CONTAINS END SUBROUTINE pair_potential_gal_clean ! ************************************************************************************************** -!> \brief Destroys the GAL17 potential +!> \brief Destroys the GAL19 potential !> (??) !> \param gal ... ! ************************************************************************************************** @@ -2099,5 +2196,128 @@ CONTAINS NULLIFY (gal) END SUBROUTINE pair_potential_gal_release +! ************************************************************************************************** +!> \brief Creates the GAL21 potential type +!> (??) +!> \param gal21 ... +! ************************************************************************************************** + SUBROUTINE pair_potential_gal21_create(gal21) + TYPE(gal21_pot_type), POINTER :: gal21 + + CPASSERT(.NOT. ASSOCIATED(gal21)) + ALLOCATE (gal21) + CALL pair_potential_gal21_clean(gal21) + END SUBROUTINE pair_potential_gal21_create + +! ************************************************************************************************** +!> \brief Copy two GAL21 potential type +!> (??) +!> \param gal21_source ... +!> \param gal21_dest ... +! ************************************************************************************************** + SUBROUTINE pair_potential_gal21_copy(gal21_source, gal21_dest) + TYPE(gal21_pot_type), POINTER :: gal21_source, gal21_dest + + IF (.NOT. ASSOCIATED(gal21_source)) RETURN + IF (ASSOCIATED(gal21_dest)) CALL pair_potential_gal21_release(gal21_dest) + CALL pair_potential_gal21_create(gal21_dest) + gal21_dest%met1 = gal21_source%met1 + gal21_dest%met2 = gal21_source%met2 + gal21_dest%epsilon1 = gal21_source%epsilon1 + gal21_dest%epsilon2 = gal21_source%epsilon2 + gal21_dest%epsilon3 = gal21_source%epsilon3 + gal21_dest%bxy1 = gal21_source%bxy1 + gal21_dest%bxy2 = gal21_source%bxy2 + gal21_dest%bz1 = gal21_source%bz1 + gal21_dest%bz2 = gal21_source%bz2 + gal21_dest%r1 = gal21_source%r1 + gal21_dest%r2 = gal21_source%r2 + gal21_dest%a11 = gal21_source%a11 + gal21_dest%a12 = gal21_source%a12 + gal21_dest%a13 = gal21_source%a13 + gal21_dest%a21 = gal21_source%a21 + gal21_dest%a22 = gal21_source%a22 + gal21_dest%a23 = gal21_source%a23 + gal21_dest%a31 = gal21_source%a31 + gal21_dest%a32 = gal21_source%a32 + gal21_dest%a33 = gal21_source%a33 + gal21_dest%a41 = gal21_source%a41 + gal21_dest%a42 = gal21_source%a42 + gal21_dest%a43 = gal21_source%a43 + gal21_dest%AO1 = gal21_source%AO1 + gal21_dest%AO2 = gal21_source%AO2 + gal21_dest%BO1 = gal21_source%BO1 + gal21_dest%BO2 = gal21_source%BO2 + gal21_dest%c = gal21_source%c + gal21_dest%AH1 = gal21_source%AH1 + gal21_dest%AH2 = gal21_source%AH2 + gal21_dest%BH1 = gal21_source%BH1 + gal21_dest%BH2 = gal21_source%BH2 + ALLOCATE (gal21_dest%gcn(SIZE(gal21_source%gcn))) + gal21_dest%gcn = gal21_source%gcn + gal21_dest%express = gal21_source%express + gal21_dest%rcutsq = gal21_source%rcutsq + + END SUBROUTINE pair_potential_gal21_copy + +! ************************************************************************************************** +!> \brief Creates the GAL21 potential type +!> (??) +!> \param gal21 ... +! ************************************************************************************************** + SUBROUTINE pair_potential_gal21_clean(gal21) + TYPE(gal21_pot_type), POINTER :: gal21 + + IF (.NOT. ASSOCIATED(gal21)) RETURN + gal21%epsilon1 = 0.0_dp + gal21%epsilon2 = 0.0_dp + gal21%epsilon3 = 0.0_dp + gal21%bxy1 = 0.0_dp + gal21%bxy2 = 0.0_dp + gal21%bz1 = 0.0_dp + gal21%bz2 = 0.0_dp + gal21%r1 = 0.0_dp + gal21%r2 = 0.0_dp + gal21%a11 = 0.0_dp + gal21%a12 = 0.0_dp + gal21%a13 = 0.0_dp + gal21%a21 = 0.0_dp + gal21%a22 = 0.0_dp + gal21%a23 = 0.0_dp + gal21%a31 = 0.0_dp + gal21%a32 = 0.0_dp + gal21%a33 = 0.0_dp + gal21%a41 = 0.0_dp + gal21%a42 = 0.0_dp + gal21%a43 = 0.0_dp + gal21%AO1 = 0.0_dp + gal21%AO2 = 0.0_dp + gal21%BO1 = 0.0_dp + gal21%BO2 = 0.0_dp + gal21%c = 0.0_dp + gal21%AH1 = 0.0_dp + gal21%AH2 = 0.0_dp + gal21%BH1 = 0.0_dp + gal21%BH2 = 0.0_dp + gal21%rcutsq = 0.0_dp + gal21%express = .FALSE. + + END SUBROUTINE pair_potential_gal21_clean + +! ************************************************************************************************** +!> \brief Destroys the GAL21 potential +!> (??) +!> \param gal21 ... +! ************************************************************************************************** + SUBROUTINE pair_potential_gal21_release(gal21) + TYPE(gal21_pot_type), POINTER :: gal21 + + IF (ASSOCIATED(gal21)) THEN + DEALLOCATE (gal21%gcn) + DEALLOCATE (gal21) + END IF + NULLIFY (gal21) + END SUBROUTINE pair_potential_gal21_release + END MODULE pair_potential_types diff --git a/tests/Fist/regtest-gal21/TEST_FILES b/tests/Fist/regtest-gal21/TEST_FILES new file mode 100644 index 0000000000..9bd9f4f35d --- /dev/null +++ b/tests/Fist/regtest-gal21/TEST_FILES @@ -0,0 +1,7 @@ +# runs are executed in the same order as in this file +# the second field tells which test should be run in order to compare with the last available output +# e.g. 0 means do not compare anything, running is enough +# 1 compares the last total energy in the file +# for details see cp2k/tools/do_regtest +test_gal21.inp 11 8e-12 -2.729695004946153 +#EOF diff --git a/tests/Fist/regtest-gal21/coord.inc b/tests/Fist/regtest-gal21/coord.inc new file mode 100644 index 0000000000..73e3608ef2 --- /dev/null +++ b/tests/Fist/regtest-gal21/coord.inc @@ -0,0 +1,870 @@ +Pt -7.791 -9.089 -3.372 Pt +Pt -4.979 -9.089 -3.372 Pt +Pt -2.167 -9.089 -3.372 Pt +Pt 0.645 -9.089 -3.372 Pt +Pt 3.457 -9.089 -3.372 Pt +Pt 6.269 -9.089 -3.372 Pt +Pt -6.385 -6.654 -3.372 Pt +Pt -3.573 -6.654 -3.372 Pt +Pt -0.761 -6.654 -3.372 Pt +Pt 2.051 -6.654 -3.372 Pt +Pt 4.863 -6.654 -3.372 Pt +Pt 7.675 -6.654 -3.372 Pt +Pt -7.791 -4.218 -3.372 Pt +Pt -4.979 -4.218 -3.372 Pt +Pt -2.167 -4.218 -3.372 Pt +Pt 0.645 -4.218 -3.372 Pt +Pt 3.457 -4.218 -3.372 Pt +Pt 6.269 -4.218 -3.372 Pt +Pt -6.385 -1.783 -3.372 Pt +Pt -3.573 -1.783 -3.372 Pt +Pt -0.761 -1.783 -3.372 Pt +Pt 2.051 -1.783 -3.372 Pt +Pt 4.863 -1.783 -3.372 Pt +Pt 7.675 -1.783 -3.372 Pt +Pt -7.791 0.652 -3.372 Pt +Pt -4.979 0.652 -3.372 Pt +Pt -2.167 0.652 -3.372 Pt +Pt 0.645 0.652 -3.372 Pt +Pt 3.457 0.652 -3.372 Pt +Pt 6.269 0.652 -3.372 Pt +Pt -6.385 3.087 -3.372 Pt +Pt -3.573 3.087 -3.372 Pt +Pt -0.761 3.087 -3.372 Pt +Pt 2.051 3.087 -3.372 Pt +Pt 4.863 3.087 -3.372 Pt +Pt 7.675 3.087 -3.372 Pt +Pt -7.791 5.523 -3.372 Pt +Pt -4.979 5.523 -3.372 Pt +Pt -2.167 5.523 -3.372 Pt +Pt 0.645 5.523 -3.372 Pt +Pt 3.457 5.523 -3.372 Pt +Pt 6.269 5.523 -3.372 Pt +Pt -6.385 7.958 -3.372 Pt +Pt -3.573 7.958 -3.372 Pt +Pt -0.761 7.958 -3.372 Pt +Pt 2.051 7.958 -3.372 Pt +Pt 4.863 7.958 -3.372 Pt +Pt 7.675 7.958 -3.372 Pt +Pt -6.385 -8.277 -1.076 Pt +Pt -3.573 -8.277 -1.076 Pt +Pt -0.761 -8.277 -1.076 Pt +Pt 2.051 -8.277 -1.076 Pt +Pt 4.863 -8.277 -1.076 Pt +Pt 7.675 -8.277 -1.076 Pt +Pt -7.791 -5.842 -1.076 Pt +Pt -4.979 -5.842 -1.076 Pt +Pt -2.167 -5.842 -1.076 Pt +Pt 0.645 -5.842 -1.076 Pt +Pt 3.457 -5.842 -1.076 Pt +Pt 6.269 -5.842 -1.076 Pt +Pt -6.385 -3.407 -1.076 Pt +Pt -3.573 -3.407 -1.076 Pt +Pt -0.761 -3.407 -1.076 Pt +Pt 2.051 -3.407 -1.076 Pt +Pt 4.863 -3.407 -1.076 Pt +Pt 7.675 -3.407 -1.076 Pt +Pt -7.791 -0.971 -1.076 Pt +Pt -4.979 -0.971 -1.076 Pt +Pt -2.167 -0.971 -1.076 Pt +Pt 0.645 -0.971 -1.076 Pt +Pt 3.457 -0.971 -1.076 Pt +Pt 6.269 -0.971 -1.076 Pt +Pt -6.385 1.464 -1.076 Pt +Pt -3.573 1.464 -1.076 Pt +Pt -0.761 1.464 -1.076 Pt +Pt 2.051 1.464 -1.076 Pt +Pt 4.863 1.464 -1.076 Pt +Pt 7.675 1.464 -1.076 Pt +Pt -7.791 3.899 -1.076 Pt +Pt -4.979 3.899 -1.076 Pt +Pt -2.167 3.899 -1.076 Pt +Pt 0.645 3.899 -1.076 Pt +Pt 3.457 3.899 -1.076 Pt +Pt 6.269 3.899 -1.076 Pt +Pt -6.385 6.334 -1.076 Pt +Pt -3.573 6.334 -1.076 Pt +Pt -0.761 6.334 -1.076 Pt +Pt 2.051 6.334 -1.076 Pt +Pt 4.863 6.334 -1.076 Pt +Pt 7.675 6.334 -1.076 Pt +Pt -7.791 8.770 -1.076 Pt +Pt -4.979 8.770 -1.076 Pt +Pt -2.167 8.770 -1.076 Pt +Pt 0.645 8.770 -1.076 Pt +Pt 3.457 8.770 -1.076 Pt +Pt 6.269 8.770 -1.076 Pt +Pt -7.791 -7.465 1.220 Pt +Pt -4.979 -7.465 1.220 Pt +Pt -2.167 -7.465 1.220 Pt +Pt 0.645 -7.465 1.220 Pt +Pt 3.457 -7.465 1.220 Pt +Pt 6.269 -7.465 1.220 Pt +Pt -6.385 -5.030 1.220 Pt +Pt -3.573 -5.030 1.220 Pt +Pt -0.761 -5.030 1.220 Pt +Pt 2.051 -5.030 1.220 Pt +Pt 4.863 -5.030 1.220 Pt +Pt 7.675 -5.030 1.220 Pt +Pt -7.791 -2.595 1.220 Pt +Pt -4.979 -2.595 1.220 Pt +Pt -2.167 -2.595 1.220 Pt +Pt 0.645 -2.595 1.220 Pt +Pt 3.457 -2.595 1.220 Pt +Pt 6.269 -2.595 1.220 Pt +Pt -6.385 -0.160 1.220 Pt +Pt -3.573 -0.160 1.220 Pt +Pt -0.761 -0.160 1.220 Pt +Pt 2.051 -0.160 1.220 Pt +Pt 4.863 -0.160 1.220 Pt +Pt 7.675 -0.160 1.220 Pt +Pt -7.791 2.276 1.220 Pt +Pt -4.979 2.276 1.220 Pt +Pt -2.167 2.276 1.220 Pt +Pt 0.645 2.276 1.220 Pt +Pt 3.457 2.276 1.220 Pt +Pt 6.269 2.276 1.220 Pt +Pt -6.385 4.711 1.220 Pt +Pt -3.573 4.711 1.220 Pt +Pt -0.761 4.711 1.220 Pt +Pt 2.051 4.711 1.220 Pt +Pt 4.863 4.711 1.220 Pt +Pt 7.675 4.711 1.220 Pt +Pt -7.791 7.146 1.220 Pt +Pt -4.979 7.146 1.220 Pt +Pt -2.167 7.146 1.220 Pt +Pt 0.645 7.146 1.220 Pt +Pt 3.457 7.146 1.220 Pt +Pt 6.269 7.146 1.220 Pt +Pt -6.385 9.582 1.220 Pt +Pt -3.573 9.582 1.220 Pt +Pt -0.761 9.582 1.220 Pt +Pt 2.051 9.582 1.220 Pt +Pt 4.863 9.582 1.220 Pt +Pt 7.675 9.582 1.220 Pt +Pt -7.791 -9.089 3.516 Pt +Pt -4.979 -9.089 3.516 Pt +Pt -2.167 -9.089 3.516 Pt +Pt 0.645 -9.089 3.516 Pt +Pt 3.457 -9.089 3.516 Pt +Pt 6.269 -9.089 3.516 Pt +Pt -6.385 -6.654 3.516 Pt +Pt -3.573 -6.654 3.516 Pt +Pt -0.761 -6.654 3.516 Pt +Pt 2.051 -6.654 3.516 Pt +Pt 4.863 -6.654 3.516 Pt +Pt 7.675 -6.654 3.516 Pt +Pt -7.791 -4.218 3.516 Pt +Pt -4.979 -4.218 3.516 Pt +Pt -2.167 -4.218 3.516 Pt +Pt 0.645 -4.218 3.516 Pt +Pt 3.457 -4.218 3.516 Pt +Pt 6.269 -4.218 3.516 Pt +Pt -6.385 -1.783 3.516 Pt +Pt -3.573 -1.783 3.516 Pt +Pt -0.761 -1.783 3.516 Pt +Pt 2.051 -1.783 3.516 Pt +Pt 4.863 -1.783 3.516 Pt +Pt 7.675 -1.783 3.516 Pt +Pt -7.791 0.652 3.516 Pt +Pt -4.979 0.652 3.516 Pt +Pt -2.167 0.652 3.516 Pt +Pt 0.645 0.652 3.516 Pt +Pt 3.457 0.652 3.516 Pt +Pt 6.269 0.652 3.516 Pt +Pt -6.385 3.087 3.516 Pt +Pt -3.573 3.087 3.516 Pt +Pt -0.761 3.087 3.516 Pt +Pt 2.051 3.087 3.516 Pt +Pt 4.863 3.087 3.516 Pt +Pt 7.675 3.087 3.516 Pt +Pt -7.791 5.523 3.516 Pt +Pt -4.979 5.523 3.516 Pt +Pt -2.167 5.523 3.516 Pt +Pt 0.645 5.523 3.516 Pt +Pt 3.457 5.523 3.516 Pt +Pt 6.269 5.523 3.516 Pt +Pt -6.385 7.958 3.516 Pt +Pt -3.573 7.958 3.516 Pt +Pt -0.761 7.958 3.516 Pt +Pt 2.051 7.958 3.516 Pt +Pt 4.863 7.958 3.516 Pt +Pt 7.675 7.958 3.516 Pt +O 2.615 5.885 23.170 WAT +H 1.727 6.227 23.068 WAT +H 2.899 5.685 22.278 WAT +O 0.579 8.486 29.054 WAT +H 1.150 7.718 29.054 WAT +H 0.579 8.785 28.144 WAT +O 4.791 4.469 24.203 WAT +H 4.898 5.369 24.510 WAT +H 3.881 4.420 23.913 WAT +O 6.632 7.437 28.536 WAT +H 7.140 7.935 29.176 WAT +H 7.163 6.661 28.360 WAT +O 8.280 1.276 27.134 WAT +H 8.588 0.721 26.417 WAT +H 7.350 1.411 26.951 WAT +O 6.274 9.868 26.652 WAT +H 5.870 9.155 26.156 WAT +H 6.360 9.523 27.541 WAT +O 2.974 8.762 24.152 WAT +H 2.854 9.488 23.541 WAT +H 2.664 7.992 23.675 WAT +O 5.441 2.106 27.535 WAT +H 4.997 2.887 27.203 WAT +H 5.999 2.429 28.242 WAT +O 8.557 1.727 22.415 WAT +H 7.932 1.723 23.141 WAT +H 8.470 0.860 22.019 WAT +O 7.478 2.650 29.394 WAT +H 7.775 3.557 29.313 WAT +H 8.027 2.162 28.780 WAT +O 5.232 7.229 25.263 WAT +H 5.006 7.013 26.168 WAT +H 4.451 7.662 24.917 WAT +O 0.105 2.462 28.115 WAT +H -0.184 3.032 28.828 WAT +H 0.849 1.982 28.478 WAT +O 0.292 6.404 19.326 WAT +H 0.209 7.347 19.467 WAT +H -0.060 6.261 18.447 WAT +O 8.092 5.138 28.393 WAT +H 8.095 4.697 27.543 WAT +H 8.987 5.459 28.500 WAT +O 3.942 5.020 27.352 WAT +H 3.602 5.498 28.109 WAT +H 3.180 4.567 26.992 WAT +O 8.154 6.139 21.383 WAT +H 8.275 6.098 22.332 WAT +H 7.768 7.001 21.228 WAT +O 7.166 8.449 20.782 WAT +H 6.597 8.001 20.155 WAT +H 7.875 8.808 20.248 WAT +O 1.256 3.942 26.268 WAT +H 0.971 3.578 27.106 WAT +H 1.194 3.210 25.655 WAT +O 0.749 9.078 20.228 WAT +H 1.619 9.361 20.512 WAT +H 0.259 9.892 20.113 WAT +O 3.814 0.809 24.821 WAT +H 3.686 0.801 25.770 WAT +H 2.990 1.148 24.470 WAT +O 7.642 6.295 24.257 WAT +H 6.934 6.764 24.698 WAT +H 7.926 5.637 24.891 WAT +O 5.846 4.603 21.632 WAT +H 6.695 5.040 21.566 WAT +H 5.479 4.915 22.459 WAT +O 3.439 8.289 20.789 WAT +H 3.520 7.342 20.900 WAT +H 4.331 8.620 20.895 WAT +O 5.672 6.926 19.025 WAT +H 5.471 6.063 19.387 WAT +H 5.502 6.840 18.087 WAT +O 4.471 4.027 19.013 WAT +H 4.773 4.033 18.105 WAT +H 5.210 3.676 19.510 WAT +O 7.689 2.970 19.999 WAT +H 7.341 2.292 19.419 WAT +H 7.891 2.509 20.813 WAT +O 1.030 9.499 26.395 WAT +H 0.375 9.230 25.751 WAT +H 1.868 9.389 25.945 WAT +O 0.244 6.709 22.338 WAT +H -0.660 6.972 22.163 WAT +H 0.728 6.980 21.558 WAT +O 2.786 5.765 20.247 WAT +H 3.495 5.344 19.762 WAT +H 2.070 5.833 19.615 WAT +O 2.434 1.819 19.287 WAT +H 3.157 2.360 19.606 WAT +H 2.706 1.545 18.412 WAT +O 2.123 1.070 29.190 WAT +H 2.771 0.519 28.749 WAT +H 2.543 1.324 30.011 WAT +O 1.858 1.858 22.881 WAT +H 2.616 1.971 22.307 WAT +H 1.396 2.695 22.835 WAT +O 4.231 2.385 21.421 WAT +H 4.760 1.598 21.292 WAT +H 4.864 3.064 21.655 WAT +O 6.425 2.122 24.476 WAT +H 5.795 1.421 24.643 WAT +H 5.885 2.883 24.259 WAT +O 0.579 8.486 10.279 WAT +H 1.150 7.718 10.280 WAT +H 0.579 8.785 9.370 WAT +O 1.821 3.894 12.066 WAT +H 1.641 2.997 12.346 WAT +H 1.419 4.442 12.741 WAT +O 2.440 6.391 10.674 WAT +H 3.121 6.654 11.293 WAT +H 2.124 5.551 11.008 WAT +O 6.632 7.437 9.762 WAT +H 7.140 7.935 10.402 WAT +H 7.163 6.661 9.585 WAT +O 4.916 2.077 14.028 WAT +H 4.918 2.594 13.222 WAT +H 4.086 1.600 14.005 WAT +O 4.597 3.788 12.035 WAT +H 4.939 4.632 11.740 WAT +H 3.659 3.936 12.151 WAT +O 8.280 1.276 8.359 WAT +H 8.588 0.721 7.643 WAT +H 7.350 1.411 8.177 WAT +O 6.274 9.868 7.877 WAT +H 5.870 9.155 7.382 WAT +H 6.360 9.523 8.766 WAT +O 4.987 6.642 11.657 WAT +H 5.538 6.752 10.883 WAT +H 5.501 7.016 12.373 WAT +O 5.441 2.106 8.760 WAT +H 4.997 2.887 8.429 WAT +H 5.999 2.429 9.468 WAT +O 7.478 2.650 10.619 WAT +H 7.775 3.557 10.539 WAT +H 8.027 2.162 10.005 WAT +O 5.232 7.329 6.489 WAT +H 5.006 7.113 7.393 WAT +H 4.451 7.762 6.143 WAT +O 7.672 8.128 17.155 WAT +H 7.576 7.187 17.005 WAT +H 8.515 8.349 16.759 WAT +O 2.803 4.539 15.308 WAT +H 2.409 5.386 15.515 WAT +H 2.115 3.902 15.499 WAT +O 0.105 2.462 9.341 WAT +H -0.184 3.032 10.053 WAT +H 0.849 1.982 9.704 WAT +O 8.092 5.138 9.618 WAT +H 8.095 4.697 8.769 WAT +H 8.987 5.459 9.725 WAT +O 3.942 5.020 8.578 WAT +H 3.602 5.498 9.334 WAT +H 3.180 4.567 8.218 WAT +O 8.392 9.198 11.365 WAT +H 8.759 8.724 12.111 WAT +H 8.203 10.072 11.706 WAT +O 1.256 3.942 7.494 WAT +H 0.971 3.578 8.332 WAT +H 1.194 3.210 6.880 WAT +O 6.182 8.239 13.683 WAT +H 5.878 9.118 13.908 WAT +H 7.137 8.299 13.719 WAT +O 2.239 9.079 14.830 WAT +H 1.880 9.744 15.418 WAT +H 1.944 8.247 15.202 WAT +O 5.440 9.537 17.714 WAT +H 5.140 9.683 16.817 WAT +H 6.211 8.978 17.622 WAT +O 0.814 2.568 15.508 WAT +H 0.099 2.768 16.113 WAT +H 1.233 1.792 15.879 WAT +O 2.999 8.254 18.145 WAT +H 3.039 7.955 19.053 WAT +H 3.853 8.654 17.986 WAT +O 4.939 6.833 16.112 WAT +H 4.409 6.342 15.485 WAT +H 5.608 7.261 15.578 WAT +O 5.427 3.799 16.192 WAT +H 5.332 3.181 15.467 WAT +H 4.846 4.526 15.965 WAT +O 8.193 2.356 13.225 WAT +H 7.744 2.379 12.379 WAT +H 7.559 1.958 13.822 WAT +O 7.702 5.526 16.487 WAT +H 8.292 5.217 17.175 WAT +H 6.884 5.054 16.643 WAT +O 1.030 9.499 7.621 WAT +H 0.375 9.230 6.977 WAT +H 1.868 9.389 7.171 WAT +O 2.067 9.884 12.161 WAT +H 2.435 9.797 13.041 WAT +H 1.498 9.121 12.062 WAT +O 1.619 6.795 16.067 WAT +H 0.725 6.762 16.406 WAT +H 2.159 6.994 16.832 WAT +O 2.123 1.070 10.415 WAT +H 2.771 0.519 9.975 WAT +H 2.543 1.324 11.237 WAT +O 1.876 1.225 13.158 WAT +H 1.528 1.682 13.924 WAT +H 1.545 0.331 13.240 WAT +O 2.548 0.753 16.433 WAT +H 3.097 0.572 15.670 WAT +H 2.614 -0.038 16.969 WAT +O 8.400 3.972 7.097 WAT +H 7.587 3.598 6.756 WAT +H 8.904 3.217 7.399 WAT +O 6.550 1.401 18.137 WAT +H 6.269 1.904 17.373 WAT +H 6.932 0.603 17.771 WAT +O 4.186 -3.815 29.581 WAT +H 4.489 -4.705 29.402 WAT +H 4.931 -3.258 29.353 WAT +O 0.671 -5.900 21.675 WAT +H 0.294 -5.411 20.944 WAT +H 0.928 -5.227 22.306 WAT +O 6.742 -5.220 28.816 WAT +H 6.445 -5.825 28.136 WAT +H 7.692 -5.182 28.705 WAT +O 1.858 -5.141 26.370 WAT +H 1.450 -4.999 25.515 WAT +H 2.749 -5.425 26.167 WAT +O 6.130 -0.926 24.482 WAT +H 5.954 -1.843 24.269 WAT +H 5.275 -0.567 24.722 WAT +O 4.710 -6.329 23.054 WAT +H 3.937 -6.870 22.894 WAT +H 5.440 -6.849 22.716 WAT +O 4.754 -6.661 26.692 WAT +H 5.272 -7.465 26.638 WAT +H 4.303 -6.724 27.534 WAT +O 7.913 -8.309 24.349 WAT +H 8.369 -7.498 24.575 WAT +H 7.562 -8.627 25.181 WAT +O 1.117 -0.483 20.794 WAT +H 1.476 0.115 20.140 WAT +H 0.976 0.063 21.568 WAT +O 7.133 -3.100 21.612 WAT +H 7.124 -2.367 20.996 WAT +H 6.213 -3.347 21.704 WAT +O 2.734 -7.879 22.356 WAT +H 2.127 -7.139 22.366 WAT +H 3.209 -7.789 21.530 WAT +O 8.761 -0.519 24.969 WAT +H 9.015 -0.557 24.047 WAT +H 7.812 -0.646 24.962 WAT +O 4.745 -7.918 20.105 WAT +H 5.511 -8.030 20.668 WAT +H 4.822 -8.614 19.453 WAT +O 4.241 -0.336 27.892 WAT +H 4.977 -0.884 28.163 WAT +H 4.625 0.530 27.753 WAT +O 1.125 -4.268 23.766 WAT +H 0.199 -4.026 23.757 WAT +H 1.591 -3.441 23.645 WAT +O 7.351 -1.163 19.715 WAT +H 8.144 -0.874 19.263 WAT +H 6.831 -1.593 19.036 WAT +O 7.869 -5.942 19.579 WAT +H 7.764 -6.557 20.305 WAT +H 7.701 -6.466 18.795 WAT +O 6.777 -7.921 22.079 WAT +H 7.198 -8.087 22.923 WAT +H 6.902 -8.731 21.585 WAT +O 3.494 -5.444 19.846 WAT +H 3.931 -4.901 20.502 WAT +H 3.979 -6.269 19.856 WAT +O 0.708 -2.655 27.298 WAT +H 1.039 -3.523 27.066 WAT +H 1.137 -2.449 28.129 WAT +O 5.305 -4.324 24.958 WAT +H 5.004 -4.447 25.859 WAT +H 5.198 -5.183 24.550 WAT +O 6.557 -1.851 28.953 WAT +H 7.074 -2.462 28.428 WAT +H 7.207 -1.348 29.444 WAT +O 1.236 -0.748 25.378 WAT +H 1.068 -1.381 26.076 WAT +H 0.369 -0.542 25.027 WAT +O 2.842 -2.054 22.823 WAT +H 2.317 -1.673 22.119 WAT +H 2.766 -1.429 23.544 WAT +O 5.441 -0.106 21.698 WAT +H 6.173 -0.390 21.150 WAT +H 5.701 -0.345 22.588 WAT +O 4.603 -3.990 21.702 WAT +H 3.974 -3.521 22.251 WAT +H 4.717 -4.835 22.137 WAT +O 4.206 -7.849 10.611 WAT +H 3.375 -8.104 11.011 WAT +H 4.814 -7.778 11.347 WAT +O 4.186 -3.815 10.807 WAT +H 4.489 -4.705 10.628 WAT +H 4.931 -3.258 10.579 WAT +O 6.166 -7.058 12.205 WAT +H 5.833 -6.431 12.847 WAT +H 6.343 -6.532 11.425 WAT +O 5.378 -5.090 14.020 WAT +H 5.995 -4.415 13.739 WAT +H 4.517 -4.681 13.929 WAT +O 6.742 -5.220 10.041 WAT +H 6.445 -5.825 9.362 WAT +H 7.692 -5.182 9.930 WAT +O 5.834 -5.482 16.692 WAT +H 4.937 -5.258 16.939 WAT +H 5.859 -5.358 15.743 WAT +O 4.661 -8.610 15.268 WAT +H 4.522 -7.738 15.637 WAT +H 3.779 -8.953 15.123 WAT +O 6.459 -2.668 17.476 WAT +H 6.431 -2.264 16.609 WAT +H 6.295 -3.597 17.313 WAT +O 0.613 -6.841 11.185 WAT +H 0.691 -6.706 10.240 WAT +H 1.227 -7.549 11.378 WAT +O 1.858 -5.141 7.596 WAT +H 1.450 -4.999 6.741 WAT +H 2.749 -5.425 7.393 WAT +O 0.669 -5.799 13.724 WAT +H 0.832 -6.269 12.906 WAT +H 1.134 -6.307 14.389 WAT +O 2.243 -6.424 15.884 WAT +H 1.858 -7.065 16.482 WAT +H 2.554 -5.720 16.453 WAT +O 4.754 -6.661 7.918 WAT +H 5.272 -7.465 7.863 WAT +H 4.303 -6.724 8.759 WAT +O 1.101 -2.218 16.494 WAT +H 1.210 -3.158 16.354 WAT +H 0.449 -1.952 15.846 WAT +O 8.761 -0.519 6.195 WAT +H 9.015 -0.557 5.273 WAT +H 7.812 -0.646 6.187 WAT +O 4.241 -0.336 9.118 WAT +H 4.977 -0.884 9.388 WAT +H 4.625 0.530 8.978 WAT +O 3.247 -1.151 18.306 WAT +H 4.048 -0.830 18.720 WAT +H 2.732 -1.512 19.028 WAT +O 8.020 -7.257 16.783 WAT +H 7.255 -6.688 16.699 WAT +H 8.202 -7.549 15.890 WAT +O 8.599 -7.523 13.846 WAT +H 9.017 -6.679 13.674 WAT +H 7.680 -7.385 13.617 WAT +O 6.449 -1.217 14.962 WAT +H 5.911 -0.493 14.643 WAT +H 6.703 -1.694 14.171 WAT +O 0.708 -2.655 8.524 WAT +H 1.039 -3.523 8.292 WAT +H 1.137 -2.449 9.354 WAT +O 1.197 -1.788 11.262 WAT +H 0.760 -1.608 12.094 WAT +H 1.765 -1.030 11.120 WAT +O 3.483 -4.377 17.323 WAT +H 3.440 -4.748 18.205 WAT +H 3.614 -3.439 17.463 WAT +O 7.183 -2.821 12.924 WAT +H 6.942 -2.554 12.037 WAT +H 8.123 -2.651 12.979 WAT +O 3.321 -1.431 14.443 WAT +H 2.727 -1.562 15.182 WAT +H 3.314 -2.267 13.977 WAT +O 2.953 -4.065 13.394 WAT +H 2.226 -4.684 13.326 WAT +H 3.014 -3.671 12.524 WAT +O 6.557 -1.851 10.179 WAT +H 7.074 -2.462 9.654 WAT +H 7.207 -1.348 10.670 WAT +O 1.236 -0.748 6.603 WAT +H 1.068 -1.381 7.302 WAT +H 0.369 -0.542 6.253 WAT +O -6.717 8.835 26.627 WAT +H -6.797 7.921 26.356 WAT +H -7.130 9.332 25.921 WAT +O -4.331 3.591 25.674 WAT +H -3.510 3.101 25.717 WAT +H -4.298 4.181 26.428 WAT +O -4.095 5.531 28.447 WAT +H -4.799 5.820 29.028 WAT +H -3.734 6.341 28.087 WAT +O -1.026 4.977 28.838 WAT +H -1.974 5.010 28.970 WAT +H -0.701 5.790 29.223 WAT +O -2.318 1.396 19.695 WAT +H -2.631 1.489 20.595 WAT +H -1.557 1.974 19.644 WAT +O -6.313 7.235 23.862 WAT +H -6.912 7.959 23.678 WAT +H -5.765 7.173 23.080 WAT +O -7.706 5.872 25.892 WAT +H -7.313 6.189 25.079 WAT +H -7.344 4.993 26.007 WAT +O -5.129 7.207 21.286 WAT +H -4.241 7.519 21.111 WAT +H -5.159 6.334 20.894 WAT +O -2.858 8.518 20.754 WAT +H -2.941 8.950 19.904 WAT +H -2.138 8.978 21.187 WAT +O -3.101 4.016 22.841 WAT +H -3.729 4.121 22.127 WAT +H -3.596 4.236 23.630 WAT +O -7.101 3.020 25.696 WAT +H -7.224 3.039 24.747 WAT +H -6.152 3.068 25.816 WAT +O -4.008 9.175 26.545 WAT +H -4.949 9.190 26.720 WAT +H -3.693 8.388 26.991 WAT +O -5.623 1.737 21.554 WAT +H -5.343 0.894 21.910 WAT +H -6.396 1.970 22.068 WAT +O -4.796 4.293 20.617 WAT +H -4.377 4.106 19.777 WAT +H -5.102 3.441 20.927 WAT +O -1.297 9.205 24.315 WAT +H -1.672 8.438 23.883 WAT +H -2.026 9.588 24.802 WAT +O -2.021 7.431 27.615 WAT +H -1.498 6.937 26.984 WAT +H -1.435 7.566 28.360 WAT +O -2.531 6.688 23.555 WAT +H -2.036 6.455 24.339 WAT +H -2.730 5.850 23.138 WAT +O -1.103 5.324 26.088 WAT +H -1.213 5.170 27.027 WAT +H -0.342 4.797 25.846 WAT +O -0.214 3.155 19.329 WAT +H 0.711 3.023 19.121 WAT +H -0.236 3.235 20.283 WAT +O -0.140 3.926 21.955 WAT +H -0.935 3.861 22.484 WAT +H 0.064 4.861 21.943 WAT +O -3.708 2.256 28.404 WAT +H -3.267 2.201 27.557 WAT +H -3.977 3.172 28.476 WAT +O -1.844 1.888 26.415 WAT +H -1.905 0.944 26.268 WAT +H -1.065 1.996 26.961 WAT +O -5.872 0.637 27.109 WAT +H -5.269 1.209 27.584 WAT +H -6.613 1.201 26.887 WAT +O -2.729 1.039 22.384 WAT +H -2.534 1.937 22.650 WAT +H -3.516 0.807 22.877 WAT +O -6.717 8.835 7.852 WAT +H -6.797 7.921 7.582 WAT +H -7.130 9.332 7.146 WAT +O -4.331 3.591 6.900 WAT +H -3.510 3.101 6.943 WAT +H -4.298 4.181 7.653 WAT +O -0.041 5.608 13.464 WAT +H -0.746 5.996 12.946 WAT +H 0.225 6.303 14.067 WAT +O -4.095 5.531 9.673 WAT +H -4.799 5.820 10.253 WAT +H -3.734 6.341 9.313 WAT +O -1.026 4.977 10.064 WAT +H -1.974 5.010 10.195 WAT +H -0.701 5.790 10.449 WAT +O -4.634 7.082 13.016 WAT +H -5.047 7.907 13.270 WAT +H -5.366 6.496 12.823 WAT +O -1.108 6.133 16.898 WAT +H -1.724 6.641 16.369 WAT +H -1.429 5.233 16.840 WAT +O -4.046 6.193 15.494 WAT +H -4.897 5.894 15.815 WAT +H -4.206 6.438 14.582 WAT +O -6.951 7.694 17.593 WAT +H -6.626 6.860 17.254 WAT +H -6.163 8.216 17.741 WAT +O -7.706 5.872 7.118 WAT +H -7.313 6.189 6.305 WAT +H -7.344 4.993 7.232 WAT +O -6.658 3.916 10.625 WAT +H -6.641 4.367 11.469 WAT +H -6.906 3.016 10.839 WAT +O -3.231 1.329 16.464 WAT +H -2.906 0.598 16.990 WAT +H -3.195 1.010 15.562 WAT +O -6.467 5.061 13.258 WAT +H -6.099 4.240 13.586 WAT +H -7.327 5.121 13.673 WAT +O -2.300 3.931 14.324 WAT +H -2.803 4.672 14.664 WAT +H -1.484 4.320 14.008 WAT +O -4.008 9.175 7.770 WAT +H -4.949 9.190 7.945 WAT +H -3.693 8.388 8.217 WAT +O -4.465 8.639 17.880 WAT +H -3.883 7.885 17.792 WAT +H -4.289 9.170 17.103 WAT +O -5.501 9.537 13.806 WAT +H -6.027 10.061 13.203 WAT +H -5.157 10.170 14.437 WAT +O -6.253 5.142 16.721 WAT +H -6.942 4.480 16.656 WAT +H -5.532 4.696 17.165 WAT +O -1.907 8.352 15.174 WAT +H -2.138 8.124 14.274 WAT +H -1.646 9.272 15.129 WAT +O -3.468 3.578 17.952 WAT +H -2.591 3.586 18.336 WAT +H -3.456 2.841 17.341 WAT +O -2.021 7.431 8.841 WAT +H -1.498 6.937 8.210 WAT +H -1.435 7.566 9.586 WAT +O -1.842 7.515 12.363 WAT +H -2.792 7.481 12.477 WAT +H -1.661 8.425 12.128 WAT +O -1.103 5.324 7.314 WAT +H -1.213 5.170 8.252 WAT +H -0.342 4.797 7.072 WAT +O -2.756 1.264 13.772 WAT +H -2.119 1.877 14.138 WAT +H -2.609 1.300 12.827 WAT +O -6.748 2.074 16.319 WAT +H -6.134 2.391 15.657 WAT +H -6.203 1.588 16.939 WAT +O -6.007 0.637 12.646 WAT +H -6.881 0.970 12.443 WAT +H -6.116 0.164 13.471 WAT +O -5.100 2.844 14.030 WAT +H -5.429 2.029 13.650 WAT +H -4.152 2.801 13.907 WAT +O -3.708 2.256 9.630 WAT +H -3.267 2.201 8.782 WAT +H -3.977 3.172 9.701 WAT +O -1.844 1.888 7.641 WAT +H -1.905 0.944 7.494 WAT +H -1.065 1.996 8.186 WAT +O -5.872 0.637 8.335 WAT +H -5.269 1.209 8.810 WAT +H -6.613 1.201 8.113 WAT +O -5.623 0.816 18.715 WAT +H -5.736 1.136 19.610 WAT +H -4.806 1.217 18.418 WAT +O -4.564 -3.755 27.288 WAT +H -5.026 -4.526 27.616 WAT +H -3.752 -3.728 27.793 WAT +O -4.943 -1.781 29.424 WAT +H -4.930 -2.185 28.557 WAT +H -5.495 -1.007 29.321 WAT +O -3.648 -5.060 24.514 WAT +H -3.759 -6.002 24.638 WAT +H -3.651 -4.698 25.400 WAT +O -4.076 -7.871 24.373 WAT +H -4.695 -8.209 23.726 WAT +H -4.094 -8.515 25.081 WAT +O -5.123 -0.706 22.705 WAT +H -5.283 -0.996 23.603 WAT +H -5.289 -1.482 22.170 WAT +O -5.421 -3.197 21.265 WAT +H -6.114 -3.199 20.605 WAT +H -4.679 -3.621 20.834 WAT +O -5.856 -1.571 25.340 WAT +H -5.400 -2.352 25.653 WAT +H -5.517 -0.860 25.884 WAT +O -3.305 -7.333 28.847 WAT +H -2.585 -6.837 28.457 WAT +H -3.507 -8.011 28.202 WAT +O -1.829 -3.402 23.482 WAT +H -2.471 -4.056 23.757 WAT +H -2.216 -2.996 22.706 WAT +O -6.250 -6.073 27.329 WAT +H -5.898 -6.670 26.670 WAT +H -7.134 -6.399 27.500 WAT +O -3.439 -4.510 20.020 WAT +H -3.666 -5.227 20.612 WAT +H -3.830 -4.759 19.182 WAT +O -3.858 -6.593 21.681 WAT +H -4.146 -6.132 22.469 WAT +H -3.062 -7.053 21.950 WAT +O -2.005 -2.429 28.220 WAT +H -1.881 -1.606 28.693 WAT +H -1.139 -2.640 27.873 WAT +O -2.370 -1.936 21.138 WAT +H -2.656 -2.349 20.324 WAT +H -2.568 -1.007 21.018 WAT +O -0.701 -4.544 19.661 WAT +H -1.531 -4.476 20.132 WAT +H -0.523 -3.654 19.358 WAT +O -1.245 -7.799 22.164 WAT +H -0.483 -7.230 22.270 WAT +H -1.230 -8.367 22.935 WAT +O -0.239 -6.944 27.009 WAT +H 0.132 -7.728 26.602 WAT +H 0.401 -6.255 26.830 WAT +O -1.193 -7.506 19.322 WAT +H -0.918 -6.598 19.194 WAT +H -1.382 -7.570 20.258 WAT +O -3.198 -0.809 25.948 WAT +H -2.602 -1.130 25.272 WAT +H -2.940 -1.282 26.740 WAT +O -2.647 -0.159 29.778 WAT +H -3.106 0.592 29.404 WAT +H -3.276 -0.879 29.720 WAT +O -1.176 -0.495 24.144 WAT +H -1.001 -1.182 23.500 WAT +H -1.874 0.029 23.752 WAT +O -4.564 -3.755 8.513 WAT +H -5.026 -4.526 8.842 WAT +H -3.752 -3.728 9.019 WAT +O -2.479 -3.626 11.936 WAT +H -3.162 -2.956 11.927 WAT +H -1.939 -3.428 11.170 WAT +O -4.943 -1.781 10.650 WAT +H -4.930 -2.185 9.782 WAT +H -5.495 -1.007 10.546 WAT +O -5.544 -6.878 11.978 WAT +H -5.127 -7.092 11.144 WAT +H -4.820 -6.621 12.549 WAT +O -1.598 -8.480 12.070 WAT +H -2.426 -8.410 11.595 WAT +H -1.003 -7.894 11.602 WAT +O -3.944 -1.647 13.704 WAT +H -4.421 -0.840 13.511 WAT +H -4.628 -2.299 13.856 WAT +O -0.364 -1.886 18.773 WAT +H 0.112 -1.407 19.452 WAT +H 0.257 -1.954 18.047 WAT +O -0.961 -4.694 15.933 WAT +H -0.834 -4.933 15.014 WAT +H -0.918 -5.526 16.405 WAT +O -5.856 -1.571 6.565 WAT +H -5.400 -2.352 6.879 WAT +H -5.517 -0.860 7.109 WAT +O -3.305 -7.333 10.073 WAT +H -2.585 -6.837 9.682 WAT +H -3.507 -8.011 9.428 WAT +O -1.780 -7.799 15.057 WAT +H -2.506 -7.206 14.860 WAT +H -1.253 -7.807 14.258 WAT +O -6.250 -6.073 8.555 WAT +H -5.898 -6.670 7.895 WAT +H -7.134 -6.399 8.726 WAT +O -7.285 -4.795 11.522 WAT +H -6.643 -5.496 11.635 WAT +H -6.760 -4.015 11.341 WAT +O -5.511 -3.996 14.725 WAT +H -5.176 -3.409 15.403 WAT +H -6.431 -4.126 14.955 WAT +O -2.005 -2.429 9.446 WAT +H -1.881 -1.606 9.918 WAT +H -1.139 -2.640 9.098 WAT +O -6.691 -5.043 18.254 WAT +H -7.297 -5.020 17.513 WAT +H -6.996 -5.775 18.790 WAT +O -5.261 -1.315 16.865 WAT +H -4.309 -1.399 16.815 WAT +H -5.408 -0.651 17.538 WAT +O -3.612 -5.526 17.341 WAT +H -2.871 -5.189 16.838 WAT +H -4.375 -5.357 16.788 WAT +O -4.937 -7.475 15.876 WAT +H -4.523 -7.320 16.725 WAT +H -5.872 -7.349 16.039 WAT +O -1.941 -2.286 15.550 WAT +H -2.541 -2.227 14.806 WAT +H -1.811 -3.225 15.679 WAT +O -0.239 -6.944 8.234 WAT +H 0.132 -7.728 7.828 WAT +H 0.401 -6.255 8.056 WAT +O -3.256 -5.561 13.647 WAT +H -4.070 -5.131 13.910 WAT +H -2.848 -4.949 13.034 WAT +O -3.198 -0.809 7.174 WAT +H -2.602 -1.130 6.498 WAT +H -2.940 -1.282 7.965 WAT +O -0.331 -0.802 13.725 WAT +H -0.825 -0.405 13.008 WAT +H -0.939 -1.427 14.119 WAT +O -2.647 -0.159 11.004 WAT +H -3.106 0.592 10.630 WAT +H -3.276 -0.879 10.945 WAT +O -2.843 -0.811 18.196 WAT +H -2.050 -1.306 17.989 WAT +H -2.549 -0.132 18.804 WAT diff --git a/tests/Fist/regtest-gal21/gcn.inc b/tests/Fist/regtest-gal21/gcn.inc new file mode 100644 index 0000000000..1d28a01e90 --- /dev/null +++ b/tests/Fist/regtest-gal21/gcn.inc @@ -0,0 +1,870 @@ +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +11.25 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +7.5 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 diff --git a/tests/Fist/regtest-gal21/test_gal21.inp b/tests/Fist/regtest-gal21/test_gal21.inp new file mode 100644 index 0000000000..3f98eae268 --- /dev/null +++ b/tests/Fist/regtest-gal21/test_gal21.inp @@ -0,0 +1,148 @@ +&GLOBAL + PROJECT md_gal21 + RUN_TYPE energy +&END GLOBAL + +&MOTION + &PRINT + &TRAJECTORY + &EACH + MD 1 + &END EACH + &END TRAJECTORY + &END PRINT + &CONSTRAINT + &FIXED_ATOMS + LIST 1..192 !Pt atoms + &END FIXED_ATOMS + &G3X3 + DISTANCES 1.808845716 1.808845716 2.860856174 + MOLNAME WAT !To be specified in the coord.inc file + ATOMS 1 2 3 + &END G3X3 + &END CONSTRAINT + &MD + ENSEMBLE NVT + STEPS 10 + TIMESTEP 2 + &THERMOSTAT + &END THERMOSTAT + TEMPERATURE 298 + &END MD +&END MOTION + +&MULTIPLE_FORCE_EVALS + FORCE_EVAL_ORDER 1 2 3 +&END MULTIPLE_FORCE_EVALS + +&FORCE_EVAL + METHOD MIXED + &MIXED + GROUP_PARTITION 1 1 + # MIXING_TYPE LINEAR_COMBINATION + MIXING_TYPE GENMIX + &GENERIC + # X: Energy force_eval 2 + # Y: Energy force_eval 3 + MIXING_FUNCTION X+Y + VARIABLES X Y + &END + &END + &SUBSYS + &CELL + ABC 16.8721334845359721 19.4822616182002832 75 + ALPHA_BETA_GAMMA 90 90 90 + PERIODIC XYZ + &END CELL + &COORD + @INCLUDE "coord.inc" + &END COORD + &TOPOLOGY + &GENERATE + BONDPARM_FACTOR 0.95 + &END GENERATE + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL + +&FORCE_EVAL + METHOD FIST + &MM + &FORCEFIELD + MULTIPLE_POTENTIAL .TRUE. + IGNORE_MISSING_CRITICAL_PARAMS T + &SPLINE + EMAX_SPLINE 10 + &END SPLINE + &CHARGE + ATOM O + CHARGE -0.830 + &END CHARGE + &CHARGE + ATOM Pt + CHARGE 0.0 + &END CHARGE + &CHARGE + ATOM H + CHARGE 0.415 + &END CHARGE + &NONBONDED + &LENNARD-JONES + atoms O O + EPSILON [kcalmol] 0.102 + SIGMA [angstrom] 3.188 + RCUT 10 + &END LENNARD-JONES + &END NONBONDED + &END FORCEFIELD + &POISSON + &EWALD + EWALD_TYPE spme + GMAX 17 20 75 !About box dimensions + &END EWALD + &END POISSON + &END MM +&END FORCE_EVAL + + +&FORCE_EVAL + METHOD FIST + &MM + &FORCEFIELD + &SPLINE + EMAX_SPLINE 10 + &END SPLINE + IGNORE_MISSING_CRITICAL_PARAMS T + &NONBONDED + &GAL21 + ATOMS O Pt + METALS Pt Cu + RCUT 10 + epsilon -0.9358113879523 9.09851523151946 -5.424038077984534 + bxy .6451191 .0010102 + bz .1214708 .0003720 + r 1.9000000 4.0 + a1 -45.03531580496108 2.2912089535351723 0.9952972316569941 + a2 35.580119531527174 -9.91621141032723 0.7804880180017064 + a3 -106.48217251468913 29.06440141041192 -1.8553892798403497 + a4 -14.13835574840837 -1.9537083062459715e-10 0.19129810724417187 + A -48399.51277185474 18486.892036112804 + B 3.4989779 .0988938 + C 917.455856665695 + AH 35728.40150389513 -3751.73899676257 + BH 4.5412501 -.1265630 + Fit_express FALSE + &GCN + @INCLUDE "gcn.inc" + &END GCN + &END GAL21 + &END NONBONDED + &END FORCEFIELD + &POISSON + &EWALD + EWALD_TYPE spme + GMAX 17 20 75 !About box size + &END EWALD + &END POISSON + &END MM +&END FORCE_EVAL diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index c90de5bfe9..9468af631d 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -270,6 +270,7 @@ Fist/regtest-field Fist/regtest-gauss QS/regtest-elpa-2 elpa Fist/regtest-gal19 +Fist/regtest-gal21 QS/regtest-fftw fftw3 QS/regtest-elpa elpa mpiranks==1||mpiranks==2||mpiranks==4||mpiranks==6||mpiranks==8||mpiranks==16 ATOM/regtest-libxc libxc