Skip to content

GitLab

S
ScalFMM
Commit 4c09a276 authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files
Options

Improvemenst in Dlpoly comparaisons

parent cd7a1e72
Hide whitespace changes
Inline Side-by-side
Showing with 264 additions and 170 deletions
Doc/Src_tex/ScalFmm-PeriodicModel/periodicmodel.tex
View file @ 4c09a276
......@@ -448,9 +448,9 @@ The first test consists in a small crystal $4\times 4\times 4$ of NaCl. It is co
Crystal & Cube size & Cut off & reciprocal lattice vector &Energy \\
& $\AA$ & $\AA$ & & $kcal\, mol^{-1}$ \\
\hline
Small & 16 & 8 & (9,9,9) & $ -2.346403 \,10^{4} $ \\
Small & 16 & 8 & (19,19,19) & $ -2.346403 \,10^{4} $ \\
\hline
Big & 40 & 12 & (15,15,15) & $-3.666255 \,10^{5} $ \\
Big & 40 & 12 & (25,25,25) & $-3.666255 \,10^{5} $ \\
\hline
\end{tabular}
\end{center}
......
Src/Files/FDlpolyLoader.hpp
View file @ 4c09a276
......@@ -13,8 +13,8 @@
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
#ifndef FEWALLOADER_HPP
#define FEWALLOADER_HPP
#ifndef FDlpolyLOADER_HPP
#define FDlpolyLOADER_HPP
#include <iostream>
......@@ -25,22 +25,15 @@
#include "FAbstractLoader.hpp"
#include "../Utils/FPoint.hpp"
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class FEwalLoader
* Please read the license
* Particle has to extend {FExtendPhysicalValue,FExtendPosition}
*/
class FEwalLoader : public FAbstractLoader {
class FDlpolyLoader : public FAbstractLoader {
protected:
std::ifstream file; //< The file to read
// std::ifstream file; //< The file to read
FPoint centerOfBox; //< The center of box read from file
FReal boxWidth; //< the box width read from file
int nbParticles; //< the number of particles read from file
int levcfg ; //< DL_POLY CONFIG file key. 0,1 or 2
FReal energy ;
public:
protected:
enum Type{
OW,
HW,
......@@ -48,7 +41,113 @@ public:
Cl,
Undefined,
};
public:
virtual ~FDlpolyLoader()
{}
virtual bool isOpen() const =0;
virtual void fillParticle(FPoint* inPosition, FReal inForces[3], FReal* inPhysicalValue, int* inIndex)=0;
/**
* To get the number of particles from this loader
* @param the number of particles the loader can fill
*/
FSize getNumberOfParticles() const{
return FSize(this->nbParticles);
}
/**
* The center of the box from the simulation file opened by the loader
* @return box center
*/
FPoint getCenterOfBox() const{
return this->centerOfBox;
}
/**
* The box width from the simulation file opened by the loader
* @return box width
*/
FReal getBoxWidth() const{
return this->boxWidth;
}
FReal getEnergy() const{
return this->energy;
}
void getPhysicalValue(std::string &type,FReal & outPhysicalValue, int & outIndex) const{
if( strncmp(type.c_str(), "OW", 2) == 0){
outPhysicalValue = FReal(-0.82);
outIndex = OW;
}
else if( strncmp(type.c_str(), "HW", 2) == 0){
outPhysicalValue = FReal(-0.41);
outIndex = HW;
}
else if( strncmp(type.c_str(), "Na", 2) == 0){
outPhysicalValue = FReal(1.0);
outIndex = Na;
}
else if( strncmp(type.c_str(), "Cl", 2) == 0){
outPhysicalValue = FReal(-1.0);
outIndex = Cl;
}
else{
std::cerr << "Atom type not defined "<< type << std::endl;
exit(-1);
}
// std::cout << "Atom type " << type << " "<< outPhysicalValue << " " << outIndex <<std::endl;
}
void getPhysicalValue(char type[2],FReal & outPhysicalValue, int & outIndex) const{
if( strncmp(type, "OW", 2) == 0){
outPhysicalValue = FReal(-0.82);
outIndex = OW;
}
else if( strncmp(type, "HW", 2) == 0){
outPhysicalValue = FReal(-0.41);
outIndex = HW;
}
else if( strncmp(type, "Na", 2) == 0){
outPhysicalValue = FReal(1.0);
outIndex = Na;
}
else if( strncmp(type, "Cl", 2) == 0){
outPhysicalValue = FReal(-1.0);
outIndex = Cl;
}
else{
std::cerr << "Atom type not defined "<< type << std::endl;
exit(-1);
}
// std::cout << "Atom type " << type << " "<< outPhysicalValue << " " << outIndex <<std::endl;
}
} ;
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class FDlpolyLoader
* Please read the license
* Particle has to extend {FExtendPhysicalValue,FExtendPosition}
*/
class FDlpolyAsciiLoader : public FDlpolyLoader {
private:
std::ifstream file; //< The file to read
// FPoint centerOfBox; //< The center of box read from file
// FReal boxWidth; //< the box width read from file
// int nbParticles; //< the number of particles read from file
// int levcfg ; //< DL_POLY CONFIG file key. 0,1 or 2
// FReal energy ;
//public:
// enum Type{
// OW,
// HW,
// Na,
// Cl,
// Undefined,
// };
//
/**
* The constructor need the file name
* @param filename the name of the file to open
......@@ -59,13 +158,15 @@ public:
0.000000000000 17.200000000000 0.000000000000
0.000000000000 0.000000000000 17.200000000000
*/
FEwalLoader(const char* const filename): file(filename,std::ifstream::in){
public:
FDlpolyAsciiLoader(const char* const filename): file(filename,std::ifstream::in){
// test if open
if(this->file.is_open()){
const int bufferSize = 512;
char buffer[bufferSize];
file.getline(buffer, bufferSize);
int imcon ;
//int tempi(0);
FReal tempf(0);
......@@ -89,11 +190,13 @@ public:
this->boxWidth = 0;
this->nbParticles = 0;
}
std::cout << "ASCII LOADER "<< this->nbParticles<< " "<< this->boxWidth<<std::endl;
}
/**
* Default destructor, simply close the file
*/
virtual ~FEwalLoader(){
virtual ~FDlpolyAsciiLoader(){
file.close();
}
......@@ -105,33 +208,6 @@ public:
return this->file.is_open() && !this->file.eof();
}
/**
* To get the number of particles from this loader
* @param the number of particles the loader can fill
*/
FSize getNumberOfParticles() const{
return FSize(this->nbParticles);
}
/**
* The center of the box from the simulation file opened by the loader
* @return box center
*/
FPoint getCenterOfBox() const{
return this->centerOfBox;
}
/**
* The box width from the simulation file opened by the loader
* @return box width
*/
FReal getBoxWidth() const{
return this->boxWidth;
}
FReal getEnergy() const{
return this->energy;
}
/**
* Fill a particle
......@@ -145,13 +221,17 @@ public:
void fillParticle(FPoint* inPosition, FReal inForces[3], FReal* inPhysicalValue, int* inIndex){
FReal x, y, z, fx, fy, fz, vx, vy, vz;
int index;
char type[2];
std::string line;
file.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
file.read(type, 2);
this->file.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
std::string atomType, line;
this->file >> atomType;
//
file >> index;
std::getline(file, line); // needed to skip the end of the line in non periodic case
std::cout << "line: " << line << std::endl;
if ( levcfg == 0) {
file >> x >> y >> z;
}else if ( levcfg == 1) {
......@@ -170,39 +250,19 @@ public:
inForces[1] = fy;
inForces[2] = fz;
if( strncmp(type, "OW", 2) == 0){
*inPhysicalValue = FReal(-0.82);
*inIndex = OW;
}
else if( strncmp(type, "HW", 2) == 0){
*inPhysicalValue = FReal(-0.41);
*inIndex = HW;
}
else if( strncmp(type, "Na", 2) == 0){
*inPhysicalValue = FReal(1.0);
*inIndex = Na;
}
else if( strncmp(type, "Cl", 2) == 0){
*inPhysicalValue = FReal(-1.0);
*inIndex = Cl;
}
else{
*inPhysicalValue = FReal(-0.41);
*inIndex = HW;
std::cerr << "Atom type not defined "<< type << std::endl;
exit(-1);
}
getPhysicalValue(atomType, *inPhysicalValue,*inIndex );
}
};
class FEwalBinLoader : public FAbstractLoader {
class FDlpolyBinLoader : public FDlpolyLoader {
protected:
FILE* const file; //< The file to read
FPoint centerOfBox; //< The center of box read from file
double boxWidth; //< the box width read from file
int nbParticles; //< the number of particles read from file
double energy;
// FPoint centerOfBox; //< The center of box read from file
// double boxWidth; //< the box width read from file
// int nbParticles; //< the number of particles read from file
// double energy;
size_t removeWarning;
template<class Type>
......@@ -212,8 +272,8 @@ protected:
removeWarning = fread(&sizeBefore, sizeof(int), 1, file);
removeWarning = fread(&value, sizeof(Type), 1, file);
removeWarning = fread(&sizeAfter, sizeof(int), 1, file);
if( sizeBefore != sizeof(Type) ) printf("Error in loader ewal Size before %d should be %lu\n", sizeBefore, sizeof(Type));
if( sizeAfter != sizeof(Type) ) printf("Error in loader ewal Size after %d should be %lu\n", sizeAfter, sizeof(Type));
if( sizeBefore != sizeof(Type) ) printf("Error in loader Dlpoly Size before %d should be %lu\n", sizeBefore, sizeof(Type));
if( sizeAfter != sizeof(Type) ) printf("Error in loader Dlpoly Size after %d should be %lu\n", sizeAfter, sizeof(Type));
return value;
}
......@@ -223,8 +283,8 @@ protected:
removeWarning = fread(&sizeBefore, sizeof(int), 1, file);
removeWarning = fread(array, sizeof(Type), size, file);
removeWarning = fread(&sizeAfter, sizeof(int), 1, file);
if( sizeBefore != int(sizeof(Type) * size) ) printf("Error in loader ewal Size before %d should be %lu\n", sizeBefore, size*sizeof(Type));
if( sizeAfter != int(sizeof(Type) * size) ) printf("Error in loader ewal Size after %d should be %lu\n", sizeAfter, size*sizeof(Type));
if( sizeBefore != int(sizeof(Type) * size) ) printf("Error in loader Dlpoly Size before %d should be %lu\n", sizeBefore, size*sizeof(Type));
if( sizeAfter != int(sizeof(Type) * size) ) printf("Error in loader Dlpoly Size after %d should be %lu\n", sizeAfter, size*sizeof(Type));
return array;
}
......@@ -237,7 +297,7 @@ public:
[index charge x y z fx fy fz]
int double double ...
*/
FEwalBinLoader(const char* const filename): file(fopen(filename, "rb")) {
FDlpolyBinLoader(const char* const filename): file(fopen(filename, "rb")) {
// test if open
if(this->file != NULL){
energy = readValue<double>();
......@@ -255,7 +315,7 @@ public:
/**
* Default destructor, simply close the file
*/
virtual ~FEwalBinLoader(){
virtual ~FDlpolyBinLoader(){
fclose(file);
}
......@@ -267,33 +327,33 @@ public:
return this->file != NULL;
}
/**
* To get the number of particles from this loader
* @param the number of particles the loader can fill
*/
FSize getNumberOfParticles() const{
return FSize(this->nbParticles);
}
/**
* The center of the box from the simulation file opened by the loader
* @return box center
*/
FPoint getCenterOfBox() const{
return this->centerOfBox;
}
/**
* The box width from the simulation file opened by the loader
* @return box width
*/
FReal getBoxWidth() const{
return static_cast<FReal>(this->boxWidth);
}
FReal getEnergy() const{
return static_cast<FReal>(this->energy);
}
// /**
// * To get the number of particles from this loader
// * @param the number of particles the loader can fill
// */
// FSize getNumberOfParticles() const{
// return FSize(this->nbParticles);
// }
//
// /**
// * The center of the box from the simulation file opened by the loader
// * @return box center
// */
// FPoint getCenterOfBox() const{
// return this->centerOfBox;
// }
//
// /**
// * The box width from the simulation file opened by the loader
// * @return box width
// */
// FReal getBoxWidth() const{
// return static_cast<FReal>(this->boxWidth);
// }
//
// FReal getEnergy() const{
// return static_cast<FReal>(this->energy);
// }
/**
* Fill a particle
......@@ -307,7 +367,7 @@ public:
int size;
removeWarning = fread(&size, sizeof(int), 1, file);
if(size != 60) printf("Error in loader ewal Size %d should be %d\n", size, 60);
if(size != 60) printf("Error in loader Dlpoly Size %d should be %d\n", size, 60);
removeWarning = fread(&index, sizeof(int), 1, file);
removeWarning = fread(&charge, sizeof(double), 1, file);
......@@ -321,7 +381,7 @@ public:
removeWarning = fread(&fz, sizeof(double), 1, file);
removeWarning = fread(&size, sizeof(int), 1, file);
if(size != 60) printf("Error in loader ewal Size %d should be %d\n", size, 60);
if(size != 60) printf("Error in loader Dlpoly Size %d should be %d\n", size, 60);
inPosition->setPosition( static_cast<FReal>(x), static_cast<FReal>(y) ,static_cast<FReal>(z));
inForces[0] = static_cast<FReal>(fx);
......@@ -334,6 +394,6 @@ public:
};
#endif //FEwalLoader_HPP
#endif //FDlpolyLoader_HPP
Tests/Kernels/testSphericalDlpolyAlgorithm.cpp
View file @ 4c09a276
// ===================================================================================
// 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.
//
......@@ -35,7 +35,7 @@
#include "../../Src/Core/FFmmAlgorithm.hpp"
#include "../../Src/Files/FEwalLoader.hpp"
#include "../../Src/Files/FDlpolyLoader.hpp"
#include "../../Src/Components/FSimpleLeaf.hpp"
#include "../../Src/Kernels/P2P/FP2PParticleContainerIndexed.hpp"
......@@ -70,9 +70,9 @@ int main(int argc, char ** argv){
typedef FSimpleLeaf< ContainerClass > LeafClass;
#ifdef ScalFMM_USE_BLAS
// begin Chebyshef kernel
// begin Chebyshev kernel
// accuracy
const unsigned int ORDER = 8;
const unsigned int ORDER = 10;
// typedefs
typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FChebCell<ORDER> CellClass;
......@@ -106,16 +106,23 @@ int main(int argc, char ** argv){
const FReal r4pie0 = FReal(138935.4835);
FReal scaleEnergy, scaleForce ;
// -----------------------------------------------------
// LOADER
// -----------------------------------------------------
std::cout << "Opening : " << filename << "\n";
FEwalLoader loader(filename);
// FEwalBinLoader loader(filename);
if(!loader.isOpen()){
FDlpolyLoader *loader = nullptr ;
if(FParameters::existParameter(argc, argv, "-bin")){
loader = new FDlpolyBinLoader(filename);
}
else {
loader = new FDlpolyAsciiLoader(filename);
}
if(! loader->isOpen()){
std::cout << "Loader Error, " << filename << " is missing\n";
return 1;
}
// ---------------------------------------------------
// DL_POLY
// DL_POLY CONSTANT
// ---------------------------------------------------
scaleEnergy = r4pie0 / 418.4 ; // kcal mol^{-1}
scaleForce = -r4pie0 ; // 10 J mol^{-1} A^{-1}
......@@ -127,22 +134,24 @@ int main(int argc, char ** argv){
#else
std::cout << " $$$$$$$$$$$$$$$ CHEBYCHEV VERSION $$$$$$$$$$$$" <<std::endl;
#endif
OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
OctreeClass tree(NbLevels, SizeSubLevels, loader->getBoxWidth(), loader->getCenterOfBox());
// ---------------------------------------------------------------------------------
// Insert particles in the Octree
// --------------------------------------------------------------------------------- std::cout << "Creating & Inserting " << loader.getNumberOfParticles() << " particles ..." << std::endl;
std::cout << "\tWidth : " << loader.getBoxWidth() << " \t center x : " << loader.getCenterOfBox().getX()
<< " y : " << loader.getCenterOfBox().getY() << " z : " << loader.getCenterOfBox().getZ() << std::endl;
// --------------------------------------------------------------------------------- std::cout << "Creating & Inserting " << loader->getNumberOfParticles() << " particles ..." << std::endl;
std::cout << "\tWidth : " << loader->getBoxWidth() << " \t center x : " << loader->getCenterOfBox().getX()
<< " y : " << loader->getCenterOfBox().getY() << " z : " << loader->getCenterOfBox().getZ() << std::endl;
std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
counter.tic();
FPoint electricMoment(0.0,0.0,0.0) ;
EwalParticle * const particles = new EwalParticle[loader.getNumberOfParticles()];
memset(particles, 0, sizeof(EwalParticle) * loader.getNumberOfParticles());
EwalParticle * const particles = new EwalParticle[loader->getNumberOfParticles()];
memset(particles, 0, sizeof(EwalParticle) * loader->getNumberOfParticles());
double totalCharge = 0.0;
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
loader.fillParticle(&particles[idxPart].position, particles[idxPart].forces,
for(int idxPart = 0 ; idxPart < loader->getNumberOfParticles() ; ++idxPart){
//
loader->fillParticle(&particles[idxPart].position, particles[idxPart].forces,
&particles[idxPart].physicalValue,&particles[idxPart].index);
//
totalCharge += particles[idxPart].physicalValue ;
electricMoment.incX(particles[idxPart].physicalValue*particles[idxPart].position.getX() );
electricMoment.incY(particles[idxPart].physicalValue*particles[idxPart].position.getY() );
......@@ -153,13 +162,13 @@ int main(int argc, char ** argv){
counter.tac();
double dipoleNorm = electricMoment.norm2() ;
double volume = loader.getBoxWidth()*loader.getBoxWidth()*loader.getBoxWidth() ;
double volume = loader->getBoxWidth()*loader->getBoxWidth()*loader->getBoxWidth() ;
double coeffCorrectionDLPOLY = 2.0*FMath::FPi/volume/3.0 ;
std::cout << std::endl;
std::cout << "Total Charge = "<< totalCharge <<std::endl;
std::cout << "Electric Moment = "<< electricMoment <<std::endl;
std::cout << "electric Moment = "<< dipoleNorm <<std::endl;
std::cout << "Total Charge = "<< totalCharge <<std::endl;
std::cout << "Electric Moment = "<< electricMoment <<std::endl;
std::cout << "Electric Moment norm = "<< dipoleNorm <<std::endl;
std::cout << std::endl;
std::cout << std::endl;
......@@ -172,8 +181,8 @@ int main(int argc, char ** argv){
typename OctreeClass::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
int inTreeHeight = NbLevels ;
double inBoxWidth = loader.getBoxWidth() ;
FPoint inBoxCenter(loader.getCenterOfBox()) ;
double inBoxWidth = loader->getBoxWidth() ;
FPoint inBoxCenter(loader->getCenterOfBox()) ;
double widthAtLeafLevel(inBoxWidth/FReal(1 << (inTreeHeight-1))) , widthAtLeafLevelDiv2 = widthAtLeafLevel/2;
FPoint boxCorner(inBoxCenter.getX()-(inBoxWidth/2),inBoxCenter.getY()-(inBoxWidth/2),
inBoxCenter.getZ()-(inBoxWidth/2));
......@@ -205,9 +214,9 @@ int main(int argc, char ** argv){
if( FParameters::existParameter(argc, argv, "-noper") ){
#ifndef ScalFMM_USE_BLAS
KernelClass kernels( DevP, NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
KernelClass kernels( DevP, NbLevels, loader->getBoxWidth(), loader->getCenterOfBox());
#else
KernelClass kernels( NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
KernelClass kernels( NbLevels, loader->getBoxWidth(), loader->getCenterOfBox());
#endif
FmmClassNoPer algo(&tree,&kernels);
algo.execute();
......@@ -249,7 +258,7 @@ int main(int argc, char ** argv){
printf("Box [%d;%d][%d;%d][%d;%d]\n", min.getX(), max.getX(), min.getY(),
max.getY(), min.getZ(), max.getZ());
particlesDirect = new EwalParticle[loader.getNumberOfParticles()];
particlesDirect = new EwalParticle[loader->getNumberOfParticles()];
FReal denergy = 0.0;
FMath::FAccurater dfx, dfy, dfz ;
......@@ -259,7 +268,7 @@ int main(int argc, char ** argv){
// particlesDirect : Direct results
//
counter.tic();
for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
for(int idxTarget = 0 ; idxTarget < loader->getNumberOfParticles() ; ++idxTarget){
particlesDirect[idxTarget] = particles[idxTarget];
EwalParticle & part = particlesDirect[idxTarget];
part.forces[0] = part.forces[1] = part.forces[2] = 0.0;
......@@ -269,7 +278,7 @@ int main(int argc, char ** argv){
//
// Compute force and potential between particles[idxTarget] and particles inside the box
//
for(int idxOther = 0; idxOther < loader.getNumberOfParticles() ; ++idxOther){
for(int idxOther = 0; idxOther < loader->getNumberOfParticles() ; ++idxOther){
if( idxOther != idxTarget ){
FP2P::NonMutualParticles(
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
......@@ -286,19 +295,19 @@ int main(int argc, char ** argv){
for(int idxX = min.getX() ; idxX <= max.getX() ; ++idxX){
for(int idxY = min.getY() ; idxY <= max.getY() ; ++idxY){
for(int idxZ = min.getZ() ; idxZ <= max.getZ() ; ++idxZ){
int nL = PeriodicDeep ;
{
//
//int nL = PeriodicDeep ;
// for(int idxX = -nL ; idxX <= nL -1; ++idxX){
// for(int idxY = -nL ; idxY <=nL-1 ; ++idxY){
// for(int idxZ = -nL ; idxZ <= nL-1; ++idxZ){
if(idxX == 0 && idxY == 0 && idxZ == 0) continue;
const FPoint offset(loader.getBoxWidth() * FReal(idxX),
loader.getBoxWidth() * FReal(idxY),
loader.getBoxWidth() * FReal(idxZ));
const FPoint offset(loader->getBoxWidth() * FReal(idxX),
loader->getBoxWidth() * FReal(idxY),
loader->getBoxWidth() * FReal(idxZ));
// std::cout <<" ( "<< idxX<<" , "<<idxY << " , "<< idxZ << " ) "<< offset <<std::endl;
for(int idxSource = 0 ; idxSource < loader.getNumberOfParticles() ; ++idxSource){
for(int idxSource = 0 ; idxSource < loader->getNumberOfParticles() ; ++idxSource){
EwalParticle source = particles[idxSource];
source.position += offset;
// std::cout << "Part "<<idxSource<< " " <<source.position.getX()<< " " << source.position.getY()<< " " <<source.position.getZ()<< " " <<source.physicalValue <<std::endl ;
......@@ -360,10 +369,10 @@ int main(int argc, char ** argv){
printf("DirectEwald InfNorm %e\n",dfz.getRelativeInfNorm());
double L2error = (dfx.getRelativeL2Norm()*dfx.getRelativeL2Norm() + dfy.getRelativeL2Norm()*dfy.getRelativeL2Norm() + dfz.getRelativeL2Norm() *dfz.getRelativeL2Norm() );
printf("DirectEwald L2 Force Error= %e\n",FMath::Sqrt(L2error)) ;
printf("DirectEwald Energy DIRECT= %e\n",denergy);
printf("DirectEwald Energy EWALD = %e\n",loader.getEnergy());
printf("DirectEwald Energy EWALD - Energy DIRECT = %e\n", loader.getEnergy()-denergy );
printf("DirectEwald |Energy EWALD - Energy DIRECT|/directEnergy= %e\n",FMath::Abs((loader.getEnergy()-denergy)/loader.getEnergy()));
printf("DirectEwald Energy DIRECT= %.12e\n",denergy);
printf("DirectEwald Energy EWALD = %.12e\n",loader->getEnergy());
printf("DirectEwald Energy EWALD - Energy DIRECT = %e\n", loader->getEnergy()-denergy );
printf("DirectEwald |Energy EWALD - Energy DIRECT|/directEnergy= %e\n",FMath::Abs((loader->getEnergy()-denergy)/loader->getEnergy()));
//
if(FParameters::existParameter(argc, argv, "-saveError")){
......@@ -383,8 +392,8 @@ int main(int argc, char ** argv){
// particlesDirect : Direct results
//
FReal potential = 0.0, energy = 0.0;
FMath::FAccurater fx, fy, fz, fmmdfx, fmmdfy, fmmdfz;
FReal energy = 0.0;
FMath::FAccurater fx, fy, fz, fmmdfx, fmmdfy, fmmdfz,fmmpot;
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const potentials = leaf->getTargets()->getPotentials();
......@@ -418,6 +427,7 @@ int main(int argc, char ** argv){
fmmdfx.add(particlesDirect[indexPartOrig].forces[0],forcesX[idxPart]);
fmmdfy.add(particlesDirect[indexPartOrig].forces[1],forcesY[idxPart]);
fmmdfz.add(pa