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

CELL_OPT: implement optimization with fixed cell volume (#5086)

Browse source
This commit is contained in:
HE Zilong 2026-04-22 20:17:54 +08:00 • committed by GitHub
parent d79ec3fd85
commit d749a202a6
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
8 changed files with 306 additions and 25 deletions
Download patch file Download diff file Expand all files Collapse all files

11
src/motion/cell_opt_types.F
View file

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

12
src/motion/gopt_f_methods.F
View file

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

67
src/start/input_cp2k_motion.F
View file

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

1
tests/QS/regtest-spgr/TEST_FILES.toml
View file

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

82
tests/QS/regtest-spgr/quartz-kv-ksg.inp Normal file
View file

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

2
tests/xTB/regtest-1/TEST_FILES.toml
View file

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

79
tests/xTB/regtest-1/quartz-kv-ks.inp Normal file
View file

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

77
tests/xTB/regtest-1/quartz-kv.inp Normal file
View file

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