OOPSE/libmdtools/SimInfo.cpp

#include <cstdlib>
#include <cstring>
#include <cmath>

#include <iostream>
using namespace std;

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

#include "fortranWrappers.hpp"

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

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

SimInfo* currentInfo;

SimInfo::SimInfo(){
  excludes = NULL;
  n_constraints = 0;
  n_oriented = 0;
  n_dipoles = 0;
  ndf = 0;
  ndfRaw = 0;
  the_integrator = NULL;
  setTemp = 0;
  thermalTime = 0.0;
  rCut = 0.0;
  ecr = 0.0;
  est = 0.0;
  oldEcr = 0.0;
  oldRcut = 0.0;

  haveOrigRcut = 0;
  haveOrigEcr = 0;
  boxIsInit = 0;
  
  
  usePBC = 0;
  useLJ = 0; 
  useSticky = 0;
  useDipole = 0;
  useReactionField = 0;
  useGB = 0;
  useEAM = 0;

  wrapMeSimInfo( this );
}

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, status;
  double smallestBoxL, maxCutoff;
  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 i,j;
  double smallDiag;
  double tol;
  double sanity[3][3];

  invertMat3( Hmat, HmatInv );

  // Check the inverse to make sure it is sane:

  matMul3( Hmat, HmatInv, sanity );
    
  // check to see if Hmat is orthorhombic
  
  smallDiag = Hmat[0][0];
  if(smallDiag > Hmat[1][1]) smallDiag = Hmat[1][1];
  if(smallDiag > Hmat[2][2]) smallDiag = Hmat[2][2];
  tol = smallDiag * 1E-6;

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

double SimInfo::matDet3(double a[3][3]) {
  int i, j, k;
  double determinant;

  determinant = 0.0;

  for(i = 0; i < 3; i++) {
    j = (i+1)%3;
    k = (i+2)%3;

    determinant += a[0][i] * (a[1][j]*a[2][k] - a[1][k]*a[2][j]);
  }

  return determinant;
}

void SimInfo::invertMat3(double a[3][3], double b[3][3]) {
  
  int  i, j, k, l, m, n;
  double determinant;

  determinant = matDet3( a );

  if (determinant == 0.0) {
    sprintf( painCave.errMsg,
             "Can't invert a matrix with a zero determinant!\n");
    painCave.isFatal = 1;
    simError();
  }

  for (i=0; i < 3; i++) {
    j = (i+1)%3;
    k = (i+2)%3;
    for(l = 0; l < 3; l++) {
      m = (l+1)%3;
      n = (l+2)%3;
      
      b[l][i] = (a[j][m]*a[k][n] - a[j][n]*a[k][m]) / determinant;
    }
  }
}

void SimInfo::matMul3(double a[3][3], double b[3][3], double c[3][3]) {
  double r00, r01, r02, r10, r11, r12, r20, r21, r22;

  r00 = a[0][0]*b[0][0] + a[0][1]*b[1][0] + a[0][2]*b[2][0];
  r01 = a[0][0]*b[0][1] + a[0][1]*b[1][1] + a[0][2]*b[2][1];
  r02 = a[0][0]*b[0][2] + a[0][1]*b[1][2] + a[0][2]*b[2][2];
  
  r10 = a[1][0]*b[0][0] + a[1][1]*b[1][0] + a[1][2]*b[2][0];
  r11 = a[1][0]*b[0][1] + a[1][1]*b[1][1] + a[1][2]*b[2][1];
  r12 = a[1][0]*b[0][2] + a[1][1]*b[1][2] + a[1][2]*b[2][2];
  
  r20 = a[2][0]*b[0][0] + a[2][1]*b[1][0] + a[2][2]*b[2][0];
  r21 = a[2][0]*b[0][1] + a[2][1]*b[1][1] + a[2][2]*b[2][1];
  r22 = a[2][0]*b[0][2] + a[2][1]*b[1][2] + a[2][2]*b[2][2];
  
  c[0][0] = r00; c[0][1] = r01; c[0][2] = r02;
  c[1][0] = r10; c[1][1] = r11; c[1][2] = r12;
  c[2][0] = r20; c[2][1] = r21; c[2][2] = r22;
}

