CP2K+SMEAGOL interface for electron transport with non-equilibrium Gr…

…een’s functions Co-authored-by: Sergey Chulkov <[email protected]> Co-authored-by: Clotilde Cucinotta <[email protected]>
hfp · Oct 21, 2024 · a1c7dd2 · a1c7dd2
1 parent 8e29cc6
commit a1c7dd2
Show file tree

Hide file tree

Showing 34 changed files with 5,641 additions and 31 deletions.
diff --git a/INSTALL.md b/INSTALL.md
@@ -466,7 +466,7 @@ Calls to `offload_dgemm` also accept pointers on GPU or a combination of them.
 
 ### 2y. DeePMD-kit (optional, wider range of interaction potentials)
 
-DeePMD-kit - Deep Potential Molecular Dyanmics. Support for DeePMD-kit can be enabled via the flag
+DeePMD-kit - Deep Potential Molecular Dynamics. Support for DeePMD-kit can be enabled via the flag
 `-D__DEEPMD`.
 
 - DeePMD-kit C interface can be downloaded from
@@ -480,6 +480,17 @@ DeePMD-kit - Deep Potential Molecular Dyanmics. Support for DeePMD-kit can be en
 - Add `-D__DFTD4` to DFLAGS, `-ldftd4 -lmstore -lmulticharge -lmctc-lib` to LIBS and
   `-I'${DFTD4_DFTD4}/../..' -I'${DFTD4_DFTD4}' -I'${DFTD4_MCTC}'` to CFLAGS
 
