diff --git a/src/USER-EFF/README b/src/USER-EFF/README index 6bae419b1..aaa11c349 100644 --- a/src/USER-EFF/README +++ b/src/USER-EFF/README @@ -1,76 +1,91 @@ This package contains a LAMMPS implementation of the electron Force Field (eFF) currently under development at Caltech, as described in A. Jaramillo-Botero, J. Su, Q. An, and W.A. Goddard III, JCC, 2010. The eFF potential was first introduced by Su and Goddard, in 2007. eFF can be viewed as an approximation to QM wave packet dynamics and Fermionic molecular dynamics, combining the ability of electronic structure methods to describe atomic structure, bonding, and chemistry in materials, and of plasma methods to describe nonequilibrium dynamics of large systems with a large number of highly excited electrons. We classify it as a mixed QM-classical approach rather than a conventional force field method, which introduces QM-based terms (a spin-dependent repulsion term to account for the Pauli exclusion principle and the electron wavefunction kinetic energy associated with the Heisenberg principle) that reduce, along with classical electrostatic terms between nuclei and electrons, to the sum of a set of effective pairwise potentials. This makes eFF uniquely suited to simulate materials over a wide range of temperatures and pressures where electronically excited and ionized states of matter can occur and coexist. The necessary customizations to the LAMMPS core are in place to enable the correct handling of explicit electron properties during minimization and dynamics. See the doc page for the pair_style eff/cut command to get started. There are example scripts for using this package in examples/USER/eff. There are auxiliary tools for using this package in tools/eff. The person who created this package is Andres Jaramillo-Botero at CalTech (ajaramil at wag.caltech.edu). Contact him directly if you have questions. ------------------------- AUTHOR INFORMATION: Andres Jaramillo-Botero California Institute of Technology (Caltech) Chemistry and Chemical Engineering, 139-74 1200 E. California Blvd., Pasadena, CA 91125 Phone: (626) 395-3591 e-mail: ajaramil@wag.caltech.edu Co-Authors: Julius Su (jsu@wag.caltech.edu) William A. Goddard III (wag@wag.caltech.edu) 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. + +12/2011: Added support for "zero" option in fix langevin/eff (see doc), and +adjusted fix_langevin_eff.cpp to correctly thermostat between nuclear and electronic dof +(required additional scaling of friction term in the Langevin equations of motion). +Radial electron mass now scales as a function of system dimension. + +Bug fixes: +(10-2011): Thanks to Christian Chenard-Lemire (U Montreal) for reporting a bug in the +fixed pair_eff_cut.cpp fixedcore-pseudocore interactions (an incorrect index and a missing +elec-core call to account for the 2 electrons from the fixed core.) +(12-2011): Corrected undefined natom variable in fix_langevin_eff (recent changes in +main fix_langevin class caused compilation error in user-eff). +Corrected thermostat in fix langevin/eff as described in version 12/2011. + + diff --git a/src/USER-EFF/compute_ke_atom_eff.cpp b/src/USER-EFF/compute_ke_atom_eff.cpp index 1b68b60d2..695cdc887 100644 --- a/src/USER-EFF/compute_ke_atom_eff.cpp +++ b/src/USER-EFF/compute_ke_atom_eff.cpp @@ -1,113 +1,115 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Andres Jaramillo-Botero ------------------------------------------------------------------------- */ #include "math.h" #include "string.h" #include "compute_ke_atom_eff.h" #include "atom.h" #include "update.h" #include "modify.h" #include "comm.h" #include "force.h" #include "memory.h" #include "error.h" +#include "domain.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeKEAtomEff::ComputeKEAtomEff(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 3) error->all(FLERR,"Illegal compute ke/atom/eff command"); peratom_flag = 1; size_peratom_cols = 0; nmax = 0; ke = NULL; // error check if (!atom->electron_flag) error->all(FLERR,"Compute ke/atom/eff requires atom style electron"); } /* ---------------------------------------------------------------------- */ ComputeKEAtomEff::~ComputeKEAtomEff() { memory->destroy(ke); } /* ---------------------------------------------------------------------- */ void ComputeKEAtomEff::init() { int count = 0; for (int i = 0; i < modify->ncompute; i++) if (strcmp(modify->compute[i]->style,"ke/atom/eff") == 0) count++; if (count > 1 && comm->me == 0) error->warning(FLERR,"More than one compute ke/atom/eff"); } /* ---------------------------------------------------------------------- */ void ComputeKEAtomEff::compute_peratom() { invoked_peratom = update->ntimestep; // grow ke array if necessary if (atom->nlocal > nmax) { memory->destroy(ke); nmax = atom->nmax; memory->create(ke,nmax,"compute/ke/atom/eff:ke"); vector_atom = ke; } // compute kinetic energy for each atom in group double mvv2e = force->mvv2e; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *mask = atom->mask; int *type = atom->type; int nlocal = atom->nlocal; + double mefactor = domain->dimension/4.0; 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 (fabs(spin[i])==1) - ke[i] += 0.5 * mvv2e * mass[type[i]] * ervel[i]*ervel[i] * 0.75; + ke[i] += 0.5 * mvv2e * mefactor * mass[type[i]] * ervel[i]*ervel[i]; } else ke[i] = 0.0; } } /* ---------------------------------------------------------------------- memory usage of local atom-based array ------------------------------------------------------------------------- */ double ComputeKEAtomEff::memory_usage() { double bytes = nmax * sizeof(double); return bytes; } diff --git a/src/USER-EFF/compute_ke_eff.cpp b/src/USER-EFF/compute_ke_eff.cpp index 43a2a952f..5da7d2574 100644 --- a/src/USER-EFF/compute_ke_eff.cpp +++ b/src/USER-EFF/compute_ke_eff.cpp @@ -1,81 +1,82 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Andres Jaramillo-Botero ------------------------------------------------------------------------- */ #include "mpi.h" #include "math.h" #include "compute_ke_eff.h" #include "atom.h" #include "update.h" #include "force.h" #include "domain.h" #include "group.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeKEEff::ComputeKEEff(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 3) error->all(FLERR,"Illegal compute ke/eff command"); scalar_flag = 1; extscalar = 1; // error check if (!atom->electron_flag) error->all(FLERR,"Compute ke/eff requires atom style electron"); } /* ---------------------------------------------------------------------- */ void ComputeKEEff::init() { pfactor = 0.5 * force->mvv2e; } /* ---------------------------------------------------------------------- */ double ComputeKEEff::compute_scalar() { invoked_scalar = update->ntimestep; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *mask = atom->mask; int *type = atom->type; int nlocal = atom->nlocal; + double mefactor = domain->dimension/4.0; double ke = 0.0; 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 (fabs(spin[i])==1) ke += 0.75*mass[type[i]]*ervel[i]*ervel[i]; + if (fabs(spin[i])==1) ke += mefactor*mass[type[i]]*ervel[i]*ervel[i]; } } MPI_Allreduce(&ke,&scalar,1,MPI_DOUBLE,MPI_SUM,world); scalar *= pfactor; return scalar; } diff --git a/src/USER-EFF/compute_temp_deform_eff.cpp b/src/USER-EFF/compute_temp_deform_eff.cpp index ccba4d782..68e73a83d 100644 --- a/src/USER-EFF/compute_temp_deform_eff.cpp +++ b/src/USER-EFF/compute_temp_deform_eff.cpp @@ -1,316 +1,319 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Andres Jaramillo-Botero (Caltech) ------------------------------------------------------------------------- */ #include "mpi.h" #include "math.h" #include "string.h" #include "compute_temp_deform_eff.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 /* ---------------------------------------------------------------------- */ ComputeTempDeformEff::ComputeTempDeformEff(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 3) error->all(FLERR,"Illegal compute temp/deform/eff command"); if (!atom->electron_flag) error->all(FLERR,"Compute temp/deform/eff requires atom style electron"); 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(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(FLERR,"Using compute temp/deform/eff with inconsistent " "fix deform remap option"); break; } if (i == modify->nfix && comm->me == 0) error->warning(FLERR,"Using compute temp/deform/eff with no fix deform defined"); } /* ---------------------------------------------------------------------- */ void ComputeTempDeformEff::dof_compute() { double natoms = group->count(igroup); dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; // just include nuclear dof int *spin = atom->spin; int *mask = atom->mask; int nlocal = atom->nlocal; int one = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (fabs(spin[i]) == 1) one++; } int nelectrons; MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world); + // Assume 3/2 k T per nucleus dof -= domain->dimension * nelectrons; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempDeformEff::compute_scalar() { double lamda[3],vstream[3],vthermal[3]; invoked_scalar = update->ntimestep; double **x = atom->x; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; + double mefactor = domain->dimension/4.0; // lamda = 0-1 triclinic lamda coords // vstream = streaming velocity = Hrate*lamda + Hratelo // vthermal = thermal velocity = v - vstream double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; double t = 0.0; for (int i = 0; i < nlocal; i++) 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]; 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 (mass) { t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2])* mass[type[i]]; - if (fabs(spin[i])==1) t += 0.75*mass[type[i]]*ervel[i]*ervel[i]; + if (fabs(spin[i])==1) t += mefactor*mass[type[i]]*ervel[i]*ervel[i]; } } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempDeformEff::compute_vector() { double lamda[3],vstream[3],vthermal[3]; invoked_vector = update->ntimestep; double **x = atom->x; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; - + double mefactor = domain->dimension/4.0; + double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; double massone,t[6]; for (int i = 0; i < 6; i++) t[i] = 0.0; for (int i = 0; i < nlocal; i++) 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]; 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]; 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 (fabs(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]; + t[0] += mefactor * massone * ervel[i]*ervel[i]; + t[1] += mefactor * massone * ervel[i]*ervel[i]; + t[2] += mefactor * massone * ervel[i]*ervel[i]; } } 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(vbiasall); maxbias = atom->nmax; memory->create(vbiasall,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_eff.cpp b/src/USER-EFF/compute_temp_eff.cpp index e2089a0d6..cf406c53e 100644 --- a/src/USER-EFF/compute_temp_eff.cpp +++ b/src/USER-EFF/compute_temp_eff.cpp @@ -1,164 +1,166 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Andres Jaramillo-Botero (Caltech) ------------------------------------------------------------------------- */ #include "mpi.h" #include "math.h" #include "string.h" #include "compute_temp_eff.h" #include "atom.h" #include "update.h" #include "force.h" #include "domain.h" #include "modify.h" #include "fix.h" #include "group.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTempEff::ComputeTempEff(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (!atom->electron_flag) error->all(FLERR,"Compute temp/eff requires atom style electron"); 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(); } /* ---------------------------------------------------------------------- */ void ComputeTempEff::dof_compute() { double natoms = group->count(igroup); dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; int *spin = atom->spin; int *mask = atom->mask; int nlocal = atom->nlocal; int one = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (fabs(spin[i])==1) one++; } int nelectrons; MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world); - // average over nuclear dof only + // Assume 3/2 k T per nucleus dof -= domain->dimension * nelectrons; if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempEff::compute_scalar() { invoked_scalar = update->ntimestep; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; + double mefactor = domain->dimension/4.0; double t = 0.0; 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 (fabs(spin[i])==1) t += 0.75*mass[type[i]]*ervel[i]*ervel[i]; + if (fabs(spin[i])==1) t += mefactor*mass[type[i]]*ervel[i]*ervel[i]; } } } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempEff::compute_vector() { int i; invoked_vector = update->ntimestep; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; + double mefactor = domain->dimension/4.0; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { 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 (fabs(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]; + t[0] += mefactor*massone*ervel[i]*ervel[i]; + t[1] += mefactor*massone*ervel[i]*ervel[i]; + t[2] += mefactor*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; } diff --git a/src/USER-EFF/compute_temp_region_eff.cpp b/src/USER-EFF/compute_temp_region_eff.cpp index f43e92095..0dd2bf3cf 100644 --- a/src/USER-EFF/compute_temp_region_eff.cpp +++ b/src/USER-EFF/compute_temp_region_eff.cpp @@ -1,276 +1,277 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Andres Jaramillo-Botero (Caltech) ------------------------------------------------------------------------- */ #include "mpi.h" #include "math.h" #include "string.h" #include "compute_temp_region_eff.h" #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) : Compute(lmp, narg, arg) { if (!atom->electron_flag) error->all(FLERR,"Compute temp/region/eff requires atom style electron"); if (narg != 4) error->all(FLERR,"Illegal compute temp/region/eff command"); iregion = domain->find_region(arg[3]); if (iregion == -1) error->all(FLERR,"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(vbiasall); delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempRegionEff::init() { // set index and check validity of region iregion = domain->find_region(idregion); if (iregion == -1) error->all(FLERR,"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; } /* ---------------------------------------------------------------------- */ double ComputeTempRegionEff::compute_scalar() { invoked_scalar = update->ntimestep; double **x = atom->x; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; + double mefactor = domain->dimension/4.0; Region *region = domain->regions[iregion]; int count = 0; + int ecount = 0; double t = 0.0; 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]) * mass[type[i]]; - if (fabs(spin[i])==1) t += 0.75*mass[type[i]]*ervel[i]*ervel[i]; + if (fabs(spin[i])==1) { + t += mefactor*mass[type[i]]*ervel[i]*ervel[i]; + ecount++; + } } } double tarray[2],tarray_all[2]; - tarray[0] = count; + // Assume 3/2 k T per nucleus + tarray[0] = count-ecount; tarray[1] = t; MPI_Allreduce(tarray,tarray_all,2,MPI_DOUBLE,MPI_SUM,world); dof = domain->dimension * tarray_all[0] - extra_dof; 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 (fabs(spin[i])==1) one++; } - int nelectrons_region; - MPI_Allreduce(&one,&nelectrons_region,1,MPI_INT,MPI_SUM,world); - - // average over nuclear dof only - - dof -= domain->dimension * nelectrons_region ; if (dof > 0) scalar = force->mvv2e * tarray_all[1] / (dof * force->boltz); else scalar = 0.0; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempRegionEff::compute_vector() { int i; invoked_vector = update->ntimestep; double **x = atom->x; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; + double mefactor = domain->dimension/4.0; Region *region = domain->regions[iregion]; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { 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 (fabs(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]; + t[0] += mefactor * massone * ervel[i]*ervel[i]; + t[1] += mefactor * massone * ervel[i]*ervel[i]; + t[2] += mefactor * 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(vbiasall); maxbias = atom->nmax; memory->create(vbiasall,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/fix_langevin_eff.cpp b/src/USER-EFF/fix_langevin_eff.cpp index 7328fdb23..faff9909d 100644 --- a/src/USER-EFF/fix_langevin_eff.cpp +++ b/src/USER-EFF/fix_langevin_eff.cpp @@ -1,251 +1,435 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Andres Jaramillo-Botero ------------------------------------------------------------------------- */ #include "mpi.h" #include "math.h" #include "string.h" #include "stdlib.h" #include "fix_langevin_eff.h" +#include "math_extra.h" #include "atom.h" #include "force.h" #include "update.h" #include "modify.h" #include "compute.h" #include "domain.h" #include "region.h" #include "respa.h" #include "comm.h" +#include "input.h" +#include "variable.h" #include "random_mars.h" #include "memory.h" #include "error.h" +#include "group.h" using namespace LAMMPS_NS; enum{NOBIAS,BIAS}; +enum{CONSTANT,EQUAL,ATOM}; + +#define SINERTIA 0.4 // moment of inertia prefactor for sphere +#define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid /* ---------------------------------------------------------------------- */ FixLangevinEff::FixLangevinEff(LAMMPS *lmp, int narg, char **arg) : FixLangevin(lmp, narg, arg) { erforcelangevin = NULL; } /* ---------------------------------------------------------------------- */ FixLangevinEff::~FixLangevinEff() { memory->destroy(erforcelangevin); } /* ---------------------------------------------------------------------- */ void FixLangevinEff::post_force_no_tally() { - double gamma1,gamma2; + double gamma1,gamma2,t_target; double **v = atom->v; double **f = atom->f; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - double *ervel = atom->ervel; double *erforce = atom->erforce; int *spin = atom->spin; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + double mefactor = domain->dimension/4.0; + double sqrtmefactor = sqrt(mefactor); double delta = update->ntimestep - update->beginstep; delta /= update->endstep - update->beginstep; - double t_target = t_start + delta * (t_stop-t_start); - double tsqrt = sqrt(t_target); + + // set current t_target and t_sqrt + // if variable temp, evaluate variable, wrap with clear/add + // reallocate tforce array if necessary + + if (tstyle == CONSTANT) { + t_target = t_start + delta * (t_stop-t_start); + tsqrt = sqrt(t_target); + } else { + modify->clearstep_compute(); + if (tstyle == EQUAL) { + t_target = input->variable->compute_equal(tvar); + if (t_target < 0.0) + error->one(FLERR,"Fix langevin/eff variable returned negative temperature"); + tsqrt = sqrt(t_target); + } else { + if (nlocal > maxatom2) { + maxatom2 = atom->nmax; + memory->destroy(tforce); + memory->create(tforce,maxatom2,"langevin/eff:tforce"); + } + input->variable->compute_atom(tvar,igroup,tforce,1,0); + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (tforce[i] < 0.0) + error->one(FLERR, + "Fix langevin/eff variable returned negative temperature"); + } + modify->addstep_compute(update->ntimestep + 1); + } // apply damping and thermostat to atoms in group // for BIAS: // calculate temperature since some computes require temp // computed on current nlocal atoms to remove bias // test v = 0 since some computes mask non-participating atoms via v = 0 // and added force has extra term not multiplied by v = 0 + // for ZEROFLAG: + // sum random force over all atoms in group + // subtract sum/particles from each atom in group + + double fran[4],fsum[4],fsumall[4]; + fsum[0] = fsum[1] = fsum[2] = fsum[3] = 0.0; + + double boltz = force->boltz; + double dt = update->dt; + double mvv2e = force->mvv2e; + double ftm2v = force->ftm2v; + + int particles = group->count(igroup); + if (zeroflag) { + if (particles == 0) + error->all(FLERR,"Cannot zero Langevin force of 0 atoms/electrons"); + } + + // find number of electrons in group + int dof,fix_dof; + dof = domain->dimension * particles; + fix_dof = 0; + for (int i = 0; i < modify->nfix; i++) + fix_dof += modify->fix[i]->dof(igroup); + + // extra_dof = domain->dimension + dof -= domain->dimension + fix_dof; + + int one = 0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (fabs(spin[i])==1) one++; + } + int nelectrons, dofelectrons, dofnuclei; + MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world); + dofelectrons = domain->dimension*nelectrons; + dofnuclei = dof-dofelectrons; + + // thermal partitioning factor between nuclei and electrons + // extra dof from electron size + double gfactor3=(double) (dof+nelectrons)/dofnuclei; if (which == NOBIAS) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - gamma1 = gfactor1[type[i]]; + if (tstyle == ATOM) tsqrt = sqrt(tforce[i]); + gamma1 = gfactor1[type[i]] * gfactor3; 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); + fran[0] = gamma2*(random->uniform()-0.5); + fran[1] = gamma2*(random->uniform()-0.5); + fran[2] = gamma2*(random->uniform()-0.5); + f[i][0] += gamma1*v[i][0] + fran[0]; + f[i][1] += gamma1*v[i][1] + fran[1]; + f[i][2] += gamma1*v[i][2] + fran[2]; + fsum[0] += fran[0]; + fsum[1] += fran[1]; + fsum[2] += fran[2]; + if (fabs(spin[i])==1) { + fran[3] = sqrtmefactor*gamma2*(random->uniform()-0.5); + erforce[i] += mefactor*gamma1*ervel[i]+fran[3]; + fsum[3] += fran[3]; + } } } } else if (which == BIAS) { - double tmp = temperature->compute_scalar(); + temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - gamma1 = gfactor1[type[i]]; + if (tstyle == ATOM) tsqrt = sqrt(tforce[i]); + gamma1 = gfactor1[type[i]] * gfactor3; gamma2 = gfactor2[type[i]] * tsqrt; temperature->remove_bias(i,v[i]); + fran[0] = gamma2*(random->uniform()-0.5); + fran[1] = gamma2*(random->uniform()-0.5); + fran[2] = gamma2*(random->uniform()-0.5); if (v[i][0] != 0.0) - f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5); + f[i][0] += gamma1*v[i][0] + fran[0]; if (v[i][1] != 0.0) - f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5); + f[i][1] += gamma1*v[i][1] + fran[1]; if (v[i][2] != 0.0) - f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5); - if (abs(spin[i])==1 && ervel[i] != 0.0) - erforce[i] += 0.75*gamma1*ervel[i] + - 0.866025404*gamma2*(random->uniform()-0.5); + f[i][2] += gamma1*v[i][2] + fran[2]; + fsum[0] += fran[0]; + fsum[1] += fran[1]; + fsum[2] += fran[2]; + if (fabs(spin[i])==1) { + fran[3] = sqrtmefactor*gamma2*(random->uniform()-0.5); + if (ervel[i] != 0.0) erforce[i] += mefactor*gamma1*ervel[i]+fran[3]; + fsum[3] += fran[3]; + } temperature->restore_bias(i,v[i]); } } } + + // set total force to zero + + if (zeroflag) { + MPI_Allreduce(fsum,fsumall,3,MPI_DOUBLE,MPI_SUM,world); + fsumall[0] /= particles; + fsumall[1] /= particles; + fsumall[2] /= particles; + fsumall[3] /= nelectrons; + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + f[i][0] -= fsumall[0]; + f[i][1] -= fsumall[1]; + f[i][2] -= fsumall[2]; + if (fabs(spin[i])==1) erforce[i] -= fsumall[3]; + } + } + } } /* ---------------------------------------------------------------------- */ void FixLangevinEff::post_force_tally() { - double gamma1,gamma2; + double gamma1,gamma2,t_target; - // reallocate flangevin if necessary + // reallocate flangevin and erforcelangevin if necessary - if (atom->nmax > nmax) { + if (atom->nlocal > maxatom1) { memory->destroy(flangevin); memory->destroy(erforcelangevin); - nmax = atom->nmax; - memory->create(flangevin,nmax,3,"langevin:flangevin"); - memory->create(erforcelangevin,nmax,"langevin/eff:erforcelangevin"); + maxatom1 = atom->nmax; + memory->create(flangevin,maxatom1,3,"langevin/eff:flangevin"); + memory->create(erforcelangevin,maxatom1,"langevin/eff:erforcelangevin"); } double **v = atom->v; double **f = atom->f; - 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 mefactor = domain->dimension/4.0; + double sqrtmefactor = sqrt(mefactor); + + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; double delta = update->ntimestep - update->beginstep; delta /= update->endstep - update->beginstep; - double t_target = t_start + delta * (t_stop-t_start); - double tsqrt = sqrt(t_target); + + // set current t_target and t_sqrt + // if variable temp, evaluate variable, wrap with clear/add + // reallocate tforce array if necessary + + if (tstyle == CONSTANT) { + t_target = t_start + delta * (t_stop-t_start); + tsqrt = sqrt(t_target); + } else { + modify->clearstep_compute(); + if (tstyle == EQUAL) { + t_target = input->variable->compute_equal(tvar); + if (t_target < 0.0) + error->one(FLERR,"Fix langevin/eff variable returned negative temperature"); + tsqrt = sqrt(t_target); + } else { + if (nlocal > maxatom2) { + maxatom2 = atom->nmax; + memory->destroy(tforce); + memory->create(tforce,maxatom2,"langevin/eff:tforce"); + } + input->variable->compute_atom(tvar,igroup,tforce,1,0); + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (tforce[i] < 0.0) + error->one(FLERR, + "Fix langevin/eff variable returned negative temperature"); + } + modify->addstep_compute(update->ntimestep + 1); + } // apply damping and thermostat to appropriate atoms // for BIAS: // calculate temperature since some computes require temp // computed on current nlocal atoms to remove bias // test v = 0 since some computes mask non-participating atoms via v = 0 // and added force has extra term not multiplied by v = 0 + double boltz = force->boltz; + double dt = update->dt; + double mvv2e = force->mvv2e; + double ftm2v = force->ftm2v; + + int particles = group->count(igroup); + if (zeroflag) { + if (particles == 0) + error->all(FLERR,"Cannot zero Langevin force of 0 atoms/electrons"); + } + + // find number of electrons in group + int dof,fix_dof; + dof = domain->dimension * particles; + fix_dof = 0; + for (int i = 0; i < modify->nfix; i++) + fix_dof += modify->fix[i]->dof(igroup); + + // extra_dof = domain->dimension + dof -= domain->dimension + fix_dof; + + int one = 0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (fabs(spin[i])==1) one++; + } + int nelectrons, dofelectrons, dofnuclei; + MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world); + dofelectrons = domain->dimension*nelectrons; + dofnuclei = dof-dofelectrons; + + // thermal partitioning factor between nuclei and electrons + // extra dof from electron size + double gfactor3=(double) (dof+nelectrons)/dofnuclei; + if (which == NOBIAS) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - gamma1 = gfactor1[type[i]]; + if (tstyle == ATOM) tsqrt = sqrt(tforce[i]); + gamma1 = gfactor1[type[i]] * gfactor3; 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]; + if (fabs(spin[i])==1) { + erforcelangevin[i] = mefactor*gamma1*ervel[i]+sqrtmefactor*gamma2*(random->uniform()-0.5); + erforce[i] += erforcelangevin[i]; + } } } } else if (which == BIAS) { - double tmp = temperature->compute_scalar(); + temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - gamma1 = gfactor1[type[i]]; + if (tstyle == ATOM) tsqrt = sqrt(tforce[i]); + gamma1 = gfactor1[type[i]] * gfactor3; 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]; + if (fabs(spin[i])==1) { + erforcelangevin[i] = mefactor*gamma1*ervel[i]+sqrtmefactor*gamma2*(random->uniform()-0.5); + if (ervel[i] != 0.0) erforce[i] += erforcelangevin[i]; + else erforcelangevin[i] = 0.0; + } temperature->restore_bias(i,v[i]); } } } } /* ---------------------------------------------------------------------- tally energy transfer to thermal reservoir ------------------------------------------------------------------------- */ void FixLangevinEff::end_of_step() { if (!tally) return; 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) { energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] + - flangevin[i][2]*v[i][2]; - if (abs(spin[i])==1) energy_onestep += erforcelangevin[i]; + flangevin[i][2]*v[i][2]; + if (fabs(spin[i])==1) energy_onestep += erforcelangevin[i]; } energy += energy_onestep*update->dt; } /* ---------------------------------------------------------------------- */ double FixLangevinEff::compute_scalar() { - if (!tally) return 0.0; + if (!tally || flangevin == NULL || erforcelangevin == NULL) return 0.0; // capture the very first energy transfer to thermal reservoir 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) { energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] + flangevin[i][2]*v[i][2]; - if (abs(spin[i])==1) energy_onestep += erforcelangevin[i]; + if (fabs(spin[i])==1) energy_onestep += erforcelangevin[i]; } energy = 0.5*energy_onestep*update->dt; } double energy_me = energy - 0.5*energy_onestep*update->dt; double energy_all; MPI_Allreduce(&energy_me,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); return -energy_all; } + diff --git a/src/USER-EFF/fix_nh_eff.cpp b/src/USER-EFF/fix_nh_eff.cpp index e3dde8417..12bcb9076 100644 --- a/src/USER-EFF/fix_nh_eff.cpp +++ b/src/USER-EFF/fix_nh_eff.cpp @@ -1,110 +1,112 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Andres Jaramillo-Botero (Caltech) ------------------------------------------------------------------------- */ #include "math.h" #include "fix_nh_eff.h" #include "atom.h" #include "atom_vec.h" #include "group.h" #include "error.h" +#include "domain.h" using namespace LAMMPS_NS; enum{NOBIAS,BIAS}; /* ---------------------------------------------------------------------- */ FixNHEff::FixNHEff(LAMMPS *lmp, int narg, char **arg) : FixNH(lmp, narg, arg) { if (!atom->electron_flag) error->all(FLERR,"Fix nvt/nph/npt/eff requires atom style electron"); } /* ---------------------------------------------------------------------- perform half-step update of electron radial velocities -----------------------------------------------------------------------*/ void FixNHEff::nve_v() { // standard nve_v velocity update FixNH::nve_v(); double *erforce = atom->erforce; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; + double mefactor = domain->dimension/4.0; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; int itype; double dtfm; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { if (fabs(spin[i])==1) { dtfm = dtf / mass[type[i]]; - ervel[i] = dtfm * erforce[i] / 0.75; + ervel[i] = dtfm * erforce[i] / mefactor; } } } } /* ---------------------------------------------------------------------- perform full-step update of electron radii -----------------------------------------------------------------------*/ void FixNHEff::nve_x() { // standard nve_x position update FixNH::nve_x(); double *eradius = atom->eradius; double *ervel = atom->ervel; int *spin = atom->spin; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) if (fabs(spin[i])==1) eradius[i] += dtv * ervel[i]; } /* ---------------------------------------------------------------------- perform half-step scaling of electron radial velocities -----------------------------------------------------------------------*/ void FixNHEff::nh_v_temp() { // standard nh_v_temp velocity scaling FixNH::nh_v_temp(); double *ervel = atom->ervel; int *spin = atom->spin; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) if (fabs(spin[i])==1) ervel[i] *= factor_eta; } diff --git a/src/USER-EFF/fix_nve_eff.cpp b/src/USER-EFF/fix_nve_eff.cpp index 1466c777d..9b4613df3 100644 --- a/src/USER-EFF/fix_nve_eff.cpp +++ b/src/USER-EFF/fix_nve_eff.cpp @@ -1,171 +1,174 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Andres Jaramillo-Botero (Caltech) ------------------------------------------------------------------------- */ #include "math.h" #include "stdio.h" #include "string.h" #include "fix_nve_eff.h" #include "atom.h" #include "force.h" #include "update.h" #include "respa.h" #include "error.h" #include "math.h" +#include "domain.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ FixNVEEff::FixNVEEff(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { if (!atom->electron_flag) error->all(FLERR,"Fix nve/eff requires atom style electron"); 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 (strstr(update->integrate_style,"respa")) step_respa = ((Respa *) update->integrate)->step; } /* ---------------------------------------------------------------------- allow for both per-type and per-atom mass ------------------------------------------------------------------------- */ void FixNVEEff::initial_integrate(int vflag) { double dtfm; // update v,vr and x,radius of atoms in group double **x = atom->x; double *eradius = atom->eradius; double **v = atom->v; double *ervel = atom->ervel; double **f = atom->f; double *erforce = atom->erforce; double *mass = atom->mass; int *spin = atom->spin; + double mefactor = domain->dimension/4.0; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; // x + dt * [v + 0.5 * dt * (f / m)]; if (mass) { 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 (fabs(spin[i])==1) { - ervel[i] += dtfm * erforce[i] / 0.75; + ervel[i] += dtfm * erforce[i] / mefactor; eradius[i] += dtv * ervel[i]; } } } } } /* ---------------------------------------------------------------------- */ void FixNVEEff::final_integrate() { double dtfm; double **v = atom->v; double *ervel = atom->ervel; double *erforce = atom->erforce; double **f = atom->f; double *mass = atom->mass; int *spin = atom->spin; + double mefactor = domain->dimension/4.0; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; // dyn_v[i] += m * dt * dyn_f[i]; if (mass) { 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 (fabs(spin[i])==1) - ervel[i] += dtfm * erforce[i] / 0.75; + ervel[i] += dtfm * erforce[i] / mefactor; } } } } /* ---------------------------------------------------------------------- */ 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/pair_eff_cut.cpp b/src/USER-EFF/pair_eff_cut.cpp index 8f9913423..ee191a0d1 100644 --- a/src/USER-EFF/pair_eff_cut.cpp +++ b/src/USER-EFF/pair_eff_cut.cpp @@ -1,954 +1,964 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Andres Jaramillo-Botero ------------------------------------------------------------------------- */ #include "math.h" #include "stdio.h" #include "stdlib.h" #include "string.h" #include "pair_eff_cut.h" #include "pair_eff_inline.h" #include "atom.h" #include "update.h" #include "min.h" #include "domain.h" #include "comm.h" #include "force.h" #include "neighbor.h" #include "neigh_list.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ PairEffCut::PairEffCut(LAMMPS *lmp) : Pair(lmp) { single_enable = 0; nmax = 0; min_eradius = NULL; min_erforce = NULL; nextra = 4; pvector = new double[nextra]; } /* ---------------------------------------------------------------------- */ PairEffCut::~PairEffCut() { delete [] pvector; memory->destroy(min_eradius); memory->destroy(min_erforce); if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); memory->destroy(cut); } } /* ---------------------------------------------------------------------- */ void PairEffCut::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; 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; 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; double **x = atom->x; double **f = atom->f; double *q = atom->q; double *erforce = atom->erforce; double *eradius = atom->eradius; int *spin = atom->spin; int *type = atom->type; int nlocal = atom->nlocal; int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; jlist = firstneigh[i]; jnum = numneigh[i]; // add electron wavefuntion kinetic energy (not pairwise) 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; // Tally energy and compute radial atomic virial contribution if (evflag) { ev_tally_eff(i,i,nlocal,newton_pair,energy,0.0); if (flexible_pressure_flag) // iff flexible pressure flag on 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]; j &= NEIGHMASK; 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); jtype = type[j]; if (rsq < cutsq[itype][jtype]) { 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) { double qxq = q[i]*q[j]; 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; // 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); } } // electron (i) - nucleus (j) Coul interaction 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; // 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 (i) - electron (j) interactions else if (abs(spin[i]) == 1 && abs(spin[j]) == 1) { e1rforce = e2rforce = 0.0; s_e1rforce = s_e2rforce = 0.0; 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); // Apply conversion factor epauli *= hhmss2e; s_fpair *= hhmss2e; e1rforce = spline * (qqrd2e * e1rforce + hhmss2e * s_e1rforce); erforce[i] += e1rforce; e2rforce = spline * (qqrd2e * e2rforce + hhmss2e * s_e2rforce); erforce[j] += e2rforce; // 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); + ElecCoreElec(q[i],rc,eradius[i],eradius[j],&ecoul, + &fpair,&e2rforce); + PauliCoreElec(rc,eradius[j],&ecp_epauli,&ecp_fpair, + &ecp_e2rforce,PAULI_CORE_A, PAULI_CORE_B, + PAULI_CORE_C); PauliCoreElec(rc,eradius[j],&ecp_epauli,&ecp_fpair, &ecp_e2rforce,PAULI_CORE_A, PAULI_CORE_B, PAULI_CORE_C); } // 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) { + if (abs(spin[i]) == 1) { ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul, &fpair,&e1rforce); PauliCoreElec(rc,eradius[i],&ecp_epauli,&ecp_fpair, &ecp_e1rforce,PAULI_CORE_A,PAULI_CORE_B, PAULI_CORE_C); } else { double qxq = q[i]*q[j]; ElecCoreNuc(qxq,rc,eradius[i],&ecoul,&fpair); ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul, &fpair,&e1rforce); + ElecCoreElec(q[j],rc,eradius[j],eradius[i],&ecoul, + &fpair,&e1rforce); + PauliCoreElec(rc,eradius[i],&ecp_epauli,&ecp_fpair, + &ecp_e1rforce,PAULI_CORE_A, PAULI_CORE_B, + PAULI_CORE_C); PauliCoreElec(rc,eradius[i],&ecp_epauli,&ecp_fpair, &ecp_e1rforce,PAULI_CORE_A, PAULI_CORE_B, PAULI_CORE_C); } // 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 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 (domain->xperiodic == 1 || domain->yperiodic == 1 || domain->zperiodic == 1) { delx = domain->boxhi[0]-domain->boxlo[0]; dely = domain->boxhi[1]-domain->boxlo[1]; delz = domain->boxhi[2]-domain->boxlo[2]; half_box_length = 0.5 * MIN(delx, MIN(dely, delz)); if (eradius[i] > half_box_length) { dr = eradius[i]-half_box_length; errestrain=0.5*kfactor*dr*dr; e1rforce=-kfactor*dr; if (eflag_global) pvector[3] += errestrain; erforce[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_fdotr_compute(); if (flexible_pressure_flag) virial_eff_compute(); } } /* ---------------------------------------------------------------------- eff-specific contribution to global virial ------------------------------------------------------------------------- */ void PairEffCut::virial_eff_compute() { double *eradius = atom->eradius; double *erforce = atom->erforce; double e_virial; int *spin = atom->spin; // sum over force on all particles including ghosts if (neighbor->includegroup == 0) { int nall = atom->nlocal + atom->nghost; for (int i = 0; i < nall; i++) { if (spin[i]) { e_virial = erforce[i]*eradius[i]/3; virial[0] += e_virial; virial[1] += e_virial; virial[2] += e_virial; } } // neighbor includegroup flag is set // sum over force on initial nfirst particles and ghosts } else { int nall = atom->nfirst; for (int i = 0; i < nall; i++) { if (spin[i]) { e_virial = erforce[i]*eradius[i]/3; virial[0] += e_virial; virial[1] += e_virial; virial[2] += e_virial; } } nall = atom->nlocal + atom->nghost; for (int i = atom->nlocal; i < nall; i++) { if (spin[i]) { e_virial = erforce[i]*eradius[i]/3; virial[0] += e_virial; virial[1] += e_virial; virial[2] += e_virial; } } } } /* ---------------------------------------------------------------------- tally eng_vdwl and virial into per-atom accumulators for virial radial electronic contributions ------------------------------------------------------------------------- */ void PairEffCut::ev_tally_eff(int i, int j, int nlocal, int newton_pair, double energy, double e_virial) { double energyhalf; double partial_evirial = e_virial/3.0; double half_partial_evirial = partial_evirial/2; int *spin = atom->spin; if (eflag_either) { if (eflag_global) { 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) { 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] += half_partial_evirial; virial[1] += half_partial_evirial; virial[2] += half_partial_evirial; } if (spin[j] && j < nlocal) { 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] += 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] += half_partial_evirial; vatom[j][1] += half_partial_evirial; vatom[j][2] += half_partial_evirial; } } } } } /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ void PairEffCut::allocate() { allocated = 1; int n = atom->ntypes; memory->create(setflag,n+1,n+1,"pair:setflag"); for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; memory->create(cutsq,n+1,n+1,"pair:cutsq"); memory->create(cut,n+1,n+1,"pair:cut"); } /* --------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ void PairEffCut::settings(int narg, char **arg) { if (narg != 1 && narg != 3 && narg != 4 && narg != 7) error->all(FLERR,"Illegal pair_style command"); // Defaults ECP parameters for Si PAULI_CORE_A = 0.320852; PAULI_CORE_B = 2.283269; PAULI_CORE_C = 0.814857; if (narg == 1) { cut_global = force->numeric(arg[0]); limit_size_flag = 0; flexible_pressure_flag = 0; } else if (narg == 3) { cut_global = force->numeric(arg[0]); limit_size_flag = force->inumeric(arg[1]); flexible_pressure_flag = force->inumeric(arg[2]); } else if (narg == 4) { cut_global = force->numeric(arg[0]); limit_size_flag = 0; flexible_pressure_flag = 0; if (strcmp(arg[1],"ecp") != 0) error->all(FLERR,"Illegal pair_style command"); else { PAULI_CORE_A = force->numeric(arg[2]); PAULI_CORE_B = force->numeric(arg[3]); PAULI_CORE_C = force->numeric(arg[4]); } } else if (narg == 7) { cut_global = force->numeric(arg[0]); limit_size_flag = force->inumeric(arg[1]); flexible_pressure_flag = force->inumeric(arg[2]); if (strcmp(arg[3],"ecp") != 0) error->all(FLERR,"Illegal pair_style command"); else { PAULI_CORE_A = force->numeric(arg[4]); PAULI_CORE_B = force->numeric(arg[5]); PAULI_CORE_C = force->numeric(arg[6]); } } // 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(FLERR,"Check your units"); // reset cutoffs that have been explicitly set if (allocated) { int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i+1; j <= atom->ntypes; j++) if (setflag[i][j]) cut[i][j] = cut_global; } } /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ void PairEffCut::coeff(int narg, char **arg) { if (narg < 2 || narg > 3) error->all(FLERR,"Incorrect args for pair coefficients"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; force->bounds(arg[0],atom->ntypes,ilo,ihi); force->bounds(arg[1],atom->ntypes,jlo,jhi); double cut_one = cut_global; if (narg == 3) cut_one = atof(arg[2]); int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { cut[i][j] = cut_one; setflag[i][j] = 1; count++; } } if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ void PairEffCut::init_style() { // error and warning checks if (!atom->q_flag || !atom->spin_flag || !atom->eradius_flag || !atom->erforce_flag) error->all(FLERR,"Pair eff/cut requires atom attributes " "q, spin, eradius, erforce"); // add hook to minimizer for eradius and erforce 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(FLERR,"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); } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairEffCut::init_one(int i, int j) { if (setflag[i][j] == 0) cut[i][j] = mix_distance(cut[i][i],cut[j][j]); return cut[i][j]; } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairEffCut::write_restart(FILE *fp) { write_restart_settings(fp); int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { fwrite(&setflag[i][j],sizeof(int),1,fp); if (setflag[i][j]) fwrite(&cut[i][j],sizeof(double),1,fp); } } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairEffCut::read_restart(FILE *fp) { read_restart_settings(fp); allocate(); int i,j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); if (setflag[i][j]) { if (me == 0) fread(&cut[i][j],sizeof(double),1,fp); MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); } } } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairEffCut::write_restart_settings(FILE *fp) { fwrite(&cut_global,sizeof(double),1,fp); fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairEffCut::read_restart_settings(FILE *fp) { if (comm->me == 0) { fread(&cut_global,sizeof(double),1,fp); fread(&offset_flag,sizeof(int),1,fp); fread(&mix_flag,sizeof(int),1,fp); } MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } /* ---------------------------------------------------------------------- returns pointers to the log() of electron radius and corresponding force minimizer operates on log(radius) so radius never goes negative these arrays are stored locally by pair style ------------------------------------------------------------------------- */ void PairEffCut::min_xf_pointers(int ignore, double **xextra, double **fextra) { // grow arrays if necessary // need to be atom->nmax in length if (atom->nmax > nmax) { memory->destroy(min_eradius); memory->destroy(min_erforce); nmax = atom->nmax; memory->create(min_eradius,nmax,"pair:min_eradius"); memory->create(min_erforce,nmax,"pair:min_erforce"); } *xextra = min_eradius; *fextra = min_erforce; } /* ---------------------------------------------------------------------- minimizer requests the log() of electron radius and corresponding force calculate and store in min_eradius and min_erforce ------------------------------------------------------------------------- */ void PairEffCut::min_xf_get(int ignore) { double *eradius = atom->eradius; double *erforce = atom->erforce; int *spin = atom->spin; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (spin[i]) { min_eradius[i] = log(eradius[i]); min_erforce[i] = eradius[i]*erforce[i]; } else min_eradius[i] = min_erforce[i] = 0.0; } /* ---------------------------------------------------------------------- minimizer has changed the log() of electron radius propagate the change back to eradius ------------------------------------------------------------------------- */ void PairEffCut::min_x_set(int ignore) { double *eradius = atom->eradius; int *spin = atom->spin; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (spin[i]) eradius[i] = exp(min_eradius[i]); } /* ---------------------------------------------------------------------- memory usage of local atom-based arrays ------------------------------------------------------------------------- */ double PairEffCut::memory_usage() { double bytes = maxeatom * sizeof(double); bytes += maxvatom*6 * sizeof(double); bytes += 2 * nmax * sizeof(double); return bytes; }