ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
sjplimp
2013-01-02 16:53:23 +00:00
parent 8106d6ed04
commit 0a4cdacb67
4 changed files with 266 additions and 114 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/USER-EFF/atom_vec_electron.cpp
View File

@ -49,6 +49,8 @@ AtomVecElectron::AtomVecElectron(LAMMPS *lmp, int narg, char **arg) :
size_data_vel = 5;
xcol_data = 6;
atom->ecp_flag = 0;
atom->electron_flag = 1;
atom->q_flag = atom->spin_flag = atom->eradius_flag =
atom->ervel_flag = atom->erforce_flag = 1;
@ -766,6 +768,8 @@ void AtomVecElectron::data_atom(double *coord, tagint imagetmp, char **values)
q[nlocal] = atof(values[2]);
spin[nlocal] = atoi(values[3]);
if (spin[nlocal] == 3) atom->ecp_flag = 1;
eradius[nlocal] = atof(values[4]);
x[nlocal][0] = coord[0];

335
src/USER-EFF/pair_eff_cut.cpp
View File

@ -31,6 +31,7 @@
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "atom_vec_electron.h"
using namespace LAMMPS_NS;
@ -76,7 +77,7 @@ void PairEffCut::compute(int eflag, int vflag)
double rsq,rc;
int *ilist,*jlist,*numneigh,**firstneigh;
energy = eke = epauli = ecoul = errestrain = 0.0;
energy = eke = epauli = ecp_epauli = ecoul = errestrain = 0.0;
// pvector = [KE, Pauli, ecoul, radial_restraint]
for (i=0; i<4; i++) pvector[i] = 0.0;
@ -92,6 +93,8 @@ void PairEffCut::compute(int eflag, int vflag)
int *type = atom->type;
int nlocal = atom->nlocal;
int *id = atom->tag;
int newton_pair = force->newton_pair;
double qqrd2e = force->qqrd2e;
@ -101,7 +104,6 @@ void PairEffCut::compute(int eflag, int vflag)
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
@ -154,7 +156,7 @@ void PairEffCut::compute(int eflag, int vflag)
// Tally energy and compute radial atomic virial contribution
if (evflag) {
ev_tally_eff(i,i,nlocal,newton_pair,energy,0.0);
if (flexible_pressure_flag) // iff flexible pressure flag on
if (pressure_with_evirials_flag) // iff flexible pressure flag on
ev_tally_eff(i,i,nlocal,newton_pair,0.0,e1rforce*eradius[i]);
}
if (eflag_global) {
@ -221,6 +223,12 @@ void PairEffCut::compute(int eflag, int vflag)
ElecCoreNuc(qxq, rc, eradius[i], &ecoul, &fpair);
}
else if (spin[i] == 4 && spin[j] == 0) {
double qxq = q[i]*q[j];
ElecCoreNuc(qxq, rc, eradius[i], &ecoul, &fpair);
}
// nucleus (i) - pseudo-core nucleus (j) interaction
else if (spin[i] == 0 && spin[j] == 3) {
double qxq = q[i]*q[j];
@ -228,6 +236,12 @@ void PairEffCut::compute(int eflag, int vflag)
ElecCoreNuc(qxq, rc, eradius[j], &ecoul, &fpair);
}
else if (spin[i] == 0 && spin[j] == 4) {
double qxq = q[i]*q[j];
ElecCoreNuc(qxq, rc, eradius[j], &ecoul, &fpair);
}
// nucleus (i) - electron (j) Coul interaction
else if (spin[i] == 0 && abs(spin[j]) == 1) {
@ -239,7 +253,7 @@ void PairEffCut::compute(int eflag, int vflag)
erforce[j] += e1rforce;
// Radial electron virial, iff flexible pressure flag set
if (evflag && flexible_pressure_flag) {
if (evflag && pressure_with_evirials_flag) {
e1rvirial = eradius[j] * e1rforce;
ev_tally_eff(j,j,nlocal,newton_pair,0.0,e1rvirial);
}
@ -256,7 +270,7 @@ void PairEffCut::compute(int eflag, int vflag)
erforce[i] += e1rforce;
// Radial electron virial, iff flexible pressure flag set
if (evflag && flexible_pressure_flag) {
if (evflag && pressure_with_evirials_flag) {
e1rvirial = eradius[i] * e1rforce;
ev_tally_eff(i,i,nlocal,newton_pair,0.0,e1rvirial);
}
@ -283,7 +297,7 @@ void PairEffCut::compute(int eflag, int vflag)
erforce[j] += e2rforce;
// Radial electron virial, iff flexible pressure flag set
if (evflag && flexible_pressure_flag) {
if (evflag && pressure_with_evirials_flag) {
e1rvirial = eradius[i] * e1rforce;
e2rvirial = eradius[j] * e2rforce;
ev_tally_eff(i,j,nlocal,newton_pair,0.0,e1rvirial+e2rvirial);
@ -315,7 +329,7 @@ void PairEffCut::compute(int eflag, int vflag)
erforce[j] += e2rforce;
// Radial electron virial, iff flexible pressure flag set
if (evflag && flexible_pressure_flag) {
if (evflag && pressure_with_evirials_flag) {
e2rvirial = eradius[j] * e2rforce;
ev_tally_eff(j,j,nlocal,newton_pair,0.0,e2rvirial);
}
@ -347,7 +361,7 @@ void PairEffCut::compute(int eflag, int vflag)
erforce[i] += e2rforce;
// add radial atomic virial, iff flexible pressure flag set
if (evflag && flexible_pressure_flag) {
if (evflag && pressure_with_evirials_flag) {
e2rvirial = eradius[i] * e2rforce;
ev_tally_eff(i,i,nlocal,newton_pair,0.0,e2rvirial);
}
@ -391,25 +405,48 @@ void PairEffCut::compute(int eflag, int vflag)
else if (spin[i] == 3 && (abs(spin[j]) == 1 || spin[j] == 2)) {
e2rforce = ecp_e2rforce = 0.0;
if (abs(spin[j]) == 1) {
ElecCoreElec(q[i],rc,eradius[i],eradius[j],&ecoul,
&fpair,&e2rforce);
PauliCoreElec(rc,eradius[j],&ecp_epauli,&ecp_fpair,
&ecp_e2rforce,PAULI_CORE_A, PAULI_CORE_B,
PAULI_CORE_C);
} else { // add second s electron contribution from fixed-core
double qxq = q[i]*q[j];
ElecCoreNuc(qxq, rc, eradius[j], &ecoul, &fpair);
ElecCoreElec(q[i],rc,eradius[i],eradius[j],&ecoul,
&fpair,&e2rforce);
ElecCoreElec(q[i],rc,eradius[i],eradius[j],&ecoul,
&fpair,&e2rforce);
PauliCoreElec(rc,eradius[j],&ecp_epauli,&ecp_fpair,
&ecp_e2rforce,PAULI_CORE_A, PAULI_CORE_B,
PAULI_CORE_C);
PauliCoreElec(rc,eradius[j],&ecp_epauli,&ecp_fpair,
&ecp_e2rforce,PAULI_CORE_A, PAULI_CORE_B,
PAULI_CORE_C);
if (((PAULI_CORE_D[ecp_type[itype]]) == 0.0) && ((PAULI_CORE_E[ecp_type[itype]]) == 0.0)) {
if (abs(spin[j]) == 1) {
ElecCoreElec(q[i],rc,eradius[i],eradius[j],&ecoul,
&fpair,&e2rforce);
PauliCoreElec(rc,eradius[j],&ecp_epauli,&ecp_fpair,
&ecp_e2rforce,PAULI_CORE_A[ecp_type[itype]], PAULI_CORE_B[ecp_type[itype]],
PAULI_CORE_C[ecp_type[itype]]);
} else { // add second s electron contribution from fixed-core
double qxq = q[i]*q[j];
ElecCoreNuc(qxq, rc, eradius[j], &ecoul, &fpair);
ElecCoreElec(q[i],rc,eradius[i],eradius[j],&ecoul,
&fpair,&e2rforce);
ElecCoreElec(q[i],rc,eradius[i],eradius[j],&ecoul,
&fpair,&e2rforce);
PauliCoreElec(rc,eradius[j],&ecp_epauli,&ecp_fpair,
&ecp_e2rforce,PAULI_CORE_A[ecp_type[itype]], PAULI_CORE_B[ecp_type[itype]],
PAULI_CORE_C[ecp_type[itype]]);
PauliCoreElec(rc,eradius[j],&ecp_epauli,&ecp_fpair,
&ecp_e2rforce,PAULI_CORE_A[ecp_type[itype]], PAULI_CORE_B[ecp_type[itype]],
PAULI_CORE_C[ecp_type[itype]]);
}
} else {
if (abs(spin[j]) == 1) {
ElecCoreElec(q[i],rc,eradius[i],eradius[j],&ecoul,
&fpair,&e2rforce);
PauliCorePElec(rc,eradius[j],&ecp_epauli,&ecp_fpair,
&ecp_e2rforce,PAULI_CORE_A[ecp_type[itype]],PAULI_CORE_B[ecp_type[itype]],
PAULI_CORE_C[ecp_type[itype]],PAULI_CORE_D[ecp_type[itype]],PAULI_CORE_E[ecp_type[itype]]);
} else { // add second s electron contribution from fixed-core
double qxq = q[i]*q[j];
ElecCoreNuc(qxq, rc, eradius[j], &ecoul, &fpair);
ElecCoreElec(q[i],rc,eradius[i],eradius[j],&ecoul,
&fpair,&e2rforce);
ElecCoreElec(q[i],rc,eradius[i],eradius[j],&ecoul,
&fpair,&e2rforce);
PauliCorePElec(rc,eradius[j],&ecp_epauli,&ecp_fpair,
&ecp_e2rforce,PAULI_CORE_A[ecp_type[itype]], PAULI_CORE_B[ecp_type[itype]],
PAULI_CORE_C[ecp_type[itype]],PAULI_CORE_D[ecp_type[itype]],PAULI_CORE_E[ecp_type[itype]]);
PauliCorePElec(rc,eradius[j],&ecp_epauli,&ecp_fpair,
&ecp_e2rforce,PAULI_CORE_A[ecp_type[itype]], PAULI_CORE_B[ecp_type[itype]],
PAULI_CORE_C[ecp_type[itype]],PAULI_CORE_D[ecp_type[itype]],PAULI_CORE_E[ecp_type[itype]]);
}
}
// Apply conversion factor from Hartree to kcal/mol
@ -421,7 +458,7 @@ void PairEffCut::compute(int eflag, int vflag)
erforce[j] += e2rforce;
// add radial atomic virial, iff flexible pressure flag set
if (evflag && flexible_pressure_flag) {
if (evflag && pressure_with_evirials_flag) {
e2rvirial = eradius[j] * e2rforce;
ev_tally_eff(j,j,nlocal,newton_pair,0.0,e2rvirial);
}
@ -432,25 +469,48 @@ void PairEffCut::compute(int eflag, int vflag)
else if ((abs(spin[i]) == 1 || spin[i] == 2) && spin[j] == 3) {
e1rforce = ecp_e1rforce = 0.0;
if (abs(spin[i]) == 1) {
ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul,
&fpair,&e1rforce);
PauliCoreElec(rc,eradius[i],&ecp_epauli,&ecp_fpair,
&ecp_e1rforce,PAULI_CORE_A,PAULI_CORE_B,
PAULI_CORE_C);
if (((PAULI_CORE_D[ecp_type[jtype]]) == 0.0) && ((PAULI_CORE_E[ecp_type[jtype]]) == 0.0)) {
if (abs(spin[i]) == 1) {
ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul,
&fpair,&e1rforce);
PauliCoreElec(rc,eradius[i],&ecp_epauli,&ecp_fpair,
&ecp_e1rforce,PAULI_CORE_A[ecp_type[jtype]],PAULI_CORE_B[ecp_type[jtype]],
PAULI_CORE_C[ecp_type[jtype]]);
} else {
double qxq = q[i]*q[j];
ElecCoreNuc(qxq,rc,eradius[i],&ecoul,&fpair);
ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul,
&fpair,&e1rforce);
ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul,
&fpair,&e1rforce);
PauliCoreElec(rc,eradius[i],&ecp_epauli,&ecp_fpair,
&ecp_e1rforce,PAULI_CORE_A[ecp_type[jtype]], PAULI_CORE_B[ecp_type[jtype]],
PAULI_CORE_C[ecp_type[jtype]]);
PauliCoreElec(rc,eradius[i],&ecp_epauli,&ecp_fpair,
&ecp_e1rforce,PAULI_CORE_A[ecp_type[jtype]], PAULI_CORE_B[ecp_type[jtype]],
PAULI_CORE_C[ecp_type[jtype]]);
}
} else {
double qxq = q[i]*q[j];
ElecCoreNuc(qxq,rc,eradius[i],&ecoul,&fpair);
ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul,
&fpair,&e1rforce);
ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul,
&fpair,&e1rforce);
PauliCoreElec(rc,eradius[i],&ecp_epauli,&ecp_fpair,
&ecp_e1rforce,PAULI_CORE_A, PAULI_CORE_B,
PAULI_CORE_C);
PauliCoreElec(rc,eradius[i],&ecp_epauli,&ecp_fpair,
&ecp_e1rforce,PAULI_CORE_A, PAULI_CORE_B,
PAULI_CORE_C);
if (abs(spin[i]) == 1) {
ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul,
&fpair,&e1rforce);
PauliCorePElec(rc,eradius[i],&ecp_epauli,&ecp_fpair,
&ecp_e1rforce,PAULI_CORE_A[ecp_type[jtype]],PAULI_CORE_B[ecp_type[jtype]],
PAULI_CORE_C[ecp_type[jtype]],PAULI_CORE_D[ecp_type[jtype]],PAULI_CORE_E[ecp_type[jtype]]);
} else {
double qxq = q[i]*q[j];
ElecCoreNuc(qxq,rc,eradius[i],&ecoul,&fpair);
ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul,
&fpair,&e1rforce);
ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul,
&fpair,&e1rforce);
PauliCorePElec(rc,eradius[i],&ecp_epauli,&ecp_fpair,
&ecp_e1rforce,PAULI_CORE_A[ecp_type[jtype]], PAULI_CORE_B[ecp_type[jtype]],
PAULI_CORE_C[ecp_type[jtype]],PAULI_CORE_D[ecp_type[jtype]],PAULI_CORE_E[ecp_type[jtype]]);
PauliCorePElec(rc,eradius[i],&ecp_epauli,&ecp_fpair,
&ecp_e1rforce,PAULI_CORE_A[ecp_type[jtype]], PAULI_CORE_B[ecp_type[jtype]],
PAULI_CORE_C[ecp_type[jtype]],PAULI_CORE_D[ecp_type[jtype]],PAULI_CORE_E[ecp_type[jtype]]);
}
}
// Apply conversion factor from Hartree to kcal/mol
@ -462,7 +522,7 @@ void PairEffCut::compute(int eflag, int vflag)
erforce[i] += e1rforce;
// add radial atomic virial, iff flexible pressure flag set
if (evflag && flexible_pressure_flag) {
if (evflag && pressure_with_evirials_flag) {
e1rvirial = eradius[i] * e1rforce;
ev_tally_eff(i,i,nlocal,newton_pair,0.0,e1rvirial);
}
@ -470,7 +530,7 @@ void PairEffCut::compute(int eflag, int vflag)
// pseudo-core (i) - pseudo-core (j) interactions
else if (spin[i] == 3 && abs(spin[j]) == 3) {
else if (spin[i] == 3 && spin[j] == 3) {
double qxq = q[i]*q[j];
ElecCoreCore(qxq,rc,eradius[i],eradius[j],&ecoul,&fpair);
@ -527,7 +587,7 @@ void PairEffCut::compute(int eflag, int vflag)
// limit electron stifness (size) for periodic systems, to max=half-box-size
if (abs(spin[i]) == 1 && limit_size_flag) {
if (abs(spin[i]) == 1 && limit_eradius_flag) {
double half_box_length=0, dr, kfactor=hhmss2e*1.0;
e1rforce = errestrain = 0.0;
@ -548,7 +608,7 @@ void PairEffCut::compute(int eflag, int vflag)
// Tally radial restrain energy and add radial restrain virial
if (evflag) {
ev_tally_eff(i,i,nlocal,newton_pair,errestrain,0.0);
if (flexible_pressure_flag) // flexible electron pressure
if (pressure_with_evirials_flag) // flexible electron pressure
ev_tally_eff(i,i,nlocal,newton_pair,0.0,eradius[i]*e1rforce);
}
}
@ -558,7 +618,7 @@ void PairEffCut::compute(int eflag, int vflag)
}
if (vflag_fdotr) {
virial_fdotr_compute();
if (flexible_pressure_flag) virial_eff_compute();
if (pressure_with_evirials_flag) virial_eff_compute();
}
}
@ -700,46 +760,80 @@ void PairEffCut::allocate()
void PairEffCut::settings(int narg, char **arg)
{
if (narg != 1 && narg != 3 && narg != 5 && narg != 7)
if (narg < 1)
error->all(FLERR,"Illegal pair_style command");
// Defaults ECP parameters for Si
PAULI_CORE_A = 0.320852;
PAULI_CORE_B = 2.283269;
PAULI_CORE_C = 0.814857;
// Defaults ECP parameters for C (radius=0.154)
PAULI_CORE_A[6] = 22.721015;
PAULI_CORE_B[6] = 0.728733;
PAULI_CORE_C[6] = 1.103199;
PAULI_CORE_D[6] = 17.695345;
PAULI_CORE_E[6] = 6.693621;
if (narg == 1) {
cut_global = force->numeric(arg[0]);
limit_size_flag = 0;
flexible_pressure_flag = 0;
} else if (narg == 3) {
cut_global = force->numeric(arg[0]);
limit_size_flag = force->inumeric(arg[1]);
flexible_pressure_flag = force->inumeric(arg[2]);
} else if (narg == 5) {
cut_global = force->numeric(arg[0]);
limit_size_flag = 0;
flexible_pressure_flag = 0;
if (strcmp(arg[1],"ecp") != 0)
error->all(FLERR,"Illegal pair_style command");
else {
PAULI_CORE_A = force->numeric(arg[2]);
PAULI_CORE_B = force->numeric(arg[3]);
PAULI_CORE_C = force->numeric(arg[4]);
// Defaults ECP parameters for N (radius=0.394732)
PAULI_CORE_A[7] = 16.242367;
PAULI_CORE_B[7] = 0.602818;
PAULI_CORE_C[7] = 1.081856;
PAULI_CORE_D[7] = 7.150803;
PAULI_CORE_E[7] = 5.351936;
// Defaults p-element ECP parameters for Oxygen (radius=0.15)
PAULI_CORE_A[8] = 29.5185;
PAULI_CORE_B[8] = 0.32995;
PAULI_CORE_C[8] = 1.21676;
PAULI_CORE_D[8] = 11.98757;
PAULI_CORE_E[8] = 3.073417;
// Defaults ECP parameters for Al (radius=1.660)
PAULI_CORE_A[13] = 0.486;
PAULI_CORE_B[13] = 1.049;
PAULI_CORE_C[13] = 0.207;
PAULI_CORE_D[13] = 0.0;
PAULI_CORE_E[13] = 0.0;
// Defaults ECP parameters for Si (radius=1.691)
PAULI_CORE_A[14] = 0.320852;
PAULI_CORE_B[14] = 2.283269;
PAULI_CORE_C[14] = 0.814857;
PAULI_CORE_D[14] = 0.0;
PAULI_CORE_E[14] = 0.0;
cut_global = force->numeric(arg[0]);
limit_eradius_flag = 0;
pressure_with_evirials_flag = 0;
int atype;
int iarg = 1;
int ecp_found = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"limit_eradius") == 0) {
limit_eradius_flag = 1;
iarg += 1;
}
} else if (narg == 7) {
cut_global = force->numeric(arg[0]);
limit_size_flag = force->inumeric(arg[1]);
flexible_pressure_flag = force->inumeric(arg[2]);
if (strcmp(arg[3],"ecp") != 0)
error->all(FLERR,"Illegal pair_style command");
else {
PAULI_CORE_A = force->numeric(arg[4]);
PAULI_CORE_B = force->numeric(arg[5]);
PAULI_CORE_C = force->numeric(arg[6]);
else if (strcmp(arg[iarg],"pressure_with_evirials") == 0) {
pressure_with_evirials_flag = 1;
iarg += 1;
}
else if (strcmp(arg[iarg],"ecp") == 0) {
iarg += 1;
while (iarg < narg) {
atype = force->inumeric(arg[iarg]);
if (strcmp(arg[iarg+1],"C") == 0) ecp_type[atype] = 6;
else if (strcmp(arg[iarg+1],"N") == 0) ecp_type[atype] = 7;
else if (strcmp(arg[iarg+1],"O") == 0) ecp_type[atype] = 8;
else if (strcmp(arg[iarg+1],"Al") == 0) ecp_type[atype] = 13;
else if (strcmp(arg[iarg+1],"Si") == 0) ecp_type[atype] = 14;
else error->all(FLERR, "Note: there are no default parameters for this atom ECP\n");
iarg += 2;
ecp_found = 1;
}
}
}
if (!ecp_found && atom->ecp_flag)
error->all(FLERR,"Need to specify ECP type on pair_style command");
// Need to introduce 2 new constants w/out changing update.cpp
if (force->qqr2e==332.06371) { // i.e. Real units chosen
h2e = 627.509; // hartree->kcal/mol
@ -759,34 +853,6 @@ void PairEffCut::settings(int narg, char **arg)
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairEffCut::coeff(int narg, char **arg)
{
if (narg < 2 || narg > 3) error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
force->bounds(arg[0],atom->ntypes,ilo,ihi);
force->bounds(arg[1],atom->ntypes,jlo,jhi);
double cut_one = cut_global;
if (narg == 3) cut_one = atof(arg[2]);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
cut[i][j] = cut_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
@ -817,6 +883,51 @@ void PairEffCut::init_style()
int irequest = neighbor->request(this);
}
/* ----------------------------------------------------------------------
set coeffs for one or more type electron pairs (ECP-only)
------------------------------------------------------------------------- */
void PairEffCut::coeff(int narg, char **arg)
{
if (!allocated) allocate();
if ((strcmp(arg[0],"*") == 0) || (strcmp(arg[1],"*") == 0)) {
int ilo,ihi,jlo,jhi;
force->bounds(arg[0],atom->ntypes,ilo,ihi);
force->bounds(arg[1],atom->ntypes,jlo,jhi);
double cut_one = cut_global;
if (narg == 3) cut_one = atof(arg[2]);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
cut[i][j] = cut_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
} else {
int ecp;
ecp = force->inumeric(arg[0]);
if (strcmp(arg[1],"s") ==0) {
PAULI_CORE_A[ecp_type[ecp]] = force->numeric(arg[2]);
PAULI_CORE_B[ecp_type[ecp]] = force->numeric(arg[3]);
PAULI_CORE_C[ecp_type[ecp]] = force->numeric(arg[4]);
PAULI_CORE_D[ecp_type[ecp]] = 0.0;
PAULI_CORE_E[ecp_type[ecp]] = 0.0;
} else if (strcmp(arg[1],"p") ==0) {
PAULI_CORE_A[ecp_type[ecp]] = force->numeric(arg[2]);
PAULI_CORE_B[ecp_type[ecp]] = force->numeric(arg[3]);
PAULI_CORE_C[ecp_type[ecp]] = force->numeric(arg[4]);
PAULI_CORE_D[ecp_type[ecp]] = force->numeric(arg[5]);
PAULI_CORE_E[ecp_type[ecp]] = force->numeric(arg[6]);
} else error->all(FLERR,"Illegal pair_coeff command");
}
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */

5
src/USER-EFF/pair_eff_cut.h
View File

@ -45,10 +45,11 @@ class PairEffCut : public Pair {
double memory_usage();
private:
int limit_size_flag, flexible_pressure_flag;
int limit_eradius_flag, pressure_with_evirials_flag;
double cut_global;
double **cut;
double PAULI_CORE_A, PAULI_CORE_B, PAULI_CORE_C;
int ecp_type[100];
double PAULI_CORE_A[100], PAULI_CORE_B[100], PAULI_CORE_C[100], PAULI_CORE_D[100], PAULI_CORE_E[100];
double hhmss2e, h2e;
int nmax;

36
src/USER-EFF/pair_eff_inline.h
View File

@ -494,6 +494,42 @@ inline void PauliCoreElec(double rc, double re2, double *epauli, double *frc,
/* ---------------------------------------------------------------------- */
inline void PauliCorePElec(double rc, double re2, double *epauli, double *frc,
double *fre2, double PAULI_CORE_P_A,
double PAULI_CORE_P_B, double PAULI_CORE_P_C,
double PAULI_CORE_P_D, double PAULI_CORE_P_E)
{
double E, dEdrc, dEdre2;
E = PAULI_CORE_P_A *
pow((2.0 / (PAULI_CORE_P_B / re2 + re2 / PAULI_CORE_P_B)), 5.0) *
pow((rc - PAULI_CORE_P_C * re2), 2.0) * exp(-PAULI_CORE_P_D *
pow((rc - PAULI_CORE_P_C * re2), 2.0) / (PAULI_CORE_P_E + re2 * re2));
dEdrc = PAULI_CORE_P_A *
pow((2.0 / (PAULI_CORE_P_B / re2 + re2 / PAULI_CORE_P_B)), 5.0) *
2.0 * (rc - PAULI_CORE_P_C * re2) * exp(-PAULI_CORE_P_D *
pow((rc - PAULI_CORE_P_C * re2), 2.0) / (PAULI_CORE_P_E + re2 * re2))
+ E * (-PAULI_CORE_P_D * 2.0 * (rc - PAULI_CORE_P_C * re2) /
(PAULI_CORE_P_E + re2 * re2));
dEdre2 = E * (-5.0 / (PAULI_CORE_P_B / re2 + re2 / PAULI_CORE_P_B) *
(-PAULI_CORE_P_B / (re2 * re2) + 1.0 / PAULI_CORE_P_B))
+ PAULI_CORE_P_A *
pow((2.0 / (PAULI_CORE_P_B / re2 + re2 / PAULI_CORE_P_B)), 5.0) *
2.0 * (rc - PAULI_CORE_P_C * re2) * (-PAULI_CORE_P_C) *
exp(-PAULI_CORE_P_D * pow((rc - PAULI_CORE_P_C * re2), 2.0) /
(PAULI_CORE_P_E + re2 * re2)) + E * (2.0 * PAULI_CORE_P_D *
(rc - PAULI_CORE_P_C * re2) * (PAULI_CORE_P_C * PAULI_CORE_P_E +
rc * re2) / pow((PAULI_CORE_P_E + re2 * re2), 2.0));
*epauli += E;
*frc -= dEdrc;
*fre2 -= dEdre2;
}
/* ---------------------------------------------------------------------- */
inline void RForce(double dx, double dy, double dz,
double rc, double force, double *fx, double *fy, double *fz)
{