ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

undo cosmetic differences between LAMMPS-ICMS and LAMMPS

Browse Source
This commit is contained in:
Axel Kohlmeyer
2015-09-03 06:51:55 -04:00
parent 31a8e6220b
commit c503b8736d
42 changed files with 268 additions and 288 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/MPIIO/restart_mpiio.cpp
View File

@ -16,7 +16,6 @@
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#include "mpi.h" #include "mpi.h"
#include "mpiio.h"
#include "restart_mpiio.h" #include "restart_mpiio.h"
#include "error.h" #include "error.h"
#include "limits.h" #include "limits.h"

2
src/REPLICA/verlet_split.cpp
View File

@ -194,7 +194,7 @@ VerletSplit::VerletSplit(LAMMPS *lmp, int narg, char **arg) :
// f_kspace = Rspace copy of Kspace forces // f_kspace = Rspace copy of Kspace forces
// allocate dummy version for Kspace partition // allocate dummy version for Kspace partition
maxatom = -1; maxatom = 0;
f_kspace = NULL; f_kspace = NULL;
if (!master) memory->create(f_kspace,1,1,"verlet/split:f_kspace"); if (!master) memory->create(f_kspace,1,1,"verlet/split:f_kspace");
} }

2
src/RIGID/fix_rigid_nh.h
View File

