Add explanation and error message about CELL/SYMMETRY (#4974)

This commit is contained in:
HE Zilong 2026-03-22 05:41:24 +08:00 committed by GitHub
parent 7765b1ffe2
commit 34a1bdb69b
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
3 changed files with 147 additions and 98 deletions

View file

@ -284,9 +284,10 @@ CONTAINS
LOGICAL, INTENT(IN), OPTIONAL :: check_for_ref
TYPE(mp_para_env_type), POINTER :: para_env
REAL(KIND=dp), PARAMETER :: eps = EPSILON(0.0_dp)
REAL(KIND=dp), PARAMETER :: eps = 1.0E-14_dp
CHARACTER(LEN=default_path_length) :: cell_file_name, coord_file_name
CHARACTER(LEN=default_path_length) :: cell_file_name, coord_file_name, &
error_msg
INTEGER :: cell_file_format, coord_file_format, &
my_per
INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell
@ -300,9 +301,12 @@ CONTAINS
my_check_ref = .TRUE.
NULLIFY (cell_ref_section, cell_par, cell_tmp, multiple_unit_cell)
! cell_tmp is only for transferring matrix of cell vectors from
! individual file parser subroutines to read_mat here, assuming
! that unit conversion has been done in those subroutines.
! cell_tmp has two purposes:
! 1. for transferring matrix of cell vectors from individual
! file parser subroutines to read_mat here, assuming that
! unit conversion has been done in those subroutines;
! 2. for testing whether enforcing symmetry makes a new set
! of cell vectors significantly different from parsed input
CALL cell_create(cell_tmp)
IF (.NOT. ASSOCIATED(cell)) CALL cell_create(cell, tag="CELL")
IF (.NOT. ASSOCIATED(cell_ref)) CALL cell_create(cell_ref, tag="CELL_REF")
@ -442,7 +446,6 @@ CONTAINS
END IF
END IF
END IF
CALL cell_release(cell_tmp)
! Convert read_mat or read_len and read_ang to actual cell%hmat
IF (ANY(read_mat(:, :) > eps)) THEN
@ -497,9 +500,22 @@ CONTAINS
! Load requested cell symmetry
CALL section_vals_val_get(cell_section, "SYMMETRY", i_val=cell%symmetry_id)
! Initialize cell
CALL init_cell(cell)
! Try enforcing symmetry by initializing a temporary copy of cell
! and see if the resulting cell matrix differ significantly
CALL cell_clone(cell, cell_tmp)
CALL init_cell(cell_tmp)
IF (ANY(ABS(cell_tmp%hmat - cell%hmat) > eps)) THEN
WRITE (UNIT=error_msg, FMT="(A)") &
"When initializing cell vectors with requested symmetry, one "// &
"or more elements of the cell matrix has varied significantly. "// &
"The input parameters are either deviating from the symmetry, "// &
"or not conforming to the program convention that cell matrix "// &
"is a lower triangle. The symmetrized cell vectors will be used "// &
"anyway with the input atomic coordinates."
CALL cp_warn(__LOCATION__, error_msg)
END IF
CALL cell_clone(cell_tmp, cell)
CALL cell_release(cell_tmp)
IF (my_check_ref) THEN
! Recursive check for reference cell requested
@ -639,12 +655,12 @@ CONTAINS
!> \param para_env ...
!> \date 12.2008
!> \par Format Information implemented:
!> _cell_length_a
!> _cell_length_b
!> _cell_length_c
!> _cell_angle_alpha
!> _cell_angle_beta
!> _cell_angle_gamma
!> _cell_length_a (_cell.length_a)
!> _cell_length_b (_cell.length_b)
!> _cell_length_c (_cell.length_c)
!> _cell_angle_alpha (_cell.length_alpha)
!> _cell_angle_beta (_cell.length_beta)
!> _cell_angle_gamma (_cell.length_gamma)
!>
!> \author Teodoro Laino [tlaino]
!> moved from topology_cif (1/2019 JHU)
@ -670,51 +686,75 @@ CONTAINS
! Parsing cell infos
periodic = 1
! Check for _cell_length_a
! Check for _cell_length_a or _cell.length_a
CALL parser_search_string(parser, "_cell_length_a", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
CPABORT("The field (_cell_length_a) was not found in CIF file! ")
IF (.NOT. found) THEN
CALL parser_search_string(parser, "_cell.length_a", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
CPABORT("The field _cell_length_a or _cell.length_a was not found in CIF file! ")
END IF
CALL cif_get_real(parser, cell_lengths(1))
cell_lengths(1) = cp_unit_to_cp2k(cell_lengths(1), "angstrom")
! Check for _cell_length_b
! Check for _cell_length_b or _cell.length_b
CALL parser_search_string(parser, "_cell_length_b", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
CPABORT("The field (_cell_length_b) was not found in CIF file! ")
IF (.NOT. found) THEN
CALL parser_search_string(parser, "_cell.length_b", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
CPABORT("The field _cell_length_b or _cell.length_b was not found in CIF file! ")
END IF
CALL cif_get_real(parser, cell_lengths(2))
cell_lengths(2) = cp_unit_to_cp2k(cell_lengths(2), "angstrom")
! Check for _cell_length_c
! Check for _cell_length_c or _cell.length_c
CALL parser_search_string(parser, "_cell_length_c", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
CPABORT("The field (_cell_length_c) was not found in CIF file! ")
IF (.NOT. found) THEN
CALL parser_search_string(parser, "_cell.length_c", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
CPABORT("The field _cell_length_c or _cell.length_c was not found in CIF file! ")
END IF
CALL cif_get_real(parser, cell_lengths(3))
cell_lengths(3) = cp_unit_to_cp2k(cell_lengths(3), "angstrom")
! Check for _cell_angle_alpha
! Check for _cell_angle_alpha or _cell.angle_alpha
CALL parser_search_string(parser, "_cell_angle_alpha", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
CPABORT("The field (_cell_angle_alpha) was not found in CIF file! ")
IF (.NOT. found) THEN
CALL parser_search_string(parser, "_cell.angle_alpha", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
CPABORT("The field _cell_angle_alpha or _cell.angle_alpha was not found in CIF file! ")
END IF
CALL cif_get_real(parser, cell_angles(1))
cell_angles(1) = cp_unit_to_cp2k(cell_angles(1), "deg")
! Check for _cell_angle_beta
! Check for _cell_angle_beta or _cell.angle_beta
CALL parser_search_string(parser, "_cell_angle_beta", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
CPABORT("The field (_cell_angle_beta) was not found in CIF file! ")
IF (.NOT. found) THEN
CALL parser_search_string(parser, "_cell.angle_beta", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
CPABORT("The field _cell_angle_beta or _cell.angle_beta was not found in CIF file! ")
END IF
CALL cif_get_real(parser, cell_angles(2))
cell_angles(2) = cp_unit_to_cp2k(cell_angles(2), "deg")
! Check for _cell_angle_gamma
! Check for _cell_angle_gamma or _cell.angle_gamma
CALL parser_search_string(parser, "_cell_angle_gamma", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
CPABORT("The field (_cell_angle_gamma) was not found in CIF file! ")
IF (.NOT. found) THEN
CALL parser_search_string(parser, "_cell.angle_gamma", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
CPABORT("The field _cell_angle_gamma or _cell.angle_gamma was not found in CIF file! ")
END IF
CALL cif_get_real(parser, cell_angles(3))
cell_angles(3) = cp_unit_to_cp2k(cell_angles(3), "deg")

View file

@ -12,28 +12,28 @@
!> \author Manuel Guidon
! **************************************************************************************************
MODULE input_cp2k_hfx
USE bibliography, ONLY: Guidon2008,&
Guidon2009
USE cp_output_handling, ONLY: add_last_numeric,&
cp_print_key_section_create,&
high_print_level,&
medium_print_level
USE input_constants, ONLY: &
do_potential_coulomb, do_potential_gaussian, do_potential_id, do_potential_long, &
do_potential_mix_cl, do_potential_mix_cl_trunc, do_potential_mix_lg, do_potential_short, &
do_potential_truncated, ehrenfest, gaussian, hfx_ri_do_2c_cholesky, hfx_ri_do_2c_diag, &
hfx_ri_do_2c_iter
USE input_keyword_types, ONLY: keyword_create,&
keyword_release,&
keyword_type
USE input_section_types, ONLY: section_add_keyword,&
section_add_subsection,&
section_create,&
section_release,&
section_type
USE input_val_types, ONLY: real_t
USE kinds, ONLY: dp
USE string_utilities, ONLY: s2a
USE bibliography, ONLY: Guidon2008, &
Guidon2009
USE cp_output_handling, ONLY: add_last_numeric, &
cp_print_key_section_create, &
high_print_level, &
medium_print_level
USE input_constants, ONLY: &
do_potential_coulomb, do_potential_gaussian, do_potential_id, do_potential_long, &
do_potential_mix_cl, do_potential_mix_cl_trunc, do_potential_mix_lg, do_potential_short, &
do_potential_truncated, ehrenfest, gaussian, hfx_ri_do_2c_cholesky, hfx_ri_do_2c_diag, &
hfx_ri_do_2c_iter
USE input_keyword_types, ONLY: keyword_create, &
keyword_release, &
keyword_type
USE input_section_types, ONLY: section_add_keyword, &
section_add_subsection, &
section_create, &
section_release, &
section_type
USE input_val_types, ONLY: real_t
USE kinds, ONLY: dp
USE string_utilities, ONLY: s2a
#include "./base/base_uses.f90"
IMPLICIT NONE
@ -229,12 +229,13 @@ CONTAINS
enum_i_vals=[do_potential_coulomb, do_potential_short, do_potential_long, &
do_potential_mix_cl, do_potential_gaussian, do_potential_mix_lg, &
do_potential_id, do_potential_truncated, do_potential_mix_cl_trunc], &
enum_desc=s2a("Coulomb potential: 1/r", &
"Shortrange potential: erfc(omega*r)/r", &
"Longrange potential: erf(omega*r)/r", &
"Mix coulomb and longrange potential: 1/r + erf(omega*r)/r", &
"Damped Gaussian potential: exp(-omega^2*r^2)", &
"Mix Gaussian and longrange potential: erf(omega*r)/r + exp(-omega^2*r^2)", &
enum_desc=s2a("Coulomb potential: $\frac{1}{r}$", &
"Shortrange potential: $\frac{\mathrm{erfc}(\omega \cdot r)}{r}$", &
"Longrange potential: $\frac{\mathrm{erf}(\omega \cdot r)}{r}$", &
"Mix coulomb and longrange potential: $\frac{1}{r} + \frac{\mathrm{erf}(\omega \cdot r)}{r}$", &
"Damped Gaussian potential: $\exp{(-\omega^2 \cdot r^2)}$", &
"Mix Gaussian and longrange potential: "// &
"$\frac{\mathrm{erf}(\omega \cdot r)}{r} + \exp{(-\omega^2 \cdot r^2)}$", &
"Overlap", &
"Truncated coulomb potential: if (r < R_c) 1/r else 0", &
"Truncated Mix coulomb and longrange potential, assumes/requires that the erf has fully decayed at R_c"), &
@ -246,7 +247,7 @@ CONTAINS
CALL keyword_create( &
keyword, __LOCATION__, &
name="OMEGA", &
description="Parameter for short/longrange interaction", &
description="Parameter $\omega$ for short/longrange interaction", &
usage="OMEGA 0.5", &
default_r_val=0.0_dp)
CALL section_add_keyword(section, keyword)
@ -274,10 +275,11 @@ CONTAINS
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="CUTOFF_RADIUS", &
description="Determines cutoff radius (in Angstroms) for the truncated 1/r "// &
"potential or the shortrange erfc(omega*r)/r potential. Default value "// &
"for shortrange potential is solved from erfc(omega*r)/r = EPS_SCHWARZ "// &
"by Newton-Raphson method with SCREENING/EPS_SCHWARZ", &
description="Determines cutoff radius (in Angstroms) for the truncated $\frac{1}{r}$ "// &
"potential or the shortrange $\frac{\mathrm{erfc}(\omega \cdot r)}{r}$ potential. Default value "// &
"for shortrange potential when this keyword is omitted is solved from "// &
"$\frac{\mathrm{erfc}(\omega \cdot r)}{r} = \epsilon_{\mathrm{schwarz}}$ "// &
"by Newton-Raphson method, with $\epsilon_{\mathrm{schwarz}}$ set by SCREENING/EPS_SCHWARZ", &
usage="CUTOFF_RADIUS 10.0", type_of_var=real_t, & ! default_r_val=10.0_dp,&
unit_str="angstrom")
CALL section_add_keyword(section, keyword)

View file

@ -25,36 +25,29 @@ MODULE input_cp2k_subsys
cell_sym_rhombohedral, cell_sym_tetragonal_ab, cell_sym_tetragonal_ac, &
cell_sym_tetragonal_bc, cell_sym_triclinic, use_perd_none, use_perd_x, use_perd_xy, &
use_perd_xyz, use_perd_xz, use_perd_y, use_perd_yz, use_perd_z
USE cp_output_handling, ONLY: cp_print_key_section_create, &
debug_print_level, &
high_print_level, &
medium_print_level
USE cp_output_handling, ONLY: cp_print_key_section_create, debug_print_level, &
high_print_level, medium_print_level
USE cp_units, ONLY: cp_unit_to_cp2k
USE input_constants, ONLY: &
do_add, do_bondparm_covalent, do_bondparm_vdw, &
do_cell_cif, do_cell_cp2k, do_cell_xsc, do_cell_extxyz, do_cell_pdb, &
do_conn_amb7, do_conn_g87, do_conn_g96, do_conn_generate, do_conn_mol_set, do_conn_off, &
do_conn_psf, do_conn_psf_u, do_conn_user, do_coord_cif, do_coord_cp2k, do_coord_crd, &
do_coord_g96, do_coord_off, do_coord_pdb, do_coord_xtl, do_coord_xyz, do_remove, &
do_skip_11, do_skip_12, do_skip_13, do_skip_14, dump_pdb, gaussian
USE input_constants, ONLY: do_add, do_bondparm_covalent, do_bondparm_vdw, &
do_cell_cif, do_cell_cp2k, do_cell_xsc, &
do_cell_extxyz, do_cell_pdb, &
do_conn_amb7, do_conn_g87, do_conn_g96, &
do_conn_generate, do_conn_mol_set, do_conn_off, &
do_conn_psf, do_conn_psf_u, do_conn_user, &
do_coord_cif, do_coord_cp2k, do_coord_crd, &
do_coord_g96, do_coord_off, do_coord_pdb, &
do_coord_xtl, do_coord_xyz, do_remove, &
do_skip_11, do_skip_12, do_skip_13, do_skip_14, &
dump_pdb, gaussian
USE input_cp2k_colvar, ONLY: create_colvar_section
USE input_cp2k_mm, ONLY: create_neighbor_lists_section
USE input_keyword_types, ONLY: keyword_create, &
keyword_release, &
keyword_type
USE input_section_types, ONLY: section_add_keyword, &
section_add_subsection, &
section_create, &
section_release, &
section_type
USE input_val_types, ONLY: char_t, &
integer_t, &
lchar_t, &
real_t
USE input_keyword_types, ONLY: keyword_create, keyword_release, keyword_type
USE input_section_types, ONLY: section_add_keyword, section_add_subsection, &
section_create, section_release, section_type
USE input_val_types, ONLY: char_t, integer_t, lchar_t, real_t
USE kinds, ONLY: dp
USE physcon, ONLY: bohr
USE string_utilities, ONLY: newline, &
s2a
USE string_utilities, ONLY: newline, s2a
#include "./base/base_uses.f90"
IMPLICIT NONE
@ -170,7 +163,7 @@ CONTAINS
description="Specify the angles between the vectors A, B and C when using the ABC keyword. "// &
"The convention is that A lies along the X-axis, B is in the XY plane. "// &
"ALPHA is the angle between B and C, BETA is the angle between A and C and "// &
"GAMMA the angle between A and B.", &
"GAMMA is the angle between A and B.", &
usage="ALPHA_BETA_GAMMA [deg] 90.0 90.0 120.0", unit_str="deg", &
n_var=3, default_r_vals=[cp_unit_to_cp2k(value=90.0_dp, unit_str="deg"), &
cp_unit_to_cp2k(value=90.0_dp, unit_str="deg"), &
@ -180,23 +173,29 @@ CONTAINS
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="CELL_FILE_NAME", &
description="Possibility to read the cell from an external file ", &
description="The external file from which cell is parsed ", &
repeats=.FALSE., usage="CELL_FILE_NAME <CHARACTER>", &
type_of_var=lchar_t)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="CELL_FILE_FORMAT", &
description="Specify the format of the cell file (if used)", &
description="Format of the external file from which "// &
"cell is parsed. If the format specifies a cell by "// &
"lengths and angles of three vectors, then a cell "// &
"matrix is constructed with the convention that A "// &
"lies along the X-axis, B is in the XY plane. ALPHA "// &
"is the angle between B and C, BETA is the angle "// &
"between A and C, and GAMMA is the angle between A and B.", &
usage="CELL_FILE_FORMAT (CP2K|CIF|XSC|EXTXYZ|XYZ|PDB)", &
enum_c_vals=s2a("CP2K", "CIF", "XSC", "EXTXYZ", "XYZ", "PDB"), &
enum_i_vals=[do_cell_cp2k, do_cell_cif, do_cell_xsc, do_cell_extxyz, do_cell_extxyz, do_cell_pdb], &
enum_desc=s2a("Cell info in the CP2K native format", &
"Cell info from CIF file", &
"Cell info from CIF file (from fields `_cell_length_a` or `_cell.length_a`, etc)", &
"Cell info in the XSC format (NAMD)", &
"Cell info in the comment line of Extended XYZ format", &
"Cell info as `lattice=...` field in the comment line of Extended XYZ format", &
"Alias for Extended XYZ", &
"Cell info in the CRYST1 record of PDB format"), &
"Cell info in the `CRYST1` record of PDB format"), &
default_i_val=do_cell_cp2k)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
@ -227,7 +226,15 @@ CONTAINS
CALL keyword_create( &
keyword, __LOCATION__, name="SYMMETRY", &
description="Imposes an initial cell symmetry.", &
description="Imposes an initial cell symmetry, according to the convention "// &
"that A lies along the X-axis, B is in the XY plane. After the "// &
"input cell information is parsed, the symmetry is enforced by "// &
"reconstructing the cell matrix from lengths and angles of the "// &
"cell vectors, taking averages if necessary. This process does "// &
"not affect input atomic coordinates; in case a space group is "// &
"to be detected and preserved for an optimization task, atomic "// &
"coordinates should correspond to cell vectors already obeying "// &
"the convention mentioned above.", &
usage="SYMMETRY monoclinic", &
enum_desc=s2a("No cell symmetry", &
"Triclinic (a &ne; b &ne; c &ne; a, &alpha; &ne; &beta; &ne; &gamma; &ne; &alpha; &ne; 90&deg;)", &