diff --git a/src/COLLOID/pair_lubricate.cpp b/src/COLLOID/pair_lubricate.cpp index 9bfe41f487..7b4aac5a9f 100644 --- a/src/COLLOID/pair_lubricate.cpp +++ b/src/COLLOID/pair_lubricate.cpp @@ -28,8 +28,8 @@ #include "neigh_request.h" #include "update.h" #include "memory.h" +#include "random_mars.h" #include "error.h" -#include "domain.h" using namespace LAMMPS_NS; @@ -41,6 +41,8 @@ using namespace LAMMPS_NS; PairLubricate::PairLubricate(LAMMPS *lmp) : Pair(lmp) { single_enable = 0; + + random = NULL; } /* ---------------------------------------------------------------------- */ @@ -54,6 +56,7 @@ PairLubricate::~PairLubricate() memory->destroy_2d_double_array(cut); memory->destroy_2d_double_array(cut_inner); } + delete random; } /* ---------------------------------------------------------------------- */ @@ -61,8 +64,8 @@ PairLubricate::~PairLubricate() void PairLubricate::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,fpair,fx,fy,fz; - double rsq,r,h_sep,radi; + double xtmp,ytmp,ztmp,delx,dely,delz,fpair,fx,fy,fz,tx,ty,tz; + double rsq,r,h_sep,radi,tfmag; double vr1,vr2,vr3,vnnr,vn1,vn2,vn3; double vt1,vt2,vt3,w1,w2,w3,v_shear1,v_shear2,v_shear3; double omega_t_1,omega_t_2,omega_t_3; @@ -88,6 +91,9 @@ void PairLubricate::compute(int eflag, int vflag) int omega_flag = atom->omega_flag; double vxmu2f = force->vxmu2f; + double prethermostat = sqrt(2.0 * force->boltz * t_target / update->dt); + prethermostat *= sqrt(force->vxmu2f/force->ftm2v/force->mvv2e); + inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; @@ -167,7 +173,7 @@ void PairLubricate::compute(int eflag, int vflag) omega_t_2 = w2 - dely*(dely*w2) / rsq; omega_t_3 = w3 - delz*(delz*w3) / rsq; - // n x omega_t + // n X omega_t n_cross_omega_t_1 = (dely*omega_t_3 - delz*omega_t_2) / r; n_cross_omega_t_2 = -(delx*omega_t_3 - delz*omega_t_1) / r; @@ -186,9 +192,9 @@ void PairLubricate::compute(int eflag, int vflag) } wnnr = wr1*delx + wr2*dely + wr3*delz; - wn1 = delx*wnnr / rsq; - wn2 = dely*wnnr / rsq; - wn3 = delz*wnnr / rsq; + wn1 = delx*wnnr / rsq; + wn2 = dely*wnnr / rsq; + wn3 = delz*wnnr / rsq; P_dot_wrel_1 = wr1 - delx*(delx*wr1)/rsq; P_dot_wrel_2 = wr2 - dely*(dely*wr2)/rsq; @@ -199,17 +205,16 @@ void PairLubricate::compute(int eflag, int vflag) h_sep = r - 2.0*radi; if (flag1) - a_squeeze = vxmu2f * - (3.0*PI*mu*2.0*radi/2.0) * (2.0*radi/4.0/h_sep); + a_squeeze = (3.0*PI*mu*2.0*radi/2.0) * (2.0*radi/4.0/h_sep); if (flag2) - a_shear = vxmu2f * (PI*mu*2.*radi/2.0) * + a_shear = (PI*mu*2.*radi/2.0) * log(2.0*radi/2.0/h_sep)*(2.0*radi+h_sep)*(2.0*radi+h_sep)/4.0; if (flag3) - a_pump = vxmu2f * (PI*mu*pow(2.0*radi,4)/8.0) * + a_pump = (PI*mu*pow(2.0*radi,4)/8.0) * ((3.0/20.0) * log(2.0*radi/2.0/h_sep) + (63.0/250.0) * (h_sep/2.0/radi) * log(2.0*radi/2.0/h_sep)); if (flag4) - a_twist = vxmu2f * (PI*mu*pow(2.0*radi,4)/4.0) * + a_twist = (PI*mu*pow(2.0*radi,4)/4.0) * (h_sep/2.0/radi) * log(2.0/(2.0*h_sep)); if (h_sep >= cut_inner[itype][jtype]) { @@ -219,19 +224,40 @@ void PairLubricate::compute(int eflag, int vflag) (2.0/r)*a_shear*n_cross_omega_t_2; fz = -a_squeeze*vn3 - a_shear*(2.0/r)*(2.0/r)*vt3 + (2.0/r)*a_shear*n_cross_omega_t_3; + fx *= vxmu2f; + fy *= vxmu2f; + fz *= vxmu2f; + + // add in thermostat force + + tfmag = prethermostat*sqrt(a_squeeze)*(random->uniform()-0.5); + fx -= tfmag * delx/r; + fy -= tfmag * dely/r; + fz -= tfmag * delz/r; - torque[i][0] += -(2.0/r)*a_shear*v_shear1 - a_shear*omega_t_1 - + tx = -(2.0/r)*a_shear*v_shear1 - a_shear*omega_t_1 - a_pump*P_dot_wrel_1 - a_twist*wn1; - torque[i][1] += -(2.0/r)*a_shear*v_shear2 - a_shear*omega_t_2 - + ty = -(2.0/r)*a_shear*v_shear2 - a_shear*omega_t_2 - a_pump*P_dot_wrel_2 - a_twist*wn2; - torque[i][2] += -(2.0/r)*a_shear*v_shear3 - a_shear*omega_t_3 - + tz = -(2.0/r)*a_shear*v_shear3 - a_shear*omega_t_3 - a_pump*P_dot_wrel_3 - a_twist*wn3; + torque[i][0] += vxmu2f * tx; + torque[i][1] += vxmu2f * ty; + torque[i][2] += vxmu2f * tz; } else { - fpair = -vnnr*(3.0*PI*mu*radi/2.0)*radi/4.0/cut_inner[itype][jtype]; - fx = fpair*delx/r; - fy = fpair*dely/r; - fz = fpair*delz/r; + a_squeeze = (3.0*PI*mu*2.0*radi/2.0) * + (2.0*radi/4.0/cut_inner[itype][jtype]); + fpair = -a_squeeze*vnnr; + fpair *= vxmu2f; + + // add in thermostat force + + fpair -= prethermostat*sqrt(a_squeeze)*(random->uniform()-0.5); + + fx = fpair * delx/r; + fy = fpair * dely/r; + fz = fpair * delz/r; } f[i][0] += fx; @@ -244,12 +270,15 @@ void PairLubricate::compute(int eflag, int vflag) f[j][2] -= fz; if (h_sep >= cut_inner[itype][jtype]) { - torque[j][0] += -(2.0/r)*a_shear*v_shear1 - a_shear*omega_t_1 + + tx = -(2.0/r)*a_shear*v_shear1 - a_shear*omega_t_1 + a_pump*P_dot_wrel_1 + a_twist*wn1; - torque[j][1] += -(2.0/r)*a_shear*v_shear2 - a_shear*omega_t_2 + + ty = -(2.0/r)*a_shear*v_shear2 - a_shear*omega_t_2 + a_pump*P_dot_wrel_2 + a_twist*wn2; - torque[j][2] += -(2.0/r)*a_shear*v_shear3 - a_shear*omega_t_3 + + tz = -(2.0/r)*a_shear*v_shear3 - a_shear*omega_t_3 + a_pump*P_dot_wrel_3 + a_twist*wn3; + torque[j][0] += vxmu2f * tx; + torque[j][1] += vxmu2f * ty; + torque[j][2] += vxmu2f * tz; } } @@ -288,7 +317,7 @@ void PairLubricate::allocate() void PairLubricate::settings(int narg, char **arg) { - if (narg != 7) error->all("Illegal pair_style command"); + if (narg != 9) error->all("Illegal pair_style command"); mu = atof(arg[0]); flag1 = atoi(arg[1]); @@ -297,6 +326,8 @@ void PairLubricate::settings(int narg, char **arg) flag4 = atoi(arg[4]); cut_inner_global = atof(arg[5]); cut_global = atof(arg[6]); + t_target = atof(arg[7]); + seed = atoi(arg[8]); // reset cutoffs that have been explicitly set @@ -309,6 +340,11 @@ void PairLubricate::settings(int narg, char **arg) cut[i][j] = cut_global; } } + + // initialize Marsaglia RNG with processor-unique seed + + delete random; + random = new RanMars(lmp,seed + comm->me); } /* ---------------------------------------------------------------------- @@ -446,6 +482,8 @@ void PairLubricate::write_restart_settings(FILE *fp) fwrite(&flag4,sizeof(int),1,fp); fwrite(&cut_inner_global,sizeof(double),1,fp); fwrite(&cut_global,sizeof(double),1,fp); + fwrite(&t_target,sizeof(double),1,fp); + fwrite(&seed,sizeof(int),1,fp); fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); } @@ -465,6 +503,8 @@ void PairLubricate::read_restart_settings(FILE *fp) fread(&flag4,sizeof(int),1,fp); fread(&cut_inner_global,sizeof(double),1,fp); fread(&cut_global,sizeof(double),1,fp); + fread(&t_target, sizeof(double),1,fp); + fread(&seed, sizeof(int),1,fp); fread(&offset_flag,sizeof(int),1,fp); fread(&mix_flag,sizeof(int),1,fp); } @@ -475,6 +515,13 @@ void PairLubricate::read_restart_settings(FILE *fp) MPI_Bcast(&flag4,1,MPI_INT,0,world); MPI_Bcast(&cut_inner_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&t_target,1,MPI_DOUBLE,0,world); + MPI_Bcast(&seed,1,MPI_INT,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); + + // additional setup based on restart parameters + + delete random; + random = new RanMars(lmp,seed + comm->me); } diff --git a/src/COLLOID/pair_lubricate.h b/src/COLLOID/pair_lubricate.h index b67f8bddf0..e7db8ab576 100644 --- a/src/COLLOID/pair_lubricate.h +++ b/src/COLLOID/pair_lubricate.h @@ -34,9 +34,12 @@ class PairLubricate : public Pair { protected: double cut_inner_global,cut_global; - double **cut_inner,**cut; - double mu; + double t_target,mu; int flag1,flag2,flag3,flag4; + int seed; + double **cut_inner,**cut; + + class RanMars *random; void allocate(); }; @@ -44,4 +47,3 @@ class PairLubricate : public Pair { } #endif - diff --git a/src/fix_npt.cpp b/src/fix_npt.cpp index 04d3efb30d..c7ede62785 100644 --- a/src/fix_npt.cpp +++ b/src/fix_npt.cpp @@ -63,7 +63,6 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) : double p_period[3]; if (strcmp(arg[6],"xyz") == 0) { if (narg < 10) error->all("Illegal fix npt command"); - press_couple = XYZ; p_start[0] = p_start[1] = p_start[2] = atof(arg[7]); p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[8]);