Skip to content
GitLab
Commit 217c8c55 authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files

BPC: Check with Dlpoly works not yet with Stamp

parent 3349dbe1
Hide whitespace changes
Inline Side-by-side
Tests/noDist/ChebyshevPeriodic.cpp
View file @ 217c8c55
......@@ -87,10 +87,11 @@ int main(int argc, char* argv[])
#endif
//
std::cout << "Parameters "<< std::endl
<< " Octree Depth "<< TreeHeight <<std::endl
<< " SubOctree depth " << SubTreeHeight <<std::endl
<< " Input file name: " <<filename <<std::endl
<< " Thread number: " << NbThreads <<std::endl
<< "\t Octree Depth \t"<< TreeHeight <<std::endl
<< "\t SubOctree depth \t" << SubTreeHeight <<std::endl
<< "\t Periodic depth \t" << PeriodicDeep <<std::endl
<< "\t Input file name: \t" <<filename <<std::endl
<< "\t Thread number: \t " << NbThreads <<std::endl
<<std::endl;
//
// init timer
......@@ -128,7 +129,6 @@ int main(int argc, char* argv[])
{ // -----------------------------------------------------
std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
<< " particles ..." << std::endl;
std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
time.tic();
//
FPoint position;
......
Tests/noDist/compare2Files.cpp → Utils/noDist/compare2Files.cpp
View file @ 217c8c55
......@@ -16,7 +16,7 @@
#include "Files/FFmaGenericLoader.hpp"
#include "Utils/FParameters.hpp"
#include "Utils/FCompareResults.hpp"
#include "../../Src/Utils/FParameterNames.hpp"
#include "Utils/FParameterNames.hpp"
//
/// \file compare2files.cpp
//!
......@@ -58,7 +58,7 @@ int main(int argc, char ** argv){
//
// Allocation
//
FSize nbParticles = loader1.getNumberOfParticles();
FSize nbParticles = loader1.getNumberOfParticles();
const unsigned int nbData = loader1.getNbRecordPerline() ;
if(nbParticles != loader2.getNumberOfParticles()){
std::cerr << "Number of points is different in the two files."<<std::endl ;
......@@ -75,13 +75,8 @@ int main(int argc, char ** argv){
//
loader1.fillParticle(particles1,nbParticles);
loader2.fillParticle(particles2,nbParticles);
if(FParameters::existParameter(argc, argv, LocalParameterEwald.options) ) {
FReal volume =1.0 ;
// double volume = boxsize[0] *boxsize[1]*boxsize[2] ;
removeFirstMoment( "DLPOLY", nbParticles, particles2, volume) ;
FPoint FirstMoment ;
}
//
const int error = compareTwoArrays("TAG", nbParticles, particles1, particles2);
//
......
Utils/noDist/removeMoment.cpp
View file @ 217c8c55
......@@ -49,6 +49,7 @@ void genDistusage() {
<< " -fin fileIn.{bin, txt) input file name "<< std::endl
<< " -fout fileout.{bfma, fma) output file name "<< std::endl
<< " -Ewald2FMM to add the Electric Moment in the potential and the force " << std::endl
<< " -energy Value of the energy (Ewald)" << std::endl
<< " -FMM2Ewald to remove the Electric Moment in the potential and the force " << std::endl ;
std::cout << " -verbose : print index x y z fx fy fy Q and V" << std::endl;
exit(-1);
......@@ -56,12 +57,15 @@ void genDistusage() {
// Simply create particles and try the kernels
int main(int argc, char ** argv){
//
bool useDLPOLY=false,useFMA=false,useSTAMP=false;
//
FReal scaleEnergy, scaleForce ;
const FReal r4pie0 = FReal(138935.4835);
const double q0 = 1.6021892e-19;
const double e0 = 8.854187e-12;
const double ang = 1.0e-10;
const double Unsur4pie0 =8.98755179e+9; // 1.0/(4*FMath::FPi*e0); //8.98755179e+9 ;
double D0 = 1.0;
std::cout << "Unsur4pie0: " << Unsur4pie0 << " 8.98755179e+9 Diff " << Unsur4pie0-8.98755179e+9 <<std::endl;
......@@ -86,11 +90,15 @@ int main(int argc, char ** argv){
if(FParameters::existParameter(argc, argv, "-dlpoly")){
// scale = false ;
scaleEnergy = r4pie0 / 418.4 ; // kcal mol^{-1}
scaleForce = -r4pie0 ; // 10 J mol^{-1} A^{-1}
scaleForce = -r4pie0 ; // 10 J mol^{-1} A^{-1}
D0 = 1.0;
useDLPOLY = true ;
}
else if(FParameters::existParameter(argc, argv, "-stamp")){
scaleEnergy = q0*q0*Unsur4pie0 /ang;
scaleForce = -scaleEnergy/ang;
D0 = q0/ang/ang;
useSTAMP = true;
}
std::cout <<"Scaling factor " <<std::endl
<< " Energy factor " << scaleEnergy<<std::endl
......@@ -106,20 +114,18 @@ int main(int argc, char ** argv){
std::cout << "Opening : " << filenameIn << "\n";
FDlpolyLoader *loaderDlpoly = nullptr ;
FFmaGenericLoader *loaderFMA = nullptr ;
bool useDLPOLY=false,useFMA=false;
std::string ext(".bin");
std::string ext1(".fma"),ext2(".bfma");
// open particle file
if(filenameIn.find(ext) != std::string::npos) {
loaderDlpoly = new FDlpolyBinLoader(filenameIn.c_str());
useDLPOLY = true ;
}
else if(filenameIn.find(".txt")!=std::string::npos ) {
loaderDlpoly = new FDlpolyAsciiLoader(filenameIn.c_str());
useDLPOLY = true ;
}
else if(filenameIn.find(ext1)!=std::string::npos ) {
std::cout << "FMA read " << std::endl;
loaderFMA = new FFmaGenericLoader(filenameIn.c_str());
useFMA = true ;
} else if(filenameIn.find(ext2)!=std::string::npos ) {
......@@ -133,14 +139,18 @@ int main(int argc, char ** argv){
double boxsize[3] ;
FPoint centre ;
FSize numberofParticles ;
double energy = 0.0 ;
FmaRWParticle<8,8>* particlesIn ; //= new FmaRWParticle<8,8>[numberofParticles];
if( FParameters::existParameter(argc, argv, "-energy") ){
energy = FParameters::getValue(argc,argv,"-energy", 0.0);
if (useDLPOLY) {
}
if (loaderDlpoly) {
// -----------------------------------------------------
// LOADER EWALD PARTICLES FROM DLPOLY
// -----------------------------------------------------
if(! loaderDlpoly->isOpen()){
std::cout << "Loader Error, " << filenameIn << " is missing\n";
std::cout << "DLPOLYLoader Error, " << filenameIn << " is missing\n";
exit(EXIT_FAILURE);
}
boxsize[0] = boxsize[1]= boxsize[2]=loaderDlpoly->getBoxWidth() ;
......@@ -157,6 +167,7 @@ int main(int argc, char ** argv){
FReal potential=0.0;
int index ;
loaderDlpoly->fillParticle(&position, forces,&physicalValue,&index);
index = index-1 ;
particlesIn[index].setPosition(position) ;
particlesIn[index].setPhysicalValue(physicalValue) ;
particlesIn[index].setPotential(potential) ;
......@@ -164,55 +175,51 @@ int main(int argc, char ** argv){
}
}
else if (useFMA) {
// -----------------------------------------------------
// -----------------------------------------------------
// LOADER EWALD PARTICLES FROM STAMP
// -----------------------------------------------------
if(! loaderFMA->isOpen()){
std::cout << "Loader Error, " << filenameIn << " is missing\n";
std::cout << "FMALoader Error, " << filenameIn << " is missing\n";
exit(EXIT_FAILURE);
}
boxsize[0] = boxsize[1]= boxsize[2] = loaderFMA->getBoxWidth()/ang;
boxsize[0] = boxsize[1]= boxsize[2] = loaderFMA->getBoxWidth();
// scaleEnergy /= boxsize[0] ;
// scaleForce = -scaleEnergy/boxsize[0] ;
// boxsize[0] = boxsize[1]= boxsize[2] = 1.0 ;
numberofParticles = loaderFMA->getNumberOfParticles() ;
particlesIn = new FmaRWParticle<8,8>[numberofParticles];
memset(particlesIn, 0, sizeof( FmaRWParticle<8,8>) * numberofParticles) ;
numberofParticles = loaderFMA->getNumberOfParticles() ;
std::cout << " numberofParticles "<< std::endl;
particlesIn = new FmaRWParticle<8,8>[numberofParticles];
memset(particlesIn, 0, sizeof( FmaRWParticle<8,8>) * numberofParticles) ;
loaderFMA->fillParticle(particlesIn,numberofParticles);
centre = loaderFMA->getCenterOfBox();
centre *= boxsize[0];
for(int idxPart = 0 ; idxPart < numberofParticles; ++idxPart){
//
FPoint PP = particlesIn[idxPart].getPosition() ;
PP *= boxsize[0];
particlesIn[idxPart].setPosition(PP) ;
//
}
}
counter.tic();
FPoint electricMoment(0.0,0.0,0.0) ;
// const --> then shared
//
double totalCharge = 0.0;
double totalCharge = 0.0, readEnergy;
//
for(int idxPart = 0 ; idxPart < numberofParticles; ++idxPart){
//
totalCharge += particlesIn[idxPart].getPhysicalValue() ;
readEnergy += particlesIn[idxPart].getPhysicalValue() *particlesIn[idxPart].getPotential() ;
//
electricMoment.incX(particlesIn[idxPart].getPhysicalValue()*particlesIn[idxPart].getPosition().getX() );
electricMoment.incY(particlesIn[idxPart].getPhysicalValue()*particlesIn[idxPart].getPosition().getY() );
electricMoment.incZ(particlesIn[idxPart].getPhysicalValue()*particlesIn[idxPart].getPosition().getZ() );
}
readEnergy *= 0.5 ;
counter.tac();
std::cout << std::endl;
std::cout << "Total Charge = "<< totalCharge <<std::endl;
std::cout << "Energy = "<< readEnergy <<std::endl;
std::cout << "Electric Moment = "<< electricMoment <<std::endl;
std::cout << "Electric Moment norm = "<< electricMoment.norm2() <<std::endl;
std::cout << "----------------------------------------------------"<<std::endl;
......@@ -230,16 +237,17 @@ int main(int argc, char ** argv){
double preScaleEnergy = 1.0, postScaleEnergy = 1.0, preScaleForce = 1.0, postScaleForce = 1.0 ;
if( FParameters::existParameter(argc, argv, "-Ewald2FMM") ){
preScaleEnergy = 1.0/scaleEnergy ; postScaleEnergy = 1.0 ; preScaleForce = 1.0/scaleForce; postScaleForce = 1.0 ;
coeffCorrection = -coeffCorrection;
}
else if( FParameters::existParameter(argc, argv, "-FMM2Ewald") ){
preScaleEnergy = 1.0 ; postScaleEnergy = scaleEnergy ; preScaleForce = 1.0/scaleForce; postScaleForce = scaleForce ;
coeffCorrection = -coeffCorrection;
preScaleEnergy = 1.0 ; postScaleEnergy = scaleEnergy ; preScaleForce = 1.0; postScaleForce = scaleForce ;
coeffCorrection = coeffCorrection;
}
else {
std::cout << " -Ewald2FMM ou -FMM2Ewald should be set"<<std::endl;
std::cout << " -Ewald2FMM or -FMM2Ewald should be set"<<std::endl;
exit(EXIT_FAILURE);
}
std::cout << "coeffCorrection: "<<coeffCorrection<<std::endl;
std::cout << "coeffCorrection: "<<coeffCorrection<<" " <<coeffCorrection*D0<< std::endl;
std::cout << "preScaleEnergy: "<<preScaleEnergy<<std::endl;
std::cout << "postScaleEnergy: "<<postScaleEnergy<<std::endl;
std::cout << "preScaleForce: "<<preScaleForce<<std::endl;
......@@ -249,15 +257,18 @@ int main(int argc, char ** argv){
for(int idx = 0 ; idx < numberofParticles ; ++idx){
tmp = particlesIn[idx].getPosition().getX()*electricMoment.getX() + particlesIn[idx].getPosition().getY()*electricMoment.getY()
+ particlesIn[idx].getPosition().getZ()*electricMoment.getZ() ;
//
// Potential rescaling ( Not checked)
double P = particlesIn[idx].getPotential();
particlesIn[idx].setPotential( P + 2.0*tmp*coeffCorrection);
particlesIn[idx].setPotential( P + tmp*coeffCorrection);
//
// Force rescaling
double * forces ;
forces = particlesIn[idx].getForces() ;
forces[0] = preScaleForce*forces[0] - coeffCorrection*particlesIn[idx].getPhysicalValue()*electricMoment.getX() ;
forces[1] = preScaleForce*forces[1] - coeffCorrection*particlesIn[idx].getPhysicalValue()*electricMoment.getY() ;
forces[2] = preScaleForce*forces[2] - coeffCorrection*particlesIn[idx].getPhysicalValue()*electricMoment.getZ() ;
double coef = coeffCorrection*particlesIn[idx].getPhysicalValue() ;
//
forces[0] = preScaleForce*forces[0] - coef*electricMoment.getX() ;
forces[1] = preScaleForce*forces[1] - coef*electricMoment.getY() ;
forces[2] = preScaleForce*forces[2] - coef*electricMoment.getZ() ;
//
forces[0] *= postScaleForce;
forces[1] *= postScaleForce;
......@@ -270,17 +281,16 @@ int main(int argc, char ** argv){
std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
std::cout << std::scientific;
std::cout.precision(10) ;
if (useDLPOLY){
std::cout << " " << coeffCorrection*electricMoment.norm2() <<std::endl;
std::cout << " " << preScaleEnergy*loaderDlpoly->getEnergy() <<std::endl;
newEnergy = preScaleEnergy*loaderDlpoly->getEnergy() + coeffCorrection*electricMoment.norm2() ;
newEnergy *= postScaleEnergy ;
//
std::cout << "Energy EWALD = "<< loaderDlpoly->getEnergy() <<std::endl ;
//
if (loaderDlpoly){
energy = loaderDlpoly->getEnergy() ;
}
std::cout << "Energy New = "<<newEnergy<<std::endl ;
std::cout <<std::endl;
std::cout << " Ener --> " <<preScaleEnergy <<" "<< readEnergy*preScaleEnergy<< " " << coeffCorrection*electricMoment.norm2() <<std::endl ;
//
newEnergy = preScaleEnergy*readEnergy + 0.5*coeffCorrection*electricMoment.norm2() ;
newEnergy *= postScaleEnergy ;
std::cout << " Cmp Energy " << newEnergy << " " << energy << " " <<FMath::Abs(newEnergy -energy ) << std::endl;
//
// Save data in FMAFormat
//
......
Supports Markdown
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