Merge branch 'master' into lammps-icms

doc/fix_ave_histo.html

@@ -295,11 +295,11 @@ values can only be accessed on timesteps that are multiples of <I>Nfreq</I>
 since that is when a histogram is generated.
 The global vector has 4 values:
 </P>
-<P>1 = total counts in the histogram
-2 = values that were not histogrammed (see <I>beyond</I> keyword)
-3 = min value of all input values, including ones not histogrammed
-4 = max value of all input values, including ones not histogrammed
-</P>
+<UL><LI>1 = total counts in the histogram
+<LI>2 = values that were not histogrammed (see <I>beyond</I> keyword)
+<LI>3 = min value of all input values, including ones not histogrammed
+<LI>4 = max value of all input values, including ones not histogrammed 
+</UL>
 <P>The global array has # of rows = Nbins and # of columns = 3.  The
 first column has the bin coordinate, the 2nd column has the count of
 values in that histogram bin, and the 3rd column has the bin count

doc/fix_ave_histo.txt

@@ -283,7 +283,7 @@ The global vector has 4 values:
 1 = total counts in the histogram
 2 = values that were not histogrammed (see {beyond} keyword)
 3 = min value of all input values, including ones not histogrammed
-4 = max value of all input values, including ones not histogrammed
+4 = max value of all input values, including ones not histogrammed :ul
 
 The global array has # of rows = Nbins and # of columns = 3.  The
 first column has the bin coordinate, the 2nd column has the count of

examples/ELASTIC/displace.mod

@@ -4,17 +4,17 @@
 
 # Find which reference length to use
 
-if ${dir} == 1 then &
+if "${dir} == 1" then &
    "variable len0 equal ${lx0}" 
-if ${dir} == 2 then &
+if "${dir} == 2" then &
    "variable len0 equal ${ly0}" 
-if ${dir} == 3 then &
+if "${dir} == 3" then &
    "variable len0 equal ${lz0}" 
-if ${dir} == 4 then &
+if "${dir} == 4" then &
    "variable len0 equal ${lz0}" 
-if ${dir} == 5 then &
+if "${dir} == 5" then &
    "variable len0 equal ${lz0}" 
-if ${dir} == 6 then &
+if "${dir} == 6" then &
    "variable len0 equal ${ly0}" 
 
 # Reset box and simulation parameters
@@ -24,17 +24,17 @@ include potential.mod
 
 # Negative deformation
 variable delta equal -${up}*${len0}
-if ${dir} == 1 then &
+if "${dir} == 1" then &
    "displace_box all x delta 0 ${delta} units box"
-if ${dir} == 2 then &
+if "${dir} == 2" then &
    "displace_box all y delta 0 ${delta} units box"
-if ${dir} == 3 then &
+if "${dir} == 3" then &
    "displace_box all z delta 0 ${delta} units box"
-if ${dir} == 4 then &
+if "${dir} == 4" then &
    "displace_box all yz delta ${delta} units box"
-if ${dir} == 5 then &
+if "${dir} == 5" then &
    "displace_box all xz delta ${delta} units box"
-if ${dir} == 6 then &
+if "${dir} == 6" then &
    "displace_box all xy delta ${delta} units box"
 
 # Relax atoms positions
@@ -72,17 +72,17 @@ include potential.mod
 
 # Positive deformation
 variable delta equal ${up}*${len0}
-if ${dir} == 1 then &
+if "${dir} == 1" then &
    "displace_box all x delta 0 ${delta} units box"
-if ${dir} == 2 then &
+if "${dir} == 2" then &
    "displace_box all y delta 0 ${delta} units box"
-if ${dir} == 3 then &
+if "${dir} == 3" then &
    "displace_box all z delta 0 ${delta} units box"
-if ${dir} == 4 then &
+if "${dir} == 4" then &
    "displace_box all yz delta ${delta} units box"
-if ${dir} == 5 then &
+if "${dir} == 5" then &
    "displace_box all xz delta ${delta} units box"
-if ${dir} == 6 then &
+if "${dir} == 6" then &
    "displace_box all xy delta ${delta} units box"
 
 # Relax atoms positions

examples/reax/AB/README

