diff --git a/gromos++/programs/0index.doxy b/gromos++/programs/0index.doxy
index 1ad263ce..7193f1df 100644
--- a/gromos++/programs/0index.doxy
+++ b/gromos++/programs/0index.doxy
@@ -44,6 +44,7 @@
  * - Drazen Petrov @anchor dp (dp)
  * - Anita de Ruiter @anchor adr (adr)
  * - Martina Setz @anchor ms (ms)
+ * - Stefan Bachmann @anchor sb (sb)
  * </p>
  *
  * <h2>Available Programs:</h2>
@@ -60,6 +61,8 @@
  * - @ref com_top combine molecular topology files into one
  * - @ref con_top convert topology to different force field version
  * - @ref copy_box repeats a simulation box in a given direction
+ * - @ref cos_dipole calculate molecular dipoles including COS charges
+ * - @ref cos_epsilon calculate relative permittivity and box dipole moment autocorrelations 
  * - @ref cry perform symmetry operations on molecules
  * - @ref cry_rms perform RMSD and RMSF calculations based on crystallographic symmetry.
  * - @ref dfmult calculates multiple free energy differences
diff --git a/gromos++/programs/Makefile.am b/gromos++/programs/Makefile.am
index fd74cd02..18b9318f 100644
--- a/gromos++/programs/Makefile.am
+++ b/gromos++/programs/Makefile.am
@@ -105,7 +105,9 @@ bin_PROGRAMS = rmsd\
     trs_ana\
     dGslv_pbsolv\
     dfgrid\
-    contactnum
+    contactnum\
+    cos_dipole\
+    cos_epsilon
 
 nhoparam_SOURCES = nhoparam.cc
 frameout_SOURCES = frameout.cc
@@ -207,6 +209,9 @@ trs_ana_SOURCES = trs_ana.cc
 dGslv_pbsolv_SOURCES = dGslv_pbsolv.cc
 dfgrid_SOURCES = dfgrid.cc
 contactnum_SOURCES = contactnum.cc
+cos_dipole_SOURCES = cos_dipole.cc
+cos_epsilon_SOURCES = cos_epsilon.cc
+
 
 LDADD = $(top_builddir)/src/libgromos.la
 
