diff --git a/doc/html/boundary.html b/doc/html/boundary.html index 75b29e927c..4166285302 100644 --- a/doc/html/boundary.html +++ b/doc/html/boundary.html @@ -173,7 +173,15 @@ you have set the

For style s, the position of the face is set so as to encompass the atoms in that dimension (shrink-wrapping), no matter how far they -move.

+move. Note that when the difference between the current box dimensions +and the shrink-wrap box dimensions is large, this can lead to lost +atoms at the beginning of a run when running in parallel. This is due +to the large change in the (global) box dimensions also causing +significant changes in the individual sub-domain sizes. If these +changes are farther than the communication cutoff, atoms will be lost. +This is best addressed by setting initial box dimensions to match the +shrink-wrapped dimensions more closely, by using m style boundaries +(see below).

For style m, shrink-wrapping occurs, but is bounded by the value specified in the data or restart file or set by the create_box command. For example, if the upper z diff --git a/doc/html/compute_rdf.html b/doc/html/compute_rdf.html index ef1d1d8cfc..7565e4888f 100644 --- a/doc/html/compute_rdf.html +++ b/doc/html/compute_rdf.html @@ -141,13 +141,13 @@

Examples

-
compute 1 all rdf 100
-compute 1 all rdf 100 1 1
-compute 1 all rdf 100 * 3
-compute 1 fluid rdf 500 1 1 1 2 2 1 2 2
-compute 1 fluid rdf 500 1*3 2 5 *10
-
-
+
+compute 1 all rdf 100
+compute 1 all rdf 100 1 1
+compute 1 all rdf 100 * 3
+compute 1 fluid rdf 500 1 1 1 2 2 1 2 2
+compute 1 fluid rdf 500 1*3 2 5 *10
+

Description

@@ -184,7 +184,7 @@ listed, then a separate histogram is generated for each

The itypeN and jtypeN settings can be specified in one of two ways. An explicit numeric value can be used, as in the 4th example above. Or a wild-card asterisk can be used to specify a range of atom -types. This takes the form “*” or “n” or “n” or “m*n”. If N = the +types. This takes the form “*” or “*n” or “n*” or “m*n”. If N = the number of atom types, then an asterisk with no numeric values means all types from 1 to N. A leading asterisk means all types from 1 to n (inclusive). A trailing asterisk means all types from n to N @@ -226,10 +226,10 @@ atoms of type jtypeN.

The simplest way to output the results of the compute rdf calculation to a file is to use the fix ave/time command, for example:

-
compute myRDF all rdf 50
-fix 1 all ave/time 100 1 100 c_myRDF[*] file tmp.rdf mode vector
-
-
+
+compute myRDF all rdf 50
+fix 1 all ave/time 100 1 100 c_myRDF[*] file tmp.rdf mode vector
+

Output info:

This compute calculates a global array with the number of rows = Nbins, and the number of columns = 1 + 2*Npairs, where Npairs is the diff --git a/doc/html/fix_adapt.html b/doc/html/fix_adapt.html index 2f5eece4a1..45eee9b4af 100644 --- a/doc/html/fix_adapt.html +++ b/doc/html/fix_adapt.html @@ -165,12 +165,12 @@

Examples

-
fix 1 all adapt 1 pair soft a 1 1 v_prefactor
-fix 1 all adapt 1 pair soft a 2* 3 v_prefactor
-fix 1 all adapt 1 pair lj/cut epsilon * * v_scale1 coul/cut scale 3 3 v_scale2 scale yes reset yes
-fix 1 all adapt 10 atom diameter v_size
-
-
+
+fix 1 all adapt 1 pair soft a 1 1 v_prefactor
+fix 1 all adapt 1 pair soft a 2* 3 v_prefactor
+fix 1 all adapt 1 pair lj/cut epsilon * * v_scale1 coul/cut scale 3 3 v_scale2 scale yes reset yes
+fix 1 all adapt 10 atom diameter v_size
+

Description

@@ -225,8 +225,8 @@ meaning of these parameters:

--++ @@ -257,23 +257,27 @@ meaning of these parameters:

- + + + + + - + - + - + - + @@ -315,7 +319,7 @@ each, as in the 1st example above. I <= J is required. LAMMPS sets the coefficients for the symmetric J,I interaction to the same values.

