From d749a202a69c8b7a7ce4137b3ebb279304492fcb Mon Sep 17 00:00:00 2001 From: HE Zilong <159878975+he-zilong@users.noreply.github.com> Date: Wed, 22 Apr 2026 20:17:54 +0800 Subject: [PATCH] CELL_OPT: implement optimization with fixed cell volume (#5086) --- src/motion/cell_opt_types.F | 11 +++- src/motion/gopt_f_methods.F | 12 +++- src/start/input_cp2k_motion.F | 67 +++++++++++++------- tests/QS/regtest-spgr/TEST_FILES.toml | 1 + tests/QS/regtest-spgr/quartz-kv-ksg.inp | 82 +++++++++++++++++++++++++ tests/xTB/regtest-1/TEST_FILES.toml | 2 + tests/xTB/regtest-1/quartz-kv-ks.inp | 79 ++++++++++++++++++++++++ tests/xTB/regtest-1/quartz-kv.inp | 77 +++++++++++++++++++++++ 8 files changed, 306 insertions(+), 25 deletions(-) create mode 100644 tests/QS/regtest-spgr/quartz-kv-ksg.inp create mode 100644 tests/xTB/regtest-1/quartz-kv-ks.inp create mode 100644 tests/xTB/regtest-1/quartz-kv.inp diff --git a/src/motion/cell_opt_types.F b/src/motion/cell_opt_types.F index cd7bb08e2d..d2bd73dd6d 100644 --- a/src/motion/cell_opt_types.F +++ b/src/motion/cell_opt_types.F @@ -60,7 +60,8 @@ MODULE cell_opt_types ! Simulation cell optimization parameters INTEGER :: constraint_id = fix_none LOGICAL :: keep_angles = .FALSE., & - keep_symmetry = .FALSE. + keep_symmetry = .FALSE., & + keep_volume = .FALSE. REAL(KIND=dp) :: pres_ext = 0.0_dp, pres_int = 0.0_dp, pres_tol = 0.0_dp, pres_constr = 0.0_dp REAL(KIND=dp), DIMENSION(3, 3) :: mtrx = 0.0_dp REAL(KIND=dp), DIMENSION(3, 3) :: rot_matrix = 0.0_dp @@ -95,6 +96,7 @@ CONTAINS CALL force_env_get(force_env, cell=cell, subsys=subsys) CALL cell_create(cell_env%ref_cell) CALL cell_clone(cell, cell_env%ref_cell, tag="REF_CELL_OPT") + CALL section_vals_val_get(geo_section, "KEEP_VOLUME", l_val=cell_env%keep_volume) CALL section_vals_val_get(geo_section, "KEEP_ANGLES", l_val=cell_env%keep_angles) CALL section_vals_val_get(geo_section, "KEEP_SYMMETRY", l_val=cell_env%keep_symmetry) CALL section_vals_val_get(geo_section, "PRESSURE_TOLERANCE", r_val=cell_env%pres_tol) @@ -121,6 +123,13 @@ CONTAINS IF (output_unit > 0) THEN WRITE (UNIT=output_unit, FMT="(/,T2,A,T61,F20.1)") & "CELL_OPT| Pressure tolerance [bar]: ", cp_unit_from_cp2k(cell_env%pres_tol, "bar") + IF (cell_env%keep_volume) THEN + WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") & + "CELL_OPT| Keep volume of cell: ", "YES" + ELSE + WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") & + "CELL_OPT| Keep volume of cell: ", " NO" + END IF IF (cell_env%keep_angles) THEN WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") & "CELL_OPT| Keep angles between the cell vectors: ", "YES" diff --git a/src/motion/gopt_f_methods.F b/src/motion/gopt_f_methods.F index f8152f93a3..65d02f8084 100644 --- a/src/motion/gopt_f_methods.F +++ b/src/motion/gopt_f_methods.F @@ -55,6 +55,7 @@ MODULE gopt_f_methods dp, & int_8 USE machine, ONLY: m_flush + USE mathlib, ONLY: det_3x3 USE md_energies, ONLY: sample_memory USE message_passing, ONLY: mp_para_env_type USE motion_utils, ONLY: write_simulation_cell, & @@ -1081,7 +1082,7 @@ CONTAINS INTEGER :: i, iatom, idg, j, natom, nparticle, & shell_index - REAL(KIND=dp) :: fc, fs, mass + REAL(KIND=dp) :: fc, fs, mass, ratio REAL(KIND=dp), DIMENSION(3) :: s TYPE(cell_type), POINTER :: cell_ref TYPE(cp_subsys_type), POINTER :: subsys @@ -1136,6 +1137,15 @@ CONTAINS cell%hmat(j, i) = x(idg) END DO END DO + IF (gopt_env%cell_env%keep_volume) THEN + ratio = cell%deth/ABS(det_3x3(cell%hmat)) + ratio = ratio**(1.0_dp/3.0_dp) ! no need to use EXP(LOG(ratio)/3) for cube root + DO i = 1, 3 + DO j = 1, i + cell%hmat(j, i) = cell%hmat(j, i)*ratio + END DO + END DO + END IF CALL init_cell(cell) CALL cp_subsys_set(subsys, cell=cell) diff --git a/src/start/input_cp2k_motion.F b/src/start/input_cp2k_motion.F index 46bac74be3..8c41eec28f 100644 --- a/src/start/input_cp2k_motion.F +++ b/src/start/input_cp2k_motion.F @@ -1013,26 +1013,32 @@ CONTAINS TYPE(section_type), POINTER :: print_key, subsection CALL create_geoopt_section(section, __LOCATION__, label="CELL_OPT", & - description="This section sets the environment for the optimization of the simulation cell."// & - " Two possible schemes are available: (1) Zero temperature optimization;"// & - " (2) Finite temperature optimization.", & - just_optimizers=.TRUE., & - use_model_hessian=.FALSE.) + description="This section sets the environment for the optimization "// & + "of the simulation cell. As is noted in FORCE_EVAL/SUBSYS/CELL, the "// & + "program convention is that the first cell vector A lies along the "// & + "X-axis and the second cell vector B is in the XY plane, such that "// & + "the cell vector matrix is a lower triangle. The algorithm support "// & + "for updating the three upper triangular components during a cell "// & + "optimization is not complete or tested, so the input structure has "// & + "to be prepared accordingly with these three components precisely 0 "// & + "even for cases like the primitive rhombohedral cell of the FCC lattice.", & + just_optimizers=.TRUE., use_model_hessian=.FALSE.) NULLIFY (keyword, print_key, subsection) CALL keyword_create( & keyword, __LOCATION__, name="TYPE", & description="Specify which kind of method to use for the optimization of the simulation cell", & - usage="TYPE (GEO_OPT|MD|DIRECT_CELL_OPT)", & - enum_c_vals=s2a("GEO_OPT", "MD", "DIRECT_CELL_OPT"), & + usage="TYPE (DIRECT_CELL_OPT|GEO_OPT|MD)", & + enum_c_vals=s2a("DIRECT_CELL_OPT", "GEO_OPT", "MD"), & enum_desc=s2a( & - "Performs a geometry optimization (the GEO_OPT section must be defined) between cell optimization steps."// & - " The stress tensor is computed at the optimized geometry.", & - "Performs a molecular dynamics run (the MD section needs must defined) for computing the stress tensor"// & - " used for the cell optimization.", & - "Performs a geometry and cell optimization at the same time."// & - " The stress tensor is computed at every step"), & - enum_i_vals=[default_cell_geo_opt_id, default_cell_md_id, default_cell_direct_id], & + "Performs a geometry and cell optimization at the same time. "// & + "The stress tensor is computed at every step", & + "Performs a geometry optimization between cell optimization steps. "// & + "The MOTION/GEO_OPT section must be defined. "// & + "The stress tensor is computed at the optimized geometry.", & + "Performs a molecular dynamics run for computing the stress tensor "// & + "used for the cell optimization. The MOTION/MD section must be defined."), & + enum_i_vals=[default_cell_direct_id, default_cell_geo_opt_id, default_cell_md_id], & default_i_val=default_cell_direct_id) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -1048,19 +1054,34 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create( & - keyword, __LOCATION__, name="KEEP_ANGLES", & - description="Keep angles between the cell vectors constant, but allow the lengths of the"// & - " cell vectors to change independently."// & - " Albeit general, this is most useful for triclinic cells, to enforce higher symmetry, see KEEP_SYMMETRY.", & - usage="KEEP_ANGLES TRUE", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL keyword_create(keyword, __LOCATION__, name="KEEP_VOLUME", & + description="Keep the volume of the cell constant during cell optimization. "// & + "This is implemented by comparing the cell volumes and scaling the new "// & + "cell vectors just before updating the cell information, and can be "// & + "used together with KEEP_ANGLES or KEEP_SYMMETRY.", & + usage="KEEP_VOLUME TRUE", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="KEEP_ANGLES", & + description="Keep angles between the cell vectors constant, but "// & + "allow the lengths of the cell vectors to change independently "// & + "during cell optimization. This is implemented by projecting out "// & + "the components of angles in the cell gradient before the cell "// & + "is updated. Albeit general, this is most useful for triclinic "// & + "cells; to enforce higher symmetry, see KEEP_SYMMETRY.", & + usage="KEEP_ANGLES TRUE", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) CALL keyword_create(keyword, __LOCATION__, name="KEEP_SYMMETRY", & - description="Keep the requested initial cell symmetry (e.g. during a cell optimisation). "// & - "The initial symmetry must be specified in the &CELL section.", & - usage="KEEP_SYMMETRY yes", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + description="Keep the requested initial cell symmetry as specified "// & + "in the FORCE_EVAL/SUBSYS/CELL section during cell optimization. "// & + "This is implemented by removing symmetry-breaking components and "// & + "taking averages of components if necessary in the cell gradient "// & + "before the cell is updated. To enforce the space group (which "// & + "requires spglib package), see KEEP_SPACE_GROUP.", & + usage="KEEP_SYMMETRY TRUE", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) diff --git a/tests/QS/regtest-spgr/TEST_FILES.toml b/tests/QS/regtest-spgr/TEST_FILES.toml index 4217f5d036..b33f2c50e1 100644 --- a/tests/QS/regtest-spgr/TEST_FILES.toml +++ b/tests/QS/regtest-spgr/TEST_FILES.toml @@ -5,4 +5,5 @@ # for details see cp2k/tools/do_regtest "copper.inp" = [{matcher="M093", tol=0, ref=225}] "quartz.inp" = [{matcher="M093", tol=0, ref=154}] +"quartz-kv-ksg.inp" = [{matcher="E_total", tol=1.0E-12, ref=-34.85757998070905}] #EOF diff --git a/tests/QS/regtest-spgr/quartz-kv-ksg.inp b/tests/QS/regtest-spgr/quartz-kv-ksg.inp new file mode 100644 index 0000000000..d66479ce59 --- /dev/null +++ b/tests/QS/regtest-spgr/quartz-kv-ksg.inp @@ -0,0 +1,82 @@ +&GLOBAL + PRINT_LEVEL MEDIUM + PROJECT_NAME quartz-kv-ksg + RUN_TYPE CELL_OPT +&END GLOBAL + +&MOTION + &CELL_OPT + EPS_SYMMETRY 1e-4 + KEEP_ANGLES T + KEEP_SPACE_GROUP T + KEEP_SYMMETRY T + KEEP_VOLUME T + MAX_ITER 5 + OPTIMIZER BFGS + &END CELL_OPT + &PRINT + &CELL ON + &END CELL + &RESTART_HISTORY OFF + &END RESTART_HISTORY + &END PRINT +&END MOTION + +&FORCE_EVAL + METHOD QS + STRESS_TENSOR ANALYTICAL + &DFT + CHARGE 0 + MULTIPLICITY 1 + &POISSON + PERIODIC XYZ + POISSON_SOLVER PERIODIC + &END POISSON + &QS + METHOD XTB + &XTB + CHECK_ATOMIC_CHARGES F + DO_EWALD T + &END XTB + &END QS + &SCF + EPS_SCF 1.E-5 + MAX_SCF 10 + SCF_GUESS ATOMIC + &OT T + MINIMIZER CG + PRECONDITIONER FULL_ALL + &END OT + &OUTER_SCF T + EPS_SCF 1.E-5 + MAX_SCF 10 + &END OUTER_SCF + &END SCF + &END DFT + &SUBSYS + &CELL + ABC 4.9159998894 4.9159998894 5.4053997993 + ALPHA_BETA_GAMMA 90.00 90.00 120.00 + MULTIPLE_UNIT_CELL 1 1 1 + PERIODIC XYZ + SYMMETRY MONOCLINIC_GAMMA_AB + &END CELL + &COORD + Si 0.469699985 0.000000000 0.000000000 + Si 0.000000018 0.469700021 0.666666696 + Si 0.530300031 0.530300048 0.333333348 + O 0.413500032 0.266900007 0.119099996 + O 0.266900018 0.413500038 0.547566700 + O 0.733100016 0.146600017 0.785766692 + O 0.586499976 0.853400010 0.214233330 + O 0.853400027 0.586500031 0.452433344 + O 0.146600008 0.733100006 0.880900004 + UNIT angstrom + SCALED T + &END COORD + &TOPOLOGY + MULTIPLE_UNIT_CELL 1 1 1 + NUMBER_OF_ATOMS 9 + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/xTB/regtest-1/TEST_FILES.toml b/tests/xTB/regtest-1/TEST_FILES.toml index dc3c0813d0..77d44f4aef 100644 --- a/tests/xTB/regtest-1/TEST_FILES.toml +++ b/tests/xTB/regtest-1/TEST_FILES.toml @@ -20,4 +20,6 @@ "h2o_dimer.inp" = [{matcher="E_total", tol=1.0E-12, ref=-11.54506384435821}] "ghost.inp" = [{matcher="E_total", tol=1.0E-12, ref=-1.03458221876736}] "tcif.inp" = [{matcher="E_total", tol=1.0E-12, ref=-27.71649842848894}] +"quartz-kv.inp" = [{matcher="E_total", tol=1.0E-12, ref=-34.8491812462054}] +"quartz-kv-ks.inp" = [{matcher="E_total", tol=1.0E-12, ref=-34.85757891946725}] #EOF diff --git a/tests/xTB/regtest-1/quartz-kv-ks.inp b/tests/xTB/regtest-1/quartz-kv-ks.inp new file mode 100644 index 0000000000..ef80cc2876 --- /dev/null +++ b/tests/xTB/regtest-1/quartz-kv-ks.inp @@ -0,0 +1,79 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT_NAME quartz-kv-ks + RUN_TYPE CELL_OPT +&END GLOBAL + +&MOTION + &CELL_OPT + KEEP_SYMMETRY T + KEEP_VOLUME T + MAX_ITER 5 + OPTIMIZER BFGS + &END CELL_OPT + &PRINT + &CELL ON + &END CELL + &RESTART_HISTORY OFF + &END RESTART_HISTORY + &END PRINT +&END MOTION + +&FORCE_EVAL + METHOD QS + STRESS_TENSOR ANALYTICAL + &DFT + CHARGE 0 + MULTIPLICITY 1 + &POISSON + PERIODIC XYZ + POISSON_SOLVER PERIODIC + &END POISSON + &QS + METHOD XTB + &XTB + CHECK_ATOMIC_CHARGES F + DO_EWALD T + &END XTB + &END QS + &SCF + EPS_SCF 1.E-5 + MAX_SCF 10 + SCF_GUESS ATOMIC + &OT T + MINIMIZER CG + PRECONDITIONER FULL_ALL + &END OT + &OUTER_SCF T + EPS_SCF 1.E-5 + MAX_SCF 10 + &END OUTER_SCF + &END SCF + &END DFT + &SUBSYS + &CELL + ABC 4.9159998894 4.9159998894 5.4053997993 + ALPHA_BETA_GAMMA 90.00 90.00 120.00 + MULTIPLE_UNIT_CELL 1 1 1 + PERIODIC XYZ + SYMMETRY HEXAGONAL_GAMMA_120 + &END CELL + &COORD + Si 0.469699985 0.000000000 0.000000000 + Si 0.000000018 0.469700021 0.666666696 + Si 0.530300031 0.530300048 0.333333348 + O 0.413500032 0.266900007 0.119099996 + O 0.266900018 0.413500038 0.547566700 + O 0.733100016 0.146600017 0.785766692 + O 0.586499976 0.853400010 0.214233330 + O 0.853400027 0.586500031 0.452433344 + O 0.146600008 0.733100006 0.880900004 + UNIT angstrom + SCALED T + &END COORD + &TOPOLOGY + MULTIPLE_UNIT_CELL 1 1 1 + NUMBER_OF_ATOMS 9 + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/xTB/regtest-1/quartz-kv.inp b/tests/xTB/regtest-1/quartz-kv.inp new file mode 100644 index 0000000000..6300cac587 --- /dev/null +++ b/tests/xTB/regtest-1/quartz-kv.inp @@ -0,0 +1,77 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT_NAME quartz-kv + RUN_TYPE CELL_OPT +&END GLOBAL + +&MOTION + &CELL_OPT + KEEP_VOLUME T + MAX_ITER 5 + OPTIMIZER BFGS + &END CELL_OPT + &PRINT + &CELL ON + &END CELL + &RESTART_HISTORY OFF + &END RESTART_HISTORY + &END PRINT +&END MOTION + +&FORCE_EVAL + METHOD QS + STRESS_TENSOR ANALYTICAL + &DFT + CHARGE 0 + MULTIPLICITY 1 + &POISSON + PERIODIC XYZ + POISSON_SOLVER PERIODIC + &END POISSON + &QS + METHOD XTB + &XTB + CHECK_ATOMIC_CHARGES F + DO_EWALD T + &END XTB + &END QS + &SCF + EPS_SCF 1.E-5 + MAX_SCF 10 + SCF_GUESS ATOMIC + &OT T + MINIMIZER CG + PRECONDITIONER FULL_ALL + &END OT + &OUTER_SCF T + EPS_SCF 1.E-5 + MAX_SCF 10 + &END OUTER_SCF + &END SCF + &END DFT + &SUBSYS + &CELL + ABC 4.9159998894 4.9159998894 5.4053997993 + ALPHA_BETA_GAMMA 90.00 90.00 120.00 + MULTIPLE_UNIT_CELL 1 1 1 + PERIODIC XYZ + &END CELL + &COORD + Si 0.469699985 0.000000000 0.000000000 + Si 0.000000018 0.469700021 0.666666696 + Si 0.530300031 0.530300048 0.333333348 + O 0.413500032 0.266900007 0.119099996 + O 0.266900018 0.413500038 0.547566700 + O 0.733100016 0.146600017 0.785766692 + O 0.586499976 0.853400010 0.214233330 + O 0.853400027 0.586500031 0.452433344 + O 0.146600008 0.733100006 0.880900004 + UNIT angstrom + SCALED T + &END COORD + &TOPOLOGY + MULTIPLE_UNIT_CELL 1 1 1 + NUMBER_OF_ATOMS 9 + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL