From 5a536fbc1c7f40a6c4774b9f08c2af48e003f474 Mon Sep 17 00:00:00 2001 From: David Bowler Date: Mon, 20 Jan 2025 17:11:01 +0900 Subject: [PATCH 1/2] Updates to cell constraints and how they are applied Specifically related to constraining ratios e.g. c/a --- src/force_module.f90 | 26 +++++++++++++++++++++++++- src/initial_read_module.f90 | 8 ++++---- src/move_atoms.module.f90 | 20 +++++++++++++++++--- 3 files changed, 46 insertions(+), 8 deletions(-) diff --git a/src/force_module.f90 b/src/force_module.f90 index 1a0534099..ef752c71b 100644 --- a/src/force_module.f90 +++ b/src/force_module.f90 @@ -225,6 +225,8 @@ module force_module !! Add printint forces/stress in ASE output !! 2023/01/10 18:41 lionel !! Secure ASE printing when using ordern + !! 2025/01/20 17:07 dave + !! Add constraints on stress when constraining cell optimisation !! SOURCE !! subroutine force(fixed_potential, vary_mu, n_cg_L_iterations, & @@ -288,7 +290,7 @@ subroutine force(fixed_potential, vary_mu, n_cg_L_iterations, & ! Local variables integer :: i, j, ii, stat, max_atom, max_compt, ispin, & direction, dir1, dir2, counter - real(double) :: max_force, volume, scale, g0 + real(double) :: max_force, volume, scale, g0, scaleC type(cq_timer) :: tmr_l_tmp1 type(cq_timer) :: backtrace_timer integer :: backtrace_level @@ -740,6 +742,28 @@ subroutine force(fixed_potential, vary_mu, n_cg_L_iterations, & else if (leqi(cell_constraint_flag, 'b c') .or. leqi(cell_constraint_flag, 'c b')) then stress(2,2) = zero stress(3,3) = zero + ! These ensure that the stresses maintain the desired ratio while keeping the average fixed + else if (leqi(cell_constraint_flag,'a/b') .or. leqi(cell_constraint_flag,'b/a')) then + call print_stress(trim(prefix)//" Orig stress: ", stress, -2, write_ase) ! Force output + ! Desired ratio + scaleC = rcelly/rcellx + ! Average x-y stress + stress(1,1) = (stress(1,1) + stress(2,2))/(one + scaleC) + stress(2,2) = scaleC*stress(1,1) + else if (leqi(cell_constraint_flag,'a/c') .or. leqi(cell_constraint_flag,'c/a')) then + call print_stress(trim(prefix)//" Orig stress: ", stress, -2, write_ase) ! Force output + ! Desired ratio + scaleC = rcellz/rcellx + ! Average x-z stress + stress(1,1) = (stress(1,1) + stress(3,3))/(one + scaleC) + stress(3,3) = scaleC*stress(1,1) + else if (leqi(cell_constraint_flag,'c/b') .or. leqi(cell_constraint_flag,'b/c')) then + call print_stress(trim(prefix)//" Orig stress: ", stress, -2, write_ase) ! Force output + ! Desired ratio + scaleC = rcelly/rcellz + ! Average y-z stress + stress(3,3) = (stress(3,3) + stress(2,2))/(one + scaleC) + stress(2,2) = scaleC*stress(3,3) end if ! Output if (inode == ionode.AND.iprint_MD + min_layer>2) then diff --git a/src/initial_read_module.f90 b/src/initial_read_module.f90 index 38b367936..dc3db8d17 100644 --- a/src/initial_read_module.f90 +++ b/src/initial_read_module.f90 @@ -1632,10 +1632,10 @@ subroutine read_input(start, start_L, titles, vary_mu,& optcell_method = fdf_integer('AtomMove.OptCellMethod', 1) cell_constraint_flag = fdf_string(20,'AtomMove.OptCell.Constraint','none') ! Warn user if applying constraints with OptCellMethod 3 - if(optcell_method==3.and.(.not.leqi(cell_constraint_flag,'none'))) then - call cq_warn(sub_name,"Cell constraints NOT applied for OptCellMethod 3") - cell_constraint_flag = 'none' - end if + !if(optcell_method==3.and.(.not.leqi(cell_constraint_flag,'none'))) then + ! call cq_warn(sub_name,"Cell constraints NOT applied for OptCellMethod 3") + ! cell_constraint_flag = 'none' + !end if cell_en_tol = fdf_double('AtomMove.OptCell.EnTol',0.00001_double) ! It makes sense to use GPa here so I'm changing the default to 0.1GPa cell_stress_tol = fdf_double('AtomMove.StressTolerance',0.1_double) !005_double) diff --git a/src/move_atoms.module.f90 b/src/move_atoms.module.f90 index fabc9cb9c..764104b1d 100644 --- a/src/move_atoms.module.f90 +++ b/src/move_atoms.module.f90 @@ -5183,12 +5183,15 @@ end subroutine update_pos_and_matrices !! CREATION DATE !! 2019/02/08 !! MODIFICATION HISTORY - !! + !! 2025/01/20 13:45 dave + !! Added conditions for fixed cell side ratios (average stresses) !! SOURCE subroutine propagate_vector(force, config, config_new, cell_ref, k) use GenComms, only: inode, ionode - use global_module, only: iprint_MD, ni_in_cell, id_glob + use global_module, only: iprint_MD, ni_in_cell, id_glob, cell_constraint_flag + use numbers, only: half + use input_module, only: leqi implicit none @@ -5210,7 +5213,18 @@ subroutine propagate_vector(force, config, config_new, cell_ref, k) do j=1,3 config_new(j,i) = config(j,i) + k*force(j,i) end do - end do + end do + ! To maintain cell ratios the strains must be equal, so we average them (the most general way) + if (leqi(cell_constraint_flag, 'a/b') .OR. leqi(cell_constraint_flag, 'b/a')) then + config_new(1,ni_in_cell+1) = half*(config_new(1,ni_in_cell+1) + config_new(2,ni_in_cell+1)) + config_new(2,ni_in_cell+1) = config_new(1,ni_in_cell+1) + else if (leqi(cell_constraint_flag, 'a/c') .OR. leqi(cell_constraint_flag, 'c/a')) then + config_new(1,ni_in_cell+1) = half*(config_new(1,ni_in_cell+1) + config_new(3,ni_in_cell+1)) + config_new(3,ni_in_cell+1) = config_new(1,ni_in_cell+1) + else if (leqi(cell_constraint_flag, 'c/b') .OR. leqi(cell_constraint_flag, 'b/c')) then + config_new(3,ni_in_cell+1) = half*(config_new(3,ni_in_cell+1) + config_new(2,ni_in_cell+1)) + config_new(2,ni_in_cell+1) = config_new(3,ni_in_cell+1) + end if end subroutine propagate_vector !!*** From ac5aeed01aae56f7e42507f4dbb7f71dbc969f79 Mon Sep 17 00:00:00 2001 From: David Bowler Date: Tue, 21 Jan 2025 11:22:15 +0900 Subject: [PATCH 2/2] Updated documentation for constrained cell optimisation --- docs/input_tags.rst | 17 +++++++++-------- docs/strucrelax.rst | 14 +++++++++++++- 2 files changed, 22 insertions(+), 9 deletions(-) diff --git a/docs/input_tags.rst b/docs/input_tags.rst index b3180c33c..62479fd57 100644 --- a/docs/input_tags.rst +++ b/docs/input_tags.rst @@ -951,25 +951,26 @@ AtomMove.OptCell.Constraint (*string*) *Fixing a single cell dimension:* - a: Fix the x-dimension of the simulation box + ``a``: Fix the x-dimension of the simulation cell - b: Fix the y-dimension of the simulation box + ``b``: Fix the y-dimension of the simulation cell - c: Fix the z-dimension of the simulation box + ``c``: Fix the z-dimension of the simulation cell *Fixing multiple cell dimensions:* - any combination of the above separated by a space character. e.g: "a b" fixes - both the x and y dimensions of the simulation box + + Any combination of the above separated by a space character. e.g: ``a b`` fixes + both the x- and y-dimensions of the simulation cell. *Fixing Ratios:* - Any combination of a, b or c separated by a "/" character. e.g "c/a" fixes - the initial ratio of the z-dimension to the x-direction. + Any combination of a, b or c separated by a "/" character, e.g ``c/a`` fixes + the initial ratio of the z-dimension to the x-dimension. *Global scaling factor:* - volume: minimize the total energy by scaling each simulation box dimension by + ``volume``: minimize the total energy by scaling each simulation cell dimension by the same global scaling factor. Search directions are set by the mean stress. AtomMove.TestSpecificForce (*integer*) diff --git a/docs/strucrelax.rst b/docs/strucrelax.rst index 2f0832f47..3b89d843e 100644 --- a/docs/strucrelax.rst +++ b/docs/strucrelax.rst @@ -123,6 +123,15 @@ coordinates* (``AtomMove.OptCellMethod 1``) using the following input: Note that stress is in GPa and enthalpy is in Ha by default. +It is possible to apply constraints to the cell when optimising it, using the +flag ``AtomMove.OptCell.Constraint``. The constraint takes three different possible +forms: fixing one or two of the cell lengths (e.g. ``AtomMove.OptCell.Constraint a`` or +``AtomMove.OptCell.Constraint a b``); fixing the ratio between two cell lengths +(e.g. ``AtomMove.OptCell.Constraint c/a``); and varying only the volume but not the +cell shape (``AtomMove.OptCell.Constraint volume``). Fixing the ratio between two +cell lengths does not determine the minimisation fully: we choose to maintain the +*average* stress in the two directions as well as the ratio. + Go to :ref:`top `. .. _sr_both: @@ -147,7 +156,10 @@ allows *orthorhombic* unit cells). This can be done by setting AtomMove.EnthalpyTolerance 1E-5 AtomMove.StressTolerance 0.1 -Note that stress is in GPa and enthalpy is in Ha by default. +Note that stress is in GPa and enthalpy is in Ha by default. It is possible +to apply constraints to the simulation cell as described above, but this will +require extra care from the user to ensure that the simulation proceeds as +desired. The enthalpy will generally converge much more rapidly than the force and stress, and that it may be necessary to tighten ``minE.SCTolerance``