jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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

This commit is contained in:
sjplimp 2010-09-01 18:51:00 +00:00
parent b6ce16eee4
commit 333c561d22
35 changed files with 4471 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/Makefile
View File

@ -16,7 +16,7 @@ OBJ = $(SRC:.cpp=.o)
PACKAGE = asphere class2 colloid dipole dsmc gpu granular \
kspace manybody meam molecule opt peri poems prd reax shock xtc
PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm \
PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm user-eff \
user-ewaldn user-imd user-smd
PACKALL = $(PACKAGE) $(PACKUSER)

73
src/USER-EFF/Install.sh Normal file
View File

@ -0,0 +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_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_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

59
src/USER-EFF/README Normal file
View File

@ -0,0 +1,59 @@
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
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.

798
src/USER-EFF/atom_vec_electron.cpp Normal file
View File

@ -0,0 +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 "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;
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;
}
/* ---------------------------------------------------------------------- */
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<int> (buf[m++]);
type[i] = static_cast<int> (buf[m++]);
mask[i] = static_cast<int> (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<int> (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<int> (buf[m++]);
type[i] = static_cast<int> (buf[m++]);
mask[i] = static_cast<int> (buf[m++]);
q[i] = buf[m++];
spin[i] = static_cast<int> (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<int> (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<int> (buf[m++]);
type[nlocal] = static_cast<int> (buf[m++]);
mask[nlocal] = static_cast<int> (buf[m++]);
image[nlocal] = static_cast<int> (buf[m++]);
q[nlocal] = buf[m++];
spin[nlocal] = static_cast<int> (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<int> (buf[m++]);
type[nlocal] = static_cast<int> (buf[m++]);
mask[nlocal] = static_cast<int> (buf[m++]);
image[nlocal] = static_cast<int> (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<int> (buf[m++]);
eradius[nlocal] = buf[m++];
ervel[nlocal] = buf[m++];
double **extra = atom->extra;
if (atom->nextra_store) {
int size = static_cast<int> (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;
}

74
src/USER-EFF/atom_vec_electron.h Normal file
View File

@ -0,0 +1,74 @@
/* ----------------------------------------------------------------------
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 ATOM_CLASS
AtomStyle(electron,AtomVecElectron)
#else
#ifndef LMP_ATOM_VEC_ELECTRON_H
#define LMP_ATOM_VEC_ELECTRON_H
#include "atom_vec.h"
namespace LAMMPS_NS {
class AtomVecElectron : public AtomVec {
public:
AtomVecElectron(class LAMMPS *, int, char **);
~AtomVecElectron() {}
void grow(int);
void grow_reset();
void copy(int, int);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
int pack_comm_one(int, double *);
void unpack_comm(int, int, double *);
void unpack_comm_vel(int, int, double *);
int unpack_comm_one(int, double *);
int pack_reverse(int, int, double *);
int pack_reverse_one(int, double *);
void unpack_reverse(int, int *, double *);
int unpack_reverse_one(int, double *);
int pack_border(int, int *, double *, int, int *);
int pack_border_vel(int, int *, double *, int, int *);
int pack_border_one(int, double *);
void unpack_border(int, int, double *);
void unpack_border_vel(int, int, double *);
int unpack_border_one(int, double *);
int pack_exchange(int, double *);
int unpack_exchange(double *);
int size_restart();
int pack_restart(int, double *);
int unpack_restart(double *);
void create_atom(int, double *);
void data_atom(double *, int, char **);
int data_atom_hybrid(int, char **);
void data_vel(int, char **);
int data_vel_hybrid(int, char **);
double memory_usage();
private:
double PI;
int *tag,*type,*mask,*image;
double **x,**v,**f;
int *spin;
double *q,*eradius,*ervel,*erforce;
};
}
#endif
#endif

124
src/USER-EFF/compute_ke_atom_eff.cpp Normal file
View File

@ -0,0 +1,124 @@
/* ----------------------------------------------------------------------
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 "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
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])
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;
}

44
src/USER-EFF/compute_ke_atom_eff.h Normal file
View File

@ -0,0 +1,44 @@
/* ----------------------------------------------------------------------
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(ke/atom/eff,ComputeKEAtomEff)
#else
#ifndef LMP_COMPUTE_KE_ATOM_EFF_H
#define LMP_COMPUTE_KE_ATOM_EFF_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeKEAtomEff : public Compute {
public:
ComputeKEAtomEff(class LAMMPS *, int, char **);
~ComputeKEAtomEff();
void init();
void compute_peratom();
double memory_usage();
private:
int nmax;
double *ke;
};
}
#endif
#endif

88
src/USER-EFF/compute_ke_eff.cpp Normal file
View File

@ -0,0 +1,88 @@
/* ----------------------------------------------------------------------
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 "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 {
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];
}
}
MPI_Allreduce(&ke,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
scalar *= pfactor;
return scalar;
}

41
src/USER-EFF/compute_ke_eff.h Normal file
View File

@ -0,0 +1,41 @@
/* ----------------------------------------------------------------------
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(ke/eff,ComputeKEEff)
#else
#ifndef LMP_COMPUTE_KE_EFF_H
#define LMP_COMPUTE_KE_EFF_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeKEEff : public Compute {
public:
ComputeKEEff(class LAMMPS *, int, char **);
~ComputeKEEff() {};
void init();
double compute_scalar();
private:
double pfactor;
};
}
#endif
#endif

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

@ -0,0 +1,180 @@
/* ----------------------------------------------------------------------
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 "compute_temp_deform_eff.h"
#include "update.h"
#include "atom.h"
#include "domain.h"
#include "force.h"
#include "modify.h"
#include "group.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)
{
// error check
if (!atom->spin_flag || !atom->ervel_flag)
error->all("Compute temp/deform/eff requires atom attributes spin, ervel");
}
/* ---------------------------------------------------------------------- */
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++;
}
int nelectrons;
MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world);
// the -3 recovers an extra_dof taken out because it's used by eradius
dof -= domain->dimension * nelectrons - 3;
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];
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];
}
}
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];
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]];
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]) {
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;
}

