Add Voronoi and Maxwell Eq. interfaces

Co-authored-by: Jürg Hutter <hutter@chem.uzh.ch>

INSTALL.md

@@ -303,6 +303,21 @@ SIRIUS is a domain specific library for electronic structure calculations.
   library should be placed in front of the scalapack library.
 - see <https://github.com/eth-cscs/COSMA> for more information.
 
+### 2t. LibVori (Voronoi Integration for Electrostatic Properties from Electron Density)
+
+- LibVori is a library which enables to compute electrostatic properties
+  (charge, dipole vector, quadrupole tensor, ...) via integration of the total
+  electron density in the Voronoi cell of each atom.
+- Add `-D__LIBVORI` to DFLAGS to enable support for LibVori.
+- see <https://brehm-research.de/voronoi> for more information.
+
+### 2u. LibMaxwell (External Maxwell Solver)
+
+- LibMaxwell is a library to solve the time-dependent Maxwell equations
+  and use the resulting electric field in MD runs or real-time propagation.
+- Add `-D__LIBMAXWELL` to DFLAGS to enable support for LibMaxwell.
+- see <https://brehm-research.de> for more information.
+
 ## 3. Compile
 
 ### 3a. ARCH files

src/common/bibliography.F

@@ -90,7 +90,8 @@ MODULE bibliography
                             Ren2011, Ren2013, Cohen2000, Rogers2002, Filippetti2000, Paziani2006, &
                             Toulouse2004, Limpanuparb2011, Martin2003, Yin2017, Goerigk2017, &
                             Wilhelm2016a, Wilhelm2016b, Wilhelm2017, Wilhelm2018, Lass2018, cp2kqs2020, &
-                            Behler2007, Behler2011, Schran2020a, Schran2020b
+                            Behler2007, Behler2011, Schran2020a, Schran2020b, &
+                            Rycroft2009, Thomas2015, Brehm2018, Brehm2020
 
 CONTAINS
 
@@ -4592,6 +4593,88 @@ SUBROUTINE add_all_references()
                          "ER"), &
                          DOI="10.1063/5.0007045")
 
+      CALL add_reference(key=Rycroft2009, ISI_record=s2a( &
+                         "TY JOUR", &
+                         "PT J", &
+                         "AU Rycroft, Chris H", &
+                         "TI VORO++: A three-dimensional Voronoi cell library in C++", &
+                         "SO Chaos: An Interdisciplinary Journal of Nonlinear Science", &
+                         "JF Chaos: An Interdisciplinary Journal of Nonlinear Science", &
+                         "JO Chaos", &
+                         "SN 1054-1500", &
+                         "IS 4", &
+                         "PY 2009", &
+                         "VL 19", &
+                         "SP 041111", &
+                         "DI 10.1063/1.3215722", &
+                         "PB American Institute of Physics", &
+                         "ER"), &
+                         DOI="10.1063/1.3215722")
+
+      CALL add_reference(key=Thomas2015, ISI_record=s2a( &
+                         "TY JOUR", &
+                         "PT J", &
+                         "AU Thomas, Martin", &
+                         "   Brehm, Martin", &
+                         "   Kirchner, Barbara", &
+                         "TI Voronoi dipole moments for the simulation of bulk phase vibrational", &
+                         "   spectra", &
+                         "SO Physical Chemistry Chemical Physics", &
+                         "JF Physical Chemistry Chemical Physics", &
+                         "JO Phys. Chem. Chem. Phys.", &
+                         "SN 1463-9076", &
+                         "IS 5", &
+                         "PY 2015", &
+                         "VL 17", &
+                         "SP 3207", &
+                         "EP 3213", &
+                         "DI 10.1039/C4CP05272B", &
+                         "PB The Royal Society of Chemistry", &
+                         "ER"), &
+                         DOI="10.1039/C4CP05272B")
+
+      CALL add_reference(key=Brehm2018, ISI_record=s2a( &
+                         "TY JOUR", &
+                         "PT J", &
+                         "AU Brehm, Martin", &
+                         "   Thomas, Martin", &
+                         "TI An Efficient Lossless Compression Algorithm for Trajectories of Atom", &
+                         "   Positions and Volumetric Data", &
+                         "SO Journal of Chemical Information and Modeling", &
+                         "JF Journal of Chemical Information and Modeling", &
+                         "JO J. Chem. Inf. Model.", &
+                         "SN 1549-9596", &
+                         "IS 10", &
+                         "PY 2018", &
+                         "VL 58", &
+                         "SP 2092", &
+                         "EP 2107", &
+                         "DI 10.1021/acs.jcim.8b00501", &
+                         "PB American Chemical Society", &
+                         "ER"), &
+                         DOI="10.1021/acs.jcim.8b00501")
+
+      CALL add_reference(key=Brehm2020, ISI_record=s2a( &
+                         "TY JOUR", &
+                         "PT J", &
+                         "AU Brehm, Martin", &
+                         "   Thomas, Martin", &
+                         "   Gehrke, Sascha", &
+                         "   Kirchner, Barbara", &
+                         "TI TRAVIS - A free analyzer for trajectories from molecular simulation", &
+                         "SO The Journal of Chemical Physics", &
+                         "JF The Journal of Chemical Physics", &
+                         "JO J. Chem. Phys.", &
+                         "SN 0021-9606", &
+                         "IS 16", &
+                         "PY 2020", &
+                         "VL 152", &
+                         "SP 164105", &
+                         "DI 10.1063/5.0005078", &
+                         "PB American Institute of Physics", &
+                         "ER"), &
+                         DOI="10.1063/5.0005078")
+
    END SUBROUTINE add_all_references
 
 END MODULE bibliography

