/*
* Project: MoleCuilder
* Description: creates and alters molecular systems
* Copyright (C) 2010-2012 University of Bonn. All rights reserved.
*
*
* This file is part of MoleCuilder.
*
* MoleCuilder is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 2 of the License, or
* (at your option) any later version.
*
* MoleCuilder 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 License for more details.
*
* You should have received a copy of the GNU General Public License
* along with MoleCuilder. If not, see .
*/
/** \file datacreator.cpp
*
* Declarations of assisting functions in creating data and plot files.
*
*/
//============================ INCLUDES ===========================
// include config.h
#ifdef HAVE_CONFIG_H
#include
#endif
#include "CodePatterns/MemDebug.hpp"
#include
#include
#include
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "LinearAlgebra/defs.hpp"
#include "Helpers/defs.hpp"
#include "datacreator.hpp"
#include "Fragmentation/EnergyMatrix.hpp"
#include "Fragmentation/ForceMatrix.hpp"
#include "Fragmentation/HessianMatrix.hpp"
#include "Fragmentation/KeySetsContainer.hpp"
using namespace std;
//=========================== FUNCTIONS============================
/** Opens a file with \a *filename in \a *dir.
* \param output file handle on return
* \param *dir directory
* \param *filename name of file
* \return true if file has been opened
*/
bool OpenOutputFile(ofstream &output, const char *dir, const char *filename)
{
stringstream name;
name << dir << "/" << filename;
output.open(name.str().c_str(), ios::out);
if (output == NULL) {
LOG(0, "Unable to open " << name.str() << " for writing, is directory correct?");
return false;
}
return true;
};
/** Opens a file for appending with \a *filename in \a *dir.
* \param output file handle on return
* \param *dir directory
* \param *filename name of file
* \return true if file has been opened
*/
bool AppendOutputFile(ofstream &output, const char *dir, const char *filename)
{
stringstream name;
name << dir << "/" << filename;
output.open(name.str().c_str(), ios::app);
if (output == NULL) {
LOG(0, "Unable to open " << name.str() << " for writing, is directory correct?");
return false;
}
return true;
};
/** Plots an energy vs. order.
* \param &Fragments EnergyMatrix class containing matrix values
* \param KeySets KeySetContainer class holding bond KeySetContainer::Order
* \param *prefix prefix in filename (without ending)
* \param *msg message to be place in first line as a comment
* \return true if file was written successfully
*/
bool CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum)
{
stringstream filename;
ofstream output;
filename << prefix << ".dat";
if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
LOG(0, msg);
output << "# " << msg << ", created on " << datum;
output << "#Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
for (int BondOrder=0;BondOrder fabs(Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][1])) {
for (int k=Fragments.ColumnCounter[ Fragments.MatrixCounter ];k--;)
Fragments.Matrix[ Fragments.MatrixCounter ][j][k] = Fragments.Matrix[ KeySets.OrderSet[BondOrder][i] ][j][k];
}
}
}
};
/** Plot matrix vs. fragment.
*/
bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySets, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateFragmentOrder)(class MatrixContainer &, class KeySetsContainer &, int))
{
stringstream filename;
ofstream output;
filename << prefix << ".dat";
if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
LOG(0, msg);
output << "# " << msg << ", created on " << datum;
output << "#Order\tFrag.No.\t" << Fragment.Header[ Fragment.MatrixCounter ] << endl;
// max
for (int BondOrder=0;BondOrder MYEPSILON) && (tmp < stored)) { // current force is greater than stored
for (int m=NDIM;m--;)
Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l+m] = Force.Matrix[MatrixNumber][ k ][l+m];
stored = tmp;
}
}
}
};
/** Scans forces for the mean in magnitude.
* Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
* \param Force ForceMatrix class containing matrix values
* \param MatrixNumber the index for the ForceMatrix::matrix array
*/
void CreateMeanForce(class MatrixContainer &Force, int MatrixNumber)
{
int divisor = 0;
for (int l=Force.ColumnCounter[MatrixNumber];l--;)
Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
for (int l=5;l MYEPSILON) divisor++;
}
tmp /= (double)divisor;
Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = tmp;
}
};
/** Scans forces for the maximum in magnitude.
* Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
* \param Force ForceMatrix class containing matrix values
* \param MatrixNumber the index for the ForceMatrix::matrix array
*/
void CreateMaximumForce(class MatrixContainer &Force, int MatrixNumber)
{
for (int l=5;l stored) { // current force is greater than stored
for (int m=NDIM;m--;)
Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l+m] = Force.Matrix[MatrixNumber][ k ][l+m];
stored = tmp;
}
}
}
};
/** Leaves the Force.Matrix as it is.
* \param Force ForceMatrix class containing matrix values
* \param MatrixNumber the index for the ForceMatrix::matrix array
*/
void CreateSameForce(class MatrixContainer &Force, int MatrixNumber)
{
// does nothing
};
/** Adds vectorwise all forces.
* Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
* \param Force ForceMatrix class containing matrix values
* \param MatrixNumber the index for the ForceMatrix::matrix array
*/
void CreateVectorSumForce(class MatrixContainer &Force, int MatrixNumber)
{
for (int l=Force.ColumnCounter[MatrixNumber];l--;)
Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
for (int l=0;l