diff --git a/examples/vashishta/in.indiumphosphide-tabulated b/examples/vashishta/in.indiumphosphide-tabulated new file mode 100644 index 0000000000..3224b9fb75 --- /dev/null +++ b/examples/vashishta/in.indiumphosphide-tabulated @@ -0,0 +1,76 @@ +# calculate the energy volume curve for InP zincblende + +# define volume range and filename + +variable ndelta equal 100 +variable volatom_min equal 20.0 +variable volatom_max equal 29.0 +variable evsvolfile string evsvol.dat + +# set up cell + +units metal + +boundary p p p + +# setup loop variables for box volume + +variable amin equal ${volatom_min}^(1/3)*2 +variable delta equal (${volatom_max}-${volatom_min})/${ndelta} +variable scale equal (${delta}/v_volatom+1)^(1/3) + +# set up 8 atom InP zincblende unit cell + +lattice diamond ${amin} + +region box prism & + 0 1 & + 0 1 & + 0 1 & + 0 0 0 + +create_box 2 box + +create_atoms 1 box & + basis 5 2 & + basis 6 2 & + basis 7 2 & + basis 8 2 + +mass 1 114.76 +mass 2 30.98 + +# choose potential + +pair_style vashishta +pair_coeff * * InP.vashishta In P +pair_modify table 16 + +# setup neighbor style + +neighbor 1.0 nsq +neigh_modify once no every 1 delay 0 check yes + +# setup output + +thermo_style custom step temp pe press vol +thermo_modify norm no +variable volatom equal vol/atoms +variable eatom equal pe/atoms +print "# Volume [A^3/atom] Energy [eV/atom]" file ${evsvolfile} + +# loop over range of volumes + +label loop +variable i loop ${ndelta} + +change_box all x scale ${scale} y scale ${scale} z scale ${scale} remap + +# calculate energy +# no energy minimization needed for zincblende + +run 0 +print "${volatom} ${eatom}" append ${evsvolfile} + +next i +jump SELF loop diff --git a/examples/vashishta/in.sio2-tabulated b/examples/vashishta/in.sio2-tabulated new file mode 100644 index 0000000000..dc1381af2b --- /dev/null +++ b/examples/vashishta/in.sio2-tabulated @@ -0,0 +1,29 @@ +# test Vashishta potential for quartz + +units metal +boundary p p p + +atom_style atomic + +read_data data.quartz + +replicate 4 4 4 +velocity all create 2000.0 277387 mom yes +displace_atoms all move 0.05 0.9 0.4 units box + +pair_style vashishta +pair_coeff * * SiO.1990.vashishta Si O +pair_modify table 16 + +neighbor 0.3 bin +neigh_modify delay 10 + +fix 1 all nve +thermo 10 +timestep 0.001 + +#dump 1 all cfg 10 *.cfg mass type xs ys zs vx vy vz fx fy fz +#dump_modify 1 element Si O + +run 100 + diff --git a/src/MANYBODY/pair_vashishta.cpp b/src/MANYBODY/pair_vashishta.cpp index 0a0ea75d2a..a451e28aea 100755 --- a/src/MANYBODY/pair_vashishta.cpp +++ b/src/MANYBODY/pair_vashishta.cpp @@ -97,11 +97,10 @@ void PairVashishta::modify_params(int narg, char **arg) if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); nTablebits = force->inumeric(FLERR,arg[iarg+1]); - if (nTablebits > sizeof(float)*CHAR_BIT) { - error->all(FLERR,"Too many total bits for bitmapped lookup table"); - } + if ((nTablebits < 0) || (nTablebits > sizeof(int)*CHAR_BIT-1)) + error->all(FLERR,"Illegal interpolation table size"); - if(nTablebits == 0) { + if (nTablebits == 0) { useTable = false; } else { useTable = true; @@ -112,11 +111,12 @@ void PairVashishta::modify_params(int narg, char **arg) if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); tabinner = force->numeric(FLERR,arg[iarg+1]); + if (tabinner <= 0.0) + error->all(FLERR,"Illegal inner cutoff for tabulation"); iarg += 2; } } - - createTable(); + updateTables(); } void PairVashishta::validateNeigh3Body() { @@ -129,40 +129,34 @@ void PairVashishta::validateNeigh3Body() { } } -void PairVashishta::createTable() +void PairVashishta::updateTables() { - int ntypes = atom->ntypes+1; + memory->destroy(forceTable); + memory->destroy(potentialTable); + forceTable = NULL; + potentialTable = NULL; + + if (!useTable) return; + + int ntable = 1<destroy(forceTable); - memory->destroy(potentialTable); - memory->create(forceTable,ntypes,ntypes,ntable+1,"pair:vashishta:forceTable"); - memory->create(potentialTable,ntypes,ntypes,ntable+1,"pair:vashishta:potentialTable"); + memory->create(forceTable,nelements,nelements,ntable+1,"pair:vashishta:forceTable"); + memory->create(potentialTable,nelements,nelements,ntable+1,"pair:vashishta:potentialTable"); - for (int ii = 0; ii < ntypes; ii++) { - int i = map[ii]; - for (int jj = 0; jj < ntypes; jj++) { - int j = map[jj]; - if (i < 0 || j < 0 || ii == 0 || jj == 0) { - for(int tableIndex=0; tableIndex<=ntable; tableIndex++) { - forceTable[ii][jj][tableIndex] = 0; - potentialTable[ii][jj][tableIndex] = 0; - } - } else { - int ijparam = elem2param[i][j][j]; - for(int tableIndex=0; tableIndex<=ntable; tableIndex++) { - double rsq = tabinnersq + tableIndex*deltaR2; - double fpair, eng; - twobody(¶ms[ijparam], rsq, fpair, 1, eng, false /* Don't ask for tabulated since we are generating it now*/ ); - forceTable[ii][jj][tableIndex] = fpair; - potentialTable[ii][jj][tableIndex] = eng; - } + for (int i = 0; i < nelements; i++) { + for (int j = 0; j < nelements; j++) { + int ijparam = elem2param[i][j][j]; + for(int idx=0; idx <= ntable; idx++) { + double rsq = tabinnersq + idx*deltaR2; + double fpair, eng; + // call analytical version of two-body term function + twobody(¶ms[ijparam], rsq, fpair, 1, eng, false); + forceTable[i][j][idx] = fpair; + potentialTable[i][j][idx] = eng; } } } @@ -642,7 +636,7 @@ void PairVashishta::setup_params() if (params[m].r0 > cutmax) cutmax = params[m].r0; } - createTable(); + updateTables(); } /* ---------------------------------------------------------------------- */ @@ -660,12 +654,12 @@ void PairVashishta::twobody(Param *param, double rsq, double &fforce, int tableIndex = (rsq - tabinnersq)*oneOverDeltaR2; double fraction = (rsq - tabinnersq)*oneOverDeltaR2 - tableIndex; // double - int will only keep the 0.xxxx part - double force0 = forceTable[param->ielement+1][param->jelement+1][tableIndex]; - double force1 = forceTable[param->ielement+1][param->jelement+1][tableIndex+1]; + double force0 = forceTable[param->ielement][param->jelement][tableIndex]; + double force1 = forceTable[param->ielement][param->jelement][tableIndex+1]; fforce = (1.0 - fraction)*force0 + fraction*force1; // force is linearly interpolated between the two values if(evflag) { - double energy0 = potentialTable[param->ielement+1][param->jelement+1][tableIndex]; - double energy1 = potentialTable[param->ielement+1][param->jelement+1][tableIndex+1]; + double energy0 = potentialTable[param->ielement][param->jelement][tableIndex]; + double energy1 = potentialTable[param->ielement][param->jelement][tableIndex+1]; eng = (1.0 - fraction)*energy0 + fraction*energy1; } } else { @@ -747,3 +741,15 @@ void PairVashishta::threebody(Param *paramij, Param *paramik, Param *paramijk, if (eflag) eng = facrad; } + +/* ---------------------------------------------------------------------- + * add memory usage for tabulation + * ---------------------------------------------------------------------- */ + +double PairVashishta::memory_usage() +{ + double bytes = Pair::memory_usage(); + if (useTable) + bytes += 2*nelements*nelements*sizeof(double)*(1<