forked from gromacs/gromacs
-
Notifications
You must be signed in to change notification settings - Fork 1
/
patch-v2021.2
316 lines (305 loc) · 14.4 KB
/
patch-v2021.2
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
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
patch -p 1 << EOF
diff --git a/src/gromacs/fileio/tpxio.cpp b/src/gromacs/fileio/tpxio.cpp
index 2d87034417..a6c62c2b8f 100644
--- a/src/gromacs/fileio/tpxio.cpp
+++ b/src/gromacs/fileio/tpxio.cpp
@@ -1339,9 +1339,13 @@ static void do_inputrec(gmx::ISerializer* serializer, t_inputrec* ir, int file_v
serializer->doRvec(&ir->compress[XX]);
serializer->doRvec(&ir->compress[YY]);
serializer->doRvec(&ir->compress[ZZ]);
+ serializer->doRvec(&ir->h_start[XX]);
+ serializer->doRvec(&ir->h_start[YY]);
+ serializer->doRvec(&ir->h_start[ZZ]);
serializer->doInt(&ir->refcoord_scaling);
serializer->doRvec(&ir->posres_com);
serializer->doRvec(&ir->posres_comB);
+ serializer->doBool(&ir->deviatoric_crescale);
if (file_version < 79)
{
diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp
index 2ca9db0586..d18139e9e0 100644
--- a/src/gromacs/gmxpreprocess/readir.cpp
+++ b/src/gromacs/gmxpreprocess/readir.cpp
@@ -1915,8 +1915,8 @@ void get_ir(const char* mdparin,
WriteMdpHeader writeMdpHeader,
warninp_t wi)
{
- char* dumstr[2];
- double dumdub[2][6];
+ char* dumstr[3];
+ double dumdub[3][6];
int i, j, m;
char warn_buf[STRLEN];
t_lambda* fep = ir->fepvals;
@@ -1930,6 +1930,7 @@ void get_ir(const char* mdparin,
snew(dumstr[0], STRLEN);
snew(dumstr[1], STRLEN);
+ snew(dumstr[2], STRLEN);
/* ignore the following deprecated commands */
replace_inp_entry(inp, "title", nullptr);
@@ -2148,6 +2149,8 @@ void get_ir(const char* mdparin,
ir->tau_p = get_ereal(&inp, "tau-p", 1.0, wi);
setStringEntry(&inp, "compressibility", dumstr[0], nullptr);
setStringEntry(&inp, "ref-p", dumstr[1], nullptr);
+ printStringNoNewline(&inp, "Initial box (nm)");
+ setStringEntry(&inp, "h-start", dumstr[2], nullptr);
printStringNoNewline(&inp, "Scaling of reference coordinates, No, All or COM");
ir->refcoord_scaling = get_eeenum(&inp, "refcoord-scaling", erefscaling_names, wi);
@@ -2502,6 +2505,9 @@ void get_ir(const char* mdparin,
}
/* Process options if necessary */
+ ir->deviatoric_crescale = false;
+ dumdub[2][XX] = dumdub[2][YY] = dumdub[2][ZZ] = 0.0;
+ dumdub[2][3] = dumdub[2][4] = dumdub[2][5] = 0.0;
for (m = 0; m < 2; m++)
{
for (i = 0; i < 2 * DIM; i++)
@@ -2540,6 +2546,20 @@ void get_ir(const char* mdparin,
wi,
"Pressure coupling incorrect number of values (I need exactly 6)");
}
+ // Deviatoric stress with C-rescale
+ if (m == 1 && (ir->epc == epcCRESCALE) && (dumdub[1][XX] != dumdub[1][YY] ||
+ dumdub[1][XX] != dumdub[1][ZZ] || dumdub[1][YY] != dumdub[1][ZZ] ||
+ dumdub[1][3] != 0.0 || dumdub[1][4] != 0.0 || dumdub[1][5] != 0.0))
+ {
+ ir->deviatoric_crescale = true;
+ if (sscanf(dumstr[2], "%lf%lf%lf%lf%lf%lf", &(dumdub[2][XX]), &(dumdub[2][YY]),
+ &(dumdub[2][ZZ]), &(dumdub[2][3]), &(dumdub[2][4]), &(dumdub[2][5])) != 6)
+ {
+ warning_error(
+ wi,
+ "h-start incorrect number of values (I need exactly 6)");
+ }
+ }
break;
default:
gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
@@ -2549,10 +2569,12 @@ void get_ir(const char* mdparin,
}
clear_mat(ir->ref_p);
clear_mat(ir->compress);
+ clear_mat(ir->h_start);
for (i = 0; i < DIM; i++)
{
ir->ref_p[i][i] = dumdub[1][i];
ir->compress[i][i] = dumdub[0][i];
+ ir->h_start[i][i] = dumdub[2][i];
}
if (ir->epct == epctANISOTROPIC)
{
@@ -2568,12 +2590,16 @@ void get_ir(const char* mdparin,
ir->compress[XX][YY] = dumdub[0][3];
ir->compress[XX][ZZ] = dumdub[0][4];
ir->compress[YY][ZZ] = dumdub[0][5];
+ ir->h_start[YY][XX] = dumdub[2][3];
+ ir->h_start[ZZ][XX] = dumdub[2][4];
+ ir->h_start[ZZ][YY] = dumdub[2][5];
for (i = 0; i < DIM; i++)
{
for (m = 0; m < i; m++)
{
ir->ref_p[i][m] = ir->ref_p[m][i];
ir->compress[i][m] = ir->compress[m][i];
+ ir->h_start[m][i] = 0.0;
}
}
}
diff --git a/src/gromacs/mdlib/coupling.cpp b/src/gromacs/mdlib/coupling.cpp
index 504ff512ed..25943df3b5 100644
--- a/src/gromacs/mdlib/coupling.cpp
+++ b/src/gromacs/mdlib/coupling.cpp
@@ -41,7 +41,6 @@
#include <cassert>
#include <cmath>
-
#include <algorithm>
#include "gromacs/domdec/domdec_struct.h"
@@ -952,9 +951,11 @@ void crescale_pcoupl(FILE* fplog,
*/
real scalar_pressure = 0;
real xy_pressure = 0;
+ real ref_p_hydro = 0;
for (int d = 0; d < DIM; d++)
{
scalar_pressure += pres[d][d] / DIM;
+ ref_p_hydro += ir->ref_p[d][d] / DIM;
if (d != ZZ)
{
xy_pressure += pres[d][d] / (DIM - 1);
@@ -1028,20 +1029,74 @@ void crescale_pcoupl(FILE* fplog,
+ std::sqrt(2.0 / 3.0 * kt * compressibilityFactor * PRESFAC / vol) * gauss2);
}
break;
+ case epctANISOTROPIC:
+ for (int d1 = 0; d1 < DIM; d1++)
+ {
+ const real compressibilityFactor = ir->compress[d1][d1] * dt / (DIM * ir->tau_p);
+ mu[d1][d1] = 1.0 - compressibilityFactor * (ref_p_hydro - PRESFAC * kt / vol);
+ for (int d2 = 0; d2 < DIM; d2++)
+ {
+ real gauss = normalDist(rng);
+ mu[d1][d2] += compressibilityFactor * pres[d1][d2]
+ + std::sqrt(2.0 * kt * compressibilityFactor * PRESFAC / vol) * gauss;
+ }
+ }
+ if (ir->deviatoric_crescale)
+ {
+ real vol_start = 1.0;
+ for (int d = 0; d < DIM; d++)
+ {
+ vol_start *= ir->h_start[d][d];
+ }
+ matrix h_start_inv;
+ matrix deviatoric_stress;
+ matrix buffer_tensor1;
+ matrix buffer_tensor2;
+ matrix pres_shear;
+ gmx::invertBoxMatrix(ir->h_start, h_start_inv);
+ copy_mat(ir->ref_p, deviatoric_stress);
+ for (int d = 0; d < DIM; d++)
+ {
+ deviatoric_stress[d][d] -= ref_p_hydro;
+ }
+ mmul(deviatoric_stress, h_start_inv, buffer_tensor1);
+ tmmul(h_start_inv, buffer_tensor1, buffer_tensor2);
+ mmul(buffer_tensor2, box, buffer_tensor1);
+ tmmul(box, buffer_tensor1, pres_shear);
+ for (int d1 = 0; d1 < DIM; d1++)
+ {
+ for (int d2 = 0; d2 < DIM; d2++)
+ {
+ const real compressibilityFactor = ir->compress[d1][d2] * dt / (DIM * ir->tau_p);
+ mu[d1][d2] -= compressibilityFactor * pres_shear[d1][d2] * vol_start / vol;
+ }
+ }
+ }
+ break;
default:
gmx_fatal(FARGS, "C-rescale pressure coupling type %s not supported yet\n",
EPCOUPLTYPETYPE(ir->epct));
}
/* To fullfill the orientation restrictions on triclinic boxes
- * we will set mu_yx, mu_zx and mu_zy to 0 and correct
- * the other elements of mu to first order.
+ * we rotate mu applying Gram-Schmidt process on its raws
*/
- mu[YY][XX] += mu[XX][YY];
- mu[ZZ][XX] += mu[XX][ZZ];
- mu[ZZ][YY] += mu[YY][ZZ];
- mu[XX][YY] = 0;
- mu[XX][ZZ] = 0;
- mu[YY][ZZ] = 0;
+ matrix mu_buffer;
+ copy_mat(mu, mu_buffer);
+ mu[XX][XX] = std::sqrt(mu_buffer[XX][XX]*mu_buffer[XX][XX] + mu_buffer[XX][YY]*mu_buffer[XX][YY]
+ + mu_buffer[XX][ZZ]*mu_buffer[XX][ZZ]);
+ mu[YY][XX] = (mu_buffer[XX][XX]*mu_buffer[YY][XX] + mu_buffer[XX][YY]*mu_buffer[YY][YY]
+ + mu_buffer[XX][ZZ]*mu_buffer[YY][ZZ]) / mu[XX][XX];
+ mu[YY][YY] = std::sqrt(mu_buffer[YY][XX]*mu_buffer[YY][XX] + mu_buffer[YY][YY]*mu_buffer[YY][YY]
+ + mu_buffer[YY][ZZ]*mu_buffer[YY][ZZ] - mu[YY][XX]*mu[YY][XX]);
+ mu[ZZ][XX] = (mu_buffer[XX][XX]*mu_buffer[ZZ][XX] + mu_buffer[XX][YY]*mu_buffer[ZZ][YY]
+ + mu_buffer[XX][ZZ]*mu_buffer[ZZ][ZZ]) / mu[XX][XX];
+ mu[ZZ][YY] = (mu_buffer[YY][XX]*mu_buffer[ZZ][XX] + mu_buffer[YY][YY]*mu_buffer[ZZ][YY]
+ + mu_buffer[YY][ZZ]*mu_buffer[ZZ][ZZ] - mu[ZZ][XX]*mu[YY][XX]) / mu[YY][YY];
+ mu[ZZ][ZZ] = std::sqrt(mu_buffer[ZZ][XX]*mu_buffer[ZZ][XX] + mu_buffer[ZZ][YY]*mu_buffer[ZZ][YY]
+ + mu_buffer[ZZ][ZZ]*mu_buffer[ZZ][ZZ] - mu[ZZ][XX]*mu[ZZ][XX] - mu[ZZ][YY]*mu[ZZ][YY]);
+ mu[XX][YY] = 0.0;
+ mu[XX][ZZ] = 0.0;
+ mu[YY][ZZ] = 0.0;
/* Keep track of the work the barostat applies on the system.
* Without constraints force_vir tells us how Epot changes when scaling.
@@ -1127,13 +1182,13 @@ void crescale_pscale(const t_inputrec* ir,
if (!nFreeze[g][XX])
{
x[n][XX] = mu[XX][XX] * x[n][XX] + mu[YY][XX] * x[n][YY] + mu[ZZ][XX] * x[n][ZZ];
- v[n][XX] = inv_mu[XX][XX] * v[n][XX] + inv_mu[YY][XX] * v[n][YY]
- + inv_mu[ZZ][XX] * v[n][ZZ];
+ v[n][XX] = inv_mu[XX][XX] * v[n][XX] + inv_mu[XX][YY] * v[n][YY]
+ + inv_mu[XX][ZZ] * v[n][ZZ];
}
if (!nFreeze[g][YY])
{
x[n][YY] = mu[YY][YY] * x[n][YY] + mu[ZZ][YY] * x[n][ZZ];
- v[n][YY] = inv_mu[YY][YY] * v[n][YY] + inv_mu[ZZ][YY] * v[n][ZZ];
+ v[n][YY] = inv_mu[YY][YY] * v[n][YY] + inv_mu[YY][ZZ] * v[n][ZZ];
}
if (!nFreeze[g][ZZ])
{
diff --git a/src/gromacs/mdlib/update_constrain_gpu_impl.cu b/src/gromacs/mdlib/update_constrain_gpu_impl.cu
index 825890ce82..271e615d45 100644
--- a/src/gromacs/mdlib/update_constrain_gpu_impl.cu
+++ b/src/gromacs/mdlib/update_constrain_gpu_impl.cu
@@ -106,6 +106,24 @@ __launch_bounds__(c_maxThreadsPerBlock) __global__
}
}
+__launch_bounds__(c_maxThreadsPerBlock) __global__
+ static void scaleVelocities_kernel(const int numAtoms,
+ float3* __restrict__ gm_v,
+ const ScalingMatrix scalingMatrix)
+{
+ int threadIndex = blockIdx.x * blockDim.x + threadIdx.x;
+ if (threadIndex < numAtoms)
+ {
+ float3 v = gm_v[threadIndex];
+
+ v.x = scalingMatrix.xx * v.x;
+ v.y = scalingMatrix.yx * v.x + scalingMatrix.yy * v.y;
+ v.z = scalingMatrix.zx * v.x + scalingMatrix.zy * v.y + scalingMatrix.zz * v.z;
+
+ gm_v[threadIndex] = v;
+ }
+}
+
void UpdateConstrainGpu::Impl::integrate(GpuEventSynchronizer* fReadyOnDevice,
const real dt,
const bool updateVelocities,
@@ -194,9 +212,9 @@ void UpdateConstrainGpu::Impl::scaleVelocities(const matrix scalingMatrix)
mu.zy = scalingMatrix[ZZ][YY];
const auto kernelArgs = prepareGpuKernelArguments(
- scaleCoordinates_kernel, coordinateScalingKernelLaunchConfig_, &numAtoms_, &d_v_, &mu);
- launchGpuKernel(scaleCoordinates_kernel, coordinateScalingKernelLaunchConfig_, deviceStream_,
- nullptr, "scaleCoordinates_kernel", kernelArgs);
+ scaleVelocities_kernel, coordinateScalingKernelLaunchConfig_, &numAtoms_, &d_v_, &mu);
+ launchGpuKernel(scaleVelocities_kernel, coordinateScalingKernelLaunchConfig_, deviceStream_,
+ nullptr, "scaleVelocities_kernel", kernelArgs);
// TODO: Although this only happens on the pressure coupling steps, this synchronization
// can affect the performance if nstpcouple is small.
deviceStream_.synchronize();
diff --git a/src/gromacs/mdtypes/inputrec.cpp b/src/gromacs/mdtypes/inputrec.cpp
index 05a09523fb..f86cb9b9aa 100644
--- a/src/gromacs/mdtypes/inputrec.cpp
+++ b/src/gromacs/mdtypes/inputrec.cpp
@@ -922,6 +922,7 @@ void pr_inputrec(FILE* fp, int indent, const char* title, const t_inputrec* ir,
PR("tau-p", ir->tau_p);
pr_matrix(fp, indent, "compressibility", ir->compress, bMDPformat);
pr_matrix(fp, indent, "ref-p", ir->ref_p, bMDPformat);
+ pr_matrix(fp, indent, "h-start", ir->h_start, bMDPformat);
PS("refcoord-scaling", EREFSCALINGTYPE(ir->refcoord_scaling));
if (bMDPformat)
@@ -1353,6 +1354,9 @@ void cmp_inputrec(FILE* fp, const t_inputrec* ir1, const t_inputrec* ir2, real f
cmp_rvec(fp, "inputrec->compress(x)", -1, ir1->compress[XX], ir2->compress[XX], ftol, abstol);
cmp_rvec(fp, "inputrec->compress(y)", -1, ir1->compress[YY], ir2->compress[YY], ftol, abstol);
cmp_rvec(fp, "inputrec->compress(z)", -1, ir1->compress[ZZ], ir2->compress[ZZ], ftol, abstol);
+ cmp_rvec(fp, "inputrec->h_start(x)", -1, ir1->h_start[XX], ir2->h_start[XX], ftol, abstol);
+ cmp_rvec(fp, "inputrec->h_start(y)", -1, ir1->h_start[YY], ir2->h_start[YY], ftol, abstol);
+ cmp_rvec(fp, "inputrec->h_start(z)", -1, ir1->h_start[ZZ], ir2->h_start[ZZ], ftol, abstol);
cmp_int(fp, "refcoord_scaling", -1, ir1->refcoord_scaling, ir2->refcoord_scaling);
cmp_rvec(fp, "inputrec->posres_com", -1, ir1->posres_com, ir2->posres_com, ftol, abstol);
cmp_rvec(fp, "inputrec->posres_comB", -1, ir1->posres_comB, ir2->posres_comB, ftol, abstol);
diff --git a/src/gromacs/mdtypes/inputrec.h b/src/gromacs/mdtypes/inputrec.h
index 6e4ee727ab..fb0532c8f4 100644
--- a/src/gromacs/mdtypes/inputrec.h
+++ b/src/gromacs/mdtypes/inputrec.h
@@ -400,6 +400,10 @@ struct t_inputrec // NOLINT (clang-analyzer-optin.performance.Padding)
tensor ref_p;
//! Compressibility ((mol nm^3)/kJ)
tensor compress;
+ //! Initial box (nm)
+ tensor h_start;
+ //! Whether deviatoric stress is present in anisotropic C-rescale
+ bool deviatoric_crescale;
//! How to scale absolute reference coordinates
int refcoord_scaling;
//! The COM of the posres atoms
EOF