src/cp2k_info.F

@@ -273,6 +273,13 @@ FUNCTION cp2k_flags() RESULT(flags)
 #if defined __CHECK_DIAG
       flags = TRIM(flags)//" check_diag"
 #endif
+#if defined __LIBVORI
+      flags = TRIM(flags)//" libvori"
+      flags = TRIM(flags)//" libbqb"
+#endif
+#if defined __LIBMAXWELL
+      flags = TRIM(flags)//" libmaxwell"
+#endif
 
 #if defined __GRID_CUDA
       flags = TRIM(flags)//" grid_cuda"

src/cp_control_types.F

@@ -268,6 +268,25 @@ MODULE cp_control_types
          ALLOCATABLE                       :: blocks
    END TYPE admm_control_type
 
+! **************************************************************************************************
+! \brief Parameters for external potential
+! **************************************************************************************************
+   TYPE expot_control_type
+      LOGICAL                              :: read_from_cube
+      LOGICAL                              :: maxwell_solver
+      LOGICAL                              :: static
+      REAL(KIND=dp)                        :: scaling_factor
+   END TYPE expot_control_type
+
+! **************************************************************************************************
+! \brief Parameters useful for Maxwell equation evaluation of external potential
+! **************************************************************************************************
+   TYPE maxwell_control_type
+      LOGICAL                              :: log_test
+      INTEGER                              :: int_test
+      REAL(KIND=dp)                        :: real_test
+   END TYPE maxwell_control_type
+
 ! **************************************************************************************************
 ! \brief Control parameters for a QUICKSTEP and KIM-GORDON calculation ***
 !        eps_pgf_orb: Cutoff value for the interaction of the primitive
@@ -489,6 +508,8 @@ MODULE cp_control_types
       TYPE(tddfpt_control_type), POINTER   :: tddfpt_control
       TYPE(tddfpt2_control_type), POINTER  :: tddfpt2_control
       TYPE(xas_control_type), POINTER      :: xas_control
+      TYPE(expot_control_type), POINTER    :: expot_control
+      TYPE(maxwell_control_type), POINTER  :: maxwell_control
       TYPE(efield_p_type), POINTER, &
          DIMENSION(:)                      :: efield_fields
       INTEGER                              :: nspins, &
