possible attributes = patom1 patom2 +possible attributes = natom1 natom2 + patom1 patom2 batom1 batom2 btype aatom1 aatom2 aatom3 atype datom1 datom2 datom3 dtype iatom1 iatom2 iatom3 itype-patom1, patom2 = IDs of 2 atoms in each pair +natom1, natom2 = IDs of 2 atoms in each pair (within neighbor cutoff) + patom1, patom2 = IDs of 2 atoms in each pair (within force cutoff) batom1, batom2 = IDs of 2 atoms in each bond btype = bond type of each bond aatom1, aatom2, aatom3 = IDs of 3 atoms in each angle @@ -63,10 +65,13 @@ property/local command.If the inputs are pair attributes, the local data is generated by looping over the pairwise neighbor list. Info about an individual pairwise interaction will only be included if both atoms in the pair -are in the specified compute group, and if the current pairwise -distance is less than the force cutoff distance for that interaction, -as defined by the pair_style and -pair_coeff commands. +are in the specified compute group. For natom1 and natom2, all +atom pairs in the neighbor list are considered (out to the neighbor +cutoff = force cutoff + neighbor skin). For patom1 +and patom2, the distance between the atoms must be less than the +force cutoff distance for that pair to be included, as defined by the +pair_style and pair_coeff +commands.
If the inputs are bond, angle, etc attributes, the local data is generated by looping over all the atoms owned on a processor and @@ -84,9 +89,9 @@ bond/local command can be combined with bond atom indices from this command and output by the dump local command in a consistent way.
-The patom1 and patom2 attributes refer to the atom IDs of the 2 -atoms in each pairwise interaction computed by the -pair_style command. +
The natom1 and natom2, or patom1 and patom2 attributes refer +to the atom IDs of the 2 atoms in each pairwise interaction computed +by the pair_style command.
IMPORTANT NOTE: For pairs, if two atoms I,J are involved in 1-2, 1-3, 1-4 interactions within the molecular topology, their pairwise diff --git a/doc/compute_property_local.txt b/doc/compute_property_local.txt index d0409b31fd..1baaedd069 100644 --- a/doc/compute_property_local.txt +++ b/doc/compute_property_local.txt @@ -15,13 +15,15 @@ compute ID group-ID property/local input1 input2 ... :pre ID, group-ID are documented in "compute"_compute.html command :ulb,l property/local = style name of this compute command :l input = one or more attributes :l - possible attributes = patom1 patom2 + possible attributes = natom1 natom2 + patom1 patom2 batom1 batom2 btype aatom1 aatom2 aatom3 atype datom1 datom2 datom3 dtype iatom1 iatom2 iatom3 itype :pre - patom1, patom2 = IDs of 2 atoms in each pair + natom1, natom2 = IDs of 2 atoms in each pair (within neighbor cutoff) + patom1, patom2 = IDs of 2 atoms in each pair (within force cutoff) batom1, batom2 = IDs of 2 atoms in each bond btype = bond type of each bond aatom1, aatom2, aatom3 = IDs of 3 atoms in each angle @@ -56,10 +58,13 @@ property/local command. If the inputs are pair attributes, the local data is generated by looping over the pairwise neighbor list. Info about an individual pairwise interaction will only be included if both atoms in the pair -are in the specified compute group, and if the current pairwise -distance is less than the force cutoff distance for that interaction, -as defined by the "pair_style"_pair_style.html and -"pair_coeff"_pair_coeff.html commands. +are in the specified compute group. For {natom1} and {natom2}, all +atom pairs in the neighbor list are considered (out to the neighbor +cutoff = force cutoff + "neighbor skin"_neighbor.html). For {patom1} +and {patom2}, the distance between the atoms must be less than the +force cutoff distance for that pair to be included, as defined by the +"pair_style"_pair_style.html and "pair_coeff"_pair_coeff.html +commands. If the inputs are bond, angle, etc attributes, the local data is generated by looping over all the atoms owned on a processor and @@ -77,9 +82,9 @@ bond/local"_compute_bond_local.html command can be combined with bond atom indices from this command and output by the "dump local"_dump.html command in a consistent way. -The {patom1} and {patom2} attributes refer to the atom IDs of the 2 -atoms in each pairwise interaction computed by the -"pair_style"_pair_style.html command. +The {natom1} and {natom2}, or {patom1} and {patom2} attributes refer +to the atom IDs of the 2 atoms in each pairwise interaction computed +by the "pair_style"_pair_style.html command. IMPORTANT NOTE: For pairs, if two atoms I,J are involved in 1-2, 1-3, 1-4 interactions within the molecular topology, their pairwise diff --git a/doc/fix.html b/doc/fix.html index 71c151b381..fc6e8828cc 100644 --- a/doc/fix.html +++ b/doc/fix.html @@ -160,7 +160,8 @@ for further info. arguments and what it does, as listed below. Here is an alphabetic list of fix styles available in LAMMPS:
-
- addforce - add a force to each atom
+
- adapt - change a simulation parameter over time +
- addforce - add a force to each atom
- aveforce - add an averaged force to each atom
- ave/atom - compute per-atom time-averaged quantities
- ave/histo - compute/output time-averaged histograms
diff --git a/doc/fix.txt b/doc/fix.txt
index a9bbf59241..0ec33264f2 100644
--- a/doc/fix.txt
+++ b/doc/fix.txt
@@ -155,6 +155,7 @@ Each fix style has its own documentation page which describes its
arguments and what it does, as listed below. Here is an alphabetic
list of fix styles available in LAMMPS:
+"adapt"_fix_adapt.html - change a simulation parameter over time
"addforce"_fix_addforce.html - add a force to each atom
"aveforce"_fix_aveforce.html - add an averaged force to each atom
"ave/atom"_fix_ave_atom.html - compute per-atom time-averaged quantities
diff --git a/doc/pair_lubricate.html b/doc/pair_lubricate.html
index 6f47436ec4..fbd8aa7fee 100644
--- a/doc/pair_lubricate.html
+++ b/doc/pair_lubricate.html
@@ -43,6 +43,10 @@ their relative velocities in the presence of a background solvent with
viscosity mu. Note that this is dynamic viscosity which has units of
mass/distance/time, not kinematic viscosity.
+
The viscosity mu can be varied in a time-dependent manner over the +course of a simluation. See the fix adapt command +for details. +
Rc is the outer cutoff specified in the pair_style command, the translational velocities of the 2 particles are v1 and v2, the angular velocities are w1 and w2, and n is the unit vector in the direction diff --git a/doc/pair_lubricate.txt b/doc/pair_lubricate.txt index aa1caae94b..b793a7f546 100644 --- a/doc/pair_lubricate.txt +++ b/doc/pair_lubricate.txt @@ -40,6 +40,10 @@ their relative velocities in the presence of a background solvent with viscosity mu. Note that this is dynamic viscosity which has units of mass/distance/time, not kinematic viscosity. +The viscosity mu can be varied in a time-dependent manner over the +course of a simluation. See the "fix adapt"_fix_adapt.html command +for details. + Rc is the outer cutoff specified in the pair_style command, the translational velocities of the 2 particles are v1 and v2, the angular velocities are w1 and w2, and n is the unit vector in the direction diff --git a/doc/pair_soft.html b/doc/pair_soft.html index bcacfea72f..d2aab51c72 100644 --- a/doc/pair_soft.html +++ b/doc/pair_soft.html @@ -20,8 +20,8 @@
Examples:
pair_style soft 2.5 -pair_coeff * * 0.0 60.0 -pair_coeff 1 1 0.0 60.0 3.0 +pair_coeff * * 10.0 +pair_coeff 1 1 10.0 3.0
Description:
@@ -30,11 +30,12 @@ pair_coeff 1 1 0.0 60.0 3.0
It is useful for pushing apart overlapping atoms, since it does not -blow up as r goes to 0. A is a pre-factor that varies in time from -the start to the end of the run. The run command documents -how to make the ramping take place across multiple runs. Rc is the -cutoff. See the fix nve/limit command for -another way to push apart overlapping atoms. +blow up as r goes to 0. A is a pre-factor that can be made to vary in +time from the start to the end of the run (see discussion below), +e.g. to start with a very soft potential and slowly harden the +interactions over time. Rc is the cutoff. See the fix +nve/limit command for another way to push apart +overlapping atoms.
The following coefficients must be defined for each pair of atom types via the pair_coeff command as in the examples above, @@ -42,28 +43,34 @@ or in the data file or restart files read by the read_data or read_restart commands, or by mixing as described below:
-- Astart (energy units) -
- Astop (energy units)
+
- A (energy units)
- cutoff (distance units)
Astart and Astop are the values of the prefactor at the start and end -of the next run. At intermediate times the value of A will be ramped -between these 2 values. Note that before performing a 2nd run, you -will want to adjust the values of Astart and Astop for all type pairs, -or switch to a new pair style. -
The last coefficient is optional. If not specified, the global soft cutoff is used.
+The fix adapt command can be used to vary A over the +course of a simulation. For example these commands will vary the +prefactor A for all pairwise interactions from 0.0 at the beginning to +30.0 at the end of a 10,000 step run: +
+variable prefactor equal 30.0*elapsed/10000 +fix 1 all adapt 1 pair soft a * * prefactor +
+Note that a formula defined by an equal-style variable +can use the current timestep, elapsed time in the current run, elapsed +time since the beginning of a series of runs, as well as access other +variables. +
Mixing, shift, table, tail correction, restart, rRESPA info:
-For atom type pairs I,J and I != J, the Astart, Astop coefficients and -cutoff distance for this pair style can be mixed. Astart and Atop are -always mixed via a geometric rule. The cutoff is mixed according to -the pair_modify mix value. The default mix value is geometric. See -the "pair_modify" command for details. +
For atom type pairs I,J and I != J, the A coefficient and cutoff +distance for this pair style can be mixed. A is always mixed via a +geometric rule. The cutoff is mixed according to the pair_modify +mix value. The default mix value is geometric. See the +"pair_modify" command for details.
This pair style does not support the pair_modify shift option, since the pair interaction goes to 0.0 at the cutoff. @@ -85,7 +92,8 @@ to be specified in an input script that reads a restart file.
Related commands:
-pair_coeff, fix nve/limit, fix +adapt
Default: none
diff --git a/doc/pair_soft.txt b/doc/pair_soft.txt index a5be96d086..86425b648f 100644 --- a/doc/pair_soft.txt +++ b/doc/pair_soft.txt @@ -17,8 +17,8 @@ cutoff = global cutoff for soft interactions (distance units) :ul [Examples:] pair_style soft 2.5 -pair_coeff * * 0.0 60.0 -pair_coeff 1 1 0.0 60.0 3.0 :pre +pair_coeff * * 10.0 +pair_coeff 1 1 10.0 3.0 :pre [Description:] @@ -27,11 +27,12 @@ Style {soft} computes pairwise interactions with the formula :c,image(Eqs/pair_soft.jpg) It is useful for pushing apart overlapping atoms, since it does not -blow up as r goes to 0. A is a pre-factor that varies in time from -the start to the end of the run. The "run"_run.html command documents -how to make the ramping take place across multiple runs. Rc is the -cutoff. See the "fix nve/limit"_fix_nve_limit.html command for -another way to push apart overlapping atoms. +blow up as r goes to 0. A is a pre-factor that can be made to vary in +time from the start to the end of the run (see discussion below), +e.g. to start with a very soft potential and slowly harden the +interactions over time. Rc is the cutoff. See the "fix +nve/limit"_fix_nve_limit.html command for another way to push apart +overlapping atoms. The following coefficients must be defined for each pair of atom types via the "pair_coeff"_pair_coeff.html command as in the examples above, @@ -39,28 +40,34 @@ or in the data file or restart files read by the "read_data"_read_data.html or "read_restart"_read_restart.html commands, or by mixing as described below: -Astart (energy units) -Astop (energy units) +A (energy units) cutoff (distance units) :ul -Astart and Astop are the values of the prefactor at the start and end -of the next run. At intermediate times the value of A will be ramped -between these 2 values. Note that before performing a 2nd run, you -will want to adjust the values of Astart and Astop for all type pairs, -or switch to a new pair style. - The last coefficient is optional. If not specified, the global soft cutoff is used. +The "fix adapt"_fix_adapt.html command can be used to vary A over the +course of a simulation. For example these commands will vary the +prefactor A for all pairwise interactions from 0.0 at the beginning to +30.0 at the end of a 10,000 step run: + +variable prefactor equal 30.0*elapsed/10000 +fix 1 all adapt 1 pair soft a * * prefactor :pre + +Note that a formula defined by an "equal-style variable"_variable.html +can use the current timestep, elapsed time in the current run, elapsed +time since the beginning of a series of runs, as well as access other +variables. + :line [Mixing, shift, table, tail correction, restart, rRESPA info]: -For atom type pairs I,J and I != J, the Astart, Astop coefficients and -cutoff distance for this pair style can be mixed. Astart and Atop are -always mixed via a {geometric} rule. The cutoff is mixed according to -the pair_modify mix value. The default mix value is {geometric}. See -the "pair_modify" command for details. +For atom type pairs I,J and I != J, the A coefficient and cutoff +distance for this pair style can be mixed. A is always mixed via a +{geometric} rule. The cutoff is mixed according to the pair_modify +mix value. The default mix value is {geometric}. See the +"pair_modify" command for details. This pair style does not support the "pair_modify"_pair_modify.html shift option, since the pair interaction goes to 0.0 at the cutoff. @@ -82,6 +89,7 @@ This pair style can only be used via the {pair} keyword of the [Related commands:] -"pair_coeff"_pair_coeff.html, "fix nve/limit"_fix_nve_limit.html +"pair_coeff"_pair_coeff.html, "fix nve/limit"_fix_nve_limit.html, "fix +adapt"_fix_adapt.html [Default:] none diff --git a/doc/run.html b/doc/run.html index 05524ee43b..a93d102276 100644 --- a/doc/run.html +++ b/doc/run.html @@ -69,8 +69,7 @@ performed and you want a fix command that changes some value over time (e.g. temperature) to make the change across the entire set of runs and not just a single run. See the doc page for individual fixes to see which ones can be used with the start/stop -keywords. The pair_style soft potential also -changes its pair potential coefficients in this manner. +keywords.For example, consider this fix followed by 10 run commands:
diff --git a/doc/run.txt b/doc/run.txt index c80f308dc6..56f47be800 100644 --- a/doc/run.txt +++ b/doc/run.txt @@ -62,8 +62,7 @@ performed and you want a "fix"_fix.html command that changes some value over time (e.g. temperature) to make the change across the entire set of runs and not just a single run. See the doc page for individual fixes to see which ones can be used with the {start/stop} -keywords. The "pair_style soft"_pair_style.html potential also -changes its pair potential coefficients in this manner. +keywords. For example, consider this fix followed by 10 run commands: diff --git a/examples/micelle/in.micelle b/examples/micelle/in.micelle index 128e4539ef..841e55ab2c 100644 --- a/examples/micelle/in.micelle +++ b/examples/micelle/in.micelle @@ -13,20 +13,25 @@ read_data data.micelle special_bonds fene pair_style soft 1.12246 -pair_coeff * * 1.0 20.0 1.12246 +pair_coeff * * 0.0 1.12246 bond_style harmonic bond_coeff 1 50.0 0.75 velocity all create 0.45 2349852 +variable prefactor equal 1.0+elapsed*(20.0-1.0)/1000 + fix 1 all nve fix 2 all temp/rescale 100 0.45 0.45 0.02 1.0 -fix 3 all enforce2d +fix 3 all adapt 1 pair soft a * * prefactor +fix 4 all enforce2d thermo 50 run 1000 +unfix 3 + # Main run pair_style lj/cut 2.5 diff --git a/src/COLLOID/pair_colloid.cpp b/src/COLLOID/pair_colloid.cpp index 09dbab7de7..ef9e8aa83b 100644 --- a/src/COLLOID/pair_colloid.cpp +++ b/src/COLLOID/pair_colloid.cpp @@ -471,7 +471,7 @@ double PairColloid::single(int i, int j, int itype, int jtype, double rsq, c2 = a2[itype][jtype]; K[1] = c2*c2; K[2] = rsq; - K[0] = K[1] - rsq; + K[0] = K[1] - rsq; K[4] = rsq*rsq; K[3] = K[1] - K[2]; K[3] *= K[3]*K[3]; diff --git a/src/COLLOID/pair_lubricate.cpp b/src/COLLOID/pair_lubricate.cpp index b17be877ec..69a207430b 100644 --- a/src/COLLOID/pair_lubricate.cpp +++ b/src/COLLOID/pair_lubricate.cpp @@ -533,3 +533,27 @@ void PairLubricate::read_restart_settings(FILE *fp) delete random; random = new RanMars(lmp,seed + comm->me); } + +/* ---------------------------------------------------------------------- + check if name is recognized, return integer index for that name + if name not recognized, return -1 + if type pair setting, return -2 if no type pairs are set +------------------------------------------------------------------------- */ + +int PairLubricate::pre_adapt(char *name, int ilo, int ihi, int jlo, int jhi) +{ + if (strcmp(name,"mu") == 0) return 0; + return -1; +} + +/* ---------------------------------------------------------------------- + adapt parameter indexed by which + change all pair variables affected by the reset parameter + if type pair setting, set I-J and J-I coeffs +------------------------------------------------------------------------- */ + +void PairLubricate::adapt(int which, int ilo, int ihi, int jlo, int jhi, + double value) +{ + mu = value; +} diff --git a/src/COLLOID/pair_lubricate.h b/src/COLLOID/pair_lubricate.h index b4e76c2b6a..5206c6d5fc 100644 --- a/src/COLLOID/pair_lubricate.h +++ b/src/COLLOID/pair_lubricate.h @@ -37,6 +37,8 @@ class PairLubricate : public Pair { void read_restart(FILE *); void write_restart_settings(FILE *); void read_restart_settings(FILE *); + int pre_adapt(char *, int, int, int, int); + void adapt(int, int, int, int, int, double); protected: double cut_inner_global,cut_global; diff --git a/src/MOLECULE/bond_table.cpp b/src/MOLECULE/bond_table.cpp index d404cede0e..28e28f17f2 100644 --- a/src/MOLECULE/bond_table.cpp +++ b/src/MOLECULE/bond_table.cpp @@ -375,7 +375,7 @@ void BondTable::compute_table(Table *tb) // delta = table spacing for N-1 bins int tlm1 = tablength-1; - tb->delta = (tb->hi - tb->lo)/ nm1; + tb->delta = (tb->hi - tb->lo)/ tlm1; tb->invdelta = 1.0/tb->delta; tb->deltasq6 = tb->delta*tb->delta / 6.0; diff --git a/src/Makefile b/src/Makefile index d172e4dc14..3d4843a1bd 100755 --- a/src/Makefile +++ b/src/Makefile @@ -14,7 +14,7 @@ OBJ = $(SRC:.cpp=.o) # Package variables PACKAGE = asphere class2 colloid dipole dsmc gpu granular \ - kspace manybody meam molecule opt peri poems prd reax xtc + kspace manybody meam molecule opt peri poems prd reax shock xtc PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm user-ewaldn \ user-imd user-smd diff --git a/src/compute_dihedral_local.cpp b/src/compute_dihedral_local.cpp index 1aa41e9a01..a8f2eddc50 100644 --- a/src/compute_dihedral_local.cpp +++ b/src/compute_dihedral_local.cpp @@ -109,9 +109,8 @@ int ComputeDihedralLocal::compute_dihedrals(int flag) { int i,m,n,atom1,atom2,atom3,atom4; double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; - double sb1,sb2,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2; - double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2; - double c2mag,sin2,sc1,sc2,s1,s2,s12,c; + double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv; + double s,c; double *pbuf; double **x = atom->x; @@ -148,71 +147,53 @@ int ComputeDihedralLocal::compute_dihedrals(int flag) if (flag) { - // phi calculation from dihedral style OPLS + // phi calculation from dihedral style harmonic if (pflag >= 0) { vb1x = x[atom1][0] - x[atom2][0]; vb1y = x[atom1][1] - x[atom2][1]; vb1z = x[atom1][2] - x[atom2][2]; domain->minimum_image(vb1x,vb1y,vb1z); - + vb2x = x[atom3][0] - x[atom2][0]; vb2y = x[atom3][1] - x[atom2][1]; vb2z = x[atom3][2] - x[atom2][2]; domain->minimum_image(vb2x,vb2y,vb2z); - + vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; domain->minimum_image(vb2xm,vb2ym,vb2zm); - + vb3x = x[atom4][0] - x[atom3][0]; vb3y = x[atom4][1] - x[atom3][1]; vb3z = x[atom4][2] - x[atom3][2]; domain->minimum_image(vb3x,vb3y,vb3z); - - sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); - sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); - sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); - - rb1 = sqrt(sb1); - rb3 = sqrt(sb3); - - c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3; - - b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; - b1mag = sqrt(b1mag2); - b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; - b2mag = sqrt(b2mag2); - b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; - b3mag = sqrt(b3mag2); - ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; - r12c1 = 1.0 / (b1mag*b2mag); - c1mag = ctmp * r12c1; + ax = vb1y*vb2zm - vb1z*vb2ym; + ay = vb1z*vb2xm - vb1x*vb2zm; + az = vb1x*vb2ym - vb1y*vb2xm; + bx = vb3y*vb2zm - vb3z*vb2ym; + by = vb3z*vb2xm - vb3x*vb2zm; + bz = vb3x*vb2ym - vb3y*vb2xm; - ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; - r12c2 = 1.0 / (b2mag*b3mag); - c2mag = ctmp * r12c2; + rasq = ax*ax + ay*ay + az*az; + rbsq = bx*bx + by*by + bz*bz; + rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm; + rg = sqrt(rgsq); + + rginv = ra2inv = rb2inv = 0.0; + if (rg > 0) rginv = 1.0/rg; + if (rasq > 0) ra2inv = 1.0/rasq; + if (rbsq > 0) rb2inv = 1.0/rbsq; + rabinv = sqrt(ra2inv*rb2inv); - sin2 = MAX(1.0 - c1mag*c1mag,0.0); - sc1 = sqrt(sin2); - if (sc1 < SMALL) sc1 = SMALL; - sc1 = 1.0/sc1; - - sin2 = MAX(1.0 - c2mag*c2mag,0.0); - sc2 = sqrt(sin2); - if (sc2 < SMALL) sc2 = SMALL; - sc2 = 1.0/sc2; - - s1 = sc1 * sc1; - s2 = sc2 * sc2; - s12 = sc1 * sc2; - c = (c0 + c1mag*c2mag) * s12; + c = (ax*bx + ay*by + az*bz)*rabinv; + s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z); if (c > 1.0) c = 1.0; if (c < -1.0) c = -1.0; - pbuf[n] = 180.0*acos(c)/PI; + pbuf[n] = 180.0*atan2(s,c)/PI; } n += nvalues; } diff --git a/src/compute_erotate_sphere.cpp b/src/compute_erotate_sphere.cpp index 2ddcb65641..1558196cce 100644 --- a/src/compute_erotate_sphere.cpp +++ b/src/compute_erotate_sphere.cpp @@ -109,6 +109,7 @@ double ComputeERotateSphere::compute_scalar() radius[i]*radius[i]*mass[itype]; } } + } else { if (rmass) { for (i = 0; i < nlocal; i++) diff --git a/src/compute_property_local.cpp b/src/compute_property_local.cpp index bf10dc7101..d7c814eafb 100644 --- a/src/compute_property_local.cpp +++ b/src/compute_property_local.cpp @@ -26,7 +26,7 @@ using namespace LAMMPS_NS; -enum{NONE,PAIR,BOND,ANGLE,DIHEDRAL,IMPROPER}; +enum{NONE,NEIGH,PAIR,BOND,ANGLE,DIHEDRAL,IMPROPER}; #define DELTA 10000 @@ -50,7 +50,18 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) : for (int iarg = 3; iarg < narg; iarg++) { i = iarg-3; - if (strcmp(arg[iarg],"patom1") == 0) { + if (strcmp(arg[iarg],"natom1") == 0) { + pack_choice[i] = &ComputePropertyLocal::pack_patom1; + if (kindflag != NONE && kindflag != NEIGH) + error->all("Compute property/local cannot use these inputs together"); + kindflag = NEIGH; + } else if (strcmp(arg[iarg],"natom2") == 0) { + pack_choice[i] = &ComputePropertyLocal::pack_patom2; + if (kindflag != NONE && kindflag != NEIGH) + error->all("Compute property/local cannot use these inputs together"); + kindflag = NEIGH; + + } else if (strcmp(arg[iarg],"patom1") == 0) { pack_choice[i] = &ComputePropertyLocal::pack_patom1; if (kindflag != NONE && kindflag != PAIR) error->all("Compute property/local cannot use these inputs together"); @@ -184,16 +195,16 @@ ComputePropertyLocal::~ComputePropertyLocal() void ComputePropertyLocal::init() { - if (kindflag == PAIR) { + if (kindflag == NEIGH || kindflag == PAIR) { if (force->pair == NULL) error->all("No pair style is defined for compute property/local"); if (force->pair->single_enable == 0) error->all("Pair style does not support compute property/local"); } - // for PAIR< need an occasional half neighbor list + // for NEIGH/PAIR need an occasional half neighbor list - if (kindflag == PAIR) { + if (kindflag == NEIGH || kindflag == PAIR) { int irequest = neighbor->request((void *) this); neighbor->requests[irequest]->pair = 0; neighbor->requests[irequest]->compute = 1; @@ -201,9 +212,10 @@ void ComputePropertyLocal::init() } // do initial memory allocation so that memory_usage() is correct - // cannot be done yet for PAIR, since neigh list does not exist + // cannot be done yet for NEIGH/PAIR, since neigh list does not exist - if (kindflag == PAIR) ncount = 0; + if (kindflag == NEIGH) ncount = 0; + else if (kindflag == PAIR) ncount = 0; else if (kindflag == BOND) ncount = count_bonds(0); else if (kindflag == ANGLE) ncount = count_angles(0); else if (kindflag == DIHEDRAL) ncount = count_dihedrals(0); @@ -228,7 +240,8 @@ void ComputePropertyLocal::compute_local() // count local entries and generate list of indices - if (kindflag == PAIR) ncount = count_pairs(0); + if (kindflag == NEIGH) ncount = count_pairs(0,0); + else if (kindflag == PAIR) ncount = count_pairs(0,1); else if (kindflag == BOND) ncount = count_bonds(0); else if (kindflag == ANGLE) ncount = count_angles(0); else if (kindflag == DIHEDRAL) ncount = count_dihedrals(0); @@ -237,7 +250,8 @@ void ComputePropertyLocal::compute_local() if (ncount > nmax) reallocate(ncount); size_local_rows = ncount; - if (kindflag == PAIR) ncount = count_pairs(1); + if (kindflag == NEIGH) ncount = count_pairs(1,0); + else if (kindflag == PAIR) ncount = count_pairs(1,1); else if (kindflag == BOND) ncount = count_bonds(1); else if (kindflag == ANGLE) ncount = count_angles(1); else if (kindflag == DIHEDRAL) ncount = count_dihedrals(1); @@ -259,10 +273,11 @@ void ComputePropertyLocal::compute_local() count pairs and compute pair info on this proc only count pair once if newton_pair is off both atom I,J must be in group - if flag is set, compute requested info about pair + if allflag is set, compute requested info about pair + if forceflag = 1, pair must be within force cutoff, else neighbor cutoff ------------------------------------------------------------------------- */ -int ComputePropertyLocal::count_pairs(int flag) +int ComputePropertyLocal::count_pairs(int allflag, int forceflag) { int i,j,m,n,ii,jj,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz; @@ -281,7 +296,7 @@ int ComputePropertyLocal::count_pairs(int flag) // invoke half neighbor list (will copy or build if necessary) - if (flag == 0) neighbor->build_one(list->index); + if (allflag == 0) neighbor->build_one(list->index); inum = list->inum; ilist = list->ilist; @@ -323,9 +338,9 @@ int ComputePropertyLocal::count_pairs(int flag) delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; jtype = type[j]; - if (rsq >= cutsq[itype][jtype]) continue; + if (forceflag && rsq >= cutsq[itype][jtype]) continue; - if (flag) { + if (allflag) { indices[m][0] = tag[i]; indices[m][1] = tag[j]; } diff --git a/src/compute_property_local.h b/src/compute_property_local.h index 59935acfe1..7517c674bc 100644 --- a/src/compute_property_local.h +++ b/src/compute_property_local.h @@ -46,7 +46,7 @@ class ComputePropertyLocal : public Compute { int ncount; int **indices; - int count_pairs(int); + int count_pairs(int, int); int count_bonds(int); int count_angles(int); int count_dihedrals(int); diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp index 6fc6f94f8c..32ff36a97a 100644 --- a/src/compute_reduce.cpp +++ b/src/compute_reduce.cpp @@ -488,8 +488,8 @@ double ComputeReduce::compute_one(int m, int flag) // only include atoms in group for atom properties and per-atom quantities index = -1; - int v_idx = value2index[m]; - int a_idx = argindex[m]; + int vidx = value2index[m]; + int aidx = argindex[m]; int *mask = atom->mask; int nlocal = atom->nlocal; @@ -504,28 +504,28 @@ double ComputeReduce::compute_one(int m, int flag) double **x = atom->x; if (flag < 0) { for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) combine(one,x[i][a_idx],i); - } else one = x[flag][a_idx]; + if (mask[i] & groupbit) combine(one,x[i][aidx],i); + } else one = x[flag][aidx]; } else if (which[m] == V) { double **v = atom->v; if (flag < 0) { for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) combine(one,v[i][a_idx],i); - } else one = v[flag][a_idx]; + if (mask[i] & groupbit) combine(one,v[i][aidx],i); + } else one = v[flag][aidx]; } else if (which[m] == F) { double **f = atom->f; if (flag < 0) { for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) combine(one,f[i][a_idx],i); - } else one = f[flag][a_idx]; + if (mask[i] & groupbit) combine(one,f[i][aidx],i); + } else one = f[flag][aidx]; // invoke compute if not previously invoked } else if (which[m] == COMPUTE) { - Compute *compute = modify->compute[v_idx]; + Compute *compute = modify->compute[vidx]; if (flavor[m] == GLOBAL) { - if (a_idx == 0) { + if (aidx == 0) { if (!(compute->invoked_flag & INVOKED_VECTOR)) { compute->compute_vector(); compute->invoked_flag |= INVOKED_VECTOR; @@ -541,13 +541,13 @@ double ComputeReduce::compute_one(int m, int flag) compute->compute_array(); compute->invoked_flag |= INVOKED_ARRAY; } - double **comp_arr = compute->array; + double **carray = compute->array; int n = compute->size_array_rows; - int a_idx_1 = a_idx - 1; + int aidx_1 = aidx - 1; if (flag < 0) for (i = 0; i < n; i++) - combine(one,comp_arr[i][a_idx_1],i); - else one = comp_arr[flag][a_idx_1]; + combine(one,carray[i][aidx_1],i); + else one = carray[flag][aidx_1]; } } else if (flavor[m] == PERATOM) { @@ -556,7 +556,7 @@ double ComputeReduce::compute_one(int m, int flag) compute->invoked_flag |= INVOKED_PERATOM; } - if (a_idx == 0) { + if (aidx == 0) { double *comp_vec = compute->vector_atom; int n = nlocal; if (flag < 0) { @@ -564,13 +564,13 @@ double ComputeReduce::compute_one(int m, int flag) if (mask[i] & groupbit) combine(one,comp_vec[i],i); } else one = comp_vec[flag]; } else { - double **comp_arr = compute->array_atom; + double **carray_atom = compute->array_atom; int n = nlocal; - int a_idxm1 = a_idx - 1; + int aidxm1 = aidx - 1; if (flag < 0) { for (i = 0; i < n; i++) - if (mask[i] & groupbit) combine(one,comp_arr[i][a_idxm1],i); - } else one = comp_arr[flag][a_idxm1]; + if (mask[i] & groupbit) combine(one,carray_atom[i][aidxm1],i); + } else one = carray_atom[flag][aidxm1]; } } else if (flavor[m] == LOCAL) { @@ -579,7 +579,7 @@ double ComputeReduce::compute_one(int m, int flag) compute->invoked_flag |= INVOKED_LOCAL; } - if (a_idx == 0) { + if (aidx == 0) { double *comp_vec = compute->vector_local; int n = compute->size_local_rows; if (flag < 0) @@ -587,40 +587,40 @@ double ComputeReduce::compute_one(int m, int flag) combine(one,comp_vec[i],i); else one = comp_vec[flag]; } else { - double **comp_arr = compute->array_local; + double **carray_local = compute->array_local; int n = compute->size_local_rows; - int a_idxm1 = a_idx - 1; + int aidxm1 = aidx - 1; if (flag < 0) for (i = 0; i < n; i++) - combine(one,comp_arr[i][a_idxm1],i); - else one = comp_arr[flag][a_idxm1]; + combine(one,carray_local[i][aidxm1],i); + else one = carray_local[flag][aidxm1]; } } // access fix fields, check if fix frequency is a match } else if (which[m] == FIX) { - if (update->ntimestep % modify->fix[v_idx]->peratom_freq) + if (update->ntimestep % modify->fix[vidx]->peratom_freq) error->all("Fix used in compute reduce not computed at compatible time"); - Fix *fix = modify->fix[v_idx]; + Fix *fix = modify->fix[vidx]; if (flavor[m] == GLOBAL) { - if (a_idx == 0) { + if (aidx == 0) { int n = fix->size_vector; if (flag < 0) for (i = 0; i < n; i++) combine(one,fix->compute_vector(i),i); else one = fix->compute_vector(flag); } else { - int a_idxm1 = a_idx - 1; + int aidxm1 = aidx - 1; if (flag < 0) for (i = 0; i < nlocal; i++) - combine(one,fix->compute_array(i,a_idxm1),i); - else one = fix->compute_array(flag,a_idxm1); + combine(one,fix->compute_array(i,aidxm1),i); + else one = fix->compute_array(flag,aidxm1); } } else if (flavor[m] == PERATOM) { - if (a_idx == 0) { + if (aidx == 0) { double *fix_vector = fix->vector_atom; int n = nlocal; if (flag < 0) { @@ -629,15 +629,15 @@ double ComputeReduce::compute_one(int m, int flag) } else one = fix_vector[flag]; } else { double **fix_array = fix->array_atom; - int a_idxm1 = a_idx - 1; + int aidxm1 = aidx - 1; if (flag < 0) { for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) combine(one,fix_array[i][a_idxm1],i); - } else one = fix_array[flag][a_idxm1]; + if (mask[i] & groupbit) combine(one,fix_array[i][aidxm1],i); + } else one = fix_array[flag][aidxm1]; } } else if (flavor[m] == LOCAL) { - if (a_idx == 0) { + if (aidx == 0) { double *fix_vector = fix->vector_local; int n = fix->size_local_rows; if (flag < 0) @@ -647,11 +647,11 @@ double ComputeReduce::compute_one(int m, int flag) } else { double **fix_array = fix->array_local; int n = fix->size_local_rows; - int a_idxm1 = a_idx - 1; + int aidxm1 = aidx - 1; if (flag < 0) for (i = 0; i < n; i++) - combine(one,fix_array[i][a_idxm1],i); - else one = fix_array[flag][a_idxm1]; + combine(one,fix_array[i][aidxm1],i); + else one = fix_array[flag][aidxm1]; } } @@ -665,7 +665,7 @@ double ComputeReduce::compute_one(int m, int flag) memory->smalloc(maxatom*sizeof(double),"reduce:varatom"); } - input->variable->compute_atom(v_idx,igroup,varatom,1,0); + input->variable->compute_atom(vidx,igroup,varatom,1,0); if (flag < 0) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) combine(one,varatom[i],i); @@ -679,15 +679,15 @@ double ComputeReduce::compute_one(int m, int flag) double ComputeReduce::count(int m) { - int v_idx = value2index[m]; - int a_idx = argindex[m]; + int vidx = value2index[m]; + int aidx = argindex[m]; if (which[m] == X || which[m] == V || which[m] == F) return group->count(igroup); else if (which[m] == COMPUTE) { - Compute *compute = modify->compute[v_idx]; + Compute *compute = modify->compute[vidx]; if (flavor[m] == GLOBAL) { - if (a_idx == 0) return(1.0*compute->size_vector); + if (aidx == 0) return(1.0*compute->size_vector); else return(1.0*compute->size_array_rows); } else if (flavor[m] == PERATOM) { return group->count(igroup); @@ -698,9 +698,9 @@ double ComputeReduce::count(int m) return ncountall; } } else if (which[m] == FIX) { - Fix *fix = modify->fix[v_idx]; + Fix *fix = modify->fix[vidx]; if (flavor[m] == GLOBAL) { - if (a_idx == 0) return(1.0*fix->size_vector); + if (aidx == 0) return(1.0*fix->size_vector); else return(1.0*fix->size_array_rows); } else if (flavor[m] == PERATOM) { return group->count(igroup); diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp new file mode 100644 index 0000000000..f9a5ba5899 --- /dev/null +++ b/src/fix_adapt.cpp @@ -0,0 +1,212 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "math.h" +#include "string.h" +#include "stdlib.h" +#include "fix_adapt.h" +#include "atom.h" +#include "force.h" +#include "pair.h" +#include "input.h" +#include "variable.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{PAIR,ATOM}; +enum{DIAMETER}; + +/* ---------------------------------------------------------------------- */ + +FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) +{ + if (narg < 4) error->all("Illegal fix adapt command"); + nevery = atoi(arg[3]); + if (nevery <= 0) error->all("Illegal fix adapt command"); + + // count # of adaptations + + nadapt = 0; + + int iarg = 4; + while (iarg < narg) { + if (strcmp(arg[iarg],"pair") == 0) { + if (iarg+6 > narg) error->all("Illegal fix adapt command"); + nadapt++; + iarg += 6; + } else if (strcmp(arg[iarg],"atom") == 0) { + if (iarg+3 > narg) error->all("Illegal fix adapt command"); + nadapt++; + iarg += 3; + } else error->all("Illegal fix adapt command"); + } + + // allocate per-adapt vectors + + which = new int[nadapt]; + pair = new char*[nadapt]; + param = new char*[nadapt]; + ilo = new int[nadapt]; + ihi = new int[nadapt]; + jlo = new int[nadapt]; + jhi = new int[nadapt]; + var = new char*[nadapt]; + ivar = new int[nadapt]; + pairptr = new Pair*[nadapt]; + pairindex = new int[nadapt]; + awhich = new int[nadapt]; + + // parse keywords + + nadapt = 0; + + iarg = 4; + while (iarg < narg) { + if (strcmp(arg[iarg],"pair") == 0) { + if (iarg+6 > narg) error->all("Illegal fix adapt command"); + which[nadapt] = PAIR; + int n = strlen(arg[iarg+1]) + 1; + pair[nadapt] = new char[n]; + strcpy(pair[nadapt],arg[iarg+1]); + n = strlen(arg[iarg+2]) + 1; + param[nadapt] = new char[n]; + strcpy(param[nadapt],arg[iarg+2]); + force->bounds(arg[iarg+3],atom->ntypes,ilo[nadapt],ihi[nadapt]); + force->bounds(arg[iarg+4],atom->ntypes,jlo[nadapt],jhi[nadapt]); + n = strlen(arg[iarg+5]) + 1; + var[nadapt] = new char[n]; + strcpy(var[nadapt],arg[iarg+5]); + nadapt++; + iarg += 6; + } else if (strcmp(arg[iarg],"atom") == 0) { + if (iarg+3 > narg) error->all("Illegal fix adapt command"); + which[nadapt] = ATOM; + int n = strlen(arg[iarg+1]) + 1; + param[nadapt] = new char[n]; + strcpy(param[nadapt],arg[iarg+1]); + n = strlen(arg[iarg+2]) + 1; + var[nadapt] = new char[n]; + strcpy(var[nadapt],arg[iarg+2]); + nadapt++; + iarg += 3; + } else error->all("Illegal fix adapt command"); + } +} + +/* ---------------------------------------------------------------------- */ + +FixAdapt::~FixAdapt() +{ + for (int i = 0; i < nadapt; i++) { + if (which[i] == PAIR) delete [] pair[i]; + delete [] param[i]; + delete [] var[i]; + } + delete [] which; + delete [] pair; + delete [] param; + delete [] ilo; + delete [] ihi; + delete [] jlo; + delete [] jhi; + delete [] var; + delete [] ivar; + delete [] pairptr; + delete [] pairindex; + delete [] awhich; +} + +/* ---------------------------------------------------------------------- */ + +int FixAdapt::setmask() +{ + int mask = 0; + mask |= PRE_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixAdapt::init() +{ + // error checks + + for (int m = 0; m < nadapt; m++) { + if (which[m] == PAIR) { + pairptr[m] = force->pair_match(pair[m],1); + if (pairptr[m] == NULL) + error->all("Fix adapt pair style does not exist"); + pairindex[m] = + pairptr[m]->pre_adapt(param[m],ilo[m],ihi[m],jlo[m],jhi[m]); + if (pairindex[m] == -1) + error->all("Fix adapt pair parameter is not recognized"); + if (pairindex[m] == -2) + error->all("Fix adapt pair types are not valid"); + + } else if (which[m] == ATOM) { + if (strcmp(param[m],"diameter") == 0) { + awhich[m] = DIAMETER; + if (!atom->radius_flag) + error->all("Fix adapt requires atom attribute radius"); + } else error->all("Fix adapt atom attribute is not recognized"); + } + + ivar[m] = input->variable->find(var[m]); + if (ivar[m] < 0) error->all("Variable name for fix adapt does not exist"); + if (!input->variable->equalstyle(ivar[m])) + error->all("Variable for fix adapt is not equal style"); + } + + // set params to values for initial force calculation + // needs to happen here in init() instead of setup() + // because modify->setup() is called after pre-Verlet forces are computed + + pre_force(0); +} + +/* ---------------------------------------------------------------------- */ + +void FixAdapt::pre_force(int vflag) +{ + for (int m = 0; m < nadapt; m++) { + double value = input->variable->compute_equal(ivar[m]); + + if (which[m] == PAIR) + pairptr[m]->adapt(pairindex[m],ilo[m],ihi[m],jlo[m],jhi[m],value); + + else if (which[m] == ATOM) { + + // set radius from diameter + // set rmass if both rmass and density are defined + + if (awhich[m] == DIAMETER) { + int mflag = 0; + if (atom->rmass_flag && atom->density_flag) mflag = 1; + double PI = 4.0*atan(1.0); + + double *radius = atom->radius; + double *rmass = atom->rmass; + double *density = atom->density; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + radius[i] = 0.5*value; + rmass[i] = 4.0*PI/3.0 * radius[i]*radius[i]*radius[i] * density[i]; + } + } + } + } +} diff --git a/src/fix_adapt.h b/src/fix_adapt.h new file mode 100644 index 0000000000..1411708e71 --- /dev/null +++ b/src/fix_adapt.h @@ -0,0 +1,50 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(adapt,FixAdapt) + +#else + +#ifndef LMP_FIX_ADAPT_H +#define LMP_FIX_ADAPT_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixAdapt : public Fix { + public: + FixAdapt(class LAMMPS *, int, char **); + ~FixAdapt(); + int setmask(); + void init(); + void pre_force(int); + + private: + int nadapt; + int *which; + char **pair,**param,**var; + int *ilo,*ihi,*jlo,*jhi; + + int *ivar; + class Pair **pairptr; + int *pairindex; + int *awhich; +}; + +} + +#endif +#endif diff --git a/src/fix_aveforce.cpp b/src/fix_aveforce.cpp index 913b82a4d2..7abece0d37 100644 --- a/src/fix_aveforce.cpp +++ b/src/fix_aveforce.cpp @@ -40,6 +40,8 @@ FixAveForce::FixAveForce(LAMMPS *lmp, int narg, char **arg) : global_freq = 1; extvector = 1; + xstr = ystr = zstr = NULL; + if (strstr(arg[3],"v_") == arg[3]) { int n = strlen(&arg[3][2]) + 1; xstr = new char[n]; diff --git a/src/fix_wall_reflect.cpp b/src/fix_wall_reflect.cpp index 7e93c2da0e..26d35980e6 100644 --- a/src/fix_wall_reflect.cpp +++ b/src/fix_wall_reflect.cpp @@ -74,7 +74,7 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) : if (iarg+2 > narg) error->all("Illegal fix wall/reflect command"); if (strcmp(arg[iarg+1],"EDGE") == 0) { yloflag = EDGE; - ylo = domain->boxlo[0]; + ylo = domain->boxlo[1]; } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { yloflag = VARIABLE; int n = strlen(&arg[iarg+1][2]) + 1; @@ -89,7 +89,7 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) : if (iarg+2 > narg) error->all("Illegal fix wall/reflect command"); if (strcmp(arg[iarg+1],"EDGE") == 0) { yhiflag = EDGE; - yhi = domain->boxhi[0]; + yhi = domain->boxhi[1]; } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { yhiflag = VARIABLE; int n = strlen(&arg[iarg+1][2]) + 1; @@ -104,7 +104,7 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) : if (iarg+2 > narg) error->all("Illegal fix wall/reflect command"); if (strcmp(arg[iarg+1],"EDGE") == 0) { zloflag = EDGE; - zlo = domain->boxlo[0]; + zlo = domain->boxlo[2]; } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { zloflag = VARIABLE; int n = strlen(&arg[iarg+1][2]) + 1; @@ -119,7 +119,7 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) : if (iarg+2 > narg) error->all("Illegal fix wall/reflect command"); if (strcmp(arg[iarg+1],"EDGE") == 0) { zhiflag = EDGE; - zhi = domain->boxhi[0]; + zhi = domain->boxhi[2]; } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { zhiflag = VARIABLE; int n = strlen(&arg[iarg+1][2]) + 1; diff --git a/src/pair.cpp b/src/pair.cpp index d72566b4e2..e108c712a1 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -23,7 +23,6 @@ #include "stdlib.h" #include "string.h" #include "pair.h" -#include "pair_soft.h" #include "atom.h" #include "neighbor.h" #include "neigh_list.h" @@ -895,13 +894,6 @@ void Pair::write_file(int narg, char **arg) force->init(); - // if pair style = soft, set prefactor to final value - - Pair *spair = force->pair_match("soft",1); - if (spair) - ((PairSoft *) spair)->prefactor[itype][jtype] = - ((PairSoft *) spair)->prestop[itype][jtype]; - // if pair style = any of EAM, swap in dummy fp vector double eamfp[2]; diff --git a/src/pair.h b/src/pair.h index 8bb36b13c8..3b3e3f9f36 100644 --- a/src/pair.h +++ b/src/pair.h @@ -97,6 +97,8 @@ class Pair : protected Pointers { virtual void swap_eam(double *, double **) {} virtual void reset_dt() {} virtual void min_pointers(double **, double **) {} + virtual int pre_adapt(char *, int, int, int, int) {return -1;} + virtual void adapt(int, int, int, int, int, double) {} protected: int allocated; // 0/1 = whether arrays are allocated diff --git a/src/pair_hybrid.h b/src/pair_hybrid.h index 6b66ccd5ec..9554d0595d 100644 --- a/src/pair_hybrid.h +++ b/src/pair_hybrid.h @@ -32,7 +32,7 @@ class PairHybrid : public Pair { char **keywords; // style name of each Pair style PairHybrid(class LAMMPS *); - ~PairHybrid(); + virtual ~PairHybrid(); void compute(int, int); void settings(int, char **); virtual void coeff(int, char **); diff --git a/src/pair_hybrid_overlay.h b/src/pair_hybrid_overlay.h index 770e0e066c..dea84bae4c 100644 --- a/src/pair_hybrid_overlay.h +++ b/src/pair_hybrid_overlay.h @@ -27,6 +27,7 @@ namespace LAMMPS_NS { class PairHybridOverlay : public PairHybrid { public: PairHybridOverlay(class LAMMPS *); + ~PairHybridOverlay() {} void coeff(int, char **); private: diff --git a/src/pair_soft.cpp b/src/pair_soft.cpp index 92fab52cf6..8add30e678 100644 --- a/src/pair_soft.cpp +++ b/src/pair_soft.cpp @@ -14,6 +14,7 @@ #include "math.h" #include "stdio.h" #include "stdlib.h" +#include "string.h" #include "pair_soft.h" #include "atom.h" #include "comm.h" @@ -43,8 +44,6 @@ PairSoft::~PairSoft() memory->destroy_2d_int_array(setflag); memory->destroy_2d_double_array(cutsq); - memory->destroy_2d_double_array(prestart); - memory->destroy_2d_double_array(prestop); memory->destroy_2d_double_array(prefactor); memory->destroy_2d_double_array(cut); } @@ -63,21 +62,6 @@ void PairSoft::compute(int eflag, int vflag) if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; - // set current prefactor - // for minimization, set to prestop - // for dynamics, ramp from prestart to prestop - // for 0-step dynamics, set to prestart - - double delta = update->ntimestep - update->beginstep; - if (update->whichflag == 2) delta = 1.0; - else if (update->nsteps) delta /= update->endstep - update->beginstep; - else delta = 0.0; - int ntypes = atom->ntypes; - for (i = 1; i <= ntypes; i++) - for (j = 1; j <= ntypes; j++) - prefactor[i][j] = prestart[i][j] + - delta * (prestop[i][j] - prestart[i][j]); - double **x = atom->x; double **f = atom->f; int *type = atom->type; @@ -161,20 +145,8 @@ void PairSoft::allocate() cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); - prestart = memory->create_2d_double_array(n+1,n+1,"pair:prestart"); - prestop = memory->create_2d_double_array(n+1,n+1,"pair:prestop"); prefactor = memory->create_2d_double_array(n+1,n+1,"pair:prefactor"); cut = memory->create_2d_double_array(n+1,n+1,"pair:cut"); - - // init prestart and prestop to 0.0 - // since pair_hybrid can use all types even if pair_soft sub-class - // never sets them - - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) { - prestart[i][j] = 0.0; - prestop[i][j] = 0.0; - } } /* ---------------------------------------------------------------------- @@ -203,24 +175,22 @@ void PairSoft::settings(int narg, char **arg) void PairSoft::coeff(int narg, char **arg) { - if (narg < 4 || narg > 5) error->all("Incorrect args for pair coefficients"); + if (narg < 3 || narg > 4) error->all("Incorrect args for pair coefficients"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; force->bounds(arg[0],atom->ntypes,ilo,ihi); force->bounds(arg[1],atom->ntypes,jlo,jhi); - double prestart_one = force->numeric(arg[2]); - double prestop_one = force->numeric(arg[3]); + double prefactor_one = force->numeric(arg[2]); double cut_one = cut_global; - if (narg == 5) cut_one = force->numeric(arg[4]); + if (narg == 4) cut_one = force->numeric(arg[3]); int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { - prestart[i][j] = prestart_one; - prestop[i][j] = prestop_one; + prefactor[i][j] = prefactor_one; cut[i][j] = cut_one; setflag[i][j] = 1; count++; @@ -239,13 +209,11 @@ double PairSoft::init_one(int i, int j) // always mix prefactors geometrically if (setflag[i][j] == 0) { - prestart[i][j] = sqrt(prestart[i][i]*prestart[j][j]); - prestop[i][j] = sqrt(prestop[i][i]*prestop[j][j]); + prefactor[i][j] = sqrt(prefactor[i][i]*prefactor[j][j]); cut[i][j] = mix_distance(cut[i][i],cut[j][j]); } - prestart[j][i] = prestart[i][j]; - prestop[j][i] = prestop[i][j]; + prefactor[j][i] = prefactor[i][j]; cut[j][i] = cut[i][j]; return cut[i][j]; @@ -264,8 +232,7 @@ void PairSoft::write_restart(FILE *fp) for (j = i; j <= atom->ntypes; j++) { fwrite(&setflag[i][j],sizeof(int),1,fp); if (setflag[i][j]) { - fwrite(&prestart[i][j],sizeof(double),1,fp); - fwrite(&prestop[i][j],sizeof(double),1,fp); + fwrite(&prefactor[i][j],sizeof(double),1,fp); fwrite(&cut[i][j],sizeof(double),1,fp); } } @@ -289,12 +256,10 @@ void PairSoft::read_restart(FILE *fp) MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); if (setflag[i][j]) { if (me == 0) { - fread(&prestart[i][j],sizeof(double),1,fp); - fread(&prestop[i][j],sizeof(double),1,fp); + fread(&prefactor[i][j],sizeof(double),1,fp); fread(&cut[i][j],sizeof(double),1,fp); } - MPI_Bcast(&prestart[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&prestop[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&prefactor[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); } } @@ -324,6 +289,38 @@ void PairSoft::read_restart_settings(FILE *fp) MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } +/* ---------------------------------------------------------------------- + check if name is recognized, return integer index for that name + if name not recognized, return -1 + if type pair setting, return -2 if no type pairs are set +------------------------------------------------------------------------- */ + +int PairSoft::pre_adapt(char *name, int ilo, int ihi, int jlo, int jhi) +{ + int count = 0; + for (int i = ilo; i <= ihi; i++) + for (int j = MAX(jlo,i); j <= jhi; j++) + count++; + if (count == 0) return -2; + + if (strcmp(name,"a") == 0) return 0; + return -1; +} + +/* ---------------------------------------------------------------------- + adapt parameter indexed by which + change all pair variables affected by the reset parameter + if type pair setting, set I-J and J-I coeffs +------------------------------------------------------------------------- */ + +void PairSoft::adapt(int which, int ilo, int ihi, int jlo, int jhi, + double value) +{ + for (int i = ilo; i <= ihi; i++) + for (int j = MAX(jlo,i); j <= jhi; j++) + prefactor[i][j] = prefactor[j][i] = value; +} + /* ---------------------------------------------------------------------- */ double PairSoft::single(int i, int j, int itype, int jtype, double rsq, diff --git a/src/pair_soft.h b/src/pair_soft.h index 026c351f64..eceb98ec6f 100644 --- a/src/pair_soft.h +++ b/src/pair_soft.h @@ -38,12 +38,13 @@ class PairSoft : public Pair { void read_restart(FILE *); void write_restart_settings(FILE *); void read_restart_settings(FILE *); + int pre_adapt(char *, int, int, int, int); + void adapt(int, int, int, int, int, double); double single(int, int, int, int, double, double, double, double &); private: double PI; double cut_global; - double **prestart,**prestop; double **prefactor; double **cut; diff --git a/src/pair_table.cpp b/src/pair_table.cpp index d3966d6c75..a116857108 100644 --- a/src/pair_table.cpp +++ b/src/pair_table.cpp @@ -287,15 +287,13 @@ void PairTable::coeff(int narg, char **arg) // match = 1 if don't need to spline read-in tables // this is only the case if r values needed by final tables // exactly match r values read from file + // for tabstyle SPLINE, always need to build spline tables tb->match = 0; if (tabstyle == LINEAR && tb->ninput == tablength && tb->rflag == RSQ && tb->rhi == tb->cut) tb->match = 1; - if (tabstyle == SPLINE && tb->ninput == tablength && - tb->rflag == RSQ && tb->rhi == tb->cut) tb->match = 1; if (tabstyle == BITMAP && tb->ninput == 1 << tablength && tb->rflag == BMP && tb->rhi == tb->cut) tb->match = 1; - if (tb->rflag == BMP && tb->match == 0) error->all("Bitmapped table in file does not match requested table"); diff --git a/src/region_cone.cpp b/src/region_cone.cpp index b8a82914fd..79c003b8e2 100644 --- a/src/region_cone.cpp +++ b/src/region_cone.cpp @@ -436,7 +436,7 @@ int RegCone::surface_exterior(double *x, double cutoff) // y is exterior to cone or on its surface // corner1234 = 4 corner pts of half trapezoid = cone surf in plane of y - // project y to 3 line segments in half trapezoid (4th is axis of cone) + // project x to 3 line segments in half trapezoid (4th is axis of cone) // nearest = point on surface of cone that y is closest to // could be edge of cone // do not add contact point if r >= cutoff @@ -481,7 +481,7 @@ int RegCone::surface_exterior(double *x, double cutoff) // z is exterior to cone or on its surface // corner1234 = 4 corner pts of half trapezoid = cone surf in plane of z - // project z to 3 line segments in half trapezoid (4th is axis of cone) + // project x to 3 line segments in half trapezoid (4th is axis of cone) // nearest = point on surface of cone that z is closest to // could be edge of cone // do not add contact point if r >= cutoff diff --git a/src/region_cylinder.cpp b/src/region_cylinder.cpp index 541b71d5e7..58bb51f212 100644 --- a/src/region_cylinder.cpp +++ b/src/region_cylinder.cpp @@ -352,7 +352,7 @@ int RegCylinder::surface_exterior(double *x, double cutoff) if (r < radius && x[1] > lo && x[1] < hi) return 0; // y is exterior to cylinder or on its surface - // xp,yp,zp = point on surface of cylinder that y is closest to + // xp,yp,zp = point on surface of cylinder that x is closest to // could be edge of cylinder // do not add contact point if r >= cutoff @@ -383,7 +383,7 @@ int RegCylinder::surface_exterior(double *x, double cutoff) if (r < radius && x[2] > lo && x[2] < hi) return 0; // z is exterior to cylinder or on its surface - // xp,yp,zp = point on surface of cylinder that z is closest to + // xp,yp,zp = point on surface of cylinder that x is closest to // could be edge of cylinder // do not add contact point if r >= cutoff diff --git a/src/version.h b/src/version.h index 9aa8998568..62e2d576b7 100644 --- a/src/version.h +++ b/src/version.h @@ -1 +1 @@ -#define LAMMPS_VERSION "19 Jun 2010" +#define LAMMPS_VERSION "24 Jun 2010"