MODIFY FIX EFIELD SOURCE CODE

Hello everyone , i want to add the region command to the fix efield i did modify the .h and .cpp files in src folder but i still can not use an electrostatic field in a specified region with at an angle needed for my simulations i hope you can help me and thanks

this is the .h code

#ifdef FIX_CLASS

FixStyle(efield,FixEfield)

#else

#ifndef LMP_FIX_EFIELD_H
#define LMP_FIX_EFIELD_H

#include "fix.h"

namespace LAMMPS_NS {

class FixEfield : public Fix {
public:
FixEfield(class LAMMPS *, int, char **);
~FixEfield();
int setmask();
void init();
void setup(int);
void post_force(int);
void post_force_respa(int, int, int);
double memory_usage();

private:
double ex,ey,ez;
int varflag,iregion;
char *xstr,*ystr,*zstr,*estr;
char *idregion;
int xvar,yvar,zvar,evar,xstyle,ystyle,zstyle,estyle;
double foriginal[4],foriginal_all[4];
int nlevels_respa;
int force_flag;
double qe2f;

int maxatom;
double **efield;
};

}

#endif
#endif

this is my .cpp code

#include
#include
#include
#include "fix_efield.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "domain.h"
#include "region.h"
#include "force.h"
#include "respa.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"

using namespace LAMMPS_NS;
using namespace FixConst;

enum{NONE,CONSTANT,EQUAL,ATOM};

/* ---------------------------------------------------------------------- */

FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg != 6) error->all(FLERR,"Illegal fix efield command");

qe2f = force->qe2f;
xstr = ystr = zstr = NULL;

if (strstr(arg[3],"v_") == arg[3]) {
int n = strlen(&arg[3][2]) + 1;
xstr = new char[n];
strcpy(xstr,&arg[3][2]);
} else {
ex = qe2f * atof(arg[3]);
xstyle = CONSTANT;
}

if (strstr(arg[4],"v_") == arg[4]) {
int n = strlen(&arg[4][2]) + 1;
ystr = new char[n];
strcpy(ystr,&arg[4][2]);
} else {
ey = qe2f * atof(arg[4]);
ystyle = CONSTANT;
}

if (strstr(arg[5],"v_") == arg[5]) {
int n = strlen(&arg[5][2]) + 1;
zstr = new char[n];
strcpy(zstr,&arg[5][2]);
} else {
ez = qe2f * atof(arg[5]);
zstyle = CONSTANT;
}

// optional args

iregion = -1;
idregion = NULL;
estr = NULL;

int iarg = 6;
while (iarg < narg) {
if (strcmp(arg[iarg],"region") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix efield command");
iregion = domain->find_region(arg[iarg+1]);
if (iregion == -1)
error->all(FLERR,"Region ID for fix efield does not exist");
int n = strlen(arg[iarg+1]) + 1;
idregion = new char[n];
strcpy(idregion,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"energy") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix efield command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
int n = strlen(&arg[iarg+1][2]) + 1;
estr = new char[n];
strcpy(estr,&arg[iarg+1][2]);
} else error->all(FLERR,"Illegal fix efield command");
iarg += 2;
} else error->all(FLERR,"Illegal fix efield command");
}

force_flag = 0;
foriginal[0] = foriginal[1] = foriginal[2] = foriginal[3] = 0.0;

maxatom = 0;
efield = NULL;
}

/* ---------------------------------------------------------------------- */

FixEfield::~FixEfield()
{
delete [] xstr;
delete [] ystr;
delete [] zstr;
delete [] estr;
delete [] idregion;
memory->destroy(efield);
}

/* ---------------------------------------------------------------------- */

int FixEfield::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
return mask;
}

/* ---------------------------------------------------------------------- */

