Skip to content

Commit

Permalink
Merge branch 'master' into develop
Browse files Browse the repository at this point in the history
  • Loading branch information
@ttadano
ttadano committed Jun 26, 2015
2 parents d2199fb + 4ad1574 commit d1be0b4
Show file tree
Hide file tree
Showing 27 changed files with 431 additions and 172 deletions.
2 changes: 1 addition & 1 deletion alm/alamode.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ ALM::ALM(int narg, char **arg)
{
std::cout << " +------------------------------------------------------------+" << std::endl;
std::cout << " + Program ALM +" << std::endl;
std::cout << " + Ver. 0.9.4 +" << std::endl;
std::cout << " + Ver. 0.9.5 +" << std::endl;
std::cout << " +------------------------------------------------------------+" << std::endl;
std::cout << std::endl;

Expand Down
109 changes: 102 additions & 7 deletions alm/constraint.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -232,7 +232,7 @@ void Constraint::setup()
memory->allocate(const_relate, maxorder);
memory->allocate(index_bimap, maxorder);

get_mapping_constraint(maxorder, const_self, const_fix, const_relate, index_bimap);
get_mapping_constraint(maxorder, const_self, const_fix, const_relate, index_bimap, false);

for (order = 0; order < maxorder; ++order) {
std::cout << " Number of free" << std::setw(9) << interaction->str_order[order] << " FCs : " << index_bimap[order].size() << std::endl;
Expand Down Expand Up @@ -425,7 +425,7 @@ void Constraint::calc_constraint_matrix(const int N, int &P)
void Constraint::get_mapping_constraint(const int nmax, std::set<ConstraintClass> *const_in,
std::vector<ConstraintTypeFix> *const_fix_out,
std::vector<ConstraintTypeRelate> *const_relate_out,
boost::bimap<int, int> *index_bimap_out)
boost::bimap<int, int> *index_bimap_out, const bool is_suggest_mode)
{
int order;
unsigned int i;
Expand All @@ -441,10 +441,13 @@ void Constraint::get_mapping_constraint(const int nmax, std::set<ConstraintClass
fix_forceconstant[i] = false;
file_forceconstant[i] = "";
}
fix_forceconstant[0] = fix_harmonic;
fix_forceconstant[1] = fix_cubic;
if (fix_forceconstant[0]) file_forceconstant[0] = fc2_file;
if (fix_forceconstant[1]) file_forceconstant[1] = fc3_file;

if (!is_suggest_mode) {
fix_forceconstant[0] = fix_harmonic;
fix_forceconstant[1] = fix_cubic;
if (fix_forceconstant[0]) file_forceconstant[0] = fc2_file;
if (fix_forceconstant[1]) file_forceconstant[1] = fc3_file;
}

int nparam;
for (order = 0; order < nmax; ++order) {
Expand Down Expand Up @@ -701,6 +704,7 @@ void Constraint::translational_invariance()

int i, j;
int iat, jat, icrd, jcrd;
int idata;
int order;
int maxorder = interaction->maxorder;

Expand All @@ -720,6 +724,7 @@ void Constraint::translational_invariance()
std::vector<int> intlist, data;
std::set<FcProperty> list_found;
std::set<FcProperty>::iterator iter_found;
std::vector<std::vector<int> > data_vec;

std::cout << " Generating constraints for translational invariance ..." << std::endl;

Expand Down Expand Up @@ -799,6 +804,91 @@ void Constraint::translational_invariance()
}
std::sort(intlist.begin(), intlist.end());

data_vec.clear();

CombinationWithRepetition<int> g2(intlist.begin(), intlist.end(), order);
do {
data = g2.now();

intarr[0] = iat;
for (isize = 0; isize < data.size(); ++isize) {
intarr[isize + 1] = data[isize];
}

// Skip if the atoms don't interact with each other.
if (!interaction->is_incutoff(order + 1, intarr)) continue;

data_vec.push_back(data);
} while(g2.next());

int ndata = data_vec.size();

#pragma omp parallel
{
int *intarr_omp, *intarr_copy_omp;
double *arr_constraint_omp;

memory->allocate(intarr_omp, order + 2);
memory->allocate(intarr_copy_omp, order + 2);
memory->allocate(arr_constraint_omp, nparams);

std::vector<int> data_omp;

#pragma omp for private(isize, ixyz, jcrd, j, jat, iter_found), schedule(guided)
for (idata = 0; idata < ndata; ++idata) {

data_omp = data_vec[idata];

intarr_omp[0] = iat;
for (isize = 0; isize < data_omp.size(); ++isize) {
intarr_omp[isize + 1] = data_omp[isize];
}
if (!interaction->is_incutoff(order + 1, intarr_omp)) continue;


// Loop for xyz component
for (ixyz = 0; ixyz < nxyz; ++ixyz) {
for (jcrd = 0; jcrd < 3; ++jcrd) {

// Reset the temporary array for another constraint
for (j = 0; j < nparams; ++j) arr_constraint_omp[j] = 0.0;

for (jat = 0; jat < 3 * nat; jat += 3) {
intarr_omp[order + 1] = jat / 3;

if (!interaction->is_incutoff(order + 2, intarr_omp)) continue;

for (j = 0; j < order + 1; ++j) {
intarr_copy_omp[j] = 3 * intarr_omp[j] + xyzcomponent[ixyz][j];
}
intarr_copy_omp[order + 1] = jat + jcrd;

fcs->sort_tail(order + 2, intarr_copy_omp);

iter_found = list_found.find(FcProperty(order + 2, 1.0, intarr_copy_omp, 1));
if (iter_found != list_found.end()) {
FcProperty arrtmp = *iter_found;
arr_constraint_omp[arrtmp.mother] += arrtmp.coef;
}

}
if (!is_allzero(nparams,arr_constraint_omp)) {
#pragma omp critical
const_translation[order].insert(ConstraintClass(nparams, arr_constraint_omp));
}
}
}


}

memory->deallocate(intarr_omp);
memory->deallocate(intarr_copy_omp);
memory->deallocate(arr_constraint_omp);
}

/*
CombinationWithRepetition<int> g(intlist.begin(), intlist.end(), order);
do {
data = g.now();
Expand Down Expand Up @@ -842,6 +932,7 @@ void Constraint::translational_invariance()
}
} while(g.next());
*/
intlist.clear();
}
}
Expand Down Expand Up @@ -1325,7 +1416,11 @@ void Constraint::rotational_invariance()

