Reax ff

.gitignore

@@ -17,4 +17,7 @@ log.lammps
 
 crates.d/config.toml
 
+log.cite
+etc/profile.d/lammps.csh
+
 venv

build.rs

@@ -1,17 +1,17 @@
-extern crate vergen;
-
-use vergen::{ConstantsFlags, Result, Vergen};
-
-fn main() {
-    gen_constants().expect("Unable to generate vergen constants!");
-}
-
-fn gen_constants() -> Result<()> {
-    let vergen = Vergen::new(ConstantsFlags::all())?;
-
-    for (k, v) in vergen.build_info() {
-        println!("cargo:rustc-env={}={}", k.name(), v);
-    }
-
-    Ok(())
-}
+extern crate vergen;
+
+use vergen::{ConstantsFlags, Result, Vergen};
+
+fn main() {
+    gen_constants().expect("Unable to generate vergen constants!");
+}
+
+fn gen_constants() -> Result<()> {
+    let vergen = Vergen::new(ConstantsFlags::all())?;
+
+    for (k, v) in vergen.build_info() {
+        println!("cargo:rustc-env={}={}", k.name(), v);
+    }
+
+    Ok(())
+}

etc/profile.d/lammps.csh

@@ -0,0 +1,4 @@
+# set environment for LAMMPS and msi2lmp executables
+# to find potential and force field files
+if ( "$?LAMMPS_POTENTIALS" == 0 ) setenv LAMMPS_POTENTIALS /home/brandon/.local/share/lammps/potentials
+if ( "$?MSI2LMP_LIBRARY" == 0 ) setenv MSI2LMP_LIBRARY /home/brandon/.local/share/lammps/frc_files

etc/profile.d/lammps.sh

@@ -0,0 +1,5 @@
+# set environment for LAMMPS and msi2lmp executables
+# to find potential and force field files
+LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS-/home/brandon/.local/share/lammps/potentials}
+MSI2LMP_LIBRARY=${MSI2LMP_LIBRARY-/home/brandon/.local/share/lammps/frc_files}
+export LAMMPS_POTENTIALS MSI2LMP_LIBRARY

log.cite

@@ -0,0 +1,26 @@
+This LAMMPS simulation made specific use of work described in the
+following references.  See http://lammps.sandia.gov/cite.html
+for details.
+
+pair reax/c command:
+
+@Article{Aktulga12,
+ author = {H. M. Aktulga, J. C. Fogarty, S. A. Pandit, A. Y. Grama},
+ title = {Parallel reactive molecular dynamics: Numerical methods and algorithmic techniques},
+ journal = {Parallel Computing},
+ year =    2012,
+ volume =  38,
+ pages =   {245--259}
+}
+
+fix qeq/reax command:
+
+@Article{Aktulga12,
+ author = {H. M. Aktulga, J. C. Fogarty, S. A. Pandit, A. Y. Grama},
+ title = {Parallel reactive molecular dynamics: Numerical methods and algorithmic techniques},
+ journal = {Parallel Computing},
+ year =    2012,
+ volume =  38,
+ pages =   {245--259}
+}
+

requirements.txt

@@ -1,5 +1,5 @@
-numpy==1.20.1
-scipy==1.6.1
+numpy==1.19.5
+scipy==1.3.0
 networkx==2.5
 spglib==1.16.1
 pymatgen==2021.2.16

src/io/lammps/lib.rs

@@ -621,6 +621,8 @@ impl<P: Potential> Lammps<P>
             )?;
         }
 
+        let InitInfo { masses, pair_style, pair_coeffs } = init_info;
+
         lmp.command(format!("package omp {}", builder.openmp_threads.unwrap_or(0)))?;
         lmp.command(format!("processors {}", {
             builder.processors.iter()
@@ -631,14 +633,27 @@ impl<P: Potential> Lammps<P>
                 .collect::<Vec<_>>().join(" ")
         }))?;
 
