← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/pair_coul_streitz.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/pair_coul_streitz.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 author: Ray Shan (Sandia)
 
16
------------------------------------------------------------------------- */
 
17
 
 
18
#include "math.h"
 
19
#include "stdio.h"
 
20
#include "stdlib.h"
 
21
#include "string.h"
 
22
#include "pair_coul_streitz.h"
 
23
#include "atom.h"
 
24
#include "comm.h"
 
25
#include "force.h"
 
26
#include "kspace.h"
 
27
#include "neighbor.h"
 
28
#include "neigh_list.h"
 
29
#include "neigh_request.h"
 
30
#include "group.h"
 
31
#include "update.h"
 
32
#include "math_const.h"
 
33
#include "memory.h"
 
34
#include "error.h"
 
35
 
 
36
using namespace LAMMPS_NS;
 
37
using namespace MathConst;
 
38
 
 
39
#define MAXLINE 1024
 
40
#define DELTA 4
 
41
#define PGDELTA 1
 
42
#define MAXNEIGH 24
 
43
 
 
44
/* ---------------------------------------------------------------------- */
 
45
 
 
46
PairCoulStreitz::PairCoulStreitz(LAMMPS *lmp) : Pair(lmp)
 
47
{
 
48
  single_enable = 0;
 
49
  restartinfo = 0;
 
50
  one_coeff = 1;
 
51
  nmax = 0;
 
52
  nelements = 0;
 
53
 
 
54
  elements = NULL;
 
55
  nparams = 0;
 
56
  maxparam = 0;
 
57
  params = NULL;
 
58
  elem2param = NULL;
 
59
}
 
60
 
 
61
/* ----------------------------------------------------------------------
 
62
   check if allocated, since class can be destructed when incomplete
 
63
------------------------------------------------------------------------- */
 
64
 
 
65
PairCoulStreitz::~PairCoulStreitz()
 
66
{
 
67
  if (elements)
 
68
    for (int i = 0; i < nelements; i++) delete [] elements[i];
 
69
 
 
70
  delete [] elements;
 
71
  memory->sfree(params);
 
72
  memory->destroy(elem2param);
 
73
 
 
74
  if (allocated) {
 
75
    memory->destroy(setflag);
 
76
    memory->destroy(cutsq);
 
77
    memory->destroy(scale);
 
78
    memory->destroy(qeq_x);
 
79
    memory->destroy(qeq_j);
 
80
    memory->destroy(qeq_g);
 
81
    memory->destroy(qeq_z);
 
82
    memory->destroy(qeq_c);
 
83
    delete [] map;
 
84
  }
 
85
}
 
86
 
 
87
/* ---------------------------------------------------------------------- */
 
88
 
 
89
void PairCoulStreitz::allocate()
 
90
{
 
91
 allocated = 1;
 
92
 int n = atom->ntypes;
 
93
 
 
94
 memory->create(setflag,n+1,n+1,"pair:setflag");
 
95
 memory->create(cutsq,n+1,n+1,"pair:cutsq");
 
96
 memory->create(scale,n+1,n+1,"pair:scale");
 
97
 memory->create(qeq_x,n+1,"pair:qeq_x");
 
98
 memory->create(qeq_j,n+1,"pair:qeq_j");
 
99
 memory->create(qeq_g,n+1,"pair:qeq_g");
 
100
 memory->create(qeq_z,n+1,"pair:qeq_z");
 
101
 memory->create(qeq_c,n+1,"pair:qeq_c");
 
102
 
 
103
 map = new int[n+1];
 
104
}
 
105
 
 
106
/* ----------------------------------------------------------------------
 
107
   global settings
 
108
------------------------------------------------------------------------- */
 
109
 
 
110
void PairCoulStreitz::settings(int narg, char **arg)
 
111
{
 
112
  if (narg < 2) error->all(FLERR,"Illegal pair_style command");
 
113
 
 
114
  cut_coul = force->numeric(FLERR,arg[0]);
 
115
 
 
116
  if (strcmp(arg[1],"wolf") == 0){
 
117
    kspacetype = 1;
 
118
    g_wolf = force->numeric(FLERR,arg[2]);
 
119
  } else if (strcmp(arg[1],"ewald") == 0){
 
120
    ewaldflag = pppmflag = 1;
 
121
    kspacetype = 2;
 
122
  } else {
 
123
    error->all(FLERR,"Illegal pair_style command");
 
124
  }
 
125
}
 
126
 
 
127
/* ----------------------------------------------------------------------
 
128
   set coeffs for one or more type pairs
 
129
------------------------------------------------------------------------- */
 
130
 
 
131
void PairCoulStreitz::coeff(int narg, char **arg)
 
