MADNESS  0.10.1
potentialmanager.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 MADNESS_CHEM_POTENTIALMANAGER_H__INCLUDED
34 #define MADNESS_CHEM_POTENTIALMANAGER_H__INCLUDED
35 
36 /// \file moldft/potentialmanager.h
37 /// \brief Declaration of molecule related classes and functions
38 
42 #include <vector>
43 #include <string>
44 #include <iostream>
45 #include <fstream>
46 #include <sstream>
47 #include <algorithm>
48 #include <ctype.h>
49 #include <cmath>
50 #include <madness/tensor/tensor.h>
51 #include <madness/misc/misc.h>
52 #include <madness/mra/mra.h>
53 
54 namespace madness {
56 private:
58 public:
60  : molecule(molecule) {}
61 
62  double operator()(const coord_3d& x) const {
63  return molecule.nuclear_attraction_potential(x[0], x[1], x[2]);
64  }
65 
66  std::vector<coord_3d> special_points() const {return molecule.get_all_coords_vec();}
67 };
68 
70 private:
72 public:
74  : molecule(molecule) {}
75 
76  double operator()(const coord_3d& x) const {
77  return molecule.molecular_core_potential(x[0], x[1], x[2]);
78  }
79 
80  std::vector<coord_3d> special_points() const {return molecule.get_all_coords_vec();}
81 };
82 
83 class CoreOrbitalFunctor : public FunctionFunctorInterface<double,3> {
85  const int atom;
86  const unsigned int core;
87  const int m;
88 public:
89  CoreOrbitalFunctor(Molecule& molecule, int atom, unsigned int core, int m)
90  : molecule(molecule), atom(atom), core(core), m(m) {};
91  double operator()(const coord_3d& r) const {
92  return molecule.core_eval(atom, core, m, r[0], r[1], r[2]);
93  };
94 };
95 
98  const int atom, axis;
99  const unsigned int core;
100  const int m;
101 public:
103  : molecule(molecule), atom(atom), axis(axis), core(core), m(m) {};
104  double operator()(const coord_3d& r) const {
105  return molecule.core_derivative(atom, axis, core, m, r[0], r[1], r[2]);
106  };
107 };
108 
109 
111 private:
114 std::string core_type_;
115 
116 public:
117  PotentialManager(const Molecule& molecule, const std::string& core_type)
119 
120  const Molecule& molecule() const {
121  return this->mol;
122  }
123 
124  const std::string& core_type() const {
125  return this->core_type_;
126  }
127 
129  return vnuc;
130  }
131 
132  vector_real_function_3d core_projection(World & world, const vector_real_function_3d& psi, const bool include_Bc = true)
133  {
134  int npsi = psi.size();
135  if (npsi == 0) return psi;
136  int natom = mol.natom();
137  vector_real_function_3d proj = zero_functions_compressed<double,3>(world, npsi);
138  real_tensor overlap_sum(static_cast<long>(npsi));
139 
140  for (int i=0; i<natom; ++i) {
141  Atom at = mol.get_atom(i);
142  unsigned int atn = at.atomic_number;
143  unsigned int nshell = mol.n_core_orb(atn);
144  if (nshell == 0) continue;
145  for (unsigned int c=0; c<nshell; ++c) {
146  unsigned int l = mol.get_core_l(atn, c);
147  int max_m = (l+1)*(l+2)/2;
148  nshell -= max_m - 1;
149  for (int m=0; m<max_m; ++m) {
150  real_function_3d core = real_factory_3d(world).functor(real_functor_3d(new CoreOrbitalFunctor(mol, i, c, m)));
151  real_tensor overlap = inner(world, core, psi);
152  overlap_sum += overlap;
153  for (int j=0; j<npsi; ++j) {
154  if (include_Bc) overlap[j] *= mol.get_core_bc(atn, c);
155  proj[j] += core.scale(overlap[j]);
156  }
157  }
158  }
159  world.gop.fence();
160  }
161  if (world.rank() == 0) print("sum_k <core_k|psi_i>:", overlap_sum);
162  return proj;
163  }
164 
165  double core_projector_derivative(World & world, const vector_real_function_3d& mo, const real_tensor& occ, int atom, int axis)
166  {
167  vector_real_function_3d cores, dcores;
168  std::vector<double> bc;
169  unsigned int atn = mol.get_atom(atom).atomic_number;
170  unsigned int ncore = mol.n_core_orb(atn);
171 
172  // projecting core & d/dx core
173  for (unsigned int c=0; c<ncore; ++c) {
174  unsigned int l = mol.get_core_l(atn, c);
175  int max_m = (l+1)*(l+2)/2;
176  for (int m=0; m<max_m; ++m) {
178  cores.push_back(real_function_3d(real_factory_3d(world).functor(func).truncate_on_project()));
180  dcores.push_back(real_function_3d(real_factory_3d(world).functor(func).truncate_on_project()));
181  bc.push_back(mol.get_core_bc(atn, c));
182  }
183  }
184 
185  // calc \sum_i occ_i <psi_i|(\sum_c Bc d/dx |core><core|)|psi_i>
186  double r = 0.0;
187  for (unsigned int c=0; c<cores.size(); ++c) {
188  double rcore= 0.0;
189  real_tensor rcores = inner(world, cores[c], mo);
190  real_tensor rdcores = inner(world, dcores[c], mo);
191  for (unsigned int i=0; i<mo.size(); ++i) {
192  rcore += rdcores[i] * rcores[i] * occ[i];
193  }
194  r += 2.0 * bc[c] * rcore;
195  }
196 
197  return r;
198  }
199 
201  if (core_type_.substr(0,3) == "mcp") {
202  // START_TIMER(world);
203  gaxpy(world, 1.0, Vpsi, 1.0, core_projection(world, amo));
204  // END_TIMER(world, "MCP Core Projector");
205  }
206  }
207 
209  double safety = 0.1;
210  double vtol = FunctionDefaults<3>::get_thresh() * safety;
211  vnuc = real_factory_3d(world).functor(real_functor_3d(new MolecularPotentialFunctor(mol))).thresh(vtol).truncate_on_project();
213  vnuc.reconstruct();
214  // "" is legacy core_type value for all-electron (also be used by CorePotentialManager)
215  // "none" is current core_type value for all-electron
216  if (core_type_ != "" && core_type_ != "none") {
217  real_function_3d c_pot = real_factory_3d(world).functor(real_functor_3d(new MolecularCorePotentialFunctor(mol))).thresh(vtol).initial_level(4);
219  c_pot.reconstruct();
220  vnuc += c_pot;
221  vnuc.truncate();
222  }
223  }
224 };
225 }
226 
227 #endif
Declaration of utility class and functions for atom.
Definition: molecule.h:58
unsigned int atomic_number
Atomic number.
Definition: molecule.h:61
Definition: potentialmanager.h:96
double operator()(const coord_3d &r) const
Definition: potentialmanager.h:104
CoreOrbitalDerivativeFunctor(Molecule &molecule, int atom, int axis, unsigned int core, int m)
Definition: potentialmanager.h:102
const Molecule molecule
Definition: potentialmanager.h:97
const int axis
Definition: potentialmanager.h:98
const int atom
Definition: potentialmanager.h:98
const unsigned int core
Definition: potentialmanager.h:99
const int m
Definition: potentialmanager.h:100
Definition: potentialmanager.h:83
CoreOrbitalFunctor(Molecule &molecule, int atom, unsigned int core, int m)
Definition: potentialmanager.h:89
double operator()(const coord_3d &r) const
Definition: potentialmanager.h:91
const int m
Definition: potentialmanager.h:87
const unsigned int core
Definition: potentialmanager.h:86
const int atom
Definition: potentialmanager.h:85
const Molecule molecule
Definition: potentialmanager.h:84
FunctionDefaults holds default paramaters as static class members.
Definition: funcdefaults.h:204
static const double & get_thresh()
Returns the default threshold.
Definition: funcdefaults.h:279
Abstract base class interface required for functors used as input to Functions.
Definition: function_interface.h:68
void set_thresh(double value, bool fence=true)
Sets the value of the truncation threshold. Optional global fence.
Definition: mra.h:577
Function< T, NDIM > & scale(const Q q, bool fence=true)
Inplace, scale the function by a constant. No communication except for optional fence.
Definition: mra.h:953
Function< T, NDIM > & truncate(double tol=0.0, bool fence=true)
Truncate the function with optional fence. Compresses with fence if not compressed.
Definition: mra.h:602
const Function< T, NDIM > & reconstruct(bool fence=true) const
Reconstructs the function, transforming into scaling function basis. Possible non-blocking comm.
Definition: mra.h:775
Definition: potentialmanager.h:69
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: potentialmanager.h:80
MolecularCorePotentialFunctor(const Molecule &molecule)
Definition: potentialmanager.h:73
double operator()(const coord_3d &x) const
Definition: potentialmanager.h:76
const Molecule & molecule
Definition: potentialmanager.h:71
Definition: potentialmanager.h:55
MolecularPotentialFunctor(const Molecule &molecule)
Definition: potentialmanager.h:59
const Molecule & molecule
Definition: potentialmanager.h:57
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: potentialmanager.h:66
double operator()(const coord_3d &x) const
Definition: potentialmanager.h:62
Definition: molecule.h:124
std::vector< madness::Vector< double, 3 > > get_all_coords_vec() const
Definition: molecule.cc:408
double get_core_bc(unsigned int atn, unsigned int c) const
Definition: molecule.h:371
const Atom & get_atom(unsigned int i) const
Definition: molecule.cc:447
double core_derivative(int atom, int axis, unsigned int core, int m, double x, double y, double z) const
Definition: molecule.cc:1103
double molecular_core_potential(double x, double y, double z) const
Definition: molecule.cc:1116
double nuclear_attraction_potential(double x, double y, double z) const
nuclear attraction potential for the whole molecule
Definition: molecule.cc:971
size_t natom() const
Definition: molecule.h:387
double core_eval(int atom, unsigned int core, int m, double x, double y, double z) const
Definition: molecule.cc:1094
unsigned int get_core_l(unsigned int atn, unsigned int c) const
Definition: molecule.h:367
unsigned int n_core_orb(unsigned int atn) const
Definition: molecule.h:360
Definition: potentialmanager.h:110
Molecule mol
Definition: potentialmanager.h:112
const Molecule & molecule() const
Definition: potentialmanager.h:120
void apply_nonlocal_potential(World &world, const vector_real_function_3d &amo, vector_real_function_3d Vpsi)
Definition: potentialmanager.h:200
const real_function_3d & vnuclear()
Definition: potentialmanager.h:128
PotentialManager(const Molecule &molecule, const std::string &core_type)
Definition: potentialmanager.h:117
const std::string & core_type() const
Definition: potentialmanager.h:124
real_function_3d vnuc
Definition: potentialmanager.h:113
vector_real_function_3d core_projection(World &world, const vector_real_function_3d &psi, const bool include_Bc=true)
Definition: potentialmanager.h:132
double core_projector_derivative(World &world, const vector_real_function_3d &mo, const real_tensor &occ, int atom, int axis)
Definition: potentialmanager.h:165
void make_nuclear_potential(World &world)
Definition: potentialmanager.h:208
std::string core_type_
Definition: potentialmanager.h:114
void fence(bool debug=false)
Synchronizes all processes in communicator AND globally ensures no pending AM or tasks.
Definition: worldgop.cc:161
A parallel world class.
Definition: world.h:132
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition: world.h:318
WorldGopInterface & gop
Global operations.
Definition: world.h:205
Declaration of core potential related class.
const double m
Definition: gfit.cc:199
double psi(const Vector< double, 3 > &r)
Definition: hatom_energy.cc:78
Header to declare stuff which has not yet found a home.
Main include file for MADNESS and defines Function interface.
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
std::shared_ptr< FunctionFunctorInterface< double, 3 > > real_functor_3d
Definition: functypedefs.h:107
std::vector< real_function_3d > vector_real_function_3d
Definition: functypedefs.h:79
std::shared_ptr< FunctionFunctorInterface< double, 3 > > func(new opT(g))
FunctionFactory< double, 3 > real_factory_3d
Definition: functypedefs.h:93
void print(const T &t, const Ts &... ts)
Print items to std::cout (items separated by spaces) and terminate with a new line.
Definition: print.h:225
Function< double, 3 > real_function_3d
Definition: functypedefs.h:65
double inner(response_space &a, response_space &b)
Definition: response_functions.h:442
void gaxpy(const double a, ScalarResult< T > &left, const double b, const T &right, const bool fence=true)
the result type of a macrotask must implement gaxpy
Definition: macrotaskq.h:140
static const double c
Definition: relops.cc:10
Defines and implements most of Tensor.
std::size_t axis
Definition: testpdiff.cc:59