SIRIUS: Update to version 6.4.1, add stress tensor, SpFFT, and more t…

…ests.

src/force_env_methods.F

@@ -260,7 +260,11 @@ RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &
          CASE (use_qs_force)
             CALL qs_calc_energy_force(force_env%qs_env, calculate_forces, energy_consistency, linres_run)
          CASE (use_pwdft_force)
-            CALL pwdft_calc_energy_force(force_env%pwdft_env, calculate_forces)
+            IF (virial%pv_availability .AND. calculate_stress_tensor) THEN
+               CALL pwdft_calc_energy_force(force_env%pwdft_env, calculate_forces,.NOT. virial%pv_numer)
+            ELSE
+               CALL pwdft_calc_energy_force(force_env%pwdft_env, calculate_forces, .FALSE.)
+            END IF
          CASE (use_eip_force)
             IF (force_env%eip_env%eip_model == use_lenosky_eip) THEN
                CALL eip_lenosky(force_env%eip_env)

src/pwdft_environment.F

@@ -46,6 +46,7 @@ MODULE pwdft_environment
    USE sirius_interface,                ONLY: cp_sirius_create_env,&
                                               cp_sirius_energy_force,&
                                               cp_sirius_update_context
+   USE virial_types,                    ONLY: virial_type
 #include "./base/base_uses.f90"
 
    IMPLICIT NONE
@@ -185,14 +186,14 @@ SUBROUTINE pwdft_init(pwdft_env, root_section, para_env, force_env_section, subs
       CALL write_symmetry(particle_set, my_cell, subsys_section)
 
       IF (iw > 0) THEN
-         WRITE (iw, '(A,A,A)') " ==================================", " SIRIUS INIT", &
+         WRITE (iw, '(A,A,A)') " ==================================", " SIRIUS INIT ", &
             "================================="
       END IF
       ! Sirius initialization
       CALL cp_sirius_create_env(pwdft_env)
 
       IF (iw > 0) THEN
-         WRITE (iw, '(A,A,A)') " ==================================", " SIRIUS INIT FINISHED ", &
+         WRITE (iw, '(A,A,A)') " =========================", " SIRIUS INIT FINISHED ", &
             "================================="
       END IF
       IF (iw > 0) CALL m_flush(iw)
@@ -205,22 +206,26 @@ END SUBROUTINE pwdft_init
 !> \brief Calculate energy and forces within the PWDFT/SIRIUS code
 !> \param pwdft_env The pwdft environment to retain
 !> \param calculate_forces ...
+!> \param calculate_stress ...
 !> \par History
 !>      03.2018 initial create
 !> \author JHU
 ! **************************************************************************************************
-   SUBROUTINE pwdft_calc_energy_force(pwdft_env, calculate_forces)
+   SUBROUTINE pwdft_calc_energy_force(pwdft_env, calculate_forces, calculate_stress)
       TYPE(pwdft_environment_type), POINTER              :: pwdft_env
-      LOGICAL, INTENT(IN)                                :: calculate_forces
+      LOGICAL, INTENT(IN)                                :: calculate_forces, calculate_stress
 
       CHARACTER(len=*), PARAMETER :: routineN = 'pwdft_calc_energy_force', &
          routineP = moduleN//':'//routineN
 
       INTEGER                                            :: handle, iatom, iw, natom
+      REAL(KIND=dp), DIMENSION(1:3, 1:3)                 :: stress
       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: force
+      TYPE(cell_type), POINTER                           :: my_cell
       TYPE(cp_logger_type), POINTER                      :: logger
       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
       TYPE(qs_subsys_type), POINTER                      :: qs_subsys
+      TYPE(virial_type), POINTER                         :: virial
 
       CALL timeset(routineN, handle)
 
@@ -239,13 +244,13 @@ SUBROUTINE pwdft_calc_energy_force(pwdft_env, calculate_forces)
       CALL cp_sirius_update_context(pwdft_env)
 
       IF (iw > 0) THEN
-         WRITE (iw, '(A,A,A)') " =============================", " SIRIUS UPDATE CONTEXT FINISHED", &
+         WRITE (iw, '(A,A,A)') " ====================", " SIRIUS UPDATE CONTEXT FINISHED", &
             "============================"
       END IF
       IF (iw > 0) CALL m_flush(iw)
 
       ! calculate energy and forces/stress
-      CALL cp_sirius_energy_force(pwdft_env, calculate_forces)
+      CALL cp_sirius_energy_force(pwdft_env, calculate_forces, calculate_stress)
 
       IF (calculate_forces) THEN
          CALL pwdft_env_get(pwdft_env=pwdft_env, qs_subsys=qs_subsys)
@@ -256,6 +261,13 @@ SUBROUTINE pwdft_calc_energy_force(pwdft_env, calculate_forces)
          END DO
       END IF
 
+      IF (calculate_stress) THEN
+         ! i need to retrieve the volume of the unit cell for the stress tensor
+         CALL qs_subsys_get(qs_subsys, cell=my_cell, virial=virial)
+         CALL pwdft_env_get(pwdft_env=pwdft_env, stress=stress)
+         virial%pv_virial(1:3, 1:3) = -stress(1:3, 1:3)*my_cell%deth
+      END IF
+
       CALL timestop(handle)
 
    END SUBROUTINE pwdft_calc_energy_force

src/sirius_interface.F

@@ -190,12 +190,12 @@ SUBROUTINE cp_sirius_create_env(pwdft_env)
       pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "parameters")
       IF (ASSOCIATED(pwdft_sub_section)) THEN
          CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, string("parameters"))