132
{
 
133
  int i,j,n;
 
134
 
 
135
  if (!allocated) allocate();
 
136
 
 
137
  if (narg != 3 + atom->ntypes)
 
138
    error->all(FLERR,"Incorrect args for pair coefficients");
 
139
 
 
140
  // insure I,J args are * *
 
141
 
 
142
  if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
 
143
    error->all(FLERR,"Incorrect args for pair coefficients");
 
144
 
 
145
  // read args that map atom types to elements in potential file
 
146
  // map[i] = which element the Ith atom type is, -1 if NULL
 
147
  // nelements = # of unique elements
 
148
  // elements = list of element names
 
149
 
 
150
  if (elements) {
 
151
    for (i = 0; i < nelements; i++) delete [] elements[i];
 
152
    delete [] elements;
 
153
  }
 
154
  elements = new char*[atom->ntypes];
 
155
  for (i = 0; i < atom->ntypes; i++) elements[i] = NULL;
 
156
 
 
157
  nelements = 0;
 
158
  for (i = 3; i < narg; i++) {
 
159
    if (strcmp(arg[i],"NULL") == 0) {
 
160
      map[i-2] = -1;
 
161
      continue;
 
162
    }
 
163
    for (j = 0; j < nelements; j++)
 
164
      if (strcmp(arg[i],elements[j]) == 0) break;
 
165
    map[i-2] = j;
 
166
    if (j == nelements) {
 
167
      n = strlen(arg[i]) + 1;
 
168
      elements[j] = new char[n];
 
169
      strcpy(elements[j],arg[i]);
 
170
      nelements++;
 
171
    }
 
172
  }
 
173
 
 
174
  // read potential file and initialize potential parameters
 
175
 
 
176
  read_file(arg[2]);
 
177
  setup();
 
178
  n = atom->ntypes;
 
179
 
 
180
  // clear setflag since coeff() called once with I,J = * *
 
181
 
 
182
  for (int i = 1; i <= n; i++)
 
183
    for (int j = i; j <= n; j++)
 
184
      setflag[i][j] = 0;
 
185
 
 
186
  // set setflag i,j for type pairs where both are mapped to elements
 
187
 
 
188
  int count = 0;
 
189
  for (int i = 1; i <= n; i++)
 
190
    for (int j = i; j <= n; j++)
 
191
      if (map[i] >= 0 && map[j] >= 0) {
 
192
        scale[i][j] = 1.0;
 
193
        setflag[i][j] = 1;
 
194
        count++;
 
195
      }
 
196
 
 
197
  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
 
198
}
 
199
 
 
200
/* ----------------------------------------------------------------------
 
201
   init specific to this pair style
 
202
------------------------------------------------------------------------- */
 
203
 
 
204
void PairCoulStreitz::init_style()
 
205
{
 
206
  if (!atom->q_flag)
 
207
    error->all(FLERR,"Pair style coul/sm requires atom attribute q");
 
208
 
 
209
  //neighbor->request(this);
 
210
  int irequest = neighbor->request(this,instance_me);
 
211
  neighbor->requests[irequest]->half = 0;
 
212
  neighbor->requests[irequest]->full = 1;
 
213
 
 
214
  cut_coulsq = cut_coul * cut_coul;
 
215
 
 
216
  // insure use of KSpace long-range solver when ewald specified, set g_ewald
 
217
 
 
218
  if (ewaldflag) {
 
219
    if (force->kspace == NULL)
 
220
      error->all(FLERR,"Pair style requires KSpace style ewald");
 
221
    g_ewald = force->kspace->g_ewald;
 
222
  }
 
223
 
 
224
  // ptr to QEQ fix
 
225
  //for (i = 0; i < modify->nfix; i++)
 
226
  //  if (strcmp(modify->fix[i]->style,"qeq") == 0) break;
 
227
  //if (i < modify->nfix) fixqeq = (FixQEQ *) modify->fix[i];
 
228
  //else fixqeq = NULL;
 
229
}
 
230
 
 
231
/* ----------------------------------------------------------------------
 
232
   init for one type pair i,j and corresponding j,i
 
233
------------------------------------------------------------------------- */
 
234
 
 
235
double PairCoulStreitz::init_one(int i, int j)
 
236
{
 
237
  scale[j][i] = scale[i][j];
 
238
 
 
239
  if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
 
240
  return cut_coul;
 
241
}
 
242
 
 
243
/* ---------------------------------------------------------------------- */
 
244
 
 
245
void PairCoulStreitz::read_file(char *file)
 