@@ -0,0 +1,30 @@
+Disclaimer:  Using these force fields for systems they 
+have not been explicitly trained against may produce
+unrealistic results.  Please see the README file in 
+each subdirectory for more detailed information.
+
+Ammonia Borane:
+
+     The follow information is reproduced from:
+ 
+     "Weismiller, M.R.; van Duin, A.C.T.; Lee, J.; 
+     Yetter, R.A. J. Phys. Chem. A 2010, 114, 5485-5492"
+
+   - QM data were generated describing the single and 
+   (if relevant) double and triple bond dissociation 
+   for all B/N/O/H combinations. These data were used 
+   to derive initial ReaxFF bond parameters, and all 
+   calculations were performed using DFT with the B3LYP 
+   functional and the Pople 6-311G** basis set.
+
+   - The training set was then extended with QM data 
+   describing angular distortions in a set of small 
+   AB-related (AB = H3N-BH3) molecules. These data 
+   were used to derive the initial ReaxFF angular 
+   parameters.
+
+   - The training set was extended with reaction barriers 
+   for key reaction steps such as H2 release 
+   from AB, dimerization of H2B-NH2 and reaction 
+   energies associated with H2 release from AB and with AB 
+   oxidation.

examples/reax/Au_O/README

@@ -0,0 +1,18 @@
+Disclaimer:  Using these force fields for systems they 
+have not been explicitly trained against may produce
+unrealistic results.  Please see the README file in 
+each subdirectory for more detailed information.
+
+Au/O
+
+     The follow information is reproduced from
+     "Keith, J. A.;  Fantauzzi, D.; Jacob, T.; 
+     van Duin, A. C. T. Phys Rev B 2010, 81, 235404"
+
+   - The force field optimization involved parameterization 
+   of the fcc, bcc, ideal-hcp, sc, diamond, and a15 bulk phases of Au.  
+   For these phases the ReaxFF force field gives good 
+   agreement for the binding energy, volume of minimum 
+   energy, and curvature of the binding well around the 
+   minimum compared to the QM calculations used 
+   for parameterization.

examples/reax/CHO/README

@@ -0,0 +1,25 @@
+Disclaimer:  Using these force fields for systems they 
+have not been explicitly trained against may produce
+unrealistic results.  Please see the README file in 
+each subdirectory for more detailed information.
+
+Hydrocarbon oxidation C/H/O:
+
+     The follow information is reproduced from:
+
+     "Chenoweth, K.; van Duin, A.C.T.; Goddard, W.A. 
+     J. Phys. Chem. A 2008, 112, 1040-1053."
+
+     - To obtain the H/C/O compound data required to 
+     extend the hydrocarbon-training set, DFT 
+     calculations were performed on the following systems: 
+     (a) dissociation energies for various bonds 
+     containing carbon, oxygen, and hydrogen.  The 
+     ground state structure was obtained through 
+     full geometry optimization.  Dissociation curves 
+     were calculated by constraining only the bond length of 
+     interest and re-optimization of the remaining 
+     internal coordinates. Optimization was also performed 
+     for the various angles and torsions associated 
+     with C/H/O interactions.
+

examples/reax/Fe_O_H/README

@@ -0,0 +1,29 @@
+Disclaimer:  Using these force fields for systems they 
+have not been explicitly trained against may produce
+unrealistic results.  Please see the README file in 
+each subdirectory for more detailed information.
+
+Fe/O/H
+
+     The follow information is reproduced from:
+
+     "Aryanpour, M.; van Duin, A. C. T.; Kubicki, 
+     J. D. J. Phys. Chem. A 2010, 114, 6298-6307"
+
+     - The initial force field parameters for the 
+     Fe-Fe parameters were taken from an 
+     earlier force field development 
+     project on bulk-iron metal, based on 
+     DFT-calculations on antiferromagnetic 
+     BCC and FCC. The DFT data can 
+     be found in Ref 31 of the above-mentioned manuscript. 
+     The O/H parameters were taken from the ReaxFF bulk 
+     water description. The Fe/Fe and O/H parameters were 
+     kept fixed to these initial values,
+     whereas the Fe/O parameters were reoptimized 
+     against the quantum mechanical results 
+     presented in the above-mentioned manuscript.
+
+     - Detailed information on the force field parameters 
+     is given in the supporting information of the above mentioned
+     manuscript.

