option to remove molecular internal interaction added

Project.toml

@@ -1,7 +1,7 @@
 name = "ACEds"
 uuid = "286ca7f0-9926-4809-930a-857c1b5bb5cd"
 authors = ["Christoph Ortner <[email protected]> and Matthias Sachs <[email protected]>"]
-version = "0.1.6"
+version = "0.1.7"
 
 [deps]
 ACE = "3e8ccfd2-c8b0-11ea-32f1-f3a5990fd77a"

examples/fit-ACEFrictionModel-new-ac.jl

@@ -35,24 +35,25 @@ data = Dict("train" => rdata[1:n_train], "test"=> rdata[n_train+1:end]);
 
 species_friction = [:H]
 species_env = [:Cu]
+species_mol = [:H]
 rcut = 8.0
 coupling= RowCoupling()
-m_inv = ac_matrixmodel(ACE.Invariant(),species_friction,species_env, coupling; n_rep = 2, rcut_on = rcut, rcut_off = rcut, maxorder_on=2, maxdeg_on=5,
+m_inv = ac_matrixmodel(ACE.Invariant(),species_friction,species_env, coupling,species_mol; n_rep = 2, rcut_on = rcut, rcut_off = rcut, maxorder_on=2, maxdeg_on=5,
         species_maxorder_dict_on = Dict( :H => 1), 
         species_weight_cat_on = Dict(:H => .75, :Cu=> 1.0),
         species_maxorder_dict_off = Dict( :H => 0), 
         species_weight_cat_off = Dict(:H => 1.0, :Cu=> 1.0),
         bond_weight = .5
     );
-m_cov = ac_matrixmodel(ACE.EuclideanVector(Float64),species_friction,species_env, coupling; n_rep=3, rcut_on = rcut, rcut_off = rcut, maxorder_on=2, maxdeg_on=5,
+m_cov = ac_matrixmodel(ACE.EuclideanVector(Float64),species_friction,species_env, coupling,species_mol; n_rep=3, rcut_on = rcut, rcut_off = rcut, maxorder_on=2, maxdeg_on=5,
         species_maxorder_dict_on = Dict( :H => 1), 
         species_weight_cat_on = Dict(:H => .75, :Cu=> 1.0),
         species_maxorder_dict_off = Dict( :H => 0), 
         species_weight_cat_off = Dict(:H => 1.0, :Cu=> 1.0),
         bond_weight = .5
     );
 
-m_equ = ac_matrixmodel(ACE.EuclideanMatrix(Float64),species_friction,species_env, coupling; n_rep=2, rcut_on = rcut, rcut_off = rcut, maxorder_on=2, maxdeg_on=5,
+m_equ = ac_matrixmodel(ACE.EuclideanMatrix(Float64),species_friction,species_env, coupling,species_mol; n_rep=2, rcut_on = rcut, rcut_off = rcut, maxorder_on=2, maxdeg_on=5,
         species_maxorder_dict_on = Dict( :H => 1), 
         species_weight_cat_on = Dict(:H => .75, :Cu=> 1.0),
         species_maxorder_dict_off = Dict( :H => 0), 

examples/fit-ACEFrictionModel-new-pw-ellipsoid.jl

@@ -38,11 +38,12 @@ using ACEds.AtomCutoffs: SphericalCutoff
 using ACEds.MatrixModels: NoZ2Sym, SpeciesUnCoupled
 species_friction = [:H]
 species_env = [:Cu]
+species_mol = [:H]
 rcut = 8.0
 z2sym= NoZ2Sym()
 speciescoupling = SpeciesUnCoupled()
 cutoff = EllipsoidCutoff(3.5,4.0,6.0)
-m_inv = pwc_matrixmodel(ACE.Invariant(),species_friction,species_env, z2sym,  speciescoupling;
+m_inv = pwc_matrixmodel(ACE.Invariant(),species_friction,species_env, z2sym,  speciescoupling, species_mol;
         n_rep = 3,
         maxorder_off=2, 
         maxdeg_off=5, 
@@ -53,7 +54,7 @@ m_inv = pwc_matrixmodel(ACE.Invariant(),species_friction,species_env, z2sym,  sp
         species_weight_cat_off = Dict(:H => 1.0, :Cu=> 1.0),
         bond_weight = .25
     );
-m_cov = pwc_matrixmodel(ACE.EuclideanVector(Float64),species_friction,species_env, z2sym,  speciescoupling;
+m_cov = pwc_matrixmodel(ACE.EuclideanVector(Float64),species_friction,species_env, z2sym,  speciescoupling, species_mol;
         n_rep = 3,
         maxorder_off=2, 
         maxdeg_off=5, 
@@ -64,7 +65,7 @@ m_cov = pwc_matrixmodel(ACE.EuclideanVector(Float64),species_friction,species_en
         species_weight_cat_off = Dict(:H => 1.0, :Cu=> 1.0),
         bond_weight = .25
     );
-m_equ = pwc_matrixmodel(ACE.EuclideanMatrix(Float64),species_friction,species_env, z2sym,  speciescoupling;
+m_equ = pwc_matrixmodel(ACE.EuclideanMatrix(Float64),species_friction,species_env, z2sym,  speciescoupling, species_mol;
         n_rep = 3,
         maxorder_off=2, 
         maxdeg_off=5, 
@@ -76,15 +77,15 @@ m_equ = pwc_matrixmodel(ACE.EuclideanMatrix(Float64),species_friction,species_en
         bond_weight = .25
     );
 
-m_inv0 = onsiteonly_matrixmodel(ACE.Invariant(), species_friction, species_env; id=:inv0, n_rep = 3, rcut_on = rcut, maxorder_on=2, maxdeg_on=5,
+m_inv0 = onsiteonly_matrixmodel(ACE.Invariant(), species_friction, species_env, species_mol; id=:inv0, n_rep = 3, rcut_on = rcut, maxorder_on=2, maxdeg_on=5,
     species_maxorder_dict_on = Dict( :H => 1), 
     species_weight_cat_on = Dict(:H => .75, :Cu=> 1.0)
     );
-m_cov0 = onsiteonly_matrixmodel(ACE.EuclideanVector(Float64), species_friction, species_env; id=:cov0, n_rep = 3, rcut_on = rcut, maxorder_on=2, maxdeg_on=5,
+m_cov0 = onsiteonly_matrixmodel(ACE.EuclideanVector(Float64), species_friction, species_env, species_mol; id=:cov0, n_rep = 3, rcut_on = rcut, maxorder_on=2, maxdeg_on=5,
     species_maxorder_dict_on = Dict( :H => 1), 
     species_weight_cat_on = Dict(:H => .75, :Cu=> 1.0)
     );
-m_equ0 = onsiteonly_matrixmodel(ACE.EuclideanMatrix(Float64), species_friction, species_env; id=:equ0, n_rep = 3, rcut_on = rcut, maxorder_on=2, maxdeg_on=5,
+m_equ0 = onsiteonly_matrixmodel(ACE.EuclideanMatrix(Float64), species_friction, species_env, species_mol; id=:equ0, n_rep = 3, rcut_on = rcut, maxorder_on=2, maxdeg_on=5,
     species_maxorder_dict_on = Dict( :H => 1), 
     species_weight_cat_on = Dict(:H => .75, :Cu=> 1.0)
     );

examples/fit-ACEFrictionModel-new-pw.jl

@@ -37,10 +37,11 @@ using ACEds.AtomCutoffs: SphericalCutoff
 using ACEds.MatrixModels: NoZ2Sym, SpeciesUnCoupled
 species_friction = [:H]
 species_env = [:Cu]
+species_mol = [:H]
 rcut = 8.0
 z2sym= NoZ2Sym()
 speciescoupling = SpeciesUnCoupled()
-m_inv = pwc_matrixmodel(ACE.Invariant(),species_friction,species_env, z2sym,  speciescoupling;
+m_inv = pwc_matrixmodel(ACE.Invariant(),species_friction, species_env, z2sym,  speciescoupling, species_mol;
         n_rep = 3,
         maxorder_off=2, 
         maxdeg_off=5, 
@@ -51,7 +52,7 @@ m_inv = pwc_matrixmodel(ACE.Invariant(),species_friction,species_env, z2sym,  sp
         species_weight_cat_off = Dict(:H => 1.0, :Cu=> 1.0),
         bond_weight = .5
     );
-m_cov = pwc_matrixmodel(ACE.EuclideanVector(Float64),species_friction,species_env, z2sym,  speciescoupling;
+m_cov = pwc_matrixmodel(ACE.EuclideanVector(Float64),species_friction, species_env, z2sym,  speciescoupling, species_mol;
         n_rep = 3,
         maxorder_off=2, 
         maxdeg_off=5, 
@@ -62,7 +63,7 @@ m_cov = pwc_matrixmodel(ACE.EuclideanVector(Float64),species_friction,species_en
         species_weight_cat_off = Dict(:H => 1.0, :Cu=> 1.0),
         bond_weight = .5
     );
-m_equ = pwc_matrixmodel(ACE.EuclideanMatrix(Float64),species_friction,species_env, z2sym,  speciescoupling;
+m_equ = pwc_matrixmodel(ACE.EuclideanMatrix(Float64),species_friction, species_env, z2sym, speciescoupling, species_mol;
         n_rep = 3,
         maxorder_off=2, 
         maxdeg_off=5, 
@@ -74,15 +75,15 @@ m_equ = pwc_matrixmodel(ACE.EuclideanMatrix(Float64),species_friction,species_en
         bond_weight = .5
     );
 
-m_inv0 = onsiteonly_matrixmodel(ACE.Invariant(), species_friction, species_env; id=:inv0, n_rep = 3, rcut_on = rcut, maxorder_on=2, maxdeg_on=5,
+m_inv0 = onsiteonly_matrixmodel(ACE.Invariant(), species_friction, species_env, species_mol; id=:inv0, n_rep = 3, rcut_on = rcut, maxorder_on=2, maxdeg_on=5,
     species_maxorder_dict_on = Dict( :H => 1), 
     species_weight_cat_on = Dict(:H => .75, :Cu=> 1.0)
     );
-m_cov0 = onsiteonly_matrixmodel(ACE.EuclideanVector(Float64), species_friction, species_env; id=:cov0, n_rep = 3, rcut_on = rcut, maxorder_on=2, maxdeg_on=5,
+m_cov0 = onsiteonly_matrixmodel(ACE.EuclideanVector(Float64), species_friction, species_env, species_mol; id=:cov0, n_rep = 3, rcut_on = rcut, maxorder_on=2, maxdeg_on=5,
     species_maxorder_dict_on = Dict( :H => 1), 
     species_weight_cat_on = Dict(:H => .75, :Cu=> 1.0)
     );
-m_equ0 = onsiteonly_matrixmodel(ACE.EuclideanMatrix(Float64), species_friction, species_env; id=:equ0, n_rep = 3, rcut_on = rcut, maxorder_on=2, maxdeg_on=5,
+m_equ0 = onsiteonly_matrixmodel(ACE.EuclideanMatrix(Float64), species_friction, species_env, species_mol; id=:equ0, n_rep = 3, rcut_on = rcut, maxorder_on=2, maxdeg_on=5,
     species_maxorder_dict_on = Dict( :H => 1), 
     species_weight_cat_on = Dict(:H => .75, :Cu=> 1.0)
     );

src/matrixmodels/matrixmodels.jl

@@ -35,7 +35,8 @@ using ACEbonds.BondCutoffs: AbstractBondCutoff
 using ACEds.AtomCutoffs: SphericalCutoff
 using ACE
 using ACEds.MatrixModels
-import ACEbonds: SymmetricEllipsoidBondBasis
+#import ACEbonds: SymmetricEllipsoidBondBasis
+include("../patches/acebonds_basisselectors.jl")
 using ACEds
 using JuLIP: AtomicNumber
 
@@ -399,6 +400,23 @@ _z2couplingToString(::Odd) = "Covariant"
 _cutoff(cutoff::SphericalCutoff) = cutoff.r_cut
 _cutoff(cutoff::EllipsoidCutoff) = cutoff.r_cut
 
+"""
+`NoMolOnly`: selects all basis functions which model interactions between atoms of the molecule only. Use this filter if the molecule feels only
+friction if in contact to the substrat.   
+"""
+struct NoMolOnly
+      isym::Symbol
+      categories
+end
+
+function (f::NoMolOnly)(bb) 
+      if isempty(bb)
+            return true
+      else
+            return !all([getproperty(b, f.isym) in f.categories for b in bb])
+      end
+end
+
 function offsite_linbasis(property,species;
     z2symmetry = NoZ2Sym(), 
     maxorder = 2,
@@ -415,9 +433,16 @@ function offsite_linbasis(property,species;
     species_minorder_dict = Dict{Any, Float64}(),
     species_maxorder_dict = Dict{Any, Float64}(),
     species_weight_cat = Dict(c => 1.0 for c in species),
+    molspecies = []
     )
+    if isempty(molspecies)
+        filterfun = _ -> true
+    else
+        filterfun = NoMolOnly(:mube, vcat(molspecies,:bond))
+    end 
+
     @info "Generate offsite basis"
-    @time offsite = SymmetricEllipsoidBondBasis(property; 
+    @time offsite = SymmetricEllipsoidBondBasis2(property; 
                 r0 = r0_ratio, 
                 rin = rin_ratio, 
                 pcut = pcut, 
@@ -433,7 +458,8 @@ function offsite_linbasis(property,species;
                 bondsymmetry=_z2couplingToString(z2symmetry),
                 species=species, 
                 isym=isym, 
-                bond_weight = bond_weight,  
+                bond_weight = bond_weight,
+                filterfun = filterfun
     )
     @info "Size of offsite basis elements: $(length(offsite))"
     return BondBasis(offsite,z2symmetry)
@@ -446,7 +472,8 @@ function onsite_linbasis(property,species;
     species_minorder_dict = Dict{Any, Float64}(),
     species_maxorder_dict = Dict{Any, Float64}(),
     weight = Dict(:l => 1.0, :n => 1.0), 
-    species_weight_cat = Dict(c => 1.0 for c in species)    
+    species_weight_cat = Dict(c => 1.0 for c in species),
+    molspecies = []
     )
     @info "Generate onsite basis"
     Bsel = ACE.SparseBasis(; maxorder=maxorder, p = p_sel, default_maxdeg = maxdeg, weight=weight ) 
@@ -470,8 +497,12 @@ function onsite_linbasis(property,species;
         maxorder_dict = species_maxorder_dict, 
         weight_cat = species_weight_cat
     )
-
-    @time onsite = ACE.SymmetricBasis(property, RnYlm * Zk, Bselcat;);
+    if isempty(molspecies)
+        filter = _ -> true
+    else
+        filter = NoMolOnly(:mu, molspecies)
+    end
+    @time onsite = ACE.SymmetricBasis(property, RnYlm * Zk, Bselcat; filterfun=filter);
     @info "Size of onsite basis elements: $(length(onsite))"
     return onsite
 end

src/matrixmodelsutils.jl

@@ -2,7 +2,6 @@ using ACEds.AtomCutoffs: SphericalCutoff
 #using ACEds.Utils: SymmetricBondSpecies_basis
 using ACE
 using ACEds.MatrixModels
-import ACEbonds: SymmetricEllipsoidBondBasis
 using ACEds
 using JuLIP: AtomicNumber
 using ACEds.MatrixModels: _o3symmetry
@@ -11,7 +10,7 @@ using ACEds.PWMatrix: _msort
 using ACEds.MatrixModels: _default_id
 export ac_matrixmodel, mbdpd_matrixmodel, pwc_matrixmodel, onsiteonly_matrixmodel
 
-function ac_matrixmodel(property, species_friction, species_env, coupling=RowCoupling(); 
+function ac_matrixmodel(property, species_friction, species_env, coupling=RowCoupling(), species_mol=[]; 
     id=nothing, 
     n_rep = 3, 
     maxorder_on=2, 
@@ -52,7 +51,8 @@ function ac_matrixmodel(property, species_friction, species_env, coupling=RowCou
         species_minorder_dict = species_minorder_dict_on,
         species_maxorder_dict = species_maxorder_dict_on,
         weight = weight_on, 
-        species_weight_cat = species_weight_cat_on  
+        species_weight_cat = species_weight_cat_on,
+        molspecies = species_mol  
     )
     #@info "Size of onsite basis elements: $(length(onsitebasis))"
 
@@ -73,7 +73,8 @@ function ac_matrixmodel(property, species_friction, species_env, coupling=RowCou
         bond_weight = bond_weight,
         species_minorder_dict = species_minorder_dict_off,
         species_maxorder_dict = species_maxorder_dict_off,
-        species_weight_cat = species_weight_cat_off
+        species_weight_cat = species_weight_cat_off,
+        molspecies = species_mol  
     )
     @info "Size of offsite basis elements: $(length(offsitebasis))"
 
@@ -86,7 +87,7 @@ function ac_matrixmodel(property, species_friction, species_env, coupling=RowCou
     return ACMatrixModel(onsitemodels, offsitemodels, id, coupling)
 end
 
-function onsiteonly_matrixmodel(property, species_friction, species_env; 
+function onsiteonly_matrixmodel(property, species_friction, species_env, species_mol=[]; 
     id=nothing, 
     n_rep = 3, 
     maxorder_on=2, 
@@ -120,7 +121,8 @@ function onsiteonly_matrixmodel(property, species_friction, species_env;
         species_minorder_dict = species_minorder_dict_on,
         species_maxorder_dict = species_maxorder_dict_on,
         weight = weight_on, 
-        species_weight_cat = species_weight_cat_on  
+        species_weight_cat = species_weight_cat_on,
+        molspecies = species_mol  
     )
     #@info "Size of onsite basis elements: $(length(onsitebasis))"
 
@@ -133,7 +135,7 @@ function onsiteonly_matrixmodel(property, species_friction, species_env;
 end
 
 
-function mbdpd_matrixmodel(property, species_friction, species_env; 
+function mbdpd_matrixmodel(property, species_friction, species_env, species_mol=[]; 
     id=nothing, 
     n_rep = 3, 
     maxorder_off=2, 
@@ -169,7 +171,8 @@ function mbdpd_matrixmodel(property, species_friction, species_env;
         bond_weight = bond_weight,
         species_minorder_dict = species_minorder_dict_off,
         species_maxorder_dict = species_maxorder_dict_off,
-        species_weight_cat = species_weight_cat_off
+        species_weight_cat = species_weight_cat_off,
+        molspecies = species_mol
     )
 
     if typeof(cutoff_off)<:AbstractBondCutoff
@@ -184,7 +187,7 @@ function mbdpd_matrixmodel(property, species_friction, species_env;
     return MBDPDMatrixModel(offsitemodels, id)
 end
 
-function pwc_matrixmodel(property, species_friction, species_env, z2sym, speciescoupling; 
+function pwc_matrixmodel(property, species_friction, species_env, z2sym, speciescoupling,species_mol=[]; 
     id=nothing, 
     n_rep = 3, 
     maxorder_off=2, 
@@ -220,7 +223,8 @@ function pwc_matrixmodel(property, species_friction, species_env, z2sym, species
         bond_weight = bond_weight,
         species_minorder_dict = species_minorder_dict_off,
         species_maxorder_dict = species_maxorder_dict_off,
-        species_weight_cat = species_weight_cat_off
+        species_weight_cat = species_weight_cat_off,
+        molspecies = species_mol
     )
 
     if typeof(speciescoupling)<:SpeciesUnCoupled

src/patches/ace_basisselectors.jl

@@ -0,0 +1,17 @@
+
+
+"""
+`EvenL`: selects all basis functions where the sum `L = sum_i l_i` of the degrees `l_i` of the spherical harmonics is even.   
+"""
+struct NoMolOnly
+      isym::Symbol
+      categories
+end
+
+function (f::NoMolPW)(bb) 
+      if isempty(bb)
+            return true
+      else
+            return length(bb)= sum( [getproperty(b, f.isym) in categories for b in bb])
+      end
+end