<commit_before>#include "mpi_common.h"
#include <iostream>
#include <iomanip>
#include <fstream>
#include "dynamical.h"
#include "error.h"
#include "fcs_phonon.h"
#include "gruneisen.h"
#include "pointers.h"
#include "kpoint.h"
#include "memory.h"
#include "system.h"
#include "parsephon.h"
#include "write_phonons.h"
#include "../alm_c++/mathfunctions.h"
#include "relaxation.h"

using namespace PHON_NS;

Gruneisen::Gruneisen(PHON *phon): Pointers(phon){};
Gruneisen::~Gruneisen(){};

void Gruneisen::setup()
{
	//	MPI_Bcast(&delta_a, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

	//	if (mympi->my_rank == 0) {
	std::cout << std::endl; 
	std::cout << "Change in cell parameters : " << delta_a << std::endl;
	//	}

	memory->allocate(gruneisen, kpoint->nk, dynamical->neval);

	prepare_delta_fc2();
	prepare_newfc2();

	if (relaxation->quartic_mode) {
		prepare_newfc3();
		write_newinfo_all();
	}
}

void Gruneisen::calc_gruneisen()
{
	std::cout << "Calculating Gruneisen parameters ..." << std::endl;

	unsigned int nk = kpoint->nk;
	unsigned int ns = dynamical->neval;

	unsigned int i;
	unsigned int ik, is;

	double *eval_orig;
	double **eval_plus;
	double **eval_minus;

	double xk_tmp[3];


	std::complex<double> **evec_tmp;

	memory->allocate(evec_tmp, 1, 1); // dummy allocation
	memory->allocate(eval_orig, ns);
	memory->allocate(eval_plus, nk, ns);
	memory->allocate(eval_minus, nk, ns);

	for (ik = 0; ik < nk; ++ik){

		for (i = 0; i < 3; ++i) xk_tmp[i] = kpoint->xk[ik][i];

		if (fcs_phonon->is_fc2_ext) {
			dynamical->eval_k(xk_tmp, dynamical->kvec_na[ik], fcs_phonon->fc2_ext, eval_orig, evec_tmp, false);
			dynamical->eval_k(xk_tmp, dynamical->kvec_na[ik], fc2_plus_ext, eval_plus[ik], evec_tmp, false);
			dynamical->eval_k(xk_tmp, dynamical->kvec_na[ik], fc2_minus_ext, eval_minus[ik], evec_tmp, false);
		} else {
			dynamical->eval_k(xk_tmp, dynamical->kvec_na[ik], fcs_phonon->fc2, eval_orig, evec_tmp, false);
			dynamical->eval_k(xk_tmp, dynamical->kvec_na[ik], fc2_plus, eval_plus[ik], evec_tmp, false);
			dynamical->eval_k(xk_tmp, dynamical->kvec_na[ik], fc2_minus, eval_minus[ik], evec_tmp, false);
		}


		for (is = 0; is < ns; ++is) {
			gruneisen[ik][is] = (eval_plus[ik][is] - eval_minus[ik][is]) / (2.0 * delta_a) / (-6.0 * eval_orig[is]);
		}
	}

	for (ik = 0; ik < nk; ++ik) {
		for (is = 0; is < ns; ++is) {
			eval_plus[ik][is] = dynamical->freq(eval_plus[ik][is]);
			eval_minus[ik][is] = dynamical->freq(eval_minus[ik][is]);
		}
	}

	if (kpoint->kpoint_mode == 1) {
		std::string file_band_plus, file_band_minus;

		file_band_plus = input->job_title + ".band_+";
		file_band_minus = input->job_title + ".band_-";

		std::ofstream ofs_plus, ofs_minus;

		ofs_plus.open(file_band_plus.c_str(), std::ios::out);
		if (!ofs_plus) error->exit("calc_gruneisen", "Could not create band_plus file.");
		ofs_minus.open(file_band_minus.c_str(), std::ios::out);
		if (!ofs_minus) error->exit("calc_gruneisen", "Could not create band_minus file.");

		ofs_plus << "# Phonon energy (cm^-1) of the system expanded by " << std::setw(10) << delta_a * 100 << " %." << std::endl;
		ofs_minus << "# Phonon energy (cm^-1) of the system compressed by " << std::setw(10) << delta_a * 100 << " %." << std::endl;

		for (ik = 0; ik < nk; ++ik){
			ofs_plus << std::setw(8) << std::fixed << kpoint->kaxis[ik];
			ofs_minus << std::setw(8) << std::fixed << kpoint->kaxis[ik];

			for (is = 0; is < ns; ++is){
				ofs_plus << std::setw(15) << std::scientific << writes->in_kayser(eval_plus[ik][is]);
				ofs_minus << std::setw(15) << std::scientific << writes->in_kayser(eval_minus[ik][is]);
			}
			ofs_plus << std::endl;
			ofs_minus << std::endl;
		}

		ofs_plus.close();
		ofs_minus.close();
	}


	memory->deallocate(evec_tmp);
	memory->deallocate(eval_orig);
	memory->deallocate(eval_plus);
	memory->deallocate(eval_minus);

	std::cout << "done !" << std::endl;
}

