← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to src/USER-MISC/fix_pimd.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/USER-MISC/fix_pimd.cpp
 
1
/* ----------------------------------------------------------------------
 
2
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
 
3
   http://lammps.sandia.gov, Sandia National Laboratories
 
4
   Steve Plimpton, sjplimp@sandia.gov
 
5
 
 
6
   Copyright (2003) Sandia Corporation.  Under the terms of Contract
 
7
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
 
8
   certain rights in this software.  This software is distributed under
 
9
   the GNU General Public License.
 
10
 
 
11
   See the README file in the top-level LAMMPS directory.
 
12
------------------------------------------------------------------------- */
 
13
 
 
14
/* ----------------------------------------------------------------------
 
15
   Package      FixPIMD
 
16
   Purpose      Quantum Path Integral Algorithm for Quantum Chemistry
 
17
   Copyright    Voth Group @ University of Chicago
 
18
   Authors      Chris Knight & Yuxing Peng (yuxing at uchicago.edu)
 
19
   
 
20
   Updated      Oct-01-2011
 
21
   Version      1.0
 
22
------------------------------------------------------------------------- */
 
23
 
 
24
#include "math.h"
 
25
#include "string.h"
 
26
#include "stdlib.h"
 
27
#include "fix_pimd.h"
 
28
#include "universe.h"
 
29
#include "comm.h"
 
30
#include "force.h"
 
31
#include "atom.h"
 
32
#include "domain.h"
 
33
#include "update.h"
 
34
#include "memory.h"
 
35
#include "error.h"
 
36
 
 
37
using namespace LAMMPS_NS;
 
38
using namespace FixConst;
 
39
 
 
40
enum{PIMD,NMPIMD,CMD};
 
41
 
 
42
/* ---------------------------------------------------------------------- */
 
43
 
 
44
FixPIMD::FixPIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
 
45
{ 
 
46
  method     = PIMD;
 
47
  fmass      = 1.0;
 
48
  nhc_temp   = 298.15;
 
49
  nhc_nchain = 2;
 
50
  sp         = 1.0;
 
51
 
 
52
  for(int i=3; i<narg-1; i+=2)
 
53
  {
 
54
    if(strcmp(arg[i],"method")==0)
 
55
    {
 
56
      if(strcmp(arg[i+1],"pimd")==0) method=PIMD;
 
57
      else if(strcmp(arg[i+1],"nmpimd")==0) method=NMPIMD;
 
58
      else if(strcmp(arg[i+1],"cmd")==0) method=CMD;
 
59
      else error->universe_all(FLERR,"Unkown method parameter for fix pimd");
 
60
    }
 
61
    else if(strcmp(arg[i],"fmass")==0)
 
62
    {
 
63
      fmass = atof(arg[i+1]);
 
64
      if(fmass<0.0 || fmass>1.0) error->universe_all(FLERR,"Invalid fmass value for fix pimd");
 
65
    }
 
66
    else if(strcmp(arg[i],"sp")==0)
 
67
    {
 
68
      sp = atof(arg[i+1]);
 
69
      if(fmass<0.0) error->universe_all(FLERR,"Invalid sp value for fix pimd");
 
70
    }
 
71
    else if(strcmp(arg[i],"temp")==0)
 
72
    {
 
73
      nhc_temp = atof(arg[i+1]);
 
74
      if(nhc_temp<0.0) error->universe_all(FLERR,"Invalid temp value for fix pimd");
 
75
    }
 
76
    else if(strcmp(arg[i],"nhc")==0)
 
77
    {
 
78
      nhc_nchain = atoi(arg[i+1]);
 
79
      if(nhc_nchain<2) error->universe_all(FLERR,"Invalid nhc value for fix pimd");
 
80
    }
 
81
    else error->universe_all(arg[i],i+1,"Unkown keyword for fix pimd");
 
82
  }
 
83
  
 
84
  /* Initiation */
 
85
  
 
86
  max_nsend = 0;
 
87
  tag_send = NULL;
 
88
  buf_send = NULL;
 
89
  
 
90
  max_nlocal = 0;
 
91
  buf_recv = NULL;
 
92
  buf_beads = NULL;
 
93
  
 
94
  size_plan = 0;
 
95
  plan_send = plan_recv = NULL;
 
96
  
 
97
  M_x2xp = M_xp2x = M_f2fp = M_fp2f = NULL;
 
98
  lam = NULL;
 
99
  mode_index = NULL;
 
100
 
 
101
  mass = NULL;
 
102
  
 
103
  array_atom = NULL;
 
104
  nhc_eta = NULL;
 
105
  nhc_eta_dot = NULL;
 
106
  nhc_eta_dotdot = NULL;
 
107
  nhc_eta_mass = NULL;
 
108
  
 
109
  size_peratom_cols = 12 * nhc_nchain + 3;
 
110
  
 
111
  nhc_offset_one_1 = 3 * nhc_nchain;
 
112
  nhc_offset_one_2 = 3 * nhc_nchain +3;
 
113
  nhc_size_one_1 = sizeof(double) * nhc_offset_one_1;
 
114
  nhc_size_one_2 = sizeof(double) * nhc_offset_one_2;
 
115
  
 
116
  restart_peratom = 1;
 
117
  peratom_flag    = 1;
 
118
  peratom_freq    = 1;
 
119
  
 
120
  global_freq = 1;
 
121
  thermo_energy = 1;
 
122
  vector_flag = 1;
 
123
  size_vector = 2;
 
124
  extvector   = 1;
 
125
  comm_forward = 3;
 
126
  
 
127
  atom->add_callback(0); // Call LAMMPS to allocate memory for per-atom array
 
128
  atom->add_callback(1); // Call LAMMPS to re-assign restart-data for per-atom array
 
129
  
 
130
  grow_arrays(atom->nmax);
 
131
  
 
132
  // some initilizations
 
133
  
 
134
  nhc_ready = false;
 
135
}
 
