/*----------------------------------------------------------------------
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
the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
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 )
{
double c_alpha, c_beta, c_gamma, s_gamma, zi;
c_alpha = cos(DEG2RAD(alpha));
c_beta = cos(DEG2RAD(beta));
c_gamma = cos(DEG2RAD(gamma));
s_gamma = sin(DEG2RAD(gamma));
zi = (c_alpha - c_beta * c_gamma) / s_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 );
#if defined(DEBUG_FOCUS)
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 )
{
int i, j;
for (i = 0; i < 3; i++)
for (j = 0; j < 3; j++)
box->box[i][j] = box_tensor[i][j];
Make_Consistent( box );
}
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)
{
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,"%lf\t",box->trans[i][j]);
// fprintf(stderr,"\n");
// }
// fprintf(stderr,"\n");
// for (i=0; i < 3; i++)
// {
// for (j=0; j < 3; j++)
// fprintf(stderr,"%lf\t",box->trans_inv[i][j]);
// fprintf(stderr,"\n");
// }
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;
}
}
// printf(">x2: (%lf, %lf, %lf)\n", x2[0], x2[1], x2[2]);
}
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 )
{
int i;
real tmp;
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;
}
}
return norm;
}
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]);
Sq_Distance_on_T3( xa, xb, box, ra );
Transform( ra, box, 1, r );
}
void Inc_on_T3_Gen( rvec x, rvec dx, simulation_box* box )
{
rvec xa, dxa;
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]);
Inc_on_T3( xa, dxa, box );
//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 )
{
int i, j;
real dist = 0.0, tmp;
for ( i = 0; i < 3; i++ )
{
tmp = 0.0;
for ( j = 0; j < 3; j++ )
tmp += box->g[i][j] * x2[j];
dist += x1[i] * tmp;
}
return dist;
}
int Are_Far_Neighbors( rvec x1, rvec x2, simulation_box *box,
real cutoff, far_neighbor_data *data )
{
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 )
{
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 );
}
else *count = 0;
}
/* Finds all periodic far neighborhoods between x1 and x2
((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 )
{
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 )
{
return (9 * ix + 3 * iy + iz + 13);
// 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;
map = 9 * ix + 3 * iy + iz + 13;
// 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;
map = 9 * ix + 3 * iy + iz + 13;
// 13 is to adjust -1,-1,-1 correspond to index 0
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 )
{
int i, j;
fprintf( out, "box: {" );
for ( i = 0; i < 3; ++i )
{
fprintf( out, "{" );
for ( j = 0; j < 3; ++j )
fprintf( out, "%8.3f ", box->box[i][j] );
fprintf( out, "}" );
}
fprintf( out, "}\n" );
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, "}\n" );
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, "}" );
}
fprintf( out, "}\n" );
}