1 // clang-format off
2 /* ----------------------------------------------------------------------
3    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4    https://www.lammps.org/, Sandia National Laboratories
5    Steve Plimpton, sjplimp@sandia.gov
6 
7    Copyright (2003) Sandia Corporation.  Under the terms of Contract
8    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9    certain rights in this software.  This software is distributed under
10    the GNU General Public License.
11 
12    See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14 
15 #include "compute_temp_drude.h"
16 
17 #include "atom.h"
18 #include "comm.h"
19 #include "domain.h"
20 #include "error.h"
21 #include "fix_drude.h"
22 #include "force.h"
23 #include "modify.h"
24 #include "update.h"
25 
26 #include <cstring>
27 
28 using namespace LAMMPS_NS;
29 
30 /* ---------------------------------------------------------------------- */
31 
ComputeTempDrude(LAMMPS * lmp,int narg,char ** arg)32 ComputeTempDrude::ComputeTempDrude(LAMMPS *lmp, int narg, char **arg) :
33   Compute(lmp, narg, arg)
34 {
35   if (narg != 3) error->all(FLERR,"Illegal compute temp command");
36 
37   vector_flag = 1;
38   scalar_flag = 1;
39   size_vector = 6;
40   extscalar = 0;
41   extvector = -1;
42   extlist = new int[6];
43   extlist[0] = extlist[1] = 0;
44   extlist[2] = extlist[3] = extlist[4] = extlist[5] = 1;
45   tempflag = 0; // because does not compute a single temperature (scalar and vector)
46 
47   vector = new double[size_vector];
48   fix_drude = nullptr;
49   id_temp = nullptr;
50   temperature = nullptr;
51 }
52 
53 /* ---------------------------------------------------------------------- */
54 
~ComputeTempDrude()55 ComputeTempDrude::~ComputeTempDrude()
56 {
57   delete [] vector;
58   delete [] extlist;
59   delete [] id_temp;
60 }
61 
62 /* ---------------------------------------------------------------------- */
63 
init()64 void ComputeTempDrude::init()
65 {
66   int ifix;
67   for (ifix = 0; ifix < modify->nfix; ifix++)
68     if (strcmp(modify->fix[ifix]->style,"drude") == 0) break;
69   if (ifix == modify->nfix) error->all(FLERR, "compute temp/drude requires fix drude");
70   fix_drude = (FixDrude *) modify->fix[ifix];
71 
72   if (!comm->ghost_velocity)
73     error->all(FLERR,"compute temp/drude requires ghost velocities. Use comm_modify vel yes");
74 }
75 
76 /* ---------------------------------------------------------------------- */
77 
setup()78 void ComputeTempDrude::setup()
79 {
80   dof_compute();
81 }
82 
83 /* ---------------------------------------------------------------------- */
84 
dof_compute()85 void ComputeTempDrude::dof_compute()
86 {
87   int nlocal = atom->nlocal;
88   int *type = atom->type;
89   int dim = domain->dimension;
90   int *drudetype = fix_drude->drudetype;
91 
92   fix_dof = 0;
93   for (int i = 0; i < modify->nfix; i++)
94     fix_dof += modify->fix[i]->dof(igroup);
95 
96   bigint dof_core_loc = 0, dof_drude_loc = 0;
97   for (int i = 0; i < nlocal; i++) {
98     if (atom->mask[i] & groupbit) {
99       if (drudetype[type[i]] == DRUDE_TYPE) // Non-polarizable atom
100           dof_drude_loc++;
101       else
102           dof_core_loc++;
103     }
104   }
105   dof_core_loc *= dim;
106   dof_drude_loc *= dim;
107   MPI_Allreduce(&dof_core_loc,  &dof_core,  1, MPI_LMP_BIGINT, MPI_SUM, world);
108   MPI_Allreduce(&dof_drude_loc, &dof_drude, 1, MPI_LMP_BIGINT, MPI_SUM, world);
109   dof_core -= fix_dof;
110   vector[2] = dof_core;
111   vector[3] = dof_drude;
112 }
113 
114 /* ---------------------------------------------------------------------- */
115 
modify_param(int narg,char ** arg)116 int ComputeTempDrude::modify_param(int narg, char **arg)
117 {
118   if (strcmp(arg[0],"temp") == 0) {
119     if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
120     delete [] id_temp;
121     id_temp = utils::strdup(arg[1]);
122 
123     int icompute = modify->find_compute(id_temp);
124     if (icompute < 0)
125       error->all(FLERR,"Could not find fix_modify temperature ID");
126     temperature = modify->compute[icompute];
127 
128     if (temperature->tempflag == 0)
129       error->all(FLERR,
130                  "Fix_modify temperature ID does not compute temperature");
131     if (temperature->igroup != igroup && comm->me == 0)
132       error->warning(FLERR,"Group for fix_modify temp != fix group");
133     return 2;
134   }
135   return 0;
136 }
137 
138 /* ---------------------------------------------------------------------- */
139 
compute_vector()140 void ComputeTempDrude::compute_vector()
141 {
142     invoked_vector = update->ntimestep;
143 
144     int nlocal = atom->nlocal;
145     int *mask = atom->mask;
146     int *type = atom->type;
147     double *rmass = atom->rmass, *mass = atom->mass;
148     double **v = atom->v;
149     tagint *drudeid = fix_drude->drudeid;
150     int *drudetype = fix_drude->drudetype;
151     int dim = domain->dimension;
152     double mvv2e = force->mvv2e, kb = force->boltz;
153 
154     double mcore, mdrude;
155     double ecore, edrude;
156     double *vcore, *vdrude;
157     double kineng_core_loc = 0., kineng_drude_loc = 0.;
158     for (int i=0; i<nlocal; i++) {
159         if (groupbit & mask[i] && drudetype[type[i]] != DRUDE_TYPE) {
160             if (drudetype[type[i]] == NOPOL_TYPE) {
161                 ecore = 0.;
162                 vcore = v[i];
163                 if (temperature) temperature->remove_bias(i, vcore);
164                 for (int k=0; k<dim; k++) ecore += vcore[k]*vcore[k];
165                 if (temperature) temperature->restore_bias(i, vcore);
166                 if (rmass) mcore = rmass[i];
167                 else mcore = mass[type[i]];
168                 kineng_core_loc += mcore * ecore;
169             } else { // CORE_TYPE
170                 int j = atom->map(drudeid[i]);
171                 if (rmass) {
172                     mcore = rmass[i];
173                     mdrude = rmass[j];
174                 } else {
175                     mcore = mass[type[i]];
176                     mdrude = mass[type[j]];
177                 }
178                 double mtot_inv = 1. / (mcore + mdrude);
179                 ecore = 0.;
180                 edrude = 0.;
181                 vcore = v[i];
182                 vdrude = v[j];
183                 if (temperature) {
184                     temperature->remove_bias(i, vcore);
185                     temperature->remove_bias(j, vdrude);
186                 }
187                 for (int k=0; k<dim; k++) {
188                     double v1 = mdrude * vdrude[k] + mcore * vcore[k];
189                     ecore += v1 * v1;
190                     double v2 = vdrude[k] - vcore[k];
191                     edrude += v2 * v2;
192                 }
193                 if (temperature) {
194                     temperature->restore_bias(i, vcore);
195                     temperature->restore_bias(j, vdrude);
196                 }
197                 kineng_core_loc += mtot_inv * ecore;
198                 kineng_drude_loc += mtot_inv * mcore * mdrude * edrude;
199             }
200         }
201     }
202 
203     if (dynamic) dof_compute();
204     kineng_core_loc *= 0.5 * mvv2e;
205     kineng_drude_loc *= 0.5 * mvv2e;
206     MPI_Allreduce(&kineng_core_loc,&kineng_core,1,MPI_DOUBLE,MPI_SUM,world);
207     MPI_Allreduce(&kineng_drude_loc,&kineng_drude,1,MPI_DOUBLE,MPI_SUM,world);
208     temp_core = 2.0 * kineng_core / (dof_core * kb);
209     temp_drude = 2.0 * kineng_drude / (dof_drude * kb);
210     vector[0] = temp_core;
211     vector[1] = temp_drude;
212     vector[4] = kineng_core;
213     vector[5] = kineng_drude;
214 }
215 
compute_scalar()216 double ComputeTempDrude::compute_scalar() {
217     compute_vector();
218     scalar = vector[0];
219     return scalar;
220 }
221 
222