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@5716 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2011-02-22 22:36:09 +00:00
parent 1d4a6440d1
commit 691894ab34
25 changed files with 1319 additions and 624 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/USER-EFF/Install.sh
View File

@ -44,7 +44,7 @@ elif (test $1 = 0) then
rm ../compute_temp_deform_eff.cpp
rm ../compute_temp_eff.cpp
rm ../compute_temp_region_eff.cpp
rm ../fix_langevin_eff.cpp
rm ../fix_langevin_eff.cpp ..
rm ../fix_nh_eff.cpp
rm ../fix_nve_eff.cpp
rm ../fix_nvt_eff.cpp
@ -61,7 +61,7 @@ elif (test $1 = 0) then
rm ../compute_temp_deform_eff.h
rm ../compute_temp_eff.h
rm ../compute_temp_region_eff.h
rm ../fix_langevin_eff.h
rm ../fix_langevin_eff.h ..
rm ../fix_nh_eff.h
rm ../fix_nve_eff.h
rm ../fix_nvt_eff.h

21
src/USER-EFF/README
View File

@ -54,6 +54,27 @@ OTHERS FILES INCLUDED:
User examples are under examples/USER/eff
eFF tools are under tools/eff
ACKNOWLEDGMENTS:
Thanks to Steve Plimpton and Aidan Thompson for their input on the
LAMMPS architecture and for their help in customizing some of the
required LAMMPS core modules.
Version 01/2010: Special thanks to:
- Hai Xiao (Caltech) for reviewing the fixed-core implementation and
providing useful insights to improve it, and for his work on the effective core pseudopotential.
- Vaclav Cvicek (Caltech) for thoroughly revising the units, for finding a bug in the
fix_langevin_eff radial scaling factors, and for suggesting changes to clean-up the code.
- Patrick Theofanis (Caltech) for providing an optimized set of parameters for the Si ECP
(default) and for providing basic cases.
- Qi An (Caltech) for providing feedback on usage, application cases, and testing.
VERSION NOTES:
01/2010: Added support for fixed-core and effective core pseudopotentials [ECP]
(useful for C, Al, Si, O and other elements). Cleaned up the code to make it
easier to maintain, revised support for real units, upgraded post-processing
and visualization tools, added support for "compute pair eff" to allow thermo
prints with the different eFF energy components (eke, epauli, ecoul and errestrain),
fixed radial scaling factors in the eff langevin thermostat.

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

@ -15,7 +15,7 @@
Contributing author: Andres Jaramillo-Botero (Caltech)
------------------------------------------------------------------------- */
#include "lmptype.h"
#include "math.h"
#include "stdlib.h"
#include "atom_vec_electron.h"
#include "atom.h"
@ -63,9 +63,7 @@ void AtomVecElectron::grow(int n)
if (n == 0) nmax += DELTA;
else nmax = n;
atom->nmax = nmax;
if (nmax < 0 || nmax > MAXSMALLINT)
error->one("Per-processor system is too big");
tag = atom->tag = (int *)
memory->srealloc(atom->tag,nmax*sizeof(int),"atom:tag");
type = atom->type = (int *)
@ -104,8 +102,7 @@ void AtomVecElectron::grow_reset()
tag = atom->tag; type = atom->type;
mask = atom->mask; image = atom->image;
x = atom->x; v = atom->v; f = atom->f;
q = atom->q;
eradius = atom->eradius; ervel = atom->ervel; erforce = atom->erforce;
q = atom->q; eradius = atom->eradius; ervel = atom->ervel; erforce = atom->erforce;
}
/* ---------------------------------------------------------------------- */

21
src/USER-EFF/compute_ke_atom_eff.cpp
View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Andres Jaramillo-Botero (Caltech)
Contributing author: Andres Jaramillo-Botero
------------------------------------------------------------------------- */
#include "string.h"
@ -86,28 +86,17 @@ void ComputeKEAtomEff::compute_peratom()
double **v = atom->v;
double *ervel = atom->ervel;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *spin = atom->spin;
int *mask = atom->mask;
int *type = atom->type;
int nlocal = atom->nlocal;
if (rmass)
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
ke[i] = 0.5 * mvv2e * rmass[i] *
(v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
if (spin[i])
ke[i] += 0.5 * mvv2e * rmass[i] * ervel[i]*ervel[i] * 0.75;
} else ke[i] = 0.0;
}
else
if (mass)
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
ke[i] = 0.5 * mvv2e * mass[type[i]] *
(v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
if (spin[i])
(v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
if (abs(spin[i])==1)
ke[i] += 0.5 * mvv2e * mass[type[i]] * ervel[i]*ervel[i] * 0.75;
} else ke[i] = 0.0;
}

16
src/USER-EFF/compute_ke_eff.cpp
View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Andres Jaramillo-Botero (Caltech)
Contributing author: Andres Jaramillo-Botero
------------------------------------------------------------------------- */
#include "mpi.h"
@ -57,7 +57,6 @@ double ComputeKEEff::compute_scalar()
double **v = atom->v;
double *ervel = atom->ervel;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *spin = atom->spin;
int *mask = atom->mask;
@ -66,19 +65,12 @@ double ComputeKEEff::compute_scalar()
double ke = 0.0;
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
ke += rmass[i] * (v[i][0]*v[i][0] + v[i][1]*v[i][1] +
v[i][2]*v[i][2]);
if (spin[i]) ke += 0.75*rmass[i]*ervel[i]*ervel[i];
}
} else {
if (mass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
ke += mass[type[i]]*(v[i][0]*v[i][0] + v[i][1]*v[i][1] +
v[i][2]*v[i][2]);
if (spin[i]) ke += 0.75*mass[type[i]]*ervel[i]*ervel[i];
v[i][2]*v[i][2]);
if (abs(spin[i])==1) ke += 0.75*mass[type[i]]*ervel[i]*ervel[i];
}
}

180
src/USER-EFF/compute_temp_deform_eff.cpp
View File

