removeMoment.cpp 13.4 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
// ===================================================================================
// ===================================================================================
// 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.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.  
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info". 
// "http://www.gnu.org/licenses".
// ===================================================================================
//

//
//
//
// Tests/Release/removeMoment  -Ewald2FMM -fin ../Data/forceNacl_2000_dlpolyPer.bin -dlpoly -fout ../Data/forceNacl_2000.fma


#include <iostream>
#include <iomanip>

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <algorithm>

#include  "ScalFmmConfig.h"
#include "Utils/FTic.hpp"
#include "Utils/FParameters.hpp"


#include "Files/FDlpolyLoader.hpp"
#include "Files/FFmaGenericLoader.hpp"


/** DLpoly particle is used in the gadget program
 * here we try to make the same simulation
 */

COULAUD Olivier's avatar
COULAUD Olivier committed
46 47 48 49 50 51 52 53 54 55
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);
}
56 57 58 59 60
// Simply create particles and try the kernels
int main(int argc, char ** argv){
	//
	FReal scaleEnergy, scaleForce  ;
	const FReal r4pie0 = FReal(138935.4835);
COULAUD Olivier's avatar
COULAUD Olivier committed
61 62 63 64 65 66
	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 <<std::endl;
67

COULAUD Olivier's avatar
COULAUD Olivier committed
68 69 70
	if( argc == 1 ){
		genDistusage() ;
	}
71 72
	///////////////////////What we do/////////////////////////////
	if( FParameters::existParameter(argc, argv, "-help") || FParameters::existParameter(argc, argv, "-h") ){
COULAUD Olivier's avatar
COULAUD Olivier committed
73
		genDistusage() ;
74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
	}
	if( FParameters::existParameter(argc, argv, "-Ewald2FMM") &&FParameters::existParameter(argc, argv, "-FMM2Ewald") ){
		std::cerr << " Only -Ewald2FMM or FDMM2Ewald have to be set " <<std::endl;
		exit(EXIT_FAILURE);
	}
	///
	// ---------------------------------------------------
	//        DL_POLY CONSTANT
	//  ---------------------------------------------------
	scaleEnergy =  1.0;
	scaleForce   = 1.0 ;
	//	bool scale = true ;
	if(FParameters::existParameter(argc, argv, "-dlpoly")){
		//		scale = false ;
		scaleEnergy =  r4pie0 / 418.4 ;   // kcal mol^{-1}
		scaleForce  = -r4pie0 ;           // 10 J mol^{-1} A^{-1}
	}
COULAUD Olivier's avatar
COULAUD Olivier committed
91 92 93 94
	else 	if(FParameters::existParameter(argc, argv, "-stamp")){
		scaleEnergy =  q0*q0*Unsur4pie0/ang ;
		scaleForce  = -scaleEnergy/ang;
	}
95 96 97 98 99 100 101 102 103 104 105 106
	std::cout <<"Scaling factor " <<std::endl
			<<      "Energy factor " << scaleEnergy<<std::endl
			<<      "  Force factor " << scaleForce<<std::endl ;
	//////////////////////////////////////////////////////////////
	//
	const std::string    filenameIn (FParameters::getStr(argc,argv,"-fin", "../Data/forceNacl_128_dlpolyPer.bin"));
	const std::string    filenameOut (FParameters::getStr(argc,argv,"-fout", "result.fma"));
	//  file for -saveError option

	FTic counter;

	std::cout << "Opening : " << filenameIn << "\n";
COULAUD Olivier's avatar
COULAUD Olivier committed
107 108 109
	FDlpolyLoader          *loaderDlpoly = nullptr ;
	FFmaGenericLoader  *loaderFMA = nullptr ;
	bool useDLPOLY=false,useFMA=false;
110
	std::string ext(".bin");
COULAUD Olivier's avatar
COULAUD Olivier committed
111 112
	std::string ext1(".fma"),ext2(".bfma");

113 114
	// open particle file
	if(filenameIn.find(ext) != std::string::npos) {
COULAUD Olivier's avatar
COULAUD Olivier committed
115 116
		loaderDlpoly  = new FDlpolyBinLoader(filenameIn.c_str());
		useDLPOLY = true ;
117 118
	}
	else if(filenameIn.find(".txt")!=std::string::npos ) {
COULAUD Olivier's avatar
COULAUD Olivier committed
119 120
		loaderDlpoly  = new FDlpolyAsciiLoader(filenameIn.c_str());
		useDLPOLY = true ;
121
	}
COULAUD Olivier's avatar
COULAUD Olivier committed
122 123 124 125 126 127 128 129
	else if(filenameIn.find(ext1)!=std::string::npos ) {
		loaderFMA  = new FFmaGenericLoader(filenameIn.c_str());
		useFMA = true ;
	}  else if(filenameIn.find(ext2)!=std::string::npos ) {
		loaderFMA  = new FFmaGenericLoader(filenameIn.c_str());
		useFMA    = true ;
	}  else
	{
130 131 132
		std::cout << "Input file not allowed only .bin or .txt extensions" <<std::endl;
		std::exit ( EXIT_FAILURE) ;
	}
COULAUD Olivier's avatar
COULAUD Olivier committed
133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
	double boxsize[3] ;
	FPoint centre ;
	FSize numberofParticles ;
	FmaRWParticle<8,8>*  particlesIn  ; //= new FmaRWParticle<8,8>[numberofParticles];

	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]) ;
		}
	}
	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 ;
		centre = loaderFMA->getCenterOfBox();
		centre *= boxsize[0];
		numberofParticles                    = loaderFMA->getNumberOfParticles() ;
		particlesIn  = new FmaRWParticle<8,8>[numberofParticles];
		memset(particlesIn, 0, sizeof( FmaRWParticle<8,8>) * numberofParticles) ;
