Random issue with `symmetrizer.F90`

The assertion all(nint(srcpoint) >= 0) fails, but when using print, these assertions no longer fail. I did not test on different compilers, optimization levels etc.

gs.inp

## Calculation parameters
CalculationMode = gs
FromScratch = False
ExperimentalFeatures = yes
Debug = no
DebugTrapSignals = yes

## Paralleliztaion parameters
ParStates = auto
ParDomains = no
ParKPoints = auto

## Physical parameters
Dimensions = 3
PeriodicDimensions = 2

# Real-space parameters
BoxShape = parallelepiped
a_par = 4.7621094
Spacing = 0.38

Lz = 25  # fixed: Vacuum

%LatticeParameters
  a_par | a_par | 2 * Lz
%

%LatticeVectors
   1   | 0         | 0
  -1/2 | sqrt(3)/2 | 0
   0   | 0         | 1
%

# Reciprocal-space parameters
nk = 36
%KPointsGrid
 nk | nk | 1
%

KPointsUseSymmetries = no
KPointsUseTimeReversal = no
SymmetrizeDensity = no

# Species parameters
PseudopotentialSet = hgh_lda
% ReducedCoordinates
 "N" | 0 | 0 | 0 | no # yes
 "B" | 1/3 | 2/3 | 0 | no #yes
%

# Occupation and extra states
ExtraStates = 0

## CalcMode parameters

## Extension parameters

## Output parameters
OutputFormat = netcdf

td.inp

## Calculation parameters
CalculationMode = td
FromScratch = True
ExperimentalFeatures = yes
Debug = no
DebugTrapSignals = yes

## Paralleliztaion parameters
ParStates = auto
ParDomains = no
ParKPoints = auto

## Physical parameters
Dimensions = 3
PeriodicDimensions = 2

# Real-space parameters
BoxShape = parallelepiped
a_par = 4.7621094
Spacing = 0.38

Lz = 25  # fixed: Vacuum

%LatticeParameters
  a_par | a_par | 2 * Lz
%

%LatticeVectors
   1   | 0         | 0
  -1/2 | sqrt(3)/2 | 0
   0   | 0         | 1
%

# Reciprocal-space parameters
nk = 36
%KPointsGrid
 nk | nk | 1
%

KPointsUseSymmetries = no
KPointsUseTimeReversal = no
SymmetrizeDensity = no

# Species parameters
PseudopotentialSet = hgh_lda
% ReducedCoordinates
 "N" | 0 | 0 | 0 | no # yes
 "B" | 1/3 | 2/3 | 0 | no #yes
%

# Occupation and extra states
ExtraStates = 0

## CalcMode parameters
# Propagator parameters
TDPropagator = aetrs
TDSystemPropagator = prop_aetrs  # Dummy value because we are using the legacy code
TDEXPOrder = 16
TDExponentialMethod = lanczos

# Physical propagation model
MoveIons = no
TDFreezeHXC = no

# Absorber parameters (for non-periodic part)
wid = 12  # width of absorber counted from Lz down.
Lmin = Lz - wid;
AbsorbingBoundaries = cap
ABCapHeight = -1.0
%ABShape
  Lmin | Lz | "abs(z)"
%

# laser details:
I0 = 1000000000000.0
lam = 1600

# resulting qunatities in a.u.:
light_speed = 137.03640  # in a.u. - 1 over the f.s.c.
amp = 5.33685e-9 * sqrt(I0)  # derived in a.u from I0.
om = 2 * pi * (1 / 7.2973525664e-3) / (10 * lam / 0.52917725)  # the resulting frequency in atomic units.
t_cyc = 2 * pi / om
n_cyc = 8
T_pulse = n_cyc * t_cyc

TDPropagationTime = T_pulse
TDTimeStep = 0.29

%TDExternalFields
  vector_potential | -i | 0 | 0 | om | "envelope_super_sin"
%

sigma = 0.75; # sigma parameter for super-sin.
%TDFunctions
  "envelope_super_sin"   | tdf_from_expr | "-(amp/om) * light_speed * (sin(pi*t/T_pulse)^(abs(pi*t/T_pulse-pi/2)/sigma))*(1-step(t-T_pulse))*(step(t))"
  "envelope_sine_square" | tdf_from_expr | "-(amp/om) * light_speed * sin(pi*t/T_pulse)^2"
%

## Output parameters
OutputInterval = 50
OutputFormat = netcdf

%TDOutput
  laser
  total_current

%
TDMultipoleLmax = 0