@ -12,19 +12,25 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Andres Jaramillo-Botero (Caltech))
Contributing author: Andres Jaramillo-Botero (Caltech)
------------------------------------------------------------------------- */
#include "mpi.h"
#include "compute_temp_deform_eff.h"
#include "update.h"
#include "atom.h"
#include "string.h"
#include "domain.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "modify.h"
#include "fix.h"
#include "fix_deform.h"
#include "group.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp
@ -32,12 +38,56 @@ enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp
/* ---------------------------------------------------------------------- */
ComputeTempDeformEff::ComputeTempDeformEff(LAMMPS *lmp, int narg, char **arg) :
ComputeTempDeform(lmp, narg, arg)
Compute(lmp, narg, arg)
{
// error check
if (narg != 3) error->all("Illegal compute temp/deform/eff command");
if (!atom->spin_flag || !atom->ervel_flag)
error->all("Compute temp/deform/eff requires atom attributes spin, ervel");
scalar_flag = vector_flag = 1;
size_vector = 6;
extscalar = 0;
extvector = 1;
tempflag = 1;
tempbias = 1;
maxbias = 0;
vbiasall = NULL;
vector = new double[6];
}
/* ---------------------------------------------------------------------- */
ComputeTempDeformEff::~ComputeTempDeformEff()
{
memory->destroy_2d_double_array(vbiasall);
delete [] vector;
}
/* ---------------------------------------------------------------------- */
void ComputeTempDeformEff::init()
{
int i;
fix_dof = 0;
for (i = 0; i < modify->nfix; i++)
fix_dof += modify->fix[i]->dof(igroup);
dof_compute();
// check fix deform remap settings
for (i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"deform") == 0) {
if (((FixDeform *) modify->fix[i])->remapflag == X_REMAP &&
comm->me == 0)
error->warning("Using compute temp/deform/eff with inconsistent "
"fix deform remap option");
break;
}
if (i == modify->nfix && comm->me == 0)
error->warning("Using compute temp/deform/eff with no fix deform defined");
}
/* ---------------------------------------------------------------------- */
@ -57,7 +107,7 @@ void ComputeTempDeformEff::dof_compute()
int one = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (spin[i]) one++;
if (abs(spin[i])==1) one++;
}
int nelectrons;
MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world);
@ -80,7 +130,6 @@ double ComputeTempDeformEff::compute_scalar()
double **v = atom->v;
double *ervel = atom->ervel;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *spin = atom->spin;
int *type = atom->type;
int *mask = atom->mask;
@ -99,21 +148,17 @@ double ComputeTempDeformEff::compute_scalar()
if (mask[i] & groupbit) {
domain->x2lamda(x[i],lamda);
vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] +
h_rate[4]*lamda[2] + h_ratelo[0];
h_rate[4]*lamda[2] + h_ratelo[0];
vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2];
vthermal[0] = v[i][0] - vstream[0];
vthermal[1] = v[i][1] - vstream[1];
vthermal[2] = v[i][2] - vstream[2];
if (rmass) {
t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] +
vthermal[2]*vthermal[2]) * rmass[i];
if (spin[i]) t += 0.75*rmass[i]*ervel[i]*ervel[i];
} else {
t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] +
vthermal[2]*vthermal[2])* mass[type[i]];
if (spin[i]) t += 0.75*mass[type[i]]*ervel[i]*ervel[i];
if (mass) {
t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] +
vthermal[2]*vthermal[2])* mass[type[i]];
if (abs(spin[i])==1) t += 0.75*mass[type[i]]*ervel[i]*ervel[i];
}
}
@ -135,7 +180,6 @@ void ComputeTempDeformEff::compute_vector()
double **v = atom->v;
double *ervel = atom->ervel;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *spin = atom->spin;
int *type = atom->type;
int *mask = atom->mask;
@ -151,22 +195,21 @@ void ComputeTempDeformEff::compute_vector()
if (mask[i] & groupbit) {
domain->x2lamda(x[i],lamda);
vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] +
h_rate[4]*lamda[2] + h_ratelo[0];
h_rate[4]*lamda[2] + h_ratelo[0];
vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2];
vthermal[0] = v[i][0] - vstream[0];
vthermal[1] = v[i][1] - vstream[1];
vthermal[2] = v[i][2] - vstream[2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
massone = mass[type[i]];
t[0] += massone * vthermal[0]*vthermal[0];
t[1] += massone * vthermal[1]*vthermal[1];
t[2] += massone * vthermal[2]*vthermal[2];
t[3] += massone * vthermal[0]*vthermal[1];
t[4] += massone * vthermal[0]*vthermal[2];
t[5] += massone * vthermal[1]*vthermal[2];
if (spin[i]) {
if (abs(spin[i])==1) {
t[0] += 0.75 * massone * ervel[i]*ervel[i];
t[1] += 0.75 * massone * ervel[i]*ervel[i];
t[2] += 0.75 * massone * ervel[i]*ervel[i];
@ -176,3 +219,98 @@ void ComputeTempDeformEff::compute_vector()
MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e;
}
/* ----------------------------------------------------------------------
remove velocity bias from atom I to leave thermal velocity
------------------------------------------------------------------------- */
void ComputeTempDeformEff::remove_bias(int i, double *v)
{
double lamda[3];
double *h_rate = domain->h_rate;
double *h_ratelo = domain->h_ratelo;
domain->x2lamda(atom->x[i],lamda);
vbias[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] +
h_rate[4]*lamda[2] + h_ratelo[0];
vbias[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
vbias[2] = h_rate[2]*lamda[2] + h_ratelo[2];
v[0] -= vbias[0];
v[1] -= vbias[1];
v[2] -= vbias[2];
}
/* ----------------------------------------------------------------------
remove velocity bias from all atoms to leave thermal velocity
NOTE: only removes translational velocity bias from electrons
------------------------------------------------------------------------- */
void ComputeTempDeformEff::remove_bias_all()
{
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (nlocal > maxbias) {
memory->destroy_2d_double_array(vbiasall);
maxbias = atom->nmax;
vbiasall = memory->create_2d_double_array(maxbias,3,
"temp/deform/eff:vbiasall");
}
double lamda[3];
double *h_rate = domain->h_rate;
double *h_ratelo = domain->h_ratelo;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
domain->x2lamda(atom->x[i],lamda);
vbiasall[i][0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] +
h_rate[4]*lamda[2] + h_ratelo[0];
vbiasall[i][1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
vbiasall[i][2] = h_rate[2]*lamda[2] + h_ratelo[2];
v[i][0] -= vbiasall[i][0];
v[i][1] -= vbiasall[i][1];
v[i][2] -= vbiasall[i][2];
}
}
/* ----------------------------------------------------------------------
add back in velocity bias to atom I removed by remove_bias()
assume remove_bias() was previously called
------------------------------------------------------------------------- */
void ComputeTempDeformEff::restore_bias(int i, double *v)
{
v[0] += vbias[0];
v[1] += vbias[1];
v[2] += vbias[2];
}
/* ----------------------------------------------------------------------
add back in velocity bias to all atoms removed by remove_bias_all()
assume remove_bias_all() was previously called
------------------------------------------------------------------------- */
void ComputeTempDeformEff::restore_bias_all()
{
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
v[i][0] += vbiasall[i][0];
v[i][1] += vbiasall[i][1];
v[i][2] += vbiasall[i][2];
}
}
/* ---------------------------------------------------------------------- */
double ComputeTempDeformEff::memory_usage()
{
double bytes = maxbias * sizeof(double);
return bytes;
}

27
src/USER-EFF/compute_temp_deform_eff.h
View File

@ -20,19 +20,32 @@ ComputeStyle(temp/deform/eff,ComputeTempDeformEff)
#ifndef LMP_COMPUTE_TEMP_DEFORM_EFF_H
#define LMP_COMPUTE_TEMP_DEFORM_EFF_H
#include "compute_temp_deform.h"
#include "compute.h"
namespace LAMMPS_NS {
class ComputeTempDeformEff : public ComputeTempDeform {
class ComputeTempDeformEff : public Compute {
public:
ComputeTempDeformEff(class LAMMPS *, int, char **);
~ComputeTempDeformEff() {}
double compute_scalar();
void compute_vector();
virtual ~ComputeTempDeformEff();
void init();
virtual double compute_scalar();
virtual void compute_vector();
private:
void dof_compute();
void remove_bias(int, double *);
void remove_bias_all();
void restore_bias(int, double *);
void restore_bias_all();
double memory_usage();
protected:
int fix_dof;
double tfactor;
double vbias[3]; // stored velocity bias for one atom
double **vbiasall; // stored velocity bias for all atoms
int maxbias; // size of vbiasall array
virtual void dof_compute();
};
}

53
src/USER-EFF/compute_temp_eff.cpp
View File

@ -32,12 +32,37 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeTempEff::ComputeTempEff(LAMMPS *lmp, int narg, char **arg) :
ComputeTemp(lmp, narg, arg)
Compute(lmp, narg, arg)
{
// error check
if (!atom->spin_flag || !atom->ervel_flag)
error->all("Compute temp/eff requires atom attributes spin, ervel");
scalar_flag = vector_flag = 1;
size_vector = 6;
extscalar = 0;
extvector = 1;
tempflag = 1;
vector = new double[size_vector];
}
/* ---------------------------------------------------------------------- */
ComputeTempEff::~ComputeTempEff()
{
delete [] vector;
}
/* ---------------------------------------------------------------------- */
void ComputeTempEff::init()
{
fix_dof = 0;
for (int i = 0; i < modify->nfix; i++)
fix_dof += modify->fix[i]->dof(igroup);
dof_compute();
}
/* ---------------------------------------------------------------------- */
@ -55,7 +80,7 @@ void ComputeTempEff::dof_compute()
int one = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (spin[i]) one++;
if (abs(spin[i])==1) one++;
}
int nelectrons;
MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world);
@ -77,29 +102,19 @@ double ComputeTempEff::compute_scalar()
double **v = atom->v;
double *ervel = atom->ervel;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *spin = atom->spin;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double massone;
double t = 0.0;
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2])*rmass[i];
if (spin[i]) t += 0.75*rmass[i]*ervel[i]*ervel[i];
}
}
} else {
if (mass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
mass[type[i]];
if (spin[i]) t += 0.75*mass[type[i]]*ervel[i]*ervel[i];
mass[type[i]];
if (abs(spin[i])==1) t += 0.75*mass[type[i]]*ervel[i]*ervel[i];
}
}
}
@ -121,26 +136,24 @@ void ComputeTempEff::compute_vector()
double **v = atom->v;
double *ervel = atom->ervel;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *spin = atom->spin;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double massone,ervel_adj,t[6];
double massone,t[6];
for (i = 0; i < 6; i++) t[i] = 0.0;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
massone = mass[type[i]];
t[0] += massone * v[i][0]*v[i][0];
t[1] += massone * v[i][1]*v[i][1];
t[2] += massone * v[i][2]*v[i][2];
t[3] += massone * v[i][0]*v[i][1];
t[4] += massone * v[i][0]*v[i][2];
t[5] += massone * v[i][1]*v[i][2];
if (spin[i]) {
if (abs(spin[i])==1) {
t[0] += 0.75*massone*ervel[i]*ervel[i];
t[1] += 0.75*massone*ervel[i]*ervel[i];
t[2] += 0.75*massone*ervel[i]*ervel[i];

11
src/USER-EFF/compute_temp_eff.h
View File

@ -20,18 +20,23 @@ ComputeStyle(temp/eff,ComputeTempEff)
#ifndef LMP_COMPUTE_TEMP_EFF_H
#define LMP_COMPUTE_TEMP_EFF_H
#include "compute_temp.h"
#include "compute.h"
namespace LAMMPS_NS {
class ComputeTempEff : public ComputeTemp {
class ComputeTempEff : public Compute {
public:
ComputeTempEff(class LAMMPS *, int, char **);
~ComputeTempEff() {}
virtual ~ComputeTempEff();
void init();
double compute_scalar();
void compute_vector();
private:
int fix_dof;
double tfactor;
double *inertia;
void dof_compute();
};

179
src/USER-EFF/compute_temp_region_eff.cpp
View File

@ -21,22 +21,72 @@
#include "atom.h"
#include "update.h"
#include "force.h"
#include "modify.h"
#include "domain.h"
#include "region.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeTempRegionEff::
ComputeTempRegionEff(LAMMPS *lmp, int narg, char **arg) :
ComputeTempRegion(lmp, narg, arg)
ComputeTempRegionEff::ComputeTempRegionEff(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
// error check
if (!atom->spin_flag || !atom->ervel_flag)
error->all("Compute temp/region/eff requires atom attributes spin, ervel");
if (narg != 4) error->all("Illegal compute temp/region/eff command");
iregion = domain->find_region(arg[3]);
if (iregion == -1)
error->all("Region ID for compute temp/region/eff does not exist");
int n = strlen(arg[3]) + 1;
idregion = new char[n];
strcpy(idregion,arg[3]);
scalar_flag = vector_flag = 1;
size_vector = 6;
extscalar = 0;
extvector = 1;
tempflag = 1;
tempbias = 1;
maxbias = 0;
vbiasall = NULL;
vector = new double[6];
}
/* ---------------------------------------------------------------------- */
ComputeTempRegionEff::~ComputeTempRegionEff()
{
delete [] idregion;
memory->destroy_2d_double_array(vbiasall);
delete [] vector;
}
/* ---------------------------------------------------------------------- */
void ComputeTempRegionEff::init()
{
// set index and check validity of region
iregion = domain->find_region(idregion);
if (iregion == -1)
error->all("Region ID for compute temp/region/eff does not exist");
dof = 0.0;
}
/* ---------------------------------------------------------------------- */
int ComputeTempRegionEff::dof_remove(int i)
{
double *x = atom->x[i];
if (domain->regions[iregion]->match(x[0],x[1],x[2])) return 0;
return 1;
}
/* ---------------------------------------------------------------------- */
@ -49,7 +99,6 @@ double ComputeTempRegionEff::compute_scalar()
double **v = atom->v;
double *ervel = atom->ervel;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *spin = atom->spin;
int *type = atom->type;
int *mask = atom->mask;
@ -59,20 +108,13 @@ double ComputeTempRegionEff::compute_scalar()
int count = 0;
double t = 0.0;
if (rmass) {
if (mass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
count++;
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2])*rmass[i];
if (spin[i]) t += 0.75*rmass[i]*ervel[i]*ervel[i];
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
count++;
count++;
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
mass[type[i]];
if (spin[i]) t += 0.75*mass[type[i]]*ervel[i]*ervel[i];
mass[type[i]];
if (abs(spin[i])==1) t += 0.75*mass[type[i]]*ervel[i]*ervel[i];
}
}
@ -85,7 +127,7 @@ double ComputeTempRegionEff::compute_scalar()
int one = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
if (spin[i]) one++;
if (abs(spin[i])==1) one++;
}
int nelectrons_region;
MPI_Allreduce(&one,&nelectrons_region,1,MPI_INT,MPI_SUM,world);
@ -111,7 +153,6 @@ void ComputeTempRegionEff::compute_vector()
double **v = atom->v;
double *ervel = atom->ervel;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *spin = atom->spin;
int *type = atom->type;
int *mask = atom->mask;
@ -123,8 +164,7 @@ void ComputeTempRegionEff::compute_vector()
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
massone = mass[type[i]];
t[0] += massone * v[i][0]*v[i][0];
t[1] += massone * v[i][1]*v[i][1];
@ -132,14 +172,105 @@ void ComputeTempRegionEff::compute_vector()
t[3] += massone * v[i][0]*v[i][1];
t[4] += massone * v[i][0]*v[i][2];
t[5] += massone * v[i][1]*v[i][2];
if (spin[i]) {
if (abs(spin[i])==1) {
t[0] += 0.75 * massone * ervel[i]*ervel[i];
t[1] += 0.75 * massone * ervel[i]*ervel[i];
t[2] += 0.75 * massone * ervel[i]*ervel[i];
}
}
}
MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < 6; i++) vector[i] *= force->mvv2e;
}
/* ----------------------------------------------------------------------
remove velocity bias from atom I to leave thermal velocity
NOTE: the following commands do not bias the radial electron velocities
------------------------------------------------------------------------- */
void ComputeTempRegionEff::remove_bias(int i, double *v)
{
double *x = atom->x[i];
if (domain->regions[iregion]->match(x[0],x[1],x[2]))
vbias[0] = vbias[1] = vbias[2] = 0.0;
else {
vbias[0] = v[0];
vbias[1] = v[1];
vbias[2] = v[2];
v[0] = v[1] = v[2] = 0.0;
}
}
/* ----------------------------------------------------------------------
remove velocity bias from all atoms to leave thermal velocity
------------------------------------------------------------------------- */
void ComputeTempRegionEff::remove_bias_all()
{
double **x = atom->x;
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (nlocal > maxbias) {
memory->destroy_2d_double_array(vbiasall);
maxbias = atom->nmax;
vbiasall = memory->create_2d_double_array(maxbias,3,
"temp/region:vbiasall");
}
Region *region = domain->regions[iregion];
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (region->match(x[i][0],x[i][1],x[i][2]))
vbiasall[i][0] = vbiasall[i][1] = vbiasall[i][2] = 0.0;
else {
vbiasall[i][0] = v[i][0];
vbiasall[i][1] = v[i][1];
vbiasall[i][2] = v[i][2];
v[i][0] = v[i][1] = v[i][2] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
add back in velocity bias to atom I removed by remove_bias()
assume remove_bias() was previously called
------------------------------------------------------------------------- */
void ComputeTempRegionEff::restore_bias(int i, double *v)
{
v[0] += vbias[0];
v[1] += vbias[1];
v[2] += vbias[2];
}
/* ----------------------------------------------------------------------
add back in velocity bias to all atoms removed by remove_bias_all()
assume remove_bias_all() was previously called
------------------------------------------------------------------------- */
void ComputeTempRegionEff::restore_bias_all()
{
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
v[i][0] += vbiasall[i][0];
v[i][1] += vbiasall[i][1];
v[i][2] += vbiasall[i][2];
}
}
/* ---------------------------------------------------------------------- */
double ComputeTempRegionEff::memory_usage()
{
double bytes = maxbias * sizeof(double);
return bytes;
}

22
src/USER-EFF/compute_temp_region_eff.h
View File

@ -20,16 +20,28 @@ ComputeStyle(temp/region/eff,ComputeTempRegionEff)
#ifndef LMP_COMPUTE_TEMP_REGION_EFF_H
#define LMP_COMPUTE_TEMP_REGION_EFF_H
#include "compute_temp_region.h"
#include "compute.h"
namespace LAMMPS_NS {
class ComputeTempRegionEff : public ComputeTempRegion {
class ComputeTempRegionEff : public Compute {
public:
ComputeTempRegionEff(class LAMMPS *, int, char **);
~ComputeTempRegionEff() {}
double compute_scalar();
void compute_vector();
virtual ~ComputeTempRegionEff();
void init();
virtual double compute_scalar();
virtual void compute_vector();
int dof_remove(int);
void remove_bias(int, double *);
void remove_bias_all();
void restore_bias(int, double *);
void restore_bias_all();
double memory_usage();
protected:
int iregion;
char *idregion;
};
}

224
src/USER-EFF/fix_langevin_eff.cpp
View File

@ -60,13 +60,13 @@ void FixLangevinEff::post_force_no_tally()
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double *ervel = atom->ervel;
double *erforce = atom->erforce;
int *spin = atom->spin;
double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep;
@ -80,78 +80,33 @@ void FixLangevinEff::post_force_no_tally()
// test v = 0 since some computes mask non-participating atoms via v = 0
// and added force has extra term not multiplied by v = 0
if (rmass) {
double boltz = force->boltz;
double dt = update->dt;
double mvv2e = force->mvv2e;
double ftm2v = force->ftm2v;
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
gamma1 = -rmass[i] / t_period / ftm2v;
gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
gamma1 *= 1.0/ratio[type[i]];
gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
erforce[i] += gamma1*ervel[i] + gamma2*(random->uniform()-0.5);
}
}
} else if (which == BIAS) {
double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
gamma1 = -rmass[i] / t_period / ftm2v;
gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
gamma1 *= 1.0/ratio[type[i]];
gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
temperature->remove_bias(i,v[i]);
if (v[i][0] != 0.0)
f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
if (v[i][1] != 0.0)
f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
if (v[i][2] != 0.0)
f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
if (ervel[i] != 0.0)
erforce[i] += gamma1*ervel[i] + gamma2*(random->uniform()-0.5);
temperature->restore_bias(i,v[i]);
}
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
if (abs(spin[i])==1) erforce[i] += 0.75*gamma1*ervel[i] + 0.866025404*gamma2*(random->uniform()-0.5);
}
}
} else {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
} else if (which == BIAS) {
double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
temperature->remove_bias(i,v[i]);
if (v[i][0] != 0.0)
f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
if (v[i][1] != 0.0)
f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
if (v[i][2] != 0.0)
f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
erforce[i] += gamma1*ervel[i] + gamma2*(random->uniform()-0.5);
}
}
} else if (which == BIAS) {
double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
temperature->remove_bias(i,v[i]);
if (v[i][0] != 0.0)
f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
if (v[i][1] != 0.0)
f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
if (v[i][2] != 0.0)
f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
if (ervel[i] != 0.0)
erforce[i] += gamma1*ervel[i] + gamma2*(random->uniform()-0.5);
temperature->restore_bias(i,v[i]);
}
if (abs(spin[i])==1 && ervel[i] != 0.0)
erforce[i] += 0.75*gamma1*ervel[i] + 0.866025404*gamma2*(random->uniform()-0.5);
temperature->restore_bias(i,v[i]);
}
}
}
@ -176,13 +131,13 @@ void FixLangevinEff::post_force_tally()
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double *erforce = atom->erforce;
double *ervel = atom->ervel;
int *spin = atom->spin;
double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep;
@ -196,92 +151,40 @@ void FixLangevinEff::post_force_tally()
// test v = 0 since some computes mask non-participating atoms via v = 0
// and added force has extra term not multiplied by v = 0
if (rmass) {
double boltz = force->boltz;
double dt = update->dt;
double mvv2e = force->mvv2e;
double ftm2v = force->ftm2v;
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
gamma1 = -rmass[i] / t_period / ftm2v;
gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
gamma1 *= 1.0/ratio[type[i]];
gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
erforcelangevin[i] = gamma1*ervel[i]+gamma2*(random->uniform()-0.5);
f[i][0] += flangevin[i][0];
f[i][1] += flangevin[i][1];
f[i][2] += flangevin[i][2];
erforce[i] += erforcelangevin[i];
}
}
} else if (which == BIAS) {
double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
gamma1 = -rmass[i] / t_period / ftm2v;
gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
gamma1 *= 1.0/ratio[type[i]];
gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
temperature->remove_bias(i,v[i]);
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
erforcelangevin[i] = gamma1*ervel[i]+gamma2*(random->uniform()-0.5);
if (v[i][0] != 0.0) f[i][0] += flangevin[i][0];
else flangevin[i][0] = 0;
if (v[i][1] != 0.0) f[i][1] += flangevin[i][1];
else flangevin[i][1] = 0;
if (v[i][2] != 0.0) f[i][2] += flangevin[i][2];
else flangevin[i][2] = 0;
if (ervel[i] != 0.0) erforce[i] += erforcelangevin[i];
temperature->restore_bias(i,v[i]);
}
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
erforcelangevin[i] = 0.75*gamma1*ervel[i]+0.866025404*gamma2*(random->uniform()-0.5);
f[i][0] += flangevin[i][0];
f[i][1] += flangevin[i][1];
f[i][2] += flangevin[i][2];
if (abs(spin[i])==1) erforce[i] += erforcelangevin[i];
}
}
} else {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
erforcelangevin[i] = gamma1*ervel[i]+gamma2*(random->uniform()-0.5);
f[i][0] += flangevin[i][0];
f[i][1] += flangevin[i][1];
f[i][2] += flangevin[i][2];
erforce[i] += erforcelangevin[i];
}
}
} else if (which == BIAS) {
double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
temperature->remove_bias(i,v[i]);
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
erforcelangevin[i] = gamma1*ervel[i]+gamma2*(random->uniform()-0.5);
if (v[i][0] != 0.0) f[i][0] += flangevin[i][0];
else flangevin[i][0] = 0.0;
if (v[i][1] != 0.0) f[i][1] += flangevin[i][1];
else flangevin[i][1] = 0.0;
if (v[i][2] != 0.0) f[i][2] += flangevin[i][2];
else flangevin[i][2] = 0.0;
if (ervel[i] != 0.0) erforce[i] += erforcelangevin[i];
temperature->restore_bias(i,v[i]);
}
} else if (which == BIAS) {
double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
temperature->remove_bias(i,v[i]);
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
erforcelangevin[i] = 0.75*gamma1*ervel[i]+0.866025404*gamma2*(random->uniform()-0.5);
if (v[i][0] != 0.0) f[i][0] += flangevin[i][0];
else flangevin[i][0] = 0.0;
if (v[i][1] != 0.0) f[i][1] += flangevin[i][1];
else flangevin[i][1] = 0.0;
if (v[i][2] != 0.0) f[i][2] += flangevin[i][2];
else flangevin[i][2] = 0.0;
if (abs(spin[i])==1 && ervel[i] != 0.0) erforce[i] += erforcelangevin[i];
temperature->restore_bias(i,v[i]);
}
}
}
@ -298,14 +201,16 @@ void FixLangevinEff::end_of_step()
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int *spin = atom->spin;
energy_onestep = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (mask[i] & groupbit) {
energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
flangevin[i][2]*v[i][2] + erforcelangevin[i];
flangevin[i][2]*v[i][2];
if (abs(spin[i])==1) energy_onestep += erforcelangevin[i];
}
energy += energy_onestep*update->dt;
}
@ -320,13 +225,16 @@ double FixLangevinEff::compute_scalar()
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int *spin = atom->spin;
if (update->ntimestep == update->beginstep) {
energy_onestep = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (mask[i] & groupbit) {
energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
flangevin[i][2]*v[i][2] + erforcelangevin[i];
flangevin[i][2]*v[i][2];
if (abs(spin[i])==1) energy_onestep += erforcelangevin[i];
}
energy = 0.5*energy_onestep*update->dt;
}

