jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5099 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2010-10-22 21:25:05 +00:00
parent 1d39d319af
commit 8372b98a29
14 changed files with 82 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
src/KSPACE/ewald.cpp
View File

@ -18,6 +18,7 @@
#include "mpi.h"
#include "stdlib.h"
#include "stdio.h"
#include "string.h"
#include "math.h"
#include "ewald.h"
#include "atom.h"
@ -97,10 +98,12 @@ void Ewald::init()
// extract short-range Coulombic cutoff from pair style
qqrd2e = force->qqrd2e;
scale = 1.0;
if (force->pair == NULL)
error->all("KSpace style is incompatible with Pair style");
double *p_cutoff = (double *) force->pair->extract("cut_coul");
int itmp;
double *p_cutoff = (double *) force->pair->extract("cut_coul",itmp);
if (p_cutoff == NULL)
error->all("KSpace style is incompatible with Pair style");
double cutoff = *p_cutoff;
@ -270,9 +273,9 @@ void Ewald::compute(int eflag, int vflag)
// convert E-field to force
for (i = 0; i < nlocal; i++) {
f[i][0] += qqrd2e*q[i]*ek[i][0];
f[i][1] += qqrd2e*q[i]*ek[i][1];
f[i][2] += qqrd2e*q[i]*ek[i][2];
f[i][0] += qqrd2e*scale * q[i]*ek[i][0];
f[i][1] += qqrd2e*scale * q[i]*ek[i][1];
f[i][2] += qqrd2e*scale * q[i]*ek[i][2];
}
// energy if requested
@ -284,7 +287,7 @@ void Ewald::compute(int eflag, int vflag)
PI = 4.0*atan(1.0);
energy -= g_ewald*qsqsum/1.772453851 +
0.5*PI*qsum*qsum / (g_ewald*g_ewald*volume);
energy *= qqrd2e;
energy *= qqrd2e*scale;
}
// virial if requested
@ -295,11 +298,10 @@ void Ewald::compute(int eflag, int vflag)
uk = ug[k] * (sfacrl_all[k]*sfacrl_all[k] + sfacim_all[k]*sfacim_all[k]);
for (n = 0; n < 6; n++) virial[n] += uk*vg[k][n];
}
for (n = 0; n < 6; n++) virial[n] *= qqrd2e;
for (n = 0; n < 6; n++) virial[n] *= qqrd2e*scale;
}
if (slabflag) slabcorr(eflag);
}
/* ---------------------------------------------------------------------- */
@ -820,14 +822,14 @@ void Ewald::slabcorr(int eflag)
double e_slabcorr = 2.0*PI*dipole_all*dipole_all/volume;
if (eflag) energy += qqrd2e*e_slabcorr;
if (eflag) energy += qqrd2e*scale * e_slabcorr;
// add on force corrections
double ffact = -4.0*PI*dipole_all/volume;
double **f = atom->f;
for (int i = 0; i < nlocal; i++) f[i][2] += qqrd2e*q[i]*ffact;
for (int i = 0; i < nlocal; i++) f[i][2] += qqrd2e*scale * q[i]*ffact;
}
/* ----------------------------------------------------------------------

3
src/KSPACE/pair_born_coul_long.cpp
View File

@ -497,8 +497,9 @@ double PairBornCoulLong::single(int i, int j, int itype, int jtype,
/* ---------------------------------------------------------------------- */
void *PairBornCoulLong::extract(char *str)
void *PairBornCoulLong::extract(char *str, int &dim)
{
dim = 0;
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
return NULL;
}

2
src/KSPACE/pair_born_coul_long.h
View File

