git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7454 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/USER-REAXC/Install.sh

@@ -5,10 +5,14 @@ if (test $1 = 1) then
     cp -p pair_reax_c.cpp ..
     cp -p fix_qeq_reax.cpp ..
     cp -p fix_reax_c.cpp ..
+    cp -p fix_reaxc_bonds.cpp ..
+    cp -p fix_reaxc_species.cpp ..
 
     cp -p pair_reax_c.h ..
     cp -p fix_qeq_reax.h ..
     cp -p fix_reax_c.h ..
+    cp -p fix_reaxc_bonds.h ..
+    cp -p fix_reaxc_species.h ..
 
     cp -p reaxc_allocate.cpp ..
     cp -p reaxc_basic_comm.cpp ..
@@ -61,10 +65,14 @@ elif (test $1 = 0) then
     rm -f ../pair_reax_c.cpp
     rm -f ../fix_qeq_reax.cpp
     rm -f ../fix_reax_c.cpp
+    rm -f ../fix_reaxc_bonds.cpp 
+    rm -f ../fix_reaxc_species.cpp 
 
     rm -f ../pair_reax_c.h
     rm -f ../fix_qeq_reax.h
     rm -f ../fix_reax_c.h
+    rm -f ../fix_reaxc_bonds.h
+    rm -f ../fix_reaxc_species.h
 
     rm -f ../reaxc_allocate.cpp
     rm -f ../reaxc_basic_comm.cpp

src/USER-REAXC/pair_reax_c.cpp

@@ -47,6 +47,8 @@
 #include "reaxc_reset_tools.h"
 #include "reaxc_traj.h"
 #include "reaxc_vector.h"
+#include "fix_reaxc_bonds.h"
+#include "fix_reaxc_species.h"
 
 using namespace LAMMPS_NS;
 
@@ -167,7 +169,7 @@ void PairReaxC::allocate( )
 
 void PairReaxC::settings(int narg, char **arg)
 {
-  if (narg != 1 && narg != 3) error->all(FLERR,"Illegal pair_style command");
+  if (narg < 1) error->all(FLERR,"Illegal pair_style command");
 
   // read name of control file or use default controls
 
@@ -194,6 +196,7 @@ void PairReaxC::settings(int narg, char **arg)
   // default values
 
   qeqflag = 1;
+  control->lgflag = 0;
 
   // process optional keywords
 
@@ -206,6 +209,12 @@ void PairReaxC::settings(int narg, char **arg)
       else if (strcmp(arg[iarg+1],"no") == 0) qeqflag = 0;
       else error->all(FLERR,"Illegal pair_style reax/c command");
       iarg += 2;
+    } else if (strcmp(arg[iarg],"lgvdw") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal pair_style reax/c command");
+      if (strcmp(arg[iarg+1],"yes") == 0) control->lgflag = 1;
+      else if (strcmp(arg[iarg+1],"no") == 0) control->lgflag = 0;
+      else error->all(FLERR,"Illegal pair_style reax/c command");
+      iarg += 2;
     } else error->all(FLERR,"Illegal pair_style reax/c command");
   }
 
@@ -400,6 +409,13 @@ void PairReaxC::compute(int eflag, int vflag)
   if (eflag || vflag) ev_setup(eflag,vflag);
   else evflag = vflag_fdotr = eflag_global = vflag_global = 0;
 
+/*  if ((eflag_atom || vflag_atom) && firstwarn) {
+    firstwarn = 0;
+    if (comm->me == 0) 
+      error->warning(FLERR,"Pair reax/c cannot yet compute "
+		     "per-atom energy or stress");
+  } */
+
   if (vflag_global) control->virial = 1;
   else control->virial = 0;
 
@@ -411,10 +427,6 @@ void PairReaxC::compute(int eflag, int vflag)
   system->big_box.box_norms[0] = 0; 
   system->big_box.box_norms[1] = 0; 
   system->big_box.box_norms[2] = 0;
-
-  system->evflag = evflag;
-  system->vflag_atom = vflag_atom;
-
   if( comm->me == 0 ) t_start = MPI_Wtime();
 
   // setup data structures
