Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F68827259
fix_pimd.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Jun 29, 02:33
Size
22 KB
Mime Type
text/x-c
Expires
Mon, Jul 1, 02:33 (2 d)
Engine
blob
Format
Raw Data
Handle
18323497
Attached To
rLAMMPS lammps
fix_pimd.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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Package FixPIMD
Purpose Quantum Path Integral Algorithm for Quantum Chemistry
Copyright Voth Group @ University of Chicago
Authors Chris Knight & Yuxing Peng (yuxing at uchicago.edu)
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 (atom->map_style == 0)
error->all(FLERR,"Fix pimd requires an atom map, see atom_modify");
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;
if((method==CMD || method==NMPIMD) && universe->iworld==0) ; else 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;
}
/* ----------------------------------------------------------------------
Normal Mode PIMD
------------------------------------------------------------------------- */
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)
{
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, MPI_STATUS_IGNORE);
// 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, MPI_STATUS_IGNORE);
// 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]);
if(index<0)
{
char error_line[256];
sprintf(error_line, "Atom " TAGINT_FORMAT " is missing at world [%d] "
"rank [%d] required by rank [%d] (" TAGINT_FORMAT ", "
TAGINT_FORMAT ", " TAGINT_FORMAT ").\n",tag_send[i],
universe->iworld, comm->me, plan_recv[iplan],
atom->tag[0], atom->tag[1], atom->tag[2]);
error->universe_one(FLERR,error_line);
}
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, MPI_STATUS_IGNORE);
// copy x
memcpy(buf_beads[mode_index[iplan]], buf_recv, sizeof(double)*nlocal*3);
}
}
/* ---------------------------------------------------------------------- */
int FixPIMD::pack_forward_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 m;
}
/* ---------------------------------------------------------------------- */
void FixPIMD::unpack_forward_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;
}
Event Timeline
Log In to Comment