diff --git a/doc/variable.txt b/doc/variable.txt index 8208c8a90..7fd1a27c5 100644 --- a/doc/variable.txt +++ b/doc/variable.txt @@ -46,7 +46,8 @@ style = {delete} or {index} or {loop} or {world} or {universe} or {uloop} or {st math functions = sqrt(x), exp(x), ln(x), log(x), abs(x), sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x), random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x) - ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z) + ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z), + seed(n) group functions = count(group), mass(group), charge(group), xcm(group,dim), vcm(group,dim), fcm(group,dim), bound(group,xmin), gyration(group), ke(group), @@ -553,6 +554,12 @@ the {start} keyword of the "run"_run.html command. See the "thermo_style"_thermo_style.html keyword elaplong = timestep-startstep. +The seed(n) function generates a prime number between 10000 and +INT_MAX (2147483647). The output of this function can be used to be the seed of +particletemplate, particledistribution, etc. +If the argument is 0, a random prime number is generated, else the n-th greater, +than 10000 prime (n = 1 -> 10007, n = 2 -> 10009, etc) is the output. + :line Group and Region Functions :h4 diff --git a/src/lammps.cpp b/src/lammps.cpp index 1c347ad37..95e8f9c35 100644 --- a/src/lammps.cpp +++ b/src/lammps.cpp @@ -47,6 +47,8 @@ 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. + + Tóth János (MATE, Gödöllő) ------------------------------------------------------------------------- */ #include @@ -86,6 +88,9 @@ #include "error.h" #include "granular_styles.h" +#include +#include + using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- @@ -559,6 +564,12 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) input->one(cmd); error->done(); } + + // used in Variable::random_seed() + if(me == 0){ + // NOTE: better random generator for random_seed()? + srand(time(NULL)); + } } /* ---------------------------------------------------------------------- diff --git a/src/variable.cpp b/src/variable.cpp index 6e44a8f5a..dd1e4a9ef 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -41,6 +41,8 @@ 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. + + Tóth János (MATE, Gödöllő) ------------------------------------------------------------------------- */ #include @@ -69,6 +71,7 @@ #include "memory.h" #include "error.h" #include "force.h" +#include "math_extra_liggghts.h" #if defined(_WIN32) || defined(_WIN64) #include @@ -99,7 +102,7 @@ enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,MODULO,UNARY, SQRT,EXP,LN,LOG,ABS,SIN,COS,TAN,ASIN,ACOS,ATAN,ATAN2, RANDOM,NORMAL,CEIL,FLOOR,ROUND,RAMP,STAGGER,LOGFREQ,STRIDE, VDISPLACE,SWIGGLE,CWIGGLE,GMASK,RMASK,GRMASK, - VALUE,ATOMARRAY,TYPEARRAY,INTARRAY}; + VALUE,ATOMARRAY,TYPEARRAY,INTARRAY,SEED}; // customize by adding a special function @@ -1673,7 +1676,8 @@ double Variable::evaluate(char *str, Tree **tree) atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(), ramp(x,y),stagger(x,y),logfreq(x,y,z),stride(x,y,z), vdisplace(x,y),swiggle(x,y,z),cwiggle(x,y,z), - gmask(x),rmask(x),grmask(x,y) + gmask(x),rmask(x),grmask(x,y), + seed(n) ---------------------------------------------------------------------- */ double Variable::collapse_tree(Tree *tree) @@ -2107,6 +2111,17 @@ double Variable::collapse_tree(Tree *tree) return tree->value; } + if (tree->type == SEED) { + int ivalue1 = static_cast (collapse_tree(tree->left)); + + if (ivalue1 < 0) + error->one(FLERR,"seed() must be called with a non-negative argument (>= 0)"); + + tree->type = VALUE; + tree->value = generate_seed(ivalue1); + return tree->value; + } + // mask functions do not become a single collapsed value if (tree->type == GMASK) return 0.0; @@ -2124,7 +2139,8 @@ double Variable::collapse_tree(Tree *tree) atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(), ramp(x,y),stagger(x,y),logfreq(x,y,z),stride(x,y,z), vdisplace(x,y),swiggle(x,y,z),cwiggle(x,y,z), - gmask(x),rmask(x),grmask(x,y) + gmask(x),rmask(x),grmask(x,y), + seed(n) ---------------------------------------------------------------------- */ double Variable::eval_tree(Tree *tree, int i) @@ -2368,6 +2384,15 @@ double Variable::eval_tree(Tree *tree, int i) return arg; } + if (tree->type == SEED){ + int ivalue1 = static_cast (eval_tree(tree->left,i)); + + if (ivalue1 < 0) + error->one(FLERR,"seed() must be called with a non-negative argument (>= 0)"); + + return generate_seed(ivalue1); + } + if (tree->type == GMASK) { if (atom->mask[i] & tree->ivalue1) return 1.0; else return 0.0; @@ -2496,7 +2521,8 @@ int Variable::math_function(char *word, char *contents, Tree **tree, strcmp(word,"ramp") && strcmp(word,"stagger") && strcmp(word,"logfreq") && strcmp(word,"stride") && strcmp(word,"vdisplace") && - strcmp(word,"swiggle") && strcmp(word,"cwiggle")) + strcmp(word,"swiggle") && strcmp(word,"cwiggle") && + strcmp(word,"seed")) return 0; // parse contents for arg1,arg2,arg3 separated by commas @@ -2819,6 +2845,19 @@ int Variable::math_function(char *word, char *contents, Tree **tree, double value = value1 + value2*(1.0-cos(omega*delta*update->dt)); argstack[nargstack++] = value; } + + } else if (strcmp(word,"seed") == 0) { + if (narg != 1) + error->all(FLERR,"seed() must be called with 1 argument"); + if (tree) newtree->type = SEED; + else { + int ivalue1 = static_cast (value1); + + if(ivalue1 < 0) + error->all(FLERR,"seed() must be called with a non-negative argument (>= 0)"); + + argstack[nargstack++] = generate_seed(ivalue1); + } } delete [] arg1; @@ -3931,6 +3970,114 @@ unsigned int Variable::data_mask(char *str) return datamask; } +/* ---------------------------------------------------------------------- + generates a random prime number between 10'000 and INT_MAX + in a form of 6k +- 1 +------------------------------------------------------------------------- */ +static int generate_random_seed(){ + static const int MIN_K = 1678; // 6*1678-1 = 10067, prime + static const int MAX_K = 357913904; // 6*357913904-1 = 2147483423, prime + + static const int RANGE_K = MAX_K - MIN_K + 1; + + int p = 0; + + // it will find a prime + while(true){ + int k = (rand() % RANGE_K) + MIN_K; + p = 6 * k - 1; + + if(MathExtraLiggghts::isPrime(p)){ + break; + } + + p += 2; + + if(MathExtraLiggghts::isPrime(p)){ + break; + } + } + + return p; +} + +/* + first 200 prime greater than 10'000: + */ +static const int PRIME_TABLE_COUNT = 200; +static const int PRIME_TABLE[PRIME_TABLE_COUNT] = { + 10007, 10009, 10037, 10039, 10061, 10067, 10069, 10079, 10091, 10093, 10099, + 10103, 10111, 10133, 10139, 10141, 10151, 10159, 10163, 10169, 10177, 10181, + 10193, 10211, 10223, 10243, 10247, 10253, 10259, 10267, 10271, 10273, 10289, + 10301, 10303, 10313, 10321, 10331, 10333, 10337, 10343, 10357, 10369, 10391, + 10399, 10427, 10429, 10433, 10453, 10457, 10459, 10463, 10477, 10487, 10499, + 10501, 10513, 10529, 10531, 10559, 10567, 10589, 10597, 10601, 10607, 10613, + 10627, 10631, 10639, 10651, 10657, 10663, 10667, 10687, 10691, 10709, 10711, + 10723, 10729, 10733, 10739, 10753, 10771, 10781, 10789, 10799, 10831, 10837, + 10847, 10853, 10859, 10861, 10867, 10883, 10889, 10891, 10903, 10909, 10937, + 10939, 10949, 10957, 10973, 10979, 10987, 10993, 11003, 11027, 11047, 11057, + 11059, 11069, 11071, 11083, 11087, 11093, 11113, 11117, 11119, 11131, 11149, + 11159, 11161, 11171, 11173, 11177, 11197, 11213, 11239, 11243, 11251, 11257, + 11261, 11273, 11279, 11287, 11299, 11311, 11317, 11321, 11329, 11351, 11353, + 11369, 11383, 11393, 11399, 11411, 11423, 11437, 11443, 11447, 11467, 11471, + 11483, 11489, 11491, 11497, 11503, 11519, 11527, 11549, 11551, 11579, 11587, + 11593, 11597, 11617, 11621, 11633, 11657, 11677, 11681, 11689, 11699, 11701, + 11717, 11719, 11731, 11743, 11777, 11779, 11783, 11789, 11801, 11807, 11813, + 11821, 11827, 11831, 11833, 11839, 11863, 11867, 11887, 11897, 11903, 11909, + 11923, 11927 +}; + +static int generate_nth_prime_impl(int n){ + int p = PRIME_TABLE[PRIME_TABLE_COUNT - 1]; + n -= PRIME_TABLE_COUNT; + + while(n > 0){ + p += 2; + + if(MathExtraLiggghts::isPrime(p)){ + n -= 1; + } + } + + return p; +} + +/* ---------------------------------------------------------------------- + gives the n-th prime number between 10'000 and INT_MAX + + n p + 1 10007 + 2 10009 + 3 10037 + ... +------------------------------------------------------------------------- */ +static int generate_nth_prime(int n){ + // fast + if(n <= PRIME_TABLE_COUNT){ + return PRIME_TABLE[n-1]; + // small + } else { + return generate_nth_prime_impl(n); + } +} + +int Variable::generate_seed(int n){ + + int seed = 0; + + if(me == 0){ + if (n == 0){ + seed = generate_random_seed(); + } else { + seed = generate_nth_prime(n); + } + } + + MPI_Bcast(&seed,1,MPI_INT,0,world); + + return seed; +} + /* ---------------------------------------------------------------------- class to read variable values from a file for flag = SCALARFILE, reads one value per line diff --git a/src/variable.h b/src/variable.h index 328bfbafa..c538a770e 100644 --- a/src/variable.h +++ b/src/variable.h @@ -41,6 +41,8 @@ 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. + + Tóth János (MATE, Gödöllő) ------------------------------------------------------------------------- */ #ifndef LMP_VARIABLE_H @@ -128,6 +130,8 @@ class Variable : protected Pointers { int inumeric(char *); char *find_next_comma(char *); void print_tree(Tree *, int); + + int generate_seed(int); }; class VarReader : protected Pointers {