41
src/USER-EFF/compute_temp_deform_eff.h Normal file
View File

@ -0,0 +1,41 @@
/* ----------------------------------------------------------------------
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"
namespace LAMMPS_NS {
class ComputeTempDeformEff : public ComputeTempDeform {
public:
ComputeTempDeformEff(class LAMMPS *, int, char **);
~ComputeTempDeformEff() {}
double compute_scalar();
void compute_vector();
private:
void dof_compute();
};
}
#endif
#endif

152
src/USER-EFF/compute_temp_eff.cpp Normal file
View File

@ -0,0 +1,152 @@
/* ----------------------------------------------------------------------
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)
{
// error check
if (!atom->spin_flag || !atom->ervel_flag)
error->all("Compute temp/eff requires atom attributes spin, ervel");
}
/* ---------------------------------------------------------------------- */
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++;
}
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 {
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];
}
}
}
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];
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]];
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]) {
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;
}

41
src/USER-EFF/compute_temp_eff.h Normal file
View File

@ -0,0 +1,41 @@
/* ----------------------------------------------------------------------
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"
namespace LAMMPS_NS {
class ComputeTempEff : public ComputeTemp {
public:
ComputeTempEff(class LAMMPS *, int, char **);
~ComputeTempEff() {}
double compute_scalar();
void compute_vector();
private:
void dof_compute();
};
}
#endif
#endif

145
src/USER-EFF/compute_temp_region_eff.cpp Normal file
View File

@ -0,0 +1,145 @@
/* ----------------------------------------------------------------------
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 "domain.h"
#include "region.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeTempRegionEff::
ComputeTempRegionEff(LAMMPS *lmp, int narg, char **arg) :
ComputeTempRegion(lmp, narg, arg)
{
// error check
if (!atom->spin_flag || !atom->ervel_flag)
error->all("Compute temp/region/eff requires atom attributes spin, ervel");
}
/* ---------------------------------------------------------------------- */
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 {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
count++;
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
mass[type[i]];
if (spin[i]) 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++;
}
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]];
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]) {
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;
}

38
src/USER-EFF/compute_temp_region_eff.h Normal file
View File

