Feature force lever

corelib/src/libs/SireCAS/lambdaschedule.h

@@ -79,6 +79,9 @@ namespace SireCAS
         static LambdaSchedule standard_morph();
         static LambdaSchedule charge_scaled_morph(double scale = 0.2);
 
+        static LambdaSchedule standard_decouple(bool perturbed_is_decoupled = true);
+        static LambdaSchedule charge_scaled_decouple(double scale = 0.2, bool perturbed_is_decoupled = true);
+
         static SireCAS::Symbol lam();
         static SireCAS::Symbol initial();
         static SireCAS::Symbol final();
@@ -95,6 +98,16 @@ namespace SireCAS
 
         QStringList getLevers() const;
 
+        void addForce(const QString &force);
+        void addForces(const QStringList &forces);
+
+        void removeForce(const QString &force);
+        void removeForces(const QStringList &forces);
+
+        int nForces() const;
+
+        QStringList getForces() const;
+
         int nStages() const;
 
         QStringList getStages() const;
@@ -116,6 +129,8 @@ namespace SireCAS
                          const QString &stage,
                          const SireCAS::Expression &equation);
 
+        void removeStage(const QString &stage);
+
         void addMorphStage();
         void addMorphStage(const QString &name);
 
@@ -124,20 +139,39 @@ namespace SireCAS
                                   const QString &recharge_name,
                                   double scale = 0.2);
 
-        void setEquation(const QString &stage,
-                         const QString &lever,
-                         const SireCAS::Expression &equation);
+        void addDecoupleStage(bool perturbed_is_decoupled = true);
+        void addDecoupleStage(const QString &name, bool perturbed_is_decoupled = true);
+
+        void setDefaultStageEquation(const QString &stage,
+                                     const SireCAS::Expression &equation);
 
-        void setDefaultEquation(const QString &stage,
-                                const SireCAS::Expression &equation);
+        void setEquation(const QString &stage = "*",
+                         const QString &force = "*",
+                         const QString &lever = "*",
+                         const SireCAS::Expression &equation = SireCAS::Expression());
 
-        void removeEquation(const QString &stage,
-                            const QString &lever);
+        void removeEquation(const QString &stage = "*",
+                            const QString &force = "*",
+                            const QString &lever = "*");
 
-        SireCAS::Expression getEquation(const QString &stage) const;
+        bool hasForceSpecificEquation(const QString &stage = "*",
+                                      const QString &force = "*",
+                                      const QString &lever = "*") const;
 
