diff --git a/Data/periodicCube20000.bfma b/Data/periodicCube20000.bfma new file mode 100644 index 0000000000000000000000000000000000000000..e4d65cdec1e22513d8b8159020be7c49d33bd818 Binary files /dev/null and b/Data/periodicCube20000.bfma differ diff --git a/Src/Files/FFmaGenericLoader.hpp b/Src/Files/FFmaGenericLoader.hpp index 6df9315833226025211f690247f75d32754b5384..1d73eef20e8c3790b97f394e42a738548d91c9c0 100755 --- a/Src/Files/FFmaGenericLoader.hpp +++ b/Src/Files/FFmaGenericLoader.hpp @@ -507,7 +507,7 @@ public: << " asked in structure " <<(*dataToRead).getReadDataNumber() <read((char*)((*dataToRead).getPtrFirstData()), sizeof(FReal)*(N*(*dataToRead).getReadDataNumber())); } diff --git a/Utils/noDist/removeMoment.cpp b/Utils/noDist/removeMoment.cpp index e5f82aeeca9319bb68e18530e83d45713923d976..e0e432e0c003acf0d68a8d0a8281f7a8fbaae7f0 100755 --- a/Utils/noDist/removeMoment.cpp +++ b/Utils/noDist/removeMoment.cpp @@ -42,30 +42,35 @@ /** DLpoly particle is used in the gadget program * here we try to make the same simulation */ -struct MDParticle { - FPoint position; - FReal forces[3]; - FReal physicalValue; - FReal potential; - int index; -}; +void genDistusage() { + std::cout << ">> Remove or Add First moment This executable has to be used to compute direct interaction either for periodic or non periodic system.\n"; + std::cout << " Options " << std::endl + << " -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 + << " -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); +} // Simply create particles and try the kernels int main(int argc, char ** argv){ // 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 ; + + std::cout << "Unsur4pie0: " << Unsur4pie0 << " 8.98755179e+9 Diff " << Unsur4pie0-8.98755179e+9 <> Remove or Add First moment This executable has to be used to compute direct interaction either for periodic or non periodic system.\n"; - std::cout << " Options " << std::endl - << " -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 - << " -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); + genDistusage() ; } if( FParameters::existParameter(argc, argv, "-Ewald2FMM") &&FParameters::existParameter(argc, argv, "-FMM2Ewald") ){ std::cerr << " Only -Ewald2FMM or FDMM2Ewald have to be set " <* particlesIn ; //= new FmaRWParticle<8,8>[numberofParticles]; - if(! loader->isOpen()){ - std::cout << "Loader Error, " << filenameIn << " is missing\n"; - exit(EXIT_FAILURE); + if (useDLPOLY) { + // ----------------------------------------------------- + // LOADER EWALD PARTICLES FROM DLPOLY + // ----------------------------------------------------- + if(! loaderDlpoly->isOpen()){ + std::cout << "Loader Error, " << filenameIn << " is missing\n"; + exit(EXIT_FAILURE); + } + boxsize[0] = boxsize[1]= boxsize[2]=loaderDlpoly->getBoxWidth() ; + centre = loaderDlpoly->getCenterOfBox(); + numberofParticles = loaderDlpoly->getNumberOfParticles() ; + particlesIn = new FmaRWParticle<8,8>[numberofParticles]; + memset(particlesIn, 0, sizeof( FmaRWParticle<8,8>) * numberofParticles) ; + + for(int idxPart = 0 ; idxPart < numberofParticles; ++idxPart){ + // + FPoint position; + FReal forces[3]; + FReal physicalValue; + FReal potential=0.0; + int index ; + loaderDlpoly->fillParticle(&position, forces,&physicalValue,&index); + particlesIn[index].setPosition(position) ; + particlesIn[index].setPhysicalValue(physicalValue) ; + particlesIn[index].setPotential(potential) ; + particlesIn[index].setForces(forces[0],forces[1],forces[2]) ; + } } - // --------------------------------------------------------------------------------- - // Read particles - // --------------------------------------------------------------------------------- - double boxsize[3] ; - boxsize[0] = boxsize[1]= boxsize[2]=loader->getBoxWidth() ; + else if (useFMA) { + // ----------------------------------------------------- + // LOADER EWALD PARTICLES FROM STAMP + // ----------------------------------------------------- + if(! loaderFMA->isOpen()){ + std::cout << "Loader Error, " << filenameIn << " is missing\n"; + exit(EXIT_FAILURE); + } + boxsize[0] = boxsize[1]= boxsize[2] = loaderFMA->getBoxWidth()/ang; + + // scaleEnergy /= boxsize[0] ; + // scaleForce = -scaleEnergy/boxsize[0] ; + // boxsize[0] = boxsize[1]= boxsize[2] = 1.0 ; + - FSize numberofParticles = loader->getNumberOfParticles() ; - std::cout << "Creating & Inserting " << numberofParticles<< " particles ..." << std::endl; - std::cout << "\tWidth : " << loader->getBoxWidth() << " \t center x : " << loader->getCenterOfBox().getX() - << " y : " << loader->getCenterOfBox().getY() << " z : " << loader->getCenterOfBox().getZ() << std::endl; + numberofParticles = loaderFMA->getNumberOfParticles() ; + 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 - MDParticle * const particles = new MDParticle[numberofParticles]; - memset(particles, 0, sizeof(MDParticle) * numberofParticles) ; // double totalCharge = 0.0; // for(int idxPart = 0 ; idxPart < numberofParticles; ++idxPart){ // - loader->fillParticle(&particles[idxPart].position, particles[idxPart].forces, - &particles[idxPart].physicalValue,&particles[idxPart].index); + totalCharge += particlesIn[idxPart].getPhysicalValue() ; // - totalCharge += particles[idxPart].physicalValue ; - // - electricMoment.incX(particles[idxPart].physicalValue*particles[idxPart].position.getX() ); - electricMoment.incY(particles[idxPart].physicalValue*particles[idxPart].position.getY() ); - electricMoment.incZ(particles[idxPart].physicalValue*particles[idxPart].position.getZ() ); + 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() ); } counter.tac(); @@ -159,20 +220,13 @@ int main(int argc, char ** argv){ std::cout << "Done " << "(@Reading Ewald file = " << counter.elapsed() << " s)." << std::endl; // - // - // --------------------------------------------------------------------------------- - // READ DIRECT COMPUTATION - // ---------------------------------------------------------------- - // Save direct computation in binary format - // write binary output file - // // ---------------------------------------------------------------------------- // Remove or Add Electric Moment //----------------------------------------------------------------------------- // remove polarization component // double volume = boxsize[0] *boxsize[1]*boxsize[2] ; - double coeffCorrection = 2.0*FMath::FPi/volume/3.0 ; + double coeffCorrection = 4.0*FMath::FPi/volume/3.0 ; 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 ; @@ -193,49 +247,42 @@ int main(int argc, char ** argv){ // double tmp, newEnergy =0.0; for(int idx = 0 ; idx < numberofParticles ; ++idx){ - // std::cout << " Pos " << particles[idx].position.getX()<< " "<< particles[idx].position.getY()<< " "<< particles[idx].position.getZ()<< std::endl - // << " F ori " << particles[idx].forces[0]<< " "<< particles[idx].forces[1]<< " "<< particles[idx].forces[2]<< std::endl; - - tmp = particles[idx].position.getX()*electricMoment.getX() + particles[idx].position.getY()*electricMoment.getY() - + particles[idx].position.getZ()*electricMoment.getZ() ; + tmp = particlesIn[idx].getPosition().getX()*electricMoment.getX() + particlesIn[idx].getPosition().getY()*electricMoment.getY() + + particlesIn[idx].getPosition().getZ()*electricMoment.getZ() ; // - particles[idx].potential += 2.0*tmp*coeffCorrection; + double P = particlesIn[idx].getPotential(); + particlesIn[idx].setPotential( P + 2.0*tmp*coeffCorrection); // - particles[idx].forces[0] = preScaleForce*particles[idx].forces[0] + 2.0*particles[idx].physicalValue*coeffCorrection*electricMoment.getX() ; - particles[idx].forces[1] = preScaleForce*particles[idx].forces[1] + 2.0*particles[idx].physicalValue*coeffCorrection*electricMoment.getY() ; - particles[idx].forces[2] = preScaleForce*particles[idx].forces[2] + 2.0*particles[idx].physicalValue*coeffCorrection*electricMoment.getZ() ; + 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() ; // - particles[idx].forces[0] *= postScaleForce; - particles[idx].forces[1] *= postScaleForce; - particles[idx].forces[2] *= postScaleForce; - // std::cout << " F new " << particles[idx].forces[0]<< " "<< particles[idx].forces[1]<< " "<< particles[idx].forces[2]<< std::endl<< std::endl; - + forces[0] *= postScaleForce; + forces[1] *= postScaleForce; + forces[2] *= postScaleForce; + particlesIn[idx].setForces(forces[0],forces[1],forces[2]) ; // } + std::cout.precision(10) ; + std::cout <getEnergy() <getEnergy() + coeffCorrection*electricMoment.norm2() ; - newEnergy *= postScaleEnergy ; + if (useDLPOLY){ + std::cout << " " << coeffCorrection*electricMoment.norm2() <getEnergy() <getEnergy() + coeffCorrection*electricMoment.norm2() ; + newEnergy *= postScaleEnergy ; - // - std::cout.precision(10) ; - std::cout << "Energy EWALD = "<< loader->getEnergy() <getCenterOfBox(), loader->getBoxWidth() , numberofParticles,*particlesOut) ; - writer.writeArrayOfParticles(particlesOut, numberofParticles); + writer.writeHeader(centre, boxsize[0] , numberofParticles,*particlesIn) ; + writer.writeArrayOfParticles(particlesIn, numberofParticles); // - delete[] particles; - + delete[] particlesIn; return 0; } diff --git a/Utils/python/manipFMA.py b/Utils/python/manipFMA.py new file mode 100644 index 0000000000000000000000000000000000000000..7cbf244307ae7d2f85b7dc11f3305de74b368562 --- /dev/null +++ b/Utils/python/manipFMA.py @@ -0,0 +1,87 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon Sep 8 17:25:37 2014 + +@author: coulaud +""" +#Read file comming from Stamp + +import numpy +#import sys + +#Filename='/Users/coulaud/Dev/src/ScalFMM/scalfmmT/Data/cea-200-per.xyzqf' + +# +# Read Header +# +#global nbRecordPerline,dataType,NbParticles,BoxWidth,Centre +dataType = 8 +nbRecordPerline = 4 +NbParticles = 0 +BoxWidth = 1.0 +Centre= (0.0,0.0,0.0) +# +def readFMAHeader(Filename): + """ + This function read the header of the FMA file Filename + """ + global nbRecordPerline,dataType,NbParticles,BoxWidth,Centre + f = open(Filename, 'r') + header1 = f.readline() + header2 = f.readline() + s1,s2=(item.strip() for item in header1.split()) + dataType = int(s1) + nbRecordPerline = int(s2) + s1,s2,s3,s4,s5 = (item.strip() for item in header2.split()) + NbParticles=int(s1) + BoxWidth=float(s2)*2 + Centre = [float(s3),float(s4),float(s5)] + f.close() + return dataType,nbRecordPerline,NbParticles,BoxWidth,Centre +# +def readFMAfile(fileName): + """ + Read data from an ascci file il the FMA format + Return data (x,y,z,q) or (x,y,z,q,p,fx,fy,fz) + """ + global nbRecordPerline,dataType,NbParticles,BoxWidth,Centre + readFMAHeader(fileName) + print NbParticles + # Read data in a non compact mode + # + data = numpy.loadtxt(fileName,skiprows=2) + +# if nbRecordPerline == 4 or nbRecordPerline == 8 : + data = numpy.loadtxt(fileName,skiprows=2) +# else: + #print " wrong data nbRecordPerline (",nbRecordPerline,"). Should be 4 or 8" + # os._exit() + return data,BoxWidth,Centre + +# +def writeFMAfilefmt(fileName,data,boxWidth, centre,myformat): + """ + write data from an ascci file il the FMA format + Return data (x,y,z,q) or (x,y,z,q,p,fx,fy,fz) + """ + # f = open(fileName, "w") + s = ' 8 '+ str(numpy.shape(data)[1]) + '\n' + # f.write(s) + s= s + str(numpy.shape(data)[0])+ ' '+ str(boxWidth/2)+ ' '+ str(Centre[0])+ ' '+ str(Centre[1])+ ' '+ str(Centre[2]) + # f.write(s) + # f.close() + numpy.savetxt(fileName,data,delimiter=' ',header=s,comments=' ',fmt=myformat) + +def writeFMAfile(fileName,data,boxWidth, centre): + """ + write data from an ascci file il the FMA format + Return data (x,y,z,q) or (x,y,z,q,p,fx,fy,fz) + """ + # f = open(fileName, "w") + s = ' 8 '+ str(numpy.shape(data)[1]) + '\n' + # f.write(s) + s= s + str(numpy.shape(data)[0])+ ' '+ str(boxWidth/2)+ ' '+ str(Centre[0])+ ' '+ str(Centre[1])+ ' '+ str(Centre[2]) + # f.write(s) + # f.close() + numpy.savetxt(fileName,data,delimiter=' ',header=s,comments=' ') + diff --git a/Utils/scripts/manipFMA.py b/Utils/scripts/manipFMA.py deleted file mode 100644 index e46eda09e2e6c1f7de6ec0e170d949965d1cb0ab..0000000000000000000000000000000000000000 --- a/Utils/scripts/manipFMA.py +++ /dev/null @@ -1,58 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Sun May 11 16:32:08 2014 - -@author: coulaud -""" - -import numpy as np -from random import * -import matplotlib.pyplot as plt -# -# Manipulation de fichier FMA -# -#Npart =0 - -#def readFMA(inputFile,N) -def readFMA(N, BoxSize) : -# -#part BoxSize,BoxCentre): -# - print __doc__ - inputFile="/Users/coulaud/Dev/src/ScalFMM/scalfmmT/Data/prolate50.fma" - fichier = open(inputFile, "r") - A= fichier.readline().split() - N = int(A[0]) - BoxSize = float(A[1]) - fichier.close() - return N - -def writeFMA(N, xyzq, outFile) : -# -#part BoxSize,BoxCentre): -# - print __doc__ -# outFile="/Users/coulaud/Dev/src/ScalFMM/scalfmmT/Data/prolate50.fma" - fichier = open(outFile, "r") - # A= fichier.readline().split() - # N = int(A[0]) - # BoxSize = float(A[1]) - fichier.close() - - -def build1DDistribution(a, alpha,N): - xyzq = np.zeros((N, 4)) - x = np.zeros(N) - y = np.zeros(N) - z = np.random.random_sample(N) - - sa = sin(alpha) - ca = cos(alpha) - print a,alpha,ca,sa - for i in range(0, N): - x[i] = a*sa*z[i]*cos(z[i]) - y[i] = a*sa*z[i]*sin(z[i]) - z[i] = a*ca*z[i] - print i, x[i],y[i], z[i] - return x,y,z - diff --git a/email-sorted b/email-sorted index fe0a2716e8fa6b1dd4c13537486b336ea4b1cd6b..4f0b77c281443248a6ea2a01333f72e7a1acab39 100644 --- a/email-sorted +++ b/email-sorted @@ -10,6 +10,7 @@ ccowingzitron@ucsd.edu cheekubadshah@sify.com dhairyamal@gmail.com fabien.casenave@gmail.com +guillaume.fuhr@univ-amu.fr jean-christophe.toussaint@grenoble.cnrs.fr j8asic@gmail.com joseluiscasalssainz@gmail.com @@ -29,6 +30,7 @@ ronojoy.adhikari@gmail.com rxrjj@126.com samuel.thibault@labri.fr shujin.cao@163.com +t.ben.thompson@gmail.com t.betcke@ucl.ac.uk vitoreafeliciano@msn.com weygand@kit.edu