#1293 and #1327 (#1331)

* Initial guess based on standard core definitions.

* Code prettify

* Possible fixes for:
calculate_vxc_atom references unassociated pointer drho_h (#1289)
calculate_xc_2nd_deriv_atom references disassociated pointer (#1292)
qs_ks_build_kohn_sham_matrix references undefined variable virial%pv_xc (#1293)

* Fix code mess up

* Remove debug code

* Conditional initialization of drho and tau variables. (#1327)
Unconditional initialization of virial_xc (intent(out) variable) (#1293)

src/qs_vxc_atom.F

@@ -493,6 +493,7 @@ SUBROUTINE calculate_xc_2nd_deriv_atom(local_rho_set, qs_env, xc_section, para_e
                                                             paw_atom, tau_f
       REAL(KIND=dp)                                      :: density_cut, gradient_cut, rtot, tau_cut
       REAL(KIND=dp), DIMENSION(1, 1)                     :: rtau
+      REAL(KIND=dp), DIMENSION(1, 1, 1, 1)               :: rrho
       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rho1_h, rho1_s, rho_h, rho_s, tau1_h, &
                                                             tau1_s, tau_h, tau_s, weight
       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: vxc_h, vxc_s
@@ -639,6 +640,9 @@ SUBROUTINE calculate_xc_2nd_deriv_atom(local_rho_set, qs_env, xc_section, para_e
             ALLOCATE (vxg_h(1:3, 1:na, 1:nr, 1:nspins), vxg_s(1:3, 1:na, 1:nr, 1:nspins))
             vxg_h = 0.0_dp
             vxg_s = 0.0_dp
+         ELSE
+            NULLIFY (vxg_h, vxg_s)
+            rrho = 0.0_dp
          END IF
 
          ! parallelization
@@ -690,11 +694,16 @@ SUBROUTINE calculate_xc_2nd_deriv_atom(local_rho_set, qs_env, xc_section, para_e
                   CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, tau1_h, na, ir)
                   CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
                   CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, tau1_s, na, ir)
-               ELSE
+               ELSE IF (gradient_functional) THEN
                   CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, rtau, na, ir)
                   CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, rtau, na, ir)
                   CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, rtau, na, ir)
                   CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, rtau, na, ir)
+               ELSE
+                  CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rrho, rtau, na, ir)
+                  CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, rrho, rtau, na, ir)
+                  CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rrho, rtau, na, ir)
+                  CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, rrho, rtau, na, ir)
                END IF
 
             END DO

src/xc/xc.F

@@ -129,6 +129,8 @@ SUBROUTINE xc_vxc_pw_create1(vxc_rho, vxc_tau, rho_r, rho_g, tau, exc, xc_sectio
       CPASSERT(.NOT. ASSOCIATED(vxc_rho))
       CPASSERT(.NOT. ASSOCIATED(vxc_tau))
 
+      virial_xc = 0.0_dp
+
       CALL section_vals_val_get(xc_section, "FUNCTIONAL_ROUTINE", &
                                 i_val=f_routine)
       SELECT CASE (f_routine)
@@ -1116,6 +1118,7 @@ SUBROUTINE xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau, xc_section
       END IF
 
       use_virial = compute_virial
+      virial_xc = 0.0_dp
 
       bo = rho_r(1)%pw%pw_grid%bounds_local
       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)