A wild-card asterisk can be used in place of or in conjunction with the I,J arguments to set the coefficients for multiple pairs of atom -types. This takes the form “*” or “n” or “n” or “m*n”. If N = the +types. This takes the form “*” or “*n” or “n*” or “m*n”. If N = the number of atom types, then an asterisk with no numeric values means all types from 1 to N. A leading asterisk means all types from 1 to n (inclusive). A trailing asterisk means all types from n to N @@ -341,10 +345,10 @@ details.

For example, these commands would change the prefactor coefficient of the pair_style soft potential from 10.0 to 30.0 in a linear fashion over the course of a simulation:

-
variable prefactor equal ramp(10,30)
-fix 1 all adapt 1 pair soft a * * v_prefactor
-
-
+
+variable prefactor equal ramp(10,30)
+fix 1 all adapt 1 pair soft a * * v_prefactor
+

The kspace keyword used the specified variable as a scale factor on the energy, forces, virial calculated by whatever K-Space solver is @@ -380,14 +384,12 @@ constant).

For example, these commands would shrink the diameter of all granular particles in the “center” group from 1.0 to 0.1 in a linear fashion over the course of a 1000-step simulation:

-
variable size equal ramp(1.0,0.1)
-fix 1 center adapt 10 atom diameter v_size
-
-
- +
+variable size equal ramp(1.0,0.1)
+fix 1 center adapt 10 atom diameter v_size
+

-
-

Restart, fix_modify, output, run start/stop, minimize info

+

Restart, fix_modify, output, run start/stop, minimize info:

No information about this fix is written to binary restart files. None of the fix_modify options are relevant to this fix. No global or per-atom quantities are stored by this fix for access by various output commands. No parameter of this fix can diff --git a/doc/html/pair_dipole.html b/doc/html/pair_dipole.html index a8b5437180..6b903c7a92 100644 --- a/doc/html/pair_dipole.html +++ b/doc/html/pair_dipole.html @@ -149,12 +149,12 @@

pair_style lj/long/dipole/long command

Syntax

-
pair_style lj/cut/dipole/cut cutoff (cutoff2)
-pair_style lj/sf/dipole/sf cutoff (cutoff2)
-pair_style lj/cut/dipole/long cutoff (cutoff2)
-pair_style lj/long/dipole/long flag_lj flag_coul cutoff (cutoff2)
-
-
+
+pair_style lj/cut/dipole/cut cutoff (cutoff2)
+pair_style lj/sf/dipole/sf cutoff (cutoff2)
+pair_style lj/cut/dipole/long cutoff (cutoff2)
+pair_style lj/long/dipole/long flag_lj flag_coul cutoff (cutoff2)
+
  • cutoff = global cutoff LJ (and Coulombic if only 1 arg) (distance units)
  • cutoff2 = global cutoff for Coulombic and dipole (optional) (distance units)
  • @@ -175,26 +175,27 @@

Examples

-
pair_style lj/cut/dipole/cut 10.0
-pair_coeff * * 1.0 1.0
-pair_coeff 2 3 1.0 1.0 2.5 4.0
-
-
-
pair_style lj/sf/dipole/sf 9.0
-pair_coeff * * 1.0 1.0
-pair_coeff 2 3 1.0 1.0 2.5 4.0
-
-
-
pair_style lj/cut/dipole/long 10.0
-pair_coeff * * 1.0 1.0
-pair_coeff 2 3 1.0 1.0 2.5 4.0
-
-
-
pair_style lj/long/dipole/long long long 3.5 10.0
-pair_coeff * * 1.0 1.0
-pair_coeff 2 3 1.0 1.0 2.5 4.0
-
-
+
+pair_style lj/cut/dipole/cut 10.0
+pair_coeff * * 1.0 1.0
+pair_coeff 2 3 1.0 1.0 2.5 4.0
+
+
+pair_style lj/sf/dipole/sf 9.0
+pair_coeff * * 1.0 1.0
+pair_coeff 2 3 1.0 1.0 2.5 4.0 scale 0.5
+pair_coeff 2 3 1.0 1.0 2.5 4.0
+
+
+pair_style lj/cut/dipole/long 10.0
+pair_coeff * * 1.0 1.0
+pair_coeff 2 3 1.0 1.0 2.5 4.0
+
+
+pair_style lj/long/dipole/long long long 3.5 10.0
+pair_coeff * * 1.0 1.0
+pair_coeff 2 3 1.0 1.0 2.5 4.0
+

Description

