MADNESS  0.10.1
poperator.h
Go to the documentation of this file.
1 /*
2  This file is part of MADNESS.
3 
4  Copyright (C) 2007,2010 Oak Ridge National Laboratory
5 
6  This program is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2 of the License, or
9  (at your option) any later version.
10 
11  This program is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with this program; if not, write to the Free Software
18  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20  For more information please contact:
21 
22  Robert J. Harrison
23  Oak Ridge National Laboratory
24  One Bethel Valley Road
25  P.O. Box 2008, MS-6367
26 
27  email: harrisonrj@ornl.gov
28  tel: 865-241-3937
29  fax: 865-572-0680
30 
31  $Id$
32 */
33 #ifndef POPERATOR_H_
34 #define POPERATOR_H_
35 
36 #include <madness/constants.h>
37 #include <madness/mra/operator.h>
38 
39 #define WST_PI madness::constants::pi
40 //#define WST_PI 3.14159265358979323846264338328
41 
42 namespace madness
43 {
44  //***************************************************************************
45  const double acut1e_6 = 0.25; //0.6450626287524907;
46  //***************************************************************************
47 
48 // //***************************************************************************
49 // template<typename Q, int NDIM>
50 // SeparatedConvolution<Q, NDIM> PeriodicHFExchangeOp(World& world, long k,
51 // double lo, double eps, Tensor<double> L, Vector<double,3> q)
52 // {
53 // // bsh_fit generates representation for 1/4Pir but we want 1/r
54 // // so have to scale eps by 1/4Pi
55 // Tensor<double> coeff, expnt;
56 // //if (mu==0) eps /= 4.0*pi;
57 // bsh_fit(0.0, lo, 100.0*L[0], eps/(4.0 * WST_PI), &coeff, &expnt, false); //eps /(4.0*pi)
58 //
59 // // Scale coefficients according to the dimensionality and add to the list of operators
60 // std::vector< std::shared_ptr< Convolution1D<Q> > > ops;
61 // for (int i=0; i < coeff.dim(0); ++i)
62 // {
63 // if (expnt[i]*L[0]*L[0] > acut1e_6)
64 // {
65 // double c = pow(4 * WST_PI * coeff[i], 1.0/double(NDIM));
66 // ops.push_back(std::shared_ptr< Convolution1D<Q> >(new GaussianConvolution1D<Q>(k,
67 // c*L[0], expnt[i]*L[0]*L[0], 1.0, 0, true, q)));
68 // }
69 // }
70 //
71 // return SeparatedConvolution<Q, NDIM>(world, ops);
72 // }
73 // //***************************************************************************
74 
75  //***************************************************************************
76  template<typename Q, int NDIM>
78  {
79  // bsh_fit generates representation for 1/4Pir but we want 1/r
80  // so have to scale eps by 1/4Pi
81  Tensor<double> coeff, expnt;
82  //if (mu==0) eps /= 4.0*pi;
83  bsh_fit(0.0, lo, 100.0*L[0], eps/(4.0 * WST_PI), &coeff, &expnt, false); //eps /(4.0*pi)
84 
85  // Scale coefficients according to the dimensionality and add to the list of operators
86  std::vector< std::shared_ptr< Convolution1D<Q> > > ops;
87  for (int i=0; i < coeff.dim(0); ++i)
88  {
89  if (expnt[i]*L[0]*L[0] > acut1e_6)
90  {
91  double c = pow(4 * WST_PI * coeff[i], 1.0/double(NDIM));
92  ops.push_back(std::shared_ptr< Convolution1D<Q> >(new GaussianConvolution1D<Q>(k, c*L[0], expnt[i]*L[0]*L[0], 1.0, 0, true)));
93  }
94  }
95 
96  return SeparatedConvolution<Q, NDIM>(world, ops);
97  }
98  //***************************************************************************
99 
100  //***************************************************************************
101  template<typename Q, int NDIM>
103  {
104  // bsh_fit generates representation for 1/4Pir but we want 1/r
105  // so have to scale eps by 1/4Pi
106  Tensor<double> coeff, expnt;
107  //if (mu==0) eps /= 4.0*pi;
108  bsh_fit(0.0, lo, 100.0*L[0], eps/(4.0 * WST_PI), &coeff, &expnt, false); //eps /(4.0*pi)
109 
110  // Scale coefficients according to the dimensionality and add to the list of operators
111  std::vector< std::shared_ptr< Convolution1D<Q> > > ops;
112  for (int i=0; i < coeff.dim(0); ++i)
113  {
114  if (expnt[i]*L[0]*L[0] > acut1e_6)
115  {
116  double c = pow(4 * WST_PI * coeff[i], 1.0/double(NDIM));
117  ops.push_back(std::shared_ptr< Convolution1D<Q> >(new GaussianConvolution1D<double>(k, c*L[0], expnt[i]*L[0]*L[0], 1.0, 0, true)));
118  }
119  }
120 
121  return new SeparatedConvolution<Q, NDIM>(world, ops);
122  }
123  //***************************************************************************
124 
125  //***************************************************************************
126  template<typename Q, int NDIM>
127  SeparatedConvolution<Q, NDIM> PeriodicBSHOp(World& world, double mu, long k, double lo, double eps, Tensor<double> L)
128  {
129  // bsh_fit generates representation for 1/4Pir but we want 1/r
130  // so have to scale eps by 1/4Pi
131  Tensor<double> coeff, expnt;
132  //if (mu==0) eps /= 4.0*pi;
133  bsh_fit(mu, lo, 10.0*L[0], eps, &coeff, &expnt, false); //eps /(4.0*pi)
134 
135  // Scale coefficients according to the dimensionality and add to the list of operators
136  std::vector< std::shared_ptr< Convolution1D<Q> > > ops;
137  for (int i=0; (i < coeff.dim(0)); ++i)
138  {
139  double c = pow(coeff[i], 1.0/double(NDIM));
140  ops.push_back(std::shared_ptr< Convolution1D<Q> >(new GaussianConvolution1D<double>(k, c*L[0], expnt[i]*L[0]*L[0], 1.0, 0, true)));
141  }
142 
143  return SeparatedConvolution<Q, NDIM>(world, ops);
144  }
145  //***************************************************************************
146 
147  //***************************************************************************
148  template<typename Q, int NDIM>
149  SeparatedConvolution<Q, NDIM>* PeriodicBSHOpPtr(World& world, double mu, long k, double lo, double eps, Tensor<double> L)
150  {
151  // bsh_fit generates representation for 1/4Pir but we want 1/r
152  // so have to scale eps by 1/4Pi
153  Tensor<double> coeff, expnt;
154  //if (mu==0) eps /= 4.0*pi;
155  bsh_fit(mu, lo, 10.0*L[0], eps, &coeff, &expnt, false); //eps /(4.0*pi)
156 
157  // Scale coefficients according to the dimensionality and add to the list of operators
158  std::vector< std::shared_ptr< Convolution1D<Q> > > ops;
159  for (int i=0; (i < coeff.dim(0)); ++i)
160  {
161  double c = pow(coeff[i], 1.0/double(NDIM));
162  ops.push_back(std::shared_ptr< Convolution1D<Q> >(new GaussianConvolution1D<double>(k, c*L[0], expnt[i]*L[0]*L[0], 1.0, 0, true)));
163  }
164 
165  return new SeparatedConvolution<Q, NDIM>(world, ops);
166  }
167  //***************************************************************************
168 
169 };
170 
171 #endif
long dim(int i) const
Returns the size of dimension i.
Definition: basetensor.h:147
Provides the common functionality/interface of all 1D convolutions.
Definition: convolution1d.h:257
1D convolution with (derivative) Gaussian; coeff and expnt given in simulation coordinates [0,...
Definition: convolution1d.h:683
Convolutions in separated form (including Gaussian)
Definition: operator.h:136
A parallel world class.
Definition: world.h:132
Defines common mathematical and physical constants.
static double lo
Definition: dirac-hatom.cc:23
static double pow(const double *a, const double *b)
Definition: lda.h:74
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
SeparatedConvolution< Q, NDIM > PeriodicCoulombOp(World &world, long k, double lo, double eps, Tensor< double > L)
Definition: poperator.h:77
SeparatedConvolution< Q, NDIM > PeriodicBSHOp(World &world, double mu, long k, double lo, double eps, Tensor< double > L)
Definition: poperator.h:127
SeparatedConvolution< Q, NDIM > * PeriodicBSHOpPtr(World &world, double mu, long k, double lo, double eps, Tensor< double > L)
Definition: poperator.h:149
const double acut1e_6
Definition: poperator.h:45
SeparatedConvolution< Q, NDIM > * PeriodicCoulombOpPtr(World &world, long k, double lo, double eps, Tensor< double > L)
Definition: poperator.h:102
const double mu
Definition: navstokes_cosines.cc:95
Implements most functionality of separated operators.
#define WST_PI
Definition: poperator.h:39
static const double c
Definition: relops.cc:10
static const double L
Definition: rk.cc:46
static const long k
Definition: rk.cc:44
static const std::size_t NDIM
Definition: testpdiff.cc:42