136
 
 
137
/* ---------------------------------------------------------------------- */
 
138
 
 
139
int FixPIMD::setmask()
 
140
{
 
141
  int mask = 0;
 
142
  mask |= POST_FORCE;
 
143
  mask |= INITIAL_INTEGRATE;
 
144
  mask |= FINAL_INTEGRATE;
 
145
  return mask;
 
146
}
 
147
 
 
148
/* ---------------------------------------------------------------------- */
 
149
 
 
150
void FixPIMD::init()
 
151
{
 
152
  if(universe->me==0 && screen) fprintf(screen,"Fix pimd initializing Path-Integral ...\n");
 
153
  
 
154
  // prepare the constants
 
155
  
 
156
  np = universe->nworlds;
 
157
  inverse_np = 1.0 / np;
 
158
  
 
159
  /* The first solution for the force constant, using SI units 
 
160
  
 
161
  const double Boltzmann = 1.3806488E-23;    // SI unit: J/K
 
162
  const double Plank     = 6.6260755E-34;    // SI unit: m^2 kg / s
 
163
  
 
164
  double hbar = Plank / ( 2.0 * M_PI ) * sp;
 
165
  double beta = 1.0 / ( Boltzmann * input.nh_temp);
 
166
  
 
167
  // - P / ( beta^2 * hbar^2)   SI unit: s^-2
 
168
  double _fbond = -1.0 / (beta*beta*hbar*hbar) * input.nbeads;
 
169
  
 
170
  // convert the units: s^-2 -> (kcal/mol) / (g/mol) / (A^2)
 
171
  fbond = _fbond * 4.184E+26;
 
172
  
 
173
  */
 
174
  
 
175
  /* The current solution, using LAMMPS internal real units */
 
176
  
 
177
  const double Boltzmann = force->boltz; 
 
178
  const double Plank     = force->hplanck; 
 
179
  
 
180
  double hbar   = Plank / ( 2.0 * M_PI );
 
181
  double beta   = 1.0 / (Boltzmann * nhc_temp);  
 
182
  double _fbond = 1.0 * np / (beta*beta*hbar*hbar) ;
 
183
  
 
184
  omega_np = sqrt(np) / (hbar * beta) * sqrt(force->mvv2e);
 
185
  fbond = - _fbond * force->mvv2e;
 
186
  
 
187
  if(universe->me==0)
 
188
    printf("Fix pimd -P/(beta^2 * hbar^2) = %20.7lE (kcal/mol/A^2)\n\n", fbond);
 
189
  
 
190
  dtv = update->dt;
 
191
  dtf = 0.5 * update->dt * force->ftm2v;
 
192
  
 
193
  comm_init();
 
194
  
 
195
  mass = new double [atom->ntypes+1];
 
196
  
 
197
  if(method==CMD || method==NMPIMD) nmpimd_init();
 
198
  else for(int i=1; i<=atom->ntypes; i++) mass[i] = atom->mass[i] / np * fmass;
 
199
  
 
200
  if(!nhc_ready) nhc_init();
 
201
}
 
202
 
 
203
/* ---------------------------------------------------------------------- */
 
204
 
 
205
void FixPIMD::setup(int vflag)
 
206
{
 
207
  if(universe->me==0 && screen) fprintf(screen,"Setting up Path-Integral ...\n");
 
208
  
 
209
  post_force(vflag);
 
210
}
 
211
 
 
212
/* ---------------------------------------------------------------------- */
 
213
 
 
214
void FixPIMD::initial_integrate(int vflag)
 
