91 |
|
} |
92 |
|
|
93 |
|
void NPAT::scaleSimBox(){ |
94 |
– |
|
95 |
– |
int i; |
96 |
– |
int j; |
97 |
– |
int k; |
94 |
|
Mat3x3d scaleMat; |
99 |
– |
double eta2ij; |
100 |
– |
double bigScale, smallScale, offDiagMax; |
101 |
– |
Mat3x3d hm; |
102 |
– |
Mat3x3d hmnew; |
95 |
|
|
96 |
< |
|
97 |
< |
|
98 |
< |
// Scale the box after all the positions have been moved: |
99 |
< |
|
100 |
< |
// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat) |
101 |
< |
// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2) |
110 |
< |
|
111 |
< |
bigScale = 1.0; |
112 |
< |
smallScale = 1.0; |
113 |
< |
offDiagMax = 0.0; |
114 |
< |
|
115 |
< |
for(i=0; i<3; i++){ |
116 |
< |
for(j=0; j<3; j++){ |
117 |
< |
|
118 |
< |
// Calculate the matrix Product of the eta array (we only need |
119 |
< |
// the ij element right now): |
120 |
< |
|
121 |
< |
eta2ij = 0.0; |
122 |
< |
for(k=0; k<3; k++){ |
123 |
< |
eta2ij += eta(i, k) * eta(k, j); |
124 |
< |
} |
125 |
< |
|
126 |
< |
scaleMat(i, j) = 0.0; |
127 |
< |
// identity matrix (see above): |
128 |
< |
if (i == j) scaleMat(i, j) = 1.0; |
129 |
< |
// Taylor expansion for the exponential truncated at second order: |
130 |
< |
scaleMat(i, j) += dt*eta(i, j) + 0.5*dt*dt*eta2ij; |
131 |
< |
|
132 |
< |
|
133 |
< |
if (i != j) |
134 |
< |
if (fabs(scaleMat(i, j)) > offDiagMax) |
135 |
< |
offDiagMax = fabs(scaleMat(i, j)); |
96 |
> |
for(int i=0; i<3; i++){ |
97 |
> |
for(int j=0; j<3; j++){ |
98 |
> |
scaleMat(i, j) = 0.0; |
99 |
> |
if(i==j) { |
100 |
> |
scaleMat(i, j) = 1.0; |
101 |
> |
} |
102 |
|
} |
137 |
– |
|
138 |
– |
if (scaleMat(i, i) > bigScale) bigScale = scaleMat(i, i); |
139 |
– |
if (scaleMat(i, i) < smallScale) smallScale = scaleMat(i, i); |
103 |
|
} |
104 |
< |
|
105 |
< |
if ((bigScale > 1.01) || (smallScale < 0.99)) { |
106 |
< |
sprintf( painCave.errMsg, |
107 |
< |
"NPAT error: Attempting a Box scaling of more than 1 percent.\n" |
108 |
< |
" Check your tauBarostat, as it is probably too small!\n\n" |
146 |
< |
" scaleMat = [%lf\t%lf\t%lf]\n" |
147 |
< |
" [%lf\t%lf\t%lf]\n" |
148 |
< |
" [%lf\t%lf\t%lf]\n" |
149 |
< |
" eta = [%lf\t%lf\t%lf]\n" |
150 |
< |
" [%lf\t%lf\t%lf]\n" |
151 |
< |
" [%lf\t%lf\t%lf]\n", |
152 |
< |
scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2), |
153 |
< |
scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2), |
154 |
< |
scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2), |
155 |
< |
eta(0, 0),eta(0, 1),eta(0, 2), |
156 |
< |
eta(1, 0),eta(1, 1),eta(1, 2), |
157 |
< |
eta(2, 0),eta(2, 1),eta(2, 2)); |
158 |
< |
painCave.isFatal = 1; |
159 |
< |
simError(); |
160 |
< |
} else if (offDiagMax > 0.01) { |
161 |
< |
sprintf( painCave.errMsg, |
162 |
< |
"NPAT error: Attempting an off-diagonal Box scaling of more than 1 percent.\n" |
163 |
< |
" Check your tauBarostat, as it is probably too small!\n\n" |
164 |
< |
" scaleMat = [%lf\t%lf\t%lf]\n" |
165 |
< |
" [%lf\t%lf\t%lf]\n" |
166 |
< |
" [%lf\t%lf\t%lf]\n" |
167 |
< |
" eta = [%lf\t%lf\t%lf]\n" |
168 |
< |
" [%lf\t%lf\t%lf]\n" |
169 |
< |
" [%lf\t%lf\t%lf]\n", |
170 |
< |
scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2), |
171 |
< |
scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2), |
172 |
< |
scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2), |
173 |
< |
eta(0, 0),eta(0, 1),eta(0, 2), |
174 |
< |
eta(1, 0),eta(1, 1),eta(1, 2), |
175 |
< |
eta(2, 0),eta(2, 1),eta(2, 2)); |
176 |
< |
painCave.isFatal = 1; |
177 |
< |
simError(); |
178 |
< |
} else { |
179 |
< |
|
180 |
< |
Mat3x3d hmat = currentSnapshot_->getHmat(); |
181 |
< |
hmat = hmat *scaleMat; |
182 |
< |
currentSnapshot_->setHmat(hmat); |
183 |
< |
|
184 |
< |
} |
104 |
> |
|
105 |
> |
scaleMat(2, 2) = exp(dt*eta(2, 2)); |
106 |
> |
Mat3x3d hmat = currentSnapshot_->getHmat(); |
107 |
> |
hmat = hmat *scaleMat; |
108 |
> |
currentSnapshot_->setHmat(hmat); |
109 |
|
} |
110 |
|
|
111 |
|
bool NPAT::etaConverged() { |