pele
Python energy landscape explorer
/home/js850/projects/pele/source/pele/atomlist_potential.h
00001 #ifndef _ATOMLIST_POTENTIAL_H_
00002 #define _ATOMLIST_POTENTIAL_H_
00003 
00004 #include "array.h"
00005 #include "base_potential.h"
00006 #include <iostream>
00007 
00008 namespace pele {
00009 
00010 template<typename pairwise_interaction, typename distance_policy>
00011 class AtomListPotential : public BasePotential
00012 {
00013 protected:
00014     std::shared_ptr<pairwise_interaction> _interaction;
00015     std::shared_ptr<distance_policy> _dist;
00016     Array<size_t> _atoms1;
00017     Array<size_t> _atoms2;
00018     bool _one_list;
00019     static const size_t _ndim = distance_policy::_ndim;
00020 
00021     AtomListPotential(
00022             std::shared_ptr<pairwise_interaction> interaction,
00023             std::shared_ptr<distance_policy> dist,
00024             Array<size_t> & atoms1, Array<size_t> & atoms2) 
00025         : _interaction(interaction),
00026           _dist(dist),
00027           _atoms1(atoms1.copy()),
00028           _atoms2(atoms2.copy()),
00029           _one_list(false)
00030     {}
00031 
00032     AtomListPotential( std::shared_ptr<pairwise_interaction> interaction,
00033             std::shared_ptr<distance_policy> dist, Array<size_t> & atoms1) 
00034         : _interaction(interaction),
00035           _dist(dist),
00036           _atoms1(atoms1.copy()),
00037           _atoms2(_atoms1),
00038           _one_list(true)
00039     {}
00040 
00041 
00042 public:
00043 
00044     virtual inline double get_energy(Array<double> x)
00045     {
00046         double e=0.;
00047         size_t jstart = 0;
00048         double dr[_ndim];
00049 
00050         for(size_t i=0; i<_atoms1.size(); ++i) {
00051             size_t atom1 = _atoms1[i];
00052             size_t i1 = _ndim * atom1;
00053             if (_one_list)
00054                 jstart = i+1;
00055             for(size_t j=jstart; j<_atoms2.size(); ++j) {
00056                 size_t atom2 = _atoms2[j];
00057                 size_t i2 = _ndim * atom2;
00058 
00059                 _dist->get_rij(dr, &x[i1], &x[i2]);
00060                 double r2 = 0;
00061                 for (size_t k=0;k<_ndim;++k){r2 += dr[k]*dr[k];}
00062                 e += _interaction->energy(r2, atom1, atom2);
00063             }
00064         }
00065 
00066         return e;
00067 
00068     }
00069 
00070     virtual inline double add_energy_gradient(Array<double> x, Array<double> grad)
00071     {
00072         if (x.size() != grad.size()) {
00073             throw std::invalid_argument("the gradient has the wrong size");
00074         }
00075         double e=0.;
00076         double gij;
00077         size_t jstart = 0;
00078         double dr[_ndim];;
00079 
00080         for(size_t i=0; i<_atoms1.size(); ++i) {
00081             size_t atom1 = _atoms1[i];
00082             size_t i1 = _ndim * atom1;
00083             if (_one_list){
00084                 jstart = i+1;
00085             }
00086             for(size_t j=jstart; j<_atoms2.size(); ++j) {
00087                 size_t atom2 = _atoms2[j];
00088                 size_t i2 = _ndim * atom2;
00089 
00090                 _dist->get_rij(dr, &x[i1], &x[i2]);
00091 
00092                 double r2 = 0;
00093                 for (size_t k=0;k<_ndim;++k){r2 += dr[k]*dr[k];}
00094 
00095                 e += _interaction->energy_gradient(r2, &gij, atom1, atom2);
00096                 for(size_t k=0; k<_ndim; ++k)
00097                     grad[i1+k] -= gij * dr[k];
00098                 for(size_t k=0; k<_ndim; ++k)
00099                     grad[i2+k] += gij * dr[k];
00100             }
00101         }
00102 
00103         return e;
00104     }
00105 
00106     virtual inline double add_energy_gradient_hessian(Array<double> x,
00107             Array<double> grad, Array<double> hess)
00108     {
00109         if (x.size() != grad.size()) {
00110             throw std::invalid_argument("the gradient has the wrong size");
00111         }
00112         if (hess.size() != x.size() * x.size()) {
00113             throw std::invalid_argument("the Hessian has the wrong size");
00114         }
00115 
00116         double e=0.;
00117         double hij, gij;
00118         size_t jstart = 0;
00119         double dr[_ndim];
00120         const size_t N = x.size();
00121 
00122         for(size_t i=0; i<_atoms1.size(); ++i) {
00123             size_t atom1 = _atoms1[i];
00124             size_t i1 = _ndim * atom1;
00125 
00126             if (_one_list){
00127                 jstart = i+1;
00128             }
00129             for(size_t j=jstart; j<_atoms2.size(); ++j) {
00130                 size_t atom2 = _atoms2[j];
00131                 size_t i2 = _ndim * atom2;
00132 
00133                 _dist->get_rij(dr, &x[i1], &x[i2]);
00134                 double r2 = 0;
00135                 for (size_t k=0;k<_ndim;++k){r2 += dr[k]*dr[k];}
00136 
00137                 e += _interaction->energy_gradient_hessian(r2, &gij, &hij, atom1, atom2);
00138                 for(size_t k=0; k<_ndim; ++k)
00139                     grad[i1+k] -= gij * dr[k];
00140                 for(size_t k=0; k<_ndim; ++k)
00141                     grad[i2+k] += gij * dr[k];
00142 
00143 
00144                 for (size_t k=0; k<_ndim; ++k){
00145                     //diagonal block - diagonal terms
00146                     double Hii_diag = (hij+gij)*dr[k]*dr[k]/r2 - gij;
00147                     hess[N*(i1+k)+i1+k] += Hii_diag;
00148                     hess[N*(i2+k)+i2+k] += Hii_diag;
00149                     //off diagonal block - diagonal terms
00150                     double Hij_diag = -Hii_diag;
00151                     hess[N*(i1+k)+i2+k] = Hij_diag;
00152                     hess[N*(i2+k)+i1+k] = Hij_diag;
00153                     for (size_t l = k+1; l<_ndim; ++l){
00154                         //diagonal block - off diagonal terms
00155                         double Hii_off = (hij+gij)*dr[k]*dr[l]/r2;
00156                         hess[N*(i1+k)+i1+l] += Hii_off;
00157                         hess[N*(i1+l)+i1+k] += Hii_off;
00158                         hess[N*(i2+k)+i2+l] += Hii_off;
00159                         hess[N*(i2+l)+i2+k] += Hii_off;
00160                         //off diagonal block - off diagonal terms
00161                         double Hij_off = -Hii_off;
00162                         hess[N*(i1+k)+i2+l] = Hij_off;
00163                         hess[N*(i1+l)+i2+k] = Hij_off;
00164                         hess[N*(i2+k)+i1+l] = Hij_off;
00165                         hess[N*(i2+l)+i1+k] = Hij_off;
00166                     }
00167                 }
00168             }
00169         }
00170 
00171         return e;
00172     }
00173 
00174 };
00175 }
00176 
00177 #endif
 All Classes Namespaces Functions Variables Typedefs