We present a novel direct three-dimensional meshless modeling algorithm for geophysical electromagnetic problems. This formulation uses the well known Coulomb gauged, coupled, magnetic vector A and electric scalar Φ decomposition. Assuming a divergence free magnetic vector potential, we recast the commonly used curl-curlequation into a set of four weakly coupled scalar differential equations. These fields are smooththroughout the entire domain, even in the presence of material discontinuities. This makes the decomposition well suited for a meshless discretization approach.To derive a discretization on given point clouds, we apply a novel direct meshless approach. The partial differential equations are written as a set of linear functionals and approximated utilizing the well-known moving least square procedure. Here, the amount of work is shifted to simple polynomials. Thus, no shape functions are needed to create a discretization. We use this kind of approximation to discretize the strong formulation of partial differential equations on multiple processors. To achieve a good load balance and also a suitable domain decomposition, we apply a space filling curve approach. Here, we use the three-dimensional Hilbert curve to derive a decomposition that always leads to connected subdomains. The resulting unsymmetrical sparse linear system is solved with a parallel transpose free QMR implementation. The whole modeling approach is implemented in a fortran coarray structure which runs on machines with several hundreds of cores. We present the mathematical formulation and show example calculations. We also discuss the scalability and performance of this formulation.