246
{
 
247
  int params_per_line = 6;
 
248
  char **words = new char*[params_per_line+1];
 
249
 
 
250
  memory->sfree(params);
 
251
  params = NULL;
 
252
  nparams = 0;
 
253
  maxparam = 0;
 
254
 
 
255
  // open file on proc 0
 
256
 
 
257
  FILE *fp;
 
258
  if (comm->me == 0) {
 
259
    fp = fopen(file,"r");
 
260
    if (fp == NULL) {
 
261
      char str[128];
 
262
      sprintf(str,"Cannot open coul/Streitz potential file %s",file);
 
263
      error->one(FLERR,str);
 
264
    }
 
265
  }
 
266
 
 
267
  // read each line out of file, skipping blank lines or leading '#'
 
268
  // store line of params if all 3 element tags are in element list
 
269
 
 
270
  int n,nwords,ielement;
 
271
  char line[MAXLINE],*ptr;
 
272
  int eof = 0;
 
273
 
 
274
  while (1) {
 
275
    if (comm->me == 0) {
 
276
      ptr = fgets(line,MAXLINE,fp);
 
277
      if (ptr == NULL) {
 
278
        eof = 1;
 
279
        fclose(fp);
 
280
      } else n = strlen(line) + 1;
 
281
    }
 
282
    MPI_Bcast(&eof,1,MPI_INT,0,world);
 
283
    if (eof) break;
 
284
    MPI_Bcast(&n,1,MPI_INT,0,world);
 
285
    MPI_Bcast(line,n,MPI_CHAR,0,world);
 
286
 
 
287
    // strip comment, skip line if blank
 
288
 
 
289
    if ((ptr = strchr(line,'#'))) *ptr = '\0';
 
290
    nwords = atom->count_words(line);
 
291
    if (nwords == 0) continue;
 
292
 
 
293
    // concatenate additional lines until have params_per_line words
 
294
 
 
295
    while (nwords < params_per_line) {
 
296
      n = strlen(line);
 
297
      if (comm->me == 0) {
 
298
        ptr = fgets(&line[n],MAXLINE-n,fp);
 
299
        if (ptr == NULL) {
 
300
          eof = 1;
 
301
          fclose(fp);
 
302
        } else n = strlen(line) + 1;
 
303
      }
 
304
      MPI_Bcast(&eof,1,MPI_INT,0,world);
 
305
      if (eof) break;
 
306
      MPI_Bcast(&n,1,MPI_INT,0,world);
 
307
      MPI_Bcast(line,n,MPI_CHAR,0,world);
 
308
      if ((ptr = strchr(line,'#'))) *ptr = '\0';
 
309
      nwords = atom->count_words(line);
 
310
    }
 
311
 
 
312
    if (nwords != params_per_line)
 
313
      error->all(FLERR,"Incorrect format in coul/Streitz potential file");
 
314
 
 
315
    // words = ptrs to all words in line
 
316
 
 
317
    nwords = 0;
 
318
    words[nwords++] = strtok(line," \t\n\r\f");
 
319
    while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
 
320
 
 
321
    // ielement = 1st args
 
322
 
 
323
    for (ielement = 0; ielement < nelements; ielement++)
 
324
      if (strcmp(words[0],elements[ielement]) == 0) break;
 
325
    if (ielement == nelements) continue;
 
326
 
 
327
    // load up parameter settings and error check their values
 
328
 
 
329
    if (nparams == maxparam) {
 
330
      maxparam += DELTA;
 
331
      params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
 
332
                                          "pair:params");
 
333
    }
 
334
 
 
335
    params[nparams].ielement = ielement;
 
336
    params[nparams].chi = atof(words[1]);
 
337
    params[nparams].eta = atof(words[2]);
 
338
    params[nparams].gamma = atof(words[3]);
 
339
    params[nparams].zeta = atof(words[4]);
 
340
    params[nparams].zcore = atof(words[5]);
 
341
 
 
342
    // parameter sanity check
 
343
    
 
344
    if (params[nparams].eta < 0.0 || params[nparams].zeta < 0.0 || 
 
345
        params[nparams].zcore < 0.0 || params[nparams].gamma != 0.0 )
 
346
      error->all(FLERR,"Illegal coul/Streitz parameter");
 
347
 
 
348
    nparams++;
 
349
  }
 
350
 
 
351
  delete [] words;
 
352
}
 
353
 
 
354
/* ---------------------------------------------------------------------- */
 
355
 
 
356
void PairCoulStreitz::setup()
 
357
{
 
358
  int i,m,n;
 
359
 
 
360
  // set elem2param 
 
361
 
 
362
  memory->destroy(elem2param);
 
363
  memory->create(elem2param,nelements,"pair:elem2param");
 
364
 
 
365
  for (i = 0; i < nelements; i++) {
 
366
    n = -1;
 
367
    for (m = 0; m < nparams; m++) {
 
368
      if (i == params[m].ielement ) {
 
369
        if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
 
370
        n = m;
 
371
      }
 
372
    }
 
373
    if (n < 0) error->all(FLERR,"Potential file is missing an entry");
 
374
    elem2param[i] = n;
 
375
  }
 
376
 
 
377
  // Wolf sum self energy
 
378
 
 
379
  if (kspacetype == 1) {
 
380
    double a = g_wolf;
 
381
    double r = cut_coul;
 
382
    double ar = a*r;
 
383
 
 
384
    woself = 0.50*erfc(ar)/r + a/MY_PIS;  // kc constant not yet multiplied
 
385
    dwoself = -(erfc(ar)/r/r + 2.0*a/MY_PIS*exp(-ar*ar)/r);
 
386
  }
 
387
}
 
