00001
00002 #include "colvarmodule.h"
00003 #include "colvar.h"
00004 #include "colvarcomp.h"
00005
00006 #include <cmath>
00007
00008
00009
00010 colvar::angle::angle (std::string const &conf)
00011 : cvc (conf)
00012 {
00013 function_type = "angle";
00014 b_inverse_gradients = true;
00015 b_Jacobian_derivative = true;
00016 parse_group (conf, "group1", group1);
00017 parse_group (conf, "group2", group2);
00018 parse_group (conf, "group3", group3);
00019 atom_groups.push_back (&group1);
00020 atom_groups.push_back (&group2);
00021 atom_groups.push_back (&group3);
00022 if (get_keyval (conf, "oneSiteSystemForce", b_1site_force, false)) {
00023 cvm::log ("Computing system force on group 1 only");
00024 }
00025 x.type (colvarvalue::type_scalar);
00026 }
00027
00028
00029 colvar::angle::angle (cvm::atom const &a1,
00030 cvm::atom const &a2,
00031 cvm::atom const &a3)
00032 : group1 (std::vector<cvm::atom> (1, a1)),
00033 group2 (std::vector<cvm::atom> (1, a2)),
00034 group3 (std::vector<cvm::atom> (1, a3))
00035 {
00036 function_type = "angle";
00037 b_inverse_gradients = true;
00038 b_Jacobian_derivative = true;
00039 b_1site_force = false;
00040 atom_groups.push_back (&group1);
00041 atom_groups.push_back (&group2);
00042 atom_groups.push_back (&group3);
00043
00044 x.type (colvarvalue::type_scalar);
00045 }
00046
00047
00048 colvar::angle::angle()
00049 {
00050 function_type = "angle";
00051 x.type (colvarvalue::type_scalar);
00052 }
00053
00054
00055 void colvar::angle::calc_value()
00056 {
00057 group1.read_positions();
00058 group2.read_positions();
00059 group3.read_positions();
00060
00061 cvm::atom_pos const g1_pos = group1.center_of_mass();
00062 cvm::atom_pos const g2_pos = group2.center_of_mass();
00063 cvm::atom_pos const g3_pos = group3.center_of_mass();
00064
00065 r21 = cvm::position_distance (g2_pos, g1_pos);
00066 r21l = r21.norm();
00067 r23 = cvm::position_distance (g2_pos, g3_pos);
00068 r23l = r23.norm();
00069
00070 cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
00071
00072 x.real_value = (180.0/PI) * std::acos (cos_theta);
00073 }
00074
00075
00076 void colvar::angle::calc_gradients()
00077 {
00078 cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
00079 cvm::real const dxdcos = -1.0 / std::sqrt (1.0 - cos_theta*cos_theta);
00080
00081 dxdr1 = (180.0/PI) * dxdcos *
00082 (1.0/r21l) * ( r23/r23l + (-1.0) * cos_theta * r21/r21l );
00083
00084 dxdr3 = (180.0/PI) * dxdcos *
00085 (1.0/r23l) * ( r21/r21l + (-1.0) * cos_theta * r23/r23l );
00086
00087 for (size_t i = 0; i < group1.size(); i++) {
00088 group1[i].grad = (group1[i].mass/group1.total_mass) *
00089 (dxdr1);
00090 }
00091
00092 for (size_t i = 0; i < group2.size(); i++) {
00093 group2[i].grad = (group2[i].mass/group2.total_mass) *
00094 (dxdr1 + dxdr3) * (-1.0);
00095 }
00096
00097 for (size_t i = 0; i < group3.size(); i++) {
00098 group3[i].grad = (group3[i].mass/group3.total_mass) *
00099 (dxdr3);
00100 }
00101 }
00102
00103 void colvar::angle::calc_force_invgrads()
00104 {
00105
00106
00107
00108
00109
00110
00111 if (b_1site_force) {
00112 group1.read_system_forces();
00113 cvm::real norm_fact = 1.0 / dxdr1.norm2();
00114 ft.real_value = norm_fact * dxdr1 * group1.system_force();
00115 } else {
00116 group1.read_system_forces();
00117 group3.read_system_forces();
00118 cvm::real norm_fact = 1.0 / (dxdr1.norm2() + dxdr3.norm2());
00119 ft.real_value = norm_fact * ( dxdr1 * group1.system_force()
00120 + dxdr3 * group3.system_force());
00121 }
00122 return;
00123 }
00124
00125 void colvar::angle::calc_Jacobian_derivative()
00126 {
00127
00128
00129 const cvm::real theta = x.real_value * PI / 180.0;
00130 jd = PI / 180.0 * (theta != 0.0 ? std::cos(theta) / std::sin(theta) : 0.0);
00131 }
00132
00133
00134 void colvar::angle::apply_force (colvarvalue const &force)
00135 {
00136 if (!group1.noforce)
00137 group1.apply_colvar_force (force.real_value);
00138
00139 if (!group2.noforce)
00140 group2.apply_colvar_force (force.real_value);
00141
00142 if (!group3.noforce)
00143 group3.apply_colvar_force (force.real_value);
00144 }
00145
00146
00147
00148
00149 colvar::dihedral::dihedral (std::string const &conf)
00150 : cvc (conf)
00151 {
00152 function_type = "dihedral";
00153 period = 360.0;
00154 b_periodic = true;
00155 b_inverse_gradients = true;
00156 b_Jacobian_derivative = true;
00157 if (get_keyval (conf, "oneSiteSystemForce", b_1site_force, false)) {
00158 cvm::log ("Computing system force on group 1 only");
00159 }
00160 parse_group (conf, "group1", group1);
00161 parse_group (conf, "group2", group2);
00162 parse_group (conf, "group3", group3);
00163 parse_group (conf, "group4", group4);
00164 atom_groups.push_back (&group1);
00165 atom_groups.push_back (&group2);
00166 atom_groups.push_back (&group3);
00167 atom_groups.push_back (&group4);
00168
00169 x.type (colvarvalue::type_scalar);
00170 }
00171
00172
00173 colvar::dihedral::dihedral (cvm::atom const &a1,
00174 cvm::atom const &a2,
00175 cvm::atom const &a3,
00176 cvm::atom const &a4)
00177 : group1 (std::vector<cvm::atom> (1, a1)),
00178 group2 (std::vector<cvm::atom> (1, a2)),
00179 group3 (std::vector<cvm::atom> (1, a3)),
00180 group4 (std::vector<cvm::atom> (1, a4))
00181 {
00182 if (cvm::debug())
00183 cvm::log ("Initializing dihedral object from atom groups.\n");
00184
00185 function_type = "dihedral";
00186 period = 360.0;
00187 b_periodic = true;
00188 b_inverse_gradients = true;
00189 b_Jacobian_derivative = true;
00190 b_1site_force = false;
00191
00192 atom_groups.push_back (&group1);
00193 atom_groups.push_back (&group2);
00194 atom_groups.push_back (&group3);
00195 atom_groups.push_back (&group4);
00196
00197 x.type (colvarvalue::type_scalar);
00198
00199 if (cvm::debug())
00200 cvm::log ("Done initializing dihedral object from atom groups.\n");
00201 }
00202
00203
00204 colvar::dihedral::dihedral()
00205 {
00206 function_type = "dihedral";
00207 period = 360.0;
00208 b_periodic = true;
00209 b_inverse_gradients = true;
00210 b_Jacobian_derivative = true;
00211 x.type (colvarvalue::type_scalar);
00212 }
00213
00214
00215 void colvar::dihedral::calc_value()
00216 {
00217 group1.read_positions();
00218 group2.read_positions();
00219 group3.read_positions();
00220 group4.read_positions();
00221
00222 cvm::atom_pos const g1_pos = group1.center_of_mass();
00223 cvm::atom_pos const g2_pos = group2.center_of_mass();
00224 cvm::atom_pos const g3_pos = group3.center_of_mass();
00225 cvm::atom_pos const g4_pos = group4.center_of_mass();
00226
00227
00228 r12 = cvm::position_distance (g1_pos, g2_pos);
00229 r23 = cvm::position_distance (g2_pos, g3_pos);
00230 r34 = cvm::position_distance (g3_pos, g4_pos);
00231
00232 cvm::rvector const n1 = cvm::rvector::outer (r12, r23);
00233 cvm::rvector const n2 = cvm::rvector::outer (r23, r34);
00234
00235 cvm::real const cos_phi = n1 * n2;
00236 cvm::real const sin_phi = n1 * r34 * r23.norm();
00237
00238 x.real_value = (180.0/PI) * std::atan2 (sin_phi, cos_phi);
00239 this->wrap (x);
00240 }
00241
00242
00243 void colvar::dihedral::calc_gradients()
00244 {
00245 cvm::rvector A = cvm::rvector::outer (r12, r23);
00246 cvm::real rA = A.norm();
00247 cvm::rvector B = cvm::rvector::outer (r23, r34);
00248 cvm::real rB = B.norm();
00249 cvm::rvector C = cvm::rvector::outer (r23, A);
00250 cvm::real rC = C.norm();
00251
00252 cvm::real const cos_phi = (A*B)/(rA*rB);
00253 cvm::real const sin_phi = (C*B)/(rC*rB);
00254
00255 cvm::rvector f1, f2, f3;
00256
00257 rB = 1.0/rB;
00258 B *= rB;
00259
00260 if (std::fabs (sin_phi) > 0.1) {
00261 rA = 1.0/rA;
00262 A *= rA;
00263 cvm::rvector const dcosdA = rA*(cos_phi*A-B);
00264 cvm::rvector const dcosdB = rB*(cos_phi*B-A);
00265 rA = 1.0;
00266
00267 cvm::real const K = (1.0/sin_phi) * (180.0/PI);
00268
00269 f1 = K * cvm::rvector::outer (r23, dcosdA);
00270 f3 = K * cvm::rvector::outer (dcosdB, r23);
00271 f2 = K * (cvm::rvector::outer (dcosdA, r12)
00272 + cvm::rvector::outer (r34, dcosdB));
00273 }
00274 else {
00275 rC = 1.0/rC;
00276 C *= rC;
00277 cvm::rvector const dsindC = rC*(sin_phi*C-B);
00278 cvm::rvector const dsindB = rB*(sin_phi*B-C);
00279 rC = 1.0;
00280
00281 cvm::real const K = (-1.0/cos_phi) * (180.0/PI);
00282
00283 f1.x = K*((r23.y*r23.y + r23.z*r23.z)*dsindC.x
00284 - r23.x*r23.y*dsindC.y
00285 - r23.x*r23.z*dsindC.z);
00286 f1.y = K*((r23.z*r23.z + r23.x*r23.x)*dsindC.y
00287 - r23.y*r23.z*dsindC.z
00288 - r23.y*r23.x*dsindC.x);
00289 f1.z = K*((r23.x*r23.x + r23.y*r23.y)*dsindC.z
00290 - r23.z*r23.x*dsindC.x
00291 - r23.z*r23.y*dsindC.y);
00292
00293 f3 = cvm::rvector::outer (dsindB, r23);
00294 f3 *= K;
00295
00296 f2.x = K*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x
00297 +(2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y
00298 +(2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z
00299 +dsindB.z*r34.y - dsindB.y*r34.z);
00300 f2.y = K*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y
00301 +(2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z
00302 +(2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x
00303 +dsindB.x*r34.z - dsindB.z*r34.x);
00304 f2.z = K*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z
00305 +(2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x
00306 +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y
00307 +dsindB.y*r34.x - dsindB.x*r34.y);
00308 }
00309
00310 for (size_t i = 0; i < group1.size(); i++)
00311 group1[i].grad = (group1[i].mass/group1.total_mass) * (-f1);
00312
00313 for (size_t i = 0; i < group2.size(); i++)
00314 group2[i].grad = (group2[i].mass/group2.total_mass) * (-f2 + f1);
00315
00316 for (size_t i = 0; i < group3.size(); i++)
00317 group3[i].grad = (group3[i].mass/group3.total_mass) * (-f3 + f2);
00318
00319 for (size_t i = 0; i < group4.size(); i++)
00320 group4[i].grad = (group4[i].mass/group4.total_mass) * (f3);
00321 }
00322
00323
00324 void colvar::dihedral::calc_force_invgrads()
00325 {
00326 cvm::rvector const u12 = r12.unit();
00327 cvm::rvector const u23 = r23.unit();
00328 cvm::rvector const u34 = r34.unit();
00329
00330 cvm::real const d12 = r12.norm();
00331 cvm::real const d34 = r34.norm();
00332
00333 cvm::rvector const cross1 = (cvm::rvector::outer (u23, u12)).unit();
00334 cvm::rvector const cross4 = (cvm::rvector::outer (u23, u34)).unit();
00335
00336 cvm::real const dot1 = u23 * u12;
00337 cvm::real const dot4 = u23 * u34;
00338
00339 cvm::real const fact1 = d12 * std::sqrt (1.0 - dot1 * dot1);
00340 cvm::real const fact4 = d34 * std::sqrt (1.0 - dot4 * dot4);
00341
00342 group1.read_system_forces();
00343 if ( b_1site_force ) {
00344
00345 ft.real_value = PI/180.0 * fact1 * (cross1 * group1.system_force());
00346 } else {
00347
00348 group4.read_system_forces();
00349 ft.real_value = PI/180.0 * 0.5 * (fact1 * (cross1 * group1.system_force())
00350 + fact4 * (cross4 * group4.system_force()));
00351 }
00352 }
00353
00354
00355 void colvar::dihedral::calc_Jacobian_derivative()
00356 {
00357
00358 jd = 0.0;
00359 }
00360
00361
00362 void colvar::dihedral::apply_force (colvarvalue const &force)
00363 {
00364 if (!group1.noforce)
00365 group1.apply_colvar_force (force.real_value);
00366
00367 if (!group2.noforce)
00368 group2.apply_colvar_force (force.real_value);
00369
00370 if (!group3.noforce)
00371 group3.apply_colvar_force (force.real_value);
00372
00373 if (!group4.noforce)
00374 group4.apply_colvar_force (force.real_value);
00375 }
00376
00377