Skip to content

Commit

Permalink
Refactored susceptibility calculation into statistics module
Browse files Browse the repository at this point in the history
  • Loading branch information
richard-evans committed Dec 1, 2014
1 parent 45dd902 commit 9a866b0
Show file tree
Hide file tree
Showing 9 changed files with 231 additions and 85 deletions.
39 changes: 34 additions & 5 deletions hdr/stats.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -52,10 +52,6 @@ namespace stats

extern double torque_data_counter;

extern double mean_susceptibility[4];
extern double mean_susceptibility_squared[4];
extern bool calculate_susceptibility;

extern bool calculate_energy;

/// Statistics energy types
Expand All @@ -82,15 +78,21 @@ namespace stats
extern bool calculate_material_magnetization;
extern bool calculate_height_magnetization;
extern bool calculate_material_height_magnetization;
extern bool calculate_system_susceptibility;

class susceptibility_statistic_t;

//----------------------------------
// Magnetization Class definition
//----------------------------------
class magnetization_statistic_t{

friend class susceptibility_statistic_t;

public:
//magnetization_statistic_t (const int in_mask_size, std::vector<int> in_mask);
magnetization_statistic_t ();
bool is_initialized();
void set_mask(const int mask_size, std::vector<int> inmask, const std::vector<double>& mm);
void calculate_magnetization(const std::vector<double>& sx, const std::vector<double>& sy, const std::vector<double>& sz, const std::vector<double>& mm);
void reset_magnetization_averages();
Expand All @@ -103,7 +105,7 @@ namespace stats
std::string output_normalized_magnetization_dot_product(const std::vector<double>& vec);

private:
bool is_initialized;
bool initialized;
int num_atoms;
int mask_size;
double mean_counter;
Expand All @@ -115,12 +117,39 @@ namespace stats

};

//----------------------------------
// Susceptibility Class definition
//----------------------------------
class susceptibility_statistic_t{

public:
susceptibility_statistic_t ();
void initialize(magnetization_statistic_t& mag_stat);
void calculate(const std::vector<double>& magnetization);
void reset_averages();
std::string output_mean_susceptibility(const double temperature);
//std::string output_mean_absolute_susceptibility();

private:
bool initialized;
int num_elements;
double mean_counter;
std::vector<double> mean_susceptibility;
std::vector<double> mean_susceptibility_squared;
std::vector<double> mean_absolute_susceptibility;
std::vector<double> mean_absolute_susceptibility_squared;
std::vector<double> saturation;

};

// Statistics classes
extern magnetization_statistic_t system_magnetization;
extern magnetization_statistic_t material_magnetization;
extern magnetization_statistic_t height_magnetization;
extern magnetization_statistic_t material_height_magnetization;

extern susceptibility_statistic_t system_susceptibility;
//extern susceptibility_statistic_t material_susceptibility;

}

Expand Down
1 change: 1 addition & 0 deletions makefile
Original file line number Diff line number Diff line change
Expand Up @@ -120,6 +120,7 @@ obj/statistics/data.o \
obj/statistics/initialize.o \
obj/statistics/magnetization.o \
obj/statistics/statistics.o \
obj/statistics/susceptibility.o \
obj/utility/checkpoint.o \
obj/utility/errors.o \
obj/utility/statistics.o \
Expand Down
11 changes: 7 additions & 4 deletions src/statistics/data.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -14,16 +14,19 @@

