lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

Final changes to fix_box_relax and fix_nh 

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@3961 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
athomps 2010-04-07 18:34:00 +00:00
parent 28bbe502e5
commit 20cfa91d13
3 changed files with 57 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
src/fix_box_relax.cpp
View File

@ -44,7 +44,7 @@ enum{ISO,ANISO,TRICLINIC};
FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 4) error->all("Illegal fix box/relax command");
if (narg < 5) error->all("Illegal fix box/relax command");
scalar_flag = 1;
extscalar = 1;
@ -62,6 +62,8 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) :
p_target[0] = p_target[1] = p_target[2] =
p_target[3] = p_target[4] = p_target[5] = 0.0;
p_flag[0] = p_flag[1] = p_flag[2] =
p_flag[3] = p_flag[4] = p_flag[5] = 0;
// process keywords
@ -71,7 +73,7 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) :
while (iarg < narg) {
if (strcmp(arg[iarg],"iso") == 0) {
if (iarg+4 > narg) error->all("Illegal fix box/relax command");
if (iarg+2 > narg) error->all("Illegal fix box/relax command");
pcouple = XYZ;
p_target[0] = p_target[1] = p_target[2] = atof(arg[iarg+1]);
p_flag[0] = p_flag[1] = p_flag[2] = 1;
@ -79,9 +81,9 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) :
p_target[2] = 0.0;
p_flag[2] = 0;
}
iarg += 4;
iarg += 2;
} else if (strcmp(arg[iarg],"aniso") == 0) {
if (iarg+4 > narg) error->all("Illegal fix box/relax command");
if (iarg+2 > narg) error->all("Illegal fix box/relax command");
pcouple = NONE;
p_target[0] = p_target[1] = p_target[2] = atof(arg[iarg+1]);
p_flag[0] = p_flag[1] = p_flag[2] = 1;
@ -89,9 +91,9 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) :
p_target[2] = 0.0;
p_flag[2] = 0;
}
iarg += 4;
iarg += 2;
} else if (strcmp(arg[iarg],"tri") == 0) {
if (iarg+4 > narg) error->all("Illegal fix box/relax command");
if (iarg+2 > narg) error->all("Illegal fix box/relax command");
pcouple = NONE;
p_target[0] = p_target[1] = p_target[2] = atof(arg[iarg+1]);
p_flag[0] = p_flag[1] = p_flag[2] = 1;
@ -101,51 +103,51 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) :
p_target[2] = p_target[3] = p_target[4] = 0.0;
p_flag[2] = p_flag[3] = p_flag[4] = 0;
}
iarg += 4;
iarg += 2;
} else if (strcmp(arg[iarg],"x") == 0) {
if (iarg+4 > narg) error->all("Illegal fix box/relax command");
if (iarg+2 > narg) error->all("Illegal fix box/relax command");
p_target[0] = atof(arg[iarg+1]);
p_flag[0] = 1;
deviatoric_flag = 1;
iarg += 4;
iarg += 2;
} else if (strcmp(arg[iarg],"y") == 0) {
if (iarg+4 > narg) error->all("Illegal fix box/relax command");
if (iarg+2 > narg) error->all("Illegal fix box/relax command");
p_target[1] = atof(arg[iarg+1]);
p_flag[1] = 1;
deviatoric_flag = 1;
iarg += 4;
iarg += 2;
} else if (strcmp(arg[iarg],"z") == 0) {
if (iarg+4 > narg) error->all("Illegal fix box/relax command");
if (iarg+2 > narg) error->all("Illegal fix box/relax command");
p_target[2] = atof(arg[iarg+1]);
p_flag[2] = 1;
deviatoric_flag = 1;
iarg += 4;
iarg += 2;
if (dimension == 2)
error->all("Invalid fix box/relax command for a 2d simulation");
} else if (strcmp(arg[iarg],"yz") == 0) {
if (iarg+4 > narg) error->all("Illegal fix box/relax command");
if (iarg+2 > narg) error->all("Illegal fix box/relax command");
p_target[3] = atof(arg[iarg+1]);
p_flag[3] = 1;
deviatoric_flag = 1;
iarg += 4;
iarg += 2;
if (dimension == 2)
error->all("Invalid fix box/relax command for a 2d simulation");
} else if (strcmp(arg[iarg],"xz") == 0) {
if (iarg+4 > narg) error->all("Illegal fix box/relax command");
if (iarg+2 > narg) error->all("Illegal fix box/relax command");
p_target[4] = atof(arg[iarg+1]);
p_flag[4] = 1;
deviatoric_flag = 1;
iarg += 4;
iarg += 2;
if (dimension == 2)
error->all("Invalid fix box/relax command for a 2d simulation");
} else if (strcmp(arg[iarg],"xy") == 0) {
if (iarg+4 > narg) error->all("Illegal fix box/relax command");
if (iarg+2 > narg) error->all("Illegal fix box/relax command");
p_target[5] = atof(arg[iarg+1]);
p_flag[5] = 1;
deviatoric_flag = 1;
iarg += 4;
iarg += 2;
} else if (strcmp(arg[iarg],"couple") == 0) {
if (iarg+2 > narg) error->all("Illegal fix box/relax command");

59
src/fix_nh.cpp
View File

@ -493,6 +493,7 @@ void FixNH::init()
dthalf = 0.5 * update->dt;
dt4 = 0.25 * update->dt;
dt8 = 0.125 * update->dt;
dto = dthalf;
p_freq_max = 0.0;
if (pstat_flag) {
@ -503,11 +504,11 @@ void FixNH::init()
p_freq_max = MAX(p_freq_max,p_freq[4]);
p_freq_max = MAX(p_freq_max,p_freq[5]);
}
drag_factor = 1.0 - (update->dt * p_freq_max * drag);
pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain);
}
if (tstat_flag)
drag_factor = 1.0 - (update->dt * MAX(p_freq_max,t_freq) * drag);
tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain);
// tally the number of dimensions that are barostatted
// also compute the initial volume and reference cell
@ -536,6 +537,7 @@ void FixNH::init()
if (strcmp(update->integrate_style,"respa") == 0) {
nlevels_respa = ((Respa *) update->integrate)->nlevels;
step_respa = ((Respa *) update->integrate)->step;
dto = 0.5*step_respa[0];
}
// detect if any rigid fixes exist so rigid bodies move when box is remapped
@ -709,6 +711,7 @@ void FixNH::initial_integrate(int vflag)
void FixNH::final_integrate()
{
nve_v();
if (pstat_flag) nh_v_press();
@ -779,14 +782,14 @@ void FixNH::initial_integrate_respa(int vflag, int ilevel, int iloop)
temperature->compute_scalar();
double tmp = pressure->compute_scalar();
} else {
temperature->compute_vector();
temperature->compute_vector();
pressure->compute_vector();
}
couple();
pressure->addstep(update->ntimestep+1);
if (mtk_flag) couple_ke();
}
if (pstat_flag) {
compute_press_target();
nh_omega_dot();
@ -795,24 +798,22 @@ void FixNH::initial_integrate_respa(int vflag, int ilevel, int iloop)
nve_v();
// remap simulation box by 1/2 step
if (pstat_flag) remap();
} else nve_v();
// innermost level - also update x only for atoms in group
// if iloop = 0, then is 1st call at innermost level from rRESPA
// perform 2nd half of box remap
// redo KSpace coeffs since volume has changed
// if barostat, perform 1/2 step remap before and after
if (ilevel == 0) {
if (pstat_flag) remap();
nve_x();
if (iloop == 0) {
remap();
if (kspace_flag) force->kspace->setup();
}
if (pstat_flag) remap();
}
// if barostat, redo KSpace coeffs at outermost level,
// since volume has changed
if (ilevel == nlevels_respa-1 && kspace_flag && pstat_flag) force->kspace->setup();
}
/* ---------------------------------------------------------------------- */
@ -1339,7 +1340,7 @@ double FixNH::compute_vector(int n)
if (pstyle == ISO) {
ilen = 1;
if (n < ilen)
return 0.5*omega_dot[n]*omega_dot[n]*omega_mass[n];
return pdim*0.5*omega_dot[n]*omega_dot[n]*omega_mass[n];
n -= ilen;
} else if (pstyle == ANISO) {
ilen = 3;
@ -1396,6 +1397,12 @@ void FixNH::reset_dt()
dthalf = 0.5 * update->dt;
dt4 = 0.25 * update->dt;
dt8 = 0.125 * update->dt;
dto = dthalf;
// If using respa, then remap is performed in innermost level
if (strcmp(update->integrate_style,"respa") == 0)
dto = 0.5*step_respa[0];
p_freq_max = 0.0;
if (pstat_flag) {
@ -1406,11 +1413,11 @@ void FixNH::reset_dt()
p_freq_max = MAX(p_freq_max,p_freq[4]);
p_freq_max = MAX(p_freq_max,p_freq[5]);
}
drag_factor = 1.0 - (update->dt * p_freq_max * drag);
pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain);
}
if (tstat_flag)
drag_factor = 1.0 - (update->dt * MAX(p_freq_max,t_freq) * drag);
tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain);
}
/* ----------------------------------------------------------------------
@ -1433,14 +1440,14 @@ void FixNH::nhc_temp_integrate()
expfac = exp(-ncfac*dt8*eta_dot[ich+1]);
eta_dot[ich] *= expfac;
eta_dot[ich] += eta_dotdot[ich] * ncfac*dt4;
eta_dot[ich] *= drag_factor;
eta_dot[ich] *= tdrag_factor;
eta_dot[ich] *= expfac;
}
expfac = exp(-ncfac*dt8*eta_dot[1]);
eta_dot[0] *= expfac;
eta_dot[0] += eta_dotdot[0] * ncfac*dt4;
eta_dot[0] *= drag_factor;
eta_dot[0] *= tdrag_factor;
eta_dot[0] *= expfac;
factor_eta = exp(-ncfac*dthalf*eta_dot[0]);
@ -1501,14 +1508,14 @@ void FixNH::nhc_press_integrate()
expfac = exp(-ncfac*dt8*etap_dot[ich+1]);
etap_dot[ich] *= expfac;
etap_dot[ich] += etap_dotdot[ich] * ncfac*dt4;
etap_dot[ich] *= drag_factor;
etap_dot[ich] *= pdrag_factor;
etap_dot[ich] *= expfac;
}
expfac = exp(-ncfac*dt8*etap_dot[1]);
etap_dot[0] *= expfac;
etap_dot[0] += etap_dotdot[0] * ncfac*dt4;
etap_dot[0] *= drag_factor;
etap_dot[0] *= pdrag_factor;
etap_dot[0] *= expfac;
for (ich = 0; ich < mpchain; ich++)
@ -1854,11 +1861,11 @@ void FixNH::nh_omega_dot()
(omega_mass[i] * nktv2p);
if (deviatoric_flag) f_omega -= fdev[i]/(omega_mass[i] * nktv2p);
omega_dot[i] += f_omega*dthalf;
omega_dot[i] *= drag_factor;
omega_dot[i] *= pdrag_factor;
}
omega[i] += dthalf*omega_dot[i];
factor[i] = exp(-dthalf*omega_dot[i]*mtk_factor);
dilation[i] = exp(dthalf*omega_dot[i]);
dilation[i] = exp(dto*omega_dot[i]);
}
if (pstyle == TRICLINIC) {
@ -1868,11 +1875,11 @@ void FixNH::nh_omega_dot()
if (deviatoric_flag)
f_omega -= fdev[i]/(omega_mass[i] * nktv2p);
omega_dot[i] += f_omega*dthalf;
omega_dot[i] *= drag_factor;
omega_dot[i] *= pdrag_factor;
}
omega[i] += dthalf*omega_dot[i];
factor[i] = -dthalf*omega_dot[i];
dilation[i] = dthalf*omega_dot[i];
dilation[i] = dto*omega_dot[i];
}
}
}

5
src/fix_nh.h
View File

@ -38,7 +38,7 @@ class FixNH : public Fix {
protected:
int dimension,which;
double dtv,dtf,dthalf,dt4,dt8;
double dtv,dtf,dthalf,dt4,dt8,dto;
double boltz,nktv2p,tdof;
double vol0,t0;
@ -56,7 +56,8 @@ class FixNH : public Fix {
double omega[6],omega_dot[6];
double omega_mass[6];
double p_current[6],dilation[6];
double drag,drag_factor; // drag factor on particle thermostat
double drag,tdrag_factor; // drag factor on particle thermostat
double pdrag_factor; // drag factor on barostat
double factor[6]; // velocity scaling due to barostat
int kspace_flag; // 1 if KSpace invoked, 0 if not
int nrigid; // number of rigid fixes