Add k-point grid to mbd

Modules/mbdlib.f90

@@ -17,8 +17,8 @@ MODULE libmbd_interface
   USE kinds,            ONLY : DP
   USE io_global,        ONLY : stdout 
   USE tsvdw_module,     ONLY : veff_pub, vfree_pub, vdw_isolated
-  USE ions_base,        ONLY : nat, atm, zv, tau, taui, ntyp => nsp, ityp
-  USE cell_base,        ONLY : alat, at, bg, omega, ainv
+  USE ions_base,        ONLY : nat, atm, tau, ityp
+  USE cell_base,        ONLY : alat, at, ainv
   USE funct,            ONLY : get_dft_short
   USE control_flags,    ONLY : conv_elec
   USE constants,        ONLY : ry_kbar
@@ -35,6 +35,7 @@ MODULE libmbd_interface
 
 
   INTEGER:: na
+  LOGICAL:: do_gradients
   TYPE(mbd_input_t):: inp
   TYPE(mbd_calc_t):: calc
   REAL(dp), DIMENSION(:), ALLOCATABLE:: ratios
@@ -47,15 +48,25 @@ MODULE libmbd_interface
 !#############################################################
 ! This subroutine sets up the library before the first call
 !#############################################################
-  SUBROUTINE init_mbd()
-  IMPLICIT NONE
-
+SUBROUTINE init_mbd ( nks_start, nk1, nk2, nk3, k1, k2, k3 )
+  !
+  INTEGER, INTENT(IN) :: nks_start, nk1, nk2, nk3, k1, k2, k3
   !
   ! Allocation of variables that depend on the number of atoms
   !
   ALLOCATE(inp%atom_types(nat))
-  IF(.NOT.ALLOCATED(mbd_gradient)) ALLOCATE(mbd_gradient(3, nat))
-  IF(tprnfor .OR. tstress .AND. .NOT.ALLOCATED(FmbdvdW)) ALLOCATE(FmbdvdW(3, nat))
+  !
+  EmbdvdW  = 0.0_dp
+  !
+  do_gradients = tprnfor .OR. tstress
+  IF ( do_gradients ) THEN
+     !
+     IF(.NOT.ALLOCATED(mbd_gradient)) ALLOCATE(mbd_gradient(3, nat))
+     !
+     IF(.NOT.ALLOCATED(FmbdvdW)) ALLOCATE(FmbdvdW(3, nat))
+     !
+  END IF
+  !
   ALLOCATE(ratios(nat))
 
   inp%log_level=1
@@ -66,16 +77,31 @@ SUBROUTINE init_mbd()
     inp%atom_types(na) = trim(atm(ityp(na)))
   ENDDO
   inp%coords = tau*alat  ! HK-TODO: this one works for PW (check if it is for CP)
-
   !
   ! If we pass lattice vectors to the library, it uses the algorithm for
   ! periodic system automatically
   !
   IF( .NOT.vdw_isolated ) THEN
     inp%lattice_vectors = at*alat  ! Lattice vector in real space
-    inp%k_grid = [1, 1, 1] !the k points grid is not needed to be the same as for the PW calculation, but it would help the convergence (TODO)
+    !
+    IF ( nks_start == 0 ) THEN
+      ! K-point mesh
+      inp%k_grid = [nk1, nk2, nk3]
+      !
+      inp%k_grid_shift = 0.5_DP
+      IF (k1 .EQ. k2 .AND. k2 .EQ. k3 .AND. k3 .EQ. 0) &
+        CALL infomsg('mbdlib','k-point shift ignored')
+      !
+    ELSE
+      inp%k_grid = [1, 1, 1] !set default k points grid
+      inp%k_grid_shift = 0.5_DP ! set default shift
+    ENDIF
+    !
   ENDIF
-
+  !
+  WRITE(stdout, '(5x,"mbdlib: K-point grid set to ",3I3,", shift: ",F4.2)') &
+        inp%k_grid, inp%k_grid_shift
+  !
   select case (TRIM(get_dft_short()))  ! An empirical factor needs to be set based on the functiona
   CASE ('PBE')
     inp%xc = 'pbe'
@@ -125,15 +151,15 @@ SUBROUTINE mbd_interface()
   ENDIF
 
   CALL calc%evaluate_vdw_method(EmbdvdW)  !MBD energy
-  IF ( tprnfor .OR. tstress) THEN
-    CALL calc%get_gradients(mbd_gradient)  
+  IF ( do_gradients ) THEN
+    CALL calc%get_gradients(mbd_gradient)
     FmbdvdW = -mbd_gradient  ! Ionic forces with correct sign
   ENDIF
   !
 
-  IF( tprnfor .OR. tstress .AND. .NOT.vdw_isolated ) THEN
+  IF( do_gradients .AND. .NOT.vdw_isolated ) THEN
     CALL calc%get_lattice_stress(cell_derivs)
-    HmbdvdW=MATMUL(cell_derivs, TRANSPOSE(ainv)) 
+    HmbdvdW=MATMUL(cell_derivs, TRANSPOSE(ainv))
   ENDIF
 
   RETURN

PW/src/init_run.f90

@@ -10,6 +10,7 @@ SUBROUTINE init_run()
   !----------------------------------------------------------------------------
   !
   USE klist,              ONLY : nkstot
+  USE start_k,            ONLY : nks_start, nk1, nk2, nk3, k1, k2, k3
   USE symme,              ONLY : sym_rho_init
   USE wvfct,              ONLY : nbnd, et, wg, btype
   USE control_flags,      ONLY : lmd, gamma_only, smallmem, ts_vdw, mbd_vdw
@@ -126,7 +127,7 @@ SUBROUTINE init_run()
      CALL set_h_ainv()
   END IF
   IF (mbd_vdw) THEN
-     CALL init_mbd()
+     CALL init_mbd( nks_start, nk1, nk2, nk3, k1, k2, k3 )
   END IF
   !
   CALL allocate_wfc_k()