388
 
 
389
/* ---------------------------------------------------------------------- */
 
390
 
 
391
void PairCoulStreitz::compute(int eflag, int vflag)
 
392
{
 
393
  int i, j, ii, jj, inum, jnum;
 
394
  int itype, jtype, iparam_i,iparam_j;
 
395
  int *ilist, *jlist, *numneigh, **firstneigh;
 
396
 
 
397
  int *type = atom->type;
 
398
  int nlocal = atom->nlocal;
 
399
  int newton_pair = force->newton_pair;
 
400
 
 
401
  double xtmp, ytmp, ztmp, ecoul, fpair;
 
402
  double qi, qj, selfion, r, rsq, delr[3];
 
403
  double zei, zej, zj, ci_jfi, dci_jfi, ci_fifj, dci_fifj; 
 
404
  double forcecoul, factor_coul;
 
405
 
 
406
  double **x = atom->x;
 
407
  double **f = atom->f;
 
408
  double *q = atom->q;
 
409
  double *special_coul = force->special_coul;
 
410
 
 
411
  ecoul = 0.0;
 
412
  selfion = fpair = 0.0;
 
413
  ci_jfi = ci_fifj = dci_jfi = dci_fifj = 0.0;
 
414
  forcecoul = 0.0;
 
415
 
 
416
  if (eflag || vflag) ev_setup(eflag,vflag);
 
417
  else evflag = vflag_fdotr = vflag_atom = 0;
 
418
 
 
419
  inum = list->inum;
 
420
  ilist = list->ilist;
 
421
  numneigh = list->numneigh;
 
422
  firstneigh = list->firstneigh;
 
423
 
 
424
  // Wolf sum
 
425
  
 
426
  if (kspacetype == 1) {
 
427
 
 
428
  for (ii = 0; ii < inum; ii++) {
 
429
    i = ilist[ii];
 
430
    xtmp = x[i][0];
 
431
    ytmp = x[i][1];
 
432
    ztmp = x[i][2];
 
433
    itype = map[type[i]];
 
434
    iparam_i = elem2param[itype];
 
435
    qi = q[i];
 
436
    zei = params[iparam_i].zeta;
 
437
 
 
438
    // self energy: ionization + wolf sum
 
439
 
 
440
    selfion = self(&params[iparam_i],qi);
 
441
 
 
442
    if (evflag) ev_tally(i,i,nlocal,0,0.0,selfion,0.0,0.0,0.0,0.0);
 
443
 
 
444
    // two-body interaction
 
445
 
 
446
    jlist = firstneigh[i];
 
447
    jnum = numneigh[i];
 
448
 
 
449
    for (jj = 0; jj < jnum; jj++) {
 
450
      j = jlist[jj];
 
451
      j &= NEIGHMASK;
 
452
 
 
453
      jtype = map[type[j]];
 
454
      iparam_j = elem2param[jtype];
 
455
      qj = q[j];
 
456
      zej = params[iparam_j].zeta;
 
457
      zj = params[iparam_j].zcore;
 
458
      factor_coul = special_coul[sbmask(j)];
 
459
 
 
460
      delr[0] = xtmp - x[j][0];
 
461
      delr[1] = ytmp - x[j][1];
 
462
      delr[2] = ztmp - x[j][2];
 
463
      rsq = delr[0]*delr[0] + delr[1]*delr[1] + delr[2]*delr[2];
 
464
 
 
465
      if (rsq > cut_coulsq) continue;
 
466
 
 
467
      r = sqrt(rsq);
 
468
 
 
469
      // Streitz-Mintmire Coulomb integrals
 
470
 
 
471
      coulomb_integral_wolf(zei, zej, r, ci_jfi, dci_jfi, ci_fifj, dci_fifj);
 
472
 
 
473
      // Wolf Sum
 
474
 
 
475
      wolf_sum(qi, qj, zj, r, ci_jfi, dci_jfi, ci_fifj, dci_fifj, 
 
476
               ecoul, forcecoul);
 
477
      
 
478
      // Forces
 
479
 
 
480
      fpair = -forcecoul / r;
 
481
 
 
482
      f[i][0] += delr[0]*fpair;
 
483
      f[i][1] += delr[1]*fpair;
 
484
      f[i][2] += delr[2]*fpair;
 
485
      if (newton_pair || j < nlocal) {
 
486
        f[j][0] -= delr[0]*fpair;
 
487
        f[j][1] -= delr[1]*fpair;
 
488
        f[j][2] -= delr[2]*fpair;
 
489
      }
 
490
 
 
491
      if (evflag) ev_tally(i,j,nlocal,newton_pair,
 
492
                           0.0,ecoul,fpair,delr[0],delr[1],delr[2]);
 
493
    }
 
494
 
 
495
  }
 
496
 
 
497
  // Ewald Sum
 
498
  
 
499
  } else if (kspacetype == 2) {
 
500
 
 
501
  for (ii = 0; ii < inum; ii++) {
 
502
    i = ilist[ii];
 
503
    xtmp = x[i][0];
 
504
    ytmp = x[i][1];
 
505
    ztmp = x[i][2];
 
506
    itype = map[type[i]];
 
507
    iparam_i = elem2param[itype];
 
508
    qi = q[i];
 
509
    zei = params[iparam_i].zeta;
 
510
 
 
511
    // self ionizition energy, only on i atom
 
512
 
 
513
    selfion = self(&params[iparam_i],qi);
 
514
 
 
515
    if (evflag) ev_tally(i,i,nlocal,0,0.0,selfion,0.0,0.0,0.0,0.0);
 
516
 
 
517
    // two-body interaction
 
518
 
 
519
    jlist = firstneigh[i];
 
520
    jnum = numneigh[i];
 
521
 
 
522
    for (jj = 0; jj < jnum; jj++) {
 
523
      j = jlist[jj];
 
524
      j &= NEIGHMASK;
 
525
      jtype = map[type[j]];
 
526
      iparam_j = elem2param[jtype];
 
527
      qj = q[j];
 
528
      zej = params[iparam_j].zeta;
 
529
      zj = params[iparam_j].zcore;
 
530
      factor_coul = special_coul[sbmask(j)];
 
531
 
 
532
      delr[0] = xtmp - x[j][0];
 
533
      delr[1] = ytmp - x[j][1];
 
534
      delr[2] = ztmp - x[j][2];
 
535
      rsq = delr[0]*delr[0] + delr[1]*delr[1] + delr[2]*delr[2];
 
536
 
 
537
      if (rsq > cut_coulsq) continue;
 
538
 
 
539
      r = sqrt(rsq);
 
540
 
 
541
      // Streitz-Mintmire Coulomb integrals
 
542
 
 
543
      coulomb_integral_ewald(zei, zej, r, ci_jfi, dci_jfi, ci_fifj, dci_fifj);
 
544
 
 
545
      // Ewald: real-space
 
546
      
 
547
      ewald_sum(qi, qj, zj, r, ci_jfi, dci_jfi, ci_fifj, dci_fifj, 
 
548
                      ecoul, forcecoul, factor_coul);
 
549
 
 
550
      // Forces
 
551
 
 
552
      fpair = -forcecoul / r;
 
553
 
 
554
      f[i][0] += delr[0]*fpair;
 
555
      f[i][1] += delr[1]*fpair;
 
556
      f[i][2] += delr[2]*fpair;
 
557
      if (newton_pair || j < nlocal) {
 
558
        f[j][0] -= delr[0]*fpair;
 
559
        f[j][1] -= delr[1]*fpair;
 
560
        f[j][2] -= delr[2]*fpair;
 
561
      }
 
562
 
 
563
      if (evflag) ev_tally(i,j,nlocal,newton_pair,
 
564
                           0.0,ecoul,fpair,delr[0],delr[1],delr[2]);
 
565
    }
 
566
  }
 
567
  }
 
568
 
 
569
  if (vflag_fdotr) virial_fdotr_compute();
 
570
}
 
