Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F66891352
fix_ave_correlate_long.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
Tue, Jun 18, 17:34
Size
24 KB
Mime Type
text/x-c
Expires
Thu, Jun 20, 17:34 (2 d)
Engine
blob
Format
Raw Data
Handle
18311946
Attached To
rLAMMPS lammps
fix_ave_correlate_long.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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Jorge Ramirez (jorge.ramirez@upm.es, Universidad Politecnica de Madrid),
Alexei Likhtman (University of Reading)
Structure and syntax of fix inspired by fix_ave_correlate
Scalar Correlator f(tau)=<A(t)A(t+tau)> and
Cross-correlator f(tau)=<A(t)B(t+tau)>
see J. Chem. Phys. 133, 154103 (2010)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include "fix_ave_correlate_long.h"
#include "update.h"
#include "modify.h"
#include "compute.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"
#include "force.h"
using
namespace
LAMMPS_NS
;
using
namespace
FixConst
;
enum
{
COMPUTE
,
FIX
,
VARIABLE
};
enum
{
AUTO
,
UPPER
,
LOWER
,
AUTOUPPER
,
AUTOLOWER
,
FULL
};
#define INVOKED_SCALAR 1
#define INVOKED_VECTOR 2
#define INVOKED_ARRAY 4
/* ---------------------------------------------------------------------- */
FixAveCorrelateLong
::
FixAveCorrelateLong
(
LAMMPS
*
lmp
,
int
narg
,
char
**
arg
)
:
Fix
(
lmp
,
narg
,
arg
)
{
// At least nevery nfrez and one value are needed
if
(
narg
<
6
)
error
->
all
(
FLERR
,
"Illegal fix ave/correlate/long command"
);
MPI_Comm_rank
(
world
,
&
me
);
nevery
=
force
->
inumeric
(
FLERR
,
arg
[
3
]);
nfreq
=
force
->
inumeric
(
FLERR
,
arg
[
4
]);
restart_global
=
1
;
global_freq
=
nfreq
;
// parse values until one isn't recognized
which
=
new
int
[
narg
-
5
];
argindex
=
new
int
[
narg
-
5
];
ids
=
new
char
*
[
narg
-
5
];
value2index
=
new
int
[
narg
-
5
];
nvalues
=
0
;
int
iarg
=
5
;
while
(
iarg
<
narg
)
{
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/correlate/long 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
);
delete
[]
suffix
;
nvalues
++
;
iarg
++
;
}
else
break
;
}
// optional args
type
=
AUTO
;
startstep
=
0
;
fp
=
NULL
;
overwrite
=
0
;
numcorrelators
=
20
;
p
=
16
;
m
=
2
;
char
*
title1
=
NULL
;
char
*
title2
=
NULL
;
while
(
iarg
<
narg
)
{
if
(
strcmp
(
arg
[
iarg
],
"type"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/correlate/long command"
);
if
(
strcmp
(
arg
[
iarg
+
1
],
"auto"
)
==
0
)
type
=
AUTO
;
else
if
(
strcmp
(
arg
[
iarg
+
1
],
"upper"
)
==
0
)
type
=
UPPER
;
else
if
(
strcmp
(
arg
[
iarg
+
1
],
"lower"
)
==
0
)
type
=
LOWER
;
else
if
(
strcmp
(
arg
[
iarg
+
1
],
"auto/upper"
)
==
0
)
type
=
AUTOUPPER
;
else
if
(
strcmp
(
arg
[
iarg
+
1
],
"auto/lower"
)
==
0
)
type
=
AUTOLOWER
;
else
if
(
strcmp
(
arg
[
iarg
+
1
],
"full"
)
==
0
)
type
=
FULL
;
else
error
->
all
(
FLERR
,
"Illegal fix ave/correlate/long command"
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"start"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/correlate/long command"
);
startstep
=
force
->
inumeric
(
FLERR
,
arg
[
iarg
+
1
]);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"ncorr"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/correlate/long command"
);
numcorrelators
=
force
->
inumeric
(
FLERR
,
arg
[
iarg
+
1
]);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"nlen"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/correlate/long command"
);
p
=
force
->
inumeric
(
FLERR
,
arg
[
iarg
+
1
]);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"ncount"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/correlate/long command"
);
m
=
force
->
inumeric
(
FLERR
,
arg
[
iarg
+
1
]);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"file"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/correlate/long command"
);
if
(
me
==
0
)
{
fp
=
fopen
(
arg
[
iarg
+
1
],
"w"
);
if
(
fp
==
NULL
)
{
char
str
[
128
];
sprintf
(
str
,
"Cannot open fix ave/correlate/long 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
],
"title1"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix ave/correlate/long 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/correlate/long command"
);
delete
[]
title2
;
int
n
=
strlen
(
arg
[
iarg
+
1
])
+
1
;
title2
=
new
char
[
n
];
strcpy
(
title2
,
arg
[
iarg
+
1
]);
iarg
+=
2
;
}
else
error
->
all
(
FLERR
,
"Illegal fix ave/correlate/long command"
);
}
if
(
p
%
m
!=
0
)
error
->
all
(
FLERR
,
"fix_correlator: p mod m must be 0"
);
dmin
=
p
/
m
;
length
=
numcorrelators
*
p
;
npcorr
=
0
;
kmax
=
0
;
// setup and error check
// for fix inputs, check that fix frequency is acceptable
if
(
nevery
<=
0
||
nfreq
<=
0
)
error
->
all
(
FLERR
,
"Illegal fix ave/correlate/long command"
);
if
(
nfreq
%
nevery
)
error
->
all
(
FLERR
,
"Illegal fix ave/correlate/long command"
);
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/correlate/long does not exist"
);
if
(
argindex
[
i
]
==
0
&&
modify
->
compute
[
icompute
]
->
scalar_flag
==
0
)
error
->
all
(
FLERR
,
"Fix ave/correlate/long compute does not calculate a scalar"
);
if
(
argindex
[
i
]
&&
modify
->
compute
[
icompute
]
->
vector_flag
==
0
)
error
->
all
(
FLERR
,
"Fix ave/correlate/long compute does not calculate a vector"
);
if
(
argindex
[
i
]
&&
argindex
[
i
]
>
modify
->
compute
[
icompute
]
->
size_vector
)
error
->
all
(
FLERR
,
"Fix ave/correlate/long 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/correlate/long does not exist"
);
if
(
argindex
[
i
]
==
0
&&
modify
->
fix
[
ifix
]
->
scalar_flag
==
0
)
error
->
all
(
FLERR
,
"Fix ave/correlate/long fix does not calculate a scalar"
);
if
(
argindex
[
i
]
&&
modify
->
fix
[
ifix
]
->
vector_flag
==
0
)
error
->
all
(
FLERR
,
"Fix ave/correlate/long fix does not calculate a vector"
);
if
(
argindex
[
i
]
&&
argindex
[
i
]
>
modify
->
fix
[
ifix
]
->
size_vector
)
error
->
all
(
FLERR
,
"Fix ave/correlate/long fix vector is accessed out-of-range"
);
if
(
nevery
%
modify
->
fix
[
ifix
]
->
global_freq
)
error
->
all
(
FLERR
,
"Fix for fix ave/correlate/long "
"not computed at compatible time"
);
}
else
if
(
which
[
i
]
==
VARIABLE
)
{
int
ivariable
=
input
->
variable
->
find
(
ids
[
i
]);
if
(
ivariable
<
0
)
error
->
all
(
FLERR
,
"Variable name for fix ave/correlate/long does not exist"
);
if
(
input
->
variable
->
equalstyle
(
ivariable
)
==
0
)
error
->
all
(
FLERR
,
"Fix ave/correlate/long variable is not equal-style variable"
);
}
}
// npair = # of correlation pairs to calculate
if
(
type
==
AUTO
)
npair
=
nvalues
;
if
(
type
==
UPPER
||
type
==
LOWER
)
npair
=
nvalues
*
(
nvalues
-
1
)
/
2
;
if
(
type
==
AUTOUPPER
||
type
==
AUTOLOWER
)
npair
=
nvalues
*
(
nvalues
+
1
)
/
2
;
if
(
type
==
FULL
)
npair
=
nvalues
*
nvalues
;
// print file comment lines
if
(
fp
&&
me
==
0
)
{
if
(
title1
)
fprintf
(
fp
,
"%s
\n
"
,
title1
);
else
fprintf
(
fp
,
"# Time-correlated data for fix %s
\n
"
,
id
);
if
(
title2
)
fprintf
(
fp
,
"%s
\n
"
,
title2
);
else
{
fprintf
(
fp
,
"# Time"
);
if
(
type
==
AUTO
)
for
(
int
i
=
0
;
i
<
nvalues
;
i
++
)
fprintf
(
fp
,
" %s*%s"
,
arg
[
5
+
i
],
arg
[
5
+
i
]);
else
if
(
type
==
UPPER
)
for
(
int
i
=
0
;
i
<
nvalues
;
i
++
)
for
(
int
j
=
i
+
1
;
j
<
nvalues
;
j
++
)
fprintf
(
fp
,
" %s*%s"
,
arg
[
5
+
i
],
arg
[
5
+
j
]);
else
if
(
type
==
LOWER
)
for
(
int
i
=
0
;
i
<
nvalues
;
i
++
)
for
(
int
j
=
0
;
j
<
i
-
1
;
j
++
)
fprintf
(
fp
,
" %s*%s"
,
arg
[
5
+
i
],
arg
[
5
+
j
]);
else
if
(
type
==
AUTOUPPER
)
for
(
int
i
=
0
;
i
<
nvalues
;
i
++
)
for
(
int
j
=
i
;
j
<
nvalues
;
j
++
)
fprintf
(
fp
,
" %s*%s"
,
arg
[
5
+
i
],
arg
[
5
+
j
]);
else
if
(
type
==
AUTOLOWER
)
for
(
int
i
=
0
;
i
<
nvalues
;
i
++
)
for
(
int
j
=
0
;
j
<
i
;
j
++
)
fprintf
(
fp
,
" %s*%s"
,
arg
[
5
+
i
],
arg
[
5
+
j
]);
else
if
(
type
==
FULL
)
for
(
int
i
=
0
;
i
<
nvalues
;
i
++
)
for
(
int
j
=
0
;
j
<
nvalues
;
j
++
)
fprintf
(
fp
,
" %s*%s"
,
arg
[
5
+
i
],
arg
[
5
+
j
]);
fprintf
(
fp
,
"
\n
"
);
}
filepos
=
ftell
(
fp
);
}
delete
[]
title1
;
delete
[]
title2
;
// allocate and initialize memory for calculated values and correlators
memory
->
create
(
values
,
nvalues
,
"correlator:values"
);
memory
->
create
(
shift
,
npair
,
numcorrelators
,
p
,
"correlator:shift"
);
memory
->
create
(
shift2
,
npair
,
numcorrelators
,
p
,
"correlator:shift2"
);
//NOT OPTMAL
memory
->
create
(
correlation
,
npair
,
numcorrelators
,
p
,
"correlator:correlation"
);
memory
->
create
(
accumulator
,
npair
,
numcorrelators
,
"correlator:accumulator"
);
memory
->
create
(
accumulator2
,
npair
,
numcorrelators
,
"correlator:accumulator2"
);
// NOT OPTIMAL
memory
->
create
(
ncorrelation
,
numcorrelators
,
p
,
"correlator:ncorrelation"
);
memory
->
create
(
naccumulator
,
numcorrelators
,
"correlator:naccumulator"
);
memory
->
create
(
insertindex
,
numcorrelators
,
"correlator:insertindex"
);
memory
->
create
(
t
,
length
,
"correlator:t"
);
memory
->
create
(
f
,
npair
,
length
,
"correlator:f"
);
for
(
int
i
=
0
;
i
<
npair
;
i
++
)
for
(
int
j
=
0
;
j
<
numcorrelators
;
j
++
)
{
for
(
int
k
=
0
;
k
<
p
;
k
++
)
{
shift
[
i
][
j
][
k
]
=-
2E10
;
shift2
[
i
][
j
][
k
]
=
0.0
;
correlation
[
i
][
j
][
k
]
=
0.0
;
}
accumulator
[
i
][
j
]
=
0.0
;
accumulator2
[
i
][
j
]
=
0.0
;
}
for
(
int
i
=
0
;
i
<
numcorrelators
;
i
++
)
{
for
(
int
j
=
0
;
j
<
p
;
j
++
)
ncorrelation
[
i
][
j
]
=
0
;
naccumulator
[
i
]
=
0
;
insertindex
[
i
]
=
0
;
}
for
(
int
i
=
0
;
i
<
length
;
i
++
)
t
[
i
]
=
0.0
;
for
(
int
i
=
0
;
i
<
npair
;
i
++
)
for
(
int
j
=
0
;
j
<
length
;
j
++
)
f
[
i
][
j
]
=
0.0
;
// 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
);
}
/* ---------------------------------------------------------------------- */
FixAveCorrelateLong
::~
FixAveCorrelateLong
()
{
delete
[]
which
;
delete
[]
argindex
;
delete
[]
value2index
;
for
(
int
i
=
0
;
i
<
nvalues
;
i
++
)
delete
[]
ids
[
i
];
delete
[]
ids
;
memory
->
destroy
(
values
);
memory
->
destroy
(
shift
);
memory
->
destroy
(
shift2
);
memory
->
destroy
(
correlation
);
memory
->
destroy
(
accumulator
);
memory
->
destroy
(
accumulator2
);
memory
->
destroy
(
ncorrelation
);
memory
->
destroy
(
naccumulator
);
memory
->
destroy
(
insertindex
);
memory
->
destroy
(
t
);
memory
->
destroy
(
f
);
if
(
fp
&&
me
==
0
)
fclose
(
fp
);
}
/* ---------------------------------------------------------------------- */
int
FixAveCorrelateLong
::
setmask
()
{
int
mask
=
0
;
mask
|=
END_OF_STEP
;
return
mask
;
}
/* ---------------------------------------------------------------------- */
void
FixAveCorrelateLong
::
init
()
{
// set current indices for all computes,fixes,variables
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/correlate/long does not exist"
);
value2index
[
i
]
=
icompute
;
}
else
if
(
which
[
i
]
==
FIX
)
{
int
ifix
=
modify
->
find_fix
(
ids
[
i
]);
if
(
ifix
<
0
)
error
->
all
(
FLERR
,
"Fix ID for fix ave/correlate/long does not exist"
);
value2index
[
i
]
=
ifix
;
}
else
if
(
which
[
i
]
==
VARIABLE
)
{
int
ivariable
=
input
->
variable
->
find
(
ids
[
i
]);
if
(
ivariable
<
0
)
error
->
all
(
FLERR
,
"Variable name for fix ave/correlate/long does not exist"
);
value2index
[
i
]
=
ivariable
;
}
}
// need to reset nvalid if nvalid < ntimestep b/c minimize was performed
if
(
nvalid
<
update
->
ntimestep
)
{
nvalid
=
nextvalid
();
modify
->
addstep_compute_all
(
nvalid
);
}
}
/* ----------------------------------------------------------------------
only does something if nvalid = current timestep
------------------------------------------------------------------------- */
void
FixAveCorrelateLong
::
setup
(
int
vflag
)
{
end_of_step
();
}
/* ---------------------------------------------------------------------- */
void
FixAveCorrelateLong
::
end_of_step
()
{
int
i
,
m
;
double
scalar
;
// 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/correlate/long"
);
if
(
ntimestep
!=
nvalid
)
return
;
nvalid_last
=
nvalid
;
// accumulate results of computes,fixes,variables to origin
// compute/fix/variable may invoke computes so wrap with clear/add
modify
->
clearstep_compute
();
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
{
m
=
value2index
[
i
];
scalar
=
0.0
;
// invoke compute if not previously invoked
if
(
which
[
i
]
==
COMPUTE
)
{
Compute
*
compute
=
modify
->
compute
[
m
];
if
(
argindex
[
i
]
==
0
)
{
if
(
!
(
compute
->
invoked_flag
&
INVOKED_SCALAR
))
{
compute
->
compute_scalar
();
compute
->
invoked_flag
|=
INVOKED_SCALAR
;
}
scalar
=
compute
->
scalar
;
}
else
{
if
(
!
(
compute
->
invoked_flag
&
INVOKED_VECTOR
))
{
compute
->
compute_vector
();
compute
->
invoked_flag
|=
INVOKED_VECTOR
;
}
scalar
=
compute
->
vector
[
argindex
[
i
]
-
1
];
}
// access fix fields, guaranteed to be ready
}
else
if
(
which
[
i
]
==
FIX
)
{
if
(
argindex
[
i
]
==
0
)
scalar
=
modify
->
fix
[
m
]
->
compute_scalar
();
else
scalar
=
modify
->
fix
[
m
]
->
compute_vector
(
argindex
[
i
]
-
1
);
// evaluate equal-style variable
}
else
if
(
which
[
i
]
==
VARIABLE
)
scalar
=
input
->
variable
->
compute_equal
(
m
);
values
[
i
]
=
scalar
;
}
// fistindex = index in values ring of earliest time sample
nvalid
+=
nevery
;
modify
->
addstep_compute
(
nvalid
);
// calculate all Cij() enabled by latest values
accumulate
();
if
(
ntimestep
%
nfreq
)
return
;
// output result to file
evaluate
();
if
(
fp
&&
me
==
0
)
{
if
(
overwrite
)
fseek
(
fp
,
filepos
,
SEEK_SET
);
fprintf
(
fp
,
"# Timestep: "
BIGINT_FORMAT
"
\n
"
,
ntimestep
);
for
(
unsigned
int
i
=
0
;
i
<
npcorr
;
++
i
)
{
fprintf
(
fp
,
"%lg "
,
t
[
i
]
*
update
->
dt
);
for
(
unsigned
int
j
=
0
;
j
<
npair
;
++
j
)
{
fprintf
(
fp
,
"%lg "
,
f
[
j
][
i
]);
}
fprintf
(
fp
,
"
\n
"
);
}
fflush
(
fp
);
if
(
overwrite
)
{
long
fileend
=
ftell
(
fp
);
if
(
fileend
>
0
)
ftruncate
(
fileno
(
fp
),
fileend
);
}
}
return
;
}
void
FixAveCorrelateLong
::
evaluate
()
{
unsigned
int
jm
=
0
;
// First correlator
for
(
unsigned
int
j
=
0
;
j
<
p
;
++
j
)
{
if
(
ncorrelation
[
0
][
j
]
>
0
)
{
t
[
jm
]
=
j
;
for
(
int
i
=
0
;
i
<
npair
;
++
i
)
f
[
i
][
jm
]
=
correlation
[
i
][
0
][
j
]
/
ncorrelation
[
0
][
j
];
++
jm
;
}
}
// Subsequent correlators
for
(
int
k
=
1
;
k
<
kmax
;
++
k
)
{
for
(
int
j
=
dmin
;
j
<
p
;
++
j
)
{
if
(
ncorrelation
[
k
][
j
]
>
0
)
{
t
[
jm
]
=
j
*
pow
((
double
)
m
,
k
);
for
(
int
i
=
0
;
i
<
npair
;
++
i
)
f
[
i
][
jm
]
=
correlation
[
i
][
k
][
j
]
/
ncorrelation
[
k
][
j
];
++
jm
;
}
}
}
npcorr
=
jm
;
}
/* ----------------------------------------------------------------------
accumulate correlation data using more recently added values
------------------------------------------------------------------------- */
void
FixAveCorrelateLong
::
accumulate
()
{
int
i
,
j
,
ipair
;
//printf("DEBUG %i %i\n", nvalues, npair);
if
(
type
==
AUTO
)
{
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
add
(
i
,
values
[
i
]);
}
else
if
(
type
==
UPPER
)
{
ipair
=
0
;
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
for
(
j
=
i
+
1
;
j
<
nvalues
;
j
++
)
add
(
ipair
++
,
values
[
i
],
values
[
j
]);
}
else
if
(
type
==
LOWER
)
{
ipair
=
0
;
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
for
(
j
=
0
;
j
<
i
;
j
++
)
add
(
ipair
++
,
values
[
i
],
values
[
j
]);
}
else
if
(
type
==
AUTOUPPER
)
{
ipair
=
0
;
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
for
(
j
=
i
;
j
<
nvalues
;
j
++
)
{
if
(
i
==
j
)
add
(
ipair
++
,
values
[
i
]);
else
add
(
ipair
++
,
values
[
i
],
values
[
j
]);
}
}
else
if
(
type
==
AUTOLOWER
)
{
ipair
=
0
;
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
for
(
j
=
0
;
j
<=
i
;
j
++
)
{
if
(
i
==
j
)
add
(
ipair
++
,
values
[
i
]);
else
add
(
ipair
++
,
values
[
i
],
values
[
j
]);
}
}
else
if
(
type
==
FULL
)
{
ipair
=
0
;
for
(
i
=
0
;
i
<
nvalues
;
i
++
)
for
(
j
=
0
;
j
<
nvalues
;
j
++
)
{
if
(
i
==
j
)
add
(
ipair
++
,
values
[
i
]);
else
add
(
ipair
++
,
values
[
i
],
values
[
j
]);
}
}
}
/* ----------------------------------------------------------------------
Add a scalar value to the autocorrelator k of pair i
------------------------------------------------------------------------- */
void
FixAveCorrelateLong
::
add
(
const
int
i
,
const
double
w
,
const
unsigned
int
k
){
// If we exceed the correlator side, the value is discarded
if
(
k
==
numcorrelators
)
return
;
if
(
k
>
kmax
)
kmax
=
k
;
// Insert new value in shift array
shift
[
i
][
k
][
insertindex
[
k
]]
=
w
;
// Add to accumulator and, if needed, add to next correlator
accumulator
[
i
][
k
]
+=
w
;
if
(
i
==
0
)
++
naccumulator
[
k
];
if
(
naccumulator
[
k
]
==
m
)
{
add
(
i
,
accumulator
[
i
][
k
]
/
m
,
k
+
1
);
accumulator
[
i
][
k
]
=
0
;
if
(
i
==
npair
-
1
)
naccumulator
[
k
]
=
0
;
}
// Calculate correlation function
unsigned
int
ind1
=
insertindex
[
k
];
if
(
k
==
0
)
{
// First correlator is different
int
ind2
=
ind1
;
for
(
unsigned
int
j
=
0
;
j
<
p
;
++
j
)
{
if
(
shift
[
i
][
k
][
ind2
]
>
-
1e10
)
{
correlation
[
i
][
k
][
j
]
+=
shift
[
i
][
k
][
ind1
]
*
shift
[
i
][
k
][
ind2
];
if
(
i
==
0
)
++
ncorrelation
[
k
][
j
];
}
--
ind2
;
if
(
ind2
<
0
)
ind2
+=
p
;
}
}
else
{
int
ind2
=
ind1
-
dmin
;
for
(
unsigned
int
j
=
dmin
;
j
<
p
;
++
j
)
{
if
(
ind2
<
0
)
ind2
+=
p
;
if
(
shift
[
i
][
k
][
ind2
]
>
-
1e10
)
{
correlation
[
i
][
k
][
j
]
+=
shift
[
i
][
k
][
ind1
]
*
shift
[
i
][
k
][
ind2
];
if
(
i
==
0
)
++
ncorrelation
[
k
][
j
];
}
--
ind2
;
}
}
if
(
i
==
npair
-
1
)
{
++
insertindex
[
k
];
if
(
insertindex
[
k
]
==
p
)
insertindex
[
k
]
=
0
;
}
}
/* ----------------------------------------------------------------------
Add 2 scalar values to the cross-correlator k of pair i
------------------------------------------------------------------------- */
void
FixAveCorrelateLong
::
add
(
const
int
i
,
const
double
wA
,
const
double
wB
,
const
unsigned
int
k
)
{
if
(
k
==
numcorrelators
)
return
;
if
(
k
>
kmax
)
kmax
=
k
;
shift
[
i
][
k
][
insertindex
[
k
]]
=
wA
;
shift2
[
i
][
k
][
insertindex
[
k
]]
=
wB
;
accumulator
[
i
][
k
]
+=
wA
;
accumulator2
[
i
][
k
]
+=
wB
;
if
(
i
==
0
)
++
naccumulator
[
k
];
if
(
naccumulator
[
k
]
==
m
)
{
add
(
i
,
accumulator
[
i
][
k
]
/
m
,
accumulator2
[
i
][
k
]
/
m
,
k
+
1
);
accumulator
[
i
][
k
]
=
0
;
accumulator2
[
i
][
k
]
=
0
;
if
(
i
==
npair
-
1
)
naccumulator
[
k
]
=
0
;
}
unsigned
int
ind1
=
insertindex
[
k
];
if
(
k
==
0
)
{
int
ind2
=
ind1
;
for
(
unsigned
int
j
=
0
;
j
<
p
;
++
j
)
{
if
(
shift
[
i
][
k
][
ind2
]
>
-
1e10
)
{
correlation
[
i
][
k
][
j
]
+=
shift
[
i
][
k
][
ind1
]
*
shift2
[
i
][
k
][
ind2
];
if
(
i
==
0
)
++
ncorrelation
[
k
][
j
];
}
--
ind2
;
if
(
ind2
<
0
)
ind2
+=
p
;
}
}
else
{
int
ind2
=
ind1
-
dmin
;
for
(
unsigned
int
j
=
dmin
;
j
<
p
;
++
j
)
{
if
(
ind2
<
0
)
ind2
+=
p
;
if
(
shift
[
i
][
k
][
ind2
]
>
-
1e10
)
{
correlation
[
i
][
k
][
j
]
+=
shift
[
i
][
k
][
ind1
]
*
shift2
[
i
][
k
][
ind2
];
if
(
i
==
0
)
++
ncorrelation
[
k
][
j
];
}
--
ind2
;
}
}
if
(
i
==
npair
-
1
)
{
++
insertindex
[
k
];
if
(
insertindex
[
k
]
==
p
)
insertindex
[
k
]
=
0
;
}
}
/* ----------------------------------------------------------------------
nvalid = next step on which end_of_step does something
this step if multiple of nevery, else next multiple
startstep is lower bound
------------------------------------------------------------------------- */
bigint
FixAveCorrelateLong
::
nextvalid
()
{
bigint
nvalid
=
update
->
ntimestep
;
if
(
startstep
>
nvalid
)
nvalid
=
startstep
;
if
(
nvalid
%
nevery
)
nvalid
=
(
nvalid
/
nevery
)
*
nevery
+
nevery
;
return
nvalid
;
}
/* ----------------------------------------------------------------------
memory_usage
------------------------------------------------------------------------- */
double
FixAveCorrelateLong
::
memory_usage
()
{
// shift: npair x numcorrelators x p
// shift2: npair x numcorrelators x p
// correlation: npair x numcorrelators x p
// accumulator: npair x numcorrelators
// accumulator2: npair x numcorrelators
// ncorrelation: numcorrelators x p
// naccumulator: numcorrelators
// insertindex: numcorrelators
// t: numcorrelators x p
// f: npair x numcorrelators x p
double
bytes
=
(
4
*
npair
*
numcorrelators
*
p
+
2
*
npair
*
numcorrelators
+
numcorrelators
*
p
)
*
sizeof
(
double
)
+
numcorrelators
*
p
*
sizeof
(
unsigned
long
int
)
+
2
*
numcorrelators
*
sizeof
(
unsigned
int
);
return
bytes
;
}
/* ----------------------------------------------------------------------
Write Restart data to restart file
------------------------------------------------------------------------- */
// Save everything except t and f
void
FixAveCorrelateLong
::
write_restart
(
FILE
*
fp
)
{
if
(
me
==
0
)
{
int
nsize
=
3
*
npair
*
numcorrelators
*
p
+
2
*
npair
*
numcorrelators
+
numcorrelators
*
p
+
2
*
numcorrelators
+
6
;
int
n
=
0
;
double
*
list
;
memory
->
create
(
list
,
nsize
,
"correlator:list"
);
list
[
n
++
]
=
npair
;
list
[
n
++
]
=
numcorrelators
;
list
[
n
++
]
=
p
;
list
[
n
++
]
=
m
;
list
[
n
++
]
=
nvalid
;
list
[
n
++
]
=
nvalid_last
;
for
(
int
i
=
0
;
i
<
npair
;
i
++
)
for
(
int
j
=
0
;
j
<
numcorrelators
;
j
++
)
{
for
(
int
k
=
0
;
k
<
p
;
k
++
)
{
list
[
n
++
]
=
shift
[
i
][
j
][
k
];
list
[
n
++
]
=
shift2
[
i
][
j
][
k
];
list
[
n
++
]
=
correlation
[
i
][
j
][
k
];
}
list
[
n
++
]
=
accumulator
[
i
][
j
];
list
[
n
++
]
=
accumulator2
[
i
][
j
];
}
for
(
int
i
=
0
;
i
<
numcorrelators
;
i
++
)
{
for
(
int
j
=
0
;
j
<
p
;
j
++
)
list
[
n
++
]
=
ncorrelation
[
i
][
j
];
list
[
n
++
]
=
naccumulator
[
i
];
list
[
n
++
]
=
insertindex
[
i
];
}
int
size
=
n
*
sizeof
(
double
);
fwrite
(
&
size
,
sizeof
(
int
),
1
,
fp
);
fwrite
(
list
,
sizeof
(
double
),
n
,
fp
);
memory
->
destroy
(
list
);
}
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void
FixAveCorrelateLong
::
restart
(
char
*
buf
)
{
int
n
=
0
;
double
*
list
=
(
double
*
)
buf
;
int
npairin
=
static_cast
<
int
>
(
list
[
n
++
]);
int
numcorrelatorsin
=
static_cast
<
int
>
(
list
[
n
++
]);
int
pin
=
static_cast
<
int
>
(
list
[
n
++
]);
int
min
=
static_cast
<
int
>
(
list
[
n
++
]);
nvalid
=
static_cast
<
int
>
(
list
[
n
++
]);
nvalid_last
=
static_cast
<
int
>
(
list
[
n
++
]);
if
((
npairin
!=
npair
)
||
(
numcorrelatorsin
!=
numcorrelators
)
||
(
pin
!=
p
)
||
(
min
!=
m
))
error
->
all
(
FLERR
,
"Fix ave/correlate/long: restart and input data are different"
);
for
(
int
i
=
0
;
i
<
npair
;
i
++
)
for
(
int
j
=
0
;
j
<
numcorrelators
;
j
++
)
{
for
(
int
k
=
0
;
k
<
p
;
k
++
)
{
shift
[
i
][
j
][
k
]
=
list
[
n
++
];
shift2
[
i
][
j
][
k
]
=
list
[
n
++
];
correlation
[
i
][
j
][
k
]
=
list
[
n
++
];
}
accumulator
[
i
][
j
]
=
list
[
n
++
];
accumulator2
[
i
][
j
]
=
list
[
n
++
];
}
for
(
int
i
=
0
;
i
<
numcorrelators
;
i
++
)
{
for
(
int
j
=
0
;
j
<
p
;
j
++
)
ncorrelation
[
i
][
j
]
=
static_cast
<
unsigned
long
int
>
(
list
[
n
++
]);
naccumulator
[
i
]
=
static_cast
<
unsigned
int
>
(
list
[
n
++
]);
insertindex
[
i
]
=
static_cast
<
unsigned
int
>
(
list
[
n
++
]);
}
}
Event Timeline
Log In to Comment