Newer
Older
/*----------------------------------------------------------------------
SerialReax - Reax Force Field Simulator
Copyright (2010) Purdue University
Hasan Metin Aktulga, haktulga@cs.purdue.edu
Joseph Fogarty, jcfogart@mail.usf.edu
Sagar Pandit, pandit@usf.edu
Ananth Y Grama, ayg@cs.purdue.edu
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2 of
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
See the GNU General Public License for more details:
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "box.h"
#include "vector.h"
void Init_Box_From_CRYST(real a, real b, real c,
real alpha, real beta, real gamma,
simulation_box* box )
c_alpha = cos(DEG2RAD(alpha));
c_beta = cos(DEG2RAD(beta));
c_gamma = cos(DEG2RAD(gamma));
s_gamma = sin(DEG2RAD(gamma));
box->box[0][0] = a;
box->box[0][1] = 0.0;
box->box[0][2] = 0.0;
box->box[1][0] = b * c_gamma;
box->box[1][1] = b * s_gamma;
box->box[1][2] = 0.0;
box->box[2][0] = c * c_beta;
box->box[2][1] = c * zi;
box->box[2][2] = c * SQRT(1.0 - SQR(c_beta) - SQR(zi));
Make_Consistent( box );
fprintf( stderr, "box is %8.2f x %8.2f x %8.2f\n",
box->box[0][0], box->box[1][1], box->box[2][2] );
#endif
}
void Update_Box( rtensor box_tensor, simulation_box* box )
{
for (i = 0; i < 3; i++)
for (j = 0; j < 3; j++)
box->box[i][j] = box_tensor[i][j];
}
void Update_Box_Isotropic( simulation_box *box, real mu )
{
/*box->box[0][0] =
POW( V_new / ( box->side_prop[1] * box->side_prop[2] ), 1.0/3.0 );
box->box[1][1] = box->box[0][0] * box->side_prop[1];
box->box[2][2] = box->box[0][0] * box->side_prop[2];
*/
rtensor_Copy( box->old_box, box->box );
box->box[0][0] *= mu;
box->box[1][1] *= mu;
box->box[2][2] *= mu;
box->volume = box->box[0][0] * box->box[1][1] * box->box[2][2];
Make_Consistent(box/*, periodic*/);
}
void Update_Box_SemiIsotropic( simulation_box *box, rvec mu )
{
/*box->box[0][0] =
POW( V_new / ( box->side_prop[1] * box->side_prop[2] ), 1.0/3.0 );
box->box[1][1] = box->box[0][0] * box->side_prop[1];
box->box[2][2] = box->box[0][0] * box->side_prop[2]; */
rtensor_Copy( box->old_box, box->box );
box->box[0][0] *= mu[0];
box->box[1][1] *= mu[1];
box->box[2][2] *= mu[2];
box->volume = box->box[0][0] * box->box[1][1] * box->box[2][2];
Make_Consistent(box);
}
void Make_Consistent(simulation_box* box)
{
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
real one_vol;
box->volume =
box->box[0][0] * (box->box[1][1] * box->box[2][2] -
box->box[2][1] * box->box[2][1]) +
box->box[0][1] * (box->box[2][0] * box->box[1][2] -
box->box[1][0] * box->box[2][2]) +
box->box[0][2] * (box->box[1][0] * box->box[2][1] -
box->box[2][0] * box->box[1][1]);
one_vol = 1.0 / box->volume;
box->box_inv[0][0] = (box->box[1][1] * box->box[2][2] -
box->box[1][2] * box->box[2][1]) * one_vol;
box->box_inv[0][1] = (box->box[0][2] * box->box[2][1] -
box->box[0][1] * box->box[2][2]) * one_vol;
box->box_inv[0][2] = (box->box[0][1] * box->box[1][2] -
box->box[0][2] * box->box[1][1]) * one_vol;
box->box_inv[1][0] = (box->box[1][2] * box->box[2][0] -
box->box[1][0] * box->box[2][2]) * one_vol;
box->box_inv[1][1] = (box->box[0][0] * box->box[2][2] -
box->box[0][2] * box->box[2][0]) * one_vol;
box->box_inv[1][2] = (box->box[0][2] * box->box[1][0] -
box->box[0][0] * box->box[1][2]) * one_vol;
box->box_inv[2][0] = (box->box[1][0] * box->box[2][1] -
box->box[1][1] * box->box[2][0]) * one_vol;
box->box_inv[2][1] = (box->box[0][1] * box->box[2][0] -
box->box[0][0] * box->box[2][1]) * one_vol;
box->box_inv[2][2] = (box->box[0][0] * box->box[1][1] -
box->box[0][1] * box->box[1][0]) * one_vol;
box->box_norms[0] = SQRT( SQR(box->box[0][0]) +
SQR(box->box[0][1]) +
SQR(box->box[0][2]) );
box->box_norms[1] = SQRT( SQR(box->box[1][0]) +
SQR(box->box[1][1]) +
SQR(box->box[1][2]) );
box->box_norms[2] = SQRT( SQR(box->box[2][0]) +
SQR(box->box[2][1]) +
SQR(box->box[2][2]) );
box->trans[0][0] = box->box[0][0] / box->box_norms[0];
box->trans[0][1] = box->box[1][0] / box->box_norms[0];
box->trans[0][2] = box->box[2][0] / box->box_norms[0];
box->trans[1][0] = box->box[0][1] / box->box_norms[1];
box->trans[1][1] = box->box[1][1] / box->box_norms[1];
box->trans[1][2] = box->box[2][1] / box->box_norms[1];
box->trans[2][0] = box->box[0][2] / box->box_norms[2];
box->trans[2][1] = box->box[1][2] / box->box_norms[2];
box->trans[2][2] = box->box[2][2] / box->box_norms[2];
one_vol = box->box_norms[0] * box->box_norms[1] * box->box_norms[2] * one_vol;
box->trans_inv[0][0] = (box->trans[1][1] * box->trans[2][2] -
box->trans[1][2] * box->trans[2][1]) * one_vol;
box->trans_inv[0][1] = (box->trans[0][2] * box->trans[2][1] -
box->trans[0][1] * box->trans[2][2]) * one_vol;
box->trans_inv[0][2] = (box->trans[0][1] * box->trans[1][2] -
box->trans[0][2] * box->trans[1][1]) * one_vol;
box->trans_inv[1][0] = (box->trans[1][2] * box->trans[2][0] -
box->trans[1][0] * box->trans[2][2]) * one_vol;
box->trans_inv[1][1] = (box->trans[0][0] * box->trans[2][2] -
box->trans[0][2] * box->trans[2][0]) * one_vol;
box->trans_inv[1][2] = (box->trans[0][2] * box->trans[1][0] -
box->trans[0][0] * box->trans[1][2]) * one_vol;
box->trans_inv[2][0] = (box->trans[1][0] * box->trans[2][1] -
box->trans[1][1] * box->trans[2][0]) * one_vol;
box->trans_inv[2][1] = (box->trans[0][1] * box->trans[2][0] -
box->trans[0][0] * box->trans[2][1]) * one_vol;
box->trans_inv[2][2] = (box->trans[0][0] * box->trans[1][1] -
box->trans[0][1] * box->trans[1][0]) * one_vol;
// for (i=0; i < 3; i++)
// {
// for (j=0; j < 3; j++)
// fprintf(stderr,"\n");
// }
// fprintf(stderr,"\n");
// for (i=0; i < 3; i++)
// {
// for (j=0; j < 3; j++)
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
box->g[0][0] = box->box[0][0] * box->box[0][0] +
box->box[0][1] * box->box[0][1] +
box->box[0][2] * box->box[0][2];
box->g[1][0] =
box->g[0][1] = box->box[0][0] * box->box[1][0] +
box->box[0][1] * box->box[1][1] +
box->box[0][2] * box->box[1][2];
box->g[2][0] =
box->g[0][2] = box->box[0][0] * box->box[2][0] +
box->box[0][1] * box->box[2][1] +
box->box[0][2] * box->box[2][2];
box->g[1][1] = box->box[1][0] * box->box[1][0] +
box->box[1][1] * box->box[1][1] +
box->box[1][2] * box->box[1][2];
box->g[1][2] =
box->g[2][1] = box->box[1][0] * box->box[2][0] +
box->box[1][1] * box->box[2][1] +
box->box[1][2] * box->box[2][2];
box->g[2][2] = box->box[2][0] * box->box[2][0] +
box->box[2][1] * box->box[2][1] +
box->box[2][2] * box->box[2][2];
// These proportions are only used for isotropic_NPT!
box->side_prop[0] = box->box[0][0] / box->box[0][0];
box->side_prop[1] = box->box[1][1] / box->box[0][0];
box->side_prop[2] = box->box[2][2] / box->box[0][0];
}
void Transform( rvec x1, simulation_box *box, char flag, rvec x2 )
{
int i, j;
real tmp;
// printf(">x1: (%lf, %lf, %lf)\n",x1[0],x1[1],x1[2]);
if (flag > 0)
{
for (i = 0; i < 3; i++)
{
tmp = 0.0;
for (j = 0; j < 3; j++)
tmp += box->trans[i][j] * x1[j];
x2[i] = tmp;
}
else
{
for (i = 0; i < 3; i++)
{
tmp = 0.0;
for (j = 0; j < 3; j++)
tmp += box->trans_inv[i][j] * x1[j];
x2[i] = tmp;
}
}
void Transform_to_UnitBox( rvec x1, simulation_box *box, char flag, rvec x2 )
{
Transform( x1, box, flag, x2 );
x2[0] /= box->box_norms[0];
x2[1] /= box->box_norms[1];
x2[2] /= box->box_norms[2];
}
void Inc_on_T3( rvec x, rvec dx, simulation_box *box )
{
for (i = 0; i < 3; i++)
{
tmp = x[i] + dx[i];
if ( tmp <= -box->box_norms[i] || tmp >= box->box_norms[i] )
tmp = fmod( tmp, box->box_norms[i] );
if ( tmp < 0 ) tmp += box->box_norms[i];
x[i] = tmp;
}
}
real Sq_Distance_on_T3(rvec x1, rvec x2, simulation_box* box, rvec r)
{
real norm = 0.0;
real d, tmp;
int i;
for (i = 0; i < 3; i++)
{
d = x2[i] - x1[i];
tmp = SQR(d);
if ( tmp >= SQR( box->box_norms[i] / 2.0 ) )
{
if (x2[i] > x1[i])
d -= box->box_norms[i];
else
d += box->box_norms[i];
r[i] = d;
norm += SQR(d);
}
else
{
r[i] = d;
norm += tmp;
}
}
void Distance_on_T3_Gen( rvec x1, rvec x2, simulation_box* box, rvec r )
{
rvec xa, xb, ra;
Transform( x1, box, -1, xa );
Transform( x2, box, -1, xb );
//printf(">xa: (%lf, %lf, %lf)\n",xa[0],xa[1],xa[2]);
//printf(">xb: (%lf, %lf, %lf)\n",xb[0],xb[1],xb[2]);
}
void Inc_on_T3_Gen( rvec x, rvec dx, simulation_box* box )
{
Transform( x, box, -1, xa );
Transform( dx, box, -1, dxa );
//printf(">xa: (%lf, %lf, %lf)\n",xa[0],xa[1],xa[2]);
//printf(">dxa: (%lf, %lf, %lf)\n",dxa[0],dxa[1],dxa[2]);
//printf(">new_xa: (%lf, %lf, %lf)\n",xa[0],xa[1],xa[2]);
Transform( xa, box, 1, x );
}
real Metric_Product( rvec x1, rvec x2, simulation_box* box )
{
tmp = 0.0;
for ( j = 0; j < 3; j++ )
tmp += box->g[i][j] * x2[j];
dist += x1[i] * tmp;
int Are_Far_Neighbors( rvec x1, rvec x2, simulation_box *box,
real cutoff, far_neighbor_data *data )
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
real norm_sqr, d, tmp;
int i;
norm_sqr = 0;
for ( i = 0; i < 3; i++ )
{
d = x2[i] - x1[i];
tmp = SQR(d);
if ( tmp >= SQR( box->box_norms[i] / 2.0 ) )
{
if ( x2[i] > x1[i] )
{
d -= box->box_norms[i];
data->rel_box[i] = -1;
}
else
{
d += box->box_norms[i];
data->rel_box[i] = +1;
}
data->dvec[i] = d;
norm_sqr += SQR(d);
}
else
{
data->dvec[i] = d;
norm_sqr += tmp;
data->rel_box[i] = 0;
}
if ( norm_sqr <= SQR(cutoff) )
{
data->d = sqrt(norm_sqr);
return 1;
}
return 0;
/* Determines if the distance between x1 and x2 is < vlist_cut.
If so, this neighborhood is added to the list of far neighbors.
Periodic boundary conditions do not apply. */
void Get_NonPeriodic_Far_Neighbors( rvec x1, rvec x2, simulation_box *box,
control_params *control,
far_neighbor_data *new_nbrs, int *count )
real norm_sqr;
rvec_ScaledSum( new_nbrs[0].dvec, 1.0, x2, -1.0, x1 );
norm_sqr = rvec_Norm_Sqr( new_nbrs[0].dvec );
if ( norm_sqr <= SQR( control->vlist_cut ) )
{
*count = 1;
new_nbrs[0].d = SQRT( norm_sqr );
ivec_MakeZero( new_nbrs[0].rel_box );
// rvec_MakeZero( new_nbrs[0].ext_factor );
}
else *count = 0;
}
/* Finds periodic neighbors in a 'big_box'. Here 'big_box' means:
the current simulation box has all dimensions > 2 *vlist_cut.
If the periodic distance between x1 and x2 is than vlist_cut, this
neighborhood is added to the list of far neighbors. */
void Get_Periodic_Far_Neighbors_Big_Box( rvec x1, rvec x2, simulation_box *box,
control_params *control,
far_neighbor_data *periodic_nbrs,
int *count )
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
real norm_sqr, d, tmp;
int i;
norm_sqr = 0;
for ( i = 0; i < 3; i++ )
{
d = x2[i] - x1[i];
tmp = SQR(d);
// fprintf(out,"Inside Sq_Distance_on_T3, %d, %lf, %lf\n",
// i,tmp,SQR(box->box_norms[i]/2.0));
if ( tmp >= SQR( box->box_norms[i] / 2.0 ) )
{
if ( x2[i] > x1[i] )
{
d -= box->box_norms[i];
periodic_nbrs[0].rel_box[i] = -1;
// periodic_nbrs[0].ext_factor[i] = +1;
}
else
{
d += box->box_norms[i];
periodic_nbrs[0].rel_box[i] = +1;
// periodic_nbrs[0].ext_factor[i] = -1;
}
periodic_nbrs[0].dvec[i] = d;
norm_sqr += SQR(d);
}
else
{
periodic_nbrs[0].dvec[i] = d;
norm_sqr += tmp;
periodic_nbrs[0].rel_box[i] = 0;
// periodic_nbrs[0].ext_factor[i] = 0;
}
}
if ( norm_sqr <= SQR( control->vlist_cut ) )
{
*count = 1;
periodic_nbrs[0].d = SQRT( norm_sqr );
((dist(x1, x2') < vlist_cut, periodic images of x2 are also considered).
Here the box is 'small' meaning that at least one dimension is < 2*vlist_cut.
IMPORTANT: This part might need some improvement. In NPT, the simulation box
might get too small (such as <5 A!). In this case we have to consider the
periodic images of x2 that are two boxs away!!!
*/
void Get_Periodic_Far_Neighbors_Small_Box( rvec x1, rvec x2, simulation_box *box,
control_params *control,
far_neighbor_data *periodic_nbrs,
int *count )
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
int i, j, k;
int imax, jmax, kmax;
real sqr_norm, d_i, d_j, d_k;
*count = 0;
/* determine the max stretch of imaginary boxs in each direction
to handle periodic boundary conditions correctly. */
imax = (int)(control->vlist_cut / box->box_norms[0] + 1);
jmax = (int)(control->vlist_cut / box->box_norms[1] + 1);
kmax = (int)(control->vlist_cut / box->box_norms[2] + 1);
/*if( imax > 1 || jmax > 1 || kmax > 1 )
fprintf( stderr, "box %8.3f x %8.3f x %8.3f --> %2d %2d %2d\n",
box->box_norms[0], box->box_norms[1], box->box_norms[2],
imax, jmax, kmax ); */
for ( i = -imax; i <= imax; ++i )
if (fabs(d_i = ((x2[0] + i * box->box_norms[0]) - x1[0])) <= control->vlist_cut)
{
for ( j = -jmax; j <= jmax; ++j )
if (fabs(d_j = ((x2[1] + j * box->box_norms[1]) - x1[1])) <= control->vlist_cut)
{
for ( k = -kmax; k <= kmax; ++k )
if (fabs(d_k = ((x2[2] + k * box->box_norms[2]) - x1[2])) <= control->vlist_cut)
{
sqr_norm = SQR(d_i) + SQR(d_j) + SQR(d_k);
if ( sqr_norm <= SQR(control->vlist_cut) )
{
periodic_nbrs[ *count ].d = SQRT( sqr_norm );
periodic_nbrs[ *count ].dvec[0] = d_i;
periodic_nbrs[ *count ].dvec[1] = d_j;
periodic_nbrs[ *count ].dvec[2] = d_k;
periodic_nbrs[ *count ].rel_box[0] = i;
periodic_nbrs[ *count ].rel_box[1] = j;
periodic_nbrs[ *count ].rel_box[2] = k;
/* if( i || j || k ) {
fprintf(stderr, "x1: %.2f %.2f %.2f\n", x1[0], x1[1], x1[2]);
fprintf(stderr, "x2: %.2f %.2f %.2f\n", x2[0], x2[1], x2[2]);
fprintf( stderr, "d : %8.2f%8.2f%8.2f\n\n", d_i, d_j, d_k );
} */
/* if(i) periodic_nbrs[*count].ext_factor[0] = (real)i/-abs(i);
else periodic_nbrs[*count].ext_factor[0] = 0;
if(j) periodic_nbrs[*count].ext_factor[1] = (real)j/-abs(j);
else periodic_nbrs[*count].ext_factor[1] = 0;
if(k) periodic_nbrs[*count].ext_factor[2] = (real)k/-abs(k);
else periodic_nbrs[*count].ext_factor[2] = 0; */
/* if( i == 0 && j == 0 && k == 0 )
* periodic_nbrs[ *count ].imaginary = 0;
* else periodic_nbrs[ *count ].imaginary = 1;
*/
++(*count);
}
}
}
}
/* Returns the mapping for the neighbor box pointed by (ix,iy,iz) */
/*int Get_Nbr_Box( simulation_box *box, int ix, int iy, int iz )
{
// 13 is to handle negative indexes properly
}*/
/* Returns total pressure vector for the neighbor box pointed by (ix,iy,iz) */
/*rvec Get_Nbr_Box_Press( simulation_box *box, int ix, int iy, int iz )
{
int map;
// 13 is to adjust -1,-1,-1 correspond to index 0
return box->nbr_box_press[map];
}*/
/* Increments total pressure vector for the nbr box pointed by (ix,iy,iz) */
/*void Inc_Nbr_Box_Press( simulation_box *box, int ix, int iy, int iz, rvec v )
{
int map;
rvec_Add( box->nbr_box_press[map], v );
}*/
/* Increments the total pressure vector for the neighbor box mapped to 'map' */
/*void Inc_Nbr_Box_Press( simulation_box *box, int map, rvec v )
{
rvec_Add( box->nbr_box_press[map], v );
}*/
void Print_Box_Information( simulation_box* box, FILE *out )
{
fprintf( out, "{" );
for ( j = 0; j < 3; ++j )
fprintf( out, "%8.3f ", box->box[i][j] );
fprintf( out, "}" );
fprintf( out, "V: %8.3f\tdims: {%8.3f, %8.3f, %8.3f}\n",
box->volume,
box->box_norms[0], box->box_norms[1], box->box_norms[2] );
fprintf( out, "box_trans: {" );
for ( i = 0; i < 3; ++i )
fprintf( out, "{" );
for ( j = 0; j < 3; ++j )
fprintf( out, "%8.3f ", box->trans[i][j] );
fprintf( out, "}" );
fprintf( out, "box_trinv: {" );
for ( i = 0; i < 3; ++i )
fprintf( out, "{" );
for ( j = 0; j < 3; ++j )
fprintf( out, "%8.3f ", box->trans_inv[i][j] );
fprintf( out, "}" );