ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4028 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2010-04-23 14:24:09 +00:00
parent 91c72cea03
commit 7dbb6d500e
9 changed files with 579 additions and 231 deletions
Download Patch File Download Diff File Expand all files Collapse all files

155
src/fix_addforce.cpp
View File

@ -19,10 +19,15 @@
#include "domain.h"
#include "region.h"
#include "respa.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
enum{NONE,CONSTANT,EQUAL,ATOM};
/* ---------------------------------------------------------------------- */
FixAddForce::FixAddForce(LAMMPS *lmp, int narg, char **arg) :
@ -37,13 +42,38 @@ FixAddForce::FixAddForce(LAMMPS *lmp, int narg, char **arg) :
extscalar = 1;
extvector = 1;
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 {
xvalue = 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 {
yvalue = 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 {
zvalue = atof(arg[5]);
zstyle = CONSTANT;
}
// optional args
iregion = -1;
evar = NONE;
estr = NULL;
int iarg = 6;
while (iarg < narg) {
@ -52,11 +82,31 @@ FixAddForce::FixAddForce(LAMMPS *lmp, int narg, char **arg) :
iregion = domain->find_region(arg[iarg+1]);
if (iregion == -1) error->all("Fix addforce region ID does not exist");
iarg += 2;
} else if (strcmp(arg[iarg],"energy") == 0) {
if (iarg+2 > narg) error->all("Illegal fix addforce command");
int n = strlen(arg[iarg+1]) + 1;
estr = new char[n];
strcpy(estr,arg[iarg+1]);
iarg += 2;
} else error->all("Illegal fix addforce command");
}
force_flag = 0;
foriginal[0] = foriginal[1] = foriginal[2] = foriginal[3] = 0.0;
maxatom = 0;
sforce = NULL;
}
/* ---------------------------------------------------------------------- */
FixAddForce::~FixAddForce()
{
delete [] xstr;
delete [] ystr;
delete [] zstr;
delete [] estr;
memory->destroy_2d_double_array(sforce);
}
/* ---------------------------------------------------------------------- */
@ -75,6 +125,49 @@ int FixAddForce::setmask()
void FixAddForce::init()
{
// check variables
if (xstr) {
xvar = input->variable->find(xstr);
if (xvar < 0) error->all("Variable name for fix addforce does not exist");
if (input->variable->equalstyle(xvar)) xstyle = EQUAL;
else if (input->variable->atomstyle(xvar)) xstyle = ATOM;
else error->all("Variable for fix setforce is invalid style");
}
if (ystr) {
yvar = input->variable->find(ystr);
if (yvar < 0) error->all("Variable name for fix addforce does not exist");
if (input->variable->equalstyle(yvar)) ystyle = EQUAL;
else if (input->variable->atomstyle(yvar)) ystyle = ATOM;
else error->all("Variable for fix setforce is invalid style");
}
if (zstr) {
zvar = input->variable->find(zstr);
if (zvar < 0) error->all("Variable name for fix addforce does not exist");
if (input->variable->equalstyle(zvar)) zstyle = EQUAL;
else if (input->variable->atomstyle(zvar)) zstyle = ATOM;
else error->all("Variable for fix setforce is invalid style");
}
if (estr) {
evar = input->variable->find(estr);
if (evar < 0) error->all("Variable name for fix addforce does not exist");
if (input->variable->atomstyle(evar)) estyle = ATOM;
else error->all("Variable for fix setforce is invalid style");
}
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("Cannot use variable energy with "
"constant force in fix addforce");
if ((varflag == EQUAL || varflag == ATOM) &&
update->whichflag == 2 && estyle == NONE)
error->all("Must use variable energy with fix addforce");
if (strcmp(update->integrate_style,"respa") == 0)
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
@ -108,12 +201,24 @@ void FixAddForce::post_force(int vflag)
int *mask = atom->mask;
int nlocal = atom->nlocal;
// reallocate sforce array if necessary
if ((varflag == ATOM || estyle == ATOM) && nlocal > maxatom) {
maxatom = atom->nmax;
memory->destroy_2d_double_array(sforce);
sforce = memory->create_2d_double_array(maxatom,4,"addforce:sforce");
}
// foriginal[0] = "potential energy" for added force
// foriginal[123] = force on atoms before extra force added
foriginal[0] = foriginal[1] = foriginal[2] = foriginal[3] = 0.0;
force_flag = 0;
// foriginal[0] = - x dot f = "potential" for added force
// foriginal[123] = force on atoms before extra force added
// constant force
// potential energy = - x dot f
if (varflag == CONSTANT) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (iregion >= 0 &&
@ -128,6 +233,41 @@ void FixAddForce::post_force(int vflag)
f[i][1] += yvalue;
f[i][2] += zvalue;
}
// variable force
// potential energy = evar if defined, else 0.0
} else {
if (xstyle == EQUAL) xvalue = input->variable->compute_equal(xvar);
else if (xstyle == ATOM)
input->variable->compute_atom(xvar,igroup,&sforce[0][0],4,0);
if (ystyle == EQUAL) yvalue = input->variable->compute_equal(yvar);
else if (ystyle == ATOM)
input->variable->compute_atom(yvar,igroup,&sforce[0][1],4,0);
if (zstyle == EQUAL) zvalue = input->variable->compute_equal(zvar);
else if (zstyle == ATOM)
input->variable->compute_atom(zvar,igroup,&sforce[0][2],4,0);
if (estyle == ATOM)
input->variable->compute_atom(evar,igroup,&sforce[0][3],4,0);
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 (estyle == ATOM) foriginal[0] += sforce[i][3];
foriginal[1] += f[i][0];
foriginal[2] += f[i][1];
foriginal[3] += f[i][2];
if (xstyle == ATOM) f[i][0] += sforce[i][0];
else if (xstyle) f[i][0] += xvalue;
if (ystyle == ATOM) f[i][1] += sforce[i][1];
else if (ystyle) f[i][1] += yvalue;
if (zstyle == ATOM) f[i][2] += sforce[i][2];
else if (zstyle) f[i][2] += zvalue;
}
}
}
/* ---------------------------------------------------------------------- */
@ -173,3 +313,14 @@ double FixAddForce::compute_vector(int n)
}
return foriginal_all[n+1];
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double FixAddForce::memory_usage()
{
double bytes = 0.0;
if (varflag == ATOM) bytes = atom->nmax*4 * sizeof(double);
return bytes;
}