@ -0,0 +1,38 @@
/* ----------------------------------------------------------------------
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"
namespace LAMMPS_NS {
class ComputeTempRegionEff : public ComputeTempRegion {
public:
ComputeTempRegionEff(class LAMMPS *, int, char **);
~ComputeTempRegionEff() {}
double compute_scalar();
void compute_vector();
};
}
#endif
#endif

338
src/USER-EFF/fix_langevin_eff.cpp Normal file
View File

@ -0,0 +1,338 @@
/* ----------------------------------------------------------------------
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;
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]);
}
}
}
} else {
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);
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]);
}
}
}
}
}
/* ---------------------------------------------------------------------- */
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;
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]);
}
}
}
} 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]);
}
}
}
}
}
/* ----------------------------------------------------------------------
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;
energy_onestep = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
flangevin[i][2]*v[i][2] + 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;
if (update->ntimestep == update->beginstep) {
energy_onestep = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
flangevin[i][2]*v[i][2] + 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;
}

44
src/USER-EFF/fix_langevin_eff.h Normal file
View File

@ -0,0 +1,44 @@
/* ----------------------------------------------------------------------
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(langevin/eff,FixLangevinEff)
#else
#ifndef LMP_FIX_LANGEVIN_EFF_H
#define LMP_FIX_LANGEVIN_EFF_H
#include "fix_langevin.h"
namespace LAMMPS_NS {
class FixLangevinEff : public FixLangevin {
public:
FixLangevinEff(class LAMMPS *, int, char **);
~FixLangevinEff();
void end_of_step();
double compute_scalar();
private:
double *erforcelangevin;
void post_force_no_tally();
void post_force_tally();
};
}
#endif
#endif

126
src/USER-EFF/fix_nh_eff.cpp Normal file
View File

@ -0,0 +1,126 @@
/* ----------------------------------------------------------------------
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 "
"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;
}
}
}
}
}
/* ----------------------------------------------------------------------
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];
}
/* ----------------------------------------------------------------------
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;
}

34
src/USER-EFF/fix_nh_eff.h Normal file
View File

@ -0,0 +1,34 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifndef LMP_FIX_NH_EFF_H
#define LMP_FIX_NH_EFF_H
#include "fix_nh.h"
namespace LAMMPS_NS {
class FixNHEff : public FixNH {
public:
FixNHEff(class LAMMPS *, int, char **);
virtual ~FixNHEff() {}
protected:
void nve_v();
void nve_x();
virtual void nh_v_temp();
};
}
#endif

67
src/USER-EFF/fix_nph_eff.cpp Normal file
View File

@ -0,0 +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 "fix_nph_eff.h"
#include "modify.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
FixNPHEff::FixNPHEff(LAMMPS *lmp, int narg, char **arg) :
FixNHEff(lmp, narg, arg)
{
if (tstat_flag)
error->all("Temperature control can not be used with fix nph/eff");
if (!pstat_flag)
error->all("Pressure control must be used with fix nph/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;
}

36
src/USER-EFF/fix_nph_eff.h Normal file
View File

@ -0,0 +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)
#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

78
src/USER-EFF/fix_npt_eff.cpp Normal file
View File

@ -0,0 +1,78 @@
/* ----------------------------------------------------------------------
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;
}

36
src/USER-EFF/fix_npt_eff.h Normal file
View File

@ -0,0 +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(npt/eff,FixNPTEff)
#else
#ifndef LMP_FIX_NPT_EFF_H
#define LMP_FIX_NPT_EFF_H
#include "fix_nh_eff.h"
namespace LAMMPS_NS {
class FixNPTEff : public FixNHEff {
public:
FixNPTEff(class LAMMPS *, int, char **);
~FixNPTEff() {}
};
}
#endif
#endif

149
src/USER-EFF/fix_nve_eff.cpp Normal file
View File

@ -0,0 +1,149 @@
/* ----------------------------------------------------------------------
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)
{
// 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");
}
/* ----------------------------------------------------------------------
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) {
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];
}
}
}
}
}
/* ---------------------------------------------------------------------- */
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 {
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;
}
}
}
}

37
src/USER-EFF/fix_nve_eff.h Normal file
View File

@ -0,0 +1,37 @@
/* ----------------------------------------------------------------------
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"
namespace LAMMPS_NS {
class FixNVEEff : public FixNVE {
public:
FixNVEEff(class LAMMPS *, int, char **);
void initial_integrate(int);
void final_integrate();
};
}
#endif
#endif

49
src/USER-EFF/fix_nvt_eff.cpp Normal file
View File

@ -0,0 +1,49 @@
/* ----------------------------------------------------------------------
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 "math.h"
#include "fix_nvt_eff.h"
#include "group.h"
#include "modify.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
FixNVTEff::FixNVTEff(LAMMPS *lmp, int narg, char **arg) :
FixNHEff(lmp, narg, arg)
{
if (!tstat_flag)
error->all("Temperature control must be used with fix nvt/eff");
if (pstat_flag)
error->all("Pressure control can not be used with fix nvt/eff");
// 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/eff";
modify->add_compute(3,newarg);
delete [] newarg;
tflag = 1;
}

36
src/USER-EFF/fix_nvt_eff.h Normal file
View File

@ -0,0 +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() {}
};
}
#endif
#endif

127
src/USER-EFF/fix_nvt_sllod_eff.cpp Normal file
View File

@ -0,0 +1,127 @@
/* ----------------------------------------------------------------------
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)
{
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");
// 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();
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])
ervel[i] = ervel[i]*factor_eta -
dthalf*sqrt(pow(vdelu[0],2)+pow(vdelu[1],2)+pow(vdelu[2],2));
}
}
}

42
src/USER-EFF/fix_nvt_sllod_eff.h Normal file
View File

@ -0,0 +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"
namespace LAMMPS_NS {
class FixNVTSllodEff : public FixNVTEff {
public:
FixNVTSllodEff(class LAMMPS *, int, char **);
~FixNVTSllodEff() {}
void init();
private:
int nondeformbias;
void nh_v_temp();
};
}
#endif
#endif

109
src/USER-EFF/fix_temp_rescale_eff.cpp Normal file
View File

@ -0,0 +1,109 @@
/* ----------------------------------------------------------------------
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 "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)
{
// create a new compute temp/eff, wiping out one parent class just created
// id = fix-ID + temp, compute group = fix group
modify->delete_compute(id_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;
}
/* ---------------------------------------------------------------------- */
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])
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])
ervel[i] *= factor;
temperature->restore_bias(i,v[i]);
}
}
}
}
}

