00001 #include <cmath>
00002
00003 #include "colvarmodule.h"
00004 #include "colvarvalue.h"
00005 #include "colvarparse.h"
00006 #include "colvar.h"
00007 #include "colvarcomp.h"
00008
00009
00011
00013
00014 colvar::alpha_angles::alpha_angles (std::string const &conf)
00015 : cvc (conf)
00016 {
00017 if (cvm::debug())
00018 cvm::log ("Initializing alpha_angles object.\n");
00019
00020 function_type = "alpha_angles";
00021 x.type (colvarvalue::type_scalar);
00022
00023 std::string segment_id;
00024 get_keyval (conf, "psfSegID", segment_id, std::string ("MAIN"));
00025
00026 std::vector<int> residues;
00027 {
00028 std::string residues_conf = "";
00029 key_lookup (conf, "residueRange", residues_conf);
00030 if (residues_conf.size()) {
00031 std::istringstream is (residues_conf);
00032 int initial, final;
00033 char dash;
00034 if ( (is >> initial) && (initial > 0) &&
00035 (is >> dash) && (dash == '-') &&
00036 (is >> final) && (final > 0) ) {
00037 for (int rnum = initial; rnum <= final; rnum++) {
00038 residues.push_back (rnum);
00039 }
00040 }
00041 } else {
00042 cvm::fatal_error ("Error: no residues defined in \"residueRange\".\n");
00043 }
00044 }
00045
00046 if (residues.size() < 5) {
00047 cvm::fatal_error ("Error: not enough residues defined in \"residueRange\".\n");
00048 }
00049
00050 std::string const &sid = segment_id;
00051 std::vector<int> const &r = residues;
00052
00053
00054 get_keyval (conf, "hBondCoeff", hb_coeff, 0.5);
00055 if ( (hb_coeff < 0.0) || (hb_coeff > 1.0) ) {
00056 cvm::fatal_error ("Error: hBondCoeff must be defined between 0 and 1.\n");
00057 }
00058
00059
00060 get_keyval (conf, "angleRef", theta_ref, 88.0, parse_silent);
00061 get_keyval (conf, "angleTol", theta_tol, 15.0, parse_silent);
00062
00063 if (hb_coeff < 1.0) {
00064
00065 for (size_t i = 0; i < residues.size()-2; i++) {
00066 theta.push_back (new colvar::angle (cvm::atom (r[i ], "CA", sid),
00067 cvm::atom (r[i+1], "CA", sid),
00068 cvm::atom (r[i+2], "CA", sid)));
00069 }
00070
00071 } else {
00072 cvm::log ("The hBondCoeff specified will disable the Calpha-Calpha-Calpha angle terms.\n");
00073 }
00074
00075 {
00076 cvm::real r0;
00077 size_t en, ed;
00078 get_keyval (conf, "hBondCutoff", r0, (3.3 * cvm::unit_angstrom()));
00079 get_keyval (conf, "hBondExpNumer", en, 6);
00080 get_keyval (conf, "hBondExpDenom", ed, 8);
00081
00082 if (hb_coeff > 0.0) {
00083
00084 for (size_t i = 0; i < residues.size()-4; i++) {
00085 hb.push_back (new colvar::h_bond (cvm::atom (r[i ], "O", sid),
00086 cvm::atom (r[i+4], "N", sid),
00087 r0, en, ed));
00088 }
00089
00090 } else {
00091 cvm::log ("The hBondCoeff specified will disable the hydrogen bond terms.\n");
00092 }
00093 }
00094
00095 if (cvm::debug())
00096 cvm::log ("Done initializing alpha_angles object.\n");
00097 }
00098
00099
00100 colvar::alpha_angles::alpha_angles()
00101 : cvc ()
00102 {
00103 function_type = "alpha_angles";
00104 x.type (colvarvalue::type_scalar);
00105 }
00106
00107
00108 void colvar::alpha_angles::calc_value()
00109 {
00110 x.real_value = 0.0;
00111
00112 if (theta.size()) {
00113
00114 cvm::real const theta_norm =
00115 (1.0-hb_coeff) / cvm::real (theta.size());
00116
00117 for (size_t i = 0; i < theta.size(); i++) {
00118
00119 (theta[i])->calc_value();
00120
00121 cvm::real const t = ((theta[i])->value().real_value-theta_ref)/theta_tol;
00122 cvm::real const f = ( (1.0 - ::pow (t, (int) 2)) /
00123 (1.0 - ::pow (t, (int) 4)) );
00124
00125 x.real_value += theta_norm * f;
00126
00127 if (cvm::debug())
00128 cvm::log ("Calpha-Calpha angle no. "+cvm::to_str (i+1)+" in \""+
00129 this->name+"\" has a value of "+
00130 (cvm::to_str ((theta[i])->value().real_value))+
00131 " degrees, f = "+cvm::to_str (f)+".\n");
00132 }
00133 }
00134
00135 if (hb.size()) {
00136
00137 cvm::real const hb_norm =
00138 hb_coeff / cvm::real (hb.size());
00139
00140 for (size_t i = 0; i < hb.size(); i++) {
00141 (hb[i])->calc_value();
00142 x.real_value += hb_norm * (hb[i])->value().real_value;
00143 if (cvm::debug())
00144 cvm::log ("Hydrogen bond no. "+cvm::to_str (i+1)+" in \""+
00145 this->name+"\" has a value of "+
00146 (cvm::to_str ((hb[i])->value().real_value))+".\n");
00147 }
00148 }
00149 }
00150
00151
00152 void colvar::alpha_angles::calc_gradients()
00153 {
00154 for (size_t i = 0; i < theta.size(); i++)
00155 (theta[i])->calc_gradients();
00156
00157 for (size_t i = 0; i < hb.size(); i++)
00158 (hb[i])->calc_gradients();
00159 }
00160
00161
00162 void colvar::alpha_angles::apply_force (colvarvalue const &force)
00163 {
00164
00165 if (theta.size()) {
00166
00167 cvm::real const theta_norm =
00168 (1.0-hb_coeff) / cvm::real (theta.size());
00169
00170 for (size_t i = 0; i < theta.size(); i++) {
00171
00172 cvm::real const t = ((theta[i])->value().real_value-theta_ref)/theta_tol;
00173 cvm::real const f = ( (1.0 - ::pow (t, (int) 2)) /
00174 (1.0 - ::pow (t, (int) 4)) );
00175
00176 cvm::real const dfdt =
00177 1.0/(1.0 - ::pow (t, (int) 4)) *
00178 ( (-2.0 * t) + (-1.0*f)*(-4.0 * ::pow (t, (int) 3)) );
00179
00180 (theta[i])->apply_force (theta_norm *
00181 dfdt * (1.0/theta_tol) *
00182 force.real_value );
00183 }
00184 }
00185
00186 if (hb.size()) {
00187
00188 cvm::real const hb_norm =
00189 hb_coeff / cvm::real (hb.size());
00190
00191 for (size_t i = 0; i < hb.size(); i++) {
00192 (hb[i])->apply_force (0.5 * hb_norm * force.real_value);
00193 }
00194 }
00195 }
00196
00197
00198
00200
00202
00203
00204
00205 inline cvm::real dih_func (cvm::real const &w,
00206 cvm::real const &w0)
00207 {
00208 return 0.5*(1.0 + ::cos (w-w0));
00209 }
00210
00211
00212 inline cvm::real dih_deriv (cvm::real const &w,
00213 cvm::real const &w0)
00214 {
00215 return (-0.5 * ::sin (w-w0));
00216 }
00217
00218
00219 colvar::alpha_dihedrals::alpha_dihedrals (std::string const &conf)
00220 : cvc (conf)
00221 {
00222 if (cvm::debug())
00223 cvm::log ("Initializing alpha_dihedrals object.\n");
00224
00225 function_type = "alpha_dihedrals";
00226 x.type (colvarvalue::type_scalar);
00227
00228 std::vector<int> residues;
00229
00230 get_keyval (conf, "residues", residues, std::vector<int>());
00231
00232 phi.reserve (residues.size());
00233 psi.reserve (residues.size());
00234 hb.reserve (residues.size());
00235
00236 get_keyval (conf, "phi_ref", phi_ref, -57.8, parse_silent);
00237 get_keyval (conf, "psi_ref", psi_ref, -47.0, parse_silent);
00238
00239 if (residues.size() < 5) {
00240 cvm::fatal_error ("Error: not enough residues defined.\n");
00241 }
00242
00243 std::string segment_id;
00244 get_keyval (conf, "segment_id", segment_id, std::string ("MAIN"));
00245
00246 std::string const &sid = segment_id;
00247 std::vector<int> const &r = residues;
00248 for (size_t i = 0; i < residues.size()-1; i++) {
00249
00250 phi.push_back (new colvar::dihedral (cvm::atom (r[i ], "C", sid),
00251 cvm::atom (r[i+1], "N", sid),
00252 cvm::atom (r[i+1], "CA", sid),
00253 cvm::atom (r[i+1], "C", sid)));
00254
00255 psi.push_back (new colvar::dihedral (cvm::atom (r[i ], "N", sid),
00256 cvm::atom (r[i ], "CA", sid),
00257 cvm::atom (r[i ], "C", sid),
00258 cvm::atom (r[i+1], "N", sid)));
00259 }
00260
00261 for (size_t i = 0; i < residues.size()-4; i++) {
00262 hb.push_back (new colvar::h_bond (cvm::atom (r[i ], "O", sid),
00263 cvm::atom (r[i+4], "N", sid),
00264 3.3 * cvm::unit_angstrom(),
00265 6, 8));
00266 }
00267
00268 if (cvm::debug())
00269 cvm::log ("Done initializing alpha_dihedrals object.\n");
00270 }
00271
00272
00273 colvar::alpha_dihedrals::alpha_dihedrals()
00274 : cvc ()
00275 {
00276 function_type = "alpha_dihedrals";
00277 x.type (colvarvalue::type_scalar);
00278 }
00279
00280
00281 void colvar::alpha_dihedrals::calc_value()
00282 {
00283 x.real_value = 0.0;
00284
00285 for (size_t i = 0; i < phi.size(); i++) {
00286
00287 (phi[i])->calc_value();
00288 (psi[i])->calc_value();
00289
00290 x.real_value += 0.5 *
00291 dih_func (((phi[i])->value()).real_value, phi_ref) *
00292 dih_func (((psi[i])->value()).real_value, psi_ref);
00293
00294 if (cvm::debug())
00295 cvm::log ("Phi dihedral no. "+cvm::to_str (i+1)+" in \""+
00296 this->name+"\" has a value of "+
00297 (cvm::to_str ((phi[i])->value().real_value))+
00298 " degrees.\n");
00299
00300 if (cvm::debug())
00301 cvm::log ("Psi dihedral no. "+cvm::to_str (i+1)+" in \""+
00302 this->name+"\" has a value of "+
00303 (cvm::to_str ((psi[i])->value().real_value))+
00304 " degrees.\n");
00305 }
00306
00307 for (size_t i = 0; i < hb.size(); i++) {
00308 (hb[i])->calc_value();
00309 x.real_value += 0.5 * (hb[i])->value().real_value;
00310 if (cvm::debug())
00311 cvm::log ("Hydrogen bond no. "+cvm::to_str (i+1)+" in \""+
00312 this->name+"\" has a value of "+
00313 (cvm::to_str ((hb[i])->value().real_value))+".\n");
00314 }
00315 }
00316
00317
00318 void colvar::alpha_dihedrals::calc_gradients()
00319 {
00320 for (size_t i = 0; i < phi.size(); i++) {
00321 (phi[i])->calc_gradients();
00322 (psi[i])->calc_gradients();
00323 }
00324 for (size_t i = 0; i < hb.size(); i++) {
00325 (hb[i])->calc_gradients();
00326 }
00327 }
00328
00329
00330 void colvar::alpha_dihedrals::apply_force (colvarvalue const &force)
00331 {
00332 for (size_t i = 0; i < phi.size(); i++) {
00333
00334 (phi[i])->apply_force ( 0.5 *
00335 dih_func (((psi[i])->value()).real_value, psi_ref) *
00336 dih_deriv (((phi[i])->value()).real_value, phi_ref) *
00337 force.real_value );
00338
00339 (psi[i])->apply_force ( 0.5 *
00340 dih_func (((phi[i])->value()).real_value, phi_ref) *
00341 dih_deriv (((psi[i])->value()).real_value, psi_ref) *
00342 force.real_value );
00343 }
00344
00345 for (size_t i = 0; i < hb.size(); i++) {
00346 (hb[i])->apply_force (0.5 * force.real_value);
00347 }
00348 }
00349
00350
00351
