diff --git a/lib/atc/LinearSolver.cpp b/lib/atc/LinearSolver.cpp index ea31a64a7..655d6130c 100644 --- a/lib/atc/LinearSolver.cpp +++ b/lib/atc/LinearSolver.cpp @@ -1,564 +1,564 @@ // Header file for this class #include "LinearSolver.h" #include using std::stringstream; using std::set; namespace ATC { const double kPenalty = 1.0e4; const double kTol = 1.0e-8; const int kMaxDirect = 1000; // ==================================================================== // LinearSolver // ==================================================================== LinearSolver::LinearSolver( const SPAR_MAT & A, const BC_SET & bcs, const int solverType, const int constraintHandlerType, bool parallel ) : solverType_(solverType), constraintHandlerType_(constraintHandlerType), nVariables_(0), initialized_(false), initializedMatrix_(false), initializedInverse_(false), matrixModified_(false), allowReinitialization_(false), homogeneousBCs_(false), bcs_(&bcs), rhs_(NULL), rhsDense_(), b_(NULL), matrix_(A), matrixDense_(), matrixFreeFree_(), matrixFreeFixed_(),matrixInverse_(), penalty_(1), maxIterations_(0), maxRestarts_(0), tol_(0), parallel_(parallel) { // deep copy matrixCopy_ = A; matrixSparse_ = &matrixCopy_; setup(); } LinearSolver::LinearSolver( const SPAR_MAT & A, const int solverType, bool parallel ) : solverType_(solverType), constraintHandlerType_(NO_CONSTRAINTS), nVariables_(0), initialized_(false), initializedMatrix_(true), initializedInverse_(false), matrixModified_(false), allowReinitialization_(false), homogeneousBCs_(false), - bcs_(NULL), // null implies no contraints will be added later + bcs_(NULL), // null implies no constraints will be added later rhs_(NULL), rhsDense_(), b_(NULL), matrix_(A), matrixDense_(), matrixFreeFree_(), matrixFreeFixed_(),matrixInverse_(), penalty_(1), maxIterations_(0), maxRestarts_(0), tol_(0), parallel_(parallel) { // shallow copy matrixSparse_ = &A; setup(); } // -------------------------------------------------------------------- // Setup // -------------------------------------------------------------------- void LinearSolver::setup(void) { tol_ = kTol; nVariables_ = matrix_.nRows(); maxIterations_=2*nVariables_; maxRestarts_=nVariables_; // switch method based on size if (solverType_ < 0) { if (nVariables_ > kMaxDirect ) { solverType_ = ITERATIVE_SOLVE_SYMMETRIC; constraintHandlerType_ = PENALIZE_CONSTRAINTS; } else { solverType_ = DIRECT_SOLVE; } } if (constraintHandlerType_ < 0) { constraintHandlerType_ = PENALIZE_CONSTRAINTS; if (solverType_ == DIRECT_SOLVE) constraintHandlerType_ = CONDENSE_CONSTRAINTS; } if ( solverType_ == DIRECT_SOLVE && constraintHandlerType_ == CONDENSE_CONSTRAINTS ) allowReinitialization_ = true; if ( solverType_ == ITERATIVE_SOLVE_SYMMETRIC && constraintHandlerType_ == CONDENSE_CONSTRAINTS ) { throw ATC_Error("LinearSolver::unimplemented method"); } } // -------------------------------------------------------------------- // Initialize // -------------------------------------------------------------------- void LinearSolver::allow_reinitialization(void) { if (constraintHandlerType_ == PENALIZE_CONSTRAINTS) { if (matrixModified_ ) throw ATC_Error("LinearSolver: can't allow reinitialization after matrix has been modified"); matrixOriginal_ = *matrixSparse_; } allowReinitialization_ = true; } void LinearSolver::initialize(const BC_SET * bcs) { if (bcs) { if (! allowReinitialization_ ) throw ATC_Error("LinearSolver: reinitialization not allowed"); //if (! bcs_ ) throw ATC_Error("LinearSolver: adding constraints after constructing without constraints is not allowed"); // shallow --> deep copy if (! bcs_ ) { // constraintHandlerType_ == NO_CONSTRAINTS if (matrixModified_) { throw ATC_Error("LinearSolver: adding constraints after constructing without constraints is not allowed if matrix has been modified"); } else { matrixCopy_ = *matrixSparse_; matrixSparse_ = &matrixCopy_; constraintHandlerType_ = -1; setup(); } } bcs_ = bcs; initializedMatrix_ = false; initializedInverse_ = false; if (matrixModified_) { matrixCopy_ = matrixOriginal_; matrixSparse_ = &matrixCopy_; } } initialize_matrix(); initialize_inverse(); initialize_rhs(); initialized_ = true; } // -------------------------------------------------------------------- // initialize_matrix // -------------------------------------------------------------------- void LinearSolver::initialize_matrix(void) { if ( initializedMatrix_ ) return; if (constraintHandlerType_ == PENALIZE_CONSTRAINTS) { add_matrix_penalty(); } else if (constraintHandlerType_ == CONDENSE_CONSTRAINTS) { partition_matrix(); } initializedMatrix_ = true; } // -------------------------------------------------------------------- // initialize_inverse // -------------------------------------------------------------------- void LinearSolver::initialize_inverse(void) { if ( initializedInverse_ ) return; if (solverType_ == ITERATIVE_SOLVE_SYMMETRIC || solverType_ == ITERATIVE_SOLVE ) { matrixDiagonal_ = matrixSparse_->diag(); // preconditioner } else { // DIRECT_SOLVE if (constraintHandlerType_ == CONDENSE_CONSTRAINTS) { if( num_unknowns() > 0 ) { matrixInverse_ = inv(matrixFreeFree_); } } else { // NO_CONSTRAINTS || PENALIZE_CONSTRAINTS matrixDense_ = matrixSparse_->dense_copy(); // need dense for lapack matrixInverse_ = inv(matrixDense_); } } initializedInverse_ = true; } // -------------------------------------------------------------------- // initialize_rhs // -------------------------------------------------------------------- void LinearSolver::initialize_rhs(void) { if (! rhs_ ) return; if (! bcs_ ) { b_ = rhs_; return; } if (constraintHandlerType_ == PENALIZE_CONSTRAINTS) { add_rhs_penalty(); } else if (constraintHandlerType_ == CONDENSE_CONSTRAINTS) { add_rhs_influence(); } } // -------------------------------------------------------------------- // add matrix penalty // - change matrix for Dirichlet conditions: add penalty // -------------------------------------------------------------------- void LinearSolver::add_matrix_penalty(void) { penalty_ = kPenalty; // relative to matrix diagonal SPAR_MAT & A = matrixCopy_; penalty_ *= (A.diag()).maxabs(); BC_SET::const_iterator itr; for (itr = bcs_->begin(); itr != bcs_->end(); itr++) { int i = itr->first; A.add(i,i,penalty_); // modifies matrix } A.compress(); matrixModified_ = true; } // -------------------------------------------------------------------- // partition matrix // - partition matrix based on Dirichlet constraints // -------------------------------------------------------------------- void LinearSolver::partition_matrix(void) { fixedSet_.clear(); BC_SET::const_iterator itr; for (itr = bcs_->begin(); itr != bcs_->end(); itr++) { int i = itr->first; fixedSet_.insert(i); } freeSet_.clear(); freeGlobalToCondensedMap_.clear(); int j = 0; // local index for (int i = 0; i < nVariables_; i++) { if (fixedSet_.find(i) == fixedSet_.end() ) { freeSet_.insert(i); freeGlobalToCondensedMap_[i] = j++; } } if (matrixDense_.nRows() == 0) matrixDense_ =matrixSparse_->dense_copy(); DENS_MAT & K = matrixDense_; K.row_partition(freeSet_,matrixFreeFree_,matrixFreeFixed_); } // -------------------------------------------------------------------- // add_rhs_penalty // -------------------------------------------------------------------- void LinearSolver::add_rhs_penalty() { // deep copy VECTOR & b = rhsDense_; const VECTOR & r = *rhs_; int size = r.nRows(); b.reset(size); for (int i = 0; i < size; i++) { b(i) = r(i); } if ( ! homogeneousBCs_ ){ BC_SET::const_iterator itr; for (itr = bcs_->begin(); itr != bcs_->end(); itr++) { int i = itr->first; double v = itr->second; b(i) += penalty_ * v; } } b_ = &rhsDense_; } // -------------------------------------------------------------------- // add_rhs_influence // -------------------------------------------------------------------- void LinearSolver::add_rhs_influence() { if (! initializedMatrix_ ) partition_matrix(); // rhs = rhs + K_free,fixed * x_fixed int nbcs = bcs_->size(); if (nbcs == 0) { // no bcs to handle b_ = rhs_; } else { DENS_VEC & b = rhsDense_; if ( ! homogeneousBCs_ ){ DENS_VEC xFixed(nbcs); BC_SET::const_iterator itr; int i = 0; for (itr = bcs_->begin(); itr != bcs_->end(); itr++,i++) { double v = itr->second; xFixed(i,0) = -v; } b = matrixFreeFixed_*xFixed; // matrix and bcs have same ordering } else { b.reset(matrixFreeFixed_.nRows()); } const VECTOR & r = *rhs_; set::const_iterator iter; int i = 0; for (iter = freeSet_.begin(); iter != freeSet_.end(); iter++,i++) { b(i) += r(*iter); } b_ = &rhsDense_; } } // -------------------------------------------------------------------- // set fixed values // - {x_i = y_i} // -------------------------------------------------------------------- void LinearSolver::set_fixed_values(VECTOR & X) { BC_SET::const_iterator itr; for (itr = bcs_->begin(); itr != bcs_->end(); itr++) { int i = itr->first; double v = 0; if ( ! homogeneousBCs_ ) v = itr->second; X(i) = v; } } // -------------------------------------------------------------------- // Eigensystem // -------------------------------------------------------------------- // calls lapack void LinearSolver::eigen_system( DENS_MAT & eigenvalues, DENS_MAT & eigenvectors, const DENS_MAT * M) /* const */ { initialize_matrix(); // no inverse needed const DENS_MAT * Kp = NULL; const DENS_MAT * Mp =M; DENS_MAT MM; DENS_MAT KM; if (constraintHandlerType_ == CONDENSE_CONSTRAINTS) { Kp = &matrixFreeFree_; if (M) { DENS_MAT MfreeFixed; // not used M->row_partition(freeSet_,MM,MfreeFixed); Mp = &MM; } } else { if (matrixDense_.nRows() == 0) matrixDense_ =matrixSparse_->dense_copy(); Kp = &matrixDense_; } if (!M) { MM.identity(Kp->nRows()); Mp = &MM; } DENS_MAT eVecs, eVals; eVecs = eigensystem(*Kp,*Mp,eVals); eigenvalues.reset(nVariables_,1); eigenvectors.reset(nVariables_,nVariables_); set::const_iterator itr; for (int i = 0; i < Kp->nRows(); i++) { // ordering is by energy not node eigenvalues(i,0) = eVals(i,0); int j = 0; for (itr = freeSet_.begin(); itr != freeSet_.end(); itr++,j++) { int jj = *itr; eigenvectors(jj,i) = eVecs(j,i); // transpose } } } // -------------------------------------------------------------------- // solve // - solves A x = b // - if a "b" is provided it is used as the new rhs // -------------------------------------------------------------------- bool LinearSolver::solve(VECTOR & x, const VECTOR & b) { SPAR_MAT * A = NULL; rhs_ = &b; initialized_ = false; initialize(); if (num_unknowns() == 0) { set_fixed_values(x); return true; } const VECTOR & r = *b_; if (solverType_ == ITERATIVE_SOLVE_SYMMETRIC) { if (parallel_) { A = new PAR_SPAR_MAT(LammpsInterface::instance()->world(), *matrixSparse_); } else { A = new SPAR_MAT(*matrixSparse_); } DIAG_MAT & PC = matrixDiagonal_; int niter = maxIterations_; double tol = tol_; int convergence = CG(*A, x, r, PC, niter, tol);// CG changes niter, tol if (convergence>0) { stringstream ss; ss << "CG solve did not converge,"; ss << " iterations: " << niter; ss << " residual: " << tol; throw ATC_Error(ss.str()); } } else if (solverType_ == ITERATIVE_SOLVE) { if (parallel_) { A = new PAR_SPAR_MAT(LammpsInterface::instance()->world(), *matrixSparse_); } else { A = new SPAR_MAT(*matrixSparse_); } const DIAG_MAT & PC = matrixDiagonal_; int iterations = maxIterations_; int restarts = maxRestarts_; double tol = tol_; DENS_MAT H(maxRestarts_+1, maxRestarts_); DENS_VEC xx(nVariables_); DENS_VEC bb; bb = b; int convergence = GMRES(*A, xx, bb, PC, H, restarts, iterations, tol); if (convergence>0) { stringstream ss; ss << "GMRES greens_function solve did not converge,"; ss << " iterations: " << iterations; ss << " residual: " << tol; throw ATC_Error(ss.str()); } x.copy(xx.ptr(),xx.nRows()); } else { // DIRECT_SOLVE const DENS_MAT & invA = matrixInverse_; if (constraintHandlerType_ == CONDENSE_CONSTRAINTS) { DENS_MAT xx = invA*r; int i = 0; set::const_iterator itr; for (itr = freeSet_.begin(); itr != freeSet_.end(); itr++,i++) { int ii = *itr; x(ii) = xx(i,0); } set_fixed_values(x); } else { DENS_VEC xx = invA*r; for (int i = 0; i < xx.nRows(); i++) { x(i) = xx(i); } } } delete A; return true; } // -------------------------------------------------------------------- // greens function // - returns the solution to a Kronecker delta rhs b = {0 0 .. 1 .. 0 0} // and with homogeneous constraints {x_i = 0} // -------------------------------------------------------------------- void LinearSolver::greens_function(int I, VECTOR & G_I) { SPAR_MAT * A = NULL; initialize_matrix(); initialize_inverse(); G_I.reset(nVariables_); VECTOR & x = G_I; if (solverType_ == ITERATIVE_SOLVE_SYMMETRIC) { DENS_VEC b(nVariables_); b = 0.0; b(I) = 1.0; if (parallel_) { A = new PAR_SPAR_MAT(LammpsInterface::instance()->world(), *matrixSparse_); } else { A = new SPAR_MAT(*matrixSparse_); } const DIAG_MAT & PC = matrixDiagonal_; int niter = maxIterations_; double tol = tol_; int convergence = CG(*A, x, b, PC, niter, tol); if (convergence>0) { stringstream ss; ss << "CG greens_function solve did not converge,"; ss << " iterations: " << niter; ss << " residual: " << tol; throw ATC_Error(ss.str()); } } else if (solverType_ == ITERATIVE_SOLVE) { DENS_VEC b(nVariables_); b = 0.0; b(I) = 1.0; // VECTOR & bb = b; if (parallel_) { A = new PAR_SPAR_MAT(LammpsInterface::instance()->world(), *matrixSparse_); } else { A = new SPAR_MAT(*matrixSparse_); } // const DENS_MAT A = matrixSparse_->dense_copy(); const DIAG_MAT & PC = matrixDiagonal_; int iterations = maxIterations_; int restarts = maxRestarts_; double tol = tol_; DENS_MAT H(maxRestarts_+1, maxRestarts_); DENS_VEC xx(nVariables_); int convergence = GMRES(*A, xx, b, PC, H, restarts, iterations, tol); if (convergence>0) { stringstream ss; ss << "GMRES greens_function solve did not converge,"; ss << " iterations: " << iterations; ss << " residual: " << tol; throw ATC_Error(ss.str()); } x.copy(xx.ptr(),xx.nRows()); } else { const DENS_MAT & invA = matrixInverse_; if (constraintHandlerType_ == CONDENSE_CONSTRAINTS) { set::const_iterator itr; for (itr = fixedSet_.begin(); itr != fixedSet_.end(); itr++) { int ii = *itr; x(ii) = 0; } itr = freeSet_.find(I); if (itr !=freeSet_.end() ) { int j = freeGlobalToCondensedMap_[I]; int i = 0; for (itr = freeSet_.begin(); itr != freeSet_.end(); itr++,i++) { int ii = *itr; x(ii) = invA(j,i); } } } else { for (int i = 0; i < nVariables_; ++i) x(i) = invA(I,i); } } delete A; } } // namespace ATC diff --git a/python/README b/python/README index f0e95d3c4..86e599a40 100644 --- a/python/README +++ b/python/README @@ -1,130 +1,130 @@ This directory contains Python code which wraps LAMMPS as a library and allows the LAMMPS library interface to be invoked from Python, either from a Python script or using Python interactively. Details on the Python interface to LAMMPS and how to build LAMMPS as a shared library, for use with Python, are given in doc/Section_python.html and in doc/Section_start.html#start_5. Basically you need to follow these steps in the src directory: % make g++ mode=shlib # build for whatever machine target you wish % make install-python # may need to do this via sudo You can replace the last step by a one-time setting of environment variables in your shell script. Or you can run the python/install.py script directly to give you more control over where the two relevant files are installed. See doc/Section_python.html for details. You should then be able to launch Python and instantiate an instance of LAMMPS: % python >>> from lammps import lammps >>> lmp = lammps() -If that gives no errors, you have succesfully wrapped LAMMPS with +If that gives no errors, you have successfully wrapped LAMMPS with Python. See doc/Section_python.html#py_7 for tests you can then use to run LAMMPS both in serial or parallel thru Python. Note that you can also invoke Python code from within a LAMMPS input script, using the "python" command. See the doc/python.html doc page for details. The Python code you invoke can also call back to LAMMPS using the same interface described here for wrapping LAMMPS. ------------------------------------------------------------------- Once you have successfully wrapped LAMMPS, you can run the Python scripts in the examples sub-directory: trivial.py read/run a LAMMPS input script thru Python demo.py invoke various LAMMPS library interface routines simple.py parallel example, mimicing examples/COUPLE/simple/simple.cpp split.py parallel example mc.py Monte Carlo energy relaxation wrapper on LAMMPS gui.py GUI go/stop/temperature-slider to control LAMMPS plot.py real-time temeperature plot with GnuPlot via Pizza.py viz_tool.py real-time viz via some viz package vizplotgui_tool.py combination of viz.py and plot.py and gui.py For the viz_tool.py and vizplotgui_tool.py commands, replace "tool" with "gl" or "atomeye" or "pymol", depending on what visualization package you have installed. We hope to add a VMD option soon. Note that for GL, you need to be able to run the Pizza.py GL tool, which is included in the pizza sub-directory. See the Pizza.py doc pages for more info: http://www.sandia.gov/~sjplimp/pizza.html Note that for AtomEye, you need version 3, and their is a line in the scripts that specifies the path and name of the executable. See the AtomEye WWW pages for more details: http://mt.seas.upenn.edu/Archive/Graphics/A http://mt.seas.upenn.edu/Archive/Graphics/A3/A3.html The latter link is to AtomEye 3 which has the scriping capability needed by these Python scripts. Note that for PyMol, you need to have built and installed the open-source version of PyMol in your Python, so that you can import it from a Python script. See the PyMol WWW pages for more details: http://www.pymol.org http://sourceforge.net/scm/?type=svn&group_id=4546 The latter link is to the open-source version. ------------------------------------------------------------------- Each example script has more documentation in the file that explains how to use it and what it is doing. You can run a particular script in either of the following ways: % trivial.py in.trivial % python -i trivial.py in.trivial The former assumes that you have changed the first line of the script to point to the Python installed on your box and made the script exectable (e.g. chmod +x trivial.py). The example scripts take the following arguments. The in.* args are LAMPS input scripts. trivial.py in.trivial demo.py simple.py in.simple # can run in parallel (see below) split.py in.simple # can run in parallel (see below) gui.py in.gui 100 plot.py in.plot 10 1000 thermo_temp viz_tool.py in.viz 100 5000 vizplotgui_tool.py in.viz 100 thermo_temp To run LAMMPS in parallel from Python, so something like this: % mpirun -np 4 simple.py in.simple % mpirun -np 4 python split.py in.simple If you run simple.py as-is, this will invoke P instances of a one-processor run, where both Python and LAMMPS will run on single processors. Each running job will read the same input file, and write to same log.lammps file, which isn't too useful. However, if you have either the Pypar or mpi4py packages installed in your Python, and uncomment the Pypar or mpi4py code in simple.py, then the above commands will invoke 1 instance of a P-processor run. Both Python and LAMMPS will run on P processors. The job will read the input file and write a single log.lammps file. The split.py script can also be run in parallel. It uses mpi4py version 2.0.0 (or later), which makes it possible to pass a communicator when creating the LAMMPS object and thus run multiple instances of LAMMPS at the same time, each on a different subset of MPI ranks. Or run LAMMPS on one subset and some other program on the rest of the MPI ranks, concurrently. See comments in the split.py script for more details. diff --git a/src/ASPHERE/compute_erotate_asphere.h b/src/ASPHERE/compute_erotate_asphere.h index be1c53492..23e30d7c9 100644 --- a/src/ASPHERE/compute_erotate_asphere.h +++ b/src/ASPHERE/compute_erotate_asphere.h @@ -1,61 +1,61 @@ /* -*- c++ -*- ---------------------------------------------------------- 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 COMPUTE_CLASS ComputeStyle(erotate/asphere,ComputeERotateAsphere) #else #ifndef LMP_COMPUTE_EROTATE_ASPHERE_H #define LMP_COMPUTE_EROTATE_ASPHERE_H #include "compute.h" namespace LAMMPS_NS { class ComputeERotateAsphere : public Compute { public: ComputeERotateAsphere(class LAMMPS *, int, char **); void init(); double compute_scalar(); private: double pfactor; class AtomVecEllipsoid *avec_ellipsoid; class AtomVecLine *avec_line; class AtomVecTri *avec_tri; }; } #endif #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Compute erotate/asphere requires atom style ellipsoid or line or tri Self-explanatory. E: Compute erotate/asphere requires extended particles -This compute cannot be used with point paritlces. +This compute cannot be used with point particles. */ diff --git a/src/ASPHERE/compute_temp_asphere.h b/src/ASPHERE/compute_temp_asphere.h index 4688a5612..d1cce3802 100644 --- a/src/ASPHERE/compute_temp_asphere.h +++ b/src/ASPHERE/compute_temp_asphere.h @@ -1,92 +1,92 @@ /* -*- c++ -*- ---------------------------------------------------------- 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 COMPUTE_CLASS ComputeStyle(temp/asphere,ComputeTempAsphere) #else #ifndef LMP_COMPUTE_TEMP_ASPHERE_H #define LMP_COMPUTE_TEMP_ASPHERE_H #include "compute.h" namespace LAMMPS_NS { class ComputeTempAsphere : public Compute { public: ComputeTempAsphere(class LAMMPS *, int, char **); ~ComputeTempAsphere(); void init(); void setup(); double compute_scalar(); void compute_vector(); void remove_bias(int, double *); void restore_bias(int, double *); private: int mode; double tfactor; char *id_bias; class Compute *tbias; // ptr to additional bias compute class AtomVecEllipsoid *avec; void dof_compute(); }; } #endif #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Compute temp/asphere requires atom style ellipsoid Self-explanatory. E: Compute temp/asphere requires extended particles -This compute cannot be used with point paritlces. +This compute cannot be used with point particles. E: Could not find compute ID for temperature bias Self-explanatory. E: Bias compute does not calculate temperature The specified compute must compute temperature. E: Bias compute does not calculate a velocity bias The specified compute must compute a bias for temperature. E: Bias compute group does not match compute group The specified compute must operate on the same group as the parent compute. E: Temperature compute degrees of freedom < 0 This should not happen if you are calculating the temperature on a valid set of atoms. */ diff --git a/src/PYTHON/python.h b/src/PYTHON/python.h index 61e5f015e..5f65e3970 100644 --- a/src/PYTHON/python.h +++ b/src/PYTHON/python.h @@ -1,130 +1,130 @@ /* -*- c++ -*- ---------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifndef LMP_PYTHON_H #define LMP_PYTHON_H #include "pointers.h" namespace LAMMPS_NS { class Python : protected Pointers { public: int python_exists; Python(class LAMMPS *); ~Python(); void command(int, char **); void invoke_function(int, char *); int find(char *); int variable_match(char *, char *, int); char *long_string(int); private: int ninput,noutput,length_longstr; char **istr; char *ostr,*format; void *pyMain; struct PyFunc { char *name; int ninput,noutput; int *itype,*ivarflag; int *ivalue; double *dvalue; char **svalue; int otype; char *ovarname; char *longstr; int length_longstr; void *pFunc; }; PyFunc *pfuncs; int nfunc; int create_entry(char *); void deallocate(int); }; } #endif /* ERROR/WARNING messages: E: Invalid python command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Python invoke of undefined function Cannot invoke a function that has not been previously defined. E: Python variable does not match Python function This matching is defined by the python-style variable and the python command. E: Cannot embed Python when also extending Python with LAMMPS When running LAMMPS via Python through the LAMMPS library interface you cannot also user the input script python command. E: Could not initialize embedded Python The main module in Python was not accessible. E: Could not open Python file The specified file of Python code cannot be opened. Check that the path and name are correct. E: Could not process Python file -The Python code in the specified file was not run sucessfully by +The Python code in the specified file was not run successfully by Python, probably due to errors in the Python code. E: Could not process Python string -The Python code in the here string was not run sucessfully by Python, +The Python code in the here string was not run successfully by Python, probably due to errors in the Python code. E: Could not find Python function The provided Python code was run successfully, but it not define a callable function with the required name. E: Python function is not callable The provided Python code was run successfully, but it not define a callable function with the required name. E: Could not create Python function arguments This is an internal Python error, possibly because the number of inputs to the function is too large. E: Could not evaluate Python function input variable Self-explanatory. E: Python function evaluation failed -The Python function did not run succesfully and/or did not return a +The Python function did not run successfully and/or did not return a value (if it is supposed to return a value). This is probably due to some error condition in the function. */ diff --git a/src/SHOCK/fix_append_atoms.h b/src/SHOCK/fix_append_atoms.h index d9885998b..cc26acc8b 100644 --- a/src/SHOCK/fix_append_atoms.h +++ b/src/SHOCK/fix_append_atoms.h @@ -1,108 +1,108 @@ /* -*- c++ -*- ---------------------------------------------------------- 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(append/atoms,FixAppendAtoms) #else #ifndef FIX_APPEND_ATOMS_H #define FIX_APPEND_ATOMS_H #include "fix.h" namespace LAMMPS_NS { class FixAppendAtoms : public Fix { public: FixAppendAtoms(class LAMMPS *, int, char **); ~FixAppendAtoms(); int setmask(); void setup(int); void pre_exchange(); void initial_integrate(int); void post_force(int); private: int get_spatial(); int spatflag, xloflag, xhiflag, yloflag, yhiflag, zloflag, zhiflag; int ranflag, tempflag, xseed, tseed; double ranx, rany, ranz, t_target, t_period, t_extent; class RanMars *randomx; class RanMars *randomt; int scaleflag, freq; int nbasis; int *basistype; int advance,advance_sum; double size,spatlead; char *spatialid; double tfactor; double *gfactor1,*gfactor2; }; } #endif #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Fix append/atoms requires a lattice be defined Use the lattice command for this purpose. E: Only zhi currently implemented for fix append/atoms Self-explanatory. E: Append boundary must be shrink/minimum The boundary style of the face where atoms are added must be of type m (shrink/minimum). E: Bad fix ID in fix append/atoms command The value of the fix_id for keyword spatial must start with 'f_'. E: Invalid basis setting in fix append/atoms command The basis index must be between 1 to N where N is the number of basis atoms in the lattice. The type index must be between 1 to N where N is the number of atom types. E: Cannot use append/atoms in periodic dimension The boundary style of the face where atoms are added can not be of type p (periodic). E: Cannot append atoms to a triclinic box -The simulation box must be defined with edges alligned with the +The simulation box must be defined with edges aligned with the Cartesian axes. E: Fix ID for fix ave/spatial does not exist Self-explanatory. E: Too many total atoms See the setting for bigint in the src/lmptype.h file. */ diff --git a/src/SRD/fix_srd.h b/src/SRD/fix_srd.h index fae9780b8..e14652e52 100644 --- a/src/SRD/fix_srd.h +++ b/src/SRD/fix_srd.h @@ -1,434 +1,434 @@ /* -*- c++ -*- ---------------------------------------------------------- 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(srd,FixSRD) #else #ifndef LMP_FIX_SRD_H #define LMP_FIX_SRD_H #include "fix.h" namespace LAMMPS_NS { class FixSRD : public Fix { public: FixSRD(class LAMMPS *, int, char **); ~FixSRD(); int setmask(); void init(); void setup(int); void pre_neighbor(); void post_force(int); double compute_vector(int); double memory_usage(); int pack_reverse_comm(int, int, double *); void unpack_reverse_comm(int, int *, double *); private: int me,nprocs; int bigexist,biggroup,biggroupbit; int collidestyle,lamdaflag,overlap,insideflag,exactflag,maxbounceallow; int cubicflag,shiftuser,shiftseed,shiftflag,tstat; int rescale_rotate,rescale_collide; double gridsrd,gridsearch,lamda,radfactor,cubictol; int triclinic,change_size,change_shape,deformflag; double dt_big,dt_srd; double mass_big,mass_srd; double temperature_srd; double sigma; double srd_per_cell; double dmax,vmax,vmaxsq; double maxbigdiam,minbigdiam; double dist_ghost,dist_srd,dist_srd_reneigh; // explained in code int wallexist,nwall,wallvarflag; class FixWallSRD *wallfix; int *wallwhich; double *xwall,*xwallhold,*vwall; double **fwall; double walltrigger; class AtomVecEllipsoid *avec_ellipsoid; class AtomVecLine *avec_line; class AtomVecTri *avec_tri; // for orthogonal box, these are in box units // for triclinic box, these are in lamda units double srdlo[3],srdhi[3]; // SRDs must stay inside double srdlo_reneigh[3],srdhi_reneigh[3]; // SRDs trigger a reneigh int dimension; int initflag,setupflag,reneighflag; class RanMars *random; class RanPark *randomshift; // stats int ncheck,ncollide,ninside,nrescale,reneighcount; int nbounce,bouncemaxnum,bouncemax; int stats_flag; int srd_bin_count; double srd_bin_temp; double stats[12],stats_all[12]; double **flocal; // local ptrs to atom force and torque double **tlocal; // info to store for each owned and ghost big particle and wall struct Big { int index; // local index of particle/wall int type; // SPHERE or ELLIPSOID or LINE or TRI or WALL double radius,radsq; // radius of sphere double aradsqinv; // 3 ellipsoid radii double bradsqinv; double cradsqinv; double length; // length of line segment double normbody[3]; // normal of tri in body-frame double cutbinsq; // add big to bin if within this distance double omega[3]; // current omega for sphere/ellipsoid/tri/line double ex[3],ey[3],ez[3]; // current orientation vecs for ellipsoid/tri double norm[3]; // current unit normal of tri in space-frame double theta; // current orientation of line }; Big *biglist; // list of info for each owned & ghost big and wall int torqueflag; // 1 if any big particle is torqued // current size of particle-based arrays int nbig; // # of owned/ghost big particles and walls int maxbig; // max number of owned/ghost big particles and walls int nmax; // max number of SRD particles // bins for SRD velocity remap, shifting and communication // binsize and inv are in lamda units for triclinic int nbins1,nbin1x,nbin1y,nbin1z; double binsize1x,binsize1y,binsize1z; double bininv1x,bininv1y,bininv1z; struct BinAve { int owner; // 1 if I am owner of this bin, 0 if not int n; // # of SRD particles in bin double xctr[3]; // center point of bin, only used for triclinic double vsum[3]; // sum of v components for SRD particles in bin double random; // random value if I am owner double value[12]; // extra per-bin values }; struct BinComm { int nsend,nrecv; // # of bins to send/recv int sendproc,recvproc; // who to send/recv to/from int *sendlist,*recvlist; // list of bins to send/recv }; struct BinShift { int commflag; // 1 if this shift requires any comm int nbins,nbinx,nbiny,nbinz; // extent of my bins int maxbinsq,maxvbin; int binlo[3],binhi[3]; // extent of my bins in global array double corner[3]; // lower,left corner to offset from // corner is in lamda units for triclinic BinAve *vbin; // my bins BinComm bcomm[6]; // bin communication pattern for overlaps }; BinShift shifts[2]; // 0 = no shift, 1 = shift int maxbin1; int *binhead; // 1st SRD particle in each bin int *binnext; // next SRD particle in same bin int maxbuf; double *sbuf1,*sbuf2; // buffers for send/recv of velocity bin data double *rbuf1,*rbuf2; // bins and stencil for collision searching for SRDs & big particles int nbins2,nbin2x,nbin2y,nbin2z; int maxbin2; double binsize2x,binsize2y,binsize2z; double bininv2x,bininv2y,bininv2z; double xblo2,yblo2,zblo2; int *nbinbig; // # of big particles overlapping each bin int **binbig; // indices of big particles overlapping each bin int *binsrd; // which bin each SRD particle is in int nstencil; // # of bins in stencil int maxstencil; // max # of bins stencil array can hold int **stencil; // list of 3d bin offsets a big particle can overlap // persistent data for line/tri collision calculations double tfraction,theta0,theta1; double xs0[3],xs1[3],xsc[3]; double xb0[3],xb1[3],xbc[3]; double nbc[3]; // shared data for triangle collision calculations // private functions void reset_velocities(); void vbin_comm(int); void vbin_pack(BinAve *, int, int *, double *); void vbin_unpack(double *, BinAve *, int, int *); void xbin_comm(int, int); void xbin_pack(BinAve *, int, int *, double *, int); void xbin_unpack(double *, BinAve *, int, int *, int); void collisions_single(); void collisions_multi(); int inside_sphere(double *, double *, Big *); int inside_ellipsoid(double *, double *, Big *); int inside_line(double *, double *, double *, double *, Big *, double); int inside_tri(double *, double *, double *, double *, Big *, double); int inside_wall(double *, int); double collision_sphere_exact(double *, double *, double *, double *, Big *, double *, double *, double *); void collision_sphere_inexact(double *, double *, Big *, double *, double *, double *); double collision_ellipsoid_exact(double *, double *, double *, double *, Big *, double *, double *, double *); void collision_ellipsoid_inexact(double *, double *, Big *, double *, double *, double *); double collision_line_exact(double *, double *, double *, double *, Big *, double, double *, double *, double *); double collision_tri_exact(double *, double *, double *, double *, Big *, double, double *, double *, double *); double collision_wall_exact(double *, int, double *, double *, double *, double *); void collision_wall_inexact(double *, int, double *, double *, double *); void slip(double *, double *, double *, Big *, double *, double *, double *); void slip_wall(double *, int, double *, double *); void noslip(double *, double *, double *, Big *, int, double *, double *, double *); void force_torque(double *, double *, double *, double *, double *, double *); void force_wall(double *, double *, int); int update_srd(int, double, double *, double *, double *, double *); void parameterize(); void setup_bounds(); void setup_velocity_bins(); void setup_velocity_shift(int, int); void setup_search_bins(); void setup_search_stencil(); void big_static(); void big_dynamic(); double point_bin_distance(double *, int, int, int); double bin_bin_distance(int, int, int); void velocity_stats(int); double newton_raphson(double, double); void lineside(double, double &, double &); void triside(double, double &, double &); double distance(int, int); void print_collision(int, int, int, double, double, double *, double *, double *, int); }; } #endif #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Could not find fix srd group ID Self-explanatory. E: Fix srd requires newton pair on Self-explanatory. E: Fix srd requires ghost atoms store velocity Use the comm_modify vel yes command to enable this. E: Fix srd no-slip requires atom attribute torque This is because the SRD collisions will impart torque to the solute particles. E: Cannot change timestep once fix srd is setup This is because various SRD properties depend on the timestep size. E: Fix srd can only currently be used with comm_style brick This is a current restriction in LAMMPS. E: Cannot use fix wall/srd more than once Nor is their a need to since multiple walls can be specified in one command. W: Fix SRD walls overlap but fix srd overlap not set You likely want to set this in your input script. E: Using fix srd with inconsistent fix deform remap option When shearing the box in an SRD simulation, the remap v option for fix deform needs to be used. W: Using fix srd with box deformation but no SRD thermostat The deformation will heat the SRD particles so this can be dangerous. W: Fix srd SRD moves may trigger frequent reneighboring This is because the SRD particles may move long distances. E: Fix SRD: bad search bin assignment Something has gone wrong in your SRD model; try using more conservative settings. E: Fix SRD: bad stencil bin for big particle Something has gone wrong in your SRD model; try using more conservative settings. E: Fix SRD: too many big particles in bin Reset the ATOMPERBIN parameter at the top of fix_srd.cpp to a larger value, and re-compile the code. E: Fix SRD: too many walls in bin This should not happen unless your system has been setup incorrectly. E: Fix SRD: bad bin assignment for SRD advection Something has gone wrong in your SRD model; try using more conservative settings. E: SRD particle %d started inside big particle %d on step %ld bounce %d See the inside keyword if you want this message to be an error vs warning. W: SRD particle %d started inside big particle %d on step %ld bounce %d See the inside keyword if you want this message to be an error vs warning. E: SRD particle %d started inside wall %d on step %ld bounce %d See the inside keyword if you want this message to be an error vs warning. W: SRD particle %d started inside wall %d on step %ld bounce %d See the inside keyword if you want this message to be an error vs warning. E: Bad quadratic solve for particle/line collision -This is an internal error. It should nornally not occur. +This is an internal error. It should normally not occur. E: Bad quadratic solve for particle/tri collision -This is an internal error. It should nornally not occur. +This is an internal error. It should normally not occur. W: Fix srd particle moved outside valid domain This may indicate a problem with your simulation parameters. E: Big particle in fix srd cannot be point particle Big particles must be extended spheriods or ellipsoids. E: Cannot use lines with fix srd unless overlap is set -This is because line segements are connected to each other. +This is because line segments are connected to each other. E: Cannot use tris with fix srd unless overlap is set This is because triangles are connected to each other. E: Fix srd requires SRD particles all have same mass Self-explanatory. E: Fewer SRD bins than processors in some dimension This is not allowed. Make your SRD bin size smaller. E: SRD bins for fix srd are not cubic enough The bin shape is not within tolerance of cubic. See the cubic keyword if you want this message to be an error vs warning. W: SRD bins for fix srd are not cubic enough The bin shape is not within tolerance of cubic. See the cubic keyword if you want this message to be an error vs warning. E: SRD bin size for fix srd differs from user request Fix SRD had to adjust the bin size to fit the simulation box. See the cubic keyword if you want this message to be an error vs warning. W: SRD bin size for fix srd differs from user request Fix SRD had to adjust the bin size to fit the simulation box. See the cubic keyword if you want this message to be an error vs warning. E: Fix srd lamda must be >= 0.6 of SRD grid size This is a requirement for accuracy reasons. W: SRD bin shifting turned on due to small lamda This is done to try to preserve accuracy. W: Fix srd grid size > 1/4 of big particle diameter This may cause accuracy problems. W: Fix srd viscosity < 0.0 due to low SRD density This may cause accuracy problems. W: Fix srd particles may move > big particle diameter This may cause accuracy problems. */ diff --git a/src/angle.h b/src/angle.h index 119f206c3..8c8129bbb 100644 --- a/src/angle.h +++ b/src/angle.h @@ -1,82 +1,82 @@ /* -*- c++ -*- ---------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifndef LMP_ANGLE_H #define LMP_ANGLE_H #include #include "pointers.h" namespace LAMMPS_NS { class Angle : protected Pointers { friend class ThrOMP; friend class FixOMP; public: int allocated; int *setflag; int writedata; // 1 if writes coeffs to data file double energy; // accumulated energies - double virial[6]; // accumlated virial + double virial[6]; // accumulated virial double *eatom,**vatom; // accumulated per-atom energy/virial // KOKKOS host/device flag and data masks ExecutionSpace execution_space; unsigned int datamask_read,datamask_modify; int copymode; Angle(class LAMMPS *); virtual ~Angle(); virtual void init(); virtual void compute(int, int) = 0; virtual void settings(int, char **) {} virtual void coeff(int, char **) = 0; virtual void init_style() {}; virtual double equilibrium_angle(int) = 0; virtual void write_restart(FILE *) = 0; virtual void read_restart(FILE *) = 0; virtual void write_data(FILE *) {} virtual double single(int, int, int, int) = 0; virtual double memory_usage(); protected: int suffix_flag; // suffix compatibility flag int evflag; int eflag_either,eflag_global,eflag_atom; int vflag_either,vflag_global,vflag_atom; int maxeatom,maxvatom; void ev_setup(int, int); void ev_tally(int, int, int, int, int, double, double *, double *, double, double, double, double, double, double); }; } #endif /* ERROR/WARNING messages: E: Angle coeffs are not set No angle coefficients have been assigned in the data file or via the angle_coeff command. E: All angle coeffs are not set All angle coefficients must be set in the data file or by the angle_coeff command before running a simulation. */ diff --git a/src/bond.h b/src/bond.h index 06e05a726..ffa1f5fb0 100644 --- a/src/bond.h +++ b/src/bond.h @@ -1,83 +1,83 @@ /* -*- c++ -*- ---------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifndef LMP_BOND_H #define LMP_BOND_H #include #include "pointers.h" namespace LAMMPS_NS { class Bond : protected Pointers { friend class ThrOMP; friend class FixOMP; public: int allocated; int *setflag; int writedata; // 1 if writes coeffs to data file double energy; // accumulated energies - double virial[6]; // accumlated virial + double virial[6]; // accumulated virial double *eatom,**vatom; // accumulated per-atom energy/virial // KOKKOS host/device flag and data masks ExecutionSpace execution_space; unsigned int datamask_read,datamask_modify; int copymode; Bond(class LAMMPS *); virtual ~Bond(); virtual void init(); virtual void init_style() {} virtual void compute(int, int) = 0; virtual void settings(int, char **) {} virtual void coeff(int, char **) = 0; virtual double equilibrium_distance(int) = 0; virtual void write_restart(FILE *) = 0; virtual void read_restart(FILE *) = 0; virtual void write_data(FILE *) {} virtual double single(int, double, int, int, double &) = 0; virtual double memory_usage(); void write_file(int, char**); protected: int suffix_flag; // suffix compatibility flag int evflag; int eflag_either,eflag_global,eflag_atom; int vflag_either,vflag_global,vflag_atom; int maxeatom,maxvatom; void ev_setup(int, int); void ev_tally(int, int, int, int, double, double, double, double, double); }; } #endif /* ERROR/WARNING messages: E: Bond coeffs are not set No bond coefficients have been assigned in the data file or via the bond_coeff command. E: All bond coeffs are not set All bond coefficients must be set in the data file or by the bond_coeff command before running a simulation. */ diff --git a/src/citeme.h b/src/citeme.h index 80b642ab6..18e9a712d 100644 --- a/src/citeme.h +++ b/src/citeme.h @@ -1,48 +1,48 @@ /* -*- c++ -*- ---------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifndef LMP_CITEME_H #define LMP_CITEME_H #include "pointers.h" #include #include namespace LAMMPS_NS { class CiteMe : protected Pointers { public: CiteMe(class LAMMPS *); virtual ~CiteMe(); void add(const char *); // print out and register publication private: FILE *fp; // opaque pointer to log.cite file object typedef std::set citeset; citeset *cs; // registered set of publications }; } #endif /* ERROR/WARNING messages: E: Cannot open log.cite file This file is created when you use some LAMMPS features, to indicate what paper you should cite on behalf of those who implemented -the feature. Check that you have write priveleges into the directory +the feature. Check that you have write privileges into the directory you are running in. */ diff --git a/src/dihedral.h b/src/dihedral.h index 68167eb86..5f3909244 100644 --- a/src/dihedral.h +++ b/src/dihedral.h @@ -1,81 +1,81 @@ /* -*- c++ -*- ---------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifndef LMP_DIHEDRAL_H #define LMP_DIHEDRAL_H #include #include "pointers.h" namespace LAMMPS_NS { class Dihedral : protected Pointers { friend class ThrOMP; friend class FixOMP; public: int allocated; int *setflag; int writedata; // 1 if writes coeffs to data file double energy; // accumulated energy - double virial[6]; // accumlated virial + double virial[6]; // accumulated virial double *eatom,**vatom; // accumulated per-atom energy/virial // KOKKOS host/device flag and data masks ExecutionSpace execution_space; unsigned int datamask_read,datamask_modify; int copymode; Dihedral(class LAMMPS *); virtual ~Dihedral(); virtual void init(); virtual void init_style() {} virtual void compute(int, int) = 0; virtual void settings(int, char **) {} virtual void coeff(int, char **) = 0; virtual void write_restart(FILE *) = 0; virtual void read_restart(FILE *) = 0; virtual void write_data(FILE *) {} virtual double memory_usage(); protected: int suffix_flag; // suffix compatibility flag int evflag; int eflag_either,eflag_global,eflag_atom; int vflag_either,vflag_global,vflag_atom; int maxeatom,maxvatom; void ev_setup(int, int); void ev_tally(int, int, int, int, int, int, double, double *, double *, double *, double, double, double, double, double, double, double, double, double); }; } #endif /* ERROR/WARNING messages: E: Dihedral coeffs are not set No dihedral coefficients have been assigned in the data file or via the dihedral_coeff command. E: All dihedral coeffs are not set All dihedral coefficients must be set in the data file or by the dihedral_coeff command before running a simulation. */ diff --git a/src/fix.h b/src/fix.h index bd3189afe..8005da1ad 100644 --- a/src/fix.h +++ b/src/fix.h @@ -1,283 +1,283 @@ /* -*- c++ -*- ---------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifndef LMP_FIX_H #define LMP_FIX_H #include "pointers.h" namespace LAMMPS_NS { class Fix : protected Pointers { public: static int instance_total; // # of Fix classes ever instantiated char *id,*style; int igroup,groupbit; int restart_global; // 1 if Fix saves global state, 0 if not int restart_peratom; // 1 if Fix saves peratom state, 0 if not int restart_file; // 1 if Fix writes own restart file, 0 if not int force_reneighbor; // 1 if Fix forces reneighboring, 0 if not int box_change_size; // 1 if Fix changes box size, 0 if not int box_change_shape; // 1 if Fix changes box shape, 0 if not int box_change_domain; // 1 if Fix changes proc sub-domains, 0 if not bigint next_reneighbor; // next timestep to force a reneighboring int thermo_energy; // 1 if fix_modify enabled ThEng, 0 if not int nevery; // how often to call an end_of_step fix int rigid_flag; // 1 if Fix integrates rigid bodies, 0 if not int peatom_flag; // 1 if Fix contributes per-atom eng, 0 if not int virial_flag; // 1 if Fix contributes to virial, 0 if not int no_change_box; // 1 if cannot swap ortho <-> triclinic int time_integrate; // 1 if fix performs time integration, 0 if no int time_depend; // 1 if requires continuous timestepping int create_attribute; // 1 if fix stores attributes that need // setting when a new atom is created int restart_pbc; // 1 if fix moves atoms (except integrate) // so write_restart must remap to PBC int wd_header; // # of header values fix writes to data file int wd_section; // # of sections fix writes to data file int dynamic_group_allow; // 1 if can be used with dynamic group, else 0 int dof_flag; // 1 if has dof() method (not min_dof()) int special_alter_flag; // 1 if has special_alter() meth for spec lists int enforce2d_flag; // 1 if has enforce2d method int respa_level_support; // 1 if fix supports fix_modify respa int respa_level; // which respa level to apply fix (1-Nrespa) int scalar_flag; // 0/1 if compute_scalar() function exists int vector_flag; // 0/1 if compute_vector() function exists int array_flag; // 0/1 if compute_array() function exists int size_vector; // length of global vector int size_array_rows; // rows in global array int size_array_cols; // columns in global array int size_vector_variable; // 1 if vec length is unknown in advance int size_array_rows_variable; // 1 if array rows is unknown in advance int global_freq; // frequency s/v data is available at int peratom_flag; // 0/1 if per-atom data is stored int size_peratom_cols; // 0 = vector, N = columns in peratom array int peratom_freq; // frequency per-atom data is available at int local_flag; // 0/1 if local data is stored int size_local_rows; // rows in local vector or array int size_local_cols; // 0 = vector, N = columns in local array int local_freq; // frequency local data is available at int extscalar; // 0/1 if global scalar is intensive/extensive int extvector; // 0/1/-1 if global vector is all int/ext/extlist int *extlist; // list of 0/1 int/ext for each vec component int extarray; // 0/1 if global array is intensive/extensive double *vector_atom; // computed per-atom vector double **array_atom; // computed per-atom array double *vector_local; // computed local vector double **array_local; // computed local array int comm_forward; // size of forward communication (0 if none) int comm_reverse; // size of reverse communication (0 if none) int comm_border; // size of border communication (0 if none) - double virial[6]; // accumlated virial + double virial[6]; // accumulated virial double *eatom,**vatom; // accumulated per-atom energy/virial int restart_reset; // 1 if restart just re-initialized fix // KOKKOS host/device flag and data masks int kokkosable; // 1 if Kokkos fix ExecutionSpace execution_space; unsigned int datamask_read,datamask_modify; Fix(class LAMMPS *, int, char **); virtual ~Fix(); void modify_params(int, char **); virtual int setmask() = 0; virtual void post_constructor() {} virtual void init() {} virtual void init_list(int, class NeighList *) {} virtual void setup(int) {} virtual void setup_pre_exchange() {} virtual void setup_pre_neighbor() {} virtual void setup_pre_force(int) {} virtual void setup_pre_reverse(int, int) {} virtual void min_setup(int) {} virtual void initial_integrate(int) {} virtual void post_integrate() {} virtual void pre_exchange() {} virtual void pre_neighbor() {} virtual void pre_force(int) {} virtual void pre_reverse(int,int) {} virtual void post_force(int) {} virtual void final_integrate() {} virtual void end_of_step() {} virtual void post_run() {} virtual void write_restart(FILE *) {} virtual void write_restart_file(char *) {} virtual void restart(char *) {} virtual void grow_arrays(int) {} virtual void copy_arrays(int, int, int) {} virtual void set_arrays(int) {} virtual void update_arrays(int, int) {} virtual void set_molecule(int, tagint, int, double *, double *, double *) {} virtual void clear_bonus() {} virtual int pack_border(int, int *, double *) {return 0;} virtual int unpack_border(int, int, double *) {return 0;} virtual int pack_exchange(int, double *) {return 0;} virtual int unpack_exchange(int, double *) {return 0;} virtual int pack_restart(int, double *) {return 0;} virtual void unpack_restart(int, int) {} virtual int size_restart(int) {return 0;} virtual int maxsize_restart() {return 0;} virtual void setup_pre_force_respa(int, int) {} virtual void initial_integrate_respa(int, int, int) {} virtual void post_integrate_respa(int, int) {} virtual void pre_force_respa(int, int, int) {} virtual void post_force_respa(int, int, int) {} virtual void final_integrate_respa(int, int) {} virtual void min_pre_exchange() {} virtual void min_pre_neighbor() {} virtual void min_pre_force(int) {} virtual void min_pre_reverse(int,int) {} virtual void min_post_force(int) {} virtual double min_energy(double *) {return 0.0;} virtual void min_store() {} virtual void min_clearstore() {} virtual void min_pushstore() {} virtual void min_popstore() {} virtual int min_reset_ref() {return 0;} virtual void min_step(double, double *) {} virtual double max_alpha(double *) {return 0.0;} virtual int min_dof() {return 0;} virtual int pack_forward_comm(int, int *, double *, int, int *) {return 0;} virtual void unpack_forward_comm(int, int, double *) {} virtual int pack_reverse_comm_size(int, int) {return 0;} virtual int pack_reverse_comm(int, int, double *) {return 0;} virtual void unpack_reverse_comm(int, int *, double *) {} virtual double compute_scalar() {return 0.0;} virtual double compute_vector(int) {return 0.0;} virtual double compute_array(int,int) {return 0.0;} virtual int dof(int) {return 0;} virtual void deform(int) {} virtual void reset_target(double) {} virtual void reset_dt() {} virtual void enforce2d() {} virtual void read_data_header(char *) {} virtual void read_data_section(char *, int, char *, tagint) {} virtual bigint read_data_skip_lines(char *) {return 0;} virtual void write_data_header(FILE *, int) {} virtual void write_data_section_size(int, int &, int &) {} virtual void write_data_section_pack(int, double **) {} virtual void write_data_section_keyword(int, FILE *) {} virtual void write_data_section(int, FILE *, int, double **, int) {} virtual void zero_momentum() {} virtual void zero_rotation() {} virtual void rebuild_special() {} virtual int image(int *&, double **&) {return 0;} virtual int modify_param(int, char **) {return 0;} virtual void *extract(const char *, int &) {return NULL;} virtual double memory_usage() {return 0.0;} protected: int instance_me; // which Fix class instantiation I am int evflag; int eflag_either,eflag_global,eflag_atom; int vflag_either,vflag_global,vflag_atom; int maxeatom,maxvatom; int copymode; // if set, do not deallocate during destruction // required when classes are used as functors by Kokkos void ev_setup(int, int); void ev_tally(int, int *, double, double, double *); void v_setup(int); void v_tally(int, int *, double, double *); // union data struct for packing 32-bit and 64-bit ints into double bufs // see atom_vec.h for documentation union ubuf { double d; int64_t i; ubuf(double arg) : d(arg) {} ubuf(int64_t arg) : i(arg) {} ubuf(int arg) : i(arg) {} }; }; namespace FixConst { static const int INITIAL_INTEGRATE = 1<<0; static const int POST_INTEGRATE = 1<<1; static const int PRE_EXCHANGE = 1<<2; static const int PRE_NEIGHBOR = 1<<3; static const int PRE_FORCE = 1<<4; static const int PRE_REVERSE = 1<<5; static const int POST_FORCE = 1<<6; static const int FINAL_INTEGRATE = 1<<7; static const int END_OF_STEP = 1<<8; static const int POST_RUN = 1<<9; static const int THERMO_ENERGY = 1<<10; static const int INITIAL_INTEGRATE_RESPA = 1<<11; static const int POST_INTEGRATE_RESPA = 1<<12; static const int PRE_FORCE_RESPA = 1<<13; static const int POST_FORCE_RESPA = 1<<14; static const int FINAL_INTEGRATE_RESPA = 1<<15; static const int MIN_PRE_EXCHANGE = 1<<16; static const int MIN_PRE_NEIGHBOR = 1<<17; static const int MIN_PRE_FORCE = 1<<18; static const int MIN_PRE_REVERSE = 1<<19; static const int MIN_POST_FORCE = 1<<20; static const int MIN_ENERGY = 1<<21; static const int FIX_CONST_LAST = 1<<22; } } #endif /* ERROR/WARNING messages: E: Fix ID must be alphanumeric or underscore characters Self-explanatory. E: Could not find fix group ID A group ID used in the fix command does not exist. E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. */ diff --git a/src/fix_ave_chunk.h b/src/fix_ave_chunk.h index 18622699d..eb51043e0 100644 --- a/src/fix_ave_chunk.h +++ b/src/fix_ave_chunk.h @@ -1,185 +1,185 @@ /* -*- c++ -*- ---------------------------------------------------------- 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(ave/chunk,FixAveChunk) #else #ifndef LMP_FIX_AVE_CHUNK_H #define LMP_FIX_AVE_CHUNK_H #include #include "fix.h" namespace LAMMPS_NS { class FixAveChunk : public Fix { public: FixAveChunk(class LAMMPS *, int, char **); ~FixAveChunk(); int setmask(); void init(); void setup(int); void end_of_step(); double compute_array(int,int); double memory_usage(); private: int me,nvalues; int nrepeat,nfreq,irepeat; int normflag,scaleflag,overwrite,biasflag,colextra; bigint nvalid,nvalid_last; double adof,cdof; char *format,*format_user; char *tstring,*sstring,*id_bias; int *which,*argindex,*value2index; char **ids; class Compute *tbias; // ptr to additional bias compute FILE *fp; int ave,nwindow; int normcount,iwindow,window_limit; int nchunk,maxchunk; char *idchunk; class ComputeChunkAtom *cchunk; int lockforever; long filepos; int maxvar; double *varatom; // one,many,sum vecs/arrays are used with a single Nfreq epoch // total,list vecs/arrays are used across epochs double *count_one,*count_many,*count_sum; double **values_one,**values_many,**values_sum; double *count_total,**count_list; double **values_total,***values_list; void allocate(); bigint nextvalid(); }; } #endif #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: No values in fix ave/chunk command Self-explanatory. E: Cannot open fix ave/chunk file %s The specified file cannot be opened. Check that the path and name are correct. E: Could not find compute ID for temperature bias Self-explanatory. E: Bias compute does not calculate temperature The specified compute must compute temperature. E: Bias compute does not calculate a velocity bias The specified compute must compute a bias for temperature. E: Compute ID for fix ave/chunk does not exist Self-explanatory. E: Fix ave/chunk compute does not calculate per-atom values Self-explanatory. E: Fix ave/chunk compute does not calculate a per-atom vector Self-explanatory. E: Fix ave/chunk compute does not calculate a per-atom array Self-explanatory. E: Fix ave/chunk compute vector is accessed out-of-range Self-explanatory. E: Fix ID for fix ave/chunk does not exist Self-explanatory. E: Fix ave/chunk fix does not calculate per-atom values Self-explanatory. E: Fix ave/chunk fix does not calculate a per-atom vector Self-explanatory. E: Fix ave/chunk fix does not calculate a per-atom array Self-explanatory. E: Fix ave/chunk fix vector is accessed out-of-range Self-explanatory. E: Variable name for fix ave/chunk does not exist Self-explanatory. E: Fix ave/chunk variable is not atom-style variable Self-explanatory. E: Chunk/atom compute does not exist for fix ave/chunk Self-explanatory. E: Fix ave/chunk does not use chunk/atom compute -The specified conpute is not for a compute chunk/atom command. +The specified compute is not for a compute chunk/atom command. E: Error writing file header Something in the output to the file triggered an error. E: Fix for fix ave/chunk not computed at compatible time Fixes generate their values on specific timesteps. Fix ave/chunk is requesting a value on a non-allowed timestep. E: Invalid timestep reset for fix ave/chunk Resetting the timestep has invalidated the sequence of timesteps this fix needs to process. E: Error writing averaged chunk data Something in the output to the file triggered an error. */ diff --git a/src/fix_langevin.h b/src/fix_langevin.h index 863c0f554..34eec6398 100644 --- a/src/fix_langevin.h +++ b/src/fix_langevin.h @@ -1,160 +1,160 @@ /* -*- c++ -*- ---------------------------------------------------------- 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(langevin,FixLangevin) #else #ifndef LMP_FIX_LANGEVIN_H #define LMP_FIX_LANGEVIN_H #include "fix.h" namespace LAMMPS_NS { class FixLangevin : public Fix { public: FixLangevin(class LAMMPS *, int, char **); virtual ~FixLangevin(); int setmask(); void init(); void setup(int); virtual void post_force(int); void post_force_respa(int, int, int); virtual void end_of_step(); void reset_target(double); void reset_dt(); int modify_param(int, char **); virtual double compute_scalar(); double memory_usage(); virtual void *extract(const char *, int &); void grow_arrays(int); void copy_arrays(int, int, int); int pack_exchange(int, double *); int unpack_exchange(int, double *); protected: int gjfflag,oflag,tallyflag,zeroflag,tbiasflag; int flangevin_allocated; double ascale; double t_start,t_stop,t_period,t_target; double *gfactor1,*gfactor2,*ratio; double energy,energy_onestep; double tsqrt; int tstyle,tvar; double gjffac; char *tstr; class AtomVecEllipsoid *avec; int maxatom1,maxatom2; double **flangevin; double *tforce; double **franprev; int nvalues; char *id_temp; class Compute *temperature; int nlevels_respa; class RanMars *random; int seed; // comment next line to turn off templating #define TEMPLATED_FIX_LANGEVIN #ifdef TEMPLATED_FIX_LANGEVIN template < int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY, int Tp_BIAS, int Tp_RMASS, int Tp_ZERO > void post_force_templated(); #else void post_force_untemplated(int, int, int, int, int, int); #endif void omega_thermostat(); void angmom_thermostat(); void compute_target(); }; } #endif #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Fix langevin period must be > 0.0 The time window for temperature relaxation must be > 0 W: Energy tally does not account for 'zero yes' The energy removed by using the 'zero yes' flag is not accounted for in the energy tally and thus energy conservation cannot be monitored in this case. E: Fix langevin omega requires atom style sphere Self-explanatory. E: Fix langevin angmom requires atom style ellipsoid Self-explanatory. E: Variable name for fix langevin does not exist Self-explanatory. E: Variable for fix langevin is invalid style It must be an equal-style variable. E: Fix langevin omega requires extended particles One of the particles has radius 0.0. E: Fix langevin angmom requires extended particles -This fix option cannot be used with point paritlces. +This fix option cannot be used with point particles. E: Cannot zero Langevin force of 0 atoms The group has zero atoms, so you cannot request its force be zeroed. E: Fix langevin variable returned negative temperature Self-explanatory. E: Could not find fix_modify temperature ID The compute ID for computing temperature does not exist. E: Fix_modify temperature ID does not compute temperature The compute ID assigned to the fix must compute temperature. W: Group for fix_modify temp != fix group The fix_modify command is specifying a temperature computation that computes a temperature on a different group of atoms than the fix itself operates on. This is probably not what you want to do. */ diff --git a/src/improper.h b/src/improper.h index b20bc732d..be78b6acf 100644 --- a/src/improper.h +++ b/src/improper.h @@ -1,81 +1,81 @@ /* -*- c++ -*- ---------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifndef LMP_IMPROPER_H #define LMP_IMPROPER_H #include #include "pointers.h" namespace LAMMPS_NS { class Improper : protected Pointers { friend class ThrOMP; friend class FixOMP; public: int allocated; int *setflag; int writedata; // 1 if writes coeffs to data file double energy; // accumulated energies - double virial[6]; // accumlated virial + double virial[6]; // accumulated virial double *eatom,**vatom; // accumulated per-atom energy/virial // KOKKOS host/device flag and data masks ExecutionSpace execution_space; unsigned int datamask_read,datamask_modify; int copymode; Improper(class LAMMPS *); virtual ~Improper(); virtual void init(); virtual void init_style() {} virtual void compute(int, int) = 0; virtual void settings(int, char **) {} virtual void coeff(int, char **) = 0; virtual void write_restart(FILE *) = 0; virtual void read_restart(FILE *) = 0; virtual void write_data(FILE *) {} virtual double memory_usage(); protected: int suffix_flag; // suffix compatibility flag int evflag; int eflag_either,eflag_global,eflag_atom; int vflag_either,vflag_global,vflag_atom; int maxeatom,maxvatom; void ev_setup(int, int); void ev_tally(int, int, int, int, int, int, double, double *, double *, double *, double, double, double, double, double, double, double, double, double); }; } #endif /* ERROR/WARNING messages: E: Improper coeffs are not set No improper coefficients have been assigned in the data file or via the improper_coeff command. E: All improper coeffs are not set All improper coefficients must be set in the data file or by the improper_coeff command before running a simulation. */ diff --git a/src/kspace.h b/src/kspace.h index 95fb6ffaf..c51451619 100644 --- a/src/kspace.h +++ b/src/kspace.h @@ -1,266 +1,266 @@ /* -*- c++ -*- ---------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifndef LMP_KSPACE_H #define LMP_KSPACE_H #include "pointers.h" #include "accelerator_kokkos.h" #ifdef FFT_SINGLE typedef float FFT_SCALAR; #define MPI_FFT_SCALAR MPI_FLOAT #else typedef double FFT_SCALAR; #define MPI_FFT_SCALAR MPI_DOUBLE #endif namespace LAMMPS_NS { class KSpace : protected Pointers { friend class ThrOMP; friend class FixOMP; public: double energy; // accumulated energies double energy_1,energy_6; - double virial[6]; // accumlated virial + double virial[6]; // accumulated virial double *eatom,**vatom; // accumulated per-atom energy/virial double e2group; // accumulated group-group energy double f2group[3]; // accumulated group-group force int triclinic_support; // 1 if supports triclinic geometries int ewaldflag; // 1 if a Ewald solver int pppmflag; // 1 if a PPPM solver int msmflag; // 1 if a MSM solver int dispersionflag; // 1 if a LJ/dispersion solver int tip4pflag; // 1 if a TIP4P solver int dipoleflag; // 1 if a dipole solver int differentiation_flag; int neighrequest_flag; // used to avoid obsolete construction // of neighbor lists int mixflag; // 1 if geometric mixing rules are enforced // for LJ coefficients int slabflag; int scalar_pressure_flag; // 1 if using MSM fast scalar pressure double slab_volfactor; int warn_nonneutral; // 0 = error if non-neutral system // 1 = warn once if non-neutral system // 2 = warn, but already warned int warn_nocharge; // 0 = already warned // 1 = warn if zero charge int order,order_6,order_allocated; double accuracy; // accuracy of KSpace solver (force units) double accuracy_absolute; // user-specifed accuracy in force units double accuracy_relative; // user-specified dimensionless accuracy // accurary = acc_rel * two_charge_force double accuracy_real_6; // real space accuracy for // dispersion solver (force units) double accuracy_kspace_6; // reciprocal space accuracy for // dispersion solver (force units) int auto_disp_flag; // use automatic paramter generation for pppm/disp double two_charge_force; // force in user units of two point // charges separated by 1 Angstrom double g_ewald,g_ewald_6; int nx_pppm,ny_pppm,nz_pppm; // global FFT grid for Coulombics int nx_pppm_6,ny_pppm_6,nz_pppm_6; // global FFT grid for dispersion int nx_msm_max,ny_msm_max,nz_msm_max; int group_group_enable; // 1 if style supports group/group calculation // KOKKOS host/device flag and data masks ExecutionSpace execution_space; unsigned int datamask_read,datamask_modify; int copymode; int compute_flag; // 0 if skip compute() int fftbench; // 0 if skip FFT timing int collective_flag; // 1 if use MPI collectives for FFT/remap int stagger_flag; // 1 if using staggered PPPM grids double splittol; // tolerance for when to truncate splitting KSpace(class LAMMPS *, int, char **); virtual ~KSpace(); void triclinic_check(); void modify_params(int, char **); void *extract(const char *); void compute_dummy(int, int); // triclinic void x2lamdaT(double *, double *); void lamda2xT(double *, double *); void lamda2xvector(double *, double *); void kspacebbox(double, double *); // public so can be called by commands that change charge void qsum_qsq(); // general child-class methods virtual void init() = 0; virtual void setup() = 0; virtual void setup_grid() {}; virtual void compute(int, int) = 0; virtual void compute_group_group(int, int, int) {}; virtual void pack_forward(int, FFT_SCALAR *, int, int *) {}; virtual void unpack_forward(int, FFT_SCALAR *, int, int *) {}; virtual void pack_reverse(int, FFT_SCALAR *, int, int *) {}; virtual void unpack_reverse(int, FFT_SCALAR *, int, int *) {}; virtual void pack_forward_kokkos(int, DAT::tdual_FFT_SCALAR_1d &, int, DAT::tdual_int_2d &, int) {}; virtual void unpack_forward_kokkos(int, DAT::tdual_FFT_SCALAR_1d &, int, DAT::tdual_int_2d &, int) {}; virtual void pack_reverse_kokkos(int, DAT::tdual_FFT_SCALAR_1d &, int, DAT::tdual_int_2d &, int) {}; virtual void unpack_reverse_kokkos(int, DAT::tdual_FFT_SCALAR_1d &, int, DAT::tdual_int_2d &, int) {}; virtual int timing(int, double &, double &) {return 0;} virtual int timing_1d(int, double &) {return 0;} virtual int timing_3d(int, double &) {return 0;} virtual double memory_usage() {return 0.0;} /* ---------------------------------------------------------------------- compute gamma for MSM and pair styles see Eq 4 from Parallel Computing 35 (2009) 164–177 ------------------------------------------------------------------------- */ double gamma(const double &rho) const { if (rho <= 1.0) { const int split_order = order/2; const double rho2 = rho*rho; double g = gcons[split_order][0]; double rho_n = rho2; for (int n = 1; n <= split_order; n++) { g += gcons[split_order][n]*rho_n; rho_n *= rho2; } return g; } else return (1.0/rho); } /* ---------------------------------------------------------------------- compute the derivative of gamma for MSM and pair styles see Eq 4 from Parallel Computing 35 (2009) 164-177 ------------------------------------------------------------------------- */ double dgamma(const double &rho) const { if (rho <= 1.0) { const int split_order = order/2; const double rho2 = rho*rho; double dg = dgcons[split_order][0]*rho; double rho_n = rho*rho2; for (int n = 1; n < split_order; n++) { dg += dgcons[split_order][n]*rho_n; rho_n *= rho2; } return dg; } else return (-1.0/rho/rho); } double **get_gcons() { return gcons; } double **get_dgcons() { return dgcons; } protected: int gridflag,gridflag_6; int gewaldflag,gewaldflag_6; int minorder,overlap_allowed; int adjust_cutoff_flag; int suffix_flag; // suffix compatibility flag bigint natoms_original; double scale,qqrd2e; double qsum,qsqsum,q2; double **gcons,**dgcons; // accumulated per-atom energy/virial int evflag,evflag_atom; int eflag_either,eflag_global,eflag_atom; int vflag_either,vflag_global,vflag_atom; int maxeatom,maxvatom; int kewaldflag; // 1 if kspace range set for Ewald sum int kx_ewald,ky_ewald,kz_ewald; // kspace settings for Ewald sum void pair_check(); void ev_setup(int, int); double estimate_table_accuracy(double, double); }; } #endif /* ERROR/WARNING messages: E: KSpace style does not yet support triclinic geometries The specified kspace style does not allow for non-orthogonal simulation boxes. E: KSpace solver requires a pair style No pair style is defined. E: KSpace style is incompatible with Pair style Setting a kspace style requires that a pair style with matching long-range Coulombic or dispersion components be used. W: Using kspace solver on system with no charge Self-explanatory. E: System is not charge neutral, net charge = %g The total charge on all atoms on the system is not 0.0. For some KSpace solvers this is an error. W: System is not charge neutral, net charge = %g The total charge on all atoms on the system is not 0.0. For some KSpace solvers this is only a warning. W: For better accuracy use 'pair_modify table 0' The user-specified force accuracy cannot be achieved unless the table feature is disabled by using 'pair_modify table 0'. E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Bad kspace_modify slab parameter Kspace_modify value for the slab/volume keyword must be >= 2.0. W: Kspace_modify slab param < 2.0 may cause unphysical behavior The kspace_modify slab parameter should be larger to insure periodic grids padded with empty space do not overlap. E: Bad kspace_modify kmax/ewald parameter Kspace_modify values for the kmax/ewald keyword must be integers > 0 E: Kspace_modify eigtol must be smaller than one Self-explanatory. */ diff --git a/src/procmap.cpp b/src/procmap.cpp index 088939b21..f54ac1915 100644 --- a/src/procmap.cpp +++ b/src/procmap.cpp @@ -1,895 +1,895 @@ /* ---------------------------------------------------------------------- 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 author (NUMA option) : Mike Brown (ORNL) ------------------------------------------------------------------------- */ #include "procmap.h" #include "universe.h" #include "domain.h" #include "math_extra.h" #include "memory.h" #include "error.h" #include #include using namespace LAMMPS_NS; #define MAXLINE 128 enum{MULTIPLE}; // same as in Comm /* ---------------------------------------------------------------------- */ ProcMap::ProcMap(LAMMPS *lmp) : Pointers(lmp) {} /* ---------------------------------------------------------------------- create a one-level 3d grid of procs ------------------------------------------------------------------------- */ void ProcMap::onelevel_grid(int nprocs, int *user_procgrid, int *procgrid, int otherflag, int other_style, int *other_procgrid, int *other_coregrid) { int **factors; // factors = list of all possible 3 factors of processor count int npossible = factor(nprocs,NULL); memory->create(factors,npossible,3,"procmap:factors"); npossible = factor(nprocs,factors); // constrain by 2d, user request, other partition if (domain->dimension == 2) npossible = cull_2d(npossible,factors,3); npossible = cull_user(npossible,factors,3,user_procgrid); if (otherflag) npossible = cull_other(npossible,factors,3, other_style,other_procgrid, other_coregrid); // user/other constraints make failure possible if (npossible == 0) error->all(FLERR,"Could not create 3d grid of processors"); // select best set of 3 factors based on surface area of proc sub-domains best_factors(npossible,factors,procgrid,1,1,1); // clean-up memory->destroy(factors); } /* ---------------------------------------------------------------------- create a two-level 3d grid of procs ------------------------------------------------------------------------- */ void ProcMap::twolevel_grid(int nprocs, int *user_procgrid, int *procgrid, int ncores, int *user_coregrid, int *coregrid, int otherflag, int other_style, int *other_procgrid, int *other_coregrid) { int **nfactors,**cfactors,**factors; if (nprocs % ncores) error->all(FLERR,"Processors twogrid requires proc count " "be a multiple of core count"); // nfactors = list of all possible 3 factors of node count // constrain by 2d int nnpossible = factor(nprocs/ncores,NULL); memory->create(nfactors,nnpossible,3,"procmap:nfactors"); nnpossible = factor(nprocs/ncores,nfactors); if (domain->dimension == 2) nnpossible = cull_2d(nnpossible,nfactors,3); // cfactors = list of all possible 3 factors of core count // constrain by 2d int ncpossible = factor(ncores,NULL); memory->create(cfactors,ncpossible,3,"procmap:cfactors"); ncpossible = factor(ncores,cfactors); if (domain->dimension == 2) ncpossible = cull_2d(ncpossible,cfactors,3); ncpossible = cull_user(ncpossible,cfactors,3,user_coregrid); // factors = all combinations of nfactors and cfactors // factors stores additional index pointing to corresponding cfactors // constrain by user request, other partition int npossible = nnpossible * ncpossible; memory->create(factors,npossible,4,"procmap:factors"); npossible = combine_factors(nnpossible,nfactors,ncpossible,cfactors,factors); npossible = cull_user(npossible,factors,4,user_procgrid); if (otherflag) npossible = cull_other(npossible,factors,4, other_style,other_procgrid, other_coregrid); // user/other constraints make failure possible if (npossible == 0) error->all(FLERR,"Could not create twolevel 3d grid of processors"); // select best set of 3 factors based on surface area of proc sub-domains // index points to corresponding core factorization int index = best_factors(npossible,factors,procgrid,1,1,1); coregrid[0] = cfactors[factors[index][3]][0]; coregrid[1] = cfactors[factors[index][3]][1]; coregrid[2] = cfactors[factors[index][3]][2]; // clean-up memory->destroy(nfactors); memory->destroy(cfactors); memory->destroy(factors); } /* ---------------------------------------------------------------------- create a 3d grid of procs that does a 2-level hierarchy within a node auto-detects NUMA sockets within a multi-core node ------------------------------------------------------------------------- */ void ProcMap::numa_grid(int nprocs, int *user_procgrid, int *procgrid, int *numagrid) { // hardwire this for now int numa_nodes = 1; // get names of all nodes int name_length; char node_name[MPI_MAX_PROCESSOR_NAME]; MPI_Get_processor_name(node_name,&name_length); node_name[name_length] = '\0'; char *node_names = new char[MPI_MAX_PROCESSOR_NAME*nprocs]; MPI_Allgather(node_name,MPI_MAX_PROCESSOR_NAME,MPI_CHAR,node_names, MPI_MAX_PROCESSOR_NAME,MPI_CHAR,world); std::string node_string = std::string(node_name); // get number of procs per node // NOTE: could do this without STL map std::map name_map; std::map::iterator np; for (int i = 0; i < nprocs; i++) { std::string i_string = std::string(&node_names[i*MPI_MAX_PROCESSOR_NAME]); np = name_map.find(i_string); if (np == name_map.end()) name_map[i_string] = 1; else np->second++; } procs_per_node = name_map.begin()->second; procs_per_numa = procs_per_node / numa_nodes; delete [] node_names; // error if any of these conditions met if (nprocs % procs_per_numa || // total procs not a multiple of node user_procgrid[0] > 1 || // user specified grid > 1 in any dim user_procgrid[1] > 1 || user_procgrid[2] > 1) error->all(FLERR,"Could not create numa grid of processors"); // user settings for the factorization per numa node // currently not user settable // if user specifies 1 for a proc grid dimension, // also use 1 for the numa grid dimension int user_numagrid[3]; user_numagrid[0] = user_numagrid[1] = user_numagrid[2] = 0; if (user_procgrid[0] == 1) user_numagrid[0] = 1; if (user_procgrid[1] == 1) user_numagrid[1] = 1; if (user_procgrid[2] == 1) user_numagrid[2] = 1; // initial factorization within NUMA node int **numafactors; int numapossible = factor(procs_per_numa,NULL); memory->create(numafactors,numapossible,3,"procmap:numafactors"); numapossible = factor(procs_per_numa,numafactors); if (domain->dimension == 2) numapossible = cull_2d(numapossible,numafactors,3); numapossible = cull_user(numapossible,numafactors,3,user_numagrid); if (numapossible == 0) error->all(FLERR,"Could not create numa grid of processors"); best_factors(numapossible,numafactors,numagrid,1,1,1); - // user_nodegrid = implied user contraints on nodes + // user_nodegrid = implied user constraints on nodes int user_nodegrid[3]; user_nodegrid[0] = user_procgrid[0] / numagrid[0]; user_nodegrid[1] = user_procgrid[1] / numagrid[1]; user_nodegrid[2] = user_procgrid[2] / numagrid[2]; // factorization for the grid of NUMA nodes int node_count = nprocs / procs_per_numa; int **nodefactors; int nodepossible = factor(node_count,NULL); memory->create(nodefactors,nodepossible,3,"procmap:nodefactors"); nodepossible = factor(node_count,nodefactors); if (domain->dimension == 2) nodepossible = cull_2d(nodepossible,nodefactors,3); nodepossible = cull_user(nodepossible,nodefactors,3,user_nodegrid); if (nodepossible == 0) error->all(FLERR,"Could not create numa grid of processors"); best_factors(nodepossible,nodefactors,nodegrid, numagrid[0],numagrid[1],numagrid[2]); // repeat NUMA node factorization using subdomain sizes // refines the factorization if the user specified the node layout // NOTE: this will not re-enforce user-procgrid constraint will it? best_factors(numapossible,numafactors,numagrid, nodegrid[0],nodegrid[1],nodegrid[2]); memory->destroy(numafactors); memory->destroy(nodefactors); // assign a unique id to each node node_id = 0; int node_num = 0; for (np = name_map.begin(); np != name_map.end(); ++np) { if (np->first == node_string) node_id = node_num; node_num++; } // return the proc-level factorization procgrid[0] = nodegrid[0] * numagrid[0]; procgrid[1] = nodegrid[1] * numagrid[1]; procgrid[2] = nodegrid[2] * numagrid[2]; } /* ---------------------------------------------------------------------- define a 3d grid from a custom file ------------------------------------------------------------------------- */ void ProcMap::custom_grid(char *cfile, int nprocs, int *user_procgrid, int *procgrid) { int me; MPI_Comm_rank(world,&me); char line[MAXLINE]; FILE *fp = NULL; if (me == 0) { fp = fopen(cfile,"r"); if (fp == NULL) error->one(FLERR,"Cannot open custom file"); // skip header = blank and comment lines char *ptr; if (!fgets(line,MAXLINE,fp)) error->one(FLERR,"Unexpected end of custom file"); while (1) { if ((ptr = strchr(line,'#'))) *ptr = '\0'; if (strspn(line," \t\n\r") != strlen(line)) break; if (!fgets(line,MAXLINE,fp)) error->one(FLERR,"Unexpected end of custom file"); } } int n = strlen(line) + 1; MPI_Bcast(&n,1,MPI_INT,0,world); MPI_Bcast(line,n,MPI_CHAR,0,world); sscanf(line,"%d %d %d",&procgrid[0],&procgrid[1],&procgrid[2]); int flag = 0; if (procgrid[0]*procgrid[1]*procgrid[2] != nprocs) flag = 1; if (user_procgrid[0] && procgrid[0] != user_procgrid[0]) flag = 1; if (user_procgrid[1] && procgrid[1] != user_procgrid[1]) flag = 1; if (user_procgrid[2] && procgrid[2] != user_procgrid[2]) flag = 1; if (flag) error->all(FLERR,"Processors custom grid file is inconsistent"); // cmap = map of procs to grid // store for use in custom_map() memory->create(cmap,nprocs,4,"procmap:cmap"); for (int i = 0; i < nprocs; i++) cmap[i][0] = -1; if (me == 0) { for (int i = 0; i < nprocs; i++) { if (!fgets(line,MAXLINE,fp)) error->one(FLERR,"Unexpected end of custom file"); sscanf(line,"%d %d %d %d", &cmap[i][0],&cmap[i][1],&cmap[i][2],&cmap[i][3]); } fclose(fp); } MPI_Bcast(&cmap[0][0],nprocs*4,MPI_INT,0,world); // error check on cmap values flag = 0; for (int i = 0; i < nprocs; i++) { if (cmap[i][0] == -1) flag = 1; else { if (cmap[i][1] <= 0 || cmap[i][1] > procgrid[0]) flag = 1; if (cmap[i][2] <= 0 || cmap[i][2] > procgrid[1]) flag = 1; if (cmap[i][3] <= 0 || cmap[i][3] > procgrid[2]) flag = 1; } } if (flag) error->all(FLERR,"Processors custom grid file is inconsistent"); } /* ---------------------------------------------------------------------- map processors to 3d grid via MPI_Cart routines MPI may do layout in machine-optimized fashion ------------------------------------------------------------------------- */ void ProcMap::cart_map(int reorder, int *procgrid, int *myloc, int procneigh[3][2], int ***grid2proc) { int periods[3]; periods[0] = periods[1] = periods[2] = 1; MPI_Comm cartesian; MPI_Cart_create(world,3,procgrid,periods,reorder,&cartesian); MPI_Cart_get(cartesian,3,procgrid,periods,myloc); MPI_Cart_shift(cartesian,0,1,&procneigh[0][0],&procneigh[0][1]); MPI_Cart_shift(cartesian,1,1,&procneigh[1][0],&procneigh[1][1]); MPI_Cart_shift(cartesian,2,1,&procneigh[2][0],&procneigh[2][1]); int coords[3]; int i,j,k; for (i = 0; i < procgrid[0]; i++) for (j = 0; j < procgrid[1]; j++) for (k = 0; k < procgrid[2]; k++) { coords[0] = i; coords[1] = j; coords[2] = k; MPI_Cart_rank(cartesian,coords,&grid2proc[i][j][k]); } MPI_Comm_free(&cartesian); } /* ---------------------------------------------------------------------- map processors to 3d grid via MPI_Cart routines respect sub-grid of cores within each node MPI may do layout in machine-optimized fashion ------------------------------------------------------------------------- */ void ProcMap::cart_map(int reorder, int *procgrid, int ncores, int *coregrid, int *myloc, int procneigh[3][2], int ***grid2proc) { // setup NUMA params that numa_grid() sets up int me; MPI_Comm_rank(world,&me); procs_per_node = ncores; procs_per_numa = ncores; node_id = me/ncores; nodegrid[0] = procgrid[0] / coregrid[0]; nodegrid[1] = procgrid[1] / coregrid[1]; nodegrid[2] = procgrid[2] / coregrid[2]; // now can use numa_map() to perform mapping numa_map(reorder,coregrid,myloc,procneigh,grid2proc); } /* ---------------------------------------------------------------------- map processors to 3d grid in XYZ order called by onelevel ------------------------------------------------------------------------- */ void ProcMap::xyz_map(char *xyz, int *procgrid, int *myloc, int procneigh[3][2], int ***grid2proc) { int me; MPI_Comm_rank(world,&me); int i,j,k; for (i = 0; i < procgrid[0]; i++) for (j = 0; j < procgrid[1]; j++) for (k = 0; k < procgrid[2]; k++) { if (xyz[0] == 'x' && xyz[1] == 'y' && xyz[2] == 'z') grid2proc[i][j][k] = k*procgrid[1]*procgrid[0] + j*procgrid[0] + i; else if (xyz[0] == 'x' && xyz[1] == 'z' && xyz[2] == 'y') grid2proc[i][j][k] = j*procgrid[2]*procgrid[0] + k*procgrid[0] + i; else if (xyz[0] == 'y' && xyz[1] == 'x' && xyz[2] == 'z') grid2proc[i][j][k] = k*procgrid[0]*procgrid[1] + i*procgrid[1] + j; else if (xyz[0] == 'y' && xyz[1] == 'z' && xyz[2] == 'x') grid2proc[i][j][k] = i*procgrid[2]*procgrid[1] + k*procgrid[1] + j; else if (xyz[0] == 'z' && xyz[1] == 'x' && xyz[2] == 'y') grid2proc[i][j][k] = j*procgrid[0]*procgrid[2] + i*procgrid[2] + k; else if (xyz[0] == 'z' && xyz[1] == 'y' && xyz[2] == 'x') grid2proc[i][j][k] = i*procgrid[1]*procgrid[2] + j*procgrid[2] + k; if (grid2proc[i][j][k] == me) { myloc[0] = i; myloc[1] = j, myloc[2] = k; } } // proc IDs of neighbors int minus,plus; grid_shift(myloc[0],procgrid[0],minus,plus); procneigh[0][0] = grid2proc[minus][myloc[1]][myloc[2]]; procneigh[0][1] = grid2proc[plus][myloc[1]][myloc[2]]; grid_shift(myloc[1],procgrid[1],minus,plus); procneigh[1][0] = grid2proc[myloc[0]][minus][myloc[2]]; procneigh[1][1] = grid2proc[myloc[0]][plus][myloc[2]]; grid_shift(myloc[2],procgrid[2],minus,plus); procneigh[2][0] = grid2proc[myloc[0]][myloc[1]][minus]; procneigh[2][1] = grid2proc[myloc[0]][myloc[1]][plus]; } /* ---------------------------------------------------------------------- map processors to 3d grid in XYZ order respect sub-grid of cores within each node called by twolevel ------------------------------------------------------------------------- */ void ProcMap::xyz_map(char *xyz, int *procgrid, int ncores, int *coregrid, int *myloc, int procneigh[3][2], int ***grid2proc) { int me; MPI_Comm_rank(world,&me); nodegrid[0] = procgrid[0] / coregrid[0]; nodegrid[1] = procgrid[1] / coregrid[1]; nodegrid[2] = procgrid[2] / coregrid[2]; int i,j,k,inode,jnode,knode,icore,jcore,kcore; for (i = 0; i < procgrid[0]; i++) for (j = 0; j < procgrid[1]; j++) for (k = 0; k < procgrid[2]; k++) { inode = i/coregrid[0]; jnode = j/coregrid[1]; knode = k/coregrid[2]; icore = i % coregrid[0]; jcore = j % coregrid[1]; kcore = k % coregrid[2]; if (xyz[0] == 'x' && xyz[1] == 'y' && xyz[2] == 'z') { grid2proc[i][j][k] = ncores * (knode*nodegrid[1]*nodegrid[0] + jnode*nodegrid[0] + inode) + (kcore*coregrid[1]*coregrid[0] + jcore*coregrid[0] + icore); } else if (xyz[0] == 'x' && xyz[1] == 'z' && xyz[2] == 'y') grid2proc[i][j][k] = ncores * (jnode*nodegrid[2]*nodegrid[0] + knode*nodegrid[0] + inode) + (jcore*coregrid[2]*coregrid[0] + kcore*coregrid[0] + icore); else if (xyz[0] == 'y' && xyz[1] == 'x' && xyz[2] == 'z') grid2proc[i][j][k] = ncores * (knode*nodegrid[0]*nodegrid[1] + inode*nodegrid[1] + jnode) + (kcore*coregrid[0]*coregrid[1] + icore*coregrid[1] + jcore); else if (xyz[0] == 'y' && xyz[1] == 'z' && xyz[2] == 'x') grid2proc[i][j][k] = ncores * (inode*nodegrid[2]*nodegrid[1] + knode*nodegrid[1] + jnode) + (icore*coregrid[2]*coregrid[1] + kcore*coregrid[1] + jcore); else if (xyz[0] == 'z' && xyz[1] == 'x' && xyz[2] == 'y') grid2proc[i][j][k] = ncores * (jnode*nodegrid[0]*nodegrid[2] + inode*nodegrid[2] + knode) + (jcore*coregrid[0]*coregrid[2] + icore*coregrid[2] + kcore); else if (xyz[0] == 'z' && xyz[1] == 'y' && xyz[2] == 'x') grid2proc[i][j][k] = ncores * (inode*nodegrid[1]*nodegrid[2] + jnode*nodegrid[2] + knode) + (icore*coregrid[1]*coregrid[2] + jcore*coregrid[2] + kcore); if (grid2proc[i][j][k] == me) { myloc[0] = i; myloc[1] = j, myloc[2] = k; } } // proc IDs of neighbors int minus,plus; grid_shift(myloc[0],procgrid[0],minus,plus); procneigh[0][0] = grid2proc[minus][myloc[1]][myloc[2]]; procneigh[0][1] = grid2proc[plus][myloc[1]][myloc[2]]; grid_shift(myloc[1],procgrid[1],minus,plus); procneigh[1][0] = grid2proc[myloc[0]][minus][myloc[2]]; procneigh[1][1] = grid2proc[myloc[0]][plus][myloc[2]]; grid_shift(myloc[2],procgrid[2],minus,plus); procneigh[2][0] = grid2proc[myloc[0]][myloc[1]][minus]; procneigh[2][1] = grid2proc[myloc[0]][myloc[1]][plus]; } /* ---------------------------------------------------------------------- map processors to 3d grid in 2-level NUMA ordering ------------------------------------------------------------------------- */ void ProcMap::numa_map(int reorder, int *numagrid, int *myloc, int procneigh[3][2], int ***grid2proc) { // setup a per node communicator and find rank within MPI_Comm node_comm; MPI_Comm_split(world,node_id,0,&node_comm); int node_rank; MPI_Comm_rank(node_comm,&node_rank); // setup a per numa communicator and find rank within MPI_Comm numa_comm; int local_numa = node_rank / procs_per_numa; MPI_Comm_split(node_comm,local_numa,0,&numa_comm); int numa_rank; MPI_Comm_rank(numa_comm,&numa_rank); // setup a communicator with the rank 0 procs from each numa node MPI_Comm numa_leaders; MPI_Comm_split(world,numa_rank,0,&numa_leaders); // use the MPI Cartesian routines to map the nodes to the grid int periods[3]; periods[0] = periods[1] = periods[2] = 1; MPI_Comm cartesian; if (numa_rank == 0) { MPI_Cart_create(numa_leaders,3,nodegrid,periods,reorder,&cartesian); MPI_Cart_get(cartesian,3,nodegrid,periods,myloc); } // broadcast numa node location in grid to other procs in numa node MPI_Bcast(myloc,3,MPI_INT,0,numa_comm); // compute my location within the node grid int z_offset = numa_rank / (numagrid[0] * numagrid[1]); int y_offset = (numa_rank % (numagrid[0] * numagrid[1]))/numagrid[0]; int x_offset = numa_rank % numagrid[0]; myloc[0] = myloc[0] * numagrid[0] + x_offset; myloc[1] = myloc[1] * numagrid[1] + y_offset; myloc[2] = myloc[2] * numagrid[2] + z_offset; // allgather of myloc into gridi to fill grid2proc int nprocs; MPI_Comm_size(world,&nprocs); int **gridi; memory->create(gridi,nprocs,3,"comm:gridi"); MPI_Allgather(myloc,3,MPI_INT,gridi[0],3,MPI_INT,world); for (int i = 0; i < nprocs; i++) grid2proc[gridi[i][0]][gridi[i][1]][gridi[i][2]] = i; memory->destroy(gridi); // proc IDs of neighbors int minus,plus; grid_shift(myloc[0],nodegrid[0]*numagrid[0],minus,plus); procneigh[0][0] = grid2proc[minus][myloc[1]][myloc[2]]; procneigh[0][1] = grid2proc[plus][myloc[1]][myloc[2]]; grid_shift(myloc[1],nodegrid[1]*numagrid[1],minus,plus); procneigh[1][0] = grid2proc[myloc[0]][minus][myloc[2]]; procneigh[1][1] = grid2proc[myloc[0]][plus][myloc[2]]; grid_shift(myloc[2],nodegrid[2]*numagrid[2],minus,plus); procneigh[2][0] = grid2proc[myloc[0]][myloc[1]][minus]; procneigh[2][1] = grid2proc[myloc[0]][myloc[1]][plus]; // clean-up if (numa_rank == 0) MPI_Comm_free(&cartesian); MPI_Comm_free(&numa_leaders); MPI_Comm_free(&numa_comm); MPI_Comm_free(&node_comm); } /* ---------------------------------------------------------------------- map processors to 3d grid in custom ordering ------------------------------------------------------------------------- */ void ProcMap::custom_map(int *procgrid, int *myloc, int procneigh[3][2], int ***grid2proc) { int me,nprocs; MPI_Comm_rank(world,&me); MPI_Comm_size(world,&nprocs); for (int i = 0; i < nprocs; i++) { grid2proc[cmap[i][1]-1][cmap[i][2]-1][cmap[i][3]-1] = cmap[i][0]; if (cmap[i][0] == me) { myloc[0] = cmap[i][1] - 1; myloc[1] = cmap[i][2] - 1; myloc[2] = cmap[i][3] - 1; } } // proc IDs of neighbors int minus,plus; grid_shift(myloc[0],procgrid[0],minus,plus); procneigh[0][0] = grid2proc[minus][myloc[1]][myloc[2]]; procneigh[0][1] = grid2proc[plus][myloc[1]][myloc[2]]; grid_shift(myloc[1],procgrid[1],minus,plus); procneigh[1][0] = grid2proc[myloc[0]][minus][myloc[2]]; procneigh[1][1] = grid2proc[myloc[0]][plus][myloc[2]]; grid_shift(myloc[2],procgrid[2],minus,plus); procneigh[2][0] = grid2proc[myloc[0]][myloc[1]][minus]; procneigh[2][1] = grid2proc[myloc[0]][myloc[1]][plus]; memory->destroy(cmap); } /* ---------------------------------------------------------------------- output mapping of processors to 3d grid to file ------------------------------------------------------------------------- */ void ProcMap::output(char *file, int *procgrid, int ***grid2proc) { int me,nprocs; MPI_Comm_rank(world,&me); MPI_Comm_size(world,&nprocs); FILE *fp; if (me == 0) { fp = fopen(file,"w"); if (fp == NULL) error->one(FLERR,"Cannot open processors output file"); fprintf(fp,"LAMMPS mapping of processors to 3d grid\n"); fprintf(fp,"partition = %d\n",universe->iworld+1); fprintf(fp,"Px Py Pz = %d %d %d\n",procgrid[0],procgrid[1],procgrid[2]); fprintf(fp,"world-ID universe-ID original-ID: I J K: name\n\n"); } // find me in the grid int ime,jme,kme; for (int i = 0; i < procgrid[0]; i++) for (int j = 0; j < procgrid[1]; j++) for (int k = 0; k < procgrid[2]; k++) if (grid2proc[i][j][k] == me) { ime = i; jme = j; kme = k; } // polled comm of grid mapping info from each proc to proc 0 int tmp; int vec[6]; char procname[MPI_MAX_PROCESSOR_NAME+1]; vec[0] = me; vec[1] = universe->me; MPI_Comm_rank(universe->uorig,&vec[2]); vec[3] = ime + 1; vec[4] = jme + 1; vec[5] = kme + 1; int len; MPI_Get_processor_name(procname,&len); procname[len] = '\0'; if (me == 0) { for (int iproc = 0; iproc < nprocs; iproc++) { if (iproc) { MPI_Send(&tmp,0,MPI_INT,iproc,0,world); MPI_Recv(vec,6,MPI_INT,iproc,0,world,MPI_STATUS_IGNORE); MPI_Recv(procname,MPI_MAX_PROCESSOR_NAME+1,MPI_CHAR, iproc,0,world,MPI_STATUS_IGNORE); } fprintf(fp,"%d %d %d: %d %d %d: %s\n", vec[0],vec[1],vec[2],vec[3],vec[4],vec[5],procname); } } else { MPI_Recv(&tmp,0,MPI_INT,0,0,world,MPI_STATUS_IGNORE); MPI_Send(vec,6,MPI_INT,0,0,world); MPI_Send(procname,strlen(procname)+1,MPI_CHAR,0,0,world); } // close output file if (me == 0) fclose(fp); } /* ---------------------------------------------------------------------- generate all possible 3-integer factorizations of N store them in factors if non-NULL return # of factorizations ------------------------------------------------------------------------- */ int ProcMap::factor(int n, int **factors) { int i,j,nyz; int m = 0; for (i = 1; i <= n; i++) { if (n % i) continue; nyz = n/i; for (j = 1; j <= nyz; j++) { if (nyz % j) continue; if (factors) { factors[m][0] = i; factors[m][1] = j; factors[m][2] = nyz/j; } m++; } } return m; } /* ---------------------------------------------------------------------- create N1*N2 new factors (procs) from factors1 (nodes) and factors2 (cores) store index of corresponding core factors in factors[][3] ------------------------------------------------------------------------- */ int ProcMap::combine_factors(int n1, int **factors1, int n2, int **factors2, int **factors) { int m = 0; for (int i = 0; i < n1; i++) for (int j = 0; j < n2; j++) { factors[m][0] = factors1[i][0]*factors2[j][0]; factors[m][1] = factors1[i][1]*factors2[j][1]; factors[m][2] = factors1[i][2]*factors2[j][2]; factors[m][3] = j; m++; } return n1*n2; } /* ---------------------------------------------------------------------- remove any factors where Pz != 1 for 2d ------------------------------------------------------------------------- */ int ProcMap::cull_2d(int n, int **factors, int m) { int i = 0; while (i < n) { if (factors[i][2] != 1) { for (int j = 0; j < m; j++) factors[i][j] = factors[n-1][j]; n--; } else i++; } return n; } /* ---------------------------------------------------------------------- remove any factors that do not match non-zero user_factors Px,Py,Pz ------------------------------------------------------------------------- */ int ProcMap::cull_user(int n, int **factors, int m, int *user_factors) { int i = 0; while (i < n) { int flag = 0; if (user_factors[0] && factors[i][0] != user_factors[0]) flag = 1; if (user_factors[1] && factors[i][1] != user_factors[1]) flag = 1; if (user_factors[2] && factors[i][2] != user_factors[2]) flag = 1; if (flag) { for (int j = 0; j < m; j++) factors[i][j] = factors[n-1][j]; n--; } else i++; } return n; } /* ---------------------------------------------------------------------- remove any factors that do not match settings from other partition MULTIPLE = other Nx,Ny,Nz must be multiple of my Px,Py,Pz where Nx,Ny,Nz = node grid = procgrid/coregrid ------------------------------------------------------------------------- */ int ProcMap::cull_other(int n, int **factors, int m, int other_style, int *other_procgrid, int *other_coregrid) { int i = 0; while (i < n) { if (other_style == MULTIPLE) { int flag = 0; if ((other_procgrid[0]/other_coregrid[0]) % factors[i][0]) flag = 1; if ((other_procgrid[1]/other_coregrid[1]) % factors[i][1]) flag = 1; if ((other_procgrid[2]/other_coregrid[2]) % factors[i][2]) flag = 1; if (flag) { for (int j = 0; j < m; j++) factors[i][j] = factors[n-1][j]; n--; } else i++; } } return n; } /* ---------------------------------------------------------------------- choose best factors from list of Npossible factors best = minimal surface area of sub-domain return best = 3 factors return index of best factors in factors ------------------------------------------------------------------------- */ int ProcMap::best_factors(int npossible, int **factors, int *best, const int sx, const int sy, const int sz) { // determine cross-sectional areas for orthogonal and triclinic boxes // for triclinic, area = cross product of 2 edge vectors stored in h matrix // area[3] = surface area 3 box faces divided by sx,sy,sz // area[0] = xy, area[1] = xz, area[2] = yz double area[3]; if (domain->triclinic == 0) { area[0] = domain->xprd * domain->yprd / (sx*sy); area[1] = domain->xprd * domain->zprd / (sx*sz); area[2] = domain->yprd * domain->zprd / (sy*sz); } else { double *h = domain->h; double a[3],b[3],c[3]; a[0] = h[0]; a[1] = 0.0; a[2] = 0.0; b[0] = h[5]; b[1] = h[1]; b[2] = 0.0; MathExtra::cross3(a,b,c); area[0] = sqrt(c[0]*c[0] + c[1]*c[1] + c[2]*c[2]) / (sx*sy); a[0] = h[0]; a[1] = 0.0; a[2] = 0.0; b[0] = h[4]; b[1] = h[3]; b[2] = h[2]; MathExtra::cross3(a,b,c); area[1] = sqrt(c[0]*c[0] + c[1]*c[1] + c[2]*c[2]) / (sx*sz); a[0] = h[5]; a[1] = h[1]; a[2] = 0.0; b[0] = h[4]; b[1] = h[3]; b[2] = h[2]; MathExtra::cross3(a,b,c); area[2] = sqrt(c[0]*c[0] + c[1]*c[1] + c[2]*c[2]) / (sy*sz); } int index; double surf; double bestsurf = 2.0 * (area[0]+area[1]+area[2]); for (int m = 0; m < npossible; m++) { surf = area[0]/factors[m][0]/factors[m][1] + area[1]/factors[m][0]/factors[m][2] + area[2]/factors[m][1]/factors[m][2]; if (surf < bestsurf) { bestsurf = surf; best[0] = factors[m][0]; best[1] = factors[m][1]; best[2] = factors[m][2]; index = m; } } return index; } /* ---------------------------------------------------------------------- minus,plus = indices of neighboring processors in a dimension ------------------------------------------------------------------------- */ void ProcMap::grid_shift(int myloc, int nprocs, int &minus, int &plus) { minus = myloc - 1; if (minus < 0) minus = nprocs - 1; plus = myloc + 1; if (plus == nprocs) plus = 0; } diff --git a/src/region_cylinder.h b/src/region_cylinder.h index bc5255b34..e065c7d3f 100644 --- a/src/region_cylinder.h +++ b/src/region_cylinder.h @@ -1,83 +1,83 @@ /* -*- c++ -*- ---------------------------------------------------------- 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 REGION_CLASS RegionStyle(cylinder,RegCylinder) #else #ifndef LMP_REGION_CYLINDER_H #define LMP_REGION_CYLINDER_H #include "region.h" namespace LAMMPS_NS { class RegCylinder : public Region { friend class FixPour; public: RegCylinder(class LAMMPS *, int, char **); ~RegCylinder(); void init(); int inside(double, double, double); int surface_interior(double *, double); int surface_exterior(double *, double); void shape_update(); void set_velocity_shape(); void velocity_contact_shape(double *, double *); private: char axis; double c1,c2; double radius; double lo,hi; int rstyle,rvar; char *rstr; void variable_check(); }; } #endif #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Cannot use region INF or EDGE when box does not exist Regions that extend to the box boundaries can only be used after the create_box command has been used. E: Variable evaluation in region gave bad value Variable returned a radius < 0.0. E: Variable name for region cylinder does not exist Self-explanatory. E: Variable for region cylinder is invalid style -Only equal-style varaibles are allowed. +Only equal-style variables are allowed. */ diff --git a/src/region_sphere.h b/src/region_sphere.h index 9e1fb4e73..c003a9165 100644 --- a/src/region_sphere.h +++ b/src/region_sphere.h @@ -1,74 +1,74 @@ /* -*- c++ -*- ---------------------------------------------------------- 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 REGION_CLASS RegionStyle(sphere,RegSphere) #else #ifndef LMP_REGION_SPHERE_H #define LMP_REGION_SPHERE_H #include "region.h" namespace LAMMPS_NS { class RegSphere : public Region { public: RegSphere(class LAMMPS *, int, char **); ~RegSphere(); void init(); int inside(double, double, double); int surface_interior(double *, double); int surface_exterior(double *, double); void shape_update(); void set_velocity_shape(); void velocity_contact_shape(double *, double *); private: double xc,yc,zc; double radius; int rstyle,rvar; char *rstr; void variable_check(); }; } #endif #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Variable evaluation in region gave bad value Variable returned a radius < 0.0. E: Variable name for region sphere does not exist Self-explanatory. E: Variable for region sphere is invalid style -Only equal-style varaibles are allowed. +Only equal-style variables are allowed. */ diff --git a/src/thermo.h b/src/thermo.h index ddbe3cc95..d6b3822ef 100644 --- a/src/thermo.h +++ b/src/thermo.h @@ -1,399 +1,399 @@ /* -*- c++ -*- ---------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifndef LMP_THERMO_H #define LMP_THERMO_H #include "pointers.h" namespace LAMMPS_NS { class Thermo : protected Pointers { friend class MinCG; // accesses compute_pe public: char *style; int normflag; // 0 if do not normalize by atoms, 1 if normalize int modified; // 1 if thermo_modify has been used, else 0 int lostflag; // IGNORE,WARN,ERROR int lostbond; // ditto for atoms in bonds Thermo(class LAMMPS *, int, char **); ~Thermo(); void init(); bigint lost_check(); void modify_params(int, char **); void header(); void compute(int); int evaluate_keyword(char *, double *); private: char *line; char **keyword; int *vtype; int nfield,nfield_initial; int me; char **format; char *format_line_user; char *format_float_user,*format_int_user,*format_bigint_user; char **format_column_user; char *format_float_one_def,*format_float_multi_def; char *format_int_one_def,*format_int_multi_def; char format_multi[128]; char format_bigint_one_def[8],format_bigint_multi_def[8]; int normvalue; // use this for normflag unless natoms = 0 int normuserflag; // 0 if user has not set, 1 if has int normuser; int firststep; int lostbefore; int flushflag,lineflag; double last_tpcpu,last_spcpu; double last_time; bigint last_step; bigint natoms; // data used by routines that compute single values int ivalue; // integer value to print double dvalue; // double value to print bigint bivalue; // big integer value to print int ifield; // which field in thermo output is being computed int *field2index; // which compute,fix,variable calcs this field int *argindex1; // indices into compute,fix scalar,vector int *argindex2; // data for keyword-specific Compute objects // index = where they are in computes list // id = ID of Compute objects // Compute * = ptrs to the Compute objects int index_temp,index_press_scalar,index_press_vector,index_pe; char *id_temp,*id_press,*id_pe; class Compute *temperature,*pressure,*pe; int ncompute; // # of Compute objects called by thermo char **id_compute; // their IDs int *compute_which; // 0/1/2 if should call scalar,vector,array class Compute **computes; // list of ptrs to the Compute objects int nfix; // # of Fix objects called by thermo char **id_fix; // their IDs class Fix **fixes; // list of ptrs to the Fix objects int nvariable; // # of variables evaulated by thermo char **id_variable; // list of variable names int *variables; // list of Variable indices // private methods void allocate(); void deallocate(); void parse_fields(char *); int add_compute(const char *, int); int add_fix(const char *); int add_variable(const char *); typedef void (Thermo::*FnPtr)(); void addfield(const char *, FnPtr, int); FnPtr *vfunc; // list of ptrs to functions void compute_compute(); // functions that compute a single value void compute_fix(); // via calls to Compute,Fix,Variable classes void compute_variable(); // functions that compute a single value // customize a new keyword by adding a method prototype void compute_step(); void compute_elapsed(); void compute_elapsed_long(); void compute_dt(); void compute_time(); void compute_cpu(); void compute_tpcpu(); void compute_spcpu(); void compute_cpuremain(); void compute_part(); void compute_timeremain(); void compute_atoms(); void compute_temp(); void compute_press(); void compute_pe(); void compute_ke(); void compute_etotal(); void compute_enthalpy(); void compute_evdwl(); void compute_ecoul(); void compute_epair(); void compute_ebond(); void compute_eangle(); void compute_edihed(); void compute_eimp(); void compute_emol(); void compute_elong(); void compute_etail(); void compute_vol(); void compute_density(); void compute_lx(); void compute_ly(); void compute_lz(); void compute_xlo(); void compute_xhi(); void compute_ylo(); void compute_yhi(); void compute_zlo(); void compute_zhi(); void compute_xy(); void compute_xz(); void compute_yz(); void compute_xlat(); void compute_ylat(); void compute_zlat(); void compute_bonds(); void compute_angles(); void compute_dihedrals(); void compute_impropers(); void compute_pxx(); void compute_pyy(); void compute_pzz(); void compute_pxy(); void compute_pyz(); void compute_pxz(); void compute_fmax(); void compute_fnorm(); void compute_nbuild(); void compute_ndanger(); void compute_cella(); void compute_cellb(); void compute_cellc(); void compute_cellalpha(); void compute_cellbeta(); void compute_cellgamma(); }; } #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Could not find thermo compute ID Compute ID specified in thermo_style command does not exist. E: Could not find thermo fix ID Fix ID specified in thermo_style command does not exist. E: Thermo and fix not computed at compatible times Fixes generate values on specific timesteps. The thermo output does not match these timesteps. E: Could not find thermo variable name Self-explanatory. E: Too many total atoms See the setting for bigint in the src/lmptype.h file. E: Lost atoms: original %ld current %ld Lost atoms are checked for each time thermo output is done. See the thermo_modify lost command for options. Lost atoms usually indicate bad dynamics, e.g. atoms have been blown far out of the simulation -box, or moved futher than one processor's sub-domain away before +box, or moved further than one processor's sub-domain away before reneighboring. W: Lost atoms: original %ld current %ld Lost atoms are checked for each time thermo output is done. See the thermo_modify lost command for options. Lost atoms usually indicate bad dynamics, e.g. atoms have been blown far out of the simulation -box, or moved futher than one processor's sub-domain away before +box, or moved further than one processor's sub-domain away before reneighboring. E: Thermo style does not use temp Cannot use thermo_modify to set this parameter since the thermo_style is not computing this quantity. E: Could not find thermo_modify temperature ID The compute ID needed by thermo style custom to compute temperature does not exist. E: Thermo_modify temperature ID does not compute temperature The specified compute ID does not compute temperature. W: Temperature for thermo pressure is not for group all User-assigned temperature to thermo via the thermo_modify command does not compute temperature for all atoms. Since thermo computes a global pressure, the kinetic energy contribution from the temperature is assumed to also be for all atoms. Thus the pressure printed by thermo could be inaccurate. E: Pressure ID for thermo does not exist The compute ID needed to compute pressure for thermodynamics does not exist. E: Thermo style does not use press Cannot use thermo_modify to set this parameter since the thermo_style is not computing this quantity. E: Could not find thermo_modify pressure ID The compute ID needed by thermo style custom to compute pressure does not exist. E: Thermo_modify pressure ID does not compute pressure The specified compute ID does not compute pressure. E: Thermo_modify int format does not contain d character Self-explanatory. E: Could not find thermo custom compute ID The compute ID needed by thermo style custom to compute a requested quantity does not exist. E: Thermo compute does not compute scalar Self-explanatory. E: Thermo compute does not compute vector Self-explanatory. E: Thermo compute vector is accessed out-of-range Self-explanatory. E: Thermo compute does not compute array Self-explanatory. E: Thermo compute array is accessed out-of-range Self-explanatory. E: Could not find thermo custom fix ID The fix ID needed by thermo style custom to compute a requested quantity does not exist. E: Thermo fix does not compute scalar Self-explanatory. E: Thermo fix does not compute vector Self-explanatory. E: Thermo fix vector is accessed out-of-range Self-explanatory. E: Thermo fix does not compute array Self-explanatory. E: Thermo fix array is accessed out-of-range Self-explanatory. E: Could not find thermo custom variable name Self-explanatory. E: Thermo custom variable is not equal-style variable Only equal-style variables can be output with thermodynamics, not atom-style variables. E: Thermo custom variable cannot be indexed Self-explanatory. E: Unknown keyword in thermo_style custom command One or more specified keywords are not recognized. E: This variable thermo keyword cannot be used between runs Keywords that refer to time (such as cpu, elapsed) do not make sense in between runs. E: Thermo keyword in variable requires thermo to use/init temp You are using a thermo keyword in a variable that requires temperature to be calculated, but your thermo output does not use it. Add it to your thermo output. E: Compute used in variable thermo keyword between runs is not current Some thermo keywords rely on a compute to calculate their value(s). Computes cannot be invoked by a variable in between runs. Thus they must have been evaluated on the last timestep of the previous run in order for their value(s) to be accessed. See the doc page for the variable command for more info. E: Thermo keyword in variable requires thermo to use/init press You are using a thermo keyword in a variable that requires pressure to be calculated, but your thermo output does not use it. Add it to your thermo output. E: Thermo keyword in variable requires thermo to use/init pe You are using a thermo keyword in a variable that requires potential energy to be calculated, but your thermo output does not use it. Add it to your thermo output. E: Energy was not tallied on needed timestep You are using a thermo keyword that requires potentials to have tallied energy, but they didn't on this timestep. See the variable doc page for ideas on how to make this work. */ diff --git a/tools/eff/cfg2lammps.py b/tools/eff/cfg2lammps.py index 1bf364aa9..49830dbb0 100644 --- a/tools/eff/cfg2lammps.py +++ b/tools/eff/cfg2lammps.py @@ -1,377 +1,377 @@ #!/usr/local/bin/python-2.5/bin/python Info=""" Module name: cfg2lammps.py Author: (c) Andres Jaramillo-Botero California Institute of Technology ajaramil@wag.caltech.edu Project: pEFF Version: August 2009 Reads in an eff .cfg file and produces the corresponding lammps data and input files NOTE: Unsupported functions will be reported in the output log 12/2010: Added support for fixed-core and pseudo-core structures """ # import essentials: import sys, os from math import log10 from shutil import rmtree from getopt import gnu_getopt as getopt import numpy def printHelp(): print Info print "Usage: python cfg2lammps cfgfile\n" return general=""" # Created %s # General parameters variable sname index %s log ${sname}.log units electron newton on boundary %s atom_style electron read_data data.${sname} pair_style eff/cut %s pair_coeff * * compute energies all pair eff/cut variable eke equal c_energies[1] variable epauli equal c_energies[2] variable estatics equal c_energies[3] variable errestrain equal c_energies[4] communicate single vel yes compute peratom all stress/atom compute p all reduce sum c_peratom[1] c_peratom[2] c_peratom[3] variable press equal -(c_p[1]+c_p[2]+c_p[3])/(3*vol) compute effTemp all temp/eff compute effPress all pressure effTemp thermo %s thermo_style custom step etotal pe ke v_eke v_epauli v_estatics v_errestrain temp press v_press thermo_modify temp effTemp press effPress """ #%(date,name,boundary,cutoff,period) minimize=""" # Minimization min_style cg dump 1 %s xyz %s ${sname}.min.xyz dump 2 %s custom %s ${sname}.min.lammpstrj id type q spin eradius x y z fx fy fz erforce min_modify line quadratic minimize 0 1.0e-5 %s %s undump 1 undump 2 """ #%(group,period,group,period,iterations,fcalls) single_pt=""" # Single point energy run 0 """ dynamics=""" # %s Dynamics timestep %s fix %s dump 1 %s custom %s ${sname}.%s.lammpstrj id type q spin eradius x y z dump 2 %s custom %s ${sname}.%s.xyz run %s unfix 1 undump 1 undump 2 """ task={'single_pt':single_pt,'minimize':minimize,'dynamics':dynamics} q2m={1:'1.007940',2:'4.002602',3:'6.941000',4:'9.012182',5:'10.811000',6:'12.010700',7:'14.006700',8:'15.999400', 9:'18.9984032',10:'20.179700',11:'22.98976928',12:'24.305000',13:'26.9815386',14:'28.085500',15:'30.973762', 16:'32.065000',17:'35.453000',18:'39.948000'} def generate_lammps_input(infile): # Defaults values ensemble={"nve":"1 %s nve/eff",'nvt':"1 %s nvt/eff %s %s %s %s",'npt':"1 %s npt/eff %s %s %s %s %s %s"} boundary="f f f" xbound="-1000.000 1000.0 xlo xhi\n" ybound="-1000.000 1000.0 ylo yhi\n" zbound="-1000.000 1000.0 zlo zhi\n" cutoff=1000.0 period="1" emass=0 vels="" datafile=open("data."+infile[:-4],'w') scriptfile=open("in."+infile[:-4],'w') print "Reading %s ... [WAIT]"%infile, fin = open(infile,'r') lines = fin.xreadlines() print 7*"\b"+"[DONE]" numcores=0 numnuclei=0 numelec=0 cores={} nuclei={} electrons={} masses=[] massstr="Masses\n\n" types=1 q2type={} Tflag=False # Default ensemble is NVE steps='1000' print "Extracting run parameters from %s ... "%(infile), for line in lines: # 1st level keywords if line.find("@params")==0: flag='params' continue elif line.find("@cores")==0: flag='cores' continue elif line.find("@nuclei")==0: flag='nuclei' continue elif line.find("@electrons")==0: flag='electrons' continue elif line.find("@nuc_velocities")==0: flag='n_vels' continue elif line.find("@elec_velocities")==0: flag='e_vels' continue elif line.find("@nuc_masses")==0: flag='n_mass' continue elif line.find("@elec_masses")==0: flag='e_mass' continue elif line.find("@restraints")==0: flag='restraints' continue # 2nd level keywords if flag=='params': if line.find("calc")>=0: op=line.split()[2] if line.find("print_every")>=0: period=line.split()[2] if line.find("num_steps")>=0: steps=line.split()[2] if line.find("min_freeze")>=0: setforce="velocity\t% set 0.0 0.0 0.0\nfix\tfreeze %s setforce 0.0 0.0 0.0"%(line.split()[2],line.split()[2]) if line.find("thermostat")>=0: tflag=True #ensemble="fix\t1 all nvt/eff " if line.find("start_temperature")>=0: Tstart=line.split()[2] #ensemble+=Tstart if line.find("end_temperature")>=0: Tstop=line.split()[2] #ensemble+=Tstop if line.find("andersen_coupling")>=0 or line.find("nose_hoover_coupling")>=0: Tdamp=line.split()[2] #ensemble+=Tdamp if line.find("dt")>=0: dt=line.split()[2] if line.find("electron_mass")>=0: emass=line.split()[2] if line.find("adaptive_step_size")>=0: continue if line.find("adaptive_energy")>=0: continue if line.find("e_field_freq")>=0: continue if line.find("e_field_packet_duration")>=0: continue if line.find("e_field")>=0: field=line.split()[2:5] efield="fix\field all efield %s %s %s"%(field[0],field[1],field[2]) if line.find("e_field_packet_duration")>=0: continue if line.find("set_limit")>=0: - continue # need to add this contraint + continue # need to add this constraint if line.find("set_limit_stiffness")>=0: continue if line.find("output_position")>=0: dump_pos="dump\t1 all custom %s ${sname}.lammpstrj id type q spin eradius x y z "%(period) if line.find("output_velocities")>=0: dump_pos+="vx vy vz " if line.find("output_energy_forces")>=0: dump_pos="compute\tenergy all pe/atom\n"+dump_pos dump_pos+="c_energy fx fy fz\n" if line.find("output_restart")>=0: restart="restart\t%s ${sname}.restart1 ${sname}.restart2"%(period) if line.find("output_restraints")>=0: continue if line.find("ewald_re_cutoff")>=0 or line.find("ewald_autoset")>=0 or line.find("ewald_log_precision")>=0 or line.find("ewald_max_re")>=0 or \ line.find("ewald_r_cutoff")>=0 or line.find("ewald_k_cutoff")>=0 or line.find("ewald_nuc_r")>=0: continue if line.find("periodic")>=0: bounds=line.split()[2] if bounds=="True": boundary="p p p" elif bounds=="minimage_x": boundary="p f f" elif bounds=="minimage_xy": boundary="p p f" elif bounds=="minimage_y": boundary="f p f" elif bounds=="minimage_xyz": boundary="p p p" elif bounds=="minimage_z": boundary="f f p" if line.find("x_bound")>=0: xbnds=line.split()[2:4] xbound="%s %s xlo xhi\n"%(xbnds[0],xbnds[1]) if line.find("y_bound")>=0: ybnds=line.split()[2:4] ybound="%s %s ylo yhi\n"%(ybnds[0],ybnds[1]) if line.find("z_bound")>=0: zbnds=line.split()[2:4] zbound="%s %s zlo zhi\n"%(zbnds[0],zbnds[1]) if line.find("taper_cutoff")>=0: cutoff=line.split()[2] continue if flag=='cores' and len(line)>1: numcores+=1 ln=line.split() nc=' '.join(ln[0:3]) q=ln[3] spin='3' radius=ln[4] m=q2m[int(float(q))] if m not in masses: masses.append(m) massstr+="%d %s\n"%(types,m) q2type[q]=types types+=1 cores[numcores]=[nc,q,spin,radius] continue if flag=='nuclei' and len(line)>1: numnuclei+=1 ln=line.split() np=' '.join(ln[0:3]) q=ln[3] m=q2m[int(float(q))] if m not in masses: masses.append(m) massstr+="%d %s\n"%(types,m) q2type[q]=types types+=1 nuclei[numnuclei]=[np,q] continue if flag=='electrons' and len(line)>1: numelec+=1 ln=line.split() ep=' '.join(ln[0:3]) spin=ln[3] radius=ln[4] electrons[numelec]=[ep,spin,radius] if numelec==1: if emass!=0: massstr+="%d %s\n\n"%(types,emass) # electron mass=1 else: massstr+="%d 1.000000\n\n"%(types) continue if flag=='n_vels' and len(line)>1: vels+=line+" 0.0" continue if flag=='e_vels' and len(line)>1: ln=line.split() ln[0]=ln[0]+numnuclei vels+=ln[0]+" "+ln[1]+" "+ln[2]+" "+ln[3]+" "+ln[4]+"\n" continue if flag=='n_mass' and len(line)>1: print "Setting nuclear masses is unsupported\n" continue if flag=='e_mass' and len(line)>1: print "Setting electron masses is unsupported\n" continue print "\bDone" # Build data file print "Writing datafile to %s ... "%('data.'+infile), sys.stdout.flush() print "\b"*19+"General section ", datafile.writelines("Created using cfg2lammps (c) AJB-2009\n\n%d atoms\n%d atom types\n\n%s%s%s\n"%(numcores+numnuclei+numelec,types,xbound,ybound,zbound)) print "\b"*19+"Masses section ", datafile.writelines(massstr) print "\b"*19+"Atoms section ", datafile.writelines("Atoms\n\n") for n in range(numcores): datafile.writelines("%d %d %2.2f %s %s %s\n"%(n+1,q2type[cores[n+1][1]],float(cores[n+1][1]),cores[n+1][2],cores[n+1][3],cores[n+1][0])) for n in range(numnuclei): datafile.writelines("%d %d %2.2f 0 0.0 %s\n"%(n+numcores+1,q2type[nuclei[n+1][1]],float(nuclei[n+1][1]),nuclei[n+1][0])) for e in range(numelec): datafile.write("%d %d 0.0 %s %s %s\n"%(e+numnuclei+numcores+1,types,electrons[e+1][1],electrons[e+1][2],electrons[e+1][0])) print "\b"*19+"Velocities section\n", datafile.writelines(vels) datafile.writelines("\n") print "DONE .... GOODBYE !!" datafile.close() # Build input script import datetime scriptfile.writelines(general%(datetime.date.today(),infile[:-4],boundary,cutoff,period)) if op=='minimize': scriptfile.writelines(minimize%('all',period,'all',period,steps,'10000')) #%(group,period,group,period,iterations,fcalls) elif op=='single_pt': scriptfile.writelines(single_pt%()) elif op=='dynamics': if Tflag==True: scriptfile.writelines(dynamics%('NVT',dt,ensemble['nvt']%('all',Tstart,Tstop,Tdamp,''),'all',period,'nvt','all',period,'nve',steps)) #%(ensemble,dt,group,ensemble%(group,tstart,tstop,tdamp,options)) else: scriptfile.writelines(dynamics%('NVE',dt,ensemble['nve']%('all'),'all',period,'nve','all',period,'nve',steps)) #%(ensemble,dt,group,ensemble%(group)) scriptfile.writelines("\n") if __name__ == '__main__': # set defaults # check for input: opts, argv = getopt(sys.argv[1:], 'h') # if no input, print help and exit if len(argv) != 1: printHelp() sys.exit(1) else: infile=argv[0] # read options for opt, arg in opts: if opt == '-h': # -h: print help printHelp() generate_lammps_input(infile) diff --git a/tools/i-pi/ipi/engine/initializer.py b/tools/i-pi/ipi/engine/initializer.py index 466677938..de629c7f9 100644 --- a/tools/i-pi/ipi/engine/initializer.py +++ b/tools/i-pi/ipi/engine/initializer.py @@ -1,549 +1,549 @@ """Contains the classes that are used to initialize data in the simulation. Copyright (C) 2013, Joshua More and Michele Ceriotti This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . These classes can either be used to restart a simulation with some different data or used to start a calculation. Any data given in these classes will overwrite data given elsewhere. Classes: Initializer: Holds the functions that are required to initialize objects in the code. Data can be initialized from a file, or according to a particular parameter. An example of the former would be initializing the configurations from a xyz file, an example of the latter would be initializing the velocities according to the physical temperature. InitBase: Simple class that reads data from a string or file. InitIndexed: The same as init base, but can also optionally hold information about which atom or bead to initialize from. Functions: init_xyz: Reads beads data from a xyz file. init_pdb: Reads beads and cell data from a pdb file. init_chk: Reads beads, cell and thermostat data from a checkpoint file. init_beads: Initializes a beads object from an Initializer object. init_vector: Initializes a vector from an Initializer object. set_vector: Initializes a vector from another vector. """ import numpy as np from ipi.engine.beads import Beads from ipi.engine.cell import Cell from ipi.engine.normalmodes import NormalModes from ipi.engine.ensembles import Ensemble from ipi.utils.io.io_xyz import read_xyz from ipi.utils.io.io_pdb import read_pdb from ipi.utils.io.io_xml import xml_parse_file from ipi.utils.depend import dobject from ipi.utils.units import Constants, unit_to_internal from ipi.utils.nmtransform import nm_rescale from ipi.utils.messages import verbosity, warning, info __all__ = ['Initializer', 'InitBase', 'InitIndexed'] class InitBase(dobject): """Base class for initializer objects. Attributes: value: A duck-typed stored value. mode: A string that determines how the value is to be interpreted. units: A string giving which unit the value is in. """ def __init__(self, value="", mode="", units="", **others): """Initializes InitFile. Args: value: A string which specifies what value to initialize the simulation property to. - mode: A string specifiying what style of initialization should be + mode: A string specifying what style of initialization should be used to read the data. units: A string giving which unit the value is in. """ self.value = value self.mode = mode self.units = units for (o, v) in others.items(): self.__dict__[o] = v class InitIndexed(InitBase): """Class to initialize objects which can be set for a particular bead. Attributes: index: Which atom to initialize the value of. bead: Which bead to initialize the value of. """ def __init__(self, value="", mode="", units="", index=-1, bead=-1): """Initializes InitFile. Args: value: A string which specifies what value to initialize the simulation property to. - mode: A string specifiying what style of initialization should be + mode: A string specifying what style of initialization should be used to read the data. units: A string giving which unit the value is in. index: Which atom to initialize the value of. bead: Which bead to initialize the value of. """ super(InitIndexed,self).__init__(value=value, mode=mode, units=units, index=index, bead=bead) def init_xyz(filename): """Reads an xyz file and returns the data contained in it. Args: filename: A string giving the name of the xyz file to be read from. Returns: A list of Atoms objects as read from each frame of the xyz file. """ rfile = open(filename,"r") ratoms = [] while True: #while loop, so that more than one configuration can be given #so multiple beads can be initialized at once. try: myatoms = read_xyz(rfile) except EOFError: break ratoms.append(myatoms) return ratoms def init_pdb(filename): """Reads an pdb file and returns the data contained in it. Args: filename: A string giving the name of the pdb file to be read from. Returns: A list of Atoms objects as read from each frame of the pdb file, and a Cell object as read from the final pdb frame. """ rfile = open(filename,"r") ratoms = [] while True: #while loop, so that more than one configuration can be given #so multiple beads can be initialized at once. try: myatoms, rcell = read_pdb(rfile) except EOFError: break ratoms.append(myatoms) return ( ratoms, rcell ) # if multiple frames, the last cell is returned def init_chk(filename): """Reads an checkpoint file and returns the data contained in it. Args: filename: A string giving the name of the checkpoint file to be read from. Returns: A Beads object, Cell object and Thermostat object as read from the checkpoint file. """ # reads configuration from a checkpoint file rfile = open(filename,"r") xmlchk = xml_parse_file(rfile) # Parses the file. from ipi.inputs.simulation import InputSimulation simchk = InputSimulation() simchk.parse(xmlchk.fields[0][1]) rcell = simchk.cell.fetch() rbeads = simchk.beads.fetch() rthermo = simchk.ensemble.thermostat.fetch() return (rbeads, rcell, rthermo) def init_beads(iif, nbeads): """A file to initialize a beads object from an appropriate initializer object. Args: iif: An Initializer object which has information on the bead positions. nbeads: The number of beads. Raises: ValueError: If called using an Initializer object with a 'manual' mode. """ mode = iif.mode; value = iif.value if mode == "xyz" or mode == "pdb": if mode == "xyz": ratoms = init_xyz(value) if mode == "pdb": ratoms = init_pdb(value)[0] rbeads = Beads(ratoms[0].natoms,len(ratoms)) for i in range(len(ratoms)): rbeads[i] = ratoms[i] elif mode == "chk": rbeads = init_chk(value)[0] elif mode == "manual": raise ValueError("Cannot initialize manually a whole beads object.") return rbeads def init_vector(iif, nbeads, momenta=False): """A file to initialize a vector from an appropriate initializer object. Args: iif: An Initializer object specifying the value of a vector. nbeads: The number of beads. momenta: If bead momenta rather than positions are being initialized from a checkpoint file, this is set to True. """ mode = iif.mode; value = iif.value if mode == "xyz" or mode == "pdb": rq = init_beads(iif, nbeads).q elif mode == "chk": if momenta: rq = init_beads(iif, nbeads).p else: rq = init_beads(iif, nbeads).q elif mode == "manual": rq = value # determines the size of the input data if mode == "manual": if iif.bead >= 0: # if there is a bead specifier then we return a single bead slice nbeads = 1 natoms = len(rq)/nbeads/3 rq.shape = (nbeads,3*natoms) return rq def set_vector(iif, dq, rq): """A file to initialize a vector from an another vector. If the first dimension is different, i.e. the two vectors correspond to a different number of beads, then the ring polymer contraction/expansion is used to rescale the original vector to the one used in the simulation, as described in the paper T. E. Markland and D. E. Manolopoulos, J. Chem. Phys. 129, 024105, (2008). Args: iif: An Initializer object specifying the value of a vector. dq: The vector to be initialized. rq: The vector to initialize from. """ (nbeads, natoms) = rq.shape; natoms /= 3 (dbeads, datoms) = dq.shape; datoms /= 3 # Check that indices make sense if iif.index < 0 and natoms != datoms: raise ValueError("Initialization tries to mix up structures with different atom numbers.") if iif.index >= datoms: raise ValueError("Cannot initialize single atom as atom index %d is larger than the number of atoms" % iif.index) if iif.bead >= dbeads: raise ValueError("Cannot initialize single bead as bead index %d is larger than the number of beads" % iif.bead) if iif.bead < 0: # we are initializing the path res = nm_rescale(nbeads,dbeads) # path rescaler if nbeads != dbeads: info(" # Initialize is rescaling from %5d beads to %5d beads" % (nbeads, dbeads), verbosity.low) if iif.index < 0: dq[:] = res.b1tob2(rq) else: # we are initializing a specific atom dq[:,3*iif.index:3*(iif.index+1)] = res.b1tob2(rq) else: # we are initializing a specific bead if iif.index < 0: dq[iif.bead] = rq else: dq[iif.bead,3*iif.index:3*(iif.index+1)] = rq class Initializer(dobject): """Class that deals with the initialization of data. This can either be used to initialize the atom positions and the cell data from a file, or to initialize them from a beads, atoms or cell object. Currently, we use a ring polymer contraction scheme to create a new beads object from one given in initialize if they have different numbers of beads, as described in the paper T. E. Markland and D. E. Manolopoulos, J. Chem. Phys. 129, 024105, (2008). If the new beads object has more beads than the beads object it was initialized from, we set the higher ring polymer normal modes to zero. Attributes: queue: A list of things to initialize. Each member of the list is a tuple of the form ('type', 'object'), where 'type' specifies what kind of initialization is being done, and 'object' gives the data to initialize it from. """ def __init__(self, nbeads=0, queue=None): """Initializes Initializer. Arguments: nbeads: The number of beads that we need in the simulation. Not necessarily the same as the number of beads of the objects we are initializing the data from. queue: A list of things to initialize. Each member of the list is a tuple of the form ('type', 'object'), where 'type' specifies what kind of initialization is being done, and 'object' gives the data to initialize it from. """ self.nbeads = nbeads if queue is None: self.queue = [] else: self.queue = queue def init_stage1(self, simul): """Initializes the simulation -- first stage. Takes a simulation object, and uses all the data in the initialization queue to fill up the beads and cell data needed to run the simulation. Args: simul: A simulation object to be initialized. Raises: ValueError: Raised if there is a problem with the initialization, if something that should have been has not been, or if the objects that have been specified are not compatible with each other. """ if simul.beads.nbeads == 0: fpos = fmom = fmass = flab = fcell = False # we don't have an explicitly defined beads object yet else: fpos = fmom = fmass = flab = fcell = True for (k,v) in self.queue: info(" # Initializer (stage 1) parsing " + str(k) + " object.", verbosity.high) if k == "cell": if fcell : warning("Overwriting previous cell parameters", verbosity.medium) if v.mode == "pdb": rh = init_pdb(v.value)[1].h elif v.mode == "chk": rh = init_chk(v.value)[1].h else: rh = v.value.reshape((3,3)) rh *= unit_to_internal("length",v.units,1.0) simul.cell.h = rh if simul.cell.V == 0.0: ValueError("Cell provided has zero volume") fcell = True elif k == "masses": if simul.beads.nbeads == 0: raise ValueError("Cannot initialize the masses before the size of the system is known") if fmass: warning("Overwriting previous atomic masses", verbosity.medium) if v.mode == "manual": rm = v.value else: rm = init_beads(v, self.nbeads).m rm *= unit_to_internal("mass",v.units,1.0) if v.bead < 0: # we are initializing the path if (fmom and fmass): warning("Rescaling momenta to make up for changed mass", verbosity.medium) simul.beads.p /= simul.beads.sm3 # go to mass-scaled momenta, that are mass-invariant if v.index < 0: simul.beads.m = rm else: # we are initializing a specific atom simul.beads.m[v.index:v.index+1] = rm if (fmom and fmass): # finishes correcting the momenta simul.beads.p *= simul.beads.sm3 # back to normal momenta else: raise ValueError("Cannot change the mass of a single bead") fmass = True elif k == "labels": if simul.beads.nbeads == 0: raise ValueError("Cannot initialize the labels before the size of the system is known") if flab: warning("Overwriting previous atomic labels", verbosity.medium) if v.mode == "manual": rn = v.value else: rn = init_beads(v, self.nbeads).names if v.bead < 0: # we are initializing the path if v.index < 0: simul.beads.names = rn else: # we are initializing a specific atom simul.beads.names[v.index:v.index+1] = rn else: raise ValueError("Cannot change the label of a single bead") flab = True elif k == "positions": if fpos: warning("Overwriting previous atomic positions", verbosity.medium) # read the atomic positions as a vector rq = init_vector(v, self.nbeads) rq *= unit_to_internal("length",v.units,1.0) (nbeads, natoms) = rq.shape; natoms /= 3 # check if we must initialize the simulation beads if simul.beads.nbeads == 0: if v.index >= 0: raise ValueError("Cannot initialize single atoms before the size of the system is known") simul.beads.resize(natoms,self.nbeads) set_vector(v, simul.beads.q, rq) fpos = True elif (k == "velocities" or k == "momenta") and v.mode == "thermal" : # intercept here thermal initialization, so we don't need to check further down if fmom: warning("Overwriting previous atomic momenta", verbosity.medium) if simul.beads.natoms == 0: raise ValueError("Cannot initialize momenta before the size of the system is known.") if not fmass: raise ValueError("Trying to resample velocities before having masses.") rtemp = v.value * unit_to_internal("temperature",v.units,1.0) if rtemp <= 0: warning("Using the simulation temperature to resample velocities", verbosity.low) rtemp = simul.ensemble.temp else: info(" # Resampling velocities at temperature %s %s" % (v.value, v.units), verbosity.low) # pull together a mock initialization to get NM masses right #without too much code duplication if v.bead >= 0: raise ValueError("Cannot thermalize a single bead") if v.index >= 0: rnatoms = 1 else: rnatoms = simul.beads.natoms rbeads = Beads(rnatoms, simul.beads.nbeads) if v.index < 0: rbeads.m[:] = simul.beads.m else: rbeads.m[:] = simul.beads.m[v.index] rnm = NormalModes(mode=simul.nm.mode, transform_method=simul.nm.transform_method, freqs=simul.nm.nm_freqs) rens = Ensemble(dt=simul.ensemble.dt, temp=simul.ensemble.temp) rnm.bind(rbeads,rens) # then we exploit the sync magic to do a complicated initialization # in the NM representation # with (possibly) shifted-frequencies NM rnm.pnm = simul.prng.gvec((rbeads.nbeads,3*rbeads.natoms))*np.sqrt(rnm.dynm3)*np.sqrt(rbeads.nbeads*rtemp*Constants.kb) if v.index < 0: simul.beads.p = rbeads.p else: simul.beads.p[:,3*v.index:3*(v.index+1)] = rbeads.p fmom = True elif k == "momenta": if fmom: warning("Overwriting previous atomic momenta", verbosity.medium) # read the atomic momenta as a vector rp = init_vector(v, self.nbeads, momenta = True) rp *= unit_to_internal("momentum",v.units,1.0) (nbeads, natoms) = rp.shape; natoms /= 3 # checks if we must initialize the simulation beads if simul.beads.nbeads == 0: if v.index >= 0 : raise ValueError("Cannot initialize single atoms before the size of the system is known") simul.beads.resize(natoms,self.nbeads) rp *= np.sqrt(self.nbeads/nbeads) set_vector(v, simul.beads.p, rp) fmom = True elif k == "velocities": if fmom: warning("Overwriting previous atomic momenta", verbosity.medium) # read the atomic velocities as a vector rv = init_vector(v, self.nbeads) rv *= unit_to_internal("velocity",v.units,1.0) (nbeads, natoms) = rv.shape; natoms /= 3 # checks if we must initialize the simulation beads if simul.beads.nbeads == 0 or not fmass: ValueError("Cannot initialize velocities before the masses of the atoms are known") simul.beads.resize(natoms,self.nbeads) warning("Initializing from velocities uses the previously defined masses -- not the masses inferred from the file -- to build momenta", verbosity.low) if v.index >= 0: rv *= simul.beads.m[v.index] elif v.bead >= 0: rv *= simul.beads.m3[0] else: rv *= simul.beads.m3 rv *= np.sqrt(self.nbeads/nbeads) set_vector(v, simul.beads.p, rv) fmom = True elif k == "thermostat": pass # thermostats must be initialised in a second stage if simul.beads.natoms == 0: raise ValueError("Initializer could not initialize the atomic positions") if simul.cell.V == 0: raise ValueError("Initializer could not initialize the cell") for i in range(simul.beads.natoms): if simul.beads.m[i] <= 0: raise ValueError("Initializer could not initialize the masses") if simul.beads.names[i] == "": raise ValueError("Initializer could not initialize the atom labels") if not fmom: warning("Momenta not specified in initialize. Will start with zero velocity if they are not specified in beads.", verbosity.low) def init_stage2(self, simul): """Initializes the simulation -- second stage. Takes a simulation object which has been fully generated, and restarts additional information such as the thermostat internal state. Args: simul: A simulation object to be initialized. Raises: ValueError: Raised if there is a problem with the initialization, if something that should have been has not been, or if the objects that have been specified are not compatible with each other. """ for (k,v) in self.queue: info(" # Initializer (stage 2) parsing " + str(k) + " object.", verbosity.high) if k == "gle": # read thermostat parameters from file if not ( hasattr(simul.ensemble, "thermostat") ): raise ValueError("Ensemble does not have a thermostat to initialize") if not ( hasattr(simul.ensemble.thermostat, "s") ): raise ValueError("There is nothing to initialize in non-GLE thermostats") ssimul = simul.ensemble.thermostat.s if v.mode == "manual": sinput = v.value.copy() if (sinput.size() != ssimul.size() ): raise ValueError("Size mismatch in thermostat initialization data") sinput.shape = ssimul.shape elif v.mode == "chk": rthermo = init_chk(v.value)[2] if not hasattr(rthermo,"s"): raise ValueError("Checkpoint file does not contain usable thermostat data") sinput = rthermo.s.copy() if sinput.shape != ssimul.shape : raise ValueError("Shape mismatch in thermostat initialization data") # if all the preliminary checks are good, we can initialize the s's ssimul[:] = sinput