diff -Naur lammps-16Jan08/doc/Section_errors.html lammps-17Jan08/doc/Section_errors.html --- lammps-16Jan08/doc/Section_errors.html 2007-11-05 16:34:17.000000000 -0700 +++ lammps-17Jan08/doc/Section_errors.html 2008-01-17 10:02:19.000000000 -0700 @@ -222,6 +222,10 @@
No angle coefficients have been assigned in the data file or via the angle_coeff command. +
Angle potential must be defined for SHAKE + +
When shaking angles, an angle_style potential must be used. +
Angle style hybrid cannot have hybrid as an argument
Self-explanatory. @@ -294,6 +298,16 @@
Self-explanatory. +
Atom vector in equal-style variable formula + +
Atom vectors generate one value per atom which is not allowed +in an equal-style variable. + +
Atom-style variable in equal-style variable formula + +
Atom-style variables generate one value per atom which is not allowed +in an equal-style variable. +
Atom_modify command after simulation box is defined
The atom_modify command cannot be used after a read_data, @@ -380,11 +394,6 @@
Cannot use fix shake unless bond potential is defined. -
Bond style does not support computing per-atom bond energy - -
The bond style does not have a single() function, so it can not be -used to compute per-atom bond energy. -
Bond style hybrid cannot have hybrid as an argument
Self-explanatory. @@ -449,10 +458,6 @@
This feature is not yet supported. -
Cannot build parse tree for variable that is not atom style - -
Only atom style variables can be evaluated once per atom. -
Cannot change box to orthogonal when tilt is non-zero
Self-explanatory @@ -470,23 +475,16 @@
The frequency of writing dump dcd snapshots cannot be changed. +
Cannot change timestep with fix pour + +
This fix pre-computes some values based on the timestep, so it cannot +be changed during a simulation run. +
Cannot compute PPPM G
LAMMPS failed to compute a valid approximation for the PPPM g_ewald factor that partitions the computation between real space and k-space. -
Cannot compute PPPM X grid spacing - -
LAMMPS failed to compute a valid PPPM grid spacing in the x dimension. - -
Cannot compute PPPM Y grid spacing - -
LAMMPS failed to compute a valid PPPM grid spacing in the y dimension. - -
Cannot compute PPPM Z grid spacing - -
LAMMPS failed to compute a valid PPPM grid spacing in the z dimension. -
Cannot create an atom map unless atoms have IDs
The simulation requires a mapping from global atom IDs to local atoms, @@ -510,10 +508,6 @@
Use dump custom with x,y,z instead. -
Cannot evaluate variable - -
Self-explanatory. -
Cannot find delete_bonds group ID
Group ID used in the delete_bonds command does not exist. @@ -578,11 +572,6 @@
The output file for the fix com command cannot be opened. Check that the path and name are correct. -
Cannot open fix gran/diag file %s - -
The output file for the fix gran/diag command cannot be opened. Check -that the path and name are correct. -
Cannot open fix gyration file %s
Self-explanatory. @@ -597,6 +586,10 @@
The specified file cannot be opened. Check that the path and name are correct. +
Cannot open fix print file %s + +
The output file generated by the fix print command cannot be opened +
Cannot open fix rdf file %s
The output file for the fix rdf command cannot be opened. Check that @@ -690,6 +683,17 @@ You should use the replicate command before fixes are applied to the system. +
Cannot reset timestep with dump file already written to + +
Changing the timestep will confuse when a dump file is written. Use +the undump command, then restart the dump file. + +
Cannot reset timestep with restart file already written + +
Changing the timestep will confuse when a restart file is written. +Use the "restart 0" command to turn off restarts, then start them +again. +
Cannot run 2d simulation with nonperiodic Z dimension
Use the boundary command to make the z dimension periodic in order to @@ -745,20 +749,6 @@
The kspace style pppm cannot be used in 2d simulations. You can use 2d PPPM in a 3d simulation; see the kspace_modify command. -
Cannot use atom vector in variable unless atom map exists - -
Atom vectors require the ability to lookup an atom index, which -is provided by an atom map. An atom map does not exist (by default) -for non-molecular problems. Using the atom_modify map command will -force an atom map to be created. - -
Cannot use fix tmd unless atom map exists - -
TMD requires the ability to lookup an atom index, which -is provided by an atom map. An atom map does not exist (by default) -for non-molecular problems. Using the atom_modify map command will -force an atom map to be created. -
Cannot use both region, partial options in fix temp/rescale
Self-explanatory. @@ -776,6 +766,13 @@
Your choice of atom style does not have bonds. +
Cannot use fix TMD unless atom map exists + +
Using this fix requires the ability to lookup an atom index, which is +provided by an atom map. An atom map does not exist (by default) for +non-molecular problems. Using the atom_modify map command will force +an atom map to be created. +
Cannot use fix deform trate on a box with zero tilt
The trate style alters the current strain. @@ -797,6 +794,11 @@
The defined atom style uses per-atom mass, not per-type mass. +
Cannot use fix npt and fix deform on same dimension + +
These commands both change the box size/shape, so you cannot use both +together. +
Cannot use fix npt on a non-periodic dimension
Only periodic dimensions can be controlled with this fix. @@ -884,6 +886,10 @@
Cannot change orthogonal box to orthogonal or a triclinic box to triclinic. +
Compute ID for compute sum does not exist + +
Self-explanatory. +
Compute ID for fix ave/atom does not exist
Self-explanatory. @@ -896,14 +902,29 @@
Self-explanatory. +
Compute ID must be alphanumeric or underscore characters + +
Self-explanatory. +
Compute coord/atom cutoff is longer than pairwise cutoff
Cannot compute coordination at distances longer than the pair cutoff, since those atoms are not in the neighbor list. +
Compute in variable formula before initial run + +
Calculating this compute before the first run is not allowed because +various quantities may not yet be initialized. + +
Compute pe must use group all + +
Energies computed by potentials (pair, bond, etc) are computed on all +atoms. +
Compute pressure must use group all -
Self-explanatory. +
Virial contributions computed by potentials (pair, bond, etc) are +computed on all atoms.
Compute pressure temp ID does not compute temperature @@ -914,9 +935,33 @@
The atom style defined does not have these attributes. -
Compute sum/atom compute does not compute vector per atom +
Compute sum compute does not calculate a per-atom scalar + +
A compute accessed by compute sum must produce per-atom values. + +
Compute sum compute does not calculate a per-atom vector + +
A compute accessed by compute sum must produce per-atom values. + +
Compute sum compute does not calculate per-atom values + +
A compute accessed by compute sum must produce per-atom values. + +
Compute sum fix does not calculate a per-atom scalar -
If one compute in list calculates a vector, all of them must. +
A fix accessed by compute sum must produce per-atom values. + +
Compute sum fix does not calculate a per-atom vector + +
A fix accessed by compute sum must produce per-atom values. + +
Compute sum fix does not calculate per-atom values + +
A fix accessed by compute sum must produce per-atom values. + +
Compute sum variable is not atom-style variable + +
A variable accessed by compute sum must produce per-atom values.
Compute temp/asphere requires atom attributes quat, angmom @@ -926,6 +971,10 @@
The atom style defined does not have these attributes. +
Compute vector in variable formula is too small + +
The index to the vector is out of bounds. +
Could not create 3d FFT plan
The FFT setup in pppm failed. @@ -938,6 +987,10 @@
Self-explanatory. +
Could not find compute displace/atom fix ID + +
Self-explanatory. +
Could not find compute group ID
Self-explanatory. @@ -946,19 +999,6 @@
The compute ID for calculating temperature does not exist. -
Could not find compute sum/atom pre-compute ID - -
The compute ID for calculating the requested per-atom quantity does -not exist. - -
Could not find compute variable name - -
The variable being used by a compute is not defined. - -
Could not find compute variable name - -
The variable name accessed by compute variable/atom does not exist. -
Could not find compute_modify ID
Self-explanatory. @@ -988,6 +1028,10 @@
Self-explanatory. +
Could not find dump custom variable name + +
Self-explanatory. +
Could not find dump group ID
A group ID used in the dump command does not exist. @@ -1046,10 +1090,9 @@
The fix ID needed by thermo style custom to compute a requested quantity does not exist. -
Could not find thermo custom variable ID +
Could not find thermo custom variable name -
The variable name needed by thermo style custom to compute a requested -quantity does not exist. +
Self-explanatory.
Could not find thermo fix ID @@ -1088,11 +1131,6 @@
The compute ID needed by the velocity command to compute temperature does not exist. -
Could not pre-compute in variable - -
A compute required to evaulate a variable does not have its pre-compute -defined. -
Coulomb cutoffs of pair hybrid sub-styles do not match
If using a Kspace solver, all Coulomb cutoffs of long pair styles must @@ -1267,6 +1305,10 @@
Cannot use tilt factors unless the simulation box is non-orthogonal. +
Divide by 0 in variable formula + +
Self-explanatory. +
Domain too large for neighbor bins
The domain has become extremely large so that neighbor bins cannot be @@ -1313,6 +1355,11 @@
The dump custom command is attempting to access an out-of-range vector value. +
Dump custom variable is not atom-style variable + +
Only atom-style variables generate per-atom quantities, needed for +dump output. +
Dump dcd must use group all
Self-explanatory. @@ -1338,6 +1385,10 @@
The chosen atom style does not define the per-atom vector being dumped. +
Empty brackets in variable formula + +
A value inside the brackets is required for this formula element. +
Failed to allocate %d bytes for array %s
Your LAMMPS simulation has run out of memory. You need to run a @@ -1352,46 +1403,107 @@
Self-explanatory. +
Fix ID for compute sum does not exist + +
Self-explanatory. + +
Fix ID for fix ave/atom does not exist + +
Self-explanatory. +
Fix ID for fix ave/spatial does not exist
Self-explanatory. -
Fix ave/atom cannot be started on this timestep +
Fix ID for fix ave/time does not exist -
The first timestep on which data needs to be accumulated for -the average is before the current timestep. +
Self-explanatory. -
Fix ave/atom compute does not calculate per-atom info +
Fix ID must be alphanumeric or underscore characters -
The specified compute does not calculate per-atom values. +
Self-explanatory. -
Fix ave/spatial and fix not computed at compatible times +
Fix ave/atom compute does not calculate a per-atom scalar -
The fix must produce per-atom quantities on timesteps that fix -ave/spatial needs them. +
A compute used by fix ave/atom must generate per-atom values. -
Fix ave/spatial cannot be started on this timestep +
Fix ave/atom compute does not calculate a per-atom vector -
The first timestep on which data needs to be accumulated for -the average is before the current timestep. +
A compute used by fix ave/atom must generate per-atom values. -
Fix ave/spatial compute does not calculate per-atom info +
Fix ave/atom compute does not calculate per-atom values -
Only computes that calculate a per-atom quantity (not a scalar or -vector quantity can be used with fix ave/spatial. +
A compute used by fix ave/atom must generate per-atom values. -
Fix ave/spatial fix does not calculate per-atom info +
Fix ave/atom compute vector is accessed out-of-range -
The specified fix does not calculate per-atom values. +
The index for the vector is out of bounds. + +
Fix ave/atom fix does not calculate a per-atom scalar + +
A fix used by fix ave/atom must generate per-atom values. + +
Fix ave/atom fix does not calculate a per-atom vector + +
A fix used by fix ave/atom must generate per-atom values. + +
Fix ave/atom fix does not calculate per-atom values + +
A fix used by fix ave/atom must generate per-atom values. + +
Fix ave/atom fix vector is accessed out-of-range + +
The index for the vector is out of bounds. + +
Fix ave/atom variable is not atom-style variable + +
A variable used by fix ave/atom must generate per-atom values. + +
Fix ave/spatial compute does not calculate a per-atom scalar + +
A compute used by fix ave/spatial must generate per-atom values. + +
Fix ave/spatial compute does not calculate a per-atom vector + +
A compute used by fix ave/spatial must generate per-atom values. + +
Fix ave/spatial compute does not calculate per-atom values + +
A compute used by fix ave/spatial must generate per-atom values. + +
Fix ave/spatial compute vector is accessed out-of-range + +
The index for the vector is out of bounds. + +
Fix ave/spatial fix does not calculate a per-atom scalar + +
A fix used by fix ave/spatial must generate per-atom values. + +
Fix ave/spatial fix does not calculate a per-atom vector + +
A fix used by fix ave/spatial must generate per-atom values. + +
Fix ave/spatial fix does not calculate per-atom values + +
A fix used by fix ave/spatial must generate per-atom values. + +
Fix ave/spatial fix vector is accessed out-of-range + +
The index for the vector is out of bounds.
Fix ave/spatial for triclinic boxes requires units reduced
Self-explanatory. -
Fix ave/time cannot be started on this timestep +
Fix ave/spatial settings invalid with changing box + +
If the ave setting is "running" or "window" and the box size/shape +changes during the simulation, then the units setting must be +"reduced", else the number of bins may change. + +
Fix ave/spatial variable is not atom-style variable -
The first timestep on which data needs to be accumulated for -the average is before the current timestep. +
A variable used by fix ave/spatial must generate per-atom values.
Fix ave/time compute does not calculate a scalar @@ -1403,6 +1515,26 @@
Only computes that calculate a scalar or vector quantity (not a per-atom quantity) can be used with fix ave/time. +
Fix ave/time compute vector is accessed out-of-range + +
The index for the vector is out of bounds. + +
Fix ave/time fix does not calculate a scalar + +
A fix used by fix ave/time must generate global values. + +
Fix ave/time fix does not calculate a vector + +
A fix used by fix ave/time must generate global values. + +
Fix ave/time fix vector is accessed out-of-range + +
The index for the vector is out of bounds. + +
Fix ave/time variable is not equal-style variable + +
A variable used by fix ave/time must generate a global value. +
Fix command before simulation box is defined
The fix command cannot be used before a read_data, read_restart, or @@ -1427,22 +1559,39 @@
Self-explanatory -
Fix freeze requires atom attribute torque +
Fix for fix ave/atom not computed at compatible time -
The atom style defined does not have this attribute. +
Fixes generate their values on specific timesteps. Fix ave/atom is +requesting a value on a non-allowed timestep. -
Fix gran/diag is incompatible with Pair style +
Fix for fix ave/spatial not computed at compatible time -
Must use atom style granular. +
Fixes generate their values on specific timesteps. Fix ave/spatial is +requesting a value on a non-allowed timestep. -
Fix gran/diag requires atom attributes radius, rmass, omega +
Fix for fix ave/time not computed at compatible time -
The atom style defined does not have these attributes. +
Fixes generate their values on specific timesteps. Fix ave/time +is requesting a value on a non-allowed timestep. + +
Fix freeze requires atom attribute torque + +
The atom style defined does not have this attribute.
Fix heat group has no atoms
Self-explanatory. +
Fix in variable formula before initial run + +
Calculating this fix before the first run is not allowed because +various quantities may not yet be initialized. + +
Fix in variable not computed at compatible time + +
Fixes generate their values on specific timesteps. The variable is +requesting the values on a non-allowed timestep. +
Fix langevin period must be > 0.0
The time window for temperature relaxation must be > 0 @@ -1518,6 +1667,15 @@ integration fixes (nve, nvt, npt). See the fix tmd documentation for details. +
Fix used in compute sum not computed at compatible time + +
Fixes generate their values on specific timesteps. Compute sum is +requesting a value on a non-allowed timestep. + +
Fix vector in variable formula is too small + +
Index into vector is out of bounds. +
Fix wall/gran is incompatible with Pair style
Must use a granular pair style to define the parameters needed for @@ -1555,6 +1713,10 @@
A group ID used in the group command does not exist. +
Group ID in variable formula does not exist + +
Self-explanatory. +
Group command before simulation box is defined
The group command cannot be used before a read_data, read_restart, or @@ -1566,8 +1728,7 @@
Illegal ... command -
** DELETE_POSSIBLE -Self-explanatory. Check the input script syntax and compare to the +
Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. @@ -1655,10 +1816,6 @@
The data file header lists improper but no improper types. -
Inconsistent sizes of compute sum/atom compute vectors - -
All vectors of the computes in the list must be the same length. -
Incorrect args for angle coefficients
Self-explanatory. Check the input script or data file. @@ -1735,6 +1892,13 @@
Self-explanatory. Check the input script or data file. +
Indexed per-atom vector in variable formula without atom map + +
Accessing a value from an atom vector requires the ability to lookup +an atom index, which is provided by an atom map. An atom map does not +exist (by default) for non-molecular problems. Using the atom_modify +map command will force an atom map to be created. +
Induced tilt by displace_box is too large
The final tilt value must be between -1/2 and 1/2 of the perpendicular @@ -1836,9 +2000,9 @@
Atom types must range from 1 to Ntypes inclusive. -
Invalid atom vector in variable +
Invalid atom vector in variable formula -
An atom vector specified in a variable definition is not recognized. +
The atom vector is not recognized.
Invalid bond style @@ -1870,14 +2034,9 @@
One or more command-line arguments is invalid. Check the syntax of the command you are using to launch LAMMPS. -
Invalid compute ID in variable - -
A compute specified in a variable definition is not defined. +
Invalid compute ID in variable formula -
Invalid compute ID index in variable - -
The argument index of a compute specified in a variable definition is -not valid. +
The compute is not recognized.
Invalid compute style @@ -1991,6 +2150,10 @@
Operator keyword used for threshhold specification in not recognized. +
Invalid fix ID in variable formula + +
The fix is not recognized. +
Invalid fix nph command for a 2d simulation
Cannot use style xy, yz, or xz for a 2d simulation. @@ -1999,6 +2162,10 @@
Cannot use style xy, yz, or xz for a 2d simulation. +
Invalid fix style used in compute displace/atom command + +
Only a fix of style coord/original can be used with this compute. +
Invalid fix style
The choice of fix style is unknown. @@ -2023,6 +2190,10 @@
A group ID used in the neigh_modify command does not exist. +
Invalid group function in variable formula + +
Group function is not recognized. +
Invalid improper style
The choice of improper style is unknown. @@ -2032,6 +2203,10 @@
Improper type must be positive integer and within range of specified improper types. +
Invalid index in variable formula + +
The index between brackets is less than or equal to 0. +
Invalid keyword in dump custom command
One or more attribute keywords are not recognized. @@ -2052,9 +2227,9 @@
Self-explanatory. -
Invalid math/group function in variable +
Invalid math or group function in variable formula -
Self-explanatory. +
The math or group function is not recognized.
Invalid natoms for dump dcd @@ -2101,9 +2276,15 @@
The choice of region style is unknown. +
Invalid seed for Marsaglia random # generator + +
The initial seed for this random number generator must be a positive +integer less than or equal to 900 million. +
Invalid seed for Park random # generator -
The random number generator cannot be given a seed <= 0. +
The initial seed for this random number generator must be a positive +integer.
Invalid shape line in data file @@ -2117,10 +2298,14 @@
Self-explanatory. Check the input script. -
Invalid thermo keyword in variable +
Invalid syntax in variable formula
Self-explanatory. +
Invalid thermo keyword in variable formula + +
The keyword is not recognized. +
Invalid type for dipole set
Dipole command must set a type from 1-N where N is the number of atom @@ -2140,13 +2325,17 @@
The value specified for the setting is invalid, likely because it is too small or too large. +
Invalid variable evaluation in variable formula + +
A variable used in a formula could not be evaluated. +
Invalid variable in next command
Next command in input script must set variables from "a" to "z". -
Invalid variable name in variable +
Invalid variable name in variable formula -
Self-explanatory. +
Variable name is not recognized.
Invalid variable name @@ -2213,6 +2402,10 @@
2d simulation can use sq, sq2, or hex lattice. 3d simulation can use sc, bcc, or fcc lattice. +
Log of zero/negative in variable formula + +
Self-explanatory. +
Lost atoms via displace_atoms: original %.15g current %.15g
The displace_atoms command lost one or more atoms. @@ -2229,11 +2422,6 @@
A call to the MEAM Fortran library returned an error. -
Marsaglia RNG cannot use 0 seed - -
The random number generator use for the fix langevin command cannot -use 0 as an initial seed. -
Mass command before simulation box is defined
The mass command cannot be used before a read_data, read_restart, or @@ -2244,11 +2432,32 @@
The min_style command cannot be used before a read_data, read_restart, or create_box command. +
Minimization could not find thermo_pe compute + +
This compute is created by the thermo command. It must have been +explicitly deleted by a uncompute command. +
Minimize command before simulation box is defined
The minimize command cannot be used before a read_data, read_restart, or create_box command. +
Mismatched compute in variable formula + +
A compute is referenced incorrectly or a compute that produces per-atom +values is used in an equal-style variable formula. + +
Mismatched fix in variable formula + +
A fix is referenced incorrectly or a fix that produces per-atom +values is used in an equal-style variable formula. + +
Mismatched variable in variable formula + +
A variable is referenced incorrectly or an atom-style variable that +produces per-atom values is used in an equal-style variable +formula. +
More than one fix deform
Only one fix deform can be defined at a time. @@ -2380,6 +2589,11 @@
Self-explanatory. +
Must use a bond style with TIP4P potential + +
TIP4P potentials assume bond lengths in water are constrained +by a fix shake command. +
Must use a molecular atom style with fix poems molecule
Self-explanatory. @@ -2393,6 +2607,11 @@
The axis of the cylinder region used with the fix insert command must be oriented along the z dimension. +
Must use an angle style with TIP4P potential + +
TIP4P potentials assume angles in water are constrained by a fix shake +command. +
Must use charged atom style with fix efield
The atom style being used does not allow atoms to have assigned @@ -2447,6 +2666,10 @@
Self-explanatory. +
Neighbor page size must be >= 10x the one atom setting + +
This is required to prevent wasting too much memory. +
Newton bond change after simulation box is defined
The newton command cannot be used to change the newton bond value @@ -2577,6 +2800,34 @@
The specified cutoffs for the pair style are inconsistent. +
Pair lubricate only available for 3d + +
Self-explanatory. + +
Pair lubricate requires atom attributes torque, shape + +
Use a different atom style. + +
Pair lubricate requires mono-disperse particles + +
This is a current restriction of this pair style. + +
Pair lubricate requires spherical particles + +
This is a current restriction of this pair style. + +
Pair resquared epsilon a,b,c coeffs are not all set + +
Self-explanatory. + +
Pair resquared epsilon and sigma coeffs are not all set + +
Self-explanatory. + +
Pair resquared requires atom attributes quat, torque + +
Use a different atom style. +
Pair style AIREBO requires atom IDs
Self-explanatory. @@ -2625,16 +2876,6 @@
The pair style does not have a single() function, so it can not be invoked by bond_style quartic. -
Pair style does not support computing per-atom energy - -
The pair style does not have a single() function, so it can -not be used to compute per-atom energy. - -
Pair style does not support computing per-atom stress - -
The pair style does not have a single() function, so it can -not be used to compute per-atom stress. -
Pair style does not support pair_write
The pair style does not have a single() function, so it can @@ -2738,6 +2979,14 @@
Self-explanatory. +
Per-atom compute in equal-style variable formula + +
Equal-style variables cannot use per-atom quantities. + +
Per-atom fix in equal-style variable formula + +
Equal-style variables cannot use per-atom quantities. +
Potential file has duplicate entry
The potential file for a SW or Tersoff potential has more than @@ -2753,20 +3002,9 @@
Granular potentials that include shear history effects can only be run with a newton setting where pairwise newton is "off". -
Precompute ID for fix ave/spatial does not exist - -
The compute used by fix ave/spatial requires a second pre-computation -compute, which isn't defined. - -
Precompute ID for fix ave/time does not exist - -
The compute used by fix ave/time requires a second pre-computation -compute, which isn't defined. - -
Precompute ID of compute for fix ave/atom does not exist +
Power by 0 in variable formula -
The compute specified for fix ave/atom has other pre-computed computes -defined. One of them does not exist. +
Self-explanatory.
Press ID for fix nph does not exist @@ -2803,11 +3041,6 @@
A numeric error occurred in the creation of a rigid body by the fix rigid command. -
Quotes in a single arg - -
A single word should not be quoted in the input script; only a set of -words with intervening spaces should be quoted. -
R0 < 0 for fix spring command
Equilibrium spring length is invalid. @@ -2958,10 +3191,14 @@
Self-explanatory. -
Substitution for undefined variable +
Sqrt of negative in variable formula -
The variable specified with a $ symbol in an input script command has -not been previously defined with a variable command. +
Self-explanatory. + +
Substitution for illegal variable + +
Input script line contained a variable that could not be substituted +for.
TIP4P hydrogen has incorrect atom type @@ -3006,16 +3243,6 @@
The compute ID needed to compute temperature for the fix does not exist. -
Temp ID of press ID for fix nph does not exist - -
The compute ID needed to compute temperature within the pressure -compute ID for the fix does not exist. - -
Temp ID of press ID for fix npt does not exist - -
The compute ID needed to compute temperature within the pressure -compute ID for the fix does not exist. -
Temper command before simulation box is defined
The temper command cannot be used before a read_data, read_restart, or @@ -3025,6 +3252,11 @@
The region ID specified in the temperature command does not exist. +
Tempering could not find thermo_pe compute + +
This compute is created by the thermo command. It must have been +explicitly deleted by a uncompute command. +
Tempering fix ID is not defined
The fix ID specified by the temper command does not exist. @@ -3034,6 +3266,11 @@
The fix specified by the temper command is not one that controls temperature (nvt or langevin). +
Thermo and fix not computed at compatible times + +
Fixes generate values on specific timesteps. The thermo output +does not match these timesteps. +
Thermo compute ID does not compute scalar info
The specified compute ID does not compute a scalar quantity @@ -3049,6 +3286,30 @@
The specified compute ID does not compute a large enough vector quantity for the requested index. +
Thermo custom variable is not equal-style variable + +
Only equal-style variables can be output with thermodynamics, not +atom-style variables. + +
Thermo fix ID does not compute scalar info + +
Only fixes that compute global values can be output with +thermodynamics. + +
Thermo fix ID does not compute vector info + +
Only fixes that compute global values can be output with +thermodynamics. + +
Thermo fix ID vector is not large enough + +
Index into vector is out of bounds. + +
Thermo keyword in variable formula before initial run + +
Calculating a thermo keyword before the first run is not allowed +because various quantities may not yet be initialized. +
Thermo style does not use drot
Cannot use thermo_modify to set this parameter since the thermo_style @@ -3082,16 +3343,6 @@
The thermo_style command cannot be used before a read_data, read_restart, or create_box command. -
Thermodynamics must compute PE for temper command - -
The thermo style must insure that thermodynamics computations include -potential energy when tempering is performed. - -
Thermodynamics not computed on tempering swap steps - -
The thermo command must insure that thermodynamics (including energy) -is computed on the timesteps that tempering swaps are attempted. -
Timestep must be >= 0
Specified timestep size is invalid. @@ -3154,6 +3405,10 @@ length in that dimension. E.g. the xy tilt must be between -half and +half of the x box length. +
Two groups cannot be the same in fix spring couple + +
Self-explanatory. +
Unbalanced quotes in input line
No matching end double quote was found following a leading double @@ -3229,6 +3484,11 @@
Must use lattice command with compute fix deposit command if units option is set to lattice. +
Use of fix dt/reset with undefined lattice + +
Must use lattice command with fix dt/reset command if units option is +set to lattice. +
Use of fix indent with undefined lattice
The lattice command must be used to define a lattice before using the @@ -3258,27 +3518,32 @@
Self-explanatory. -
Variable compute ID does not compute scalar info +
Variable name for compute sum does not exist -
The specified compute ID does not compute a scalar quantity -as requested. +
Self-explanatory. -
Variable compute ID vector is not large enough +
Variable name for fix ave/atom does not exist -
The specified compute ID does not compute a large enough vector -quantity for the requested index. +
Self-explanatory. -
Variable equal keyword used before initial run +
Variable name for fix ave/spatial does not exist -
Cannot evaluate the variable at this stage of input script. +
Self-explanatory. -
Variable equal keyword used before simulation box defined +
Variable name for fix ave/time does not exist -
Cannot evaluate the variable at this stage of input script. +
Self-explanatory. -
Variable group ID does not exist +
Variable name must be alphanumeric or underscore characters -
A group specified in a variable definition does not exist. +
Self-explanatory. + +
Variable uses compute via thermo keyword that thermo does not + +
A thermo keyword used in the variable requires a particular compute be +invoked. However, the thermodynamic output does not use this compute, +so it has not been initialized. Include the keyword in thermodynamic +output.
Velocity command before simulation box is defined @@ -3318,6 +3583,11 @@
+
Dump dcd/xtc timestamp may be wrong with fix dt/reset + +
If the fix changes the timestep, the dump dcd file will not +reflect the change. +
FENE bond too long: %d %d %d %g
A FENE bond has stretched dangerously far. It's interaction strength @@ -3368,22 +3638,10 @@
It is not efficient to use compute coord/atom more than once. -
More than one compute ebond/atom - -
It is not efficient to use compute ebond/atom more than once. - -
More than one compute epair/atom - -
It is not efficient to use compute epair/atom more than once. -
More than one compute ke/atom
It is not efficient to use compute ke/atom more than once. -
More than one compute stress/atom - -
It is not efficient to use compute stress/atom more than once. -
More than one fix msd
It is not efficient to use fix msd more than once. @@ -3566,10 +3824,6 @@ computed by integrating the density of a periodic system out to infinity. -
Variable equal keyword used with non-current thermo - -
The evaluation of the variable may be inaccurate as a result. -
diff -Naur lammps-16Jan08/doc/Section_errors.txt lammps-17Jan08/doc/Section_errors.txt --- lammps-16Jan08/doc/Section_errors.txt 2007-11-05 16:34:17.000000000 -0700 +++ lammps-17Jan08/doc/Section_errors.txt 2008-01-17 10:02:19.000000000 -0700 @@ -219,6 +219,10 @@ No angle coefficients have been assigned in the data file or via the angle_coeff command. :dd +{Angle potential must be defined for SHAKE} :dt + +When shaking angles, an angle_style potential must be used. :dd + {Angle style hybrid cannot have hybrid as an argument} :dt Self-explanatory. :dd @@ -291,6 +295,16 @@ Self-explanatory. :dd +{Atom vector in equal-style variable formula} :dt + +Atom vectors generate one value per atom which is not allowed +in an equal-style variable. :dd + +{Atom-style variable in equal-style variable formula} :dt + +Atom-style variables generate one value per atom which is not allowed +in an equal-style variable. :dd + {Atom_modify command after simulation box is defined} :dt The atom_modify command cannot be used after a read_data, @@ -377,11 +391,6 @@ Cannot use fix shake unless bond potential is defined. :dd -{Bond style does not support computing per-atom bond energy} :dt - -The bond style does not have a single() function, so it can not be -used to compute per-atom bond energy. :dd - {Bond style hybrid cannot have hybrid as an argument} :dt Self-explanatory. :dd @@ -446,10 +455,6 @@ This feature is not yet supported. :dd -{Cannot build parse tree for variable that is not atom style} :dt - -Only atom style variables can be evaluated once per atom. :dd - {Cannot change box to orthogonal when tilt is non-zero} :dt Self-explanatory :dd @@ -467,23 +472,16 @@ The frequency of writing dump dcd snapshots cannot be changed. :dd +{Cannot change timestep with fix pour} :dt + +This fix pre-computes some values based on the timestep, so it cannot +be changed during a simulation run. :dd + {Cannot compute PPPM G} :dt LAMMPS failed to compute a valid approximation for the PPPM g_ewald factor that partitions the computation between real space and k-space. :dd -{Cannot compute PPPM X grid spacing} :dt - -LAMMPS failed to compute a valid PPPM grid spacing in the x dimension. :dd - -{Cannot compute PPPM Y grid spacing} :dt - -LAMMPS failed to compute a valid PPPM grid spacing in the y dimension. :dd - -{Cannot compute PPPM Z grid spacing} :dt - -LAMMPS failed to compute a valid PPPM grid spacing in the z dimension. :dd - {Cannot create an atom map unless atoms have IDs} :dt The simulation requires a mapping from global atom IDs to local atoms, @@ -507,10 +505,6 @@ Use dump custom with x,y,z instead. :dd -{Cannot evaluate variable} :dt - -Self-explanatory. :dd - {Cannot find delete_bonds group ID} :dt Group ID used in the delete_bonds command does not exist. :dd @@ -575,11 +569,6 @@ The output file for the fix com command cannot be opened. Check that the path and name are correct. :dd -{Cannot open fix gran/diag file %s} :dt - -The output file for the fix gran/diag command cannot be opened. Check -that the path and name are correct. :dd - {Cannot open fix gyration file %s} :dt Self-explanatory. :dd @@ -594,6 +583,10 @@ The specified file cannot be opened. Check that the path and name are correct. :dd +{Cannot open fix print file %s} :dt + +The output file generated by the fix print command cannot be opened :dd + {Cannot open fix rdf file %s} :dt The output file for the fix rdf command cannot be opened. Check that @@ -687,6 +680,17 @@ You should use the replicate command before fixes are applied to the system. :dd +{Cannot reset timestep with dump file already written to} :dt + +Changing the timestep will confuse when a dump file is written. Use +the undump command, then restart the dump file. :dd + +{Cannot reset timestep with restart file already written} :dt + +Changing the timestep will confuse when a restart file is written. +Use the "restart 0" command to turn off restarts, then start them +again. :dd + {Cannot run 2d simulation with nonperiodic Z dimension} :dt Use the boundary command to make the z dimension periodic in order to @@ -742,20 +746,6 @@ The kspace style pppm cannot be used in 2d simulations. You can use 2d PPPM in a 3d simulation; see the kspace_modify command. :dd -{Cannot use atom vector in variable unless atom map exists} :dt - -Atom vectors require the ability to lookup an atom index, which -is provided by an atom map. An atom map does not exist (by default) -for non-molecular problems. Using the atom_modify map command will -force an atom map to be created. :dd - -{Cannot use fix tmd unless atom map exists} :dt - -TMD requires the ability to lookup an atom index, which -is provided by an atom map. An atom map does not exist (by default) -for non-molecular problems. Using the atom_modify map command will -force an atom map to be created. :dd - {Cannot use both region, partial options in fix temp/rescale} :dt Self-explanatory. :dd @@ -773,6 +763,13 @@ Your choice of atom style does not have bonds. :dd +{Cannot use fix TMD unless atom map exists} :dt + +Using this fix requires the ability to lookup an atom index, which is +provided by an atom map. An atom map does not exist (by default) for +non-molecular problems. Using the atom_modify map command will force +an atom map to be created. :dd + {Cannot use fix deform trate on a box with zero tilt} :dt The trate style alters the current strain. :dd @@ -794,6 +791,11 @@ The defined atom style uses per-atom mass, not per-type mass. :dd +{Cannot use fix npt and fix deform on same dimension} :dt + +These commands both change the box size/shape, so you cannot use both +together. :dd + {Cannot use fix npt on a non-periodic dimension} :dt Only periodic dimensions can be controlled with this fix. :dd @@ -881,6 +883,10 @@ Cannot change orthogonal box to orthogonal or a triclinic box to triclinic. :dd +{Compute ID for compute sum does not exist} :dt + +Self-explanatory. :dd + {Compute ID for fix ave/atom does not exist} :dt Self-explanatory. :dd @@ -893,14 +899,29 @@ Self-explanatory. :dd +{Compute ID must be alphanumeric or underscore characters} :dt + +Self-explanatory. :dd + {Compute coord/atom cutoff is longer than pairwise cutoff} :dt Cannot compute coordination at distances longer than the pair cutoff, since those atoms are not in the neighbor list. :dd +{Compute in variable formula before initial run} :dt + +Calculating this compute before the first run is not allowed because +various quantities may not yet be initialized. :dd + +{Compute pe must use group all} :dt + +Energies computed by potentials (pair, bond, etc) are computed on all +atoms. :dd + {Compute pressure must use group all} :dt -Self-explanatory. :dd +Virial contributions computed by potentials (pair, bond, etc) are +computed on all atoms. :dd {Compute pressure temp ID does not compute temperature} :dt @@ -911,9 +932,33 @@ The atom style defined does not have these attributes. :dd -{Compute sum/atom compute does not compute vector per atom} :dt +{Compute sum compute does not calculate a per-atom scalar} :dt + +A compute accessed by compute sum must produce per-atom values. :dd + +{Compute sum compute does not calculate a per-atom vector} :dt + +A compute accessed by compute sum must produce per-atom values. :dd + +{Compute sum compute does not calculate per-atom values} :dt + +A compute accessed by compute sum must produce per-atom values. :dd + +{Compute sum fix does not calculate a per-atom scalar} :dt -If one compute in list calculates a vector, all of them must. :dd +A fix accessed by compute sum must produce per-atom values. :dd + +{Compute sum fix does not calculate a per-atom vector} :dt + +A fix accessed by compute sum must produce per-atom values. :dd + +{Compute sum fix does not calculate per-atom values} :dt + +A fix accessed by compute sum must produce per-atom values. :dd + +{Compute sum variable is not atom-style variable} :dt + +A variable accessed by compute sum must produce per-atom values. :dd {Compute temp/asphere requires atom attributes quat, angmom} :dt @@ -923,6 +968,10 @@ The atom style defined does not have these attributes. :dd +{Compute vector in variable formula is too small} :dt + +The index to the vector is out of bounds. :dd + {Could not create 3d FFT plan} :dt The FFT setup in pppm failed. :dd @@ -935,6 +984,10 @@ Self-explanatory. :dd +{Could not find compute displace/atom fix ID} :dt + +Self-explanatory. :dd + {Could not find compute group ID} :dt Self-explanatory. :dd @@ -943,19 +996,6 @@ The compute ID for calculating temperature does not exist. :dd -{Could not find compute sum/atom pre-compute ID} :dt - -The compute ID for calculating the requested per-atom quantity does -not exist. :dd - -{Could not find compute variable name} :dt - -The variable being used by a compute is not defined. :dd - -{Could not find compute variable name} :dt - -The variable name accessed by compute variable/atom does not exist. :dd - {Could not find compute_modify ID} :dt Self-explanatory. :dd @@ -985,6 +1025,10 @@ Self-explanatory. :dd +{Could not find dump custom variable name} :dt + +Self-explanatory. :dd + {Could not find dump group ID} :dt A group ID used in the dump command does not exist. :dd @@ -1043,10 +1087,9 @@ The fix ID needed by thermo style custom to compute a requested quantity does not exist. :dd -{Could not find thermo custom variable ID} :dt +{Could not find thermo custom variable name} :dt -The variable name needed by thermo style custom to compute a requested -quantity does not exist. :dd +Self-explanatory. :dd {Could not find thermo fix ID} :dt @@ -1085,11 +1128,6 @@ The compute ID needed by the velocity command to compute temperature does not exist. :dd -{Could not pre-compute in variable} :dt - -A compute required to evaulate a variable does not have its pre-compute -defined. :dd - {Coulomb cutoffs of pair hybrid sub-styles do not match} :dt If using a Kspace solver, all Coulomb cutoffs of long pair styles must @@ -1264,6 +1302,10 @@ Cannot use tilt factors unless the simulation box is non-orthogonal. :dd +{Divide by 0 in variable formula} :dt + +Self-explanatory. :dd + {Domain too large for neighbor bins} :dt The domain has become extremely large so that neighbor bins cannot be @@ -1310,6 +1352,11 @@ The dump custom command is attempting to access an out-of-range vector value. :dd +{Dump custom variable is not atom-style variable} :dt + +Only atom-style variables generate per-atom quantities, needed for +dump output. :dd + {Dump dcd must use group all} :dt Self-explanatory. :dd @@ -1335,6 +1382,10 @@ The chosen atom style does not define the per-atom vector being dumped. :dd +{Empty brackets in variable formula} :dt + +A value inside the brackets is required for this formula element. :dd + {Failed to allocate %d bytes for array %s} :dt Your LAMMPS simulation has run out of memory. You need to run a @@ -1349,46 +1400,107 @@ Self-explanatory. :dd +{Fix ID for compute sum does not exist} :dt + +Self-explanatory. :dd + +{Fix ID for fix ave/atom does not exist} :dt + +Self-explanatory. :dd + {Fix ID for fix ave/spatial does not exist} :dt Self-explanatory. :dd -{Fix ave/atom cannot be started on this timestep} :dt +{Fix ID for fix ave/time does not exist} :dt -The first timestep on which data needs to be accumulated for -the average is before the current timestep. :dd +Self-explanatory. :dd -{Fix ave/atom compute does not calculate per-atom info} :dt +{Fix ID must be alphanumeric or underscore characters} :dt -The specified compute does not calculate per-atom values. :dd +Self-explanatory. :dd -{Fix ave/spatial and fix not computed at compatible times} :dt +{Fix ave/atom compute does not calculate a per-atom scalar} :dt -The fix must produce per-atom quantities on timesteps that fix -ave/spatial needs them. :dd +A compute used by fix ave/atom must generate per-atom values. :dd -{Fix ave/spatial cannot be started on this timestep} :dt +{Fix ave/atom compute does not calculate a per-atom vector} :dt -The first timestep on which data needs to be accumulated for -the average is before the current timestep. :dd +A compute used by fix ave/atom must generate per-atom values. :dd -{Fix ave/spatial compute does not calculate per-atom info} :dt +{Fix ave/atom compute does not calculate per-atom values} :dt -Only computes that calculate a per-atom quantity (not a scalar or -vector quantity can be used with fix ave/spatial. :dd +A compute used by fix ave/atom must generate per-atom values. :dd -{Fix ave/spatial fix does not calculate per-atom info} :dt +{Fix ave/atom compute vector is accessed out-of-range} :dt -The specified fix does not calculate per-atom values. :dd +The index for the vector is out of bounds. :dd + +{Fix ave/atom fix does not calculate a per-atom scalar} :dt + +A fix used by fix ave/atom must generate per-atom values. :dd + +{Fix ave/atom fix does not calculate a per-atom vector} :dt + +A fix used by fix ave/atom must generate per-atom values. :dd + +{Fix ave/atom fix does not calculate per-atom values} :dt + +A fix used by fix ave/atom must generate per-atom values. :dd + +{Fix ave/atom fix vector is accessed out-of-range} :dt + +The index for the vector is out of bounds. :dd + +{Fix ave/atom variable is not atom-style variable} :dt + +A variable used by fix ave/atom must generate per-atom values. :dd + +{Fix ave/spatial compute does not calculate a per-atom scalar} :dt + +A compute used by fix ave/spatial must generate per-atom values. :dd + +{Fix ave/spatial compute does not calculate a per-atom vector} :dt + +A compute used by fix ave/spatial must generate per-atom values. :dd + +{Fix ave/spatial compute does not calculate per-atom values} :dt + +A compute used by fix ave/spatial must generate per-atom values. :dd + +{Fix ave/spatial compute vector is accessed out-of-range} :dt + +The index for the vector is out of bounds. :dd + +{Fix ave/spatial fix does not calculate a per-atom scalar} :dt + +A fix used by fix ave/spatial must generate per-atom values. :dd + +{Fix ave/spatial fix does not calculate a per-atom vector} :dt + +A fix used by fix ave/spatial must generate per-atom values. :dd + +{Fix ave/spatial fix does not calculate per-atom values} :dt + +A fix used by fix ave/spatial must generate per-atom values. :dd + +{Fix ave/spatial fix vector is accessed out-of-range} :dt + +The index for the vector is out of bounds. :dd {Fix ave/spatial for triclinic boxes requires units reduced} :dt Self-explanatory. :dd -{Fix ave/time cannot be started on this timestep} :dt +{Fix ave/spatial settings invalid with changing box} :dt + +If the ave setting is "running" or "window" and the box size/shape +changes during the simulation, then the units setting must be +"reduced", else the number of bins may change. :dd + +{Fix ave/spatial variable is not atom-style variable} :dt -The first timestep on which data needs to be accumulated for -the average is before the current timestep. :dd +A variable used by fix ave/spatial must generate per-atom values. :dd {Fix ave/time compute does not calculate a scalar} :dt @@ -1400,6 +1512,26 @@ Only computes that calculate a scalar or vector quantity (not a per-atom quantity) can be used with fix ave/time. :dd +{Fix ave/time compute vector is accessed out-of-range} :dt + +The index for the vector is out of bounds. :dd + +{Fix ave/time fix does not calculate a scalar} :dt + +A fix used by fix ave/time must generate global values. :dd + +{Fix ave/time fix does not calculate a vector} :dt + +A fix used by fix ave/time must generate global values. :dd + +{Fix ave/time fix vector is accessed out-of-range} :dt + +The index for the vector is out of bounds. :dd + +{Fix ave/time variable is not equal-style variable} :dt + +A variable used by fix ave/time must generate a global value. :dd + {Fix command before simulation box is defined} :dt The fix command cannot be used before a read_data, read_restart, or @@ -1424,22 +1556,39 @@ Self-explanatory :dd -{Fix freeze requires atom attribute torque} :dt +{Fix for fix ave/atom not computed at compatible time} :dt -The atom style defined does not have this attribute. :dd +Fixes generate their values on specific timesteps. Fix ave/atom is +requesting a value on a non-allowed timestep. :dd -{Fix gran/diag is incompatible with Pair style} :dt +{Fix for fix ave/spatial not computed at compatible time} :dt -Must use atom style granular. :dd +Fixes generate their values on specific timesteps. Fix ave/spatial is +requesting a value on a non-allowed timestep. :dd -{Fix gran/diag requires atom attributes radius, rmass, omega} :dt +{Fix for fix ave/time not computed at compatible time} :dt -The atom style defined does not have these attributes. :dd +Fixes generate their values on specific timesteps. Fix ave/time +is requesting a value on a non-allowed timestep. :dd + +{Fix freeze requires atom attribute torque} :dt + +The atom style defined does not have this attribute. :dd {Fix heat group has no atoms} :dt Self-explanatory. :dd +{Fix in variable formula before initial run} :dt + +Calculating this fix before the first run is not allowed because +various quantities may not yet be initialized. :dd + +{Fix in variable not computed at compatible time} :dt + +Fixes generate their values on specific timesteps. The variable is +requesting the values on a non-allowed timestep. :dd + {Fix langevin period must be > 0.0} :dt The time window for temperature relaxation must be > 0 :dd @@ -1515,6 +1664,15 @@ integration fixes (nve, nvt, npt). See the fix tmd documentation for details. :dd +{Fix used in compute sum not computed at compatible time} :dt + +Fixes generate their values on specific timesteps. Compute sum is +requesting a value on a non-allowed timestep. :dd + +{Fix vector in variable formula is too small} :dt + +Index into vector is out of bounds. :dd + {Fix wall/gran is incompatible with Pair style} :dt Must use a granular pair style to define the parameters needed for @@ -1552,6 +1710,10 @@ A group ID used in the group command does not exist. :dd +{Group ID in variable formula does not exist} :dt + +Self-explanatory. :dd + {Group command before simulation box is defined} :dt The group command cannot be used before a read_data, read_restart, or @@ -1563,7 +1725,6 @@ {Illegal ... command} :dt -** DELETE_POSSIBLE Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. :dd @@ -1652,10 +1813,6 @@ The data file header lists improper but no improper types. :dd -{Inconsistent sizes of compute sum/atom compute vectors} :dt - -All vectors of the computes in the list must be the same length. :dd - {Incorrect args for angle coefficients} :dt Self-explanatory. Check the input script or data file. :dd @@ -1732,6 +1889,13 @@ Self-explanatory. Check the input script or data file. :dd +{Indexed per-atom vector in variable formula without atom map} :dt + +Accessing a value from an atom vector requires the ability to lookup +an atom index, which is provided by an atom map. An atom map does not +exist (by default) for non-molecular problems. Using the atom_modify +map command will force an atom map to be created. :dd + {Induced tilt by displace_box is too large} :dt The final tilt value must be between -1/2 and 1/2 of the perpendicular @@ -1833,9 +1997,9 @@ Atom types must range from 1 to Ntypes inclusive. :dd -{Invalid atom vector in variable} :dt +{Invalid atom vector in variable formula} :dt -An atom vector specified in a variable definition is not recognized. :dd +The atom vector is not recognized. :dd {Invalid bond style} :dt @@ -1867,14 +2031,9 @@ One or more command-line arguments is invalid. Check the syntax of the command you are using to launch LAMMPS. :dd -{Invalid compute ID in variable} :dt - -A compute specified in a variable definition is not defined. :dd +{Invalid compute ID in variable formula} :dt -{Invalid compute ID index in variable} :dt - -The argument index of a compute specified in a variable definition is -not valid. :dd +The compute is not recognized. :dd {Invalid compute style} :dt @@ -1988,6 +2147,10 @@ Operator keyword used for threshhold specification in not recognized. :dd +{Invalid fix ID in variable formula} :dt + +The fix is not recognized. :dd + {Invalid fix nph command for a 2d simulation} :dt Cannot use style xy, yz, or xz for a 2d simulation. :dd @@ -1996,6 +2159,10 @@ Cannot use style xy, yz, or xz for a 2d simulation. :dd +{Invalid fix style used in compute displace/atom command} :dt + +Only a fix of style coord/original can be used with this compute. :dd + {Invalid fix style} :dt The choice of fix style is unknown. :dd @@ -2020,6 +2187,10 @@ A group ID used in the neigh_modify command does not exist. :dd +{Invalid group function in variable formula} :dt + +Group function is not recognized. :dd + {Invalid improper style} :dt The choice of improper style is unknown. :dd @@ -2029,6 +2200,10 @@ Improper type must be positive integer and within range of specified improper types. :dd +{Invalid index in variable formula} :dt + +The index between brackets is less than or equal to 0. :dd + {Invalid keyword in dump custom command} :dt One or more attribute keywords are not recognized. :dd @@ -2049,9 +2224,9 @@ Self-explanatory. :dd -{Invalid math/group function in variable} :dt +{Invalid math or group function in variable formula} :dt -Self-explanatory. :dd +The math or group function is not recognized. :dd {Invalid natoms for dump dcd} :dt @@ -2098,9 +2273,15 @@ The choice of region style is unknown. :dd +{Invalid seed for Marsaglia random # generator} :dt + +The initial seed for this random number generator must be a positive +integer less than or equal to 900 million. :dd + {Invalid seed for Park random # generator} :dt -The random number generator cannot be given a seed <= 0. :dd +The initial seed for this random number generator must be a positive +integer. :dd {Invalid shape line in data file} :dt @@ -2114,10 +2295,14 @@ Self-explanatory. Check the input script. :dd -{Invalid thermo keyword in variable} :dt +{Invalid syntax in variable formula} :dt Self-explanatory. :dd +{Invalid thermo keyword in variable formula} :dt + +The keyword is not recognized. :dd + {Invalid type for dipole set} :dt Dipole command must set a type from 1-N where N is the number of atom @@ -2137,13 +2322,17 @@ The value specified for the setting is invalid, likely because it is too small or too large. :dd +{Invalid variable evaluation in variable formula} :dt + +A variable used in a formula could not be evaluated. :dd + {Invalid variable in next command} :dt Next command in input script must set variables from "a" to "z". :dd -{Invalid variable name in variable} :dt +{Invalid variable name in variable formula} :dt -Self-explanatory. :dd +Variable name is not recognized. :dd {Invalid variable name} :dt @@ -2210,6 +2399,10 @@ 2d simulation can use sq, sq2, or hex lattice. 3d simulation can use sc, bcc, or fcc lattice. :dd +{Log of zero/negative in variable formula} :dt + +Self-explanatory. :dd + {Lost atoms via displace_atoms: original %.15g current %.15g} :dt The displace_atoms command lost one or more atoms. :dd @@ -2226,11 +2419,6 @@ A call to the MEAM Fortran library returned an error. :dd -{Marsaglia RNG cannot use 0 seed} :dt - -The random number generator use for the fix langevin command cannot -use 0 as an initial seed. :dd - {Mass command before simulation box is defined} :dt The mass command cannot be used before a read_data, read_restart, or @@ -2241,11 +2429,32 @@ The min_style command cannot be used before a read_data, read_restart, or create_box command. :dd +{Minimization could not find thermo_pe compute} :dt + +This compute is created by the thermo command. It must have been +explicitly deleted by a uncompute command. :dd + {Minimize command before simulation box is defined} :dt The minimize command cannot be used before a read_data, read_restart, or create_box command. :dd +{Mismatched compute in variable formula} :dt + +A compute is referenced incorrectly or a compute that produces per-atom +values is used in an equal-style variable formula. :dd + +{Mismatched fix in variable formula} :dt + +A fix is referenced incorrectly or a fix that produces per-atom +values is used in an equal-style variable formula. :dd + +{Mismatched variable in variable formula} :dt + +A variable is referenced incorrectly or an atom-style variable that +produces per-atom values is used in an equal-style variable +formula. :dd + {More than one fix deform} :dt Only one fix deform can be defined at a time. :dd @@ -2377,6 +2586,11 @@ Self-explanatory. :dd +{Must use a bond style with TIP4P potential} :dt + +TIP4P potentials assume bond lengths in water are constrained +by a fix shake command. :dd + {Must use a molecular atom style with fix poems molecule} :dt Self-explanatory. :dd @@ -2390,6 +2604,11 @@ The axis of the cylinder region used with the fix insert command must be oriented along the z dimension. :dd +{Must use an angle style with TIP4P potential} :dt + +TIP4P potentials assume angles in water are constrained by a fix shake +command. :dd + {Must use charged atom style with fix efield} :dt The atom style being used does not allow atoms to have assigned @@ -2444,6 +2663,10 @@ Self-explanatory. :dd +{Neighbor page size must be >= 10x the one atom setting} :dt + +This is required to prevent wasting too much memory. :dd + {Newton bond change after simulation box is defined} :dt The newton command cannot be used to change the newton bond value @@ -2574,6 +2797,34 @@ The specified cutoffs for the pair style are inconsistent. :dd +{Pair lubricate only available for 3d} :dt + +Self-explanatory. :dd + +{Pair lubricate requires atom attributes torque, shape} :dt + +Use a different atom style. :dd + +{Pair lubricate requires mono-disperse particles} :dt + +This is a current restriction of this pair style. :dd + +{Pair lubricate requires spherical particles} :dt + +This is a current restriction of this pair style. :dd + +{Pair resquared epsilon a,b,c coeffs are not all set} :dt + +Self-explanatory. :dd + +{Pair resquared epsilon and sigma coeffs are not all set} :dt + +Self-explanatory. :dd + +{Pair resquared requires atom attributes quat, torque} :dt + +Use a different atom style. :dd + {Pair style AIREBO requires atom IDs} :dt Self-explanatory. :dd @@ -2622,16 +2873,6 @@ The pair style does not have a single() function, so it can not be invoked by bond_style quartic. :dd -{Pair style does not support computing per-atom energy} :dt - -The pair style does not have a single() function, so it can -not be used to compute per-atom energy. :dd - -{Pair style does not support computing per-atom stress} :dt - -The pair style does not have a single() function, so it can -not be used to compute per-atom stress. :dd - {Pair style does not support pair_write} :dt The pair style does not have a single() function, so it can @@ -2735,6 +2976,14 @@ Self-explanatory. :dd +{Per-atom compute in equal-style variable formula} :dt + +Equal-style variables cannot use per-atom quantities. :dd + +{Per-atom fix in equal-style variable formula} :dt + +Equal-style variables cannot use per-atom quantities. :dd + {Potential file has duplicate entry} :dt The potential file for a SW or Tersoff potential has more than @@ -2750,20 +2999,9 @@ Granular potentials that include shear history effects can only be run with a newton setting where pairwise newton is "off". :dd -{Precompute ID for fix ave/spatial does not exist} :dt - -The compute used by fix ave/spatial requires a second pre-computation -compute, which isn't defined. :dd - -{Precompute ID for fix ave/time does not exist} :dt - -The compute used by fix ave/time requires a second pre-computation -compute, which isn't defined. :dd - -{Precompute ID of compute for fix ave/atom does not exist} :dt +{Power by 0 in variable formula} :dt -The compute specified for fix ave/atom has other pre-computed computes -defined. One of them does not exist. :dd +Self-explanatory. :dd {Press ID for fix nph does not exist} :dt @@ -2800,11 +3038,6 @@ A numeric error occurred in the creation of a rigid body by the fix rigid command. :dd -{Quotes in a single arg} :dt - -A single word should not be quoted in the input script; only a set of -words with intervening spaces should be quoted. :dd - {R0 < 0 for fix spring command} :dt Equilibrium spring length is invalid. :dd @@ -2955,10 +3188,14 @@ Self-explanatory. :dd -{Substitution for undefined variable} :dt +{Sqrt of negative in variable formula} :dt -The variable specified with a $ symbol in an input script command has -not been previously defined with a variable command. :dd +Self-explanatory. :dd + +{Substitution for illegal variable} :dt + +Input script line contained a variable that could not be substituted +for. :dd {TIP4P hydrogen has incorrect atom type} :dt @@ -3003,16 +3240,6 @@ The compute ID needed to compute temperature for the fix does not exist. :dd -{Temp ID of press ID for fix nph does not exist} :dt - -The compute ID needed to compute temperature within the pressure -compute ID for the fix does not exist. :dd - -{Temp ID of press ID for fix npt does not exist} :dt - -The compute ID needed to compute temperature within the pressure -compute ID for the fix does not exist. :dd - {Temper command before simulation box is defined} :dt The temper command cannot be used before a read_data, read_restart, or @@ -3022,6 +3249,11 @@ The region ID specified in the temperature command does not exist. :dd +{Tempering could not find thermo_pe compute} :dt + +This compute is created by the thermo command. It must have been +explicitly deleted by a uncompute command. :dd + {Tempering fix ID is not defined} :dt The fix ID specified by the temper command does not exist. :dd @@ -3031,6 +3263,11 @@ The fix specified by the temper command is not one that controls temperature (nvt or langevin). :dd +{Thermo and fix not computed at compatible times} :dt + +Fixes generate values on specific timesteps. The thermo output +does not match these timesteps. :dd + {Thermo compute ID does not compute scalar info} :dt The specified compute ID does not compute a scalar quantity @@ -3046,6 +3283,30 @@ The specified compute ID does not compute a large enough vector quantity for the requested index. :dd +{Thermo custom variable is not equal-style variable} :dt + +Only equal-style variables can be output with thermodynamics, not +atom-style variables. :dd + +{Thermo fix ID does not compute scalar info} :dt + +Only fixes that compute global values can be output with +thermodynamics. :dd + +{Thermo fix ID does not compute vector info} :dt + +Only fixes that compute global values can be output with +thermodynamics. :dd + +{Thermo fix ID vector is not large enough} :dt + +Index into vector is out of bounds. :dd + +{Thermo keyword in variable formula before initial run} :dt + +Calculating a thermo keyword before the first run is not allowed +because various quantities may not yet be initialized. :dd + {Thermo style does not use drot} :dt Cannot use thermo_modify to set this parameter since the thermo_style @@ -3079,16 +3340,6 @@ The thermo_style command cannot be used before a read_data, read_restart, or create_box command. :dd -{Thermodynamics must compute PE for temper command} :dt - -The thermo style must insure that thermodynamics computations include -potential energy when tempering is performed. :dd - -{Thermodynamics not computed on tempering swap steps} :dt - -The thermo command must insure that thermodynamics (including energy) -is computed on the timesteps that tempering swaps are attempted. :dd - {Timestep must be >= 0} :dt Specified timestep size is invalid. :dd @@ -3151,6 +3402,10 @@ length in that dimension. E.g. the xy tilt must be between -half and +half of the x box length. :dd +{Two groups cannot be the same in fix spring couple} :dt + +Self-explanatory. :dd + {Unbalanced quotes in input line} :dt No matching end double quote was found following a leading double @@ -3226,6 +3481,11 @@ Must use lattice command with compute fix deposit command if units option is set to lattice. :dd +{Use of fix dt/reset with undefined lattice} :dt + +Must use lattice command with fix dt/reset command if units option is +set to lattice. :dd + {Use of fix indent with undefined lattice} :dt The lattice command must be used to define a lattice before using the @@ -3255,27 +3515,32 @@ Self-explanatory. :dd -{Variable compute ID does not compute scalar info} :dt +{Variable name for compute sum does not exist} :dt -The specified compute ID does not compute a scalar quantity -as requested. :dd +Self-explanatory. :dd -{Variable compute ID vector is not large enough} :dt +{Variable name for fix ave/atom does not exist} :dt -The specified compute ID does not compute a large enough vector -quantity for the requested index. :dd +Self-explanatory. :dd -{Variable equal keyword used before initial run} :dt +{Variable name for fix ave/spatial does not exist} :dt -Cannot evaluate the variable at this stage of input script. :dd +Self-explanatory. :dd -{Variable equal keyword used before simulation box defined} :dt +{Variable name for fix ave/time does not exist} :dt -Cannot evaluate the variable at this stage of input script. :dd +Self-explanatory. :dd -{Variable group ID does not exist} :dt +{Variable name must be alphanumeric or underscore characters} :dt -A group specified in a variable definition does not exist. :dd +Self-explanatory. :dd + +{Variable uses compute via thermo keyword that thermo does not} :dt + +A thermo keyword used in the variable requires a particular compute be +invoked. However, the thermodynamic output does not use this compute, +so it has not been initialized. Include the keyword in thermodynamic +output. :dd {Velocity command before simulation box is defined} :dt @@ -3315,6 +3580,11 @@ :dlb +{Dump dcd/xtc timestamp may be wrong with fix dt/reset} :dt + +If the fix changes the timestep, the dump dcd file will not +reflect the change. :dd + {FENE bond too long: %d %d %d %g} :dt A FENE bond has stretched dangerously far. It's interaction strength @@ -3365,22 +3635,10 @@ It is not efficient to use compute coord/atom more than once. :dd -{More than one compute ebond/atom} :dt - -It is not efficient to use compute ebond/atom more than once. :dd - -{More than one compute epair/atom} :dt - -It is not efficient to use compute epair/atom more than once. :dd - {More than one compute ke/atom} :dt It is not efficient to use compute ke/atom more than once. :dd -{More than one compute stress/atom} :dt - -It is not efficient to use compute stress/atom more than once. :dd - {More than one fix msd} :dt It is not efficient to use fix msd more than once. :dd @@ -3563,8 +3821,4 @@ computed by integrating the density of a periodic system out to infinity. :dd -{Variable equal keyword used with non-current thermo} :dt - -The evaluation of the variable may be inaccurate as a result. :dd - :dle diff -Naur lammps-16Jan08/doc/fix_indent.html lammps-17Jan08/doc/fix_indent.html --- lammps-16Jan08/doc/fix_indent.html 2008-01-02 12:25:15.000000000 -0700 +++ lammps-17Jan08/doc/fix_indent.html 2008-01-17 10:02:19.000000000 -0700 @@ -94,6 +94,10 @@ velocity and force constant since they are defined in terms of distance as well.

+

IMPORTANT NOTE: You should insure the indenter's extent does not +overlap a periodic boundary, either for a fixed indenter, or one that +moves. No check for such overlaps is performed by the code. +

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

No information about this fix is written to binary restart diff -Naur lammps-16Jan08/doc/fix_indent.txt lammps-17Jan08/doc/fix_indent.txt --- lammps-16Jan08/doc/fix_indent.txt 2008-01-02 12:25:15.000000000 -0700 +++ lammps-17Jan08/doc/fix_indent.txt 2008-01-17 10:02:19.000000000 -0700 @@ -85,6 +85,10 @@ velocity and force constant since they are defined in terms of distance as well. +IMPORTANT NOTE: You should insure the indenter's extent does not +overlap a periodic boundary, either for a fixed indenter, or one that +moves. No check for such overlaps is performed by the code. + [Restart, fix_modify, output, run start/stop, minimize info:] No information about this fix is written to "binary restart diff -Naur lammps-16Jan08/lib/meam/meam_dens_final.F lammps-17Jan08/lib/meam/meam_dens_final.F --- lammps-16Jan08/lib/meam/meam_dens_final.F 2007-01-30 14:53:58.000000000 -0700 +++ lammps-17Jan08/lib/meam/meam_dens_final.F 2008-01-17 10:04:27.000000000 -0700 @@ -1,19 +1,22 @@ c Extern "C" declaration has the form: c -c void meam_dens_final_(int *, int *, int *, double *, int *, int *, int *, +c void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *, +c int *, int *, int *, c double *, double *, double *, double *, double *, double *, c double *, double *, double *, double *, double *, double *, c double *, double *, double *, double *, int *); c c Call from pair_meam.cpp has the form: c -c meam_dens_final_(&nlocal,&nmax,&eflag,&eng_vdwl,ntype,type,fmap, +c meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom, +c &eng_vdwl,eatom,ntype,type,fmap, c &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0], c &arho3b[0][0],&t_ave[0][0],gamma,dgamma1, c dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag); c - subroutine meam_dens_final(nlocal, nmax, eflag, eng_vdwl, + subroutine meam_dens_final(nlocal, nmax, + $ eflag_either, eflag_global, eflag_atom, eng_vdwl, eatom, $ ntype, type, fmap, $ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, $ Gamma, dGamma1, dGamma2, dGamma3, @@ -22,8 +25,9 @@ use meam_data implicit none - integer nlocal, nmax, eflag, ntype, type, fmap - real*8 eng_vdwl, Arho1, Arho2 + integer nlocal, nmax, eflag_either, eflag_global, eflag_atom + integer ntype, type, fmap + real*8 eng_vdwl, eatom, Arho1, Arho2 real*8 Arho2b, Arho3, Arho3b real*8 t_ave real*8 Gamma, dGamma1, dGamma2, dGamma3 @@ -31,6 +35,7 @@ real*8 fp integer errorflag + dimension eatom(nmax) dimension type(nmax), fmap(ntype) dimension Arho1(3,nmax), Arho2(6,nmax), Arho2b(nmax) dimension Arho3(10,nmax), Arho3b(3,nmax), t_ave(3,nmax) @@ -48,80 +53,87 @@ do i = 1,nlocal - elti = fmap(type(i)) - rho1(i) = 0.d0 - rho2(i) = -1.d0/3.d0*Arho2b(i)*Arho2b(i) - rho3(i) = 0.d0 - do m = 1,3 - rho1(i) = rho1(i) + Arho1(m,i)*Arho1(m,i) - rho3(i) = rho3(i) - 3.d0/5.d0*Arho3b(m,i)*Arho3b(m,i) - enddo - do m = 1,6 - rho2(i) = rho2(i) + v2D(m)*Arho2(m,i)*Arho2(m,i) - enddo - do m = 1,10 - rho3(i) = rho3(i) + v3D(m)*Arho3(m,i)*Arho3(m,i) - enddo - - if( rho0(i) .gt. 0.0 ) then - t_ave(1,i) = t_ave(1,i)/rho0(i) - t_ave(2,i) = t_ave(2,i)/rho0(i) - t_ave(3,i) = t_ave(3,i)/rho0(i) - endif - - Gamma(i) = t_ave(1,i)*rho1(i) - $ + t_ave(2,i)*rho2(i) + t_ave(3,i)*rho3(i) - - if( rho0(i) .gt. 0.0 ) then - Gamma(i) = Gamma(i)/(rho0(i)*rho0(i)) - end if - - call G_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,errorflag) - if (errorflag.ne.0) return - if (ibar_meam(elti).le.0) then - Gbar = 1.d0 - else - call get_shpfcn(shp,lattce_meam(elti,elti)) - gam = (t_ave(1,i)*shp(1)+t_ave(2,i)*shp(2) - $ +t_ave(3,i)*shp(3))/(Z_meam(elti)**2) - call G_gam(gam,ibar_meam(elti),gsmooth_factor, - $ Gbar,errorflag) - endif - rho(i) = rho0(i) * G - rhob = rho(i)/(rho0_meam(elti)*Z_meam(elti)*Gbar) - - Z = Z_meam(elti) - call dG_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,dG) - if (ibar_meam(elti).le.0) then - Gbar = 1.d0 - dGbar = 0.d0 - else - call get_shpfcn(shpi,lattce_meam(elti,elti)) - gam = (t_ave(1,i)*shpi(1)+t_ave(2,i)*shpi(2) - $ +t_ave(3,i)*shpi(3))/(Z*Z) - call dG_gam(gam,ibar_meam(elti),gsmooth_factor, - $ Gbar,dGbar) - endif - denom = 1.d0/(rho0_meam(elti)*Z*Gbar) - dGamma1(i) = (G - 2*dG*Gamma(i))*denom - - if( rho0(i) .ne. 0.d0 ) then - dGamma2(i) = (dG/rho0(i))*denom - else - dGamma2(i) = 0.d0 - end if + elti = fmap(type(i)) + rho1(i) = 0.d0 + rho2(i) = -1.d0/3.d0*Arho2b(i)*Arho2b(i) + rho3(i) = 0.d0 + do m = 1,3 + rho1(i) = rho1(i) + Arho1(m,i)*Arho1(m,i) + rho3(i) = rho3(i) - 3.d0/5.d0*Arho3b(m,i)*Arho3b(m,i) + enddo + do m = 1,6 + rho2(i) = rho2(i) + v2D(m)*Arho2(m,i)*Arho2(m,i) + enddo + do m = 1,10 + rho3(i) = rho3(i) + v3D(m)*Arho3(m,i)*Arho3(m,i) + enddo + + if( rho0(i) .gt. 0.0 ) then + t_ave(1,i) = t_ave(1,i)/rho0(i) + t_ave(2,i) = t_ave(2,i)/rho0(i) + t_ave(3,i) = t_ave(3,i)/rho0(i) + endif - dGamma3(i) = rho0(i)*G*dGbar/(Gbar*Z*Z)*denom + Gamma(i) = t_ave(1,i)*rho1(i) + $ + t_ave(2,i)*rho2(i) + t_ave(3,i)*rho3(i) + + if( rho0(i) .gt. 0.0 ) then + Gamma(i) = Gamma(i)/(rho0(i)*rho0(i)) + end if + + call G_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,errorflag) + if (errorflag.ne.0) return + if (ibar_meam(elti).le.0) then + Gbar = 1.d0 + else + call get_shpfcn(shp,lattce_meam(elti,elti)) + gam = (t_ave(1,i)*shp(1)+t_ave(2,i)*shp(2) + $ +t_ave(3,i)*shp(3))/(Z_meam(elti)**2) + call G_gam(gam,ibar_meam(elti),gsmooth_factor, + $ Gbar,errorflag) + endif + rho(i) = rho0(i) * G + rhob = rho(i)/(rho0_meam(elti)*Z_meam(elti)*Gbar) - B = A_meam(elti)*Ec_meam(elti,elti) + Z = Z_meam(elti) + call dG_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,dG) + if (ibar_meam(elti).le.0) then + Gbar = 1.d0 + dGbar = 0.d0 + else + call get_shpfcn(shpi,lattce_meam(elti,elti)) + gam = (t_ave(1,i)*shpi(1)+t_ave(2,i)*shpi(2) + $ +t_ave(3,i)*shpi(3))/(Z*Z) + call dG_gam(gam,ibar_meam(elti),gsmooth_factor, + $ Gbar,dGbar) + endif + denom = 1.d0/(rho0_meam(elti)*Z*Gbar) + dGamma1(i) = (G - 2*dG*Gamma(i))*denom - if( rhob .ne. 0.d0 ) then - fp(i) = B*(log(rhob)+1.d0) - if (eflag.eq.1) eng_vdwl = eng_vdwl + B*rhob*log(rhob) - else - fp(i) = B - endif - + if( rho0(i) .ne. 0.d0 ) then + dGamma2(i) = (dG/rho0(i))*denom + else + dGamma2(i) = 0.d0 + end if + + dGamma3(i) = rho0(i)*G*dGbar/(Gbar*Z*Z)*denom + + B = A_meam(elti)*Ec_meam(elti,elti) + + if( rhob .ne. 0.d0 ) then + fp(i) = B*(log(rhob)+1.d0) + if (eflag_either.ne.0) then + if (eflag_global.ne.0) then + eng_vdwl = eng_vdwl + B*rhob*log(rhob) + endif + if (eflag_atom.ne.0) then + eatom(i) = eatom(i) + B*rhob*log(rhob) + endif + endif + else + fp(i) = B + endif + enddo return diff -Naur lammps-16Jan08/lib/meam/meam_dens_init.F lammps-17Jan08/lib/meam/meam_dens_init.F --- lammps-16Jan08/lib/meam/meam_dens_init.F 2007-01-30 14:53:58.000000000 -0700 +++ lammps-17Jan08/lib/meam/meam_dens_init.F 2008-01-17 10:04:27.000000000 -0700 @@ -8,14 +8,14 @@ c c Call from pair_meam.cpp has the form: c -c meam_dens_init_(&i,&nmax,&eflag,&eng_vdwl,ntype,type,fmap,&x[0][0], +c meam_dens_init_(&i,&nmax,ntype,type,fmap,&x[0][0], c &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i], c &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], c rho0,&arho1[0][0],&arho2[0][0],arho2b, c &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&errorflag); c - subroutine meam_dens_init(i, nmax, eflag, eng_vdwl, + subroutine meam_dens_init(i, nmax, $ ntype, type, fmap, x, $ numneigh, firstneigh, $ numneigh_full, firstneigh_full, @@ -25,8 +25,7 @@ use meam_data implicit none - integer i, nmax, eflag, ntype, type, fmap - real*8 eng_vdwl + integer i, nmax, ntype, type, fmap real*8 x integer numneigh, firstneigh, numneigh_full, firstneigh_full real*8 scrfcn, dscrfcn, fcpair diff -Naur lammps-16Jan08/lib/meam/meam_force.F lammps-17Jan08/lib/meam/meam_force.F --- lammps-16Jan08/lib/meam/meam_force.F 2007-01-30 14:53:58.000000000 -0700 +++ lammps-17Jan08/lib/meam/meam_force.F 2008-01-17 10:04:27.000000000 -0700 @@ -8,36 +8,40 @@ c c Call from pair_meam.cpp has the form: c -c meam_force_(&i,&nmax,&eflag,&eng_vdwl,&ntype,type,fmap,&x[0][0], +c meam_force_(&i,&nmax,&eflag_either,&eflag_global,&eflag_atom,&vflag_atom, +c &eng_vdwl,eatom,&ntype,type,fmap,&x[0][0], c &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i], c &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], c dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop, c &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0], -c &t_ave[0][0],&f[0][0],&strssa[0][0][0],&errorflag); +c &t_ave[0][0],&f[0][0],&vatom[0][0],&errorflag); c - subroutine meam_force(i, nmax, eflag, eng_vdwl, - $ ntype, type, fmap, x, + subroutine meam_force(i, nmax, + $ eflag_either, eflag_global, eflag_atom, vflag_atom, + $ eng_vdwl, eatom, ntype, type, fmap, x, $ numneigh, firstneigh, numneigh_full, firstneigh_full, $ scrfcn, dscrfcn, fcpair, $ dGamma1, dGamma2, dGamma3, rho0, rho1, rho2, rho3, fp, - $ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, f, strssa, + $ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, f, vatom, $ errorflag) use meam_data implicit none - integer eflag, nmax, ntype, type, fmap - real*8 eng_vdwl, x + integer eflag_either, eflag_global, eflag_atom, vflag_atom + integer nmax, ntype, type, fmap + real*8 eng_vdwl, eatom, x integer numneigh, firstneigh, numneigh_full, firstneigh_full real*8 scrfcn, dscrfcn, fcpair real*8 dGamma1, dGamma2, dGamma3 real*8 rho0, rho1, rho2, rho3, fp real*8 Arho1, Arho2, Arho2b real*8 Arho3, Arho3b - real*8 t_ave, f, strssa + real*8 t_ave, f, vatom integer errorflag + dimension eatom(nmax) dimension type(nmax), fmap(ntype) dimension x(3,nmax) dimension firstneigh(numneigh), firstneigh_full(numneigh_full) @@ -46,13 +50,13 @@ dimension rho0(nmax), rho1(nmax), rho2(nmax), rho3(nmax), fp(nmax) dimension Arho1(3,nmax), Arho2(6,nmax), Arho2b(nmax) dimension Arho3(10,nmax), Arho3b(3,nmax) - dimension t_ave(3,nmax), f(3,nmax), strssa(3,3,nmax) + dimension t_ave(3,nmax), f(3,nmax), vatom(6,nmax) - integer i,j,jn,k,kn,kk,m,n,p,q,strscomp + integer i,j,jn,k,kn,kk,m,n,p,q integer nv2,nv3,elti,eltj,ind real*8 xitmp,yitmp,zitmp,delij(3),delref(3),rij2,rij,rij3 - real*8 delik(3),deljk(3) - real*8 Eu,astar,astarp + real*8 delik(3),deljk(3),v(6),fi(3),fj(3) + real*8 Eu,astar,astarp,third,sixth real*8 pp,phiforce,dUdrij,dUdsij,dUdrijm(3),force,forcem real*8 B,r,recip,phi,phip,rhop,a real*8 sij,fcij,dfcij,ds(3) @@ -82,8 +86,8 @@ real*8 t1i,t2i,t3i,t1j,t2j,t3j errorflag = 0 - - strscomp = 0 + third = 1.0/3.0 + sixth = 1.0/6.0 c Compute forces atom i @@ -122,7 +126,13 @@ $ + phirar4(kk,ind) recip = 1.0d0/r - if (eflag.eq.1) eng_vdwl = eng_vdwl + phi*sij + if (eflag_either.ne.0) then + if (eflag_global.ne.0) eng_vdwl = eng_vdwl + phi*sij + if (eflag_atom.ne.0) then + eatom(i) = eatom(i) + 0.5*phi*sij + eatom(j) = eatom(j) + 0.5*phi*sij + endif + endif c write(1,*) "force_meamf: phi: ",phi c write(1,*) "force_meamf: phip: ",phip @@ -371,23 +381,43 @@ enddo c Add the part of the force due to dUdrij and dUdsij + force = dUdrij*recip + dUdsij*dscrfcn(jn) do m = 1,3 forcem = delij(m)*force + dUdrijm(m) f(m,i) = f(m,i) + forcem f(m,j) = f(m,j) - forcem - if (strscomp>0) then - do n = 1,3 - arg = delij(n)*forcem - strssa(n,m,i) = strssa(n,m,i) + arg - strssa(n,m,j) = strssa(n,m,j) + arg - enddo - endif enddo -c write(1,*)"forcei1: ",f(1,i)," ",f(2,i)," ",f(3,i) -c write(1,*)"forcej1: ",f(1,j)," ",f(2,j)," ",f(3,j) + +c Tabulate per-atom virial as symmetrized stress tensor + + if (vflag_atom.ne.0) then + fi(1) = delij(1)*force + dUdrijm(1) + fi(2) = delij(2)*force + dUdrijm(2) + fi(3) = delij(3)*force + dUdrijm(3) + v(1) = -0.5 * (delij(1) * fi(1)) + v(2) = -0.5 * (delij(2) * fi(2)) + v(3) = -0.5 * (delij(3) * fi(3)) + v(4) = -0.25 * (delij(1)*fi(2) + delij(2)*fi(1)) + v(5) = -0.25 * (delij(1)*fi(3) + delij(3)*fi(1)) + v(6) = -0.25 * (delij(2)*fi(3) + delij(3)*fi(2)) + + vatom(1,i) = vatom(1,i) + v(1) + vatom(2,i) = vatom(2,i) + v(2) + vatom(3,i) = vatom(3,i) + v(3) + vatom(4,i) = vatom(4,i) + v(4) + vatom(5,i) = vatom(5,i) + v(5) + vatom(6,i) = vatom(6,i) + v(6) + vatom(1,j) = vatom(1,j) + v(1) + vatom(2,j) = vatom(2,j) + v(2) + vatom(3,j) = vatom(3,j) + v(3) + vatom(4,j) = vatom(4,j) + v(4) + vatom(5,j) = vatom(5,j) + v(5) + vatom(6,j) = vatom(6,j) + v(6) + endif c Now compute forces on other atoms k due to change in sij + if (sij.eq.0.d0.or.sij.eq.1.d0) goto 100 do kn = 1,numneigh_full k = firstneigh_full(kn) @@ -406,24 +436,47 @@ f(m,j) = f(m,j) + force2*deljk(m) f(m,k) = f(m,k) - force1*delik(m) $ - force2*deljk(m) - if (strscomp>0) then - do n = m,3 - arg1 = force1*delik(m)*delik(n) - arg2 = force2*deljk(m)*deljk(n) - strssa(m,n,i) = strssa(m,n,i) + arg1 - strssa(m,n,j) = strssa(m,n,j) + arg2 - strssa(m,n,k) = strssa(m,n,k) + - $ arg1 + arg2 - if (n.ne.m) then - strssa(n,m,i) = strssa(n,m,i) + arg1 - strssa(n,m,j) = strssa(n,m,j) + arg2 - strssa(n,m,k) = strssa(n,m,k) - $ + arg1 + arg2 - endif - enddo - endif enddo +c Tabulate per-atom virial as symmetrized stress tensor + + if (vflag_atom.ne.0) then + fi(1) = force1*delik(1) + fi(2) = force1*delik(2) + fi(3) = force1*delik(3) + fj(1) = force2*deljk(1) + fj(2) = force2*deljk(2) + fj(3) = force2*deljk(3) + v(1) = -third * (delik(1)*fi(1) + deljk(1)*fj(1)) + v(2) = -third * (delik(2)*fi(2) + deljk(2)*fj(2)) + v(3) = -third * (delik(3)*fi(3) + deljk(3)*fj(3)) + v(4) = -sixth * (delik(1)*fi(2) + deljk(1)*fj(2) + + $ delik(2)*fi(1) + deljk(2)*fj(1)) + v(5) = -sixth * (delik(1)*fi(3) + deljk(1)*fj(3) + + $ delik(3)*fi(1) + deljk(3)*fj(1)) + v(6) = -sixth * (delik(2)*fi(3) + deljk(2)*fj(3) + + $ delik(3)*fi(2) + deljk(3)*fj(2)) + + vatom(1,i) = vatom(1,i) + v(1) + vatom(2,i) = vatom(2,i) + v(2) + vatom(3,i) = vatom(3,i) + v(3) + vatom(4,i) = vatom(4,i) + v(4) + vatom(5,i) = vatom(5,i) + v(5) + vatom(6,i) = vatom(6,i) + v(6) + vatom(1,j) = vatom(1,j) + v(1) + vatom(2,j) = vatom(2,j) + v(2) + vatom(3,j) = vatom(3,j) + v(3) + vatom(4,j) = vatom(4,j) + v(4) + vatom(5,j) = vatom(5,j) + v(5) + vatom(6,j) = vatom(6,j) + v(6) + vatom(1,k) = vatom(1,k) + v(1) + vatom(2,k) = vatom(2,k) + v(2) + vatom(3,k) = vatom(3,k) + v(3) + vatom(4,k) = vatom(4,k) + v(4) + vatom(5,k) = vatom(5,k) + v(5) + vatom(6,k) = vatom(6,k) + v(6) + endif + endif endif c end of k loop diff -Naur lammps-16Jan08/potentials/Ni_smf7.eam lammps-17Jan08/potentials/Ni_smf7.eam --- lammps-16Jan08/potentials/Ni_smf7.eam 2007-06-11 10:54:20.000000000 -0600 +++ lammps-17Jan08/potentials/Ni_smf7.eam 2008-01-17 10:04:27.000000000 -0700 @@ -1,4 +1,4 @@ -Ni functions for NiPd alloy (exponential Z) +Ni functions for NiPd alloy (exponential Z) - S.M. Foiles, PRB, 32, 7685 (1985) 28 58.710 3.5200 FCC 500 4.0080160320641114e-04 500 9.6192384769538952e-03 4.8000000000000114e+00 0. -3.7126554861002070e-01 -6.0479555138022789e-01 -7.9881657783159099e-01 -9.7033414815713570e-01 diff -Naur lammps-16Jan08/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp lammps-17Jan08/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp --- lammps-16Jan08/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp 2007-12-21 08:05:18.000000000 -0700 +++ lammps-17Jan08/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp 2008-01-17 10:01:57.000000000 -0700 @@ -71,28 +71,29 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype,itable; + int n,vlist[6]; + int iH1,iH2,jH1,jH2; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; double fraction,table; double delx1,dely1,delz1,delx2,dely2,delz2,delx3,dely3,delz3; double r,r2inv,r6inv,forcecoul,forcelj,cforce,negforce; double factor_coul,factor_lj; double grij,expm2,prefactor,t,erfc; - int iH1,iH2,jH1,jH2; - double xiM[3],xjM[3]; + double xiM[3],xjM[3],fO[3],fH[3],v[6]; double *x1,*x2; - double fO[3],fH[3]; int *ilist,*jlist,*numneigh,**firstneigh; float rsq; int *int_rsq = (int *) &rsq; + // if vflag_global = 2, reset vflag as if vflag_global = 1 + // necessary since TIP4P cannot compute virial as F dot r + // due to find_M() finding bonded H atoms which are not near O atom + + if (vflag % 4 == 2) vflag = 1 + vflag/4 * 4; + evdwl = ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = 0; - - // error check (for now) - - if (eflag_atom || vflag_atom) - error->all("Pair style lj/cut/coul/long/tip4p does not yet support peratom energy/virial"); + else evflag = 0; double **f = atom->f; double **x = atom->x; @@ -105,23 +106,6 @@ int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - // grow and zero temporary force array if necessary - - if (vflag && atom->nmax > nmax) { - memory->destroy_2d_double_array(ftmp); - nmax = atom->nmax; - ftmp = memory->create_2d_double_array(nmax,3,"pair:ftmp"); - } - - if (vflag) { - for (i = 0; i < 6; i++) virialtmp[i] = 0.0; - for (i = 0; i < nall; i++) { - ftmp[i][0] = 0.0; - ftmp[i][1] = 0.0; - ftmp[i][2] = 0.0; - } - } - inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; @@ -179,8 +163,11 @@ evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; - } - } else evdwl = 0.0; + } else evdwl = 0.0; + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,forcelj,delx,dely,delz); + } // adjust rsq for off-site O charge(s) @@ -227,27 +214,25 @@ cforce = forcecoul * r2inv; // if i,j are not O atoms, force is applied directly - // if i or j are O atoms, force is on fictitious atoms - // spread force to all 3 atoms in water molecule + // if i or j are O atoms, force is on fictitious atom + // and is thus spread to all 3 atoms in water molecule // formulas due to Feenstra et al, J Comp Chem, 20, 786 (1999) + n = 0; + if (itype != typeO) { - if (vflag == 0) { - f[i][0] += delx * cforce; - f[i][1] += dely * cforce; - f[i][2] += delz * cforce; - - } else { - ftmp[i][0] += delx * cforce; - ftmp[i][1] += dely * cforce; - ftmp[i][2] += delz * cforce; - - virialtmp[0] += 0.5 * delx * delx * cforce; - virialtmp[1] += 0.5 * dely * dely * cforce; - virialtmp[2] += 0.5 * delz * delz * cforce; - virialtmp[3] += 0.5 * dely * delx * cforce; - virialtmp[4] += 0.5 * delz * delx * cforce; - virialtmp[5] += 0.5 * delz * dely * cforce; + f[i][0] += delx * cforce; + f[i][1] += dely * cforce; + f[i][2] += delz * cforce; + + if (vflag) { + v[0] = 0.5 * delx * delx * cforce; + v[1] = 0.5 * dely * dely * cforce; + v[2] = 0.5 * delz * delz * cforce; + v[3] = 0.5 * delx * dely * cforce; + v[4] = 0.5 * delx * delz * cforce; + v[5] = 0.5 * dely * delz * cforce; + vlist[n++] = i; } } else { @@ -259,32 +244,19 @@ fH[1] = alpha * (dely*cforce); fH[2] = alpha * (delz*cforce); - if (vflag == 0) { - f[i][0] += fO[0]; - f[i][1] += fO[1]; - f[i][2] += fO[2]; - - f[iH1][0] += fH[0]; - f[iH1][1] += fH[1]; - f[iH1][2] += fH[2]; + f[i][0] += fO[0]; + f[i][1] += fO[1]; + f[i][2] += fO[2]; + + f[iH1][0] += fH[0]; + f[iH1][1] += fH[1]; + f[iH1][2] += fH[2]; - f[iH2][0] += fH[0]; - f[iH2][1] += fH[1]; - f[iH2][2] += fH[2]; - - } else { - ftmp[i][0] += fO[0]; - ftmp[i][1] += fO[1]; - ftmp[i][2] += fO[2]; - - ftmp[iH1][0] += fH[0]; - ftmp[iH1][1] += fH[1]; - ftmp[iH1][2] += fH[2]; - - ftmp[iH2][0] += fH[0]; - ftmp[iH2][1] += fH[1]; - ftmp[iH2][2] += fH[2]; + f[iH2][0] += fH[0]; + f[iH2][1] += fH[1]; + f[iH2][2] += fH[2]; + if (vflag) { delx1 = x[i][0] - x2[0]; dely1 = x[i][1] - x2[1]; delz1 = x[i][2] - x2[2]; @@ -300,32 +272,33 @@ delz3 = x[iH2][2] - x2[2]; domain->minimum_image(delx3,dely3,delz3); - virialtmp[0] += 0.5 * (delx1 * fO[0] + (delx2 + delx3) * fH[0]); - virialtmp[1] += 0.5 * (dely1 * fO[1] + (dely2 + dely3) * fH[1]); - virialtmp[2] += 0.5 * (delz1 * fO[2] + (delz2 + delz3) * fH[2]); - virialtmp[3] += 0.5 * (dely1 * fO[0] + (dely2 + dely3) * fH[0]); - virialtmp[4] += 0.5 * (delz1 * fO[0] + (delz2 + delz3) * fH[0]); - virialtmp[5] += 0.5 * (delz1 * fO[1] + (delz2 + delz3) * fH[1]); + v[0] = 0.5 * (delx1 * fO[0] + (delx2 + delx3) * fH[0]); + v[1] = 0.5 * (dely1 * fO[1] + (dely2 + dely3) * fH[1]); + v[2] = 0.5 * (delz1 * fO[2] + (delz2 + delz3) * fH[2]); + v[3] = 0.5 * (dely1 * fO[0] + (dely2 + dely3) * fH[0]); + v[4] = 0.5 * (delz1 * fO[0] + (delz2 + delz3) * fH[0]); + v[5] = 0.5 * (delz1 * fO[1] + (delz2 + delz3) * fH[1]); + + vlist[n++] = i; + vlist[n++] = iH1; + vlist[n++] = iH2; } } if (jtype != typeO) { - if (vflag == 0) { - f[j][0] -= delx * cforce; - f[j][1] -= dely * cforce; - f[j][2] -= delz * cforce; - - } else { - ftmp[j][0] -= delx * cforce; - ftmp[j][1] -= dely * cforce; - ftmp[j][2] -= delz * cforce; - - virialtmp[0] += 0.5 * delx * delx * cforce; - virialtmp[1] += 0.5 * dely * dely * cforce; - virialtmp[2] += 0.5 * delz * delz * cforce; - virialtmp[3] += 0.5 * dely * delx * cforce; - virialtmp[4] += 0.5 * delz * delx * cforce; - virialtmp[5] += 0.5 * delz * dely * cforce; + f[j][0] -= delx * cforce; + f[j][1] -= dely * cforce; + f[j][2] -= delz * cforce; + + if (vflag) { + v[0] += 0.5 * delx * delx * cforce; + v[1] += 0.5 * dely * dely * cforce; + v[2] += 0.5 * delz * delz * cforce; + v[3] += 0.5 * delx * dely * cforce; + v[4] += 0.5 * delx * delz * cforce; + v[5] += 0.5 * dely * delz * cforce; + + vlist[n++] = j; } } else { @@ -339,32 +312,19 @@ fH[1] = alpha * (dely*negforce); fH[2] = alpha * (delz*negforce); - if (vflag != 2) { - f[j][0] += fO[0]; - f[j][1] += fO[1]; - f[j][2] += fO[2]; + f[j][0] += fO[0]; + f[j][1] += fO[1]; + f[j][2] += fO[2]; - f[jH1][0] += fH[0]; - f[jH1][1] += fH[1]; - f[jH1][2] += fH[2]; - - f[jH2][0] += fH[0]; - f[jH2][1] += fH[1]; - f[jH2][2] += fH[2]; - - } else { - ftmp[j][0] += fO[0]; - ftmp[j][1] += fO[1]; - ftmp[j][2] += fO[2]; - - ftmp[jH1][0] += fH[0]; - ftmp[jH1][1] += fH[1]; - ftmp[jH1][2] += fH[2]; - - ftmp[jH2][0] += fH[0]; - ftmp[jH2][1] += fH[1]; - ftmp[jH2][2] += fH[2]; + f[jH1][0] += fH[0]; + f[jH1][1] += fH[1]; + f[jH1][2] += fH[2]; + + f[jH2][0] += fH[0]; + f[jH2][1] += fH[1]; + f[jH2][2] += fH[2]; + if (vflag) { delx1 = x[j][0] - x1[0]; dely1 = x[j][1] - x1[1]; delz1 = x[j][2] - x1[2]; @@ -380,12 +340,16 @@ delz3 = x[jH2][2] - x1[2]; domain->minimum_image(delx3,dely3,delz3); - virialtmp[0] += 0.5 * (delx1 * fO[0] + (delx2 + delx3) * fH[0]); - virialtmp[1] += 0.5 * (dely1 * fO[1] + (dely2 + dely3) * fH[1]); - virialtmp[2] += 0.5 * (delz1 * fO[2] + (delz2 + delz3) * fH[2]); - virialtmp[3] += 0.5 * (dely1 * fO[0] + (dely2 + dely3) * fH[0]); - virialtmp[4] += 0.5 * (delz1 * fO[0] + (delz2 + delz3) * fH[0]); - virialtmp[5] += 0.5 * (delz1 * fO[1] + (delz2 + delz3) * fH[1]); + v[0] += 0.5 * (delx1 * fO[0] + (delx2 + delx3) * fH[0]); + v[1] += 0.5 * (dely1 * fO[1] + (dely2 + dely3) * fH[1]); + v[2] += 0.5 * (delz1 * fO[2] + (delz2 + delz3) * fH[2]); + v[3] += 0.5 * (dely1 * fO[0] + (dely2 + dely3) * fH[0]); + v[4] += 0.5 * (delz1 * fO[0] + (delz2 + delz3) * fH[0]); + v[5] += 0.5 * (delz1 * fO[1] + (delz2 + delz3) * fH[1]); + + vlist[n++] = j; + vlist[n++] = jH1; + vlist[n++] = jH2; } } @@ -397,24 +361,13 @@ ecoul = qtmp*q[j] * table; } if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; - } - } else ecoul = 0.0; + } else ecoul = 0.0; - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,ecoul,fpair,delx,dely,delz); + if (evflag) ev_tally_list(n,vlist,ecoul,v); + } } } } - - if (vflag_fdotr) { - virial_compute(); - for (int i = 0; i < 6; i++) virial[i] += virialtmp[i]; - for (int i = 0; i < nall; i++) { - f[i][0] += ftmp[i][0]; - f[i][1] += ftmp[i][1]; - f[i][2] += ftmp[i][2]; - } - } } /* ---------------------------------------------------------------------- @@ -459,68 +412,14 @@ error->all("Pair style lj/cut/coul/long/tip4p requires newton pair on"); if (!atom->q_flag) error->all("Pair style lj/cut/coul/long/tip4p requires atom attribute q"); - - // request regular or rRESPA neighbor lists - - int irequest; - - if (update->whichflag == 0 && strcmp(update->integrate_style,"respa") == 0) { - int respa = 0; - if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; - if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; - - if (respa == 0) irequest = neighbor->request(this); - else if (respa == 1) { - irequest = neighbor->request(this); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->respaouter = 1; - } else { - irequest = neighbor->request(this); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this); - neighbor->requests[irequest]->id = 2; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->respamiddle = 1; - irequest = neighbor->request(this); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->respaouter = 1; - } - - } else irequest = neighbor->request(this); - - cut_coulsq = cut_coul * cut_coul; - - // set & error check interior rRESPA cutoffs - - if (strcmp(update->integrate_style,"respa") == 0) { - if (((Respa *) update->integrate)->level_inner >= 0) { - cut_respa = ((Respa *) update->integrate)->cutoff; - for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) - if (MIN(cut_lj[i][j],cut_coul) < cut_respa[3]) - error->all("Pair cutoff < Respa interior cutoff"); - } - } else cut_respa = NULL; - - // insure use of correct KSpace long-range solver, set g_ewald - - if (force->kspace == NULL) + if (strcmp(force->kspace_style,"pppm/tip4p") != 0) error->all("Pair style is incompatible with KSpace style"); - if (strcmp(force->kspace_style,"pppm/tip4p") == 0) - g_ewald = force->kspace->g_ewald; - else error->all("Pair style is incompatible with KSpace style"); - - // setup force tables + if (force->bond == NULL) + error->all("Must use a bond style with TIP4P potential"); + if (force->angle == NULL) + error->all("Must use an angle style with TIP4P potential"); - if (ncoultablebits) init_tables(); + PairLJCutCoulLong::init_style(); // set alpha parameter diff -Naur lammps-16Jan08/src/MANYBODY/pair_airebo.cpp lammps-17Jan08/src/MANYBODY/pair_airebo.cpp --- lammps-16Jan08/src/MANYBODY/pair_airebo.cpp 2008-01-08 16:34:19.000000000 -0700 +++ lammps-17Jan08/src/MANYBODY/pair_airebo.cpp 2008-01-17 10:01:57.000000000 -0700 @@ -90,11 +90,6 @@ if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; - // error check (for now) - - if (eflag_atom || vflag_atom) - error->all("Pair style airebo does not yet support peratom energy/virial"); - REBO_neigh(); FREBO(eflag,vflag); if (ljflag) FLJ(eflag,vflag); @@ -507,7 +502,7 @@ del[0] = delx; del[1] = dely; del[2] = delz; - bij = bondorder(i,j,del,rij,VA,f); + bij = bondorder(i,j,del,rij,VA,f,vflag_atom); dVAdi = bij*dVA; fpair = -(dVRdi+dVAdi) / rij; @@ -537,15 +532,16 @@ int testpath,npath,done; double evdwl,fpair; double rsq,best,wik,wkm,cij,rij,dwij,dwik,dwkj,dwkm,dwmj; - double delij[3],rijsq,delik[3],rik,delkj[3]; + double delij[3],rijsq,delik[3],rik,deljk[3]; double rkj,wkj,dC,VLJ,dVLJ,VA,Str,dStr,Stb; - double delkm[3],rkm,delmj[3],rmj,wmj,r2inv,r6inv,scale,delscale[3]; + double delkm[3],rkm,deljm[3],rmj,wmj,r2inv,r6inv,scale,delscale[3]; int *ilist,*jlist,*numneigh,**firstneigh; int *REBO_neighs_i,*REBO_neighs_k; - double delikS[3],delkjS[3],delkmS[3],delmjS[3]; + double delikS[3],deljkS[3],delkmS[3],deljmS[3],delimS[3]; double rikS,rkjS,rkmS,rmjS,wikS,dwikS; double wkjS,dwkjS,wkmS,dwkmS,wmjS,dwmjS; double fpair1,fpair2,fpair3; + double fi[3],fj[3],fk[3],fm[3]; // I-J interaction from full neighbor list // skip 1/2 of interactions since only consider each pair once @@ -629,10 +625,10 @@ } else wik = 0.0; if (wik > best) { - delkj[0] = x[k][0] - x[j][0]; - delkj[1] = x[k][1] - x[j][1]; - delkj[2] = x[k][2] - x[j][2]; - rsq = delkj[0]*delkj[0] + delkj[1]*delkj[1] + delkj[2]*delkj[2]; + deljk[0] = x[j][0] - x[k][0]; + deljk[1] = x[j][1] - x[k][1]; + deljk[2] = x[j][2] - x[k][2]; + rsq = deljk[0]*deljk[0] + deljk[1]*deljk[1] + deljk[2]*deljk[2]; if (rsq < rcmaxsq[ktype][jtype]) { rkj = sqrt(rsq); wkj = Sp(rkj,rcmin[ktype][jtype],rcmax[ktype][jtype],dwkj); @@ -646,9 +642,9 @@ rikS = rik; wikS = wik; dwikS = dwik; - delkjS[0] = delkj[0]; - delkjS[1] = delkj[1]; - delkjS[2] = delkj[2]; + deljkS[0] = deljk[0]; + deljkS[1] = deljk[1]; + deljkS[2] = deljk[2]; rkjS = rkj; wkjS = wkj; dwkjS = dwkj; @@ -679,11 +675,11 @@ } else wkm = 0.0; if (wik*wkm > best) { - delmj[0] = x[m][0] - x[j][0]; - delmj[1] = x[m][1] - x[j][1]; - delmj[2] = x[m][2] - x[j][2]; - rsq = delmj[0]*delmj[0] + delmj[1]*delmj[1] + - delmj[2]*delmj[2]; + deljm[0] = x[j][0] - x[m][0]; + deljm[1] = x[j][1] - x[m][1]; + deljm[2] = x[j][2] - x[m][2]; + rsq = deljm[0]*deljm[0] + deljm[1]*deljm[1] + + deljm[2]*deljm[2]; if (rsq < rcmaxsq[mtype][jtype]) { rmj = sqrt(rsq); wmj = Sp(rmj,rcmin[mtype][jtype],rcmax[mtype][jtype],dwmj); @@ -704,9 +700,9 @@ rkmS = rkm; wkmS = wkm; dwkmS = dwkm; - delmjS[0] = delmj[0]; - delmjS[1] = delmj[1]; - delmjS[2] = delmj[2]; + deljmS[0] = deljm[0]; + deljmS[1] = deljm[1]; + deljmS[2] = deljm[2]; rmjS = rmj; wmjS = wmj; dwmjS = dwmj; @@ -739,7 +735,8 @@ delscale[0] = scale * delij[0]; delscale[1] = scale * delij[1]; delscale[2] = scale * delij[2]; - Stb = bondorderLJ(i,j,delscale,rcmin[itype][jtype],VA,delij,rij,f); + Stb = bondorderLJ(i,j,delscale,rcmin[itype][jtype],VA, + delij,rij,f,vflag_atom); } else Stb = 0.0; fpair = -(dStr * (Stb*cij*VLJ - cij*VLJ) + @@ -766,50 +763,71 @@ f[atomj][0] -= delij[0]*fpair; f[atomj][1] -= delij[1]*fpair; f[atomj][2] -= delij[2]*fpair; - if (vflag_atom) ev_tally(atomi,atomj,nlocal,newton_pair, - 0.0,0.0,fpair,delij[0],delij[1],delij[2]); + + if (vflag_atom) v_tally2(atomi,atomj,fpair,delij); + } else if (npath == 3) { fpair1 = dC*dwikS*wkjS / rikS; - f[atomi][0] += delikS[0]*fpair1; - f[atomi][1] += delikS[1]*fpair1; - f[atomi][2] += delikS[2]*fpair1; - f[atomk][0] -= delikS[0]*fpair1; - f[atomk][1] -= delikS[1]*fpair1; - f[atomk][2] -= delikS[2]*fpair1; + fi[0] = delikS[0]*fpair1; + fi[1] = delikS[1]*fpair1; + fi[2] = delikS[2]*fpair1; fpair2 = dC*wikS*dwkjS / rkjS; - f[atomk][0] += delkjS[0]*fpair2; - f[atomk][1] += delkjS[1]*fpair2; - f[atomk][2] += delkjS[2]*fpair2; - f[atomj][0] -= delkjS[0]*fpair2; - f[atomj][1] -= delkjS[1]*fpair2; - f[atomj][2] -= delkjS[2]*fpair2; + fj[0] = deljkS[0]*fpair2; + fj[1] = deljkS[1]*fpair2; + fj[2] = deljkS[2]*fpair2; + + f[atomi][0] += fi[0]; + f[atomi][1] += fi[1]; + f[atomi][2] += fi[2]; + f[atomj][0] += fj[0]; + f[atomj][1] += fj[1]; + f[atomj][2] += fj[2]; + f[atomk][0] -= fi[0] + fj[0]; + f[atomk][1] -= fi[1] + fj[1]; + f[atomk][2] -= fi[2] + fj[2]; + if (vflag_atom) - v_tally3(atomi,atomj,atomk,fpair1,fpair2,delikS,delkjS); + v_tally3(atomi,atomj,atomk,fi,fj,delikS,deljkS); + } else { fpair1 = dC*dwikS*wkmS*wmjS / rikS; - f[atomi][0] += delikS[0]*fpair1; - f[atomi][1] += delikS[1]*fpair1; - f[atomi][2] += delikS[2]*fpair1; - f[atomk][0] -= delikS[0]*fpair1; - f[atomk][1] -= delikS[1]*fpair1; - f[atomk][2] -= delikS[2]*fpair1; + fi[0] = delikS[0]*fpair1; + fi[1] = delikS[1]*fpair1; + fi[2] = delikS[2]*fpair1; + fpair2 = dC*wikS*dwkmS*wmjS / rkmS; - f[atomk][0] += delkmS[0]*fpair2; - f[atomk][1] += delkmS[1]*fpair2; - f[atomk][2] += delkmS[2]*fpair2; - f[atomm][0] -= delkmS[0]*fpair2; - f[atomm][1] -= delkmS[1]*fpair2; - f[atomm][2] -= delkmS[2]*fpair2; + fk[0] = delkmS[0]*fpair2 - fi[0]; + fk[1] = delkmS[1]*fpair2 - fi[1]; + fk[2] = delkmS[2]*fpair2 - fi[2]; + fpair3 = dC*wikS*wkmS*dwmjS / rmjS; - f[atomm][0] += delmjS[0]*fpair3; - f[atomm][1] += delmjS[1]*fpair3; - f[atomm][2] += delmjS[2]*fpair3; - f[atomj][0] -= delmjS[0]*fpair3; - f[atomj][1] -= delmjS[1]*fpair3; - f[atomj][2] -= delmjS[2]*fpair3; - if (vflag_atom) - v_tally4(atomi,atomj,atomk,atomm, - fpair1,fpair2,fpair3,delikS,delkmS,delmjS); + fj[0] = deljmS[0]*fpair3; + fj[1] = deljmS[1]*fpair3; + fj[2] = deljmS[2]*fpair3; + + fm[0] = -delkmS[0]*fpair2 - fj[0]; + fm[1] = -delkmS[1]*fpair2 - fj[1]; + fm[2] = -delkmS[2]*fpair2 - fj[2]; + + f[atomi][0] += fi[0]; + f[atomi][1] += fi[1]; + f[atomi][2] += fi[2]; + f[atomj][0] += fj[0]; + f[atomj][1] += fj[1]; + f[atomj][2] += fj[2]; + f[atomk][0] += fk[0]; + f[atomk][1] += fk[1]; + f[atomk][2] += fk[2]; + f[atomm][0] += fm[0]; + f[atomm][1] += fm[1]; + f[atomm][2] += fm[2]; + + if (vflag_atom) { + delimS[0] = delikS[0] + delkmS[0]; + delimS[1] = delikS[1] + delkmS[1]; + delimS[2] = delikS[2] + delkmS[2]; + v_tally4(atomi,atomj,atomk,atomm,fi,fj,fk,delimS,deljmS,delkmS); + } } } } @@ -839,13 +857,11 @@ double dxidij,dxidik,dxidjk,dxjdji,dxjdjl,dxjdil; double ddndij,ddndik,ddndjk,ddndjl,ddndil,dcwddn,dcwdn,dvpdcw,Ftmp[3]; double del32[3],rsq,r32,del23[3],del21[3],r21; - double deljk[3],del34[3],delil[3],r23,r34; + double deljk[3],del34[3],delil[3],delkl[3],r23,r34; double fi[3],fj[3],fk[3],fl[3]; int itype,jtype,ktype,ltype,kk,ll,jj; int *REBO_neighs_i,*REBO_neighs_j; - evdwl = 0.0; - double **x = atom->x; double **f = atom->f; int *type = atom->type; @@ -1146,7 +1162,12 @@ f[k][0] += fk[0]; f[k][1] += fk[1]; f[k][2] += fk[2]; f[l][0] += fl[0]; f[l][1] += fl[1]; f[l][2] += fl[2]; - //if (evflag) ev_tally4(i,j,k,l,fi,fj,fk,fl); + if (evflag) { + delkl[0] = delil[0] - del21[0]; + delkl[1] = delil[1] - del21[1]; + delkl[2] = delil[2] - del21[2]; + ev_tally4(i,j,k,l,evdwl,fi,fj,fk,delil,del34,delkl); + } } } } @@ -1212,11 +1233,12 @@ ------------------------------------------------------------------------- */ double PairAIREBO::bondorder(int i, int j, double rij[3], - double rijmag, double VA, double **f) + double rijmag, double VA, + double **f, int vflag_atom) { int atomi,atomj,k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4; int itype,jtype,ktype,ltype,ntype; - double rik[3], rjl[3], rkn[3],rknmag,dNki,dwjl,bij; + double rik[3],rjl[3],rkn[3],rji[3],rki[3],rlj[3],rknmag,dNki,dwjl,bij; double NijC,NijH,NjiC,NjiH,wik,dwik,dwkn,wjl; double rikmag,rjlmag,cosjik,cosijl,g,tmp2,tmp3; double Etmp,pij,tmp,wij,dwij,NconjtmpI,NconjtmpJ,Nki,Nlj,dS; @@ -1225,7 +1247,7 @@ double dN2[2],dN3[3]; double dcosjikdrj[3],dcosijldrj[3],dcosijldrl[3]; double Tij; - double r32[3],r32mag,cos321; + double r32[3],r32mag,cos321,r43[3],r13[3]; double dNlj; double om1234,rln[3]; double rlnmag,dwln,r23[3],r23mag,r21[3],r21mag; @@ -1236,7 +1258,8 @@ double sin321,sin234,rr,rijrik,rijrjl,rjk2,rik2,ril2,rjl2; double dctik,dctjk,dctjl,dctij,dctji,dctil,rik2i,rjl2i,sink2i,sinl2i; double rjk[3],ril[3],dt1dik,dt1djk,dt1djl,dt1dil,dt1dij; - double F23[3],F12[3],F34[3],F31[3],F24[3]; + double F23[3],F12[3],F34[3],F31[3],F24[3],fi[3],fj[3],fk[3],fl[3]; + double f1[3],f2[3],f3[3],f4[4]; double dcut321,PijS,PjiS; double rij2,tspjik,dtsjik,tspijl,dtsijl,costmp; int *REBO_neighs,*REBO_neighs_i,*REBO_neighs_j,*REBO_neighs_k,*REBO_neighs_l; @@ -1338,58 +1361,68 @@ g = gSpline(cosjik,(NijC+NijH),itype,&dgdc,&dgdN); tmp2 = VA*.5*(tmp*wik*dgdc*exp(lamdajik)); - f[atomj][0] -= tmp2*dcosjikdrj[0]; - f[atomj][1] -= tmp2*dcosjikdrj[1]; - f[atomj][2] -= tmp2*dcosjikdrj[2]; - f[atomi][0] -= tmp2*dcosjikdri[0]; - f[atomi][1] -= tmp2*dcosjikdri[1]; - f[atomi][2] -= tmp2*dcosjikdri[2]; - f[atomk][0] -= tmp2*dcosjikdrk[0]; - f[atomk][1] -= tmp2*dcosjikdrk[1]; - f[atomk][2] -= tmp2*dcosjikdrk[2]; + fj[0] = -tmp2*dcosjikdrj[0]; + fj[1] = -tmp2*dcosjikdrj[1]; + fj[2] = -tmp2*dcosjikdrj[2]; + fi[0] = -tmp2*dcosjikdri[0]; + fi[1] = -tmp2*dcosjikdri[1]; + fi[2] = -tmp2*dcosjikdri[2]; + fk[0] = -tmp2*dcosjikdrk[0]; + fk[1] = -tmp2*dcosjikdrk[1]; + fk[2] = -tmp2*dcosjikdrk[2]; tmp2 = VA*.5*(tmp*wik*g*exp(lamdajik)*4.0*kronecker(itype,1)); - f[atomj][0] -= tmp2*(-rij[0]/rijmag); - f[atomj][1] -= tmp2*(-rij[1]/rijmag); - f[atomj][2] -= tmp2*(-rij[2]/rijmag); - f[atomi][0] -= tmp2*((-rik[0]/rikmag)+(rij[0]/rijmag)); - f[atomi][1] -= tmp2*((-rik[1]/rikmag)+(rij[1]/rijmag)); - f[atomi][2] -= tmp2*((-rik[2]/rikmag)+(rij[2]/rijmag)); - f[atomk][0] -= tmp2*(rik[0]/rikmag); - f[atomk][1] -= tmp2*(rik[1]/rikmag); - f[atomk][2] -= tmp2*(rik[2]/rikmag); + fj[0] -= tmp2*(-rij[0]/rijmag); + fj[1] -= tmp2*(-rij[1]/rijmag); + fj[2] -= tmp2*(-rij[2]/rijmag); + fi[0] -= tmp2*((-rik[0]/rikmag)+(rij[0]/rijmag)); + fi[1] -= tmp2*((-rik[1]/rikmag)+(rij[1]/rijmag)); + fi[2] -= tmp2*((-rik[2]/rikmag)+(rij[2]/rijmag)); + fk[0] -= tmp2*(rik[0]/rikmag); + fk[1] -= tmp2*(rik[1]/rikmag); + fk[2] -= tmp2*(rik[2]/rikmag); // coordination forces // dwik forces tmp2 = VA*.5*(tmp*dwik*g*exp(lamdajik))/rikmag; - f[atomi][0] -= tmp2*rik[0]; - f[atomi][1] -= tmp2*rik[1]; - f[atomi][2] -= tmp2*rik[2]; - f[atomk][0] += tmp2*rik[0]; - f[atomk][1] += tmp2*rik[1]; - f[atomk][2] += tmp2*rik[2]; + fi[0] -= tmp2*rik[0]; + fi[1] -= tmp2*rik[1]; + fi[2] -= tmp2*rik[2]; + fk[0] += tmp2*rik[0]; + fk[1] += tmp2*rik[1]; + fk[2] += tmp2*rik[2]; // PIJ forces tmp2 = VA*.5*(tmp*dN2[ktype]*dwik)/rikmag; - f[atomi][0] -= tmp2*rik[0]; - f[atomi][1] -= tmp2*rik[1]; - f[atomi][2] -= tmp2*rik[2]; - f[atomk][0] += tmp2*rik[0]; - f[atomk][1] += tmp2*rik[1]; - f[atomk][2] += tmp2*rik[2]; + fi[0] -= tmp2*rik[0]; + fi[1] -= tmp2*rik[1]; + fi[2] -= tmp2*rik[2]; + fk[0] += tmp2*rik[0]; + fk[1] += tmp2*rik[1]; + fk[2] += tmp2*rik[2]; // dgdN forces tmp2 = VA*.5*(tmp*tmp3*dwik)/rikmag; - f[atomi][0] -= tmp2*rik[0]; - f[atomi][1] -= tmp2*rik[1]; - f[atomi][2] -= tmp2*rik[2]; - f[atomk][0] += tmp2*rik[0]; - f[atomk][1] += tmp2*rik[1]; - f[atomk][2] += tmp2*rik[2]; + fi[0] -= tmp2*rik[0]; + fi[1] -= tmp2*rik[1]; + fi[2] -= tmp2*rik[2]; + fk[0] += tmp2*rik[0]; + fk[1] += tmp2*rik[1]; + fk[2] += tmp2*rik[2]; + + f[atomi][0] += fi[0]; f[atomi][1] += fi[1]; f[atomi][2] += fi[2]; + f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; + f[atomk][0] += fk[0]; f[atomk][1] += fk[1]; f[atomk][2] += fk[2]; + + if (vflag_atom) { + rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2]; + rki[0] = -rik[0]; rki[1] = -rik[1]; rki[2] = -rik[2]; + v_tally3(atomi,atomj,atomk,fj,fk,rji,rki); + } } } @@ -1470,58 +1503,67 @@ g = gSpline(cosijl,NjiC+NjiH,jtype,&dgdc,&dgdN); tmp2 = VA*.5*(tmp*wjl*dgdc*exp(lamdaijl)); - f[atomi][0] -= tmp2*dcosijldri[0]; - f[atomi][1] -= tmp2*dcosijldri[1]; - f[atomi][2] -= tmp2*dcosijldri[2]; - f[atomj][0] -= tmp2*dcosijldrj[0]; - f[atomj][1] -= tmp2*dcosijldrj[1]; - f[atomj][2] -= tmp2*dcosijldrj[2]; - f[atoml][0] -= tmp2*dcosijldrl[0]; - f[atoml][1] -= tmp2*dcosijldrl[1]; - f[atoml][2] -= tmp2*dcosijldrl[2]; + fi[0] = -tmp2*dcosijldri[0]; + fi[1] = -tmp2*dcosijldri[1]; + fi[2] = -tmp2*dcosijldri[2]; + fj[0] = -tmp2*dcosijldrj[0]; + fj[1] = -tmp2*dcosijldrj[1]; + fj[2] = -tmp2*dcosijldrj[2]; + fl[0] = -tmp2*dcosijldrl[0]; + fl[1] = -tmp2*dcosijldrl[1]; + fl[2] = -tmp2*dcosijldrl[2]; tmp2 = VA*.5*(tmp*wjl*g*exp(lamdaijl)*4.0*kronecker(jtype,1)); - f[atomi][0] -= tmp2*(rij[0]/rijmag); - f[atomi][1] -= tmp2*(rij[1]/rijmag); - f[atomi][2] -= tmp2*(rij[2]/rijmag); - f[atomj][0] -= tmp2*((-rjl[0]/rjlmag)-(rij[0]/rijmag)); - f[atomj][1] -= tmp2*((-rjl[1]/rjlmag)-(rij[1]/rijmag)); - f[atomj][2] -= tmp2*((-rjl[2]/rjlmag)-(rij[2]/rijmag)); - f[atoml][0] -= tmp2*(rjl[0]/rjlmag); - f[atoml][1] -= tmp2*(rjl[1]/rjlmag); - f[atoml][2] -= tmp2*(rjl[2]/rjlmag); + fi[0] -= tmp2*(rij[0]/rijmag); + fi[1] -= tmp2*(rij[1]/rijmag); + fi[2] -= tmp2*(rij[2]/rijmag); + fj[0] -= tmp2*((-rjl[0]/rjlmag)-(rij[0]/rijmag)); + fj[1] -= tmp2*((-rjl[1]/rjlmag)-(rij[1]/rijmag)); + fj[2] -= tmp2*((-rjl[2]/rjlmag)-(rij[2]/rijmag)); + fl[0] -= tmp2*(rjl[0]/rjlmag); + fl[1] -= tmp2*(rjl[1]/rjlmag); + fl[2] -= tmp2*(rjl[2]/rjlmag); // coordination forces // dwik forces tmp2 = VA*.5*(tmp*dwjl*g*exp(lamdaijl))/rjlmag; - f[atomj][0] -= tmp2*rjl[0]; - f[atomj][1] -= tmp2*rjl[1]; - f[atomj][2] -= tmp2*rjl[2]; - f[atoml][0] += tmp2*rjl[0]; - f[atoml][1] += tmp2*rjl[1]; - f[atoml][2] += tmp2*rjl[2]; + fj[0] -= tmp2*rjl[0]; + fj[1] -= tmp2*rjl[1]; + fj[2] -= tmp2*rjl[2]; + fl[0] += tmp2*rjl[0]; + fl[1] += tmp2*rjl[1]; + fl[2] += tmp2*rjl[2]; // PIJ forces tmp2 = VA*.5*(tmp*dN2[ltype]*dwjl)/rjlmag; - f[atomj][0] -= tmp2*rjl[0]; - f[atomj][1] -= tmp2*rjl[1]; - f[atomj][2] -= tmp2*rjl[2]; - f[atoml][0] += tmp2*rjl[0]; - f[atoml][1] += tmp2*rjl[1]; - f[atoml][2] += tmp2*rjl[2]; + fj[0] -= tmp2*rjl[0]; + fj[1] -= tmp2*rjl[1]; + fj[2] -= tmp2*rjl[2]; + fl[0] += tmp2*rjl[0]; + fl[1] += tmp2*rjl[1]; + fl[2] += tmp2*rjl[2]; // dgdN forces tmp2 = VA*.5*(tmp*tmp3*dwjl)/rjlmag; - f[atomj][0] -= tmp2*rjl[0]; - f[atomj][1] -= tmp2*rjl[1]; - f[atomj][2] -= tmp2*rjl[2]; - f[atoml][0] += tmp2*rjl[0]; - f[atoml][1] += tmp2*rjl[1]; - f[atoml][2] += tmp2*rjl[2]; + fj[0] -= tmp2*rjl[0]; + fj[1] -= tmp2*rjl[1]; + fj[2] -= tmp2*rjl[2]; + fl[0] += tmp2*rjl[0]; + fl[1] += tmp2*rjl[1]; + fl[2] += tmp2*rjl[2]; + + f[atomi][0] += fi[0]; f[atomi][1] += fi[1]; f[atomi][2] += fi[2]; + f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; + f[atoml][0] += fl[0]; f[atoml][1] += fl[1]; f[atoml][2] += fl[2]; + + if (vflag_atom) { + rlj[0] = -rjl[0]; rlj[1] = -rjl[1]; rlj[2] = -rjl[2]; + v_tally3(atomi,atomj,atoml,fi,fl,rij,rlj); + } } } @@ -1554,6 +1596,8 @@ f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; + if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + tmp2 = VA*dN3[2]*(2.0*NconjtmpI*dwik*SpN)/rikmag; f[atomi][0] -= tmp2*rik[0]; f[atomi][1] -= tmp2*rik[1]; @@ -1562,6 +1606,8 @@ f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; + if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + if (fabs(dNki) > TOL) { REBO_neighs_k = REBO_firstneigh[atomk]; for (n = 0; n < REBO_numneigh[atomk]; n++) { @@ -1581,6 +1627,8 @@ f[atomn][0] += tmp2*rkn[0]; f[atomn][1] += tmp2*rkn[1]; f[atomn][2] += tmp2*rkn[2]; + + if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn); } } } @@ -1611,6 +1659,8 @@ f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; + if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + tmp2 = VA*dN3[2]*(2.0*NconjtmpJ*dwjl*SpN)/rjlmag; f[atomj][0] -= tmp2*rjl[0]; f[atomj][1] -= tmp2*rjl[1]; @@ -1619,6 +1669,8 @@ f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; + if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + if (fabs(dNlj) > TOL) { REBO_neighs_l = REBO_firstneigh[atoml]; for (n = 0; n < REBO_numneigh[atoml]; n++) { @@ -1638,6 +1690,8 @@ f[atomn][0] += tmp2*rln[0]; f[atomn][1] += tmp2*rln[1]; f[atomn][2] += tmp2*rln[2]; + + if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln); } } } @@ -1809,38 +1863,53 @@ F24[1] = (fcilpc*ril[1]); F24[2] = (fcilpc*ril[2]); - f[atom1][0] +=-F12[0]-F31[0]; - f[atom1][1] +=-F12[1]-F31[1]; - f[atom1][2] +=-F12[2]-F31[2]; - f[atom2][0] += F23[0]+F12[0]+F24[0]; - f[atom2][1] += F23[1]+F12[1]+F24[1]; - f[atom2][2] += F23[2]+F12[2]+F24[2]; - f[atom3][0] += -F23[0]+F34[0]+F31[0]; - f[atom3][1] += -F23[1]+F34[1]+F31[1]; - f[atom3][2] += -F23[2]+F34[2]+F31[2]; - f[atom4][0] += -F34[0]-F24[0]; - f[atom4][1] += -F34[1]-F24[1]; - f[atom4][2] += -F34[2]-F24[2]; + f1[0] = -F12[0]-F31[0]; + f1[1] = -F12[1]-F31[1]; + f1[2] = -F12[2]-F31[2]; + f2[0] = F23[0]+F12[0]+F24[0]; + f2[1] = F23[1]+F12[1]+F24[1]; + f2[2] = F23[2]+F12[2]+F24[2]; + f3[0] = -F23[0]+F34[0]+F31[0]; + f3[1] = -F23[1]+F34[1]+F31[1]; + f3[2] = -F23[2]+F34[2]+F31[2]; + f4[0] = -F34[0]-F24[0]; + f4[1] = -F34[1]-F24[1]; + f4[2] = -F34[2]-F24[2]; // coordination forces tmp2 = VA*Tij*((1.0-(om1234*om1234))) * (1.0-tspjik)*(1.0-tspijl)*dw21*w34/r21mag; - f[atom2][0] -= tmp2*r21[0]; - f[atom2][1] -= tmp2*r21[1]; - f[atom2][2] -= tmp2*r21[2]; - f[atom1][0] += tmp2*r21[0]; - f[atom1][1] += tmp2*r21[1]; - f[atom1][2] += tmp2*r21[2]; + f2[0] -= tmp2*r21[0]; + f2[1] -= tmp2*r21[1]; + f2[2] -= tmp2*r21[2]; + f1[0] += tmp2*r21[0]; + f1[1] += tmp2*r21[1]; + f1[2] += tmp2*r21[2]; tmp2 = VA*Tij*((1.0-(om1234*om1234))) * (1.0-tspjik)*(1.0-tspijl)*w21*dw34/r34mag; - f[atom3][0] -= tmp2*r34[0]; - f[atom3][1] -= tmp2*r34[1]; - f[atom3][2] -= tmp2*r34[2]; - f[atom4][0] += tmp2*r34[0]; - f[atom4][1] += tmp2*r34[1]; - f[atom4][2] += tmp2*r34[2]; + f3[0] -= tmp2*r34[0]; + f3[1] -= tmp2*r34[1]; + f3[2] -= tmp2*r34[2]; + f4[0] += tmp2*r34[0]; + f4[1] += tmp2*r34[1]; + f4[2] += tmp2*r34[2]; + + f[atom1][0] += f1[0]; f[atom1][1] += f1[1]; + f[atom1][2] += f1[2]; + f[atom2][0] += f2[0]; f[atom2][1] += f2[1]; + f[atom2][2] += f2[2]; + f[atom3][0] += f3[0]; f[atom3][1] += f3[1]; + f[atom3][2] += f3[2]; + f[atom4][0] += f4[0]; f[atom4][1] += f4[1]; + f[atom4][2] += f4[2]; + + if (vflag_atom) { + r13[0] = -rjk[0]; r13[1] = -rjk[1]; r13[2] = -rjk[2]; + r43[0] = -r34[0]; r43[1] = -r34[1]; r43[2] = -r34[2]; + v_tally4(atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43); + } } } } @@ -1872,6 +1941,8 @@ f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; + if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + tmp2 = VA*dN3[2]*(2.0*NconjtmpI*dwik*SpN)*Etmp/rikmag; f[atomi][0] -= tmp2*rik[0]; f[atomi][1] -= tmp2*rik[1]; @@ -1880,6 +1951,8 @@ f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; + if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + if (fabs(dNki) > TOL) { REBO_neighs_k = REBO_firstneigh[atomk]; for (n = 0; n < REBO_numneigh[atomk]; n++) { @@ -1899,6 +1972,8 @@ f[atomn][0] += tmp2*rkn[0]; f[atomn][1] += tmp2*rkn[1]; f[atomn][2] += tmp2*rkn[2]; + + if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn); } } } @@ -1929,6 +2004,8 @@ f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; + if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + tmp2 = VA*dN3[2]*(2.0*NconjtmpJ*dwjl*SpN)*Etmp/rjlmag; f[atomj][0] -= tmp2*rjl[0]; f[atomj][1] -= tmp2*rjl[1]; @@ -1937,6 +2014,8 @@ f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; + if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + if (fabs(dNlj) > TOL) { REBO_neighs_l = REBO_firstneigh[atoml]; for (n = 0; n < REBO_numneigh[atoml]; n++) { @@ -1956,6 +2035,8 @@ f[atomn][0] += tmp2*rln[0]; f[atomn][1] += tmp2*rln[1]; f[atomn][2] += tmp2*rln[2]; + + if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln); } } } @@ -1973,7 +2054,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, double VA, double rij0[3], double rij0mag, - double **f) + double **f, int vflag_atom) { int k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4; int atomi,atomj,itype,jtype,ktype,ltype,ntype; @@ -1987,11 +2068,9 @@ double dcosijldrj[3],dcosijldrl[3],dcosjikdrj[3],dwjl; double Tij,crosskij[3],crosskijmag; double crossijl[3],crossijlmag,omkijl; - double tmppij,tmppji,dN2PIJ[2],dN2PJI[2],dN3piRC[3],dN3Tij[3]; double bij,tmp3pij,tmp3pji,Stb,dStb; double r32[3],r32mag,cos321; - double om1234,rln[3]; double rlnmag,dwln,r23[3],r23mag,r21[3],r21mag; double w21,dw21,r34[3],r34mag,cos234,w34,dw34; @@ -2006,6 +2085,8 @@ double rij2,tspjik,dtsjik,tspijl,dtsijl,costmp; int *REBO_neighs,*REBO_neighs_i,*REBO_neighs_j,*REBO_neighs_k,*REBO_neighs_l; double F12[3],F23[3],F34[3],F31[3],F24[3]; + double fi[3],fj[3],fk[3],fl[3],f1[3],f2[3],f3[3],f4[4]; + double rji[3],rki[3],rlj[3],r13[3],r43[3]; double **x = atom->x; int *type = atom->type; @@ -2251,58 +2332,68 @@ g = gSpline(cosjik,(NijC+NijH),itype,&dgdc,&dgdN); tmp2 = VA*.5*(tmp*wik*dgdc*exp(lamdajik)); - f[atomj][0] -= tmp2*dcosjikdrj[0]; - f[atomj][1] -= tmp2*dcosjikdrj[1]; - f[atomj][2] -= tmp2*dcosjikdrj[2]; - f[atomi][0] -= tmp2*dcosjikdri[0]; - f[atomi][1] -= tmp2*dcosjikdri[1]; - f[atomi][2] -= tmp2*dcosjikdri[2]; - f[atomk][0] -= tmp2*dcosjikdrk[0]; - f[atomk][1] -= tmp2*dcosjikdrk[1]; - f[atomk][2] -= tmp2*dcosjikdrk[2]; + fj[0] = -tmp2*dcosjikdrj[0]; + fj[1] = -tmp2*dcosjikdrj[1]; + fj[2] = -tmp2*dcosjikdrj[2]; + fi[0] = -tmp2*dcosjikdri[0]; + fi[1] = -tmp2*dcosjikdri[1]; + fi[2] = -tmp2*dcosjikdri[2]; + fk[0] = -tmp2*dcosjikdrk[0]; + fk[1] = -tmp2*dcosjikdrk[1]; + fk[2] = -tmp2*dcosjikdrk[2]; tmp2 = VA*.5*(tmp*wik*g*exp(lamdajik)*4.0*kronecker(itype,1)); - f[atomj][0] -= tmp2*(-rij[0]/rijmag); - f[atomj][1] -= tmp2*(-rij[1]/rijmag); - f[atomj][2] -= tmp2*(-rij[2]/rijmag); - f[atomi][0] -= tmp2*((-rik[0]/rikmag)+(rij[0]/rijmag)); - f[atomi][1] -= tmp2*((-rik[1]/rikmag)+(rij[1]/rijmag)); - f[atomi][2] -= tmp2*((-rik[2]/rikmag)+(rij[2]/rijmag)); - f[atomk][0] -= tmp2*(rik[0]/rikmag); - f[atomk][1] -= tmp2*(rik[1]/rikmag); - f[atomk][2] -= tmp2*(rik[2]/rikmag); + fj[0] -= tmp2*(-rij[0]/rijmag); + fj[1] -= tmp2*(-rij[1]/rijmag); + fj[2] -= tmp2*(-rij[2]/rijmag); + fi[0] -= tmp2*((-rik[0]/rikmag)+(rij[0]/rijmag)); + fi[1] -= tmp2*((-rik[1]/rikmag)+(rij[1]/rijmag)); + fi[2] -= tmp2*((-rik[2]/rikmag)+(rij[2]/rijmag)); + fk[0] -= tmp2*(rik[0]/rikmag); + fk[1] -= tmp2*(rik[1]/rikmag); + fk[2] -= tmp2*(rik[2]/rikmag); // coordination forces // dwik forces tmp2 = VA*.5*(tmp*dwik*g*exp(lamdajik))/rikmag; - f[atomi][0] -= tmp2*rik[0]; - f[atomi][1] -= tmp2*rik[1]; - f[atomi][2] -= tmp2*rik[2]; - f[atomk][0] += tmp2*rik[0]; - f[atomk][1] += tmp2*rik[1]; - f[atomk][2] += tmp2*rik[2]; + fi[0] -= tmp2*rik[0]; + fi[1] -= tmp2*rik[1]; + fi[2] -= tmp2*rik[2]; + fk[0] += tmp2*rik[0]; + fk[1] += tmp2*rik[1]; + fk[2] += tmp2*rik[2]; // PIJ forces tmp2 = VA*.5*(tmp*dN2[ktype]*dwik)/rikmag; - f[atomi][0] -= tmp2*rik[0]; - f[atomi][1] -= tmp2*rik[1]; - f[atomi][2] -= tmp2*rik[2]; - f[atomk][0] += tmp2*rik[0]; - f[atomk][1] += tmp2*rik[1]; - f[atomk][2] += tmp2*rik[2]; + fi[0] -= tmp2*rik[0]; + fi[1] -= tmp2*rik[1]; + fi[2] -= tmp2*rik[2]; + fk[0] += tmp2*rik[0]; + fk[1] += tmp2*rik[1]; + fk[2] += tmp2*rik[2]; // dgdN forces tmp2 = VA*.5*(tmp*tmp3*dwik)/rikmag; - f[atomi][0] -= tmp2*rik[0]; - f[atomi][1] -= tmp2*rik[1]; - f[atomi][2] -= tmp2*rik[2]; - f[atomk][0] += tmp2*rik[0]; - f[atomk][1] += tmp2*rik[1]; - f[atomk][2] += tmp2*rik[2]; + fi[0] -= tmp2*rik[0]; + fi[1] -= tmp2*rik[1]; + fi[2] -= tmp2*rik[2]; + fk[0] += tmp2*rik[0]; + fk[1] += tmp2*rik[1]; + fk[2] += tmp2*rik[2]; + + f[atomi][0] += fi[0]; f[atomi][1] += fi[1]; f[atomi][2] += fi[2]; + f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; + f[atomk][0] += fk[0]; f[atomk][1] += fk[1]; f[atomk][2] += fk[2]; + + if (vflag_atom) { + rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2]; + rki[0] = -rik[0]; rki[1] = -rik[1]; rki[2] = -rik[2]; + v_tally3(atomi,atomj,atomk,fj,fk,rji,rki); + } } } @@ -2350,57 +2441,66 @@ g = gSpline(cosijl,NjiC+NjiH,jtype,&dgdc,&dgdN); tmp2 = VA*.5*(tmp*wjl*dgdc*exp(lamdaijl)); - f[atomi][0] -= tmp2*dcosijldri[0]; - f[atomi][1] -= tmp2*dcosijldri[1]; - f[atomi][2] -= tmp2*dcosijldri[2]; - f[atomj][0] -= tmp2*dcosijldrj[0]; - f[atomj][1] -= tmp2*dcosijldrj[1]; - f[atomj][2] -= tmp2*dcosijldrj[2]; - f[atoml][0] -= tmp2*dcosijldrl[0]; - f[atoml][1] -= tmp2*dcosijldrl[1]; - f[atoml][2] -= tmp2*dcosijldrl[2]; + fi[0] = -tmp2*dcosijldri[0]; + fi[1] = -tmp2*dcosijldri[1]; + fi[2] = -tmp2*dcosijldri[2]; + fj[0] = -tmp2*dcosijldrj[0]; + fj[1] = -tmp2*dcosijldrj[1]; + fj[2] = -tmp2*dcosijldrj[2]; + fl[0] = -tmp2*dcosijldrl[0]; + fl[1] = -tmp2*dcosijldrl[1]; + fl[2] = -tmp2*dcosijldrl[2]; tmp2 = VA*.5*(tmp*wjl*g*exp(lamdaijl)*4.0*kronecker(jtype,1)); - f[atomi][0] -= tmp2*(rij[0]/rijmag); - f[atomi][1] -= tmp2*(rij[1]/rijmag); - f[atomi][2] -= tmp2*(rij[2]/rijmag); - f[atomj][0] -= tmp2*((-rjl[0]/rjlmag)-(rij[0]/rijmag)); - f[atomj][1] -= tmp2*((-rjl[1]/rjlmag)-(rij[1]/rijmag)); - f[atomj][2] -= tmp2*((-rjl[2]/rjlmag)-(rij[2]/rijmag)); - f[atoml][0] -= tmp2*(rjl[0]/rjlmag); - f[atoml][1] -= tmp2*(rjl[1]/rjlmag); - f[atoml][2] -= tmp2*(rjl[2]/rjlmag); + fi[0] -= tmp2*(rij[0]/rijmag); + fi[1] -= tmp2*(rij[1]/rijmag); + fi[2] -= tmp2*(rij[2]/rijmag); + fj[0] -= tmp2*((-rjl[0]/rjlmag)-(rij[0]/rijmag)); + fj[1] -= tmp2*((-rjl[1]/rjlmag)-(rij[1]/rijmag)); + fj[2] -= tmp2*((-rjl[2]/rjlmag)-(rij[2]/rijmag)); + fl[0] -= tmp2*(rjl[0]/rjlmag); + fl[1] -= tmp2*(rjl[1]/rjlmag); + fl[2] -= tmp2*(rjl[2]/rjlmag); // coordination forces // dwik forces tmp2 = VA*.5*(tmp*dwjl*g*exp(lamdaijl))/rjlmag; - f[atomj][0] -= tmp2*rjl[0]; - f[atomj][1] -= tmp2*rjl[1]; - f[atomj][2] -= tmp2*rjl[2]; - f[atoml][0] += tmp2*rjl[0]; - f[atoml][1] += tmp2*rjl[1]; - f[atoml][2] += tmp2*rjl[2]; + fj[0] -= tmp2*rjl[0]; + fj[1] -= tmp2*rjl[1]; + fj[2] -= tmp2*rjl[2]; + fl[0] += tmp2*rjl[0]; + fl[1] += tmp2*rjl[1]; + fl[2] += tmp2*rjl[2]; // PIJ forces tmp2 = VA*.5*(tmp*dN2[ltype]*dwjl)/rjlmag; - f[atomj][0] -= tmp2*rjl[0]; - f[atomj][1] -= tmp2*rjl[1]; - f[atomj][2] -= tmp2*rjl[2]; - f[atoml][0] += tmp2*rjl[0]; - f[atoml][1] += tmp2*rjl[1]; - f[atoml][2] += tmp2*rjl[2]; + fj[0] -= tmp2*rjl[0]; + fj[1] -= tmp2*rjl[1]; + fj[2] -= tmp2*rjl[2]; + fl[0] += tmp2*rjl[0]; + fl[1] += tmp2*rjl[1]; + fl[2] += tmp2*rjl[2]; // dgdN forces tmp2=VA*.5*(tmp*tmp3*dwjl)/rjlmag; - f[atomj][0] -= tmp2*rjl[0]; - f[atomj][1] -= tmp2*rjl[1]; - f[atomj][2] -= tmp2*rjl[2]; - f[atoml][0] += tmp2*rjl[0]; - f[atoml][1] += tmp2*rjl[1]; - f[atoml][2] += tmp2*rjl[2]; + fj[0] -= tmp2*rjl[0]; + fj[1] -= tmp2*rjl[1]; + fj[2] -= tmp2*rjl[2]; + fl[0] += tmp2*rjl[0]; + fl[1] += tmp2*rjl[1]; + fl[2] += tmp2*rjl[2]; + + f[atomi][0] += fi[0]; f[atomi][1] += fi[1]; f[atomi][2] += fi[2]; + f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; + f[atoml][0] += fl[0]; f[atoml][1] += fl[1]; f[atoml][2] += fl[2]; + + if (vflag_atom) { + rlj[0] = -rjl[0]; rlj[1] = -rjl[1]; rlj[2] = -rjl[2]; + v_tally3(atomi,atomj,atoml,fi,fl,rij,rlj); + } } } @@ -2432,6 +2532,8 @@ f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; + if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + tmp2 = VA*dN3[2]*(2.0*NconjtmpI*dwik*SpN)/rikmag; f[atomi][0] -= tmp2*rik[0]; f[atomi][1] -= tmp2*rik[1]; @@ -2440,6 +2542,8 @@ f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; + if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + if (fabs(dNki) > TOL) { REBO_neighs_k = REBO_firstneigh[atomk]; for (n = 0; n < REBO_numneigh[atomk]; n++) { @@ -2459,6 +2563,8 @@ f[atomn][0] += tmp2*rkn[0]; f[atomn][1] += tmp2*rkn[1]; f[atomn][2] += tmp2*rkn[2]; + + if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn); } } } @@ -2489,6 +2595,8 @@ f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; + if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + tmp2 = VA*dN3[2]*(2.0*NconjtmpJ*dwjl*SpN)/rjlmag; f[atomj][0] -= tmp2*rjl[0]; f[atomj][1] -= tmp2*rjl[1]; @@ -2497,6 +2605,8 @@ f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; + if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + if (fabs(dNlj) > TOL) { REBO_neighs_l = REBO_firstneigh[atoml]; for (n = 0; n < REBO_numneigh[atoml]; n++) { @@ -2516,6 +2626,8 @@ f[atomn][0] += tmp2*rln[0]; f[atomn][1] += tmp2*rln[1]; f[atomn][2] += tmp2*rln[2]; + + if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln); } } } @@ -2686,38 +2798,53 @@ F24[1] = (fcilpc*ril[1]); F24[2] = (fcilpc*ril[2]); - f[atom1][0] +=-F12[0]-F31[0]; - f[atom1][1] +=-F12[1]-F31[1]; - f[atom1][2] +=-F12[2]-F31[2]; - f[atom2][0] += F23[0]+F12[0]+F24[0]; - f[atom2][1] += F23[1]+F12[1]+F24[1]; - f[atom2][2] += F23[2]+F12[2]+F24[2]; - f[atom3][0] += -F23[0]+F34[0]+F31[0]; - f[atom3][1] += -F23[1]+F34[1]+F31[1]; - f[atom3][2] += -F23[2]+F34[2]+F31[2]; - f[atom4][0] += -F34[0]-F24[0]; - f[atom4][1] += -F34[1]-F24[1]; - f[atom4][2] += -F34[2]-F24[2]; + f1[0] = -F12[0]-F31[0]; + f1[1] = -F12[1]-F31[1]; + f1[2] = -F12[2]-F31[2]; + f2[0] = F23[0]+F12[0]+F24[0]; + f2[1] = F23[1]+F12[1]+F24[1]; + f2[2] = F23[2]+F12[2]+F24[2]; + f3[0] = -F23[0]+F34[0]+F31[0]; + f3[1] = -F23[1]+F34[1]+F31[1]; + f3[2] = -F23[2]+F34[2]+F31[2]; + f4[0] = -F34[0]-F24[0]; + f4[1] = -F34[1]-F24[1]; + f4[2] = -F34[2]-F24[2]; // coordination forces tmp2 = VA*Tij*((1.0-(om1234*om1234))) * (1.0-tspjik)*(1.0-tspijl)*dw21*w34/r21mag; - f[atom2][0] -= tmp2*r21[0]; - f[atom2][1] -= tmp2*r21[1]; - f[atom2][2] -= tmp2*r21[2]; - f[atom1][0] += tmp2*r21[0]; - f[atom1][1] += tmp2*r21[1]; - f[atom1][2] += tmp2*r21[2]; + f2[0] -= tmp2*r21[0]; + f2[1] -= tmp2*r21[1]; + f2[2] -= tmp2*r21[2]; + f1[0] += tmp2*r21[0]; + f1[1] += tmp2*r21[1]; + f1[2] += tmp2*r21[2]; tmp2 = VA*Tij*((1.0-(om1234*om1234))) * (1.0-tspjik)*(1.0-tspijl)*w21*dw34/r34mag; - f[atom3][0] -= tmp2*r34[0]; - f[atom3][1] -= tmp2*r34[1]; - f[atom3][2] -= tmp2*r34[2]; - f[atom4][0] += tmp2*r34[0]; - f[atom4][1] += tmp2*r34[1]; - f[atom4][2] += tmp2*r34[2]; + f3[0] -= tmp2*r34[0]; + f3[1] -= tmp2*r34[1]; + f3[2] -= tmp2*r34[2]; + f4[0] += tmp2*r34[0]; + f4[1] += tmp2*r34[1]; + f4[2] += tmp2*r34[2]; + + f[atom1][0] += f1[0]; f[atom1][1] += f1[1]; + f[atom1][2] += f1[2]; + f[atom2][0] += f2[0]; f[atom2][1] += f2[1]; + f[atom2][2] += f2[2]; + f[atom3][0] += f3[0]; f[atom3][1] += f3[1]; + f[atom3][2] += f3[2]; + f[atom4][0] += f4[0]; f[atom4][1] += f4[1]; + f[atom4][2] += f4[2]; + + if (vflag_atom) { + r13[0] = -rjk[0]; r13[1] = -rjk[1]; r13[2] = -rjk[2]; + r43[0] = -r34[0]; r43[1] = -r34[1]; r43[2] = -r34[2]; + v_tally4(atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43); + } } } } @@ -2747,6 +2874,8 @@ f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; + if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + tmp2 = VA*dN3[2]*(2.0*NconjtmpI*dwik*SpN)*Etmp/rikmag; f[atomi][0] -= tmp2*rik[0]; f[atomi][1] -= tmp2*rik[1]; @@ -2755,6 +2884,8 @@ f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; + if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + if (fabs(dNki) >TOL) { REBO_neighs_k = REBO_firstneigh[atomk]; for (n = 0; n < REBO_numneigh[atomk]; n++) { @@ -2774,6 +2905,8 @@ f[atomn][0] += tmp2*rkn[0]; f[atomn][1] += tmp2*rkn[1]; f[atomn][2] += tmp2*rkn[2]; + + if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn); } } } @@ -2804,6 +2937,8 @@ f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; + if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + tmp2 = VA*dN3[2]*(2.0*NconjtmpJ*dwjl*SpN)*Etmp/rjlmag; f[atomj][0] -= tmp2*rjl[0]; f[atomj][1] -= tmp2*rjl[1]; @@ -2812,6 +2947,8 @@ f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; + if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + if (fabs(dNlj) > TOL) { REBO_neighs_l = REBO_firstneigh[atoml]; for (n = 0; n < REBO_numneigh[atoml]; n++) { @@ -2831,6 +2968,8 @@ f[atomn][0] += tmp2*rln[0]; f[atomn][1] += tmp2*rln[1]; f[atomn][2] += tmp2*rln[2]; + + if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln); } } } diff -Naur lammps-16Jan08/src/MANYBODY/pair_airebo.h lammps-17Jan08/src/MANYBODY/pair_airebo.h --- lammps-16Jan08/src/MANYBODY/pair_airebo.h 2008-01-08 16:34:19.000000000 -0700 +++ lammps-17Jan08/src/MANYBODY/pair_airebo.h 2008-01-17 10:01:57.000000000 -0700 @@ -82,9 +82,9 @@ void FLJ(int, int); void TORSION(int, int); - double bondorder(int, int, double *, double, double, double **); + double bondorder(int, int, double *, double, double, double **, int); double bondorderLJ(int, int, double *, double, double, - double *, double, double **); + double *, double, double **, int); double Sp(double, double, double, double &); double Sp2(double, double, double, double &); diff -Naur lammps-16Jan08/src/MANYBODY/pair_tersoff.cpp lammps-17Jan08/src/MANYBODY/pair_tersoff.cpp --- lammps-16Jan08/src/MANYBODY/pair_tersoff.cpp 2007-11-30 18:01:35.000000000 -0700 +++ lammps-17Jan08/src/MANYBODY/pair_tersoff.cpp 2008-01-17 10:01:57.000000000 -0700 @@ -226,7 +226,7 @@ f[k][1] += fk[1]; f[k][2] += fk[2]; - if (evflag) ev_tally3(i,j,k,0.0,0.0,fj,fk,delr1,delr2); + if (vflag_atom) v_tally3(i,j,k,fj,fk,delr1,delr2); } } } diff -Naur lammps-16Jan08/src/MEAM/pair_meam.cpp lammps-17Jan08/src/MEAM/pair_meam.cpp --- lammps-16Jan08/src/MEAM/pair_meam.cpp 2008-01-08 16:34:19.000000000 -0700 +++ lammps-17Jan08/src/MEAM/pair_meam.cpp 2008-01-17 10:01:57.000000000 -0700 @@ -54,7 +54,6 @@ rho = rho0 = rho1 = rho2 = rho3 = frhop = NULL; gamma = dgamma1 = dgamma2 = dgamma3 = arho2b = NULL; arho1 = arho2 = arho3 = arho3b = t_ave = NULL; - strssa = NULL; maxneigh = 0; scrfcn = dscrfcn = fcpair = NULL; @@ -96,8 +95,6 @@ memory->destroy_2d_double_array(arho3b); memory->destroy_2d_double_array(t_ave); - memory->destroy_3d_double_array(strssa); - memory->sfree(scrfcn); memory->sfree(dscrfcn); memory->sfree(fcpair); @@ -127,11 +124,6 @@ if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; - // error check (for now) - - if (eflag_atom || vflag_atom) - error->all("Pair style meam does not yet support peratom energy/virial"); - int newton_pair = force->newton_pair; // grow local arrays if necessary @@ -153,7 +145,6 @@ memory->destroy_2d_double_array(arho3); memory->destroy_2d_double_array(arho3b); memory->destroy_2d_double_array(t_ave); - memory->destroy_3d_double_array(strssa); nmax = atom->nmax; @@ -173,7 +164,6 @@ arho3 = memory->create_2d_double_array(nmax,10,"pair:arho3"); arho3b = memory->create_2d_double_array(nmax,3,"pair:arho3b"); t_ave = memory->create_2d_double_array(nmax,3,"pair:t_ave"); - strssa = memory->create_3d_double_array(nmax,3,3,"pair:strssa"); } // neighbor list info @@ -238,7 +228,7 @@ for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; ifort = i+1; - meam_dens_init_(&ifort,&nmax,&eflag,&evdwl,&ntype,type,fmap,&x[0][0], + meam_dens_init_(&ifort,&nmax,&ntype,type,fmap,&x[0][0], &numneigh_half[i],firstneigh_half[i], &numneigh_full[i],firstneigh_full[i], &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], @@ -252,10 +242,10 @@ offset += numneigh_half[i]; } - reverse_flag = 0; comm->reverse_comm_pair(this); - meam_dens_final_(&nlocal,&nmax,&eflag,&evdwl,&ntype,type,fmap, + meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom, + &eng_vdwl,eatom,&ntype,type,fmap, &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0], &arho3b[0][0],&t_ave[0][0],gamma,dgamma1, dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag); @@ -269,16 +259,24 @@ offset = 0; + // vptr is first value in vatom if it will be used by meam_force() + // else vatom may not exist, so pass dummy ptr + + double *vptr; + if (vflag_atom) vptr = &vatom[0][0]; + else vptr = &cutmax; + for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; ifort = i+1; - meam_force_(&ifort,&nmax,&eflag,&evdwl,&ntype,type,fmap,&x[0][0], + meam_force_(&ifort,&nmax,&eflag_either,&eflag_global,&eflag_atom, + &vflag_atom,&eng_vdwl,eatom,&ntype,type,fmap,&x[0][0], &numneigh_half[i],firstneigh_half[i], &numneigh_full[i],firstneigh_full[i], &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop, &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0], - &t_ave[0][0],&f[0][0],&strssa[0][0][0],&errorflag); + &t_ave[0][0],&f[0][0],vptr,&errorflag); if (errorflag) { char str[128]; sprintf(str,"MEAM library error %d",errorflag); @@ -287,17 +285,11 @@ offset += numneigh_half[i]; } - reverse_flag = 1; - comm->reverse_comm_pair(this); - // change neighbor list indices back to C indexing neigh_f2c(inum_half,ilist_half,numneigh_half,firstneigh_half); neigh_f2c(inum_half,ilist_half,numneigh_full,firstneigh_full); - // just sum global energy (for now) - - if (evflag) ev_tally(0,0,nlocal,newton_pair,evdwl,0.0,0.0,0.0,0.0,0.0); if (vflag_fdotr) virial_compute(); } @@ -811,45 +803,28 @@ m = 0; last = first + n; - if (reverse_flag == 0) { - for (i = first; i < last; i++) { - buf[m++] = rho0[i]; - buf[m++] = arho2b[i]; - buf[m++] = arho1[i][0]; - buf[m++] = arho1[i][1]; - buf[m++] = arho1[i][2]; - buf[m++] = arho2[i][0]; - buf[m++] = arho2[i][1]; - buf[m++] = arho2[i][2]; - buf[m++] = arho2[i][3]; - buf[m++] = arho2[i][4]; - buf[m++] = arho2[i][5]; - for (k = 0; k < 10; k++) buf[m++] = arho3[i][k]; - buf[m++] = arho3b[i][0]; - buf[m++] = arho3b[i][1]; - buf[m++] = arho3b[i][2]; - buf[m++] = t_ave[i][0]; - buf[m++] = t_ave[i][1]; - buf[m++] = t_ave[i][2]; - } - size = 27; - - } else { - for (i = first; i < last; i++) { - buf[m++] = strssa[i][0][0]; - buf[m++] = strssa[i][0][1]; - buf[m++] = strssa[i][0][2]; - buf[m++] = strssa[i][1][0]; - buf[m++] = strssa[i][1][1]; - buf[m++] = strssa[i][1][2]; - buf[m++] = strssa[i][2][0]; - buf[m++] = strssa[i][2][1]; - buf[m++] = strssa[i][2][2]; - } - size = 9; + for (i = first; i < last; i++) { + buf[m++] = rho0[i]; + buf[m++] = arho2b[i]; + buf[m++] = arho1[i][0]; + buf[m++] = arho1[i][1]; + buf[m++] = arho1[i][2]; + buf[m++] = arho2[i][0]; + buf[m++] = arho2[i][1]; + buf[m++] = arho2[i][2]; + buf[m++] = arho2[i][3]; + buf[m++] = arho2[i][4]; + buf[m++] = arho2[i][5]; + for (k = 0; k < 10; k++) buf[m++] = arho3[i][k]; + buf[m++] = arho3b[i][0]; + buf[m++] = arho3b[i][1]; + buf[m++] = arho3b[i][2]; + buf[m++] = t_ave[i][0]; + buf[m++] = t_ave[i][1]; + buf[m++] = t_ave[i][2]; } - return size; + return m; } /* ---------------------------------------------------------------------- */ @@ -859,42 +834,26 @@ int i,j,k,m; m = 0; - if (reverse_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - rho0[j] += buf[m++]; - arho2b[j] += buf[m++]; - arho1[j][0] += buf[m++]; - arho1[j][1] += buf[m++]; - arho1[j][2] += buf[m++]; - arho2[j][0] += buf[m++]; - arho2[j][1] += buf[m++]; - arho2[j][2] += buf[m++]; - arho2[j][3] += buf[m++]; - arho2[j][4] += buf[m++]; - arho2[j][5] += buf[m++]; - for (k = 0; k < 10; k++) arho3[j][k] += buf[m++]; - arho3b[j][0] += buf[m++]; - arho3b[j][1] += buf[m++]; - arho3b[j][2] += buf[m++]; - t_ave[j][0] += buf[m++]; - t_ave[j][1] += buf[m++]; - t_ave[j][2] += buf[m++]; - } - - } else { - for (i = 0; i < n; i++) { - j = list[i]; - strssa[j][0][0] += buf[m++]; - strssa[j][0][1] += buf[m++]; - strssa[j][0][2] += buf[m++]; - strssa[j][1][0] += buf[m++]; - strssa[j][1][1] += buf[m++]; - strssa[j][1][2] += buf[m++]; - strssa[j][2][0] += buf[m++]; - strssa[j][2][1] += buf[m++]; - strssa[j][2][2] += buf[m++]; - } + for (i = 0; i < n; i++) { + j = list[i]; + rho0[j] += buf[m++]; + arho2b[j] += buf[m++]; + arho1[j][0] += buf[m++]; + arho1[j][1] += buf[m++]; + arho1[j][2] += buf[m++]; + arho2[j][0] += buf[m++]; + arho2[j][1] += buf[m++]; + arho2[j][2] += buf[m++]; + arho2[j][3] += buf[m++]; + arho2[j][4] += buf[m++]; + arho2[j][5] += buf[m++]; + for (k = 0; k < 10; k++) arho3[j][k] += buf[m++]; + arho3b[j][0] += buf[m++]; + arho3b[j][1] += buf[m++]; + arho3b[j][2] += buf[m++]; + t_ave[j][0] += buf[m++]; + t_ave[j][1] += buf[m++]; + t_ave[j][2] += buf[m++]; } } @@ -906,7 +865,6 @@ { double bytes = 11 * nmax * sizeof(double); bytes += (3 + 6 + 10 + 3 + 3) * nmax * sizeof(double); - bytes += 3*3 * nmax * sizeof(double); bytes += 3 * maxneigh * sizeof(double); return bytes; } diff -Naur lammps-16Jan08/src/MEAM/pair_meam.h lammps-17Jan08/src/MEAM/pair_meam.h --- lammps-16Jan08/src/MEAM/pair_meam.h 2007-10-04 18:32:35.000000000 -0600 +++ lammps-17Jan08/src/MEAM/pair_meam.h 2008-01-17 10:01:57.000000000 -0700 @@ -22,20 +22,22 @@ void meam_setup_param_(int *, double *, int *, int *, int *); void meam_setup_done_(double *); - void meam_dens_init_(int *, int *, int *, double *, int *, int *, int *, + void meam_dens_init_(int *, int *, int *, int *, int *, double *, int *, int *, int *, int *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, int *); - void meam_dens_final_(int *, int *, int *, double *, int *, int *, int *, + void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *, + int *, int *, int *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, int *); - void meam_force_(int *, int *, int *, double *, int *, int *, int *, + void meam_force_(int *, int *, int *, int *, int *, int *, + double *, double *, int *, int *, int *, double *, int *, int *, int *, int *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, @@ -71,7 +73,6 @@ int nelements; // # of unique elements char **elements; // names of unique elements double *mass; // mass of each element - int reverse_flag; // which pass of reverse comm is being done int *map; // mapping from atom types to elements int *fmap; // Fortran version of map array for MEAM lib @@ -83,7 +84,6 @@ double *rho,*rho0,*rho1,*rho2,*rho3,*frhop; double *gamma,*dgamma1,*dgamma2,*dgamma3,*arho2b; double **arho1,**arho2,**arho3,**arho3b,**t_ave; - double ***strssa; void allocate(); void read_files(char *, char *); diff -Naur lammps-16Jan08/src/POEMS/fix_poems.cpp lammps-17Jan08/src/POEMS/fix_poems.cpp --- lammps-16Jan08/src/POEMS/fix_poems.cpp 2007-10-10 16:22:16.000000000 -0600 +++ lammps-17Jan08/src/POEMS/fix_poems.cpp 2008-01-17 10:01:57.000000000 -0700 @@ -57,6 +57,7 @@ rigid_flag = 1; virial_flag = 1; + MPI_Comm_rank(world,&me); // perform initial allocation of atom-based arrays @@ -254,10 +255,6 @@ fprintf(logfile,"%d clusters, %d bodies, %d joints, %d atoms\n", ncluster,nbody,njoint,nsum); } - - // zero fix_poems virial in case pressure uses it before 1st fix_poems call - - for (int n = 0; n < 6; n++) virial[n] = 0.0; } /* ---------------------------------------------------------------------- @@ -346,14 +343,6 @@ error->all("POEMS fix must come before NPT/NPH fix"); } - // compute poems contribution to virial every step if fix NPT,NPH exists - - pressure_flag = 0; - for (int i = 0; i < modify->nfix; i++) { - if (strcmp(modify->fix[i]->style,"npt") == 0) pressure_flag = 1; - if (strcmp(modify->fix[i]->style,"nph") == 0) pressure_flag = 1; - } - // timestep info dtv = update->dt; @@ -597,9 +586,9 @@ make setup call to POEMS ------------------------------------------------------------------------- */ -void FixPOEMS::setup() +void FixPOEMS::setup(int vflag) { - int i,j,ibody; + int i,j,n,ibody; // vcm = velocity of center-of-mass of each rigid body // angmom = angular momentum of each rigid body @@ -654,29 +643,31 @@ angmom[ibody][2] = all[ibody][5]; } - // set omega from angmom + // virial setup before call to set_v + + if (vflag) v_setup(vflag); + else evflag = 0; + + // set velocities from angmom & omega for (ibody = 0; ibody < nbody; ibody++) omega_from_mq(angmom[ibody],ex_space[ibody],ey_space[ibody], ez_space[ibody],inertia[ibody],omega[ibody]); + set_v(); - // reset velocities from omega // guestimate virial as 2x the set_v contribution - // since post_force doesn't compute virial - - for (int n = 0; n < 6; n++) virial[n] = 0.0; - set_v(1); - for (int n = 0; n < 6; n++) virial[n] *= 2.0; - double ke = 0.0; - for (i = 0; i < nlocal; i++) - if (natom2body[i]) - ke += mass[atom->type[i]] * - (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); + if (vflag_global) + for (n = 0; n < 6; n++) virial[n] *= 2.0; + if (vflag_atom) { + for (i = 0; i < nlocal; i++) + for (n = 0; n < 6; n++) + vatom[i][n] *= 2.0; + } // use post_force() to compute initial fcm & torque - post_force(1); + post_force(vflag); // setup for POEMS @@ -691,17 +682,20 @@ set x,v of body atoms accordingly /* ---------------------------------------------------------------------- */ -void FixPOEMS::initial_integrate() +void FixPOEMS::initial_integrate(int vflag) { // perform POEMS integration poems->LobattoOne(xcm,vcm,omega,torque,fcm,ex_space,ey_space,ez_space); - // set coords and velocities of atoms in rigid bodies + // virial setup before call to set_xv + + if (vflag) v_setup(vflag); + else evflag = 0; - int vflag = 0; - if (pressure_flag || output->next_thermo == update->ntimestep) vflag = 1; - set_xv(vflag); + // set coords and velocities of atoms in rigid bodies + + set_xv(); } /* ---------------------------------------------------------------------- @@ -712,6 +706,8 @@ void FixPOEMS::post_force(int vflag) { int i,j,ibody; + int xbox,ybox,zbox; + double dx,dy,dz; int *image = atom->image; double **x = atom->x; @@ -722,8 +718,6 @@ double yprd = domain->yprd; double zprd = domain->zprd; - int xbox,ybox,zbox; - double dx,dy,dz; for (ibody = 0; ibody < nbody; ibody++) for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; @@ -772,15 +766,14 @@ poems->LobattoTwo(vcm,omega,torque,fcm); // set velocities of atoms in rigid bodies + // virial is already setup from initial_integrate - int vflag = 0; - if (pressure_flag || output->next_thermo == update->ntimestep) vflag = 1; - set_v(vflag); + set_v(); } /* ---------------------------------------------------------------------- */ -void FixPOEMS::initial_integrate_respa(int ilevel, int flag) +void FixPOEMS::initial_integrate_respa(int vflag, int ilevel, int flag) { if (flag) return; // only used by NPT,NPH @@ -788,7 +781,7 @@ dtf = 0.5 * step_respa[ilevel] * force->ftm2v; dthalf = 0.5 * step_respa[ilevel]; - if (ilevel == 0) initial_integrate(); + if (ilevel == 0) initial_integrate(vflag); else final_integrate(); } @@ -1344,30 +1337,26 @@ v = Vcm + (W cross (x - Xcm)) ------------------------------------------------------------------------- */ -void FixPOEMS::set_xv(int vflag) +void FixPOEMS::set_xv() { + int ibody; + int xbox,ybox,zbox; + double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone; + double vr[6]; + int *image = atom->image; double **x = atom->x; double **v = atom->v; + double **f = atom->f; + double *mass = atom->mass; + int *type = atom->type; int nlocal = atom->nlocal; double xprd = domain->xprd; double yprd = domain->yprd; double zprd = domain->zprd; - int ibody; - int xbox,ybox,zbox; - - double vold0,vold1,vold2,fc0,fc1,fc2,massone,x0,x1,x2; - double *mass = atom->mass; - double **f = atom->f; - int *type = atom->type; - - // zero out fix_poems virial - - if (vflag) for (int n = 0; n < 6; n++) virial[n] = 0.0; - - // set x and v + // set x and v of each atom // only set joint atoms for 1st rigid body they belong to for (int i = 0; i < nlocal; i++) { @@ -1378,18 +1367,21 @@ ybox = (image[i] >> 10 & 1023) - 512; zbox = (image[i] >> 20) - 512; - // save old positions and velocities for virial contribution + // save old positions and velocities for virial - if (vflag) { + if (evflag) { x0 = x[i][0] + xbox*xprd; x1 = x[i][1] + ybox*yprd; x2 = x[i][2] + zbox*zprd; - vold0 = v[i][0]; - vold1 = v[i][1]; - vold2 = v[i][2]; + v0 = v[i][0]; + v1 = v[i][1]; + v2 = v[i][2]; } + // x = displacement from center-of-mass, based on body orientation + // v = vcm + omega around center-of-mass + x[i][0] = ex_space[ibody][0]*displace[i][0] + ey_space[ibody][0]*displace[i][1] + ez_space[ibody][0]*displace[i][2]; @@ -1407,24 +1399,33 @@ v[i][2] = omega[ibody][0]*x[i][1] - omega[ibody][1]*x[i][0] + vcm[ibody][2]; + // add center of mass to displacement + // map back into periodic box via xbox,ybox,zbox + x[i][0] += xcm[ibody][0] - xbox*xprd; x[i][1] += xcm[ibody][1] - ybox*yprd; x[i][2] += xcm[ibody][2] - zbox*zprd; - // compute body constraint forces for virial + // virial = unwrapped coords dotted into body constraint force + // body constraint force = implied force due to v change minus f external + // assume f does not include forces internal to body + // 1/2 factor b/c final_integrate contributes other half + // assume per-atom contribution is due to constraint force on that atom - if (vflag) { + if (evflag) { massone = mass[type[i]]; - fc0 = massone*(v[i][0] - vold0)/dtf - f[i][0]; - fc1 = massone*(v[i][1] - vold1)/dtf - f[i][1]; - fc2 = massone*(v[i][2] - vold2)/dtf - f[i][2]; - - virial[0] += 0.5*fc0*x0; - virial[1] += 0.5*fc1*x1; - virial[2] += 0.5*fc2*x2; - virial[3] += 0.5*fc1*x0; - virial[4] += 0.5*fc2*x0; - virial[5] += 0.5*fc2*x1; + fc0 = massone*(v[i][0] - v0)/dtf - f[i][0]; + fc1 = massone*(v[i][1] - v1)/dtf - f[i][1]; + fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; + + vr[0] = 0.5*fc0*x0; + vr[1] = 0.5*fc1*x1; + vr[2] = 0.5*fc2*x2; + vr[3] = 0.5*fc1*x0; + vr[4] = 0.5*fc2*x0; + vr[5] = 0.5*fc2*x1; + + v_tally(1,&i,1.0,vr); } } } @@ -1434,26 +1435,27 @@ v = Vcm + (W cross (x - Xcm)) ------------------------------------------------------------------------- */ -void FixPOEMS::set_v(int vflag) +void FixPOEMS::set_v() { - double **v = atom->v; - int nlocal = atom->nlocal; - int ibody; + int xbox,ybox,zbox; double dx,dy,dz; + double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone; + double vr[6]; - double vold0,vold1,vold2,fc0,fc1,fc2,massone,x0,x1,x2; double *mass = atom->mass; double **f = atom->f; double **x = atom->x; + double **v = atom->v; int *type = atom->type; int *image = atom->image; + int nlocal = atom->nlocal; + double xprd = domain->xprd; double yprd = domain->yprd; double zprd = domain->zprd; - int xbox,ybox,zbox; - // set v + // set v of each atom // only set joint atoms for 1st rigid body they belong to for (int i = 0; i < nlocal; i++) { @@ -1472,24 +1474,27 @@ // save old velocities for virial - if (vflag) { - vold0 = v[i][0]; - vold1 = v[i][1]; - vold2 = v[i][2]; + if (evflag) { + v0 = v[i][0]; + v1 = v[i][1]; + v2 = v[i][2]; } v[i][0] = omega[ibody][1]*dz - omega[ibody][2]*dy + vcm[ibody][0]; v[i][1] = omega[ibody][2]*dx - omega[ibody][0]*dz + vcm[ibody][1]; v[i][2] = omega[ibody][0]*dy - omega[ibody][1]*dx + vcm[ibody][2]; - // compute body constraint forces for virial - // use unwrapped atom positions + // virial = unwrapped coords dotted into body constraint force + // body constraint force = implied force due to v change minus f external + // assume f does not include forces internal to body + // 1/2 factor b/c initial_integrate contributes other half + // assume per-atom contribution is due to constraint force on that atom - if (vflag) { + if (evflag) { massone = mass[type[i]]; - fc0 = massone*(v[i][0] - vold0)/dtf - f[i][0]; - fc1 = massone*(v[i][1] - vold1)/dtf - f[i][1]; - fc2 = massone*(v[i][2] - vold2)/dtf - f[i][2]; + fc0 = massone*(v[i][0] - v0)/dtf - f[i][0]; + fc1 = massone*(v[i][1] - v1)/dtf - f[i][1]; + fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; xbox = (image[i] & 1023) - 512; ybox = (image[i] >> 10 & 1023) - 512; @@ -1499,12 +1504,14 @@ x1 = x[i][1] + ybox*yprd; x2 = x[i][2] + zbox*zprd; - virial[0] += 0.5*fc0*x0; - virial[1] += 0.5*fc1*x1; - virial[2] += 0.5*fc2*x2; - virial[3] += 0.5*fc1*x0; - virial[4] += 0.5*fc2*x0; - virial[5] += 0.5*fc2*x1; + vr[0] = 0.5*fc0*x0; + vr[1] = 0.5*fc1*x1; + vr[2] = 0.5*fc2*x2; + vr[3] = 0.5*fc1*x0; + vr[4] = 0.5*fc2*x0; + vr[5] = 0.5*fc2*x1; + + v_tally(1,&i,1.0,vr); } } } diff -Naur lammps-16Jan08/src/POEMS/fix_poems.h lammps-17Jan08/src/POEMS/fix_poems.h --- lammps-16Jan08/src/POEMS/fix_poems.h 2007-10-10 16:22:16.000000000 -0600 +++ lammps-17Jan08/src/POEMS/fix_poems.h 2008-01-17 10:01:57.000000000 -0700 @@ -24,11 +24,11 @@ ~FixPOEMS(); int setmask(); void init(); - void setup(); - void initial_integrate(); + void setup(int); + void initial_integrate(int); void post_force(int); void final_integrate(); - void initial_integrate_respa(int, int); + void initial_integrate_respa(int, int, int); void post_force_respa(int, int, int); void final_integrate_respa(int); @@ -48,7 +48,6 @@ double dtv,dtf,dthalf; double *step_respa; int nlevels_respa; - int pressure_flag; double total_ke; // atom assignment to rigid bodies @@ -100,8 +99,8 @@ void rotate(double **, int, int, int, int, double, double); void omega_from_mq(double *, double *, double *, double *, double *, double *); - void set_v(int); - void set_xv(int); + void set_v(); + void set_xv(); }; } diff -Naur lammps-16Jan08/src/compute_stress_atom.cpp lammps-17Jan08/src/compute_stress_atom.cpp --- lammps-16Jan08/src/compute_stress_atom.cpp 2008-01-09 14:56:57.000000000 -0700 +++ lammps-17Jan08/src/compute_stress_atom.cpp 2008-01-17 10:01:57.000000000 -0700 @@ -117,7 +117,7 @@ stress[i][j] = 0.0; // add in per-atom contributions from each force - + if (pairflag && force->pair) { double **vatom = force->pair->vatom; for (i = 0; i < npair; i++) diff -Naur lammps-16Jan08/src/fix_rigid.cpp lammps-17Jan08/src/fix_rigid.cpp --- lammps-16Jan08/src/fix_rigid.cpp 2008-01-09 14:56:57.000000000 -0700 +++ lammps-17Jan08/src/fix_rigid.cpp 2008-01-17 10:01:57.000000000 -0700 @@ -296,14 +296,6 @@ error->all("Rigid fix must come before NPT/NPH fix"); } - // compute rigid contribution to virial every step if fix NPT,NPH exists - - pressure_flag = 0; - for (int i = 0; i < modify->nfix; i++) { - if (strcmp(modify->fix[i]->style,"npt") == 0) pressure_flag = 1; - if (strcmp(modify->fix[i]->style,"nph") == 0) pressure_flag = 1; - } - // timestep info dtv = update->dt; @@ -720,7 +712,7 @@ if (vflag) v_setup(vflag); else evflag = 0; - // set coords and velocities if atoms in rigid bodies + // set coords and velocities of atoms in rigid bodies // from quarternion and omega set_xv(); @@ -1240,6 +1232,9 @@ int *image = atom->image; double **x = atom->x; double **v = atom->v; + double **f = atom->f; + double *mass = atom->mass; + int *type = atom->type; int nlocal = atom->nlocal; double xprd = domain->xprd; @@ -1251,10 +1246,8 @@ xz = domain->xz; yz = domain->yz; } - - double *mass = atom->mass; - double **f = atom->f; - int *type = atom->type; + + // set x and v of each atom for (int i = 0; i < nlocal; i++) { if (body[i] < 0) continue; @@ -1340,7 +1333,7 @@ } /* ---------------------------------------------------------------------- - set space-frame velocity of each atom in rigid body + set space-frame velocity of each atom in a rigid body v = Vcm + (W cross (x - Xcm)) ------------------------------------------------------------------------- */ @@ -1348,7 +1341,7 @@ { int ibody; int xbox,ybox,zbox; - double xunwrap,yunwrap,zunwrap,dx,dy,dz; + double dx,dy,dz; double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone; double xy,xz,yz; double vr[6]; @@ -1370,6 +1363,8 @@ yz = domain->yz; } + // set v of each atom + for (int i = 0; i < nlocal; i++) { if (body[i] < 0) continue; ibody = body[i]; diff -Naur lammps-16Jan08/src/fix_rigid.h lammps-17Jan08/src/fix_rigid.h --- lammps-16Jan08/src/fix_rigid.h 2008-01-09 14:56:57.000000000 -0700 +++ lammps-17Jan08/src/fix_rigid.h 2008-01-17 10:01:57.000000000 -0700 @@ -44,7 +44,6 @@ private: double dtv,dtf,dtq; double *step_respa; - int pressure_flag; int triclinic; int nbody; // # of rigid bodies diff -Naur lammps-16Jan08/src/lammps.cpp lammps-17Jan08/src/lammps.cpp --- lammps-16Jan08/src/lammps.cpp 2007-10-16 12:08:31.000000000 -0600 +++ lammps-17Jan08/src/lammps.cpp 2008-01-17 10:01:57.000000000 -0700 @@ -44,6 +44,7 @@ memory = new Memory(this); error = new Error(this); universe = new Universe(this,communicator); + output = NULL; // parse input switches diff -Naur lammps-16Jan08/src/pair.cpp lammps-17Jan08/src/pair.cpp --- lammps-16Jan08/src/pair.cpp 2008-01-02 12:24:46.000000000 -0700 +++ lammps-17Jan08/src/pair.cpp 2008-01-17 10:01:57.000000000 -0700 @@ -485,13 +485,12 @@ /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global and per-atom accumulators - called by SW and Tersoff potentials, newton_pair is always on - virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drij*fj + drik*fk - could just pass ri,rj,rk since coords are already unwrapped by PBC + called by SW potential, newton_pair is always on + virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drji*fj + drki*fk ------------------------------------------------------------------------- */ void Pair::ev_tally3(int i, int j, int k, double evdwl, double ecoul, - double *fj, double *fk, double *drij, double *drik) + double *fj, double *fk, double *drji, double *drki) { double epairthird,v[6]; @@ -509,12 +508,12 @@ } if (vflag_atom) { - v[0] = THIRD * (drij[0]*fj[0] + drik[0]*fk[0]); - v[1] = THIRD * (drij[1]*fj[1] + drik[1]*fk[1]); - v[2] = THIRD * (drij[2]*fj[2] + drik[2]*fk[2]); - v[3] = THIRD * (drij[0]*fj[1] + drik[0]*fk[1]); - v[4] = THIRD * (drij[0]*fj[2] + drik[0]*fk[2]); - v[5] = THIRD * (drij[1]*fj[2] + drik[1]*fk[2]); + v[0] = THIRD * (drji[0]*fj[0] + drki[0]*fk[0]); + v[1] = THIRD * (drji[1]*fj[1] + drki[1]*fk[1]); + v[2] = THIRD * (drji[2]*fj[2] + drki[2]*fk[2]); + v[3] = THIRD * (drji[0]*fj[1] + drki[0]*fk[1]); + v[4] = THIRD * (drji[0]*fj[2] + drki[0]*fk[2]); + v[5] = THIRD * (drji[1]*fj[2] + drki[1]*fk[2]); vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; @@ -526,22 +525,133 @@ } /* ---------------------------------------------------------------------- + tally eng_vdwl and virial into global and per-atom accumulators + called by AIREBO potential, newton_pair is always on + ------------------------------------------------------------------------- */ + +void Pair::ev_tally4(int i, int j, int k, int m, double evdwl, + double *fi, double *fj, double *fk, + double *drim, double *drjm, double *drkm) +{ + double epairfourth,v[6]; + + if (eflag_either) { + if (eflag_global) eng_vdwl += evdwl; + if (eflag_atom) { + epairfourth = 0.25 * evdwl; + eatom[i] += epairfourth; + eatom[j] += epairfourth; + eatom[k] += epairfourth; + eatom[m] += epairfourth; + } + } + + if (vflag_atom) { + v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); + v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); + v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); + v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); + v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); + v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); + + vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; + vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; + vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; + vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; + vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2]; + vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5]; + vatom[m][0] += v[0]; vatom[m][1] += v[1]; vatom[m][2] += v[2]; + vatom[m][3] += v[3]; vatom[m][4] += v[4]; vatom[m][5] += v[5]; + } +} + +/* ---------------------------------------------------------------------- + tally ecoul and virial into each of n atoms in list + called by TIP4P potential, newton_pair is always on + changes v values by dividing by n + ------------------------------------------------------------------------- */ + +void Pair::ev_tally_list(int n, int *list, double ecoul, double *v) +{ + int i,j; + + if (eflag_either) { + if (eflag_global) eng_coul += ecoul; + if (eflag_atom) { + double epairatom = ecoul/n; + for (i = 0; i < n; i++) eatom[list[i]] += epairatom; + } + } + + if (vflag_either) { + if (vflag_global) { + virial[0] += v[0]; + virial[1] += v[1]; + virial[2] += v[2]; + virial[3] += v[3]; + virial[4] += v[4]; + virial[5] += v[5]; + } + + if (vflag_atom) { + v[0] /= n; + v[1] /= n; + v[2] /= n; + v[3] /= n; + v[4] /= n; + v[5] /= n; + for (i = 0; i < n; i++) { + j = list[i]; + vatom[j][0] += v[0]; + vatom[j][1] += v[1]; + vatom[j][2] += v[2]; + vatom[j][3] += v[3]; + vatom[j][4] += v[4]; + vatom[j][5] += v[5]; + } + } + } +} + +/* ---------------------------------------------------------------------- + tally virial into per-atom accumulators + called by AIREBO potential, newton_pair is always on + fpair is magnitude of force on atom I +------------------------------------------------------------------------- */ + +void Pair::v_tally2(int i, int j, double fpair, double *drij) +{ + double v[6]; + + v[0] = 0.5 * drij[0]*drij[0]*fpair; + v[1] = 0.5 * drij[1]*drij[1]*fpair; + v[2] = 0.5 * drij[2]*drij[2]*fpair; + v[3] = 0.5 * drij[0]*drij[1]*fpair; + v[4] = 0.5 * drij[0]*drij[2]*fpair; + v[5] = 0.5 * drij[1]*drij[2]*fpair; + + vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; + vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; + vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; + vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; +} + +/* ---------------------------------------------------------------------- tally virial into per-atom accumulators - called by airebo potential, newton_pair is always on - could just pass ri,rj,rk since coords are already unwrapped by PBC + called by AIREBO and Tersoff potential, newton_pair is always on ------------------------------------------------------------------------- */ void Pair::v_tally3(int i, int j, int k, - double f1, double f2, double *drik, double *drkj) + double *fi, double *fj, double *drik, double *drjk) { double v[6]; - v[0] = 0.0; - v[1] = 0.0; - v[2] = 0.0; - v[3] = 0.0; - v[4] = 0.0; - v[5] = 0.0; + v[0] = THIRD * (drik[0]*fi[0] + drjk[0]*fj[0]); + v[1] = THIRD * (drik[1]*fi[1] + drjk[1]*fj[1]); + v[2] = THIRD * (drik[2]*fi[2] + drjk[2]*fj[2]); + v[3] = THIRD * (drik[0]*fi[1] + drjk[0]*fj[1]); + v[4] = THIRD * (drik[0]*fi[2] + drjk[0]*fj[2]); + v[5] = THIRD * (drik[1]*fi[2] + drjk[1]*fj[2]); vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; @@ -553,22 +663,21 @@ /* ---------------------------------------------------------------------- tally virial into per-atom accumulators - called by airebo potential, newton_pair is always on - could just pass ri,rj,rk,rm since coords are already unwrapped by PBC + called by AIREBO potential, newton_pair is always on ------------------------------------------------------------------------- */ void Pair::v_tally4(int i, int j, int k, int m, - double f1, double f2, double f3, - double *drik, double *drkm, double *drmj) + double *fi, double *fj, double *fk, + double *drim, double *drjm, double *drkm) { double v[6]; - v[0] = 0.0; - v[1] = 0.0; - v[2] = 0.0; - v[3] = 0.0; - v[4] = 0.0; - v[5] = 0.0; + v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); + v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); + v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); + v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); + v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); + v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; diff -Naur lammps-16Jan08/src/pair.h lammps-17Jan08/src/pair.h --- lammps-16Jan08/src/pair.h 2008-01-02 12:24:46.000000000 -0700 +++ lammps-17Jan08/src/pair.h 2008-01-17 10:01:57.000000000 -0700 @@ -119,8 +119,12 @@ double, double, double, double, double, double); void ev_tally3(int, int, int, double, double, double *, double *, double *, double *); - void v_tally3(int, int, int, double, double, double *, double *); - void v_tally4(int, int, int, int, double, double, double, + void ev_tally4(int, int, int, int, double, + double *, double *, double *, double *, double *, double *); + void ev_tally_list(int, int *, double, double *); + void v_tally2(int, int, double, double *); + void v_tally3(int, int, int, double *, double *, double *, double *); + void v_tally4(int, int, int, int, double *, double *, double *, double *, double *, double *); void virial_compute(); }; diff -Naur lammps-16Jan08/src/thermo.cpp lammps-17Jan08/src/thermo.cpp --- lammps-16Jan08/src/thermo.cpp 2008-01-02 12:24:46.000000000 -0700 +++ lammps-17Jan08/src/thermo.cpp 2008-01-17 10:01:57.000000000 -0700 @@ -256,7 +256,7 @@ int ivariable; for (i = 0; i < nvariable; i++) { ivariable = input->variable->find(id_variable[i]); - if (ivariable < 0) error->all("Could not find thermo variable name"); + if (ivariable < 0) error->all("Could not find thermo custom variable name"); variables[i] = ivariable; } diff -Naur lammps-16Jan08/src/variable.cpp lammps-17Jan08/src/variable.cpp --- lammps-16Jan08/src/variable.cpp 2008-01-09 08:39:05.000000000 -0700 +++ lammps-17Jan08/src/variable.cpp 2008-01-17 10:01:57.000000000 -0700 @@ -1094,7 +1094,7 @@ return exp(eval_tree(tree->left,i)); if (tree->type == LN) { double arg = eval_tree(tree->left,i); - if (arg <= 0.0) error->all("Log of negative/zero in variable formula"); + if (arg <= 0.0) error->all("Log of zero/negative in variable formula"); return log(arg); }