ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
sjplimp
2015-07-22 15:49:58 +00:00
committed by Pierre de Buyl
parent 27def043d9
commit 161f0205b2
5 changed files with 1876 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

56
src/USER-QTB/README Normal file
View File

@ -0,0 +1,56 @@
This package contains 2 fix commands, "fix qtb" and
"fix qbmsst" which involve quantum nuclear effects.
What is quantum nuclear effects (in molecular dynamics
simulation)?
Quantum treatment of the vibrational modes will
introduce zero point energy into the system, alter
the energy power spectrum and bias the heat capacity
from the classical limit. classical MD leaves out
these effects completely and thus need to be corrected.
When should I consider quantum nuclear effects in
my simulation?
(1) When you want to model systems at temperatures
lower than their classical limits. This is especially
important for materials with a large population of
hydrogen atoms and thus higher classical limits.
(2) In MD simulations when temperatures ramp up
across (or starting from somewhere close to) the
classical limits, tiny differences between the
classical and quantum heat capacity can accumulate
and cause the final state to deviate (could be as
large as 30 %).
What command should I use regarding these two cases?
(1) "fix qtb" provides quantum nulcear correction
through a colored thermostat and can be used with
other time integration schemes like fix nve or
fix nph. In this case, you tell "fix qtb" the temperature
and the other time integration scheme either the
pressure ("fix nph") or the volume ("fix nve").
(2) "fix qbmsst" deals with shock MD simulations
in which cases the temperature may increase a lot.
It enables quantum nuclear correction of a multi-scale
shock technique simulation by coupling the quantum
thermal bath with the shocked system. It is an independent
command from "fix msst" and does not assume that
the SHOCK package has been installed.
See the doc page for the fix qtb and fix qbmsst
command for detailed usage instructions.
Where can I find examples of these two commands?
There are example scripts for using this package in
examples/USER/qtb, including one "fix qtb" and one
"fix qbmsst" example for each of alpha quartz and
methane. Running the alpha quartz example requires
installation of the kspace package while the methane
example requires the REAX package for the force field.
Authors Information and Contact
The person who created this package is Yuan Shen
(sy0302 at stanford.edu) at Stanford University.
Contact him directly if you have questions.

1195
src/USER-QTB/fix_qbmsst.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

124
src/USER-QTB/fix_qbmsst.h Normal file
View File

@ -0,0 +1,124 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Shen,Yuan, Qi,Tingting, and Reed,Evan
Implementation of the Multi-Scale Shock Method with quantum nuclear effects
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(qbmsst,FixQBMSST)
#else
#ifndef FIX_QBMSST_H
#define FIX_QBMSST_H
#include "fix.h"
namespace LAMMPS_NS {
class FixQBMSST : public Fix {
public:
FixQBMSST(class LAMMPS *, int, char **);
~FixQBMSST();
int setmask();
void init();
void setup(int);
void initial_integrate(int);
void final_integrate();
double compute_scalar();
double compute_vector(int);
void write_restart(FILE *);
void restart(char *);
int modify_param(int, char **);
double memory_usage();
void grow_arrays(int);
void copy_arrays(int, int, int);
int pack_exchange(int, double *);
int unpack_exchange(int, double *);
private:
// msst parameters
int direction; // Direction of shock
double velocity; // Velocity of the shock
double qmass; // Effective cell mass
double mu; // Effective cell viscosity
double e0; // Initial energy
double v0; // Initial volume
double p0; // Initial pressure
int p0_set; // Is pressure set
int v0_set; // Is volume set
int e0_set; // Is energy set
double tscale; // Converts thermal energy to compressive strain ke at simulation start
char *id_temp,*id_press, *id_pe; // Strings with identifiers of created computes.
int tflag,pflag,vsflag,peflag; // Flags to keep track of computes that were created.
double dtv; // update->dt
double dtf; // Full step size
double dthalf; // half step size
bigint ntotal; // atom->natoms
double boltz,nktv2p, mvv2e; // Boltzmann factor and unit conversions
class Compute *temperature; // Computes created to evaluate
class Compute *pressure; // thermodynamic quantities.
class Compute *pe;
double **old_velocity; // Saved velocities.
int atoms_allocated; // size of old_velocity
double dilation[3];
double omega[3]; // Time derivative of the volume.
double total_mass; // Mass of the computational cell
int kspace_flag; // 1 if KSpace invoked, 0 if not
int nrigid; // number of rigid fixes
int *rfix; // indices of rigid fixes
double p_current[3]; // pressure
double velocity_sum; // Sum of the velocities squared.
double lagrangian_position; // Lagrangian location of computational cell
// qbmsst parameters
double t_period, fric_coef; // friction coefficient
int seed; // seed for the random number generator
double f_max; // frequency cutoff
int N_f; // number of frequency grid
double eta; // coupling coefficient bewteen shock and the qtb
int beta; // average beta steps before updating the qtb temperature
double t_init; // initial qtb temperature
int qtb_set; // 1 if its a restarting qbmsst, 0 if not
int counter_l, counter_mu; // counter l and mu
double t_current; // qtb temperature
double h_timestep; // time step to update the random forces
int alpha; // number of time steps to update the random forces
class RanMars *random; // random number generator
double *gfactor; // factors of random forces
double *omega_H,*time_H; // H gives the desired power spectrum
double **random_array_0, **random_array_1, **random_array_2; // random number arrays give independence between atoms and directions
double **fran; // random forces
double old_eavg; // average previous energies
// functions
void couple();
void remap(int);
void check_alloc(int n);
double compute_etotal();
double compute_egrand();
double compute_vol();
double compute_vsum();
double compute_hugoniot();
double compute_rayleigh();
double compute_lagrangian_speed();
double compute_lagrangian_position();
};
}
#endif
#endif

