← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/USER-MISC/pair_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/pair_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
This pair style srp command calculates a segmental repulsive force
 
18
between bonds. This is useful for preventing the crossing of bonds if
 
19
soft non-bonded potentials are used, such as DPD polymer chains.
 
20
 
 
21
See the doc page for pair_style srp command for usage instructions.
 
22
 
 
23
There is an example script for this package in examples/USER/srp.
 
24
 
 
25
Please contact Timothy Sirk for questions (tim.sirk@us.army.mil).
 
26
------------------------------------------------------------------------- */
 
27
 
 
28
#include "stdlib.h" 
 
29
#include "pair_srp.h" 
 
30
#include "atom.h" 
 
31
#include "comm.h" 
 
32
#include "force.h" 
 
33
#include "neighbor.h" 
 
34
#include "neigh_list.h" 
 
35
#include "memory.h" 
 
36
#include "error.h" 
 
37
#include "domain.h" 
 
38
#include "modify.h"
 
39
#include "fix.h"
 
40
#include "fix_srp.h"
 
41
#include "thermo.h"
 
42
#include "output.h"
 
43
#include "string.h"
 
44
#include "citeme.h"
 
45
 
 
46
using namespace LAMMPS_NS; 
 
47
 
 
48
#define SMALL 1.0e-10
 
49
#define BIG 1e10
 
50
#define ONETWOBIT 0x40000000
 
51
 
 
52
static const char cite_srp[] =
 
53
  "@Article{Sirk2012\n"
 
54
  " author = {T. Sirk and Y. Sliozberg and J. Brennan and M. Lisal and J. Andzelm},\n"
 
55
  " title = {An enhanced entangled polymer model for dissipative particle dynamics},\n"
 
56
  " journal = {J.~Chem.~Phys.},\n"
 
57
  " year =    2012,\n"
 
58
  " volume =  136,\n"
 
59
  " pages =   {134903}\n"
 
60
  "}\n\n";
 
61
 
 
62
/* ----------------------------------------------------------------------
 
63
 set size of pair comms in constructor
 
64
 ---------------------------------------------------------------------- */
 
65
 
 
66
PairSRP::PairSRP(LAMMPS *lmp) : Pair(lmp) 
 
67
{
 
68
  writedata = 1; 
 
69
 
 
70
  if (lmp->citeme) lmp->citeme->add(cite_srp);
 
71
 
 
72
  nextra = 1;
 
73
  segment = NULL;
 
74
} 
 
75
 
 
76
/* ----------------------------------------------------------------------
 
77
 allocate all arrays
 
78
 ------------------------------------------------------------------------- */
 
79
 
 
80
void PairSRP::allocate()
 
81
{
 
82
    allocated = 1;
 
83
    // particles of bptype inserted by fix srp
 
84
    // bptype is the highest numbered atom type
 
85
    int n = bptype;
 
86
    memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
 
87
    memory->create(cut, n + 1, n + 1, "pair:cut");
 
88
    memory->create(a0, n + 1, n + 1, "pair:a0");
 
89
 
 
90
    // setflag for atom types
 
91
    memory->create(setflag,n+1,n+1,"pair:setflag");
 
92
    for (int i = 1; i <= n; i++)
 
93
        for (int j = i; j <= n; j++)
 
94
            setflag[i][j] = 0;
 
95
 
 
96
    maxcount = 0;
 
97
}
 
98
 
 
99
/* ---------------------------------------------------------------------- 
 
100
 free 
 
101
 ------------------------------------------------------------------------- */ 
 
102
 
 
103
PairSRP::~PairSRP() 
 
104
{ 
 
105
    if (allocated) 
 
106
    { 
 
107
        memory->destroy(setflag); 
 
108
        memory->destroy(cutsq); 
 
109
        memory->destroy(cut); 
 
110
        memory->destroy(a0); 
 
111
        memory->destroy(segment); 
 
112
    }
 
113
 
 
114
  // check nfix in case all fixes have already been deleted
 
115
  if (modify->nfix) modify->delete_fix("mysrp");
 
116
} 
 
117
 
 
118
/* ---------------------------------------------------------------------- 
 
119
 compute bond-bond repulsions 
 
120
 ------------------------------------------------------------------------- */ 
 
121
 
 
122
void PairSRP::compute(int eflag, int vflag) 
 