@@ -562,6 +583,8 @@ MODULE cp_control_types
              semi_empirical_control_type, &
              s2_restraint_type, &
              admm_control_type, &
+             maxwell_control_type, &
+             expot_control_type, &
              rtp_control_type, &
              sccs_control_type, &
              stda_control_type
@@ -573,6 +596,8 @@ MODULE cp_control_types
              dft_control_create, &
              tddfpt_control_create, &
              admm_control_create, &
+             maxwell_control_create, &
+             expot_control_create, &
              ddapc_control_create, &
              sccs_control_create
 
@@ -771,6 +796,8 @@ SUBROUTINE dft_control_create(dft_control)
       NULLIFY (dft_control%efield_fields)
       NULLIFY (dft_control%period_efield)
       NULLIFY (dft_control%admm_control)
+      NULLIFY (dft_control%expot_control)
+      NULLIFY (dft_control%maxwell_control)
       NULLIFY (dft_control%rtp_control)
       NULLIFY (dft_control%sccs_control)
       dft_control%do_sccs = .FALSE.
@@ -814,6 +841,8 @@ SUBROUTINE dft_control_release(dft_control)
             CALL tddfpt2_control_release(dft_control%tddfpt2_control)
             CALL xas_control_release(dft_control%xas_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 efield_fields_release(dft_control%efield_fields)
             CALL sccs_control_release(dft_control%sccs_control)
             IF (ASSOCIATED(dft_control%period_efield)) THEN
@@ -1187,4 +1216,58 @@ SUBROUTINE sccs_control_release(sccs_control)
 
    END SUBROUTINE sccs_control_release
 
+! **************************************************************************************************
+!> \brief ...
+!> \param expot_control ...
+! **************************************************************************************************
+   SUBROUTINE expot_control_create(expot_control)
+      TYPE(expot_control_type), POINTER                  :: expot_control
+
+      CPASSERT(.NOT. ASSOCIATED(expot_control))
+      ALLOCATE (expot_control)
+      expot_control%read_from_cube = .FALSE.
+      expot_control%maxwell_solver = .FALSE.
+      expot_control%static = .TRUE.
+      expot_control%scaling_factor = 1.0_dp
+
+   END SUBROUTINE expot_control_create
+
+! **************************************************************************************************
+!> \brief ...
+!> \param expot_control ...
+! **************************************************************************************************
+   SUBROUTINE expot_control_release(expot_control)
+      TYPE(expot_control_type), POINTER                  :: expot_control
+
+      IF (ASSOCIATED(expot_control)) THEN
+         DEALLOCATE (expot_control)
+      END IF
+
+   END SUBROUTINE expot_control_release
+
+! **************************************************************************************************
+!> \brief ...
+!> \param maxwell_control ...
+! **************************************************************************************************
+   SUBROUTINE maxwell_control_create(maxwell_control)
+      TYPE(maxwell_control_type), POINTER                :: maxwell_control
+
+      CPASSERT(.NOT. ASSOCIATED(maxwell_control))
+      ALLOCATE (maxwell_control)
+
+   END SUBROUTINE maxwell_control_create
+
+! **************************************************************************************************
+!> \brief ...
+!> \param maxwell_control ...
+! **************************************************************************************************
+   SUBROUTINE maxwell_control_release(maxwell_control)
+      TYPE(maxwell_control_type), POINTER                :: maxwell_control
+
+      IF (ASSOCIATED(maxwell_control)) THEN
+         DEALLOCATE (maxwell_control)
+      END IF
+
+   END SUBROUTINE maxwell_control_release
+
 END MODULE cp_control_types

src/cp_control_utils.F

@@ -16,8 +16,9 @@ MODULE cp_control_utils
         VanVoorhis2015, VandeVondele2005a, VandeVondele2005b, Yin2017, Zhechkov2005, cite_reference
    USE cp_control_types,                ONLY: &
         admm_control_create, admm_control_type, ddapc_control_create, ddapc_restraint_type, &
-        dft_control_create, dft_control_type, efield_type, qs_control_type, sccs_control_create, &
-        tddfpt2_control_type, tddfpt_control_create, tddfpt_control_type, xtb_control_type
+        dft_control_create, dft_control_type, efield_type, expot_control_create, &
+        maxwell_control_create, qs_control_type, sccs_control_create, tddfpt2_control_type, &
+        tddfpt_control_create, tddfpt_control_type, xtb_control_type
    USE cp_files,                        ONLY: close_file,&
                                               open_file
    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
@@ -102,7 +103,8 @@ SUBROUTINE read_dft_control(dft_control, dft_section)
       REAL(KIND=dp)                                      :: density_cut, gradient_cut, tau_cut
       REAL(KIND=dp), DIMENSION(:), POINTER               :: pol
       TYPE(cp_logger_type), POINTER                      :: logger
-      TYPE(section_vals_type), POINTER                   :: sccs_section, scf_section, tmp_section, &
+      TYPE(section_vals_type), POINTER                   :: maxwell_section, sccs_section, &
+                                                            scf_section, tmp_section, &
                                                             xc_fun_section, xc_section
       TYPE(xc_rho_cflags_type)                           :: needs
 
@@ -436,6 +438,31 @@ SUBROUTINE read_dft_control(dft_control, dft_section)
       ! Read the external potential input section
       tmp_section => section_vals_get_subs_vals(dft_section, "EXTERNAL_POTENTIAL")
       CALL section_vals_get(tmp_section, explicit=dft_control%apply_external_potential)
+      IF (dft_control%apply_external_potential) THEN
+         CALL expot_control_create(dft_control%expot_control)
+         CALL section_vals_val_get(tmp_section, "READ_FROM_CUBE", &
+                                   l_val=dft_control%expot_control%read_from_cube)
+         CALL section_vals_val_get(tmp_section, "STATIC", &
+                                   l_val=dft_control%expot_control%static)
+         CALL section_vals_val_get(tmp_section, "SCALING_FACTOR", &
+                                   r_val=dft_control%expot_control%scaling_factor)
+         ! External potential using Maxwell equation
+         maxwell_section => section_vals_get_subs_vals(tmp_section, "MAXWELL")
+         CALL section_vals_get(maxwell_section, explicit=is_present)
+         IF (is_present) THEN
+            dft_control%expot_control%maxwell_solver = .TRUE.
+            CALL maxwell_control_create(dft_control%maxwell_control)
+            ! read the input values from Maxwell section
+            CALL section_vals_val_get(maxwell_section, "TEST_REAL", &
+                                      r_val=dft_control%maxwell_control%real_test)
+            CALL section_vals_val_get(maxwell_section, "TEST_INTEGER", &
+                                      i_val=dft_control%maxwell_control%int_test)
+            CALL section_vals_val_get(maxwell_section, "TEST_LOGICAL", &
+                                      l_val=dft_control%maxwell_control%log_test)
+         ELSE
+            dft_control%expot_control%maxwell_solver = .FALSE.
+         END IF
+      END IF
 
       ! Read the SCCS input section if present
       sccs_section => section_vals_get_subs_vals(dft_section, "SCCS")

src/input_constants.F

@@ -608,6 +608,19 @@ MODULE input_constants
                                                radius_vdw = 13, &
                                                radius_default = 14
 
+   ! Voronoi integration
+   INTEGER, PARAMETER, PUBLIC               :: voro_radii_unity = 100, &
+                                               voro_radii_vdw = 101, &
+                                               voro_radii_cov = 102, &
+                                               voro_radii_user = 103
+
+   ! BQB Format
+   INTEGER, PARAMETER, PUBLIC               :: bqb_opt_off = 100, &
+                                               bqb_opt_quick = 101, &
+                                               bqb_opt_normal = 102, &
+                                               bqb_opt_patient = 103, &
+                                               bqb_opt_exhaustive = 104
+
    ! MAO
    INTEGER, PARAMETER, PUBLIC               :: mao_basis_orb = 2000, &
                                                mao_basis_prim = 2001, &