30
src/USER-EFF/fix_nh_eff.cpp
View File

@ -32,7 +32,7 @@ FixNHEff::FixNHEff(LAMMPS *lmp, int narg, char **arg) : FixNH(lmp, narg, arg)
{
if (!atom->spin_flag || !atom->eradius_flag ||
!atom->ervel_flag || !atom->erforce_flag)
error->all("Fix nvt/nph/npt eff requires atom attributes "
error->all("Fix nvt/nph/npt/eff requires atom attributes "
"spin, eradius, ervel, erforce");
}
@ -48,7 +48,6 @@ void FixNHEff::nve_v()
double *erforce = atom->erforce;
double *ervel = atom->ervel;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *spin = atom->spin;
int *type = atom->type;
@ -56,27 +55,14 @@ void FixNHEff::nve_v()
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
// 2 cases depending on rmass vs mass
int itype;
double dtfm;
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (spin[i]) {
dtfm = dtf / rmass[i];
ervel[i] = dtfm * erforce[i] / 0.75;
}
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (spin[i]) {
dtfm = dtf / mass[type[i]];
ervel[i] = dtfm * erforce[i] / 0.75;
}
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (abs(spin[i])==1) {
dtfm = dtf / mass[type[i]];
ervel[i] = dtfm * erforce[i] / 0.75;
}
}
}
@ -101,7 +87,7 @@ void FixNHEff::nve_x()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (spin[i]) eradius[i] += dtv * ervel[i];
if (abs(spin[i])==1) eradius[i] += dtv * ervel[i];
}
/* ----------------------------------------------------------------------
@ -122,5 +108,5 @@ void FixNHEff::nh_v_temp()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (spin[i]) ervel[i] *= factor_eta;
if (abs(spin[i])==1) ervel[i] *= factor_eta;
}

2
src/USER-EFF/fix_nph_eff.h
View File

@ -13,7 +13,7 @@
#ifdef FIX_CLASS
FixStyle(nph/sphere,FixNPHEff)
FixStyle(nph/eff,FixNPHEff)
#else

11
src/USER-EFF/fix_npt_eff.cpp
View File

@ -12,19 +12,8 @@
------------------------------------------------------------------------- */
#include "string.h"
#include "stdlib.h"
#include "math.h"
#include "fix_npt_eff.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "modify.h"
#include "fix_deform.h"
#include "compute.h"
#include "kspace.h"
#include "update.h"
#include "respa.h"
#include "domain.h"
#include "error.h"
using namespace LAMMPS_NS;