123
 
 
124
{
 
125
    // setup energy and virial 
 
126
    if (eflag || vflag) 
 
127
        ev_setup(eflag, vflag); 
 
128
    else 
 
129
        evflag = vflag_fdotr = 0; 
 
130
 
 
131
    double **x = atom->x; 
 
132
    double **f = atom->f; 
 
133
    int nlocal = atom->nlocal; 
 
134
    int nall = nlocal + atom->nghost; 
 
135
    int i0, i1, j0, j1; 
 
136
    int i,j,ii,jj,inum,jnum;
 
137
    double dijsq, dij;
 
138
 
 
139
    int *ilist,*jlist,*numneigh,**firstneigh;
 
140
    inum = list->inum;
 
141
    ilist = list->ilist;
 
142
    numneigh = list->numneigh;
 
143
    firstneigh = list->firstneigh;
 
144
 
 
145
    double dx,dy,dz,ti,tj;
 
146
    double wd, lever0, lever1, evdwl, fpair;
 
147
    double fxlever0, fylever0, fzlever0, fxlever1, fylever1, fzlever1;
 
148
    double fx, fy, fz;
 
149
 
 
150
    // mapping global to local for atoms inside bond particles 
 
151
    // exclude 1-2 neighs if requested 
 
152
    if (neighbor->ago == 0){
 
153
      remapBonds(nall);
 
154
      if(exclude) onetwoexclude(ilist, inum, jlist, numneigh, firstneigh);
 
155
    }
 
156
 
 
157
  // this pair style only used with hybrid
 
158
  // due to exclusions
 
159
  // each atom i is type bptype  
 
160
  // each neigh j is type bptype 
 
161
  
 
162
  // using midpoint distance option
 
163
  if(midpoint){
 
164
 
 
165
    for (ii = 0; ii < inum; ii++) {
 
166
 
 
167
      i = ilist[ii];
 
168
      jnum = numneigh[i];
 
169
      // two atoms inside bond particle
 
170
      i0 = segment[i][0];
 
171
      j0 = segment[i][1];
 
172
 
 
173
      for (jj = 0; jj < jnum; jj++) {
 
174
 
 
175
        jlist = firstneigh[i];
 
176
        j = jlist[jj];
 
177
 
 
178
        // enforce 1-2 exclusions
 
179
        if( (sbmask(j) & exclude) )
 
180
          continue;
 
181
 
 
182
        j &= NEIGHMASK;
 
183
        //retrieve atoms from bond particle
 
184
        i1 = segment[j][0];
 
185
        j1 = segment[j][1];
 
186
 
 
187
        // midpt dist bond 0 and 1
 
188
        dx = 0.5*(x[i0][0] - x[i1][0] + x[j0][0] - x[j1][0]);
 
189
        dy = 0.5*(x[i0][1] - x[i1][1] + x[j0][1] - x[j1][1]);  
 
190
        dz = 0.5*(x[i0][2] - x[i1][2] + x[j0][2] - x[j1][2]);  
 
191
        dijsq = dx*dx + dy*dy + dz*dz;
 
192
 
 
193
        if (dijsq < cutsq[bptype][bptype]){
 
194
        dij = sqrt(dijsq);
 
195
 
 
196
        if (dij < SMALL) 
 
197
          continue;     // dij can be 0.0 with soft potentials
 
198
 
 
199
        wd = 1.0 - dij / cut[bptype][bptype];
 
200
        fpair = 0.5 * a0[bptype][bptype] * wd / dij; // 0.5 factor for lever rule 
 
201
 
 
202
        // force for bond 0, beads 0,1
 
203
        //force between bonds
 
204
        fx = fpair * dx;
 
205
        fy = fpair * dy;
 
206
        fz = fpair * dz;
 
207
 
 
208
        f[i0][0] += fx; //keep force sign for bond 0
 
209
        f[i0][1] += fy;
 
210
        f[i0][2] += fz;
 
211
 
 
212
        f[j0][0] += fx;
 
213
        f[j0][1] += fy;
 
214
        f[j0][2] += fz;
 
215
 
 
216
        f[i1][0] -= fx; //flip force sign for bond 1
 
217
        f[i1][1] -= fy;
 
218
        f[i1][2] -= fz;
 
219
 
 
220
        f[j1][0] -= fx;
 
221
        f[j1][1] -= fy;
 
222
        f[j1][2] -= fz;
 
223
 
 
224
        // ************************************************* //
 
225
 
 
226
        if (eflag){
 
227
          evdwl = 0.5 * a0[bptype][bptype] * cut[bptype][bptype] * wd * wd;
 
228
        }
 
229
 
 
230
        if (evflag){
 
231
          ev_tally(i0,i1,nlocal,1,0.5*evdwl,0.0,fpair,dx,dy,dz);
 
232
          ev_tally(j0,j1,nlocal,1,0.5*evdwl,0.0,fpair,dx,dy,dz);
 
233
        }
 
234
 
 
235
        if (vflag_fdotr) virial_fdotr_compute();
 
236
 
 
237
        }
 
238
      }
 
239
   }
 
240
 } 
 
241
  else{
 
242
  // using min distance option
 
243
  
 
244
    for (ii = 0; ii < inum; ii++) {
 
245
 
 
246
      i = ilist[ii];
 
247
      jnum = numneigh[i];
 
248
      i0 = segment[i][0];
 
249
      j0 = segment[i][1];
 
250
 
 
251
      for (jj = 0; jj < jnum; jj++) {
 
252
 
 
253
        jlist = firstneigh[i];
 
254
        j = jlist[jj];
 
255
 
 
256
        // enforce 1-2 exclusions
 
257
        if( (sbmask(j) & exclude) )
 
258
          continue;
 
259
 
 
260
        j &= NEIGHMASK;
 
261
 
 
262
        i1 = segment[j][0];
 
263
        j1 = segment[j][1];
 
264
 
 
265
        getMinDist(x, dx, dy, dz, ti, tj, i0, j0, i1, j1);
 
266
        dijsq = dx*dx + dy*dy + dz*dz;
 
267
 
 
268
        if (dijsq < cutsq[bptype][bptype]){
 
269
 
 
270
        dij = sqrt(dijsq); 
 
271
 
 
272
        if (dij < SMALL)
 
273
          continue;     // dij can be 0.0 with soft potentials
 
274
 
 
275
        wd = 1.0 - dij / cut[bptype][bptype];
 
276
        fpair = a0[bptype][bptype] * wd / dij; 
 
277
 
 
278
        // force for bond 0, beads 0,1
 
279
        lever0 = 0.5 + ti; // assign force according to lever rule
 
280
        lever1 = 0.5 + tj; // assign force according to lever rule
 
281
        //force between bonds
 
282
        fx = fpair * dx;
 
283
        fy = fpair * dy;
 
284
        fz = fpair * dz;
 
285
 
 
286
        //decompose onto atoms
 
287
        fxlever0 = fx * lever0;
 
288
        fylever0 = fy * lever0;
 
289
        fzlever0 = fz * lever0;
 
290
        fxlever1 = fx * lever1;
 
291
        fylever1 = fy * lever1;
 
292
        fzlever1 = fz * lever1;
 
293
 
 
294
        f[i0][0] += fxlever0; //keep force sign for bond 0
 
295
        f[i0][1] += fylever0;
 
296
        f[i0][2] += fzlever0;
 
297
 
 
298
        f[j0][0] += (fx - fxlever0);
 
299
        f[j0][1] += (fy - fylever0);
 
300
        f[j0][2] += (fz - fzlever0);
 
301
 
 
302
        f[i1][0] -= fxlever1; //flip force sign for bond 1
 
303
        f[i1][1] -= fylever1;
 
304
        f[i1][2] -= fzlever1;
 
305
 
 
306
        f[j1][0] -= (fx - fxlever1);
 
307
        f[j1][1] -= (fy - fylever1);
 
308
        f[j1][2] -= (fz - fzlever1);
 
309
 
 
310
        // ************************************************* //
 
311
 
 
312
        if (eflag){
 
313
          evdwl = 0.5 * a0[bptype][bptype] * cut[bptype][bptype] * wd * wd;
 
314
        }
 
315
 
 
316
        if (evflag){
 
317
          ev_tally(i0,i1,nlocal,1,0.5*evdwl,0.0,0.5*fpair,dx,dy,dz);
 
318
          ev_tally(j0,j1,nlocal,1,0.5*evdwl,0.0,0.5*fpair,dx,dy,dz);
 
319
        }
 
320
 
 
321
       if (vflag_fdotr) virial_fdotr_compute();
 
322
 
 
323
      }
 
324
    }
 
325
  }
 
326
 }
 
327
}
 