namespace stats{

bool calculate_system_magnetization=false;
bool calculate_material_magnetization=false;
bool calculate_height_magnetization=false;
bool calculate_material_height_magnetization=false;
bool calculate_system_magnetization = false;
bool calculate_material_magnetization = false;
bool calculate_height_magnetization = false;
bool calculate_material_height_magnetization = false;
bool calculate_system_susceptibility = false;

magnetization_statistic_t system_magnetization;
magnetization_statistic_t material_magnetization;
magnetization_statistic_t height_magnetization;
magnetization_statistic_t material_height_magnetization;

susceptibility_statistic_t system_susceptibility;

//-----------------------------------------------------------------------------
// Shared variables used for statistics calculation
//-----------------------------------------------------------------------------
Expand Down
3 changes: 3 additions & 0 deletions src/statistics/initialize.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,9 @@ namespace stats{
stats::material_height_magnetization.set_mask(num_materials*(max_height+1),mask,magnetic_moment_array);
}

// system susceptibility
stats::system_susceptibility.initialize(stats::system_magnetization);

return;
}
} // end of namespace stats
12 changes: 11 additions & 1 deletion src/statistics/magnetization.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,14 @@ namespace stats{
//------------------------------------------------------------------------------------------------------
// Constructor to initialize data structures
//------------------------------------------------------------------------------------------------------
magnetization_statistic_t::magnetization_statistic_t (): is_initialized(false){}
magnetization_statistic_t::magnetization_statistic_t (): initialized(false){}

//------------------------------------------------------------------------------------------------------
// Function to determine if class is properly initialized
//------------------------------------------------------------------------------------------------------
bool magnetization_statistic_t::is_initialized(){
return initialized;
}

//------------------------------------------------------------------------------------------------------
// Function to initialize mask
Expand Down Expand Up @@ -85,6 +92,9 @@ void magnetization_statistic_t::set_mask(const int in_mask_size, std::vector<int
}
}

// Set flag indicating correct initialization
initialized=true;

return;

}
Expand Down
6 changes: 6 additions & 0 deletions src/statistics/statistics.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,9 @@ namespace stats{
if(stats::calculate_height_magnetization) stats::height_magnetization.calculate_magnetization(sx,sy,sz,mm);
if(stats::calculate_material_height_magnetization) stats::material_height_magnetization.calculate_magnetization(sx,sy,sz,mm);

// update susceptibility statistics
if(stats::calculate_system_susceptibility) stats::system_susceptibility.calculate(stats::system_magnetization.get_magnetization());

return;

}
Expand All @@ -48,6 +51,9 @@ namespace stats{
if(stats::calculate_height_magnetization) stats::height_magnetization.reset_magnetization_averages();
if(stats::calculate_material_height_magnetization) stats::material_height_magnetization.reset_magnetization_averages();

// reset susceptibility statistics
if(stats::calculate_system_susceptibility) stats::system_susceptibility.reset_averages();

return;

}
Expand Down
165 changes: 165 additions & 0 deletions src/statistics/susceptibility.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,165 @@
//-----------------------------------------------------------------------------
//
// This source file is part of the VAMPIRE open source package under the
// GNU GPL (version 2) licence (see licence file for details).
//
// (c) R F L Evans 2014. All rights reserved.
//
//-----------------------------------------------------------------------------

// C++ standard library headers
#include <algorithm>
#include <cmath>
#include <iostream>
#include <sstream>

// Vampire headers
#include "errors.hpp"
#include "stats.hpp"
#include "vmpi.hpp"
#include "vio.hpp"