215
{
 
216
  nhc_update_v();
 
217
  nhc_update_x();
 
218
}
 
219
 
 
220
/* ---------------------------------------------------------------------- */
 
221
 
 
222
void FixPIMD::final_integrate()
 
223
{
 
224
  nhc_update_v();
 
225
}
 
226
 
 
227
/* ---------------------------------------------------------------------- */
 
228
 
 
229
void FixPIMD::post_force(int flag)
 
230
{
 
231
  for(int i=0; i<atom->nlocal; i++) for(int j=0; j<3; j++) atom->f[i][j] /= np;
 
232
  
 
233
  comm_exec(atom->x);
 
234
  spring_force();
 
235
  
 
236
  if(method==CMD || method==NMPIMD)
 
237
  {
 
238
    /* forward comm for the force on ghost atoms */
 
239
    
 
240
    nmpimd_fill(atom->f);    
 
241
  
 
242
    /* inter-partition comm */
 
243
    
 
244
    comm_exec(atom->f);
 
245
    
 
246
    /* normal-mode transform */
 
247
    
 
248
    nmpimd_transform(buf_beads, atom->f, M_f2fp[universe->iworld]);
 
249
  }
 
250
}
 
251
 
 
252
/* ----------------------------------------------------------------------
 
253
   Nose-Hoover Chains
 
254
------------------------------------------------------------------------- */
 
255
 
 
256
void FixPIMD::nhc_init()
 
257
{
 
258
  double tau = 1.0 / omega_np;
 
259
  double KT  = force->boltz * nhc_temp;
 
260
 
 
261
  double mass0 = KT * tau * tau;
 
262
  int max = 3 * atom->nlocal;
 
263
    
 
264
  for(int i=0; i<max; i++) 
 
265
  {
 
266
    for(int ichain=0; ichain<nhc_nchain; ichain++) 
 
267
    {
 
268
      nhc_eta[i][ichain]        = 0.0;
 
269
      nhc_eta_dot[i][ichain]    = 0.0;
 
270
      nhc_eta_dot[i][ichain]    = 0.0;
 
271
      nhc_eta_dotdot[i][ichain] = 0.0;
 
272
      nhc_eta_mass[i][ichain]   = mass0;
 
273
      nhc_eta_mass[i][ichain] *= fmass; 
 
274
    }
 
275
      
 
276
    nhc_eta_dot[i][nhc_nchain]    = 0.0;
 
277
    
 
278
    for(int ichain=1; ichain<nhc_nchain; ichain++)
 
279
      nhc_eta_dotdot[i][ichain] = (nhc_eta_mass[i][ichain-1] * nhc_eta_dot[i][ichain-1] 
 
280
        * nhc_eta_dot[i][ichain-1] * force->mvv2e - KT) / nhc_eta_mass[i][ichain];
 
281
  }
 
282
 
 
283
  // Zero NH acceleration for CMD
 
284
  
 
285
  if(method==CMD && universe->iworld==0) for(int i=0; i<max; i++) 
 
286
    for(int ichain=0; ichain<nhc_nchain; ichain++) nhc_eta_dotdot[i][ichain] = 0.0;
 
287
  
 
288
  nhc_ready = true;
 
289
}
 
290
 
 
291
/* ---------------------------------------------------------------------- */
 
292
 
 
293
void FixPIMD::nhc_update_x()
 
294
{ 
 
295
  int n = atom->nlocal;
 
296
  double **x = atom->x;
 
297
  double **v = atom->v;
 
298
 
 
299
  if(method==CMD || method==NMPIMD)
 
300
  {
 
301
    nmpimd_fill(atom->v);     
 
302
    comm_exec(atom->v);
 
303
    
 
304
    /* borrow the space of atom->f to store v in cartisian */
 
305
    
 
306
    v = atom->f;      
 
307
    nmpimd_transform(buf_beads, v, M_xp2x[universe->iworld]);
 
308
  }
 
309
 
 
310
  for(int i=0; i<n; i++)
 
311
  {
 
312
    x[i][0] += dtv * v[i][0];
 
313
    x[i][1] += dtv * v[i][1];
 
314
    x[i][2] += dtv * v[i][2];
 
315
  }
 
316
}
 
317
 
 
318
/* ---------------------------------------------------------------------- */
 
319
 
 
320
void FixPIMD::nhc_update_v()
 
