Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80121949
fix_ave_chunk.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Aug 30, 00:24
Size
35 KB
Mime Type
text/x-c
Expires
Sun, Sep 1, 00:24 (2 d)
Engine
blob
Format
Raw Data
Handle
20301216
Attached To
rLAMMPS lammps
fix_ave_chunk.cpp
View Options
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include "fix_ave_chunk.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "domain.h"
#include "modify.h"
#include "compute.h"
#include "compute_chunk_atom.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"
using
namespace
LAMMPS_NS
;
using
namespace
FixConst
;
enum
{
V
,
F
,
DENSITY_NUMBER
,
DENSITY_MASS
,
MASS
,
TEMPERATURE
,
COMPUTE
,
FIX
,
VARIABLE
};
enum
{
SAMPLE
,
ALL
};
enum
{
NOSCALE
,
ATOM
};
enum
{
ONE
,
RUNNING
,
WINDOW
};
#define INVOKED_PERATOM 8
/* ---------------------------------------------------------------------- */
FixAveChunk
::
FixAveChunk
(
LAMMPS
*
lmp
,
int
narg
,
char
**
arg
)
:
Fix
(
lmp
,
narg
,
arg
),
nvalues
(
0
),
nrepeat
(
0
),
which
(
NULL
),
argindex
(
NULL
),
value2index
(
NULL
),
ids
(
NULL
),
fp
(
NULL
),
idchunk
(
NULL
),
varatom
(
NULL
),
count_one
(
NULL
),
count_many
(
NULL
),
count_sum
(
NULL
),
values_one
(
NULL
),
values_many
(
NULL
),
values_sum
(
NULL
),
count_total
(
NULL
),
count_list
(
NULL
),
values_total
(
NULL
),
values_list
(
NULL
)
{
if
(
narg
<
7
)
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
MPI_Comm_rank
(
world
,
&
me
);
nevery
=
force
->
inumeric
(
FLERR
,
arg
[
3
]);
nrepeat
=
force
->
inumeric
(
FLERR
,
arg
[
4
]);
nfreq
=
force
->
inumeric
(
FLERR
,
arg
[
5
]);
int
n
=
strlen
(
arg
[
6
])
+
1
;
idchunk
=
new
char
[
n
];
strcpy
(
idchunk
,
arg
[
6
]);
global_freq
=
nfreq
;
no_change_box
=
1
;
// expand args if any have wildcard character "*"
int
expand
=
0
;
char
**
earg
;
int
nargnew
=
input
->
expand_args
(
narg
-
7
,
&
arg
[
7
],
1
,
earg
);
if
(
earg
!=
&
arg
[
7
])
expand
=
1
;
arg
=
earg
;
// parse values until one isn't recognized
which
=
new
int
[
nargnew
];
argindex
=
new
int
[
nargnew
];
ids
=
new
char
*
[
nargnew
];
value2index
=
new
int
[
nargnew
];
int
iarg
=
0
;
while
(
iarg
<
nargnew
)
{
ids
[
nvalues
]
=
NULL
;
if
(
strcmp
(
arg
[
iarg
],
"vx"
)
==
0
)
{
which
[
nvalues
]
=
V
;
argindex
[
nvalues
++
]
=
0
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"vy"
)
==
0
)
{
which
[
nvalues
]
=
V
;
argindex
[
nvalues
++
]
=
1
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"vz"
)
==
0
)
{
which
[
nvalues
]
=
V
;
argindex
[
nvalues
++
]
=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"fx"
)
==
0
)
{
which
[
nvalues
]
=
F
;
argindex
[
nvalues
++
]
=
0
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"fy"
)
==
0
)
{
which
[
nvalues
]
=
F
;
argindex
[
nvalues
++
]
=
1
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"fz"
)
==
0
)
{
which
[
nvalues
]
=
F
;
argindex
[
nvalues
++
]
=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"density/number"
)
==
0
)
{
which
[
nvalues
]
=
DENSITY_NUMBER
;
argindex
[
nvalues
++
]
=
0
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"density/mass"
)
==
0
)
{
which
[
nvalues
]
=
DENSITY_MASS
;
argindex
[
nvalues
++
]
=
0
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"mass"
)
==
0
)
{
which
[
nvalues
]
=
MASS
;
argindex
[
nvalues
++
]
=
0
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"temp"
)
==
0
)
{
which
[
nvalues
]
=
TEMPERATURE
;
argindex
[
nvalues
++
]
=
0
;
}
else
if
(
strncmp
(
arg
[
iarg
],
"c_"
,
2
)
==
0
||
strncmp
(
arg
[
iarg
],
"f_"
,
2
)
==
0
||
strncmp
(
arg
[
iarg
],
"v_"
,
2
)
==
0
)
{
if
(
arg
[
iarg
][
0
]
==
'c'
)
which
[
nvalues
]
=
COMPUTE
;
else
if
(
arg
[
iarg
][
0
]
==
'f'
)
which
[
nvalues
]
=
FIX
;
else
if
(
arg
[
iarg
][
0
]
==
'v'
)
which
[
nvalues
]
=
VARIABLE
;
int
n
=
strlen
(
arg
[
iarg
]);
char
*
suffix
=
new
char
[
n
];
strcpy
(
suffix
,
&
arg
[
iarg
][
2
]);
char
*
ptr
=
strchr
(
suffix
,
'['
);
if
(
ptr
)
{
if
(
suffix
[
strlen
(
suffix
)
-
1
]
!=
']'
)
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
argindex
[
nvalues
]
=
atoi
(
ptr
+
1
);
*
ptr
=
'\0'
;
}
else
argindex
[
nvalues
]
=
0
;
n
=
strlen
(
suffix
)
+
1
;
ids
[
nvalues
]
=
new
char
[
n
];
strcpy
(
ids
[
nvalues
],
suffix
);
nvalues
++
;
delete
[]
suffix
;
}
else
break
;
iarg
++
;
}
if
(
nvalues
==
0
)
error
->
all
(
FLERR
,
"No values in fix ave/chunk command"
);
// optional args
normflag
=
ALL
;
scaleflag
=
ATOM
;
ave
=
ONE
;
nwindow
=
0
;
biasflag
=
0
;
id_bias
=
NULL
;
adof
=
domain
->
dimension
;
cdof
=
0.0
;
overwrite
=
0
;
format_user
=
NULL
;
format
=
(
char
*
)
" %g"
;
char
*
title1
=
NULL
;
char
*
title2
=
NULL
;
char
*
title3
=
NULL
;
while
(
iarg
<
nargnew
)
{
if
(
strcmp
(
arg
[
iarg
],
"norm"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
if
(
strcmp
(
arg
[
iarg
+
1
],
"all"
)
==
0
)
{
normflag
=
ALL
;
scaleflag
=
ATOM
;
}
else
if
(
strcmp
(
arg
[
iarg
+
1
],
"sample"
)
==
0
)
{
normflag
=
SAMPLE
;
scaleflag
=
ATOM
;
}
else
if
(
strcmp
(
arg
[
iarg
+
1
],
"none"
)
==
0
)
{
normflag
=
SAMPLE
;
scaleflag
=
NOSCALE
;
}
else
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"ave"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
if
(
strcmp
(
arg
[
iarg
+
1
],
"one"
)
==
0
)
ave
=
ONE
;
else
if
(
strcmp
(
arg
[
iarg
+
1
],
"running"
)
==
0
)
ave
=
RUNNING
;
else
if
(
strcmp
(
arg
[
iarg
+
1
],
"window"
)
==
0
)
ave
=
WINDOW
;
else
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
if
(
ave
==
WINDOW
)
{
if
(
iarg
+
3
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
nwindow
=
force
->
inumeric
(
FLERR
,
arg
[
iarg
+
2
]);
if
(
nwindow
<=
0
)
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
}
iarg
+=
2
;
if
(
ave
==
WINDOW
)
iarg
++
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"bias"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
biasflag
=
1
;
int
n
=
strlen
(
arg
[
iarg
+
1
])
+
1
;
id_bias
=
new
char
[
n
];
strcpy
(
id_bias
,
arg
[
iarg
+
1
]);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"adof"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
adof
=
force
->
numeric
(
FLERR
,
arg
[
iarg
+
1
]);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"cdof"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
cdof
=
force
->
numeric
(
FLERR
,
arg
[
iarg
+
1
]);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"file"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
if
(
me
==
0
)
{
fp
=
fopen
(
arg
[
iarg
+
1
],
"w"
);
if
(
fp
==
NULL
)
{
char
str
[
128
];
sprintf
(
str
,
"Cannot open fix ave/chunk file %s"
,
arg
[
iarg
+
1
]);
error
->
one
(
FLERR
,
str
);
}
}
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"overwrite"
)
==
0
)
{
overwrite
=
1
;
iarg
+=
1
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"format"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
delete
[]
format_user
;
int
n
=
strlen
(
arg
[
iarg
+
1
])
+
2
;
format_user
=
new
char
[
n
];
sprintf
(
format_user
,
" %s"
,
arg
[
iarg
+
1
]);
format
=
format_user
;
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"title1"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
delete
[]
title1
;
int
n
=
strlen
(
arg
[
iarg
+
1
])
+
1
;
title1
=
new
char
[
n
];
strcpy
(
title1
,
arg
[
iarg
+
1
]);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"title2"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
delete
[]
title2
;
int
n
=
strlen
(
arg
[
iarg
+
1
])
+
1
;
title2
=
new
char
[
n
];
strcpy
(
title2
,
arg
[
iarg
+
1
]);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"title3"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
delete
[]
title3
;
int
n
=
strlen
(
arg
[
iarg
+
1
])
+
1
;
title3
=
new
char
[
n
];
strcpy
(
title3
,
arg
[
iarg
+
1
]);
iarg
+=
2
;
}
else
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
}
// setup and error check
if
(
nevery
<=
0
||
nrepeat
<=
0
||
nfreq
<=
0
)
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
if
(
nfreq
%
nevery
||
nrepeat
*
nevery
>
nfreq
)
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
if
(
ave
!=
RUNNING
&&
overwrite
)
error
->
all
(
FLERR
,
"Illegal fix ave/chunk command"
);
if
(
biasflag
)
{
int
i
=
modify
->
find_compute
(
id_bias
);
if
(
i
<
0
)
error
->
all
(
FLERR
,
"Could not find compute ID for temperature bias"
);
tbias
=
modify
->
compute
[
i
];
if
(
tbias
->
tempflag
==
0
)
error
->
all
(
FLERR
,
"Bias compute does not calculate temperature"
);
if
(
tbias
->
tempbias
==
0
)
error
->
all
(
FLERR
,
"Bias compute does not calculate a velocity bias"
);
}
for
(
int
i
=
0
;
i
<
nvalues
;
i
++
)
{
if
(
which
[
i
]
==
COMPUTE
)
{
int
icompute
=
modify
->
find_compute
(
ids
[
i
]);
if
(
icompute
<
0
)
error
->
all
(
FLERR
,
"Compute ID for fix ave/chunk does not exist"
);
if
(
modify
->
compute
[
icompute
]
->
peratom_flag
==
0
)
error
->
all
(
FLERR
,
"Fix ave/chunk compute does not "
"calculate per-atom values"
);
if
(
argindex
[
i
]
==
0
&&
modify
->
compute
[
icompute
]
->
size_peratom_cols
!=
0
)
error
->
all
(
FLERR
,
"Fix ave/chunk compute does not "
"calculate a per-atom vector"
);
if
(
argindex
[
i
]
&&
modify
->
compute
[
icompute
]
->
size_peratom_cols
==
0
)
error
->
all
(
FLERR
,
"Fix ave/chunk compute does not "
"calculate a per-atom array"
);
if
(
argindex
[
i
]
&&
argindex
[
i
]
>
modify
->
compute
[
icompute
]
->
size_peratom_cols
)
error
->
all
(
FLERR
,
"Fix ave/chunk compute vector is accessed out-of-range"
);
}
else
if
(
which
[
i
]
==
FIX
)
{
int
ifix
=
modify
->
find_fix
(
ids
[
i
]);
if
(
ifix
<
0
)
error
->
all
(
FLERR
,
"Fix ID for fix ave/chunk does not exist"
);
if
(
modify
->
fix
[
ifix
]
->
peratom_flag
==
0
)
error
->
all
(
FLERR
,
"Fix ave/chunk fix does not calculate per-atom values"
);
if
(
argindex
[
i
]
==
0
&&
modify
->
fix
[
ifix
]
->
size_peratom_cols
!=
0
)
error
->
all
(
FLERR
,
"Fix ave/chunk fix does not calculate a per-atom vector"
);
if
(
argindex
[
i
]
&&
modify
->
fix
[
ifix
]
->
size_peratom_cols
==
0
)
error
->
all
(
FLERR
,
"Fix ave/chunk fix does not calculate a per-atom array"
);
if
(
argindex
[
i
]
&&
argindex
[
i
]
>
modify
->
fix
[
ifix
]
->
size_peratom_cols
)
error
->
all
(
FLERR
,
"Fix ave/chunk fix vector is accessed out-of-range"
);
}
else
if
(
which
[
i
]
==
VARIABLE
)
{
int
ivariable
=
input
->
variable
->
find
(
ids
[
i
]);
if
(
ivariable
<
0
)
error
->
all
(
FLERR
,
"Variable name for fix ave/chunk does not exist"
);
if
(
input
->
variable
->
atomstyle
(
ivariable
)
==
0
)
error
->
all
(
FLERR
,
"Fix ave/chunk variable is not atom-style variable"
);
}
}
// increment lock counter in compute chunk/atom
// only if nrepeat > 1 or ave = RUNNING/WINDOW,
// so that locking spans multiple timesteps
int
icompute
=
modify
->
find_compute
(
idchunk
);
if
(
icompute
<
0
)
error
->
all
(
FLERR
,
"Chunk/atom compute does not exist for fix ave/chunk"
);
cchunk
=
(
ComputeChunkAtom
*
)
modify
->
compute
[
icompute
];
if
(
strcmp
(
cchunk
->
style
,
"chunk/atom"
)
!=
0
)
error
->
all
(
FLERR
,
"Fix ave/chunk does not use chunk/atom compute"
);
if
(
nrepeat
>
1
||
ave
==
RUNNING
||
ave
==
WINDOW
)
cchunk
->
lockcount
++
;
lockforever
=
0
;
// print file comment lines
if
(
fp
&&
me
==
0
)
{
clearerr
(
fp
);
if
(
title1
)
fprintf
(
fp
,
"%s
\n
"
,
title1
);
else
fprintf
(
fp
,
"# Chunk-averaged data for fix %s and group %s
\n
"
,
id
,
arg
[
1
]);
if
(
title2
)
fprintf
(
fp
,
"%s
\n
"
,
title2
);
else
fprintf
(
fp
,
"# Timestep Number-of-chunks Total-count
\n
"
);
if
(
title3
)
fprintf
(
fp
,
"%s
\n
"
,
title3
);
else
{
int
compress
=
cchunk
->
compress
;
int
ncoord
=
cchunk
->
ncoord
;
if
(
!
compress
)
{
if
(
ncoord
==
0
)
fprintf
(
fp
,
"# Chunk Ncount"
);
else
if
(
ncoord
==
1
)
fprintf
(
fp
,
"# Chunk Coord1 Ncount"
);
else
if
(
ncoord
==
2
)
fprintf
(
fp
,
"# Chunk Coord1 Coord2 Ncount"
);
else
if
(
ncoord
==
3
)
fprintf
(
fp
,
"# Chunk Coord1 Coord2 Coord3 Ncount"
);
}
else
{
if
(
ncoord
==
0
)
fprintf
(
fp
,
"# Chunk OrigID Ncount"
);
else
if
(
ncoord
==
1
)
fprintf
(
fp
,
"# Chunk OrigID Coord1 Ncount"
);
else
if
(
ncoord
==
2
)
fprintf
(
fp
,
"# Chunk OrigID Coord1 Coord2 Ncount"
);
else
if
(
ncoord
==
3
)
fprintf
(
fp
,
"# Chunk OrigID Coord1 Coord2 Coord3 Ncount"
);
}
for
(
int
i
=
0
;
i
<
nvalues
;
i
++
)
fprintf
(
fp
,
" %s"
,
earg
[
i
]);
fprintf
(
fp
,
"
\n
"
);
}
if
(
ferror
(
fp
))
error
->
one
(
FLERR
,
"Error writing file header"
);
filepos
=
ftell
(
fp
);
}
delete
[]
title1
;
delete
[]
title2
;
delete
[]
title3
;
// if wildcard expansion occurred, free earg memory from expand_args()
// wait to do this until after file comment lines are printed
if
(
expand
)
{
for
(
int
i
=
0
;
i
<
nargnew
;
i
++
)
delete
[]
earg
[
i
];
memory
->
sfree
(
earg
);
}
// this fix produces a global array
// size_array_rows is variable and set by allocate()
int
compress
=
cchunk
->
compress
;
int
ncoord
=
cchunk
->
ncoord
;
colextra
=
compress
+
ncoord
;
array_flag
=
1
;
size_array_cols
=
colextra
+
1
+
nvalues
;
size_array_rows_variable
=
1
;
extarray
=
0
;
// initializations
irepeat
=
0
;
iwindow
=
window_limit
=
0
;
normcount
=
0
;
maxvar
=
0
;
varatom
=
NULL
;
count_one
=
count_many
=
count_sum
=
count_total
=
NULL
;
count_list
=
NULL
;
values_one
=
values_many
=
values_sum
=
values_total
=
NULL
;
values_list
=
NULL
;
maxchunk
=
0
;
nchunk
=
1
;
allocate
();
// nvalid = next step on which end_of_step does something
// add nvalid to all computes that store invocation times
// since don't know a priori which are invoked by this fix
// once in end_of_step() can set timestep for ones actually invoked
nvalid_last
=
-
1
;
nvalid
=
nextvalid
();
modify
->
addstep_compute_all
(
nvalid
);
}
/* ---------------------------------------------------------------------- */
FixAveChunk
::~
FixAveChunk
()
{
delete
[]
which
;
delete
[]
argindex
;
for
(
int
i
=
0
;
i
<
nvalues
;
i
++
)
delete
[]
ids
[
i
];
delete
[]
ids
;
delete
[]
value2index
;
if
(
fp
&&
me
==
0
)
fclose
(
fp
);
memory
->
destroy
(
varatom
);
memory
->
destroy
(
count_one
);
memory
->
destroy
(
count_many
);
memory
->
destroy
(
count_sum
);
memory
->
destroy
(
count_total
);
memory
->
destroy
(
count_list
);
memory
->
destroy
(
values_one
);
memory
->
destroy
(
values_many
);
memory
->
destroy
(
values_sum
);
memory
->
destroy
(
values_total
);
memory
->
destroy
(
values_list
);
// decrement lock counter in compute chunk/atom, it if still exists
if
(
nrepeat
>
1
||
ave
==
RUNNING
||
ave
==
WINDOW
)
{
int
icompute
=
modify
->
find_compute
(
idchunk
);
if
(
icompute
>=
0
)
{
cchunk
=
(
ComputeChunkAtom
*
)
modify
->
compute
[
icompute
];
if
(
ave
==
RUNNING
||
ave
==
WINDOW
)
cchunk
->
unlock
(
this
);
cchunk
->
lockcount
--
;
}
}
delete
[]
idchunk
;
which
=
NULL
;
argindex
=
NULL
;
ids
=
NULL
;
value2index
=
NULL
;
fp
=
NULL
;
varatom
=
NULL
;
count_one
=
NULL
;
count_many
=
NULL
;
count_sum
=
NULL
;
count_total
=
NULL
;
count_list
=
NULL
;
values_one
=
NULL
;
values_many
=
NULL
;
values_sum
=
NULL
;
values_total
=
NULL
;
values_list
=
NULL
;
idchunk
=
NULL
;
cchunk
=
NULL
;
}
/* ---------------------------------------------------------------------- */
int
FixAveChunk
::
setmask
()
{
int
mask
=
0
;
mask
|=
END_OF_STEP
;
return
mask
;
}
/* ---------------------------------------------------------------------- */
void
FixAveChunk
::
init
()
{
// set indices and check validity of all computes,fixes,variables
// check that fix frequency is acceptable
int
icompute
=
modify
->
find_compute
(
idchunk
);
if
(
icompute
<
0
)
error
->
all
(
FLERR
,
"Chunk/atom compute does not exist for fix ave/chunk"
);
cchunk
=
(
ComputeChunkAtom
*
)
modify
->
compute
[
icompute
];
if
(
biasflag
)
{
int
i
=
modify
->
find_compute
(
id_bias
);
if
(
i
<
0
)
error
->
all
(
FLERR
,
"Could not find compute ID for temperature bias"
);
tbias
=
modify
->
compute
[
i
];
}
for
(
int
m
=
0
;
m
<
nvalues
;
m
++
)
{
if
(
which
[
m
]
==
COMPUTE
)
{
int
icompute
=
modify
->
find_compute
(
ids
[
m
]);
if
(
icompute
<
0
)
error
->
all
(
FLERR
,
"Compute ID for fix ave/chunk does not exist"
);
value2index
[
m
]
=
icompute
;
}
else
if
(
which
[
m
]
==
FIX
)
{
int
ifix
=
modify
->
find_fix
(
ids
[
m
]);
if
(
ifix
<
0
)
error
->
all
(
FLERR
,
"Fix ID for fix ave/chunk does not exist"
);
value2index
[
m
]
=
ifix
;
if
(
nevery
%
modify
->
fix
[
ifix
]
->
peratom_freq
)
error
->
all
(
FLERR
,
"Fix for fix ave/chunk not computed at compatible time"
);
}
else
if
(
which
[
m
]
==
VARIABLE
)
{
int
ivariable
=
input
->
variable
->
find
(
ids
[
m
]);
if
(
ivariable
<
0
)
error
->
all
(
FLERR
,
"Variable name for fix ave/chunk does not exist"
);
value2index
[
m
]
=
ivariable
;
}
else
value2index
[
m
]
=
-
1
;
}
// need to reset nvalid if nvalid < ntimestep b/c minimize was performed
if
(
nvalid
<
update
->
ntimestep
)
{
irepeat
=
0
;
nvalid
=
nextvalid
();
modify
->
addstep_compute_all
(
nvalid
);
}
}
/* ----------------------------------------------------------------------
only does averaging if nvalid = current timestep
do not call setup_chunks(), even though fix ave/spatial called setup_bins()
b/c could cause nchunk to change if Nfreq epoch crosses 2 runs
does mean that if change_box is used between runs to change box size,
that nchunk may not track it
------------------------------------------------------------------------- */
void
FixAveChunk
::
setup
(
int
vflag
)
{
end_of_step
();
}
/* ---------------------------------------------------------------------- */
void
FixAveChunk
::
end_of_step
()
{
int
i
,
j
,
m
,
n
,
index
;
// skip if not step which requires doing something
// error check if timestep was reset in an invalid manner
bigint
ntimestep
=
update
->
ntimestep
;
if
(
ntimestep
<
nvalid_last
||
ntimestep
>
nvalid
)
error
->
all
(
FLERR
,
"Invalid timestep reset for fix ave/chunk"
);
if
(
ntimestep
!=
nvalid
)
return
;
nvalid_last
=
nvalid
;
// first sample within single Nfreq epoch
// zero out arrays that accumulate over many samples, but not across epochs
// invoke setup_chunks() to determine current nchunk
// re-allocate per-chunk arrays if needed
// invoke lock() in two cases:
// if nrepeat > 1: so nchunk cannot change until Nfreq epoch is over,
// will be unlocked on last repeat of this Nfreq
// if ave = RUNNING/WINDOW and not yet locked:
// set forever, will be unlocked in fix destructor
// wrap setup_chunks in clearstep/addstep b/c it may invoke computes
// both nevery and nfreq are future steps,
// since call below to cchunk->ichunk()
// does not re-invoke internal cchunk compute on this same step
if
(
irepeat
==
0
)
{
if
(
cchunk
->
computeflag
)
modify
->
clearstep_compute
();
nchunk
=
cchunk
->
setup_chunks
();
if
(
cchunk
->
computeflag
)
{
modify
->
addstep_compute
(
ntimestep
+
nevery
);
modify
->
addstep_compute
(
ntimestep
+
nfreq
);
}
allocate
();
if
(
nrepeat
>
1
&&
ave
==
ONE
)
cchunk
->
lock
(
this
,
ntimestep
,
ntimestep
+
(
nrepeat
-
1
)
*
nevery
);
else
if
((
ave
==
RUNNING
||
ave
==
WINDOW
)
&&
!
lockforever
)
{
cchunk
->
lock
(
this
,
update
->
ntimestep
,
-
1
);
lockforever
=
1
;
}
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
count_many
[
m
]
=
count_sum
[
m
]
=
0.0
;
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
values_many
[
m
][
i
]
=
0.0
;
}
}
// zero out arrays for one sample
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
count_one
[
m
]
=
0.0
;
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
values_one
[
m
][
i
]
=
0.0
;
}
// compute chunk/atom assigns atoms to chunk IDs
// extract ichunk index vector from compute
// ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms
// wrap compute_ichunk in clearstep/addstep b/c it may invoke computes
if
(
cchunk
->
computeflag
)
modify
->
clearstep_compute
();
cchunk
->
compute_ichunk
();
int
*
ichunk
=
cchunk
->
ichunk
;
if
(
cchunk
->
computeflag
)
modify
->
addstep_compute
(
ntimestep
+
nevery
);
// perform the computation for one sample
// count # of atoms in each bin
// accumulate results of attributes,computes,fixes,variables to local copy
// sum within each chunk, only include atoms in fix group
// compute/fix/variable may invoke computes so wrap with clear/add
int
*
mask
=
atom
->
mask
;
int
nlocal
=
atom
->
nlocal
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
ichunk
[
i
]
>
0
)
count_one
[
ichunk
[
i
]
-
1
]
++
;
modify
->
clearstep_compute
();
for
(
m
=
0
;
m
<
nvalues
;
m
++
)
{
n
=
value2index
[
m
];
j
=
argindex
[
m
];
// V,F adds velocities,forces to values
if
(
which
[
m
]
==
V
||
which
[
m
]
==
F
)
{
double
**
attribute
;
if
(
which
[
m
]
==
V
)
attribute
=
atom
->
v
;
else
attribute
=
atom
->
f
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
ichunk
[
i
]
>
0
)
{
index
=
ichunk
[
i
]
-
1
;
values_one
[
index
][
m
]
+=
attribute
[
i
][
j
];
}
// DENSITY_NUMBER adds 1 to values
}
else
if
(
which
[
m
]
==
DENSITY_NUMBER
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
ichunk
[
i
]
>
0
)
{
index
=
ichunk
[
i
]
-
1
;
values_one
[
index
][
m
]
+=
1.0
;
}
// DENSITY_MASS or MASS adds mass to values
}
else
if
(
which
[
m
]
==
DENSITY_MASS
||
which
[
m
]
==
MASS
)
{
int
*
type
=
atom
->
type
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
if
(
rmass
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
ichunk
[
i
]
>
0
)
{
index
=
ichunk
[
i
]
-
1
;
values_one
[
index
][
m
]
+=
rmass
[
i
];
}
}
else
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
ichunk
[
i
]
>
0
)
{
index
=
ichunk
[
i
]
-
1
;
values_one
[
index
][
m
]
+=
mass
[
type
[
i
]];
}
}
// TEMPERATURE adds KE to values
// subtract and restore velocity bias if requested
}
else
if
(
which
[
m
]
==
TEMPERATURE
)
{
if
(
biasflag
)
{
if
(
tbias
->
invoked_scalar
!=
ntimestep
)
tbias
->
compute_scalar
();
tbias
->
remove_bias_all
();
}
double
**
v
=
atom
->
v
;
int
*
type
=
atom
->
type
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
if
(
rmass
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
ichunk
[
i
]
>
0
)
{
index
=
ichunk
[
i
]
-
1
;
values_one
[
index
][
m
]
+=
(
v
[
i
][
0
]
*
v
[
i
][
0
]
+
v
[
i
][
1
]
*
v
[
i
][
1
]
+
v
[
i
][
2
]
*
v
[
i
][
2
])
*
rmass
[
i
];
}
}
else
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
ichunk
[
i
]
>
0
)
{
index
=
ichunk
[
i
]
-
1
;
values_one
[
index
][
m
]
+=
(
v
[
i
][
0
]
*
v
[
i
][
0
]
+
v
[
i
][
1
]
*
v
[
i
][
1
]
+
v
[
i
][
2
]
*
v
[
i
][
2
])
*
mass
[
type
[
i
]];
}
}
if
(
biasflag
)
tbias
->
restore_bias_all
();
// COMPUTE adds its scalar or vector component to values
// invoke compute if not previously invoked
}
else
if
(
which
[
m
]
==
COMPUTE
)
{
Compute
*
compute
=
modify
->
compute
[
n
];
if
(
!
(
compute
->
invoked_flag
&
INVOKED_PERATOM
))
{
compute
->
compute_peratom
();
compute
->
invoked_flag
|=
INVOKED_PERATOM
;
}
double
*
vector
=
compute
->
vector_atom
;
double
**
array
=
compute
->
array_atom
;
int
jm1
=
j
-
1
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
ichunk
[
i
]
>
0
)
{
index
=
ichunk
[
i
]
-
1
;
if
(
j
==
0
)
values_one
[
index
][
m
]
+=
vector
[
i
];
else
values_one
[
index
][
m
]
+=
array
[
i
][
jm1
];
}
// FIX adds its scalar or vector component to values
// access fix fields, guaranteed to be ready
}
else
if
(
which
[
m
]
==
FIX
)
{
double
*
vector
=
modify
->
fix
[
n
]
->
vector_atom
;
double
**
array
=
modify
->
fix
[
n
]
->
array_atom
;
int
jm1
=
j
-
1
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
ichunk
[
i
]
>
0
)
{
index
=
ichunk
[
i
]
-
1
;
if
(
j
==
0
)
values_one
[
index
][
m
]
+=
vector
[
i
];
else
values_one
[
index
][
m
]
+=
array
[
i
][
jm1
];
}
// VARIABLE adds its per-atom quantities to values
// evaluate atom-style variable
}
else
if
(
which
[
m
]
==
VARIABLE
)
{
if
(
atom
->
nmax
>
maxvar
)
{
maxvar
=
atom
->
nmax
;
memory
->
destroy
(
varatom
);
memory
->
create
(
varatom
,
maxvar
,
"ave/chunk:varatom"
);
}
input
->
variable
->
compute_atom
(
n
,
igroup
,
varatom
,
1
,
0
);
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
&&
ichunk
[
i
]
>
0
)
{
index
=
ichunk
[
i
]
-
1
;
values_one
[
index
][
m
]
+=
varatom
[
i
];
}
}
}
// process the current sample
// if normflag = ALL, accumulate values,count separately to many
// if normflag = SAMPLE, one = value/count, accumulate one to many
// count is MPI summed here, value is MPI summed below across samples
// exception is TEMPERATURE: normalize by DOF
// exception is DENSITYs: no normalize by atom count
// exception is scaleflag = NOSCALE : no normalize by atom count
// check last so other options can take precedence
double
mvv2e
=
force
->
mvv2e
;
double
boltz
=
force
->
boltz
;
if
(
normflag
==
ALL
)
{
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
count_many
[
m
]
+=
count_one
[
m
];
for
(
j
=
0
;
j
<
nvalues
;
j
++
)
values_many
[
m
][
j
]
+=
values_one
[
m
][
j
];
}
}
else
if
(
normflag
==
SAMPLE
)
{
MPI_Allreduce
(
count_one
,
count_many
,
nchunk
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
if
(
count_many
[
m
]
>
0.0
)
for
(
j
=
0
;
j
<
nvalues
;
j
++
)
{
if
(
which
[
j
]
==
TEMPERATURE
)
values_many
[
m
][
j
]
+=
mvv2e
*
values_one
[
m
][
j
]
/
((
cdof
+
adof
*
count_many
[
m
])
*
boltz
);
else
if
(
which
[
j
]
==
DENSITY_NUMBER
||
which
[
j
]
==
DENSITY_MASS
||
scaleflag
==
NOSCALE
)
values_many
[
m
][
j
]
+=
values_one
[
m
][
j
];
else
values_many
[
m
][
j
]
+=
values_one
[
m
][
j
]
/
count_many
[
m
];
}
count_sum
[
m
]
+=
count_many
[
m
];
}
}
// done if irepeat < nrepeat
// else reset irepeat and nvalid
irepeat
++
;
if
(
irepeat
<
nrepeat
)
{
nvalid
+=
nevery
;
modify
->
addstep_compute
(
nvalid
);
return
;
}
irepeat
=
0
;
nvalid
=
ntimestep
+
nfreq
-
(
nrepeat
-
1
)
*
nevery
;
modify
->
addstep_compute
(
nvalid
);
// unlock compute chunk/atom at end of Nfreq epoch
// do not unlock if ave = RUNNING or WINDOW
if
(
nrepeat
>
1
&&
ave
==
ONE
)
cchunk
->
unlock
(
this
);
// time average across samples
// if normflag = ALL, final is total value / total count
// exception is TEMPERATURE: normalize by DOF for total count
// exception is DENSITYs: normalize by repeat, not total count
// exception is scaleflag == NOSCALE: normalize by repeat, not total count
// check last so other options can take precedence
// if normflag = SAMPLE, final is sum of ave / repeat
double
repeat
=
nrepeat
;
double
mv2d
=
force
->
mv2d
;
if
(
normflag
==
ALL
)
{
MPI_Allreduce
(
count_many
,
count_sum
,
nchunk
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
MPI_Allreduce
(
&
values_many
[
0
][
0
],
&
values_sum
[
0
][
0
],
nchunk
*
nvalues
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
if
(
count_sum
[
m
]
>
0.0
)
for
(
j
=
0
;
j
<
nvalues
;
j
++
)
{
if
(
which
[
j
]
==
TEMPERATURE
)
values_sum
[
m
][
j
]
*=
mvv2e
/
((
cdof
+
adof
*
count_sum
[
m
])
*
boltz
);
else
if
(
which
[
j
]
==
DENSITY_MASS
)
values_sum
[
m
][
j
]
*=
mv2d
/
repeat
;
else
if
(
which
[
j
]
==
DENSITY_NUMBER
||
scaleflag
==
NOSCALE
)
values_sum
[
m
][
j
]
/=
repeat
;
else
values_sum
[
m
][
j
]
/=
count_sum
[
m
];
}
count_sum
[
m
]
/=
repeat
;
}
}
else
if
(
normflag
==
SAMPLE
)
{
MPI_Allreduce
(
&
values_many
[
0
][
0
],
&
values_sum
[
0
][
0
],
nchunk
*
nvalues
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
for
(
j
=
0
;
j
<
nvalues
;
j
++
)
values_sum
[
m
][
j
]
/=
repeat
;
count_sum
[
m
]
/=
repeat
;
}
}
// DENSITYs are additionally normalized by chunk volume
// use scalar or vector values for volume(s)
// if chunks are not spatial bins, chunk_volume_scalar = 1.0
for
(
j
=
0
;
j
<
nvalues
;
j
++
)
if
(
which
[
j
]
==
DENSITY_NUMBER
||
which
[
j
]
==
DENSITY_MASS
)
{
if
(
cchunk
->
chunk_volume_vec
)
{
double
*
chunk_volume_vec
=
cchunk
->
chunk_volume_vec
;
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
values_sum
[
m
][
j
]
/=
chunk_volume_vec
[
m
];
}
else
{
double
chunk_volume_scalar
=
cchunk
->
chunk_volume_scalar
;
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
values_sum
[
m
][
j
]
/=
chunk_volume_scalar
;
}
}
// if ave = ONE, only single Nfreq timestep value is needed
// if ave = RUNNING, combine with all previous Nfreq timestep values
// if ave = WINDOW, comine with nwindow most recent Nfreq timestep values
if
(
ave
==
ONE
)
{
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
values_total
[
m
][
i
]
=
values_sum
[
m
][
i
];
count_total
[
m
]
=
count_sum
[
m
];
}
normcount
=
1
;
}
else
if
(
ave
==
RUNNING
)
{
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
values_total
[
m
][
i
]
+=
values_sum
[
m
][
i
];
count_total
[
m
]
+=
count_sum
[
m
];
}
normcount
++
;
}
else
if
(
ave
==
WINDOW
)
{
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
{
values_total
[
m
][
i
]
+=
values_sum
[
m
][
i
];
if
(
window_limit
)
values_total
[
m
][
i
]
-=
values_list
[
iwindow
][
m
][
i
];
values_list
[
iwindow
][
m
][
i
]
=
values_sum
[
m
][
i
];
}
count_total
[
m
]
+=
count_sum
[
m
];
if
(
window_limit
)
count_total
[
m
]
-=
count_list
[
iwindow
][
m
];
count_list
[
iwindow
][
m
]
=
count_sum
[
m
];
}
iwindow
++
;
if
(
iwindow
==
nwindow
)
{
iwindow
=
0
;
window_limit
=
1
;
}
if
(
window_limit
)
normcount
=
nwindow
;
else
normcount
=
iwindow
;
}
// output result to file
if
(
fp
&&
me
==
0
)
{
clearerr
(
fp
);
if
(
overwrite
)
fseek
(
fp
,
filepos
,
SEEK_SET
);
double
count
=
0.0
;
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
count
+=
count_total
[
m
];
fprintf
(
fp
,
BIGINT_FORMAT
" %d %g
\n
"
,
ntimestep
,
nchunk
,
count
);
int
compress
=
cchunk
->
compress
;
int
*
chunkID
=
cchunk
->
chunkID
;
int
ncoord
=
cchunk
->
ncoord
;
double
**
coord
=
cchunk
->
coord
;
if
(
!
compress
)
{
if
(
ncoord
==
0
)
{
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
fprintf
(
fp
,
" %d %g"
,
m
+
1
,
count_total
[
m
]
/
normcount
);
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
fprintf
(
fp
,
format
,
values_total
[
m
][
i
]
/
normcount
);
fprintf
(
fp
,
"
\n
"
);
}
}
else
if
(
ncoord
==
1
)
{
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
fprintf
(
fp
,
" %d %g %g"
,
m
+
1
,
coord
[
m
][
0
],
count_total
[
m
]
/
normcount
);
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
fprintf
(
fp
,
format
,
values_total
[
m
][
i
]
/
normcount
);
fprintf
(
fp
,
"
\n
"
);
}
}
else
if
(
ncoord
==
2
)
{
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
fprintf
(
fp
,
" %d %g %g %g"
,
m
+
1
,
coord
[
m
][
0
],
coord
[
m
][
1
],
count_total
[
m
]
/
normcount
);
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
fprintf
(
fp
,
format
,
values_total
[
m
][
i
]
/
normcount
);
fprintf
(
fp
,
"
\n
"
);
}
}
else
if
(
ncoord
==
3
)
{
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
fprintf
(
fp
,
" %d %g %g %g %g"
,
m
+
1
,
coord
[
m
][
0
],
coord
[
m
][
1
],
coord
[
m
][
2
],
count_total
[
m
]
/
normcount
);
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
fprintf
(
fp
,
format
,
values_total
[
m
][
i
]
/
normcount
);
fprintf
(
fp
,
"
\n
"
);
}
}
}
else
{
int
j
;
if
(
ncoord
==
0
)
{
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
fprintf
(
fp
,
" %d %d %g"
,
m
+
1
,
chunkID
[
m
],
count_total
[
m
]
/
normcount
);
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
fprintf
(
fp
,
format
,
values_total
[
m
][
i
]
/
normcount
);
fprintf
(
fp
,
"
\n
"
);
}
}
else
if
(
ncoord
==
1
)
{
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
j
=
chunkID
[
m
];
fprintf
(
fp
,
" %d %d %g %g"
,
m
+
1
,
j
,
coord
[
j
-
1
][
0
],
count_total
[
m
]
/
normcount
);
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
fprintf
(
fp
,
format
,
values_total
[
m
][
i
]
/
normcount
);
fprintf
(
fp
,
"
\n
"
);
}
}
else
if
(
ncoord
==
2
)
{
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
j
=
chunkID
[
m
];
fprintf
(
fp
,
" %d %d %g %g %g"
,
m
+
1
,
j
,
coord
[
j
-
1
][
0
],
coord
[
j
-
1
][
1
],
count_total
[
m
]
/
normcount
);
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
fprintf
(
fp
,
format
,
values_total
[
m
][
i
]
/
normcount
);
fprintf
(
fp
,
"
\n
"
);
}
}
else
if
(
ncoord
==
3
)
{
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
j
=
chunkID
[
m
];
fprintf
(
fp
,
" %d %d %g %g %g %g"
,
m
+
1
,
j
,
coord
[
j
-
1
][
0
],
coord
[
j
-
1
][
1
],
coord
[
j
-
1
][
2
],
count_total
[
m
]
/
normcount
);
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
fprintf
(
fp
,
format
,
values_total
[
m
][
i
]
/
normcount
);
fprintf
(
fp
,
"
\n
"
);
}
}
}
if
(
ferror
(
fp
))
error
->
one
(
FLERR
,
"Error writing averaged chunk data"
);
fflush
(
fp
);
if
(
overwrite
)
{
long
fileend
=
ftell
(
fp
);
if
(
fileend
>
0
)
ftruncate
(
fileno
(
fp
),
fileend
);
}
}
}
/* ----------------------------------------------------------------------
allocate all per-chunk vectors
------------------------------------------------------------------------- */
void
FixAveChunk
::
allocate
()
{
size_array_rows
=
nchunk
;
// reallocate chunk arrays if needed
if
(
nchunk
>
maxchunk
)
{
maxchunk
=
nchunk
;
memory
->
grow
(
count_one
,
nchunk
,
"ave/chunk:count_one"
);
memory
->
grow
(
count_many
,
nchunk
,
"ave/chunk:count_many"
);
memory
->
grow
(
count_sum
,
nchunk
,
"ave/chunk:count_sum"
);
memory
->
grow
(
count_total
,
nchunk
,
"ave/chunk:count_total"
);
memory
->
grow
(
values_one
,
nchunk
,
nvalues
,
"ave/chunk:values_one"
);
memory
->
grow
(
values_many
,
nchunk
,
nvalues
,
"ave/chunk:values_many"
);
memory
->
grow
(
values_sum
,
nchunk
,
nvalues
,
"ave/chunk:values_sum"
);
memory
->
grow
(
values_total
,
nchunk
,
nvalues
,
"ave/chunk:values_total"
);
// only allocate count and values list for ave = WINDOW
if
(
ave
==
WINDOW
)
{
memory
->
create
(
count_list
,
nwindow
,
nchunk
,
"ave/chunk:count_list"
);
memory
->
create
(
values_list
,
nwindow
,
nchunk
,
nvalues
,
"ave/chunk:values_list"
);
}
// reinitialize regrown count/values total since they accumulate
int
i
,
m
;
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
values_total
[
m
][
i
]
=
0.0
;
count_total
[
m
]
=
0.0
;
}
}
}
/* ----------------------------------------------------------------------
return I,J array value
if I exceeds current nchunks, return 0.0 instead of generating an error
columns 1 to colextra = chunkID + ncoord
next column = count, remaining columns = Nvalues
------------------------------------------------------------------------- */
double
FixAveChunk
::
compute_array
(
int
i
,
int
j
)
{
if
(
values_total
==
NULL
)
return
0.0
;
if
(
i
>=
nchunk
)
return
0.0
;
if
(
j
<
colextra
)
{
if
(
cchunk
->
compress
)
{
if
(
j
==
0
)
return
(
double
)
cchunk
->
chunkID
[
i
];
return
cchunk
->
coord
[
i
][
j
-
1
];
}
else
return
cchunk
->
coord
[
i
][
j
];
}
j
-=
colextra
+
1
;
if
(
!
normcount
)
return
0.0
;
if
(
j
<
0
)
return
count_total
[
i
]
/
normcount
;
return
values_total
[
i
][
j
]
/
normcount
;
}
/* ----------------------------------------------------------------------
calculate nvalid = next step on which end_of_step does something
can be this timestep if multiple of nfreq and nrepeat = 1
else backup from next multiple of nfreq
------------------------------------------------------------------------- */
bigint
FixAveChunk
::
nextvalid
()
{
bigint
nvalid
=
(
update
->
ntimestep
/
nfreq
)
*
nfreq
+
nfreq
;
if
(
nvalid
-
nfreq
==
update
->
ntimestep
&&
nrepeat
==
1
)
nvalid
=
update
->
ntimestep
;
else
nvalid
-=
(
nrepeat
-
1
)
*
nevery
;
if
(
nvalid
<
update
->
ntimestep
)
nvalid
+=
nfreq
;
return
nvalid
;
}
/* ----------------------------------------------------------------------
memory usage of varatom and bins
------------------------------------------------------------------------- */
double
FixAveChunk
::
memory_usage
()
{
double
bytes
=
maxvar
*
sizeof
(
double
);
// varatom
bytes
+=
4
*
maxchunk
*
sizeof
(
double
);
// count one,many,sum,total
bytes
+=
nvalues
*
maxchunk
*
sizeof
(
double
);
// values one,many,sum,total
bytes
+=
nwindow
*
maxchunk
*
sizeof
(
double
);
// count_list
bytes
+=
nwindow
*
maxchunk
*
nvalues
*
sizeof
(
double
);
// values_list
return
bytes
;
}
Event Timeline
Log In to Comment