-         CALL section_vals_val_get(pwdft_sub_section, "k_grid", i_vals=kk)
+         CALL section_vals_val_get(pwdft_sub_section, "ngridk", i_vals=kk)
          k_grid(1) = kk(1)
          k_grid(2) = kk(2)
          k_grid(3) = kk(3)
 
-         CALL section_vals_val_get(pwdft_sub_section, "k_shift", i_vals=kk)
+         CALL section_vals_val_get(pwdft_sub_section, "shiftk", i_vals=kk)
          k_shift(1) = kk(1)
          k_shift(2) = kk(2)
          k_shift(3) = kk(3)
@@ -557,15 +557,16 @@ END SUBROUTINE cp_sirius_fill_in_section
 !> \brief ...
 !> \param pwdft_env ...
 !> \param calculate_forces ...
+!> \param calculate_stress_tensor ...
 !> \param
 !> \par History
 !>      07.2018 start the Sirius library
 !> \author JHU
 ! **************************************************************************************************
-   SUBROUTINE cp_sirius_energy_force(pwdft_env, calculate_forces)
+   SUBROUTINE cp_sirius_energy_force(pwdft_env, calculate_forces, calculate_stress_tensor)
       TYPE(pwdft_environment_type), INTENT(INOUT), &
          POINTER                                         :: pwdft_env
-      LOGICAL, INTENT(IN)                                :: calculate_forces
+      LOGICAL, INTENT(IN)                                :: calculate_forces, calculate_stress_tensor
 
       INTEGER                                            :: n1, n2
       REAL(KIND=C_DOUBLE)                                :: etotal
@@ -603,12 +604,15 @@ SUBROUTINE cp_sirius_energy_force(pwdft_env, calculate_forces)
          ! cp2k : forces(atoms, coordinates)
          forces = -TRANSPOSE(cforces(:, :))
          DEALLOCATE (cforces)
+      ENDIF
 
+      IF (calculate_stress_tensor) THEN
          cstress = 0.0_C_DOUBLE
          CALL sirius_get_stress_tensor(gs_handler, string('total'), cstress(1, 1))
          stress(1:3, 1:3) = cstress(1:3, 1:3)
          CALL pwdft_env_set(pwdft_env=pwdft_env, stress=stress)
       ENDIF
+
    END SUBROUTINE cp_sirius_energy_force
 
 END MODULE sirius_interface
@@ -657,13 +661,15 @@ END SUBROUTINE cp_sirius_create_env
 !> \brief Empty implementation in case SIRIUS is not compiled in.
 !> \param pwdft_env ...
 !> \param calculate_forces ...
+!> \param calculate_stress ...
 ! **************************************************************************************************
-   SUBROUTINE cp_sirius_energy_force(pwdft_env, calculate_forces)
+   SUBROUTINE cp_sirius_energy_force(pwdft_env, calculate_forces, calculate_stress)
       TYPE(pwdft_environment_type), POINTER              :: pwdft_env
-      LOGICAL                                            :: calculate_forces
+      LOGICAL                                            :: calculate_forces, calculate_stress
 
       MARK_USED(pwdft_env)
       MARK_USED(calculate_forces)
+      MARK_USED(calculate_stress)
       CPABORT("Sirius library is missing")
    END SUBROUTINE cp_sirius_energy_force
 

tests/SIRIUS/regtest-1/Au.inp

@@ -11,15 +11,15 @@
     &PARAMETERS
        ELECTRONIC_STRUCTURE_METHOD  pseudopotential
        SMEARING_WIDTH 0.01
-       USE_SYMMETRY true
+       USE_SYMMETRY false
        NUM_MAG_DIMS 3
        SPIN_ORBIT true
        GK_CUTOFF 6.0
        PW_CUTOFF 20.00
        ENERGY_TOL 1e-10
        POTENTIAL_TOL 1e-8
        NUM_DFT_ITER 100
-       K_GRID 2 2 2
+       NGRIDK 2 2 2
     &END PARAMETERS
     &ITERATIVE_SOLVER
        ENERGY_TOLERANCE 1e-5
@@ -31,7 +31,8 @@
     &MIXER
        TYPE broyden1
        MAX_HISTORY 8
-    &END MIXER
+       BETA 0.6
+  &END MIXER
   &END PW_DFT
     &DFT
       &XC
