comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 22:25:32 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Fist: Add GAL21 forcefield

Browse source
This commit is contained in:
pclabaut 2021-07-09 14:51:05 +02:00 • committed by GitHub
parent 7a55edd99b
commit c923b2c4bf
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
15 changed files with 3615 additions and 63 deletions
Download patch file Download diff file Expand all files Collapse all files

14
src/common/bibliography.F
View file

@ -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", &

6
src/fist_neighbor_list_control.F
View file

@ -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

7
src/fist_nonbond_env_types.F
View file

@ -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.

15
src/fist_nonbond_force.F
View file

@ -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

145
src/force_fields_input.F
View file

@ -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 ...

157
src/input_cp2k_mm.F
View file

@ -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 ...

937
src/manybody_gal21.F Normal file
View file

@ -0,0 +1,937 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2021 CP2K developers group <https://cp2k.org> !
! !
! 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

186
src/manybody_potential.F
View file

@ -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

43
src/pair_potential.F
View file

@ -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

272
src/pair_potential_types.F
View file

@ -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

7
tests/Fist/regtest-gal21/TEST_FILES Normal file
View file

@ -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

870
tests/Fist/regtest-gal21/coord.inc Normal file
View file

@ -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

870
tests/Fist/regtest-gal21/gcn.inc Normal file
View file

@ -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

148
tests/Fist/regtest-gal21/test_gal21.inp Normal file
View file

@ -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

1
tests/TEST_DIRS
View file

@ -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