libhoomd/computes/PPPMForceCompute.h

   1 /*
   2 Highly Optimized Object-oriented Many-particle Dynamics -- Blue Edition
   3 (HOOMD-blue) Open Source Software License Copyright 2009-2014 The Regents of
   4 the University of Michigan All rights reserved.
   5
   6 HOOMD-blue may contain modifications ("Contributions") provided, and to which
   7 copyright is held, by various Contributors who have granted The Regents of the
   8 University of Michigan the right to modify and/or distribute such Contributions.
   9
  10 You may redistribute, use, and create derivate works of HOOMD-blue, in source
  11 and binary forms, provided you abide by the following conditions:
  12
  13 * Redistributions of source code must retain the above copyright notice, this
  14 list of conditions, and the following disclaimer both in the code and
  15 prominently in any materials provided with the distribution.
  16
  17 * Redistributions in binary form must reproduce the above copyright notice, this
  18 list of conditions, and the following disclaimer in the documentation and/or
  19 other materials provided with the distribution.
  20
  21 * All publications and presentations based on HOOMD-blue, including any reports
  22 or published results obtained, in whole or in part, with HOOMD-blue, will
  23 acknowledge its use according to the terms posted at the time of submission on:
  24 http://codeblue.umich.edu/hoomd-blue/citations.html
  25
  26 * Any electronic documents citing HOOMD-Blue will link to the HOOMD-Blue website:
  27 http://codeblue.umich.edu/hoomd-blue/
  28
  29 * Apart from the above required attributions, neither the name of the copyright
  30 holder nor the names of HOOMD-blue's contributors may be used to endorse or
  31 promote products derived from this software without specific prior written
  32 permission.
  33
  34 Disclaimer
  35
  36 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS ``AS IS'' AND
  37 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  38 WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, AND/OR ANY
  39 WARRANTIES THAT THIS SOFTWARE IS FREE OF INFRINGEMENT ARE DISCLAIMED.
  40
  41 IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
  42 INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  43 BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  44 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  45 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
  46 OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
  47 ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  48 */
  49
  50 // Maintainer: sbarr
  51
  52 #include <boost/shared_ptr.hpp>
  53 #include <boost/signals2.hpp>
  54
  55 #include "ForceCompute.h"
  56 #include "NeighborList.h"
  57 #include "ParticleGroup.h"
  58
  59 #include <vector>
  60
  61 #ifdef ENABLE_CUDA
  62 #include <cufft.h>
  63 #endif
  64
  65 /*! \file PPPMForceCompute.h
  66     \brief Declares a class for computing harmonic bonds
  67 */
  68
  69 #ifdef NVCC
  70 #error This header cannot be compiled by nvcc
  71 #endif
  72
  73 // CUFFTCOMPLEX is cufftDoubleComplex when run in double precision with CUDA enabled, and cufftComplex otherwise
  74 #ifndef SINGLE_PRECISION
  75 #define CUFFTCOMPLEX cufftDoubleComplex
  76 #else
  77 #define CUFFTCOMPLEX cufftComplex
  78 #endif
  79
  80 #ifndef __PPPMFORCECOMPUTE_H__
  81 #define __PPPMFORCECOMPUTE_H__
  82
  83 // slave KISS data type to HOOMD Scalar
  84 #ifndef kiss_fft_scalar
  85 #define kiss_fft_scalar Scalar
  86 #endif
  87 #include "HOOMDMath.h"
  88 #include "kiss_fftnd.h"
  89
  90
  91 // MAX gives the larger of two values
  92 #define MAX(a,b) ((a) > (b) ? (a) : (b))
  93 // MIN gives the lesser of two values
  94 #define MIN(a,b) ((a) < (b) ? (a) : (b))
  95 // MaxOrder is the largest allowed value of the interpolation order
  96 #define MaxOrder 7
  97 // ConstSize is used to make sure the rho_coeff will fit into memory on the GPU
  98 #define CONSTANT_SIZE 2048
  99 // EPS_HOC is used to calculate the Green's function
 100 #define EPS_HOC 1.0e-7
 101
 102
 103 //! Computes the long ranged part of the electrostatic forces on each particle
 104 /*! PPPM forces are computed on every particle in the simulation.
 105
 106 */
 107 class PPPMForceCompute : public ForceCompute
 108     {
 109     public:
 110         //! Constructs the compute
 111         PPPMForceCompute(boost::shared_ptr<SystemDefinition> sysdef,
 112                          boost::shared_ptr<NeighborList> nlist,
 113                          boost::shared_ptr<ParticleGroup> group);
 114
 115         //! Destructor
 116         ~PPPMForceCompute();
 117
 118         //! Set the parameters
 119         virtual void setParams(int Nx, int Ny, int Nz, int order, Scalar kappa, Scalar rcut);
 120
 121         //! Returns a list of log quantities this compute calculates
 122         virtual std::vector< std::string > getProvidedLogQuantities();
 123
 124         //! Calculates the requested log value and returns it
 125         virtual Scalar getLogValue(const std::string& quantity, unsigned int timestep);
 126
 127         //! Notification of a box size change
 128         void slotBoxChanged()
 129             {
 130             m_box_changed = true;
 131             }
 132
 133         //! root mean square error in force calculation
 134         Scalar rms(Scalar h, Scalar prd, Scalar natoms);
 135         //! computes coefficients for assigning charges to grid points
 136         void compute_rho_coeff();
 137         //! computes coefficients for the Green's function
 138         void compute_gf_denom();
 139         //! computes coefficients for the Green's function
 140         Scalar gf_denom(Scalar x, Scalar y, Scalar z);
 141         //! resets kvec, Green's function and virial coefficients if the box size changes
 142         void reset_kvec_green_hat_cpu();
 143         //! assigns charges to grid points
 144         void assign_charges_to_grid();
 145         //! multiply Green's function by charge density to get electric field
 146         void combined_green_e();
 147         //! Do the final force calculation
 148         void calculate_forces();
 149         //! fix the force due to excluded particles
 150         void fix_exclusions_cpu();
 151         //! fix the energy and virial thermodynamic quantities
 152         virtual void fix_thermo_quantities();
 153
 154     protected:
 155         GPUArray<Scalar>m_vg;                    //!< Virial coefficient
 156         Scalar m_thermo_data[7];                 //!< PPPM contribution to energy and virial
 157         bool m_params_set;                       //!< Set to true when the parameters are set
 158         int m_Nx;                                //!< Number of grid points in x direction
 159         int m_Ny;                                //!< Number of grid points in y direction
 160         int m_Nz;                                //!< Number of grid points in z direction
 161         int m_order;                             //!< Interpolation order
 162         Scalar m_kappa;                          //!< screening parameter for erfc(kappa*r)
 163         Scalar m_rcut;                           //!< Real space cutoff
 164         Scalar m_q;                              //!< Total system charge
 165         Scalar m_q2;                             //!< Sum(q_i*q_i), where q_i is the charge of each particle
 166         Scalar m_energy_virial_factor;           //!< Multiplication factor for energy and virial
 167         bool m_box_changed;                      //!< Set to ttrue when the box size has changed
 168         GPUArray<Scalar3> m_kvec;                //!< k-vectors for each grid point
 169         GPUArray<CUFFTCOMPLEX> m_rho_real_space; //!< x component of the grid based electric field
 170         GPUArray<CUFFTCOMPLEX> m_Ex;             //!< x component of the grid based electric field
 171         GPUArray<CUFFTCOMPLEX> m_Ey;             //!< y component of the grid based electric field
 172         GPUArray<CUFFTCOMPLEX> m_Ez;             //!< z component of the grid based electric field
 173         GPUArray<Scalar3>m_field;                //!< grid based Electric field, combined
 174         GPUArray<Scalar> m_rho_coeff;            //!< Coefficients for computing the grid based charge density
 175         GPUArray<Scalar> m_gf_b;                 //!< Used to compute the grid based Green's function
 176         GPUArray<Scalar> m_green_hat;            //!< Modified Hockney-Eastwood Green's function
 177         GPUArray<Scalar> o_data;                 //!< Used to quickly sum grid points for pressure and energy calcuation (output)
 178         GPUArray<Scalar> m_energy_sum;           //!< Used to quickly sum grid points for pressure and energy calcuation (input)
 179         GPUArray<Scalar> m_v_xx_sum;             //!< Used to quickoy sum grid points for virial_xx
 180         GPUArray<Scalar> m_v_xy_sum;             //!< Used to quickoy sum grid points for virial_xy
 181         GPUArray<Scalar> m_v_xz_sum;             //!< Used to quickoy sum grid points for virial_xz
 182         GPUArray<Scalar> m_v_yy_sum;             //!< Used to quickoy sum grid points for virial_yy
 183         GPUArray<Scalar> m_v_yz_sum;             //!< Used to quickoy sum grid points for virial_yz
 184         GPUArray<Scalar> m_v_zz_sum;             //!< Used to quickoy sum grid points for virial_zz
 185         boost::signals2::connection m_boxchange_connection;   //!< Connection to the ParticleData box size change signal
 186         boost::shared_ptr<NeighborList> m_nlist; //!< The neighborlist to use for the computation
 187         boost::shared_ptr<ParticleGroup> m_group;//!< Group to compute properties for
 188         kiss_fft_cpx *fft_in;                    //!< For FFTs on CPU rho_real_space
 189         kiss_fft_cpx *fft_ex;                    //!< For FFTs on CPU E-field x component
 190         kiss_fft_cpx *fft_ey;                    //!< For FFTs on CPU E-field y component
 191         kiss_fft_cpx *fft_ez;                    //!< For FFTs on CPU E-field z component
 192         kiss_fftnd_cfg fft_forward;              //!< Forward FFT on CPU
 193         kiss_fftnd_cfg fft_inverse;              //!< Inverse FFT on CPU
 194         int first_run;                           //!< flag for allocating arrays
 195
 196         //! Actually compute the forces
 197         virtual void computeForces(unsigned int timestep);
 198     };
 199
 200
 201 //! Exports the PPPMForceCompute class to python
 202 void export_PPPMForceCompute();
 203
 204
 205
 206 #endif