Rearrange nonbond evaluation

Separate out 1-4 interactions into a separate _Terms14 vector.
Add code to rebuild the pair list in _Terms

Author: Christian Schafmeisterr

include/cando/chem/energyFunction.h

@@ -263,16 +263,16 @@ namespace chem {
     EnergyAtom*     getEnergyAtomPointer(Atom_sp a);
 
     void assignAtomTypes(Matter_sp matter, bool show_progress);
-    void defineForMatter(Matter_sp agg, bool useExcludedAtoms, core::T_sp keepInteractionFactory=nil<core::T_O>(), bool assign_types=true );
-    void defineForMatterWithAtomTypes(Matter_sp matter, bool useExcludedAtoms, core::T_sp keepInteractionFactory, core::T_sp cip_priorities, core::HashTable_sp atomTypes );
+    void defineForMatter(Matter_sp agg, bool useExcludedAtoms, core::T_sp keepInteractionFactory=nil<core::T_O>(), bool assign_types=true, core::T_sp coordinates=nil<core::T_O>() );
+    void defineForMatterWithAtomTypes(Matter_sp matter, bool useExcludedAtoms, core::T_sp keepInteractionFactory, core::T_sp cip_priorities, core::HashTable_sp atomTypes, core::T_sp coordinates );
     void generateStandardEnergyFunctionTables(Matter_sp mol,
                                               FFStretchDb_sp stretchDb,
                                               FFAngleDb_sp angleDb,
                                               FFPtorDb_sp ptorDb,
                                               FFItorDb_sp itorDb,
                                               core::T_sp keepInteractionFactory,
                                               core::HashTable_sp atomTypes);
-    void generateNonbondEnergyFunctionTables(bool useExcludedAtoms, Matter_sp agg, core::T_sp forceField, core::T_sp keepInteractionFactory, core::HashTable_sp atomTypes );
+    void generateNonbondEnergyFunctionTables(bool useExcludedAtoms, Matter_sp agg, core::T_sp forceField, core::T_sp keepInteractionFactory, core::HashTable_sp atomTypes, core::T_sp coordinates );
     void generateRestraintEnergyFunctionTables(Matter_sp agg, core::T_sp nonbonds, core::T_sp keepInteractionFactory, core::T_sp cip_priorities, core::HashTable_sp atomTypes );
 
     /*! Add the restraints to the energy function.

include/cando/chem/energyNonbond.h

@@ -74,7 +74,6 @@ class EnergyNonbond : public EnergyTerm {
   Atom_sp      	_Atom1_enb;
   Atom_sp      	_Atom2_enb;
   TermNonBond	term;
-  bool		_Is14;
 #if TURN_ENERGY_FUNCTION_DEBUG_ON
   bool		_calcForce;
   bool		_calcDiagonalHessian;
@@ -171,13 +170,17 @@ class EnergyNonbond_O : public EnergyComponent_O
   double              _Nonbond_invdd; // Distance dependent dielectric reciprocal
     size_t                              _NonbondsKept;
     size_t                              _NonbondsDiscarded;
-    // Original way of defining nonbonds with list of nonbond terms
+  gctools::Vec0<TermType>	_Terms14;
+  // Pairlist
+  AtomTable_sp                  _AtomTable;
+  core::T_sp                    _NonbondForceField;
+  core::HashTable_sp            _AtomTypes;
+  core::T_sp                    _KeepInteractionFactory;
   gctools::Vec0<TermType>	_Terms;
   gctools::Vec0<TermType>	_BeyondThresholdTerms;
     // Correct way of defining nonbonds using excluded atom indexes
   // FIXME:  Get rid of _FFNonbondDb - and the old code that uses it.
   core::T_sp        _FFNonbondDb;
-  AtomTable_sp          _AtomTable;
   //  Tables for nonbonded calculation using excluded atoms and the Amber way
   size_t                     _ntypes;          // ntypes
   core::SimpleVector_sp      _atom_name_vector;  // atom-name-vector
@@ -200,12 +203,12 @@ class EnergyNonbond_O : public EnergyComponent_O
 
  public:
   typedef gctools::Vec0<TermType>::iterator iterator;
-  iterator begin() { return this->_Terms.begin(); };
-  iterator end() { return this->_Terms.end(); };
+  //  iterator begin() { return this->_Terms.begin(); };
+  //  iterator end() { return this->_Terms.end(); };
 //added by G 7.19.2011
  public:
 //  core::List_sp termAtIndex(size_t index) const;
-  virtual size_t numberOfTerms() { return this->_Terms.size();};
+  virtual size_t numberOfTerms() { return this->_Terms.size(); };
   void callForEachTerm(core::Function_sp callback);
  public:
 
@@ -216,6 +219,7 @@ class EnergyNonbond_O : public EnergyComponent_O
   CL_DEFMETHOD core::SimpleVector_int32_t_sp excluded_atom_list() const { return this->_ExcludedAtomIndexes;}
  public:
   void constructNonbondTermsBetweenResidues(Residue_sp res1, Residue_sp res2, core::T_sp nbForceField, core::HashTable_sp atomTypes );
+  void addTerm14(const TermType& term);
   void addTerm(const TermType& term);
   virtual void dumpTerms(core::HashTable_sp atomTypes);
 
@@ -269,7 +273,8 @@ class EnergyNonbond_O : public EnergyComponent_O
 
   core::T_sp getFFNonbondDb();
 
-  void constructNonbondTermsFromAtomTable(bool ignore14s, AtomTable_sp atomTable, core::T_sp nbforceField, core::HashTable_sp atomTypes, core::T_sp keepInteractionFactory );
+  void updatePairList(core::T_sp tcoordinates);
+  void constructNonbondTermsFromAtomTable( AtomTable_sp atomTable, core::T_sp nbforceField, core::HashTable_sp atomTypes, core::T_sp keepInteractionFactory, core::T_sp coordinates);
   void constructNonbondTermsBetweenMatters( Matter_sp matter1, Matter_sp matter2, EnergyFunction_sp energyFunction, core::T_sp keepInteractionFactory );
   void construct14InteractionTerms(AtomTable_sp atomTable, Matter_sp matter, core::T_sp nbforceField, core::T_sp keepInteractionFactory, core::HashTable_sp atomTypes );
   void constructExcludedAtomListFromAtomTable(AtomTable_sp atomTable, core::T_sp nbforceField, core::T_sp keepInteractionFactory );

src/chem/energyFunction.cc

@@ -1292,103 +1292,92 @@ __END_DOC
   SYMBOL_EXPORT_SC_(ChemPkg,STARcurrent_aromaticity_informationSTAR);
 
   CL_LISPIFY_NAME("defineForMatter");
-  CL_LAMBDA((energy-function chem:energy-function) matter &key use-excluded-atoms (keep-interaction-factory t) (assign-types t));
-  CL_DEFMETHOD void EnergyFunction_O::defineForMatter(Matter_sp matter, bool useExcludedAtoms, core::T_sp keepInteractionFactory, bool assign_types )
-  {
-    if ( !(matter.isA<Aggregate_O>() || matter.isA<Molecule_O>() ) )
-      {
-        SIMPLE_ERROR("You can only define energy functions for Aggregates or Molecules");
-      }
-
-    core::DynamicScopeManager scope(_sym_STARparameter_warningsSTAR,nil<core::T_O>());
-    //
-    // Identify rings
-    //
-    if (chem__verbose(0)) core::clasp_write_string("Searching for rings.\n");
-
-    // If *current-rings* is bound then reuse it, otherwise, calculate rings and bind those
-    core::T_sp rings = unbound<core::T_sp>();
-    if (_sym_STARcurrent_ringsSTAR->boundP()) {
-      rings = _sym_STARcurrent_ringsSTAR->symbolValue();
-    } else if (keepInteractionFactory.notnilp()) {
-      rings = RingFinder_O::identifyRings(matter);
-    } else {
-      rings = nil<core::T_O>();
+CL_LAMBDA((energy-function chem:energy-function) matter &key use-excluded-atoms (keep-interaction-factory t) (assign-types t) (coordinates nil));
+CL_DEFMETHOD void EnergyFunction_O::defineForMatter(Matter_sp matter, bool useExcludedAtoms, core::T_sp keepInteractionFactory, bool assign_types, core::T_sp coordinates )
+{
+  if ( !(matter.isA<Aggregate_O>() || matter.isA<Molecule_O>() ) )
+    {
+      SIMPLE_ERROR("You can only define energy functions for Aggregates or Molecules");
     }
-    core::DynamicScopeManager ring_scope(_sym_STARcurrent_ringsSTAR, rings );
 
-    //
-    // Assign relative Cahn-Ingold-Preylog priorities
-    //
-    if (chem__verbose(0)) core::clasp_write_string("Assigning CIP priorities.\n");
-    core::HashTable_sp cip = unbound<core::HashTable_O>();
-    if (keepInteractionFactory.notnilp()) {
-      cip = CipPrioritizer_O::assignPrioritiesHashTable(matter);
-    }
+  core::DynamicScopeManager scope(_sym_STARparameter_warningsSTAR,nil<core::T_O>());
+  //
+  // Identify rings
+  //
+  if (chem__verbose(0)) core::clasp_write_string("Searching for rings.\n");
+
+  // If *current-rings* is bound then reuse it, otherwise, calculate rings and bind those
+  core::T_sp rings = unbound<core::T_sp>();
+  if (_sym_STARcurrent_ringsSTAR->boundP()) {
+    rings = _sym_STARcurrent_ringsSTAR->symbolValue();
+  } else if (keepInteractionFactory.notnilp()) {
+    rings = RingFinder_O::identifyRings(matter);
+  } else {
+    rings = nil<core::T_O>();
+  }
+  core::DynamicScopeManager ring_scope(_sym_STARcurrent_ringsSTAR, rings );
 
-    //
-    // Assign atom types for each molecule
-    //
+  //
+  // Assign relative Cahn-Ingold-Preylog priorities
+  //
+  if (chem__verbose(0)) core::clasp_write_string("Assigning CIP priorities.\n");
+  core::HashTable_sp cip = unbound<core::HashTable_O>();
+  if (keepInteractionFactory.notnilp()) {
+    cip = CipPrioritizer_O::assignPrioritiesHashTable(matter);
+  }
 
-    core::HashTable_sp force_fields = core::HashTable_O::createEq();
-    Loop moleculeLoop;
+  //
+  // Assign atom types for each molecule
+  //
+
+  core::HashTable_sp force_fields = core::HashTable_O::createEq();
+  Loop moleculeLoop;
+  moleculeLoop.loopTopGoal(matter,MOLECULES);
+  while (moleculeLoop.advanceLoopAndProcess() ) {
+    Molecule_sp molecule = moleculeLoop.getMolecule();
+    core::T_sp force_field_name = molecule->force_field_name();
+    if (force_fields->gethash(force_field_name).nilp()) {
+      core::T_sp combined_force_field = core::eval::funcall(_sym_find_force_field,force_field_name);
+      force_fields->setf_gethash(force_field_name,combined_force_field);
+    }
+  }
+  core::HashTable_sp atomTypes = core::HashTable_O::createEq();
+  // Assign given atom-types
+  Loop atom_loop;
+  atom_loop.loopTopGoal(matter,ATOMS);
+  while (atom_loop.advanceLoopAndProcess()) {
+    Atom_sp atom = atom_loop.getAtom();
+    core::T_sp type = atom->getPropertyOrDefault( kw::_sym_given_atom_type, nil<core::T_O>() );
+    if (type.notnilp()) {
+      atomTypes->setf_gethash(atom,type);
+    }
+  }
+  this->_AtomTypes = atomTypes;
+  if (assign_types) {
+    if (chem__verbose(0)) core::clasp_write_string("Assigning atom types.\n");
     moleculeLoop.loopTopGoal(matter,MOLECULES);
     while (moleculeLoop.advanceLoopAndProcess() ) {
       Molecule_sp molecule = moleculeLoop.getMolecule();
       core::T_sp force_field_name = molecule->force_field_name();
-      if (force_fields->gethash(force_field_name).nilp()) {
-        core::T_sp combined_force_field = core::eval::funcall(_sym_find_force_field,force_field_name);
-        force_fields->setf_gethash(force_field_name,combined_force_field);
-      }
-    }
-    core::HashTable_sp atomTypes = core::HashTable_O::createEq();
-    // Assign given atom-types
-    Loop atom_loop;
-    atom_loop.loopTopGoal(matter,ATOMS);
-    while (atom_loop.advanceLoopAndProcess()) {
-      Atom_sp atom = atom_loop.getAtom();
-      core::T_sp type = atom->getPropertyOrDefault( kw::_sym_given_atom_type, nil<core::T_O>() );
-      if (type.notnilp()) {
-        atomTypes->setf_gethash(atom,type);
-      }
-    }
-    this->_AtomTypes = atomTypes;
-    if (assign_types) {
-      if (chem__verbose(0)) core::clasp_write_string("Assigning atom types.\n");
-      moleculeLoop.loopTopGoal(matter,MOLECULES);
-      while (moleculeLoop.advanceLoopAndProcess() ) {
-        Molecule_sp molecule = moleculeLoop.getMolecule();
-        core::T_sp force_field_name = molecule->force_field_name();
-        [[maybe_unused]]core::T_sp use_given_types = molecule->force_field_use_given_types();
-        //      if (use_given_types.nilp()) {
-        core::T_sp combined_force_field = force_fields->gethash(force_field_name);
-        if (chem__verbose(0)) core::clasp_write_string(fmt::format("Assigning atom types for molecule {} using {}.\n" , _rep_(molecule->getName()) , _rep_(force_field_name)));
-        //
-        // Calculate aromaticity using the rings we just calculated
-        //
-        core::T_sp aromaticity_info = core::eval::funcall(_sym_identify_aromatic_rings,matter,force_field_name);
-        if (aromaticity_info.nilp()) SIMPLE_ERROR("The aromaticity-info was NIL when about to assign force field types - it should not be");
-        core::DynamicScopeManager aromaticity_scope(_sym_STARcurrent_aromaticity_informationSTAR,aromaticity_info);
-        core::eval::funcall(_sym_assign_force_field_types,combined_force_field,molecule,atomTypes);
-#if 0
-      } else {
-        if (chem__verbose(0)) core::clasp_write_string(fmt::format("Assigning atom types for molecule {} using given-types for {}.\n" , _rep_(molecule->getName()) , _rep_(force_field_name)));
-      Loop atom_loop;
-      atom_loop.loopTopGoal(molecule,ATOMS);
-      while (atom_loop.advanceLoopAndProcess()) {
-        Atom_sp atom = atom_loop.getAtom();
-        atomTypes->setf_gethash(atom,atom->atomType());
-      }
+      [[maybe_unused]]core::T_sp use_given_types = molecule->force_field_use_given_types();
+      //      if (use_given_types.nilp()) {
+      core::T_sp combined_force_field = force_fields->gethash(force_field_name);
+      if (chem__verbose(0)) core::clasp_write_string(fmt::format("Assigning atom types for molecule {} using {}.\n" , _rep_(molecule->getName()) , _rep_(force_field_name)));
+      //
+      // Calculate aromaticity using the rings we just calculated
+      //
+      core::T_sp aromaticity_info = core::eval::funcall(_sym_identify_aromatic_rings,matter,force_field_name);
+      if (aromaticity_info.nilp()) SIMPLE_ERROR("The aromaticity-info was NIL when about to assign force field types - it should not be");
+      core::DynamicScopeManager aromaticity_scope(_sym_STARcurrent_aromaticity_informationSTAR,aromaticity_info);
+      core::eval::funcall(_sym_assign_force_field_types,combined_force_field,molecule,atomTypes);
     }
-#endif
   }
-}
-  this->defineForMatterWithAtomTypes(matter,useExcludedAtoms,keepInteractionFactory,cip,atomTypes);
+  this->defineForMatterWithAtomTypes(matter,useExcludedAtoms,keepInteractionFactory,cip,atomTypes,coordinates);
 }
 
 
-CL_LAMBDA((energy-function chem:energy-function) matter &key use-excluded-atoms (keep-interaction-factory t) cip-priorities atom-types);
-CL_DEFMETHOD void EnergyFunction_O::defineForMatterWithAtomTypes(Matter_sp matter, bool useExcludedAtoms, core::T_sp keepInteractionFactory, core::T_sp cip_priorities, core::HashTable_sp atomTypes )
+CL_LAMBDA((energy-function chem:energy-function) matter &key use-excluded-atoms (keep-interaction-factory t) cip-priorities atom-types coordinates);
+CL_DEFMETHOD void EnergyFunction_O::defineForMatterWithAtomTypes(Matter_sp matter, bool useExcludedAtoms, core::T_sp keepInteractionFactory, core::T_sp cip_priorities, core::HashTable_sp atomTypes, core::T_sp coordinates )
 {
   if (keepInteractionFactory.notnilp() && !gc::IsA<core::HashTable_sp>(cip_priorities)) {
     SIMPLE_ERROR("You need to provide a hash-table of atoms to relative CIP priorities - see CipPrioritizer_O::assignPrioritiesHashTable(matter)");
@@ -1507,7 +1496,7 @@ CL_DEFMETHOD void EnergyFunction_O::defineForMatterWithAtomTypes(Matter_sp matte
   if (keepInteractionFactory.notnilp()) {
     if (chem__verbose(1)) core::clasp_write_string("About to calculate nonbond and restraint terms");
     core::T_sp nonbondForceField = this->_AtomTable->nonbondForceFieldForAggregate();
-    this->generateNonbondEnergyFunctionTables(useExcludedAtoms,matter,nonbondForceField,keepInteractionFactory,atomTypes);
+    this->generateNonbondEnergyFunctionTables(useExcludedAtoms,matter,nonbondForceField,keepInteractionFactory,atomTypes,coordinates);
     this->generateRestraintEnergyFunctionTables(matter,nonbondForceField,keepInteractionFactory,cip_priorities,atomTypes);
   }
   core::eval::funcall(_sym_report_parameter_warnings);
@@ -1775,7 +1764,7 @@ CL_DEFMETHOD void EnergyFunction_O::generateStandardEnergyFunctionTables(Matter_
 SYMBOL_EXPORT_SC_(ChemPkg,prepare_amber_energy_nonbond);
 
 CL_DOCSTRING(R"dx(Generate the nonbond energy function tables. The atom types, and CIP priorities need to be precalculated.)dx");
-CL_DEFMETHOD void EnergyFunction_O::generateNonbondEnergyFunctionTables(bool useExcludedAtoms, Matter_sp matter, core::T_sp nonbondForceField, core::T_sp keepInteractionFactory, core::HashTable_sp atomTypes )
+CL_DEFMETHOD void EnergyFunction_O::generateNonbondEnergyFunctionTables(bool useExcludedAtoms, Matter_sp matter, core::T_sp nonbondForceField, core::T_sp keepInteractionFactory, core::HashTable_sp atomTypes, core::T_sp coordinates )
 {
   if (keepInteractionFactory.nilp()) return;
   if (chem__verbose(0))
@@ -1792,7 +1781,9 @@ CL_DEFMETHOD void EnergyFunction_O::generateNonbondEnergyFunctionTables(bool use
     this->_Nonbond->constructExcludedAtomListFromAtomTable(this->_AtomTable, nonbondForceField,keepInteractionFactory);
     this->_Nonbond->construct14InteractionTerms(this->_AtomTable,matter,nonbondForceField,keepInteractionFactory,atomTypes);
   } else {
-    this->_Nonbond->constructNonbondTermsFromAtomTable(false,this->_AtomTable, nonbondForceField,atomTypes, keepInteractionFactory );
+    this->_Nonbond->constructNonbondTermsFromAtomTable(this->_AtomTable, nonbondForceField,atomTypes,
+                                                       keepInteractionFactory, coordinates );
+    this->_Nonbond->construct14InteractionTerms(this->_AtomTable,matter,nonbondForceField,keepInteractionFactory,atomTypes);
   }
   if (chem__verbose(0)) core::clasp_write_string(fmt::format("Built nonbond table for {} terms\n" , this->_Nonbond->numberOfTerms()));
 }