pele
Python energy landscape explorer
/home/js850/projects/pele/source/aatopology.cpp
00001 #include "pele/rotations.h"
00002 #include "pele/aatopology.h"
00003 #include "pele/vecn.h"
00004 
00005 namespace pele{
00006 
00007 pele::Array<double>
00008 pele::RigidFragment::to_atomistic(pele::Array<double> const com,
00009         pele::VecN<3> const & p)
00010 {
00011     assert(com.size() == _ndim);
00012     assert(p.size() == 3);
00013     auto rmat = pele::aa_to_rot_mat(p);
00014     Array<double> pos(_atom_positions.size());
00015     MatrixAdapter<double> mpos(pos, _ndim);
00016 
00017     // in python this is:
00018     //      return com + np.dot(R, np.transpose(self.atom_positions)).transpose()
00019     for (size_t atom = 0; atom<_natoms; ++atom) {
00020         for (size_t j = 0; j<_ndim; ++j) {
00021             double val = com[j];
00022             for (size_t k = 0; k<_ndim; ++k) {
00023                 val += rmat(j,k) * _atom_positions_matrix(atom,k);
00024             }
00025             mpos(atom, j) = val;
00026         }
00027     }
00028     return pos;
00029 }
00030 
00031 void
00032 pele::RigidFragment::transform_grad(
00033         pele::VecN<3> const & p,
00034         pele::Array<double> const g,
00035         pele::VecN<3> & g_com,
00036         pele::VecN<3> & g_rot
00037         )
00038 {
00039     assert(g.size() == natoms() * 3);
00040     // view the array as a matrix
00041     MatrixAdapter<double> gmat(g, 3);
00042 
00043     // compute the rotation matrix and derivatives
00044     pele::MatrixNM<3,3> rmat;
00045     pele::MatrixNM<3,3> drm1;
00046     pele::MatrixNM<3,3> drm2;
00047     pele::MatrixNM<3,3> drm3;
00048     rot_mat_derivatives(p, rmat, drm1, drm2, drm3);
00049 
00050     // do the center of mass coordinates
00051     for (size_t k=0; k<3; ++k) {
00052         double val = 0;
00053         for (size_t atom=0; atom < _natoms; ++atom) {
00054             val += gmat(atom,k);
00055         }
00056         g_com[k] = val;
00057     }
00058 
00059     // now do the rotations
00060     g_rot.assign(0);
00061     for (size_t atom=0; atom < _natoms; ++atom) {
00062         double val1 = 0;
00063         double val2 = 0;
00064         double val3 = 0;
00065         for (size_t i=0; i<3; ++i) {
00066             for (size_t j=0; j<3; ++j) {
00067                 val1 += gmat(atom,i) * drm1(i,j) * _atom_positions_matrix(atom,j);
00068                 val2 += gmat(atom,i) * drm2(i,j) * _atom_positions_matrix(atom,j);
00069                 val3 += gmat(atom,i) * drm3(i,j) * _atom_positions_matrix(atom,j);
00070             }
00071         }
00072         g_rot[0] += val1;
00073         g_rot[1] += val2;
00074         g_rot[2] += val3;
00075     }
00076 }
00077 
00078 double
00079 pele::RigidFragment::distance_squared(pele::VecN<3> const & com1, pele::VecN<3> const & p1,
00080         pele::VecN<3> const & com2, pele::VecN<3> const & p2) const
00081 {
00082     VecN<3> drij = get_smallest_rij(com1, com2);
00083     pele::MatrixNM<3,3> R1 = pele::aa_to_rot_mat(p1);
00084     pele::MatrixNM<3,3> R2 = pele::aa_to_rot_mat(p2);
00085 
00086     MatrixNM<3,3> dR = R2 - R1;  
00087 
00088     double d_M = m_W * dot<3>(drij, drij);
00089     // we only need the trace, so this can be sped up
00090     double d_P = dot<3,3,3>(dR, dot<3,3,3>(m_S, transpose<3>(dR))).trace();
00091     double d_mix = 2. * m_W * dot<3>(drij, dot<3,3>(dR, m_cog));
00092 
00093     double dist2 = d_M + d_P + d_mix;
00094     return dist2;
00095 }
00096 
00097 void
00098 pele::RigidFragment::distance_squared_grad(pele::VecN<3> const & com1, pele::VecN<3> const & p1,
00099         pele::VecN<3> const & com2, pele::VecN<3> const & p2,
00100         VecN<3> & g_M, VecN<3> & g_P
00101         ) const
00102 {
00103     VecN<3> drij = get_smallest_rij(com1, com2);
00104     auto R2 = pele::aa_to_rot_mat(p2);
00105     MatrixNM<3,3> R1, R11, R12, R13;
00106     pele::rot_mat_derivatives(p1, R1, R11, R12, R13);
00107 
00108     auto dR = R2 - R1;
00109 
00110     g_M = drij;
00111     g_M *= -2. * m_W;
00112 
00113     // this linear algebra can be done more efficiently
00114     auto dRT = pele::transpose(dR);
00115     g_P[0] = -2. * dot<3,3,3>(R11, dot<3,3,3>(m_S, dRT)).trace();
00116     g_P[1] = -2. * dot<3,3,3>(R12, dot<3,3,3>(m_S, dRT)).trace();
00117     g_P[2] = -2. * dot<3,3,3>(R13, dot<3,3,3>(m_S, dRT)).trace();
00118 
00119     // this can also be done more efficiently
00120     auto temp = dot<3,3>(dR, m_cog);
00121     temp *= 2. * m_W;
00122     g_M -= temp;
00123     g_P[0] -= 2. * m_W * dot<3>(drij, dot<3,3>(R11, m_cog));
00124     g_P[1] -= 2. * m_W * dot<3>(drij, dot<3,3>(R12, m_cog));
00125     g_P[2] -= 2. * m_W * dot<3>(drij, dot<3,3>(R13, m_cog));
00126 }
00127 
00128 void pele::MeasureAngleAxisCluster::
00129 align(pele::Array<double> const x1, pele::Array<double> x2)
00130 {
00131     auto c1 = m_topology->get_coords_adaptor(x1);
00132     auto c2 = m_topology->get_coords_adaptor(x2);
00133 
00134     // now account for the symmetries
00135     for (size_t isite = 0; isite < m_topology->nrigid(); ++isite) {
00136         auto const & rotations = m_topology->get_sites()[isite].get_symmetry_rotations();
00137         auto p1 = c1.get_rb_rotation(isite);
00138         auto p2 = c2.get_rb_rotation(isite);
00139         auto mx2 = pele::aa_to_rot_mat(p2);
00140         auto mx1 = pele::transpose(pele::aa_to_rot_mat(p1));
00141         auto mx =  pele::dot(mx1, mx2);
00142         double theta_min = 10.;
00143         MatrixNM<3,3> rot_best = pele::identity<3>();
00144         for (auto const & rot : rotations){
00145             auto mx_diff = dot(mx, rot);
00146             double theta = norm<3>(rot_mat_to_aa(mx_diff));
00147             theta -= int(theta / (2. * M_PI)) * 2. * M_PI;
00148             if (theta < theta_min) {
00149                 theta_min = theta;
00150                 rot_best = rot;
00151             }
00152         }
00153         auto newp2 = rotate_aa(rot_mat_to_aa(rot_best), p2);
00154         std::copy(newp2.begin(), newp2.end(), p2.begin());
00155     }
00156 }
00157 
00158 Array<double>
00159 pele::RBTopology::to_atomistic(Array<double> rbcoords)
00160 {
00161     if ( rbcoords.size() != nrigid() * 6 ) {
00162         throw std::invalid_argument("rbcoords has the wrong size");
00163     }
00164 
00165     size_t const nrigid = _sites.size();
00166     auto ca = get_coords_adaptor(rbcoords);
00167     Array<double> atomistic(3 * natoms_total());
00168     // view the atomistic coords as a matrix
00169     size_t istart = 0;
00170     for (size_t isite=0; isite<nrigid; ++isite) {
00171         VecN<3> psite = ca.get_rb_rotation(isite);
00172         auto site_atom_positions = _sites[isite].to_atomistic(
00173                 ca.get_rb_position(isite),
00174                 psite);
00175         Array<double> atomistic_view(atomistic.view(istart, istart + site_atom_positions.size()));
00176         atomistic_view.assign(site_atom_positions);
00177 
00178         istart += site_atom_positions.size();
00179     }
00180     assert(istart == natoms_total() * 3);
00181     return atomistic;
00182 }
00183 
00184 void
00185 pele::RBTopology::transform_gradient(pele::Array<double> rbcoords,
00186         pele::Array<double> grad, pele::Array<double> rbgrad)
00187 {
00188     if ( rbcoords.size() != nrigid() * 6 ) {
00189         throw std::invalid_argument("rbcoords has the wrong size");
00190     }
00191     if (grad.size() != natoms_total() * 3) {
00192         throw std::invalid_argument("grad has the wrong size");
00193     }
00194     if (rbgrad.size() != rbcoords.size()) {
00195         throw std::invalid_argument("rbgrad has the wrong size");
00196     }
00197 
00198     CoordsAdaptor ca(nrigid(), 0, rbcoords);
00199     pele::Array<double> coords_rot(ca.get_rb_rotations());
00200 //        pele::Array<double> rbgrad(rbcoords.size());
00201     CoordsAdaptor rbgrad_ca(nrigid(), 0, rbgrad);
00202 
00203     size_t istart = 0;
00204     for (size_t isite=0; isite<nrigid(); ++isite) {
00205         size_t const site_ndof = _sites[isite].natoms() * 3;
00206 //            std::cout << grad.size() << " " << istart << " " << site_ndof << " " << istart + site_ndof << "\n";
00207         Array<double> g_site     = grad.view(istart, istart + site_ndof);
00208         Array<double> p          = ca.get_rb_rotation(isite);
00209         Array<double> g_com_site = rbgrad_ca.get_rb_position(isite);
00210         Array<double> g_rot_site = rbgrad_ca.get_rb_rotation(isite);
00211         _sites[isite].transform_grad(
00212                 ca.get_rb_rotation(isite),
00213                 grad.view(istart, istart + site_ndof),
00214                 rbgrad_ca.get_rb_position(isite),
00215                 rbgrad_ca.get_rb_rotation(isite));
00216 //                p, g_site, g_com_site, g_rot_site);
00217         istart += site_ndof;
00218     }
00219 }
00220 
00221 pele::VecN<3>
00222 pele::RBTopology::align_angle_axis_vectors(pele::VecN<3> const & p1,
00223         pele::VecN<3> const & p2in)
00224 {
00225     pele::VecN<3> p2 = p2in;
00226     pele::VecN<3> n2, p2n;
00227     if (norm<3>(p2) < 1e-6) {
00228         if (norm<3>(p1) < 1e-6) {
00229             return p2;
00230         }
00231         n2 = p1;
00232         n2 *= 2. * M_PI / norm<3>(p1);
00233     } else {
00234         n2 = p2;
00235         n2 *= 2. * M_PI / norm<3>(p2);
00236     }
00237 
00238     while (true) {
00239         p2n = p2;
00240         p2n += n2;
00241         if (norm<3>(p2n - p1) > norm<3>(p2 - p1)) {
00242             break;
00243         }
00244         p2 = p2n;
00245     }
00246 
00247     while (true) {
00248         p2n = p2;
00249         p2n -= n2;
00250         if (norm<3>(p2n - p1) > norm<3>(p2 - p1)) {
00251             break;
00252         }
00253         p2 = p2n;
00254     }
00255     return p2;
00256 }
00257 
00258 void
00259 pele::RBTopology::align_all_angle_axis_vectors(pele::Array<double> x1,
00260         pele::Array<double> x2)
00261 {
00262     auto c1 = get_coords_adaptor(x1);
00263     auto c2 = get_coords_adaptor(x2);
00264     for (size_t isite = 0; isite < nrigid(); ++isite) {
00265         VecN<3> p1 = c1.get_rb_rotation(isite);
00266         pele::Array<double> p2 = c2.get_rb_rotation(isite);
00267         auto p2new = align_angle_axis_vectors(p1, p2);
00268         std::copy(p2new.begin(), p2new.end(), p2.begin());
00269     }
00270 }
00271 
00272 void
00273 pele::RBTopology::align_path(std::list<pele::Array<double> > path)
00274 {
00275     auto iter1 = path.begin();
00276     auto iter2 = path.begin();
00277     iter2++;
00278     while (iter2 != path.end()) {
00279         align_all_angle_axis_vectors(*iter1, *iter2);
00280         ++iter1;
00281         ++iter2;
00282     }
00283 }
00284 
00285 void
00286 pele::TransformAACluster::rotate(pele::Array<double> x,
00287         pele::MatrixNM<3,3> const & mx)
00288 {
00289     auto ca = m_topology->get_coords_adaptor(x);
00290     if(m_topology->nrigid() > 0) {
00291         // rotate the center of mass positions by mx
00292         pele::MatrixAdapter<double> rb_pos(ca.get_rb_positions(), 3);
00293         // make a MatrixAdapter view of the transposed rotation matrix
00294         auto mxT = pele::transpose(mx);
00295         pele::MatrixAdapter<double> mxT_view(mxT.data(), 3, 3);
00296         assert(mxT_view(0,1) == mx(1,0));
00297         // do the multiplication
00298         auto result = hacky_mat_mul(rb_pos, mxT_view);
00299         // copy the results back into the coordinates array
00300 //        std::cout << "result " << result << std::endl;
00301         rb_pos.assign(result);
00302 
00303         // rotate each aa rotation by mx
00304         VecN<3> dp = pele::rot_mat_to_aa(mx);
00305         for (size_t isite = 0; isite < m_topology->nrigid(); ++isite) {
00306             pele::Array<double> pview = ca.get_rb_rotation(isite);
00307             VecN<3> p = pele::rotate_aa(pview, dp);
00308             // copy the vector back into pview
00309             std::copy(p.begin(), p.end(), pview.begin());
00310         }
00311     }
00312     if (m_topology->number_of_non_rigid_atoms() > 0) {
00313         throw std::runtime_error("non-rigid atoms is not yet supported");
00314 //            ca.posAtom[:] = np.dot(mx, ca.posAtom.transpose()).transpose()
00315     }
00316 }
00317 
00318 } // namespace pele
 All Classes Namespaces Functions Variables Typedefs