131
src/USER-EFF/fix_nve_eff.cpp
View File

@ -30,14 +30,37 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
FixNVEEff::FixNVEEff(LAMMPS *lmp, int narg, char **arg) :
FixNVE(lmp, narg, arg)
Fix(lmp, narg, arg)
{
// error check
if (!atom->spin_flag || !atom->eradius_flag ||
!atom->ervel_flag || !atom->erforce_flag)
error->all("Fix nve/eff requires atom attributes "
"spin, eradius, ervel, erforce");
"spin, eradius, ervel, erforce");
time_integrate = 1;
}
/* ---------------------------------------------------------------------- */
int FixNVEEff::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
mask |= INITIAL_INTEGRATE_RESPA;
mask |= FINAL_INTEGRATE_RESPA;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixNVEEff::init()
{
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
if (strcmp(update->integrate_style,"respa") == 0)
step_respa = ((Respa *) update->integrate)->step;
}
/* ----------------------------------------------------------------------
@ -57,7 +80,6 @@ void FixNVEEff::initial_integrate(int vflag)
double **f = atom->f;
double *erforce = atom->erforce;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *spin = atom->spin;
int *type = atom->type;
int *mask = atom->mask;
@ -66,37 +88,20 @@ void FixNVEEff::initial_integrate(int vflag)
// x + dt * [v + 0.5 * dt * (f / m)];
if (rmass) {
if (mass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
if (spin[i]) {
ervel[i] += dtfm * erforce[i] / 0.75;
eradius[i] += dtv * ervel[i];
}
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
if (spin[i]) {
ervel[i] += dtfm * erforce[i] / 0.75;
eradius[i] += dtv * ervel[i];
}
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
if (abs(spin[i])==1) {
ervel[i] += dtfm * erforce[i] / 0.75;
eradius[i] += dtv * ervel[i];
}
}
}
}
@ -113,7 +118,6 @@ void FixNVEEff::final_integrate()
double *erforce = atom->erforce;
double **f = atom->f;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *spin = atom->spin;
int *type = atom->type;
int *mask = atom->mask;
@ -122,28 +126,47 @@ void FixNVEEff::final_integrate()
// dyn_v[i] += m * dt * dyn_f[i];
if (rmass) {
if (mass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
if (spin[i] != 0)
ervel[i] += dtfm * erforce[i] / 0.75;
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
if (spin[i] != 0)
ervel[i] += dtfm * erforce[i] / 0.75;
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
if (abs(spin[i])==1)
ervel[i] += dtfm * erforce[i] / 0.75;
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixNVEEff::initial_integrate_respa(int vflag, int ilevel, int iloop)
{
dtv = step_respa[ilevel];
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
// innermost level - NVE update of v and x
// all other levels - NVE update of v
if (ilevel == 0) initial_integrate(vflag);
else final_integrate();
}
/* ---------------------------------------------------------------------- */
void FixNVEEff::final_integrate_respa(int ilevel, int iloop)
{
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
final_integrate();
}
/* ---------------------------------------------------------------------- */
void FixNVEEff::reset_dt()
{
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
}

18
src/USER-EFF/fix_nve_eff.h
View File

