Skip to content

Commit f25c3ac

Browse files
mbugeanumbugeanu
committed
added wem_lin_benzene file - although already there before
1 parent fa519c3 commit f25c3ac

File tree

1 file changed

+174
-0
lines changed

1 file changed

+174
-0
lines changed

‎tests/wem_lin/wem_lin_benzene.cpp‎

Lines changed: 174 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,174 @@
1+
/* pcmsolver_copyright_start */
2+
/*
3+
* PCMSolver, an API for the Polarizable Continuum Model
4+
* Copyright (C) 2013 Roberto Di Remigio, Luca Frediani and contributors
5+
*
6+
* This file is part of PCMSolver.
7+
*
8+
* PCMSolver is free software: you can redistribute it and/or modify
9+
* it under the terms of the GNU Lesser General Public License as published by
10+
* the Free Software Foundation, either version 3 of the License, or
11+
* (at your option) any later version.
12+
*
13+
* PCMSolver is distributed in the hope that it will be useful,
14+
* but WITHOUT ANY WARRANTY; without even the implied warranty of
15+
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16+
* GNU Lesser General Public License for more details.
17+
*
18+
* You should have received a copy of the GNU Lesser General Public License
19+
* along with PCMSolver. If not, see <http://www.gnu.org/licenses/>.
20+
*
21+
* For information on the complete list of contributors to the
22+
* PCMSolver API, see: <http://pcmsolver.github.io/pcmsolver-doc>
23+
*/
24+
/* pcmsolver_copyright_end */
25+
26+
#define BOOST_TEST_MODULE WEMSolverbenzene
27+
#define DEBUG
28+
29+
#include <boost/test/unit_test.hpp>
30+
#include <boost/test/floating_point_comparison.hpp>
31+
32+
#include <iostream>
33+
#include <vector>
34+
35+
#include "Config.hpp"
36+
37+
#include "CollocationIntegrator.hpp"
38+
#include "DerivativeTypes.hpp"
39+
#include "WEMSolver.hpp"
40+
#include "UniformDielectric.hpp"
41+
#include "Vacuum.hpp"
42+
#include "WaveletCavity.hpp"
43+
44+
/*! \class WEMSolver
45+
* \test \b benzene tests WEMSolver with linear ansatz functions using ammonia and a wavelet cavity
46+
*/
47+
BOOST_AUTO_TEST_CASE(benzene)
48+
{
49+
printf("STARTING TEST\n");
50+
// Set up cavity
51+
/*
52+
Eigen::Vector3d C1(-0.694303272975, -0.000000000000, -1.202568545351);
53+
Eigen::Vector3d C2( 0.694303272975, 0.000000000000, -1.202568545351);
54+
Eigen::Vector3d C3( 1.388606546154, 0.000000000000, 0.000000000000);
55+
Eigen::Vector3d C4( 0.694303272975, 0.000000000000, 1.202568545351);
56+
Eigen::Vector3d C5(-0.694303272975, -0.000000000000, 1.202568545351);
57+
Eigen::Vector3d C6(-1.388606546154, -0.000000000000, 0.000000000000);
58+
Eigen::Vector3d H1(-1.235418032354, -0.000000000000, -2.139806800843);
59+
Eigen::Vector3d H2( 1.235418032354, 0.000000000000, -2.139806800843);
60+
Eigen::Vector3d H3( 2.470836065313, 0.000000000000, 0.000000000000);
61+
Eigen::Vector3d H4( 1.235418032354, 0.000000000000, 2.139806800843);
62+
Eigen::Vector3d H5(-1.235418032354, -0.000000000000, 2.139806800843);
63+
Eigen::Vector3d H6(-2.470836065313, -0.000000000000, 0.000000000000);
64+
*/
65+
Eigen::Vector3d C1( 5.274, 1.999, -8.568);
66+
Eigen::Vector3d C2( 6.627, 2.018, -8.209);
67+
Eigen::Vector3d C3( 7.366, 0.829, -8.202);
68+
Eigen::Vector3d C4( 6.752, -0.379, -8.554);
69+
Eigen::Vector3d C5( 5.399, -0.398, -8.912);
70+
Eigen::Vector3d C6( 4.660, 0.791, -8.919);
71+
Eigen::Vector3d H1( 4.704, 2.916, -8.573);
72+
Eigen::Vector3d H2( 7.101, 2.950, -7.938);
73+
Eigen::Vector3d H3( 8.410, 0.844, -7.926);
74+
Eigen::Vector3d H4( 7.322, -1.296, -8.548);
75+
Eigen::Vector3d H5( 4.925, -1.330, -9.183);
76+
Eigen::Vector3d H6( 3.616, 0.776, -9.196);
77+
78+
std::vector<Sphere> spheres;
79+
/*
80+
Sphere sph1(C1, 3.212534412);
81+
Sphere sph2(C2, 3.212534412);
82+
Sphere sph3(C3, 3.212534412);
83+
Sphere sph4(C4, 3.212534412);
84+
Sphere sph5(C5, 3.212534412);
85+
Sphere sph6(C6, 3.212534412);
86+
87+
Sphere sph7(H1, 2.267671349);
88+
Sphere sph8(H2, 2.267671349);
89+
Sphere sph9(H3, 2.267671349);
90+
Sphere sph10(H4, 2.267671349);
91+
Sphere sph11(H5, 2.267671349);
92+
Sphere sph12(H6, 2.267671349);
93+
*/
94+
Sphere sph1(C1, 1.7*1.2);
95+
Sphere sph2(C2, 1.7*1.2);
96+
Sphere sph3(C3, 1.7*1.2);
97+
Sphere sph4(C4, 1.7*1.2);
98+
Sphere sph5(C5, 1.7*1.2);
99+
Sphere sph6(C6, 1.7*1.2);
100+
101+
Sphere sph7(H1, 1.2*1.2);
102+
Sphere sph8(H2, 1.2*1.2);
103+
Sphere sph9(H3, 1.2*1.2);
104+
Sphere sph10(H4, 1.2*1.2);
105+
Sphere sph11(H5, 1.2*1.2);
106+
Sphere sph12(H6, 1.2*1.2);
107+
108+
spheres.push_back(sph1);
109+
spheres.push_back(sph2);
110+
spheres.push_back(sph3);
111+
spheres.push_back(sph4);
112+
spheres.push_back(sph5);
113+
spheres.push_back(sph6);
114+
spheres.push_back(sph7);
115+
spheres.push_back(sph8);
116+
spheres.push_back(sph9);
117+
spheres.push_back(sph10);
118+
spheres.push_back(sph11);
119+
spheres.push_back(sph12);
120+
121+
double probeRadius = 1.385*1.2; // Probe Radius for water
122+
int patchLevel = 3;
123+
double coarsity = 0.5;
124+
printf("TEST 1\n");
125+
WaveletCavity cavity(spheres, probeRadius, patchLevel, coarsity);
126+
printf("TEST 1a\n");
127+
cavity.readCavity("molec_dyadic.dat");
128+
cavity.scaleNodePoints(0.52917721092);
129+
130+
printf("TEST 2\n");
131+
CollocationIntegrator * diag = new CollocationIntegrator();
132+
double permittivity = 78.39;
133+
Vacuum<AD_directional> * gfInside = new Vacuum<AD_directional>();
134+
UniformDielectric<AD_directional> * gfOutside = new
135+
UniformDielectric<AD_directional>(permittivity);
136+
int firstKind = 0;
137+
#ifdef DEBUG
138+
FILE* debugFile = fopen("debug.out","w");
139+
fclose(debugFile);
140+
#endif
141+
WEMSolver solver(gfInside, gfOutside, "Linear", firstKind);
142+
solver.buildSystemMatrix(cavity);
143+
cavity.uploadPoints(solver.getQuadratureLevel(), solver.getT_());
144+
145+
double Hcharge = 1.0;
146+
double Ccharge = 6.0;
147+
int size = cavity.size();
148+
Eigen::VectorXd fake_mep = Eigen::VectorXd::Zero(size);
149+
for (int i = 0; i < size; ++i) {
150+
Eigen::Vector3d center = cavity.elementCenter(i);
151+
double C1mep = Ccharge/(center - C1).norm();
152+
double C2mep = Ccharge/(center - C2).norm();
153+
double C3mep = Ccharge/(center - C3).norm();
154+
double C4mep = Ccharge/(center - C4).norm();
155+
double C5mep = Ccharge/(center - C5).norm();
156+
double C6mep = Ccharge/(center - C6).norm();
157+
158+
double H1mep = Hcharge/(center - H1).norm();
159+
double H2mep = Hcharge/(center - H2).norm();
160+
double H3mep = Hcharge/(center - H3).norm();
161+
double H4mep = Hcharge/(center - H4).norm();
162+
double H5mep = Hcharge/(center - H5).norm();
163+
double H6mep = Hcharge/(center - H6).norm();
164+
fake_mep(i) = C1mep + C2mep + C3mep + C4mep + C5mep + C6mep +
165+
H1mep + H2mep + H3mep + H4mep + H5mep + H6mep;
166+
}
167+
// The total ASC for a dielectric is -Q*[(epsilon-1)/epsilon]
168+
Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size);
169+
solver.compCharge(fake_mep, fake_asc);
170+
double totalASC = - (6 * Ccharge + 6 * Hcharge) * ( permittivity - 1) / permittivity;
171+
double totalFakeASC = fake_asc.sum();
172+
std::cout << "totalASC - totalFakeASC = " << totalASC - totalFakeASC << std::endl;
173+
BOOST_REQUIRE_CLOSE(totalASC, totalFakeASC, 3e-2);
174+
}

0 commit comments

Comments
 (0)