namespace stats{

//------------------------------------------------------------------------------------------------------
// Constructor
//------------------------------------------------------------------------------------------------------
susceptibility_statistic_t::susceptibility_statistic_t (){}

//------------------------------------------------------------------------------------------------------
// Function to initialize data structures
//------------------------------------------------------------------------------------------------------
void susceptibility_statistic_t::initialize(stats::magnetization_statistic_t& mag_stat) {

// Check that magnetization statistic is properly initialized
if(!mag_stat.is_initialized()){
terminaltextcolor(RED);
std::cerr << "Programmer Error - Uninitialized magnetization statistic passed to susceptibility statistic - please initialize first." << std::endl;
terminaltextcolor(WHITE);
zlog << zTs() << "Programmer Error - Uninitialized magnetization statistic passed to susceptibility statistic - please initialize first." << std::endl;
err::vexit();
}

// Determine number of magnetization statistics*4
std::vector<double> temp = mag_stat.get_magnetization();
num_elements = temp.size()/4;

// Now set number of susceptibility values to match
mean_susceptibility.resize(4*num_elements,0.0);
mean_susceptibility_squared.resize(4*num_elements,0.0);
mean_absolute_susceptibility.resize(4*num_elements,0.0);
mean_absolute_susceptibility_squared.resize(4*num_elements,0.0);

// copy saturation data
saturation = mag_stat.saturation;

// initialize mean counter
mean_counter = 0.0;

// Set flag indicating correct initialization
initialized=true;

}

//------------------------------------------------------------------------------------------------------
// Function to calculate susceptibility of the magnetisation and retain the mean value
//
// chi_l = sum_i mu_i
// ---------- ( <m_l^2> - <m_l>^2 )
// k_B T
//
// m_l = sum_i mu_i S_i
// --------------
// sum_i mu_i
//
//-------------------------------------------------------------------------------------------------------
void susceptibility_statistic_t::calculate(const std::vector<double>& magnetization){

// loop over all elements
for(int id=0; id< num_elements; ++id){

// copy reduced magnetisation
const double mx = magnetization[4*id + 0];
const double my = magnetization[4*id + 1];
const double mz = magnetization[4*id + 2];
const double mm = magnetization[4*id + 3];

mean_susceptibility[4*id + 0]+=mx*mm;
mean_susceptibility[4*id + 1]+=my*mm;
mean_susceptibility[4*id + 2]+=mz*mm;
mean_susceptibility[4*id + 3]+=mm;

mean_susceptibility_squared[4*id + 0]+=mx*mx*mm*mm;
mean_susceptibility_squared[4*id + 1]+=my*my*mm*mm;
mean_susceptibility_squared[4*id + 2]+=mz*mz*mm*mm;
mean_susceptibility_squared[4*id + 3]+=mm*mm;

mean_absolute_susceptibility[4*id + 0]+=fabs(mx*mm);
mean_absolute_susceptibility[4*id + 1]+=fabs(my*mm);
mean_absolute_susceptibility[4*id + 2]+=fabs(mz*mm);
mean_absolute_susceptibility[4*id + 3]+=mm;

mean_absolute_susceptibility_squared[4*id + 0]+=fabs(mx*mx*mm*mm);
mean_absolute_susceptibility_squared[4*id + 1]+=fabs(my*my*mm*mm);
mean_absolute_susceptibility_squared[4*id + 2]+=fabs(mz*mz*mm*mm);
mean_absolute_susceptibility_squared[4*id + 3]+=mm*mm;

}

mean_counter+=1.0;

return;

}

//------------------------------------------------------------------------------------------------------
// Function to reset statistical averages
//------------------------------------------------------------------------------------------------------
void susceptibility_statistic_t::reset_averages(){

// reinitialise mean magnetization to zero
std::fill(mean_susceptibility.begin(),mean_susceptibility.end(),0.0);
std::fill(mean_susceptibility_squared.begin(),mean_susceptibility_squared.end(),0.0);
std::fill(mean_absolute_susceptibility.begin(),mean_absolute_susceptibility.end(),0.0);
std::fill(mean_absolute_susceptibility_squared.begin(),mean_absolute_susceptibility_squared.end(),0.0);

// reset data counter
mean_counter = 0.0;

return;

}

//------------------------------------------------------------------------------------------------------
// Function to output mean susceptibility values as string
//------------------------------------------------------------------------------------------------------
std::string susceptibility_statistic_t::output_mean_susceptibility(const double temperature){

// result string stream
std::ostringstream result;

// determine inverse temperature mu_B/(kB T) (flushing to zero for very low temperatures)
const double itemp = temperature < 1.e-300 ? 0.0 : 9.274e-24/(1.3806503e-23*temperature);

// determine inverse mean counter and its square
const double imean_counter = 1.0/mean_counter;
const double imean_counter_sq = 1.0/(mean_counter*mean_counter);

// loop over all elements
for(int id=0; id< num_elements; ++id){

const double prefactor = itemp*saturation[id]; // in mu_B
const double sus_x = prefactor*(mean_susceptibility_squared[4*id + 0]*imean_counter-mean_susceptibility[4*id + 0]*mean_susceptibility[4*id + 0]*imean_counter_sq);
const double sus_y = prefactor*(mean_susceptibility_squared[4*id + 1]*imean_counter-mean_susceptibility[4*id + 1]*mean_susceptibility[4*id + 1]*imean_counter_sq);
const double sus_z = prefactor*(mean_susceptibility_squared[4*id + 2]*imean_counter-mean_susceptibility[4*id + 2]*mean_susceptibility[4*id + 2]*imean_counter_sq);
const double sus_m = prefactor*(mean_susceptibility_squared[4*id + 3]*imean_counter-mean_susceptibility[4*id + 3]*mean_susceptibility[4*id + 3]*imean_counter_sq);

result << sus_x << "\t" << sus_y << "\t" << sus_z << "\t" << sus_m << "\t";

}

return result.str();

}

} // end of namespace stats
Loading

0 comments on commit 9a866b0

Please sign in to comment.