OOPSE/libmdtools/SimInfo.cpp

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

#include <iostream>
using namespace std;

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

#include "fortranWrappers.hpp"

#include "MatVec3.h"

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

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

SimInfo* currentInfo;

SimInfo::SimInfo(){

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

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

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

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

  excludes = Exclude::Instance();

  myConfiguration = new SimState();

  has_minimizer = false;
  the_minimizer =NULL;

  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;
  fInfo.SIM_uses_molecular_cutoffs = useMolecularCutoffs;

  n_exclude = excludes->getSize();
  excl = excludes->getFortranArray();

#ifdef IS_MPI
  n_global = mpiSim->getTotAtoms();
#else
  n_global = n_atoms;
#endif

  isError = 0;

  setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl, 
                  &nGlobalExcludes, globalExcludes, molMembershipArray, 
                  &isError );

  if( isError ){

    sprintf( painCave.errMsg,
             "There was an error setting the simulation information in fortran.\n" );
    painCave.isFatal = 1;
    simError();
  }

#ifdef IS_MPI
  sprintf( checkPointMsg,
           "succesfully sent the simulation information to fortran.\n");
  MPIcheckPoint();
#endif // is_mpi

  this->ndf = this->getNDF();
  this->ndfRaw = this->getNDFraw();
  this->ndfTrans = this->getNDFtranslational();
}

void SimInfo::setDefaultRcut( double theRcut ){

  haveRcut = 1;
  rCut = theRcut;

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

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

void SimInfo::setDefaultEcr( double theEcr ){

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

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

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

  est = theEst;
  setDefaultEcr( theEcr );
}


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

void SimInfo::addProperty(GenericData* prop){

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

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

  }
    
}

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

