<commit_before>/*******************************************************************************
GPU OPTIMIZED MONTE CARLO (GOMC) 2.1
Copyright (C) 2016  GOMC Group
A copy of the GNU General Public License can be found in the COPYRIGHT.txt
along with this program, also can be found at <http://www.gnu.org/licenses/>.
********************************************************************************/
#include "Coordinates.h"
#include "TransformMatrix.h"
#include <algorithm>          //For copy
#include <cmath>
#include <cassert>

void Coordinates::InitFromPDB(pdb_setup::Atoms const& atoms)
{
  //Allocate master array and push stuff to it.
  XYZArray::Init(atoms.beta.size());
  //Transfer without creating a new array.
  std::copy(atoms.x.begin(), atoms.x.end(), x);
  std::copy(atoms.y.begin(), atoms.y.end(), y);
  std::copy(atoms.z.begin(), atoms.z.end(), z);

  CheckCoordinate();
  comRef.CalcCOM();
}

void Coordinates::CheckCoordinate()
{
  int p, start, atom, length, stRange, endRange;
  XYZ min, max, diffV;

  for (uint b = 0; b < BOX_TOTAL; b++)
  {
    MoleculeLookup::box_iterator thisMol = molLookRef.BoxBegin(b),
      end = molLookRef.BoxEnd(b), endc = molLookRef.BoxEnd(b);
    //find the min and max coordinate
    stRange = molRef.MolStart(*thisMol);
    --endc;
    endRange = molRef.MolStart(*endc) + molRef.GetKind(*endc).NumAtoms();

    min.x = *std::min_element(x + stRange, x + endRange);
    max.x = *std::max_element(x + stRange, x + endRange);
    min.y = *std::min_element(y + stRange, y + endRange);
    max.y = *std::max_element(y + stRange, y + endRange);
    min.z = *std::min_element(z + stRange, z + endRange);
    max.z = *std::max_element(z + stRange, z + endRange);
    
    //printf("start: %d, end: %d\n", stRange, endRange);
    printf("Minimum coordinates in box %d: x = %8.3f, y = %8.3f, z = %8.3f\n",
	   b+1, min.x, min.y, min.z);
    printf("Maximum coordinates in box %d: x = %8.3f, y = %8.3f, z = %8.3f\n",
	   b+1, max.x, max.y, max.z);

    min = boxDimRef.TransformUnSlant(min, b);
    max = boxDimRef.TransformUnSlant(max, b);
    diffV = max - min;
    //check to see if molecules are in the box or not
    if( diffV.x > 1.5 * boxDimRef.axis.Get(b).x ||
	diffV.y > 1.5 * boxDimRef.axis.Get(b).y ||
	diffV.z > 1.5 * boxDimRef.axis.Get(b).z)
    {
      printf("Molecules are packed outside of the defined box dimension.\n");
      exit(0);
    }

    printf("Wrapping molecules inside the simulation box %d:\n", b+1);
    while (thisMol != end)
    {
       start = molRef.MolStart(*thisMol);
       MoleculeKind const& thisKind = molRef.GetKind(*thisMol);

       for (p = 0; p < thisKind.NumAtoms(); p++)
       {
	  atom = start + p;
	  boxDimRef.WrapPBC(x[atom], y[atom], z[atom], b);
       }
       ++thisMol;
    }

  }

}
<commit_msg>Fix the logic to see if molecule is inside box or not for non-orthogonal box<commit_after>/*******************************************************************************
GPU OPTIMIZED MONTE CARLO (GOMC) 2.1
Copyright (C) 2016  GOMC Group
A copy of the GNU General Public License can be found in the COPYRIGHT.txt
along with this program, also can be found at <http://www.gnu.org/licenses/>.
********************************************************************************/
#include "Coordinates.h"
#include "TransformMatrix.h"
#include <algorithm>          //For copy
#include <cmath>
#include <cassert>

void Coordinates::InitFromPDB(pdb_setup::Atoms const& atoms)
{
  //Allocate master array and push stuff to it.
  XYZArray::Init(atoms.beta.size());
  //Transfer without creating a new array.
  std::copy(atoms.x.begin(), atoms.x.end(), x);
  std::copy(atoms.y.begin(), atoms.y.end(), y);
  std::copy(atoms.z.begin(), atoms.z.end(), z);

  CheckCoordinate();
  comRef.CalcCOM();
}

void Coordinates::CheckCoordinate()
{
  int p, start, atom, length, stRange, endRange;
  XYZ min, max, diffV, diffX, diffY, diffZ;

  for (uint b = 0; b < BOX_TOTAL; b++)
  {
    MoleculeLookup::box_iterator thisMol = molLookRef.BoxBegin(b),
      end = molLookRef.BoxEnd(b), endc = molLookRef.BoxEnd(b);
    //find the min and max coordinate
    stRange = molRef.MolStart(*thisMol);
    --endc;
    endRange = molRef.MolStart(*endc) + molRef.GetKind(*endc).NumAtoms();

    min.x = *std::min_element(x + stRange, x + endRange);
    max.x = *std::max_element(x + stRange, x + endRange);
    min.y = *std::min_element(y + stRange, y + endRange);
    max.y = *std::max_element(y + stRange, y + endRange);
    min.z = *std::min_element(z + stRange, z + endRange);
    max.z = *std::max_element(z + stRange, z + endRange);
    
    //printf("start: %d, end: %d\n", stRange, endRange);
    printf("Minimum coordinates in box %d: x = %8.3f, y = %8.3f, z = %8.3f\n",
	   b+1, min.x, min.y, min.z);
    printf("Maximum coordinates in box %d: x = %8.3f, y = %8.3f, z = %8.3f\n",
	   b+1, max.x, max.y, max.z);

    printf("Wrapping molecules inside the simulation box %d:\n", b+1);
    while (thisMol != end)
    {
       start = molRef.MolStart(*thisMol);
       MoleculeKind const& thisKind = molRef.GetKind(*thisMol);

       for (p = 0; p < thisKind.NumAtoms(); p++)
       {
	  atom = start + p;
	  boxDimRef.WrapPBC(x[atom], y[atom], z[atom], b);
       }

       diffV = comRef.Get(*thisMol);
       diffV = boxDimRef.TransformUnSlant(diffV, b);
       if( diffV.x > 1.5 * boxDimRef.axis.Get(b).x ||
	   diffV.y > 1.5 * boxDimRef.axis.Get(b).y ||
	   diffV.z > 1.5 * boxDimRef.axis.Get(b).z)
       {
	 printf("Molecules %d is packed outside of the defined box dimension.\n", *thisMol);
	 exit(0);
       }

       ++thisMol;
    }

  }

}

{'fraction_non_alphanumeric': 0.21412057620487285, 'fraction_numerical': 0.010670460608216255, 'mean_word_length': 3.356312935708753, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1262133', 'n_tokens_mistral': 2112, 'n_tokens_neox': 2023, 'n_words': 512}