Merge branch 'center-phases' into v2.5

CHANGES/v2.5.md

@@ -69,6 +69,7 @@ Changes from version 2.4 which are relevant for users:
   - \ref EXTERNAL can now SCALE the input grid. This allows for more flexibility without modifying the grid file.
   - \ref ALPHABETA can now combine dihedrals with different coefficients
   - \ref INCLUDE can now be used also before setup actions.
+  - \ref CENTER can now be computed using trigonometric functions (PHASES) to simplify its calculation with PBCs.
   - Libmatheval is not used anymore. \ref MATHEVAL (and \ref CUSTOM) are still available
     but employ an internal implementation of the lepton library.
     Functions available in libmatheval and absent in the original lepton library have been added so as to have backward compatbility.

regtest/basic/rt71c/Makefile

+include ../../scripts/test.make

regtest/basic/rt71c/centers.reference

+8
+ 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
+X 0.000000 0.000000 0.000000
+X 1.000000 0.000000 0.000000
+X 0.000000 0.000000 -1.000000
+X 0.333333 0.000000 -0.333333
+X 0.333333 0.000000 -0.333333
+X 0.333333 0.000000 -0.333333
+X 0.333333 0.000000 0.000000
+X 0.333333 0.000000 0.000000
+8
+ 4.000000 4.000000 4.000000
+X 10.000000 10.000000 10.000000
+X 11.000000 10.000000 10.000000
+X 10.000000 10.000000 9.000000
+X 10.333333 10.000000 9.666667
+X 10.333333 10.000000 9.666667
+X -1.704833 2.000000 1.704833
+X -1.704833 2.000000 2.000000
+X -1.704833 2.000000 2.000000
+8
+ 100.000000 1.000000 1.000000 1.000000 100.000000 1.000000 1.000000 1.000000 100.000000
+X 1.000000 2.000000 3.000000
+X 4.000000 7.000000 2.000000
+X 1.000000 9.000000 13.000000
+X 2.000000 6.000000 6.000000
+X 2.000000 6.000000 6.000000
+X 1.998141 6.007338 5.943370
+X 1.998671 3.660581 2.666655
+X 1.998671 3.660581 2.666655
+8
+ 100.000000 9.000000 8.000000 7.000000 40.000000 3.000000 2.000000 1.000000 30.000000
+X 1.000000 2.000000 3.000000
+X 4.000000 7.000000 2.000000
+X 1.000000 9.000000 13.000000
+X 2.000000 6.000000 6.000000
+X 2.000000 6.000000 6.000000
+X 1.945398 6.017210 5.062714
+X 1.993191 3.630205 2.665897
+X 1.993191 3.630205 2.665897
+8
+ 1.000000 1.000000 1.000000
+X 0.100000 0.100000 -99.900000
+X 0.200000 0.200000 23.200000
+X 0.300000 1.300000 33.300000
+X 0.200000 0.533333 -14.466667
+X 0.200000 0.200000 -99.800000
+X 0.200000 0.200000 0.200000
+X 0.132829 0.132829 0.132829
+X 0.132829 0.132829 0.132829

regtest/basic/rt71c/config

+type=driver
+arg="--plumed plumed.dat --ixyz trajectory.xyz --dump-forces ff --dump-forces-fmt=%10.6f --dump-full-virial --debug-forces debugforces"

regtest/basic/rt71c/ff.reference

+3
+  0.235702   0.000000  -0.235702   0.000000   0.000000   0.000000  -0.235702   0.000000   0.235702
+X   0.471405   0.000000  -0.471405
+X  -0.235702   0.000000   0.235702
+X  -0.235702   0.000000   0.235702
+3
+  0.208715   0.000000  -0.208715   0.000000   0.000000   0.000000  -0.208715   0.000000   0.208715
+X   0.424264   0.000000  -0.424264
+X  -0.141421   0.000000   0.282843
+X  -0.282843   0.000000   0.141421
+3
+  0.196454   0.788725   0.579314   0.788725   3.166572   2.325831   0.579314   2.325831   1.708312
+X   0.131057   0.530627   0.380547
+X  -0.065339  -0.267311  -0.196897
+X  -0.065718  -0.263316  -0.183650
+3
+  0.193721   0.823164   0.422670   0.823164   3.497813   1.796020   0.422670   1.796020   0.922201
+X   0.120758   0.607062   0.254492
+X  -0.083411  -0.320661  -0.159031
+X  -0.083505  -0.290919   0.016328
+3
+  0.057735   0.057735   0.057735   0.057735   0.057735   0.057735   0.057735   0.057735   0.057735
+X   0.398939   0.398939   0.398939
+X  -0.220528  -0.220528  -0.220528
+X  -0.178411  -0.178411  -0.178411

