improving the energy computation model in SPIN/compute_spin.cpp

This commit is contained in:
Julien Tranchida
2022-09-21 10:06:41 +02:00
parent 7e8e40fefa
commit 1ef2c8c5dc
2 changed files with 44 additions and 6 deletions

View File

@ -43,7 +43,8 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */
ComputeSpin::ComputeSpin(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), pair(nullptr), spin_pairs(nullptr), lockprecessionspin(nullptr)
Compute(lmp, narg, arg), pair(nullptr), spin_pairs(nullptr),
lockprecessionspin(nullptr)
{
if ((narg != 3) && (narg != 4)) error->all(FLERR,"Illegal compute compute/spin command");
@ -136,15 +137,49 @@ void ComputeSpin::init()
}
}
// ptrs FixPrecessionSpin classes
// set ptrs for fix precession/spin styles
// loop 1: obtain # of fix precession/spin styles
int iforce;
nprecspin = 0;
for (iforce = 0; iforce < modify->nfix; iforce++) {
if (utils::strmatch(modify->fix[iforce]->style,"^precession/spin")) {
precession_spin_flag = 1;
lockprecessionspin = dynamic_cast<FixPrecessionSpin *>(modify->fix[iforce]);
nprecspin++;
}
}
// init length of vector of ptrs to precession/spin styles
if (nprecspin > 0) {
lockprecessionspin = new FixPrecessionSpin*[nprecspin];
}
// loop 2: fill vector with ptrs to precession/spin styles
int count2 = 0;
if (nprecspin > 0) {
for (iforce = 0; iforce < modify->nfix; iforce++) {
if (utils::strmatch(modify->fix[iforce]->style,"^precession/spin")) {
precession_spin_flag = 1;
lockprecessionspin[count2] = dynamic_cast<FixPrecessionSpin *>(modify->fix[iforce]);
count2++;
}
}
}
if (count2 != nprecspin)
error->all(FLERR,"Incorrect number of precession/spin fixes");
// // ptrs FixPrecessionSpin classes
// int iforce;
// for (iforce = 0; iforce < modify->nfix; iforce++) {
// if (utils::strmatch(modify->fix[iforce]->style,"^precession/spin")) {
// precession_spin_flag = 1;
// lockprecessionspin = dynamic_cast<FixPrecessionSpin *>(modify->fix[iforce]);
// }
// }
}
/* ---------------------------------------------------------------------- */
@ -192,7 +227,9 @@ void ComputeSpin::compute_vector()
// update magnetic precession energies
if (precession_spin_flag) {
magenergy += lockprecessionspin->emag[i];
for (int k = 0; k < nprecspin; k++) {
magenergy += lockprecessionspin[k]->emag[i];
}
}
// update magnetic pair interactions

View File

@ -40,7 +40,8 @@ class ComputeSpin : public Compute {
// pointers to magnetic fixes
class FixPrecessionSpin *lockprecessionspin;
int nprecspin;
class FixPrecessionSpin **lockprecessionspin;
// pointers to magnetic pair styles