diff --git a/src/madness/chem/BSHApply.h b/src/madness/chem/BSHApply.h
index 141b3d6e21b..137481911f3 100644
--- a/src/madness/chem/BSHApply.h
+++ b/src/madness/chem/BSHApply.h
@@ -85,7 +85,7 @@ class BSHApply {
 
 	    Tensor<T> in=inner(world,Vpsi,bra_res);	// no shift here!
 	    Tensor<double> delta_eps(psi.size());
-	    for (int i=0; i<psi.size(); ++i) delta_eps(i)=std::real(in(i))/(norms[i]*norms[i]);
+	    for (size_t i=0; i<psi.size(); ++i) delta_eps(i)=std::real(in(i))/(norms[i]*norms[i]);
 
 	    if (printme) print("orbital energy update",delta_eps);
 	    double cpu1=cpu_time();
@@ -124,8 +124,8 @@ class BSHApply {
 			const Tensor<T> fock1) const {
 
 		// check dimensions
-		bool consistent=(psi.size()==fock1.dim(0));
-		if ((fock1.ndim()==2) and not (psi.size()==fock1.dim(1))) consistent=false;
+   	        bool consistent=(psi.size()==size_t(fock1.dim(0)));
+		if ((fock1.ndim()==2) and not (psi.size()==size_t(fock1.dim(1)))) consistent=false;
 
 		if (not consistent) {
 			print("Fock matrix dimensions",fock1.ndim(), fock1.dim(0), fock1.dim(1));
diff --git a/src/madness/chem/CC2.cc b/src/madness/chem/CC2.cc
index af9db385497..0ffa7e5e10a 100644
--- a/src/madness/chem/CC2.cc
+++ b/src/madness/chem/CC2.cc
@@ -318,7 +318,7 @@ Tensor<double> CC2::enforce_core_valence_separation(const Tensor<double>& fmat)
     if (nemo->get_param().localize_method()=="canon") {
         auto nmo=nemo->get_calc()->amo.size();
         Tensor<double> fmat1(nmo,nmo);
-        for (int i=0; i<nmo; ++i) fmat1(i,i)=nemo->get_calc()->aeps(i);
+        for (size_t i=0; i<nmo; ++i) fmat1(i,i)=nemo->get_calc()->aeps(i);
         return fmat1;
     }
 
@@ -334,7 +334,7 @@ Tensor<double> CC2::enforce_core_valence_separation(const Tensor<double>& fmat)
     //hf->reset_orbitals(lmo);
     nemo->get_calc()->amo=lmo.get_mos();
     nemo->get_calc()->aeps=lmo.get_eps();
-    MADNESS_CHECK(nemo->get_calc()->aeps.size()==nemo->get_calc()->amo.size());
+    MADNESS_CHECK(size_t(nemo->get_calc()->aeps.size())==nemo->get_calc()->amo.size());
     //orbitals_ = nemo->R*nemo->get_calc()->amo;
     //R2orbitals_ = nemo->ncf->square()*nemo->get_calc()->amo;
 
@@ -415,7 +415,7 @@ double CC2::solve_mp2_coupled(Pairs<CCPair>& doubles) {
 //   std::vector<real_function_6d> coupling_constant_term_vec=Pairs<real_function_6d>::pairs2vector(coupling_constant_term,triangular_map);
 
     // transform vector back to Pairs structure
-    for (int i = 0; i < pair_vec.size(); i++) {
+    for (size_t i = 0; i < pair_vec.size(); i++) {
         pair_vec[i].constant_part = result_vec[i];// - coupling_constant_term_vec[i];
         //save(pair_vec[i].constant_part, pair_vec[i].name() + "_const");
         pair_vec[i].constant_part.truncate().reduce_rank();
@@ -503,20 +503,20 @@ double CC2::solve_mp2_coupled(Pairs<CCPair>& doubles) {
             std::vector<real_function_6d> u;
             for (auto p : pair_vec) u.push_back(p.function());
             std::vector<real_function_6d> kain_update = copy(world,solver.update(u, u_update));
-            for (int i=0; i<pair_vec.size(); ++i) {
+            for (size_t i=0; i<pair_vec.size(); ++i) {
                 kain_update[i].truncate().reduce_rank();
                 kain_update[i].print_size("Kain-Update-Function");
                 pair_vec[i].update_u(copy(kain_update[i]));
             }
         } else {
             if (world.rank()==0) std::cout << "Update without KAIN" << std::endl;
-            for (int i=0; i<pair_vec.size(); ++i) {
+            for (size_t i=0; i<pair_vec.size(); ++i) {
                 pair_vec[i].update_u(pair_vec[i].function() - u_update[i]);
             }
         }
 
         // calculate energy and error and update pairs
-        for (int i = 0; i < pair_vec.size(); i++) {
+        for (size_t i = 0; i < pair_vec.size(); i++) {
             //const real_function_6d residue = pair_vec[i].function() - u_update[i];
             const double error = u_update[i].norm2();
             if (world.rank()==0) std::cout << "residual " << pair_vec[i].i << " " << pair_vec[i].j << " " << error << std::endl;
diff --git a/src/madness/chem/CCPotentials.cc b/src/madness/chem/CCPotentials.cc
index dbb56b27814..478e489d44e 100644
--- a/src/madness/chem/CCPotentials.cc
+++ b/src/madness/chem/CCPotentials.cc
@@ -3075,7 +3075,7 @@ CC_vecfunction CCPotentials::make_full_t_intermediate(const CC_vecfunction& tau)
 
     CC_vecfunction result(returntype);
     for (size_t i = 0; i < mo_ket_.size(); i++) {
-        if (i < parameters.freeze()) {
+        if (int(i) < parameters.freeze()) {
             result.insert(i, mo_ket_(i));
         } else {
             CCFunction t(mo_ket_(i).function + tau(i).function, i, MIXED);
diff --git a/src/madness/chem/CCStructures.h b/src/madness/chem/CCStructures.h
index 73857b64b9e..e4eb3ba4cd1 100644
--- a/src/madness/chem/CCStructures.h
+++ b/src/madness/chem/CCStructures.h
@@ -462,7 +462,7 @@ struct Pairs {
 
     static std::vector<T> pairs2vector(const Pairs<T>& argument, const PairVectorMap map) {
         std::vector<T> vector;
-        for (int i=0; i<argument.allpairs.size(); ++i) {
+        for (size_t i=0; i<argument.allpairs.size(); ++i) {
             vector.push_back(argument(map.map[i].first,map.map[i].second));
         }
         return vector;
diff --git a/src/madness/chem/SCF.cc b/src/madness/chem/SCF.cc
index 23523eed637..4c2ed155da2 100644
--- a/src/madness/chem/SCF.cc
+++ b/src/madness/chem/SCF.cc
@@ -1013,7 +1013,7 @@ void SCF::initial_guess(World& world) {
 
             // Get center of charge for madness
             std::vector<double> mad_coc(3, 0);
-            for (int i = 0; i < molecule.natom(); ++i) {
+            for (size_t i = 0; i < molecule.natom(); ++i) {
                 const Atom& atom = molecule.get_atom(i);
                 int charge = atom.atomic_number;
                 mad_coc[0] += atom.x * charge;
@@ -1060,7 +1060,7 @@ void SCF::initial_guess(World& world) {
             // And rotate
             molecule.rotate(q(Slice(0, 2), Slice(0, 2)));
             if (world.rank() == 0 && param.print_level() > 3) print("New MADNESS coordinates:");
-            for (int i = 0; i < molecule.natom(); ++i) {
+            for (size_t i = 0; i < molecule.natom(); ++i) {
                 const Atom& atom = molecule.get_atom(i);
                 if (world.rank() == 0 && param.print_level() > 3)
                     print(atomic_number_to_symbol(atom.atomic_number), atom.x, atom.y, atom.z);
@@ -2038,7 +2038,7 @@ void SCF::solve(World& world) {
     const double dconv = std::max(FunctionDefaults<3>::get_thresh(),
                                   param.dconv());
     const double trantol = vtol / std::min(30.0, double(amo.size()));
-    const double tolloc = 1e-6; // was std::min(1e-6,0.01*dconv) but now trying to avoid unnecessary change // moved to localizer.h
+    //const double tolloc = 1e-6; // was std::min(1e-6,0.01*dconv) but now trying to avoid unnecessary change // moved to localizer.h
     double update_residual = 0.0, bsh_residual = 0.0;
     subspaceT subspace;
     tensorT Q;
diff --git a/src/madness/chem/TDHF.cc b/src/madness/chem/TDHF.cc
index 63f2ef49c89..c892840ea3e 100644
--- a/src/madness/chem/TDHF.cc
+++ b/src/madness/chem/TDHF.cc
@@ -144,7 +144,7 @@ void TDHF::prepare_calculation() {
     } else {
         std::size_t nmo=get_calc()->aeps.size();
         fmat=Tensor<double>(nmo,nmo);
-        for (int i=0; i<nmo; ++i) fmat(i,i)= get_calc()->aeps(i);
+        for (size_t i=0; i<nmo; ++i) fmat(i,i)= get_calc()->aeps(i);
     }
 
     std::size_t nfrozen=Localizer::determine_frozen_orbitals(fmat);
@@ -302,7 +302,7 @@ std::vector<CC_vecfunction> TDHF::solve_cis() const {
     if (parameters.restart()=="iterate" or parameters.restart()=="no_compute") {
         auto excitations_list=parameters.excitations();
         if (excitations_list.empty()) {
-            for (int i=0; i<parameters.nexcitations(); ++i) excitations_list.push_back(i);
+            for (size_t i=0; i<parameters.nexcitations(); ++i) excitations_list.push_back(i);
         }
         for (auto ex : excitations_list) {
             std::string filename= filename_for_roots(ex);
diff --git a/src/madness/chem/mp2.h b/src/madness/chem/mp2.h
index d9dcaee4484..72e13dce929 100644
--- a/src/madness/chem/mp2.h
+++ b/src/madness/chem/mp2.h
@@ -122,7 +122,7 @@ class HartreeFock {
     void reset_orbitals(const MolecularOrbitals<double,3>& mos) {
         nemo_ptr->get_calc()->amo=mos.get_mos();
         nemo_ptr->get_calc()->aeps=mos.get_eps();
-        MADNESS_CHECK(nemo_ptr->get_calc()->aeps.size()==nemo_ptr->get_calc()->amo.size());
+        MADNESS_CHECK(size_t(nemo_ptr->get_calc()->aeps.size())==nemo_ptr->get_calc()->amo.size());
         orbitals_ = nemo_ptr->R*nemo_ptr->get_calc()->amo;
         R2orbitals_ = nemo_ptr->ncf->square()*nemo_ptr->get_calc()->amo;
     }
diff --git a/src/madness/chem/oep.cc b/src/madness/chem/oep.cc
index 1ebc2baf7fd..f9adb8d0b38 100644
--- a/src/madness/chem/oep.cc
+++ b/src/madness/chem/oep.cc
@@ -239,10 +239,10 @@ double OEP::iterate(const std::string model, const vecfuncT& HF_nemo, const tens
 	bool converged=false;
 
 	timer timer1(world,param.print_level()>=3);
-	for (int iter = 0; iter < oep_param.maxiter(); ++iter) {
+	for (size_t iter = 0; iter < oep_param.maxiter(); ++iter) {
 
 	    if (param.do_localize()) {
-	    	for (int i=0; i<KS_nemo.size(); ++i) calc->aeps(i)=KS_Fock(i,i);
+	    	for (size_t i=0; i<KS_nemo.size(); ++i) calc->aeps(i)=KS_Fock(i,i);
 	    	KS_nemo=localize(KS_nemo,param.econv(),iter==0);
 	    	if (param.print_level()>=10) calc->analyze_vectors(world,KS_nemo,calc->aocc,tensorT(),calc->aset);
 	    }
@@ -372,13 +372,13 @@ double OEP::iterate(const std::string model, const vecfuncT& HF_nemo, const tens
 		KS_nemo = nemo_new;
 		timer1.tag("post-process");
 		if (oep_param.saving_amount() >= 3)
-			for (int n=0; n<KS_nemo.size(); ++n) save(KS_nemo[n], "KS_nemo"+stringify(n)+"iter"+stringify(iter));
+			for (size_t n=0; n<KS_nemo.size(); ++n) save(KS_nemo[n], "KS_nemo"+stringify(n)+"iter"+stringify(iter));
 
 		double deltadensity=0.0;
 		converged=check_convergence(energies,oldenergies,bshnorm,deltadensity,param,
 				param.econv(),param.dconv());
 		if (world.rank() == 0) {
-			printf("\nfinished %s iteration %2d at time %8.1fs with energy %12.8f; residual %12.8f\n\n",
+			printf("\nfinished %s iteration %2lu at time %8.1fs with energy %12.8f; residual %12.8f\n\n",
 					model.c_str(), iter, wall_time(), energy,bshnorm);
 		}
 
diff --git a/src/madness/mra/funcimpl.h b/src/madness/mra/funcimpl.h
index 24b04861b1d..399f01d8fb1 100644
--- a/src/madness/mra/funcimpl.h
+++ b/src/madness/mra/funcimpl.h
@@ -2269,6 +2269,8 @@ template<size_t NDIM>
             /// return the norm of the difference of this node and its "mirror" node
             double operator()(typename rangeT::iterator& it) const {
 
+	      // Temporary fix to GCC whining about out of range access for NDIM!=6
+  	      if constexpr(NDIM==6) {
                 const keyT& key = it->first;
                 const nodeT& fnode = it->second;
 
@@ -2318,6 +2320,10 @@ template<size_t NDIM>
                     norm=fnode.coeff().normf();
                 }
                 return norm*norm;
+	      }
+	      else {
+		throw "ONLY FOR DIM 6!";
+	      }
             }
 
             double operator()(double a, double b) const {
diff --git a/src/madness/mra/mra.h b/src/madness/mra/mra.h
index f5ec8c9fbd6..79d3be60031 100644
--- a/src/madness/mra/mra.h
+++ b/src/madness/mra/mra.h
@@ -2116,7 +2116,7 @@ namespace madness {
     	op.reset_timer();
 
 		// will fence here
-        for (int i=0; i<f1.size(); ++i)
+        for (size_t i=0; i<f1.size(); ++i)
             result.get_impl()->recursive_apply(op, f1[i].get_impl().get(),f2[i].get_impl().get(),false);
         world.gop.fence();
 
diff --git a/src/madness/mra/mraimpl.h b/src/madness/mra/mraimpl.h
index d0e120ccc7c..a3e9f8bd49b 100644
--- a/src/madness/mra/mraimpl.h
+++ b/src/madness/mra/mraimpl.h
@@ -402,7 +402,7 @@ namespace madness {
     Tensor<double> FunctionImpl<T,NDIM>::print_plane_local(const int xaxis, const int yaxis, const coordT& el2) {
         coordT x_sim;
         user_to_sim<NDIM>(el2,x_sim);
-        x_sim[2]+=1.e-10;
+        x_sim[0]+=1.e-10;
 
         // dimensions are: (# boxes)(hue, x lo left, y lo left, x hi right, y hi right)
         Tensor<double> plotinfo(coeffs.size(),5);