src/brains/SimInfo.cpp

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

#include <iostream>
using namespace std;

#include "brains/SimInfo.hpp"
#define __C
#include "brains/fSimulation.h"
#include "utils/simError.h"
#include "UseTheForce/DarkSide/simulation_interface.h"
#include "UseTheForce/notifyCutoffs_interface.h"

//#include "UseTheForce/fortranWrappers.hpp"

#include "math/MatVec3.h"

#ifdef IS_MPI
#include "brains/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;
  useDirectionalAtoms = 0;
  useLennardJones = 0; 
  useElectrostatics = 0;
  useCharges = 0;
  useDipoles = 0;
  useSticky = 0;
  useGayBerne = 0;
  useEAM = 0;
  useShapes = 0;
  useFLARB = 0;

  useSolidThermInt = 0;
  useLiquidThermInt = 0;

  haveCutoffGroups = false;

  excludes = Exclude::Instance();

  myConfiguration = new SimState();

  has_minimizer = false;
  the_minimizer =NULL;

  ngroup = 0;

}


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];
    }
  }

  setFortranBox(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);
      painCave.severity = OOPSE_INFO;
      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);
      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;

  if (useSticky || useDipoles || useGayBerne || useShapes) {
    useDirectionalAtoms = 1;
    fInfo.SIM_uses_DirectionalAtoms = useDirectionalAtoms;
  }

  fInfo.SIM_uses_LennardJones = useLennardJones;

  if (useCharges || useDipoles) {
    useElectrostatics = 1;
    fInfo.SIM_uses_Electrostatics = useElectrostatics;
  }

  fInfo.SIM_uses_Charges = useCharges;
  fInfo.SIM_uses_Dipoles = useDipoles;
  fInfo.SIM_uses_Sticky = useSticky;
  fInfo.SIM_uses_GayBerne = useGayBerne;
  fInfo.SIM_uses_EAM = useEAM;
  fInfo.SIM_uses_Shapes = useShapes;
  fInfo.SIM_uses_FLARB = useFLARB;
  fInfo.SIM_uses_RF = useReactionField;

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