diff --git a/diis.h b/diis.h
index e2fb3c5..80832c7 100644
--- a/diis.h
+++ b/diis.h
@@ -35,17 +35,17 @@ public:
 	DIIS(const int n, const bool core=1);
 	void setup(const int n, const bool core=1);
 	~DIIS();
-	typename LA_traits<U>::normtype extrapolate(T &vec, const U &errvec); //vec is input/output; returns square residual norm
+	typename LA_traits<U>::normtype extrapolate(T &vec, const U &errvec, bool verbose=false); //vec is input/output; returns square residual norm
 	};
 
 template<typename T, typename U>
-DIIS<T,U>::DIIS(const int n, const bool core) : dim(n), incore(core), bmat(n,n)
+DIIS<T,U>::DIIS(const int n, const bool core) : dim(n), incore(core), bmat(n+1,n+1)
 {
 st=incore?NULL: new AuxStorage<Te>;
 errst=incore?NULL: new AuxStorage<Ue>;
 stor= incore? new T[dim] : NULL;
 errstor= incore? new U[dim] : NULL;
-bmat= (Ue)0; for(int i=1; i<n; ++i) bmat(0,i) = (Ue)-1; 
+bmat= (Ue)0; for(int i=1; i<=n; ++i) bmat(0,i) = (Ue)-1; 
 aktdim=cyclicshift=0;
 }
 
@@ -54,12 +54,12 @@ void DIIS<T,U>::setup(const int n, const bool core)
 {
 dim=n;
 incore=core;
-bmat.resize(n);
+bmat.resize(n+1);
 st=incore?NULL: new AuxStorage<Te>;
 errst=incore?NULL: new AuxStorage<Ue>;
 stor= incore? new T[dim] : NULL;
 errstor= incore? new U[dim] : NULL;
-bmat= (Ue)0; for(int i=1; i<n; ++i) bmat(0,i) = (Ue)-1;
+bmat= (Ue)0; for(int i=1; i<=n; ++i) bmat(0,i) = (Ue)-1;
 aktdim=cyclicshift=0;
 }
 
@@ -76,14 +76,14 @@ if(errstor) delete[] errstor;
 
 
 template<typename T, typename U>
-typename LA_traits<U>::normtype DIIS<T,U>::extrapolate(T &vec, const U &errvec)
+typename LA_traits<U>::normtype DIIS<T,U>::extrapolate(T &vec, const U &errvec, bool verbose)
 {
 if(!dim) laerror("attempt to extrapolate from uninitialized DIIS");
 //if dim exceeded, shift 
 if(aktdim==dim)
 	{
 	cyclicshift=(cyclicshift+1)%dim;
-	for(int i=1; i<dim-1; ++i)
+	for(int i=1; i<dim; ++i)
 		for(int j=1; j<=i; ++j)
 			bmat(i,j)=bmat(i+1,j+1);
 	}
@@ -93,8 +93,8 @@ else
 //store vector
 if(incore) 
 	{
-	stor[(aktdim-1+cyclicshift)%dim]=vec;
-	errstor[(aktdim-1+cyclicshift)%dim]=errvec;
+	stor[(aktdim-1+cyclicshift)%dim]|=vec;
+	errstor[(aktdim-1+cyclicshift)%dim]|=errvec;
 	}
 else 
 	{
@@ -107,23 +107,23 @@ if(aktdim==1) return (typename LA_traits<T>::normtype)1000000000;
 
 //calculate overlaps of the new error with old ones
 typename LA_traits<T>::normtype norm=errvec.norm();
-bmat(aktdim-1,aktdim-1)= norm*norm + DIISEPS;
+bmat(aktdim,aktdim)= norm*norm + DIISEPS;
 if(incore)
-	for(int i=1; i<aktdim-1; ++i) 
-		bmat(i,aktdim-1)=errvec.dot(errstor[(i+cyclicshift)%dim]);
+	for(int i=1; i<aktdim; ++i) 
+		bmat(i,aktdim)=errvec.dot(errstor[(i+cyclicshift-1)%dim]);
 else
 	{
 	U tmp = errvec; tmp.copyonwrite(); //copy dimensions
-	for(int i=1; i<aktdim-1; ++i) 
+	for(int i=1; i<aktdim; ++i) 
 		{
-		errst->get(tmp,(i+cyclicshift)%dim);
-		bmat(i,aktdim-1)= errvec.dot(tmp);
+		errst->get(tmp,(i-1+cyclicshift)%dim);
+		bmat(i,aktdim)= errvec.dot(tmp);
 		}
 	}
 
 
 //prepare rhs-solution vector
-NRVec<Ue> rhs(dim);
+NRVec<Ue> rhs(dim+1);
 rhs= (Ue)0; rhs[0]= (Ue)-1;
 
 //solve for coefficients
@@ -131,19 +131,21 @@ rhs= (Ue)0; rhs[0]= (Ue)-1;
 //@@@ explicit solution - cf. remarks in Pulay memorial book
 {
 NRSMat<Te> amat=bmat;
-linear_solve(amat,rhs,NULL,aktdim);
+linear_solve(amat,rhs,NULL,aktdim+1);
 }
+if(verbose) cout <<"DIIS coefficients: "<<rhs<<endl;
 
 //build the new linear combination
 vec.clear();
 if(incore)
-	for(int i=1; i<aktdim; ++i) vec.axpy(rhs[i],stor[(i+cyclicshift)%dim]);
+	for(int i=1; i<=aktdim; ++i) 
+		vec.axpy(rhs[i],stor[(i-1+cyclicshift)%dim]);
 else
         {
 	T tmp=vec; //copy dimensions
-	for(int i=1; i<aktdim; ++i)
+	for(int i=1; i<=aktdim; ++i)
                 {
-                st->get(tmp,(i+cyclicshift)%dim);
+                st->get(tmp,(i-1+cyclicshift)%dim);
 		vec.axpy(rhs[i],tmp);
                 }
         }