examples/reax/RDX/README

@@ -0,0 +1,40 @@
+Disclaimer:  Using these force fields for systems they 
+have not been explicitly trained against may produce
+unrealistic results.  Please see the README file in 
+each subdirectory for more detailed information.
+
+RDX/High Energy C/H/O/N:
+
+     The follow information is reproduced from:
+
+     "Strachan, A.; van Duin, A. C. T.; Chakraborty, D.; 
+     Dasgupta, S.; Goddard III, W. A. Phys. Rev. Lett. 
+     2003, 91, 098301"
+
+     "Strachan, A.; Kober, E.; van Duin, A. C. T.; 
+     Oxgaard, J.; Goddard III, W. A. J. Chem. Phys. 
+     2005, 122, 054502"
+
+     "Zhang, L.; van Duin, A.C.T.; Zybin, S.; 
+     Goddard, W.A. J. Phys. Chem. B 2009, 113, 10770-10778"
+
+     "Zhang, L.; Zybin, S.; van Duin, A.C.T.; Dasgupta, S.; 
+     Goddard, W.A.; Kober, E. J. Phys. Chem. A. 2009, 
+     113, 10619-10640"
+
+     - The parameters of the nitramine ReaxFF are based 
+     on a large number of ab initio QM calculations. 
+     Over 40 reactions and over 1600 equilibrated molecules 
+     have been used; they are designed to characterize 
+     the atomic interactions under various environments 
+     likely and unlikely high energy each atom can encounter. 
+     The training set contains bond breaking and compression 
+     curves for all possible bonds, angle and torsion bending 
+     data for all possible cases, as well as crystal 
+     data.
+
+     - Please see the supplemental material from 
+     Phys. Rev. Lett. 2003, 91, 098301 
+     for a detailed description of the 
+     parameterization of this force field.
+

examples/reax/README

@@ -9,17 +9,21 @@ contains validated parameter files for a particular
 published version of ReaxFF. In all cases, the examples
 use pair_style reax/c.
 
-CHO: Reactive MD-force field c/h/o combustion force field
-This version requires tripflag = 0 and hbnewflag = 1
+Disclaimer:  Using these force fields for systems they 
+have not been explicitly trained against may produce
+unrealistic results.  Please see the README file in 
+each subdirectory for more detailed information.
 
-AB: Reactive MD-force field for Ammonia Borane: Weismiller, M.R.; van Duin, A.C.T.; Lee, J.; Yetter, R.A. J. Phys. Chem. A 2010, 114, 5485-5492
+AB (Ammonia Borane): Weismiller, M.R.; van Duin, A.C.T.; Lee, J.; Yetter, R.A. J. Phys. Chem. A 2010, 114, 5485-5492
 
-Au_O: Reactive MD-force field: Au-Au: Keith, J. A. et al. Phys Rev B 2010, 81, 235404; Joshi, K. et al., J. Mat. Chem 2010, 20, 10431-10437
+Au_O (Gold-Oxygen): Keith, J. A. et al. Phys Rev B 2010, 81, 235404; Joshi, K. et al., J. Mat. Chem 2010, 20, 10431-10437
 
-Fe_O_C_H: Reactive MD-force field: Fe/O/water: Aryanpour, M.; van Duin, A. C. T.; Kubicki, J. D. J. Phys. Chem. A 2010, 114, 6298-6307
+CHO (C/H/O combustion): Chenoweth, K.; van Duin, A.C.T.; Goddard, W.A. J. Phys. Chem. A 2008, 112, 1040-1053.
 
-RDX: Reactive MD-force field nitramines (RDX/HMX/TATB/PETN): Strachan, A.; et al. Phys. Rev. Lett. 2003, 91, 098301
+Fe_O_C_H (Fe/O/water): Aryanpour, M.; van Duin, A. C. T.; Kubicki, J. D. J. Phys. Chem. A 2010, 114, 6298-6307
 
-V_O_H: Reactive MD-force field: V/O/C/H potential: Chenoweth, K.; et al. J. Phys. Chem. C, 2008, 112, 14645-14654. 
+RDX (Nitramines RDX/HMX/TATB/PETN): Strachan, A.; et al. Phys. Rev. Lett. 2003, 91, 098301
 
