diff --git a/src/fix_respa.cpp b/src/fix_respa.cpp index 03c227fd9c..6d3a948a95 100644 --- a/src/fix_respa.cpp +++ b/src/fix_respa.cpp @@ -12,6 +12,7 @@ ------------------------------------------------------------------------- */ #include "stdlib.h" +#include "string.h" #include "fix_respa.h" #include "atom.h" #include "force.h" @@ -29,10 +30,18 @@ FixRespa::FixRespa(LAMMPS *lmp, int narg, char **arg) : nlevels = force->inumeric(FLERR,arg[3]); + // optional arguments + store_torque = 0; + for (int iarg=4; iarg < narg; ++iarg) { + if (strcmp(arg[iarg],"torque") == 0) + store_torque = 1; + } + // perform initial allocation of atom-based arrays // register with Atom class f_level = NULL; + t_level = NULL; grow_arrays(atom->nmax); atom->add_callback(0); } @@ -48,6 +57,7 @@ FixRespa::~FixRespa() // delete locally stored arrays memory->destroy(f_level); + if (store_torque) memory->destroy(t_level); } /* ---------------------------------------------------------------------- */ @@ -64,6 +74,7 @@ int FixRespa::setmask() double FixRespa::memory_usage() { double bytes = atom->nmax*nlevels*3 * sizeof(double); + if (store_torque) bytes += atom->nmax*nlevels*3 * sizeof(double); return bytes; } @@ -74,6 +85,7 @@ double FixRespa::memory_usage() void FixRespa::grow_arrays(int nmax) { memory->grow(f_level,nmax,nlevels,3,"fix_respa:f_level"); + if (store_torque) memory->grow(t_level,nmax,nlevels,3,"fix_respa:t_level"); } /* ---------------------------------------------------------------------- @@ -87,6 +99,13 @@ void FixRespa::copy_arrays(int i, int j, int delflag) f_level[j][k][1] = f_level[i][k][1]; f_level[j][k][2] = f_level[i][k][2]; } + if (store_torque) { + for (int k = 0; k < nlevels; k++) { + t_level[j][k][0] = t_level[i][k][0]; + t_level[j][k][1] = t_level[i][k][1]; + t_level[j][k][2] = t_level[i][k][2]; + } + } } /* ---------------------------------------------------------------------- @@ -101,6 +120,13 @@ int FixRespa::pack_exchange(int i, double *buf) buf[m++] = f_level[i][k][1]; buf[m++] = f_level[i][k][2]; } + if (store_torque) { + for (int k = 0; k < nlevels; k++) { + buf[m++] = t_level[i][k][0]; + buf[m++] = t_level[i][k][1]; + buf[m++] = t_level[i][k][2]; + } + } return m; } @@ -116,5 +142,12 @@ int FixRespa::unpack_exchange(int nlocal, double *buf) f_level[nlocal][k][1] = buf[m++]; f_level[nlocal][k][2] = buf[m++]; } + if (store_torque) { + for (int k = 0; k < nlevels; k++) { + t_level[nlocal][k][0] = buf[m++]; + t_level[nlocal][k][1] = buf[m++]; + t_level[nlocal][k][2] = buf[m++]; + } + } return m; } diff --git a/src/fix_respa.h b/src/fix_respa.h index 7c7bd53a61..96b5b5531b 100644 --- a/src/fix_respa.h +++ b/src/fix_respa.h @@ -43,7 +43,9 @@ class FixRespa : public Fix { private: int nlevels; - double ***f_level; // force at each rRESPA level + int store_torque; // 1 if torques should be stored in addition to forces + double ***f_level; // force at each rRESPA level + double ***t_level; // torque at each rRESPA level }; } diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index e24f9788ee..dce3d21978 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -25,6 +25,7 @@ #include "comm.h" #include "memory.h" #include "error.h" +#include "respa.h" using namespace LAMMPS_NS; @@ -36,8 +37,11 @@ PairHybrid::PairHybrid(LAMMPS *lmp) : Pair(lmp) styles = NULL; keywords = NULL; multiple = NULL; + special_lj = NULL; + special_coul = NULL; outerflag = 0; + respaflag = 0; } /* ---------------------------------------------------------------------- */ @@ -45,13 +49,20 @@ PairHybrid::PairHybrid(LAMMPS *lmp) : Pair(lmp) PairHybrid::~PairHybrid() { if (nstyles) { - for (int m = 0; m < nstyles; m++) delete styles[m]; - for (int m = 0; m < nstyles; m++) delete [] keywords[m]; + for (int m = 0; m < nstyles; m++) { + delete styles[m]; + delete [] keywords[m]; + if (special_lj[m]) delete [] special_lj[m]; + if (special_coul[m]) delete [] special_coul[m]; + } } delete [] styles; delete [] keywords; delete [] multiple; + delete [] special_lj; + delete [] special_coul; + delete [] svector; if (allocated) { @@ -98,15 +109,36 @@ void PairHybrid::compute(int eflag, int vflag) if (vflag % 4 == 2) vflag_substyle = vflag/4 * 4; else vflag_substyle = vflag; + double *saved_special = save_special(); + + // check if we are running with r-RESPA using the hybrid keyword + + Respa *respa = NULL; + if (strstr(update->integrate_style,"respa")) { + respa = (Respa *) update->integrate; + if (respa->nhybrid_styles > 0) respaflag = 1; + } + for (m = 0; m < nstyles; m++) { - // invoke compute() unless compute flag is turned off or - // outerflag is set and sub-style has a compute_outer() method + set_special(m); - if (styles[m]->compute_flag == 0) continue; - if (outerflag && styles[m]->respa_enable) - styles[m]->compute_outer(eflag,vflag_substyle); - else styles[m]->compute(eflag,vflag_substyle); + if (!respaflag || (respaflag && respa->hybrid_compute[m])) { + + // invoke compute() unless compute flag is turned off or + // outerflag is set and sub-style has a compute_outer() method + + if (styles[m]->compute_flag == 0) continue; + if (outerflag && styles[m]->respa_enable) + styles[m]->compute_outer(eflag,vflag_substyle); + else styles[m]->compute(eflag,vflag_substyle); + } + + restore_special(saved_special); + + // jump to next sub-style if r-RESPA does not want global accumulated data + + if (respaflag && !respa->tally_global) continue; if (eflag_global) { eng_vdwl += styles[m]->eng_vdwl; @@ -131,6 +163,8 @@ void PairHybrid::compute(int eflag, int vflag) } } + delete [] saved_special; + if (vflag_fdotr) virial_fdotr_compute(); } @@ -215,6 +249,9 @@ void PairHybrid::settings(int narg, char **arg) keywords = new char*[narg]; multiple = new int[narg]; + special_lj = new double*[narg]; + special_coul = new double*[narg]; + // allocate each sub-style // allocate uses suffix, but don't store suffix version in keywords, // else syntax in coeff() will not match @@ -233,6 +270,7 @@ void PairHybrid::settings(int narg, char **arg) styles[nstyles] = force->new_pair(arg[iarg],1,dummy); force->store_style(keywords[nstyles],arg[iarg],0); + special_lj[nstyles] = special_coul[nstyles] = NULL; jarg = iarg + 1; while (jarg < narg && !force->pair_map->count(arg[jarg])) jarg++; @@ -421,6 +459,28 @@ void PairHybrid::init_style() if (used == 0) error->all(FLERR,"Pair hybrid sub-style is not used"); } + // check if special_lj/special_coul overrides are compatible + for (istyle = 0; istyle < nstyles; istyle++) { + if (special_lj[istyle]) { + for (i = 1; i < 4; ++i) { + if (((force->special_lj[i] == 0.0) || (force->special_lj[i] == 1.0)) + && (force->special_lj[i] != special_lj[istyle][i])) + error->all(FLERR,"Pair_modify special setting incompatible with" + " global special_bonds setting"); + } + } + + if (special_coul[istyle]) { + for (i = 1; i < 4; ++i) { + if (((force->special_coul[i] == 0.0) + || (force->special_coul[i] == 1.0)) + && (force->special_coul[i] != special_coul[istyle][i])) + error->all(FLERR,"Pair_modify special setting incompatible with" + "global special_bonds setting"); + } + } + } + // each sub-style makes its neighbor list request(s) for (istyle = 0; istyle < nstyles; istyle++) styles[istyle]->init_style(); @@ -610,6 +670,13 @@ void PairHybrid::write_restart(FILE *fp) fwrite(&n,sizeof(int),1,fp); fwrite(keywords[m],sizeof(char),n,fp); styles[m]->write_restart_settings(fp); + // write out per style special settings, if present + n = (special_lj[m] == NULL) ? 0 : 1; + fwrite(&n,sizeof(int),1,fp); + if (n) fwrite(special_lj[m],sizeof(double),4,fp); + n = (special_coul[m] == NULL) ? 0 : 1; + fwrite(&n,sizeof(int),1,fp); + if (n) fwrite(special_coul[m],sizeof(double),4,fp); } } @@ -629,6 +696,9 @@ void PairHybrid::read_restart(FILE *fp) keywords = new char*[nstyles]; multiple = new int[nstyles]; + special_lj = new double*[nstyles]; + special_coul = new double*[nstyles]; + // each sub-style is created via new_pair() // each reads its settings, but no coeff info @@ -641,6 +711,22 @@ void PairHybrid::read_restart(FILE *fp) MPI_Bcast(keywords[m],n,MPI_CHAR,0,world); styles[m] = force->new_pair(keywords[m],0,dummy); styles[m]->read_restart_settings(fp); + // read back per style special settings, if present + special_lj[m] = special_coul[m] = NULL; + if (me == 0) fread(&n,sizeof(int),1,fp); + MPI_Bcast(&n,1,MPI_INT,0,world); + if (n > 0 ) { + special_lj[m] = new double[4]; + if (me == 0) fread(special_lj[m],sizeof(double),4,fp); + MPI_Bcast(special_lj[m],4,MPI_DOUBLE,0,world); + } + if (me == 0) fread(&n,sizeof(int),1,fp); + MPI_Bcast(&n,1,MPI_INT,0,world); + if (n > 0 ) { + special_coul[m] = new double[4]; + if (me == 0) fread(special_coul[m],sizeof(double),4,fp); + MPI_Bcast(special_coul[m],4,MPI_DOUBLE,0,world); + } } // multiple[i] = 1 to M if sub-style used multiple times, else 0 @@ -681,6 +767,11 @@ double PairHybrid::single(int i, int j, int itype, int jtype, if (styles[map[itype][jtype][m]]->single_enable == 0) error->one(FLERR,"Pair hybrid sub-style does not support single call"); + if ((special_lj[map[itype][jtype][m]] != NULL) || + (special_coul[map[itype][jtype][m]] != NULL)) + error->one(FLERR,"Pair hybrid single calls do not support" + " per sub-style special bond values"); + esum += styles[map[itype][jtype][m]]-> single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fone); fforce += fone; @@ -716,8 +807,14 @@ void PairHybrid::modify_params(int narg, char **arg) if (strcmp(arg[1],keywords[m]) == 0) break; if (m == nstyles) error->all(FLERR,"Unknown pair_modify hybrid sub-style"); if (multiple[m] == 0) { - Pair::modify_params(narg-2,&arg[2]); - styles[m]->modify_params(narg-2,&arg[2]); + // augment special settings for this one pair style + if (strcmp(arg[2],"special") == 0) { + if (narg < 7) error->all(FLERR,"Illegal pair_modify special command"); + modify_special(m,narg-3,&arg[3]); + } else { + Pair::modify_params(narg-2,&arg[2]); + styles[m]->modify_params(narg-2,&arg[2]); + } } else { if (narg < 3) error->all(FLERR,"Illegal pair_modify command"); int multiflag = force->inumeric(FLERR,arg[2]); @@ -725,13 +822,98 @@ void PairHybrid::modify_params(int narg, char **arg) if (strcmp(arg[1],keywords[m]) == 0 && multiflag == multiple[m]) break; if (m == nstyles) error->all(FLERR,"Unknown pair_modify hybrid sub-style"); - Pair::modify_params(narg-3,&arg[3]); - styles[m]->modify_params(narg-3,&arg[3]); + // augment special settings for this one pair style + if (strcmp(arg[3],"special") == 0) { + if (narg < 8) error->all(FLERR,"Illegal pair_modify special command"); + modify_special(m,narg-4,&arg[4]); + } else { + Pair::modify_params(narg-3,&arg[3]); + styles[m]->modify_params(narg-3,&arg[3]); + } } } else { - Pair::modify_params(narg,arg); - for (int m = 0; m < nstyles; m++) styles[m]->modify_params(narg,arg); + // augment special settings for all pair styles + if (strcmp(arg[0],"special") == 0) { + if (narg < 5) error->all(FLERR,"Illegal pair_modify special command"); + for (int m = 0; m < nstyles; m++) modify_special(m,narg-1,&arg[1]); + } else { + Pair::modify_params(narg,arg); + for (int m = 0; m < nstyles; m++) styles[m]->modify_params(narg,arg); + } + } +} + +/* ---------------------------------------------------------------------- + store a local per pair style override for special_lj and special_coul +------------------------------------------------------------------------- */ + +void PairHybrid::modify_special(int m, int narg, char **arg) +{ + double special[4]; + int i; + + special[0] = 1.0; + special[1] = force->numeric(FLERR,arg[1]); + special[2] = force->numeric(FLERR,arg[2]); + special[3] = force->numeric(FLERR,arg[3]); + + if (strcmp(arg[0],"lj/coul") == 0) { + if (!special_lj[m]) special_lj[m] = new double[4]; + if (!special_coul[m]) special_coul[m] = new double[4]; + for (i = 0; i < 4; ++i) + special_lj[m][i] = special_coul[m][i] = special[i]; + + } else if (strcmp(arg[0],"lj") == 0) { + if (!special_lj[m]) special_lj[m] = new double[4]; + for (i = 0; i < 4; ++i) + special_lj[m][i] = special[i]; + + } else if (strcmp(arg[0],"coul") == 0) { + if (!special_coul[m]) special_coul[m] = new double[4]; + for (i = 0; i < 4; ++i) + special_coul[m][i] = special[i]; + + } else error->all(FLERR,"Illegal pair_modify special command"); +} + +/* ---------------------------------------------------------------------- + override global special bonds settings with per substyle values +------------------------------------------------------------------------- */ + +void PairHybrid::set_special(int m) +{ + int i; + if (special_lj[m]) + for (i = 0; i < 4; ++i) force->special_lj[i] = special_lj[m][i]; + if (special_coul[m]) + for (i = 0; i < 4; ++i) force->special_coul[i] = special_coul[m][i]; +} + +/* ---------------------------------------------------------------------- + store global special settings +------------------------------------------------------------------------- */ + +double * PairHybrid::save_special() +{ + double *saved = new double[8]; + + for (int i = 0; i < 4; ++i) { + saved[i] = force->special_lj[i]; + saved[i+4] = force->special_coul[i]; + } + return saved; +} + +/* ---------------------------------------------------------------------- + restore global special settings from saved data +------------------------------------------------------------------------- */ + +void PairHybrid::restore_special(double *saved) +{ + for (int i = 0; i < 4; ++i) { + force->special_lj[i] = saved[i]; + force->special_coul[i] = saved[i+4]; } } diff --git a/src/pair_hybrid.h b/src/pair_hybrid.h index 9fd29673ee..bc41b993e2 100644 --- a/src/pair_hybrid.h +++ b/src/pair_hybrid.h @@ -26,12 +26,11 @@ PairStyle(hybrid,PairHybrid) namespace LAMMPS_NS { class PairHybrid : public Pair { + friend class FixGPU; + friend class FixOMP; + friend class Force; + friend class Respa; public: - int nstyles; // # of sub-styles - Pair **styles; // list of Pair style classes - char **keywords; // style name of each Pair style - int *multiple; // 0 if style used once, else Mth instance - PairHybrid(class LAMMPS *); virtual ~PairHybrid(); void compute(int, int); @@ -54,13 +53,28 @@ class PairHybrid : public Pair { int check_ijtype(int, int, char *); protected: + int nstyles; // # of sub-styles + Pair **styles; // list of Pair style classes + char **keywords; // style name of each Pair style + int *multiple; // 0 if style used once, else Mth instance + int outerflag; // toggle compute() when invoked by outer() + int respaflag; // 1 if different substyles are assigned to + // different r-RESPA levels int **nmap; // # of sub-styles itype,jtype points to int ***map; // list of sub-styles itype,jtype points to + double **special_lj; // list of per style LJ exclusion factors + double **special_coul; // list of per style Coulomb exclusion factors void allocate(); void flags(); + + void modify_special(int, int, char**); + double *save_special(); + void set_special(int); + void restore_special(double *); + virtual void modify_requests(); }; diff --git a/src/respa.cpp b/src/respa.cpp index 549415f16f..3af1dafcfb 100644 --- a/src/respa.cpp +++ b/src/respa.cpp @@ -38,6 +38,7 @@ #include "timer.h" #include "memory.h" #include "error.h" +#include "pair_hybrid.h" using namespace LAMMPS_NS; @@ -65,6 +66,11 @@ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg) level_pair = level_kspace = -1; level_inner = level_middle = level_outer = -1; + // defaults for hybrid pair styles + nhybrid_styles = 0; + tally_global = 1; + pair_compute = 1; + int iarg = nlevels; while (iarg < narg) { if (strcmp(arg[iarg],"bond") == 0) { @@ -107,6 +113,22 @@ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg) if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command"); level_kspace = force->inumeric(FLERR,arg[iarg+1]) - 1; iarg += 2; + } else if (strcmp(arg[iarg],"hybrid") == 0) { + // the hybrid keyword requires a hybrid pair style + if (!strstr(force->pair_style,"hybrid")) + error->all(FLERR,"Illegal run_style respa command"); + PairHybrid *hybrid = (PairHybrid *) force->pair; + nhybrid_styles = hybrid->nstyles; + // each hybrid sub-style needs to be assigned to a respa level + if (iarg+nhybrid_styles > narg) + error->all(FLERR,"Illegal run_style respa command"); + hybrid_level = new int[nhybrid_styles]; + hybrid_compute = new int[nhybrid_styles]; + for (int i=0; i < nhybrid_styles; ++i) { + ++iarg; + hybrid_level[i] = force->inumeric(FLERR,arg[iarg])-1; + } + ++iarg; } else error->all(FLERR,"Illegal run_style respa command"); } @@ -127,6 +149,10 @@ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg) if (level_middle >= 0 && level_inner == -1) error->all(FLERR,"Cannot set respa middle without inner/outer"); + // cannot combine hybrid with any of pair/inner/middle/outer + if ((nhybrid_styles > 0) && (level_pair >= 0 || level_inner >= 0 + || level_middle >= 0 || level_outer >= 0)) + error->all(FLERR,"Cannot set respa hybrid and any of pair/inner/middle/outer"); // set defaults if user did not specify level // bond to innermost level // angle same as bond, dihedral same as angle, improper same as dihedral @@ -138,9 +164,23 @@ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg) if (level_angle == -1) level_angle = level_bond; if (level_dihedral == -1) level_dihedral = level_angle; if (level_improper == -1) level_improper = level_dihedral; - if (level_pair == -1 && level_inner == -1) level_pair = nlevels-1; + + if (level_pair == -1 && level_inner == -1 && nhybrid_styles < 1) + level_pair = nlevels-1; + if (level_kspace == -1 && level_pair >= 0) level_kspace = level_pair; - if (level_kspace == -1 && level_pair == -1) level_kspace = level_outer; + if (level_kspace == -1 && level_pair == -1) { + if (nhybrid_styles < 1) { + level_kspace = level_outer; + } else { + int max_hybrid_level = -1; + for (int i=0; i < nhybrid_styles; ++i) { + if (max_hybrid_level < hybrid_level[i]) + max_hybrid_level = hybrid_level[i]; + } + level_kspace = max_hybrid_level; + } + } // print respa levels @@ -157,6 +197,9 @@ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg) if (level_inner == i) fprintf(screen," pair-inner"); if (level_middle == i) fprintf(screen," pair-middle"); if (level_outer == i) fprintf(screen," pair-outer"); + for (int j=0;jall(FLERR,"Invalid order of forces within respa levels"); } - if (level_pair == -1 && level_middle == -1) { + if (level_pair == -1 && level_middle == -1 && nhybrid_styles < 1) { if (level_inner < level_improper || level_outer < level_inner || level_kspace < level_outer) error->all(FLERR,"Invalid order of forces within respa levels"); @@ -223,6 +269,12 @@ Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg) cutoff[3] = cutoff[1]; } + // ensure that pair->compute() is run properly when the "hybrid" keyword is not used. + if (nhybrid_styles < 1) { + pair_compute = 1; + tally_global = 1; + } + // allocate other needed arrays newton = new int[nlevels]; @@ -236,6 +288,10 @@ Respa::~Respa() delete [] loop; delete [] newton; delete [] step; + if (nhybrid_styles > 0) { + delete [] hybrid_level; + delete [] hybrid_compute; + } } /* ---------------------------------------------------------------------- @@ -254,13 +310,19 @@ void Respa::init() // create fix needed for storing atom-based respa level forces // will delete it at end of run - char **fixarg = new char*[4]; + char **fixarg = new char*[5]; fixarg[0] = (char *) "RESPA"; fixarg[1] = (char *) "all"; fixarg[2] = (char *) "RESPA"; fixarg[3] = new char[8]; sprintf(fixarg[3],"%d",nlevels); - modify->add_fix(4,fixarg); + // if supported, we also store torques on a per-level basis + if (atom->torque_flag) { + fixarg[4] = (char *) "torque"; + modify->add_fix(5,fixarg); + } else { + modify->add_fix(4,fixarg); + } delete [] fixarg[3]; delete [] fixarg; fix_respa = (FixRespa *) modify->fix[modify->nfix-1]; @@ -309,6 +371,11 @@ void Respa::init() if (level_pair == ilevel || level_inner == ilevel || level_middle == ilevel || level_outer == ilevel) newton[ilevel] = 1; + + if (nhybrid_styles > 0) { + set_compute_flags(ilevel); + if (pair_compute) newton[ilevel] = 1; + } } } @@ -360,6 +427,11 @@ void Respa::setup() for (int ilevel = 0; ilevel < nlevels; ilevel++) { force_clear(newton[ilevel]); modify->setup_pre_force_respa(vflag,ilevel); + + if (nhybrid_styles > 0) { + set_compute_flags(ilevel); + force->pair->compute(eflag,vflag); + } if (level_pair == ilevel && pair_compute_flag) force->pair->compute(eflag,vflag); if (level_inner == ilevel && pair_compute_flag) @@ -428,6 +500,12 @@ void Respa::setup_minimal(int flag) for (int ilevel = 0; ilevel < nlevels; ilevel++) { force_clear(newton[ilevel]); modify->setup_pre_force_respa(vflag,ilevel); + + if (nhybrid_styles > 0) { + set_compute_flags(ilevel); + force->pair->compute(eflag,vflag); + } + if (level_pair == ilevel && pair_compute_flag) force->pair->compute(eflag,vflag); if (level_inner == ilevel && pair_compute_flag) @@ -574,6 +652,11 @@ void Respa::recurse(int ilevel) modify->pre_force_respa(vflag,ilevel,iloop); timer->stamp(); + if (nhybrid_styles > 0) { + set_compute_flags(ilevel); + force->pair->compute(eflag,vflag); + timer->stamp(TIME_PAIR); + } if (level_pair == ilevel && pair_compute_flag) { force->pair->compute(eflag,vflag); timer->stamp(TIME_PAIR); @@ -654,12 +737,19 @@ void Respa::copy_f_flevel(int ilevel) { double ***f_level = fix_respa->f_level; double **f = atom->f; + double ***t_level = fix_respa->t_level; + double **t = atom->torque; int n = atom->nlocal; for (int i = 0; i < n; i++) { f_level[i][ilevel][0] = f[i][0]; f_level[i][ilevel][1] = f[i][1]; f_level[i][ilevel][2] = f[i][2]; + if (fix_respa->store_torque) { + t_level[i][ilevel][0] = t[i][0]; + t_level[i][ilevel][1] = t[i][1]; + t_level[i][ilevel][2] = t[i][2]; + } } } @@ -671,12 +761,19 @@ void Respa::copy_flevel_f(int ilevel) { double ***f_level = fix_respa->f_level; double **f = atom->f; + double ***t_level = fix_respa->t_level; + double **t = atom->torque; int n = atom->nlocal; for (int i = 0; i < n; i++) { f[i][0] = f_level[i][ilevel][0]; f[i][1] = f_level[i][ilevel][1]; f[i][2] = f_level[i][ilevel][2]; + if (fix_respa->store_torque) { + t[i][0] = t_level[i][ilevel][0]; + t[i][1] = t_level[i][ilevel][1]; + t[i][2] = t_level[i][ilevel][2]; + } } } @@ -690,6 +787,8 @@ void Respa::sum_flevel_f() double ***f_level = fix_respa->f_level; double **f = atom->f; + double ***t_level = fix_respa->t_level; + double **t = atom->torque; int n = atom->nlocal; for (int ilevel = 1; ilevel < nlevels; ilevel++) { @@ -697,6 +796,28 @@ void Respa::sum_flevel_f() f[i][0] += f_level[i][ilevel][0]; f[i][1] += f_level[i][ilevel][1]; f[i][2] += f_level[i][ilevel][2]; + if (fix_respa->store_torque) { + t[i][0] += t_level[i][ilevel][0]; + t[i][1] += t_level[i][ilevel][1]; + t[i][2] += t_level[i][ilevel][2]; + } } } } + +/*----------------------------------------------------------------------- + set flags for when some hybrid forces should be computed +------------------------------------------------------------------------- */ + +void Respa::set_compute_flags(int ilevel) +{ + + if (nhybrid_styles < 1) return; + + pair_compute = 0; + for (int i=0; i