Homec4science
Diffusion lammps (lammm-master-deprecated) Browse

lammps/src/USER-QTB/fix_qtb.cpp
1ba7f726cce2lammm-master-deprecated

Last Change
Raw File

fix_qtb.cpp
View Options

/* ----------------------------------------------------------------------
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;
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;
}
c4science · Help