321
{
 
322
  int n = atom->nlocal;
 
323
  int *type = atom->type;
 
324
  double **v = atom->v;
 
325
  double **f = atom->f;
 
326
  
 
327
  for(int i=0; i<n; i++)
 
328
  {
 
329
    double dtfm = dtf / mass[type[i]];
 
330
    v[i][0] += dtfm * f[i][0];
 
331
    v[i][1] += dtfm * f[i][1];
 
332
    v[i][2] += dtfm * f[i][2];
 
333
  }
 
334
  
 
335
  t_sys = 0.0;
 
336
  if(method==CMD && universe->iworld==0) return;
 
337
  
 
338
  double expfac;
 
339
  int nmax = 3 * atom->nlocal;
 
340
  double KT = force->boltz * nhc_temp;
 
341
  double kecurrent, t_current;
 
342
  
 
343
  double dthalf = 0.5   * update->dt;
 
344
  double dt4    = 0.25  * update->dt;
 
345
  double dt8    = 0.125 * update->dt;
 
346
  
 
347
  for(int i=0; i<nmax; i++) 
 
348
  {
 
349
    int iatm = i/3;
 
350
    int idim = i%3;
 
351
      
 
352
    double *vv = v[iatm];      
 
353
      
 
354
    kecurrent = mass[type[iatm]] * vv[idim]* vv[idim] * force->mvv2e;
 
355
    t_current = kecurrent / force->boltz;
 
356
      
 
357
    double *eta = nhc_eta[i];
 
358
    double *eta_dot = nhc_eta_dot[i];
 
359
    double *eta_dotdot = nhc_eta_dotdot[i];
 
360
      
 
361
    eta_dotdot[0] = (kecurrent - KT) / nhc_eta_mass[i][0];
 
362
 
 
363
    for(int ichain=nhc_nchain-1; ichain>0; ichain--) 
 
364
    {
 
365
      expfac = exp(-dt8 * eta_dot[ichain+1]);
 
366
      eta_dot[ichain] *= expfac;
 
367
      eta_dot[ichain] += eta_dotdot[ichain] * dt4;
 
368
      eta_dot[ichain] *= expfac;
 
369
    }
 
370
 
 
371
    expfac = exp(-dt8 * eta_dot[1]);
 
372
    eta_dot[0] *= expfac;
 
373
    eta_dot[0] += eta_dotdot[0] * dt4;
 
374
    eta_dot[0] *= expfac;
 
375
 
 
376
    // Update particle velocities half-step
 
377
 
 
378
    double factor_eta = exp(-dthalf * eta_dot[0]);
 
379
    vv[idim] *= factor_eta;
 
380
 
 
381
    t_current *= (factor_eta * factor_eta);
 
382
    kecurrent = force->boltz * t_current;
 
383
    eta_dotdot[0] = (kecurrent - KT) / nhc_eta_mass[i][0];
 
384
 
 
385
    for(int ichain=0; ichain<nhc_nchain; ichain++) 
 
386
      eta[ichain] += dthalf * eta_dot[ichain];
 
387
 
 
388
    eta_dot[0] *= expfac;
 
389
    eta_dot[0] += eta_dotdot[0] * dt4;
 
390
    eta_dot[0] *= expfac;
 
391
 
 
392
    for(int ichain=1; ichain<nhc_nchain; ichain++) 
 
393
    {
 
394
      expfac = exp(-dt8 * eta_dot[ichain+1]);
 
395
      eta_dot[ichain] *= expfac;
 
396
      eta_dotdot[ichain] = (nhc_eta_mass[i][ichain-1] * eta_dot[ichain-1] * eta_dot[ichain-1]
 
397
                           - KT) / nhc_eta_mass[i][ichain];
 
398
      eta_dot[ichain] += eta_dotdot[ichain] * dt4;
 
399
      eta_dot[ichain] *= expfac;
 
400
    }
 
401
    
 
402
    t_sys += t_current;    
 
403
  }
 
404
  
 
405
  t_sys /= nmax;
 
406
}
 
407
 
 
408
/* ----------------------------------------------------------------------
 
409
   Normal Mode PIMD
 
410
------------------------------------------------------------------------- */
 
411
 
 
412
void FixPIMD::nmpimd_init()
 
