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:29:23 -07:00
parent c16caea13b
commit 3286fcb13a
2 changed files with 6 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
examples/numdiff/in.numdiff
View File

@ -31,9 +31,8 @@ dump errors all custom 100 force_error.dump v_errx v_erry v_errz
variable ferrsq atom (fx-f_numforce[1])^2+(fy-f_numforce[2])^2+(fz-f_numforce[3])^2 variable ferrsq atom (fx-f_numforce[1])^2+(fy-f_numforce[2])^2+(fz-f_numforce[3])^2
compute faverrsq all reduce ave v_ferrsq compute faverrsq all reduce ave v_ferrsq
fix avfaverrsq all ave/time 100 1 100 c_faverrsq ave running
compute myvirial all pressure NULL virial compute myvirial all pressure NULL virial
variable ratio11 equal f_numstress[1]/c_myvirial[1]
thermo_style custom step temp pe press c_faverrsq f_avfaverrsq c_myvirial[1] f_numstress[1] thermo_style custom step temp pe press c_faverrsq c_myvirial[*] f_numstress[*] v_ratio11
thermo 100 thermo 100
run 500 run 500

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

@ -217,7 +217,10 @@ void FixNumDiffStress::calculate_stress()
// compute a finite difference force in each dimension // compute a finite difference force in each dimension
int flag,allflag; int flag,allflag;
double denominator = 0.5 / delta; double nktv2p = force->nktv2p;
double inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);
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);