@@ -494,6 +506,12 @@ void PairReaxC::compute(int eflag, int vflag)
 
   Output_Results( system, control, data, &lists, out_control, mpi_data );
 
+  if(fixbond_flag) 
+	  fixbond( system, control, data, &lists, out_control, mpi_data );
+
+  if(fixspecies_flag) 
+	  fixspecies( system, control, data, &lists, out_control, mpi_data );
+
 }
 
 /* ---------------------------------------------------------------------- */
@@ -745,7 +763,7 @@ void PairReaxC::read_reax_forces()
 
 /* ---------------------------------------------------------------------- */
 
-void *PairReaxC::extract(const char *str, int &dim)
+void *PairReaxC::extract(char *str, int &dim)
 {
   dim = 1;
   if (strcmp(str,"chi") == 0 && chi) {
@@ -768,3 +786,5 @@ void *PairReaxC::extract(const char *str, int &dim)
   }
   return NULL;
 }
+
+/* ---------------------------------------------------------------------- */

src/USER-REAXC/pair_reax_c.h

@@ -45,16 +45,17 @@ class PairReaxC : public Pair {
   void init_style();
   double init_one(int, int);
   void *extract(char *, int &);
+  int fixbond_flag, fixspecies_flag;
 
- private:
-  reax_system *system;
   control_params *control;
+  reax_system *system;
+  output_controls *out_control;
   simulation_data *data;
   storage *workspace;
   reax_list *lists;
-  output_controls *out_control;
   mpi_datatypes *mpi_data;
 
+ private:
   double cutmax;
   int *map;
   class FixReaxC *fix_reax;

src/USER-REAXC/reaxc_allocate.cpp

@@ -92,6 +92,7 @@ inline void reax_atom_Copy( reax_atom *dest, reax_atom *src )
   dest->Hindex = src->Hindex;
   dest->num_bonds = src->num_bonds;
   dest->num_hbonds = src->num_hbonds;
+  dest->numbonds = src->numbonds;
 }
 
 

src/USER-REAXC/reaxc_bond_orders.cpp

