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>
51#include <madness/misc/misc.h>
52#include <madness/mra/mra.h>
53
54namespace madness {
56private:
58public:
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
70private:
72public:
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
85 const int atom;
86 const unsigned int core;
87 const int m;
88public:
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;
101public:
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
111private:
114std::string core_type_;
115
116public:
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();
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
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
double thresh() const
Returns value of truncation threshold. No communication.
Definition mra.h:567
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 > & 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
void apply_nonlocal_potential(World &world, const vector_real_function_3d &amo, vector_real_function_3d Vpsi)
Definition potentialmanager.h:200
const Molecule & molecule() const
Definition potentialmanager.h:120
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.
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.
Namespace for all elements and tools of MADNESS.
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
static const double m
Definition relops.cc:9
Defines and implements most of Tensor.
std::size_t axis
Definition testpdiff.cc:59