diff --git a/src/USER-SMTBQ/pair_smtbq.cpp b/src/USER-SMTBQ/pair_smtbq.cpp
index b2d17109be..f8097a78b4 100644
--- a/src/USER-SMTBQ/pair_smtbq.cpp
+++ b/src/USER-SMTBQ/pair_smtbq.cpp
@@ -38,10 +38,10 @@
.
------------------------------------------------------------------------- */
-#include "math.h"
-#include "stdio.h"
-#include "stdlib.h"
-#include "string.h"
+#include
+#include
+#include
+#include
#include "pair_smtbq.h"
#include "atom.h"
#include "comm.h"
@@ -147,7 +147,6 @@ PairSMTBQ::PairSMTBQ(LAMMPS *lmp) : Pair(lmp)
fct = NULL;
-
maxpage = 0;
// set comm size needed by this Pair
@@ -258,6 +257,9 @@ void PairSMTBQ::coeff(int narg, char **arg)
if (!allocated) allocate();
+ if (strstr(force->pair_style,"hybrid"))
+ error->all(FLERR,"Pair style SMTBQ is not compatible with hybrid styles");
+
if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
@@ -460,11 +462,11 @@ void PairSMTBQ::read_file(char *file)
for (i=1; i <= maxintparam; i++)
intparams[i].mode = (char*) malloc(sizeof(char)*6);
- QEqMode = (char*) malloc(sizeof(char)*18);
- Bavard = (char*) malloc(sizeof(char)*5);
- QInitMode = (char*) malloc(sizeof(char)*18);
- writepot = (char*) malloc(sizeof(char)*5);
- writeenerg = (char*) malloc(sizeof(char)*5);
+ QEqMode = (char*) malloc(sizeof(char)*19);
+ Bavard = (char*) malloc(sizeof(char)*6);
+ QInitMode = (char*) malloc(sizeof(char)*19);
+ writepot = (char*) malloc(sizeof(char)*6);
+ writeenerg = (char*) malloc(sizeof(char)*6);
// Little loop for ion's parameters
@@ -844,7 +846,7 @@ void PairSMTBQ::compute(int eflag, int vflag)
int *ilist,*jlist,*numneigh,**firstneigh;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
- double rsq,iq,jq,Eself,natom;
+ double rsq,iq,jq,Eself;
double ecovtot,ErepOO,ErepMO,Eion,Ecoh;
double **tmp,**tmpAll,*nmol;
double dq,dqcov;
@@ -891,13 +893,7 @@ void PairSMTBQ::compute(int eflag, int vflag)
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
-
-#ifdef LAMMPS_BIGBIG
- long int *tag = atom->tag;
-#else
tagint *tag = atom->tag;
-#endif
-
int *type = atom->type;
int newton_pair = force->newton_pair;
int nlocal = atom->nlocal;
@@ -950,10 +946,6 @@ void PairSMTBQ::compute(int eflag, int vflag)
iq = q[i];
gp = flag_QEq[i];
- if (gp == 0 && itype > 0) natom += 1.0;
- // if (gp == 0 && itype > 0) nmol[gp] += 1.0;
-
-
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
@@ -999,6 +991,7 @@ void PairSMTBQ::compute(int eflag, int vflag)
for (jj = 0; jj < jnum; jj++) {
// ===============================
j = jlist[jj];
+ j &= NEIGHMASK;
jtype = map[type[j]];
jtag = tag[j]; jq = q[j];
@@ -2345,7 +2338,9 @@ void PairSMTBQ::QForce_charge(int loop)
if (loop == 0) {
// ==================
-
+ memset(sbcov,0,sizeof(double)*atom->nmax);
+ memset(coord,0,sizeof(double)*atom->nmax);
+ memset(sbmet,0,sizeof(double)*atom->nmax);
for (ii = 0; ii < inum; ii ++) {
//--------------------------------
@@ -2354,8 +2349,6 @@ void PairSMTBQ::QForce_charge(int loop)
gp = flag_QEq[i];
- sbcov[i] =coord[i]= sbmet[i] = 0.0;
-
itype = map[type[i]];
xtmp = x[i][0];
ytmp = x[i][1];
@@ -2454,6 +2447,7 @@ void PairSMTBQ::QForce_charge(int loop)
for (jj = 0; jj < jnum; jj++)
{
j = jlist[jj];
+ j &= NEIGHMASK;
jtype = map[type[j]];
m = intype[itype][jtype];
jq = q[j];
@@ -2897,16 +2891,17 @@ void PairSMTBQ::groupQEqAllParallel_QEq()
{
int ii,i,jj,j,kk,k,itype,jtype,ktype,jnum,m,gp,zz,z,kgp;
int iproc,team_elt[10][nproc],team_QEq[10][nproc][5];
- int *ilist,*jlist,*numneigh,**firstneigh,ngp,igp,nboite;
+ int *ilist,*jlist,*numneigh,**firstneigh,ngp,igp;
double delr[3],xtmp,ytmp,ztmp,rsq;
int **flag_gp, *nelt, **tab_gp;
int QEq,QEqall[nproc];
double **x = atom->x;
int *type = atom->type;
+ const int nlocal = atom->nlocal;
+ const int nghost = atom->nghost;
+ const int nall = nlocal + nghost;
int inum = list->inum;
- int nlocal = atom->nlocal;
- int nghost = atom->nghost;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
@@ -2916,7 +2911,6 @@ void PairSMTBQ::groupQEqAllParallel_QEq()
// On declare et initialise nos p'tits tableaux
// +++++++++++++++++++++++++++++++++++++++++++++++++
- nboite = nlocal + nghost;
int **tabtemp,**Alltabtemp, *gptmp, *Allgptmp;
memory->create(tabtemp,10*nproc+10,nproc,"pair:tabtemp");
@@ -2924,19 +2918,19 @@ void PairSMTBQ::groupQEqAllParallel_QEq()
memory->create(gptmp,10*nproc+10,"pair:gptmp");
memory->create(Allgptmp,10*nproc+10,"pair:Allgptmp");
- memory->create(flag_gp,nproc,nboite,"pair:flag_gp");
- memory->create(nelt,nboite,"pair:nelt");
- memory->create(tab_gp,10,nboite,"pair:flag_gp");
+ memory->create(flag_gp,nproc,nall,"pair:flag_gp");
+ memory->create(nelt,nall,"pair:nelt");
+ memory->create(tab_gp,10,nall,"pair:flag_gp");
- for (i = 0; i < nlocal+nghost ; i++) { flag_QEq[i] = 0; }
+ for (i = 0; i < nall ; i++) { flag_QEq[i] = 0; }
for (i = 0; i < 10*nproc; i++) {
gptmp[i] = 0; Allgptmp[i] = 0;
for (j=0;jmolecule_flag is set
+ // reset new molecule bond,angle,etc and special values if defined
// send atoms to new owning procs via irregular comm
// since not all atoms I created will be within my sub-domain
// perform special list build if needed
if (mode == MOLECULE) {
+ int molecule_flag = atom->molecule_flag;
+ int molecular = atom->molecular;
+ tagint *molecule = atom->molecule;
+
// molcreate = # of molecules I created
int molcreate = (atom->nlocal - nlocal_previous) / onemol->natoms;
@@ -444,8 +449,7 @@ void CreateAtoms::command(int narg, char **arg)
// including molecules existing before this creation
tagint moloffset;
- tagint *molecule = atom->molecule;
- if (molecule) {
+ if (molecule_flag) {
tagint max = 0;
for (int i = 0; i < nlocal_previous; i++) max = MAX(max,molecule[i]);
tagint maxmol;
@@ -480,13 +484,12 @@ void CreateAtoms::command(int narg, char **arg)
tagint **improper_atom4 = atom->improper_atom4;
int **nspecial = atom->nspecial;
tagint **special = atom->special;
- int molecular = atom->molecular;
int ilocal = nlocal_previous;
for (int i = 0; i < molcreate; i++) {
if (tag) offset = tag[ilocal]-1;
for (int m = 0; m < natoms; m++) {
- if (molecular) molecule[ilocal] = moloffset + i+1;
+ if (molecule_flag) molecule[ilocal] = moloffset + i+1;
if (molecular == 2) {
atom->molindex[ilocal] = 0;
atom->molatom[ilocal] = m;
diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp
index 5f7bec6580..89bd562e8e 100644
--- a/src/fix_nh.cpp
+++ b/src/fix_nh.cpp
@@ -78,6 +78,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
omega_mass_flag = 0;
etap_mass_flag = 0;
flipflag = 1;
+ dipole_flag = 0;
// turn on tilt factor scaling, whenever applicable
@@ -321,6 +322,11 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
else if (strcmp(arg[iarg+1],"no") == 0) flipflag = 0;
else error->all(FLERR,"Illegal fix nvt/npt/nph command");
iarg += 2;
+ } else if (strcmp(arg[iarg],"update") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
+ if (strcmp(arg[iarg+1],"dipole") == 0) dipole_flag = 1;
+ else error->all(FLERR,"Illegal fix nvt/npt/nph command");
+ iarg += 2;
} else if (strcmp(arg[iarg],"fixedpoint") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
fixedpoint[0] = force->numeric(FLERR,arg[iarg+1]);
@@ -417,6 +423,13 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
p_period[0] != p_period[2]))
error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings");
+ if (dipole_flag) {
+ if (!atom->sphere_flag)
+ error->all(FLERR,"Using update dipole flag requires atom style sphere");
+ if (!atom->mu_flag)
+ error->all(FLERR,"Using update dipole flag requires atom attribute mu");
+ }
+
if ((tstat_flag && t_period <= 0.0) ||
(p_flag[0] && p_period[0] <= 0.0) ||
(p_flag[1] && p_period[1] <= 0.0) ||
diff --git a/src/fix_nh.h b/src/fix_nh.h
index be645d0be8..d2201e1612 100644
--- a/src/fix_nh.h
+++ b/src/fix_nh.h
@@ -109,6 +109,7 @@ class FixNH : public Fix {
int eta_mass_flag; // 1 if eta_mass updated, 0 if not.
int omega_mass_flag; // 1 if omega_mass updated, 0 if not.
int etap_mass_flag; // 1 if etap_mass updated, 0 if not.
+ int dipole_flag; // 1 if dipole is updated, 0 if not.
int scaleyz; // 1 if yz scaled with lz
int scalexz; // 1 if xz scaled with lz
diff --git a/src/fix_nh_sphere.cpp b/src/fix_nh_sphere.cpp
index 9db1e34e39..f978b67660 100644
--- a/src/fix_nh_sphere.cpp
+++ b/src/fix_nh_sphere.cpp
@@ -91,6 +91,42 @@ void FixNHSphere::nve_v()
}
}
+/* ----------------------------------------------------------------------
+ perform full-step update of position with dipole orientation, if requested
+-----------------------------------------------------------------------*/
+
+void FixNHSphere::nve_x()
+{
+ // standard nve_x position update
+
+ FixNH::nve_x();
+
+ // update mu for dipoles
+ // d_mu/dt = omega cross mu
+ // renormalize mu to dipole length
+
+ if (dipole_flag) {
+ double msq,scale,g[3];
+ double **mu = atom->mu;
+ double **omega = atom->omega;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit)
+ if (mu[i][3] > 0.0) {
+ g[0] = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]);
+ g[1] = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]);
+ g[2] = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]);
+ msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
+ scale = mu[i][3]/sqrt(msq);
+ mu[i][0] = g[0]*scale;
+ mu[i][1] = g[1]*scale;
+ mu[i][2] = g[2]*scale;
+ }
+ }
+}
+
/* ----------------------------------------------------------------------
perform half-step scaling of rotatonal velocities
-----------------------------------------------------------------------*/
diff --git a/src/fix_nh_sphere.h b/src/fix_nh_sphere.h
index a98356055d..0897f4a93b 100644
--- a/src/fix_nh_sphere.h
+++ b/src/fix_nh_sphere.h
@@ -26,6 +26,7 @@ class FixNHSphere : public FixNH {
protected:
void nve_v();
+ void nve_x();
void nh_v_temp();
};
diff --git a/src/fix_property_atom.cpp b/src/fix_property_atom.cpp
index 10f73aa53b..9b281938c8 100644
--- a/src/fix_property_atom.cpp
+++ b/src/fix_property_atom.cpp
@@ -222,8 +222,21 @@ void FixPropertyAtom::read_data_section(char *keyword, int n, char *buf,
next = strchr(buf,'\n');
values[0] = strtok(buf," \t\n\r\f");
- for (j = 1; j < nwords; j++)
+ if (values[0] == NULL) {
+ char str[128];
+ sprintf(str,"Too few lines in %s section of data file",keyword);
+ error->one(FLERR,str);
+ }
+ int format_ok = 1;
+ for (j = 1; j < nwords; j++) {
values[j] = strtok(NULL," \t\n\r\f");
+ if (values[j] == NULL) format_ok = 0;
+ }
+ if (!format_ok) {
+ char str[128];
+ sprintf(str,"Incorrect %s format in data file",keyword);
+ error->all(FLERR,str);
+ }
itag = ATOTAGINT(values[0]) + id_offset;
if (itag <= 0 || itag > map_tag_max) {