diff --git a/src/MANYBODY/fix_qeq_comb.cpp b/src/MANYBODY/fix_qeq_comb.cpp index b7c87e3407..dcc9c15dbc 100644 --- a/src/MANYBODY/fix_qeq_comb.cpp +++ b/src/MANYBODY/fix_qeq_comb.cpp @@ -110,7 +110,7 @@ int FixQEQComb::setmask() void FixQEQComb::init() { if (!atom->q_flag) - error->all("Atoms must have charge to use fix qeq/comb"); + error->all("Fix qeq/comb requires atom attribute q"); comb = (PairComb *) force->pair_match("comb",1); if (comb == NULL) error->all("Must use pair_style comb with fix qeq/comb"); diff --git a/src/MANYBODY/pair_comb.cpp b/src/MANYBODY/pair_comb.cpp index 0f1b5f1a0c..9ed7d0c93a 100644 --- a/src/MANYBODY/pair_comb.cpp +++ b/src/MANYBODY/pair_comb.cpp @@ -14,7 +14,7 @@ /* ---------------------------------------------------------------------- Contributing author: Tzu-Ray Shan (U Florida, rayshan@ufl.edu) LAMMPS implementation of the Charge-optimized many-body (COMB) potential - based on HELL MD program (Prof Simon R. Phillpot, UF, sphil@mse.ufl.edu) + based on the HELL MD program (Prof Simon Phillpot, UF, sphil@mse.ufl.edu) and Aidan Thompson's Tersoff code in LAMMPS ------------------------------------------------------------------------- */ diff --git a/src/Makefile b/src/Makefile index 3d4843a1bd..aa39dff1c6 100755 --- a/src/Makefile +++ b/src/Makefile @@ -16,8 +16,8 @@ OBJ = $(SRC:.cpp=.o) PACKAGE = asphere class2 colloid dipole dsmc gpu granular \ kspace manybody meam molecule opt peri poems prd reax shock xtc -PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm user-ewaldn \ - user-imd user-smd +PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm \ + user-ewaldn user-imd user-smd PACKALL = $(PACKAGE) $(PACKUSER) diff --git a/src/atom.cpp b/src/atom.cpp index 00334ca3b5..49e95007cc 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -73,6 +73,9 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) vfrac = s0 = NULL; x0 = NULL; + spin = NULL; + eradius = evel = eforce = NULL; + maxspecial = 1; nspecial = NULL; special = NULL; @@ -98,6 +101,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) q_flag = mu_flag = 0; quat_flag = omega_flag = angmom_flag = torque_flag = 0; radius_flag = density_flag = rmass_flag = vfrac_flag = 0; + spin_flag = eradius_flag = evel_flag = eforce_flag = 0; // ntype-length arrays @@ -176,6 +180,11 @@ Atom::~Atom() memory->sfree(s0); memory->destroy_2d_double_array(x0); + memory->sfree(spin); + memory->sfree(eradius); + memory->sfree(evel); + memory->sfree(eforce); + memory->sfree(molecule); memory->destroy_2d_int_array(nspecial); diff --git a/src/atom.h b/src/atom.h index 39f567f883..87d75e3b45 100644 --- a/src/atom.h +++ b/src/atom.h @@ -52,6 +52,9 @@ class Atom : protected Pointers { double *radius,*density,*rmass,*vfrac,*s0; double **x0; + int *spin; + double *eradius,*evel,*eforce; + int **nspecial; // 0,1,2 = cummulative # of 1-2,1-3,1-4 neighs int **special; // IDs of 1-2,1-3,1-4 neighs of each atom int maxspecial; // special[nlocal][maxspecial] @@ -80,6 +83,7 @@ class Atom : protected Pointers { int q_flag,mu_flag; int quat_flag,omega_flag,angmom_flag,torque_flag; int radius_flag,density_flag,rmass_flag,vfrac_flag; + int spin_flag,eradius_flag,evel_flag,eforce_flag; // extra peratom info in restart file destined for fix & diag diff --git a/src/compute_centro_atom.h b/src/compute_centro_atom.h index 889c925d64..44553aec65 100644 --- a/src/compute_centro_atom.h +++ b/src/compute_centro_atom.h @@ -1,50 +1,50 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -#ifdef COMPUTE_CLASS - -ComputeStyle(centro/atom,ComputeCentroAtom) - -#else - -#ifndef COMPUTE_CENTRO_ATOM_H -#define COMPUTE_CENTRO_ATOM_H - -#include "compute.h" - -namespace LAMMPS_NS { - -class ComputeCentroAtom : public Compute { - public: - ComputeCentroAtom(class LAMMPS *, int, char **); - ~ComputeCentroAtom(); - void init(); - void init_list(int, class NeighList *); - void compute_peratom(); - double memory_usage(); - - private: - int nmax,maxneigh,nnn; - double *distsq; - int *nearest; - class NeighList *list; - double *centro; - - void select(int, int, double *); - void select2(int, int, double *, int *); -}; - -} - -#endif -#endif +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS + +ComputeStyle(centro/atom,ComputeCentroAtom) + +#else + +#ifndef COMPUTE_CENTRO_ATOM_H +#define COMPUTE_CENTRO_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeCentroAtom : public Compute { + public: + ComputeCentroAtom(class LAMMPS *, int, char **); + ~ComputeCentroAtom(); + void init(); + void init_list(int, class NeighList *); + void compute_peratom(); + double memory_usage(); + + private: + int nmax,maxneigh,nnn; + double *distsq; + int *nearest; + class NeighList *list; + double *centro; + + void select(int, int, double *); + void select2(int, int, double *, int *); +}; + +} + +#endif +#endif diff --git a/src/compute_temp.h b/src/compute_temp.h index 5bbe0ba7f4..3d3c030b67 100644 --- a/src/compute_temp.h +++ b/src/compute_temp.h @@ -27,16 +27,16 @@ namespace LAMMPS_NS { class ComputeTemp : public Compute { public: ComputeTemp(class LAMMPS *, int, char **); - ~ComputeTemp(); + virtual ~ComputeTemp(); void init(); double compute_scalar(); void compute_vector(); - private: + protected: int fix_dof; double tfactor; - void dof_compute(); + virtual void dof_compute(); }; } diff --git a/src/compute_temp_deform.h b/src/compute_temp_deform.h index 0b69b73c8f..ae066856b8 100644 --- a/src/compute_temp_deform.h +++ b/src/compute_temp_deform.h @@ -27,10 +27,10 @@ namespace LAMMPS_NS { class ComputeTempDeform : public Compute { public: ComputeTempDeform(class LAMMPS *, int, char **); - ~ComputeTempDeform(); + virtual ~ComputeTempDeform(); void init(); - double compute_scalar(); - void compute_vector(); + virtual double compute_scalar(); + virtual void compute_vector(); void remove_bias(int, double *); void remove_bias_all(); @@ -38,14 +38,14 @@ class ComputeTempDeform : public Compute { void restore_bias_all(); double memory_usage(); - private: + protected: int fix_dof; double tfactor; double vbias[3]; // stored velocity bias for one atom double **vbiasall; // stored velocity bias for all atoms int maxbias; // size of vbiasall array - void dof_compute(); + virtual void dof_compute(); }; } diff --git a/src/compute_temp_region.h b/src/compute_temp_region.h index 7c2f7e9b43..d886da5ceb 100644 --- a/src/compute_temp_region.h +++ b/src/compute_temp_region.h @@ -27,10 +27,10 @@ namespace LAMMPS_NS { class ComputeTempRegion : public Compute { public: ComputeTempRegion(class LAMMPS *, int, char **); - ~ComputeTempRegion(); + virtual ~ComputeTempRegion(); void init(); - double compute_scalar(); - void compute_vector(); + virtual double compute_scalar(); + virtual void compute_vector(); int dof_remove(int); void remove_bias(int, double *); @@ -39,7 +39,7 @@ class ComputeTempRegion : public Compute { void restore_bias_all(); double memory_usage(); - private: + protected: int iregion; char *idregion; }; diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 4e8783f199..20bbe1f9dc 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -11,6 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +#include "math.h" #include "stdlib.h" #include "string.h" #include "dump_custom.h" @@ -37,7 +38,7 @@ enum{ID,MOL,TYPE,MASS, X,Y,Z,XS,YS,ZS,XSTRI,YSTRI,ZSTRI,XU,YU,ZU,XUTRI,YUTRI,ZUTRI,IX,IY,IZ, VX,VY,VZ,FX,FY,FZ, Q,MUX,MUY,MUZ,RADIUS,OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ, - QUATW,QUATI,QUATJ,QUATK,TQX,TQY,TQZ, + QUATW,QUATI,QUATJ,QUATK,TQX,TQY,TQZ,SPIN,ERADIUS,EVEL,EFORCE, COMPUTE,FIX,VARIABLE}; enum{LT,LE,GT,GE,EQ,NEQ}; enum{INT,DOUBLE}; @@ -586,6 +587,7 @@ int DumpCustom::count() } else if (thresh_array[ithresh] == FZ) { ptr = &atom->f[0][2]; nstride = 3; + } else if (thresh_array[ithresh] == Q) { if (!atom->q_flag) error->all("Threshhold for an atom property that isn't allocated"); @@ -606,6 +608,7 @@ int DumpCustom::count() error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->mu[0][2]; nstride = 3; + } else if (thresh_array[ithresh] == RADIUS) { if (!atom->radius_flag) error->all("Threshhold for an atom property that isn't allocated"); @@ -641,6 +644,7 @@ int DumpCustom::count() error->all("Threshhold for an atom property that isn't allocated"); ptr = &atom->angmom[0][2]; nstride = 3; + } else if (thresh_array[ithresh] == QUATW) { if (!atom->quat_flag) error->all("Threshhold for an atom property that isn't allocated"); @@ -677,6 +681,29 @@ int DumpCustom::count() ptr = &atom->torque[0][2]; nstride = 3; + } else if (thresh_array[ithresh] == SPIN) { + if (!atom->spin_flag) + error->all("Threshhold for an atom property that isn't allocated"); + int *spin = atom->spin; + for (i = 0; i < nlocal; i++) dchoose[i] = spin[i]; + ptr = dchoose; + nstride = 1; + } else if (thresh_array[ithresh] == ERADIUS) { + if (!atom->eradius_flag) + error->all("Threshhold for an atom property that isn't allocated"); + ptr = atom->eradius; + nstride = 1; + } else if (thresh_array[ithresh] == EVEL) { + if (!atom->evel_flag) + error->all("Threshhold for an atom property that isn't allocated"); + ptr = atom->evel; + nstride = 1; + } else if (thresh_array[ithresh] == EFORCE) { + if (!atom->eforce_flag) + error->all("Threshhold for an atom property that isn't allocated"); + ptr = atom->eforce; + nstride = 1; + } else if (thresh_array[ithresh] == COMPUTE) { i = nfield + ithresh; if (argindex[i] == 0) { @@ -971,6 +998,27 @@ void DumpCustom::parse_fields(int narg, char **arg) pack_choice[i] = &DumpCustom::pack_tqz; vtype[i] = DOUBLE; + } else if (strcmp(arg[iarg],"spin") == 0) { + if (!atom->spin_flag) + error->all("Dumping an atom quantity that isn't allocated"); + pack_choice[i] = &DumpCustom::pack_spin; + vtype[i] = INT; + } else if (strcmp(arg[iarg],"eradius") == 0) { + if (!atom->eradius_flag) + error->all("Dumping an atom quantity that isn't allocated"); + pack_choice[i] = &DumpCustom::pack_eradius; + vtype[i] = DOUBLE; + } else if (strcmp(arg[iarg],"evel") == 0) { + if (!atom->evel_flag) + error->all("Dumping an atom quantity that isn't allocated"); + pack_choice[i] = &DumpCustom::pack_evel; + vtype[i] = DOUBLE; + } else if (strcmp(arg[iarg],"eforce") == 0) { + if (!atom->evel_flag) + error->all("Dumping an atom quantity that isn't allocated"); + pack_choice[i] = &DumpCustom::pack_eforce; + vtype[i] = DOUBLE; + // compute value = c_ID // if no trailing [], then arg is set to 0, else arg is int between [] @@ -1234,10 +1282,12 @@ int DumpCustom::modify_param(int narg, char **arg) else if (strcmp(arg[1],"fx") == 0) thresh_array[nthresh] = FX; else if (strcmp(arg[1],"fy") == 0) thresh_array[nthresh] = FY; else if (strcmp(arg[1],"fz") == 0) thresh_array[nthresh] = FZ; + else if (strcmp(arg[1],"q") == 0) thresh_array[nthresh] = Q; else if (strcmp(arg[1],"mux") == 0) thresh_array[nthresh] = MUX; else if (strcmp(arg[1],"muy") == 0) thresh_array[nthresh] = MUY; else if (strcmp(arg[1],"muz") == 0) thresh_array[nthresh] = MUZ; + else if (strcmp(arg[1],"radius") == 0) thresh_array[nthresh] = RADIUS; else if (strcmp(arg[1],"omegax") == 0) thresh_array[nthresh] = OMEGAX; else if (strcmp(arg[1],"omegay") == 0) thresh_array[nthresh] = OMEGAY; @@ -1245,6 +1295,7 @@ int DumpCustom::modify_param(int narg, char **arg) else if (strcmp(arg[1],"angmomx") == 0) thresh_array[nthresh] = ANGMOMX; else if (strcmp(arg[1],"angmomy") == 0) thresh_array[nthresh] = ANGMOMY; else if (strcmp(arg[1],"angmomz") == 0) thresh_array[nthresh] = ANGMOMZ; + else if (strcmp(arg[1],"quatw") == 0) thresh_array[nthresh] = QUATW; else if (strcmp(arg[1],"quati") == 0) thresh_array[nthresh] = QUATI; else if (strcmp(arg[1],"quatj") == 0) thresh_array[nthresh] = QUATJ; @@ -1252,7 +1303,12 @@ int DumpCustom::modify_param(int narg, char **arg) else if (strcmp(arg[1],"tqx") == 0) thresh_array[nthresh] = TQX; else if (strcmp(arg[1],"tqy") == 0) thresh_array[nthresh] = TQY; else if (strcmp(arg[1],"tqz") == 0) thresh_array[nthresh] = TQZ; - + + else if (strcmp(arg[1],"spin") == 0) thresh_array[nthresh] = SPIN; + else if (strcmp(arg[1],"eradius") == 0) thresh_array[nthresh] = ERADIUS; + else if (strcmp(arg[1],"evel") == 0) thresh_array[nthresh] = EVEL; + else if (strcmp(arg[1],"eforce") == 0) thresh_array[nthresh] = EFORCE; + // compute value = c_ID // if no trailing [], then arg is set to 0, else arg is between [] // must grow field2index and argindex arrays, since access is beyond nfield @@ -2175,3 +2231,86 @@ void DumpCustom::pack_tqz(int n) n += size_one; } } + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_spin(int n) +{ + int *spin = atom->spin; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (choose[i]) { + buf[n] = spin[i]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- + different interpretation of electron radius + depending on dynamics vs minimization +------------------------------------------------------------------------- */ + +void DumpCustom::pack_eradius(int n) +{ + double *eradius = atom->eradius; + int *spin = atom->spin; + int nlocal = atom->nlocal; + + if (update->whichflag == 1) { + for (int i = 0; i < nlocal; i++) + if (choose[i]) { + buf[n] = eradius[i]; + n += size_one; + } + } else { + for (int i = 0; i < nlocal; i++) + if (choose[i]) { + if (spin[i]) buf[n] = exp(eradius[i]); + else buf[n] = 0.0; + n += size_one; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_evel(int n) +{ + double *evel = atom->evel; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (choose[i]) { + buf[n] = evel[i]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- + different interpretation of electron radial force + depending on dynamics vs minimization +------------------------------------------------------------------------- */ + +void DumpCustom::pack_eforce(int n) +{ + double *eradius = atom->eradius; + double *eforce = atom->eforce; + int *spin = atom->spin; + int nlocal = atom->nlocal; + + if (update->whichflag == 1) { + for (int i = 0; i < nlocal; i++) + if (choose[i]) { + buf[n] = eforce[i]; + n += size_one; + } + } else { + for (int i = 0; i < nlocal; i++) + if (choose[i]) { + if (spin[i]) buf[n] = eforce[i]/exp(eradius[i]); + else buf[n] = 0.0; + n += size_one; + } + } +} diff --git a/src/dump_custom.h b/src/dump_custom.h index 4cdea81412..03828adb63 100644 --- a/src/dump_custom.h +++ b/src/dump_custom.h @@ -148,6 +148,10 @@ class DumpCustom : public Dump { void pack_tqx(int); void pack_tqy(int); void pack_tqz(int); + void pack_spin(int); + void pack_eradius(int); + void pack_evel(int); + void pack_eforce(int); void pack_compute(int); void pack_fix(int); diff --git a/src/min.cpp b/src/min.cpp index 4fb96e7744..bc4bfc3c79 100644 --- a/src/min.cpp +++ b/src/min.cpp @@ -135,10 +135,12 @@ void Min::init() ev_setup(); // set flags for what arrays to clear in force_clear() - // clear torques if array exists + // need to clear torques,eforce if arrays exists torqueflag = 0; - if (atom->torque) torqueflag = 1; + if (atom->torque_flag) torqueflag = 1; + eforceflag = 0; + if (atom->eforce_flag) eforceflag = 1; // orthogonal vs triclinic simulation box @@ -496,6 +498,8 @@ double Min::energy_force(int resetflag) void Min::force_clear() { + int i; + // clear global force array // nall includes ghosts only if either newton flag is set @@ -504,7 +508,7 @@ void Min::force_clear() else nall = atom->nlocal; double **f = atom->f; - for (int i = 0; i < nall; i++) { + for (i = 0; i < nall; i++) { f[i][0] = 0.0; f[i][1] = 0.0; f[i][2] = 0.0; @@ -512,12 +516,18 @@ void Min::force_clear() if (torqueflag) { double **torque = atom->torque; - for (int i = 0; i < nall; i++) { + for (i = 0; i < nall; i++) { torque[i][0] = 0.0; torque[i][1] = 0.0; torque[i][2] = 0.0; } } + + if (eforceflag) { + double *eforce = atom->eforce; + for (i = 0; i < nall; i++) + eforce[i] = 0.0; + } } /* ---------------------------------------------------------------------- diff --git a/src/min.h b/src/min.h index dcea3d85c6..53457acc85 100644 --- a/src/min.h +++ b/src/min.h @@ -57,7 +57,8 @@ class Min : protected Pointers { class Compute **vlist_global; class Compute **vlist_atom; - int pairflag,torqueflag; + int pairflag; + int torqueflag,eforceflag; int triclinic; // 0 if domain is orthog, 1 if triclinic int narray; // # of arrays stored by fix_minimize diff --git a/src/min_hftn.cpp b/src/min_hftn.cpp index 81b1cdb5b7..d4322ef06a 100644 --- a/src/min_hftn.cpp +++ b/src/min_hftn.cpp @@ -180,7 +180,7 @@ void MinHFTN::reset_vectors() int n = NUM_HFTN_ATOM_BASED_VECTORS; for (int m = 0; m < nextra_atom; m++) { extra_nlen[m] = extra_peratom[m] * atom->nlocal; - requestor[m]->min_pointers (&xextra_atom[m], &fextra_atom[m]); + requestor[m]->min_pointers(&xextra_atom[m], &fextra_atom[m]); for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) _daExtraAtom[i][m] = fix_minimize->request_vector (n++); } diff --git a/src/update.cpp b/src/update.cpp index 502aa869f2..f1a9fd2cdb 100644 --- a/src/update.cpp +++ b/src/update.cpp @@ -166,6 +166,19 @@ void Update::set_units(const char *style) dt = 1.0e-8; neighbor->skin = 0.1; + } else if (strcmp(style,"electron") == 0) { + force->boltz = 3.16679e-6; + force->mvv2e = 1.0; + force->ftm2v = 1 / 1.0327499; + force->mv2d = 1.0; + force->nktv2p = 2.9421e13; + force->qqr2e = 1.0; + force->qe2f = 1.0; + force->vxmu2f = 0.6241509647; + force->xxt2kmu = 1.0e-4; + dt = 0.001; + neighbor->skin = 2.0; + } else error->all("Illegal units command"); delete [] unit_style; diff --git a/src/verlet.cpp b/src/verlet.cpp index c454bdd60b..523a03839a 100644 --- a/src/verlet.cpp +++ b/src/verlet.cpp @@ -63,10 +63,12 @@ void Verlet::init() ev_setup(); // set flags for what arrays to clear in force_clear() - // need to clear torques if array exists + // need to clear torques,eforce if arrays exists torqueflag = 0; if (atom->torque_flag) torqueflag = 1; + eforceflag = 0; + if (atom->eforce_flag) eforceflag = 1; // orthogonal vs triclinic simulation box @@ -317,6 +319,12 @@ void Verlet::force_clear() } } + if (eforceflag) { + double *eforce = atom->eforce; + for (i = 0; i < nall; i++) + eforce[i] = 0.0; + } + // neighbor includegroup flag is set // clear force only on initial nfirst particles // if either newton flag is set, also include ghosts @@ -340,6 +348,12 @@ void Verlet::force_clear() } } + if (eforceflag) { + double *eforce = atom->eforce; + for (i = 0; i < nall; i++) + eforce[i] = 0.0; + } + if (force->newton) { nall = atom->nlocal + atom->nghost; @@ -357,6 +371,12 @@ void Verlet::force_clear() torque[i][2] = 0.0; } } + + if (eforceflag) { + double *eforce = atom->eforce; + for (i = atom->nlocal; i < nall; i++) + eforce[i] = 0.0; + } } } } diff --git a/src/verlet.h b/src/verlet.h index de57b4dbc0..9257691640 100644 --- a/src/verlet.h +++ b/src/verlet.h @@ -35,7 +35,8 @@ class Verlet : public Integrate { private: int triclinic; // 0 if domain is orthog, 1 if triclinic - int torqueflag; // zero out array every step + int torqueflag; // zero out arrays every step + int eforceflag; void force_clear(); };