Clear["Global`*"];
listofplots = {};
n = 50; a =1.5;
list = {}; listdens = {}; d2 = .5; radius = 1;
f[r1_,r2_] := 1/Norm[r1-r2];
g[x_,y_,z_] = x^2+y^2+z^2;
U[x_,y_,z_] := Sum[f[points[[i]],{x,y,z}],{i,1,n}];
(*points = Table[{RandomReal[{1-.4,1+.4}],RandomReal[{Pi/2-.1,Pi/2+.1}],RandomReal[{0,Pi}]},{i,1,n}];*)
points = Table[{RandomReal[{1-.4,1+.4}],RandomReal[{Pi/2-.1,Pi/2+.1}],RandomReal[{0,Pi}]},{i,1,n}];
points = Map[CoordinateTransformData["Spherical"->"Cartesian","Mapping"],points]//N;
plot1= ListPointPlot3D[points,{PlotRange->{{-a,a},{-a,a},{-a,a}},BoxRatios->{1,1,1},PlotStyle->PointSize[0.02]},Boxed->False,Axes->False];
plot2 = ContourPlot3D[{g[x,y,z]==d2^2},{x,-a,a},{y,-a,a},{z,-a,a},ContourStyle->Directive[Opacity[0.4],Specularity[White,10],Hue[.5,.5,1]],Mesh->None,Boxed->False,Axes->False];
before =Show[plot1,plot2]

time1 = AbsoluteTime[];
iterations = 400; 
small = 0.005; large = 10;
For[i=1, i<iterations, i=i+1, 

For[k=1,k<n+1,k=k+1,
dist = large; index =0;
For[j=1,j<n+1,j=j+1,
If[k==j,Continue[],
If[dist > EuclideanDistance[points[[k]],points[[j]]],dist = EuclideanDistance[points[[k]],points[[j]]]; index = j]
];
];
displace = small*(points[[k]]-points[[index]])/Norm[points[[k]]-points[[index]]];

np = points[[k]]+displace;
condition1 = g[np[[1]],np[[2]],np[[3]]]> radius;
condition2 = g[np[[1]],np[[2]],np[[3]]]< d2;

If[condition1 || condition2, 
perp = Grad[g[x,y,z],{x,y,z}]/.{x->points[[k,1]],y->points[[k,2]],z->points[[k,3]]};
displace = displace - Dot[displace,perp]*perp/Norm[perp];
points[[k]]=points[[k]]+displace;
(*displace = small*(points[[k]]-points[[index]])/Norm[points[[k]]-points[[index]]];
displace = displace +(radius-Norm[points[[k]]])*perp/Norm[perp];
points[[k]]=points[[k]]+displace;*)
];

If[!(condition1 && condition2),
points[[k]]=points[[k]]+displace;
];
];


If[Mod[i,5]==0,plot1= ListPointPlot3D[points,{PlotRange->{{-a,a},{-a,a},{-a,a}},BoxRatios->{1,1,1},PlotStyle->PointSize[0.02]},Boxed->False,Axes->False];
after =Show[plot1,plot2];
list=Append[list,after];
];
];



plot1= ListPointPlot3D[points,{PlotRange->{{-a,a},{-a,a},{-a,a}},BoxRatios->{1,1,1},PlotStyle->PointSize[0.02]},Boxed->False,Axes->False];
time2=AbsoluteTime[];
Print[time2-time1];
after =Show[plot1,plot2]

ListAnimate[list]

Export["torus.gif",list]