@@ -242,7 +243,10 @@ C.3 of If one cutoff is specified in the pair_style command, it is used for both the LJ and Coulombic (q,p) terms. If two cutoffs are specified, they are used as cutoffs for the LJ and Coulombic (q,p) terms -respectively.

+respectively. This pair style also supports an optional scale keyword +as part of a pair_coeff statement, where the interactions can be +scaled according to this factor. This scale factor is also made available +for use with fix adapt.

Style lj/cut/dipole/long computes long-range point-dipole interactions as discussed in (Toukmaji). Dipole-dipole, dipole-charge, and charge-charge interactions are all supported, along diff --git a/doc/html/variable.html b/doc/html/variable.html index f5eb258506..2b9cdb1f79 100644 --- a/doc/html/variable.html +++ b/doc/html/variable.html @@ -965,11 +965,24 @@ if $(is_active(pair,respa)) then "run_style respa 4 3 2 2 improper 1 inner

The is_defined() function allows to query categories like compute, dump, fix, group, region, and variable whether an entry with the provided name or id is defined.

-

The is_available() function allows to query whether a specific -optional feature is available, i.e. compiled in. This currently -works for the following categories: command, compute, fix, -and pair_style. For all categories except command also appending -active suffixes is tried before reporting failure.

+

The is_available(category,name) function allows to query whether +a specific optional feature is available, i.e. compiled in. +This currently works for the following categories: command, +compute, fix, pair_style and feature. For all categories +except command and feature also appending active suffixes is +tried before reporting failure.

+

The feature category is used to check the availability of compiled in +features such as GZIP support, PNG support, JPEG support, FFMPEG support, +and C++ exceptions for error handling. Corresponding values for name are +gzip, png, jpeg, ffmpeg and exceptions.

+

This enables writing input scripts which only dump using a given format if +the compiled binary supports it.

+
+if "$(is_available(feature,png))" then "print 'PNG supported'" else "print 'PNG not supported'"
+
+
+if "$(is_available(feature,ffmpeg)" then "dump 3 all movie 25 movie.mp4 type type zoom 1.6 adiam 1.0"
+