430
src/USER-QTB/fix_qtb.cpp Normal file
View File

@ -0,0 +1,430 @@
/* ----------------------------------------------------------------------
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: Shen,Yuan, Qi,Tingting, and Reed,Evan
Implementation of the colored thermostat for quantum nuclear effects
------------------------------------------------------------------------- */
#include "mpi.h"
#include "math.h"
#include "string.h"
#include "stdlib.h"
#include "fix_qtb.h"
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
#include "force.h"
#include "update.h"
#include "modify.h"
#include "compute.h"
#include "domain.h"
#include "region.h"
#include "respa.h"
#include "comm.h"
#include "input.h"
#include "variable.h"
#include "random_mars.h"
#include "memory.h"
#include "error.h"
#include "group.h"
using namespace LAMMPS_NS;
using namespace FixConst;
/* ----------------------------------------------------------------------
read parameters
------------------------------------------------------------------------- */
FixQTB::FixQTB(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 3) error->all(FLERR,"Illegal fix qtb command");
// default parameters
global_freq = 1;
extscalar = 1;
nevery = 1;
t_target = 300.0;
t_period = 1.0;
fric_coef = 1/t_period;
seed = 880302;
f_max = 200.0;
N_f = 100;
// reading parameters
int iarg = 3;
while (iarg < narg) {
if (strcmp(arg[iarg],"temp") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix qtb command");
t_target = atof(arg[iarg+1]); if (t_target < 0.0) error->all(FLERR,"Fix qtb temp must be >= 0.0");
iarg += 2;
} else if (strcmp(arg[iarg],"damp") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix qtb command");
t_period = atof(arg[iarg+1]); if (t_period <= 0.0) error->all(FLERR,"Fix qtb damp must be > 0.0"); fric_coef = 1/t_period;
iarg += 2;
} else if (strcmp(arg[iarg],"seed") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix qtb command");
seed = atof(arg[iarg+1]); if (seed <= 0) error->all(FLERR,"Illegal fix qtb command");
iarg += 2;
} else if (strcmp(arg[iarg],"f_max") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix qtb command");
f_max = atof(arg[iarg+1]); if (f_max <= 0) error->all(FLERR,"Illegal fix qtb command");
iarg += 2;
} else if (strcmp(arg[iarg],"N_f") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix qtb command");
N_f = atof(arg[iarg+1]); if (N_f <= 0) error->all(FLERR,"Illegal fix qtb command");
iarg += 2;
} else error->all(FLERR,"Illegal fix qtb command");
}
// allocate qtb
gfactor1 = NULL;
gfactor3 = NULL;
omega_H = NULL;
time_H = NULL;
random_array_0 = NULL;
random_array_1 = NULL;
random_array_2 = NULL;
fran = NULL;
id_temp = NULL;
temperature = NULL;
// initialize Marsaglia RNG with processor-unique seed
random = new RanMars(lmp,seed + comm->me);
// allocate per-type arrays for force prefactors
gfactor1 = new double[atom->ntypes+1];
gfactor3 = new double[atom->ntypes+1];
// allocate random-arrays and fran
grow_arrays(atom->nmax);
atom->add_callback(0);
// memory->create(random_array_0,atom->nmax+300,2*N_f,"qtb:random_array_0");
// memory->create(random_array_1,atom->nmax+300,2*N_f,"qtb:random_array_1");
// memory->create(random_array_2,atom->nmax+300,2*N_f,"qtb:random_array_2");
// memory->create(fran,atom->nmax+300,3,"qtb:fran");
// allocate omega_H and time_H
memory->create(omega_H,2*N_f,"qtb:omega_H");
memory->create(time_H,2*N_f,"qtb:time_H");
}
/* ----------------------------------------------------------------------
release memories
------------------------------------------------------------------------- */
FixQTB::~FixQTB()
{
delete random;
delete [] gfactor1;
delete [] gfactor3;
delete [] id_temp;
memory->destroy(fran);
memory->destroy(random_array_0);
memory->destroy(random_array_1);
memory->destroy(random_array_2);
memory->destroy(omega_H);
memory->destroy(time_H);
atom->delete_callback(id,0);
}
/* ----------------------------------------------------------------------
setmask
------------------------------------------------------------------------- */
int FixQTB::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= THERMO_ENERGY;
return mask;
}
/* ----------------------------------------------------------------------
fix initiation
------------------------------------------------------------------------- */
void FixQTB::init()
{
// copy parameters from other classes
double dtv = update->dt;
if (atom->mass == NULL)
error->all(FLERR,"Cannot use fix msst without per-type mass defined");
//initiate the counter \mu
counter_mu=0;
//set up the h time step for updating the random force \delta{}h=\frac{\pi}{\Omega_{max}}
if (int(1.0/(2*f_max*dtv)) == 0) {
if (comm->me == 0) printf ("Warning: Either f_max is too high or the time step is too big, setting f_max to be 1/timestep!\n");
h_timestep=dtv;
alpha=1;
} else {
alpha=int(1.0/(2*f_max*dtv));
h_timestep=alpha*dtv;
}
if (comm->me == 0) printf ("The effective maximum frequncy is now %f inverse time unit with alpha value as %d!\n", 0.5/h_timestep, alpha);
// set force prefactors
if (!atom->rmass) {
for (int i = 1; i <= atom->ntypes; i++) {
//gfactor1 is the friction force \gamma{}m_{i}\frac{dv}{dt}
gfactor1[i] = (atom->mass[i]*fric_coef) / force->ftm2v;
//gfactor3 is the random force \sqrt{\frac{2\gamma{}m_{i}}{\alpha*\delta{}t}}, \sqrt{12} makes the random array variance equal to unit.
gfactor3[i] = sqrt(2*fric_coef*atom->mass[i])*sqrt(force->mvv2e)*sqrt(12/h_timestep); //this still leaves a square energy term from the power spectrum H.
}
}
// generate random number array with zero mean and variance equal 1/12.
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
for (int m = 0; m < 2*N_f; m++) {
random_array_0[i][m] = random->uniform()-0.5;
random_array_1[i][m] = random->uniform()-0.5;
random_array_2[i][m] = random->uniform()-0.5;
}
}
// load omega_H with calculated spectrum at a specific temperature (corrected spectrum), omega_H is the Fourier transformation of time_H
for (int k = 0; k < 2*N_f; k++) {
double f_k=(k-N_f)/(2*N_f*h_timestep); //\omega_k=\frac{2\pi}{\delta{}h}\frac{k}{2N_f} for k from -N_f to N_f-1
if(k == N_f) {
omega_H[k]=sqrt(force->boltz * t_target);
} else {
double energy_k= force->hplanck * fabs(f_k);
omega_H[k]=sqrt( energy_k * (0.5+1.0/( exp(energy_k/(force->boltz * t_target)) - 1.0 )) );
omega_H[k]*=alpha*sin((k-N_f)*M_PI/(2*alpha*N_f))/sin((k-N_f)*M_PI/(2*N_f));
}
}
// construct the signal filter H, filter has the unit of of sqrt(energy) \sqrt{2N_f}^{-1}H\left(t_n\right)
for (int n = 0; n < 2*N_f; n++) {
time_H[n] = 0;
double t_n=(n-N_f);
for (int k = 0; k < 2*N_f; k++) {
double omega_k=(k-N_f)*M_PI/N_f;
time_H[n] += omega_H[k]*(cos(omega_k*t_n));
}
time_H[n]/=(2.0*N_f);
}
// respa
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ----------------------------------------------------------------------
no MD, so setup returns post force
------------------------------------------------------------------------- */
void FixQTB::setup(int vflag)
{
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
}
}
/* ----------------------------------------------------------------------
post_force
------------------------------------------------------------------------- */
void FixQTB::post_force(int vflag)
{
double gamma1,gamma3;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
int *type = atom->type;
int *mask = atom->mask;
bigint nlocal = atom->nlocal;
bigint ntotal = atom->natoms;
//update the colored random force every alpha MD steps
if (counter_mu == alpha) {
//propagate h_timestep ahead
for (int j = 0; j < nlocal; j++) {
//update random array
for (int m = 0; m < 2*N_f-1; m++) {
random_array_0[j][m] = random_array_0[j][m+1];
random_array_1[j][m] = random_array_1[j][m+1];
random_array_2[j][m] = random_array_2[j][m+1];
}
random_array_0[j][2*N_f-1] = random->uniform()-0.5;
random_array_1[j][2*N_f-1] = random->uniform()-0.5;
random_array_2[j][2*N_f-1] = random->uniform()-0.5;
}
//reset counter \mu
counter_mu=0;
}
if (counter_mu == 0) {
for (int j = 0; j < nlocal; j++) {
fran[j][0] = 0.0;
fran[j][1] = 0.0;
fran[j][2] = 0.0;
//reset random force
if (mask[j] & groupbit) {
gamma3 = gfactor3[type[j]];
for (int m = 0; m < 2*N_f; m++) {
fran[j][0] += time_H[m] * random_array_0[j][2*N_f-m-1];
fran[j][1] += time_H[m] * random_array_1[j][2*N_f-m-1];
fran[j][2] += time_H[m] * random_array_2[j][2*N_f-m-1];
}
fran[j][0] = fran[j][0]*gamma3;
fran[j][1] = fran[j][1]*gamma3;
fran[j][2] = fran[j][2]*gamma3;
}
}
}
//reset all the force sums
fsum[0]=0.0; fsumall[0]=0.0;
fsum[1]=0.0; fsumall[1]=0.0;
fsum[2]=0.0; fsumall[2]=0.0;
for (int j = 0; j < nlocal; j++) {
//sum over each atom
if (mask[j] & groupbit) {
gamma1 = gfactor1[type[j]];
fsum[0]+=fran[j][0]-gamma1*v[j][0];
fsum[1]+=fran[j][1]-gamma1*v[j][1];
fsum[2]+=fran[j][2]-gamma1*v[j][2];
}
}
//compute force sums
MPI_Allreduce(fsum,fsumall,3,MPI_DOUBLE,MPI_SUM,world);
//implement random forces
for (int j = 0; j < nlocal; j++) {
//make sure there is no net force on the system
f[j][0] -= fsumall[0]/ntotal;
f[j][1] -= fsumall[1]/ntotal;
f[j][2] -= fsumall[2]/ntotal;
if (mask[j] & groupbit) {
gamma1 = gfactor1[type[j]];
f[j][0]+=fran[j][0]-gamma1*v[j][0];
f[j][1]+=fran[j][1]-gamma1*v[j][1];
f[j][2]+=fran[j][2]-gamma1*v[j][2];
}
}
//move 1 step forward
counter_mu++;
}
/* ----------------------------------------------------------------------
post_force_respa
------------------------------------------------------------------------- */
void FixQTB::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
}
/* ----------------------------------------------------------------------
modifications of fix qtb
------------------------------------------------------------------------- */
int FixQTB::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"temp") == 0) {
if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
int icompute = modify->find_compute(id_temp);
if (icompute < 0) error->all(FLERR,"Could not find fix_modify temperature ID");
temperature = modify->compute[icompute];
if (temperature->tempflag == 0)
error->all(FLERR,"Fix_modify temperature ID does not compute temperature");
if (temperature->igroup != igroup && comm->me == 0)
error->warning(FLERR,"Group for fix_modify temp != fix group");
return 2;
}
return 0;
}
/* ----------------------------------------------------------------------
memory usage of fix qtb
------------------------------------------------------------------------- */
double FixQTB::memory_usage()
{
double bytes = 0.0;
// random_arrays memory usage
bytes += (atom->nmax* 6*N_f * sizeof(double));
// fran memory usage
bytes += (atom->nmax* 3 * sizeof(double));
bytes += (4*N_f * sizeof(double));
return bytes;
}
/* ----------------------------------------------------------------------
allocate atom-based array for fran and random_array
------------------------------------------------------------------------- */
void FixQTB::grow_arrays(int nmax)
{
memory->grow(random_array_0,nmax,2*N_f,"qtb:random_array_0");
memory->grow(random_array_1,nmax,2*N_f,"qtb:random_array_1");
memory->grow(random_array_2,nmax,2*N_f,"qtb:random_array_2");
memory->grow(fran,nmax,3,"qtb:fran");
}
/* ----------------------------------------------------------------------
copy values within local atom-based array
------------------------------------------------------------------------- */
void FixQTB::copy_arrays(int i, int j, int delflag)
{
for (int m = 0; m < 2*N_f; m++) {
random_array_0[j][m] = random_array_0[i][m];
random_array_1[j][m] = random_array_1[i][m];
random_array_2[j][m] = random_array_2[i][m];
}
for (int m = 0; m < 3; m++)
fran[j][m] = fran[i][m];
}
/* ----------------------------------------------------------------------
pack values in local atom-based array for exchange with another proc
------------------------------------------------------------------------- */
int FixQTB::pack_exchange(int i, double *buf)
{
for (int m = 0; m < 2*N_f; m++) buf[m] = random_array_0[i][m];
for (int m = 0; m < 2*N_f; m++) buf[m+2*N_f] = random_array_1[i][m];
for (int m = 0; m < 2*N_f; m++) buf[m+4*N_f] = random_array_2[i][m];
for (int m = 0; m < 3; m++) buf[m+6*N_f] = fran[i][m];
return 6*N_f+3;
}
/* ----------------------------------------------------------------------
unpack values in local atom-based array from exchange with another proc
------------------------------------------------------------------------- */
int FixQTB::unpack_exchange(int nlocal, double *buf)
{
for (int m = 0; m < 2*N_f; m++) random_array_0[nlocal][m] = buf[m];
for (int m = 0; m < 2*N_f; m++) random_array_1[nlocal][m] = buf[m+2*N_f];
for (int m = 0; m < 2*N_f; m++) random_array_2[nlocal][m] = buf[m+4*N_f];
for (int m = 0; m < 3; m++) fran[nlocal][m] = buf[m+6*N_f];
return 6*N_f+3;
}