diff --git a/gromos++/programs/atominfo.cc b/gromos++/programs/atominfo.cc
index 6bfe9673..8ee2ce3d 100644
--- a/gromos++/programs/atominfo.cc
+++ b/gromos++/programs/atominfo.cc
@@ -155,7 +155,7 @@ int main(int argc, char **argv){
 	 << setw(12) << "Atom Code"
 	 << endl;
     
-    for(int i=0; i < as.size(); ++i){
+    for(unsigned int i=0; i < as.size(); ++i){
       if(as.atom()[i]->type()==utils::spec_virtual){
 	utils::AtomSpecifier conf=as.atom()[i]->conf();
 	cout << "----------------------------------------"
@@ -202,7 +202,7 @@ int main(int argc, char **argv){
 	    
 	}
 	
-	for(int j=0; j< conf.size(); ++j){
+	for(unsigned int j=0; j< conf.size(); ++j){
 	  cout << setw(13) << conf.toString(j)
 	       << setw(10) << conf.gromosAtom(j)+1
 	       << setw(10) << conf.resnum(j)+1
diff --git a/gromos++/programs/bilayer_dist.cc b/gromos++/programs/bilayer_dist.cc
index c52b7d14..1e133541 100644
--- a/gromos++/programs/bilayer_dist.cc
+++ b/gromos++/programs/bilayer_dist.cc
@@ -197,7 +197,7 @@ int main(int argc, char** argv) {
         min_d = -sys.box().M().abs()/2; // before was: cm[2]-sys.box()[2]/2;
         max_d = sys.box().M().abs()/2;
         
-        for(int i = 0; i < selection.size(); i++) {
+        for(unsigned int i = 0; i < selection.size(); i++) {
           Vec & vector1 = selection.pos(i);
           Vec dist = vector1 - pbc->nearestImage(vector1, cm, sys.box());
           zcoord.push_back(dist[2]);
diff --git a/gromos++/programs/check_box.cc b/gromos++/programs/check_box.cc
index e105c762..a85fce96 100644
--- a/gromos++/programs/check_box.cc
+++ b/gromos++/programs/check_box.cc
@@ -149,7 +149,7 @@ int main(int argc, char **argv){
 
     //init variables
     double overall_min_dist2=1e4;
-    int limit = 1;
+    unsigned int limit = 1;
     double cutoff = 1.4;
     double cutoff2 = cutoff*cutoff;
     bool no_cutoff=false; //do not report all distances
@@ -359,12 +359,12 @@ int main(int argc, char **argv){
                     //if 2 neighbours have the same atoms (=there is only 1 cube total in this direction), only half of the distances must be calculated:
                     bool skip=cubes.same_cube(c);
 
-                    for(int i=0; i < cubes.cube_i_atomlist(c).size(); ++i ){
+                    for(unsigned int i=0; i < cubes.cube_i_atomlist(c).size(); ++i ){
 
                         int atom_i = cubes.cube_i_atomlist(c)[i];
                         Vec atom_i_pos = atoms.pos(atom_i);
 
-                        int j=0;
+                        unsigned int j=0;
                         if(skip)
                             j=i+1;
 
@@ -438,7 +438,7 @@ int main(int argc, char **argv){
     cout << "###" << endl;
     cout << "# Time     Distance [nm]     Atom A - Atom B      Mol:Residue Atom A -  Mol:Resiude Atom B" << endl;
     double itime=0;
-    for (int ii =0; ii < output.size(); ii++) {
+    for (unsigned int ii =0; ii < output.size(); ii++) {
       for(vec_it = output[ii].begin(); vec_it != output[ii].end(); ++vec_it){
 
         if (read_time) itime=vec_it->time;
diff --git a/gromos++/programs/close_pair.cc b/gromos++/programs/close_pair.cc
index d4e76451..1e74c60a 100644
--- a/gromos++/programs/close_pair.cc
+++ b/gromos++/programs/close_pair.cc
@@ -148,10 +148,10 @@ int main(int argc, char **argv) {
       debug = atoi(args["debug"].c_str());
     }
     if (debug == 1) {
-      for (int i = 0; i < groupA.size(); ++i) {
+      for (unsigned int i = 0; i < groupA.size(); ++i) {
         cout << "# groupA mol " << groupA.mol(i) << " size " << groupA.size() << endl;
       }
-      for (int j = 0; j < groupB.size(); ++j) {
+      for (unsigned int j = 0; j < groupB.size(); ++j) {
         cout << "# groupB mol " << groupB.mol(j) << " size " << groupB.size() << endl;
       }
     }
@@ -224,7 +224,7 @@ int main(int argc, char **argv) {
 
         int pass = 1000;
 
-        for (int i = 0; i < groupA.size(); ++i) {
+        for (unsigned int i = 0; i < groupA.size(); ++i) {
           // skip hydrogen atoms
           int m = groupA.mol(i);
 
@@ -236,7 +236,7 @@ int main(int argc, char **argv) {
                   "atom name not specified properly.");
           if (m != pass && m > 0)
             if (groupA.mass(i) != 1.008) {
-              for (int j = 0; j < groupB.size(); ++j) {
+              for (unsigned int j = 0; j < groupB.size(); ++j) {
                 if (debug >= 3) {
                   cout << "# groupB size " << groupB.size() << endl;
                   cout << "# molA " << groupA.mol(i)
diff --git a/gromos++/programs/cog.cc b/gromos++/programs/cog.cc
index 04785cef..7ff25a96 100644
--- a/gromos++/programs/cog.cc
+++ b/gromos++/programs/cog.cc
@@ -196,7 +196,7 @@ int main(int argc, char **argv){
       }
     }
     
-    for(int a=0; a<atomspec.size(); a++){
+    for(unsigned int a=0; a<atomspec.size(); a++){
 	  osys.addMolecule(mt_spec);
 	  osys.mol(molnum).initPos();
 	  molnum++;
@@ -266,10 +266,10 @@ int main(int argc, char **argv){
 	  }
 	  }
 	  
-	  for (int i=0; i < atomspec.size(); i++) {
+	  for (unsigned int i=0; i < atomspec.size(); i++) {
 	    Vec cog(0.0,0.0,0.0);
 	    double molmass=0;
-	    for (int a=0; a < atomspec[i].size(); a++) {
+	    for (unsigned int a=0; a < atomspec[i].size(); a++) {
 	      if(inc=="cog"){
 		cog+=atomspec[i].pos(a);
 	      }
diff --git a/gromos++/programs/contactnum.cc b/gromos++/programs/contactnum.cc
index c39bf75f..eabb57bb 100644
--- a/gromos++/programs/contactnum.cc
+++ b/gromos++/programs/contactnum.cc
@@ -91,7 +91,7 @@ double fastpow(double base, int exp);
 int main(int argc, char **argv){
   Argument_List knowns; 
   knowns << "topo" << "traj" << "atoms" << "cutoffs" << "pbc" << "exp"
-         << "time" << "debug" << "printpairs" << "excludemass";
+         << "time" << "printpairs" << "excludemass";
 
   string usage = "# " + string(argv[0]);
   usage += "\n\t@topo         <molecular topology file>\n";
@@ -120,10 +120,6 @@ int main(int argc, char **argv){
     InTopology it(args["topo"]);
     System sys(it.system());
     
-    int debug=0;
-    if(args.count("debug")>0)
-        debug=1;
-        
     if (args.count("pbc") <=0) {
       cout << "# WARNING: no @pbc argument found - no gathering will be done!";
     }
@@ -189,7 +185,7 @@ int main(int argc, char **argv){
     
     // exclude atoms with given mass
     for (int i=atomgroup1.size()-1; i>=0 ; i--) {
-      for(int j=0; j<excludemass.size(); j++) {
+      for(unsigned int j=0; j<excludemass.size(); j++) {
         if (atomgroup1.mass(i)== excludemass[j]) {
           //cerr << "atom " << atomgroup1.atom(i) << " mass " << atomgroup1.mass(i)<<endl;
           atomgroup1.removeAtom(i);
@@ -197,7 +193,7 @@ int main(int argc, char **argv){
       }
     }
     for (int i=atomgroup2.size()-1; i>=0 ; i--) {
-      for(int j=0; j<excludemass.size(); j++) {
+      for(unsigned int j=0; j<excludemass.size(); j++) {
         if (atomgroup2.mass(i)== excludemass[j]){
           //cerr << "atom " << atomgroup2.atom(i) << " mass " << atomgroup2.mass(i) << endl;
           atomgroup2.removeAtom(i);
@@ -213,7 +209,7 @@ int main(int argc, char **argv){
     // Parse boundary conditions for sys
     Boundary *pbc = BoundaryParser::boundary(sys, args);
     
-	for (int cut=0; cut<cutoffs.size(); cut++) {
+	for (unsigned int cut=0; cut<cutoffs.size(); cut++) {
 	   cout <<setw(15) << cutoffs[cut];
 	}
 	cout << endl;
@@ -240,10 +236,10 @@ int main(int argc, char **argv){
 	
         cout.precision(2);
         cout << time;
-	    for (int cut=0; cut<cutoffs.size(); cut++) {
+	    for (unsigned int cut=0; cut<cutoffs.size(); cut++) {
 	      double ct=0;
-          for (int i=0; i< atomgroup1.size(); i++) {
-            for (int j=0; j < atomgroup2.size(); j++) {
+          for (unsigned int i=0; i< atomgroup1.size(); i++) {
+            for (unsigned int j=0; j < atomgroup2.size(); j++) {
               double func;
               gmath::Vec v = atomgroup1.pos(i) - pbc->nearestImage(atomgroup1.pos(i), atomgroup2.pos(j), sys.box());
               double rdist=v.abs();
@@ -263,7 +259,7 @@ int main(int argc, char **argv){
       ic.close();
     }  
     cout <<"# averages     ";
-    for (int i=0; i<cutoffs.size(); i++) {
+    for (unsigned int i=0; i<cutoffs.size(); i++) {
       cout << setw(15) << ct_ave[i]/numFrames; 
     }
     cout << endl;
diff --git a/gromos++/programs/cos_dipole.cc b/gromos++/programs/cos_dipole.cc
new file mode 100644
index 00000000..bd3ecb4f
--- /dev/null
+++ b/gromos++/programs/cos_dipole.cc
@@ -0,0 +1,534 @@
+/**
+ * @file cos_dipole.cc
+ * Calculate molecular dipoles including COS charges
+ */
+
+/**
+ * @page programs Program Documentation
+ *
+ * @anchor cos_dipole
+ * @section cos_dipole Calculate molecular dipoles including COS charges
+ * @author @ref sb @ref mp
+ * @date 11-01-2018
+ *
+ * Program cos_dipole calculates the average dipole moments over a selected set of molecules, taking into account also polarizable sites.
+ * Standardly it outputs the magnitude of the average total, fixed and induced molecular dipoles, but if needed, the x-, y- and z-components can be written to a file by specifying the \@xyz flag.
+ * 
+ * Note that the total dipole moment is only well-defined for 
+ * systems consisting of neutral molecules. If the system carries a net-charge,
+ * the dipole moment will depend on the position of the origin. In cases where
+ * the overall system is neutral but contains ions, the dipole moment will
+ * depend on which periodic copy of the ions is taken. In these cases, the program
+ * issues a warning that results will critically depend on the choice of
+ * gathering method.
+ *
+ * <b>arguments:</b>
+ * <table border=0 cellpadding=0>
+ * <tr><td> \@topo</td><td>&lt;molecular topology file&gt; </td></tr>
+ * <tr><td> \@pbc</td><td>&lt;boundary type&gt; [&lt;gather method&gt;] </td></tr>
+ * <tr><td> [\@time</td><td>&lt;@ref utils::Time "time and dt"&gt;] </td></tr>
+ * <tr><td> [\@fac</td><td>&lt;conversion factor for the unit of the dipole, default: 1; use 48.032045 to convert from e*nm to Debye&gt;] </td></tr>
+ * <tr><td> [\@molecules</td><td>&lt;solute molecules to average over, e.g. 1-5,37,100-101&gt;] </td></tr>
+ * <tr><td> [\@solv</td><td>&lt;include solvent&gt;]</td></tr>
+ * <tr><td> [\@xyz</td><td>&lt;filename for writing out dipole x-,y-,z-components, Default: Mxyz.out&gt;]</td></tr>
+ * <tr><td> \@traj</td><td>&lt;trajectory files&gt; </td></tr>
+ * <tr><td> \@trs</td><td>&lt;special trajectory files with COS displacements&gt; </td></tr>
+ * </table>
+ *
+ *
+ * Example:
+ * @verbatim
+  cos_dipole
+    @topo  ex.top
+    @pbc   r
+    @fac 48.032045
+    @molecules 1-5,8-10
+    @solv
+    [@time  0 0.2]
+    @traj  ex.trc
+    @trs   ex.trs
+ @endverbatim
+ *
+ * <hr>
+ */
+
+#include <vector>
+#include <iomanip>
+#include <iostream>
+#include <fstream>
+#include <cassert>
+#include <locale>
+
+#include "../src/args/Arguments.h"
+#include "../src/args/BoundaryParser.h"
+#include "../src/args/GatherParser.h"
+#include "../src/fit/Reference.h"
+#include "../src/gio/InG96.h"
+#include "../src/gcore/System.h"
+#include "../src/gcore/Molecule.h"
+#include "../src/gcore/Box.h"
+#include "../src/gio/InTopology.h"
+#include "../src/bound/Boundary.h"
+#include "../src/fit/PositionUtils.h"
+#include "../src/gmath/Vec.h"
+#include "../src/gmath/Correlation.h"
+#include "../src/gmath/Physics.h"
+#include "../src/utils/AtomSpecifier.h"
+#include "../src/utils/groTime.h"
+#include "../src/utils/parse.h"
+
+using namespace std;
+using namespace fit;
+using namespace gcore;
+using namespace gio;
+using namespace bound;
+using namespace args;
+
+void get_atom_dipole(const AtomTopology &current_atom, Vec &atom_pos,
+                     Vec &offsite_i, Vec &offsite_j, Vec &cosdispl,
+                     int &pol_count,
+                     Vec &mol_dip, Vec &fix_dip, Vec &ind_dip,
+                     double &nchargepa)
+{
+  Vec pol_site_pos(0, 0, 0);
+  double pol_site_q = 0;
+  if (current_atom.isPolarisable())
+  {
+    pol_count++;
+    if (current_atom.poloffsiteGamma() > 0.000001)
+    {
+      double off_set_gamma;
+      off_set_gamma = current_atom.poloffsiteGamma();
+
+      //position of the offset atom
+      pol_site_pos = atom_pos + off_set_gamma *
+                                    (offsite_i + offsite_j - 2 * atom_pos);
+    }
+    else
+    {
+      pol_site_pos = atom_pos;
+    }
+    pol_site_q = current_atom.charge() + (current_atom.cosCharge() * -1);
+    Vec cos_pos = pol_site_pos + cosdispl;
+
+    // contribution of the cos charges with cos position
+    mol_dip += cos_pos * current_atom.cosCharge();
+    ind_dip += cosdispl * current_atom.cosCharge();
+  }
+  else
+  {
+    // case no polarisation nor offset
+    pol_site_pos = atom_pos;
+    pol_site_q = current_atom.charge();
+  }
+  mol_dip += (pol_site_q - nchargepa) * pol_site_pos;
+  fix_dip += (current_atom.charge() - nchargepa) * pol_site_pos;
+}
+
+int main(int argc, char **argv)
+{
+
+  Argument_List knowns;
+  knowns << "topo"
+         << "pbc"
+         << "time"
+         << "molecules"
+         << "fac"
+         << "solv"
+         << "xyz"
+         << "traj"
+         << "trs";
+
+  string usage = "# " + string(argv[0]);
+  usage += "\n\t@topo <molecular topology file>\n";
+  usage += "\t@pbc    <boundary type> [<gather method>]\n";
+  usage += "\t[@time   <time and dt>]\n";
+  usage += "\t[@fac   <conversion factor for the unit of the dipole, default: 1; use 48.032045 to convert from e*nm to Debye>]\n";
+  usage += "\t[@molecules <solute molecules to average over, e.g. 1-5,37,100-101>]\n";
+  usage += "\t[@solv <include solvent>]\n";
+  usage += "\t[@xyz <filename for writing out dipole x-,y-,z-components, Default: Mxyz.out>]\n";
+  usage += "\t@traj   <trajectory files>\n";
+  usage += "\t@trs    <special traj with cosdisplacement>\n";
+
+  try
+  {
+    Arguments args(argc, argv, knowns, usage);
+
+    // get the @time argument
+    utils::Time time(args), time_trs(args);
+
+    //  read topology
+    InTopology it(args["topo"]);
+    System sys(it.system());
+
+    System refSys(it.system());
+
+    // include solvent
+    bool include_solvent = false;
+    if (args.count("solv") > -1)
+      include_solvent = true;
+
+    // write x-, y- z-components?
+    std::string fname_xyz;
+    bool write_xyz = false;
+    if (args.count("xyz") == 0)
+    {
+      fname_xyz = "Mxyz.out";
+      write_xyz = true;
+    }
+    else if (args.count("xyz") > 0)
+    {
+      fname_xyz = args.getValue<std::string>("xyz", false, "Mxyz.out");
+      write_xyz = true;
+    }
+
+    // create an atomspecifier that contains all atoms (including the solvent)
+    // or just the ones that got specified by @molecules
+    utils::AtomSpecifier atoms(sys);
+
+    for (int m = 0; m < sys.numMolecules(); ++m)
+      for (int a = 0; a < sys.mol(m).numAtoms(); ++a)
+        atoms.addAtom(m, a);
+
+    //determine net charge
+    double ncharge = 0;
+    double nchargepa = 0;
+    for (unsigned int i = 0; i < atoms.size(); i++)
+    {
+      ncharge += atoms.charge(i);
+    }
+    // atoms.size() does not include the solvent, so for solvent-only systems:
+    int num_solv_mol = 0;
+    bool charged_solvent = false;
+    if (include_solvent)
+    {
+      // read first frame to get number of solvent atoms
+      InG96 ic;
+      if (args.count("traj") > 0)
+      {
+        ic.open(args.lower_bound("traj")->second);
+
+        ic.select("ALL");
+        ic >> sys;
+        ic.close();
+      }
+      else
+      {
+        throw gromos::Exception("dipole", "no trajectory specified (@traj)");
+      }
+      atoms.addSpecifier("s:a");
+      for (int i = 0; i < sys.numSolvents(); i++)
+      {
+        double solv_charge_mol = 0;
+        for (int a = 0; a < sys.sol(i).topology().numAtoms(); a++)
+          solv_charge_mol += sys.sol(i).topology().atom(a).charge();
+        if (solv_charge_mol)
+          charged_solvent = true;
+
+        int _num_solv_mol = sys.sol(i).numAtoms() / sys.sol(i).topology().numAtoms();
+        num_solv_mol += _num_solv_mol;
+        ncharge += num_solv_mol * solv_charge_mol;
+      }
+    }
+
+    if (atoms.size() + atoms.numSolventAtoms())
+      nchargepa = ncharge / (atoms.size() + atoms.numSolventAtoms());
+    else
+      nchargepa = 0;
+
+    int num_mol = 0;
+    std::vector<int> molecules;
+    std::string mol_parse_string;
+    if (args.count("molecules") == -1)
+    {
+      mol_parse_string = "all";
+      for (int i = 0; i < sys.numMolecules(); i++)
+      {
+        molecules.push_back(i);
+      }
+      num_mol = sys.numMolecules();
+      cout << "# average dipole and quadrupole moments over all " << num_mol << " solute molecules" << endl;
+    }
+    else
+    {
+
+      mol_parse_string = args.getValue<std::string>("molecules", false, "1-2");
+      utils::parse_range(mol_parse_string, molecules);
+      num_mol = molecules.size();
+
+      for (unsigned int i = 0; i < molecules.size(); i++)
+      {
+        if (molecules[i] > sys.numMolecules())
+        {
+          std::ostringstream oss;
+          oss << "flag @molecules: " << mol_parse_string << " , but there are only " << sys.numMolecules() << " molecules in the system";
+          throw gromos::Exception("cos_dipole", oss.str());
+        }
+      }
+      cout << "# average dipole and quadrupole moments of molecules " << mol_parse_string << endl;
+    }
+
+    if (include_solvent)
+      cout << "# including " << num_solv_mol << " solvent molecules " << std::endl;
+
+    // read conversion factor
+    double conversion_fac = args.getValue<double>("fac", false, 1);
+
+    // parse boundary conditions
+    Boundary *pbc = BoundaryParser::boundary(sys, args);
+    // parse gather method
+    Boundary::MemPtr gathmethod = args::GatherParser::parse(sys, refSys, args);
+
+    // write a title
+    cout << "#\n";
+    if (ncharge != 0.0)
+    {
+      cout << "# WARNING the system carries a net charge ( "
+           << ncharge << ")\n"
+           << "#         this means that the dipole depends on the origin!\n"
+           << "#\n";
+    }
+    else
+    {
+      // we also want to know if there are any ions
+      int ion_count = 0;
+      for (int m = 0; m < sys.numMolecules(); ++m)
+      {
+        double tc = 0.0;
+        for (int a = 0; a < sys.mol(m).numAtoms(); ++a)
+          tc += sys.mol(m).topology().atom(a).charge();
+        if (tc != 0.0)
+          ion_count++;
+      }
+      if (ion_count != 0 || charged_solvent)
+      {
+        cout << "# WARNING although the system is neutral overall,\n";
+        if (ion_count != 0)
+          cout << "#        - there are" << ion_count << " molecules that carry a charge\n";
+        if (charged_solvent)
+          cout << "#        - there is charged solvent\n";
+        cout << "#\n"
+             << "# -> the system-dipole will depend on their positions in "
+             << "the periodic box\n"
+             << "#         this is likely to give random results\n";
+      }
+    }
+
+    cout << "#" << setw(14) << " time"
+         << setw(15) << " dipole_tot"
+         << setw(15) << " dipole_fix"
+         << setw(15) << " dipole_induced"
+         << setw(15) << " Q_xx"
+         << setw(15) << " Q_yy"
+         << setw(15) << " Q_zz\n";
+
+    ofstream os(fname_xyz.c_str());
+    if (write_xyz)
+    {
+      os << "# x-, y- and z- components of the total, fixed and reduced dipole moments\n";
+      os << "# averages over the selected solute molecules: " << mol_parse_string << "\n";
+      if (include_solvent)
+        os << "# including " << num_solv_mol << " solvent molecules " << std::endl;
+      os << "#\n";
+      os << "#" << setw(9) << "Time"
+         << setw(12) << "Mtot_x" << setw(12) << "Mtot_y" << setw(12) << "Mtot_z"
+         << setw(12) << "Mfix_x" << setw(12) << "Mfix_y" << setw(12) << "Mfix_z"
+         << setw(12) << "Mind_x" << setw(12) << "Mind_y" << setw(12) << "Mind_z"
+         << std::endl;
+    }
+
+    int numFrames = 0;
+
+    // sum over the timeseries
+    double tot_sum_mol_dip = 0, tot_sum_fix_dip = 0, tot_sum_ind_dip_frame = 0;
+
+    vector<Vec> dipole_vector;
+
+    // define input coordinate
+    InG96 ic, is;
+
+    // loop over all trajectories
+    for (Arguments::const_iterator
+             iter = args.lower_bound("traj"),
+             to = args.upper_bound("traj"),
+             siter = args.lower_bound("trs"),
+             sto = args.upper_bound("trs");
+         iter != to || siter != sto; ++iter, ++siter)
+    {
+
+      // open file
+      ic.open((iter->second).c_str());
+      ic.select("ALL");
+
+      //open special traj
+      is.open((siter->second).c_str());
+      is.select("ALL");
+
+      // loop over single trajectory
+      while (!ic.eof() && !is.eof())
+      {
+        is >> sys >> time_trs;
+
+        ic >> sys >> time;
+
+        // check if the times are the same
+        if (time_trs.time() != time.time())
+        {
+          cerr << "Time is not corresponding between trajectories\n"
+               << "Special traj:\t" << time_trs.time() << "\n"
+               << "Position traj:\t" << time.time() << endl;
+          return -1;
+        }
+
+        // we have to reconnect the molecules
+        // in the case of neutral molecules, this will be sufficient
+        // non-neutral molecules will behave weirdly.
+        (*pbc.*gathmethod)();
+
+        Vec tot_dipole(0.0, 0.0, 0.0), fix_dipole(0.0, 0.0, 0.0), ind_dipole(0.0, 0.0, 0.0);
+        double sum_mol_dip_frame = 0;
+        double sum_fix_dip_frame = 0;
+        double sum_ind_dip_frame = 0;
+
+        int pol_count = 0;
+
+//calculate molecular dipoles for selected solute molecules
+        for (unsigned int m = 0; m < molecules.size(); m++)
+        {
+          Molecule molecule = sys.mol(molecules[m] - 1);
+
+          Vec mol_com(0, 0, 0);
+          double mass = 0;
+          for (int a = 0; a < molecule.numAtoms(); a++)
+          {
+            mol_com += molecule.topology().atom(a).mass() * molecule.pos(a);
+            mass += molecule.topology().atom(a).mass();
+          }
+          mol_com /= mass;
+          //cerr << "mass " << mass << " com " <<  v2s(mol_com) << endl;
+
+          //cout << "id" << omp_get_thread_num() << " mol" << m << " " <<std::endl;
+          Vec mol_dip(0, 0, 0);
+          Vec fix_dip(0, 0, 0);
+          Vec ind_dip(0, 0, 0);
+
+          Vec offsite_i(0, 0, 0), offsite_j(0, 0, 0);
+
+          //cerr << "#Amount of atoms in molecule:\t" << sys.mol(m).numAtoms() << endl;
+          for (int a = 0; a < molecule.numAtoms(); a++)
+          {
+            const AtomTopology &current_atom = molecule.topology().atom(a);
+
+            Vec atom_pos = molecule.pos(a) - mol_com;
+            if (current_atom.isPolarisable() && current_atom.poloffsiteGamma() > 0.000001)
+            {
+              offsite_i = molecule.pos(current_atom.poloffsiteI()) - mol_com;
+              offsite_j = molecule.pos(current_atom.poloffsiteJ()) - mol_com;
+            }
+            Vec cosdispl = molecule.cosDisplacement(a);
+
+            get_atom_dipole(current_atom, atom_pos, offsite_i, offsite_j, cosdispl, pol_count, mol_dip, fix_dip, ind_dip, nchargepa);
+          }
+          //cerr << m << ":\t" << mol_dip.abs() << endl;
+          tot_dipole += mol_dip;
+          fix_dipole += fix_dip;
+          ind_dipole += ind_dip;
+          sum_mol_dip_frame += mol_dip.abs();
+          sum_fix_dip_frame += fix_dip.abs();
+          sum_ind_dip_frame += ind_dip.abs();
+        }
+
+        // calculate molecular dipoles for solvent molecules
+        if (include_solvent)
+        {
+          for (int s = 0; s < sys.numSolvents(); s++)
+          {
+            int num_solvent_atoms = sys.sol(s).topology().numAtoms();
+            for (int i = 0; i < sys.sol(s).numAtoms(); i += num_solvent_atoms)
+            {
+              Vec mol_dip(0, 0, 0);
+              Vec fix_dip(0, 0, 0);
+              Vec ind_dip(0, 0, 0);
+
+              Vec offsite_i(0, 0, 0), offsite_j(0, 0, 0);
+
+              //cerr << "#Amount of atoms in molecule:\t" << sys.mol(m).numAtoms() << endl;
+              for (int a = 0; a < num_solvent_atoms; a++)
+              {
+                const AtomTopology &current_atom = sys.sol(s).topology().atom(a);
+
+                Vec atom_pos = sys.sol(s).pos(i + a);
+
+                if (current_atom.isPolarisable() && current_atom.poloffsiteGamma() > 0.000001)
+                {
+                  offsite_i = sys.sol(s).pos(i + a + current_atom.poloffsiteI());
+                  offsite_j = sys.sol(s).pos(i + a + current_atom.poloffsiteJ());
+                }
+                Vec cosdispl = sys.sol(s).cosDisplacement(i + a);
+
+                get_atom_dipole(current_atom, atom_pos, offsite_i, offsite_j, cosdispl, pol_count, mol_dip, fix_dip, ind_dip, nchargepa);
+              }
+              //cerr << m << ":\t" << mol_dip.abs() << endl;
+              tot_dipole += mol_dip;
+              fix_dipole += fix_dip;
+              ind_dipole += ind_dip;
+              sum_mol_dip_frame += mol_dip.abs();
+              sum_fix_dip_frame += fix_dip.abs();
+              sum_ind_dip_frame += ind_dip.abs();
+            }
+          }
+        }
+
+        int tot_num_mol = num_mol + num_solv_mol;
+        tot_sum_mol_dip += sum_mol_dip_frame / tot_num_mol;
+        tot_sum_fix_dip += sum_fix_dip_frame / tot_num_mol;
+
+        double ave_ind_dip_frame;
+        if (pol_count != 0)
+        {
+          tot_sum_ind_dip_frame += sum_ind_dip_frame / pol_count;
+          ave_ind_dip_frame = sum_ind_dip_frame * conversion_fac / pol_count;
+        }
+        else
+        {
+          ave_ind_dip_frame = 0;
+        }
+        // do some bookkeeping
+        numFrames++;
+
+        cout << time
+             << setw(15) << setprecision(8) << sum_mol_dip_frame * conversion_fac / tot_num_mol
+             << setw(15) << setprecision(8) << sum_fix_dip_frame * conversion_fac / tot_num_mol
+             << setw(15) << setprecision(8) << ave_ind_dip_frame
+             << endl;
+        if (write_xyz)
+        {
+          os << setw(10) << setprecision(5) << time.time()
+             << setw(12) << tot_dipole[0] * conversion_fac / tot_num_mol << setw(12) << tot_dipole[1] * conversion_fac / tot_num_mol << setw(12) << tot_dipole[2] * conversion_fac / tot_num_mol
+             << setw(12) << fix_dipole[0] * conversion_fac / tot_num_mol << setw(12) << fix_dipole[1] * conversion_fac / tot_num_mol << setw(12) << fix_dipole[2] * conversion_fac / tot_num_mol
+             << setw(12) << ind_dipole[0] * conversion_fac / tot_num_mol << setw(12) << ind_dipole[1] * conversion_fac / tot_num_mol << setw(12) << ind_dipole[2] * conversion_fac / tot_num_mol
+             << std::endl;
+        }
+      }
+      ic.close();
+      is.close();
+    }
+
+    delete pbc;
+    os.close();
+
+    //now print average over the systems for all three values
+    cout << "#\n#" << setw(14) << "Averages:"
+         << setw(15) << setprecision(8) << tot_sum_mol_dip * conversion_fac / numFrames
+         << setw(15) << setprecision(8) << tot_sum_fix_dip * conversion_fac / numFrames
+         << setw(15) << setprecision(8) << tot_sum_ind_dip_frame * conversion_fac / numFrames
+         << endl;
+  }
+
+  catch (const gromos::Exception &e)
+  {
+    cerr << e.what() << endl;
+    exit(1);
+  }
+  return 0;
+}
diff --git a/gromos++/programs/cos_epsilon.cc b/gromos++/programs/cos_epsilon.cc
new file mode 100644
index 00000000..655fcfd9
--- /dev/null
+++ b/gromos++/programs/cos_epsilon.cc
@@ -0,0 +1,508 @@
+/**
+ * @file cos_epsilon.cc
+ * Calculate relative permittivity and box dipole moment autocorrelations
+ */
+
+/**
+ * @page programs Program Documentation
+ *
+ * @anchor cos_epsilon
+ * @section cos_epsilon Calculate relative permittivity and box dipole moment autocorrelations
+ * @author @ref sb @ref mp
+ * @date 11-01-2018
+ *
+ * Program cos\_epsilon estimates the relative dielectric permittivity, 
+ * @f$\epsilon(0)@f$, of a simulation box from a Kirkwood - Fr&ouml;hlich type of
+ * equation, as derived by Neumann [Mol. Phys. 50, 841 (1983)], taking into account also polarizable centers.
+ *
+ * @f[ (\epsilon(0) - 1) \frac{2\epsilon_{RF}+1}{2\epsilon_{RF}+\epsilon(0)} = \frac{<\vec{M}^2> - <\vec{M}>^2}{3\epsilon_0 Vk_BT} @f]
+ *
+ * where @f$\vec{M}@f$ is the total dipole moment of the system, 
+ * @f$\epsilon_0@f$ is the dielectric permittivity of vacuum,
+ * @f$\epsilon_{RF}@f$ is a reaction-field epsilon value, @f$V@f$ is the volume
+ * and @f$k_BT@f$ is the absolute temperature multiplied by the Boltzmann 
+ * constant. 
+ * 
+ * If the \@autocorr flag is given, it also writes out the normalized autocorrelation function @f$\Theta(\tau)@f$ of the total dipole moment of the system: 
+ * @f[ \Theta(\tau) = \frac{<\vec{M}(t)\cdot\vec{M}(t+\tau)>_t}{<M(t)^2>_t} = exp(-\frac{t}{\tau_\phi})@f]
+ * from which the Debye relaxation time @f$\tau_D@f$ can be calculated [J. Chem. Phys. 82, 5663 (1985)]:
+ * @f[ \tau_D = \frac{2\epsilon_{rf} + \epsilon(0)}{2\epsilon_{rf} + 1} \tau_{\phi}@f] 
+ * Using the \@truncate flag, the autocorrelation output can be truncated when the number of independent contributing frames drops below a given threshold.
+ *
+ * Note that the total dipole moment of the system is only well-defined for 
+ * systems consisting of neutral molecules. If the system carries a net-charge,
+ * the dipole moment will depend on the position of the origin. In cases where
+ * the overall system is neutral but contains ions, the dipole moment will
+ * depend on which periodic copy of the ions is taken. In these cases, cos\_epsilon
+ * issues a warning that results will critically depend on the choice of
+ * gathering method.
+ *
+ * <b>arguments:</b>
+ * <table border=0 cellpadding=0>
+ * <tr><td> \@topo</td><td>&lt;molecular topology file&gt; </td></tr>
+ * <tr><td> \@pbc</td><td>&lt;boundary type&gt; [&lt;gather method&gt;] </td></tr>
+ * <tr><td> [\@time</td><td>&lt;@ref utils::Time "time and dt"&gt;] </td></tr>
+ * <tr><td> [\@e_rf</td><td>&lt;reaction field epsilon&gt;] </td></tr>
+ * <tr><td> [\@fac</td><td>&lt;conversion factor for the unit of the dipole, default: 1; use 48.032045 to convert from e*nm to Debye&gt;] </td></tr>
+ * <tr><td> \@temp</td><td>&lt;temperature&gt; </td></tr>
+ * <tr><td> [\@autocorr</td><td>&lt;filename for storing time autocorrelation&gt;] </td></tr>
+ * <tr><td> [\@truncate</td><td>&lt;minimum number of independent contributing frames after which to truncate the correlation function&gt;] </td></tr>
+ * <tr><td> \@traj</td><td>&lt;trajectory files&gt; </td></tr>
+ * <tr><td> \@trs</td><td>&lt;special trajectory files with COS displacements&gt; </td></tr>
+ * </table>
+ *
+ *
+ * Example:
+ * @verbatim
+  cos_epsilon
+    @topo  ex.top
+    @pbc   r
+    [@time  0 0.2]
+    @e_rf  61
+    @fac 48.032045
+    @temp  300
+    @traj  ex.trc
+    @trs   ex.trs
+ @endverbatim
+ *
+ * <hr>
+ */
+
+#include <vector>
+#include <iomanip>
+#include <iostream>
+#include <fstream>
+#include <cassert>
+#include <locale>
+
+#include "../src/args/Arguments.h"
+#include "../src/args/BoundaryParser.h"
+#include "../src/args/GatherParser.h"
+#include "../src/fit/Reference.h"
+#include "../src/gio/InG96.h"
+#include "../src/gcore/System.h"
+#include "../src/gcore/Molecule.h"
+#include "../src/gcore/Box.h"
+#include "../src/gio/InTopology.h"
+#include "../src/bound/Boundary.h"
+#include "../src/fit/PositionUtils.h"
+#include "../src/gmath/Vec.h"
+#include "../src/gmath/Physics.h"
+#include "../src/gmath/Correlation.h"
+#include "../src/utils/AtomSpecifier.h"
+#include "../src/utils/groTime.h"
+
+using namespace std;
+using namespace fit;
+using namespace gcore;
+using namespace gio;
+using namespace bound;
+using namespace args;
+
+void get_atom_dipole(const AtomTopology &current_atom, Vec &atom_pos,
+ Vec &offsite_i, Vec &offsite_j, Vec &cosdispl,
+ Vec &mol_dip,
+ double &nchargepa)
+{
+  Vec pol_site_pos(0, 0, 0);
+  double pol_site_q = 0;
+  if (current_atom.isPolarisable())
+  {
+    if (current_atom.poloffsiteGamma() > 0.000001)
+    {
+      double off_set_gamma;
+      off_set_gamma = current_atom.poloffsiteGamma();
+
+      //position of the offset atom
+      pol_site_pos = atom_pos + off_set_gamma *
+                                    (offsite_i + offsite_j - 2 * atom_pos);
+    }
+    else
+    {
+      pol_site_pos = atom_pos;
+    }
+    pol_site_q = current_atom.charge() + (current_atom.cosCharge() * -1);
+    Vec cos_pos = pol_site_pos + cosdispl;
+
+    // contribution of the cos charges with cos position
+    mol_dip += cos_pos * current_atom.cosCharge();
+  }
+  else
+  {
+    // case no polarisation nor offset
+    pol_site_pos = atom_pos;
+    pol_site_q = current_atom.charge();
+  }
+  mol_dip += (pol_site_q - nchargepa) * pol_site_pos;
+}
+
+int main(int argc, char **argv)
+{
+
+  Argument_List knowns;
+  knowns << "topo"
+         << "pbc"
+         << "time"
+         << "e_rf"
+         << "fac"
+         << "temp"
+         << "autocorr"
+         << "truncate"
+         << "traj"
+         << "trs";
+
+  string usage = "# " + string(argv[0]);
+  usage += "\n\t@topo <molecular topology file>\n";
+  usage += "\t@pbc    <boundary type> [<gather method>]\n";
+  usage += "\t[@time   <time and dt>]\n";
+  usage += "\t@e_rf   <reaction field epsilon>\n";
+  usage += "\t[@fac   <conversion factor for the unit of the dipole, default: 1; use 48.032045 to convert from e*nm to Debye>]\n";
+  usage += "\t@temp   <temperature>\n";
+  usage += "\t[@autocorr    <filename for storing time autocorrelation, if no name given: Mxyz.out>]\n";
+  usage += "\t[@truncate    <minimum number of independent contributing frames after which to truncate the correlation function>]\n";
+  usage += "\t@traj   <trajectory files>\n";
+  usage += "\t@trs    <special traj with cosdisplacement>\n";
+
+  try
+  {
+    Arguments args(argc, argv, knowns, usage);
+
+    // get the @time argument
+    utils::Time time_trj(args), time_trs(args);
+    std::vector<double> time;
+
+    // read the temperature
+    double temp = args.getValue<double>("temp");
+
+    //  read topology
+    InTopology it(args["topo"]);
+    System sys(it.system());
+
+    System refSys(it.system());
+
+    unsigned int truncate_at = args.getValue<unsigned int>("truncate", false, 0);
+
+    bool do_autocorr = false;
+    if (args.count("autocorr") >= 0)
+      do_autocorr = true;
+
+    std::string autocorr_fname;
+    if (args.count("autocorr") == 0)
+    {
+      autocorr_fname = "Mcorr.out";
+    }
+    else if (args.count("autocorr") > 0)
+    {
+      autocorr_fname = args.getValue<std::string>("autocorr", false, "Mcorr.out");
+    }
+
+    // create an atomspecifier that contains all atoms (including the solvent)
+    utils::AtomSpecifier atoms(sys);
+    for (int m = 0; m < sys.numMolecules(); ++m)
+      for (int a = 0; a < sys.mol(m).numAtoms(); ++a)
+        atoms.addAtom(m, a);
+    atoms.addSpecifier("s:a");
+
+    //determine net charge
+    double ncharge = 0;
+    double nchargepa = 0;
+    for (unsigned int i = 0; i < atoms.size(); i++)
+    {
+      ncharge += atoms.charge(i);
+    }
+
+    // atoms.size() does not include the solvent, so for solvent-only systems:
+    int num_solv_mol = 0;
+    bool charged_solvent = false;
+    // read first frame to get number of solvent atoms
+    InG96 ic;
+    if (args.count("traj") > 0)
+    {
+      ic.open(args.lower_bound("traj")->second);
+
+      ic.select("ALL");
+      ic >> sys;
+      ic.close();
+    }
+    else
+    {
+      throw gromos::Exception("cos_epsilon", "no trajectory specified (@traj)");
+    }
+    
+    for (int i = 0; i < sys.numSolvents(); i++)
+    {
+      double solv_charge_mol = 0;
+      for (int a = 0; a < sys.sol(i).topology().numAtoms(); a++)
+        solv_charge_mol += sys.sol(i).topology().atom(a).charge();
+      if (solv_charge_mol)
+        charged_solvent = true;
+
+      int _num_solv_mol = sys.sol(i).numAtoms() / sys.sol(i).topology().numAtoms();
+      num_solv_mol += _num_solv_mol;
+      ncharge += num_solv_mol * solv_charge_mol;
+    }
+
+    if (atoms.size() + atoms.numSolventAtoms())
+      nchargepa = ncharge / (atoms.size() + atoms.numSolventAtoms());
+    else
+      nchargepa = 0;
+
+    // read e_rf
+    double e_rf = args.getValue<double>("e_rf", true);
+    // read conversion factor
+    double conversion_fac = args.getValue<double>("fac", false, 1);
+
+    // parse boundary conditions
+    Boundary *pbc = BoundaryParser::boundary(sys, args);
+    // parse gather method
+    Boundary::MemPtr gathmethod = args::GatherParser::parse(sys, refSys, args);
+
+    // write a title
+    if (ncharge != 0.0)
+    {
+      cout << "# WARNING the system carries a net charge ( "
+           << ncharge << ")\n"
+           << "#         this means that the dipole depends on the origin\n"
+           << "#         and your estimate of eps is probably wrong\n"
+           << "#\n";
+    }
+    else
+    {
+      // we also want to know if there are any ions
+      int ion_count = 0;
+      for (int m = 0; m < sys.numMolecules(); ++m)
+      {
+        double tc = 0.0;
+        for (int a = 0; a < sys.mol(m).numAtoms(); ++a)
+          tc += sys.mol(m).topology().atom(a).charge();
+        if (tc != 0.0)
+          ion_count++;
+      }
+      if (ion_count != 0 || charged_solvent)
+      {
+        cout << "# WARNING although the system is neutral overall,\n";
+        if (ion_count != 0)
+          cout << "#        - there are" << ion_count << " molecules that carry a charge\n";
+        if (charged_solvent)
+          cout << "#        - there is charged solvent\n";
+        cout << "#\n"
+             << "# -> the system-dipole will depend on their positions in "
+             << "the periodic box\n"
+             << "#         this is likely to give random results\n";
+      }
+    }
+
+    cout << "#" << setw(14) << " time"
+         << setw(15) << " boxdipole"
+         << setw(15) << " <mol.dipole>"
+         << setw(15) << " epsilon\n";
+
+    // prepare the calculation of the average volume
+    double vol = 0, sum_vol = 0, vcorr = 1.0;
+    Arguments::const_iterator iter = args.lower_bound("pbc");
+    if (iter != args.upper_bound("pbc"))
+      if (iter->second[0] == 't')
+        vcorr = 0.5;
+
+    // and values to calculate the dipole fluctuation
+    Vec sum_dip(0.0, 0.0, 0.0);
+    vector<Vec> dipole_vector;
+    double sum_dip2 = 0.0, fluc = 0.0;
+    int numFrames = 0;
+
+    // define input coordinate
+    InG96 is;
+
+    // we also need these factors
+    double fac, a, b, eps;
+    double f = 3.0 * gmath::physConst.get_eps0() * gmath::physConst.get_boltzmann() * temp * (2 * e_rf + 1.0);
+
+    // loop over all trajectories
+    for (Arguments::const_iterator
+             iter = args.lower_bound("traj"),
+             to = args.upper_bound("traj"),
+             siter = args.lower_bound("trs"),
+             sto = args.upper_bound("trs");
+         iter != to || siter != sto; ++iter, ++siter)
+    {
+
+      // open file
+      ic.open((iter->second).c_str());
+      ic.select("ALL");
+
+      //open special traj
+      is.open((siter->second).c_str());
+      is.select("ALL");
+
+      // loop over single trajectory
+      while (!ic.eof() && !is.eof())
+      {
+        is >> sys >> time_trs;
+
+        ic >> sys >> time_trj;
+
+        // check if the times are the same
+        if (time_trs.time() != time_trj.time())
+        {
+          cerr << "Time is not corresponding between trajectories\n"
+               << "Special traj:\t" << time_trs.time() << "\n"
+               << "Position traj:\t" << time_trj.time() << endl;
+          return -1;
+        }
+
+        // I store it because even when no @time flag is given I want to determine a timestep
+        // as the difference between the first and second time point for the correlation
+        time.push_back(time_trs.time());
+
+        (*pbc.*gathmethod)();
+
+        // calculate the volume
+        sys.box().update_triclinic();
+        vol = vcorr * sys.box().K_L_M();
+        sum_vol += vol;
+
+        Vec tot_dipole(0.0, 0.0, 0.0);
+        double sum_mol_dip_frame = 0;
+        for (int m = 0; m < sys.numMolecules(); m++)
+        {
+          Vec mol_dip(0, 0, 0);
+
+          Vec offsite_i(0, 0, 0), offsite_j(0, 0, 0);
+
+          Molecule molecule = sys.mol(m);
+
+          for (int a = 0; a < molecule.numAtoms(); a++)
+          {
+            const AtomTopology &current_atom = molecule.topology().atom(a);
+
+            Vec atom_pos = molecule.pos(a);
+            if (current_atom.isPolarisable() && current_atom.poloffsiteGamma() > 0.000001)
+            {
+              offsite_i = molecule.pos(current_atom.poloffsiteI());
+              offsite_j = molecule.pos(current_atom.poloffsiteJ());
+            }
+            Vec cosdispl = molecule.cosDisplacement(a);
+
+            get_atom_dipole(current_atom, atom_pos, offsite_i, offsite_j, cosdispl, mol_dip, nchargepa);
+          }
+          tot_dipole += mol_dip;
+          sum_mol_dip_frame += mol_dip.abs();
+        }
+
+        uint num_solv_molecules = 0;
+        for (int s = 0; s < sys.numSolvents(); s++)
+        {
+          int num_solvent_atoms = sys.sol(s).topology().numAtoms();
+          //for (unsigned int i=0; i < sys.sol(s).numAtoms(); i++)
+          //unsigned int i=0;
+
+          for (int i = 0; i < sys.sol(s).numAtoms(); i += num_solvent_atoms)
+          {
+            num_solv_molecules++;
+            Vec mol_dip(0, 0, 0);
+
+            Vec offsite_i(0, 0, 0), offsite_j(0, 0, 0);
+
+            //cerr << "#Amount of atoms in molecule:\t" << sys.mol(m).numAtoms() << endl;
+            for (int a = 0; a < num_solvent_atoms; a++)
+            {
+              const AtomTopology &current_atom = sys.sol(s).topology().atom(a);
+
+              Vec atom_pos = sys.sol(s).pos(i + a);
+
+              if (current_atom.isPolarisable() && current_atom.poloffsiteGamma() > 0.000001)
+              {
+                offsite_i = sys.sol(s).pos(i + a + current_atom.poloffsiteI());
+                offsite_j = sys.sol(s).pos(i + a + current_atom.poloffsiteJ());
+              }
+              Vec cosdispl = sys.sol(s).cosDisplacement(i + a);
+
+              get_atom_dipole(current_atom, atom_pos, offsite_i, offsite_j, cosdispl, mol_dip, nchargepa);
+            }
+
+            tot_dipole += mol_dip;
+            sum_mol_dip_frame += mol_dip.abs();
+          }
+        }
+        sum_mol_dip_frame = sum_mol_dip_frame / (sys.numMolecules() + num_solv_molecules);
+
+        // store for autocorrelation
+        dipole_vector.push_back(tot_dipole);
+
+        // do some bookkeeping
+        numFrames++;
+
+        sum_dip2 += tot_dipole.abs2();
+        sum_dip += tot_dipole;
+
+        // calculate the fluctuations of the dipole moment
+        fluc = sum_dip2 / numFrames - (sum_dip / numFrames).abs2();
+
+        // calculate the current estimate of eps
+        fac = f * sum_vol / numFrames;
+        a = 2 * e_rf * fluc + fac;
+        b = -fluc + fac;
+        eps = a / b;
+
+        cout << time_trj
+             << setw(15) << setprecision(8) << tot_dipole.abs() * conversion_fac
+             << setw(15) << setprecision(8) << sum_mol_dip_frame * conversion_fac
+             << setw(15) << setprecision(8) << eps
+             << endl;
+      }
+      ic.close();
+      is.close();
+    }
+
+    sum_dip2 /= numFrames;
+
+    double dt;
+    if (time.size() > 1)
+      dt = time[1] - time[0];
+    else
+      dt = 1;
+
+    if (do_autocorr)
+    {
+      //now doing the auto correlation calculation
+
+      ofstream os(autocorr_fname.c_str());
+      os << "#\n";
+      os << "# normalized autocorrelation function C(t) of the box dipole\n";
+      os << "#" << setw(11) << "time" << setw(20) << "C(t)" << setw(15) << "indep.frames\n";
+
+      gmath::Correlation *corr;
+      // several cases
+      corr = new gmath::Correlation(dipole_vector, dipole_vector);
+      corr->calc_direct();
+
+      double tau = 0;
+      for (unsigned int i = 0; i < corr->size(); i++, tau += dt)
+      {
+        uint independent_contrib = 0;
+        if (i != 0)
+        {
+          independent_contrib = corr->size() / i;
+          if (independent_contrib < truncate_at)
+            break;
+        }
+
+        os << setw(12) << tau
+           << setw(20) << (*corr)[i] / sum_dip2
+           << setw(15) << independent_contrib
+           << std::endl;
+      }
+
+      delete pbc;
+      delete corr;
+      os.close();
+    }
+  }
+
+  catch (const gromos::Exception &e)
+  {
+    cerr << e.what() << endl;
+    exit(1);
+  }
+  return 0;
+}
diff --git a/gromos++/programs/cry_rms.cc b/gromos++/programs/cry_rms.cc
index b265c574..f6aaa70b 100644
--- a/gromos++/programs/cry_rms.cc
+++ b/gromos++/programs/cry_rms.cc
@@ -214,7 +214,7 @@ int main(int argc, char **argv) {
         asu_pointer.addSpecifier(spec);
       }
     }
-    if (asu_pointer.size() != int(num_symop)) {
+    if (asu_pointer.size() != (unsigned int)(num_symop)) {
       ostringstream os;
       os << "Please give the first atom of every ASU for @asuspec. Got "
               << asu_pointer.size() << " atoms but need " << num_symop << ":" << endl;
@@ -241,7 +241,7 @@ int main(int argc, char **argv) {
         atomsrmsd[0].addSpecifier(spec);
       }
       for (unsigned int i = 1; i < num_symop; ++i) {
-        for (int j = 0; j < atomsrmsd[0].size(); ++j) {
+        for (unsigned int j = 0; j < atomsrmsd[0].size(); ++j) {
           int rel_pointer = atomsrmsd[0].gromosAtom(j) - asu_pointer.gromosAtom(0);
           int atom = asu_pointer.gromosAtom(i) + rel_pointer;
           if (rel_pointer < 0 || atom >= num_atoms) {
@@ -264,7 +264,7 @@ int main(int argc, char **argv) {
       }
       atomsrmsf[0].sort();
       for (unsigned int i = 1; i < num_symop; ++i) {
-        for (int j = 0; j < atomsrmsf[0].size(); ++j) {
+        for (unsigned int j = 0; j < atomsrmsf[0].size(); ++j) {
           int rel_pointer = atomsrmsf[0].gromosAtom(j) - asu_pointer.gromosAtom(0);
           int atom = asu_pointer.gromosAtom(i) + rel_pointer;
           if (rel_pointer < 0 || atom >= num_atoms) {
@@ -376,7 +376,7 @@ int main(int argc, char **argv) {
             for (unsigned int j = i + 1; j < num_symop; ++j) {
               // calculate RMSD between each pair.
               double sum = 0.0;
-              for (int a = 0; a < atomsrmsd[i].size(); ++a) {
+              for (unsigned int a = 0; a < atomsrmsd[i].size(); ++a) {
                 Vec pos_i(inv_rotation[i] * (atomsrmsd[i].pos(a) - translation[i]));
                 Vec pos_j(inv_rotation[j] * (atomsrmsd[j].pos(a) - translation[j]));
                 Vec d(pos_i - pbc->nearestImage(pos_i, pos_j, sys.box()));
diff --git a/gromos++/programs/dGslv_pbsolv.cc b/gromos++/programs/dGslv_pbsolv.cc
index dd9435af..63ecbf57 100644
--- a/gromos++/programs/dGslv_pbsolv.cc
+++ b/gromos++/programs/dGslv_pbsolv.cc
@@ -207,7 +207,7 @@ int main(int argc, char **argv){
     if (atoms_to_charge.size()==0)  throw gromos::Exception("dGslv_pbsolv","No atoms to charge specified. Exiting ...");
 
    std::cout << "# READ: atoms_to_charge " << endl;
-   for (int i=0;i< atoms_to_charge.size();i++){
+   for (unsigned int i=0;i< atoms_to_charge.size();i++){
        std::cout << "# READ: mol " <<  atoms_to_charge.mol(i) << " " << atoms_to_charge.name(i) << endl;
    }
 
@@ -460,7 +460,7 @@ int main(int argc, char **argv){
 
       // loop over atoms and set H radii to hradius_param
       // also, if sigma is used, adapt appropriately
-      for (int i=0;i<atoms.size();i++){
+      for (unsigned int i=0;i<atoms.size();i++){
   
           if (fabs(atoms.radius(i))< ppp.tiny_real){
      //     sys.mol(atoms.mol(i)).topology().atom(atoms.atom(i)).setradius( ppp.get_hradius() );
@@ -484,7 +484,7 @@ int main(int argc, char **argv){
           }
           std::cout << "# radius of atom " << i << " : rad(i) = " << atoms.radius(i)  << endl;
       }
-       for (int i=0;i<atoms_to_charge.size();i++){
+       for (unsigned int i=0;i<atoms_to_charge.size();i++){
           std::cout << "# radius of atom_to_charge " << i << " : rad(i) = " << atoms_to_charge.radius(i)  << endl;
       }
 
diff --git a/gromos++/programs/dfgrid.cc b/gromos++/programs/dfgrid.cc
index 9f2b2734..252ff2fe 100644
--- a/gromos++/programs/dfgrid.cc
+++ b/gromos++/programs/dfgrid.cc
@@ -242,15 +242,6 @@ int main(int argc, char **argv) {
         }
     }
 
-      
-
-    // get simulation time either from the user or from the files
-    bool usertime=false;
-    
-    if (args.count("time") > 0) {
-      usertime=true;
-    }
-
     // parse boundary conditions
     Boundary *pbc = BoundaryParser::boundary(sys, args);
     //parse gather method
@@ -345,7 +336,7 @@ int main(int argc, char **argv) {
     
     if (args.count("distatoms") >0)  {
       std::cout << "# ";
-      for (int i=0; i<distatoms.size(); i++) {
+      for (unsigned int i=0; i<distatoms.size(); i++) {
         std::cout << distatoms.toString(i)<<" ";
       }
       std::cout << "# in protein" << std::endl;
@@ -425,7 +416,7 @@ int main(int argc, char **argv) {
               distance[j] = 4*(boxK+boxL+boxM);
 
           std::set<int> protein;
-          for (int a = 0; a < proteinatoms.size(); ++a){
+          for (unsigned int a = 0; a < proteinatoms.size(); ++a){
             gmath::Vec ppos = pbc->nearestImage(origin, proteinatoms.pos(a), box);
             int nx_min=int(-( proteincutoff - ppos[0] - boxK/2)/gridspacing);
             int nx_max=int(-(-proteincutoff - ppos[0] - boxK/2)/gridspacing)+1;
@@ -552,7 +543,7 @@ int main(int argc, char **argv) {
             if (args.count("distatoms")>0) {
               std::cout << numFrames << " ";
               ostringstream ostr;
-              for (int a=0; a < distatoms.size(); a++) {
+              for (unsigned int a=0; a < distatoms.size(); a++) {
                 gmath::Vec pos = pbc->nearestImage(origin,distatoms.pos(a), box);
                 int x = int(-(-pos[0] - boxK/2)/gridspacing+0.5);
                 int y = int(-(-pos[1] - boxL/2)/gridspacing+0.5);
@@ -583,7 +574,6 @@ int main(int argc, char **argv) {
 
             // add atoms to the new molecule's topology
             // VA atom for center
-            int numva = distatoms.size()+1;
             int resnum = 1;
             int atomnum = 0;
             AtomTopology at_va;
@@ -597,7 +587,7 @@ int main(int argc, char **argv) {
             // DA atoms for distatoms, each as separate residue
             AtomTopology at_da;
             at_da.setName("DA");
-            for (int i=0; i<distatoms.size();i++) {
+            for (unsigned int i=0; i<distatoms.size();i++) {
               mt.addAtom(at_da);
               mt.setResNum(atomnum,resnum);
               mt.setResName(resnum,"DA");
@@ -608,8 +598,8 @@ int main(int argc, char **argv) {
             // PA: path atoms
             AtomTopology at_pa;
             at_pa.setName("PA");
-            for (int i=0; i<distatoms.size();i++) {
-              for (int j=0; j<minpaths[i].size();j++) {
+            for (unsigned int i=0; i<distatoms.size();i++) {
+              for (unsigned int j=0; j<minpaths[i].size();j++) {
                 mt.addAtom(at_pa);
                 if (j!=0) {
                   Bond bond(atomnum-1,atomnum);
@@ -645,14 +635,14 @@ int main(int argc, char **argv) {
             atomnum++;
 
             // add distatom coordinates        
-            for (int i=0; i<distatoms.size();i++) {
+            for (unsigned int i=0; i<distatoms.size();i++) {
               gridsys.mol(nummol).pos(atomnum) = distatoms.pos(i);
               atomnum++;
             }        
 
             // add path coordinates
-            for (int i=0; i<distatoms.size();i++) {
-              for (int j=0; j<minpaths[i].size();j++) {
+            for (unsigned int i=0; i<distatoms.size();i++) {
+              for (unsigned int j=0; j<minpaths[i].size();j++) {
                 int idx = minpaths[i][j];
                 gmath::Vec gpos = get_coords(idx, ngrid,gridspacing,boxK,boxL,boxM);
                 gridsys.mol(nummol).pos(atomnum) = gpos;
diff --git a/gromos++/programs/dg_ener.cc b/gromos++/programs/dg_ener.cc
index 71670878..fe258ea0 100644
--- a/gromos++/programs/dg_ener.cc
+++ b/gromos++/programs/dg_ener.cc
@@ -118,7 +118,7 @@ int main(int argc, char **argv){
     vector<double> delta_v;
     vector<double> t;
 
-    double covA, nbA, totA, covB, nbB, totB;
+    double totA, totB; //covA, nbA, covB, nbB
 
     while(true){
 
diff --git a/gromos++/programs/diffus.cc b/gromos++/programs/diffus.cc
index a71c30a1..564ee2f5 100644
--- a/gromos++/programs/diffus.cc
+++ b/gromos++/programs/diffus.cc
@@ -272,7 +272,7 @@ int main(int argc, char **argv) {
     
     // calculate the com of the reference state
     Vec com0(0.0, 0.0, 0.0);
-    for (int i = 0; i < ref_at.size(); i++)
+    for (unsigned int i = 0; i < ref_at.size(); i++)
       com0 += ref_at.pos(i);
     com0 /= ref_at.size();
 
@@ -312,7 +312,7 @@ int main(int argc, char **argv) {
         Stat<double> disp_data;
         //loop over the molecules
 
-        for(int i = 0; i < at.size(); i++) {
+        for(unsigned int i = 0; i < at.size(); i++) {
           at.pos(i) = pbc->nearestImage(old_at.pos(i), at.pos(i), sys.box());
           position[ i ].push_back(at.pos(i));
           old_at.pos(i) = at.pos(i);
@@ -338,7 +338,7 @@ int main(int argc, char **argv) {
       //for(int j = it; j < it+window; j++) {
         counter++;
         //cout << "  j : " << j << "  it : " << it << endl;
-        for(int i = 0; i < at.size(); i++) {
+        for(unsigned int i = 0; i < at.size(); i++) {
           square_d = 0;
 
           for(int k = 0; k < ndim; k++) {
diff --git a/gromos++/programs/dipole.cc b/gromos++/programs/dipole.cc
index fc97abb1..d4f651e6 100644
--- a/gromos++/programs/dipole.cc
+++ b/gromos++/programs/dipole.cc
@@ -110,7 +110,7 @@ int main(int argc, char **argv){
     //determine net charge
     double ncharge=0;
     double nchargepa=0;
-    for(int i=0; i < atoms.size(); i++){
+    for(unsigned int i=0; i < atoms.size(); i++){
       ncharge+=atoms.charge(i);
     }
     nchargepa=ncharge/atoms.size();
@@ -162,12 +162,12 @@ int main(int argc, char **argv){
 	(*pbc.*gathmethod)();
 	Vec cog(0.0,0.0,0.0);
 	if(lcog){
-	  for(int i=0; i<atoms.size(); ++i)
+	  for(unsigned int i=0; i<atoms.size(); ++i)
 	    cog+=atoms.pos(i);
 	  cog/=atoms.size();
 	}
 	Vec dipole(0.0,0.0,0.0);
-	for(int i=0; i<atoms.size(); ++i){
+	for(unsigned int i=0; i<atoms.size(); ++i){
 	  dipole += (atoms.charge(i)-nchargepa)*(atoms.pos(i)-cog);
 	}
 	
diff --git a/gromos++/programs/disicl.cc b/gromos++/programs/disicl.cc
index 5b16df17..84c5bc8d 100644
--- a/gromos++/programs/disicl.cc
+++ b/gromos++/programs/disicl.cc
@@ -210,7 +210,6 @@ int main(int argc, char **argv) {
     }
     
     // define input coordinate
-    bool do_bfac=true;
     InG96 ic(skip,stride);
     
     dscl.initSummary();
diff --git a/gromos++/programs/edyn.cc b/gromos++/programs/edyn.cc
index 084bde1b..3acc9b07 100644
--- a/gromos++/programs/edyn.cc
+++ b/gromos++/programs/edyn.cc
@@ -558,7 +558,7 @@ void writePdb(const char *name, AtomSpecifier &atoms){
   out.setf(ios::fixed, ios::floatfield);
   out.setf(ios::unitbuf);
   out.precision(3);
-  for (int h=0; h<atoms.size(); ++h){
+  for (unsigned int h=0; h<atoms.size(); ++h){
     int res=atoms.resnum(h);
     out << "ATOM";
     out.setf(ios::right, ios::adjustfield);
diff --git a/gromos++/programs/ene_ana.cc b/gromos++/programs/ene_ana.cc
index b9531a77..660465d3 100644
--- a/gromos++/programs/ene_ana.cc
+++ b/gromos++/programs/ene_ana.cc
@@ -195,20 +195,20 @@ int main(int argc, char **argv){
       to_en=args.upper_bound("en_files"),
       it_fr=args.lower_bound("fr_files"),
       to_fr=args.upper_bound("fr_files");
-    int cont=0, en_cont=0, fr_cont=0;
+    //int en_cont=0, fr_cont=0, cont=0;
     if(do_energy_files) {
       gin_en.open(it_en->second.c_str()); 
       ++it_en; 
-      en_cont=1;
+      //en_cont=1;
     }
     
     if(do_free_energy_files) {
       gin_fr.open(it_fr->second.c_str());
       ++it_fr;
-      fr_cont=1;
+      //fr_cont=1;
     }
     
-    cont=en_cont+fr_cont;
+    //cont=en_cont+fr_cont;
     
     bool version_checked = false;
     while (true) {
diff --git a/gromos++/programs/epath.cc b/gromos++/programs/epath.cc
index 49c17c00..34d93262 100644
--- a/gromos++/programs/epath.cc
+++ b/gromos++/programs/epath.cc
@@ -443,7 +443,7 @@ int main(int argc, char **argv) {
             withH.addSpecifier("a:a");
 
             // remove H
-            for (int i = 0; i < withH.size(); i++) {
+            for (unsigned int i = 0; i < withH.size(); i++) {
                 if (!sys.mol(withH.mol(i)).topology().atom(withH.atom(i)).isH()) {
                     allatoms.addAtom(withH.mol(i), withH.atom(i));
                 }
@@ -547,7 +547,7 @@ int main(int argc, char **argv) {
                 }
 
                 bool domore = true;
-                while ((int)visited.size() != allatoms.size() && domore) {
+                while (visited.size() != allatoms.size() && domore) {
                     ///////////////////
 
                     // create specifier pairlist
@@ -567,7 +567,7 @@ int main(int argc, char **argv) {
 
                     Neighbours nb(sys, current.mol(0), current.atom(0));
                     // put the pairlist into the pair specifier
-                    for (int i = 0; i < around.size(); ++i) {
+                    for (unsigned int i = 0; i < around.size(); ++i) {
                         // check if atom is visited
                         if (!sys.mol(around.mol(i)).topology().atom(around.atom(i)).isH()
                                 && !visited.count(allatoms.findAtom(around.mol(i), around.atom(i)))) {
@@ -592,7 +592,7 @@ int main(int argc, char **argv) {
                     for (unsigned int i = 0; i < nb.size(); i++) {
                         if (sys.mol(current.mol(0)).topology().atom(nb[i]).isH()) {
                             // now search in the pairs list for a suitable acceptor
-                            for (int j = 0; j < pairs.size(); j++) {
+                            for (unsigned int j = 0; j < pairs.size(); j++) {
                                 // loop over masses
                                 bool couldbeh = false;
                                 for (unsigned int h = 0; h < masses.size(); h++) {
@@ -625,7 +625,7 @@ int main(int argc, char **argv) {
                         }
                     }
                     if (couldbeh) {
-                        for (int i = 0; i < pairs.size(); i++) {
+                        for (unsigned int i = 0; i < pairs.size(); i++) {
                             Neighbours pnb(sys, pairs.mol(i), pairs.atom(i));
                             for (unsigned int j = 0; j < pnb.size(); j++) {
                                 if (sys.mol(pairs.mol(i)).topology().atom(pnb[j]).isH()) {
@@ -646,7 +646,7 @@ int main(int argc, char **argv) {
                     }
 
                     // remove the hb from the pairs
-                    for (int i = 0; i < hb.size(); i++) {
+                    for (unsigned int i = 0; i < hb.size(); i++) {
                         pairs.removeAtom(hb.mol(i), hb.atom(i));
                     }
 
@@ -671,7 +671,7 @@ int main(int argc, char **argv) {
                     int jtype;
 
                     // covalent
-                    for (int i = 0; i < cov.size(); i++) {
+                    for (unsigned int i = 0; i < cov.size(); i++) {
                         // current in the donor or acceptor range?
                         if (((start.findAtom(cov.mol(i), cov.atom(i)) != -1) && (start.findAtom(current.mol(0), current.atom(0)) != -1))|| ((stop.findAtom(cov.mol(i), cov.atom(i)) != -1) && (stop.findAtom(current.mol(0), current.atom(0)) != -1))) {
                             ec = currentdecay * 1.00;
@@ -696,7 +696,7 @@ int main(int argc, char **argv) {
                     }
 
                     // hbond
-                    for (int i = 0; i < hb.size(); i++) {
+                    for (unsigned int i = 0; i < hb.size(); i++) {
                         // current in the donor or acceptor range?
                         if (((start.findAtom(hb.mol(i), hb.atom(i)) != -1) && (start.findAtom(current.mol(0), current.atom(0)) != -1))|| ((stop.findAtom(hb.mol(i), hb.atom(i)) != -1) && (stop.findAtom(current.mol(0), current.atom(0)) != -1))) {
                             eh = currentdecay * 1.00;
@@ -721,7 +721,7 @@ int main(int argc, char **argv) {
                     }
 
                     // space
-                    for (int i = 0; i < pairs.size(); i++) {
+                    for (unsigned int i = 0; i < pairs.size(); i++) {
                         // current in the donor or acceptor range?
                         if (((start.findAtom(pairs.mol(i), pairs.atom(i)) != -1) && (start.findAtom(current.mol(0), current.atom(0)) != -1))|| ((stop.findAtom(pairs.mol(i), pairs.atom(i)) != -1) && (stop.findAtom(current.mol(0), current.atom(0)) != -1))) {
                             es = currentdecay * 1.00;
@@ -951,7 +951,7 @@ int main(int argc, char **argv) {
                     int atomoffset = 0;
                     int m = 0;
                     int atomsoff = 0;
-                    for (int i = 0; i < allatoms.size(); i++) {
+                    for (unsigned int i = 0; i < allatoms.size(); i++) {
 
                         if (allatoms.mol(i) > m) {
                             resoffset += sys.mol(allatoms.mol(i) - 1).topology().numRes();
@@ -987,7 +987,7 @@ int main(int argc, char **argv) {
 
 
                     // Chain B
-                    for (int i = 0; i < path.size(); i++) {
+                    for (unsigned int i = 0; i < path.size(); i++) {
                         atomsoff = 0;
                         // complicated way of doing the resoffset
                         resoffset = 0;
@@ -1052,7 +1052,7 @@ int main(int argc, char **argv) {
                     int atomoffset_connect1plus = 0;
                     int atomoffset_plus1 = 0;
 
-                    for (int i = 0; i < path.size() - 1; i++) {
+                    for (unsigned int i = 0; i < path.size() - 1; i++) {
                         // complicated way of doing the offset
                         // reset
                         atomoffset_connect1 = 0;
diff --git a/gromos++/programs/eps_field.cc b/gromos++/programs/eps_field.cc
index 688c09d0..fbe6e785 100644
--- a/gromos++/programs/eps_field.cc
+++ b/gromos++/programs/eps_field.cc
@@ -118,12 +118,12 @@ int main(int argc, char **argv) {
     cout << "#Note: Only the dipole moment in z-direction was considered for calculation." << endl;
 
     // parse boundary conditions
-    Boundary *pbc = BoundaryParser::boundary(sys, args);
+    //Boundary *pbc = BoundaryParser::boundary(sys, args);
 
     // prepare the calculation of the average volume
-    double vcorr = 1.0;
-    if (pbc->type() == 't')
-      vcorr = 0.5;
+    //double vcorr = 1.0;
+    //if (pbc->type() == 't')
+    //  vcorr = 0.5;
 
     // and values to calculate the dipole fluctuation
     double sum_pz = 0.0, sum_px = 0.0, sum_py = 0.0;
diff --git a/gromos++/programs/epsilon.cc b/gromos++/programs/epsilon.cc
index 0e567e19..7836aec3 100644
--- a/gromos++/programs/epsilon.cc
+++ b/gromos++/programs/epsilon.cc
@@ -166,7 +166,7 @@ void calc_update(System sys,utils::Time time,utils::AtomSpecifier neut_mol){
 	Vec Mj(0.0,0.0,0.0);
 
 	//loop over neutral molecules
-	for (int a=0;a<neut_mol.size();a++)
+	for (unsigned int a=0;a<neut_mol.size();a++)
 		Md+=(neut_mol.pos(a))*neut_mol.charge(a);
 
 	// loop over ions
@@ -174,12 +174,12 @@ void calc_update(System sys,utils::Time time,utils::AtomSpecifier neut_mol){
 	Vec temp_com;
 	// solute
 	int m;
-        for(int ch_m=0; ch_m<mol_ions.size(); ++ch_m) {
+        for(unsigned int ch_m=0; ch_m<mol_ions.size(); ++ch_m) {
 		m=mol_ions[ch_m];
 		temp_atoms.clear();
 		temp_atoms.addMolecule(m);
 		temp_com=PositionUtils::com(sys, temp_atoms);
-		for(int a=0; a < temp_atoms.size(); a++) {
+		for(unsigned int a=0; a < temp_atoms.size(); a++) {
 			Md+=(temp_atoms.pos(a)-temp_com)*temp_atoms.charge(a);
 		}
 		Mj+=temp_com*mol_net_ch[ch_m];
@@ -189,7 +189,7 @@ void calc_update(System sys,utils::Time time,utils::AtomSpecifier neut_mol){
 	int temp_sol_atom=0;
 	string temp_spec;
 	ostringstream ostrs;
-	int ch_s;
+	unsigned int ch_s;
         for(int s=0; s<sys.numSolvents(); ++s) {
 		ch_s=find(sol_ions.begin(), sol_ions.end(), s)-sol_ions.begin();
 		if (ch_s != sol_ions.size()){
@@ -204,7 +204,7 @@ void calc_update(System sys,utils::Time time,utils::AtomSpecifier neut_mol){
 				temp_spec+=ostrs.str();
 				temp_atoms.addSpecifier(temp_spec);
 				temp_com=PositionUtils::com(sys, temp_atoms);
-				for(int a=0; a < temp_atoms.size(); a++) {
+				for(unsigned int a=0; a < temp_atoms.size(); a++) {
 					Md+=(temp_atoms.pos(a)-temp_com)*temp_atoms.charge(a);
 				}
 				Mj+=temp_com*sol_net_ch[ch_s];
@@ -220,7 +220,7 @@ void calc_update(System sys,utils::Time time,utils::AtomSpecifier neut_mol){
         sum_Md2+= Md.abs2();
 
 	// calculate the average sum_Md2
-	double Md2_avg=sum_Md2 / numFrames;
+	//double Md2_avg=sum_Md2 / numFrames;
 
 	// calculate the current estimate of eps
 	sum_Md+=Md;
@@ -247,7 +247,7 @@ void calc_update_Mj_p(System sys,utils::Time time,utils::AtomSpecifier neut_mol)
 	Vec temp_com_v;
 	// solute
 	int m;
-        for(int ch_m=0; ch_m<mol_ions.size(); ++ch_m) {
+        for(unsigned int ch_m=0; ch_m<mol_ions.size(); ++ch_m) {
 		m=mol_ions[ch_m];
 		temp_atoms.clear();
 		temp_atoms.addMolecule(m);
@@ -262,7 +262,7 @@ void calc_update_Mj_p(System sys,utils::Time time,utils::AtomSpecifier neut_mol)
 	int ch_s;
         for(int s=0; s<sys.numSolvents(); ++s) {
 		ch_s=find(sol_ions.begin(), sol_ions.end(), s)-sol_ions.begin();
-		if (ch_s != sol_ions.size()){
+		if (ch_s != (int)sol_ions.size()){
 			for (int m=0; m<(sys.sol(s).numAtoms()/sys.sol(s).topology().numAtoms()); ++m) {
 				temp_atoms.clear();
 				ostrs.str("");
@@ -433,7 +433,7 @@ int main(int argc, char **argv){
     // loop over all trajectories
     iter=args.lower_bound("traj");
     Arguments::const_iterator to=args.upper_bound("traj");
-    for(iter; iter!=to; ++iter){
+    for(; iter!=to; ++iter){
       
       // open file
       ic.open((iter->second).c_str());
@@ -470,7 +470,7 @@ int main(int argc, char **argv){
 	InG96 ic_v;
 	Arguments::const_iterator iter=args.lower_bound("traj_vel");
 	Arguments::const_iterator to=args.upper_bound("traj_vel");
-	for(iter; iter!=to; ++iter){
+	for(; iter!=to; ++iter){
 	// open file
 		ic_v.open((iter->second).c_str());
 		ic_v.select("ALL");
diff --git a/gromos++/programs/ext_ti_ana.cc b/gromos++/programs/ext_ti_ana.cc
index b8176106..c7e2a79e 100644
--- a/gromos++/programs/ext_ti_ana.cc
+++ b/gromos++/programs/ext_ti_ana.cc
@@ -296,7 +296,6 @@ int main(int argc, char **argv){
   try{
     #ifdef OMP
     double start_total = omp_get_wtime();
-    double start;
     #endif
     Arguments args(argc, argv, knowns, usage);
     // check whether we have both energy and free energy files
@@ -335,7 +334,7 @@ int main(int argc, char **argv){
     //we will do some numerical comparisons between lambda values that may lead to problems
     double lam_epsilon = pow(10, -(lam_precision+2));
     
-    int num_cpus = args.getValue<int>("cpus",false, 1);
+    unsigned int num_cpus = args.getValue<unsigned int>("cpus",false, 1);
 
     // which contributions do we include?
     bool no_lj = true;
@@ -382,7 +381,7 @@ int main(int argc, char **argv){
     std::string tmp[] = {"kin", "bond", "ang", "dih", "impr", "lj", "slj", "crf", "scrf"}; 
     std::vector<std::string> lambdas_types (tmp, tmp+(sizeof(tmp)/sizeof(std::string)));
     
-    for (int i=0; i<lambdas_types.size(); i++) {
+    for (unsigned int i=0; i<lambdas_types.size(); i++) {
        // ref: by default 0 0 0 1 0
        vector<double> coeffs_sim = args.getValues<double>("lam"+lambdas_types[i]+"_sim",5,false, Arguments::Default<double>() << 0 << 0 << 0 << 1 << 0);
        abcde_sim[lambdas_types[i]] = coeffs_sim;
@@ -497,7 +496,7 @@ int main(int argc, char **argv){
     // prepare Lambdas
     std::map<std::string, LambdaStruct> reflambdas;
     std::map<std::string, std::vector<LambdaStruct> > lambdas;
-    for (int i=0; i<lambdas_types.size(); i++) {
+    for (unsigned int i=0; i<lambdas_types.size(); i++) {
        std::string lambda_type = lambdas_types[i];
        reflambdas[lambda_type] = getSingleLambda(abcde_sim[lambda_type], minlam, lam_step, nrlambdas, slam, p_LAM, lambda_type);
        // if no flag given, by default the lambda coefficients for 
@@ -598,20 +597,20 @@ int main(int argc, char **argv){
       to_fr=args.upper_bound("fr_files");
 
 
-    int cont=0, en_cont=0, fr_cont=0;
+    //int cont=0, en_cont=0, fr_cont=0;
     if(do_energy_files) {
       gin_en.open(it_en->second.c_str()); 
       ++it_en; 
-      en_cont=1;
+      //en_cont=1;
     }
     
     if(do_free_energy_files) {
       gin_fr.open(it_fr->second.c_str());
       ++it_fr;
-      fr_cont=1;
+      //fr_cont=1;
     }
 
-    cont=en_cont+fr_cont;
+    //cont=en_cont+fr_cont;
     bool version_checked = false;
     
     // starting with the calculations    
@@ -747,8 +746,8 @@ int main(int argc, char **argv){
         //calculate lj and crf contributions for all slj and scrf parameterizations
         double Emult;
         double dEmult;
-        for(int i=0; i<SLJ.size(); i++){
-          for(int j=0; j<SCRF.size(); j++){
+        for(unsigned int i=0; i<SLJ.size(); i++){
+          for(unsigned int j=0; j<SCRF.size(); j++){
           // calculate energies (with interpolation between lambda_i)
           Emult = E;
           Emult += pow1_lj * (SLJ[i][p].wlow*se[SLJ[i][p].idx_low][0] + SLJ[i][p].whigh*se[SLJ[i][p].idx_high][0])
@@ -790,7 +789,7 @@ int main(int argc, char **argv){
         
     #ifdef OMP
         if (num_cpus > SLJ.size()){
-            if(SLJ.size() > omp_get_max_threads())
+            if(int(SLJ.size()) > omp_get_max_threads())
                 num_cpus = omp_get_max_threads();
             else
                 num_cpus = SLJ.size();
@@ -798,7 +797,7 @@ int main(int argc, char **argv){
                  << num_cpus << " threads." << endl;
         }
 
-        if (num_cpus > omp_get_max_threads()){
+        if ((int)num_cpus > omp_get_max_threads()){
             cerr << "# You specified " << num_cpus << " number of threads. There are only " 
                  << omp_get_max_threads() << " threads available." << endl;
             num_cpus = omp_get_max_threads();
@@ -813,13 +812,12 @@ int main(int argc, char **argv){
     //initialize random number generator for bootstraps
     int ti = time(NULL);
     #ifdef OMP
-    double prep_time = omp_get_wtime() - start_total;
     double totaltime=0;
     double bootstrap_time=0;
     #pragma omp parallel for reduction(+:totaltime,bootstrap_time) //firstprivate(sys, time) 
     #endif
     // loop over all trajectories
-    for(int i=0; i<SLJ.size(); i++){ 
+    for(unsigned int i=0; i<SLJ.size(); i++){ 
         //#ifdef OMP
         //#pragma omp critical
         //#endif
@@ -831,7 +829,7 @@ int main(int argc, char **argv){
         unsigned int seed=ti+i;
       
         abcde["slj"] = SLJ_abcde[i];
-        for(int j=0; j<SCRF.size(); j++){
+        for(unsigned int j=0; j<SCRF.size(); j++){
             abcde["scrf"] = SCRF_abcde[j];
 
             ostringstream oss;
@@ -843,7 +841,7 @@ int main(int argc, char **argv){
             
                 oss << "#----- lambda coefficients: ------------\n";
                 oss << "#        coeffs_sim    coeffs_prediction" << endl;
-                for (int ii=0; ii<lambdas_types.size(); ii++) {
+                for (unsigned int ii=0; ii<lambdas_types.size(); ii++) {
                     oss << "# " << setw(4) <<lambdas_types[ii]<< ": "
                         << abcde_sim[lambdas_types[ii]] << " | " 
                         << abcde[lambdas_types[ii]] << std::endl;
@@ -931,7 +929,6 @@ int main(int argc, char **argv){
                             #endif
                             r=rand_r(&seed);
                             int index = int((r*length)/RAND_MAX);
-                            int sign = 0;
                             x_new.addval(x[i][j][p].val(index));
                             vyvr_new.addval(vyvr[i][j][p].val(index));
                         }
@@ -1040,7 +1037,7 @@ vector <LambdaStruct> getLambdas(vector<double> abcde, double minlam, double lam
 {
   vector<LambdaStruct> newlambdas(p_LAM.size());  
 
-  for(int p=0; p<p_LAM.size(); p++){
+  for(unsigned int p=0; p<p_LAM.size(); p++){
     newlambdas[p]=getSingleLambda(abcde, minlam, lam_step, nrlambdas, p_LAM[p], p_LAM, lambda_type);
   }
 
@@ -1230,17 +1227,17 @@ std::map<std::string, vector<double> > read_lambdas(std::vector<ilambdas::lambin
   std::map<std::string, vector<double> > lambdas_map;
   bool warn=false;
   int counts[lambdas_types.size()];
-  for (int t = 0; t < lambdas_types.size(); t++) {
+  for (unsigned int t = 0; t < lambdas_types.size(); t++) {
     counts[t]=0;
   }
   
-  for (int i=0; i < lambints.size(); i++) {
+  for (unsigned int i=0; i < lambints.size(); i++) {
     double tmp[] = {lambints[i].ali, lambints[i].bli, lambints[i].cli, lambints[i].dli, lambints[i].eli};
     std::vector<double> coeffs (tmp, tmp + sizeof(tmp)/sizeof(double));
 
     // t=0 is the kinetic part, which is not defined in gromosXX LAMBDAS block
-    for (int t = 1; t < lambdas_types.size(); t++) {
-      if (lambints[i].ntli == t) {
+    for (unsigned int t = 1; t < lambdas_types.size(); t++) {
+      if (lambints[i].ntli == int(t)) {
         counts[t]++;
 
         // check if the coefficients are the same in all occurrences of this lambdas type
diff --git a/gromos++/programs/ext_ti_merge.cc b/gromos++/programs/ext_ti_merge.cc
index f1975bfe..5bcb7241 100644
--- a/gromos++/programs/ext_ti_merge.cc
+++ b/gromos++/programs/ext_ti_merge.cc
@@ -220,7 +220,7 @@ int main(int argc, char **argv) {
     if (!no_errors) interr=integrate(prediction[0],prediction[2]);
     
     
-    for (int i=0; i < prediction[0].size(); i++) {
+    for (unsigned int i=0; i < prediction[0].size(); i++) {
       cout << setw(4)<< setprecision(2) << fixed<< prediction[0][i]<< " " << setprecision(6)<< fixed << prediction[1][i];
       if (!no_errors) cout << " " <<setprecision(6)<< fixed <<  prediction[2][i];
       cout << endl;
@@ -243,7 +243,7 @@ double integrate(vector<double> x, vector<double> y) {
   
   if (x.size() != y.size()) throw gromos::Exception("integrate", "Vectors are not the same length!");
   
-  for (int i=0; i<x.size()-1; i++) {
+  for (unsigned int i=0; i<x.size()-1; i++) {
     integral+=0.5 * (y[i+1]+y[i])*(x[i+1]-x[i]);
   }
   return integral;
@@ -255,7 +255,7 @@ vector<vector<double> > weight_function(ostream &os, vector<double> &slam, map<d
     sort(slam.begin(), slam.end());
     vector<vector<double> > pred_out(3);
     
-    for (int i=0; i < slam.size()-1; i++) {
+    for (unsigned int i=0; i < slam.size()-1; i++) {
       double sl1=slam[i], sl2=slam[i+1];
       if (sl1==sl2)
         throw gromos::Exception("ext_ti_merge", "Two slam values are the same!\n");      
@@ -263,7 +263,7 @@ vector<vector<double> > weight_function(ostream &os, vector<double> &slam, map<d
       // get the relevant predicted lambdas from i and i+1 in the range between i and i+1
       vector<vector<double> > plam_tmp(2), dxdl_tmp(2), dxdl_err_tmp(2);
       for (int s=0; s<=1; s++) {
-        for (int j=0; j<vplam[slam[i+s]].size(); j++) {
+        for (unsigned int j=0; j<vplam[slam[i+s]].size(); j++) {
           if (vplam[slam[i+s]][j]>=sl1 && vplam[slam[i+s]][j] <sl2) {
             plam_tmp[s].push_back(vplam[slam[i+s]][j]);
             dxdl_tmp[s].push_back(vdxdl[slam[i+s]][j]);
@@ -278,7 +278,7 @@ vector<vector<double> > weight_function(ostream &os, vector<double> &slam, map<d
         cerr << "ext_ti_merge: no predicted lambdas between " << sl1 << " and " << sl2 << ".\n";
         continue;
       } else {
-        for (int j=0; j<plam_tmp[0].size(); j++) {
+        for (unsigned int j=0; j<plam_tmp[0].size(); j++) {
           double ws1 = (plam_tmp[0][j]-sl2)/(sl1-sl2);
           double ws2 = (plam_tmp[0][j]-sl1)/(sl2-sl1);
           double predict = dxdl_tmp[0][j]*ws1 + dxdl_tmp[1][j]*ws2;
diff --git a/gromos++/programs/filter.cc b/gromos++/programs/filter.cc
index 449da981..094624bc 100644
--- a/gromos++/programs/filter.cc
+++ b/gromos++/programs/filter.cc
@@ -132,13 +132,6 @@ int main(int argc, char **argv) {
     if (args.count("notimeblock") >= 0)
       notimeblock = true;
 
-    // get simulation time either from the user or from the files
-    bool usertime=false;
-    
-    if (args.count("time") > 0) {
-      usertime=true;
-    }
-
     // parse boundary conditions
     Boundary *pbc = BoundaryParser::boundary(sys, args);
     //parse gather method
@@ -170,7 +163,7 @@ int main(int argc, char **argv) {
     }
     // check that there are no doubles in ls and rej
     bool warn = false;
-    for (int i = 0; i < rej.size(); i++) {
+    for (unsigned int i = 0; i < rej.size(); i++) {
       if (ls.findAtom(rej.mol(i), rej.atom(i)) != -1) {
         warn = true;
       }
@@ -273,7 +266,7 @@ int main(int argc, char **argv) {
 #ifdef OMP
 #pragma omp parallel for 
 #endif
-          for (int i = 0; i < ref.size(); i++) {
+          for (unsigned int i = 0; i < ref.size(); i++) {
             utils::SimplePairlist spl(sys, *pbc, cut);
             spl.setAtom(*ref.atom()[i]);
             spl.setType(t);
@@ -291,7 +284,7 @@ int main(int argc, char **argv) {
           }
 
           // remove atoms that are to be rejected
-          for (int i = 0; i < rej.size(); i++) {
+          for (unsigned int i = 0; i < rej.size(); i++) {
             rls.removeAtom(rej.mol(i), rej.atom(i));
           }
           rls.sort();
diff --git a/gromos++/programs/follow.cc b/gromos++/programs/follow.cc
index 3d2a86ed..b5ab44a0 100644
--- a/gromos++/programs/follow.cc
+++ b/gromos++/programs/follow.cc
@@ -146,13 +146,13 @@ try{
   // we need to store the old coordinates of the atoms to follow
   vector<Vec> oldpos;
   
-  for(int i=0; i<at.size(); i++){
+  for(unsigned int i=0; i<at.size(); i++){
     oldpos.push_back(Vec(0.0,0.0,0.0));
   }
   // open at.size() files to write the trajectories to pdb
   vector<ofstream *> opdb(at.size());
 
-  for(int i=0; i<at.size(); i++){
+  for(unsigned int i=0; i<at.size(); i++){
     stringstream os;
     string s=molecule(at.mol(i), at.atom(i));
     s=s.substr(0,s.find(' '));
@@ -167,11 +167,11 @@ try{
     
   // print title
   cout << setw(6) << "#     ";
-  for(int i=0; i<at.size(); i++)
+  for(unsigned int i=0; i<at.size(); i++)
     cout << setw(ndim*12) << molecule(at.mol(i), at.atom(i));
   cout << endl;
   cout << setw(6) << "# time";
-  for(int i=0; i<at.size(); i++)
+  for(unsigned int i=0; i<at.size(); i++)
     for(int k=0; k<ndim; k++)
       cout << setw(12) << dimension(dim[k]);
   cout << endl;
@@ -205,7 +205,7 @@ try{
       cout << time;
       
       // loop over the atoms to consider
-      for(int i=0; i<at.size(); i++){
+      for(unsigned int i=0; i<at.size(); i++){
 	// gather with respect to its old position
 	*at.coord(i)=
 	  pbc->nearestImage(oldpos[i], *at.coord(i), sys.box());
@@ -227,7 +227,7 @@ try{
     
     ic.close();
   }
-  for(int i=0; i<at.size(); i++){
+  for(unsigned int i=0; i<at.size(); i++){
     
     writeCON(*opdb[i], count, at.size(), i);
 
diff --git a/gromos++/programs/frameout.cc b/gromos++/programs/frameout.cc
index 848d926d..2a47618b 100644
--- a/gromos++/programs/frameout.cc
+++ b/gromos++/programs/frameout.cc
@@ -149,11 +149,11 @@ int main(int argc, char **argv) {
       notimeblock = true;
 
     // get simulation time either from the user or from the files
-    bool usertime=false;
+    //bool usertime=false;
     
-    if (args.count("time") > 0) {
-      usertime=true;
-    }
+    //if (args.count("time") > 0) {
+    //  usertime=true;
+    //}
 
     // do we want to fit to a reference structure
     bool fit = false;
diff --git a/gromos++/programs/gca.cc b/gromos++/programs/gca.cc
index 2bcae660..9d6f59c9 100644
--- a/gromos++/programs/gca.cc
+++ b/gromos++/programs/gca.cc
@@ -327,7 +327,7 @@ int main(int argc, char **argv){
 int in_property(Property &p, int i)
 {
   // checks whether the atom i, is part of the property definition
-  for(int j=0; j< p.atoms().size(); j++){
+  for(unsigned int j=0; j< p.atoms().size(); j++){
     if(i==p.atoms().atom(j)) return 1;
   }
   return 0;
@@ -339,7 +339,7 @@ int in_atoms(AtomSpecifier &as, int i)
   // As our properties are always real bonds, angles or dihedrals, we do
   // not care about the molecule number here.
   // The AtomSpecifier should maybe get a function like this on its own.
-  for(int j=0; j< as.size(); j++){
+  for(unsigned int j=0; j< as.size(); j++){
     if(i==as.atom(j)) return 1;
   }
   return 0;
@@ -393,7 +393,7 @@ void move_atoms(System &sys, AtomSpecifier &as, Vec v)
   // Move the atoms in the atom-specifier by the vector v
   int m, a;
   
-  for(int i=0; i<as.size(); i++){
+  for(unsigned int i=0; i<as.size(); i++){
     m=as.mol(i);
     a=as.atom(i);
     sys.mol(m).pos(a) += v;
@@ -406,7 +406,7 @@ void rotate_atoms(System &sys, AtomSpecifier &as, gmath::Matrix rot, Vec v)
   // rotation the atoms are moved back by v again.
   int m, a;
   Vec t;
-  for(int i=0; i<as.size(); i++){
+  for(unsigned int i=0; i<as.size(); i++){
     m = as.mol(i);
     a = as.atom(i);
     t=sys.mol(m).pos(a) - v;
diff --git a/gromos++/programs/inbox.cc b/gromos++/programs/inbox.cc
index 74aa2fa2..edf730b6 100644
--- a/gromos++/programs/inbox.cc
+++ b/gromos++/programs/inbox.cc
@@ -180,7 +180,7 @@ int main(int argc, char **argv) {
         const Vec centre = (sys.box().K() * 0.5) + (sys.box().L() * 0.5) +
                 (sys.box().M() * 0.5);
 
-        for(int i = 0; i < atoms.size(); ++i) {
+        for(unsigned int i = 0; i < atoms.size(); ++i) {
           if (shift) {
             if (boxlength) {
               gmath::Vec boxshift=(sys.box().K() * shifts[0]) + (sys.box().L() * shifts[1]) +
diff --git a/gromos++/programs/int_ener.cc b/gromos++/programs/int_ener.cc
index 46ebd769..e8c2a146 100644
--- a/gromos++/programs/int_ener.cc
+++ b/gromos++/programs/int_ener.cc
@@ -275,7 +275,7 @@ int main(int argc, char **argv) {
             // calculate the energies (for atomsA)
             en.calc();
             // loop over selected atoms
-            for (int i = 0; i < atomsA.size(); ++i) {
+            for (unsigned int i = 0; i < atomsA.size(); ++i) {
               vdw += en.vdw_s(i);
               elec += en.el_s(i);
             }
@@ -284,10 +284,10 @@ int main(int argc, char **argv) {
           else {
 
             // loop over first set of selected atoms (A)
-            for (int i = 0; i < atomsA.size(); ++i) {
+            for (unsigned int i = 0; i < atomsA.size(); ++i) {
 
 	      // loop over the combined set of atoms (AB)
-              for (int j = atomsA.size(); j < atomsAB.size(); ++j) {
+              for (unsigned int j = atomsA.size(); j < atomsAB.size(); ++j) {
 
 		// set totals to zero
                 double vdw_ij = 0.0;
diff --git a/gromos++/programs/ion.cc b/gromos++/programs/ion.cc
index cb58ac89..232e1195 100644
--- a/gromos++/programs/ion.cc
+++ b/gromos++/programs/ion.cc
@@ -95,7 +95,7 @@ using namespace args;
 using namespace utils;
 
 vector<int> selectPotential(utils::AtomSpecifier &sol, 
-			    int num[2], double mindist,
+			    unsigned int num[2], double mindist,
 			    utils::Energy &en, bound::Boundary *pbc);
 
 vector<int> selectRandom(utils::AtomSpecifier &sol, 
@@ -140,7 +140,7 @@ try{
   const int positive=0;
   const int negative=1;
   
-  int num_ions[2]={0,0};
+  unsigned int num_ions[2]={0,0};
   std::string ion_names[2], res_names[2];
   {
     Arguments::const_iterator iter=args.lower_bound("positive"), 
@@ -224,7 +224,7 @@ try{
   sol.addSolventType(sys.sol(0).topology().atom(0).name());
 
   // remove the excluded solvents
-  for(int i=0; i<exclude.size(); i++)
+  for(unsigned int i=0; i<exclude.size(); i++)
     sol.removeAtom(exclude.mol(i), exclude.atom(i));
   
   if(sol.size()<num_ions[positive]+num_ions[negative])
@@ -260,7 +260,7 @@ try{
   
   // now, add the ions
   for(int t=positive; t<=negative; t++){
-    for(int i=0; i<num_ions[t]; i++){
+    for(unsigned int i=0; i<num_ions[t]; i++){
       Molecule mol(mt[t]);
       // initialize the coordinates
       mol.initPos();
@@ -325,7 +325,7 @@ try{
 
 
 vector<int> selectPotential(utils::AtomSpecifier &sol, 
-			 int num[2], double mindist, 
+			 unsigned int num[2], double mindist, 
 			 utils::Energy &en, bound::Boundary *pbc)
 {
   vector<int> selected;
@@ -340,11 +340,11 @@ vector<int> selectPotential(utils::AtomSpecifier &sol,
   
   for(int t=0; t<2; t++){
     double sign=charge*s[t];
-    for(int i=0; i<num[t]; i++){
+    for(unsigned int i=0; i<num[t]; i++){
       double min=1e6;
       int min_index=-1;
     
-      for(int j=0; j<sol.size(); j++){
+      for(unsigned int j=0; j<sol.size(); j++){
 	if(en.el(j) / sign < min && not_allowed.count(j)==0){
 	  min=en.el(j)/sign;
 	  min_index=j;
@@ -357,7 +357,7 @@ vector<int> selectPotential(utils::AtomSpecifier &sol,
 
       v1 = sol.pos(min_index);
       // now put all solvents within mindist into not_allowed.
-      for(int j=0; j<sol.size(); j++){
+      for(unsigned int j=0; j<sol.size(); j++){
 	v2 = pbc->nearestImage(sol.pos(j), v1, sol.sys()->box());
 	if((v2 - sol.pos(j)).abs2() < mindist2){
 	  not_allowed.insert(j);
@@ -377,7 +377,7 @@ vector<int> selectRandom(utils::AtomSpecifier &sol,
   double mindist2=mindist*mindist;
   gmath::Vec v1, v2;
   
-  for(int i=0; i<sol.size(); i++) allowed.push_back(i);
+  for(unsigned int i=0; i<sol.size(); i++) allowed.push_back(i);
   
   srand(seed);
   for(int i=0; i<num; i++){
@@ -388,7 +388,7 @@ vector<int> selectRandom(utils::AtomSpecifier &sol,
     selected.push_back(remove);
     gmath::Vec v1=*sol.coord(remove);
     // now remove all solvents within mindist from allowed.
-    for(int j=0; j<sol.size(); j++){
+    for(unsigned int j=0; j<sol.size(); j++){
       v2 = pbc->nearestImage(*sol.coord(j), v1, sol.sys()->box());
       if((v2 - *sol.coord(j)).abs2() < mindist2)
 	allowed.erase(allowed.begin() + index);
diff --git a/gromos++/programs/iondens.cc b/gromos++/programs/iondens.cc
index 81455797..7d5b526d 100644
--- a/gromos++/programs/iondens.cc
+++ b/gromos++/programs/iondens.cc
@@ -290,7 +290,7 @@ int main(int argc, char **argv){
 	cog=PositionUtils::cog(sys, ref);
         
 	// calculate nim with respect to cog from above
-	for (int i=0; i < ions.size(); ++i){
+	for (unsigned int i=0; i < ions.size(); ++i){
 	  ions.pos(i) = pbc->nearestImage(cog,ions.pos(i),sys.box());
 	}
 	
diff --git a/gromos++/programs/link_top.cc b/gromos++/programs/link_top.cc
index b544bbdf..769dce5c 100644
--- a/gromos++/programs/link_top.cc
+++ b/gromos++/programs/link_top.cc
@@ -199,7 +199,7 @@ int main(int argc, char *argv[]){
         if(link.atom(j).iac() < 0){
           // we will need resiue res[i][j]  
           //std::cerr << "searching for residue " << res[i][link.linkRes(j)] +1 << " and atom " << link.atom(j).name() << std::endl;
-          for(int k=0; k < lt.atoms().size(); k++){
+          for(unsigned int k=0; k < lt.atoms().size(); k++){
             if(lt.resMap()[k] == res[i][link.linkRes(j)] && lt.atoms()[k].name() == link.atom(j).name()){
               lt.atoms()[k].setIac(-1);
             }
@@ -211,10 +211,10 @@ int main(int argc, char *argv[]){
       // now we make a map of the numbers in the LinkingBlock to the current
       // atom numbers in the reduced topology
       std::map<int, int> atommap;
-      for(int j=0; j < link.numAtoms(); j++){
+      for(unsigned int j=0; j < (unsigned int)link.numAtoms(); j++){
         if(link.atom(j).iac() >= 0){
           bool found=false;
-	  for(int k=0; k < lt.atoms().size(); k++){
+	  for(unsigned int k=0; k < lt.atoms().size(); k++){
             if(lt.resMap()[k] == res[i][link.linkRes(j)] && lt.atoms()[k].name() == link.atom(j).name()){
               //std::cerr << "found atom " << link.linkRes(j) << " " << link.atom(j).name() << std::endl;
               atommap[j]=k;
@@ -452,7 +452,7 @@ int main(int argc, char *argv[]){
     os << "Modified topology based on " << args["topo"] << endl;
     os << "Applied following crosslinks according to definitions in " 
        << args["links"] << endl;
-    for(int i=0; i< name.size(); i++){
+    for(unsigned int i=0; i< name.size(); i++){
       os << "- residues " << res[i][0]+1 << " to " << res[i][1]+1 << " according to " 
          << name[i] << std::endl;
     }
diff --git a/gromos++/programs/m_widom.cc b/gromos++/programs/m_widom.cc
index 2475856b..19d0d6be 100644
--- a/gromos++/programs/m_widom.cc
+++ b/gromos++/programs/m_widom.cc
@@ -365,7 +365,7 @@ int main(int argc, char **argv) {
               rdf[j].clear();
 
               // loop over the 'with' atoms
-              for (int i = 0; i < rdfatoms[j].size(); i++) {
+              for (unsigned int i = 0; i < rdfatoms[j].size(); i++) {
                 Vec tmp;
                 tmp = pbc->nearestImage(move, *rdfatoms[j].coord(i),
                         sys.box());
diff --git a/gromos++/programs/mdf.cc b/gromos++/programs/mdf.cc
index e06e0baf..b91a558c 100644
--- a/gromos++/programs/mdf.cc
+++ b/gromos++/programs/mdf.cc
@@ -133,7 +133,7 @@ try{
   // open centre.size() files
   vector<ofstream *> fout(centre.size());
   
-  for(int i=0; i<centre.size(); i++){
+  for(unsigned int i=0; i<centre.size(); i++){
     ostringstream os;
     if(centre.mol(i)!=-3)
       os << "MIN_" << centre.mol(i)+1 << ":" << centre.atom(i)+1 
@@ -165,13 +165,13 @@ try{
       int minat=0;
       
       // now really loop over the centre atoms
-      for(int i=start;i<centre.size();i++){
+      for(unsigned int i=start;i<centre.size();i++){
 //	double min2=sys.box().K().abs()*sys.box().K().abs();
 //      The box dimension might not be the appropriate initial distance.
 //      Therefore simply use the maximum finite representable floating-point number:
         double min2=DBL_MAX;
 	//loop over the atoms to consider
-        for(int j=0;j<with.size();j++){
+        for(unsigned int j=0;j<with.size();j++){
           //calculate distance only if this atom is not the current centre
           if(!(with.mol(j)==centre.mol(i)&&with.atom(j)==centre.atom(i))){
 	    Vec tmp=pbc->nearestImage(centre.pos(i),
@@ -196,7 +196,7 @@ try{
   }
   
   //close the files
-  for(int i=0;i<centre.size();i++){
+  for(unsigned int i=0;i<centre.size();i++){
     fout[i]->close();
     delete fout[i];
     fout[i] = NULL;
diff --git a/gromos++/programs/nhoparam.cc b/gromos++/programs/nhoparam.cc
index 4a2f1e69..0ba10f12 100644
--- a/gromos++/programs/nhoparam.cc
+++ b/gromos++/programs/nhoparam.cc
@@ -209,7 +209,7 @@ int main(int argc, char **argv) {
 
     //see whether no atoms are supplied
     vector<set<int> > hydrogens(atoms.size(), set<int>());
-    for (int i = 0; i < atoms.size(); ++i) {
+    for (unsigned int i = 0; i < atoms.size(); ++i) {
       if (atoms.mass(i) != nitrogen_mass) {
         ostringstream msg;
         msg << "Atom " << atoms.toString(i) << " is not a nitrogen.";
@@ -254,7 +254,7 @@ int main(int argc, char **argv) {
     // print gromos numbers to title
     ts << "#" << setw(14) << "time";
     tsw << "#" << setw(14) << "time";
-    for (int i = 0; i < atoms.size(); ++i) {
+    for (unsigned int i = 0; i < atoms.size(); ++i) {
       ts << setw(10) << (atoms.gromosAtom(i) + 1);
       tsw << setw(10) << (atoms.gromosAtom(i) + 1);
     }
@@ -283,7 +283,7 @@ int main(int argc, char **argv) {
           tsw << time;
 
         // calculate the z-vector between atom i-1 and i+1, normalize
-        for (int i = 0; i < atoms.size(); i++) {
+        for (unsigned int i = 0; i < atoms.size(); i++) {
           // holds the N-H bond
           gmath::Vec nh; 
           set<int>::const_iterator it = hydrogens[i].begin();
@@ -372,7 +372,7 @@ int main(int argc, char **argv) {
             << setw(8) << "RESNAME" << setw(10) << "S2" << setw(10) << "WINAV"
             << setw(10) << "WINRMSD" << setw(10) << "WINEE" << endl;
 
-    for (int i = 0; i < atoms.size(); ++i) {
+    for (unsigned int i = 0; i < atoms.size(); ++i) {
       cout.setf(ios::floatfield, ios::fixed);
       cout.setf(ios::right, ios::adjustfield);
       cout.precision(4);
diff --git a/gromos++/programs/pairlist.cc b/gromos++/programs/pairlist.cc
index a7fae55e..f0bef56d 100644
--- a/gromos++/programs/pairlist.cc
+++ b/gromos++/programs/pairlist.cc
@@ -210,7 +210,7 @@ int main(int argc, char **argv){
 	pl.addAtom(ref_as.mol(0), ref_as.atom(0));
     }
     if(do_diff){
-      for(int i=0; i<pl.size(); i++)
+      for(unsigned int i=0; i<pl.size(); i++)
 	if(pp.findAtom(pl.mol(i), pl.atom(i))==-1)
 	  diff.addAtom(pl.mol(i), pl.atom(i));
     }
@@ -257,7 +257,7 @@ int main(int argc, char **argv){
 	     << setw(12) << "Atom Code"
 	     << endl;
 	
-	for(int i=0; i < pp.size(); ++i){
+	for(unsigned int i=0; i < pp.size(); ++i){
 	  cout << setw(13) << pp.toString(i)
 	       << setw(10) << pp.gromosAtom(i)+1
 	       << setw(10) << pp.resnum(i)+1
@@ -305,7 +305,7 @@ int main(int argc, char **argv){
 	     << setw(12) << "Atom Code"
 	     << endl;
 
-	for(int i=0; i < pl.size(); ++i){
+	for(unsigned int i=0; i < pl.size(); ++i){
 	  cout << setw(13) << pl.toString(i)
 	       << setw(10) << pl.gromosAtom(i)+1
 	       << setw(10) << pl.resnum(i)+1
@@ -354,7 +354,7 @@ int main(int argc, char **argv){
 	     << endl;
 
 
-	for(int i=0; i < diff.size(); ++i){
+	for(unsigned int i=0; i < diff.size(); ++i){
 	  cout << setw(13) << diff.toString(i)
 	       << setw(10) << diff.gromosAtom(i)+1
 	       << setw(10) << diff.resnum(i)+1
diff --git a/gromos++/programs/pert_top.cc b/gromos++/programs/pert_top.cc
index c6fab4a6..7bcadff4 100644
--- a/gromos++/programs/pert_top.cc
+++ b/gromos++/programs/pert_top.cc
@@ -154,11 +154,11 @@ int main(int argc, char *argv[]){
     // set the title
     ostringstream title;
     title << "Perturbation of atoms: ";
-    for (int i = 0; i < at.size(); i++) title << at.toString(i) << " ";
+    for (unsigned int i = 0; i < at.size(); i++) title << at.toString(i) << " ";
     title << endl << "From topology: " << args["topo"];
     out_pt.setTitle(title.str());
     
-    for (int i = 0; i < at.size(); i++) {
+    for (unsigned int i = 0; i < at.size(); i++) {
       // set state A
       pttop.setAtomNum(i, at.gromosAtom(i));
       pttop.setAtomName(i, at.name(i));
@@ -175,7 +175,7 @@ int main(int argc, char *argv[]){
         iac = at.iac(i);
       } else {
         // take the right value or the last one
-        iac = i < int(iacs.size()) ? iacs[i] : iacs.back();
+        iac = i < (unsigned int)(iacs.size()) ? iacs[i] : iacs.back();
       }
       pttop.setIac(i, 1, iac);
       
@@ -185,7 +185,7 @@ int main(int argc, char *argv[]){
         mass = at.mass(i);
       } else {
         // take the right value or the last one
-        mass = i < int(masses.size()) ? masses[i] : masses.back();
+        mass = i < (unsigned int)(masses.size()) ? masses[i] : masses.back();
       }
       pttop.setMass(i, 1, mass);
       
@@ -195,7 +195,7 @@ int main(int argc, char *argv[]){
         charge = at.charge(i);
       } else {
         // take the right value or the last one
-        charge = i < int(charges.size()) ? charges[i] : charges.back();
+        charge = i < (unsigned int)(charges.size()) ? charges[i] : charges.back();
       }
       pttop.setCharge(i, 1, charge);
     }
diff --git a/gromos++/programs/predict_noe.cc b/gromos++/programs/predict_noe.cc
index 610b4a7f..3fc4cdff 100644
--- a/gromos++/programs/predict_noe.cc
+++ b/gromos++/programs/predict_noe.cc
@@ -156,7 +156,7 @@ int main(int argc, char *argv[]) {
     
     vector<VirtualAtom*> atoms;
     vector<pair<string, unsigned int> > atom_description;
-    for(int i = 0; i < specatoms.size(); ++i) {
+    for(unsigned int i = 0; i < specatoms.size(); ++i) {
       vector<Noelib>::const_iterator it = noelib.begin(),
               to = noelib.end();
       for(; it != to; ++it) {
@@ -245,4 +245,4 @@ int main(int argc, char *argv[]) {
     exit(1);
   }
   return 0;
-}
\ No newline at end of file
+}
diff --git a/gromos++/programs/prep_noe.cc b/gromos++/programs/prep_noe.cc
index 2cde2a2b..c73040e8 100644
--- a/gromos++/programs/prep_noe.cc
+++ b/gromos++/programs/prep_noe.cc
@@ -306,7 +306,7 @@ int main(int argc, char *argv[]) {
       // to have it in the file for comparision with the output of post_noe
       int a = atoi(tokens[1].c_str());
       int b = atoi(tokens[3].c_str());
-      int noe_number = atoi(tokens[0].c_str());
+      unsigned int noe_number = atoi(tokens[0].c_str());
       double d = atof(tokens[5].c_str());
       double e = atof(tokens[6].c_str());
       double f = atof(tokens[7].c_str());
@@ -319,7 +319,7 @@ int main(int argc, char *argv[]) {
       }
       // the noe numbers must be consecutive due to downstream programs
       if (!noevec.empty()) {
-          int last_noe_num = noevec.rbegin()->first;
+          unsigned int last_noe_num = noevec.rbegin()->first;
           if (noe_number <= last_noe_num){ // check if the last element of this ordered map is bigger than the current noe number
             stringstream msg;
             msg <<"NOE number " << noe_number << " is smaller than its predecessor " << last_noe_num << ". This is not allowed";
@@ -676,7 +676,7 @@ int main(int argc, char *argv[]) {
         int mol = VA.conf().mol(0);
         for (int l = 0; l < mol; ++l) offsetA += sys.mol(l).numAtoms();
 
-        for (int aa = 0; aa < 4; ++aa) {
+        for (unsigned int aa = 0; aa < 4; ++aa) {
           int att;
           if (VA.conf().size() > aa)
             att = VA.conf().atom(aa);
@@ -692,7 +692,7 @@ int main(int argc, char *argv[]) {
           int mol = VB.conf().mol(0);
           for (int l = 0; l < mol; ++l) offsetB += sys.mol(l).numAtoms();
 
-          for (int bb = 0; bb < 4; ++bb) {
+          for (unsigned int bb = 0; bb < 4; ++bb) {
             int att;
             if (VB.conf().size() > bb)
               att = VB.conf().atom(bb);
diff --git a/gromos++/programs/prep_xray_le.cc b/gromos++/programs/prep_xray_le.cc
index 31e81ac8..98d5d26a 100644
--- a/gromos++/programs/prep_xray_le.cc
+++ b/gromos++/programs/prep_xray_le.cc
@@ -129,7 +129,7 @@ int main(int argc, char *argv[]) {
     atoms.sort();
     { // loop over the atoms and include the first atom of the residue
       set<pair<int,int> > res;
-      for(int i = 0; i < atoms.size(); ++i)
+      for(unsigned int i = 0; i < atoms.size(); ++i)
         res.insert(pair<int,int>(atoms.mol(i), atoms.resnum(i)));
       atoms.clear();
       for(set<pair<int,int> >::const_iterator it = res.begin(), to = res.end();
@@ -217,7 +217,7 @@ int main(int argc, char *argv[]) {
     vector<string> res_desc;
     vector<vector<int> > res_indices;
     vector<set<int> > res_unique_indices;
-    for (int r = 0; r < atoms.size(); ++r) {
+    for (unsigned int r = 0; r < atoms.size(); ++r) {
       string resname = atoms.resname(r);
       // skip residue if not present in lib
       if (sidechain_spec.find(resname) == sidechain_spec.end())
diff --git a/gromos++/programs/ran_box.cc b/gromos++/programs/ran_box.cc
index d0126043..e768a8c7 100644
--- a/gromos++/programs/ran_box.cc
+++ b/gromos++/programs/ran_box.cc
@@ -118,7 +118,7 @@ int main(int argc, char **argv) {
     Arguments args(argc, argv, knowns, usage);
 
     // define some variables
-    const double pi = acos(-1.0);
+    //const double pi = acos(-1.0);
     const double nano_i = 1 / gmath::physConst.get_nano();
     // 1.66054 * 10^-27 * 10^27 = 1.66054
     const double fac_amu2kg = gmath::physConst.get_atomic_mass_unit() * nano_i * nano_i * nano_i ;
diff --git a/gromos++/programs/ran_solvation.cc b/gromos++/programs/ran_solvation.cc
index 65962ef2..e76d2c53 100644
--- a/gromos++/programs/ran_solvation.cc
+++ b/gromos++/programs/ran_solvation.cc
@@ -204,7 +204,7 @@ int main(int argc, char **argv) {
     //    solute is rotated and the specified numbers are added to the
     //    maximum distances in the resulting dimensions
     vector<double> minwall;
-    int calculate_dimensions = 0;
+    //int calculate_dimensions = 0;
     {
       Arguments::const_iterator iter = args.lower_bound("minwall"),
               to = args.upper_bound("minwall");
@@ -212,7 +212,7 @@ int main(int argc, char **argv) {
         minwall.push_back(atof(iter->second.c_str()));
         ++iter;
       }
-      if (minwall.size()) calculate_dimensions = 1;
+      //if (minwall.size()) calculate_dimensions = 1;
     }
 
     // read the boxsize
@@ -241,13 +241,13 @@ int main(int argc, char **argv) {
     System solu(it.system());
     Box box;
     Boundary *pbc = BoundaryParser::boundary(solu, args);
-    int truncoct = 0, rectbox = 0, cubic = 0;
+    //int truncoct = 0, rectbox = 0, cubic = 0;
     double size_corr = 1.0;
     double fac_vol = 1.0;
     switch (pbc->type()) {
       case('t'):
-        truncoct = 1;
-        rectbox = 0;
+        //truncoct = 1;
+        //rectbox = 0;
         size_corr = 2.0 * sqrt(3.0) / 3.0;
         fac_vol = 0.5;
         if (minwall.size() > 1)
@@ -260,13 +260,13 @@ int main(int argc, char **argv) {
         box.setNtb(gcore::Box::truncoct);
         break;
       case('r'):
-        truncoct = 0;
-        rectbox = 1;
+        //truncoct = 0;
+        //rectbox = 1;
         fac_vol = 1.0;
-        if (minwall.size() == 1) cubic = 1;
-        else if (minwall.size() == 0 &&
-                (boxsize[0] == boxsize[1] && boxsize[0] == boxsize[2]))
-          cubic = 1;
+        //if (minwall.size() == 1) cubic = 1;
+        //else if (minwall.size() == 0 &&
+        //        (boxsize[0] == boxsize[1] && boxsize[0] == boxsize[2]))
+          //cubic = 1;
         box.setNtb(gcore::Box::rectangular);
         break;
       case('v'):
diff --git a/gromos++/programs/rep_ana.cc b/gromos++/programs/rep_ana.cc
index 88d08076..466a3867 100644
--- a/gromos++/programs/rep_ana.cc
+++ b/gromos++/programs/rep_ana.cc
@@ -177,7 +177,7 @@ int main(int argc, char *argv[]) {
         first_run = r.run;  // so this also works when user used equilibration steps and the first run != 1
       }
 
-      if (r.run == first_run) {
+      if (int(r.run) == first_run) {
         // lambda
         if (r.s == 1) { // exchange successfull
           lam_series[r.ID - 1].push_back(r.li);
diff --git a/gromos++/programs/rgyr.cc b/gromos++/programs/rgyr.cc
index eedc391f..95b27bf0 100644
--- a/gromos++/programs/rgyr.cc
+++ b/gromos++/programs/rgyr.cc
@@ -126,7 +126,7 @@ int main(int argc, char **argv){
 
     // calculate total mass
     double totalMass=0.0;
-    for(int i=0; i< atom.size(); ++i)
+    for(unsigned int i=0; i< atom.size(); ++i)
       totalMass+=atom.mass(i);
     
     bool massweighted = false;
@@ -151,7 +151,7 @@ int main(int argc, char **argv){
 	
 	//calculate cm, rgyr
 	Vec cm(0.0,0.0,0.0);
-	for (int i=0;i < atom.size(); i++) {
+	for (unsigned int i=0;i < atom.size(); i++) {
 	  cm += atom.pos(i) * atom.mass(i);
 	}
 	cm /= totalMass;
@@ -159,13 +159,13 @@ int main(int argc, char **argv){
 	double rg=0; 
 
 	if(massweighted) {
-	  for(int i=0; i < atom.size(); ++i){
+	  for(unsigned int i=0; i < atom.size(); ++i){
 	    rg += atom.mass(i)*(atom.pos(i) - cm).abs2();
 	  }
 	  rg = sqrt(rg/totalMass);
 	}
 	else {
-	  for (int i=0;i < atom.size(); i++) {
+	  for (unsigned int i=0;i < atom.size(); i++) {
 	    // should we correct for periodicity?
 	    // Only if atom contains different molecules.
 	    // But then cm would be wrong already. The user should just use 
diff --git a/gromos++/programs/rmsd.cc b/gromos++/programs/rmsd.cc
index 7746e3d1..e205de00 100644
--- a/gromos++/programs/rmsd.cc
+++ b/gromos++/programs/rmsd.cc
@@ -115,7 +115,7 @@ double rmsdproperty(const utils::PropertyContainer &prop_ref,
 int main(int argc, char **argv){
   Argument_List knowns; 
   knowns << "topo" << "traj" << "atomsfit" << "atomsrmsd" << "prop" << "pbc" << "ref"
-         << "time" << "debug" << "fit" << "reftopo" << "refpbc" << "printatoms";
+         << "time"  << "reftopo" << "refpbc" << "printatoms";
 
   string usage = "# " + string(argv[0]);
   usage += "\n\t@topo       <molecular topology file>\n";
@@ -190,15 +190,6 @@ int main(int argc, char **argv){
     } else {
       cerr << "# Warning: @reftopo but no @refpbc : reference will not be gathered\n";
     }
-    
-
-    int debug=0;
-    if(args.count("debug")>0)
-        debug=1;
-    int fit=0;
-    if(args.count("fit")>0)
-        //fit=1;
-        fit=atoi(args.lower_bound("fit")->second.c_str());
 
     // this always goes wrong. check that there is a box block in the refSys
     //if(refSys.hasBox == false && pbc->type()!='v')
@@ -247,11 +238,11 @@ int main(int argc, char **argv){
     }
     
     // optionally print atomspecs to check them
-    int maxsize=0;
+    unsigned int maxsize=0;
     if (fitatoms.size() > rmsdatoms.size()) maxsize=fitatoms.size();
     else maxsize=rmsdatoms.size();
     if (printatoms) {
-        for (int i=0; i < maxsize; i++) {
+        for (unsigned int i=0; i < maxsize; i++) {
               cerr << "# " ;
             if (i<fitatoms.size()) {
               cerr << fitatoms.toString(i) << " "<< fitatoms.mol(i)+1 <<":"<< fitatoms.resname(i) << fitatoms.resnum(i)+1<<":" << fitatoms.name(i);
@@ -328,7 +319,7 @@ int main(int argc, char **argv){
         Vec cogfitref (0,0,0);
         Vec cogfitsys (0,0,0);
         if (fitatomsref.size()){
-          for(int i=0;i<fitatomsref.size();++i) {
+          for(unsigned int i=0;i<fitatomsref.size();++i) {
             cogfitref+=fitatomsref.pos(i);
             cogfitsys+=fitatoms.pos(i);
           }
@@ -341,11 +332,11 @@ int main(int argc, char **argv){
         std::vector<gmath::Vec> rmsdref;
         std::vector<gmath::Vec> rmsdsys;
         if (fitatomsref.size()){
-          for(int i=0;i<fitatomsref.size();++i) {
+          for(unsigned int i=0;i<fitatomsref.size();++i) {
             fitref.push_back(fitatomsref.pos(i)-cogfitref);
             fitsys.push_back(fitatoms.pos(i)-cogfitsys);
           }
-          for(int i=0;i<rmsdatomsref.size();++i) {
+          for(unsigned int i=0;i<rmsdatomsref.size();++i) {
             rmsdref.push_back(rmsdatomsref.pos(i)-cogfitref);
             rmsdsys.push_back(rmsdatoms.pos(i)-cogfitsys);
           }
diff --git a/gromos++/programs/rmsdmat.cc b/gromos++/programs/rmsdmat.cc
index dfc26002..337d9803 100644
--- a/gromos++/programs/rmsdmat.cc
+++ b/gromos++/programs/rmsdmat.cc
@@ -119,7 +119,7 @@ using namespace fit;
 
 double props_rmsd(const vector<vector<Value> > &props, int i, int j) {
   double rmsd = 0.0;
-  for (int k=0; k< props[i].size(); k++) {
+  for (unsigned int k=0; k< props[i].size(); k++) {
     // get nearestimagedistance for dihedrals+angles
     // assuming distances or Hbonds will never be >|180|, they should be ok
     double diff = props[j][k].scalar() - props[i][k].scalar(); 
@@ -141,8 +141,8 @@ void update_bins(map<int, int > & bins, double bin_size, double x) {
 }
 
 void write_bins(map<int, int > bins, double bin_size) {
-  int maxbin=bins.rbegin()->first;
-  int minbin=bins.begin()->first;
+  //int maxbin=bins.rbegin()->first;
+  //int minbin=bins.begin()->first;
   ofstream distfile;
   distfile.open("distr.out");
   distfile << "# rmsd distribution, bin size " <<bin_size <<" nm" << endl;
@@ -318,10 +318,10 @@ int main(int argc, char **argv) {
     // get atom specifications for refsys from args
     if (do_atomrmsd) 
         getAtomSpecs(refSys, args, atomspecs, fitatoms, rmsdatoms);
-    int atomspecnum=atomspecs.size();
+    unsigned int atomspecnum=atomspecs.size();
     
     if (printatoms) {
-        for (int i=0; i < atomspecs.size(); i++) {
+        for (unsigned int i=0; i < atomspecs.size(); i++) {
             cout << atomspecs.toString(i) << " "<< atomspecs.mol(i)+1 <<":"<< atomspecs.resname(i) << atomspecs.resnum(i)+1<<":" << atomspecs.name(i) <<  endl;
         }
     }
@@ -331,7 +331,7 @@ int main(int argc, char **argv) {
     if (atomspecs.size() != fitatoms.size() || atomspecs.size() != rmsdatoms.size()) {
       fit_spec.resize(atomspecs.size(), true);
       rmsd_spec.resize(atomspecs.size(), true);
-      for (int i = 0; i < atomspecs.size(); ++i) {
+      for (unsigned int i = 0; i < atomspecs.size(); ++i) {
         if (fitatoms.findAtom(atomspecs.mol(i), atomspecs.atom(i)) < 0)
           fit_spec[i] = false;
         if (rmsdatoms.findAtom(atomspecs.mol(i), atomspecs.atom(i)) < 0)
@@ -367,14 +367,14 @@ int main(int argc, char **argv) {
 
     // calculate the centre of geometry of the relevant atoms
     Vec cog;
-    for (int i = 0; i < fitatoms.size(); ++i) {
+    for (unsigned int i = 0; i < fitatoms.size(); ++i) {
       cog += *fitatoms.coord(i);
     }
     cog /= fitatoms.size();
 
     // put it in the trajectory
     vector< Vec > frame(atomspecs.size());
-    for (int i = 0; i < atomspecs.size(); ++i) {
+    for (unsigned int i = 0; i < atomspecs.size(); ++i) {
       frame[i] = *atomspecs.coord(i) - cog;
     }
     traj.push_back(frame);
@@ -414,7 +414,7 @@ int main(int argc, char **argv) {
     if (do_atomrmsd)
       getAtomSpecs(sys, args, atomspecs, fitatoms, rmsdatoms);
     if (printatoms) {
-        for (int i=0; i < atomspecs.size(); i++) {
+        for (unsigned int i=0; i < atomspecs.size(); i++) {
             cout << atomspecs.toString(i) << " "<< atomspecs.mol(i)+1 <<":"
                  << atomspecs.resname(i) << atomspecs.resnum(i)+1<<":" 
                  << atomspecs.name(i) <<  endl;
@@ -451,7 +451,7 @@ int main(int argc, char **argv) {
         if (do_atomrmsd) getAtomSpecs(sys, args, atomspecs, fitatoms, rmsdatoms);
         
         if (printatoms) {
-          for (int i=0; i < atomspecs.size(); i++) {
+          for (unsigned int i=0; i < atomspecs.size(); i++) {
             cout << atomspecs.toString(i) << " "<< atomspecs.mol(i)+1 <<":"<< atomspecs.resname(i) << atomspecs.resnum(i)+1<<":" << atomspecs.name(i) <<  endl;
           }
         }
@@ -488,11 +488,11 @@ int main(int argc, char **argv) {
             propcnt++;
           } else {
             Vec cog;
-            for (int i = 0; i < fitatoms.size(); ++i) {
+            for (unsigned int i = 0; i < fitatoms.size(); ++i) {
               cog += *fitatoms.coord(i);
             }
             cog /= fitatoms.size();
-            for (int i = 0; i < atomspecs.size(); ++i) {
+            for (unsigned int i = 0; i < atomspecs.size(); ++i) {
               frame[i] = *atomspecs.coord(i) - cog;
             }
             int err = frf.fit(traj[0], frame);
@@ -608,8 +608,8 @@ int main(int argc, char **argv) {
 
               rmsd = frf.rmsd(rot, traj[i], traj[j]);
             }
-          if (do_dist)
-            update_bins(bins, bin_size, rmsd);
+            if (do_dist)
+              update_bins(bins, bin_size, rmsd);
 
             rmsd *= precision;
             if (rmsd > std::numeric_limits<unsigned short>::max()) {
@@ -638,8 +638,8 @@ int main(int argc, char **argv) {
 
               rmsd = frf.rmsd(rot, traj[i], traj[j]);
             }
-          if (do_dist)
-            update_bins(bins, bin_size, rmsd);
+            if (do_dist)
+              update_bins(bins, bin_size, rmsd);
 
             rmsd *= precision;
             if (rmsd > std::numeric_limits<unsigned int>::max())
diff --git a/gromos++/programs/rmsf.cc b/gromos++/programs/rmsf.cc
index 89fa8aaa..5e589491 100644
--- a/gromos++/programs/rmsf.cc
+++ b/gromos++/programs/rmsf.cc
@@ -221,7 +221,7 @@ int main(int argc, char **argv){
           apos.resize(rmsfatoms.size(), Vec(0.0, 0.0, 0.0));
           apos2.resize(rmsfatoms.size(), 0.0);
           rmsf.resize(rmsfatoms.size(), 0.0);
-          for(int i = 0; i < rmsfatoms.size(); ++i) {
+          for(unsigned int i = 0; i < rmsfatoms.size(); ++i) {
             firstpos.push_back(rmsfatoms.pos(i));
           }
         }
@@ -233,7 +233,7 @@ int main(int argc, char **argv){
 	  rf->fit(&sys);
 	
 	// calculate <r> and <r^2>
-	for(int i=0; i< rmsfatoms.size(); ++i){
+	for(unsigned int i=0; i< rmsfatoms.size(); ++i){
           const Vec & gathpos = rmsfatoms.pos(i); // pbc->nearestImage(firstpos[i], rmsfatoms.pos(i), sys.box());
 	  apos[i] += gathpos;
 	  apos2[i] += gathpos.abs2();
@@ -244,7 +244,7 @@ int main(int argc, char **argv){
     } //end loop over trajectory
     
     // calculate the rmsf's
-    for(int i=0; i < rmsfatoms.size(); ++i){
+    for(unsigned int i=0; i < rmsfatoms.size(); ++i){
       apos2[i]/=numFrames;
       apos[i]/=numFrames;
       
@@ -254,7 +254,7 @@ int main(int argc, char **argv){
     //spit out results
     cout << "#\n#  at          rmsf name\n";
      
-    for (int i=0; i < rmsfatoms.size(); ++i) {
+    for (unsigned int i=0; i < rmsfatoms.size(); ++i) {
       cout.precision(8);
       cout << setw(5) << i+1
 	   << setw(14) << rmsf[i]
@@ -270,7 +270,7 @@ int main(int argc, char **argv){
     
     if (writepdb) {
     // get a system
-      for(int i=0; i< rmsfatoms.size(); ++i){
+      for(unsigned int i=0; i< rmsfatoms.size(); ++i){
         rmsfatoms.pos(i) = apos[i];
       }
       ofstream fout(outpdb.c_str());
diff --git a/gromos++/programs/sasa.cc b/gromos++/programs/sasa.cc
index ee74f065..dddea08b 100644
--- a/gromos++/programs/sasa.cc
+++ b/gromos++/programs/sasa.cc
@@ -176,7 +176,7 @@ int main(int argc, char **argv) {
         atoms.addSpecifier(spec);
       }
       // check that they are all in 'roll atoms'
-      for(int i=0; i < atoms.size(); ++i){
+      for(unsigned int i=0; i < atoms.size(); ++i){
 	if(rollatoms.findAtom(atoms.mol(i), atoms.atom(i))==-1){
 	  ostringstream os;
 	  os << "Error: not all selected atoms are part of the sasaatoms: " << atoms.toString(i);
@@ -201,7 +201,7 @@ int main(int argc, char **argv) {
     utils::compute_atomic_radii_vdw(probe_iac, probe, sys, it.forceField());
 
     //get all heavy atoms...
-    for(int i = 0; i < rollatoms.size(); ++i) {
+    for(unsigned int i = 0; i < rollatoms.size(); ++i) {
       if (!sys.mol(rollatoms.mol(i)).topology().atom(rollatoms.atom(i)).isH()) {
 	heavyatoms.addAtom(rollatoms.mol(i), rollatoms.atom(i));
 	double rad = rollatoms.radius(i);
@@ -259,7 +259,7 @@ int main(int argc, char **argv) {
 	  
 	  utils::AtomSpecifier pairlist(sys);
 	  std::vector<int> neighbour;
-	  for(int i=0; i < heavyatoms.size(); ++i){
+	  for(unsigned int i=0; i < heavyatoms.size(); ++i){
 	    if(i!=ir){
 		
 	      Vec dist=heavyatoms.pos(i) - heavyatoms.pos(ir);
@@ -325,7 +325,7 @@ int main(int argc, char **argv) {
               vector<double> arci;
 
               // inner loop over neighbors
-              for (int j = 0; j < pairlist.size(); ++j) {
+              for (unsigned int j = 0; j < pairlist.size(); ++j) {
                 
 		// calculate some distances
 		Vec tmp = heavyatoms.pos(neighbour[j]);
@@ -502,7 +502,7 @@ int main(int argc, char **argv) {
 
       double sumaccs = 0.0, sumper = 0.0, sumpera = 0.0;
 
-      for (int i = 0; i < atoms.size(); ++i) {
+      for (unsigned int i = 0; i < atoms.size(); ++i) {
         int index = heavyatoms.findAtom(atoms.mol(i), atoms.atom(i));
         if (index >= 0) {
           cout << "# "
diff --git a/gromos++/programs/shake_analysis.cc b/gromos++/programs/shake_analysis.cc
index 63e36eec..9f8ed606 100644
--- a/gromos++/programs/shake_analysis.cc
+++ b/gromos++/programs/shake_analysis.cc
@@ -241,7 +241,7 @@ int main(int argc, char **argv) {
                 << endl;
 
         cout << "Covalent interactions involving atoms ";
-        for (int i = 0; i < atoms.size(); i++) {
+        for (unsigned int i = 0; i < atoms.size(); i++) {
           if (atoms.mol(i) < 0) cout << "s";
           else cout << atoms.mol(i) + 1;
           cout << ":" << atoms.atom(i) + 1 << " ";
@@ -418,7 +418,7 @@ int main(int argc, char **argv) {
     // double max_e_vdw = 0.0, max_e_crf = 0.0, max_e_tot = 0.0;
 
     // loop over the relevant atoms
-    for (int i = 0; i < atoms.size(); i++) {
+    for (unsigned int i = 0; i < atoms.size(); i++) {
 
       //make a pairlist
       SimplePairlist pl(sys, *pbc, cut);
diff --git a/gromos++/programs/solute_entropy.cc b/gromos++/programs/solute_entropy.cc
index e06972c6..c75eaed0 100644
--- a/gromos++/programs/solute_entropy.cc
+++ b/gromos++/programs/solute_entropy.cc
@@ -175,7 +175,7 @@ int main(int argc, char **argv) {
                     string spec=iter->second.c_str();
                     gathlist.addSpecifierStrict(spec);
                 }
-                for(int j=0;j<gathlist.size()/2;++j){
+                for(unsigned int j=0;j<gathlist.size()/2;++j){
                     int i=2*j;
                     sys.primlist[gathlist.mol(i)][0]=gathlist.atom(i);
                     sys.primlist[gathlist.mol(i)][1]=gathlist.mol(i+1);
@@ -486,7 +486,7 @@ double freq_etc(gsl_matrix * cov, gsl_vector * av, gsl_vector * mass,
   // define necessary constants
   //Boltzmann constant in J/K
   const double KB = gmath::physConst.get_boltzmann() * 1000 / gmath::physConst.get_avogadro();
-  const double E = gmath::physConst.get_euler(); /// Euler number
+  //const double E = gmath::physConst.get_euler(); /// Euler number
   const double HBAR = gmath::physConst.get_hbar() / (1e9 * gmath::physConst.get_avogadro()); /// Plank constant over 2 Pi, SI units
   const double MU = gmath::physConst.get_atomic_mass_unit(); /// Atomic mass unit
   const double NA = gmath::physConst.get_avogadro(); /// Avogadros number
diff --git a/gromos++/programs/tcf.cc b/gromos++/programs/tcf.cc
index 6e5b3912..9d8cce3e 100644
--- a/gromos++/programs/tcf.cc
+++ b/gromos++/programs/tcf.cc
@@ -436,7 +436,7 @@ int main(int argc, char **argv) {
       }
       cout << "\n#        t          C(t)\n";
 
-      for (int i = 0; i < corr->size(); i++, time += dt) {
+      for (unsigned int i = 0; i < corr->size(); i++, time += dt) {
         cout << setw(10) << time
                 << setw(14) << (*corr)[i] << "\n";
       }
diff --git a/gromos++/src/args/BoundaryParser.cc b/gromos++/src/args/BoundaryParser.cc
index 9f41ac56..6893cc1a 100644
--- a/gromos++/src/args/BoundaryParser.cc
+++ b/gromos++/src/args/BoundaryParser.cc
@@ -99,7 +99,7 @@ bound::Boundary *BoundaryParser::boundary(gcore::System &sys,
       ++it;
         std::vector<int> molecules;
         utils::parse_range(it->second.c_str(), molecules);
-        for (int i=0; i<molecules.size(); i++) {
+        for (unsigned int i=0; i<molecules.size(); i++) {
           pbc->addRefMol(molecules[i]);
           //std::cout << "# Gathering molecule " << molecules[i] << " to reference" << std::endl;
         }
diff --git a/gromos++/src/args/GatherParser.cc b/gromos++/src/args/GatherParser.cc
index e75434de..889117cd 100644
--- a/gromos++/src/args/GatherParser.cc
+++ b/gromos++/src/args/GatherParser.cc
@@ -158,7 +158,7 @@ bound::Boundary::MemPtr GatherParser::parse(gcore::System &sys, gcore::System &r
       if (uselist) {
         if (list_available) {
             //the block of primlist
-              for (int j = 0; j < gathlist.size() / 2; ++j) {
+              for (unsigned int j = 0; j < gathlist.size() / 2; ++j) {
                 int i = 2 * j;
                 sys.primlist[gathlist.mol(i)][0] = gathlist.atom(i);
                 sys.primlist[gathlist.mol(i)][1] = gathlist.mol(i + 1);
diff --git a/gromos++/src/bound/Boundary.cc b/gromos++/src/bound/Boundary.cc
index dfa3184a..cb31152b 100644
--- a/gromos++/src/bound/Boundary.cc
+++ b/gromos++/src/bound/Boundary.cc
@@ -150,6 +150,10 @@ void Boundary::gatherlist() {
     throw gromos::Exception("Gather problem",
           "Box block contains element(s) of value 0.0! Abort!");
 
+  if (sys().numMolecules() < 1)
+    throw gromos::Exception("Gather problem",
+          "the glist gather method requires at least one solute molecule");
+
   // gather the first molecule
   Molecule &mol = sys().mol(0);
   mol.pos(0) = nearestImage(reference(0), mol.pos(0), sys().box());
@@ -278,6 +282,10 @@ void Boundary::gatherltime() {
     throw gromos::Exception("Gather problem",
           "Box block contains element(s) of value 0.0! Abort!");
 
+  if (sys().numMolecules() < 1)
+    throw gromos::Exception("Gather problem",
+          "the gltime gather method requires at least one solute molecule");
+
   if (sys().primlist[0][0] == 31415926) {
     for (int i = 0; i < sys().numMolecules(); ++i) {
       Molecule &mol = sys().mol(i);
@@ -657,6 +665,10 @@ void Boundary::coggather() {
     throw gromos::Exception("Gather problem",
           "Box block contains element(s) of value 0.0! Abort!");
 
+  if (sys().numMolecules() < 1)
+    throw gromos::Exception("Gather problem",
+          "the cog gather method requires at least one solute molecule");
+
   Molecule &mol = sys().mol(0);
   Solvent &sol = sys().sol(0);
 
@@ -747,7 +759,7 @@ void Boundary::gfitgather() {
   
   Vec cog(0.0, 0.0, 0.0);
   int atoms = 0;
-  for (int i=0; i<d_refmol.size(); i++) {
+  for (unsigned int i=0; i<d_refmol.size(); i++) {
     int molnum=d_refmol[i];
 
     Molecule &mol = sys().mol(molnum);
diff --git a/gromos++/src/fit/FastRotationalFit.cc b/gromos++/src/fit/FastRotationalFit.cc
index f6136e54..d1085660 100644
--- a/gromos++/src/fit/FastRotationalFit.cc
+++ b/gromos++/src/fit/FastRotationalFit.cc
@@ -81,7 +81,7 @@ int FastRotationalFit::fit(utils::AtomSpecifier & ref_spec,
   std::vector<Vec> v_ref(ref_spec.size()), 
     v_sys(sys_spec.size());
 
-  for(int i=0; i<ref_spec.size(); ++i){
+  for(unsigned int i=0; i<ref_spec.size(); ++i){
     v_ref[i] = ref_spec.pos(i);
     v_sys[i] = sys_spec.pos(i);
   }
diff --git a/gromos++/src/fit/PositionUtils.cc b/gromos++/src/fit/PositionUtils.cc
index 001b855c..0d7eb2eb 100644
--- a/gromos++/src/fit/PositionUtils.cc
+++ b/gromos++/src/fit/PositionUtils.cc
@@ -77,7 +77,7 @@ Vec PositionUtils::com(const System &sys, utils::AtomSpecifier & atoms){
   double totalMass=0;
   Vec cm;
 
-  for(int a=0; a<atoms.size(); ++a){
+  for(unsigned int a=0; a<atoms.size(); ++a){
     cm += 
       atoms.mass(a) * atoms.pos(a) ;
     
@@ -93,7 +93,7 @@ Vec PositionUtils::com_v(const System &sys, utils::AtomSpecifier & atoms){
   double totalMass=0;
   Vec cm_v;
 
-  for(int a=0; a<atoms.size(); ++a){
+  for(unsigned int a=0; a<atoms.size(); ++a){
     cm_v += 
       atoms.mass(a) * atoms.vel(a) ;
     
@@ -122,7 +122,7 @@ Vec PositionUtils::cog(const System &sys, const Reference &ref){
 Vec PositionUtils::cog(const System &sys, utils::AtomSpecifier & atoms){
   Vec cg;
 
-  for(int a=0;a<atoms.size();++a)
+  for(unsigned int a=0;a<atoms.size();++a)
     cg += atoms.pos(a);
   
   cg /= atoms.size();
diff --git a/gromos++/src/fit/Reference.cc b/gromos++/src/fit/Reference.cc
index 352ead8a..b9a27e73 100644
--- a/gromos++/src/fit/Reference.cc
+++ b/gromos++/src/fit/Reference.cc
@@ -56,7 +56,7 @@ void Reference::addAtom(int m, int i){
 }
 
 void Reference::addAtomSpecifier(utils::AtomSpecifier as) {
-  for (int i=0; i < as.size(); ++i) addAtom(as.mol(i), as.atom(i));
+  for (unsigned int i=0; i < as.size(); ++i) addAtom(as.mol(i), as.atom(i));
 } 
 
 
diff --git a/gromos++/src/gcore/BbLink.cc b/gromos++/src/gcore/BbLink.cc
index d6c0d1b6..12032750 100644
--- a/gromos++/src/gcore/BbLink.cc
+++ b/gromos++/src/gcore/BbLink.cc
@@ -69,7 +69,7 @@ BbLink &BbLink::operator=(const BbLink &mt){
   return *this;
 }
 
-void BbLink::setLinkRes(const int a, const int i){
+void BbLink::setLinkRes(const unsigned int a, const unsigned int i){
   if(d_linkres.size() <= a)
    d_linkres.resize(a+1); 
   d_linkres[a]=i;
diff --git a/gromos++/src/gcore/BbLink.h b/gromos++/src/gcore/BbLink.h
index 66330104..6d920e19 100644
--- a/gromos++/src/gcore/BbLink.h
+++ b/gromos++/src/gcore/BbLink.h
@@ -48,7 +48,7 @@ namespace gcore{
      * @param a atom
      * @param i residue
      */
-    void setLinkRes(const int a, const int i);
+    void setLinkRes(const unsigned int a, const unsigned int i);
 
     // Accessors
     /** 
diff --git a/gromos++/src/gcore/LinearTopology.cc b/gromos++/src/gcore/LinearTopology.cc
index 884b799d..9ffec0ae 100644
--- a/gromos++/src/gcore/LinearTopology.cc
+++ b/gromos++/src/gcore/LinearTopology.cc
@@ -189,7 +189,16 @@ void LinearTopology::parse(gcore::System &sys)
 
     // add Atoms
     for(; int(atomCounter) < lastAtom; atomCounter++){
-      // std::cerr << "adding atom " << atomCounter << std::endl;
+
+      //std::cerr << "adding atom " << atomCounter << std::endl;
+      // convert poloffsite atom numbers to be relative to the molecule
+      if (d_atom[atomCounter].isPolarisable()) {
+        int poloffsite_i=d_atom[atomCounter].poloffsiteI();
+        int poloffsite_j=d_atom[atomCounter].poloffsiteJ();
+        d_atom[atomCounter].setPoloffsiteI(poloffsite_i-prevMol);
+        d_atom[atomCounter].setPoloffsiteJ(poloffsite_j-prevMol);
+      }
+
       mt->addAtom(d_atom[atomCounter]);
       
       // adapt exclusions:
@@ -537,9 +546,8 @@ void LinearTopology::moveAtoms(std::vector<std::pair<int, int> > moveatoms) {
   std::vector<std::pair<int, int> >::iterator p_it = moveatoms.begin(),
                                               p_to = moveatoms.end();
   for (;p_it != p_to;p_it++) {
-    int move_from=p_it->first-1;
-    int move_to=p_it->second-1;
-    int lower=std::min(move_from,move_to);
+    unsigned int move_from=p_it->first-1;
+    unsigned int move_to=p_it->second-1;
     int upper=std::max(move_from,move_to);
     
     if (d_resmap[move_from] != d_resmap[move_to])  {
@@ -557,13 +565,13 @@ void LinearTopology::moveAtoms(std::vector<std::pair<int, int> > moveatoms) {
     // insert/remove atoms and fill map that connects old and new indices
     change_map[move_from]=move_to;
     if (move_from > move_to) {
-        for (int ii=move_to; ii < move_from; ii++) {
+        for (unsigned int ii=move_to; ii < move_from; ii++) {
           change_map[ii]=ii+1;
         }
         d_atom.erase(firstatom+move_from);
         d_atom.insert(firstatom+move_to, at);
     } else if (move_from < move_to) {
-        for (int ii=move_from+1; ii <= move_to; ii++) {
+        for (unsigned int ii=move_from+1; ii <= move_to; ii++) {
           change_map[ii]=ii-1;
         }
         d_atom.insert(firstatom+move_to, at);
diff --git a/gromos++/src/gcore/MoleculeTopology.cc b/gromos++/src/gcore/MoleculeTopology.cc
index 724ec7ce..269a0549 100644
--- a/gromos++/src/gcore/MoleculeTopology.cc
+++ b/gromos++/src/gcore/MoleculeTopology.cc
@@ -18,6 +18,7 @@
 
 using namespace std;
 using gcore::MoleculeTopology;
+using gcore::AtomIterator;
 using gcore::BondIterator;
 using gcore::BondDipoleIterator;
 using gcore::AngleIterator;
@@ -36,6 +37,7 @@ using gcore::AtomTopology;
 class gcore::MoleculeTopology_i{
 
   friend class gcore::MoleculeTopology;
+  friend class gcore::AtomIterator;
   friend class gcore::BondIterator;
   friend class gcore::BondDipoleIterator;
   friend class gcore::AngleIterator;
@@ -217,6 +219,56 @@ const string &MoleculeTopology::resName(int i)const{
   return d_this->d_resNames[i];
 }
 
+class gcore::AtomIterator_i{
+  friend class gcore::AtomIterator;
+  vector<AtomTopology>::iterator d_it;
+  const MoleculeTopology *d_mt;
+  // not implemented
+  AtomIterator_i(const AtomIterator_i&);
+  AtomIterator_i &operator=(const AtomIterator_i &);
+public:
+  AtomIterator_i():
+    d_it(){d_mt=0;}
+};
+
+gcore::AtomIterator::AtomIterator(const MoleculeTopology &mt):
+  d_this(new AtomIterator_i())
+{
+  d_this->d_it=mt.d_this->d_atoms.begin();
+  d_this->d_mt=&mt;
+}
+
+AtomIterator::~AtomIterator(){delete d_this;}
+
+void AtomIterator::operator++(){
+  ++(d_this->d_it);
+}
+
+const AtomTopology &AtomIterator::operator()()const{
+  return *(d_this->d_it);
+}
+
+AtomTopology &AtomIterator::operator()(){
+  return const_cast<AtomTopology&>(*(d_this->d_it));
+}
+
+AtomIterator::operator bool()const{
+  return d_this->d_it != d_this->d_mt->d_this->d_atoms.end();
+}
+
+bool AtomIterator::last() const{
+  // these iterator are not random access iterators and so they don't support
+  // end()-1. Thus we have to copy the iterator and advance it an check whether
+  // it's at the end.
+  vector<AtomTopology>::iterator it(d_this->d_it);
+  ++it;
+  return it == d_this->d_mt->d_this->d_atoms.end();
+}
+
+bool AtomIterator::first() const{
+  return d_this->d_it == d_this->d_mt->d_this->d_atoms.begin();
+}
+
 class gcore::BondIterator_i{
   friend class gcore::BondIterator;
   set<Bond>::iterator d_it;
diff --git a/gromos++/src/gcore/MoleculeTopology.h b/gromos++/src/gcore/MoleculeTopology.h
index ae9c331c..f9912a79 100644
--- a/gromos++/src/gcore/MoleculeTopology.h
+++ b/gromos++/src/gcore/MoleculeTopology.h
@@ -46,6 +46,19 @@ namespace gcore{
   class MoleculeTopology{
     MoleculeTopology_i *d_this;
     // This class contains all topological information
+    /**
+     * Atom Iterator for a MoleculeTopology
+     *
+     * The MoleculeTopology Atom iterator is used to loop over the AtomTopologies 
+     * in a MoleculeTopology. 
+     * It is constructed with the MoleculeTopology as an argument. Use the 
+     * ++ operator to move to the next atom topology. The () operator returns the 
+     * current atom topology. 
+     * This can also be used as a boolean: the bool() returns 1 as long as 
+     * the iterator is not at the end of the atom list.
+     * @author R. Buergi
+     */
+    friend class AtomIterator;
     /**
      * Bond Iterator for the a MoleculeTopology
      *
@@ -298,7 +311,7 @@ namespace gcore{
 
   }; /* class MoleculeTopology */
 
-
+  class AtomIterator_i;
   class BondIterator_i;
   class BondDipoleIterator_i;
   class AngleIterator_i;
@@ -307,6 +320,24 @@ namespace gcore{
   class CrossDihedralIterator_i;
   class LJExceptionIterator_i;
 
+
+  class AtomIterator{
+    AtomIterator_i *d_this;
+    // not implemented
+    AtomIterator();
+    AtomIterator(const AtomIterator&);
+    AtomIterator &operator=(const AtomIterator &);
+  public:
+    AtomIterator(const MoleculeTopology &mt);
+    ~AtomIterator();
+    void operator++();
+    const AtomTopology &operator()()const;
+    AtomTopology & operator()();
+    operator bool()const;
+    bool last()const;
+    bool first()const;
+  };
+
   class BondIterator{
     BondIterator_i *d_this;
     // not implemented
diff --git a/gromos++/src/gcore/System.h b/gromos++/src/gcore/System.h
index 16e1be49..51f38ab3 100644
--- a/gromos++/src/gcore/System.h
+++ b/gromos++/src/gcore/System.h
@@ -210,12 +210,12 @@ namespace gcore{
   };
 
   inline const Molecule &System::mol(int i)const{
-    assert (i < this->numMolecules());
+    assert (i < int(this->numMolecules()));
     return *d_mol[i];
   }
 
   inline Molecule &System::mol(int i){
-    assert (i < this->numMolecules());
+    assert (i < int(this->numMolecules()));
     return *d_mol[i];
   }
 
diff --git a/gromos++/src/gio/InBuildingBlock.cc b/gromos++/src/gio/InBuildingBlock.cc
index 2d77e91a..16482555 100644
--- a/gromos++/src/gio/InBuildingBlock.cc
+++ b/gromos++/src/gio/InBuildingBlock.cc
@@ -820,7 +820,6 @@ void gio::InBuildingBlock_i::readDihedrals(BbSolute &bb, std::string bname, std:
 void gio::InBuildingBlock_i::readLJExceptions(BbSolute &bb, std::string bname, std::string resname) {
   int i[4], num;
   std::string s;
-  double d[2];
 
   _lineStream >> num;
   if (_lineStream.fail())
diff --git a/gromos++/src/gio/InLinkingBlock.cc b/gromos++/src/gio/InLinkingBlock.cc
index 51da3690..7bbe2dca 100644
--- a/gromos++/src/gio/InLinkingBlock.cc
+++ b/gromos++/src/gio/InLinkingBlock.cc
@@ -162,7 +162,7 @@ void gio::InLinkingBlock_i::readLink(std::vector<std::string> &buffer) {
 
 void gio::InLinkingBlock_i::readAtoms(BbLink &bb, std::string bname, std::string resname) {
   double d[3];
-  int i[5], na, npe;
+  int i[5], na;
   std::string s;
 
   _lineStream >> na;
diff --git a/gromos++/src/gio/InParameter.cc b/gromos++/src/gio/InParameter.cc
index 7a65250d..1795766c 100644
--- a/gromos++/src/gio/InParameter.cc
+++ b/gromos++/src/gio/InParameter.cc
@@ -86,7 +86,7 @@ int gio::InParameter_i::_initBlock(std::vector<std::string> &buffer,
     throw InParameter::Exception("Parameter file "+name()+
 				" is corrupted. No "+blockname+
 				" block!");
-    if(buffer[buffer.size()-1].find("END")!=0)
+  if(buffer[buffer.size()-1].find("END")!=0)
       throw InParameter::Exception("Topology file " + name() +
 				       " is corrupted. No END in "+blockname+
 				       " block. Got\n"
diff --git a/gromos++/src/gio/InTopology.cc b/gromos++/src/gio/InTopology.cc
index 5a9bd9a3..2a0a0e05 100644
--- a/gromos++/src/gio/InTopology.cc
+++ b/gromos++/src/gio/InTopology.cc
@@ -375,7 +375,7 @@ void gio::InTopology_i::parseForceField() {
           throw InTopology::Exception("Different number of bond types in BONDTYPE and HARMBONDTYPE block.");
     }
     if (quart_k.size() && harm_k.size()){
-        for (int i=0; i < quart_k.size(); i++) {  
+        for (unsigned int i=0; i < quart_k.size(); i++) {  
           if ( quart_d0[i] !=  harm_d0[i] ) {
             ostringstream msg;
             msg << "Bond type "<< i <<" in BONDTYPE (b0="<<quart_d0[i]<<") and HARMBONDTYPE (b0="<<harm_d0[i]<<") block do not agree.";
@@ -385,11 +385,11 @@ void gio::InTopology_i::parseForceField() {
           d_gff.addBondType(BondType(i, quart_k[i], harm_k[i], quart_d0[i]));
         }
     } else if (quart_k.size()) {  
-          for (int i=0; i < quart_k.size(); i++) { 
+          for (unsigned int i=0; i < quart_k.size(); i++) { 
             d_gff.addBondType(BondType(i, quart_k[i], quart_d0[i], true));
           }
     } else if (harm_k.size()) {
-          for (int i=0; i < harm_k.size(); i++) { 
+          for (unsigned int i=0; i < harm_k.size(); i++) { 
             d_gff.addBondType(BondType(i, harm_k[i], harm_d0[i], false));
           }
     }
@@ -492,7 +492,7 @@ void gio::InTopology_i::parseForceField() {
           throw InTopology::Exception("Different number of angle types in BONDANGLETYPE and HARMBONDANGLETYPE block.");
       }
       if (cosharm_k.size() && harm_k.size()){
-        for (int i=0; i < cosharm_k.size(); i++) {  
+        for (unsigned int i=0; i < cosharm_k.size(); i++) {  
           if ( cosharm_th[i] !=  harm_th[i] ) {
             ostringstream msg;
             msg << "Bond angle "<< i <<" in BONDANGLETYPE (theta="<<cosharm_th[i]<<") and HARMBONDANGLETYPE (theta="<<harm_th[i]<<") block do not agree.";
@@ -502,7 +502,7 @@ void gio::InTopology_i::parseForceField() {
           d_gff.addAngleType(AngleType(i, cosharm_k[i], harm_k[i], cosharm_th[i]));
         }
       } else if (cosharm_k.size()) {  
-          for (int i=0; i < cosharm_k.size(); i++) {          
+          for (unsigned int i=0; i < cosharm_k.size(); i++) {          
             try {
                 d_gff.addAngleType(AngleType(i, cosharm_k[i], cosharm_th[i]));
             } catch (gromos::Exception & exp) {
@@ -515,7 +515,7 @@ void gio::InTopology_i::parseForceField() {
           }
       } else if (harm_k.size()) {
             std::cerr << "# Warning: Angle types: only harmonic force constants (HARMBONDANGLETYPE block) found, setting cosine harmonic force constant to -1.0." << std::endl;
-          for (int i=0; i < harm_k.size(); i++) { 
+          for (unsigned int i=0; i < harm_k.size(); i++) { 
             d_gff.addAngleType(AngleType(i, -1, harm_k[i], harm_th[i]));
           }
       }
diff --git a/gromos++/src/gio/OutG96.cc b/gromos++/src/gio/OutG96.cc
index dfe30ea3..a92853a9 100644
--- a/gromos++/src/gio/OutG96.cc
+++ b/gromos++/src/gio/OutG96.cc
@@ -324,7 +324,7 @@ void gio::OutG96_i::writeAtomSpecifier(const AtomSpecifier& atoms) {
   d_os.setf(ios::fixed, ios::floatfield);
   d_os.precision(9);
   d_os << "# selected " << atoms.size() << " atoms" << endl;
-  for (int i = 0; i < atoms.size(); ++i) {
+  for (unsigned int i = 0; i < atoms.size(); ++i) {
     d_os << setw(15) << atoms.pos(i)[0]
             << setw(15) << atoms.pos(i)[1]
             << setw(15) << atoms.pos(i)[2]
diff --git a/gromos++/src/gio/OutG96S.cc b/gromos++/src/gio/OutG96S.cc
index a4c2d321..5812ae1a 100644
--- a/gromos++/src/gio/OutG96S.cc
+++ b/gromos++/src/gio/OutG96S.cc
@@ -494,7 +494,7 @@ void gio::OutG96S_i::writeAtomSpecifier(const AtomSpecifier & atoms, bool vel) {
   d_os << "# selected " << atoms.size() << " atoms" << endl;
   d_os.setf(ios::unitbuf);
 
-  for (int i = 0; i < atoms.size(); ++i) {
+  for (unsigned int i = 0; i < atoms.size(); ++i) {
     d_os.setf(ios::right, ios::adjustfield);
     int offset = 1;
     if (atoms.mol(i) >= 0) {
diff --git a/gromos++/src/gio/OutPdb.cc b/gromos++/src/gio/OutPdb.cc
index 86498ba9..617ec00c 100644
--- a/gromos++/src/gio/OutPdb.cc
+++ b/gromos++/src/gio/OutPdb.cc
@@ -266,7 +266,7 @@ void gio::OutPdb_i::writeAtomSpecifier(const AtomSpecifier& atoms) {
   int resoff = 0;
   char chain = 'A';
 
-  for (int i = 0; i < atoms.size(); ++i) {
+  for (unsigned int i = 0; i < atoms.size(); ++i) {
     int maxmol = atoms.mol(i);
     if (maxmol < 0) maxmol = sys.numMolecules();
     count = atoms.atom(i);
diff --git a/gromos++/src/gio/Outvmdam.cc b/gromos++/src/gio/Outvmdam.cc
index 29cb631c..98c4c49e 100644
--- a/gromos++/src/gio/Outvmdam.cc
+++ b/gromos++/src/gio/Outvmdam.cc
@@ -169,7 +169,7 @@ void gio::Outvmdam_i::writeAtomSpecifier(const AtomSpecifier& atoms) {
   d_os.setf(ios::unitbuf);
   d_os.precision(3);
 
-  for (int i = 0; i < atoms.size(); ++i) {
+  for (unsigned int i = 0; i < atoms.size(); ++i) {
     ++d_count;
 
     d_os.setf(ios::right, ios::adjustfield);
diff --git a/gromos++/src/gmath/Correlation.cc b/gromos++/src/gmath/Correlation.cc
index 87ad36e2..4a95754d 100644
--- a/gromos++/src/gmath/Correlation.cc
+++ b/gromos++/src/gmath/Correlation.cc
@@ -51,22 +51,22 @@ namespace gmath
   void Correlation::calc_direct()
   {
     int num=d_f.size();
-    for(int i=0; i< num; i++) d_f[i]=0.0;
+    for(unsigned int i=0; i< num; i++) d_f[i]=0.0;
     if(d_vec){
-      for(int i=0; i<num; i++){
-        for(int j=i+1; j<num; j++){
+      for(unsigned int i=0; i<num; i++){
+        for(unsigned int j=i; j<num; j++){
           d_f[j-i]+=(*d_va)[i].dot((*d_vb)[j]);
         }
       }
     }
     else{
-      for(int i=0; i<num; i++){
-        for(int j=i+1; j<num; j++){
+      for(unsigned int i=0; i<num; i++){
+        for(unsigned int j=i; j<num; j++){
           d_f[j-i]+=(*d_a)[i] * (*d_b)[j];
         }
       }
     }
-    for(int i=1; i<num; i++){
+    for(unsigned int i=0; i<num; i++){
       d_f[i]/=(num-i);
     }
     d_calc=true;
@@ -172,10 +172,10 @@ namespace gmath
       
       // do a (direct) calculation
       int num=d_f.size();
-      for(int i=0; i<num; i++) d_f[i]=0.0;
-      for(int i=0; i<num; i++){
+      for(unsigned int i=0; i<num; i++) d_f[i]=0.0;
+      for(unsigned int i=0; i<num; i++){
 	v[0]=(*d_a)[i];
-        for(int j=i+1; j<num; j++){
+        for(unsigned int j=i; j<num; j++){
 	  v[1]=(*d_b)[j];
 	  
 	  e.setValues(v);
@@ -183,7 +183,7 @@ namespace gmath
         }
       }
 
-      for(int i=1; i<num; i++){
+      for(unsigned int i=0; i<num; i++){
 	d_f[i]/=(num-i);
       }
       d_calc=true;
@@ -193,7 +193,7 @@ namespace gmath
     assert(i < int(d_f.size()));
     return d_f[i];
   }
-  int Correlation::size(){
+  unsigned int Correlation::size(){
     return d_f.size();
   }
 
diff --git a/gromos++/src/gmath/Correlation.h b/gromos++/src/gmath/Correlation.h
index 71020848..3820a844 100644
--- a/gromos++/src/gmath/Correlation.h
+++ b/gromos++/src/gmath/Correlation.h
@@ -81,7 +81,7 @@ namespace gmath
       /**
        * Accessor to the size of the correlation function
        */
-      int size();
+      unsigned int size();
       /**
        * Method to calculate the power spectrum of the correlation function
        * using fft
diff --git a/gromos++/src/gmath/Matrix.cc b/gromos++/src/gmath/Matrix.cc
index 633c98ba..76911826 100644
--- a/gromos++/src/gmath/Matrix.cc
+++ b/gromos++/src/gmath/Matrix.cc
@@ -15,9 +15,9 @@ namespace gmath{
   Matrix::Matrix(int rows, int columns, double value){
     d_rows=rows;
     d_columns=columns;
-    d_val=new (double* [d_rows]);
+    d_val=new double* [d_rows];
     for (int i=0;i<d_rows;++i){
-      d_val[i]=new (double [d_columns]);
+      d_val[i]=new double [d_columns];
       for(int j=0;j<d_columns;++j){
 	d_val[i][j]=value;
       }
@@ -27,9 +27,9 @@ namespace gmath{
   Matrix::Matrix(const Matrix &mat){
     d_rows=mat.d_rows;
     d_columns = mat.d_columns;
-    d_val=new (double* [d_rows]);
+    d_val=new double* [d_rows];
     for (int i=0;i<d_rows;++i){
-      d_val[i]=new (double [d_columns]);
+      d_val[i]=new double [d_columns];
       for(int j=0;j<d_columns;++j){
 	d_val[i][j]=mat.d_val[i][j];
       }
@@ -39,9 +39,9 @@ namespace gmath{
   Matrix::Matrix(const Vec &v, const Vec &w){
     d_rows=3;
     d_columns=3;
-    d_val=new (double* [d_rows]);
+    d_val=new double* [d_rows];
     for (int i=0;i<d_rows;++i){
-      d_val[i]=new (double [d_columns]);
+      d_val[i]=new double [d_columns];
       for(int j=0;j<d_columns;++j){
 	d_val[i][j]=v[i]*w[j];
       }
@@ -51,9 +51,9 @@ namespace gmath{
   Matrix::Matrix(const Vec &u, const Vec &v, const Vec &w){
     d_rows=3;
     d_columns=3;
-    d_val=new (double* [d_rows]);
+    d_val=new double* [d_rows];
     for (int i=0;i<d_rows;++i){
-      d_val[i]=new (double [d_columns]);
+      d_val[i]=new double [d_columns];
       d_val[i][0]=u[i];
       d_val[i][1]=v[i];
       d_val[i][2]=w[i];
diff --git a/gromos++/src/pb/Ewald_edir.cc b/gromos++/src/pb/Ewald_edir.cc
index 14716047..549b66f6 100644
--- a/gromos++/src/pb/Ewald_edir.cc
+++ b/gromos++/src/pb/Ewald_edir.cc
@@ -68,7 +68,7 @@ using namespace std;
      eir.resize(kmax);
      for (int i=0; i<kmax; i++){
          eir[i].resize(atoms.size());
-            for (int j=0; j<atoms.size(); j++){
+            for (unsigned int j=0; j<atoms.size(); j++){
                  eir[i][j].resize(3);
              }}
      
@@ -227,8 +227,8 @@ using namespace std;
 
 		double energy = 0.0;
 
-		for (int i=0; i < atoms.size()-1; ++i) {
-			for (int j=i+1; j < atoms.size(); ++j) {
+		for (unsigned int i=0; i < atoms.size()-1; ++i) {
+			for (unsigned int j=i+1; j < atoms.size(); ++j) {
                       
                             
 
@@ -238,20 +238,20 @@ using namespace std;
                           double distance=dd.abs();
 
 
-                          bool excluded_atoms;
+                          //bool excluded_atoms;
 
 
-                          int im=atoms.mol(i);
-                          int ia=atoms.atom(i);
-                          int jm=atoms.mol(j);
-                          int ja=atoms.atom(j);
+                          //int im=atoms.mol(i);
+                          //int ia=atoms.atom(i);
+                          //int jm=atoms.mol(j);
+                          //int ja=atoms.atom(j);
 
                       
 
 
-                          gcore::System *sss= atoms.sys();
+                          //gcore::System *sss= atoms.sys();
                           
-                         excluded_atoms = false;
+                         //excluded_atoms = false;
 
                         //  if (im==jm && sss->mol(im).topology().atom(ia).exclusion().contains(ja)){
                         //      excluded_atoms=true;
@@ -316,7 +316,8 @@ using namespace std;
 	void Ewald_edir::tabulate_eir(double (& lll)[3]) {
 
             
-	int  i,j,m;
+	unsigned int  i;
+        int j,m;
 
               try{
 		if (kmax < 1) {
@@ -373,7 +374,7 @@ using namespace std;
 
 	
 			for (int i=0; i < kmax; ++i) {
-				for (int ii=0; ii < atoms.size(); ++ii) {
+				for (unsigned int ii=0; ii < atoms.size(); ++ii) {
 					for (int iii=0; iii < 3; ++iii) {
 						 eir[i][ii][iii]= complex<double> (0.0, 0.0);
                                               					}
@@ -394,7 +395,8 @@ using namespace std;
 		tabulate_eir(lll);
 
 
-		int  ix,iy,iz,n;
+		int  ix,iy,iz;
+                unsigned int n;
 		double tmp,cs,ss,ak,akv,mx,my,mz,m2;
 		double energy=0.0;
 		int lowiy=0;
@@ -481,7 +483,7 @@ using namespace std;
 
 
 			for (int i=0; i < kmax; ++i) {
-				for (int ii=0; ii < atoms.size(); ++ii) {
+				for (unsigned int ii=0; ii < atoms.size(); ++ii) {
 					for (int iii=0; iii < 3; ++iii) {
 						 eir[i][ii][iii]= complex<double> (0.0, 0.0);
                                               					}
@@ -502,7 +504,8 @@ using namespace std;
 		tabulate_eir(lll);
 
 
-		int  ix,iy,iz,n;
+		int  ix,iy,iz;
+                unsigned int n;
 		double tmp,cs,ss,ak,akv,mx,my,mz,m2;
 		double energy=0.0;
 
@@ -647,8 +650,8 @@ using namespace std;
 		//tabulate_eir(lll);
 
 
-		int  ix,iy,iz,n;
-		double tmp,cs,ss,ak,akv,mx,my,mz,m2;
+		int  ix,iy,iz;
+		double tmp,ak,mx,my,mz,m2;
 	
 
 
@@ -686,7 +689,7 @@ using namespace std;
 
                 // now get stildesquare
                 double q2sum=0;
-               for (int i=0; i < atoms.size();++i){
+               for (unsigned int i=0; i < atoms.size();++i){
                     q2sum += atoms.charge(i) * atoms.charge(i);
                }
                 std::cout << "a2sum " << a2sum << "\n";
@@ -704,7 +707,7 @@ using namespace std;
 			double ewaldcoeff) {
 
 		double q2_sum = 0.0;
-		for (int i=0; i < atoms.size();++i){
+		for (unsigned int i=0; i < atoms.size();++i){
                     q2_sum += atoms.charge(i) * atoms.charge(i);
                 }
 		return ewaldcoeff * ppp.getFPEPSI() * q2_sum/(sqrt(ppp.getPI()));
@@ -715,7 +718,7 @@ using namespace std;
 	double Ewald_edir::minusA1timesSSquare(
 			double ewaldcoeff) {
 		double q_sum = 0.0;
-		for (int i=0; i < atoms.size();++i) {q_sum += atoms.charge(i);};
+		for (unsigned int i=0; i < atoms.size();++i) {q_sum += atoms.charge(i);};
 	
 		double vol = box[0]*box[1]*box[2];
 		return((q_sum*q_sum * ppp.getPI() * ppp.getFPEPSI())/(2.0*vol*ewaldcoeff*ewaldcoeff));
@@ -724,7 +727,7 @@ using namespace std;
         double Ewald_edir::A1timesStildeSquare(
 			double ewaldcoeff) {
 		double q_sum = 0.0;
-		for (int i=0; i < atoms.size();++i){ q_sum += (atoms.charge(i) * atoms.charge(i));};
+		for (unsigned int i=0; i < atoms.size();++i){ q_sum += (atoms.charge(i) * atoms.charge(i));};
 
 		double vol = box[0]*box[1]*box[2];
 		return( -1.0 * (q_sum * ppp.getPI() * ppp.getFPEPSI())/(2.0*vol*ewaldcoeff*ewaldcoeff));
@@ -740,7 +743,7 @@ using namespace std;
 
         double Ewald_edir::self_other(){
             double qsum = 0;
-            for (int i=0; i < atoms.size();++i) qsum += atoms.charge(i)*atoms.charge(i);
+            for (unsigned int i=0; i < atoms.size();++i) qsum += atoms.charge(i)*atoms.charge(i);
             double res= ppp.getFPEPSI()/2.0 * qsum / box[0] * (ppp.get_xiew()) ;
             return res;
         }
@@ -779,8 +782,8 @@ using namespace std;
 
 		double energy = 0.0;
 
-		for (int i=0; i < atoms.size()-1; ++i) {
-			for (int j=i+1; j < atoms.size(); ++j) {
+		for (unsigned int i=0; i < atoms.size()-1; ++i) {
+			for (unsigned int j=i+1; j < atoms.size(); ++j) {
 
 
 
@@ -790,20 +793,20 @@ using namespace std;
                           double distance=dd.abs();
 
 
-                          bool excluded_atoms;
+                          //bool excluded_atoms;
 
 
-                          int im=atoms.mol(i);
-                          int ia=atoms.atom(i);
-                          int jm=atoms.mol(j);
-                          int ja=atoms.atom(j);
+                          //int im=atoms.mol(i);
+                          //int ia=atoms.atom(i);
+                          //int jm=atoms.mol(j);
+                          //int ja=atoms.atom(j);
 
 
 
 
-                          gcore::System *sss= atoms.sys();
+                          //gcore::System *sss= atoms.sys();
 
-                          excluded_atoms = false;
+                          //excluded_atoms = false;
 
                         //  if (im==jm && sss->mol(im).topology().atom(ia).exclusion().contains(ja)){
                         //      excluded_atoms=true;
@@ -833,29 +836,29 @@ using namespace std;
 
 		double energy = 0.0;
 
-		for (int i=0; i < atoms.size()-1; ++i) {
-			for (int j=i+1; j < atoms.size(); ++j) {
+		for (unsigned int i=0; i < atoms.size()-1; ++i) {
+			for (unsigned int j=i+1; j < atoms.size(); ++j) {
 
 
 
                         // get the distance between the atoms
                           gmath::Vec jimage = (*pbc).nearestImage(atoms.pos(i),atoms.pos(j), thebox);
                           gmath::Vec dd = atoms.pos(i)-jimage;
-                          double distance=dd.abs();
+                          //double distance=dd.abs();
 
 
-                          bool excluded_atoms;
+                          //bool excluded_atoms;
 
 
-                          int im=atoms.mol(i);
-                          int ia=atoms.atom(i);
-                          int jm=atoms.mol(j);
-                          int ja=atoms.atom(j);
+                          //int im=atoms.mol(i);
+                          //int ia=atoms.atom(i);
+                          //int jm=atoms.mol(j);
+                          //int ja=atoms.atom(j);
 
 
 
 
-                          gcore::System *sss= atoms.sys();
+                          //gcore::System *sss= atoms.sys();
 
                        //   excluded_atoms = false;
 
diff --git a/gromos++/src/pb/FDPoissonBoltzmann.cc b/gromos++/src/pb/FDPoissonBoltzmann.cc
index d512a6a1..99c9ac0b 100755
--- a/gromos++/src/pb/FDPoissonBoltzmann.cc
+++ b/gromos++/src/pb/FDPoissonBoltzmann.cc
@@ -52,7 +52,7 @@ using pb::PB_Parameters;
     std::cout << "# GPY " << GPY << endl;
     std::cout << "# GPZ " << GPZ << endl;
     std::cout << "# epssolvent " << epssolvent << endl;
-    for (int i=0; i<atoms_to_charge.size(); i++){
+    for (unsigned int i=0; i<atoms_to_charge.size(); i++){
         std::cout << "# atom " << i << " charge " << atoms_to_charge.charge(i) << endl;
         std::cout << "# atom " << i << " radius " << atoms_to_charge.radius(i) << endl;
     }
@@ -635,10 +635,10 @@ std::cout << "# iter anorm acceptance*znorm " << iter << " " << anorm << " " <<
         double FDPoissonBoltzmann::dGelec(){
 
 
-                double potential = 0;
+                //double potential = 0;
                 double potential_rest = 0;
 
-                int I, J, K;
+                //int I, J, K;
 
 
                 double phitot = 0.0;
@@ -959,7 +959,7 @@ int size=atoms.size();
         double range1 = (atoms.radius(i))/gridspacing;
 
          double ijk[3];
-         double fraction[3];
+         //double fraction[3];
 
          double fuzoff = 0.5;
          int I, J, K;
@@ -967,7 +967,7 @@ int size=atoms.size();
          for (int j=0; j < 3; ++j) {
           ijk[j] = (  (atoms.pos(i))[j] - gridstart[j])/gridspacing;
 	      //get the fractional remainder
-          fraction[j] = ijk[j] - (int) ijk[j];
+          //fraction[j] = ijk[j] - (int) ijk[j];
          }
 
          int low1 = (int) (ijk[2] - range1 + 1.0);
@@ -1051,7 +1051,7 @@ int size=atoms.size();
           double range1 = (atoms.radius(i))/gridspacing;
 
          double ijk[3];
-         double fraction[3];
+         //double fraction[3];
 
          double fuzoff = 0.5;
 
@@ -1061,7 +1061,7 @@ int size=atoms.size();
          for (int j=0; j < 3; ++j) {
              ijk[j] = ((atoms.pos(i))[j] - gridstart[j])/gridspacing;
 	    //get the fractional remainder
-             fraction[j] = ijk[j] - (int) ijk[j];
+            // fraction[j] = ijk[j] - (int) ijk[j];
 	    }
 
          int low1 = (int) (ijk[2] - range1 + 1.0);
@@ -1135,7 +1135,7 @@ int size=atoms.size();
      double range1 = atoms.radius(i) /gridspacing;
 
          double ijk[3];
-         double fraction[3];
+         //double fraction[3];
 
          double fuzoff = 0.5;
 
@@ -1144,7 +1144,7 @@ int size=atoms.size();
          for (int j=0; j < 3; ++j) {
              ijk[j] = (   (atoms.pos(i))[j] - gridstart[j])/gridspacing;
 	    //get the fractional remainder
-             fraction[j] = ijk[j] - (int) ijk[j];
+            // fraction[j] = ijk[j] - (int) ijk[j];
 	    }
 
          int low1 = (int) (ijk[1] - range1 + 1.0);
diff --git a/gromos++/src/pb/FFTChargeDipole_RF.cc b/gromos++/src/pb/FFTChargeDipole_RF.cc
index 1fe531ee..b17ecfef 100755
--- a/gromos++/src/pb/FFTChargeDipole_RF.cc
+++ b/gromos++/src/pb/FFTChargeDipole_RF.cc
@@ -71,7 +71,7 @@ FFTChargeDipole_RF::FFTChargeDipole_RF(
 	/* kr the product of the norm of the k-vector and the cutoff distance */
 
 	double FFTChargeDipole_RF::fkr(double kr) {
-            double res;
+            //double res;
             //return FFTDipoleDipole_RF.fkr(kr);
 
 
diff --git a/gromos++/src/pb/FFTChargeShapingFunction.h b/gromos++/src/pb/FFTChargeShapingFunction.h
index 5bb57304..248166bb 100755
--- a/gromos++/src/pb/FFTChargeShapingFunction.h
+++ b/gromos++/src/pb/FFTChargeShapingFunction.h
@@ -28,7 +28,7 @@ class FFTChargeShapingFunction{
     //methods
 	/* FFTChargeShapingFunction getCSF(int type);*/
 	virtual double calc(double xx, double yy, double zz,
-		            double alpha, double eps0){}
+		            double alpha, double eps0) {}
 
 
 
diff --git a/gromos++/src/pb/FFTInsideOutside.cc b/gromos++/src/pb/FFTInsideOutside.cc
index c47cb26e..c9a65dfb 100755
--- a/gromos++/src/pb/FFTInsideOutside.cc
+++ b/gromos++/src/pb/FFTInsideOutside.cc
@@ -35,7 +35,7 @@ void FFTInsideOutside::inside_sol(
 			 std::vector<double> & in) {
 		
 		
-        	int numGridPoints = gridN[0] * gridN[1] * gridN[2];
+        	//int numGridPoints = gridN[0] * gridN[1] * gridN[2];
 	
 	//	int fullyInside = 0;
 		
@@ -60,8 +60,8 @@ void FFTInsideOutside::inside_sol(
 		int numAtoms = atoms.size();
 	
 		
-		int cell[3];
-		int nim[3];
+		//int cell[3];
+		//int nim[3];
 
 
 
@@ -338,7 +338,7 @@ void FFTInsideOutside::inside_sol(
 
 
 
-                int index;
+                unsigned int index;
                 double integral=0.0;
                 int inPoints_noB=inside.size(); /* without boundary */
 
@@ -382,4 +382,4 @@ void FFTInsideOutside::inside_sol(
          //       }
 
 
-        }
\ No newline at end of file
+        }
diff --git a/gromos++/src/pb/FFTPoisson.cc b/gromos++/src/pb/FFTPoisson.cc
index 35439698..34188e23 100755
--- a/gromos++/src/pb/FFTPoisson.cc
+++ b/gromos++/src/pb/FFTPoisson.cc
@@ -58,7 +58,7 @@ FFTPoisson::FFTPoisson(utils::AtomSpecifier atoms,utils::AtomSpecifier atoms_to_
 
 
         std::cout << "# epssolvent " << epssolvent << endl;
-        for (int i=0; i<atoms_to_charge.size();i++){
+        for (unsigned int i=0; i<atoms_to_charge.size();i++){
         std::cout << "# atom " << i << " charge " << atoms_to_charge.charge(i) << endl;
         std::cout << "# atom " << i << " radius " << atoms_to_charge.radius(i) << endl;
          }
@@ -388,7 +388,7 @@ FFTPoisson::FFTPoisson(utils::AtomSpecifier atoms,utils::AtomSpecifier atoms_to_
     // now shift
 
 
-    for (int i=0;i<atoms.size(); i++){
+    for (unsigned int i=0;i<atoms.size(); i++){
         atoms.pos(i)[0] += tvec[0];
         atoms.pos(i)[1] += tvec[1];
         atoms.pos(i)[2] += tvec[2];
diff --git a/gromos++/src/pb/FFTPoissonIterator.cc b/gromos++/src/pb/FFTPoissonIterator.cc
index 1c310ad7..f01bca8b 100755
--- a/gromos++/src/pb/FFTPoissonIterator.cc
+++ b/gromos++/src/pb/FFTPoissonIterator.cc
@@ -333,7 +333,7 @@ std::cout << "# deltaSigma = " << deltaSigma << endl;
 	
 	void FFTPoissonIterator::enforceZeroImaginaryComponent(std::vector<double> & complexVector) {
 		
-		for (int ii = 1; ii < complexVector.size(); ii += 2)
+		for (unsigned int ii = 1; ii < complexVector.size(); ii += 2)
 			complexVector[ii] = 0.0;
 	}
 	
@@ -606,7 +606,7 @@ std::cout << "# deltaSigma = " << deltaSigma << endl;
 		
 		double dg;
 		int i,j,k;
-		int ion;
+		unsigned int ion;
 		double p;
 		double fx,fy,fz;
 		
@@ -648,7 +648,7 @@ std::cout << "# deltaSigma = " << deltaSigma << endl;
 			
 			//System.out.println(i + " " + j + " " + k);
 			
-		                                                                                                                                                                                                                        +     fx *     fy *     fz *pot[k+1+gt.ngrdz*(j+1+gt.ngrdy*(i+1)) + shift7];	
+		                                                                                                                                  // +     fx *     fy *     fz *pot[k+1+gt.ngrdz*(j+1+gt.ngrdy*(i+1)) + shift7];	
 			
 			// we multiply the shifts by two, because pot represents a complex vector
 			p = 
@@ -677,7 +677,7 @@ std::cout << "# deltaSigma = " << deltaSigma << endl;
 
 		double dg;
 		int i,j,k;
-		int ion;
+		unsigned int ion;
 		double p;
 		double fx,fy,fz;
 
@@ -718,7 +718,7 @@ std::cout << "# deltaSigma = " << deltaSigma << endl;
 
 			//System.out.println(i + " " + j + " " + k);
 
-		                                                                                                                                                                                                                        +     fx *     fy *     fz *pot[k+1+gt.ngrdz*(j+1+gt.ngrdy*(i+1)) + shift7];
+                                                                                                                                                  //+     fx *     fy *     fz *pot[k+1+gt.ngrdz*(j+1+gt.ngrdy*(i+1)) + shift7];
 
 			// we multiply the shifts by two, because pot represents a complex vector
 			p =
diff --git a/gromos++/src/pb/FFTVacuumField_LS.cc b/gromos++/src/pb/FFTVacuumField_LS.cc
index 6c38be81..23d876b5 100755
--- a/gromos++/src/pb/FFTVacuumField_LS.cc
+++ b/gromos++/src/pb/FFTVacuumField_LS.cc
@@ -96,7 +96,7 @@ using pb::FFTChargeShapingFunction;
         void FFTVacuumField_LS::complexFromDouble(
 			std::vector <double> & doubleArray, std::vector <double> & complexArray) {
 
-		for (int ii = 0; ii < doubleArray.size(); ii++) {
+		for (unsigned int ii = 0; ii < doubleArray.size(); ii++) {
 			complexArray[2 * ii] =      doubleArray[ii];
 			complexArray[2 * ii + 1] =  doubleArray[ii] ; // 0;  // Maria, why not 0??
 		}
@@ -111,7 +111,7 @@ using pb::FFTChargeShapingFunction;
 
             try{
 		ion_count1 = 0;ion_count2 = 0;
-		for (int ion=0;ion< atoms.size();ion++) {
+		for (unsigned int ion=0;ion< atoms.size();ion++) {
 			if ((atoms.radius(ion) > (bc.alpha1+tinynum))){
 
 
@@ -305,15 +305,15 @@ using pb::FFTChargeShapingFunction;
 				ion_list1, ion_count1);		
 		
                 complexFromDouble(Vx1_store, complexTmp);
-		for (int i=0;i<destX.size();i++){
+		for (unsigned int i=0;i<destX.size();i++){
                     destX[i]+= complexTmp[i]*multipliers[i];
                     }
                 complexFromDouble(Vy1_store, complexTmp);
-                for (int i=0;i<destY.size();i++){
+                for (unsigned int i=0;i<destY.size();i++){
                     destY[i]+=complexTmp[i]*multipliers[i];
                     }
                 complexFromDouble(Vz1_store, complexTmp);
-                for (int i=0;i<destZ.size();i++){
+                for (unsigned int i=0;i<destZ.size();i++){
                     destZ[i]+=complexTmp[i]*multipliers[i];
                     }
 		
@@ -321,15 +321,15 @@ using pb::FFTChargeShapingFunction;
 				ion_list2, ion_count2);		
 
                 complexFromDouble(Vx2_store, complexTmp);
-                for (int i=0;i<destX.size();i++){
+                for (unsigned int i=0;i<destX.size();i++){
                    destX[i]+= complexTmp[i]*multipliers[i];
                     }
                 complexFromDouble(Vy2_store, complexTmp);
-                for (int i=0;i<destY.size();i++){
+                for (unsigned int i=0;i<destY.size();i++){
                    destY[i]+=complexTmp[i]*multipliers[i];
                     }
                 complexFromDouble(Vz2_store, complexTmp);
-                for (int i=0;i<destZ.size();i++){
+                for (unsigned int i=0;i<destZ.size();i++){
                    destZ[i]+=complexTmp[i]*multipliers[i];
                     }
 
diff --git a/gromos++/src/pb/FFTVacuumField_RF.cc b/gromos++/src/pb/FFTVacuumField_RF.cc
index a5c5b7cf..08b355ed 100755
--- a/gromos++/src/pb/FFTVacuumField_RF.cc
+++ b/gromos++/src/pb/FFTVacuumField_RF.cc
@@ -117,11 +117,11 @@ FFTVacuumField_RF::FFTVacuumField_RF(utils::AtomSpecifier atoms, FFTGridType gt,
               //  eCC[0]=(int) (rcut / gridD[0]) + 1;
               //  eCC[1]=(int) (int) (rcut / gridD[1]) + 1;
               //  eCC[2]=(int) (rcut / gridD[2]) + 1;
-              std:vector <double> scaledCharges;
+              std::vector <double> scaledCharges;
                 scaledCharges.resize(atoms.size());
 	
 		
-		for (int ii = 0; ii < atoms.size(); ii++) {
+		for (unsigned int ii = 0; ii < atoms.size(); ii++) {
 			scaledCharges[ii] = fpepsi * atoms.charge(ii);
 	    }
 
diff --git a/gromos++/src/utils/AtomSpecifier.cc b/gromos++/src/utils/AtomSpecifier.cc
index 9b95d27d..6e1fd744 100644
--- a/gromos++/src/utils/AtomSpecifier.cc
+++ b/gromos++/src/utils/AtomSpecifier.cc
@@ -240,7 +240,7 @@ int utils::AtomSpecifier::addGromosAtom(int a)
 
 int utils::AtomSpecifier::addMolecule(int m)
 {
-  if(m >= d_sys->numMolecules())
+  if(m >= int(d_sys->numMolecules()))
     throw utils::AtomSpecifier::Exception(" molecule number out of range.\n");
 
   for(int i=0; i < d_sys->mol(m).numAtoms(); ++i)
@@ -252,7 +252,7 @@ int utils::AtomSpecifier::addMolecule(int m)
 
 int utils::AtomSpecifier::addAtomStrict(int m, int a)
 {
-  if(m >= d_sys->numMolecules())
+  if(m >= int(d_sys->numMolecules()))
     throw utils::AtomSpecifier::Exception(" molecule number out of range.\n");
   if(m >= 0){
     if(a >= d_sys->mol(m).topology().numAtoms())
@@ -617,7 +617,7 @@ bool utils::AtomSpecifier::empty()const
 	  d_solventType.size() > 0);
 }
 
-int utils::AtomSpecifier::size()const
+unsigned int utils::AtomSpecifier::size()const
 {
   if(_expand()) _expandSolvent();
   return d_specatom.size();
diff --git a/gromos++/src/utils/AtomSpecifier.h b/gromos++/src/utils/AtomSpecifier.h
index 0903e32e..32aa2505 100644
--- a/gromos++/src/utils/AtomSpecifier.h
+++ b/gromos++/src/utils/AtomSpecifier.h
@@ -667,7 +667,7 @@ namespace utils
     /**
      * Accessor, returns the number of atoms in the AtomSpecifier
      */
-    int size()const;
+    unsigned int size()const;
     /**
      * Accessor, returns if empty
      */
diff --git a/gromos++/src/utils/AtomSpecifier_parse.cc b/gromos++/src/utils/AtomSpecifier_parse.cc
index 3b283a11..30846f50 100644
--- a/gromos++/src/utils/AtomSpecifier_parse.cc
+++ b/gromos++/src/utils/AtomSpecifier_parse.cc
@@ -49,7 +49,7 @@ void AtomSpecifier::parse(std::string s, int x)
 
   if (d_not_atoms != NULL) {
     // remove atoms
-    for (int i = 0; i < d_not_atoms->size(); ++i) {
+    for (unsigned int i = 0; i < d_not_atoms->size(); ++i) {
       removeAtom(d_not_atoms->mol(i), d_not_atoms->atom(i));
     }
     // and solvent types
@@ -62,7 +62,7 @@ void AtomSpecifier::parse(std::string s, int x)
     }
     d_solventType = new_solv;
     // check for specific atoms violation this
-    for(int i = 0; i < size(); i++) {
+    for(unsigned int i = 0; i < size(); i++) {
       if (mol(i) >= 0) continue;
 
       if (std::find(d_not_atoms->d_solventType.begin(),
@@ -324,7 +324,7 @@ void AtomSpecifier::parse_va(std::string s, int x)
 
 void AtomSpecifier::parse_minus(std::string s, int x) {
   AtomSpecifier minus(*d_sys, s, x);
-  for(int i = 0; i < minus.size(); ++i) {
+  for(unsigned int i = 0; i < minus.size(); ++i) {
     if (minus.mol(i) < 0)
       throw Exception("Solvent is not supported in the minus() atom specifier. Use not() instead.");
     removeAtom(minus.mol(i), minus.atom(i));
diff --git a/gromos++/src/utils/CubeSystem.hcc b/gromos++/src/utils/CubeSystem.hcc
index ea2092c9..251e924b 100644
--- a/gromos++/src/utils/CubeSystem.hcc
+++ b/gromos++/src/utils/CubeSystem.hcc
@@ -708,7 +708,7 @@ void CubeSystem<Type>::assign_outer_neighbours(const int my_index){
 template <typename Type>
 void CubeSystem<Type>::assign_inner_neighbours(){
 
-    for(int i=0; i < cube_list.size(); ++i){
+    for(unsigned int i=0; i < cube_list.size(); ++i){
 
         const Vec my_pos=calc_position(i);
 
diff --git a/gromos++/src/utils/Disicl.cc b/gromos++/src/utils/Disicl.cc
index 4ccee3b3..72fdc9ee 100644
--- a/gromos++/src/utils/Disicl.cc
+++ b/gromos++/src/utils/Disicl.cc
@@ -34,7 +34,6 @@ vector<string> Dscl::split(const string &s, char delim) {
 void Dscl::readLibrary(gio::Ginstream &lib) {  
   map<std::string, std::vector<std::string> > libContent;
   std::vector<std::string> buffer;
-  typedef multimap<string,vector<double> >::value_type mapTypeVec;
   
   // check for valid block structure and read lib file contents
   while(!lib.stream().eof()) {
@@ -92,7 +91,7 @@ void Dscl::readLibAng(std::vector<std::string> &buffer) {
     numAng=libAngNames.size();
     
     // check for complex atom specifications
-    for (int i=0; i< angAtomsUniq.size(); i++) {
+    for (unsigned int i=0; i< angAtomsUniq.size(); i++) {
       string &s=angAtomsUniq[i];
       multAtom ma;
 
@@ -108,7 +107,7 @@ void Dscl::readLibAng(std::vector<std::string> &buffer) {
         ma.defaultAtom=substrings[0];
         
         // divide residue and atomnames
-        for (int j=1; j<substrings.size(); j++) {
+        for (unsigned int j=1; j<substrings.size(); j++) {
           vector<string> sstring=split(substrings[j],':');
           if (sstring.size()<2) {
             throw gromos::Exception("disicl","ERROR in library angle specifications!");
@@ -117,7 +116,7 @@ void Dscl::readLibAng(std::vector<std::string> &buffer) {
           // there could be more than one residue for a given atomname
           vector<string> ssstring=split(sstring[0],',');
           if (ssstring.size() >1) {
-            for (int k=0; k<ssstring.size(); k++) {
+            for (unsigned int k=0; k<ssstring.size(); k++) {
               ma.atoms.push_back(sstring[1]);
               ma.residues.push_back(ssstring[k]);
             }
@@ -148,17 +147,16 @@ void Dscl::readLibAng(std::vector<std::string> &buffer) {
 
     //some feedback to see if multi-atom input was read correctly
     cerr<< "# Atoms in angle definitions:\n#\t";
-    for (int i=0; i< angAtomsUniq.size(); i++) {
+    for (unsigned int i=0; i< angAtomsUniq.size(); i++) {
       cerr << angAtomsUniq[i] << "  ";
     }
     cerr << "\n";
     if (multAtoms.size()>0) {
       //cerr << "# Complex atom specifications:\n";
-      for (int j=0; j<multAtoms.size();j++) {
+      for (unsigned int j=0; j<multAtoms.size();j++) {
         cerr << "#\t" << multAtoms[j].id << " -> using ";
-        for (int i=0; i< multAtoms[j].residues.size(); i++) {
+        for (unsigned int i=0; i< multAtoms[j].residues.size(); i++) {
           cerr << multAtoms[j].residues[i] <<":" << multAtoms[j].atoms[i] << ", ";
-          cerr ;
         }
         cerr << " for all other residues: "<< multAtoms[j].defaultAtom << endl;
       }
@@ -179,7 +177,7 @@ void Dscl::readLibReg(std::vector<std::string> &buffer) {
 	    std::istringstream linestream(buffer[i]);
 	    linestream >> regName;
 	    bool outside=false;
-	    for (int j = 0; j < 2*numAng; ++j) { 
+	    for (unsigned int j = 0; j < 2*numAng; ++j) { 
 	      linestream >> regLimits[j];
 	      if (regLimits[j] < periodic[0] || regLimits[j] > periodic[1]) {
 	          outside=true;
@@ -187,7 +185,7 @@ void Dscl::readLibReg(std::vector<std::string> &buffer) {
 	    }
 	    if (outside) {
   	      cerr << "WARNING: region "<< regName << "( ";
-  	      for (int i=0; i < regLimits.size(); i++) {
+  	      for (unsigned int i=0; i < regLimits.size(); i++) {
   	        cerr << regLimits[i] << " ";
   	      }
   	      cerr  <<") outside periodic range ("
@@ -199,8 +197,8 @@ void Dscl::readLibReg(std::vector<std::string> &buffer) {
       }
       
     // check region limits
-    for (int i=0; i< libRegions.size(); i++) {
-      for (int j=0; j<numAng; j=j+2) {
+    for (unsigned int i=0; i< libRegions.size(); i++) {
+      for (unsigned int j=0; j<numAng; j=j+2) {
         if (libRegions[i][j*2] >= libRegions[i][j*2+1]) { 
             cerr << "ERROR: something is wrong with the limits of region "<< 
                      libRegionNames[i]<<"\n";
@@ -214,10 +212,10 @@ void Dscl::readLibReg(std::vector<std::string> &buffer) {
     numReg=libRegions.size();
     
     // test for overlapping region limits
-    for (int i=0; i < numReg; i++) {
-      for (int j=i+1; j < numReg; j++) {
-        int overlap=0;
-        for (int ang=0; ang < numAng; ang++) {
+    for (unsigned int i=0; i < numReg; i++) {
+      for (unsigned int j=i+1; j < numReg; j++) {
+        unsigned int overlap=0;
+        for (unsigned int ang=0; ang < numAng; ang++) {
           double min_i=libRegions[i][ang*2];
           double max_i=libRegions[i][ang*2+1];
           double min_j=libRegions[j][ang*2];
@@ -231,11 +229,11 @@ void Dscl::readLibReg(std::vector<std::string> &buffer) {
         if (overlap == numAng) {
           stringstream errstring;
           errstring << libRegionNames[i] << " and " << libRegionNames[j] << "\n";
-          for (int ang=0; ang < numAng; ang++) {
+          for (unsigned int ang=0; ang < numAng; ang++) {
             errstring<<"|"<<libRegions[i][2*ang]<<" "<<libRegions[i][2*ang+1]<<"|";
           }
           errstring<<"\n";
-          for (int ang=0; ang < numAng; ang++) {
+          for (unsigned int ang=0; ang < numAng; ang++) {
             errstring<<"|"<<libRegions[j][2*ang]<<" "<<libRegions[j][2*ang+1]<<"|";
           }
           errstring<<"\n";
@@ -296,7 +294,7 @@ void Dscl::readLibClass(std::vector<std::string> &buffer) {
 
 void Dscl::print_angles() {
   cout << "# ***ANGLES***" << endl;
-  for (int i=0; i < numAng; i++) {
+  for (unsigned int i=0; i < numAng; i++) {
     cout << "# " << libAngNames[i] << " ";
     for (int j=0; j < 4; j++) { cout << libAngAtoms[i][j]  << " "; }
     for (int j=0; j < 4; j++) { cout << libAngResShifts[i][j] << " "; }
@@ -306,9 +304,9 @@ void Dscl::print_angles() {
 
 void Dscl::print_regions() {
   cout << "# ***REGIONS***" << endl;
-  for ( int i=0; i < numReg ; i++) {
+  for (unsigned int i=0; i < numReg ; i++) {
     cout  << "# " << libRegionNames[i];
-    for(int j=0; j < numAng*2; j++) {
+    for(unsigned int j=0; j < numAng*2; j++) {
       cout << "\t" << libRegions[i][j];
     }
     cout << endl;
@@ -330,7 +328,7 @@ void Dscl::writeHeader() {
   for (int i = 0; i < 87; i++) { dihts << "-"; }
   dihts << "\n#" << setw(9) << "Time";
   dihts << setw(8) << "Residue";
-  for(int i = 0; i < numAng; i++)
+  for(unsigned int i = 0; i < numAng; i++)
     dihts << setw(12) << libAngNames[i];
   dihts << "    # " << setw(12) << "Region" << setw(30) << "Class" << endl;
   dihts << "# ";
@@ -359,10 +357,10 @@ void Dscl::determineAtoms(utils::AtomSpecifier &atomSelection, gcore::System &sy
     if (sys.mol(m).topology().numRes()>(maxResShift-minResShift)) {
       vector<string> atomList (sys.mol(m).topology().numRes(), "");
   
-      for (int i=0; i<angAtomsUniq.size(); i++) {
+      for (unsigned int i=0; i<angAtomsUniq.size(); i++) {
         atomExists.insert(mapTypeString(angAtomsUniq[i],atomList));
       }
-      for (int i=0; i<multAtoms.size(); i++) {
+      for (unsigned int i=0; i<multAtoms.size(); i++) {
         atomExists.insert(mapTypeString(multAtoms[i].id,atomList));
       }
     } 
@@ -375,22 +373,22 @@ void Dscl::determineAtoms(utils::AtomSpecifier &atomSelection, gcore::System &sy
   maxResSel.resize(sys.numMolecules(), 0);
   int prevmol = -1;
   
-  for (int i=0; i < atomSelection.size(); i++) {
+  for (unsigned int i=0; i < atomSelection.size(); i++) {
     // do only if molecule has enough residues to calculate all angles
     if (sys.mol(atomSelection.mol(i)).topology().numRes()>(maxResShift-minResShift)) {
       bool found=false;
       // for single atom specs
-      for (int j=0; j < angAtomsUniq.size(); j++) {
+      for (unsigned int j=0; j < angAtomsUniq.size(); j++) {
         if (atomSelection.name(i) == angAtomsUniq[j]) {
           vecAtomExists[atomSelection.mol(i)][angAtomsUniq[j]][atomSelection.resnum(i)]=atomSelection.name(i);
           found=true;
         }
       }
       // for multiatom specs (alternative atoms for different residues)
-      for (int j=0; j < multAtoms.size(); j++) {
+      for (unsigned int j=0; j < multAtoms.size(); j++) {
         multAtom &ma = multAtoms[j];
         bool resfound=false;
-        for (int r=0; r < ma.residues.size(); r++) {
+        for (unsigned int r=0; r < ma.residues.size(); r++) {
           if (atomSelection.resname(i) == ma.residues[r]) {
             resfound=true;
             if (atomSelection.name(i) == ma.atoms[r]) {
@@ -458,8 +456,8 @@ void Dscl::getPropSpec(gcore::System &sys, bound::Boundary * pbc) {
       
       // for each dihedral read from library, test if all atoms of
       // all angles needed for classification exist
-      for (int ang=0; ang < numAng; ang++) {
-        for (int atom=0; atom < libAngAtoms[ang].size(); atom++) {
+      for (unsigned int ang=0; ang < numAng; ang++) {
+        for (unsigned int atom=0; atom < libAngAtoms[ang].size(); atom++) {
           string &atomname = libAngAtoms[ang][atom];
           int const &shiftedResNum = i+libAngResShifts[ang][atom];
           if ((molecule[atomname][shiftedResNum]).empty()) {
@@ -486,7 +484,7 @@ void Dscl::getPropSpec(gcore::System &sys, bound::Boundary * pbc) {
         fragments.back().resNums.push_back(i+1);
         fragments.back().molNums.push_back(molnum+1);
         
-        for (int ang=0; ang < numAng; ang++) {
+        for (unsigned int ang=0; ang < numAng; ang++) {
         vector<int> shifts = libAngResShifts[ang];
         vector<string> atomnames = libAngAtoms[ang];
         // adding +1 because internal numbering of molecules and residues starts at 0
@@ -504,7 +502,7 @@ void Dscl::getPropSpec(gcore::System &sys, bound::Boundary * pbc) {
   }
   numFrags=fragments.size();
   numResTot=0;
-  for (int f =0; f<numFrags; f++ ) {
+  for (unsigned int f =0; f<numFrags; f++ ) {
     fragments[f].numRes=fragments[f].resIds.size();
     numResTot+=fragments[f].numRes;
   }
@@ -516,8 +514,8 @@ void Dscl::getPropSpec(gcore::System &sys, bound::Boundary * pbc) {
 
 
 void Dscl::getAtoms(gcore::System &sys) {
-  for (int f=0; f<fragments.size(); f++) {
-    for (int i=0; i<fragments[f].numRes; i++) {
+  for (unsigned int f=0; f<fragments.size(); f++) {
+    for (unsigned int i=0; i<fragments[f].numRes; i++) {
       AtomSpecifier as(sys);
       stringstream ss;
       ss << fragments[f].molNums[i] << ":res(" << fragments[f].resNums[i] << ":a)";
@@ -530,8 +528,8 @@ void Dscl::getAtoms(gcore::System &sys) {
 }
 
 void Dscl::calcDih() {
- for (int f=0; f<fragments.size(); f++) {
-  for (int ang=0; ang < numAng; ang++) {
+ for (unsigned int f=0; f<fragments.size(); f++) {
+  for (unsigned int ang=0; ang < numAng; ang++) {
     fragments[f].propStore[ang].calc();
   }  
  }
@@ -553,12 +551,12 @@ double Dscl::modulo(double value, int min, int max) {
 
 
 void Dscl::classifyRegions() {
- for (int f = 0; f<fragments.size(); f++) {
+ for (unsigned int f = 0; f<fragments.size(); f++) {
   // convert property-objects to Double to be able to compare them to 
   //library values and move to the given period  (default -180to180)
   fragments[f].torsions.clear();
   fragments[f].torsions.resize(numAng);
-  for (int ang=0; ang < numAng; ang++) {
+  for (unsigned int ang=0; ang < numAng; ang++) {
     for (unsigned int i=0; i < fragments[f].numRes; i++) {
       fragments[f].torsions[ang].push_back(modulo(fragments[f].propStore[ang][i]->getValue().scalar(),periodic[0],periodic[1]));
     }
@@ -570,12 +568,12 @@ void Dscl::classifyRegions() {
   vector<string> regionstmp(fragments[f].numRes);
   string region;
   
-  for ( int i=0; i < fragments[f].numRes; i++) {
+  for (unsigned  int i=0; i < fragments[f].numRes; i++) {
     regionstmp[i]=unclassString;
-    for (int regcnt=0; regcnt < numReg; regcnt++) {
+    for (unsigned int regcnt=0; regcnt < numReg; regcnt++) {
       vector<double> limits=libRegions[regcnt];  
       bool regcheck=true;
-      for (int ang=0; ang < numAng; ang++) {
+      for (unsigned int ang=0; ang < numAng; ang++) {
         if (!(fragments[f].torsions[ang][i] >= limits[2*ang] && fragments[f].torsions[ang][i] < limits[2*ang+1])) {
           regcheck=false;
           break;
@@ -593,12 +591,12 @@ void Dscl::classifyRegions() {
 
 
 void Dscl::classifyClasses(double const &time) {
- for (int f = 0; f<fragments.size(); f++) {
+ for (unsigned int f = 0; f<fragments.size(); f++) {
   fragments[f].classes.clear();
   fragments[f].bfactors.clear();
   string classname;
   // the last residue can never be classified because we need the +1 region
-  for (int i=0; i < fragments[f].numRes-1; i++) {
+  for (unsigned int i=0; i < fragments[f].numRes-1; i++) {
       if (libClassNames[fragments[f].regions[i]+"-"+fragments[f].regions[i+1]].size()) {
         classname=libClassNames[fragments[f].regions[i]+"-"+fragments[f].regions[i+1]][1];
       } else {
@@ -615,9 +613,9 @@ void Dscl::classifyClasses(double const &time) {
 }
 
 void Dscl::getBfactorValues(gcore::System &sys) {
- for (int f = 0; f<fragments.size(); f++) {
-  for (int i=0; i<fragments[f].numRes; i++) {
-    for (int j=0; j<fragments[f].resAtoms[i].size(); j++) {
+ for (unsigned int f = 0; f<fragments.size(); f++) {
+  for (unsigned int i=0; i<fragments[f].numRes; i++) {
+    for (unsigned int j=0; j<fragments[f].resAtoms[i].size(); j++) {
       sys.mol(fragments[f].molNums[i]-1).setBfac(fragments[f].resAtoms[i].atom(j),fragments[f].bfactors[i]);
       
     }
@@ -627,11 +625,11 @@ void Dscl::getBfactorValues(gcore::System &sys) {
 
 
 void Dscl::writeDihTs(double const &time) {
- for (int f = 0; f<fragments.size(); f++) {
-  for (int i=0; i< fragments[f].numRes; i++) {
+ for (unsigned int f = 0; f<fragments.size(); f++) {
+  for (unsigned int i=0; i< fragments[f].numRes; i++) {
     dihts << setw(10) << std::setprecision(2) << time;
     dihts << setw(8) << fragments[f].resIds[i];
-    for (int ang=0; ang<numAng; ang++) {
+    for (unsigned int ang=0; ang<numAng; ang++) {
       dihts << setw(12) << std::setprecision(4) << fragments[f].torsions[ang][i];
     }
     dihts << "    # " << setw(12) << fragments[f].regions[i] << setw(30) << fragments[f].classes[i]<< endl;
@@ -643,7 +641,7 @@ void Dscl::writeDihTs(double const &time) {
 
 void Dscl::initTimeseries() {
   string classname;
-  for (int i=0; i<numClass; i++) { 
+  for (unsigned int i=0; i<numClass; i++) { 
     classname = classShortnUniq[i];
     ofstream* of = new ofstream;
     classts.insert(map<string,ofstream* >::value_type(classname,of));
@@ -656,7 +654,7 @@ void Dscl::initTimeseries() {
 
 void Dscl::closeTimeseries() {
   string classname;
-  for (int i=0; i<numClass; i++) {
+  for (unsigned int i=0; i<numClass; i++) {
     classname = classShortnUniq[i];
     classts[classname]->close();
     //deallocate
@@ -668,7 +666,7 @@ void Dscl::closeTimeseries() {
 void Dscl::initSummary() {
   typedef multimap<string,vector<int> >::value_type mapTypeInt;
   vector<int> counts (numResTot, 0);
-  for (int i=0; i < classShortnUniq.size(); i++) {
+  for (unsigned int i=0; i < classShortnUniq.size(); i++) {
     summary.insert(mapTypeInt(classShortnUniq[i],counts));
   }
   summary.insert(mapTypeInt(unclassString,counts));
@@ -677,7 +675,7 @@ void Dscl::initSummary() {
 string Dscl::pdbTitle() {
   stringstream ss;
   ss << "DISICL out pdb with class information in the B-factor column\n";
-  for (int i=0; i<classNumMap.size(); i++) {
+  for (unsigned int i=0; i<classNumMap.size(); i++) {
     ss << classNumMap[classShortnUniq[i]] << " - " << classShortnUniq[i]  << ";  ";
     
     if ((i % 7) == 6) ss << endl;
@@ -687,8 +685,8 @@ string Dscl::pdbTitle() {
 
 void Dscl::keepStatistics() {
 int rescounter=0;
-for (int f = 0; f<fragments.size(); f++) {
-  for (int i=0; i<fragments[f].numRes-1; i++) {
+for (unsigned int f = 0; f<fragments.size(); f++) {
+  for (unsigned int i=0; i<fragments[f].numRes-1; i++) {
     ++summary[fragments[f].classes[i]][rescounter];
     rescounter++;
   }
@@ -703,19 +701,19 @@ void Dscl::writeStatistics(unsigned int  frameNum, bool do_tser) {
   stats.precision(1);
   
   stats << "# Number of frames: " <<  frameNum << "\n";
-  for (int j=0; j<classShortnUniq.size()-1; j++) {
+  for (unsigned int j=0; j<classShortnUniq.size()-1; j++) {
     stats << "# "<< classShortnUniq[j] << " ... " << classNameMap[classShortnUniq[j]]<< "\n";
   }
   stats << "# " << unclassString<<" ... unclassified\n";
   stats << "# \n";
   
   stats <<"#" << setw(5) << "resi"<< setw(4) << "mol";
-  for (int j=0; j<classShortnUniq.size(); j++) {
+  for (unsigned int j=0; j<classShortnUniq.size(); j++) {
     stats << setw(13) << classShortnUniq[j];
   }
   stats << "\n";
   stats  << "#         ";
-  for (int j=0; j<classShortnUniq.size(); j++) {
+  for (unsigned int j=0; j<classShortnUniq.size(); j++) {
     stats << "      #     %";
   }
   stats <<"\n";
@@ -723,10 +721,10 @@ void Dscl::writeStatistics(unsigned int  frameNum, bool do_tser) {
   // I exclude the last residue of each molecule, as I can not classify a 
   // class without the next residue's region anyways
   int rescounter=0;
-  for (int f = 0; f<fragments.size(); f++) {
-    for (int i=0; i<fragments[f].numRes-1; i++) { 
+  for (unsigned int f = 0; f<fragments.size(); f++) {
+    for (unsigned int i=0; i<fragments[f].numRes-1; i++) { 
       stats << setw(7) << fragments[f].resIds[i] << setw(3) << fragments[f].molNums[i];
-      for (int j=0; j<classShortnUniq.size(); j++) {
+      for (unsigned int j=0; j<classShortnUniq.size(); j++) {
         entry = &summary[classShortnUniq[j]][rescounter];
         stats << setw(7) << *entry << setw(6) << 100*double(*entry)/ frameNum;
         sum[j]+=*entry;
@@ -738,13 +736,13 @@ void Dscl::writeStatistics(unsigned int  frameNum, bool do_tser) {
     
   stats << "# \n";
   stats << "#         ";
-  for (int j=0; j<classShortnUniq.size(); j++) {
+  for (unsigned int j=0; j<classShortnUniq.size(); j++) {
     stats << setw(13) << classShortnUniq[j];
   }
   stats << "\n";
   stats << setw(7) << "# summary:";
   
-  for (int j=0; j<sum.size(); j++) {
+  for (unsigned int j=0; j<sum.size(); j++) {
     stats << setw(7) << sum[j] << setw(6) << 100*double(sum[j])/ frameNum/(rescounter);
   }
   stats << "\n";
@@ -752,18 +750,18 @@ void Dscl::writeStatistics(unsigned int  frameNum, bool do_tser) {
   // averages at the end of the dihedral timeseries file
   if (do_tser) {
     dihts << "# Averages over run: ";
-    for (int i=0; i<numAng; i++) {
+    for (unsigned int i=0; i<numAng; i++) {
       dihts << libAngNames[i] << "(<average> <rmsd> <error estimate>)  " ;
     }
     dihts << "\n"; 
-    for (int f=0; f<fragments.size(); f++) { 
-    for (int i=0; i<fragments[f].numRes; i++) {
+    for (unsigned int f=0; f<fragments.size(); f++) { 
+    for (unsigned int i=0; i<fragments[f].numRes; i++) {
       dihts << "# Residue " << setw(4) << fragments[f].resIds[i] << ": ";
-      for (int j=0; j<numAng; j++) {
+      for (unsigned int j=0; j<numAng; j++) {
         dihts << fragments[f].propNames[j][i] << " ";
       }
       dihts << endl << "# averages:     ";
-      for (int j=0; j<numAng; j++) {
+      for (unsigned int j=0; j<numAng; j++) {
         dihts << fragments[f].propStore[j][i]->average() << " ";
       } 
       dihts << endl;
@@ -781,7 +779,7 @@ void Dscl::writePdbColorLegend() {
   os << "REMARK residue names correspond to class names\n";
   os << "REMARK to be loaded into a visualization program together\n";
   os << "REMARK with the DISICL output pdbs and colored by b-factor\n";
-  for (int i=0; i < classNumMap.size(); i++) {
+  for (unsigned int i=0; i < classNumMap.size(); i++) {
     string resname=classShortnUniq[i].substr(0,3);
     string atomname;
     if (classShortnUniq[i].size() >3) atomname=classShortnUniq[i].substr(3,6);
diff --git a/gromos++/src/utils/Disicl.h b/gromos++/src/utils/Disicl.h
index 48b3873e..13b81fad 100644
--- a/gromos++/src/utils/Disicl.h
+++ b/gromos++/src/utils/Disicl.h
@@ -87,7 +87,7 @@ namespace utils
    */
    
   class Dscl {
-    int numAng, numReg, numClass, numResTot, numFrags;
+    unsigned int numAng, numReg, numClass, numResTot, numFrags;
     int minResShift, maxResShift;
     string const unclassString;
     
@@ -120,7 +120,7 @@ namespace utils
     // a fragment contains the information for each set of consecutive residues 
     // from a single molecule in the atom selection
     struct fragment {
-      int numRes;
+      unsigned int numRes;
       vector<PropertyContainer> propStore;
       vector<vector<string> > propNames;
       vector<int> resIds;
@@ -163,8 +163,8 @@ namespace utils
      * default constructor
      */
     Dscl() : numAng(0), numReg(0), numClass(0), minResShift(9999), maxResShift(-9999), 
-             tsFile("ts_disicl.dat"), statFile("stat_disicl.out"), 
-             unclassString("UC") {} 
+             unclassString("UC"), tsFile("ts_disicl.dat"), statFile("stat_disicl.out")
+             {} 
        
     /**
      * Destructor.
diff --git a/gromos++/src/utils/Dssp.cc b/gromos++/src/utils/Dssp.cc
index 4fd9abb5..089c9d03 100644
--- a/gromos++/src/utils/Dssp.cc
+++ b/gromos++/src/utils/Dssp.cc
@@ -38,7 +38,7 @@ using utils::AtomSpecifier;
 
 void Dssp::determineAtoms(utils::AtomSpecifier &protein) {
   protein.sort();
-  for (int m = 1; m < protein.size(); m++) {
+  for (unsigned int m = 1; m < protein.size(); m++) {
     if (protein.mol(m) == protein.mol(m - 1)) {
       if (protein.name(m - 1) == "N" && protein.name(m) == "H") {
 	d_H.addAtom(protein.mol(m), protein.atom(m));
@@ -58,7 +58,7 @@ void Dssp::calcHintra_init(utils::AtomSpecifier &protein)
   d_pbc = BoundaryParser::boundary(*d_sys, *d_args);  
   //this gather call does not do anything, 'cause we dont have coords...
   //d_pbc -> gather();
-  for(int m=0; m<protein.size(); m++) {
+  for(unsigned int m=0; m<protein.size(); m++) {
     if(protein.name(m)=="CA") {
       d_CA.addAtom(protein.mol(m), protein.atom(m));
     }
@@ -339,9 +339,9 @@ void Dssp::filter_SecStruct()
   classes.push_back(&turn);
   classes.push_back(&Bend);
   
-  for (int c = 0; c < classes.size()-1; c++) {
-    for (int i = 0; i < (int) classes[c]->size(); ++i) {
-      for (int j = c+1; j < classes.size(); j++) {
+  for (unsigned int c = 0; c < classes.size()-1; c++) {
+    for (unsigned int i = 0; i < classes[c]->size(); ++i) {
+      for (unsigned int j = c+1; j < classes.size(); j++) {
         for (iter=classes[j]->begin(); iter!=classes[j]->end(); ++iter){
           if (*iter == (*classes[c])[i]) {
 	        (*classes[j]).erase(iter);
@@ -552,7 +552,7 @@ void Dssp::calcnumres(utils::AtomSpecifier &protein, const System &sys)
     d_resOffSets.push_back(currentResNums);
   }
 
-  for(int i=0, j=i+1; i<protein.size()-1; i++, j++) {
+  for(unsigned int i=0, j=i+1; i<protein.size()-1; i++, j++) {
 
     while (protein.resnum(i) + d_resOffSets[protein.mol(i)]
 	== protein.resnum(j) + d_resOffSets[protein.mol(j)]
diff --git a/gromos++/src/utils/Energy.cc b/gromos++/src/utils/Energy.cc
index ec418c06..b5acd1d8 100644
--- a/gromos++/src/utils/Energy.cc
+++ b/gromos++/src/utils/Energy.cc
@@ -83,7 +83,7 @@ namespace utils {
         ((1 + 2 * d_eps)*(1 + d_kap * d_cut) + d_eps * (d_kap * d_kap * d_cut * d_cut));
 
     // loop over the atoms
-    for (int i = 0; i < d_as->size(); i++) {
+    for (unsigned int i = 0; i < d_as->size(); i++) {
 
       double qi = d_as->charge(i);
       gmath::Vec vi = d_as->pos(i);
@@ -93,7 +93,7 @@ namespace utils {
       d_f_el_s[i] = Vec(0.0, 0.0, 0.0);
 
       // now, loop over the pairlist
-      for (int j = 0; j < d_pl[i].size(); j++) {
+      for (unsigned int j = 0; j < d_pl[i].size(); j++) {
 
         // determine parameters
 
@@ -134,7 +134,7 @@ namespace utils {
 #ifdef OMP
 #pragma omp parallel for reduction(+ : tmp_el, tmp_vdw)
 #endif
-    for (int i = 0; i < d_as->size(); i++) {
+    for (unsigned int i = 0; i < d_as->size(); i++) {
       int mi = d_as->mol(i);
       int ai = d_as->atom(i);
       int gi = d_as->gromosAtom(i);
@@ -153,7 +153,7 @@ namespace utils {
       d_el_s[i] = 0.0;
 
       // now, loop over the pairlist
-      for (int j = 0; j < d_pl[i].size(); j++) {
+      for (unsigned int j = 0; j < d_pl[i].size(); j++) {
         int mj = d_pl[i].mol(j);
         int aj = d_pl[i].atom(j);
         int gj = d_pl[i].gromosAtom(j);
@@ -375,7 +375,7 @@ namespace utils {
     // for all specified atoms, determine all excluded atoms and all third 
     // neighbours
 
-    for (int i = 0; i < d_as->size(); i++) {
+    for (unsigned int i = 0; i < d_as->size(); i++) {
       if (d_as->atom()[i]->type() == spec_virtual)
         throw gromos::Exception("Energy", "Cannot calculate energy for a virtual atom");
 
@@ -694,70 +694,70 @@ namespace utils {
 
   double Energy::vdw() const{
     double e = 0.0;
-    for (int i = 0; i < d_as->size(); i++)
+    for (unsigned int i = 0; i < d_as->size(); i++)
       e += this->vdw(i);
     return e - d_p_vdw;
   }
 
   double Energy::el() const{
     double e = 0.0;
-    for (int i = 0; i < d_as->size(); i++)
+    for (unsigned int i = 0; i < d_as->size(); i++)
       e += this->el(i);
     return e - d_p_el;
   }
 
-  double Energy::vdw(int i) const{
+  double Energy::vdw(unsigned int i) const{
     assert(i < d_as->size());
     return d_vdw_m[i] + d_vdw_s[i];
   }
 
-  double Energy::el(int i) const{
+  double Energy::el(unsigned int i) const{
     assert(i < d_as->size());
     return d_el_m[i] + d_el_s[i];
   }
 
-  double Energy::cov(int i) const{
-    assert(i < int(d_pc->size()));
+  double Energy::cov(unsigned int i) const{
+    assert(i < d_pc->size());
     return d_cov[i];
   }
 
-  double Energy::vdw_m(int i) const{
+  double Energy::vdw_m(unsigned int i) const{
     assert(i < d_as->size());
     return d_vdw_m[i];
   }
 
-  double Energy::vdw_s(int i) const{
+  double Energy::vdw_s(unsigned int i) const{
     assert(i < d_as->size());
     return d_vdw_s[i];
   }
 
-  double Energy::el_m(int i) const{
+  double Energy::el_m(unsigned int i) const{
     assert(i < d_as->size());
     return d_el_m[i];
   }
 
-  double Energy::el_s(int i) const{
+  double Energy::el_s(unsigned int i) const{
     assert(i < d_as->size());
     return d_el_s[i];
   }
 
-  gmath::Vec Energy::f_el(int i) const{
+  gmath::Vec Energy::f_el(unsigned int i) const{
     assert(i < d_as->size());
     return d_f_el_s[i] + d_f_el_m[i];
   }
 
-  gmath::Vec Energy::f_el_m(int i) const{
+  gmath::Vec Energy::f_el_m(unsigned int i) const{
     assert(i < d_as->size());
     return d_f_el_m[i];
   }
 
-  gmath::Vec Energy::f_el_s(int i) const{
+  gmath::Vec Energy::f_el_s(unsigned int i) const{
     assert(i < d_as->size());
     return d_f_el_s[i];
   }
 
   void Energy::calcPairlist() {
-    if (int(d_pl.size()) != d_as->size())
+    if (d_pl.size() != d_as->size())
       throw Energy::Exception(
         " Cannot calculate pairlist without setting atoms first");
     const int size = d_pl.size();
@@ -780,7 +780,7 @@ namespace utils {
   }
 
   void Energy::setPairlistType(string t) {
-    if (int(d_pl.size()) != d_as->size() || d_pl.size() == 0)
+    if (d_pl.size() != d_as->size() || d_pl.size() == 0)
       throw Energy::Exception(
         " Cannot set pairlist type, without setting atoms first");
     for (unsigned int i = 0; i < d_pl.size(); ++i)
diff --git a/gromos++/src/utils/Energy.h b/gromos++/src/utils/Energy.h
index 62324ea4..e717ece6 100644
--- a/gromos++/src/utils/Energy.h
+++ b/gromos++/src/utils/Energy.h
@@ -258,67 +258,67 @@ namespace utils
      * in the AtomSpecifier
      * @param i The i-th atom in the AtomSpecifier
      */
-    double vdw(int i) const;
+    double vdw(unsigned int i) const;
 
     /**
      * Accessor, returns the total electrostatic energy of the i-th atom
      * in the AtomSpecifier
      * @param i The i-th atom in the AtomSpecifier
      */
-    double el(int i) const;
+    double el(unsigned int i) const;
 
     /**
      * Accessor, returns the potential energy corresponding to the i-th 
      * property in the PropertyContainer
      * @param i The i-th property in the PropertyContainer
      */
-    double cov(int i) const;
+    double cov(unsigned int i) const;
 
     /**
      * Accessor, returns the total Vanderwaals interaction of the i-th atom
      * in the AtomSpecifier with other Solute atoms
      * @param i The i-th atom in the AtomSpecifier
      */
-    double vdw_m(int i) const;
+    double vdw_m(unsigned int i) const;
 
     /**
      * Accessor, returns the total Vanderwaals interaction of the i-th atom
      * in the AtomSpecifier with Solvent atoms
      * @param i The i-th atom in the AtomSpecifier
      */
-    double vdw_s(int i) const;
+    double vdw_s(unsigned int i) const;
 
     /**
      * Accessor, returns the total Electrostatic interaction of the i-the atom
      * in the AtomSpecifier with other Soluter atoms
      * @param i The i-th atom in the AtomSpecifier
      */
-    double el_m(int i) const;
+    double el_m(unsigned int i) const;
 
     /**
      * Accessor, returns the total Electrostratic interaction of the i-th atom
      * in the AtomSpecifier with Solvent atoms
      * @param i The i-th atom in the AtomSpecifier
      */
-    double el_s(int i) const;
+    double el_s(unsigned int i) const;
     /**
      * Accessor, returns the total electrostatic force of the i-th atom
      * i n the AtomSpecifier
      * @param i The i-th atom in the AtomSpecifier
      */
-    gmath::Vec f_el(int i) const;
+    gmath::Vec f_el(unsigned int i) const;
     /**
      * Accessor, returns the total electrostatic force of the i-th atom
      * i n the AtomSpecifier with Solute atoms
      * @param i The i-th atom in the AtomSpecifier
      */
-    gmath::Vec f_el_m(int i) const;
+    gmath::Vec f_el_m(unsigned int i) const;
     /**
      * Accessor, returns the total electrostatic force of the i-th atom
      * i n the AtomSpecifier with Solvent atoms
      * @param i The i-th atom in the AtomSpecifier
      */
-    gmath::Vec f_el_s(int i) const;
+    gmath::Vec f_el_s(unsigned int i) const;
     
     // Exception
     /**
diff --git a/gromos++/src/utils/FfExpertGraph.cc b/gromos++/src/utils/FfExpertGraph.cc
index 3a703d1d..a77a4aba 100644
--- a/gromos++/src/utils/FfExpertGraph.cc
+++ b/gromos++/src/utils/FfExpertGraph.cc
@@ -25,10 +25,11 @@
 
 #include "FfExpertGraph.h"
 
+namespace utils {
 /**
  * write a vertex
  */
-std::ostream & utils::operator<<(std::ostream &os, const utils::Vertex & v) {
+std::ostream & operator<<(std::ostream &os, const utils::Vertex & v) {
   os << "V(" << v.id << "," << v.name << ") ";
   return os;
 }
@@ -37,10 +38,11 @@ std::ostream & utils::operator<<(std::ostream &os, const utils::Vertex & v) {
  * write the graph in the graphviz DOT format. 
  * You can polt it with: neato -Tpng file.dot > file.png
  */
-std::ostream & utils::operator<<(std::ostream &os, const utils::FfExpertGraph & g) {
+std::ostream & operator<<(std::ostream &os, const utils::FfExpertGraph & g) {
   g.print(os);
   return os;
 }
+}
 
 void utils::FfExpertGraph::print(std::ostream &os, const utils::Vertex::equality_enum criterion) const {
   os << "graph G  {" << std::endl;
@@ -203,7 +205,6 @@ bool utils::FfExpertGraph::match_neighbors(
  */
 std::vector<utils::Vertex> utils::FfExpertGraph::equal_subgraph(const utils::FfExpertGraph & g, 
         const utils::Vertex::equality_enum criterion, size_t radius) const {
-  typedef std::set<utils::Vertex> v_set;
   // two graphs are equal when all nodes are physically equal and all
   // edges have the same weights
   
@@ -212,7 +213,6 @@ std::vector<utils::Vertex> utils::FfExpertGraph::equal_subgraph(const utils::FfE
   // first loop over the whole graph and try to find the first vertex
   for(std::vector<utils::Vertex>::const_iterator it = vertices().begin(),
           to = vertices().end(); it != to; ++it) {
-    typedef std::set<utils::Vertex> v_set;
     if (utils::FfExpertGraph::match_neighbors(g.vertices()[0], g, *it, *this, criterion, radius, g.vertices()[0], *it)) {
         solutions.push_back(*it);
     }
diff --git a/gromos++/src/utils/FfExpertGraph.h b/gromos++/src/utils/FfExpertGraph.h
index eb15fec2..faeb8902 100644
--- a/gromos++/src/utils/FfExpertGraph.h
+++ b/gromos++/src/utils/FfExpertGraph.h
@@ -200,7 +200,7 @@ namespace utils {
         if (m_vertices[i] == v) {
             a=i;
         }
-        return a;
+      return a;
     }
 
     /**
diff --git a/gromos++/src/utils/Hbond_calc.cc b/gromos++/src/utils/Hbond_calc.cc
index 178ec7ee..0057614b 100644
--- a/gromos++/src/utils/Hbond_calc.cc
+++ b/gromos++/src/utils/Hbond_calc.cc
@@ -135,7 +135,7 @@ void HB_calc::determineAtoms() {
   acceptors.sort();
 
   int m, a;
-  for (int i = 0; i < donors.size(); ++i) {
+  for (unsigned int i = 0; i < donors.size(); ++i) {
     m = donors.mol(i);
     a = donors.atom(i);
 
@@ -161,14 +161,14 @@ void HB_calc::determineAtoms() {
 
   // if there is no B specified, we take A
   if (args->count("DonorAtomsB") <= 0 && args->count("AcceptorAtomsB") <= 0) {
-    for (int i = 0; i < num_A_donors; i++) {
+    for (unsigned int i = 0; i < num_A_donors; i++) {
       donors.addAtomStrict(donors.mol(i), donors.atom(i));
       bound.addAtomStrict(bound.mol(i), bound.atom(i));
       if(donors.mol(i) < 0){
         donors.addSolventType();
       }
     }
-    for (int i = 0; i < num_A_acceptors; i++) {
+    for (unsigned int i = 0; i < num_A_acceptors; i++) {
       acceptors.addAtomStrict(acceptors.mol(i), acceptors.atom(i));
       if(acceptors.mol(i) < 0)
         acceptors.addSolventType();
@@ -197,7 +197,7 @@ void HB_calc::determineAtoms() {
     donor_B.sort();
     acceptor_B.sort();
 
-    for (int i = 0; i < donor_B.size(); ++i) {
+    for (unsigned int i = 0; i < donor_B.size(); ++i) {
       m = donor_B.mol(i);
       a = donor_B.atom(i);
       if (m < 0) {
@@ -215,18 +215,18 @@ void HB_calc::determineAtoms() {
     }
 
     // copy them into the d_donors, d_bound, d_acceptors
-    for (int i = 0; i < donor_B.size(); ++i) {
+    for (unsigned int i = 0; i < donor_B.size(); ++i) {
       donors.addAtomStrict(donor_B.mol(i), donor_B.atom(i));
       if(donors.mol(i) < 0)
         donors.addSolventType();
     }
-    for (int i =0 ;i < bound_B.size() ;++i ){
+    for (unsigned int i =0 ;i < bound_B.size() ;++i ){
       bound.addAtomStrict(bound_B.mol(i), bound_B.atom(i));
       //if(bound.mol(i) < 0)
         //bound.addSolventType();
     }
 
-    for (int i = 0; i < acceptor_B.size(); ++i) {
+    for (unsigned int i = 0; i < acceptor_B.size(); ++i) {
       acceptors.addAtomStrict(acceptor_B.mol(i), acceptor_B.atom(i));
       if(acceptors.mol(i) < 0)
         acceptors.addSolventType();
@@ -247,8 +247,8 @@ void HB_calc::set_subatomspecs(){
     (2) donA & (accB + accAB)
     (3) donB & (accA + accAB)
     */
-    int i=0;
-    for(int j=num_A_donors; j<donors.size(); ++j){ //go through all atoms in donorsB. search for donAB or donB atoms
+    unsigned int i=0;
+    for(unsigned int j=num_A_donors; j<donors.size(); ++j){ //go through all atoms in donorsB. search for donAB or donB atoms
         for(; i<num_A_donors; ++i){ //and search all atoms in donA
             if(donors.atom(i) == donors.atom(j) && donors.mol(i) == donors.mol(j))
                 break;
@@ -263,7 +263,7 @@ void HB_calc::set_subatomspecs(){
         }
         //increase j by 1 and start again
     }
-    for(int j=0; j<num_A_donors; ++j){ //search donA for donA atoms
+    for(unsigned int j=0; j<num_A_donors; ++j){ //search donA for donA atoms
     	for(i=0; i<donAB.size(); ++i) //all atoms that are present in donorsA and donorsB are already in donAB, so every donorA atom that is not in donAB is stored on donA
     		if(donors.atom(donAB[i]) == donors.atom(j) && donors.mol(donAB[i]) == donors.mol(j))
     		     break;
@@ -273,7 +273,7 @@ void HB_calc::set_subatomspecs(){
     }
     i=0;
     //do the same for acceptors
-    for(int j=num_A_acceptors; j<acceptors.size(); ++j){
+    for(unsigned int j=num_A_acceptors; j<acceptors.size(); ++j){
         for(; i<num_A_acceptors; ++i){ //an search all atoms in accA
             if(acceptors.atom(i) == acceptors.atom(j) && acceptors.mol(i) == acceptors.mol(j))
                 break;
@@ -287,7 +287,7 @@ void HB_calc::set_subatomspecs(){
             ++i; //increase by one, since the atomspec is ordererd, the next atom in accB could be this in accA
         }
     }
-    for(int j=0; j<num_A_acceptors; ++j){ //search accA for accA atoms
+    for(unsigned int j=0; j<num_A_acceptors; ++j){ //search accA for accA atoms
         for(i=0; i<accAB.size(); ++i)
             if(acceptors.atom(accAB[i]) == acceptors.atom(j) && acceptors.mol(accAB[i]) == acceptors.mol(j))
                  break;
@@ -324,7 +324,7 @@ void HB_calc::set_subatomspecs(){
 void HB_calc::determineAtomsbymass() {
   bool keep = false;
   //donors
-  for (int i = 0; i < donors.size(); ++i) {
+  for (unsigned int i = 0; i < donors.size(); ++i) {
     keep = false;
     for (unsigned int j = 0; j < mass_hydrogens.size(); ++j) {
       if (donors.mass(i) == mass_hydrogens[j]) {
@@ -353,7 +353,7 @@ void HB_calc::determineAtomsbymass() {
     throw gromos::Exception("hbond","DonorAtomsA and DonorAtomsB do not contain any atoms. Please check your mass file.");
 
   // acceptors
-  for (int i = 0; i < acceptors.size(); ++i) {
+  for (unsigned int i = 0; i < acceptors.size(); ++i) {
     keep = false;
     for (unsigned int j = 0; j < mass_acceptors.size(); ++j) {
       if (acceptors.mass(i) == mass_acceptors[j]) {
@@ -392,7 +392,8 @@ void HB_calc::readframe() {
 
 void HB_calc::set_reduce(){
   if(reduce){
-      int i, size, start;
+      unsigned int i;
+      int size, start;
       if(donors.numSolventAtoms()){ //if there was a solvent specified
           for(i=0; i < donors.size(); ++i) //find the first atom in the atomspecifier that denotes a solvent
             if(donors.mol(i) < 0)
diff --git a/gromos++/src/utils/Hbond_calc.h b/gromos++/src/utils/Hbond_calc.h
index bc0209ec..955a6e72 100644
--- a/gromos++/src/utils/Hbond_calc.h
+++ b/gromos++/src/utils/Hbond_calc.h
@@ -314,7 +314,7 @@ namespace utils {
     utils::AtomSpecifier donors, bound, acceptors;
     std::vector<int> donAB,donA,donB, accAB,accA,accB;
     std::vector<double> mass_hydrogens, mass_acceptors, mass_donors;
-    int frames, num_A_donors, num_A_acceptors, numHb;
+    unsigned int frames, num_A_donors, num_A_acceptors, numHb;
     std::ofstream timeseriesHB, timeseriesHBtot;
     bool reduce, read_time;
     std::vector<int> solv_donor, solv_acc;
@@ -471,7 +471,7 @@ AcceptorAtomsB=1:a @endverbatim
     /**
      * Method which returns the number of atoms in DonorsA and AcceptorsA.
      */
-    void get_num_A(int& a_don, int& a_acc) const{
+    void get_num_A(unsigned int& a_don, unsigned int& a_acc) const{
         a_don = num_A_donors;
         a_acc = num_A_acceptors;
     }
diff --git a/gromos++/src/utils/Hbond_calc_2c.cc b/gromos++/src/utils/Hbond_calc_2c.cc
index 11ad6765..555ddb1a 100644
--- a/gromos++/src/utils/Hbond_calc_2c.cc
+++ b/gromos++/src/utils/Hbond_calc_2c.cc
@@ -76,13 +76,13 @@ void HB2c_calc::calc(int i, int j) {
         if(reduce){
             if(donors.mol(i) < 0 && !solv_donor.empty() ){ //if atom i is from a solvent molecule
                 i = donors.atom(i) % donors.numSolventAtoms(); //get the solvent atom number (0,1,2,...)
-                int t;
+                unsigned int t;
                 for(t = 0; t < solv_donor.size()-1 && i != solv_donor[t]; ++t); //find the atom number of i: this is stored in solv_donor[t]. t = the position in donors
                 i = solv_donor.back() + t; //solv_donor.back stores where the first solvent starts in donors: add position of the solvent in donors to where the first solvent starts in donors
             }
             if(acceptors.mol(j) < 0 && !solv_acc.empty()){ //same for acceptors
                 j = acceptors.atom(j) % acceptors.numSolventAtoms();
-                int t;
+                unsigned int t;
                 for(t = 0; t < solv_acc.size()-1 && j != solv_acc[t]; ++t);
                 j = solv_acc.back() + t;
             }
@@ -115,7 +115,7 @@ void HB2c_calc::merge(utils::HB2c_calc& input, int traj_num){
 }
 
 void HB2c_calc::go_through_cubes(CubeSystem<int>& cubes_donors, CubeSystem<int>& cubes_acceptors){
-    for(int c=0; c<cubes_donors.size(); ++c){ //go through all cubes with donor atoms from A and all cubes with acceptor atoms from B. the size() is the same for cubes_donors and cubes_acceptors
+    for(unsigned int c=0; c<cubes_donors.size(); ++c){ //go through all cubes with donor atoms from A and all cubes with acceptor atoms from B. the size() is the same for cubes_donors and cubes_acceptors
 
         int don_atoms_i=cubes_donors.cube_i_atomlist(c).size(), //donor atoms my cube
             don_atoms_j=cubes_donors.cube_j_atomlist(c).size(), //donor atoms neighbour cube
@@ -149,17 +149,17 @@ void HB2c_calc::calc_hb(CubeSystem<int>& cubes_donors, CubeSystem<int>& cubes_ac
     //donAB & (accA + accB + accAB)
     if(!donAB.empty() && !(accA.empty() && accB.empty() && accAB.empty() ) ){
 		//assign donAB atoms to cubes
-		for (int i = 0; i < donAB.size(); ++i)
+		for (unsigned int i = 0; i < donAB.size(); ++i)
 		    cubes_donors.assign_atom(donAB[i], donors.pos(donAB[i]));
 
 		//assign accA + accB + accAB atoms to cubes
-		for (int i = 0; i < accA.size(); ++i)
+		for (unsigned int i = 0; i < accA.size(); ++i)
 		    cubes_acceptors.assign_atom(accA[i], acceptors.pos(accA[i]));
 
-        for (int i = 0; i < accB.size(); ++i)
+        for (unsigned int i = 0; i < accB.size(); ++i)
 		    cubes_acceptors.assign_atom(accB[i], acceptors.pos(accB[i]));
 
-        for (int i = 0; i < accAB.size(); ++i)
+        for (unsigned int i = 0; i < accAB.size(); ++i)
 		    cubes_acceptors.assign_atom(accAB[i], acceptors.pos(accAB[i]));
 
         go_through_cubes(cubes_donors, cubes_acceptors);
@@ -172,15 +172,15 @@ void HB2c_calc::calc_hb(CubeSystem<int>& cubes_donors, CubeSystem<int>& cubes_ac
         cubes_acceptors.delete_atoms();
 
         //assign donB
-        for (int i = 0; i < donA.size(); ++i)
+        for (unsigned int i = 0; i < donA.size(); ++i)
             cubes_donors.assign_atom(donA[i], donors.pos(donA[i]));
 
         //assign accB atoms to cubes
-        for (int i = 0; i < accB.size(); ++i)
+        for (unsigned int i = 0; i < accB.size(); ++i)
             cubes_acceptors.assign_atom(accB[i], acceptors.pos(accB[i]));
 
         //assign accAB atoms to cubes
-		for (int i = 0; i < accAB.size(); ++i)
+		for (unsigned int i = 0; i < accAB.size(); ++i)
 		    cubes_acceptors.assign_atom(accAB[i], acceptors.pos(accAB[i]));
 
         go_through_cubes(cubes_donors, cubes_acceptors);
@@ -192,14 +192,14 @@ void HB2c_calc::calc_hb(CubeSystem<int>& cubes_donors, CubeSystem<int>& cubes_ac
         cubes_acceptors.delete_atoms();
 
         //assign donB
-        for (int i = 0; i < donB.size(); ++i)
+        for (unsigned int i = 0; i < donB.size(); ++i)
             cubes_donors.assign_atom(donB[i], donors.pos(donB[i]));
 
         //assign accA atoms to cubes
-        for (int i = 0; i < accA.size(); ++i)
+        for (unsigned int i = 0; i < accA.size(); ++i)
             cubes_acceptors.assign_atom(accA[i], acceptors.pos(accA[i]));
         //assign accAB atoms to cubes
-        for (int i = 0; i < accAB.size(); ++i)
+        for (unsigned int i = 0; i < accAB.size(); ++i)
             cubes_acceptors.assign_atom(accAB[i], acceptors.pos(accAB[i]));
 
         go_through_cubes(cubes_donors, cubes_acceptors);
@@ -214,37 +214,37 @@ void HB2c_calc::calc_vac(){
 
     //donAB & (accA + accB + accAB)
     if(!donAB.empty() && !(accA.empty() && accB.empty() && accAB.empty() ) )
-        for (int i = 0; i < donAB.size(); ++i){
+        for (unsigned int i = 0; i < donAB.size(); ++i){
 
-            for (int j = 0; j < accA.size(); ++j)
+            for (unsigned int j = 0; j < accA.size(); ++j)
                 calc(donAB[i], accA[j]);
 
-            for (int j = 0; j < accB.size(); ++j)
+            for (unsigned int j = 0; j < accB.size(); ++j)
                 calc(donAB[i], accB[j]);
 
-            for (int j = 0; j < accAB.size(); ++j)
+            for (unsigned int j = 0; j < accAB.size(); ++j)
                 calc(donAB[i], accAB[j]);
         }
 
     //donA & (accB + accAB)
     if(!donA.empty() && !(accB.empty() && accAB.empty() ))
-        for (int i = 0; i < donA.size(); ++i){
+        for (unsigned int i = 0; i < donA.size(); ++i){
 
-            for (int j = 0; j < accB.size(); ++j)
+            for (unsigned int j = 0; j < accB.size(); ++j)
                 calc(donA[i], accB[j]);
 
-			for (int j = 0; j < accAB.size(); ++j)
+			for (unsigned int j = 0; j < accAB.size(); ++j)
                 calc(donA[i], accAB[j]);
 		}
 
 	//donB & (accA + accAB)
     if(!donB.empty() && !(accA.empty() && accAB.empty() ) )
-        for (int i = 0; i < donB.size(); ++i){
+        for (unsigned int i = 0; i < donB.size(); ++i){
 
-            for (int j = 0; j < accA.size(); ++j)
+            for (unsigned int j = 0; j < accA.size(); ++j)
                 calc(donB[i], accA[j]);
 
-            for (int j = 0; j < accAB.size(); ++j)
+            for (unsigned int j = 0; j < accAB.size(); ++j)
                 calc(donB[i], accAB[j]);
         }
     ts.back().sort_keys();
@@ -271,7 +271,7 @@ void HB2c_calc::printstatistics(bool sort_occ, double higher){
     //write timeseries - numhb and   timeseries- hbindex
 
     double print_time=time_start-time_dt;
-    for (int traj = 0; traj < traj_map.size(); traj++) {
+    for (unsigned int traj = 0; traj < traj_map.size(); traj++) {
       for(TimeseriesContainer::const_iterator it = traj_map[traj].begin(); it != traj_map[traj].end(); ++it){
         if (!read_time) print_time = print_time+time_dt;
         else print_time = it->time();
@@ -301,7 +301,7 @@ void HB2c_calc::printstatistics(bool sort_occ, double higher){
 		for(HB2cContainerIteratorList::const_iterator it = hb_vec.begin(); it!= hb_vec.end(); ++it){
             if((**it).second.num()/(double)frames * 100 < higher) // as soon as the occurence drops below higher value: stop
                 break;
-        	print((**it).first);
+            print((**it).first);
         }
 
     }
diff --git a/gromos++/src/utils/Hbond_calc_3c.cc b/gromos++/src/utils/Hbond_calc_3c.cc
index 5536e2f9..a4b5f057 100644
--- a/gromos++/src/utils/Hbond_calc_3c.cc
+++ b/gromos++/src/utils/Hbond_calc_3c.cc
@@ -111,19 +111,19 @@ void HB3c_calc::calc(int i, int j, const std::vector<int>& k_atoms) {
 
                 if(donors.mol(i) < 0){ //if atom i is from a solvent molecule
                     i = donors.atom(i) % donors.numSolventAtoms(); //get the solvent atom number (0,1,2,...)
-                    int t;
+                    unsigned int t;
                     for(t = 0; t < solv_donor.size()-1 && i != solv_donor[t]; ++t); //find the atom number of i: this is stored in solv_donor[t]. t = the position in donors
                     i = solv_donor.back() + t; //solv_donor.back stores where the first solvent starts in donors: add position of the solvent in donors to where the first solvent starts in donors
                 }
                 if(acceptors.mol(j) < 0){ //same for acceptors
                     j = acceptors.atom(j) % acceptors.numSolventAtoms();
-                    int t;
+                    unsigned int t;
                     for(t = 0; t < solv_acc.size()-1 && j != solv_acc[t]; ++t);
                     j = solv_acc.back() + t;
                 }
                 if(acceptors.mol(k) < 0){
                     k = acceptors.atom(k) % acceptors.numSolventAtoms();
-                    int t;
+                    unsigned int t;
                     for(t = 0; t < solv_acc.size()-1 && k != solv_acc[t]; ++t);
                     k = solv_acc.back() + t;
                 }
@@ -146,7 +146,7 @@ void HB3c_calc::calc(int i, int j, const std::vector<int>& k_atoms) {
 
 void HB3c_calc::go_through_cubes(CubeSystem<int>& cubes_donors, CubeSystem<int>& cubes_acceptors){
 
-    for(int c=0; c<cubes_donors.size(); ++c){ //go through all cubes with donor atoms from A and all cubes with acceptor atoms from B. the size() is the same for cubes_donors and cubes_acceptors
+    for(unsigned int c=0; c<cubes_donors.size(); ++c){ //go through all cubes with donor atoms from A and all cubes with acceptor atoms from B. the size() is the same for cubes_donors and cubes_acceptors
 
         int don_atoms_i=cubes_donors.cube_i_atomlist(c).size(), //donor atoms my cube
             don_atoms_j=cubes_donors.cube_j_atomlist(c).size(), //donor atoms neighbour cube
@@ -154,7 +154,7 @@ void HB3c_calc::go_through_cubes(CubeSystem<int>& cubes_donors, CubeSystem<int>&
             acc_atoms_j=cubes_acceptors.cube_j_atomlist(c).size(); //acceptor atoms neighbour cube
 
         std::vector<int> k_list;
-        if(don_atoms_i && acc_atoms_j || don_atoms_j && acc_atoms_i) //only update if needed
+        if((don_atoms_i && acc_atoms_j) || (don_atoms_j && acc_atoms_i)) //only update if needed
             k_list = cubes_acceptors.neighbour_atomlist(c); //cubes_k_list atoms k list
 
         if(don_atoms_i && acc_atoms_j) {
@@ -191,17 +191,17 @@ void HB3c_calc::calc_hb(CubeSystem<int>& cubes_donors, CubeSystem<int>& cubes_ac
         cubes_acceptors.delete_atoms();
 
 		//assign donAB atoms to cubes
-		for (int i = 0; i < donAB.size(); ++i)
+		for (unsigned int i = 0; i < donAB.size(); ++i)
 		    cubes_donors.assign_atom(donAB[i], donors.pos(donAB[i]));
 
 		//assign accA + accB + accAB atoms to cubes
-		for (int i = 0; i < accA.size(); ++i)
+		for (unsigned int i = 0; i < accA.size(); ++i)
 		    cubes_acceptors.assign_atom(accA[i], acceptors.pos(accA[i]));
 
-        for (int i = 0; i < accB.size(); ++i)
+        for (unsigned int i = 0; i < accB.size(); ++i)
 		    cubes_acceptors.assign_atom(accB[i], acceptors.pos(accB[i]));
 
-        for (int i = 0; i < accAB.size(); ++i)
+        for (unsigned int i = 0; i < accAB.size(); ++i)
 		    cubes_acceptors.assign_atom(accAB[i], acceptors.pos(accAB[i]));
 
         go_through_cubes(cubes_donors, cubes_acceptors);
@@ -214,15 +214,15 @@ void HB3c_calc::calc_hb(CubeSystem<int>& cubes_donors, CubeSystem<int>& cubes_ac
         cubes_acceptors.delete_atoms();
 
         //assign donB
-        for (int i = 0; i < donA.size(); ++i)
+        for (unsigned int i = 0; i < donA.size(); ++i)
             cubes_donors.assign_atom(donA[i], donors.pos(donA[i]));
 
         //assign accB atoms to cubes
-        for (int i = 0; i < accB.size(); ++i)
+        for (unsigned int i = 0; i < accB.size(); ++i)
             cubes_acceptors.assign_atom(accB[i], acceptors.pos(accB[i]));
 
         //assign accAB atoms to cubes
-		for (int i = 0; i < accAB.size(); ++i)
+		for (unsigned int i = 0; i < accAB.size(); ++i)
 		    cubes_acceptors.assign_atom(accAB[i], acceptors.pos(accAB[i]));
 
         go_through_cubes(cubes_donors, cubes_acceptors);
@@ -234,14 +234,14 @@ void HB3c_calc::calc_hb(CubeSystem<int>& cubes_donors, CubeSystem<int>& cubes_ac
         cubes_acceptors.delete_atoms();
 
         //assign donB
-        for (int i = 0; i < donB.size(); ++i)
+        for (unsigned int i = 0; i < donB.size(); ++i)
             cubes_donors.assign_atom(donB[i], donors.pos(donB[i]));
 
         //assign accA atoms to cubes
-        for (int i = 0; i < accA.size(); ++i)
+        for (unsigned int i = 0; i < accA.size(); ++i)
             cubes_acceptors.assign_atom(accA[i], acceptors.pos(accA[i]));
         //assign accAB atoms to cubes
-        for (int i = 0; i < accAB.size(); ++i)
+        for (unsigned int i = 0; i < accAB.size(); ++i)
             cubes_acceptors.assign_atom(accAB[i], acceptors.pos(accAB[i]));
 
         go_through_cubes(cubes_donors, cubes_acceptors);
@@ -263,15 +263,15 @@ void HB3c_calc::calc_vac(){
         k_list.insert(k_list.end(), accB.begin(), accB.end());
         k_list.insert(k_list.end(), accAB.begin(), accAB.end());
 
-        for (int i = 0; i < donAB.size(); ++i){
+        for (unsigned int i = 0; i < donAB.size(); ++i){
 
-            for (int j = 0; j < accA.size(); ++j)
+            for (unsigned int j = 0; j < accA.size(); ++j)
                 calc(donAB[i], accA[j], k_list);
 
-            for (int j = 0; j < accB.size(); ++j)
+            for (unsigned int j = 0; j < accB.size(); ++j)
                 calc(donAB[i], accB[j], k_list);
 
-            for (int j = 0; j < accAB.size(); ++j)
+            for (unsigned int j = 0; j < accAB.size(); ++j)
                 calc(donAB[i], accAB[j], k_list);
         }
     }
@@ -281,12 +281,12 @@ void HB3c_calc::calc_vac(){
     if(!donA.empty() && !(accB.empty() && accAB.empty() )){
         k_list = accB; //assignment deletes everything that was in the vector before
         k_list.insert(k_list.end(), accAB.begin(), accAB.end());
-        for (int i = 0; i < donA.size(); ++i){
+        for (unsigned int i = 0; i < donA.size(); ++i){
 
-            for (int j = 0; j < accB.size(); ++j)
+            for (unsigned int j = 0; j < accB.size(); ++j)
                 calc(donA[i], accB[j], k_list);
 
-			for (int j = 0; j < accAB.size(); ++j)
+			for (unsigned int j = 0; j < accAB.size(); ++j)
                 calc(donA[i], accAB[j], k_list);
 		}
     }
@@ -295,12 +295,12 @@ void HB3c_calc::calc_vac(){
     if(!donB.empty() && !(accA.empty() && accAB.empty() ) ){
         k_list = accA; //assignment deletes everything that was in the vector before
         k_list.insert(k_list.end(), accAB.begin(), accAB.end());
-        for (int i = 0; i < donB.size(); ++i){
+        for (unsigned int i = 0; i < donB.size(); ++i){
 
-            for (int j = 0; j < accA.size(); ++j)
+            for (unsigned int j = 0; j < accA.size(); ++j)
                 calc(donB[i], accA[j], k_list);
 
-            for (int j = 0; j < accAB.size(); ++j)
+            for (unsigned int j = 0; j < accAB.size(); ++j)
                 calc(donB[i], accAB[j], k_list);
         }
     }
@@ -324,7 +324,7 @@ void HB3c_calc::printstatistics(bool sort_occ, double higher){
 
     //write timeseries - numhb and   timeseries- hbindex
     double print_time=time_start-time_dt;
-    for (int traj = 0; traj < traj_map.size(); traj++) {
+    for (unsigned int traj = 0; traj < traj_map.size(); traj++) {
       for(TimeseriesContainer::const_iterator it = traj_map[traj].begin(); it != traj_map[traj].end(); ++it){
         if (!read_time) {
           print_time = print_time+time_dt;
@@ -356,7 +356,7 @@ void HB3c_calc::printstatistics(bool sort_occ, double higher){
 		for(HB3cContainerIteratorList::const_iterator it = hb_vec.begin(); it!= hb_vec.end(); ++it){
             if((**it).second.num()/(double)frames * 100 < higher) // as soon as the occurence drops below higher value: stop
                 break;
-        	print((**it).first);
+            print((**it).first);
         }
     }
 
@@ -408,7 +408,6 @@ void HB3c_calc::print(const Key3c& key){
     int i_d, i_a1, i_a2,i_d2;
     key.get_atom_num(i_d,i_a1,i_d2,i_a2);
     const HB3c& hb3cprint = hb3cc[key]; //= HB3c
-    int occur = hb3cprint.num(); //how many hbonds of this index
 
     cout << setw(7) << hb3cc[key].id() << " ";
 
diff --git a/gromos++/src/utils/Hbond_calc_bridges.cc b/gromos++/src/utils/Hbond_calc_bridges.cc
index ef2f1001..86c3024d 100644
--- a/gromos++/src/utils/Hbond_calc_bridges.cc
+++ b/gromos++/src/utils/Hbond_calc_bridges.cc
@@ -148,38 +148,38 @@ void HB_bridges::calc(const Key2c& key_left, const Key2c& key_right){
 
             //only take solute - solvent - solute bridges:
     if( //A1..D1/D2..A2 and A1..D1-B1/B2-D2..A2
-        bound.atom(d_left) == bound.atom(d_right) && bound.mol(d_left) < 0 && bound.mol(d_right) < 0 && acceptors.mol(a_left) >= 0 && acceptors.mol(a_right) >= 0  || //&& donors.atom(d_left) != donors.atom(d_right) || //only a bound atom (or maybe also the bound donor) acts as bridge
+        (bound.atom(d_left) == bound.atom(d_right) && bound.mol(d_left) < 0 && bound.mol(d_right) < 0 && acceptors.mol(a_left) >= 0 && acceptors.mol(a_right) >= 0)  || //&& donors.atom(d_left) != donors.atom(d_right) || //only a bound atom (or maybe also the bound donor) acts as bridge
         // B1-D1..A1/A2..D2-B2
-        acceptors.atom(a_left) == acceptors.atom(a_right) && acceptors.mol(a_left) < 0 && acceptors.mol(a_right) < 0 && donors.mol(d_left) >= 0 && donors.mol(d_right) >= 0 || // an acceptor atom acts as bridge
+        (acceptors.atom(a_left) == acceptors.atom(a_right) && acceptors.mol(a_left) < 0 && acceptors.mol(a_right) < 0 && donors.mol(d_left) >= 0 && donors.mol(d_right) >= 0) || // an acceptor atom acts as bridge
         // B2-D2..A2/B1-D1..A1
-        bound.atom(d_left) == acceptors.atom(a_right) && bound.mol(d_left) < 0 && acceptors.mol(a_right) < 0 && acceptors.mol(a_left) >= 0 && donors.mol(d_right) >= 0 || //either outer or inner loop donor is the same atom as the inner/outer acceptor
+        (bound.atom(d_left) == acceptors.atom(a_right) && bound.mol(d_left) < 0 && acceptors.mol(a_right) < 0 && acceptors.mol(a_left) >= 0 && donors.mol(d_right) >= 0) || //either outer or inner loop donor is the same atom as the inner/outer acceptor
         // B1-D1..A1/B2-D2..A2
-        bound.atom(d_right) == acceptors.atom(a_left) && bound.mol(d_right) < 0 && acceptors.mol(a_left) < 0 && acceptors.mol(a_right) >= 0 && donors.mol(d_left) >= 0){
+        (bound.atom(d_right) == acceptors.atom(a_left) && bound.mol(d_right) < 0 && acceptors.mol(a_left) < 0 && acceptors.mol(a_right) >= 0 && donors.mol(d_left) >= 0)){
 
             Key3c key(key_left, key_right);
 
             if(reduce){
                 if(donors.mol(d_left) < 0 && !solv_donor.empty() ){ //if atom i is from a solvent molecule
                     d_left = donors.atom(d_left) % donors.numSolventAtoms(); //get the solvent atom number (0,1,2,...)
-                    int t;
+                    unsigned int t;
                     for(t = 0; t < solv_donor.size()-1 && d_left != solv_donor[t]; ++t); //find the atom number of i: this is stored in solv_donor[t]. t = the position in donors
                     d_left = solv_donor.back() + t; //solv_donor.back stores where the first solvent starts in donors: add position of the solvent in donors to where the first solvent starts in donors
                 }
                 if(donors.mol(d_right) < 0 && !solv_donor.empty() ){ //if atom i is from a solvent molecule
                     d_right = donors.atom(d_right) % donors.numSolventAtoms(); //get the solvent atom number (0,1,2,...)
-                    int t;
+                    unsigned int t;
                     for(t = 0; t < solv_donor.size()-1 && d_right != solv_donor[t]; ++t); //find the atom number of i: this is stored in solv_donor[t]. t = the position in donors
                     d_right = solv_donor.back() + t; //solv_donor.back stores where the first solvent starts in donors: add position of the solvent in donors to where the first solvent starts in donors
                 }
                 if(acceptors.mol(a_left) < 0 && !solv_acc.empty()){ //same for acceptors
                     a_left = acceptors.atom(a_left) % acceptors.numSolventAtoms();
-                    int t;
+                    unsigned int t;
                     for(t = 0; t < solv_acc.size()-1 && a_left != solv_acc[t]; ++t);
                     a_left = solv_acc.back() + t;
                 }
                 if(acceptors.mol(a_right) < 0 && !solv_acc.empty()){ //same for acceptors
                     a_right = acceptors.atom(a_right) % acceptors.numSolventAtoms();
-                    int t;
+                    unsigned int t;
                     for(t = 0; t < solv_acc.size()-1 && a_right != solv_acc[t]; ++t);
                     a_right = solv_acc.back() + t;
                 }
@@ -213,7 +213,7 @@ void HB_bridges::printstatistics(bool sort_occ, double higher){
 
     //write timeseries - numhb and   timeseries- hbindex
     double print_time=time_start-time_dt;
-    for (int traj = 0; traj < traj_map.size(); traj++) {
+    for (unsigned int traj = 0; traj < traj_map.size(); traj++) {
       for(TimeseriesContainer::const_iterator it = traj_map[traj].begin(); it != traj_map[traj].end(); ++it){
         if (!read_time)  print_time = print_time+time_dt;
         else print_time = it->time();
@@ -243,7 +243,7 @@ void HB_bridges::printstatistics(bool sort_occ, double higher){
 		for(BridgeContainerIteratorList::const_iterator it = hb_vec.begin(); it!= hb_vec.end(); ++it){
             if((**it).second.num()/(double)frames * 100 < higher) // as soon as the occurence drops below higher value: stop
                 break;
-        	print((**it).first);
+            print((**it).first);
         }
     }
 
diff --git a/gromos++/src/utils/NeutronScattering.cc b/gromos++/src/utils/NeutronScattering.cc
index 96451f40..2184d3da 100644
--- a/gromos++/src/utils/NeutronScattering.cc
+++ b/gromos++/src/utils/NeutronScattering.cc
@@ -106,8 +106,8 @@ namespace utils {
   int NS::getCombinations(void) {
     // make sure there are no (old) combinations in the list now
     d_this->d_comb.clear();
-    for(int c = 0; c < d_this->d_atoms.size(); ++c) {
-      for(int w = 0; w < d_this->d_atoms.size(); ++w) {
+    for(unsigned int c = 0; c < d_this->d_atoms.size(); ++c) {
+      for(unsigned int w = 0; w < d_this->d_atoms.size(); ++w) {
         int iacc = d_this->d_atoms.iac(c); // has to be in this loop since it
                                            // need to be reset in case of
                                            // iacw > iacc, which is switched
@@ -169,7 +169,7 @@ namespace utils {
     map<int, set<int> >::iterator it;
     for(it = comb.begin(); it != comb.end(); ++it) {
       int num = 0;
-      for(int a = 0; a < d_this->d_atoms.size(); ++a) {
+      for(unsigned int a = 0; a < d_this->d_atoms.size(); ++a) {
         if(it->first == d_this->d_atoms.iac(a)) {
           num++;
         }
@@ -198,7 +198,7 @@ namespace utils {
     int i = 0;
     multimap<int, int>::iterator it;
     for(it = d_this->d_comb.begin(); it != d_this->d_comb.end(); ++it) {
-      for(int asn = 0; asn < d_this->d_atoms.size(); ++asn) {
+      for(unsigned int asn = 0; asn < d_this->d_atoms.size(); ++asn) {
         if (d_this->d_atoms.iac(asn) == it->first) {
           int m = d_this->d_atoms.mol(asn);
           int a = d_this->d_atoms.atom(asn);
@@ -243,15 +243,15 @@ namespace utils {
   void NS::getWeights(void) {
     assert(d_this != NULL);
     assert(d_this->d_atoms.size() != 0);
-    int i = 0;
+    unsigned int i = 0;
     multimap<int, int>::iterator it;
     for (it = d_this->d_comb.begin(); it != d_this->d_comb.end(); ++it) {
       // make sure there are no old weights
       d_this->d_weightIntra[i] = 0;
       d_this->d_weightInter[i] = 0;
       // get the intramolecular weights
-      for (int c = 0; c < d_this->d_atoms.size(); ++c) {
-        for (int w = 0; w < d_this->d_atoms.size(); ++w) {
+      for (unsigned int c = 0; c < d_this->d_atoms.size(); ++c) {
+        for (unsigned int w = 0; w < d_this->d_atoms.size(); ++w) {
           int iacc = d_this->d_atoms.iac(c); // has to be in this loop since it
                                              // need to be reset in case of
                                              // iacw > iacc, which is switched
diff --git a/gromos++/src/utils/Property.cc b/gromos++/src/utils/Property.cc
index e5cc7930..31e6c86a 100644
--- a/gromos++/src/utils/Property.cc
+++ b/gromos++/src/utils/Property.cc
@@ -88,7 +88,7 @@ namespace utils {
     if (atoms().size()) {
       int the_mol = atoms().mol(0);
 
-      for (int m = 1; m < atoms().size(); ++m)
+      for (unsigned int m = 1; m < atoms().size(); ++m)
         if (atoms().mol(m) != the_mol) return -1;
 
       return findTopologyType(sys.mol(the_mol).topology());
@@ -1030,7 +1030,7 @@ namespace utils {
     
     //displacement vector in z direction and total puckering amplitude Q
     std::vector<double> zeta(6);
-    double Q = util::cremerpople::calcZeta(R, RpxRdp, zeta);
+    //double Q = util::cremerpople::calcZeta(R, RpxRdp, zeta);
   
     //theta
     double theta = util::cremerpople::calcTheta(zeta);
@@ -1075,7 +1075,7 @@ namespace utils {
 
     //displacement vector in z direction and total puckering amplitude Q
     std::vector<double> zeta(6);
-    double Q = util::cremerpople::calcZeta(R, RpxRdp, zeta);
+    //double Q = util::cremerpople::calcZeta(R, RpxRdp, zeta);
   
     //phi
     double phi = util::cremerpople::calcPhi(zeta);
@@ -1308,7 +1308,7 @@ namespace utils {
     }
     // add the plane to the atoms: this is not really needed but the user may
     // expect something useful if accessing atoms()
-    for (int i = 0; i < atoms1().size(); ++i)
+    for (unsigned int i = 0; i < atoms1().size(); ++i)
       atoms().addAtom(atoms1().mol(i), atoms1().atom(i));
 
     // add the second arguments to the second plane - with subst.
@@ -1320,7 +1320,7 @@ namespace utils {
         msg << " " << atoms2().name(i);
       throw Exception(msg.str());
     }
-    for (int i = 0; i < atoms2().size(); ++i)
+    for (unsigned int i = 0; i < atoms2().size(); ++i)
       atoms().addAtom(atoms2().mol(i), atoms2().atom(i));
 
     if (int(arguments.size()) == REQUIREDARGUMENTS) {
diff --git a/gromos++/src/utils/RDF.cc b/gromos++/src/utils/RDF.cc
index f660eea3..9295e076 100644
--- a/gromos++/src/utils/RDF.cc
+++ b/gromos++/src/utils/RDF.cc
@@ -225,7 +225,7 @@ namespace utils {
 #ifdef OMP
 #pragma omp parallel for
 #endif
-        for (int c = 0; c < d_this->d_centre.size(); c++) {
+        for (unsigned int c = 0; c < d_this->d_centre.size(); c++) {
 
           // the distribution array
           gmath::Distribution dist(0, d_this->d_cut, d_this->d_grid);
@@ -238,7 +238,7 @@ namespace utils {
           if (d_this->d_with.findAtom(d_this->d_centre.mol(c), d_this->d_centre.atom(c))>-1) inwith = 1;
 
           // loop over the width atoms
-          for (int w = 0; w < d_this->d_with.size(); w++) {
+          for (unsigned int w = 0; w < d_this->d_with.size(); w++) {
 
             // only do the calculations if the centre and with atom are not identical
             if (!(d_this->d_with.mol(w) == d_this->d_centre.mol(c) && d_this->d_with.atom(w) == d_this->d_centre.atom(c))) {
@@ -335,7 +335,7 @@ namespace utils {
 #ifdef OMP
 #pragma omp parallel for
 #endif
-        for (int c = 0; c < d_this->d_centre.size(); c++) {
+        for (unsigned int c = 0; c < d_this->d_centre.size(); c++) {
 
           // the distribution array
           gmath::Distribution dist(0, d_this->d_cut, d_this->d_grid);
@@ -348,7 +348,7 @@ namespace utils {
           if (d_this->d_with.findAtom(d_this->d_centre.mol(c), d_this->d_centre.atom(c))>-1) inwith = 1;
 
           // loop over the width atoms
-          for (int w = 0; w < d_this->d_with.size(); w++) {
+          for (unsigned int w = 0; w < d_this->d_with.size(); w++) {
 
             // only do the calculations if the centre and with atom are within different molecules
             if (!(d_this->d_with.mol(w) == d_this->d_centre.mol(c))) {
@@ -432,7 +432,7 @@ namespace utils {
 #ifdef OMP
 #pragma omp parallel for
 #endif
-        for (int c = 0; c < d_this->d_centre.size(); c++) {
+        for (unsigned int c = 0; c < d_this->d_centre.size(); c++) {
 
           // the distribution array
           gmath::Distribution dist(0, d_this->d_cut, d_this->d_grid);
@@ -442,7 +442,7 @@ namespace utils {
           const gmath::Vec & centre_coord = *(d_this->d_centre.coord(c));
 
           // loop over the with atoms
-          for (int w = 0; w < d_this->d_with.size(); w++) {
+          for (unsigned int w = 0; w < d_this->d_with.size(); w++) {
             // only do the calculations if the center and with atom are not identical
             if (!(d_this->d_with.mol(w) == d_this->d_centre.mol(c) && d_this->d_with.atom(w) == d_this->d_centre.atom(c))) {
               const Vec & tmp = pbc->nearestImage(centre_coord, *(d_this->d_with.coord(w)), d_this->d_sys->box());
@@ -451,7 +451,7 @@ namespace utils {
           } /* end of loop over with atoms */
           
           // loop over the center atoms
-          for (int w = 0; w < d_this->d_centre.size(); w++) {
+          for (unsigned int w = 0; w < d_this->d_centre.size(); w++) {
             // only do the calculations if the center and with atom are not identical
             if (!(d_this->d_centre.mol(w) == d_this->d_centre.mol(c) && d_this->d_centre.atom(w) == d_this->d_centre.atom(c))) {
               const Vec & tmp = pbc->nearestImage(centre_coord, *(d_this->d_centre.coord(w)), d_this->d_sys->box());
@@ -549,7 +549,7 @@ namespace utils {
 #ifdef OMP
 #pragma omp parallel for
 #endif
-        for (int c = 0; c < d_this->d_centre.size(); c++) {
+        for (unsigned int c = 0; c < d_this->d_centre.size(); c++) {
 
           // the distribution array
           gmath::Distribution dist(0, d_this->d_cut, d_this->d_grid);
@@ -562,7 +562,7 @@ namespace utils {
           if (d_this->d_with.findAtom(d_this->d_centre.mol(c), d_this->d_centre.atom(c))>-1) inwith = 1;
 
           // loop over the width atoms
-          for (int w = 0; w < d_this->d_with.size(); w++) {
+          for (unsigned int w = 0; w < d_this->d_with.size(); w++) {
 
             // only do the calculations if the centre and with atom are within different molecules
             if (!(d_this->d_with.mol(w) == d_this->d_centre.mol(c))) {
diff --git a/gromos++/src/utils/RdcFuncs.cc b/gromos++/src/utils/RdcFuncs.cc
index 681f9eba..7c2bf31e 100644
--- a/gromos++/src/utils/RdcFuncs.cc
+++ b/gromos++/src/utils/RdcFuncs.cc
@@ -85,7 +85,7 @@ utils::rdcdata_t utils::read_rdc(const vector<string> &buffer, const System &sys
 
     // first get atom numbers
     istringstream is(buffer[jj]);
-    int i, j, k, l;
+    unsigned int i, j, k, l;
     is >> i >> j >> k >> l;
     DEBUG(5, "i,j,k,l:" << i << ", " << j << ", " << k << ", " << l)
     if (i < 1 || j < 1) {
diff --git a/gromos++/src/utils/SimplePairlist.cc b/gromos++/src/utils/SimplePairlist.cc
index fd37c7bb..1c39a0e1 100644
--- a/gromos++/src/utils/SimplePairlist.cc
+++ b/gromos++/src/utils/SimplePairlist.cc
@@ -172,7 +172,7 @@ namespace utils{
     double d2;
     double d2min = cutmin * cutmin;
     gmath::Vec v;
-    for (int b = 0; b < B.size(); b++) {
+    for (unsigned int b = 0; b < B.size(); b++) {
       v = d_pbc->nearestImage(pos_i, B.pos(b), sys()->box());
       d2 = (pos_i - v).abs2();
       if (d2 < d_cut2 && d2 >= d2min) {
@@ -196,7 +196,7 @@ namespace utils{
     double d2;
     double d2min = cutmin * cutmin;
     gmath::Vec v;
-    for (int b = 0; b < B.size(); b++) {
+    for (unsigned int b = 0; b < B.size(); b++) {
       int molNumB = B.mol(b);
       int atomNumB = B.atom(b);
       v = chargeGroupPosition(molNumB, atomNumB);
diff --git a/gromos++/src/utils/SoluteWeightedDistance.cc b/gromos++/src/utils/SoluteWeightedDistance.cc
index e096f3ad..e1454311 100644
--- a/gromos++/src/utils/SoluteWeightedDistance.cc
+++ b/gromos++/src/utils/SoluteWeightedDistance.cc
@@ -146,13 +146,13 @@ namespace utils {
     Distances * const ds[] = {&_fgDistances, &_cgDistances};
     for (unsigned int k = 0; k < 2; k++) {
       ds[k]->clear();
-      for (int j = 0; j < as[k]->size(); j++) {
+      for (unsigned int j = 0; j < as[k]->size(); j++) {
         double sum = 0.0;
         if (_withMeasures || _withWeights) {
           _weights.clear();
         }
         double minDist;
-        for (int i = 0; i < _solute.size(); i++) {
+        for (unsigned int i = 0; i < _solute.size(); i++) {
           gmath::Vec nim_sol = _pbc->nearestImage(as[k]->pos(j), _solute.pos(i), _sys.box());
           gmath::Vec dist = as[k]->pos(j) - nim_sol;
           double r_ij = dist.abs();
@@ -242,7 +242,7 @@ namespace utils {
     
     AtomSpecifier * const as[] = {&_fgSolvent, &_cgSolvent};
     for (unsigned int k = 0; k < 2; k++) {
-      for (int i = 0; i < as[k]->size(); i++) {
+      for (unsigned int i = 0; i < as[k]->size(); i++) {
         title += " ";
         title += as[k]->toString(i);
       }
diff --git a/gromos++/src/utils/TrajArray.cc b/gromos++/src/utils/TrajArray.cc
index 4036517e..557d1476 100644
--- a/gromos++/src/utils/TrajArray.cc
+++ b/gromos++/src/utils/TrajArray.cc
@@ -36,7 +36,8 @@ void TrajArray::store(const gcore::System &sys,
   const unsigned int frameIndex){
 
   int i;
-  int nAtomsMol, molAtomIndex, molIndex;
+  int nAtomsMol, molAtomIndex;
+  int molIndex;
   unsigned int nAtomsSystem = 0;
   double *framePointer;
 
diff --git a/gromos++/src/utils/VirtualAtom.cc b/gromos++/src/utils/VirtualAtom.cc
index 12b38758..9b16d0f5 100644
--- a/gromos++/src/utils/VirtualAtom.cc
+++ b/gromos++/src/utils/VirtualAtom.cc
@@ -41,7 +41,7 @@ class utils::VirtualAtom_i {
 
   int d_orient;
 
-  int d_required_atoms;
+  unsigned int d_required_atoms;
 
     /**
    * calculates the required atoms
@@ -263,7 +263,7 @@ Vec VirtualAtom::pos()const {
     {
       gmath::Vec v = 0;
       
-      for (int i = 0; i < spec.size(); ++i) {
+      for (unsigned int i = 0; i < spec.size(); ++i) {
         
         v += spec.pos(i);
       }
@@ -275,7 +275,7 @@ Vec VirtualAtom::pos()const {
       double m = 0;
       gmath::Vec v = 0;
 
-      for (int i = 0; i < spec.size(); ++i) {
+      for (unsigned int i = 0; i < spec.size(); ++i) {
         v += spec.mass(i) * spec.pos(i);
         m += spec.mass(i);
       }
diff --git a/gromos++/src/utils/make_top.h b/gromos++/src/utils/make_top.h
index 3f13945b..a7480cce 100644
--- a/gromos++/src/utils/make_top.h
+++ b/gromos++/src/utils/make_top.h
@@ -112,8 +112,8 @@ void addSolute(gcore::LinearTopology &lt,
       lt.atoms()[strt+i].setDampingLevel(bb.atom(i).dampingLevel());
       lt.atoms()[strt+i].setDampingPower(bb.atom(i).dampingPower());
       lt.atoms()[strt+i].setPoloffsiteGamma(bb.atom(i).poloffsiteGamma());
-      lt.atoms()[strt+i].setPoloffsiteI(bb.atom(i).poloffsiteI());
-      lt.atoms()[strt+i].setPoloffsiteJ(bb.atom(i).poloffsiteJ());
+      lt.atoms()[strt+i].setPoloffsiteI(bb.atom(i).poloffsiteI()+strt);
+      lt.atoms()[strt+i].setPoloffsiteJ(bb.atom(i).poloffsiteJ()+strt);
     }
 
     //now, bonded interactions
@@ -718,7 +718,7 @@ void setCysteines(gcore::LinearTopology &lt,
         throw gromos::Exception("setCysteines", "Angle connecting the building blocks was not found");
 
     //three dihedrals
-    int removed_d1 = 0, removed_d2 = 0, removed_d3 = 0;
+    int removed_d1 = 0, removed_d2 = 0; // removed_d3 = 0;
     Dihedral di(a, a + 1, a + 2, b + 2);
     for (std::set<gcore::Dihedral>::iterator iter = lt.dihedrals().begin();
             iter != lt.dihedrals().end(); ++iter) {
@@ -755,7 +755,7 @@ void setCysteines(gcore::LinearTopology &lt,
         if ((*iter)[0] == a + 2 && (*iter)[1] == a - 8 && (*iter)[2] == a - 7 && (*iter)[3] == a - 6) {
             di.setType(iter->type());
             lt.dihedrals().erase(iter);
-            removed_d3 = 1;
+            //removed_d3 = 1;
             break;
         }
     }