-        SireCAS::Expression getEquation(const QString &stage,
-                                        const QString &lever) const;
+        SireCAS::Expression getEquation(const QString &stage = "*",
+                                        const QString &force = "*",
+                                        const QString &lever = "*") const;
+
+        void setMoleculeSchedule(int pert_mol_id,
+                                 const LambdaSchedule &schedule);
+
+        bool hasMoleculeSchedule(int pert_mol_id) const;
+
+        void removeMoleculeSchedule(int pert_mol_id);
+
+        LambdaSchedule takeMoleculeSchedule(int pert_mol_id);
+
+        const LambdaSchedule &getMoleculeSchedule(int pert_mol_id) const;
 
         QHash<QString, QVector<double>> getLeverValues(const QVector<double> &lambda_values,
                                                        double initial = 1.0,
@@ -160,18 +194,20 @@ namespace SireCAS
 
         SireCAS::Symbol getConstantSymbol(const QString &constant) const;
 
-        double morph(const QString &lever,
-                     double initial, double final, double lambda_value) const;
+        double morph(const QString &force = "*", const QString &lever = "*",
+                     double initial = 0, double final = 1, double lambda_value = 0) const;
 
-        QVector<double> morph(const QString &lever,
-                              const QVector<double> &initial,
-                              const QVector<double> &final,
-                              double lambda_value) const;
+        QVector<double> morph(const QString &force = "*",
+                              const QString &lever = "*",
+                              const QVector<double> &initial = QVector<double>(),
+                              const QVector<double> &final = QVector<double>(),
+                              double lambda_value = 0.0) const;
 
-        QVector<int> morph(const QString &lever,
-                           const QVector<int> &initial,
-                           const QVector<int> &final,
-                           double lambda_value) const;
+        QVector<int> morph(const QString &force = "*",
+                           const QString &lever = "*",
+                           const QVector<int> &initial = QVector<int>(),
+                           const QVector<int> &final = QVector<int>(),
+                           double lambda_value = 0.0) const;
 
         double clamp(double lambda_value) const;
 
@@ -180,9 +216,19 @@ namespace SireCAS
 
         std::tuple<int, double> resolve_lambda(double lambda) const;
 
+        SireCAS::Expression _getEquation(int stage, const QString &force, const QString &lever) const;
+
+        /** Additional schedules for extra molecules, i.e. that
+         *  run in parallel alongside the default schedule
+         */
+        QHash<qint32, LambdaSchedule> mol_schedules;
+
         /** The set of all constants used across all stages */
         SireCAS::Values constant_values;
 
+        /** The names of all of the forces */
+        QStringList force_names;
+
         /** The names of all of the levers provided by the forcefields */
         QStringList lever_names;
 

corelib/src/libs/SireIO/amberprm.cpp

@@ -562,7 +562,12 @@ void AmberPrm::rebuildLJParameters()
 
                 bool is_exception = false;
 
-                if (std::abs(lj_ij.epsilon().value() - expect.epsilon().value()) <= 1e-6)
+                if (std::abs(epsilon) < 1e-6 and std::abs(expect.epsilon().value()) < 1e-6)
+                {
+                    // this is a LJ pair that involves a ghost or dummy atom
+                    // It should not impact exceptions or combining rules
+                }
+                else if (std::abs(lj_ij.epsilon().value() - expect.epsilon().value()) <= 1e-6)
                 {
                     if (std::abs(lj_ij.sigma().value() - expect.sigma().value()) > 1e-6)
                     {

corelib/src/libs/SireMol/moleculedata.cpp

@@ -363,7 +363,8 @@ bool MoleculeData::operator!=(const MoleculeData &other) const
 
 /** Return a new MoleculeData that contains only the passed selected
     atoms. This allows parts of the molecule to be pulled out and used independently */
-MoleculeData MoleculeData::extract(const AtomSelection &selected_atoms) const
+MoleculeData MoleculeData::extract(const AtomSelection &selected_atoms,
+                                   bool to_same_molecule) const
 {
     selected_atoms.assertCompatibleWith(*this);
 
@@ -433,7 +434,8 @@ MoleculeData MoleculeData::extract(const AtomSelection &selected_atoms) const
 
     // renumber the molecule to remove confusion as to why
     // extracted molecules cannot be combined
-    editor.renumber();
+    if (not to_same_molecule)
+        editor.renumber();
 
     return editor.commit().data();
 }

corelib/src/libs/SireMol/moleculedata.h

@@ -188,7 +188,8 @@ namespace SireMol
             return props;
         }
 
-        MoleculeData extract(const AtomSelection &selected_atoms) const;
+        MoleculeData extract(const AtomSelection &selected_atoms,
+                             bool to_same_molecule = false) const;
 
         QStringList propertyKeys() const;
 

corelib/src/libs/SireMol/moleculeview.cpp

@@ -476,20 +476,6 @@ void MoleculeView::update(const MoleculeData &moldata)
 {
     if (moldata.number() == this->data().number())
     {
-        if (moldata.info().UID() != this->data().info().UID())
-        {
-            throw SireError::incompatible_error(QObject::tr(
-                                                    "Cannot update molecule %1 because the layout for the new "
-                                                    "version of the molecule (%2) is different. This is likely "
-                                                    "because the molecule has been edited and the layout of "
-                                                    "atoms, residues etc has been changed. To update this molecule "
-                                                    "you will need to replace it in the container.")
-                                                    .arg(this->molecule().toString())
-                                                    .arg(Molecule(moldata).toString()),
-                                                CODELOC);
-        }
-
-        // this is the same molecule with the same molecular layout
         d = moldata;
     }
 }
@@ -1934,7 +1920,7 @@ Selector<Segment> MoleculeView::selectAllSegments() const
     selected atoms. This allows the used to pull out parts of a larger molecule,
     e.g. if they want to have only selected residues in a protein and do not
     want to have to store or manipulate the larger protein molecule */
-Molecule MoleculeView::extract() const
+Molecule MoleculeView::extract(bool to_same_molecule) const
 {
     if (d->info().nAtoms() == 0 or this->isEmpty())
         return Molecule();
@@ -1944,7 +1930,7 @@ Molecule MoleculeView::extract() const
 
     else
     {
-        return Molecule(d->extract(this->selection()));
+        return Molecule(d->extract(this->selection(), to_same_molecule));
     }
 }
 

corelib/src/libs/SireMol/moleculeview.h

@@ -278,7 +278,7 @@ namespace SireMol
 
         virtual bool isSelector() const;
 
-        Molecule extract() const;
+        Molecule extract(bool to_same_molecule = false) const;
 
         void update(const MoleculeData &moldata);
         void update(const MoleculeView &molview);

corelib/src/libs/SireMol/selectorm.hpp

@@ -65,7 +65,7 @@ namespace SireMol
     {
 
         friend SIREMOL_EXPORT QDataStream & ::operator<< <>(QDataStream &, const SelectorM<T> &);
-        friend SIREMOL_EXPORT QDataStream & ::operator>><>(QDataStream &, SelectorM<T> &);
+        friend SIREMOL_EXPORT QDataStream & ::operator>> <>(QDataStream &, SelectorM<T> &);
 
     public:
         typedef typename QList<Selector<T>>::const_iterator iterator;
@@ -93,7 +93,7 @@ namespace SireMol
         template <class U>
         SelectorM(const SelectorM<U> &other, const QList<qint64> &idxs);
         template <class U>
-        SelectorM(const SelectorM<U> &other, const QString &name);
+        SelectorM(const SelectorM<U> &other, const QString &name, const SireBase::PropertyMap &map = SireBase::PropertyMap());
         template <class U>
         SelectorM(const SelectorM<U> &other, const typename T::ID &id);
 
@@ -1072,7 +1072,9 @@ namespace SireMol
 
     template <class T>
     template <class U>
-    SIRE_OUTOFLINE_TEMPLATE SelectorM<T>::SelectorM(const SelectorM<U> &other, const QString &name)
+    SIRE_OUTOFLINE_TEMPLATE SelectorM<T>::SelectorM(const SelectorM<U> &other,
+                                                    const QString &name,
+                                                    const SireBase::PropertyMap &map)
         : SireBase::ConcreteProperty<SelectorM<T>, SireBase::Property>()
     {
         for (const auto &view : other)
@@ -1093,7 +1095,7 @@ namespace SireMol
             // try a search
             try
             {
-                this->operator=(SelectorM<T>(other.search(name)));
+                this->operator=(SelectorM<T>(other.search(name, map)));
             }
             catch (...)
             {