Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88007234
colvar.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
Wed, Oct 16, 06:59
Size
60 KB
Mime Type
text/x-c
Expires
Fri, Oct 18, 06:59 (2 d)
Engine
blob
Format
Raw Data
Handle
21699438
Attached To
rLAMMPS lammps
colvar.cpp
View Options
// -*- c++ -*-
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
#include "colvar.h"
#include "colvarcomp.h"
#include "colvarscript.h"
#include <algorithm>
/// Compare two cvcs using their names
/// Used to sort CVC array in scripted coordinates
bool
compare
(
colvar
::
cvc
*
i
,
colvar
::
cvc
*
j
)
{
return
i
->
name
<
j
->
name
;
}
colvar
::
colvar
(
std
::
string
const
&
conf
)
:
colvarparse
(
conf
)
{
cvm
::
log
(
"Initializing a new collective variable.
\n
"
);
int
error_code
=
COLVARS_OK
;
get_keyval
(
conf
,
"name"
,
this
->
name
,
(
std
::
string
(
"colvar"
)
+
cvm
::
to_str
(
cvm
::
colvars
.
size
()
+
1
)));
if
(
cvm
::
colvar_by_name
(
this
->
name
)
!=
NULL
)
{
cvm
::
error
(
"Error: this colvar cannot have the same name,
\"
"
+
this
->
name
+
"
\"
, as another colvar.
\n
"
,
INPUT_ERROR
);
return
;
}
// Initialize dependency members
// Could be a function defined in a different source file, for space?
this
->
description
=
"colvar "
+
this
->
name
;
// Initialize static array once and for all
init_cv_requires
();
kinetic_energy
=
0.0
;
potential_energy
=
0.0
;
cvm
::
combine_errors
(
error_code
,
init_components
(
conf
));
size_t
i
;
// Setup colvar as scripted function of components
if
(
get_keyval
(
conf
,
"scriptedFunction"
,
scripted_function
,
""
,
colvarparse
::
parse_silent
))
{
// Make feature available only on user request
provide
(
f_cv_scripted
);
enable
(
f_cv_scripted
);
cvm
::
log
(
"This colvar uses scripted function
\"
"
+
scripted_function
+
"
\"
."
);
std
::
string
type_str
;
get_keyval
(
conf
,
"scriptedFunctionType"
,
type_str
,
"scalar"
);
x
.
type
(
colvarvalue
::
type_notset
);
int
t
;
for
(
t
=
0
;
t
<
colvarvalue
::
type_all
;
t
++
)
{
if
(
type_str
==
colvarvalue
::
type_keyword
(
colvarvalue
::
Type
(
t
)))
{
x
.
type
(
colvarvalue
::
Type
(
t
));
break
;
}
}
if
(
x
.
type
()
==
colvarvalue
::
type_notset
)
{
cvm
::
error
(
"Could not parse scripted colvar type."
);
return
;
}
cvm
::
log
(
std
::
string
(
"Expecting colvar value of type "
)
+
colvarvalue
::
type_desc
(
x
.
type
()));
if
(
x
.
type
()
==
colvarvalue
::
type_vector
)
{
int
size
;
if
(
!
get_keyval
(
conf
,
"scriptedFunctionVectorSize"
,
size
))
{
cvm
::
error
(
"Error: no size specified for vector scripted function."
);
return
;
}
x
.
vector1d_value
.
resize
(
size
);
}
x_reported
.
type
(
x
);
// Sort array of cvcs based on their names
// Note: default CVC names are in input order for same type of CVC
std
::
sort
(
cvcs
.
begin
(),
cvcs
.
end
(),
compare
);
if
(
cvcs
.
size
()
>
1
)
{
cvm
::
log
(
"Sorted list of components for this scripted colvar:"
);
for
(
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
cvm
::
log
(
cvm
::
to_str
(
i
+
1
)
+
" "
+
cvcs
[
i
]
->
name
);
}
}
// Build ordered list of component values that will be
// passed to the script
for
(
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
sorted_cvc_values
.
push_back
(
&
(
cvcs
[
i
]
->
value
()));
}
}
if
(
!
is_enabled
(
f_cv_scripted
))
{
colvarvalue
const
&
cvc_value
=
(
cvcs
[
0
])
->
value
();
if
(
cvm
::
debug
())
cvm
::
log
(
"This collective variable is a "
+
colvarvalue
::
type_desc
(
cvc_value
.
type
())
+
((
cvc_value
.
size
()
>
1
)
?
" with "
+
cvm
::
to_str
(
cvc_value
.
size
())
+
" individual components.
\n
"
:
".
\n
"
));
x
.
type
(
cvc_value
);
x_reported
.
type
(
cvc_value
);
}
// If using scripted biases, any colvar may receive bias forces
// and will need its gradient
if
(
cvm
::
scripted_forces
())
{
enable
(
f_cv_gradient
);
}
// check for linear combinations
{
bool
lin
=
!
is_enabled
(
f_cv_scripted
);
for
(
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
// FIXME this is a reverse dependency, ie. cv feature depends on cvc flag
// need to clarify this case
// if ((cvcs[i])->b_debug_gradients)
// enable(task_gradients);
if
((
cvcs
[
i
])
->
sup_np
!=
1
)
{
if
(
cvm
::
debug
()
&&
lin
)
cvm
::
log
(
"Warning: You are using a non-linear polynomial "
"combination to define this collective variable, "
"some biasing methods may be unavailable.
\n
"
);
lin
=
false
;
if
((
cvcs
[
i
])
->
sup_np
<
0
)
{
cvm
::
log
(
"Warning: you chose a negative exponent in the combination; "
"if you apply forces, the simulation may become unstable "
"when the component
\"
"
+
(
cvcs
[
i
])
->
function_type
+
"
\"
approaches zero.
\n
"
);
}
}
}
feature_states
[
f_cv_linear
]
->
enabled
=
lin
;
}
// Colvar is homogeneous iff:
// - it is linear (hence not scripted)
// - all cvcs have coefficient 1 or -1
// i.e. sum or difference of cvcs
{
bool
homogeneous
=
is_enabled
(
f_cv_linear
);
for
(
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
if
((
std
::
fabs
(
cvcs
[
i
]
->
sup_coeff
)
-
1.0
)
>
1.0e-10
)
{
homogeneous
=
false
;
}
}
feature_states
[
f_cv_homogeneous
]
->
enabled
=
homogeneous
;
}
// Colvar is deemed periodic iff:
// - it is homogeneous
// - all cvcs are periodic
// - all cvcs have the same period
b_periodic
=
cvcs
[
0
]
->
b_periodic
&&
is_enabled
(
f_cv_homogeneous
);
period
=
cvcs
[
0
]
->
period
;
for
(
i
=
1
;
i
<
cvcs
.
size
();
i
++
)
{
if
(
!
cvcs
[
i
]
->
b_periodic
||
cvcs
[
i
]
->
period
!=
period
)
{
b_periodic
=
false
;
period
=
0.0
;
}
}
feature_states
[
f_cv_periodic
]
->
enabled
=
b_periodic
;
// check that cvcs are compatible
for
(
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
if
((
cvcs
[
i
])
->
b_periodic
&&
!
b_periodic
)
{
cvm
::
log
(
"Warning: although this component is periodic, the colvar will "
"not be treated as periodic, either because the exponent is not "
"1, or because multiple components are present. Make sure that "
"you know what you are doing!"
);
}
// components may have different types only for scripted functions
if
(
!
is_enabled
(
f_cv_scripted
)
&&
(
colvarvalue
::
check_types
(
cvcs
[
i
]
->
value
(),
cvcs
[
0
]
->
value
()))
)
{
cvm
::
error
(
"ERROR: you are definining this collective variable "
"by using components of different types. "
"You must use the same type in order to "
" sum them together.
\n
"
,
INPUT_ERROR
);
return
;
}
}
active_cvc_square_norm
=
0.
;
for
(
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
active_cvc_square_norm
+=
cvcs
[
i
]
->
sup_coeff
*
cvcs
[
i
]
->
sup_coeff
;
}
// at this point, the colvar's type is defined
f
.
type
(
value
());
f_accumulated
.
type
(
value
());
fb
.
type
(
value
());
get_keyval
(
conf
,
"width"
,
width
,
1.0
);
if
(
width
<=
0.0
)
{
cvm
::
error
(
"Error:
\"
width
\"
must be positive.
\n
"
,
INPUT_ERROR
);
}
// NOTE: not porting wall stuff to new deps, as this will change to a separate bias
// the grid functions will wait a little as well
lower_boundary
.
type
(
value
());
lower_wall
.
type
(
value
());
upper_boundary
.
type
(
value
());
upper_wall
.
type
(
value
());
feature_states
[
f_cv_scalar
]
->
enabled
=
(
value
().
type
()
==
colvarvalue
::
type_scalar
);
if
(
is_enabled
(
f_cv_scalar
))
{
if
(
get_keyval
(
conf
,
"lowerBoundary"
,
lower_boundary
,
lower_boundary
))
{
provide
(
f_cv_lower_boundary
);
enable
(
f_cv_lower_boundary
);
}
get_keyval
(
conf
,
"lowerWallConstant"
,
lower_wall_k
,
0.0
);
if
(
lower_wall_k
>
0.0
)
{
get_keyval
(
conf
,
"lowerWall"
,
lower_wall
,
lower_boundary
);
enable
(
f_cv_lower_wall
);
}
if
(
get_keyval
(
conf
,
"upperBoundary"
,
upper_boundary
,
upper_boundary
))
{
provide
(
f_cv_upper_boundary
);
enable
(
f_cv_upper_boundary
);
}
get_keyval
(
conf
,
"upperWallConstant"
,
upper_wall_k
,
0.0
);
if
(
upper_wall_k
>
0.0
)
{
get_keyval
(
conf
,
"upperWall"
,
upper_wall
,
upper_boundary
);
enable
(
f_cv_upper_wall
);
}
}
if
(
is_enabled
(
f_cv_lower_boundary
))
{
get_keyval
(
conf
,
"hardLowerBoundary"
,
hard_lower_boundary
,
false
);
}
if
(
is_enabled
(
f_cv_upper_boundary
))
{
get_keyval
(
conf
,
"hardUpperBoundary"
,
hard_upper_boundary
,
false
);
}
// consistency checks for boundaries and walls
if
(
is_enabled
(
f_cv_lower_boundary
)
&&
is_enabled
(
f_cv_upper_boundary
))
{
if
(
lower_boundary
>=
upper_boundary
)
{
cvm
::
error
(
"Error: the upper boundary, "
+
cvm
::
to_str
(
upper_boundary
)
+
", is not higher than the lower boundary, "
+
cvm
::
to_str
(
lower_boundary
)
+
".
\n
"
,
INPUT_ERROR
);
}
}
if
(
is_enabled
(
f_cv_lower_wall
)
&&
is_enabled
(
f_cv_upper_wall
))
{
if
(
lower_wall
>=
upper_wall
)
{
cvm
::
error
(
"Error: the upper wall, "
+
cvm
::
to_str
(
upper_wall
)
+
", is not higher than the lower wall, "
+
cvm
::
to_str
(
lower_wall
)
+
".
\n
"
,
INPUT_ERROR
);
}
}
get_keyval
(
conf
,
"expandBoundaries"
,
expand_boundaries
,
false
);
if
(
expand_boundaries
&&
periodic_boundaries
())
{
cvm
::
error
(
"Error: trying to expand boundaries that already "
"cover a whole period of a periodic colvar.
\n
"
,
INPUT_ERROR
);
}
if
(
expand_boundaries
&&
hard_lower_boundary
&&
hard_upper_boundary
)
{
cvm
::
error
(
"Error: inconsistent configuration "
"(trying to expand boundaries with both "
"hardLowerBoundary and hardUpperBoundary enabled).
\n
"
,
INPUT_ERROR
);
}
get_keyval
(
conf
,
"timeStepFactor"
,
time_step_factor
,
1
);
{
bool
b_extended_Lagrangian
;
get_keyval
(
conf
,
"extendedLagrangian"
,
b_extended_Lagrangian
,
false
);
if
(
b_extended_Lagrangian
)
{
cvm
::
real
temp
,
tolerance
,
period
;
cvm
::
log
(
"Enabling the extended Lagrangian term for colvar
\"
"
+
this
->
name
+
"
\"
.
\n
"
);
// Make feature available only on user request
provide
(
f_cv_extended_Lagrangian
);
enable
(
f_cv_extended_Lagrangian
);
provide
(
f_cv_Langevin
);
xr
.
type
(
value
());
vr
.
type
(
value
());
fr
.
type
(
value
());
const
bool
found
=
get_keyval
(
conf
,
"extendedTemp"
,
temp
,
cvm
::
temperature
());
if
(
temp
<=
0.0
)
{
if
(
found
)
cvm
::
error
(
"Error:
\"
extendedTemp
\"
must be positive.
\n
"
,
INPUT_ERROR
);
else
cvm
::
error
(
"Error: a positive temperature must be provided, either "
"by enabling a thermostat, or through
\"
extendedTemp
\"
.
\n
"
,
INPUT_ERROR
);
}
get_keyval
(
conf
,
"extendedFluctuation"
,
tolerance
);
if
(
tolerance
<=
0.0
)
{
cvm
::
error
(
"Error:
\"
extendedFluctuation
\"
must be positive.
\n
"
,
INPUT_ERROR
);
}
ext_force_k
=
cvm
::
boltzmann
()
*
temp
/
(
tolerance
*
tolerance
);
cvm
::
log
(
"Computed extended system force constant: "
+
cvm
::
to_str
(
ext_force_k
)
+
" kcal/mol/U^2"
);
get_keyval
(
conf
,
"extendedTimeConstant"
,
period
,
200.0
);
if
(
period
<=
0.0
)
{
cvm
::
error
(
"Error:
\"
extendedTimeConstant
\"
must be positive.
\n
"
,
INPUT_ERROR
);
}
ext_mass
=
(
cvm
::
boltzmann
()
*
temp
*
period
*
period
)
/
(
4.0
*
PI
*
PI
*
tolerance
*
tolerance
);
cvm
::
log
(
"Computed fictitious mass: "
+
cvm
::
to_str
(
ext_mass
)
+
" kcal/mol/(U/fs)^2 (U: colvar unit)"
);
{
bool
b_output_energy
;
get_keyval
(
conf
,
"outputEnergy"
,
b_output_energy
,
false
);
if
(
b_output_energy
)
{
enable
(
f_cv_output_energy
);
}
}
get_keyval
(
conf
,
"extendedLangevinDamping"
,
ext_gamma
,
1.0
);
if
(
ext_gamma
<
0.0
)
{
cvm
::
error
(
"Error:
\"
extendedLangevinDamping
\"
may not be negative.
\n
"
,
INPUT_ERROR
);
}
if
(
ext_gamma
!=
0.0
)
{
enable
(
f_cv_Langevin
);
ext_gamma
*=
1.0e-3
;
// convert from ps-1 to fs-1
ext_sigma
=
std
::
sqrt
(
2.0
*
cvm
::
boltzmann
()
*
temp
*
ext_gamma
*
ext_mass
/
cvm
::
dt
());
}
}
}
{
bool
b_output_value
;
get_keyval
(
conf
,
"outputValue"
,
b_output_value
,
true
);
if
(
b_output_value
)
{
enable
(
f_cv_output_value
);
}
}
{
bool
b_output_velocity
;
get_keyval
(
conf
,
"outputVelocity"
,
b_output_velocity
,
false
);
if
(
b_output_velocity
)
{
enable
(
f_cv_output_velocity
);
}
}
{
bool
b_output_system_force
;
get_keyval
(
conf
,
"outputSystemForce"
,
b_output_system_force
,
false
);
if
(
b_output_system_force
)
{
enable
(
f_cv_output_system_force
);
}
}
{
bool
b_output_applied_force
;
get_keyval
(
conf
,
"outputAppliedForce"
,
b_output_applied_force
,
false
);
if
(
b_output_applied_force
)
{
enable
(
f_cv_output_applied_force
);
}
}
// Start in active state by default
enable
(
f_cv_active
);
// Make sure dependency side-effects are correct
refresh_deps
();
x_old
.
type
(
value
());
v_fdiff
.
type
(
value
());
v_reported
.
type
(
value
());
fj
.
type
(
value
());
ft
.
type
(
value
());
ft_reported
.
type
(
value
());
if
(
cvm
::
b_analysis
)
parse_analysis
(
conf
);
if
(
cvm
::
debug
())
cvm
::
log
(
"Done initializing collective variable
\"
"
+
this
->
name
+
"
\"
.
\n
"
);
}
// read the configuration and set up corresponding instances, for
// each type of component implemented
template
<
typename
def_class_name
>
int
colvar
::
init_components_type
(
std
::
string
const
&
conf
,
char
const
*
def_desc
,
char
const
*
def_config_key
)
{
size_t
def_count
=
0
;
std
::
string
def_conf
=
""
;
size_t
pos
=
0
;
while
(
this
->
key_lookup
(
conf
,
def_config_key
,
def_conf
,
pos
)
)
{
if
(
!
def_conf
.
size
())
continue
;
cvm
::
log
(
"Initializing "
"a new
\"
"
+
std
::
string
(
def_config_key
)
+
"
\"
component"
+
(
cvm
::
debug
()
?
", with configuration:
\n
"
+
def_conf
:
".
\n
"
));
cvm
::
increase_depth
();
cvc
*
cvcp
=
new
def_class_name
(
def_conf
);
if
(
cvcp
!=
NULL
)
{
cvcs
.
push_back
(
cvcp
);
cvcp
->
check_keywords
(
def_conf
,
def_config_key
);
if
(
cvm
::
get_error
())
{
cvm
::
error
(
"Error: in setting up component
\"
"
+
std
::
string
(
def_config_key
)
+
"
\"
.
\n
"
,
INPUT_ERROR
);
return
INPUT_ERROR
;
}
cvm
::
decrease_depth
();
}
else
{
cvm
::
error
(
"Error: in allocating component
\"
"
+
std
::
string
(
def_config_key
)
+
"
\"
.
\n
"
,
MEMORY_ERROR
);
return
MEMORY_ERROR
;
}
if
(
(
cvcp
->
period
!=
0.0
)
||
(
cvcp
->
wrap_center
!=
0.0
)
)
{
if
(
(
cvcp
->
function_type
!=
std
::
string
(
"distance_z"
))
&&
(
cvcp
->
function_type
!=
std
::
string
(
"dihedral"
))
&&
(
cvcp
->
function_type
!=
std
::
string
(
"spin_angle"
))
)
{
cvm
::
error
(
"Error: invalid use of period and/or "
"wrapAround in a
\"
"
+
std
::
string
(
def_config_key
)
+
"
\"
component.
\n
"
+
"Period: "
+
cvm
::
to_str
(
cvcp
->
period
)
+
" wrapAround: "
+
cvm
::
to_str
(
cvcp
->
wrap_center
),
INPUT_ERROR
);
return
INPUT_ERROR
;
}
}
if
(
!
cvcs
.
back
()
->
name
.
size
())
{
std
::
ostringstream
s
;
s
<<
def_config_key
<<
std
::
setfill
(
'0'
)
<<
std
::
setw
(
4
)
<<
++
def_count
;
cvcs
.
back
()
->
name
=
s
.
str
();
/* pad cvc number for correct ordering when sorting by name */
}
cvcs
.
back
()
->
setup
();
if
(
cvm
::
debug
())
{
cvm
::
log
(
"Done initializing a
\"
"
+
std
::
string
(
def_config_key
)
+
"
\"
component"
+
(
cvm
::
debug
()
?
", named
\"
"
+
cvcs
.
back
()
->
name
+
"
\"
"
:
""
)
+
".
\n
"
);
}
def_conf
=
""
;
if
(
cvm
::
debug
())
{
cvm
::
log
(
"Parsed "
+
cvm
::
to_str
(
cvcs
.
size
())
+
" components at this time.
\n
"
);
}
}
return
COLVARS_OK
;
}
int
colvar
::
init_components
(
std
::
string
const
&
conf
)
{
int
error_code
=
COLVARS_OK
;
cvm
::
combine_errors
(
error_code
,
init_components_type
<
distance
>
(
conf
,
"distance"
,
"distance"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
distance_vec
>
(
conf
,
"distance vector"
,
"distanceVec"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
cartesian
>
(
conf
,
"Cartesian coordinates"
,
"cartesian"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
distance_dir
>
(
conf
,
"distance vector "
"direction"
,
"distanceDir"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
distance_z
>
(
conf
,
"distance projection "
"on an axis"
,
"distanceZ"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
distance_xy
>
(
conf
,
"distance projection "
"on a plane"
,
"distanceXY"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
distance_inv
>
(
conf
,
"average distance "
"weighted by inverse power"
,
"distanceInv"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
distance_pairs
>
(
conf
,
"N1xN2-long vector "
"of pairwise distances"
,
"distancePairs"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
coordnum
>
(
conf
,
"coordination "
"number"
,
"coordNum"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
selfcoordnum
>
(
conf
,
"self-coordination "
"number"
,
"selfCoordNum"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
angle
>
(
conf
,
"angle"
,
"angle"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
dipole_angle
>
(
conf
,
"dipole angle"
,
"dipoleAngle"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
dihedral
>
(
conf
,
"dihedral"
,
"dihedral"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
h_bond
>
(
conf
,
"hydrogen bond"
,
"hBond"
));
// cvm::combine_errors(error_code, init_components_type<alpha_dihedrals>(conf, "alpha helix", "alphaDihedrals"));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
alpha_angles
>
(
conf
,
"alpha helix"
,
"alpha"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
dihedPC
>
(
conf
,
"dihedral "
"principal component"
,
"dihedralPC"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
orientation
>
(
conf
,
"orientation"
,
"orientation"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
orientation_angle
>
(
conf
,
"orientation "
"angle"
,
"orientationAngle"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
orientation_proj
>
(
conf
,
"orientation "
"projection"
,
"orientationProj"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
tilt
>
(
conf
,
"tilt"
,
"tilt"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
spin_angle
>
(
conf
,
"spin angle"
,
"spinAngle"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
rmsd
>
(
conf
,
"RMSD"
,
"rmsd"
));
// cvm::combine_errors(error_code, init_components_type <logmsd>(conf,"logarithm of MSD", "logmsd"));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
gyration
>
(
conf
,
"radius of "
"gyration"
,
"gyration"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
inertia
>
(
conf
,
"moment of "
"inertia"
,
"inertia"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
inertia_z
>
(
conf
,
"moment of inertia around an axis"
,
"inertiaZ"
));
cvm
::
combine_errors
(
error_code
,
init_components_type
<
eigenvector
>
(
conf
,
"eigenvector"
,
"eigenvector"
));
if
(
!
cvcs
.
size
()
||
(
error_code
!=
COLVARS_OK
))
{
cvm
::
error
(
"Error: no valid components were provided "
"for this collective variable.
\n
"
,
INPUT_ERROR
);
return
INPUT_ERROR
;
}
n_active_cvcs
=
cvcs
.
size
();
cvm
::
log
(
"All components initialized.
\n
"
);
// Store list of children cvcs for dependency checking purposes
for
(
size_t
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
add_child
(
cvcs
[
i
]);
}
return
COLVARS_OK
;
}
int
colvar
::
refresh_deps
()
{
// If enabled features are changed upstream, the features below should be refreshed
if
(
is_enabled
(
f_cv_system_force_calc
))
{
cvm
::
request_system_force
();
}
if
(
is_enabled
(
f_cv_collect_gradient
)
&&
atom_ids
.
size
()
==
0
)
{
build_atom_list
();
}
return
COLVARS_OK
;
}
void
colvar
::
build_atom_list
(
void
)
{
// If atomic gradients are requested, build full list of atom ids from all cvcs
std
::
list
<
int
>
temp_id_list
;
for
(
size_t
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
for
(
size_t
j
=
0
;
j
<
cvcs
[
i
]
->
atom_groups
.
size
();
j
++
)
{
cvm
::
atom_group
&
ag
=
*
(
cvcs
[
i
]
->
atom_groups
[
j
]);
for
(
size_t
k
=
0
;
k
<
ag
.
size
();
k
++
)
{
temp_id_list
.
push_back
(
ag
[
k
].
id
);
}
}
}
temp_id_list
.
sort
();
temp_id_list
.
unique
();
// atom_ids = std::vector<int> (temp_id_list.begin(), temp_id_list.end());
unsigned
int
id_i
=
0
;
std
::
list
<
int
>::
iterator
li
;
for
(
li
=
temp_id_list
.
begin
();
li
!=
temp_id_list
.
end
();
++
li
)
{
atom_ids
[
id_i
]
=
*
li
;
id_i
++
;
}
temp_id_list
.
clear
();
atomic_gradients
.
resize
(
atom_ids
.
size
());
if
(
atom_ids
.
size
())
{
if
(
cvm
::
debug
())
cvm
::
log
(
"Colvar: created atom list with "
+
cvm
::
to_str
(
atom_ids
.
size
())
+
" atoms.
\n
"
);
}
else
{
cvm
::
log
(
"Warning: colvar components communicated no atom IDs.
\n
"
);
}
}
int
colvar
::
parse_analysis
(
std
::
string
const
&
conf
)
{
// if (cvm::debug())
// cvm::log ("Parsing analysis flags for collective variable \""+
// this->name+"\".\n");
runave_length
=
0
;
bool
b_runave
=
false
;
if
(
get_keyval
(
conf
,
"runAve"
,
b_runave
)
&&
b_runave
)
{
enable
(
f_cv_runave
);
get_keyval
(
conf
,
"runAveLength"
,
runave_length
,
1000
);
get_keyval
(
conf
,
"runAveStride"
,
runave_stride
,
1
);
if
((
cvm
::
restart_out_freq
%
runave_stride
)
!=
0
)
{
cvm
::
error
(
"Error: runAveStride must be commensurate with the restart frequency.
\n
"
,
INPUT_ERROR
);
}
std
::
string
runave_outfile
;
get_keyval
(
conf
,
"runAveOutputFile"
,
runave_outfile
,
std
::
string
(
cvm
::
output_prefix
+
"."
+
this
->
name
+
".runave.traj"
));
size_t
const
this_cv_width
=
x
.
output_width
(
cvm
::
cv_width
);
cvm
::
backup_file
(
runave_outfile
.
c_str
());
runave_os
.
open
(
runave_outfile
.
c_str
());
runave_os
<<
"# "
<<
cvm
::
wrap_string
(
"step"
,
cvm
::
it_width
-
2
)
<<
" "
<<
cvm
::
wrap_string
(
"running average"
,
this_cv_width
)
<<
" "
<<
cvm
::
wrap_string
(
"running stddev"
,
this_cv_width
)
<<
"
\n
"
;
}
acf_length
=
0
;
bool
b_acf
=
false
;
if
(
get_keyval
(
conf
,
"corrFunc"
,
b_acf
)
&&
b_acf
)
{
enable
(
f_cv_corrfunc
);
std
::
string
acf_colvar_name
;
get_keyval
(
conf
,
"corrFuncWithColvar"
,
acf_colvar_name
,
this
->
name
);
if
(
acf_colvar_name
==
this
->
name
)
{
cvm
::
log
(
"Calculating auto-correlation function.
\n
"
);
}
else
{
cvm
::
log
(
"Calculating correlation function with
\"
"
+
this
->
name
+
"
\"
.
\n
"
);
}
std
::
string
acf_type_str
;
get_keyval
(
conf
,
"corrFuncType"
,
acf_type_str
,
to_lower_cppstr
(
std
::
string
(
"velocity"
)));
if
(
acf_type_str
==
to_lower_cppstr
(
std
::
string
(
"coordinate"
)))
{
acf_type
=
acf_coor
;
}
else
if
(
acf_type_str
==
to_lower_cppstr
(
std
::
string
(
"velocity"
)))
{
acf_type
=
acf_vel
;
enable
(
f_cv_fdiff_velocity
);
if
(
acf_colvar_name
.
size
())
(
cvm
::
colvar_by_name
(
acf_colvar_name
))
->
enable
(
f_cv_fdiff_velocity
);
}
else
if
(
acf_type_str
==
to_lower_cppstr
(
std
::
string
(
"coordinate_p2"
)))
{
acf_type
=
acf_p2coor
;
}
else
{
cvm
::
log
(
"Unknown type of correlation function,
\"
"
+
acf_type_str
+
"
\"
.
\n
"
);
cvm
::
set_error_bit
(
INPUT_ERROR
);
}
get_keyval
(
conf
,
"corrFuncOffset"
,
acf_offset
,
0
);
get_keyval
(
conf
,
"corrFuncLength"
,
acf_length
,
1000
);
get_keyval
(
conf
,
"corrFuncStride"
,
acf_stride
,
1
);
if
((
cvm
::
restart_out_freq
%
acf_stride
)
!=
0
)
{
cvm
::
error
(
"Error: corrFuncStride must be commensurate with the restart frequency.
\n
"
,
INPUT_ERROR
);
}
get_keyval
(
conf
,
"corrFuncNormalize"
,
acf_normalize
,
true
);
get_keyval
(
conf
,
"corrFuncOutputFile"
,
acf_outfile
,
std
::
string
(
cvm
::
output_prefix
+
"."
+
this
->
name
+
".corrfunc.dat"
));
}
return
(
cvm
::
get_error
()
?
COLVARS_ERROR
:
COLVARS_OK
);
}
void
colvar
::
setup
()
{
// loop over all components to reset masses of all groups
for
(
size_t
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
for
(
size_t
ig
=
0
;
ig
<
cvcs
[
i
]
->
atom_groups
.
size
();
ig
++
)
{
cvm
::
atom_group
&
atoms
=
*
(
cvcs
[
i
]
->
atom_groups
[
ig
]);
atoms
.
setup
();
atoms
.
reset_mass
(
name
,
i
,
ig
);
atoms
.
read_positions
();
}
}
}
colvar
::~
colvar
()
{
// Clear references to this colvar's cvcs as children
// for dependency purposes
remove_all_children
();
for
(
std
::
vector
<
cvc
*>::
reverse_iterator
ci
=
cvcs
.
rbegin
();
ci
!=
cvcs
.
rend
();
++
ci
)
{
// clear all children of this cvc (i.e. its atom groups)
// because the cvc base class destructor can't do it early enough
// and we don't want to have each cvc derived class do it separately
(
*
ci
)
->
remove_all_children
();
delete
*
ci
;
}
// remove reference to this colvar from the CVM
for
(
std
::
vector
<
colvar
*>::
iterator
cvi
=
cvm
::
colvars
.
begin
();
cvi
!=
cvm
::
colvars
.
end
();
++
cvi
)
{
if
(
*
cvi
==
this
)
{
cvm
::
colvars
.
erase
(
cvi
);
break
;
}
}
}
// ******************** CALC FUNCTIONS ********************
// Default schedule (everything is serialized)
int
colvar
::
calc
()
{
// Note: if anything is added here, it should be added also in the SMP block of calc_colvars()
int
error_code
=
COLVARS_OK
;
if
(
is_enabled
(
f_cv_active
))
{
cvm
::
combine_errors
(
error_code
,
update_cvc_flags
());
cvm
::
combine_errors
(
error_code
,
calc_cvcs
());
cvm
::
combine_errors
(
error_code
,
collect_cvc_data
());
}
return
error_code
;
}
int
colvar
::
calc_cvcs
(
int
first_cvc
,
size_t
num_cvcs
)
{
int
error_code
=
COLVARS_OK
;
if
(
cvm
::
debug
())
cvm
::
log
(
"Calculating colvar
\"
"
+
this
->
name
+
"
\"
, components "
+
cvm
::
to_str
(
first_cvc
)
+
" through "
+
cvm
::
to_str
(
first_cvc
+
num_cvcs
)
+
".
\n
"
);
cvm
::
combine_errors
(
error_code
,
check_cvc_range
(
first_cvc
,
num_cvcs
));
if
(
error_code
!=
COLVARS_OK
)
{
return
error_code
;
}
cvm
::
combine_errors
(
error_code
,
calc_cvc_values
(
first_cvc
,
num_cvcs
));
cvm
::
combine_errors
(
error_code
,
calc_cvc_gradients
(
first_cvc
,
num_cvcs
));
cvm
::
combine_errors
(
error_code
,
calc_cvc_sys_forces
(
first_cvc
,
num_cvcs
));
cvm
::
combine_errors
(
error_code
,
calc_cvc_Jacobians
(
first_cvc
,
num_cvcs
));
if
(
cvm
::
debug
())
cvm
::
log
(
"Done calculating colvar
\"
"
+
this
->
name
+
"
\"
.
\n
"
);
return
error_code
;
}
int
colvar
::
collect_cvc_data
()
{
if
(
cvm
::
debug
())
cvm
::
log
(
"Calculating colvar
\"
"
+
this
->
name
+
"
\"
's properties.
\n
"
);
int
error_code
=
COLVARS_OK
;
cvm
::
combine_errors
(
error_code
,
collect_cvc_values
());
cvm
::
combine_errors
(
error_code
,
collect_cvc_gradients
());
cvm
::
combine_errors
(
error_code
,
collect_cvc_sys_forces
());
cvm
::
combine_errors
(
error_code
,
collect_cvc_Jacobians
());
cvm
::
combine_errors
(
error_code
,
calc_colvar_properties
());
if
(
cvm
::
debug
())
cvm
::
log
(
"Done calculating colvar
\"
"
+
this
->
name
+
"
\"
's properties.
\n
"
);
return
error_code
;
}
int
colvar
::
check_cvc_range
(
int
first_cvc
,
size_t
num_cvcs
)
{
if
((
first_cvc
<
0
)
||
(
first_cvc
>=
((
int
)
cvcs
.
size
())))
{
cvm
::
error
(
"Error: trying to address a component outside the "
"range defined for colvar
\"
"
+
name
+
"
\"
.
\n
"
,
BUG_ERROR
);
return
BUG_ERROR
;
}
return
COLVARS_OK
;
}
int
colvar
::
calc_cvc_values
(
int
first_cvc
,
size_t
num_cvcs
)
{
size_t
const
cvc_max_count
=
num_cvcs
?
num_cvcs
:
num_active_cvcs
();
size_t
i
,
cvc_count
;
// calculate the value of the colvar
if
(
cvm
::
debug
())
cvm
::
log
(
"Calculating colvar components.
\n
"
);
// First, calculate component values
cvm
::
increase_depth
();
for
(
i
=
first_cvc
,
cvc_count
=
0
;
(
i
<
cvcs
.
size
())
&&
(
cvc_count
<
cvc_max_count
);
i
++
)
{
if
(
!
cvcs
[
i
]
->
is_enabled
())
continue
;
cvc_count
++
;
(
cvcs
[
i
])
->
read_data
();
(
cvcs
[
i
])
->
calc_value
();
if
(
cvm
::
debug
())
cvm
::
log
(
"Colvar component no. "
+
cvm
::
to_str
(
i
+
1
)
+
" within colvar
\"
"
+
this
->
name
+
"
\"
has value "
+
cvm
::
to_str
((
cvcs
[
i
])
->
value
(),
cvm
::
cv_width
,
cvm
::
cv_prec
)
+
".
\n
"
);
}
cvm
::
decrease_depth
();
return
COLVARS_OK
;
}
int
colvar
::
collect_cvc_values
()
{
x
.
reset
();
size_t
i
;
// combine them appropriately, using either a scripted function or a polynomial
if
(
is_enabled
(
f_cv_scripted
))
{
// cvcs combined by user script
int
res
=
cvm
::
proxy
->
run_colvar_callback
(
scripted_function
,
sorted_cvc_values
,
x
);
if
(
res
==
COLVARS_NOT_IMPLEMENTED
)
{
cvm
::
error
(
"Scripted colvars are not implemented."
);
return
COLVARS_NOT_IMPLEMENTED
;
}
if
(
res
!=
COLVARS_OK
)
{
cvm
::
error
(
"Error running scripted colvar"
);
return
COLVARS_OK
;
}
}
else
if
(
x
.
type
()
==
colvarvalue
::
type_scalar
)
{
// polynomial combination allowed
for
(
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
if
(
!
cvcs
[
i
]
->
is_enabled
())
continue
;
x
+=
(
cvcs
[
i
])
->
sup_coeff
*
(
((
cvcs
[
i
])
->
sup_np
!=
1
)
?
std
::
pow
((
cvcs
[
i
])
->
value
().
real_value
,
(
cvcs
[
i
])
->
sup_np
)
:
(
cvcs
[
i
])
->
value
().
real_value
);
}
}
else
{
for
(
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
if
(
!
cvcs
[
i
]
->
is_enabled
())
continue
;
x
+=
(
cvcs
[
i
])
->
sup_coeff
*
(
cvcs
[
i
])
->
value
();
}
}
if
(
cvm
::
debug
())
cvm
::
log
(
"Colvar
\"
"
+
this
->
name
+
"
\"
has value "
+
cvm
::
to_str
(
x
,
cvm
::
cv_width
,
cvm
::
cv_prec
)
+
".
\n
"
);
return
COLVARS_OK
;
}
int
colvar
::
calc_cvc_gradients
(
int
first_cvc
,
size_t
num_cvcs
)
{
size_t
const
cvc_max_count
=
num_cvcs
?
num_cvcs
:
num_active_cvcs
();
size_t
i
,
cvc_count
;
if
(
is_enabled
(
f_cv_gradient
))
{
if
(
cvm
::
debug
())
cvm
::
log
(
"Calculating gradients of colvar
\"
"
+
this
->
name
+
"
\"
.
\n
"
);
// calculate the gradients of each component
cvm
::
increase_depth
();
for
(
i
=
first_cvc
,
cvc_count
=
0
;
(
i
<
cvcs
.
size
())
&&
(
cvc_count
<
cvc_max_count
);
i
++
)
{
if
(
!
cvcs
[
i
]
->
is_enabled
())
continue
;
cvc_count
++
;
(
cvcs
[
i
])
->
calc_gradients
();
// if requested, propagate (via chain rule) the gradients above
// to the atoms used to define the roto-translation
for
(
size_t
ig
=
0
;
ig
<
cvcs
[
i
]
->
atom_groups
.
size
();
ig
++
)
{
if
(
cvcs
[
i
]
->
atom_groups
[
ig
]
->
b_fit_gradients
)
cvcs
[
i
]
->
atom_groups
[
ig
]
->
calc_fit_gradients
();
if
(
cvcs
[
i
]
->
is_enabled
(
f_cvc_debug_gradient
))
{
cvm
::
log
(
"Debugging gradients for "
+
cvcs
[
i
]
->
description
);
cvcs
[
i
]
->
debug_gradients
(
cvcs
[
i
]
->
atom_groups
[
ig
]);
}
}
}
cvm
::
decrease_depth
();
if
(
cvm
::
debug
())
cvm
::
log
(
"Done calculating gradients of colvar
\"
"
+
this
->
name
+
"
\"
.
\n
"
);
}
return
COLVARS_OK
;
}
int
colvar
::
collect_cvc_gradients
()
{
size_t
i
;
if
(
is_enabled
(
f_cv_collect_gradient
))
{
if
(
is_enabled
(
f_cv_scripted
))
{
cvm
::
error
(
"Collecting atomic gradients is not implemented for "
"scripted colvars."
,
COLVARS_NOT_IMPLEMENTED
);
return
COLVARS_NOT_IMPLEMENTED
;
}
// Collect the atomic gradients inside colvar object
for
(
unsigned
int
a
=
0
;
a
<
atomic_gradients
.
size
();
a
++
)
{
atomic_gradients
[
a
].
reset
();
}
for
(
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
if
(
!
cvcs
[
i
]
->
is_enabled
())
continue
;
// Coefficient: d(a * x^n) = a * n * x^(n-1) * dx
cvm
::
real
coeff
=
(
cvcs
[
i
])
->
sup_coeff
*
cvm
::
real
((
cvcs
[
i
])
->
sup_np
)
*
std
::
pow
((
cvcs
[
i
])
->
value
().
real_value
,
(
cvcs
[
i
])
->
sup_np
-
1
);
for
(
size_t
j
=
0
;
j
<
cvcs
[
i
]
->
atom_groups
.
size
();
j
++
)
{
cvm
::
atom_group
&
ag
=
*
(
cvcs
[
i
]
->
atom_groups
[
j
]);
// If necessary, apply inverse rotation to get atomic
// gradient in the laboratory frame
if
(
ag
.
b_rotate
)
{
cvm
::
rotation
const
rot_inv
=
ag
.
rot
.
inverse
();
for
(
size_t
k
=
0
;
k
<
ag
.
size
();
k
++
)
{
size_t
a
=
std
::
lower_bound
(
atom_ids
.
begin
(),
atom_ids
.
end
(),
ag
[
k
].
id
)
-
atom_ids
.
begin
();
atomic_gradients
[
a
]
+=
coeff
*
rot_inv
.
rotate
(
ag
[
k
].
grad
);
}
}
else
{
for
(
size_t
k
=
0
;
k
<
ag
.
size
();
k
++
)
{
size_t
a
=
std
::
lower_bound
(
atom_ids
.
begin
(),
atom_ids
.
end
(),
ag
[
k
].
id
)
-
atom_ids
.
begin
();
atomic_gradients
[
a
]
+=
coeff
*
ag
[
k
].
grad
;
}
}
}
}
}
return
COLVARS_OK
;
}
int
colvar
::
calc_cvc_sys_forces
(
int
first_cvc
,
size_t
num_cvcs
)
{
size_t
const
cvc_max_count
=
num_cvcs
?
num_cvcs
:
num_active_cvcs
();
size_t
i
,
cvc_count
;
if
(
is_enabled
(
f_cv_system_force_calc
))
{
if
(
cvm
::
debug
())
cvm
::
log
(
"Calculating system force of colvar
\"
"
+
this
->
name
+
"
\"
.
\n
"
);
// if (!tasks[task_extended_lagrangian] && (cvm::step_relative() > 0)) {
// Disabled check to allow for explicit system force calculation
// even with extended Lagrangian
if
(
cvm
::
step_relative
()
>
0
)
{
cvm
::
increase_depth
();
// get from the cvcs the system forces from the PREVIOUS step
for
(
i
=
first_cvc
,
cvc_count
=
0
;
(
i
<
cvcs
.
size
())
&&
(
cvc_count
<
cvc_max_count
);
i
++
)
{
if
(
!
cvcs
[
i
]
->
is_enabled
())
continue
;
cvc_count
++
;
(
cvcs
[
i
])
->
calc_force_invgrads
();
}
cvm
::
decrease_depth
();
}
if
(
cvm
::
debug
())
cvm
::
log
(
"Done calculating system force of colvar
\"
"
+
this
->
name
+
"
\"
.
\n
"
);
}
return
COLVARS_OK
;
}
int
colvar
::
collect_cvc_sys_forces
()
{
if
(
is_enabled
(
f_cv_system_force_calc
))
{
ft
.
reset
();
if
(
cvm
::
step_relative
()
>
0
)
{
// get from the cvcs the system forces from the PREVIOUS step
for
(
size_t
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
if
(
!
cvcs
[
i
]
->
is_enabled
())
continue
;
// linear combination is assumed
ft
+=
(
cvcs
[
i
])
->
system_force
()
*
(
cvcs
[
i
])
->
sup_coeff
/
active_cvc_square_norm
;
}
}
if
(
!
is_enabled
(
f_cv_hide_Jacobian
))
{
// add the Jacobian force to the system force, and don't apply any silent
// correction internally: biases such as colvarbias_abf will handle it
ft
+=
fj
;
}
}
return
COLVARS_OK
;
}
int
colvar
::
calc_cvc_Jacobians
(
int
first_cvc
,
size_t
num_cvcs
)
{
size_t
const
cvc_max_count
=
num_cvcs
?
num_cvcs
:
num_active_cvcs
();
if
(
is_enabled
(
f_cv_Jacobian
))
{
cvm
::
increase_depth
();
size_t
i
,
cvc_count
;
for
(
i
=
first_cvc
,
cvc_count
=
0
;
(
i
<
cvcs
.
size
())
&&
(
cvc_count
<
cvc_max_count
);
i
++
)
{
if
(
!
cvcs
[
i
]
->
is_enabled
())
continue
;
cvc_count
++
;
(
cvcs
[
i
])
->
calc_Jacobian_derivative
();
}
cvm
::
decrease_depth
();
}
return
COLVARS_OK
;
}
int
colvar
::
collect_cvc_Jacobians
()
{
if
(
is_enabled
(
f_cv_Jacobian
))
{
fj
.
reset
();
for
(
size_t
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
if
(
!
cvcs
[
i
]
->
is_enabled
())
continue
;
// linear combination is assumed
fj
+=
(
cvcs
[
i
])
->
Jacobian_derivative
()
*
(
cvcs
[
i
])
->
sup_coeff
/
active_cvc_square_norm
;
}
fj
*=
cvm
::
boltzmann
()
*
cvm
::
temperature
();
}
return
COLVARS_OK
;
}
int
colvar
::
calc_colvar_properties
()
{
if
(
is_enabled
(
f_cv_fdiff_velocity
))
{
// calculate the velocity by finite differences
if
(
cvm
::
step_relative
()
==
0
)
x_old
=
x
;
else
{
v_fdiff
=
fdiff_velocity
(
x_old
,
x
);
v_reported
=
v_fdiff
;
}
}
if
(
is_enabled
(
f_cv_extended_Lagrangian
))
{
// initialize the restraint center in the first step to the value
// just calculated from the cvcs
// TODO: put it in the restart information
if
(
cvm
::
step_relative
()
==
0
)
{
xr
=
x
;
vr
=
0.0
;
// (already 0; added for clarity)
}
// report the restraint center as "value"
x_reported
=
xr
;
v_reported
=
vr
;
// the "system force" with the extended Lagrangian is just the
// harmonic term acting on the extended coordinate
// Note: this is the force for current timestep
ft_reported
=
(
-
0.5
*
ext_force_k
)
*
this
->
dist2_lgrad
(
xr
,
x
);
}
else
{
x_reported
=
x
;
ft_reported
=
ft
;
}
return
COLVARS_OK
;
}
cvm
::
real
colvar
::
update_forces_energy
()
{
if
(
cvm
::
debug
())
cvm
::
log
(
"Updating colvar
\"
"
+
this
->
name
+
"
\"
.
\n
"
);
// set to zero the applied force
f
.
type
(
value
());
f
.
reset
();
// add the biases' force, which at this point should already have
// been summed over each bias using this colvar
f
+=
fb
;
if
(
is_enabled
(
f_cv_Jacobian
))
{
// the instantaneous Jacobian force was not included in the reported system force;
// instead, it is subtracted from the applied force (silent Jacobian correction)
if
(
is_enabled
(
f_cv_hide_Jacobian
))
f
-=
fj
;
}
if
(
is_enabled
(
f_cv_extended_Lagrangian
))
{
cvm
::
real
dt
=
cvm
::
dt
();
cvm
::
real
f_ext
;
// the total force is applied to the fictitious mass, while the
// atoms only feel the harmonic force
// fr: bias force on extended coordinate (without harmonic spring), for output in trajectory
// f_ext: total force on extended coordinate (including harmonic spring)
// f: - initially, external biasing force
// - after this code block, colvar force to be applied to atomic coordinates, ie. spring force
// (note: wall potential is added to f after this block)
fr
=
f
;
f_ext
=
f
+
(
-
0.5
*
ext_force_k
)
*
this
->
dist2_lgrad
(
xr
,
x
);
f
=
(
-
0.5
*
ext_force_k
)
*
this
->
dist2_rgrad
(
xr
,
x
);
// leapfrog: starting from x_i, f_i, v_(i-1/2)
vr
+=
(
0.5
*
dt
)
*
f_ext
/
ext_mass
;
// Because of leapfrog, kinetic energy at time i is approximate
kinetic_energy
=
0.5
*
ext_mass
*
vr
*
vr
;
potential_energy
=
0.5
*
ext_force_k
*
this
->
dist2
(
xr
,
x
);
// leap to v_(i+1/2)
if
(
is_enabled
(
f_cv_Langevin
))
{
vr
-=
dt
*
ext_gamma
*
vr
.
real_value
;
vr
+=
dt
*
ext_sigma
*
cvm
::
rand_gaussian
()
/
ext_mass
;
}
vr
+=
(
0.5
*
dt
)
*
f_ext
/
ext_mass
;
xr
+=
dt
*
vr
;
xr
.
apply_constraints
();
if
(
this
->
b_periodic
)
this
->
wrap
(
xr
);
}
// Adding wall potential to "true" colvar force, whether or not an extended coordinate is in use
if
(
is_enabled
(
f_cv_lower_wall
)
||
is_enabled
(
f_cv_upper_wall
))
{
// Wall force
colvarvalue
fw
(
x
);
fw
.
reset
();
if
(
cvm
::
debug
())
cvm
::
log
(
"Calculating wall forces for colvar
\"
"
+
this
->
name
+
"
\"
.
\n
"
);
// For a periodic colvar, both walls may be applicable at the same time
// in which case we pick the closer one
if
(
(
!
is_enabled
(
f_cv_upper_wall
))
||
(
this
->
dist2
(
x
,
lower_wall
)
<
this
->
dist2
(
x
,
upper_wall
))
)
{
cvm
::
real
const
grad
=
this
->
dist2_lgrad
(
x
,
lower_wall
);
if
(
grad
<
0.0
)
{
fw
=
-
0.5
*
lower_wall_k
*
grad
;
f
+=
fw
;
if
(
cvm
::
debug
())
cvm
::
log
(
"Applying a lower wall force("
+
cvm
::
to_str
(
fw
)
+
") to
\"
"
+
this
->
name
+
"
\"
.
\n
"
);
}
}
else
{
cvm
::
real
const
grad
=
this
->
dist2_lgrad
(
x
,
upper_wall
);
if
(
grad
>
0.0
)
{
fw
=
-
0.5
*
upper_wall_k
*
grad
;
f
+=
fw
;
if
(
cvm
::
debug
())
cvm
::
log
(
"Applying an upper wall force("
+
cvm
::
to_str
(
fw
)
+
") to
\"
"
+
this
->
name
+
"
\"
.
\n
"
);
}
}
}
f_accumulated
+=
f
;
if
(
is_enabled
(
f_cv_fdiff_velocity
))
{
// set it for the next step
x_old
=
x
;
}
if
(
cvm
::
debug
())
cvm
::
log
(
"Done updating colvar
\"
"
+
this
->
name
+
"
\"
.
\n
"
);
return
(
potential_energy
+
kinetic_energy
);
}
void
colvar
::
communicate_forces
()
{
size_t
i
;
if
(
cvm
::
debug
())
cvm
::
log
(
"Communicating forces from colvar
\"
"
+
this
->
name
+
"
\"
.
\n
"
);
if
(
is_enabled
(
f_cv_scripted
))
{
std
::
vector
<
cvm
::
matrix2d
<
cvm
::
real
>
>
func_grads
;
func_grads
.
reserve
(
cvcs
.
size
());
for
(
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
if
(
!
cvcs
[
i
]
->
is_enabled
())
continue
;
func_grads
.
push_back
(
cvm
::
matrix2d
<
cvm
::
real
>
(
x
.
size
(),
cvcs
[
i
]
->
value
().
size
()));
}
int
res
=
cvm
::
proxy
->
run_colvar_gradient_callback
(
scripted_function
,
sorted_cvc_values
,
func_grads
);
if
(
res
!=
COLVARS_OK
)
{
if
(
res
==
COLVARS_NOT_IMPLEMENTED
)
{
cvm
::
error
(
"Colvar gradient scripts are not implemented."
,
COLVARS_NOT_IMPLEMENTED
);
}
else
{
cvm
::
error
(
"Error running colvar gradient script"
);
}
return
;
}
int
grad_index
=
0
;
// index in the scripted gradients, to account for some components being disabled
for
(
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
if
(
!
cvcs
[
i
]
->
is_enabled
())
continue
;
// cvc force is colvar force times colvar/cvc Jacobian
// (vector-matrix product)
(
cvcs
[
i
])
->
apply_force
(
colvarvalue
(
f_accumulated
.
as_vector
()
*
func_grads
[
grad_index
++
],
cvcs
[
i
]
->
value
().
type
()));
}
}
else
if
(
x
.
type
()
==
colvarvalue
::
type_scalar
)
{
for
(
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
if
(
!
cvcs
[
i
]
->
is_enabled
())
continue
;
(
cvcs
[
i
])
->
apply_force
(
f_accumulated
*
(
cvcs
[
i
])
->
sup_coeff
*
cvm
::
real
((
cvcs
[
i
])
->
sup_np
)
*
(
std
::
pow
((
cvcs
[
i
])
->
value
().
real_value
,
(
cvcs
[
i
])
->
sup_np
-
1
))
);
}
}
else
{
for
(
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
if
(
!
cvcs
[
i
]
->
is_enabled
())
continue
;
(
cvcs
[
i
])
->
apply_force
(
f_accumulated
*
(
cvcs
[
i
])
->
sup_coeff
);
}
}
// Accumulated forces have been applied, impulse-style
// Reset to start accumulating again
f_accumulated
.
reset
();
if
(
cvm
::
debug
())
cvm
::
log
(
"Done communicating forces from colvar
\"
"
+
this
->
name
+
"
\"
.
\n
"
);
}
int
colvar
::
set_cvc_flags
(
std
::
vector
<
bool
>
const
&
flags
)
{
if
(
flags
.
size
()
!=
cvcs
.
size
())
{
cvm
::
error
(
"ERROR: Wrong number of CVC flags provided."
);
return
COLVARS_ERROR
;
}
// We cannot enable or disable cvcs in the middle of a timestep or colvar evaluation sequence
// so we store the flags that will be enforced at the next call to calc()
cvc_flags
=
flags
;
return
COLVARS_OK
;
}
int
colvar
::
update_cvc_flags
()
{
// Update the enabled/disabled status of cvcs if necessary
if
(
cvc_flags
.
size
())
{
n_active_cvcs
=
0
;
active_cvc_square_norm
=
0.
;
for
(
size_t
i
=
0
;
i
<
cvcs
.
size
();
i
++
)
{
cvcs
[
i
]
->
feature_states
[
f_cvc_active
]
->
enabled
=
cvc_flags
[
i
];
if
(
cvcs
[
i
]
->
is_enabled
())
{
n_active_cvcs
++
;
active_cvc_square_norm
+=
cvcs
[
i
]
->
sup_coeff
*
cvcs
[
i
]
->
sup_coeff
;
}
}
if
(
!
n_active_cvcs
)
{
cvm
::
error
(
"ERROR: All CVCs are disabled for colvar "
+
this
->
name
+
"
\n
"
);
return
COLVARS_ERROR
;
}
cvc_flags
.
resize
(
0
);
}
return
COLVARS_OK
;
}
// ******************** METRIC FUNCTIONS ********************
// Use the metrics defined by \link cvc \endlink objects
bool
colvar
::
periodic_boundaries
(
colvarvalue
const
&
lb
,
colvarvalue
const
&
ub
)
const
{
if
(
(
!
is_enabled
(
f_cv_lower_boundary
))
||
(
!
is_enabled
(
f_cv_upper_boundary
))
)
{
cvm
::
log
(
"Error: checking periodicity for collective variable
\"
"
+
this
->
name
+
"
\"
"
"requires lower and upper boundaries to be defined.
\n
"
);
cvm
::
set_error_bit
(
INPUT_ERROR
);
}
if
(
period
>
0.0
)
{
if
(
((
std
::
sqrt
(
this
->
dist2
(
lb
,
ub
)))
/
this
->
width
)
<
1.0E-10
)
{
return
true
;
}
}
return
false
;
}
bool
colvar
::
periodic_boundaries
()
const
{
if
(
(
!
is_enabled
(
f_cv_lower_boundary
))
||
(
!
is_enabled
(
f_cv_upper_boundary
))
)
{
cvm
::
log
(
"Error: checking periodicity for collective variable
\"
"
+
this
->
name
+
"
\"
"
"requires lower and upper boundaries to be defined.
\n
"
);
}
return
periodic_boundaries
(
lower_boundary
,
upper_boundary
);
}
cvm
::
real
colvar
::
dist2
(
colvarvalue
const
&
x1
,
colvarvalue
const
&
x2
)
const
{
if
(
is_enabled
(
f_cv_homogeneous
))
{
return
(
cvcs
[
0
])
->
dist2
(
x1
,
x2
);
}
else
{
return
x1
.
dist2
(
x2
);
}
}
colvarvalue
colvar
::
dist2_lgrad
(
colvarvalue
const
&
x1
,
colvarvalue
const
&
x2
)
const
{
if
(
is_enabled
(
f_cv_homogeneous
))
{
return
(
cvcs
[
0
])
->
dist2_lgrad
(
x1
,
x2
);
}
else
{
return
x1
.
dist2_grad
(
x2
);
}
}
colvarvalue
colvar
::
dist2_rgrad
(
colvarvalue
const
&
x1
,
colvarvalue
const
&
x2
)
const
{
if
(
is_enabled
(
f_cv_homogeneous
))
{
return
(
cvcs
[
0
])
->
dist2_rgrad
(
x1
,
x2
);
}
else
{
return
x2
.
dist2_grad
(
x1
);
}
}
void
colvar
::
wrap
(
colvarvalue
&
x
)
const
{
if
(
is_enabled
(
f_cv_homogeneous
))
{
(
cvcs
[
0
])
->
wrap
(
x
);
}
return
;
}
// ******************** INPUT FUNCTIONS ********************
std
::
istream
&
colvar
::
read_restart
(
std
::
istream
&
is
)
{
size_t
const
start_pos
=
is
.
tellg
();
std
::
string
conf
;
if
(
!
(
is
>>
colvarparse
::
read_block
(
"colvar"
,
conf
))
)
{
// this is not a colvar block
is
.
clear
();
is
.
seekg
(
start_pos
,
std
::
ios
::
beg
);
is
.
setstate
(
std
::
ios
::
failbit
);
return
is
;
}
{
std
::
string
check_name
=
""
;
if
(
(
get_keyval
(
conf
,
"name"
,
check_name
,
std
::
string
(
""
),
colvarparse
::
parse_silent
))
&&
(
check_name
!=
name
)
)
{
cvm
::
error
(
"Error: the state file does not match the "
"configuration file, at colvar
\"
"
+
name
+
"
\"
.
\n
"
);
}
if
(
check_name
.
size
()
==
0
)
{
cvm
::
error
(
"Error: Collective variable in the "
"restart file without any identifier.
\n
"
);
}
}
if
(
!
(
get_keyval
(
conf
,
"x"
,
x
,
colvarvalue
(
x
.
type
()),
colvarparse
::
parse_silent
))
)
{
cvm
::
log
(
"Error: restart file does not contain "
"the value of the colvar
\"
"
+
name
+
"
\"
.
\n
"
);
}
else
{
cvm
::
log
(
"Restarting collective variable
\"
"
+
name
+
"
\"
from value: "
+
cvm
::
to_str
(
x
)
+
"
\n
"
);
}
if
(
is_enabled
(
f_cv_extended_Lagrangian
))
{
if
(
!
(
get_keyval
(
conf
,
"extended_x"
,
xr
,
colvarvalue
(
x
.
type
()),
colvarparse
::
parse_silent
))
&&
!
(
get_keyval
(
conf
,
"extended_v"
,
vr
,
colvarvalue
(
x
.
type
()),
colvarparse
::
parse_silent
))
)
{
cvm
::
log
(
"Error: restart file does not contain "
"
\"
extended_x
\"
or
\"
extended_v
\"
for the colvar
\"
"
+
name
+
"
\"
, but you requested
\"
extendedLagrangian
\"
.
\n
"
);
}
x_reported
=
xr
;
}
else
{
x_reported
=
x
;
}
if
(
is_enabled
(
f_cv_output_velocity
))
{
if
(
!
(
get_keyval
(
conf
,
"v"
,
v_fdiff
,
colvarvalue
(
x
.
type
()),
colvarparse
::
parse_silent
))
)
{
cvm
::
log
(
"Error: restart file does not contain "
"the velocity for the colvar
\"
"
+
name
+
"
\"
, but you requested
\"
outputVelocity
\"
.
\n
"
);
}
if
(
is_enabled
(
f_cv_extended_Lagrangian
))
{
v_reported
=
vr
;
}
else
{
v_reported
=
v_fdiff
;
}
}
return
is
;
}
std
::
istream
&
colvar
::
read_traj
(
std
::
istream
&
is
)
{
size_t
const
start_pos
=
is
.
tellg
();
if
(
is_enabled
(
f_cv_output_value
))
{
if
(
!
(
is
>>
x
))
{
cvm
::
log
(
"Error: in reading the value of colvar
\"
"
+
this
->
name
+
"
\"
from trajectory.
\n
"
);
is
.
clear
();
is
.
seekg
(
start_pos
,
std
::
ios
::
beg
);
is
.
setstate
(
std
::
ios
::
failbit
);
return
is
;
}
if
(
is_enabled
(
f_cv_extended_Lagrangian
))
{
is
>>
xr
;
x_reported
=
xr
;
}
else
{
x_reported
=
x
;
}
}
if
(
is_enabled
(
f_cv_output_velocity
))
{
is
>>
v_fdiff
;
if
(
is_enabled
(
f_cv_extended_Lagrangian
))
{
is
>>
vr
;
v_reported
=
vr
;
}
else
{
v_reported
=
v_fdiff
;
}
}
if
(
is_enabled
(
f_cv_output_system_force
))
{
is
>>
ft
;
ft_reported
=
ft
;
}
if
(
is_enabled
(
f_cv_output_applied_force
))
{
is
>>
f
;
}
return
is
;
}
// ******************** OUTPUT FUNCTIONS ********************
std
::
ostream
&
colvar
::
write_restart
(
std
::
ostream
&
os
)
{
os
<<
"colvar {
\n
"
<<
" name "
<<
name
<<
"
\n
"
<<
" x "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
x
<<
"
\n
"
;
if
(
is_enabled
(
f_cv_output_velocity
))
{
os
<<
" v "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
v_reported
<<
"
\n
"
;
}
if
(
is_enabled
(
f_cv_extended_Lagrangian
))
{
os
<<
" extended_x "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
xr
<<
"
\n
"
<<
" extended_v "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
vr
<<
"
\n
"
;
}
os
<<
"}
\n\n
"
;
return
os
;
}
std
::
ostream
&
colvar
::
write_traj_label
(
std
::
ostream
&
os
)
{
size_t
const
this_cv_width
=
x
.
output_width
(
cvm
::
cv_width
);
os
<<
" "
;
if
(
is_enabled
(
f_cv_output_value
))
{
os
<<
" "
<<
cvm
::
wrap_string
(
this
->
name
,
this_cv_width
);
if
(
is_enabled
(
f_cv_extended_Lagrangian
))
{
// extended DOF
os
<<
" r_"
<<
cvm
::
wrap_string
(
this
->
name
,
this_cv_width
-
2
);
}
}
if
(
is_enabled
(
f_cv_output_velocity
))
{
os
<<
" v_"
<<
cvm
::
wrap_string
(
this
->
name
,
this_cv_width
-
2
);
if
(
is_enabled
(
f_cv_extended_Lagrangian
))
{
// extended DOF
os
<<
" vr_"
<<
cvm
::
wrap_string
(
this
->
name
,
this_cv_width
-
3
);
}
}
if
(
is_enabled
(
f_cv_output_energy
))
{
os
<<
" Ep_"
<<
cvm
::
wrap_string
(
this
->
name
,
this_cv_width
-
3
)
<<
" Ek_"
<<
cvm
::
wrap_string
(
this
->
name
,
this_cv_width
-
3
);
}
if
(
is_enabled
(
f_cv_output_system_force
))
{
os
<<
" fs_"
<<
cvm
::
wrap_string
(
this
->
name
,
this_cv_width
-
3
);
}
if
(
is_enabled
(
f_cv_output_applied_force
))
{
os
<<
" fa_"
<<
cvm
::
wrap_string
(
this
->
name
,
this_cv_width
-
3
);
}
return
os
;
}
std
::
ostream
&
colvar
::
write_traj
(
std
::
ostream
&
os
)
{
os
<<
" "
;
if
(
is_enabled
(
f_cv_output_value
))
{
if
(
is_enabled
(
f_cv_extended_Lagrangian
))
{
os
<<
" "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
x
;
}
os
<<
" "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
x_reported
;
}
if
(
is_enabled
(
f_cv_output_velocity
))
{
if
(
is_enabled
(
f_cv_extended_Lagrangian
))
{
os
<<
" "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
v_fdiff
;
}
os
<<
" "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
v_reported
;
}
if
(
is_enabled
(
f_cv_output_energy
))
{
os
<<
" "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
potential_energy
<<
" "
<<
kinetic_energy
;
}
if
(
is_enabled
(
f_cv_output_system_force
))
{
os
<<
" "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
ft_reported
;
}
if
(
is_enabled
(
f_cv_output_applied_force
))
{
if
(
is_enabled
(
f_cv_extended_Lagrangian
))
{
os
<<
" "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
fr
;
}
else
{
os
<<
" "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
f
;
}
}
return
os
;
}
int
colvar
::
write_output_files
()
{
if
(
cvm
::
b_analysis
)
{
if
(
acf
.
size
())
{
cvm
::
log
(
"Writing acf to file
\"
"
+
acf_outfile
+
"
\"
.
\n
"
);
cvm
::
backup_file
(
acf_outfile
.
c_str
());
cvm
::
ofstream
acf_os
(
acf_outfile
.
c_str
());
if
(
!
acf_os
.
is_open
())
{
cvm
::
error
(
"Cannot open file
\"
"
+
acf_outfile
+
"
\"
.
\n
"
,
FILE_ERROR
);
}
write_acf
(
acf_os
);
acf_os
.
close
();
}
if
(
runave_os
.
is_open
())
{
runave_os
.
close
();
}
}
return
(
cvm
::
get_error
()
?
COLVARS_ERROR
:
COLVARS_OK
);
}
// ******************** ANALYSIS FUNCTIONS ********************
void
colvar
::
analyze
()
{
if
(
is_enabled
(
f_cv_runave
))
{
calc_runave
();
}
if
(
is_enabled
(
f_cv_corrfunc
))
{
calc_acf
();
}
}
inline
void
history_add_value
(
size_t
const
&
history_length
,
std
::
list
<
colvarvalue
>
&
history
,
colvarvalue
const
&
new_value
)
{
history
.
push_front
(
new_value
);
if
(
history
.
size
()
>
history_length
)
history
.
pop_back
();
}
inline
void
history_incr
(
std
::
list
<
std
::
list
<
colvarvalue
>
>
&
history
,
std
::
list
<
std
::
list
<
colvarvalue
>
>::
iterator
&
history_p
)
{
if
((
++
history_p
)
==
history
.
end
())
history_p
=
history
.
begin
();
}
int
colvar
::
calc_acf
()
{
// using here an acf_stride-long list of vectors for either
// coordinates(acf_x_history) or velocities (acf_v_history); each vector can
// contain up to acf_length values, which are contiguous in memory
// representation but separated by acf_stride in the time series;
// the pointer to each vector is changed at every step
if
(
acf_x_history
.
empty
()
&&
acf_v_history
.
empty
())
{
// first-step operations
colvar
*
cfcv
=
(
acf_colvar_name
.
size
()
?
cvm
::
colvar_by_name
(
acf_colvar_name
)
:
this
);
if
(
colvarvalue
::
check_types
(
cfcv
->
value
(),
value
()))
{
cvm
::
error
(
"Error: correlation function between
\"
"
+
cfcv
->
name
+
"
\"
and
\"
"
+
this
->
name
+
"
\"
cannot be calculated, "
"because their value types are different.
\n
"
,
INPUT_ERROR
);
}
acf_nframes
=
0
;
cvm
::
log
(
"Colvar
\"
"
+
this
->
name
+
"
\"
: initializing ACF calculation.
\n
"
);
if
(
acf
.
size
()
<
acf_length
+
1
)
acf
.
resize
(
acf_length
+
1
,
0.0
);
size_t
i
;
switch
(
acf_type
)
{
case
acf_vel:
// allocate space for the velocities history
for
(
i
=
0
;
i
<
acf_stride
;
i
++
)
{
acf_v_history
.
push_back
(
std
::
list
<
colvarvalue
>
());
}
acf_v_history_p
=
acf_v_history
.
begin
();
break
;
case
acf_coor:
case
acf_p2coor:
// allocate space for the coordinates history
for
(
i
=
0
;
i
<
acf_stride
;
i
++
)
{
acf_x_history
.
push_back
(
std
::
list
<
colvarvalue
>
());
}
acf_x_history_p
=
acf_x_history
.
begin
();
break
;
default
:
break
;
}
}
else
{
colvar
*
cfcv
=
(
acf_colvar_name
.
size
()
?
cvm
::
colvar_by_name
(
acf_colvar_name
)
:
this
);
switch
(
acf_type
)
{
case
acf_vel:
if
(
is_enabled
(
f_cv_fdiff_velocity
))
{
// calc() should do this already, but this only happens in a
// simulation; better do it again in case a trajectory is
// being read
v_reported
=
v_fdiff
=
fdiff_velocity
(
x_old
,
cfcv
->
value
());
}
calc_vel_acf
((
*
acf_v_history_p
),
cfcv
->
velocity
());
// store this value in the history
history_add_value
(
acf_length
+
acf_offset
,
*
acf_v_history_p
,
cfcv
->
velocity
());
// if stride is larger than one, cycle among different histories
history_incr
(
acf_v_history
,
acf_v_history_p
);
break
;
case
acf_coor:
calc_coor_acf
((
*
acf_x_history_p
),
cfcv
->
value
());
history_add_value
(
acf_length
+
acf_offset
,
*
acf_x_history_p
,
cfcv
->
value
());
history_incr
(
acf_x_history
,
acf_x_history_p
);
break
;
case
acf_p2coor:
calc_p2coor_acf
((
*
acf_x_history_p
),
cfcv
->
value
());
history_add_value
(
acf_length
+
acf_offset
,
*
acf_x_history_p
,
cfcv
->
value
());
history_incr
(
acf_x_history
,
acf_x_history_p
);
break
;
default
:
break
;
}
}
if
(
is_enabled
(
f_cv_fdiff_velocity
))
{
// set it for the next step
x_old
=
x
;
}
return
(
cvm
::
get_error
()
?
COLVARS_ERROR
:
COLVARS_OK
);
}
int
colvar
::
calc_vel_acf
(
std
::
list
<
colvarvalue
>
&
v_list
,
colvarvalue
const
&
v
)
{
// loop over stored velocities and add to the ACF, but only the
// length is sufficient to hold an entire row of ACF values
if
(
v_list
.
size
()
>=
acf_length
+
acf_offset
)
{
std
::
list
<
colvarvalue
>::
iterator
vs_i
=
v_list
.
begin
();
std
::
vector
<
cvm
::
real
>::
iterator
acf_i
=
acf
.
begin
();
for
(
size_t
i
=
0
;
i
<
acf_offset
;
i
++
)
++
vs_i
;
// current vel with itself
*
(
acf_i
)
+=
v
.
norm2
();
++
acf_i
;
// inner products of previous velocities with current (acf_i and
// vs_i are updated)
colvarvalue
::
inner_opt
(
v
,
vs_i
,
v_list
.
end
(),
acf_i
);
acf_nframes
++
;
}
return
(
cvm
::
get_error
()
?
COLVARS_ERROR
:
COLVARS_OK
);
}
void
colvar
::
calc_coor_acf
(
std
::
list
<
colvarvalue
>
&
x_list
,
colvarvalue
const
&
x
)
{
// same as above but for coordinates
if
(
x_list
.
size
()
>=
acf_length
+
acf_offset
)
{
std
::
list
<
colvarvalue
>::
iterator
xs_i
=
x_list
.
begin
();
std
::
vector
<
cvm
::
real
>::
iterator
acf_i
=
acf
.
begin
();
for
(
size_t
i
=
0
;
i
<
acf_offset
;
i
++
)
++
xs_i
;
*
(
acf_i
++
)
+=
x
.
norm2
();
colvarvalue
::
inner_opt
(
x
,
xs_i
,
x_list
.
end
(),
acf_i
);
acf_nframes
++
;
}
}
void
colvar
::
calc_p2coor_acf
(
std
::
list
<
colvarvalue
>
&
x_list
,
colvarvalue
const
&
x
)
{
// same as above but with second order Legendre polynomial instead
// of just the scalar product
if
(
x_list
.
size
()
>=
acf_length
+
acf_offset
)
{
std
::
list
<
colvarvalue
>::
iterator
xs_i
=
x_list
.
begin
();
std
::
vector
<
cvm
::
real
>::
iterator
acf_i
=
acf
.
begin
();
for
(
size_t
i
=
0
;
i
<
acf_offset
;
i
++
)
++
xs_i
;
// value of P2(0) = 1
*
(
acf_i
++
)
+=
1.0
;
colvarvalue
::
p2leg_opt
(
x
,
xs_i
,
x_list
.
end
(),
acf_i
);
acf_nframes
++
;
}
}
void
colvar
::
write_acf
(
std
::
ostream
&
os
)
{
if
(
!
acf_nframes
)
cvm
::
log
(
"Warning: ACF was not calculated (insufficient frames).
\n
"
);
os
.
setf
(
std
::
ios
::
scientific
,
std
::
ios
::
floatfield
);
os
<<
"# Autocorrelation function for collective variable
\"
"
<<
this
->
name
<<
"
\"\n
"
;
// one frame is used for normalization, the statistical sample is
// hence decreased
os
<<
"# nframes = "
<<
(
acf_normalize
?
acf_nframes
-
1
:
acf_nframes
)
<<
"
\n
"
;
cvm
::
real
const
acf_norm
=
acf
.
front
()
/
cvm
::
real
(
acf_nframes
);
std
::
vector
<
cvm
::
real
>::
iterator
acf_i
;
size_t
it
=
acf_offset
;
for
(
acf_i
=
acf
.
begin
();
acf_i
!=
acf
.
end
();
++
acf_i
)
{
os
<<
std
::
setw
(
cvm
::
it_width
)
<<
acf_stride
*
(
it
++
)
<<
" "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
(
acf_normalize
?
(
*
acf_i
)
/
(
acf_norm
*
cvm
::
real
(
acf_nframes
))
:
(
*
acf_i
)
/
(
cvm
::
real
(
acf_nframes
))
)
<<
"
\n
"
;
}
}
void
colvar
::
calc_runave
()
{
if
(
x_history
.
empty
())
{
runave
.
type
(
value
().
type
());
runave
.
reset
();
// first-step operations
if
(
cvm
::
debug
())
cvm
::
log
(
"Colvar
\"
"
+
this
->
name
+
"
\"
: initializing running average calculation.
\n
"
);
acf_nframes
=
0
;
x_history
.
push_back
(
std
::
list
<
colvarvalue
>
());
x_history_p
=
x_history
.
begin
();
}
else
{
if
(
(
cvm
::
step_relative
()
%
runave_stride
)
==
0
)
{
if
((
*
x_history_p
).
size
()
>=
runave_length
-
1
)
{
runave
=
x
;
std
::
list
<
colvarvalue
>::
iterator
xs_i
;
for
(
xs_i
=
(
*
x_history_p
).
begin
();
xs_i
!=
(
*
x_history_p
).
end
();
++
xs_i
)
{
runave
+=
(
*
xs_i
);
}
runave
*=
1.0
/
cvm
::
real
(
runave_length
);
runave
.
apply_constraints
();
runave_variance
=
0.0
;
runave_variance
+=
this
->
dist2
(
x
,
runave
);
for
(
xs_i
=
(
*
x_history_p
).
begin
();
xs_i
!=
(
*
x_history_p
).
end
();
++
xs_i
)
{
runave_variance
+=
this
->
dist2
(
x
,
(
*
xs_i
));
}
runave_variance
*=
1.0
/
cvm
::
real
(
runave_length
-
1
);
runave_os
<<
std
::
setw
(
cvm
::
it_width
)
<<
cvm
::
step_relative
()
<<
" "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
runave
<<
" "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
std
::
sqrt
(
runave_variance
)
<<
"
\n
"
;
}
history_add_value
(
runave_length
,
*
x_history_p
,
x
);
}
}
}
// Static members
std
::
vector
<
cvm
::
deps
::
feature
*>
colvar
::
cv_features
;
Event Timeline
Log In to Comment