10
src/fix_addforce.h
View File

@ -27,6 +27,7 @@ namespace LAMMPS_NS {
class FixAddForce : public Fix {
public:
FixAddForce(class LAMMPS *, int, char **);
~FixAddForce();
int setmask();
void init();
void setup(int);
@ -36,13 +37,20 @@ class FixAddForce : public Fix {
void min_post_force(int);
double compute_scalar();
double compute_vector(int);
double memory_usage();
private:
int iregion;
double xvalue,yvalue,zvalue;
int varflag;
char *xstr,*ystr,*zstr,*estr;
int xvar,yvar,zvar,evar,xstyle,ystyle,zstyle,estyle;
int iregion;
double foriginal[4],foriginal_all[4];
int force_flag;
int nlevels_respa;
int maxatom;
double **sforce;
};
}

94
src/fix_aveforce.cpp
View File

@ -20,10 +20,14 @@
#include "domain.h"
#include "region.h"
#include "respa.h"
#include "input.h"
#include "variable.h"
#include "error.h"
using namespace LAMMPS_NS;
enum{NONE,CONSTANT,EQUAL};
/* ---------------------------------------------------------------------- */
FixAveForce::FixAveForce(LAMMPS *lmp, int narg, char **arg) :
@ -36,13 +40,36 @@ FixAveForce::FixAveForce(LAMMPS *lmp, int narg, char **arg) :
global_freq = 1;
extvector = 1;
xflag = yflag = zflag = 1;
if (strcmp(arg[3],"NULL") == 0) xflag = 0;
else xvalue = atof(arg[3]);
if (strcmp(arg[4],"NULL") == 0) yflag = 0;
else yvalue = atof(arg[4]);
if (strcmp(arg[5],"NULL") == 0) zflag = 0;
else zvalue = atof(arg[5]);
if (strstr(arg[3],"v_") == arg[3]) {
int n = strlen(&arg[3][2]) + 1;
xstr = new char[n];
strcpy(xstr,&arg[3][2]);
} else if (strcmp(arg[3],"NULL") == 0) {
xstyle = NONE;
} else {
xvalue = 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 if (strcmp(arg[4],"NULL") == 0) {
ystyle = NONE;
} else {
yvalue = 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 if (strcmp(arg[5],"NULL") == 0) {
zstyle = NONE;
} else {
zvalue = atof(arg[5]);
zstyle = CONSTANT;
}
// optional args
@ -65,6 +92,15 @@ FixAveForce::FixAveForce(LAMMPS *lmp, int narg, char **arg) :
/* ---------------------------------------------------------------------- */
FixAveForce::~FixAveForce()
{
delete [] xstr;
delete [] ystr;
delete [] zstr;
}
/* ---------------------------------------------------------------------- */
int FixAveForce::setmask()
{
int mask = 0;
@ -78,6 +114,30 @@ int FixAveForce::setmask()
void FixAveForce::init()
{
// check variables
if (xstr) {
xvar = input->variable->find(xstr);
if (xvar < 0) error->all("Variable name for fix aveforce does not exist");
if (input->variable->equalstyle(xvar)) xstyle = EQUAL;
else error->all("Variable for fix aveforce is invalid style");
}
if (ystr) {
yvar = input->variable->find(ystr);
if (yvar < 0) error->all("Variable name for fix aveforce does not exist");
if (input->variable->equalstyle(yvar)) ystyle = EQUAL;
else error->all("Variable for fix aveforce is invalid style");
}
if (zstr) {
zvar = input->variable->find(zstr);
if (zvar < 0) error->all("Variable name for fix aveforce does not exist");
if (input->variable->equalstyle(zvar)) zstyle = EQUAL;
else error->all("Variable for fix aveforce is invalid style");
}
if (xstyle == EQUAL || ystyle == EQUAL || zstyle == EQUAL) varflag = EQUAL;
else varflag = CONSTANT;
if (strcmp(update->integrate_style,"respa") == 0)
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
@ -130,13 +190,19 @@ void FixAveForce::post_force(int vflag)
}
// average the force on participating atoms
// add in requested amount
// add in requested amount, computed via variable evaluation if necessary
MPI_Allreduce(foriginal,foriginal_all,4,MPI_DOUBLE,MPI_SUM,world);
int ncount = static_cast<int> (foriginal_all[3]);
if (ncount == 0) return;
if (varflag == EQUAL) {
if (xstyle == EQUAL) xvalue = input->variable->compute_equal(xvar);
if (ystyle == EQUAL) yvalue = input->variable->compute_equal(yvar);
if (zstyle == EQUAL) zvalue = input->variable->compute_equal(zvar);
}
double fave[3];
fave[0] = foriginal_all[0]/ncount + xvalue;
fave[1] = foriginal_all[1]/ncount + yvalue;
@ -151,9 +217,9 @@ void FixAveForce::post_force(int vflag)
!domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]))
continue;
if (xflag) f[i][0] = fave[0];
if (yflag) f[i][1] = fave[1];
if (zflag) f[i][2] = fave[2];
if (xstyle) f[i][0] = fave[0];
if (ystyle) f[i][1] = fave[1];
if (zstyle) f[i][2] = fave[2];
}
}
@ -202,9 +268,9 @@ void FixAveForce::post_force_respa(int vflag, int ilevel, int iloop)
!domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]))
continue;
if (xflag) f[i][0] = fave[0];
if (yflag) f[i][1] = fave[1];
if (zflag) f[i][2] = fave[2];
if (xstyle) f[i][0] = fave[0];
if (ystyle) f[i][1] = fave[1];
if (zstyle) f[i][2] = fave[2];
}
}
}

