@@ -759,6 +759,7 @@ namespace Faunus {
759759 Average<double > msqd, msqd_angle, N;
760760 double thresholdsq=0 , dptrans=0 , dprot=0 , angle=0 , _bias=0 ;
761761 size_t bias_rejected=0 ;
762+ bool rotate; // true if cluster should be rotated
762763 Point dir={1 ,1 ,1 }, dp;
763764 std::vector<std::string> names; // names of molecules to be considered
764765 std::vector<int > ids; // molecule id's of molecules to be considered
@@ -804,7 +805,7 @@ namespace Faunus {
804805 * @param first Index of initial molecule (randomly selected)
805806 * @param index w. all molecules clustered around first (first included)
806807 */
807- void findCluster (Tspace &spc, size_t first, std::set<size_t >& cluster) const {
808+ void findCluster (Tspace &spc, size_t first, std::set<size_t >& cluster) {
808809 assert (first < spc.p .size ());
809810 std::set<size_t > pool (index.begin (), index.end ());
810811 assert (pool.count (first)>0 );
@@ -838,11 +839,13 @@ namespace Faunus {
838839 for (auto j : cluster)
839840 if (j>i)
840841 if (spc.geo .sqdist (spc.groups .at (i).cm , spc.groups .at (j).cm )>=max*max)
841- throw std::runtime_error (name+ " : cluster larger than half box length" );
842+ rotate= false ; // skip rotation if cluster larger than half the box length
842843 }
843844
844845 void _move (Change &change) override {
845- if (thresholdsq>0 && !index.empty ()) {
846+ _bias=0 ;
847+ rotate=true ;
848+ if (thresholdsq>0 and not index.empty ()) {
846849 std::set<size_t > cluster; // all group index in cluster
847850 size_t first = *slump.sample (index.begin (), index.end ()); // random molecule (nuclei)
848851 findCluster (spc, first, cluster); // find cluster around first
@@ -852,6 +855,8 @@ namespace Faunus {
852855 d.all =true ;
853856 dp = 0.5 *ranunit (slump).cwiseProduct (dir) * dptrans;
854857 angle = dprot * (slump ()-0.5 );
858+ if (not rotate)
859+ angle=0 ;
855860
856861 Point COM = Geometry::trigoCom (spc, cluster); // cluster center
857862 Eigen::Quaterniond Q;
@@ -860,11 +865,13 @@ namespace Faunus {
860865 for (auto i : cluster) { // loop over molecules in cluster
861866 auto &g = spc.groups [i];
862867
863- Geometry::rotate (g.begin (), g.end (), Q, spc.geo .boundaryFunc , -COM);
864- g.cm = g.cm -COM;
865- spc.geo .boundary (g.cm );
866- g.cm = Q*g.cm +COM;
867- spc.geo .boundary (g.cm );
868+ if (rotate) {
869+ Geometry::rotate (g.begin (), g.end (), Q, spc.geo .boundaryFunc , -COM);
870+ g.cm = g.cm -COM;
871+ spc.geo .boundary (g.cm );
872+ g.cm = Q*g.cm +COM;
873+ spc.geo .boundary (g.cm );
874+ }
868875
869876 g.translate ( dp, spc.geo .boundaryFunc );
870877 d.index =i;
0 commit comments