diff --git a/src/respa.cpp b/src/respa.cpp index 7f09b03f1a..96d947a509 100644 --- a/src/respa.cpp +++ b/src/respa.cpp @@ -22,7 +22,6 @@ #include "atom.h" #include "domain.h" #include "comm.h" -#include "atom.h" #include "force.h" #include "pair.h" #include "bond.h" @@ -148,7 +147,7 @@ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg) if (screen) { fprintf(screen,"Respa levels:\n"); for (int i = 0; i < nlevels; i++) { - fprintf(screen," %d =",i); + fprintf(screen," %d =",i+1); if (level_bond == i) fprintf(screen," bond"); if (level_angle == i) fprintf(screen," angle"); if (level_dihedral == i) fprintf(screen," dihedral"); @@ -164,7 +163,7 @@ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg) if (logfile) { fprintf(logfile,"Respa levels:\n"); for (int i = 0; i < nlevels; i++) { - fprintf(logfile," %d =",i); + fprintf(logfile," %d =",i+1); if (level_bond == i) fprintf(logfile," bond"); if (level_angle == i) fprintf(logfile," angle"); if (level_dihedral == i) fprintf(logfile," dihedral"); @@ -287,6 +286,10 @@ void Respa::init() ev_setup(); + // detect if fix omp is present and will clear force arrays for us + int ifix = modify->find_fix("package_omp"); + if (ifix >= 0) external_force_clear = 1; + // set flags for what arrays to clear in force_clear() // need to clear additionals arrays if they exist @@ -615,6 +618,8 @@ void Respa::recurse(int ilevel) void Respa::force_clear(int newtonflag) { + if (external_force_clear) return; + int i; if (external_force_clear) return; @@ -625,37 +630,15 @@ void Respa::force_clear(int newtonflag) int nall; if (newtonflag) nall = atom->nlocal + atom->nghost; else nall = atom->nlocal; - int ntot = nall * comm->nthreads; - double **f = atom->f; - for (i = 0; i < ntot; i++) { - f[i][0] = 0.0; - f[i][1] = 0.0; - f[i][2] = 0.0; - } + size_t nbytes = sizeof(double) * nall; - if (torqueflag) { - double **torque = atom->torque; - for (i = 0; i < nall; i++) { - torque[i][0] = 0.0; - torque[i][1] = 0.0; - torque[i][2] = 0.0; - } - } - - if (erforceflag) { - double *erforce = atom->erforce; - for (i = 0; i < nall; i++) erforce[i] = 0.0; - } - - if (e_flag) { - double *de = atom->de; - for (i = 0; i < nall; i++) de[i] = 0.0; - } - - if (rho_flag) { - double *drho = atom->drho; - for (i = 0; i < nall; i++) drho[i] = 0.0; + if (nbytes > 0 ) { + memset(&(atom->f[0][0]),0,3*nbytes); + if (torqueflag) memset(&(atom->torque[0][0]),0,3*nbytes); + if (erforceflag) memset(&(atom->erforce[0]), 0, nbytes); + if (e_flag) memset(&(atom->de[0]), 0, nbytes); + if (rho_flag) memset(&(atom->drho[0]), 0, nbytes); } } diff --git a/src/respa.h b/src/respa.h index e806dcc92b..df2d55a0b0 100644 --- a/src/respa.h +++ b/src/respa.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov @@ -40,18 +40,18 @@ class Respa : public Integrate { int level_inner,level_middle,level_outer; Respa(class LAMMPS *, int, char **); - ~Respa(); - void init(); - void setup(); - void setup_minimal(int); - void run(int); - void cleanup(); - void reset_dt(); + virtual ~Respa(); + virtual void init(); + virtual void setup(); + virtual void setup_minimal(int); + virtual void run(int); + virtual void cleanup(); + virtual void reset_dt(); void copy_f_flevel(int); void copy_flevel_f(int); - private: + protected: int triclinic; // 0 if domain is orthog, 1 if triclinic int torqueflag,erforceflag; int e_flag,rho_flag; @@ -59,7 +59,7 @@ class Respa : public Integrate { int *newton; // newton flag at each level class FixRespa *fix_respa; // Fix to store the force level array - void recurse(int); + virtual void recurse(int); void force_clear(int); void sum_flevel_f(); };