diff --git a/src/Makefile b/src/Makefile index aa39dff1c..7b32442c3 100755 --- a/src/Makefile +++ b/src/Makefile @@ -1,185 +1,185 @@ # LAMMPS multiple-machine Makefile SHELL = /bin/sh #.IGNORE: # Definitions ROOT = lmp EXE = $(ROOT)_$@ SRC = $(wildcard *.cpp) INC = $(wildcard *.h) OBJ = $(SRC:.cpp=.o) # Package variables 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) PACKAGEUC = $(shell perl -e 'printf("%s", uc("$(PACKAGE)"));') PACKUSERUC = $(shell perl -e 'printf("%s", uc("$(PACKUSER)"));') YESDIR = $(shell perl -e 'printf("%s", uc("$(@:yes-%=%)"));') NODIR = $(shell perl -e 'printf("%s", uc("$(@:no-%=%)"));') # List of all targets help: @echo '' @echo 'make clean-all delete all object files' @echo 'make clean-machine delete object files for one machine' @echo 'make tar lmp_src.tar.gz of src dir and packages' @echo 'make makelib update Makefile.lib for library build' @echo 'make makelist update Makefile.list used by old makes' @echo '' @echo 'make package list available packages' @echo 'make package-status status of all packages' @echo 'make yes-package install a single package in src dir' @echo 'make no-package remove a single package from src dir' @echo 'make yes-all install all packages in src dir' @echo 'make no-all remove all packages from src dir' @echo 'make yes-standard install all standard packages' @echo 'make no-standard remove all standard packages' @echo 'make yes-user install all user packages' @echo 'make no-user remove all user packages' @echo '' @echo 'make package-update replace src files with package files' @echo 'make package-overwrite replace package files with src files' @echo '' @echo 'make machine build LAMMPS where machine is one of:' @echo '' @files="`ls MAKE/Makefile.*`"; \ for file in $$files; do head -1 $$file; done @echo '' # Build the code .DEFAULT: @test -f MAKE/Makefile.$@ @if [ ! -d Obj_$@ ]; then mkdir Obj_$@; fi @/bin/sh Make.sh style @cp -p *.cpp *.h Obj_$@ @cp MAKE/Makefile.$@ Obj_$@/Makefile @if [ ! -e Makefile.package ]; then make package-regenerate; fi @cp Makefile.package Obj_$@ @cd Obj_$@; \ $(MAKE) $(MFLAGS) "OBJ = $(OBJ)" "INC = $(INC)" "EXE = ../$(EXE)" ../$(EXE) @if [ -d Obj_$@ ]; then cd Obj_$@; rm -f $(SRC) $(INC) Makefile*; fi # Remove machine-specific object files clean: @echo 'make clean-all delete all object files' @echo 'make clean-machine delete object files for one machine' clean-all: rm -rf Obj_* clean-%: rm -rf Obj_$(@:clean-%=%) # Create a tarball of src dir and packages tar: @cd STUBS; make clean @cd ..; tar cvzf src/$(ROOT)_src.tar.gz \ src/Make* src/Package.sh src/MAKE src/*.cpp src/*.h src/STUBS \ $(patsubst %,src/%,$(PACKAGEUC)) $(patsubst %,src/%,$(PACKUSERUC)) \ --exclude=*/.svn @cd STUBS; make @echo "Created $(ROOT)_src.tar.gz" # Update Makefile.lib and Makefile.list makelib: @/bin/sh Make.sh style @/bin/sh Make.sh Makefile.lib makelist: @/bin/sh Make.sh style @/bin/sh Make.sh Makefile.list # Package management package: @echo 'Standard packages:' $(PACKAGE) @echo '' @echo 'User-contributed packages:' $(PACKUSER) @echo '' @echo 'make package list available packages' @echo 'make package-status status of all packages' @echo 'make yes-package install a single package in src dir' @echo 'make no-package remove a single package from src dir' @echo 'make yes-all install all packages in src dir' @echo 'make no-all remove all packages from src dir' @echo 'make yes-standard install all standard packages' @echo 'make no-standard remove all standard packages' @echo 'make yes-user install all user packages' @echo 'make no-user remove all user packages' @echo '' @echo 'make package-update replace src files with package files' @echo 'make package-overwrite replace package files with src files' yes-all: @for p in $(PACKALL); do $(MAKE) yes-$$p; done no-all: @for p in $(PACKALL); do $(MAKE) no-$$p; done yes-standard: @for p in $(PACKAGE); do $(MAKE) yes-$$p; done no-standard: @for p in $(PACKAGE); do $(MAKE) no-$$p; done yes-user: @for p in $(PACKUSER); do $(MAKE) yes-$$p; done no-user: @for p in $(PACKUSER); do $(MAKE) no-$$p; done yes-%: @if [ ! -e $(YESDIR) ]; then \ echo "Package $(@:yes-%=%) does not exist"; \ else \ echo "Installing package $(@:yes-%=%)"; \ cd $(YESDIR); /bin/sh Install.sh 1; \ fi; no-%: @if [ ! -e $(NODIR) ]; then \ echo "Package $(@:no-%=%) does not exist"; \ else \ echo "Uninstalling package $(@:no-%=%), ignore errors"; \ cd $(NODIR); /bin/sh Install.sh 0; cd ..; \ fi; # status = list differences between src and package files # update = replace src files with newer package files # overwrite = overwrite package files with newer src files # regenerate = regenerate Makefile.package from Makefile.package.empty package-status: @for p in $(PACKAGEUC); do /bin/sh Package.sh $$p status; done @echo '' @for p in $(PACKUSERUC); do /bin/sh Package.sh $$p status; done package-update: @for p in $(PACKAGEUC); do /bin/sh Package.sh $$p update; done @echo '' @for p in $(PACKUSERUC); do /bin/sh Package.sh $$p update; done package-overwrite: @for p in $(PACKAGEUC); do /bin/sh Package.sh $$p overwrite; done @echo '' @for p in $(PACKUSERUC); do /bin/sh Package.sh $$p overwrite; done package-regenerate: @cp Makefile.package.empty Makefile.package @for p in $(PACKAGEUC); do /bin/sh Package.sh $$p regenerate; done @for p in $(PACKUSERUC); do /bin/sh Package.sh $$p regenerate; done diff --git a/src/USER-EFF/Install.sh b/src/USER-EFF/Install.sh new file mode 100644 index 000000000..407168c56 --- /dev/null +++ b/src/USER-EFF/Install.sh @@ -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 diff --git a/src/USER-EFF/README b/src/USER-EFF/README new file mode 100644 index 000000000..2e500bafc --- /dev/null +++ b/src/USER-EFF/README @@ -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. diff --git a/src/USER-EFF/atom_vec_electron.cpp b/src/USER-EFF/atom_vec_electron.cpp new file mode 100644 index 000000000..f3a399768 --- /dev/null +++ b/src/USER-EFF/atom_vec_electron.cpp @@ -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 (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + q[i] = buf[m++]; + spin[i] = static_cast (buf[m++]); + eradius[i] = buf[m++]; + ervel[i] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecElectron::pack_border(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + double dx,dy,dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = q[j]; + buf[m++] = spin[j]; + buf[m++] = eradius[j]; + buf[m++] = ervel[j]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0]*domain->xprd; + dy = pbc[1]*domain->yprd; + dz = pbc[2]*domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = q[j]; + buf[m++] = spin[j]; + buf[m++] = eradius[j]; + buf[m++] = ervel[j]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecElectron::pack_border_one(int i, double *buf) +{ + buf[0] = q[i]; + buf[1] = spin[i]; + buf[2] = eradius[i]; + buf[3] = ervel[i]; + + return 4; +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecElectron::unpack_border(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = static_cast (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + q[i] = buf[m++]; + spin[i] = static_cast (buf[m++]); + eradius[i] = buf[m++]; + ervel[i] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecElectron::unpack_border_one(int i, double *buf) +{ + q[i] = buf[0]; + spin[i] = static_cast (buf[1]); + eradius[i] = buf[2]; + ervel[i] = buf[3]; + + return 4; +} + +/* ---------------------------------------------------------------------- + pack data for atom I for sending to another proc + xyz must be 1st 3 values, so comm::exchange() can test on them +------------------------------------------------------------------------- */ + +int AtomVecElectron::pack_exchange(int i, double *buf) +{ + int m = 1; + buf[m++] = x[i][0]; + buf[m++] = x[i][1]; + buf[m++] = x[i][2]; + buf[m++] = v[i][0]; + buf[m++] = v[i][1]; + buf[m++] = v[i][2]; + buf[m++] = tag[i]; + buf[m++] = type[i]; + buf[m++] = mask[i]; + buf[m++] = image[i]; + buf[m++] = q[i]; + buf[m++] = spin[i]; + buf[m++] = eradius[i]; + buf[m++] = ervel[i]; + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]); + + buf[0] = m; + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecElectron::unpack_exchange(double *buf) +{ + int nlocal = atom->nlocal; + if (nlocal == nmax) grow(0); + + int m = 1; + x[nlocal][0] = buf[m++]; + x[nlocal][1] = buf[m++]; + x[nlocal][2] = buf[m++]; + v[nlocal][0] = buf[m++]; + v[nlocal][1] = buf[m++]; + v[nlocal][2] = buf[m++]; + tag[nlocal] = static_cast (buf[m++]); + type[nlocal] = static_cast (buf[m++]); + mask[nlocal] = static_cast (buf[m++]); + image[nlocal] = static_cast (buf[m++]); + q[nlocal] = buf[m++]; + spin[nlocal] = static_cast (buf[m++]); + eradius[nlocal] = buf[m++]; + ervel[nlocal] = buf[m++]; + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + m += modify->fix[atom->extra_grow[iextra]]-> + unpack_exchange(nlocal,&buf[m]); + + atom->nlocal++; + return m; +} + +/* ---------------------------------------------------------------------- + size of restart data for all atoms owned by this proc + include extra data stored by fixes +------------------------------------------------------------------------- */ + +int AtomVecElectron::size_restart() +{ + int i; + + int nlocal = atom->nlocal; + int n = 15 * nlocal; // Associated with pack_restart + + if (atom->nextra_restart) + for (int iextra = 0; iextra < atom->nextra_restart; iextra++) + for (i = 0; i < nlocal; i++) + n += modify->fix[atom->extra_restart[iextra]]->size_restart(i); + + return n; +} + +/* ---------------------------------------------------------------------- + pack atom I's data for restart file including extra quantities + xyz must be 1st 3 values, so that read_restart can test on them + molecular types may be negative, but write as positive +------------------------------------------------------------------------- */ + +int AtomVecElectron::pack_restart(int i, double *buf) +{ + int m = 1; + buf[m++] = x[i][0]; + buf[m++] = x[i][1]; + buf[m++] = x[i][2]; + buf[m++] = tag[i]; + buf[m++] = type[i]; + buf[m++] = mask[i]; + buf[m++] = image[i]; + buf[m++] = v[i][0]; + buf[m++] = v[i][1]; + buf[m++] = v[i][2]; + + buf[m++] = q[i]; + buf[m++] = spin[i]; + buf[m++] = eradius[i]; + buf[m++] = ervel[i]; + + if (atom->nextra_restart) + for (int iextra = 0; iextra < atom->nextra_restart; iextra++) + m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]); + + buf[0] = m; + return m; +} + +/* ---------------------------------------------------------------------- + unpack data for one atom from restart file including extra quantities +------------------------------------------------------------------------- */ + +int AtomVecElectron::unpack_restart(double *buf) +{ + int nlocal = atom->nlocal; + if (nlocal == nmax) { + grow(0); + if (atom->nextra_store) + atom->extra = memory->grow_2d_double_array(atom->extra,nmax, + atom->nextra_store, + "atom:extra"); + } + + int m = 1; + x[nlocal][0] = buf[m++]; + x[nlocal][1] = buf[m++]; + x[nlocal][2] = buf[m++]; + tag[nlocal] = static_cast (buf[m++]); + type[nlocal] = static_cast (buf[m++]); + mask[nlocal] = static_cast (buf[m++]); + image[nlocal] = static_cast (buf[m++]); + v[nlocal][0] = buf[m++]; + v[nlocal][1] = buf[m++]; + v[nlocal][2] = buf[m++]; + + q[nlocal] = buf[m++]; + spin[nlocal] = static_cast (buf[m++]); + eradius[nlocal] = buf[m++]; + ervel[nlocal] = buf[m++]; + + double **extra = atom->extra; + if (atom->nextra_store) { + int size = static_cast (buf[0]) - m; + for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++]; + } + + atom->nlocal++; + return m; +} + +/* ---------------------------------------------------------------------- + create one atom of itype at coord + set other values to defaults +------------------------------------------------------------------------- */ + +void AtomVecElectron::create_atom(int itype, double *coord) +{ + int nlocal = atom->nlocal; + if (nlocal == nmax) grow(0); + + tag[nlocal] = 0; + type[nlocal] = itype; + x[nlocal][0] = coord[0]; + x[nlocal][1] = coord[1]; + x[nlocal][2] = coord[2]; + mask[nlocal] = 1; + image[nlocal] = (512 << 20) | (512 << 10) | 512; + v[nlocal][0] = 0.0; + v[nlocal][1] = 0.0; + v[nlocal][2] = 0.0; + + q[nlocal] = 0.0; + spin[nlocal] = 1; + eradius[nlocal] = 1.0; + ervel[nlocal] = 0.0; + + atom->nlocal++; +} + +/* ---------------------------------------------------------------------- + unpack one line from Atoms section of data file + initialize other atom quantities +------------------------------------------------------------------------- */ + +void AtomVecElectron::data_atom(double *coord, int imagetmp, char **values) +{ + int nlocal = atom->nlocal; + + if (nlocal == nmax) grow(0); + + tag[nlocal] = atoi(values[0]); + if (tag[nlocal] <= 0) + error->one("Invalid atom ID in Atoms section of data file"); + + type[nlocal] = atoi(values[1]); + if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes) + error->one("Invalid atom type in Atoms section of data file"); + + q[nlocal] = atof(values[2]); + spin[nlocal] = atoi(values[3]); + eradius[nlocal] = atof(values[4]); + + x[nlocal][0] = coord[0]; + x[nlocal][1] = coord[1]; + x[nlocal][2] = coord[2]; + + image[nlocal] = imagetmp; + + mask[nlocal] = 1; + v[nlocal][0] = 0.0; + v[nlocal][1] = 0.0; + v[nlocal][2] = 0.0; + ervel[nlocal] = 0.0; + + atom->nlocal++; +} + +/* ---------------------------------------------------------------------- + unpack hybrid quantities from one line in Atoms section of data file + initialize other atom quantities for this sub-style +------------------------------------------------------------------------- */ + +int AtomVecElectron::data_atom_hybrid(int nlocal, char **values) +{ + q[nlocal] = atof(values[0]); + spin[nlocal] = atoi(values[1]); + eradius[nlocal] = atof(values[2]); + if (eradius[nlocal] < 0.0) + error->one("Invalid eradius in Atoms section of data file"); + + v[nlocal][0] = 0.0; + v[nlocal][1] = 0.0; + v[nlocal][2] = 0.0; + ervel[nlocal] = 0.0; + + return 2; +} + +/* ---------------------------------------------------------------------- + unpack one line from Velocities section of data file +------------------------------------------------------------------------- */ + +void AtomVecElectron::data_vel(int m, char **values) +{ + v[m][0] = atof(values[0]); + v[m][1] = atof(values[1]); + v[m][2] = atof(values[2]); + ervel[m] = atof(values[3]); +} + +/* ---------------------------------------------------------------------- + unpack hybrid quantities from one line in Velocities section of data file +------------------------------------------------------------------------- */ + +int AtomVecElectron::data_vel_hybrid(int m, char **values) +{ + ervel[m] = atof(values[0]); + + return 1; +} + +/* ---------------------------------------------------------------------- + return # of bytes of allocated memory +------------------------------------------------------------------------- */ + +double AtomVecElectron::memory_usage() +{ + double bytes = 0.0; + + if (atom->memcheck("tag")) bytes += nmax * sizeof(int); + if (atom->memcheck("type")) bytes += nmax * sizeof(int); + if (atom->memcheck("mask")) bytes += nmax * sizeof(int); + if (atom->memcheck("image")) bytes += nmax * sizeof(int); + if (atom->memcheck("x")) bytes += nmax*3 * sizeof(double); + if (atom->memcheck("v")) bytes += nmax*3 * sizeof(double); + if (atom->memcheck("f")) bytes += nmax*3 * sizeof(double); + + if (atom->memcheck("q")) bytes += nmax * sizeof(double); + if (atom->memcheck("spin")) bytes += nmax * sizeof(int); + if (atom->memcheck("eradius")) bytes += nmax * sizeof(double); + if (atom->memcheck("ervel")) bytes += nmax * sizeof(double); + if (atom->memcheck("erforce")) bytes += nmax * sizeof(double); + + return bytes; +} diff --git a/src/USER-EFF/atom_vec_electron.h b/src/USER-EFF/atom_vec_electron.h new file mode 100644 index 000000000..704a8378f --- /dev/null +++ b/src/USER-EFF/atom_vec_electron.h @@ -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 + diff --git a/src/USER-EFF/compute_ke_atom_eff.cpp b/src/USER-EFF/compute_ke_atom_eff.cpp new file mode 100644 index 000000000..f9dba6e5d --- /dev/null +++ b/src/USER-EFF/compute_ke_atom_eff.cpp @@ -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; +} diff --git a/src/fix_temp_rescale.h b/src/USER-EFF/compute_ke_atom_eff.h similarity index 58% copy from src/fix_temp_rescale.h copy to src/USER-EFF/compute_ke_atom_eff.h index f710bc7cd..1f04412e5 100644 --- a/src/fix_temp_rescale.h +++ b/src/USER-EFF/compute_ke_atom_eff.h @@ -1,51 +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(temp/rescale,FixTempRescale) +#ifdef COMPUTE_CLASS + +ComputeStyle(ke/atom/eff,ComputeKEAtomEff) #else -#ifndef LMP_FIX_TEMP_RESCALE_H -#define LMP_FIX_TEMP_RESCALE_H +#ifndef LMP_COMPUTE_KE_ATOM_EFF_H +#define LMP_COMPUTE_KE_ATOM_EFF_H -#include "fix.h" +#include "compute.h" namespace LAMMPS_NS { -class FixTempRescale : public Fix { +class ComputeKEAtomEff : public Compute { public: - FixTempRescale(class LAMMPS *, int, char **); - ~FixTempRescale(); - int setmask(); + ComputeKEAtomEff(class LAMMPS *, int, char **); + ~ComputeKEAtomEff(); void init(); - void end_of_step(); - int modify_param(int, char **); - void reset_target(double); - double compute_scalar(); - + void compute_peratom(); + double memory_usage(); + private: - int which; - double t_start,t_stop,t_window; - double fraction,energy,efactor; - - char *id_temp; - class Compute *temperature; - int tflag; + int nmax; + double *ke; }; } #endif #endif diff --git a/src/USER-EFF/compute_ke_eff.cpp b/src/USER-EFF/compute_ke_eff.cpp new file mode 100644 index 000000000..600f0ebae --- /dev/null +++ b/src/USER-EFF/compute_ke_eff.cpp @@ -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; +} diff --git a/src/fix_temp_rescale.h b/src/USER-EFF/compute_ke_eff.h similarity index 60% copy from src/fix_temp_rescale.h copy to src/USER-EFF/compute_ke_eff.h index f710bc7cd..e8a31db55 100644 --- a/src/fix_temp_rescale.h +++ b/src/USER-EFF/compute_ke_eff.h @@ -1,51 +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 FIX_CLASS +#ifdef COMPUTE_CLASS -FixStyle(temp/rescale,FixTempRescale) +ComputeStyle(ke/eff,ComputeKEEff) #else -#ifndef LMP_FIX_TEMP_RESCALE_H -#define LMP_FIX_TEMP_RESCALE_H +#ifndef LMP_COMPUTE_KE_EFF_H +#define LMP_COMPUTE_KE_EFF_H -#include "fix.h" +#include "compute.h" namespace LAMMPS_NS { -class FixTempRescale : public Fix { +class ComputeKEEff : public Compute { public: - FixTempRescale(class LAMMPS *, int, char **); - ~FixTempRescale(); - int setmask(); + ComputeKEEff(class LAMMPS *, int, char **); + ~ComputeKEEff() {}; void init(); - void end_of_step(); - int modify_param(int, char **); - void reset_target(double); double compute_scalar(); - + private: - int which; - double t_start,t_stop,t_window; - double fraction,energy,efactor; - - char *id_temp; - class Compute *temperature; - int tflag; + double pfactor; }; - + } #endif #endif diff --git a/src/USER-EFF/compute_temp_deform_eff.cpp b/src/USER-EFF/compute_temp_deform_eff.cpp new file mode 100644 index 000000000..e05e4d131 --- /dev/null +++ b/src/USER-EFF/compute_temp_deform_eff.cpp @@ -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; +} diff --git a/src/fix_temp_rescale.h b/src/USER-EFF/compute_temp_deform_eff.h similarity index 59% copy from src/fix_temp_rescale.h copy to src/USER-EFF/compute_temp_deform_eff.h index f710bc7cd..32bf165c4 100644 --- a/src/fix_temp_rescale.h +++ b/src/USER-EFF/compute_temp_deform_eff.h @@ -1,51 +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 FIX_CLASS +#ifdef COMPUTE_CLASS -FixStyle(temp/rescale,FixTempRescale) +ComputeStyle(temp/deform/eff,ComputeTempDeformEff) #else -#ifndef LMP_FIX_TEMP_RESCALE_H -#define LMP_FIX_TEMP_RESCALE_H +#ifndef LMP_COMPUTE_TEMP_DEFORM_EFF_H +#define LMP_COMPUTE_TEMP_DEFORM_EFF_H -#include "fix.h" +#include "compute_temp_deform.h" namespace LAMMPS_NS { -class FixTempRescale : public Fix { +class ComputeTempDeformEff : public ComputeTempDeform { public: - FixTempRescale(class LAMMPS *, int, char **); - ~FixTempRescale(); - int setmask(); - void init(); - void end_of_step(); - int modify_param(int, char **); - void reset_target(double); + ComputeTempDeformEff(class LAMMPS *, int, char **); + ~ComputeTempDeformEff() {} double compute_scalar(); + void compute_vector(); private: - int which; - double t_start,t_stop,t_window; - double fraction,energy,efactor; - - char *id_temp; - class Compute *temperature; - int tflag; + void dof_compute(); }; } #endif #endif diff --git a/src/USER-EFF/compute_temp_eff.cpp b/src/USER-EFF/compute_temp_eff.cpp new file mode 100644 index 000000000..26f6539e6 --- /dev/null +++ b/src/USER-EFF/compute_temp_eff.cpp @@ -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; +} diff --git a/src/fix_temp_rescale.h b/src/USER-EFF/compute_temp_eff.h similarity index 59% copy from src/fix_temp_rescale.h copy to src/USER-EFF/compute_temp_eff.h index f710bc7cd..683b73a29 100644 --- a/src/fix_temp_rescale.h +++ b/src/USER-EFF/compute_temp_eff.h @@ -1,51 +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 FIX_CLASS +#ifdef COMPUTE_CLASS -FixStyle(temp/rescale,FixTempRescale) +ComputeStyle(temp/eff,ComputeTempEff) #else -#ifndef LMP_FIX_TEMP_RESCALE_H -#define LMP_FIX_TEMP_RESCALE_H +#ifndef LMP_COMPUTE_TEMP_EFF_H +#define LMP_COMPUTE_TEMP_EFF_H -#include "fix.h" +#include "compute_temp.h" namespace LAMMPS_NS { - -class FixTempRescale : public Fix { + +class ComputeTempEff : public ComputeTemp { public: - FixTempRescale(class LAMMPS *, int, char **); - ~FixTempRescale(); - int setmask(); - void init(); - void end_of_step(); - int modify_param(int, char **); - void reset_target(double); + ComputeTempEff(class LAMMPS *, int, char **); + ~ComputeTempEff() {} double compute_scalar(); - + void compute_vector(); + private: - int which; - double t_start,t_stop,t_window; - double fraction,energy,efactor; - - char *id_temp; - class Compute *temperature; - int tflag; + void dof_compute(); }; - + } #endif #endif diff --git a/src/USER-EFF/compute_temp_region_eff.cpp b/src/USER-EFF/compute_temp_region_eff.cpp new file mode 100644 index 000000000..be6bee419 --- /dev/null +++ b/src/USER-EFF/compute_temp_region_eff.cpp @@ -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; +} diff --git a/src/fix_temp_rescale.h b/src/USER-EFF/compute_temp_region_eff.h similarity index 58% copy from src/fix_temp_rescale.h copy to src/USER-EFF/compute_temp_region_eff.h index f710bc7cd..e2d7add69 100644 --- a/src/fix_temp_rescale.h +++ b/src/USER-EFF/compute_temp_region_eff.h @@ -1,51 +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 FIX_CLASS +#ifdef COMPUTE_CLASS -FixStyle(temp/rescale,FixTempRescale) +ComputeStyle(temp/region/eff,ComputeTempRegionEff) #else -#ifndef LMP_FIX_TEMP_RESCALE_H -#define LMP_FIX_TEMP_RESCALE_H +#ifndef LMP_COMPUTE_TEMP_REGION_EFF_H +#define LMP_COMPUTE_TEMP_REGION_EFF_H -#include "fix.h" +#include "compute_temp_region.h" namespace LAMMPS_NS { -class FixTempRescale : public Fix { +class ComputeTempRegionEff : public ComputeTempRegion { public: - FixTempRescale(class LAMMPS *, int, char **); - ~FixTempRescale(); - int setmask(); - void init(); - void end_of_step(); - int modify_param(int, char **); - void reset_target(double); + ComputeTempRegionEff(class LAMMPS *, int, char **); + ~ComputeTempRegionEff() {} double compute_scalar(); - - private: - int which; - double t_start,t_stop,t_window; - double fraction,energy,efactor; - - char *id_temp; - class Compute *temperature; - int tflag; + void compute_vector(); }; } #endif #endif diff --git a/src/USER-EFF/fix_langevin_eff.cpp b/src/USER-EFF/fix_langevin_eff.cpp new file mode 100644 index 000000000..7aae42b52 --- /dev/null +++ b/src/USER-EFF/fix_langevin_eff.cpp @@ -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; +} diff --git a/src/fix_temp_rescale.h b/src/USER-EFF/fix_langevin_eff.h similarity index 62% copy from src/fix_temp_rescale.h copy to src/USER-EFF/fix_langevin_eff.h index f710bc7cd..6eb78ff9b 100644 --- a/src/fix_temp_rescale.h +++ b/src/USER-EFF/fix_langevin_eff.h @@ -1,51 +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(temp/rescale,FixTempRescale) +FixStyle(langevin/eff,FixLangevinEff) #else -#ifndef LMP_FIX_TEMP_RESCALE_H -#define LMP_FIX_TEMP_RESCALE_H +#ifndef LMP_FIX_LANGEVIN_EFF_H +#define LMP_FIX_LANGEVIN_EFF_H -#include "fix.h" +#include "fix_langevin.h" namespace LAMMPS_NS { -class FixTempRescale : public Fix { +class FixLangevinEff : public FixLangevin { public: - FixTempRescale(class LAMMPS *, int, char **); - ~FixTempRescale(); - int setmask(); - void init(); + FixLangevinEff(class LAMMPS *, int, char **); + ~FixLangevinEff(); void end_of_step(); - int modify_param(int, char **); - void reset_target(double); double compute_scalar(); private: - int which; - double t_start,t_stop,t_window; - double fraction,energy,efactor; + double *erforcelangevin; - char *id_temp; - class Compute *temperature; - int tflag; + void post_force_no_tally(); + void post_force_tally(); }; } #endif #endif diff --git a/src/USER-EFF/fix_nh_eff.cpp b/src/USER-EFF/fix_nh_eff.cpp new file mode 100644 index 000000000..0a30b02ee --- /dev/null +++ b/src/USER-EFF/fix_nh_eff.cpp @@ -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; +} diff --git a/src/fix_temp_rescale.h b/src/USER-EFF/fix_nh_eff.h similarity index 55% copy from src/fix_temp_rescale.h copy to src/USER-EFF/fix_nh_eff.h index f710bc7cd..c51b6be01 100644 --- a/src/fix_temp_rescale.h +++ b/src/USER-EFF/fix_nh_eff.h @@ -1,51 +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. ------------------------------------------------------------------------- */ -#ifdef FIX_CLASS +#ifndef LMP_FIX_NH_EFF_H +#define LMP_FIX_NH_EFF_H -FixStyle(temp/rescale,FixTempRescale) - -#else - -#ifndef LMP_FIX_TEMP_RESCALE_H -#define LMP_FIX_TEMP_RESCALE_H - -#include "fix.h" +#include "fix_nh.h" namespace LAMMPS_NS { -class FixTempRescale : public Fix { +class FixNHEff : public FixNH { public: - FixTempRescale(class LAMMPS *, int, char **); - ~FixTempRescale(); - int setmask(); - void init(); - void end_of_step(); - int modify_param(int, char **); - void reset_target(double); - double compute_scalar(); - - private: - int which; - double t_start,t_stop,t_window; - double fraction,energy,efactor; - - char *id_temp; - class Compute *temperature; - int tflag; + FixNHEff(class LAMMPS *, int, char **); + virtual ~FixNHEff() {} + + protected: + void nve_v(); + void nve_x(); + virtual void nh_v_temp(); }; } #endif -#endif diff --git a/src/USER-EFF/fix_nph_eff.cpp b/src/USER-EFF/fix_nph_eff.cpp new file mode 100644 index 000000000..98f76c346 --- /dev/null +++ b/src/USER-EFF/fix_nph_eff.cpp @@ -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; +} diff --git a/src/fix_temp_rescale.h b/src/USER-EFF/fix_nph_eff.h similarity index 57% copy from src/fix_temp_rescale.h copy to src/USER-EFF/fix_nph_eff.h index f710bc7cd..04d265e7b 100644 --- a/src/fix_temp_rescale.h +++ b/src/USER-EFF/fix_nph_eff.h @@ -1,51 +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(temp/rescale,FixTempRescale) +FixStyle(nph/sphere,FixNPHEff) #else -#ifndef LMP_FIX_TEMP_RESCALE_H -#define LMP_FIX_TEMP_RESCALE_H +#ifndef LMP_FIX_NPH_EFF_H +#define LMP_FIX_NPH_EFF_H -#include "fix.h" +#include "fix_nh_eff.h" namespace LAMMPS_NS { -class FixTempRescale : public Fix { +class FixNPHEff : public FixNHEff { public: - FixTempRescale(class LAMMPS *, int, char **); - ~FixTempRescale(); - int setmask(); - void init(); - void end_of_step(); - int modify_param(int, char **); - void reset_target(double); - double compute_scalar(); - - private: - int which; - double t_start,t_stop,t_window; - double fraction,energy,efactor; - - char *id_temp; - class Compute *temperature; - int tflag; + FixNPHEff(class LAMMPS *, int, char **); + ~FixNPHEff() {} }; } #endif #endif diff --git a/src/USER-EFF/fix_npt_eff.cpp b/src/USER-EFF/fix_npt_eff.cpp new file mode 100644 index 000000000..4c2b381fa --- /dev/null +++ b/src/USER-EFF/fix_npt_eff.cpp @@ -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; +} diff --git a/src/fix_temp_rescale.h b/src/USER-EFF/fix_npt_eff.h similarity index 57% copy from src/fix_temp_rescale.h copy to src/USER-EFF/fix_npt_eff.h index f710bc7cd..ace91a596 100644 --- a/src/fix_temp_rescale.h +++ b/src/USER-EFF/fix_npt_eff.h @@ -1,51 +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(temp/rescale,FixTempRescale) +FixStyle(npt/eff,FixNPTEff) #else -#ifndef LMP_FIX_TEMP_RESCALE_H -#define LMP_FIX_TEMP_RESCALE_H +#ifndef LMP_FIX_NPT_EFF_H +#define LMP_FIX_NPT_EFF_H -#include "fix.h" +#include "fix_nh_eff.h" namespace LAMMPS_NS { -class FixTempRescale : public Fix { +class FixNPTEff : public FixNHEff { public: - FixTempRescale(class LAMMPS *, int, char **); - ~FixTempRescale(); - int setmask(); - void init(); - void end_of_step(); - int modify_param(int, char **); - void reset_target(double); - double compute_scalar(); - - private: - int which; - double t_start,t_stop,t_window; - double fraction,energy,efactor; - - char *id_temp; - class Compute *temperature; - int tflag; + FixNPTEff(class LAMMPS *, int, char **); + ~FixNPTEff() {} }; } #endif #endif diff --git a/src/USER-EFF/fix_nve_eff.cpp b/src/USER-EFF/fix_nve_eff.cpp new file mode 100644 index 000000000..0f20bec33 --- /dev/null +++ b/src/USER-EFF/fix_nve_eff.cpp @@ -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; + } + } + } +} diff --git a/src/fix_temp_rescale.h b/src/USER-EFF/fix_nve_eff.h similarity index 57% copy from src/fix_temp_rescale.h copy to src/USER-EFF/fix_nve_eff.h index f710bc7cd..33d437961 100644 --- a/src/fix_temp_rescale.h +++ b/src/USER-EFF/fix_nve_eff.h @@ -1,51 +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,FixTempRescale) +FixStyle(nve/eff,FixNVEEff) #else -#ifndef LMP_FIX_TEMP_RESCALE_H -#define LMP_FIX_TEMP_RESCALE_H +#ifndef LMP_FIX_NVE_EFF_H +#define LMP_FIX_NVE_EFF_H -#include "fix.h" +#include "fix_nve.h" namespace LAMMPS_NS { -class FixTempRescale : public Fix { +class FixNVEEff : public FixNVE { public: - FixTempRescale(class LAMMPS *, int, char **); - ~FixTempRescale(); - int setmask(); - void init(); - void end_of_step(); - int modify_param(int, char **); - void reset_target(double); - double compute_scalar(); - - private: - int which; - double t_start,t_stop,t_window; - double fraction,energy,efactor; - - char *id_temp; - class Compute *temperature; - int tflag; + FixNVEEff(class LAMMPS *, int, char **); + void initial_integrate(int); + void final_integrate(); }; } #endif #endif diff --git a/src/USER-EFF/fix_nvt_eff.cpp b/src/USER-EFF/fix_nvt_eff.cpp new file mode 100644 index 000000000..b742b66ee --- /dev/null +++ b/src/USER-EFF/fix_nvt_eff.cpp @@ -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; +} diff --git a/src/fix_temp_rescale.h b/src/USER-EFF/fix_nvt_eff.h similarity index 57% copy from src/fix_temp_rescale.h copy to src/USER-EFF/fix_nvt_eff.h index f710bc7cd..3e832094c 100644 --- a/src/fix_temp_rescale.h +++ b/src/USER-EFF/fix_nvt_eff.h @@ -1,51 +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(temp/rescale,FixTempRescale) +FixStyle(nvt/eff,FixNVTEff) #else -#ifndef LMP_FIX_TEMP_RESCALE_H -#define LMP_FIX_TEMP_RESCALE_H +#ifndef LMP_FIX_NVT_EFF_H +#define LMP_FIX_NVT_EFF_H -#include "fix.h" +#include "fix_nh_eff.h" namespace LAMMPS_NS { -class FixTempRescale : public Fix { +class FixNVTEff : public FixNHEff { public: - FixTempRescale(class LAMMPS *, int, char **); - ~FixTempRescale(); - int setmask(); - void init(); - void end_of_step(); - int modify_param(int, char **); - void reset_target(double); - double compute_scalar(); - - private: - int which; - double t_start,t_stop,t_window; - double fraction,energy,efactor; - - char *id_temp; - class Compute *temperature; - int tflag; + FixNVTEff(class LAMMPS *, int, char **); + virtual ~FixNVTEff() {} }; } #endif #endif diff --git a/src/USER-EFF/fix_nvt_sllod_eff.cpp b/src/USER-EFF/fix_nvt_sllod_eff.cpp new file mode 100644 index 000000000..792208cdc --- /dev/null +++ b/src/USER-EFF/fix_nvt_sllod_eff.cpp @@ -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)); + } + } +} diff --git a/src/fix_temp_rescale.h b/src/USER-EFF/fix_nvt_sllod_eff.h similarity index 59% copy from src/fix_temp_rescale.h copy to src/USER-EFF/fix_nvt_sllod_eff.h index f710bc7cd..e8f9007cf 100644 --- a/src/fix_temp_rescale.h +++ b/src/USER-EFF/fix_nvt_sllod_eff.h @@ -1,51 +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(temp/rescale,FixTempRescale) +FixStyle(nvt/sllod/eff,FixNVTSllodEff) #else -#ifndef LMP_FIX_TEMP_RESCALE_H -#define LMP_FIX_TEMP_RESCALE_H +#ifndef LMP_FIX_NVT_SLODD_EFF_H +#define LMP_FIX_NVT_SLODD_EFF_H -#include "fix.h" +#include "fix_nvt_eff.h" namespace LAMMPS_NS { -class FixTempRescale : public Fix { +class FixNVTSllodEff : public FixNVTEff { public: - FixTempRescale(class LAMMPS *, int, char **); - ~FixTempRescale(); - int setmask(); + FixNVTSllodEff(class LAMMPS *, int, char **); + ~FixNVTSllodEff() {} void init(); - void end_of_step(); - int modify_param(int, char **); - void reset_target(double); - double compute_scalar(); private: - int which; - double t_start,t_stop,t_window; - double fraction,energy,efactor; + int nondeformbias; - char *id_temp; - class Compute *temperature; - int tflag; + void nh_v_temp(); }; } #endif #endif diff --git a/src/USER-EFF/fix_temp_rescale_eff.cpp b/src/USER-EFF/fix_temp_rescale_eff.cpp new file mode 100644 index 000000000..9a31b1023 --- /dev/null +++ b/src/USER-EFF/fix_temp_rescale_eff.cpp @@ -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]); + } + } + } + + } +} diff --git a/src/fix_temp_rescale.h b/src/USER-EFF/fix_temp_rescale_eff.h similarity index 59% copy from src/fix_temp_rescale.h copy to src/USER-EFF/fix_temp_rescale_eff.h index f710bc7cd..76e7700a6 100644 --- a/src/fix_temp_rescale.h +++ b/src/USER-EFF/fix_temp_rescale_eff.h @@ -1,51 +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,FixTempRescale) +FixStyle(temp/rescale/eff,FixTempRescaleEff) #else -#ifndef LMP_FIX_TEMP_RESCALE_H -#define LMP_FIX_TEMP_RESCALE_H +#ifndef LMP_FIX_TEMP_RESCALE_EFF_H +#define LMP_FIX_TEMP_RESCALE_EFF_H -#include "fix.h" +#include "fix_temp_rescale.h" namespace LAMMPS_NS { -class FixTempRescale : public Fix { +class FixTempRescaleEff : public FixTempRescale { public: - FixTempRescale(class LAMMPS *, int, char **); - ~FixTempRescale(); - int setmask(); - void init(); + FixTempRescaleEff(class LAMMPS *, int, char **); + ~FixTempRescaleEff() {} void end_of_step(); - int modify_param(int, char **); - void reset_target(double); - double compute_scalar(); - - private: - int which; - double t_start,t_stop,t_window; - double fraction,energy,efactor; - - char *id_temp; - class Compute *temperature; - int tflag; }; } #endif #endif diff --git a/src/USER-EFF/pair_eff_cut.cpp b/src/USER-EFF/pair_eff_cut.cpp new file mode 100644 index 000000000..013f8ecec --- /dev/null +++ b/src/USER-EFF/pair_eff_cut.cpp @@ -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; +} diff --git a/src/USER-EFF/pair_eff_cut.h b/src/USER-EFF/pair_eff_cut.h new file mode 100644 index 000000000..e47a13f82 --- /dev/null +++ b/src/USER-EFF/pair_eff_cut.h @@ -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 diff --git a/src/USER-EFF/pair_eff_inline.h b/src/USER-EFF/pair_eff_inline.h new file mode 100644 index 000000000..b21e35188 --- /dev/null +++ b/src/USER-EFF/pair_eff_inline.h @@ -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; +} + +} diff --git a/src/fix_temp_rescale.h b/src/fix_temp_rescale.h index f710bc7cd..bd87e516b 100644 --- a/src/fix_temp_rescale.h +++ b/src/fix_temp_rescale.h @@ -1,51 +1,51 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #ifdef FIX_CLASS FixStyle(temp/rescale,FixTempRescale) #else #ifndef LMP_FIX_TEMP_RESCALE_H #define LMP_FIX_TEMP_RESCALE_H #include "fix.h" 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; char *id_temp; class Compute *temperature; int tflag; }; } #endif #endif