OOPSE/libmdtools/SimInfo.cpp

#include <stdlib.h>
#include <string.h>
#include <math.h>

#include <iostream>
using namespace std;

#include "SimInfo.hpp"
#define __C
#include "fSimulation.h"
#include "simError.h"

#include "fortranWrappers.hpp"

#include "MatVec3.h"

#ifdef IS_MPI
#include "mpiSimulation.hpp"
#endif

inline double roundMe( double x ){
  return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
}
          
inline double min( double a, double b ){
  return (a < b ) ? a : b;
}

SimInfo* currentInfo;

SimInfo::SimInfo(){

  n_constraints = 0;
  nZconstraints = 0;
  n_oriented = 0;
  n_dipoles = 0;
  ndf = 0;
  ndfRaw = 0;
  nZconstraints = 0;
  the_integrator = NULL;
  setTemp = 0;
  thermalTime = 0.0;
  currentTime = 0.0;
  rCut = 0.0;
  ecr = 0.0;
  est = 0.0;

  haveRcut = 0;
  haveEcr = 0;
  boxIsInit = 0;
  
  resetTime = 1e99;

  orthoRhombic = 0;
  orthoTolerance = 1E-6;
  useInitXSstate = true;

  usePBC = 0;
  useLJ = 0; 
  useSticky = 0;
  useCharges = 0;
  useDipoles = 0;
  useReactionField = 0;
  useGB = 0;
  useEAM = 0;
  useMolecularCutoffs = 0;

  excludes = Exclude::Instance();

  myConfiguration = new SimState();

  has_minimizer = false;
  the_minimizer =NULL;

  ngroup = 0;

  wrapMeSimInfo( this );
}


SimInfo::~SimInfo(){

  delete myConfiguration;

  map<string, GenericData*>::iterator i;
  
  for(i = properties.begin(); i != properties.end(); i++)
    delete (*i).second;
  
}

void SimInfo::setBox(double newBox[3]) {
  
  int i, j;
  double tempMat[3][3];

  for(i=0; i<3; i++) 
    for (j=0; j<3; j++) tempMat[i][j] = 0.0;;

  tempMat[0][0] = newBox[0];
  tempMat[1][1] = newBox[1];
  tempMat[2][2] = newBox[2];

  setBoxM( tempMat );

}

void SimInfo::setBoxM( double theBox[3][3] ){
  
  int i, j;
  double FortranHmat[9]; // to preserve compatibility with Fortran the
                         // ordering in the array is as follows:
                         // [ 0 3 6 ]
                         // [ 1 4 7 ]
                         // [ 2 5 8 ]
  double FortranHmatInv[9]; // the inverted Hmat (for Fortran);

  if( !boxIsInit ) boxIsInit = 1;

  for(i=0; i < 3; i++) 
    for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
  
  calcBoxL();
  calcHmatInv();

  for(i=0; i < 3; i++) {
    for (j=0; j < 3; j++) {
      FortranHmat[3*j + i] = Hmat[i][j];
      FortranHmatInv[3*j + i] = HmatInv[i][j];
    }
  }

  setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic);
 
}
 

void SimInfo::getBoxM (double theBox[3][3]) {

  int i, j;
  for(i=0; i<3; i++) 
    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j];
}


void SimInfo::scaleBox(double scale) {
  double theBox[3][3];
  int i, j;

  // cerr << "Scaling box by " << scale << "\n";

  for(i=0; i<3; i++) 
    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;

  setBoxM(theBox);

}

void SimInfo::calcHmatInv( void ) {
  
  int oldOrtho;
  int i,j;
  double smallDiag;
  double tol;
  double sanity[3][3];

  invertMat3( Hmat, HmatInv );

  // check to see if Hmat is orthorhombic
  
  oldOrtho = orthoRhombic;

  smallDiag = fabs(Hmat[0][0]);
  if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]);
  if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]);
  tol = smallDiag * orthoTolerance;

  orthoRhombic = 1;
  
  for (i = 0; i < 3; i++ ) {
    for (j = 0 ; j < 3; j++) {
      if (i != j) {
        if (orthoRhombic) {
          if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0;
        }        
      }
    }
  }

  if( oldOrtho != orthoRhombic ){
    
    if( orthoRhombic ){
      sprintf( painCave.errMsg,
               "OOPSE is switching from the default Non-Orthorhombic\n"
               "\tto the faster Orthorhombic periodic boundary computations.\n"
               "\tThis is usually a good thing, but if you wan't the\n"
               "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
               "\tvariable ( currently set to %G ) smaller.\n",
               orthoTolerance);
      simError();
    }
    else {
      sprintf( painCave.errMsg,
               "OOPSE is switching from the faster Orthorhombic to the more\n"
               "\tflexible Non-Orthorhombic periodic boundary computations.\n"
               "\tThis is usually because the box has deformed under\n"
               "\tNPTf integration. If you wan't to live on the edge with\n"
               "\tthe Orthorhombic computations, make the orthoBoxTolerance\n"
               "\tvariable ( currently set to %G ) larger.\n",
               orthoTolerance);
      simError();
    }
  }
}