+### 2y. libsmeagol (optional, electron transport calculation with current-induced forces)
+
+- libsmeagol is an external library to compute electron transport properties using Non-Equilibrium
+  Green Functions (NEGF) method. The library can be downloaded from
+  <https://github.com/StefanoSanvitoGroup/libsmeagol>.
+- libsmeagol depends on an MPI library and can only be linked with MPI parallel CP2K binaries.
+- During the installation, the directories `$(LIBSMEAGOL_DIR)/lib` and `$(LIBGRPP_DIR)/obj` are
+  created.
+- Add `-D__SMEAGOL` to DFLAGS, `-I$(LIBSMEAGOL_DIR)/obj` to FCFLAGS and
+  `-L$(LIBSMEAGOL_DIR)/lib -lsmeagol` to LIBS
+
 ## 3. Compile
 
 ### 3a. ARCH files

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
@@ -286,6 +286,7 @@ list(
   input_cp2k_restarts_util.F
   input_cp2k_rsgrid.F
   input_cp2k_subsys.F
+  input_cp2k_smeagol.F
   input_cp2k_tb.F
   input_cp2k_thermostats.F
   input_cp2k_transport.F
@@ -777,6 +778,10 @@ list(
   semi_empirical_utils.F
   simpar_types.F
   sirius_interface.F
+  smeagol_control_types.F
+  smeagol_emtoptions.F
+  smeagol_interface.F
+  smeagol_matrix_utils.F
   spglib_f08.F
   splines_methods.F
   splines_types.F

diff --git a/src/common/bibliography.F b/src/common/bibliography.F
@@ -88,7 +88,7 @@ MODULE bibliography
                             Caldeweyher2017, Caldeweyher2019, Caldeweyher2020, Freeman1977, Gruneis2009, &
                             Stein2022, Stein2024, &
                             Blase2018, Blase2020, Bruneval2015, Golze2019, Gui2018, Jacquemin2017, Liu2020, &
-                            Sander2015, Schreiber2008, vanSetten2015, Setyawan2010
+                            Sander2015, Schreiber2008, vanSetten2015, Setyawan2010, Ahart2024
 
 CONTAINS
 
@@ -1877,6 +1877,13 @@ SUBROUTINE add_all_references()
                          source="Comput. Mater. Sci.", volume="49", pages="299-312", &
                          year=2010, doi="10.1016/j.commatsci.2010.05.010")
 
+      CALL add_reference(key=Ahart2024, &
+                         authors=s2a("C. Ahart", "S. Chulkov", "C. Cucinotta"), &
+                         title="Enabling Ab Initio Molecular Dynamics under Bias: The CP2K+SMEAGOL Interface "// &
+                         "for Integrating Density Functional Theory and Non-Equilibrium Green Functions", &
+                         source="J. Chem. Theory Comput.", volume="20", pages="6772-6780", &
+                         year=2024, doi="10.1021/acs.jctc.4c00371")
+
    END SUBROUTINE add_all_references
 
 END MODULE bibliography
diff --git a/src/cp2k_info.F b/src/cp2k_info.F
@@ -271,6 +271,10 @@ FUNCTION cp2k_flags() RESULT(flags)
       flags = TRIM(flags)//" offload_unified_memory"
 #endif
 
+#if defined(__SMEAGOL)
+      flags = TRIM(flags)//" libsmeagol"
+#endif
+
    END FUNCTION cp2k_flags
 
 ! **************************************************************************************************

diff --git a/src/cp_control_types.F b/src/cp_control_types.F
@@ -23,6 +23,9 @@ MODULE cp_control_types
    USE qs_cdft_types,                   ONLY: cdft_control_create,&
                                               cdft_control_release,&
                                               cdft_control_type
+   USE smeagol_control_types,           ONLY: smeagol_control_create,&
+                                              smeagol_control_release,&
+                                              smeagol_control_type
    USE xas_control,                     ONLY: xas_control_release,&
                                               xas_control_type
 #include "./base/base_uses.f90"
@@ -571,6 +574,7 @@ MODULE cp_control_types
       TYPE(xas_control_type), POINTER      :: xas_control => NULL()
       TYPE(expot_control_type), POINTER    :: expot_control => NULL()
       TYPE(maxwell_control_type), POINTER  :: maxwell_control => NULL()
+      TYPE(smeagol_control_type), POINTER  :: smeagol_control => NULL()
       TYPE(efield_p_type), POINTER, &
          DIMENSION(:)                      :: efield_fields => NULL()
       INTEGER                              :: nspins = 0, &
@@ -779,13 +783,15 @@ SUBROUTINE dft_control_create(dft_control)
       NULLIFY (dft_control%admm_control)
       NULLIFY (dft_control%expot_control)
       NULLIFY (dft_control%maxwell_control)
+      NULLIFY (dft_control%smeagol_control)
       NULLIFY (dft_control%rtp_control)
       NULLIFY (dft_control%sccs_control)
       dft_control%do_sccs = .FALSE.
       dft_control%apply_embed_pot = .FALSE.
       dft_control%apply_dmfet_pot = .FALSE.
       CALL qs_control_create(dft_control%qs_control)
       CALL tddfpt2_control_create(dft_control%tddfpt2_control)
+      CALL smeagol_control_create(dft_control%smeagol_control)
    END SUBROUTINE dft_control_create
 
 ! **************************************************************************************************
@@ -808,6 +814,7 @@ SUBROUTINE dft_control_release(dft_control)
       CALL admm_control_release(dft_control%admm_control)
       CALL expot_control_release(dft_control%expot_control)
       CALL maxwell_control_release(dft_control%maxwell_control)
+      CALL smeagol_control_release(dft_control%smeagol_control)
       CALL efield_fields_release(dft_control%efield_fields)
       IF (ASSOCIATED(dft_control%sccs_control)) DEALLOCATE (dft_control%sccs_control)
       IF (ASSOCIATED(dft_control%period_efield)) THEN

diff --git a/src/cp_control_utils.F b/src/cp_control_utils.F
@@ -67,6 +67,7 @@ MODULE cp_control_utils
    USE pair_potential_types,            ONLY: pair_potential_reallocate
    USE periodic_table,                  ONLY: get_ptable_info
    USE qs_cdft_utils,                   ONLY: read_cdft_control_section
+   USE smeagol_control_types,           ONLY: read_smeagol_control
    USE string_utilities,                ONLY: uppercase
    USE util,                            ONLY: sort
    USE xc,                              ONLY: xc_uses_kinetic_energy_density,&
@@ -614,6 +615,10 @@ SUBROUTINE read_dft_control(dft_control, dft_section)
       tmp_section => section_vals_get_subs_vals(dft_section, "EXTERNAL_VXC")
       CALL section_vals_get(tmp_section, explicit=dft_control%apply_external_vxc)
 
+      ! SMEAGOL interface
+      tmp_section => section_vals_get_subs_vals(dft_section, "SMEAGOL")
+      CALL read_smeagol_control(dft_control%smeagol_control, tmp_section)
+
    END SUBROUTINE read_dft_control
 
 ! **************************************************************************************************

diff --git a/src/input/cp_output_handling.F b/src/input/cp_output_handling.F
@@ -62,7 +62,7 @@ MODULE cp_output_handling
 
    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_output_handling'
-   PUBLIC :: cp_print_key_should_output, cp_iterate, cp_add_iter_level, cp_rm_iter_level
+   PUBLIC :: cp_print_key_should_output, cp_iterate, cp_get_iter_level_by_name, cp_get_iter_nr, cp_add_iter_level, cp_rm_iter_level
    PUBLIC :: cp_iter_string, cp_print_key_section_create
    PUBLIC :: cp_print_key_unit_nr, cp_print_key_finished_output
    PUBLIC :: cp_print_key_generate_filename, cp_printkey_is_on
@@ -600,6 +600,56 @@ SUBROUTINE cp_iterate(iteration_info, last, iter_nr, increment, iter_nr_out)
       iteration_info%last_iter(iteration_info%n_rlevel) = my_last
    END SUBROUTINE cp_iterate
 
+! **************************************************************************************************
+!> \brief Return the index of an iteration level by its name.
+!> \param iteration_info the iteration info to query
+!> \param level_name     level name to query.
+!> \return iteration level index or 0 if there is no such level
+!> \author Sergey Chulkov
+! **************************************************************************************************
+   FUNCTION cp_get_iter_level_by_name(iteration_info, level_name) RESULT(rlevel)
+      TYPE(cp_iteration_info_type), INTENT(IN), POINTER  :: iteration_info
+      CHARACTER(LEN=*), INTENT(IN)                       :: level_name
+      INTEGER                                            :: rlevel
+
+      CPASSERT(ASSOCIATED(iteration_info))
+      CPASSERT(iteration_info%ref_count > 0)
+      DO rlevel = iteration_info%n_rlevel, 1, -1
+         IF (iteration_info%level_name(rlevel) == level_name) EXIT
+      END DO
+
+   END FUNCTION cp_get_iter_level_by_name
+
+! **************************************************************************************************
+!> \brief Return the current iteration number at a given level.
+!> \param iteration_info the iteration info to query
+!> \param rlevel         index of the iteration level. Use the level on top of the stack,
+!>                       if it is not given
+!> \param iter_nr        iteration number [out]
+!> \param last_iter      last iteration flag [out]
+!> \author Sergey Chulkov
+! **************************************************************************************************
+   SUBROUTINE cp_get_iter_nr(iteration_info, rlevel, iter_nr, last_iter)
+      TYPE(cp_iteration_info_type), INTENT(IN), POINTER  :: iteration_info
+      INTEGER, INTENT(IN), OPTIONAL                      :: rlevel
+      INTEGER, INTENT(OUT), OPTIONAL                     :: iter_nr
+      LOGICAL, INTENT(OUT), OPTIONAL                     :: last_iter
+
+      INTEGER                                            :: ilevel
+
+      CPASSERT(ASSOCIATED(iteration_info))
+      CPASSERT(iteration_info%ref_count > 0)
+      IF (PRESENT(rlevel)) THEN
+         CPASSERT(rlevel > 0 .AND. rlevel <= iteration_info%n_rlevel)
+         ilevel = rlevel
+      ELSE
+         ilevel = iteration_info%n_rlevel
+      END IF
+
+      IF (PRESENT(iter_nr)) iter_nr = iteration_info%iteration(ilevel)
+      IF (PRESENT(last_iter)) last_iter = iteration_info%last_iter(ilevel)
+   END SUBROUTINE cp_get_iter_nr
+
 ! **************************************************************************************************
 !> \brief Adds an iteration level
 !> \param iteration_info the iteration info to which an iteration level has

diff --git a/src/input_constants.F b/src/input_constants.F
@@ -1327,4 +1327,18 @@ MODULE input_constants
                                                sigma_PBE0_S2 = 3, &
                                                sigma_PBE_S1 = 4
 
+   ! SMEAGOL interface
+   INTEGER, PARAMETER, PUBLIC               :: smeagol_bulklead_leftright = 1, &
+                                               smeagol_bulklead_left = 2, &
+                                               smeagol_bulklead_right = 3
+
+   INTEGER, PARAMETER, PUBLIC               :: smeagol_gridmethod_traditional = 0, &
+                                               smeagol_gridmethod_adaptive = 1
+
+   INTEGER, PARAMETER, PUBLIC               :: smeagol_integraltype_gauss_legendre = 1, &
+                                               smeagol_integraltype_gauss_chebyshev = 2
+
+   INTEGER, PARAMETER, PUBLIC               :: smeagol_runtype_bulktransport = 0, &
+                                               smeagol_runtype_emtransport = 1
+
 END MODULE input_constants
diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F
@@ -71,6 +71,7 @@ MODULE input_cp2k_dft
                                               create_qs_section
    USE input_cp2k_rsgrid,               ONLY: create_rsgrid_section
    USE input_cp2k_scf,                  ONLY: create_scf_section
+   USE input_cp2k_smeagol,              ONLY: create_dft_smeagol_section
    USE input_cp2k_transport,            ONLY: create_transport_section
    USE input_cp2k_xas,                  ONLY: create_xas_section,&
                                               create_xas_tdp_section
@@ -452,6 +453,10 @@ SUBROUTINE create_dft_section(section)
       CALL section_add_subsection(section, subsection)
       CALL section_release(subsection)
 
+      CALL create_dft_smeagol_section(subsection)
+      CALL section_add_subsection(section, subsection)
+      CALL section_release(subsection)
+
    END SUBROUTINE create_dft_section
 
 ! **************************************************************************************************