Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F82555043
colvaratoms.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
Thu, Sep 12, 02:22
Size
25 KB
Mime Type
text/x-c
Expires
Sat, Sep 14, 02:22 (2 d)
Engine
blob
Format
Raw Data
Handle
20600885
Attached To
rLAMMPS lammps
colvaratoms.cpp
View Options
/// -*- c++ -*-
#include "colvarmodule.h"
#include "colvarparse.h"
#include "colvaratoms.h"
// member functions of the "atom" class depend tightly on the MD interface, and are
// thus defined in colvarproxy_xxx.cpp
// in this file only atom_group functions are defined
// Note: "conf" is the configuration of the cvc who is using this atom group;
// "key" is the name of the atom group (e.g. "atoms", "group1", "group2", ...)
cvm
::
atom_group
::
atom_group
(
std
::
string
const
&
conf
,
char
const
*
key
)
:
b_center
(
false
),
b_rotate
(
false
),
b_user_defined_fit
(
false
),
b_fit_gradients
(
false
),
ref_pos_group
(
NULL
),
noforce
(
false
)
{
cvm
::
log
(
"Defining atom group
\"
"
+
std
::
string
(
key
)
+
"
\"
.
\n
"
);
// real work is done by parse
parse
(
conf
,
key
);
}
cvm
::
atom_group
::
atom_group
(
std
::
vector
<
cvm
::
atom
>
const
&
atoms
)
:
b_dummy
(
false
),
b_center
(
false
),
b_rotate
(
false
),
b_fit_gradients
(
false
),
ref_pos_group
(
NULL
),
noforce
(
false
)
{
this
->
reserve
(
atoms
.
size
());
for
(
size_t
i
=
0
;
i
<
atoms
.
size
();
i
++
)
{
this
->
push_back
(
atoms
[
i
]);
}
total_mass
=
0.0
;
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
total_mass
+=
ai
->
mass
;
}
}
cvm
::
atom_group
::
atom_group
()
:
b_dummy
(
false
),
b_center
(
false
),
b_rotate
(
false
),
b_user_defined_fit
(
false
),
b_fit_gradients
(
false
),
ref_pos_group
(
NULL
),
noforce
(
false
)
{
total_mass
=
0.0
;
}
cvm
::
atom_group
::~
atom_group
()
{
if
(
ref_pos_group
)
{
delete
ref_pos_group
;
ref_pos_group
=
NULL
;
}
}
void
cvm
::
atom_group
::
add_atom
(
cvm
::
atom
const
&
a
)
{
if
(
b_dummy
)
{
cvm
::
error
(
"Error: cannot add atoms to a dummy group.
\n
"
,
INPUT_ERROR
);
}
else
{
this
->
push_back
(
a
);
total_mass
+=
a
.
mass
;
}
}
void
cvm
::
atom_group
::
reset_mass
(
std
::
string
&
name
,
int
i
,
int
j
)
{
total_mass
=
0.0
;
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
total_mass
+=
ai
->
mass
;
}
cvm
::
log
(
"Re-initialized atom group "
+
name
+
":"
+
cvm
::
to_str
(
i
)
+
"/"
+
cvm
::
to_str
(
j
)
+
". "
+
cvm
::
to_str
(
this
->
size
())
+
" atoms: total mass = "
+
cvm
::
to_str
(
this
->
total_mass
)
+
".
\n
"
);
}
int
cvm
::
atom_group
::
parse
(
std
::
string
const
&
conf
,
char
const
*
key
)
{
std
::
string
group_conf
;
// save_delimiters is set to false for this call, because "conf" is
// not the config string of this group, but of its parent object
// (which has already taken care of the delimiters)
save_delimiters
=
false
;
key_lookup
(
conf
,
key
,
group_conf
,
dummy_pos
);
// restoring the normal value, because we do want keywords checked
// inside "group_conf"
save_delimiters
=
true
;
if
(
group_conf
.
size
()
==
0
)
{
cvm
::
error
(
"Error: atom group
\"
"
+
std
::
string
(
key
)
+
"
\"
is set, but has no definition.
\n
"
,
INPUT_ERROR
);
return
COLVARS_ERROR
;
}
cvm
::
increase_depth
();
cvm
::
log
(
"Initializing atom group
\"
"
+
std
::
string
(
key
)
+
"
\"
.
\n
"
);
// whether or not to include messages in the log
// colvarparse::Parse_Mode mode = parse_silent;
// {
// bool b_verbose;
// get_keyval (group_conf, "verboseOutput", b_verbose, false, parse_silent);
// if (b_verbose) mode = parse_normal;
// }
colvarparse
::
Parse_Mode
mode
=
parse_normal
;
{
// std::vector<int> atom_indexes;
std
::
string
numbers_conf
=
""
;
size_t
pos
=
0
;
std
::
vector
<
int
>
atom_indexes
;
while
(
key_lookup
(
group_conf
,
"atomNumbers"
,
numbers_conf
,
pos
))
{
if
(
numbers_conf
.
size
())
{
std
::
istringstream
is
(
numbers_conf
);
int
ia
;
while
(
is
>>
ia
)
{
atom_indexes
.
push_back
(
ia
);
}
}
if
(
atom_indexes
.
size
())
{
this
->
reserve
(
this
->
size
()
+
atom_indexes
.
size
());
for
(
size_t
i
=
0
;
i
<
atom_indexes
.
size
();
i
++
)
{
this
->
push_back
(
cvm
::
atom
(
atom_indexes
[
i
]));
}
if
(
cvm
::
get_error
())
return
COLVARS_ERROR
;
}
else
{
cvm
::
error
(
"Error: no numbers provided for
\"
"
"atomNumbers
\"
.
\n
"
,
INPUT_ERROR
);
return
COLVARS_ERROR
;
}
atom_indexes
.
clear
();
}
std
::
string
index_group_name
;
if
(
get_keyval
(
group_conf
,
"indexGroup"
,
index_group_name
))
{
// use an index group from the index file read globally
std
::
list
<
std
::
string
>::
iterator
names_i
=
cvm
::
index_group_names
.
begin
();
std
::
list
<
std
::
vector
<
int
>
>::
iterator
index_groups_i
=
cvm
::
index_groups
.
begin
();
for
(
;
names_i
!=
cvm
::
index_group_names
.
end
()
;
++
names_i
,
++
index_groups_i
)
{
if
(
*
names_i
==
index_group_name
)
break
;
}
if
(
names_i
==
cvm
::
index_group_names
.
end
())
{
cvm
::
error
(
"Error: could not find index group "
+
index_group_name
+
" among those provided by the index file.
\n
"
,
INPUT_ERROR
);
return
COLVARS_ERROR
;
}
this
->
reserve
(
index_groups_i
->
size
());
for
(
size_t
i
=
0
;
i
<
index_groups_i
->
size
();
i
++
)
{
this
->
push_back
(
cvm
::
atom
((
*
index_groups_i
)[
i
]));
}
if
(
cvm
::
get_error
())
return
COLVARS_ERROR
;
}
}
{
std
::
string
range_conf
=
""
;
size_t
pos
=
0
;
while
(
key_lookup
(
group_conf
,
"atomNumbersRange"
,
range_conf
,
pos
))
{
if
(
range_conf
.
size
())
{
std
::
istringstream
is
(
range_conf
);
int
initial
,
final
;
char
dash
;
if
(
(
is
>>
initial
)
&&
(
initial
>
0
)
&&
(
is
>>
dash
)
&&
(
dash
==
'-'
)
&&
(
is
>>
final
)
&&
(
final
>
0
)
)
{
for
(
int
anum
=
initial
;
anum
<=
final
;
anum
++
)
{
this
->
push_back
(
cvm
::
atom
(
anum
));
}
if
(
cvm
::
get_error
())
return
COLVARS_ERROR
;
range_conf
=
""
;
continue
;
}
}
cvm
::
error
(
"Error: no valid definition for
\"
atomNumbersRange
\"
,
\"
"
+
range_conf
+
"
\"
.
\n
"
,
INPUT_ERROR
);
}
}
std
::
vector
<
std
::
string
>
psf_segids
;
get_keyval
(
group_conf
,
"psfSegID"
,
psf_segids
,
std
::
vector
<
std
::
string
>
(),
mode
);
for
(
std
::
vector
<
std
::
string
>::
iterator
psii
=
psf_segids
.
begin
();
psii
<
psf_segids
.
end
();
++
psii
)
{
if
(
(
psii
->
size
()
==
0
)
||
(
psii
->
size
()
>
4
)
)
{
cvm
::
error
(
"Error: invalid segmend identifier provided,
\"
"
+
(
*
psii
)
+
"
\"
.
\n
"
,
INPUT_ERROR
);
}
}
{
std
::
string
range_conf
=
""
;
size_t
pos
=
0
;
size_t
range_count
=
0
;
std
::
vector
<
std
::
string
>::
iterator
psii
=
psf_segids
.
begin
();
while
(
key_lookup
(
group_conf
,
"atomNameResidueRange"
,
range_conf
,
pos
))
{
range_count
++
;
if
(
range_count
>
psf_segids
.
size
())
{
cvm
::
error
(
"Error: more instances of
\"
atomNameResidueRange
\"
than "
"values of
\"
psfSegID
\"
.
\n
"
,
INPUT_ERROR
);
}
std
::
string
const
&
psf_segid
=
psf_segids
.
size
()
?
*
psii
:
std
::
string
(
""
);
if
(
range_conf
.
size
())
{
std
::
istringstream
is
(
range_conf
);
std
::
string
atom_name
;
int
initial
,
final
;
char
dash
;
if
(
(
is
>>
atom_name
)
&&
(
atom_name
.
size
())
&&
(
is
>>
initial
)
&&
(
initial
>
0
)
&&
(
is
>>
dash
)
&&
(
dash
==
'-'
)
&&
(
is
>>
final
)
&&
(
final
>
0
)
)
{
for
(
int
resid
=
initial
;
resid
<=
final
;
resid
++
)
{
this
->
push_back
(
cvm
::
atom
(
resid
,
atom_name
,
psf_segid
));
}
if
(
cvm
::
get_error
())
return
COLVARS_ERROR
;
range_conf
=
""
;
}
else
{
cvm
::
error
(
"Error: cannot parse definition for
\"
"
"atomNameResidueRange
\"
,
\"
"
+
range_conf
+
"
\"
.
\n
"
);
}
}
else
{
cvm
::
error
(
"Error: atomNameResidueRange with empty definition.
\n
"
);
}
if
(
psf_segid
.
size
())
++
psii
;
}
}
{
// read the atoms from a file
std
::
string
atoms_file_name
;
if
(
get_keyval
(
group_conf
,
"atomsFile"
,
atoms_file_name
,
std
::
string
(
""
),
mode
))
{
std
::
string
atoms_col
;
if
(
!
get_keyval
(
group_conf
,
"atomsCol"
,
atoms_col
,
std
::
string
(
""
),
mode
))
{
cvm
::
error
(
"Error: parameter atomsCol is required if atomsFile is set.
\n
"
,
INPUT_ERROR
);
}
double
atoms_col_value
;
bool
const
atoms_col_value_defined
=
get_keyval
(
group_conf
,
"atomsColValue"
,
atoms_col_value
,
0.0
,
mode
);
if
(
atoms_col_value_defined
&&
(
!
atoms_col_value
))
{
cvm
::
error
(
"Error: atomsColValue, if provided, must be non-zero.
\n
"
,
INPUT_ERROR
);
}
cvm
::
load_atoms
(
atoms_file_name
.
c_str
(),
*
this
,
atoms_col
,
atoms_col_value
);
}
}
// Catch any errors from all the initialization steps above
if
(
cvm
::
get_error
())
return
COLVARS_ERROR
;
for
(
std
::
vector
<
cvm
::
atom
>::
iterator
a1
=
this
->
begin
();
a1
!=
this
->
end
();
++
a1
)
{
std
::
vector
<
cvm
::
atom
>::
iterator
a2
=
a1
;
++
a2
;
for
(
;
a2
!=
this
->
end
();
++
a2
)
{
if
(
a1
->
id
==
a2
->
id
)
{
if
(
cvm
::
debug
())
cvm
::
log
(
"Discarding doubly counted atom with number "
+
cvm
::
to_str
(
a1
->
id
+
1
)
+
".
\n
"
);
a2
=
this
->
erase
(
a2
);
if
(
a2
==
this
->
end
())
break
;
}
}
}
if
(
get_keyval
(
group_conf
,
"dummyAtom"
,
dummy_atom_pos
,
cvm
::
atom_pos
(),
mode
))
{
b_dummy
=
true
;
this
->
total_mass
=
1.0
;
}
else
b_dummy
=
false
;
if
(
b_dummy
&&
(
this
->
size
()))
{
cvm
::
error
(
"Error: cannot set up group
\"
"
+
std
::
string
(
key
)
+
"
\"
as a dummy atom "
"and provide it with atom definitions.
\n
"
,
INPUT_ERROR
);
}
#if(! defined(COLVARS_STANDALONE))
if
(
(
!
b_dummy
)
&&
(
!
cvm
::
b_analysis
)
&&
(
!
(
this
->
size
()))
)
{
cvm
::
error
(
"Error: no atoms defined for atom group
\"
"
+
std
::
string
(
key
)
+
"
\"
.
\n
"
);
}
#endif
if
(
!
b_dummy
)
{
// calculate total mass (TODO: this is the step that most needs deferred re-initialization)
this
->
total_mass
=
0.0
;
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
this
->
total_mass
+=
ai
->
mass
;
}
// whether these atoms will ever receive forces or not
bool
enable_forces
=
true
;
// disableForces is deprecated
if
(
get_keyval
(
group_conf
,
"enableForces"
,
enable_forces
,
true
,
mode
))
{
noforce
=
!
enable_forces
;
}
else
{
get_keyval
(
group_conf
,
"disableForces"
,
noforce
,
false
,
colvarparse
::
parse_silent
);
}
}
// FITTING OPTIONS
bool
b_defined_center
=
get_keyval
(
group_conf
,
"centerReference"
,
b_center
,
false
,
mode
);
bool
b_defined_rotate
=
get_keyval
(
group_conf
,
"rotateReference"
,
b_rotate
,
false
,
mode
);
// is the user setting explicit options?
b_user_defined_fit
=
b_defined_center
||
b_defined_rotate
;
get_keyval
(
group_conf
,
"enableFitGradients"
,
b_fit_gradients
,
true
,
mode
);
if
(
b_center
||
b_rotate
)
{
if
(
b_dummy
)
cvm
::
error
(
"Error: centerReference or rotateReference "
"cannot be defined for a dummy atom.
\n
"
);
if
(
key_lookup
(
group_conf
,
"refPositionsGroup"
))
{
// instead of this group, define another group to compute the fit
if
(
ref_pos_group
)
{
cvm
::
error
(
"Error: the atom group
\"
"
+
std
::
string
(
key
)
+
"
\"
has already a reference group "
"for the rototranslational fit, which was communicated by the "
"colvar component. You should not use refPositionsGroup "
"in this case.
\n
"
);
}
cvm
::
log
(
"Within atom group
\"
"
+
std
::
string
(
key
)
+
"
\"
:
\n
"
);
ref_pos_group
=
new
atom_group
(
group_conf
,
"refPositionsGroup"
);
// regardless of the configuration, fit gradients must be calculated by refPositionsGroup
ref_pos_group
->
b_fit_gradients
=
this
->
b_fit_gradients
;
this
->
b_fit_gradients
=
false
;
}
atom_group
*
group_for_fit
=
ref_pos_group
?
ref_pos_group
:
this
;
get_keyval
(
group_conf
,
"refPositions"
,
ref_pos
,
ref_pos
,
mode
);
std
::
string
ref_pos_file
;
if
(
get_keyval
(
group_conf
,
"refPositionsFile"
,
ref_pos_file
,
std
::
string
(
""
),
mode
))
{
if
(
ref_pos
.
size
())
{
cvm
::
error
(
"Error: cannot specify
\"
refPositionsFile
\"
and "
"
\"
refPositions
\"
at the same time.
\n
"
);
}
std
::
string
ref_pos_col
;
double
ref_pos_col_value
=
0.0
;
if
(
get_keyval
(
group_conf
,
"refPositionsCol"
,
ref_pos_col
,
std
::
string
(
""
),
mode
))
{
// if provided, use PDB column to select coordinates
bool
found
=
get_keyval
(
group_conf
,
"refPositionsColValue"
,
ref_pos_col_value
,
0.0
,
mode
);
if
(
found
&&
ref_pos_col_value
==
0.0
)
cvm
::
error
(
"Error: refPositionsColValue, "
"if provided, must be non-zero.
\n
"
);
}
else
{
// if not, rely on existing atom indices for the group
group_for_fit
->
create_sorted_ids
();
ref_pos
.
resize
(
group_for_fit
->
size
());
}
cvm
::
load_coords
(
ref_pos_file
.
c_str
(),
ref_pos
,
group_for_fit
->
sorted_ids
,
ref_pos_col
,
ref_pos_col_value
);
}
if
(
ref_pos
.
size
())
{
if
(
b_rotate
)
{
if
(
ref_pos
.
size
()
!=
group_for_fit
->
size
())
cvm
::
error
(
"Error: the number of reference positions provided("
+
cvm
::
to_str
(
ref_pos
.
size
())
+
") does not match the number of atoms within
\"
"
+
std
::
string
(
key
)
+
"
\"
("
+
cvm
::
to_str
(
group_for_fit
->
size
())
+
"): to perform a rotational fit, "
+
"these numbers should be equal.
\n
"
,
INPUT_ERROR
);
}
// save the center of geometry of ref_pos and subtract it
center_ref_pos
();
}
else
{
#if(! defined(COLVARS_STANDALONE))
cvm
::
error
(
"Error: no reference positions provided.
\n
"
);
#endif
}
if
(
b_fit_gradients
)
{
group_for_fit
->
fit_gradients
.
assign
(
group_for_fit
->
size
(),
cvm
::
atom_pos
(
0.0
,
0.0
,
0.0
));
rot
.
request_group1_gradients
(
group_for_fit
->
size
());
}
if
(
b_rotate
&&
!
noforce
)
{
cvm
::
log
(
"Warning: atom group
\"
"
+
std
::
string
(
key
)
+
"
\"
will be aligned to a fixed orientation given by the reference positions provided. "
"If the internal structure of the group changes too much (i.e. its RMSD is comparable "
"to its radius of gyration), the optimal rotation and its gradients may become discontinuous. "
"If that happens, use refPositionsGroup (or a different definition for it if already defined) "
"to align the coordinates.
\n
"
);
// initialize rot member data
rot
.
request_group1_gradients
(
this
->
size
());
}
}
if
(
cvm
::
debug
())
cvm
::
log
(
"Done initializing atom group with name
\"
"
+
std
::
string
(
key
)
+
"
\"
.
\n
"
);
this
->
check_keywords
(
group_conf
,
key
);
if
(
cvm
::
get_error
())
{
cvm
::
error
(
"Error setting up atom group
\"
"
+
std
::
string
(
key
)
+
"
\"
."
);
return
COLVARS_ERROR
;
}
cvm
::
log
(
"Atom group
\"
"
+
std
::
string
(
key
)
+
"
\"
defined, "
+
cvm
::
to_str
(
this
->
size
())
+
" atoms initialized: total mass = "
+
cvm
::
to_str
(
this
->
total_mass
)
+
".
\n
"
);
cvm
::
decrease_depth
();
return
(
cvm
::
get_error
()
?
COLVARS_ERROR
:
COLVARS_OK
);
}
int
cvm
::
atom_group
::
create_sorted_ids
(
void
)
{
// Only do the work if the vector is not yet populated
if
(
sorted_ids
.
size
())
return
COLVARS_OK
;
std
::
list
<
int
>
temp_id_list
;
cvm
::
atom_iter
ai
;
for
(
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
temp_id_list
.
push_back
(
ai
->
id
);
}
temp_id_list
.
sort
();
temp_id_list
.
unique
();
if
(
temp_id_list
.
size
()
!=
this
->
size
())
{
cvm
::
error
(
"Error: duplicate atom IDs in atom group? (found "
+
cvm
::
to_str
(
temp_id_list
.
size
())
+
" unique atom IDs instead of"
+
cvm
::
to_str
(
this
->
size
())
+
").
\n
"
);
return
COLVARS_ERROR
;
}
sorted_ids
=
std
::
vector
<
int
>
(
temp_id_list
.
size
());
unsigned
int
id_i
=
0
;
std
::
list
<
int
>::
iterator
li
;
for
(
li
=
temp_id_list
.
begin
();
li
!=
temp_id_list
.
end
();
++
li
)
{
sorted_ids
[
id_i
]
=
*
li
;
id_i
++
;
}
return
(
cvm
::
get_error
()
?
COLVARS_ERROR
:
COLVARS_OK
);
}
void
cvm
::
atom_group
::
center_ref_pos
()
{
// save the center of geometry of ref_pos and then subtract it from
// them; in this way it will be possible to use ref_pos also for
// the rotational fit
// This is called either by atom_group::parse or by CVCs that set
// reference positions (eg. RMSD, eigenvector)
ref_pos_cog
=
cvm
::
atom_pos
(
0.0
,
0.0
,
0.0
);
std
::
vector
<
cvm
::
atom_pos
>::
iterator
pi
;
for
(
pi
=
ref_pos
.
begin
();
pi
!=
ref_pos
.
end
();
++
pi
)
{
ref_pos_cog
+=
*
pi
;
}
ref_pos_cog
/=
(
cvm
::
real
)
ref_pos
.
size
();
for
(
pi
=
ref_pos
.
begin
();
pi
!=
ref_pos
.
end
();
++
pi
)
{
(
*
pi
)
-=
ref_pos_cog
;
}
}
void
cvm
::
atom_group
::
read_positions
()
{
if
(
b_dummy
)
return
;
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
ai
->
read_position
();
}
if
(
ref_pos_group
)
ref_pos_group
->
read_positions
();
}
void
cvm
::
atom_group
::
calc_apply_roto_translation
()
{
atom_group
*
fit_group
=
ref_pos_group
?
ref_pos_group
:
this
;
if
(
b_center
)
{
// center on the origin first
cvm
::
atom_pos
const
cog
=
fit_group
->
center_of_geometry
();
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
ai
->
pos
-=
cog
;
}
}
if
(
b_rotate
)
{
// rotate the group (around the center of geometry if b_center is
// true, around the origin otherwise)
rot
.
calc_optimal_rotation
(
fit_group
->
positions
(),
ref_pos
);
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
ai
->
pos
=
rot
.
rotate
(
ai
->
pos
);
}
}
if
(
b_center
)
{
// align with the center of geometry of ref_pos
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
ai
->
pos
+=
ref_pos_cog
;
}
}
}
void
cvm
::
atom_group
::
apply_translation
(
cvm
::
rvector
const
&
t
)
{
if
(
b_dummy
)
return
;
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
ai
->
pos
+=
t
;
}
}
void
cvm
::
atom_group
::
apply_rotation
(
cvm
::
rotation
const
&
rot
)
{
if
(
b_dummy
)
return
;
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
ai
->
pos
=
rot
.
rotate
(
ai
->
pos
);
}
}
void
cvm
::
atom_group
::
read_velocities
()
{
if
(
b_dummy
)
return
;
if
(
b_rotate
)
{
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
ai
->
read_velocity
();
ai
->
vel
=
rot
.
rotate
(
ai
->
vel
);
}
}
else
{
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
ai
->
read_velocity
();
}
}
}
void
cvm
::
atom_group
::
read_system_forces
()
{
if
(
b_dummy
)
return
;
if
(
b_rotate
)
{
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
ai
->
read_system_force
();
ai
->
system_force
=
rot
.
rotate
(
ai
->
system_force
);
}
}
else
{
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
ai
->
read_system_force
();
}
}
}
cvm
::
atom_pos
cvm
::
atom_group
::
center_of_geometry
()
const
{
if
(
b_dummy
)
return
dummy_atom_pos
;
cvm
::
atom_pos
cog
(
0.0
,
0.0
,
0.0
);
for
(
cvm
::
atom_const_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
cog
+=
ai
->
pos
;
}
cog
/=
this
->
size
();
return
cog
;
}
cvm
::
atom_pos
cvm
::
atom_group
::
center_of_mass
()
const
{
if
(
b_dummy
)
return
dummy_atom_pos
;
cvm
::
atom_pos
com
(
0.0
,
0.0
,
0.0
);
for
(
cvm
::
atom_const_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
com
+=
ai
->
mass
*
ai
->
pos
;
}
com
/=
this
->
total_mass
;
return
com
;
}
void
cvm
::
atom_group
::
set_weighted_gradient
(
cvm
::
rvector
const
&
grad
)
{
if
(
b_dummy
)
return
;
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
ai
->
grad
=
(
ai
->
mass
/
this
->
total_mass
)
*
grad
;
}
}
void
cvm
::
atom_group
::
calc_fit_gradients
()
{
if
(
b_dummy
)
return
;
if
((
!
b_center
)
&&
(
!
b_rotate
))
return
;
// no fit
if
(
cvm
::
debug
())
cvm
::
log
(
"Calculating fit gradients.
\n
"
);
atom_group
*
group_for_fit
=
ref_pos_group
?
ref_pos_group
:
this
;
group_for_fit
->
fit_gradients
.
assign
(
group_for_fit
->
size
(),
cvm
::
rvector
(
0.0
,
0.0
,
0.0
));
if
(
b_center
)
{
// add the center of geometry contribution to the gradients
for
(
size_t
i
=
0
;
i
<
this
->
size
();
i
++
)
{
// need to bring the gradients in original frame first
cvm
::
rvector
const
atom_grad
=
b_rotate
?
(
rot
.
inverse
()).
rotate
((
*
this
)[
i
].
grad
)
:
(
*
this
)[
i
].
grad
;
for
(
size_t
j
=
0
;
j
<
group_for_fit
->
size
();
j
++
)
{
group_for_fit
->
fit_gradients
[
j
]
+=
(
-
1.0
)
/
(
cvm
::
real
(
group_for_fit
->
size
()))
*
atom_grad
;
}
}
}
if
(
b_rotate
)
{
// add the rotation matrix contribution to the gradients
cvm
::
rotation
const
rot_inv
=
rot
.
inverse
();
cvm
::
atom_pos
const
cog
=
this
->
center_of_geometry
();
for
(
size_t
i
=
0
;
i
<
this
->
size
();
i
++
)
{
cvm
::
atom_pos
const
pos_orig
=
rot_inv
.
rotate
((
b_center
?
((
*
this
)[
i
].
pos
-
cog
)
:
((
*
this
)[
i
].
pos
)));
for
(
size_t
j
=
0
;
j
<
group_for_fit
->
size
();
j
++
)
{
// calculate \partial(R(q) \vec{x}_i)/\partial q) \cdot \partial\xi/\partial\vec{x}_i
cvm
::
quaternion
const
dxdq
=
rot
.
q
.
position_derivative_inner
(
pos_orig
,
(
*
this
)[
i
].
grad
);
// multiply by \cdot {\partial q}/\partial\vec{x}_j and add it to the fit gradients
for
(
size_t
iq
=
0
;
iq
<
4
;
iq
++
)
{
group_for_fit
->
fit_gradients
[
j
]
+=
dxdq
[
iq
]
*
rot
.
dQ0_1
[
j
][
iq
];
}
}
}
}
if
(
cvm
::
debug
())
cvm
::
log
(
"Done calculating fit gradients.
\n
"
);
}
std
::
vector
<
cvm
::
atom_pos
>
cvm
::
atom_group
::
positions
()
const
{
if
(
b_dummy
)
{
cvm
::
error
(
"Error: positions are not available "
"from a dummy atom group.
\n
"
);
}
std
::
vector
<
cvm
::
atom_pos
>
x
(
this
->
size
(),
0.0
);
cvm
::
atom_const_iter
ai
=
this
->
begin
();
std
::
vector
<
cvm
::
atom_pos
>::
iterator
xi
=
x
.
begin
();
for
(
;
ai
!=
this
->
end
();
++
xi
,
++
ai
)
{
*
xi
=
ai
->
pos
;
}
return
x
;
}
std
::
vector
<
cvm
::
atom_pos
>
cvm
::
atom_group
::
positions_shifted
(
cvm
::
rvector
const
&
shift
)
const
{
if
(
b_dummy
)
{
cvm
::
error
(
"Error: positions are not available "
"from a dummy atom group.
\n
"
);
}
std
::
vector
<
cvm
::
atom_pos
>
x
(
this
->
size
(),
0.0
);
cvm
::
atom_const_iter
ai
=
this
->
begin
();
std
::
vector
<
cvm
::
atom_pos
>::
iterator
xi
=
x
.
begin
();
for
(
;
ai
!=
this
->
end
();
++
xi
,
++
ai
)
{
*
xi
=
(
ai
->
pos
+
shift
);
}
return
x
;
}
std
::
vector
<
cvm
::
rvector
>
cvm
::
atom_group
::
velocities
()
const
{
if
(
b_dummy
)
{
cvm
::
error
(
"Error: velocities are not available "
"from a dummy atom group.
\n
"
);
}
std
::
vector
<
cvm
::
rvector
>
v
(
this
->
size
(),
0.0
);
cvm
::
atom_const_iter
ai
=
this
->
begin
();
std
::
vector
<
cvm
::
atom_pos
>::
iterator
vi
=
v
.
begin
();
for
(
;
ai
!=
this
->
end
();
vi
++
,
ai
++
)
{
*
vi
=
ai
->
vel
;
}
return
v
;
}
std
::
vector
<
cvm
::
rvector
>
cvm
::
atom_group
::
system_forces
()
const
{
if
(
b_dummy
)
{
cvm
::
error
(
"Error: system forces are not available "
"from a dummy atom group.
\n
"
);
}
std
::
vector
<
cvm
::
rvector
>
f
(
this
->
size
(),
0.0
);
cvm
::
atom_const_iter
ai
=
this
->
begin
();
std
::
vector
<
cvm
::
atom_pos
>::
iterator
fi
=
f
.
begin
();
for
(
;
ai
!=
this
->
end
();
++
fi
,
++
ai
)
{
*
fi
=
ai
->
system_force
;
}
return
f
;
}
cvm
::
rvector
cvm
::
atom_group
::
system_force
()
const
{
if
(
b_dummy
)
{
cvm
::
error
(
"Error: system forces are not available "
"from a dummy atom group.
\n
"
);
}
cvm
::
rvector
f
(
0.0
);
for
(
cvm
::
atom_const_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
f
+=
ai
->
system_force
;
}
return
f
;
}
void
cvm
::
atom_group
::
apply_colvar_force
(
cvm
::
real
const
&
force
)
{
if
(
b_dummy
)
return
;
if
(
noforce
)
{
cvm
::
error
(
"Error: sending a force to a group that has "
"
\"
enableForces
\"
set to off.
\n
"
);
return
;
}
if
(
b_rotate
)
{
// rotate forces back to the original frame
cvm
::
rotation
const
rot_inv
=
rot
.
inverse
();
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
ai
->
apply_force
(
rot_inv
.
rotate
(
force
*
ai
->
grad
));
}
}
else
{
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
ai
->
apply_force
(
force
*
ai
->
grad
);
}
}
if
((
b_center
||
b_rotate
)
&&
b_fit_gradients
)
{
atom_group
*
group_for_fit
=
ref_pos_group
?
ref_pos_group
:
this
;
// add the contribution from the roto-translational fit to the gradients
if
(
b_rotate
)
{
// rotate forces back to the original frame
cvm
::
rotation
const
rot_inv
=
rot
.
inverse
();
for
(
size_t
j
=
0
;
j
<
group_for_fit
->
size
();
j
++
)
{
(
*
group_for_fit
)[
j
].
apply_force
(
rot_inv
.
rotate
(
force
*
group_for_fit
->
fit_gradients
[
j
]));
}
}
else
{
for
(
size_t
j
=
0
;
j
<
group_for_fit
->
size
();
j
++
)
{
(
*
group_for_fit
)[
j
].
apply_force
(
force
*
group_for_fit
->
fit_gradients
[
j
]);
}
}
}
}
void
cvm
::
atom_group
::
apply_force
(
cvm
::
rvector
const
&
force
)
{
if
(
b_dummy
)
return
;
if
(
noforce
)
{
cvm
::
error
(
"Error: sending a force to a group that has "
"
\"
disableForces
\"
defined.
\n
"
);
return
;
}
if
(
b_rotate
)
{
cvm
::
rotation
const
rot_inv
=
rot
.
inverse
();
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
ai
->
apply_force
(
rot_inv
.
rotate
((
ai
->
mass
/
this
->
total_mass
)
*
force
));
}
}
else
{
for
(
cvm
::
atom_iter
ai
=
this
->
begin
();
ai
!=
this
->
end
();
ai
++
)
{
ai
->
apply_force
((
ai
->
mass
/
this
->
total_mass
)
*
force
);
}
}
}
void
cvm
::
atom_group
::
apply_forces
(
std
::
vector
<
cvm
::
rvector
>
const
&
forces
)
{
if
(
b_dummy
)
return
;
if
(
noforce
)
cvm
::
error
(
"Error: sending a force to a group that has "
"
\"
disableForces
\"
defined.
\n
"
);
if
(
forces
.
size
()
!=
this
->
size
())
{
cvm
::
error
(
"Error: trying to apply an array of forces to an atom "
"group which does not have the same length.
\n
"
);
}
if
(
b_rotate
)
{
cvm
::
rotation
const
rot_inv
=
rot
.
inverse
();
cvm
::
atom_iter
ai
=
this
->
begin
();
std
::
vector
<
cvm
::
rvector
>::
const_iterator
fi
=
forces
.
begin
();
for
(
;
ai
!=
this
->
end
();
++
fi
,
++
ai
)
{
ai
->
apply_force
(
rot_inv
.
rotate
(
*
fi
));
}
}
else
{
cvm
::
atom_iter
ai
=
this
->
begin
();
std
::
vector
<
cvm
::
rvector
>::
const_iterator
fi
=
forces
.
begin
();
for
(
;
ai
!=
this
->
end
();
++
fi
,
++
ai
)
{
ai
->
apply_force
(
*
fi
);
}
}
}
Event Timeline
Log In to Comment