From bb253b67e7b6f3782cc0b0b080ddc7f262315329 Mon Sep 17 00:00:00 2001
From: Jana Pazurikova <pazurikova@ics.muni.cz>
Date: Thu, 2 Nov 2017 11:14:31 +0100
Subject: [PATCH] Change intel pragmas to OpenMP v4+ pragmas

---
 src/colvar/PathMSDBase.cpp | 20 ++++++++---------
 src/tools/RMSD.cpp         | 45 +++++++++++++++++++++-----------------
 2 files changed, 35 insertions(+), 30 deletions(-)

diff --git a/src/colvar/PathMSDBase.cpp b/src/colvar/PathMSDBase.cpp
index c244708d5..b8746221e 100644
--- a/src/colvar/PathMSDBase.cpp
+++ b/src/colvar/PathMSDBase.cpp
@@ -139,7 +139,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
+    #pragma omp simd
     for(unsigned i=0; i<nframes; i++) {
       imgVec[i].property=indexvec[i];
       imgVec[i].index=i;
@@ -170,7 +170,6 @@ void PathMSDBase::calculate() {
       // and we need to accurately recalculate for all reference structures
       computeRefClose = true;
       imgVec.resize(nframes);
-#pragma simd
       for(unsigned i=0; i<nframes; i++) {
         imgVec[i].property=indexvec[i];
         imgVec[i].index=i;
@@ -197,7 +196,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
+        #pragma omp simd
         for(unsigned j=0; j<nat; j++) tmp_derivs2[i*nat+j]=tmp_derivs[j];
       }
     }
@@ -206,7 +205,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
+        #pragma omp simd
         for(unsigned j=0; j<nat; j++) tmp_derivs2[i*nat+j]=tmp_derivs[j];
         if (debugClose) {
           double withclose = tmp_distances[i];
@@ -226,7 +225,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
+      #pragma omp simd
       for(unsigned j=0; j<nat; j++) tmp_derivs2[i*nat+j]=tmp_derivs[j];
     }
   }
@@ -275,17 +274,18 @@ 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
+    #pragma omp 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
+      #pragma omp simd
       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;}}
+      if(j==0) {
+        #pragma omp simd
+        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++) {
       setAtomsDerivatives (val_s_path[j],i,derivs_s[i]);
diff --git a/src/tools/RMSD.cpp b/src/tools/RMSD.cpp
index ddbec1b25..c719a9121 100644
--- a/src/tools/RMSD.cpp
+++ b/src/tools/RMSD.cpp
@@ -97,9 +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
+  #pragma omp simd
   for(unsigned i=0; i<n; i++) reference_center+=this->reference[i]*align[i];
-#pragma ivdep
+  #pragma omp simd
   for(unsigned i=0; i<n; i++) this->reference[i]-=reference_center;
   reference_center_is_calculated=true;
   reference_center_is_removed=true;
@@ -116,10 +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
+    #pragma omp simd
     for(unsigned i=0; i<n; i++) w+=this->align[i];
     double inv=1.0/w;
-#pragma ivdep
+    #pragma omp simd
     for(unsigned i=0; i<n; i++) this->align[i]*=inv;
   }
   // recalculate the center anyway
@@ -149,11 +149,13 @@ 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
+  #pragma omp simd
   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;}
+  if(normalize_weights) {
+    #pragma omp simd
+    for(unsigned i=0; i<n; i++) this->displace[i]*=inv;
+  }
 }
 std::vector<double> RMSD::getDisplace() {
   return displace;
@@ -971,7 +973,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
+  #pragma omp simd
   for(unsigned iat=0; iat<n; iat++) {
     double w=align[iat];
     rr00+=dotProduct(positions[iat]-cp,positions[iat]-cp)*w;
@@ -1082,15 +1084,16 @@ 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
+  #pragma omp simd
   for(unsigned iat=0; iat<n; iat++) {
     // components differences: this is useful externally
     d[iat]=positions[iat]-cp - matmul(rotation,reference[iat]-cr);
     //cerr<<"D "<<iat<<" "<<d[iat][0]<<" "<<d[iat][1]<<" "<<d[iat][2]<<"\n";
-
-    // ddist_drotation if needed
-    if(!alEqDis or !only_rotation) ddist_drotation+=-2*displace[iat]*extProduct(d[iat],reference[iat]-cr);
   }
+  // ddist_drotation if needed
+  if(!alEqDis or !only_rotation)
+    for (unsigned iat=0; iat<n; iat++)
+      ddist_drotation+=-2*displace[iat]*extProduct(d[iat],reference[iat]-cr);
 
   if(!alEqDis or !only_rotation) {
     ddist_drr01.zero();
@@ -1112,7 +1115,7 @@ double RMSDCoreData::getDistance( bool squared) {
     dist=eigenvals[0]+rr00+rr11;
   const unsigned n=static_cast<unsigned int>(reference.size());
   double localDist = 0;
-#pragma simd reduction(+:localDist)
+  #pragma omp simd reduction(+:localDist)
   for(unsigned iat=0; iat<n; iat++) {
     if(alEqDis) {
       if(safe) localDist+=align[iat]*modulo2(d[iat]);
@@ -1145,7 +1148,7 @@ void RMSDCoreData::doCoreCalcWithCloseStructure(bool safe,bool alEqDis, Tensor &
 
   Tensor rotation = matmul(rotationPosClose, rotationRefClose);
 
-#pragma simd
+  #pragma omp simd
   for (unsigned iat=0; iat<natoms; iat++) {
     d[iat] = positions[iat] - cp - matmul(rotation, reference[iat]-cr);
   }
@@ -1178,7 +1181,7 @@ std::vector<Vector> RMSDCoreData::getDDistanceDPositions() {
   if(!isInitialized)plumed_merror("getDPositionsDerivatives needs to initialize the coreData first!");
   Vector csum;
   Vector tmp1,tmp2;
-#pragma ivdep
+  #pragma omp 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
@@ -1197,8 +1200,9 @@ 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]); }
+  if(!alEqDis)
+    #pragma omp simd
+    for(unsigned iat=0; iat<n; iat++) {derivatives[iat]= prefactor*(derivatives[iat]+(ddist_dcpositions-csum)*align[iat]); }
 
   return derivatives;
 }
@@ -1220,7 +1224,7 @@ std::vector<Vector>  RMSDCoreData::getDDistanceDReference() {
   Tensor t_ddist_drr01=ddist_drr01.transpose();
 
 // third expensive loop: derivatives
-#pragma simd
+  #pragma omp 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
@@ -1240,8 +1244,9 @@ 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]);}
+  if(!alEqDis)
+    #pragma omp simd
+    for(unsigned iat=0; iat<n; iat++) {derivatives[iat]= prefactor*(derivatives[iat]+(ddist_dcreference-csum)*align[iat]);}
 
   return derivatives;
 }
-- 
GitLab