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@12788 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2014-11-24 20:57:12 +00:00
parent 0e4d24ade8
commit 1f84af1b91
4 changed files with 977 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
src/Makefile
View File

@ -19,12 +19,10 @@ PACKAGE = asphere body class2 colloid dipole fld gpu granular kim \
PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars \
user-cuda user-eff user-fep user-intel user-lb user-misc \
user-molfile user-omp user-phonon user-pimd \
user-qmmm user-reaxc user-sph
user-molfile user-omp user-phonon user-qmmm user-reaxc user-sph
PACKLIB = gpu kim meam poems reax voronoi \
user-atc user-awpmd user-colvars user-pimd \
user-qmmm user-cuda user-molfile
user-atc user-awpmd user-colvars user-qmmm user-cuda user-molfile
PACKALL = $(PACKAGE) $(PACKUSER)

1
src/USER-MISC/README
View File

@ -39,6 +39,7 @@ fix ave/spatial/sphere, Niall Jackson (Imperial College), niall.jackson at gmail
fix gle, Michele Ceriotti (EPFL Lausanne), michele.ceriotti at gmail.com, 24 Nov 2014
fix imd, Axel Kohlmeyer, akohlmey at gmail.com, 9 Nov 2009
fix ipi, Michele Ceriotti (EPFL Lausanne), michele.ceriotti at gmail.com, 24 Nov 2014
fix pimd, Yuxing Peng (U Chicago), yuxing at uchicago.edu, 24 Nov 2014
fix smd, Axel Kohlmeyer, akohlmey at gmail.com, 19 May 2008
fix ti/rs, Rodrigo Freitas (Unicamp/Brazil), rodrigohb at gmail.com, 7 Nov 2013
fix ti/spring, Rodrigo Freitas (Unicamp/Brazil), rodrigohb at gmail.com, 7 Nov 2013

859
src/USER-MISC/fix_pimd.cpp Normal file
View File