regtest/basic/rt71c/plumed.dat

+# dummy stuff using during optimization of the code:
+# DEBUG DETAILED_TIMERS
+# gg: GROUP ATOMS=1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3
+# g: GROUP ATOMS=gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg,gg
+
+## NOTE
+## I used this with --debug-forces to debug forces on CENTER PHASES
+## since phases are disabled when PBC matrix is null
+## they clearly cannot agree with numerical difference in that case.
+## GB
+g: GROUP ATOMS=1,2,3
+ce1: CENTER ATOMS=g NOPBC
+ce2: CENTER ATOMS=g 
+ce3: CENTER ATOMS=g PHASES
+ce4: CENTER ATOMS=1,2,3 PHASES WEIGHTS=2,1,0
+ce5: CENTER ATOMS=1,1,2 PHASES
+DUMPATOMS ATOMS=1,2,3,ce1,ce2,ce3,ce4,ce5 FILE=centers
+de:  DISTANCE ATOMS=1,ce3
+dev: BIASVALUE ARG=de
+
+

regtest/basic/rt71c/trajectory.xyz

+3
+0 0 0 
+X 0 0 0
+X 1 0 0
+X 0 0 -1
+3
+4 4 4
+X 10 10 10
+X 11 10 10
+X 10 10 9
+3
+100 1 1 1 100 1 1 1 100
+X 1 2 3
+X 4 7 2
+X 1 9 13
+3
+100 9 8 7 40 3 2 1 30
+X 1 2 3
+X 4 7 2
+X 1 9 13
+3
+1 0 0 0 1 0 0 0 1
+X 0.1 0.1 -99.9
+X 0.2 0.2 23.2
+X 0.3 1.3 33.3

src/vatom/Center.cpp

@@ -51,6 +51,15 @@ In case you want to recover the old behavior you should use the NOPBC flag.
 In that case you need to take care that atoms are in the correct
 periodic image.
 
+\note As an experimental feature, CENTER also supports a keyword PHASES.
+This keyword solves PBCs by computing scaled coordinates and average
+trigonometric functions, similarly to \ref CENTER_OF_MULTICOLVAR.
+Notice that by construction this center position is
+not invariant with respect to rotations of the atoms at fixed cell lattice.
+In addition, for symmetric Bravais lattices, it is not invariant with respect
+to special symmetries. E.g., if you have an hexagonal cell, the center will
+not be invariant with respect to rotations of 120 degrees.
+On the other hand, it might make the treatment of PBC easier in difficult cases.
 
 \par Examples
 
