diff --git a/doc/src/atom_modify.txt b/doc/src/atom_modify.txt index d5c82f16a..1dc0fa6bf 100644 --- a/doc/src/atom_modify.txt +++ b/doc/src/atom_modify.txt @@ -1,170 +1,174 @@ "LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c :link(lws,http://lammps.sandia.gov) :link(ld,Manual.html) :link(lc,Section_commands.html#comm) :line atom_modify command :h3 [Syntax:] atom_modify keyword values ... :pre one or more keyword/value pairs may be appended :ulb,l keyword = {id} or {map} or {first} or {sort} :l {id} value = {yes} or {no} - {map} value = {array} or {hash} + {map} value = {yes} or {array} or {hash} {first} value = group-ID = group whose atoms will appear first in internal atom lists {sort} values = Nfreq binsize Nfreq = sort atoms spatially every this many time steps binsize = bin size for spatial sorting (distance units) :pre :ule [Examples:] -atom_modify map hash -atom_modify map array sort 10000 2.0 +atom_modify map yes +atom_modify map hash sort 10000 2.0 atom_modify first colloid :pre [Description:] Modify certain attributes of atoms defined and stored within LAMMPS, in addition to what is specified by the "atom_style"_atom_style.html command. The {id} and {map} keywords must be specified before a simulation box is defined; other keywords can be specified any time. The {id} keyword determines whether non-zero atom IDs can be assigned to each atom. If the value is {yes}, which is the default, IDs are assigned, whether you use the "create atoms"_create_atoms.html or "read_data"_read_data.html or "read_restart"_read_restart.html commands to initialize atoms. If the value is {no} the IDs for all atoms are assumed to be 0. If atom IDs are used, they must all be positive integers. They should also be unique, though LAMMPS does not check for this. Typically they should also be consecutively numbered (from 1 to Natoms), though this is not required. Molecular "atom styles"_atom_style.html are those that store bond topology information (styles bond, angle, molecular, full). These styles require atom IDs since the IDs are used to encode the topology. Some other LAMMPS commands also require the use of atom IDs. E.g. some many-body pair styles use them to avoid double computation of the I-J interaction between two atoms. The only reason not to use atom IDs is if you are running an atomic simulation so large that IDs cannot be uniquely assigned. For a default LAMMPS build this limit is 2^31 or about 2 billion atoms. However, even in this case, you can use 64-bit atom IDs, allowing 2^63 or about 9e18 atoms, if you build LAMMPS with the - DLAMMPS_BIGBIG switch. This is described in "Section 2.2"_Section_start.html#start_2 of the manual. If atom IDs are not used, they must be specified as 0 for all atoms, e.g. in a data or restart file. -The {map} keyword determines how atom ID lookup is done for molecular -atom styles. Lookups are performed by bond (angle, etc) routines in -LAMMPS to find the local atom index associated with a global atom ID. - -When the {array} value is used, each processor stores a lookup table -of length N, where N is the largest atom ID in the system. This is a +The {map} keyword determines how atoms with specific IDs are found +when required. An example are the bond (angle, etc) methods which +need to find the local index of an atom with a specific global ID +which is a bond (angle, etc) partner. LAMMPS performs this operation +efficiently by creating a "map", which is either an {array} or {hash} +table, as descibed below. + +When the {map} keyword is not specified in your input script, LAMMPS +only creates a map for "atom_styles"_atom_style.html for molecular +systems which have permanent bonds (angles, etc). No map is created +for atomic systems, since it is normally not needed. However some +LAMMPS commands require a map, even for atomic systems, and will +generate an error if one does not exist. The {map} keyword thus +allows you to force the creation of a map. The {yes} value will +create either an {array} or {hash} style map, as explained in the next +paragraph. The {array} and {hash} values create an atom-style or +hash-style map respectively. + +For an {array}-style map, each processor stores a lookup table of +length N, where N is the largest atom ID in the system. This is a fast, simple method for many simulations, but requires too much memory -for large simulations. The {hash} value uses a hash table to perform -the lookups. This can be slightly slower than the {array} method, but -its memory cost is proportional to the number of atoms owned by a -processor, i.e. N/P when N is the total number of atoms in the system -and P is the number of processors. - -When this setting is not specified in your input script, LAMMPS -creates a map, if one is needed, as an array or hash. See the -discussion of default values below for how LAMMPS chooses which kind -of map to build. Note that atomic systems do not normally need to -create a map. However, even in this case some LAMMPS commands will -create a map to find atoms (and then destroy it), or require a -permanent map. An example of the former is the "velocity loop -all"_velocity.html command, which uses a map when looping over all -atoms and insuring the same velocity values are assigned to an atom -ID, no matter which processor owns it. +for large simulations. For a {hash}-style map, a hash table is +created on each processor, which finds an atom ID in constant time +(independent of the global number of atom IDs). It can be slightly +slower than the {array} map, but its memory cost is proportional to +the number of atoms owned by a processor, i.e. N/P when N is the total +number of atoms in the system and P is the number of processors. The {first} keyword allows a "group"_group.html to be specified whose atoms will be maintained as the first atoms in each processor's list of owned atoms. This in only useful when the specified group is a small fraction of all the atoms, and there are other operations LAMMPS is performing that will be sped-up significantly by being able to loop over the smaller set of atoms. Otherwise the reordering required by this option will be a net slow-down. The "neigh_modify include"_neigh_modify.html and "comm_modify group"_comm_modify.html commands are two examples of commands that require this setting to work efficiently. Several "fixes"_fix.html, most notably time integration fixes like "fix nve"_fix_nve.html, also take advantage of this setting if the group they operate on is the group specified by this command. Note that specifying "all" as the group-ID effectively turns off the {first} option. It is OK to use the {first} keyword with a group that has not yet been defined, e.g. to use the atom_modify first command at the beginning of your input script. LAMMPS does not use the group until a simulation is run. The {sort} keyword turns on a spatial sorting or reordering of atoms within each processor's sub-domain every {Nfreq} timesteps. If {Nfreq} is set to 0, then sorting is turned off. Sorting can improve cache performance and thus speed-up a LAMMPS simulation, as discussed in a paper by "(Meloni)"_#Meloni. Its efficacy depends on the problem size (atoms/processor), how quickly the system becomes disordered, and various other factors. As a general rule, sorting is typically more effective at speeding up simulations of liquids as opposed to solids. In tests we have done, the speed-up can range from zero to 3-4x. Reordering is performed every {Nfreq} timesteps during a dynamics run or iterations during a minimization. More precisely, reordering occurs at the first reneighboring that occurs after the target timestep. The reordering is performed locally by each processor, using bins of the specified {binsize}. If {binsize} is set to 0.0, then a binsize equal to half the "neighbor"_neighbor.html cutoff distance (force cutoff plus skin distance) is used, which is a reasonable value. After the atoms have been binned, they are reordered so that atoms in the same bin are adjacent to each other in the processor's 1d list of atoms. The goal of this procedure is for atoms to put atoms close to each other in the processor's one-dimensional list of atoms that are also near to each other spatially. This can improve cache performance when pairwise interactions and neighbor lists are computed. Note that if bins are too small, there will be few atoms/bin. Likewise if bins are too large, there will be many atoms/bin. In both cases, the goal of cache locality will be undermined. NOTE: Running a simulation with sorting on versus off should not change the simulation results in a statistical sense. However, a different ordering will induce round-off differences, which will lead to diverging trajectories over time when comparing two simulations. Various commands, particularly those which use random numbers (e.g. "velocity create"_velocity.html, and "fix langevin"_fix_langevin.html), may generate (statistically identical) results which depend on the order in which atoms are processed. The order of atoms in a "dump"_dump.html file will also typically change if sorting is enabled. [Restrictions:] The {first} and {sort} options cannot be used together. Since sorting is on by default, it will be turned off if the {first} keyword is used with a group-ID that is not "all". [Related commands:] none [Default:] By default, {id} is yes. By default, atomic systems (no bond topology info) do not use a map. For molecular systems (with bond topology info), a map is used. The default map style is array if no atom ID is larger than 1 million, otherwise the default is hash. By default, a "first" group is not defined. By default, sorting is enabled with a frequency of 1000 and a binsize of 0.0, which means the neighbor cutoff will be used to set the bin size. :line :link(Meloni) [(Meloni)] Meloni, Rosati and Colombo, J Chem Phys, 126, 121102 (2007). diff --git a/doc/src/fix_nh.txt b/doc/src/fix_nh.txt index 8fa30ac22..41d0e6438 100644 --- a/doc/src/fix_nh.txt +++ b/doc/src/fix_nh.txt @@ -1,642 +1,646 @@ <"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c :link(lws,http://lammps.sandia.gov) :link(ld,Manual.html) :link(lc,Section_commands.html#comm) :line fix nvt command :h3 fix nvt/intel command :h3 fix nvt/kk command :h3 fix nvt/omp command :h3 fix npt command :h3 fix npt/intel command :h3 fix npt/kk command :h3 fix npt/omp command :h3 fix nph command :h3 fix nph/kk command :h3 fix nph/omp command :h3 [Syntax:] fix ID group-ID style_name keyword value ... :pre ID, group-ID are documented in "fix"_fix.html command :ulb,l style_name = {nvt} or {npt} or {nph} :l one or more keyword/value pairs may be appended :l keyword = {temp} or {iso} or {aniso} or {tri} or {x} or {y} or {z} or {xy} or {yz} or {xz} or {couple} or {tchain} or {pchain} or {mtk} or {tloop} or {ploop} or {nreset} or {drag} or {dilate} or {scalexy} or {scaleyz} or {scalexz} or {flip} or {fixedpoint} or {update} {temp} values = Tstart Tstop Tdamp Tstart,Tstop = external temperature at start/end of run Tdamp = temperature damping parameter (time units) {iso} or {aniso} or {tri} values = Pstart Pstop Pdamp Pstart,Pstop = scalar external pressure at start/end of run (pressure units) Pdamp = pressure damping parameter (time units) {x} or {y} or {z} or {xy} or {yz} or {xz} values = Pstart Pstop Pdamp Pstart,Pstop = external stress tensor component at start/end of run (pressure units) Pdamp = stress damping parameter (time units) {couple} = {none} or {xyz} or {xy} or {yz} or {xz} {tchain} value = N N = length of thermostat chain (1 = single thermostat) {pchain} values = N N length of thermostat chain on barostat (0 = no thermostat) {mtk} value = {yes} or {no} = add in MTK adjustment term or not {tloop} value = M M = number of sub-cycles to perform on thermostat {ploop} value = M M = number of sub-cycles to perform on barostat thermostat {nreset} value = reset reference cell every this many timesteps {drag} value = Df Df = drag factor added to barostat/thermostat (0.0 = no drag) {dilate} value = dilate-group-ID dilate-group-ID = only dilate atoms in this group due to barostat volume changes {scalexy} value = {yes} or {no} = scale xy with ly {scaleyz} value = {yes} or {no} = scale yz with lz {scalexz} value = {yes} or {no} = scale xz with lz {flip} value = {yes} or {no} = allow or disallow box flips when it becomes highly skewed {fixedpoint} values = x y z x,y,z = perform barostat dilation/contraction around this point (distance units) {update} value = {dipole} or {dipole/dlm} dipole = update dipole orientation (only for sphere variants) dipole/dlm = use DLM integrator to update dipole orientation (only for sphere variants) :pre :ule [Examples:] fix 1 all nvt temp 300.0 300.0 100.0 fix 1 water npt temp 300.0 300.0 100.0 iso 0.0 0.0 1000.0 fix 2 jello npt temp 300.0 300.0 100.0 tri 5.0 5.0 1000.0 fix 2 ice nph x 1.0 1.0 0.5 y 2.0 2.0 0.5 z 3.0 3.0 0.5 yz 0.1 0.1 0.5 xz 0.2 0.2 0.5 xy 0.3 0.3 0.5 nreset 1000 :pre [Description:] These commands perform time integration on Nose-Hoover style non-Hamiltonian equations of motion which are designed to generate positions and velocities sampled from the canonical (nvt), isothermal-isobaric (npt), and isenthalpic (nph) ensembles. This updates the position and velocity for atoms in the group each timestep. The thermostatting and barostatting is achieved by adding some dynamic variables which are coupled to the particle velocities (thermostatting) and simulation domain dimensions (barostatting). In addition to basic thermostatting and barostatting, these fixes can also create a chain of thermostats coupled to the particle thermostat, and another chain of thermostats coupled to the barostat variables. The barostat can be coupled to the overall box volume, or to individual dimensions, including the {xy}, {xz} and {yz} tilt dimensions. The external pressure of the barostat can be specified as either a scalar pressure (isobaric ensemble) or as components of a symmetric stress tensor (constant stress ensemble). When used correctly, the time-averaged temperature and stress tensor of the particles will match the target values specified by Tstart/Tstop and Pstart/Pstop. The equations of motion used are those of Shinoda et al in "(Shinoda)"_#nh-Shinoda, which combine the hydrostatic equations of Martyna, Tobias and Klein in "(Martyna)"_#nh-Martyna with the strain energy proposed by Parrinello and Rahman in "(Parrinello)"_#nh-Parrinello. The time integration schemes closely follow the time-reversible measure-preserving Verlet and rRESPA integrators derived by Tuckerman et al in "(Tuckerman)"_#nh-Tuckerman. :line The thermostat parameters for fix styles {nvt} and {npt} is specified using the {temp} keyword. Other thermostat-related keywords are {tchain}, {tloop} and {drag}, which are discussed below. The thermostat is applied to only the translational degrees of freedom for the particles. The translational degrees of freedom can also have a bias velocity removed before thermostatting takes place; see the description below. The desired temperature at each timestep is a ramped value during the run from {Tstart} to {Tstop}. The {Tdamp} parameter is specified in time units and determines how rapidly the temperature is relaxed. For example, a value of 10.0 means to relax the temperature in a timespan of (roughly) 10 time units (e.g. tau or fmsec or psec - see the "units"_units.html command). The atoms in the fix group are the only ones whose velocities and positions are updated by the velocity/position update portion of the integration. NOTE: A Nose-Hoover thermostat will not work well for arbitrary values of {Tdamp}. If {Tdamp} is too small, the temperature can fluctuate wildly; if it is too large, the temperature will take a very long time to equilibrate. A good choice for many models is a {Tdamp} of around 100 timesteps. Note that this is NOT the same as 100 time units for most "units"_units.html settings. :line The barostat parameters for fix styles {npt} and {nph} is specified using one or more of the {iso}, {aniso}, {tri}, {x}, {y}, {z}, {xy}, {xz}, {yz}, and {couple} keywords. These keywords give you the ability to specify all 6 components of an external stress tensor, and to couple various of these components together so that the dimensions they represent are varied together during a constant-pressure simulation. Other barostat-related keywords are {pchain}, {mtk}, {ploop}, {nreset}, {drag}, and {dilate}, which are discussed below. Orthogonal simulation boxes have 3 adjustable dimensions (x,y,z). Triclinic (non-orthogonal) simulation boxes have 6 adjustable dimensions (x,y,z,xy,xz,yz). The "create_box"_create_box.html, "read data"_read_data.html, and "read_restart"_read_restart.html commands specify whether the simulation box is orthogonal or non-orthogonal (triclinic) and explain the meaning of the xy,xz,yz tilt factors. The target pressures for each of the 6 components of the stress tensor can be specified independently via the {x}, {y}, {z}, {xy}, {xz}, {yz} keywords, which correspond to the 6 simulation box dimensions. For each component, the external pressure or tensor component at each timestep is a ramped value during the run from {Pstart} to {Pstop}. If a target pressure is specified for a component, then the corresponding box dimension will change during a simulation. For example, if the {y} keyword is used, the y-box length will change. If the {xy} keyword is used, the xy tilt factor will change. A box dimension will not change if that component is not specified, although you have the option to change that dimension via the "fix deform"_fix_deform.html command. Note that in order to use the {xy}, {xz}, or {yz} keywords, the simulation box must be triclinic, even if its initial tilt factors are 0.0. For all barostat keywords, the {Pdamp} parameter operates like the {Tdamp} parameter, determining the time scale on which pressure is relaxed. For example, a value of 10.0 means to relax the pressure in a timespan of (roughly) 10 time units (e.g. tau or fmsec or psec - see the "units"_units.html command). NOTE: A Nose-Hoover barostat will not work well for arbitrary values of {Pdamp}. If {Pdamp} is too small, the pressure and volume can fluctuate wildly; if it is too large, the pressure will take a very long time to equilibrate. A good choice for many models is a {Pdamp} of around 1000 timesteps. However, note that {Pdamp} is specified in time units, and that timesteps are NOT the same as time units for most "units"_units.html settings. Regardless of what atoms are in the fix group (the only atoms which are time integrated), a global pressure or stress tensor is computed for all atoms. Similarly, when the size of the simulation box is changed, all atoms are re-scaled to new positions, unless the keyword {dilate} is specified with a {dilate-group-ID} for a group that represents a subset of the atoms. This can be useful, for example, to leave the coordinates of atoms in a solid substrate unchanged and controlling the pressure of a surrounding fluid. This option should be used with care, since it can be unphysical to dilate some atoms and not others, because it can introduce large, instantaneous displacements between a pair of atoms (one dilated, one not) that are far from the dilation origin. Also note that for atoms not in the fix group, a separate time integration fix like "fix nve"_fix_nve.html or "fix nvt"_fix_nh.html can be used on them, independent of whether they are dilated or not. :line The {couple} keyword allows two or three of the diagonal components of the pressure tensor to be "coupled" together. The value specified with the keyword determines which are coupled. For example, {xz} means the {Pxx} and {Pzz} components of the stress tensor are coupled. {Xyz} means all 3 diagonal components are coupled. Coupling means two things: the instantaneous stress will be computed as an average of the corresponding diagonal components, and the coupled box dimensions will be changed together in lockstep, meaning coupled dimensions will be dilated or contracted by the same percentage every timestep. The {Pstart}, {Pstop}, {Pdamp} parameters for any coupled dimensions must be identical. {Couple xyz} can be used for a 2d simulation; the {z} dimension is simply ignored. :line The {iso}, {aniso}, and {tri} keywords are simply shortcuts that are equivalent to specifying several other keywords together. The keyword {iso} means couple all 3 diagonal components together when pressure is computed (hydrostatic pressure), and dilate/contract the dimensions together. Using "iso Pstart Pstop Pdamp" is the same as specifying these 4 keywords: x Pstart Pstop Pdamp y Pstart Pstop Pdamp z Pstart Pstop Pdamp couple xyz :pre The keyword {aniso} means {x}, {y}, and {z} dimensions are controlled independently using the {Pxx}, {Pyy}, and {Pzz} components of the stress tensor as the driving forces, and the specified scalar external pressure. Using "aniso Pstart Pstop Pdamp" is the same as specifying these 4 keywords: x Pstart Pstop Pdamp y Pstart Pstop Pdamp z Pstart Pstop Pdamp couple none :pre The keyword {tri} means {x}, {y}, {z}, {xy}, {xz}, and {yz} dimensions are controlled independently using their individual stress components as the driving forces, and the specified scalar pressure as the external normal stress. Using "tri Pstart Pstop Pdamp" is the same as specifying these 7 keywords: x Pstart Pstop Pdamp y Pstart Pstop Pdamp z Pstart Pstop Pdamp xy 0.0 0.0 Pdamp yz 0.0 0.0 Pdamp xz 0.0 0.0 Pdamp couple none :pre :line In some cases (e.g. for solids) the pressure (volume) and/or temperature of the system can oscillate undesirably when a Nose/Hoover barostat and thermostat is applied. The optional {drag} keyword will damp these oscillations, although it alters the Nose/Hoover equations. A value of 0.0 (no drag) leaves the Nose/Hoover formalism unchanged. A non-zero value adds a drag term; the larger the value specified, the greater the damping effect. Performing a short run and monitoring the pressure and temperature is the best way to determine if the drag term is working. Typically a value between 0.2 to 2.0 is sufficient to damp oscillations after a few periods. Note that use of the drag keyword will interfere with energy conservation and will also change the distribution of positions and velocities so that they do not correspond to the nominal NVT, NPT, or NPH ensembles. An alternative way to control initial oscillations is to use chain thermostats. The keyword {tchain} determines the number of thermostats in the particle thermostat. A value of 1 corresponds to the original Nose-Hoover thermostat. The keyword {pchain} specifies the number of thermostats in the chain thermostatting the barostat degrees of freedom. A value of 0 corresponds to no thermostatting of the barostat variables. The {mtk} keyword controls whether or not the correction terms due to Martyna, Tuckerman, and Klein are included in the equations of motion "(Martyna)"_#nh-Martyna. Specifying {no} reproduces the original Hoover barostat, whose volume probability distribution function differs from the true NPT and NPH ensembles by a factor of 1/V. Hence using {yes} is more correct, but in many cases the difference is negligible. The keyword {tloop} can be used to improve the accuracy of integration scheme at little extra cost. The initial and final updates of the thermostat variables are broken up into {tloop} substeps, each of length {dt}/{tloop}. This corresponds to using a first-order Suzuki-Yoshida scheme "(Tuckerman)"_#nh-Tuckerman. The keyword {ploop} does the same thing for the barostat thermostat. The keyword {nreset} controls how often the reference dimensions used to define the strain energy are reset. If this keyword is not used, or is given a value of zero, then the reference dimensions are set to those of the initial simulation domain and are never changed. If the simulation domain changes significantly during the simulation, then the final average pressure tensor will differ significantly from the specified values of the external stress tensor. A value of {nstep} means that every {nstep} timesteps, the reference dimensions are set to those of the current simulation domain. The {scaleyz}, {scalexz}, and {scalexy} keywords control whether or not the corresponding tilt factors are scaled with the associated box dimensions when barostatting triclinic periodic cells. The default values {yes} will turn on scaling, which corresponds to adjusting the linear dimensions of the cell while preserving its shape. Choosing {no} ensures that the tilt factors are not scaled with the box dimensions. See below for restrictions and default values in different situations. In older versions of LAMMPS, scaling of tilt factors was not performed. The old behavior can be recovered by setting all three scale keywords to {no}. The {flip} keyword allows the tilt factors for a triclinic box to exceed half the distance of the parallel box length, as discussed below. If the {flip} value is set to {yes}, the bound is enforced by flipping the box when it is exceeded. If the {flip} value is set to {no}, the tilt will continue to change without flipping. Note that if applied stress induces large deformations (e.g. in a liquid), this means the box shape can tilt dramatically and LAMMPS will run less efficiently, due to the large volume of communication needed to acquire ghost atoms around a processor's irregular-shaped sub-domain. For extreme values of tilt, LAMMPS may also lose atoms and generate an error. The {fixedpoint} keyword specifies the fixed point for barostat volume changes. By default, it is the center of the box. Whatever point is chosen will not move during the simulation. For example, if the lower periodic boundaries pass through (0,0,0), and this point is provided to {fixedpoint}, then the lower periodic boundaries will remain at (0,0,0), while the upper periodic boundaries will move twice as far. In all cases, the particle trajectories are unaffected by the chosen value, except for a time-dependent constant translation of positions. If the {update} keyword is used with the {dipole} value, then the orientation of the dipole moment of each particle is also updated during the time integration. This option should be used for models where a dipole moment is assigned to finite-size particles, e.g. spheroids via use of the "atom_style hybrid sphere dipole"_atom_style.html command. The default dipole orientation integrator can be changed to the Dullweber-Leimkuhler-McLachlan integration scheme "(Dullweber)"_#nh-Dullweber when using {update} with the value {dipole/dlm}. This integrator is symplectic and time-reversible, giving better energy conservation and allows slightly longer timesteps at only a small additional computational cost. :line NOTE: Using a barostat coupled to tilt dimensions {xy}, {xz}, {yz} can sometimes result in arbitrarily large values of the tilt dimensions, i.e. a dramatically deformed simulation box. LAMMPS allows the tilt factors to grow a small amount beyond the normal limit of half the box length (0.6 times the box length), and then performs a box "flip" to an equivalent periodic cell. See the discussion of the {flip} keyword above, to allow this bound to be exceeded, if desired. The flip operation is described in more detail in the doc page for "fix deform"_fix_deform.html. Both the barostat dynamics and the atom trajectories are unaffected by this operation. However, if a tilt factor is incremented by a large amount (1.5 times the box length) on a single timestep, LAMMPS can not accomodate this event and will terminate the simulation with an error. This error typically indicates that there is something badly wrong with how the simulation was constructed, such as specifying values of {Pstart} that are too far from the current stress value, or specifying a timestep that is too large. Triclinic barostatting should be used with care. This also is true for other barostat styles, although they tend to be more forgiving of insults. In particular, it is important to recognize that equilibrium liquids can not support a shear stress and that equilibrium solids can not support shear stresses that exceed the yield stress. One exception to this rule is if the 1st dimension in the tilt factor (x for xy) is non-periodic. In that case, the limits on the tilt factor are not enforced, since flipping the box in that dimension does not change the atom positions due to non-periodicity. In this mode, if you tilt the system to extreme angles, the simulation will simply become inefficient due to the highly skewed simulation box. NOTE: Unlike the "fix temp/berendsen"_fix_temp_berendsen.html command which performs thermostatting but NO time integration, these fixes perform thermostatting/barostatting AND time integration. Thus you should not use any other time integration fix, such as "fix nve"_fix_nve.html on atoms to which this fix is applied. Likewise, fix nvt and fix npt should not normally be used on atoms that also have their temperature controlled by another fix - e.g. by "fix langevin"_fix_nh.html or "fix temp/rescale"_fix_temp_rescale.html commands. See "this howto section"_Section_howto.html#howto_16 of the manual for a discussion of different ways to compute temperature and perform thermostatting and barostatting. :line These fixes compute a temperature and pressure each timestep. To do -this, the fix creates its own computes of style "temp" and "pressure", -as if one of these two sets of commands had been issued: +this, the thermostat and barostat fixes create their own computes of +style "temp" and "pressure", as if one of these sets of commands had +been issued: +For fix nvt: compute fix-ID_temp group-ID temp -compute fix-ID_press group-ID pressure fix-ID_temp :pre +For fix npt and fix nph: compute fix-ID_temp all temp compute fix-ID_press all pressure fix-ID_temp :pre -See the "compute temp"_compute_temp.html and "compute -pressure"_compute_pressure.html commands for details. Note that the -IDs of the new computes are the fix-ID + underscore + "temp" or fix_ID -+ underscore + "press". For fix nvt, the group for the new computes -is the same as the fix group. For fix nph and fix npt, the group for -the new computes is "all" since pressure is computed for the entire -system. +For fix nvt, the group for the new temperature compute is the same as +the fix group. For fix npt and fix nph, the group for both the new +temperature and pressure compute is "all" since pressure is computed +for the entire system. In the case of fix nph, the temperature +compute is not used for thermostatting, but just for a kinetic-energy +contribution to the pressure. See the "compute +temp"_compute_temp.html and "compute pressure"_compute_pressure.html +commands for details. Note that the IDs of the new computes are the +fix-ID + underscore + "temp" or fix_ID + underscore + "press". Note that these are NOT the computes used by thermodynamic output (see the "thermo_style"_thermo_style.html command) with ID = {thermo_temp} -and {thermo_press}. This means you can change the attributes of this +and {thermo_press}. This means you can change the attributes of these fix's temperature or pressure via the -"compute_modify"_compute_modify.html command or print this temperature -or pressure during thermodynamic output via the "thermo_style -custom"_thermo_style.html command using the appropriate compute-ID. -It also means that changing attributes of {thermo_temp} or -{thermo_press} will have no effect on this fix. +"compute_modify"_compute_modify.html command. Or you can print this +temperature or pressure during thermodynamic output via the +"thermo_style custom"_thermo_style.html command using the appropriate +compute-ID. It also means that changing attributes of {thermo_temp} +or {thermo_press} will have no effect on this fix. Like other fixes that perform thermostatting, fix nvt and fix npt can be used with "compute commands"_compute.html that calculate a temperature after removing a "bias" from the atom velocities. E.g. removing the center-of-mass velocity from a group of atoms or only calculating temperature on the x-component of velocity or only calculating temperature for atoms in a geometric region. This is not done by default, but only if the "fix_modify"_fix_modify.html command is used to assign a temperature compute to this fix that includes such a bias term. See the doc pages for individual "compute commands"_compute.html to determine which ones include a bias. In this case, the thermostat works in the following manner: the current temperature is calculated taking the bias into account, bias is removed from each atom, thermostatting is performed on the remaining thermal degrees of freedom, and the bias is added back in. :line These fixes can be used with either the {verlet} or {respa} "integrators"_run_style.html. When using one of the barostat fixes with {respa}, LAMMPS uses an integrator constructed according to the following factorization of the Liouville propagator (for two rRESPA levels): :c,image(Eqs/fix_nh1.jpg) This factorization differs somewhat from that of Tuckerman et al, in that the barostat is only updated at the outermost rRESPA level, whereas Tuckerman's factorization requires splitting the pressure into pieces corresponding to the forces computed at each rRESPA level. In theory, the latter method will exhibit better numerical stability. In practice, because Pdamp is normally chosen to be a large multiple of the outermost rRESPA timestep, the barostat dynamics are not the limiting factor for numerical stability. Both factorizations are time-reversible and can be shown to preserve the phase space measure of the underlying non-Hamiltonian equations of motion. NOTE: This implementation has been shown to conserve linear momentum up to machine precision under NVT dynamics. Under NPT dynamics, for a system with zero initial total linear momentum, the total momentum fluctuates close to zero. It may occasionally undergo brief excursions to non-negligible values, before returning close to zero. Over long simulations, this has the effect of causing the center-of-mass to undergo a slow random walk. This can be mitigated by resetting the momentum at infrequent intervals using the "fix momentum"_fix_momentum.html command. :line The fix npt and fix nph commands can be used with rigid bodies or mixtures of rigid bodies and non-rigid particles (e.g. solvent). But there are also "fix rigid/npt"_fix_rigid.html and "fix rigid/nph"_fix_rigid.html commands, which are typically a more natural choice. See the doc page for those commands for more discussion of the various ways to do this. :line Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are functionally the same as the corresponding style without the suffix. They have been optimized to run faster, depending on your available hardware, as discussed in "Section 5"_Section_accelerate.html of the manual. The accelerated styles take the same arguments and should produce the same results, except for round-off and precision issues. These accelerated styles are part of the GPU, USER-INTEL, KOKKOS, USER-OMP and OPT packages, respectively. They are only enabled if LAMMPS was built with those packages. See the "Making LAMMPS"_Section_start.html#start_3 section for more info. You can specify the accelerated styles explicitly in your input script by including their suffix, or you can use the "-suffix command-line switch"_Section_start.html#start_6 when you invoke LAMMPS, or you can use the "suffix"_suffix.html command in your input script. See "Section 5"_Section_accelerate.html of the manual for more instructions on how to use the accelerated styles effectively. :line [Restart, fix_modify, output, run start/stop, minimize info:] These fixes writes the state of all the thermostat and barostat variables to "binary restart files"_restart.html. See the "read_restart"_read_restart.html command for info on how to re-specify a fix in an input script that reads a restart file, so that the operation of the fix continues in an uninterrupted fashion. The "fix_modify"_fix_modify.html {temp} and {press} options are supported by these fixes. You can use them to assign a "compute"_compute.html you have defined to this fix which will be used in its thermostatting or barostatting procedure, as described above. If you do this, note that the kinetic energy derived from the compute temperature should be consistent with the virial term computed using all atoms for the pressure. LAMMPS will warn you if you choose to compute temperature on a subset of atoms. NOTE: If both the {temp} and {press} keywords are used in a single thermo_modify command (or in two separate commands), then the order in which the keywords are specified is important. Note that a "pressure compute"_compute_pressure.html defines its own temperature compute as an argument when it is specified. The {temp} keyword will override this (for the pressure compute being used by fix npt), but only if the {temp} keyword comes after the {press} keyword. If the {temp} keyword comes before the {press} keyword, then the new pressure compute specified by the {press} keyword will be unaffected by the {temp} setting. The "fix_modify"_fix_modify.html {energy} option is supported by these fixes to add the energy change induced by Nose/Hoover thermostatting and barostatting to the system's potential energy as part of "thermodynamic output"_thermo_style.html. These fixes compute a global scalar and a global vector of quantities, which can be accessed by various "output commands"_Section_howto.html#howto_15. The scalar value calculated by these fixes is "extensive"; the vector values are "intensive". The scalar is the cumulative energy change due to the fix. The vector stores internal Nose/Hoover thermostat and barostat variables. The number and meaning of the vector values depends on which fix is used and the settings for keywords {tchain} and {pchain}, which specify the number of Nose/Hoover chains for the thermostat and barostat. If no thermostatting is done, then {tchain} is 0. If no barostatting is done, then {pchain} is 0. In the following list, "ndof" is 0, 1, 3, or 6, and is the number of degrees of freedom in the barostat. Its value is 0 if no barostat is used, else its value is 6 if any off-diagonal stress tensor component is barostatted, else its value is 1 if {couple xyz} is used or {couple xy} for a 2d simulation, otherwise its value is 3. The order of values in the global vector and their meaning is as follows. The notation means there are tchain values for eta, followed by tchain for eta_dot, followed by ndof for omega, etc: eta\[tchain\] = particle thermostat displacements (unitless) eta_dot\[tchain\] = particle thermostat velocities (1/time units) omega\[ndof\] = barostat displacements (unitless) omega_dot\[ndof\] = barostat velocities (1/time units) etap\[pchain\] = barostat thermostat displacements (unitless) etap_dot\[pchain\] = barostat thermostat velocities (1/time units) PE_eta\[tchain\] = potential energy of each particle thermostat displacement (energy units) KE_eta_dot\[tchain\] = kinetic energy of each particle thermostat velocity (energy units) PE_omega\[ndof\] = potential energy of each barostat displacement (energy units) KE_omega_dot\[ndof\] = kinetic energy of each barostat velocity (energy units) PE_etap\[pchain\] = potential energy of each barostat thermostat displacement (energy units) KE_etap_dot\[pchain\] = kinetic energy of each barostat thermostat velocity (energy units) PE_strain\[1\] = scalar strain energy (energy units) :ul These fixes can ramp their external temperature and pressure over multiple runs, using the {start} and {stop} keywords of the "run"_run.html command. See the "run"_run.html command for details of how to do this. These fixes are not invoked during "energy minimization"_minimize.html. :line [Restrictions:] {X}, {y}, {z} cannot be barostatted if the associated dimension is not periodic. {Xy}, {xz}, and {yz} can only be barostatted if the simulation domain is triclinic and the 2nd dimension in the keyword ({y} dimension in {xy}) is periodic. {Z}, {xz}, and {yz}, cannot be barostatted for 2D simulations. The "create_box"_create_box.html, "read data"_read_data.html, and "read_restart"_read_restart.html commands specify whether the simulation box is orthogonal or non-orthogonal (triclinic) and explain the meaning of the xy,xz,yz tilt factors. For the {temp} keyword, the final Tstop cannot be 0.0 since it would make the external T = 0.0 at some timestep during the simulation which is not allowed in the Nose/Hoover formulation. The {scaleyz yes} and {scalexz yes} keyword/value pairs can not be used for 2D simulations. {scaleyz yes}, {scalexz yes}, and {scalexy yes} options can only be used if the 2nd dimension in the keyword is periodic, and if the tilt factor is not coupled to the barostat via keywords {tri}, {yz}, {xz}, and {xy}. These fixes can be used with dynamic groups as defined by the "group"_group.html command. Likewise they can be used with groups to which atoms are added or deleted over time, e.g. a deposition simulation. However, the conservation properties of the thermostat and barostat are defined for systems with a static set of atoms. You may observe odd behavior if the atoms in a group vary dramatically over time or the atom count becomes very small. [Related commands:] "fix nve"_fix_nve.html, "fix_modify"_fix_modify.html, "run_style"_run_style.html [Default:] The keyword defaults are tchain = 3, pchain = 3, mtk = yes, tloop = ploop = 1, nreset = 0, drag = 0.0, dilate = all, couple = none, scaleyz = scalexz = scalexy = yes if periodic in 2nd dimension and not coupled to barostat, otherwise no. :line :link(nh-Martyna) [(Martyna)] Martyna, Tobias and Klein, J Chem Phys, 101, 4177 (1994). :link(nh-Parrinello) [(Parrinello)] Parrinello and Rahman, J Appl Phys, 52, 7182 (1981). :link(nh-Tuckerman) [(Tuckerman)] Tuckerman, Alejandre, Lopez-Rendon, Jochim, and Martyna, J Phys A: Math Gen, 39, 5629 (2006). :link(nh-Shinoda) [(Shinoda)] Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004). :link(nh-Dullweber) [(Dullweber)] Dullweber, Leimkuhler and McLachlan, J Chem Phys, 107, 5840 (1997). diff --git a/src/atom.cpp b/src/atom.cpp index 1191f0f2b..7d343a080 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -1,2268 +1,2269 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include #include #include #include #include #include #include "atom.h" #include "style_atom.h" #include "atom_vec.h" #include "atom_vec_ellipsoid.h" #include "comm.h" #include "neighbor.h" #include "force.h" #include "modify.h" #include "fix.h" #include "compute.h" #include "output.h" #include "thermo.h" #include "update.h" #include "domain.h" #include "group.h" #include "input.h" #include "variable.h" #include "molecule.h" #include "atom_masks.h" #include "math_const.h" #include "memory.h" #include "error.h" #ifdef LMP_USER_INTEL #include "neigh_request.h" #endif using namespace LAMMPS_NS; using namespace MathConst; #define DELTA 1 #define DELTA_MEMSTR 1024 #define EPSILON 1.0e-6 enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED}; // several files /* ---------------------------------------------------------------------- */ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) { natoms = 0; nlocal = nghost = nmax = 0; ntypes = 0; nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0; nbonds = nangles = ndihedrals = nimpropers = 0; firstgroupname = NULL; sortfreq = 1000; nextsort = 0; userbinsize = 0.0; maxbin = maxnext = 0; binhead = NULL; next = permute = NULL; // initialize atom arrays // customize by adding new array tag = NULL; type = mask = NULL; image = NULL; x = v = f = NULL; molecule = NULL; molindex = molatom = NULL; q = NULL; mu = NULL; omega = angmom = torque = NULL; radius = rmass = NULL; ellipsoid = line = tri = body = NULL; vfrac = s0 = NULL; x0 = NULL; spin = NULL; eradius = ervel = erforce = NULL; cs = csforce = vforce = ervelforce = NULL; etag = NULL; rho = drho = e = de = cv = NULL; vest = NULL; // USER-DPD uCond = uMech = uChem = uCG = uCGnew = NULL; duChem = NULL; dpdTheta = NULL; // USER-MESO cc = cc_flux = NULL; edpd_temp = edpd_flux = edpd_cv = NULL; // USER-SMD contact_radius = NULL; smd_data_9 = NULL; smd_stress = NULL; eff_plastic_strain = NULL; eff_plastic_strain_rate = NULL; damage = NULL; // molecular info bond_per_atom = extra_bond_per_atom = 0; num_bond = NULL; bond_type = NULL; bond_atom = NULL; angle_per_atom = extra_angle_per_atom = 0; num_angle = NULL; angle_type = NULL; angle_atom1 = angle_atom2 = angle_atom3 = NULL; dihedral_per_atom = extra_dihedral_per_atom = 0; num_dihedral = NULL; dihedral_type = NULL; dihedral_atom1 = dihedral_atom2 = dihedral_atom3 = dihedral_atom4 = NULL; improper_per_atom = extra_improper_per_atom = 0; num_improper = NULL; improper_type = NULL; improper_atom1 = improper_atom2 = improper_atom3 = improper_atom4 = NULL; maxspecial = 1; nspecial = NULL; special = NULL; // user-defined molecules nmolecule = 0; molecules = NULL; // custom atom arrays nivector = ndvector = 0; ivector = NULL; dvector = NULL; iname = dname = NULL; // initialize atom style and array existence flags // customize by adding new flag sphere_flag = peri_flag = electron_flag = 0; wavepacket_flag = sph_flag = 0; molecule_flag = 0; q_flag = mu_flag = 0; omega_flag = torque_flag = angmom_flag = 0; radius_flag = rmass_flag = 0; ellipsoid_flag = line_flag = tri_flag = body_flag = 0; vfrac_flag = 0; spin_flag = eradius_flag = ervel_flag = erforce_flag = ervelforce_flag = 0; cs_flag = csforce_flag = vforce_flag = etag_flag = 0; rho_flag = e_flag = cv_flag = vest_flag = 0; dpd_flag = edpd_flag = tdpd_flag = 0; // USER-SMD smd_flag = 0; contact_radius_flag = 0; smd_data_9_flag = 0; smd_stress_flag = 0; x0_flag = 0; eff_plastic_strain_flag = 0; eff_plastic_strain_rate_flag = 0; damage_flag = 0; // Peridynamic scale factor pdscale = 1.0; // ntype-length arrays mass = NULL; mass_setflag = NULL; // callback lists & extra restart info nextra_grow = nextra_restart = nextra_border = 0; extra_grow = extra_restart = extra_border = NULL; nextra_grow_max = nextra_restart_max = nextra_border_max = 0; nextra_store = 0; extra = NULL; // default atom ID and mapping values tag_enable = 1; map_style = map_user = 0; map_tag_max = -1; map_maxarray = map_nhash = -1; max_same = 0; sametag = NULL; map_array = NULL; map_bucket = NULL; map_hash = NULL; atom_style = NULL; avec = NULL; avec_map = new AtomVecCreatorMap(); #define ATOM_CLASS #define AtomStyle(key,Class) \ (*avec_map)[#key] = &avec_creator; #include "style_atom.h" #undef AtomStyle #undef ATOM_CLASS } /* ---------------------------------------------------------------------- */ Atom::~Atom() { delete [] atom_style; delete avec; delete avec_map; delete [] firstgroupname; memory->destroy(binhead); memory->destroy(next); memory->destroy(permute); // delete atom arrays // customize by adding new array memory->destroy(tag); memory->destroy(type); memory->destroy(mask); memory->destroy(image); memory->destroy(x); memory->destroy(v); memory->destroy(f); memory->destroy(molecule); memory->destroy(molindex); memory->destroy(molatom); memory->destroy(q); memory->destroy(mu); memory->destroy(omega); memory->destroy(angmom); memory->destroy(torque); memory->destroy(radius); memory->destroy(rmass); memory->destroy(ellipsoid); memory->destroy(line); memory->destroy(tri); memory->destroy(body); memory->destroy(vfrac); memory->destroy(s0); memory->destroy(x0); memory->destroy(spin); memory->destroy(eradius); memory->destroy(ervel); memory->destroy(erforce); memory->destroy(ervelforce); memory->destroy(cs); memory->destroy(csforce); memory->destroy(vforce); memory->destroy(etag); memory->destroy(rho); memory->destroy(drho); memory->destroy(e); memory->destroy(de); memory->destroy(cv); memory->destroy(vest); memory->destroy(contact_radius); memory->destroy(smd_data_9); memory->destroy(smd_stress); memory->destroy(eff_plastic_strain); memory->destroy(eff_plastic_strain_rate); memory->destroy(damage); memory->destroy(dpdTheta); memory->destroy(uCond); memory->destroy(uMech); memory->destroy(uChem); memory->destroy(uCG); memory->destroy(uCGnew); memory->destroy(duChem); memory->destroy(cc); memory->destroy(cc_flux); memory->destroy(edpd_temp); memory->destroy(edpd_flux); memory->destroy(edpd_cv); memory->destroy(nspecial); memory->destroy(special); memory->destroy(num_bond); memory->destroy(bond_type); memory->destroy(bond_atom); memory->destroy(num_angle); memory->destroy(angle_type); memory->destroy(angle_atom1); memory->destroy(angle_atom2); memory->destroy(angle_atom3); memory->destroy(num_dihedral); memory->destroy(dihedral_type); memory->destroy(dihedral_atom1); memory->destroy(dihedral_atom2); memory->destroy(dihedral_atom3); memory->destroy(dihedral_atom4); memory->destroy(num_improper); memory->destroy(improper_type); memory->destroy(improper_atom1); memory->destroy(improper_atom2); memory->destroy(improper_atom3); memory->destroy(improper_atom4); // delete custom atom arrays for (int i = 0; i < nivector; i++) { delete [] iname[i]; memory->destroy(ivector[i]); } if (dvector != NULL) { for (int i = 0; i < ndvector; i++) { delete [] dname[i]; memory->destroy(dvector[i]); } } memory->sfree(iname); memory->sfree(dname); memory->sfree(ivector); memory->sfree(dvector); // delete user-defined molecules for (int i = 0; i < nmolecule; i++) delete molecules[i]; memory->sfree(molecules); // delete per-type arrays delete [] mass; delete [] mass_setflag; // delete extra arrays memory->destroy(extra_grow); memory->destroy(extra_restart); memory->destroy(extra_border); memory->destroy(extra); // delete mapping data structures map_delete(); } /* ---------------------------------------------------------------------- copy modify settings from old Atom class to current Atom class ------------------------------------------------------------------------- */ void Atom::settings(Atom *old) { tag_enable = old->tag_enable; map_user = old->map_user; map_style = old->map_style; sortfreq = old->sortfreq; userbinsize = old->userbinsize; if (old->firstgroupname) { int n = strlen(old->firstgroupname) + 1; firstgroupname = new char[n]; strcpy(firstgroupname,old->firstgroupname); } } /* ---------------------------------------------------------------------- create an AtomVec style called from lammps.cpp, input script, restart file, replicate ------------------------------------------------------------------------- */ void Atom::create_avec(const char *style, int narg, char **arg, int trysuffix) { delete [] atom_style; if (avec) delete avec; atom_style = NULL; avec = NULL; // unset atom style and array existence flags // may have been set by old avec // customize by adding new flag sphere_flag = peri_flag = electron_flag = 0; wavepacket_flag = sph_flag = 0; molecule_flag = 0; q_flag = mu_flag = 0; omega_flag = torque_flag = angmom_flag = 0; radius_flag = rmass_flag = 0; ellipsoid_flag = line_flag = tri_flag = body_flag = 0; vfrac_flag = 0; spin_flag = eradius_flag = ervel_flag = erforce_flag = ervelforce_flag = 0; cs_flag = csforce_flag = vforce_flag = etag_flag = 0; rho_flag = e_flag = cv_flag = vest_flag = 0; // create instance of AtomVec // use grow() to initialize atom-based arrays to length 1 // so that x[0][0] can always be referenced even if proc has no atoms int sflag; avec = new_avec(style,trysuffix,sflag); avec->store_args(narg,arg); avec->process_args(narg,arg); avec->grow(1); if (sflag) { char estyle[256]; if (sflag == 1) sprintf(estyle,"%s/%s",style,lmp->suffix); else sprintf(estyle,"%s/%s",style,lmp->suffix2); int n = strlen(estyle) + 1; atom_style = new char[n]; strcpy(atom_style,estyle); } else { int n = strlen(style) + 1; atom_style = new char[n]; strcpy(atom_style,style); } // if molecular system: // atom IDs must be defined // force atom map to be created - // map style may be reset by map_init() and its call to map_style_set() + // map style will be reset to array vs hash to by map_init() molecular = avec->molecular; if (molecular && tag_enable == 0) error->all(FLERR,"Atom IDs must be used for molecular systems"); - if (molecular) map_style = 1; + if (molecular) map_style = 3; } /* ---------------------------------------------------------------------- generate an AtomVec class, first with suffix appended ------------------------------------------------------------------------- */ AtomVec *Atom::new_avec(const char *style, int trysuffix, int &sflag) { if (trysuffix && lmp->suffix_enable) { if (lmp->suffix) { sflag = 1; char estyle[256]; sprintf(estyle,"%s/%s",style,lmp->suffix); if (avec_map->find(estyle) != avec_map->end()) { AtomVecCreator avec_creator = (*avec_map)[estyle]; return avec_creator(lmp); } } if (lmp->suffix2) { sflag = 2; char estyle[256]; sprintf(estyle,"%s/%s",style,lmp->suffix2); if (avec_map->find(estyle) != avec_map->end()) { AtomVecCreator avec_creator = (*avec_map)[estyle]; return avec_creator(lmp); } } } sflag = 0; if (avec_map->find(style) != avec_map->end()) { AtomVecCreator avec_creator = (*avec_map)[style]; return avec_creator(lmp); } error->all(FLERR,"Unknown atom style"); return NULL; } /* ---------------------------------------------------------------------- one instance per AtomVec style in style_atom.h ------------------------------------------------------------------------- */ template AtomVec *Atom::avec_creator(LAMMPS *lmp) { return new T(lmp); } /* ---------------------------------------------------------------------- */ void Atom::init() { // delete extra array since it doesn't persist past first run if (nextra_store) { memory->destroy(extra); extra = NULL; nextra_store = 0; } // check arrays that are atom type in length check_mass(FLERR); // setup of firstgroup if (firstgroupname) { firstgroup = group->find(firstgroupname); if (firstgroup < 0) error->all(FLERR,"Could not find atom_modify first group ID"); } else firstgroup = -1; // init AtomVec avec->init(); } /* ---------------------------------------------------------------------- */ void Atom::setup() { // setup bins for sorting // cannot do this in init() because uses neighbor cutoff if (sortfreq > 0) setup_sort_bins(); } /* ---------------------------------------------------------------------- return ptr to AtomVec class if matches style or to matching hybrid sub-class return NULL if no match ------------------------------------------------------------------------- */ AtomVec *Atom::style_match(const char *style) { if (strcmp(atom_style,style) == 0) return avec; else if (strcmp(atom_style,"hybrid") == 0) { AtomVecHybrid *avec_hybrid = (AtomVecHybrid *) avec; for (int i = 0; i < avec_hybrid->nstyles; i++) if (strcmp(avec_hybrid->keywords[i],style) == 0) return avec_hybrid->styles[i]; } return NULL; } /* ---------------------------------------------------------------------- modify parameters of the atom style some options can only be invoked before simulation box is defined first and sort options cannot be used together ------------------------------------------------------------------------- */ void Atom::modify_params(int narg, char **arg) { if (narg == 0) error->all(FLERR,"Illegal atom_modify command"); int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg],"id") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command"); if (domain->box_exist) error->all(FLERR, "Atom_modify id command after simulation box is defined"); if (strcmp(arg[iarg+1],"yes") == 0) tag_enable = 1; else if (strcmp(arg[iarg+1],"no") == 0) tag_enable = 0; else error->all(FLERR,"Illegal atom_modify command"); iarg += 2; } else if (strcmp(arg[iarg],"map") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command"); if (domain->box_exist) error->all(FLERR, "Atom_modify map command after simulation box is defined"); if (strcmp(arg[iarg+1],"array") == 0) map_user = 1; else if (strcmp(arg[iarg+1],"hash") == 0) map_user = 2; + else if (strcmp(arg[iarg+1],"yes") == 0) map_user = 3; else error->all(FLERR,"Illegal atom_modify command"); map_style = map_user; iarg += 2; } else if (strcmp(arg[iarg],"first") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command"); if (strcmp(arg[iarg+1],"all") == 0) { delete [] firstgroupname; firstgroupname = NULL; } else { int n = strlen(arg[iarg+1]) + 1; firstgroupname = new char[n]; strcpy(firstgroupname,arg[iarg+1]); sortfreq = 0; } iarg += 2; } else if (strcmp(arg[iarg],"sort") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal atom_modify command"); sortfreq = force->inumeric(FLERR,arg[iarg+1]); userbinsize = force->numeric(FLERR,arg[iarg+2]); if (sortfreq < 0 || userbinsize < 0.0) error->all(FLERR,"Illegal atom_modify command"); if (sortfreq >= 0 && firstgroupname) error->all(FLERR,"Atom_modify sort and first options " "cannot be used together"); iarg += 3; } else error->all(FLERR,"Illegal atom_modify command"); } } /* ---------------------------------------------------------------------- check that atom IDs are valid error if any atom ID < 0 or atom ID = MAXTAGINT if any atom ID > 0, error if any atom ID == 0 if any atom ID > 0, error if tag_enable = 0 if all atom IDs = 0, tag_enable must be 0 if max atom IDs < natoms, must be duplicates OK if max atom IDs > natoms NOTE: not fully checking that atom IDs are unique ------------------------------------------------------------------------- */ void Atom::tag_check() { tagint min = MAXTAGINT; tagint max = 0; for (int i = 0; i < nlocal; i++) { min = MIN(min,tag[i]); max = MAX(max,tag[i]); } tagint minall,maxall; MPI_Allreduce(&min,&minall,1,MPI_LMP_TAGINT,MPI_MIN,world); MPI_Allreduce(&max,&maxall,1,MPI_LMP_TAGINT,MPI_MAX,world); if (minall < 0) error->all(FLERR,"One or more Atom IDs is negative"); if (maxall >= MAXTAGINT) error->all(FLERR,"One or more atom IDs is too big"); if (maxall > 0 && minall == 0) error->all(FLERR,"One or more atom IDs is zero"); if (maxall > 0 && tag_enable == 0) error->all(FLERR,"Non-zero atom IDs with atom_modify id = no"); if (maxall == 0 && natoms && tag_enable) error->all(FLERR,"All atom IDs = 0 but atom_modify id = yes"); if (tag_enable && maxall < natoms) error->all(FLERR,"Duplicate atom IDs exist"); } /* ---------------------------------------------------------------------- add unique tags to any atoms with tag = 0 new tags are grouped by proc and start after max current tag called after creating new atoms error if new tags will exceed MAXTAGINT ------------------------------------------------------------------------- */ void Atom::tag_extend() { // maxtag_all = max tag for all atoms tagint maxtag = 0; for (int i = 0; i < nlocal; i++) maxtag = MAX(maxtag,tag[i]); tagint maxtag_all; MPI_Allreduce(&maxtag,&maxtag_all,1,MPI_LMP_TAGINT,MPI_MAX,world); // DEBUG: useful for generating 64-bit IDs even for small systems // use only when LAMMPS is compiled with BIGBIG //maxtag_all += 1000000000000; // notag = # of atoms I own with no tag (tag = 0) // notag_sum = # of total atoms on procs <= me with no tag bigint notag = 0; for (int i = 0; i < nlocal; i++) if (tag[i] == 0) notag++; bigint notag_total; MPI_Allreduce(¬ag,¬ag_total,1,MPI_LMP_BIGINT,MPI_SUM,world); if (notag_total >= MAXTAGINT) error->all(FLERR,"New atom IDs exceed maximum allowed ID"); bigint notag_sum; MPI_Scan(¬ag,¬ag_sum,1,MPI_LMP_BIGINT,MPI_SUM,world); // itag = 1st new tag that my untagged atoms should use tagint itag = maxtag_all + notag_sum - notag + 1; for (int i = 0; i < nlocal; i++) if (tag[i] == 0) tag[i] = itag++; } /* ---------------------------------------------------------------------- check that atom IDs span range from 1 to Natoms inclusive return 0 if mintag != 1 or maxtag != Natoms return 1 if OK doesn't actually check if all tag values are used ------------------------------------------------------------------------- */ int Atom::tag_consecutive() { tagint idmin = MAXTAGINT; tagint idmax = 0; for (int i = 0; i < nlocal; i++) { idmin = MIN(idmin,tag[i]); idmax = MAX(idmax,tag[i]); } tagint idminall,idmaxall; MPI_Allreduce(&idmin,&idminall,1,MPI_LMP_TAGINT,MPI_MIN,world); MPI_Allreduce(&idmax,&idmaxall,1,MPI_LMP_TAGINT,MPI_MAX,world); if (idminall != 1 || idmaxall != natoms) return 0; return 1; } /* ---------------------------------------------------------------------- count and return words in a single line make copy of line before using strtok so as not to change line trim anything from '#' onward ------------------------------------------------------------------------- */ int Atom::count_words(const char *line) { int n = strlen(line) + 1; char *copy; memory->create(copy,n,"atom:copy"); strcpy(copy,line); char *ptr; if ((ptr = strchr(copy,'#'))) *ptr = '\0'; if (strtok(copy," \t\n\r\f") == NULL) { memory->destroy(copy); return 0; } n = 1; while (strtok(NULL," \t\n\r\f")) n++; memory->destroy(copy); return n; } /* ---------------------------------------------------------------------- count and return words in a single line using provided copy buf make copy of line before using strtok so as not to change line trim anything from '#' onward ------------------------------------------------------------------------- */ int Atom::count_words(const char *line, char *copy) { strcpy(copy,line); char *ptr; if ((ptr = strchr(copy,'#'))) *ptr = '\0'; if (strtok(copy," \t\n\r\f") == NULL) { memory->destroy(copy); return 0; } int n = 1; while (strtok(NULL," \t\n\r\f")) n++; return n; } /* ---------------------------------------------------------------------- deallocate molecular topology arrays done before realloc with (possibly) new 2nd dimension set to correctly initialized per-atom values, e.g. bond_per_atom needs to be called whenever 2nd dimensions are changed and these arrays are already pre-allocated, e.g. due to grow(1) in create_avec() ------------------------------------------------------------------------- */ void Atom::deallocate_topology() { memory->destroy(atom->bond_type); memory->destroy(atom->bond_atom); atom->bond_type = NULL; atom->bond_atom = NULL; memory->destroy(atom->angle_type); memory->destroy(atom->angle_atom1); memory->destroy(atom->angle_atom2); memory->destroy(atom->angle_atom3); atom->angle_type = NULL; atom->angle_atom1 = atom->angle_atom2 = atom->angle_atom3 = NULL; memory->destroy(atom->dihedral_type); memory->destroy(atom->dihedral_atom1); memory->destroy(atom->dihedral_atom2); memory->destroy(atom->dihedral_atom3); memory->destroy(atom->dihedral_atom4); atom->dihedral_type = NULL; atom->dihedral_atom1 = atom->dihedral_atom2 = atom->dihedral_atom3 = atom->dihedral_atom4 = NULL; memory->destroy(atom->improper_type); memory->destroy(atom->improper_atom1); memory->destroy(atom->improper_atom2); memory->destroy(atom->improper_atom3); memory->destroy(atom->improper_atom4); atom->improper_type = NULL; atom->improper_atom1 = atom->improper_atom2 = atom->improper_atom3 = atom->improper_atom4 = NULL; } /* ---------------------------------------------------------------------- unpack N lines from Atom section of data file call style-specific routine to parse line ------------------------------------------------------------------------- */ void Atom::data_atoms(int n, char *buf, tagint id_offset, int type_offset, int shiftflag, double *shift) { int m,xptr,iptr; imageint imagedata; double xdata[3],lamda[3]; double *coord; char *next; next = strchr(buf,'\n'); *next = '\0'; int nwords = count_words(buf); *next = '\n'; if (nwords != avec->size_data_atom && nwords != avec->size_data_atom + 3) error->all(FLERR,"Incorrect atom format in data file"); char **values = new char*[nwords]; // set bounds for my proc // if periodic and I am lo/hi proc, adjust bounds by EPSILON // insures all data atoms will be owned even with round-off int triclinic = domain->triclinic; double epsilon[3]; if (triclinic) epsilon[0] = epsilon[1] = epsilon[2] = EPSILON; else { epsilon[0] = domain->prd[0] * EPSILON; epsilon[1] = domain->prd[1] * EPSILON; epsilon[2] = domain->prd[2] * EPSILON; } double sublo[3],subhi[3]; if (triclinic == 0) { sublo[0] = domain->sublo[0]; subhi[0] = domain->subhi[0]; sublo[1] = domain->sublo[1]; subhi[1] = domain->subhi[1]; sublo[2] = domain->sublo[2]; subhi[2] = domain->subhi[2]; } else { sublo[0] = domain->sublo_lamda[0]; subhi[0] = domain->subhi_lamda[0]; sublo[1] = domain->sublo_lamda[1]; subhi[1] = domain->subhi_lamda[1]; sublo[2] = domain->sublo_lamda[2]; subhi[2] = domain->subhi_lamda[2]; } if (comm->layout != LAYOUT_TILED) { if (domain->xperiodic) { if (comm->myloc[0] == 0) sublo[0] -= epsilon[0]; if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] += epsilon[0]; } if (domain->yperiodic) { if (comm->myloc[1] == 0) sublo[1] -= epsilon[1]; if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] += epsilon[1]; } if (domain->zperiodic) { if (comm->myloc[2] == 0) sublo[2] -= epsilon[2]; if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] += epsilon[2]; } } else { if (domain->xperiodic) { if (comm->mysplit[0][0] == 0.0) sublo[0] -= epsilon[0]; if (comm->mysplit[0][1] == 1.0) subhi[0] += epsilon[0]; } if (domain->yperiodic) { if (comm->mysplit[1][0] == 0.0) sublo[1] -= epsilon[1]; if (comm->mysplit[1][1] == 1.0) subhi[1] += epsilon[1]; } if (domain->zperiodic) { if (comm->mysplit[2][0] == 0.0) sublo[2] -= epsilon[2]; if (comm->mysplit[2][1] == 1.0) subhi[2] += epsilon[2]; } } // xptr = which word in line starts xyz coords // iptr = which word in line starts ix,iy,iz image flags xptr = avec->xcol_data - 1; int imageflag = 0; if (nwords > avec->size_data_atom) imageflag = 1; if (imageflag) iptr = nwords - 3; // loop over lines of atom data // tokenize the line into values // extract xyz coords and image flags // remap atom into simulation box // if atom is in my sub-domain, unpack its values for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); values[0] = strtok(buf," \t\n\r\f"); if (values[0] == NULL) error->all(FLERR,"Incorrect atom format in data file"); for (m = 1; m < nwords; m++) { values[m] = strtok(NULL," \t\n\r\f"); if (values[m] == NULL) error->all(FLERR,"Incorrect atom format in data file"); } if (imageflag) imagedata = ((imageint) (atoi(values[iptr]) + IMGMAX) & IMGMASK) | (((imageint) (atoi(values[iptr+1]) + IMGMAX) & IMGMASK) << IMGBITS) | (((imageint) (atoi(values[iptr+2]) + IMGMAX) & IMGMASK) << IMG2BITS); else imagedata = ((imageint) IMGMAX << IMG2BITS) | ((imageint) IMGMAX << IMGBITS) | IMGMAX; xdata[0] = atof(values[xptr]); xdata[1] = atof(values[xptr+1]); xdata[2] = atof(values[xptr+2]); if (shiftflag) { xdata[0] += shift[0]; xdata[1] += shift[1]; xdata[2] += shift[2]; } domain->remap(xdata,imagedata); if (triclinic) { domain->x2lamda(xdata,lamda); coord = lamda; } else coord = xdata; if (coord[0] >= sublo[0] && coord[0] < subhi[0] && coord[1] >= sublo[1] && coord[1] < subhi[1] && coord[2] >= sublo[2] && coord[2] < subhi[2]) { avec->data_atom(xdata,imagedata,values); if (id_offset) tag[nlocal-1] += id_offset; if (type_offset) { type[nlocal-1] += type_offset; if (type[nlocal-1] > ntypes) error->one(FLERR,"Invalid atom type in Atoms section of data file"); } } buf = next + 1; } delete [] values; } /* ---------------------------------------------------------------------- unpack N lines from Velocity section of data file check that atom IDs are > 0 and <= map_tag_max call style-specific routine to parse line ------------------------------------------------------------------------- */ void Atom::data_vels(int n, char *buf, tagint id_offset) { int j,m; tagint tagdata; char *next; next = strchr(buf,'\n'); *next = '\0'; int nwords = count_words(buf); *next = '\n'; if (nwords != avec->size_data_vel) error->all(FLERR,"Incorrect velocity format in data file"); char **values = new char*[nwords]; // loop over lines of atom velocities // tokenize the line into values // if I own atom tag, unpack its values for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); values[0] = strtok(buf," \t\n\r\f"); for (j = 1; j < nwords; j++) values[j] = strtok(NULL," \t\n\r\f"); tagdata = ATOTAGINT(values[0]) + id_offset; if (tagdata <= 0 || tagdata > map_tag_max) error->one(FLERR,"Invalid atom ID in Velocities section of data file"); if ((m = map(tagdata)) >= 0) avec->data_vel(m,&values[1]); buf = next + 1; } delete [] values; } /* ---------------------------------------------------------------------- process N bonds read into buf from data files if count is non-NULL, just count bonds per atom else store them with atoms check that atom IDs are > 0 and <= map_tag_max ------------------------------------------------------------------------- */ void Atom::data_bonds(int n, char *buf, int *count, tagint id_offset, int type_offset) { int m,tmp,itype; tagint atom1,atom2; char *next; int newton_bond = force->newton_bond; for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); *next = '\0'; sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&itype,&atom1,&atom2); if (id_offset) { atom1 += id_offset; atom2 += id_offset; } itype += type_offset; if (atom1 <= 0 || atom1 > map_tag_max || atom2 <= 0 || atom2 > map_tag_max) error->one(FLERR,"Invalid atom ID in Bonds section of data file"); if (itype <= 0 || itype > nbondtypes) error->one(FLERR,"Invalid bond type in Bonds section of data file"); if ((m = map(atom1)) >= 0) { if (count) count[m]++; else { bond_type[m][num_bond[m]] = itype; bond_atom[m][num_bond[m]] = atom2; num_bond[m]++; } } if (newton_bond == 0) { if ((m = map(atom2)) >= 0) { if (count) count[m]++; else { bond_type[m][num_bond[m]] = itype; bond_atom[m][num_bond[m]] = atom1; num_bond[m]++; } } } buf = next + 1; } } /* ---------------------------------------------------------------------- process N angles read into buf from data files if count is non-NULL, just count angles per atom else store them with atoms check that atom IDs are > 0 and <= map_tag_max ------------------------------------------------------------------------- */ void Atom::data_angles(int n, char *buf, int *count, tagint id_offset, int type_offset) { int m,tmp,itype; tagint atom1,atom2,atom3; char *next; int newton_bond = force->newton_bond; for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); *next = '\0'; sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&itype,&atom1,&atom2,&atom3); if (id_offset) { atom1 += id_offset; atom2 += id_offset; atom3 += id_offset; } itype += type_offset; if (atom1 <= 0 || atom1 > map_tag_max || atom2 <= 0 || atom2 > map_tag_max || atom3 <= 0 || atom3 > map_tag_max) error->one(FLERR,"Invalid atom ID in Angles section of data file"); if (itype <= 0 || itype > nangletypes) error->one(FLERR,"Invalid angle type in Angles section of data file"); if ((m = map(atom2)) >= 0) { if (count) count[m]++; else { angle_type[m][num_angle[m]] = itype; angle_atom1[m][num_angle[m]] = atom1; angle_atom2[m][num_angle[m]] = atom2; angle_atom3[m][num_angle[m]] = atom3; num_angle[m]++; } } if (newton_bond == 0) { if ((m = map(atom1)) >= 0) { if (count) count[m]++; else { angle_type[m][num_angle[m]] = itype; angle_atom1[m][num_angle[m]] = atom1; angle_atom2[m][num_angle[m]] = atom2; angle_atom3[m][num_angle[m]] = atom3; num_angle[m]++; } } if ((m = map(atom3)) >= 0) { if (count) count[m]++; else { angle_type[m][num_angle[m]] = itype; angle_atom1[m][num_angle[m]] = atom1; angle_atom2[m][num_angle[m]] = atom2; angle_atom3[m][num_angle[m]] = atom3; num_angle[m]++; } } } buf = next + 1; } } /* ---------------------------------------------------------------------- process N dihedrals read into buf from data files if count is non-NULL, just count diihedrals per atom else store them with atoms check that atom IDs are > 0 and <= map_tag_max ------------------------------------------------------------------------- */ void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset, int type_offset) { int m,tmp,itype; tagint atom1,atom2,atom3,atom4; char *next; int newton_bond = force->newton_bond; for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); *next = '\0'; sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&itype,&atom1,&atom2,&atom3,&atom4); if (id_offset) { atom1 += id_offset; atom2 += id_offset; atom3 += id_offset; atom4 += id_offset; } itype += type_offset; if (atom1 <= 0 || atom1 > map_tag_max || atom2 <= 0 || atom2 > map_tag_max || atom3 <= 0 || atom3 > map_tag_max || atom4 <= 0 || atom4 > map_tag_max) error->one(FLERR,"Invalid atom ID in Dihedrals section of data file"); if (itype <= 0 || itype > ndihedraltypes) error->one(FLERR, "Invalid dihedral type in Dihedrals section of data file"); if ((m = map(atom2)) >= 0) { if (count) count[m]++; else { dihedral_type[m][num_dihedral[m]] = itype; dihedral_atom1[m][num_dihedral[m]] = atom1; dihedral_atom2[m][num_dihedral[m]] = atom2; dihedral_atom3[m][num_dihedral[m]] = atom3; dihedral_atom4[m][num_dihedral[m]] = atom4; num_dihedral[m]++; } } if (newton_bond == 0) { if ((m = map(atom1)) >= 0) { if (count) count[m]++; else { dihedral_type[m][num_dihedral[m]] = itype; dihedral_atom1[m][num_dihedral[m]] = atom1; dihedral_atom2[m][num_dihedral[m]] = atom2; dihedral_atom3[m][num_dihedral[m]] = atom3; dihedral_atom4[m][num_dihedral[m]] = atom4; num_dihedral[m]++; } } if ((m = map(atom3)) >= 0) { if (count) count[m]++; else { dihedral_type[m][num_dihedral[m]] = itype; dihedral_atom1[m][num_dihedral[m]] = atom1; dihedral_atom2[m][num_dihedral[m]] = atom2; dihedral_atom3[m][num_dihedral[m]] = atom3; dihedral_atom4[m][num_dihedral[m]] = atom4; num_dihedral[m]++; } } if ((m = map(atom4)) >= 0) { if (count) count[m]++; else { dihedral_type[m][num_dihedral[m]] = itype; dihedral_atom1[m][num_dihedral[m]] = atom1; dihedral_atom2[m][num_dihedral[m]] = atom2; dihedral_atom3[m][num_dihedral[m]] = atom3; dihedral_atom4[m][num_dihedral[m]] = atom4; num_dihedral[m]++; } } } buf = next + 1; } } /* ---------------------------------------------------------------------- process N impropers read into buf from data files if count is non-NULL, just count impropers per atom else store them with atoms check that atom IDs are > 0 and <= map_tag_max ------------------------------------------------------------------------- */ void Atom::data_impropers(int n, char *buf, int *count, tagint id_offset, int type_offset) { int m,tmp,itype; tagint atom1,atom2,atom3,atom4; char *next; int newton_bond = force->newton_bond; for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); *next = '\0'; sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&itype,&atom1,&atom2,&atom3,&atom4); if (id_offset) { atom1 += id_offset; atom2 += id_offset; atom3 += id_offset; atom4 += id_offset; } itype += type_offset; if (atom1 <= 0 || atom1 > map_tag_max || atom2 <= 0 || atom2 > map_tag_max || atom3 <= 0 || atom3 > map_tag_max || atom4 <= 0 || atom4 > map_tag_max) error->one(FLERR,"Invalid atom ID in Impropers section of data file"); if (itype <= 0 || itype > nimpropertypes) error->one(FLERR, "Invalid improper type in Impropers section of data file"); if ((m = map(atom2)) >= 0) { if (count) count[m]++; else { improper_type[m][num_improper[m]] = itype; improper_atom1[m][num_improper[m]] = atom1; improper_atom2[m][num_improper[m]] = atom2; improper_atom3[m][num_improper[m]] = atom3; improper_atom4[m][num_improper[m]] = atom4; num_improper[m]++; } } if (newton_bond == 0) { if ((m = map(atom1)) >= 0) { if (count) count[m]++; else { improper_type[m][num_improper[m]] = itype; improper_atom1[m][num_improper[m]] = atom1; improper_atom2[m][num_improper[m]] = atom2; improper_atom3[m][num_improper[m]] = atom3; improper_atom4[m][num_improper[m]] = atom4; num_improper[m]++; } } if ((m = map(atom3)) >= 0) { if (count) count[m]++; else { improper_type[m][num_improper[m]] = itype; improper_atom1[m][num_improper[m]] = atom1; improper_atom2[m][num_improper[m]] = atom2; improper_atom3[m][num_improper[m]] = atom3; improper_atom4[m][num_improper[m]] = atom4; num_improper[m]++; } } if ((m = map(atom4)) >= 0) { if (count) count[m]++; else { improper_type[m][num_improper[m]] = itype; improper_atom1[m][num_improper[m]] = atom1; improper_atom2[m][num_improper[m]] = atom2; improper_atom3[m][num_improper[m]] = atom3; improper_atom4[m][num_improper[m]] = atom4; num_improper[m]++; } } } buf = next + 1; } } /* ---------------------------------------------------------------------- unpack N lines from atom-style specific bonus section of data file check that atom IDs are > 0 and <= map_tag_max call style-specific routine to parse line ------------------------------------------------------------------------- */ void Atom::data_bonus(int n, char *buf, AtomVec *avec_bonus, tagint id_offset) { int j,m,tagdata; char *next; next = strchr(buf,'\n'); *next = '\0'; int nwords = count_words(buf); *next = '\n'; if (nwords != avec_bonus->size_data_bonus) error->all(FLERR,"Incorrect bonus data format in data file"); char **values = new char*[nwords]; // loop over lines of bonus atom data // tokenize the line into values // if I own atom tag, unpack its values for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); values[0] = strtok(buf," \t\n\r\f"); for (j = 1; j < nwords; j++) values[j] = strtok(NULL," \t\n\r\f"); tagdata = ATOTAGINT(values[0]) + id_offset; if (tagdata <= 0 || tagdata > map_tag_max) error->one(FLERR,"Invalid atom ID in Bonus section of data file"); // ok to call child's data_atom_bonus() method thru parent avec_bonus, // since data_bonus() was called with child ptr, and method is virtual if ((m = map(tagdata)) >= 0) avec_bonus->data_atom_bonus(m,&values[1]); buf = next + 1; } delete [] values; } /* ---------------------------------------------------------------------- unpack N bodies from Bodies section of data file each body spans multiple lines check that atom IDs are > 0 and <= map_tag_max call style-specific routine to parse line ------------------------------------------------------------------------- */ void Atom::data_bodies(int n, char *buf, AtomVecBody *avec_body, tagint id_offset) { int j,m,nvalues,tagdata,ninteger,ndouble; int maxint = 0; int maxdouble = 0; int *ivalues = NULL; double *dvalues = NULL; // loop over lines of body data // if I own atom tag, tokenize lines into ivalues/dvalues, call data_body() // else skip values for (int i = 0; i < n; i++) { if (i == 0) tagdata = ATOTAGINT(strtok(buf," \t\n\r\f")) + id_offset; else tagdata = ATOTAGINT(strtok(NULL," \t\n\r\f")) + id_offset; if (tagdata <= 0 || tagdata > map_tag_max) error->one(FLERR,"Invalid atom ID in Bodies section of data file"); ninteger = force->inumeric(FLERR,strtok(NULL," \t\n\r\f")); ndouble = force->inumeric(FLERR,strtok(NULL," \t\n\r\f")); if ((m = map(tagdata)) >= 0) { if (ninteger > maxint) { delete [] ivalues; maxint = ninteger; ivalues = new int[maxint]; } if (ndouble > maxdouble) { delete [] dvalues; maxdouble = ndouble; dvalues = new double[maxdouble]; } for (j = 0; j < ninteger; j++) ivalues[j] = force->inumeric(FLERR,strtok(NULL," \t\n\r\f")); for (j = 0; j < ndouble; j++) dvalues[j] = force->numeric(FLERR,strtok(NULL," \t\n\r\f")); avec_body->data_body(m,ninteger,ndouble,ivalues,dvalues); } else { nvalues = ninteger + ndouble; // number of values to skip for (j = 0; j < nvalues; j++) strtok(NULL," \t\n\r\f"); } } delete [] ivalues; delete [] dvalues; } /* ---------------------------------------------------------------------- init per-atom fix/compute/variable values for newly created atoms called from create_atoms, read_data, read_dump, lib::lammps_create_atoms() fixes, computes, variables may or may not exist when called ------------------------------------------------------------------------- */ void Atom::data_fix_compute_variable(int nprev, int nnew) { for (int m = 0; m < modify->nfix; m++) { Fix *fix = modify->fix[m]; if (fix->create_attribute) for (int i = nprev; i < nnew; i++) fix->set_arrays(i); } for (int m = 0; m < modify->ncompute; m++) { Compute *compute = modify->compute[m]; if (compute->create_attribute) for (int i = nprev; i < nnew; i++) compute->set_arrays(i); } for (int i = nprev; i < nnew; i++) input->variable->set_arrays(i); } /* ---------------------------------------------------------------------- allocate arrays of length ntypes only done after ntypes is set ------------------------------------------------------------------------- */ void Atom::allocate_type_arrays() { if (avec->mass_type) { mass = new double[ntypes+1]; mass_setflag = new int[ntypes+1]; for (int itype = 1; itype <= ntypes; itype++) mass_setflag[itype] = 0; } } /* ---------------------------------------------------------------------- set a mass and flag it as set called from reading of data file type_offset may be used when reading multiple data files ------------------------------------------------------------------------- */ void Atom::set_mass(const char *file, int line, const char *str, int type_offset) { if (mass == NULL) error->all(file,line,"Cannot set mass for this atom style"); int itype; double mass_one; int n = sscanf(str,"%d %lg",&itype,&mass_one); if (n != 2) error->all(file,line,"Invalid mass line in data file"); itype += type_offset; if (itype < 1 || itype > ntypes) error->all(file,line,"Invalid type for mass set"); mass[itype] = mass_one; mass_setflag[itype] = 1; if (mass[itype] <= 0.0) error->all(file,line,"Invalid mass value"); } /* ---------------------------------------------------------------------- set a mass and flag it as set called from EAM pair routine ------------------------------------------------------------------------- */ void Atom::set_mass(const char *file, int line, int itype, double value) { if (mass == NULL) error->all(file,line,"Cannot set mass for this atom style"); if (itype < 1 || itype > ntypes) error->all(file,line,"Invalid type for mass set"); mass[itype] = value; mass_setflag[itype] = 1; if (mass[itype] <= 0.0) error->all(file,line,"Invalid mass value"); } /* ---------------------------------------------------------------------- set one or more masses and flag them as set called from reading of input script ------------------------------------------------------------------------- */ void Atom::set_mass(const char *file, int line, int narg, char **arg) { if (mass == NULL) error->all(file,line,"Cannot set mass for this atom style"); int lo,hi; force->bounds(file,line,arg[0],ntypes,lo,hi); if (lo < 1 || hi > ntypes) error->all(file,line,"Invalid type for mass set"); for (int itype = lo; itype <= hi; itype++) { mass[itype] = atof(arg[1]); mass_setflag[itype] = 1; if (mass[itype] <= 0.0) error->all(file,line,"Invalid mass value"); } } /* ---------------------------------------------------------------------- set all masses as read in from restart file ------------------------------------------------------------------------- */ void Atom::set_mass(double *values) { for (int itype = 1; itype <= ntypes; itype++) { mass[itype] = values[itype]; mass_setflag[itype] = 1; } } /* ---------------------------------------------------------------------- check that all per-atom-type masses have been set ------------------------------------------------------------------------- */ void Atom::check_mass(const char *file, int line) { if (mass == NULL) return; if (rmass_flag) return; for (int itype = 1; itype <= ntypes; itype++) if (mass_setflag[itype] == 0) error->all(file,line,"Not all per-type masses are set"); } /* ---------------------------------------------------------------------- check that radii of all particles of itype are the same return 1 if true, else return 0 also return the radius value for that type ------------------------------------------------------------------------- */ int Atom::radius_consistency(int itype, double &rad) { double value = -1.0; int flag = 0; for (int i = 0; i < nlocal; i++) { if (type[i] != itype) continue; if (value < 0.0) value = radius[i]; else if (value != radius[i]) flag = 1; } int flagall; MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); if (flagall) return 0; MPI_Allreduce(&value,&rad,1,MPI_DOUBLE,MPI_MAX,world); return 1; } /* ---------------------------------------------------------------------- check that shape of all particles of itype are the same return 1 if true, else return 0 also return the 3 shape params for itype ------------------------------------------------------------------------- */ int Atom::shape_consistency(int itype, double &shapex, double &shapey, double &shapez) { double zero[3] = {0.0, 0.0, 0.0}; double one[3] = {-1.0, -1.0, -1.0}; double *shape; AtomVecEllipsoid *avec_ellipsoid = (AtomVecEllipsoid *) style_match("ellipsoid"); AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus; int flag = 0; for (int i = 0; i < nlocal; i++) { if (type[i] != itype) continue; if (ellipsoid[i] < 0) shape = zero; else shape = bonus[ellipsoid[i]].shape; if (one[0] < 0.0) { one[0] = shape[0]; one[1] = shape[1]; one[2] = shape[2]; } else if (one[0] != shape[0] || one[1] != shape[1] || one[2] != shape[2]) flag = 1; } int flagall; MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); if (flagall) return 0; double oneall[3]; MPI_Allreduce(one,oneall,3,MPI_DOUBLE,MPI_MAX,world); shapex = oneall[0]; shapey = oneall[1]; shapez = oneall[2]; return 1; } /* ---------------------------------------------------------------------- add a new molecule template = set of molecules ------------------------------------------------------------------------- */ void Atom::add_molecule(int narg, char **arg) { if (narg < 1) error->all(FLERR,"Illegal molecule command"); if (find_molecule(arg[0]) >= 0) error->all(FLERR,"Reuse of molecule template ID"); // 1st molecule in set stores nset = # of mols, others store nset = 0 // ifile = count of molecules in set // index = argument index where next molecule starts, updated by constructor int ifile = 1; int index = 1; while (1) { molecules = (Molecule **) memory->srealloc(molecules,(nmolecule+1)*sizeof(Molecule *), "atom::molecules"); molecules[nmolecule] = new Molecule(lmp,narg,arg,index); molecules[nmolecule]->nset = 0; molecules[nmolecule-ifile+1]->nset++; nmolecule++; if (molecules[nmolecule-1]->last) break; ifile++; } } /* ---------------------------------------------------------------------- find first molecule in set with template ID return -1 if does not exist ------------------------------------------------------------------------- */ int Atom::find_molecule(char *id) { if(id == NULL) return -1; int imol; for (imol = 0; imol < nmolecule; imol++) if (strcmp(id,molecules[imol]->id) == 0) return imol; return -1; } /* ---------------------------------------------------------------------- add info to current atom ilocal from molecule template onemol and its iatom offset = atom ID preceding IDs of atoms in this molecule called by fixes and commands that add molecules ------------------------------------------------------------------------- */ void Atom::add_molecule_atom(Molecule *onemol, int iatom, int ilocal, tagint offset) { if (onemol->qflag && q_flag) q[ilocal] = onemol->q[iatom]; if (onemol->radiusflag && radius_flag) radius[ilocal] = onemol->radius[iatom]; if (onemol->rmassflag && rmass_flag) rmass[ilocal] = onemol->rmass[iatom]; else if (rmass_flag) rmass[ilocal] = 4.0*MY_PI/3.0 * radius[ilocal]*radius[ilocal]*radius[ilocal]; if (onemol->bodyflag) { body[ilocal] = 0; // as if a body read from data file onemol->avec_body->data_body(ilocal,onemol->nibody,onemol->ndbody, onemol->ibodyparams,onemol->dbodyparams); onemol->avec_body->set_quat(ilocal,onemol->quat_external); } if (molecular != 1) return; // add bond topology info // for molecular atom styles, but not atom style template if (avec->bonds_allow) { num_bond[ilocal] = onemol->num_bond[iatom]; for (int i = 0; i < num_bond[ilocal]; i++) { bond_type[ilocal][i] = onemol->bond_type[iatom][i]; bond_atom[ilocal][i] = onemol->bond_atom[iatom][i] + offset; } } if (avec->angles_allow) { num_angle[ilocal] = onemol->num_angle[iatom]; for (int i = 0; i < num_angle[ilocal]; i++) { angle_type[ilocal][i] = onemol->angle_type[iatom][i]; angle_atom1[ilocal][i] = onemol->angle_atom1[iatom][i] + offset; angle_atom2[ilocal][i] = onemol->angle_atom2[iatom][i] + offset; angle_atom3[ilocal][i] = onemol->angle_atom3[iatom][i] + offset; } } if (avec->dihedrals_allow) { num_dihedral[ilocal] = onemol->num_dihedral[iatom]; for (int i = 0; i < num_dihedral[ilocal]; i++) { dihedral_type[ilocal][i] = onemol->dihedral_type[iatom][i]; dihedral_atom1[ilocal][i] = onemol->dihedral_atom1[iatom][i] + offset; dihedral_atom2[ilocal][i] = onemol->dihedral_atom2[iatom][i] + offset; dihedral_atom3[ilocal][i] = onemol->dihedral_atom3[iatom][i] + offset; dihedral_atom4[ilocal][i] = onemol->dihedral_atom4[iatom][i] + offset; } } if (avec->impropers_allow) { num_improper[ilocal] = onemol->num_improper[iatom]; for (int i = 0; i < num_improper[ilocal]; i++) { improper_type[ilocal][i] = onemol->improper_type[iatom][i]; improper_atom1[ilocal][i] = onemol->improper_atom1[iatom][i] + offset; improper_atom2[ilocal][i] = onemol->improper_atom2[iatom][i] + offset; improper_atom3[ilocal][i] = onemol->improper_atom3[iatom][i] + offset; improper_atom4[ilocal][i] = onemol->improper_atom4[iatom][i] + offset; } } if (onemol->specialflag) { nspecial[ilocal][0] = onemol->nspecial[iatom][0]; nspecial[ilocal][1] = onemol->nspecial[iatom][1]; int n = nspecial[ilocal][2] = onemol->nspecial[iatom][2]; for (int i = 0; i < n; i++) special[ilocal][i] = onemol->special[iatom][i] + offset; } } /* ---------------------------------------------------------------------- reorder owned atoms so those in firstgroup appear first called by comm->exchange() if atom_modify first group is set only owned atoms exist at this point, no ghost atoms ------------------------------------------------------------------------- */ void Atom::first_reorder() { // insure there is one extra atom location at end of arrays for swaps if (nlocal == nmax) avec->grow(0); // loop over owned atoms // nfirst = index of first atom not in firstgroup // when find firstgroup atom out of place, swap it with atom nfirst int bitmask = group->bitmask[firstgroup]; nfirst = 0; while (nfirst < nlocal && mask[nfirst] & bitmask) nfirst++; for (int i = 0; i < nlocal; i++) { if (mask[i] & bitmask && i > nfirst) { avec->copy(i,nlocal,0); avec->copy(nfirst,i,0); avec->copy(nlocal,nfirst,0); while (nfirst < nlocal && mask[nfirst] & bitmask) nfirst++; } } } /* ---------------------------------------------------------------------- perform spatial sort of atoms within my sub-domain always called between comm->exchange() and comm->borders() don't have to worry about clearing/setting atom->map since done in comm ------------------------------------------------------------------------- */ void Atom::sort() { int i,m,n,ix,iy,iz,ibin,empty; // set next timestep for sorting to take place nextsort = (update->ntimestep/sortfreq)*sortfreq + sortfreq; // re-setup sort bins if needed if (domain->box_change) setup_sort_bins(); if (nbins == 1) return; // reallocate per-atom vectors if needed if (nlocal > maxnext) { memory->destroy(next); memory->destroy(permute); maxnext = atom->nmax; memory->create(next,maxnext,"atom:next"); memory->create(permute,maxnext,"atom:permute"); } // insure there is one extra atom location at end of arrays for swaps if (nlocal == nmax) avec->grow(0); // bin atoms in reverse order so linked list will be in forward order for (i = 0; i < nbins; i++) binhead[i] = -1; for (i = nlocal-1; i >= 0; i--) { ix = static_cast ((x[i][0]-bboxlo[0])*bininvx); iy = static_cast ((x[i][1]-bboxlo[1])*bininvy); iz = static_cast ((x[i][2]-bboxlo[2])*bininvz); ix = MAX(ix,0); iy = MAX(iy,0); iz = MAX(iz,0); ix = MIN(ix,nbinx-1); iy = MIN(iy,nbiny-1); iz = MIN(iz,nbinz-1); ibin = iz*nbiny*nbinx + iy*nbinx + ix; next[i] = binhead[ibin]; binhead[ibin] = i; } // permute = desired permutation of atoms // permute[I] = J means Ith new atom will be Jth old atom n = 0; for (m = 0; m < nbins; m++) { i = binhead[m]; while (i >= 0) { permute[n++] = i; i = next[i]; } } // current = current permutation, just reuse next vector // current[I] = J means Ith current atom is Jth old atom int *current = next; for (i = 0; i < nlocal; i++) current[i] = i; // reorder local atom list, when done, current = permute // perform "in place" using copy() to extra atom location at end of list // inner while loop processes one cycle of the permutation // copy before inner-loop moves an atom to end of atom list // copy after inner-loop moves atom at end of list back into list // empty = location in atom list that is currently empty for (i = 0; i < nlocal; i++) { if (current[i] == permute[i]) continue; avec->copy(i,nlocal,0); empty = i; while (permute[empty] != i) { avec->copy(permute[empty],empty,0); empty = current[empty] = permute[empty]; } avec->copy(nlocal,empty,0); current[empty] = permute[empty]; } // sanity check that current = permute //int flag = 0; //for (i = 0; i < nlocal; i++) // if (current[i] != permute[i]) flag = 1; //int flagall; //MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); //if (flagall) error->all(FLERR,"Atom sort did not operate correctly"); } /* ---------------------------------------------------------------------- setup bins for spatial sorting of atoms ------------------------------------------------------------------------- */ void Atom::setup_sort_bins() { // binsize: // user setting if explicitly set // default = 1/2 of neighbor cutoff // check if neighbor cutoff = 0.0 double binsize; if (userbinsize > 0.0) binsize = userbinsize; else binsize = 0.5 * neighbor->cutneighmax; if (binsize == 0.0) error->all(FLERR,"Atom sorting has bin size = 0.0"); double bininv = 1.0/binsize; // nbin xyz = local bins // bbox lo/hi = bounding box of my sub-domain if (domain->triclinic) domain->bbox(domain->sublo_lamda,domain->subhi_lamda,bboxlo,bboxhi); else { bboxlo[0] = domain->sublo[0]; bboxlo[1] = domain->sublo[1]; bboxlo[2] = domain->sublo[2]; bboxhi[0] = domain->subhi[0]; bboxhi[1] = domain->subhi[1]; bboxhi[2] = domain->subhi[2]; } nbinx = static_cast ((bboxhi[0]-bboxlo[0]) * bininv); nbiny = static_cast ((bboxhi[1]-bboxlo[1]) * bininv); nbinz = static_cast ((bboxhi[2]-bboxlo[2]) * bininv); if (domain->dimension == 2) nbinz = 1; if (nbinx == 0) nbinx = 1; if (nbiny == 0) nbiny = 1; if (nbinz == 0) nbinz = 1; bininvx = nbinx / (bboxhi[0]-bboxlo[0]); bininvy = nbiny / (bboxhi[1]-bboxlo[1]); bininvz = nbinz / (bboxhi[2]-bboxlo[2]); #ifdef LMP_USER_INTEL int intel_neigh = 0; if (neighbor->nrequest) { if (neighbor->requests[0]->intel) intel_neigh = 1; } else if (neighbor->old_nrequest) if (neighbor->old_requests[0]->intel) intel_neigh = 1; if (intel_neigh && userbinsize == 0.0) { if (neighbor->binsizeflag) bininv = 1.0/neighbor->binsize_user; double nx_low = neighbor->bboxlo[0]; double ny_low = neighbor->bboxlo[1]; double nz_low = neighbor->bboxlo[2]; double nxbbox = neighbor->bboxhi[0] - nx_low; double nybbox = neighbor->bboxhi[1] - ny_low; double nzbbox = neighbor->bboxhi[2] - nz_low; int nnbinx = static_cast (nxbbox * bininv); int nnbiny = static_cast (nybbox * bininv); int nnbinz = static_cast (nzbbox * bininv); if (domain->dimension == 2) nnbinz = 1; if (nnbinx == 0) nnbinx = 1; if (nnbiny == 0) nnbiny = 1; if (nnbinz == 0) nnbinz = 1; double binsizex = nxbbox/nnbinx; double binsizey = nybbox/nnbiny; double binsizez = nzbbox/nnbinz; bininvx = 1.0 / binsizex; bininvy = 1.0 / binsizey; bininvz = 1.0 / binsizez; int lxo = (bboxlo[0] - nx_low) * bininvx; int lyo = (bboxlo[1] - ny_low) * bininvy; int lzo = (bboxlo[2] - nz_low) * bininvz; bboxlo[0] = nx_low + static_cast(lxo) / bininvx; bboxlo[1] = ny_low + static_cast(lyo) / bininvy; bboxlo[2] = nz_low + static_cast(lzo) / bininvz; nbinx = static_cast((bboxhi[0] - bboxlo[0]) * bininvx) + 1; nbiny = static_cast((bboxhi[1] - bboxlo[1]) * bininvy) + 1; nbinz = static_cast((bboxhi[2] - bboxlo[2]) * bininvz) + 1; bboxhi[0] = bboxlo[0] + static_cast(nbinx) / bininvx; bboxhi[1] = bboxlo[1] + static_cast(nbiny) / bininvy; bboxhi[2] = bboxlo[2] + static_cast(nbinz) / bininvz; } #endif if (1.0*nbinx*nbiny*nbinz > INT_MAX) error->one(FLERR,"Too many atom sorting bins"); nbins = nbinx*nbiny*nbinz; // reallocate per-bin memory if needed if (nbins > maxbin) { memory->destroy(binhead); maxbin = nbins; memory->create(binhead,maxbin,"atom:binhead"); } } /* ---------------------------------------------------------------------- register a callback to a fix so it can manage atom-based arrays happens when fix is created flag = 0 for grow, 1 for restart, 2 for border comm ------------------------------------------------------------------------- */ void Atom::add_callback(int flag) { int ifix; // find the fix // if find NULL ptr: // it's this one, since it is being replaced and has just been deleted // at this point in re-creation // if don't find NULL ptr: // i is set to nfix = new one currently being added at end of list for (ifix = 0; ifix < modify->nfix; ifix++) if (modify->fix[ifix] == NULL) break; // add callback to lists, reallocating if necessary if (flag == 0) { if (nextra_grow == nextra_grow_max) { nextra_grow_max += DELTA; memory->grow(extra_grow,nextra_grow_max,"atom:extra_grow"); } extra_grow[nextra_grow] = ifix; nextra_grow++; } else if (flag == 1) { if (nextra_restart == nextra_restart_max) { nextra_restart_max += DELTA; memory->grow(extra_restart,nextra_restart_max,"atom:extra_restart"); } extra_restart[nextra_restart] = ifix; nextra_restart++; } else if (flag == 2) { if (nextra_border == nextra_border_max) { nextra_border_max += DELTA; memory->grow(extra_border,nextra_border_max,"atom:extra_border"); } extra_border[nextra_border] = ifix; nextra_border++; } } /* ---------------------------------------------------------------------- unregister a callback to a fix happens when fix is deleted, called by its destructor flag = 0 for grow, 1 for restart ------------------------------------------------------------------------- */ void Atom::delete_callback(const char *id, int flag) { if (id == NULL) return; int ifix; for (ifix = 0; ifix < modify->nfix; ifix++) if (strcmp(id,modify->fix[ifix]->id) == 0) break; // compact the list of callbacks if (flag == 0) { int match; for (match = 0; match < nextra_grow; match++) if (extra_grow[match] == ifix) break; for (int i = match; i < nextra_grow-1; i++) extra_grow[i] = extra_grow[i+1]; nextra_grow--; } else if (flag == 1) { int match; for (match = 0; match < nextra_restart; match++) if (extra_restart[match] == ifix) break; for (int i = match; i < nextra_restart-1; i++) extra_restart[i] = extra_restart[i+1]; nextra_restart--; } else if (flag == 2) { int match; for (match = 0; match < nextra_border; match++) if (extra_border[match] == ifix) break; for (int i = match; i < nextra_border-1; i++) extra_border[i] = extra_border[i+1]; nextra_border--; } } /* ---------------------------------------------------------------------- decrement ptrs in callback lists to fixes beyond the deleted ifix happens after fix is deleted ------------------------------------------------------------------------- */ void Atom::update_callback(int ifix) { for (int i = 0; i < nextra_grow; i++) if (extra_grow[i] > ifix) extra_grow[i]--; for (int i = 0; i < nextra_restart; i++) if (extra_restart[i] > ifix) extra_restart[i]--; for (int i = 0; i < nextra_border; i++) if (extra_border[i] > ifix) extra_border[i]--; } /* ---------------------------------------------------------------------- find custom per-atom vector with name return index if found, and flag = 0/1 for int/double return -1 if not found ------------------------------------------------------------------------- */ int Atom::find_custom(const char *name, int &flag) { if(name == NULL) return -1; for (int i = 0; i < nivector; i++) if (iname[i] && strcmp(iname[i],name) == 0) { flag = 0; return i; } for (int i = 0; i < ndvector; i++) if (dname[i] && strcmp(dname[i],name) == 0) { flag = 1; return i; } return -1; } /* ---------------------------------------------------------------------- add a custom variable with name of type flag = 0/1 for int/double assumes name does not already exist return index in ivector or dvector of its location ------------------------------------------------------------------------- */ int Atom::add_custom(const char *name, int flag) { int index; if (flag == 0) { index = nivector; nivector++; iname = (char **) memory->srealloc(iname,nivector*sizeof(char *), "atom:iname"); int n = strlen(name) + 1; iname[index] = new char[n]; strcpy(iname[index],name); ivector = (int **) memory->srealloc(ivector,nivector*sizeof(int *), "atom:ivector"); memory->create(ivector[index],nmax,"atom:ivector"); } else { index = ndvector; ndvector++; dname = (char **) memory->srealloc(dname,ndvector*sizeof(char *), "atom:dname"); int n = strlen(name) + 1; dname[index] = new char[n]; strcpy(dname[index],name); dvector = (double **) memory->srealloc(dvector,ndvector*sizeof(double *), "atom:dvector"); memory->create(dvector[index],nmax,"atom:dvector"); } return index; } /* ---------------------------------------------------------------------- remove a custom variable of type flag = 0/1 for int/double at index free memory for vector and name and set ptrs to NULL ivector/dvector and iname/dname lists never shrink ------------------------------------------------------------------------- */ void Atom::remove_custom(int flag, int index) { if (flag == 0) { memory->destroy(ivector[index]); ivector[index] = NULL; delete [] iname[index]; iname[index] = NULL; } else { memory->destroy(dvector[index]); dvector[index] = NULL; delete [] dname[index]; dname[index] = NULL; } } /* ---------------------------------------------------------------------- return a pointer to a named internal variable if don't recognize name, return NULL customize by adding names ------------------------------------------------------------------------- */ void *Atom::extract(char *name) { if (strcmp(name,"mass") == 0) return (void *) mass; if (strcmp(name,"id") == 0) return (void *) tag; if (strcmp(name,"type") == 0) return (void *) type; if (strcmp(name,"mask") == 0) return (void *) mask; if (strcmp(name,"image") == 0) return (void *) image; if (strcmp(name,"x") == 0) return (void *) x; if (strcmp(name,"v") == 0) return (void *) v; if (strcmp(name,"f") == 0) return (void *) f; if (strcmp(name,"molecule") == 0) return (void *) molecule; if (strcmp(name,"q") == 0) return (void *) q; if (strcmp(name,"mu") == 0) return (void *) mu; if (strcmp(name,"omega") == 0) return (void *) omega; if (strcmp(name,"angmom") == 0) return (void *) angmom; if (strcmp(name,"torque") == 0) return (void *) torque; if (strcmp(name,"radius") == 0) return (void *) radius; if (strcmp(name,"rmass") == 0) return (void *) rmass; if (strcmp(name,"ellipsoid") == 0) return (void *) ellipsoid; if (strcmp(name,"line") == 0) return (void *) line; if (strcmp(name,"tri") == 0) return (void *) tri; if (strcmp(name,"vfrac") == 0) return (void *) vfrac; if (strcmp(name,"s0") == 0) return (void *) s0; if (strcmp(name,"x0") == 0) return (void *) x0; if (strcmp(name,"spin") == 0) return (void *) spin; if (strcmp(name,"eradius") == 0) return (void *) eradius; if (strcmp(name,"ervel") == 0) return (void *) ervel; if (strcmp(name,"erforce") == 0) return (void *) erforce; if (strcmp(name,"ervelforce") == 0) return (void *) ervelforce; if (strcmp(name,"cs") == 0) return (void *) cs; if (strcmp(name,"csforce") == 0) return (void *) csforce; if (strcmp(name,"vforce") == 0) return (void *) vforce; if (strcmp(name,"etag") == 0) return (void *) etag; if (strcmp(name,"rho") == 0) return (void *) rho; if (strcmp(name,"drho") == 0) return (void *) drho; if (strcmp(name,"e") == 0) return (void *) e; if (strcmp(name,"de") == 0) return (void *) de; if (strcmp(name,"cv") == 0) return (void *) cv; if (strcmp(name,"vest") == 0) return (void *) vest; if (strcmp(name, "contact_radius") == 0) return (void *) contact_radius; if (strcmp(name, "smd_data_9") == 0) return (void *) smd_data_9; if (strcmp(name, "smd_stress") == 0) return (void *) smd_stress; if (strcmp(name, "eff_plastic_strain") == 0) return (void *) eff_plastic_strain; if (strcmp(name, "eff_plastic_strain_rate") == 0) return (void *) eff_plastic_strain_rate; if (strcmp(name, "damage") == 0) return (void *) damage; if (strcmp(name,"dpdTheta") == 0) return (void *) dpdTheta; if (strcmp(name,"edpd_temp") == 0) return (void *) edpd_temp; return NULL; } /* ---------------------------------------------------------------------- return # of bytes of allocated memory call to avec tallies per-atom vectors add in global to local mapping storage ------------------------------------------------------------------------- */ bigint Atom::memory_usage() { memlength = DELTA_MEMSTR; memory->create(memstr,memlength,"atom:memstr"); memstr[0] = '\0'; bigint bytes = avec->memory_usage(); memory->destroy(memstr); bytes += max_same*sizeof(int); if (map_style == 1) bytes += memory->usage(map_array,map_maxarray); else if (map_style == 2) { bytes += map_nbucket*sizeof(int); bytes += map_nhash*sizeof(HashElem); } if (maxnext) { bytes += memory->usage(next,maxnext); bytes += memory->usage(permute,maxnext); } return bytes; } /* ---------------------------------------------------------------------- accumulate per-atom vec names in memstr, padded by spaces return 1 if padded str is not already in memlist, else 0 ------------------------------------------------------------------------- */ int Atom::memcheck(const char *str) { int n = strlen(str) + 3; char *padded = new char[n]; strcpy(padded," "); strcat(padded,str); strcat(padded," "); if (strstr(memstr,padded)) { delete [] padded; return 0; } if (strlen(memstr) + n >= memlength) { memlength += DELTA_MEMSTR; memory->grow(memstr,memlength,"atom:memstr"); } strcat(memstr,padded); delete [] padded; return 1; } diff --git a/src/atom_map.cpp b/src/atom_map.cpp index bbfe014de..9d257d99d 100644 --- a/src/atom_map.cpp +++ b/src/atom_map.cpp @@ -1,380 +1,380 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include #include "atom.h" #include "comm.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; #define EXTRA 1000 /* ---------------------------------------------------------------------- allocate and initialize array or hash table for global -> local map for array option: array length = 1 to map_tag_max set entire array to -1 as initial values for hash option: map_nhash = length of hash table map_nbucket = # of hash buckets, prime larger than map_nhash * 2 so buckets will only be filled with 0 or 1 atoms on average ------------------------------------------------------------------------- */ void Atom::map_init(int check) { // check for new map style if max atomID changed (check = 1 = default) // recreate = 1 if must delete old map and create new map // recreate = 0 if can re-use old map w/out realloc and just adjust settings // map_maxarray/map_nhash initially -1, to force recreate even when no atoms int recreate = 0; if (check) recreate = map_style_set(); if (map_style == 1 && map_tag_max > map_maxarray) recreate = 1; else if (map_style == 2 && nlocal+nghost > map_nhash) recreate = 1; // if not recreating: // for array, initialize current map_tag_max values // for hash, set all buckets to empty, put all entries in free list if (!recreate) { if (map_style == 1) { for (int i = 0; i <= map_tag_max; i++) map_array[i] = -1; } else { for (int i = 0; i < map_nbucket; i++) map_bucket[i] = -1; map_nused = 0; map_free = 0; for (int i = 0; i < map_nhash; i++) map_hash[i].next = i+1; map_hash[map_nhash-1].next = -1; } // recreating: delete old map and create new one for array or hash } else { map_delete(); if (map_style == 1) { map_maxarray = map_tag_max; memory->create(map_array,map_maxarray+1,"atom:map_array"); for (int i = 0; i <= map_tag_max; i++) map_array[i] = -1; } else { // map_nhash = max # of atoms that can be hashed on this proc // set to max of ave atoms/proc or atoms I can store // multiply by 2, require at least 1000 // doubling means hash table will need to be re-init only rarely int nper = static_cast (natoms/comm->nprocs); map_nhash = MAX(nper,nmax); map_nhash *= 2; map_nhash = MAX(map_nhash,1000); // map_nbucket = prime just larger than map_nhash // next_prime() should be fast enough, // about 10% of odd integers are prime above 1M map_nbucket = next_prime(map_nhash); // set all buckets to empty // set hash to map_nhash in length // put all hash entries in free list and point them to each other map_bucket = new int[map_nbucket]; for (int i = 0; i < map_nbucket; i++) map_bucket[i] = -1; map_hash = new HashElem[map_nhash]; map_nused = 0; map_free = 0; for (int i = 0; i < map_nhash; i++) map_hash[i].next = i+1; map_hash[map_nhash-1].next = -1; } } } /* ---------------------------------------------------------------------- clear global -> local map for all of my own and ghost atoms for hash table option: global ID may not be in table if image atom was already cleared ------------------------------------------------------------------------- */ void Atom::map_clear() { if (map_style == 1) { int nall = nlocal + nghost; for (int i = 0; i < nall; i++) { sametag[i] = -1; map_array[tag[i]] = -1; } } else { int previous,ibucket,index; tagint global; int nall = nlocal + nghost; for (int i = 0; i < nall; i++) { sametag[i] = -1; // search for key // if don't find it, done previous = -1; global = tag[i]; ibucket = global % map_nbucket; index = map_bucket[ibucket]; while (index > -1) { if (map_hash[index].global == global) break; previous = index; index = map_hash[index].next; } if (index == -1) continue; // delete the hash entry and add it to free list // special logic if entry is 1st in the bucket if (previous == -1) map_bucket[ibucket] = map_hash[index].next; else map_hash[previous].next = map_hash[index].next; map_hash[index].next = map_free; map_free = index; map_nused--; } } } /* ---------------------------------------------------------------------- set global -> local map for all of my own and ghost atoms loop in reverse order so that nearby images take precedence over far ones and owned atoms take precedence over images this enables valid lookups of bond topology atoms for hash table option: if hash table too small, re-init global ID may already be in table if image atom was set ------------------------------------------------------------------------- */ void Atom::map_set() { int nall = nlocal + nghost; if (map_style == 1) { // possible reallocation of sametag must come before loop over atoms // since loop sets sametag if (nall > max_same) { max_same = nall + EXTRA; memory->destroy(sametag); memory->create(sametag,max_same,"atom:sametag"); } for (int i = nall-1; i >= 0 ; i--) { sametag[i] = map_array[tag[i]]; map_array[tag[i]] = i; } } else { // if this proc has more atoms than hash table size, call map_init() // call with 0 since max atomID in system has not changed // possible reallocation of sametag must come after map_init(), // b/c map_init() may invoke map_delete(), whacking sametag if (nall > map_nhash) map_init(0); if (nall > max_same) { max_same = nall + EXTRA; memory->destroy(sametag); memory->create(sametag,max_same,"atom:sametag"); } int previous,ibucket,index; tagint global; for (int i = nall-1; i >= 0 ; i--) { sametag[i] = map_find_hash(tag[i]); // search for key // if found it, just overwrite local value with index previous = -1; global = tag[i]; ibucket = global % map_nbucket; index = map_bucket[ibucket]; while (index > -1) { if (map_hash[index].global == global) break; previous = index; index = map_hash[index].next; } if (index > -1) { map_hash[index].local = i; continue; } // take one entry from free list // add the new global/local pair as entry at end of bucket list // special logic if this entry is 1st in bucket index = map_free; map_free = map_hash[map_free].next; if (previous == -1) map_bucket[ibucket] = index; else map_hash[previous].next = index; map_hash[index].global = global; map_hash[index].local = i; map_hash[index].next = -1; map_nused++; } } } /* ---------------------------------------------------------------------- set global to local map for one atom for hash table option: global ID may already be in table if atom was already set called by Special class ------------------------------------------------------------------------- */ void Atom::map_one(tagint global, int local) { if (map_style == 1) map_array[global] = local; else { // search for key // if found it, just overwrite local value with index int previous = -1; int ibucket = global % map_nbucket; int index = map_bucket[ibucket]; while (index > -1) { if (map_hash[index].global == global) break; previous = index; index = map_hash[index].next; } if (index > -1) { map_hash[index].local = local; return; } // take one entry from free list // add the new global/local pair as entry at end of bucket list // special logic if this entry is 1st in bucket index = map_free; map_free = map_hash[map_free].next; if (previous == -1) map_bucket[ibucket] = index; else map_hash[previous].next = index; map_hash[index].global = global; map_hash[index].local = local; map_hash[index].next = -1; map_nused++; } } /* ---------------------------------------------------------------------- set map style to array or hash based on user request or max atomID set map_tag_max = max atom ID (may be larger than natoms) called whenever map_init() called with new total # of atoms return 1 if map_style changed, else 0 ------------------------------------------------------------------------- */ int Atom::map_style_set() { if (tag_enable == 0) error->all(FLERR,"Cannot create an atom map unless atoms have IDs"); // map_tag_max = max ID of any atom that will be in new map // map_tag_max = -1 if no atoms tagint max = -1; for (int i = 0; i < nlocal; i++) max = MAX(max,tag[i]); MPI_Allreduce(&max,&map_tag_max,1,MPI_LMP_TAGINT,MPI_MAX,world); // set map_style for new map - // if user-selected, use that setting + // if user-selected to array/hash, use that setting // else if map_tag_max > 1M, use hash // else use array int map_style_old = map_style; - if (map_user) map_style = map_user; + if (map_user == 1 || map_user == 2) map_style = map_user; else if (map_tag_max > 1000000) map_style = 2; else map_style = 1; // recreate = 1 if must create new map b/c map_style changed int recreate = 0; if (map_style != map_style_old) recreate = 1; return recreate; } /* ---------------------------------------------------------------------- free the array or hash table for global to local mapping ------------------------------------------------------------------------- */ void Atom::map_delete() { memory->destroy(sametag); sametag = NULL; max_same = 0; if (map_style == 1) { memory->destroy(map_array); map_array = NULL; } else { if (map_nhash) { delete [] map_bucket; delete [] map_hash; map_bucket = NULL; map_hash = NULL; } map_nhash = 0; } } /* ---------------------------------------------------------------------- lookup global ID in hash table, return local index called by map() in atom.h ------------------------------------------------------------------------- */ int Atom::map_find_hash(tagint global) { int local = -1; int index = map_bucket[global % map_nbucket]; while (index > -1) { if (map_hash[index].global == global) { local = map_hash[index].local; break; } index = map_hash[index].next; } return local; } /* ---------------------------------------------------------------------- return next prime larger than n ------------------------------------------------------------------------- */ int Atom::next_prime(int n) { int factor; int nprime = n+1; if (nprime % 2 == 0) nprime++; int root = static_cast (sqrt(1.0*n)) + 2; while (nprime < MAXSMALLINT) { for (factor = 3; factor < root; factor++) if (nprime % factor == 0) break; if (factor == root) return nprime; nprime += 2; } return MAXSMALLINT; } diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index 04a2df91f..992049a81 100644 --- a/src/create_atoms.cpp +++ b/src/create_atoms.cpp @@ -1,838 +1,850 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include #include #include #include "create_atoms.h" #include "atom.h" #include "atom_vec.h" #include "molecule.h" #include "comm.h" #include "irregular.h" #include "modify.h" #include "force.h" #include "special.h" #include "fix.h" #include "compute.h" #include "domain.h" #include "lattice.h" #include "region.h" #include "input.h" #include "variable.h" #include "random_park.h" #include "random_mars.h" #include "math_extra.h" #include "math_const.h" #include "error.h" using namespace LAMMPS_NS; using namespace MathConst; #define BIG 1.0e30 #define EPSILON 1.0e-6 enum{BOX,REGION,SINGLE,RANDOM}; enum{ATOM,MOLECULE}; enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED}; // several files /* ---------------------------------------------------------------------- */ CreateAtoms::CreateAtoms(LAMMPS *lmp) : Pointers(lmp) {} /* ---------------------------------------------------------------------- */ void CreateAtoms::command(int narg, char **arg) { if (domain->box_exist == 0) error->all(FLERR,"Create_atoms command before simulation box is defined"); if (modify->nfix_restart_peratom) error->all(FLERR,"Cannot create_atoms after " "reading restart file with per-atom info"); // parse arguments if (narg < 2) error->all(FLERR,"Illegal create_atoms command"); ntype = force->inumeric(FLERR,arg[0]); int iarg; if (strcmp(arg[1],"box") == 0) { style = BOX; iarg = 2; } else if (strcmp(arg[1],"region") == 0) { style = REGION; if (narg < 3) error->all(FLERR,"Illegal create_atoms command"); nregion = domain->find_region(arg[2]); if (nregion == -1) error->all(FLERR, "Create_atoms region ID does not exist"); domain->regions[nregion]->init(); domain->regions[nregion]->prematch(); iarg = 3;; } else if (strcmp(arg[1],"single") == 0) { style = SINGLE; if (narg < 5) error->all(FLERR,"Illegal create_atoms command"); xone[0] = force->numeric(FLERR,arg[2]); xone[1] = force->numeric(FLERR,arg[3]); xone[2] = force->numeric(FLERR,arg[4]); iarg = 5; } else if (strcmp(arg[1],"random") == 0) { style = RANDOM; if (narg < 5) error->all(FLERR,"Illegal create_atoms command"); nrandom = force->inumeric(FLERR,arg[2]); seed = force->inumeric(FLERR,arg[3]); if (strcmp(arg[4],"NULL") == 0) nregion = -1; else { nregion = domain->find_region(arg[4]); if (nregion == -1) error->all(FLERR, "Create_atoms region ID does not exist"); domain->regions[nregion]->init(); domain->regions[nregion]->prematch(); } iarg = 5; } else error->all(FLERR,"Illegal create_atoms command"); // process optional keywords int scaleflag = 1; remapflag = 0; mode = ATOM; int molseed; varflag = 0; vstr = xstr = ystr = zstr = NULL; quatone[0] = quatone[1] = quatone[2] = 0.0; nbasis = domain->lattice->nbasis; basistype = new int[nbasis]; for (int i = 0; i < nbasis; i++) basistype[i] = ntype; while (iarg < narg) { if (strcmp(arg[iarg],"basis") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal create_atoms command"); int ibasis = force->inumeric(FLERR,arg[iarg+1]); int itype = force->inumeric(FLERR,arg[iarg+2]); if (ibasis <= 0 || ibasis > nbasis || itype <= 0 || itype > atom->ntypes) error->all(FLERR,"Invalid basis setting in create_atoms command"); basistype[ibasis-1] = itype; iarg += 3; } else if (strcmp(arg[iarg],"remap") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal create_atoms command"); if (strcmp(arg[iarg+1],"yes") == 0) remapflag = 1; else if (strcmp(arg[iarg+1],"no") == 0) remapflag = 0; else error->all(FLERR,"Illegal create_atoms command"); iarg += 2; } else if (strcmp(arg[iarg],"mol") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal create_atoms command"); int imol = atom->find_molecule(arg[iarg+1]); if (imol == -1) error->all(FLERR,"Molecule template ID for " "create_atoms does not exist"); if (atom->molecules[imol]->nset > 1 && comm->me == 0) error->warning(FLERR,"Molecule template for " "create_atoms has multiple molecules"); mode = MOLECULE; onemol = atom->molecules[imol]; molseed = force->inumeric(FLERR,arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg],"units") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal create_atoms command"); if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; else error->all(FLERR,"Illegal create_atoms command"); iarg += 2; } else if (strcmp(arg[iarg],"var") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal create_atoms command"); delete [] vstr; int n = strlen(arg[iarg+1]) + 1; vstr = new char[n]; strcpy(vstr,arg[iarg+1]); varflag = 1; iarg += 2; } else if (strcmp(arg[iarg],"set") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal create_atoms command"); if (strcmp(arg[iarg+1],"x") == 0) { delete [] xstr; int n = strlen(arg[iarg+2]) + 1; xstr = new char[n]; strcpy(xstr,arg[iarg+2]); } else if (strcmp(arg[iarg+1],"y") == 0) { delete [] ystr; int n = strlen(arg[iarg+2]) + 1; ystr = new char[n]; strcpy(ystr,arg[iarg+2]); } else if (strcmp(arg[iarg+1],"z") == 0) { delete [] zstr; int n = strlen(arg[iarg+2]) + 1; zstr = new char[n]; strcpy(zstr,arg[iarg+2]); } else error->all(FLERR,"Illegal create_atoms command"); iarg += 3; } else if (strcmp(arg[iarg],"rotate") == 0) { if (style != SINGLE) error->all(FLERR,"Cannot use create_atoms rotate unless single style"); if (iarg+5 > narg) error->all(FLERR,"Illegal create_atoms command"); double thetaone; double axisone[3]; thetaone = force->numeric(FLERR,arg[iarg+1]); axisone[0] = force->numeric(FLERR,arg[iarg+2]); axisone[1] = force->numeric(FLERR,arg[iarg+3]); axisone[2] = force->numeric(FLERR,arg[iarg+4]); if (axisone[0] == 0.0 && axisone[1] == 0.0 && axisone[2] == 0.0) error->all(FLERR,"Illegal create_atoms command"); if (domain->dimension == 2 && (axisone[0] != 0.0 || axisone[1] != 0.0)) error->all(FLERR,"Invalid create_atoms rotation vector for 2d model"); MathExtra::norm3(axisone); MathExtra::axisangle_to_quat(axisone,thetaone,quatone); iarg += 5; } else error->all(FLERR,"Illegal create_atoms command"); } // error checks if (mode == ATOM && (ntype <= 0 || ntype > atom->ntypes)) error->all(FLERR,"Invalid atom type in create_atoms command"); if (style == RANDOM) { if (nrandom < 0) error->all(FLERR,"Illegal create_atoms command"); if (seed <= 0) error->all(FLERR,"Illegal create_atoms command"); } // error check and further setup for mode = MOLECULE ranmol = NULL; if (mode == MOLECULE) { if (onemol->xflag == 0) error->all(FLERR,"Create_atoms molecule must have coordinates"); if (onemol->typeflag == 0) error->all(FLERR,"Create_atoms molecule must have atom types"); if (ntype+onemol->ntypes <= 0 || ntype+onemol->ntypes > atom->ntypes) error->all(FLERR,"Invalid atom type in create_atoms mol command"); if (onemol->tag_require && !atom->tag_enable) error->all(FLERR, "Create_atoms molecule has atom IDs, but system does not"); onemol->check_attributes(0); // create_atoms uses geoemetric center of molecule for insertion onemol->compute_center(); // molecule random number generator, different for each proc ranmol = new RanMars(lmp,molseed+comm->me); } // error check and further setup for variable test if (!vstr && (xstr || ystr || zstr)) error->all(FLERR,"Incomplete use of variables in create_atoms command"); if (vstr && (!xstr && !ystr && !zstr)) error->all(FLERR,"Incomplete use of variables in create_atoms command"); if (varflag) { vvar = input->variable->find(vstr); if (vvar < 0) error->all(FLERR,"Variable name for create_atoms does not exist"); if (!input->variable->equalstyle(vvar)) error->all(FLERR,"Variable for create_atoms is invalid style"); if (xstr) { xvar = input->variable->find(xstr); if (xvar < 0) error->all(FLERR,"Variable name for create_atoms does not exist"); if (!input->variable->internalstyle(xvar)) error->all(FLERR,"Variable for create_atoms is invalid style"); } if (ystr) { yvar = input->variable->find(ystr); if (yvar < 0) error->all(FLERR,"Variable name for create_atoms does not exist"); if (!input->variable->internalstyle(yvar)) error->all(FLERR,"Variable for create_atoms is invalid style"); } if (zstr) { zvar = input->variable->find(zstr); if (zvar < 0) error->all(FLERR,"Variable name for create_atoms does not exist"); if (!input->variable->internalstyle(zvar)) error->all(FLERR,"Variable for create_atoms is invalid style"); } } // demand non-none lattice be defined for BOX and REGION // else setup scaling for SINGLE and RANDOM // could use domain->lattice->lattice2box() to do conversion of // lattice to box, but not consistent with other uses of units=lattice // triclinic remapping occurs in add_single() if (style == BOX || style == REGION) { if (nbasis == 0) error->all(FLERR,"Cannot create atoms with undefined lattice"); } else if (scaleflag == 1) { xone[0] *= domain->lattice->xlattice; xone[1] *= domain->lattice->ylattice; xone[2] *= domain->lattice->zlattice; } // set bounds for my proc in sublo[3] & subhi[3] // if periodic and style = BOX or REGION, i.e. using lattice: // should create exactly 1 atom when 2 images are both "on" the boundary // either image may be slightly inside/outside true box due to round-off // if I am lo proc, decrement lower bound by EPSILON // this will insure lo image is created // if I am hi proc, decrement upper bound by 2.0*EPSILON // this will insure hi image is not created // thus insertion box is EPSILON smaller than true box // and is shifted away from true boundary // which is where atoms are likely to be generated triclinic = domain->triclinic; double epsilon[3]; if (triclinic) epsilon[0] = epsilon[1] = epsilon[2] = EPSILON; else { epsilon[0] = domain->prd[0] * EPSILON; epsilon[1] = domain->prd[1] * EPSILON; epsilon[2] = domain->prd[2] * EPSILON; } if (triclinic == 0) { sublo[0] = domain->sublo[0]; subhi[0] = domain->subhi[0]; sublo[1] = domain->sublo[1]; subhi[1] = domain->subhi[1]; sublo[2] = domain->sublo[2]; subhi[2] = domain->subhi[2]; } else { sublo[0] = domain->sublo_lamda[0]; subhi[0] = domain->subhi_lamda[0]; sublo[1] = domain->sublo_lamda[1]; subhi[1] = domain->subhi_lamda[1]; sublo[2] = domain->sublo_lamda[2]; subhi[2] = domain->subhi_lamda[2]; } if (style == BOX || style == REGION) { if (comm->layout != LAYOUT_TILED) { if (domain->xperiodic) { if (comm->myloc[0] == 0) sublo[0] -= epsilon[0]; if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] -= 2.0*epsilon[0]; } if (domain->yperiodic) { if (comm->myloc[1] == 0) sublo[1] -= epsilon[1]; if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] -= 2.0*epsilon[1]; } if (domain->zperiodic) { if (comm->myloc[2] == 0) sublo[2] -= epsilon[2]; if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] -= 2.0*epsilon[2]; } } else { if (domain->xperiodic) { if (comm->mysplit[0][0] == 0.0) sublo[0] -= epsilon[0]; if (comm->mysplit[0][1] == 1.0) subhi[0] -= 2.0*epsilon[0]; } if (domain->yperiodic) { if (comm->mysplit[1][0] == 0.0) sublo[1] -= epsilon[1]; if (comm->mysplit[1][1] == 1.0) subhi[1] -= 2.0*epsilon[1]; } if (domain->zperiodic) { if (comm->mysplit[2][0] == 0.0) sublo[2] -= epsilon[2]; if (comm->mysplit[2][1] == 1.0) subhi[2] -= 2.0*epsilon[2]; } } } + // Record wall time for atom creation + + MPI_Barrier(world); + double time1 = MPI_Wtime(); + // clear ghost count and any ghost bonus data internal to AtomVec // same logic as beginning of Comm::exchange() // do it now b/c creating atoms will overwrite ghost atoms atom->nghost = 0; atom->avec->clear_bonus(); // add atoms/molecules in one of 3 ways bigint natoms_previous = atom->natoms; int nlocal_previous = atom->nlocal; if (style == SINGLE) add_single(); else if (style == RANDOM) add_random(); else add_lattice(); // init per-atom fix/compute/variable values for created atoms atom->data_fix_compute_variable(nlocal_previous,atom->nlocal); // set new total # of atoms and error check bigint nblocal = atom->nlocal; MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world); if (atom->natoms < 0 || atom->natoms >= MAXBIGINT) error->all(FLERR,"Too many total atoms"); // add IDs for newly created atoms // check that atom IDs are valid if (atom->tag_enable) atom->tag_extend(); atom->tag_check(); // if global map exists, reset it // invoke map_init() b/c atom count has grown if (atom->map_style) { atom->map_init(); atom->map_set(); } // for MOLECULE mode: // molecule can mean just a mol ID or bonds/angles/etc or mol templates // set molecule IDs for created atoms if atom->molecule_flag is set // reset new molecule bond,angle,etc and special values if defined // send atoms to new owning procs via irregular comm // since not all atoms I created will be within my sub-domain // perform special list build if needed if (mode == MOLECULE) { int molecule_flag = atom->molecule_flag; int molecular = atom->molecular; tagint *molecule = atom->molecule; // molcreate = # of molecules I created int molcreate = (atom->nlocal - nlocal_previous) / onemol->natoms; // increment total bonds,angles,etc bigint nmolme = molcreate; bigint nmoltotal; MPI_Allreduce(&nmolme,&nmoltotal,1,MPI_LMP_BIGINT,MPI_SUM,world); atom->nbonds += nmoltotal * onemol->nbonds; atom->nangles += nmoltotal * onemol->nangles; atom->ndihedrals += nmoltotal * onemol->ndihedrals; atom->nimpropers += nmoltotal * onemol->nimpropers; // if atom style template // maxmol = max molecule ID across all procs, for previous atoms // moloffset = max molecule ID for all molecules owned by previous procs // including molecules existing before this creation tagint moloffset; if (molecule_flag) { tagint max = 0; for (int i = 0; i < nlocal_previous; i++) max = MAX(max,molecule[i]); tagint maxmol; MPI_Allreduce(&max,&maxmol,1,MPI_LMP_TAGINT,MPI_MAX,world); MPI_Scan(&molcreate,&moloffset,1,MPI_INT,MPI_SUM,world); moloffset = moloffset - molcreate + maxmol; } // loop over molecules I created // set their molecule ID // reset their bond,angle,etc and special values int natoms = onemol->natoms; tagint offset = 0; tagint *tag = atom->tag; int *num_bond = atom->num_bond; int *num_angle = atom->num_angle; int *num_dihedral = atom->num_dihedral; int *num_improper = atom->num_improper; tagint **bond_atom = atom->bond_atom; tagint **angle_atom1 = atom->angle_atom1; tagint **angle_atom2 = atom->angle_atom2; tagint **angle_atom3 = atom->angle_atom3; tagint **dihedral_atom1 = atom->dihedral_atom1; tagint **dihedral_atom2 = atom->dihedral_atom2; tagint **dihedral_atom3 = atom->dihedral_atom3; tagint **dihedral_atom4 = atom->dihedral_atom4; tagint **improper_atom1 = atom->improper_atom1; tagint **improper_atom2 = atom->improper_atom2; tagint **improper_atom3 = atom->improper_atom3; tagint **improper_atom4 = atom->improper_atom4; int **nspecial = atom->nspecial; tagint **special = atom->special; int ilocal = nlocal_previous; for (int i = 0; i < molcreate; i++) { if (tag) offset = tag[ilocal]-1; for (int m = 0; m < natoms; m++) { if (molecule_flag) molecule[ilocal] = moloffset + i+1; if (molecular == 2) { atom->molindex[ilocal] = 0; atom->molatom[ilocal] = m; } else if (molecular) { if (onemol->bondflag) for (int j = 0; j < num_bond[ilocal]; j++) bond_atom[ilocal][j] += offset; if (onemol->angleflag) for (int j = 0; j < num_angle[ilocal]; j++) { angle_atom1[ilocal][j] += offset; angle_atom2[ilocal][j] += offset; angle_atom3[ilocal][j] += offset; } if (onemol->dihedralflag) for (int j = 0; j < num_dihedral[ilocal]; j++) { dihedral_atom1[ilocal][j] += offset; dihedral_atom2[ilocal][j] += offset; dihedral_atom3[ilocal][j] += offset; dihedral_atom4[ilocal][j] += offset; } if (onemol->improperflag) for (int j = 0; j < num_improper[ilocal]; j++) { improper_atom1[ilocal][j] += offset; improper_atom2[ilocal][j] += offset; improper_atom3[ilocal][j] += offset; improper_atom4[ilocal][j] += offset; } if (onemol->specialflag) for (int j = 0; j < nspecial[ilocal][2]; j++) special[ilocal][j] += offset; } ilocal++; } } // perform irregular comm to migrate atoms to new owning procs double **x = atom->x; imageint *image = atom->image; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]); if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->reset_box(); Irregular *irregular = new Irregular(lmp); irregular->migrate_atoms(1); delete irregular; if (domain->triclinic) domain->lamda2x(atom->nlocal); } + MPI_Barrier(world); + double time2 = MPI_Wtime(); + // clean up delete ranmol; if (domain->lattice) delete [] basistype; delete [] vstr; delete [] xstr; delete [] ystr; delete [] zstr; // print status if (comm->me == 0) { - if (screen) + if (screen) { fprintf(screen,"Created " BIGINT_FORMAT " atoms\n", atom->natoms-natoms_previous); - if (logfile) + fprintf(screen," Time spent = %g secs\n",time2-time1); + } + if (logfile) { fprintf(logfile,"Created " BIGINT_FORMAT " atoms\n", atom->natoms-natoms_previous); + fprintf(logfile," Time spent = %g secs\n",time2-time1); + } } // for MOLECULE mode: // create special bond lists for molecular systems, // but not for atom style template // only if onemol added bonds but not special info if (mode == MOLECULE) { if (atom->molecular == 1 && onemol->bondflag && !onemol->specialflag) { Special special(lmp); special.build(); } } } /* ---------------------------------------------------------------------- add single atom with coords at xone if it's in my sub-box if triclinic, xone is in lamda coords ------------------------------------------------------------------------- */ void CreateAtoms::add_single() { // remap atom if requested if (remapflag) { imageint imagetmp = ((imageint) IMGMAX << IMG2BITS) | ((imageint) IMGMAX << IMGBITS) | IMGMAX; domain->remap(xone,imagetmp); } // if triclinic, convert to lamda coords (0-1) double lamda[3],*coord; if (triclinic) { domain->x2lamda(xone,lamda); coord = lamda; } else coord = xone; // if atom/molecule is in my subbox, create it if (coord[0] >= sublo[0] && coord[0] < subhi[0] && coord[1] >= sublo[1] && coord[1] < subhi[1] && coord[2] >= sublo[2] && coord[2] < subhi[2]) { if (mode == ATOM) atom->avec->create_atom(ntype,xone); else if (quatone[0] == 0.0 && quatone[1] == 0.0 && quatone[2] == 0.0) add_molecule(xone); else add_molecule(xone,quatone); } } /* ---------------------------------------------------------------------- add Nrandom atoms at random locations ------------------------------------------------------------------------- */ void CreateAtoms::add_random() { double xlo,ylo,zlo,xhi,yhi,zhi,zmid; double lamda[3],*coord; double *boxlo,*boxhi; // random number generator, same for all procs RanPark *random = new RanPark(lmp,seed); // bounding box for atom creation // in real units, even if triclinic // only limit bbox by region if its bboxflag is set (interior region) if (triclinic == 0) { xlo = domain->boxlo[0]; xhi = domain->boxhi[0]; ylo = domain->boxlo[1]; yhi = domain->boxhi[1]; zlo = domain->boxlo[2]; zhi = domain->boxhi[2]; zmid = zlo + 0.5*(zhi-zlo); } else { xlo = domain->boxlo_bound[0]; xhi = domain->boxhi_bound[0]; ylo = domain->boxlo_bound[1]; yhi = domain->boxhi_bound[1]; zlo = domain->boxlo_bound[2]; zhi = domain->boxhi_bound[2]; zmid = zlo + 0.5*(zhi-zlo); boxlo = domain->boxlo_lamda; boxhi = domain->boxhi_lamda; } if (nregion >= 0 && domain->regions[nregion]->bboxflag) { xlo = MAX(xlo,domain->regions[nregion]->extent_xlo); xhi = MIN(xhi,domain->regions[nregion]->extent_xhi); ylo = MAX(ylo,domain->regions[nregion]->extent_ylo); yhi = MIN(yhi,domain->regions[nregion]->extent_yhi); zlo = MAX(zlo,domain->regions[nregion]->extent_zlo); zhi = MIN(zhi,domain->regions[nregion]->extent_zhi); } // generate random positions for each new atom/molecule within bounding box // iterate until atom is within region, variable, and triclinic simulation box // if final atom position is in my subbox, create it if (xlo > xhi || ylo > yhi || zlo > zhi) error->all(FLERR,"No overlap of box and region for create_atoms"); int valid; for (int i = 0; i < nrandom; i++) { while (1) { xone[0] = xlo + random->uniform() * (xhi-xlo); xone[1] = ylo + random->uniform() * (yhi-ylo); xone[2] = zlo + random->uniform() * (zhi-zlo); if (domain->dimension == 2) xone[2] = zmid; valid = 1; if (nregion >= 0 && domain->regions[nregion]->match(xone[0],xone[1],xone[2]) == 0) valid = 0; if (varflag && vartest(xone) == 0) valid = 0; if (triclinic) { domain->x2lamda(xone,lamda); coord = lamda; if (coord[0] < boxlo[0] || coord[0] >= boxhi[0] || coord[1] < boxlo[1] || coord[1] >= boxhi[1] || coord[2] < boxlo[2] || coord[2] >= boxhi[2]) valid = 0; } else coord = xone; if (valid) break; } // if triclinic, coord is now in lamda units if (coord[0] >= sublo[0] && coord[0] < subhi[0] && coord[1] >= sublo[1] && coord[1] < subhi[1] && coord[2] >= sublo[2] && coord[2] < subhi[2]) { if (mode == ATOM) atom->avec->create_atom(ntype,xone); else add_molecule(xone); } } // clean-up delete random; } /* ---------------------------------------------------------------------- add many atoms by looping over lattice ------------------------------------------------------------------------- */ void CreateAtoms::add_lattice() { // convert 8 corners of my subdomain from box coords to lattice coords // for orthogonal, use corner pts of my subbox // for triclinic, use bounding box of my subbox // xyz min to max = bounding box around the domain corners in lattice space double bboxlo[3],bboxhi[3]; if (triclinic == 0) { bboxlo[0] = domain->sublo[0]; bboxhi[0] = domain->subhi[0]; bboxlo[1] = domain->sublo[1]; bboxhi[1] = domain->subhi[1]; bboxlo[2] = domain->sublo[2]; bboxhi[2] = domain->subhi[2]; } else domain->bbox(domain->sublo_lamda,domain->subhi_lamda,bboxlo,bboxhi); double xmin,ymin,zmin,xmax,ymax,zmax; xmin = ymin = zmin = BIG; xmax = ymax = zmax = -BIG; domain->lattice->bbox(1,bboxlo[0],bboxlo[1],bboxlo[2], xmin,ymin,zmin,xmax,ymax,zmax); domain->lattice->bbox(1,bboxhi[0],bboxlo[1],bboxlo[2], xmin,ymin,zmin,xmax,ymax,zmax); domain->lattice->bbox(1,bboxlo[0],bboxhi[1],bboxlo[2], xmin,ymin,zmin,xmax,ymax,zmax); domain->lattice->bbox(1,bboxhi[0],bboxhi[1],bboxlo[2], xmin,ymin,zmin,xmax,ymax,zmax); domain->lattice->bbox(1,bboxlo[0],bboxlo[1],bboxhi[2], xmin,ymin,zmin,xmax,ymax,zmax); domain->lattice->bbox(1,bboxhi[0],bboxlo[1],bboxhi[2], xmin,ymin,zmin,xmax,ymax,zmax); domain->lattice->bbox(1,bboxlo[0],bboxhi[1],bboxhi[2], xmin,ymin,zmin,xmax,ymax,zmax); domain->lattice->bbox(1,bboxhi[0],bboxhi[1],bboxhi[2], xmin,ymin,zmin,xmax,ymax,zmax); // ilo:ihi,jlo:jhi,klo:khi = loop bounds for lattice overlap of my subbox // overlap = any part of a unit cell (face,edge,pt) in common with my subbox // in lattice space, subbox is a tilted box // but bbox of subbox is aligned with lattice axes // so ilo:khi unit cells should completely tile bounding box // decrement lo, increment hi to avoid round-off issues in lattice->bbox(), // which can lead to missing atoms in rare cases // extra decrement of lo if min < 0, since static_cast(-1.5) = -1 int ilo,ihi,jlo,jhi,klo,khi; ilo = static_cast (xmin) - 1; jlo = static_cast (ymin) - 1; klo = static_cast (zmin) - 1; ihi = static_cast (xmax) + 1; jhi = static_cast (ymax) + 1; khi = static_cast (zmax) + 1; if (xmin < 0.0) ilo--; if (ymin < 0.0) jlo--; if (zmin < 0.0) klo--; // iterate on 3d periodic lattice of unit cells using loop bounds // iterate on nbasis atoms in each unit cell // convert lattice coords to box coords // add atom or molecule (on each basis point) if it meets all criteria double **basis = domain->lattice->basis; double x[3],lamda[3]; double *coord; int i,j,k,m; for (k = klo; k <= khi; k++) for (j = jlo; j <= jhi; j++) for (i = ilo; i <= ihi; i++) for (m = 0; m < nbasis; m++) { x[0] = i + basis[m][0]; x[1] = j + basis[m][1]; x[2] = k + basis[m][2]; // convert from lattice coords to box coords domain->lattice->lattice2box(x[0],x[1],x[2]); // if a region was specified, test if atom is in it if (style == REGION) if (!domain->regions[nregion]->match(x[0],x[1],x[2])) continue; // if variable test specified, eval variable if (varflag && vartest(x) == 0) continue; // test if atom/molecule position is in my subbox if (triclinic) { domain->x2lamda(x,lamda); coord = lamda; } else coord = x; if (coord[0] < sublo[0] || coord[0] >= subhi[0] || coord[1] < sublo[1] || coord[1] >= subhi[1] || coord[2] < sublo[2] || coord[2] >= subhi[2]) continue; // add the atom or entire molecule to my list of atoms if (mode == ATOM) atom->avec->create_atom(basistype[m],x); else add_molecule(x); } } /* ---------------------------------------------------------------------- add a randomly rotated molecule with its center at center if quat_user set, perform requested rotation ------------------------------------------------------------------------- */ void CreateAtoms::add_molecule(double *center, double *quat_user) { int n; double r[3],rotmat[3][3],quat[4],xnew[3]; if (quat_user) { quat[0] = quat_user[0]; quat[1] = quat_user[1]; quat[2] = quat_user[2]; quat[3] = quat_user[3]; } else { if (domain->dimension == 3) { r[0] = ranmol->uniform() - 0.5; r[1] = ranmol->uniform() - 0.5; r[2] = ranmol->uniform() - 0.5; } else { r[0] = r[1] = 0.0; r[2] = 1.0; } MathExtra::norm3(r); double theta = ranmol->uniform() * MY_2PI; MathExtra::axisangle_to_quat(r,theta,quat); } MathExtra::quat_to_mat(quat,rotmat); onemol->quat_external = quat; // create atoms in molecule with atom ID = 0 and mol ID = 0 // reset in caller after all molecules created by all procs // pass add_molecule_atom an offset of 0 since don't know // max tag of atoms in previous molecules at this point int natoms = onemol->natoms; for (int m = 0; m < natoms; m++) { MathExtra::matvec(rotmat,onemol->dx[m],xnew); MathExtra::add3(xnew,center,xnew); atom->avec->create_atom(ntype+onemol->type[m],xnew); n = atom->nlocal - 1; atom->add_molecule_atom(onemol,m,n,0); } } /* ---------------------------------------------------------------------- test a generated atom position against variable evaluation first set x,y,z values in internal variables ------------------------------------------------------------------------- */ int CreateAtoms::vartest(double *x) { if (xstr) input->variable->internal_set(xvar,x[0]); if (ystr) input->variable->internal_set(yvar,x[1]); if (zstr) input->variable->internal_set(zvar,x[2]); double value = input->variable->compute_equal(vvar); if (value == 0.0) return 0; return 1; } diff --git a/src/output.cpp b/src/output.cpp index ce7fcb7cc..ce593ec6a 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -1,857 +1,857 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include #include #include #include "output.h" #include "style_dump.h" #include "atom.h" #include "neighbor.h" #include "input.h" #include "variable.h" #include "comm.h" #include "update.h" #include "group.h" #include "domain.h" #include "thermo.h" #include "modify.h" #include "compute.h" #include "force.h" #include "dump.h" #include "write_restart.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; #define DELTA 1 /* ---------------------------------------------------------------------- initialize all output ------------------------------------------------------------------------- */ Output::Output(LAMMPS *lmp) : Pointers(lmp) { // create default computes for temp,pressure,pe char **newarg = new char*[4]; newarg[0] = (char *) "thermo_temp"; newarg[1] = (char *) "all"; newarg[2] = (char *) "temp"; modify->add_compute(3,newarg); newarg[0] = (char *) "thermo_press"; newarg[1] = (char *) "all"; newarg[2] = (char *) "pressure"; newarg[3] = (char *) "thermo_temp"; modify->add_compute(4,newarg); newarg[0] = (char *) "thermo_pe"; newarg[1] = (char *) "all"; newarg[2] = (char *) "pe"; modify->add_compute(3,newarg); delete [] newarg; // create default Thermo class newarg = new char*[1]; newarg[0] = (char *) "one"; thermo = new Thermo(lmp,1,newarg); delete [] newarg; thermo_every = 0; var_thermo = NULL; ndump = 0; max_dump = 0; every_dump = NULL; next_dump = NULL; last_dump = NULL; var_dump = NULL; ivar_dump = NULL; dump = NULL; restart_flag = restart_flag_single = restart_flag_double = 0; restart_every_single = restart_every_double = 0; last_restart = -1; restart1 = restart2a = restart2b = NULL; var_restart_single = var_restart_double = NULL; restart = NULL; dump_map = new DumpCreatorMap(); #define DUMP_CLASS #define DumpStyle(key,Class) \ (*dump_map)[#key] = &dump_creator; #include "style_dump.h" #undef DumpStyle #undef DUMP_CLASS } /* ---------------------------------------------------------------------- free all memory ------------------------------------------------------------------------- */ Output::~Output() { if (thermo) delete thermo; delete [] var_thermo; memory->destroy(every_dump); memory->destroy(next_dump); memory->destroy(last_dump); for (int i = 0; i < ndump; i++) delete [] var_dump[i]; memory->sfree(var_dump); memory->destroy(ivar_dump); for (int i = 0; i < ndump; i++) delete dump[i]; memory->sfree(dump); delete [] restart1; delete [] restart2a; delete [] restart2b; delete [] var_restart_single; delete [] var_restart_double; delete restart; delete dump_map; } /* ---------------------------------------------------------------------- */ void Output::init() { thermo->init(); if (var_thermo) { ivar_thermo = input->variable->find(var_thermo); if (ivar_thermo < 0) error->all(FLERR,"Variable name for thermo every does not exist"); if (!input->variable->equalstyle(ivar_thermo)) error->all(FLERR,"Variable for thermo every is invalid style"); } for (int i = 0; i < ndump; i++) dump[i]->init(); for (int i = 0; i < ndump; i++) if (every_dump[i] == 0) { ivar_dump[i] = input->variable->find(var_dump[i]); if (ivar_dump[i] < 0) error->all(FLERR,"Variable name for dump every does not exist"); if (!input->variable->equalstyle(ivar_dump[i])) error->all(FLERR,"Variable for dump every is invalid style"); } if (restart_flag_single && restart_every_single == 0) { ivar_restart_single = input->variable->find(var_restart_single); if (ivar_restart_single < 0) error->all(FLERR,"Variable name for restart does not exist"); if (!input->variable->equalstyle(ivar_restart_single)) error->all(FLERR,"Variable for restart is invalid style"); } if (restart_flag_double && restart_every_double == 0) { ivar_restart_double = input->variable->find(var_restart_double); if (ivar_restart_double < 0) error->all(FLERR,"Variable name for restart does not exist"); if (!input->variable->equalstyle(ivar_restart_double)) error->all(FLERR,"Variable for restart is invalid style"); } } /* ---------------------------------------------------------------------- perform output for setup of run/min do dump first, so memory_usage will include dump allocation do thermo last, so will print after memory_usage memflag = 0/1 for printing out memory usage ------------------------------------------------------------------------- */ void Output::setup(int memflag) { bigint ntimestep = update->ntimestep; // perform dump at start of run only if: // current timestep is multiple of every and last dump not >= this step // this is first run after dump created and firstflag is set // note that variable freq will not write unless triggered by firstflag // set next_dump to multiple of every or variable value // set next_dump_any to smallest next_dump // wrap dumps that invoke computes and variable eval with clear/add // if dump not written now, use addstep_compute_all() since don't know // what computes the dump write would invoke // if no dumps, set next_dump_any to last+1 so will not influence next int writeflag; if (ndump && update->restrict_output == 0) { for (int idump = 0; idump < ndump; idump++) { if (dump[idump]->clearstep || every_dump[idump] == 0) modify->clearstep_compute(); writeflag = 0; if (every_dump[idump] && ntimestep % every_dump[idump] == 0 && last_dump[idump] != ntimestep) writeflag = 1; if (last_dump[idump] < 0 && dump[idump]->first_flag == 1) writeflag = 1; if (writeflag) { dump[idump]->write(); last_dump[idump] = ntimestep; } if (every_dump[idump]) next_dump[idump] = (ntimestep/every_dump[idump])*every_dump[idump] + every_dump[idump]; else { bigint nextdump = static_cast (input->variable->compute_equal(ivar_dump[idump])); if (nextdump <= ntimestep) error->all(FLERR,"Dump every variable returned a bad timestep"); next_dump[idump] = nextdump; } if (dump[idump]->clearstep || every_dump[idump] == 0) { if (writeflag) modify->addstep_compute(next_dump[idump]); else modify->addstep_compute_all(next_dump[idump]); } if (idump) next_dump_any = MIN(next_dump_any,next_dump[idump]); else next_dump_any = next_dump[0]; } } else next_dump_any = update->laststep + 1; // do not write restart files at start of run // set next_restart values to multiple of every or variable value // wrap variable eval with clear/add // if no restarts, set next_restart to last+1 so will not influence next if (restart_flag && update->restrict_output == 0) { if (restart_flag_single) { if (restart_every_single) next_restart_single = (ntimestep/restart_every_single)*restart_every_single + restart_every_single; else { bigint nextrestart = static_cast (input->variable->compute_equal(ivar_restart_single)); if (nextrestart <= ntimestep) error->all(FLERR,"Restart variable returned a bad timestep"); next_restart_single = nextrestart; } } else next_restart_single = update->laststep + 1; if (restart_flag_double) { if (restart_every_double) next_restart_double = (ntimestep/restart_every_double)*restart_every_double + restart_every_double; else { bigint nextrestart = static_cast (input->variable->compute_equal(ivar_restart_double)); if (nextrestart <= ntimestep) error->all(FLERR,"Restart variable returned a bad timestep"); next_restart_double = nextrestart; } } else next_restart_double = update->laststep + 1; next_restart = MIN(next_restart_single,next_restart_double); } else next_restart = update->laststep + 1; // print memory usage unless being called between multiple runs if (memflag) memory_usage(); // set next_thermo to multiple of every or variable eval if var defined // insure thermo output on last step of run // thermo may invoke computes so wrap with clear/add modify->clearstep_compute(); thermo->header(); thermo->compute(0); last_thermo = ntimestep; if (var_thermo) { next_thermo = static_cast (input->variable->compute_equal(ivar_thermo)); if (next_thermo <= ntimestep) error->all(FLERR,"Thermo every variable returned a bad timestep"); } else if (thermo_every) { next_thermo = (ntimestep/thermo_every)*thermo_every + thermo_every; next_thermo = MIN(next_thermo,update->laststep); } else next_thermo = update->laststep; modify->addstep_compute(next_thermo); // next = next timestep any output will be done next = MIN(next_dump_any,next_restart); next = MIN(next,next_thermo); } /* ---------------------------------------------------------------------- perform all output for this timestep only perform output if next matches current step and last output doesn't do dump/restart before thermo so thermo CPU time will include them ------------------------------------------------------------------------- */ void Output::write(bigint ntimestep) { // next_dump does not force output on last step of run // wrap dumps that invoke computes or eval of variable with clear/add if (next_dump_any == ntimestep) { for (int idump = 0; idump < ndump; idump++) { if (next_dump[idump] == ntimestep) { if (dump[idump]->clearstep || every_dump[idump] == 0) modify->clearstep_compute(); if (last_dump[idump] != ntimestep) { dump[idump]->write(); last_dump[idump] = ntimestep; } if (every_dump[idump]) next_dump[idump] += every_dump[idump]; else { bigint nextdump = static_cast (input->variable->compute_equal(ivar_dump[idump])); if (nextdump <= ntimestep) error->all(FLERR,"Dump every variable returned a bad timestep"); next_dump[idump] = nextdump; } if (dump[idump]->clearstep || every_dump[idump] == 0) modify->addstep_compute(next_dump[idump]); } if (idump) next_dump_any = MIN(next_dump_any,next_dump[idump]); else next_dump_any = next_dump[0]; } } // next_restart does not force output on last step of run // for toggle = 0, replace "*" with current timestep in restart filename // eval of variable may invoke computes so wrap with clear/add if (next_restart == ntimestep) { if (next_restart_single == ntimestep) { char *file = new char[strlen(restart1) + 16]; char *ptr = strchr(restart1,'*'); *ptr = '\0'; sprintf(file,"%s" BIGINT_FORMAT "%s",restart1,ntimestep,ptr+1); *ptr = '*'; if (last_restart != ntimestep) restart->write(file); delete [] file; if (restart_every_single) next_restart_single += restart_every_single; else { modify->clearstep_compute(); bigint nextrestart = static_cast (input->variable->compute_equal(ivar_restart_single)); if (nextrestart <= ntimestep) error->all(FLERR,"Restart variable returned a bad timestep"); next_restart_single = nextrestart; modify->addstep_compute(next_restart_single); } } if (next_restart_double == ntimestep) { if (last_restart != ntimestep) { if (restart_toggle == 0) { restart->write(restart2a); restart_toggle = 1; } else { restart->write(restart2b); restart_toggle = 0; } } if (restart_every_double) next_restart_double += restart_every_double; else { modify->clearstep_compute(); bigint nextrestart = static_cast (input->variable->compute_equal(ivar_restart_double)); if (nextrestart <= ntimestep) error->all(FLERR,"Restart variable returned a bad timestep"); next_restart_double = nextrestart; modify->addstep_compute(next_restart_double); } } last_restart = ntimestep; next_restart = MIN(next_restart_single,next_restart_double); } // insure next_thermo forces output on last step of run // thermo may invoke computes so wrap with clear/add if (next_thermo == ntimestep) { modify->clearstep_compute(); if (last_thermo != ntimestep) thermo->compute(1); last_thermo = ntimestep; if (var_thermo) { next_thermo = static_cast (input->variable->compute_equal(ivar_thermo)); if (next_thermo <= ntimestep) error->all(FLERR,"Thermo every variable returned a bad timestep"); } else if (thermo_every) next_thermo += thermo_every; else next_thermo = update->laststep; next_thermo = MIN(next_thermo,update->laststep); modify->addstep_compute(next_thermo); } // next = next timestep any output will be done next = MIN(next_dump_any,next_restart); next = MIN(next,next_thermo); } /* ---------------------------------------------------------------------- force a snapshot to be written for all dumps called from PRD and TAD ------------------------------------------------------------------------- */ void Output::write_dump(bigint ntimestep) { for (int idump = 0; idump < ndump; idump++) { dump[idump]->write(); last_dump[idump] = ntimestep; } } /* ---------------------------------------------------------------------- force restart file(s) to be written called from PRD and TAD ------------------------------------------------------------------------- */ void Output::write_restart(bigint ntimestep) { if (restart_flag_single) { char *file = new char[strlen(restart1) + 16]; char *ptr = strchr(restart1,'*'); *ptr = '\0'; sprintf(file,"%s" BIGINT_FORMAT "%s",restart1,ntimestep,ptr+1); *ptr = '*'; restart->write(file); delete [] file; } if (restart_flag_double) { if (restart_toggle == 0) { restart->write(restart2a); restart_toggle = 1; } else { restart->write(restart2b); restart_toggle = 0; } } last_restart = ntimestep; } /* ---------------------------------------------------------------------- timestep is being changed, called by update->reset_timestep() reset next timestep values for dumps, restart, thermo output reset to smallest value >= new timestep if next timestep set by variable evaluation, eval for ntimestep-1, so current ntimestep can be returned if needed no guarantee that variable can be evaluated for ntimestep-1 if it depends on computes, but live with that rare case for now ------------------------------------------------------------------------- */ void Output::reset_timestep(bigint ntimestep) { next_dump_any = MAXBIGINT; for (int idump = 0; idump < ndump; idump++) { if (every_dump[idump]) { next_dump[idump] = (ntimestep/every_dump[idump])*every_dump[idump]; if (next_dump[idump] < ntimestep) next_dump[idump] += every_dump[idump]; } else { // ivar_dump may not be initialized if (ivar_dump[idump] < 0) { ivar_dump[idump] = input->variable->find(var_dump[idump]); if (ivar_dump[idump] < 0) error->all(FLERR,"Variable name for dump every does not exist"); if (!input->variable->equalstyle(ivar_dump[idump])) error->all(FLERR,"Variable for dump every is invalid style"); } modify->clearstep_compute(); update->ntimestep--; bigint nextdump = static_cast (input->variable->compute_equal(ivar_dump[idump])); if (nextdump < ntimestep) error->all(FLERR,"Dump every variable returned a bad timestep"); update->ntimestep++; next_dump[idump] = nextdump; modify->addstep_compute(next_dump[idump]); } next_dump_any = MIN(next_dump_any,next_dump[idump]); } if (restart_flag_single) { if (restart_every_single) { next_restart_single = (ntimestep/restart_every_single)*restart_every_single; if (next_restart_single < ntimestep) next_restart_single += restart_every_single; } else { modify->clearstep_compute(); update->ntimestep--; bigint nextrestart = static_cast (input->variable->compute_equal(ivar_restart_single)); if (nextrestart < ntimestep) error->all(FLERR,"Restart variable returned a bad timestep"); update->ntimestep++; next_restart_single = nextrestart; modify->addstep_compute(next_restart_single); } } else next_restart_single = update->laststep + 1; if (restart_flag_double) { if (restart_every_double) { next_restart_double = (ntimestep/restart_every_double)*restart_every_double; if (next_restart_double < ntimestep) next_restart_double += restart_every_double; } else { modify->clearstep_compute(); update->ntimestep--; bigint nextrestart = static_cast (input->variable->compute_equal(ivar_restart_double)); if (nextrestart < ntimestep) error->all(FLERR,"Restart variable returned a bad timestep"); update->ntimestep++; next_restart_double = nextrestart; modify->addstep_compute(next_restart_double); } } else next_restart_double = update->laststep + 1; next_restart = MIN(next_restart_single,next_restart_double); if (var_thermo) { modify->clearstep_compute(); update->ntimestep--; next_thermo = static_cast (input->variable->compute_equal(ivar_thermo)); if (next_thermo < ntimestep) error->all(FLERR,"Thermo_modify every variable returned a bad timestep"); update->ntimestep++; next_thermo = MIN(next_thermo,update->laststep); modify->addstep_compute(next_thermo); } else if (thermo_every) { next_thermo = (ntimestep/thermo_every)*thermo_every; if (next_thermo < ntimestep) next_thermo += thermo_every; next_thermo = MIN(next_thermo,update->laststep); } else next_thermo = update->laststep; next = MIN(next_dump_any,next_restart); next = MIN(next,next_thermo); } /* ---------------------------------------------------------------------- add a Dump to list of Dumps ------------------------------------------------------------------------- */ void Output::add_dump(int narg, char **arg) { if (narg < 5) error->all(FLERR,"Illegal dump command"); // error checks for (int idump = 0; idump < ndump; idump++) if (strcmp(arg[0],dump[idump]->id) == 0) error->all(FLERR,"Reuse of dump ID"); int igroup = group->find(arg[1]); if (igroup == -1) error->all(FLERR,"Could not find dump group ID"); if (force->inumeric(FLERR,arg[3]) <= 0) error->all(FLERR,"Invalid dump frequency"); // extend Dump list if necessary if (ndump == max_dump) { max_dump += DELTA; dump = (Dump **) memory->srealloc(dump,max_dump*sizeof(Dump *),"output:dump"); memory->grow(every_dump,max_dump,"output:every_dump"); memory->grow(next_dump,max_dump,"output:next_dump"); memory->grow(last_dump,max_dump,"output:last_dump"); var_dump = (char **) memory->srealloc(var_dump,max_dump*sizeof(char *),"output:var_dump"); memory->grow(ivar_dump,max_dump,"output:ivar_dump"); } // initialize per-dump data to suitable default values every_dump[ndump] = 0; last_dump[ndump] = -1; var_dump[ndump] = NULL; ivar_dump[ndump] = -1; // create the Dump if (dump_map->find(arg[2]) != dump_map->end()) { DumpCreator dump_creator = (*dump_map)[arg[2]]; dump[ndump] = dump_creator(lmp, narg, arg); } else error->all(FLERR,"Unknown dump style"); every_dump[ndump] = force->inumeric(FLERR,arg[3]); if (every_dump[ndump] <= 0) error->all(FLERR,"Illegal dump command"); last_dump[ndump] = -1; var_dump[ndump] = NULL; ndump++; } /* ---------------------------------------------------------------------- one instance per dump style in style_dump.h ------------------------------------------------------------------------- */ template Dump *Output::dump_creator(LAMMPS *lmp, int narg, char ** arg) { return new T(lmp, narg, arg); } /* ---------------------------------------------------------------------- modify parameters of a Dump ------------------------------------------------------------------------- */ void Output::modify_dump(int narg, char **arg) { if (narg < 1) error->all(FLERR,"Illegal dump_modify command"); // find which dump it is int idump; for (idump = 0; idump < ndump; idump++) if (strcmp(arg[0],dump[idump]->id) == 0) break; if (idump == ndump) error->all(FLERR,"Cound not find dump_modify ID"); dump[idump]->modify_params(narg-1,&arg[1]); } /* ---------------------------------------------------------------------- delete a Dump from list of Dumps ------------------------------------------------------------------------- */ void Output::delete_dump(char *id) { // find which dump it is and delete it int idump; for (idump = 0; idump < ndump; idump++) if (strcmp(id,dump[idump]->id) == 0) break; if (idump == ndump) error->all(FLERR,"Could not find undump ID"); delete dump[idump]; delete [] var_dump[idump]; // move other dumps down in list one slot for (int i = idump+1; i < ndump; i++) { dump[i-1] = dump[i]; every_dump[i-1] = every_dump[i]; next_dump[i-1] = next_dump[i]; last_dump[i-1] = last_dump[i]; var_dump[i-1] = var_dump[i]; ivar_dump[i-1] = ivar_dump[i]; } ndump--; } /* ---------------------------------------------------------------------- find a dump by ID return index of dump or -1 if not found ------------------------------------------------------------------------- */ int Output::find_dump(const char *id) { if (id == NULL) return -1; int idump; for (idump = 0; idump < ndump; idump++) if (strcmp(id,dump[idump]->id) == 0) break; if (idump == ndump) return -1; return idump; } /* ---------------------------------------------------------------------- set thermo output frequency from input script ------------------------------------------------------------------------- */ void Output::set_thermo(int narg, char **arg) { if (narg != 1) error->all(FLERR,"Illegal thermo command"); if (strstr(arg[0],"v_") == arg[0]) { delete [] var_thermo; int n = strlen(&arg[0][2]) + 1; var_thermo = new char[n]; strcpy(var_thermo,&arg[0][2]); } else { thermo_every = force->inumeric(FLERR,arg[0]); if (thermo_every < 0) error->all(FLERR,"Illegal thermo command"); } } /* ---------------------------------------------------------------------- new Thermo style ------------------------------------------------------------------------- */ void Output::create_thermo(int narg, char **arg) { if (narg < 1) error->all(FLERR,"Illegal thermo_style command"); // don't allow this so that dipole style can safely allocate inertia vector if (domain->box_exist == 0) error->all(FLERR,"Thermo_style command before simulation box is defined"); // warn if previous thermo had been modified via thermo_modify command if (thermo->modified && comm->me == 0) error->warning(FLERR,"New thermo_style command, " "previous thermo_modify settings will be lost"); // set thermo = NULL in case new Thermo throws an error delete thermo; thermo = NULL; thermo = new Thermo(lmp,narg,arg); } /* ---------------------------------------------------------------------- setup restart capability for single or double output files if only one filename and it contains no "*", then append ".*" ------------------------------------------------------------------------- */ void Output::create_restart(int narg, char **arg) { if (narg < 1) error->all(FLERR,"Illegal restart command"); int every = 0; int varflag = 0; if (strstr(arg[0],"v_") == arg[0]) varflag = 1; else every = force->inumeric(FLERR,arg[0]); if (!varflag && every == 0) { if (narg != 1) error->all(FLERR,"Illegal restart command"); restart_flag = restart_flag_single = restart_flag_double = 0; last_restart = -1; delete restart; restart = NULL; delete [] restart1; delete [] restart2a; delete [] restart2b; restart1 = restart2a = restart2b = NULL; delete [] var_restart_single; delete [] var_restart_double; var_restart_single = var_restart_double = NULL; return; } if (narg < 2) error->all(FLERR,"Illegal restart command"); int nfile = 0; if (narg % 2 == 0) nfile = 1; else nfile = 2; if (nfile == 1) { restart_flag = restart_flag_single = 1; if (varflag) { delete [] var_restart_single; int n = strlen(&arg[0][2]) + 1; var_restart_single = new char[n]; strcpy(var_restart_single,&arg[0][2]); restart_every_single = 0; } else restart_every_single = every; int n = strlen(arg[1]) + 3; delete [] restart1; restart1 = new char[n]; strcpy(restart1,arg[1]); if (strchr(restart1,'*') == NULL) strcat(restart1,".*"); } if (nfile == 2) { restart_flag = restart_flag_double = 1; if (varflag) { delete [] var_restart_double; int n = strlen(&arg[0][2]) + 1; var_restart_double = new char[n]; strcpy(var_restart_double,&arg[0][2]); restart_every_double = 0; } else restart_every_double = every; delete [] restart2a; delete [] restart2b; restart_toggle = 0; int n = strlen(arg[1]) + 3; restart2a = new char[n]; strcpy(restart2a,arg[1]); n = strlen(arg[2]) + 1; restart2b = new char[n]; strcpy(restart2b,arg[2]); } // check for multiproc output and an MPI-IO filename // if 2 filenames, must be consistent int multiproc; if (strchr(arg[1],'%')) multiproc = comm->nprocs; else multiproc = 0; if (nfile == 2) { if (multiproc && !strchr(arg[2],'%')) error->all(FLERR,"Both restart files must use % or neither"); if (!multiproc && strchr(arg[2],'%')) error->all(FLERR,"Both restart files must use % or neither"); } int mpiioflag; if (strstr(arg[1],".mpi")) mpiioflag = 1; else mpiioflag = 0; if (nfile == 2) { if (mpiioflag && !strstr(arg[2],".mpi")) error->all(FLERR,"Both restart files must use MPI-IO or neither"); if (!mpiioflag && strstr(arg[2],".mpi")) error->all(FLERR,"Both restart files must use MPI-IO or neither"); } // setup output style and process optional args delete restart; restart = new WriteRestart(lmp); int iarg = nfile+1; restart->multiproc_options(multiproc,mpiioflag,narg-iarg,&arg[iarg]); } /* ---------------------------------------------------------------------- sum and print memory usage result is only memory on proc 0, not averaged across procs ------------------------------------------------------------------------- */ + void Output::memory_usage() { - bigint bytes = 0; bytes += atom->memory_usage(); bytes += neighbor->memory_usage(); bytes += comm->memory_usage(); bytes += update->memory_usage(); bytes += force->memory_usage(); bytes += modify->memory_usage(); for (int i = 0; i < ndump; i++) bytes += dump[i]->memory_usage(); double mbytes = bytes/1024.0/1024.0; double mbavg,mbmin,mbmax; MPI_Reduce(&mbytes,&mbavg,1,MPI_DOUBLE,MPI_SUM,0,world); MPI_Reduce(&mbytes,&mbmin,1,MPI_DOUBLE,MPI_MIN,0,world); MPI_Reduce(&mbytes,&mbmax,1,MPI_DOUBLE,MPI_MAX,0,world); mbavg /= comm->nprocs; if (comm->me == 0) { if (screen) fprintf(screen,"Per MPI rank memory allocation (min/avg/max) = " "%.4g | %.4g | %.4g Mbytes\n",mbmin,mbavg,mbmax); if (logfile) fprintf(logfile,"Per MPI rank memory allocation (min/avg/max) = " "%.4g | %.4g | %.4g Mbytes\n",mbmin,mbavg,mbmax); } } diff --git a/src/replicate.cpp b/src/replicate.cpp index e2ed718f6..f3d196416 100644 --- a/src/replicate.cpp +++ b/src/replicate.cpp @@ -1,427 +1,444 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include #include #include "replicate.h" #include "atom.h" #include "atom_vec.h" #include "atom_vec_hybrid.h" #include "force.h" #include "domain.h" #include "comm.h" #include "special.h" #include "accelerator_kokkos.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; #define LB_FACTOR 1.1 #define EPSILON 1.0e-6 enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED}; // several files /* ---------------------------------------------------------------------- */ Replicate::Replicate(LAMMPS *lmp) : Pointers(lmp) {} /* ---------------------------------------------------------------------- */ void Replicate::command(int narg, char **arg) { int i,j,m,n; if (domain->box_exist == 0) error->all(FLERR,"Replicate command before simulation box is defined"); if (narg != 3) error->all(FLERR,"Illegal replicate command"); int me = comm->me; int nprocs = comm->nprocs; if (me == 0 && screen) fprintf(screen,"Replicating atoms ...\n"); // nrep = total # of replications int nx = force->inumeric(FLERR,arg[0]); int ny = force->inumeric(FLERR,arg[1]); int nz = force->inumeric(FLERR,arg[2]); int nrep = nx*ny*nz; // error and warning checks if (nx <= 0 || ny <= 0 || nz <= 0) error->all(FLERR,"Illegal replicate command"); if (domain->dimension == 2 && nz != 1) error->all(FLERR,"Cannot replicate 2d simulation in z dimension"); if ((nx > 1 && domain->xperiodic == 0) || (ny > 1 && domain->yperiodic == 0) || (nz > 1 && domain->zperiodic == 0)) { if (comm->me == 0) error->warning(FLERR,"Replicating in a non-periodic dimension"); } if (atom->nextra_grow || atom->nextra_restart || atom->nextra_store) error->all(FLERR,"Cannot replicate with fixes that store atom quantities"); + // Record wall time for atom replication + + MPI_Barrier(world); + double time1 = MPI_Wtime(); + // maxtag = largest atom tag across all existing atoms tagint maxtag = 0; if (atom->tag_enable) { for (i = 0; i < atom->nlocal; i++) maxtag = MAX(atom->tag[i],maxtag); tagint maxtag_all; MPI_Allreduce(&maxtag,&maxtag_all,1,MPI_LMP_TAGINT,MPI_MAX,world); maxtag = maxtag_all; } // maxmol = largest molecule tag across all existing atoms tagint maxmol = 0; if (atom->molecule_flag) { for (i = 0; i < atom->nlocal; i++) maxmol = MAX(atom->molecule[i],maxmol); tagint maxmol_all; MPI_Allreduce(&maxmol,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world); maxmol = maxmol_all; } // unmap existing atoms via image flags for (i = 0; i < atom->nlocal; i++) domain->unmap(atom->x[i],atom->image[i]); // communication buffer for all my atom's info // max_size = largest buffer needed by any proc // must do before new Atom class created, // since size_restart() uses atom->nlocal int max_size; int send_size = atom->avec->size_restart(); MPI_Allreduce(&send_size,&max_size,1,MPI_INT,MPI_MAX,world); double *buf; memory->create(buf,max_size,"replicate:buf"); // old = original atom class // atom = new replicated atom class // also set atomKK for Kokkos version of Atom class Atom *old = atom; atomKK = NULL; if (lmp->kokkos) atom = atomKK = new AtomKokkos(lmp); else atom = new Atom(lmp); atom->settings(old); atom->create_avec(old->atom_style,old->avec->nargcopy,old->avec->argcopy,0); // check that new system will not be too large // new tags cannot exceed MAXTAGINT // new system sizes cannot exceed MAXBIGINT if (atom->tag_enable) { bigint maxnewtag = maxtag + (nrep-1)*old->natoms; if (maxnewtag < 0 || maxnewtag >= MAXTAGINT) error->all(FLERR,"Replicated system atom IDs are too big"); } if (nrep*old->natoms < 0 || nrep*old->natoms >= MAXBIGINT || nrep*old->nbonds < 0 || nrep*old->nbonds >= MAXBIGINT || nrep*old->nangles < 0 || nrep*old->nangles >= MAXBIGINT || nrep*old->ndihedrals < 0 || nrep*old->ndihedrals >= MAXBIGINT || nrep*old->nimpropers < 0 || nrep*old->nimpropers >= MAXBIGINT) error->all(FLERR,"Replicated system is too big"); // assign atom and topology counts in new class from old one atom->natoms = old->natoms * nrep; atom->nbonds = old->nbonds * nrep; atom->nangles = old->nangles * nrep; atom->ndihedrals = old->ndihedrals * nrep; atom->nimpropers = old->nimpropers * nrep; atom->ntypes = old->ntypes; atom->nbondtypes = old->nbondtypes; atom->nangletypes = old->nangletypes; atom->ndihedraltypes = old->ndihedraltypes; atom->nimpropertypes = old->nimpropertypes; atom->bond_per_atom = old->bond_per_atom; atom->angle_per_atom = old->angle_per_atom; atom->dihedral_per_atom = old->dihedral_per_atom; atom->improper_per_atom = old->improper_per_atom; // store old simulation box int triclinic = domain->triclinic; double old_xprd = domain->xprd; double old_yprd = domain->yprd; double old_zprd = domain->zprd; double old_xy = domain->xy; double old_xz = domain->xz; double old_yz = domain->yz; // setup new simulation box domain->boxhi[0] = domain->boxlo[0] + nx*old_xprd; domain->boxhi[1] = domain->boxlo[1] + ny*old_yprd; domain->boxhi[2] = domain->boxlo[2] + nz*old_zprd; if (triclinic) { domain->xy *= ny; domain->xz *= nz; domain->yz *= nz; } // new problem setup using new box boundaries if (nprocs == 1) n = static_cast (atom->natoms); else n = static_cast (LB_FACTOR * atom->natoms / nprocs); atom->allocate_type_arrays(); atom->avec->grow(n); n = atom->nmax; domain->print_box(" "); domain->set_initial_box(); domain->set_global_box(); comm->set_proc_grid(); domain->set_local_box(); // copy type arrays to new atom class if (atom->mass) { for (int itype = 1; itype <= atom->ntypes; itype++) { atom->mass_setflag[itype] = old->mass_setflag[itype]; if (atom->mass_setflag[itype]) atom->mass[itype] = old->mass[itype]; } } // set bounds for my proc // if periodic and I am lo/hi proc, adjust bounds by EPSILON // insures all replicated atoms will be owned even with round-off double epsilon[3]; if (triclinic) epsilon[0] = epsilon[1] = epsilon[2] = EPSILON; else { epsilon[0] = domain->prd[0] * EPSILON; epsilon[1] = domain->prd[1] * EPSILON; epsilon[2] = domain->prd[2] * EPSILON; } double sublo[3],subhi[3]; if (triclinic == 0) { sublo[0] = domain->sublo[0]; subhi[0] = domain->subhi[0]; sublo[1] = domain->sublo[1]; subhi[1] = domain->subhi[1]; sublo[2] = domain->sublo[2]; subhi[2] = domain->subhi[2]; } else { sublo[0] = domain->sublo_lamda[0]; subhi[0] = domain->subhi_lamda[0]; sublo[1] = domain->sublo_lamda[1]; subhi[1] = domain->subhi_lamda[1]; sublo[2] = domain->sublo_lamda[2]; subhi[2] = domain->subhi_lamda[2]; } if (comm->layout != LAYOUT_TILED) { if (domain->xperiodic) { if (comm->myloc[0] == 0) sublo[0] -= epsilon[0]; if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] += epsilon[0]; } if (domain->yperiodic) { if (comm->myloc[1] == 0) sublo[1] -= epsilon[1]; if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] += epsilon[1]; } if (domain->zperiodic) { if (comm->myloc[2] == 0) sublo[2] -= epsilon[2]; if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] += epsilon[2]; } } else { if (domain->xperiodic) { if (comm->mysplit[0][0] == 0.0) sublo[0] -= epsilon[0]; if (comm->mysplit[0][1] == 1.0) subhi[0] += epsilon[0]; } if (domain->yperiodic) { if (comm->mysplit[1][0] == 0.0) sublo[1] -= epsilon[1]; if (comm->mysplit[1][1] == 1.0) subhi[1] += epsilon[1]; } if (domain->zperiodic) { if (comm->mysplit[2][0] == 0.0) sublo[2] -= epsilon[2]; if (comm->mysplit[2][1] == 1.0) subhi[2] += epsilon[2]; } } // loop over all procs // if this iteration of loop is me: // pack my unmapped atom data into buf // bcast it to all other procs // performs 3d replicate loop with while loop over atoms in buf // x = new replicated position, remapped into simulation box // unpack atom into new atom class from buf if I own it // adjust tag, mol #, coord, topology info as needed AtomVec *old_avec = old->avec; AtomVec *avec = atom->avec; int ix,iy,iz; tagint atom_offset,mol_offset; imageint image; double x[3],lamda[3]; double *coord; int tag_enable = atom->tag_enable; for (int iproc = 0; iproc < nprocs; iproc++) { if (me == iproc) { n = 0; for (i = 0; i < old->nlocal; i++) n += old_avec->pack_restart(i,&buf[n]); } MPI_Bcast(&n,1,MPI_INT,iproc,world); MPI_Bcast(buf,n,MPI_DOUBLE,iproc,world); for (ix = 0; ix < nx; ix++) { for (iy = 0; iy < ny; iy++) { for (iz = 0; iz < nz; iz++) { // while loop over one proc's atom list m = 0; while (m < n) { image = ((imageint) IMGMAX << IMG2BITS) | ((imageint) IMGMAX << IMGBITS) | IMGMAX; if (triclinic == 0) { x[0] = buf[m+1] + ix*old_xprd; x[1] = buf[m+2] + iy*old_yprd; x[2] = buf[m+3] + iz*old_zprd; } else { x[0] = buf[m+1] + ix*old_xprd + iy*old_xy + iz*old_xz; x[1] = buf[m+2] + iy*old_yprd + iz*old_yz; x[2] = buf[m+3] + iz*old_zprd; } domain->remap(x,image); if (triclinic) { domain->x2lamda(x,lamda); coord = lamda; } else coord = x; if (coord[0] >= sublo[0] && coord[0] < subhi[0] && coord[1] >= sublo[1] && coord[1] < subhi[1] && coord[2] >= sublo[2] && coord[2] < subhi[2]) { m += avec->unpack_restart(&buf[m]); i = atom->nlocal - 1; if (tag_enable) atom_offset = iz*ny*nx*maxtag + iy*nx*maxtag + ix*maxtag; else atom_offset = 0; mol_offset = iz*ny*nx*maxmol + iy*nx*maxmol + ix*maxmol; atom->x[i][0] = x[0]; atom->x[i][1] = x[1]; atom->x[i][2] = x[2]; atom->tag[i] += atom_offset; atom->image[i] = image; if (atom->molecular) { if (atom->molecule[i] > 0) atom->molecule[i] += mol_offset; if (atom->molecular == 1) { if (atom->avec->bonds_allow) for (j = 0; j < atom->num_bond[i]; j++) atom->bond_atom[i][j] += atom_offset; if (atom->avec->angles_allow) for (j = 0; j < atom->num_angle[i]; j++) { atom->angle_atom1[i][j] += atom_offset; atom->angle_atom2[i][j] += atom_offset; atom->angle_atom3[i][j] += atom_offset; } if (atom->avec->dihedrals_allow) for (j = 0; j < atom->num_dihedral[i]; j++) { atom->dihedral_atom1[i][j] += atom_offset; atom->dihedral_atom2[i][j] += atom_offset; atom->dihedral_atom3[i][j] += atom_offset; atom->dihedral_atom4[i][j] += atom_offset; } if (atom->avec->impropers_allow) for (j = 0; j < atom->num_improper[i]; j++) { atom->improper_atom1[i][j] += atom_offset; atom->improper_atom2[i][j] += atom_offset; atom->improper_atom3[i][j] += atom_offset; atom->improper_atom4[i][j] += atom_offset; } } } } else m += static_cast (buf[m]); } } } } } // free communication buffer and old atom class memory->destroy(buf); delete old; // check that all atoms were assigned to procs bigint natoms; bigint nblocal = atom->nlocal; MPI_Allreduce(&nblocal,&natoms,1,MPI_LMP_BIGINT,MPI_SUM,world); if (me == 0) { if (screen) fprintf(screen," " BIGINT_FORMAT " atoms\n",natoms); if (logfile) fprintf(logfile," " BIGINT_FORMAT " atoms\n",natoms); } if (natoms != atom->natoms) error->all(FLERR,"Replicate did not assign all atoms correctly"); if (me == 0) { if (atom->nbonds) { if (screen) fprintf(screen," " BIGINT_FORMAT " bonds\n",atom->nbonds); if (logfile) fprintf(logfile," " BIGINT_FORMAT " bonds\n",atom->nbonds); } if (atom->nangles) { if (screen) fprintf(screen," " BIGINT_FORMAT " angles\n", atom->nangles); if (logfile) fprintf(logfile," " BIGINT_FORMAT " angles\n", atom->nangles); } if (atom->ndihedrals) { if (screen) fprintf(screen," " BIGINT_FORMAT " dihedrals\n", atom->ndihedrals); if (logfile) fprintf(logfile," " BIGINT_FORMAT " dihedrals\n", atom->ndihedrals); } if (atom->nimpropers) { if (screen) fprintf(screen," " BIGINT_FORMAT " impropers\n", atom->nimpropers); if (logfile) fprintf(logfile," " BIGINT_FORMAT " impropers\n", atom->nimpropers); } } // check that atom IDs are valid atom->tag_check(); // create global mapping of atoms if (atom->map_style) { atom->map_init(); atom->map_set(); } // create special bond lists for molecular systems if (atom->molecular == 1) { Special special(lmp); special.build(); } + + // Wall time + + MPI_Barrier(world); + double time2 = MPI_Wtime(); + + if (me == 0) { + if (screen) + fprintf(screen," Time spent = %g secs\n",time2-time1); + if (logfile) + fprintf(logfile," Time spent = %g secs\n",time2-time1); + } }