Calculating the potential of a matrix of size (424600, 55)
Your code loops N x M x N times, where N is currently 424600 and M is 55. I don't think it's going to finish. 

I tried a couple of variations of your code, and the results are below. The version marked "original" is pretty true to your original code, except that I returned the results instead of storing them in a file. (So I could compare them against the original version.)

The version marked "fortran" is an attempt to change the internal arrangement of data in the array, to try to get better cache performance out of consecutive reads. It didn't work, but I may not have done it right, or I maybe should have done it on the numpyier version. (It's late, I'm tired, sorry.)

The "cached" version is the original code, but I tried to cache everything I could to eliminate repeated reads of the same items. I also tried to eliminate repeated sub-expressions, and I eliminated the sigma variable, since it was 1 and never changed.

The "numpyier" version uses numpy operations in the inner loop. There may be a way to do the same thing with the columns, doing a 2-d computation and eliminating the middle loop as well. Someone better at numpy than I am may be able to provide that. At any rate, this version does not produce quite the same results, but they're close and it runs a lot faster. 

import sys
import csv
from timeit import timeit

import numpy as np
#from normalization import normalized_matrix
#MAKE_ROWS = 424600
MAKE_ROWS = 424
MAKE_COLS = 55
normalized_matrix = np.random.rand(MAKE_ROWS, MAKE_COLS)

def do_original():
    print(f"Original method", file=sys.stderr)
    n = len(normalized_matrix)
    sigma = 1
    #potential_writer = csv.writer(open('pot_orig.csv', 'w'))
    writer = []

    for l in range(n):
        elem = []
        for j in range(len(normalized_matrix[0])):
            GWF = 0
            PotSecPart = 0
            for i in range(n):
                if (i!=l):
                    PotSecPart = (PotSecPart + (((normalized_matrix[l][j] - normalized_matrix[i][j])**2)*(np.exp((-(normalized_matrix[l][j] - normalized_matrix[i][j])**2)/(2*(sigma**2))))))
                    GWF = GWF + (np.exp((-(normalized_matrix[l][j] - normalized_matrix[i][j])**2)/(2*(sigma**2))))
                    PotFirPart = (2*(sigma**2))* GWF
            elem.append(np.round(PotSecPart/PotFirPart, 2))
        #potential_writer.writerows([elem])
        writer.append(elem)
    return writer

def do_cached():
    print(f"Cached method", file=sys.stderr)
    num_rows = len(normalized_matrix)
    num_cols = len(normalized_matrix[0])

    sigma = 1
    #potential_writer = csv.writer(open('pot_orig.csv', 'w'))
    writer = []

    exp = np.exp

    for l in range(num_rows):
        elem = []
        for j in range(num_cols):
            GWF = 0
            PotSecPart = 0

            nm_lj = normalized_matrix[l][j]

            for i in range(num_rows):
                if i != l:
                    nm_ij = normalized_matrix[i][j]
                    diff_sq = (nm_lj - nm_ij) ** 2
                    exp_diff_sq_2 = exp(-diff_sq/2)

                    PotSecPart += diff_sq * exp_diff_sq_2
                    GWF += exp_diff_sq_2
                    PotFirPart = 2 * GWF
            elem.append(np.round(PotSecPart/PotFirPart, 2))
        #potential_writer.writerows([elem])
        writer.append(elem)
    return writer

def do_fortran():
    print(f"Fortran (cached) method", file=sys.stderr)
    matrix = normalized_matrix
    shape = matrix.shape
    matrix = np.reshape(matrix, shape, order='F')
    num_rows, num_cols = shape

    sigma = 1
    #potential_writer = csv.writer(open('pot_orig.csv', 'w'))
    writer = []

    exp = np.exp

    for l in range(num_rows):
        elem = []
        for j in range(num_cols):
            GWF = 0
            PotSecPart = 0

            nm_lj = matrix[l][j]

            for i in range(num_rows):
                if i != l:
                    nm_ij = matrix[i][j]
                    diff_sq = (nm_lj - nm_ij) ** 2
                    exp_diff_sq_2 = exp(-diff_sq/2)

                    PotSecPart += diff_sq * exp_diff_sq_2
                    GWF += exp_diff_sq_2
                    PotFirPart = 2 * GWF
            elem.append(np.round(PotSecPart/PotFirPart, 2))
        #potential_writer.writerows([elem])
        writer.append(elem)
    return writer

def do_numpyier():
    print(f"Numpy-ier method", file=sys.stderr)
    matrix = normalized_matrix
    num_rows, num_cols = normalized_matrix.shape

    sigma = 1
    #potential_writer = csv.writer(open('pot_orig.csv', 'w'))
    writer = []

    exp = np.exp

    for l in range(num_rows):
        elem = []
        for j in range(num_cols):
            GWF = 0
            PotSecPart = 0

            nm_lj = normalized_matrix[l][j]

            diff_sq = (nm_lj - matrix[:, j]) ** 2
            np.delete(diff_sq, l)
            exp_diff_sq_2 = exp(-diff_sq / 2)
            np.delete(exp_diff_sq_2, l)

            PotSecPart = np.sum(diff_sq * exp_diff_sq_2)
            GWF = np.sum(exp_diff_sq_2)
            PotFirPart = 2 * GWF

            elem.append(np.round(PotSecPart/PotFirPart, 2))
        #potential_writer.writerows([elem])
        writer.append(elem)
    return writer

if __name__ == '__main__':
    #orig = do_original()
    #print(timeit('do_original()', globals=globals(), number=1))
    #cached = do_cached()
    #assert cached == orig
    print(timeit('do_cached()', globals=globals(), number=1))
    #fortran = do_fortran()
    #assert fortran == orig
    #print(timeit('do_fortran()', globals=globals(), number=1))
    #numpyier = do_numpyier()
    #assert numpyier == cached
    print(timeit('do_numpyier()', globals=globals(), number=1))


The cached version was a little better than 2x faster than the original version, when I ran it on a 42x55 array. After a few passes, I stopped running the original.

The fortran version surprised me, because I expected to see better performance. Its results were consistently a smidge faster (.01-.02 sec) than the cached version, but that's it.

The numpyier version uses numpy operations to do the individual element calculations. This effectively replaces the inner loop with C. 

Please Note: The numpyier version gets different results than your python-only version. The differences are typically +/- 0.01 in your rounded output. But they're different, and I haven't looked into why.

I did most of my work at a shape of 42 x 55. That was in the 3-second range for me, about what I was willing to run repeatedly. After I got the numpyier version sort-of working, and satisfied myself that the differences were small, I ran it again with a shape of 424 x 55. That's 10 times bigger, which means 100x slower since this is an \$n^2\$ operation. The results were:

Cached method
114.699519728
Numpy-ier method
4.552203466999998


This seems intuitively correct: the numpy-ier version eliminates the python inner loop in favor of C looping operations, reducing the runtime by a factor of 424. The results are 23x faster, which is on the order of 100x. 

Note that this suggests that increasing the shape by a factor of 1000 (424600 instead of 424) would increase the runtime by somewhere between 1,000x and 1,000,000x. The numpyier version would require somewhere between 4000 seconds (1 hr, 20 mins) and 4,000,000 seconds (50 days).

I am working on a 10+ year old, 32-bit laptop. Your mileage is very likely going to be better.