diff --git a/src/colvar/PathMSDBase.cpp b/src/colvar/PathMSDBase.cpp
index f908ae4a47962b7dcca78f6965164f55a63d3793..3f3c852f4fdddd0ff9de301f4c99b6c5b96a6081 100644
--- a/src/colvar/PathMSDBase.cpp
+++ b/src/colvar/PathMSDBase.cpp
@@ -142,6 +142,7 @@ void PathMSDBase::calculate(){
// resize the list to full
if(imgVec.empty()){ // this is the signal that means: recalculate all
imgVec.resize(nframes);
+#pragma simd
for(unsigned i=0;i<nframes;i++){
imgVec[i].property=indexvec[i];
imgVec[i].index=i;
@@ -199,6 +200,7 @@ void PathMSDBase::calculate(){
for(unsigned i=rank;i<imgVec.size();i+=stride){
tmp_distances[i] = msdv[imgVec[i].index].calc_Rot(getPositions(), tmp_derivs, tmp_rotationRefClose[imgVec[i].index], true);
plumed_assert(tmp_derivs.size()==nat);
+#pragma simd
for(unsigned j=0;j<nat;j++) tmp_derivs2[i*nat+j]=tmp_derivs[j];
}
}
@@ -207,6 +209,7 @@ void PathMSDBase::calculate(){
for(unsigned i=rank;i<imgVec.size();i+=stride){
tmp_distances[i] = msdv[imgVec[i].index].calculateWithCloseStructure(getPositions(), tmp_derivs, rotationPosClose, rotationRefClose[imgVec[i].index], drotationPosCloseDrr01, true);
plumed_assert(tmp_derivs.size()==nat);
+#pragma simd
for(unsigned j=0;j<nat;j++) tmp_derivs2[i*nat+j]=tmp_derivs[j];
if (debugClose){
double withclose = tmp_distances[i];
@@ -226,6 +229,7 @@ void PathMSDBase::calculate(){
for(unsigned i=rank;i<imgVec.size();i+=stride){
tmp_distances[i]=msdv[imgVec[i].index].calculate(getPositions(),tmp_derivs,true);
plumed_assert(tmp_derivs.size()==nat);
+#pragma simd
for(unsigned j=0;j<nat;j++) tmp_derivs2[i*nat+j]=tmp_derivs[j];
}
}
@@ -275,12 +279,16 @@ void PathMSDBase::calculate(){
val_z_path->set(-(1./lambda)*std::log(partition));
for(unsigned j=0;j<s_path.size();j++){
// clean up
+#pragma simd
for(unsigned i=0;i< derivs_s.size();i++){derivs_s[i].zero();}
// do the derivative
+#pragma simd
for(const auto & it : imgVec){
double expval=it.similarity;
tmp=lambda*expval*(s_path[j]-it.property[j])/partition;
+#pragma ivdep
for(unsigned i=0;i< derivs_s.size();i++){ derivs_s[i]+=tmp*it.distder[i] ;}
+#pragma ivdep
if(j==0){for(unsigned i=0;i< derivs_z.size();i++){ derivs_z[i]+=it.distder[i]*expval/partition;}}
}
for(unsigned i=0;i< derivs_s.size();i++){
diff --git a/src/tools/RMSD.cpp b/src/tools/RMSD.cpp
index 6d589438bb3b80da2c46a864bdd7c339096460c5..fd2a9cae9eb101c5e59255e4f56f152868322dc6 100644
--- a/src/tools/RMSD.cpp
+++ b/src/tools/RMSD.cpp
@@ -97,7 +97,9 @@ void RMSD::setReference(const vector<Vector> & reference){
plumed_massert(displace.empty(),"you should first clear() an RMSD object, then set a new reference");
align.resize(n,1.0/n);
displace.resize(n,1.0/n);
+#pragma ivdep
for(unsigned i=0;i<n;i++) reference_center+=this->reference[i]*align[i];
+#pragma ivdep
for(unsigned i=0;i<n;i++) this->reference[i]-=reference_center;
reference_center_is_calculated=true;
reference_center_is_removed=true;
@@ -114,8 +116,10 @@ void RMSD::setAlign(const vector<double> & align, bool normalize_weights, bool r
this->align=align;
if(normalize_weights){
double w=0.0;
+#pragma ivdep
for(unsigned i=0;i<n;i++) w+=this->align[i];
double inv=1.0/w;
+#pragma ivdep
for(unsigned i=0;i<n;i++) this->align[i]*=inv;
}
// recalculate the center anyway
@@ -145,8 +149,10 @@ void RMSD::setDisplace(const vector<double> & displace, bool normalize_weights){
plumed_massert(this->displace.size()==displace.size(),"mismatch in dimension of align/displace arrays");
this->displace=displace;
double w=0.0;
+#pragma ivdep
for(unsigned i=0;i<n;i++) w+=this->displace[i];
double inv=1.0/w;
+#pragma ivdep
if(normalize_weights){for(unsigned i=0;i<n;i++) this->displace[i]*=inv;}
}
std::vector<double> RMSD::getDisplace(){
@@ -965,6 +971,7 @@ void RMSDCoreData::doCoreCalc(bool safe,bool alEqDis, bool only_rotation){
Vector cp; cp.zero(); if(!cpositions_is_removed)cp=cpositions;
Vector cr; cr.zero(); if(!creference_is_removed)cr=creference;
// second expensive loop: compute second moments wrt centers
+#pragma ivdep
for(unsigned iat=0;iat<n;iat++){
double w=align[iat];
rr00+=dotProduct(positions[iat]-cp,positions[iat]-cp)*w;
@@ -1075,6 +1082,7 @@ void RMSDCoreData::doCoreCalc(bool safe,bool alEqDis, bool only_rotation){
// calculate rotation matrix derivatives and components distances needed for components only when align!=displacement
if(!alEqDis)ddist_drotation.zero();
+#pragma ivdep
for(unsigned iat=0;iat<n;iat++){
// components differences: this is useful externally
d[iat]=positions[iat]-cp - matmul(rotation,reference[iat]-cr);
@@ -1105,17 +1113,20 @@ double RMSDCoreData::getDistance( bool squared){
else
dist=eigenvals[0]+rr00+rr11;
const unsigned n=static_cast<unsigned int>(reference.size());
+ double localDist = 0;
+#pragma simd reduction(+:localDist)
for(unsigned iat=0;iat<n;iat++){
if(alEqDis){
- if(safe) dist+=align[iat]*modulo2(d[iat]);
+ if(safe) localDist+=align[iat]*modulo2(d[iat]);
} else {
- dist+=displace[iat]*modulo2(d[iat]);
+ localDist+=displace[iat]*modulo2(d[iat]);
}
}
if(!squared){
- dist=sqrt(dist);
+ dist=sqrt(localDist);
distanceIsMSD=false;
}else{
+ dist=localDist;
distanceIsMSD=true;
}
hasDistance=true;
@@ -1136,10 +1147,13 @@ void RMSDCoreData::doCoreCalcWithCloseStructure(bool safe,bool alEqDis, Tensor &
Tensor rotation = matmul(rotationPosClose, rotationRefClose);
+#pragma simd
for (unsigned iat=0; iat<natoms; iat++){
d[iat] = positions[iat] - cp - matmul(rotation, reference[iat]-cr);
+ }
+ if (!alEqDis){
+ for (unsigned iat=0; iat<natoms; iat++) {
//dist = \sum w_i(x_i - cpos - R_xy * R_ay * a_i)
- if (!alEqDis){
ddist_drxy += -2*displace[iat]*extProduct(matmul(d[iat], rotationRefClose), reference[iat]-cr);
}
}
@@ -1166,6 +1180,7 @@ std::vector<Vector> RMSDCoreData::getDDistanceDPositions(){
if(!isInitialized)plumed_merror("getDPositionsDerivatives needs to initialize the coreData first!");
Vector csum;
Vector tmp1,tmp2;
+#pragma ivdep
for(unsigned iat=0;iat<n;iat++){
if(alEqDis){
// there is no need for derivatives of rotation and shift here as it is by construction zero
@@ -1184,6 +1199,7 @@ std::vector<Vector> RMSDCoreData::getDDistanceDPositions(){
}
}
+#pragma ivdep
if(!alEqDis) for(unsigned iat=0;iat<n;iat++){derivatives[iat]= prefactor*(derivatives[iat]+(ddist_dcpositions-csum)*align[iat]); }
return derivatives;
@@ -1206,6 +1222,7 @@ std::vector<Vector> RMSDCoreData::getDDistanceDReference(){
Tensor t_ddist_drr01=ddist_drr01.transpose();
// third expensive loop: derivatives
+#pragma simd
for(unsigned iat=0;iat<n;iat++){
if(alEqDis){
// there is no need for derivatives of rotation and shift here as it is by construction zero
@@ -1225,6 +1242,7 @@ std::vector<Vector> RMSDCoreData::getDDistanceDReference(){
}
}
+#pragma ivdep
if(!alEqDis) for(unsigned iat=0;iat<n;iat++){derivatives[iat]= prefactor*(derivatives[iat]+(ddist_dcreference-csum)*align[iat]);}
return derivatives;