1 |
< |
#include <cstdlib> |
2 |
< |
#include <cstring> |
3 |
< |
#include <cmath> |
1 |
> |
#include <stdlib.h> |
2 |
> |
#include <string.h> |
3 |
> |
#include <math.h> |
4 |
|
|
5 |
|
#include <iostream> |
6 |
|
using namespace std; |
49 |
|
boxIsInit = 0; |
50 |
|
|
51 |
|
resetTime = 1e99; |
52 |
– |
|
52 |
|
|
53 |
+ |
orthoTolerance = 1E-6; |
54 |
+ |
useInitXSstate = true; |
55 |
+ |
|
56 |
|
usePBC = 0; |
57 |
|
useLJ = 0; |
58 |
|
useSticky = 0; |
96 |
|
|
97 |
|
void SimInfo::setBoxM( double theBox[3][3] ){ |
98 |
|
|
99 |
< |
int i, j, status; |
98 |
< |
double smallestBoxL, maxCutoff; |
99 |
> |
int i, j; |
100 |
|
double FortranHmat[9]; // to preserve compatibility with Fortran the |
101 |
|
// ordering in the array is as follows: |
102 |
|
// [ 0 3 6 ] |
148 |
|
|
149 |
|
void SimInfo::calcHmatInv( void ) { |
150 |
|
|
151 |
+ |
int oldOrtho; |
152 |
|
int i,j; |
153 |
|
double smallDiag; |
154 |
|
double tol; |
156 |
|
|
157 |
|
invertMat3( Hmat, HmatInv ); |
158 |
|
|
157 |
– |
// Check the inverse to make sure it is sane: |
158 |
– |
|
159 |
– |
matMul3( Hmat, HmatInv, sanity ); |
160 |
– |
|
159 |
|
// check to see if Hmat is orthorhombic |
160 |
|
|
161 |
< |
smallDiag = Hmat[0][0]; |
164 |
< |
if(smallDiag > Hmat[1][1]) smallDiag = Hmat[1][1]; |
165 |
< |
if(smallDiag > Hmat[2][2]) smallDiag = Hmat[2][2]; |
166 |
< |
tol = smallDiag * 1E-6; |
161 |
> |
oldOrtho = orthoRhombic; |
162 |
|
|
163 |
+ |
smallDiag = fabs(Hmat[0][0]); |
164 |
+ |
if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]); |
165 |
+ |
if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]); |
166 |
+ |
tol = smallDiag * orthoTolerance; |
167 |
+ |
|
168 |
|
orthoRhombic = 1; |
169 |
|
|
170 |
|
for (i = 0; i < 3; i++ ) { |
171 |
|
for (j = 0 ; j < 3; j++) { |
172 |
|
if (i != j) { |
173 |
|
if (orthoRhombic) { |
174 |
< |
if (Hmat[i][j] >= tol) orthoRhombic = 0; |
174 |
> |
if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0; |
175 |
|
} |
176 |
|
} |
177 |
+ |
} |
178 |
+ |
} |
179 |
+ |
|
180 |
+ |
if( oldOrtho != orthoRhombic ){ |
181 |
+ |
|
182 |
+ |
if( orthoRhombic ){ |
183 |
+ |
sprintf( painCave.errMsg, |
184 |
+ |
"Hmat is switching from Non-Orthorhombic to OrthoRhombic\n" |
185 |
+ |
" If this is a bad thing, change the orthoBoxTolerance( currently %G ).\n", |
186 |
+ |
orthoTolerance); |
187 |
+ |
simError(); |
188 |
|
} |
189 |
+ |
else { |
190 |
+ |
sprintf( painCave.errMsg, |
191 |
+ |
"Hmat is switching from Orthorhombic to Non-OrthoRhombic\n" |
192 |
+ |
" If this is a bad thing, change the orthoBoxTolerance( currently %G ).\n", |
193 |
+ |
orthoTolerance); |
194 |
+ |
simError(); |
195 |
+ |
} |
196 |
|
} |
197 |
|
} |
198 |
|
|
320 |
|
void SimInfo::calcBoxL( void ){ |
321 |
|
|
322 |
|
double dx, dy, dz, dsq; |
305 |
– |
int i; |
323 |
|
|
324 |
|
// boxVol = Determinant of Hmat |
325 |
|
|
389 |
|
|
390 |
|
void SimInfo::wrapVector( double thePos[3] ){ |
391 |
|
|
392 |
< |
int i, j, k; |
392 |
> |
int i; |
393 |
|
double scaled[3]; |
394 |
|
|
395 |
|
if( !orthoRhombic ){ |
427 |
|
|
428 |
|
|
429 |
|
int SimInfo::getNDF(){ |
430 |
< |
int ndf_local, ndf; |
430 |
> |
int ndf_local; |
431 |
|
|
432 |
|
ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints; |
433 |
|
|
443 |
|
} |
444 |
|
|
445 |
|
int SimInfo::getNDFraw() { |
446 |
< |
int ndfRaw_local, ndfRaw; |
446 |
> |
int ndfRaw_local; |
447 |
|
|
448 |
|
// Raw degrees of freedom that we have to set |
449 |
|
ndfRaw_local = 3 * n_atoms + 3 * n_oriented; |
458 |
|
} |
459 |
|
|
460 |
|
int SimInfo::getNDFtranslational() { |
461 |
< |
int ndfTrans_local, ndfTrans; |
461 |
> |
int ndfTrans_local; |
462 |
|
|
463 |
|
ndfTrans_local = 3 * n_atoms - n_constraints; |
464 |
|
|
534 |
|
|
535 |
|
void SimInfo::setRcut( double theRcut ){ |
536 |
|
|
520 |
– |
if( !haveOrigRcut ){ |
521 |
– |
haveOrigRcut = 1; |
522 |
– |
origRcut = theRcut; |
523 |
– |
} |
524 |
– |
|
537 |
|
rCut = theRcut; |
538 |
|
checkCutOffs(); |
539 |
|
} |
540 |
|
|
541 |
< |
void SimInfo::setEcr( double theEcr ){ |
541 |
> |
void SimInfo::setDefaultRcut( double theRcut ){ |
542 |
|
|
543 |
< |
if( !haveOrigEcr ){ |
544 |
< |
haveOrigEcr = 1; |
545 |
< |
origEcr = theEcr; |
534 |
< |
} |
543 |
> |
haveOrigRcut = 1; |
544 |
> |
origRcut = theRcut; |
545 |
> |
rCut = theRcut; |
546 |
|
|
547 |
+ |
( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0; |
548 |
+ |
|
549 |
+ |
notifyFortranCutOffs( &rCut, &rList, &ecr, &est ); |
550 |
+ |
} |
551 |
+ |
|
552 |
+ |
void SimInfo::setEcr( double theEcr ){ |
553 |
+ |
|
554 |
|
ecr = theEcr; |
555 |
|
checkCutOffs(); |
556 |
|
} |
557 |
|
|
558 |
+ |
void SimInfo::setDefaultEcr( double theEcr ){ |
559 |
+ |
|
560 |
+ |
haveOrigEcr = 1; |
561 |
+ |
origEcr = theEcr; |
562 |
+ |
|
563 |
+ |
( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0; |
564 |
+ |
|
565 |
+ |
ecr = theEcr; |
566 |
+ |
|
567 |
+ |
notifyFortranCutOffs( &rCut, &rList, &ecr, &est ); |
568 |
+ |
} |
569 |
+ |
|
570 |
|
void SimInfo::setEcr( double theEcr, double theEst ){ |
571 |
|
|
572 |
|
est = theEst; |
573 |
|
setEcr( theEcr ); |
574 |
|
} |
575 |
|
|
576 |
+ |
void SimInfo::setDefaultEcr( double theEcr, double theEst ){ |
577 |
|
|
578 |
+ |
est = theEst; |
579 |
+ |
setDefaultEcr( theEcr ); |
580 |
+ |
} |
581 |
+ |
|
582 |
+ |
|
583 |
|
void SimInfo::checkCutOffs( void ){ |
584 |
|
|
585 |
|
int cutChanged = 0; |
646 |
|
( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0; |
647 |
|
|
648 |
|
if( cutChanged ){ |
613 |
– |
|
649 |
|
notifyFortranCutOffs( &rCut, &rList, &ecr, &est ); |
650 |
|
} |
651 |
|
|