571
 
 
572
/* ---------------------------------------------------------------------- */
 
573
 
 
574
double PairCoulStreitz::self(Param *param, double qi)
 
575
{
 
576
 double s1=param->chi, s2=param->eta;
 
577
 double qqrd2e = force->qqrd2e;
 
578
 
 
579
 if (kspacetype == 1) return 1.0*qi*(s1+qi*(0.50*s2 - qqrd2e*woself));
 
580
 
 
581
 if (kspacetype == 2) return 1.0*qi*(s1+qi*(0.50*s2));
 
582
 
 
583
 return 0.0;
 
584
}
 
585
 
 
586
/* ---------------------------------------------------------------------- */
 
587
 
 
588
void PairCoulStreitz::coulomb_integral_wolf(double zei, double zej, double r,
 
589
                  double &ci_jfi, double &dci_jfi, double &ci_fifj, 
 
590
                  double &dci_fifj)
 
591
{
 
592
  double rinv = 1.0/r;
 
593
  double rinv2 = rinv*rinv;
 
594
 
 
595
  double exp2zir = exp(-2.0*zei*r);
 
596
  double zei2 = zei*zei;
 
597
  double zei4 = zei2*zei2;
 
598
  double zei6 = zei2*zei4;
 
599
 
 
600
  double exp2zjr = exp(-2.0*zej*r);
 
601
  double zej2 = zej*zej;
 
602
  double zej4 = zej2*zej2;
 
603
  double zej6 = zej2*zej4;
 
604
 
 
605
  double sm1 = 11.0/8.0;
 
606
  double sm2 = 3.00/4.0;
 
607
  double sm3 = 1.00/6.0;
 
608
  double e1, e2, e3, e4;
 
609
 
 
610
  double rc = cut_coul;
 
611
  double rcinv = 1.0/rc;
 
612
  double rcinv2 = rcinv*rcinv;
 
613
  double exp2zirsh = exp(-2.0*zei*rc);
 
614
  double exp2zjrsh = exp(-2.0*zej*rc);
 
615
  double eshift, fshift;
 
616
 
 
617
  e1 = e2 = e3 = e4 = 0.0;
 
618
 
 
619
  eshift = -zei*exp2zirsh - rcinv*exp2zirsh;
 
620
  fshift = 2.0*zei2*exp2zirsh + rcinv2*exp2zirsh + 2.0*zei*rcinv*exp2zirsh;
 
621
 
 
622
  ci_jfi = -zei*exp2zir - rinv*exp2zir - eshift - (r-rc)*fshift;
 
623
  dci_jfi = 2.0*zei2*exp2zir + rinv2*exp2zir + 2.0*zei*rinv*exp2zir - fshift;
 
624
 
 
625
  if (zei == zej) {
 
626
 
 
627
    eshift = -exp2zirsh*(rcinv + zei*(sm1 + sm2*zei*rc + sm3*zei2*rc*rc));
 
628
    fshift =  exp2zirsh*(rcinv2 + 2.0*zei*rcinv + zei2*
 
629
                (2.0 + 7.0/6.0*zei*rc + 1.0/3.0*zei2*rc*rc));
 
630
 
 
631
    ci_fifj = -exp2zir*(rinv + zei*(sm1 + sm2*zei*r + sm3*zei2*r*r))
 
632
              - eshift - (r-rc)*fshift;
 
633
    dci_fifj = exp2zir*(rinv2 + 2.0*zei*rinv + zei2*
 
634
                (2.0 + 7.0/6.0*zei*r + 1.0/3.0*zei2*r*r)) - fshift;
 
635
 
 
636
  } else {
 
637
 
 
638
    e1 = zei*zej4/((zei+zej)*(zei+zej)*(zei-zej)*(zei-zej));
 
639
    e2 = zej*zei4/((zei+zej)*(zei+zej)*(zej-zei)*(zej-zei));
 
640
    e3 = (3.0*zei2*zej4-zej6) /
 
641
         ((zei+zej)*(zei+zej)*(zei+zej)*(zei-zej)*(zei-zej)*(zei-zej));
 
642
    e4 = (3.0*zej2*zei4-zei6) /
 
643
         ((zei+zej)*(zei+zej)*(zei+zej)*(zej-zei)*(zej-zei)*(zej-zei));
 
644
 
 
645
    eshift = -exp2zirsh*(e1+e3/rc) - exp2zjrsh*(e2+e4/rc);
 
646
    fshift = (exp2zirsh*(2.0*zei*(e1+e3/rc) + e3*rcinv2)
 
647
            + exp2zjrsh*(2.0*zej*(e2+e4/rc) + e4*rcinv2));
 
648
 
 
649
    ci_fifj = -exp2zir*(e1+e3/r) - exp2zjr*(e2+e4/r) 
 
650
              - eshift - (r-rc)*fshift;
 
651
    dci_fifj = (exp2zir*(2.0*zei*(e1+e3/r) + e3*rinv2) + 
 
652
                exp2zjr*(2.0*zej*(e2+e4/r) + e4*rinv2)) - fshift;
 
653
  }
 
654
}
 
