Range is working slow. So for solution on python requre more advanced approaches. Also you can reduce search diapason, enough iterate one coordinate quartal for each vertex.

Greetings. Your solution is 2000**2 + 2000*2 which is 8*10^6 operations. Python can hardly make more than 7*10^6 operations per second, so your solution must be just above the threshold of 1 second.

You can simply iterate over  < 0, a >  not  <  - a, a >  and if you get match, add (i, j) and (i,  - j) Also, count once and keep a * a, i * iandb * b in separate variables.