413
{
 
414
  memory->create(M_x2xp, np, np, "fix_feynman:M_x2xp");
 
415
  memory->create(M_xp2x, np, np, "fix_feynman:M_xp2x");
 
416
  memory->create(M_f2fp, np, np, "fix_feynman:M_f2fp");
 
417
  memory->create(M_fp2f, np, np, "fix_feynman:M_fp2f");
 
418
  
 
419
  lam = (double*) memory->smalloc(sizeof(double)*np, "FixPIMD::lam");
 
420
 
 
421
  // Set up  eigenvalues
 
422
  
 
423
  lam[0] = 0.0;
 
424
  if(np%2==0) lam[np-1] = 4.0 * np;
 
425
  
 
426
  for(int i=2; i<=np/2; i++)
 
427
  {
 
428
    lam[2*i-3] = lam[2*i-2] = 2.0 * np * (1.0 - 1.0 *cos(2.0*M_PI*(i-1)/np));
 
429
  }
 
430
  
 
431
  // Set up eigenvectors for non-degenerated modes
 
432
  
 
433
  for(int i=0; i<np; i++)
 
434
  {
 
435
    M_x2xp[0][i] = 1.0 / np;
 
436
    if(np%2==0) M_x2xp[np-1][i] = 1.0 / np * pow(-1.0, i);
 
437
  }
 
438
  
 
439
  // Set up eigenvectors for degenerated modes
 
440
  
 
441
  for(int i=0; i<(np-1)/2; i++) for(int j=0; j<np; j++)
 
442
  {
 
443
    M_x2xp[2*i+1][j] =   sqrt(2.0) * cos ( 2.0 * M_PI * (i+1) * j / np) / np;
 
444
    M_x2xp[2*i+2][j] = - sqrt(2.0) * sin ( 2.0 * M_PI * (i+1) * j / np) / np;
 
445
  }
 
446
  
 
447
  // Set up Ut
 
448
  
 
449
  for(int i=0; i<np; i++) 
 
450
    for(int j=0; j<np; j++)
 
451
    {
 
452
      M_xp2x[i][j] = M_x2xp[j][i] * np; 
 
453
      M_f2fp[i][j] = M_x2xp[i][j] * np;
 
454
      M_fp2f[i][j] = M_xp2x[i][j];
 
455
    }
 
456
  
 
457
  // Set up masses
 
458
  
 
459
  int iworld = universe->iworld;
 
460
  
 
461
  for(int i=1; i<=atom->ntypes; i++)
 
462
  {
 
463
    mass[i] = atom->mass[i];
 
464
 
 
465
    if(iworld) 
 
466
    {
 
467
      mass[i] *= lam[iworld];
 
468
      mass[i] *= fmass;
 
469
    }
 
470
  }
 
471
}
 
472
 
 
473
/* ---------------------------------------------------------------------- */
 
474
 
 
475
void FixPIMD::nmpimd_fill(double **ptr)
 
476
{
 
477
  comm_ptr = ptr;
 
478
  comm->forward_comm_fix(this);
 
479
}
 
480
 
 
481
/* ---------------------------------------------------------------------- */
 
482
 
 
483
void FixPIMD::nmpimd_transform(double** src, double** des, double *vector)
 
484
{
 
485
  int n = atom->nlocal;
 
486
  int m = 0;
 
487
  
 
488
  for(int i=0; i<n; i++) for(int d=0; d<3; d++)
 
489
  {
 
490
    des[i][d] = 0.0;   
 
491
    for(int j=0; j<np; j++) { des[i][d] += (src[j][m] * vector[j]); }
 
492
    m++;
 
493
  }
 
494
}
 
495
 
 
496
/* ---------------------------------------------------------------------- */
 
497
 
 
498
void FixPIMD::spring_force()
 
499
{
 
500
  spring_energy = 0.0;
 
501
 
 
502
  double **x = atom->x;
 
503
  double **f = atom->f;
 
504
  double* _mass = atom->mass;
 
505
  int* type = atom->type;
 
506
  int nlocal = atom->nlocal;
 
507
 
 
508
  double* xlast = buf_beads[x_last];  
 
509
  double* xnext = buf_beads[x_next];
 
510
  
 
511
  for(int i=0; i<nlocal; i++)
 
512
  {
 
513
    double delx1 = xlast[0] - x[i][0];
 
514
    double dely1 = xlast[1] - x[i][1];
 
515
    double delz1 = xlast[2] - x[i][2];
 
516
    xlast += 3;
 
517
    domain->minimum_image(delx1, dely1, delz1);
 
518
 
 
519
    double delx2 = xnext[0] - x[i][0];
 
520
    double dely2 = xnext[1] - x[i][1];
 
521
    double delz2 = xnext[2] - x[i][2];
 
522
    xnext += 3;
 
523
    domain->minimum_image(delx2, dely2, delz2);
 
524
    
 
525
    double ff = fbond * _mass[type[i]];
 
526
    
 
527
    double dx = delx1+delx2;
 
528
    double dy = dely1+dely2;
 
529
    double dz = delz1+delz2;
 
530
    
 
531
    f[i][0] -= (dx) * ff;
 
532
    f[i][1] -= (dy) * ff;
 
533
    f[i][2] -= (dz) * ff;
 
534
    
 
535
    spring_energy += (dx*dx+dy*dy+dz*dz);
 
536
  }
 
537
}
 
538
 
 
539
/* ----------------------------------------------------------------------
 
540
   Comm operations
 
541
------------------------------------------------------------------------- */
 