328
 
 
329
/* ----------------------------------------------------------------------
 
330
 global settings
 
331
 ------------------------------------------------------------------------- */
 
332
 
 
333
void PairSRP::settings(int narg, char **arg)
 
334
{
 
335
    if (narg < 3 || narg > 5)
 
336
        error->all(FLERR,"Illegal pair_style command");
 
337
 
 
338
    cut_global = force->numeric(FLERR,arg[0]);
 
339
    btype = force->inumeric(FLERR,arg[1]);
 
340
    if (btype > atom->nbondtypes) error->all(FLERR,"Illegal pair_style command");
 
341
 
 
342
    // settings
 
343
    midpoint = 0;
 
344
    min = 0;
 
345
 
 
346
  if (strcmp(arg[2],"min") == 0) min = 1;
 
347
  else if (strcmp(arg[2],"mid") == 0) midpoint = 1;
 
348
  else
 
349
     error->all(FLERR,"Illegal pair_style command");
 
350
 
 
351
  int iarg = 3;
 
352
  // default exclude 1-2
 
353
  // scaling for 1-2, etc not supported
 
354
  exclude = 1; 
 
355
 
 
356
  while (iarg < narg) {
 
357
    if (strcmp(arg[iarg],"exclude") == 0) {
 
358
      if (iarg+2 > narg) error->all(FLERR,"Illegal pair srp command");
 
359
      if (strcmp(arg[iarg+1],"yes") == 0)
 
360
        exclude = 1;
 
361
      if (strcmp(arg[iarg+1],"no") == 0){
 
362
        if (min) error->all(FLERR,"Illegal exclude option in pair srp command");
 
363
        exclude = 0;
 
364
      }
 
365
      iarg += 2;
 
366
    } else error->all(FLERR,"Illegal pair srp command");
 
367
  }
 
368
 
 
369
  // highest atom type is saved for bond particles
 
370
  bptype = atom->ntypes;
 
371
 
 
372
  // reset cutoffs if explicitly set
 
373
  if (allocated) {
 
374
    int i,j;
 
375
    for (i = 1; i <= bptype; i++)
 
376
      for (j = i+1; j <= bptype; j++)
 
377
        if (setflag[i][j]) cut[i][j] = cut_global;
 
378
  }
 
379
}
 
