Commit 9a09208c authored by COULAUD Olivier's avatar COULAUD Olivier

periodic improvement

parent bd6e637b
This diff is collapsed.
NaCl (10x10x10 unit cell) crystal
2 1 2000 -0.4269175E+06 3
2 1 2000 -0.4269175E+06
40.00
Na 1
-0.3909861360E-25 -0.3909861360E-25 -0.3909861360E-25
......
......@@ -382,17 +382,25 @@ The kernel should be able to proceed usual FMM operator in a tree of height of s
\section{Results}
\subsection{Ewald summation comparisons}
The energy computed by molecular dynamics codes is given by
The energy computed by molecular dynamics codes for a finite system is given by
$$
U = \frac{1}{4 \pi\epsilon_0}\sum_{i=0}^{N}{\sum_{j<i}{\frac{q_i q_j}{\|x_i-x_j\|}}}
$$
and the force on atom $x_i$
and the force on atom $x_i$ is
$$
f(x_i) = \frac{q_i }{4 \pi\epsilon_0}\sum_{j=0,i\neq j}^{N}{q_j\frac{x_i-x_j}{\|x_i-x_j\|^3}}
$$
Now consider a periodic system, the energy computed by molecular dynamics codes atypically by Ewald method is
$$
U = \frac{1}{8 \pi\epsilon_0}\sum_{|n|=0}^{\infty}{\sum_{i,j}^{N}{\sum_{j<i}{'\frac{q_i q_j}{\|x_i-x_j + n L \| }}}}
$$
where $L$ is the box size and ${'}$ specifies that for $n=0$ the omit i equal $j$ in the sum. The force on atom $x_i$ is
$$
f(x_i) = \frac{q_i}{8 \pi\epsilon_0}\sum_{|n|=0}^{\infty}{\sum_{i,j}^{N}{'q_j\frac{x_i-x_j}{\|x_i-x_j + n L\|^3}}}
$$
In the Ewald method, to original box is replicated up to the infinity, and then the method is equivalent to consider a surrounding conductor ($\epsilon=\infty$).
The Ewald method is equivalent to consider a surrounding conductor ($\epsilon=\infty$). The method described in this document and implemented in ScalFMM consider a finite set of the original box surrounded by a vacuum with dielectric constant $\epsilon=1$. The results by these two approaches are different \cite{DeLeeuw1980, Heyes1981} and for a very large sphere surrounding all boxes we have the following equation
The method described in this document and implemented in ScalFMM consider a finite set of the original box surrounded by a vacuum with dielectric constant $\epsilon=1$. The results by these two approaches are different \cite{DeLeeuw1980, Heyes1981} and for a very large sphere surrounding all boxes we have the following equation
\begin{equation}
E_{Ewald} = E_{FMM} -\frac{2\pi}{3V}|\sum_{i=1}^{N}{q_i x_i}|^2.
\end{equation}
......@@ -400,6 +408,20 @@ By tacking the gradient of this equation, the relationship between FMM forces an
\begin{equation}
F_{Ewald}(x_k) = -\nabla_k E_{Ewald} = F_{FMM}(x_k) + \frac{4\pi}{3V} q_k \sum_{i=1}^{N}{q_i x_i}.
\end{equation}
If we denote by $\mathbf{D}= \sum_{i=1}^{N}{q_i x_i}$ the dipole moment of the system then the periodic energy is
\begin{equation}
E_{Ewald} = E_{FMM} - \frac{2\pi}{3V}\mathbf{D}.\mathbf{D} =\sum_{i=1}^{N}{q_i V_i^{Ewald}}
\end{equation}
the potential
\begin{equation}
V^{Ewald}(r) = V^{FMM}(r) - \frac{4\pi}{3V} r.\mathbf{D}
\end{equation}
and the forces are
\begin{equation}
F_{Ewald}(x_k) = -\nabla_k E_{Ewald} = F_{FMM}(x_k) + \frac{4\pi}{3V} q_k \mathbf{D}.
\end{equation}
\subsubsection{DL\_Poly comparisons}
DL\_POLY\_2 uses the following internal molecular units \\
......@@ -502,14 +524,14 @@ struct MDParticle {
\end{verbatim}
We save the direct evaluation of the forces in a binary file. To read the file you should proceed as follow
\begin{verbatim}
std::ifstream fileout(filenameOut,
std::ifstream filein(filename,
std::ifstream::in| std::ios::binary);
file.seekg (std::ios::beg);
1 file.read( (char*)(&nbParticles), sizeof(int));
2 fileout.read( (char*)(&boxsize[0]), sizeof(double)*3);
3 fileout.read((char*)&PER[0],sizeof(int)*4);
4 fileout.read( (char*)(&denergy), sizeof(denergy));
filein.seekg (std::ios::beg);
1 filein.read( (char*)(&nbParticles), sizeof(int));
2 filein.read( (char*)(&boxsize[0]), sizeof(double)*3);
3 filein.read((char*)&PER[0],sizeof(int)*4);
4 filein.read( (char*)(&denergy), sizeof(denergy));
//
particles = new MDParticle[nbParticles*2];
fileout.read((char*)&particles[0], sizeof(MDParticle)*nbParticles);
......
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
......@@ -131,10 +131,12 @@ class TestSphericalDirectPeriodic : public FUTester<TestSphericalDirectPeriodic>
Print("Compute Diff...");
FMath::FAccurater potentialDiff;
FMath::FAccurater fx, fy, fz;
FReal energy = 0.0 ;
{ // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
......@@ -147,32 +149,39 @@ class TestSphericalDirectPeriodic : public FUTester<TestSphericalDirectPeriodic>
fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
energy += potentials[idxPart]*physicalValues[idxPart];
}
});
}
Print("Potential diff is = ");
Print(potentialDiff.getL2Norm());
Print(potentialDiff.getInfNorm());
printf(" L2Norm %e\n",potentialDiff.getRelativeL2Norm());
printf(" RMSError %e\n",potentialDiff.getRMSError());
Print("Fx diff is = ");
Print(fx.getL2Norm());
Print(fx.getInfNorm());
printf(" L2Norm %e\n",fx.getRelativeL2Norm());
printf(" RMSError %e\n",fx.getRMSError());
Print(fx.getRelativeL2Norm());
Print(fx.getRelativeInfNorm());
Print("Fy diff is = ");
Print(fy.getL2Norm());
Print(fy.getInfNorm());
printf(" L2Norm %e\n",fy.getRelativeL2Norm());
printf(" RMSError %e\n",fy.getRMSError());
Print("Fz diff is = ");
Print(fz.getL2Norm());
Print(fz.getInfNorm());
const FReal MaximumDiff = FReal(0.0001);
uassert(potentialDiff.getL2Norm() < MaximumDiff);
uassert(potentialDiff.getInfNorm() < MaximumDiff);
uassert(fx.getL2Norm() < MaximumDiff);
uassert(fx.getInfNorm() < MaximumDiff);
uassert(fy.getL2Norm() < MaximumDiff);
uassert(fy.getInfNorm() < MaximumDiff);
uassert(fz.getL2Norm() < MaximumDiff);
uassert(fz.getInfNorm() < MaximumDiff);
printf(" L2Norm %e\n",fz.getRelativeL2Norm());
printf(" RMSError %e\n",fz.getRMSError());
FReal L2error = (fx.getRelativeL2Norm()*fx.getRelativeL2Norm() + fy.getRelativeL2Norm()*fy.getRelativeL2Norm() + fz.getRelativeL2Norm() *fz.getRelativeL2Norm() );
printf(" L2 Force Error= %e\n",FMath::Sqrt(L2error)) ;
printf(" Energy = %.12e\n",energy);
const FReal MaximumDiff = FReal(0.001);
uassert(potentialDiff.getRelativeL2Norm() < MaximumDiff); // 1
uassert(potentialDiff.getRMSError() < MaximumDiff); // 2
uassert(fx.getRelativeL2Norm() < MaximumDiff); // 3
uassert(fx.getRMSError() < MaximumDiff); // 4
uassert(fy.getRelativeL2Norm() < MaximumDiff); // 5
uassert(fy.getRMSError() < MaximumDiff); // 6
uassert(fz.getRelativeL2Norm() < MaximumDiff); // 7
uassert(fz.getRMSError() < MaximumDiff); // 8
delete[] particles;
}
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment