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