diff --git a/.gitignore b/.gitignore index b426b3bb..8b2aab34 100644 --- a/.gitignore +++ b/.gitignore @@ -17,4 +17,7 @@ log.lammps crates.d/config.toml +log.cite +etc/profile.d/lammps.csh + venv diff --git a/build.rs b/build.rs index 1cc09974..b65ba533 100644 --- a/build.rs +++ b/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(()) +} diff --git a/etc/profile.d/lammps.csh b/etc/profile.d/lammps.csh new file mode 100644 index 00000000..ef47dbe5 --- /dev/null +++ b/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 diff --git a/etc/profile.d/lammps.sh b/etc/profile.d/lammps.sh new file mode 100644 index 00000000..01d01bfd --- /dev/null +++ b/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 diff --git a/log.cite b/log.cite new file mode 100644 index 00000000..9aa702cc --- /dev/null +++ b/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} +} + diff --git a/requirements.txt b/requirements.txt index 5fc4f362..50876811 100644 --- a/requirements.txt +++ b/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 diff --git a/src/io/lammps/lib.rs b/src/io/lammps/lib.rs index ff332d2c..6228ac74 100644 --- a/src/io/lammps/lib.rs +++ b/src/io/lammps/lib.rs @@ -621,6 +621,8 @@ impl Lammps

)?; } + 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 Lammps

