25
25
26
26
/* * The MRChem sandbox */
27
27
28
- #include " MRCPP/Parallel"
29
- #include < MRCPP/MWOperators>
30
28
#include < MRCPP/Printer>
31
29
#include < MRCPP/Timer>
32
-
33
- #include " chemistry/Nucleus.h"
34
- #include " chemistry/PeriodicTable.h"
35
- #include " environment/Cavity.h"
36
- #include " environment/DHScreening.h"
37
- #include " environment/GPESolver.h"
38
- #include " environment/Permittivity.h"
39
-
40
- #include " qmfunctions/qmfunction_fwd.h"
41
- #include " qmoperators/two_electron/ReactionOperator.h"
30
+ #include " MRCPP/Parallel"
42
31
43
32
#include " mrchem.h"
44
33
#include " mrenv.h"
45
34
46
35
// Initializing global variables
47
36
mrcpp::MultiResolutionAnalysis<3 > *mrchem::MRA;
48
37
49
- using mrcpp::ABGVOperator;
50
-
51
38
using json = nlohmann::json;
52
39
using Timer = mrcpp::Timer;
53
40
using namespace mrchem ;
54
41
55
- using DerivativeOperator = mrcpp::DerivativeOperator<3 >;
56
- using DerivativeOperator_p = std::shared_ptr<mrcpp::DerivativeOperator<3 >>;
57
-
58
- using PoissonOperator = mrcpp::PoissonOperator;
59
- using PoissonOperator_p = std::shared_ptr<mrcpp::PoissonOperator>;
60
-
61
42
int main (int argc, char **argv) {
62
43
mrcpp::mpi::initialize ();
63
44
const auto json_inp = mrenv::fetch_json (argc, argv);
@@ -66,72 +47,6 @@ int main(int argc, char **argv) {
66
47
Timer timer;
67
48
68
49
// Do your stuff here
69
- bool acc_pot = true ;
70
- bool dyn_thrs = false ;
71
- int kain = 7 ;
72
- int max_iter = 200 ;
73
- double poisson_prec = 1.0e-5 ;
74
- double eps_in = 1.0 ;
75
- double eps_out = 78.54 ;
76
- double kappa_out = 0.054995 ;
77
- double slope = 0.2 ;
78
-
79
- std::cout << " setting up the cavity" << std::endl;
80
-
81
- std::vector<double > R = {3.7794522509156563 }; // must be 2 A
82
- auto sph_coords = std::vector<mrcpp::Coord<3 >>({{0.0 , 0.0 , 0.0 }});
83
- Cavity C (sph_coords, R, slope);
84
-
85
- std::cout << " setting up permittivity and kappa" << std::endl;
86
- Permittivity eps (C, eps_in, eps_out, " exponential" );
87
- DHScreening kappa_sq (C, kappa_out, " continuous" );
88
-
89
- std::cout << " setting up the operators" << std::endl;
90
- PoissonOperator_p P_p = std::make_shared<PoissonOperator>(*MRA, poisson_prec);
91
- DerivativeOperator_p D_p = std::make_shared<ABGVOperator<3 >>(*MRA, 0.0 , 0.0 );
92
-
93
- std::cout << " setting up the nuclei" << std::endl;
94
- PeriodicTable P;
95
- auto q_coords = std::vector<mrcpp::Coord<3 >>({{0.7558904498503081 , 0.0 , 0.0 },
96
- {0.0 , 1.5117808997006161 , 0.0 },
97
- {0.0 , 0.0 , 2.267671349550924 },
98
- {0.0 , 0.0 , -0.7558904498503081 },
99
- {-1.5117808997006161 , 0.0 , 0.0 },
100
- {0.0 , -2.267671349550924 , 0.0 }});
101
-
102
- std::cout << " please work" ;
103
- Nucleus Q1 (P.getElement (0 ), q_coords[0 ]);
104
- Nucleus Q2 (P.getElement (0 ), q_coords[1 ]);
105
- Nucleus Q3 (P.getElement (0 ), q_coords[2 ]);
106
- Nucleus Q4 (P.getElement (0 ), q_coords[3 ]);
107
- Nucleus Q5 (P.getElement (0 ), q_coords[4 ]);
108
- Nucleus Q6 (P.getElement (0 ), q_coords[5 ]);
109
-
110
- Nuclei N;
111
- N.push_back (Q1);
112
- N.push_back (Q2);
113
- N.push_back (Q3);
114
- N.push_back (Q4);
115
- N.push_back (Q5);
116
- N.push_back (Q6);
117
-
118
- auto ORB = std::make_shared<mrchem::OrbitalVector>();
119
- ORB->push_back (Orbital (SPIN::Paired));
120
-
121
- std::cout << " setting up the solver" << std::endl;
122
- auto PB_solver = std::make_unique<GPESolver>(eps, N, P_p, D_p, poisson_prec, kain, max_iter, acc_pot, dyn_thrs, " nuclear" );
123
- // auto PB_solver = std::make_unique<GPESolver>(eps, kappa_sq, N, P_p, D_p, poisson_prec, kain, max_iter, acc_pot, dyn_thrs, "nuclear");
124
- std::cout << " Computing the potential standard: " << std::endl;
125
-
126
- PB_solver->setConvergenceThreshold (poisson_prec);
127
- // PB_solver->setStaticSalt();
128
- PB_solver->solveEquation (poisson_prec, ORB);
129
- double total_energy = PB_solver->getTotalEnergy () / 2 ;
130
- PB_solver->clear ();
131
-
132
- total_energy = PB_solver->getTotalEnergy () / 2.0 ;
133
- std::cout << " total electrostatic potential linearized pb: " << total_energy << std::endl;
134
-
135
50
println (0 , json_inp.dump (2 ));
136
51
137
52
timer.stop ();
0 commit comments