Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85872485
colvarcomp_protein.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 2, 16:40
Size
11 KB
Mime Type
text/x-c
Expires
Fri, Oct 4, 16:40 (2 d)
Engine
blob
Format
Raw Data
Handle
21287715
Attached To
rLAMMPS lammps
colvarcomp_protein.cpp
View Options
// -*- c++ -*-
#include <cmath>
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
#include "colvar.h"
#include "colvarcomp.h"
//////////////////////////////////////////////////////////////////////
// alpha component
//////////////////////////////////////////////////////////////////////
// FIXME: this will not make collect_gradients work
// because gradients in individual atom groups
// are those of the sub-cvcs (angle, hb), not those
// of this cvc (alpha)
// This is true of all cvcs with sub-cvcs, and those
// that do not calculate explicit gradients
// SO: we need a flag giving the availability of
// atomic gradients
colvar
::
alpha_angles
::
alpha_angles
(
std
::
string
const
&
conf
)
:
cvc
(
conf
)
{
if
(
cvm
::
debug
())
cvm
::
log
(
"Initializing alpha_angles object.
\n
"
);
function_type
=
"alpha_angles"
;
x
.
type
(
colvarvalue
::
type_scalar
);
std
::
string
segment_id
;
get_keyval
(
conf
,
"psfSegID"
,
segment_id
,
std
::
string
(
"MAIN"
));
std
::
vector
<
int
>
residues
;
{
std
::
string
residues_conf
=
""
;
key_lookup
(
conf
,
"residueRange"
,
residues_conf
);
if
(
residues_conf
.
size
())
{
std
::
istringstream
is
(
residues_conf
);
int
initial
,
final
;
char
dash
;
if
(
(
is
>>
initial
)
&&
(
initial
>
0
)
&&
(
is
>>
dash
)
&&
(
dash
==
'-'
)
&&
(
is
>>
final
)
&&
(
final
>
0
)
)
{
for
(
int
rnum
=
initial
;
rnum
<=
final
;
rnum
++
)
{
residues
.
push_back
(
rnum
);
}
}
}
else
{
cvm
::
fatal_error
(
"Error: no residues defined in
\"
residueRange
\"
.
\n
"
);
}
}
if
(
residues
.
size
()
<
5
)
{
cvm
::
fatal_error
(
"Error: not enough residues defined in
\"
residueRange
\"
.
\n
"
);
}
std
::
string
const
&
sid
=
segment_id
;
std
::
vector
<
int
>
const
&
r
=
residues
;
get_keyval
(
conf
,
"hBondCoeff"
,
hb_coeff
,
0.5
);
if
(
(
hb_coeff
<
0.0
)
||
(
hb_coeff
>
1.0
)
)
{
cvm
::
fatal_error
(
"Error: hBondCoeff must be defined between 0 and 1.
\n
"
);
}
get_keyval
(
conf
,
"angleRef"
,
theta_ref
,
88.0
);
get_keyval
(
conf
,
"angleTol"
,
theta_tol
,
15.0
);
if
(
hb_coeff
<
1.0
)
{
for
(
size_t
i
=
0
;
i
<
residues
.
size
()
-
2
;
i
++
)
{
theta
.
push_back
(
new
colvar
::
angle
(
cvm
::
atom
(
r
[
i
],
"CA"
,
sid
),
cvm
::
atom
(
r
[
i
+
1
],
"CA"
,
sid
),
cvm
::
atom
(
r
[
i
+
2
],
"CA"
,
sid
)));
atom_groups
.
push_back
(
theta
.
back
()
->
atom_groups
[
0
]);
atom_groups
.
push_back
(
theta
.
back
()
->
atom_groups
[
1
]);
atom_groups
.
push_back
(
theta
.
back
()
->
atom_groups
[
2
]);
}
}
else
{
cvm
::
log
(
"The hBondCoeff specified will disable the Calpha-Calpha-Calpha angle terms.
\n
"
);
}
{
cvm
::
real
r0
;
size_t
en
,
ed
;
get_keyval
(
conf
,
"hBondCutoff"
,
r0
,
(
3.3
*
cvm
::
unit_angstrom
()));
get_keyval
(
conf
,
"hBondExpNumer"
,
en
,
6
);
get_keyval
(
conf
,
"hBondExpDenom"
,
ed
,
8
);
if
(
hb_coeff
>
0.0
)
{
for
(
size_t
i
=
0
;
i
<
residues
.
size
()
-
4
;
i
++
)
{
hb
.
push_back
(
new
colvar
::
h_bond
(
cvm
::
atom
(
r
[
i
],
"O"
,
sid
),
cvm
::
atom
(
r
[
i
+
4
],
"N"
,
sid
),
r0
,
en
,
ed
));
atom_groups
.
push_back
(
hb
.
back
()
->
atom_groups
[
0
]);
}
}
else
{
cvm
::
log
(
"The hBondCoeff specified will disable the hydrogen bond terms.
\n
"
);
}
}
if
(
cvm
::
debug
())
cvm
::
log
(
"Done initializing alpha_angles object.
\n
"
);
}
colvar
::
alpha_angles
::
alpha_angles
()
:
cvc
()
{
function_type
=
"alpha_angles"
;
x
.
type
(
colvarvalue
::
type_scalar
);
}
colvar
::
alpha_angles
::~
alpha_angles
()
{
while
(
theta
.
size
()
!=
0
)
{
delete
theta
.
back
();
theta
.
pop_back
();
}
while
(
hb
.
size
()
!=
0
)
{
delete
hb
.
back
();
hb
.
pop_back
();
}
}
void
colvar
::
alpha_angles
::
calc_value
()
{
x
.
real_value
=
0.0
;
if
(
theta
.
size
())
{
cvm
::
real
const
theta_norm
=
(
1.0
-
hb_coeff
)
/
cvm
::
real
(
theta
.
size
());
for
(
size_t
i
=
0
;
i
<
theta
.
size
();
i
++
)
{
(
theta
[
i
])
->
calc_value
();
cvm
::
real
const
t
=
((
theta
[
i
])
->
value
().
real_value
-
theta_ref
)
/
theta_tol
;
cvm
::
real
const
f
=
(
(
1.0
-
std
::
pow
(
t
,
(
int
)
2
))
/
(
1.0
-
std
::
pow
(
t
,
(
int
)
4
))
);
x
.
real_value
+=
theta_norm
*
f
;
if
(
cvm
::
debug
())
cvm
::
log
(
"Calpha-Calpha angle no. "
+
cvm
::
to_str
(
i
+
1
)
+
" in
\"
"
+
this
->
name
+
"
\"
has a value of "
+
(
cvm
::
to_str
((
theta
[
i
])
->
value
().
real_value
))
+
" degrees, f = "
+
cvm
::
to_str
(
f
)
+
".
\n
"
);
}
}
if
(
hb
.
size
())
{
cvm
::
real
const
hb_norm
=
hb_coeff
/
cvm
::
real
(
hb
.
size
());
for
(
size_t
i
=
0
;
i
<
hb
.
size
();
i
++
)
{
(
hb
[
i
])
->
calc_value
();
x
.
real_value
+=
hb_norm
*
(
hb
[
i
])
->
value
().
real_value
;
if
(
cvm
::
debug
())
cvm
::
log
(
"Hydrogen bond no. "
+
cvm
::
to_str
(
i
+
1
)
+
" in
\"
"
+
this
->
name
+
"
\"
has a value of "
+
(
cvm
::
to_str
((
hb
[
i
])
->
value
().
real_value
))
+
".
\n
"
);
}
}
}
void
colvar
::
alpha_angles
::
calc_gradients
()
{
size_t
i
;
for
(
i
=
0
;
i
<
theta
.
size
();
i
++
)
(
theta
[
i
])
->
calc_gradients
();
for
(
i
=
0
;
i
<
hb
.
size
();
i
++
)
(
hb
[
i
])
->
calc_gradients
();
}
void
colvar
::
alpha_angles
::
apply_force
(
colvarvalue
const
&
force
)
{
if
(
theta
.
size
())
{
cvm
::
real
const
theta_norm
=
(
1.0
-
hb_coeff
)
/
cvm
::
real
(
theta
.
size
());
for
(
size_t
i
=
0
;
i
<
theta
.
size
();
i
++
)
{
cvm
::
real
const
t
=
((
theta
[
i
])
->
value
().
real_value
-
theta_ref
)
/
theta_tol
;
cvm
::
real
const
f
=
(
(
1.0
-
std
::
pow
(
t
,
(
int
)
2
))
/
(
1.0
-
std
::
pow
(
t
,
(
int
)
4
))
);
cvm
::
real
const
dfdt
=
1.0
/
(
1.0
-
std
::
pow
(
t
,
(
int
)
4
))
*
(
(
-
2.0
*
t
)
+
(
-
1.0
*
f
)
*
(
-
4.0
*
std
::
pow
(
t
,
(
int
)
3
))
);
(
theta
[
i
])
->
apply_force
(
theta_norm
*
dfdt
*
(
1.0
/
theta_tol
)
*
force
.
real_value
);
}
}
if
(
hb
.
size
())
{
cvm
::
real
const
hb_norm
=
hb_coeff
/
cvm
::
real
(
hb
.
size
());
for
(
size_t
i
=
0
;
i
<
hb
.
size
();
i
++
)
{
(
hb
[
i
])
->
apply_force
(
0.5
*
hb_norm
*
force
.
real_value
);
}
}
}
//////////////////////////////////////////////////////////////////////
// dihedral principal component
//////////////////////////////////////////////////////////////////////
// FIXME: this will not make collect_gradients work
// because gradients in individual atom groups
// are those of the sub-cvcs (dihedral), not those
// of this cvc
// This is true of all cvcs with sub-cvcs, and those
// that do not calculate explicit gradients
// SO: we need a flag giving the availability of
// atomic gradients
colvar
::
dihedPC
::
dihedPC
(
std
::
string
const
&
conf
)
:
cvc
(
conf
)
{
if
(
cvm
::
debug
())
cvm
::
log
(
"Initializing dihedral PC object.
\n
"
);
function_type
=
"dihedPC"
;
x
.
type
(
colvarvalue
::
type_scalar
);
std
::
string
segment_id
;
get_keyval
(
conf
,
"psfSegID"
,
segment_id
,
std
::
string
(
"MAIN"
));
std
::
vector
<
int
>
residues
;
{
std
::
string
residues_conf
=
""
;
key_lookup
(
conf
,
"residueRange"
,
residues_conf
);
if
(
residues_conf
.
size
())
{
std
::
istringstream
is
(
residues_conf
);
int
initial
,
final
;
char
dash
;
if
(
(
is
>>
initial
)
&&
(
initial
>
0
)
&&
(
is
>>
dash
)
&&
(
dash
==
'-'
)
&&
(
is
>>
final
)
&&
(
final
>
0
)
)
{
for
(
int
rnum
=
initial
;
rnum
<=
final
;
rnum
++
)
{
residues
.
push_back
(
rnum
);
}
}
}
else
{
cvm
::
fatal_error
(
"Error: no residues defined in
\"
residueRange
\"
.
\n
"
);
}
}
if
(
residues
.
size
()
<
2
)
{
cvm
::
fatal_error
(
"Error: dihedralPC requires at least two residues.
\n
"
);
}
std
::
string
const
&
sid
=
segment_id
;
std
::
vector
<
int
>
const
&
r
=
residues
;
std
::
string
vecFileName
;
int
vecNumber
;
if
(
get_keyval
(
conf
,
"vectorFile"
,
vecFileName
,
vecFileName
))
{
get_keyval
(
conf
,
"vectorNumber"
,
vecNumber
,
0
);
if
(
vecNumber
<
1
)
cvm
::
fatal_error
(
"A positive value of vectorNumber is required."
);
std
::
ifstream
vecFile
;
vecFile
.
open
(
vecFileName
.
c_str
());
if
(
!
vecFile
.
good
())
cvm
::
fatal_error
(
"Error opening dihedral PCA vector file "
+
vecFileName
+
" for reading"
);
// TODO: adapt to different formats by setting this flag
bool
eigenvectors_as_columns
=
true
;
if
(
eigenvectors_as_columns
)
{
// Carma-style dPCA file
std
::
string
line
;
cvm
::
real
c
;
while
(
vecFile
.
good
())
{
getline
(
vecFile
,
line
);
if
(
line
.
length
()
<
2
)
break
;
std
::
istringstream
ls
(
line
);
for
(
int
i
=
0
;
i
<
vecNumber
;
i
++
)
ls
>>
c
;
coeffs
.
push_back
(
c
);
}
}
/* TODO Uncomment this when different formats are recognized
else {
// Eigenvectors as lines
// Skip to the right line
for (int i = 1; i<vecNumber; i++)
vecFile.ignore(999999, '\n');
if (!vecFile.good())
cvm::fatal_error("Error reading dihedral PCA vector file " + vecFileName);
std::string line;
getline(vecFile, line);
std::istringstream ls(line);
cvm::real c;
while (ls.good()) {
ls >> c;
coeffs.push_back(c);
}
}
*/
vecFile
.
close
();
}
else
{
get_keyval
(
conf
,
"vector"
,
coeffs
,
coeffs
);
}
if
(
coeffs
.
size
()
!=
4
*
(
residues
.
size
()
-
1
))
{
cvm
::
fatal_error
(
"Error: wrong number of coefficients: "
+
cvm
::
to_str
(
coeffs
.
size
())
+
". Expected "
+
cvm
::
to_str
(
4
*
(
residues
.
size
()
-
1
))
+
" (4 coeffs per residue, minus one residue).
\n
"
);
}
for
(
size_t
i
=
0
;
i
<
residues
.
size
()
-
1
;
i
++
)
{
// Psi
theta
.
push_back
(
new
colvar
::
dihedral
(
cvm
::
atom
(
r
[
i
],
"N"
,
sid
),
cvm
::
atom
(
r
[
i
],
"CA"
,
sid
),
cvm
::
atom
(
r
[
i
],
"C"
,
sid
),
cvm
::
atom
(
r
[
i
+
1
],
"N"
,
sid
)));
// Phi (next res)
theta
.
push_back
(
new
colvar
::
dihedral
(
cvm
::
atom
(
r
[
i
],
"C"
,
sid
),
cvm
::
atom
(
r
[
i
+
1
],
"N"
,
sid
),
cvm
::
atom
(
r
[
i
+
1
],
"CA"
,
sid
),
cvm
::
atom
(
r
[
i
+
1
],
"C"
,
sid
)));
}
if
(
cvm
::
debug
())
cvm
::
log
(
"Done initializing dihedPC object.
\n
"
);
}
colvar
::
dihedPC
::
dihedPC
()
:
cvc
()
{
function_type
=
"dihedPC"
;
x
.
type
(
colvarvalue
::
type_scalar
);
}
colvar
::
dihedPC
::~
dihedPC
()
{
while
(
theta
.
size
()
!=
0
)
{
delete
theta
.
back
();
theta
.
pop_back
();
}
}
void
colvar
::
dihedPC
::
calc_value
()
{
x
.
real_value
=
0.0
;
for
(
size_t
i
=
0
;
i
<
theta
.
size
();
i
++
)
{
theta
[
i
]
->
calc_value
();
cvm
::
real
const
t
=
(
PI
/
180.
)
*
theta
[
i
]
->
value
().
real_value
;
x
.
real_value
+=
coeffs
[
2
*
i
]
*
std
::
cos
(
t
)
+
coeffs
[
2
*
i
+
1
]
*
std
::
sin
(
t
);
}
}
void
colvar
::
dihedPC
::
calc_gradients
()
{
for
(
size_t
i
=
0
;
i
<
theta
.
size
();
i
++
)
{
theta
[
i
]
->
calc_gradients
();
}
}
void
colvar
::
dihedPC
::
apply_force
(
colvarvalue
const
&
force
)
{
for
(
size_t
i
=
0
;
i
<
theta
.
size
();
i
++
)
{
cvm
::
real
const
t
=
(
PI
/
180.
)
*
theta
[
i
]
->
value
().
real_value
;
cvm
::
real
const
dcosdt
=
-
(
PI
/
180.
)
*
std
::
sin
(
t
);
cvm
::
real
const
dsindt
=
(
PI
/
180.
)
*
std
::
cos
(
t
);
theta
[
i
]
->
apply_force
((
coeffs
[
2
*
i
]
*
dcosdt
+
coeffs
[
2
*
i
+
1
]
*
dsindt
)
*
force
);
}
}
Event Timeline
Log In to Comment