655
 
 
656
/* ---------------------------------------------------------------------- */
 
657
 
 
658
void PairCoulStreitz::wolf_sum(double qi, double qj, double zj, double r, 
 
659
                double ci_jfi, double dci_jfi, double ci_fifj, 
 
660
                double dci_fifj, double &etmp, double &ftmp)
 
661
{
 
662
  double a = g_wolf;
 
663
  double rc = cut_coul;
 
664
  double qqrd2e = force->qqrd2e;
 
665
 
 
666
  double erfcr = erfc(a*r);
 
667
  double derfcr = exp(-a*a*r*r);
 
668
  double erfcrc = erfc(a*rc);
 
669
 
 
670
  double etmp1, etmp2, etmp3;
 
671
  double ftmp1, ftmp2, ftmp3;
 
672
 
 
673
  etmp = etmp1 = etmp2 = etmp3 = 0.0;
 
674
  ftmp = ftmp1 = ftmp2 = ftmp3 = 0.0;
 
675
 
 
676
  etmp1 = erfcr/r - erfcrc/rc;
 
677
  etmp2 = qi * zj * (ci_jfi - ci_fifj);
 
678
  etmp3 = qi * qj * 0.50 * (etmp1 + ci_fifj);
 
679
  
 
680
  ftmp1 = -erfcr/r/r - 2.0*a/MY_PIS*derfcr/r - dwoself;
 
681
  ftmp2 = qi * zj * (dci_jfi - dci_fifj);
 
682
  ftmp3 = qi * qj * 0.50 * (ftmp1 + dci_fifj);
 
683
 
 
684
  etmp = qqrd2e * (etmp2 + etmp3);
 
685
  ftmp = qqrd2e * (ftmp2 + ftmp3);
 
686
 
 
687
}
 