542
 
 
543
void FixPIMD::comm_init()
 
544
{
 
545
  if(size_plan)
 
546
  {
 
547
    delete [] plan_send;
 
548
    delete [] plan_recv;
 
549
  }
 
550
  
 
551
  if(method == PIMD) 
 
552
  {
 
553
    size_plan = 2;
 
554
    plan_send = new int [2];
 
555
    plan_recv = new int [2];
 
556
    mode_index = new int [2];
 
557
    
 
558
    int rank_last = universe->me - comm->nprocs;
 
559
    int rank_next = universe->me + comm->nprocs;    
 
560
    if(rank_last<0) rank_last += universe->nprocs;
 
561
    if(rank_next>=universe->nprocs) rank_next -= universe->nprocs;
 
562
    
 
563
    plan_send[0] = rank_next; plan_send[1] = rank_last;
 
564
    plan_recv[0] = rank_last; plan_recv[1] = rank_next;
 
565
    
 
566
    mode_index[0] = 0; mode_index[1] = 1;
 
567
    x_last = 1; x_next = 0;
 
568
  }
 
569
  else 
 
570
  {
 
571
    size_plan = np - 1;
 
572
    plan_send = new int [size_plan];
 
573
    plan_recv = new int [size_plan];
 
574
    mode_index = new int [size_plan];
 
575
 
 
576
    for(int i=0; i<size_plan; i++)
 
577
    {
 
578
      plan_send[i] = universe->me + comm->nprocs * (i+1);
 
579
      if(plan_send[i]>=universe->nprocs) plan_send[i] -= universe->nprocs;
 
580
      
 
581
      plan_recv[i] = universe->me - comm->nprocs * (i+1);
 
582
      if(plan_recv[i]<0) plan_recv[i] += universe->nprocs;
 
583
 
 
584
      mode_index[i]=(universe->iworld+i+1)%(universe->nworlds);
 
585
    }
 
586
    
 
587
    x_next = (universe->iworld+1+universe->nworlds)%(universe->nworlds);
 
588
    x_last = (universe->iworld-1+universe->nworlds)%(universe->nworlds);
 
589
  }
 
590
  
 
591
  if(buf_beads)
 
592
  {
 
593
    for(int i=0; i<np; i++) if(buf_beads[i]) delete [] buf_beads[i];
 
594
    delete [] buf_beads;
 
595
  }
 
596
  
 
597
  buf_beads = new double* [np];
 
598
  for(int i=0; i<np; i++) buf_beads[i] = NULL;
 
599
}
 
600
 
 
601
/* ---------------------------------------------------------------------- */
 
602
 
 
603
void FixPIMD::comm_exec(double **ptr)
 
604
{
 
605
  int nlocal = atom->nlocal;
 
606
  
 
607
  if(nlocal > max_nlocal)
 
608
  {
 
609
    max_nlocal = nlocal+200;
 
610
    int size = sizeof(double) * max_nlocal * 3;
 
611
    buf_recv = (double*) memory->srealloc(buf_recv, size, "FixPIMD:x_recv");
 
612
    
 
613
    for(int i=0; i<np; i++) 
 
614
      buf_beads[i] = (double*) memory->srealloc(buf_beads[i], size, "FixPIMD:x_beads[i]");
 
615
  }
 
616
  
 
617
  // copy local positions
 
618
  
 
619
  memcpy(buf_beads[universe->iworld], &(ptr[0][0]), sizeof(double)*nlocal*3);
 
620
  
 
621
  // go over comm plans
 
622
  
 
623
  for(int iplan = 0; iplan<size_plan; iplan++)
 
624
  {  
 
625
    // sendrecv nlocal
 
626
   
 
627
    int nsend;
 
628
    
 
629
    MPI_Sendrecv( &(nlocal), 1, MPI_INT, plan_send[iplan], 0, 
 
630
                  &(nsend),  1, MPI_INT, plan_recv[iplan], 0, universe->uworld, MPI_STATUS_IGNORE);
 
631
  
 
632
    // allocate arrays
 
633
    
 
634
    if(nsend > max_nsend)
 
635
    {
 
636
      max_nsend = nsend+200;
 
637
      tag_send = (int*) memory->srealloc(tag_send, sizeof(int)*max_nsend, "FixPIMD:tag_send");
 
638
      buf_send = (double*) memory->srealloc(buf_send, sizeof(double)*max_nsend*3, "FixPIMD:x_send");
 
639
    }
 
640
  
 
641
    // send tags 
 
642
  
 
643
    MPI_Sendrecv( atom->tag, nlocal, MPI_INT, plan_send[iplan], 0, 
 
644
                  tag_send,  nsend,  MPI_INT, plan_recv[iplan], 0, universe->uworld, MPI_STATUS_IGNORE);
 
645
 
 
646
    // wrap positions
 
647
  
 
648
    double *wrap_ptr = buf_send;
 
649
    int ncpy = sizeof(double)*3;
 
650
    
 
651
    for(int i=0; i<nsend; i++)
 
652
    {
 
653
      int index = atom->map(tag_send[i]);
 
654
 
 
655
      if(index<0)
 
656
      {
 
657
        char error_line[256];
 
658
      
 
659
        sprintf(error_line, "Atom %d is missing at world [%d] rank [%d] required by  rank [%d].\n",
 
660
          tag_send[i], universe->iworld, comm->me, plan_recv[iplan]);
 
661
          
 
662
        error->universe_one(FLERR,error_line);
 
663
      }
 
664
    
 
665
      memcpy(wrap_ptr, ptr[index], ncpy);
 
666
      wrap_ptr += 3;
 
667
    }
 
668
 
 
669
    // sendrecv x 
 
670
  
 
671
    MPI_Sendrecv( buf_send, nsend*3,  MPI_DOUBLE, plan_recv[iplan], 0, 
 
672
                  buf_recv, nlocal*3, MPI_DOUBLE, plan_send[iplan], 0, universe->uworld, MPI_STATUS_IGNORE);
 
673
  
 
674
    // copy x
 
675
    
 
676
    memcpy(buf_beads[mode_index[iplan]], buf_recv, sizeof(double)*nlocal*3);
 
677
  }
 
678
}
 
