diff --git a/include/tadah/mlip/models/m_blr.h b/include/tadah/mlip/models/m_blr.h
index eafc784ea1734ddbc6028c99637067fcdfc39b1f..82ebbffa15b34b1fb2d404a98dada21d796f5b89 100644
--- a/include/tadah/mlip/models/m_blr.h
+++ b/include/tadah/mlip/models/m_blr.h
@@ -71,9 +71,9 @@ template
<class BF=DM_Function_Base&>
class M_BLR: public M_Tadah_Base, public M_BLR_Train<BF> {
- public:
+public:
- /** This constructor will preapare this object for either training
+ /** This constructor will preapare this object for either training
* or prediction (if potential is provides as a Config)
*
* Usage example:
@@ -84,16 +84,16 @@ class M_BLR: public M_Tadah_Base, public M_BLR_Train<BF> {
* \endcode
*
*/
- using M_BLR_Train<BF>::config;
- using M_BLR_Train<BF>::bf;
- M_BLR(Config &c):
- M_BLR_Train<BF>(c),
- desmat(M_BLR_Train<BF>::bf,c)
+ using M_BLR_Train<BF>::config;
+ using M_BLR_Train<BF>::bf;
+ M_BLR(Config &c):
+ M_BLR_Train<BF>(c),
+ desmat(M_BLR_Train<BF>::bf,c)
{
norm = Normaliser(c);
}
- /** This constructor will preapare this object for either training
+ /** This constructor will preapare this object for either training
* or prediction (if potential is provides as a Config)
*
* Usage example:
@@ -105,260 +105,260 @@ class M_BLR: public M_Tadah_Base, public M_BLR_Train<BF> {
* \endcode
*
*/
- M_BLR(BF &bf, Config &c):
- M_BLR_Train<BF>(bf,c),
- desmat(bf,c)
+ M_BLR(BF &bf, Config &c):
+ M_BLR_Train<BF>(bf,c),
+ desmat(bf,c)
{
norm = Normaliser(c);
}
- double epredict(const aed_type2 &aed) const{
- return bf.epredict(weights,aed);
- };
+ double epredict(const aed_type2 &aed) const{
+ return bf.epredict(weights,aed);
+ };
- double fpredict(const fd_type &fdij, const aed_type2 &aedi, const size_t k) const{
- return bf.fpredict(weights,fdij,aedi,k);
- }
-
- force_type fpredict(const fd_type &fdij, const aed_type2 &aedi) const{
- return bf.fpredict(weights,fdij,aedi);
- }
+ double fpredict(const fd_type &fdij, const aed_type2 &aedi, const size_t k) const{
+ return bf.fpredict(weights,fdij,aedi,k);
+ }
- void train(StDescriptorsDB &st_desc_db, const StructureDB &stdb) {
+ force_type fpredict(const fd_type &fdij, const aed_type2 &aedi) const{
+ return bf.fpredict(weights,fdij,aedi);
+ }
- if(config.template get<bool>("NORM"))
- norm = Normaliser(config,st_desc_db);
+ void train(StDescriptorsDB &st_desc_db, const StructureDB &stdb) {
- desmat.build(st_desc_db,stdb);
- train(desmat);
- }
+ if(config.template get<bool>("NORM"))
+ norm = Normaliser(config,st_desc_db);
- void train(StructureDB &stdb, DC_Base &dc) {
+ desmat.build(st_desc_db,stdb);
+ train(desmat);
+ }
- if(config.template get<bool>("NORM")) {
+ void train(StructureDB &stdb, DC_Base &dc) {
- std::string force=config.template get<std::string>("FORCE");
- std::string stress=config.template get<std::string>("STRESS");
+ if(config.template get<bool>("NORM")) {
- config.remove("FORCE");
- config.remove("STRESS");
- config.add("FORCE", "false");
- config.add("STRESS", "false");
+ std::string force=config.template get<std::string>("FORCE");
+ std::string stress=config.template get<std::string>("STRESS");
- StDescriptorsDB st_desc_db_temp = dc.calc(stdb);
+ config.remove("FORCE");
+ config.remove("STRESS");
+ config.add("FORCE", "false");
+ config.add("STRESS", "false");
- if(config.template get<bool>("NORM")) {
- norm = Normaliser(config);
- norm.learn(st_desc_db_temp);
- // norm.normalise(st_desc_db_temp);
- }
+ StDescriptorsDB st_desc_db_temp = dc.calc(stdb);
- config.remove("FORCE");
- config.remove("STRESS");
- config.add("FORCE", force);
- config.add("STRESS", stress);
+ if(config.template get<bool>("NORM")) {
+ norm = Normaliser(config);
+ norm.learn(st_desc_db_temp);
+ // norm.normalise(st_desc_db_temp);
}
- desmat.build(stdb,norm,dc);
- train(desmat);
+ config.remove("FORCE");
+ config.remove("STRESS");
+ config.add("FORCE", force);
+ config.add("STRESS", stress);
}
- Structure predict(const Config &c, StDescriptors &std, const Structure &st) {
- if(config.template get<bool>("NORM") && !std.normalised && bf.get_label()!="BF_Linear")
- norm.normalise(std);
- return M_Tadah_Base::predict(c,std,st);
- }
+ desmat.build(stdb,norm,dc);
+ train(desmat);
+ }
- StructureDB predict(Config &c, const StructureDB &stdb, DC_Base &dc) {
- return M_Tadah_Base::predict(c,stdb,dc);
+ Structure predict(const Config &c, StDescriptors &std, const Structure &st) {
+ if(config.template get<bool>("NORM") && !std.normalised && bf.get_label()!="BF_Linear")
+ norm.normalise(std);
+ return M_Tadah_Base::predict(c,std,st);
+ }
+
+ StructureDB predict(Config &c, const StructureDB &stdb, DC_Base &dc) {
+ return M_Tadah_Base::predict(c,stdb,dc);
+ }
+
+ Config get_param_file() {
+ Config c = config;
+ //c.remove("ALPHA");
+ //c.remove("BETA");
+ c.remove("DBFILE");
+ c.remove("FORCE");
+ c.remove("STRESS");
+ c.remove("VERBOSE");
+ c.add("VERBOSE", 0);
+
+ c.clear_internal_keys();
+ c.remove("MODEL");
+ c.add("MODEL", label);
+ c.add("MODEL", bf.get_label());
+
+ for (size_t i=0;i<weights.size();++i) {
+ c.add("WEIGHTS", weights(i));
}
- Config get_param_file() {
- Config c = config;
- //c.remove("ALPHA");
- //c.remove("BETA");
- c.remove("DBFILE");
- c.remove("FORCE");
- c.remove("STRESS");
- c.remove("VERBOSE");
- c.add("VERBOSE", 0);
-
- c.clear_internal_keys();
- c.remove("MODEL");
- c.add("MODEL", label);
- c.add("MODEL", bf.get_label());
-
- for (size_t i=0;i<weights.size();++i) {
- c.add("WEIGHTS", weights(i));
+ if(config.template get<bool>("NORM")) {
+ for (size_t i=0;i<norm.mean.size();++i) {
+ c.add("NMEAN", norm.mean[i]);
}
-
- if(config.template get<bool>("NORM")) {
- for (size_t i=0;i<norm.mean.size();++i) {
- c.add("NMEAN", norm.mean[i]);
- }
- for (size_t i=0;i<norm.std_dev.size();++i) {
- c.add("NSTDEV", norm.std_dev[i]);
- }
+ for (size_t i=0;i<norm.std_dev.size();++i) {
+ c.add("NSTDEV", norm.std_dev[i]);
}
- return c;
}
- StructureDB predict(Config config_pred, StructureDB &stdb, DC_Base &dc,
- aed_type2 &predicted_error) {
-
- LinearRegressor::read_sigma(config_pred,Sigma);
- DesignMatrix<BF> dm(bf,config_pred);
- dm.scale=false; // do not scale energy, forces and stresses
- dm.build(stdb,norm,dc);
-
- predicted_error = T_MDMT_diag(dm.Phi, Sigma);
- double pmean = sqrt(predicted_error.mean());
-
- // compute energy, forces and stresses
- aed_type2 Tpred = T_dgemv(dm.Phi, weights);
-
- // Construct StructureDB object with predicted values
- StructureDB stdb_;
- stdb_.structures.resize(stdb.size());
- size_t i=0;
- for (size_t s=0; s<stdb.size(); ++s) {
- stdb_(s) = Structure(stdb(s));
-
- predicted_error(i) = (sqrt(predicted_error(i))-pmean)/stdb(s).natoms();
- stdb_(s).energy = Tpred(i++);
- if (config_pred.get<bool>("FORCE")) {
- for (size_t a=0; a<stdb(s).natoms(); ++a) {
- for (size_t k=0; k<3; ++k) {
- stdb_(s).atoms[a].force[k] = Tpred(i++);
- predicted_error(i) = (sqrt(predicted_error(i))-pmean);
- }
+ return c;
+ }
+ StructureDB predict(Config config_pred, StructureDB &stdb, DC_Base &dc,
+ aed_type2 &predicted_error) {
+
+ LinearRegressor::read_sigma(config_pred,Sigma);
+ DesignMatrix<BF> dm(bf,config_pred);
+ dm.scale=false; // do not scale energy, forces and stresses
+ dm.build(stdb,norm,dc);
+
+ predicted_error = T_MDMT_diag(dm.Phi, Sigma);
+ double pmean = sqrt(predicted_error.mean());
+
+ // compute energy, forces and stresses
+ aed_type2 Tpred = T_dgemv(dm.Phi, weights);
+
+ // Construct StructureDB object with predicted values
+ StructureDB stdb_;
+ stdb_.structures.resize(stdb.size());
+ size_t i=0;
+ for (size_t s=0; s<stdb.size(); ++s) {
+ stdb_(s) = Structure(stdb(s));
+
+ predicted_error(i) = (sqrt(predicted_error(i))-pmean)/stdb(s).natoms();
+ stdb_(s).energy = Tpred(i++);
+ if (config_pred.get<bool>("FORCE")) {
+ for (size_t a=0; a<stdb(s).natoms(); ++a) {
+ for (size_t k=0; k<3; ++k) {
+ stdb_(s).atoms[a].force[k] = Tpred(i++);
+ predicted_error(i) = (sqrt(predicted_error(i))-pmean);
}
}
- if (config_pred.get<bool>("STRESS")) {
- for (size_t x=0; x<3; ++x) {
- for (size_t y=x; y<3; ++y) {
- stdb_(s).stress(x,y) = Tpred(i++);
- predicted_error(i) = (sqrt(predicted_error(i))-pmean);
- if (x!=y)
- stdb_(s).stress(y,x) = stdb_(s).stress(x,y);
- }
+ }
+ if (config_pred.get<bool>("STRESS")) {
+ for (size_t x=0; x<3; ++x) {
+ for (size_t y=x; y<3; ++y) {
+ stdb_(s).stress(x,y) = Tpred(i++);
+ predicted_error(i) = (sqrt(predicted_error(i))-pmean);
+ if (x!=y)
+ stdb_(s).stress(y,x) = stdb_(s).stress(x,y);
}
}
}
- return stdb_;
}
- StructureDB predict(StructureDB &stdb) {
- if(!trained) throw std::runtime_error("This object is not trained!\n\
- Hint: check different predict() methods.");
-
- phi_type &Phi = desmat.Phi;
- //std::cout << Phi.row(0) << std::endl;
-
- // compute energy, forces and stresses
- aed_type2 Tpred = T_dgemv(Phi, weights);
-
- double eweightglob=config.template get<double>("EWEIGHT");
- double fweightglob=config.template get<double>("FWEIGHT");
- double sweightglob=config.template get<double>("SWEIGHT");
-
- // Construct StructureDB object with predicted values
- StructureDB stdb_;
- stdb_.structures.resize(stdb.size());
- size_t s=0;
- size_t i=0;
- while (i<Phi.rows()) {
-
- stdb_(s).energy = Tpred(i++)*stdb(s).natoms()/eweightglob/stdb(s).eweight;
- if (config.template get<bool>("FORCE")) {
- stdb_(s).atoms.resize(stdb(s).natoms());
- for (size_t a=0; a<stdb(s).natoms(); ++a) {
- for (size_t k=0; k<3; ++k) {
- stdb_(s).atoms[a].force[k] = Tpred(i++)/fweightglob/stdb(s).fweight;
- }
+ return stdb_;
+ }
+ StructureDB predict(StructureDB &stdb) {
+ if(!trained) throw std::runtime_error("This object is not trained!\n\
+Hint: check different predict() methods.");
+
+ phi_type &Phi = desmat.Phi;
+ //std::cout << Phi.row(0) << std::endl;
+
+ // compute energy, forces and stresses
+ aed_type2 Tpred = T_dgemv(Phi, weights);
+
+ double eweightglob=config.template get<double>("EWEIGHT");
+ double fweightglob=config.template get<double>("FWEIGHT");
+ double sweightglob=config.template get<double>("SWEIGHT");
+
+ // Construct StructureDB object with predicted values
+ StructureDB stdb_;
+ stdb_.structures.resize(stdb.size());
+ size_t s=0;
+ size_t i=0;
+ while (i<Phi.rows()) {
+
+ stdb_(s).energy = Tpred(i++)*stdb(s).natoms()/eweightglob/stdb(s).eweight;
+ if (config.template get<bool>("FORCE")) {
+ stdb_(s).atoms.resize(stdb(s).natoms());
+ for (size_t a=0; a<stdb(s).natoms(); ++a) {
+ for (size_t k=0; k<3; ++k) {
+ stdb_(s).atoms[a].force[k] = Tpred(i++)/fweightglob/stdb(s).fweight;
}
}
- if (config.template get<bool>("STRESS")) {
- for (size_t x=0; x<3; ++x) {
- for (size_t y=x; y<3; ++y) {
- stdb_(s).stress(x,y) = Tpred(i++)/sweightglob/stdb(s).sweight;
- if (x!=y)
- stdb_(s).stress(y,x) = stdb_(s).stress(x,y);
- }
+ }
+ if (config.template get<bool>("STRESS")) {
+ for (size_t x=0; x<3; ++x) {
+ for (size_t y=x; y<3; ++y) {
+ stdb_(s).stress(x,y) = Tpred(i++)/sweightglob/stdb(s).sweight;
+ if (x!=y)
+ stdb_(s).stress(y,x) = stdb_(s).stress(x,y);
}
}
- s++;
}
- return stdb_;
+ s++;
}
+ return stdb_;
+ }
- private:
- std::string label="M_BLR";
- DesignMatrix<BF> desmat;
+private:
+ std::string label="M_BLR";
+ DesignMatrix<BF> desmat;
- // normalise weights such that when predict is called
- // we can supply it with a non-normalised descriptor
- t_type convert_to_nweights(const t_type &weights) const {
- if(bf.get_label()!="BF_Linear") {
- throw std::runtime_error("Cannot convert weights to nweights for\n\
- non linear basis function\n");
- }
- t_type nw(weights.rows());
- nw(0) = weights(0);
- for (size_t i=1; i<weights.size(); ++i) {
+ // normalise weights such that when predict is called
+ // we can supply it with a non-normalised descriptor
+ t_type convert_to_nweights(const t_type &weights) const {
+ if(bf.get_label()!="BF_Linear") {
+ throw std::runtime_error("Cannot convert weights to nweights for\n\
+non linear basis function\n");
+ }
+ t_type nw(weights.rows());
+ nw(0) = weights(0);
+ for (size_t i=1; i<weights.size(); ++i) {
- if (norm.std_dev[i] > std::numeric_limits<double>::min())
- nw(i) = weights(i) / norm.std_dev[i];
- else
- nw(i) = weights(i);
+ if (norm.std_dev[i] > std::numeric_limits<double>::min())
+ nw(i) = weights(i) / norm.std_dev[i];
+ else
+ nw(i) = weights(i);
- nw(0) -= norm.mean[i]*nw(i);
+ nw(0) -= norm.mean[i]*nw(i);
- }
- return nw;
}
- // The opposite of convert_to_nweights()
- t_type convert_to_weights(const t_type &nw) const {
- if(bf.get_label()!="BF_Linear") {
- throw std::runtime_error("Cannot convert nweights to weights for\n\
- non linear basis function\n");
- }
- // convert normalised weights back to "normal"
- t_type w(nw.rows());
- w(0) = nw(0);
- for (size_t i=1; i<nw.size(); ++i) {
- if (norm.std_dev[i] > std::numeric_limits<double>::min())
- w(i) = nw(i) * norm.std_dev[i];
- else
- w(i) = nw(i);
-
- w(0) += nw(i)*norm.mean[i];
- }
- return w;
+ return nw;
+ }
+ // The opposite of convert_to_nweights()
+ t_type convert_to_weights(const t_type &nw) const {
+ if(bf.get_label()!="BF_Linear") {
+ throw std::runtime_error("Cannot convert nweights to weights for\n\
+non linear basis function\n");
}
+ // convert normalised weights back to "normal"
+ t_type w(nw.rows());
+ w(0) = nw(0);
+ for (size_t i=1; i<nw.size(); ++i) {
+ if (norm.std_dev[i] > std::numeric_limits<double>::min())
+ w(i) = nw(i) * norm.std_dev[i];
+ else
+ w(i) = nw(i);
+
+ w(0) += nw(i)*norm.mean[i];
+ }
+ return w;
+ }
- template <typename D>
- void train(D &desmat) {
- // TODO
- // the train method destroys the Phi matrix
- // In consequence, we cannot use it for quick prediction
- // The simple solution, for now, is to make a copy of the Phi matrix
- //phi_type &Phi = desmat.Phi;
- phi_type Phi = desmat.Phi;
- t_type T = desmat.T;
- //t_type &T = desmat.T;
- M_BLR_Train<BF>::train(Phi,T);
-
- if (config.template get<bool>("NORM") &&
- bf.get_label()=="BF_Linear") {
- weights = convert_to_nweights(weights);
- }
- }
+ template <typename D>
+ void train(D &desmat) {
+ // TODO
+ // the train method destroys the Phi matrix
+ // In consequence, we cannot use it for quick prediction
+ // The simple solution, for now, is to make a copy of the Phi matrix
+ //phi_type &Phi = desmat.Phi;
+ phi_type Phi = desmat.Phi;
+ t_type T = desmat.T;
+ //t_type &T = desmat.T;
+ M_BLR_Train<BF>::train(Phi,T);
+
+ if (config.template get<bool>("NORM") &&
+ bf.get_label()=="BF_Linear") {
+ weights = convert_to_nweights(weights);
+ }
+ }
- // Do we want to confuse user with those and make them public?
- // Either way they must be stated as below to silence clang warning
- using M_BLR_Train<BF>::predict;
- using M_BLR_Train<BF>::train;
- using M_BLR_Train<BF>::trained;
- using M_BLR_Train<BF>::weights;
- using M_BLR_Train<BF>::Sigma;
+ // Do we want to confuse user with those and make them public?
+ // Either way they must be stated as below to silence clang warning
+ using M_BLR_Train<BF>::predict;
+ using M_BLR_Train<BF>::train;
+ using M_BLR_Train<BF>::trained;
+ using M_BLR_Train<BF>::weights;
+ using M_BLR_Train<BF>::Sigma;
};
#endif
diff --git a/include/tadah/mlip/models/m_krr.h b/include/tadah/mlip/models/m_krr.h
index 9062779b400a239743eb54f244205176b392d46f..6e79ae3363f5079b1be9c7cb9a2a7c68a002541a 100644
--- a/include/tadah/mlip/models/m_krr.h
+++ b/include/tadah/mlip/models/m_krr.h
@@ -49,12 +49,12 @@ template
<class K=DM_Function_Base&>
class M_KRR: public M_Tadah_Base,
public M_KRR_Train<K>
- //public M_KRR_Predict<K>
+//public M_KRR_Predict<K>
{
- public:
+public:
- /** This constructor will preapare this object for either training
+ /** This constructor will preapare this object for either training
* or prediction (if potential is provides as a Config)
*
* Usage example:
@@ -65,15 +65,15 @@ class M_KRR: public M_Tadah_Base,
* \endcode
*
*/
- M_KRR(Config &c):
- M_KRR_Train<K>(c),
- basis(c),
- desmat(kernel,c)
+ M_KRR(Config &c):
+ M_KRR_Train<K>(c),
+ basis(c),
+ desmat(kernel,c)
{
norm = Normaliser(c);
}
- /** This constructor will preapare this object for either training
+ /** This constructor will preapare this object for either training
* or prediction (if potential is provides as a Config)
*
* Usage example:
@@ -85,330 +85,330 @@ class M_KRR: public M_Tadah_Base,
* \endcode
*
*/
- M_KRR(K &kernel, Config &c):
- M_KRR_Train<K>(kernel,c),
- basis(c),
- desmat(kernel,c)
+ M_KRR(K &kernel, Config &c):
+ M_KRR_Train<K>(kernel,c),
+ basis(c),
+ desmat(kernel,c)
{
norm = Normaliser(c);
}
- double epredict(const aed_type2 &aed) const {
- return kernel.epredict(weights,aed);
- };
+ double epredict(const aed_type2 &aed) const {
+ return kernel.epredict(weights,aed);
+ };
- double fpredict(const fd_type &fdij, const aed_type2 &aedi, const size_t k) const {
- return kernel.fpredict(weights,fdij,aedi,k);
- }
+ double fpredict(const fd_type &fdij, const aed_type2 &aedi, const size_t k) const {
+ return kernel.fpredict(weights,fdij,aedi,k);
+ }
- force_type fpredict(const fd_type &fdij, const aed_type2 &aedi) const {
- return kernel.fpredict(weights,fdij,aedi);
- }
+ force_type fpredict(const fd_type &fdij, const aed_type2 &aedi) const {
+ return kernel.fpredict(weights,fdij,aedi);
+ }
- void train(StDescriptorsDB &st_desc_db, const StructureDB &stdb) {
+ void train(StDescriptorsDB &st_desc_db, const StructureDB &stdb) {
- if(config.template get<bool>("NORM"))
- norm = Normaliser(config,st_desc_db);
+ if(config.template get<bool>("NORM"))
+ norm = Normaliser(config,st_desc_db);
- desmat.build(st_desc_db,stdb);
- train(desmat);
- }
+ desmat.build(st_desc_db,stdb);
+ train(desmat);
+ }
- void train(StructureDB &stdb, DC_Base &dc) {
- int modelN;
- try {
- modelN=config.template get<int>("MODEL",2);
- }
- catch(std::runtime_error &e) {
- // use default model
- modelN=1;
- config.add("MODEL", modelN);
- }
- if (modelN==1)
- train1(stdb,dc);
- else if (modelN==2)
- train2(stdb,dc);
- else
- throw
- std::runtime_error(
- "This KRR implementation does exist: "\
- + std::to_string(modelN)+"\n");
+ void train(StructureDB &stdb, DC_Base &dc) {
+ int modelN;
+ try {
+ modelN=config.template get<int>("MODEL",2);
}
- void train1(StructureDB &stdb, DC_Base &dc) {
- // KRR implementation using EKM
- if(config.template get<bool>("NORM") || kernel.get_label()!="Kern_Linear") {
- // either build basis or prep normaliser
- std::string force=config.template get<std::string>("FORCE");
- std::string stress=config.template get<std::string>("STRESS");
-
- config.remove("FORCE");
- config.remove("STRESS");
- config.add("FORCE", "false");
- config.add("STRESS", "false");
-
- StDescriptorsDB st_desc_db_temp = dc.calc(stdb);
-
- if(config.template get<bool>("NORM")) {
- norm = Normaliser(config);
- norm.learn(st_desc_db_temp);
- norm.normalise(st_desc_db_temp);
- }
-
- config.remove("FORCE");
- config.remove("STRESS");
- config.add("FORCE", force);
- config.add("STRESS", stress);
-
- if (kernel.get_label()!="Kern_Linear") {
- basis.build_random_basis(config.template get<size_t>("SBASIS"),st_desc_db_temp);
- desmat.f.set_basis(basis.b);
- kernel.set_basis(basis.b);
- // to configure ekm, we need basis vectors
- ekm.configure(basis.b);
- }
- }
- desmat.build(stdb,norm,dc);
- train(desmat);
+ catch(std::runtime_error &e) {
+ // use default model
+ modelN=1;
+ config.add("MODEL", modelN);
}
-
- void train2(StructureDB &stdb, DC_Base &dc) {
-
- // NEW (BASIC) IMPLEMENTATION OF KRR //
+ if (modelN==1)
+ train1(stdb,dc);
+ else if (modelN==2)
+ train2(stdb,dc);
+ else
+ throw
+ std::runtime_error(
+ "This KRR implementation does exist: "\
+ + std::to_string(modelN)+"\n");
+ }
+ void train1(StructureDB &stdb, DC_Base &dc) {
+ // KRR implementation using EKM
+ if(config.template get<bool>("NORM") || kernel.get_label()!="Kern_Linear") {
+ // either build basis or prep normaliser
std::string force=config.template get<std::string>("FORCE");
std::string stress=config.template get<std::string>("STRESS");
+
config.remove("FORCE");
config.remove("STRESS");
config.add("FORCE", "false");
config.add("STRESS", "false");
+
StDescriptorsDB st_desc_db_temp = dc.calc(stdb);
+
if(config.template get<bool>("NORM")) {
norm = Normaliser(config);
norm.learn(st_desc_db_temp);
norm.normalise(st_desc_db_temp);
}
+
config.remove("FORCE");
config.remove("STRESS");
config.add("FORCE", force);
config.add("STRESS", stress);
- basis.prep_basis_for_krr(st_desc_db_temp,stdb);
- kernel.set_basis(basis.b);
- M_KRR_Train<K>::train2(basis.b, basis.T);
+
+ if (kernel.get_label()!="Kern_Linear") {
+ basis.build_random_basis(config.template get<size_t>("SBASIS"),st_desc_db_temp);
+ desmat.f.set_basis(basis.b);
+ kernel.set_basis(basis.b);
+ // to configure ekm, we need basis vectors
+ ekm.configure(basis.b);
+ }
}
+ desmat.build(stdb,norm,dc);
+ train(desmat);
+ }
- Structure predict(const Config &c, StDescriptors &std, const Structure &st) {
- if(config.template get<bool>("NORM") && !std.normalised && kernel.get_label()!="Kern_Linear")
- norm.normalise(std);
- return M_Tadah_Base::predict(c,std,st);
+ void train2(StructureDB &stdb, DC_Base &dc) {
+
+ // NEW (BASIC) IMPLEMENTATION OF KRR //
+ std::string force=config.template get<std::string>("FORCE");
+ std::string stress=config.template get<std::string>("STRESS");
+ config.remove("FORCE");
+ config.remove("STRESS");
+ config.add("FORCE", "false");
+ config.add("STRESS", "false");
+ StDescriptorsDB st_desc_db_temp = dc.calc(stdb);
+ if(config.template get<bool>("NORM")) {
+ norm = Normaliser(config);
+ norm.learn(st_desc_db_temp);
+ norm.normalise(st_desc_db_temp);
}
+ config.remove("FORCE");
+ config.remove("STRESS");
+ config.add("FORCE", force);
+ config.add("STRESS", stress);
+ basis.prep_basis_for_krr(st_desc_db_temp,stdb);
+ kernel.set_basis(basis.b);
+ M_KRR_Train<K>::train2(basis.b, basis.T);
+ }
+
+ Structure predict(const Config &c, StDescriptors &std, const Structure &st) {
+ if(config.template get<bool>("NORM") && !std.normalised && kernel.get_label()!="Kern_Linear")
+ norm.normalise(std);
+ return M_Tadah_Base::predict(c,std,st);
+ }
+
+ StructureDB predict(Config &c, const StructureDB &stdb, DC_Base &dc) {
+ return M_Tadah_Base::predict(c,stdb,dc);
+ }
- StructureDB predict(Config &c, const StructureDB &stdb, DC_Base &dc) {
- return M_Tadah_Base::predict(c,stdb,dc);
+ Config get_param_file() {
+ Config c = config;
+ c.remove("ALPHA");
+ c.remove("BETA");
+ c.remove("DBFILE");
+ c.remove("FORCE");
+ c.remove("STRESS");
+ c.remove("VERBOSE");
+ c.add("VERBOSE", 0);
+
+ c.clear_internal_keys();
+ c.remove("MODEL");
+ c.add("MODEL", label);
+ c.add("MODEL", kernel.get_label());
+ int modelN=config.template get<int>("MODEL",2);
+ c.add("MODEL", modelN);
+
+ for (size_t i=0;i<weights.size();++i) {
+ c.add("WEIGHTS", weights(i));
}
- Config get_param_file() {
- Config c = config;
- c.remove("ALPHA");
- c.remove("BETA");
- c.remove("DBFILE");
- c.remove("FORCE");
- c.remove("STRESS");
- c.remove("VERBOSE");
- c.add("VERBOSE", 0);
-
- c.clear_internal_keys();
- c.remove("MODEL");
- c.add("MODEL", label);
- c.add("MODEL", kernel.get_label());
- int modelN=config.template get<int>("MODEL",2);
- c.add("MODEL", modelN);
-
- for (size_t i=0;i<weights.size();++i) {
- c.add("WEIGHTS", weights(i));
+ if(config.template get<bool>("NORM")) {
+ for (size_t i=0;i<norm.mean.size();++i) {
+ c.add("NMEAN", norm.mean[i]);
}
-
- if(config.template get<bool>("NORM")) {
- for (size_t i=0;i<norm.mean.size();++i) {
- c.add("NMEAN", norm.mean[i]);
- }
- for (size_t i=0;i<norm.std_dev.size();++i) {
- c.add("NSTDEV", norm.std_dev[i]);
- }
+ for (size_t i=0;i<norm.std_dev.size();++i) {
+ c.add("NSTDEV", norm.std_dev[i]);
}
- if (kernel.get_label()!="Kern_Linear") {
- // dump basis to the config file file
- // make sure keys are not accidently assigned
- if (c.exist("SBASIS"))
- c.remove("SBASIS");
- if (c.exist("BASIS"))
- c.remove("BASIS");
- c.add("SBASIS", basis.b.cols());
- for (size_t i=0;i<basis.b.cols();++i) {
- for (size_t j=0;j<basis.b.rows();++j) {
- c.add("BASIS", basis.b(j,i));
- }
+ }
+ if (kernel.get_label()!="Kern_Linear") {
+ // dump basis to the config file file
+ // make sure keys are not accidently assigned
+ if (c.exist("SBASIS"))
+ c.remove("SBASIS");
+ if (c.exist("BASIS"))
+ c.remove("BASIS");
+ c.add("SBASIS", basis.b.cols());
+ for (size_t i=0;i<basis.b.cols();++i) {
+ for (size_t j=0;j<basis.b.rows();++j) {
+ c.add("BASIS", basis.b(j,i));
}
}
- return c;
}
- StructureDB predict(Config config_pred, StructureDB &stdb, DC_Base &dc,
- aed_type2 &predicted_error) {
-
- LinearRegressor::read_sigma(config_pred,Sigma);
- DesignMatrix<K> dm(kernel,config_pred);
- dm.scale=false; // do not scale energy, forces and stresses
- dm.build(stdb,norm,dc);
-
- // compute error
- predicted_error = T_MDMT_diag(dm.Phi, Sigma);
- double pmean = sqrt(predicted_error.mean());
-
- // compute energy, forces and stresses
- aed_type2 Tpred = T_dgemv(dm.Phi, weights);
-
- // Construct StructureDB object with predicted values
- StructureDB stdb_;
- stdb_.structures.resize(stdb.size());
- size_t i=0;
- for (size_t s=0; s<stdb.size(); ++s) {
- stdb_(s) = Structure(stdb(s));
-
- stdb_(s).energy = Tpred(i++);
- predicted_error(i) = (sqrt(predicted_error(i))-pmean)/stdb(s).natoms();
- if (config_pred.get<bool>("FORCE")) {
- for (size_t a=0; a<stdb(s).natoms(); ++a) {
- for (size_t k=0; k<3; ++k) {
- stdb_(s).atoms[a].force[k] = Tpred(i++);
- predicted_error(i) = (sqrt(predicted_error(i))-pmean);
- }
+ return c;
+ }
+ StructureDB predict(Config config_pred, StructureDB &stdb, DC_Base &dc,
+ aed_type2 &predicted_error) {
+
+ LinearRegressor::read_sigma(config_pred,Sigma);
+ DesignMatrix<K> dm(kernel,config_pred);
+ dm.scale=false; // do not scale energy, forces and stresses
+ dm.build(stdb,norm,dc);
+
+ // compute error
+ predicted_error = T_MDMT_diag(dm.Phi, Sigma);
+ double pmean = sqrt(predicted_error.mean());
+
+ // compute energy, forces and stresses
+ aed_type2 Tpred = T_dgemv(dm.Phi, weights);
+
+ // Construct StructureDB object with predicted values
+ StructureDB stdb_;
+ stdb_.structures.resize(stdb.size());
+ size_t i=0;
+ for (size_t s=0; s<stdb.size(); ++s) {
+ stdb_(s) = Structure(stdb(s));
+
+ stdb_(s).energy = Tpred(i++);
+ predicted_error(i) = (sqrt(predicted_error(i))-pmean)/stdb(s).natoms();
+ if (config_pred.get<bool>("FORCE")) {
+ for (size_t a=0; a<stdb(s).natoms(); ++a) {
+ for (size_t k=0; k<3; ++k) {
+ stdb_(s).atoms[a].force[k] = Tpred(i++);
+ predicted_error(i) = (sqrt(predicted_error(i))-pmean);
}
}
- if (config_pred.get<bool>("STRESS")) {
- for (size_t x=0; x<3; ++x) {
- for (size_t y=x; y<3; ++y) {
- predicted_error(i) = (sqrt(predicted_error(i))-pmean);
- stdb_(s).stress(x,y) = Tpred(i++);
- if (x!=y)
- stdb_(s).stress(y,x) = stdb_(s).stress(x,y);
- }
+ }
+ if (config_pred.get<bool>("STRESS")) {
+ for (size_t x=0; x<3; ++x) {
+ for (size_t y=x; y<3; ++y) {
+ predicted_error(i) = (sqrt(predicted_error(i))-pmean);
+ stdb_(s).stress(x,y) = Tpred(i++);
+ if (x!=y)
+ stdb_(s).stress(y,x) = stdb_(s).stress(x,y);
}
}
}
- return stdb_;
}
- StructureDB predict(StructureDB &stdb) {
- if(!trained) throw std::runtime_error("This object is not trained!\n\
- Hint: check different predict() methods.");
-
- phi_type &Phi = desmat.Phi;
-
- // compute energy, forces and stresses
- aed_type2 Tpred = T_dgemv(Phi, weights);
-
- double eweightglob=config.template get<double>("EWEIGHT");
- double fweightglob=config.template get<double>("FWEIGHT");
- double sweightglob=config.template get<double>("SWEIGHT");
-
- // Construct StructureDB object with predicted values
- StructureDB stdb_;
- stdb_.structures.resize(stdb.size());
- size_t s=0;
- size_t i=0;
- while (i<Phi.rows()) {
-
- stdb_(s).energy = Tpred(i++)*stdb(s).natoms()/eweightglob/stdb(s).eweight;
- if (config.template get<bool>("FORCE")) {
- stdb_(s).atoms.resize(stdb(s).natoms());
- for (size_t a=0; a<stdb(s).natoms(); ++a) {
- for (size_t k=0; k<3; ++k) {
- stdb_(s).atoms[a].force[k] = Tpred(i++)/fweightglob/stdb(s).fweight;
- }
+ return stdb_;
+ }
+ StructureDB predict(StructureDB &stdb) {
+ if(!trained) throw std::runtime_error("This object is not trained!\n\
+Hint: check different predict() methods.");
+
+ phi_type &Phi = desmat.Phi;
+
+ // compute energy, forces and stresses
+ aed_type2 Tpred = T_dgemv(Phi, weights);
+
+ double eweightglob=config.template get<double>("EWEIGHT");
+ double fweightglob=config.template get<double>("FWEIGHT");
+ double sweightglob=config.template get<double>("SWEIGHT");
+
+ // Construct StructureDB object with predicted values
+ StructureDB stdb_;
+ stdb_.structures.resize(stdb.size());
+ size_t s=0;
+ size_t i=0;
+ while (i<Phi.rows()) {
+
+ stdb_(s).energy = Tpred(i++)*stdb(s).natoms()/eweightglob/stdb(s).eweight;
+ if (config.template get<bool>("FORCE")) {
+ stdb_(s).atoms.resize(stdb(s).natoms());
+ for (size_t a=0; a<stdb(s).natoms(); ++a) {
+ for (size_t k=0; k<3; ++k) {
+ stdb_(s).atoms[a].force[k] = Tpred(i++)/fweightglob/stdb(s).fweight;
}
}
- if (config.template get<bool>("STRESS")) {
- for (size_t x=0; x<3; ++x) {
- for (size_t y=x; y<3; ++y) {
- stdb_(s).stress(x,y) = Tpred(i++)/sweightglob/stdb(s).sweight;
- if (x!=y)
- stdb_(s).stress(y,x) = stdb_(s).stress(x,y);
- }
+ }
+ if (config.template get<bool>("STRESS")) {
+ for (size_t x=0; x<3; ++x) {
+ for (size_t y=x; y<3; ++y) {
+ stdb_(s).stress(x,y) = Tpred(i++)/sweightglob/stdb(s).sweight;
+ if (x!=y)
+ stdb_(s).stress(y,x) = stdb_(s).stress(x,y);
}
}
- s++;
}
- return stdb_;
+ s++;
}
+ return stdb_;
+ }
- private:
- std::string label="M_KRR";
- Basis<K> basis;
- DesignMatrix<K> desmat;
+private:
+ std::string label="M_KRR";
+ Basis<K> basis;
+ DesignMatrix<K> desmat;
- t_type convert_to_nweights(const t_type &weights) const {
- if(kernel.get_label()!="Kern_Linear") {
- throw std::runtime_error("Cannot convert weights to nweights for\n\
- non linear kernel\n");
- }
- t_type kw(weights.rows());
- if(config.template get<bool>("NORM") && kernel.get_label()=="Kern_Linear") {
- // normalise weights such that when predict is called
- // we can supply it with a non-normalised descriptor
- kw.resize(weights.rows());
- kw(0) = weights(0);
- for (size_t i=1; i<weights.size(); ++i) {
+ t_type convert_to_nweights(const t_type &weights) const {
+ if(kernel.get_label()!="Kern_Linear") {
+ throw std::runtime_error("Cannot convert weights to nweights for\n\
+non linear kernel\n");
+ }
+ t_type kw(weights.rows());
+ if(config.template get<bool>("NORM") && kernel.get_label()=="Kern_Linear") {
+ // normalise weights such that when predict is called
+ // we can supply it with a non-normalised descriptor
+ kw.resize(weights.rows());
+ kw(0) = weights(0);
+ for (size_t i=1; i<weights.size(); ++i) {
- if (norm.std_dev[i] > std::numeric_limits<double>::min())
- kw(i) = weights(i) / norm.std_dev[i];
- else
- kw(i) = weights(i);
+ if (norm.std_dev[i] > std::numeric_limits<double>::min())
+ kw(i) = weights(i) / norm.std_dev[i];
+ else
+ kw(i) = weights(i);
- kw(0) -= norm.mean[i]*kw(i);
+ kw(0) -= norm.mean[i]*kw(i);
- }
}
- return kw;
}
- // The opposite of convert_to_nweights()
- t_type convert_to_weights(const t_type &kw) const {
- if(kernel.get_label()!="Kern_Linear") {
- throw std::runtime_error("Cannot convert nweights to weights for\n\
- non linear kernel\n");
- }
- // convert normalised weights back to "normal"
- t_type w(kw.rows());
- w(0) = kw(0);
- for (size_t i=1; i<kw.size(); ++i) {
- if (norm.std_dev[i] > std::numeric_limits<double>::min())
- w(i) = kw(i) * norm.std_dev[i];
- else
- w(i) = kw(i);
+ return kw;
+ }
+ // The opposite of convert_to_nweights()
+ t_type convert_to_weights(const t_type &kw) const {
+ if(kernel.get_label()!="Kern_Linear") {
+ throw std::runtime_error("Cannot convert nweights to weights for\n\
+non linear kernel\n");
+ }
+ // convert normalised weights back to "normal"
+ t_type w(kw.rows());
+ w(0) = kw(0);
+ for (size_t i=1; i<kw.size(); ++i) {
+ if (norm.std_dev[i] > std::numeric_limits<double>::min())
+ w(i) = kw(i) * norm.std_dev[i];
+ else
+ w(i) = kw(i);
- w(0) += kw(i)*norm.mean[i];
- }
- return w;
+ w(0) += kw(i)*norm.mean[i];
}
+ return w;
+ }
- template <typename D>
- void train(D &desmat) {
- // TODO see comments in M_BLR
- phi_type Phi = desmat.Phi;
- t_type T = desmat.T;
- M_KRR_Train<K>::train(Phi,T);
+ template <typename D>
+ void train(D &desmat) {
+ // TODO see comments in M_BLR
+ phi_type Phi = desmat.Phi;
+ t_type T = desmat.T;
+ M_KRR_Train<K>::train(Phi,T);
- if (config.template get<bool>("NORM") &&
- kernel.get_label()=="Kern_Linear") {
- weights = convert_to_nweights(weights);
- }
- }
+ if (config.template get<bool>("NORM") &&
+ kernel.get_label()=="Kern_Linear") {
+ weights = convert_to_nweights(weights);
+ }
+ }
- // Do we want to confuse user with those and make them public?
- // Either way they must be stated as below to silence clang warning
- using M_KRR_Train<K>::predict;
- using M_KRR_Train<K>::train;
- using M_KRR_Train<K>::trained;
- using M_KRR_Train<K>::weights;
- using M_KRR_Train<K>::Sigma;
- using M_KRR_Train<K>::config;
- using M_KRR_Train<K>::kernel;
- using M_KRR_Train<K>::ekm;
+ // Do we want to confuse user with those and make them public?
+ // Either way they must be stated as below to silence clang warning
+ using M_KRR_Train<K>::predict;
+ using M_KRR_Train<K>::train;
+ using M_KRR_Train<K>::trained;
+ using M_KRR_Train<K>::weights;
+ using M_KRR_Train<K>::Sigma;
+ using M_KRR_Train<K>::config;
+ using M_KRR_Train<K>::kernel;
+ using M_KRR_Train<K>::ekm;
};
#endif
diff --git a/src/dm_bf_linear.cpp b/src/dm_bf_linear.cpp
index f04f66e4921900f8ec3dea20e5c5adae8bc3fd1d..ab467e8c4ba7eb31f43ab3cdc26ec860061b312a 100644
--- a/src/dm_bf_linear.cpp
+++ b/src/dm_bf_linear.cpp
@@ -44,6 +44,7 @@ void DM_BF_Linear::calc_phi_force_rows(phi_type &Phi, size_t &row,
void DM_BF_Linear::calc_phi_stress_rows(phi_type &Phi, size_t &row,
const double fac[6], const Structure &st, const StDescriptors &st_d)
{
+ double V_inv = 1/st.get_volume();
for (size_t i=0; i<st.natoms(); ++i) {
const Vec3d &ri = st(i).position;
for (size_t jj=0; jj<st_d.fd[i].size(); ++jj) {
@@ -56,7 +57,7 @@ void DM_BF_Linear::calc_phi_stress_rows(phi_type &Phi, size_t &row,
for (size_t x=0; x<3; ++x) {
for (size_t y=x; y<3; ++y) {
for (size_t d=0; d<fdij.rows(); ++d) {
- Phi(row+mn,d) += (fdij(d,y)-fdji(d,y))*0.5*fac[mn]*(ri(x)-rj(x));
+ Phi(row+mn,d) += V_inv*(fdij(d,y)-fdji(d,y))*0.5*fac[mn]*(ri(x)-rj(x));
}
mn++;
}
diff --git a/src/dm_bf_polynomial2.cpp b/src/dm_bf_polynomial2.cpp
index 51c3f561ce32df7c67017f379ec9f76a5eb3f123..6638fdb19ce46998b511e14e90b8b0b0f2f51c0f 100644
--- a/src/dm_bf_polynomial2.cpp
+++ b/src/dm_bf_polynomial2.cpp
@@ -64,6 +64,7 @@ void DM_BF_Polynomial2::calc_phi_stress_rows(phi_type &Phi,
const Structure &st,
const StDescriptors &st_d)
{
+ double V_inv = 1/st.get_volume();
for (size_t a=0; a<st.natoms(); ++a) {
const Vec3d &ri = st(a).position;
const aed_type2& aedi = st_d.get_aed(a);
@@ -80,7 +81,7 @@ void DM_BF_Polynomial2::calc_phi_stress_rows(phi_type &Phi,
size_t b=0;
for (size_t i=0; i<fdij.rows(); ++i) {
for (size_t ii=i; ii<fdij.rows(); ++ii) {
- Phi(row+mn,b) += 0.5*fac[mn]*(ri(x)-rj(x))
+ Phi(row+mn,b) += V_inv*0.5*fac[mn]*(ri(x)-rj(x))
*(fdij(i,y)*aedi(ii) + fdij(ii,y)*aedi(i)
- fdji(i,y)*aedj(ii) - fdji(ii,y)*aedj(i));
b++;
diff --git a/src/dm_kern_base.cpp b/src/dm_kern_base.cpp
index 92d62af48cb6227225fd049ada42970b2bfed683..ecf679c39ad18a36219dc36aa0aa3de34cb80c95 100644
--- a/src/dm_kern_base.cpp
+++ b/src/dm_kern_base.cpp
@@ -3,64 +3,65 @@
DM_Kern_Base::~DM_Kern_Base() {}
size_t DM_Kern_Base::get_phi_cols(const Config &)
{
- return basis.cols();
+ return basis.cols();
}
void DM_Kern_Base::calc_phi_energy_row(phi_type &Phi, size_t &row, const double fac,
- const Structure &, const StDescriptors &st_d)
+ const Structure &, const StDescriptors &st_d)
{
- for (size_t a=0; a<st_d.naed();++a) {
- for (size_t b=0; b<basis.cols(); ++b) {
- Phi(row,b) += (*this)(basis.col(b),st_d.get_aed(a))*fac;
- }
+ for (size_t a=0; a<st_d.naed();++a) {
+ for (size_t b=0; b<basis.cols(); ++b) {
+ Phi(row,b) += (*this)(basis.col(b),st_d.get_aed(a))*fac;
}
- row++;
- //Phi.row(row++) *= fac;
+ }
+ row++;
+ //Phi.row(row++) *= fac;
}
void DM_Kern_Base::calc_phi_force_rows(phi_type &Phi, size_t &row, const double fac,
- const Structure &st, const StDescriptors &st_d) {
- for (size_t i=0; i<st.natoms(); ++i) {
- const aed_type2& aedi = st_d.get_aed(i);
- for (size_t jj=0; jj<st_d.fd[i].size(); ++jj) {
- size_t j=st.near_neigh_idx[i][jj];
- size_t ii = st.get_nn_iindex(i,j,jj);
- const fd_type &fdji = st_d.fd[j][ii];
- const fd_type &fdij = st_d.fd[i][jj];
- const aed_type2& aedj = st_d.get_aed(j);
- for (size_t b=0; b<basis.cols(); ++b) {
- for (size_t k=0; k<3; ++k) {
- Phi(row+k,b) -= fac*((*this).prime(basis.col(b), aedi,fdij(k)) -
- (*this).prime(basis.col(b),aedj,fdji(k)));
- }
- }
+ const Structure &st, const StDescriptors &st_d) {
+ for (size_t i=0; i<st.natoms(); ++i) {
+ const aed_type2& aedi = st_d.get_aed(i);
+ for (size_t jj=0; jj<st_d.fd[i].size(); ++jj) {
+ size_t j=st.near_neigh_idx[i][jj];
+ size_t ii = st.get_nn_iindex(i,j,jj);
+ const fd_type &fdji = st_d.fd[j][ii];
+ const fd_type &fdij = st_d.fd[i][jj];
+ const aed_type2& aedj = st_d.get_aed(j);
+ for (size_t b=0; b<basis.cols(); ++b) {
+ for (size_t k=0; k<3; ++k) {
+ Phi(row+k,b) -= fac*((*this).prime(basis.col(b), aedi,fdij(k)) -
+ (*this).prime(basis.col(b),aedj,fdji(k)));
}
- row+=3;
+ }
}
+ row+=3;
+ }
}
void DM_Kern_Base::calc_phi_stress_rows(phi_type &Phi, size_t &row, const double fac[6],
- const Structure &st, const StDescriptors &st_d)
+ const Structure &st, const StDescriptors &st_d)
{
- for (size_t i=0; i<st.natoms(); ++i) {
- const Vec3d &ri = st(i).position;
- const aed_type2& aedi = st_d.get_aed(i);
- for (size_t jj=0; jj<st_d.fd[i].size(); ++jj) {
- size_t j=st.near_neigh_idx[i][jj];
- size_t ii = st.get_nn_iindex(i,j,jj);
- const fd_type &fdji = st_d.fd[j][ii];
- const fd_type &fdij = st_d.fd[i][jj];
- const aed_type2& aedj = st_d.get_aed(j);
- const Vec3d &rj = st.nn_pos(i,jj);
- size_t mn=0;
- for (size_t x=0; x<3; ++x) {
- for (size_t y=x; y<3; ++y) {
- for (size_t b=0; b<basis.cols(); ++b) {
- Phi(row+mn,b) += 0.5*fac[mn]*(ri(x)-rj(x))*
- ((*this).prime(basis.col(b),aedi,fdij(y)) -
- (*this).prime(basis.col(b),aedj,fdji(y)));
- }
- mn++;
- }
- }
+ double V_inv = 1/st.get_volume();
+ for (size_t i=0; i<st.natoms(); ++i) {
+ const Vec3d &ri = st(i).position;
+ const aed_type2& aedi = st_d.get_aed(i);
+ for (size_t jj=0; jj<st_d.fd[i].size(); ++jj) {
+ size_t j=st.near_neigh_idx[i][jj];
+ size_t ii = st.get_nn_iindex(i,j,jj);
+ const fd_type &fdji = st_d.fd[j][ii];
+ const fd_type &fdij = st_d.fd[i][jj];
+ const aed_type2& aedj = st_d.get_aed(j);
+ const Vec3d &rj = st.nn_pos(i,jj);
+ size_t mn=0;
+ for (size_t x=0; x<3; ++x) {
+ for (size_t y=x; y<3; ++y) {
+ for (size_t b=0; b<basis.cols(); ++b) {
+ Phi(row+mn,b) += V_inv*0.5*fac[mn]*(ri(x)-rj(x))*
+ ((*this).prime(basis.col(b),aedi,fdij(y)) -
+ (*this).prime(basis.col(b),aedj,fdji(y)));
+ }
+ mn++;
}
+ }
}
- row+=6;
+ }
+ row+=6;
}
diff --git a/src/dm_kern_linear.cpp b/src/dm_kern_linear.cpp
index c04a361d072d474389ef8b0bcbdee9cc25340ba2..90f97a7eda24374a7bad06a72f926a97c67370fe 100644
--- a/src/dm_kern_linear.cpp
+++ b/src/dm_kern_linear.cpp
@@ -8,61 +8,62 @@ DM_Kern_Linear::DM_Kern_Linear (const Config &c): Kern_Linear(c)
{}
size_t DM_Kern_Linear::get_phi_cols(const Config &config)
{
- size_t cols = config.get<size_t>("DSIZE");
- return cols;
+ size_t cols = config.get<size_t>("DSIZE");
+ return cols;
}
void DM_Kern_Linear::calc_phi_energy_row(phi_type &Phi, size_t &row,
- const double fac, const Structure &, const StDescriptors &st_d)
+ const double fac, const Structure &, const StDescriptors &st_d)
{
- for (size_t a=0; a<st_d.naed();++a) {
- const aed_type2 &aed = st_d.get_aed(a); // TODO
- for (size_t j=0; j<aed.size(); ++j) {
- Phi(row,j)+=aed[j]*fac;
- }
+ for (size_t a=0; a<st_d.naed();++a) {
+ const aed_type2 &aed = st_d.get_aed(a); // TODO
+ for (size_t j=0; j<aed.size(); ++j) {
+ Phi(row,j)+=aed[j]*fac;
}
- row++;
+ }
+ row++;
}
void DM_Kern_Linear::calc_phi_force_rows(phi_type &Phi, size_t &row,
- const double fac, const Structure &st, const StDescriptors &st_d)
+ const double fac, const Structure &st, const StDescriptors &st_d)
{
- for (size_t a=0; a<st.natoms(); ++a) {
- for (size_t jj=0; jj<st_d.fd[a].size(); ++jj) {
- const size_t j=st.near_neigh_idx[a][jj];
- const size_t aa = st.get_nn_iindex(a,j,jj);
- for (size_t k=0; k<3; ++k) {
- aed_type2 temp = (st_d.fd[a][jj](k)-
- st_d.fd[j][aa](k))*fac;
- for (size_t d=0; d<temp.size(); ++d) {
- Phi(row+k,d) -= temp[d];
- }
- }
+ for (size_t a=0; a<st.natoms(); ++a) {
+ for (size_t jj=0; jj<st_d.fd[a].size(); ++jj) {
+ const size_t j=st.near_neigh_idx[a][jj];
+ const size_t aa = st.get_nn_iindex(a,j,jj);
+ for (size_t k=0; k<3; ++k) {
+ aed_type2 temp = (st_d.fd[a][jj](k)-
+ st_d.fd[j][aa](k))*fac;
+ for (size_t d=0; d<temp.size(); ++d) {
+ Phi(row+k,d) -= temp[d];
}
- row+=3;
+ }
}
+ row+=3;
+ }
}
void DM_Kern_Linear::calc_phi_stress_rows(phi_type &Phi, size_t &row,
- const double fac[6], const Structure &st, const StDescriptors &st_d)
+ const double fac[6], const Structure &st, const StDescriptors &st_d)
{
- for (size_t i=0; i<st.natoms(); ++i) {
- const Vec3d &ri = st(i).position;
- for (size_t jj=0; jj<st_d.fd[i].size(); ++jj) {
- const size_t j=st.near_neigh_idx[i][jj];
- const size_t ii = st.get_nn_iindex(i,j,jj);
- const fd_type &fdij = st_d.fd[i][jj];
- const fd_type &fdji = st_d.fd[j][ii];
- const Vec3d &rj = st.nn_pos(i,jj);
- size_t mn=0;
- for (size_t x=0; x<3; ++x) {
- for (size_t y=x; y<3; ++y) {
- aed_type2 temp = (fdij(y)-fdji(y))*0.5*fac[mn]*(ri(x)-rj(x));
- for (size_t d=0; d<temp.size(); ++d) {
- Phi(row+mn,d) += temp[d];
- }
- mn++;
- }
- }
+ double V_inv = 1/st.get_volume();
+ for (size_t i=0; i<st.natoms(); ++i) {
+ const Vec3d &ri = st(i).position;
+ for (size_t jj=0; jj<st_d.fd[i].size(); ++jj) {
+ const size_t j=st.near_neigh_idx[i][jj];
+ const size_t ii = st.get_nn_iindex(i,j,jj);
+ const fd_type &fdij = st_d.fd[i][jj];
+ const fd_type &fdji = st_d.fd[j][ii];
+ const Vec3d &rj = st.nn_pos(i,jj);
+ size_t mn=0;
+ for (size_t x=0; x<3; ++x) {
+ for (size_t y=x; y<3; ++y) {
+ aed_type2 temp = V_inv*(fdij(y)-fdji(y))*0.5*fac[mn]*(ri(x)-rj(x));
+ for (size_t d=0; d<temp.size(); ++d) {
+ Phi(row+mn,d) += temp[d];
+ }
+ mn++;
}
+ }
}
- row += 6;
+ }
+ row += 6;
}
diff --git a/src/dm_kern_lq.cpp b/src/dm_kern_lq.cpp
index 0333db06a7ce4a9c9113e4b484bc1dc007d43db4..3c14dfaad5465facd3ef725b13d7a323fd29f7b0 100644
--- a/src/dm_kern_lq.cpp
+++ b/src/dm_kern_lq.cpp
@@ -4,8 +4,8 @@
//CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_Kern_LQ> DM_Kern_LQ_2( "Kern_LQ" );
DM_Kern_LQ::DM_Kern_LQ():
- Kern_LQ()
+ Kern_LQ()
{}
DM_Kern_LQ::DM_Kern_LQ(const Config &c):
- Kern_LQ(c)
+ Kern_LQ(c)
{}
diff --git a/src/dm_kern_polynomial.cpp b/src/dm_kern_polynomial.cpp
index b5fccbd4ab82ae3c397a002abcdadddeacb582c8..b60e9d913d808a63f48a6bb76cbb2547888b3500 100644
--- a/src/dm_kern_polynomial.cpp
+++ b/src/dm_kern_polynomial.cpp
@@ -4,8 +4,8 @@
//CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_Kern_Polynomial> DM_Kern_Polynomial_2( "Kern_Polynomial" );
DM_Kern_Polynomial::DM_Kern_Polynomial():
- Kern_Polynomial()
+ Kern_Polynomial()
{}
DM_Kern_Polynomial::DM_Kern_Polynomial(const Config &c):
- Kern_Polynomial(c)
+ Kern_Polynomial(c)
{}
diff --git a/src/dm_kern_quadratic.cpp b/src/dm_kern_quadratic.cpp
index 402c7771f57933a289a08802eddc6de1ad79ddfb..83a71981b1312f04f7c713d21f1e4ac546402d08 100644
--- a/src/dm_kern_quadratic.cpp
+++ b/src/dm_kern_quadratic.cpp
@@ -4,8 +4,8 @@
//CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_Kern_Quadratic> DM_Kern_Quadratic_2( "Kern_Quadratic" );
DM_Kern_Quadratic::DM_Kern_Quadratic():
- Kern_Quadratic()
+ Kern_Quadratic()
{}
DM_Kern_Quadratic::DM_Kern_Quadratic(const Config &c):
- Kern_Quadratic(c)
+ Kern_Quadratic(c)
{}
diff --git a/src/dm_kern_rbf.cpp b/src/dm_kern_rbf.cpp
index 7ff05ef7af2d432ef6a82e33d1118c1d1ae31374..1870190d591d2c0df7b88acadbcd00f9b5839326 100644
--- a/src/dm_kern_rbf.cpp
+++ b/src/dm_kern_rbf.cpp
@@ -4,8 +4,8 @@
//CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_Kern_RBF> DM_Kern_RBF_2( "Kern_RBF" );
DM_Kern_RBF::DM_Kern_RBF():
- Kern_RBF()
+ Kern_RBF()
{}
DM_Kern_RBF::DM_Kern_RBF(const Config &c):
- Kern_RBF(c)
+ Kern_RBF(c)
{}
diff --git a/src/dm_kern_sigmoid.cpp b/src/dm_kern_sigmoid.cpp
index d1fb598553daa2c9fb695f2da49c3393e246ed6e..2da37c905a1be2f801beafd7ed9bbcbd83f99295 100644
--- a/src/dm_kern_sigmoid.cpp
+++ b/src/dm_kern_sigmoid.cpp
@@ -4,8 +4,8 @@
//CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_Kern_Sigmoid> DM_Kern_Sigmoid_2( "Kern_Sigmoid" );
DM_Kern_Sigmoid::DM_Kern_Sigmoid():
- Kern_Sigmoid()
+ Kern_Sigmoid()
{}
DM_Kern_Sigmoid::DM_Kern_Sigmoid(const Config &c):
- Kern_Sigmoid(c)
+ Kern_Sigmoid(c)
{}
diff --git a/src/m_tadah_base.cpp b/src/m_tadah_base.cpp
index 11cac065b7dff94a3c0eb954d44c95b9c88ab74c..d78f49839c6cea61e2c57595059d130f0bfb280a 100644
--- a/src/m_tadah_base.cpp
+++ b/src/m_tadah_base.cpp
@@ -105,6 +105,7 @@ void M_Tadah_Base::
spredict(const size_t a, stress_type &s,
const StDescriptors &std, const Structure &st)
{
+ double V_inv = 1/st.get_volume();
const Vec3d &ri = st.atoms[a].position;
const aed_type2 &aedi = std.get_aed(a);
for (size_t jj=0; jj<st.nn_size(a); ++jj) {
@@ -122,7 +123,7 @@ spredict(const size_t a, stress_type &s,
// resutling from the pairwise i-j interaction.
for (size_t x=0; x<3; ++x) {
for (size_t y=x; y<3; ++y) {
- s(x,y) -= 0.5*(ri[x]-rj[x])*(fij[y]-fji[y]);
+ s(x,y) -= V_inv*0.5*(ri[x]-rj[x])*(fij[y]-fji[y]);
}
}
}
diff --git a/tests/test_structure.cpp b/tests/test_structure.cpp
index 1a0292c725f931e458256a6f5717bd6b58ead8a5..e1b3446f0f29818e78b1554a3c49f1dc9195e460 100644
--- a/tests/test_structure.cpp
+++ b/tests/test_structure.cpp
@@ -1,5 +1,4 @@
#include <string>
-#include <filesystem>
#include "catch2/catch.hpp"
#include <tadah/mlip/atom.h>
#include <tadah/mlip/structure.h>