-        lmp.commands(&[
-            "units metal",                  // Angstroms, picoseconds, eV
-            "neigh_modify delay 0",         // disable delay for a safer `run pre no`
-            "atom_style atomic",            // attributes to store per-atom
-            "thermo_modify lost error",     // don't let atoms disappear without telling us
-            "atom_modify map array",        // store all positions in an array
-            "atom_modify sort 0 0.0",       // don't reorder atoms during simulation
-        ])?;
+        if &(pair_style.to_string())[..17] == "pair_style reax/c"{
+            lmp.commands(&[
+                "units real",                   // Angstroms, picoseconds, Kcal/mol
+                "neigh_modify delay 0",         // disable delay for a safer `run pre no`
+                "atom_style charge",            // attributes to store per-atom and required q value
+                "thermo_modify lost error",     // don't let atoms disappear without telling us
+                "atom_modify map array",        // store all positions in an array
+                "atom_modify sort 0 0.0",       // don't reorder atoms during simulation
+                "neighbor 6 bin"                // give more room to the neighbor lists for dynamic charge equalibriations
+            ])?;
+        }
+        else{
+            lmp.commands(&[
+                "units metal",                   // Angstroms, picoseconds, eV
+                "neigh_modify delay 0",         // disable delay for a safer `run pre no`
+                "atom_style atomic",            // attributes to store per-atom
+                "thermo_modify lost error",     // don't let atoms disappear without telling us
+                "atom_modify map array",        // store all positions in an array
+                "atom_modify sort 0 0.0",       // don't reorder atoms during simulation
+            ])?;
+        }
 
         if let Some(_) = molecule_ids {
             lmp.commands(&["fix RSP2_HasMolIds all property/atom mol ghost yes"])?;
@@ -652,7 +667,6 @@ impl<P: Potential> Lammps<P>
         ])?;
 
         {
-            let InitInfo { masses, pair_style, pair_coeffs } = init_info;
 
             lmp.command(format!("create_box {} sim", masses.len()))?;
             for (i, mass) in (1..).zip(masses) {
@@ -661,6 +675,9 @@ impl<P: Potential> Lammps<P>
 
             lmp.command(pair_style.to_string())?;
             lmp.commands(pair_coeffs)?;
+            if &(pair_style.to_string())[..17] == "pair_style reax/c"{
+                lmp.command("compute reax all pair reax/c".to_string())?;
+            }
         }
 
         // HACK:
@@ -699,12 +716,20 @@ impl<P: Potential> Lammps<P>
             }
         }
 
+        // FIXME: the charges of 4 and -2 are only valid for TMD's. 
+        // If other ReaxFF approaches are used in the future, these will need to be read element-wise.
+        if &(pair_style.to_string())[..17] == "pair_style reax/c"{
+            lmp.command("set type 1 charge 4".to_string())?;
+            lmp.command("set type 2 charge -2".to_string())?;
+            lmp.command("fix 1 all qeq/reax 1 0 10 1.0e-6 reax/c".to_string())?;
+        }
+
         // set up computes
         lmp.commands(&[
             &format!("compute RSP2_PE all pe"),
             &format!("compute RSP2_Pressure all pressure NULL virial"),
         ])?;
-
+        
         lmp
     })}
 }

src/io/structure/assemble.rs

@@ -11,7 +11,7 @@
 
 use crate::{FailResult};
 
-use rsp2_structure::{CoordsKind, Lattice, Coords};
+use rsp2_structure::{CoordsKind, Lattice, Coords, Element};
 
 use rsp2_soa_ops::{Part, Perm, Permute};
 use rsp2_array_types::{M33, V3, dot};
@@ -41,6 +41,8 @@ pub struct Assemble {
     /// between the first layer encountered on either side of the boundary)
     pub vacuum_sep: f64,
     layer_seps: Vec<f64>, // [layer] (length: nlayer - 1)
+    //elements of the lattice
+    pub elements: Vec<Element>
 }
 
 impl Assemble {
@@ -163,14 +165,18 @@ pub struct RawAssemble {
     /// This is to help catch bugs related to accidentally wrapping `carts_along_normal`
     /// across a periodic boundary.
     pub check_intralayer_distance: Option<f64>,
+
+    //elements of the lattice
+    pub elements: Vec<Element>
+
 }
 
 impl Assemble {
     pub fn from_raw(raw: RawAssemble) -> FailResult<Self> {
         let RawAssemble {
             normal_axis, lattice, fracs_in_plane, carts_along_normal,
             initial_vacuum_sep, initial_layer_seps, initial_scale,
-            check_intralayer_distance, part,
+            check_intralayer_distance, part, elements
         } = raw;
         assert!(normal_axis < 3);
 
@@ -259,7 +265,7 @@ impl Assemble {
         let layer_seps = initial_layer_seps;
         Ok(Assemble {
             normal_axis, scale, lattice, fracs_in_plane, carts_along_normal,
-            vacuum_sep, layer_seps, perm,
+            vacuum_sep, layer_seps, perm, elements
         })
     }
 }

src/io/structure/layers_yaml.rs

@@ -15,7 +15,7 @@
 use crate::{FailResult, FailOk};
 use crate::assemble::{RawAssemble, Assemble};
 
-use rsp2_structure::{CoordsKind, Lattice, Coords};
+use rsp2_structure::{CoordsKind, Lattice, Coords, Element};
 use std::io::Read;
 
 use rsp2_array_types::{M22, M33, mat, V2, V3, inv, Unvee};
@@ -79,6 +79,7 @@ mod cereal {
         // Number of unique images to generate along each layer lattice vector
         #[serde(default = "defaults::layer::repeat")]
         pub repeat: [u32; 2],
+        pub elements: Vec<String>,
         // Common translation for all positions in layer.
         // NOTE: units of layer lattice
         #[serde(default = "defaults::layer::shift")]
@@ -90,7 +91,7 @@ mod cereal {
 
         pub fn a() -> f64 { 1.0 }
         pub fn vacuum_sep() -> f64 { 10.0 }
-        pub fn layer_sep() -> Either<f64, Vec<f64>> { Either::A(1.0) }
+        pub fn layer_sep() -> Either<f64, Vec<f64>> { Either::A(1.5) }
 
         pub mod layer {
             use super::*;
@@ -123,6 +124,7 @@ mod middle {
         pub transform: M33,
         pub repeat: [u32; 3],
         pub shift: V3,
+        pub elements: Vec<Element>,
     }
 }
 
@@ -148,7 +150,7 @@ fn interpret_cereal(cereal: self::cereal::Root) -> FailResult<middle::Layers>
         let self::cereal::Layer {
             frac_lattice, frac_sites,
             cart_lattice, cart_sites,
-            transform, repeat, shift,
+            transform, repeat, shift, elements
         } = layer;
 
         #[derive(Debug, Copy, Clone, PartialEq, Eq)]
@@ -181,11 +183,15 @@ fn interpret_cereal(cereal: self::cereal::Root) -> FailResult<middle::Layers>
                 }
             }
         };
-
+        let mut els = vec![];
+        for atom in elements{
+            els.push(Element::from_symbol(&atom).unwrap());
+        }
+        let elements = els;
         let transform = m22_to_m33(&transform);
         let shift = V3([shift[0], shift[1], 0.0]);
         let repeat = [repeat[0], repeat[1], 1];
-        middle::Layer { cart_lattice, frac_lattice, cart_sites, transform, repeat, shift }
+        middle::Layer { cart_lattice, frac_lattice, cart_sites, transform, repeat, elements, shift }
     })}).collect::<FailResult<Vec<_>>>()?;
 
     middle::Layers { lattice_a, full_lattice, layers, layer_seps, vacuum_sep }