@@ -596,7 +596,7 @@ void Add_dBond_to_Forces( reax_system *system, int i, int pj,
     rvec_ScaledAdd( workspace->f[k], -coef.C3dbopi2, nbr_k->bo_data.dBOp );
 
     // tally into per-atom virial
-    if( system->vflag_atom ) {
+    if( system->pair_ptr->vflag_atom ) {
       f_scaler = -(coef.C2dbo + coef.C2dDelta + coef.C3dbopi + coef.C3dbopi2);
       rvec_Scale( fk_tmp, -f_scaler, nbr_k->bo_data.dBOp);
       system->pair_ptr->v_tally(k, fk_tmp);
@@ -642,7 +642,7 @@ void Add_dBond_to_Forces( reax_system *system, int i, int pj,
     rvec_ScaledAdd( workspace->f[k], -coef.C4dbopi2, nbr_k->bo_data.dBOp );
 
     // tally into per-atom virial
-    if( system->vflag_atom ) {
+    if( system->pair_ptr->vflag_atom ) {
       f_scaler = -(coef.C3dbo + coef.C3dDelta + coef.C4dbopi + coef.C4dbopi2);
       rvec_Scale( fk_tmp, -f_scaler, nbr_k->bo_data.dBOp);
       system->pair_ptr->v_tally(k, fk_tmp);
@@ -673,7 +673,7 @@ void Add_dBond_to_Forces( reax_system *system, int i, int pj,
   /*3rd, dBOpi2*/
   rvec_ScaledAdd( workspace->f[j], coef.C4dbopi2, workspace->dDeltap_self[j] );
 
-  if( system->vflag_atom) {
+  if( system->pair_ptr->vflag_atom) {
 
     // forces on i 
     f_scaler = coef.C1dbo + coef.C1dDelta + coef.C2dbopi + coef.C2dbopi2;

src/USER-REAXC/reaxc_bonds.cpp

@@ -93,7 +93,7 @@ void Bonds( reax_system *system, control_params *control,
 	  -twbp->De_pp * bo_ij->BO_pi2;
 	
 	/* tally into per-atom energy */
-        if( system->evflag)
+        if( system->pair_ptr->evflag)
           system->pair_ptr->ev_tally(i,j,natoms,1,ebond,0.0,0.0,0.0,0.0,0.0);
 	
 	/* calculate derivatives of Bond Orders */
@@ -136,7 +136,7 @@ void Bonds( reax_system *system, control_params *control,
 	      (gp3*exphub1 + 25.0*gp4*exphuov*hulpov*(exphua1+exphub1));
 	    
 	    /* tally into per-atom energy */
-	    if( system->evflag)
+	    if( system->pair_ptr->evflag)
 	      system->pair_ptr->ev_tally(i,j,natoms,1,estriph,0.0,0.0,0.0,0.0,0.0);
 	    
 	    bo_ij->Cdbo += decobdbo;

src/USER-REAXC/reaxc_defs.h

@@ -104,6 +104,7 @@
 #define MAX_ITR 10
 #define RESTART 30
 
+#define MAX_BOND 20
 
 
 /******************* ENUMERATIONS *************************/

src/USER-REAXC/reaxc_ffield.cpp

@@ -25,6 +25,7 @@
   ----------------------------------------------------------------------*/
 
 #include "pair_reax_c.h"
+#include "error.h"
 #if defined(PURE_REAX)
 #include "ffield.h"
 #include "tool_box.h"
@@ -42,6 +43,8 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
   char   **tmp;
   char ****tor_flag;
   int      c, i, j, k, l, m, n, o, p, cnt;
+  int lgflag = control->lgflag;
+  int errorflag = 1;
   real     val;
   MPI_Comm comm;
 
@@ -161,9 +164,10 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
   reax->gp.vdw_type = 0;
   
   /* reading single atom parameters */
-  /* there are 4 lines of each single atom parameters in ff files. these
+  /* there are 4 or 5 lines of each single atom parameters in ff files, 
-     parameters later determine some of the pair and triplet parameters using 
+     depending on using lgvdw or not. These parameters later determine 
-     combination rules. */
+     some of the pair and triplet parameters using combination rules. */
+
   for( i = 0; i < reax->num_atom_types; i++ ) {
     /* line one */
     fgets( s, MAX_LINE, fp );
@@ -212,6 +216,12 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
     fgets( s, MAX_LINE, fp );
     c = Tokenize( s, &tmp );
 
+    /* Sanity check */
+    if (c < 3) {
+      fprintf(stderr, "Inconsistent ffield file (reaxc_ffield.cpp) \n");
+      MPI_Abort( comm, FILE_NOT_FOUND );
+    }
+
     val = atof(tmp[0]); reax->sbp[i].p_ovun2    = val;
     val = atof(tmp[1]); reax->sbp[i].p_val3     = val;
     val = atof(tmp[2]); 
@@ -221,17 +231,34 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
     val = atof(tmp[6]); reax->sbp[i].ecore2     = val;
     val = atof(tmp[7]); reax->sbp[i].acore2     = val;
 
+    /* line 5, only if lgvdw is yes */
+    if (lgflag) {
+      fgets( s, MAX_LINE, fp );
+      c = Tokenize( s, &tmp );
+
+      /* Sanity check */
+      if (c > 3) {
+        fprintf(stderr, "Inconsistent ffield file (reaxc_ffield.cpp) \n");
+        MPI_Abort( comm, FILE_NOT_FOUND );
+      }
+
+      val = atof(tmp[0]); reax->sbp[i].lgcij   	= val;
+      val = atof(tmp[1]); reax->sbp[i].lgre   	= val;
+    }
 
     if( reax->sbp[i].rcore2>0.01 && reax->sbp[i].acore2>0.01 ){ // Inner-wall
       if( reax->sbp[i].gamma_w>0.5 ){ // Shielding vdWaals
-	if( reax->gp.vdw_type != 0 && reax->gp.vdw_type != 3 )
+	if( reax->gp.vdw_type != 0 && reax->gp.vdw_type != 3 ) {
-	  fprintf( stderr, "Warning: inconsistent vdWaals-parameters\n"	\
+	  if (errorflag) 
+	    fprintf( stderr, "Warning: inconsistent vdWaals-parameters\n"	\
 		   "Force field parameters for element %s\n"		\
 		   "indicate inner wall+shielding, but earlier\n"	\
 		   "atoms indicate different vdWaals-method.\n"		\
 		   "This may cause division-by-zero errors.\n"		\
 		   "Keeping vdWaals-setting for earlier atoms.\n", 
 		   reax->sbp[i].name );
+	  errorflag = 0;
+	}
 	else{
 	  reax->gp.vdw_type = 3;
 #if defined(DEBUG)
@@ -439,6 +466,7 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
 	    reax->sbp[i].gamma,-1.5);
 
       // additions for additional vdWaals interaction types - inner core
+      
       reax->tbp[i][j].rcore = reax->tbp[j][i].rcore =
 	sqrt( reax->sbp[i].rcore2 * reax->sbp[j].rcore2 );
 
@@ -447,6 +475,15 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
 
       reax->tbp[i][j].acore = reax->tbp[j][i].acore =
 	sqrt( reax->sbp[i].acore2 * reax->sbp[j].acore2 );
+
+      // additions for additional vdWalls interaction types lg correction 
+
+      reax->tbp[i][j].lgcij = reax->tbp[j][i].lgcij =
+	sqrt( reax->sbp[i].lgcij * reax->sbp[j].lgcij );
+
+      reax->tbp[i][j].lgre = reax->tbp[j][i].lgre = 2.0 *	
+        sqrt( reax->sbp[i].lgre*reax->sbp[j].lgre );
+ 
     }
 
 
@@ -500,6 +537,12 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
 	reax->tbp[j][k].r_pp = val;
 	reax->tbp[k][j].r_pp = val;
       }
+
+      val = atof(tmp[8]); 
+      if (val >= 0.0) {
+	reax->tbp[j][k].lgcij = val;
+	reax->tbp[k][j].lgcij = val;
+      }
     }
   }
   

src/USER-REAXC/reaxc_hydrogen_bonds.cpp

@@ -146,8 +146,6 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
 	    CEhb3 = -hbp->p_hb3 * 
 	      (-hbp->r0_hb / SQR(r_jk) + 1.0 / hbp->r0_hb) * e_hb;
 
-	    fprintf(stderr,"	H bond called \n");
-		    
 	    /*fprintf( stdout, 
 	      "%6d%6d%6d%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n",
 	      system->my_atoms[i].orig_id, system->my_atoms[j].orig_id, 
@@ -192,7 +190,7 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
 	    }
 
 	    /* tally into per-atom virials */
-	    if (system->vflag_atom || system->evflag) {
+	    if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) {
 	      rvec_ScaledSum( delij, 1., system->my_atoms[i].x,
 				    -1., system->my_atoms[j].x );
 	      rvec_ScaledSum( delkj, 1., system->my_atoms[k].x,

src/USER-REAXC/reaxc_io_tools.cpp

@@ -25,6 +25,7 @@
   ----------------------------------------------------------------------*/
 
 #include "pair_reax_c.h"
+#include "update.h"
 #if defined(PURE_REAX)
 #include "io_tools.h"
 #include "basic_comm.h"
@@ -1030,6 +1031,70 @@ void Print_Total_Force( reax_system *system, simulation_data *data,
 	     workspace->f[i][0], workspace->f[i][1], workspace->f[i][2] );
 }
 
+void fixbond( reax_system *system, control_params *control, 
+		     simulation_data *data, reax_list **lists, 
+		     output_controls *out_control, mpi_datatypes *mpi_data )
+{
+  // count the number of bonds around each atom, for fix reax/c/bond
+  int i, j, pj, my_bonds_0, i_id, j_id;
+  int my_bonds_max = 0;
+  double BO_tmp;
+
+  control->bg_cut = 0.3;  // this values will not change with control file
+  reax_list *bonds = (*lists) + BONDS;
+
+  for( i=0; i < system->n; ++i ) {
+    my_bonds_0 = 0;
+    i_id = system->my_atoms[i].orig_id;  // orig_id is atom->tag
+    for( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj ) {
+      j = bonds->select.bond_list[pj].nbr;
+      j_id = system->my_atoms[j].orig_id;
+      BO_tmp = bonds->select.bond_list[pj].bo_data.BO;
+      if( i_id != j_id && BO_tmp >= control->bg_cut ) {
+	++my_bonds_0;
+	system->my_atoms[i].nbr_id[my_bonds_0] = j_id;
+	system->my_atoms[i].nbr_bo[my_bonds_0] = BO_tmp;
+      }
+    }
+    my_bonds_max = MAX(my_bonds_0, my_bonds_max);
+    system->my_atoms[i].numbonds = my_bonds_0;
+    system->my_bonds = my_bonds_max;
+  }
+}
+
+
+void fixspecies( reax_system *system, control_params *control, 
+		     simulation_data *data, reax_list **lists, 
+		     output_controls *out_control, mpi_datatypes *mpi_data )
+{
+  // count the number of bonds around each atom, for fix reax/c/bond
+  int i, j, pj, my_bonds_0, i_id, j_id;
+  int my_bonds_max = 0;
+  double BO_tmp;
+
+  control->bg_cut = 0.3;  // this values will not change with control file
+  reax_list *bonds = (*lists) + BONDS;
+
+  for( i=0; i < system->n; ++i ) {
+    my_bonds_0 = 0;
+    i_id = system->my_atoms[i].orig_id;
+    for( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj ) {
+      j = bonds->select.bond_list[pj].nbr;
+      j_id = system->my_atoms[j].orig_id;
+      BO_tmp = bonds->select.bond_list[pj].bo_data.BO;
+      if( i_id != j_id && BO_tmp >= control->bg_cut ) {
+	++my_bonds_0;
+	system->my_atoms[i].nbr_id[my_bonds_0] = j_id;
+	system->my_atoms[i].nbr_bo[my_bonds_0] = BO_tmp;
+      }
+    }
+    my_bonds_max = MAX(my_bonds_0, my_bonds_max);
+    system->my_atoms[i].numbonds = my_bonds_0;
+    system->my_bonds = my_bonds_max;
+  }
+}
+
+
 void Output_Results( reax_system *system, control_params *control, 
 		     simulation_data *data, reax_list **lists, 
 		     output_controls *out_control, mpi_datatypes *mpi_data )

src/USER-REAXC/reaxc_io_tools.h

@@ -56,6 +56,10 @@ void  Print_Bond_List2( reax_system*, reax_list*, char* );
 void  Print_Total_Force( reax_system*, simulation_data*, storage* );
 void  Output_Results( reax_system*, control_params*, simulation_data*, 
 		      reax_list**, output_controls*, mpi_datatypes* );
+void  fixbond( reax_system*, control_params*, simulation_data*, 
+		      reax_list**, output_controls*, mpi_datatypes* );
+void  fixspecies( reax_system*, control_params*, simulation_data*, 
+		      reax_list**, output_controls*, mpi_datatypes* );
 
 #if defined(DEBUG_FOCUS) || defined(TEST_FORCES) || defined(TEST_ENERGY) 
 void Debug_Marker_Bonded( output_controls*, int );

src/USER-REAXC/reaxc_lookup.cpp

@@ -266,7 +266,7 @@ int Init_Lookup_Tables( reax_system *system, control_params *control,
 		     "lookup:LR[i,j].CEclmb", comm );
 	  
 	  for( r = 1; r <= control->tabulate; ++r ) {
-	    LR_vdW_Coulomb( system, workspace, i, j, r * dr, &(LR[i][j].y[r]) );
+	    LR_vdW_Coulomb( system, workspace, control, i, j, r * dr, &(LR[i][j].y[r]) );
 	    h[r] = LR[i][j].dx;
 	    fh[r] = LR[i][j].y[r].H;
 	    fvdw[r] = LR[i][j].y[r].e_vdW;

src/USER-REAXC/reaxc_multi_body.cpp

@@ -91,7 +91,7 @@ void Atom_Energy( reax_system *system, control_params *control,
     workspace->CdDelta[i] += CElp;  // lp - 1st term  
 
     /* tally into per-atom energy */
-    if( system->evflag)
+    if( system->pair_ptr->evflag)
       system->pair_ptr->ev_tally(i,i,system->n,1,e_lp,0.0,0.0,0.0,0.0,0.0);
 
 #ifdef TEST_ENERGY
@@ -128,7 +128,7 @@ void Atom_Energy( reax_system *system, control_params *control,
 	    workspace->CdDelta[i] += deahu2dsbo;
 
 	    /* tally into per-atom energy */
-	    if( system->evflag)
+	    if( system->pair_ptr->evflag)
               system->pair_ptr->ev_tally(i,j,system->n,1,e_lph,0.0,0.0,0.0,0.0,0.0);
 		
 #ifdef TEST_ENERGY
@@ -219,7 +219,7 @@ void Atom_Energy( reax_system *system, control_params *control,
       p_ovun4 * exp_ovun1 * SQR(inv_exp_ovun1) + CEunder2;
     
     /* tally into per-atom energy */
-    if( system->evflag) {
+    if( system->pair_ptr->evflag) {
       eng_tmp = e_ov + e_un;
       f_tmp = CEover3 + CEunder3;
       system->pair_ptr->ev_tally(i,i,system->n,1,eng_tmp,0.0,0.0,0.0,0.0,0.0);

src/USER-REAXC/reaxc_nonbonded.cpp

@@ -50,6 +50,7 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
   real Tap, dTap, dfn13, CEvd, CEclmb, de_core;
   real dr3gamij_1, dr3gamij_3;
   real e_ele, e_vdW, e_core, SMALL = 0.0001;
+  real e_lg, de_lg, r_ij5, r_ij6, re6;
   rvec temp, ext_press;
   two_body_parameters *twbp;
   far_neighbor_data *nbr_pj;
@@ -65,6 +66,7 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
   p_vdW1i = 1.0 / p_vdW1;
   e_core = 0;
   e_vdW = 0;
+  e_lg = de_lg = 0.0;
 
   for( i = 0; i < natoms; ++i ) {
     start_i = Start_Index(i, far_nbrs);
@@ -155,6 +157,19 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
 
 	  de_core = -(twbp->acore/twbp->rcore) * e_core;
 	  CEvd += dTap * e_core + Tap * de_core / r_ij;
+
+          //  lg correction, only if lgvdw is yes
+	  if (control->lgflag) {
+            r_ij5 = pow( r_ij, 5.0 ); 
+            r_ij6 = pow( r_ij, 6.0 );
+            re6 = pow( twbp->lgre, 6.0 );
+            e_lg = -(twbp->lgcij/( r_ij6 + re6 )); 
+            data->my_en.e_vdW += Tap * e_lg;
+
+            de_lg = -6.0 * e_lg *  r_ij5 / ( r_ij6 + re6 ) ;
+            CEvd += dTap * e_lg + Tap * de_lg / r_ij;
+	  }
+
 	}
 
       /*Coulomb Calculations*/
@@ -169,8 +184,8 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
 	( dTap -  Tap * r_ij / dr3gamij_1 ) / dr3gamij_3;
 
       /* tally into per-atom energy */
-      if( system->evflag || system->vflag_atom) {
+      if( system->pair_ptr->evflag || system->pair_ptr->vflag_atom) {
-        pe_vdw = Tap * (e_vdW + e_core);
+        pe_vdw = Tap * (e_vdW + e_core + e_lg);
 	rvec_ScaledSum( delij, 1., system->my_atoms[i].x,
 			      -1., system->my_atoms[j].x );
 	f_tmp = -(CEvd + CEclmb);
@@ -334,7 +349,7 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control,
 
 
       /* tally into per-atom energy */
-      if( system->evflag || system->vflag_atom) {
+      if( system->pair_ptr->evflag || system->pair_ptr->vflag_atom) {
 	rvec_ScaledSum( delij, 1., system->my_atoms[i].x,
 			      -1., system->my_atoms[j].x );
 	f_tmp = -(CEvd + CEclmb);
@@ -399,14 +414,13 @@ void Compute_Polarization_Energy( reax_system *system, simulation_data *data )
     data->my_en.e_pol += en_tmp;
 
     /* tally into per-atom energy */
-    if( system->evflag) 
+    if( system->pair_ptr->evflag) 
       system->pair_ptr->ev_tally(i,i,system->n,1,0.0,en_tmp,0.0,0.0,0.0,0.0);
   }
 }
 
-
 void LR_vdW_Coulomb( reax_system *system, storage *workspace, 
-		     int i, int j, real r_ij, LR_data *lr )
+	control_params *control, int i, int j, real r_ij, LR_data *lr )
 {
   real p_vdW1 = system->reax_param.gp.l[28];
   real p_vdW1i = 1.0 / p_vdW1;
@@ -415,11 +429,13 @@ void LR_vdW_Coulomb( reax_system *system, storage *workspace,
   real Tap, dTap, dfn13;
   real dr3gamij_1, dr3gamij_3;
   real e_core, de_core;
+  real e_lg, de_lg, r_ij5, r_ij6, re6;
   two_body_parameters *twbp;
 
   twbp = &(system->reax_param.tbp[i][j]);
   e_core = 0;
   de_core = 0;
+  e_lg = de_lg = 0.0;
 
   /* calculate taper and its derivative */
   Tap = workspace->Tap[7] * r_ij + workspace->Tap[6];
@@ -470,6 +486,19 @@ void LR_vdW_Coulomb( reax_system *system, storage *workspace,
       
       de_core = -(twbp->acore/twbp->rcore) * e_core;
       lr->CEvd += dTap * e_core + Tap * de_core / r_ij;
+
+      //  lg correction, only if lgvdw is yes
+      if (control->lgflag) {
+        r_ij5 = pow( r_ij, 5.0 ); 
+        r_ij6 = pow( r_ij, 6.0 );
+        re6 = pow( twbp->lgre, 6.0 );
+        e_lg = -(twbp->lgcij/( r_ij6 + re6 )); 
+        lr->e_vdW += Tap * e_lg;
+
+        de_lg = -6.0 * e_lg *  r_ij5 / ( r_ij6 + re6 ) ;
+        lr->CEvd += dTap * e_lg + Tap * de_lg/r_ij;
+      }
+
     }
 
   

src/USER-REAXC/reaxc_nonbonded.h

@@ -38,5 +38,6 @@ void Tabulated_vdW_Coulomb_Energy( reax_system*, control_params*,
 
 void Compute_Polarization_Energy( reax_system*, simulation_data* );
 
-void LR_vdW_Coulomb( reax_system*, storage*, int, int, real, LR_data* );
+void LR_vdW_Coulomb( reax_system*, storage*, control_params*,
+		int, int, real, LR_data* );
 #endif

src/USER-REAXC/reaxc_torsion_angles.cpp

@@ -484,7 +484,7 @@ void Torsion_Angles( reax_system *system, control_params *control,
 		  }
 
 		  /* tally into per-atom virials */
-		  if( system->vflag_atom || system->evflag) {
+		  if( system->pair_ptr->vflag_atom || system->pair_ptr->evflag) {
 
 		    // acquire vectors
 		    rvec_ScaledSum( delil, 1., system->my_atoms[l].x,
@@ -509,9 +509,9 @@ void Torsion_Angles( reax_system *system, control_params *control,
 
 		    // tally
 		    eng_tmp = e_tor + e_con;
-		    if( system->evflag)
+		    if( system->pair_ptr->evflag)
 			    system->pair_ptr->ev_tally(j,k,natoms,1,eng_tmp,0.0,0.0,0.0,0.0,0.0);
-		    if( system->vflag_atom)
+		    if( system->pair_ptr->vflag_atom)
 			    system->pair_ptr->v_tally4(i,j,k,l,fi_tmp,fj_tmp,fk_tmp,delil,deljl,delkl);
 		  }
 		  

src/USER-REAXC/reaxc_traj.cpp

@@ -815,6 +815,7 @@ int Write_Bonds(reax_system *system, control_params *control, reax_list *bonds,
 	  bonds->select.bond_list[pj].bo_data.BO >= control->bg_cut )
 	++my_bonds;
     }
+
   /* allreduce - total number of bonds */
   MPI_Allreduce( &my_bonds, &num_bonds, 1, MPI_INT, MPI_SUM, mpi_data->world );
 
@@ -828,10 +829,11 @@ int Write_Bonds(reax_system *system, control_params *control, reax_list *bonds,
     Reallocate_Output_Buffer( out_control, buffer_req, mpi_data->world );
   
   /* fill in the buffer */
-  my_bonds = 0;
   out_control->line[0] = 0;
   out_control->buffer[0] = 0;
-  for( i=0; i < system->n; ++i )
+
+  my_bonds = 0;
+  for( i=0; i < system->n; ++i ) {
     for( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj ) {
       bo_ij = &( bonds->select.bond_list[pj] );
       j = bo_ij->nbr;
@@ -854,12 +856,13 @@ int Write_Bonds(reax_system *system, control_params *control, reax_list *bonds,
 	  fprintf(stderr, "write_traj_bonds: FATAL! invalid bond_info option");
 	  MPI_Abort( mpi_data->world, UNKNOWN_OPTION );
 	}
-	    
 	strncpy( out_control->buffer + my_bonds*line_len, 
 		 out_control->line, line_len+1 );
 	++my_bonds;
       }
     }
+  }
+
 
 #if defined(PURE_REAX)
   if( out_control->traj_method == MPI_TRAJ ) {

src/USER-REAXC/reaxc_types.h

@@ -57,6 +57,7 @@
 #define REAX_MAX_4BODY_PARAM    5
 #define REAX_MAX_ATOM_TYPES     25
 #define REAX_MAX_MOLECULE_SIZE  20
+#define MAX_BOND    		20  // same as reaxc_defs.h
 
 /********************** TYPE DEFINITIONS ********************/
 typedef int  ivec[3];
@@ -246,6 +247,11 @@ typedef struct
   real rcore2;
   real ecore2;
   real acore2;
+
+  /* Line five in the ffield file, only for lgvdw yes */
+  real lgcij;
+  real lgre;
+
 } single_body_parameters;
 
 
@@ -270,6 +276,7 @@ typedef struct {
   real r_vdW;
   real gamma_w;
   real rcore, ecore, acore;
+  real lgcij, lgre;
 
   /* electrostatic parameters */
   real gamma; // note: this parameter is gamma^-3 and not gamma.
@@ -359,6 +366,11 @@ typedef struct
   int num_bonds;
   int num_hbonds;
   int renumber;
+
+  int numbonds;  		// true number of bonds around atoms
+  int nbr_id[MAX_BOND];  	// ids of neighbors around atoms
+  double nbr_bo[MAX_BOND];  	// BO values of bond between i and nbr
+  double sum_bo, no_lp;   	// sum of BO values and no. of lone pairs
 } reax_atom;
 
 
@@ -472,12 +484,11 @@ typedef struct
   simulation_box   big_box, my_box, my_ext_box;
   grid             my_grid;
   boundary_cutoff  bndry_cuts;
-
   reax_atom       *my_atoms;
 
-  int evflag;
-  int vflag_atom;
   class Pair *pair_ptr;
+  int my_bonds;
+
 } reax_system;
 
 
@@ -543,6 +554,9 @@ typedef struct
   int  diffusion_coef;
   int  freq_diffusion_coef;
   int  restrict_type;
+
+  int lgflag;
+
 } control_params;
 
 

src/USER-REAXC/reaxc_valence_angles.cpp

@@ -416,7 +416,7 @@ void Valence_Angles( reax_system *system, control_params *control,
 		}
 
 		/* tally into per-atom virials */
-		if( system->vflag_atom || system->evflag) {
+		if( system->pair_ptr->vflag_atom || system->pair_ptr->evflag) {
 			
 		  /* Acquire vectors */
 		  rvec_ScaledSum( delij, 1., system->my_atoms[i].x,
@@ -430,9 +430,9 @@ void Valence_Angles( reax_system *system, control_params *control,
 		  
 		  eng_tmp = e_ang + e_pen + e_coa;
 
-		  if( system->evflag)
+		  if( system->pair_ptr->evflag)
 			  system->pair_ptr->ev_tally(j,j,system->N,1,eng_tmp,0.0,0.0,0.0,0.0,0.0);
-		  if( system->vflag_atom)
+		  if( system->pair_ptr->vflag_atom)
 			  system->pair_ptr->v_tally3(i,j,k,fi_tmp,fk_tmp,delij,delkj);
 		}