380
 
 
381
/* ----------------------------------------------------------------------
 
382
 set coeffs 
 
383
 ------------------------------------------------------------------------- */
 
384
 
 
385
void PairSRP::coeff(int narg, char **arg)
 
386
{
 
387
    if (narg < 3 || narg > 4)
 
388
        error->all(FLERR,"PairSRP: Incorrect args for pair coeff");
 
389
    if (!allocated) allocate();
 
390
 
 
391
    // set ij bond-bond cutoffs
 
392
    int ilo, ihi, jlo, jhi;
 
393
    force->bounds(arg[0], bptype, ilo, ihi);
 
394
    force->bounds(arg[1], bptype, jlo, jhi);
 
395
 
 
396
    double a0_one = force->numeric(FLERR,arg[2]);
 
397
    double cut_one = cut_global;
 
398
    if (narg == 4)  cut_one = force->numeric(FLERR,arg[3]);
 
399
 
 
400
    int count = 0;
 
401
    for (int i = ilo; i <= ihi; i++)
 
402
    {
 
403
        for (int j = MAX(jlo,i); j <= jhi; j++)
 
404
        {
 
405
            a0[i][j] = a0_one;
 
406
            cut[i][j] = cut_one;
 
407
            cutsq[i][j] = cut_one * cut_one;
 
408
            setflag[i][j] = 1;
 
409
            count++;
 
410
        }
 
411
    }
 
412
 
 
413
    if (count == 0) error->warning(FLERR,"PairSRP: No pair coefficients were set");
 
414
}
 
415
 
 
416
/* ----------------------------------------------------------------------
 
417
 init specific to this pair style
 
418
 ------------------------------------------------------------------------- */
 
419
 
 
420
void PairSRP::init_style()
 
