diff --git a/doc/tests/scripts/sidechain_steps.py b/doc/tests/scripts/sidechain_steps.py
index 854ff709f9b88187d7db704bc681a2132740a6df..d73a626dd775fd254e62cd1e085f355a278f7d01 100644
--- a/doc/tests/scripts/sidechain_steps.py
+++ b/doc/tests/scripts/sidechain_steps.py
@@ -7,7 +7,7 @@ prot = io.LoadPDB('data/1CRN.pdb')
 library = sidechain.LoadBBDepLib()
 # we need a rotamer constructor to create any rotamers or 
 # frame residues 
-rot_constructor = sidechain.SCWRLRotamerConstructor(False)
+rot_constructor = sidechain.SCWRL4RotamerConstructor(False)
 
 # create new entity from protein only containing the amino acids
 prot = mol.CreateEntityFromView(prot.Select("peptide=true"), True)
diff --git a/modelling/pymod/_reconstruct_sidechains.py b/modelling/pymod/_reconstruct_sidechains.py
index d9a7dbd5541b9e7ca3102bf9e8cc78bda4b8f1e6..969df78d40bcda96d881abf800892d8de6a0fd3a 100644
--- a/modelling/pymod/_reconstruct_sidechains.py
+++ b/modelling/pymod/_reconstruct_sidechains.py
@@ -339,17 +339,26 @@ def _GetDisulfidBridges(frame_residues, keep_sidechains, cystein_indices,
         if keep_sidechains:
             sg = r.FindAtom("SG")
             if sg.IsValid():
-                particle = sidechain.Particle(sidechain.SidechainParticle.SParticle, 
-                                              sg.GetPos(),
-                                              -0.19, "SG")   
-                particle_list = [particle]
-                # The temperature and self internal_e_prefactor parameter have 
-                # been copied from the SCWRLRotamerConstructor. Hacky...
-                frm_rotamer = sidechain.FRMRotamer(particle_list, 
-                                                   1.69, 1.0, 4.07)
-                frm_rotamer.AddSubrotamerDefinition([0])               
-                frm_rotamer.SetInternalEnergy(0.0)
-                rot_group = sidechain.FRMRotamerGroup([frm_rotamer], i)
+                # we're constructing the required particle through a frame 
+                # residue in the rotamer constructor
+                cys_frame_res = \
+                rotamer_constructor.ConstructSidechainFrameResidue(r.handle, 
+                                                                   sidechain.CYD, 
+                                                                   0)
+                for j in range(len(cys_frame_res)):
+                    if cys_frame_res[j].GetName() == "SG":
+                        particle_list = [cys_frame_res[j]]
+                        # The temperature and self internal_e_prefactor 
+                        # parameter have been copied from the 
+                        # SCWRL4RotamerConstructor. Hacky...
+                        frm_rotamer = sidechain.FRMRotamer(particle_list, 
+                                                           1.69, 1.0, 4.07)
+                        frm_rotamer.AddSubrotamerDefinition([0])               
+                        frm_rotamer.SetInternalEnergy(0.0)
+                        rot_group = sidechain.FRMRotamerGroup([frm_rotamer], i)
+                        break
+                if rot_group == None:
+                    raise RuntimeError("Could not find SG atom in CYS frame res")
 
         if rot_group == None:
             rot_group = _GetRotamerGroup(r.handle, sidechain.CYD, i, 
@@ -447,7 +456,7 @@ def ReconstructSidechains(ent, keep_sidechains=False, build_disulfids=True,
     if type(rotamer_library) is sidechain.BBDepRotamerLib:
         bbdep = True
 
-    rotamer_constructor = sidechain.SCWRLRotamerConstructor(False)
+    rotamer_constructor = sidechain.SCWRL4RotamerConstructor(False)
     
     # take out ligand chain and any non-peptides
     prot = ent.Select("peptide=true and cname!='_'")
diff --git a/modelling/src/sidechain_env_listener.hh b/modelling/src/sidechain_env_listener.hh
index 309930e65e0b0db0469a36c1338e7c39db22a118..cc8c2709a42e01b56f41cabd54a994460acd5de5 100644
--- a/modelling/src/sidechain_env_listener.hh
+++ b/modelling/src/sidechain_env_listener.hh
@@ -24,7 +24,7 @@
 #include <promod3/sidechain/rotamer_id.hh>
 #include <promod3/sidechain/rotamer_lib.hh>
 #include <promod3/sidechain/bb_dep_rotamer_lib.hh>
-#include <promod3/sidechain/scwrl_rotamer_constructor.hh>
+#include <promod3/sidechain/scwrl4_rotamer_constructor.hh>
 
 namespace promod3 { namespace modelling {
 
@@ -182,7 +182,7 @@ private:
   bool add_cyd_rotamers_;  // build extra rotamers with r_id = CYD for cysteins
   sidechain::RotamerLibPtr library_;
   sidechain::BBDepRotamerLibPtr bbdep_library_;
-  sidechain::SCWRLRotamerConstructor rot_constructor_;
+  sidechain::SCWRL4RotamerConstructor rot_constructor_;
   ScCBetaSpatialOrganizer env_;
 
   // fixed data per residue (same numbering as base_env.GetAllPosData())
diff --git a/modelling/src/sidechain_reconstructor.cc b/modelling/src/sidechain_reconstructor.cc
index 3e63c3221e189cb5cca392298fcc34545581e971..008c5a41694820a3047ec325fc58831c03646fcf 100644
--- a/modelling/src/sidechain_reconstructor.cc
+++ b/modelling/src/sidechain_reconstructor.cc
@@ -275,14 +275,21 @@ void SidechainReconstructor::BuildDisulfids_(
           // We want to keep existing sidechains (keep_sidechains_ == true) 
           // and everything seems to be present since the sidechain frame 
           // residue exists. We therefore replace the extracted FRMRotamerGroup 
-          // (p) with position data from the frame. 
+          // (p) with the SG position from the frame. 
           // The reason for still creating a FRMRotamerGroup is that it can 
           // enter the ResolveCystein evaluation.
           std::vector<sidechain::Particle> particles;
-          particles.push_back(sidechain::Particle((*f_res)[0].GetParticleType(),
-                                                  (*f_res)[0].GetPos(),
-                                                  (*f_res)[0].GetCharge(),
-                                                  (*f_res)[0].GetName()));
+          for(uint j = 0; j < f_res->size(); ++j) {
+            if((*f_res)[j].GetName() == "SG") {
+              particles.push_back((*f_res)[j]);
+              break;
+            }
+          }
+          if(particles.empty()) {
+            throw promod3::Error("Could not find SG atom in CYS frame residue");
+          }
+          // The temperature and self internal_e_prefactor parameter have been 
+          // copied from the SCWRL4RotamerConstructor. Hacky...
           RotamerPtr r(new Rotamer(particles, (*p)[0]->GetTemperature(), 1.0,
                                    (*p)[0]->GetInternalEnergyPrefactor()));
           std::vector<int> subrotamer_definition(1,0);
diff --git a/sidechain/doc/rotamer.rst b/sidechain/doc/rotamer.rst
index 30b3c7ca3b06602a1fbade371a3abdbc5a4e9029..d3ae44ebe66b3f2139cd7a5ba04e29de80efc9ac 100644
--- a/sidechain/doc/rotamer.rst
+++ b/sidechain/doc/rotamer.rst
@@ -28,104 +28,111 @@ Pairwise interactions between particles give raise to pairwise energies between
 rotamers. Nevertheless, the energy calculation itself happens on the level
 of RotamerGroups and is mostly hidden away in the construction of the
 the :class:`RotamerGraph`. If you're too lazy to build up your rotamers
-by hand, you might be interested in the :class:`SCWRLRotamerConstructor`.
+by hand, you might be interested in the :class:`RotamerConstructor`.
 
 
 
 The Smallest Building Block - The Particle
 --------------------------------------------------------------------------------
 
-Particles give raise to more complex objects such as rotamers and frame residues.
-They are the basis for calculating interaction energies based on a Lennard-Jones
-like term and an approximation of an hbond energy term.
-The Lennard-Jones term gets mainly controlled by the type of the particle, that
-has to be defined at initialization.
+Particles give raise to more complex objects such as rotamers and frame 
+residues. They contain all data required to calculate pairwise energies.
+For every energy function available in ProMod3, there's a particle creation
+function.
 
-.. class:: SidechainParticle
+.. class:: PScoringFunction
 
-  Enumerates the types of particle. Possible values:
+  The available scoring functions between :class:`Particle` objects
 
-  * HParticle   - represents hydrogen
-  * CParticle   - default representation of a carbon
-  * CH1Particle - represents carbon bound to 1 hydrogen
-  * CH2Particle - represents carbon bound to 2 hydrogen
-  * CH3Particle - represents carbon bound to 3 hydrogen
-  * NParticle   - represents nitrogen
-  * OParticle   - default representation of oxygen
-  * OCParticle  - represents carbonyl-oxygen for ASP/GLU
-  * SParticle   - represents sulfur
+  * SCWRL4
 
-.. class:: Particle(type, pos, charge, name)
+.. class:: Particle
 
-  :param type:          Type of particle
-  :param pos:           Positions of particle
-  :param charge:        Charge of particle, will be used in case H-Bonds
-  :param name:          Name of particle. This name will be given to an actual
-                        :class:`ost.mol.AtomHandle` if a rotamer gets applied
-                        to a :class:`ost.mol.ResidueHandle`
+  The particle class. There's no constructor. You can either use the
+  :class:`RotamerConstructor` to create full :class:`RotamerGroup` objects
+  with all underlying particles or the energy function specific creation 
+  functions.
 
-  :type type:           :class:`SidechainParticle`
-  :type pos:            :class:`ost.geom.Vec3`
-  :type charge:         :class:`float`
-  :type name:           :class:`str`
-
-  .. method:: PairwiseEnergy(other)
-
-    Calculates the interaction energy between two particles.
-
-    :param other:       Interaction partner
+  .. method:: PairwiseScore(other_particle)
 
-    :type other:        :class:`Particle`
-
-    :returns:           the interaction energy
-
-
-  .. method:: GetPos()
+    Calculates score between the two particles
 
-    :returns:           The Particle position
+    :param other_particle:  The interacting particle
+    :type other_particle:  :class:`Particle`
 
+    :returns:           The score
+    :rtype:             :class:`float`
 
-  .. method:: GetCharge()
+    :raises:  :exc:`~exceptions.RuntimeError` if the scoring function 
+              parametrization of the two particles is inconsistent         
 
-    :returns:           Charge of particle used for HBond term
-
-  .. method:: GetParticleType()
+  .. method:: GetName()
 
-    :returns:           Internal type of particle
+    :returns:           The name of the particle, which corresponds to the
+                        atom name
+    :rtype:             :class:`str`
 
+  .. method:: GetCollisionDistance()
 
-  .. method:: GetName()
+    :returns:           The "collision distance" all pairs of particles with 
+                        their distance below the sum of their collision 
+                        distances are considered as interacting and thus
+                        evaluated by the underlying scoring function.
 
-    :returns:           Name of particle
+    :rtype:             :class:`float`
 
+  .. method:: GetPos()
 
-  .. method:: IsHBondDonor()
-
-    :returns:           Whether particle has defined polar direction
+    :returns:           The position of the particle
+    :rtype:             :class:`ost.geom.Vec3`
 
+  .. method:: GetScoringFunction()
 
-  .. method:: IsHBondAcceptor()
+    :returns:           The underlying scoring function
+    :rtype:             :class:`PScoringFunction`
 
-    :returns:           Whether particle has defined lone pairs
 
 
-  .. method:: AddLonePair(lone_pair)
 
-    Adds lone pair, that gets evaluated in the HBond term.
+The SCWRL4 scoring function
+--------------------------------------------------------------------------------
 
-    :param lone_pair:   Lone pair direction
+.. class:: SCWRL4ParticleType
 
-    :type lone_pair:    :class:`ost.geom.Vec3`
+  The SCWRL4 energy function differentiates between following particle types
 
+  * HParticle   - represents hydrogen
+  * CParticle   - default representation of a carbon
+  * CH1Particle - represents carbon bound to 1 hydrogen
+  * CH2Particle - represents carbon bound to 2 hydrogen
+  * CH3Particle - represents carbon bound to 3 hydrogen
+  * NParticle   - represents nitrogen
+  * OParticle   - default representation of oxygen
+  * OCParticle  - represents carbonyl-oxygen for ASP/GLU
+  * SParticle   - represents sulfur
 
-  .. method:: SetPolarDirection(polar_direction)
+.. method:: CreateSCWRL4Particle(name, particle_type, pos, [charge, \
+                                 lone_pairs=None, polar_direction=None])
 
-    Sets the polar direction of particle. In case of backbone
-    hydrogen, this would be the direction of backbone-N to backbone-H.
+  Creates and returns a :class:`Particle` that can evaluate the SCWRL4 scoring 
+  function
 
-    :param polar_direction: Polar direction of particle
+  :param name:          The name of the particle
+  :param particle_type: The type of the particle
+  :param pos:           The position of the particle
+  :param charge:        The charge of the particle, relevant for the hydrogen 
+                        bond term
+  :param lone_pairs:    Direction of all possible lone pairs of the particle,
+                        relevant for the hydrogen bond term
+  :param polar_direction: The polar direction of the particle,
+                          relevant for the hdrogen bond term
 
-    :type polar_direction:  :class:`ost.geom.Vec3` 
+  :type name:           :class:`str`
+  :type particle_type:  :class:`SCWRL4ParticleType`
+  :type pos:            :class:`ost.geom.Vec3`
+  :type charge:         :class:`float`
+  :type lone_pairs:     :class:`ost.geom.Vec3List`
+  :type polar_direction: :class:`ost.geom.Vec3`
 
 
 Rotamers
diff --git a/sidechain/doc/rotamer_constructor.rst b/sidechain/doc/rotamer_constructor.rst
index 4b1bf20bd00c0813ac54895f2153a77367a19ea1..b36d7d88619892d046e8ebe09aeaca7f71713046 100644
--- a/sidechain/doc/rotamer_constructor.rst
+++ b/sidechain/doc/rotamer_constructor.rst
@@ -31,7 +31,7 @@ Constructing Rotamers and Frame Residues
 
   Abstract base class that cannot be initialized from Python. It builds 
   an interface implemented by energy function specific constructors 
-  (e.g. :class:`SCWRLRotamerConstructor`). 
+  (e.g. :class:`SCWRL4RotamerConstructor`). 
 
   .. method:: ConstructRRMRotamerGroup(res, id, residue_index, rot_lib,\
                                        [probability_cutoff = 0.98])
@@ -181,7 +181,7 @@ Constructing Rotamers and Frame Residues
 
 
 
-.. class:: SCWRLRotamerConstructor(cb_in_sidechain)
+.. class:: SCWRL4RotamerConstructor(cb_in_sidechain)
 
   This object implements the full interface defined in 
   :class:`RotamerConstructor` and constructs rotamers and frame residues that 
diff --git a/sidechain/pymod/__init__.py b/sidechain/pymod/__init__.py
index 311426d9e13dc020e82e4f64a0afad86803bc11a..286993f2ec39f36072634dab870cf384e0c4e14e 100644
--- a/sidechain/pymod/__init__.py
+++ b/sidechain/pymod/__init__.py
@@ -13,6 +13,5 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-
+from ost import geom
 from _sidechain import *
-
diff --git a/sidechain/pymod/export_particle.cc b/sidechain/pymod/export_particle.cc
index 64eb1a0660bcf678d29ebde03afe41d2e3e93377..7d0d472fe650daacff93f13669fa7476c0ed318f 100644
--- a/sidechain/pymod/export_particle.cc
+++ b/sidechain/pymod/export_particle.cc
@@ -15,24 +15,53 @@
 
 
 #include <boost/python.hpp>
-#include <boost/python/register_ptr_to_python.hpp>
 #include <promod3/sidechain/particle.hh>
-
+#include <promod3/sidechain/scwrl4_particle_scoring.hh>
 
 using namespace boost::python;
 using namespace promod3::sidechain;
 
-
 namespace{
-  Real WrapPairwiseEnergy(ParticlePtr p, ParticlePtr other){
-    return p->PairwiseEnergy(other);
+
+ParticlePtr WrapCreateSCWRL4Particle(const String& name, 
+                                     SCWRL4ParticleType particle_type,
+                                     const geom::Vec3& pos,
+                                     Real charge,
+                                     const geom::Vec3List& lone_pairs,
+                                     const geom::Vec3& polar_direction) {
+  SCWRL4Param* p = new SCWRL4Param(particle_type, pos, charge);
+  if(geom::Length2(polar_direction) > Real(0.000001)) {
+    p->SetPolarDirection(polar_direction);
+  }
+  for(uint i = 0; i < lone_pairs.size(); ++i) {
+    p->AddLonePair(lone_pairs[i]);
   }
+  ParticlePtr return_ptr(new Particle(name, p));
+  return return_ptr;
+}
 
 }
 
-void export_Particle()
-{
-  enum_<SidechainParticle>("SidechainParticle")
+
+void export_Particle() {
+  
+  enum_<PScoringFunction>("PScoringFunction")
+    .value("SCWRL4", SCWRL4)
+  ;
+
+  class_<Particle, ParticlePtr>("Particle", no_init)
+    .def("PairwiseScore", &Particle::PairwiseScore, (arg("other_particle")))
+    .def("GetName", &Particle::GetName, 
+         return_value_policy<copy_const_reference>())
+    .def("GetPos", &Particle::GetPos, 
+         return_value_policy<copy_const_reference>())
+    .def("GetCollisionDistance", &Particle::GetCollisionDistance)
+    .def("GetScoringFunction", &Particle::GetScoringFunction)
+  ;  
+
+
+  // SCWRL4 specific stuff
+  enum_<SCWRL4ParticleType>("SCWRL4ParticleType")
     .value("HParticle", HParticle)
     .value("CParticle", CParticle)
     .value("CH1Particle", CH1Particle)
@@ -44,20 +73,12 @@ void export_Particle()
     .value("SParticle", SParticle)
   ;
 
-  class_<Particle>
-    ("Particle", init<SidechainParticle,const geom::Vec3&,Real,const String&>())
-    .def("PairwiseEnergy", &WrapPairwiseEnergy, (arg("other")))
-    .def("GetPos", &Particle::GetPos,
-         return_value_policy<copy_const_reference>())
-    .def("GetCharge", &Particle::GetCharge)
-    .def("GetParticleType", &Particle::GetParticleType)
-    .def("GetName", &Particle::GetName,
-         return_value_policy<copy_const_reference>())
-    .def("IsHBondDonor", &Particle::IsHBondDonor)
-    .def("IsHBondAcceptor", &Particle::IsHBondAcceptor)
-    .def("AddLonePair", &Particle::AddLonePair,(arg("lone_pair")))
-    .def("SetPolarDirection", &Particle::SetPolarDirection,(arg("polar_direction")))
-  ;  
+  def("CreateSCWRL4Particle", 
+      &WrapCreateSCWRL4Particle, (arg("name"), 
+                                  arg("particle_type"),
+                                  arg("pos"),  
+                                  arg("charge") = 0.0, 
+                                  arg("lone_pairs") = geom::Vec3List(),
+                                  arg("polar_direction") = geom::Vec3()));
 
-  register_ptr_to_python<ParticlePtr>();
 }
diff --git a/sidechain/pymod/export_rotamer_constructor.cc b/sidechain/pymod/export_rotamer_constructor.cc
index 11829eadee51a756d7dd8b368acf9b1e17ee2894..8b6b6890e0a65cda168fc55e89001f8e95ac6fe6 100644
--- a/sidechain/pymod/export_rotamer_constructor.cc
+++ b/sidechain/pymod/export_rotamer_constructor.cc
@@ -20,7 +20,7 @@
 
 #include <promod3/sidechain/rotamer.hh>
 #include <promod3/sidechain/rotamer_group.hh>
-#include <promod3/sidechain/scwrl_rotamer_constructor.hh>
+#include <promod3/sidechain/scwrl4_rotamer_constructor.hh>
 
 using namespace boost::python;
 using namespace promod3::sidechain;
@@ -29,8 +29,8 @@ using namespace promod3;
 
 namespace{
 
-SCWRLRotamerConstructorPtr WrapSCWRLRotamerConstructorInit(bool cb_in_sidechain) {
-  SCWRLRotamerConstructorPtr p(new SCWRLRotamerConstructor(cb_in_sidechain));
+SCWRL4RotamerConstructorPtr WrapSCWRL4RotamerConstructorInit(bool cb_in_sidechain) {
+  SCWRL4RotamerConstructorPtr p(new SCWRL4RotamerConstructor(cb_in_sidechain));
   return p;
 }
 
@@ -311,38 +311,38 @@ void export_RotamerConstructor(){
                                                                arg("entries"),
                                                                arg("probability_cutoff") = 0.98))
     .def("ConstructBackboneFrameResidue", &WrapBBFrame_res,(arg("residue"),
-                                                                          arg("rotamer_id"),
-                                                                          arg("residue_index"),
-                                                                          arg("phi") = -1.0472,
-                                                                          arg("n_ter")=false,
-                                                                          arg("c_ter")=false))
+                                                            arg("rotamer_id"),
+                                                            arg("residue_index"),
+                                                            arg("phi") = -1.0472,
+                                                            arg("n_ter")=false,
+                                                            arg("c_ter")=false))
     .def("ConstructBackboneFrameResidue", &WrapBBFrame_aa,(arg("all_atom_pos"),
-                                                                         arg("all_atom_pos_idx"),
-                                                                         arg("rotamer_id"),
-                                                                         arg("residue_index"),
-                                                                         arg("phi") = -1.0472,
-                                                                         arg("n_ter")=false,
-                                                                         arg("c_ter")=false))
+                                                           arg("all_atom_pos_idx"),
+                                                           arg("rotamer_id"),
+                                                           arg("residue_index"),
+                                                           arg("phi") = -1.0472,
+                                                           arg("n_ter")=false,
+                                                           arg("c_ter")=false))
     .def("ConstructSidechainFrameResidue", &WrapSCFrame_res,(arg("residue"),
-                                                                           arg("rotamer_id"),
-                                                                           arg("residue_index")))
+                                                             arg("rotamer_id"),
+                                                             arg("residue_index")))
     .def("ConstructSidechainFrameResidue", &WrapSCFrame_aa,(arg("all_atom_pos"),
-                                                                          arg("all_atom_pos_idx"),
-                                                                          arg("rotamer_id"),
-                                                                          arg("residue_index")))
+                                                            arg("all_atom_pos_idx"),
+                                                            arg("rotamer_id"),
+                                                            arg("residue_index")))
     .def("AssignInternalEnergies", &AssignInternalEnergiesRRM)
     .def("AssignInternalEnergies", &AssignInternalEnergiesFRM)
   ;
 
 
-  class_<SCWRLRotamerConstructor, SCWRLRotamerConstructorPtr, 
-         bases<RotamerConstructor> >("SCWRLRotamerConstructor", no_init)
-    .def("__init__", boost::python::make_constructor(&WrapSCWRLRotamerConstructorInit))
+  class_<SCWRL4RotamerConstructor, SCWRL4RotamerConstructorPtr, 
+         bases<RotamerConstructor> >("SCWRL4RotamerConstructor", no_init)
+    .def("__init__", boost::python::make_constructor(&WrapSCWRL4RotamerConstructorInit))
     .def("ConstructFrameResidue",
-         &SCWRLRotamerConstructor::ConstructFrameResidue, (arg("res"),
+         &SCWRL4RotamerConstructor::ConstructFrameResidue, (arg("res"),
                                                            arg("residue_index")))
     .def("ConstructFrameResidueHeuristic",
-         &SCWRLRotamerConstructor::ConstructFrameResidueHeuristic, (arg("res"),
+         &SCWRL4RotamerConstructor::ConstructFrameResidueHeuristic, (arg("res"),
                                                                     arg("residue_index"),
                                                                     arg("comp_lib")))
   ;
diff --git a/sidechain/src/CMakeLists.txt b/sidechain/src/CMakeLists.txt
index aca188a6a4e1198a1583f1b7d01267ca628b6992..ac3b4ed867e19ba8c82c75c02dbc454206686938 100644
--- a/sidechain/src/CMakeLists.txt
+++ b/sidechain/src/CMakeLists.txt
@@ -4,6 +4,9 @@ set(SIDECHAIN_HEADERS
   frame.hh
   rotamer_graph.hh
   particle.hh
+  particle_scoring_base.hh
+  particle_scoring.hh
+  scwrl4_particle_scoring.hh
   rotamer.hh
   rotamer_cruncher.hh
   rotamer_density.hh
@@ -15,7 +18,7 @@ set(SIDECHAIN_HEADERS
   sidechain_connector.hh
   sidechain_object_loader.hh
   rotamer_constructor.hh
-  scwrl_rotamer_constructor.hh
+  scwrl4_rotamer_constructor.hh
   subrotamer_optimizer.hh
 )
 
@@ -23,7 +26,8 @@ set(SIDECHAIN_SOURCES
   bb_dep_rotamer_lib.cc
   disulfid.cc
   frame.cc
-  particle.cc
+  particle_scoring.cc
+  scwrl4_particle_scoring.cc
   rotamer.cc
   rotamer_cruncher.cc
   rotamer_density.cc
@@ -35,7 +39,7 @@ set(SIDECHAIN_SOURCES
   sidechain_connector.cc
   sidechain_object_loader.cc
   rotamer_constructor.cc
-  scwrl_rotamer_constructor.cc
+  scwrl4_rotamer_constructor.cc
   subrotamer_optimizer.cc
 )
 
diff --git a/sidechain/src/frame.cc b/sidechain/src/frame.cc
index 891e2959674a32fe2bfae76ef33fd636e5cddd42..b6f7dbba2119e94bdf154b0053fa66be9656d377 100644
--- a/sidechain/src/frame.cc
+++ b/sidechain/src/frame.cc
@@ -47,7 +47,7 @@ Frame::Frame(const std::vector<FrameResiduePtr>& frame_residues):
       i != frame_residues_.end(); ++i){
     overall_size += (*i)->size();
   }
-  frame_particles_.resize(overall_size);
+  scoring_parameters_.clear();
   residue_indices_.resize(overall_size);
 
   std::vector<Real> x_particle_positions(overall_size);
@@ -64,7 +64,7 @@ Frame::Frame(const std::vector<FrameResiduePtr>& frame_residues):
     actual_index = (*i)->GetResidueIndex();
     for(uint j = 0; j < (*i)->size(); ++j){
       Particle& p = (**i)[j];
-      frame_particles_[particle_counter] = &p;
+      scoring_parameters_.push_back(p.GetSParam());
       residue_indices_[particle_counter] = actual_index;
       pos = p.GetPos();
       x_particle_positions[particle_counter] = pos[0];
diff --git a/sidechain/src/frame.hh b/sidechain/src/frame.hh
index b073bb0989e60eb581a17f33d299de0bdc201368..324680c5653e661eb84dda1c501e406e8982e6da 100644
--- a/sidechain/src/frame.hh
+++ b/sidechain/src/frame.hh
@@ -83,13 +83,13 @@ public:
     return collision_tree_; 
   }
 
-  const std::vector<Particle*>& GetParticles() const {
-    return frame_particles_;
+  const PScoringParamContainer& GetScoringParameters() const {
+    return scoring_parameters_;
   }
 
 private:
 
-  std::vector<Particle*> frame_particles_;
+  PScoringParamContainer scoring_parameters_;
   std::vector<FrameResiduePtr> frame_residues_;
   std::vector<uint> residue_indices_;
   promod3::core::TetrahedralPolytopeTree collision_tree_;
diff --git a/sidechain/src/particle.hh b/sidechain/src/particle.hh
index 9675e0ccf95afc8bd9682641c0e8673e0bcae711..16829fd763eb8666e5dcc8f99d54c887b81dc741 100644
--- a/sidechain/src/particle.hh
+++ b/sidechain/src/particle.hh
@@ -17,31 +17,12 @@
 #ifndef PROMOD3_PARTICLE_HH
 #define PROMOD3_PARTICLE_HH
 
-#include <math.h>
 #include <boost/shared_ptr.hpp>
-#include <ost/geom/vecmat3_op.hh>
-
-namespace{
-
-const Real _max_radii[9] = {1.0300, 2.2400, 2.2400, 2.2400, 
-                            2.2400, 1.7067, 1.7200,
-                            1.7200, 2.2800};
-}
+#include <promod3/sidechain/particle_scoring.hh>
+#include <promod3/core/message.hh>
 
 namespace promod3{ namespace sidechain{
 
-enum SidechainParticle {
-  HParticle,
-  CParticle,
-  CH1Particle,
-  CH2Particle,
-  CH3Particle,
-  NParticle,
-  OParticle,
-  OCParticle,
-  SParticle
-};
-
 class Particle;
 typedef boost::shared_ptr<Particle> ParticlePtr;
 
@@ -49,63 +30,108 @@ class Particle {
 
 public:
 
-  Particle() { }
+  Particle(): name_(""), s_param_(NULL) { }
+
+  Particle(const String& name, PScoringParam* s_param): name_(name),
+                                                        s_param_(s_param) { }
 
-  Particle(SidechainParticle p, const geom::Vec3& pos, Real charge,
-           String name = ""):
-              particle_id_(p), 
-              pos_(pos), charge_(charge), name_(name),
-              is_hbond_donor_(false), is_hbond_acceptor_(false) { }
+  Particle(const Particle& other) {
+    name_ = other.name_;
+    if(other.s_param_ != NULL) {
+      s_param_ = other.s_param_->Copy();
+    } else {
+      s_param_ = NULL;
+    }
+  }
 
-  bool operator==(const Particle& particle) const;   
-  bool operator!=(const Particle& particle) const {     
-    return !this->operator==(particle);   
+  ~Particle() {
+    // particle takes the ownership of the param ptr
+    if(s_param_ != NULL) {
+      delete s_param_;
+    }
   }
 
-  Real GetCollisionDistance() const { return _max_radii[particle_id_]; }
+  Particle& operator= (const Particle& other) {
+    if(this != &other) {      
+      name_ = other.name_;        
+      PScoringParam* s_param_2 = NULL;
+      if(other.s_param_ != NULL) {
+        s_param_2 = other.s_param_->Copy();
+      }    
+      if(s_param_ != NULL) {
+        delete s_param_;      
+      }             
+      s_param_ = s_param_2;
+    }
+    return *this;
+  }
 
-  Real PairwiseEnergy(const Particle* other) const;
+  bool operator==(const Particle& other) const {
+    if(name_ != other.name_) return false;
+    if(s_param_ == NULL && other.s_param_ == NULL) return true;
+    if(s_param_ != NULL) return s_param_->EqualTo(other.s_param_);
+    return false;
+  }
 
-  Real PairwiseEnergy(const ParticlePtr other) const {
-  	return this->PairwiseEnergy(other.get());
+  bool operator!=(const Particle& other) const {
+     return !(*this == other);
   }
 
-  const geom::Vec3& GetPos() const { return pos_; }
+  Real PairwiseScore(const Particle& other) const {
+    if(s_param_ == NULL || other.s_param_ == NULL) {
+      throw promod3::Error("Try to get pairwise score of particles that are "
+                           "not fully initialized");
+    }
+    return PairwiseParticleScore(s_param_, other.s_param_);
+  }
 
-  Real GetCharge() { return charge_; }
+  Real GetCollisionDistance() const { 
+    if(s_param_ != NULL) {
+      return s_param_->GetCollisionDistance(); 
+    } else {
+      throw promod3::Error("Try to get scoring parameters from a particle "
+                           "that is not fully initialized");
+    }
+  }
 
-  SidechainParticle GetParticleType() { return particle_id_; }
+  const geom::Vec3& GetPos() const { 
+    if(s_param_ != NULL) {
+      return s_param_->GetPos(); 
+    } else {
+      throw promod3::Error("Try to get scoring parameters from a particle "
+                           "that is not fully initialized");
+    }
+  } 
 
   const String& GetName() const { return name_; }
 
-  bool IsHBondDonor() const { return is_hbond_donor_; }
-
-  bool IsHBondAcceptor() const { return is_hbond_acceptor_; }
+  void SetName(const String& name) { name_ = name; }
 
-  const std::vector<geom::Vec3>& GetLonePairs() const { return lone_pairs_; }
+  PScoringParam* GetSParam() const { return s_param_; }
 
-  const geom::Vec3& GetPolarDirection() const { return polar_direction_; }
+  void SetSParam(PScoringParam* s_param) { s_param_ = s_param; }
 
-  void AddLonePair(const geom::Vec3& lone_pair){
-    is_hbond_acceptor_ = true;
-    lone_pairs_.push_back(geom::Normalize(lone_pair));
+  PScoringFunction GetScoringFunction() const {
+    if(s_param_ != NULL) {
+      return s_param_->GetScoringFunction(); 
+    } else {
+      throw promod3::Error("Try to get scoring function from a particle "
+                           "that is not fully initialized");
+    }    
   }
 
-  void SetPolarDirection(const geom::Vec3& dir){
-    is_hbond_donor_ = true;
-    polar_direction_ = geom::Normalize(dir);
+  void ApplyTransform(geom::Mat4& transform) {
+    if(s_param_ != NULL) {
+      return s_param_->ApplyTransform(transform); 
+    } else {
+      throw promod3::Error("Try to transform scoring parameters of a particle "
+                           "that is not fully initialized");
+    }    
   }
 
 private:
-
-  SidechainParticle particle_id_;
-  geom::Vec3 pos_;
-  Real charge_;
-  String name_; 
-  bool is_hbond_donor_;
-  bool is_hbond_acceptor_; 
-  std::vector<geom::Vec3> lone_pairs_;
-  geom::Vec3 polar_direction_;
+  String name_;
+  PScoringParam* s_param_;
 };
 
 }}//ns
diff --git a/sidechain/src/particle_scoring.cc b/sidechain/src/particle_scoring.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c95fc241162e7acf48234422e7bf71537cb484ee
--- /dev/null
+++ b/sidechain/src/particle_scoring.cc
@@ -0,0 +1,97 @@
+// Copyright (c) 2013-2018, SIB - Swiss Institute of Bioinformatics and 
+//                          Biozentrum - University of Basel
+// 
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+// 
+//   http://www.apache.org/licenses/LICENSE-2.0
+// 
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include <promod3/core/message.hh>
+#include <promod3/sidechain/particle_scoring.hh>
+#include <promod3/sidechain/scwrl4_particle_scoring.hh>
+
+namespace promod3{ namespace sidechain{
+
+Real PairwiseParticleScore(PScoringParam* p_one, PScoringParam* p_two) {
+
+  if(p_one->GetScoringFunction() != p_two->GetScoringFunction()) {
+    throw promod3::Error("Try to calculate pairwise particle energy with "
+                         "parameters of different type!");
+  }
+
+  switch(p_one->GetScoringFunction()) {
+    case SCWRL4: {
+      return SCWRL4PairwiseScore(reinterpret_cast<SCWRL4Param*>(p_one),
+                                 reinterpret_cast<SCWRL4Param*>(p_two));
+    }
+    default: {
+      throw promod3::Error("Don't know what scoring function to call between "
+                           "particles... maybe a lazy developer forgot to "
+                           "update the PairwiseParticleEnergy function after "
+                           "implementing a new score");
+    }
+  }
+}
+
+
+void PScoringParamContainer::push_back(PScoringParam* p) {
+  if(params_.empty()) {
+    scoring_function_ = p->GetScoringFunction();
+  } else{
+    if(scoring_function_ != p->GetScoringFunction()) {
+      throw promod3::Error("Tried to setup PScoringParamContainer with "
+                            "inconsistent scoring functions");
+    }
+  }
+  params_.push_back(p);
+}
+
+
+void PScoringParamContainer::clear() {
+  scoring_function_ = INVALID_PSCORING_FUNCTION;
+  params_.clear();
+}
+
+
+void PScoringParamContainer::PairwiseScores(const PScoringParamContainer& other,
+                                   const std::vector<std::pair<int, int> >& idx, 
+                                   std::vector<Real>& scores) const {
+
+  if(idx.empty()) {
+    scores.clear();
+    return;
+  }
+
+  if(scoring_function_ != other.GetScoringFunction()) {
+    throw promod3::Error("Tried to get pairwise scores from two "
+                         "PScoringParamContainers with inconsistent energy "
+                         "functions");
+  }
+
+  switch(scoring_function_) {
+    case SCWRL4: {
+      scores.resize(idx.size());
+      for(uint i = 0; i < idx.size(); ++i) {
+        scores[i] = 
+        SCWRL4PairwiseScore(reinterpret_cast<SCWRL4Param*>(params_[idx[i].first]),
+                            reinterpret_cast<SCWRL4Param*>(other.params_[idx[i].second]));
+      }
+      break;
+    }
+    default: {
+      throw promod3::Error("Don't know what scoring function to call between "
+                           "particles... maybe a lazy developer forgot to "
+                           "update the PairwiseParticleEnergy function after "
+                           "implementing a new score");
+    }
+  }
+}
+
+}} //ns
diff --git a/sidechain/src/particle_scoring.hh b/sidechain/src/particle_scoring.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b04fc8c04dde1628576f92e1018c54a0c8480dd8
--- /dev/null
+++ b/sidechain/src/particle_scoring.hh
@@ -0,0 +1,60 @@
+// Copyright (c) 2013-2018, SIB - Swiss Institute of Bioinformatics and 
+//                          Biozentrum - University of Basel
+// 
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+// 
+//   http://www.apache.org/licenses/LICENSE-2.0
+// 
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef PROMOD3_PARTICLE_SCORING_HH
+#define PROMOD3_PARTICLE_SCORING_HH
+
+#include <vector>
+#include <promod3/sidechain/particle_scoring_base.hh>
+
+namespace promod3{ namespace sidechain{
+
+// The generic pairwise energy function that calls the right pairwise energy 
+// function. Be aware of the virtual function calls to check what scoring 
+// function should be used. If you already know what you're dealing with
+// you should directly call the specific scoring function. reinterpret_cast the 
+// ptrs if necessary but only do that if you know what you're doing...
+Real PairwiseParticleScore(PScoringParam* p_one, PScoringParam* p_two);
+
+// this class is for internal use only and serves the RotamerGroup objects to
+// avoid too many virtual function calls. Its a container for pointeres to
+// PScoringParam objects. It guarantees that they're the same type. Virtual 
+// function calls to check for the underlying scoring function are therefore
+// not necessary anymore.
+class PScoringParamContainer{
+
+public:
+
+  PScoringParamContainer(): scoring_function_(INVALID_PSCORING_FUNCTION) { }
+
+  void push_back(PScoringParam* p);
+
+  void clear();
+
+  PScoringFunction GetScoringFunction() const { return scoring_function_; }
+
+  void PairwiseScores(const PScoringParamContainer& other,
+                      const std::vector<std::pair<int, int> >& idx, 
+                      std::vector<Real>& scores) const;
+
+private:
+
+  PScoringFunction scoring_function_;
+  std::vector<PScoringParam*> params_;
+};
+
+}} //ns
+
+#endif
diff --git a/sidechain/src/particle_scoring_base.hh b/sidechain/src/particle_scoring_base.hh
new file mode 100644
index 0000000000000000000000000000000000000000..96cdb6d97b80353e3a5331d394d6e707db9ae501
--- /dev/null
+++ b/sidechain/src/particle_scoring_base.hh
@@ -0,0 +1,53 @@
+// Copyright (c) 2013-2018, SIB - Swiss Institute of Bioinformatics and 
+//                          Biozentrum - University of Basel
+// 
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+// 
+//   http://www.apache.org/licenses/LICENSE-2.0
+// 
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef PROMOD3_PARTICLE_SCORING_BASE_HH
+#define PROMOD3_PARTICLE_SCORING_BASE_HH
+
+#include <ost/geom/vec3.hh>
+#include <ost/geom/mat4.hh>
+
+
+namespace promod3{ namespace sidechain{
+
+// the available particle scoring functions
+enum PScoringFunction{
+SCWRL4,
+INVALID_PSCORING_FUNCTION
+};
+
+
+// abstract base struct to be implemented by all particle scoring function 
+// specific parameter classes
+struct PScoringParam {
+
+  virtual ~PScoringParam() { }
+
+  virtual const geom::Vec3& GetPos() const = 0;
+
+  virtual Real GetCollisionDistance() const = 0;
+
+  virtual void ApplyTransform(const geom::Mat4& transform) = 0;
+
+  virtual PScoringParam* Copy() const = 0;
+
+  virtual bool EqualTo(PScoringParam* other) const = 0;
+  
+  virtual PScoringFunction GetScoringFunction() const =  0;
+};
+
+}} //ns
+
+#endif
diff --git a/sidechain/src/rotamer_constructor.cc b/sidechain/src/rotamer_constructor.cc
index 7967eecddb251c94c790684c73343ede4c16d684..3dfe130fcb331240a44ae2ab2b691af6dd4a95af 100644
--- a/sidechain/src/rotamer_constructor.cc
+++ b/sidechain/src/rotamer_constructor.cc
@@ -246,7 +246,7 @@ void RotamerConstructor::SetPosBuffer(const ost::mol::ResidueHandle& res,
                                            RotamerID id) {
 
   core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
-          "SCWRLRotamerConstructor::SetPosBuffer_residue", 2);
+          "SCWRL4RotamerConstructor::SetPosBuffer_residue", 2);
 
   ost::mol::AtomHandle n = res.FindAtom("N");
   ost::mol::AtomHandle ca = res.FindAtom("CA");
@@ -277,7 +277,7 @@ void RotamerConstructor::SetPosBuffer(const promod3::loop::AllAtomPositions& all
                                            uint all_atom_idx, RotamerID id) {
 
   core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
-          "SCWRLRotamerConstructor::SetPosBuffer_all_atom_pos", 2);
+          "SCWRL4RotamerConstructor::SetPosBuffer_all_atom_pos", 2);
 
 
   if(!(all_atom_pos.IsSet(all_atom_idx, promod3::loop::BB_N_INDEX) &&
diff --git a/sidechain/src/rotamer_constructor.hh b/sidechain/src/rotamer_constructor.hh
index 86b0c6f8cb2ce4d5d8c7f9fcda4e382b954a5edf..4a496ef0e13e68be29730a3be01c04022e18c2ec 100644
--- a/sidechain/src/rotamer_constructor.hh
+++ b/sidechain/src/rotamer_constructor.hh
@@ -18,14 +18,11 @@
 
 #include <vector>
 #include <ost/base.hh>
-//#include <ost/conop/amino_acids.hh>
 #include <ost/conop/compound_lib.hh>
 #include <promod3/loop/amino_acid_atoms.hh>
 #include <promod3/loop/all_atom_positions.hh>
 #include <promod3/loop/hydrogen_constructor.hh>
-#include <promod3/sidechain/particle.hh>
 #include <promod3/sidechain/rotamer_id.hh>
-#include <promod3/sidechain/rotamer.hh>
 #include <promod3/sidechain/rotamer_group.hh>
 #include <promod3/sidechain/frame.hh>
 #include <promod3/sidechain/rotamer_lib.hh>
diff --git a/sidechain/src/rotamer_group.cc b/sidechain/src/rotamer_group.cc
index ed6fde6c96b03d07a152f5211bf6f9102ca55e05..e963a685c7e86e6319a71f81df8b6a7f922dd018 100644
--- a/sidechain/src/rotamer_group.cc
+++ b/sidechain/src/rotamer_group.cc
@@ -16,6 +16,7 @@
 
 #include <promod3/sidechain/rotamer_group.hh>
 #include <promod3/core/runtime_profiling.hh>
+#include <promod3/sidechain/particle_scoring.hh>
 
 namespace{
 
@@ -110,17 +111,18 @@ void RRMRotamerGroup::AddFrameEnergy(FramePtr frame) {
   collision_tree_.OverlappingPolytopes(frame->GetCollisionTree(), 
                                        overlapping_particles);
 
+  if(overlapping_particles.empty()) return;
+
   const std::vector<uint>& frame_residue_indices = frame->GetResidueIndices();
-  const std::vector<Particle*>& frame_particles = frame->GetParticles();
   std::vector<Real> frame_energies(this->size(), 0.0);
-
-  for(std::vector<std::pair<int,int> >::iterator i = 
-      overlapping_particles.begin(),
-      e = overlapping_particles.end(); i != e; ++i) {
-    if(frame_residue_indices[i->second] != residue_index_) {
-      Real energy = 
-      particles_[i->first]->PairwiseEnergy(frame_particles[i->second]);
-      frame_energies[rotamer_indices_[i->first]] += energy;
+  std::vector<Real> scores;
+  scoring_parameters_.PairwiseScores(frame->GetScoringParameters(), 
+                                     overlapping_particles, scores);
+
+  for(uint i = 0; i < scores.size(); ++i) {
+    const std::pair<int, int>& idx = overlapping_particles[i];
+    if(frame_residue_indices[idx.second] != residue_index_) {
+      frame_energies[rotamer_indices_[idx.first]] += scores[i];
     }
   }
 
@@ -148,17 +150,15 @@ bool RRMRotamerGroup::CalculatePairwiseEnergies(RRMRotamerGroupPtr other,
   if(overlapping_particles.empty()) return false;
   
   emat = promod3::core::EMatXX::Zero(this->size(), other->size());
-
-  for(std::vector<std::pair<int,int> >::iterator i = 
-      overlapping_particles.begin(),
-      e = overlapping_particles.end(); i!= e; ++i) {
-
-    Real energy = 
-    particles_[i->first]->PairwiseEnergy(other->particles_[i->second]);
-
-    if(energy != Real(0.0)) {
-      emat(rotamer_indices_[i->first], other->rotamer_indices_[i->second]) += 
-      energy;
+  std::vector<Real> scores;
+  scoring_parameters_.PairwiseScores(other->scoring_parameters_, 
+                                     overlapping_particles, scores);
+
+  for(uint i = 0; i < scores.size(); ++i) {
+    if(scores[i] != 0.0) {
+      const std::pair<int, int>& idx = overlapping_particles[i];
+      emat(rotamer_indices_[idx.first], other->rotamer_indices_[idx.second]) += 
+      scores[i];
     }
   }
 
@@ -226,7 +226,7 @@ void RRMRotamerGroup::Init() {
     overall_size += (*i)->size();
   }
 
-  particles_.resize(overall_size);
+  scoring_parameters_.clear();
   rotamer_indices_.resize(overall_size);
 
   std::vector<Real> x_particle_positions(overall_size);
@@ -243,7 +243,7 @@ void RRMRotamerGroup::Init() {
 
     for(uint j = 0; j < (*i)->size(); ++j) {
       Particle& p = (**i)[j];
-      particles_[actual_position] = &p;
+      scoring_parameters_.push_back(p.GetSParam());
       rotamer_indices_[actual_position] = actual_rotamer;
       pos = p.GetPos();
       x_particle_positions[actual_position] = pos[0];
@@ -331,8 +331,9 @@ void FRMRotamerGroup::AddFrameEnergy(FramePtr frame) {
   collision_tree_.OverlappingPolytopes(frame->GetCollisionTree(),
                                        overlapping_particles);
 
+  if(overlapping_particles.empty()) return;
+
   const std::vector<uint>& frame_residue_indices = frame->GetResidueIndices();
-  const std::vector<Particle*>& frame_particles = frame->GetParticles();
 
   std::vector<int> subrotamer_sizes(this->size());
   int max_subrotamer_size = 0;
@@ -354,20 +355,18 @@ void FRMRotamerGroup::AddFrameEnergy(FramePtr frame) {
     }
   }
 
-  for(std::vector<std::pair<int,int> >::iterator i = 
-      overlapping_particles.begin(),
-      e = overlapping_particles.end(); i != e ; ++i) {
-
-    if(frame_residue_indices[i->second] != residue_index_) {
-
-      Real energy = 
-      particles_[i->first]->PairwiseEnergy(frame_particles[i->second]);
+  std::vector<Real> scores;
+  scoring_parameters_.PairwiseScores(frame->GetScoringParameters(), 
+                                     overlapping_particles, scores);
 
-      if(energy != Real(0.0)) {
-        int row_idx = rotamer_indices_[i->first];
+  for(uint i = 0; i < scores.size(); ++i) {
+    const std::pair<int, int>& idx = overlapping_particles[i];
+    if(frame_residue_indices[idx.second] != residue_index_) {
+      if(scores[i] != 0.0) {
+        int row_idx = rotamer_indices_[idx.first];
         for(FRMRotamer::subrotamer_association_iterator j = 
-            ass_ind_b_[i->first]; j != ass_ind_e_[i->first]; ++j) {
-          frame_energies(row_idx, *j) += energy;
+            ass_ind_b_[idx.first]; j != ass_ind_e_[idx.first]; ++j) {
+          frame_energies(row_idx, *j) += scores[i];
         }
       }
     }
@@ -450,36 +449,34 @@ bool FRMRotamerGroup::CalculatePairwiseEnergies(FRMRotamerGroupPtr other,
   promod3::core::EMatXX::Zero(rotamer_combinations, 
                               max_subrotamer_combinations);
 
-  Real energy;
   int r_index_j;
   int e_row_idx, e_col_idx;
   int other_subrotamer_size;
   FRMRotamer::subrotamer_association_iterator subrot_i_b, subrot_i_e,
                                               subrot_j_b, subrot_j_e, temp_it;
 
-  for(std::vector<std::pair<int,int> >::iterator i = 
-      overlapping_particles.begin(),
-      e = overlapping_particles.end(); i != e; ++i) {
-    
-    energy = particles_[i->first]->PairwiseEnergy(other->particles_[i->second]);
-
-    if(energy != Real(0.0)) {
+  std::vector<Real> scores;
+  scoring_parameters_.PairwiseScores(other->scoring_parameters_, 
+                                     overlapping_particles, scores);
 
+  for(uint i = 0; i < scores.size(); ++i) {
+    if(scores[i] != 0.0) {
+      const std::pair<int, int>& idx = overlapping_particles[i];
       // we have to add this energy to all combinations of subrotamers where
       // this pair of particles is involved in...
-      r_index_j = other->rotamer_indices_[i->second];
-      e_row_idx = rotamer_indices_[i->first] * other_size + r_index_j;
+      r_index_j = other->rotamer_indices_[idx.second];
+      e_row_idx = rotamer_indices_[idx.first] * other_size + r_index_j;
       other_subrotamer_size = other_subrotamer_sizes[r_index_j];
-      subrot_i_b = ass_ind_b_[i->first];
-      subrot_i_e = ass_ind_e_[i->first];
-      temp_it = other->ass_ind_b_[i->second];
-      subrot_j_e = other->ass_ind_e_[i->second];
+      subrot_i_b = ass_ind_b_[idx.first];
+      subrot_i_e = ass_ind_e_[idx.first];
+      temp_it = other->ass_ind_b_[idx.second];
+      subrot_j_e = other->ass_ind_e_[idx.second];
 
       for(;subrot_i_b != subrot_i_e; ++subrot_i_b) {
         subrot_j_b = temp_it;
         int temp = (*subrot_i_b) * other_subrotamer_size;
         for(;subrot_j_b != subrot_j_e; ++subrot_j_b) {
-          subrotamer_energies(e_row_idx, temp + (*subrot_j_b)) += energy;
+          subrotamer_energies(e_row_idx, temp + (*subrot_j_b)) += scores[i];
         }
       }
     }
@@ -593,13 +590,12 @@ void FRMRotamerGroup::ApplySelfEnergyThresh(Real self_energy_thresh) {
 void FRMRotamerGroup::Init() {
 
   int overall_size = 0;
-
   for(std::vector<FRMRotamerPtr>::iterator i = rotamers_.begin(); 
       i != rotamers_.end(); ++i) {
     overall_size += (*i)->size();
   }
 
-  particles_.resize(overall_size);
+  scoring_parameters_.clear();
   rotamer_indices_.resize(overall_size);
 
   ass_ind_b_.resize(overall_size);
@@ -622,7 +618,7 @@ void FRMRotamerGroup::Init() {
 
     for(uint j = 0; j < (*i)->size(); ++j) {
       Particle& p = (**i)[j];
-      particles_[actual_position] = &p;
+      scoring_parameters_.push_back(p.GetSParam());
       rotamer_indices_[actual_position] = actual_rotamer;
       ass_ind_b_[actual_position] = 
       (*i)->subrotamer_association_begin(actual_internal_index);
diff --git a/sidechain/src/rotamer_group.hh b/sidechain/src/rotamer_group.hh
index b2556becc273619a187b085275e7fde33395dc08..eae0f2a5c17a68d34435cfaaf5e52feeead3bed3 100644
--- a/sidechain/src/rotamer_group.hh
+++ b/sidechain/src/rotamer_group.hh
@@ -25,6 +25,7 @@
 #include <promod3/core/eigen_types.hh>
 #include <promod3/sidechain/rotamer.hh>
 #include <promod3/sidechain/frame.hh>
+#include <promod3/sidechain/particle_scoring.hh>
 
 namespace promod3 { namespace sidechain {
 
@@ -95,7 +96,7 @@ private:
   void Init();
 
   std::vector<RRMRotamerPtr> rotamers_;
-  std::vector<Particle*> particles_;
+  PScoringParamContainer scoring_parameters_;
   std::vector<int> rotamer_indices_;
   promod3::core::TetrahedralPolytopeTree collision_tree_;
   uint residue_index_;
@@ -165,7 +166,7 @@ private:
   void Init();
   
   std::vector<FRMRotamerPtr> rotamers_;
-  std::vector<Particle*> particles_;
+  PScoringParamContainer scoring_parameters_;
   std::vector<int> rotamer_indices_;
   std::vector<FRMRotamer::subrotamer_association_iterator> ass_ind_b_;
   std::vector<FRMRotamer::subrotamer_association_iterator> ass_ind_e_;
diff --git a/sidechain/src/particle.cc b/sidechain/src/scwrl4_particle_scoring.cc
similarity index 97%
rename from sidechain/src/particle.cc
rename to sidechain/src/scwrl4_particle_scoring.cc
index e4c79f01a01b4216a97b30cc3c53b0dedb2dc785..1e93f0b01c9bd296e0ba3f99c55415716c1f4cab 100644
--- a/sidechain/src/particle.cc
+++ b/sidechain/src/scwrl4_particle_scoring.cc
@@ -14,8 +14,7 @@
 // limitations under the License.
 
 
-#include <promod3/sidechain/particle.hh>
-
+#include <promod3/sidechain/scwrl4_particle_scoring.hh>
 
 namespace{
 
@@ -4595,83 +4594,131 @@ namespace{
     int bin = std::min(499,static_cast<int>(d*100));
     return _energies[a][b][bin];
   }
+} // anon ns
+
+
+namespace promod3 { namespace sidechain {
+
+void SCWRL4Param::ApplyTransform(const geom::Mat4& t) {
+
+  // get transformed pos
+  Real a = t(0,0)*pos_[0] + t(0,1)*pos_[1] + t(0,2)*pos_[2] + t(0,3);
+  Real b = t(1,0)*pos_[0] + t(1,1)*pos_[1] + t(1,2)*pos_[2] + t(1,3);
+  Real c = t(2,0)*pos_[0] + t(2,1)*pos_[1] + t(2,2)*pos_[2] + t(2,3);
+  geom::Vec3 transformed_pos(a, b, c);
+
+  if(is_hbond_donor_) {
+    geom::Vec3 pos = pos_ + polar_direction_;
+    Real a = t(0,0)*pos[0] + t(0,1)*pos[1] + t(0,2)*pos[2] + t(0,3);
+    Real b = t(1,0)*pos[0] + t(1,1)*pos[1] + t(1,2)*pos[2] + t(1,3);
+    Real c = t(2,0)*pos[0] + t(2,1)*pos[1] + t(2,2)*pos[2] + t(2,3);
+    polar_direction_ = geom::Normalize(geom::Vec3(a,b,c) - transformed_pos);
+  }
+
+  if(is_hbond_acceptor_) {
+    for(uint i = 0; i < lone_pairs_.size(); ++i) {
+      geom::Vec3 pos = pos_ + lone_pairs_[i];
+      Real a = t(0,0)*pos[0] + t(0,1)*pos[1] + t(0,2)*pos[2] + t(0,3);
+      Real b = t(1,0)*pos[0] + t(1,1)*pos[1] + t(1,2)*pos[2] + t(1,3);
+      Real c = t(2,0)*pos[0] + t(2,1)*pos[1] + t(2,2)*pos[2] + t(2,3);
+      lone_pairs_[i] = geom::Normalize(geom::Vec3(a,b,c) - transformed_pos);
+    }
+  }
+
+  // do the position update
+  pos_ = transformed_pos;
+}
+
+
+SCWRL4Param* SCWRL4Param::Copy() const {
+  SCWRL4Param* return_ptr = new SCWRL4Param(particle_type_, pos_, charge_);
+  return_ptr->is_hbond_donor_ = is_hbond_donor_;
+  return_ptr->is_hbond_acceptor_ = is_hbond_acceptor_;
+  return_ptr->lone_pairs_ = lone_pairs_;
+  return_ptr->polar_direction_ = polar_direction_;
+  return return_ptr;
 }
 
 
+bool SCWRL4Param::EqualTo(PScoringParam* other) const {
 
+  if(other == NULL) return false;
+  if(other->GetScoringFunction() != SCWRL4) return false;
 
-namespace promod3 { namespace sidechain {
+  // as the other scoring function is also SCWRL4, we can assume that
+  // the following dynamic cast is fine...
+  SCWRL4Param* p = dynamic_cast<SCWRL4Param*>(other);
 
-bool Particle::operator==(const Particle& particle) const {   
-  return (particle_id_ == particle.particle_id_ && 
-          pos_ == particle.pos_ && 
-          charge_ == particle.charge_ && 
-          name_ == particle.name_ && 
-          is_hbond_donor_ == particle.is_hbond_donor_ && 
-          is_hbond_acceptor_ == particle.is_hbond_acceptor_ && 
-          lone_pairs_ == particle.lone_pairs_ && 
-          polar_direction_ == particle.polar_direction_); 
+  return pos_ == p->pos_ &&
+         collision_distance_ == p->collision_distance_ &&
+         particle_type_ == p->particle_type_ &&
+         charge_ == p->charge_ &&
+         is_hbond_donor_ == p->is_hbond_donor_ &&
+         is_hbond_acceptor_ == p->is_hbond_acceptor_ &&
+         lone_pairs_ == p->lone_pairs_ &&
+         polar_direction_ == p->polar_direction_;
 }
 
-Real Particle::PairwiseEnergy(const Particle* other) const{
 
+Real SCWRL4PairwiseScore(SCWRL4Param* p_one, SCWRL4Param* p_two) {
 
-  Real d = geom::Distance(pos_,other->pos_);
-  Real e = _GetEnergy(particle_id_, other->particle_id_, d);
-  if (is_hbond_donor_) {
-    if (other->is_hbond_acceptor_) {
-      Real distance_term = Real(0.44890) - (d-Real(2.08))*(d-Real(2.08));
+  Real d = geom::Distance(p_one->pos_, p_two->pos_);
+  Real e = _GetEnergy(p_one->particle_type_, p_two->particle_type_, d);
+  if (p_one->is_hbond_donor_) {
+    if (p_two->is_hbond_acceptor_) {
+      Real distance_term = Real(0.44890) - (d - Real(2.08))*(d - Real(2.08));
       if (distance_term < 0) return e;
       d = Real(1.0) / d;
-      Real n[] = {(pos_[0]-other->pos_[0])*d,
-                  (pos_[1]-other->pos_[1])*d,
-                  (pos_[2]-other->pos_[2])*d};
-      Real cos_a_term = -n[0]*polar_direction_[0] - n[1]*polar_direction_[1]
-                        -n[2]*polar_direction_[2] - Real(0.79864);
+      Real n[] = {(p_one->pos_[0] - p_two->pos_[0])*d,
+                  (p_one->pos_[1] - p_two->pos_[1])*d,
+                  (p_one->pos_[2] - p_two->pos_[2])*d};
+      Real cos_a_term = -n[0]*p_one->polar_direction_[0] 
+                        -n[1]*p_one->polar_direction_[1]
+                        -n[2]*p_one->polar_direction_[2] - Real(0.79864);
       if (cos_a_term < 0) return e;
       // find lone pair, that gives maximal value to the final nominator term
-      Real cos_b_term = n[0]*other->lone_pairs_[0][0]
-                      + n[1]*other->lone_pairs_[0][1]
-                      + n[2]*other->lone_pairs_[0][2] - Real(0.65606);
-      for (uint i = 1; i < other->lone_pairs_.size(); ++i) {
-        d = n[0]*other->lone_pairs_[i][0] + n[1]*other->lone_pairs_[i][1]
-          + n[2]*other->lone_pairs_[i][2] - Real(0.65606);
+      Real cos_b_term = n[0]*p_two->lone_pairs_[0][0]
+                       +n[1]*p_two->lone_pairs_[0][1]
+                       +n[2]*p_two->lone_pairs_[0][2] - Real(0.65606);
+      for (uint i = 1; i < p_two->lone_pairs_.size(); ++i) {
+        d = n[0]*p_two->lone_pairs_[i][0] + n[1]*p_two->lone_pairs_[i][1]
+          + n[2]*p_two->lone_pairs_[i][2] - Real(0.65606);
         cos_b_term = (cos_b_term < d) ? d : cos_b_term;
       }
       if (cos_b_term < 0) return e;
       d = std::sqrt(distance_term*cos_a_term*cos_b_term) * Real(6.5904);
-      return (Real(1)-d)*e + d*Real(35)*charge_*other->charge_;
+      return (Real(1)-d)*e + d*Real(35)*p_one->charge_*p_two->charge_;
     }
     else return e;
   }
-  if (is_hbond_acceptor_) {
-    if (other->is_hbond_donor_) {
-      Real distance_term = Real(0.44890) - (d-Real(2.08))*(d-Real(2.08));
+  if (p_one->is_hbond_acceptor_) {
+    if (p_two->is_hbond_donor_) {
+      Real distance_term = Real(0.44890) - (d - Real(2.08))*(d - Real(2.08));
       if (distance_term < 0) return e;
       d = Real(1.0) / d;
-      Real n[] = {(other->pos_[0]-pos_[0])*d,
-                  (other->pos_[1]-pos_[1])*d,
-                  (other->pos_[2]-pos_[2])*d};
-      Real cos_a_term = -n[0]*other->polar_direction_[0]
-                      - n[1]*other->polar_direction_[1]
-                      - n[2]*other->polar_direction_[2] - Real(0.79864);
+      Real n[] = {(p_two->pos_[0] - p_one->pos_[0])*d,
+                  (p_two->pos_[1] - p_one->pos_[1])*d,
+                  (p_two->pos_[2] - p_one->pos_[2])*d};
+      Real cos_a_term = -n[0]*p_two->polar_direction_[0]
+                        -n[1]*p_two->polar_direction_[1]
+                        -n[2]*p_two->polar_direction_[2] - Real(0.79864);
       if (cos_a_term < 0) return e;
       //find lone pair, that gives maximal value to the final nominator term
-      Real cos_b_term = n[0]*lone_pairs_[0][0] + n[1]*lone_pairs_[0][1]
-                      + n[2]*lone_pairs_[0][2] - Real(0.65606);
-      for (uint i = 1; i < lone_pairs_.size(); ++i) {
-        d = n[0]*lone_pairs_[i][0] + n[1]*lone_pairs_[i][1]
-          + n[2]*lone_pairs_[i][2] - Real(0.65606);
+      Real cos_b_term = n[0]*p_one->lone_pairs_[0][0] 
+                       +n[1]*p_one->lone_pairs_[0][1]
+                       +n[2]*p_one->lone_pairs_[0][2] - Real(0.65606);
+      for (uint i = 1; i < p_one->lone_pairs_.size(); ++i) {
+        d = n[0]*p_one->lone_pairs_[i][0] + n[1]*p_one->lone_pairs_[i][1]
+          + n[2]*p_one->lone_pairs_[i][2] - Real(0.65606);
         cos_b_term = (cos_b_term < d) ? d : cos_b_term;
       }
       if (cos_b_term < 0) return e;
       d = std::sqrt(distance_term*cos_a_term*cos_b_term) * Real(6.5904);
-      return (Real(1)-d)*e + d*Real(35)*charge_*other->charge_;
+      return (Real(1)-d)*e + d*Real(35)*p_one->charge_*p_two->charge_;
     }
     else return e;
   }
   return e;
-
 }
 
 }}//ns
diff --git a/sidechain/src/scwrl4_particle_scoring.hh b/sidechain/src/scwrl4_particle_scoring.hh
new file mode 100644
index 0000000000000000000000000000000000000000..cb3af990d2bb3973773c4a45b0c9a90cd1245c72
--- /dev/null
+++ b/sidechain/src/scwrl4_particle_scoring.hh
@@ -0,0 +1,101 @@
+// Copyright (c) 2013-2018, SIB - Swiss Institute of Bioinformatics and 
+//                          Biozentrum - University of Basel
+// 
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+// 
+//   http://www.apache.org/licenses/LICENSE-2.0
+// 
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef PROMOD3_SCWRL4_PARTICLE_SCORING_HH
+#define PROMOD3_SCWRL4_PARTICLE_SCORING_HH
+
+#include <ost/geom/mat4.hh>
+#include <ost/geom/vecmat3_op.hh>
+#include <promod3/sidechain/particle_scoring_base.hh>
+
+namespace promod3{ namespace sidechain{
+
+enum SCWRL4ParticleType {
+  HParticle,
+  CParticle,
+  CH1Particle,
+  CH2Particle,
+  CH3Particle,
+  NParticle,
+  OParticle,
+  OCParticle,
+  SParticle,
+};
+
+const Real SCWRL4CollisionDistances[9] = {1.0300, 2.2400, 2.2400, 2.2400, 
+                                          2.2400, 1.7067, 1.7200,
+                                          1.7200, 2.2800};
+
+struct SCWRL4Param : public PScoringParam {
+
+  SCWRL4Param(SCWRL4ParticleType p, const geom::Vec3& pos, 
+              Real charge): pos_(pos), 
+                            collision_distance_(SCWRL4CollisionDistances[p]),
+                            particle_type_(p), charge_(charge),
+                            is_hbond_donor_(false), 
+                            is_hbond_acceptor_(false) { }
+
+  virtual ~SCWRL4Param() { }
+
+  virtual const geom::Vec3& GetPos() const { return pos_; }
+
+  virtual Real GetCollisionDistance() const { return collision_distance_; }
+
+  virtual void ApplyTransform(const geom::Mat4& t);
+
+  virtual SCWRL4Param* Copy() const;
+
+  virtual bool EqualTo(PScoringParam* other) const;
+
+  virtual PScoringFunction GetScoringFunction() const { return SCWRL4; }
+
+  Real GetCharge() { return charge_; }
+
+  SCWRL4ParticleType GetParticleType() { return particle_type_; }
+
+  bool IsHBondDonor() const { return is_hbond_donor_; }
+
+  bool IsHBondAcceptor() const { return is_hbond_acceptor_; }
+
+  const std::vector<geom::Vec3>& GetLonePairs() const { return lone_pairs_; }
+
+  const geom::Vec3& GetPolarDirection() const { return polar_direction_; }
+
+  void AddLonePair(const geom::Vec3& lone_pair){
+    is_hbond_acceptor_ = true;
+    lone_pairs_.push_back(geom::Normalize(lone_pair));
+  }
+
+  void SetPolarDirection(const geom::Vec3& dir){
+    is_hbond_donor_ = true;
+    polar_direction_ = geom::Normalize(dir);
+  }
+
+  geom::Vec3 pos_;
+  Real collision_distance_;
+  SCWRL4ParticleType particle_type_;
+  Real charge_;
+  bool is_hbond_donor_;
+  bool is_hbond_acceptor_; 
+  std::vector<geom::Vec3> lone_pairs_;
+  geom::Vec3 polar_direction_;
+};
+
+
+Real SCWRL4PairwiseScore(SCWRL4Param* p_one, SCWRL4Param* p_two);
+
+}} //ns
+
+#endif
diff --git a/sidechain/src/scwrl_rotamer_constructor.cc b/sidechain/src/scwrl4_rotamer_constructor.cc
similarity index 92%
rename from sidechain/src/scwrl_rotamer_constructor.cc
rename to sidechain/src/scwrl4_rotamer_constructor.cc
index b0e1129429f3f2c869e114a49d9d0c0c101a9c02..7619d2ec8cef748701cf4ba50b046dd8cdc53ee4 100644
--- a/sidechain/src/scwrl_rotamer_constructor.cc
+++ b/sidechain/src/scwrl4_rotamer_constructor.cc
@@ -14,7 +14,8 @@
 // limitations under the License.
 
 
-#include <promod3/sidechain/scwrl_rotamer_constructor.hh>
+#include <promod3/sidechain/scwrl4_rotamer_constructor.hh>
+#include <promod3/sidechain/scwrl4_particle_scoring.hh>
 #include <promod3/core/runtime_profiling.hh>
 #include <promod3/core/message.hh>
 #include <promod3/loop/sidechain_atom_constructor.hh>
@@ -26,20 +27,20 @@ namespace{
 // add lone pairs for carbonyl (X-C-O) to O particle p
 void _AddLonePairsCarbonyl(const geom::Vec3& x_pos, const geom::Vec3& c_pos,
                            const geom::Vec3& o_pos, 
-                           promod3::sidechain::Particle& p) {
+                           promod3::sidechain::SCWRL4Param* p) {
   geom::Vec3 lone_pair_center_one, lone_pair_center_two;
   promod3::core::ConstructAtomPos(x_pos, c_pos, o_pos, 1.00, 2.0944,
                                   M_PI, lone_pair_center_one);
   promod3::core::ConstructAtomPos(x_pos, c_pos, o_pos, 1.00, 2.0944,
                                   0.0, lone_pair_center_two);
-  p.AddLonePair(lone_pair_center_one - o_pos);
-  p.AddLonePair(lone_pair_center_two - o_pos);
+  p->AddLonePair(lone_pair_center_one - o_pos);
+  p->AddLonePair(lone_pair_center_two - o_pos);
 }
 
 // add lone pairs for C-O-H combos (SER / THR / TYR) to O particle p
 void _AddLonePairsCOH(const geom::Vec3& c_pos, const geom::Vec3& o_pos,
                       const geom::Vec3& h_pos, 
-                      promod3::sidechain::Particle& p) {
+                      promod3::sidechain::SCWRL4Param* p) {
   // create a center point between the two lone pair clouds with distance to O
   // equals 1
   const geom::Vec3 center_point
@@ -53,19 +54,19 @@ void _AddLonePairsCOH(const geom::Vec3& c_pos, const geom::Vec3& o_pos,
                                   M_PI/2, lone_pair_center_one);
   promod3::core::ConstructAtomPos(c_pos, o_pos, center_point, 1.4150, M_PI/2,
                                   -M_PI/2, lone_pair_center_two);
-  p.AddLonePair(lone_pair_center_one - o_pos);
-  p.AddLonePair(lone_pair_center_two - o_pos);
+  p->AddLonePair(lone_pair_center_one - o_pos);
+  p->AddLonePair(lone_pair_center_two - o_pos);
 }
 
 // add lone pairs for C-N-C combos (HIS) to N particle p
 void _AddLonePairsCNC(const geom::Vec3& c1_pos, const geom::Vec3& n_pos,
                       const geom::Vec3& c2_pos, 
-                      promod3::sidechain::Particle& p) {
-  p.AddLonePair(  geom::Normalize(n_pos - c1_pos)
-                + geom::Normalize(n_pos - c2_pos));
+                      promod3::sidechain::SCWRL4Param* p) {
+  p->AddLonePair(  geom::Normalize(n_pos - c1_pos)
+                 + geom::Normalize(n_pos - c2_pos));
 }
 
-void EvalLonePairRule(promod3::sidechain::Particle& p,
+void EvalLonePairRule(promod3::sidechain::SCWRL4Param* p,
                       promod3::sidechain::SCWRLLPRule rule,
                       const geom::Vec3& a,
                       const geom::Vec3& b,
@@ -1398,18 +1399,18 @@ SCWRLRotamerLookup::SCWRLRotamerLookup(bool cb_in_sidechain) {
   }
 }
 
-SCWRLRotamerConstructor::SCWRLRotamerConstructor(bool cb_in_sidechain): 
+SCWRL4RotamerConstructor::SCWRL4RotamerConstructor(bool cb_in_sidechain): 
                                 rotamer_lookup_(cb_in_sidechain) {
   hydrogen_buffer_ = boost::make_shared<promod3::loop::HydrogenStorage>();
 }
 
 
-FrameResiduePtr SCWRLRotamerConstructor::ConstructBackboneFrameResidue(
+FrameResiduePtr SCWRL4RotamerConstructor::ConstructBackboneFrameResidue(
           const ost::mol::ResidueHandle& res, RotamerID id, uint residue_index,
           Real phi, bool n_ter, bool c_ter) {
 
   core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
-          "SCWRLRotamerConstructor::ConstructBackboneFrameResidue", 2);
+          "SCWRL4RotamerConstructor::ConstructBackboneFrameResidue", 2);
 
   if(id == XXX) {
     throw promod3::Error("Cannot construct FrameResidue for unknown "
@@ -1446,13 +1447,13 @@ FrameResiduePtr SCWRLRotamerConstructor::ConstructBackboneFrameResidue(
   return p;
 }
 
-FrameResiduePtr SCWRLRotamerConstructor::ConstructBackboneFrameResidue(
+FrameResiduePtr SCWRL4RotamerConstructor::ConstructBackboneFrameResidue(
           const promod3::loop::AllAtomPositions& all_atom, uint aa_res_idx, 
           RotamerID id, uint residue_index,
           Real phi, bool n_ter, bool c_ter) {
 
   core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
-          "SCWRLRotamerConstructor::ConstructBackboneFrameResidue", 2);
+          "SCWRL4RotamerConstructor::ConstructBackboneFrameResidue", 2);
 
   if(id == XXX) {
     throw promod3::Error("Cannot construct FrameResidue for unknown "
@@ -1491,11 +1492,11 @@ FrameResiduePtr SCWRLRotamerConstructor::ConstructBackboneFrameResidue(
   return p; 
 }
 
-FrameResiduePtr SCWRLRotamerConstructor::ConstructSidechainFrameResidue(
+FrameResiduePtr SCWRL4RotamerConstructor::ConstructSidechainFrameResidue(
         const ost::mol::ResidueHandle& res, RotamerID id, uint residue_index) {
 
   core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
-          "SCWRLRotamerConstructor::ConstructSidechainFrameResidue", 2);
+          "SCWRL4RotamerConstructor::ConstructSidechainFrameResidue", 2);
 
   if(id == XXX) {
     throw promod3::Error("Cannot construct FrameResidue for unknown "
@@ -1546,12 +1547,12 @@ FrameResiduePtr SCWRLRotamerConstructor::ConstructSidechainFrameResidue(
   return p;
 }
 
-FrameResiduePtr SCWRLRotamerConstructor::ConstructSidechainFrameResidue(
+FrameResiduePtr SCWRL4RotamerConstructor::ConstructSidechainFrameResidue(
           const promod3::loop::AllAtomPositions& all_atom, uint aa_res_idx, 
           RotamerID id, uint residue_index) {
 
   core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
-          "SCWRLRotamerConstructor::ConstructSidechainFrameResidue", 2);
+          "SCWRL4RotamerConstructor::ConstructSidechainFrameResidue", 2);
 
   if(id == XXX) {
     throw promod3::Error("Cannot construct FrameResidue for unknown "
@@ -1602,7 +1603,7 @@ FrameResiduePtr SCWRLRotamerConstructor::ConstructSidechainFrameResidue(
   return p; 
 }
 
-FrameResiduePtr SCWRLRotamerConstructor::ConstructFrameResidue(
+FrameResiduePtr SCWRL4RotamerConstructor::ConstructFrameResidue(
           const ost::mol::ResidueHandle& res, uint residue_index) {
 
   // get atoms and count non-hydrogens (= num. particles to add)
@@ -1611,8 +1612,9 @@ FrameResiduePtr SCWRLRotamerConstructor::ConstructFrameResidue(
   for (uint i = 0; i < atom_list.size(); ++i) {
     const String& ele = atom_list[i].GetElement();
     if(_IsHydrogen(ele)) continue;
-    Particle p(CParticle, atom_list[i].GetPos(), 0.0,
-               atom_list[i].GetName());
+    Particle p(atom_list[i].GetName(), new SCWRL4Param(CParticle, 
+                                                       atom_list[i].GetPos(), 
+                                                       0.0));
     particles.push_back(p);
   }
 
@@ -1622,7 +1624,7 @@ FrameResiduePtr SCWRLRotamerConstructor::ConstructFrameResidue(
 
 // NOTE: this only deals with few possible ligand cases and has not been
 //       tested extensively!
-FrameResiduePtr SCWRLRotamerConstructor::ConstructFrameResidueHeuristic(
+FrameResiduePtr SCWRL4RotamerConstructor::ConstructFrameResidueHeuristic(
           const ost::mol::ResidueHandle& res, uint residue_index,
           ost::conop::CompoundLibPtr comp_lib) {
 
@@ -1655,20 +1657,22 @@ FrameResiduePtr SCWRLRotamerConstructor::ConstructFrameResidueHeuristic(
       if (idx >= 0) {
         // create particle according to atom's element
         if (atom_elem == "C") {
-          SidechainParticle p_id = CParticle;
+          SCWRL4ParticleType p_id = CParticle;
           if      (atom_infos[idx].num_hydrogens == 1) p_id = CH1Particle;
           else if (atom_infos[idx].num_hydrogens == 2) p_id = CH2Particle;
           else if (atom_infos[idx].num_hydrogens == 3) p_id = CH3Particle;
-          Particle p(p_id, atom_pos, 0.0, atom_name);
+
+          Particle p(atom_name, new SCWRL4Param(p_id, atom_pos, 0.0));
           particles.push_back(p);
         } else if (atom_elem == "N") {
-          particles.push_back(Particle(NParticle, atom_pos, 0.0, atom_name));
+          Particle p(atom_name, new SCWRL4Param(NParticle, atom_pos, 0.0));
+          particles.push_back(p);
         } else if (atom_elem == "O") {
           // carbonyl?
           _AtomSpecCPList carbonyl_atoms = _GetCarboxylAtoms(atom_infos[idx]);
           if (carbonyl_atoms.size() == 2) {
             // charged particle needed for hydrogen bonds
-            Particle p(OParticle, atom_pos, -0.6, atom_name);
+            Particle p(atom_name, new SCWRL4Param(OParticle, atom_pos, -0.6));
             particles.push_back(p);
             // try to get lone pair directions (try name and alt. name)
             ost::mol::AtomHandle c = res.FindAtom(carbonyl_atoms[0]->name);
@@ -1677,22 +1681,22 @@ FrameResiduePtr SCWRLRotamerConstructor::ConstructFrameResidueHeuristic(
             if (!x.IsValid()) x = res.FindAtom(carbonyl_atoms[1]->alt_name);
             if (c.IsValid() && x.IsValid()) {
               _AddLonePairsCarbonyl(x.GetPos(), c.GetPos(), atom_pos,
-                                    particles[p_idx]);
+                  reinterpret_cast<SCWRL4Param*>(particles[p_idx].GetSParam()));
             }
           } else {
-            Particle p(OParticle, atom_pos, 0.0, atom_name);
+            Particle p(atom_name, new SCWRL4Param(OParticle, atom_pos, 0.0));
             particles.push_back(p);
           }
         } else if (atom_elem == "S") {
-          Particle p(SParticle, atom_pos, 0.0, atom_name);
+          Particle p(atom_name, new SCWRL4Param(SParticle, atom_pos, 0.0));
           particles.push_back(p);
         } else {
-          Particle p(CParticle, atom_pos, 0.0, atom_name);
+          Particle p(atom_name, new SCWRL4Param(CParticle, atom_pos, 0.0));
           particles.push_back(p);
         }
       } else {
         // fallback -> treat as uncharged carbon
-        Particle p(CParticle, atom_pos, 0.0, atom_name);
+        Particle p(atom_name, new SCWRL4Param(CParticle, atom_pos, 0.0));
         particles.push_back(p);
       }
       // particle was created anyways...
@@ -1705,10 +1709,10 @@ FrameResiduePtr SCWRLRotamerConstructor::ConstructFrameResidueHeuristic(
 
 
 
-void SCWRLRotamerConstructor::AssignInternalEnergies(RRMRotamerGroupPtr group) {
+void SCWRL4RotamerConstructor::AssignInternalEnergies(RRMRotamerGroupPtr group) {
 
   core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
-          "SCWRLRotamerConstructor::AssignInternalEnergies", 2);
+          "SCWRL4RotamerConstructor::AssignInternalEnergies", 2);
 
   Real max_p = group->GetMaxP();
   for(uint i = 0; i < group->size(); ++i){
@@ -1720,10 +1724,10 @@ void SCWRLRotamerConstructor::AssignInternalEnergies(RRMRotamerGroupPtr group) {
   }
 }
 
-void SCWRLRotamerConstructor::AssignInternalEnergies(FRMRotamerGroupPtr group) {
+void SCWRL4RotamerConstructor::AssignInternalEnergies(FRMRotamerGroupPtr group) {
 
   core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
-          "SCWRLRotamerConstructor::AssignInternalEnergies", 2);
+          "SCWRL4RotamerConstructor::AssignInternalEnergies", 2);
 
   Real max_p = group->GetMaxP();
   for(uint i = 0; i < group->size(); ++i){
@@ -1737,13 +1741,13 @@ void SCWRLRotamerConstructor::AssignInternalEnergies(FRMRotamerGroupPtr group) {
 
 // INTERNAL FUNCTIONS
 
-RRMRotamerGroupPtr SCWRLRotamerConstructor::ConstructRRMGroup(
+RRMRotamerGroupPtr SCWRL4RotamerConstructor::ConstructRRMGroup(
           RotamerID id, uint residue_index, 
           std::pair<RotamerLibEntry*,uint> lib_entries, 
           Real probability_cutoff) {
 
   core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
-          "SCWRLRotamerConstructor::ConstructRRMGroup", 2);
+          "SCWRL4RotamerConstructor::ConstructRRMGroup", 2);
 
   std::vector<RRMRotamerPtr> rotamers;
   Real summed_prob = 0.0;
@@ -1823,13 +1827,13 @@ RRMRotamerGroupPtr SCWRLRotamerConstructor::ConstructRRMGroup(
   return group;
 }
 
-FRMRotamerGroupPtr SCWRLRotamerConstructor::ConstructFRMGroup(
+FRMRotamerGroupPtr SCWRL4RotamerConstructor::ConstructFRMGroup(
           RotamerID id, uint residue_index, 
           std::pair<RotamerLibEntry*,uint> lib_entries, 
           Real probability_cutoff) {
 
   core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
-          "SCWRLRotamerConstructor::ConstructFRMGroup", 2);
+          "SCWRL4RotamerConstructor::ConstructFRMGroup", 2);
 
   std::vector<FRMRotamerPtr> rotamers;
   Real summed_prob = 0.0;
@@ -1913,11 +1917,11 @@ FRMRotamerGroupPtr SCWRLRotamerConstructor::ConstructFRMGroup(
   return group;
 }
 
-RRMRotamerPtr SCWRLRotamerConstructor::ConstructRRM(RotamerID id, 
+RRMRotamerPtr SCWRL4RotamerConstructor::ConstructRRM(RotamerID id, 
                                                     uint residue_idx) {
 
   core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
-          "SCWRLRotamerConstructor::ConstructRRM", 2);
+          "SCWRL4RotamerConstructor::ConstructRRM", 2);
 
   const SCWRLRotamerInfo& info = rotamer_lookup_.GetSidechainInfo(id);
   int num_particles = rotamer_lookup_.GetNumSidechainParticles(id);
@@ -1959,11 +1963,11 @@ RRMRotamerPtr SCWRLRotamerConstructor::ConstructRRM(RotamerID id,
   return rot;  
 }
 
-FRMRotamerPtr SCWRLRotamerConstructor::ConstructFRM(RotamerID id,
+FRMRotamerPtr SCWRL4RotamerConstructor::ConstructFRM(RotamerID id,
                                                     uint residue_idx) {
 
   core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
-          "SCWRLRotamerConstructor::ConstructFRM", 2);
+          "SCWRL4RotamerConstructor::ConstructFRM", 2);
 
   const SCWRLRotamerInfo& info = rotamer_lookup_.GetSidechainInfo(id);
   int num_rrm_particles = rotamer_lookup_.GetNumSidechainParticles(id);
@@ -2005,9 +2009,9 @@ FRMRotamerPtr SCWRLRotamerConstructor::ConstructFRM(RotamerID id,
   subrotamer_definitions.push_back(actual_definition);
 
   const std::vector<SCWRLFRMRule>& frm_rules = rotamer_lookup_.GetFRMRules(id);
-  geom::Vec3 rotation_axis, rotation_anchor;
+  geom::Vec3 rot_axis, rot_anchor;
   geom::Vec3 orig_pos, rot_pos;
-  geom::Mat4 rot;
+  geom::Mat4 transform;
   int base_idx = num_rrm_particles;
 
   for(uint frm_rule_idx = 0; frm_rule_idx < frm_rules.size(); ++frm_rule_idx){
@@ -2015,24 +2019,24 @@ FRMRotamerPtr SCWRLRotamerConstructor::ConstructFRM(RotamerID id,
     int num_fix_particles = frm_rule.fix_particles.size();
     int num_rotating_particles = frm_rule.rotating_particles.size();
 
-    // Update the subrotamer definition... all particles that are fixed are taken
-    // from the first subrotamer.
+    // Update the subrotamer definition... all particles that are fixed are 
+    // taken from the first subrotamer.
     for(int i = 0; i < num_fix_particles; ++i) {
       actual_definition[i] = frm_rule.fix_particles[i];
     }
 
     // The data required for the rotation around the specified axis
-    rotation_anchor = pos_buffer_->GetPos(id, frm_rule.anchor_idx_two);
-    rotation_axis = rotation_anchor - 
-                    pos_buffer_->GetPos(id, frm_rule.anchor_idx_one); 
+    rot_anchor = pos_buffer_->GetPos(id, frm_rule.anchor_idx_two);
+    rot_axis = rot_anchor - pos_buffer_->GetPos(id, frm_rule.anchor_idx_one); 
 
     // Every rotation results in a new subrotamer
-    for(uint prefac_idx = 0; prefac_idx < frm_rule.prefactors.size(); ++prefac_idx){
+    for(uint prefac_idx = 0; prefac_idx < frm_rule.prefactors.size(); 
+        ++prefac_idx){
 
       // Get the rotation matrix
-      Real rotation_angle = frm_rule.prefactors[prefac_idx] * chi_dev_[frm_rule_idx];
-      rot = promod3::core::RotationAroundLine(rotation_axis, rotation_anchor,
-                                              rotation_angle);
+      Real rot_angle = frm_rule.prefactors[prefac_idx] * chi_dev_[frm_rule_idx];
+      transform = promod3::core::RotationAroundLine(rot_axis, rot_anchor,
+                                                    rot_angle);
 
       for(int i = 0; i < num_rotating_particles; ++i){
         int orig_particle_idx = frm_rule.rotating_particles[i];
@@ -2041,57 +2045,9 @@ FRMRotamerPtr SCWRLRotamerConstructor::ConstructFRM(RotamerID id,
         // replace the old with the new index...
         actual_definition[orig_particle_idx] = new_particle_idx;
 
-        // get the new position
-        orig_pos = particles[orig_particle_idx].GetPos();
-        rot_pos[0] = rot(0,0) * orig_pos[0] + rot(0,1) * orig_pos[1] + 
-                     rot(0,2) * orig_pos[2] + rot(0,3);
-        rot_pos[1] = rot(1,0) * orig_pos[0] + rot(1,1) * orig_pos[1] + 
-                     rot(1,2) * orig_pos[2] + rot(1,3);
-        rot_pos[2] = rot(2,0) * orig_pos[0] + rot(2,1) * orig_pos[1] + 
-                     rot(2,2) * orig_pos[2] + rot(2,3);
-
-        particles[new_particle_idx] = 
-        Particle(info.particles[orig_particle_idx].type, rot_pos, 
-                 info.particles[orig_particle_idx].charge, 
-                 info.particles[orig_particle_idx].name);
-      }
-
-      // apply the polar directions and lone pairs
-      // instead of building them from scratch, they get extracted from the
-      // initial particles and transformed
-      for(int i = 0; i < num_rotating_particles; ++i){
-        int orig_particle_idx = frm_rule.rotating_particles[i];
-        int new_particle_idx = base_idx + i;
-
-        if(particles[orig_particle_idx].IsHBondDonor()){
-          orig_pos = particles[orig_particle_idx].GetPos() +
-                     particles[orig_particle_idx].GetPolarDirection();
-          rot_pos[0] = rot(0,0) * orig_pos[0] + rot(0,1) * orig_pos[1] + 
-                       rot(0,2) * orig_pos[2] + rot(0,3);
-          rot_pos[1] = rot(1,0) * orig_pos[0] + rot(1,1) * orig_pos[1] + 
-                       rot(1,2) * orig_pos[2] + rot(1,3);
-          rot_pos[2] = rot(2,0) * orig_pos[0] + rot(2,1) * orig_pos[1] + 
-                       rot(2,2) * orig_pos[2] + rot(2,3);
-          geom::Vec3 polar_direction = rot_pos - 
-                                       particles[new_particle_idx].GetPos();
-          particles[new_particle_idx].SetPolarDirection(polar_direction);
-        }
-
-        if(particles[orig_particle_idx].IsHBondAcceptor()){
-          const std::vector<geom::Vec3>& lp = 
-          particles[orig_particle_idx].GetLonePairs();
-          for(uint j = 0; j < lp.size(); ++j){
-            orig_pos = particles[orig_particle_idx].GetPos() + lp[j];
-            rot_pos[0] = rot(0,0) * orig_pos[0] + rot(0,1) * orig_pos[1] + 
-                         rot(0,2) * orig_pos[2] + rot(0,3);
-            rot_pos[1] = rot(1,0) * orig_pos[0] + rot(1,1) * orig_pos[1] + 
-                         rot(1,2) * orig_pos[2] + rot(1,3);
-            rot_pos[2] = rot(2,0) * orig_pos[0] + rot(2,1) * orig_pos[1] + 
-                         rot(2,2) * orig_pos[2] + rot(2,3);
-            geom::Vec3 new_lp = rot_pos - particles[new_particle_idx].GetPos();
-            particles[new_particle_idx].AddLonePair(new_lp);
-          }
-        }
+        // construct the rotated particle
+        particles[new_particle_idx] = particles[orig_particle_idx];
+        particles[new_particle_idx].ApplyTransform(transform);
       }
 
       subrotamer_definitions.push_back(actual_definition);
@@ -2112,7 +2068,7 @@ FRMRotamerPtr SCWRLRotamerConstructor::ConstructFRM(RotamerID id,
   return r;
 }
 
-void SCWRLRotamerConstructor::MVBBPosBuffer(RotamerID from, RotamerID to) {
+void SCWRL4RotamerConstructor::MVBBPosBuffer(RotamerID from, RotamerID to) {
   pos_buffer_->SetPos(to, promod3::loop::BB_N_INDEX, 
                       pos_buffer_->GetPos(from, promod3::loop::BB_N_INDEX));
   pos_buffer_->SetPos(to, promod3::loop::BB_CA_INDEX, 
@@ -2128,12 +2084,12 @@ void SCWRLRotamerConstructor::MVBBPosBuffer(RotamerID from, RotamerID to) {
                       pos_buffer_->GetPos(from, promod3::loop::BB_CB_INDEX));
 }
 
-void SCWRLRotamerConstructor::ConstructBaseParticles(const SCWRLRotamerInfo& info,
+void SCWRL4RotamerConstructor::ConstructBaseParticles(const SCWRLRotamerInfo& info,
                                                      RotamerID id,
                                                      std::vector<Particle>& particles) {
 
   core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
-          "SCWRLRotamerConstructor::ConstructBaseParticles", 2);
+          "SCWRL4RotamerConstructor::ConstructBaseParticles", 2);
 
   geom::Vec3 a, b, c;
 
@@ -2144,14 +2100,16 @@ void SCWRLRotamerConstructor::ConstructBaseParticles(const SCWRLRotamerInfo& inf
       a = hydrogen_buffer_->GetPos(pinfo.atom_idx);
     }
     else a = pos_buffer_->GetPos(id, pinfo.atom_idx);
-    particles[i] = Particle(pinfo.type, a, pinfo.charge, pinfo.name);
+    particles[i] = Particle(pinfo.name, new SCWRL4Param(pinfo.type, a, 
+                                                        pinfo.charge));
   }
 
   // set the polar directions
   for(uint i = 0; i < info.polar_directions.size(); ++i){
     a = hydrogen_buffer_->GetPos(info.polar_directions[i].polar_idx);
     b = pos_buffer_->GetPos(id, info.polar_directions[i].anchor_idx);
-    particles[info.polar_directions[i].p_idx].SetPolarDirection(a - b);
+    PScoringParam* p = particles[info.polar_directions[i].p_idx].GetSParam();
+    reinterpret_cast<SCWRL4Param*>(p)->SetPolarDirection(a - b);
   }
 
   // set the lone pairs
@@ -2167,17 +2125,18 @@ void SCWRLRotamerConstructor::ConstructBaseParticles(const SCWRLRotamerInfo& inf
     if(lp.C_is_hydrogen) c = hydrogen_buffer_->GetPos(lp.index_C);
     else c = pos_buffer_->GetPos(id, lp.index_C);
     
-    EvalLonePairRule(particles[lp.p_idx], lp.rule, a, b, c);
+    EvalLonePairRule(reinterpret_cast<SCWRL4Param*>(particles[lp.p_idx].GetSParam()), 
+                     lp.rule, a, b, c);
   }
 }
 
-void SCWRLRotamerConstructor::ConstructBBFrameParticles(
+void SCWRL4RotamerConstructor::ConstructBBFrameParticles(
                                  const SCWRLRotamerInfo& info, RotamerID id,
                                  bool n_ter, bool c_ter,
                                  std::vector<Particle>& particles) {
 
   core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
-          "SCWRLRotamerConstructor::ConstructBBFrameParticles", 2);
+          "SCWRL4RotamerConstructor::ConstructBBFrameParticles", 2);
 
   ConstructBaseParticles(info, id, particles);
   if(!(n_ter || c_ter)) return; // we're already done
@@ -2192,9 +2151,11 @@ void SCWRLRotamerConstructor::ConstructBBFrameParticles(
     geom::Vec3 o_pos = pos_buffer_->GetPos(id, promod3::loop::BB_O_INDEX);
     promod3::core::ConstructCTerminalOxygens(c_pos, ca_pos, n_pos,
                                              o_pos, oxt_pos);
-    //manually add oxygen without a rule...
-    particles[p_idx] = Particle(OParticle, oxt_pos, -0.55, "OT");
-    _AddLonePairsCarbonyl(ca_pos, c_pos, oxt_pos, particles[p_idx]);
+    //manually add oxygen without a rule...    
+    particles[p_idx] = Particle("OT", new SCWRL4Param(OParticle, 
+                                                      oxt_pos, -0.55));
+    _AddLonePairsCarbonyl(ca_pos, c_pos, oxt_pos, 
+      reinterpret_cast<SCWRL4Param*>(particles[p_idx].GetSParam()));
     ++p_idx;
   }
 
@@ -2208,10 +2169,15 @@ void SCWRLRotamerConstructor::ConstructBBFrameParticles(
                                     1.0472, ht1_pos);
     promod3::core::ConstructAtomPos(c_pos,ca_pos,n_pos,0.997, 2.042, 
                                     -1.0472, ht2_pos);
-    particles[p_idx] = Particle(HParticle, ht1_pos, 0.35, "HT1");
-    particles[p_idx].SetPolarDirection(ht1_pos - n_pos);
-    particles[p_idx + 1] = Particle(HParticle, ht2_pos, 0.35, "HT2");
-    particles[p_idx + 1].SetPolarDirection(ht2_pos - n_pos);
+    particles[p_idx] = Particle("HT1", new SCWRL4Param(HParticle, 
+                                                       ht1_pos, 0.35));
+    SCWRL4Param* p = reinterpret_cast<SCWRL4Param*>(particles[p_idx].GetSParam());
+    p->SetPolarDirection(ht1_pos - n_pos);
+
+    particles[p_idx + 1] = Particle("HT2", new SCWRL4Param(HParticle, 
+                                                           ht2_pos, 0.35));
+    p = reinterpret_cast<SCWRL4Param*>(particles[p_idx + 1].GetSParam());
+    p->SetPolarDirection(ht2_pos - n_pos);
 
     if(id == PRO || id == TPR || id == CPR) {
       // there are only two hydrogens... please note, that they are not
@@ -2227,8 +2193,10 @@ void SCWRLRotamerConstructor::ConstructBBFrameParticles(
     // and replace it with ht3
     for(uint i = 0; i < info.particles.size(); ++i){
       if(info.particles[i].is_hydrogen){
-        particles[i] = Particle(HParticle, ht3_pos, 0.35, "HT3");
-        particles[i].SetPolarDirection(ht3_pos - n_pos);
+        particles[i] = Particle("HT3", new SCWRL4Param(HParticle, 
+                                                       ht3_pos, 0.35));
+        SCWRL4Param* p = reinterpret_cast<SCWRL4Param*>(particles[i].GetSParam());
+        p->SetPolarDirection(ht3_pos - n_pos);
         break;
       }
     }
diff --git a/sidechain/src/scwrl_rotamer_constructor.hh b/sidechain/src/scwrl4_rotamer_constructor.hh
similarity index 94%
rename from sidechain/src/scwrl_rotamer_constructor.hh
rename to sidechain/src/scwrl4_rotamer_constructor.hh
index 21c31228cf1fad17f10a2d01e98b1b34aaa3edd9..887ea3c39f71690123fd5c81f0e2def69455d225 100644
--- a/sidechain/src/scwrl_rotamer_constructor.hh
+++ b/sidechain/src/scwrl4_rotamer_constructor.hh
@@ -18,11 +18,12 @@
 #define PROMOD3_SCWRL_ROTAMER_CONSTRUCTOR_HH
 
 #include <promod3/sidechain/rotamer_constructor.hh>
+#include <promod3/sidechain/scwrl4_particle_scoring.hh>
 
 namespace promod3 { namespace sidechain {
 
-class SCWRLRotamerConstructor;
-typedef boost::shared_ptr<SCWRLRotamerConstructor> SCWRLRotamerConstructorPtr;
+class SCWRL4RotamerConstructor;
+typedef boost::shared_ptr<SCWRL4RotamerConstructor> SCWRL4RotamerConstructorPtr;
 
 /// \brief Types of lone pair construction
 enum SCWRLLPRule {
@@ -72,11 +73,11 @@ struct SCWRLPDInfo {
 /// \brief Info for particle construction
 struct SCWRLPInfo {
   SCWRLPInfo() { }
-  SCWRLPInfo(SidechainParticle t, Real c, 
+  SCWRLPInfo(SCWRL4ParticleType t, Real c, 
              const String& n, int idx, bool ih) : type(t), charge(c), 
                                                   name(n), atom_idx(idx), 
                                                   is_hydrogen(ih) { }
-  SidechainParticle type;
+  SCWRL4ParticleType type;
   Real charge;
   String name;   // PDB naming
   int atom_idx;  // for AllAtomPositions 
@@ -91,7 +92,7 @@ struct SCWRLPInfo {
 ///        This rule defines how to construct them by using
 ///        three indices of heavy atoms, a bond length and an angle. 
 ///        The dihedral is defined as an index and is set globally in
-///        SCWRLRotamerConstructor. The default hydrogen construction
+///        SCWRL4RotamerConstructor. The default hydrogen construction
 ///        gets overriden for cases defined with such a SCWRLCustomHydrogenInfo
 struct SCWRLCustomHydrogenInfo {
   SCWRLCustomHydrogenInfo() { }
@@ -171,7 +172,7 @@ public:
 
 private:
 
-  void AddInfo(RotamerID id, SidechainParticle p_type, Real charge, 
+  void AddInfo(RotamerID id, SCWRL4ParticleType p_type, Real charge, 
                const String& name, int idx, bool is_h){
     SCWRLPInfo p(p_type, charge, name, idx, is_h);
     sidechain_infos_[id].particles.push_back(p);
@@ -185,7 +186,7 @@ private:
     sidechain_infos_[id].custom_hydrogens.push_back(i);
   }
 
-  void AddBBInfo(RotamerID id, SidechainParticle p_type, Real charge, 
+  void AddBBInfo(RotamerID id, SCWRL4ParticleType p_type, Real charge, 
                  const String& name, int idx, bool is_h){
     SCWRLPInfo p(p_type, charge, name, idx, is_h);
     backbone_infos_[id].particles.push_back(p);
@@ -246,13 +247,13 @@ private:
 };
 
 
-class SCWRLRotamerConstructor : public RotamerConstructor{
+class SCWRL4RotamerConstructor : public RotamerConstructor{
 
 public:
 
-  SCWRLRotamerConstructor(bool cb_in_sidechain = false);
+  SCWRL4RotamerConstructor(bool cb_in_sidechain = false);
 
-  virtual ~SCWRLRotamerConstructor() { }
+  virtual ~SCWRL4RotamerConstructor() { }
 
   // Construct frame residues
   virtual FrameResiduePtr ConstructBackboneFrameResidue(
diff --git a/sidechain/tests/test_frame_construction.cc b/sidechain/tests/test_frame_construction.cc
index 0817985b9a9d2708693957562ac4f196fb4efeac..893b1f130db866c4120b140e17f79e841ead705d 100644
--- a/sidechain/tests/test_frame_construction.cc
+++ b/sidechain/tests/test_frame_construction.cc
@@ -15,7 +15,7 @@
 
 
 
-#include <promod3/sidechain/scwrl_rotamer_constructor.hh>
+#include <promod3/sidechain/scwrl4_rotamer_constructor.hh>
 #include <promod3/core/message.hh>
 #define BOOST_TEST_DYN_LINK
 #include <boost/test/unit_test.hpp>
@@ -64,7 +64,7 @@ void CheckBackboneFrameConstruction(ost::mol::ResidueHandleList& res_list,
   BOOST_CHECK_EQUAL(res_list.size(), all_atoms.GetNumResidues());
   FrameResiduePtr fr1, fr2;
   geom::Vec3 cb_pos;
-  SCWRLRotamerConstructor rc;
+  SCWRL4RotamerConstructor rc;
   for (uint i = 0; i < res_list.size(); ++i) {
     const RotamerID r_id = TLCToRotID(res_list[i].GetName());
     // terminal?
@@ -83,7 +83,7 @@ void CheckBackboneFrameConstruction(ost::mol::ResidueHandleList& res_list,
 void CheckSidechainFrameConstruction(ost::mol::ResidueHandleList& res_list,
                                      loop::AllAtomPositions& all_atoms) {
   // check all frames
-  SCWRLRotamerConstructor rc;
+  SCWRL4RotamerConstructor rc;
   BOOST_CHECK_EQUAL(res_list.size(), all_atoms.GetNumResidues());
   FrameResiduePtr fr1, fr2;
   for (uint i = 0; i < res_list.size(); ++i) {
diff --git a/sidechain/tests/test_rotamers.cc b/sidechain/tests/test_rotamers.cc
index 11143888c2c85c246ee6c46c170560e5cfd045a0..7524af71093b8d86dd8270d475bc94b4e13a3543 100644
--- a/sidechain/tests/test_rotamers.cc
+++ b/sidechain/tests/test_rotamers.cc
@@ -15,7 +15,7 @@
 
 
 
-#include <promod3/sidechain/scwrl_rotamer_constructor.hh>
+#include <promod3/sidechain/scwrl4_rotamer_constructor.hh>
 #include <promod3/core/message.hh>
 #define BOOST_TEST_DYN_LINK
 #include <boost/test/unit_test.hpp>
@@ -114,7 +114,7 @@ void CheckRotamerGroupConstruction(ost::mol::ResidueHandleList& res_list,
   BOOST_CHECK_EQUAL(res_list.size(), all_atoms.GetNumResidues());
   RRMRotamerGroupPtr rrm1, rrm2;
   FRMRotamerGroupPtr frm1, frm2;
-  SCWRLRotamerConstructor rot_constructor;
+  SCWRL4RotamerConstructor rot_constructor;
   for (uint i = 0; i < res_list.size(); ++i) {
     const RotamerID r_id = TLCToRotID(res_list[i].GetName());
     // terminal?
@@ -170,7 +170,7 @@ void CheckRotamerApplyOnResidue(ost::mol::ResidueHandleList& res_list,
   BOOST_CHECK_EQUAL(res_list.size(), all_atoms.GetNumResidues());
   RRMRotamerGroupPtr rrm;
   FRMRotamerGroupPtr frm;
-  SCWRLRotamerConstructor rc;
+  SCWRL4RotamerConstructor rc;
   for (uint i = 0; i < res_list.size(); ++i) {
     // do it
     const RotamerID r_id = TLCToRotID(res_list[i].GetName());