void SimInfo::calcBoxL( void ){

  double dx, dy, dz, dsq;

  // boxVol = Determinant of Hmat

  boxVol = matDet3( Hmat );

  // boxLx
  
  dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
  dsq = dx*dx + dy*dy + dz*dz;
  boxL[0] = sqrt( dsq );
  //maxCutoff = 0.5 * boxL[0];

  // boxLy
  
  dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
  dsq = dx*dx + dy*dy + dz*dz;
  boxL[1] = sqrt( dsq );
  //if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];


  // boxLz
  
  dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
  dsq = dx*dx + dy*dy + dz*dz;
  boxL[2] = sqrt( dsq );
  //if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];

  //calculate the max cutoff
  maxCutoff =  calcMaxCutOff(); 
  
  checkCutOffs();

}


double SimInfo::calcMaxCutOff(){

  double ri[3], rj[3], rk[3];
  double rij[3], rjk[3], rki[3];
  double minDist;

  ri[0] = Hmat[0][0];
  ri[1] = Hmat[1][0];
  ri[2] = Hmat[2][0];

  rj[0] = Hmat[0][1];
  rj[1] = Hmat[1][1];
  rj[2] = Hmat[2][1];

  rk[0] = Hmat[0][2];
  rk[1] = Hmat[1][2];
  rk[2] = Hmat[2][2];
    
  crossProduct3(ri, rj, rij);
  distXY = dotProduct3(rk,rij) / norm3(rij);

  crossProduct3(rj,rk, rjk);
  distYZ = dotProduct3(ri,rjk) / norm3(rjk);

  crossProduct3(rk,ri, rki);
  distZX = dotProduct3(rj,rki) / norm3(rki);

  minDist = min(min(distXY, distYZ), distZX);
  return minDist/2;
  
}

void SimInfo::wrapVector( double thePos[3] ){

  int i;
  double scaled[3];

  if( !orthoRhombic ){
    // calc the scaled coordinates.
  

    matVecMul3(HmatInv, thePos, scaled);
    
    for(i=0; i<3; i++)
      scaled[i] -= roundMe(scaled[i]);
    
    // calc the wrapped real coordinates from the wrapped scaled coordinates
    
    matVecMul3(Hmat, scaled, thePos);

  }
  else{
    // calc the scaled coordinates.
    
    for(i=0; i<3; i++)
      scaled[i] = thePos[i]*HmatInv[i][i];
    
    // wrap the scaled coordinates
    
    for(i=0; i<3; i++)
      scaled[i] -= roundMe(scaled[i]);
    
    // calc the wrapped real coordinates from the wrapped scaled coordinates
    
    for(i=0; i<3; i++)
      thePos[i] = scaled[i]*Hmat[i][i];
  }
    
}


int SimInfo::getNDF(){
  int ndf_local;

  ndf_local = 0;
  
  for(int i = 0; i < integrableObjects.size(); i++){
    ndf_local += 3;
    if (integrableObjects[i]->isDirectional()) {
      if (integrableObjects[i]->isLinear())
        ndf_local += 2;
      else
        ndf_local += 3;
    }
  }

  // n_constraints is local, so subtract them on each processor:

  ndf_local -= n_constraints;

#ifdef IS_MPI
  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
#else
  ndf = ndf_local;
#endif

  // nZconstraints is global, as are the 3 COM translations for the 
  // entire system:

  ndf = ndf - 3 - nZconstraints;

  return ndf;
}

int SimInfo::getNDFraw() {
  int ndfRaw_local;

  // Raw degrees of freedom that we have to set
  ndfRaw_local = 0;

  for(int i = 0; i < integrableObjects.size(); i++){
    ndfRaw_local += 3;
    if (integrableObjects[i]->isDirectional()) {
       if (integrableObjects[i]->isLinear())
        ndfRaw_local += 2;
      else
        ndfRaw_local += 3;
    }
  }
    
#ifdef IS_MPI
  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
#else
  ndfRaw = ndfRaw_local;
#endif

  return ndfRaw;
}

int SimInfo::getNDFtranslational() {
  int ndfTrans_local;

  ndfTrans_local = 3 * integrableObjects.size() - n_constraints;


#ifdef IS_MPI
  MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
#else
  ndfTrans = ndfTrans_local;
#endif

  ndfTrans = ndfTrans - 3 - nZconstraints;

  return ndfTrans;
}

int SimInfo::getTotIntegrableObjects() {
  int nObjs_local;
  int nObjs;

  nObjs_local =  integrableObjects.size();


#ifdef IS_MPI
  MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
#else
  nObjs = nObjs_local;
#endif


  return nObjs;
}