688
 
 
689
/* ---------------------------------------------------------------------- */
 
690
 
 
691
void PairCoulStreitz::coulomb_integral_ewald(double zei, double zej, double r,
 
692
                  double &ci_jfi, double &dci_jfi, double &ci_fifj, 
 
693
                  double &dci_fifj)
 
694
{
 
695
  double rinv = 1.0/r;
 
696
  double rinv2 = rinv*rinv;
 
697
 
 
698
  double exp2zir = exp(-2.0*zei*r);
 
699
  double zei2 = zei*zei;
 
700
  double zei4 = zei2*zei2;
 
701
  double zei6 = zei2*zei4;
 
702
 
 
703
  double exp2zjr = exp(-2.0*zej*r);
 
704
  double zej2 = zej*zej;
 
705
  double zej4 = zej2*zej2;
 
706
  double zej6 = zej2*zej4;
 
707
 
 
708
  double sm1 = 11.0/8.0;
 
709
  double sm2 = 3.00/4.0;
 
710
  double sm3 = 1.00/6.0;
 
711
 
 
712
  double e1, e2, e3, e4;
 
713
  e1 = e2 = e3 = e4 = 0.0;
 
714
 
 
715
  ci_jfi = -zei*exp2zir - rinv*exp2zir;
 
716
  dci_jfi = 2.0*zei2*exp2zir + rinv2*exp2zir + 2.0*zei*rinv*exp2zir;
 
717
 
 
718
  if (zei == zej) {
 
719
 
 
720
    ci_fifj = -exp2zir*(rinv + zei*(sm1 + sm2*zei*r + sm3*zei2*r*r));
 
721
    dci_fifj = exp2zir*(rinv2 + 2.0*zei*rinv + 
 
722
                zei2*(2.0 + 7.0/6.0*zei*r + 1.0/3.0*zei2*r*r));
 
723
 
 
724
  } else {
 
725
 
 
726
    e1 = zei*zej4/((zei+zej)*(zei+zej)*(zei-zej)*(zei-zej));
 
727
    e2 = zej*zei4/((zei+zej)*(zei+zej)*(zej-zei)*(zej-zei));
 
728
    e3 = (3.0*zei2*zej4-zej6) /
 
729
         ((zei+zej)*(zei+zej)*(zei+zej)*(zei-zej)*(zei-zej)*(zei-zej));
 
730
    e4 = (3.0*zej2*zei4-zei6) /
 
731
         ((zei+zej)*(zei+zej)*(zei+zej)*(zej-zei)*(zej-zei)*(zej-zei));
 
732
 
 
733
    ci_fifj = -exp2zir*(e1+e3/r) - exp2zjr*(e2+e4/r);
 
734
    dci_fifj = (exp2zir*(2.0*zei*(e1+e3/r) + e3*rinv2)
 
735
              + exp2zjr*(2.0*zej*(e2+e4/r) + e4*rinv2));
 
736
  }
 
737
 
 
738
}
 
739
 
 
740
/* ---------------------------------------------------------------------- */
 
741
 
 
742
void PairCoulStreitz::ewald_sum(double qi, double qj, double zj, double r, 
 
743
                double ci_jfi, double dci_jfi, double ci_fifj, 
 
744
                double dci_fifj, double &etmp, double &ftmp, double fac)
 
