modifications in fix adapt.cpp

This commit is contained in:
Evangelos Voyiatzis
2025-04-07 11:58:24 +02:00
committed by GitHub
parent 2b38b6d1c0
commit 239eb64cbb

View File

@ -19,6 +19,8 @@
#include "atom.h"
#include "bond.h"
#include "bond_hybrid.h"
#include "dihedral.h"
#include "dihedral_hybrid.h"
#include "domain.h"
#include "error.h"
#include "fix_store_atom.h"
@ -43,14 +45,14 @@ using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
enum{PAIR, KSPACE, ATOM, BOND, ANGLE, IMPROPER};
enum{PAIR, KSPACE, ATOM, BOND, ANGLE, DIHEDRAL, IMPROPER};
enum{DIAMETER, CHARGE};
/* ---------------------------------------------------------------------- */
FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), nadapt(0), anypair(0), anybond(0), anyangle(0),
anyimproper(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr)
anydihedral(0), anyimproper(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr)
{
if (narg < 5) utils::missing_cmd_args(FLERR,"fix adapt", error);
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
@ -85,6 +87,10 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) :
if (iarg+5 > narg) utils::missing_cmd_args(FLERR,"fix adapt angle", error);
nadapt++;
iarg += 5;
} else if (strcmp(arg[iarg],"dihedral") == 0) {
if (iarg+5 > narg) utils::missing_cmd_args(FLERR,"fix adapt dihedral", error);
nadapt++;
iarg += 5;
} else if (strcmp(arg[iarg],"improper") == 0) {
if (iarg+5 > narg) utils::missing_cmd_args(FLERR,"fix adapt improper", error);
nadapt++;
@ -155,6 +161,19 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) :
nadapt++;
iarg += 5;
} else if (strcmp(arg[iarg],"dihedral") == 0) {
adapt[nadapt].which = DIHEDRAL;
adapt[nadapt].dihedral = nullptr;
adapt[nadapt].dstyle = utils::strdup(arg[iarg+1]);
adapt[nadapt].dparam = utils::strdup(arg[iarg+2]);
utils::bounds_typelabel(FLERR, arg[iarg+3], 1, atom->ndihedraltypes,
adapt[nadapt].ilo, adapt[nadapt].ihi, lmp, Atom::DIHEDRAL);
if (utils::strmatch(arg[iarg+4],"^v_")) {
adapt[nadapt].var = utils::strdup(arg[iarg+4]+2);
} else error->all(FLERR,"Argument #{} must be variable not {}", iarg+5, arg[iarg+4]);
nadapt++;
iarg += 5;
} else if (strcmp(arg[iarg],"improper") == 0) {
adapt[nadapt].which = IMPROPER;
adapt[nadapt].improper = nullptr;
@ -246,6 +265,12 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) :
// allocate improper style arrays:
n = atom->ndihedraltypes;
for (int m = 0; m < nadapt; ++m)
if (adapt[m].which == DIHEDRAL) memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig");
// allocate improper style arrays:
n = atom->nimpropertypes;
for (int m = 0; m < nadapt; ++m)
if (adapt[m].which == IMPROPER) memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig");
@ -270,6 +295,10 @@ FixAdapt::~FixAdapt()
delete[] adapt[m].astyle;
delete[] adapt[m].aparam;
memory->destroy(adapt[m].vector_orig);
} else if (adapt[m].which == DIHEDRAL) {
delete[] adapt[m].dstyle;
delete[] adapt[m].dparam;
memory->destroy(adapt[m].vector_orig);
} else if (adapt[m].which == IMPROPER) {
delete[] adapt[m].istyle;
delete[] adapt[m].iparam;
@ -368,6 +397,7 @@ void FixAdapt::init()
anypair = 0;
anybond = 0;
anyangle = 0;
anydihedral = 0;
anyimproper = 0;
for (int m = 0; m < nadapt; m++) {
@ -509,6 +539,40 @@ void FixAdapt::init()
delete[] astyle;
} else if (ad->which == DIHEDRAL) {
ad->dihedral = nullptr;
anydihedral = 1;
char *dstyle = utils::strdup(ad->dstyle);
if (lmp->suffix_enable)
ad->dihedral = force->dihedral_match(fmt::format("{}/{}",dstyle,lmp->suffix));
if (ad->dihedral == nullptr) ad->dihedral = force->dihedral_match(dstyle);
if (ad->dihedral == nullptr )
error->all(FLERR,"Fix adapt dihedral style {} does not exist", dstyle);
void *ptr = ad->dihedral->extract(ad->dparam,ad->ddim);
if (ptr == nullptr)
error->all(FLERR,"Fix adapt dihedral style parameter {} not supported", ad->dparam);
// for dihedral styles, use a vector
if (ad->ddim == 1) ad->vector = (double *) ptr;
if (utils::strmatch(force->dihedral_style,"^hybrid")) {
auto dihedral = dynamic_cast<DihedralHybrid *>(force->dihedral);
if (dihedral) {
for (i = ad->ilo; i <= ad->ihi; i++) {
if (!dihedral->check_itype(i,dstyle))
error->all(FLERR,"Fix adapt type dihedral range is not valid "
"for dihedral hybrid sub-style {}", dstyle);
}
}
}
delete[] dstyle;
} else if (ad->which == IMPROPER) {
ad->improper = nullptr;
anyimproper = 1;
@ -569,7 +633,7 @@ void FixAdapt::init()
if (restart_reset) restart_reset = 0;
// make copy of original pair/bond/angle/improper array values
// make copy of original pair/bond/angle/dihedral/improper array values
for (int m = 0; m < nadapt; m++) {
Adapt *ad = &adapt[m];
@ -588,6 +652,10 @@ void FixAdapt::init()
for (i = ad->ilo; i <= ad->ihi; ++i )
ad->vector_orig[i] = ad->vector[i];
} else if (ad->which == DIHEDRAL && ad->ddim == 1) {
for (i = ad->ilo; i <= ad->ihi; ++i )
ad->vector_orig[i] = ad->vector[i];
} else if (ad->which == IMPROPER && ad->idim == 1) {
for (i = ad->ilo; i <= ad->ihi; ++i )
ad->vector_orig[i] = ad->vector[i];
@ -710,6 +778,18 @@ void FixAdapt::change_settings()
ad->vector[i] = value;
}
// set dihedral type array values:
} else if (ad->which == DIHEDRAL) {
if (ad->ddim == 1) {
if (scaleflag)
for (i = ad->ilo; i <= ad->ihi; ++i )
ad->vector[i] = value*ad->vector_orig[i];
else
for (i = ad->ilo; i <= ad->ihi; ++i )
ad->vector[i] = value;
}
// set improper type array values:
} else if (ad->which == IMPROPER) {
@ -797,6 +877,7 @@ void FixAdapt::change_settings()
if (anypair) force->pair->reinit();
if (anybond) force->bond->reinit();
if (anyangle) force->angle->reinit();
if (anydihedral) force->dihedral->reinit();
if (anyimproper) force->improper->reinit();
// reset KSpace charges if charges have changed
@ -832,6 +913,12 @@ void FixAdapt::restore_settings()
ad->vector[i] = ad->vector_orig[i];
}
} else if (ad->which == DIHEDRAL) {
if (ad->ddim == 1) {
for (int i = ad->ilo; i <= ad->ihi; i++)
ad->vector[i] = ad->vector_orig[i];
}
} else if (ad->which == IMPROPER) {
if (ad->idim == 1) {
for (int i = ad->ilo; i <= ad->ihi; i++)
@ -878,6 +965,7 @@ void FixAdapt::restore_settings()
if (anypair) force->pair->reinit();
if (anybond) force->bond->reinit();
if (anyangle) force->angle->reinit();
if (anydihedral) force->dihedral->reinit();
if (anyimproper) force->improper->reinit();
if (chgflag && force->kspace) force->kspace->qsum_qsq();
}