diff --git a/doc/src/boundary.txt b/doc/src/boundary.txt index e832a74d24..91096589a7 100644 --- a/doc/src/boundary.txt +++ b/doc/src/boundary.txt @@ -54,7 +54,15 @@ allow for lost atoms. For style {s}, the position of the face is set so as to encompass the atoms in that dimension (shrink-wrapping), no matter how far they -move. +move. Note that when the difference between the current box dimensions +and the shrink-wrap box dimensions is large, this can lead to lost +atoms at the beginning of a run when running in parallel. This is due +to the large change in the (global) box dimensions also causing +significant changes in the individual sub-domain sizes. If these +changes are farther than the communication cutoff, atoms will be lost. +This is best addressed by setting initial box dimensions to match the +shrink-wrapped dimensions more closely, by using {m} style boundaries +(see below). For style {m}, shrink-wrapping occurs, but is bounded by the value specified in the data or restart file or set by the diff --git a/doc/src/fix_adapt.txt b/doc/src/fix_adapt.txt index e13f877748..712626fd25 100644 --- a/doc/src/fix_adapt.txt +++ b/doc/src/fix_adapt.txt @@ -110,6 +110,7 @@ meaning of these parameters: "coul/long"_pair_coul.html: scale: type pairs: "lj/cut"_pair_lj.html: epsilon,sigma: type pairs: "lj/expand"_pair_lj_expand.html: epsilon,sigma,delta: type pairs: +"lj/sf/dipole/sf"_pair_dipole.html: epsilon,sigma,scale: type pairs: "lubricate"_pair_lubricate.html: mu: global: "gauss"_pair_gauss.html: a: type pairs: "morse"_pair_morse.html: d0,r0,alpha: type pairs: diff --git a/doc/src/pair_dipole.txt b/doc/src/pair_dipole.txt index be6608d57b..b3eefa63f2 100755 --- a/doc/src/pair_dipole.txt +++ b/doc/src/pair_dipole.txt @@ -41,6 +41,7 @@ pair_coeff 2 3 1.0 1.0 2.5 4.0 :pre pair_style lj/sf/dipole/sf 9.0 pair_coeff * * 1.0 1.0 +pair_coeff 2 3 1.0 1.0 2.5 4.0 scale 0.5 pair_coeff 2 3 1.0 1.0 2.5 4.0 :pre pair_style lj/cut/dipole/long 10.0 @@ -103,7 +104,10 @@ C.3 of "(Allen)"_#Allen. If one cutoff is specified in the pair_style command, it is used for both the LJ and Coulombic (q,p) terms. If two cutoffs are specified, they are used as cutoffs for the LJ and Coulombic (q,p) terms -respectively. +respectively. This pair style also supports an optional {scale} keyword +as part of a pair_coeff statement, where the interactions can be +scaled according to this factor. This scale factor is also made available +for use with fix adapt. Style {lj/cut/dipole/long} computes long-range point-dipole interactions as discussed in "(Toukmaji)"_#Toukmaji. Dipole-dipole, diff --git a/doc/src/variable.txt b/doc/src/variable.txt index 9f40beacb4..7940f2b5e9 100644 --- a/doc/src/variable.txt +++ b/doc/src/variable.txt @@ -894,11 +894,25 @@ The {is_defined()} function allows to query categories like {compute}, {dump}, {fix}, {group}, {region}, and {variable} whether an entry with the provided name or id is defined. -The {is_available()} function allows to query whether a specific -optional feature is available, i.e. compiled in. This currently -works for the following categories: {command}, {compute}, {fix}, -and {pair_style}. For all categories except {command} also appending -active suffixes is tried before reporting failure. +The {is_available(category,name)} function allows to query whether +a specific optional feature is available, i.e. compiled in. +This currently works for the following categories: {command}, +{compute}, {fix}, {pair_style} and {feature}. For all categories +except {command} and {feature} also appending active suffixes is +tried before reporting failure. + +The {feature} category is used to check the availability of compiled in +features such as GZIP support, PNG support, JPEG support, FFMPEG support, +and C++ exceptions for error handling. Corresponding values for name are +{gzip}, {png}, {jpeg}, {ffmpeg} and {exceptions}. + +This enables writing input scripts which only dump using a given format if +the compiled binary supports it. + +if "$(is_available(feature,png))" then "print 'PNG supported'" else "print 'PNG not supported'" :pre + +if "$(is_available(feature,ffmpeg)" then "dump 3 all movie 25 movie.mp4 type type zoom 1.6 adiam 1.0" :pre + :line diff --git a/src/MEAM/pair_meam.cpp b/src/MEAM/pair_meam.cpp index 2d14ffa0cf..bc386d0320 100644 --- a/src/MEAM/pair_meam.cpp +++ b/src/MEAM/pair_meam.cpp @@ -57,6 +57,7 @@ PairMEAM::PairMEAM(LAMMPS *lmp) : Pair(lmp) arho1 = arho2 = arho3 = arho3b = t_ave = tsq_ave = NULL; maxneigh = 0; + allocated = 0; scrfcn = dscrfcn = fcpair = NULL; nelements = 0; diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index a963e07a40..2f2f84570e 100644 --- a/src/RIGID/fix_rigid.cpp +++ b/src/RIGID/fix_rigid.cpp @@ -2006,7 +2006,7 @@ void FixRigid::setup_bodies_static() MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); - // error check that re-computed momemts of inertia match diagonalized ones + // error check that re-computed moments of inertia match diagonalized ones // do not do test for bodies with params read from infile double norm; diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index 36dfc1bdea..bc2cd04bcc 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -2192,7 +2192,7 @@ void FixRigidSmall::setup_bodies_static() commflag = ITENSOR; comm->reverse_comm_fix(this,6); - // error check that re-computed momemts of inertia match diagonalized ones + // error check that re-computed moments of inertia match diagonalized ones // do not do test for bodies with params read from infile double norm; diff --git a/src/USER-COLVARS/fix_colvars.cpp b/src/USER-COLVARS/fix_colvars.cpp index 4f85c882d4..9977625444 100644 --- a/src/USER-COLVARS/fix_colvars.cpp +++ b/src/USER-COLVARS/fix_colvars.cpp @@ -383,7 +383,7 @@ void FixColvars::init() error->all(FLERR,"Cannot use fix colvars without atom IDs"); if (atom->map_style == 0) - error->all(FLERR,"Fix colvars requires an atom map"); + error->all(FLERR,"Fix colvars requires an atom map, see atom_modify"); if ((me == 0) && (update->whichflag == 2)) error->warning(FLERR,"Using fix colvars with minimization"); diff --git a/src/USER-COLVARS/fix_colvars.h b/src/USER-COLVARS/fix_colvars.h index b87c4c3ed5..a1dae0d757 100644 --- a/src/USER-COLVARS/fix_colvars.h +++ b/src/USER-COLVARS/fix_colvars.h @@ -128,7 +128,7 @@ E: Cannot use fix colvars without atom IDs Atom IDs are not defined, but fix colvars needs them to identify an atom. -E: Fix colvars requires an atom map +E: Fix colvars requires an atom map, see atom_modify Use the atom_modify command to create an atom map. diff --git a/src/USER-MISC/pair_lj_sf_dipole_sf.cpp b/src/USER-MISC/pair_lj_sf_dipole_sf.cpp index 0a4ddd7908..072a6c9943 100755 --- a/src/USER-MISC/pair_lj_sf_dipole_sf.cpp +++ b/src/USER-MISC/pair_lj_sf_dipole_sf.cpp @@ -12,7 +12,8 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Mario Orsi (U Southampton), orsimario@gmail.com + Contributing authors: Mario Orsi (QMUL), m.orsi@qmul.ac.uk + Samuel Genheden (University of Southampton) ------------------------------------------------------------------------- */ #include @@ -34,7 +35,6 @@ using namespace LAMMPS_NS; PairLJSFDipoleSF::PairLJSFDipoleSF(LAMMPS *lmp) : Pair(lmp) { - single_enable = 0; } /* ---------------------------------------------------------------------- */ @@ -55,6 +55,7 @@ PairLJSFDipoleSF::~PairLJSFDipoleSF() memory->destroy(lj2); memory->destroy(lj3); memory->destroy(lj4); + memory->destroy(scale); } } @@ -86,6 +87,7 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag) double **torque = atom->torque; int *type = atom->type; int nlocal = atom->nlocal; +// The global scaling parameters aren't used anymore double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -234,7 +236,7 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag) // total force - fq = factor_coul*qqrd2e; + fq = factor_coul*qqrd2e*scale[itype][jtype]; fx = fq*forcecoulx + delx*forcelj; fy = fq*forcecouly + dely*forcelj; fz = fq*forcecoulz + delz*forcelj; @@ -268,7 +270,7 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag) ecoul += -q[j] * r3inv * pqfac * pidotr; if (mu[j][3] > 0.0 && qtmp != 0.0) ecoul += qtmp * r3inv * qpfac * pjdotr; - ecoul *= factor_coul*qqrd2e; + ecoul *= factor_coul*qqrd2e*scale[itype][jtype]; } else ecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { @@ -315,6 +317,7 @@ void PairLJSFDipoleSF::allocate() memory->create(lj2,n+1,n+1,"pair:lj2"); memory->create(lj3,n+1,n+1,"pair:lj3"); memory->create(lj4,n+1,n+1,"pair:lj4"); + memory->create(scale,n+1,n+1,"pair:scale"); } /* ---------------------------------------------------------------------- @@ -352,7 +355,7 @@ void PairLJSFDipoleSF::settings(int narg, char **arg) void PairLJSFDipoleSF::coeff(int narg, char **arg) { - if (narg < 4 || narg > 6) + if (narg < 4 || narg > 7) error->all(FLERR,"Incorrect args for pair coefficients"); if (!allocated) allocate(); @@ -363,10 +366,13 @@ void PairLJSFDipoleSF::coeff(int narg, char **arg) double epsilon_one = force->numeric(FLERR,arg[2]); double sigma_one = force->numeric(FLERR,arg[3]); + double scale_one = 1.0; + if (narg >= 5) scale_one = force->numeric(FLERR,arg[4]); + double cut_lj_one = cut_lj_global; double cut_coul_one = cut_coul_global; - if (narg >= 5) cut_coul_one = cut_lj_one = force->numeric(FLERR,arg[4]); - if (narg == 6) cut_coul_one = force->numeric(FLERR,arg[5]); + if (narg >= 6) cut_coul_one = cut_lj_one = force->numeric(FLERR,arg[5]); + if (narg == 7) cut_coul_one = force->numeric(FLERR,arg[6]); int count = 0; for (int i = ilo; i <= ihi; i++) { @@ -376,6 +382,7 @@ void PairLJSFDipoleSF::coeff(int narg, char **arg) cut_lj[i][j] = cut_lj_one; cut_coul[i][j] = cut_coul_one; setflag[i][j] = 1; + scale[i][j] = scale_one; count++; } } @@ -424,6 +431,7 @@ double PairLJSFDipoleSF::init_one(int i, int j) lj2[j][i] = lj2[i][j]; lj3[j][i] = lj3[i][j]; lj4[j][i] = lj4[i][j]; + scale[j][i] = scale[i][j]; return cut; } @@ -445,6 +453,7 @@ void PairLJSFDipoleSF::write_restart(FILE *fp) fwrite(&sigma[i][j],sizeof(double),1,fp); fwrite(&cut_lj[i][j],sizeof(double),1,fp); fwrite(&cut_coul[i][j],sizeof(double),1,fp); + fwrite(&scale[i][j],sizeof(double),1,fp); } } } @@ -471,11 +480,13 @@ void PairLJSFDipoleSF::read_restart(FILE *fp) fread(&sigma[i][j],sizeof(double),1,fp); fread(&cut_lj[i][j],sizeof(double),1,fp); fread(&cut_coul[i][j],sizeof(double),1,fp); + fread(&scale[i][j],sizeof(double),1,fp); } MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_coul[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&scale[i][j],1,MPI_DOUBLE,0,world); } } } @@ -506,3 +517,105 @@ void PairLJSFDipoleSF::read_restart_settings(FILE *fp) MPI_Bcast(&cut_coul_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } + +// PairLJSFDipoleSF: calculation of force is missing (to be implemented) +double PairLJSFDipoleSF::single(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double r2inv,r6inv; + double pdotp,pidotr,pjdotr,pre1,delx,dely,delz; + double rinv, r3inv,r5inv, rcutlj2inv, rcutcoul2inv,rcutlj6inv; + double qtmp,xtmp,ytmp,ztmp,bfac,pqfac,qpfac, ecoul, evdwl; + + double **x = atom->x; + double *q = atom->q; + double **mu = atom->mu; + + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + r2inv = 1.0/rsq; + rinv = sqrt(r2inv); + fforce = 0.0; + + if (rsq < cut_coulsq[itype][jtype]) { + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + // if (qtmp != 0.0 && q[j] != 0.0) { + // pre1 = qtmp*q[j]*rinv*(r2inv-1.0/cut_coulsq[itype][jtype]); + // forcecoulx += pre1*delx; + // forcecouly += pre1*dely; + // forcecoulz += pre1*delz; + // } + if (mu[i][3] > 0.0 && mu[j][3] > 0.0) { + r3inv = r2inv*rinv; + r5inv = r3inv*r2inv; + rcutcoul2inv=1.0/cut_coulsq[itype][jtype]; + pdotp = mu[i][0]*mu[j][0] + mu[i][1]*mu[j][1] + mu[i][2]*mu[j][2]; + pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz; + pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz; + bfac = 1.0 - 4.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv) + + 3.0*rsq*rsq*rcutcoul2inv*rcutcoul2inv; + } + if (mu[i][3] > 0.0 && q[j] != 0.0) { + r3inv = r2inv*rinv; + r5inv = r3inv*r2inv; + pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz; + rcutcoul2inv=1.0/cut_coulsq[itype][jtype]; + pqfac = 1.0 - 3.0*rsq*rcutcoul2inv + + 2.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv); + } + if (mu[j][3] > 0.0 && qtmp != 0.0) { + r3inv = r2inv*rinv; + r5inv = r3inv*r2inv; + pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz; + rcutcoul2inv=1.0/cut_coulsq[itype][jtype]; + qpfac = 1.0 - 3.0*rsq*rcutcoul2inv + + 2.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv); + } + } + if (rsq < cut_ljsq[itype][jtype]) { + r6inv = r2inv*r2inv*r2inv; + rcutlj2inv = 1.0 / cut_ljsq[itype][jtype]; + rcutlj6inv = rcutlj2inv * rcutlj2inv * rcutlj2inv; + } + + double eng = 0.0; + if (rsq < cut_coulsq[itype][jtype]) { + ecoul = (1.0-sqrt(rsq)/sqrt(cut_coulsq[itype][jtype])); + ecoul *= ecoul; + ecoul *= qtmp * q[j] * rinv; + if (mu[i][3] > 0.0 && mu[j][3] > 0.0) + ecoul += bfac * (r3inv*pdotp - 3.0*r5inv*pidotr*pjdotr); + if (mu[i][3] > 0.0 && q[j] != 0.0) + ecoul += -q[j] * r3inv * pqfac * pidotr; + if (mu[j][3] > 0.0 && qtmp != 0.0) + ecoul += qtmp * r3inv * qpfac * pjdotr; + ecoul *= factor_coul*force->qqrd2e*scale[itype][jtype]; + eng += ecoul; + } + if (rsq < cut_ljsq[itype][jtype]) { + evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype])+ + rcutlj6inv*(6*lj3[itype][jtype]*rcutlj6inv-3*lj4[itype][jtype])* + rsq*rcutlj2inv+ + rcutlj6inv*(-7*lj3[itype][jtype]*rcutlj6inv+4*lj4[itype][jtype]); + eng += evdwl*factor_lj; + } + + return eng; +} + +/* ---------------------------------------------------------------------- */ + +void *PairLJSFDipoleSF::extract(const char *str, int &dim) +{ + dim = 2; + if (strcmp(str,"epsilon") == 0) return (void *) epsilon; + if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str,"scale") == 0) return (void *) scale; + return NULL; +} diff --git a/src/USER-MISC/pair_lj_sf_dipole_sf.h b/src/USER-MISC/pair_lj_sf_dipole_sf.h index 1a6c57b2b1..76cc058e69 100755 --- a/src/USER-MISC/pair_lj_sf_dipole_sf.h +++ b/src/USER-MISC/pair_lj_sf_dipole_sf.h @@ -37,6 +37,8 @@ class PairLJSFDipoleSF : public Pair { void read_restart(FILE *); void write_restart_settings(FILE *); void read_restart_settings(FILE *); + virtual double single(int, int, int, int, double, double, double, double &); + void *extract(const char *, int &); protected: double cut_lj_global,cut_coul_global; @@ -44,6 +46,7 @@ class PairLJSFDipoleSF : public Pair { double **cut_coul,**cut_coulsq; double **epsilon,**sigma; double **lj1,**lj2,**lj3,**lj4; + double **scale; void allocate(); }; diff --git a/src/USER-PHONON/fix_phonon.cpp b/src/USER-PHONON/fix_phonon.cpp index d675736179..1470bc3eef 100644 --- a/src/USER-PHONON/fix_phonon.cpp +++ b/src/USER-PHONON/fix_phonon.cpp @@ -190,7 +190,7 @@ FixPhonon::FixPhonon(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) sprintf(str,"Can not open output file %s",logfile); error->one(FLERR,str); } - for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); + fprintf(flog,"############################################################\n"); fprintf(flog,"# group name of the atoms under study : %s\n", group->names[igroup]); fprintf(flog,"# total number of atoms in the group : %d\n", ngroup); fprintf(flog,"# dimension of the system : %d D\n", sysdim); @@ -200,7 +200,7 @@ FixPhonon::FixPhonon(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) fprintf(flog,"# frequency of the measurement : %d\n", nevery); fprintf(flog,"# output result after this many measurement: %d\n", nfreq); fprintf(flog,"# number of processors used by this run : %d\n", nprocs); - for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); + fprintf(flog,"############################################################\n"); fprintf(flog,"# mapping information between lattice indices and atom id\n"); fprintf(flog,"# nx ny nz nucell\n"); fprintf(flog,"%d %d %d %d\n", nx, ny, nz, nucell); @@ -214,7 +214,7 @@ FixPhonon::FixPhonon(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) ix = (idx/(nucell*nz*ny))%nx; fprintf(flog,"%d %d %d %d " TAGINT_FORMAT "\n", ix, iy, iz, iu, itag); } - for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); + fprintf(flog,"############################################################\n"); fflush(flog); } surf2tag.clear(); @@ -734,16 +734,16 @@ void FixPhonon::postprocess( ) fclose(fp_bin); // write log file, here however, it is the dynamical matrix that is written - for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); - fprintf(flog, "# Current time step : " BIGINT_FORMAT "\n", update->ntimestep); - fprintf(flog, "# Total number of measurements : %d\n", neval); - fprintf(flog, "# Average temperature of the measurement : %lg\n", TempAve); - fprintf(flog, "# Boltzmann constant under current units : %lg\n", boltz); - fprintf(flog, "# basis vector A1 = [%lg %lg %lg]\n", basevec[0], basevec[1], basevec[2]); - fprintf(flog, "# basis vector A2 = [%lg %lg %lg]\n", basevec[3], basevec[4], basevec[5]); - fprintf(flog, "# basis vector A3 = [%lg %lg %lg]\n", basevec[6], basevec[7], basevec[8]); - for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); - fprintf(flog, "# qx\t qy \t qz \t\t Phi(q)\n"); + fprintf(flog,"############################################################\n"); + fprintf(flog,"# Current time step : " BIGINT_FORMAT "\n", update->ntimestep); + fprintf(flog,"# Total number of measurements : %d\n", neval); + fprintf(flog,"# Average temperature of the measurement : %lg\n", TempAve); + fprintf(flog,"# Boltzmann constant under current units : %lg\n", boltz); + fprintf(flog,"# basis vector A1 = [%lg %lg %lg]\n", basevec[0], basevec[1], basevec[2]); + fprintf(flog,"# basis vector A2 = [%lg %lg %lg]\n", basevec[3], basevec[4], basevec[5]); + fprintf(flog,"# basis vector A3 = [%lg %lg %lg]\n", basevec[6], basevec[7], basevec[8]); + fprintf(flog,"############################################################\n"); + fprintf(flog,"# qx\t qy \t qz \t\t Phi(q)\n"); EnforceASR(); diff --git a/src/VORONOI/compute_voronoi_atom.cpp b/src/VORONOI/compute_voronoi_atom.cpp index 5239a27ecc..bf58cdd347 100644 --- a/src/VORONOI/compute_voronoi_atom.cpp +++ b/src/VORONOI/compute_voronoi_atom.cpp @@ -127,6 +127,9 @@ ComputeVoronoi::ComputeVoronoi(LAMMPS *lmp, int narg, char **arg) : if (occupation && ( surface!=VOROSURF_NONE || maxedge>0 ) ) error->all(FLERR,"Illegal compute voronoi/atom command (occupation and (surface or edges))"); + if (occupation && (atom->map_style == 0)) + error->all(FLERR,"Compute voronoi/atom occupation requires an atom map, see atom_modify"); + nmax = rmax = 0; edge = rfield = sendvector = NULL; voro = NULL; diff --git a/src/info.cpp b/src/info.cpp index 175ea7c94f..30c696393f 100644 --- a/src/info.cpp +++ b/src/info.cpp @@ -925,6 +925,18 @@ bool Info::is_available(const char *category, const char *name) delete[] name_w_suffix; } } + } else if (strcmp(category,"feature") == 0) { + if (strcmp(name,"gzip") == 0) { + return has_gzip_support(); + } else if (strcmp(name,"png") == 0) { + return has_png_support(); + } else if (strcmp(name,"jpeg") == 0) { + return has_jpeg_support(); + } else if (strcmp(name,"ffmpeg") == 0) { + return has_ffmpeg_support(); + } else if (strcmp(name,"exceptions") == 0) { + return has_exceptions(); + } } else error->all(FLERR,"Unknown category for info is_available()"); return match ? true : false; @@ -1027,3 +1039,43 @@ static void print_columns(FILE* fp, vector & styles) } } } + +bool Info::has_gzip_support() const { +#ifdef LAMMPS_GZIP + return true; +#else + return false; +#endif +} + +bool Info::has_png_support() const { +#ifdef LAMMPS_PNG + return true; +#else + return false; +#endif +} + +bool Info::has_jpeg_support() const { +#ifdef LAMMPS_JPEG + return true; +#else + return false; +#endif +} + +bool Info::has_ffmpeg_support() const { +#ifdef LAMMPS_FFMPEG + return true; +#else + return false; +#endif +} + +bool Info::has_exceptions() const { +#ifdef LAMMPS_EXCEPTIONS + return true; +#else + return false; +#endif +} diff --git a/src/info.h b/src/info.h index f4549badf7..8e7a96d15a 100644 --- a/src/info.h +++ b/src/info.h @@ -33,6 +33,12 @@ class Info : protected Pointers { bool is_defined(const char *, const char *); bool is_available(const char *, const char *); + bool has_gzip_support() const; + bool has_png_support() const; + bool has_jpeg_support() const; + bool has_ffmpeg_support() const; + bool has_exceptions() const; + private: void available_styles(FILE * out, int flags);
born epsilon,sigma,delta type pairs
lubricate
lj/sf/dipole/sfepsilon,sigma,scaletype pairs
lubricate mu global
gauss
gauss a type pairs
morse
morse d0,r0,alpha type pairs
soft
soft a type pairs
kim
kim PARAM_FREE_*&#58i,j,... global