ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added correct units for pressure, still getting wrong answer

Browse Source
This commit is contained in:
Aidan Thompson
2022-01-21 19:46:54 -07:00
parent 3286fcb13a
commit 3e3cd4c94d
2 changed files with 10 additions and 16 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
examples/numdiff/in.numdiff
View File

@ -20,7 +20,7 @@ neighbor 1.0 bin
timestep 0.001 timestep 0.001
fix numforce all numdiff 100 0.0001 fix numforce all numdiff 100 0.0001
fix numstress all numdiff/stress 100 0.0001 fix numstress all numdiff/stress 100 1.0e-2
fix nve all nve fix nve all nve
variable errx atom fx-f_numforce[1] variable errx atom fx-f_numforce[1]

24
src/EXTRA-FIX/fix_numdiff_stress.cpp
View File

@ -70,8 +70,6 @@ FixNumDiffStress::FixNumDiffStress(LAMMPS *lmp, int narg, char **arg) :
// zero numdiff_forces since a variable may access it before first run // zero numdiff_forces since a variable may access it before first run
reallocate(); reallocate();
stress_clear();
// set fixed-point to default = center of cell // set fixed-point to default = center of cell
@ -204,15 +202,11 @@ void FixNumDiffStress::calculate_stress()
int nall = nlocal + atom->nghost; int nall = nlocal + atom->nghost;
for (int i = 0; i < nall; i++) for (int i = 0; i < nall; i++)
for (int j = 0; j < 3; j++) { for (int k = 0; k < 3; k++) {
temp_x[i][j] = x[i][j]; temp_x[i][k] = x[i][k];
temp_f[i][j] = f[i][j]; temp_f[i][k] = f[i][k];
} }
// initialize numerical stress to zero
stress_clear();
// loop over 6 strain directions // loop over 6 strain directions
// compute a finite difference force in each dimension // compute a finite difference force in each dimension
@ -220,12 +214,12 @@ void FixNumDiffStress::calculate_stress()
double nktv2p = force->nktv2p; double nktv2p = force->nktv2p;
double inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd); double inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);
double denominator = 0.5 / delta * inv_volume * nktv2p; double denominator = -0.5 / delta * inv_volume * nktv2p;
for (int idir = 0; idir < NDIR_STRESS; idir++) { for (int idir = 0; idir < NDIR_STRESS; idir++) {
displace_atoms(nall, idir, 1.0); displace_atoms(nall, idir, 1.0);
energy = update_energy(); energy = update_energy();
stress[idir] += energy; stress[idir] = energy;
displace_atoms(nall, idir,-1.0); displace_atoms(nall, idir,-1.0);
energy = update_energy(); energy = update_energy();
stress[idir] -= energy; stress[idir] -= energy;
@ -236,8 +230,8 @@ void FixNumDiffStress::calculate_stress()
// restore original forces for owned and ghost atoms // restore original forces for owned and ghost atoms
for (int i = 0; i < nall; i++) for (int i = 0; i < nall; i++)
for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++)
f[i][j] = temp_f[i][j]; f[i][k] = temp_f[i][k];
} }
@ -262,9 +256,9 @@ void FixNumDiffStress::displace_atoms(int nall, int idir, double magnitude)
void FixNumDiffStress::restore_atoms(int nall, int idir) void FixNumDiffStress::restore_atoms(int nall, int idir)
{ {
double **x = atom->x; double **x = atom->x;
int j = dirlist[idir][0]; int k = dirlist[idir][0];
for (int i = 0; i < nall; i++) { for (int i = 0; i < nall; i++) {
x[i][j] = temp_x[i][j]; x[i][k] = temp_x[i][k];
} }
} }