Add stress tensor terms with 2nd order Laplace functionals (#2037)

* Add stress tensors with analytical 2nd derivatives for Laplace dependent functionals * Refactor weights for numerical derivatives * Add virial Laplace terms first order * Continue with Laplace numer 2nd * Fix leaks and add regtests
cp2k · Mar 29, 2022 · c73f66f · c73f66f
1 parent 721f5f9
commit c73f66f
Show file tree

Hide file tree

Showing 19 changed files with 813 additions and 405 deletions.
diff --git a/src/input_cp2k_xc.F b/src/input_cp2k_xc.F
@@ -1282,6 +1282,13 @@ SUBROUTINE create_xc_section(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword, __LOCATION__, name="NSTEPS", &
+                          description="Number of steps to consider in each direction for the numerical "// &
+                          "evaluation of XC derivatives. Must be a value from 1 to 4 (Default: 3).", &
+                          usage="NSTEPS 4", default_i_val=3)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       CALL section_add_subsection(section, subsection)
       CALL section_release(subsection)
 

diff --git a/src/xc/xc.F b/src/xc/xc.F
diff --git a/src/xc/xc.fypp b/src/xc/xc.fypp
@@ -62,6 +62,19 @@
 !$OMP       END PARALLEL WORKSHARE
             END IF
          #:endif
+
+         #! Extra code of the Laplacian of the densities
+         #:if arg1.startswith("laplace") and inner_arguments == []
+            IF (my_compute_virial) THEN
+               deriv_att => xc_dset_get_derivative(deriv_set, [deriv_${arg1}$])
+               IF (ASSOCIATED(deriv_att)) THEN
+                  CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
+
+                  virial_pw%cr3d(:, :, :) = -rho1${appendix1}$ (:, :, :)
+                  CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
+               END IF
+            END IF ! my_compute_virial
+         #:endif
       #:endfor
 
    #:enddef

diff --git a/tests/QS/regtest-double-hybrid-grad-numer-meta/TEST_FILES b/tests/QS/regtest-double-hybrid-grad-numer-meta/TEST_FILES
@@ -1,3 +1,3 @@
-H2O_tpss_mp2.inp                                      11      4e-09            -17.159885409905911
+H2O_tpss_mp2.inp                                      11      4e-09            -17.159885409843870
 CH_tpss_mp2.inp                                       11      5e-09             -7.570243568596605
 #EOF
diff --git a/tests/QS/regtest-double-hybrid-grad-numer/H2O_pbe_mp2.inp b/tests/QS/regtest-double-hybrid-grad-numer/H2O_pbe_mp2.inp
@@ -28,6 +28,7 @@
     &END SCF
     &XC
       2ND_DERIV_ANALYTICAL .FALSE.
+      NSTEPS 4
       # Not a very useful functional, just for testing
       &XC_FUNCTIONAL PBE
       &END XC_FUNCTIONAL

diff --git a/tests/QS/regtest-double-hybrid-stress-laplace/CH3_br89_mp2_an.inp b/tests/QS/regtest-double-hybrid-stress-laplace/CH3_br89_mp2_an.inp
@@ -0,0 +1,86 @@
+&GLOBAL                                                                                           
+  PROJECT  CH3_br89_mp2
+  PRINT_LEVEL LOW
+  RUN_TYPE ENERGY_FORCE
+  &TIMINGS
+     THRESHOLD 0.01
+  &END
+&END GLOBAL
+&FORCE_EVAL
+  METHOD Quickstep
+  STRESS_TENSOR  ANALYTICAL
+  &PRINT
+    &FORCES
+    &END FORCES
+    &STRESS_TENSOR
+     COMPONENTS
+    &END STRESS_TENSOR
+  &END
+  &DFT
+    BASIS_SET_FILE_NAME  GTH_BASIS_SETS
+    BASIS_SET_FILE_NAME  HFX_BASIS
+    POTENTIAL_FILE_NAME  POTENTIAL
+    &MGRID
+      CUTOFF  50
+      REL_CUTOFF  20
+    &END MGRID
+    &QS
+      METHOD GPW
+      EPS_DEFAULT 1.0E-10
+    &END QS
+    &SCF
+      SCF_GUESS ATOMIC
+      EPS_SCF 1.0E-5
+      MAX_SCF 100
+    &END SCF
+    &XC
+      # Not a very useful functional, just for testing
+      &XC_FUNCTIONAL
+        &MGGA_X_BR89
+        &END
+      &END XC_FUNCTIONAL
+      &WF_CORRELATION
+        &RI_MP2
+          BLOCK_SIZE  1
+          EPS_CANONICAL 0.1
+          FREE_HFX_BUFFER .TRUE.
+        &END RI_MP2
+        &INTEGRALS
+        &WFC_GPW
+          CUTOFF  20
+          REL_CUTOFF 6
+        &END WFC_GPW
+        &END INTEGRALS
+        MEMORY  500.0
+        NUMBER_PROC  1
+      &END
+    &END XC
+    UKS
+    MULTIPLICITY  2
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC [angstrom] 5.0 5.0 5.0
+    &END CELL
+    &KIND H
+      BASIS_SET  SZV-GTH
+      BASIS_SET RI_AUX  RI_DZVP-GTH
+      POTENTIAL  GTH-HF-q1
+    &END KIND
+    &KIND C
+      BASIS_SET  SZV-GTH
+      BASIS_SET RI_AUX  RI_DZVP-GTH
+      POTENTIAL  GTH-HF-q4
+    &END KIND
+    &COORD 
+        C     0.000     0.000     0.300
+        H     0.000     1.173     0.000
+        H     0.929    -0.536     0.000
+        H    -0.929    -0.536     0.000
+    &END COORD
+    &TOPOLOGY
+      &CENTER_COORDINATES
+      &END
+    &END TOPOLOGY
+  &END SUBSYS
+&END FORCE_EVAL
diff --git a/tests/QS/regtest-double-hybrid-stress-laplace/H2O_br89_mp2_an.inp b/tests/QS/regtest-double-hybrid-stress-laplace/H2O_br89_mp2_an.inp
@@ -0,0 +1,90 @@
+&GLOBAL                                                                                           
+  PROJECT     H2O_br89_mp2
+  PRINT_LEVEL LOW
+  RUN_TYPE    ENERGY_FORCE
+&END GLOBAL
+&FORCE_EVAL
+  METHOD         Quickstep
+  STRESS_TENSOR  ANALYTICAL
+  &PRINT
+    &FORCES
+    &END FORCES
+    &STRESS_TENSOR
+     COMPONENTS
+    &END STRESS_TENSOR
+  &END
+  &DFT
+    BASIS_SET_FILE_NAME  GTH_BASIS_SETS
+    BASIS_SET_FILE_NAME  HFX_BASIS
+    POTENTIAL_FILE_NAME  POTENTIAL
+    &MGRID
+      CUTOFF      200
+      REL_CUTOFF   30
+    &END MGRID
+    &QS
+      METHOD GPW
+      EPS_DEFAULT 1.0E-10
+    &END QS
+    &SCF
+      SCF_GUESS   ATOMIC
+      EPS_SCF     1.0E-6
+      MAX_SCF     100
+      &PRINT
+        &RESTART OFF
+        &END
+      &END
+    &END SCF
+    &XC
+      # Not a very useful functional, just for testing
+      &XC_FUNCTIONAL
+        &MGGA_X_BR89
+        &END
+      &END XC_FUNCTIONAL
+      &WF_CORRELATION
+        &RI_MP2
+          BLOCK_SIZE  1
+          EPS_CANONICAL 0.0001
+          FREE_HFX_BUFFER .TRUE.
+          &CPHF
+            EPS_CONV  1.0E-4
+            MAX_ITER  30
+          &END
+        &END
+        &INTEGRALS
+          &WFC_GPW
+            CUTOFF      25
+            REL_CUTOFF  10
+            EPS_FILTER  1.0E-6
+            EPS_GRID    1.0E-6
+          &END WFC_GPW
+        &END INTEGRALS
+        MEMORY       1.00
+        NUMBER_PROC  1
+      &END
+    &END XC
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC [angstrom]  3.5 3.5 3.5
+    &END CELL
+    &KIND H
+      BASIS_SET         SZV-GTH
+      BASIS_SET RI_AUX  RI_DZVP-GTH
+      POTENTIAL         GTH-PBE-q1
+    &END KIND
+    &KIND O
+      BASIS_SET         SZV-GTH
+      BASIS_SET RI_AUX  RI_DZVP-GTH
+      POTENTIAL         GTH-PBE-q6
+    &END KIND
+    &COORD
+      O       0.000000    0.000000    -0.211000
+      H       0.000000   -0.844000     0.495000
+      H       0.000000    0.744000     0.495000
+    &END
+    &TOPOLOGY
+      &CENTER_COORDINATES
+      &END
+    &END TOPOLOGY
+  &END SUBSYS
+&END FORCE_EVAL
diff --git a/tests/QS/regtest-double-hybrid-stress-laplace/TEST_FILES b/tests/QS/regtest-double-hybrid-stress-laplace/TEST_FILES
@@ -0,0 +1,3 @@
+H2O_br89_mp2_an.inp                                     31       2E-5                 2.55476623536E+01
+CH3_br89_mp2_an.inp                                     31      3e-04                 -2.03435661393E+00
+#EOF
diff --git a/tests/QS/regtest-double-hybrid-stress-numer-laplace/CH3_br89_mp2_numer.inp b/tests/QS/regtest-double-hybrid-stress-numer-laplace/CH3_br89_mp2_numer.inp
@@ -0,0 +1,89 @@
+&GLOBAL                                                                                           
+  PROJECT  CH3_br89_mp2_numer
+  PRINT_LEVEL LOW
+  RUN_TYPE ENERGY_FORCE
+  &TIMINGS
+     THRESHOLD 0.01
+  &END
+&END GLOBAL
+&FORCE_EVAL
+  METHOD Quickstep
+  STRESS_TENSOR  ANALYTICAL
+  &PRINT
+    &FORCES
+    &END FORCES
+    &STRESS_TENSOR
+     COMPONENTS
+    &END STRESS_TENSOR
+  &END
+  &DFT
+    BASIS_SET_FILE_NAME  GTH_BASIS_SETS
+    BASIS_SET_FILE_NAME  HFX_BASIS
+    POTENTIAL_FILE_NAME  POTENTIAL
+    &MGRID
+      CUTOFF  50
+      REL_CUTOFF  20
+    &END MGRID
+    &QS
+      METHOD GPW
+      EPS_DEFAULT 1.0E-10
+    &END QS
+    &SCF
+      SCF_GUESS ATOMIC
+      EPS_SCF 1.0E-5
+      MAX_SCF 100
+    &END SCF
+    &XC
+      2ND_DERIV_ANALYTICAL .FALSE.
+      STEP_SIZE 1.0E-2
+      NSTEPS 1
+      # Not a very useful functional, just for testing
+      &XC_FUNCTIONAL
+        &MGGA_X_BR89
+        &END 
+      &END XC_FUNCTIONAL
+      &WF_CORRELATION
+        &RI_MP2
+          BLOCK_SIZE  1
+          EPS_CANONICAL 0.1
+          FREE_HFX_BUFFER .TRUE.
+        &END RI_MP2
+        &INTEGRALS
+        &WFC_GPW
+          CUTOFF  20
+          REL_CUTOFF 6
+        &END WFC_GPW
+        &END INTEGRALS
+        MEMORY  500.0
+        NUMBER_PROC  1
+      &END
+    &END XC
+    UKS
+    MULTIPLICITY  2
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC [angstrom] 5.0 5.0 5.0
+    &END CELL
+    &KIND H
+      BASIS_SET  SZV-GTH
+      BASIS_SET RI_AUX  RI_DZVP-GTH
+      POTENTIAL  GTH-HF-q1
+    &END KIND
+    &KIND C
+      BASIS_SET  SZV-GTH
+      BASIS_SET RI_AUX  RI_DZVP-GTH
+      POTENTIAL  GTH-HF-q4
+    &END KIND
+    &COORD 
+        C     0.000     0.000     0.300
+        H     0.000     1.173     0.000
+        H     0.929    -0.536     0.000
+        H    -0.929    -0.536     0.000
+    &END COORD
+    &TOPOLOGY
+      &CENTER_COORDINATES
+      &END
+    &END TOPOLOGY
+  &END SUBSYS
+&END FORCE_EVAL