421
{
 
422
  if (!force->newton_pair)
 
423
    error->all(FLERR,"PairSRP: Pair srp requires newton pair on");
 
424
 
 
425
  // need fix srp
 
426
  // invoke here instead of constructor
 
427
  // to make restart possible
 
428
  char **fixarg = new char*[3];
 
429
  fixarg[0] = (char *) "mysrp";
 
430
  fixarg[1] = (char *) "all";
 
431
  fixarg[2] = (char *) "SRP";
 
432
  modify->add_fix(3,fixarg);
 
433
  f_srp = (FixSRP *) modify->fix[modify->nfix-1];
 
434
  delete [] fixarg;
 
435
 
 
436
  // set bond type in fix srp
 
437
  // bonds of this type will be represented by bond particles
 
438
  // btype = bond type
 
439
  // bptype = bond particle type
 
440
  char c0[20];
 
441
  char* arg0[2];
 
442
  sprintf(c0, "%d", btype);
 
443
  arg0[0] = (char *) "btype";
 
444
  arg0[1] = c0;
 
445
  f_srp->modify_params(2, arg0);
 
446
 
 
447
  // bond particles do not contribute to energy or virial
 
448
  // bond particles do not belong to group all
 
449
  // but thermo normalization is by nall 
 
450
  // therefore should turn off normalization
 
451
  int me;
 
452
  MPI_Comm_rank(world,&me);
 
453
  char *arg1[2];
 
454
  arg1[0] = (char *) "norm";
 
455
  arg1[1] = (char *) "no";
 
456
  output->thermo->modify_params(2, arg1);
 
457
  if (me == 0) 
 
458
    error->message(FLERR,"Thermo normalization turned off by pair srp");
 
459
 
 
460
  neighbor->request(this,instance_me);
 
461
}
 
462
 
 
463
/* ----------------------------------------------------------------------
 
464
 init for one type pair i,j and corresponding j,i
 
465
 ------------------------------------------------------------------------- */
 
466
 
 
467
double PairSRP::init_one(int i, int j)
 
468
{
 
469
 
 
470
 if (setflag[i][j] == 0) error->all(FLERR,"PairSRP: All pair coeffs are not set");
 
471
 
 
472
  cut[j][i] = cut[i][j];
 
473
  a0[j][i] = a0[i][j];
 
474
 
 
475
  return cut[i][j];
 
476
}
 
477
 
 
478
/* ---------------------------------------------------------------------- 
 
479
 find min distance for bonds i0/j0 and i1/j1
 
480
 ------------------------------------------------------------------------- */ 
 
481
inline void PairSRP::getMinDist(double** &x, double &dx, double &dy, double &dz, double &ti, double &tj, int &i0, int &j0, int &i1, int &j1)
 
482
{ 
 
483
    // move these outside the loop 
 
484
    double diffx0, diffy0, diffz0, diffx1, diffy1, diffz1, dPx, dPy, dPz, RiRi, RiRj, RjRj;
 
485
    double denom, termx0, termy0, termz0, num0, termx1, termy1, termz1, num1;
 
486
 
 
487
    // compute midpt dist from 1st atom, 1st bond
 
488
    diffx0 = x[j0][0] - x[i0][0]; // x,y,z from bond 0
 
489
    diffy0 = x[j0][1] - x[i0][1];
 
490
    diffz0 = x[j0][2] - x[i0][2];
 
491
 
 
492
    // compute midpt dist from 1st atom, 2nd bond
 
493
    diffx1 = x[j1][0] - x[i1][0]; 
 
494
    diffy1 = x[j1][1] - x[i1][1];
 
495
    diffz1 = x[j1][2] - x[i1][2];
 
496
 
 
497
    // midpoint distance
 
498
    dPx = 0.5*(diffx0-diffx1) + x[i0][0]-x[i1][0];
 
499
    dPy = 0.5*(diffy0-diffy1) + x[i0][1]-x[i1][1];
 
500
    dPz = 0.5*(diffz0-diffz1) + x[i0][2]-x[i1][2];
 
501
 
 
502
    // Ri^2 Rj^2
 
503
    RiRi = diffx0*diffx0 + diffy0*diffy0 + diffz0*diffz0;
 
504
    RiRj = diffx0*diffx1 + diffy0*diffy1 + diffz0*diffz1;
 
505
    RjRj = diffx1*diffx1 + diffy1*diffy1 + diffz1*diffz1;
 
506
    denom = RiRj*RiRj - RiRi*RjRj;
 
507
 
 
508
    // handle case of parallel lines
 
509
    // reduce to midpt distance 
 
510
    if (fabs(denom) < SMALL){ 
 
511
        if(denom < 0) denom = -BIG;
 
512
        else denom = BIG;
 
513
    } 
 
514
 
 
515
    // calc ti  
 
516
    termx0 = RiRj*diffx1 - RjRj*diffx0;
 
517
    termy0 = RiRj*diffy1 - RjRj*diffy0;
 
518
    termz0 = RiRj*diffz1 - RjRj*diffz0;
 
519
    num0 = dPx*termx0 + dPy*termy0 + dPz*termz0;
 
520
    ti = num0 / denom;
 
521
    if (ti > 0.5) ti = 0.5; 
 
522
    if (ti < -0.5) ti = -0.5; 
 
523
 
 
524
    // calc tj  
 
525
    termx1 = RiRj*diffx0 - RiRi*diffx1;
 
526
    termy1 = RiRj*diffy0 - RiRi*diffy1;
 
527
    termz1 = RiRj*diffz0 - RiRi*diffz1;
 
528
    num1 = dPx*termx1 + dPy*termy1 + dPz*termz1;
 
529
    tj = -num1/ denom;
 
530
    if (tj > 0.5)  tj = 0.5;
 
531
    if (tj < -0.5) tj = -0.5;
 
532
 
 
533
    // min dist 
 
534
    dx = dPx - ti*diffx0 + tj*diffx1;
 
535
    dy = dPy - ti*diffy0 + tj*diffy1;
 
536
    dz = dPz - ti*diffz0 + tj*diffz1;
 
537
} 
 
