Skip to content

Commit

Permalink
Fix CELL_OPT ignoring CONSTRAINT for convergence criterion
Browse files Browse the repository at this point in the history
The internal pressure is calculated including dimensions fixed
with the CONSTRAINT keyword. This can result in CELL_OPT
calculations that never converge because the direction with
remaining stress is not allowed to change. With this patch
only directions that are not constrained contribute to the
internal pressure.
  • Loading branch information
ducryf authored and mkrack committed Jun 30, 2020
1 parent 3b42c2a commit c405d2f
Showing 1 changed file with 17 additions and 5 deletions.
22 changes: 17 additions & 5 deletions src/motion/cell_opt_utils.F
Original file line number Diff line number Diff line change
Expand Up @@ -272,7 +272,7 @@ SUBROUTINE get_dg_dh(gradient, av_ptens, pres_ext, cell, mtrx, keep_angles, &

CHARACTER(LEN=*), PARAMETER :: routineN = 'get_dg_dh', routineP = moduleN//':'//routineN

INTEGER :: i, my_constraint_id
INTEGER :: my_constraint_id
LOGICAL :: my_keep_angles, my_keep_symmetry
REAL(KIND=dp), DIMENSION(3, 3) :: correction, pten_hinv_old, ptens

Expand All @@ -292,10 +292,22 @@ SUBROUTINE get_dg_dh(gradient, av_ptens, pres_ext, cell, mtrx, keep_angles, &
ptens = av_ptens

! Evaluating the internal pressure
pres_int = 0.0_dp
DO i = 1, 3
pres_int = pres_int + ptens(i, i)
ENDDO
SELECT CASE (my_constraint_id)
CASE (fix_x)
pres_int = ptens(2, 2) + ptens(3, 3)
CASE (fix_y)
pres_int = ptens(1, 1) + ptens(3, 3)
CASE (fix_z)
pres_int = ptens(1, 1) + ptens(2, 2)
CASE (fix_xy)
pres_int = ptens(3, 3)
CASE (fix_xz)
pres_int = ptens(2, 2)
CASE (fix_yz)
pres_int = ptens(1, 1)
CASE (fix_none)
pres_int = ptens(1, 1) + ptens(2, 2) + ptens(3, 3)
END SELECT
pres_int = pres_int/3.0_dp

ptens(1, 1) = av_ptens(1, 1) - pres_ext
Expand Down

0 comments on commit c405d2f

Please sign in to comment.