void Constraint::remove_redundant_rows(const int n, std::set<ConstraintClass> &Constraint_Set, const double tolerance)
{
#ifdef _USE_EIGEN
#ifdef _USE_EIGEN_DISABLE

// This function doesn't make the reduced row echelon form of the constraint matrix.
// It just returns the image of the matrix, though they are similar.

using namespace Eigen;

int nrow = n;
Expand Down
9 changes: 5 additions & 4 deletions alm/constraint.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,11 +90,14 @@ namespace ALM_NS
std::string rotation_axis;
std::set<ConstraintClass> *const_symmetry;

void constraint_from_symmetry(std::set<ConstraintClass> *);
std::vector<ConstraintTypeFix> *const_fix;
std::vector<ConstraintTypeRelate> *const_relate;
boost::bimap<int, int> *index_bimap;

void constraint_from_symmetry(std::set<ConstraintClass> *);
void get_mapping_constraint(const int, std::set<ConstraintClass> *,
std::vector<ConstraintTypeFix> *, std::vector<ConstraintTypeRelate> *,
boost::bimap<int, int> *, const bool);

private:

Expand All @@ -113,9 +116,7 @@ namespace ALM_NS
void translational_invariance();
void rotational_invariance();
void calc_constraint_matrix(const int, int &);
void get_mapping_constraint(const int, std::set<ConstraintClass> *,
std::vector<ConstraintTypeFix> *, std::vector<ConstraintTypeRelate> *,
boost::bimap<int, int> *);

void setup_rotation_axis(bool [3][3]);
bool is_allzero(const int, const double *, const int nshift = 0);
void remove_redundant_rows(const int, std::set<ConstraintClass> &, const double tolerance = eps12);
Expand Down
5 changes: 1 addition & 4 deletions alm/fcs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,8 +88,6 @@ void Fcs::generate_fclists(int maxorder)
int nxyz;
unsigned int isym;

IntList list_tmp;

double c_tmp;

int **xyzcomponent;
Expand Down Expand Up @@ -127,8 +125,7 @@ void Fcs::generate_fclists(int maxorder)

for (std::set<IntList>::iterator iter = interaction->pairs[order].begin(); iter != interaction->pairs[order].end(); ++iter) {

IntList list_tmp = *iter;
for (i = 0; i < order + 2; ++i) atmn[i] = list_tmp.iarray[i];
for (i = 0; i < order + 2; ++i) atmn[i] = (*iter).iarray[i];

for (i1 = 0; i1 < nxyz; ++i1) {
for (i = 0; i < order + 2; ++i) ind[i] = 3 * atmn[i] + xyzcomponent[i1][i];
Expand Down
34 changes: 27 additions & 7 deletions alm/fitting.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,8 @@ void Fitting::fitmain()
std::cout << " NSTART = " << nstart << "; NEND = " << nend << std::endl;
std::cout << " " << nend - nstart + 1 << " entries will be used for fitting." << std::endl << std::endl;

data_multiplier(nat, ndata, nstart, nend, ndata_used, nmulti, symmetry->multiply_data, u, f);
data_multiplier(nat, ndata, nstart, nend, ndata_used, nmulti, symmetry->multiply_data, u, f,
files->file_disp, files->file_force);

N = 0;
for (i = 0; i < maxorder; ++i) {
Expand Down Expand Up @@ -208,7 +209,8 @@ void Fitting::fitmain()
}

void Fitting::data_multiplier(const int nat, const int ndata, const int nstart, const int nend,
const int ndata_used, int &nmulti, const int multiply_data, double **&u, double **&f)
const int ndata_used, int &nmulti, const int multiply_data, double **&u, double **&f,
const std::string file_disp, const std::string file_force)
{
int i, j, k;
int idata, itran, isym;
Expand All @@ -220,6 +222,13 @@ void Fitting::data_multiplier(const int nat, const int ndata, const int nstart,
unsigned int nline_f, nline_u;
unsigned int nreq;

std::ifstream ifs_disp, ifs_force;

ifs_disp.open(file_disp.c_str(), std::ios::in);
if (!ifs_disp) error->exit("openfiles", "cannot open disp file");
ifs_force.open(file_force.c_str(), std::ios::in);
if (!ifs_force) error->exit("openfiles", "cannot open force file");

nreq = 3 * nat * ndata;

memory->allocate(u_tmp, nreq);
Expand All @@ -228,7 +237,7 @@ void Fitting::data_multiplier(const int nat, const int ndata, const int nstart,
// Read displacements from DFILE

nline_u = 0;
while (files->ifs_disp >> u_in) {
while (ifs_disp >> u_in) {
u_tmp[nline_u++] = u_in;
if (nline_u == nreq) break;
}
Expand All @@ -238,7 +247,7 @@ void Fitting::data_multiplier(const int nat, const int ndata, const int nstart,
// Read forces from FFILE

nline_f = 0;
while (files->ifs_force >> f_in) {
while (ifs_force >> f_in) {
f_tmp[nline_f++] = f_in;
if (nline_f == nreq) break;
}
Expand Down Expand Up @@ -344,6 +353,9 @@ void Fitting::data_multiplier(const int nat, const int ndata, const int nstart,

memory->deallocate(u_tmp);
memory->deallocate(f_tmp);

ifs_disp.close();
ifs_force.close();
}

void Fitting::fit_without_constraints(int N, int M_Start, int M_End, double **amat, double *bvec)
Expand Down Expand Up @@ -969,10 +981,14 @@ void Fitting::calc_matrix_elements_algebraic_constraint(const int M, const int N
int irow;
int ncycle;


std::cout << " Calculation of matrix elements for direct fitting started ... ";

ncycle = ndata_fit * nmulti;


#ifdef _OPENMP
#pragma omp parallel for private(j)
#endif
for (i = 0; i < M; ++i) {
for (j = 0; j < N_new; ++j) {
amat[i][j] = 0.0;
Expand All @@ -981,8 +997,9 @@ void Fitting::calc_matrix_elements_algebraic_constraint(const int M, const int N
bvec_orig[i] = 0.0;
}

ncycle = ndata_fit * nmulti;

#ifdef _OPENMP
#pragma omp parallel private(irow, i, j)
#endif
{
int *ind;
int mm, order, iat, k;
Expand All @@ -997,6 +1014,9 @@ void Fitting::calc_matrix_elements_algebraic_constraint(const int M, const int N
memory->allocate(amat_orig, 3 * natmin, N);
memory->allocate(amat_mod, 3 * natmin, N_new);

#ifdef _OPENMP
#pragma omp for schedule(guided)
#endif
for (irow = 0; irow < ncycle; ++irow) {

// generate r.h.s vector B
Expand Down
13 changes: 9 additions & 4 deletions alm/fitting.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,16 +31,22 @@ namespace ALM_NS {
unsigned int nboot;
unsigned int seed;

void data_multiplier(const int, const int, const int, const int, const int, int &, const int, double **&, double **&,
const std::string, const std::string);
void calc_matrix_elements_algebraic_constraint(const int, const int, const int, const int,
const int, const int, const int, const int, double **, double **, double **, double *, double *);

#ifdef _VSL
VSLStreamStatePtr stream;
int brng;
#endif

double gamma(const int, const int *);

private:

int inprim_index(const int);
void wrtfcs(const double *);
void data_multiplier(const int, const int, const int, const int, const int, int &, const int, double **&, double **&);
void fit_without_constraints(int, int, int, double **, double *);
void fit_algebraic_constraints(int, int, int, double **, double *, double *, const int);

Expand All @@ -49,10 +55,9 @@ namespace ALM_NS {
const int, const int, double **, double *, double **, double *);
void calc_matrix_elements(const int, const int, const int,
const int, const int, const int, const int, double **, double **, double **, double *);
void calc_matrix_elements_algebraic_constraint(const int, const int, const int, const int,
const int, const int, const int, const int, double **, double **, double **, double *, double *);

void fit_bootstrap(int, int, int, int, int, double **, double *, double **, double *);
double gamma(const int, const int *);
// double gamma(const int, const int *);
int factorial(const int);
int rankSVD(const int, const int, double *, const double);
int rankQRD(const int, const int, double *, const double);
Expand Down

0 comments on commit d1be0b4

Please sign in to comment.