Merge pull request #395 from OrderN/f-cell-opt-constraint

Updates to cell constraints and how they are applied

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

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 <strucrelax>`.
 
 .. _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``

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

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)

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
   !!***