@@ -199,9 +205,12 @@ fn assemble_from_cereal(cereal: self::cereal::Root) -> FailResult<Assemble>
     } = interpret_cereal(cereal)?;
 
     let mut fracs_in_plane = vec![];
+    let mut elements = vec![];
+
     for layer in layers.into_iter() {
         let lattice = layer.cart_lattice.clone();
         let sites = layer.cart_sites.clone();
+        let atoms = layer.elements.clone();
 
         let mut structure = Coords::new(lattice, CoordsKind::Carts(sites));
         structure.translate_frac(&layer.shift);
@@ -217,18 +226,21 @@ fn assemble_from_cereal(cereal: self::cereal::Root) -> FailResult<Assemble>
         //       to get the integer sc matrices, and somehow verify that the 'repeat' field
         //       is correct.  Or just ignore the 'repeat' field and do HNF reduction to find
         //       a set of periods (but that feels wasteful).
-        let (superstructure, _) = rsp2_structure::supercell::diagonal(layer.repeat).build(&structure);
+        let (superstructure, sc) = rsp2_structure::supercell::diagonal(layer.repeat).build(&structure);
 
         // put them in frac coords for the full lattice
         let mut superstructure = Coords::new(
             Lattice::new(&full_lattice),
             CoordsKind::Carts(superstructure.to_carts()),
         );
+
         // FIXME this reduction is just a bandaid for the above-mentioned issue.
         //       (taking unique positions in the diagonal layer supercells and mapping
         //        them into the cell that we generally use for the structure)
         superstructure.reduce_positions();
         fracs_in_plane.push(superstructure.to_fracs());
+        elements.extend(sc.replicate(&atoms));
+
     }
 
     let carts_along_normal = {
@@ -247,6 +259,7 @@ fn assemble_from_cereal(cereal: self::cereal::Root) -> FailResult<Assemble>
         initial_vacuum_sep: vacuum_sep,
         check_intralayer_distance: None,
         part: None,
+        elements: elements,
     };
 
     // FIXME: some of the possible errors produced by `from_raw` here are

src/potentials/rebo/nonreactive.rs

@@ -159,10 +159,10 @@ pub const SITE_MAX_BONDS: usize = 4;
 #[derive(Debug, Copy, Clone, PartialEq, Eq, PartialOrd, Ord, Hash, enum_map::Enum)]
 pub enum AtomType { Carbon, Hydrogen }
 impl AtomType {
-    pub fn char(self) -> char {
+    pub fn element(self) -> Element {
         match self {
-            AtomType::Carbon => 'C',
-            AtomType::Hydrogen => 'H',
+            AtomType::Carbon => Element::from_symbol("C").unwrap(),
+            AtomType::Hydrogen => Element::from_symbol("H").unwrap(),
         }
     }