tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'iss2345' of github.com:eagunn/lammps into iss2345

Browse Source
This commit is contained in:
Axel Kohlmeyer
2020-09-15 14:04:16 -04:00
parent 93a0d4c096 6315b277c8
commit 5cbb54b844
20 changed files with 245 additions and 374 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
doc/src/fix_qeq_reax.rst
View File

@ -20,7 +20,14 @@ Syntax
* cutlo,cuthi = lo and hi cutoff for Taper radius
* tolerance = precision to which charges will be equilibrated
* params = reax/c or a filename
* args = *dual* (optional)
* one or more keywords or keyword/value pairs may be appended
.. parsed-literal::
keyword = *dual* or *maxiter*
*dual* = process S and T matrix in parallel (only for qeq/reax/omp)
*maxiter* N = limit the number of iterations to *N*
Examples
""""""""
@ -28,7 +35,7 @@ Examples
.. code-block:: LAMMPS
fix 1 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c
fix 1 all qeq/reax 1 0.0 10.0 1.0e-6 param.qeq
fix 1 all qeq/reax 1 0.0 10.0 1.0e-6 param.qeq maxiter 500
Description
"""""""""""
@ -69,6 +76,9 @@ The optional *dual* keyword allows to perform the optimization
of the S and T matrices in parallel. This is only supported for
the *qeq/reax/omp* style. Otherwise they are processed separately.
The optional *maxiter* keyword allows changing the max number
of iterations in the linear solver. The default value is 200.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
@ -102,7 +112,7 @@ Related commands
Default
"""""""
none
maxiter 200
----------

8
lib/awpmd/ivutils/include/logexc.h
View File

@ -129,7 +129,7 @@ public:
message_logger(const string &descriptor_="", int out_level=vblALLBAD|vblMESS1,
int stop_level=vblFATAL, int throw_exceptions=0, int use_globally=0)
:descriptor(descriptor_),prev(NULL),next(NULL){
:descriptor(descriptor_),prev(nullptr),next(nullptr){
set_throw(throw_exceptions);
set_levels(out_level,stop_level);
extra_levels(0,0);
@ -157,8 +157,8 @@ public:
return -1;
glogp=prev;
if(glogp)
glogp->next=NULL;
prev=NULL;
glogp->next=nullptr;
prev=nullptr;
}
return 1;
}
@ -244,7 +244,7 @@ public:
FILE *out=stdout, FILE *err=stderr,
int out_level=vblALLBAD|vblMESS1,int stop_level=vblFATAL,
int use_globally=0)
: message_logger(descriptor_,out_level,stop_level,throw_exceptions,use_globally),fout(NULL), ferr(NULL){
: message_logger(descriptor_,out_level,stop_level,throw_exceptions,use_globally),fout(nullptr), ferr(nullptr){
set_out(out);
set_err(err);
}

28
lib/awpmd/ivutils/include/pairhash.h
View File

@ -146,7 +146,7 @@ public:
public:
iterator(const iterator &other):ptr(other.ptr),incr(other.incr){
}
iterator():ptr(NULL),incr(0){}
iterator():ptr(nullptr),incr(0){}
iterator &operator++(){ // prefix
ptr+=incr;
return *this;
@ -173,13 +173,13 @@ public:
size_t sizex, sizey;
//e default constructor
recmatrix(): parr(NULL,1) {
recmatrix(): parr(nullptr,1) {
sizey=sizex=0;
arr=NULL;
arr=nullptr;
}
//e copy constructor: makes a managed copy
recmatrix(const recmatrix &other):sizex(0),sizey(0),arr(NULL){
recmatrix(const recmatrix &other):sizex(0),sizey(0),arr(nullptr){
*this=other;
}
@ -222,7 +222,7 @@ public:
virtual int init(size_t nx, size_t ny, int smanaged=-1){
int managed=parr.managed();
if(managed && (sizex!=nx || sizey!=ny)){
parr.reset(NULL,0);
parr.reset(nullptr,0);
}
if(smanaged>=0){ // for changing the managed flag?
parr.reset(parr.ptr(),smanaged ? smanaged|0x8 : 0 );
@ -355,7 +355,7 @@ public:
public:
iterator(const iterator &other):ptr(other.ptr),incr(other.incr){
}
iterator():ptr(NULL),incr(0){}
iterator():ptr(nullptr),incr(0){}
iterator &operator++(){ // prefix
ptr+=incr;
return *this;
@ -382,13 +382,13 @@ public:
size_t size;
//e default constructor
sqmatrix(): parr(NULL,1) {
sqmatrix(): parr(nullptr,1) {
size=0;
arr=NULL;
arr=nullptr;
}
//e copy constructor: makes a managed copy
sqmatrix(const sqmatrix &other):size(0),arr(NULL){
sqmatrix(const sqmatrix &other):size(0),arr(nullptr){
*this=other;
}
@ -430,7 +430,7 @@ public:
virtual int init(size_t n, int smanaged=-1){
int managed=parr.managed();
if(managed && size!=n){
parr.reset(NULL,0);
parr.reset(nullptr,0);
}
if(smanaged>=0){ // for changing the managed flag?
parr.reset(parr.ptr(),smanaged ? smanaged|0x8 : 0 );
@ -600,9 +600,9 @@ class PairHash{
public:
//e find the value with indexes i, j
//e @return 0 if not found, 1 otherwise
//e if retval is not NULL, puts the found value there
virtual int Find(long i, long j, T *retval=NULL)=0;
virtual int Find(long i, long j, T **retval=NULL)=0;
//e if retval is not a null pointer, puts the found value there
virtual int Find(long i, long j, T *retval=nullptr)=0;
virtual int Find(long i, long j, T **retval=nullptr)=0;
virtual int Del(long i, long j)=0;
virtual int Put(long i, long j, const T *value)=0;
virtual int Put(long i, long j, const T& value)=0;
@ -621,7 +621,7 @@ public:
indm.Set(-1);
arr= new T[n*(n+1)/2];
}
int Find(long i, long j, T *retval=NULL){
int Find(long i, long j, T *retval=nullptr){
long ind=indm(i,j);
if(ind>=0){
if(retval){

20
lib/awpmd/ivutils/include/refobj.h
View File

@ -158,7 +158,7 @@ public:
using base_t::second;
using base_t::first;
mngptr(T* ptr=NULL, int managed=0): pair<T*,int>(ptr,managed){
mngptr(T* ptr=nullptr, int managed=0): pair<T*,int>(ptr,managed){
//if(managed==2)ptr= new T(*ptr);
}
mngptr(const mngarg<T> &arg): pair<T*,int>(arg.first,arg.second){}
@ -166,7 +166,7 @@ public:
reset(arg.first,arg.second);
return *this;
}
void reset(T* ptr=NULL, int managed=0){
void reset(T* ptr=nullptr, int managed=0){
if(second && first && first!=ptr){
if(second&0x8)delete [] first;
else delete first;
@ -313,7 +313,7 @@ template<class T, class delete_t=delete_ptr<T> >
class shptr{
template<class Y, class Z> friend class shptr;
T *p;
int *num; //if num==NULL than p is not managed (as in mngptr)
int *num; //if num==nullptr than p is not managed (as in mngptr)
void set(T *p_, int managed){
p=p_;
@ -321,7 +321,7 @@ class shptr{
num=new int;
*num=1;
}
else num=NULL;
else num=nullptr;
}
template<class Y>
void set(const Y &other){
@ -330,7 +330,7 @@ class shptr{
num=other.num;
if(num)(*num)++;
}
else num=NULL;
else num=nullptr;
}
void set(const shptr &other){
p=other.p;
@ -338,11 +338,11 @@ class shptr{
num=other.num;
if(num)(*num)++;
}
else num=NULL;
else num=nullptr;
}
public:
shptr(T* p=NULL, int managed=1){
shptr(T* p=nullptr, int managed=1){
set(p,managed);
}
shptr(const mngarg<T> &arg){
@ -398,14 +398,14 @@ public:
delete_t()(p);
delete num;
}
num=NULL;
num=nullptr;
}
p=NULL;
p=nullptr;
}
}
bool valid() const {
return p!=NULL;
return p!=nullptr;
}
T* ptr() const {

14
lib/awpmd/ivutils/include/vector_3.h
View File

@ -374,7 +374,7 @@ struct Vector_Nt {
}
T maxcoord(int *ind=NULL) const {
T maxcoord(int *ind=nullptr) const {
int im=0;
T vv=v[0];
for (int i=1; i<N; i++) {
@ -389,7 +389,7 @@ struct Vector_Nt {
//e returns the corrd having maximal absolute value
T maxabscoord(int *ind=NULL) const {
T maxabscoord(int *ind=nullptr) const {
int im=0;
T vv=fabs(v[0]);
for (int i=1; i<N; i++) {
@ -403,7 +403,7 @@ struct Vector_Nt {
}
//e returns the corrd having minimal absolute value
T minabscoord(int *ind=NULL) const {
T minabscoord(int *ind=nullptr) const {
int im=0;
T vv=fabs(v[0]);
for (int i=1; i<N; i++) {
@ -417,7 +417,7 @@ struct Vector_Nt {
}
T mincoord(int *ind=NULL) const {
T mincoord(int *ind=nullptr) const {
int im=0;
T vv=v[0];
for (int i=1; i<N; i++) {
@ -485,7 +485,7 @@ vec_type dist_av(Vector_3 *va1,Vector_3 *va2,int n);
//e finds the average difference norm between two vector sets of the same length
/*e optionally gives the indexes for maximal and minimal difference
va2 can be NULL, then the norm of va1 is used */
va2 can be nullptr, then the norm of va1 is used */
vec_type diff_av(Vector_3 *va1,Vector_3 *va2,int n, int *minind=0, int *maxind=0);
@ -615,9 +615,9 @@ inline Vector_3 randdir(){
///\en Calculates extent of the vector container.
/// \return the center of the vector set, optionally
/// (if arguments are not NULL) fills the bounding box in \a box_min, \a box_max.
/// (if arguments are not null pointers) fills the bounding box in \a box_min, \a box_max.
template<class vec_inp_it>
Vector_3 get_extent(vec_inp_it beg,vec_inp_it end, Vector_3* box_min=NULL,Vector_3* box_max=NULL){
Vector_3 get_extent(vec_inp_it beg,vec_inp_it end, Vector_3* box_min=nullptr,Vector_3* box_max=nullptr){
if(beg==end)
return Vector_3();
Vector_3 center(*beg++);

2
lib/awpmd/ivutils/src/logexc.cpp
View File

@ -13,7 +13,7 @@ message_logger &message_logger::global(){
return *glogp;
}
message_logger *message_logger::glogp=NULL;
message_logger *message_logger::glogp=nullptr;
stdfile_logger default_log("",0,stdout,stderr,vblALLBAD|vblMESS1,vblFATAL,1);
const char *logfmt(const char *format,...){

20
lib/awpmd/systems/interact/TCP/wpmd.h
View File

@ -255,7 +255,7 @@ public:
}
//e Create NormDeriv object and calculate the derivatived for the given WP
void set2(const WavePacket& w2_, const cdouble *I0_=NULL){
void set2(const WavePacket& w2_, const cdouble *I0_=nullptr){
w2=w2_;
d2.set(w2);
w12=conj(w1)*w2;
@ -558,13 +558,13 @@ public:
//e if PBCs are used, the corresponding coordinates of electrons and ions
//e in periodic directions must be within the range [0, cell[per_dir])
//e @returns 1 if OK
int set_pbc(const Vector_3P pcell=NULL, int pbc_=0x7);
int set_pbc(const Vector_3P pcell=nullptr, int pbc_=0x7);
///\en Setup electrons: forms internal wave packet representations.
/// If PBCs are used the coords must be within a range [0, cell).
/// Default electron mass is AWPMD::me.
/// Default (q=NULL )electron charges are -1.
int set_electrons(int spin, int n, Vector_3P x, Vector_3P v, double* w, double* pw, double mass=-1, double *q=NULL);
/// Default (q=nullptr )electron charges are -1.
int set_electrons(int spin, int n, Vector_3P x, Vector_3P v, double* w, double* pw, double mass=-1, double *q=nullptr);
//e setup ion charges and coordinates
//e if PBCs are used the coords must be within a range [0, cell)
@ -593,16 +593,16 @@ public:
// 0x2 -- add ion forces to the existing set
// 0x4 -- calculate derivatives for electronic time step (NOT IMPLEMENTED)
//e if PBCs are used the coords must be within a range [0, cell)
virtual int interaction(int flag=0, Vector_3P fi=NULL, Vector_3P fe_x=NULL,
Vector_3P fe_p=NULL, double *fe_w=NULL, double *fe_pw=NULL, Vector_2P fe_c=NULL);
virtual int interaction(int flag=0, Vector_3P fi=nullptr, Vector_3P fe_x=nullptr,
Vector_3P fe_p=nullptr, double *fe_w=nullptr, double *fe_pw=nullptr, Vector_2P fe_c=nullptr);
//e same as interaction, but using Hartee factorization (no antisymmetrization)
virtual int interaction_hartree(int flag=0, Vector_3P fi=NULL, Vector_3P fe_x=NULL,
Vector_3P fe_p=NULL, double *fe_w=NULL, double *fe_pw=NULL, Vector_2P fe_c=NULL);
virtual int interaction_hartree(int flag=0, Vector_3P fi=nullptr, Vector_3P fe_x=nullptr,
Vector_3P fe_p=nullptr, double *fe_w=nullptr, double *fe_pw=nullptr, Vector_2P fe_c=nullptr);
///\en Calculates ion-ion interactions and updates Eii and ion forces if requested. This function
/// is called form intaraction() and interaction_hartree if calc_ii is set.
virtual int interaction_ii(int flag,Vector_3P fi=NULL);
virtual int interaction_ii(int flag,Vector_3P fi=nullptr);
//e Calculates Norm matrix
//e The result is saved in AWPMD::Norm[s]
@ -643,7 +643,7 @@ public:
///\en Prepares force arrays according to \a flag setting for interaction()
virtual void clear_forces(int flagi,Vector_3P fi, Vector_3P fe_x,
Vector_3P fe_p, double *fe_w, double *fe_pw, Vector_2P fe_c=NULL);
Vector_3P fe_p, double *fe_w, double *fe_pw, Vector_2P fe_c=nullptr);
///\en Creates wave packet according to the given physical parameters.

4
lib/awpmd/systems/interact/TCP/wpmd_split.cpp
View File

@ -312,7 +312,7 @@ void AWPMD_split::clear_forces(int flag,Vector_3P fi, Vector_3P fe_x,
void AWPMD_split::get_el_forces(int flag, Vector_3P fe_x,
Vector_3P fe_p, double *fe_w, double *fe_pw, Vector_2P fe_c){
if(flag&0x4) //need to replace the forces
clear_forces(0x4,NULL,fe_x,fe_p,fe_w,fe_pw,fe_c);
clear_forces(0x4,nullptr,fe_x,fe_p,fe_w,fe_pw,fe_c);
// recalculating derivatives
if(flag&(0x8|0x4)){ //electron forces needed
@ -622,7 +622,7 @@ void AWPMD_split::y_deriv(cdouble v,int s,int c2, int c1){
/// 0x4 -- calculate electronic forces \n
/// 0x8 -- add electronic forces to the existing arrays \n
/// 0x10 -- calculate internal electronic derivatives only: \n
/// will not update electronic force arrays, which may be NULL, \n
/// will not update electronic force arrays, which may be null pointers, \n
/// the forces may be obtained then using \ref get_el_forces() for all WPs \n
/// or separately for each WP using \ref get_wp_force()
/// if PBCs are used the coords must be within a range [0, cell)

16
lib/awpmd/systems/interact/TCP/wpmd_split.h
View File

@ -116,8 +116,8 @@ public:
/// \a n is the number of electrons of a given spin component
/// Electron velocity v is multiplied by mass to obtain momentum.
/// Default mass (-1) means me.
/// Electronic charges q are -1 by default (when q=NULL), otherwise the charges are assigned for each split
int set_electrons(int s, int nel, Vector_3P x, Vector_3P v, double* w, double* pw, Vector_2 *c, int *splits, double mass=-1, double *q=NULL);
/// Electronic charges q are -1 by default (when q=nullptr), otherwise the charges are assigned for each split
int set_electrons(int s, int nel, Vector_3P x, Vector_3P v, double* w, double* pw, Vector_2 *c, int *splits, double mass=-1, double *q=nullptr);
///\en Starts adding new electron: continue with \ref add_split functions.
@ -141,7 +141,7 @@ public:
///\en gets current electronic coordinates, and (optionally) number of wave packets for each electron
int get_electrons(int spin, Vector_3P x, Vector_3P v, double* w, double* pw, cdouble *c, int *splits=NULL, double mass=-1);
int get_electrons(int spin, Vector_3P x, Vector_3P v, double* w, double* pw, cdouble *c, int *splits=nullptr, double mass=-1);
void eterm_deriv(int ic1,int s1, int c1,int k1,int ic2,int s2, int c2,int j2,cdouble pref,
@ -164,8 +164,8 @@ public:
cdouble overlap(int ic1, int s1, int c1,int ic2, int s2, int c2);
//e same as interaction, but using Hartee factorization (no antisymmetrization)
int interaction_hartree(int flag=0, Vector_3P fi=NULL, Vector_3P fe_x=NULL,
Vector_3P fe_p=NULL, double *fe_w=NULL, double *fe_pw=NULL, Vector_2P fe_c=NULL);
int interaction_hartree(int flag=0, Vector_3P fi=nullptr, Vector_3P fe_x=nullptr,
Vector_3P fe_p=nullptr, double *fe_w=nullptr, double *fe_pw=nullptr, Vector_2P fe_c=nullptr);
///\en Calculates interaction in the system of ni ions + electrons
/// the electonic subsystem must be previously setup by set_electrons, ionic by set_ions
@ -174,12 +174,12 @@ public:
/// 0x4 -- calculate electronic forces \n
/// 0x8 -- add electronic forces to the existing arrays \n
/// 0x10 -- calculate internal electronic derivatives only: \n
/// will not update electronic force arrays, which may be NULL, \n
/// will not update electronic force arrays, which may be null pointers, \n
/// the forces may be obtained then using \ref get_el_forces() for all WPs \n
/// or separately for each WP using \ref get_wp_force()
/// if PBCs are used the coords must be within a range [0, cell)
int interaction(int flag=0, Vector_3P fi=NULL, Vector_3P fe_x=NULL,
Vector_3P fe_p=NULL, double *fe_w=NULL, double *fe_pw=NULL, Vector_2P fe_c=NULL);
int interaction(int flag=0, Vector_3P fi=nullptr, Vector_3P fe_x=nullptr,
Vector_3P fe_p=nullptr, double *fe_w=nullptr, double *fe_pw=nullptr, Vector_2P fe_c=nullptr);
///\en Get electronic forcess in the arrays provided, using calculated internal representation
/// Valid flag settings are:\n

30
lib/poems/SystemProcessor.h
View File

@ -110,21 +110,21 @@ void SystemProcessor::processArray(int** links, int numLinks)
do
{
currentNode = findSingleLink(temp); //find the start of the next available chain
if(currentNode != NULL)
if(currentNode != nullptr)
{
headsOfSystems.Append(AddNewChain(currentNode)); //and add it to the headsOfSystems list of chains
}
}
while(currentNode != NULL); //repeat this until all chains have been added
while(currentNode != nullptr); //repeat this until all chains have been added
}
POEMSChain * SystemProcessor::AddNewChain(POEMSNode * currentNode){
if(currentNode == NULL) //Termination condition; if the currentNode is null, then return null
if(currentNode == nullptr) //Termination condition; if the currentNode is null, then return null
{
return NULL;
return nullptr;
}
int * tmp;
POEMSNode * nextNode = NULL; //nextNode stores the proposed next node to add to the chain. this will be checked to make sure no backtracking is occurring before being assigned as the current node.
POEMSNode * nextNode = nullptr; //nextNode stores the proposed next node to add to the chain. this will be checked to make sure no backtracking is occurring before being assigned as the current node.
POEMSChain * newChain = new POEMSChain; //make a new POEMSChain object. This will be the object returned
if(currentNode->links.GetNumElements() == 0) //if we have no links from this node, then the whole chain is only one node. Add this node to the chain and return it; mark node as visited for future reference
@ -168,8 +168,8 @@ POEMSChain * SystemProcessor::AddNewChain(POEMSNode * currentNode){
newChain->listOfNodes.Append(tmp); //append the last node before branch (node shared jointly with branch chains)
//re-mark as visited, just to make sure
ListElement<POEMSNode> * tempNode = currentNode->links.GetHeadElement(); //go through all of the links, one at a time that branch
POEMSChain * tempChain = NULL; //temporary variable to hold data
while(tempNode != NULL) //when we have followed all links, stop
POEMSChain * tempChain = nullptr; //temporary variable to hold data
while(tempNode != nullptr) //when we have followed all links, stop
{
if(setLinkVisited(tempNode->value, currentNode)) //dont backtrack, or create closed loops
{
@ -187,12 +187,12 @@ POEMSNode * SystemProcessor::findSingleLink(TreeNode * aNode)
//This function takes the root of a search tree containing POEMSNodes and returns a POEMSNode corresponding to the start of a chain in the
//system. It finds a node that has not been visited before, and only has one link; this node will be used as the head of the chain.
{
if(aNode == NULL)
if(aNode == nullptr)
{
return NULL;
return nullptr;
}
POEMSNode * returnVal = (POEMSNode *)aNode->GetAuxData(); //get the poemsnode data out of the treenode
POEMSNode * detectLoneLoops = NULL; //is used to handle a loop that has no protruding chains
POEMSNode * detectLoneLoops = nullptr; //is used to handle a loop that has no protruding chains
if(returnVal->visited == false)
{
detectLoneLoops = returnVal; //if we find any node that has not been visited yet, save it
@ -202,15 +202,15 @@ POEMSNode * SystemProcessor::findSingleLink(TreeNode * aNode)
return returnVal; //return the node is it meets this criteria
}
returnVal = findSingleLink(aNode->Left()); //otherwise, check the left subtree
if(returnVal == NULL) //and if we find nothing...
if(returnVal == nullptr) //and if we find nothing...
{
returnVal = findSingleLink(aNode->Right()); //check the right subtree
}
if(returnVal == NULL) //if we could not find any chains
if(returnVal == nullptr) //if we could not find any chains
{
returnVal = detectLoneLoops; //see if we found any nodes at all that havent been processed
}
return returnVal; //return what we find (will be NULL if no new chains are
return returnVal; //return what we find (will be a null pointer if no new chains are
//found)
}
@ -226,7 +226,7 @@ bool SystemProcessor::setLinkVisited(POEMSNode * firstNode, POEMSNode * secondNo
//cout << "Checking link between nodes " << firstNode->idNumber << " and " << secondNode->idNumber << "... ";
ListElement<POEMSNode> * tmp = firstNode->links.GetHeadElement(); //get the head element of the list of pointers for node 1
ListElement<bool> * tmp2 = firstNode->taken.GetHeadElement(); //get the head element of the list of bool isVisited flags for node 1
while(tmp->value != NULL || tmp2->value != NULL) //go through until we reach the end of the lists
while(tmp->value != nullptr || tmp2->value != nullptr) //go through until we reach the end of the lists
{
if(tmp->value == secondNode) //if we find the link to the other node
{
@ -248,7 +248,7 @@ bool SystemProcessor::setLinkVisited(POEMSNode * firstNode, POEMSNode * secondNo
tmp = secondNode->links.GetHeadElement(); //now, if the link was unvisited, we need to go set the other node's list such that
//it also knows this link is being visited
tmp2 = secondNode->taken.GetHeadElement();
while(tmp->value != NULL || tmp2->value != NULL) //go through the list
while(tmp->value != nullptr || tmp2->value != nullptr) //go through the list
{
if(tmp->value == firstNode) //if we find the link
{

20
lib/poems/poemsnodelib.h
View File

@ -23,7 +23,7 @@
using namespace std;
TreeNode *GetTreeNode(int item,TreeNode *lptr = NULL,TreeNode *rptr =NULL);
TreeNode *GetTreeNode(int item,TreeNode *lptr = nullptr,TreeNode *rptr =nullptr);
void FreeTreeNode(TreeNode *p);
@ -46,7 +46,7 @@ void PrintTree (TreeNode *t, int level);
void Postorder (TreeNode *t, void visit(TreeNode* &t))
{
// the recursive scan terminates on a empty subtree
if (t != NULL)
if (t != nullptr)
{
Postorder(t->Left(), visit); // descend left
Postorder(t->Right(), visit); // descend right
@ -59,7 +59,7 @@ void Postorder (TreeNode *t, void visit(TreeNode* &t))
void Preorder (TreeNode *t, void visit(TreeNode* &t))
{
// the recursive scan terminates on a empty subtree
if (t != NULL)
if (t != nullptr)
{
visit(t); // visit the node
Preorder(t->Left(), visit); // descend left
@ -69,7 +69,7 @@ void Preorder (TreeNode *t, void visit(TreeNode* &t))
//create TreeNode object with pointer fields lptr and rptr
// The pointers have default value NULL
// The pointers have default value nullptr
TreeNode *GetTreeNode(int item,TreeNode *lptr,TreeNode *rptr)
{
TreeNode *p;
@ -79,7 +79,7 @@ TreeNode *GetTreeNode(int item,TreeNode *lptr,TreeNode *rptr)
p = new TreeNode(item, lptr, rptr);
// if insufficient memory, terminatewith an error message
if (p == NULL)
if (p == nullptr)
{
cerr << "Memory allocation failure!\n";
exit(1);
@ -103,14 +103,14 @@ void FreeTreeNode(TreeNode *p)
void CountLeaf (TreeNode *t, int& count)
{
//use postorder descent
if(t !=NULL)
if(t !=nullptr)
{
CountLeaf(t->Left(), count); // descend left
CountLeaf(t->Right(), count); // descend right
// check if node t is a leaf node (no descendants)
// if so, increment the variable count
if (t->Left() == NULL && t->Right() == NULL)
if (t->Left() == nullptr && t->Right() == nullptr)
count++;
}
}
@ -124,7 +124,7 @@ int Depth (TreeNode *t)
{
int depthLeft, depthRight, depthval;
if (t == NULL)
if (t == nullptr)
depthval = -1;
else
{
@ -145,8 +145,8 @@ void IndentBlanks(int num)
void PrintTree (TreeNode *t, int level)
{
//print tree with root t, as long as t!=NULL
if (t != NULL)
//print tree with root t, as long as t!=nullptr
if (t != nullptr)
{
int indentUnit = 5;
// print right branch of tree t

64
lib/poems/poemstree.h
View File

@ -82,7 +82,7 @@ public:
DeleteAuxData = callback;
}
void Insert(const int& item, const int& data, void * AuxData = NULL);
void Insert(const int& item, const int& data, void * AuxData = nullptr);
void Delete(const int& item);
void AVLInsert(TreeNode* &tree, TreeNode* newNode, int &reviseBalanceFactor);
void ClearList(void);
@ -99,7 +99,7 @@ Tree::Tree(void)
root = 0;
current = 0;
size = 0;
DeleteAuxData = NULL;
DeleteAuxData = nullptr;
}
@ -131,18 +131,18 @@ Tree& Tree::operator = (const Tree& rhs)
}
// search for data item in the tree. if found, return its node
// address and a pointer to its parent; otherwise, return NULL
// address and a pointer to its parent; otherwise, return a null pointer
TreeNode *Tree::FindNode(const int& item,
TreeNode* & parent) const
{
// cycle t through the tree starting with root
TreeNode *t = root;
// the parent of the root is NULL
parent = NULL;
// the parent of the root is a null pointer
parent = nullptr;
// terminate on empty subtree
while(t != NULL)
while(t != nullptr)
{
// stop on a match
if (item == t->data)
@ -158,7 +158,7 @@ TreeNode *Tree::FindNode(const int& item,
}
}
// return pointer to node; NULL if not found
// return pointer to node; a null pointer if not found
return t;
}
@ -172,14 +172,14 @@ void * Tree::Find(int& item)
current = FindNode (item, parent);
// if item found, assign data to item and return True
if (current != NULL)
if (current != nullptr)
{
item = current->data;
return current->GetAuxData();
}
else
// item not found in the tree. return False
return NULL;
// item not found in the tree. return a null pointer
return nullptr;
}
@ -194,7 +194,7 @@ void Tree::Insert(const int& item, const int& data, void * AuxData)
int reviseBalanceFactor = 0;
// get a new AVL tree node with empty pointer fields
newNode = GetTreeNode(item,NULL,NULL);
newNode = GetTreeNode(item,nullptr,nullptr);
newNode->data = data;
newNode->SetAuxData(AuxData);
// call recursive routine to actually insert the element
@ -213,7 +213,7 @@ void Tree::AVLInsert(TreeNode *&tree, TreeNode *newNode, int &reviseBalanceFacto
int rebalanceCurrNode;
// scan reaches an empty tree; time to insert the new node
if (tree == NULL)
if (tree == nullptr)
{
// update the parent to point at newNode
tree = newNode;
@ -456,16 +456,16 @@ void Tree::Delete(const int& item)
// search for a node containing data value item. obtain its
// node address and that of its parent
if ((DNodePtr = FindNode (item, PNodePtr)) == NULL)
if ((DNodePtr = FindNode (item, PNodePtr)) == nullptr)
return;
// If D has NULL pointer, the
// If D has null pointer, the
// replacement node is the one on the other branch
if (DNodePtr->right == NULL)
if (DNodePtr->right == nullptr)
RNodePtr = DNodePtr->left;
else if (DNodePtr->left == NULL)
else if (DNodePtr->left == nullptr)
RNodePtr = DNodePtr->right;
// Both pointers of DNodePtr are non-NULL
// Both pointers of DNodePtr are non-null pointers
else
{
// Find and unlink replacement node for D
@ -483,7 +483,7 @@ void Tree::Delete(const int& item)
// descend down right subtree of the left child of D
// keeping a record of current node and its parent.
// when we stop, we have found the replacement
while (RNodePtr->right != NULL)
while (RNodePtr->right != nullptr)
{
PofRNodePtr = RNodePtr;
RNodePtr = RNodePtr;
@ -508,7 +508,7 @@ void Tree::Delete(const int& item)
// complete the link to the parent node
// deleting the root node. assign new root
if (PNodePtr == NULL)
if (PNodePtr == nullptr)
root = RNodePtr;
// attach R to the correct branch of P
else if (DNodePtr->data < PNodePtr->data)
@ -529,7 +529,7 @@ void Tree::Delete(const int& item)
// assign node value to item; otherwise, insert item in tree
void Tree::Update(const int& item)
{
if (current !=NULL && current->data == item)
if (current !=nullptr && current->data == item)
current->data = item;
else
Insert(item, item);
@ -545,25 +545,25 @@ TreeNode *Tree::CopyTree(TreeNode *t)
TreeNode *newlptr, *newrptr, *newnode;
// stop the recursive scan when we arrive at an empty tree
if (t == NULL)
return NULL;
if (t == nullptr)
return nullptr;
// CopyTree builds a new tree by scanning the nodes of t.
// At each node in t, CopyTree checks for a left child. if
// present it makes a copy of left child or returns NULL.
// present it makes a copy of left child or returns a null pointer.
// the algorithm similarly checks for a right child.
// CopyTree builds a copy of node using GetTreeNode and
// appends copy of left and right children to node.
if (t->Left() !=NULL)
if (t->Left() !=nullptr)
newlptr = CopyTree(t->Left());
else
newlptr = NULL;
newlptr = nullptr;
if (t->Right() !=NULL)
if (t->Right() !=nullptr)
newrptr = CopyTree(t->Right());
else
newrptr = NULL;
newrptr = nullptr;
// Build new tree from the bottom up by building the two
@ -579,12 +579,12 @@ TreeNode *Tree::CopyTree(TreeNode *t)
// the tree and delete each node at the visit operation
void Tree::DeleteTree(TreeNode *t)
{
if (t != NULL) {
if (t != nullptr) {
DeleteTree(t->Left());
DeleteTree(t->Right());
void *aux = t->GetAuxData();
if (aux != NULL) {
if (DeleteAuxData != NULL) {
if (aux != nullptr) {
if (DeleteAuxData != nullptr) {
(*DeleteAuxData)(aux);
} else {
delete (TreeNode *) aux;
@ -595,11 +595,11 @@ void Tree::DeleteTree(TreeNode *t)
}
// call the function DeleteTree to deallocate the nodes. then
// set the root pointer back to NULL
// set the root pointer back to a null pointer
void Tree::ClearTree(TreeNode * &t)
{
DeleteTree(t);
t = NULL; // root now NULL
t = nullptr; // root now a null pointer
}
// delete all nodes in list

2
lib/poems/poemstreenode.cpp
View File

@ -18,7 +18,7 @@
#include "poemstreenode.h"
// constructor; initialize the data and pointer fields
// The pointer value NULL assigns a empty subtree
// A null pointer value assigns a empty subtree
TreeNode::TreeNode (const int & item, TreeNode *lptr,TreeNode *rptr,
int balfac):data(item), left(lptr), right(rptr), balanceFactor(balfac)
{

2
lib/poems/solver.cpp
View File

@ -45,7 +45,7 @@ Solver * Solver::GetSolver(SolverType solverToMake) //returning a pointer to a n
switch((int)solverToMake)
{
case ONSOLVER: return new OnSolver();
default: return NULL;
default: return nullptr;
}
}

10
lib/poems/system.cpp
View File

@ -23,7 +23,7 @@
System::System(){
mappings = NULL;
mappings = nullptr;
}
System::~System(){
@ -238,7 +238,7 @@ void System::Create_DegenerateSystem(int& nfree, int*freelist, double *&masstota
//-------------------------------------------------------------------------//
// Declaring Temporary Entities
//-------------------------------------------------------------------------//
Body* body = NULL;
Body* body = nullptr;
Body* prev;
Body* Inertial;
Point* origin;
@ -246,7 +246,7 @@ void System::Create_DegenerateSystem(int& nfree, int*freelist, double *&masstota
Point* point_CM;
Point* point_p;
Point* point_k;
Point* point_ch = NULL;
Point* point_ch = nullptr;
Vect3 r1,r2,r3,v1,v2,v3;
Mat3x3 IM, N, PKCK,PKCN ;
ColMatrix qo, uo, q, qdot,w;
@ -391,7 +391,7 @@ void System::Create_System_LAMMPS(int numbodies, double *mass,double **inertia,
// Declaring Temporary Entities
//-------------------------------------------------------------------------//
Body* body = NULL;
Body* body = nullptr;
Body* prev;
Body* Inertial;
Point* origin;
@ -399,7 +399,7 @@ void System::Create_System_LAMMPS(int numbodies, double *mass,double **inertia,
Point* point_CM;
Point* point_p;
Point* point_k;
Point* point_ch = NULL;
Point* point_ch = nullptr;
Vect3 r1,r2,r3,v1,v2,v3;
Mat3x3 IM, N, PKCK,PKCN ;
ColMatrix qo, uo, q, qdot,w;

6
lib/poems/workspace.cpp
View File

@ -52,7 +52,7 @@ void Workspace::allocateNewSystem() {
Workspace::Workspace(){
currentIndex = -1;
maxAlloc = 0;
system = NULL;
system = nullptr;
}
Workspace::~Workspace(){
@ -133,7 +133,7 @@ if(njoint){
int ttk = 0;
while(NodeValue != NULL) {
while(NodeValue != nullptr) {
array = new int[NodeValue->value->listOfNodes.GetNumElements()];
arrayFromChain = NodeValue->value->listOfNodes.CreateArray();
numElementsInSystem = NodeValue->value->listOfNodes.GetNumElements();
@ -200,7 +200,7 @@ System* Workspace::GetSystem(int index){
}
}
else{
return NULL;
return nullptr;
}
}

31
src/KOKKOS/fix_qeq_reax_kokkos.cpp
View File

@ -112,7 +112,7 @@ void FixQEqReaxKokkos<DeviceType>::init()
("FixQEqReax::params",ntypes+1);
params = k_params.template view<DeviceType>();
for (n = 1; n <= ntypes; n++) {
for (int n = 1; n <= ntypes; n++) {
k_params.h_view(n).chi = chi[n];
k_params.h_view(n).eta = eta[n];
k_params.h_view(n).gamma = gamma[n];
@ -351,34 +351,35 @@ void FixQEqReaxKokkos<DeviceType>::allocate_array()
if (atom->nmax > nmax) {
nmax = atom->nmax;
k_o = DAT::tdual_ffloat_1d("qeq/kk:h_o",nmax);
k_o = DAT::tdual_ffloat_1d("qeq/kk:o",nmax);
d_o = k_o.template view<DeviceType>();
h_o = k_o.h_view;
d_Hdia_inv = typename AT::t_ffloat_1d("qeq/kk:h_Hdia_inv",nmax);
d_Hdia_inv = typename AT::t_ffloat_1d("qeq/kk:Hdia_inv",nmax);
d_b_s = typename AT::t_ffloat_1d("qeq/kk:h_b_s",nmax);
d_b_s = typename AT::t_ffloat_1d("qeq/kk:b_s",nmax);
d_b_t = typename AT::t_ffloat_1d("qeq/kk:h_b_t",nmax);
d_b_t = typename AT::t_ffloat_1d("qeq/kk:b_t",nmax);
k_s = DAT::tdual_ffloat_1d("qeq/kk:h_s",nmax);
k_s = DAT::tdual_ffloat_1d("qeq/kk:s",nmax);
d_s = k_s.template view<DeviceType>();
h_s = k_s.h_view;
k_t = DAT::tdual_ffloat_1d("qeq/kk:h_t",nmax);
k_t = DAT::tdual_ffloat_1d("qeq/kk:t",nmax);
d_t = k_t.template view<DeviceType>();
h_t = k_t.h_view;
d_p = typename AT::t_ffloat_1d("qeq/kk:h_p",nmax);
d_p = typename AT::t_ffloat_1d("qeq/kk:p",nmax);
d_r = typename AT::t_ffloat_1d("qeq/kk:h_r",nmax);
d_r = typename AT::t_ffloat_1d("qeq/kk:r",nmax);
k_d = DAT::tdual_ffloat_1d("qeq/kk:h_d",nmax);
k_d = DAT::tdual_ffloat_1d("qeq/kk:d",nmax);
d_d = k_d.template view<DeviceType>();
h_d = k_d.h_view;
}
// init_storage
FixQEqReaxKokkosZeroFunctor<DeviceType> zero_functor(this);
Kokkos::parallel_for(ignum,zero_functor);
@ -779,8 +780,7 @@ void FixQEqReaxKokkos<DeviceType>::cg_solve1()
F_FLOAT sig_new = dot_sqr;
int loop;
const int loopmax = 200;
for (loop = 1; (loop < loopmax) && (sqrt(sig_new)/b_norm > tolerance); loop++) {
for (loop = 1; (loop < imax) && (sqrt(sig_new)/b_norm > tolerance); loop++) {
// comm->forward_comm_fix(this); //Dist_vector( d );
pack_flag = 1;
@ -848,7 +848,7 @@ void FixQEqReaxKokkos<DeviceType>::cg_solve1()
Kokkos::parallel_for(inum,vecsum2_functor);
}
if (loop >= loopmax && comm->me == 0) {
if (loop >= imax && comm->me == 0) {
char str[128];
sprintf(str,"Fix qeq/reax cg_solve1 convergence failed after %d iterations "
"at " BIGINT_FORMAT " step: %f",loop,update->ntimestep,sqrt(sig_new)/b_norm);
@ -918,8 +918,7 @@ void FixQEqReaxKokkos<DeviceType>::cg_solve2()
F_FLOAT sig_new = dot_sqr;
int loop;
const int loopmax = 200;
for (loop = 1; (loop < loopmax) && (sqrt(sig_new)/b_norm > tolerance); loop++) {
for (loop = 1; (loop < imax) && (sqrt(sig_new)/b_norm > tolerance); loop++) {
// comm->forward_comm_fix(this); //Dist_vector( d );
pack_flag = 1;
@ -987,7 +986,7 @@ void FixQEqReaxKokkos<DeviceType>::cg_solve2()
Kokkos::parallel_for(inum,vecsum2_functor);
}
if (loop >= loopmax && comm->me == 0) {
if (loop >= imax && comm->me == 0) {
char str[128];
sprintf(str,"Fix qeq/reax cg_solve2 convergence failed after %d iterations "
"at " BIGINT_FORMAT " step: %f",loop,update->ntimestep,sqrt(sig_new)/b_norm);

144
src/USER-OMP/fix_qeq_reax_omp.cpp
View File

@ -63,8 +63,6 @@ using namespace FixConst;
FixQEqReaxOMP::FixQEqReaxOMP(LAMMPS *lmp, int narg, char **arg) :
FixQEqReax(lmp, narg, arg)
{
if (narg<8 || narg>9) error->all(FLERR,"Illegal fix qeq/reax/omp command");
b_temp = nullptr;
// ASPC: Kolafa, J. Comp. Chem., 25(3), 335 (2003)
@ -85,6 +83,11 @@ FixQEqReaxOMP::~FixQEqReaxOMP()
void FixQEqReaxOMP::post_constructor()
{
grow_arrays(atom->nmax);
for (int i = 0; i < atom->nmax; i++)
for (int j = 0; j < nprev; ++j)
s_hist[i][j] = t_hist[i][j] = 0;
pertype_parameters(pertype_option);
}
@ -148,7 +151,6 @@ void FixQEqReaxOMP::init()
void FixQEqReaxOMP::compute_H()
{
int inum, *ilist, *numneigh, **firstneigh;
double SMALL = 0.0001;
int *type = atom->type;
@ -156,17 +158,6 @@ void FixQEqReaxOMP::compute_H()
double **x = atom->x;
int *mask = atom->mask;
if (reaxc) {
inum = reaxc->list->inum;
ilist = reaxc->list->ilist;
numneigh = reaxc->list->numneigh;
firstneigh = reaxc->list->firstneigh;
} else {
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
}
int ai, num_nbrs;
// sumscan of the number of neighbors per atom to determine the offsets
@ -175,7 +166,7 @@ void FixQEqReaxOMP::compute_H()
num_nbrs = 0;
for (int itr_i = 0; itr_i < inum; ++itr_i) {
for (int itr_i = 0; itr_i < nn; ++itr_i) {
ai = ilist[itr_i];
H.firstnbr[ai] = num_nbrs;
num_nbrs += numneigh[ai];
@ -197,7 +188,7 @@ void FixQEqReaxOMP::compute_H()
#if defined(_OPENMP)
#pragma omp for schedule(guided)
#endif
for (int ii = 0; ii < inum; ii++) {
for (int ii = 0; ii < nn; ii++) {
int i = ilist[ii];
if (mask[i] & groupbit) {
jlist = firstneigh[i];
@ -214,7 +205,7 @@ void FixQEqReaxOMP::compute_H()
flag = 0;
if (r_sqr <= SQR(swb)) {
if (j < n) flag = 1;
if (j < atom->nlocal) flag = 1;
else if (tag[i] < tag[j]) flag = 1;
else if (tag[i] == tag[j]) {
if (dz > SMALL) flag = 1;
@ -251,11 +242,6 @@ void FixQEqReaxOMP::compute_H()
void FixQEqReaxOMP::init_storage()
{
int NN;
if (reaxc) NN = reaxc->list->inum + reaxc->list->gnum;
else NN = list->inum + list->gnum;
#if defined(_OPENMP)
#pragma omp parallel for schedule(static)
#endif
@ -284,8 +270,21 @@ void FixQEqReaxOMP::pre_force(int /* vflag */)
if (update->ntimestep % nevery) return;
if (comm->me == 0) t_start = MPI_Wtime();
n = atom->nlocal;
N = atom->nlocal + atom->nghost;
int n = atom->nlocal;
if (reaxc) {
nn = reaxc->list->inum;
NN = reaxc->list->inum + reaxc->list->gnum;
ilist = reaxc->list->ilist;
numneigh = reaxc->list->numneigh;
firstneigh = reaxc->list->firstneigh;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
}
// grow arrays if necessary
// need to be atom->nmax in length
@ -365,16 +364,7 @@ void FixQEqReaxOMP::init_matvec()
/* fill-in H matrix */
compute_H();
int nn,i;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
ilist = list->ilist;
}
int i;
// Should really be more careful with initialization and first (aspc_order+2) MD steps
if (do_aspc) {
@ -450,24 +440,12 @@ void FixQEqReaxOMP::init_matvec()
int FixQEqReaxOMP::CG( double *b, double *x)
{
int i, imax;
int i;
double alpha, beta, b_norm;
double sig_old, sig_new;
double my_buf[2], buf[2];
int nn;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
ilist = list->ilist;
}
imax = 200;
pack_flag = 1;
sparse_matvec( &H, x, q );
comm->reverse_comm_fix( this); //Coll_Vector( q );
@ -579,8 +557,7 @@ void FixQEqReaxOMP::sparse_matvec( sparse_matrix *A, double *x, double *b)
#endif
{
int i, j, itr_j;
int nn, NN, ii;
int *ilist;
int ii;
int nthreads = comm->nthreads;
#if defined(_OPENMP)
int tid = omp_get_thread_num();
@ -588,16 +565,6 @@ void FixQEqReaxOMP::sparse_matvec( sparse_matrix *A, double *x, double *b)
int tid = 0;
#endif
if (reaxc) {
nn = reaxc->list->inum;
NN = reaxc->list->inum + reaxc->list->gnum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
}
#if defined(_OPENMP)
#pragma omp for schedule(dynamic,50)
#endif
@ -655,17 +622,6 @@ void FixQEqReaxOMP::calculate_Q()
int i;
double *q = atom->q;
int nn;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
ilist = list->ilist;
}
double tmp1, tmp2;
tmp1 = tmp2 = 0.0;
#if defined(_OPENMP)
@ -718,10 +674,6 @@ void FixQEqReaxOMP::vector_sum( double* dest, double c, double* v,
double d, double* y, int k)
{
int i;
int *ilist;
if (reaxc) ilist = reaxc->list->ilist;
else ilist = list->ilist;
#if defined(_OPENMP)
#pragma omp parallel for schedule(static) private(i)
@ -737,10 +689,6 @@ void FixQEqReaxOMP::vector_sum( double* dest, double c, double* v,
void FixQEqReaxOMP::vector_add( double* dest, double c, double* v, int k)
{
int i;
int *ilist;
if (reaxc) ilist = reaxc->list->ilist;
else ilist = list->ilist;
#if defined(_OPENMP)
#pragma omp parallel for schedule(static) private(i)
@ -765,24 +713,12 @@ int FixQEqReaxOMP::dual_CG( double *b1, double *b2, double *x1, double *x2)
startTimeBase = MPI_Wtime();
#endif
int i, imax;
int i;
double alpha_s, alpha_t, beta_s, beta_t, b_norm_s, b_norm_t;
double sig_old_s, sig_old_t, sig_new_s, sig_new_t;
double my_buf[4], buf[4];
int nn;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
ilist = list->ilist;
}
imax = 200;
pack_flag = 5; // forward 2x d and reverse 2x q
dual_sparse_matvec( &H, x1, x2, q );
comm->reverse_comm_fix(this); //Coll_Vector( q );
@ -975,8 +911,7 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x1, double *x2
#endif
{
int i, j, itr_j;
int nn, NN, ii;
int *ilist;
int ii;
int indxI, indxJ;
int nthreads = comm->nthreads;
@ -986,16 +921,6 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x1, double *x2
int tid = 0;
#endif
if (reaxc) {
nn = reaxc->list->inum;
NN = reaxc->list->inum + reaxc->list->gnum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
}
#if defined(_OPENMP)
#pragma omp for schedule(dynamic,50)
#endif
@ -1077,8 +1002,7 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x, double *b )
#endif
{
int i, j, itr_j;
int nn, NN, ii;
int *ilist;
int ii;
int indxI, indxJ;
int nthreads = comm->nthreads;
@ -1088,16 +1012,6 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x, double *b )
int tid = 0;
#endif
if (reaxc) {
nn = reaxc->list->inum;
NN = reaxc->list->inum + reaxc->list->gnum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
}
#if defined(_OPENMP)
#pragma omp for schedule(dynamic,50)
#endif

178
src/USER-REAXC/fix_qeq_reax.cpp
View File

@ -62,9 +62,7 @@ static const char cite_fix_qeq_reax[] =
FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), pertype_option(nullptr)
{
if (lmp->citeme) lmp->citeme->add(cite_fix_qeq_reax);
if (narg<8 || narg>9) error->all(FLERR,"Illegal fix qeq/reax command");
if (narg<8 || narg>11) error->all(FLERR,"Illegal fix qeq/reax command");
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
if (nevery <= 0) error->all(FLERR,"Illegal fix qeq/reax command");
@ -79,14 +77,23 @@ FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) :
// dual CG support only available for USER-OMP variant
// check for compatibility is in Fix::post_constructor()
dual_enabled = 0;
if (narg == 9) {
if (strcmp(arg[8],"dual") == 0) dual_enabled = 1;
else error->all(FLERR,"Illegal fix qeq/reax command");
imax = 200;
int iarg = 8;
while (iarg < narg) {
if (strcmp(arg[iarg],"dual") == 0) dual_enabled = 1;
else if (strcmp(arg[iarg],"maxiter") == 0) {
if (iarg+1 > narg-1)
error->all(FLERR,"Illegal fix qeq/reax command");
imax = utils::numeric(FLERR,arg[iarg+1],false,lmp);
iarg++;
} else error->all(FLERR,"Illegal fix qeq/reax command");
iarg++;
}
shld = nullptr;
n = n_cap = 0;
N = nmax = 0;
nn = n_cap = 0;
NN = nmax = 0;
m_fill = m_cap = 0;
pack_flag = 0;
s = nullptr;
@ -123,11 +130,7 @@ FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) :
reaxc = (PairReaxC *) force->pair_match("^reax/c",0);
s_hist = t_hist = nullptr;
grow_arrays(atom->nmax);
atom->add_callback(0);
for (int i = 0; i < atom->nmax; i++)
for (int j = 0; j < nprev; ++j)
s_hist[i][j] = t_hist[i][j] = 0;
}
/* ---------------------------------------------------------------------- */
@ -161,6 +164,13 @@ FixQEqReax::~FixQEqReax()
void FixQEqReax::post_constructor()
{
if (lmp->citeme) lmp->citeme->add(cite_fix_qeq_reax);
grow_arrays(atom->nmax);
for (int i = 0; i < atom->nmax; i++)
for (int j = 0; j < nprev; ++j)
s_hist[i][j] = t_hist[i][j] = 0;
pertype_parameters(pertype_option);
if (dual_enabled)
error->all(FLERR,"Dual keyword only supported with fix qeq/reax/omp");
@ -287,8 +297,7 @@ void FixQEqReax::reallocate_storage()
void FixQEqReax::allocate_matrix()
{
int i,ii,inum,m;
int *ilist, *numneigh;
int i,ii,n,m;
int mincap;
double safezone;
@ -306,18 +315,8 @@ void FixQEqReax::allocate_matrix()
// determine the total space for the H matrix
if (reaxc) {
inum = reaxc->list->inum;
ilist = reaxc->list->ilist;
numneigh = reaxc->list->numneigh;
} else {
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
}
m = 0;
for (ii = 0; ii < inum; ii++) {
for (ii = 0; ii < nn; ii++) {
i = ilist[ii];
m += numneigh[i];
}
@ -432,6 +431,20 @@ void FixQEqReax::init_taper()
void FixQEqReax::setup_pre_force(int vflag)
{
if (reaxc) {
nn = reaxc->list->inum;
NN = reaxc->list->inum + reaxc->list->gnum;
ilist = reaxc->list->ilist;
numneigh = reaxc->list->numneigh;
firstneigh = reaxc->list->firstneigh;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
}
deallocate_storage();
allocate_storage();
@ -495,8 +508,21 @@ void FixQEqReax::pre_force(int /*vflag*/)
if (update->ntimestep % nevery) return;
if (comm->me == 0) t_start = MPI_Wtime();
n = atom->nlocal;
N = atom->nlocal + atom->nghost;
int n = atom->nlocal;
if (reaxc) {
nn = reaxc->list->inum;
NN = reaxc->list->inum + reaxc->list->gnum;
ilist = reaxc->list->ilist;
numneigh = reaxc->list->numneigh;
firstneigh = reaxc->list->firstneigh;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
}
// grow arrays if necessary
// need to be atom->nmax in length
@ -540,16 +566,7 @@ void FixQEqReax::init_matvec()
/* fill-in H matrix */
compute_H();
int nn, ii, i;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
ilist = list->ilist;
}
int ii, i;
for (ii = 0; ii < nn; ++ii) {
i = ilist[ii];
@ -578,7 +595,7 @@ void FixQEqReax::init_matvec()
void FixQEqReax::compute_H()
{
int inum, jnum, *ilist, *jlist, *numneigh, **firstneigh;
int jnum;
int i, j, ii, jj, flag;
double dx, dy, dz, r_sqr;
const double SMALL = 0.0001;
@ -588,22 +605,10 @@ void FixQEqReax::compute_H()
double **x = atom->x;
int *mask = atom->mask;
if (reaxc) {
inum = reaxc->list->inum;
ilist = reaxc->list->ilist;
numneigh = reaxc->list->numneigh;
firstneigh = reaxc->list->firstneigh;
} else {
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
}
// fill in the H matrix
m_fill = 0;
r_sqr = 0;
for (ii = 0; ii < inum; ii++) {
for (ii = 0; ii < nn; ii++) {
i = ilist[ii];
if (mask[i] & groupbit) {
jlist = firstneigh[i];
@ -621,7 +626,7 @@ void FixQEqReax::compute_H()
flag = 0;
if (r_sqr <= SQR(swb)) {
if (j < n) flag = 1;
if (j < atom->nlocal) flag = 1;
else if (tag[i] < tag[j]) flag = 1;
else if (tag[i] == tag[j]) {
if (dz > SMALL) flag = 1;
@ -676,21 +681,11 @@ double FixQEqReax::calculate_H( double r, double gamma)
int FixQEqReax::CG( double *b, double *x)
{
int i, j, imax;
int i, j;
double tmp, alpha, beta, b_norm;
double sig_old, sig_new;
int nn, jj;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
ilist = list->ilist;
}
imax = 200;
int jj;
pack_flag = 1;
sparse_matvec( &H, x, q);
@ -748,18 +743,7 @@ int FixQEqReax::CG( double *b, double *x)
void FixQEqReax::sparse_matvec( sparse_matrix *A, double *x, double *b)
{
int i, j, itr_j;
int nn, NN, ii;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
NN = reaxc->list->inum + reaxc->list->gnum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
}
int ii;
for (ii = 0; ii < nn; ++ii) {
i = ilist[ii];
@ -794,16 +778,7 @@ void FixQEqReax::calculate_Q()
double u, s_sum, t_sum;
double *q = atom->q;
int nn, ii;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
ilist = list->ilist;
}
int ii;
s_sum = parallel_vector_acc( s, nn);
t_sum = parallel_vector_acc( t, nn);
@ -988,12 +963,6 @@ double FixQEqReax::parallel_norm( double *v, int n)
double my_sum, norm_sqr;
int ii;
int *ilist;
if (reaxc)
ilist = reaxc->list->ilist;
else
ilist = list->ilist;
my_sum = 0.0;
norm_sqr = 0.0;
@ -1016,12 +985,6 @@ double FixQEqReax::parallel_dot( double *v1, double *v2, int n)
double my_dot, res;
int ii;
int *ilist;
if (reaxc)
ilist = reaxc->list->ilist;
else
ilist = list->ilist;
my_dot = 0.0;
res = 0.0;
@ -1044,12 +1007,6 @@ double FixQEqReax::parallel_vector_acc( double *v, int n)
double my_acc, res;
int ii;
int *ilist;
if (reaxc)
ilist = reaxc->list->ilist;
else
ilist = list->ilist;
my_acc = 0.0;
res = 0.0;
@ -1070,12 +1027,6 @@ void FixQEqReax::vector_sum( double* dest, double c, double* v,
double d, double* y, int k)
{
int kk;
int *ilist;
if (reaxc)
ilist = reaxc->list->ilist;
else
ilist = list->ilist;
for (--k; k>=0; --k) {
kk = ilist[k];
@ -1089,12 +1040,6 @@ void FixQEqReax::vector_sum( double* dest, double c, double* v,
void FixQEqReax::vector_add( double* dest, double c, double* v, int k)
{
int kk;
int *ilist;
if (reaxc)
ilist = reaxc->list->ilist;
else
ilist = list->ilist;
for (--k; k>=0; --k) {
kk = ilist[k];
@ -1102,3 +1047,4 @@ void FixQEqReax::vector_add( double* dest, double c, double* v, int k)
dest[kk] += c * v[kk];
}
}

4
src/USER-REAXC/fix_qeq_reax.h
View File

@ -57,12 +57,13 @@ class FixQEqReax : public Fix {
protected:
int nevery,reaxflag;
int n, N, m_fill;
int nn, NN, m_fill;
int n_cap, nmax, m_cap;
int pack_flag;
int nlevels_respa;
class NeighList *list;
class PairReaxC *reaxc;
int *ilist, *jlist, *numneigh, **firstneigh;
double swa, swb; // lower/upper Taper cutoff radius
double Tap[8]; // Taper function
@ -94,6 +95,7 @@ class FixQEqReax : public Fix {
//CG storage
double *p, *q, *r, *d;
int imax;
//GMRES storage
//double *g,*y;