@@ -77,8 +86,11 @@ class Center:
   public ActionWithVirtualAtom
 {
   std::vector<double> weights;
+  std::vector<Tensor> dcenter_sin;
+  std::vector<Tensor> dcenter_cos;
   bool weight_mass;
   bool nopbc;
+  bool phases;
 public:
   explicit Center(const ActionOptions&ao);
   void calculate();
@@ -92,13 +104,15 @@ void Center::registerKeywords(Keywords& keys) {
   keys.add("optional","WEIGHTS","Center is computed as a weighted average.");
   keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
   keys.addFlag("MASS",false,"If set center is mass weighted");
+  keys.addFlag("PHASES",false,"Compute center using trigonometric phases");
 }
 
 Center::Center(const ActionOptions&ao):
   Action(ao),
   ActionWithVirtualAtom(ao),
   weight_mass(false),
-  nopbc(false)
+  nopbc(false),
+  phases(false)
 {
   vector<AtomNumber> atoms;
   parseAtomList("ATOMS",atoms);
@@ -106,6 +120,7 @@ Center::Center(const ActionOptions&ao):
   parseVector("WEIGHTS",weights);
   parseFlag("MASS",weight_mass);
   parseFlag("NOPBC",nopbc);
+  parseFlag("PHASES",phases);
   checkRead();
   log.printf("  of atoms:");
   for(unsigned i=0; i<atoms.size(); ++i) {
@@ -131,7 +146,9 @@ Center::Center(const ActionOptions&ao):
       log.printf("\n");
     }
   }
-  if(nopbc) {
+  if(phases) {
+    log<<"  Phases will be used to take into account PBC\n";
+  } else if(nopbc) {
     log<<"  PBC will be ignored\n";
   } else {
     log<<"  broken molecules will be rebuilt assuming atoms are in the proper order\n";
@@ -142,7 +159,9 @@ Center::Center(const ActionOptions&ao):
 void Center::calculate() {
   Vector pos;
   double mass(0.0);
-  if(!nopbc) makeWhole();
+  const bool dophases=(getPbc().isSet() ? phases : false);
+
+  if(!nopbc && !dophases) makeWhole();
   vector<Tensor> deriv(getNumberOfAtoms());
   for(unsigned i=0; i<getNumberOfAtoms(); i++) mass+=getMass(i);
   if( plumed.getAtoms().chargesWereSet() ) {
@@ -154,16 +173,79 @@ void Center::calculate() {
   }
   double wtot=0.0;
   for(unsigned i=0; i<weights.size(); i++) wtot+=weights[i];
-  for(unsigned i=0; i<getNumberOfAtoms(); i++) {
-    double w=0;
-    if(weight_mass) w=getMass(i)/mass;
-    else w=weights[i]/wtot;
-    pos+=w*getPosition(i);
-    deriv[i]=w*Tensor::identity();
+
+  if(dophases) {
+    dcenter_sin.resize(getNumberOfAtoms());
+    dcenter_cos.resize(getNumberOfAtoms());
+    Vector center_sin;
+    Vector center_cos;
+    Tensor invbox2pi=2*pi*getPbc().getInvBox();
+    Tensor box2pi=getPbc().getBox() / (2*pi);
+    for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
+      double w=0;
+      if(weight_mass) w=getMass(i)/mass;
+      else w=weights[i]/wtot;
+
+// real to scaled
+      const Vector scaled=matmul(getPosition(i),invbox2pi);
+      const Vector ccos(
+        w*std::cos(scaled[0]),
+        w*std::cos(scaled[1]),
+        w*std::cos(scaled[2])
+      );
+      const Vector csin(
+        w*std::sin(scaled[0]),
+        w*std::sin(scaled[1]),
+        w*std::sin(scaled[2])
+      );
+      center_cos+=ccos;
+      center_sin+=csin;
+      for(unsigned l=0; l<3; l++) for(unsigned k=0; k<3; k++) {
+// k over real coordinates
+// l over scaled coordinates
+          dcenter_sin[i][l][k]=ccos[l]*invbox2pi[k][l];
+          dcenter_cos[i][l][k]=-csin[l]*invbox2pi[k][l];
+        }
+    }
+    const Vector c(
+      std::atan2(center_sin[0],center_cos[0]),
+      std::atan2(center_sin[1],center_cos[1]),
+      std::atan2(center_sin[2],center_cos[2])
+    );
+
+// normalization is convenient for doing derivatives later
+    for(unsigned l=0; l<3; l++) {
+      double norm=1.0/(center_sin[l]*center_sin[l]+center_cos[l]*center_cos[l]);
+      center_sin[l]*=norm;
+      center_cos[l]*=norm;
+    }
+
+    for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
+      Tensor dd;
+      for(unsigned l=0; l<3; l++) for(unsigned k=0; k<3; k++) {
+// k over real coordinates
+// l over scaled coordinates
+          dd[l][k]= (center_cos[l]*dcenter_sin[i][l][k] - center_sin[l]*dcenter_cos[i][l][k]);
+        }
+// scaled to real
+      deriv[i]=matmul(dd,box2pi);
+    }
+    setMass(mass);
+    setAtomsDerivatives(deriv);
+// scaled to real
+    setPosition(matmul(c,box2pi));
+  } else {
+    for(unsigned i=0; i<getNumberOfAtoms(); i++) {
+      double w=0;
+      if(weight_mass) w=getMass(i)/mass;
+      else w=weights[i]/wtot;
+      pos+=w*getPosition(i);
+      deriv[i]=w*Tensor::identity();
+    }
+    setPosition(pos);
+    setMass(mass);
+    setAtomsDerivatives(deriv);
   }
-  setPosition(pos);
-  setMass(mass);
-  setAtomsDerivatives(deriv);
 }
 
 }