diff --git a/src/core/Colvar.cpp b/src/core/Colvar.cpp index c636d89beb02ad9b28f66ddaef692c9afb39d5a0..a2b5ed623934821891991c4eb8c22e8cb90aae97 100644 --- a/src/core/Colvar.cpp +++ b/src/core/Colvar.cpp @@ -68,29 +68,28 @@ void Colvar::apply(){ if(!isEnergy){ for(int i=rank;i<getNumberOfComponents();i+=stride){ - if(getPntrToComponent(i)->applyForce( forces )){ - for(unsigned j=0;j<nat;++j){ + if(getPntrToComponent(i)->applyForce(forces)){ + for(unsigned j=0;j<nat;++j){ f[j][0]+=forces[3*j+0]; f[j][1]+=forces[3*j+1]; f[j][2]+=forces[3*j+2]; - } - v(0,0)+=forces[3*nat+0]; - v(0,1)+=forces[3*nat+1]; - v(0,2)+=forces[3*nat+2]; - v(1,0)+=forces[3*nat+3]; - v(1,1)+=forces[3*nat+4]; - v(1,2)+=forces[3*nat+5]; - v(2,0)+=forces[3*nat+6]; - v(2,1)+=forces[3*nat+7]; - v(2,2)+=forces[3*nat+8]; - } - } - comm.Sum(&f[0][0],3*nat); - comm.Sum(&v[0][0],9); - + } + v(0,0)+=forces[3*nat+0]; + v(0,1)+=forces[3*nat+1]; + v(0,2)+=forces[3*nat+2]; + v(1,0)+=forces[3*nat+3]; + v(1,1)+=forces[3*nat+4]; + v(1,2)+=forces[3*nat+5]; + v(2,0)+=forces[3*nat+6]; + v(2,1)+=forces[3*nat+7]; + v(2,2)+=forces[3*nat+8]; + } + } + comm.Sum(&f[0][0],3*f.size()); + comm.Sum(&v[0][0],9); } else if( isEnergy ){ - forces.resize(1); - if( getPntrToComponent(0)->applyForce( forces ) ) modifyForceOnEnergy()+=forces[0]; + forces.resize(1); + if(getPntrToComponent(0)->applyForce(forces)) modifyForceOnEnergy()+=forces[0]; } } diff --git a/src/tools/Matrix.h b/src/tools/Matrix.h index 24382b93cb75458a9374bd4cdfa9662b8a171b57..859eda8dba6594eb12858b76384449abe852e347 100644 --- a/src/tools/Matrix.h +++ b/src/tools/Matrix.h @@ -241,8 +241,8 @@ template <typename T> int diagMat( const Matrix<T>& A, std::vector<double>& eige // of each of them is positive // We can do it because the phase is arbitrary, and helps making // the result reproducible - for(unsigned i=0;i<n;++i) { - unsigned j; + for(int i=0;i<n;++i) { + int j; for(j=0;j<n;j++) if(eigenvecs(i,j)*eigenvecs(i,j)>1e-14) break; if(j<n) if(eigenvecs(i,j)<0.0) for(unsigned j=0;j<n;j++) eigenvecs(i,j)*=-1; } @@ -277,16 +277,16 @@ template <typename T> int pseudoInvert( const Matrix<T>& A, Matrix<double>& pseu if(info!=0) return info; // Compute the tolerance on the singular values ( machine epsilon * number of singular values * maximum singular value ) - double tol; tol=S[0]; for(unsigned i=1;i<nsv;++i){ if( S[i]>tol ){ tol=S[i]; } } tol*=nsv*epsilon; + double tol; tol=S[0]; for(int i=1;i<nsv;++i){ if( S[i]>tol ){ tol=S[i]; } } tol*=nsv*epsilon; // Get the inverses of the singlular values Matrix<double> Si( ncols, nrows ); Si=0.0; - for(unsigned i=0;i<nsv;++i){ if( S[i]>tol ){ Si(i,i)=1./S[i]; }else{ Si(i,i)=0.0; } } + for(int i=0;i<nsv;++i){ if( S[i]>tol ){ Si(i,i)=1./S[i]; }else{ Si(i,i)=0.0; } } // Now extract matrices for pseudoinverse Matrix<double> V( ncols, ncols ), UT( nrows, nrows ), tmp( ncols, nrows ); - k=0; for(unsigned i=0;i<nrows;++i){ for(unsigned j=0;j<nrows;++j){ UT(i,j)=U[k++]; } } - k=0; for(unsigned i=0;i<ncols;++i){ for(unsigned j=0;j<ncols;++j){ V(i,j)=VT[k++]; } } + k=0; for(int i=0;i<nrows;++i){ for(int j=0;j<nrows;++j){ UT(i,j)=U[k++]; } } + k=0; for(int i=0;i<ncols;++i){ for(int j=0;j<ncols;++j){ V(i,j)=VT[k++]; } } // And do matrix algebra to construct the pseudoinverse if( pseudoinverse.rw!=ncols || pseudoinverse.cl!=nrows ) pseudoinverse.resize( ncols, nrows );