@ -0,0 +1,859 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Package FixPIMD
Purpose Quantum Path Integral Algorithm for Quantum Chemistry
Copyright Voth Group @ University of Chicago
Authors Chris Knight & Yuxing Peng
Updated Oct-01-2011
Version 1.0
------------------------------------------------------------------------- */
#include "math.h"
#include "string.h"
#include "stdlib.h"
#include "fix_pimd.h"
#include "universe.h"
#include "comm.h"
#include "force.h"
#include "atom.h"
#include "domain.h"
#include "update.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum{PIMD,NMPIMD,CMD};
/* ---------------------------------------------------------------------- */
FixPIMD::FixPIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
method = PIMD;
fmass = 1.0;
nhc_temp = 298.15;
nhc_nchain = 2;
sp = 1.0;
for(int i=3; i<narg-1; i+=2)
{
if(strcmp(arg[i],"method")==0)
{
if(strcmp(arg[i+1],"pimd")==0) method=PIMD;
else if(strcmp(arg[i+1],"nmpimd")==0) method=NMPIMD;
else if(strcmp(arg[i+1],"cmd")==0) method=CMD;
else error->universe_all(FLERR,"Unkown method parameter for fix pimd");
}
else if(strcmp(arg[i],"fmass")==0)
{
fmass = atof(arg[i+1]);
if(fmass<0.0 || fmass>1.0) error->universe_all(FLERR,"Invalid fmass value for fix pimd");
}
else if(strcmp(arg[i],"sp")==0)
{
sp = atof(arg[i+1]);
if(fmass<0.0) error->universe_all(FLERR,"Invalid sp value for fix pimd");
}
else if(strcmp(arg[i],"temp")==0)
{
nhc_temp = atof(arg[i+1]);
if(nhc_temp<0.0) error->universe_all(FLERR,"Invalid temp value for fix pimd");
}
else if(strcmp(arg[i],"nhc")==0)
{
nhc_nchain = atoi(arg[i+1]);
if(nhc_nchain<2) error->universe_all(FLERR,"Invalid nhc value for fix pimd");
}
else error->universe_all(arg[i],i+1,"Unkown keyword for fix pimd");
}
/* Initiation */
max_nsend = 0;
tag_send = NULL;
buf_send = NULL;
max_nlocal = 0;
buf_recv = NULL;
buf_beads = NULL;
size_plan = 0;
plan_send = plan_recv = NULL;
M_x2xp = M_xp2x = M_f2fp = M_fp2f = NULL;
lam = NULL;
mode_index = NULL;
mass = NULL;
array_atom = NULL;
nhc_eta = NULL;
nhc_eta_dot = NULL;
nhc_eta_dotdot = NULL;
nhc_eta_mass = NULL;
size_peratom_cols = 12 * nhc_nchain + 3;
nhc_offset_one_1 = 3 * nhc_nchain;
nhc_offset_one_2 = 3 * nhc_nchain +3;
nhc_size_one_1 = sizeof(double) * nhc_offset_one_1;
nhc_size_one_2 = sizeof(double) * nhc_offset_one_2;
restart_peratom = 1;
peratom_flag = 1;
peratom_freq = 1;
global_freq = 1;
thermo_energy = 1;
vector_flag = 1;
size_vector = 2;
extvector = 1;
comm_forward = 3;
atom->add_callback(0); // Call LAMMPS to allocate memory for per-atom array
atom->add_callback(1); // Call LAMMPS to re-assign restart-data for per-atom array
grow_arrays(atom->nmax);
// some initilizations
nhc_ready = false;
}
/* ---------------------------------------------------------------------- */
int FixPIMD::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixPIMD::init()
{
if(universe->me==0 && screen) fprintf(screen,"Fix pimd initializing Path-Integral ...\n");
// prepare the constants
np = universe->nworlds;
inverse_np = 1.0 / np;
/* The first solution for the force constant, using SI units
const double Boltzmann = 1.3806488E-23; // SI unit: J/K
const double Plank = 6.6260755E-34; // SI unit: m^2 kg / s
double hbar = Plank / ( 2.0 * M_PI ) * sp;
double beta = 1.0 / ( Boltzmann * input.nh_temp);
// - P / ( beta^2 * hbar^2) SI unit: s^-2
double _fbond = -1.0 / (beta*beta*hbar*hbar) * input.nbeads;
// convert the units: s^-2 -> (kcal/mol) / (g/mol) / (A^2)
fbond = _fbond * 4.184E+26;
*/
/* The current solution, using LAMMPS internal real units */
const double Boltzmann = force->boltz;
const double Plank = force->hplanck;
double hbar = Plank / ( 2.0 * M_PI );
double beta = 1.0 / (Boltzmann * nhc_temp);
double _fbond = 1.0 * np / (beta*beta*hbar*hbar) ;
omega_np = sqrt(np) / (hbar * beta) * sqrt(force->mvv2e);
fbond = - _fbond * force->mvv2e;
if(universe->me==0)
printf("Fix pimd -P/(beta^2 * hbar^2) = %20.7lE (kcal/mol/A^2)\n\n", fbond);
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
comm_init();
mass = new double [atom->ntypes+1];
if(method==CMD || method==NMPIMD) nmpimd_init();
else for(int i=1; i<=atom->ntypes; i++) mass[i] = atom->mass[i] / np * fmass;
if(!nhc_ready) nhc_init();
}
/* ---------------------------------------------------------------------- */
void FixPIMD::setup(int vflag)
{
if(universe->me==0 && screen) fprintf(screen,"Setting up Path-Integral ...\n");
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixPIMD::initial_integrate(int vflag)
{
nhc_update_v();
nhc_update_x();
}
/* ---------------------------------------------------------------------- */
void FixPIMD::final_integrate()
{
nhc_update_v();
}
/* ---------------------------------------------------------------------- */
void FixPIMD::post_force(int flag)
{
for(int i=0; i<atom->nlocal; i++) for(int j=0; j<3; j++) atom->f[i][j] /= np;
comm_exec(atom->x);
spring_force();
if(method==CMD || method==NMPIMD)
{
/* forward comm for the force on ghost atoms */
nmpimd_fill(atom->f);
/* inter-partition comm */
comm_exec(atom->f);
/* normal-mode transform */
nmpimd_transform(buf_beads, atom->f, M_f2fp[universe->iworld]);
}
}
/* ----------------------------------------------------------------------
Nose-Hoover Chains
------------------------------------------------------------------------- */
void FixPIMD::nhc_init()
{
double tau = 1.0 / omega_np;
double KT = force->boltz * nhc_temp;
double mass0 = KT * tau * tau;
int max = 3 * atom->nlocal;
for(int i=0; i<max; i++)
{
for(int ichain=0; ichain<nhc_nchain; ichain++)
{
nhc_eta[i][ichain] = 0.0;
nhc_eta_dot[i][ichain] = 0.0;
nhc_eta_dot[i][ichain] = 0.0;
nhc_eta_dotdot[i][ichain] = 0.0;
nhc_eta_mass[i][ichain] = mass0;
nhc_eta_mass[i][ichain] *= fmass;
}
nhc_eta_dot[i][nhc_nchain] = 0.0;
for(int ichain=1; ichain<nhc_nchain; ichain++)
nhc_eta_dotdot[i][ichain] = (nhc_eta_mass[i][ichain-1] * nhc_eta_dot[i][ichain-1]
* nhc_eta_dot[i][ichain-1] * force->mvv2e - KT) / nhc_eta_mass[i][ichain];
}
// Zero NH acceleration for CMD
if(method==CMD && universe->iworld==0) for(int i=0; i<max; i++)
for(int ichain=0; ichain<nhc_nchain; ichain++) nhc_eta_dotdot[i][ichain] = 0.0;
nhc_ready = true;
}
/* ---------------------------------------------------------------------- */
void FixPIMD::nhc_update_x()
{
int n = atom->nlocal;
double **x = atom->x;
double **v = atom->v;
if(method==CMD || method==NMPIMD)
{
nmpimd_fill(atom->v);
comm_exec(atom->v);
/* borrow the space of atom->f to store v in cartisian */
v = atom->f;
nmpimd_transform(buf_beads, v, M_xp2x[universe->iworld]);
}
for(int i=0; i<n; i++)
{
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
/* ---------------------------------------------------------------------- */
void FixPIMD::nhc_update_v()
{
int n = atom->nlocal;
int *type = atom->type;
double **v = atom->v;
double **f = atom->f;
for(int i=0; i<n; i++)
{
double 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];
}
t_sys = 0.0;
if(method==CMD && universe->iworld==0) return;
double expfac;
int nmax = 3 * atom->nlocal;
double KT = force->boltz * nhc_temp;
double kecurrent, t_current;
double dthalf = 0.5 * update->dt;
double dt4 = 0.25 * update->dt;
double dt8 = 0.125 * update->dt;
for(int i=0; i<nmax; i++)
{
int iatm = i/3;
int idim = i%3;
double *vv = v[iatm];
kecurrent = mass[type[iatm]] * vv[idim]* vv[idim] * force->mvv2e;
t_current = kecurrent / force->boltz;
double *eta = nhc_eta[i];
double *eta_dot = nhc_eta_dot[i];
double *eta_dotdot = nhc_eta_dotdot[i];
eta_dotdot[0] = (kecurrent - KT) / nhc_eta_mass[i][0];
for(int ichain=nhc_nchain-1; ichain>0; ichain--)
{
expfac = exp(-dt8 * eta_dot[ichain+1]);
eta_dot[ichain] *= expfac;
eta_dot[ichain] += eta_dotdot[ichain] * dt4;
eta_dot[ichain] *= expfac;
}
expfac = exp(-dt8 * eta_dot[1]);
eta_dot[0] *= expfac;
eta_dot[0] += eta_dotdot[0] * dt4;
eta_dot[0] *= expfac;
// Update particle velocities half-step
double factor_eta = exp(-dthalf * eta_dot[0]);
vv[idim] *= factor_eta;
t_current *= (factor_eta * factor_eta);
kecurrent = force->boltz * t_current;
eta_dotdot[0] = (kecurrent - KT) / nhc_eta_mass[i][0];
for(int ichain=0; ichain<nhc_nchain; ichain++)
eta[ichain] += dthalf * eta_dot[ichain];
eta_dot[0] *= expfac;
eta_dot[0] += eta_dotdot[0] * dt4;
eta_dot[0] *= expfac;
for(int ichain=1; ichain<nhc_nchain; ichain++)
{
expfac = exp(-dt8 * eta_dot[ichain+1]);
eta_dot[ichain] *= expfac;
eta_dotdot[ichain] = (nhc_eta_mass[i][ichain-1] * eta_dot[ichain-1] * eta_dot[ichain-1]
- KT) / nhc_eta_mass[i][ichain];
eta_dot[ichain] += eta_dotdot[ichain] * dt4;
eta_dot[ichain] *= expfac;
}
t_sys += t_current;
}
t_sys /= nmax;
}
/* ----------------------------------------------------------------------
NM PIMD Operations
------------------------------------------------------------------------- */
void FixPIMD::nmpimd_init()
{
memory->create(M_x2xp, np, np, "fix_feynman:M_x2xp");
memory->create(M_xp2x, np, np, "fix_feynman:M_xp2x");
memory->create(M_f2fp, np, np, "fix_feynman:M_f2fp");
memory->create(M_fp2f, np, np, "fix_feynman:M_fp2f");
lam = (double*) memory->smalloc(sizeof(double)*np, "FixPIMD::lam");
// Set up eigenvalues
lam[0] = 0.0;
if(np%2==0) lam[np-1] = 4.0 * np;
for(int i=2; i<=np/2; i++)
{
lam[2*i-3] = lam[2*i-2] = 2.0 * np * (1.0 - 1.0 *cos(2.0*M_PI*(i-1)/np));
}
// Set up eigenvectors for non-degenerated modes
for(int i=0; i<np; i++)
{
M_x2xp[0][i] = 1.0 / np;
if(np%2==0) M_x2xp[np-1][i] = 1.0 / np * pow(-1.0, i);
}
// Set up eigenvectors for degenerated modes
for(int i=0; i<(np-1)/2; i++) for(int j=0; j<np; j++)
{
M_x2xp[2*i+1][j] = sqrt(2.0) * cos ( 2.0 * M_PI * (i+1) * j / np) / np;
M_x2xp[2*i+2][j] = - sqrt(2.0) * sin ( 2.0 * M_PI * (i+1) * j / np) / np;
}
// Set up Ut
for(int i=0; i<np; i++)
for(int j=0; j<np; j++)
{
M_xp2x[i][j] = M_x2xp[j][i] * np;
M_f2fp[i][j] = M_x2xp[i][j] * np;
M_fp2f[i][j] = M_xp2x[i][j];
}
// Set up masses
int iworld = universe->iworld;
for(int i=1; i<=atom->ntypes; i++)
{
mass[i] = atom->mass[i];
if(iworld)
{
mass[i] *= lam[iworld];
mass[i] *= fmass;
}
}
}
/* ---------------------------------------------------------------------- */
void FixPIMD::nmpimd_fill(double **ptr)
{
comm_ptr = ptr;
comm->forward_comm_fix(this);
}
/* ---------------------------------------------------------------------- */
void FixPIMD::nmpimd_transform(double** src, double** des, double *vector)
{
int n = atom->nlocal;
int m = 0;
for(int i=0; i<n; i++) for(int d=0; d<3; d++)
{
des[i][d] = 0.0;
for(int j=0; j<np; j++) { des[i][d] += (src[j][m] * vector[j]); }
m++;
}
}
/* ---------------------------------------------------------------------- */
void FixPIMD::spring_force()
{
spring_energy = 0.0;
double **x = atom->x;
double **f = atom->f;
double* _mass = atom->mass;
int* type = atom->type;
int nlocal = atom->nlocal;
double* xlast = buf_beads[x_last];
double* xnext = buf_beads[x_next];
for(int i=0; i<nlocal; i++)
{
double delx1 = xlast[0] - x[i][0];
double dely1 = xlast[1] - x[i][1];
double delz1 = xlast[2] - x[i][2];
xlast += 3;
domain->minimum_image(delx1, dely1, delz1);
double delx2 = xnext[0] - x[i][0];
double dely2 = xnext[1] - x[i][1];
double delz2 = xnext[2] - x[i][2];
xnext += 3;
domain->minimum_image(delx2, dely2, delz2);
double ff = fbond * _mass[type[i]];
double dx = delx1+delx2;
double dy = dely1+dely2;
double dz = delz1+delz2;
f[i][0] -= (dx) * ff;
f[i][1] -= (dy) * ff;
f[i][2] -= (dz) * ff;
spring_energy += (dx*dx+dy*dy+dz*dz);
}
}
/* ----------------------------------------------------------------------
Comm Operations
------------------------------------------------------------------------- */
void FixPIMD::comm_init()
{
if(size_plan)
{
delete [] plan_send;
delete [] plan_recv;
}
if(method == PIMD)
{
size_plan = 2;
plan_send = new int [2];
plan_recv = new int [2];
mode_index = new int [2];
int rank_last = universe->me - comm->nprocs;
int rank_next = universe->me + comm->nprocs;
if(rank_last<0) rank_last += universe->nprocs;
if(rank_next>=universe->nprocs) rank_next -= universe->nprocs;
plan_send[0] = rank_next; plan_send[1] = rank_last;
plan_recv[0] = rank_last; plan_recv[1] = rank_next;
mode_index[0] = 0; mode_index[1] = 1;
x_last = 1; x_next = 0;
}
else
{
size_plan = np - 1;
plan_send = new int [size_plan];
plan_recv = new int [size_plan];
mode_index = new int [size_plan];
for(int i=0; i<size_plan; i++)
{
plan_send[i] = universe->me + comm->nprocs * (i+1);
if(plan_send[i]>=universe->nprocs) plan_send[i] -= universe->nprocs;
plan_recv[i] = universe->me - comm->nprocs * (i+1);
if(plan_recv[i]<0) plan_recv[i] += universe->nprocs;
mode_index[i]=(universe->iworld+i+1)%(universe->nworlds);
}
x_next = (universe->iworld+1+universe->nworlds)%(universe->nworlds);
x_last = (universe->iworld-1+universe->nworlds)%(universe->nworlds);
}
if(buf_beads)
{
for(int i=0; i<np; i++) if(buf_beads[i]) delete [] buf_beads[i];
delete [] buf_beads;
}
buf_beads = new double* [np];
for(int i=0; i<np; i++) buf_beads[i] = NULL;
}
/* ---------------------------------------------------------------------- */
void FixPIMD::comm_exec(double **ptr)
{
MPI_Status status;
int nlocal = atom->nlocal;
if(nlocal > max_nlocal)
{
max_nlocal = nlocal+200;
int size = sizeof(double) * max_nlocal * 3;
buf_recv = (double*) memory->srealloc(buf_recv, size, "FixPIMD:x_recv");
for(int i=0; i<np; i++)
buf_beads[i] = (double*) memory->srealloc(buf_beads[i], size, "FixPIMD:x_beads[i]");
}
// copy local positions
memcpy(buf_beads[universe->iworld], &(ptr[0][0]), sizeof(double)*nlocal*3);
// go over comm plans
for(int iplan = 0; iplan<size_plan; iplan++)
{
// sendrecv nlocal
int nsend;
MPI_Sendrecv( &(nlocal), 1, MPI_INT, plan_send[iplan], 0,
&(nsend), 1, MPI_INT, plan_recv[iplan], 0, universe->uworld, &status);
// allocate arrays
if(nsend > max_nsend)
{
max_nsend = nsend+200;
tag_send = (int*) memory->srealloc(tag_send, sizeof(int)*max_nsend, "FixPIMD:tag_send");
buf_send = (double*) memory->srealloc(buf_send, sizeof(double)*max_nsend*3, "FixPIMD:x_send");
}
// send tags
MPI_Sendrecv( atom->tag, nlocal, MPI_INT, plan_send[iplan], 0,
tag_send, nsend, MPI_INT, plan_recv[iplan], 0, universe->uworld, &status);
// wrap positions
double *wrap_ptr = buf_send;
int ncpy = sizeof(double)*3;
for(int i=0; i<nsend; i++)
{
int index = atom->map(tag_send[i]);
#define ERR_HEAD "FEYNMAN|"
#define ERR_LAST "Aborted in \"fix pimd\""
#define _ECHO if(universe->me==0 && screen) fprintf(screen,
#define _END );
if(index<0)
{
_ECHO "%s Error: Atom %d is missing at world [%d] rank [%d] required by rank [%d].\n",
ERR_HEAD, tag_send[i], universe->iworld, comm->me, plan_recv[iplan] _END
/* --------------------------------------------------------------------------------------------------- */
char fname[100];
sprintf(fname,"err_%d_%d.xyz", universe->iworld, comm->me);
FILE *fp = fopen(fname,"w");
int n = atom->nlocal + atom->nghost;
fprintf(fp,"%d\n\n", n);
for(int i=0; i<atom->nlocal; i++) fprintf(fp, "O %lf %lf %lf %d\n", atom->x[i][0], atom->x[i][1], atom->x[i][2],atom->tag[i]);
for(int i=nlocal; i<n; i++) fprintf(fp, "H %lf %lf %lf %d\n", atom->x[i][0], atom->x[i][1], atom->x[i][2],atom->tag[i]);
fclose(fp);
/* --------------------------------------------------------------------------------------------------- */
error->universe_one(FLERR,ERR_LAST);
}
memcpy(wrap_ptr, ptr[index], ncpy);
wrap_ptr += 3;
}
// sendrecv x
MPI_Sendrecv( buf_send, nsend*3, MPI_DOUBLE, plan_recv[iplan], 0,
buf_recv, nlocal*3, MPI_DOUBLE, plan_send[iplan], 0, universe->uworld, &status);
// copy x
memcpy(buf_beads[mode_index[iplan]], buf_recv, sizeof(double)*nlocal*3);
}
}
/* ---------------------------------------------------------------------- */
int FixPIMD::pack_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = comm_ptr[j][0];
buf[m++] = comm_ptr[j][1];
buf[m++] = comm_ptr[j][2];
}
return 3;
}
/* ---------------------------------------------------------------------- */
void FixPIMD::unpack_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
comm_ptr[i][0] = buf[m++];
comm_ptr[i][1] = buf[m++];
comm_ptr[i][2] = buf[m++];
}
}
/* ----------------------------------------------------------------------
Memory Operations
------------------------------------------------------------------------- */
double FixPIMD::memory_usage()
{
double bytes = 0;
bytes = atom->nmax * size_peratom_cols * sizeof(double);
return bytes;
}
/* ---------------------------------------------------------------------- */
void FixPIMD::grow_arrays(int nmax)
{
if (nmax==0) return;
int count = nmax*3;
memory->grow(array_atom, nmax, size_peratom_cols, "FixPIMD::array_atom");
memory->grow(nhc_eta, count, nhc_nchain, "FixPIMD::nh_eta");
memory->grow(nhc_eta_dot, count, nhc_nchain+1, "FixPIMD::nh_eta_dot");
memory->grow(nhc_eta_dotdot, count, nhc_nchain, "FixPIMD::nh_eta_dotdot");
memory->grow(nhc_eta_mass, count, nhc_nchain, "FixPIMD::nh_eta_mass");
}
/* ---------------------------------------------------------------------- */
void FixPIMD::copy_arrays(int i, int j, int delflag)
{
int i_pos = i*3;
int j_pos = j*3;
memcpy(nhc_eta [j_pos], nhc_eta [i_pos], nhc_size_one_1);
memcpy(nhc_eta_dot [j_pos], nhc_eta_dot [i_pos], nhc_size_one_2);
memcpy(nhc_eta_dotdot[j_pos], nhc_eta_dotdot[i_pos], nhc_size_one_1);
memcpy(nhc_eta_mass [j_pos], nhc_eta_mass [i_pos], nhc_size_one_1);
}
/* ---------------------------------------------------------------------- */
int FixPIMD::pack_exchange(int i, double *buf)
{
int offset=0;
int pos = i * 3;
memcpy(buf+offset, nhc_eta[pos], nhc_size_one_1); offset += nhc_offset_one_1;
memcpy(buf+offset, nhc_eta_dot[pos], nhc_size_one_2); offset += nhc_offset_one_2;
memcpy(buf+offset, nhc_eta_dotdot[pos], nhc_size_one_1); offset += nhc_offset_one_1;
memcpy(buf+offset, nhc_eta_mass[pos], nhc_size_one_1); offset += nhc_offset_one_1;
return size_peratom_cols;
}
/* ---------------------------------------------------------------------- */
int FixPIMD::unpack_exchange(int nlocal, double *buf)
{
int offset=0;
int pos = nlocal*3;
memcpy(nhc_eta[pos], buf+offset, nhc_size_one_1); offset += nhc_offset_one_1;
memcpy(nhc_eta_dot[pos], buf+offset, nhc_size_one_2); offset += nhc_offset_one_2;
memcpy(nhc_eta_dotdot[pos], buf+offset, nhc_size_one_1); offset += nhc_offset_one_1;
memcpy(nhc_eta_mass[pos], buf+offset, nhc_size_one_1); offset += nhc_offset_one_1;
return size_peratom_cols;
}
/* ---------------------------------------------------------------------- */
int FixPIMD::pack_restart(int i, double *buf)
{
int offset=0;
int pos = i * 3;
buf[offset++] = size_peratom_cols+1;
memcpy(buf+offset, nhc_eta[pos], nhc_size_one_1); offset += nhc_offset_one_1;
memcpy(buf+offset, nhc_eta_dot[pos], nhc_size_one_2); offset += nhc_offset_one_2;
memcpy(buf+offset, nhc_eta_dotdot[pos], nhc_size_one_1); offset += nhc_offset_one_1;
memcpy(buf+offset, nhc_eta_mass[pos], nhc_size_one_1); offset += nhc_offset_one_1;
return size_peratom_cols+1;
}
/* ---------------------------------------------------------------------- */
void FixPIMD::unpack_restart(int nlocal, int nth)
{
double **extra = atom->extra;
// skip to Nth set of extra values
int m = 0;
for (int i=0; i<nth; i++) m += static_cast<int> (extra[nlocal][m]);
m++;
int pos = nlocal * 3;
memcpy(nhc_eta[pos], extra[nlocal]+m, nhc_size_one_1); m += nhc_offset_one_1;
memcpy(nhc_eta_dot[pos], extra[nlocal]+m, nhc_size_one_2); m += nhc_offset_one_2;
memcpy(nhc_eta_dotdot[pos], extra[nlocal]+m, nhc_size_one_1); m += nhc_offset_one_1;
memcpy(nhc_eta_mass[pos], extra[nlocal]+m, nhc_size_one_1); m += nhc_offset_one_1;
nhc_ready = true;
}
/* ---------------------------------------------------------------------- */
int FixPIMD::maxsize_restart()
{
return size_peratom_cols+1;
}
/* ---------------------------------------------------------------------- */
int FixPIMD::size_restart(int nlocal)
{
return size_peratom_cols+1;
}
/* ---------------------------------------------------------------------- */
double FixPIMD::compute_vector(int n)
{
if(n==0) { return spring_energy; }
if(n==1) { return t_sys; }
return 0.0;
}

