How do I align the raytraced spheres from my fragment shader with GL_POINTS?

I have a very simple shader program that takes in a bunch of position data as GL_POINTS that generate screen-aligned squares of fragments like normal with a size depending on depth, and then in the fragment shader I wanted to draw a very simple ray-traced sphere for each one with just the shadow that is on the sphere opposite to the light. I went to this shadertoy to try to figure it out on my own. I used the sphIntersect function for ray-sphere intersection, and sphNormal to get the normal vectors on the sphere for lighting. The problem is that the spheres do not align with the squares of fragments, causing them to be cut off. This is because I am not sure how to match the projections of the spheres and the vertex positions so that they line up. Can I have an explanation of how to do this?

Here is a picture for reference.

Here are my vertex and fragment shaders for reference:

//vertex shader:
#version 460

layout(location = 0) in vec4 position; // position of each point in space
layout(location = 1) in vec4 color; //color of each point in space
layout(location = 2) uniform mat4 view_matrix; // projection * camera matrix
layout(location = 6) uniform mat4 cam_matrix; //just the camera matrix
out vec4 col; // color of vertex
out vec4 posi; // position of vertex

void main() {
    vec4 p = view_matrix * vec4(position.xyz, 1.0);
    gl_PointSize = clamp(1024.0 * position.w / p.z, 0.0, 4000.0);
    gl_Position = p;
    col = color;
    posi = cam_matrix * position;
}

//fragment shader:
#version 460

in vec4 col; // color of vertex associated with this fragment
in vec4 posi; // position of the vertex associated with this fragment relative to camera

out vec4 f_color;

layout (depth_less) out float gl_FragDepth;

float sphIntersect( in vec3 ro, in vec3 rd, in vec4 sph )
{
    vec3 oc = ro - sph.xyz;
    float b = dot( oc, rd );
    float c = dot( oc, oc ) - sph.w*sph.w;
    float h = b*b - c;
    if( h<0.0 ) return -1.0;
    return -b - sqrt( h );
}

vec3 sphNormal( in vec3 pos, in vec4 sph )
{
    return normalize(pos-sph.xyz);
}

void main() {

    vec4 c = clamp(col, 0.0, 1.0);
    vec2 p = ((2.0*gl_FragCoord.xy)-vec2(1920.0, 1080.0)) / 2.0;
    
    vec3 ro = vec3(0.0, 0.0, -960.0 );
    vec3 rd = normalize(vec3(p.x, p.y,960.0));
    
    vec3 lig = normalize(vec3(0.6,0.3,0.1));

    vec4 k = vec4(posi.x, posi.y, -posi.z, 2.0*posi.w);

    float t = sphIntersect(ro, rd, k);
    vec3 ps = ro + (t * rd);
    vec3 nor = sphNormal(ps, k);
    
    if(t < 0.0) c = vec4(1.0);
    
    else c.xyz *= clamp(dot(nor,lig), 0.0, 1.0);

    f_color = c;

    gl_FragDepth = t * 0.0001;

}

Solution

To create a ray direction that matches a perspective matrix from screen space, the following ray direction formula can be used:

vec3 rd = normalize(vec3(((2.0 / screenWidth) * gl_FragCoord.xy) - vec2(aspectRatio, 1.0), -proj_matrix[1][1]));

The value of 2.0 / screenWidth can be pre-computed or the opengl built-in uniform structs can be used.

To get a bounding box or other shape for your spheres, it is very important to use camera-facing shapes, and not camera-plane-facing shapes. Use the following process where position is the incoming VBO position data, and the w-component of position is the radius:

vec4 p = vec4((cam_matrix * vec4(position.xyz, 1.0)).xyz, position.w);
o.vpos = p;

float l2 = dot(p.xyz, p.xyz);
float r2 = p.w * p.w;
float k = 1.0 - (r2/l2);
float radius = p.w * sqrt(k);
if(l2 < r2) {
    p = vec4(0.0, 0.0, -p.w * 0.49, p.w);
    radius = p.w;
    k = 0.0;
}
vec3 hx = radius * normalize(vec3(-p.z, 0.0, p.x));
vec3 hy = radius * normalize(vec3(-p.x * p.y, p.z * p.z + p.x * p.x, -p.z * p.y));
p.xyz *= k;

Then use hx and hy as basis vectors for any 2D shape that you want the billboard to be shaped like for the vertices. Don't forget later to multiply each vertex by a perspective matrix to get the final position of each vertex. Here is a visualization of the billboarding on desmos using a hexagon shape: https://www.desmos.com/calculator/yeeew6tqwx