void SimInfo::matVecMul3(double m[3][3], double inVec[3], double outVec[3]) {
  double a0, a1, a2;

  a0 = inVec[0];  a1 = inVec[1];  a2 = inVec[2];

  outVec[0] = m[0][0]*a0 + m[0][1]*a1 + m[0][2]*a2;
  outVec[1] = m[1][0]*a0 + m[1][1]*a1 + m[1][2]*a2;
  outVec[2] = m[2][0]*a0 + m[2][1]*a1 + m[2][2]*a2;
}

void SimInfo::transposeMat3(double in[3][3], double out[3][3]) {
  double temp[3][3];
  int i, j;

  for (i = 0; i < 3; i++) {
    for (j = 0; j < 3; j++) {
      temp[j][i] = in[i][j];
    }
  }
  for (i = 0; i < 3; i++) {
    for (j = 0; j < 3; j++) {
      out[i][j] = temp[i][j];
    }
  }
}
  
void SimInfo::printMat3(double A[3][3] ){

  std::cerr 
            << "[ " << A[0][0] << ", " << A[0][1] << ", " << A[0][2] << " ]\n"
            << "[ " << A[1][0] << ", " << A[1][1] << ", " << A[1][2] << " ]\n"
            << "[ " << A[2][0] << ", " << A[2][1] << ", " << A[2][2] << " ]\n";
}

void SimInfo::printMat9(double A[9] ){

  std::cerr 
            << "[ " << A[0] << ", " << A[1] << ", " << A[2] << " ]\n"
            << "[ " << A[3] << ", " << A[4] << ", " << A[5] << " ]\n"
            << "[ " << A[6] << ", " << A[7] << ", " << A[8] << " ]\n";
}

void SimInfo::calcBoxL( void ){

  double dx, dy, dz, dsq;
  int i;

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

}


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

  int i, j, k;
  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;
  
  ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints;

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

  ndf = ndf - 3;

  return ndf;
}

int SimInfo::getNDFraw() {
  int ndfRaw_local, ndfRaw;

  // Raw degrees of freedom that we have to set
  ndfRaw_local = 3 * n_atoms + 3 * n_oriented;
  
#ifdef IS_MPI
  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
#else
  ndfRaw = ndfRaw_local;
#endif

  return ndfRaw;
}
 
void SimInfo::refreshSim(){

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

  fInfo.dielect = 0.0;

  if( useDipole ){
    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_dipoles = useDipole;
  //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;

  excl = Exclude::getArray();

#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();

}


void SimInfo::setRcut( double theRcut ){

  if( !haveOrigRcut ){
    haveOrigRcut = 1;
    origRcut = theRcut;
  }

  rCut = theRcut;
  checkCutOffs();
}

void SimInfo::setEcr( double theEcr ){

  if( !haveOrigEcr ){
    haveOrigEcr = 1;
    origEcr = theEcr;
  }

  ecr = theEcr;
  checkCutOffs();
}

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

  est = theEst;
  setEcr( theEcr );
}