Revision:	1139
Committed:	Wed Apr 28 22:06:29 2004 UTC (20 years, 2 months ago) by gezelter
File size:	12744 byte(s)
Log Message:	Adding molecular cutoffs
#	User	Rev	Content
1	gezelter	829	#include <stdlib.h>
2			#include <string.h>
3			#include <math.h>
4	mmeineke	377
5	mmeineke	572	#include <iostream>
6			using namespace std;
7	mmeineke	377
8			#include "SimInfo.hpp"
9			#define __C
10			#include "fSimulation.h"
11			#include "simError.h"
12
13			#include "fortranWrappers.hpp"
14
15	gezelter	1097	#include "MatVec3.h"
16
17	gezelter	490	#ifdef IS_MPI
18			#include "mpiSimulation.hpp"
19			#endif
20
21	mmeineke	572	inline double roundMe( double x ){
22			return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
23			}
24
25	mmeineke	860	inline double min( double a, double b ){
26			return (a < b ) ? a : b;
27			}
28	mmeineke	572
29	mmeineke	377	SimInfo* currentInfo;
30
31			SimInfo::SimInfo(){
32	gezelter	1097
33	mmeineke	377	n_constraints = 0;
34	tim	699	nZconstraints = 0;
35	mmeineke	377	n_oriented = 0;
36			n_dipoles = 0;
37	gezelter	458	ndf = 0;
38			ndfRaw = 0;
39	mmeineke	674	nZconstraints = 0;
40	mmeineke	377	the_integrator = NULL;
41			setTemp = 0;
42			thermalTime = 0.0;
43	mmeineke	642	currentTime = 0.0;
44	mmeineke	420	rCut = 0.0;
45	mmeineke	618	ecr = 0.0;
46	mmeineke	619	est = 0.0;
47	mmeineke	377
48	mmeineke	859	haveRcut = 0;
49			haveEcr = 0;
50	mmeineke	626	boxIsInit = 0;
51
52	tim	781	resetTime = 1e99;
53	mmeineke	626
54	gezelter	1097	orthoRhombic = 0;
55	mmeineke	855	orthoTolerance = 1E-6;
56			useInitXSstate = true;
57
58	mmeineke	377	usePBC = 0;
59			useLJ = 0;
60			useSticky = 0;
61	gezelter	941	useCharges = 0;
62			useDipoles = 0;
63	mmeineke	377	useReactionField = 0;
64			useGB = 0;
65			useEAM = 0;
66	gezelter	1139	useMolecularCutoffs = 0;
67	mmeineke	377
68	gezelter	1097	excludes = Exclude::Instance();
69
70	mmeineke	670	myConfiguration = new SimState();
71
72	tim	1031	has_minimizer = false;
73			the_minimizer =NULL;
74
75	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	mmeineke	670
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	gezelter	1139	fInfo.SIM_uses_molecular_cutoffs = useMolecularCutoffs;
441	mmeineke	377
442	gezelter	1097	n_exclude = excludes->getSize();
443			excl = excludes->getFortranArray();
444	mmeineke	377
445	gezelter	490	#ifdef IS_MPI
446			n_global = mpiSim->getTotAtoms();
447			#else
448			n_global = n_atoms;
449			#endif
450
451	mmeineke	377	isError = 0;
452
453	gezelter	490	setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
454	gezelter	483	&nGlobalExcludes, globalExcludes, molMembershipArray,
455			&isError );
456	mmeineke	377
457			if( isError ){
458
459			sprintf( painCave.errMsg,
460			"There was an error setting the simulation information in fortran.\n" );
461			painCave.isFatal = 1;
462			simError();
463			}
464
465			#ifdef IS_MPI
466			sprintf( checkPointMsg,
467			"succesfully sent the simulation information to fortran.\n");
468			MPIcheckPoint();
469			#endif // is_mpi
470	gezelter	458
471	gezelter	474	this->ndf = this->getNDF();
472			this->ndfRaw = this->getNDFraw();
473	tim	767	this->ndfTrans = this->getNDFtranslational();
474	mmeineke	377	}
475
476	mmeineke	841	void SimInfo::setDefaultRcut( double theRcut ){
477
478	mmeineke	859	haveRcut = 1;
479	mmeineke	841	rCut = theRcut;
480	mmeineke	843
481	gezelter	845	( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
482
483	mmeineke	843	notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
484	mmeineke	841	}
485
486			void SimInfo::setDefaultEcr( double theEcr ){
487
488	mmeineke	859	haveEcr = 1;
489	chrisfen	872	ecr = theEcr;
490	mmeineke	841
491	gezelter	845	( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
492
493	mmeineke	843	notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
494	mmeineke	841	}
495
496			void SimInfo::setDefaultEcr( double theEcr, double theEst ){
497	mmeineke	626
498	mmeineke	841	est = theEst;
499			setDefaultEcr( theEcr );
500			}
501
502
503	mmeineke	626	void SimInfo::checkCutOffs( void ){
504	gezelter	770
505	mmeineke	626	if( boxIsInit ){
506
507			//we need to check cutOffs against the box
508	mmeineke	859
509			if( rCut > maxCutoff ){
510	mmeineke	626	sprintf( painCave.errMsg,
511	gezelter	1097	"LJrcut is too large for the current periodic box.\n"
512			"\tCurrent Value of LJrcut = %G at time %G\n "
513			"\tThis is larger than half of at least one of the\n"
514			"\tperiodic box vectors. Right now, the Box matrix is:\n"
515	tim	1131	"\n"
516	gezelter	965	"\t[ %G %G %G ]\n"
517			"\t[ %G %G %G ]\n"
518			"\t[ %G %G %G ]\n",
519	tim	1131	rCut, currentTime,
520	mmeineke	874	Hmat[0][0], Hmat[0][1], Hmat[0][2],
521			Hmat[1][0], Hmat[1][1], Hmat[1][2],
522			Hmat[2][0], Hmat[2][1], Hmat[2][2]);
523	mmeineke	859	painCave.isFatal = 1;
524	mmeineke	626	simError();
525			}
526	mmeineke	859
527			if( haveEcr ){
528			if( ecr > maxCutoff ){
529			sprintf( painCave.errMsg,
530	gezelter	1097	"electrostaticCutoffRadius is too large for the current\n"
531			"\tperiodic box.\n\n"
532			"\tCurrent Value of ECR = %G at time %G\n "
533			"\tThis is larger than half of at least one of the\n"
534			"\tperiodic box vectors. Right now, the Box matrix is:\n"
535			"\n"
536			"\t[ %G %G %G ]\n"
537			"\t[ %G %G %G ]\n"
538			"\t[ %G %G %G ]\n",
539	mmeineke	874	ecr, currentTime,
540			Hmat[0][0], Hmat[0][1], Hmat[0][2],
541			Hmat[1][0], Hmat[1][1], Hmat[1][2],
542			Hmat[2][0], Hmat[2][1], Hmat[2][2]);
543	mmeineke	859	painCave.isFatal = 1;
544			simError();
545	tim	781	}
546			}
547	tim	767	} else {
548			// initialize this stuff before using it, OK?
549	gezelter	770	sprintf( painCave.errMsg,
550	gezelter	965	"Trying to check cutoffs without a box.\n"
551			"\tOOPSE should have better programmers than that.\n" );
552	gezelter	770	painCave.isFatal = 1;
553			simError();
554	mmeineke	626	}
555	gezelter	770
556	mmeineke	626	}
557	tim	658
558			void SimInfo::addProperty(GenericData* prop){
559
560			map<string, GenericData*>::iterator result;
561			result = properties.find(prop->getID());
562
563			//we can't simply use properties[prop->getID()] = prop,
564			//it will cause memory leak if we already contain a propery which has the same name of prop
565
566			if(result != properties.end()){
567
568			delete (*result).second;
569			(*result).second = prop;
570
571			}
572			else{
573
574			properties[prop->getID()] = prop;
575
576			}
577
578			}
579
580			GenericData* SimInfo::getProperty(const string& propName){
581
582			map<string, GenericData*>::iterator result;
583
584			//string lowerCaseName = ();
585
586			result = properties.find(propName);
587
588			if(result != properties.end())
589			return (*result).second;
590			else
591			return NULL;
592			}
593