diff --git a/src/USER-DRUDE/compute_temp_drude.cpp b/src/USER-DRUDE/compute_temp_drude.cpp index 99578e65e0..44b727f6b7 100644 --- a/src/USER-DRUDE/compute_temp_drude.cpp +++ b/src/USER-DRUDE/compute_temp_drude.cpp @@ -37,6 +37,7 @@ ComputeTempDrude::ComputeTempDrude(LAMMPS *lmp, int narg, char **arg) : if (narg != 3) error->all(FLERR,"Illegal compute temp command"); vector_flag = 1; + scalar_flag = 1; size_vector = 6; extscalar = 0; extvector = -1; @@ -216,3 +217,9 @@ void ComputeTempDrude::compute_vector() vector[5] = kineng_drude; } +double ComputeTempDrude::compute_scalar(){ + compute_vector(); + scalar = vector[0]; + return scalar; +} + diff --git a/src/USER-DRUDE/compute_temp_drude.h b/src/USER-DRUDE/compute_temp_drude.h index 834d1d7519..88d767f8b7 100644 --- a/src/USER-DRUDE/compute_temp_drude.h +++ b/src/USER-DRUDE/compute_temp_drude.h @@ -32,6 +32,7 @@ class ComputeTempDrude : public Compute { void init(); void setup(); void compute_vector(); + double compute_scalar(); int modify_param(int, char **); private: diff --git a/src/USER-DRUDE/fix_drude.cpp b/src/USER-DRUDE/fix_drude.cpp index 7f37b89eff..77253d1bb1 100644 --- a/src/USER-DRUDE/fix_drude.cpp +++ b/src/USER-DRUDE/fix_drude.cpp @@ -19,6 +19,8 @@ #include "modify.h" #include "error.h" #include "memory.h" +#include "molecule.h" +#include "atom_vec.h" #include #include @@ -38,6 +40,7 @@ FixDrude::FixDrude(LAMMPS *lmp, int narg, char **arg) : comm_border = 1; // drudeid special_alter_flag = 1; create_attribute = 1; + is_reduced = false; memory->create(drudetype, atom->ntypes+1, "fix_drude::drudetype"); for (int i=3; i[nlocal]; // Temporary sets of bond partner tags sptr = this; - // Build list of my atoms' bond partners - for (int i=0; inum_bond[i]; k++){ - core_drude_vec.push_back(atom->tag[i]); - core_drude_vec.push_back(atom->bond_atom[i][k]); + if (atom->molecular == 1) + { + // Build list of my atoms' bond partners + for (int i=0; inum_bond[i]; k++){ + core_drude_vec.push_back(atom->tag[i]); + core_drude_vec.push_back(atom->bond_atom[i][k]); + } + } + } + else + { + // Template case + class Molecule **atommols; + atommols = atom->avec->onemols; + + // Build list of my atoms' bond partners + for (int i=0; imolindex[i]; + int iatom = atom->molatom[i]; + tagint *batom = atommols[imol]->bond_atom[iatom]; + tagint tagprev = atom->tag[i] - iatom - 1; + int nbonds = atommols[imol]->num_bond[iatom]; + + if (drudetype[type[i]] == NOPOL_TYPE) continue; + drudeid[i] = 0; + for (int k=0; ktag[i]); + core_drude_vec.push_back(batom[k]+tagprev); + } } } // Loop on procs to fill my atoms' sets of bond partners @@ -275,6 +303,9 @@ void FixDrude::rebuild_special(){ tagint **special = atom->special; int *type = atom->type; + if (atom->molecular != 1) + return; + // Make sure that drude partners know each other //build_drudeid(); diff --git a/src/USER-DRUDE/fix_langevin_drude.cpp b/src/USER-DRUDE/fix_langevin_drude.cpp index 119328a1dd..319516aa7d 100644 --- a/src/USER-DRUDE/fix_langevin_drude.cpp +++ b/src/USER-DRUDE/fix_langevin_drude.cpp @@ -44,13 +44,9 @@ FixLangevinDrude::FixLangevinDrude(LAMMPS *lmp, int narg, char **arg) : // Langevin thermostat should be applied every step nevery = 1; - vector_flag = 1; global_freq = nevery; - extvector = 0; - size_vector = 6; comm_reverse = 3; - //extscalar = 1; - + // core temperature tstr_core = NULL; if (strstr(arg[3],"v_") == arg[3]) { diff --git a/src/USER-DRUDE/pair_lj_cut_thole_long.cpp b/src/USER-DRUDE/pair_lj_cut_thole_long.cpp index ef3747e7b3..efa11c6c5e 100644 --- a/src/USER-DRUDE/pair_lj_cut_thole_long.cpp +++ b/src/USER-DRUDE/pair_lj_cut_thole_long.cpp @@ -577,7 +577,7 @@ void PairLJCutTholeLong::read_restart_settings(FILE *fp) void PairLJCutTholeLong::write_data(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) - fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]); + fprintf(fp,"%d %g %g %g %g\n",i,epsilon[i][i],sigma[i][i],polar[i][i],thole[i][i]); } /* ---------------------------------------------------------------------- @@ -588,7 +588,7 @@ void PairLJCutTholeLong::write_data_all(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) - fprintf(fp,"%d %d %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],cut_lj[i][j]); + fprintf(fp,"%d %d %g %g %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],polar[i][j],thole[i][j],cut_lj[i][j]); } /* ---------------------------------------------------------------------- */