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