@ -20,15 +20,25 @@ FixStyle(nve/eff,FixNVEEff)
#ifndef LMP_FIX_NVE_EFF_H
#define LMP_FIX_NVE_EFF_H
#include "fix_nve.h"
#include "fix.h"
namespace LAMMPS_NS {
class FixNVEEff : public FixNVE {
class FixNVEEff : public Fix {
public:
FixNVEEff(class LAMMPS *, int, char **);
void initial_integrate(int);
void final_integrate();
int setmask();
virtual void init();
virtual void initial_integrate(int);
virtual void final_integrate();
void initial_integrate_respa(int, int, int);
void final_integrate_respa(int, int);
void reset_dt();
protected:
double dtv,dtf;
double *step_respa;
int mass_require;
};
}

2
src/USER-EFF/fix_nvt_eff.h
View File

@ -27,7 +27,7 @@ namespace LAMMPS_NS {
class FixNVTEff : public FixNHEff {
public:
FixNVTEff(class LAMMPS *, int, char **);
virtual ~FixNVTEff() {}
~FixNVTEff() {}
};
}

14
src/USER-EFF/fix_nvt_sllod_eff.cpp
View File

@ -12,8 +12,6 @@
------------------------------------------------------------------------- */
#include "math.h"
#include "string.h"
#include "stdlib.h"
#include "fix_nvt_sllod_eff.h"
#include "math_extra.h"
#include "atom.h"
@ -32,13 +30,17 @@ enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp
/* ---------------------------------------------------------------------- */
FixNVTSllodEff::FixNVTSllodEff(LAMMPS *lmp, int narg, char **arg) :
FixNVTEff(lmp, narg, arg)
FixNHEff(lmp, narg, arg)
{
if (!tstat_flag)
error->all("Temperature control must be used with fix nvt/sllod/eff");
if (pstat_flag)
error->all("Pressure control can not be used with fix nvt/sllod/eff");
// default values
if (mtchain_default_flag) mtchain = 1;
// create a new compute temp style
// id = fix-ID + temp
@ -46,7 +48,7 @@ FixNVTSllodEff::FixNVTSllodEff(LAMMPS *lmp, int narg, char **arg) :
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = group->names[igroup];
@ -61,7 +63,7 @@ FixNVTSllodEff::FixNVTSllodEff(LAMMPS *lmp, int narg, char **arg) :
void FixNVTSllodEff::init()
{
FixNVTEff::init();
FixNHEff::init();
if (!temperature->tempbias)
error->all("Temperature for fix nvt/sllod/eff does not have a bias");
@ -119,7 +121,7 @@ void FixNVTSllodEff::nh_v_temp()
v[i][1] = v[i][1]*factor_eta - dthalf*vdelu[1];
v[i][2] = v[i][2]*factor_eta - dthalf*vdelu[2];
temperature->restore_bias(i,v[i]);
if (spin[i])
if (abs(spin[i])==1)
ervel[i] = ervel[i]*factor_eta -
dthalf*sqrt(pow(vdelu[0],2)+pow(vdelu[1],2)+pow(vdelu[2],2));
}

4
src/USER-EFF/fix_nvt_sllod_eff.h
View File

@ -20,11 +20,11 @@ FixStyle(nvt/sllod/eff,FixNVTSllodEff)
#ifndef LMP_FIX_NVT_SLODD_EFF_H
#define LMP_FIX_NVT_SLODD_EFF_H
#include "fix_nvt_eff.h"
#include "fix_nh_eff.h"
namespace LAMMPS_NS {
class FixNVTSllodEff : public FixNVTEff {
class FixNVTSllodEff : public FixNHEff {
public:
FixNVTSllodEff(class LAMMPS *, int, char **);
~FixNVTSllodEff() {}

107
src/USER-EFF/fix_temp_rescale_eff.cpp
View File

@ -25,6 +25,7 @@
#include "update.h"
#include "domain.h"
#include "region.h"
#include "comm.h"
#include "modify.h"
#include "compute.h"
#include "error.h"
@ -36,12 +37,29 @@ enum{NOBIAS,BIAS};
/* ---------------------------------------------------------------------- */
FixTempRescaleEff::FixTempRescaleEff(LAMMPS *lmp, int narg, char **arg) :
FixTempRescale(lmp, narg, arg)
Fix(lmp, narg, arg)
{
// create a new compute temp/eff, wiping out one parent class just created
if (narg < 8) error->all("Illegal fix temp/rescale/eff command");
nevery = atoi(arg[3]);
if (nevery <= 0) error->all("Illegal fix temp/rescale/eff command");
scalar_flag = 1;
global_freq = nevery;
extscalar = 1;
t_start = atof(arg[4]);
t_stop = atof(arg[5]);
t_window = atof(arg[6]);
fraction = atof(arg[7]);
// create a new compute temp/eff
// id = fix-ID + temp, compute group = fix group
modify->delete_compute(id_temp);
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[6];
newarg[0] = id_temp;
@ -49,6 +67,42 @@ FixTempRescaleEff::FixTempRescaleEff(LAMMPS *lmp, int narg, char **arg) :
newarg[2] = (char *) "temp/eff";
modify->add_compute(3,newarg);
delete [] newarg;
tflag = 1;
energy = 0.0;
}
/* ---------------------------------------------------------------------- */
FixTempRescaleEff::~FixTempRescaleEff()
{
// delete temperature if fix created it
if (tflag) modify->delete_compute(id_temp);
delete [] id_temp;
}
/* ---------------------------------------------------------------------- */
int FixTempRescaleEff::setmask()
{
int mask = 0;
mask |= END_OF_STEP;
mask |= THERMO_ENERGY;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixTempRescaleEff::init()
{
int icompute = modify->find_compute(id_temp);
if (icompute < 0)
error->all("Temperature ID for fix temp/rescale/eff does not exist");
temperature = modify->compute[icompute];
if (temperature->tempbias) which = BIAS;
else which = NOBIAS;
}
/* ---------------------------------------------------------------------- */
@ -86,7 +140,7 @@ void FixTempRescaleEff::end_of_step()
v[i][0] *= factor;
v[i][1] *= factor;
v[i][2] *= factor;
if (spin[i])
if (abs(spin[i])==1)
ervel[i] *= factor;
}
}
@ -98,7 +152,7 @@ void FixTempRescaleEff::end_of_step()
v[i][0] *= factor;
v[i][1] *= factor;
v[i][2] *= factor;
if (spin[i])
if (abs(spin[i])==1)
ervel[i] *= factor;
temperature->restore_bias(i,v[i]);
}
@ -107,3 +161,46 @@ void FixTempRescaleEff::end_of_step()
}
}
/* ---------------------------------------------------------------------- */
int FixTempRescaleEff::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"temp/eff") == 0) {
if (narg < 2) error->all("Illegal fix_modify command");
if (tflag) {
modify->delete_compute(id_temp);
tflag = 0;
}
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
int icompute = modify->find_compute(id_temp);
if (icompute < 0) error->all("Could not find fix_modify temperature ID");
temperature = modify->compute[icompute];
if (temperature->tempflag == 0)
error->all("Fix_modify temperature ID does not compute temperature");
if (temperature->igroup != igroup && comm->me == 0)
error->warning("Group for fix_modify temp != fix group");
return 2;
}
return 0;
}
/* ---------------------------------------------------------------------- */
void FixTempRescaleEff::reset_target(double t_new)
{
t_start = t_stop = t_new;
}
/* ---------------------------------------------------------------------- */
double FixTempRescaleEff::compute_scalar()
{
return energy;
}

22
src/USER-EFF/fix_temp_rescale_eff.h
View File