@ -38,7 +38,7 @@ class PairBornCoulLong : public Pair {
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
double single(int, int, int, int, double, double, double, double &);
void *extract(char *);
void *extract(char *, int &);
private:
double cut_lj_global;

4
src/KSPACE/pair_buck_coul_long.cpp
View File

@ -437,9 +437,9 @@ double PairBuckCoulLong::single(int i, int j, int itype, int jtype,
/* ---------------------------------------------------------------------- */
void *PairBuckCoulLong::extract(char *str)
void *PairBuckCoulLong::extract(char *str, int &dim)
{
dim = 0;
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
return NULL;
}

2
src/KSPACE/pair_buck_coul_long.h
View File

@ -38,7 +38,7 @@ class PairBuckCoulLong : public Pair {
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
double single(int, int, int, int, double, double, double, double &);
void *extract(char *);
void *extract(char *, int &);
private:
double cut_lj_global;

37
src/KSPACE/pair_coul_long.cpp
View File

@ -59,6 +59,8 @@ PairCoulLong::~PairCoulLong()
if (allocated) {
memory->destroy_2d_int_array(setflag);
memory->destroy_2d_double_array(cutsq);
memory->destroy_2d_double_array(scale);
}
if (ftable) free_tables();
}
@ -67,7 +69,7 @@ PairCoulLong::~PairCoulLong()
void PairCoulLong::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itable;
int i,j,ii,jj,inum,jnum,itable,itype,jtype;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
double fraction,table;
double r,r2inv,forcecoul,factor_coul;
@ -82,6 +84,7 @@ void PairCoulLong::compute(int eflag, int vflag)
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
int *type = atom->type;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
double *special_coul = force->special_coul;
@ -101,6 +104,7 @@ void PairCoulLong::compute(int eflag, int vflag)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
@ -117,16 +121,17 @@ void PairCoulLong::compute(int eflag, int vflag)
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cut_coulsq) {
r2inv = 1.0/rsq;
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrt(rsq);
jtype = type[j];
if (rsq < cut_coulsq) {
r2inv = 1.0/rsq;
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
prefactor = qqrd2e * qtmp*q[j]/r;
prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
} else {
@ -136,10 +141,10 @@ void PairCoulLong::compute(int eflag, int vflag)
itable >>= ncoulshiftbits;
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = qtmp*q[j] * table;
forcecoul = scale[itype][jtype] * qtmp*q[j] * table;
if (factor_coul < 1.0) {
table = ctable[itable] + fraction*dctable[itable];
prefactor = qtmp*q[j] * table;
prefactor = scale[itype][jtype] * qtmp*q[j] * table;
forcecoul -= (1.0-factor_coul)*prefactor;
}
}
@ -160,7 +165,7 @@ void PairCoulLong::compute(int eflag, int vflag)
ecoul = prefactor*erfc;
else {
table = etable[itable] + fraction*detable[itable];
ecoul = qtmp*q[j] * table;
ecoul = scale[itype][jtype] * qtmp*q[j] * table;
}
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
}
@ -189,6 +194,8 @@ void PairCoulLong::allocate()
setflag[i][j] = 0;
cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq");
scale = memory->create_2d_double_array(n+1,n+1,"pair:scale");
}
/* ----------------------------------------------------------------------
@ -218,6 +225,7 @@ void PairCoulLong::coeff(int narg, char **arg)
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
scale[i][j] = 1.0;
setflag[i][j] = 1;
count++;
}
@ -556,8 +564,15 @@ double PairCoulLong::single(int i, int j, int itype, int jtype,
/* ---------------------------------------------------------------------- */
void *PairCoulLong::extract(char *str)
void *PairCoulLong::extract(char *str, int &dim)
{
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
if (strcmp(str,"cut_coul") == 0) {
dim = 0;
return (void *) &cut_coul;
}
if (strcmp(str,"scale") == 0) {
dim = 2;
return (void *) scale;
}
return NULL;
}

3
src/KSPACE/pair_coul_long.h
View File