679
 
 
680
/* ---------------------------------------------------------------------- */
 
681
 
 
682
int FixPIMD::pack_comm(int n, int *list, double *buf,
 
683
                             int pbc_flag, int *pbc)
 
684
{
 
685
  int i,j,m;
 
686
 
 
687
  m = 0;
 
688
  
 
689
  for (i = 0; i < n; i++) {
 
690
    j = list[i];
 
691
    buf[m++] = comm_ptr[j][0];
 
692
    buf[m++] = comm_ptr[j][1];
 
693
    buf[m++] = comm_ptr[j][2];
 
694
  }
 
695
  
 
696
  return 3;
 
697
}
 
698
 
 
699
/* ---------------------------------------------------------------------- */
 
700
 
 
701
void FixPIMD::unpack_comm(int n, int first, double *buf)
 
702
{
 
703
  int i,m,last;
 
704
 
 
705
  m = 0;
 
706
  last = first + n;
 
707
  for (i = first; i < last; i++) {
 
708
    comm_ptr[i][0] = buf[m++];
 
709
    comm_ptr[i][1] = buf[m++];
 
710
    comm_ptr[i][2] = buf[m++];
 
711
  }
 
712
}
 
713
 
 
714
/* ----------------------------------------------------------------------
 
715
   Memory operations
 
716
------------------------------------------------------------------------- */
 
717
 
 
718
double FixPIMD::memory_usage()
 
719
{
 
720
  double bytes = 0;
 
721
  bytes = atom->nmax * size_peratom_cols * sizeof(double);
 
722
  return bytes;
 
723
}
 
724
 
 
725
/* ---------------------------------------------------------------------- */
 
726
 
 
727
void FixPIMD::grow_arrays(int nmax)
 
728
{
 
729
  if (nmax==0) return;
 
730
  int count = nmax*3;
 
731
 
 
732
  memory->grow(array_atom, nmax, size_peratom_cols, "FixPIMD::array_atom");
 
733
  memory->grow(nhc_eta,        count, nhc_nchain,   "FixPIMD::nh_eta");     
 
734
  memory->grow(nhc_eta_dot,    count, nhc_nchain+1, "FixPIMD::nh_eta_dot");   
 
735
  memory->grow(nhc_eta_dotdot, count, nhc_nchain,   "FixPIMD::nh_eta_dotdot");
 
736
  memory->grow(nhc_eta_mass,   count, nhc_nchain,   "FixPIMD::nh_eta_mass"); 
 
737
}
 
738
 
 
739
/* ---------------------------------------------------------------------- */
 
740
 
 
741
void FixPIMD::copy_arrays(int i, int j, int delflag)
 
742
{
 
743
  int i_pos = i*3;
 
744
  int j_pos = j*3;
 
745
  
 
746
  memcpy(nhc_eta       [j_pos], nhc_eta       [i_pos], nhc_size_one_1);
 
747
  memcpy(nhc_eta_dot   [j_pos], nhc_eta_dot   [i_pos], nhc_size_one_2);
 
748
  memcpy(nhc_eta_dotdot[j_pos], nhc_eta_dotdot[i_pos], nhc_size_one_1);
 
749
  memcpy(nhc_eta_mass  [j_pos], nhc_eta_mass  [i_pos], nhc_size_one_1);
 
750
}
 