538
 
 
539
/* -------------------------------------------------------- 
 
540
map global id of atoms in stored by each bond particle
 
541
 ------------------------------------------------------- */ 
 
542
inline void PairSRP::remapBonds(int &nall)
 
543
{
 
544
  if(nall > maxcount){
 
545
    memory->grow(segment, nall, 2, "pair:segment");    
 
546
    maxcount = nall;
 
547
  }
 
548
 
 
549
  // loop over all bond particles
 
550
  // each bond paricle holds two bond atoms
 
551
  // map global ids of bond atoms to local ids
 
552
  // might not be able to map both bond atoms of j, if j is outside neighcut  
 
553
  // these are not on neighlist, so are not used 
 
554
  int tmp;
 
555
  srp = f_srp->array_atom;
 
556
 
 
557
    for (int i = 0; i < nall; i++) {
 
558
      if(atom->type[i] == bptype){
 
559
        // tmp is local id
 
560
        // tmp == -1 is ok
 
561
        tmp = atom->map((int)srp[i][0]);
 
562
        segment[i][0] = domain->closest_image(i,tmp);
 
563
        // repeat with other id
 
564
        tmp = atom->map((int)srp[i][1]);
 
565
        segment[i][1] = domain->closest_image(i,tmp);
 
566
      }
 
567
    }
 
568
}
 
569
 
 
570
/* -------------------------------------------------------- 
 
571
add exclusions for 1-2 neighs, if requested
 
572
more complex exclusions or scaling probably not needed 
 
573
 ------------------------------------------------------- */ 
 
574
inline void PairSRP::onetwoexclude(int* &ilist, int &inum, int* &jlist, int* &numneigh, int** &firstneigh)
 
575
{
 
576
    int i0, i1, j0, j1;
 
577
    int i,j,ii,jj,jnum;
 
578
 
 
579
    // encode neighs with exclusions
 
580
    // only need 1-2 info for normal uses of srp 
 
581
    // add 1-3, etc later if ever needed 
 
582
 
 
583
    for (ii = 0; ii < inum; ii++) {
 
584
 
 
585
      i = ilist[ii];
 
586
      jnum = numneigh[i];
 
587
      // two atoms inside bond particle
 
588
      i0 = segment[i][0];
 
589
      j0 = segment[i][1];
 
590
 
 
591
      for (jj = 0; jj < jnum; jj++) {
 
592
 
 
593
        jlist = firstneigh[i];
 
594
        j = jlist[jj];
 
595
        j &= NEIGHMASK;
 
596
        //two atoms inside bond particle
 
597
        i1 = segment[j][0];
 
598
        j1 = segment[j][1];
 
599
 
 
600
        // check for a 1-2 neigh 
 
601
        if(i0 == i1 || i0 == j1 || i1 == j0 || j0 == j1){
 
602
          j |= ONETWOBIT;
 
603
          jlist[jj] = j;
 
604
        }
 
605
      }
 
606
    }
 
607
}
 
