diff --git a/src/USER-EFF/Install.sh b/src/USER-EFF/Install.sh index 64a7b71c0..407168c56 100644 --- a/src/USER-EFF/Install.sh +++ b/src/USER-EFF/Install.sh @@ -1,73 +1,73 @@ # Install/unInstall package files in LAMMPS if (test $1 = 1) then cp -p atom_vec_electron.cpp .. cp -p pair_eff_cut.cpp .. cp -p compute_ke_eff.cpp .. cp -p compute_ke_atom_eff.cpp .. cp -p compute_temp_deform_eff.cpp .. cp -p compute_temp_eff.cpp .. cp -p compute_temp_region_eff.cpp .. cp -p fix_langevin_eff.cpp .. cp -p fix_nh_eff.cpp .. cp -p fix_nve_eff.cpp .. cp -p fix_nvt_eff.cpp .. cp -p fix_nvt_sllod_eff.cpp .. cp -p fix_npt_eff.cpp .. cp -p fix_nph_eff.cpp .. cp -p fix_temp_rescale_eff.cpp .. cp -p atom_vec_electron.h .. cp -p pair_eff_cut.h .. cp -p pair_eff_inline.h .. cp -p compute_ke_eff.h .. cp -p compute_ke_atom_eff.h .. cp -p compute_temp_deform_eff.h .. cp -p compute_temp_eff.h .. cp -p compute_temp_region_eff.h .. cp -p fix_langevin_eff.h .. cp -p fix_nh_eff.h .. cp -p fix_nve_eff.h .. cp -p fix_nvt_eff.h .. cp -p fix_nvt_sllod_eff.h .. cp -p fix_npt_eff.h .. cp -p fix_nph_eff.h .. cp -p fix_temp_rescale_eff.h .. elif (test $1 = 0) then rm ../atom_vec_electron.cpp rm ../pair_eff_cut.cpp rm ../compute_ke_eff.cpp rm ../compute_ke_atom_eff.cpp 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 rm ../fix_nvt_sllod_eff.cpp rm ../fix_npt_eff.cpp rm ../fix_nph_eff.cpp rm ../fix_temp_rescale_eff.cpp rm ../atom_vec_electron.h rm ../pair_eff_cut.h rm ../pair_eff_inline.h rm ../compute_ke_eff.h rm ../compute_ke_atom_eff.h 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 rm ../fix_nvt_sllod_eff.h rm ../fix_npt_eff.h rm ../fix_nph_eff.h rm ../fix_temp_rescale_eff.h fi diff --git a/src/USER-EFF/README b/src/USER-EFF/README index 2e500bafc..c1c2453ec 100644 --- a/src/USER-EFF/README +++ b/src/USER-EFF/README @@ -1,59 +1,80 @@ The files in this directory are a user-contributed package for LAMMPS. The person who created these files is Andres Jaramillo-Botero at CalTech (ajaramil@wag.caltech.edu). Contact him directly if you have questions. -------------------------------------- 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) PACKAGE DESCRIPTION: 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. INSTALLATION: via a normal LAMMPS package installation: make yes-user-eff 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 20e33a091..f3a399768 100644 --- a/src/USER-EFF/atom_vec_electron.cpp +++ b/src/USER-EFF/atom_vec_electron.cpp @@ -1,801 +1,798 @@ /* ---------------------------------------------------------------------- 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 "lmptype.h" +#include "math.h" #include "stdlib.h" #include "atom_vec_electron.h" #include "atom.h" #include "domain.h" #include "modify.h" #include "force.h" #include "fix.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; #define DELTA 10000 /* ---------------------------------------------------------------------- */ AtomVecElectron::AtomVecElectron(LAMMPS *lmp, int narg, char **arg) : AtomVec(lmp, narg, arg) { comm_x_only = comm_f_only = 0; mass_type = 1; molecular = 0; size_forward = 5; size_reverse = 4; size_border = 10; size_velocity = 4; size_data_atom = 8; size_data_vel = 5; xcol_data = 6; atom->q_flag = atom->spin_flag = atom->eradius_flag = atom->ervel_flag = atom->erforce_flag = 1; } /* ---------------------------------------------------------------------- grow atom-electron arrays n = 0 grows arrays by DELTA n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ 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 *) memory->srealloc(atom->type,nmax*sizeof(int),"atom:type"); mask = atom->mask = (int *) memory->srealloc(atom->mask,nmax*sizeof(int),"atom:mask"); image = atom->image = (int *) memory->srealloc(atom->image,nmax*sizeof(int),"atom:image"); x = atom->x = memory->grow_2d_double_array(atom->x,nmax,3,"atom:x"); v = atom->v = memory->grow_2d_double_array(atom->v,nmax,3,"atom:v"); f = atom->f = memory->grow_2d_double_array(atom->f,nmax,3,"atom:f"); q = atom->q = (double *) memory->srealloc(atom->q,nmax*sizeof(double),"atom:q"); spin = atom->spin = (int *) memory->srealloc(atom->spin,nmax*sizeof(int),"atom:spin"); eradius = atom->eradius = (double *) memory->srealloc(atom->eradius,nmax*sizeof(double),"atom:eradius"); ervel = atom->ervel = (double *) memory->srealloc(atom->ervel,nmax*sizeof(double),"atom:ervel"); erforce = atom->erforce = (double *) memory->srealloc(atom->erforce,nmax*sizeof(double),"atom:erforce"); if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax); } /* ---------------------------------------------------------------------- reset local array ptrs ------------------------------------------------------------------------- */ 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; } /* ---------------------------------------------------------------------- */ void AtomVecElectron::copy(int i, int j) { tag[j] = tag[i]; type[j] = type[i]; mask[j] = mask[i]; image[j] = image[i]; x[j][0] = x[i][0]; x[j][1] = x[i][1]; x[j][2] = x[i][2]; v[j][0] = v[i][0]; v[j][1] = v[i][1]; v[j][2] = v[i][2]; q[j] = q[i]; spin[j] = spin[i]; eradius[j] = eradius[i]; ervel[j] = ervel[i]; if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j); } /* ---------------------------------------------------------------------- */ int AtomVecElectron::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) { int i,j,m; double dx,dy,dz; m = 0; if (pbc_flag == 0) { for (i = 0; i < n; i++) { j = list[i]; buf[m++] = x[j][0]; buf[m++] = x[j][1]; buf[m++] = x[j][2]; buf[m++] = eradius[j]; buf[m++] = ervel[j]; } } else { if (domain->triclinic == 0) { dx = pbc[0]*domain->xprd; dy = pbc[1]*domain->yprd; dz = pbc[2]*domain->zprd; } else { dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; dz = pbc[2]*domain->zprd; } for (i = 0; i < n; i++) { j = list[i]; buf[m++] = x[j][0] + dx; buf[m++] = x[j][1] + dy; buf[m++] = x[j][2] + dz; buf[m++] = eradius[j]; buf[m++] = ervel[j]; } } return m; } /* ---------------------------------------------------------------------- */ int AtomVecElectron::pack_comm_vel(int n, int *list, double *buf, int pbc_flag, int *pbc) { int i,j,m; double dx,dy,dz; m = 0; if (pbc_flag == 0) { for (i = 0; i < n; i++) { j = list[i]; buf[m++] = x[j][0]; buf[m++] = x[j][1]; buf[m++] = x[j][2]; buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; buf[m++] = eradius[j]; buf[m++] = ervel[j]; } } else { if (domain->triclinic == 0) { dx = pbc[0]*domain->xprd; dy = pbc[1]*domain->yprd; dz = pbc[2]*domain->zprd; } else { dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; dz = pbc[2]*domain->zprd; } for (i = 0; i < n; i++) { j = list[i]; buf[m++] = x[j][0] + dx; buf[m++] = x[j][1] + dy; buf[m++] = x[j][2] + dz; buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; buf[m++] = eradius[j]; buf[m++] = ervel[j]; } } return m; } /* ---------------------------------------------------------------------- */ void AtomVecElectron::unpack_comm(int n, int first, double *buf) { int i,m,last; m = 0; last = first + n; for (i = first; i < last; i++) { x[i][0] = buf[m++]; x[i][1] = buf[m++]; x[i][2] = buf[m++]; eradius[i] = buf[m++]; ervel[i] = buf[m++]; } } /* ---------------------------------------------------------------------- */ void AtomVecElectron::unpack_comm_vel(int n, int first, double *buf) { int i,m,last; m = 0; last = first + n; for (i = first; i < last; i++) { x[i][0] = buf[m++]; x[i][1] = buf[m++]; x[i][2] = buf[m++]; v[i][0] = buf[m++]; v[i][1] = buf[m++]; v[i][2] = buf[m++]; eradius[i] = buf[m++]; ervel[i] = buf[m++]; } } /* ---------------------------------------------------------------------- */ int AtomVecElectron::pack_comm_one(int i, double *buf) { buf[0] = eradius[i]; buf[1] = ervel[i]; return 2; } /* ---------------------------------------------------------------------- */ int AtomVecElectron::unpack_comm_one(int i, double *buf) { eradius[i] = buf[0]; ervel[i] = buf[1]; return 2; } /* ---------------------------------------------------------------------- */ int AtomVecElectron::pack_reverse(int n, int first, double *buf) { int i,m,last; m = 0; last = first + n; for (i = first; i < last; i++) { buf[m++] = f[i][0]; buf[m++] = f[i][1]; buf[m++] = f[i][2]; buf[m++] = erforce[i]; } return m; } /* ---------------------------------------------------------------------- */ int AtomVecElectron::pack_reverse_one(int i, double *buf) { buf[0] = erforce[i]; return 1; } /* ---------------------------------------------------------------------- */ void AtomVecElectron::unpack_reverse(int n, int *list, double *buf) { int i,j,m; m = 0; for (i = 0; i < n; i++) { j = list[i]; f[j][0] += buf[m++]; f[j][1] += buf[m++]; f[j][2] += buf[m++]; erforce[j] += buf[m++]; } } /* ---------------------------------------------------------------------- */ int AtomVecElectron::unpack_reverse_one(int i, double *buf) { erforce[i] += buf[0]; return 1; } /* ---------------------------------------------------------------------- */ int AtomVecElectron::pack_border_vel(int n, int *list, double *buf, int pbc_flag, int *pbc) { int i,j,m; double dx,dy,dz; m = 0; if (pbc_flag == 0) { for (i = 0; i < n; i++) { j = list[i]; buf[m++] = x[j][0]; buf[m++] = x[j][1]; buf[m++] = x[j][2]; buf[m++] = tag[j]; buf[m++] = type[j]; buf[m++] = mask[j]; buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; buf[m++] = q[j]; buf[m++] = spin[j]; buf[m++] = eradius[j]; buf[m++] = ervel[j]; } } else { if (domain->triclinic == 0) { dx = pbc[0]*domain->xprd; dy = pbc[1]*domain->yprd; dz = pbc[2]*domain->zprd; } else { dx = pbc[0]; dy = pbc[1]; dz = pbc[2]; } for (i = 0; i < n; i++) { j = list[i]; buf[m++] = x[j][0] + dx; buf[m++] = x[j][1] + dy; buf[m++] = x[j][2] + dz; buf[m++] = tag[j]; buf[m++] = type[j]; buf[m++] = mask[j]; buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; buf[m++] = q[j]; buf[m++] = spin[j]; buf[m++] = eradius[j]; buf[m++] = ervel[j]; } } return m; } /* ---------------------------------------------------------------------- */ void AtomVecElectron::unpack_border_vel(int n, int first, double *buf) { int i,m,last; m = 0; last = first + n; for (i = first; i < last; i++) { if (i == nmax) grow(0); x[i][0] = buf[m++]; x[i][1] = buf[m++]; x[i][2] = buf[m++]; tag[i] = static_cast (buf[m++]); type[i] = static_cast (buf[m++]); mask[i] = static_cast (buf[m++]); v[i][0] = buf[m++]; v[i][1] = buf[m++]; v[i][2] = buf[m++]; q[i] = buf[m++]; spin[i] = static_cast (buf[m++]); eradius[i] = buf[m++]; ervel[i] = buf[m++]; } } /* ---------------------------------------------------------------------- */ int AtomVecElectron::pack_border(int n, int *list, double *buf, int pbc_flag, int *pbc) { int i,j,m; double dx,dy,dz; m = 0; if (pbc_flag == 0) { for (i = 0; i < n; i++) { j = list[i]; buf[m++] = x[j][0]; buf[m++] = x[j][1]; buf[m++] = x[j][2]; buf[m++] = tag[j]; buf[m++] = type[j]; buf[m++] = mask[j]; buf[m++] = q[j]; buf[m++] = spin[j]; buf[m++] = eradius[j]; buf[m++] = ervel[j]; } } else { if (domain->triclinic == 0) { dx = pbc[0]*domain->xprd; dy = pbc[1]*domain->yprd; dz = pbc[2]*domain->zprd; } else { dx = pbc[0]; dy = pbc[1]; dz = pbc[2]; } for (i = 0; i < n; i++) { j = list[i]; buf[m++] = x[j][0] + dx; buf[m++] = x[j][1] + dy; buf[m++] = x[j][2] + dz; buf[m++] = tag[j]; buf[m++] = type[j]; buf[m++] = mask[j]; buf[m++] = q[j]; buf[m++] = spin[j]; buf[m++] = eradius[j]; buf[m++] = ervel[j]; } } return m; } /* ---------------------------------------------------------------------- */ int AtomVecElectron::pack_border_one(int i, double *buf) { buf[0] = q[i]; buf[1] = spin[i]; buf[2] = eradius[i]; buf[3] = ervel[i]; return 4; } /* ---------------------------------------------------------------------- */ void AtomVecElectron::unpack_border(int n, int first, double *buf) { int i,m,last; m = 0; last = first + n; for (i = first; i < last; i++) { if (i == nmax) grow(0); x[i][0] = buf[m++]; x[i][1] = buf[m++]; x[i][2] = buf[m++]; tag[i] = static_cast (buf[m++]); type[i] = static_cast (buf[m++]); mask[i] = static_cast (buf[m++]); q[i] = buf[m++]; spin[i] = static_cast (buf[m++]); eradius[i] = buf[m++]; ervel[i] = buf[m++]; } } /* ---------------------------------------------------------------------- */ int AtomVecElectron::unpack_border_one(int i, double *buf) { q[i] = buf[0]; spin[i] = static_cast (buf[1]); eradius[i] = buf[2]; ervel[i] = buf[3]; return 4; } /* ---------------------------------------------------------------------- pack data for atom I for sending to another proc xyz must be 1st 3 values, so comm::exchange() can test on them ------------------------------------------------------------------------- */ int AtomVecElectron::pack_exchange(int i, double *buf) { int m = 1; buf[m++] = x[i][0]; buf[m++] = x[i][1]; buf[m++] = x[i][2]; buf[m++] = v[i][0]; buf[m++] = v[i][1]; buf[m++] = v[i][2]; buf[m++] = tag[i]; buf[m++] = type[i]; buf[m++] = mask[i]; buf[m++] = image[i]; buf[m++] = q[i]; buf[m++] = spin[i]; buf[m++] = eradius[i]; buf[m++] = ervel[i]; if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]); buf[0] = m; return m; } /* ---------------------------------------------------------------------- */ int AtomVecElectron::unpack_exchange(double *buf) { int nlocal = atom->nlocal; if (nlocal == nmax) grow(0); int m = 1; x[nlocal][0] = buf[m++]; x[nlocal][1] = buf[m++]; x[nlocal][2] = buf[m++]; v[nlocal][0] = buf[m++]; v[nlocal][1] = buf[m++]; v[nlocal][2] = buf[m++]; tag[nlocal] = static_cast (buf[m++]); type[nlocal] = static_cast (buf[m++]); mask[nlocal] = static_cast (buf[m++]); image[nlocal] = static_cast (buf[m++]); q[nlocal] = buf[m++]; spin[nlocal] = static_cast (buf[m++]); eradius[nlocal] = buf[m++]; ervel[nlocal] = buf[m++]; if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) m += modify->fix[atom->extra_grow[iextra]]-> unpack_exchange(nlocal,&buf[m]); atom->nlocal++; return m; } /* ---------------------------------------------------------------------- size of restart data for all atoms owned by this proc include extra data stored by fixes ------------------------------------------------------------------------- */ int AtomVecElectron::size_restart() { int i; int nlocal = atom->nlocal; int n = 15 * nlocal; // Associated with pack_restart if (atom->nextra_restart) for (int iextra = 0; iextra < atom->nextra_restart; iextra++) for (i = 0; i < nlocal; i++) n += modify->fix[atom->extra_restart[iextra]]->size_restart(i); return n; } /* ---------------------------------------------------------------------- pack atom I's data for restart file including extra quantities xyz must be 1st 3 values, so that read_restart can test on them molecular types may be negative, but write as positive ------------------------------------------------------------------------- */ int AtomVecElectron::pack_restart(int i, double *buf) { int m = 1; buf[m++] = x[i][0]; buf[m++] = x[i][1]; buf[m++] = x[i][2]; buf[m++] = tag[i]; buf[m++] = type[i]; buf[m++] = mask[i]; buf[m++] = image[i]; buf[m++] = v[i][0]; buf[m++] = v[i][1]; buf[m++] = v[i][2]; buf[m++] = q[i]; buf[m++] = spin[i]; buf[m++] = eradius[i]; buf[m++] = ervel[i]; if (atom->nextra_restart) for (int iextra = 0; iextra < atom->nextra_restart; iextra++) m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]); buf[0] = m; return m; } /* ---------------------------------------------------------------------- unpack data for one atom from restart file including extra quantities ------------------------------------------------------------------------- */ int AtomVecElectron::unpack_restart(double *buf) { int nlocal = atom->nlocal; if (nlocal == nmax) { grow(0); if (atom->nextra_store) atom->extra = memory->grow_2d_double_array(atom->extra,nmax, atom->nextra_store, "atom:extra"); } int m = 1; x[nlocal][0] = buf[m++]; x[nlocal][1] = buf[m++]; x[nlocal][2] = buf[m++]; tag[nlocal] = static_cast (buf[m++]); type[nlocal] = static_cast (buf[m++]); mask[nlocal] = static_cast (buf[m++]); image[nlocal] = static_cast (buf[m++]); v[nlocal][0] = buf[m++]; v[nlocal][1] = buf[m++]; v[nlocal][2] = buf[m++]; q[nlocal] = buf[m++]; spin[nlocal] = static_cast (buf[m++]); eradius[nlocal] = buf[m++]; ervel[nlocal] = buf[m++]; double **extra = atom->extra; if (atom->nextra_store) { int size = static_cast (buf[0]) - m; for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++]; } atom->nlocal++; return m; } /* ---------------------------------------------------------------------- create one atom of itype at coord set other values to defaults ------------------------------------------------------------------------- */ void AtomVecElectron::create_atom(int itype, double *coord) { int nlocal = atom->nlocal; if (nlocal == nmax) grow(0); tag[nlocal] = 0; type[nlocal] = itype; x[nlocal][0] = coord[0]; x[nlocal][1] = coord[1]; x[nlocal][2] = coord[2]; mask[nlocal] = 1; image[nlocal] = (512 << 20) | (512 << 10) | 512; v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; v[nlocal][2] = 0.0; q[nlocal] = 0.0; spin[nlocal] = 1; eradius[nlocal] = 1.0; ervel[nlocal] = 0.0; atom->nlocal++; } /* ---------------------------------------------------------------------- unpack one line from Atoms section of data file initialize other atom quantities ------------------------------------------------------------------------- */ void AtomVecElectron::data_atom(double *coord, int imagetmp, char **values) { int nlocal = atom->nlocal; if (nlocal == nmax) grow(0); tag[nlocal] = atoi(values[0]); if (tag[nlocal] <= 0) error->one("Invalid atom ID in Atoms section of data file"); type[nlocal] = atoi(values[1]); if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes) error->one("Invalid atom type in Atoms section of data file"); q[nlocal] = atof(values[2]); spin[nlocal] = atoi(values[3]); eradius[nlocal] = atof(values[4]); x[nlocal][0] = coord[0]; x[nlocal][1] = coord[1]; x[nlocal][2] = coord[2]; image[nlocal] = imagetmp; mask[nlocal] = 1; v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; v[nlocal][2] = 0.0; ervel[nlocal] = 0.0; atom->nlocal++; } /* ---------------------------------------------------------------------- unpack hybrid quantities from one line in Atoms section of data file initialize other atom quantities for this sub-style ------------------------------------------------------------------------- */ int AtomVecElectron::data_atom_hybrid(int nlocal, char **values) { q[nlocal] = atof(values[0]); spin[nlocal] = atoi(values[1]); eradius[nlocal] = atof(values[2]); if (eradius[nlocal] < 0.0) error->one("Invalid eradius in Atoms section of data file"); v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; v[nlocal][2] = 0.0; ervel[nlocal] = 0.0; return 2; } /* ---------------------------------------------------------------------- unpack one line from Velocities section of data file ------------------------------------------------------------------------- */ void AtomVecElectron::data_vel(int m, char **values) { v[m][0] = atof(values[0]); v[m][1] = atof(values[1]); v[m][2] = atof(values[2]); ervel[m] = atof(values[3]); } /* ---------------------------------------------------------------------- unpack hybrid quantities from one line in Velocities section of data file ------------------------------------------------------------------------- */ int AtomVecElectron::data_vel_hybrid(int m, char **values) { ervel[m] = atof(values[0]); return 1; } /* ---------------------------------------------------------------------- return # of bytes of allocated memory ------------------------------------------------------------------------- */ double AtomVecElectron::memory_usage() { double bytes = 0.0; if (atom->memcheck("tag")) bytes += nmax * sizeof(int); if (atom->memcheck("type")) bytes += nmax * sizeof(int); if (atom->memcheck("mask")) bytes += nmax * sizeof(int); if (atom->memcheck("image")) bytes += nmax * sizeof(int); if (atom->memcheck("x")) bytes += nmax*3 * sizeof(double); if (atom->memcheck("v")) bytes += nmax*3 * sizeof(double); if (atom->memcheck("f")) bytes += nmax*3 * sizeof(double); if (atom->memcheck("q")) bytes += nmax * sizeof(double); if (atom->memcheck("spin")) bytes += nmax * sizeof(int); if (atom->memcheck("eradius")) bytes += nmax * sizeof(double); if (atom->memcheck("ervel")) bytes += nmax * sizeof(double); if (atom->memcheck("erforce")) bytes += nmax * sizeof(double); return bytes; } diff --git a/src/USER-EFF/compute_ke_atom_eff.cpp b/src/USER-EFF/compute_ke_atom_eff.cpp index f9dba6e5d..73e63313a 100644 --- a/src/USER-EFF/compute_ke_atom_eff.cpp +++ b/src/USER-EFF/compute_ke_atom_eff.cpp @@ -1,124 +1,113 @@ /* ---------------------------------------------------------------------- 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) + Contributing author: Andres Jaramillo-Botero ------------------------------------------------------------------------- */ #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" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeKEAtomEff::ComputeKEAtomEff(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 3) error->all("Illegal compute ke/atom/eff command"); peratom_flag = 1; size_peratom_cols = 0; nmax = 0; ke = NULL; // error check if (!atom->ervel_flag) error->all("Compute ke/atom/eff requires atom attribute ervel"); } /* ---------------------------------------------------------------------- */ ComputeKEAtomEff::~ComputeKEAtomEff() { memory->sfree(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("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->sfree(ke); nmax = atom->nmax; ke = (double *) memory->smalloc(nmax*sizeof(double),"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; - 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; } } /* ---------------------------------------------------------------------- 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 600f0ebae..08759c737 100644 --- a/src/USER-EFF/compute_ke_eff.cpp +++ b/src/USER-EFF/compute_ke_eff.cpp @@ -1,88 +1,80 @@ /* ---------------------------------------------------------------------- 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) + Contributing author: Andres Jaramillo-Botero ------------------------------------------------------------------------- */ #include "mpi.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("Illegal compute ke/eff command"); scalar_flag = 1; extscalar = 1; // error check if (!atom->ervel_flag) error->all("Compute ke/eff requires atom attribute ervel"); } /* ---------------------------------------------------------------------- */ 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 *rmass = atom->rmass; double *mass = atom->mass; int *spin = atom->spin; int *mask = atom->mask; int *type = atom->type; int nlocal = atom->nlocal; 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]; } } 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 ab6b51ea4..5a093b65e 100644 --- a/src/USER-EFF/compute_temp_deform_eff.cpp +++ b/src/USER-EFF/compute_temp_deform_eff.cpp @@ -1,178 +1,316 @@ /* ---------------------------------------------------------------------- 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)) + 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 /* ---------------------------------------------------------------------- */ 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"); } /* ---------------------------------------------------------------------- */ 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 (spin[i]) one++; + if (abs(spin[i])==1) one++; } int nelectrons; MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world); 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; - double *rmass = atom->rmass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; // 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]; + 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]; } } 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; - double *rmass = atom->rmass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; 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]; + 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]; } } 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 32bf165c4..47f440942 100644 --- a/src/USER-EFF/compute_temp_deform_eff.h +++ b/src/USER-EFF/compute_temp_deform_eff.h @@ -1,41 +1,54 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifdef COMPUTE_CLASS ComputeStyle(temp/deform/eff,ComputeTempDeformEff) #else #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(); - - private: - void dof_compute(); + virtual ~ComputeTempDeformEff(); + void init(); + virtual double compute_scalar(); + virtual void compute_vector(); + + 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(); }; } #endif #endif diff --git a/src/USER-EFF/compute_temp_eff.cpp b/src/USER-EFF/compute_temp_eff.cpp index 26f6539e6..d5e359b94 100644 --- a/src/USER-EFF/compute_temp_eff.cpp +++ b/src/USER-EFF/compute_temp_eff.cpp @@ -1,152 +1,165 @@ /* ---------------------------------------------------------------------- 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 "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) : - 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(); } /* ---------------------------------------------------------------------- */ 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 (spin[i]) one++; + if (abs(spin[i])==1) one++; } int nelectrons; MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world); // average over nuclear dof only 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; - 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]; } } } 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; - 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]; } } 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_eff.h b/src/USER-EFF/compute_temp_eff.h index 683b73a29..197dd4dbc 100644 --- a/src/USER-EFF/compute_temp_eff.h +++ b/src/USER-EFF/compute_temp_eff.h @@ -1,41 +1,46 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifdef COMPUTE_CLASS ComputeStyle(temp/eff,ComputeTempEff) #else #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(); }; } #endif #endif diff --git a/src/USER-EFF/compute_temp_region_eff.cpp b/src/USER-EFF/compute_temp_region_eff.cpp index be6bee419..7af81a42c 100644 --- a/src/USER-EFF/compute_temp_region_eff.cpp +++ b/src/USER-EFF/compute_temp_region_eff.cpp @@ -1,145 +1,276 @@ /* ---------------------------------------------------------------------- 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 "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) : - 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; } /* ---------------------------------------------------------------------- */ double ComputeTempRegionEff::compute_scalar() { invoked_scalar = update->ntimestep; double **x = atom->x; 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; Region *region = domain->regions[iregion]; int count = 0; double t = 0.0; - if (rmass) { - 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 { + 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++; + 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]; } } double tarray[2],tarray_all[2]; tarray[0] = count; 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 (spin[i]) one++; + if (abs(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; - double *rmass = atom->rmass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; 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])) { - 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]; - } + } } 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 e2d7add69..3340fe826 100644 --- a/src/USER-EFF/compute_temp_region_eff.h +++ b/src/USER-EFF/compute_temp_region_eff.h @@ -1,38 +1,50 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifdef COMPUTE_CLASS ComputeStyle(temp/region/eff,ComputeTempRegionEff) #else #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; }; } #endif #endif diff --git a/src/USER-EFF/fix_langevin_eff.cpp b/src/USER-EFF/fix_langevin_eff.cpp index 7aae42b52..45b9eba0f 100644 --- a/src/USER-EFF/fix_langevin_eff.cpp +++ b/src/USER-EFF/fix_langevin_eff.cpp @@ -1,338 +1,246 @@ /* ---------------------------------------------------------------------- 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 "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 "random_mars.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; enum{NOBIAS,BIAS}; /* ---------------------------------------------------------------------- */ FixLangevinEff::FixLangevinEff(LAMMPS *lmp, int narg, char **arg) : FixLangevin(lmp, narg, arg) { erforcelangevin = NULL; } /* ---------------------------------------------------------------------- */ FixLangevinEff::~FixLangevinEff() { memory->sfree(erforcelangevin); } /* ---------------------------------------------------------------------- */ void FixLangevinEff::post_force_no_tally() { double gamma1,gamma2; 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; double t_target = t_start + delta * (t_stop-t_start); double tsqrt = sqrt(t_target); // 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 - 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]); } } } } /* ---------------------------------------------------------------------- */ void FixLangevinEff::post_force_tally() { double gamma1,gamma2; // reallocate flangevin if necessary if (atom->nmax > nmax) { memory->destroy_2d_double_array(flangevin); memory->sfree(erforcelangevin); nmax = atom->nmax; flangevin = memory->create_2d_double_array(nmax,3,"langevin:flangevin"); erforcelangevin = (double *) memory->smalloc(nmax*sizeof(double),"langevin/eff:erforcelangevin"); } 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; double t_target = t_start + delta * (t_stop-t_start); double tsqrt = sqrt(t_target); // 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 - 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]); } } } } /* ---------------------------------------------------------------------- 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) + 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; } /* ---------------------------------------------------------------------- */ double FixLangevinEff::compute_scalar() { if (!tally) 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) + 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; } 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 0a30b02ee..9dafafb24 100644 --- a/src/USER-EFF/fix_nh_eff.cpp +++ b/src/USER-EFF/fix_nh_eff.cpp @@ -1,126 +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" using namespace LAMMPS_NS; enum{NOBIAS,BIAS}; /* ---------------------------------------------------------------------- */ 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"); } /* ---------------------------------------------------------------------- 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 *rmass = atom->rmass; double *mass = atom->mass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; 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; } } } } /* ---------------------------------------------------------------------- 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 (spin[i]) eradius[i] += dtv * ervel[i]; + if (abs(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 (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 04d265e7b..9862bf8d9 100644 --- a/src/USER-EFF/fix_nph_eff.h +++ b/src/USER-EFF/fix_nph_eff.h @@ -1,36 +1,36 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifdef FIX_CLASS -FixStyle(nph/sphere,FixNPHEff) +FixStyle(nph/eff,FixNPHEff) #else #ifndef LMP_FIX_NPH_EFF_H #define LMP_FIX_NPH_EFF_H #include "fix_nh_eff.h" namespace LAMMPS_NS { class FixNPHEff : public FixNHEff { public: FixNPHEff(class LAMMPS *, int, char **); ~FixNPHEff() {} }; } #endif #endif diff --git a/src/USER-EFF/fix_npt_eff.cpp b/src/USER-EFF/fix_npt_eff.cpp index 4c2b381fa..d71a25d46 100644 --- a/src/USER-EFF/fix_npt_eff.cpp +++ b/src/USER-EFF/fix_npt_eff.cpp @@ -1,78 +1,67 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #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; /* ---------------------------------------------------------------------- */ FixNPTEff::FixNPTEff(LAMMPS *lmp, int narg, char **arg) : FixNHEff(lmp, narg, arg) { if (!tstat_flag) error->all("Temperature control must be used with fix npt/eff"); if (!pstat_flag) error->all("Pressure control must be used with fix npt/eff"); // create a new compute temp style // id = fix-ID + temp // compute group = all since pressure is always global (group all) // and thus its KE/temperature contribution should use group all int n = strlen(id) + 6; id_temp = new char[n]; strcpy(id_temp,id); strcat(id_temp,"_temp"); char **newarg = new char*[3]; newarg[0] = id_temp; newarg[1] = (char *) "all"; newarg[2] = (char *) "temp/eff"; modify->add_compute(3,newarg); delete [] newarg; tflag = 1; // create a new compute pressure style // id = fix-ID + press, compute group = all // pass id_temp as 4th arg to pressure constructor n = strlen(id) + 7; id_press = new char[n]; strcpy(id_press,id); strcat(id_press,"_press"); newarg = new char*[4]; newarg[0] = id_press; newarg[1] = (char *) "all"; newarg[2] = (char *) "pressure"; newarg[3] = id_temp; modify->add_compute(4,newarg); delete [] newarg; pflag = 1; } diff --git a/src/USER-EFF/fix_nve_eff.cpp b/src/USER-EFF/fix_nve_eff.cpp index 0f20bec33..99866afe0 100644 --- a/src/USER-EFF/fix_nve_eff.cpp +++ b/src/USER-EFF/fix_nve_eff.cpp @@ -1,149 +1,172 @@ /* ---------------------------------------------------------------------- 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 "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" 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; } /* ---------------------------------------------------------------------- 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; - double *rmass = atom->rmass; int *spin = atom->spin; 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 (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]; + } } } } } /* ---------------------------------------------------------------------- */ 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; - double *rmass = atom->rmass; int *spin = atom->spin; 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 (rmass) { - 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 { + 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 (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 33d437961..68be95295 100644 --- a/src/USER-EFF/fix_nve_eff.h +++ b/src/USER-EFF/fix_nve_eff.h @@ -1,37 +1,47 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifdef FIX_CLASS FixStyle(nve/eff,FixNVEEff) #else #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; }; } #endif #endif diff --git a/src/USER-EFF/fix_nvt_eff.h b/src/USER-EFF/fix_nvt_eff.h index 3e832094c..139febcf3 100644 --- a/src/USER-EFF/fix_nvt_eff.h +++ b/src/USER-EFF/fix_nvt_eff.h @@ -1,36 +1,36 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifdef FIX_CLASS FixStyle(nvt/eff,FixNVTEff) #else #ifndef LMP_FIX_NVT_EFF_H #define LMP_FIX_NVT_EFF_H #include "fix_nh_eff.h" namespace LAMMPS_NS { class FixNVTEff : public FixNHEff { public: FixNVTEff(class LAMMPS *, int, char **); - virtual ~FixNVTEff() {} + ~FixNVTEff() {} }; } #endif #endif diff --git a/src/USER-EFF/fix_nvt_sllod_eff.cpp b/src/USER-EFF/fix_nvt_sllod_eff.cpp index 792208cdc..12bdac820 100644 --- a/src/USER-EFF/fix_nvt_sllod_eff.cpp +++ b/src/USER-EFF/fix_nvt_sllod_eff.cpp @@ -1,127 +1,129 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #include "math.h" -#include "string.h" -#include "stdlib.h" #include "fix_nvt_sllod_eff.h" #include "math_extra.h" #include "atom.h" #include "domain.h" #include "group.h" #include "modify.h" #include "fix.h" #include "fix_deform.h" #include "compute.h" #include "error.h" using namespace LAMMPS_NS; 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 int n = strlen(id) + 6; 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]; newarg[2] = (char *) "temp/deform/eff"; modify->add_compute(3,newarg); delete [] newarg; tflag = 1; } /* ---------------------------------------------------------------------- */ void FixNVTSllodEff::init() { - FixNVTEff::init(); + FixNHEff::init(); if (!temperature->tempbias) error->all("Temperature for fix nvt/sllod/eff does not have a bias"); nondeformbias = 0; if (strcmp(temperature->style,"temp/deform/eff") != 0) nondeformbias = 1; // check fix deform remap settings int i; for (i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"deform") == 0) { if (((FixDeform *) modify->fix[i])->remapflag != V_REMAP) error->all("Using fix nvt/sllod/eff with inconsistent fix deform " "remap option"); break; } if (i == modify->nfix) error->all("Using fix nvt/sllod/eff with no fix deform defined"); } /* ---------------------------------------------------------------------- perform half-step scaling of velocities -----------------------------------------------------------------------*/ void FixNVTSllodEff::nh_v_temp() { // remove and restore bias = streaming velocity = Hrate*lamda + Hratelo // thermostat thermal velocity only // vdelu = SLLOD correction = Hrate*Hinv*vthermal // for non temp/deform BIAS: // calculate temperature since some computes require temp // computed on current nlocal atoms to remove bias if (nondeformbias) double tmp = temperature->compute_scalar(); double **v = atom->v; double *ervel = atom->ervel; int *spin = atom->spin; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; double h_two[6],vdelu[3]; MathExtra::multiply_shape_shape(domain->h_rate,domain->h_inv,h_two); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); vdelu[0] = h_two[0]*v[i][0] + h_two[5]*v[i][1] + h_two[4]*v[i][2]; vdelu[1] = h_two[1]*v[i][1] + h_two[3]*v[i][2]; vdelu[2] = h_two[2]*v[i][2]; v[i][0] = v[i][0]*factor_eta - dthalf*vdelu[0]; 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 e8f9007cf..852ffb89c 100644 --- a/src/USER-EFF/fix_nvt_sllod_eff.h +++ b/src/USER-EFF/fix_nvt_sllod_eff.h @@ -1,42 +1,42 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifdef FIX_CLASS FixStyle(nvt/sllod/eff,FixNVTSllodEff) #else #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() {} void init(); private: int nondeformbias; void nh_v_temp(); }; } #endif #endif diff --git a/src/USER-EFF/fix_temp_rescale_eff.cpp b/src/USER-EFF/fix_temp_rescale_eff.cpp index 9a31b1023..e5a0d19be 100644 --- a/src/USER-EFF/fix_temp_rescale_eff.cpp +++ b/src/USER-EFF/fix_temp_rescale_eff.cpp @@ -1,109 +1,206 @@ /* ---------------------------------------------------------------------- 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 "string.h" #include "stdlib.h" #include "math.h" #include "fix_temp_rescale_eff.h" #include "atom.h" #include "force.h" #include "group.h" #include "update.h" #include "domain.h" #include "region.h" +#include "comm.h" #include "modify.h" #include "compute.h" #include "error.h" using namespace LAMMPS_NS; 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; newarg[1] = group->names[igroup]; 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; } /* ---------------------------------------------------------------------- */ void FixTempRescaleEff::end_of_step() { double t_current = temperature->compute_scalar(); if (t_current == 0.0) error->all("Computed temperature for fix temp/rescale/eff cannot be 0.0"); double delta = update->ntimestep - update->beginstep; delta /= update->endstep - update->beginstep; double t_target = t_start + delta * (t_stop-t_start); // rescale velocity of appropriate atoms if outside window // for BIAS: // temperature is current, so do not need to re-compute // OK to not test returned v = 0, since factor is multiplied by v if (fabs(t_current-t_target) > t_window) { t_target = t_current - fraction*(t_current-t_target); double factor = sqrt(t_target/t_current); double efactor = 0.5 * force->boltz * temperature->dof; double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; int *spin = atom->spin; double *ervel = atom->ervel; if (which == NOBIAS) { energy += (t_current-t_target) * efactor; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { v[i][0] *= factor; v[i][1] *= factor; v[i][2] *= factor; - if (spin[i]) + if (abs(spin[i])==1) ervel[i] *= factor; } } } else { energy += (t_current-t_target) * efactor; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); 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]); } } } } } + +/* ---------------------------------------------------------------------- */ + +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 76e7700a6..9ac4c003a 100644 --- a/src/USER-EFF/fix_temp_rescale_eff.h +++ b/src/USER-EFF/fix_temp_rescale_eff.h @@ -1,37 +1,51 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifdef FIX_CLASS FixStyle(temp/rescale/eff,FixTempRescaleEff) #else #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; }; } #endif #endif diff --git a/src/USER-EFF/pair_eff_cut.cpp b/src/USER-EFF/pair_eff_cut.cpp index 013f8ecec..9b196e98e 100644 --- a/src/USER-EFF/pair_eff_cut.cpp +++ b/src/USER-EFF/pair_eff_cut.cpp @@ -1,715 +1,960 @@ /* ---------------------------------------------------------------------- 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 authors: Andres Jaramillo-Botero and Julius Su (Caltech) + 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; #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b)) /* ---------------------------------------------------------------------- */ 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->sfree(min_eradius); memory->sfree(min_erforce); if (allocated) { memory->destroy_2d_int_array(setflag); memory->destroy_2d_double_array(cutsq); memory->destroy_2d_double_array(cut); } } /* ---------------------------------------------------------------------- */ 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; 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 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; firstneigh = list->firstneigh; // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; - qtmp = q[i]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; 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); + double qxq = q[i]*q[j]; + + ElecNucNuc(qxq, rc, &ecoul, &fpair); } - - // I is nucleus, J is electron - 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; + // 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]; 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; - 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(); } } /* ---------------------------------------------------------------------- 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 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; } } } } } /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ void PairEffCut::allocate() { allocated = 1; int n = atom->ntypes; setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); cut = memory->create_2d_double_array(n+1,n+1,"pair:cut"); } /* --------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ 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]); 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("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 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("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("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("Pair eff/cut requires atom attributes " "q, spin, eradius, erforce"); if (comm->ghost_velocity == 0) error->all("Pair eff/cut requires ghost atoms store velocity"); // 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); } /* ---------------------------------------------------------------------- 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->sfree(min_eradius); memory->sfree(min_erforce); nmax = atom->nmax; min_eradius = (double *) memory->smalloc(nmax*sizeof(double), "pair:min_eradius"); min_erforce = (double *) memory->smalloc(nmax*sizeof(double), "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; } diff --git a/src/USER-EFF/pair_eff_cut.h b/src/USER-EFF/pair_eff_cut.h index e47a13f82..7519d7335 100644 --- a/src/USER-EFF/pair_eff_cut.h +++ b/src/USER-EFF/pair_eff_cut.h @@ -1,63 +1,65 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifdef PAIR_CLASS PairStyle(eff/cut,PairEffCut) #else #ifndef LMP_PAIR_EFF_CUT_H #define LMP_PAIR_EFF_CUT_H #include "pair.h" namespace LAMMPS_NS { class PairEffCut : public Pair { public: PairEffCut(class LAMMPS *); virtual ~PairEffCut(); virtual void compute(int, int); virtual void settings(int, char **); void coeff(int, char **); void init_style(); void min_pointers(double **, double **); double init_one(int, int); void write_restart(FILE *); void read_restart(FILE *); virtual void write_restart_settings(FILE *); virtual void read_restart_settings(FILE *); void min_xf_pointers(int, double **, double **); void min_xf_get(int); void min_x_set(int); double memory_usage(); private: 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; void allocate(); void virial_eff_compute(); void ev_tally_eff(int, int, int, int, double, double); }; } #endif #endif diff --git a/src/USER-EFF/pair_eff_inline.h b/src/USER-EFF/pair_eff_inline.h index b21e35188..6f7705b19 100644 --- a/src/USER-EFF/pair_eff_inline.h +++ b/src/USER-EFF/pair_eff_inline.h @@ -1,446 +1,554 @@ /* ---------------------------------------------------------------------- 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 authors: Andres Jaramillo-Botero and Julius Su (Caltech) + Contributing authors: Andres Jaramillo-Botero, Hai Xiao, Julius Su (Caltech) ------------------------------------------------------------------------- */ namespace LAMMPS_NS { #define PAULI_RE 0.9 #define PAULI_RC 1.125 #define PAULI_RHO -0.2 #define ERF_TERMS1 12 #define ERF_TERMS2 7 #define DERF_TERMS 13 // error arrays double E1[] = { 1.483110564084803581889448079057, -3.01071073386594942470731046311E-1, 6.8994830689831566246603180718E-2, -1.3916271264722187682546525687E-2, 2.420799522433463662891678239E-3, -3.65863968584808644649382577E-4, 4.8620984432319048282887568E-5, -5.749256558035684835054215E-6, 6.11324357843476469706758E-7, -5.8991015312958434390846E-8, 5.207009092068648240455E-9, -4.23297587996554326810E-10, 3.1881135066491749748E-11, -2.236155018832684273E-12, 1.46732984799108492E-13, -9.044001985381747E-15, 5.25481371547092E-16, -2.8874261222849E-17, 1.504785187558E-18, -7.4572892821E-20, 3.522563810E-21, -1.58944644E-22, 6.864365E-24, -2.84257E-25, 1.1306E-26, -4.33E-28, 1.6E-29, -1.0E-30 }; double E2[] = { 1.077977852072383151168335910348, -2.6559890409148673372146500904E-2, -1.487073146698099509605046333E-3, -1.38040145414143859607708920E-4, -1.1280303332287491498507366E-5, -1.172869842743725224053739E-6, -1.03476150393304615537382E-7, -1.1899114085892438254447E-8, -1.016222544989498640476E-9, -1.37895716146965692169E-10, -9.369613033737303335E-12, -1.918809583959525349E-12, -3.7573017201993707E-14, -3.7053726026983357E-14, 2.627565423490371E-15, -1.121322876437933E-15, 1.84136028922538E-16, -4.9130256574886E-17, 1.0704455167373E-17, -2.671893662405E-18, 6.49326867976E-19, -1.65399353183E-19, 4.2605626604E-20, -1.1255840765E-20, 3.025617448E-21, -8.29042146E-22, 2.31049558E-22, -6.5469511E-23, 1.8842314E-23, -5.504341E-24, 1.630950E-24, -4.89860E-25, 1.49054E-25, -4.5922E-26, 1.4318E-26, -4.516E-27, 1.440E-27, -4.64E-28, 1.51E-28, -5.0E-29, 1.7E-29, -6.0E-30, 2.0E-30, -1.0E-30 }; double DE1[] = { -0.689379974848418501361491576718, 0.295939056851161774752959335568, -0.087237828075228616420029484096, 0.019959734091835509766546612696, -0.003740200486895490324750329974, 0.000593337912367800463413186784, -0.000081560801047403878256504204, 9.886099179971884018535968E-6, -1.071209234904290565745194E-6, 1.0490945447626050322784E-7, -9.370959271038746709966E-9, 7.6927263488753841874E-10, -5.8412335114551520146E-11, 4.125393291736424788E-12, -2.72304624901729048E-13, 1.6869717361387012E-14, -9.84565340276638E-16, 5.4313471880068E-17, -2.840458699772E-18, 1.4120512798E-19, -6.688772574E-21, 3.0257558E-22, -1.3097526E-23, 5.4352E-25, -2.1704E-26, 8.32E-28, -5.4E-29 }; double DE2[] = { 0.717710208167480928473053690384, -0.379868973985143305103199928808, 0.125832094465157378967135019248, -0.030917661684228839423081992424, 0.006073689914144320367855343072, -0.000996057789064916825079352632, 0.000140310790466315733723475232, -0.000017328176496070286001302184, 1.90540194670935746397168e-6, -1.8882873760163694937908e-7, 1.703176613666840587056e-8, -1.40955218086201517976e-9, 1.0776816914256065828e-10, -7.656138112778696256e-12, 5.07943557413613792e-13, -3.1608615530282912e-14, 1.852036572003432e-15, -1.02524641430496e-16, 5.37852808112e-18, -2.68128238704e-19, 1.273321788e-20, -5.77335744e-22, 2.504352e-23, -1.0446e-24, 4.16e-26, -2.808e-27 }; // inline functions for performance /* ---------------------------------------------------------------------- */ inline double ipoly02(double x) { /* P(x) in the range x > 2 */ int i; double b0, b1, b2; b1 = 0.0; b0 = 0.0; x *= 2; for (i = ERF_TERMS2; i >= 0; i--) { b2 = b1; b1 = b0; b0 = x * b1 - b2 + E2[i]; } return 0.5 * (b0 - b2); } /* ---------------------------------------------------------------------- */ inline double ipoly1(double x) { /* First derivative P'(x) in the range x < 2 */ int i; double b0, b1, b2; b1 = 0.0; b0 = 0.0; x *= 2; for (i = DERF_TERMS; i >= 0; i--) { b2 = b1; b1 = b0; b0 = x * b1 - b2 + DE1[i]; } return 0.5 * (b0 - b2); } /* ---------------------------------------------------------------------- */ inline double ipoly01(double x) { // P(x) in the range x < 2 int i; double b0, b1, b2; b1 = 0.0; b0 = 0.0; x *= 2; for (i = ERF_TERMS1; i >= 0; i--) { b2 = b1; b1 = b0; b0 = x * b1 - b2 + E1[i]; } return 0.5 * (b0 - b2); } /* ---------------------------------------------------------------------- */ inline double ierfoverx1(double x, double *df) { // Computes Erf(x)/x and its first derivative double t, f; double x2; // x squared double exp_term, recip_x; if (x < 2.0) { /* erf(x) = x * y(t) */ /* t = 2 * (x/2)^2 - 1. */ t = 0.5 * x * x - 1; f = ipoly01(t); *df = ipoly1(t) * x; } else { /* erf(x) = 1 - exp(-x^2)/x * y(t) */ /* t = (10.5 - x^2) / (2.5 + x^2) */ x2 = x * x; t = (10.5 - x2) / (2.5 + x2); exp_term = exp(-x2); recip_x = 1.0 / x; f = 1.0 / x - (exp_term / x2) * ipoly02(t); *df = (1.12837916709551257389615890312 * exp_term - f) * recip_x; } return f; } /* ---------------------------------------------------------------------- */ -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; /* sqrt(2) */ coeff_a = 1.4142135623730951; /* E = -Z/r Erf(a r / re) */ /* constants: sqrt(2), 2 / sqrt(pi) */ a = coeff_a / re1; arc = a * rc; /* Interaction between nuclear point charge and Gaussian electron */ 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; } /* ---------------------------------------------------------------------- */ 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; /* sqrt(2) */ coeff_a = 1.4142135623730951; re = sqrt(re1 * re1 + re2 * re2); double ratio; ratio = rc / re; /* constants: sqrt(2), 2 / sqrt(pi) */ a = coeff_a / re; arc = a * rc; /* V_elecelec = E * F */ /* where E = -Z/r Erf(a r / re) */ /* F = (1 - (b s + c s^2) exp(-d s^2)) */ /* and s = r / re */ double E, dEdr, dEdr1, dEdr2, f, df; f = ierfoverx1(arc, &df); dEdr = -a * a * df; fre = a * (f + arc * df) / (re * re); dEdr1 = fre * re1; dEdr2 = fre * re2; E = a * f; - *energy += E; + *ecoul += E; *frc += dEdr; *fre1 += dEdr1; *fre2 += dEdr2; } /* ---------------------------------------------------------------------- */ +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; double dSdr1, dSdr2, dSdr; double dTdr1, dTdr2, dTdr; double O, dOdS, ratio; re1 *= PAULI_RE; re2 *= PAULI_RE; rc *= PAULI_RC; ree = re1 * re1 + re2 * re2; rem = re1 * re1 - re2 * re2; S = (2.82842712474619 / pow((re2 / re1 + re1 / re2), 1.5)) * exp(-rc * rc / ree); t1 = 1.5 * (1 / (re1 * re1) + 1 / (re2 * re2)); t2 = 2.0 * (3 * ree - 2 * rc * rc) / (ree * ree); tt = t1 - t2; dSdr1 = (-1.5 / re1) * (rem / ree) + 2 * re1 * rc * rc / (ree * ree); dSdr2 = (1.5 / re2) * (rem / ree) + 2 * re2 * rc * rc / (ree * ree); dSdr = -2 * rc / ree; dTdr1 = -3 / (re1 * re1 * re1) - 12 * re1 / (ree * ree) + 8 * re1 * (-2 * rc * rc + 3 * ree) / (ree * ree * ree); dTdr2 = -3 / (re2 * re2 * re2) - 12 * re2 / (ree * ree) + 8 * re2 * (-2 * rc * rc + 3 * ree) / (ree * ree * ree); dTdr = 8 * rc / (ree * ree); if (samespin == 1) { O = S * S / (1.0 - S * S) + (1 - PAULI_RHO) * S * S / (1.0 + S * S); dOdS = 2 * S / ((1.0 - S * S) * (1.0 - S * S)) + (1 - PAULI_RHO) * 2 * S / ((1.0 + S * S) * (1.0 + S * S)); } else { O = -PAULI_RHO * S * S / (1.0 + S * S); dOdS = -PAULI_RHO * 2 * S / ((1.0 + S * S) * (1.0 + S * S)); } ratio = tt * dOdS * S; *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; } /* ---------------------------------------------------------------------- */ inline void RForce(double dx, double dy, double dz, double rc, double force, double *fx, double *fy, double *fz) { force /= rc; *fx = force * dx; *fy = force * dy; *fz = force * dz; } /* ---------------------------------------------------------------------- */ inline void SmallRForce(double dx, double dy, double dz, double rc, double force, double *fx, double *fy, double *fz) { /* Handles case where rc is small to avoid division by zero */ if (rc > 0.000001){ force /= rc; *fx = force * dx; *fy = force * dy; *fz = force * dz; } else { if (dx != 0) *fx = force / sqrt(1 + (dy * dy + dz * dz) / (dx * dx)); else *fx = 0.0; if (dy != 0) *fy = force / sqrt(1 + (dx * dx + dz * dz) / (dy * dy)); else *fy = 0.0; if (dz != 0) *fz = force / sqrt(1 + (dx * dx + dy * dy) / (dz * dz)); else *fz = 0.0; // if (dx < 0) *fx = -*fx; // if (dy < 0) *fy = -*fy; // if (dz < 0) *fz = -*fz; } } /* ---------------------------------------------------------------------- */ inline double cutoff(double x) { /* cubic: return x * x * (2.0 * x - 3.0) + 1.0; */ /* quintic: return -6 * pow(x, 5) + 15 * pow(x, 4) - 10 * pow(x, 3) + 1; */ /* Seventh order spline */ // return 20 * pow(x, 7) - 70 * pow(x, 6) + 84 * pow(x, 5) - 35 * pow(x, 4) + 1; return (((20 * x - 70) * x + 84) * x - 35) * x * x * x * x + 1; } /* ---------------------------------------------------------------------- */ inline double dcutoff(double x) { /* cubic: return (6.0 * x * x - 6.0 * x); */ /* quintic: return -30 * pow(x, 4) + 60 * pow(x, 3) - 30 * pow(x, 2); */ /* Seventh order spline */ // return 140 * pow(x, 6) - 420 * pow(x, 5) + 420 * pow(x, 4) - 140 * pow(x, 3); return (((140 * x - 420) * x + 420) * x - 140) * x * x * x; } }