diff --git a/src/cell_methods.F b/src/cell_methods.F index 24c25d20aa..8c260362ef 100644 --- a/src/cell_methods.F +++ b/src/cell_methods.F @@ -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") diff --git a/src/input_cp2k_hfx.F b/src/input_cp2k_hfx.F index 3bbfd836ec..79fdfb9297 100644 --- a/src/input_cp2k_hfx.F +++ b/src/input_cp2k_hfx.F @@ -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) diff --git a/src/input_cp2k_subsys.F b/src/input_cp2k_subsys.F index e20e6121d7..1f457e3f7b 100644 --- a/src/input_cp2k_subsys.F +++ b/src/input_cp2k_subsys.F @@ -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 ", & 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 ≠ b ≠ c ≠ a, α ≠ β ≠ γ ≠ α ≠ 90°)", &