F = Fc + Ff + Fr +F = Fc + Ff + Fr +Ff = - (m / damp) v +Fr is proportional to sqrt(Kb T m / (dt damp))Fc is the conservative force computed via the usual inter-particle interactions (pair_style, bond_style, etc).
-The Ff and Fr terms are added by this fix. Ff = - gamma v and is a -frictional drag or viscous damping term proportional to the particle's -velocity. Gamma for each atom is computed as m/damp, where m is the -mass of the particle and damp is the damping factor specified by the -user. +
The Ff and Fr terms are added by this fix. +
+Ff is a frictional drag or viscous damping term proportional to the +particle's velocity. The proportionality constant for each atom is +computed as m/damp, where m is the mass of the particle and damp is +the damping factor specified by the user.
Fr is a force due to solvent atoms at a temperature T randomly bumping into the particle. As derived from the fluctuation/dissipation -theorem, its magnitude is proportional to sqrt(T m / dt damp), where T -is the desired temperature, m is the mass of the particle, dt is the -timestep size, and damp is the damping factor. Random numbers are -used to randomize the direction and magnitude of this force as -described in (Dunweg), where a uniform random number is used -(instead of a Gaussian random number) for speed. +theorem, its magnitude as shown above is proportional to sqrt(Kb T m / +dt damp), where Kb is the Boltzmann constant, T is the desired +temperature, m is the mass of the particle, dt is the timestep size, +and damp is the damping factor. Random numbers are used to randomize +the direction and magnitude of this force as described in +(Dunweg), where a uniform random number is used (instead of +a Gaussian random number) for speed.
Note that the thermostat effect of this fix is applied to only the translational degrees of freedom for the particles, which is an diff -Naur lammps-15Apr08/doc/fix_langevin.txt lammps-16Apr08/doc/fix_langevin.txt --- lammps-15Apr08/doc/fix_langevin.txt 2008-04-09 11:29:51.000000000 -0600 +++ lammps-16Apr08/doc/fix_langevin.txt 2008-04-16 16:15:33.000000000 -0600 @@ -37,26 +37,30 @@ nve"_fix_nve.html, this command performs Brownian dynamics (BD), since the total force on each atom will have the form: -F = Fc + Ff + Fr :pre +F = Fc + Ff + Fr +Ff = - (m / damp) v +Fr is proportional to sqrt(Kb T m / (dt damp)) :pre Fc is the conservative force computed via the usual inter-particle interactions ("pair_style"_pair_style.html, "bond_style"_bond_style.html, etc). -The Ff and Fr terms are added by this fix. Ff = - gamma v and is a -frictional drag or viscous damping term proportional to the particle's -velocity. Gamma for each atom is computed as m/damp, where m is the -mass of the particle and damp is the damping factor specified by the -user. +The Ff and Fr terms are added by this fix. + +Ff is a frictional drag or viscous damping term proportional to the +particle's velocity. The proportionality constant for each atom is +computed as m/damp, where m is the mass of the particle and damp is +the damping factor specified by the user. Fr is a force due to solvent atoms at a temperature T randomly bumping into the particle. As derived from the fluctuation/dissipation -theorem, its magnitude is proportional to sqrt(T m / dt damp), where T -is the desired temperature, m is the mass of the particle, dt is the -timestep size, and damp is the damping factor. Random numbers are -used to randomize the direction and magnitude of this force as -described in "(Dunweg)"_#Dunweg, where a uniform random number is used -(instead of a Gaussian random number) for speed. +theorem, its magnitude as shown above is proportional to sqrt(Kb T m / +dt damp), where Kb is the Boltzmann constant, T is the desired +temperature, m is the mass of the particle, dt is the timestep size, +and damp is the damping factor. Random numbers are used to randomize +the direction and magnitude of this force as described in +"(Dunweg)"_#Dunweg, where a uniform random number is used (instead of +a Gaussian random number) for speed. Note that the thermostat effect of this fix is applied to only the translational degrees of freedom for the particles, which is an diff -Naur lammps-15Apr08/doc/pair_lubricate.html lammps-16Apr08/doc/pair_lubricate.html --- lammps-15Apr08/doc/pair_lubricate.html 2008-02-29 18:13:20.000000000 -0700 +++ lammps-16Apr08/doc/pair_lubricate.html 2008-04-16 16:19:08.000000000 -0600 @@ -13,7 +13,7 @@
Syntax:
-pair_style lubricate mu squeeze shear pump twist cutinner cutoff +pair_style lubricate mu squeeze shear pump twist cutinner cutoff T_target seed
Examples:
-pair_style lubricate 1.5 1 1 1 0 2.3 2.4 +pair_style lubricate 1.5 1 1 1 0 2.3 2.4 1.3 5878598 pair_coeff 1 1 1.8 2.0 pair_coeff * *@@ -59,6 +61,12 @@ outer cutoff is the separation distance beyond which the pair-wise forces are zero. +A Langevin thermostatting term is also added to the pairwise force, +similar to that provided by the fix langevin or +pair_style dpd commands. The target temperature for +the thermostat is the specified T_target. The seed is used for +the random numbers generated for the thermostat. +
The following coefficients must be defined for each pair of atoms types via the pair_coeff command as in the examples above, or in the data file or restart files read by the diff -Naur lammps-15Apr08/doc/pair_lubricate.txt lammps-16Apr08/doc/pair_lubricate.txt --- lammps-15Apr08/doc/pair_lubricate.txt 2008-02-29 18:13:20.000000000 -0700 +++ lammps-16Apr08/doc/pair_lubricate.txt 2008-04-16 16:19:08.000000000 -0600 @@ -10,7 +10,7 @@ [Syntax:] -pair_style lubricate mu squeeze shear pump twist cutinner cutoff :pre +pair_style lubricate mu squeeze shear pump twist cutinner cutoff T_target seed :pre mu = viscosity (mass/distance/time units) squeeze = 0/1 for squeeze force off/on @@ -18,11 +18,13 @@ pump = 0/1 for pump force off/on twist = 0/1 for twist force off/on cutinner = (distance units) -cutoff = outer cutoff for interactions (distance units) :ul +cutoff = outer cutoff for interactions (distance units) +T_target = desired temperature (temperature units) +seed = random number seed (positive integer) :ul [Examples:] -pair_style lubricate 1.5 1 1 1 0 2.3 2.4 +pair_style lubricate 1.5 1 1 1 0 2.3 2.4 1.3 5878598 pair_coeff 1 1 1.8 2.0 pair_coeff * * :pre @@ -56,6 +58,12 @@ outer {cutoff} is the separation distance beyond which the pair-wise forces are zero. +A Langevin thermostatting term is also added to the pairwise force, +similar to that provided by the "fix langevin"_fix_langevin.html or +"pair_style dpd"_pair_dpd.html commands. The target temperature for +the thermostat is the specified {T_target}. The {seed} is used for +the random numbers generated for the thermostat. + The following coefficients must be defined for each pair of atoms types via the "pair_coeff"_pair_coeff.html command as in the examples above, or in the data file or restart files read by the diff -Naur lammps-15Apr08/src/COLLOID/pair_lubricate.cpp lammps-16Apr08/src/COLLOID/pair_lubricate.cpp --- lammps-15Apr08/src/COLLOID/pair_lubricate.cpp 2008-04-14 16:14:13.000000000 -0600 +++ lammps-16Apr08/src/COLLOID/pair_lubricate.cpp 2008-04-16 16:00:37.000000000 -0600 @@ -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 @@ PairLubricate::PairLubricate(LAMMPS *lmp) : Pair(lmp) { single_enable = 0; + + random = NULL; } /* ---------------------------------------------------------------------- */ @@ -54,6 +56,7 @@ memory->destroy_2d_double_array(cut); memory->destroy_2d_double_array(cut_inner); } + delete random; } /* ---------------------------------------------------------------------- */ @@ -61,8 +64,8 @@ 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 @@ 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 @@ 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 @@ } 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 @@ 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 @@ (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 @@ 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::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 @@ 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 @@ 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 @@ 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 @@ 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 @@ 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 -Naur lammps-15Apr08/src/COLLOID/pair_lubricate.h lammps-16Apr08/src/COLLOID/pair_lubricate.h --- lammps-15Apr08/src/COLLOID/pair_lubricate.h 2007-10-11 17:56:12.000000000 -0600 +++ lammps-16Apr08/src/COLLOID/pair_lubricate.h 2008-04-16 16:00:37.000000000 -0600 @@ -34,9 +34,12 @@ 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 @@ } #endif - diff -Naur lammps-15Apr08/src/fix_npt.cpp lammps-16Apr08/src/fix_npt.cpp --- lammps-15Apr08/src/fix_npt.cpp 2008-04-10 14:53:37.000000000 -0600 +++ lammps-16Apr08/src/fix_npt.cpp 2008-04-16 16:00:37.000000000 -0600 @@ -63,7 +63,6 @@ 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]); diff -Naur lammps-15Apr08/src/min.h lammps-16Apr08/src/min.h --- lammps-15Apr08/src/min.h 2008-04-14 14:50:48.000000000 -0600 +++ lammps-16Apr08/src/min.h 2008-04-15 15:43:15.000000000 -0600 @@ -25,7 +25,7 @@ double alpha_final; int niter,neval; char *stopstr; - double dmin,dmax; + double dmax; int linestyle,lineiter; Min(class LAMMPS *); diff -Naur lammps-15Apr08/src/min_cg.cpp lammps-16Apr08/src/min_cg.cpp --- lammps-15Apr08/src/min_cg.cpp 2008-04-14 16:13:25.000000000 -0600 +++ lammps-16Apr08/src/min_cg.cpp 2008-04-15 15:43:15.000000000 -0600 @@ -294,8 +294,7 @@ // line minimization along direction h from current atom->x eprevious = ecurrent; - fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmin,dmax, - alpha_final,neval); + fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmax,alpha_final,neval); if (fail) return FAIL; // function evaluation criterion @@ -479,7 +478,7 @@ x = ptr to atom->x[0] as vector dir = search direction as vector eng = current energy at initial x - min/max dist = min/max distance to move any atom coord + maxdist = max distance to move any atom coord output: return 0 if successful move, non-zero alpha alpha return 1 if failed, alpha = 0.0 alpha = distance moved along dir to set x to minimun eng config @@ -502,45 +501,48 @@ ------------------------------------------------------------------------- */ int MinCG::linemin_backtrack(int n, double *x, double *dir, double eng, - double mindist, double maxdist, - double &alpha, int &nfunc) + double maxdist, double &alpha, int &nfunc) { - int i,ilist,m; - double fdotdirme,fdotdirall,fme,fmax,eoriginal,alphamax; + int i,m; // stopping criterion, must be scaled by normflag double *f = atom->f[0]; - fdotdirme = 0.0; + double fdotdirme = 0.0; for (i = 0; i < n; i++) fdotdirme += f[i]*dir[i]; + double fdotdirall; MPI_Allreduce(&fdotdirme,&fdotdirall,1,MPI_DOUBLE,MPI_SUM,world); if (output->thermo->normflag) fdotdirall /= atom->natoms; // alphamax = step that moves some atom coord by maxdist - fme = 0.0; + double fme = 0.0; for (i = 0; i < n; i++) fme = MAX(fme,fabs(dir[i])); + double fmax; MPI_Allreduce(&fme,&fmax,1,MPI_DOUBLE,MPI_MAX,world); if (fmax == 0.0) return 1; - alphamax = maxdist/fmax; + double alphamax = maxdist/fmax; - // reduce alpha by factor of ALPHA_REDUCE until energy decrease is sufficient + // reduce alpha by ALPHA_REDUCE until energy decrease is sufficient - eoriginal = eng; + double eoriginal = eng; alpha = alphamax; + double alpha_previous = 0.0; + double delta; while (1) { - for (i = 0; i < n; i++) x[i] += alpha*dir[i]; + delta = alpha - alpha_previous; + for (i = 0; i < n; i++) x[i] += delta*dir[i]; eng_force(&n,&x,&dir,&eng); nfunc++; if (eng <= eoriginal - BACKTRACK_SLOPE*alpha*fdotdirall) return 0; - for (i = 0; i < n; i++) x[i] -= alpha*dir[i]; - + alpha_previous = alpha; alpha *= ALPHA_REDUCE; if (alpha == 0.0) { + for (i = 0; i < n; i++) x[i] -= alpha_previous*dir[i]; eng_force(&n,&x,&dir,&eng); return 1; } diff -Naur lammps-15Apr08/src/min_cg.h lammps-16Apr08/src/min_cg.h --- lammps-15Apr08/src/min_cg.h 2008-04-14 14:50:48.000000000 -0600 +++ lammps-16Apr08/src/min_cg.h 2008-04-15 15:43:15.000000000 -0600 @@ -32,18 +32,20 @@ int triclinic; // 0 if domain is orthog, 1 if triclinic class FixMinimize *fix_minimize; // fix that stores gradient vecs - class Compute *pe_compute; // compute for potential energy - double ecurrent; // current potential energy + class Compute *pe_compute; // compute for potential energy + double ecurrent; // current potential energy double mindist,maxdist; // min/max dist for coord delta in line search int ndof; // # of degrees-of-freedom on this proc double *g,*h; // local portion of gradient, searchdir vectors + // ptr to linemin functions + typedef int (MinCG::*FnPtr)(int, double *, double *, double, - double, double, double &, int &); - FnPtr linemin; // ptr to linemin functions + double, double &, int &); + FnPtr linemin; int linemin_backtrack(int, double *, double *, double, - double, double, double &, int &); + double, double &, int &); void setup(); void setup_vectors(); diff -Naur lammps-15Apr08/src/min_sd.cpp lammps-16Apr08/src/min_sd.cpp --- lammps-15Apr08/src/min_sd.cpp 2008-04-14 16:13:25.000000000 -0600 +++ lammps-16Apr08/src/min_sd.cpp 2008-04-15 15:43:15.000000000 -0600 @@ -51,8 +51,7 @@ // line minimization along direction h from current atom->x eprevious = ecurrent; - fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmin,dmax, - alpha_final,neval); + fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmax,alpha_final,neval); if (fail) return FAIL; // function evaluation criterion