-Zn_O_H: Reactive MD-force field: water/zinc: Raymand, D.; van Duin A.C.T.; Baudin M.; Hermannson K. Surface Science 2008, 602, 1020-1031.
+V_O_H (Vanadium V/O/C/H potential): Chenoweth, K.; et al. J. Phys. Chem. C, 2008, 112, 14645-14654. 
+
+Zn_O_H (Water/Zinc): Raymand, D.; van Duin A.C.T.; Baudin M.; Hermannson K. Surface Science 2008, 602, 1020-1031.

examples/reax/V_O_H/README

@@ -0,0 +1,21 @@
+Disclaimer:  Using these force fields for systems they 
+have not been explicitly trained against may produce
+unrealistic results.  Please see the README file in 
+each subdirectory for more detailed information.
+
+V/O
+     The follow information is reproduced from:
+
+     "Chenoweth, K.; van Duin A.C.T.; Persson, P.; 
+     Cheng M.J.; Oxgaard, J.; Goddard W.A. 
+     J. Phys. Chem. C, 2008, 112, 14645-14654."
+
+     - The ReaxFF force field parameters have been fit 
+     to a large quantum mechanics (QM) training set containing over 700 
+     structures and energetics related to bond dissociations, angle 
+     and dihedral distortions, and reactions between hydrocarbons 
+     and vanadium oxide clusters. In addition, the training set 
+     contains charge distributions for small vanadium oxide clusters 
+     and the stabilities of condensed-phase systems 
+     including V2O5, VO2, and V2O3 in addition to metallic V (V0).
+

examples/reax/Zn_O_H/README

@@ -0,0 +1,28 @@
+Disclaimer:  Using these force fields for systems they 
+have not been explicitly trained against may produce
+unrealistic results.  Please see the README file in 
+each subdirectory for more detailed information.
+
+Zn/O/H:
+
+     The follow information is reproduced from:
+
+     "Raymand, D.; van Duin A.C.T.; Baudin M.; Hermannson K. 
+     Surface Science 2008, 602, 1020-1031."
+
+     "Raymand, D.; van Duin, A. C. T.; Spangberg, D.; 
+     Goddard, W. A.; Hermansson, K. Surface Science 
+     2010, 604, 9-10, 741-752."
+
+     - Based on QM calculations for Zn(s), ZnO(s), 
+     and Zn hydroxide clusters [Zn(OH)2 and O(ZnOH)2], 
+     ReaxFF parameters were generated for Zn-O and 
+     Zn-Zn bond energies and for Zn-O-Zn, O-Zn-O, 
+     O-Zn-Zn and Zn-O-H valence angle energies.
+
+     - QM calculations were performed for the four 
+     crystal polymorphs of the wurtzite, zincblende, 
+     rocksalt and caesium chloride structures 
+     (the structures are also referred to as h-ZnS, 
+     c-ZnS, NaCl and CsCl, respectively).
+

src/USER-REAXC/pair_reax_c.cpp

@@ -223,15 +223,12 @@ void PairReaxC::coeff( int nargs, char **args )
  
   // read args that map atom types to elements in potential file
   // map[i] = which element the Ith atom type is, -1 if NULL
-  // NOTE: for now throw an error if NULL is used to disallow use with hybrid
-  //       REAX/C lib will have to be modified to allow this
 
   int itmp;
   int nreax_types = system->reax_param.num_atom_types;
   for (int i = 3; i < nargs; i++) {
     if (strcmp(args[i],"NULL") == 0) {
       map[i-2] = -1;
-      error->all("Cannot currently use pair reax/c with pair hybrid");
       continue;
     }
 
@@ -242,6 +239,7 @@ void PairReaxC::coeff( int nargs, char **args )
 
     if (itmp < 0 || itmp >= nreax_types)
       error->all("Non-existent ReaxFF type");
+
   }
 
   int n = atom->ntypes;
@@ -292,7 +290,7 @@ void PairReaxC::init_style( )
 
   cutmax = MAX3(control->nonb_cut, control->hbond_cut, 2*control->bond_cut);
 
