error unable to open file ......(my directory) Error code=0x80070002 in visual studio 2019 for my c++ program

I wrote the following program using visual studio 2019, which for some reason shows me no compiler errors or anything remotely useful, however every time I run it I get the error: error unable to open file ......(my directory) Error code=0x80070002. The thing is if I exchange the code with something else (that works) everything is fine, so I assume it has to be code related, however I just can t figure out what the problem is. Does anyone understand the issue and could please explain. I add the whole program so that you have exactly what I have.

#include <iostream>
#include <cstdlib>
#include <cmath>
#include <fstream>
#include <ctime>
using namespace std;

float ljpot(float r, float a1, float a2) {
    float lj; float x;
    x = (a2 / r) * (a2 / r) * (a2 / r);
    x = x * x;
    lj = 4 * a1 * x * (x - 1);
    return(lj);
}
float dljpot(float r, float a1, float a2) {
    float dlj; float x;
    x = (a2 / r) * (a2 / r) * (a2 / r);
    x = x * x;
    dlj = -24 * a1 * x * (2 * x - 1);
    return(dlj);
}


float image(float x1, float x2, float l) {
    float im; float x12;
    x12 = x2 - x1;
    x12 = x12 - l * rint(x12 / l);
    im = x12;
    return(im);
}

int main()
{
    int ii;
    bool overlap;
    int eq = 1000;
    const int N = 100;
    float x[N][2], v[N][2], a[N][2], xx[2], rho, l, r, dpot[N], du;
    float dpot[N][N], der[N][N], u, koor[1000], temp, Ekin;
    float pot[N][N] = {};
    float force[N][N][2] = {};
    float acc[N][N] = {};
    ofstream box, part, ekv;
    cout << "density?        ";
    cin >> rho;
    cout << "temperature?        ";
    cin >> temp;
    float sigma = sqrt(temp);
    l = sqrt(float(N) / rho);
    float pi = acos(-1);
    float a1 = 1; // eps LJ potencial
    float a2 = 1; // sigma LJ potencial
    float interval = 0.04;
    srand((int)time(0));


    //random positions and velocities
    x[0][0] = float(rand()) / float(RAND_MAX) * l - l / 2;
    x[0][1] = float(rand()) / float(RAND_MAX) * l - l / 2;
    for (int i = 1; i < N; i++) {
        overlap = true;
        while (overlap) {
            xx[0] = float(rand()) / float(RAND_MAX) * l - l / 2;
            xx[1] = float(rand()) / float(RAND_MAX) * l - l / 2;
            overlap = false;
            for (int j = 0; j < i; j++) {
                r = sqrt(pow(image(xx[0], x[j][0], l), 2) + pow(image(xx[1], x[j][1], l), 2));
                if (r < 1)
                    overlap = true;
            }
        }
        x[i][0] = xx[0];
        x[i][1] = xx[1];
        float ran1 = float(rand()) / float(RAND_MAX);
        float ran2 = float(rand()) / float(RAND_MAX);
        v[i][0] = sigma * sqrt(-2 * log(ran1)) * cos(2 * pi * ran2);
        v[i][1] = sigma * sqrt(-2 * log(ran1)) * sin(2 * pi * ran2);
    }


    box.open("box.txt");
    box << 0 << " " << 0 << endl;
    box << 0 << " " << l << endl;
    box << l << " " << l << endl;
    box << l << " " << 0 << endl;
    box << 0 << " " << 0 << endl;
    box.close();
    part.open("part.txt");
    for (int i = 0; i < N; i++) {
        part << x[i][0] << " " << x[i][1] << " " << v[i][0] << " " << v[i][1] << endl;
    }
    part.close();


    //energy computation
    u = 0;
    for (int i = 1; i < N; i++) {
        for (int j = 0; j < i; j++) {
            if (i != j) {
                r = sqrt(pow(image(x[i][0], x[j][0], l), 2) + pow(image(x[i][1], x[j][1], l), 2));
                pot[i][j] = ljpot(r, a1, a2);
                der[i][j] = dljpot(r, a1, a2);
                pot[j][i] = pot[i][j];
                der[j][i] = der[i][j];
                force[i][j][0] = (-der[i][j] / r / r) * image(x[i][0], x[j][0], l);
                force[i][j][1] = (-der[i][j] / r / r) * image(x[i][1], x[j][1], l);
                force[j][i][0] = -force[i][j][0];
                force[j][i][1] = -force[i][j][1];
                u = u + pot[i][j];

            }
        }
    }
    cout << "inner energy " << u << endl;
    cout << "inner energy per particle " << u / N << endl;


    // equilibration
    ekv.open("ekv.txt");
    int k;
    float dt = 0.0001;
    double sum1;
    double sum2;
    Ekin = 0;
    
    for (int k = 1; k < eq; k++) {
        for (int j = 0; j < N; j++) {
            int p = 0;
            sum1 = 0;
            sum2 = 0;
            while (p < N) {
                sum1 = sum1 + force[j][p][0];
                sum2 = sum2 + force[j][p][0];
            }
            acc[j][0] = sum1;
            acc[j][1] = sum2;
            //speed Verlet alghorhythm
            v[j][0] = v[j][0] + 0.5 * dt * acc[j][0];
            v[j][1] = v[j][1] + 0.5 * dt * acc[j][1];
            //new positions
            x[j][0] = x[j][0] + 0.5 * dt * v[j][0];
            x[j][1] = x[j][1] + 0.5 * dt * v[j][1];
            x[j][0] = image(0, x[j][0], l);
            x[j][1] = image(0, x[j][1], l);
        }
        u = 0;
        for (int i = 1; i < N; i++) {
            for (int j = 0; j < i; j++) {
                if (i != j) {
                    r = sqrt(pow(image(x[i][0], x[j][0], l), 2) + pow(image(x[i][1], x[j][1], l), 2));
                    pot[i][j] = ljpot(r, a1, a2);
                    der[i][j] = dljpot(r, a1, a2);
                    pot[j][i] = pot[i][j];
                    der[j][i] = der[i][j];
                    force[i][j][0] = (-der[i][j] / r / r) * image(x[i][0], x[j][0], l);
                    force[i][j][1] = (-der[i][j] / r / r) * image(x[i][1], x[j][1], l);
                    force[j][i][0] = -force[i][j][0];
                    force[j][i][1] = -force[i][j][1];
                    u = u + pot[i][j];

                }
            }
        }
        for (int j = 0; j < N; j++) {
            int p = 0;
            sum1 = 0;
            sum2 = 0;
            while (p < N) {
                sum1 = sum1 + force[j][p][0];
                sum2 = sum2 + force[j][p][0];
            }
            acc[j][0] = sum1;
            acc[j][1] = sum2;
            v[j][0] = v[j][0] + 0.5 * dt * acc[j][0];
            v[j][1] = v[j][1] + 0.5 * dt * acc[j][1];
            Ekin = Ekin + 0.5 * (v[j][0] * v[j][0] + v[j][1] * v[j][1]);
        }
        Ekin = Ekin / float(N);
        for (int m = 0; m < N; m++) {
            v[m][0] = v[m][0] * sqrt(temp) / sqrt(Ekin);
        }
    }
    ekv.close();
    return 0;
}