void Gruneisen::finish_gruneisen()
{
	memory->deallocate(gruneisen);
	memory->deallocate(dfc2);

	if (fcs_phonon->is_fc2_ext) {
		fc2_plus_ext.clear();
		fc2_minus_ext.clear();
	} else {
		memory->deallocate(fc2_plus);
		memory->deallocate(fc2_minus);
	}

}

void Gruneisen::prepare_delta_fc2()
{
	unsigned int i;
	unsigned int iat, jat, kat;
	unsigned int icrd, jcrd;
	double coord_tmp[3];

	unsigned int nalpha[3], ncell[3], ixyz[3];

	unsigned int natmin = system->natmin;
	unsigned int nat = system->nat;

	std::cout << "Preparing delta FC2 from cubic force constants ...";

	memory->allocate(dfc2, natmin, nat, 3, 3);

	for (iat = 0; iat < natmin; ++iat){
		for (jat = 0; jat < nat; ++jat){
			for (icrd = 0; icrd < 3; ++icrd){
				for (jcrd = 0; jcrd < 3; ++jcrd){
					dfc2[iat][jat][icrd][jcrd] = 0.0;
				}
			}
		}
	}

	for (std::vector<FcsClass>::iterator it = fcs_phonon->force_constant[1].begin(); it != fcs_phonon->force_constant[1].end(); ++it) {

		FcsClass fc3_tmp = *it;

		for (i = 0; i < 3; ++i){
			nalpha[i] = fc3_tmp.elems[i].atom;
			ncell[i]  = fc3_tmp.elems[i].cell;
			ixyz[i] = fc3_tmp.elems[i].xyz;
		}

		iat = system->map_p2s[nalpha[0]][ncell[0]];
		jat = system->map_p2s[nalpha[1]][ncell[1]];
		kat = system->map_p2s[nalpha[2]][ncell[2]];

		for (i = 0; i < 3; ++i) {
			coord_tmp[i] = system->xr_s[kat][i] - system->xr_s[iat][i];
			coord_tmp[i] = dynamical->fold(coord_tmp[i]);
		}
		rotvec(coord_tmp, coord_tmp, system->lavec_s);

		dfc2[nalpha[0]][jat][icrd][jcrd] += fc3_tmp.fcs_val * coord_tmp[ixyz[2]];
	}

#ifdef _DEBUG
	for (iat = 0; iat < natmin; ++iat){
		for (jat = 0; jat < nat; ++jat){
			for (icrd = 0; icrd < 3; ++icrd){
				for (jcrd = 0; jcrd < 3; ++jcrd){
					std::cout << std::setw(15) << dfc2[iat][jat][icrd][jcrd] << std::endl;
				}
			}
		}
	}
#endif

	std::cout << "done !" << std::endl;
}<commit_msg>A critical bug in gruneisen.cpp was fixed.<commit_after>#include "mpi_common.h"
#include <iostream>
#include <iomanip>
#include <fstream>
#include "dynamical.h"
#include "error.h"
#include "fcs_phonon.h"
#include "gruneisen.h"
#include "pointers.h"
#include "kpoint.h"
#include "memory.h"
#include "system.h"
#include "parsephon.h"
#include "write_phonons.h"
#include "../alm_c++/mathfunctions.h"
#include "relaxation.h"