@@ -48,9 +49,9 @@
 
   &SUBSYS
     &CELL
-      A -3.85500000000000000000 0.0 3.85500000000000000000
-      B 0.0 3.85500000000000000000 3.85500000000000000000
-      C -3.85500000000000000000 3.85500000000000000000 0.0
+      A [bohr] -3.85500000000000000000 0.0 3.85500000000000000000
+      B [bohr] 0.0 3.85500000000000000000 3.85500000000000000000
+      C [bohr] -3.85500000000000000000 3.85500000000000000000 0.0
     &END CELL
     &COORD
       SCALED

tests/SIRIUS/regtest-1/Au_GTH.inp

@@ -19,7 +19,7 @@
        ENERGY_TOL 1e-10
        POTENTIAL_TOL 1e-8
        NUM_DFT_ITER 100
-       K_GRID 2 2 2
+       NGRIDK 2 2 2
       &END PARAMETERS
      &ITERATIVE_SOLVER
       RESIDUAL_TOLERANCE 1e-5
@@ -47,9 +47,9 @@
   &END DFT
   &SUBSYS
     &CELL
-      A -3.85500000000000000000 0.0 3.85500000000000000000
-      B 0.0 3.85500000000000000000 3.85500000000000000000
-      C -3.85500000000000000000 3.85500000000000000000 0.0
+      A [bohr] -3.85500000000000000000 0.0 3.85500000000000000000
+      B [bohr] 0.0 3.85500000000000000000 3.85500000000000000000
+      C [bohr] -3.85500000000000000000 3.85500000000000000000 0.0
     &END CELL
     &COORD
       SCALED

tests/SIRIUS/regtest-1/Fe.inp

@@ -18,7 +18,7 @@
       PW_CUTOFF 25.00
       ENERGY_TOL 1e-8
       NUM_DFT_ITER 100
-      K_GRID 2 2 2
+      NGRIDK 2 2 2
     &END PARAMETERS
     &ITERATIVE_SOLVER
       ENERGY_TOLERANCE 1e-2
@@ -48,7 +48,7 @@
   &END DFT
   &SUBSYS
     &CELL
-      ABC 5.5415878223 5.5415878223 5.5415878223
+      ABC [bohr] 5.5415878223 5.5415878223 5.5415878223
     &END CELL
     &COORD
       SCALED

tests/SIRIUS/regtest-1/He-full-potential.inp

@@ -0,0 +1,74 @@
+&FORCE_EVAL
+  METHOD SIRIUS
+  &PW_DFT
+    &CONTROL
+       PROCESSING_UNIT cpu
+       STD_EVP_SOLVER_TYPE lapack
+       GEN_EVP_SOLVER_TYPE lapack
+       VERBOSITY 2
+    &END CONTROL
+    &PARAMETERS
+       ELECTRONIC_STRUCTURE_METHOD  full_potential_lapwlo
+       SMEARING_WIDTH 0.05
+       USE_SYMMETRY true
+       NUM_MAG_DIMS 0
+       NUM_FV_STATES 40
+       AW_CUTOFF 8.0
+       PW_CUTOFF 20.00
+       AUTO_RMT 0
+       ENERGY_TOL 1e-7
+       POTENTIAL_TOL 1e-7
+       NUM_DFT_ITER 20
+       NGRIDK 1 1 1
+       MOLECULE true
+       lmax_apw 8
+       lmax_pot 8
+       lmax_rho 8
+    &END PARAMETERS
+    &ITERATIVE_SOLVER
+       ENERGY_TOLERANCE 1e-4
+       RESIDUAL_TOLERANCE 1e-5
+       NUM_STEPS 8
+       SUBSPACE_SIZE 8
+       TYPE davidson
+       CONVERGE_BY_ENERGY 0
+    &END ITERATIVE_SOLVER
+    &MIXER
+       BETA 0.95
+       TYPE broyden1
+       MAX_HISTORY 8
+    &END MIXER
+  &END PW_DFT
+    &DFT
+      &XC
+         &XC_FUNCTIONAL
+            &LIBXC
+               FUNCTIONAL XC_LDA_X
+            &END LIBXC
+            &LIBXC
+               FUNCTIONAL XC_LDA_C_PZ
+            &END LIBXC
+         &END XC_FUNCTIONAL
+      &END XC
+    &END DFT
+
+  &SUBSYS
+    &CELL
+      A [bohr] 10.0 0.0 0.0
+      B [bohr] 0 10.0 0.0
+      C [bohr] 0.0 0.0 10.0
+    &END CELL
+    &COORD
+      SCALED
+        He 0.0 0.0 0.0
+    &END COORD
+    &KIND He
+      POTENTIAL UPF "He.json"
+    &END KIND
+  &END SUBSYS
+&END FORCE_EVAL
+&GLOBAL
+  PROJECT He-full-potential
+  PRINT_LEVEL MEDIUM
+  RUN_TYPE ENERGY
+&END GLOBAL