/opt/arm/gcc-14.2.0_AmazonLinux-2023/lib/gcc/aarch64-linux-gnu/14.2.0/../../../../include/c++/14.2.0/bits/stl_vector.h: 993 - 1150
--------------------------------------------------------------------------------

993:       { return size_type(this->_M_impl._M_finish - this->_M_impl._M_start); }
[...]
1131: 	return *(this->_M_impl._M_start + __n);
[...]
1150: 	return *(this->_M_impl._M_start + __n);
/home/fmusial/MD_Benchmarks/simulation.cpp: 140 - 229
--------------------------------------------------------------------------------

140:     #pragma omp parallel for collapse(3) schedule(runtime)
141:     #endif
142:     for (int cx = 0; cx < cellList.num_cells; ++cx) {
143:         for (int cy = 0; cy < cellList.num_cells; ++cy) {
144:             for (int cz = 0; cz < cellList.num_cells; ++cz) {
145:                 int cellIndex = cx + cy * cellList.num_cells + cz * cellList.num_cells * cellList.num_cells;
[...]
158:                     for (int n = 0; n < 13; ++n) {
159: 
160:                                 int dx = neighbor_patterns[n][0];
161:                                 int dy = neighbor_patterns[n][1];
162:                                 int dz = neighbor_patterns[n][2];
163: 
164:                                 int nx = (cx + dx + cellList.num_cells) % cellList.num_cells;
165:                                 int ny = (cy + dy + cellList.num_cells) % cellList.num_cells;
166:                                 int nz = (cz + dz + cellList.num_cells) % cellList.num_cells;
167:                                 
168:                                 int neighborCellIndex = nx + ny * cellList.num_cells + nz * cellList.num_cells * cellList.num_cells;
169:                                 
170:                                 // Bounds check for neighbor cell
171:                                 if (neighborCellIndex < 0 || neighborCellIndex >= static_cast<int>(cellList.cells.size())) {
172:                                     continue;
173:                                 }
174:                                 
175:                                 // Iterate through particles in this cell
176:                                 for (size_t pi = 0; pi < cellList.cells[cellIndex].size(); ++pi) {
177:                                     int i = cellList.cells[cellIndex][pi];
178:                                     // Check particles in neighboring cell
179:                                     for (size_t pj = 0; pj < cellList.cells[neighborCellIndex].size(); ++pj) {
180:                                         int j = cellList.cells[neighborCellIndex][pj];
181:                                         
182:                                         double dx = particles.x[i] - particles.x[j];
183:                                         double dy = particles.y[i] - particles.y[j];
184:                                         double dz = particles.z[i] - particles.z[j];
185:                                         
186:                                         // Apply minimum image convention for periodic boundary conditions
187:                                         dx -= BOX_SIZE * std::round(dx / BOX_SIZE); //NOTE : precompute 1/BOX_SIZE and multiply by it ?
188:                                         dy -= BOX_SIZE * std::round(dy / BOX_SIZE);
189:                                         dz -= BOX_SIZE * std::round(dz / BOX_SIZE);
190:                                         
191:                                         double r2 = dx*dx + dy*dy + dz*dz;
192:                                                                         
193:                                         if (r2 < CUTOFF_SQUARED) {
194:                                             // Compute LJ force
195:                                             double r1i = SIGMA / sqrt(r2); //NOTE: multiply by inverse distance
196:                                             double r2i = r1i * r1i;
197:                                             double r6i = r2i * r2i * r2i;
198:                                             double force = 48.0 * EPSILON * r1i * r2i * r6i * (r6i - 0.5);
199:                                             
200:                                             double fx = force * dx;
201:                                             double fy = force * dy;
202:                                             double fz = force * dz;
[...]
208:                                             particles.fx[i] += fx;
209:                                             #ifdef USE_OPENMP
210:                                             #pragma omp atomic
211:                                             #endif
212:                                             particles.fy[i] += fy;
213:                                             #ifdef USE_OPENMP
214:                                             #pragma omp atomic
215:                                             #endif
216:                                             particles.fz[i] += fz;
217:                                             
218:                                             #ifdef USE_OPENMP
219:                                             #pragma omp atomic
220:                                             #endif
221:                                             particles.fx[j] -= fx;
222:                                             #ifdef USE_OPENMP
223:                                             #pragma omp atomic
224:                                             #endif
225:                                             particles.fy[j] -= fy;
226:                                             #ifdef USE_OPENMP
227:                                             #pragma omp atomic
228:                                             #endif
229:                                             particles.fz[j] -= fz;