181

COULAUD Olivier's avatar
COULAUD Olivier committed
182 183 184 185 186 187 188 189
		loaderFMA->fillParticle(particlesIn,numberofParticles);
		for(int idxPart = 0 ; idxPart < numberofParticles; ++idxPart){
			//
			FPoint PP = particlesIn[idxPart].getPosition() ;
			PP *=  boxsize[0];
			particlesIn[idxPart].setPosition(PP) ;
			//
		}
190
	}
COULAUD Olivier's avatar
COULAUD Olivier committed
191

192 193 194 195 196
	// ---------------------------------------------------------------------------------
	//  Read particles
	// ---------------------------------------------------------------------------------


COULAUD Olivier's avatar
COULAUD Olivier committed
197 198 199 200
	//	std::cout << "Creating & Inserting "q << 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;
	//
201 202 203 204 205 206 207 208 209

	counter.tic();
	FPoint electricMoment(0.0,0.0,0.0) ;
	// const --> then shared
	//
	double totalCharge = 0.0;
	//
	for(int idxPart = 0 ; idxPart < numberofParticles; ++idxPart){
		//
COULAUD Olivier's avatar
COULAUD Olivier committed
210
		totalCharge += particlesIn[idxPart].getPhysicalValue() ;
211
		//
COULAUD Olivier's avatar
COULAUD Olivier committed
212 213 214
		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() );
215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263
	}
	counter.tac();

	std::cout << std::endl;
	std::cout << "Total Charge               = "<< totalCharge <<std::endl;
	std::cout << "Electric Moment          = "<< electricMoment <<std::endl;
	std::cout << "Electric Moment norm = "<< electricMoment.norm2()  <<std::endl;
	std::cout << "----------------------------------------------------"<<std::endl;
	std::cout << std::endl;

	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 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 ;
	}
	else if( FParameters::existParameter(argc, argv, "-FMM2Ewald") ){
		preScaleEnergy  = 1.0 ;  postScaleEnergy = scaleEnergy ; preScaleForce = 1.0/scaleForce;  postScaleForce = scaleForce ;
		coeffCorrection  = -coeffCorrection;
	}
	else {
		std::cout << " -Ewald2FMM ou -FMM2Ewald should be set"<<std::endl;
		exit(EXIT_FAILURE);
	}
	std::cout << "coeffCorrection: "<<coeffCorrection<<std::endl;
	std::cout << "preScaleEnergy:  "<<preScaleEnergy<<std::endl;
	std::cout << "postScaleEnergy: "<<postScaleEnergy<<std::endl;
	std::cout << "preScaleForce:     "<<preScaleForce<<std::endl;
	std::cout << "postScaleForce:   "<<postScaleForce<<std::endl;
	//
	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;