.collect::>().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 Lammps

])?; { - 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 Lammps

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 Lammps

} } + // 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 })} } diff --git a/src/io/structure/assemble.rs b/src/io/structure/assemble.rs index 90eef9a7..9d8b6da2 100644 --- a/src/io/structure/assemble.rs +++ b/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, // [layer] (length: nlayer - 1) + //elements of the lattice + pub elements: Vec } impl Assemble { @@ -163,6 +165,10 @@ 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, + + //elements of the lattice + pub elements: Vec + } impl Assemble { @@ -170,7 +176,7 @@ impl Assemble { 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 }) } } diff --git a/src/io/structure/layers_yaml.rs b/src/io/structure/layers_yaml.rs index e87c7354..350ff434 100644 --- a/src/io/structure/layers_yaml.rs +++ b/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, // 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> { Either::A(1.0) } + pub fn layer_sep() -> Either> { 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, } } @@ -148,7 +150,7 @@ fn interpret_cereal(cereal: self::cereal::Root) -> FailResult 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 } } }; - + 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::>>()?; middle::Layers { lattice_a, full_lattice, layers, layer_seps, vacuum_sep } @@ -199,9 +205,12 @@ fn assemble_from_cereal(cereal: self::cereal::Root) -> FailResult } = 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 // 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 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 diff --git a/src/potentials/rebo/nonreactive.rs b/src/potentials/rebo/nonreactive.rs index 6872b7cd..ffd06af2 100644 --- a/src/potentials/rebo/nonreactive.rs +++ b/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(), } } diff --git a/src/tasks/cmd/mod.rs b/src/tasks/cmd/mod.rs index 977b2093..c5ac88b1 100644 --- a/src/tasks/cmd/mod.rs +++ b/src/tasks/cmd/mod.rs @@ -135,7 +135,7 @@ impl TrialDir { (_, Some(_)) => {}, } - let pot = PotentialBuilder::from_root_config(Some(&self), on_demand, &settings)?; + let pot = ::from_root_config(Some(&self), on_demand, &settings)?; let (optimizable_coords, mut meta) = { read_optimizable_structure( @@ -1122,7 +1122,7 @@ pub(crate) fn run_shear_plot( let EnergySurfaceArgs { density, extend_border, layer: translated_layer } = plot_args; - let pot = PotentialBuilder::from_config_parts( + let pot = ::from_config_parts( None, on_demand, &settings.threading, @@ -1250,7 +1250,7 @@ impl TrialDir { stored: StoredStructure, ) -> FailResult<()> {Ok({ - let pot = PotentialBuilder::from_root_config(Some(&self), on_demand, &settings)?; + let pot = ::from_root_config(Some(&self), on_demand, &settings)?; let phonons_settings = match &settings.phonons { Some(x) => x, @@ -1357,7 +1357,7 @@ pub(crate) fn run_plot_vdw( update_style: cfg::LammpsUpdateStyle::Fast { sync_positions_every: 1 }.into(), processor_axis_mask: [true; 3].into(), }; - let pot = PotentialBuilder::from_config_parts(None, on_demand, &threading, &lammps, pot)?; + let pot = ::from_config_parts(None, on_demand, &threading, &lammps, pot)?; let lattice = { let a = rs.iter().fold(0.0, |a, &b| f64::max(a, b)) + 20.0; @@ -1410,7 +1410,7 @@ pub(crate) fn run_converge_vdw( update_style: cfg::LammpsUpdateStyle::Fast { sync_positions_every: 1 }.into(), processor_axis_mask: [true; 3].into(), }; - let pot = PotentialBuilder::from_config_parts(None, on_demand, &threading, &lammps, pot)?; + let pot = ::from_config_parts(None, on_demand, &threading, &lammps, pot)?; let lattice = Lattice::orthorhombic(40.0, 40.0, 40.0); let direction = { @@ -1473,7 +1473,7 @@ pub(crate) fn run_single_force_computation( })).fold_ok()? }); - let pot = PotentialBuilder::from_root_config(None, on_demand, &settings)?; + let pot = ::from_root_config(None, on_demand, &settings)?; pot.one_off().compute_force(&coords, meta.sift())? })} @@ -1531,7 +1531,7 @@ pub(crate) fn run_layer_mode_frequencies( OnlyUniqueResult::NoValues => unreachable!(), }; - let pot = PotentialBuilder::from_root_config(None, on_demand, &settings)?; + let pot = ::from_root_config(None, on_demand, &settings)?; let mut diff_fn = pot.initialize_diff_fn(&original_coords, meta.sift())?; warn!("Don't quote these values, only compare them! (I'm not sure about the prefactor...)"); @@ -1587,7 +1587,7 @@ pub(crate) fn run_dynmat_at_q( qpoint_frac: V3, structure: StoredStructure, ) -> FailResult { - let pot = PotentialBuilder::from_root_config(None, on_demand, &settings)?; + let pot = ::from_root_config(None, on_demand, &settings)?; let phonons_settings = match &settings.phonons { Some(x) => x, @@ -1644,7 +1644,7 @@ impl TrialDir { None => bail!("`rsp2-run-after-diagonalization` cannot be used without a `phonons:` config section"), }; - let pot = PotentialBuilder::from_root_config(Some(&self), on_demand, &settings)?; + let pot = ::from_root_config(Some(&self), on_demand, &settings)?; let (coords, meta) = self.read_stored_structure_data(&self.structure_path(PreEvChase(prev_iteration)))?; @@ -1881,7 +1881,7 @@ pub(crate) fn read_optimizable_structure( }); out_layers = Some(layer_builder.atom_layers().into_iter().map(Layer).collect::>().into()); - out_elements = vec![CARBON; layer_builder.num_atoms()].into(); + out_elements = layer_builder.elements.clone().into(); out_masses = masses_by_config(mass_cfg, out_elements.clone())?; out_coords = ScalableCoords::KnownLayers { layer_builder }; out_bonds = None; // Determine bonds AFTER parameter optimization. diff --git a/src/tasks/cmd/param_optimization.rs b/src/tasks/cmd/param_optimization.rs index b1ebfa2e..26b4cca3 100644 --- a/src/tasks/cmd/param_optimization.rs +++ b/src/tasks/cmd/param_optimization.rs @@ -13,6 +13,7 @@ use crate::{FailResult}; use crate::potential::{CommonMeta, PotentialBuilder, DiffFn, BondGrad, BondDiffFn}; use crate::meta::{prelude::*}; +use rsp2_tasks_config::MaskBit; use rsp2_minimize::exact_ls::{Value, Golden}; use rsp2_structure::{Lattice, CoordsKind, Coords}; use rsp2_structure::layer::{LayersPerUnitCell, require_simple_axis_normal}; @@ -155,13 +156,23 @@ fn add_scalables( }, ( - &cfg::Scalable::UniformLayerSep { .. }, - &ScalableCoords::UnknownLayers { .. }, - ) | ( - &cfg::Scalable::LayerSeps { .. }, + &cfg::Scalable::UniformLayerSep { .. } | &cfg::Scalable::LayerSeps { .. }, &ScalableCoords::UnknownLayers { .. }, ) => bail!("cannot scale layer separations when layers have not been determined"), + ( + &cfg::Scalable::UniformLayerSep { ref mask, .. } | &cfg::Scalable::LayerSeps { ref mask, .. }, + &ScalableCoords::KnownLayers { ref layer_builder, .. }, + ) => { + if let Some(mask) = mask { + ensure!( + mask.len() == layer_builder.num_layer_seps(), + "a layer-sep mask has the wrong length (expected {} elements, got {:?})", + layer_builder.num_layer_seps(), mask, + ) + } + }, + _ => {}, } @@ -195,15 +206,19 @@ fn add_scalables( }); }, - &cfg::Scalable::UniformLayerSep { ref range } => { + &cfg::Scalable::UniformLayerSep { ref range, ref mask } => { // one scalable for all layers + let n_layer_seps = n_layer_seps.expect("BUG!"); + let mask = mask.clone().unwrap_or(vec![MaskBit(true); n_layer_seps]); emit(Scalable { name: format!("a uniform layer separation"), - setter: Box::new(|s, val| match s { + setter: Box::new(move |s, val| match s { ScalableCoords::UnknownLayers { .. } => unreachable!(), ScalableCoords::KnownLayers { layer_builder, .. } => { - for x in layer_builder.layer_seps() { - *x = val; + for k in 0..n_layer_seps { + if mask[k].0 { + layer_builder.layer_seps()[k] = val; + } } }, }), @@ -211,10 +226,14 @@ fn add_scalables( }); }, - &cfg::Scalable::LayerSeps { ref range } => { + &cfg::Scalable::LayerSeps { ref range, ref mask } => { // separate scalables for each layer let n_layer_seps = n_layer_seps.expect("BUG!"); + let mask = mask.clone().unwrap_or(vec![MaskBit(true); n_layer_seps]); for k in 0..n_layer_seps { + if !mask[k].0 { + continue; + } emit(Scalable { name: format!("layer separation {}", k), setter: Box::new(move |s, val| match s { @@ -303,6 +322,7 @@ impl ScalableCoords { initial_scale: Some([1.0; 3]), part: Some(layers.get_part()), check_intralayer_distance: Some(cfg.threshold), + elements: vec![], }).unwrap(); ScalableCoords::KnownLayers { layer_builder } } diff --git a/src/tasks/config/config.rs b/src/tasks/config/config.rs index 5d5a36c3..243995ce 100644 --- a/src/tasks/config/config.rs +++ b/src/tasks/config/config.rs @@ -316,7 +316,7 @@ pub enum Scalable { range: ScalableRange, }, - /// Optimize a single value shared by all layer separations. + /// Optimize a single value shared by multiple layer separations. /// /// Under certain conditions, the optimum separation IS identical for /// all layers (e.g. generated structures where all pairs of layers @@ -326,6 +326,9 @@ pub enum Scalable { /// is "good enough" that CG can be trusted to take care of the rest. #[serde(rename_all = "kebab-case")] UniformLayerSep { + /// Toggle which separations are affected. For n layers, this must have n-1 elements. + #[serde(default)] + mask: OrDefault>, #[serde(flatten)] range: ScalableRange, }, @@ -333,6 +336,9 @@ pub enum Scalable { /// Optimize each layer separation individually. Can be costly. #[serde(rename_all = "kebab-case")] LayerSeps { + /// Toggle which separations are affected. For n layers, this must have n-1 elements. + #[serde(default)] + mask: OrDefault>, #[serde(flatten)] range: ScalableRange, }, @@ -622,6 +628,7 @@ pub enum PotentialKind { /// Arranges atoms into a chain along the first lattice vector. #[serde(rename = "test-func-chainify")] TestChainify, + } #[derive(Serialize, Deserialize)] @@ -638,6 +645,7 @@ pub enum LammpsPotentialKind { #[serde(rename = "kc-z")] KolmogorovCrespiZ(LammpsPotentialKolmogorovCrespiZ), /// `pair_style hybrid rebo kolmogorov/crespi/full` #[serde(rename = "kc-full")] KolmogorovCrespiFull(LammpsPotentialKolmogorovCrespiFull), + #[serde(rename = "reaxff")] ReaxFF(LammpsPotentialReaxFF), } #[derive(Serialize, Deserialize)] @@ -662,6 +670,12 @@ pub struct LammpsPotentialAirebo { pub omp: OrDefault, } +#[derive(Serialize, Deserialize)] +#[derive(Debug, Clone, PartialEq)] +#[serde(rename_all = "kebab-case")] +pub struct LammpsPotentialReaxFF { +} + #[derive(Serialize, Deserialize)] #[derive(Debug, Clone, PartialEq)] #[serde(rename_all = "kebab-case")] diff --git a/src/tasks/filetypes/eigensols.rs b/src/tasks/filetypes/eigensols.rs index 0ab2113e..043b7b5f 100644 --- a/src/tasks/filetypes/eigensols.rs +++ b/src/tasks/filetypes/eigensols.rs @@ -16,7 +16,9 @@ use crate::traits::{Load, AsPath, save::Json}; // Conversion factor phonopy uses to scale the eigenvalues to THz angular momentum. // = sqrt(eV/amu)/angstrom/(2*pi)/THz -const SQRT_EIGENVALUE_TO_THZ: f64 = 15.6333043006705; +// 1 ev = 23.061 Kcal/mol (NIST) +// 4.80218700177 sqrt((Kcal/mol)/ev) ~= sqrt(23.061) +const SQRT_EIGENVALUE_TO_THZ: f64 = 15.6333043006705 * 4.80218700177; // = THz / (c / cm) const THZ_TO_WAVENUMBER: f64 = 33.3564095198152; const SQRT_EIGENVALUE_TO_WAVENUMBER: f64 = SQRT_EIGENVALUE_TO_THZ * THZ_TO_WAVENUMBER; diff --git a/src/tasks/lib.rs b/src/tasks/lib.rs index 4944b4b9..97ca1c0c 100644 --- a/src/tasks/lib.rs +++ b/src/tasks/lib.rs @@ -143,6 +143,8 @@ mod common { match elem { consts::HYDROGEN => 1.00794, consts::CARBON => 12.0107, + consts::MOLYBDENUM => 95.56, + consts::SULFUR => 32.06, _ => bail!("No default mass for element {}.", elem.symbol()), } }))} diff --git a/src/tasks/potential/lammps.rs b/src/tasks/potential/lammps.rs index 1b49cc46..3f95b883 100644 --- a/src/tasks/potential/lammps.rs +++ b/src/tasks/potential/lammps.rs @@ -495,6 +495,57 @@ mod airebo { } } +pub use self::reaxff::ReaxFF; +mod reaxff { + use super::*; + + #[derive(Debug, Clone)] + pub enum ReaxFF { + ReaxFF {}, + } + + impl<'a> From<&'a cfg::LammpsPotentialReaxFF> for ReaxFF { + fn from(cfg: &'a cfg::LammpsPotentialReaxFF) -> Self { + let cfg::LammpsPotentialReaxFF { .. } = *cfg; + ReaxFF::ReaxFF{} + } + } + + + impl LammpsPotential for ReaxFF { + type Meta = CommonMeta; + + fn molecule_ids(&self, _: &Coords, _: &CommonMeta) -> Option> { None } + + fn atom_types(&self, _: &Coords, meta: &CommonMeta) -> Vec + { + let elements: meta::SiteElements = meta.pick(); + elements.iter().map(|elem| match elem.symbol() { + "Mo" => AtomType::new(1), + "S" => AtomType::new(2), + sym => panic!("Unexpected element in ReaxFF: {}", sym), + }).collect() + } + + fn init_info(&self, _: &Coords, meta: &CommonMeta) -> InitInfo { + let masses = vec![ + only_unique_mass(meta, consts::MOLYBDENUM), + only_unique_mass(meta, consts::SULFUR), + ]; + let pair_style; + match *self { + ReaxFF::ReaxFF {..} => { + pair_style = PairStyle::named("reax/c NULL safezone 5 mincap 100"); //safezone and mincap extend memory segments + }, + }; + let pair_coeffs = vec![ + PairCoeff::new(.., ..).args(&["ffield.reax.mos2_old", "Mo", "S"]), + ]; + InitInfo { masses, pair_style, pair_coeffs } + } + } +} + pub use self::kc_z::KolmogorovCrespiZ; mod kc_z { use super::*; diff --git a/src/tasks/potential/mod.rs b/src/tasks/potential/mod.rs index dbe644df..9b926d54 100644 --- a/src/tasks/potential/mod.rs +++ b/src/tasks/potential/mod.rs @@ -582,9 +582,9 @@ impl dyn PotentialBuilder { cfg::PotentialKind::Lammps(_) => { assert!(!found_lammps, "(BUG!) more than one lammps potential after validation!?"); found_lammps = true; - PotentialBuilder::single_from_config_parts(trial_dir, on_demand.take(), threading, lammps, &cfg) + ::single_from_config_parts(trial_dir, on_demand.take(), threading, lammps, &cfg) }, - _ => PotentialBuilder::single_from_config_parts(trial_dir, None, threading, lammps, &cfg), + _ => ::single_from_config_parts(trial_dir, None, threading, lammps, &cfg), } }) .collect::>>()? @@ -634,6 +634,11 @@ impl dyn PotentialBuilder { let pot = self::lammps::Builder::new(trial_dir, on_demand, threading, lammps, lammps_pot)?; Ok(Box::new(pot)) }, + cfg::LammpsPotentialKind::ReaxFF(cfg) => { + let lammps_pot = self::lammps::ReaxFF::from(cfg); + let pot = self::lammps::Builder::new(trial_dir, on_demand, threading, lammps, lammps_pot)?; + Ok(Box::new(pot)) + }, }, cfg::PotentialKind::KolmogorovCrespi(cfg) => { let cfg = cfg.clone();