751
 
 
752
/* ---------------------------------------------------------------------- */
 
753
 
 
754
int FixPIMD::pack_exchange(int i, double *buf)
 
755
{
 
756
  int offset=0;  
 
757
  int pos = i * 3;
 
758
  
 
759
  memcpy(buf+offset, nhc_eta[pos],        nhc_size_one_1); offset += nhc_offset_one_1;
 
760
  memcpy(buf+offset, nhc_eta_dot[pos],    nhc_size_one_2); offset += nhc_offset_one_2;
 
761
  memcpy(buf+offset, nhc_eta_dotdot[pos], nhc_size_one_1); offset += nhc_offset_one_1;
 
762
  memcpy(buf+offset, nhc_eta_mass[pos],   nhc_size_one_1); offset += nhc_offset_one_1;
 
763
  
 
764
  return size_peratom_cols;
 
765
}
 
766
 
 
767
/* ---------------------------------------------------------------------- */
 
768
 
 
769
int FixPIMD::unpack_exchange(int nlocal, double *buf)
 
770
{
 
771
  int offset=0;
 
772
  int pos = nlocal*3;
 
773
  
 
774
  memcpy(nhc_eta[pos],        buf+offset, nhc_size_one_1); offset += nhc_offset_one_1;
 
775
  memcpy(nhc_eta_dot[pos],    buf+offset, nhc_size_one_2); offset += nhc_offset_one_2;
 
776
  memcpy(nhc_eta_dotdot[pos], buf+offset, nhc_size_one_1); offset += nhc_offset_one_1;
 
777
  memcpy(nhc_eta_mass[pos],   buf+offset, nhc_size_one_1); offset += nhc_offset_one_1;
 
778
  
 
779
  return size_peratom_cols;
 
780
}
 
781
 
 
782
/* ---------------------------------------------------------------------- */
 
783
 
 
784
int FixPIMD::pack_restart(int i, double *buf)
 
785
{
 
786
  int offset=0;
 
787
  int pos = i * 3;
 
788
  buf[offset++] = size_peratom_cols+1;
 
789
 
 
790
  memcpy(buf+offset, nhc_eta[pos],        nhc_size_one_1); offset += nhc_offset_one_1;
 
791
  memcpy(buf+offset, nhc_eta_dot[pos],    nhc_size_one_2); offset += nhc_offset_one_2;
 
792
  memcpy(buf+offset, nhc_eta_dotdot[pos], nhc_size_one_1); offset += nhc_offset_one_1;
 
793
  memcpy(buf+offset, nhc_eta_mass[pos],   nhc_size_one_1); offset += nhc_offset_one_1;
 
794
  
 
795
  return size_peratom_cols+1;  
 
796
}
 
797
 
 
798
/* ---------------------------------------------------------------------- */
 
799
 
 
800
void FixPIMD::unpack_restart(int nlocal, int nth)
 
801
{
 
802
  double **extra = atom->extra;
 
803
 
 
804
  // skip to Nth set of extra values
 
805
  
 
806
  int m = 0;
 
807
  for (int i=0; i<nth; i++) m += static_cast<int> (extra[nlocal][m]);
 
808
  m++;
 
809
  
 
810
  int pos = nlocal * 3;
 
811
    
 
812
  memcpy(nhc_eta[pos],        extra[nlocal]+m, nhc_size_one_1); m += nhc_offset_one_1;
 
813
  memcpy(nhc_eta_dot[pos],    extra[nlocal]+m, nhc_size_one_2); m += nhc_offset_one_2;
 
814
  memcpy(nhc_eta_dotdot[pos], extra[nlocal]+m, nhc_size_one_1); m += nhc_offset_one_1;
 
815
  memcpy(nhc_eta_mass[pos],   extra[nlocal]+m, nhc_size_one_1); m += nhc_offset_one_1;
 
816
  
 
817
  nhc_ready = true;
 
818
}
 
819
 
 
820
/* ---------------------------------------------------------------------- */
 
821
 
 
822
int FixPIMD::maxsize_restart()
 
823
{
 
824
  return size_peratom_cols+1;
 
825
}
 
826
 
 
827
/* ---------------------------------------------------------------------- */
 
828
 
 
829
int FixPIMD::size_restart(int nlocal)
 
830
{
 
831
  return size_peratom_cols+1;
 
832
}
 
833
 
 
834
/* ---------------------------------------------------------------------- */
 
835
 
 
836
double FixPIMD::compute_vector(int n)
 
837
{ 
 
838
  if(n==0) { return spring_energy; }
 
839
  if(n==1) { return t_sys; }
 
840
  return 0.0;
 
841
}