forked from lijiext/lammps
Made two changes:
-recomputed up-to-date pressure tensor (fixes energy conservation problem with aniso) -changed ndof for iso (fixes volume fluctuation problem with iso)
This commit is contained in:
parent
ad1b1897d8
commit
6f33ce9dba
|
@ -261,6 +261,60 @@ void ComputePressure::compute_vector()
|
|||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute pressure tensor
|
||||
assume KE tensor has already been computed
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputePressure::compute_vector_ke_scalar()
|
||||
{
|
||||
invoked_vector = update->ntimestep;
|
||||
if (update->vflag_global != invoked_vector)
|
||||
error->all(FLERR,"Virial was not tallied on needed timestep");
|
||||
|
||||
if (force->kspace && kspace_virial && force->kspace->scalar_pressure_flag)
|
||||
error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' for "
|
||||
"tensor components with kspace_style msm");
|
||||
|
||||
// invoke temperature if it hasn't been already
|
||||
|
||||
double t;
|
||||
if (keflag) {
|
||||
if (temperature->invoked_scalar != update->ntimestep)
|
||||
t = temperature->compute_scalar();
|
||||
else t = temperature->scalar;
|
||||
}
|
||||
|
||||
if (dimension == 3) {
|
||||
inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);
|
||||
virial_compute(6,3);
|
||||
if (keflag) {
|
||||
double kescalar = temperature->dof * boltz * t / 3.0;
|
||||
for (int i = 0; i < 3; i++)
|
||||
vector[i] = (kescalar + virial[i]) * inv_volume * nktv2p;
|
||||
for (int i = 3; i < 6; i++)
|
||||
vector[i] = virial[i] * inv_volume * nktv2p;
|
||||
} else
|
||||
for (int i = 0; i < 6; i++)
|
||||
vector[i] = virial[i] * inv_volume * nktv2p;
|
||||
} else {
|
||||
inv_volume = 1.0 / (domain->xprd * domain->yprd);
|
||||
virial_compute(4,2);
|
||||
if (keflag) {
|
||||
double kescalar = temperature->dof * boltz * t / 2.0;
|
||||
vector[0] = (kescalar + virial[0]) * inv_volume * nktv2p;
|
||||
vector[1] = (kescalar + virial[1]) * inv_volume * nktv2p;
|
||||
vector[3] = virial[3] * inv_volume * nktv2p;
|
||||
vector[2] = vector[4] = vector[5] = 0.0;
|
||||
} else {
|
||||
vector[0] = virial[0] * inv_volume * nktv2p;
|
||||
vector[1] = virial[1] * inv_volume * nktv2p;
|
||||
vector[3] = virial[3] * inv_volume * nktv2p;
|
||||
vector[2] = vector[4] = vector[5] = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputePressure::virial_compute(int n, int ndiag)
|
||||
|
|
|
@ -31,6 +31,7 @@ class ComputePressure : public Compute {
|
|||
virtual void init();
|
||||
virtual double compute_scalar();
|
||||
virtual void compute_vector();
|
||||
void compute_vector_ke_scalar();
|
||||
void reset_extra_compute_fix(const char *);
|
||||
|
||||
protected:
|
||||
|
|
|
@ -29,6 +29,7 @@
|
|||
#include "modify.h"
|
||||
#include "fix_deform.h"
|
||||
#include "compute.h"
|
||||
#include "compute_pressure.h"
|
||||
#include "kspace.h"
|
||||
#include "update.h"
|
||||
#include "respa.h"
|
||||
|
@ -777,7 +778,7 @@ void FixNH::setup(int /*vflag*/)
|
|||
|
||||
if (pstat_flag) {
|
||||
if (pstyle == ISO) pressure->compute_scalar();
|
||||
else pressure->compute_vector();
|
||||
else ((ComputePressure *)pressure)->compute_vector_ke_scalar();
|
||||
couple();
|
||||
pressure->addstep(update->ntimestep+1);
|
||||
}
|
||||
|
@ -850,7 +851,7 @@ void FixNH::initial_integrate(int /*vflag*/)
|
|||
pressure->compute_scalar();
|
||||
} else {
|
||||
temperature->compute_vector();
|
||||
pressure->compute_vector();
|
||||
((ComputePressure *)pressure)->compute_vector_ke_scalar();
|
||||
}
|
||||
couple();
|
||||
pressure->addstep(update->ntimestep+1);
|
||||
|
@ -904,9 +905,16 @@ void FixNH::final_integrate()
|
|||
t_current = temperature->compute_scalar();
|
||||
tdof = temperature->dof;
|
||||
|
||||
// need to recompute pressure to account for change in KE
|
||||
// t_current is up-to-date, but compute_temperature is not
|
||||
// compute appropriately coupled elements of mvv_current
|
||||
|
||||
if (pstat_flag) {
|
||||
if (pstyle == ISO) pressure->compute_scalar();
|
||||
else pressure->compute_vector();
|
||||
else {
|
||||
temperature->compute_vector();
|
||||
((ComputePressure *)pressure)->compute_vector_ke_scalar();
|
||||
}
|
||||
couple();
|
||||
pressure->addstep(update->ntimestep+1);
|
||||
}
|
||||
|
@ -957,7 +965,7 @@ void FixNH::initial_integrate_respa(int /*vflag*/, int ilevel, int /*iloop*/)
|
|||
pressure->compute_scalar();
|
||||
} else {
|
||||
temperature->compute_vector();
|
||||
pressure->compute_vector();
|
||||
((ComputePressure *)pressure)->compute_vector_ke_scalar();
|
||||
}
|
||||
couple();
|
||||
pressure->addstep(update->ntimestep+1);
|
||||
|
@ -1871,7 +1879,8 @@ void FixNH::nhc_press_integrate()
|
|||
}
|
||||
}
|
||||
|
||||
lkt_press = pdof * kt;
|
||||
if (pstyle == ISO) lkt_press = kt;
|
||||
else lkt_press = pdof * kt;
|
||||
etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0];
|
||||
|
||||
double ncfac = 1.0/nc_pchain;
|
||||
|
|
Loading…
Reference in New Issue