37
src/USER-EFF/fix_temp_rescale_eff.h Normal file
View File

@ -0,0 +1,37 @@
/* ----------------------------------------------------------------------
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"
namespace LAMMPS_NS {
class FixTempRescaleEff : public FixTempRescale {
public:
FixTempRescaleEff(class LAMMPS *, int, char **);
~FixTempRescaleEff() {}
void end_of_step();
};
}
#endif
#endif

715
src/USER-EFF/pair_eff_cut.cpp Normal file
View File

@ -0,0 +1,715 @@
/* ----------------------------------------------------------------------
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)
------------------------------------------------------------------------- */
#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;
}
/* ---------------------------------------------------------------------- */
PairEffCut::~PairEffCut()
{
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;
int *ilist,*jlist,*numneigh,**firstneigh;
ecoul = 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
if (spin[i] != 0) {
e1rforce = energy = 0;
energy = 1.5 / (eradius[i] * eradius[i]);
e1rforce = 3.0 / (eradius[i] * eradius[i] * eradius[i]);
erforce[i] += e1rforce;
// electronic ke accumulates into ecoul (pot)
if (eflag) ecoul = energy; // KE e-wavefunction
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,e1rforce*eradius[i]);
}
}
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
if (spin[i] == 0 && spin[j] == 0) {
energy = fx = fy = fz = 0;
double qxq = qqrd2e*qtmp*q[j];
forcecoul = qxq/rsq;
double dist = rc / taper;
double spline = cutoff(dist);
double dspline = dcutoff(dist) / taper;
energy = factor_coul*qxq/rc;
fpair = forcecoul*spline-energy*dspline;
fpair = qqrd2e*fpair/rc;
energy = spline*energy;
fx = delx*fpair;
fy = dely*fpair;
fz = delz*fpair;
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
if (newton_pair || j < nlocal) {
f[j][0] -= fx;
f[j][1] -= fy;
f[j][2] -= fz;
}
if (eflag) ecoul = energy; // Electrostatics:N-N
if (evflag)
ev_tally_xyz(i,j,nlocal,newton_pair,0.0,ecoul,
fx,fy,fz,delx,dely,delz);
}
// I is nucleus, J is electron
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;
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);
}
}
// I is electon, J is nucleus
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;
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);
}
}
// 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)
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
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
SmallRForce(delx, dely, delz, rc, fpair + s_fpair, &fx, &fy, &fz);
erforce[i] += e1rforce;
erforce[j] += e2rforce;
// radial virials
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);
}
}
}
}
// limit the electron size for periodic systems, to max=half-box-size
// limit_size_stiffness for electrons
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;
// constraint radial energy accumulated as ecoul
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);
}
}
}
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 ecoulhalf,epairhalf;
double partial_evirial = e_virial/3.0;
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 (eflag_atom) {
epairhalf = 0.5 * ecoul;
if (i < nlocal) eatom[i] += epairhalf;
if (j < nlocal) eatom[j] += epairhalf;
}
}
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;
}
if (spin[j] && j < nlocal) {
virial[0] += 0.5*partial_evirial;
virial[1] += 0.5*partial_evirial;
virial[2] += 0.5*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;
}
}
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;
}
}
}
}
}
/* ----------------------------------------------------------------------
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) {
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]);
}
// 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)
int ignore = update->minimize->request(this,1,0.01);
// 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;
}

63
src/USER-EFF/pair_eff_cut.h Normal file
View File

@ -0,0 +1,63 @@
/* ----------------------------------------------------------------------
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;
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

446
src/USER-EFF/pair_eff_inline.h Normal file
View File

@ -0,0 +1,446 @@
/* ----------------------------------------------------------------------
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)
------------------------------------------------------------------------- */
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)
{
*energy += 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 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;
*energy += 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 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;
*frc += dEdr;
*fre1 += dEdr1;
*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 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;
}
/* ---------------------------------------------------------------------- */
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;
}
}

6
src/fix_temp_rescale.h
View File

@ -27,15 +27,15 @@ namespace LAMMPS_NS {
class FixTempRescale : public Fix {
public:
FixTempRescale(class LAMMPS *, int, char **);
~FixTempRescale();
virtual ~FixTempRescale();
int setmask();
void init();
void end_of_step();
virtual void end_of_step();
int modify_param(int, char **);
void reset_target(double);
double compute_scalar();
private:
protected:
int which;
double t_start,t_stop,t_window;
double fraction,energy,efactor;