-  for(int i = 0; i < LIST_N; ++i )
+  for( int i = 0; i < LIST_N; ++i )
     lists[i].allocated = 0;
 
   if (fix_reax == NULL) {

src/USER-REAXC/reaxc_allocate.cpp

@@ -512,7 +512,7 @@ int Reallocate_Bonds_List( reax_system *system, reax_list *bonds,
   *total_bonds = 0;
   *est_3body = 0;
   for( i = 0; i < system->N; ++i ){
-    *est_3body += SQR( Num_Entries( i, bonds ) );
+    *est_3body += SQR(system->my_atoms[i].num_bonds);
     // commented out - already updated in validate_lists in forces.c
     // system->my_atoms[i].num_bonds = MAX( Num_Entries(i,bonds)*2, MIN_BONDS );
     *total_bonds += system->my_atoms[i].num_bonds;

src/USER-REAXC/reaxc_nonbonded.cpp

@@ -388,7 +388,7 @@ void LR_vdW_Coulomb( reax_system *system, storage *workspace,
 	
     lr->e_vdW = Tap * twbp->D * (exp1 - 2.0 * exp2);
     lr->CEvd = dTap * twbp->D * (exp1 - 2.0 * exp2) -
-      Tap * twbp->D * (twbp->alpha / twbp->r_vdW) * (exp1 - exp2);
+      Tap * twbp->D * (twbp->alpha / twbp->r_vdW) * (exp1 - exp2) / r_ij;
   }
 
   if(system->reax_param.gp.vdw_type==2 || system->reax_param.gp.vdw_type==3)
@@ -397,7 +397,7 @@ void LR_vdW_Coulomb( reax_system *system, storage *workspace,
       lr->e_vdW += Tap * e_core;
       
       de_core = -(twbp->acore/twbp->rcore) * e_core;
-      lr->CEvd += dTap * e_core + Tap * de_core;
+      lr->CEvd += dTap * e_core + Tap * de_core / r_ij;
     }
 
   

src/USER-REAXC/reaxc_valence_angles.cpp

@@ -235,7 +235,9 @@ void Valence_Angles( reax_system *system, control_params *control,
 	  
 	  
 	  if( (j < system->n) && (BOA_jk > 0.0) && 
-              (bo_ij->BO * bo_jk->BO > SQR(control->thb_cut)/*0*/) ) {
+              (bo_ij->BO > control->thb_cut) &&  
+              (bo_jk->BO > control->thb_cut) &&  
+              (bo_ij->BO * bo_jk->BO > 0.001) ) {
 	    r_jk = pbond_jk->d;		      
 	    thbh = &( system->reax_param.thbp[ type_i ][ type_j ][ type_k ] );
 	    

src/domain.cpp

@@ -1164,4 +1164,3 @@ void Domain::box_corners()
   corners[7][0] = 1.0; corners[7][1] = 1.0; corners[7][2] = 1.0;
   lamda2x(corners[7],corners[7]);
 }
-

src/fix_ave_histo.cpp

@@ -61,7 +61,7 @@ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) :
   size_vector = 4;
   extvector = 0;
   array_flag = 1;
-  size_local_cols = 3;
+  size_array_cols = 3;
   extarray = 0;
   time_depend = 1;
 
@@ -432,7 +432,7 @@ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) :
   // allocate and initialize memory for averaging
 
   if (beyond == EXTRA) nbins += 2;
-  size_local_rows = nbins;
+  size_array_rows = nbins;
 
   bin = new double[nbins];
   bin_total = new double[nbins];

src/neighbor.cpp

@@ -1101,6 +1101,7 @@ void Neighbor::build()
 	boxhi_hold[2] = bboxhi[2];
       } else {
 	domain->box_corners();
+	corners = domain->corners;
 	for (i = 0; i < 8; i++) {
 	  corners_hold[i][0] = corners[i][0];
 	  corners_hold[i][1] = corners[i][1];
@@ -1230,7 +1231,6 @@ void Neighbor::setup_bins()
     hi[1] = domain->subhi_lamda[1] + cutghost[1];
     hi[2] = domain->subhi_lamda[2] + cutghost[2];
     domain->bbox(lo,hi,bsubboxlo,bsubboxhi);
-    corners = domain->corners;
   }
 
   bbox[0] = bboxhi[0] - bboxlo[0];

src/version.h

@@ -1 +1 @@
-#define LAMMPS_VERSION "9 Feb 2011"
+#define LAMMPS_VERSION "11 Feb 2011"