diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx
index 15200ca..7f41b30 100644
--- a/qpms/qpms_c.pyx
+++ b/qpms/qpms_c.pyx
@@ -1096,10 +1096,10 @@ cdef class _ScatteringSystemAtOmegaK:
         cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
         cdef np.ndarray[complex, ndim=3] results = np.empty((evalpos_a.shape[0], fecv_size, 3), dtype=complex)
         cdef ccart3_t *res
-        res = <ccart3_t *> malloc(fecv_size*sizeof(ccart3_t))
         cdef cart3_t pos
         cdef Py_ssize_t i, j
         with nogil, wraparound(False), parallel():
+            res = <ccart3_t *> malloc(fecv_size*sizeof(ccart3_t))
             for i in prange(evalpos_a.shape[0]):
                 pos.x = evalpos_a[i,0]
                 pos.y = evalpos_a[i,1]
@@ -1109,7 +1109,7 @@ cdef class _ScatteringSystemAtOmegaK:
                     results[i,j,0] = res[j].x
                     results[i,j,1] = res[j].y
                     results[i,j,2] = res[j].z
-        free(res)
+            free(res)
         return results.reshape(evalpos.shape[:-1] + (self.fecv_size, 3))
     
     property fecv_size: 
@@ -1329,10 +1329,10 @@ cdef class _ScatteringSystemAtOmega:
         cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
         cdef np.ndarray[complex, ndim=3] results = np.empty((evalpos_a.shape[0], basissize, 3), dtype=complex)
         cdef ccart3_t *res
-        res = <ccart3_t *> malloc(basissize*sizeof(ccart3_t))
         cdef cart3_t pos
         cdef Py_ssize_t i, j
         with nogil, wraparound(False), parallel():
+            res = <ccart3_t *> malloc(basissize*sizeof(ccart3_t)) # thread-local
             for i in prange(evalpos_a.shape[0]):
                 pos.x = evalpos_a[i,0]
                 pos.y = evalpos_a[i,1]
@@ -1345,7 +1345,7 @@ cdef class _ScatteringSystemAtOmega:
                     results[i,j,0] = res[j].x
                     results[i,j,1] = res[j].y
                     results[i,j,2] = res[j].z
-        free(res)
+            free(res)
         return results.reshape(evalpos.shape[:-1] + (basissize, 3))
 
     def bspec_pi(self, pi):