6
src/fix_aveforce.h
View File

@ -27,6 +27,7 @@ namespace LAMMPS_NS {
class FixAveForce : public Fix {
public:
FixAveForce(class LAMMPS *, int, char **);
~FixAveForce();
int setmask();
void init();
void setup(int);
@ -37,8 +38,11 @@ class FixAveForce : public Fix {
double compute_vector(int);
private:
int xflag,yflag,zflag,iregion;
double xvalue,yvalue,zvalue;
int varflag;
char *xstr,*ystr,*zstr;
int xvar,yvar,zvar,xstyle,ystyle,zstyle;
int iregion;
double foriginal_all[4];
int nlevels_respa;
};

70
src/fix_efield.cpp
View File

@ -30,7 +30,7 @@
using namespace LAMMPS_NS;
enum{NONE,EQUAL,ATOM};
enum{CONSTANT,EQUAL,ATOM};
/* ---------------------------------------------------------------------- */
@ -46,19 +46,28 @@ FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) :
int n = strlen(&arg[3][2]) + 1;
xstr = new char[n];
strcpy(xstr,&arg[3][2]);
} else ex = efactor * atof(arg[3]);
} else {
ex = efactor * atof(arg[3]);
xstyle = CONSTANT;
}
if (strstr(arg[4],"v_") == arg[4]) {
int n = strlen(&arg[4][2]) + 1;
xstr = new char[n];
strcpy(xstr,&arg[4][2]);
} else ey = efactor * atof(arg[4]);
} else {
ey = efactor * atof(arg[4]);
ystyle = CONSTANT;
}
if (strstr(arg[5],"v_") == arg[5]) {
int n = strlen(&arg[5][2]) + 1;
xstr = new char[n];
strcpy(xstr,&arg[5][2]);
} else ez = efactor * atof(arg[5]);
} else {
ez = efactor * atof(arg[5]);
zstyle = CONSTANT;
}
maxatom = 0;
efield = NULL;
@ -98,30 +107,30 @@ void FixEfield::init()
if (xstr) {
xvar = input->variable->find(xstr);
if (xvar < 0) error->all("Variable name for fix efield does not exist");
if (input->variable->equalstyle(xvar)) xvarstyle = EQUAL;
else if (input->variable->atomstyle(xvar)) xvarstyle = ATOM;
if (input->variable->equalstyle(xvar)) xstyle = EQUAL;
else if (input->variable->atomstyle(xvar)) xstyle = ATOM;
else error->all("Variable for fix efield is invalid style");
} else xvarstyle = NONE;
}
if (ystr) {
yvar = input->variable->find(ystr);
if (yvar < 0) error->all("Variable name for fix efield does not exist");
if (input->variable->equalstyle(yvar)) yvarstyle = EQUAL;
else if (input->variable->atomstyle(yvar)) yvarstyle = ATOM;
if (input->variable->equalstyle(yvar)) ystyle = EQUAL;
else if (input->variable->atomstyle(yvar)) ystyle = ATOM;
else error->all("Variable for fix efield is invalid style");
} else yvarstyle = NONE;
}
if (zstr) {
zvar = input->variable->find(zstr);
if (zvar < 0) error->all("Variable name for fix efield does not exist");
if (input->variable->equalstyle(zvar)) zvarstyle = EQUAL;
else if (input->variable->atomstyle(zvar)) zvarstyle = ATOM;
if (input->variable->equalstyle(zvar)) zstyle = EQUAL;
else if (input->variable->atomstyle(zvar)) zstyle = ATOM;
else error->all("Variable for fix efield is invalid style");
} else zvarstyle = NONE;
}
if (xvarstyle == ATOM || yvarstyle == ATOM || zvarstyle == ATOM)
if (xstyle == ATOM || ystyle == ATOM || zstyle == ATOM)
varflag = ATOM;
else if (xvarstyle == EQUAL || yvarstyle == EQUAL || zvarstyle == EQUAL)
else if (xstyle == EQUAL || ystyle == EQUAL || zstyle == EQUAL)
varflag = EQUAL;
else varflag = NONE;
else varflag = CONSTANT;
if (strcmp(update->integrate_style,"respa") == 0)
nlevels_respa = ((Respa *) update->integrate)->nlevels;
@ -159,7 +168,7 @@ void FixEfield::post_force(int vflag)
efield = memory->create_2d_double_array(maxatom,3,"efield:efield");
}
if (varflag == NONE) {
if (varflag == CONSTANT) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
f[i][0] += q[i]*ex;
@ -168,32 +177,23 @@ void FixEfield::post_force(int vflag)
}
} else {
if (xstr) {
if (xvarstyle == EQUAL)
ex = efactor * input->variable->compute_equal(xvar);
else
if (xstyle == EQUAL) ex = efactor * input->variable->compute_equal(xvar);
else if (xstyle == ATOM)
input->variable->compute_atom(xvar,igroup,&efield[0][0],3,0);
}
if (ystr) {
if (yvarstyle == EQUAL)
ey = efactor * input->variable->compute_equal(yvar);
else
if (ystyle == EQUAL) ey = efactor * input->variable->compute_equal(yvar);
else if (ystyle == ATOM)
input->variable->compute_atom(yvar,igroup,&efield[0][1],3,0);
}
if (zstr) {
if (zvarstyle == EQUAL)
ez = efactor * input->variable->compute_equal(zvar);
else
if (zstyle == EQUAL) ez = efactor * input->variable->compute_equal(zvar);
else if (zstyle == ATOM)
input->variable->compute_atom(zvar,igroup,&efield[0][2],3,0);
}
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (xvarstyle == ATOM) f[i][0] += q[i]*efield[i][0];
if (xstyle == ATOM) f[i][0] += q[i]*efield[i][0];
else f[i][0] += q[i]*ex;
if (yvarstyle == ATOM) f[i][1] += q[i]*efield[i][1];
if (ystyle == ATOM) f[i][1] += q[i]*efield[i][1];
else f[i][1] += q[i]*ey;
if (zvarstyle == ATOM) f[i][2] += q[i]*efield[i][2];
if (zstyle == ATOM) f[i][2] += q[i]*efield[i][2];
else f[i][2] += q[i]*ez;
}
}

