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;

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