using namespace PHON_NS;

Gruneisen::Gruneisen(PHON *phon): Pointers(phon){};
Gruneisen::~Gruneisen(){};

void Gruneisen::setup()
{
	//	MPI_Bcast(&delta_a, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

	//	if (mympi->my_rank == 0) {
	std::cout << std::endl; 
	std::cout << "Change in cell parameters : " << delta_a << std::endl;
	//	}

	memory->allocate(gruneisen, kpoint->nk, dynamical->neval);

	prepare_delta_fc2();
	prepare_newfc2();

	if (relaxation->quartic_mode) {
		prepare_newfc3();
		write_newinfo_all();
	}
}

void Gruneisen::calc_gruneisen()
{
	std::cout << "Calculating Gruneisen parameters ..." << std::endl;

	unsigned int nk = kpoint->nk;
	unsigned int ns = dynamical->neval;

	unsigned int i;
	unsigned int ik, is;

	double *eval_orig;
	double **eval_plus;
	double **eval_minus;

	double xk_tmp[3];


	std::complex<double> **evec_tmp;

	memory->allocate(evec_tmp, 1, 1); // dummy allocation
	memory->allocate(eval_orig, ns);
	memory->allocate(eval_plus, nk, ns);
	memory->allocate(eval_minus, nk, ns);

	for (ik = 0; ik < nk; ++ik){

		for (i = 0; i < 3; ++i) xk_tmp[i] = kpoint->xk[ik][i];

		if (fcs_phonon->is_fc2_ext) {
			dynamical->eval_k(xk_tmp, dynamical->kvec_na[ik], fcs_phonon->fc2_ext, eval_orig, evec_tmp, false);
			dynamical->eval_k(xk_tmp, dynamical->kvec_na[ik], fc2_plus_ext, eval_plus[ik], evec_tmp, false);
			dynamical->eval_k(xk_tmp, dynamical->kvec_na[ik], fc2_minus_ext, eval_minus[ik], evec_tmp, false);
		} else {
			dynamical->eval_k(xk_tmp, dynamical->kvec_na[ik], fcs_phonon->fc2, eval_orig, evec_tmp, false);
			dynamical->eval_k(xk_tmp, dynamical->kvec_na[ik], fc2_plus, eval_plus[ik], evec_tmp, false);
			dynamical->eval_k(xk_tmp, dynamical->kvec_na[ik], fc2_minus, eval_minus[ik], evec_tmp, false);
		}


		for (is = 0; is < ns; ++is) {
			gruneisen[ik][is] = (eval_plus[ik][is] - eval_minus[ik][is]) / (2.0 * delta_a) / (-6.0 * eval_orig[is]);
		}
	}

	for (ik = 0; ik < nk; ++ik) {
		for (is = 0; is < ns; ++is) {
			eval_plus[ik][is] = dynamical->freq(eval_plus[ik][is]);
			eval_minus[ik][is] = dynamical->freq(eval_minus[ik][is]);
		}
	}

	if (kpoint->kpoint_mode == 1) {
		std::string file_band_plus, file_band_minus;

		file_band_plus = input->job_title + ".band_+";
		file_band_minus = input->job_title + ".band_-";

		std::ofstream ofs_plus, ofs_minus;

		ofs_plus.open(file_band_plus.c_str(), std::ios::out);
		if (!ofs_plus) error->exit("calc_gruneisen", "Could not create band_plus file.");
		ofs_minus.open(file_band_minus.c_str(), std::ios::out);
		if (!ofs_minus) error->exit("calc_gruneisen", "Could not create band_minus file.");

		ofs_plus << "# Phonon energy (cm^-1) of the system expanded by " << std::setw(10) << delta_a * 100 << " %." << std::endl;
		ofs_minus << "# Phonon energy (cm^-1) of the system compressed by " << std::setw(10) << delta_a * 100 << " %." << std::endl;

		for (ik = 0; ik < nk; ++ik){
			ofs_plus << std::setw(8) << std::fixed << kpoint->kaxis[ik];
			ofs_minus << std::setw(8) << std::fixed << kpoint->kaxis[ik];

			for (is = 0; is < ns; ++is){
				ofs_plus << std::setw(15) << std::scientific << writes->in_kayser(eval_plus[ik][is]);
				ofs_minus << std::setw(15) << std::scientific << writes->in_kayser(eval_minus[ik][is]);
			}
			ofs_plus << std::endl;
			ofs_minus << std::endl;
		}

		ofs_plus.close();
		ofs_minus.close();
	}


	memory->deallocate(evec_tmp);
	memory->deallocate(eval_orig);
	memory->deallocate(eval_plus);
	memory->deallocate(eval_minus);

	std::cout << "done !" << std::endl;
}