void SimInfo::refreshSim(){

  simtype fInfo;
  int isError;
  int n_global;
  int* excl;

  fInfo.dielect = 0.0;

  if( useDipoles ){
    if( useReactionField )fInfo.dielect = dielectric;
  }

  fInfo.SIM_uses_PBC = usePBC;
  //fInfo.SIM_uses_LJ = 0;
  fInfo.SIM_uses_LJ = useLJ;
  fInfo.SIM_uses_sticky = useSticky;
  //fInfo.SIM_uses_sticky = 0;
  fInfo.SIM_uses_charges = useCharges;
  fInfo.SIM_uses_dipoles = useDipoles;
  //fInfo.SIM_uses_dipoles = 0;
  fInfo.SIM_uses_RF = useReactionField;
  //fInfo.SIM_uses_RF = 0;
  fInfo.SIM_uses_GB = useGB;
  fInfo.SIM_uses_EAM = useEAM;

  n_exclude = excludes->getSize();
  excl = excludes->getFortranArray();
  
#ifdef IS_MPI
  n_global = mpiSim->getTotAtoms();
#else
  n_global = n_atoms;
#endif
  
  isError = 0;
  
  getFortranGroupArray(this, mfact, ngroup, groupList, groupStart);
  
  setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl, 
                  &nGlobalExcludes, globalExcludes, molMembershipArray, 
                  &mfact[0], &ngroup, &groupList[0], &groupStart[0], &isError);
  
  if( isError ){
    
    sprintf( painCave.errMsg,
             "There was an error setting the simulation information in fortran.\n" );
    painCave.isFatal = 1;
    simError();
  }
  
#ifdef IS_MPI
  sprintf( checkPointMsg,
           "succesfully sent the simulation information to fortran.\n");
  MPIcheckPoint();
#endif // is_mpi
  
  this->ndf = this->getNDF();
  this->ndfRaw = this->getNDFraw();
  this->ndfTrans = this->getNDFtranslational();
}

void SimInfo::setDefaultRcut( double theRcut ){
  
  haveRcut = 1;
  rCut = theRcut;
  
  ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
  
  notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
}

void SimInfo::setDefaultEcr( double theEcr ){

  haveEcr = 1;
  ecr = theEcr;
  
  ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;

  notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
}

void SimInfo::setDefaultEcr( double theEcr, double theEst ){

  est = theEst;
  setDefaultEcr( theEcr );
}


void SimInfo::checkCutOffs( void ){
  
  if( boxIsInit ){
    
    //we need to check cutOffs against the box
    
    if( rCut > maxCutoff ){
      sprintf( painCave.errMsg,
               "LJrcut is too large for the current periodic box.\n"
               "\tCurrent Value of LJrcut = %G at time %G\n "
               "\tThis is larger than half of at least one of the\n"
               "\tperiodic box vectors.  Right now, the Box matrix is:\n"
               "\n"
               "\t[ %G %G %G ]\n"
               "\t[ %G %G %G ]\n"
               "\t[ %G %G %G ]\n",
               rCut, currentTime,
               Hmat[0][0], Hmat[0][1], Hmat[0][2],
               Hmat[1][0], Hmat[1][1], Hmat[1][2],
               Hmat[2][0], Hmat[2][1], Hmat[2][2]);
      painCave.isFatal = 1;
      simError();
    }
    
    if( haveEcr ){
      if( ecr > maxCutoff ){
        sprintf( painCave.errMsg,
                 "electrostaticCutoffRadius is too large for the current\n"
                 "\tperiodic box.\n\n"
                 "\tCurrent Value of ECR = %G at time %G\n "
                 "\tThis is larger than half of at least one of the\n"
                 "\tperiodic box vectors.  Right now, the Box matrix is:\n"
                 "\n"
                 "\t[ %G %G %G ]\n"
                 "\t[ %G %G %G ]\n"
                 "\t[ %G %G %G ]\n",
                 ecr, currentTime,
                 Hmat[0][0], Hmat[0][1], Hmat[0][2],
                 Hmat[1][0], Hmat[1][1], Hmat[1][2],
                 Hmat[2][0], Hmat[2][1], Hmat[2][2]);
        painCave.isFatal = 1;
        simError();
      }
    }
  } else {
    // initialize this stuff before using it, OK?
    sprintf( painCave.errMsg,
             "Trying to check cutoffs without a box.\n"
             "\tOOPSE should have better programmers than that.\n" );
    painCave.isFatal = 1;
    simError();      
  }
  
}

void SimInfo::addProperty(GenericData* prop){

  map<string, GenericData*>::iterator result;
  result = properties.find(prop->getID());
  
  //we can't simply use  properties[prop->getID()] = prop,
  //it will cause memory leak if we already contain a propery which has the same name of prop
  
  if(result != properties.end()){
    
    delete (*result).second;
    (*result).second = prop;
      
  }
  else{

    properties[prop->getID()] = prop;

  }
    
}

GenericData* SimInfo::getProperty(const string& propName){
 
  map<string, GenericData*>::iterator result;
  
  //string lowerCaseName = ();
  
  result = properties.find(propName);
  
  if(result != properties.end()) 
    return (*result).second;  
  else   
    return NULL;  
}


