Use Fypp to create code for 2nd derivatives XC automatically

A commit introduced a bug with GAPW which is fixed, too.
More reasonable variable names
Switch POINTER to ALLOCATABLE

src/xc/xc.fypp

@@ -0,0 +1,67 @@
+#!-------------------------------------------------------------------------------------------------!
+#!   CP2K: A general program to perform molecular dynamics simulations                             !
+#!   Copyright 2000-2022 CP2K developers group <https://cp2k.org>                                  !
+#!                                                                                                 !
+#!   SPDX-License-Identifier: GPL-2.0-or-later                                                     !
+#!-------------------------------------------------------------------------------------------------!
+#:mute
+
+#! each list collects the different variable names
+#! element 1: name of quantity to get derivative
+#! element 2: name of spin channel (only used with density gradients)
+#! element 3: name of response quantity variable
+#! element 4: name of variable nameof the contribution to the response potential
+#! elements 5+6: first and last indices of the potential arrays corresponding to the required spin channels
+
+   #:set arguments_openshell = [("rhoa", "a", "rho1a", "v_xc", 1, 1), ("rhob", "b", "rho1b", "v_xc", 2, 2), ("norm_drho", "", "dr1dr", "v_drho", 1, 2), ("norm_drhoa", "a", "dra1dra", "v_drhoa", 1, 1), ("norm_drhob", "b", "drb1drb", "v_drhob", 2, 2)]
+   #:set arguments_triplet_outer = [("rhoa", "a", "rho1a", "v_xc", 1, 1), ("norm_drho", "", "dr1dr", "v_drho", 1, 1), ("norm_drhoa", "a", "dra1dra", "v_drhoa", 1, 1)]
+   #:set arguments_triplet_inner = arguments_triplet_outer + [("rhob", "b", "rho1b", "v_xc", 2, 2), ("norm_drhob", "b", "drb1drb", "v_drhob", 2, 2)]
+   #:set arguments_closedshell = [("rho", "", "rho1", "v_xc", 1, 1), ("norm_drho", "", "dr1dr", "v_drho", 1, 1)]
+
+   #:def add_2nd_derivative_terms(outer_arguments, inner_arguments=[])
+      #:set my_inner_arguments = outer_arguments if inner_arguments==[] else inner_arguments
+      #:for arg1, appendix1, darg1, v_arg1, start1, end1 in outer_arguments
+#! Code for the second order contributions
+         #:for arg2, _, darg2, v_arg2, _, _ in my_inner_arguments
+            deriv_att => xc_dset_get_derivative(deriv_set, "(${arg1}$)(${arg2}$)")
+            IF (ASSOCIATED(deriv_att)) THEN
+               CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
+!$OMP       PARALLEL DO PRIVATE(k,j,i) DEFAULT(NONE)&
+!$OMP       SHARED(bo,${v_arg1}$,deriv_data,${darg2}$#{if v_arg1 != v_arg2}#,${v_arg2}$#{endif}#,fac) COLLAPSE(3)
+               DO k = bo(1, 3), bo(2, 3)
+                  DO j = bo(1, 2), bo(2, 2)
+                     DO i = bo(1, 1), bo(2, 1)
+                        #:for i in range(start1, end1+1)
+                           ${v_arg1}$ (${i}$)%pw%cr3d(i, j, k) = ${v_arg1}$ (${i}$)%pw%cr3d(i, j, k) #{if arg1.startswith("norm")}#-#{else}#+#{endif}# &
+                           #{if ((outer_arguments!=my_inner_arguments) and (arg2[-1]=="b"))}#fac*#{endif}#deriv_data(i, j, k)*${darg2}$ (i, j, k)
+                        #:endfor
+                     END DO
+                  END DO
+               END DO
+            END IF
+         #:endfor
+
+#! Extra code of the gradients
+#:if arg1.startswith("norm")
+         deriv_att => xc_dset_get_derivative(deriv_set, "(${arg1}$)")
+         IF (ASSOCIATED(deriv_att)) THEN
+            CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
+            CALL xc_derivative_get(deriv_att, deriv_data=e_drho${appendix1}$)
+
+            IF (my_compute_virial) THEN
+               CALL virial_drho_drho1(virial_pw, drho${appendix1}$, drho1${appendix1}$, deriv_data, virial_xc)
+            END IF ! my_compute_virial
+
+!$OMP       PARALLEL WORKSHARE DEFAULT(NONE) SHARED(${darg1}$,gradient_cut,${arg1}$,${v_arg1}$,deriv_data)
+#:for i in range(start1, end1+1)
+            ${v_arg1}$(${i}$)%pw%cr3d(:, :, :) = ${v_arg1}$(${i}$)%pw%cr3d(:, :, :) + &
+                                          deriv_data(:, :, :)*${darg1}$(:, :, :)/MAX(gradient_cut, ${arg1}$(:, :, :))**2
+#:endfor
+!$OMP       END PARALLEL WORKSHARE
+         END IF
+#:endif
+#:endfor
+
+   #:enddef
+
+#:endmute

tests/QS/regtest-debug-1/TEST_FILES

@@ -8,5 +8,5 @@ h2o_polar.inp                                         87      1e-05
 h2o_pdip.inp                                          86      1e-05             0.961018561809E+00
 h2o_periodic.inp                                      87      1e-05             0.139741440657E+02
 h2o_gga.inp                                           87      1e-05             0.165143521103E+02
-h2o_gapw.inp                                          87      1e-05             0.170647717900E+02
+h2o_gapw.inp                                          87      1e-05             0.170220498303E+02
 #EOF