@ -20,15 +20,29 @@ FixStyle(temp/rescale/eff,FixTempRescaleEff)
#ifndef LMP_FIX_TEMP_RESCALE_EFF_H
#define LMP_FIX_TEMP_RESCALE_EFF_H
#include "fix_temp_rescale.h"
#include "fix.h"
namespace LAMMPS_NS {
class FixTempRescaleEff : public FixTempRescale {
class FixTempRescaleEff : public Fix {
public:
FixTempRescaleEff(class LAMMPS *, int, char **);
~FixTempRescaleEff() {}
void end_of_step();
virtual ~FixTempRescaleEff();
int setmask();
void init();
virtual void end_of_step();
int modify_param(int, char **);
void reset_target(double);
double compute_scalar();
protected:
int which;
double t_start,t_stop,t_window;
double fraction,energy,efactor;
char *id_temp;
class Compute *temperature;
int tflag;
};
}

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

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Andres Jaramillo-Botero and Julius Su (Caltech)
Contributing author: Andres Jaramillo-Botero
------------------------------------------------------------------------- */
#include "math.h"
@ -46,12 +46,15 @@ PairEffCut::PairEffCut(LAMMPS *lmp) : Pair(lmp)
nmax = 0;
min_eradius = NULL;
min_erforce = NULL;
nextra = 4;
pvector = new double[nextra];
}
/* ---------------------------------------------------------------------- */
PairEffCut::~PairEffCut()
{
delete [] pvector;
memory->sfree(min_eradius);
memory->sfree(min_erforce);
@ -67,12 +70,19 @@ PairEffCut::~PairEffCut()
void PairEffCut::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,energy;
double fpair,fx,fy,fz,e1rforce,e2rforce,e1rvirial,e2rvirial;
double rsq,rc,forcecoul,factor_coul;
double xtmp,ytmp,ztmp,delx,dely,delz,energy;
double eke,ecoul,epauli,errestrain,halfcoul,halfpauli;
double fpair,fx,fy,fz;
double e1rforce,e2rforce,e1rvirial,e2rvirial;
double s_fpair, s_e1rforce, s_e2rforce;
double ecp_epauli, ecp_fpair, ecp_e1rforce, ecp_e2rforce;
double rsq,rc;
int *ilist,*jlist,*numneigh,**firstneigh;
ecoul = 0.0;
energy = eke = epauli = ecoul = errestrain = 0.0;
// pvector = [KE, Pauli, ecoul, radial_restraint]
for (i=0; i<4; i++) pvector[i] = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
@ -84,12 +94,10 @@ void PairEffCut::compute(int eflag, int vflag)
int *spin = atom->spin;
int *type = atom->type;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
double *special_coul = force->special_coul;
int newton_pair = force->newton_pair;
double qqrd2e = force->qqrd2e;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
@ -99,7 +107,6 @@ void PairEffCut::compute(int eflag, int vflag)
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
qtmp = q[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
@ -107,222 +114,415 @@ void PairEffCut::compute(int eflag, int vflag)
jlist = firstneigh[i];
jnum = numneigh[i];
// add electron kinetic energy
// add electron wavefuntion kinetic energy (not pairwise)
if (spin[i] != 0) {
e1rforce = energy = 0;
energy = 1.5 / (eradius[i] * eradius[i]);
e1rforce = 3.0 / (eradius[i] * eradius[i] * eradius[i]);
if (abs(spin[i])==1 || spin[i]==2) {
// reset energy and force temp variables
eke = epauli = ecoul = 0.0;
fpair = e1rforce = e2rforce = 0.0;
s_fpair = 0.0;
KinElec(eradius[i],&eke,&e1rforce);
// Fixed-core
if (spin[i] == 2) {
// KE(2s)+Coul(1s-1s)+Coul(2s-nuclei)+Pauli(2s)
eke *= 2;
ElecNucElec(q[i],0.0,eradius[i],&ecoul,&fpair,&e1rforce);
ElecNucElec(q[i],0.0,eradius[i],&ecoul,&fpair,&e1rforce);
ElecElecElec(0.0,eradius[i],eradius[i],&ecoul,&fpair,&e1rforce,&e2rforce);
// opposite spin electron interactions
PauliElecElec(0,0.0,eradius[i],eradius[i],
&epauli,&s_fpair,&e1rforce,&e2rforce);
// fix core electron size, i.e. don't contribute to ervirial
e2rforce = e1rforce = 0.0;
}
// apply unit conversion factors
eke *= hhmss2e;
ecoul *= qqrd2e;
fpair *= qqrd2e;
epauli *= hhmss2e;
s_fpair *= hhmss2e;
e1rforce *= hhmss2e;
// Sum up contributions
energy = eke + epauli + ecoul;
fpair = fpair + s_fpair;
erforce[i] += e1rforce;
// electronic ke accumulates into ecoul (pot)
if (eflag) ecoul = energy; // KE e-wavefunction
// Tally energy and compute radial atomic virial contribution
if (evflag) {
ev_tally_eff(i,i,nlocal,newton_pair,ecoul,0.0);
if (flexible_pressure_flag) // only on electron
ev_tally_eff(i,i,nlocal,newton_pair,energy,0.0);
if (flexible_pressure_flag) // iff flexible pressure flag on
ev_tally_eff(i,i,nlocal,newton_pair,0.0,e1rforce*eradius[i]);
}
if (eflag_global) {
pvector[0] += eke;
pvector[1] += epauli;
pvector[2] += ecoul;
}
}
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
rc = sqrt(rsq);
if (j < nall) factor_coul = 1.0;
else {
factor_coul = special_coul[j/nall];
j %= nall;
}
jtype = type[j];
double taper = sqrt(cutsq[itype][jtype]);
if (rsq < cutsq[itype][jtype]) {
// nuclei-nuclei interaction
energy = ecoul = epauli = ecp_epauli = 0.0;
fx = fy = fz = fpair = s_fpair = ecp_fpair = 0.0;
double taper = sqrt(cutsq[itype][jtype]);
double dist = rc / taper;
double spline = cutoff(dist);
double dspline = dcutoff(dist) / taper;
// nucleus (i) - nucleus (j) Coul interaction
if (spin[i] == 0 && spin[j] == 0) {
energy = fx = fy = fz = 0;
double qxq = qqrd2e*qtmp*q[j];
forcecoul = qxq/rsq;
double dist = rc / taper;
double spline = cutoff(dist);
double dspline = dcutoff(dist) / taper;
energy = factor_coul*qxq/rc;
fpair = forcecoul*spline-energy*dspline;
fpair = qqrd2e*fpair/rc;
energy = spline*energy;
fx = delx*fpair;
fy = dely*fpair;
fz = delz*fpair;
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
if (newton_pair || j < nlocal) {
f[j][0] -= fx;
f[j][1] -= fy;
f[j][2] -= fz;
}
if (eflag) ecoul = energy; // Electrostatics:N-N
if (evflag)
ev_tally_xyz(i,j,nlocal,newton_pair,0.0,ecoul,
fx,fy,fz,delx,dely,delz);
}
// I is nucleus, J is electron
double qxq = q[i]*q[j];
if (spin[i] == 0 && spin[j] != 0) {
energy = fpair = e1rforce = fx = fy = fz = e1rvirial = 0;
ElecNucElec(-q[i],rc,eradius[j],&energy,&fpair,&e1rforce,i,j);
double dist = rc / taper;
double spline = cutoff(dist);
double dspline = dcutoff(dist) / taper;
fpair = qqrd2e * (fpair * spline - energy * dspline);
energy = qqrd2e * spline * energy;
e1rforce = qqrd2e * spline * e1rforce;
ElecNucNuc(qxq, rc, &ecoul, &fpair);
}
// fixed-core (i) - nucleus (j) nuclear Coul interaction
else if (spin[i] == 2 && spin[j] == 0) {
double qxq = q[i]*q[j];
e1rforce = 0.0;
ElecNucNuc(qxq, rc, &ecoul, &fpair);
ElecNucElec(q[j],rc,eradius[i],&ecoul,&fpair,&e1rforce);
ElecNucElec(q[j],rc,eradius[i],&ecoul,&fpair,&e1rforce);
}
// nucleus (i) - fixed-core (j) nuclear Coul interaction
else if (spin[i] == 0 && spin[j] == 2) {
double qxq = q[i]*q[j];
e1rforce = 0.0;
ElecNucNuc(qxq, rc, &ecoul, &fpair);
ElecNucElec(q[i],rc,eradius[j],&ecoul,&fpair,&e1rforce);
ElecNucElec(q[i],rc,eradius[j],&ecoul,&fpair,&e1rforce);
}
// pseudo-core nucleus (i) - nucleus (j) interaction
else if (spin[i] == 3 && 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];
ElecCoreNuc(qxq, rc, eradius[j], &ecoul, &fpair);
}
// nucleus (i) - electron (j) Coul interaction
else if (spin[i] == 0 && abs(spin[j]) == 1) {
e1rforce = 0.0;
ElecNucElec(q[i],rc,eradius[j],&ecoul,&fpair,&e1rforce);
e1rforce = spline * qqrd2e * e1rforce;
erforce[j] += e1rforce;
e1rvirial = eradius[j] * e1rforce;
SmallRForce(delx,dely,delz,rc,fpair,&fx,&fy,&fz);
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
if (newton_pair || j < nlocal) {
f[j][0] -= fx;
f[j][1] -= fy;
f[j][2] -= fz;
}
if (eflag) ecoul = energy; // Electrostatics:N-e
if (evflag) {
ev_tally_xyz(i,j,nlocal,newton_pair,0.0,ecoul,
fx,fy,fz,delx,dely,delz);
if (flexible_pressure_flag) // only on electron
ev_tally_eff(j,j,nlocal,newton_pair,0.0,e1rvirial);
// Radial electron virial, iff flexible pressure flag set
if (evflag && flexible_pressure_flag) {
e1rvirial = eradius[j] * e1rforce;
ev_tally_eff(j,j,nlocal,newton_pair,0.0,e1rvirial);
}
}
// I is electon, J is nucleus
// electron (i) - nucleus (j) Coul interaction
if (spin[i] != 0 && spin[j] == 0) {
energy = fpair = e1rforce = fx = fy = fz = e1rvirial = 0;
ElecNucElec(-q[j],rc,eradius[i],&energy,&fpair,&e1rforce,j,i);
double dist = rc / taper;
double spline = cutoff(dist);
double dspline = dcutoff(dist) / taper;
fpair = qqrd2e * (fpair * spline - energy * dspline);
energy = qqrd2e * spline * energy;
e1rforce = qqrd2e * spline * e1rforce;
else if (abs(spin[i]) == 1 && spin[j] == 0) {
e1rforce = 0.0;
ElecNucElec(q[j],rc,eradius[i],&ecoul,&fpair,&e1rforce);
e1rforce = spline * qqrd2e * e1rforce;
erforce[i] += e1rforce;
e1rvirial = eradius[i] * e1rforce;
SmallRForce(delx,dely,delz,rc,fpair,&fx,&fy,&fz);
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
if (newton_pair || j < nlocal) {
f[j][0] -= fx;
f[j][1] -= fy;
f[j][2] -= fz;
}
if (eflag) ecoul = energy; //Electrostatics-e-N
if (evflag) {
ev_tally_xyz(i,j,nlocal,newton_pair,0.0,ecoul,
fx,fy,fz,delx,dely,delz);
if (flexible_pressure_flag) // only on electron
ev_tally_eff(i,i,nlocal,newton_pair,0.0,e1rvirial);
// Radial electron virial, iff flexible pressure flag set
if (evflag && flexible_pressure_flag) {
e1rvirial = eradius[i] * e1rforce;
ev_tally_eff(i,i,nlocal,newton_pair,0.0,e1rvirial);
}
}
// electron-electron interaction
if (spin[i] && spin[j]) {
energy = fpair = fx = fy= fz =
e1rforce = e2rforce = e1rvirial = e2rvirial = 0.0;
ElecElecElec(rc,eradius[i],eradius[j],&energy,&fpair,
&e1rforce,&e2rforce,i,j);
double s_energy, s_fpair, s_e1rforce, s_e2rforce;
s_energy = s_fpair = s_e1rforce = s_e2rforce = 0.0;
// as with the electron ke,
// the Pauli term is also accumulated into ecoul (pot)
// electron (i) - electron (j) interactions
PauliElecElec(spin[j] == spin[i],rc,eradius[i],eradius[j],
&s_energy,&s_fpair,&s_e1rforce,&s_e2rforce,i,j);
double dist = rc / taper;
double spline = cutoff(dist);
double dspline = dcutoff(dist) / taper;
// apply spline cutoff
else if (abs(spin[i]) == 1 && abs(spin[j]) == 1) {
e1rforce = e2rforce = 0.0;
s_e1rforce = s_e2rforce = 0.0;
s_fpair = qqrd2e * (s_fpair * spline - s_energy * dspline);
s_energy = qqrd2e * spline * s_energy;
fpair = qqrd2e * (fpair * spline - energy * dspline);
energy = qqrd2e * spline * energy;
e1rforce = qqrd2e * spline * (e1rforce + s_e1rforce);
e2rforce = qqrd2e * spline * (e2rforce + s_e2rforce);
// Cartesian and radial forces
ElecElecElec(rc,eradius[i],eradius[j],&ecoul,&fpair,
&e1rforce,&e2rforce);
PauliElecElec(spin[i] == spin[j],rc,eradius[i],eradius[j],
&epauli,&s_fpair,&s_e1rforce,&s_e2rforce);
SmallRForce(delx, dely, delz, rc, fpair + s_fpair, &fx, &fy, &fz);
erforce[i] += e1rforce;
erforce[j] += e2rforce;
// Apply conversion factor
epauli *= hhmss2e;
s_fpair *= hhmss2e;
// radial virials
e1rforce = spline * (qqrd2e * e1rforce + hhmss2e * s_e1rforce);
erforce[i] += e1rforce;
e2rforce = spline * (qqrd2e * e2rforce + hhmss2e * s_e2rforce);
erforce[j] += e2rforce;
e1rvirial = eradius[i] * e1rforce;
e2rvirial = eradius[j] * e2rforce;
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
if (newton_pair || j < nlocal) {
f[j][0] -= fx;
f[j][1] -= fy;
f[j][2] -= fz;
}
if (eflag) ecoul = energy + s_energy; // Electrostatics+Pauli: e-e
if (evflag) {
ev_tally_xyz(i,j,nlocal,newton_pair,0.0,
ecoul,fx,fy,fz,delx,dely,delz);
if (flexible_pressure_flag) // on both electrons
ev_tally_eff(i,j,nlocal,newton_pair,0.0,e1rvirial+e2rvirial);
// Radial electron virial, iff flexible pressure flag set
if (evflag && flexible_pressure_flag) {
e1rvirial = eradius[i] * e1rforce;
e2rvirial = eradius[j] * e2rforce;
ev_tally_eff(i,j,nlocal,newton_pair,0.0,e1rvirial+e2rvirial);
}
}
}
// fixed-core (i) - electron (j) interactions
else if (spin[i] == 2 && abs(spin[j]) == 1) {
e1rforce = e2rforce = 0.0;
s_e1rforce = s_e2rforce = 0.0;
ElecNucElec(q[i],rc,eradius[j],&ecoul,&fpair,&e2rforce);
ElecElecElec(rc,eradius[i],eradius[j],&ecoul,&fpair,
&e1rforce,&e2rforce);
ElecElecElec(rc,eradius[i],eradius[j],&ecoul,&fpair,
&e1rforce,&e2rforce);
PauliElecElec(0,rc,eradius[i],eradius[j],&epauli,
&s_fpair,&s_e1rforce,&s_e2rforce);
PauliElecElec(1,rc,eradius[i],eradius[j],&epauli,
&s_fpair,&s_e1rforce,&s_e2rforce);
// Apply conversion factor
epauli *= hhmss2e;
s_fpair *= hhmss2e;
// only update virial for j electron
e2rforce = spline * (qqrd2e * e2rforce + hhmss2e * s_e2rforce);
erforce[j] += e2rforce;
// Radial electron virial, iff flexible pressure flag set
if (evflag && flexible_pressure_flag) {
e2rvirial = eradius[j] * e2rforce;
ev_tally_eff(j,j,nlocal,newton_pair,0.0,e2rvirial);
}
}
// electron (i) - fixed-core (j) interactions
else if (abs(spin[i]) == 1 && spin[j] == 2) {
e1rforce = e2rforce = 0.0;
s_e1rforce = s_e2rforce = 0.0;
ElecNucElec(q[j],rc,eradius[i],&ecoul,&fpair,&e2rforce);
ElecElecElec(rc,eradius[j],eradius[i],&ecoul,&fpair,
&e1rforce,&e2rforce);
ElecElecElec(rc,eradius[j],eradius[i],&ecoul,&fpair,
&e1rforce,&e2rforce);
PauliElecElec(0,rc,eradius[j],eradius[i],&epauli,
&s_fpair,&s_e1rforce,&s_e2rforce);
PauliElecElec(1,rc,eradius[j],eradius[i],&epauli,
&s_fpair,&s_e1rforce,&s_e2rforce);
// Apply conversion factor
epauli *= hhmss2e;
s_fpair *= hhmss2e;
// only update virial for i electron
e2rforce = spline * (qqrd2e * e2rforce + hhmss2e * s_e2rforce);
erforce[i] += e2rforce;
// add radial atomic virial, iff flexible pressure flag set
if (evflag && flexible_pressure_flag) {
e2rvirial = eradius[i] * e2rforce;
ev_tally_eff(i,i,nlocal,newton_pair,0.0,e2rvirial);
}
}
// fixed-core (i) - fixed-core (j) interactions
else if (spin[i] == 2 && spin[j] == 2) {
e1rforce = e2rforce = 0.0;
s_e1rforce = s_e2rforce = 0.0;
double qxq = q[i]*q[j];
ElecNucNuc(qxq, rc, &ecoul, &fpair);
ElecNucElec(q[i],rc,eradius[j],&ecoul,&fpair,&e1rforce);
ElecNucElec(q[i],rc,eradius[j],&ecoul,&fpair,&e1rforce);
ElecNucElec(q[j],rc,eradius[i],&ecoul,&fpair,&e1rforce);
ElecNucElec(q[j],rc,eradius[i],&ecoul,&fpair,&e1rforce);
ElecElecElec(rc,eradius[i],eradius[j],&ecoul,&fpair,
&e1rforce,&e2rforce);
ElecElecElec(rc,eradius[i],eradius[j],&ecoul,&fpair,
&e1rforce,&e2rforce);
ElecElecElec(rc,eradius[i],eradius[j],&ecoul,&fpair,
&e1rforce,&e2rforce);
ElecElecElec(rc,eradius[i],eradius[j],&ecoul,&fpair,
&e1rforce,&e2rforce);
PauliElecElec(0,rc,eradius[i],eradius[j],&epauli,
&s_fpair,&s_e1rforce,&s_e2rforce);
PauliElecElec(1,rc,eradius[i],eradius[j],&epauli,
&s_fpair,&s_e1rforce,&s_e2rforce);
epauli *= 2;
s_fpair *= 2;
// Apply conversion factor
epauli *= hhmss2e;
s_fpair *= hhmss2e;
}
// pseudo-core (i) - electron/fixed-core electrons (j) interactions
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);
PauliCoreElec(rc,eradius[j],&ecp_epauli,&ecp_fpair,
&ecp_e2rforce,PAULI_CORE_A, PAULI_CORE_B,
PAULI_CORE_C);
}
// Apply conversion factor from Hartree to kcal/mol
ecp_epauli *= h2e;
ecp_fpair *= h2e;
// only update virial for j electron
e2rforce = spline * (qqrd2e * e2rforce + h2e * ecp_e2rforce);
erforce[j] += e2rforce;
// add radial atomic virial, iff flexible pressure flag set
if (evflag && flexible_pressure_flag) {
e2rvirial = eradius[j] * e2rforce;
ev_tally_eff(j,j,nlocal,newton_pair,0.0,e2rvirial);
}
}
// electron/fixed-core electrons (i) - pseudo-core (j) interactions
else if ((abs(spin[i]) == 1 || spin[i] == 2) && spin[j] == 3) {
e1rforce = ecp_e1rforce = 0.0;
if (abs(spin[j]) == 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);
} 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);
PauliCoreElec(rc,eradius[i],&ecp_epauli,&ecp_fpair,
&ecp_e1rforce,PAULI_CORE_A, PAULI_CORE_B,
PAULI_CORE_C);
}
// Apply conversion factor from Hartree to kcal/mol
ecp_epauli *= h2e;
ecp_fpair *= h2e;
// only update virial for j electron
e1rforce = spline * (qqrd2e * e1rforce + h2e * ecp_e1rforce);
erforce[i] += e1rforce;
// add radial atomic virial, iff flexible pressure flag set
if (evflag && flexible_pressure_flag) {
e1rvirial = eradius[i] * e1rforce;
ev_tally_eff(i,i,nlocal,newton_pair,0.0,e1rvirial);
}
}
// pseudo-core (i) - pseudo-core (j) interactions
else if (spin[i] == 3 && abs(spin[j]) == 3) {
double qxq = q[i]*q[j];
ElecCoreCore(qxq,rc,eradius[i],eradius[j],&ecoul,&fpair);
}
// Apply Coulomb conversion factor for all cases
ecoul *= qqrd2e;
fpair *= qqrd2e;
// Sum up energy and force contributions
epauli += ecp_epauli;
energy = ecoul + epauli;
fpair = fpair + s_fpair + ecp_fpair;
// Apply cutoff spline
fpair = fpair * spline - energy * dspline;
energy = spline * energy;
// Tally cartesian forces
SmallRForce(delx,dely,delz,rc,fpair,&fx,&fy,&fz);
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
if (newton_pair || j < nlocal) {
f[j][0] -= fx;
f[j][1] -= fy;
f[j][2] -= fz;
}
// Tally energy (in ecoul) and compute normal pressure virials
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,0.0,
energy,fx,fy,fz,delx,dely,delz);
if (eflag_global) {
if (newton_pair) {
pvector[1] += spline * epauli;
pvector[2] += spline * ecoul;
}
else {
halfpauli = 0.5 * spline * epauli;
halfcoul = 0.5 * spline * ecoul;
if (i < nlocal) {
pvector[1] += halfpauli;
pvector[2] += halfcoul;
}
if (j < nlocal) {
pvector[1] += halfpauli;
pvector[2] += halfcoul;
}
}
}
}
}
// limit the electron size for periodic systems, to max=half-box-size
// limit_size_stiffness for electrons
// limit electron stifness (size) for periodic systems, to max=half-box-size
if (abs(spin[i]) == 1 && limit_size_flag) {
double half_box_length=0, dr, kfactor=hhmss2e*1.0;
e1rforce = errestrain = 0.0;
if (spin[i] && limit_size_flag) {
double half_box_length=0, dr, k=1.0;
e1rforce = energy = 0;
if (domain->xperiodic == 1 || domain->yperiodic == 1 ||
domain->zperiodic == 1) {
delx = domain->boxhi[0]-domain->boxlo[0];
@ -331,24 +531,23 @@ void PairEffCut::compute(int eflag, int vflag)
half_box_length = 0.5 * MIN(delx, MIN(dely, delz));
if (eradius[i] > half_box_length) {
dr = eradius[i]-half_box_length;
energy=0.5*k*dr*dr;
e1rforce=-k*dr;
}
}
erforce[i] += e1rforce;
errestrain=0.5*kfactor*dr*dr;
e1rforce=-kfactor*dr;
if (eflag_global) pvector[3] += errestrain;
// constraint radial energy accumulated as ecoul
erforce[i] += e1rforce;
if (eflag) ecoul = energy; // Radial constraint energy
if (evflag) {
ev_tally_eff(i,i,nlocal,newton_pair,ecoul,0.0);
if (flexible_pressure_flag) // only on electron
ev_tally_eff(i,i,nlocal,newton_pair,0.0,eradius[i]*e1rforce);
// 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
ev_tally_eff(i,i,nlocal,newton_pair,0.0,eradius[i]*e1rforce);
}
}
}
}
}
}
if (vflag_fdotr) {
virial_compute();
if (flexible_pressure_flag) virial_eff_compute();
@ -411,56 +610,58 @@ void PairEffCut::virial_eff_compute()
------------------------------------------------------------------------- */
void PairEffCut::ev_tally_eff(int i, int j, int nlocal, int newton_pair,
double ecoul, double e_virial)
double energy, double e_virial)
{
double ecoulhalf,epairhalf;
double energyhalf;
double partial_evirial = e_virial/3.0;
double half_partial_evirial = partial_evirial/2;
int *spin = atom->spin;
// accumulate electronic wavefunction ke and radial constraint as ecoul
if (eflag_either) {
if (eflag_global) {
ecoulhalf = 0.5*ecoul;
if (i < nlocal)
eng_coul += ecoulhalf;
if (j < nlocal)
eng_coul += ecoulhalf;
if (newton_pair)
eng_coul += energy;
else {
energyhalf = 0.5*energy;
if (i < nlocal)
eng_coul += energyhalf;
if (j < nlocal)
eng_coul += energyhalf;
}
}
if (eflag_atom) {
epairhalf = 0.5 * ecoul;
if (i < nlocal) eatom[i] += epairhalf;
if (j < nlocal) eatom[j] += epairhalf;
if (newton_pair || i < nlocal) eatom[i] += 0.5 * energy;
if (newton_pair || j < nlocal) eatom[j] += 0.5 * energy;
}
}
if (vflag_either) {
if (vflag_global) {
if (spin[i] && i < nlocal) {
virial[0] += 0.5*partial_evirial;
virial[1] += 0.5*partial_evirial;
virial[2] += 0.5*partial_evirial;
virial[0] += half_partial_evirial;
virial[1] += half_partial_evirial;
virial[2] += half_partial_evirial;
}
if (spin[j] && j < nlocal) {
virial[0] += 0.5*partial_evirial;
virial[1] += 0.5*partial_evirial;
virial[2] += 0.5*partial_evirial;
virial[0] += half_partial_evirial;
virial[1] += half_partial_evirial;
virial[2] += half_partial_evirial;
}
}
if (vflag_atom) {
if (spin[i]) {
if (newton_pair || i < nlocal) {
vatom[i][0] += 0.5*partial_evirial;
vatom[i][1] += 0.5*partial_evirial;
vatom[i][2] += 0.5*partial_evirial;
vatom[i][0] += half_partial_evirial;
vatom[i][1] += half_partial_evirial;
vatom[i][2] += half_partial_evirial;
}
}
if (spin[j]) {
if (newton_pair || j < nlocal) {
vatom[j][0] += 0.5*partial_evirial;
vatom[j][1] += 0.5*partial_evirial;
vatom[j][2] += 0.5*partial_evirial;
vatom[j][0] += half_partial_evirial;
vatom[j][1] += half_partial_evirial;
vatom[j][2] += half_partial_evirial;
}
}
}
@ -491,7 +692,13 @@ void PairEffCut::allocate()
void PairEffCut::settings(int narg, char **arg)
{
if (narg != 1 && narg != 3) error->all("Illegal pair_style command");
if (narg != 1 && narg != 3 && narg != 4 && narg != 7)
error->all("Illegal pair_style command");
// Defaults ECP parameters for Si
PAULI_CORE_A = 0.320852;
PAULI_CORE_B = 2.283269;
PAULI_CORE_C = 0.814857;
if (narg == 1) {
cut_global = force->numeric(arg[0]);
@ -501,7 +708,38 @@ void PairEffCut::settings(int narg, char **arg)
cut_global = force->numeric(arg[0]);
limit_size_flag = force->inumeric(arg[1]);
flexible_pressure_flag = force->inumeric(arg[2]);
} else if (narg == 4) {
cut_global = force->numeric(arg[0]);
limit_size_flag = 0;
flexible_pressure_flag = 0;
if (strcmp(arg[1],"ecp") != 0)
error->all("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]);
}
} 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("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]);
}
}
// 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
hhmss2e = 175.72044219620075; // hartree->kcal/mol * (Bohr->Angstrom)^2
} else if (force->qqr2e==1.0) { // electron units
h2e = 1.0;
hhmss2e = 1.0;
} else error->all("Check your units");
// reset cutoffs that have been explicitly set
@ -558,9 +796,16 @@ void PairEffCut::init_style()
// add hook to minimizer for eradius and erforce
if (update->whichflag == 2)
if (update->whichflag == 2)
int ignore = update->minimize->request(this,1,0.01);
// make sure to use the appropriate timestep when using real units
if (update->whichflag == 1) {
if (force->qqr2e == 332.06371 && update->dt == 1.0)
error->all("You must lower the default real units timestep for pEFF ");
}
// need a half neigh list and optionally a granular history neigh list
int irequest = neighbor->request(this);

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

@ -48,6 +48,8 @@ class PairEffCut : public Pair {
int limit_size_flag, flexible_pressure_flag;
double cut_global;
double **cut;
double PAULI_CORE_A, PAULI_CORE_B, PAULI_CORE_C;
double hhmss2e, h2e;
int nmax;
double *min_eradius,*min_erforce;

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

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Andres Jaramillo-Botero and Julius Su (Caltech)
Contributing authors: Andres Jaramillo-Botero, Hai Xiao, Julius Su (Caltech)
------------------------------------------------------------------------- */
namespace LAMMPS_NS {
@ -257,17 +257,24 @@ inline double ierfoverx1(double x, double *df)
/* ---------------------------------------------------------------------- */
inline void ElecNucNuc(double q, double rc, double *energy, double *frc)
inline void KinElec(double radius, double *eke, double *frc)
{
*energy += q / rc;
*eke += 1.5 / (radius * radius);
*frc += 3.0 / (radius * radius * radius);
}
/* ---------------------------------------------------------------------- */
inline void ElecNucNuc(double q, double rc, double *ecoul, double *frc)
{
*ecoul += q / rc;
*frc += q / (rc * rc);
}
/* ---------------------------------------------------------------------- */
inline void ElecNucElec(double q, double rc, double re1,
double *energy, double *frc, double *fre1,
int i, int j)
double *ecoul, double *frc, double *fre1)
{
double a, arc;
double coeff_a;
@ -284,11 +291,11 @@ inline void ElecNucElec(double q, double rc, double re1,
double E, dEdr, dEdr1, f, df;
f = ierfoverx1(arc, &df);
dEdr = -q * a * a * df;
dEdr1 = q * (a / re1) * (f + arc * df);
E = q * a * f;
dEdr = q * a * a * df;
dEdr1 = -q * (a / re1) * (f + arc * df);
E = -q * a * f;
*energy += E;
*ecoul += E;
*frc += dEdr;
*fre1 += dEdr1;
}
@ -296,8 +303,8 @@ inline void ElecNucElec(double q, double rc, double re1,
/* ---------------------------------------------------------------------- */
inline void ElecElecElec(double rc, double re1, double re2,
double *energy, double *frc, double *fre1,
double *fre2, int i, int j)
double *ecoul, double *frc, double *fre1,
double *fre2)
{
double a, arc, re, fre;
double coeff_a;
@ -329,7 +336,7 @@ inline void ElecElecElec(double rc, double re1, double re2,
E = a * f;
*energy += E;
*ecoul += E;
*frc += dEdr;
*fre1 += dEdr1;
*fre2 += dEdr2;
@ -337,10 +344,88 @@ inline void ElecElecElec(double rc, double re1, double re2,
/* ---------------------------------------------------------------------- */
inline void ElecCoreNuc(double q, double rc, double re1, double *ecoul, double *frc)
{
double a, arc;
double coeff_a;
double E, dEdr, df, f;
coeff_a = 1.4142135623730951; /* sqrt(2) */
a = coeff_a / re1;
arc = a * rc;
f = ierfoverx1(arc, &df);
dEdr = -q * a * a * df;
E = q * a * f;
*ecoul += E;
*frc += dEdr;
}
/* ---------------------------------------------------------------------- */
inline void ElecCoreCore(double q, double rc, double re1, double re2,
double *ecoul, double *frc)
{
double a, arc, re;
double coeff_a;
double E, dEdr, f, fre, df;
coeff_a = 1.4142135623730951;
re = sqrt(re1 * re1 + re2 * re2);
a = coeff_a / re;
arc = a * rc;
f = ierfoverx1(arc, &df);
dEdr = -q * a * a * df;
fre = q * a * (f + arc * df) / (re * re);
E = q * a * f;
*ecoul += E;
*frc += dEdr;
}
/* ---------------------------------------------------------------------- */
inline void ElecCoreElec(double q, double rc, double re1, double re2,
double *ecoul, double *frc, double *fre2)
{
double a,arc, re;
double coeff_a;
double E, dEdr, dEdr2, f, df, fre;
/* sqrt(2) */
coeff_a = 1.4142135623730951;
/*
re1: core size
re2: electron size
re3: size of the core, obtained from the electron density function rho(r) of core
e.g. rho(r) = a1*exp(-((r)/b1)^2), a1 =157.9, b1 = 0.1441 -> re3 = 0.1441 for Si4+
*/
re = sqrt(re1 * re1 + re2 * re2);
a = coeff_a / re;
arc = a * rc;
f = ierfoverx1(arc, &df);
E = -q * a * f;
dEdr = -q * a * df * a;
fre = q * a * (f + arc * df) / (re * re);
dEdr2 = fre * re2;
*ecoul += E;
*frc -= dEdr;
*fre2 -= dEdr2;
}
/* ---------------------------------------------------------------------- */
inline void PauliElecElec(int samespin, double rc,
double re1, double re2, double *energy,
double *frc, double *fre1, double *fre2,
int i, int j)
double re1, double re2, double *epauli,
double *frc, double *fre1, double *fre2)
{
double ree, rem;
double S, t1, t2, tt;
@ -381,7 +466,30 @@ inline void PauliElecElec(int samespin, double rc,
*fre1 -= PAULI_RE * (dTdr1 * O + ratio * dSdr1);
*fre2 -= PAULI_RE * (dTdr2 * O + ratio * dSdr2);
*frc -= PAULI_RC * (dTdr * O + ratio * dSdr);
*energy += tt * O;
*epauli += tt*O;
}
/* ---------------------------------------------------------------------- */
inline void PauliCoreElec(double rc, double re2, double *epauli, double *frc,
double *fre2, double PAULI_CORE_A, double PAULI_CORE_B, double PAULI_CORE_C)
{
double E, dEdrc, dEdre2, rcsq, ssq;
rcsq = rc * rc;
ssq = re2 * re2;
E = PAULI_CORE_A * exp((-PAULI_CORE_B * rcsq) / (ssq + PAULI_CORE_C));
dEdrc = -2 * PAULI_CORE_A * PAULI_CORE_B * rc * exp(-PAULI_CORE_B * rcsq /
(ssq + PAULI_CORE_C)) / (ssq + PAULI_CORE_C);
dEdre2 = 2 * PAULI_CORE_A * PAULI_CORE_B * re2 * rcsq * exp(-PAULI_CORE_B * rcsq /
(ssq + PAULI_CORE_C)) / pow((PAULI_CORE_C + ssq),2);
*epauli += E;
*frc -= dEdrc;
*fre2 -= dEdre2;
}
/* ---------------------------------------------------------------------- */