My generic normal generation function looks like this:

static void calculate_vertex_normals(vector<vec3>& vx, vector<int>& faces_tri_idx, vector<vec3>& out_normals)
{
    vector<vector<int>> adj_faces(vx.size());

    for (int i = 0; i < vx.size(); i++)
    {
	for (int j = 0; j < faces_tri_idx.size() / 3; j++)
	{
	    for (int k = 0; k < 3; k++)
	    {
	        if (faces_tri_idx[j * 3 + k] == i)
	        adj_faces[i].push_back(j);
	    }
        }
    }

    vector<vec3> face_n(faces_tri_idx.size() / 3);

    for (int i = 0; i < faces_tri_idx.size() / 3; i++)
    {
	const vec3& v1 = vx[faces_tri_idx[i * 3 + 0]];
	const vec3& v2 = vx[faces_tri_idx[i * 3 + 1]];
	const vec3& v3 = vx[faces_tri_idx[i * 3 + 2]];

	vec3 vv1 = v2 - v1;
	vec3 vv2 = v3 - v1;

	vec3 norm = normalize(cross(vv1, vv2));

	face_n[i] = norm;
    }

    for (int i = 0; i < adj_faces.size(); i++)
    {
        auto& vf = adj_faces[i];

	vec3 sum(0.0f);

	for (int j = 0; j < vf.size(); j++)
	{
	    sum += face_n[vf[j]];
	}

	if (abs(length2(sum))>1e-9)
	{
            out_normals.push_back(normalize(sum));
	}
	else
	{
	    out_normals.push_back(vec3(0.0f, 1.0f, 0.0f));
        }
    }
}

Box generation function:

static void box(vector<vec3>& out_vertices, vector<vec3>& out_normals, vec3& sz)
{
    static vec3 cube_v[] =
    {
	//top face view
	//v3----v2
	//|     |
        //v4----v1

	//y=1
	{ 1.0f, 1.0f,  1.0f },
	{ 1.0f, 1.0f, -1.0f },
	{ -1.0f, 1.0f, -1.0f },
	{ -1.0f, 1.0f,  1.0f },

	//y=-1
	{ 1.0f, -1.0f,  1.0f },
	{ 1.0f, -1.0f, -1.0f },
	{ -1.0f, -1.0f, -1.0f },
	{ -1.0f, -1.0f,  1.0f },
    };

    static int  cube_i[] =
    {
	//y=-1
	0,1,3,3,1,2,
	//y=-1
	7,6,4,4,6,5,
	//x=1
	4,5,0,0,5,1,
	//x=-1
	3,2,7,7,2,6,
	//z=1
	4,0,7,7,0,3,
	//z=-1
	1,5,2,2,5,6
    };

    vector<vec3> cube_vertices(cube_v, cube_v + 8);

    vec3 szd2 = sz / 2.0f;

    for (auto& v : cube_vertices)
	v *= szd2;

    for (int i = 0; i < 36; i += 3)
    {
	int i1 = cube_i[i];
	int i2 = cube_i[i + 1];
	int i3 = cube_i[i + 2];

	vec3& v1 = cube_vertices[i1];
	vec3& v2 = cube_vertices[i2];
	vec3& v3 = cube_vertices[i3];

	out_vertices.push_back(v1);
	out_vertices.push_back(v2);
	out_vertices.push_back(v3);
    }

    vector<vec3> v_n;
    v_n.reserve(8);

    vector<vec3> vertices(cube_v, cube_v + 8);
    vector<int> faces_idx(cube_i, cube_i + 36);

    calculate_vertex_normals(vertices, faces_idx, v_n);

    for (int i = 0; i < faces_idx.size(); i++)
    {
        out_normals.push_back(v_n[faces_idx[i]]);
    }
}

Vertex shader:

#version 330 core

layout(location = 0) in vec3 pos;
layout(location = 1) in vec3 n;

uniform vec3 lightPosition;

uniform mat4 model;
uniform mat4 view;
uniform mat4 proj;

out vec3 f_vertexPosition;
out vec3 f_lightPosition;
out vec3 f_normal;

void main()
{
	mat4 modelView = view * model;
	mat4 modelViewProj = proj * modelView;
	gl_Position = modelViewProj * vec4(pos, 1.0);

	mat3 normalMatrix = mat3(modelView);
	normalMatrix = inverse(normalMatrix);
	normalMatrix = transpose(normalMatrix);

	vec4 lpt = vec4(lightPosition,1.0f);
	
	lpt = lpt*modelView;
	
	f_vertexPosition = pos;
	f_normal = normalize(normalMatrix*n);
	f_lightPosition =  vec3(lpt);
}

Fragment shader:

#version 330 core

in vec3 f_vertexPosition;
in vec3 f_lightPosition;
in vec3 f_normal;

void main()
{
	vec3 faceColor = vec3(1.0, 0.0, 0.0);
	
	vec3 pl = f_vertexPosition-f_lightPosition;
	
	float dist = length(pl);
	
	vec3 lightDirection = normalize(pl);
	
	float magnitude = max(0.0, dot(f_normal, -lightDirection)) /(dist*dist); 
	
	gl_FragColor = vec4(faceColor*magnitude +  faceColor/3.0f, 1.0f);
}

The results I get at light position = (0.0f, 0.0f, 3.0f) by rotating box attached.

[attachment=33825:ex2.PNG]

[attachment=33826:ex3.PNG]

[attachment=33824:ex1.PNG]

I think the problem is at normal calculation step. But I cant figure out what's wrong...

Been there. You aren't normalizing your f_normal in the pixel shader. Remember, interpolating the normal between two directions can cause it to change magnitude.