1
/* ----------------------------------------------------------------------
2
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3
http://lammps.sandia.gov, Sandia National Laboratories
4
Steve Plimpton, sjplimp@sandia.gov
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.
11
See the README file in the top-level LAMMPS directory.
12
------------------------------------------------------------------------- */
14
/* ----------------------------------------------------------------------
15
Contributing author: Ray Shan (Sandia)
16
------------------------------------------------------------------------- */
22
#include "pair_coul_streitz.h"
28
#include "neigh_list.h"
29
#include "neigh_request.h"
32
#include "math_const.h"
36
using namespace LAMMPS_NS;
37
using namespace MathConst;
44
/* ---------------------------------------------------------------------- */
46
PairCoulStreitz::PairCoulStreitz(LAMMPS *lmp) : Pair(lmp)
61
/* ----------------------------------------------------------------------
62
check if allocated, since class can be destructed when incomplete
63
------------------------------------------------------------------------- */
65
PairCoulStreitz::~PairCoulStreitz()
68
for (int i = 0; i < nelements; i++) delete [] elements[i];
71
memory->sfree(params);
72
memory->destroy(elem2param);
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);
87
/* ---------------------------------------------------------------------- */
89
void PairCoulStreitz::allocate()
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");
106
/* ----------------------------------------------------------------------
108
------------------------------------------------------------------------- */
110
void PairCoulStreitz::settings(int narg, char **arg)
112
if (narg < 2) error->all(FLERR,"Illegal pair_style command");
114
cut_coul = force->numeric(FLERR,arg[0]);
116
if (strcmp(arg[1],"wolf") == 0){
118
g_wolf = force->numeric(FLERR,arg[2]);
119
} else if (strcmp(arg[1],"ewald") == 0){
120
ewaldflag = pppmflag = 1;
123
error->all(FLERR,"Illegal pair_style command");
127
/* ----------------------------------------------------------------------
128
set coeffs for one or more type pairs
129
------------------------------------------------------------------------- */
131
void PairCoulStreitz::coeff(int narg, char **arg)
135
if (!allocated) allocate();
137
if (narg != 3 + atom->ntypes)
138
error->all(FLERR,"Incorrect args for pair coefficients");
140
// insure I,J args are * *
142
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
143
error->all(FLERR,"Incorrect args for pair coefficients");
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
151
for (i = 0; i < nelements; i++) delete [] elements[i];
154
elements = new char*[atom->ntypes];
155
for (i = 0; i < atom->ntypes; i++) elements[i] = NULL;
158
for (i = 3; i < narg; i++) {
159
if (strcmp(arg[i],"NULL") == 0) {
163
for (j = 0; j < nelements; j++)
164
if (strcmp(arg[i],elements[j]) == 0) break;
166
if (j == nelements) {
167
n = strlen(arg[i]) + 1;
168
elements[j] = new char[n];
169
strcpy(elements[j],arg[i]);
174
// read potential file and initialize potential parameters
180
// clear setflag since coeff() called once with I,J = * *
182
for (int i = 1; i <= n; i++)
183
for (int j = i; j <= n; j++)
186
// set setflag i,j for type pairs where both are mapped to elements
189
for (int i = 1; i <= n; i++)
190
for (int j = i; j <= n; j++)
191
if (map[i] >= 0 && map[j] >= 0) {
197
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
200
/* ----------------------------------------------------------------------
201
init specific to this pair style
202
------------------------------------------------------------------------- */
204
void PairCoulStreitz::init_style()
207
error->all(FLERR,"Pair style coul/sm requires atom attribute q");
209
//neighbor->request(this);
210
int irequest = neighbor->request(this,instance_me);
211
neighbor->requests[irequest]->half = 0;
212
neighbor->requests[irequest]->full = 1;
214
cut_coulsq = cut_coul * cut_coul;
216
// insure use of KSpace long-range solver when ewald specified, set g_ewald
219
if (force->kspace == NULL)
220
error->all(FLERR,"Pair style requires KSpace style ewald");
221
g_ewald = force->kspace->g_ewald;
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;
231
/* ----------------------------------------------------------------------
232
init for one type pair i,j and corresponding j,i
233
------------------------------------------------------------------------- */
235
double PairCoulStreitz::init_one(int i, int j)
237
scale[j][i] = scale[i][j];
239
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
243
/* ---------------------------------------------------------------------- */
245
void PairCoulStreitz::read_file(char *file)
247
int params_per_line = 6;
248
char **words = new char*[params_per_line+1];
250
memory->sfree(params);
255
// open file on proc 0
259
fp = fopen(file,"r");
262
sprintf(str,"Cannot open coul/Streitz potential file %s",file);
263
error->one(FLERR,str);
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
270
int n,nwords,ielement;
271
char line[MAXLINE],*ptr;
276
ptr = fgets(line,MAXLINE,fp);
280
} else n = strlen(line) + 1;
282
MPI_Bcast(&eof,1,MPI_INT,0,world);
284
MPI_Bcast(&n,1,MPI_INT,0,world);
285
MPI_Bcast(line,n,MPI_CHAR,0,world);
287
// strip comment, skip line if blank
289
if ((ptr = strchr(line,'#'))) *ptr = '\0';
290
nwords = atom->count_words(line);
291
if (nwords == 0) continue;
293
// concatenate additional lines until have params_per_line words
295
while (nwords < params_per_line) {
298
ptr = fgets(&line[n],MAXLINE-n,fp);
302
} else n = strlen(line) + 1;
304
MPI_Bcast(&eof,1,MPI_INT,0,world);
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);
312
if (nwords != params_per_line)
313
error->all(FLERR,"Incorrect format in coul/Streitz potential file");
315
// words = ptrs to all words in line
318
words[nwords++] = strtok(line," \t\n\r\f");
319
while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
321
// ielement = 1st args
323
for (ielement = 0; ielement < nelements; ielement++)
324
if (strcmp(words[0],elements[ielement]) == 0) break;
325
if (ielement == nelements) continue;
327
// load up parameter settings and error check their values
329
if (nparams == maxparam) {
331
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
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]);
342
// parameter sanity check
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");
354
/* ---------------------------------------------------------------------- */
356
void PairCoulStreitz::setup()
362
memory->destroy(elem2param);
363
memory->create(elem2param,nelements,"pair:elem2param");
365
for (i = 0; i < nelements; i++) {
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");
373
if (n < 0) error->all(FLERR,"Potential file is missing an entry");
377
// Wolf sum self energy
379
if (kspacetype == 1) {
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);
389
/* ---------------------------------------------------------------------- */
391
void PairCoulStreitz::compute(int eflag, int vflag)
393
int i, j, ii, jj, inum, jnum;
394
int itype, jtype, iparam_i,iparam_j;
395
int *ilist, *jlist, *numneigh, **firstneigh;
397
int *type = atom->type;
398
int nlocal = atom->nlocal;
399
int newton_pair = force->newton_pair;
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;
406
double **x = atom->x;
407
double **f = atom->f;
409
double *special_coul = force->special_coul;
412
selfion = fpair = 0.0;
413
ci_jfi = ci_fifj = dci_jfi = dci_fifj = 0.0;
416
if (eflag || vflag) ev_setup(eflag,vflag);
417
else evflag = vflag_fdotr = vflag_atom = 0;
421
numneigh = list->numneigh;
422
firstneigh = list->firstneigh;
426
if (kspacetype == 1) {
428
for (ii = 0; ii < inum; ii++) {
433
itype = map[type[i]];
434
iparam_i = elem2param[itype];
436
zei = params[iparam_i].zeta;
438
// self energy: ionization + wolf sum
440
selfion = self(¶ms[iparam_i],qi);
442
if (evflag) ev_tally(i,i,nlocal,0,0.0,selfion,0.0,0.0,0.0,0.0);
444
// two-body interaction
446
jlist = firstneigh[i];
449
for (jj = 0; jj < jnum; jj++) {
453
jtype = map[type[j]];
454
iparam_j = elem2param[jtype];
456
zej = params[iparam_j].zeta;
457
zj = params[iparam_j].zcore;
458
factor_coul = special_coul[sbmask(j)];
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];
465
if (rsq > cut_coulsq) continue;
469
// Streitz-Mintmire Coulomb integrals
471
coulomb_integral_wolf(zei, zej, r, ci_jfi, dci_jfi, ci_fifj, dci_fifj);
475
wolf_sum(qi, qj, zj, r, ci_jfi, dci_jfi, ci_fifj, dci_fifj,
480
fpair = -forcecoul / r;
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;
491
if (evflag) ev_tally(i,j,nlocal,newton_pair,
492
0.0,ecoul,fpair,delr[0],delr[1],delr[2]);
499
} else if (kspacetype == 2) {
501
for (ii = 0; ii < inum; ii++) {
506
itype = map[type[i]];
507
iparam_i = elem2param[itype];
509
zei = params[iparam_i].zeta;
511
// self ionizition energy, only on i atom
513
selfion = self(¶ms[iparam_i],qi);
515
if (evflag) ev_tally(i,i,nlocal,0,0.0,selfion,0.0,0.0,0.0,0.0);
517
// two-body interaction
519
jlist = firstneigh[i];
522
for (jj = 0; jj < jnum; jj++) {
525
jtype = map[type[j]];
526
iparam_j = elem2param[jtype];
528
zej = params[iparam_j].zeta;
529
zj = params[iparam_j].zcore;
530
factor_coul = special_coul[sbmask(j)];
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];
537
if (rsq > cut_coulsq) continue;
541
// Streitz-Mintmire Coulomb integrals
543
coulomb_integral_ewald(zei, zej, r, ci_jfi, dci_jfi, ci_fifj, dci_fifj);
547
ewald_sum(qi, qj, zj, r, ci_jfi, dci_jfi, ci_fifj, dci_fifj,
548
ecoul, forcecoul, factor_coul);
552
fpair = -forcecoul / r;
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;
563
if (evflag) ev_tally(i,j,nlocal,newton_pair,
564
0.0,ecoul,fpair,delr[0],delr[1],delr[2]);
569
if (vflag_fdotr) virial_fdotr_compute();
572
/* ---------------------------------------------------------------------- */
574
double PairCoulStreitz::self(Param *param, double qi)
576
double s1=param->chi, s2=param->eta;
577
double qqrd2e = force->qqrd2e;
579
if (kspacetype == 1) return 1.0*qi*(s1+qi*(0.50*s2 - qqrd2e*woself));
581
if (kspacetype == 2) return 1.0*qi*(s1+qi*(0.50*s2));
586
/* ---------------------------------------------------------------------- */
588
void PairCoulStreitz::coulomb_integral_wolf(double zei, double zej, double r,
589
double &ci_jfi, double &dci_jfi, double &ci_fifj,
593
double rinv2 = rinv*rinv;
595
double exp2zir = exp(-2.0*zei*r);
596
double zei2 = zei*zei;
597
double zei4 = zei2*zei2;
598
double zei6 = zei2*zei4;
600
double exp2zjr = exp(-2.0*zej*r);
601
double zej2 = zej*zej;
602
double zej4 = zej2*zej2;
603
double zej6 = zej2*zej4;
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;
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;
617
e1 = e2 = e3 = e4 = 0.0;
619
eshift = -zei*exp2zirsh - rcinv*exp2zirsh;
620
fshift = 2.0*zei2*exp2zirsh + rcinv2*exp2zirsh + 2.0*zei*rcinv*exp2zirsh;
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;
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));
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;
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));
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));
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;
656
/* ---------------------------------------------------------------------- */
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)
663
double rc = cut_coul;
664
double qqrd2e = force->qqrd2e;
666
double erfcr = erfc(a*r);
667
double derfcr = exp(-a*a*r*r);
668
double erfcrc = erfc(a*rc);
670
double etmp1, etmp2, etmp3;
671
double ftmp1, ftmp2, ftmp3;
673
etmp = etmp1 = etmp2 = etmp3 = 0.0;
674
ftmp = ftmp1 = ftmp2 = ftmp3 = 0.0;
676
etmp1 = erfcr/r - erfcrc/rc;
677
etmp2 = qi * zj * (ci_jfi - ci_fifj);
678
etmp3 = qi * qj * 0.50 * (etmp1 + ci_fifj);
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);
684
etmp = qqrd2e * (etmp2 + etmp3);
685
ftmp = qqrd2e * (ftmp2 + ftmp3);
689
/* ---------------------------------------------------------------------- */
691
void PairCoulStreitz::coulomb_integral_ewald(double zei, double zej, double r,
692
double &ci_jfi, double &dci_jfi, double &ci_fifj,
696
double rinv2 = rinv*rinv;
698
double exp2zir = exp(-2.0*zei*r);
699
double zei2 = zei*zei;
700
double zei4 = zei2*zei2;
701
double zei6 = zei2*zei4;
703
double exp2zjr = exp(-2.0*zej*r);
704
double zej2 = zej*zej;
705
double zej4 = zej2*zej2;
706
double zej6 = zej2*zej4;
708
double sm1 = 11.0/8.0;
709
double sm2 = 3.00/4.0;
710
double sm3 = 1.00/6.0;
712
double e1, e2, e3, e4;
713
e1 = e2 = e3 = e4 = 0.0;
715
ci_jfi = -zei*exp2zir - rinv*exp2zir;
716
dci_jfi = 2.0*zei2*exp2zir + rinv2*exp2zir + 2.0*zei*rinv*exp2zir;
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));
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));
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));
740
/* ---------------------------------------------------------------------- */
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)
746
double etmp1, etmp2, etmp3, etmp4;
747
double ftmp1, ftmp2, ftmp3, ftmp4;
750
double qqrd2e = force->qqrd2e;
752
double erfcr = erfc(a*r);
753
double derfcr = exp(-a*a*r*r);
756
etmp1 = etmp2 = etmp3 = etmp4 = 0.0;
757
ftmp1 = ftmp2 = ftmp3 = ftmp4 = 0.0;
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;
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);
772
etmp -= (1.0-fac)*etmp4;
773
ftmp4 -= (1.0-fac)*etmp4;
777
ftmp = ftmp3 - ftmp4/r;
781
/* ----------------------------------------------------------------------
782
memory usage of local atom-based arrays
783
------------------------------------------------------------------------- */
785
double PairCoulStreitz::memory_usage()
787
double bytes = maxeatom * sizeof(double);
788
bytes += maxvatom*6 * sizeof(double);
789
bytes += nmax * sizeof(int);
790
bytes += MAXNEIGH * nmax * sizeof(int);
794
/* ---------------------------------------------------------------------- */
796
void *PairCoulStreitz::extract(const char *str, int &dim)
798
if (strcmp(str,"cut_coul") == 0) {
800
return (void *) &cut_coul;
802
if (strcmp(str,"scale") == 0) {
804
return (void *) scale;
806
if (strcmp(str,"chi") == 0 && qeq_x) {
808
for (int i = 1; i <= atom->ntypes; i++)
809
if (map[i] >= 0) qeq_x[i] = params[map[i]].chi;
811
return (void *) qeq_x;
813
if (strcmp(str,"eta") == 0 && qeq_j) {
815
for (int i = 1; i <= atom->ntypes; i++)
816
if (map[i] >= 0) qeq_j[i] = params[map[i]].eta;
818
return (void *) qeq_j;
820
if (strcmp(str,"gamma") == 0 && qeq_g) {
822
for (int i = 1; i <= atom->ntypes; i++)
823
if (map[i] >= 0) qeq_g[i] = params[map[i]].gamma;
825
return (void *) qeq_g;
827
if (strcmp(str,"zeta") == 0 && qeq_z) {
829
for (int i = 1; i <= atom->ntypes; i++)
830
if (map[i] >= 0) qeq_z[i] = params[map[i]].zeta;
832
return (void *) qeq_z;
834
if (strcmp(str,"zcore") == 0 && qeq_c) {
836
for (int i = 1; i <= atom->ntypes; i++)
837
if (map[i] >= 0) qeq_c[i] = params[map[i]].zcore;
839
return (void *) qeq_c;
841
if (strcmp(str,"kspacetype") == 0) {
843
return (void *) &kspacetype;
845
if (strcmp(str,"alpha") == 0) {
847
if (kspacetype == 1) return (void *) &g_wolf;
848
if (kspacetype == 2) return (void *) &g_ewald;