diff --git a/src/USER-AWPMD/Install.sh b/src/USER-AWPMD/Install.sh
new file mode 100644
index 000000000..a6e1d7a99
--- /dev/null
+++ b/src/USER-AWPMD/Install.sh
@@ -0,0 +1,37 @@
+# Install/unInstall package files in LAMMPS
+# edit Makefile.package to include/exclude ATC library
+
+if (test $1 = 1) then
+
+ if (test -e ../Makefile.package) then
+ sed -i -e 's/[^ \t]*awpmd[^ \t]* //g' ../Makefile.package
+ sed -i -e 's|^PKG_INC =[ \t]*|&-I../../lib/awpmd/ivutils/include -I../../lib/awpmd/systems/interact |' ../Makefile.package
+ sed -i -e 's|^PKG_PATH =[ \t]*|&-L../../lib/awpmd |' ../Makefile.package
+ sed -i -e 's|^PKG_LIB =[ \t]*|&-lawpmd |' ../Makefile.package
+ sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(user-awpmd_SYSPATH) |' ../Makefile.package
+ sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(user-awpmd_SYSLIB) |' ../Makefile.package
+ fi
+
+ cp atom_vec_wavepacket.cpp ..
+ cp fix_nve_awpmd.cpp ..
+ cp pair_awpmd_cut.cpp ..
+
+ cp atom_vec_wavepacket.h ..
+ cp fix_nve_awpmd.h ..
+ cp pair_awpmd_cut.h ..
+
+elif (test $1 = 0) then
+
+ if (test -e ../Makefile.package) then
+ sed -i -e 's/[^ \t]*awpmd[^ \t]* //g' ../Makefile.package
+ fi
+
+ rm -f ../atom_vec_wavepacket.cpp
+ rm -f ../fix_nve_awpmd.cpp
+ rm -f ../pair_awpmd_cut.cpp
+
+ rm -f ../atom_vec_wavepacket.h
+ rm -f ../fix_nve_awpmd.h
+ rm -f ../pair_awpmd_cut.h
+
+fi
diff --git a/src/USER-AWPMD/README b/src/USER-AWPMD/README
new file mode 100644
index 000000000..66378cd45
--- /dev/null
+++ b/src/USER-AWPMD/README
@@ -0,0 +1,24 @@
+The files in this directory are a user-contributed package for LAMMPS.
+
+The person who created these files is Ilya Valuev
+(valuev@physik.hu-berlin.de). Contact him directly if you have
+questions.
+
+PACKAGE DESCRIPTION:
+
+Contains a LAMMPS implementation of the Antisymmetrized Wave Packet
+Molecular Dynamics (AWPMD) method.
+
+INSTALLATION:
+
+- compile the awpmd library in lib/awpmd
+- follow a normal LAMMPS package installation: make yes-user-awpmd
+
+OTHERS FILES INCLUDED:
+
+User examples are under examples/USER/awpmd
+
+ACKNOWLEDGMENTS:
+
+Thanks to Steve Plimpton and Aidan Thompson for their help in
+adaptation of the AWPMD code to LAMMPS.
diff --git a/src/USER-AWPMD/atom_vec_wavepacket.cpp b/src/USER-AWPMD/atom_vec_wavepacket.cpp
new file mode 100644
index 000000000..a9edfe4e4
--- /dev/null
+++ b/src/USER-AWPMD/atom_vec_wavepacket.cpp
@@ -0,0 +1,1009 @@
+/* ----------------------------------------------------------------------
+ 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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Ilya Valuev (JIHT RAS)
+------------------------------------------------------------------------- */
+
+
+#include "math.h"
+#include "stdlib.h"
+#include "atom_vec_wavepacket.h"
+#include "atom.h"
+#include "domain.h"
+#include "modify.h"
+#include "force.h"
+#include "fix.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define DELTA 10000
+
+
+/* ---------------------------------------------------------------------- */
+
+AtomVecWavepacket::AtomVecWavepacket(LAMMPS *lmp, int narg, char **arg) :
+ AtomVec(lmp, narg, arg)
+{
+ comm_x_only = comm_f_only = 0;
+
+ mass_type = 1;
+ molecular = 0;
+
+ size_forward = 4; // coords[3]+radius[1]
+ size_reverse = 10; // force[3]+erforce[1]+ervelforce[1]+vforce[3]+csforce[2]
+ size_border = 10; // coords[3]+tag[1]+type[1]+mask[1]+q[1]+spin[1]+eradius[1]+etag[1]
+ size_velocity = 6; // +velocities[3]+ ervel[1]+cs[2]
+ size_data_atom = 11; // for input file: 1-tag 2-type 3-q 4-spin 5-eradius 6-etag 7-cs_re 8-cs_im 9-x 10-y 11-z
+ size_data_vel = 5; // for input file: vx vy vz ervel <??>
+ xcol_data = 9; // starting column for x data
+
+ atom->wavepacket_flag = 1;
+ atom->electron_flag = 1; // compatible with eff
+ atom->q_flag = atom->spin_flag = atom->eradius_flag =
+ atom->ervel_flag = atom->erforce_flag = 1;
+
+ atom->cs_flag = atom->csforce_flag = atom->vforce_flag = atom->ervelforce_flag = atom->etag_flag = 1;
+
+
+}
+
+/* ----------------------------------------------------------------------
+ grow atom-electron arrays
+ n = 0 grows arrays by DELTA
+ n > 0 allocates arrays to size n
+------------------------------------------------------------------------- */
+
+void AtomVecWavepacket::grow(int n)
+{
+ if (n == 0) nmax += DELTA;
+ else nmax = n;
+ atom->nmax = nmax;
+
+ tag = memory->grow(atom->tag,nmax,"atom:tag");
+ type = memory->grow(atom->type,nmax,"atom:type");
+ mask = memory->grow(atom->mask,nmax,"atom:mask");
+ image = memory->grow(atom->image,nmax,"atom:image");
+ x = memory->grow(atom->x,nmax,3,"atom:x");
+ v = memory->grow(atom->v,nmax,3,"atom:v");
+ f = memory->grow(atom->f,nmax,3,"atom:f");
+
+ q = memory->grow(atom->q,nmax,"atom:q");
+ spin = memory->grow(atom->spin,nmax,"atom:spin");
+ eradius = memory->grow(atom->eradius,nmax,"atom:eradius");
+ ervel = memory->grow(atom->ervel,nmax,"atom:ervel");
+ erforce = memory->grow(atom->erforce,nmax,"atom:erforce");
+
+ cs = memory->grow(atom->cs,2*nmax,"atom:cs");
+ csforce = memory->grow(atom->csforce,2*nmax,"atom:csforce");
+ vforce = memory->grow(atom->vforce,3*nmax,"atom:vforce");
+ ervelforce = memory->grow(atom->ervelforce,nmax,"atom:ervelforce");
+ etag = memory->grow(atom->etag,nmax,"atom:etag");
+
+ if (atom->nextra_grow)
+ for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
+ modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax);
+}
+
+/* ----------------------------------------------------------------------
+ reset local array ptrs
+------------------------------------------------------------------------- */
+
+void AtomVecWavepacket::grow_reset()
+{
+ tag = atom->tag; type = atom->type;
+ mask = atom->mask; image = atom->image;
+ x = atom->x; v = atom->v; f = atom->f;
+ q = atom->q;
+ eradius = atom->eradius; ervel = atom->ervel; erforce = atom->erforce;
+
+ cs = atom->cs;
+ csforce = atom->csforce;
+ vforce = atom->vforce;
+ ervelforce = atom->ervelforce;
+ etag = atom->etag;
+
+}
+
+/* ----------------------------------------------------------------------
+ copy atom I info to atom J
+------------------------------------------------------------------------- */
+
+void AtomVecWavepacket::copy(int i, int j, int delflag)
+{
+ tag[j] = tag[i];
+ type[j] = type[i];
+ mask[j] = mask[i];
+ image[j] = image[i];
+ x[j][0] = x[i][0];
+ x[j][1] = x[i][1];
+ x[j][2] = x[i][2];
+ v[j][0] = v[i][0];
+ v[j][1] = v[i][1];
+ v[j][2] = v[i][2];
+
+ q[j] = q[i];
+ spin[j] = spin[i];
+ eradius[j] = eradius[i];
+ ervel[j] = ervel[i];
+
+ cs[2*j] = cs[2*i];
+ cs[2*j+1] = cs[2*i+1];
+ etag[j] = etag[i];
+
+
+ if (atom->nextra_grow)
+ for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
+ modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j);
+}
+
+/* ---------------------------------------------------------------------- */
+// this will be used as partial pack for unsplit Hartree packets (v, ervel not regarded as separate variables)
+
+int AtomVecWavepacket::pack_comm(int n, int *list, double *buf,
+ int pbc_flag, int *pbc)
+{
+ int i,j,m;
+ double dx,dy,dz;
+
+ m = 0;
+ if (pbc_flag == 0) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0];
+ buf[m++] = x[j][1];
+ buf[m++] = x[j][2];
+ buf[m++] = eradius[j];
+ }
+ } else {
+ if (domain->triclinic == 0) {
+ dx = pbc[0]*domain->xprd;
+ dy = pbc[1]*domain->yprd;
+ dz = pbc[2]*domain->zprd;
+ } else {
+ dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
+ dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
+ dz = pbc[2]*domain->zprd;
+ }
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0] + dx;
+ buf[m++] = x[j][1] + dy;
+ buf[m++] = x[j][2] + dz;
+ buf[m++] = eradius[j];
+ }
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+// this is a complete pack of all 'position' variables of AWPMD
+
+int AtomVecWavepacket::pack_comm_vel(int n, int *list, double *buf,
+ int pbc_flag, int *pbc)
+{
+ int i,j,m;
+ double dx,dy,dz,dvx,dvy,dvz;
+
+ m = 0;
+ if (pbc_flag == 0) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0];
+ buf[m++] = x[j][1];
+ buf[m++] = x[j][2];
+ buf[m++] = eradius[j];
+ buf[m++] = v[j][0];
+ buf[m++] = v[j][1];
+ buf[m++] = v[j][2];
+
+ buf[m++] = ervel[j];
+ buf[m++] = cs[2*j];
+ buf[m++] = cs[2*j+1];
+ }
+ } else {
+ if (domain->triclinic == 0) {
+ dx = pbc[0]*domain->xprd;
+ dy = pbc[1]*domain->yprd;
+ dz = pbc[2]*domain->zprd;
+ } else {
+ dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
+ dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
+ dz = pbc[2]*domain->zprd;
+ }
+ if (!deform_vremap) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0] + dx;
+ buf[m++] = x[j][1] + dy;
+ buf[m++] = x[j][2] + dz;
+ buf[m++] = eradius[j];
+ buf[m++] = v[j][0];
+ buf[m++] = v[j][1];
+ buf[m++] = v[j][2];
+
+ buf[m++] = ervel[j];
+ buf[m++] = cs[2*j];
+ buf[m++] = cs[2*j+1];
+ }
+ } else {
+ dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
+ dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
+ dvz = pbc[2]*h_rate[2];
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0] + dx;
+ buf[m++] = x[j][1] + dy;
+ buf[m++] = x[j][2] + dz;
+ buf[m++] = eradius[j];
+ if (mask[i] & deform_groupbit) {
+ buf[m++] = v[j][0] + dvx;
+ buf[m++] = v[j][1] + dvy;
+ buf[m++] = v[j][2] + dvz;
+ } else {
+ buf[m++] = v[j][0];
+ buf[m++] = v[j][1];
+ buf[m++] = v[j][2];
+ }
+ buf[m++] = ervel[j];
+ buf[m++] = cs[2*j];
+ buf[m++] = cs[2*j+1];
+ }
+ }
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecWavepacket::pack_comm_hybrid(int n, int *list, double *buf)
+{
+ int i,j,m;
+
+ m = 0;
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = eradius[j];
+ buf[m++] = ervel[j];
+ buf[m++] = cs[2*j];
+ buf[m++] = cs[2*j+1];
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void AtomVecWavepacket::unpack_comm(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) {
+ x[i][0] = buf[m++];
+ x[i][1] = buf[m++];
+ x[i][2] = buf[m++];
+ eradius[i] = buf[m++];
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void AtomVecWavepacket::unpack_comm_vel(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) {
+ x[i][0] = buf[m++];
+ x[i][1] = buf[m++];
+ x[i][2] = buf[m++];
+ eradius[i] = buf[m++];
+ v[i][0] = buf[m++];
+ v[i][1] = buf[m++];
+ v[i][2] = buf[m++];
+
+ ervel[i] = buf[m++];
+ cs[2*i] = buf[m++];
+ cs[2*i+1] = buf[m++];
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecWavepacket::unpack_comm_hybrid(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++){
+ eradius[i] = buf[m++];
+ ervel[i] = buf[m++];
+ cs[2*i] = buf[m++];
+ cs[2*i+1] = buf[m++];
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecWavepacket::pack_reverse(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) { //10
+ buf[m++] = f[i][0];
+ buf[m++] = f[i][1];
+ buf[m++] = f[i][2];
+ buf[m++] = erforce[i];
+
+ buf[m++] = ervelforce[i];
+ buf[m++] = vforce[3*i];
+ buf[m++] = vforce[3*i+1];
+ buf[m++] = vforce[3*i+2];
+ buf[m++] = csforce[2*i];
+ buf[m++] = csforce[2*i+1];
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecWavepacket::pack_reverse_hybrid(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++){
+ buf[m++] = erforce[i];
+
+ buf[m++] = ervelforce[i];
+ buf[m++] = vforce[3*i];
+ buf[m++] = vforce[3*i+1];
+ buf[m++] = vforce[3*i+2];
+ buf[m++] = csforce[2*i];
+ buf[m++] = csforce[2*i+1];
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void AtomVecWavepacket::unpack_reverse(int n, int *list, double *buf)
+{
+ int i,j,m;
+
+ m = 0;
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ f[j][0] += buf[m++];
+ f[j][1] += buf[m++];
+ f[j][2] += buf[m++];
+ erforce[j] += buf[m++];
+
+ ervelforce[j] += buf[m++];
+ vforce[3*j] += buf[m++];
+ vforce[3*j+1] += buf[m++];
+ vforce[3*j+2] += buf[m++];
+ csforce[2*j] += buf[m++];
+ csforce[2*j+1] += buf[m++];
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecWavepacket::unpack_reverse_hybrid(int n, int *list, double *buf)
+{
+ int i,j,m;
+
+ m = 0;
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ erforce[j] += buf[m++];
+
+ ervelforce[j] += buf[m++];
+ vforce[3*j] += buf[m++];
+ vforce[3*j+1] += buf[m++];
+ vforce[3*j+2] += buf[m++];
+ csforce[2*j] += buf[m++];
+ csforce[2*j+1] += buf[m++];
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+// will be used for Hartree unsplit version (the etag is added however)
+int AtomVecWavepacket::pack_border(int n, int *list, double *buf,
+ int pbc_flag, int *pbc)
+{
+ int i,j,m;
+ double dx,dy,dz;
+
+ m = 0;
+ if (pbc_flag == 0) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0];
+ buf[m++] = x[j][1];
+ buf[m++] = x[j][2];
+ buf[m++] = tag[j];
+ buf[m++] = type[j];
+ buf[m++] = mask[j];
+ buf[m++] = q[j];
+ buf[m++] = spin[j];
+ buf[m++] = eradius[j];
+ buf[m++] = etag[j];
+ }
+ } else {
+ if (domain->triclinic == 0) {
+ dx = pbc[0]*domain->xprd;
+ dy = pbc[1]*domain->yprd;
+ dz = pbc[2]*domain->zprd;
+ } else {
+ dx = pbc[0];
+ dy = pbc[1];
+ dz = pbc[2];
+ }
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0] + dx;
+ buf[m++] = x[j][1] + dy;
+ buf[m++] = x[j][2] + dz;
+ buf[m++] = tag[j];
+ buf[m++] = type[j];
+ buf[m++] = mask[j];
+ buf[m++] = q[j];
+ buf[m++] = spin[j];
+ buf[m++] = eradius[j];
+ buf[m++] = etag[j];
+ }
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecWavepacket::pack_border_vel(int n, int *list, double *buf,
+ int pbc_flag, int *pbc)
+{
+ int i,j,m;
+ double dx,dy,dz,dvx,dvy,dvz;
+
+ m = 0;
+ if (pbc_flag == 0) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0];
+ buf[m++] = x[j][1];
+ buf[m++] = x[j][2];
+ buf[m++] = tag[j];
+ buf[m++] = type[j];
+ buf[m++] = mask[j];
+ buf[m++] = q[j];
+ buf[m++] = spin[j];
+ buf[m++] = eradius[j];
+ buf[m++] = etag[j];
+
+ buf[m++] = v[j][0];
+ buf[m++] = v[j][1];
+ buf[m++] = v[j][2];
+
+
+ buf[m++] = ervel[j];
+ buf[m++] = cs[2*j];
+ buf[m++] = cs[2*j+1];
+ }
+ } else {
+ if (domain->triclinic == 0) {
+ dx = pbc[0]*domain->xprd;
+ dy = pbc[1]*domain->yprd;
+ dz = pbc[2]*domain->zprd;
+ } else {
+ dx = pbc[0];
+ dy = pbc[1];
+ dz = pbc[2];
+ }
+ if (domain->triclinic == 0) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0] + dx;
+ buf[m++] = x[j][1] + dy;
+ buf[m++] = x[j][2] + dz;
+ buf[m++] = tag[j];
+ buf[m++] = type[j];
+ buf[m++] = mask[j];
+ buf[m++] = q[j];
+ buf[m++] = spin[j];
+ buf[m++] = eradius[j];
+ buf[m++] = etag[j];
+
+ buf[m++] = v[j][0];
+ buf[m++] = v[j][1];
+ buf[m++] = v[j][2];
+
+
+ buf[m++] = ervel[j];
+ buf[m++] = cs[2*j];
+ buf[m++] = cs[2*j+1];
+ }
+ } else {
+ dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
+ dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
+ dvz = pbc[2]*h_rate[2];
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = x[j][0] + dx;
+ buf[m++] = x[j][1] + dy;
+ buf[m++] = x[j][2] + dz;
+ buf[m++] = tag[j];
+ buf[m++] = type[j];
+ buf[m++] = mask[j];
+ buf[m++] = q[j];
+ buf[m++] = spin[j];
+ buf[m++] = eradius[j];
+ buf[m++] = etag[j];
+
+ if (mask[i] & deform_groupbit) {
+ buf[m++] = v[j][0] + dvx;
+ buf[m++] = v[j][1] + dvy;
+ buf[m++] = v[j][2] + dvz;
+ } else {
+ buf[m++] = v[j][0];
+ buf[m++] = v[j][1];
+ buf[m++] = v[j][2];
+ }
+
+ buf[m++] = ervel[j];
+ buf[m++] = cs[2*j];
+ buf[m++] = cs[2*j+1];
+ }
+ }
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecWavepacket::pack_border_hybrid(int n, int *list, double *buf)
+{
+ int i,j,m;
+
+ m = 0;
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = q[j];
+ buf[m++] = spin[j];
+ buf[m++] = eradius[j];
+
+ buf[m++] = etag[j];
+ buf[m++] = ervel[j];
+ buf[m++] = cs[2*j];
+ buf[m++] = cs[2*j+1];
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void AtomVecWavepacket::unpack_border(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) {
+ if (i == nmax) grow(0);
+ x[i][0] = buf[m++];
+ x[i][1] = buf[m++];
+ x[i][2] = buf[m++];
+ tag[i] = static_cast<int> (buf[m++]);
+ type[i] = static_cast<int> (buf[m++]);
+ mask[i] = static_cast<int> (buf[m++]);
+ q[i] = buf[m++];
+ spin[i] = static_cast<int> (buf[m++]);
+ eradius[i] = buf[m++];
+ etag[i] = (int)buf[m++];
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void AtomVecWavepacket::unpack_border_vel(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) {
+ if (i == nmax) grow(0);
+ x[i][0] = buf[m++];
+ x[i][1] = buf[m++];
+ x[i][2] = buf[m++];
+ tag[i] = static_cast<int> (buf[m++]);
+ type[i] = static_cast<int> (buf[m++]);
+ mask[i] = static_cast<int> (buf[m++]);
+ q[i] = buf[m++];
+ spin[i] = static_cast<int> (buf[m++]);
+ eradius[i] = buf[m++];
+ etag[i] = (int)buf[m++];
+
+ v[i][0] = buf[m++];
+ v[i][1] = buf[m++];
+ v[i][2] = buf[m++];
+
+
+ ervel[i] = buf[m++];
+ cs[2*i] = buf[m++];
+ cs[2*i+1] = buf[m++];
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecWavepacket::unpack_border_hybrid(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) {
+ q[i] = buf[m++];
+ spin[i] = static_cast<int> (buf[m++]);
+ eradius[i] = buf[m++];
+
+ etag[i] = (int)buf[m++];
+ ervel[i] = buf[m++];
+ cs[2*i] = buf[m++];
+ cs[2*i+1] = buf[m++];
+ }
+ return m;
+}
+
+/* ----------------------------------------------------------------------
+ pack data for atom I for sending to another proc
+ xyz must be 1st 3 values, so comm::exchange() can test on them
+------------------------------------------------------------------------- */
+
+int AtomVecWavepacket::pack_exchange(int i, double *buf)
+{
+ int m = 1;
+ buf[m++] = x[i][0];
+ buf[m++] = x[i][1];
+ buf[m++] = x[i][2];
+ buf[m++] = v[i][0];
+ buf[m++] = v[i][1];
+ buf[m++] = v[i][2];
+ buf[m++] = tag[i];
+ buf[m++] = type[i];
+ buf[m++] = mask[i];
+ buf[m++] = image[i];
+ buf[m++] = q[i];
+ buf[m++] = spin[i];
+ buf[m++] = eradius[i];
+ buf[m++] = ervel[i];
+
+ buf[m++] = etag[i];
+ buf[m++] = cs[2*i];
+ buf[m++] = cs[2*i+1];
+
+ if (atom->nextra_grow)
+ for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
+ m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]);
+
+ buf[0] = m;
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecWavepacket::unpack_exchange(double *buf)
+{
+ int nlocal = atom->nlocal;
+ if (nlocal == nmax) grow(0);
+
+ int m = 1;
+ x[nlocal][0] = buf[m++];
+ x[nlocal][1] = buf[m++];
+ x[nlocal][2] = buf[m++];
+ v[nlocal][0] = buf[m++];
+ v[nlocal][1] = buf[m++];
+ v[nlocal][2] = buf[m++];
+ tag[nlocal] = static_cast<int> (buf[m++]);
+ type[nlocal] = static_cast<int> (buf[m++]);
+ mask[nlocal] = static_cast<int> (buf[m++]);
+ image[nlocal] = static_cast<int> (buf[m++]);
+ q[nlocal] = buf[m++];
+ spin[nlocal] = static_cast<int> (buf[m++]);
+ eradius[nlocal] = buf[m++];
+ ervel[nlocal] = buf[m++];
+
+ etag[nlocal] = buf[m++];
+ cs[2*nlocal] = buf[m++];
+ cs[2*nlocal+1] = buf[m++];
+
+ if (atom->nextra_grow)
+ for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
+ m += modify->fix[atom->extra_grow[iextra]]->
+ unpack_exchange(nlocal,&buf[m]);
+
+ atom->nlocal++;
+ return m;
+}
+
+/* ----------------------------------------------------------------------
+ size of restart data for all atoms owned by this proc
+ include extra data stored by fixes
+------------------------------------------------------------------------- */
+
+int AtomVecWavepacket::size_restart()
+{
+ int i;
+
+ int nlocal = atom->nlocal;
+ int n = 18 * nlocal; // Associated with pack_restart
+
+ if (atom->nextra_restart)
+ for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
+ for (i = 0; i < nlocal; i++)
+ n += modify->fix[atom->extra_restart[iextra]]->size_restart(i);
+
+ return n;
+}
+
+/* ----------------------------------------------------------------------
+ pack atom I's data for restart file including extra quantities
+ xyz must be 1st 3 values, so that read_restart can test on them
+ molecular types may be negative, but write as positive
+------------------------------------------------------------------------- */
+
+int AtomVecWavepacket::pack_restart(int i, double *buf)
+{
+ int m = 1;
+ buf[m++] = x[i][0];
+ buf[m++] = x[i][1];
+ buf[m++] = x[i][2];
+ buf[m++] = tag[i];
+ buf[m++] = type[i];
+ buf[m++] = mask[i];
+ buf[m++] = image[i];
+ buf[m++] = v[i][0];
+ buf[m++] = v[i][1];
+ buf[m++] = v[i][2];
+
+ buf[m++] = q[i];
+ buf[m++] = spin[i];
+ buf[m++] = eradius[i];
+ buf[m++] = ervel[i];
+
+ buf[m++] = etag[i];
+ buf[m++] = cs[2*i];
+ buf[m++] = cs[2*i+1];
+
+ if (atom->nextra_restart)
+ for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
+ m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]);
+
+ buf[0] = m;
+ return m;
+}
+
+/* ----------------------------------------------------------------------
+ unpack data for one atom from restart file including extra quantities
+------------------------------------------------------------------------- */
+
+int AtomVecWavepacket::unpack_restart(double *buf)
+{
+ int nlocal = atom->nlocal;
+ if (nlocal == nmax) {
+ grow(0);
+ if (atom->nextra_store)
+ memory->grow(atom->extra,nmax,atom->nextra_store,"atom:extra");
+ }
+
+ int m = 1;
+ x[nlocal][0] = buf[m++];
+ x[nlocal][1] = buf[m++];
+ x[nlocal][2] = buf[m++];
+ tag[nlocal] = static_cast<int> (buf[m++]);
+ type[nlocal] = static_cast<int> (buf[m++]);
+ mask[nlocal] = static_cast<int> (buf[m++]);
+ image[nlocal] = static_cast<int> (buf[m++]);
+ v[nlocal][0] = buf[m++];
+ v[nlocal][1] = buf[m++];
+ v[nlocal][2] = buf[m++];
+
+ q[nlocal] = buf[m++];
+ spin[nlocal] = static_cast<int> (buf[m++]);
+ eradius[nlocal] = buf[m++];
+ ervel[nlocal] = buf[m++];
+
+ etag[nlocal] = buf[m++];
+ cs[2*nlocal] = buf[m++];
+ cs[2*nlocal+1] = buf[m++];
+
+ double **extra = atom->extra;
+ if (atom->nextra_store) {
+ int size = static_cast<int> (buf[0]) - m;
+ for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++];
+ }
+
+ atom->nlocal++;
+ return m;
+}
+
+/* ----------------------------------------------------------------------
+ create one atom of itype at coord
+ set other values to defaults
+ AWPMD: creates a proton
+------------------------------------------------------------------------- */
+
+void AtomVecWavepacket::create_atom(int itype, double *coord)
+{
+ int nlocal = atom->nlocal;
+ if (nlocal == nmax) grow(0);
+
+ tag[nlocal] = 0;
+ type[nlocal] = itype;
+ x[nlocal][0] = coord[0];
+ x[nlocal][1] = coord[1];
+ x[nlocal][2] = coord[2];
+ mask[nlocal] = 1;
+ image[nlocal] = (512 << 20) | (512 << 10) | 512;
+ v[nlocal][0] = 0.0;
+ v[nlocal][1] = 0.0;
+ v[nlocal][2] = 0.0;
+
+ q[nlocal] = 1.;
+ spin[nlocal] = 0.;
+ eradius[nlocal] = 0.0;
+ ervel[nlocal] = 0.0;
+
+ etag[nlocal]= 0.;
+ cs[2*nlocal] = 0.;
+ cs[2*nlocal+1] = 0.;
+
+ atom->nlocal++;
+}
+
+/* ----------------------------------------------------------------------
+ unpack one line from Atoms section of data file
+ initialize other atom quantities
+ AWPMD: 0-tag 1-type 2-q 3-spin 4-eradius 5-etag 6-cs_re 7-cs_im
+------------------------------------------------------------------------- */
+
+void AtomVecWavepacket::data_atom(double *coord, int imagetmp, char **values)
+{
+ int nlocal = atom->nlocal;
+
+ if (nlocal == nmax) grow(0);
+
+ tag[nlocal] = atoi(values[0]);
+ if (tag[nlocal] <= 0)
+ error->one("Invalid atom ID in Atoms section of data file (ID tag must be >0)");
+
+ type[nlocal] = atoi(values[1]);
+ if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
+ error->one("Invalid atom type in Atoms section of data file");
+
+ q[nlocal] = atof(values[2]);
+ spin[nlocal] = atoi(values[3]);
+ eradius[nlocal] = atof(values[4]);
+ if (eradius[nlocal] < 0.0)
+ error->one("Invalid eradius in Atoms section of data file");
+
+
+ etag[nlocal] = atoi(values[5]);
+ cs[2*nlocal] = atoi(values[6]);
+ cs[2*nlocal+1] = atof(values[7]);
+
+
+ x[nlocal][0] = coord[0];
+ x[nlocal][1] = coord[1];
+ x[nlocal][2] = coord[2];
+
+ image[nlocal] = imagetmp;
+
+ mask[nlocal] = 1;
+ v[nlocal][0] = 0.0;
+ v[nlocal][1] = 0.0;
+ v[nlocal][2] = 0.0;
+ ervel[nlocal] = 0.0;
+
+ atom->nlocal++;
+}
+
+/* ----------------------------------------------------------------------
+ unpack hybrid quantities from one line in Atoms section of data file
+ initialize other atom quantities for this sub-style
+------------------------------------------------------------------------- */
+
+int AtomVecWavepacket::data_atom_hybrid(int nlocal, char **values)
+{
+ q[nlocal] = atof(values[0]);
+ spin[nlocal] = atoi(values[1]);
+ eradius[nlocal] = atof(values[2]);
+ if (eradius[nlocal] < 0.0)
+ error->one("Invalid eradius in Atoms section of data file");
+
+ etag[nlocal] = atoi(values[3]);
+ cs[2*nlocal] = atoi(values[4]);
+ cs[2*nlocal+1] = atof(values[5]);
+
+
+ v[nlocal][0] = 0.0;
+ v[nlocal][1] = 0.0;
+ v[nlocal][2] = 0.0;
+ ervel[nlocal] = 0.0;
+
+ return 3;
+}
+
+/* ----------------------------------------------------------------------
+ unpack one line from Velocities section of data file
+------------------------------------------------------------------------- */
+
+void AtomVecWavepacket::data_vel(int m, char **values)
+{
+ v[m][0] = atof(values[0]);
+ v[m][1] = atof(values[1]);
+ v[m][2] = atof(values[2]);
+ ervel[m] = atof(values[3]);
+}
+
+/* ----------------------------------------------------------------------
+ unpack hybrid quantities from one line in Velocities section of data file
+------------------------------------------------------------------------- */
+
+int AtomVecWavepacket::data_vel_hybrid(int m, char **values)
+{
+ ervel[m] = atof(values[0]);
+ return 1;
+}
+
+/* ----------------------------------------------------------------------
+ return # of bytes of allocated memory
+------------------------------------------------------------------------- */
+
+bigint AtomVecWavepacket::memory_usage()
+{
+ bigint bytes = 0;
+
+ if (atom->memcheck("tag")) bytes += memory->usage(tag,nmax);
+ if (atom->memcheck("type")) bytes += memory->usage(type,nmax);
+ if (atom->memcheck("mask")) bytes += memory->usage(mask,nmax);
+ if (atom->memcheck("image")) bytes += memory->usage(image,nmax);
+ if (atom->memcheck("x")) bytes += memory->usage(x,nmax,3);
+ if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3);
+ if (atom->memcheck("f")) bytes += memory->usage(f,nmax,3);
+
+ if (atom->memcheck("q")) bytes += memory->usage(q,nmax);
+ if (atom->memcheck("spin")) bytes += memory->usage(spin,nmax);
+ if (atom->memcheck("eradius")) bytes += memory->usage(eradius,nmax);
+ if (atom->memcheck("ervel")) bytes += memory->usage(ervel,nmax);
+ if (atom->memcheck("erforce")) bytes += memory->usage(erforce,nmax);
+
+ if (atom->memcheck("ervelforce")) bytes += memory->usage(ervelforce,nmax);
+ if (atom->memcheck("cs")) bytes += memory->usage(cs,2*nmax);
+ if (atom->memcheck("csforce")) bytes += memory->usage(csforce,2*nmax);
+ if (atom->memcheck("vforce")) bytes += memory->usage(vforce,3*nmax);
+ if (atom->memcheck("etag")) bytes += memory->usage(etag,nmax);
+
+ return bytes;
+}
diff --git a/src/USER-AWPMD/atom_vec_wavepacket.h b/src/USER-AWPMD/atom_vec_wavepacket.h
new file mode 100644
index 000000000..4515dcf0f
--- /dev/null
+++ b/src/USER-AWPMD/atom_vec_wavepacket.h
@@ -0,0 +1,98 @@
+/* ----------------------------------------------------------------------
+ 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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Ilya Valuev (JIHT RAS)
+------------------------------------------------------------------------- */
+
+
+#ifdef ATOM_CLASS
+
+AtomStyle(wavepacket,AtomVecWavepacket)
+
+#else
+
+#ifndef LMP_ATOM_VEC_WAVEPACKET_H
+#define LMP_ATOM_VEC_WAVEPACKET_H
+
+#include "atom_vec.h"
+
+namespace LAMMPS_NS {
+
+class AtomVecWavepacket : public AtomVec {
+public:
+ AtomVecWavepacket(class LAMMPS *, int, char **);
+ ~AtomVecWavepacket() {}
+ void grow(int);
+ void grow_reset();
+ void copy(int, int, int);
+ int pack_comm(int, int *, double *, int, int *);
+ int pack_comm_vel(int, int *, double *, int, int *);
+ int pack_comm_hybrid(int, int *, double *);
+ void unpack_comm(int, int, double *);
+ void unpack_comm_vel(int, int, double *);
+ int unpack_comm_hybrid(int, int, double *);
+ int pack_reverse(int, int, double *);
+ int pack_reverse_hybrid(int, int, double *);
+ void unpack_reverse(int, int *, double *);
+ int unpack_reverse_hybrid(int, int *, double *);
+ int pack_border(int, int *, double *, int, int *);
+ int pack_border_vel(int, int *, double *, int, int *);
+ int pack_border_hybrid(int, int *, double *);
+ void unpack_border(int, int, double *);
+ void unpack_border_vel(int, int, double *);
+ int unpack_border_hybrid(int, int, double *);
+ int pack_exchange(int, double *);
+ int unpack_exchange(double *);
+ int size_restart();
+ int pack_restart(int, double *);
+ int unpack_restart(double *);
+ void create_atom(int, double *);
+ void data_atom(double *, int, char **);
+ int data_atom_hybrid(int, char **);
+ void data_vel(int, char **);
+ int data_vel_hybrid(int, char **);
+ bigint memory_usage();
+
+private:
+ int *tag,*type,*mask,*image;
+ double **x,**v,**f;
+
+ ///\en spin: -1 or 1 for electron, 0 for ion (compatible with eff)
+ int *spin;
+ ///\en charge: must be specified in the corresponding units (-1 for electron in real units, eff compatible)
+ double *q;
+ ///\en width of the wavepacket (compatible with eff)
+ double *eradius;
+ ///\en width velocity for the wavepacket (compatible with eff)
+ double *ervel;
+ ///\en (generalized) force on width (compatible with eff)
+ double *erforce;
+
+ // AWPMD- specific:
+ ///\en electron tag: must be the same for the WPs belonging to the same electron
+ int *etag;
+ ///\en wavepacket split coeffcients: cre, cim, size is 2*N
+ double *cs;
+ ///\en force on wavepacket split coeffcients: re, im, size is 2*N
+ double *csforce;
+ ///\en (generalized) force on velocity, size is 3*N
+ double *vforce;
+ ///\en (generalized) force on radius velocity, size is N
+ double *ervelforce;
+};
+
+}
+
+#endif
+#endif
diff --git a/src/USER-AWPMD/fix_nve_awpmd.cpp b/src/USER-AWPMD/fix_nve_awpmd.cpp
new file mode 100644
index 000000000..24af55df4
--- /dev/null
+++ b/src/USER-AWPMD/fix_nve_awpmd.cpp
@@ -0,0 +1,155 @@
+/* ----------------------------------------------------------------------
+ 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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Ilya Valuev (JIHT RAS)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdio.h"
+#include "string.h"
+#include "fix_nve_awpmd.h"
+#include "atom.h"
+#include "force.h"
+#include "update.h"
+#include "respa.h"
+#include "error.h"
+#include "math.h"
+
+#include "TCP/wpmd_split.h"
+
+
+
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+FixNVEAwpmd::FixNVEAwpmd(LAMMPS *lmp, int narg, char **arg) :
+ Fix(lmp, narg, arg)
+{
+ if (!atom->wavepacket_flag)
+ error->all("Fix nve/awpmd requires atom style wavepacket");
+ //if (!atom->mass_type != 1)
+ // error->all("Fix nve/awpmd requires per type mass");
+
+ time_integrate = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixNVEAwpmd::setmask()
+{
+ int mask = 0;
+ mask |= INITIAL_INTEGRATE;
+ mask |= FINAL_INTEGRATE;
+ mask |= INITIAL_INTEGRATE_RESPA;
+ mask |= FINAL_INTEGRATE_RESPA;
+ return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEAwpmd::init()
+{
+ dtv = update->dt;
+ dtf = 0.5 * update->dt * force->ftm2v;
+
+ if (strcmp(update->integrate_style,"respa") == 0)
+ step_respa = ((Respa *) update->integrate)->step;
+
+ awpmd_pair=(PairAWPMDCut *)force->pair;
+ awpmd_pair->wpmd->norm_needed=1;
+}
+
+/* ----------------------------------------------------------------------
+ allow for only per-type mass
+------------------------------------------------------------------------- */
+
+void FixNVEAwpmd::initial_integrate(int vflag)
+{
+
+
+ // update v,vr and x,radius of atoms in group
+
+ double **x = atom->x;
+ double *eradius = atom->eradius;
+ double **v = atom->v;
+ double *ervel = atom->ervel;
+ double **f = atom->f;
+ double *erforce = atom->erforce;
+ double *vforce=atom->vforce;
+ double *ervelforce=atom->ervelforce;
+ double *cs=atom->cs;
+ double *csforce=atom->csforce;
+
+
+ double *mass = atom->mass;
+ int *spin = atom->spin;
+ int *type = atom->type;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+ if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+ // x + dt * [v + 0.5 * dt * (f / m)];
+
+ // simple Euler update
+ for (int i = 0; i < nlocal; i++) {
+ if (mask[i] & groupbit) {
+ double dtfm = dtf / mass[type[i]];
+ double dtfmr=dtfm;
+ for(int j=0;j<3;j++){
+ x[i][j] += dtv*vforce[3*i+j];
+ v[i][j] += dtfm*f[i][j];
+ }
+ eradius[i]+= dtv*ervelforce[i];
+ ervel[i] += dtfmr*erforce[i];
+ }
+ }
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEAwpmd::final_integrate(){}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEAwpmd::initial_integrate_respa(int vflag, int ilevel, int iloop)
+{
+ dtv = step_respa[ilevel];
+ dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
+
+ // innermost level - NVE update of v and x
+ // all other levels - NVE update of v
+
+ if (ilevel == 0) initial_integrate(vflag);
+ else final_integrate();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEAwpmd::final_integrate_respa(int ilevel, int iloop)
+{
+ dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
+ final_integrate();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEAwpmd::reset_dt()
+{
+ dtv = update->dt;
+ dtf = 0.5 * update->dt * force->ftm2v;
+}
+
diff --git a/src/USER-AWPMD/fix_nve_awpmd.h b/src/USER-AWPMD/fix_nve_awpmd.h
new file mode 100644
index 000000000..9bc6d3ebb
--- /dev/null
+++ b/src/USER-AWPMD/fix_nve_awpmd.h
@@ -0,0 +1,54 @@
+/* ----------------------------------------------------------------------
+ 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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Ilya Valuev (JIHT RAS)
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(nve/awpmd,FixNVEAwpmd)
+
+#else
+
+#ifndef LMP_FIX_NVE_awpmd_H
+#define LMP_FIX_NVE_awpmd_H
+
+#include "fix.h"
+#include "pair_awpmd_cut.h"
+
+namespace LAMMPS_NS {
+
+class FixNVEAwpmd : public Fix {
+ public:
+ FixNVEAwpmd(class LAMMPS *, int, char **);
+ int setmask();
+ virtual void init();
+ virtual void initial_integrate(int);
+ virtual void final_integrate();
+ void initial_integrate_respa(int, int, int);
+ void final_integrate_respa(int, int);
+ void reset_dt();
+
+ protected:
+ double dtv,dtf;
+ double *step_respa;
+ int mass_require;
+
+ PairAWPMDCut *awpmd_pair;
+};
+
+}
+
+#endif
+#endif
diff --git a/src/USER-AWPMD/pair_awpmd_cut.cpp b/src/USER-AWPMD/pair_awpmd_cut.cpp
new file mode 100644
index 000000000..15e157bd1
--- /dev/null
+++ b/src/USER-AWPMD/pair_awpmd_cut.cpp
@@ -0,0 +1,757 @@
+/* ----------------------------------------------------------------------
+ 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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Ilya Valuev
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdio.h"
+#include "stdlib.h"
+#include "string.h"
+#include "pair_awpmd_cut.h"
+#include "atom.h"
+#include "update.h"
+#include "min.h"
+#include "domain.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "memory.h"
+#include "error.h"
+
+#include "TCP/wpmd_split.h"
+
+
+using namespace LAMMPS_NS;
+
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+
+/* ---------------------------------------------------------------------- */
+
+PairAWPMDCut::PairAWPMDCut(LAMMPS *lmp) : Pair(lmp)
+{
+ single_enable = 0;
+
+ nmax = 0;
+ min_var = NULL;
+ min_varforce = NULL;
+ nextra = 4;
+ pvector = new double[nextra];
+
+ ermscale=1.;
+ width_pbc=0.;
+ wpmd= new AWPMD_split();
+
+ half_box_length=0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairAWPMDCut::~PairAWPMDCut()
+{
+ delete [] pvector;
+ memory->destroy(min_var);
+ memory->destroy(min_varforce);
+
+ if (allocated) {
+ memory->destroy(setflag);
+ memory->destroy(cutsq);
+ memory->destroy(cut);
+ }
+
+ delete wpmd;
+}
+
+
+struct cmp_x{
+ double tol;
+ double **xx;
+ cmp_x(double **xx_=NULL, double tol_=1e-12):xx(xx_),tol(tol_){}
+ bool operator()(const pair<int,int> &left, const pair<int,int> &right) const {
+ if(left.first==right.first){
+ double d=xx[left.second][0]-xx[right.second][0];
+ if(d<-tol)
+ return true;
+ else if(d>tol)
+ return false;
+ d=xx[left.second][1]-xx[right.second][1];
+ if(d<-tol)
+ return true;
+ else if(d>tol)
+ return false;
+ d=xx[left.second][2]-xx[right.second][2];
+ if(d<-tol)
+ return true;
+ else
+ return false;
+ }
+ else
+ return left.first<right.first;
+ }
+};
+
+/* ---------------------------------------------------------------------- */
+
+void PairAWPMDCut::compute(int eflag, int vflag)
+{
+
+ // pvector = [KE, Pauli, ecoul, radial_restraint]
+ for (int i=0; i<4; i++) pvector[i] = 0.0;
+
+ if (eflag || vflag)
+ ev_setup(eflag,vflag);
+ else
+ evflag = vflag_fdotr = 0; //??
+
+ double **x = atom->x;
+ double **f = atom->f;
+ double *q = atom->q;
+ double *erforce = atom->erforce;
+ double *eradius = atom->eradius;
+ int *spin = atom->spin;
+ int *type = atom->type;
+ int *etag = atom->etag;
+ double **v = atom->v;
+
+ int nlocal = atom->nlocal;
+ int nghost = atom->nghost;
+ int ntot=nlocal+nghost;
+
+ int newton_pair = force->newton_pair;
+ double qqrd2e = force->qqrd2e;
+
+ int inum = list->inum;
+ int *ilist = list->ilist;
+ int *numneigh = list->numneigh;
+ int **firstneigh = list->firstneigh;
+
+
+
+
+ // width pbc
+ if(width_pbc<0)
+ wpmd->Lextra=2*half_box_length;
+ else
+ wpmd->Lextra=width_pbc;
+
+ wpmd->newton_pair=newton_pair;
+
+
+
+# if 1
+ // mapping of the LAMMPS numbers to the AWPMC numbers
+ vector<int> gmap(ntot,-1);
+
+ for (int ii = 0; ii < inum; ii++) {
+ int i = ilist[ii];
+ // local particles are all there
+ gmap[i]=0;
+ Vector_3 ri=Vector_3(x[i][0],x[i][1],x[i][2]);
+ int itype = type[i];
+ int *jlist = firstneigh[i];
+ int jnum = numneigh[i];
+ for (int jj = 0; jj < jnum; jj++) {
+ int j = jlist[jj];
+ j &= NEIGHMASK;
+ if(j>=nlocal){ // this is a ghost
+ Vector_3 rj=Vector_3(x[j][0],x[j][1],x[j][2]);
+ int jtype = type[j];
+ double rsq=(ri-rj).norm2();
+ if (rsq < cutsq[itype][jtype])
+ gmap[j]=0; //bingo, this ghost is really needed
+
+ }
+ }
+ }
+
+# else // old mapping
+ // mapping of the LAMMPS numbers to the AWPMC numbers
+ vector<int> gmap(ntot,-1);
+ // map for filtering the clones out: [tag,image] -> id
+ typedef map< pair<int,int>, int, cmp_x > map_t;
+ cmp_x cmp(x);
+ map_t idmap(cmp);
+ for (int ii = 0; ii < inum; ii++) {
+ int i = ilist[ii];
+ // local particles are all there
+ idmap[make_pair(atom->tag[i],i)]=i;
+ bool i_local= i<nlocal ? true : false;
+ if(i_local)
+ gmap[i]=0;
+ else if(gmap[i]==0) // this is a ghost which already has been tested
+ continue;
+ Vector_3 ri=Vector_3(x[i][0],x[i][1],x[i][2]);
+ int itype = type[i];
+ int *jlist = firstneigh[i];
+ int jnum = numneigh[i];
+ for (int jj = 0; jj < jnum; jj++) {
+ int j = jlist[jj];
+ j &= NEIGHMASK;
+
+ pair<map_t::iterator,bool> res=idmap.insert(make_pair(make_pair(atom->tag[j],j),j));
+ bool have_it=!res.second;
+ if(have_it){ // the clone of this particle is already listed
+ if(res.first->second!=j) // check that was not the very same particle
+ gmap[j]=-1; // filter out
+ continue;
+ }
+
+ bool j_local= j<nlocal ? true : false;
+ if((i_local && !j_local) || (j_local && !i_local)){ // some of them is a ghost
+ Vector_3 rj=Vector_3(x[j][0],x[j][1],x[j][2]);
+ int jtype = type[j];
+ double rsq=(ri-rj).norm2();
+ if (rsq < cutsq[itype][jtype]){
+ if(!i_local){
+ gmap[i]=0; //bingo, this ghost is really needed
+ break; // don't need to continue j loop
+ }
+ else
+ gmap[j]=0; //bingo, this ghost is really needed
+ }
+ }
+ }
+ }
+# endif
+ // prepare the solver object
+ wpmd->reset();
+
+ map<int,vector<int> > etmap;
+ // add particles to the AWPMD solver object
+ for (int i = 0; i < ntot; i++) {
+ //int i = ilist[ii];
+ if(gmap[i]<0) // this particle was filtered out
+ continue;
+ if(spin[i]==0) // this is an ion
+ gmap[i]=wpmd->add_ion(q[i], Vector_3(x[i][0],x[i][1],x[i][2]),i<nlocal ? atom->tag[i] : -atom->tag[i]);
+ else if(spin[i]==1 || spin[i]==-1){ // electron, sort them according to the tag
+ etmap[etag[i]].push_back(i);
+ }
+ else
+ error->all(fmt("Invalid spin value (%d) for particle %d !",spin[i],i));
+ }
+ // ion force vector
+ Vector_3 *fi=NULL;
+ if(wpmd->ni)
+ fi= new Vector_3[wpmd->ni];
+
+ // adding electrons
+ for(map<int,vector<int> >::iterator it=etmap.begin(); it!= etmap.end(); ++it){
+ vector<int> &el=it->second;
+ if(!el.size()) // should not happen
+ continue;
+ int s=spin[el[0]] >0 ? 0 : 1;
+ wpmd->add_electron(s); // starts adding the spits
+ for(size_t k=0;k<el.size();k++){
+ int i=el[k];
+ if(spin[el[0]]!=spin[i])
+ error->all(fmt("WP splits for one electron should have the same spin (at particles %d, %d)!",el[0],i));
+ double m= atom->mass ? atom->mass[type[i]] : force->e_mass;
+ Vector_3 xx=Vector_3(x[i][0],x[i][1],x[i][2]);
+ Vector_3 rv=m*Vector_3(v[i][0],v[i][1],v[i][2]);
+ double pv=ermscale*m*atom->ervel[i];
+ Vector_2 cc=Vector_2(atom->cs[2*i],atom->cs[2*i+1]);
+ gmap[i]=wpmd->add_split(xx,rv,atom->eradius[i],pv,cc,1.,atom->q[i],i<nlocal ? atom->tag[i] : -atom->tag[i]);
+ // resetting for the case constraints were applied
+ v[i][0]=rv[0]/m;
+ v[i][1]=rv[1]/m;
+ v[i][2]=rv[2]/m;
+ atom->ervel[i]=pv/(m*ermscale);
+ }
+ }
+ wpmd->set_pbc(NULL); // not required for LAMMPS
+ wpmd->interaction(0x1|0x4|0x10,fi);
+
+ // get forces from the AWPMD solver object
+ for (int ii = 0; ii < inum; ii++) {
+ int i = ilist[ii];
+ if(gmap[i]<0) // this particle was filtered out
+ continue;
+ if(spin[i]==0){ // this is an ion, copying forces
+ int ion=gmap[i];
+ f[i][0]=fi[ion][0];
+ f[i][0]=fi[ion][1];
+ f[i][0]=fi[ion][2];
+ }
+ else { // electron
+ int iel=gmap[i];
+ int s=spin[i] >0 ? 0 : 1;
+ wpmd->get_wp_force(s,iel,(Vector_3 *)f[i],(Vector_3 *)(atom->vforce+3*i),atom->erforce+i,atom->ervelforce+i,(Vector_2 *)(atom->csforce+2*i));
+ }
+ }
+
+ if(fi)
+ delete [] fi;
+
+ // update LAMMPS energy
+ if (eflag_either) {
+ if (eflag_global){
+ eng_coul+= wpmd->get_energy();
+ // pvector = [KE, Pauli, ecoul, radial_restraint]
+ pvector[0] = wpmd->Ee[0]+wpmd->Ee[1];
+ pvector[2] = wpmd->Eii+wpmd->Eei[0]+wpmd->Eei[1]+wpmd->Eee;
+ pvector[1] = pvector[0] + pvector[2] - wpmd->Edk - wpmd->Edc - wpmd->Eii; // All except diagonal terms
+ pvector[3] = wpmd->Ew;
+ }
+
+ if (eflag_atom) {
+ // transfer per-atom energies here
+ for (int i = 0; i < ntot; i++) {
+ if(gmap[i]<0) // this particle was filtered out
+ continue;
+ if(spin[i]==0){
+ eatom[i]=wpmd->Eiep[gmap[i]]+wpmd->Eiip[gmap[i]];
+ }
+ else {
+ int s=spin[i] >0 ? 0 : 1;
+ eatom[i]=wpmd->Eep[s][gmap[i]]+wpmd->Eeip[s][gmap[i]]+wpmd->Eeep[s][gmap[i]]+wpmd->Ewp[s][gmap[i]];
+ }
+ }
+ }
+ }
+ if (vflag_fdotr) {
+ virial_fdotr_compute();
+ if (flexible_pressure_flag)
+ virial_eradius_compute();
+ }
+}
+
+/* ----------------------------------------------------------------------
+ electron width-specific contribution to global virial
+------------------------------------------------------------------------- */
+
+void PairAWPMDCut::virial_eradius_compute()
+{
+ double *eradius = atom->eradius;
+ double *erforce = atom->erforce;
+ double e_virial;
+ int *spin = atom->spin;
+
+ // sum over force on all particles including ghosts
+
+ if (neighbor->includegroup == 0) {
+ int nall = atom->nlocal + atom->nghost;
+ for (int i = 0; i < nall; i++) {
+ if (spin[i]) {
+ e_virial = erforce[i]*eradius[i]/3;
+ virial[0] += e_virial;
+ virial[1] += e_virial;
+ virial[2] += e_virial;
+ }
+ }
+
+ // neighbor includegroup flag is set
+ // sum over force on initial nfirst particles and ghosts
+
+ } else {
+ int nall = atom->nfirst;
+ for (int i = 0; i < nall; i++) {
+ if (spin[i]) {
+ e_virial = erforce[i]*eradius[i]/3;
+ virial[0] += e_virial;
+ virial[1] += e_virial;
+ virial[2] += e_virial;
+ }
+ }
+
+ nall = atom->nlocal + atom->nghost;
+ for (int i = atom->nlocal; i < nall; i++) {
+ if (spin[i]) {
+ e_virial = erforce[i]*eradius[i]/3;
+ virial[0] += e_virial;
+ virial[1] += e_virial;
+ virial[2] += e_virial;
+ }
+ }
+ }
+}
+
+
+
+/* ----------------------------------------------------------------------
+ allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairAWPMDCut::allocate()
+{
+ allocated = 1;
+ int n = atom->ntypes;
+
+ memory->create(setflag,n+1,n+1,"pair:setflag");
+ for (int i = 1; i <= n; i++)
+ for (int j = i; j <= n; j++)
+ setflag[i][j] = 0;
+
+ memory->create(cutsq,n+1,n+1,"pair:cutsq");
+ memory->create(cut,n+1,n+1,"pair:cut");
+}
+
+/* ---------------------------------------------------------------------
+ global settings
+------------------------------------------------------------------------- */
+// the format is: pair_style awpmd/cut [<global_cutoff|-1> [command1] [command2] ...]
+// commands:
+// [hartree|dproduct|uhf] -- quantum approximation level (default is hartree)
+// [free|pbc <length|-1>|fix <w0|-1>|relax|harm <w0>] -- width restriction (default is free)
+// [ermscale <number>] -- scaling factor between electron mass and effective width mass (used for equations of motion only) (default is 1)
+// [flex_press] -- set flexible pressure flag
+// -1 for length means default setting (L/2 for cutoff and L for width PBC)
+void PairAWPMDCut::settings(int narg, char **arg){
+ cut_global=-1.;
+ ermscale=1.;
+ width_pbc=0.;
+ for(int i=0;i<narg;i++){
+ // reading commands
+ // first is cutoff value
+ if(i==0){
+ cut_global = force->numeric(arg[0]);
+ continue;
+ }
+ if(!strcmp(arg[i],"hartree"))
+ wpmd->approx=AWPMD::HARTREE;
+ else if(!strcmp(arg[i],"dproduct"))
+ wpmd->approx=AWPMD::DPRODUCT;
+ else if(!strcmp(arg[i],"uhf"))
+ wpmd->approx=AWPMD::UHF;
+ else if(!strcmp(arg[i],"free"))
+ wpmd->constraint=AWPMD::NONE;
+ else if(!strcmp(arg[i],"fix")){
+ wpmd->constraint=AWPMD::FIX;
+ i++;
+ if(i>=narg)
+ error->all("Setting 'fix' should be followed by a number in awpmd/cut");
+ wpmd->w0=force->numeric(arg[i]);
+ }
+ else if(!strcmp(arg[i],"harm")){
+ wpmd->constraint=AWPMD::HARM;
+ i++;
+ if(i>=narg)
+ error->all("Setting 'harm' should be followed by a number in awpmd/cut");
+ wpmd->w0=force->numeric(arg[i]);
+ wpmd->set_harm_constr(wpmd->w0);
+ }
+ else if(!strcmp(arg[i],"pbc")){
+ i++;
+ if(i>=narg)
+ error->all("Setting 'pbc' should be followed by a number in awpmd/cut");
+ width_pbc=force->numeric(arg[i]);
+ }
+ else if(!strcmp(arg[i],"relax"))
+ wpmd->constraint=AWPMD::RELAX;
+ else if(!strcmp(arg[i],"ermscale")){
+ i++;
+ if(i>=narg)
+ error->all("Setting 'ermscale' should be followed by a number in awpmd/cut");
+ ermscale=force->numeric(arg[i]);
+ }
+ else if(!strcmp(arg[i],"flex_press"))
+ flexible_pressure_flag = 1;
+ }
+
+
+ // reset cutoffs that have been explicitly set
+ /*
+ if (allocated) {
+ int i,j;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i+1; j <= atom->ntypes; j++)
+ if (setflag[i][j]) cut[i][j] = cut_global;
+ }*/
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+// pair settings are as usual
+void PairAWPMDCut::coeff(int narg, char **arg)
+{
+ if (narg < 2 || narg > 3) error->all("Incorrect args for pair coefficients");
+
+ /*if(domain->xperiodic == 1 || domain->yperiodic == 1 ||
+ domain->zperiodic == 1) {*/
+ double delx = domain->boxhi[0]-domain->boxlo[0];
+ double dely = domain->boxhi[1]-domain->boxlo[1];
+ double delz = domain->boxhi[2]-domain->boxlo[2];
+ half_box_length = 0.5 * MIN(delx, MIN(dely, delz));
+ //}
+ if(cut_global<0)
+ cut_global=half_box_length;
+
+ if (!allocated)
+ allocate();
+ else{
+ int i,j;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i+1; j <= atom->ntypes; j++)
+ if (setflag[i][j]) cut[i][j] = cut_global;
+ }
+
+ int ilo,ihi,jlo,jhi;
+ force->bounds(arg[0],atom->ntypes,ilo,ihi);
+ force->bounds(arg[1],atom->ntypes,jlo,jhi);
+
+ double cut_one = cut_global;
+ if (narg == 3) cut_one = atof(arg[2]);
+
+ int count = 0;
+ for (int i = ilo; i <= ihi; i++) {
+ for (int j = MAX(jlo,i); j <= jhi; j++) {
+ cut[i][j] = cut_one;
+ setflag[i][j] = 1;
+ count++;
+ }
+ }
+
+ if (count == 0) error->all("Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+ init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairAWPMDCut::init_style()
+{
+ // error and warning checks
+
+ if (!atom->q_flag || !atom->spin_flag ||
+ !atom->eradius_flag || !atom->erforce_flag ) // TO DO: adjust this to match approximation used
+ error->all("Pair awpmd/cut requires atom attributes "
+ "q, spin, eradius, erforce");
+
+ /*
+ if(vflag_atom){ // can't compute virial per atom
+ //warning->
+ error->all("Pair style awpmd can't compute per atom virials");
+ }*/
+
+ // add hook to minimizer for eradius and erforce
+
+ if (update->whichflag == 2)
+ int ignore = update->minimize->request(this,1,0.01);
+
+ // make sure to use the appropriate timestep when using real units
+
+ /*if (update->whichflag == 1) {
+ if (force->qqr2e == 332.06371 && update->dt == 1.0)
+ error->all("You must lower the default real units timestep for pEFF ");
+ }*/
+
+ // need a half neigh list and optionally a granular history neigh list
+
+ //int irequest = neighbor->request(this);
+
+ //if (atom->tag_enable == 0)
+ // error->all("Pair style reax requires atom IDs");
+
+ //if (force->newton_pair == 0)
+ //error->all("Pair style awpmd requires newton pair on");
+
+ //if (strcmp(update->unit_style,"real") != 0 && comm->me == 0)
+ //error->warning("Not using real units with pair reax");
+
+ int irequest = neighbor->request(this);
+ neighbor->requests[irequest]->newton = 2;
+
+ if(force->e_mass==0. || force->hhmrr2e==0. || force->mvh2r==0.)
+ error->all("Pair style awpmd requires e_mass and conversions hhmrr2e, mvh2r to be properly set for unit system");
+
+ wpmd->me=force->e_mass;
+ wpmd->h2_me=force->hhmrr2e/force->e_mass;
+ wpmd->one_h=force->mvh2r;
+ wpmd->coul_pref=force->qqrd2e;
+
+ wpmd->calc_ii=1;
+}
+
+/* ----------------------------------------------------------------------
+ init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairAWPMDCut::init_one(int i, int j)
+{
+ if (setflag[i][j] == 0)
+ cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
+
+ return cut[i][j];
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairAWPMDCut::write_restart(FILE *fp)
+{
+ write_restart_settings(fp);
+
+ int i,j;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i; j <= atom->ntypes; j++) {
+ fwrite(&setflag[i][j],sizeof(int),1,fp);
+ if (setflag[i][j]) fwrite(&cut[i][j],sizeof(double),1,fp);
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairAWPMDCut::read_restart(FILE *fp)
+{
+ read_restart_settings(fp);
+ allocate();
+
+ int i,j;
+ int me = comm->me;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i; j <= atom->ntypes; j++) {
+ if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
+ MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
+ if (setflag[i][j]) {
+ if (me == 0) fread(&cut[i][j],sizeof(double),1,fp);
+ MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairAWPMDCut::write_restart_settings(FILE *fp)
+{
+ fwrite(&cut_global,sizeof(double),1,fp);
+ fwrite(&offset_flag,sizeof(int),1,fp);
+ fwrite(&mix_flag,sizeof(int),1,fp);
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairAWPMDCut::read_restart_settings(FILE *fp)
+{
+ if (comm->me == 0) {
+ fread(&cut_global,sizeof(double),1,fp);
+ fread(&offset_flag,sizeof(int),1,fp);
+ fread(&mix_flag,sizeof(int),1,fp);
+ }
+ MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
+ MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
+}
+
+/* ----------------------------------------------------------------------
+ returns pointers to the log() of electron radius and corresponding force
+ minimizer operates on log(radius) so radius never goes negative
+ these arrays are stored locally by pair style
+------------------------------------------------------------------------- */
+
+void PairAWPMDCut::min_xf_pointers(int ignore, double **xextra, double **fextra)
+{
+ // grow arrays if necessary
+ // need to be atom->nmax in length
+ int nvar=atom->nmax*(3+1+1+2); // w(1), vel(3), pw(1), cs(2)
+
+ if (nvar > nmax) {
+ memory->destroy(min_var);
+ memory->destroy(min_varforce);
+ nmax = nvar;
+ memory->create(min_var,nmax,"pair:min_var");
+ memory->create(min_varforce,nmax,"pair:min_varforce");
+ }
+
+ *xextra = min_var;
+ *fextra = min_varforce;
+}
+
+/* ----------------------------------------------------------------------
+ minimizer requests the log() of electron radius and corresponding force
+ calculate and store in min_eradius and min_erforce
+------------------------------------------------------------------------- */
+
+void PairAWPMDCut::min_xf_get(int ignore)
+{
+ double *eradius = atom->eradius;
+ double *erforce = atom->erforce;
+ double **v=atom->v;
+ double *vforce=atom->vforce;
+ double *ervel=atom->ervel;
+ double *ervelforce=atom->ervelforce;
+ double *cs=atom->cs;
+ double *csforce=atom->csforce;
+
+ int *spin = atom->spin;
+ int nlocal = atom->nlocal;
+
+ for (int i = 0; i < nlocal; i++)
+ if (spin[i]) {
+ min_var[7*i] = log(eradius[i]);
+ min_varforce[7*i] = eradius[i]*erforce[i];
+ for(int j=0;j<3;j++){
+ min_var[7*i+1+3*j] = v[i][j];
+ min_varforce[7*i+1+3*j] = vforce[3*i+j];
+ }
+ min_var[7*i+4] = ervel[i];
+ min_varforce[7*i+4] = ervelforce[i];
+ min_var[7*i+5] = cs[2*i];
+ min_varforce[7*i+5] = csforce[2*i];
+ min_var[7*i+6] = cs[2*i+1];
+ min_varforce[7*i+6] = csforce[2*i+1];
+
+ } else {
+ for(int j=0;j<7;j++)
+ min_var[7*i+j] = min_varforce[7*i+j] = 0.0;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ propagate the minimizer values to the atom values
+------------------------------------------------------------------------- */
+
+void PairAWPMDCut::min_x_set(int ignore)
+{
+ double *eradius = atom->eradius;
+ double **v=atom->v;
+ double *ervel=atom->ervel;
+ double *cs=atom->cs;
+
+ int *spin = atom->spin;
+ int nlocal = atom->nlocal;
+
+ for (int i = 0; i < nlocal; i++) {
+ if (spin[i]){
+ eradius[i]=exp(min_var[7*i]);
+ for(int j=0;j<3;j++)
+ v[i][j]=min_var[7*i+1+3*j];
+ ervel[i]=min_var[7*i+4];
+ cs[2*i]=min_var[7*i+5];
+ cs[2*i+1]=min_var[7*i+6];
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ memory usage of local atom-based arrays
+------------------------------------------------------------------------- */
+
+double PairAWPMDCut::memory_usage()
+{
+ double bytes = maxeatom * sizeof(double);
+ bytes += maxvatom*6 * sizeof(double);
+ bytes += 2 * nmax * sizeof(double);
+ return bytes;
+}
diff --git a/src/USER-AWPMD/pair_awpmd_cut.h b/src/USER-AWPMD/pair_awpmd_cut.h
new file mode 100644
index 000000000..b1511d449
--- /dev/null
+++ b/src/USER-AWPMD/pair_awpmd_cut.h
@@ -0,0 +1,81 @@
+/* ----------------------------------------------------------------------
+ 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.
+ ------------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+ Contributing author: Ilya Valuev (JIHT RAS)
+------------------------------------------------------------------------- */
+
+
+#ifdef PAIR_CLASS
+
+PairStyle(awpmd/cut,PairAWPMDCut)
+
+#else
+
+#ifndef LMP_PAIR_AWPMD_CUT_H
+#define LMP_PAIR_AWPMD_CUT_H
+
+#include "pair.h"
+
+
+class AWPMD_split;
+
+
+namespace LAMMPS_NS {
+
+class PairAWPMDCut : public Pair {
+ friend class FixNVEAwpmd;
+ public:
+ PairAWPMDCut(class LAMMPS *);
+ virtual ~PairAWPMDCut();
+ virtual void compute(int, int);
+ virtual void settings(int, char **);
+ void coeff(int, char **);
+ void init_style();
+ void min_pointers(double **, double **);
+ double init_one(int, int);
+ void write_restart(FILE *);
+ void read_restart(FILE *);
+ virtual void write_restart_settings(FILE *);
+ virtual void read_restart_settings(FILE *);
+
+ void min_xf_pointers(int, double **, double **);
+ void min_xf_get(int);
+ void min_x_set(int);
+ double memory_usage();
+
+ private:
+
+
+ int flexible_pressure_flag;
+ double cut_global;
+ double **cut;
+
+
+ int nmax; // number of additional variables for minimizer
+ double *min_var,*min_varforce; // additional variables for minimizer
+
+ void allocate();
+
+ void virial_eradius_compute();
+
+
+ AWPMD_split *wpmd; // solver oblect
+ double ermscale; // scale of width mass for motion
+ double width_pbc; // setting for width pbc
+ double half_box_length; // calculated by coeff function
+};
+
+}
+
+#endif
+#endif