Feature request from Intel: constant volume variable cell relaxation

(cell_dofree = 'shape'). (D.C.)


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@6815 c92efa57-630b-4861-b058-cf58834340f0

Modules/cell_base.f90

@@ -635,7 +635,7 @@ END FUNCTION saw
      thdiag = .false.
      SELECT CASE ( TRIM( cell_dofree ) )
 
-            CASE ( 'all', 'default' )
+            CASE ( 'all', 'default', 'shape' )
               iforceh = 1
             CASE ( 'volume' )
               CALL errore(' init_dofree ', &

PW/move_ions.f90

@@ -32,6 +32,7 @@ SUBROUTINE move_ions()
   USE control_flags,          ONLY : istep, nstep, upscale, lbfgs, ldamped, &
                                      lconstrain, lcoarsegrained, conv_ions, &
                                      lmd, llang, history, tr2
+  USE input_parameters,       ONLY : cell_dofree
   USE relax,                  ONLY : epse, epsf, epsp, starting_scf_threshold
   USE lsda_mod,               ONLY : lsda, absmag
   USE metadyn_base,           ONLY : set_target, mean_force
@@ -54,6 +55,7 @@ SUBROUTINE move_ions()
   REAL(DP), ALLOCATABLE :: pos(:), grad(:)
   REAL(DP)              :: h(3,3), fcell(3,3)=0.d0, epsp1
   INTEGER,  ALLOCATABLE :: fixion(:)
+  real(dp) :: tr
   !
   !
   ! ... only one node does the calculation in the parallel case
@@ -144,6 +146,7 @@ SUBROUTINE move_ions()
         !
         IF ( lmovecell ) THEN
            ! changes needed only if cell moves
+           if (cell_dofree == 'shape') call impose_deviatoric_strain(alat*at, h)
            at = h /alat  
            CALL recips( at(1,1),at(1,2),at(1,3), bg(1,1),bg(1,2),bg(1,3) )
            CALL volume( alat, at(1,1), at(1,2), at(1,3), omega )

PW/stress.f90

@@ -29,6 +29,7 @@ subroutine stress
   USE bp,            ONLY : lelfield
   USE uspp,          ONLY : okvan
   USE london_module, ONLY : stres_london
+  USE input_parameters,       ONLY : cell_dofree
 #ifdef EXX
   USE exx,           ONLY : exx_stress
   USE funct,         ONLY : dft_is_hybrid, exx_is_active
@@ -139,11 +140,17 @@ subroutine stress
   ! Resymmetrize the total stress. This should not be strictly necessary,
   ! but prevents loss of symmetry in long vc-bfgs runs
 
-   CALL symmatrix ( sigma )
+  CALL symmatrix ( sigma )
+
+  if (cell_dofree == 'shape') then
+      WRITE(stdout,9001) (sigma(1,1) + sigma(2,2) + sigma(3,3)) * uakbar / 3d0
+      WRITE(stdout,*)
+      call impose_deviatoric_stress(sigma)
+  endif
+
   !
   ! write results in Ryd/(a.u.)^3 and in kbar
   !
-
   WRITE( stdout, 9000) (sigma(1,1) + sigma(2,2) + sigma(3,3)) * uakbar / 3d0,  &
                   (sigma(l,1), sigma(l,2), sigma(l,3),                    &
             sigma(l,1)*uakbar, sigma(l,2)*uakbar, sigma(l,3)*uakbar, l=1,3)
@@ -176,6 +183,7 @@ subroutine stress
   return
 9000 format (10x,'total   stress  (Ry/bohr**3) ',18x,'(kbar)', &
              &5x,'P=',f8.2/3 (3f13.8,4x,3f10.2/))
+9001 format (5x,'Isostatic pressure: ',f8.2,' kbar')
 9005 format &
          &  (5x,'kinetic stress (kbar)',3f10.2/2(26x,3f10.2/)/ &
          &   5x,'local   stress (kbar)',3f10.2/2(26x,3f10.2/)/ &

PW/vcsmd.f90

@@ -41,6 +41,8 @@ SUBROUTINE vcsmd()
   USE parameters,      ONLY : ntypx
   USE ener,            ONLY : etot
   USE io_files,        ONLY : prefix, delete_if_present
+  USE input_parameters,       ONLY : cell_dofree
+
   !
   IMPLICIT NONE
   !
@@ -412,6 +414,7 @@ SUBROUTINE vcsmd()
   !
   ! ... update configuration in PWSCF variables
   !
+  if (cell_dofree == 'shape') call impose_deviatoric_strain(alat*at, avec)
   at = avec / alat
   !
   CALL volume( alat, at(1,1), at(1,2), at(1,3), omega )

flib/Makefile

@@ -38,7 +38,8 @@ transto.o \
 date_and_tim.o \
 volume.o \
 dylmr2.o \
-ylmr2.o
+ylmr2.o \
+deviatoric.o
 
 POBJS = \
 flush_unit.o \

flib/deviatoric.f90

@@ -0,0 +1,54 @@
+!
+! Copyright (C) 2001 PWSCF group
+! This file is distributed under the terms of the
+! GNU General Public License. See the file `License'
+! in the root directory of the present distribution,
+! or http://www.gnu.org/copyleft/gpl.txt .
+!
+!
+!---------------------------------------------------------------------
+subroutine impose_deviatoric_strain ( at_old, at )
+  !---------------------------------------------------------------------
+  !
+  !     Impose a pure deviatoric deformation
+  !
+  use kinds, ONLY: dp
+  implicit none
+  real(dp), intent(in)    :: at_old(3,3)
+  real(dp), intent(inout) :: at(3,3)
+  real(dp) :: tr, omega, omega_old
+
+  tr = (at(1,1)+at(2,2)+at(3,3))/3.d0
+  tr = tr - (at_old(1,1)+at_old(2,2)+at_old(3,3))/3.d0
+
+  at(1,1) = at(1,1) - tr
+  at(2,2) = at(2,2) - tr
+  at(3,3) = at(3,3) - tr
+
+  call volume (1.d0, at_old(1,1), at_old(1,2), at_old(1,3), omega_old)
+  call volume (1.d0, at(1,1), at(1,2), at(1,3), omega)
+  at = at * (omega_old / omega)**(1.d0/3.d0)
+
+end subroutine impose_deviatoric_strain
+
+
+
+!---------------------------------------------------------------------
+subroutine impose_deviatoric_stress ( sigma )
+  !---------------------------------------------------------------------
+  !
+  !     Impose a pure deviatoric stress
+  !
+  use kinds, ONLY: dp
+  implicit none
+  real(dp), intent(inout) :: sigma(3,3)
+  real(dp) :: tr
+
+  tr = (sigma(1,1)+sigma(2,2)+sigma(3,3))/3.d0
+  sigma(1,1) = sigma(1,1) - tr
+  sigma(2,2) = sigma(2,2) - tr
+  sigma(3,3) = sigma(3,3) - tr
+
+end subroutine impose_deviatoric_stress
+
+