@ -72,7 +72,7 @@ class FixRigidNH : public FixRigid {
int tcomputeflag,pcomputeflag; int tcomputeflag,pcomputeflag;
void couple(); void couple();
virtual void remap(); void remap();
void nhc_temp_integrate(); void nhc_temp_integrate();
void nhc_press_integrate(); void nhc_press_integrate();

5
src/USER-MISC/compute_basal_atom.cpp
View File

@ -204,10 +204,8 @@ void ComputeBasalAtom::compute_peratom()
double bond_angle; double bond_angle;
double norm_j, norm_k; double norm_j, norm_k;
chi[0] = chi[1] = chi[2] = chi[3] = chi[4] = chi[5] = chi[6] = chi[7] = 0; chi[0] = chi[1] = chi[2] = chi[3] = chi[4] = chi[5] = chi[6] = chi[7] = 0;
double x_ij, y_ij, z_ij, x_ik, y_ik, z_ik, double x_ij, y_ij, z_ij, x_ik, y_ik, z_ik,x3[n0],y3[n0],z3[n0],
xmean5, ymean5, zmean5, xmean6, ymean6, zmean6, xmean7, ymean7, zmean7; xmean5, ymean5, zmean5, xmean6, ymean6, zmean6, xmean7, ymean7, zmean7;
double *x3,*y3,*z3;
x3 = new double[n0]; y3 = new double[n0]; z3 = new double[n0];
for (j = 0; j < n0; j++) { for (j = 0; j < n0; j++) {
x_ij = x[i][0]-x[nearest_n0[j]][0]; x_ij = x[i][0]-x[nearest_n0[j]][0];
y_ij = x[i][1]-x[nearest_n0[j]][1]; y_ij = x[i][1]-x[nearest_n0[j]][1];
@ -415,7 +413,6 @@ void ComputeBasalAtom::compute_peratom()
//if there are less than two ~180 degree bond angles, the algorithm returns null //if there are less than two ~180 degree bond angles, the algorithm returns null
else BPV[i][0] = BPV[i][1] = BPV[i][2] = 0.0; else BPV[i][0] = BPV[i][1] = BPV[i][2] = 0.0;
delete[] x3; delete[] y3; delete[] z3;
//normalize BPV: //normalize BPV:
double Mag = sqrt(BPV[i][0]*BPV[i][0] + double Mag = sqrt(BPV[i][0]*BPV[i][0] +
BPV[i][1]*BPV[i][1] + BPV[i][2]*BPV[i][2]); BPV[i][1]*BPV[i][1] + BPV[i][2]*BPV[i][2]);

4
src/USER-MISC/fix_ipi.cpp
View File

@ -380,8 +380,7 @@ void FixIPI::final_integrate()
int nat=bsize/3; int nat=bsize/3;
double **f= atom->f; double **f= atom->f;
double *lbuf; double lbuf[bsize];
lbuf = new double[bsize];
// reassembles the force vector from the local arrays // reassembles the force vector from the local arrays
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
@ -393,7 +392,6 @@ void FixIPI::final_integrate()
lbuf[3*(atom->tag[i]-1)+2]=f[i][2]*forceconv; lbuf[3*(atom->tag[i]-1)+2]=f[i][2]*forceconv;
} }
MPI_Allreduce(lbuf,buffer,bsize,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(lbuf,buffer,bsize,MPI_DOUBLE,MPI_SUM,world);
delete[] lbuf;
for (int i = 0; i < 9; ++i) vir[i]=0.0; for (int i = 0; i < 9; ++i) vir[i]=0.0;

70
src/USER-MISC/fix_srp.cpp
View File

@ -49,22 +49,22 @@ FixSRP::FixSRP(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
// per-atom array width 2 // per-atom array width 2
peratom_flag = 1; peratom_flag = 1;
size_peratom_cols = 2; size_peratom_cols = 2;
// initial allocation of atom-based array // initial allocation of atom-based array
// register with Atom class // register with Atom class
array = NULL; array = NULL;
grow_arrays(atom->nmax); grow_arrays(atom->nmax);
// extends pack_exchange() // extends pack_exchange()
atom->add_callback(0); atom->add_callback(0);
atom->add_callback(1); // restart atom->add_callback(1); // restart
atom->add_callback(2); atom->add_callback(2);
// zero // zero
for (int i = 0; i < atom->nmax; i++) for (int i = 0; i < atom->nmax; i++)
for (int m = 0; m < 3; m++) for (int m = 0; m < 3; m++)
array[i][m] = 0.0; array[i][m] = 0.0;
// assume setup of fix is needed to insert particles // assume setup of fix is needed to insert particles
// is true if reading from datafile // is true if reading from datafile
@ -72,7 +72,7 @@ FixSRP::FixSRP(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
// might be true if reading from restart // might be true if reading from restart
// a restart file written during the run has bond particles as per atom data // a restart file written during the run has bond particles as per atom data
// a restart file written after the run does not have bond particles // a restart file written after the run does not have bond particles
setup = 1; setup = 1;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -106,7 +106,7 @@ void FixSRP::init()
error->all(FLERR,"Cannot use pair srp without pair_style hybrid"); error->all(FLERR,"Cannot use pair srp without pair_style hybrid");
// the highest numbered atom type should be reserved for bond particles (bptype) // the highest numbered atom type should be reserved for bond particles (bptype)
// set bptype, unless it will be read from restart // set bptype, unless it will be read from restart
if(!restart_reset) bptype = atom->ntypes; if(!restart_reset) bptype = atom->ntypes;
// check if bptype is already in use // check if bptype is already in use
@ -118,24 +118,24 @@ void FixSRP::init()
if(!restart_reset) if(!restart_reset)
error->all(FLERR,"Fix SRP requires an extra atom type"); error->all(FLERR,"Fix SRP requires an extra atom type");
else else
setup = 0; setup = 0;
} }
// setup neigh exclusions for diff atom types // setup neigh exclusions for diff atom types
// bond particles do not interact with other types // bond particles do not interact with other types
// type bptype only interacts with itself // type bptype only interacts with itself
char* arg1[4]; char* arg1[4];
arg1[0] = (char *) "exclude"; arg1[0] = (char *) "exclude";
arg1[1] = (char *) "type"; arg1[1] = (char *) "type";
char c0[20]; char c0[20];
char c1[20]; char c1[20];
for(int z = 1; z < bptype; z++) for(int z = 1; z < bptype; z++)
{ {
sprintf(c0, "%d", z); sprintf(c0, "%d", z);
arg1[2] = c0; arg1[2] = c0;
sprintf(c1, "%d", bptype); sprintf(c1, "%d", bptype);
arg1[3] = c1; arg1[3] = c1;
neighbor->modify_params(4, arg1); neighbor->modify_params(4, arg1);
} }
} }
@ -192,19 +192,19 @@ void FixSRP::setup_pre_force(int zz)
xone[1] = (xold[i][1] + xold[j][1])*0.5; xone[1] = (xold[i][1] + xold[j][1])*0.5;
xone[2] = (xold[i][2] + xold[j][2])*0.5; xone[2] = (xold[i][2] + xold[j][2])*0.5;
// record longest bond // record longest bond
// this used to set ghost cutoff // this used to set ghost cutoff
delx = xold[j][0] - xold[i][0]; delx = xold[j][0] - xold[i][0];
dely = xold[j][1] - xold[i][1]; dely = xold[j][1] - xold[i][1];
delz = xold[j][2] - xold[i][2]; delz = xold[j][2] - xold[i][2];
rsq = delx*delx + dely*dely + delz*delz; rsq = delx*delx + dely*dely + delz*delz;
if(rsq > rsqold) rsqold = rsq; if(rsq > rsqold) rsqold = rsq;
// make one particle for each bond // make one particle for each bond
// i is local // i is local
// if newton bond, always make particle // if newton bond, always make particle
// if j is local, always make particle // if j is local, always make particle
// if j is ghost, decide from tag // if j is ghost, decide from tag
if( force->newton_bond || j < nlocal || tagold[i] > tagold[j] ){ if( force->newton_bond || j < nlocal || tagold[i] > tagold[j] ){
atom->natoms += 1; atom->natoms += 1;
@ -246,7 +246,7 @@ void FixSRP::setup_pre_force(int zz)
// ghost atoms must be present for bonds on edge of neighbor cutoff // ghost atoms must be present for bonds on edge of neighbor cutoff
// extend cutghost slightly more than half of the longest bond // extend cutghost slightly more than half of the longest bond
MPI_Allreduce(&rsqold,&rsqmax,1,MPI_DOUBLE,MPI_MAX,world); MPI_Allreduce(&rsqold,&rsqmax,1,MPI_DOUBLE,MPI_MAX,world);
rmax = sqrt(rsqmax); rmax = sqrt(rsqmax);
double cutneighmax_srp = neighbor->cutneighmax + 0.51*rmax; double cutneighmax_srp = neighbor->cutneighmax + 0.51*rmax;
// find smallest cutghost // find smallest cutghost
@ -262,7 +262,7 @@ void FixSRP::setup_pre_force(int zz)
sprintf(str, "Extending ghost comm cutoff. New %f, old %f.", cutneighmax_srp, cutghostmin); sprintf(str, "Extending ghost comm cutoff. New %f, old %f.", cutneighmax_srp, cutghostmin);
error->message(FLERR,str); error->message(FLERR,str);
} }
// cutghost updated by comm->setup // cutghost updated by comm->setup
comm->cutghostuser = cutneighmax_srp; comm->cutghostuser = cutneighmax_srp;
} }
@ -270,7 +270,7 @@ void FixSRP::setup_pre_force(int zz)
// move owned to new procs // move owned to new procs
// get ghosts // get ghosts
// build neigh lists again // build neigh lists again
// if triclinic, lambda coords needed for pbc, exchange, borders // if triclinic, lambda coords needed for pbc, exchange, borders
if (domain->triclinic) domain->x2lamda(atom->nlocal); if (domain->triclinic) domain->x2lamda(atom->nlocal);
domain->pbc(); domain->pbc();
@ -290,7 +290,7 @@ void FixSRP::setup_pre_force(int zz)
// new atom counts // new atom counts
nlocal = atom->nlocal; nlocal = atom->nlocal;
nall = atom->nlocal + atom->nghost; nall = atom->nlocal + atom->nghost;
// zero all forces // zero all forces
for(int i = 0; i < nall; i++) for(int i = 0; i < nall; i++)
for(int n = 0; n < 3; n++) for(int n = 0; n < 3; n++)
atom->f[i][n] = 0.0; atom->f[i][n] = 0.0;
@ -308,7 +308,7 @@ void FixSRP::setup_pre_force(int zz)
void FixSRP::pre_exchange() void FixSRP::pre_exchange()
{ {
// update ghosts // update ghosts
comm->forward_comm(); comm->forward_comm();
// reassign bond particle coordinates to midpoint of bonds // reassign bond particle coordinates to midpoint of bonds
@ -328,7 +328,7 @@ void FixSRP::pre_exchange()
if(j < 0) error->all(FLERR,"Fix SRP failed to map atom"); if(j < 0) error->all(FLERR,"Fix SRP failed to map atom");
j = domain->closest_image(ii,j); j = domain->closest_image(ii,j);
// position of bond particle ii // position of bond particle ii
atom->x[ii][0] = (x[i][0] + x[j][0])*0.5; atom->x[ii][0] = (x[i][0] + x[j][0])*0.5;
atom->x[ii][1] = (x[i][1] + x[j][1])*0.5; atom->x[ii][1] = (x[i][1] + x[j][1])*0.5;
atom->x[ii][2] = (x[i][2] + x[j][2])*0.5; atom->x[ii][2] = (x[i][2] + x[j][2])*0.5;
@ -374,7 +374,7 @@ void FixSRP::copy_arrays(int i, int j, int delflag)
void FixSRP::set_arrays(int i) void FixSRP::set_arrays(int i)
{ {
array[i][0] = -1; array[i][0] = -1;
array[i][1] = -1; array[i][1] = -1;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -438,9 +438,9 @@ int FixSRP::unpack_border(int n, int first, double *buf)
void FixSRP::post_run() void FixSRP::post_run()
{ {
// all bond particles are removed after each run // all bond particles are removed after each run
// useful for write_data and write_restart commands // useful for write_data and write_restart commands
// since those commands occur between runs // since those commands occur between runs
bigint natoms_previous = atom->natoms; bigint natoms_previous = atom->natoms;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int* dlist; int* dlist;
@ -508,7 +508,7 @@ void FixSRP::post_run()
comm->borders(); comm->borders();
// change back to box coordinates // change back to box coordinates
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
pack values in local atom-based arrays for restart file pack values in local atom-based arrays for restart file
@ -518,7 +518,7 @@ int FixSRP::pack_restart(int i, double *buf)
{ {
int m = 0; int m = 0;
buf[m++] = 3; buf[m++] = 3;
buf[m++] = array[i][0]; buf[m++] = array[i][0];
buf[m++] = array[i][1]; buf[m++] = array[i][1];
return m; return m;
} }
@ -527,17 +527,17 @@ int FixSRP::pack_restart(int i, double *buf)
unpack values from atom->extra array to restart the fix unpack values from atom->extra array to restart the fix
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void FixSRP::unpack_restart(int nlocal, int nth) void FixSRP::unpack_restart(int nlocal, int nth)
{ {
double **extra = atom->extra; double **extra = atom->extra;
// skip to Nth set of extra values // skip to Nth set of extra values
int m = 0; int m = 0;
for (int i = 0; i < nth; i++){ for (int i = 0; i < nth; i++){
m += static_cast<int> (extra[nlocal][m]); m += static_cast<int> (extra[nlocal][m]);
} }
m++; m++;
array[nlocal][0] = extra[nlocal][m++]; array[nlocal][0] = extra[nlocal][m++];
array[nlocal][1] = extra[nlocal][m++]; array[nlocal][1] = extra[nlocal][m++];
@ -561,7 +561,7 @@ int FixSRP::size_restart(int nlocal)
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
pack global state of Fix pack global state of Fix
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void FixSRP::write_restart(FILE *fp) void FixSRP::write_restart(FILE *fp)

2
src/USER-MISC/fix_srp.h
View File

@ -43,7 +43,7 @@ class FixSRP : public Fix {
int unpack_exchange(int, double *); int unpack_exchange(int, double *);
int pack_border(int, int *, double *); int pack_border(int, int *, double *);
int unpack_border(int, int, double *); int unpack_border(int, int, double *);
void post_run(); void post_run();
int pack_restart(int, double*); int pack_restart(int, double*);
void unpack_restart(int, int); void unpack_restart(int, int);

184
src/USER-MISC/pair_srp.cpp
View File

@ -10,7 +10,7 @@
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
Contributing authors: Timothy Sirk (ARL), Pieter in't Veld (BASF) Contributing authors: Timothy Sirk (ARL), Pieter in't Veld (BASF)
@ -25,16 +25,16 @@ There is an example script for this package in examples/USER/srp.
Please contact Timothy Sirk for questions (tim.sirk@us.army.mil). Please contact Timothy Sirk for questions (tim.sirk@us.army.mil).
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#include "stdlib.h" #include "stdlib.h"
#include "pair_srp.h" #include "pair_srp.h"
#include "atom.h" #include "atom.h"
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "neighbor.h" #include "neighbor.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
#include "domain.h" #include "domain.h"
#include "modify.h" #include "modify.h"
#include "fix.h" #include "fix.h"
#include "fix_srp.h" #include "fix_srp.h"
@ -42,9 +42,9 @@ Please contact Timothy Sirk for questions (tim.sirk@us.army.mil).
#include "output.h" #include "output.h"
#include "string.h" #include "string.h"
#include "citeme.h" #include "citeme.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
#define SMALL 1.0e-10 #define SMALL 1.0e-10
#define BIG 1e10 #define BIG 1e10
#define ONETWOBIT 0x40000000 #define ONETWOBIT 0x40000000
@ -65,10 +65,10 @@ static int srp_instance = 0;
set size of pair comms in constructor set size of pair comms in constructor
---------------------------------------------------------------------- */ ---------------------------------------------------------------------- */
PairSRP::PairSRP(LAMMPS *lmp) : Pair(lmp) PairSRP::PairSRP(LAMMPS *lmp) : Pair(lmp)
{ {
writedata = 1; writedata = 1;
if (lmp->citeme) lmp->citeme->add(cite_srp); if (lmp->citeme) lmp->citeme->add(cite_srp);
nextra = 1; nextra = 1;
@ -77,8 +77,8 @@ PairSRP::PairSRP(LAMMPS *lmp) : Pair(lmp)
fix_id[0] = '0' + srp_instance / 10; fix_id[0] = '0' + srp_instance / 10;
fix_id[1] = '0' + srp_instance % 10; fix_id[1] = '0' + srp_instance % 10;
++srp_instance; ++srp_instance;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
allocate all arrays allocate all arrays
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -102,19 +102,19 @@ void PairSRP::allocate()
maxcount = 0; maxcount = 0;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
free free
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
PairSRP::~PairSRP() PairSRP::~PairSRP()
{ {
if (allocated) if (allocated)
{ {
memory->destroy(setflag); memory->destroy(setflag);
memory->destroy(cutsq); memory->destroy(cutsq);
memory->destroy(cut); memory->destroy(cut);
memory->destroy(a0); memory->destroy(a0);
memory->destroy(segment); memory->destroy(segment);
} }
// check nfix in case all fixes have already been deleted // check nfix in case all fixes have already been deleted
@ -122,24 +122,24 @@ PairSRP::~PairSRP()
free(fix_id); free(fix_id);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
compute bond-bond repulsions compute bond-bond repulsions
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairSRP::compute(int eflag, int vflag) void PairSRP::compute(int eflag, int vflag)
{ {
// setup energy and virial // setup energy and virial
if (eflag || vflag) if (eflag || vflag)
ev_setup(eflag, vflag); ev_setup(eflag, vflag);
else else
evflag = vflag_fdotr = 0; evflag = vflag_fdotr = 0;
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost; int nall = nlocal + atom->nghost;
int i0, i1, j0, j1; int i0, i1, j0, j1;
int i,j,ii,jj,inum,jnum; int i,j,ii,jj,inum,jnum;
double dijsq, dij; double dijsq, dij;
@ -154,8 +154,8 @@ void PairSRP::compute(int eflag, int vflag)
double fxlever0, fylever0, fzlever0, fxlever1, fylever1, fzlever1; double fxlever0, fylever0, fzlever0, fxlever1, fylever1, fzlever1;
double fx, fy, fz; double fx, fy, fz;
// mapping global to local for atoms inside bond particles // mapping global to local for atoms inside bond particles
// exclude 1-2 neighs if requested // exclude 1-2 neighs if requested
if (neighbor->ago == 0){ if (neighbor->ago == 0){
remapBonds(nall); remapBonds(nall);
if(exclude) onetwoexclude(ilist, inum, jlist, numneigh, firstneigh); if(exclude) onetwoexclude(ilist, inum, jlist, numneigh, firstneigh);
@ -163,9 +163,9 @@ void PairSRP::compute(int eflag, int vflag)
// this pair style only used with hybrid // this pair style only used with hybrid
// due to exclusions // due to exclusions
// each atom i is type bptype // each atom i is type bptype
// each neigh j is type bptype // each neigh j is type bptype
// using midpoint distance option // using midpoint distance option
if(midpoint){ if(midpoint){
@ -193,18 +193,18 @@ void PairSRP::compute(int eflag, int vflag)
// midpt dist bond 0 and 1 // midpt dist bond 0 and 1
dx = 0.5*(x[i0][0] - x[i1][0] + x[j0][0] - x[j1][0]); dx = 0.5*(x[i0][0] - x[i1][0] + x[j0][0] - x[j1][0]);
dy = 0.5*(x[i0][1] - x[i1][1] + x[j0][1] - x[j1][1]); dy = 0.5*(x[i0][1] - x[i1][1] + x[j0][1] - x[j1][1]);
dz = 0.5*(x[i0][2] - x[i1][2] + x[j0][2] - x[j1][2]); dz = 0.5*(x[i0][2] - x[i1][2] + x[j0][2] - x[j1][2]);
dijsq = dx*dx + dy*dy + dz*dz; dijsq = dx*dx + dy*dy + dz*dz;
if (dijsq < cutsq[bptype][bptype]){ if (dijsq < cutsq[bptype][bptype]){
dij = sqrt(dijsq); dij = sqrt(dijsq);
if (dij < SMALL) if (dij < SMALL)
continue; // dij can be 0.0 with soft potentials continue; // dij can be 0.0 with soft potentials
wd = 1.0 - dij / cut[bptype][bptype]; wd = 1.0 - dij / cut[bptype][bptype];
fpair = 0.5 * a0[bptype][bptype] * wd / dij; // 0.5 factor for lever rule fpair = 0.5 * a0[bptype][bptype] * wd / dij; // 0.5 factor for lever rule
// force for bond 0, beads 0,1 // force for bond 0, beads 0,1
//force between bonds //force between bonds
@ -244,10 +244,10 @@ void PairSRP::compute(int eflag, int vflag)
} }
} }
} }
} }
else{ else{
// using min distance option // using min distance option
for (ii = 0; ii < inum; ii++) { for (ii = 0; ii < inum; ii++) {
i = ilist[ii]; i = ilist[ii];
@ -274,13 +274,13 @@ void PairSRP::compute(int eflag, int vflag)
if (dijsq < cutsq[bptype][bptype]){ if (dijsq < cutsq[bptype][bptype]){
dij = sqrt(dijsq); dij = sqrt(dijsq);
if (dij < SMALL) if (dij < SMALL)
continue; // dij can be 0.0 with soft potentials continue; // dij can be 0.0 with soft potentials
wd = 1.0 - dij / cut[bptype][bptype]; wd = 1.0 - dij / cut[bptype][bptype];
fpair = a0[bptype][bptype] * wd / dij; fpair = a0[bptype][bptype] * wd / dij;
// force for bond 0, beads 0,1 // force for bond 0, beads 0,1
lever0 = 0.5 + ti; // assign force according to lever rule lever0 = 0.5 + ti; // assign force according to lever rule
@ -364,7 +364,7 @@ void PairSRP::settings(int narg, char **arg)
int iarg = 3; int iarg = 3;
// default exclude 1-2 // default exclude 1-2
// scaling for 1-2, etc not supported // scaling for 1-2, etc not supported
exclude = 1; exclude = 1;
while (iarg < narg) { while (iarg < narg) {
if (strcmp(arg[iarg],"exclude") == 0) { if (strcmp(arg[iarg],"exclude") == 0) {
@ -392,7 +392,7 @@ void PairSRP::settings(int narg, char **arg)
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
set coeffs set coeffs
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairSRP::coeff(int narg, char **arg) void PairSRP::coeff(int narg, char **arg)
@ -459,7 +459,7 @@ void PairSRP::init_style()
// bond particles do not contribute to energy or virial // bond particles do not contribute to energy or virial
// bond particles do not belong to group all // bond particles do not belong to group all
// but thermo normalization is by nall // but thermo normalization is by nall
// therefore should turn off normalization // therefore should turn off normalization
int me; int me;
MPI_Comm_rank(world,&me); MPI_Comm_rank(world,&me);
@ -467,7 +467,7 @@ void PairSRP::init_style()
arg1[0] = (char *) "norm"; arg1[0] = (char *) "norm";
arg1[1] = (char *) "no"; arg1[1] = (char *) "no";
output->thermo->modify_params(2, arg1); output->thermo->modify_params(2, arg1);
if (me == 0) if (me == 0)
error->message(FLERR,"Thermo normalization turned off by pair srp"); error->message(FLERR,"Thermo normalization turned off by pair srp");
neighbor->request(this,instance_me); neighbor->request(this,instance_me);
@ -487,13 +487,13 @@ double PairSRP::init_one(int i, int j)
return cut[i][j]; return cut[i][j];
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
find min distance for bonds i0/j0 and i1/j1 find min distance for bonds i0/j0 and i1/j1
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
inline void PairSRP::getMinDist(double** &x, double &dx, double &dy, double &dz, double &ti, double &tj, int &i0, int &j0, int &i1, int &j1) inline void PairSRP::getMinDist(double** &x, double &dx, double &dy, double &dz, double &ti, double &tj, int &i0, int &j0, int &i1, int &j1)
{ {
// move these outside the loop // move these outside the loop
double diffx0, diffy0, diffz0, diffx1, diffy1, diffz1, dPx, dPy, dPz, RiRi, RiRj, RjRj; double diffx0, diffy0, diffz0, diffx1, diffy1, diffz1, dPx, dPy, dPz, RiRi, RiRj, RjRj;
double denom, termx0, termy0, termz0, num0, termx1, termy1, termz1, num1; double denom, termx0, termy0, termz0, num0, termx1, termy1, termz1, num1;
@ -503,7 +503,7 @@ inline void PairSRP::getMinDist(double** &x, double &dx, double &dy, double &dz,
diffz0 = x[j0][2] - x[i0][2]; diffz0 = x[j0][2] - x[i0][2];
// compute midpt dist from 1st atom, 2nd bond // compute midpt dist from 1st atom, 2nd bond
diffx1 = x[j1][0] - x[i1][0]; diffx1 = x[j1][0] - x[i1][0];
diffy1 = x[j1][1] - x[i1][1]; diffy1 = x[j1][1] - x[i1][1];
diffz1 = x[j1][2] - x[i1][2]; diffz1 = x[j1][2] - x[i1][2];
@ -519,22 +519,22 @@ inline void PairSRP::getMinDist(double** &x, double &dx, double &dy, double &dz,
denom = RiRj*RiRj - RiRi*RjRj; denom = RiRj*RiRj - RiRi*RjRj;
// handle case of parallel lines // handle case of parallel lines
// reduce to midpt distance // reduce to midpt distance
if (fabs(denom) < SMALL){ if (fabs(denom) < SMALL){
if(denom < 0) denom = -BIG; if(denom < 0) denom = -BIG;
else denom = BIG; else denom = BIG;
} }
// calc ti // calc ti
termx0 = RiRj*diffx1 - RjRj*diffx0; termx0 = RiRj*diffx1 - RjRj*diffx0;
termy0 = RiRj*diffy1 - RjRj*diffy0; termy0 = RiRj*diffy1 - RjRj*diffy0;
termz0 = RiRj*diffz1 - RjRj*diffz0; termz0 = RiRj*diffz1 - RjRj*diffz0;
num0 = dPx*termx0 + dPy*termy0 + dPz*termz0; num0 = dPx*termx0 + dPy*termy0 + dPz*termz0;
ti = num0 / denom; ti = num0 / denom;
if (ti > 0.5) ti = 0.5; if (ti > 0.5) ti = 0.5;
if (ti < -0.5) ti = -0.5; if (ti < -0.5) ti = -0.5;
// calc tj // calc tj
termx1 = RiRj*diffx0 - RiRi*diffx1; termx1 = RiRj*diffx0 - RiRi*diffx1;
termy1 = RiRj*diffy0 - RiRi*diffy1; termy1 = RiRj*diffy0 - RiRi*diffy1;
termz1 = RiRj*diffz0 - RiRi*diffz1; termz1 = RiRj*diffz0 - RiRi*diffz1;
@ -542,28 +542,28 @@ inline void PairSRP::getMinDist(double** &x, double &dx, double &dy, double &dz,
tj = -num1/ denom; tj = -num1/ denom;
if (tj > 0.5) tj = 0.5; if (tj > 0.5) tj = 0.5;
if (tj < -0.5) tj = -0.5; if (tj < -0.5) tj = -0.5;
// min dist // min dist
dx = dPx - ti*diffx0 + tj*diffx1; dx = dPx - ti*diffx0 + tj*diffx1;
dy = dPy - ti*diffy0 + tj*diffy1; dy = dPy - ti*diffy0 + tj*diffy1;
dz = dPz - ti*diffz0 + tj*diffz1; dz = dPz - ti*diffz0 + tj*diffz1;
} }
/* -------------------------------------------------------- /* --------------------------------------------------------
map global id of atoms in stored by each bond particle map global id of atoms in stored by each bond particle
------------------------------------------------------- */ ------------------------------------------------------- */
inline void PairSRP::remapBonds(int &nall) inline void PairSRP::remapBonds(int &nall)
{ {
if(nall > maxcount){ if(nall > maxcount){
memory->grow(segment, nall, 2, "pair:segment"); memory->grow(segment, nall, 2, "pair:segment");
maxcount = nall; maxcount = nall;
} }
// loop over all bond particles // loop over all bond particles
// each bond paricle holds two bond atoms // each bond paricle holds two bond atoms
// map global ids of bond atoms to local ids // map global ids of bond atoms to local ids
// might not be able to map both bond atoms of j, if j is outside neighcut // might not be able to map both bond atoms of j, if j is outside neighcut
// these are not on neighlist, so are not used // these are not on neighlist, so are not used
int tmp; int tmp;
srp = f_srp->array_atom; srp = f_srp->array_atom;
@ -580,18 +580,18 @@ inline void PairSRP::remapBonds(int &nall)
} }
} }
/* -------------------------------------------------------- /* --------------------------------------------------------
add exclusions for 1-2 neighs, if requested add exclusions for 1-2 neighs, if requested
more complex exclusions or scaling probably not needed more complex exclusions or scaling probably not needed
------------------------------------------------------- */ ------------------------------------------------------- */
inline void PairSRP::onetwoexclude(int* &ilist, int &inum, int* &jlist, int* &numneigh, int** &firstneigh) inline void PairSRP::onetwoexclude(int* &ilist, int &inum, int* &jlist, int* &numneigh, int** &firstneigh)
{ {
int i0, i1, j0, j1; int i0, i1, j0, j1;
int i,j,ii,jj,jnum; int i,j,ii,jj,jnum;
// encode neighs with exclusions // encode neighs with exclusions
// only need 1-2 info for normal uses of srp // only need 1-2 info for normal uses of srp
// add 1-3, etc later if ever needed // add 1-3, etc later if ever needed
for (ii = 0; ii < inum; ii++) { for (ii = 0; ii < inum; ii++) {
@ -610,7 +610,7 @@ inline void PairSRP::onetwoexclude(int* &ilist, int &inum, int* &jlist, int* &nu
i1 = segment[j][0]; i1 = segment[j][0];
j1 = segment[j][1]; j1 = segment[j][1];
// check for a 1-2 neigh // check for a 1-2 neigh
if(i0 == i1 || i0 == j1 || i1 == j0 || j0 == j1){ if(i0 == i1 || i0 == j1 || i1 == j0 || j0 == j1){
j |= ONETWOBIT; j |= ONETWOBIT;
jlist[jj] = j; jlist[jj] = j;
@ -639,7 +639,7 @@ void PairSRP::write_data_all(FILE *fp)
for (int j = i; j <= atom->ntypes; j++) for (int j = i; j <= atom->ntypes; j++)
fprintf(fp,"%d %d %g %g\n",i,j,a0[i][j],cut[i][j]); fprintf(fp,"%d %d %g %g\n",i,j,a0[i][j],cut[i][j]);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
proc 0 writes to restart file proc 0 writes to restart file
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */

4
src/USER-MISC/pair_srp.h
View File

@ -44,9 +44,9 @@ class PairSRP : public Pair {
inline void onetwoexclude(int* &, int &, int* &, int* &, int** &); inline void onetwoexclude(int* &, int &, int* &, int* &, int** &);
inline void remapBonds(int &); inline void remapBonds(int &);
void allocate(); void allocate();
void getMinDist(double** &, double &, double &, double &, double &, void getMinDist(double** &, double &, double &, double &, double &,
double &, int &, int &, int &, int &); double &, int &, int &, int &, int &);
bool min, midpoint; bool min, midpoint;
double **cut; double **cut;
double **a0; double **a0;
double **srp; double **srp;

4
src/USER-OMP/fix_nh_asphere_omp.cpp
View File

@ -79,7 +79,7 @@ void FixNHAsphereOMP::nve_v()
const int nlocal = (igroup == atom->firstgroup) ? atom->nfirst : atom->nlocal; const int nlocal = (igroup == atom->firstgroup) ? atom->nfirst : atom->nlocal;
int i; int i;
// standard nve_v velocity update. for efficiency the loop is // standard nve_v velocity update. for efficiency the loop is
// merged with FixNHOMP instead of calling it for the COM update. // merged with FixNHOMP instead of calling it for the COM update.
#if defined(_OPENMP) #if defined(_OPENMP)
@ -122,7 +122,7 @@ void FixNHAsphereOMP::nve_x()
// update quaternion a full step via Richardson iteration // update quaternion a full step via Richardson iteration
// returns new normalized quaternion // returns new normalized quaternion
// principal moments of inertia // principal moments of inertia
#if defined(_OPENMP) #if defined(_OPENMP)
#pragma omp parallel for default(none) private(i) schedule(static) #pragma omp parallel for default(none) private(i) schedule(static)
#endif #endif

2
src/USER-OMP/fix_nh_sphere_omp.cpp
View File

@ -81,7 +81,7 @@ void FixNHSphereOMP::nve_v()
const int nlocal = (igroup == atom->firstgroup) ? atom->nfirst : atom->nlocal; const int nlocal = (igroup == atom->firstgroup) ? atom->nfirst : atom->nlocal;
int i; int i;
// standard nve_v velocity update. for efficiency the loop is // standard nve_v velocity update. for efficiency the loop is
// merged with FixNHOMP instead of calling it for the COM update. // merged with FixNHOMP instead of calling it for the COM update.
// update omega for all particles // update omega for all particles

100
src/USER-OMP/fix_rigid_nh_omp.cpp
View File

@ -70,7 +70,7 @@ void FixRigidNHOMP::initial_integrate(int vflag)
scale_t[0] = scale_t[1] = scale_t[2] = tmp; scale_t[0] = scale_t[1] = scale_t[2] = tmp;
tmp = exp(-dtq * eta_dot_r[0]); tmp = exp(-dtq * eta_dot_r[0]);
scale_r = tmp; scale_r = tmp;
} }
if (pstat_flag) { if (pstat_flag) {
akin_t = akin_r = 0.0; akin_t = akin_r = 0.0;
@ -86,7 +86,7 @@ void FixRigidNHOMP::initial_integrate(int vflag)
tmp = dtq * epsilon_dot[2]; tmp = dtq * epsilon_dot[2];
scale_v[2] = dtv * exp(tmp) * maclaurin_series(tmp); scale_v[2] = dtv * exp(tmp) * maclaurin_series(tmp);
} }
// update xcm, vcm, quat, conjqm and angmom // update xcm, vcm, quat, conjqm and angmom
double akt=0.0, akr=0.0; double akt=0.0, akr=0.0;
int ibody; int ibody;
@ -97,24 +97,24 @@ void FixRigidNHOMP::initial_integrate(int vflag)
for (ibody = 0; ibody < nbody; ibody++) { for (ibody = 0; ibody < nbody; ibody++) {
double mbody[3],tbody[3],fquat[4]; double mbody[3],tbody[3],fquat[4];
const double dtf2 = dtf * 2.0; const double dtf2 = dtf * 2.0;
// step 1.1 - update vcm by 1/2 step // step 1.1 - update vcm by 1/2 step
const double dtfm = dtf / masstotal[ibody]; const double dtfm = dtf / masstotal[ibody];
vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0]; vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0];
vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1]; vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1];
vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2]; vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2];
if (tstat_flag || pstat_flag) { if (tstat_flag || pstat_flag) {
vcm[ibody][0] *= scale_t[0]; vcm[ibody][0] *= scale_t[0];
vcm[ibody][1] *= scale_t[1]; vcm[ibody][1] *= scale_t[1];
vcm[ibody][2] *= scale_t[2]; vcm[ibody][2] *= scale_t[2];
double tmp = vcm[ibody][0]*vcm[ibody][0] + vcm[ibody][1]*vcm[ibody][1] + double tmp = vcm[ibody][0]*vcm[ibody][0] + vcm[ibody][1]*vcm[ibody][1] +
vcm[ibody][2]*vcm[ibody][2]; vcm[ibody][2]*vcm[ibody][2];
akt += masstotal[ibody]*tmp; akt += masstotal[ibody]*tmp;
} }
// step 1.2 - update xcm by full step // step 1.2 - update xcm by full step
if (!pstat_flag) { if (!pstat_flag) {
@ -126,29 +126,29 @@ void FixRigidNHOMP::initial_integrate(int vflag)
xcm[ibody][1] += scale_v[1] * vcm[ibody][1]; xcm[ibody][1] += scale_v[1] * vcm[ibody][1];
xcm[ibody][2] += scale_v[2] * vcm[ibody][2]; xcm[ibody][2] += scale_v[2] * vcm[ibody][2];
} }
// step 1.3 - apply torque (body coords) to quaternion momentum // step 1.3 - apply torque (body coords) to quaternion momentum
torque[ibody][0] *= tflag[ibody][0]; torque[ibody][0] *= tflag[ibody][0];
torque[ibody][1] *= tflag[ibody][1]; torque[ibody][1] *= tflag[ibody][1];
torque[ibody][2] *= tflag[ibody][2]; torque[ibody][2] *= tflag[ibody][2];
MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody],
torque[ibody],tbody); torque[ibody],tbody);
MathExtra::quatvec(quat[ibody],tbody,fquat); MathExtra::quatvec(quat[ibody],tbody,fquat);
conjqm[ibody][0] += dtf2 * fquat[0]; conjqm[ibody][0] += dtf2 * fquat[0];
conjqm[ibody][1] += dtf2 * fquat[1]; conjqm[ibody][1] += dtf2 * fquat[1];
conjqm[ibody][2] += dtf2 * fquat[2]; conjqm[ibody][2] += dtf2 * fquat[2];
conjqm[ibody][3] += dtf2 * fquat[3]; conjqm[ibody][3] += dtf2 * fquat[3];
if (tstat_flag || pstat_flag) { if (tstat_flag || pstat_flag) {
conjqm[ibody][0] *= scale_r; conjqm[ibody][0] *= scale_r;
conjqm[ibody][1] *= scale_r; conjqm[ibody][1] *= scale_r;
conjqm[ibody][2] *= scale_r; conjqm[ibody][2] *= scale_r;
conjqm[ibody][3] *= scale_r; conjqm[ibody][3] *= scale_r;
} }
// step 1.4 to 1.13 - use no_squish rotate to update p and q // step 1.4 to 1.13 - use no_squish rotate to update p and q
MathExtra::no_squish_rotate(3,conjqm[ibody],quat[ibody],inertia[ibody],dtq); MathExtra::no_squish_rotate(3,conjqm[ibody],quat[ibody],inertia[ibody],dtq);
@ -160,22 +160,22 @@ void FixRigidNHOMP::initial_integrate(int vflag)
// update exyz_space // update exyz_space
// transform p back to angmom // transform p back to angmom
// update angular velocity // update angular velocity
MathExtra::q_to_exyz(quat[ibody],ex_space[ibody],ey_space[ibody], MathExtra::q_to_exyz(quat[ibody],ex_space[ibody],ey_space[ibody],
ez_space[ibody]); ez_space[ibody]);
MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody); MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody);
MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody],
mbody,angmom[ibody]); mbody,angmom[ibody]);
angmom[ibody][0] *= 0.5; angmom[ibody][0] *= 0.5;
angmom[ibody][1] *= 0.5; angmom[ibody][1] *= 0.5;
angmom[ibody][2] *= 0.5; angmom[ibody][2] *= 0.5;
MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody], MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
ez_space[ibody],inertia[ibody],omega[ibody]); ez_space[ibody],inertia[ibody],omega[ibody]);
if (tstat_flag || pstat_flag) { if (tstat_flag || pstat_flag) {
akr += angmom[ibody][0]*omega[ibody][0] + akr += angmom[ibody][0]*omega[ibody][0] +
angmom[ibody][1]*omega[ibody][1] + angmom[ibody][2]*omega[ibody][2]; angmom[ibody][1]*omega[ibody][1] + angmom[ibody][2]*omega[ibody][2];
} }
} // end of parallel for } // end of parallel for
@ -205,14 +205,14 @@ void FixRigidNHOMP::initial_integrate(int vflag)
if (vflag) v_setup(vflag); if (vflag) v_setup(vflag);
else evflag = 0; else evflag = 0;
// remap simulation box by 1/2 step // remap simulation box by 1/2 step
if (pstat_flag) remap(); if (pstat_flag) remap();
// set coords/orient and velocity/rotation of atoms in rigid bodies // set coords/orient and velocity/rotation of atoms in rigid bodies
// from quarternion and omega // from quarternion and omega
if (triclinic) if (triclinic)
if (evflag) if (evflag)
set_xv_thr<1,1>(); set_xv_thr<1,1>();
@ -230,7 +230,7 @@ void FixRigidNHOMP::initial_integrate(int vflag)
if (pstat_flag) { if (pstat_flag) {
remap(); remap();
if (kspace_flag) force->kspace->setup(); if (kspace_flag) force->kspace->setup();
} }
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -240,7 +240,7 @@ void FixRigidNHOMP::final_integrate()
double scale_t[3],scale_r; double scale_t[3],scale_r;
// compute scale variables // compute scale variables
scale_t[0] = scale_t[1] = scale_t[2] = 1.0; scale_t[0] = scale_t[1] = scale_t[2] = 1.0;
scale_r = 1.0; scale_r = 1.0;
@ -248,17 +248,17 @@ void FixRigidNHOMP::final_integrate()
double tmp = exp(-1.0 * dtq * eta_dot_t[0]); double tmp = exp(-1.0 * dtq * eta_dot_t[0]);
scale_t[0] = scale_t[1] = scale_t[2] = tmp; scale_t[0] = scale_t[1] = scale_t[2] = tmp;
scale_r = exp(-1.0 * dtq * eta_dot_r[0]); scale_r = exp(-1.0 * dtq * eta_dot_r[0]);
} }
if (pstat_flag) { if (pstat_flag) {
scale_t[0] *= exp(-dtq * (epsilon_dot[0] + mtk_term2)); scale_t[0] *= exp(-dtq * (epsilon_dot[0] + mtk_term2));
scale_t[1] *= exp(-dtq * (epsilon_dot[1] + mtk_term2)); scale_t[1] *= exp(-dtq * (epsilon_dot[1] + mtk_term2));
scale_t[2] *= exp(-dtq * (epsilon_dot[2] + mtk_term2)); scale_t[2] *= exp(-dtq * (epsilon_dot[2] + mtk_term2));
scale_r *= exp(-dtq * (pdim * mtk_term2)); scale_r *= exp(-dtq * (pdim * mtk_term2));
akin_t = akin_r = 0.0; akin_t = akin_r = 0.0;
} }
double * const * _noalias const x = atom->x; double * const * _noalias const x = atom->x;
const dbl3_t * _noalias const f = (dbl3_t *) atom->f[0]; const dbl3_t * _noalias const f = (dbl3_t *) atom->f[0];
const double * const * const torque_one = atom->torque; const double * const * const torque_one = atom->torque;
@ -306,7 +306,7 @@ void FixRigidNHOMP::final_integrate()
// we likely have only a rather number of groups so we loops // we likely have only a rather number of groups so we loops
// over bodies and thread over all atoms for each of them. // over bodies and thread over all atoms for each of them.
for (int ib = 0; ib < nbody; ++ib) { for (int ib = 0; ib < nbody; ++ib) {
double s0=0.0,s1=0.0,s2=0.0,s3=0.0,s4=0.0,s5=0.0; double s0=0.0,s1=0.0,s2=0.0,s3=0.0,s4=0.0,s5=0.0;
int i; int i;
@ -361,7 +361,7 @@ void FixRigidNHOMP::final_integrate()
#else #else
const int tid = 0; const int tid = 0;
#endif #endif
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
const int ibody = body[i]; const int ibody = body[i];
if ((ibody < 0) || (ibody % nthreads != tid)) continue; if ((ibody < 0) || (ibody % nthreads != tid)) continue;
@ -423,28 +423,28 @@ void FixRigidNHOMP::final_integrate()
vcm[ibody][1] *= scale_t[1]; vcm[ibody][1] *= scale_t[1];
vcm[ibody][2] *= scale_t[2]; vcm[ibody][2] *= scale_t[2];
} }
vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0]; vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0];
vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1]; vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1];
vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2]; vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2];
if (pstat_flag) { if (pstat_flag) {
double tmp = vcm[ibody][0]*vcm[ibody][0] + vcm[ibody][1]*vcm[ibody][1] + double tmp = vcm[ibody][0]*vcm[ibody][0] + vcm[ibody][1]*vcm[ibody][1] +
vcm[ibody][2]*vcm[ibody][2]; vcm[ibody][2]*vcm[ibody][2];
akt += masstotal[ibody]*tmp; akt += masstotal[ibody]*tmp;
} }
// update conjqm, then transform to angmom, set velocity again // update conjqm, then transform to angmom, set velocity again
// virial is already setup from initial_integrate // virial is already setup from initial_integrate
torque[ibody][0] *= tflag[ibody][0]; torque[ibody][0] *= tflag[ibody][0];
torque[ibody][1] *= tflag[ibody][1]; torque[ibody][1] *= tflag[ibody][1];
torque[ibody][2] *= tflag[ibody][2]; torque[ibody][2] *= tflag[ibody][2];
MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody], MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],
ez_space[ibody],torque[ibody],tbody); ez_space[ibody],torque[ibody],tbody);
MathExtra::quatvec(quat[ibody],tbody,fquat); MathExtra::quatvec(quat[ibody],tbody,fquat);
if (tstat_flag || pstat_flag) { if (tstat_flag || pstat_flag) {
conjqm[ibody][0] = scale_r * conjqm[ibody][0] + dtf2 * fquat[0]; conjqm[ibody][0] = scale_r * conjqm[ibody][0] + dtf2 * fquat[0];
conjqm[ibody][1] = scale_r * conjqm[ibody][1] + dtf2 * fquat[1]; conjqm[ibody][1] = scale_r * conjqm[ibody][1] + dtf2 * fquat[1];
@ -460,17 +460,17 @@ void FixRigidNHOMP::final_integrate()
MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody); MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody);
MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody],
mbody,angmom[ibody]); mbody,angmom[ibody]);
angmom[ibody][0] *= 0.5; angmom[ibody][0] *= 0.5;
angmom[ibody][1] *= 0.5; angmom[ibody][1] *= 0.5;
angmom[ibody][2] *= 0.5; angmom[ibody][2] *= 0.5;
MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody], MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody],
ez_space[ibody],inertia[ibody],omega[ibody]); ez_space[ibody],inertia[ibody],omega[ibody]);
if (pstat_flag) { if (pstat_flag) {
akr += angmom[ibody][0]*omega[ibody][0] + akr += angmom[ibody][0]*omega[ibody][0] +
angmom[ibody][1]*omega[ibody][1] + angmom[ibody][1]*omega[ibody][1] +
angmom[ibody][2]*omega[ibody][2]; angmom[ibody][2]*omega[ibody][2];
} }
} // end of parallel for } // end of parallel for
@ -523,11 +523,11 @@ void FixRigidNHOMP::remap()
const int nlocal = atom->nlocal; const int nlocal = atom->nlocal;
// epsilon is not used, except for book-keeping // epsilon is not used, except for book-keeping
for (int i = 0; i < 3; i++) epsilon[i] += dtq * epsilon_dot[i]; for (int i = 0; i < 3; i++) epsilon[i] += dtq * epsilon_dot[i];
// convert pertinent atoms and rigid bodies to lamda coords // convert pertinent atoms and rigid bodies to lamda coords
if (allremap) domain->x2lamda(nlocal); if (allremap) domain->x2lamda(nlocal);
else { else {
int i; int i;
@ -538,13 +538,13 @@ void FixRigidNHOMP::remap()
if (mask[i] & dilate_group_bit) if (mask[i] & dilate_group_bit)
domain->x2lamda(x[i],x[i]); domain->x2lamda(x[i],x[i]);
} }
if (nrigid) if (nrigid)
for (int i = 0; i < nrigidfix; i++) for (int i = 0; i < nrigidfix; i++)
modify->fix[rfix[i]]->deform(0); modify->fix[rfix[i]]->deform(0);
// reset global and local box to new size/shape // reset global and local box to new size/shape
for (int i = 0; i < 3; i++) { for (int i = 0; i < 3; i++) {
if (p_flag[i]) { if (p_flag[i]) {
const double oldlo = domain->boxlo[i]; const double oldlo = domain->boxlo[i];
@ -558,9 +558,9 @@ void FixRigidNHOMP::remap()
domain->set_global_box(); domain->set_global_box();
domain->set_local_box(); domain->set_local_box();
// convert pertinent atoms and rigid bodies back to box coords // convert pertinent atoms and rigid bodies back to box coords
if (allremap) domain->lamda2x(nlocal); if (allremap) domain->lamda2x(nlocal);
else { else {
int i; int i;
@ -901,7 +901,7 @@ void FixRigidNHOMP::set_v_thr()
} }
} }
// set omega, angmom of each extended particle // set omega, angmom of each extended particle
// XXX: extended particle info not yet multi-threaded // XXX: extended particle info not yet multi-threaded
if (extended) { if (extended) {

4
src/USER-OMP/fix_rigid_nh_omp.h
View File

@ -5,7 +5,7 @@
Copyright (2003) Sandia Corporation. Under the terms of Contract Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under certain rights in this software. This software is distributed under
the GNU General Public License. the GNU General Public License.
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
@ -26,7 +26,7 @@ class FixRigidNHOMP : public FixRigidNH {
virtual void initial_integrate(int); virtual void initial_integrate(int);
virtual void final_integrate(); virtual void final_integrate();
virtual void remap(); virtual void remap();
private: // copied from FixRigidOMP private: // copied from FixRigidOMP
template <int, int> void set_xv_thr(); template <int, int> void set_xv_thr();

16
src/USER-OMP/fix_rigid_nph_omp.cpp
View File

@ -5,7 +5,7 @@
Copyright (2003) Sandia Corporation. Under the terms of Contract Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under certain rights in this software. This software is distributed under
the GNU General Public License. the GNU General Public License.
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
@ -29,24 +29,24 @@ using namespace LAMMPS_NS;
FixRigidNPHOMP::FixRigidNPHOMP(LAMMPS *lmp, int narg, char **arg) : FixRigidNPHOMP::FixRigidNPHOMP(LAMMPS *lmp, int narg, char **arg) :
FixRigidNHOMP(lmp, narg, arg) FixRigidNHOMP(lmp, narg, arg)
{ {
// other setting are made by parent // other setting are made by parent
scalar_flag = 1; scalar_flag = 1;
restart_global = 1; restart_global = 1;
box_change_size = 1; box_change_size = 1;
extscalar = 1; extscalar = 1;
// error checks // error checks
if (pstat_flag == 0) if (pstat_flag == 0)
error->all(FLERR,"Pressure control must be used with fix nph/omp"); error->all(FLERR,"Pressure control must be used with fix nph/omp");
if (tstat_flag == 1) if (tstat_flag == 1)
error->all(FLERR,"Temperature control must not be used with fix nph/omp"); error->all(FLERR,"Temperature control must not be used with fix nph/omp");
if (p_start[0] < 0.0 || p_start[1] < 0.0 || p_start[2] < 0.0 || if (p_start[0] < 0.0 || p_start[1] < 0.0 || p_start[2] < 0.0 ||
p_stop[0] < 0.0 || p_stop[1] < 0.0 || p_stop[2] < 0.0) p_stop[0] < 0.0 || p_stop[1] < 0.0 || p_stop[2] < 0.0)
error->all(FLERR,"Target pressure for fix rigid/nph/omp cannot be 0.0"); error->all(FLERR,"Target pressure for fix rigid/nph/omp cannot be 0.0");
// convert input periods to frequency // convert input periods to frequency
p_freq[0] = p_freq[1] = p_freq[2] = 0.0; p_freq[0] = p_freq[1] = p_freq[2] = 0.0;
@ -67,15 +67,15 @@ FixRigidNPHOMP::FixRigidNPHOMP(LAMMPS *lmp, int narg, char **arg) :
char **newarg = new char*[3]; char **newarg = new char*[3];
newarg[0] = id_temp; newarg[0] = id_temp;
newarg[1] = (char *) "all"; newarg[1] = (char *) "all";
newarg[2] = (char *) "temp"; newarg[2] = (char *) "temp";
modify->add_compute(3,newarg); modify->add_compute(3,newarg);
delete [] newarg; delete [] newarg;
tcomputeflag = 1; tcomputeflag = 1;
// create a new compute pressure style // create a new compute pressure style
// id = fix-ID + press, compute group = all // id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor // pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7; n = strlen(id) + 7;
id_press = new char[n]; id_press = new char[n];
strcpy(id_press,id); strcpy(id_press,id);

2
src/USER-OMP/fix_rigid_nph_omp.h
View File

@ -5,7 +5,7 @@
Copyright (2003) Sandia Corporation. Under the terms of Contract Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under certain rights in this software. This software is distributed under
the GNU General Public License. the GNU General Public License.
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.

16
src/USER-OMP/fix_rigid_npt_omp.cpp
View File

@ -5,7 +5,7 @@
Copyright (2003) Sandia Corporation. Under the terms of Contract Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under certain rights in this software. This software is distributed under
the GNU General Public License. the GNU General Public License.
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
@ -29,21 +29,21 @@ using namespace LAMMPS_NS;
FixRigidNPTOMP::FixRigidNPTOMP(LAMMPS *lmp, int narg, char **arg) : FixRigidNPTOMP::FixRigidNPTOMP(LAMMPS *lmp, int narg, char **arg) :
FixRigidNHOMP(lmp, narg, arg) FixRigidNHOMP(lmp, narg, arg)
{ {
// other setting are made by parent // other setting are made by parent
scalar_flag = 1; scalar_flag = 1;
restart_global = 1; restart_global = 1;
box_change_size = 1; box_change_size = 1;
extscalar = 1; extscalar = 1;
// error checks // error checks
if (tstat_flag == 0 || pstat_flag == 0) if (tstat_flag == 0 || pstat_flag == 0)
error->all(FLERR,"Did not set temp or press for fix rigid/npt/omp"); error->all(FLERR,"Did not set temp or press for fix rigid/npt/omp");
if (t_start <= 0.0 || t_stop <= 0.0) if (t_start <= 0.0 || t_stop <= 0.0)
error->all(FLERR,"Target temperature for fix rigid/npt/omp cannot be 0.0"); error->all(FLERR,"Target temperature for fix rigid/npt/omp cannot be 0.0");
if (p_start[0] < 0.0 || p_start[1] < 0.0 || p_start[2] < 0.0 || if (p_start[0] < 0.0 || p_start[1] < 0.0 || p_start[2] < 0.0 ||
p_stop[0] < 0.0 || p_stop[1] < 0.0 || p_stop[2] < 0.0) p_stop[0] < 0.0 || p_stop[1] < 0.0 || p_stop[2] < 0.0)
error->all(FLERR,"Target pressure for fix rigid/npt/omp cannot be 0.0"); error->all(FLERR,"Target pressure for fix rigid/npt/omp cannot be 0.0");
@ -53,7 +53,7 @@ FixRigidNPTOMP::FixRigidNPTOMP(LAMMPS *lmp, int narg, char **arg) :
if (t_chain < 1) error->all(FLERR,"Illegal fix_modify command"); if (t_chain < 1) error->all(FLERR,"Illegal fix_modify command");
if (t_iter < 1) error->all(FLERR,"Illegal fix_modify command"); if (t_iter < 1) error->all(FLERR,"Illegal fix_modify command");
if (t_order != 3 && t_order != 5) if (t_order != 3 && t_order != 5)
error->all(FLERR,"Fix_modify order must be 3 or 5"); error->all(FLERR,"Fix_modify order must be 3 or 5");
// convert input periods to frequency // convert input periods to frequency
@ -79,15 +79,15 @@ FixRigidNPTOMP::FixRigidNPTOMP(LAMMPS *lmp, int narg, char **arg) :
char **newarg = new char*[3]; char **newarg = new char*[3];
newarg[0] = id_temp; newarg[0] = id_temp;
newarg[1] = (char *) "all"; newarg[1] = (char *) "all";
newarg[2] = (char *) "temp"; newarg[2] = (char *) "temp";
modify->add_compute(3,newarg); modify->add_compute(3,newarg);
delete [] newarg; delete [] newarg;
tcomputeflag = 1; tcomputeflag = 1;
// create a new compute pressure style // create a new compute pressure style
// id = fix-ID + press, compute group = all // id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor // pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7; n = strlen(id) + 7;
id_press = new char[n]; id_press = new char[n];
strcpy(id_press,id); strcpy(id_press,id);

2
src/USER-OMP/fix_rigid_npt_omp.h
View File

@ -5,7 +5,7 @@
Copyright (2003) Sandia Corporation. Under the terms of Contract Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under certain rights in this software. This software is distributed under
the GNU General Public License. the GNU General Public License.
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.

2
src/USER-OMP/fix_rigid_nve_omp.cpp
View File

@ -5,7 +5,7 @@
Copyright (2003) Sandia Corporation. Under the terms of Contract Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under certain rights in this software. This software is distributed under
the GNU General Public License. the GNU General Public License.
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.

2
src/USER-OMP/fix_rigid_nve_omp.h
View File

@ -5,7 +5,7 @@
Copyright (2003) Sandia Corporation. Under the terms of Contract Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under certain rights in this software. This software is distributed under
the GNU General Public License. the GNU General Public License.
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.

10
src/USER-OMP/fix_rigid_nvt_omp.cpp
View File

@ -5,7 +5,7 @@
Copyright (2003) Sandia Corporation. Under the terms of Contract Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under certain rights in this software. This software is distributed under
the GNU General Public License. the GNU General Public License.
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
@ -26,13 +26,13 @@ using namespace LAMMPS_NS;
FixRigidNVTOMP::FixRigidNVTOMP(LAMMPS *lmp, int narg, char **arg) : FixRigidNVTOMP::FixRigidNVTOMP(LAMMPS *lmp, int narg, char **arg) :
FixRigidNHOMP(lmp, narg, arg) FixRigidNHOMP(lmp, narg, arg)
{ {
// other settings are made by parent // other settings are made by parent
scalar_flag = 1; scalar_flag = 1;
restart_global = 1; restart_global = 1;
extscalar = 1; extscalar = 1;
// error checking // error checking
// convert input period to frequency // convert input period to frequency
@ -46,6 +46,6 @@ FixRigidNVTOMP::FixRigidNVTOMP(LAMMPS *lmp, int narg, char **arg) :
if (t_chain < 1) error->all(FLERR,"Illegal fix_modify command"); if (t_chain < 1) error->all(FLERR,"Illegal fix_modify command");
if (t_iter < 1) error->all(FLERR,"Illegal fix_modify command"); if (t_iter < 1) error->all(FLERR,"Illegal fix_modify command");
if (t_order != 3 && t_order != 5) if (t_order != 3 && t_order != 5)
error->all(FLERR,"Fix_modify order must be 3 or 5"); error->all(FLERR,"Fix_modify order must be 3 or 5");
} }

2
src/USER-OMP/fix_rigid_nvt_omp.h
View File

@ -5,7 +5,7 @@
Copyright (2003) Sandia Corporation. Under the terms of Contract Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under certain rights in this software. This software is distributed under
the GNU General Public License. the GNU General Public License.
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.

8
src/USER-OMP/fix_rigid_omp.cpp
View File

@ -155,9 +155,9 @@ void FixRigidOMP::final_integrate()
} else if (rstyle == GROUP) { } else if (rstyle == GROUP) {
// we likely have only a rather number of groups so we loop // we likely have only a rather number of groups so we loops
// over bodies and thread over all atoms for each of them. // over bodies and thread over all atoms for each of them.
for (int ib = 0; ib < nbody; ++ib) { for (int ib = 0; ib < nbody; ++ib) {
double s0=0.0,s1=0.0,s2=0.0,s3=0.0,s4=0.0,s5=0.0; double s0=0.0,s1=0.0,s2=0.0,s3=0.0,s4=0.0,s5=0.0;
int i; int i;
@ -212,7 +212,7 @@ void FixRigidOMP::final_integrate()
#else #else
const int tid = 0; const int tid = 0;
#endif #endif
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
const int ibody = body[i]; const int ibody = body[i];
if ((ibody < 0) || (ibody % nthreads != tid)) continue; if ((ibody < 0) || (ibody % nthreads != tid)) continue;
@ -616,7 +616,7 @@ void FixRigidOMP::set_v_thr()
} }
} }
// set omega, angmom of each extended particle // set omega, angmom of each extended particle
// XXX: extended particle info not yet multi-threaded // XXX: extended particle info not yet multi-threaded
if (extended) { if (extended) {

10
src/USER-OMP/fix_rigid_small_omp.cpp
View File

@ -120,7 +120,7 @@ void FixRigidSmallOMP::final_integrate()
const int nlocal = atom->nlocal; const int nlocal = atom->nlocal;
const int nthreads=comm->nthreads; const int nthreads=comm->nthreads;
int i, ibody; int i, ibody;
#if defined(_OPENMP) #if defined(_OPENMP)
#pragma omp parallel for default(none) private(ibody) schedule(static) #pragma omp parallel for default(none) private(ibody) schedule(static)
#endif #endif
@ -139,7 +139,7 @@ void FixRigidSmallOMP::final_integrate()
#if defined(_OPENMP) #if defined(_OPENMP)
#pragma omp parallel default(none) private(i,ibody) #pragma omp parallel default(none) private(i,ibody)
#endif #endif
{ {
#if defined(_OPENMP) #if defined(_OPENMP)
const int tid = omp_get_thread_num(); const int tid = omp_get_thread_num();
#else #else
@ -149,7 +149,7 @@ void FixRigidSmallOMP::final_integrate()
for (i = 0; i < nlocal; i++) { for (i = 0; i < nlocal; i++) {
ibody = atom2body[i]; ibody = atom2body[i];
if ((ibody < 0) || (ibody % nthreads != tid)) continue; if ((ibody < 0) || (ibody % nthreads != tid)) continue;
Body &b = body[ibody]; Body &b = body[ibody];
double unwrap[3]; double unwrap[3];
@ -203,7 +203,7 @@ void FixRigidSmallOMP::final_integrate()
} }
// update vcm and angmom, recompute omega // update vcm and angmom, recompute omega
#if defined(_OPENMP) #if defined(_OPENMP)
#pragma omp parallel for default(none) private(ibody) schedule(static) #pragma omp parallel for default(none) private(ibody) schedule(static)
#endif #endif
@ -562,7 +562,7 @@ void FixRigidSmallOMP::set_v_thr()
} }
} }
// set omega, angmom of each extended particle // set omega, angmom of each extended particle
// XXX: extended particle info not yet multi-threaded // XXX: extended particle info not yet multi-threaded
if (extended) { if (extended) {

1
src/USER-OMP/thr_omp.cpp
View File

@ -22,7 +22,6 @@
#include "force.h" #include "force.h"
#include "memory.h" #include "memory.h"
#include "modify.h" #include "modify.h"
#include "compute.h"
#include "neighbor.h" #include "neighbor.h"
#include "timer.h" #include "timer.h"

8
src/USER-PHONON/fix_phonon.cpp
View File

@ -666,10 +666,9 @@ void FixPhonon::postprocess( )
} }
// to get Phi = KT.G^-1; normalization of FFTW data is done here // to get Phi = KT.G^-1; normalization of FFTW data is done here
double boltz = force->boltz, *kbtsqrt, TempAve = 0.; double boltz = force->boltz, kbtsqrt[sysdim], TempAve = 0.;
double TempFac = inv_neval*inv_nTemp; double TempFac = inv_neval*inv_nTemp;
double NormFac = TempFac*double(ntotal); double NormFac = TempFac*double(ntotal);
kbtsqrt = new double[sysdim];
for (idim = 0; idim < sysdim; ++idim){ for (idim = 0; idim < sysdim; ++idim){
kbtsqrt[idim] = sqrt(TempSum[idim]*NormFac); kbtsqrt[idim] = sqrt(TempSum[idim]*NormFac);
@ -683,7 +682,6 @@ void FixPhonon::postprocess( )
for (idim = 0; idim < fft_dim; ++idim) for (idim = 0; idim < fft_dim; ++idim)
for (jdim = 0; jdim < fft_dim; ++jdim) Phi_q[idq][ndim++] *= kbtsqrt[idim%sysdim]*kbtsqrt[jdim%sysdim]; for (jdim = 0; jdim < fft_dim; ++jdim) Phi_q[idq][ndim++] *= kbtsqrt[idim%sysdim]*kbtsqrt[jdim%sysdim];
} }
delete[] kbtsqrt;
// to collect all local Phi_q to root // to collect all local Phi_q to root
displs[0]=0; displs[0]=0;
@ -692,8 +690,7 @@ void FixPhonon::postprocess( )
MPI_Gatherv(Phi_q[0],mynq*fft_dim2*2,MPI_DOUBLE,Phi_all[0],recvcnts,displs,MPI_DOUBLE,0,world); MPI_Gatherv(Phi_q[0],mynq*fft_dim2*2,MPI_DOUBLE,Phi_all[0],recvcnts,displs,MPI_DOUBLE,0,world);
// to collect all basis info and averaged it on root // to collect all basis info and averaged it on root
double *basis_root; double basis_root[fft_dim];
basis_root = new double[fft_dim];
if (fft_dim > sysdim) MPI_Reduce(&basis[1][0], &basis_root[sysdim], fft_dim-sysdim, MPI_DOUBLE, MPI_SUM, 0, world); if (fft_dim > sysdim) MPI_Reduce(&basis[1][0], &basis_root[sysdim], fft_dim-sysdim, MPI_DOUBLE, MPI_SUM, 0, world);
if (me == 0){ // output dynamic matrix by root if (me == 0){ // output dynamic matrix by root
@ -733,7 +730,6 @@ void FixPhonon::postprocess( )
fwrite(M_inv_sqrt, sizeof(double),nucell, fp_bin); fwrite(M_inv_sqrt, sizeof(double),nucell, fp_bin);
fclose(fp_bin); fclose(fp_bin);
delete[] basis_root;
// write log file, here however, it is the dynamical matrix that is written // write log file, here however, it is the dynamical matrix that is written
for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n");

16
src/atom.cpp
View File

@ -735,23 +735,23 @@ void Atom::deallocate_topology()
memory->destroy(atom->angle_atom3); memory->destroy(atom->angle_atom3);
atom->angle_type = NULL; atom->angle_type = NULL;
atom->angle_atom1 = atom->angle_atom2 = atom->angle_atom3 = NULL; atom->angle_atom1 = atom->angle_atom2 = atom->angle_atom3 = NULL;
memory->destroy(atom->dihedral_type); memory->destroy(atom->dihedral_type);
memory->destroy(atom->dihedral_atom1); memory->destroy(atom->dihedral_atom1);
memory->destroy(atom->dihedral_atom2); memory->destroy(atom->dihedral_atom2);
memory->destroy(atom->dihedral_atom3); memory->destroy(atom->dihedral_atom3);
memory->destroy(atom->dihedral_atom4); memory->destroy(atom->dihedral_atom4);
atom->dihedral_type = NULL; atom->dihedral_type = NULL;
atom->dihedral_atom1 = atom->dihedral_atom2 = atom->dihedral_atom1 = atom->dihedral_atom2 =
atom->dihedral_atom3 = atom->dihedral_atom4 = NULL; atom->dihedral_atom3 = atom->dihedral_atom4 = NULL;
memory->destroy(atom->improper_type); memory->destroy(atom->improper_type);
memory->destroy(atom->improper_atom1); memory->destroy(atom->improper_atom1);
memory->destroy(atom->improper_atom2); memory->destroy(atom->improper_atom2);
memory->destroy(atom->improper_atom3); memory->destroy(atom->improper_atom3);
memory->destroy(atom->improper_atom4); memory->destroy(atom->improper_atom4);
atom->improper_type = NULL; atom->improper_type = NULL;
atom->improper_atom1 = atom->improper_atom2 = atom->improper_atom1 = atom->improper_atom2 =
atom->improper_atom3 = atom->improper_atom4 = NULL; atom->improper_atom3 = atom->improper_atom4 = NULL;
} }
@ -865,7 +865,7 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, int type_offset,
(((imageint) (atoi(values[iptr+2]) + IMGMAX) & IMGMASK) << IMG2BITS); (((imageint) (atoi(values[iptr+2]) + IMGMAX) & IMGMASK) << IMG2BITS);
else imagedata = ((imageint) IMGMAX << IMG2BITS) | else imagedata = ((imageint) IMGMAX << IMG2BITS) |
((imageint) IMGMAX << IMGBITS) | IMGMAX; ((imageint) IMGMAX << IMGBITS) | IMGMAX;
xdata[0] = atof(values[xptr]); xdata[0] = atof(values[xptr]);
xdata[1] = atof(values[xptr+1]); xdata[1] = atof(values[xptr+1]);
xdata[2] = atof(values[xptr+2]); xdata[2] = atof(values[xptr+2]);
@ -1078,7 +1078,7 @@ void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset)
for (int i = 0; i < n; i++) { for (int i = 0; i < n; i++) {
next = strchr(buf,'\n'); next = strchr(buf,'\n');
*next = '\0'; *next = '\0';
sscanf(buf,"%d %d " sscanf(buf,"%d %d "
TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT,
&tmp,&itype,&atom1,&atom2,&atom3,&atom4); &tmp,&itype,&atom1,&atom2,&atom3,&atom4);
if (id_offset) { if (id_offset) {
@ -1163,7 +1163,7 @@ void Atom::data_impropers(int n, char *buf, int *count, tagint id_offset)
for (int i = 0; i < n; i++) { for (int i = 0; i < n; i++) {
next = strchr(buf,'\n'); next = strchr(buf,'\n');
*next = '\0'; *next = '\0';
sscanf(buf,"%d %d " sscanf(buf,"%d %d "
TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT,
&tmp,&itype,&atom1,&atom2,&atom3,&atom4); &tmp,&itype,&atom1,&atom2,&atom3,&atom4);
if (id_offset) { if (id_offset) {
@ -1540,7 +1540,7 @@ void Atom::add_molecule_atom(Molecule *onemol, int iatom,
if (onemol->qflag && q_flag) q[ilocal] = onemol->q[iatom]; if (onemol->qflag && q_flag) q[ilocal] = onemol->q[iatom];
if (onemol->radiusflag && radius_flag) radius[ilocal] = onemol->radius[iatom]; if (onemol->radiusflag && radius_flag) radius[ilocal] = onemol->radius[iatom];
if (onemol->rmassflag && rmass_flag) rmass[ilocal] = onemol->rmass[iatom]; if (onemol->rmassflag && rmass_flag) rmass[ilocal] = onemol->rmass[iatom];
else if (rmass_flag) else if (rmass_flag)
rmass[ilocal] = 4.0*MY_PI/3.0 * rmass[ilocal] = 4.0*MY_PI/3.0 *
radius[ilocal]*radius[ilocal]*radius[ilocal]; radius[ilocal]*radius[ilocal]*radius[ilocal];

36
src/atom_vec_body.cpp
View File

@ -65,8 +65,8 @@ AtomVecBody::~AtomVecBody()
{ {
int nall = nlocal_bonus + nghost_bonus; int nall = nlocal_bonus + nghost_bonus;
for (int i = 0; i < nall; i++) { for (int i = 0; i < nall; i++) {
bptr->icp->put(bonus[i].iindex); icp->put(bonus[i].iindex);
bptr->dcp->put(bonus[i].dindex); dcp->put(bonus[i].dindex);
} }
memory->sfree(bonus); memory->sfree(bonus);
@ -95,6 +95,8 @@ void AtomVecBody::process_args(int narg, char **arg)
else error->all(FLERR,"Unknown body style"); else error->all(FLERR,"Unknown body style");
bptr->avec = this; bptr->avec = this;
icp = bptr->icp;
dcp = bptr->dcp;
// max size of forward/border comm // max size of forward/border comm
// 7,16 are packed in pack_comm/pack_border // 7,16 are packed in pack_comm/pack_border
@ -189,8 +191,8 @@ void AtomVecBody::copy(int i, int j, int delflag)
// if deleting atom J via delflag and J has bonus data, then delete it // if deleting atom J via delflag and J has bonus data, then delete it
if (delflag && body[j] >= 0) { if (delflag && body[j] >= 0) {
bptr->icp->put(bonus[body[j]].iindex); icp->put(bonus[body[j]].iindex);
bptr->dcp->put(bonus[body[j]].dindex); dcp->put(bonus[body[j]].dindex);
copy_bonus(nlocal_bonus-1,body[j]); copy_bonus(nlocal_bonus-1,body[j]);
nlocal_bonus--; nlocal_bonus--;
} }
@ -226,8 +228,8 @@ void AtomVecBody::clear_bonus()
{ {
int nall = nlocal_bonus + nghost_bonus; int nall = nlocal_bonus + nghost_bonus;
for (int i = nlocal_bonus; i < nall; i++) { for (int i = nlocal_bonus; i < nall; i++) {
bptr->icp->put(bonus[i].iindex); icp->put(bonus[i].iindex);
bptr->dcp->put(bonus[i].dindex); dcp->put(bonus[i].dindex);
} }
nghost_bonus = 0; nghost_bonus = 0;
} }
@ -828,8 +830,8 @@ void AtomVecBody::unpack_border(int n, int first, double *buf)
inertia[2] = buf[m++]; inertia[2] = buf[m++];
bonus[j].ninteger = (int) ubuf(buf[m++]).i; bonus[j].ninteger = (int) ubuf(buf[m++]).i;
bonus[j].ndouble = (int) ubuf(buf[m++]).i; bonus[j].ndouble = (int) ubuf(buf[m++]).i;
bonus[j].ivalue = bptr->icp->get(bonus[j].ninteger,bonus[j].iindex); bonus[j].ivalue = icp->get(bonus[j].ninteger,bonus[j].iindex);
bonus[j].dvalue = bptr->dcp->get(bonus[j].ndouble,bonus[j].dindex); bonus[j].dvalue = dcp->get(bonus[j].ndouble,bonus[j].dindex);
m += bptr->unpack_border_body(&bonus[j],&buf[m]); m += bptr->unpack_border_body(&bonus[j],&buf[m]);
bonus[j].ilocal = i; bonus[j].ilocal = i;
body[i] = j; body[i] = j;
@ -876,8 +878,8 @@ void AtomVecBody::unpack_border_vel(int n, int first, double *buf)
inertia[2] = buf[m++]; inertia[2] = buf[m++];
bonus[j].ninteger = (int) ubuf(buf[m++]).i; bonus[j].ninteger = (int) ubuf(buf[m++]).i;
bonus[j].ndouble = (int) ubuf(buf[m++]).i; bonus[j].ndouble = (int) ubuf(buf[m++]).i;
bonus[j].ivalue = bptr->icp->get(bonus[j].ninteger,bonus[j].iindex); bonus[j].ivalue = icp->get(bonus[j].ninteger,bonus[j].iindex);
bonus[j].dvalue = bptr->dcp->get(bonus[j].ndouble,bonus[j].dindex); bonus[j].dvalue = dcp->get(bonus[j].ndouble,bonus[j].dindex);
m += bptr->unpack_border_body(&bonus[j],&buf[m]); m += bptr->unpack_border_body(&bonus[j],&buf[m]);
bonus[j].ilocal = i; bonus[j].ilocal = i;
body[i] = j; body[i] = j;
@ -923,8 +925,8 @@ int AtomVecBody::unpack_border_hybrid(int n, int first, double *buf)
inertia[2] = buf[m++]; inertia[2] = buf[m++];
bonus[j].ninteger = (int) ubuf(buf[m++]).i; bonus[j].ninteger = (int) ubuf(buf[m++]).i;
bonus[j].ndouble = (int) ubuf(buf[m++]).i; bonus[j].ndouble = (int) ubuf(buf[m++]).i;
bonus[j].ivalue = bptr->icp->get(bonus[j].ninteger,bonus[j].iindex); bonus[j].ivalue = icp->get(bonus[j].ninteger,bonus[j].iindex);
bonus[j].dvalue = bptr->dcp->get(bonus[j].ndouble,bonus[j].dindex); bonus[j].dvalue = dcp->get(bonus[j].ndouble,bonus[j].dindex);
m += bptr->unpack_border_body(&bonus[j],&buf[m]); m += bptr->unpack_border_body(&bonus[j],&buf[m]);
bonus[j].ilocal = i; bonus[j].ilocal = i;
body[i] = j; body[i] = j;
@ -1027,9 +1029,9 @@ int AtomVecBody::unpack_exchange(double *buf)
inertia[2] = buf[m++]; inertia[2] = buf[m++];
bonus[nlocal_bonus].ninteger = (int) ubuf(buf[m++]).i; bonus[nlocal_bonus].ninteger = (int) ubuf(buf[m++]).i;
bonus[nlocal_bonus].ndouble = (int) ubuf(buf[m++]).i; bonus[nlocal_bonus].ndouble = (int) ubuf(buf[m++]).i;
bonus[nlocal_bonus].ivalue = bptr->icp->get(bonus[nlocal_bonus].ninteger, bonus[nlocal_bonus].ivalue = icp->get(bonus[nlocal_bonus].ninteger,
bonus[nlocal_bonus].iindex); bonus[nlocal_bonus].iindex);
bonus[nlocal_bonus].dvalue = bptr->dcp->get(bonus[nlocal_bonus].ndouble, bonus[nlocal_bonus].dvalue = dcp->get(bonus[nlocal_bonus].ndouble,
bonus[nlocal_bonus].dindex); bonus[nlocal_bonus].dindex);
memcpy(bonus[nlocal_bonus].ivalue,&buf[m], memcpy(bonus[nlocal_bonus].ivalue,&buf[m],
bonus[nlocal_bonus].ninteger*sizeof(int)); bonus[nlocal_bonus].ninteger*sizeof(int));
@ -1179,9 +1181,9 @@ int AtomVecBody::unpack_restart(double *buf)
inertia[2] = buf[m++]; inertia[2] = buf[m++];
bonus[nlocal_bonus].ninteger = (int) ubuf(buf[m++]).i; bonus[nlocal_bonus].ninteger = (int) ubuf(buf[m++]).i;
bonus[nlocal_bonus].ndouble = (int) ubuf(buf[m++]).i; bonus[nlocal_bonus].ndouble = (int) ubuf(buf[m++]).i;
bonus[nlocal_bonus].ivalue = bptr->icp->get(bonus[nlocal_bonus].ninteger, bonus[nlocal_bonus].ivalue = icp->get(bonus[nlocal_bonus].ninteger,
bonus[nlocal_bonus].iindex); bonus[nlocal_bonus].iindex);
bonus[nlocal_bonus].dvalue = bptr->dcp->get(bonus[nlocal_bonus].ndouble, bonus[nlocal_bonus].dvalue = dcp->get(bonus[nlocal_bonus].ndouble,
bonus[nlocal_bonus].dindex); bonus[nlocal_bonus].dindex);
memcpy(bonus[nlocal_bonus].ivalue,&buf[m], memcpy(bonus[nlocal_bonus].ivalue,&buf[m],
bonus[nlocal_bonus].ninteger*sizeof(int)); bonus[nlocal_bonus].ninteger*sizeof(int));
@ -1469,7 +1471,7 @@ bigint AtomVecBody::memory_usage()
if (atom->memcheck("body")) bytes += memory->usage(body,nmax); if (atom->memcheck("body")) bytes += memory->usage(body,nmax);
bytes += nmax_bonus*sizeof(Bonus); bytes += nmax_bonus*sizeof(Bonus);
bytes += bptr->icp->size + bptr->dcp->size; bytes += icp->size + dcp->size;
int nall = nlocal_bonus + nghost_bonus; int nall = nlocal_bonus + nghost_bonus;
for (int i = 0; i < nall; i++) { for (int i = 0; i < nall; i++) {

2
src/delete_bonds.cpp
View File

@ -25,8 +25,6 @@
#include "special.h" #include "special.h"
#include "error.h" #include "error.h"
#include <stdlib.h>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
enum{MULTI,ATOM,BOND,ANGLE,DIHEDRAL,IMPROPER,STATS}; enum{MULTI,ATOM,BOND,ANGLE,DIHEDRAL,IMPROPER,STATS};

1
src/error.h
View File

@ -21,7 +21,6 @@ namespace LAMMPS_NS {
class Error : protected Pointers { class Error : protected Pointers {
public: public:
Error(class LAMMPS *); Error(class LAMMPS *);
virtual ~Error() {}
void universe_all(const char *, int, const char *); void universe_all(const char *, int, const char *);
void universe_one(const char *, int, const char *); void universe_one(const char *, int, const char *);

2
src/fix_ave_spatial.cpp
View File

@ -413,7 +413,7 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
maxvar = 0; maxvar = 0;
varatom = NULL; varatom = NULL;
maxatom = -1; maxatom = 0;
bin = NULL; bin = NULL;
nbins = maxbin = 0; nbins = maxbin = 0;

2
src/fix_heat.cpp
View File

@ -81,7 +81,7 @@ FixHeat::FixHeat(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
scale = 1.0; scale = 1.0;
maxatom = -1; maxatom = 0;
vheat = NULL; vheat = NULL;
vscale = NULL; vscale = NULL;
} }

3
src/fix_nh.h
View File

@ -120,11 +120,10 @@ class FixNH : public Fix {
double fixedpoint[3]; // location of dilation fixed-point double fixedpoint[3]; // location of dilation fixed-point
void couple(); void couple();
void remap();
void nhc_temp_integrate(); void nhc_temp_integrate();
void nhc_press_integrate(); void nhc_press_integrate();
virtual void remap();
virtual void nve_x(); // may be overwritten by child classes virtual void nve_x(); // may be overwritten by child classes
virtual void nve_v(); virtual void nve_v();
virtual void nh_v_press(); virtual void nh_v_press();

2
src/force.h
View File

@ -75,7 +75,7 @@ class Force : protected Pointers {
int special_extra; // extra space for added bonds int special_extra; // extra space for added bonds
Force(class LAMMPS *); Force(class LAMMPS *);
virtual ~Force(); ~Force();
void init(); void init();
void create_pair(const char *, int); void create_pair(const char *, int);

4
src/image.cpp
View File

@ -1065,8 +1065,7 @@ void Image::write_PNG(FILE *fp)
png_set_text(png_ptr,info_ptr,text_ptr,1); png_set_text(png_ptr,info_ptr,text_ptr,1);
png_write_info(png_ptr,info_ptr); png_write_info(png_ptr,info_ptr);
png_bytep *row_pointers; png_bytep row_pointers[height];
row_pointers = new png_bytep[height];
for (int i=0; i < height; ++i) for (int i=0; i < height; ++i)
row_pointers[i] = (png_bytep) &writeBuffer[(height-i-1)*3*width]; row_pointers[i] = (png_bytep) &writeBuffer[(height-i-1)*3*width];
@ -1074,7 +1073,6 @@ void Image::write_PNG(FILE *fp)
png_write_end(png_ptr, info_ptr); png_write_end(png_ptr, info_ptr);
png_destroy_write_struct(&png_ptr, &info_ptr); png_destroy_write_struct(&png_ptr, &info_ptr);
delete[] row_pointers;
#endif #endif
} }

2
src/lattice.cpp
View File

@ -18,8 +18,8 @@
#include "update.h" #include "update.h"
#include "domain.h" #include "domain.h"
#include "comm.h" #include "comm.h"
#include "memory.h"
#include "force.h" #include "force.h"
#include "memory.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;

1
src/memory.h
View File

@ -25,7 +25,6 @@ namespace LAMMPS_NS {
class Memory : protected Pointers { class Memory : protected Pointers {
public: public:
Memory(class LAMMPS *); Memory(class LAMMPS *);
virtual ~Memory() {}
void *smalloc(bigint n, const char *); void *smalloc(bigint n, const char *);
void *srealloc(void *, bigint n, const char *); void *srealloc(void *, bigint n, const char *);

4
src/pair.cpp
View File

@ -223,10 +223,10 @@ void Pair::init()
// I,I coeffs must be set // I,I coeffs must be set
// init_one() will check if I,J is set explicitly or inferred by mixing // init_one() will check if I,J is set explicitly or inferred by mixing
if (!allocated) error->all(FLERR,"Not all pair coeffs are set"); if (!allocated) error->all(FLERR,"All pair coeffs are not set");
for (i = 1; i <= atom->ntypes; i++) for (i = 1; i <= atom->ntypes; i++)
if (setflag[i][i] == 0) error->all(FLERR,"Not all pair coeffs are set"); if (setflag[i][i] == 0) error->all(FLERR,"All pair coeffs are not set");
// style-specific initialization // style-specific initialization

6
src/pair.h
View File

@ -277,10 +277,10 @@ This is likely not what you want to do. The exclusion settings will
eliminate neighbors in the neighbor list, which the manybody potential eliminate neighbors in the neighbor list, which the manybody potential
needs to calculated its terms correctly. needs to calculated its terms correctly.
E: Not all pair coeffs are set E: All pair coeffs are not set
Pair coefficients must be set for all pairs of atom types in either the All pair coefficients must be set in the data file or by the
data file or by the pair_coeff command before running a simulation. pair_coeff command before running a simulation.
E: Fix adapt interface to this pair style not supported E: Fix adapt interface to this pair style not supported

4
src/random_mars.cpp
View File

@ -17,7 +17,6 @@
#include "math.h" #include "math.h"
#include "random_mars.h" #include "random_mars.h"
#include "error.h" #include "error.h"
#include "math_inline.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
@ -105,8 +104,7 @@ double RanMars::gaussian()
rsq = v1*v1 + v2*v2; rsq = v1*v1 + v2*v2;
if (rsq < 1.0 && rsq != 0.0) again = 0; if (rsq < 1.0 && rsq != 0.0) again = 0;
} }
// fac = sqrt(-2.0*log(rsq)/rsq); fac = sqrt(-2.0*log(rsq)/rsq);
fac = MathInline::sqrtlgx2byx(rsq);
second = v1*fac; second = v1*fac;
first = v2*fac; first = v2*fac;
save = 1; save = 1;

4
src/random_park.cpp
View File

@ -15,7 +15,6 @@
#include "math.h" #include "math.h"
#include "random_park.h" #include "random_park.h"
#include "math_inline.h"
#include "error.h" #include "error.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
@ -65,8 +64,7 @@ double RanPark::gaussian()
rsq = v1*v1 + v2*v2; rsq = v1*v1 + v2*v2;
if (rsq < 1.0 && rsq != 0.0) again = 0; if (rsq < 1.0 && rsq != 0.0) again = 0;
} }
// fac = sqrt(-2.0*log(rsq)/rsq); fac = sqrt(-2.0*log(rsq)/rsq);
fac = MathInline::sqrtlgx2byx(rsq);
second = v1*fac; second = v1*fac;
first = v2*fac; first = v2*fac;
save = 1; save = 1;

2
src/universe.h
View File

@ -42,7 +42,7 @@ class Universe : protected Pointers {
// proc uni2orig[I] in original communicator // proc uni2orig[I] in original communicator
Universe(class LAMMPS *, MPI_Comm); Universe(class LAMMPS *, MPI_Comm);
virtual ~Universe(); ~Universe();
void reorder(char *, char *); void reorder(char *, char *);
void add_world(char *); void add_world(char *);
int consistent(); int consistent();

2
src/update.h
View File

@ -47,7 +47,7 @@ class Update : protected Pointers {
char *minimize_style; char *minimize_style;
Update(class LAMMPS *); Update(class LAMMPS *);
virtual ~Update(); ~Update();
void init(); void init();
void set_units(const char *); void set_units(const char *);
void create_integrate(int, char **, int); void create_integrate(int, char **, int);