COULAUD Olivier's avatar
COULAUD Olivier committed
264 265
		tmp = particlesIn[idx].getPosition().getX()*electricMoment.getX()  + particlesIn[idx].getPosition().getY()*electricMoment.getY()
																														+ particlesIn[idx].getPosition().getZ()*electricMoment.getZ()  ;
266
		//
COULAUD Olivier's avatar
COULAUD Olivier committed
267 268
		double P = particlesIn[idx].getPotential();
		particlesIn[idx].setPotential(  P + 2.0*tmp*coeffCorrection);
269
		//
COULAUD Olivier's avatar
COULAUD Olivier committed
270 271 272 273 274
		double * forces ;
		forces  = particlesIn[idx].getForces() ;
		forces[0] = preScaleForce*forces[0] + 2.0*particlesIn[idx].getPhysicalValue()*coeffCorrection*electricMoment.getX() ;
		forces[1] = preScaleForce*forces[1] + 2.0*particlesIn[idx].getPhysicalValue()*coeffCorrection*electricMoment.getY() ;
		forces[2] = preScaleForce*forces[2] + 2.0*particlesIn[idx].getPhysicalValue()*coeffCorrection*electricMoment.getZ() ;
275
		//
COULAUD Olivier's avatar
COULAUD Olivier committed
276 277 278 279 280
		forces[0] *= postScaleForce;
		forces[1] *= postScaleForce;
		forces[2] *= postScaleForce;
		particlesIn[idx].setForces(forces[0],forces[1],forces[2]) ;

281 282 283 284
		//	      std::cout <<  " F  new " << particles[idx].forces[0]<< "  "<< particles[idx].forces[1]<< "  "<< particles[idx].forces[2]<< std::endl<< std::endl;

		//
	}
COULAUD Olivier's avatar
COULAUD Olivier committed
285 286
	std::cout.precision(10) ;

287 288 289
	std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
	std::cout << std::scientific;
	std::cout.precision(10) ;
COULAUD Olivier's avatar
COULAUD Olivier committed
290 291 292 293 294
	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 ;
295

COULAUD Olivier's avatar
COULAUD Olivier committed
296 297 298 299

		//
		std::cout << "Energy EWALD  = "<< loaderDlpoly->getEnergy() <<std::endl ;
	}
300 301 302 303
	std::cout << "Energy New      = "<<newEnergy<<std::endl ;
	//
	// Save data in FMAFormat
	//FSize nbParticles = loader.getNumberOfParticles() ;
COULAUD Olivier's avatar
COULAUD Olivier committed
304 305 306 307 308 309 310 311 312 313
	//	FmaRWParticle<8,8>* const particlesOut = new FmaRWParticle<8,8>[numberofParticles];
	//	// remove index
	//	memset(particlesOut, 0, sizeof(FmaRWParticle<8,8>) * numberofParticles) ;
	//	for(int idx = 0 ; idx < numberofParticles ; ++idx){
	//		particlesOut[idx].setPosition(particles[idx].position) ;
	//		particlesOut[idx].setPhysicalValue(particles[idx].physicalValue) ;
	//		particlesOut[idx].setPotential (particles[idx].potential) ;
	//		particlesOut[idx].setForces(particles[idx].forces[0],particles[idx].forces[0],particles[idx].forces[2]) ;
	//	}
	//	//
314 315 316 317 318 319 320
	// ----------------------------------------------------------------
	//  Save  computation in binary format
	//
	//

	std::cout << "Generate " << filenameOut <<"  for output file" << std::endl;
	//
COULAUD Olivier's avatar
COULAUD Olivier committed
321 322 323
	//	std::cout << " numberofParticles: " << nbParticles <<"  " << sizeof(numberofParticles) <<std::endl;
	//	std::cout << " denergy: " << newEnergy <<"  " << sizeof(newEnergy) <<std::endl;
	//	std::cout << " Box size: " << loader.getBoxWidth() << "  " << sizeof(loader.getBoxWidth())<<std::endl;
324 325
	//
	FFmaGenericWriter writer(filenameOut) ;
COULAUD Olivier's avatar
COULAUD Olivier committed
326 327
	writer.writeHeader(centre, boxsize[0] , numberofParticles,*particlesIn) ;
	writer.writeArrayOfParticles(particlesIn, numberofParticles);
328 329

	//
COULAUD Olivier's avatar
COULAUD Olivier committed
330
	delete[] particlesIn;
331 332 333 334 335
	return 0;
}