2
src/fix_efield.h
View File

@ -39,7 +39,7 @@ class FixEfield : public Fix {
double ex,ey,ez;
int varflag;
char *xstr,*ystr,*zstr;
int xvar,yvar,zvar,xvarstyle,yvarstyle,zvarstyle;
int xvar,yvar,zvar,xstyle,ystyle,zstyle;
int nlevels_respa;
double efactor;

157
src/fix_set_force.cpp
View File

@ -1,157 +0,0 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "stdlib.h"
#include "fix_set_force.h"
#include "atom.h"
#include "update.h"
#include "respa.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
FixSetForce::FixSetForce(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg != 6) error->all("Illegal fix setforce command");
vector_flag = 1;
size_vector = 3;
global_freq = 1;
extvector = 1;
flagx = flagy = flagz = 1;
if (strcmp(arg[3],"NULL") == 0) flagx = 0;
else xvalue = atof(arg[3]);
if (strcmp(arg[4],"NULL") == 0) flagy = 0;
else yvalue = atof(arg[4]);
if (strcmp(arg[5],"NULL") == 0) flagz = 0;
else zvalue = atof(arg[5]);
force_flag = 0;
foriginal[0] = foriginal[1] = foriginal[2] = 0.0;
}
/* ---------------------------------------------------------------------- */
int FixSetForce::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= MIN_POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixSetForce::init()
{
if (strcmp(update->integrate_style,"respa") == 0)
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ---------------------------------------------------------------------- */
void FixSetForce::setup(int vflag)
{
if (strcmp(update->integrate_style,"verlet") == 0)
post_force(vflag);
else
for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) {
((Respa *) update->integrate)->copy_flevel_f(ilevel);
post_force_respa(vflag,ilevel,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel);
}
}
/* ---------------------------------------------------------------------- */
void FixSetForce::min_setup(int vflag)
{
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixSetForce::post_force(int vflag)
{
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
foriginal[0] = foriginal[1] = foriginal[2] = 0.0;
force_flag = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
foriginal[0] += f[i][0];
foriginal[1] += f[i][1];
foriginal[2] += f[i][2];
if (flagx) f[i][0] = xvalue;
if (flagy) f[i][1] = yvalue;
if (flagz) f[i][2] = zvalue;
}
}
/* ---------------------------------------------------------------------- */
void FixSetForce::post_force_respa(int vflag, int ilevel, int iloop)
{
// set force to desired value on outermost level, 0.0 on other levels
if (ilevel == nlevels_respa-1) post_force(vflag);
else {
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
foriginal[0] = foriginal[1] = foriginal[2] = 0.0;
force_flag = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
foriginal[0] += f[i][0];
foriginal[1] += f[i][1];
foriginal[2] += f[i][2];
if (flagx) f[i][0] = 0.0;
if (flagy) f[i][1] = 0.0;
if (flagz) f[i][2] = 0.0;
}
}
}
/* ---------------------------------------------------------------------- */
void FixSetForce::min_post_force(int vflag)
{
post_force(vflag);
}
/* ----------------------------------------------------------------------
return components of total force on fix group before force was changed
------------------------------------------------------------------------- */
double FixSetForce::compute_vector(int n)
{
// only sum across procs one time
if (force_flag == 0) {
MPI_Allreduce(foriginal,foriginal_all,3,MPI_DOUBLE,MPI_SUM,world);
force_flag = 1;
}
return foriginal_all[n];
}

