From f0ec2e3279ccb08e86e2b47950c660d3abf5cd99 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 1 Mar 2019 18:47:34 -0700 Subject: [PATCH] refactoring of global and local hyper, including amended doc pages --- doc/src/compute.txt | 7 +- doc/src/fix.txt | 6 +- doc/src/fix_ave_histo.txt | 19 +- doc/src/fix_hyper_global.txt | 45 +- doc/src/fix_hyper_local.txt | 225 ++++++---- examples/hyper/in.hyper.global | 6 +- examples/hyper/in.hyper.local | 6 + src/REPLICA/fix_hyper_global.cpp | 157 +++---- src/REPLICA/fix_hyper_global.h | 6 +- src/REPLICA/fix_hyper_local.cpp | 691 ++++++++++++++++++------------- src/REPLICA/fix_hyper_local.h | 89 ++-- src/REPLICA/hyper.cpp | 158 +++---- src/fix_ave_histo.cpp | 84 ++-- src/thermo.cpp | 3 +- 14 files changed, 868 insertions(+), 634 deletions(-) diff --git a/doc/src/compute.txt b/doc/src/compute.txt index 4886e73ab6..87dbee57d6 100644 --- a/doc/src/compute.txt +++ b/doc/src/compute.txt @@ -54,9 +54,10 @@ local quantities have the word "local" in their style, e.g. {bond/local}. Styles with neither "atom" or "local" in their style produce global quantities. -Note that a single compute produces either global or per-atom or local -quantities, but never more than one of these (with only a few -exceptions, as documented by individual compute commands). +Note that a single compute can produce either global or per-atom or +local quantities, but not both global and per-atom. It can produce +local quantities in tandem with global or per-atom quantities. The +compute doc page will explain. Global, per-atom, and local quantities each come in three kinds: a single scalar value, a vector of values, or a 2d array of values. The diff --git a/doc/src/fix.txt b/doc/src/fix.txt index 395f2ad7a9..9b48ee13bf 100644 --- a/doc/src/fix.txt +++ b/doc/src/fix.txt @@ -83,8 +83,10 @@ not in the specified fix group. Local quantities are calculated by each processor based on the atoms it owns, but there may be zero or more per atoms. -Note that a single fix may produces either global or per-atom or local -quantities (or none at all), but never more than one of these. +Note that a single fix can produce either global or per-atom or local +quantities (or none at all), but not both global and per-atom. It can +produce local quantities in tandem with global or per-atom quantities. +The fix doc page will explain. Global, per-atom, and local quantities each come in three kinds: a single scalar value, a vector of values, or a 2d array of values. The diff --git a/doc/src/fix_ave_histo.txt b/doc/src/fix_ave_histo.txt index 79f4d53481..6f99685682 100644 --- a/doc/src/fix_ave_histo.txt +++ b/doc/src/fix_ave_histo.txt @@ -35,6 +35,7 @@ keyword = {mode} or {file} or {ave} or {start} or {beyond} or {overwrite} or {ti {mode} arg = {scalar} or {vector} scalar = all input values are scalars vector = all input values are vectors + {kind} arg = {global} or {peratom} or {local} {file} arg = filename filename = name of file to output histogram(s) to {ave} args = {one} or {running} or {window} @@ -92,7 +93,8 @@ either all global, all per-atom, or all local quantities. Inputs of different kinds (e.g. global and per-atom) cannot be mixed. Atom attributes are per-atom vector values. See the doc page for individual "compute" and "fix" commands to see what kinds of -quantities they generate. +quantities they generate. See the optional {kind} keyword below for +how to force the fix ave/histo command to dis-ambiguate if necessary. Note that the output of this command is a single histogram for all input values combined together, not one histogram per input value. @@ -231,6 +233,14 @@ keyword is set to {vector}, then all input values must be global or per-atom or local vectors, or columns of global or per-atom or local arrays. +The {kind} keyword only needs to be set if a compute or fix produces +more than one kind of output (global, per-atom, local). If this is +not the case, then LAMMPS will determine what kind of input is +provided and whether all the input arguments are consistent. If a +compute or fix produces more than one kind of output, the {kind} +keyword should be used to specify which output will be used. The +remaining input arguments must still be consistent. + The {beyond} keyword determines how input values that fall outside the {lo} to {hi} bounds are treated. Values such that {lo} <= value <= {hi} are assigned to one bin. Values on a bin boundary are assigned @@ -240,7 +250,7 @@ If {beyond} is set to {end} then values < {lo} are counted in the first bin and values > {hi} are counted in the last bin. If {beyond} is set to {extend} then two extra bins are created, so that there are Nbins+2 total bins. Values < {lo} are counted in the first bin and -values > {hi} are counted in the last bin (Nbins+1). Values between +values > {hi} are counted in the last bin (Nbins+2). Values between {lo} and {hi} (inclusive) are counted in bins 2 through Nbins+1. The "coordinate" stored and printed for these two extra bins is {lo} and {hi}. @@ -354,5 +364,6 @@ ave/chunk"_fix_ave_chunk.html, "fix ave/time"_fix_ave_time.html, [Default:] none -The option defaults are mode = scalar, ave = one, start = 0, no file -output, beyond = ignore, and title 1,2,3 = strings as described above. +The option defaults are mode = scalar, kind = figured out from input +arguments, ave = one, start = 0, no file output, beyond = ignore, and +title 1,2,3 = strings as described above. diff --git a/doc/src/fix_hyper_global.txt b/doc/src/fix_hyper_global.txt index a7a938b144..81404ac6a2 100644 --- a/doc/src/fix_hyper_global.txt +++ b/doc/src/fix_hyper_global.txt @@ -102,7 +102,7 @@ Bi = exp(beta * Vij(max)) :pre where beta = 1/kTequil, and {Tequil} is the temperature of the system and an argument to this fix. Note that Bi >= 1 at every step. -NOTE: To run GHD, the input script must also use the "fix +NOTE: To run a GHD simulation, the input script must also use the "fix langevin"_fix_langevin.html command to thermostat the atoms at the same {Tequil} as specified by this fix, so that the system is running constant-temperature (NVT) dynamics. LAMMPS does not check that this @@ -166,9 +166,9 @@ correctly. There will just be fewer events because the hyper time NOTE: If you have no physical intuition as to the smallest barrier height in your system, a reasonable strategy to determine the largest -{Vmax} you can use for an LHD model, is to run a sequence of +{Vmax} you can use for a GHD model, is to run a sequence of simulations with smaller and smaller {Vmax} values, until the event -rate does not change. +rate does not change (as a function of hyper time). The {Tequil} argument is the temperature at which the system is simulated; see the comment above about the "fix @@ -177,7 +177,8 @@ beta term in the exponential factor that determines how much boost is achieved as a function of the bias potential. In general, the lower the value of {Tequil} and the higher the value -of {Vmax}, the more boost will be achievable by the GHD algorithm. +of {Vmax}, the more time boost will be achievable by the GHD +algorithm. :line @@ -190,41 +191,43 @@ The "fix_modify"_fix_modify.html {energy} option is supported by this fix to add the energy of the bias potential to the the system's potential energy as part of "thermodynamic output"_thermo_style.html. -This fix computes a global scalar and global vector of length 11, which +This fix computes a global scalar and global vector of length 12, which can be accessed by various "output commands"_Howto_output.html. The scalar is the magnitude of the bias potential (energy units) applied on the current timestep. The vector stores the following quantities: 1 = boost factor on this step (unitless) -2 = max strain Eij of any bond on this step (unitless) +2 = max strain Eij of any bond on this step (absolute value, unitless) 3 = ID of first atom in the max-strain bond 4 = ID of second atom in the max-strain bond 5 = average # of bonds/atom on this step :ul -6 = fraction of timesteps with bias = 0.0 during this run -7 = max drift distance of any atom during this run (distance units) -8 = max bond length during this run (distance units) :ul +6 = fraction of timesteps where the biased bond has bias = 0.0 during this run +7 = fraction of timesteps where the biased bond has negative strain during this run +8 = max drift distance of any atom during this run (distance units) +9 = max bond length during this run (distance units) :ul -9 = cumulative hyper time since fix was defined (time units) -10 = cumulative count of event timesteps since fix was defined -11 = cumulative count of atoms in events since fix was defined :ul +10 = cumulative hyper time since fix was defined (time units) +11 = cumulative count of event timesteps since fix was defined +12 = cumulative count of atoms in events since fix was defined :ul -The first 5 quantities are for the current timestep. Quantities 6-8 -are for the current hyper run. Quantities 9-11 are cumulative across -multiple runs (since the fix was defined in the input script). +The first 5 quantities are for the current timestep. Quantities 6-9 +are for the current hyper run. They are reset each time a new hyper +run is performed. Quantities 19-12 are cumulative across multiple +runs (since the point in the input script the fix was defined). -For value 7, drift is the distance an atom moves between timesteps -when the bond list is reset, i.e. between events. Atoms involved in -an event will typically move the greatest distance since others are -typically oscillating around their lattice site. +For value 8, drift is the distance an atom moves between two quenched +states when the second quench determines an event has occurred. Atoms +involved in an event will typically move the greatest distance since +others typically remain near their original quenched position. -For value 10, events are checked for by the "hyper"_hyper.html command +For value 11, events are checked for by the "hyper"_hyper.html command once every {Nevent} timesteps. This value is the count of those timesteps on which one (or more) events was detected. It is NOT the number of distinct events, since more than one event may occur in the same {Nevent} time window. -For value 11, each time the "hyper"_hyper.html command checks for an +For value 12, each time the "hyper"_hyper.html command checks for an event, it invokes a compute to flag zero or more atoms as participating in one or more events. E.g. atoms that have displaced more than some distance from the previous quench state. Value 11 is diff --git a/doc/src/fix_hyper_local.txt b/doc/src/fix_hyper_local.txt index c34b9ba7da..7f12e37999 100644 --- a/doc/src/fix_hyper_local.txt +++ b/doc/src/fix_hyper_local.txt @@ -22,10 +22,9 @@ Dcut = minimum distance between boosted bonds (distance units) :l alpha = boostostat relaxation time (time units) :l Btarget = desired time boost factor (unitless) :l zero or more keyword/value pairs may be appended :l -keyword = {lost} or {check/bias} or {check/coeff} - {lostbond} value = error/warn/ignore - {check/bias} values = Nevery error/warn/ignore - {check/coeff} values = Nevery error/warn/ignore :pre +keyword = {check/ghost} or {check/bias} :l + {check/ghost} values = none + {check/bias} values = Nevery error/warn/ignore :pre :ule [Examples:] @@ -65,8 +64,8 @@ To understand this description, you should first read the description of the GHD algorithm on the "fix hyper/global"_fix_hyper_global.html doc page. This description of LHD builds on the GHD description. -The definition of bonds, Eij, and Emax are the same for GHD and LHD. -The formulas for Vij(max) and Fij(max) are also the same except for a +The definition of bonds and Eij are the same for GHD and LHD. The +formulas for Vij(max) and Fij(max) are also the same except for a pre-factor Cij, explained below. The bias energy Vij applied to a bond IJ with maximum strain is @@ -117,11 +116,11 @@ where Vkl(max) is the bias energy of the maxstrain bond KL within bond IJ's neighborhood, beta = 1/kTequil, and {Tequil} is the temperature of the system and an argument to this fix. -NOTE: To run LHD, the input script must also use the "fix -langevin"_fix_langevin.html command to thermostat the atoms at the -same {Tequil} as specified by this fix, so that the system is running -constant-temperature (NVT) dynamics. LAMMPS does not check that this -is done. +NOTE: To run an LHD simulation, the input script must also use the +"fix langevin"_fix_langevin.html command to thermostat the atoms at +the same {Tequil} as specified by this fix, so that the system is +running constant-temperature (NVT) dynamics. LAMMPS does not check +that this is done. Note that if IJ = KL, then bond IJ is a biased bond on that timestep, otherwise it is not. But regardless, the boost factor Bij can be @@ -216,20 +215,20 @@ each pair. E.g. something like 2x the cutoff of the interatomic potential. In practice a {Dcut} value of ~10 Angstroms seems to work well for many solid-state systems. -NOTE: You must also insure that ghost atom communication is performed -for a distance of at least {Dcut} + {cutevent} where {cutevent} = the -distance one or more atoms move (between quenched states) to be -considered an "event". It is an argument to the "compute -event/displace" command used to detect events. By default the ghost -communication distance is set by the pair_style cutoff, which will -typically be < {Dcut}. The "comm_modify cutoff"_comm_modify.html -command can be used to set the ghost cutoff explicitly, e.g. +NOTE: You should insure that ghost atom communication is performed for +a distance of at least {Dcut} + {cutevent} = the distance one or more +atoms move (between quenched states) to be considered an "event". It +is an argument to the "compute event/displace" command used to detect +events. By default the ghost communication distance is set by the +pair_style cutoff, which will typically be < {Dcut}. The "comm_modify +cutoff"_comm_modify.html command should be used to override the ghost +cutoff explicitly, e.g. comm_modify cutoff 12.0 :pre -This fix does not know the {cutevent} parameter, but uses half the -bond length as an estimate to warn if the ghost cutoff is not long -enough. +Note that this fix does not know the {cutevent} parameter, but uses +half the {cutbond} parameter as an estimate to warn if the ghost +cutoff is not long enough. As described above the {alpha} argument is a pre-factor in the boostostat update equation for each bond's Cij prefactor. {Alpha} is @@ -269,7 +268,30 @@ NOTE: If you have no physical intuition as to the smallest barrier height in your system, a reasonable strategy to determine the largest {Btarget} you can use for an LHD model, is to run a sequence of simulations with smaller and smaller {Btarget} values, until the event -rate does not change. +rate does not change (as a function of hyper time). + +:line + +Here is additional information on the optional keywords for this fix. + +The {check/ghost} keyword turns on extra computation each timestep to +compute statistics about ghost atoms used to determine which bonds to +bias. The output of these stats are the vector values 14 and 15, +described below. If this keyword is not enabled, the output +of the stats will be zero. + +The {check/bias} keyword turns on extra computation and communcation +to check if any biased bonds are closer than {Dcut} to each other, +which should not be the case if LHD is operating correctly. Thus it +is a debugging check. The {Nevery} setting determines how often the +check is made. The {error}, {warn}, or {ignore} setting determines +what is done if the count of too-close bonds is not zero. Either the +code will exit, or issue a warning, or silently tally the count. The +count can be output as vector value 17, as described below. If this +keyword is not enabled, the output of that statistic will be 0. + +Note that both of these computations are costly, hence they are only +enabled by these keywords. :line @@ -282,95 +304,120 @@ The "fix_modify"_fix_modify.html {energy} option is supported by this fix to add the energy of the bias potential to the the system's potential energy as part of "thermodynamic output"_thermo_style.html. -This fix computes a global scalar and global vector of length 23, -which can be accessed by various "output -commands"_Howto_output.html. The scalar is the magnitude of -the bias potential (energy units) applied on the current timestep, -summed over all biased bonds. The vector stores the following -quantities: +This fix computes a global scalar and global vector of length 21, +which can be accessed by various "output commands"_Howto_output.html. +The scalar is the magnitude of the bias potential (energy units) +applied on the current timestep, summed over all biased bonds. The +vector stores the following quantities: 1 = # of biased bonds on this step -2 = max strain Eij of any bond on this step (unitless) -3 = average bias potential for all biased bonds on this step (energy units) +2 = max strain Eij of any bond on this step (absolute value, unitless) +3 = average bias coeff for all bonds on this step (unitless) 4 = average # of bonds/atom on this step 5 = average neighbor bonds/bond on this step within {Dcut} :ul -6 = fraction of steps and bonds with no bias during this run -7 = max drift distance of any atom during this run (distance units) -8 = max bond length during this run (distance units) -9 = average # of biased bonds/step during this run -10 = average bias potential for all biased bonds during this run (energy units) -11 = max bias potential for any biased bond during this run (energy units) -12 = min bias potential for any biased bond during this run (energy units) -13 = max distance from my sub-box of any ghost atom with maxstrain < qfactor during this run (distance units) -14 = max distance outside my box of any ghost atom with any maxstrain during this run (distance units) -15 = count of ghost neighbor atoms not found on reneighbor steps during this run -16 = count of lost bond partners during this run -17 = average bias coeff for lost bond partners during this run -18 = count of bias overlaps found during this run -19 = count of non-matching bias coefficients found during this run :ul +6 = max bond length during this run (distance units) +7 = average # of biased bonds/step during this run +8 = fraction of biased bonds with no bias during this run +9 = fraction of biased bonds with negative strain during this run +10 = average bias coeff for all bonds during this run (unitless) +11 = min bias coeff for any bond during this run (unitless) +12 = max bias coeff for any bond during this run (unitless) -20 = cumulative hyper time since fix created (time units) -21 = cumulative count of event timesteps since fix created -22 = cumulative count of atoms in events since fix created -23 = cumulative # of new bonds since fix created :ul +13 = max drift distance of any bond atom during this run (distance units) +14 = max distance from proc subbox of any ghost atom with maxstrain < qfactor during this run (distance units) +15 = max distance outside my box of any ghost atom with any maxstrain during this run (distance units) +16 = count of ghost atoms that could not be found on reneighbor steps during this run +17 = count of bias overlaps (< Dcut) found during this run + +18 = cumulative hyper time since fix created (time units) +19 = cumulative count of event timesteps since fix created +20 = cumulative count of atoms in events since fix created +21 = cumulative # of new bonds formed since fix created :ul The first quantities (1-5) are for the current timestep. Quantities -6-19 are for the current hyper run. They are reset each time a new -hyper run is performed. Quantities 20-23 are cumulative across -multiple runs (since the fix was defined in the input script). +6-17 are for the current hyper run. They are reset each time a new +hyper run is performed. Quantities 18-21 are cumulative across +multiple runs (since the point in the input script the fix was +defined). -For value 6, the numerator is a count of all biased bonds on every +For value 8, the numerator is a count of all biased bonds on each timestep whose bias energy = 0.0 due to Eij >= {qfactor}. The denominator is the count of all biased bonds on all timesteps. -For value 7, drift is the distance an atom moves between timesteps -when the bond list is reset, i.e. between events. Atoms involved in -an event will typically move the greatest distance since others are -typically oscillating around their lattice site. +For value 9, the numerator is a count of all biased bonds on each +timestep with negative strain. The denominator is the count of all +biased bonds on all timesteps. -For values 13 and 14, the maxstrain of a ghost atom is the maxstrain -of any bond it is part of, and it is checked for ghost atoms within -the bond neighbor cutoff. +Values 13-17 are mostly useful for debugging and diagnostic purposes. -Values 15-19 are mostly useful for debugging and diagnostic purposes. +For value 13, drift is the distance an atom moves between two quenched +states when the second quench determines an event has occurred. Atoms +involved in an event will typically move the greatest distance since +others typically remain near their original quenched position. -For values 15-17, it is possible that a ghost atom owned by another -processor will move far enough (e.g. as part of an event-in-progress) -that it will no longer be within the communication cutoff distance for -acquiring ghost atoms. Likewise it may be a ghost atom bond partner -that cannot be found because it has moved too far. These values count -those occurrences. Because they typically involve atoms that are part -of events, they do not usually indicate bad dynamics. Value 16 is the -average bias coefficient for bonds where a partner atom was lost. +For values 14-16, neighbor atoms in the full neighbor list with cutoff +{Dcut} may be ghost atoms outside a processor's sub-box. Before the +next event occurs they may move further than {Dcut} away from the +sub-box boundary. Value 14 is the furthest (from the sub-box) any +ghost atom in the neighbor list with maxstrain < {qfactor} was +accessed during the run. Value 15 is the same except that the ghost +atom's maxstrain may be >= {qfactor}, which may mean it is about to +participate in an event. Value 16 is a count of how many ghost atoms +could not be found on reneighbor steps, presumably because they moved +too far away due to their participation in an event (which will likely +be detected at the next quench). -For value 18, no two bonds should be biased if they are within a +Typical values for 14 and 15 should be slightly larger than {Dcut}, +which accounts for ghost atoms initially at a {Dcut} distance moving +thermally before the next event takes place. + +Note that for values 14 and 15 to be computed, the optional keyword +{check/ghost} must be specified. Otherwise these values will be zero. +This is because computing them incurs overhead, so the values are only +computed if requested. + +Value 16 should be zero or small. As explained above a small count +likely means some ghost atoms were participating in their own events +and moved a longer distance. If the value is large, it likely means +the communication cutoff for ghosts is too close to {Dcut} leading to +many not-found ghost atoms before the next event. This may lead to a +reduced number of bonds being selected for biasing, since the code +assumes those atoms are part of highly strained bonds. As explained +above, the "comm_modify cutoff"_comm_modify.html command can be used +to set a longer cutoff. + +For value 17, no two bonds should be biased if they are within a {Dcut} distance of each other. This value should be zero, indicating -that no pair of bonds "overlap", meaning they are closer than {Dcut} -from each other. +that no pair of biased bonds are closer than {Dcut} from each other. -For value 19, the same bias coefficient is stored by both atoms in an -IJ bond. This value should be zero, indicating that for all bonds, -each atom in the bond stores the a bias coefficient with the same -value. +Note that for values 17 to be computed, the optional keyword +{check/bias} must be specified and it determines how often this check +is performed. This is because performing the check incurs overhead, +so if only computed as often as requested. -Value 20 is simply the specified {boost} factor times the number of -timestep times the timestep size. +The result at the end of the run is the cumulative total from every +timestep the check was made. Note that the value is a count of atoms +in bonds which found other atoms in bonds too close, so it is almost +always an over-count of the number of too-close bonds. -For value 21, events are checked for by the "hyper"_hyper.html command +Value 18 is simply the specified {boost} factor times the number of +timesteps times the timestep size. + +For value 19, events are checked for by the "hyper"_hyper.html command once every {Nevent} timesteps. This value is the count of those timesteps on which one (or more) events was detected. It is NOT the number of distinct events, since more than one event may occur in the same {Nevent} time window. -For value 22, each time the "hyper"_hyper.html command checks for an +For value 20, each time the "hyper"_hyper.html command checks for an event, it invokes a compute to flag zero or more atoms as participating in one or more events. E.g. atoms that have displaced -more than some distance from the previous quench state. Value 22 is +more than some distance from the previous quench state. Value 20 is the cumulative count of the number of atoms participating in any of the events that were found. -Value 23 tallies the number of new bonds created by the bond reset +Value 21 tallies the number of new bonds created by the bond reset operation. Bonds between a specific I,J pair of atoms may persist for the entire hyperdynamics simulation if neither I or J are involved in an event. @@ -378,6 +425,16 @@ an event. The scalar and vector values calculated by this fix are all "intensive". +This fix also computes a local vector of length the number of bonds +currently in the system. The value for each bond is its Cij prefactor +(bias coefficient). These values can be can be accessed by various +"output commands"_Howto_output.html. A particularly useful one is the +"fix ave/histo"_fix_ave_histo.html command which can be used to +histogram the Cij values to see if they are distributed reasonably +close to 1.0, which indicates a good choice of {Vmax}. + +The local values calculated by this fix are unitless. + No parameter of this fix can be used with the {start/stop} keywords of the "run"_run.html command. This fix is not invoked during "energy minimization"_minimize.html. @@ -392,7 +449,9 @@ doc page for more info. "hyper"_hyper.html, "fix hyper/global"_fix_hyper_global.html -[Default:] None +[Default:] + +The check/ghost and check/bias keywords are not enabled by default. :line diff --git a/examples/hyper/in.hyper.global b/examples/hyper/in.hyper.global index 22b3b4251b..eba5c7bf89 100644 --- a/examples/hyper/in.hyper.global +++ b/examples/hyper/in.hyper.global @@ -12,6 +12,8 @@ variable cutevent index 1.1 variable steps index 100000 variable nevent index 1000 variable zoom index 1.8 +variable seed index 826626413 +variable tol index 1.0e-15 units metal atom_style atomic @@ -45,7 +47,7 @@ neighbor 0.5 bin neigh_modify every 1 delay 5 check yes fix 1 mobile nve -fix 2 mobile langevin ${Tequil} ${Tequil} 1.0 858872873 zero yes +fix 2 mobile langevin ${Tequil} ${Tequil} 1.0 ${seed} zero yes timestep 0.005 @@ -92,4 +94,4 @@ dump_modify 1 pad 6 amap 1 3 sa 1 3 blue red green # run -hyper ${steps} ${nevent} HG event min 1.0e-6 1.0e-6 100 100 dump 1 +hyper ${steps} ${nevent} HG event min ${tol} ${tol} 1000 1000 dump 1 diff --git a/examples/hyper/in.hyper.local b/examples/hyper/in.hyper.local index ef8ed4d042..cdf478ac38 100644 --- a/examples/hyper/in.hyper.local +++ b/examples/hyper/in.hyper.local @@ -107,6 +107,12 @@ dump 1 all image 10000000 local.*.jpg v_acolor type size 1024 1024 & zoom ${zoom} adiam 2.5 view 0.0 0.0 up 0 1 0 axes yes 0.9 0.01 dump_modify 1 pad 6 amap 1 3 sa 1 3 blue red green +# test of histogramming and dump output of bias coeffs + +#fix histo all ave/histo 10 100 1000 0.9 1.1 100 f_HL & +# mode vector kind local file tmp.histo +#dump 2 all local 1000 tmp.local f_HL + # run hyper ${steps} ${nevent} HL event min ${tol} ${tol} 1000 1000 dump 1 diff --git a/src/REPLICA/fix_hyper_global.cpp b/src/REPLICA/fix_hyper_global.cpp index e43f1431a9..6924fe2d93 100644 --- a/src/REPLICA/fix_hyper_global.cpp +++ b/src/REPLICA/fix_hyper_global.cpp @@ -37,7 +37,7 @@ using namespace FixConst; // possible enhancements // should there be a virial contribution from boosted bond? -// allow newton off? see Note in pre_reverse() +// allow newton off? /* ---------------------------------------------------------------------- */ @@ -52,7 +52,7 @@ FixHyperGlobal::FixHyperGlobal(LAMMPS *lmp, int narg, char **arg) : hyperflag = 1; scalar_flag = 1; vector_flag = 1; - size_vector = 11; + size_vector = 12; global_freq = 1; extscalar = 0; extvector = 0; @@ -76,6 +76,7 @@ FixHyperGlobal::FixHyperGlobal(LAMMPS *lmp, int narg, char **arg) : maxold = 0; xold = NULL; tagold = NULL; + old2now = NULL; me = comm->me; firstflag = 1; @@ -94,6 +95,7 @@ FixHyperGlobal::~FixHyperGlobal() memory->sfree(blist); memory->destroy(xold); memory->destroy(tagold); + memory->destroy(old2now); } /* ---------------------------------------------------------------------- */ @@ -114,6 +116,7 @@ void FixHyperGlobal::init_hyper() maxdriftsq = 0.0; maxbondlen = 0.0; nobias = 0; + negstrain = 0; } /* ---------------------------------------------------------------------- */ @@ -155,14 +158,16 @@ void FixHyperGlobal::setup_pre_neighbor() void FixHyperGlobal::setup_pre_reverse(int eflag, int vflag) { - // no increment in nobias or hyper time when pre-run forces are calculated + // no increment in these quantities when pre-run forces are calculated int nobias_hold = nobias; + int negstrain_hold = negstrain; double t_hyper_hold = t_hyper; pre_reverse(eflag,vflag); nobias = nobias_hold; + negstrain = negstrain_hold; t_hyper = t_hyper_hold; } @@ -171,7 +176,7 @@ void FixHyperGlobal::setup_pre_reverse(int eflag, int vflag) void FixHyperGlobal::pre_neighbor() { int i,m,iold,jold,ilocal,jlocal; - double distsq; + // double distsq; // reset local indices for owned bond atoms, since atoms have migrated // must be done after ghost atoms are setup via comm->borders() @@ -182,6 +187,7 @@ void FixHyperGlobal::pre_neighbor() // closest_image() returns the ghost atom index in that case // also compute max drift of any atom in a bond // drift = displacement from quenched coord while event has not yet occured + // NOTE: drift calc is now done in bond_build(), between 2 quenched states for (i = 0; i < nall_old; i++) old2now[i] = -1; @@ -199,8 +205,8 @@ void FixHyperGlobal::pre_neighbor() if (ilocal < 0) error->one(FLERR,"Fix hyper/global bond atom not found"); old2now[iold] = ilocal; - distsq = MathExtra::distsq3(x[ilocal],xold[iold]); - maxdriftsq = MAX(distsq,maxdriftsq); + //distsq = MathExtra::distsq3(x[ilocal],xold[iold]); + //maxdriftsq = MAX(distsq,maxdriftsq); } if (jlocal < 0) { jlocal = atom->map(tagold[jold]); @@ -208,40 +214,13 @@ void FixHyperGlobal::pre_neighbor() if (jlocal < 0) error->one(FLERR,"Fix hyper/global bond atom not found"); old2now[jold] = jlocal; - distsq = MathExtra::distsq3(x[jlocal],xold[jold]); - maxdriftsq = MAX(distsq,maxdriftsq); + //distsq = MathExtra::distsq3(x[jlocal],xold[jold]); + //maxdriftsq = MAX(distsq,maxdriftsq); } blist[m].i = ilocal; blist[m].j = jlocal; } - - /* old way - nblocal loop is re-doing index-find calculation - - // NOTE: drift may not include J atoms moving (if not themselves bond owners) - - int flag = 0; - - for (m = 0; m < nblocal; m++) { - iold = blist[m].iold; - jold = blist[m].jold; - ilocal = atom->map(tagold[iold]); - jlocal = atom->map(tagold[jold]); - ilocal = domain->closest_image(xold[iold],ilocal); - jlocal = domain->closest_image(xold[iold],jlocal); - blist[m].i = ilocal; - blist[m].j = jlocal; - - if (ilocal < 0 || jlocal < 0) flag++; - else { - distsq = MathExtra::distsq3(x[ilocal],xold[iold]); - maxdriftsq = MAX(distsq,maxdriftsq); - } - } - - if (flag) error->one(FLERR,"Fix hyper/global bond atom not found"); - - */ } /* ---------------------------------------------------------------------- */ @@ -254,12 +233,12 @@ void FixHyperGlobal::pre_reverse(int /* eflag */, int /* vflag */) double ebias,vbias,fbias,fbiasr; // compute current strain of each owned bond - // eabs_max = maximum absolute value of strain of any bond I own + // emax = maximum abs value of strain of any bond I own // imax,jmax = local indices of my 2 atoms in that bond // rmax,r0max = current and relaxed lengths of that bond double **x = atom->x; - double estrain_maxabs = 0.0; + double emax = 0.0; for (m = 0; m < nblocal; m++) { i = blist[m].i; @@ -272,8 +251,8 @@ void FixHyperGlobal::pre_reverse(int /* eflag */, int /* vflag */) r0 = blist[m].r0; estrain = fabs(r-r0) / r0; - if (estrain > estrain_maxabs) { - estrain_maxabs = estrain; + if (estrain > emax) { + emax = estrain; rmax = r; r0max = r0; imax = i; @@ -285,7 +264,7 @@ void FixHyperGlobal::pre_reverse(int /* eflag */, int /* vflag */) // finds max strain and what proc owns it // owner = proc that owns that bond - pairme.value = estrain_maxabs; + pairme.value = emax; pairme.proc = me; MPI_Allreduce(&pairme,&pairall,1,MPI_DOUBLE_INT,MPI_MAXLOC,world); owner = pairall.proc; @@ -311,16 +290,14 @@ void FixHyperGlobal::pre_reverse(int /* eflag */, int /* vflag */) // Fix = x component of force on atom I // = Fbias dEbias/dr dr/dxi, dEbias/dr = 1/r0, dr/dxi = delx/r // dt_boost = time boost factor = exp(Vbias/kT) - // NOTE: logic here would need to be different for newton off double **f = atom->f; vbias = fbias = 0.0; dt_boost = 1.0; - if (estrain_maxabs < qfactor) { - //ebias = (rmax-r0max) / r0max; - ebias = fabs(rmax-r0max) / r0max; + if (emax < qfactor) { + ebias = (rmax-r0max) / r0max; vbias = vmax * (1.0 - ebias*ebias*invqfactorsq); fbias = 2.0 * vmax * ebias * invqfactorsq; dt_boost = exp(beta*vbias); @@ -338,13 +315,15 @@ void FixHyperGlobal::pre_reverse(int /* eflag */, int /* vflag */) f[jmax][1] -= dely*fbiasr; f[jmax][2] -= delz*fbiasr; + if (ebias < 0.0) negstrain++; + } else nobias++; // output quantities outvec[0] = vbias; outvec[1] = dt_boost; - outvec[2] = ebias; + outvec[2] = emax; outvec[3] = atom->tag[imax]; outvec[4] = atom->tag[jmax]; @@ -356,8 +335,8 @@ void FixHyperGlobal::pre_reverse(int /* eflag */, int /* vflag */) void FixHyperGlobal::build_bond_list(int natom) { - int i,j,ii,jj,inum,jnum; - double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int i,j,m,ii,jj,iold,jold,ilocal,jlocal,inum,jnum; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq,distsq; int *ilist,*jlist,*numneigh,**firstneigh; if (natom) { @@ -365,6 +344,27 @@ void FixHyperGlobal::build_bond_list(int natom) nevent_atom += natom; } + // compute max distance any bond atom has moved between 2 quenched states + // xold[iold] = last quenched coord for iold + // x[ilocal] = current quenched coord for same atom + + double **x = atom->x; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + + for (m = 0; m < nblocal; m++) { + iold = blist[m].iold; + ilocal = atom->map(tagold[iold]); + ilocal = domain->closest_image(xold[iold],ilocal); + distsq = MathExtra::distsq3(x[ilocal],xold[iold]); + maxdriftsq = MAX(distsq,maxdriftsq); + jold = blist[m].jold; + jlocal = atom->map(tagold[jold]); + jlocal = domain->closest_image(xold[iold],jlocal); + distsq = MathExtra::distsq3(x[jlocal],xold[jold]); + maxdriftsq = MAX(distsq,maxdriftsq); + } + // trigger neighbor list build neighbor->build_one(list); @@ -372,7 +372,6 @@ void FixHyperGlobal::build_bond_list(int natom) // identify bonds assigned to each owned atom // do not create a bond between two non-group atoms - double **x = atom->x; int *mask = atom->mask; inum = list->inum; @@ -415,10 +414,12 @@ void FixHyperGlobal::build_bond_list(int natom) } } - // store IDs and coords for owned+ghost atoms at time of bond creation - // realloc xold and tagold as needed + // store per-atom quantities for owned+ghost atoms at time of bond creation + // nall_old = value of nall at time bonds are built - if (atom->nmax > maxold) { + tagint *tag = atom->tag; + + if (nall > maxold) { memory->destroy(xold); memory->destroy(tagold); memory->destroy(old2now); @@ -428,16 +429,11 @@ void FixHyperGlobal::build_bond_list(int natom) memory->create(old2now,maxold,"hyper/global:old2now"); } - tagint *tag = atom->tag; - int nall = atom->nlocal + atom->nghost; - nall_old = nall; + memcpy(&xold[0][0],&x[0][0],3*nall*sizeof(double)); + for (i = 0; i < nall; i++) tagold[i] = tag[i]; - for (i = 0; i < nall; i++) { - xold[i][0] = x[i][0]; - xold[i][1] = x[i][1]; - xold[i][2] = x[i][2]; - tagold[i] = tag[i]; - } + nlocal_old = nlocal; + nall_old = nall; } /* ---------------------------------------------------------------------- @@ -473,7 +469,7 @@ double FixHyperGlobal::compute_vector(int i) bcastflag = 0; } - // 11 vector outputs returned for i = 0-10 + // 12 vector outputs returned for i = 0-11 // i = 0 = boost factor on this step // i = 1 = max strain of any bond on this step (positive or negative) @@ -481,13 +477,14 @@ double FixHyperGlobal::compute_vector(int i) // i = 3 = ID of atom J in max-strain bond on this step // i = 4 = ave bonds/atom on this step - // i = 5 = fraction of steps with no bias during this run - // i = 6 = max drift of any atom during this run - // i = 7 = max bond length during this run + // i = 5 = fraction of steps where bond has no bias during this run + // i = 6 = fraction of steps where bond has negative strain during this run + // i = 7 = max drift distance of any atom during this run + // i = 8 = max bond length during this run - // i = 8 = cummulative hyper time since fix created - // i = 9 = cummulative # of event timesteps since fix created - // i = 10 = cummulative # of atoms in events since fix created + // i = 9 = cummulative hyper time since fix created + // i = 10 = cummulative # of event timesteps since fix created + // i = 11 = cummulative # of atoms in events since fix created if (i == 0) return outvec[1]; if (i == 1) return outvec[2]; @@ -509,20 +506,27 @@ double FixHyperGlobal::compute_vector(int i) } if (i == 6) { + if (update->ntimestep == update->firststep) return 0.0; + int allnegstrain; + MPI_Allreduce(&negstrain,&allnegstrain,1,MPI_INT,MPI_SUM,world); + return 1.0*allnegstrain / (update->ntimestep - update->firststep); + } + + if (i == 7) { double alldriftsq; MPI_Allreduce(&maxdriftsq,&alldriftsq,1,MPI_DOUBLE,MPI_MAX,world); return sqrt(alldriftsq); } - if (i == 7) { + if (i == 8) { double allbondlen; MPI_Allreduce(&maxbondlen,&allbondlen,1,MPI_DOUBLE,MPI_MAX,world); return allbondlen; } - if (i == 8) return t_hyper; - if (i == 9) return (double) nevent; - if (i == 10) return (double) nevent_atom; + if (i == 9) return t_hyper; + if (i == 10) return (double) nevent; + if (i == 11) return (double) nevent_atom; return 0.0; } @@ -534,13 +538,14 @@ double FixHyperGlobal::compute_vector(int i) double FixHyperGlobal::query(int i) { - if (i == 1) return compute_vector(8); // cummulative hyper time - if (i == 2) return compute_vector(9); // nevent - if (i == 3) return compute_vector(10); // nevent_atom + if (i == 1) return compute_vector(9); // cummulative hyper time + if (i == 2) return compute_vector(10); // nevent + if (i == 3) return compute_vector(11); // nevent_atom if (i == 4) return compute_vector(4); // ave bonds/atom - if (i == 5) return compute_vector(6); // maxdrift - if (i == 6) return compute_vector(7); // maxbondlen + if (i == 5) return compute_vector(7); // maxdrift + if (i == 6) return compute_vector(8); // maxbondlen if (i == 7) return compute_vector(5); // fraction with zero bias + if (i == 8) return compute_vector(6); // fraction with negative strain error->all(FLERR,"Invalid query to fix hyper/global"); diff --git a/src/REPLICA/fix_hyper_global.h b/src/REPLICA/fix_hyper_global.h index 42dd64e145..a62f80b54a 100644 --- a/src/REPLICA/fix_hyper_global.h +++ b/src/REPLICA/fix_hyper_global.h @@ -56,6 +56,7 @@ class FixHyperGlobal : public FixHyper { double maxbondlen; // max length of any bond double maxdriftsq; // max distance any atom drifts from original pos int nobias; // # of steps when bias = 0, b/c bond too long + int negstrain; // # of steps when biased bond has negative strain class NeighList *list; @@ -70,12 +71,13 @@ class FixHyperGlobal : public FixHyper { double r0; // relaxed bond length }; - struct OneBond *blist; // list of owned bonds - int nblocal; // # of owned bonds + OneBond *blist; // list of owned bonds + int nblocal; // # of owned bonds // coords and IDs of owned+ghost atoms when bonds were formed // persists on a proc from one event until the next + int nlocal_old; // nlocal for old atoms int nall_old; // nlocal+nghost for old atoms int maxold; // allocated size of old atoms diff --git a/src/REPLICA/fix_hyper_local.cpp b/src/REPLICA/fix_hyper_local.cpp index 99dd1945ad..a2af3dff6e 100644 --- a/src/REPLICA/fix_hyper_local.cpp +++ b/src/REPLICA/fix_hyper_local.cpp @@ -36,8 +36,7 @@ using namespace FixConst; #define DELTABOND 16384 #define DELTABIAS 16 #define COEFFINIT 1.0 -#define COEFFMAX 1.2 -#define MAXBONDPERATOM 30 +#define FCCBONDS 12 #define BIG 1.0e20 enum{STRAIN,STRAINDOMAIN,BIASFLAG,BIASCOEFF}; @@ -46,12 +45,11 @@ enum{IGNORE,WARN,ERROR}; /* ---------------------------------------------------------------------- */ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) : - FixHyper(lmp, narg, arg), old2now(NULL), xold(NULL), tagold(NULL), - blist(NULL), maxstrain(NULL), maxstrain_domain(NULL), - biasflag(NULL), bias(NULL) + FixHyper(lmp, narg, arg), blist(NULL), biascoeff(NULL), numbond(NULL), + maxhalf(NULL), eligible(NULL), maxhalfstrain(NULL), old2now(NULL), + tagold(NULL), xold(NULL), maxstrain(NULL), maxstrain_domain(NULL), + biasflag(NULL), bias(NULL), cpage(NULL), clist(NULL), numcoeff(NULL) { - // NOTE: need to add vecs/arrays to constructor list - // error checks if (atom->map_style == 0) @@ -64,7 +62,11 @@ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) : hyperflag = 2; scalar_flag = 1; vector_flag = 1; - size_vector = 23; + size_vector = 21; + local_flag = 1; + size_local_rows = 0; + size_local_cols = 0; + local_freq = 1; global_freq = 1; extscalar = 0; @@ -89,11 +91,15 @@ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) : // optional args + checkghost = 0; checkbias = 0; int iarg = 10; while (iarg < narg) { - if (strcmp(arg[iarg],"check/bias") == 0) { + if (strcmp(arg[iarg],"check/ghost") == 0) { + checkghost = 1; + iarg++; + } else if (strcmp(arg[iarg],"check/bias") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal fix hyper/local command"); checkbias = 1; checkbias_every = force->inumeric(FLERR,arg[iarg+1]); @@ -102,14 +108,15 @@ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) : else if (strcmp(arg[iarg+2],"ignore") == 0) checkbias_flag = IGNORE; else error->all(FLERR,"Illegal fix hyper/local command"); iarg += 3; - } else error->all(FLERR,"Illegal fix hyper/local command"); } // per-atom data structs - maxbond = 0; + maxbond = nblocal = 0; blist = NULL; + biascoeff = NULL; + allbonds = 0; maxatom = 0; maxstrain = NULL; @@ -130,24 +137,30 @@ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) : nbias = maxbias = 0; bias = NULL; + // data structs for persisting bias coeffs when bond list is reformed + // maxbondperatom = max # of bonds any atom is part of + // FCCBONDS = 12 is a good estimate for fcc lattices + // will be reset in build_bond() if necessary + maxcoeff = 0; - maxcoeffperatom = 0; + maxbondperatom = FCCBONDS; numcoeff = NULL; clist = NULL; + cpage = new MyPage; + cpage->init(maxbondperatom,1024*maxbondperatom,1); - // maxbondperatom = max # of bonds any atom is part of - // will be reset in bond_build() // set comm sizes needed by this fix - // NOTE: remove MBPA when minimize reverse Cij comm + // reverse = 2 is for sending atom index + value, though total likely < 1 + // reverse comm for bias coeffs has variable size, so not tallied here - maxbondperatom = 1; comm_forward = 1; - comm_reverse = MAXBONDPERATOM; + comm_reverse = 2; me = comm->me; firstflag = 1; - allbias = 0.0; + sumbiascoeff = 0.0; + avebiascoeff = 0.0; starttime = update->ntimestep; nostrainyet = 1; @@ -162,6 +175,7 @@ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) : FixHyperLocal::~FixHyperLocal() { memory->destroy(blist); + memory->destroy(biascoeff); memory->destroy(maxstrain); memory->destroy(maxstrain_domain); @@ -179,7 +193,8 @@ FixHyperLocal::~FixHyperLocal() memory->destroy(bias); memory->destroy(numcoeff); - memory->destroy(clist); + memory->sfree(clist); + delete cpage; } /* ---------------------------------------------------------------------- */ @@ -202,11 +217,12 @@ void FixHyperLocal::init_hyper() checkbias_count = 0; maxdriftsq = 0.0; maxbondlen = 0.0; - maxbiascoeff = 0.0; + avebiascoeff = 0.0; minbiascoeff = BIG; - sumbiascoeff = 0.0; + maxbiascoeff = 0.0; nbias_running = 0; nobias_running = 0; + negstrain_running = 0; rmaxever = 0.0; rmaxeverbig = 0.0; @@ -254,11 +270,12 @@ void FixHyperLocal::init() // need occasional full neighbor list with cutoff = Dcut // used for finding maxstrain of neighbor bonds out to Dcut // do not need to include neigh skin in cutoff, - // b/c this list will be built every time bond_build() is called + // b/c this list will be built every time build_bond() is called // NOTE: what if pair style list cutoff > Dcut // or what if neigh skin is huge? int irequest_full = neighbor->request(this,instance_me); + neighbor->requests[irequest_full]->id = 1; neighbor->requests[irequest_full]->pair = 0; neighbor->requests[irequest_full]->fix = 1; neighbor->requests[irequest_full]->half = 0; @@ -271,9 +288,10 @@ void FixHyperLocal::init() // used for building local bond list // no specified cutoff, should be longer than cutbond // this list will also be built (or derived/copied) - // every time bond_build() is called + // every time build_bond() is called int irequest_half = neighbor->request(this,instance_me); + neighbor->requests[irequest_half]->id = 2; neighbor->requests[irequest_half]->pair = 0; neighbor->requests[irequest_half]->fix = 1; neighbor->requests[irequest_half]->occasional = 1; @@ -297,7 +315,6 @@ void FixHyperLocal::init_list(int id, NeighList *ptr) void FixHyperLocal::setup_pre_neighbor() { // called for dynamics and minimization - // NOTE: check if needed for min, I think so b/c of Cij persist? pre_neighbor(); } @@ -308,7 +325,8 @@ void FixHyperLocal::setup_pre_reverse(int eflag, int vflag) { // only called for dynamics, not minimization // setupflag prevents boostostat update of bias coeffs in setup - // also prevents increments of nbias_running, nobias_running, sumbiascoeff + // also prevents increments of nbias_running, nobias_running, + // negstrain_running, sumbiascoeff setupflag = 1; pre_reverse(eflag,vflag); @@ -320,7 +338,7 @@ void FixHyperLocal::setup_pre_reverse(int eflag, int vflag) void FixHyperLocal::pre_neighbor() { int i,m,iold,jold,ilocal,jlocal; - double distsq; + // double distsq; // reset local indices for owned bond atoms, since atoms have migrated // must be done after ghost atoms are setup via comm->borders() @@ -331,6 +349,7 @@ void FixHyperLocal::pre_neighbor() // closest_image() returns the ghost atom index in that case // also compute max drift of any atom in a bond // drift = displacement from quenched coord while event has not yet occured + // NOTE: drift calc is now done in bond_build(), between 2 quenched states for (i = 0; i < nall_old; i++) old2now[i] = -1; @@ -348,17 +367,17 @@ void FixHyperLocal::pre_neighbor() if (ilocal < 0) error->one(FLERR,"Fix hyper/local bond atom not found"); old2now[iold] = ilocal; - distsq = MathExtra::distsq3(x[ilocal],xold[iold]); - maxdriftsq = MAX(distsq,maxdriftsq); + //distsq = MathExtra::distsq3(x[ilocal],xold[iold]); + //maxdriftsq = MAX(distsq,maxdriftsq); } if (jlocal < 0) { jlocal = atom->map(tagold[jold]); - jlocal = domain->closest_image(xold[iold],jlocal); // closest to iold + jlocal = domain->closest_image(xold[iold],jlocal); // close to I atom if (jlocal < 0) error->one(FLERR,"Fix hyper/local bond atom not found"); old2now[jold] = jlocal; - distsq = MathExtra::distsq3(x[jlocal],xold[jold]); - maxdriftsq = MAX(distsq,maxdriftsq); + //distsq = MathExtra::distsq3(x[jlocal],xold[jold]); + //maxdriftsq = MAX(distsq,maxdriftsq); } blist[m].i = ilocal; @@ -374,11 +393,11 @@ void FixHyperLocal::pre_neighbor() // b/c old2now is only used to access maxstrain() or biasflag() // which will be identical for every copy of the same atom ID - for (i = 0; i < nall_old; i++) { - if (old2now[i] >= 0) continue; - if (tagold[i] == 0) continue; - ilocal = atom->map(tagold[i]); - old2now[i] = ilocal; + for (iold = 0; iold < nall_old; iold++) { + if (old2now[iold] >= 0) continue; + if (tagold[iold] == 0) continue; + ilocal = atom->map(tagold[iold]); + old2now[iold] = ilocal; if (ilocal < 0) ghost_toofar++; } } @@ -389,26 +408,31 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) { int i,j,m,ii,jj,inum,jnum,iold,jold,ibond,nbond,ijhalf,ncount; double xtmp,ytmp,ztmp,delx,dely,delz; - double r,r0,estrain,emax,ebias,vbias,fbias,fbiasr,biascoeff; + double r,r0,estrain,emax,ebias,vbias,fbias,fbiasr; double halfstrain,selfstrain; int *ilist,*jlist,*numneigh,**firstneigh; //double time1,time2,time3,time4,time5,time6,time7,time8; //time1 = MPI_Wtime(); - // reallocate local vectors if necessary + nostrainyet = 0; + + // reallocate per-atom maxstrain and biasflag vectors if necessary int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; if (maxatom < nall) { + memory->destroy(maxstrain); + memory->destroy(maxstrain_domain); + if (checkbias) memory->destroy(biasflag); maxatom = atom->nmax; - memory->grow(maxstrain,maxatom,"hyper/local:maxstrain"); - memory->grow(maxstrain_domain,maxatom,"hyper/local:maxstrain_domain"); - if (checkbias) memory->grow(biasflag,maxatom,"hyper/local:biasflag"); + memory->create(maxstrain,maxatom,"hyper/local:maxstrain"); + memory->create(maxstrain_domain,maxatom,"hyper/local:maxstrain_domain"); + if (checkbias) memory->create(biasflag,maxatom,"hyper/local:biasflag"); } - // each old atom I's owned bond with max strain is eligible for biasing + // one max strain bond per old owned atom is eligible for biasing for (iold = 0; iold < nlocal_old; iold++) eligible[iold] = 1; @@ -423,7 +447,7 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) // mark atom I ineligible if it has no bonds // also store: // maxhalf = which owned bond is maxstrain for each old atom I - // maxhalfstrain = strain of that bond for each old atom I + // maxhalfstrain = abs value strain of that bond for each old atom I for (i = 0; i < nall; i++) maxstrain[i] = 0.0; @@ -479,7 +503,6 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) // ------------------------------------------------------------- // use full Dcut neighbor list to check maxstrain of all neighbor atoms - // NOTE: is II loop the same as iold over nlocal_old ?? // neighlist is from last event // has old indices for I,J (reneighboring may have occurred) // use old2now[] to convert to current indices @@ -503,8 +526,12 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) double *sublo = domain->sublo; double *subhi = domain->subhi; + // first two lines of outer loop should be identical to this: + // for (iold = 0; iold < nlocal_old; iold++) + for (ii = 0; ii < inum; ii++) { iold = ilist[ii]; + if (eligible[iold] == 0) continue; jlist = firstneigh[iold]; jnum = numneigh[iold]; @@ -533,29 +560,30 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) emax = MAX(emax,maxstrain[j]); if (selfstrain == maxstrain[j]) ncount++; - // diagnostic + // optional diagnostic // tally largest distance from subbox that a ghost atom is (rmaxbig) // and the largest distance if strain < qfactor (rmax) - // NOTE: could this be removed from loop ?? - if (j >= nlocal) { - if (x[j][0] < sublo[0]) rmaxbig = MAX(rmaxbig,sublo[0]-x[j][0]); - if (x[j][1] < sublo[1]) rmaxbig = MAX(rmaxbig,sublo[1]-x[j][1]); - if (x[j][2] < sublo[2]) rmaxbig = MAX(rmaxbig,sublo[2]-x[j][2]); - if (x[j][0] > subhi[0]) rmaxbig = MAX(rmaxbig,x[j][0]-subhi[0]); - if (x[j][1] > subhi[1]) rmaxbig = MAX(rmaxbig,x[j][1]-subhi[1]); - if (x[j][2] > subhi[2]) rmaxbig = MAX(rmaxbig,x[j][2]-subhi[2]); - if (maxstrain[j] < qfactor) { - if (x[j][0] < sublo[0]) rmax = MAX(rmax,sublo[0]-x[j][0]); - if (x[j][1] < sublo[1]) rmax = MAX(rmax,sublo[1]-x[j][1]); - if (x[j][2] < sublo[2]) rmax = MAX(rmax,sublo[2]-x[j][2]); - if (x[j][0] > subhi[0]) rmax = MAX(rmax,x[j][0]-subhi[0]); - if (x[j][1] > subhi[1]) rmax = MAX(rmax,x[j][1]-subhi[1]); - if (x[j][2] > subhi[2]) rmax = MAX(rmax,x[j][2]-subhi[2]); + if (checkghost) { + if (j >= nlocal) { + if (x[j][0] < sublo[0]) rmaxbig = MAX(rmaxbig,sublo[0]-x[j][0]); + if (x[j][1] < sublo[1]) rmaxbig = MAX(rmaxbig,sublo[1]-x[j][1]); + if (x[j][2] < sublo[2]) rmaxbig = MAX(rmaxbig,sublo[2]-x[j][2]); + if (x[j][0] > subhi[0]) rmaxbig = MAX(rmaxbig,x[j][0]-subhi[0]); + if (x[j][1] > subhi[1]) rmaxbig = MAX(rmaxbig,x[j][1]-subhi[1]); + if (x[j][2] > subhi[2]) rmaxbig = MAX(rmaxbig,x[j][2]-subhi[2]); + if (maxstrain[j] < qfactor) { + if (x[j][0] < sublo[0]) rmax = MAX(rmax,sublo[0]-x[j][0]); + if (x[j][1] < sublo[1]) rmax = MAX(rmax,sublo[1]-x[j][1]); + if (x[j][2] < sublo[2]) rmax = MAX(rmax,sublo[2]-x[j][2]); + if (x[j][0] > subhi[0]) rmax = MAX(rmax,x[j][0]-subhi[0]); + if (x[j][1] > subhi[1]) rmax = MAX(rmax,x[j][1]-subhi[1]); + if (x[j][2] > subhi[2]) rmax = MAX(rmax,x[j][2]-subhi[2]); + } } } } - + if (maxhalfstrain[iold] < selfstrain) eligible[iold] = 0; if (selfstrain < emax) eligible[iold] = 0; else if (ncount > 1) { @@ -565,15 +593,6 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) maxstrain_domain[i] = emax; } - // diagnostic // NOTE: optional, should skip - - double rmax2[2],rmax2all[2]; - rmax2[0] = rmax; - rmax2[1] = rmaxbig; - MPI_Allreduce(&rmax2,&rmax2all,2,MPI_DOUBLE,MPI_MAX,world); - rmaxever = rmax2all[0]; - rmaxeverbig = rmax2all[1]; - //time4 = MPI_Wtime(); // reverse comm to acquire maxstrain_domain from ghost atoms @@ -607,7 +626,7 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) maxbias += DELTABIAS; memory->grow(bias,maxbias,"hyper/local:bias"); } - bias[nbias++] = ibond; + bias[nbias++] = maxhalf[iold]; } //time6 = MPI_Wtime(); @@ -620,6 +639,7 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) double **f = atom->f; int nobias = 0; + int negstrain = 0; mybias = 0.0; for (int ibias = 0; ibias < nbias; ibias++) { @@ -637,11 +657,9 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) delz = x[i][2] - x[j][2]; r = sqrt(delx*delx + dely*dely + delz*delz); r0 = blist[m].r0; - //ebias = (r-r0) / r0; - ebias = fabs(r-r0) / r0; - biascoeff = blist[m].biascoeff; - vbias = biascoeff * vmax * (1.0 - ebias*ebias*invqfactorsq); - fbias = biascoeff * 2.0 * vmax * ebias * invqfactorsq; + ebias = (r-r0) / r0; + vbias = biascoeff[m] * vmax * (1.0 - ebias*ebias*invqfactorsq); + fbias = biascoeff[m] * 2.0 * vmax * ebias * invqfactorsq; fbiasr = fbias / r0 / r; f[i][0] += delx*fbiasr; @@ -652,6 +670,7 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) f[j][1] -= dely*fbiasr; f[j][2] -= delz*fbiasr; + if (ebias < 0.0) negstrain++; mybias += vbias; } @@ -662,50 +681,59 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) // ------------------------------------------------------------- // no boostostat update when pre_reverse called from setup() - // nbias_running, nobias_running, sumbiascoeff only incremented on run steps - // NOTE: maybe should also not bias any bonds on firststep of this fix + // nbias_running, nobias_running, negstrain_running only incremented + // on run steps if (setupflag) return; nbias_running += nbias; nobias_running += nobias; + negstrain_running += negstrain; // loop over bonds I own to adjust bias coeff // delta in boost coeff is function of maxboost_domain vs target boost // maxboost_domain is function of two maxstrain_domains for I,J - double mybias = 0.0; - double emaxi,emaxj,maxboost_domain; + double emaxi,emaxj,maxboost_domain,bc; + double mybiascoeff = 0.0; for (m = 0; m < nblocal; m++) { i = blist[m].i; j = blist[m].j; - emaxi = fabs(maxstrain_domain[i]); - emaxj = fabs(maxstrain_domain[j]); + emaxi = maxstrain_domain[i]; + emaxj = maxstrain_domain[j]; emax = MAX(emaxi,emaxj); if (emax < qfactor) vbias = vmax * (1.0 - emax*emax*invqfactorsq); else vbias = 0.0; - biascoeff = blist[m].biascoeff; - maxboost_domain = exp(beta * biascoeff*vbias); - biascoeff -= alpha * (maxboost_domain-boost_target) / boost_target; - blist[m].biascoeff = biascoeff; + maxboost_domain = exp(beta * biascoeff[m]*vbias); + biascoeff[m] -= alpha * (maxboost_domain-boost_target) / boost_target; // stats - mybias += biascoeff; - maxbiascoeff = MAX(maxbiascoeff,biascoeff); - minbiascoeff = MIN(minbiascoeff,biascoeff); + bc = biascoeff[m]; + mybiascoeff += bc; + minbiascoeff = MIN(minbiascoeff,bc); + maxbiascoeff = MAX(maxbiascoeff,bc); } - // running stats - - MPI_Allreduce(&mybias,&allbias,1,MPI_DOUBLE,MPI_SUM,world); - if (allbonds) sumbiascoeff += allbias/allbonds; - // ------------------------------------------------------------- - // extra diagnostics if requested + // diagnostics, some optional // ------------------------------------------------------------- + MPI_Allreduce(&mybiascoeff,&sumbiascoeff,1,MPI_DOUBLE,MPI_SUM,world); + if (allbonds) avebiascoeff += sumbiascoeff/allbonds; + + // if requested, monitor ghost distance from processor sub-boxes + + if (checkghost) { + double rmax2[2],rmax2all[2]; + rmax2[0] = rmax; + rmax2[1] = rmaxbig; + MPI_Allreduce(&rmax2,&rmax2all,2,MPI_DOUBLE,MPI_MAX,world); + rmaxever = rmax2all[0]; + rmaxeverbig = rmax2all[1]; + } + // if requsted, check for any biased bonds that are too close to each other // keep a running count for output // requires 2 additional local comm operations @@ -713,7 +741,7 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) if (checkbias && update->ntimestep % checkbias_every == 0) { // mark each atom in a biased bond with ID of partner - // nbias loop will mark some ghost atoms + // this may mark some ghost atoms for (i = 0; i < nall; i++) biasflag[i] = 0; @@ -728,13 +756,13 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) } // reverse comm to acquire biasflag from ghost atoms - // needed b/c above loop may set biasflag of ghost atoms - // forward comm to acquire biasflag of all ghost atoms + // forward comm to set biasflag for all ghost atoms commflag = BIASFLAG; comm->reverse_comm_fix(this); comm->forward_comm_fix(this); + // loop over Dcut full neighbor list // I and J may be ghost atoms // only continue if I is a biased atom // if J is unknown (drifted ghost) just ignore @@ -755,6 +783,19 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) if (biasflag[j] && biasflag[j] != tag[i]) checkbias_count++; } } + + if (checkbias_flag != IGNORE) { + int allcount; + MPI_Allreduce(&checkbias_count,&allcount,1,MPI_INT,MPI_SUM,world); + if (allcount) { + char str[128]; + sprintf(str,"Fix hyper/local biased bonds too close: " + "cumulative atom count %d",allcount); + if (checkbias_flag == WARN) { + if (me == 0) error->warning(FLERR,str); + } else error->all(FLERR,str); + } + } } } @@ -769,9 +810,9 @@ void FixHyperLocal::min_pre_neighbor() void FixHyperLocal::build_bond_list(int natom) { - int i,j,ii,jj,m,n,inum,jnum,nbond; + int i,j,ii,jj,m,n,iold,jold,ilocal,jlocal,inum,jnum,nbond; tagint itag,jtag; - double xtmp,ytmp,ztmp,delx,dely,delz,rsq,oldcoeff; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq,distsq,oldcoeff; int *ilist,*jlist,*numneigh,**firstneigh; double time1,time2; @@ -782,80 +823,117 @@ void FixHyperLocal::build_bond_list(int natom) nevent_atom += natom; } + // compute max distance any bond atom has moved between 2 quenched states + // xold[iold] = last quenched coord for iold + // x[ilocal] = current quenched coord for same atom + // use of old2now calculates distsq only once per atom + + double **x = atom->x; + + for (i = 0; i < nall_old; i++) old2now[i] = -1; + + for (m = 0; m < nblocal; m++) { + iold = blist[m].iold; + if (old2now[iold] < 0) { + ilocal = atom->map(tagold[iold]); + ilocal = domain->closest_image(xold[iold],ilocal); + if (ilocal < 0) error->one(FLERR,"Fix hyper/local bond atom not found"); + old2now[iold] = ilocal; + distsq = MathExtra::distsq3(x[ilocal],xold[iold]); + maxdriftsq = MAX(distsq,maxdriftsq); + } + jold = blist[m].jold; + if (old2now[jold] < 0) { + jold = blist[m].jold; + jlocal = atom->map(tagold[jold]); + jlocal = domain->closest_image(xold[iold],jlocal); // close to I atom + if (jlocal < 0) error->one(FLERR,"Fix hyper/local bond atom not found"); + old2now[jold] = jlocal; + distsq = MathExtra::distsq3(x[jlocal],xold[jold]); + maxdriftsq = MAX(distsq,maxdriftsq); + } + } + + // store old bond coeffs so can persist them in new blist + // while loop allows growing value of maxbondperatom + // will loop at most 2 times, stops when maxbondperatom is large enough + // requires reverse comm, no forward comm: + // b/c new coeff list is stored only by current owned atoms + + tagint *tag = atom->tag; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; - // acquire old bond coeffs so can persist them in new blist - // while loop is to allow new value of maxcoeffperatom - // will loop at most 2 times, just once when maxcoeffperatom is large enough - // just reverse comm needed, - // b/c new bond list will be bonds of current owned atoms - - tagint *tag = atom->tag; - if (maxcoeff < nall) { + memory->destroy(numcoeff); + memory->sfree(clist); maxcoeff = atom->nmax; - grow_coeff(); + memory->create(numcoeff,maxcoeff,"hyper/local:numcoeff"); + clist = (OneCoeff **) memory->smalloc(maxcoeff*sizeof(OneCoeff *), + "hyper/local:clist"); } while (1) { if (firstflag) break; for (i = 0; i < nall; i++) numcoeff[i] = 0; + for (i = 0; i < nall; i++) clist[i] = NULL; + cpage->reset(); for (m = 0; m < nblocal; m++) { i = blist[m].i; j = blist[m].j; - if (numcoeff[i] < maxcoeffperatom) { - clist[i][numcoeff[i]].biascoeff = blist[m].biascoeff; - clist[i][numcoeff[i]].jtag = tag[j]; + if (numcoeff[i] == 0) clist[i] = cpage->get(maxbondperatom); + if (numcoeff[j] == 0) clist[j] = cpage->get(maxbondperatom); + + if (numcoeff[i] < maxbondperatom) { + clist[i][numcoeff[i]].biascoeff = biascoeff[m]; + clist[i][numcoeff[i]].tag = tag[j]; } numcoeff[i]++; - if (numcoeff[j] < maxcoeffperatom) { - clist[j][numcoeff[j]].biascoeff = blist[m].biascoeff; - clist[j][numcoeff[i]].jtag = tag[i]; + if (numcoeff[j] < maxbondperatom) { + clist[j][numcoeff[j]].biascoeff = biascoeff[m]; + clist[j][numcoeff[j]].tag = tag[i]; } numcoeff[j]++; } - int maxcol = 0; - for (i = 0; i < nall; i++) maxcol = MAX(maxcol,numcoeff[i]); - int maxcolall; - MPI_Allreduce(&maxcol,&maxcolall,1,MPI_INT,MPI_MAX,world); + int mymax = 0; + for (i = 0; i < nall; i++) mymax = MAX(mymax,numcoeff[i]); + int maxcoeffall; + MPI_Allreduce(&mymax,&maxcoeffall,1,MPI_INT,MPI_MAX,world); - if (maxcolall > maxcoeffperatom) { - maxcoeffperatom = maxcolall; - grow_coeff(); - memory->destroy(clist); - memory->create(clist,maxcoeff,maxcoeffperatom,"hyper/local:clist"); + if (maxcoeffall > maxbondperatom) { + maxbondperatom = maxcoeffall; + cpage->init(maxbondperatom,1024*maxbondperatom,1); continue; } commflag = BIASCOEFF; - comm->reverse_comm_fix(this); + comm->reverse_comm_fix_variable(this); - maxcol = 0; - for (i = 0; i < nall; i++) maxcol = MAX(maxcol,numcoeff[i]); - MPI_Allreduce(&maxcol,&maxcolall,1,MPI_INT,MPI_MAX,world); - if (maxcolall <= maxcoeffperatom) break; + mymax = 0; + for (i = 0; i < nall; i++) mymax = MAX(mymax,numcoeff[i]); + MPI_Allreduce(&mymax,&maxcoeffall,1,MPI_INT,MPI_MAX,world); + if (maxcoeffall <= maxbondperatom) break; - maxcoeffperatom = maxcolall; - grow_coeff(); + maxbondperatom = maxcoeffall; + cpage->init(maxbondperatom,1024*maxbondperatom,1); } - // reallocate vectors that are maxnew xold and tagold if necessary - // initialize xold to current coords - // initialize tagold to zero, so atoms not in neighbor list will remain zero + // reallocate vectors that are maxlocal and maxall length if necessary if (nlocal > maxlocal) { memory->destroy(eligible); memory->destroy(numbond); memory->destroy(maxhalf); + memory->destroy(maxhalfstrain); maxlocal = nlocal; memory->create(eligible,maxlocal,"hyper/local:eligible"); memory->create(numbond,maxlocal,"hyper/local:numbond"); memory->create(maxhalf,maxlocal,"hyper/local:maxhalf"); + memory->create(maxhalfstrain,maxlocal,"hyper/local:maxhalfstrain"); } if (nall > maxall) { @@ -870,23 +948,26 @@ void FixHyperLocal::build_bond_list(int natom) // nlocal_old = value of nlocal at time bonds are built // nall_old = value of nall at time bonds are built - // archive current peratom info in old vecs + // archive current atom coords in xold + // tagold will be set to non-zero below for accessed atoms + // numbond will be set below nlocal_old = nlocal; nall_old = nall; - double **x = atom->x; - memcpy(&xold[0][0],&x[0][0],3*nall*sizeof(double)); for (i = 0; i < nall; i++) tagold[i] = 0; for (i = 0; i < nlocal; i++) numbond[i] = 0; - // trigger builds for both neighbor lists - // NOTE: insure the I atoms are in same order? - + // trigger neighbor list builds for both lists + // insure the I loops in both are from 1 to nlocal + neighbor->build_one(listfull); neighbor->build_one(listhalf); + if (listfull->inum != nlocal || listhalf->inum != nlocal) + error->one(FLERR,"Invalid neighbor list in fix hyper/local bond build"); + // set tagold = 1 for all J atoms used in full neighbor list // tagold remains 0 for unused atoms, skipped in pre_neighbor @@ -897,12 +978,13 @@ void FixHyperLocal::build_bond_list(int natom) for (ii = 0; ii < inum; ii++) { i = ilist[ii]; + tagold[i] = tag[i]; jlist = firstneigh[i]; jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; j &= NEIGHMASK; - tagold[j] = 1; + tagold[j] = tag[j]; } } @@ -916,7 +998,6 @@ void FixHyperLocal::build_bond_list(int natom) numneigh = listhalf->numneigh; firstneigh = listhalf->firstneigh; - bigint bondcount = 0; nblocal = 0; for (ii = 0; ii < inum; ii++) { @@ -966,16 +1047,16 @@ void FixHyperLocal::build_bond_list(int natom) jtag = tag[j]; n = numcoeff[i]; for (m = 0; m < n; m++) { - if (clist[i][m].jtag == jtag) { + if (clist[i][m].tag == jtag) { oldcoeff = clist[i][m].biascoeff; break; } } } - if (oldcoeff > 0.0) blist[nblocal].biascoeff = oldcoeff; + if (oldcoeff > 0.0) biascoeff[nblocal] = oldcoeff; else { - blist[nblocal].biascoeff = COEFFINIT; + biascoeff[nblocal] = COEFFINIT; nnewbond++; } @@ -984,9 +1065,15 @@ void FixHyperLocal::build_bond_list(int natom) } numbond[i] = nbond; - bondcount += nbond; } - + + // this fix allows access to biascoeffs as local data + + size_local_rows = nblocal; + + // allbonds = total # of bonds in system + + bigint bondcount = nblocal; MPI_Allreduce(&bondcount,&allbonds,1,MPI_LMP_BIGINT,MPI_SUM,world); time2 = MPI_Wtime(); @@ -1010,6 +1097,7 @@ int FixHyperLocal::pack_forward_comm(int n, int *list, double *buf, // STRAIN // pack maxstrain vector + // must send to all ghosts out to Dcut if (commflag == STRAIN) { for (i = 0; i < n; i++) { @@ -1019,6 +1107,8 @@ int FixHyperLocal::pack_forward_comm(int n, int *list, double *buf, // STRAINDOMAIN // pack maxstrain_domain vector + // could just send to nearby ghosts in bonds + // don't see easy way to determine precisely which atoms that is } else if (commflag == STRAINDOMAIN) { for (i = 0; i < n; i++) { @@ -1028,6 +1118,7 @@ int FixHyperLocal::pack_forward_comm(int n, int *list, double *buf, // BIASFLAG // pack biasflag vector + // must send to all ghosts out to Dcut } else if (commflag == BIASFLAG) { for (i = 0; i < n; i++) { @@ -1085,22 +1176,37 @@ int FixHyperLocal::pack_reverse_comm(int n, int first, double *buf) // STRAIN // pack maxstrain vector + // only pack for nonzero values if (commflag == STRAIN) { + int nonzero = 0; + m++; // placeholder for count of atoms for (i = first; i < last; i++) { - buf[m++] = maxstrain[i]; + if (maxstrain[i] == 0.0) continue; + nonzero++; + buf[m++] = ubuf(i-first).d; // which atom is next + buf[m++] = maxstrain[i]; // value } + buf[0] = ubuf(nonzero).d; // STRAINDOMAIN // pack maxstrain_domain vector + // only pack for nonzero values } else if (commflag == STRAINDOMAIN) { + int nonzero = 0; + m++; // placeholder for count of atoms for (i = first; i < last; i++) { - buf[m++] = maxstrain_domain[i]; + if (maxstrain[i] == 0.0) continue; + nonzero++; + buf[m++] = ubuf(i-first).d; // which atom is next + buf[m++] = maxstrain_domain[i]; // value } + buf[0] = ubuf(nonzero).d; // BIASFLAG // pack biasflag vector + // could just pack for nonzero values, like STRAIN and STRAINDOMAIN } else if (commflag == BIASFLAG) { for (i = first; i < last; i++) { @@ -1109,22 +1215,46 @@ int FixHyperLocal::pack_reverse_comm(int n, int first, double *buf) // BIASCOEFF // pack list of biascoeffs + // only pack for atoms with nonzero # of bias coeffs + // this will skip majority of ghost atoms } else if (commflag == BIASCOEFF) { int ncoeff; + int nonzero = 0; + m++; // placeholder for count of atoms for (i = first; i < last; i++) { + if (numcoeff[i] == 0) continue; + nonzero++; ncoeff = numcoeff[i]; - buf[m++] = ubuf(ncoeff).d; + buf[m++] = ubuf(i-first).d; // which atom is next + buf[m++] = ubuf(ncoeff).d; // # of bias coeffs for (j = 0; j < ncoeff; j++) { buf[m++] = clist[i][j].biascoeff; - buf[m++] = ubuf(clist[i][j].jtag).d; + buf[m++] = ubuf(clist[i][j].tag).d; } } + buf[0] = ubuf(nonzero).d; } return m; } +/* ---------------------------------------------------------------------- + callback by comm->reverse_comm_fix_variable() in build_bond() + same logic as BIASCOEFF option in pack_reverse_comm() + m = returned size of message +------------------------------------------------------------------------- */ + +int FixHyperLocal::pack_reverse_comm_size(int n, int first) +{ + int last = first + n; + int m = 1; + for (int i = first; i < last; i++) { + if (numcoeff[i]) m += 2 + 2*numcoeff[i]; + } + return m; +} + /* ---------------------------------------------------------------------- */ void FixHyperLocal::unpack_reverse_comm(int n, int *list, double *buf) @@ -1135,23 +1265,32 @@ void FixHyperLocal::unpack_reverse_comm(int n, int *list, double *buf) // STRAIN // unpack maxstrain vector + // nonzero # of entries, each has offset to which atom in receiver's list // use MAX, b/c want maximum abs value strain for each atom's bonds if (commflag == STRAIN) { - for (i = 0; i < n; i++) { - j = list[i]; + int offset; + int nonzero = (int) ubuf(buf[m++]).i; // # of atoms with values + for (int iatom = 0; iatom < nonzero; iatom++) { + offset = (int) ubuf(buf[m++]).i; // offset into list for which atom + j = list[offset]; maxstrain[j] = MAX(maxstrain[j],buf[m]); m++; } // STRAINDOMAIN // unpack maxstrain_domain vector - // use SUM, b/c exactly one ghost or owned value per atom ID is non-zero + // use MAX, b/c want maximum abs value strain for each atom's domain + // could also use SUM, b/c exactly one ghost or owned value is non-zero } else if (commflag == STRAINDOMAIN) { - for (i = 0; i < n; i++) { - j = list[i]; - maxstrain_domain[j] += buf[m++]; + int offset; + int nonzero = (int) ubuf(buf[m++]).i; // # of atoms with values + for (int iatom = 0; iatom < nonzero; iatom++) { + offset = (int) ubuf(buf[m++]).i; // offset into list for which atom + j = list[offset]; + maxstrain_domain[j] = MAX(maxstrain_domain[j],buf[m]); + m++; } // BIASFLAG @@ -1164,19 +1303,23 @@ void FixHyperLocal::unpack_reverse_comm(int n, int *list, double *buf) } // BIASCOEFF - // unpack list of biascoeffs and add to atom J's list - // protect against overflow of clist columns - // if that happens, caller will realloc clist and reverse comm again + // unpack list of biascoeffs + // nonzero # of entries, each has offset to which atom in receiver's list + // protect against overflow of clist vector + // if that happens, caller will re-setup cpage and reverse comm again - } else if (commflag == BIASFLAG) { - int ncoeff; - for (i = 0; i < n; i++) { - j = list[i]; - ncoeff = (int) ubuf(buf[m++]).i; + } else if (commflag == BIASCOEFF) { + int offset,ncoeff; + int nonzero = (int) ubuf(buf[m++]).i; // # of atoms with coeffs + for (int iatom = 0; iatom < nonzero; iatom++) { + offset = (int) ubuf(buf[m++]).i; // offset into list for which atom + j = list[offset]; + ncoeff = (int) ubuf(buf[m++]).i; // # of bias coeffs for (k = 0; k < ncoeff; k++) { - if (numcoeff[j] < maxcoeffperatom) { + if (numcoeff[j] == 0) clist[j] = cpage->get(maxbondperatom); + if (numcoeff[j] < maxbondperatom) { clist[j][numcoeff[j]].biascoeff = buf[m++]; - clist[j][numcoeff[j]].jtag = (tagint) ubuf(buf[m++]).i; + clist[j][numcoeff[j]].tag = (tagint) ubuf(buf[m++]).i; } else m += 2; numcoeff[j]++; } @@ -1185,7 +1328,7 @@ void FixHyperLocal::unpack_reverse_comm(int n, int *list, double *buf) } /* ---------------------------------------------------------------------- - grow bond list by a chunk + grow bond list and bias coeff vector by a chunk ------------------------------------------------------------------------- */ void FixHyperLocal::grow_bond() @@ -1195,16 +1338,8 @@ void FixHyperLocal::grow_bond() maxbond += DELTABOND; blist = (OneBond *) memory->srealloc(blist,maxbond*sizeof(OneBond),"hyper/local:blist"); -} - -/* ---------------------------------------------------------------------- - reallocate 2-dimensional clist -------------------------------------------------------------------------- */ - -void FixHyperLocal::grow_coeff() -{ - memory->destroy(clist); - memory->create(clist,maxcoeff,maxcoeffperatom,"hyper/local:clist"); + memory->grow(biascoeff,maxbond,"hyper/local:biascoeff"); + vector_local = biascoeff; } /* ---------------------------------------------------------------------- */ @@ -1220,41 +1355,35 @@ double FixHyperLocal::compute_scalar() double FixHyperLocal::compute_vector(int i) { - // 23 vector outputs returned for i = 0-22 + // 21 vector outputs returned for i = 0-20 // i = 0 = # of biased bonds on this step // i = 1 = max strain of any bond on this step - // i = 2 = average bias potential for all bonds on this step + // i = 2 = average bias coeff for all bonds on this step // i = 3 = ave bonds/atom on this step // i = 4 = ave neighbor bonds/bond on this step - // i = 5 = fraction of steps and bonds with no bias during this run - // i = 6 = max drift distance of any atom during this run - // i = 7 = max bond length during this run - // i = 8 = average # of biased bonds/step during this run - // i = 9 = average bias potential for all bonds during this run - // i = 10 = max bias potential for any bond during this run - // i = 11 = min bias potential for any bond during this run - // i = 12 = max dist from my box of any ghost atom with - // maxstain < qfactor during this run - // i = 13 = max dist from my box of any ghost atom with + // i = 5 = max bond length during this run + // i = 6 = average # of biased bonds/step during this run + // i = 7 = fraction of biased bonds with no bias during this run + // i = 8 = fraction of biased bonds with negative strain during this run + // i = 9 = average bias coeff for all bonds during this run + // i = 10 = min bias coeff for any bond during this run + // i = 11 = max bias coeff for any bond during this run + + // i = 12 = max drift distance of any atom during this run + // i = 13 = max distance from proc subbox of any ghost atom with + // maxstrain < qfactor during this run + // i = 14 = max distance from proc subbox of any ghost atom with // any maxstrain during this run - // i = 14 = count of ghost atoms that could not be found - // by any proc at any reneighbor step during this run + // i = 15 = count of ghost atoms that could not be found + // on reneighbor steps during this run + // i = 16 = count of bias overlaps (< Dcut) found during this run - // NOTE: these 2 are no longer relevant - // i = 15 = count of lost bond partners during this run - // i = 16 = average bias coeff for lost bond partners during this run - - // i = 17 = count of bias overlaps found during this run - - // NOTE: this is no longer relevant - // i = 18 = count of non-matching bias coefficients found during this run - - // i = 19 = cummulative hyper time - // i = 20 = cummulative # of event timesteps since fix created - // i = 21 = cummulative # of atoms in events since fix created - // i = 22 = cummulative # of new bonds formed since fix created + // i = 17 = cumulative hyper time since fix created + // i = 18 = cumulative # of event timesteps since fix created + // i = 19 = cumulative # of atoms in events since fix created + // i = 20 = cumulative # of new bonds formed since fix created if (i == 0) { int nbiasall; @@ -1274,7 +1403,7 @@ double FixHyperLocal::compute_vector(int i) } if (i == 2) { - if (allbias && allbonds) return allbias/allbonds * vmax; + if (allbonds) return sumbiascoeff/allbonds; return 1.0; } @@ -1294,6 +1423,19 @@ double FixHyperLocal::compute_vector(int i) } if (i == 5) { + double allbondlen; + MPI_Allreduce(&maxbondlen,&allbondlen,1,MPI_DOUBLE,MPI_MAX,world); + return allbondlen; + } + + if (i == 6) { + if (update->ntimestep == update->firststep) return 0.0; + int allbias_running; + MPI_Allreduce(&nbias_running,&allbias_running,1,MPI_INT,MPI_SUM,world); + return 1.0*allbias_running / (update->ntimestep - update->firststep); + } + + if (i == 7) { int allbias_running,allnobias_running; MPI_Allreduce(&nbias_running,&allbias_running,1,MPI_INT,MPI_SUM,world); MPI_Allreduce(&nobias_running,&allnobias_running,1,MPI_INT,MPI_SUM,world); @@ -1301,68 +1443,64 @@ double FixHyperLocal::compute_vector(int i) return 0.0; } - if (i == 6) { + if (i == 8) { + int allbias_running,allnegstrain_running; + MPI_Allreduce(&nbias_running,&allbias_running,1,MPI_INT,MPI_SUM,world); + MPI_Allreduce(&negstrain_running,&allnegstrain_running,1,MPI_INT, + MPI_SUM,world); + if (allbias_running) return 1.0*allnegstrain_running / allbias_running; + return 0.0; + } + + if (i == 9) { + if (update->ntimestep == update->firststep) return 0.0; + return avebiascoeff / (update->ntimestep - update->firststep); + } + + if (i == 10) { + double biascoeff; + MPI_Allreduce(&minbiascoeff,&biascoeff,1,MPI_DOUBLE,MPI_MIN,world); + return biascoeff; + } + + if (i == 11) { + double biascoeff; + MPI_Allreduce(&maxbiascoeff,&biascoeff,1,MPI_DOUBLE,MPI_MAX,world); + return biascoeff; + } + + if (i == 12) { double alldriftsq; MPI_Allreduce(&maxdriftsq,&alldriftsq,1,MPI_DOUBLE,MPI_MAX,world); return (double) sqrt(alldriftsq); } - if (i == 7) { - double allbondlen; - MPI_Allreduce(&maxbondlen,&allbondlen,1,MPI_DOUBLE,MPI_MAX,world); - return allbondlen; - } + if (i == 13) return rmaxever; + if (i == 14) return rmaxeverbig; - if (i == 8) { - if (update->ntimestep == update->firststep) return 0.0; - int allbias_running; - MPI_Allreduce(&nbias_running,&allbias_running,1,MPI_INT,MPI_SUM,world); - return 1.0*allbias_running / (update->ntimestep - update->firststep); - } - - if (i == 9) { - if (update->ntimestep == update->firststep) return 0.0; - return sumbiascoeff * vmax / (update->ntimestep - update->firststep); - } - - if (i == 10) { - double allbiascoeff; - MPI_Allreduce(&maxbiascoeff,&allbiascoeff,1,MPI_DOUBLE,MPI_MAX,world); - return allbiascoeff * vmax; - } - - if (i == 11) { - double allbiascoeff; - MPI_Allreduce(&minbiascoeff,&allbiascoeff,1,MPI_DOUBLE,MPI_MAX,world); - return allbiascoeff * vmax; - } - - if (i == 12) return rmaxever; - if (i == 13) return rmaxeverbig; - - if (i == 14) { + if (i == 15) { int allghost_toofar; MPI_Allreduce(&ghost_toofar,&allghost_toofar,1,MPI_INT,MPI_SUM,world); return 1.0*allghost_toofar; } - if (i == 17) { + if (i == 16) { int allclose; MPI_Allreduce(&checkbias_count,&allclose,1,MPI_INT,MPI_SUM,world); return 1.0*allclose; } - if (i == 19) { + if (i == 17) { return boost_target * update->dt * (update->ntimestep - starttime); } - if (i == 20) return (double) nevent; - if (i == 21) return (double) nevent_atom; + if (i == 18) return (double) nevent; + if (i == 19) return (double) nevent_atom; - if (i == 22) { - int allnew; - MPI_Allreduce(&nnewbond,&allnew,1,MPI_INT,MPI_SUM,world); - return (double) 0.5*allnew; + if (i == 20) { + int allnewbond; + MPI_Allreduce(&nnewbond,&allnewbond,1,MPI_INT,MPI_SUM,world); + return (double) allnewbond; } return 0.0; @@ -1375,32 +1513,30 @@ double FixHyperLocal::compute_vector(int i) double FixHyperLocal::query(int i) { - if (i == 1) return compute_vector(19); // cummulative hyper time - if (i == 2) return compute_vector(20); // nevent - if (i == 3) return compute_vector(21); // nevent_atom + if (i == 1) return compute_vector(17); // cummulative hyper time + if (i == 2) return compute_vector(18); // nevent + if (i == 3) return compute_vector(19); // nevent_atom if (i == 4) return compute_vector(3); // ave bonds/atom - if (i == 5) return compute_vector(6); // maxdrift - if (i == 6) return compute_vector(7); // maxbondlen - if (i == 7) return compute_vector(5); // fraction with zero bias + if (i == 5) return compute_vector(12); // maxdrift + if (i == 6) return compute_vector(5); // maxbondlen + if (i == 7) return compute_vector(7); // fraction with zero bias + if (i == 8) return compute_vector(8); // fraction with negative strain // unique to local hyper - if (i == 8) return compute_vector(22); // number of new bonds - if (i == 9) return 1.0*maxbondperatom; // max bonds/atom - if (i == 10) return compute_vector(8); // ave # of biased bonds/step - if (i == 11) return compute_vector(9); // ave bias coeff over all bonds - if (i == 12) return compute_vector(10); // max bias cooef for any bond - if (i == 13) return compute_vector(11); // max bias cooef for any bond - if (i == 14) return compute_vector(4); // neighbor bonds/bond - if (i == 15) return compute_vector(2); // ave bias coeff now - if (i == 16) return time_bondbuild; // CPU time for bond_build calls - if (i == 17) return rmaxever; // ghost atom distance for < maxstrain - if (i == 18) return rmaxeverbig; // ghost atom distance for any strain - if (i == 19) return compute_vector(14); // count of ghost atoms not found - //if (i == 20) return compute_vector(15); // count of lost bond partners - //if (i == 21) return compute_vector(16); // ave bias coeff of long bonds - if (i == 22) return compute_vector(17); // count of bias overlaps - //if (i == 23) return compute_vector(18); // count of non-matching bias coeffs + if (i == 9) return compute_vector(20); // number of new bonds + if (i == 10) return 1.0*maxbondperatom; // max bonds/atom + if (i == 11) return compute_vector(6); // ave # of biased bonds/step + if (i == 12) return compute_vector(9); // ave bias coeff over all bonds + if (i == 13) return compute_vector(10); // min bias cooef for any bond + if (i == 14) return compute_vector(11); // max bias cooef for any bond + if (i == 15) return compute_vector(4); // neighbor bonds/bond + if (i == 16) return compute_vector(2); // ave bias coeff now + if (i == 17) return time_bondbuild; // CPU time for build_bond calls + if (i == 18) return rmaxever; // ghost atom distance for < maxstrain + if (i == 19) return rmaxeverbig; // ghost atom distance for any strain + if (i == 20) return compute_vector(15); // count of ghost atoms not found + if (i == 21) return compute_vector(16); // count of bias overlaps error->all(FLERR,"Invalid query to fix hyper/local"); @@ -1408,21 +1544,22 @@ double FixHyperLocal::query(int i) } /* ---------------------------------------------------------------------- - memory usage of per-atom and per-bond arrays + memory usage of per-atom and per-bond data structs ------------------------------------------------------------------------- */ double FixHyperLocal::memory_usage() { - int nmax = atom->nmax; - double bytes = maxbond * sizeof(OneBond); // bond list + double bytes = maxbond * sizeof(OneBond); // blist + bytes = maxbond * sizeof(double); // per-bond bias coeffs bytes += 3*maxlocal * sizeof(int); // numbond,maxhalf,eligible bytes += maxlocal * sizeof(double); // maxhalfstrain bytes += maxall * sizeof(int); // old2now bytes += maxall * sizeof(tagint); // tagold bytes += 3*maxall * sizeof(double); // xold - bytes += 2*nmax * sizeof(double); // maxstrain,maxstrain_domain - if (checkbias) bytes += nmax * sizeof(tagint); // biasflag - bytes += maxcoeff*maxcoeffperatom * sizeof(OneCoeff); // clist - bytes += maxcoeff * sizeof(int); // numcoeff + bytes += 2*maxall * sizeof(double); // maxstrain,maxstrain_domain + if (checkbias) bytes += maxall * sizeof(tagint); // biasflag + bytes += maxcoeff * sizeof(int); // numcoeff + bytes += maxcoeff * sizeof(OneCoeff *); // clist + bytes += maxlocal*maxbondperatom * sizeof(OneCoeff); // cpage estimate return bytes; } diff --git a/src/REPLICA/fix_hyper_local.h b/src/REPLICA/fix_hyper_local.h index 67361ce3ac..f0075c185c 100644 --- a/src/REPLICA/fix_hyper_local.h +++ b/src/REPLICA/fix_hyper_local.h @@ -21,6 +21,7 @@ FixStyle(hyper/local,FixHyperLocal) #define LMP_FIX_HYPER_LOCAL_H #include "fix_hyper.h" +#include "my_page.h" namespace LAMMPS_NS { @@ -43,6 +44,7 @@ class FixHyperLocal : public FixHyper { int pack_forward_comm(int, int *, double *, int, int *); void unpack_forward_comm(int, int, double *); int pack_reverse_comm(int, int, double *); + int pack_reverse_comm_size(int, int); void unpack_reverse_comm(int, int *, double *); double memory_usage(); @@ -54,54 +56,72 @@ class FixHyperLocal : public FixHyper { private: int me; + + // inputs and derived quantities + double cutbond,qfactor,vmax,tequil,dcut; double alpha_user; // timescale to apply boostostat (time units) double alpha; // unitless dt/alpha_user double boost_target; // target value of boost - int checkbias,checkbias_every,checkbias_flag,checkbias_count; + int checkghost,checkbias; // flags for optional stats + + double cutbondsq,dcutsq; + double beta,invqfactorsq; + + // flags int setupflag; // 1 during setup, 0 during run int firstflag; // set for first time bond_build takes place - int nostrainyet; // 1 until maxstrain is first computed - - int nbias_running,nobias_running; - int nbondbuild; - double time_bondbuild; - bigint starttime; - double sumbiascoeff; // sum of aveboost at every timestep - bigint allbonds; // sum of bond count on this step - double allbias; // sum of biascoeff on all bonds on this step - - int nnewbond; // running tally of number of new bonds created - int maxbondperatom; // max # of bonds any atom ever has + int nostrainyet; // 1 until maxstrain is first compute + bigint starttime; // timestep when this fix was invoked int commflag; // flag for communication mode + + // stats + + int nbondbuild; // # of rebuilds of bond list + double time_bondbuild; // CPU time for bond builds + + bigint allbonds; // current total # of bonds + int nnewbond; // running tally of # of new bonds created + int maxbondperatom; // max # of bonds any atom ever has int nevent; // # of events that trigger bond rebuild int nevent_atom; // # of atoms that experienced an event - double cutbondsq,dcutsq; - double beta,invqfactorsq; - double mybias; + + int nbias_running; // running count of biased bonds + int nobias_running; // ditto for bonds with bias = 0, b/c too long + int negstrain_running; // ditto for bonds with negative strain + + double mybias; // sum of bias potentials for biased bonds double maxbondlen; // cummulative max length of any bond - double maxdriftsq; // max distance any atom drifts from original pos - double maxbiascoeff; // cummulative max bias coeff for any bond + double maxdriftsq; // max distance any bond atom drifts from quenched x + + double sumbiascoeff; // sum of all bond bias coeffs at each timestep + double avebiascoeff; // cummulative sumbiascoeff/allbonds across steps double minbiascoeff; // cummulative min bias coeff for any bond + double maxbiascoeff; // cummulative max bias coeff for any bond + double rmaxever,rmaxeverbig; - int ghost_toofar; + int ghost_toofar; // # of ghost atoms not found in Dcut neigh list + + int checkbias_every,checkbias_flag,checkbias_count; + + // 2 neighbor lists class NeighList *listfull; // full neigh list up to Dcut distance class NeighList *listhalf; // half neigh list up to pair distance // both created only when bonds are rebuilt - // list of my owned bonds + // list of my owned bonds and bias coeffs // persists on a proc from one event until the next struct OneBond { // single IJ bond, atom I is owner int i,j; // current local indices of 2 bond atoms int iold,jold; // local indices when bonds were formed double r0; // relaxed bond length - double biascoeff; // biasing coefficient = prefactor Cij }; - struct OneBond *blist; // list of owned bonds + OneBond *blist; // list of owned bonds + double *biascoeff; // biasing coefficient Cij for each bond int nblocal; // # of owned bonds int maxbond; // allocated size of blist @@ -137,24 +157,24 @@ class FixHyperLocal : public FixHyper { tagint *biasflag; // atoms in biased bonds marked with bond partner // for owned and ghost atoms - // data struct used to persist biascoeffs when bond list is re-created - - struct OneCoeff { - double biascoeff; - tagint jtag; - }; - - struct OneCoeff **clist; // list of bond coeffs for each atom's bonds - int *numcoeff; // # of coeffs per atom - int maxcoeff; // allocate size of clist - int maxcoeffperatom; // allocated # of columns in clist - // list of biased bonds this proc owns int maxbias; // allocated size of bias list int nbias; // # of biased bonds I own int *bias; // index of biased bonds in my bond list + // data structs for persisting bias coeffs when bond list is reformed + + struct OneCoeff { + double biascoeff; + tagint tag; + }; + + MyPage *cpage; // pages of OneCoeff datums for clist + OneCoeff **clist; // ptrs to vectors of bias coeffs for each atom + int *numcoeff; // # of bias coeffs per atom (one per bond) + int maxcoeff; // allocate sized of clist and numcoeff + // extra timers //double timefirst,timesecond,timethird,timefourth; @@ -163,7 +183,6 @@ class FixHyperLocal : public FixHyper { // private methods void grow_bond(); - void grow_coeff(); }; } diff --git a/src/REPLICA/hyper.cpp b/src/REPLICA/hyper.cpp index 22940de221..0d8de6d060 100644 --- a/src/REPLICA/hyper.cpp +++ b/src/REPLICA/hyper.cpp @@ -40,9 +40,7 @@ enum{NOHYPER,GLOBAL,LOCAL}; /* ---------------------------------------------------------------------- */ -Hyper::Hyper(LAMMPS *lmp) : - Pointers(lmp), dumplist(NULL) -{} +Hyper::Hyper(LAMMPS *lmp) : Pointers(lmp), dumplist(NULL) {} /* ---------------------------------------------------------------------- perform hyperdynamics simulation @@ -260,11 +258,12 @@ void Hyper::command(int narg, char **arg) double maxdrift = 0.0; double maxbondlen = 0.0; double fraczero = 1.0; + double fracneg = 1.0; - double nnewbond,avenboost,aveboostcoeff,maxboostcoeff,minboostcoeff; - double maxbondperatom,neighbondperbond,aveboostnow; + double nnewbond,avenbias,avebiascoeff,minbiascoeff,maxbiascoeff; + double maxbondperatom,neighbondperbond,avebiasnow; double tbondbuild,rmaxever,rmaxeverbig,allghost_toofar; - double lostbond,lostbondcoeff,biasoverlap,nonmatchbiascoeff; + double lostbond,lostbondcoeff,biasoverlap; if (hyperenable) { t_hyper = fix_hyper->query(1); @@ -274,115 +273,70 @@ void Hyper::command(int narg, char **arg) maxdrift = fix_hyper->query(5); maxbondlen = fix_hyper->query(6); fraczero = fix_hyper->query(7); + fracneg = fix_hyper->query(8); if (hyperstyle == LOCAL) { - nnewbond = fix_hyper->query(8); - maxbondperatom = fix_hyper->query(9); - avenboost = fix_hyper->query(10); - aveboostcoeff = fix_hyper->query(11); - maxboostcoeff = fix_hyper->query(12); - minboostcoeff = fix_hyper->query(13); - neighbondperbond = fix_hyper->query(14); - aveboostnow = fix_hyper->query(15); - tbondbuild = fix_hyper->query(16); - rmaxever = fix_hyper->query(17); - rmaxeverbig = fix_hyper->query(18); - allghost_toofar = fix_hyper->query(19); - lostbond = fix_hyper->query(20); - lostbondcoeff = fix_hyper->query(21); - biasoverlap = fix_hyper->query(22); - nonmatchbiascoeff = fix_hyper->query(23); + nnewbond = fix_hyper->query(9); + maxbondperatom = fix_hyper->query(10); + avenbias = fix_hyper->query(11); + avebiascoeff = fix_hyper->query(12); + minbiascoeff = fix_hyper->query(13); + maxbiascoeff = fix_hyper->query(14); + neighbondperbond = fix_hyper->query(15); + avebiasnow = fix_hyper->query(16); + tbondbuild = fix_hyper->query(17); + rmaxever = fix_hyper->query(18); + rmaxeverbig = fix_hyper->query(19); + allghost_toofar = fix_hyper->query(20); + biasoverlap = fix_hyper->query(21); } } if (me == 0) { - if (screen) { - fprintf(screen,"Cummulative quantities for fix hyper:\n"); - fprintf(screen," hyper time = %g\n",t_hyper); - fprintf(screen," time boost factor = %g\n",t_hyper/(nsteps*update->dt)); - fprintf(screen," event timesteps = %d\n",nevent_running); - fprintf(screen," # of atoms in events = %d\n",nevent_atoms_running); - fprintf(screen,"Quantities for this hyper run:\n"); - fprintf(screen," event timesteps = %d\n",nevent); - fprintf(screen," # of atoms in events = %d\n",nevent_atoms); - fprintf(screen," max length of any bond = %g\n",maxbondlen); - fprintf(screen," max drift distance of any atom = %g\n",maxdrift); - fprintf(screen," fraction of steps & bonds with zero bias = %g\n", - fraczero); - fprintf(screen,"Current quantities:\n"); - fprintf(screen," ave bonds/atom = %g\n",avebonds); + FILE *out; + for (int iout = 0; iout < 2; iout++) { + if (iout == 0) out = screen; + if (iout == 1) out = logfile; + if (!out) continue; + fprintf(out,"Cummulative quantities for fix hyper:\n"); + fprintf(out," hyper time = %g\n",t_hyper); + fprintf(out," time boost factor = %g\n",t_hyper/(nsteps*update->dt)); + fprintf(out," event timesteps = %d\n",nevent_running); + fprintf(out," # of atoms in events = %d\n",nevent_atoms_running); + fprintf(out,"Quantities for this hyper run:\n"); + fprintf(out," event timesteps = %d\n",nevent); + fprintf(out," # of atoms in events = %d\n",nevent_atoms); + fprintf(out," max length of any bond = %g\n",maxbondlen); + fprintf(out," max drift distance of any atom = %g\n",maxdrift); + fprintf(out," fraction of biased bonds with zero bias = %g\n",fraczero); + fprintf(out," fraction of biased bonds with negative strain = %g\n", + fracneg); + fprintf(out,"Current quantities:\n"); + fprintf(out," ave bonds/atom = %g\n",avebonds); if (hyperstyle == LOCAL) { - fprintf(screen,"Cummulative quantities specific to fix hyper/local:\n"); - fprintf(screen," # of new bonds formed = %g\n",nnewbond); - fprintf(screen," max bonds/atom = %g\n",maxbondperatom); - fprintf(screen,"Quantities for this hyper run specific to " + fprintf(out,"Cummulative quantities specific to fix hyper/local:\n"); + fprintf(out," # of new bonds formed = %g\n",nnewbond); + fprintf(out," max bonds/atom = %g\n",maxbondperatom); + fprintf(out,"Quantities for this hyper run specific to " "fix hyper/local:\n"); - fprintf(screen," ave boosted bonds/step = %g\n",avenboost); - fprintf(screen," ave boost coeff of all bonds = %g\n",aveboostcoeff); - fprintf(screen," max boost coeff of any bond = %g\n",maxboostcoeff); - fprintf(screen," min boost coeff of any bond = %g\n",minboostcoeff); - fprintf(screen," max dist from my box of any " + fprintf(out," ave biased bonds/step = %g\n",avenbias); + fprintf(out," ave bias coeff of all bonds = %g\n",avebiascoeff); + fprintf(out," min bias coeff of any bond = %g\n",minbiascoeff); + fprintf(out," max bias coeff of any bond = %g\n",maxbiascoeff); + fprintf(out," max dist from my subbox of any " "non-maxstrain bond ghost atom = %g\n",rmaxever); - fprintf(screen," max dist from my box of any bond ghost atom = %g\n", + fprintf(out," max dist from my box of any bond ghost atom = %g\n", rmaxeverbig); - fprintf(screen," count of bond ghost neighbors " + fprintf(out," count of bond ghost neighbors " "not found on reneighbor steps = %g\n",allghost_toofar); - fprintf(screen," lost bond partners = %g\n",lostbond); - fprintf(screen," ave bias coeff for lost bond partners = %g\n", - lostbondcoeff); - fprintf(screen," bias overlaps = %g\n",biasoverlap); - fprintf(screen," non-matching bias coeffs = %g\n",nonmatchbiascoeff); - fprintf(screen," CPU time for bond builds = %g\n",tbondbuild); - fprintf(screen,"Current quantities specific to fix hyper/local:\n"); - fprintf(screen," neighbor bonds/bond = %g\n",neighbondperbond); - fprintf(screen," ave boost coeff for all bonds = %g\n",aveboostnow); + fprintf(out," bias overlaps = %g\n",biasoverlap); + fprintf(out," CPU time for bond builds = %g\n",tbondbuild); + fprintf(out,"Current quantities specific to fix hyper/local:\n"); + fprintf(out," neighbor bonds/bond = %g\n",neighbondperbond); + fprintf(out," ave boost coeff for all bonds = %g\n",avebiasnow); } - fprintf(screen,"\n"); - } - - if (logfile) { - fprintf(logfile,"Cummulative quantities for fix hyper:\n"); - fprintf(logfile," hyper time = %g\n",t_hyper); - fprintf(logfile," event timesteps = %d\n",nevent_running); - fprintf(logfile," # of atoms in events = %d\n",nevent_atoms_running); - fprintf(logfile,"Quantities for this hyper run:\n"); - fprintf(logfile," event timesteps = %d\n",nevent); - fprintf(logfile," # of atoms in events = %d\n",nevent_atoms); - fprintf(logfile," max length of any bond = %g\n",maxbondlen); - fprintf(logfile," max drift distance of any atom = %g\n",maxdrift); - fprintf(logfile," fraction of steps & bonds with zero bias = %g\n", - fraczero); - fprintf(logfile,"Current quantities:\n"); - fprintf(logfile," ave bonds/atom = %g\n",avebonds); - - if (hyperstyle == LOCAL) { - fprintf(logfile,"Cummulative quantities specific tofix hyper/local:\n"); - fprintf(logfile," # of new bonds formed = %g\n",nnewbond); - fprintf(logfile," max bonds/atom = %g\n",maxbondperatom); - fprintf(logfile,"Quantities for this hyper run specific to " - "fix hyper/local:\n"); - fprintf(logfile," ave boosted bonds/step = %g\n",avenboost); - fprintf(logfile," ave boost coeff of all bonds = %g\n",aveboostcoeff); - fprintf(logfile," max boost coeff of any bond = %g\n",maxboostcoeff); - fprintf(logfile," min boost coeff of any bond = %g\n",minboostcoeff); - fprintf(logfile," max dist from my box of any " - "non-maxstrain bond ghost atom = %g\n",rmaxever); - fprintf(logfile," max dist from my box of any bond ghost atom = %g\n", - rmaxeverbig); - fprintf(logfile," count of ghost bond neighbors " - "not found on reneighbor steps = %g\n",allghost_toofar); - fprintf(logfile," lost bond partners = %g\n",lostbond); - fprintf(logfile," ave bias coeff for lost bond partners = %g\n", - lostbondcoeff); - fprintf(logfile," bias overlaps = %g\n",biasoverlap); - fprintf(logfile," non-matching bias coeffs = %g\n",nonmatchbiascoeff); - fprintf(logfile," CPU time for bond builds = %g\n",tbondbuild); - fprintf(logfile,"Current quantities specific to fix hyper/local:\n"); - fprintf(logfile," neighbor bonds/bond = %g\n",neighbondperbond); - fprintf(logfile," ave boost coeff for all bonds = %g\n",aveboostnow); - } - fprintf(logfile,"\n"); + fprintf(out,"\n"); } } diff --git a/src/fix_ave_histo.cpp b/src/fix_ave_histo.cpp index a5bf8db557..d60fe7af14 100644 --- a/src/fix_ave_histo.cpp +++ b/src/fix_ave_histo.cpp @@ -32,7 +32,7 @@ using namespace FixConst; enum{X,V,F,COMPUTE,FIX,VARIABLE}; enum{ONE,RUNNING}; enum{SCALAR,VECTOR,WINDOW}; -enum{GLOBAL,PERATOM,LOCAL}; +enum{DEFAULT,GLOBAL,PERATOM,LOCAL}; enum{IGNORE,END,EXTRA}; #define INVOKED_SCALAR 1 @@ -46,8 +46,10 @@ enum{IGNORE,END,EXTRA}; FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), - nvalues(0), which(NULL), argindex(NULL), value2index(NULL), ids(NULL), fp(NULL), stats_list(NULL), - bin(NULL), bin_total(NULL), bin_all(NULL), bin_list(NULL), coord(NULL), vector(NULL) + nvalues(0), which(NULL), argindex(NULL), value2index(NULL), + ids(NULL), fp(NULL), stats_list(NULL), + bin(NULL), bin_total(NULL), bin_all(NULL), bin_list(NULL), + coord(NULL), vector(NULL) { if (narg < 10) error->all(FLERR,"Illegal fix ave/histo command"); @@ -188,9 +190,8 @@ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) : memory->sfree(earg); } - // setup and error check - // kind = inputs are all global, or all per-atom, or all local - // for fix inputs, check that fix frequency is acceptable + // check input args for kind consistency + // all inputs must all be global, per-atom, or local if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0) error->all(FLERR,"Illegal fix ave/histo command"); @@ -201,40 +202,65 @@ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) : if (ave != RUNNING && overwrite) error->all(FLERR,"Illegal fix ave/histo command"); - int kindflag; + int kindglobal,kindperatom,kindlocal; + for (int i = 0; i < nvalues; i++) { - if (which[i] == X || which[i] == V || which[i] == F) kindflag = PERATOM; - else if (which[i] == COMPUTE) { + kindglobal = kindperatom = kindlocal = 0; + + if (which[i] == X || which[i] == V || which[i] == F) { + kindperatom = 1; + + } else if (which[i] == COMPUTE) { int c_id = modify->find_compute(ids[i]); if (c_id < 0) error->all(FLERR,"Fix ave/histo input is invalid compute"); Compute *compute = modify->compute[c_id]; + // computes can produce multiple kinds of output if (compute->scalar_flag || compute->vector_flag || compute->array_flag) - kindflag = GLOBAL; - else if (compute->peratom_flag) kindflag = PERATOM; - else if (compute->local_flag) kindflag = LOCAL; - else error->all(FLERR,"Fix ave/histo input is invalid compute"); + kindglobal = 1; + if (compute->peratom_flag) kindperatom = 1; + if (compute->local_flag) kindlocal = 1; + } else if (which[i] == FIX) { int f_id = modify->find_fix(ids[i]); if (f_id < 0) error->all(FLERR,"Fix ave/histo input is invalid fix"); Fix *fix = modify->fix[f_id]; + // fixes can produce multiple kinds of output if (fix->scalar_flag || fix->vector_flag || fix->array_flag) - kindflag = GLOBAL; - else if (fix->peratom_flag) kindflag = PERATOM; - else if (fix->local_flag) kindflag = LOCAL; - else error->all(FLERR,"Fix ave/histo input is invalid fix"); + kindglobal = 1; + if (fix->peratom_flag) kindperatom = 1; + if (fix->local_flag) kindlocal = 1; + } else if (which[i] == VARIABLE) { int ivariable = input->variable->find(ids[i]); - if (ivariable < 0) error->all(FLERR,"Fix ave/histo input is invalid variable"); - if (input->variable->equalstyle(ivariable)) kindflag = GLOBAL; - else if (input->variable->atomstyle(ivariable)) kindflag = PERATOM; - else error->all(FLERR,"Fix ave/histo input is invalid variable"); + if (ivariable < 0) + error->all(FLERR,"Fix ave/histo input is invalid variable"); + // variables only produce one kind of output + if (input->variable->equalstyle(ivariable)) kindglobal = 1; + else if (input->variable->atomstyle(ivariable)) kindperatom = 1; + else error->all(FLERR,"Fix ave/histo input is invalid kind of variable"); + } + + if (kind == DEFAULT) { + if (kindglobal + kindperatom + kindlocal > 1) + error->all(FLERR,"Fix ave/histo input kind is ambiguous"); + if (kindglobal) kind = GLOBAL; + if (kindperatom) kind = PERATOM; + if (kindlocal) kind = LOCAL; + } else if (kind == GLOBAL) { + if (!kindglobal) + error->all(FLERR,"Fix ave/histo input kind is invalid"); + } else if (kind == PERATOM) { + if (!kindperatom) + error->all(FLERR,"Fix ave/histo input kind is invalid"); + } else if (kind == LOCAL) { + if (!kindlocal) + error->all(FLERR,"Fix ave/histo input kind is invalid"); } - if (i == 0) kind = kindflag; - else if (kindflag != kind) - error->all(FLERR, - "Fix ave/histo inputs are not all global, peratom, or local"); } + // more error checks + // for fix inputs, check that fix frequency is acceptable + if (kind == PERATOM && mode == SCALAR) error->all(FLERR, "Fix ave/histo cannot input per-atom values in scalar mode"); @@ -919,6 +945,7 @@ void FixAveHisto::options(int iarg, int narg, char **arg) // option defaults fp = NULL; + kind = DEFAULT; ave = ONE; startstep = 0; mode = SCALAR; @@ -942,6 +969,13 @@ void FixAveHisto::options(int iarg, int narg, char **arg) } } iarg += 2; + } else if (strcmp(arg[iarg],"kind") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/histo command"); + if (strcmp(arg[iarg+1],"global") == 0) kind = GLOBAL; + else if (strcmp(arg[iarg+1],"peratom") == 0) kind = PERATOM; + else if (strcmp(arg[iarg+1],"local") == 0) kind = LOCAL; + else error->all(FLERR,"Illegal fix ave/histo command"); + iarg += 2; } else if (strcmp(arg[iarg],"ave") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/histo command"); if (strcmp(arg[iarg+1],"one") == 0) ave = ONE; diff --git a/src/thermo.cpp b/src/thermo.cpp index ddbbd0f496..3e777edf82 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -380,9 +380,8 @@ void Thermo::compute(int flag) loc += sprintf(&line[loc],format[ifield],dvalue); else if (vtype[ifield] == INT) loc += sprintf(&line[loc],format[ifield],ivalue); - else if (vtype[ifield] == BIGINT) { + else if (vtype[ifield] == BIGINT) loc += sprintf(&line[loc],format[ifield],bivalue); - } } // print line to screen and logfile