void Gruneisen::finish_gruneisen()
{
	memory->deallocate(gruneisen);
	memory->deallocate(dfc2);

	if (fcs_phonon->is_fc2_ext) {
		fc2_plus_ext.clear();
		fc2_minus_ext.clear();
	} else {
		memory->deallocate(fc2_plus);
		memory->deallocate(fc2_minus);
	}

}

void Gruneisen::prepare_delta_fc2()
{
	unsigned int i;
	unsigned int iat, jat, kat;
	unsigned int icrd, jcrd;
	double coord_tmp[3];

	unsigned int nalpha[3], ncell[3], ixyz[3];

	unsigned int natmin = system->natmin;
	unsigned int nat = system->nat;

	std::cout << "Preparing delta FC2 from cubic force constants ...";

	memory->allocate(dfc2, natmin, nat, 3, 3);

	for (iat = 0; iat < natmin; ++iat){
		for (jat = 0; jat < nat; ++jat){
			for (icrd = 0; icrd < 3; ++icrd){
				for (jcrd = 0; jcrd < 3; ++jcrd){
					dfc2[iat][jat][icrd][jcrd] = 0.0;
				}
			}
		}
	}

	for (std::vector<FcsClass>::iterator it = fcs_phonon->force_constant[1].begin(); it != fcs_phonon->force_constant[1].end(); ++it) {

		FcsClass fc3_tmp = *it;

		for (i = 0; i < 3; ++i){
			nalpha[i] = fc3_tmp.elems[i].atom;
			ncell[i]  = fc3_tmp.elems[i].cell;
			ixyz[i] = fc3_tmp.elems[i].xyz;
		}

		iat = system->map_p2s[nalpha[0]][ncell[0]];
		jat = system->map_p2s[nalpha[1]][ncell[1]];
		kat = system->map_p2s[nalpha[2]][ncell[2]];

		for (i = 0; i < 3; ++i) {
			coord_tmp[i] = system->xr_s[kat][i] - system->xr_s[iat][i];
			coord_tmp[i] = dynamical->fold(coord_tmp[i]);
		}
		rotvec(coord_tmp, coord_tmp, system->lavec_s);

		dfc2[nalpha[0]][jat][ixyz[0]][ixyz[1]] += fc3_tmp.fcs_val * coord_tmp[ixyz[2]];
	}

#ifdef _DEBUG
	for (iat = 0; iat < natmin; ++iat){
		for (jat = 0; jat < nat; ++jat){
			for (icrd = 0; icrd < 3; ++icrd){
				for (jcrd = 0; jcrd < 3; ++jcrd){
					std::cout << std::setw(15) << dfc2[iat][jat][icrd][jcrd] << std::endl;
				}
			}
		}
	}
#endif

	std::cout << "done !" << std::endl;
}
{'fraction_non_alphanumeric': 0.1994262366339216, 'fraction_numerical': 0.01808223941580457, 'mean_word_length': 4.485932284215546, 'pattern_counts': {'":': 0, '<': 199, '<?xml version=': 0, '>': 131, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6472657', 'n_tokens_mistral': 5403, 'n_tokens_neox': 4854, 'n_words': 914}