diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c
index 8d56866..81f1ab5 100644
--- a/qpms/scatsystem.c
+++ b/qpms/scatsystem.c
@@ -1295,7 +1295,6 @@ static inline complex double *qpms_scatsysw_scatsyswk_build_modeproblem_matrix_f
     QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double));
   complex double *tmp; // translation matrix, S or W
   QPMS_CRASHING_MALLOC(tmp, SQ(ss->max_bspecn) * sizeof(complex double));
-  memset(target, 0, SQ(full_len) * sizeof(complex double));
   const complex double zero = 0, minusone = -1;
   { // Non-diagonal part; M[piR, piC] = -T[piR] S(piR<-piC)
     size_t fullvec_offsetR = 0;
@@ -1321,7 +1320,14 @@ static inline complex double *qpms_scatsysw_scatsyswk_build_modeproblem_matrix_f
                 QPMS_EWALD_FULL, eta));
         }
 
-        cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
+        if (k == NULL && piC == piR) { 
+          // non-periodic case diagonal block: no "self-interaction", just
+          // fill with zeroes (the ones on the diagonal are added in the very end)
+          QPMS_ENSURE_SUCCESS(LAPACKE_zlaset(CblasRowMajor, 'x',
+                bspecR->n /*m*/, bspecC->n/*n*/, 0 /*alfa: offdiag*/, 0 /*beta: diag*/, 
+                target + fullvec_offsetR*full_len + fullvec_offsetC /*a*/,
+                full_len /*lda*/));
+        } else cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
             bspecR->n /*m*/, bspecC->n /*n*/, bspecR->n /*k*/, 
             &minusone/*alpha*/, tmmR/*a*/, bspecR->n/*lda*/,
             tmp/*b*/, bspecC->n/*ldb*/, &zero/*beta*/,