void SimInfo::checkCutOffs( void ){

  int cutChanged = 0;

  if( boxIsInit ){
    
    //we need to check cutOffs against the box
    
    if( maxCutoff > rCut ){
      if( rCut < origRcut ){
        rCut = origRcut;
        if (rCut > maxCutoff) rCut = maxCutoff;
        
        sprintf( painCave.errMsg,
                 "New Box size is setting the long range cutoff radius "
                 "to %lf\n",
                 rCut );
        painCave.isFatal = 0;
        simError();
      }
    }

    if( maxCutoff > ecr ){
      if( ecr < origEcr ){
        rCut = origEcr;
        if (ecr > maxCutoff) ecr = maxCutoff;
        
        sprintf( painCave.errMsg,
                 "New Box size is setting the electrostaticCutoffRadius "
                 "to %lf\n",
                 ecr );
        painCave.isFatal = 0;
        simError();
      }
    }


    if (rCut > maxCutoff) {
      sprintf( painCave.errMsg,
               "New Box size is setting the long range cutoff radius "
               "to %lf\n",
               maxCutoff );
      painCave.isFatal = 0;
      simError();
      rCut = maxCutoff;
    }

    if( ecr > maxCutoff){
      sprintf( painCave.errMsg,
               "New Box size is setting the electrostaticCutoffRadius "
               "to %lf\n",
               maxCutoff  );
      painCave.isFatal = 0;
      simError();      
      ecr = maxCutoff;
    }

    
  }
   

  if( (oldEcr != ecr) || ( oldRcut != rCut ) ) cutChanged = 1;

  // rlist is the 1.0 plus max( rcut, ecr )
  
  ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;

  if( cutChanged ){
    
    notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
  }

  oldEcr = ecr;
  oldRcut = rCut;
}
Revision:	626
Committed:	Wed Jul 16 21:30:56 2003 UTC (20 years, 11 months ago) by mmeineke
File size:	10814 byte(s)
Log Message:	Changed how cutoffs were handled from C. Now notifyCutoffs in Fortran notifies those who need the information of any changes to cutoffs.
#	User	Rev	Content
1	mmeineke	377	#include <cstdlib>
2			#include <cstring>
3	mmeineke	568	#include <cmath>
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	490	#ifdef IS_MPI
16			#include "mpiSimulation.hpp"
17			#endif
18
19	mmeineke	572	inline double roundMe( double x ){
20			return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
21			}
22
23
24	mmeineke	377	SimInfo* currentInfo;
25
26			SimInfo::SimInfo(){
27			excludes = NULL;
28			n_constraints = 0;
29			n_oriented = 0;
30			n_dipoles = 0;
31	gezelter	458	ndf = 0;
32			ndfRaw = 0;
33	mmeineke	377	the_integrator = NULL;
34			setTemp = 0;
35			thermalTime = 0.0;
36	mmeineke	420	rCut = 0.0;
37	mmeineke	618	ecr = 0.0;
38	mmeineke	619	est = 0.0;
39	mmeineke	626	oldEcr = 0.0;
40			oldRcut = 0.0;
41	mmeineke	377
42	mmeineke	626	haveOrigRcut = 0;
43			haveOrigEcr = 0;
44			boxIsInit = 0;
45
46
47
48	mmeineke	377	usePBC = 0;
49			useLJ = 0;
50			useSticky = 0;
51			useDipole = 0;
52			useReactionField = 0;
53			useGB = 0;
54			useEAM = 0;
55
56	gezelter	457	wrapMeSimInfo( this );
57			}
58	mmeineke	377
59	gezelter	457	void SimInfo::setBox(double newBox[3]) {
60	mmeineke	586
61	gezelter	588	int i, j;
62			double tempMat[3][3];
63	gezelter	463
64	gezelter	588	for(i=0; i<3; i++)
65			for (j=0; j<3; j++) tempMat[i][j] = 0.0;;
66	gezelter	463
67	gezelter	588	tempMat[0][0] = newBox[0];
68			tempMat[1][1] = newBox[1];
69			tempMat[2][2] = newBox[2];
70	gezelter	463
71	mmeineke	586	setBoxM( tempMat );
72	mmeineke	568
73	gezelter	457	}
74	mmeineke	377
75	gezelter	588	void SimInfo::setBoxM( double theBox[3][3] ){
76	mmeineke	568
77	gezelter	588	int i, j, status;
78	mmeineke	568	double smallestBoxL, maxCutoff;
79	gezelter	588	double FortranHmat[9]; // to preserve compatibility with Fortran the
80			// ordering in the array is as follows:
81			// [ 0 3 6 ]
82			// [ 1 4 7 ]
83			// [ 2 5 8 ]
84			double FortranHmatInv[9]; // the inverted Hmat (for Fortran);
85	mmeineke	568
86	mmeineke	626
87			if( !boxIsInit ) boxIsInit = 1;
88	mmeineke	586
89	gezelter	588	for(i=0; i < 3; i++)
90			for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
91
92	mmeineke	568	calcBoxL();
93	gezelter	588	calcHmatInv();
94	mmeineke	568
95	gezelter	588	for(i=0; i < 3; i++) {
96			for (j=0; j < 3; j++) {
97			FortranHmat[3*j + i] = Hmat[i][j];
98			FortranHmatInv[3*j + i] = HmatInv[i][j];
99			}
100			}
101	mmeineke	586
102	mmeineke	590	setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic);
103	mmeineke	568
104	mmeineke	377	}
105	gezelter	458
106	mmeineke	568
107	gezelter	588	void SimInfo::getBoxM (double theBox[3][3]) {
108	mmeineke	568
109	gezelter	588	int i, j;
110			for(i=0; i<3; i++)
111			for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j];
112	mmeineke	568	}
113
114	gezelter	574
115			void SimInfo::scaleBox(double scale) {
116	gezelter	588	double theBox[3][3];
117			int i, j;
118	gezelter	574
119	gezelter	617	// cerr << "Scaling box by " << scale << "\n";
120	mmeineke	586
121	gezelter	588	for(i=0; i<3; i++)
122			for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;
123	gezelter	574
124			setBoxM(theBox);
125
126			}
127
128	gezelter	588	void SimInfo::calcHmatInv( void ) {
129	mmeineke	590
130			int i,j;
131	mmeineke	569	double smallDiag;
132			double tol;
133			double sanity[3][3];
134	mmeineke	568
135	gezelter	588	invertMat3( Hmat, HmatInv );
136	mmeineke	568
137	gezelter	588	// Check the inverse to make sure it is sane:
138	mmeineke	568
139	gezelter	588	matMul3( Hmat, HmatInv, sanity );
140
141			// check to see if Hmat is orthorhombic
142	mmeineke	568
143	gezelter	588	smallDiag = Hmat[0][0];
144			if(smallDiag > Hmat[1][1]) smallDiag = Hmat[1][1];
145			if(smallDiag > Hmat[2][2]) smallDiag = Hmat[2][2];
146			tol = smallDiag * 1E-6;
147	mmeineke	568
148	gezelter	588	orthoRhombic = 1;
149	mmeineke	568
150	gezelter	588	for (i = 0; i < 3; i++ ) {
151			for (j = 0 ; j < 3; j++) {
152			if (i != j) {
153			if (orthoRhombic) {
154			if (Hmat[i][j] >= tol) orthoRhombic = 0;
155			}
156			}
157	mmeineke	568	}
158			}
159	gezelter	588	}
160	mmeineke	569
161	gezelter	588	double SimInfo::matDet3(double a[3][3]) {
162			int i, j, k;
163			double determinant;
164	mmeineke	569
165	gezelter	588	determinant = 0.0;
166
167			for(i = 0; i < 3; i++) {
168			j = (i+1)%3;
169			k = (i+2)%3;
170
171			determinant += a[0][i] * (a[1][j]a[2][k] - a[1][k]a[2][j]);
172	mmeineke	569	}
173
174	gezelter	588	return determinant;
175			}
176	mmeineke	569
177	gezelter	588	void SimInfo::invertMat3(double a[3][3], double b[3][3]) {
178	mmeineke	569
179	gezelter	588	int i, j, k, l, m, n;
180			double determinant;
181	mmeineke	569
182	gezelter	588	determinant = matDet3( a );
183
184			if (determinant == 0.0) {
185			sprintf( painCave.errMsg,
186			"Can't invert a matrix with a zero determinant!\n");
187			painCave.isFatal = 1;
188			simError();
189			}
190
191			for (i=0; i < 3; i++) {
192			j = (i+1)%3;
193			k = (i+2)%3;
194			for(l = 0; l < 3; l++) {
195			m = (l+1)%3;
196			n = (l+2)%3;
197
198			b[l][i] = (a[j][m]a[k][n] - a[j][n]a[k][m]) / determinant;
199	mmeineke	569	}
200			}
201	mmeineke	568	}
202
203	gezelter	588	void SimInfo::matMul3(double a[3][3], double b[3][3], double c[3][3]) {
204			double r00, r01, r02, r10, r11, r12, r20, r21, r22;
205
206			r00 = a[0][0]b[0][0] + a[0][1]b[1][0] + a[0][2]*b[2][0];
207			r01 = a[0][0]b[0][1] + a[0][1]b[1][1] + a[0][2]*b[2][1];
208			r02 = a[0][0]b[0][2] + a[0][1]b[1][2] + a[0][2]*b[2][2];
209
210			r10 = a[1][0]b[0][0] + a[1][1]b[1][0] + a[1][2]*b[2][0];
211			r11 = a[1][0]b[0][1] + a[1][1]b[1][1] + a[1][2]*b[2][1];
212			r12 = a[1][0]b[0][2] + a[1][1]b[1][2] + a[1][2]*b[2][2];
213
214			r20 = a[2][0]b[0][0] + a[2][1]b[1][0] + a[2][2]*b[2][0];
215			r21 = a[2][0]b[0][1] + a[2][1]b[1][1] + a[2][2]*b[2][1];
216			r22 = a[2][0]b[0][2] + a[2][1]b[1][2] + a[2][2]*b[2][2];
217
218			c[0][0] = r00; c[0][1] = r01; c[0][2] = r02;
219			c[1][0] = r10; c[1][1] = r11; c[1][2] = r12;
220			c[2][0] = r20; c[2][1] = r21; c[2][2] = r22;
221			}
222
223			void SimInfo::matVecMul3(double m[3][3], double inVec[3], double outVec[3]) {
224			double a0, a1, a2;
225
226			a0 = inVec[0]; a1 = inVec[1]; a2 = inVec[2];
227
228			outVec[0] = m[0][0]a0 + m[0][1]a1 + m[0][2]*a2;
229			outVec[1] = m[1][0]a0 + m[1][1]a1 + m[1][2]*a2;
230			outVec[2] = m[2][0]a0 + m[2][1]a1 + m[2][2]*a2;
231			}
232	mmeineke	597
233			void SimInfo::transposeMat3(double in[3][3], double out[3][3]) {
234			double temp[3][3];
235			int i, j;
236
237			for (i = 0; i < 3; i++) {
238			for (j = 0; j < 3; j++) {
239			temp[j][i] = in[i][j];
240			}
241			}
242			for (i = 0; i < 3; i++) {
243			for (j = 0; j < 3; j++) {
244			out[i][j] = temp[i][j];
245			}
246			}
247			}
248	gezelter	588
249	mmeineke	597	void SimInfo::printMat3(double A[3][3] ){
250
251			std::cerr
252			<< "[ " << A[0][0] << ", " << A[0][1] << ", " << A[0][2] << " ]\n"
253			<< "[ " << A[1][0] << ", " << A[1][1] << ", " << A[1][2] << " ]\n"
254			<< "[ " << A[2][0] << ", " << A[2][1] << ", " << A[2][2] << " ]\n";
255			}
256
257			void SimInfo::printMat9(double A[9] ){
258
259			std::cerr
260			<< "[ " << A[0] << ", " << A[1] << ", " << A[2] << " ]\n"
261			<< "[ " << A[3] << ", " << A[4] << ", " << A[5] << " ]\n"
262			<< "[ " << A[6] << ", " << A[7] << ", " << A[8] << " ]\n";
263			}
264
265	mmeineke	568	void SimInfo::calcBoxL( void ){
266
267			double dx, dy, dz, dsq;
268			int i;
269
270	gezelter	588	// boxVol = Determinant of Hmat
271	mmeineke	568
272	gezelter	588	boxVol = matDet3( Hmat );
273	mmeineke	568
274			// boxLx
275
276	gezelter	588	dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
277	mmeineke	568	dsq = dxdx + dydy + dz*dz;
278	gezelter	621	boxL[0] = sqrt( dsq );
279	mmeineke	626	maxCutoff = 0.5 * boxL[0];
280	mmeineke	568
281			// boxLy
282
283	gezelter	588	dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
284	mmeineke	568	dsq = dxdx + dydy + dz*dz;
285	gezelter	621	boxL[1] = sqrt( dsq );
286	mmeineke	626	if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];
287	mmeineke	568
288			// boxLz
289
290	gezelter	588	dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
291	mmeineke	568	dsq = dxdx + dydy + dz*dz;
292	gezelter	621	boxL[2] = sqrt( dsq );
293	mmeineke	626	if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];
294
295	mmeineke	568	}
296
297
298			void SimInfo::wrapVector( double thePos[3] ){
299
300			int i, j, k;
301			double scaled[3];
302
303	mmeineke	569	if( !orthoRhombic ){
304			// calc the scaled coordinates.
305	gezelter	588
306
307			matVecMul3(HmatInv, thePos, scaled);
308	mmeineke	569
309			for(i=0; i<3; i++)
310	mmeineke	572	scaled[i] -= roundMe(scaled[i]);
311	mmeineke	569
312			// calc the wrapped real coordinates from the wrapped scaled coordinates
313
314	gezelter	588	matVecMul3(Hmat, scaled, thePos);
315
316	mmeineke	569	}
317			else{
318			// calc the scaled coordinates.
319
320			for(i=0; i<3; i++)
321	gezelter	588	scaled[i] = thePos[i]*HmatInv[i][i];
322	mmeineke	569
323			// wrap the scaled coordinates
324
325			for(i=0; i<3; i++)
326	mmeineke	572	scaled[i] -= roundMe(scaled[i]);
327	mmeineke	569
328			// calc the wrapped real coordinates from the wrapped scaled coordinates
329
330			for(i=0; i<3; i++)
331	gezelter	588	thePos[i] = scaled[i]*Hmat[i][i];
332	mmeineke	569	}
333
334	mmeineke	568	}
335
336
337	gezelter	458	int SimInfo::getNDF(){
338			int ndf_local, ndf;
339	gezelter	457
340	gezelter	458	ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints;
341
342			#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			ndf = ndf - 3;
349
350			return ndf;
351			}
352
353			int SimInfo::getNDFraw() {
354			int ndfRaw_local, ndfRaw;
355
356			// Raw degrees of freedom that we have to set
357			ndfRaw_local = 3 * n_atoms + 3 * n_oriented;
358
359			#ifdef IS_MPI
360			MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
361			#else
362			ndfRaw = ndfRaw_local;
363			#endif
364
365			return ndfRaw;
366			}
367
368	mmeineke	377	void SimInfo::refreshSim(){
369
370			simtype fInfo;
371			int isError;
372	gezelter	490	int n_global;
373	mmeineke	424	int* excl;
374	mmeineke	626
375	mmeineke	469	fInfo.dielect = 0.0;
376	mmeineke	377
377	mmeineke	469	if( useDipole ){
378			if( useReactionField )fInfo.dielect = dielectric;
379			}
380
381	mmeineke	377	fInfo.SIM_uses_PBC = usePBC;
382	mmeineke	443	//fInfo.SIM_uses_LJ = 0;
383	chuckv	439	fInfo.SIM_uses_LJ = useLJ;
384	mmeineke	443	fInfo.SIM_uses_sticky = useSticky;
385			//fInfo.SIM_uses_sticky = 0;
386	chuckv	482	fInfo.SIM_uses_dipoles = useDipole;
387			//fInfo.SIM_uses_dipoles = 0;
388	mmeineke	443	//fInfo.SIM_uses_RF = useReactionField;
389			fInfo.SIM_uses_RF = 0;
390	mmeineke	377	fInfo.SIM_uses_GB = useGB;
391			fInfo.SIM_uses_EAM = useEAM;
392
393	mmeineke	424	excl = Exclude::getArray();
394	mmeineke	377
395	gezelter	490	#ifdef IS_MPI
396			n_global = mpiSim->getTotAtoms();
397			#else
398			n_global = n_atoms;
399			#endif
400
401	mmeineke	377	isError = 0;
402
403	gezelter	490	setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
404	gezelter	483	&nGlobalExcludes, globalExcludes, molMembershipArray,
405			&isError );
406	mmeineke	377
407			if( isError ){
408
409			sprintf( painCave.errMsg,
410			"There was an error setting the simulation information in fortran.\n" );
411			painCave.isFatal = 1;
412			simError();
413			}
414
415			#ifdef IS_MPI
416			sprintf( checkPointMsg,
417			"succesfully sent the simulation information to fortran.\n");
418			MPIcheckPoint();
419			#endif // is_mpi
420	gezelter	458
421	gezelter	474	this->ndf = this->getNDF();
422			this->ndfRaw = this->getNDFraw();
423	gezelter	458
424	mmeineke	377	}
425
426	mmeineke	626
427			void SimInfo::setRcut( double theRcut ){
428
429			if( !haveOrigRcut ){
430			haveOrigRcut = 1;
431			origRcut = theRcut;
432			}
433
434			rCut = theRcut;
435			checkCutOffs();
436			}
437
438			void SimInfo::setEcr( double theEcr ){
439
440			if( !haveOrigEcr ){
441			haveOrigEcr = 1;
442			origEcr = theEcr;
443			}
444
445			ecr = theEcr;
446			checkCutOffs();
447			}
448
449			void SimInfo::setEcr( double theEcr, double theEst ){
450
451			est = theEst;
452			setEcr( theEcr );
453			}
454
455
456			void SimInfo::checkCutOffs( void ){
457
458			int cutChanged = 0;
459
460			if( boxIsInit ){
461
462			//we need to check cutOffs against the box
463
464			if( maxCutoff > rCut ){
465			if( rCut < origRcut ){
466			rCut = origRcut;
467			if (rCut > maxCutoff) rCut = maxCutoff;
468
469			sprintf( painCave.errMsg,
470			"New Box size is setting the long range cutoff radius "
471			"to %lf\n",
472			rCut );
473			painCave.isFatal = 0;
474			simError();
475			}
476			}
477
478			if( maxCutoff > ecr ){
479			if( ecr < origEcr ){
480			rCut = origEcr;
481			if (ecr > maxCutoff) ecr = maxCutoff;
482
483			sprintf( painCave.errMsg,
484			"New Box size is setting the electrostaticCutoffRadius "
485			"to %lf\n",
486			ecr );
487			painCave.isFatal = 0;
488			simError();
489			}
490			}
491
492
493			if (rCut > maxCutoff) {
494			sprintf( painCave.errMsg,
495			"New Box size is setting the long range cutoff radius "
496			"to %lf\n",
497			maxCutoff );
498			painCave.isFatal = 0;
499			simError();
500			rCut = maxCutoff;
501			}
502
503			if( ecr > maxCutoff){
504			sprintf( painCave.errMsg,
505			"New Box size is setting the electrostaticCutoffRadius "
506			"to %lf\n",
507			maxCutoff );
508			painCave.isFatal = 0;
509			simError();
510			ecr = maxCutoff;
511			}
512
513
514			}
515
516
517			if( (oldEcr != ecr) \|\| ( oldRcut != rCut ) ) cutChanged = 1;
518
519			// rlist is the 1.0 plus max( rcut, ecr )
520
521			( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
522
523			if( cutChanged ){
524
525			notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
526			}
527
528			oldEcr = ecr;
529			oldRcut = rCut;
530			}