From 7cd9975e29568b3e9cb04c1660fcae76e86b1539 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Thu, 3 Feb 2022 17:44:47 -0700 Subject: [PATCH] Fixed one obvious error in numdiff version, matching analytic quite well for P~=0, need to fix addon term --- examples/cij_nvt/in.ljfcc | 11 +++- examples/cij_nvt/in.ljfcc_num | 15 ++++-- src/EXTRA-COMPUTE/compute_born_matrix.cpp | 64 +++++++++++++++-------- src/EXTRA-COMPUTE/compute_born_matrix.h | 1 + 4 files changed, 65 insertions(+), 26 deletions(-) diff --git a/examples/cij_nvt/in.ljfcc b/examples/cij_nvt/in.ljfcc index 71e541df24..08b0cd67dc 100644 --- a/examples/cij_nvt/in.ljfcc +++ b/examples/cij_nvt/in.ljfcc @@ -1,18 +1,25 @@ # Numerical difference calculation # of Born matrix +# Note that because of cubic symmetry and central forces, we have: +# C11, pure axial == positive mean value: 1,2,3 +# C44==C23, pure shear == positive mean value, exactly match in pairs: (4,12),(5,8),(6,7) +# C14==C56, shear/axial(normal) == zero mean, exactly match in pairs: (9,21),(14,20),(18,19) +# C15, shear/axial(in-plane) == zero mean: 10,11,13,15,16,17 + # adjustable parameters -variable nsteps index 50 # length of run +variable nsteps index 50 # length of run variable nthermo index 10 # thermo output interval variable nlat index 3 # size of box variable delta index 1.0e-6 # strain size variable temp index 0.3 # temperature +variable rho index 0.8442 # reduced density units lj atom_style atomic -lattice fcc 0.8442 +lattice fcc ${rho} region box block 0 ${nlat} 0 ${nlat} 0 ${nlat} create_box 1 box create_atoms 1 box diff --git a/examples/cij_nvt/in.ljfcc_num b/examples/cij_nvt/in.ljfcc_num index 7ef456150e..af5fa02c6e 100644 --- a/examples/cij_nvt/in.ljfcc_num +++ b/examples/cij_nvt/in.ljfcc_num @@ -1,18 +1,27 @@ # Numerical difference calculation # of Born matrix +# Note that because of cubic symmetry and central forces, we have: +# C11, pure axial == positive mean value: 1,2,3 +# C44==C23, pure shear == positive mean value, (exactly) match in pairs: (4,12),(5,8),(6,7) +# C14==C56, shear/axial(normal) == zero mean, (exactly) match in pairs: (9,21),(14,20),(18,19) +# C15, shear/axial(in-plane) == zero mean: 10,11,13,15,16,17 +# the "(exactly)" above is because these symmetries seem to be sensitive to details +# of how the finite shear deformations are applied + # adjustable parameters -variable nsteps index 50 # length of run +variable nsteps index 50 # length of run variable nthermo index 10 # thermo output interval variable nlat index 3 # size of box -variable delta index 1.0e-4 # strain size +variable delta index 1.0e-6 # strain size variable temp index 0.3 # temperature +variable rho index 0.8442 # reduced density units lj atom_style atomic -lattice fcc 0.8442 +lattice fcc ${rho} region box block 0 ${nlat} 0 ${nlat} 0 ${nlat} create_box 1 box create_atoms 1 box diff --git a/src/EXTRA-COMPUTE/compute_born_matrix.cpp b/src/EXTRA-COMPUTE/compute_born_matrix.cpp index cf2bd95f25..e943b1c3e8 100644 --- a/src/EXTRA-COMPUTE/compute_born_matrix.cpp +++ b/src/EXTRA-COMPUTE/compute_born_matrix.cpp @@ -287,7 +287,8 @@ void ComputeBornMatrix::init() // to convert to vector indexing: // matrix index to vector index, double -> single index - + // this is not used at all + voigt3MtoV[0][0]=0; voigt3MtoV[0][1]=5; voigt3MtoV[0][2]=4; voigt3MtoV[1][0]=5; voigt3MtoV[1][1]=1; voigt3MtoV[1][2]=3; voigt3MtoV[2][0]=4; voigt3MtoV[2][1]=3; voigt3MtoV[2][2]=2; @@ -299,8 +300,19 @@ void ComputeBornMatrix::init() virialMtoV[0][0]=0; virialMtoV[0][1]=3; virialMtoV[0][2]=4; virialMtoV[1][0]=3; virialMtoV[1][1]=1; virialMtoV[1][2]=5; virialMtoV[2][0]=4; virialMtoV[2][1]=5; virialMtoV[2][2]=2; - - // the following is for 6x6 matrix and vector indexing converter + + // reorder LAMMPS virial vector to Voigt order + + virialVtoV[0] = 0; + virialVtoV[1] = 1; + virialVtoV[2] = 2; + virialVtoV[3] = 5; + virialVtoV[4] = 4; + virialVtoV[5] = 3; + + // the following is for 6x6 matrix and vector indexing converter // this is clearly different order form albe[][] and revalbe[] + // should not be used + int indcounter = 0; for(int row = 0; row < NDIR_VIRIAL; row++) for(int col = row; col< NDIR_VIRIAL; col++) { @@ -351,6 +363,13 @@ void ComputeBornMatrix::compute_vector() MPI_Allreduce(values_local, values_global, nvalues, MPI_DOUBLE, MPI_SUM, world); + // convert to pressure units + + double nktv2p = force->nktv2p; + double inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd); + for (int m = 0; m < nvalues; m++) { values_global[m] *= + (nktv2p * inv_volume); } + } else { // calculate Born matrix using stress finite differences @@ -358,9 +377,9 @@ void ComputeBornMatrix::compute_vector() compute_numdiff(); } - - for (int m = 0; m < nvalues; m++) { vector[m] = values_global[m]; } - + + for (int m = 0; m < nvalues; m++) vector[m] = values_global[m]; + } /* ---------------------------------------------------------------------- @@ -443,7 +462,8 @@ void ComputeBornMatrix::compute_pairs() // Add contribution to Born tensor pair->born_matrix(i, j, itype, jtype, rsq, factor_coul, factor_lj, dupair, du2pair); - pair_pref = du2pair - dupair * rinv; + // pair_pref = du2pair - dupair * rinv; + pair_pref = du2pair; // See albemunu in compute_born_matrix.h for indices order. @@ -489,7 +509,7 @@ void ComputeBornMatrix::compute_numdiff() } // loop over 6 strain directions - // compute a finite difference force in each dimension + // compute stress finite difference in each direction int flag, allflag; @@ -499,26 +519,24 @@ void ComputeBornMatrix::compute_numdiff() update_virial(); for (int jdir = 0; jdir < NDIR_VIRIAL; jdir++) { m = revalbe[idir][jdir]; - values_global[m] = compute_virial->vector[jdir]; + values_global[m] = compute_virial->vector[virialVtoV[jdir]]; } - restore_atoms(nall, idir); + displace_atoms(nall, idir, -1.0); force_clear(nall); update_virial(); for (int jdir = 0; jdir < NDIR_VIRIAL; jdir++) { m = revalbe[idir][jdir]; - values_global[m] -= compute_virial->vector[jdir]; + values_global[m] -= compute_virial->vector[virialVtoV[jdir]]; } restore_atoms(nall, idir); } // apply derivative factor - double nktv2p = force->nktv2p; - double inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd); - // double denominator = -0.5 / numdelta * inv_volume * nktv2p; double denominator = -0.5 / numdelta; + for (int m = 0; m < nvalues; m++) values_global[m] *= denominator; // recompute virial so all virial and energy contributions are as before // also needed for virial stress addon contributions to Born matrix @@ -526,15 +544,14 @@ void ComputeBornMatrix::compute_numdiff() update_virial(); + // virial_addon(); + // restore original forces for owned and ghost atoms for (int i = 0; i < nall; i++) for (int k = 0; k < 3; k++) f[i][k] = temp_f[i][k]; - for (int m = 0; m < nvalues; m++) values_global[m] *= denominator; - - virial_addon(); } /* ---------------------------------------------------------------------- @@ -544,6 +561,7 @@ void ComputeBornMatrix::compute_numdiff() void ComputeBornMatrix::displace_atoms(int nall, int idir, double magnitude) { double **x = atom->x; + int k = dirlist[idir][0]; int l = dirlist[idir][1]; for (int i = 0; i < nall; i++) @@ -556,9 +574,13 @@ void ComputeBornMatrix::displace_atoms(int nall, int idir, double magnitude) void ComputeBornMatrix::restore_atoms(int nall, int idir) { + + // reset all coords, just to be safe, ignore idir + double **x = atom->x; - int k = dirlist[idir][0]; - for (int i = 0; i < nall; i++) { x[i][k] = temp_x[i][k]; } + for (int i = 0; i < nall; i++) + for (int k = 0; k < 3; k++) + x[i][k] = temp_x[i][k]; } /* ---------------------------------------------------------------------- @@ -585,7 +607,6 @@ void ComputeBornMatrix::update_virial() compute_virial->compute_vector(); } - /* ---------------------------------------------------------------------- calculate virial stress addon terms to the Born matrix this is based on original code of Dr. Yubao Zhen @@ -630,7 +651,8 @@ void ComputeBornMatrix::virial_addon() if(SHEAR) addon += kronecker[jd][nd]*sigv[virialMtoV[id][kd]] + kronecker[jd][kd]*sigv[virialMtoV[id][nd]]; - values_global[voigt6MtoV[ijvgt][knvgt]] += addon; + // values_global[voigt6MtoV[ijvgt][knvgt]] += addon; + values_global[revalbe[ijvgt][knvgt]] += addon; } } } diff --git a/src/EXTRA-COMPUTE/compute_born_matrix.h b/src/EXTRA-COMPUTE/compute_born_matrix.h index bdef7ef86e..ffd180d094 100644 --- a/src/EXTRA-COMPUTE/compute_born_matrix.h +++ b/src/EXTRA-COMPUTE/compute_born_matrix.h @@ -75,6 +75,7 @@ class ComputeBornMatrix : public Compute { int voigt3VtoM[NDIR_VIRIAL][2]; int voigt3MtoV[NXYZ_VIRIAL][NXYZ_VIRIAL]; int virialMtoV[NXYZ_VIRIAL][NXYZ_VIRIAL]; + int virialVtoV[NDIR_VIRIAL]; int voigt6MtoV[NDIR_VIRIAL][NDIR_VIRIAL]; int kronecker[NXYZ_VIRIAL][NXYZ_VIRIAL]; double **temp_x; // original coords