diff --git a/lib/awpmd/ivutils/include/logexc.h b/lib/awpmd/ivutils/include/logexc.h
index 7bf6125bc6..2eebfcc0ce 100644
--- a/lib/awpmd/ivutils/include/logexc.h
+++ b/lib/awpmd/ivutils/include/logexc.h
@@ -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);
   }
diff --git a/lib/awpmd/ivutils/include/pairhash.h b/lib/awpmd/ivutils/include/pairhash.h
index c29f72773f..401b7b91e3 100644
--- a/lib/awpmd/ivutils/include/pairhash.h
+++ b/lib/awpmd/ivutils/include/pairhash.h
@@ -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){
diff --git a/lib/awpmd/ivutils/include/refobj.h b/lib/awpmd/ivutils/include/refobj.h
index a1a4feb0ff..8e5fc7df45 100644
--- a/lib/awpmd/ivutils/include/refobj.h
+++ b/lib/awpmd/ivutils/include/refobj.h
@@ -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 {
diff --git a/lib/awpmd/ivutils/include/vector_3.h b/lib/awpmd/ivutils/include/vector_3.h
index 468abbda99..5711d3bf35 100644
--- a/lib/awpmd/ivutils/include/vector_3.h
+++ b/lib/awpmd/ivutils/include/vector_3.h
@@ -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++);
diff --git a/lib/awpmd/ivutils/src/logexc.cpp b/lib/awpmd/ivutils/src/logexc.cpp
index 7e82af0f65..6e809dd65c 100644
--- a/lib/awpmd/ivutils/src/logexc.cpp
+++ b/lib/awpmd/ivutils/src/logexc.cpp
@@ -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,...){
diff --git a/lib/awpmd/systems/interact/TCP/wpmd.h b/lib/awpmd/systems/interact/TCP/wpmd.h
index e118b76f36..4c65440014 100644
--- a/lib/awpmd/systems/interact/TCP/wpmd.h
+++ b/lib/awpmd/systems/interact/TCP/wpmd.h
@@ -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.
diff --git a/lib/awpmd/systems/interact/TCP/wpmd_split.cpp b/lib/awpmd/systems/interact/TCP/wpmd_split.cpp
index 2c58bdaf28..f85b1ebf59 100644
--- a/lib/awpmd/systems/interact/TCP/wpmd_split.cpp
+++ b/lib/awpmd/systems/interact/TCP/wpmd_split.cpp
@@ -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)
diff --git a/lib/awpmd/systems/interact/TCP/wpmd_split.h b/lib/awpmd/systems/interact/TCP/wpmd_split.h
index 44fb757576..10a8bda558 100644
--- a/lib/awpmd/systems/interact/TCP/wpmd_split.h
+++ b/lib/awpmd/systems/interact/TCP/wpmd_split.h
@@ -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