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@13169 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2015-03-04 15:24:51 +00:00
parent 472b4fd62b
commit ee00c0f654
7 changed files with 55 additions and 49 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/Makefile
View File

@ -13,7 +13,7 @@ OBJ = $(SRC:.cpp=.o)
# Package variables
PACKAGE = asphere body class2 colloid dipole fld gpu granular kim \
PACKAGE = asphere body class2 colloid coreshell dipole fld gpu granular kim \
kokkos kspace manybody mc meam misc molecule mpiio opt peri poems \
qeq reax replica rigid shock snap srd voronoi xtc

4
src/QEQ/fix_qeq_slater.cpp
View File

@ -57,7 +57,7 @@ FixQEqSlater::FixQEqSlater(LAMMPS *lmp, int narg, char **arg) :
} else error->all(FLERR,"Illegal fix qeq/slater command");
}
if (streitz_flag) extract();
if (streitz_flag) extract_streitz();
}
/* ---------------------------------------------------------------------- */
@ -88,7 +88,7 @@ void FixQEqSlater::init()
/* ---------------------------------------------------------------------- */
void FixQEqSlater::extract()
void FixQEqSlater::extract_streitz()
{
int ntypes = atom->ntypes;

2
src/QEQ/fix_qeq_slater.h
View File

@ -37,7 +37,7 @@ class FixQEqSlater : public FixQEq {
void compute_H();
double calculate_H(double, double, double, double, double &);
double calculate_H_wolf(double, double, double, double, double &);
void extract();
void extract_streitz();
class PairCoulStreitz *streitz;
};

4
src/compute_chunk_atom.cpp
View File

@ -373,7 +373,7 @@ ComputeChunkAtom::ComputeChunkAtom(LAMMPS *lmp, int narg, char **arg) :
id_fix = NULL;
fixstore = NULL;
if (compress) hash = new std::map<int,int>();
if (compress) hash = new std::map<tagint,int>();
else hash = NULL;
maxvar = 0;
@ -1010,7 +1010,7 @@ void ComputeChunkAtom::compress_chunk_ids()
void ComputeChunkAtom::idring(int n, char *cbuf)
{
tagint *list = (tagint *) cbuf;
std::map<int,int> *hash = cptr->hash;
std::map<tagint,int> *hash = cptr->hash;
for (int i = 0; i < n; i++) (*hash)[list[i]] = 0;
}

2
src/compute_chunk_atom.h
View File

@ -86,7 +86,7 @@ class ComputeChunkAtom : public Compute {
int molcheck; // one-time check if all molecule atoms in chunk
int *exclude; // 1 if atom is not assigned to any chunk
std::map<int,int> *hash; // store original chunks IDs before compression
std::map<tagint,int> *hash; // store original chunks IDs before compression
// static variable for ring communication callback to access class data
// callback functions for ring communication

34
src/pair_coul_streitz.cpp
View File

@ -286,7 +286,7 @@ void PairCoulStreitz::read_file(char *file)
// strip comment, skip line if blank
if (ptr = strchr(line,'#')) *ptr = '\0';
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
if (nwords == 0) continue;
@ -305,7 +305,7 @@ void PairCoulStreitz::read_file(char *file)
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
if (ptr = strchr(line,'#')) *ptr = '\0';
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
}
@ -316,7 +316,7 @@ void PairCoulStreitz::read_file(char *file)
nwords = 0;
words[nwords++] = strtok(line," \t\n\r\f");
while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue;
while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
// ielement = 1st args
@ -355,7 +355,7 @@ void PairCoulStreitz::read_file(char *file)
void PairCoulStreitz::setup()
{
int i,j,k,m,n;
int i,m,n;
// set elem2param
@ -391,29 +391,27 @@ void PairCoulStreitz::setup()
void PairCoulStreitz::compute(int eflag, int vflag)
{
int i, j, ii, jj, inum, jnum;
int itag, jtag, itype, jtype, iparam_i,iparam_j;
int itype, jtype, iparam_i,iparam_j;
int *ilist, *jlist, *numneigh, **firstneigh;
int *tag = atom->tag;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
double xtmp, ytmp, ztmp, evdwl, ecoul, fpair;
double xtmp, ytmp, ztmp, ecoul, fpair;
double qi, qj, selfion, r, rsq, delr[3];
double zei, zej, zi, zj, ci_jfi, dci_jfi, ci_fifj, dci_fifj;
double prefactor, forcecoul, factor_coul;
double zei, zej, zj, ci_jfi, dci_jfi, ci_fifj, dci_fifj;
double forcecoul, factor_coul;
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
double qqrd2e = force->qqrd2e;
double *special_coul = force->special_coul;
evdwl = ecoul = 0.0;
ecoul = 0.0;
selfion = fpair = 0.0;
ci_jfi = ci_fifj = dci_jfi = dci_fifj = 0.0;
prefactor = forcecoul = 0.0;
forcecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = vflag_atom = 0;
@ -429,7 +427,6 @@ void PairCoulStreitz::compute(int eflag, int vflag)
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itag = tag[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
@ -437,7 +434,6 @@ void PairCoulStreitz::compute(int eflag, int vflag)
iparam_i = elem2param[itype];
qi = q[i];
zei = params[iparam_i].zeta;
zi = params[iparam_i].zcore;
// self energy: ionization + wolf sum
@ -454,7 +450,6 @@ void PairCoulStreitz::compute(int eflag, int vflag)
j = jlist[jj];
j &= NEIGHMASK;
jtag = tag[j];
jtype = map[type[j]];
iparam_j = elem2param[jtype];
qj = q[j];
@ -505,7 +500,6 @@ void PairCoulStreitz::compute(int eflag, int vflag)
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itag = tag[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
@ -513,7 +507,6 @@ void PairCoulStreitz::compute(int eflag, int vflag)
iparam_i = elem2param[itype];
qi = q[i];
zei = params[iparam_i].zeta;
zi = params[iparam_i].zcore;
// self ionizition energy, only on i atom
@ -529,7 +522,6 @@ void PairCoulStreitz::compute(int eflag, int vflag)
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtag = tag[j];
jtype = map[type[j]];
iparam_j = elem2param[jtype];
qj = q[j];
@ -587,6 +579,8 @@ double PairCoulStreitz::self(Param *param, double qi)
if (kspacetype == 1) return 1.0*qi*(s1+qi*(0.50*s2 - qqrd2e*woself));
if (kspacetype == 2) return 1.0*qi*(s1+qi*(0.50*s2));
return 0.0;
}
/* ---------------------------------------------------------------------- */
@ -673,8 +667,8 @@ void PairCoulStreitz::wolf_sum(double qi, double qj, double zj, double r,
double derfcr = exp(-a*a*r*r);
double erfcrc = erfc(a*rc);
double etmp1, etmp2, etmp3, etmp4;
double ftmp1, ftmp2, ftmp3, ftmp4;
double etmp1, etmp2, etmp3;
double ftmp1, ftmp2, ftmp3;
etmp = etmp1 = etmp2 = etmp3 = 0.0;
ftmp = ftmp1 = ftmp2 = ftmp3 = 0.0;

36
src/velocity.cpp
View File

@ -98,6 +98,28 @@ void Velocity::command(int narg, char **arg)
else if (style == RAMP) options(narg-8,&arg[8]);
else if (style == ZERO) options(narg-3,&arg[3]);
// special cases where full init and border communication must be done first
// for CREATE/SET if compute temp/cs is used
// for ZERO if fix rigid/small is used
// b/c methods invoked in the compute/fix perform forward/reverse comm
int initcomm = 0;
if (style == ZERO && rfix >= 0 &&
strcmp(modify->fix[rfix]->style,"rigid/small") == 0) initcomm = 1;
if ((style == CREATE || style == SET) && temperature &&
strcmp(temperature->style,"temp/cs2") == 0) initcomm = 1;
if (initcomm) {
lmp->init();
if (domain->triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
comm->exchange();
comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
}
// initialize velocities based on style
// create() invoked differently, so can be called externally
@ -670,37 +692,27 @@ void Velocity::ramp(int narg, char **arg)
/* ----------------------------------------------------------------------
zero linear or angular momentum of a group
if using rigid/small requires init of entire system since
its methods perform forward/reverse comm,
comm::init needs neighbor::init needs pair::init needs kspace::init, etc
also requires setup_pre_neighbor call to setup bodies
------------------------------------------------------------------------- */
void Velocity::zero(int narg, char **arg)
{
if (strcmp(arg[0],"linear") == 0) {
if (rfix < 0) zero_momentum();
else {
if (strcmp(modify->fix[rfix]->style,"rigid/small") == 0) {
lmp->init();
else if (strcmp(modify->fix[rfix]->style,"rigid/small") == 0) {
modify->fix[rfix]->setup_pre_neighbor();
modify->fix[rfix]->zero_momentum();
} else if (strstr(modify->fix[rfix]->style,"rigid")) {
modify->fix[rfix]->zero_momentum();
} else error->all(FLERR,"Velocity rigid used with non-rigid fix-ID");
}
} else if (strcmp(arg[0],"angular") == 0) {
if (rfix < 0) zero_rotation();
else {
if (strcmp(modify->fix[rfix]->style,"rigid/small") == 0) {
lmp->init();
else if (strcmp(modify->fix[rfix]->style,"rigid/small") == 0) {
modify->fix[rfix]->setup_pre_neighbor();
modify->fix[rfix]->zero_rotation();
} else if (strstr(modify->fix[rfix]->style,"rigid")) {
modify->fix[rfix]->zero_rotation();
} else error->all(FLERR,"Velocity rigid used with non-rigid fix-ID");
}
} else error->all(FLERR,"Illegal velocity command");
}