void getFortranGroupArray(SimInfo* info, vector<double>& mfact, int& ngroup, 
                          vector<int>& groupList, vector<int>& groupStart){
  Molecule* mol;
  Atom** myAtoms;
  int numAtom;
  int curIndex;
  double mtot;

  mfact.clear();
  groupList.clear();
  groupStart.clear();
  
  //Be careful, fortran array begin at 1
  curIndex = 1;
  
  if(info->useMolecularCutoffs){
    
#ifdef IS_MPI
    ngroup = mpiSim->getMyNMol();
#else
    ngroup = info->n_mol;
#endif
    
    for(int i = 0; i < ngroup; i ++){
      mol = &(info->molecules[i]);
      numAtom = mol->getNAtoms();
      myAtoms = mol->getMyAtoms();
      mtot = 0.0;

      for(int j=0; j < numAtom; j++)
        mtot += myAtoms[j]->getMass();                
      
      for(int j=0; j < numAtom; j++){ 
               
        // We want the local Index:
        groupList.push_back(myAtoms[j]->getIndex() + 1);
        mfact.push_back(myAtoms[j]->getMass() / mtot);

      }
      
      groupStart.push_back(curIndex);
      curIndex += numAtom;
      
    } //end for(int i =0 ; i < ngroup; i++)    
  }
  else{
    //using atomic cutoff, every single atom is just a group
    
#ifdef IS_MPI
    ngroup = mpiSim->getMyNlocal();
#else
    ngroup = info->n_atoms;
#endif
    
    for(int i =0 ; i < ngroup; i++){
      groupStart.push_back(curIndex++);       
      groupList.push_back((info->atoms[i])->getIndex() + 1);
      mfact.push_back(1.0);
      
    }//end for(int i =0 ; i < ngroup; i++)
    
  }//end if (info->useMolecularCutoffs)
  
}
Revision:	1150
Committed:	Fri May 7 21:35:05 2004 UTC (20 years, 2 months ago) by gezelter
File size:	14399 byte(s)
Log Message:	Many changes to get group-based cutoffs to work
#	User	Rev	Content
1	gezelter	829	#include <stdlib.h>
2			#include <string.h>
3			#include <math.h>
4	mmeineke	377
5	mmeineke	572	#include <iostream>
6			using namespace std;
7	mmeineke	377
8			#include "SimInfo.hpp"
9			#define __C
10			#include "fSimulation.h"
11			#include "simError.h"
12
13			#include "fortranWrappers.hpp"
14
15	gezelter	1097	#include "MatVec3.h"
16
17	gezelter	490	#ifdef IS_MPI
18			#include "mpiSimulation.hpp"
19			#endif
20
21	mmeineke	572	inline double roundMe( double x ){
22			return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
23			}
24
25	mmeineke	860	inline double min( double a, double b ){
26			return (a < b ) ? a : b;
27			}
28	mmeineke	572
29	mmeineke	377	SimInfo* currentInfo;
30
31			SimInfo::SimInfo(){
32	gezelter	1097
33	mmeineke	377	n_constraints = 0;
34	tim	699	nZconstraints = 0;
35	mmeineke	377	n_oriented = 0;
36			n_dipoles = 0;
37	gezelter	458	ndf = 0;
38			ndfRaw = 0;
39	mmeineke	674	nZconstraints = 0;
40	mmeineke	377	the_integrator = NULL;
41			setTemp = 0;
42			thermalTime = 0.0;
43	mmeineke	642	currentTime = 0.0;
44	mmeineke	420	rCut = 0.0;
45	mmeineke	618	ecr = 0.0;
46	mmeineke	619	est = 0.0;
47	mmeineke	377
48	mmeineke	859	haveRcut = 0;
49			haveEcr = 0;
50	mmeineke	626	boxIsInit = 0;
51
52	tim	781	resetTime = 1e99;
53	mmeineke	626
54	gezelter	1097	orthoRhombic = 0;
55	mmeineke	855	orthoTolerance = 1E-6;
56			useInitXSstate = true;
57
58	mmeineke	377	usePBC = 0;
59			useLJ = 0;
60			useSticky = 0;
61	gezelter	941	useCharges = 0;
62			useDipoles = 0;
63	mmeineke	377	useReactionField = 0;
64			useGB = 0;
65			useEAM = 0;
66	gezelter	1139	useMolecularCutoffs = 0;
67	mmeineke	377
68	gezelter	1097	excludes = Exclude::Instance();
69
70	mmeineke	670	myConfiguration = new SimState();
71
72	tim	1031	has_minimizer = false;
73			the_minimizer =NULL;
74
75	tim	1144	ngroup = 0;
76
77	gezelter	457	wrapMeSimInfo( this );
78			}
79	mmeineke	377
80	mmeineke	670
81	tim	660	SimInfo::~SimInfo(){
82
83	mmeineke	670	delete myConfiguration;
84
85	tim	660	map<string, GenericData*>::iterator i;
86
87			for(i = properties.begin(); i != properties.end(); i++)
88			delete (*i).second;
89	tim	1144
90	tim	660	}
91
92	gezelter	457	void SimInfo::setBox(double newBox[3]) {
93	mmeineke	586
94	gezelter	588	int i, j;
95			double tempMat[3][3];
96	gezelter	463
97	gezelter	588	for(i=0; i<3; i++)
98			for (j=0; j<3; j++) tempMat[i][j] = 0.0;;
99	gezelter	463
100	gezelter	588	tempMat[0][0] = newBox[0];
101			tempMat[1][1] = newBox[1];
102			tempMat[2][2] = newBox[2];
103	gezelter	463
104	mmeineke	586	setBoxM( tempMat );
105	mmeineke	568
106	gezelter	457	}
107	mmeineke	377
108	gezelter	588	void SimInfo::setBoxM( double theBox[3][3] ){
109	mmeineke	568
110	mmeineke	787	int i, j;
111	gezelter	588	double FortranHmat[9]; // to preserve compatibility with Fortran the
112			// ordering in the array is as follows:
113			// [ 0 3 6 ]
114			// [ 1 4 7 ]
115			// [ 2 5 8 ]
116			double FortranHmatInv[9]; // the inverted Hmat (for Fortran);
117	mmeineke	568
118	mmeineke	626	if( !boxIsInit ) boxIsInit = 1;
119	mmeineke	586
120	gezelter	588	for(i=0; i < 3; i++)
121			for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
122
123	mmeineke	568	calcBoxL();
124	gezelter	588	calcHmatInv();
125	mmeineke	568
126	gezelter	588	for(i=0; i < 3; i++) {
127			for (j=0; j < 3; j++) {
128			FortranHmat[3*j + i] = Hmat[i][j];
129			FortranHmatInv[3*j + i] = HmatInv[i][j];
130			}
131			}
132	mmeineke	586
133	mmeineke	590	setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic);
134	mmeineke	568
135	mmeineke	377	}
136	gezelter	458
137	mmeineke	568
138	gezelter	588	void SimInfo::getBoxM (double theBox[3][3]) {
139	mmeineke	568
140	gezelter	588	int i, j;
141			for(i=0; i<3; i++)
142			for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j];
143	mmeineke	568	}
144
145	gezelter	574
146			void SimInfo::scaleBox(double scale) {
147	gezelter	588	double theBox[3][3];
148			int i, j;
149	gezelter	574
150	gezelter	617	// cerr << "Scaling box by " << scale << "\n";
151	mmeineke	586
152	gezelter	588	for(i=0; i<3; i++)
153			for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;
154	gezelter	574
155			setBoxM(theBox);
156
157			}
158
159	gezelter	588	void SimInfo::calcHmatInv( void ) {
160	mmeineke	590
161	mmeineke	853	int oldOrtho;
162	mmeineke	590	int i,j;
163	mmeineke	569	double smallDiag;
164			double tol;
165			double sanity[3][3];
166	mmeineke	568
167	gezelter	588	invertMat3( Hmat, HmatInv );
168	mmeineke	568
169	gezelter	588	// check to see if Hmat is orthorhombic
170	mmeineke	568
171	mmeineke	853	oldOrtho = orthoRhombic;
172
173			smallDiag = fabs(Hmat[0][0]);
174			if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]);
175			if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]);
176	mmeineke	855	tol = smallDiag * orthoTolerance;
177	mmeineke	568
178	gezelter	588	orthoRhombic = 1;
179	mmeineke	568
180	gezelter	588	for (i = 0; i < 3; i++ ) {
181			for (j = 0 ; j < 3; j++) {
182			if (i != j) {
183			if (orthoRhombic) {
184	mmeineke	853	if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0;
185	gezelter	588	}
186			}
187	mmeineke	568	}
188			}
189	mmeineke	853
190			if( oldOrtho != orthoRhombic ){
191
192			if( orthoRhombic ){
193			sprintf( painCave.errMsg,
194	gezelter	1097	"OOPSE is switching from the default Non-Orthorhombic\n"
195			"\tto the faster Orthorhombic periodic boundary computations.\n"
196			"\tThis is usually a good thing, but if you wan't the\n"
197			"\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
198			"\tvariable ( currently set to %G ) smaller.\n",
199	mmeineke	855	orthoTolerance);
200	mmeineke	853	simError();
201			}
202			else {
203			sprintf( painCave.errMsg,
204	gezelter	1097	"OOPSE is switching from the faster Orthorhombic to the more\n"
205			"\tflexible Non-Orthorhombic periodic boundary computations.\n"
206			"\tThis is usually because the box has deformed under\n"
207			"\tNPTf integration. If you wan't to live on the edge with\n"
208			"\tthe Orthorhombic computations, make the orthoBoxTolerance\n"
209			"\tvariable ( currently set to %G ) larger.\n",
210	mmeineke	855	orthoTolerance);
211	mmeineke	853	simError();
212			}
213			}
214	gezelter	588	}
215	mmeineke	569
216	mmeineke	568	void SimInfo::calcBoxL( void ){
217
218			double dx, dy, dz, dsq;
219
220	gezelter	588	// boxVol = Determinant of Hmat
221	mmeineke	568
222	gezelter	588	boxVol = matDet3( Hmat );
223	mmeineke	568
224			// boxLx
225
226	gezelter	588	dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
227	mmeineke	568	dsq = dxdx + dydy + dz*dz;
228	gezelter	621	boxL[0] = sqrt( dsq );
229	tim	781	//maxCutoff = 0.5 * boxL[0];
230	mmeineke	568
231			// boxLy
232
233	gezelter	588	dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
234	mmeineke	568	dsq = dxdx + dydy + dz*dz;
235	gezelter	621	boxL[1] = sqrt( dsq );
236	tim	781	//if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];
237	mmeineke	568
238	tim	781
239	mmeineke	568	// boxLz
240
241	gezelter	588	dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
242	mmeineke	568	dsq = dxdx + dydy + dz*dz;
243	gezelter	621	boxL[2] = sqrt( dsq );
244	tim	781	//if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];
245
246			//calculate the max cutoff
247			maxCutoff = calcMaxCutOff();
248	chuckv	669
249			checkCutOffs();
250	mmeineke	626
251	mmeineke	568	}
252
253
254	tim	781	double SimInfo::calcMaxCutOff(){
255
256			double ri[3], rj[3], rk[3];
257			double rij[3], rjk[3], rki[3];
258			double minDist;
259
260			ri[0] = Hmat[0][0];
261			ri[1] = Hmat[1][0];
262			ri[2] = Hmat[2][0];
263
264			rj[0] = Hmat[0][1];
265			rj[1] = Hmat[1][1];
266			rj[2] = Hmat[2][1];
267
268			rk[0] = Hmat[0][2];
269			rk[1] = Hmat[1][2];
270			rk[2] = Hmat[2][2];
271	gezelter	1097
272			crossProduct3(ri, rj, rij);
273			distXY = dotProduct3(rk,rij) / norm3(rij);
274	tim	781
275			crossProduct3(rj,rk, rjk);
276	gezelter	1097	distYZ = dotProduct3(ri,rjk) / norm3(rjk);
277	tim	781
278			crossProduct3(rk,ri, rki);
279	gezelter	1097	distZX = dotProduct3(rj,rki) / norm3(rki);
280	tim	781
281			minDist = min(min(distXY, distYZ), distZX);
282			return minDist/2;
283
284			}
285
286	mmeineke	568	void SimInfo::wrapVector( double thePos[3] ){
287
288	mmeineke	787	int i;
289	mmeineke	568	double scaled[3];
290
291	mmeineke	569	if( !orthoRhombic ){
292			// calc the scaled coordinates.
293	gezelter	588
294
295			matVecMul3(HmatInv, thePos, scaled);
296	mmeineke	569
297			for(i=0; i<3; i++)
298	mmeineke	572	scaled[i] -= roundMe(scaled[i]);
299	mmeineke	569
300			// calc the wrapped real coordinates from the wrapped scaled coordinates
301
302	gezelter	588	matVecMul3(Hmat, scaled, thePos);
303
304	mmeineke	569	}
305			else{
306			// calc the scaled coordinates.
307
308			for(i=0; i<3; i++)
309	gezelter	588	scaled[i] = thePos[i]*HmatInv[i][i];
310	mmeineke	569
311			// wrap the scaled coordinates
312
313			for(i=0; i<3; i++)
314	mmeineke	572	scaled[i] -= roundMe(scaled[i]);
315	mmeineke	569
316			// calc the wrapped real coordinates from the wrapped scaled coordinates
317
318			for(i=0; i<3; i++)
319	gezelter	588	thePos[i] = scaled[i]*Hmat[i][i];
320	mmeineke	569	}
321
322	mmeineke	568	}
323
324
325	gezelter	458	int SimInfo::getNDF(){
326	mmeineke	790	int ndf_local;
327	gezelter	458
328	tim	1113	ndf_local = 0;
329
330	gezelter	1097	for(int i = 0; i < integrableObjects.size(); i++){
331			ndf_local += 3;
332	tim	1118	if (integrableObjects[i]->isDirectional()) {
333			if (integrableObjects[i]->isLinear())
334			ndf_local += 2;
335			else
336			ndf_local += 3;
337			}
338	gezelter	1097	}
339
340			// n_constraints is local, so subtract them on each processor:
341
342			ndf_local -= n_constraints;
343
344	gezelter	458	#ifdef IS_MPI
345			MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
346			#else
347			ndf = ndf_local;
348			#endif
349
350	gezelter	1097	// nZconstraints is global, as are the 3 COM translations for the
351			// entire system:
352
353	mmeineke	674	ndf = ndf - 3 - nZconstraints;
354	gezelter	458
355			return ndf;
356			}
357
358			int SimInfo::getNDFraw() {
359	mmeineke	790	int ndfRaw_local;
360	gezelter	458
361			// Raw degrees of freedom that we have to set
362	tim	1113	ndfRaw_local = 0;
363	gezelter	1097
364			for(int i = 0; i < integrableObjects.size(); i++){
365			ndfRaw_local += 3;
366	tim	1118	if (integrableObjects[i]->isDirectional()) {
367			if (integrableObjects[i]->isLinear())
368			ndfRaw_local += 2;
369			else
370			ndfRaw_local += 3;
371			}
372	gezelter	1097	}
373
374	gezelter	458	#ifdef IS_MPI
375			MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
376			#else
377			ndfRaw = ndfRaw_local;
378			#endif
379
380			return ndfRaw;
381			}
382	tim	767
383			int SimInfo::getNDFtranslational() {
384	mmeineke	790	int ndfTrans_local;
385	tim	767
386	gezelter	1097	ndfTrans_local = 3 * integrableObjects.size() - n_constraints;
387	tim	767
388	gezelter	1097
389	tim	767	#ifdef IS_MPI
390			MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
391			#else
392			ndfTrans = ndfTrans_local;
393			#endif
394
395			ndfTrans = ndfTrans - 3 - nZconstraints;
396
397			return ndfTrans;
398			}
399
400	tim	1108	int SimInfo::getTotIntegrableObjects() {
401			int nObjs_local;
402			int nObjs;
403
404			nObjs_local = integrableObjects.size();
405
406
407			#ifdef IS_MPI
408			MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
409			#else
410			nObjs = nObjs_local;
411			#endif
412
413
414			return nObjs;
415			}
416
417	mmeineke	377	void SimInfo::refreshSim(){
418
419			simtype fInfo;
420			int isError;
421	gezelter	490	int n_global;
422	mmeineke	424	int* excl;
423	mmeineke	626
424	mmeineke	469	fInfo.dielect = 0.0;
425	mmeineke	377
426	gezelter	941	if( useDipoles ){
427	mmeineke	469	if( useReactionField )fInfo.dielect = dielectric;
428			}
429
430	mmeineke	377	fInfo.SIM_uses_PBC = usePBC;
431	mmeineke	443	//fInfo.SIM_uses_LJ = 0;
432	chuckv	439	fInfo.SIM_uses_LJ = useLJ;
433	mmeineke	443	fInfo.SIM_uses_sticky = useSticky;
434			//fInfo.SIM_uses_sticky = 0;
435	gezelter	941	fInfo.SIM_uses_charges = useCharges;
436			fInfo.SIM_uses_dipoles = useDipoles;
437	chuckv	482	//fInfo.SIM_uses_dipoles = 0;
438	chrisfen	999	fInfo.SIM_uses_RF = useReactionField;
439			//fInfo.SIM_uses_RF = 0;
440	mmeineke	377	fInfo.SIM_uses_GB = useGB;
441			fInfo.SIM_uses_EAM = useEAM;
442
443	gezelter	1097	n_exclude = excludes->getSize();
444			excl = excludes->getFortranArray();
445	tim	1144
446	gezelter	490	#ifdef IS_MPI
447			n_global = mpiSim->getTotAtoms();
448			#else
449			n_global = n_atoms;
450			#endif
451	gezelter	1150
452	mmeineke	377	isError = 0;
453	gezelter	1150
454			getFortranGroupArray(this, mfact, ngroup, groupList, groupStart);
455
456	gezelter	490	setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
457	gezelter	1150	&nGlobalExcludes, globalExcludes, molMembershipArray,
458			&mfact[0], &ngroup, &groupList[0], &groupStart[0], &isError);
459
460	mmeineke	377	if( isError ){
461	gezelter	1150
462	mmeineke	377	sprintf( painCave.errMsg,
463	gezelter	1150	"There was an error setting the simulation information in fortran.\n" );
464	mmeineke	377	painCave.isFatal = 1;
465			simError();
466			}
467	gezelter	1150
468	mmeineke	377	#ifdef IS_MPI
469			sprintf( checkPointMsg,
470			"succesfully sent the simulation information to fortran.\n");
471			MPIcheckPoint();
472			#endif // is_mpi
473	gezelter	1150
474	gezelter	474	this->ndf = this->getNDF();
475			this->ndfRaw = this->getNDFraw();
476	tim	767	this->ndfTrans = this->getNDFtranslational();
477	mmeineke	377	}
478
479	mmeineke	841	void SimInfo::setDefaultRcut( double theRcut ){
480	gezelter	1150
481	mmeineke	859	haveRcut = 1;
482	mmeineke	841	rCut = theRcut;
483	gezelter	1150
484	gezelter	845	( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
485	gezelter	1150
486	mmeineke	843	notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
487	mmeineke	841	}
488
489			void SimInfo::setDefaultEcr( double theEcr ){
490
491	mmeineke	859	haveEcr = 1;
492	chrisfen	872	ecr = theEcr;
493	mmeineke	841
494	gezelter	845	( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
495
496	mmeineke	843	notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
497	mmeineke	841	}
498
499			void SimInfo::setDefaultEcr( double theEcr, double theEst ){
500	mmeineke	626
501	mmeineke	841	est = theEst;
502			setDefaultEcr( theEcr );
503			}
504
505
506	mmeineke	626	void SimInfo::checkCutOffs( void ){
507	gezelter	770
508	mmeineke	626	if( boxIsInit ){
509
510			//we need to check cutOffs against the box
511	mmeineke	859
512			if( rCut > maxCutoff ){
513	mmeineke	626	sprintf( painCave.errMsg,
514	gezelter	1097	"LJrcut is too large for the current periodic box.\n"
515			"\tCurrent Value of LJrcut = %G at time %G\n "
516			"\tThis is larger than half of at least one of the\n"
517			"\tperiodic box vectors. Right now, the Box matrix is:\n"
518	tim	1131	"\n"
519	gezelter	965	"\t[ %G %G %G ]\n"
520			"\t[ %G %G %G ]\n"
521			"\t[ %G %G %G ]\n",
522	tim	1131	rCut, currentTime,
523	mmeineke	874	Hmat[0][0], Hmat[0][1], Hmat[0][2],
524			Hmat[1][0], Hmat[1][1], Hmat[1][2],
525			Hmat[2][0], Hmat[2][1], Hmat[2][2]);
526	mmeineke	859	painCave.isFatal = 1;
527	mmeineke	626	simError();
528			}
529	mmeineke	859
530			if( haveEcr ){
531			if( ecr > maxCutoff ){
532			sprintf( painCave.errMsg,
533	gezelter	1097	"electrostaticCutoffRadius is too large for the current\n"
534			"\tperiodic box.\n\n"
535			"\tCurrent Value of ECR = %G at time %G\n "
536			"\tThis is larger than half of at least one of the\n"
537			"\tperiodic box vectors. Right now, the Box matrix is:\n"
538			"\n"
539			"\t[ %G %G %G ]\n"
540			"\t[ %G %G %G ]\n"
541			"\t[ %G %G %G ]\n",
542	mmeineke	874	ecr, currentTime,
543			Hmat[0][0], Hmat[0][1], Hmat[0][2],
544			Hmat[1][0], Hmat[1][1], Hmat[1][2],
545			Hmat[2][0], Hmat[2][1], Hmat[2][2]);
546	mmeineke	859	painCave.isFatal = 1;
547			simError();
548	tim	781	}
549			}
550	tim	767	} else {
551			// initialize this stuff before using it, OK?
552	gezelter	770	sprintf( painCave.errMsg,
553	gezelter	965	"Trying to check cutoffs without a box.\n"
554			"\tOOPSE should have better programmers than that.\n" );
555	gezelter	770	painCave.isFatal = 1;
556			simError();
557	mmeineke	626	}
558	gezelter	770
559	mmeineke	626	}
560	tim	658
561			void SimInfo::addProperty(GenericData* prop){
562
563			map<string, GenericData*>::iterator result;
564			result = properties.find(prop->getID());
565
566			//we can't simply use properties[prop->getID()] = prop,
567			//it will cause memory leak if we already contain a propery which has the same name of prop
568
569			if(result != properties.end()){
570
571			delete (*result).second;
572			(*result).second = prop;
573
574			}
575			else{
576
577			properties[prop->getID()] = prop;
578
579			}
580
581			}
582
583			GenericData* SimInfo::getProperty(const string& propName){
584
585			map<string, GenericData*>::iterator result;
586
587			//string lowerCaseName = ();
588
589			result = properties.find(propName);
590
591			if(result != properties.end())
592			return (*result).second;
593			else
594			return NULL;
595			}
596
597	tim	1144
598			void getFortranGroupArray(SimInfo* info, vector<double>& mfact, int& ngroup,
599	gezelter	1150	vector<int>& groupList, vector<int>& groupStart){
600	tim	1144	Molecule* mol;
601	gezelter	1150	Atom** myAtoms;
602	tim	1144	int numAtom;
603			int curIndex;
604	gezelter	1150	double mtot;
605	tim	1144
606			mfact.clear();
607			groupList.clear();
608			groupStart.clear();
609	gezelter	1150
610	tim	1144	//Be careful, fortran array begin at 1
611			curIndex = 1;
612	gezelter	1150
613			if(info->useMolecularCutoffs){
614	tim	1144
615	gezelter	1150	#ifdef IS_MPI
616			ngroup = mpiSim->getMyNMol();
617			#else
618	tim	1144	ngroup = info->n_mol;
619	gezelter	1150	#endif
620
621	tim	1144	for(int i = 0; i < ngroup; i ++){
622			mol = &(info->molecules[i]);
623			numAtom = mol->getNAtoms();
624	gezelter	1150	myAtoms = mol->getMyAtoms();
625			mtot = 0.0;
626
627			for(int j=0; j < numAtom; j++)
628			mtot += myAtoms[j]->getMass();
629	tim	1144
630	gezelter	1150	for(int j=0; j < numAtom; j++){
631
632			// We want the local Index:
633			groupList.push_back(myAtoms[j]->getIndex() + 1);
634			mfact.push_back(myAtoms[j]->getMass() / mtot);
635
636			}
637
638	tim	1144	groupStart.push_back(curIndex);
639			curIndex += numAtom;
640
641	gezelter	1150	} //end for(int i =0 ; i < ngroup; i++)
642	tim	1144	}
643			else{
644			//using atomic cutoff, every single atom is just a group
645	gezelter	1150
646			#ifdef IS_MPI
647			ngroup = mpiSim->getMyNlocal();
648			#else
649	tim	1144	ngroup = info->n_atoms;
650	gezelter	1150	#endif
651
652	tim	1144	for(int i =0 ; i < ngroup; i++){
653	gezelter	1150	groupStart.push_back(curIndex++);
654	tim	1144	groupList.push_back((info->atoms[i])->getIndex() + 1);
655	gezelter	1150	mfact.push_back(1.0);
656
657	tim	1144	}//end for(int i =0 ; i < ngroup; i++)
658	gezelter	1150
659	tim	1144	}//end if (info->useMolecularCutoffs)
660	gezelter	1150
661	tim	1144	}