← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/USER-MISC/fix_srp.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_srp.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
   Contributing authors: Timothy Sirk (ARL), Pieter in't Veld (BASF)
 
16
------------------------------------------------------------------------- */
 
17
 
 
18
#include "string.h"
 
19
#include "stdlib.h"
 
20
#include "fix_srp.h"
 
21
#include "atom.h"
 
22
#include "force.h"
 
23
#include "domain.h"
 
24
#include "comm.h"
 
25
#include "memory.h"
 
26
#include "error.h"
 
27
#include "neighbor.h"
 
28
#include "atom_vec.h"
 
29
#include "modify.h"
 
30
 
 
31
using namespace LAMMPS_NS;
 
32
using namespace FixConst;
 
33
 
 
34
/* ---------------------------------------------------------------------- */
 
35
 
 
36
FixSRP::FixSRP(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
 
37
{
 
38
  // settings
 
39
  nevery=1;
 
40
  peratom_freq = 1;
 
41
  time_integrate = 0;
 
42
  create_attribute = 0;
 
43
  comm_border = 2;
 
44
 
 
45
  // restart settings
 
46
  restart_global = 1;
 
47
  restart_peratom = 1;
 
48
  restart_pbc = 1;
 
49
 
 
50
  // per-atom array width 2
 
51
  peratom_flag = 1;
 
52
  size_peratom_cols = 2;                                                  
 
53
 
 
54
  // initial allocation of atom-based array                      
 
55
  // register with Atom class
 
56
  array = NULL;
 
57
  grow_arrays(atom->nmax);                                                
 
58
 
 
59
  // extends pack_exchange() 
 
60
  atom->add_callback(0);
 
61
  atom->add_callback(1); // restart
 
62
  atom->add_callback(2); 
 
63
 
 
64
  // zero 
 
65
  for (int i = 0; i < atom->nmax; i++) 
 
66
    for (int m = 0; m < 3; m++) 
 
67
      array[i][m] = 0.0; 
 
68
 
 
69
  // assume setup of fix is needed to insert particles
 
70
  // is true if reading from datafile
 
71
  // since a datafile cannot contain bond particles
 
72
  // might be true if reading from restart
 
73
  // a restart file written during the run has bond particles as per atom data
 
74
  // a restart file written after the run does not have bond particles
 
75
  setup = 1; 
 
76
}
 
77
 
 
78
/* ---------------------------------------------------------------------- */
 
79
 
 
80
FixSRP::~FixSRP()
 
81
{
 
82
  // unregister callbacks to this fix from Atom class
 
83
  atom->delete_callback(id,0);
 
84
  atom->delete_callback(id,1);
 
85
  atom->delete_callback(id,2);
 
86
  memory->destroy(array);
 
87
}
 
88
 
 
89
/* ---------------------------------------------------------------------- */
 
90
 
 
91
int FixSRP::setmask()
 
92
{
 
93
  int mask = 0;
 
94
  mask |= PRE_FORCE;
 
95
  mask |= PRE_EXCHANGE;
 
96
  mask |= POST_RUN;
 
97
 
 
98
  return mask;
 
99
}
 
100
 
 
101
/* ---------------------------------------------------------------------- */
 
102
 
 
103
void FixSRP::init()
 
104
{
 
105
  if (force->pair_match("hybrid",1) == NULL)
 
106
    error->all(FLERR,"Cannot use pair srp without pair_style hybrid");
 
107
 
 
108
  // the highest numbered atom type should be reserved for bond particles (bptype)
 
109
  // set bptype, unless it will be read from restart 
 
110
  if(!restart_reset) bptype = atom->ntypes;
 
111
 
 
112
  // check if bptype is already in use
 
113
   for(int i=0; i < atom->nlocal; i++)
 
114
      if(atom->type[i] == bptype){
 
115
        // error if bptype is already in use
 
116
        // unless starting from a rst file
 
117
        // since rst can contain the bond particles as per atom data
 
118
        if(!restart_reset)
 
119
          error->all(FLERR,"Fix SRP requires an extra atom type");
 
120
        else
 
121
          setup = 0; 
 
122
      }
 
123
 
 
124
  // setup neigh exclusions for diff atom types
 
125
  // bond particles do not interact with other types
 
126
  // type bptype only interacts with itself
 
127
  char* arg1[4];
 
128
  arg1[0] = (char *) "exclude"; 
 
129
  arg1[1] = (char *) "type"; 
 
130
  char c0[20];
 
131
  char c1[20];
 
132
 
 
133
  for(int z = 1; z < bptype; z++)
 
134
  {
 
135
   sprintf(c0, "%d", z);
 
136
   arg1[2] = c0; 
 
137
   sprintf(c1, "%d", bptype);
 
138
   arg1[3] = c1; 
 
139
   neighbor->modify_params(4, arg1);
 
140
  }
 
141
}
 
142
 
 
143
/* ----------------------------------------------------------------------
 
144
   insert bond particles
 
145
------------------------------------------------------------------------- */
 
146
 
 
147
void FixSRP::setup_pre_force(int zz)
 
148
{
 
149
  if(!setup) return;
 
150
 
 
151
  double rsqold = 0.0;
 
152
  double delx, dely, delz, rmax, rsq, rsqmax;
 
153
  double **x = atom->x;
 
154
  bigint nall = atom->nlocal + atom->nghost;
 
155
  double xold[nall][3];
 
156
 
 
157
  // make a copy of all coordinates and tags
 
158
  // need this because create_atom overwrites ghost atoms
 
159
  for(int i = 0; i < nall; i++){
 
160
    xold[i][0] = x[i][0];
 
161
    xold[i][1] = x[i][1];
 
162
    xold[i][2] = x[i][2];
 
163
  }
 
164
 
 
165
  tagint *tag = atom->tag;
 
166
  tagint tagold[nall];
 
167
 
 
168
  for(int i = 0; i < nall; i++){
 
169
    tagold[i]=tag[i];
 
170
  }
 
171
 
 
172
  int nlocal = atom->nlocal;
 
173
  int **bondlist = neighbor->bondlist;
 
174
  int nbondlist = neighbor->nbondlist;
 
175
  int nadd = 0;
 
176
  int i,j;
 
177
 
 
178
  for (int n = 0; n < nbondlist; n++) {
 
179
 
 
180
   // consider only the user defined bond type
 
181
   if(bondlist[n][2] != btype) continue;
 
182
 
 
183
   i = bondlist[n][0];
 
184
   j = bondlist[n][1];
 
185
 
 
186
   // position of bond i
 
187
   xone[0] = (xold[i][0] + xold[j][0])*0.5;
 
188
   xone[1] = (xold[i][1] + xold[j][1])*0.5;
 
189
   xone[2] = (xold[i][2] + xold[j][2])*0.5;
 
190
 
 
191
   // record longest bond 
 
192
   // this used to set ghost cutoff
 
193
    delx = xold[j][0] - xold[i][0];
 
194
    dely = xold[j][1] - xold[i][1];
 
195
    delz = xold[j][2] - xold[i][2];
 
196
    rsq = delx*delx + dely*dely + delz*delz;
 
197
    if(rsq > rsqold) rsqold = rsq; 
 
198
 
 
199
   // make one particle for each bond
 
200
   // i is local
 
201
   // if newton bond, always make particle
 
202
   // if j is local, always make particle
 
203
   // if j is ghost, decide from tag  
 
204
 
 
205
   if( force->newton_bond || j < nlocal || tagold[i] > tagold[j] ){
 
206
        atom->natoms += 1;
 
207
        atom->avec->create_atom(bptype,xone);
 
208
        // pack tag i/j into buffer for comm
 
209
        array[atom->nlocal-1][0] =  (double)tagold[i];
 
210
        array[atom->nlocal-1][1] =  (double)tagold[j];
 
211
        nadd++;
 
212
    }
 
213
  } 
 
214
 
 
215
  // new total # of atoms
 
216
  bigint nblocal = atom->nlocal;
 
217
  MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
 
218
 
 
219
  // assign tags for new atoms, update map
 
220
  if (atom->tag_enable) atom->tag_extend();
 
221
 
 
222
  if (atom->map_style) {
 
223
    atom->nghost = 0;
 
224
    atom->map_init();
 
225
    atom->map_set();
 
226
  }
 
227
 
 
228
  char str[128];
 
229
  int Nadd = 0;
 
230
  MPI_Allreduce(&nadd,&Nadd,1,MPI_INT,MPI_SUM,world);
 
231
  if(comm->me == 0){
 
232
    sprintf(str, "Inserted %d bond particles.", Nadd);
 
233
    error->message(FLERR,str);
 
234
  }
 
235
 
 
236
  // check ghost comm distances
 
237
  // warn and change if shorter from estimate
 
238
  // ghost atoms must be present for bonds on edge of neighbor cutoff
 
239
  // extend cutghost slightly more than half of the longest bond
 
240
  MPI_Allreduce(&rsqold,&rsqmax,1,MPI_DOUBLE,MPI_MAX,world);
 
241
  rmax = sqrt(rsqmax); 
 
242
  double cutneighmax_srp = neighbor->cutneighmax + 0.51*rmax;
 
243
 
 
244
  // find smallest cutghost
 
245
  double cutghostmin = comm->cutghost[0];
 
246
  if (cutghostmin > comm->cutghost[1])
 
247
    cutghostmin = comm->cutghost[1];
 
248
  if (cutghostmin > comm->cutghost[2])
 
249
    cutghostmin = comm->cutghost[2];
 
250
 
 
251
  // reset cutghost if needed
 
252
  if (cutneighmax_srp > cutghostmin){
 
253
    if(comm->me == 0){
 
254
      sprintf(str, "Extending ghost comm cutoff. New %f, old %f.", cutneighmax_srp, cutghostmin);
 
255
      error->message(FLERR,str);
 
256
    }
 
257
    // cutghost updated by comm->setup 
 
258
    comm->cutghostuser = cutneighmax_srp;
 
259
  }
 
260
 
 
261
  // put new particles in the box before exchange
 
262
  // move owned to new procs
 
263
  // get ghosts
 
264
  // build neigh lists again
 
265
  
 
266
  // if triclinic, lambda coords needed for pbc, exchange, borders
 
267
  if (domain->triclinic) domain->x2lamda(atom->nlocal);
 
268
  domain->pbc();
 
269
  comm->setup();
 
270
  if (neighbor->style) neighbor->setup_bins();
 
271
  comm->exchange();
 
272
  if (atom->sortfreq > 0) atom->sort();
 
273
  comm->borders();
 
274
  // back to box coords
 
275
  if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
 
276
  domain->image_check();
 
277
  domain->box_too_small_check();
 
278
  modify->setup_pre_neighbor();
 
279
  neighbor->build();
 
280
  neighbor->ncalls = 0;
 
281
 
 
282
  // new atom counts
 
283
  nlocal = atom->nlocal;
 
284
  nall = atom->nlocal + atom->nghost;
 
285
  // zero all forces 
 
286
  for(int i = 0; i < nall; i++)
 
287
    for(int n = 0; n < 3; n++)
 
288
      atom->f[i][n] = 0.0;
 
289
 
 
290
  // do not include bond particles in thermo output
 
291
  // remove them from all groups
 
292
  for(int i=0; i< nlocal; i++)
 
293
    if(atom->type[i] == bptype)
 
294
      atom->mask[i] = 0;
 
295
}
 
296
 
 
297
/* ----------------------------------------------------------------------
 
298
   set position of bond particles
 
299
------------------------------------------------------------------------- */
 
300
 
 
301
void FixSRP::pre_exchange()
 
302
{
 
303
  // update ghosts 
 
304
  comm->forward_comm();
 
305
 
 
306
  // reassign bond particle coordinates to midpoint of bonds
 
307
  // only need to do this before neigh rebuild
 
308
  double **x=atom->x;
 
309
  int i,j;
 
310
  int nlocal = atom->nlocal;
 
311
 
 
312
  for(int ii = 0; ii < nlocal; ii++){
 
313
    if(atom->type[ii] != bptype) continue;
 
314
 
 
315
    i = atom->map((int)array[ii][0]);
 
316
    if(i < 0) error->all(FLERR,"Fix SRP failed to map atom");
 
317
    i = domain->closest_image(ii,i);
 
318
 
 
319
    j = atom->map((int)array[ii][1]);
 
320
    if(j < 0) error->all(FLERR,"Fix SRP failed to map atom");
 
321
    j = domain->closest_image(ii,j);
 
322
 
 
323
    // position of bond particle ii 
 
324
    atom->x[ii][0] = (x[i][0] + x[j][0])*0.5;
 
325
    atom->x[ii][1] = (x[i][1] + x[j][1])*0.5;
 
326
    atom->x[ii][2] = (x[i][2] + x[j][2])*0.5;
 
327
  }
 
328
}
 
329
 
 
330
/* ----------------------------------------------------------------------
 
331
   memory usage of local atom-based array
 
332
------------------------------------------------------------------------- */
 
333
 
 
334
double FixSRP::memory_usage()
 
335
{
 
336
  double bytes = atom->nmax*2 * sizeof(double);
 
337
  return bytes;
 
338
}
 
339
 
 
340
/* ----------------------------------------------------------------------
 
341
   allocate atom-based array
 
342
------------------------------------------------------------------------- */
 
343
 
 
344
void FixSRP::grow_arrays(int nmax)
 
345
{
 
346
  memory->grow(array,nmax,2,"fix_srp:array");
 
347
  array_atom = array;
 
348
}
 
349
 
 
350
/* ----------------------------------------------------------------------
 
351
   copy values within local atom-based array
 
352
   called when move to new proc
 
353
------------------------------------------------------------------------- */
 
354
 
 
355
void FixSRP::copy_arrays(int i, int j, int delflag)
 
356
{
 
357
  for (int m = 0; m < 2; m++)
 
358
    array[j][m] = array[i][m];
 
359
}
 
360
 
 
361
/* ----------------------------------------------------------------------
 
362
   initialize one atom's array values
 
363
   called when atom is created
 
364
------------------------------------------------------------------------- */
 
365
 
 
366
void FixSRP::set_arrays(int i)
 
367
{
 
368
  array[i][0] = -1;
 
369
  array[i][1] = -1; 
 
370
}
 
371
 
 
372
/* ----------------------------------------------------------------------
 
373
   pack values in local atom-based array for exchange with another proc
 
374
------------------------------------------------------------------------- */
 
375
 
 
376
int FixSRP::pack_exchange(int i, double *buf)
 
377
{
 
378
  for (int m = 0; m < 2; m++) buf[m] = array[i][m];
 
379
  return 2;
 
380
}
 
381
 
 
382
/* ----------------------------------------------------------------------
 
383
   unpack values in local atom-based array from exchange with another proc
 
384
------------------------------------------------------------------------- */
 
385
 
 
386
int FixSRP::unpack_exchange(int nlocal, double *buf)
 
387
{
 
388
  for (int m = 0; m < 2; m++) array[nlocal][m] = buf[m];
 
389
  return 2;
 
390
}
 
391
/* ----------------------------------------------------------------------
 
392
   pack values for border communication at re-neighboring
 
393
------------------------------------------------------------------------- */
 
394
 
 
395
int FixSRP::pack_border(int n, int *list, double *buf)
 
396
{
 
397
  // pack buf for border com
 
398
  int i,j;
 
399
  int m = 0;
 
400
    for (i = 0; i < n; i++) {
 
401
      j = list[i];
 
402
      buf[m++] = array[j][0];
 
403
      buf[m++] = array[j][1];
 
404
    }
 
405
  return m;
 
406
}
 
407
 
 
408
/* ----------------------------------------------------------------------
 
409
   unpack values for border communication at re-neighboring
 
410
------------------------------------------------------------------------- */
 
411
 
 
412
int FixSRP::unpack_border(int n, int first, double *buf)
 
413
{
 
414
  // unpack buf into array
 
415
  int i,last;
 
416
  int m = 0;
 
417
  last = first + n;
 
418
 
 
419
      for (i = first; i < last; i++){
 
420
        array[i][0] = static_cast<int> (buf[m++]);
 
421
        array[i][1] = static_cast<int> (buf[m++]);
 
422
      }
 
423
  return m;
 
424
}
 
425
 
 
426
/* ----------------------------------------------------------------------
 
427
   remove particles after run
 
428
------------------------------------------------------------------------- */
 
429
 
 
430
void FixSRP::post_run()
 
431
{
 
432
  // all bond particles are removed after each run
 
433
  // useful for write_data and write_restart commands 
 
434
  // since those commands occur between runs
 
435
 
 
436
  bigint natoms_previous = atom->natoms;
 
437
  int nlocal = atom->nlocal;
 
438
  int* dlist;
 
439
  memory->create(dlist,nlocal,"fix_srp:dlist");
 
440
 
 
441
  for (int i = 0; i < nlocal; i++){
 
442
    if(atom->type[i] == bptype)
 
443
      dlist[i] = 1;
 
444
    else
 
445
      dlist[i] = 0;
 
446
  }
 
447
 
 
448
  // delete local atoms flagged in dlist
 
449
  // reset nlocal
 
450
 
 
451
  AtomVec *avec = atom->avec;
 
452
 
 
453
  int i = 0;
 
454
  while (i < nlocal) {
 
455
    if (dlist[i]) {
 
456
      avec->copy(nlocal-1,i,1);
 
457
      dlist[i] = dlist[nlocal-1];
 
458
      nlocal--;
 
459
    } else i++;
 
460
  }
 
461
 
 
462
  atom->nlocal = nlocal;
 
463
  memory->destroy(dlist);
 
464
 
 
465
  // reset atom->natoms
 
466
  // reset atom->map if it exists
 
467
  // set nghost to 0 so old ghosts won't be mapped
 
468
 
 
469
  bigint nblocal = atom->nlocal;
 
470
  MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
 
471
  if (atom->map_style) {
 
472
    atom->nghost = 0;
 
473
    atom->map_init();
 
474
    atom->map_set();
 
475
  }
 
476
 
 
477
  // print before and after atom count
 
478
 
 
479
  bigint ndelete = natoms_previous - atom->natoms;
 
480
 
 
481
  if (comm->me == 0) {
 
482
    if (screen) fprintf(screen,"Deleted " BIGINT_FORMAT
 
483
                        " atoms, new total = " BIGINT_FORMAT "\n",
 
484
                        ndelete,atom->natoms);
 
485
    if (logfile) fprintf(logfile,"Deleted " BIGINT_FORMAT
 
486
                         " atoms, new total = " BIGINT_FORMAT "\n",
 
487
                         ndelete,atom->natoms);
 
488
  }
 
489
 
 
490
  // verlet calls box_too_small_check() in post_run
 
491
  // this check maps all bond partners
 
492
  // therefore need ghosts
 
493
 
 
494
  // need to convert to lambda coords before apply pbc
 
495
  if (domain->triclinic) domain->x2lamda(atom->nlocal);
 
496
  domain->pbc();
 
497
  comm->setup();
 
498
  comm->exchange();
 
499
  if (atom->sortfreq > 0) atom->sort();
 
500
  comm->borders();
 
501
  // change back to box coordinates
 
502
  if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
 
503
} 
 
504
 
 
505
/* ----------------------------------------------------------------------
 
506
   pack values in local atom-based arrays for restart file
 
507
------------------------------------------------------------------------- */
 
508
 
 
509
int FixSRP::pack_restart(int i, double *buf)
 
510
{
 
511
  int m = 0;
 
512
  buf[m++] = 3;
 
513
  buf[m++] = array[i][0]; 
 
514
  buf[m++] = array[i][1];
 
515
  return m;
 
516
}
 
517
 
 
518
/* ----------------------------------------------------------------------
 
519
   unpack values from atom->extra array to restart the fix
 
520
------------------------------------------------------------------------- */
 
521
 
 
522
void FixSRP::unpack_restart(int nlocal, int nth) 
 
523
{
 
524
  double **extra = atom->extra;
 
525
 
 
526
// skip to Nth set of extra values
 
527
  int m = 0; 
 
528
  for (int i = 0; i < nth; i++){
 
529
    m += static_cast<int> (extra[nlocal][m]);
 
530
  }
 
531
 
 
532
  m++; 
 
533
  array[nlocal][0] = extra[nlocal][m++];
 
534
  array[nlocal][1] = extra[nlocal][m++];
 
535
 
 
536
}
 
537
/* ----------------------------------------------------------------------
 
538
   maxsize of any atom's restart data
 
539
------------------------------------------------------------------------- */
 
540
 
 
541
int FixSRP::maxsize_restart()
 
542
{
 
543
  return 3;
 
544
}
 
545
 
 
546
/* ----------------------------------------------------------------------
 
547
   size of atom nlocal's restart data
 
548
------------------------------------------------------------------------- */
 
549
 
 
550
int FixSRP::size_restart(int nlocal)
 
551
{
 
552
  return 3;
 
553
}
 
554
 
 
555
/* ----------------------------------------------------------------------
 
556
   pack global state of Fix 
 
557
------------------------------------------------------------------------- */
 
558
 
 
559
void FixSRP::write_restart(FILE *fp)
 
560
{
 
561
  int n = 0;
 
562
  double list[3];
 
563
  list[n++] = comm->cutghostuser;
 
564
  list[n++] = btype;
 
565
  list[n++] = bptype;
 
566
 
 
567
  if (comm->me == 0) {
 
568
    int size = n * sizeof(double);
 
569
    fwrite(&size,sizeof(int),1,fp);
 
570
    fwrite(list,sizeof(double),n,fp);
 
571
  }
 
572
}
 
573
 
 
574
/* ----------------------------------------------------------------------
 
575
   use info from restart file to restart the Fix
 
576
------------------------------------------------------------------------- */
 
577
 
 
578
void FixSRP::restart(char *buf)
 
579
{
 
580
  int n = 0;
 
581
  double *list = (double *) buf;
 
582
 
 
583
  comm->cutghostuser = static_cast<double> (list[n++]);
 
584
  btype = static_cast<int> (list[n++]);
 
585
  bptype = static_cast<int> (list[n++]);
 
586
}
 
587
 
 
588
/* ----------------------------------------------------------------------
 
589
   interface with pair class
 
590
   pair srp sets the bond type in this fix
 
591
------------------------------------------------------------------------- */
 
592
 
 
593
int FixSRP::modify_param(int narg, char **arg)
 
594
{
 
595
  if (strcmp(arg[0],"btype") == 0) {
 
596
    btype = atoi(arg[1]);
 
597
    return 2;
 
598
  }
 
599
  return 0;
 
600
}
 
601