Hello,
I am trying to do a molecular dynamics code and I am having some trouble with the precision. I am writing an Molecular Dynamics code. I am doing everything in double precision, but it seems the program gives non-zero output for configurations which should have 0.
I am testing my code with a 43 particles arranged in a hexagonal configuration. All of them are identical and the total force on the particle in the center should be zero.
For each pair the force is computed using the following code:
double dx,dy;
dx=xip-xjp;
dx=(dx<-llx*0.5?dx+llx:dx);
dx=(dx> llx*0.5?dx-llx:dx);
dy=yip-yjp[j];
dy=(dy<-lly*0.5?dy+lly:dy);
dy=(dy> lly*0.5?dy-lly:dy);
double rij2=dx*dx+dy*dy;
if(rij2<=9.0*hsq)
{
double rij=sqrt(rij2);
double q=rij/hg;
double dqkrij=quinticfirstderivative(q,cdqk);
double wp=cgd*(1.0/(rhoip)+1.0/rhojp)*dqkrij/rij;
xmyacc-=(wp*dx);
ymyacc-=(wp*dy);
}
For the central particle the code is executed for all jp indeces.
Because the configuration is symmetrical the total force should be 0, but it is not. It 1.0e-10. This is large enough to make problems in a energy conserverving algorithm algorithm like the velocity Verlet.
The constant cgd is about 1.0e+6, the rhoip and rhojp 1.0e+3, while the factor dqkrij/rij is about 15.
The function quinticfirstderivative(q,cdqk) is a forth order polynomial a1*(3-q)^4+a2*(2-q)^4+a3*(1-q)^4.
When I put all this Fortran I get a similar behavior, however the forces are a few order of magnitude smaller.
When error the calcualtion depends on the cgd constant. The error decreases if I decrease it.
Is this kind of error normal for CUDA?
Is there anything I can do to make it of the order 1.0e-14 without changing the constants?