115
src/USER-MISC/fix_pimd.h Normal file
View File

@ -0,0 +1,115 @@
/* ----------------------------------------------------------------------
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(pimd,FixPIMD)
#else
#ifndef FIX_PIMD_H
#define FIX_PIMD_H
#include "fix.h"
namespace LAMMPS_NS {
class FixPIMD : public Fix {
public:
FixPIMD(class LAMMPS *, int, char **);
int setmask();
void init();
void setup(int);
void post_force(int);
void initial_integrate(int);
void final_integrate();
double memory_usage();
void grow_arrays(int);
void copy_arrays(int,int,int);
int pack_exchange(int,double*);
int unpack_exchange(int,double*);
int pack_restart(int,double*);
void unpack_restart(int,int);
int maxsize_restart();
int size_restart(int);
double compute_vector(int);
int pack_comm(int, int*, double *, int, int*);
void unpack_comm(int, int, double *);
int method;
int np;
double inverse_np;
/* ring-polymer model */
double omega_np, fbond, spring_energy, sp;
int x_last, x_next;
void spring_force();
/* fictious mass */
double fmass, *mass;
/* inter-partition communication */
int max_nsend;
int* tag_send;
double *buf_send;
int max_nlocal;
double *buf_recv, **buf_beads;
int size_plan;
int *plan_send, *plan_recv;
double **comm_ptr;
void comm_init();
void comm_exec(double **);
/* normal-mode operations */
double *lam, **M_x2xp, **M_xp2x, **M_f2fp, **M_fp2f;
int *mode_index;
void nmpimd_init();
void nmpimd_fill(double**);
void nmpimd_transform(double**, double**, double*);
/* Nose-hoover chain integration */
int nhc_offset_one_1, nhc_offset_one_2;
int nhc_size_one_1, nhc_size_one_2;
int nhc_nchain;
bool nhc_ready;
double nhc_temp, dtv, dtf, t_sys;
double **nhc_eta; /* coordinates of NH chains for ring-polymer beads */
double **nhc_eta_dot; /* velocities of NH chains */
double **nhc_eta_dotdot; /* acceleration of NH chains */
double **nhc_eta_mass; /* mass of NH chains */
void nhc_init();
void nhc_update_v();
void nhc_update_x();
};
}
#endif
#endif