void FixEfield::init()
{
if (!atom->q_flag) error->all(FLERR,"Fix efield requires atom attribute q");

// check variables

if (xstr) {
xvar = input->variable->find(xstr);
if (xvar < 0) error->all(FLERR,
"Variable name for fix efield does not exist");
if (input->variable->equalstyle(xvar)) xstyle = EQUAL;
else if (input->variable->atomstyle(xvar)) xstyle = ATOM;
else error->all(FLERR,"Variable for fix efield is invalid style");
}
if (ystr) {
yvar = input->variable->find(ystr);
if (yvar < 0) error->all(FLERR,
"Variable name for fix efield does not exist");
if (input->variable->equalstyle(yvar)) ystyle = EQUAL;
else if (input->variable->atomstyle(yvar)) ystyle = ATOM;
else error->all(FLERR,"Variable for fix efield is invalid style");
}
if (zstr) {
zvar = input->variable->find(zstr);
if (zvar < 0) error->all(FLERR,
"Variable name for fix efield does not exist");
if (input->variable->equalstyle(zvar)) zstyle = EQUAL;
else if (input->variable->atomstyle(zvar)) zstyle = ATOM;
else error->all(FLERR,"Variable for fix efield is invalid style");
}
if (estr) {
evar = input->variable->find(estr);
if (evar < 0)
error->all(FLERR,"Variable name for fix efield does not exist");
if (input->variable->atomstyle(evar)) estyle = ATOM;
else error->all(FLERR,"Variable for fix efield is invalid style");
} else estyle = NONE;

if (xstyle == ATOM || ystyle == ATOM || zstyle == ATOM)
varflag = ATOM;
else if (xstyle == EQUAL || ystyle == EQUAL || zstyle == EQUAL)
varflag = EQUAL;
else varflag = CONSTANT;

if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;

// set index and check validity of region

if (iregion >= 0) {
iregion = domain->find_region(idregion);
if (iregion == -1)
error->all(FLERR,"Region ID for fix efield does not exist");
}

if (xstyle == ATOM || ystyle == ATOM || zstyle == ATOM)
varflag = ATOM;
else if (xstyle == EQUAL || ystyle == EQUAL || zstyle == EQUAL)
varflag = EQUAL;
else varflag = CONSTANT;

if (varflag == CONSTANT && estyle != NONE)
error->all(FLERR,"Cannot use variable energy with "
"constant force in fix efield");
if ((varflag == EQUAL || varflag == ATOM) &&
update->whichflag == 2 && estyle == NONE)
error->all(FLERR,"Must use variable energy with fix efield");

if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;

// error checks on coarsegraining
if(force->cg_active())
error->cg(FLERR,this->style);
}

/* ---------------------------------------------------------------------- */

void FixEfield::setup(int vflag)
{
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
}
}

/* ----------------------------------------------------------------------
apply F = qE
------------------------------------------------------------------------- */

void FixEfield::post_force(int vflag)
{
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
int *mask = atom->mask;
int nlocal = atom->nlocal;

// reallocate efield array if necessary

if (varflag == ATOM && nlocal > maxatom) {
maxatom = atom->nmax;
memory->destroy(efield);
memory->create(efield,maxatom,3,"efield:efield");
}

if (varflag == CONSTANT) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (iregion >= 0 &&
!domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]))
continue;

f[i][0] += q[i]*ex;
f[i][1] += q[i]*ey;
f[i][2] += q[i]*ez;
}

// variable efield, wrap with clear/add

} else {

modify->clearstep_compute();

if (xstyle == EQUAL) ex = qe2f * input->variable->compute_equal(xvar);
else if (xstyle == ATOM && efield)
input->variable->compute_atom(xvar,igroup,&efield[0][0],3,0);
if (ystyle == EQUAL) ey = qe2f * input->variable->compute_equal(yvar);
else if (ystyle == ATOM && efield)
input->variable->compute_atom(yvar,igroup,&efield[0][1],3,0);
if (zstyle == EQUAL) ez = qe2f * input->variable->compute_equal(zvar);
else if (zstyle == ATOM && efield)
input->variable->compute_atom(zvar,igroup,&efield[0][2],3,0);

modify->addstep_compute(update->ntimestep + 1);

for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (iregion >= 0 &&
!domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]))
continue;

if (xstyle == ATOM) f[i][0] += qe2f * q[i]*efield[i][0];
else f[i][0] += q[i]*ex;
if (ystyle == ATOM) f[i][1] += qe2f * q[i]*efield[i][1];
else f[i][1] += q[i]*ey;
if (zstyle == ATOM) f[i][2] += qe2f * q[i]*efield[i][2];
else f[i][2] += q[i]*ez;
}
}
}

/* ---------------------------------------------------------------------- */

void FixEfield::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
}

/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */

double FixEfield::memory_usage()
{
double bytes = 0.0;
if (varflag == ATOM) bytes = atom->nmax*3 * sizeof(double);
return bytes;
}