From 493f9c469b7adecfc0a737facd3067c1becbbe3c Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 9 Dec 2015 21:44:39 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14331
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/USER-DRUDE/compute_temp_drude.cpp     |  7 ++++
 src/USER-DRUDE/compute_temp_drude.h       |  1 +
 src/USER-DRUDE/fix_drude.cpp              | 45 +++++++++++++++++++----
 src/USER-DRUDE/fix_langevin_drude.cpp     |  6 +--
 src/USER-DRUDE/pair_lj_cut_thole_long.cpp |  4 +-
 5 files changed, 49 insertions(+), 14 deletions(-)

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 <set>
 #include <vector>
@@ -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<narg; i++) {
@@ -109,13 +112,38 @@ void FixDrude::build_drudeid(){
   partner_set = new std::set<tagint>[nlocal]; // Temporary sets of bond partner tags
 
   sptr = this;
-  // Build list of my atoms' bond partners
-  for (int i=0; i<nlocal; i++){
-    if (drudetype[type[i]] == NOPOL_TYPE) continue;
-    drudeid[i] = 0;
-    for (int k=0; k<atom->num_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; i<nlocal; i++){
+      if (drudetype[type[i]] == NOPOL_TYPE) continue;
+      drudeid[i] = 0;
+      for (int k=0; k<atom->num_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; i<nlocal; i++){
+      int imol = atom->molindex[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; k<nbonds; k++){
+        core_drude_vec.push_back(atom->tag[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]);
 }
 
 /* ---------------------------------------------------------------------- */