608
 
 
609
/* ----------------------------------------------------------------------
 
610
proc 0 writes to data file
 
611
------------------------------------------------------------------------- */
 
612
 
 
613
void PairSRP::write_data(FILE *fp)
 
614
{
 
615
  for (int i = 1; i <= atom->ntypes; i++)
 
616
    fprintf(fp,"%d %g\n",i,a0[i][i]);
 
617
}
 
618
 
 
619
/* ----------------------------------------------------------------------
 
620
proc 0 writes all pairs to data file
 
621
------------------------------------------------------------------------- */
 
622
 
 
623
void PairSRP::write_data_all(FILE *fp)
 
624
{
 
625
  for (int i = 1; i <= atom->ntypes; i++)
 
626
    for (int j = i; j <= atom->ntypes; j++)
 
627
      fprintf(fp,"%d %d %g %g\n",i,j,a0[i][j],cut[i][j]);
 
628
}
 
629
 
 
630
/* ----------------------------------------------------------------------
 
631
   proc 0 writes to restart file
 
632
------------------------------------------------------------------------- */
 
633
 
 
634
void PairSRP::write_restart(FILE *fp)
 
635
{
 
636
  write_restart_settings(fp);
 
637
 
 
638
  int i,j;
 
639
  for (i = 1; i <= atom->ntypes; i++)
 
640
    for (j = i; j <= atom->ntypes; j++) {
 
641
      fwrite(&setflag[i][j],sizeof(int),1,fp);
 
642
      if (setflag[i][j]) {
 
643
        fwrite(&a0[i][j],sizeof(double),1,fp);
 
644
        fwrite(&cut[i][j],sizeof(double),1,fp);
 
645
      }
 
646
    }
 
647
}
 
648
 
 
649
/* ----------------------------------------------------------------------
 
650
   proc 0 reads from restart file, bcasts
 
651
------------------------------------------------------------------------- */
 
652
 
 
653
void PairSRP::read_restart(FILE *fp)
 
654
{
 
655
  read_restart_settings(fp);
 
656
  allocate();
 
657
 
 
658
  int i,j;
 
659
  int me = comm->me;
 
660
  for (i = 1; i <= atom->ntypes; i++)
 
661
    for (j = i; j <= atom->ntypes; j++) {
 
662
      if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
 
663
      MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
 
664
      if (setflag[i][j]) {
 
665
        if (me == 0) {
 
666
          printf(" i %d j %d \n",i,j);
 
667
          fread(&a0[i][j],sizeof(double),1,fp);
 
668
          fread(&cut[i][j],sizeof(double),1,fp);
 
669
        }
 
670
        MPI_Bcast(&a0[i][j],1,MPI_DOUBLE,0,world);
 
671
        MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
 
672
      }
 
673
    }
 
674
}
 
675
/* ----------------------------------------------------------------------
 
676
   proc 0 writes to restart file
 
677
------------------------------------------------------------------------- */
 
678
 
 
679
void PairSRP::write_restart_settings(FILE *fp)
 
680
{
 
681
  fwrite(&cut_global,sizeof(double),1,fp);
 
682
  fwrite(&bptype,sizeof(int),1,fp);
 
683
  fwrite(&btype,sizeof(int),1,fp);
 
684
  fwrite(&min,sizeof(int),1,fp);
 
685
  fwrite(&midpoint,sizeof(int),1,fp);
 
686
  fwrite(&exclude,sizeof(int),1,fp);
 
687
}
 
688
 
 
689
/* ----------------------------------------------------------------------
 
690
   proc 0 reads from restart file, bcasts
 
691
------------------------------------------------------------------------- */
 
692
 
 
693
void PairSRP::read_restart_settings(FILE *fp)
 
694
{
 
695
  if (comm->me == 0) {
 
696
    fread(&cut_global,sizeof(double),1,fp);
 
697
    fread(&bptype,sizeof(int),1,fp);
 
698
    fread(&btype,sizeof(int),1,fp);
 
699
    fread(&min,sizeof(int),1,fp);
 
700
    fread(&midpoint,sizeof(int),1,fp);
 
701
    fread(&exclude,sizeof(int),1,fp);
 
702
  }
 
703
  MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
 
704
}