@ -38,12 +38,13 @@ class PairCoulLong : public Pair {
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
double single(int, int, int, int, double, double, double, double &);
void *extract(char *);
void *extract(char *, int &);
private:
double cut_coul,cut_coulsq;
double *cut_respa;
double g_ewald;
double **scale;
double tabinnersq;
double *rtable,*drtable,*ftable,*dftable,*ctable,*dctable;

16
src/KSPACE/pair_lj_charmm_coul_long.cpp
View File

@ -1201,13 +1201,17 @@ double PairLJCharmmCoulLong::single(int i, int j, int itype, int jtype,
/* ---------------------------------------------------------------------- */
void *PairLJCharmmCoulLong::extract(char *str)
void *PairLJCharmmCoulLong::extract(char *str, int &dim)
{
dim = 2;
if (strcmp(str,"lj14_1") == 0) return (void *) lj14_1;
else if (strcmp(str,"lj14_2") == 0) return (void *) lj14_2;
else if (strcmp(str,"lj14_3") == 0) return (void *) lj14_3;
else if (strcmp(str,"lj14_4") == 0) return (void *) lj14_4;
else if (strcmp(str,"implicit") == 0) return (void *) &implicit;
else if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
if (strcmp(str,"lj14_2") == 0) return (void *) lj14_2;
if (strcmp(str,"lj14_3") == 0) return (void *) lj14_3;
if (strcmp(str,"lj14_4") == 0) return (void *) lj14_4;
dim = 0;
if (strcmp(str,"implicit") == 0) return (void *) &implicit;
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
return NULL;
}

2
src/KSPACE/pair_lj_charmm_coul_long.h
View File

@ -43,7 +43,7 @@ class PairLJCharmmCoulLong : public Pair {
void compute_inner();
void compute_middle();
void compute_outer(int, int);
void *extract(char *);
void *extract(char *, int &);
protected:
int implicit;

3
src/KSPACE/pair_lj_cut_coul_long.cpp
View File

@ -1144,8 +1144,9 @@ double PairLJCutCoulLong::single(int i, int j, int itype, int jtype,
/* ---------------------------------------------------------------------- */
void *PairLJCutCoulLong::extract(char *str)
void *PairLJCutCoulLong::extract(char *str, int &dim)
{
dim = 0;
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
return NULL;
}

2
src/KSPACE/pair_lj_cut_coul_long.h
View File

@ -43,7 +43,7 @@ class PairLJCutCoulLong : public Pair {
void compute_inner();
void compute_middle();
void compute_outer(int, int);
void *extract(char *);
void *extract(char *, int &);
protected:
double cut_lj_global;

13
src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp
View File

@ -496,13 +496,14 @@ void PairLJCutCoulLongTIP4P::find_M(int i, int &iH1, int &iH2, double *xM)
/* ---------------------------------------------------------------------- */
void *PairLJCutCoulLongTIP4P::extract(char *str)
void *PairLJCutCoulLongTIP4P::extract(char *str, int &dim)
{
dim = 0;
if (strcmp(str,"qdist") == 0) return (void *) &qdist;
else if (strcmp(str,"typeO") == 0) return (void *) &typeO;
else if (strcmp(str,"typeH") == 0) return (void *) &typeH;
else if (strcmp(str,"typeA") == 0) return (void *) &typeA;
else if (strcmp(str,"typeB") == 0) return (void *) &typeB;
else if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
if (strcmp(str,"typeO") == 0) return (void *) &typeO;
if (strcmp(str,"typeH") == 0) return (void *) &typeH;
if (strcmp(str,"typeA") == 0) return (void *) &typeA;
if (strcmp(str,"typeB") == 0) return (void *) &typeB;
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
return NULL;
}

2
src/KSPACE/pair_lj_cut_coul_long_tip4p.h
View File

@ -32,7 +32,7 @@ class PairLJCutCoulLongTIP4P : public PairLJCutCoulLong {
void init_style();
void write_restart_settings(FILE *fp);
void read_restart_settings(FILE *fp);
void *extract(char *);
void *extract(char *, int &);
private:
int typeH,typeO; // atom types of TIP4P water H and O atoms

29
src/KSPACE/pppm.cpp
View File

@ -132,10 +132,12 @@ void PPPM::init()
// extract short-range Coulombic cutoff from pair style
qqrd2e = force->qqrd2e;
scale = 1.0;
if (force->pair == NULL)
error->all("KSpace style is incompatible with Pair style");
double *p_cutoff = (double *) force->pair->extract("cut_coul");
int itmp;
double *p_cutoff = (double *) force->pair->extract("cut_coul",itmp);
if (p_cutoff == NULL)
error->all("KSpace style is incompatible with Pair style");
cutoff = *p_cutoff;
@ -147,11 +149,11 @@ void PPPM::init()
if (strcmp(force->kspace_style,"pppm/tip4p") == 0) {
if (force->pair == NULL)
error->all("KSpace style is incompatible with Pair style");
double *p_qdist = (double *) force->pair->extract("qdist");
int *p_typeO = (int *) force->pair->extract("typeO");
int *p_typeH = (int *) force->pair->extract("typeH");
int *p_typeA = (int *) force->pair->extract("typeA");
int *p_typeB = (int *) force->pair->extract("typeB");
double *p_qdist = (double *) force->pair->extract("qdist",itmp);
int *p_typeO = (int *) force->pair->extract("typeO",itmp);
int *p_typeH = (int *) force->pair->extract("typeH",itmp);
int *p_typeA = (int *) force->pair->extract("typeA",itmp);
int *p_typeB = (int *) force->pair->extract("typeB",itmp);
if (!p_qdist || !p_typeO || !p_typeH || !p_typeA || !p_typeB)
error->all("KSpace style is incompatible with Pair style");
qdist = *p_qdist;
@ -694,7 +696,7 @@ void PPPM::compute(int eflag, int vflag)
energy *= 0.5*volume;
energy -= g_ewald*qsqsum/1.772453851 +
0.5*PI*qsum*qsum / (g_ewald*g_ewald*volume);
energy *= qqrd2e;
energy *= qqrd2e*scale;
}
// sum virial across procs
@ -702,7 +704,7 @@ void PPPM::compute(int eflag, int vflag)
if (vflag) {
double virial_all[6];
MPI_Allreduce(virial,virial_all,6,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < 6; i++) virial[i] = 0.5*qqrd2e*volume*virial_all[i];
for (i = 0; i < 6; i++) virial[i] = 0.5*qqrd2e*scale*volume*virial_all[i];
}
// 2d slab correction
@ -1680,7 +1682,6 @@ void PPPM::fieldforce()
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
nx = part2grid[i][0];
ny = part2grid[i][1];
nz = part2grid[i][2];
@ -1709,9 +1710,9 @@ void PPPM::fieldforce()
// convert E-field to force
f[i][0] += qqrd2e*q[i]*ek[0];
f[i][1] += qqrd2e*q[i]*ek[1];
f[i][2] += qqrd2e*q[i]*ek[2];
f[i][0] += qqrd2e*scale * q[i]*ek[0];
f[i][1] += qqrd2e*scale * q[i]*ek[1];
f[i][2] += qqrd2e*scale * q[i]*ek[2];
}
}
@ -1854,14 +1855,14 @@ void PPPM::slabcorr(int eflag)
double e_slabcorr = 2.0*PI*dipole_all*dipole_all/volume;
if (eflag) energy += qqrd2e*e_slabcorr;
if (eflag) energy += qqrd2e*scale * e_slabcorr;
// add on force corrections
double ffact = -4.0*PI*dipole_all/volume;
double **f = atom->f;
for (int i = 0; i < nlocal; i++) f[i][2] += qqrd2e*q[i]*ffact;
for (int i = 0; i < nlocal; i++) f[i][2] += qqrd2e*scale * q[i]*ffact;
}
/* ----------------------------------------------------------------------