From 691894ab343ca7f3c3461b893f68fbd61f16d542 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 22 Feb 2011 22:36:09 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5716 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-EFF/Install.sh | 4 +- src/USER-EFF/README | 21 + src/USER-EFF/atom_vec_electron.cpp | 9 +- src/USER-EFF/compute_ke_atom_eff.cpp | 21 +- src/USER-EFF/compute_ke_eff.cpp | 16 +- src/USER-EFF/compute_temp_deform_eff.cpp | 180 +++++- src/USER-EFF/compute_temp_deform_eff.h | 27 +- src/USER-EFF/compute_temp_eff.cpp | 53 +- src/USER-EFF/compute_temp_eff.h | 11 +- src/USER-EFF/compute_temp_region_eff.cpp | 179 +++++- src/USER-EFF/compute_temp_region_eff.h | 22 +- src/USER-EFF/fix_langevin_eff.cpp | 224 +++----- src/USER-EFF/fix_nh_eff.cpp | 30 +- src/USER-EFF/fix_nph_eff.h | 2 +- src/USER-EFF/fix_npt_eff.cpp | 11 - src/USER-EFF/fix_nve_eff.cpp | 131 +++-- src/USER-EFF/fix_nve_eff.h | 18 +- src/USER-EFF/fix_nvt_eff.h | 2 +- src/USER-EFF/fix_nvt_sllod_eff.cpp | 14 +- src/USER-EFF/fix_nvt_sllod_eff.h | 4 +- src/USER-EFF/fix_temp_rescale_eff.cpp | 107 +++- src/USER-EFF/fix_temp_rescale_eff.h | 22 +- src/USER-EFF/pair_eff_cut.cpp | 693 +++++++++++++++-------- src/USER-EFF/pair_eff_cut.h | 2 + src/USER-EFF/pair_eff_inline.h | 140 ++++- 25 files changed, 1319 insertions(+), 624 deletions(-) diff --git a/src/USER-EFF/Install.sh b/src/USER-EFF/Install.sh index 64a7b71c0b..407168c564 100644 --- a/src/USER-EFF/Install.sh +++ b/src/USER-EFF/Install.sh @@ -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 diff --git a/src/USER-EFF/README b/src/USER-EFF/README index 2e500bafcc..c1c2453ecc 100644 --- a/src/USER-EFF/README +++ b/src/USER-EFF/README @@ -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. + + diff --git a/src/USER-EFF/atom_vec_electron.cpp b/src/USER-EFF/atom_vec_electron.cpp index 20e33a0919..f3a3997684 100644 --- a/src/USER-EFF/atom_vec_electron.cpp +++ b/src/USER-EFF/atom_vec_electron.cpp @@ -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; } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-EFF/compute_ke_atom_eff.cpp b/src/USER-EFF/compute_ke_atom_eff.cpp index f9dba6e5db..73e63313af 100644 --- a/src/USER-EFF/compute_ke_atom_eff.cpp +++ b/src/USER-EFF/compute_ke_atom_eff.cpp @@ -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; } diff --git a/src/USER-EFF/compute_ke_eff.cpp b/src/USER-EFF/compute_ke_eff.cpp index 600f0ebae1..08759c7379 100644 --- a/src/USER-EFF/compute_ke_eff.cpp +++ b/src/USER-EFF/compute_ke_eff.cpp @@ -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]; } } diff --git a/src/USER-EFF/compute_temp_deform_eff.cpp b/src/USER-EFF/compute_temp_deform_eff.cpp index ab6b51ea4c..5a093b65e0 100644 --- a/src/USER-EFF/compute_temp_deform_eff.cpp +++ b/src/USER-EFF/compute_temp_deform_eff.cpp @@ -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; +} + diff --git a/src/USER-EFF/compute_temp_deform_eff.h b/src/USER-EFF/compute_temp_deform_eff.h index 32bf165c47..47f4409422 100644 --- a/src/USER-EFF/compute_temp_deform_eff.h +++ b/src/USER-EFF/compute_temp_deform_eff.h @@ -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(); }; } diff --git a/src/USER-EFF/compute_temp_eff.cpp b/src/USER-EFF/compute_temp_eff.cpp index 26f6539e61..d5e359b942 100644 --- a/src/USER-EFF/compute_temp_eff.cpp +++ b/src/USER-EFF/compute_temp_eff.cpp @@ -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]; diff --git a/src/USER-EFF/compute_temp_eff.h b/src/USER-EFF/compute_temp_eff.h index 683b73a297..197dd4dbc5 100644 --- a/src/USER-EFF/compute_temp_eff.h +++ b/src/USER-EFF/compute_temp_eff.h @@ -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(); }; diff --git a/src/USER-EFF/compute_temp_region_eff.cpp b/src/USER-EFF/compute_temp_region_eff.cpp index be6bee4199..7af81a42c7 100644 --- a/src/USER-EFF/compute_temp_region_eff.cpp +++ b/src/USER-EFF/compute_temp_region_eff.cpp @@ -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; +} + diff --git a/src/USER-EFF/compute_temp_region_eff.h b/src/USER-EFF/compute_temp_region_eff.h index e2d7add698..3340fe8261 100644 --- a/src/USER-EFF/compute_temp_region_eff.h +++ b/src/USER-EFF/compute_temp_region_eff.h @@ -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; }; } diff --git a/src/USER-EFF/fix_langevin_eff.cpp b/src/USER-EFF/fix_langevin_eff.cpp index 7aae42b52d..45b9eba0f7 100644 --- a/src/USER-EFF/fix_langevin_eff.cpp +++ b/src/USER-EFF/fix_langevin_eff.cpp @@ -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; } diff --git a/src/USER-EFF/fix_nh_eff.cpp b/src/USER-EFF/fix_nh_eff.cpp index 0a30b02eee..9dafafb247 100644 --- a/src/USER-EFF/fix_nh_eff.cpp +++ b/src/USER-EFF/fix_nh_eff.cpp @@ -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; } diff --git a/src/USER-EFF/fix_nph_eff.h b/src/USER-EFF/fix_nph_eff.h index 04d265e7bb..9862bf8d99 100644 --- a/src/USER-EFF/fix_nph_eff.h +++ b/src/USER-EFF/fix_nph_eff.h @@ -13,7 +13,7 @@ #ifdef FIX_CLASS -FixStyle(nph/sphere,FixNPHEff) +FixStyle(nph/eff,FixNPHEff) #else diff --git a/src/USER-EFF/fix_npt_eff.cpp b/src/USER-EFF/fix_npt_eff.cpp index 4c2b381fa9..d71a25d461 100644 --- a/src/USER-EFF/fix_npt_eff.cpp +++ b/src/USER-EFF/fix_npt_eff.cpp @@ -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; diff --git a/src/USER-EFF/fix_nve_eff.cpp b/src/USER-EFF/fix_nve_eff.cpp index 0f20bec33a..99866afe03 100644 --- a/src/USER-EFF/fix_nve_eff.cpp +++ b/src/USER-EFF/fix_nve_eff.cpp @@ -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; +} + diff --git a/src/USER-EFF/fix_nve_eff.h b/src/USER-EFF/fix_nve_eff.h index 33d437961e..68be95295f 100644 --- a/src/USER-EFF/fix_nve_eff.h +++ b/src/USER-EFF/fix_nve_eff.h @@ -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; }; } diff --git a/src/USER-EFF/fix_nvt_eff.h b/src/USER-EFF/fix_nvt_eff.h index 3e832094c0..139febcf3d 100644 --- a/src/USER-EFF/fix_nvt_eff.h +++ b/src/USER-EFF/fix_nvt_eff.h @@ -27,7 +27,7 @@ namespace LAMMPS_NS { class FixNVTEff : public FixNHEff { public: FixNVTEff(class LAMMPS *, int, char **); - virtual ~FixNVTEff() {} + ~FixNVTEff() {} }; } diff --git a/src/USER-EFF/fix_nvt_sllod_eff.cpp b/src/USER-EFF/fix_nvt_sllod_eff.cpp index 792208cdce..12bdac820f 100644 --- a/src/USER-EFF/fix_nvt_sllod_eff.cpp +++ b/src/USER-EFF/fix_nvt_sllod_eff.cpp @@ -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)); } diff --git a/src/USER-EFF/fix_nvt_sllod_eff.h b/src/USER-EFF/fix_nvt_sllod_eff.h index e8f9007cfd..852ffb89c4 100644 --- a/src/USER-EFF/fix_nvt_sllod_eff.h +++ b/src/USER-EFF/fix_nvt_sllod_eff.h @@ -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() {} diff --git a/src/USER-EFF/fix_temp_rescale_eff.cpp b/src/USER-EFF/fix_temp_rescale_eff.cpp index 9a31b1023d..e5a0d19be7 100644 --- a/src/USER-EFF/fix_temp_rescale_eff.cpp +++ b/src/USER-EFF/fix_temp_rescale_eff.cpp @@ -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; +} + diff --git a/src/USER-EFF/fix_temp_rescale_eff.h b/src/USER-EFF/fix_temp_rescale_eff.h index 76e7700a64..9ac4c003ad 100644 --- a/src/USER-EFF/fix_temp_rescale_eff.h +++ b/src/USER-EFF/fix_temp_rescale_eff.h @@ -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; }; } diff --git a/src/USER-EFF/pair_eff_cut.cpp b/src/USER-EFF/pair_eff_cut.cpp index 013f8ececa..9b196e98eb 100644 --- a/src/USER-EFF/pair_eff_cut.cpp +++ b/src/USER-EFF/pair_eff_cut.cpp @@ -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); diff --git a/src/USER-EFF/pair_eff_cut.h b/src/USER-EFF/pair_eff_cut.h index e47a13f82b..7519d7335d 100644 --- a/src/USER-EFF/pair_eff_cut.h +++ b/src/USER-EFF/pair_eff_cut.h @@ -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; diff --git a/src/USER-EFF/pair_eff_inline.h b/src/USER-EFF/pair_eff_inline.h index b21e35188d..6f7705b193 100644 --- a/src/USER-EFF/pair_eff_inline.h +++ b/src/USER-EFF/pair_eff_inline.h @@ -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; } /* ---------------------------------------------------------------------- */