745
{
 
746
  double etmp1, etmp2, etmp3, etmp4;
 
747
  double ftmp1, ftmp2, ftmp3, ftmp4;
 
748
 
 
749
  double a = g_ewald;
 
750
  double qqrd2e = force->qqrd2e;
 
751
 
 
752
  double erfcr = erfc(a*r);
 
753
  double derfcr = exp(-a*a*r*r);
 
754
 
 
755
  etmp = ftmp = 0.0;
 
756
  etmp1 = etmp2 = etmp3 = etmp4 = 0.0;
 
757
  ftmp1 = ftmp2 = ftmp3 = ftmp4 = 0.0;
 
758
 
 
759
  etmp1 = qi * zj * (ci_jfi - ci_fifj);
 
760
  etmp2 = qi * qj * 0.50 * ci_fifj;
 
761
  etmp3 = qqrd2e * (etmp1 + etmp2);
 
762
  etmp4 = qqrd2e * 0.50*qi*qj/r;
 
763
  
 
764
  ftmp1 = qi * zj * (dci_jfi - dci_fifj);
 
765
  ftmp2 = qi * qj * 0.50 * dci_fifj;
 
766
  ftmp3 = qqrd2e * (ftmp1 + ftmp2);
 
767
  ftmp4 = etmp4 * (erfcr + 2.0/MY_PIS*a*r*derfcr);
 
768
 
 
769
  etmp = erfcr*etmp4;
 
770
 
 
771
  if (fac < 1.0) {
 
772
    etmp  -= (1.0-fac)*etmp4;
 
773
    ftmp4 -= (1.0-fac)*etmp4;
 
774
  }
 
775
 
 
776
  etmp += etmp3;
 
777
  ftmp = ftmp3 - ftmp4/r;
 
778
 
 
779
}
 
780
 
 
781
/* ----------------------------------------------------------------------
 
782
   memory usage of local atom-based arrays
 
783
------------------------------------------------------------------------- */
 
784
 
 
785
double PairCoulStreitz::memory_usage()
 
786
{
 
787
  double bytes = maxeatom * sizeof(double);
 
788
  bytes += maxvatom*6 * sizeof(double);
 
789
  bytes += nmax * sizeof(int);
 
790
  bytes += MAXNEIGH * nmax * sizeof(int);
 
791
  return bytes;
 
792
}
 
793
 
 
794
/* ---------------------------------------------------------------------- */
 
795
 
 
796
void *PairCoulStreitz::extract(const char *str, int &dim)
 
797
{
 
798
  if (strcmp(str,"cut_coul") == 0) {
 
799
    dim = 0;
 
800
    return (void *) &cut_coul;
 
801
  }
 
802
  if (strcmp(str,"scale") == 0) {
 
803
    dim = 2;
 
804
    return (void *) scale;
 
805
  }
 
806
  if (strcmp(str,"chi") == 0 && qeq_x) {
 
807
    dim = 1;
 
808
    for (int i = 1; i <= atom->ntypes; i++)
 
809
      if (map[i] >= 0) qeq_x[i] = params[map[i]].chi;
 
810
      else qeq_x[i] = 0.0;
 
811
    return (void *) qeq_x;
 
812
  }
 
813
  if (strcmp(str,"eta") == 0 && qeq_j) {
 
814
    dim = 1;
 
815
    for (int i = 1; i <= atom->ntypes; i++)
 
816
      if (map[i] >= 0) qeq_j[i] = params[map[i]].eta;
 
817
      else qeq_j[i] = 0.0;
 
818
    return (void *) qeq_j;
 
819
  }
 
820
  if (strcmp(str,"gamma") == 0 && qeq_g) {
 
821
    dim = 1;
 
822
    for (int i = 1; i <= atom->ntypes; i++)
 
823
      if (map[i] >= 0) qeq_g[i] = params[map[i]].gamma;
 
824
      else qeq_g[i] = 0.0;
 
825
    return (void *) qeq_g;
 
826
  }
 
827
  if (strcmp(str,"zeta") == 0 && qeq_z) {
 
828
    dim = 1;
 
829
    for (int i = 1; i <= atom->ntypes; i++)
 
830
      if (map[i] >= 0) qeq_z[i] = params[map[i]].zeta;
 
831
      else qeq_z[i] = 0.0;
 
832
    return (void *) qeq_z;
 
833
  }
 
834
  if (strcmp(str,"zcore") == 0 && qeq_c) {
 
835
    dim = 1;
 
836
    for (int i = 1; i <= atom->ntypes; i++)
 
837
      if (map[i] >= 0) qeq_c[i] = params[map[i]].zcore;
 
838
      else qeq_c[i] = 0.0;
 
839
    return (void *) qeq_c;
 
840
  }
 
841
  if (strcmp(str,"kspacetype") == 0) {
 
842
    dim = 0;
 
843
    return (void *) &kspacetype;
 
844
  }
 
845
  if (strcmp(str,"alpha") == 0) {
 
846
    dim = 0;
 
847
    if (kspacetype == 1) return (void *) &g_wolf;
 
848
    if (kspacetype == 2) return (void *) &g_ewald;
 
849
  }
 
850
  return NULL;
 
851
}