jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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

This commit is contained in:
sjplimp 2014-10-21 16:36:14 +00:00
parent c90746a20c
commit 4ffd3808e0
6 changed files with 139 additions and 152 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
src/QEQ/fix_qeq.cpp
View File

@ -269,6 +269,9 @@ void FixQEq::init_list(int id, NeighList *ptr)
void FixQEq::setup_pre_force(int vflag)
{
if (force->newton_pair == 0)
error->all(FLERR,"QEQ with 'newton pair off' not supported");
neighbor->build_one(list);
deallocate_storage();
@ -405,11 +408,9 @@ void FixQEq::sparse_matvec( sparse_matrix *A, double *x, double *b )
{
int i, j, itr_j;
int nn, NN;
int *ilist;
nn = atom->nlocal;
NN = atom->nlocal + atom->nghost;
ilist = list->ilist;
for( i = 0; i < nn; ++i ) {
if (atom->mask[i] & groupbit)
@ -696,7 +697,7 @@ void FixQEq::vector_add( double* dest, double c, double* v, int k )
void FixQEq::read_file(char *file)
{
int i,itype,ntypes;
int itype,ntypes;
int params_per_line = 6;
char **words = new char*[params_per_line+1];

1
src/RIGID/fix_rigid_small.cpp
View File

@ -1734,6 +1734,7 @@ void FixRigidSmall::setup_bodies_static()
// dorientflag = 1 if any particle stores dipole orientation
if (extended) {
grow_arrays(atom->nmax);
if (atom->ellipsoid_flag) orientflag = 4;
if (atom->line_flag) orientflag = 1;
if (atom->tri_flag) orientflag = 4;

255
src/USER-FEP/compute_fep.cpp
View File

@ -28,6 +28,8 @@
#include "pair_hybrid.h"
#include "kspace.h"
#include "input.h"
#include "fix.h"
#include "modify.h"
#include "variable.h"
#include "timer.h"
#include "memory.h"
@ -39,9 +41,6 @@ using namespace LAMMPS_NS;
enum{PAIR,ATOM};
enum{CHARGE};
#undef FEP_DEBUG
#undef FEP_MAXDEBUG
/* ---------------------------------------------------------------------- */
ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) :
@ -66,11 +65,13 @@ ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) :
int iarg = 4;
while (iarg < narg) {
if (strcmp(arg[iarg],"pair") == 0) {
if (iarg+6 > narg) error->all(FLERR,"Illegal pair attribute in compute fep");
if (iarg+6 > narg) error->all(FLERR,
"Illegal pair attribute in compute fep");
npert++;
iarg += 6;
} else if (strcmp(arg[iarg],"atom") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal atom attribute in compute fep");
if (iarg+4 > narg) error->all(FLERR,
"Illegal atom attribute in compute fep");
npert++;
iarg += 4;
} else break;
@ -159,6 +160,7 @@ ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) :
allocate_storage();
fixgpu = NULL;
}
/* ---------------------------------------------------------------------- */
@ -244,6 +246,11 @@ void ComputeFEP::init()
"compute tail corrections");
}
// detect if package gpu is present
int ifixgpu = modify->find_fix("package_gpu");
if (ifixgpu >= 0) fixgpu = modify->fix[ifixgpu];
if (comm->me == 0) {
if (screen) {
fprintf(screen, "FEP settings ...\n");
@ -282,6 +289,9 @@ void ComputeFEP::compute_vector()
{
double pe0,pe1;
eflag = 1;
vflag = 0;
invoked_vector = update->ntimestep;
if (atom->nmax > nmax) { // reallocate working arrays if necessary
@ -294,26 +304,36 @@ void ComputeFEP::compute_vector()
timer->stamp();
if (force->pair && force->pair->compute_flag) {
force->pair->compute(1,0);
force->pair->compute(eflag,vflag);
timer->stamp(TIME_PAIR);
}
if (force->kspace && force->kspace->compute_flag) {
force->kspace->compute(1,0);
if (chgflag && force->kspace && force->kspace->compute_flag) {
force->kspace->compute(eflag,vflag);
timer->stamp(TIME_KSPACE);
}
// accumulate force/energy/virial from /gpu pair styles
if (fixgpu) fixgpu->post_force(vflag);
pe0 = compute_epair();
perturb_params();
timer->stamp();
if (force->pair && force->pair->compute_flag) {
force->pair->compute(1,0);
force->pair->compute(eflag,vflag);
timer->stamp(TIME_PAIR);
}
if (force->kspace && force->kspace->compute_flag) {
force->kspace->compute(1,0);
if (chgflag && force->kspace && force->kspace->compute_flag) {
force->kspace->compute(eflag,vflag);
timer->stamp(TIME_KSPACE);
}
// accumulate force/energy/virial from /gpu pair styles
// this is required as to empty the answer queue,
// otherwise the force compute on the GPU in the next step would be incorrect
if (fixgpu) fixgpu->post_force(vflag);
pe1 = compute_epair();
restore_qfev(); // restore charge, force, energy, virial array values
@ -324,12 +344,6 @@ void ComputeFEP::compute_vector()
vector[2] = domain->xprd * domain->yprd * domain->zprd;
if (volumeflag)
vector[1] *= vector[2];
#ifdef FEP_DEBUG
if (comm->me == 0 && screen)
fprintf(screen, "FEP U0 = %f U1 = %f DU = %f exp(-DU/kT) = %f\n",
pe0,pe1,vector[0],vector[1]);
#endif
}
@ -351,7 +365,7 @@ double ComputeFEP::compute_epair()
eng_pair += force->pair->etail / volume;
}
if (force->kspace) eng_pair += force->kspace->energy;
if (chgflag && force->kspace) eng_pair += force->kspace->energy;
return eng_pair;
}
@ -374,18 +388,6 @@ void ComputeFEP::perturb_params()
for (i = pert->ilo; i <= pert->ihi; i++)
for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
pert->array[i][j] = pert->array_orig[i][j] + delta;
#ifdef FEP_MAXDEBUG
if (comm->me == 0 && screen) {
fprintf(screen, "###FEP change %s %s, delta = %f\n",
pert->pstyle, pert->pparam, delta);
fprintf(screen, "###FEP I J old_param new_param\n");
for (i = pert->ilo; i <= pert->ihi; i++)
for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
fprintf(screen, "###FEP %2d %2d %9.5f %9.5f\n", i, j,
pert->array_orig[i][j], pert->array[i][j]);
}
#endif
} else if (pert->which == ATOM) {
@ -400,18 +402,6 @@ void ComputeFEP::perturb_params()
if (mask[i] & groupbit)
q[i] += delta;
#ifdef FEP_MAXDEBUG
if (comm->me == 0 && screen) {
fprintf(screen, "###FEP change charge, delta = %f\n", delta);
fprintf(screen, "###FEP atom I old_q new_q\n");
for (i = 0; i < atom->nlocal; i++)
if (atype[i] >= pert->ilo && atype[i] <= pert->ihi)
if (mask[i] & groupbit)
fprintf(screen, "###FEP %5d %2d %9.5f %9.5f\n", i, atype[i],
q_orig[i], q[i]);
}
#endif
}
}
}
@ -421,6 +411,16 @@ void ComputeFEP::perturb_params()
// and also offset and tail corrections
if (pairflag) force->pair->reinit();
// when perturbing charge and using kspace,
// need to recompute additional params in kspace->setup()
// first backup the state of kspace->qsum_update_flag
if (chgflag && force->kspace) {
sys_qsum_update_flag = force->kspace->qsum_update_flag;
force->kspace->qsum_update_flag = 1;
force->kspace->setup();
}
}
@ -457,40 +457,15 @@ void ComputeFEP::restore_params()
for (i = pert->ilo; i <= pert->ihi; i++)
for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
pert->array[i][j] = pert->array_orig[i][j];
#ifdef FEP_MAXDEBUG
if (comm->me == 0 && screen) {
fprintf(screen, "###FEP restore %s %s\n", pert->pstyle, pert->pparam);
fprintf(screen, "###FEP I J param\n");
for (i = pert->ilo; i <= pert->ihi; i++)
for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
fprintf(screen, "###FEP %2d %2d %9.5f\n", i, j, pert->array[i][j]);
}
#endif
}
#ifdef FEP_MAXDEBUG
if (pert->which == ATOM) {
if (comm->me == 0 && screen) {
int *atype = atom->type;
double *q = atom->q;
int *mask = atom->mask;
int natom = atom->nlocal;
fprintf(screen, "###FEP restore charge\n");
fprintf(screen, "###FEP atom I q\n");
for (i = 0; i < natom; i++)
if (atype[i] >= pert->ilo && atype[i] <= pert->ihi)
if (mask[i] & groupbit)
fprintf(screen, "###FEP %5d %2d %9.5f\n", i, atype[i], q[i]);
}
}
#endif
}
// re-initialize pair styles if any PAIR settings were changed
// this resets other coeffs that may depend on changed values,
// and also offset and tail corrections
if (pairflag) force->pair->reinit();
if (chgflag && force->kspace) {
force->kspace->setup();
// restore kspace->qsum_update_flag to original state
force->kspace->qsum_update_flag = sys_qsum_update_flag;
}
}
@ -501,14 +476,15 @@ void ComputeFEP::restore_params()
void ComputeFEP::allocate_storage()
{
nmax = atom->nmax;
if (chgflag)
memory->create(q_orig,nmax,"fep:q_orig");
memory->create(f_orig,nmax,3,"fep:f_orig");
memory->create(peatom_orig,nmax,"fep:peatom_orig");
memory->create(pvatom_orig,nmax,6,"fep:pvatom_orig");
if (force->kspace) {
memory->create(keatom_orig,nmax,"fep:keatom_orig");
memory->create(kvatom_orig,nmax,6,"fep:kvatom_orig");
if (chgflag) {
memory->create(q_orig,nmax,"fep:q_orig");
if (force->kspace) {
memory->create(keatom_orig,nmax,"fep:keatom_orig");
memory->create(kvatom_orig,nmax,6,"fep:kvatom_orig");
}
}
}
@ -516,14 +492,15 @@ void ComputeFEP::allocate_storage()
void ComputeFEP::deallocate_storage()
{
if (chgflag)
memory->destroy(q_orig);
memory->destroy(f_orig);
memory->destroy(peatom_orig);
memory->destroy(pvatom_orig);
if (force->kspace) {
memory->destroy(keatom_orig);
memory->destroy(kvatom_orig);
if (chgflag) {
memory->destroy(q_orig);
if (force->kspace) {
memory->destroy(keatom_orig);
memory->destroy(kvatom_orig);
}
}
}
@ -541,12 +518,6 @@ void ComputeFEP::backup_qfev()
if (force->newton || force->kspace->tip4pflag)
natom += atom->nghost;
if (chgflag) {
double *q = atom->q;
for (i = 0; i < nall; i++)
q_orig[i] = q[i];
}
double **f = atom->f;
for (i = 0; i < natom; i++) {
f_orig[i][0] = f[i][0];
@ -581,29 +552,35 @@ void ComputeFEP::backup_qfev()
}
}
if (force->kspace) {
energy_orig = force->kspace->energy;
kvirial_orig[0] = force->kspace->virial[0];
kvirial_orig[1] = force->kspace->virial[1];
kvirial_orig[2] = force->kspace->virial[2];
kvirial_orig[3] = force->kspace->virial[3];
kvirial_orig[4] = force->kspace->virial[4];
kvirial_orig[5] = force->kspace->virial[5];
if (update->eflag_atom) {
double *keatom = force->kspace->eatom;
for (i = 0; i < natom; i++)
keatom_orig[i] = keatom[i];
}
if (update->vflag_atom) {
double **kvatom = force->kspace->vatom;
for (i = 0; i < natom; i++) {
kvatom_orig[i][0] = kvatom[i][0];
kvatom_orig[i][1] = kvatom[i][1];
kvatom_orig[i][2] = kvatom[i][2];
kvatom_orig[i][3] = kvatom[i][3];
kvatom_orig[i][4] = kvatom[i][4];
kvatom_orig[i][5] = kvatom[i][5];
if (chgflag) {
double *q = atom->q;
for (i = 0; i < nall; i++)
q_orig[i] = q[i];
if (force->kspace) {
energy_orig = force->kspace->energy;
kvirial_orig[0] = force->kspace->virial[0];
kvirial_orig[1] = force->kspace->virial[1];
kvirial_orig[2] = force->kspace->virial[2];
kvirial_orig[3] = force->kspace->virial[3];
kvirial_orig[4] = force->kspace->virial[4];
kvirial_orig[5] = force->kspace->virial[5];
if (update->eflag_atom) {
double *keatom = force->kspace->eatom;
for (i = 0; i < natom; i++)
keatom_orig[i] = keatom[i];
}
if (update->vflag_atom) {
double **kvatom = force->kspace->vatom;
for (i = 0; i < natom; i++) {
kvatom_orig[i][0] = kvatom[i][0];
kvatom_orig[i][1] = kvatom[i][1];
kvatom_orig[i][2] = kvatom[i][2];
kvatom_orig[i][3] = kvatom[i][3];
kvatom_orig[i][4] = kvatom[i][4];
kvatom_orig[i][5] = kvatom[i][5];
}
}
}
}
@ -620,12 +597,6 @@ void ComputeFEP::restore_qfev()
if (force->newton || force->kspace->tip4pflag)
natom += atom->nghost;
if (chgflag) {
double *q = atom->q;
for (i = 0; i < nall; i++)
q[i] = q_orig[i];
}
double **f = atom->f;
for (i = 0; i < natom; i++) {
f[i][0] = f_orig[i][0];
@ -660,29 +631,35 @@ void ComputeFEP::restore_qfev()
}
}
if (force->kspace) {
force->kspace->energy = energy_orig;
force->kspace->virial[0] = kvirial_orig[0];
force->kspace->virial[1] = kvirial_orig[1];
force->kspace->virial[2] = kvirial_orig[2];
force->kspace->virial[3] = kvirial_orig[3];
force->kspace->virial[4] = kvirial_orig[4];
force->kspace->virial[5] = kvirial_orig[5];
if (chgflag) {
double *q = atom->q;
for (i = 0; i < nall; i++)
q[i] = q_orig[i];
if (force->kspace) {
force->kspace->energy = energy_orig;
force->kspace->virial[0] = kvirial_orig[0];
force->kspace->virial[1] = kvirial_orig[1];
force->kspace->virial[2] = kvirial_orig[2];
force->kspace->virial[3] = kvirial_orig[3];
force->kspace->virial[4] = kvirial_orig[4];
force->kspace->virial[5] = kvirial_orig[5];
if (update->eflag_atom) {
double *keatom = force->kspace->eatom;
for (i = 0; i < natom; i++)
keatom[i] = keatom_orig[i];
}
if (update->vflag_atom) {
double **kvatom = force->kspace->vatom;
for (i = 0; i < natom; i++) {
kvatom[i][0] = kvatom_orig[i][0];
kvatom[i][1] = kvatom_orig[i][1];
kvatom[i][2] = kvatom_orig[i][2];
kvatom[i][3] = kvatom_orig[i][3];
kvatom[i][4] = kvatom_orig[i][4];
kvatom[i][5] = kvatom_orig[i][5];
if (update->eflag_atom) {
double *keatom = force->kspace->eatom;
for (i = 0; i < natom; i++)
keatom[i] = keatom_orig[i];
}
if (update->vflag_atom) {
double **kvatom = force->kspace->vatom;
for (i = 0; i < natom; i++) {
kvatom[i][0] = kvatom_orig[i][0];
kvatom[i][1] = kvatom_orig[i][1];
kvatom[i][2] = kvatom_orig[i][2];
kvatom[i][3] = kvatom_orig[i][3];
kvatom[i][4] = kvatom_orig[i][4];
kvatom[i][5] = kvatom_orig[i][5];
}
}
}
}

4
src/USER-FEP/compute_fep.h
View File

@ -41,6 +41,8 @@ class ComputeFEP : public Compute {
int chgflag;
int tailflag, volumeflag;
int fepinitflag;
int eflag, vflag;
int sys_qsum_update_flag;
double temp_fep;
int nmax;
@ -53,6 +55,8 @@ class ComputeFEP : public Compute {
double kvirial_orig[6];
double *keatom_orig,**kvatom_orig;
class Fix *fixgpu;
struct Perturb {
int which,ivar;
char *var;

23
src/USER-FEP/fix_adapt_fep.cpp
View File

@ -171,6 +171,11 @@ FixAdaptFEP::FixAdaptFEP(LAMMPS *lmp, int narg, char **arg) :
if (adapt[m].which == PAIR)
memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig");
// when adapting charge and using kspace,
// need to recompute additional params in kspace->setup()
if (chgflag && force->kspace) force->kspace->qsum_update_flag = 1;
id_fix_diam = id_fix_chg = NULL;
}
@ -188,7 +193,7 @@ FixAdaptFEP::~FixAdaptFEP()
}
delete [] adapt;
if (force->kspace) force->kspace->qsum_update_flag = 0;
if (chgflag && force->kspace) force->kspace->qsum_update_flag = 0;
// check nfix in case all fixes have already been deleted
@ -291,11 +296,6 @@ void FixAdaptFEP::init()
if (adapt[i].which == ATOM)
error->all(FLERR,"Cannot use dynamic group with fix adapt/fep atom");
// when using kspace, we need to recompute
// some additional parameters in kspace->setup()
if (force->kspace) force->kspace->qsum_update_flag = 1;
// setup and error checks
anypair = 0;
@ -455,6 +455,8 @@ void FixAdaptFEP::change_settings()
} else if (ad->which == KSPACE) {
*kspace_scale = value;
// set per atom values, also make changes for ghost atoms
} else if (ad->which == ATOM) {
// reset radius from diameter
@ -510,13 +512,14 @@ void FixAdaptFEP::change_settings()
if (anypair) force->pair->reinit();
// re-setup KSpace if using it, since charges may have changed
// re-setup KSpace if it exists and adapting charges
// since charges have changed
if (force->kspace) force->kspace->setup();
if (chgflag && force->kspace) force->kspace->setup();
}
/* ----------------------------------------------------------------------
restore pair,kspace,atom parameters to original values
restore pair,kspace.atom parameters to original values
------------------------------------------------------------------------- */
void FixAdaptFEP::restore_settings()
@ -565,5 +568,5 @@ void FixAdaptFEP::restore_settings()
}
if (anypair) force->pair->reinit();
if (force->kspace) force->kspace->setup();
if (chgflag && force->kspace) force->kspace->setup();
}

1
src/VORONOI/compute_voronoi_atom.cpp
View File

@ -15,6 +15,7 @@
Contributing author: Daniel Schwen
------------------------------------------------------------------------- */
#include "mpi.h"
#include "math.h"
#include "string.h"
#include "stdlib.h"