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;
  rSw = 0.0;

  haveRcut = 0;
  haveRsw = 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;
  useSolidThermInt = 0;
  useLiquidThermInt = 0;

  haveCutoffGroups = false;

  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,
               "\n\tOOPSE 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);
      painCave.severity = OOPSE_INFO;
      simError();
    }
    else {
      sprintf( painCave.errMsg,
               "\n\tOOPSE 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);
      painCave.severity = OOPSE_WARNING;
      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->getNAtomsGlobal();
#else
  n_global = n_atoms;
#endif
  
  isError = 0;
  
  getFortranGroupArrays(this, FglobalGroupMembership, mfact);
  //it may not be a good idea to pass the address of first element in vector
  //since c++ standard does not require vector to be stored continuously in meomory
  //Most of the compilers will organize the memory of vector continuously
  setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl, 
                  &nGlobalExcludes, globalExcludes, molMembershipArray, 
                  &mfact[0], &ngroup, &FglobalGroupMembership[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;
  rList = rCut + 1.0; 
  
  notifyFortranCutOffs( &rCut, &rSw, &rList );
}

void SimInfo::setDefaultRcut( double theRcut, double theRsw ){

  rSw = theRsw;
  setDefaultRcut( theRcut );
}


void SimInfo::checkCutOffs( void ){
  
  if( boxIsInit ){
    
    //we need to check cutOffs against the box
    
    if( rCut > maxCutoff ){
      sprintf( painCave.errMsg,
               "\n\tcutoffRadius is too large for the current periodic box.\n"
               "\tCurrent Value of cutoffRadius = %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.severity = OOPSE_ERROR;
      painCave.isFatal = 1;
      simError();
    }    
  } else {
    // initialize this stuff before using it, OK?
    sprintf( painCave.errMsg,
             "\n\tTrying to check cutoffs without a box.\n"
             "\tOOPSE should have better programmers than that.\n" );
    painCave.severity = OOPSE_ERROR;
    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 SimInfo::getFortranGroupArrays(SimInfo* info, 
                                    vector<int>& FglobalGroupMembership,
                                    vector<double>& mfact){
  
  Molecule* myMols;
  Atom** myAtoms;
  int numAtom;
  double mtot;
  int numMol;
  int numCutoffGroups;
  CutoffGroup* myCutoffGroup;
  vector<CutoffGroup*>::iterator iterCutoff;
  Atom* cutoffAtom;
  vector<Atom*>::iterator iterAtom;
  int atomIndex;
  double totalMass;
  
  mfact.clear();
  FglobalGroupMembership.clear();
  

  // Fix the silly fortran indexing problem
#ifdef IS_MPI
  numAtom = mpiSim->getNAtomsGlobal();
#else
  numAtom = n_atoms;
#endif
  for (int i = 0; i < numAtom; i++) 
    FglobalGroupMembership.push_back(globalGroupMembership[i] + 1);
  

  myMols = info->molecules;
  numMol = info->n_mol;
  for(int i  = 0; i < numMol; i++){
    numCutoffGroups = myMols[i].getNCutoffGroups();
    for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff); 
        myCutoffGroup != NULL; 
        myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){

      totalMass = myCutoffGroup->getMass();
      
      for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom); 
          cutoffAtom != NULL; 
          cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){
        mfact.push_back(cutoffAtom->getMass()/totalMass);
      }  
    }
  }

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