71
src/USER-QTB/fix_qtb.h Normal file
View File

@ -0,0 +1,71 @@
/* ----------------------------------------------------------------------
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: Shen,Yuan, Qi,Tingting, and Reed,Evan
Implementation of the colored thermostat for quantum nuclear effects
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(qtb,FixQTB)
#else
#ifndef LMP_FIX_QTB_H
#define LMP_FIX_QTB_H
#include "fix.h"
namespace LAMMPS_NS {
class FixQTB : public Fix {
public:
FixQTB(class LAMMPS *, int, char **);
virtual ~FixQTB();
int setmask();
void init();
void setup(int);
void post_force(int);
void post_force_respa(int, int, int);
int modify_param(int, char **);
double memory_usage();
void grow_arrays(int);
void copy_arrays(int, int, int);
int pack_exchange(int, double *);
int unpack_exchange(int, double *);
private:
// qtb parameters
int counter_mu; // counter l and mu
double t_period, fric_coef; // friction coefficient
int seed; // seed for the random number generator
double f_max; // frequency cutoff
int N_f; // number of frequency grid
double t_target; // target qtb temperature
char *id_temp;
class Compute *temperature;
double h_timestep; // time step to update the random forces
int alpha; // number of time steps to update the random forces
class RanMars *random; // random number generator
double *gfactor1,*gfactor3; // factors of frictions and random forces
double *omega_H,*time_H; // H gives the desired power spectrum
double **random_array_0, **random_array_1, **random_array_2; // random number arrays give independence between atoms and directions
int nlevels_respa;
double **fran, fsum[3], fsumall[3]; // random forces and their sums
};
}
#endif
#endif