Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65072915
colvarbias_restraint.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, May 31, 14:17
Size
40 KB
Mime Type
text/x-c
Expires
Sun, Jun 2, 14:17 (2 d)
Engine
blob
Format
Raw Data
Handle
17999118
Attached To
rLAMMPS lammps
colvarbias_restraint.cpp
View Options
// -*- c++ -*-
// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/colvars/colvars
// Please update all Colvars source files before making any changes.
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarbias_restraint.h"
colvarbias_restraint
::
colvarbias_restraint
(
char
const
*
key
)
:
colvarbias
(
key
)
{
}
int
colvarbias_restraint
::
init
(
std
::
string
const
&
conf
)
{
colvarbias
::
init
(
conf
);
enable
(
f_cvb_apply_force
);
if
(
cvm
::
debug
())
cvm
::
log
(
"Initializing a new restraint bias.
\n
"
);
return
COLVARS_OK
;
}
int
colvarbias_restraint
::
update
()
{
// Update base class (bias_energy and colvar_forces are zeroed there)
colvarbias
::
update
();
// Force and energy calculation
for
(
size_t
i
=
0
;
i
<
num_variables
();
i
++
)
{
bias_energy
+=
restraint_potential
(
i
);
colvar_forces
[
i
].
type
(
variables
(
i
)
->
value
());
colvar_forces
[
i
].
is_derivative
();
colvar_forces
[
i
]
=
restraint_force
(
i
);
}
if
(
cvm
::
debug
())
cvm
::
log
(
"Done updating the restraint bias
\"
"
+
this
->
name
+
"
\"
.
\n
"
);
if
(
cvm
::
debug
())
cvm
::
log
(
"Current forces for the restraint bias
\"
"
+
this
->
name
+
"
\"
: "
+
cvm
::
to_str
(
colvar_forces
)
+
".
\n
"
);
return
COLVARS_OK
;
}
colvarbias_restraint
::~
colvarbias_restraint
()
{
}
std
::
string
const
colvarbias_restraint
::
get_state_params
()
const
{
return
colvarbias
::
get_state_params
();
}
int
colvarbias_restraint
::
set_state_params
(
std
::
string
const
&
conf
)
{
return
colvarbias
::
set_state_params
(
conf
);
}
std
::
ostream
&
colvarbias_restraint
::
write_traj_label
(
std
::
ostream
&
os
)
{
return
colvarbias
::
write_traj_label
(
os
);
}
std
::
ostream
&
colvarbias_restraint
::
write_traj
(
std
::
ostream
&
os
)
{
return
colvarbias
::
write_traj
(
os
);
}
colvarbias_restraint_centers
::
colvarbias_restraint_centers
(
char
const
*
key
)
:
colvarbias
(
key
),
colvarbias_restraint
(
key
)
{
}
int
colvarbias_restraint_centers
::
init
(
std
::
string
const
&
conf
)
{
size_t
i
;
bool
null_centers
=
(
colvar_centers
.
size
()
==
0
);
if
(
null_centers
)
{
// try to initialize the restraint centers for the first time
colvar_centers
.
resize
(
num_variables
());
colvar_centers_raw
.
resize
(
num_variables
());
for
(
i
=
0
;
i
<
num_variables
();
i
++
)
{
colvar_centers
[
i
].
type
(
variables
(
i
)
->
value
());
colvar_centers
[
i
].
reset
();
colvar_centers_raw
[
i
].
type
(
variables
(
i
)
->
value
());
colvar_centers_raw
[
i
].
reset
();
}
}
if
(
get_keyval
(
conf
,
"centers"
,
colvar_centers
,
colvar_centers
))
{
for
(
i
=
0
;
i
<
num_variables
();
i
++
)
{
if
(
cvm
::
debug
())
{
cvm
::
log
(
"colvarbias_restraint: parsing initial centers, i = "
+
cvm
::
to_str
(
i
)
+
".
\n
"
);
}
colvar_centers_raw
[
i
]
=
colvar_centers
[
i
];
colvar_centers
[
i
].
apply_constraints
();
}
null_centers
=
false
;
}
if
(
null_centers
)
{
colvar_centers
.
clear
();
cvm
::
error
(
"Error: must define the initial centers of the restraints.
\n
"
,
INPUT_ERROR
);
return
INPUT_ERROR
;
}
if
(
colvar_centers
.
size
()
!=
num_variables
())
{
cvm
::
error
(
"Error: number of centers does not match "
"that of collective variables.
\n
"
,
INPUT_ERROR
);
return
INPUT_ERROR
;
}
return
COLVARS_OK
;
}
int
colvarbias_restraint_centers
::
change_configuration
(
std
::
string
const
&
conf
)
{
if
(
get_keyval
(
conf
,
"centers"
,
colvar_centers
,
colvar_centers
))
{
for
(
size_t
i
=
0
;
i
<
num_variables
();
i
++
)
{
colvar_centers
[
i
].
type
(
variables
(
i
)
->
value
());
colvar_centers
[
i
].
apply_constraints
();
colvar_centers_raw
[
i
].
type
(
variables
(
i
)
->
value
());
colvar_centers_raw
[
i
]
=
colvar_centers
[
i
];
}
}
return
COLVARS_OK
;
}
colvarbias_restraint_k
::
colvarbias_restraint_k
(
char
const
*
key
)
:
colvarbias
(
key
),
colvarbias_restraint
(
key
)
{
force_k
=
-
1.0
;
}
int
colvarbias_restraint_k
::
init
(
std
::
string
const
&
conf
)
{
get_keyval
(
conf
,
"forceConstant"
,
force_k
,
(
force_k
>
0.0
?
force_k
:
1.0
));
if
(
force_k
<
0.0
)
{
cvm
::
error
(
"Error: undefined or invalid force constant.
\n
"
,
INPUT_ERROR
);
return
INPUT_ERROR
;
}
return
COLVARS_OK
;
}
int
colvarbias_restraint_k
::
change_configuration
(
std
::
string
const
&
conf
)
{
get_keyval
(
conf
,
"forceConstant"
,
force_k
,
force_k
);
return
COLVARS_OK
;
}
colvarbias_restraint_moving
::
colvarbias_restraint_moving
(
char
const
*
key
)
{
target_nstages
=
0
;
target_nsteps
=
0
;
stage
=
0
;
b_chg_centers
=
false
;
b_chg_force_k
=
false
;
}
int
colvarbias_restraint_moving
::
init
(
std
::
string
const
&
conf
)
{
if
(
b_chg_centers
&&
b_chg_force_k
)
{
cvm
::
error
(
"Error: cannot specify both targetCenters and targetForceConstant.
\n
"
,
INPUT_ERROR
);
return
INPUT_ERROR
;
}
if
(
b_chg_centers
||
b_chg_force_k
)
{
get_keyval
(
conf
,
"targetNumSteps"
,
target_nsteps
,
target_nsteps
);
if
(
!
target_nsteps
)
{
cvm
::
error
(
"Error: targetNumSteps must be non-zero.
\n
"
,
INPUT_ERROR
);
return
cvm
::
get_error
();
}
if
(
get_keyval
(
conf
,
"targetNumStages"
,
target_nstages
,
target_nstages
)
&&
lambda_schedule
.
size
())
{
cvm
::
error
(
"Error: targetNumStages and lambdaSchedule are incompatible.
\n
"
,
INPUT_ERROR
);
return
cvm
::
get_error
();
}
}
return
COLVARS_OK
;
}
std
::
string
const
colvarbias_restraint_moving
::
get_state_params
()
const
{
std
::
ostringstream
os
;
os
.
setf
(
std
::
ios
::
scientific
,
std
::
ios
::
floatfield
);
if
(
b_chg_centers
||
b_chg_force_k
)
{
// TODO move this
if
(
target_nstages
)
{
os
<<
"stage "
<<
std
::
setw
(
cvm
::
it_width
)
<<
stage
<<
"
\n
"
;
}
}
return
os
.
str
();
}
int
colvarbias_restraint_moving
::
set_state_params
(
std
::
string
const
&
conf
)
{
if
(
b_chg_centers
||
b_chg_force_k
)
{
if
(
target_nstages
)
{
// cvm::log ("Reading current stage from the restart.\n");
if
(
!
get_keyval
(
conf
,
"stage"
,
stage
))
cvm
::
error
(
"Error: current stage is missing from the restart.
\n
"
);
}
}
return
COLVARS_OK
;
}
colvarbias_restraint_centers_moving
::
colvarbias_restraint_centers_moving
(
char
const
*
key
)
:
colvarbias
(
key
),
colvarbias_restraint
(
key
),
colvarbias_restraint_centers
(
key
),
colvarbias_restraint_moving
(
key
)
{
b_chg_centers
=
false
;
b_output_centers
=
false
;
b_output_acc_work
=
false
;
acc_work
=
0.0
;
}
int
colvarbias_restraint_centers_moving
::
init
(
std
::
string
const
&
conf
)
{
colvarbias_restraint_centers
::
init
(
conf
);
if
(
cvm
::
debug
())
{
cvm
::
log
(
"colvarbias_restraint: parsing target centers.
\n
"
);
}
size_t
i
;
if
(
get_keyval
(
conf
,
"targetCenters"
,
target_centers
,
colvar_centers
))
{
if
(
colvar_centers
.
size
()
!=
num_variables
())
{
cvm
::
error
(
"Error: number of target centers does not match "
"that of collective variables.
\n
"
);
}
b_chg_centers
=
true
;
for
(
i
=
0
;
i
<
target_centers
.
size
();
i
++
)
{
target_centers
[
i
].
apply_constraints
();
}
}
if
(
b_chg_centers
)
{
// parse moving restraint options
colvarbias_restraint_moving
::
init
(
conf
);
}
else
{
target_centers
.
clear
();
return
COLVARS_OK
;
}
get_keyval
(
conf
,
"outputCenters"
,
b_output_centers
,
b_output_centers
);
get_keyval
(
conf
,
"outputAccumulatedWork"
,
b_output_acc_work
,
b_output_acc_work
);
return
COLVARS_OK
;
}
int
colvarbias_restraint_centers_moving
::
update
()
{
if
(
b_chg_centers
)
{
if
(
cvm
::
debug
())
{
cvm
::
log
(
"Updating centers for the restraint bias
\"
"
+
this
->
name
+
"
\"
: "
+
cvm
::
to_str
(
colvar_centers
)
+
".
\n
"
);
}
if
(
!
centers_incr
.
size
())
{
// if this is the first calculation, calculate the advancement
// at each simulation step (or stage, if applicable)
// (take current stage into account: it can be non-zero
// if we are restarting a staged calculation)
centers_incr
.
resize
(
num_variables
());
for
(
size_t
i
=
0
;
i
<
num_variables
();
i
++
)
{
centers_incr
[
i
].
type
(
variables
(
i
)
->
value
());
centers_incr
[
i
]
=
(
target_centers
[
i
]
-
colvar_centers_raw
[
i
])
/
cvm
::
real
(
target_nstages
?
(
target_nstages
-
stage
)
:
(
target_nsteps
-
cvm
::
step_absolute
()));
}
if
(
cvm
::
debug
())
{
cvm
::
log
(
"Center increment for the restraint bias
\"
"
+
this
->
name
+
"
\"
: "
+
cvm
::
to_str
(
centers_incr
)
+
" at stage "
+
cvm
::
to_str
(
stage
)
+
".
\n
"
);
}
}
if
(
target_nstages
)
{
if
((
cvm
::
step_relative
()
>
0
)
&&
(
cvm
::
step_absolute
()
%
target_nsteps
)
==
0
&&
stage
<
target_nstages
)
{
for
(
size_t
i
=
0
;
i
<
num_variables
();
i
++
)
{
colvar_centers_raw
[
i
]
+=
centers_incr
[
i
];
colvar_centers
[
i
]
=
colvar_centers_raw
[
i
];
variables
(
i
)
->
wrap
(
colvar_centers
[
i
]);
colvar_centers
[
i
].
apply_constraints
();
}
stage
++
;
cvm
::
log
(
"Moving restraint
\"
"
+
this
->
name
+
"
\"
stage "
+
cvm
::
to_str
(
stage
)
+
" : setting centers to "
+
cvm
::
to_str
(
colvar_centers
)
+
" at step "
+
cvm
::
to_str
(
cvm
::
step_absolute
()));
}
}
else
if
((
cvm
::
step_relative
()
>
0
)
&&
(
cvm
::
step_absolute
()
<=
target_nsteps
))
{
// move the restraint centers in the direction of the targets
// (slow growth)
for
(
size_t
i
=
0
;
i
<
num_variables
();
i
++
)
{
colvar_centers_raw
[
i
]
+=
centers_incr
[
i
];
colvar_centers
[
i
]
=
colvar_centers_raw
[
i
];
variables
(
i
)
->
wrap
(
colvar_centers
[
i
]);
colvar_centers
[
i
].
apply_constraints
();
}
}
if
(
cvm
::
debug
())
{
cvm
::
log
(
"New centers for the restraint bias
\"
"
+
this
->
name
+
"
\"
: "
+
cvm
::
to_str
(
colvar_centers
)
+
".
\n
"
);
}
}
return
COLVARS_OK
;
}
int
colvarbias_restraint_centers_moving
::
update_acc_work
()
{
if
(
b_output_acc_work
)
{
if
((
cvm
::
step_relative
()
>
0
)
||
(
cvm
::
step_absolute
()
==
0
))
{
for
(
size_t
i
=
0
;
i
<
num_variables
();
i
++
)
{
// project forces on the calculated increments at this step
acc_work
+=
colvar_forces
[
i
]
*
centers_incr
[
i
];
}
}
}
return
COLVARS_OK
;
}
std
::
string
const
colvarbias_restraint_centers_moving
::
get_state_params
()
const
{
std
::
ostringstream
os
;
os
.
setf
(
std
::
ios
::
scientific
,
std
::
ios
::
floatfield
);
if
(
b_chg_centers
)
{
size_t
i
;
os
<<
"centers "
;
for
(
i
=
0
;
i
<
num_variables
();
i
++
)
{
os
<<
" "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
colvar_centers
[
i
];
}
os
<<
"
\n
"
;
os
<<
"centers_raw "
;
for
(
i
=
0
;
i
<
num_variables
();
i
++
)
{
os
<<
" "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
colvar_centers_raw
[
i
];
}
os
<<
"
\n
"
;
if
(
b_output_acc_work
)
{
os
<<
"accumulatedWork "
<<
std
::
setprecision
(
cvm
::
en_prec
)
<<
std
::
setw
(
cvm
::
en_width
)
<<
acc_work
<<
"
\n
"
;
}
}
return
colvarbias_restraint_moving
::
get_state_params
()
+
os
.
str
();
}
int
colvarbias_restraint_centers_moving
::
set_state_params
(
std
::
string
const
&
conf
)
{
colvarbias_restraint
::
set_state_params
(
conf
);
if
(
b_chg_centers
)
{
// cvm::log ("Reading the updated restraint centers from the restart.\n");
if
(
!
get_keyval
(
conf
,
"centers"
,
colvar_centers
))
cvm
::
error
(
"Error: restraint centers are missing from the restart.
\n
"
);
if
(
!
get_keyval
(
conf
,
"centers_raw"
,
colvar_centers_raw
))
cvm
::
error
(
"Error:
\"
raw
\"
restraint centers are missing from the restart.
\n
"
);
if
(
b_output_acc_work
)
{
if
(
!
get_keyval
(
conf
,
"accumulatedWork"
,
acc_work
))
cvm
::
error
(
"Error: accumulatedWork is missing from the restart.
\n
"
);
}
}
return
COLVARS_OK
;
}
std
::
ostream
&
colvarbias_restraint_centers_moving
::
write_traj_label
(
std
::
ostream
&
os
)
{
if
(
b_output_centers
)
{
for
(
size_t
i
=
0
;
i
<
num_variables
();
i
++
)
{
size_t
const
this_cv_width
=
(
variables
(
i
)
->
value
()).
output_width
(
cvm
::
cv_width
);
os
<<
" x0_"
<<
cvm
::
wrap_string
(
variables
(
i
)
->
name
,
this_cv_width
-
3
);
}
}
if
(
b_output_acc_work
)
{
os
<<
" W_"
<<
cvm
::
wrap_string
(
this
->
name
,
cvm
::
en_width
-
2
);
}
return
os
;
}
std
::
ostream
&
colvarbias_restraint_centers_moving
::
write_traj
(
std
::
ostream
&
os
)
{
if
(
b_output_centers
)
{
for
(
size_t
i
=
0
;
i
<
num_variables
();
i
++
)
{
os
<<
" "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
colvar_centers
[
i
];
}
}
if
(
b_output_acc_work
)
{
os
<<
" "
<<
std
::
setprecision
(
cvm
::
en_prec
)
<<
std
::
setw
(
cvm
::
en_width
)
<<
acc_work
;
}
return
os
;
}
colvarbias_restraint_k_moving
::
colvarbias_restraint_k_moving
(
char
const
*
key
)
:
colvarbias
(
key
),
colvarbias_restraint
(
key
),
colvarbias_restraint_k
(
key
),
colvarbias_restraint_moving
(
key
)
{
b_chg_force_k
=
false
;
target_equil_steps
=
0
;
target_force_k
=
0.0
;
starting_force_k
=
0.0
;
force_k_exp
=
1.0
;
restraint_FE
=
0.0
;
}
int
colvarbias_restraint_k_moving
::
init
(
std
::
string
const
&
conf
)
{
colvarbias_restraint_k
::
init
(
conf
);
if
(
get_keyval
(
conf
,
"targetForceConstant"
,
target_force_k
,
target_force_k
))
{
starting_force_k
=
force_k
;
b_chg_force_k
=
true
;
}
if
(
b_chg_force_k
)
{
// parse moving restraint options
colvarbias_restraint_moving
::
init
(
conf
);
}
else
{
return
COLVARS_OK
;
}
get_keyval
(
conf
,
"targetEquilSteps"
,
target_equil_steps
,
target_equil_steps
);
if
(
get_keyval
(
conf
,
"lambdaSchedule"
,
lambda_schedule
,
lambda_schedule
)
&&
target_nstages
>
0
)
{
cvm
::
error
(
"Error: targetNumStages and lambdaSchedule are incompatible.
\n
"
,
INPUT_ERROR
);
return
cvm
::
get_error
();
}
if
(
lambda_schedule
.
size
())
{
// There is one more lambda-point than stages
target_nstages
=
lambda_schedule
.
size
()
-
1
;
}
if
(
get_keyval
(
conf
,
"targetForceExponent"
,
force_k_exp
,
force_k_exp
))
{
if
(
!
b_chg_force_k
)
cvm
::
log
(
"Warning: not changing force constant: targetForceExponent will be ignored
\n
"
);
}
if
(
force_k_exp
<
1.0
)
{
cvm
::
log
(
"Warning: for all practical purposes, targetForceExponent should be 1.0 or greater.
\n
"
);
}
return
COLVARS_OK
;
}
int
colvarbias_restraint_k_moving
::
update
()
{
if
(
b_chg_force_k
)
{
cvm
::
real
lambda
;
if
(
target_nstages
)
{
if
(
cvm
::
step_absolute
()
==
0
)
{
// Setup first stage of staged variable force constant calculation
if
(
lambda_schedule
.
size
())
{
lambda
=
lambda_schedule
[
0
];
}
else
{
lambda
=
0.0
;
}
force_k
=
starting_force_k
+
(
target_force_k
-
starting_force_k
)
*
std
::
pow
(
lambda
,
force_k_exp
);
cvm
::
log
(
"Restraint "
+
this
->
name
+
", stage "
+
cvm
::
to_str
(
stage
)
+
" : lambda = "
+
cvm
::
to_str
(
lambda
)
+
", k = "
+
cvm
::
to_str
(
force_k
));
}
// TI calculation: estimate free energy derivative
// need current lambda
if
(
lambda_schedule
.
size
())
{
lambda
=
lambda_schedule
[
stage
];
}
else
{
lambda
=
cvm
::
real
(
stage
)
/
cvm
::
real
(
target_nstages
);
}
if
(
target_equil_steps
==
0
||
cvm
::
step_absolute
()
%
target_nsteps
>=
target_equil_steps
)
{
// Start averaging after equilibration period, if requested
// Square distance normalized by square colvar width
cvm
::
real
dist_sq
=
0.0
;
for
(
size_t
i
=
0
;
i
<
num_variables
();
i
++
)
{
dist_sq
+=
d_restraint_potential_dk
(
i
);
}
restraint_FE
+=
0.5
*
force_k_exp
*
std
::
pow
(
lambda
,
force_k_exp
-
1.0
)
*
(
target_force_k
-
starting_force_k
)
*
dist_sq
;
}
// Finish current stage...
if
(
cvm
::
step_absolute
()
%
target_nsteps
==
0
&&
cvm
::
step_absolute
()
>
0
)
{
cvm
::
log
(
"Restraint "
+
this
->
name
+
" Lambda= "
+
cvm
::
to_str
(
lambda
)
+
" dA/dLambda= "
+
cvm
::
to_str
(
restraint_FE
/
cvm
::
real
(
target_nsteps
-
target_equil_steps
)));
// ...and move on to the next one
if
(
stage
<
target_nstages
)
{
restraint_FE
=
0.0
;
stage
++
;
if
(
lambda_schedule
.
size
())
{
lambda
=
lambda_schedule
[
stage
];
}
else
{
lambda
=
cvm
::
real
(
stage
)
/
cvm
::
real
(
target_nstages
);
}
force_k
=
starting_force_k
+
(
target_force_k
-
starting_force_k
)
*
std
::
pow
(
lambda
,
force_k_exp
);
cvm
::
log
(
"Restraint "
+
this
->
name
+
", stage "
+
cvm
::
to_str
(
stage
)
+
" : lambda = "
+
cvm
::
to_str
(
lambda
)
+
", k = "
+
cvm
::
to_str
(
force_k
));
}
}
}
else
if
(
cvm
::
step_absolute
()
<=
target_nsteps
)
{
// update force constant (slow growth)
lambda
=
cvm
::
real
(
cvm
::
step_absolute
())
/
cvm
::
real
(
target_nsteps
);
force_k
=
starting_force_k
+
(
target_force_k
-
starting_force_k
)
*
std
::
pow
(
lambda
,
force_k_exp
);
}
}
return
COLVARS_OK
;
}
std
::
string
const
colvarbias_restraint_k_moving
::
get_state_params
()
const
{
std
::
ostringstream
os
;
os
.
setf
(
std
::
ios
::
scientific
,
std
::
ios
::
floatfield
);
if
(
b_chg_force_k
)
{
os
<<
"forceConstant "
<<
std
::
setprecision
(
cvm
::
en_prec
)
<<
std
::
setw
(
cvm
::
en_width
)
<<
force_k
<<
"
\n
"
;
}
return
colvarbias_restraint_moving
::
get_state_params
()
+
os
.
str
();
}
int
colvarbias_restraint_k_moving
::
set_state_params
(
std
::
string
const
&
conf
)
{
colvarbias_restraint
::
set_state_params
(
conf
);
if
(
b_chg_force_k
)
{
// cvm::log ("Reading the updated force constant from the restart.\n");
if
(
!
get_keyval
(
conf
,
"forceConstant"
,
force_k
,
force_k
))
cvm
::
error
(
"Error: force constant is missing from the restart.
\n
"
);
}
return
COLVARS_OK
;
}
std
::
ostream
&
colvarbias_restraint_k_moving
::
write_traj_label
(
std
::
ostream
&
os
)
{
return
os
;
}
std
::
ostream
&
colvarbias_restraint_k_moving
::
write_traj
(
std
::
ostream
&
os
)
{
return
os
;
}
// redefined due to legacy state file keyword "harmonic"
std
::
istream
&
colvarbias_restraint
::
read_state
(
std
::
istream
&
is
)
{
size_t
const
start_pos
=
is
.
tellg
();
std
::
string
key
,
brace
,
conf
;
if
(
!
(
is
>>
key
)
||
!
(
key
==
"restraint"
||
key
==
"harmonic"
)
||
!
(
is
>>
brace
)
||
!
(
brace
==
"{"
)
||
!
(
is
>>
colvarparse
::
read_block
(
"configuration"
,
conf
))
||
(
set_state_params
(
conf
)
!=
COLVARS_OK
)
)
{
cvm
::
error
(
"Error: in reading state configuration for
\"
"
+
bias_type
+
"
\"
bias
\"
"
+
this
->
name
+
"
\"
at position "
+
cvm
::
to_str
(
is
.
tellg
())
+
" in stream.
\n
"
,
INPUT_ERROR
);
is
.
clear
();
is
.
seekg
(
start_pos
,
std
::
ios
::
beg
);
is
.
setstate
(
std
::
ios
::
failbit
);
return
is
;
}
if
(
!
read_state_data
(
is
))
{
cvm
::
error
(
"Error: in reading state data for
\"
"
+
bias_type
+
"
\"
bias
\"
"
+
this
->
name
+
"
\"
at position "
+
cvm
::
to_str
(
is
.
tellg
())
+
" in stream.
\n
"
,
INPUT_ERROR
);
is
.
clear
();
is
.
seekg
(
start_pos
,
std
::
ios
::
beg
);
is
.
setstate
(
std
::
ios
::
failbit
);
}
is
>>
brace
;
if
(
brace
!=
"}"
)
{
cvm
::
log
(
"brace = "
+
brace
+
"
\n
"
);
cvm
::
error
(
"Error: corrupt restart information for
\"
"
+
bias_type
+
"
\"
bias
\"
"
+
this
->
name
+
"
\"
: no matching brace at position "
+
cvm
::
to_str
(
is
.
tellg
())
+
" in stream.
\n
"
);
is
.
setstate
(
std
::
ios
::
failbit
);
}
return
is
;
}
std
::
ostream
&
colvarbias_restraint
::
write_state
(
std
::
ostream
&
os
)
{
os
.
setf
(
std
::
ios
::
scientific
,
std
::
ios
::
floatfield
);
os
<<
"restraint {
\n
"
<<
" configuration {
\n
"
;
std
::
istringstream
is
(
get_state_params
());
std
::
string
line
;
while
(
std
::
getline
(
is
,
line
))
{
os
<<
" "
<<
line
<<
"
\n
"
;
}
os
<<
" }
\n
"
;
write_state_data
(
os
);
os
<<
"}
\n\n
"
;
return
os
;
}
colvarbias_restraint_harmonic
::
colvarbias_restraint_harmonic
(
char
const
*
key
)
:
colvarbias
(
key
),
colvarbias_restraint
(
key
),
colvarbias_restraint_centers
(
key
),
colvarbias_restraint_moving
(
key
),
colvarbias_restraint_k
(
key
),
colvarbias_restraint_centers_moving
(
key
),
colvarbias_restraint_k_moving
(
key
)
{
}
int
colvarbias_restraint_harmonic
::
init
(
std
::
string
const
&
conf
)
{
colvarbias_restraint
::
init
(
conf
);
colvarbias_restraint_moving
::
init
(
conf
);
colvarbias_restraint_centers_moving
::
init
(
conf
);
colvarbias_restraint_k_moving
::
init
(
conf
);
for
(
size_t
i
=
0
;
i
<
num_variables
();
i
++
)
{
if
(
variables
(
i
)
->
width
!=
1.0
)
cvm
::
log
(
"The force constant for colvar
\"
"
+
variables
(
i
)
->
name
+
"
\"
will be rescaled to "
+
cvm
::
to_str
(
force_k
/
(
variables
(
i
)
->
width
*
variables
(
i
)
->
width
))
+
" according to the specified width.
\n
"
);
}
return
COLVARS_OK
;
}
int
colvarbias_restraint_harmonic
::
update
()
{
// update parameters (centers or force constant)
colvarbias_restraint_centers_moving
::
update
();
colvarbias_restraint_k_moving
::
update
();
// update restraint energy and forces
colvarbias_restraint
::
update
();
// update accumulated work using the current forces
colvarbias_restraint_centers_moving
::
update_acc_work
();
return
COLVARS_OK
;
}
cvm
::
real
colvarbias_restraint_harmonic
::
restraint_potential
(
size_t
i
)
const
{
return
0.5
*
force_k
/
(
variables
(
i
)
->
width
*
variables
(
i
)
->
width
)
*
variables
(
i
)
->
dist2
(
variables
(
i
)
->
value
(),
colvar_centers
[
i
]);
}
colvarvalue
const
colvarbias_restraint_harmonic
::
restraint_force
(
size_t
i
)
const
{
return
-
0.5
*
force_k
/
(
variables
(
i
)
->
width
*
variables
(
i
)
->
width
)
*
variables
(
i
)
->
dist2_lgrad
(
variables
(
i
)
->
value
(),
colvar_centers
[
i
]);
}
cvm
::
real
colvarbias_restraint_harmonic
::
d_restraint_potential_dk
(
size_t
i
)
const
{
return
0.5
/
(
variables
(
i
)
->
width
*
variables
(
i
)
->
width
)
*
variables
(
i
)
->
dist2
(
variables
(
i
)
->
value
(),
colvar_centers
[
i
]);
}
std
::
string
const
colvarbias_restraint_harmonic
::
get_state_params
()
const
{
return
colvarbias_restraint
::
get_state_params
()
+
colvarbias_restraint_centers_moving
::
get_state_params
()
+
colvarbias_restraint_k_moving
::
get_state_params
();
}
int
colvarbias_restraint_harmonic
::
set_state_params
(
std
::
string
const
&
conf
)
{
int
error_code
=
COLVARS_OK
;
error_code
|=
colvarbias_restraint
::
set_state_params
(
conf
);
error_code
|=
colvarbias_restraint_centers_moving
::
set_state_params
(
conf
);
error_code
|=
colvarbias_restraint_k_moving
::
set_state_params
(
conf
);
return
error_code
;
}
std
::
ostream
&
colvarbias_restraint_harmonic
::
write_traj_label
(
std
::
ostream
&
os
)
{
colvarbias_restraint
::
write_traj_label
(
os
);
colvarbias_restraint_centers_moving
::
write_traj_label
(
os
);
colvarbias_restraint_k_moving
::
write_traj_label
(
os
);
return
os
;
}
std
::
ostream
&
colvarbias_restraint_harmonic
::
write_traj
(
std
::
ostream
&
os
)
{
colvarbias_restraint
::
write_traj
(
os
);
colvarbias_restraint_centers_moving
::
write_traj
(
os
);
colvarbias_restraint_k_moving
::
write_traj
(
os
);
return
os
;
}
int
colvarbias_restraint_harmonic
::
change_configuration
(
std
::
string
const
&
conf
)
{
return
colvarbias_restraint_centers
::
change_configuration
(
conf
)
|
colvarbias_restraint_k
::
change_configuration
(
conf
);
}
cvm
::
real
colvarbias_restraint_harmonic
::
energy_difference
(
std
::
string
const
&
conf
)
{
cvm
::
real
const
old_bias_energy
=
bias_energy
;
cvm
::
real
const
old_force_k
=
force_k
;
std
::
vector
<
colvarvalue
>
const
old_centers
=
colvar_centers
;
change_configuration
(
conf
);
update
();
cvm
::
real
const
result
=
(
bias_energy
-
old_bias_energy
);
bias_energy
=
old_bias_energy
;
force_k
=
old_force_k
;
colvar_centers
=
old_centers
;
return
result
;
}
colvarbias_restraint_harmonic_walls
::
colvarbias_restraint_harmonic_walls
(
char
const
*
key
)
:
colvarbias
(
key
),
colvarbias_restraint
(
key
),
colvarbias_restraint_k
(
key
),
colvarbias_restraint_moving
(
key
),
colvarbias_restraint_k_moving
(
key
)
{
lower_wall_k
=
0.0
;
upper_wall_k
=
0.0
;
}
int
colvarbias_restraint_harmonic_walls
::
init
(
std
::
string
const
&
conf
)
{
colvarbias_restraint
::
init
(
conf
);
colvarbias_restraint_moving
::
init
(
conf
);
colvarbias_restraint_k_moving
::
init
(
conf
);
get_keyval
(
conf
,
"lowerWallConstant"
,
lower_wall_k
,
(
lower_wall_k
>
0.0
)
?
lower_wall_k
:
force_k
);
get_keyval
(
conf
,
"upperWallConstant"
,
upper_wall_k
,
(
upper_wall_k
>
0.0
)
?
upper_wall_k
:
force_k
);
enable
(
f_cvb_scalar_variables
);
size_t
i
;
bool
b_null_lower_walls
=
false
;
if
(
lower_walls
.
size
()
==
0
)
{
b_null_lower_walls
=
true
;
lower_walls
.
resize
(
num_variables
());
for
(
i
=
0
;
i
<
num_variables
();
i
++
)
{
lower_walls
[
i
].
type
(
variables
(
i
)
->
value
());
lower_walls
[
i
].
reset
();
}
}
if
(
!
get_keyval
(
conf
,
"lowerWalls"
,
lower_walls
,
lower_walls
)
&&
b_null_lower_walls
)
{
cvm
::
log
(
"Lower walls were not provided.
\n
"
);
lower_walls
.
resize
(
0
);
}
bool
b_null_upper_walls
=
false
;
if
(
upper_walls
.
size
()
==
0
)
{
b_null_upper_walls
=
true
;
upper_walls
.
resize
(
num_variables
());
for
(
i
=
0
;
i
<
num_variables
();
i
++
)
{
upper_walls
[
i
].
type
(
variables
(
i
)
->
value
());
upper_walls
[
i
].
reset
();
}
}
if
(
!
get_keyval
(
conf
,
"upperWalls"
,
upper_walls
,
upper_walls
)
&&
b_null_upper_walls
)
{
cvm
::
log
(
"Upper walls were not provided.
\n
"
);
upper_walls
.
resize
(
0
);
}
if
((
lower_walls
.
size
()
==
0
)
&&
(
upper_walls
.
size
()
==
0
))
{
cvm
::
error
(
"Error: no walls provided.
\n
"
,
INPUT_ERROR
);
return
INPUT_ERROR
;
}
if
((
lower_walls
.
size
()
==
0
)
||
(
upper_walls
.
size
()
==
0
))
{
for
(
i
=
0
;
i
<
num_variables
();
i
++
)
{
if
(
variables
(
i
)
->
is_enabled
(
f_cv_periodic
))
{
cvm
::
error
(
"Error: at least one variable is periodic, "
"both walls must be provided.
\n
"
,
INPUT_ERROR
);
return
INPUT_ERROR
;
}
}
}
if
((
lower_walls
.
size
()
>
0
)
&&
(
upper_walls
.
size
()
>
0
))
{
for
(
i
=
0
;
i
<
num_variables
();
i
++
)
{
if
(
lower_walls
[
i
]
>=
upper_walls
[
i
])
{
cvm
::
error
(
"Error: one upper wall, "
+
cvm
::
to_str
(
upper_walls
[
i
])
+
", is not higher than the lower wall, "
+
cvm
::
to_str
(
lower_walls
[
i
])
+
".
\n
"
,
INPUT_ERROR
);
}
}
if
(
lower_wall_k
*
upper_wall_k
==
0.0
)
{
cvm
::
error
(
"Error: lowerWallConstant and upperWallConstant, "
"when defined, must both be positive.
\n
"
,
INPUT_ERROR
);
return
INPUT_ERROR
;
}
force_k
=
lower_wall_k
*
upper_wall_k
;
// transform the two constants to relative values
// (allow changing both via force_k)
lower_wall_k
/=
force_k
;
upper_wall_k
/=
force_k
;
}
for
(
i
=
0
;
i
<
num_variables
();
i
++
)
{
if
(
variables
(
i
)
->
width
!=
1.0
)
cvm
::
log
(
"The force constant for colvar
\"
"
+
variables
(
i
)
->
name
+
"
\"
will be rescaled to "
+
cvm
::
to_str
(
force_k
/
(
variables
(
i
)
->
width
*
variables
(
i
)
->
width
))
+
" according to the specified width.
\n
"
);
}
return
COLVARS_OK
;
}
int
colvarbias_restraint_harmonic_walls
::
update
()
{
colvarbias_restraint_k_moving
::
update
();
colvarbias_restraint
::
update
();
return
COLVARS_OK
;
}
void
colvarbias_restraint_harmonic_walls
::
communicate_forces
()
{
for
(
size_t
i
=
0
;
i
<
num_variables
();
i
++
)
{
if
(
cvm
::
debug
())
{
cvm
::
log
(
"Communicating a force to colvar
\"
"
+
variables
(
i
)
->
name
+
"
\"
.
\n
"
);
}
variables
(
i
)
->
add_bias_force_actual_value
(
colvar_forces
[
i
]);
}
}
cvm
::
real
colvarbias_restraint_harmonic_walls
::
colvar_distance
(
size_t
i
)
const
{
colvar
*
cv
=
variables
(
i
);
colvarvalue
const
&
cvv
=
variables
(
i
)
->
actual_value
();
// For a periodic colvar, both walls may be applicable at the same time
// in which case we pick the closer one
if
(
cv
->
is_enabled
(
f_cv_periodic
))
{
cvm
::
real
const
lower_wall_dist2
=
cv
->
dist2
(
cvv
,
lower_walls
[
i
]);
cvm
::
real
const
upper_wall_dist2
=
cv
->
dist2
(
cvv
,
upper_walls
[
i
]);
if
(
lower_wall_dist2
<
upper_wall_dist2
)
{
cvm
::
real
const
grad
=
cv
->
dist2_lgrad
(
cvv
,
lower_walls
[
i
]);
if
(
grad
<
0.0
)
{
return
0.5
*
grad
;
}
}
else
{
cvm
::
real
const
grad
=
cv
->
dist2_lgrad
(
cvv
,
upper_walls
[
i
]);
if
(
grad
>
0.0
)
{
return
0.5
*
grad
;
}
}
return
0.0
;
}
if
(
lower_walls
.
size
()
>
0
)
{
cvm
::
real
const
grad
=
cv
->
dist2_lgrad
(
cvv
,
lower_walls
[
i
]);
if
(
grad
<
0.0
)
{
return
0.5
*
grad
;
}
}
if
(
upper_walls
.
size
()
>
0
)
{
cvm
::
real
const
grad
=
cv
->
dist2_lgrad
(
cvv
,
upper_walls
[
i
]);
if
(
grad
>
0.0
)
{
return
0.5
*
grad
;
}
}
return
0.0
;
}
cvm
::
real
colvarbias_restraint_harmonic_walls
::
restraint_potential
(
size_t
i
)
const
{
cvm
::
real
const
dist
=
colvar_distance
(
i
);
cvm
::
real
const
scale
=
dist
>
0.0
?
upper_wall_k
:
lower_wall_k
;
return
0.5
*
force_k
*
scale
/
(
variables
(
i
)
->
width
*
variables
(
i
)
->
width
)
*
dist
*
dist
;
}
colvarvalue
const
colvarbias_restraint_harmonic_walls
::
restraint_force
(
size_t
i
)
const
{
cvm
::
real
const
dist
=
colvar_distance
(
i
);
cvm
::
real
const
scale
=
dist
>
0.0
?
upper_wall_k
:
lower_wall_k
;
return
-
force_k
*
scale
/
(
variables
(
i
)
->
width
*
variables
(
i
)
->
width
)
*
dist
;
}
cvm
::
real
colvarbias_restraint_harmonic_walls
::
d_restraint_potential_dk
(
size_t
i
)
const
{
cvm
::
real
const
dist
=
colvar_distance
(
i
);
cvm
::
real
const
scale
=
dist
>
0.0
?
upper_wall_k
:
lower_wall_k
;
return
0.5
*
scale
/
(
variables
(
i
)
->
width
*
variables
(
i
)
->
width
)
*
dist
*
dist
;
}
std
::
string
const
colvarbias_restraint_harmonic_walls
::
get_state_params
()
const
{
return
colvarbias_restraint
::
get_state_params
()
+
colvarbias_restraint_k_moving
::
get_state_params
();
}
int
colvarbias_restraint_harmonic_walls
::
set_state_params
(
std
::
string
const
&
conf
)
{
int
error_code
=
COLVARS_OK
;
error_code
|=
colvarbias_restraint
::
set_state_params
(
conf
);
error_code
|=
colvarbias_restraint_k_moving
::
set_state_params
(
conf
);
return
error_code
;
}
std
::
ostream
&
colvarbias_restraint_harmonic_walls
::
write_traj_label
(
std
::
ostream
&
os
)
{
colvarbias_restraint
::
write_traj_label
(
os
);
colvarbias_restraint_k_moving
::
write_traj_label
(
os
);
return
os
;
}
std
::
ostream
&
colvarbias_restraint_harmonic_walls
::
write_traj
(
std
::
ostream
&
os
)
{
colvarbias_restraint
::
write_traj
(
os
);
colvarbias_restraint_k_moving
::
write_traj
(
os
);
return
os
;
}
colvarbias_restraint_linear
::
colvarbias_restraint_linear
(
char
const
*
key
)
:
colvarbias
(
key
),
colvarbias_restraint
(
key
),
colvarbias_restraint_centers
(
key
),
colvarbias_restraint_moving
(
key
),
colvarbias_restraint_k
(
key
),
colvarbias_restraint_centers_moving
(
key
),
colvarbias_restraint_k_moving
(
key
)
{
}
int
colvarbias_restraint_linear
::
init
(
std
::
string
const
&
conf
)
{
colvarbias_restraint
::
init
(
conf
);
colvarbias_restraint_moving
::
init
(
conf
);
colvarbias_restraint_centers_moving
::
init
(
conf
);
colvarbias_restraint_k_moving
::
init
(
conf
);
for
(
size_t
i
=
0
;
i
<
num_variables
();
i
++
)
{
if
(
variables
(
i
)
->
is_enabled
(
f_cv_periodic
))
{
cvm
::
error
(
"Error: linear biases cannot be applied to periodic variables.
\n
"
,
INPUT_ERROR
);
return
INPUT_ERROR
;
}
if
(
variables
(
i
)
->
width
!=
1.0
)
cvm
::
log
(
"The force constant for colvar
\"
"
+
variables
(
i
)
->
name
+
"
\"
will be rescaled to "
+
cvm
::
to_str
(
force_k
/
variables
(
i
)
->
width
)
+
" according to the specified width.
\n
"
);
}
return
COLVARS_OK
;
}
int
colvarbias_restraint_linear
::
update
()
{
// update parameters (centers or force constant)
colvarbias_restraint_centers_moving
::
update
();
colvarbias_restraint_k_moving
::
update
();
// update restraint energy and forces
colvarbias_restraint
::
update
();
// update accumulated work using the current forces
colvarbias_restraint_centers_moving
::
update_acc_work
();
return
COLVARS_OK
;
}
int
colvarbias_restraint_linear
::
change_configuration
(
std
::
string
const
&
conf
)
{
// Only makes sense to change the force constant
return
colvarbias_restraint_k
::
change_configuration
(
conf
);
}
cvm
::
real
colvarbias_restraint_linear
::
energy_difference
(
std
::
string
const
&
conf
)
{
cvm
::
real
const
old_bias_energy
=
bias_energy
;
cvm
::
real
const
old_force_k
=
force_k
;
change_configuration
(
conf
);
update
();
cvm
::
real
const
result
=
(
bias_energy
-
old_bias_energy
);
bias_energy
=
old_bias_energy
;
force_k
=
old_force_k
;
return
result
;
}
cvm
::
real
colvarbias_restraint_linear
::
restraint_potential
(
size_t
i
)
const
{
return
force_k
/
variables
(
i
)
->
width
*
(
variables
(
i
)
->
value
()
-
colvar_centers
[
i
]);
}
colvarvalue
const
colvarbias_restraint_linear
::
restraint_force
(
size_t
i
)
const
{
return
-
1.0
*
force_k
/
variables
(
i
)
->
width
;
}
cvm
::
real
colvarbias_restraint_linear
::
d_restraint_potential_dk
(
size_t
i
)
const
{
return
1.0
/
variables
(
i
)
->
width
*
(
variables
(
i
)
->
value
()
-
colvar_centers
[
i
]);
}
std
::
string
const
colvarbias_restraint_linear
::
get_state_params
()
const
{
return
colvarbias_restraint
::
get_state_params
()
+
colvarbias_restraint_centers_moving
::
get_state_params
()
+
colvarbias_restraint_k_moving
::
get_state_params
();
}
int
colvarbias_restraint_linear
::
set_state_params
(
std
::
string
const
&
conf
)
{
int
error_code
=
COLVARS_OK
;
error_code
|=
colvarbias_restraint
::
set_state_params
(
conf
);
error_code
|=
colvarbias_restraint_centers_moving
::
set_state_params
(
conf
);
error_code
|=
colvarbias_restraint_k_moving
::
set_state_params
(
conf
);
return
error_code
;
}
std
::
ostream
&
colvarbias_restraint_linear
::
write_traj_label
(
std
::
ostream
&
os
)
{
colvarbias_restraint
::
write_traj_label
(
os
);
colvarbias_restraint_centers_moving
::
write_traj_label
(
os
);
colvarbias_restraint_k_moving
::
write_traj_label
(
os
);
return
os
;
}
std
::
ostream
&
colvarbias_restraint_linear
::
write_traj
(
std
::
ostream
&
os
)
{
colvarbias_restraint
::
write_traj
(
os
);
colvarbias_restraint_centers_moving
::
write_traj
(
os
);
colvarbias_restraint_k_moving
::
write_traj
(
os
);
return
os
;
}
colvarbias_restraint_histogram
::
colvarbias_restraint_histogram
(
char
const
*
key
)
:
colvarbias
(
key
)
{
lower_boundary
=
0.0
;
upper_boundary
=
0.0
;
width
=
0.0
;
gaussian_width
=
0.0
;
}
int
colvarbias_restraint_histogram
::
init
(
std
::
string
const
&
conf
)
{
colvarbias
::
init
(
conf
);
get_keyval
(
conf
,
"lowerBoundary"
,
lower_boundary
,
lower_boundary
);
get_keyval
(
conf
,
"upperBoundary"
,
upper_boundary
,
upper_boundary
);
get_keyval
(
conf
,
"width"
,
width
,
width
);
if
(
width
<=
0.0
)
{
cvm
::
error
(
"Error:
\"
width
\"
must be positive.
\n
"
,
INPUT_ERROR
);
}
get_keyval
(
conf
,
"gaussianWidth"
,
gaussian_width
,
2.0
*
width
,
colvarparse
::
parse_silent
);
get_keyval
(
conf
,
"gaussianSigma"
,
gaussian_width
,
2.0
*
width
);
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
);
}
cvm
::
real
const
nbins
=
(
upper_boundary
-
lower_boundary
)
/
width
;
int
const
nbins_round
=
(
int
)(
nbins
);
if
(
std
::
fabs
(
nbins
-
cvm
::
real
(
nbins_round
))
>
1.0E-10
)
{
cvm
::
log
(
"Warning: grid interval ("
+
cvm
::
to_str
(
lower_boundary
,
cvm
::
cv_width
,
cvm
::
cv_prec
)
+
" - "
+
cvm
::
to_str
(
upper_boundary
,
cvm
::
cv_width
,
cvm
::
cv_prec
)
+
") is not commensurate to its bin width ("
+
cvm
::
to_str
(
width
,
cvm
::
cv_width
,
cvm
::
cv_prec
)
+
").
\n
"
);
}
p
.
resize
(
nbins_round
);
ref_p
.
resize
(
nbins_round
);
p_diff
.
resize
(
nbins_round
);
bool
const
inline_ref_p
=
get_keyval
(
conf
,
"refHistogram"
,
ref_p
.
data_array
(),
ref_p
.
data_array
());
std
::
string
ref_p_file
;
get_keyval
(
conf
,
"refHistogramFile"
,
ref_p_file
,
std
::
string
(
""
));
if
(
ref_p_file
.
size
())
{
if
(
inline_ref_p
)
{
cvm
::
error
(
"Error: cannot specify both refHistogram and refHistogramFile at the same time.
\n
"
,
INPUT_ERROR
);
}
else
{
std
::
ifstream
is
(
ref_p_file
.
c_str
());
std
::
string
data_s
=
""
;
std
::
string
line
;
while
(
getline_nocomments
(
is
,
line
))
{
data_s
.
append
(
line
+
"
\n
"
);
}
if
(
data_s
.
size
()
==
0
)
{
cvm
::
error
(
"Error: file
\"
"
+
ref_p_file
+
"
\"
empty or unreadable.
\n
"
,
FILE_ERROR
);
}
is
.
close
();
cvm
::
vector1d
<
cvm
::
real
>
data
;
if
(
data
.
from_simple_string
(
data_s
)
!=
0
)
{
cvm
::
error
(
"Error: could not read histogram from file
\"
"
+
ref_p_file
+
"
\"
.
\n
"
);
}
if
(
data
.
size
()
==
2
*
ref_p
.
size
())
{
// file contains both x and p(x)
size_t
i
;
for
(
i
=
0
;
i
<
ref_p
.
size
();
i
++
)
{
ref_p
[
i
]
=
data
[
2
*
i
+
1
];
}
}
else
if
(
data
.
size
()
==
ref_p
.
size
())
{
ref_p
=
data
;
}
else
{
cvm
::
error
(
"Error: file
\"
"
+
ref_p_file
+
"
\"
contains a histogram of different length.
\n
"
,
INPUT_ERROR
);
}
}
}
cvm
::
real
const
ref_integral
=
ref_p
.
sum
()
*
width
;
if
(
std
::
fabs
(
ref_integral
-
1.0
)
>
1.0e-03
)
{
cvm
::
log
(
"Reference distribution not normalized, normalizing to unity.
\n
"
);
ref_p
/=
ref_integral
;
}
get_keyval
(
conf
,
"writeHistogram"
,
b_write_histogram
,
false
);
get_keyval
(
conf
,
"forceConstant"
,
force_k
,
1.0
);
return
COLVARS_OK
;
}
colvarbias_restraint_histogram
::~
colvarbias_restraint_histogram
()
{
p
.
resize
(
0
);
ref_p
.
resize
(
0
);
p_diff
.
resize
(
0
);
}
int
colvarbias_restraint_histogram
::
update
()
{
if
(
cvm
::
debug
())
cvm
::
log
(
"Updating the histogram restraint bias
\"
"
+
this
->
name
+
"
\"
.
\n
"
);
size_t
vector_size
=
0
;
size_t
icv
;
for
(
icv
=
0
;
icv
<
num_variables
();
icv
++
)
{
vector_size
+=
variables
(
icv
)
->
value
().
size
();
}
cvm
::
real
const
norm
=
1.0
/
(
std
::
sqrt
(
2.0
*
PI
)
*
gaussian_width
*
vector_size
);
// calculate the histogram
p
.
reset
();
for
(
icv
=
0
;
icv
<
num_variables
();
icv
++
)
{
colvarvalue
const
&
cv
=
variables
(
icv
)
->
value
();
if
(
cv
.
type
()
==
colvarvalue
::
type_scalar
)
{
cvm
::
real
const
cv_value
=
cv
.
real_value
;
size_t
igrid
;
for
(
igrid
=
0
;
igrid
<
p
.
size
();
igrid
++
)
{
cvm
::
real
const
x_grid
=
(
lower_boundary
+
(
igrid
+
0.5
)
*
width
);
p
[
igrid
]
+=
norm
*
std
::
exp
(
-
1.0
*
(
x_grid
-
cv_value
)
*
(
x_grid
-
cv_value
)
/
(
2.0
*
gaussian_width
*
gaussian_width
));
}
}
else
if
(
cv
.
type
()
==
colvarvalue
::
type_vector
)
{
size_t
idim
;
for
(
idim
=
0
;
idim
<
cv
.
vector1d_value
.
size
();
idim
++
)
{
cvm
::
real
const
cv_value
=
cv
.
vector1d_value
[
idim
];
size_t
igrid
;
for
(
igrid
=
0
;
igrid
<
p
.
size
();
igrid
++
)
{
cvm
::
real
const
x_grid
=
(
lower_boundary
+
(
igrid
+
0.5
)
*
width
);
p
[
igrid
]
+=
norm
*
std
::
exp
(
-
1.0
*
(
x_grid
-
cv_value
)
*
(
x_grid
-
cv_value
)
/
(
2.0
*
gaussian_width
*
gaussian_width
));
}
}
}
else
{
cvm
::
error
(
"Error: unsupported type for variable "
+
variables
(
icv
)
->
name
+
".
\n
"
,
COLVARS_NOT_IMPLEMENTED
);
return
COLVARS_NOT_IMPLEMENTED
;
}
}
cvm
::
real
const
force_k_cv
=
force_k
*
vector_size
;
// calculate the difference between current and reference
p_diff
=
p
-
ref_p
;
bias_energy
=
0.5
*
force_k_cv
*
p_diff
*
p_diff
;
// calculate the forces
for
(
icv
=
0
;
icv
<
num_variables
();
icv
++
)
{
colvarvalue
const
&
cv
=
variables
(
icv
)
->
value
();
colvarvalue
&
cv_force
=
colvar_forces
[
icv
];
cv_force
.
type
(
cv
);
cv_force
.
reset
();
if
(
cv
.
type
()
==
colvarvalue
::
type_scalar
)
{
cvm
::
real
const
cv_value
=
cv
.
real_value
;
cvm
::
real
&
force
=
cv_force
.
real_value
;
size_t
igrid
;
for
(
igrid
=
0
;
igrid
<
p
.
size
();
igrid
++
)
{
cvm
::
real
const
x_grid
=
(
lower_boundary
+
(
igrid
+
0.5
)
*
width
);
force
+=
force_k_cv
*
p_diff
[
igrid
]
*
norm
*
std
::
exp
(
-
1.0
*
(
x_grid
-
cv_value
)
*
(
x_grid
-
cv_value
)
/
(
2.0
*
gaussian_width
*
gaussian_width
))
*
(
-
1.0
*
(
x_grid
-
cv_value
)
/
(
gaussian_width
*
gaussian_width
));
}
}
else
if
(
cv
.
type
()
==
colvarvalue
::
type_vector
)
{
size_t
idim
;
for
(
idim
=
0
;
idim
<
cv
.
vector1d_value
.
size
();
idim
++
)
{
cvm
::
real
const
cv_value
=
cv
.
vector1d_value
[
idim
];
cvm
::
real
&
force
=
cv_force
.
vector1d_value
[
idim
];
size_t
igrid
;
for
(
igrid
=
0
;
igrid
<
p
.
size
();
igrid
++
)
{
cvm
::
real
const
x_grid
=
(
lower_boundary
+
(
igrid
+
0.5
)
*
width
);
force
+=
force_k_cv
*
p_diff
[
igrid
]
*
norm
*
std
::
exp
(
-
1.0
*
(
x_grid
-
cv_value
)
*
(
x_grid
-
cv_value
)
/
(
2.0
*
gaussian_width
*
gaussian_width
))
*
(
-
1.0
*
(
x_grid
-
cv_value
)
/
(
gaussian_width
*
gaussian_width
));
}
}
}
else
{
// TODO
}
}
return
COLVARS_OK
;
}
std
::
ostream
&
colvarbias_restraint_histogram
::
write_restart
(
std
::
ostream
&
os
)
{
if
(
b_write_histogram
)
{
std
::
string
file_name
(
cvm
::
output_prefix
()
+
"."
+
this
->
name
+
".hist.dat"
);
std
::
ofstream
os
(
file_name
.
c_str
());
os
<<
"# "
<<
cvm
::
wrap_string
(
variables
(
0
)
->
name
,
cvm
::
cv_width
)
<<
" "
<<
"p("
<<
cvm
::
wrap_string
(
variables
(
0
)
->
name
,
cvm
::
cv_width
-
3
)
<<
")
\n
"
;
size_t
igrid
;
for
(
igrid
=
0
;
igrid
<
p
.
size
();
igrid
++
)
{
cvm
::
real
const
x_grid
=
(
lower_boundary
+
(
igrid
+
1
)
*
width
);
os
<<
" "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
x_grid
<<
" "
<<
std
::
setprecision
(
cvm
::
cv_prec
)
<<
std
::
setw
(
cvm
::
cv_width
)
<<
p
[
igrid
]
<<
"
\n
"
;
}
os
.
close
();
}
return
os
;
}
std
::
istream
&
colvarbias_restraint_histogram
::
read_restart
(
std
::
istream
&
is
)
{
return
is
;
}
std
::
ostream
&
colvarbias_restraint_histogram
::
write_traj_label
(
std
::
ostream
&
os
)
{
os
<<
" "
;
if
(
b_output_energy
)
{
os
<<
" E_"
<<
cvm
::
wrap_string
(
this
->
name
,
cvm
::
en_width
-
2
);
}
return
os
;
}
std
::
ostream
&
colvarbias_restraint_histogram
::
write_traj
(
std
::
ostream
&
os
)
{
os
<<
" "
;
if
(
b_output_energy
)
{
os
<<
" "
<<
std
::
setprecision
(
cvm
::
en_prec
)
<<
std
::
setw
(
cvm
::
en_width
)
<<
bias_energy
;
}
return
os
;
}
Event Timeline
Log In to Comment