diff --git a/src/DSMC/Install.sh b/src/DSMC/Install.sh deleted file mode 100644 index 60e34982ca..0000000000 --- a/src/DSMC/Install.sh +++ /dev/null @@ -1,15 +0,0 @@ -# Install/unInstall package files in LAMMPS - -if (test $1 = 1) then - - cp pair_dsmc.cpp .. - - cp pair_dsmc.h .. - -elif (test $1 = 0) then - - rm -f ../pair_dsmc.cpp - - rm -f ../pair_dsmc.h - -fi diff --git a/src/DSMC/pair_dsmc.cpp b/src/DSMC/pair_dsmc.cpp deleted file mode 100644 index b04dd8f6d7..0000000000 --- a/src/DSMC/pair_dsmc.cpp +++ /dev/null @@ -1,525 +0,0 @@ -/* ---------------------------------------------------------------------- - 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: Paul Crozier (SNL) -------------------------------------------------------------------------- */ - -#include "math.h" -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "pair_dsmc.h" -#include "atom.h" -#include "comm.h" -#include "force.h" -#include "memory.h" -#include "error.h" -#include "domain.h" -#include "update.h" -#include "random_mars.h" -#include "limits.h" - -using namespace LAMMPS_NS; - -#define MIN(a,b) ((a) < (b) ? (a) : (b)) -#define MAX(a,b) ((a) > (b) ? (a) : (b)) - -/* ---------------------------------------------------------------------- */ - -PairDSMC::PairDSMC(LAMMPS *lmp) : Pair(lmp) -{ - single_enable = 0; - - total_number_of_collisions = 0; - max_particles = max_particle_list = 0; - next_particle = NULL; - random = NULL; -} - -/* ---------------------------------------------------------------------- */ - -PairDSMC::~PairDSMC() -{ - if (allocated) { - memory->destroy(setflag); - memory->destroy(cutsq); - memory->destroy(sigma); - memory->destroy(cut); - memory->destroy(V_sigma_max); - memory->destroy(particle_list); - memory->destroy(first); - memory->destroy(number); - } - - delete [] next_particle; - delete random; -} - -/* ---------------------------------------------------------------------- */ - -void PairDSMC::compute(int eflag, int vflag) -{ - double **x = atom->x; - double *mass = atom->mass; - int *type = atom->type; - int nlocal = atom->nlocal; - - for (int i = 1; i <= atom->ntypes; ++i) - for (int j = 0; j < total_ncells; ++j) { - first[i][j] = -1; - number[i][j] = 0; - } - - if (atom->nmax > max_particles) { - delete [] next_particle; - max_particles = atom->nmax; - next_particle = new int[max_particles]; - } - - // find each particle's cell and sort by type - // assume a constant volume and shape simulation domain - // skip particle if outside processor domain - - for (int i = 0; i < nlocal; ++i) { - int xcell = static_cast((x[i][0] - domain->boxlo[0])/cellx); - int ycell = static_cast((x[i][1] - domain->boxlo[1])/celly); - int zcell = static_cast((x[i][2] - domain->boxlo[2])/cellz); - - if ((xcell < 0) or (xcell > ncellsx-1) or - (ycell < 0) or (ycell > ncellsy-1) or - (zcell < 0) or (zcell > ncellsz-1)) continue; - - int icell = xcell + ycell*ncellsx + zcell*ncellsx*ncellsy; - itype = type[i]; - next_particle[i] = first[itype][icell]; - first[itype][icell] = i; - number[itype][icell]++; - } - - for (int icell = 0; icell < total_ncells; ++icell) { - - for (itype = 1; itype <= atom->ntypes; ++itype) { - number_of_A = number[itype][icell]; - if (number_of_A > max_particle_list) { - max_particle_list = number_of_A; - memory->grow(particle_list,atom->ntypes+1,max_particle_list, - "pair:particle_list"); - } - - int m = first[itype][icell]; - for (int k = 0; k < number_of_A; k++) { - particle_list[itype][k] = m; - m = next_particle[m]; - } - } - - for (itype = 1; itype <= atom->ntypes; ++itype) { - imass = mass[itype]; - number_of_A = number[itype][icell]; - - for (jtype = itype; jtype <= atom->ntypes; ++jtype) { - jmass = mass[jtype]; - number_of_B = number[jtype][icell]; - - reduced_mass = imass*jmass/(imass + jmass); - total_mass = imass + jmass; - jmass_tmass = jmass/total_mass; - imass_tmass = imass/total_mass; - - // if necessary, recompute V_sigma_max values - - if (recompute_vsigmamax_stride && - (update->ntimestep % recompute_vsigmamax_stride == 0)) - recompute_V_sigma_max(icell); - - // # of collisions to perform for itype-jtype pairs - - double &Vs_max = V_sigma_max[itype][jtype]; - double num_of_collisions_double = number_of_A * number_of_B * - weighting * Vs_max * update->dt / vol; - - if ((itype == jtype) and number_of_B) - num_of_collisions_double *= - 0.5 * double(number_of_B - 1) / double(number_of_B); - - int num_of_collisions = - convert_double_to_equivalent_int(num_of_collisions_double); - - if (num_of_collisions > number_of_A) - error->warning("Pair dsmc: num_of_collisions > number_of_A",0); - if (num_of_collisions > number_of_B) - error->warning("Pair dsmc: num_of_collisions > number_of_B",0); - - // perform collisions on pairs of particles in icell - - for (int k = 0; k < num_of_collisions; k++) { - if ((number_of_A < 1) or (number_of_B < 1)) break; - int ith_A = static_cast(random->uniform()*number_of_A); - int jth_B = static_cast(random->uniform()*number_of_B); - int i = particle_list[itype][ith_A]; - int j = particle_list[jtype][jth_B]; - if (i == j) { - k--; - continue; - } - double probability = V_sigma(i,j)/Vs_max; - if (probability > random->uniform()) scatter_random(i,j,icell); - } - } - } - } -} - -/* ---------------------------------------------------------------------- - allocate all arrays -------------------------------------------------------------------------- */ - -void PairDSMC::allocate() -{ - allocated = 1; - int n = atom->ntypes; - - memory->create(setflag,n+1,n+1,"pair:setflag"); - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; - - memory->create(cutsq,n+1,n+1,"pair:cutsq"); - - memory->create(cut,n+1,n+1,"pair:cut"); - memory->create(sigma,n+1,n+1,"pair:sigma"); - memory->create(V_sigma_max,n+1,n+1,"pair:V_sigma_max"); -} - -/* ---------------------------------------------------------------------- - global settings -------------------------------------------------------------------------- */ - -void PairDSMC::settings(int narg, char **arg) -{ - if (narg != 6) error->all("Illegal pair_style command"); - - cut_global = 0.0; - max_cell_size = force->numeric(arg[0]); - seed = force->inumeric(arg[1]); - weighting = force->numeric(arg[2]); - T_ref = force->numeric(arg[3]); - recompute_vsigmamax_stride = force->inumeric(arg[4]); - vsigmamax_samples = force->inumeric(arg[5]); - - // initialize Marsaglia RNG with processor-unique seed - - if (max_cell_size <= 0.0) error->all("Illegal pair_style command"); - if (seed <= 0) error->all("Illegal pair_style command"); - if (random) delete random; - random = new RanMars(lmp,seed + comm->me); - - kT_ref = force->boltz*T_ref; - - // 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 PairDSMC::coeff(int narg, char **arg) -{ - if (narg < 3 || narg > 4) 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 sigma_one = force->numeric(arg[2]); - - double cut_one = cut_global; - if (narg == 4) cut_one = force->numeric(arg[3]); - - int count = 0; - for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { - sigma[i][j] = sigma_one; - 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 PairDSMC::init_style() -{ - ncellsx = ncellsy = ncellsz = 1; - while (((domain->boxhi[0] - domain->boxlo[0])/ncellsx) > max_cell_size) - ncellsx++; - while (((domain->boxhi[1] - domain->boxlo[1])/ncellsy) > max_cell_size) - ncellsy++; - while (((domain->boxhi[2] - domain->boxlo[2])/ncellsz) > max_cell_size) - ncellsz++; - - cellx = (domain->boxhi[0] - domain->boxlo[0])/ncellsx; - celly = (domain->boxhi[1] - domain->boxlo[1])/ncellsy; - cellz = (domain->boxhi[2] - domain->boxlo[2])/ncellsz; - - if (comm->me == 0) { - if (screen) fprintf(screen,"DSMC cell size = %g x %g x %g\n", - cellx,celly,cellz); - if (logfile) fprintf(logfile,"DSMC cell size = %g x %g x %g\n", - cellx,celly,cellz); - } - - total_ncells = ncellsx*ncellsy*ncellsz; - vol = cellx*celly*cellz; - - memory->create(particle_list,atom->ntypes+1,0,"pair:particle_list"); - memory->create(first,atom->ntypes+1,total_ncells,"pair:first"); - memory->create(number,atom->ntypes+1,total_ncells,"pair:number"); - - for (int i = 1; i <= atom->ntypes; i++) - for (int j = 1; j <= atom->ntypes; j++) - V_sigma_max[i][j] = 0.0; - - two_pi = 8.0*atan(1.0); -} - -/* ---------------------------------------------------------------------- - init for one type pair i,j and corresponding j,i -------------------------------------------------------------------------- */ - -double PairDSMC::init_one(int i, int j) -{ - if (setflag[i][j] == 0) cut[i][j] = 0.0; - return cut[i][j]; -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairDSMC::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(&sigma[i][j],sizeof(double),1,fp); - fwrite(&cut[i][j],sizeof(double),1,fp); - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairDSMC::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(&sigma[i][j],sizeof(double),1,fp); - fread(&cut[i][j],sizeof(double),1,fp); - } - MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairDSMC::write_restart_settings(FILE *fp) -{ - fwrite(&cut_global,sizeof(double),1,fp); - fwrite(&max_cell_size,sizeof(double),1,fp); - fwrite(&seed,sizeof(int),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairDSMC::read_restart_settings(FILE *fp) -{ - if (comm->me == 0) { - fread(&cut_global,sizeof(double),1,fp); - fread(&max_cell_size,sizeof(double),1,fp); - fread(&seed,sizeof(int),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(&max_cell_size,1,MPI_DOUBLE,0,world); - MPI_Bcast(&seed,1,MPI_INT,0,world); - MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); - - // initialize Marsaglia RNG with processor-unique seed - // same seed that pair_style command initially specified - - if (random) delete random; - random = new RanMars(lmp,seed + comm->me); -} - -/*------------------------------------------------------------------------- - rezero and recompute the V_sigma_max values this timestep for use during - the next nrezero timesteps --------------------------------------------------------------------------*/ - -void PairDSMC::recompute_V_sigma_max(int icell) -{ - int i,j,k; - double Vsigma_max = 0; - - if (number_of_A && number_of_B) { - for (k = 0; k < vsigmamax_samples; k++) { - i = particle_list[itype] - [static_cast(random->uniform()*number_of_A)]; - j = particle_list[jtype] - [static_cast(random->uniform()*number_of_B)]; - if (i == j) continue; - Vsigma_max = MAX(Vsigma_max,V_sigma(i,j)); - } - } - V_sigma_max[itype][jtype] = Vsigma_max; -} - -/*------------------------------------------------------------------------- - VHS model - compute the velocity vector difference between i and j and multiply by - their combined collision cross section, sigma, for neutral-neutral - collisions using the Variable Hard Sphere model --------------------------------------------------------------------------*/ - -double PairDSMC::V_sigma(int i, int j) -{ - double relative_velocity_sq,relative_velocity,pair_sigma; - double delv[3]; - double *vi = atom->v[i]; - double *vj = atom->v[j]; - - subtract3d(vi,vj,delv); - relative_velocity_sq = dot3d(delv,delv); - relative_velocity = sqrt(relative_velocity_sq); - - // from Bird eq 4.63, and omega=0.67 - // (omega - 0.5) = 0.17 - // 1/GAMMA(2.5 - omega) = 1.06418029298371 - - if (relative_velocity_sq != 0.0) - pair_sigma = sigma[itype][jtype]* - pow(kT_ref/(0.5*reduced_mass*relative_velocity_sq),0.17) * - 1.06418029298371; - else - pair_sigma = 0.0; - - return relative_velocity*pair_sigma; -} - -/*------------------------------------------------------------------------- - generate new velocities for collided particles --------------------------------------------------------------------------*/ - -void PairDSMC::scatter_random(int i, int j, int icell) -{ - double mag_delv,cos_phi,cos_squared,r,theta; - double delv[3],vcm[3]; - double *vi = atom->v[i]; - double *vj = atom->v[j]; - - subtract3d(vi,vj,delv); - if (itype == jtype) mag_delv = sqrt(dot3d(delv,delv))*0.5; - else mag_delv = sqrt(dot3d(delv,delv)); - - cos_phi = 1.0 - (2.0*random->uniform()); - cos_squared = MIN(1.0,cos_phi*cos_phi); - r = sqrt(1.0 - cos_squared); - delv[0] = cos_phi*mag_delv; - theta = two_pi*random->uniform(); - delv[1] = r*mag_delv*cos(theta); - delv[2] = r*mag_delv*sin(theta); - - if (itype == jtype) { - vcm[0] = (vi[0]+vj[0])*0.5; - vcm[1] = (vi[1]+vj[1])*0.5; - vcm[2] = (vi[2]+vj[2])*0.5; - vi[0] = vcm[0] + delv[0]; - vi[1] = vcm[1] + delv[1]; - vi[2] = vcm[2] + delv[2]; - vj[0] = vcm[0] - delv[0]; - vj[1] = vcm[1] - delv[1]; - vj[2] = vcm[2] - delv[2]; - } else { - vcm[0] = vi[0]*imass_tmass + vj[0]*jmass_tmass; - vcm[1] = vi[1]*imass_tmass + vj[1]*jmass_tmass; - vcm[2] = vi[2]*imass_tmass + vj[2]*jmass_tmass; - vi[0] = vcm[0] + delv[0]*jmass_tmass; - vi[1] = vcm[1] + delv[1]*jmass_tmass; - vi[2] = vcm[2] + delv[2]*jmass_tmass; - vj[0] = vcm[0] - delv[0]*imass_tmass; - vj[1] = vcm[1] - delv[1]*imass_tmass; - vj[2] = vcm[2] - delv[2]*imass_tmass; - } - - total_number_of_collisions++; -} - -/* ---------------------------------------------------------------------- - This method converts the double supplied by the calling function into - an int, which is returned. By adding a random number between 0 and 1 - to the double before converting it to an int, we ensure that, - statistically, we round down with probability identical to the - remainder and up the rest of the time. So even though we're using an - integer, we're statistically matching the exact expression represented - by the double. -------------------------------------------------------------------------- */ - -int PairDSMC::convert_double_to_equivalent_int(double input_double) -{ - if (input_double > INT_MAX) - error->all("Tried to convert a double to int, but input_double > INT_MAX"); - - int output_int = static_cast(input_double + random->uniform()); - return output_int; -} diff --git a/src/DSMC/pair_dsmc.h b/src/DSMC/pair_dsmc.h deleted file mode 100644 index 7e2b784f73..0000000000 --- a/src/DSMC/pair_dsmc.h +++ /dev/null @@ -1,109 +0,0 @@ -/* ---------------------------------------------------------------------- - 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(dsmc,PairDSMC) - -#else - -#ifndef LMP_PAIR_DSMC_H -#define LMP_PAIR_DSMC_H - -#include "pair.h" - -namespace LAMMPS_NS { - -class PairDSMC : public Pair { - public: - PairDSMC(class LAMMPS *); - virtual ~PairDSMC(); - virtual void compute(int, int); - virtual void settings(int, char **); - void coeff(int, char **); - void init_style(); - 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 *); - - private: - double cut_global; - double **cut; - double **sigma; - - double cellx; - double celly; - double cellz; - int ncellsx; - int ncellsy; - int ncellsz; - int total_ncells; - int total_number_of_collisions; - int recompute_vsigmamax_stride; - int vsigmamax_samples; - double T_ref; - double kT_ref; - double two_pi; - double max_cell_size; - - int seed; - int number_of_A; - int number_of_B; - int max_particle_list; - - class RanMars *random; - - int **particle_list; - int **first; - int **number; - - double **V_sigma_max; - - int max_particles; - int *next_particle; - - int itype; - int jtype; - - double imass; - double jmass; - double total_mass; - double reduced_mass; - double imass_tmass; - double jmass_tmass; - double vol; - double weighting; - - void allocate(); - void recompute_V_sigma_max(int); - double V_sigma(int, int); - void scatter_random(int, int, int); - int convert_double_to_equivalent_int(double); - - inline void subtract3d(const double *v1, const double *v2, double *v3) { - v3[0] = v2[0] - v1[0]; - v3[1] = v2[1] - v1[1]; - v3[2] = v2[2] - v1[2]; - } - - inline double dot3d(const double *v1, const double *v2) { - return( v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2] ); - } -}; - -} - -#endif -#endif diff --git a/src/Makefile b/src/Makefile index 3bc766272f..d09c371239 100755 --- a/src/Makefile +++ b/src/Makefile @@ -13,8 +13,8 @@ OBJ = $(SRC:.cpp=.o) # Package variables -PACKAGE = asphere class2 colloid dipole dsmc gpu granular \ - kspace manybody meam molecule opt peri poems reax replica \ +PACKAGE = asphere class2 colloid dipole gpu granular \ + kspace manybody mc meam molecule opt peri poems reax replica \ shock srd xtc PACKUSER = user-misc user-atc user-awpmd user-cg-cmm \ diff --git a/src/fix_restrain.cpp b/src/fix_restrain.cpp index a4e040f0c3..b9e16c7a41 100644 --- a/src/fix_restrain.cpp +++ b/src/fix_restrain.cpp @@ -288,7 +288,7 @@ void FixRestrain::restrain_dihedral() MPI_Comm_rank(world,&me); if (screen) { char str[128]; - sprintf(str,"Restrain problem: %d %d %d %d %d %d", + sprintf(str,"Restrain problem: %d " BIGINT_FORMAT " %d %d %d %d", me,update->ntimestep, atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]); error->warning(str); diff --git a/src/force.h b/src/force.h index df214c974f..fc1bd0a812 100644 --- a/src/force.h +++ b/src/force.h @@ -21,6 +21,7 @@ namespace LAMMPS_NS { class Force : protected Pointers { public: double boltz; // Boltzmann constant (eng/degree-K) + double hplanck; // Planck's constant (energy-time) double mvv2e; // conversion of mv^2 to energy double ftm2v; // conversion of ft/m to velocity double mv2d; // conversion of mass/volume to density diff --git a/src/update.cpp b/src/update.cpp index 3320ce9d8a..947bba8c1e 100644 --- a/src/update.cpp +++ b/src/update.cpp @@ -118,6 +118,7 @@ void Update::set_units(const char *style) if (strcmp(style,"lj") == 0) { force->boltz = 1.0; + force->hplanck = 0.18292026; // using LJ parameters for argon force->mvv2e = 1.0; force->ftm2v = 1.0; force->mv2d = 1.0; @@ -135,6 +136,7 @@ void Update::set_units(const char *style) } else if (strcmp(style,"real") == 0) { force->boltz = 0.0019872067; + force->hplanck = 95.306976368; force->mvv2e = 48.88821291 * 48.88821291; force->ftm2v = 1.0 / 48.88821291 / 48.88821291; force->mv2d = 1.0 / 0.602214179; @@ -152,6 +154,7 @@ void Update::set_units(const char *style) } else if (strcmp(style,"metal") == 0) { force->boltz = 8.617343e-5; + force->hplanck = 4.135667403e-3; force->mvv2e = 1.0364269e-4; force->ftm2v = 1.0 / 1.0364269e-4; force->mv2d = 1.0 / 0.602214179; @@ -169,6 +172,7 @@ void Update::set_units(const char *style) } else if (strcmp(style,"si") == 0) { force->boltz = 1.3806504e-23; + force->hplanck = 6.62606896e-34; force->mvv2e = 1.0; force->ftm2v = 1.0; force->mv2d = 1.0; @@ -186,6 +190,7 @@ void Update::set_units(const char *style) } else if (strcmp(style,"cgs") == 0) { force->boltz = 1.3806504e-16; + force->hplanck = 6.62606896e-27; force->mvv2e = 1.0; force->ftm2v = 1.0; force->mv2d = 1.0; @@ -203,6 +208,7 @@ void Update::set_units(const char *style) } else if (strcmp(style,"electron") == 0) { force->boltz = 3.16681534e-6; + force->hplanck = 0.1519829846; force->mvv2e = 1.06657236; force->ftm2v = 0.937582899; force->mv2d = 1.0;