269
src/fix_setforce.cpp Normal file
View File

@ -0,0 +1,269 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "stdlib.h"
#include "fix_setforce.h"
#include "atom.h"
#include "update.h"
#include "respa.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
enum{NONE,CONSTANT,EQUAL,ATOM};
/* ---------------------------------------------------------------------- */
FixSetForce::FixSetForce(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg != 6) error->all("Illegal fix setforce command");
vector_flag = 1;
size_vector = 3;
global_freq = 1;
extvector = 1;
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 if (strcmp(arg[3],"NULL") == 0) {
xstyle = NONE;
} else {
xvalue = 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 if (strcmp(arg[4],"NULL") == 0) {
ystyle = NONE;
} else {
yvalue = 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 if (strcmp(arg[5],"NULL") == 0) {
zstyle = NONE;
} else {
zvalue = atof(arg[5]);
zstyle = CONSTANT;
}
force_flag = 0;
foriginal[0] = foriginal[1] = foriginal[2] = 0.0;
maxatom = 0;
sforce = NULL;
}
/* ---------------------------------------------------------------------- */
FixSetForce::~FixSetForce()
{
delete [] xstr;
delete [] ystr;
delete [] zstr;
memory->destroy_2d_double_array(sforce);
}
/* ---------------------------------------------------------------------- */
int FixSetForce::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= MIN_POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixSetForce::init()
{
// check variables
if (xstr) {
xvar = input->variable->find(xstr);
if (xvar < 0) error->all("Variable name for fix setforce does not exist");
if (input->variable->equalstyle(xvar)) xstyle = EQUAL;
else if (input->variable->atomstyle(xvar)) xstyle = ATOM;
else error->all("Variable for fix setforce is invalid style");
}
if (ystr) {
yvar = input->variable->find(ystr);
if (yvar < 0) error->all("Variable name for fix setforce does not exist");
if (input->variable->equalstyle(yvar)) ystyle = EQUAL;
else if (input->variable->atomstyle(yvar)) ystyle = ATOM;
else error->all("Variable for fix setforce is invalid style");
}
if (zstr) {
zvar = input->variable->find(zstr);
if (zvar < 0) error->all("Variable name for fix setforce does not exist");
if (input->variable->equalstyle(zvar)) zstyle = EQUAL;
else if (input->variable->atomstyle(zvar)) zstyle = ATOM;
else error->all("Variable for fix setforce is invalid style");
}
if (xstyle == ATOM || ystyle == ATOM || zstyle == ATOM)
varflag = ATOM;
else if (xstyle == EQUAL || ystyle == EQUAL || zstyle == EQUAL)
varflag = EQUAL;
else varflag = CONSTANT;
if (strcmp(update->integrate_style,"respa") == 0)
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ---------------------------------------------------------------------- */
void FixSetForce::setup(int vflag)
{
if (strcmp(update->integrate_style,"verlet") == 0)
post_force(vflag);
else
for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) {
((Respa *) update->integrate)->copy_flevel_f(ilevel);
post_force_respa(vflag,ilevel,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel);
}
}
/* ---------------------------------------------------------------------- */
void FixSetForce::min_setup(int vflag)
{
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixSetForce::post_force(int vflag)
{
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// reallocate sforce array if necessary
if (varflag == ATOM && nlocal > maxatom) {
maxatom = atom->nmax;
memory->destroy_2d_double_array(sforce);
sforce = memory->create_2d_double_array(maxatom,3,"setforce:sforce");
}
foriginal[0] = foriginal[1] = foriginal[2] = 0.0;
force_flag = 0;
if (varflag == CONSTANT) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
foriginal[0] += f[i][0];
foriginal[1] += f[i][1];
foriginal[2] += f[i][2];
if (xstyle) f[i][0] = xvalue;
if (ystyle) f[i][1] = yvalue;
if (zstyle) f[i][2] = zvalue;
}
} else {
if (xstyle == EQUAL) xvalue = input->variable->compute_equal(xvar);
else if (xstyle == ATOM)
input->variable->compute_atom(xvar,igroup,&sforce[0][0],3,0);
if (ystyle == EQUAL) yvalue = input->variable->compute_equal(yvar);
else if (ystyle == ATOM)
input->variable->compute_atom(yvar,igroup,&sforce[0][1],3,0);
if (zstyle == EQUAL) zvalue = input->variable->compute_equal(zvar);
else if (zstyle == ATOM)
input->variable->compute_atom(zvar,igroup,&sforce[0][2],3,0);
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
foriginal[0] += f[i][0];
foriginal[1] += f[i][1];
foriginal[2] += f[i][2];
if (xstyle == ATOM) f[i][0] = sforce[i][0];
else if (xstyle) f[i][0] = xvalue;
if (ystyle == ATOM) f[i][1] = sforce[i][1];
else if (ystyle) f[i][1] = yvalue;
if (zstyle == ATOM) f[i][2] = sforce[i][2];
else if (zstyle) f[i][2] = zvalue;
}
}
}
/* ---------------------------------------------------------------------- */
void FixSetForce::post_force_respa(int vflag, int ilevel, int iloop)
{
// set force to desired value on outermost level, 0.0 on other levels
if (ilevel == nlevels_respa-1) post_force(vflag);
else {
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (xstyle) f[i][0] = 0.0;
if (ystyle) f[i][1] = 0.0;
if (zstyle) f[i][2] = 0.0;
}
}
}
/* ---------------------------------------------------------------------- */
void FixSetForce::min_post_force(int vflag)
{
post_force(vflag);
}
/* ----------------------------------------------------------------------
return components of total force on fix group before force was changed
------------------------------------------------------------------------- */
double FixSetForce::compute_vector(int n)
{
// only sum across procs one time
if (force_flag == 0) {
MPI_Allreduce(foriginal,foriginal_all,3,MPI_DOUBLE,MPI_SUM,world);
force_flag = 1;
}
return foriginal_all[n];
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double FixSetForce::memory_usage()
{
double bytes = 0.0;
if (varflag == ATOM) bytes = atom->nmax*3 * sizeof(double);
return bytes;
}

9
src/fix_set_force.h → src/fix_setforce.h
View File

@ -27,6 +27,7 @@ namespace LAMMPS_NS {
class FixSetForce : public Fix {
public:
FixSetForce(class LAMMPS *, int, char **);
~FixSetForce();
int setmask();
void init();
void setup(int);
@ -35,13 +36,19 @@ class FixSetForce : public Fix {
void post_force_respa(int, int, int);
void min_post_force(int);
double compute_vector(int);
double memory_usage();
private:
int flagx,flagy,flagz;
double xvalue,yvalue,zvalue;
int varflag;
char *xstr,*ystr,*zstr;
int xvar,yvar,zvar,xstyle,ystyle,zstyle;
double foriginal[3],foriginal_all[3];
int force_flag;
int nlevels_respa;
int maxatom;
double **sforce;
};
}