Repository navigation
Expand file tree
/
Copy paththermo.cpp
More file actions
250 lines (205 loc) · 6.13 KB
/
Copy paththermo.cpp
File metadata and controls
250 lines (205 loc) · 6.13 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
/* ----------------------------------------------------------------------
miniMD is a simple, parallel molecular dynamics (MD) code. miniMD is
an MD microapplication in the Mantevo project at Sandia National
Laboratories ( http://www.mantevo.org ). The primary
authors of miniMD are Steve Plimpton (sjplimp@sandia.gov), Paul Crozier
(pscrozi@sandia.gov) and Christian Trott (crtrott@sandia.gov).
Copyright (2008) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This library is free software; you
can redistribute it and/or modify it under the terms of the GNU Lesser
General Public License as published by the Free Software Foundation;
either version 3 of the License, or (at your option) any later
version.
This library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this software; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
USA. See also: http://www.gnu.org/licenses/lgpl.txt .
For questions, contact Paul S. Crozier (pscrozi@sandia.gov) or
Christian Trott (crtrott@sandia.gov).
Please read the accompanying README and LICENSE files.
---------------------------------------------------------------------- */
#include "thermo.h"
#include "force.h"
#include "integrate.h"
#include "mpi.h"
#include "offload.h"
#include "stdio.h"
#include "stdlib.h"
#include "util.h"
#include <omp.h>
Thermo::Thermo() {}
Thermo::~Thermo() {}
void Thermo::setup(MMD_float rho_in, Integrate &integrate, Atom &atom, int units)
{
rho = rho_in;
ntimes = integrate.ntimes;
MMD_int maxstat;
if(nstat == 0)
{
maxstat = 2;
}
else
{
maxstat = ntimes / nstat + 2;
}
steparr = ( MMD_int * )malloc(maxstat * sizeof(MMD_int));
tmparr = ( MMD_float * )malloc(maxstat * sizeof(MMD_float));
engarr = ( MMD_float * )malloc(maxstat * sizeof(MMD_float));
prsarr = ( MMD_float * )malloc(maxstat * sizeof(MMD_float));
if(units == LJ)
{
mvv2e = 1.0;
dof_boltz = (atom.natoms * 3 - 3);
t_scale = mvv2e / dof_boltz;
p_scale = 1.0 / 3 / atom.box.xprd / atom.box.yprd / atom.box.zprd;
e_scale = 0.5;
}
else if(units == METAL)
{
mvv2e = 1.036427e-04;
dof_boltz = (atom.natoms * 3 - 3) * 8.617343e-05;
t_scale = mvv2e / dof_boltz;
p_scale = 1.602176e+06 / 3 / atom.box.xprd / atom.box.yprd / atom.box.zprd;
e_scale = 524287.985533; // 16.0;
integrate.dtforce /= mvv2e;
}
}
void Thermo::compute(MMD_int iflag, Atom &atom, Neighbor &neighbor, Force *force, Timer &timer, Comm &comm)
{
MMD_float t, eng, p;
if(iflag > 0 && iflag % nstat)
{
return;
}
if(iflag == -1 && nstat > 0 && ntimes % nstat == 0)
{
return;
}
t_act = 0;
e_act = 0;
p_act = 0;
t = temperature(atom);
eng = energy(atom, neighbor, force);
p = pressure(t, force);
MMD_int istep = iflag;
if(iflag == -1)
{
istep = ntimes;
}
if(iflag == 0)
{
mstat = 0;
}
steparr[mstat] = istep;
tmparr[mstat] = t;
engarr[mstat] = eng;
prsarr[mstat] = p;
mstat++;
double oldtime = timer.array[TIME_TOTAL];
timer.barrier_stop(TIME_TOTAL);
if(threads->mpi_me == 0)
{
fprintf(stdout, "%i %e %e %e %6.3lf\n", istep, t, eng, p, istep == 0 ? 0.0 : timer.array[TIME_TOTAL]);
}
timer.array[TIME_TOTAL] = oldtime;
}
/* reduced potential energy */
MMD_float Thermo::energy(Atom &atom, Neighbor &neighbor, Force *force)
{
e_act = force->eng_vdwl;
if(neighbor.halfneigh)
{
e_act *= 2.0;
}
e_act *= e_scale;
MMD_float eng;
if(sizeof(MMD_float) == 4)
{
MPI_Allreduce(&e_act, &eng, 1, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD);
}
else
{
MPI_Allreduce(&e_act, &eng, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
}
return eng / atom.natoms;
}
/* reduced temperature */
MMD_float Thermo::temperature(Atom &atom)
{
t_act = 0;
const int nlocal = atom.nlocal;
const MMD_float mass = atom.mass;
MMD_float * v = atom.v;
MMD_float t = 0.0;
#ifdef USE_OFFLOAD
#pragma omp target teams map(tofrom:t) thread_limit(MAX_TEAM_SIZE)
#endif
{
#ifdef USE_OFFLOAD
const int nthr = omp_get_max_threads();
MMD_float w_t[MAX_TEAM_SIZE];
for(int i = 0; i < nthr; i++)
{
w_t[i] = MMD_float(0.0);
}
#pragma omp distribute parallel for num_threads(nthr) // FIXME: this is because for x86, it appears that extra threads can run in here.
#else
#pragma omp parallel for reduction(+:t)
#endif
for(MMD_int i = 0; i < nlocal; i++)
{
const int tid = omp_get_thread_num();
const MMD_float vx = v[i * PAD + 0];
const MMD_float vy = v[i * PAD + 1];
const MMD_float vz = v[i * PAD + 2];
const MMD_float l_t = (vx * vx + vy * vy + vz * vz) * mass;
#ifdef USE_OFFLOAD
w_t[tid] += l_t;
#else
t += l_t;
#endif
}
#ifdef USE_OFFLOAD
MMD_float team_t = MMD_float(0.0);
for(int i = 0; i < nthr; i++)
{
team_t += w_t[i];
}
atomic_add(&t, team_t);
#endif
}
t_act += t;
MMD_float t1;
if(sizeof(MMD_float) == 4)
{
MPI_Allreduce(&t_act, &t1, 1, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD);
}
else
{
MPI_Allreduce(&t_act, &t1, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
}
return t1 * t_scale;
}
/* reduced pressure from virial
virial = Fi dot Ri summed over own and ghost atoms, since PBC info is
stored correctly in force array before reverse_communicate is performed */
MMD_float Thermo::pressure(MMD_float t, Force *force)
{
p_act = force->virial;
MMD_float virial = 0;
if(sizeof(MMD_float) == 4)
{
MPI_Allreduce(&p_act, &virial, 1, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD);
}
else
{
MPI_Allreduce(&p_act, &virial, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
}
// printf("Pres: %e %e %e %e\n",t,dof_boltz,virial,p_scale);
return (t * dof_boltz + virial) * p_scale;
}