jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5977 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2011-04-19 19:33:05 +00:00
parent 7d7dd36676
commit 99ed989715
3 changed files with 48 additions and 62 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
src/USER-EWALDN/ewald_n.cpp
View File

@ -353,10 +353,10 @@ void EwaldN::init_coeff_sums() // sums based on atoms
}
}
if (function[3]) { // dipole
int *type = atom->type, *ntype = type+atom->nlocal;
double *dipole = atom->dipole;
for (int *i=type; i<ntype; ++i)
sum_local[9].x2 += dipole[i[0]]*dipole[i[0]];
double **mu = atom->mu;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
sum_local[9].x2 += mu[i][3]*mu[i][3];
}
MPI_Allreduce(sum_local, sum, 2*EWALD_MAX_NSUMS, MPI_DOUBLE, MPI_SUM, world);
}
@ -413,7 +413,7 @@ void EwaldN::compute_ek()
complex *cek, zxyz, zxy, cx;
vector mui;
double *x = atom->x[0], *xn = x+3*atom->nlocal, *q = atom->q, qi, bi, ci[7];
double *dipole = atom->dipole, *mu = atom->mu ? atom->mu[0] : NULL;
double *mu = atom->mu ? atom->mu[0] : NULL;
int i, kx, ky, n = nkvec*nsums, *type = atom->type, tri = domain->triclinic;
int func[EWALD_NFUNCS];
@ -444,7 +444,8 @@ void EwaldN::compute_ek()
if (func[2]) memcpy(ci, B+7*type[0], 7*sizeof(double));
if (func[3]) {
memcpy(mui, mu, sizeof(vector)); mu += 3;
vec_scalar_mult(mui, dipole[*type]);
vec_scalar_mult(mui, mu[3]);
mu += 4;
h = hvec;
}
for (k=kvec; k<nk; ++k) { // compute rho(k)
@ -482,7 +483,7 @@ void EwaldN::compute_force()
vector sum[EWALD_MAX_NSUMS], mui;
complex *cek, zc, zx, zxy;
double *f = atom->f[0], *fn = f+3*atom->nlocal, *q = atom->q, *t = NULL;
double *dipole = atom->dipole, *mu = atom->mu ? atom->mu[0] : NULL;
double *mu = atom->mu ? atom->mu[0] : NULL;
double *ke, c[EWALD_NFUNCS] = {
8.0*M_PI*qqrd2e*scale/volume, 2.0*M_PI*sqrt(M_PI)/(12.0*volume),
2.0*M_PI*sqrt(M_PI)/(192.0*volume), 8.0*M_PI*mumurd2e/volume};
@ -500,8 +501,9 @@ void EwaldN::compute_force()
cek = cek_global;
memset(sum, 0, EWALD_MAX_NSUMS*sizeof(vector));
if (func[3]) {
register double di = dipole[*type]*c[3];
register double di = mu[3] * c[3];
mui[0] = di*(mu++)[0]; mui[1] = di*(mu++)[1]; mui[2] = di*(mu++)[2];
mu++;
}
for (nh = (h = hvec)+nkvec; h<nh; ++h, ++k) {
if (ky!=k->y) { // based on order in
@ -556,20 +558,19 @@ void EwaldN::compute_force()
}
}
void EwaldN::compute_surface()
{
if (!function[3]) return;
vector sum_local = VECTOR_NULL, sum_total;
double *mu = atom->mu ? atom->mu[0] : NULL, *dipole = atom->dipole;
int *type = atom->type, *ntype = type+atom->nlocal;
double **mu = atom->mu;
int nlocal = atom->nlocal;
for (int *i=type; i<ntype; ++i) {
register double di = dipole[i[0]];
sum_local[0] += di*(mu++)[0];
sum_local[1] += di*(mu++)[1];
sum_local[2] += di*(mu++)[2];
for (int i=0; i < nlocal; i++) {
register double di = mu[i][3];
sum_local[0] += di*mu[i][0];
sum_local[1] += di*mu[i][1];
sum_local[2] += di*mu[i][2];
}
MPI_Allreduce(sum_local, sum_total, 3, MPI_DOUBLE, MPI_SUM, world);
@ -579,7 +580,6 @@ void EwaldN::compute_surface()
energy_self[3] -= virial_self[3];
}
void EwaldN::compute_energy(int eflag)
{
energy = 0.0;

70
src/fix_rigid.cpp
View File

@ -590,7 +590,8 @@ void FixRigid::init()
// compute 6 moments of inertia of each body
// dx,dy,dz = coords relative to center-of-mass
// symmetric 3x3 inertia tensor stored in Voigt notation as 6-vector
double dx,dy,dz,rad;
for (ibody = 0; ibody < nbody; ibody++)
@ -624,16 +625,15 @@ void FixRigid::init()
sum[ibody][0] += massone * (dy*dy + dz*dz);
sum[ibody][1] += massone * (dx*dx + dz*dz);
sum[ibody][2] += massone * (dx*dx + dy*dy);
sum[ibody][3] -= massone * dx*dy;
sum[ibody][4] -= massone * dy*dz;
sum[ibody][5] -= massone * dx*dz;
sum[ibody][3] -= massone * dy*dz;
sum[ibody][4] -= massone * dx*dz;
sum[ibody][5] -= massone * dx*dy;
}
// extended particles may contribute extra terms to moments of inertia
if (extended) {
double ex[3],ey[3],ez[3],idiag[3];
double p[3][3],ptrans[3][3],ispace[3][3],itemp[3][3];
double ivec[6];
double *shape,*quatatom;
for (i = 0; i < nlocal; i++) {
@ -652,19 +652,13 @@ void FixRigid::init()
if (eflags[i] & INERTIA_ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
quatatom = ebonus[ellipsoid[i]].quat;
MathExtra::quat_to_mat(quatatom,p);
MathExtra::quat_to_mat_trans(quatatom,ptrans);
idiag[0] = 0.2*massone * (shape[1]*shape[1]+shape[2]*shape[2]);
idiag[1] = 0.2*massone * (shape[0]*shape[0]+shape[2]*shape[2]);
idiag[2] = 0.2*massone * (shape[0]*shape[0]+shape[1]*shape[1]);
MathExtra::diag_times3(idiag,ptrans,itemp);
MathExtra::times3(p,itemp,ispace);
sum[ibody][0] += ispace[0][0];
sum[ibody][1] += ispace[1][1];
sum[ibody][2] += ispace[2][2];
sum[ibody][3] += ispace[0][1];
sum[ibody][4] += ispace[1][2];
sum[ibody][5] += ispace[0][2];
MathExtra::inertia_ellipsoid(shape,quatatom,massone,ivec);
sum[ibody][0] += ivec[0];
sum[ibody][1] += ivec[1];
sum[ibody][2] += ivec[2];
sum[ibody][3] += ivec[3];
sum[ibody][4] += ivec[4];
sum[ibody][5] += ivec[5];
}
}
}
@ -674,18 +668,17 @@ void FixRigid::init()
// inertia = 3 eigenvalues = principal moments of inertia
// ex_space,ey_space,ez_space = 3 eigenvectors = principal axes of rigid body
double tensor[3][3],evectors[3][3];
int ierror;
double ez0,ez1,ez2;
double tensor[3][3],evectors[3][3];
for (ibody = 0; ibody < nbody; ibody++) {
tensor[0][0] = all[ibody][0];
tensor[1][1] = all[ibody][1];
tensor[2][2] = all[ibody][2];
tensor[0][1] = tensor[1][0] = all[ibody][3];
tensor[1][2] = tensor[2][1] = all[ibody][4];
tensor[0][2] = tensor[2][0] = all[ibody][5];
tensor[1][2] = tensor[2][1] = all[ibody][3];
tensor[0][2] = tensor[2][0] = all[ibody][4];
tensor[0][1] = tensor[1][0] = all[ibody][5];
ierror = MathExtra::jacobi(tensor,inertia[ibody],evectors);
if (ierror) error->all("Insufficient Jacobi rotations for rigid body");
@ -818,14 +811,13 @@ void FixRigid::init()
(displace[i][0]*displace[i][0] + displace[i][2]*displace[i][2]);
sum[ibody][2] += massone *
(displace[i][0]*displace[i][0] + displace[i][1]*displace[i][1]);
sum[ibody][3] -= massone * displace[i][0]*displace[i][1];
sum[ibody][4] -= massone * displace[i][1]*displace[i][2];
sum[ibody][5] -= massone * displace[i][0]*displace[i][2];
sum[ibody][3] -= massone * displace[i][1]*displace[i][2];
sum[ibody][4] -= massone * displace[i][0]*displace[i][2];
sum[ibody][5] -= massone * displace[i][0]*displace[i][1];
}
if (extended) {
double ex[3],ey[3],ez[3],idiag[3];
double p[3][3],ptrans[3][3],ispace[3][3],itemp[3][3];
double ivec[6];
double *shape;
for (i = 0; i < nlocal; i++) {
@ -843,19 +835,13 @@ void FixRigid::init()
}
if (eflags[i] & INERTIA_ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
MathExtra::quat_to_mat(qorient[i],p);
MathExtra::quat_to_mat_trans(qorient[i],ptrans);
idiag[0] = 0.2*massone * (shape[1]*shape[1]+shape[2]*shape[2]);
idiag[1] = 0.2*massone * (shape[0]*shape[0]+shape[2]*shape[2]);
idiag[2] = 0.2*massone * (shape[0]*shape[0]+shape[1]*shape[1]);
MathExtra::diag_times3(idiag,ptrans,itemp);
MathExtra::times3(p,itemp,ispace);
sum[ibody][0] += ispace[0][0];
sum[ibody][1] += ispace[1][1];
sum[ibody][2] += ispace[2][2];
sum[ibody][3] += ispace[0][1];
sum[ibody][4] += ispace[1][2];
sum[ibody][5] += ispace[0][2];
MathExtra::inertia_ellipsoid(shape,qorient[i],massone,ivec);
sum[ibody][0] += ivec[0];
sum[ibody][1] += ivec[1];
sum[ibody][2] += ivec[2];
sum[ibody][3] += ivec[3];
sum[ibody][4] += ivec[4];
sum[ibody][5] += ivec[5];
}
}
}

6
src/math_extra.cpp
View File

@ -421,9 +421,9 @@ void inertia_triangle(double *v0, double *v1, double *v2,
// NOTE: use mass
inertia[0] = inv24*a * (sum-vtsv[0][0]);
inertia[1] = inv24*a * (sum-vtsv[1][1]);
inertia[2] = inv24*a * (sum-vtsv[2][2]);
inertia[0] = inv24*a*(sum-vtsv[0][0]);
inertia[1] = inv24*a*(sum-vtsv[1][1]);
inertia[2] = inv24*a*(sum-vtsv[2][2]);
inertia[3] = -inv24*a*vtsv[1][2];
inertia[4] = -inv24*a*vtsv[0][2];
inertia[5] = -inv24*a*vtsv[0][1];