Skip to content

Commit 139c395

Browse files
author
Christian Schafmeisterr
committed
Add new kernels
1 parent 305202f commit 139c395

13 files changed

Lines changed: 19758 additions & 0 deletions

‎include/cando/chem/atomGrid.h‎

Lines changed: 97 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,97 @@
1+
/*
2+
File: atomGrid.h
3+
*/
4+
/*
5+
Open Source License
6+
Copyright (c) 2016, Christian E. Schafmeister
7+
Permission is hereby granted, free of charge, to any person obtaining a copy
8+
of this software and associated documentation files (the "Software"), to deal
9+
in the Software without restriction, including without limitation the rights
10+
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
11+
copies of the Software, and to permit persons to whom the Software is
12+
furnished to do so, subject to the following conditions:
13+
The above copyright notice and this permission notice shall be included in
14+
all copies or substantial portions of the Software.
15+
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16+
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17+
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18+
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19+
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20+
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
21+
THE SOFTWARE.
22+
23+
This is an open source license for the CANDO software from Temple University, but it is not the only one. Contact Temple University at mailto:techtransfer@temple.edu if you would like a different license.
24+
*/
25+
/* -^- */
26+
27+
28+
29+
#ifndef AtomGrid_H //[
30+
#define AtomGrid_H
31+
32+
33+
34+
#include <string>
35+
#include <vector>
36+
#include <set>
37+
#include <clasp/core/common.h>
38+
#include <clasp/core/smallMap.h>
39+
#include <clasp/core/array.h>
40+
#include <clasp/core/array_double.h>
41+
#include <cando/chem/nVector.h>
42+
43+
#include <cando/chem/chemPackage.h>
44+
45+
namespace chem {
46+
47+
48+
SMART(AtomGrid);
49+
class AtomGrid_O : public core::CxxObject_O
50+
{
51+
LISP_CLASS(chem,ChemPkg,AtomGrid_O,"AtomGrid",core::CxxObject_O);
52+
53+
public:
54+
private:
55+
double _OriginX;
56+
double _OriginY;
57+
double _OriginZ;
58+
double _CellSize;
59+
size_t _Nx;
60+
size_t _Ny;
61+
size_t _Nz;
62+
gc::Vec0<core::T_sp> _Cells;
63+
gc::Vec0<double> _CellMaxRadius;
64+
NVector_sp _Positions;
65+
core::SimpleVector_double_sp _EffRadii;
66+
core::ComplexVector_double_sp _HeapBounds;
67+
core::ComplexVector_size_t_sp _HeapCells;
68+
uint8_t _VisitedFlag;
69+
gc::Vec0<uint8_t> _Visited;
70+
public:
71+
bool fieldsp() const { return true; };
72+
void fields(core::Record_sp node);
73+
74+
static AtomGrid_sp makeAtomGrid(NVector_sp positions,
75+
core::SimpleVector_double_sp effRadii,
76+
double cellSize,
77+
double padding );
78+
79+
size_t cellIndex(size_t ix, size_t iy, size_t iz) const;
80+
void decodeIndex(size_t idx, size_t& cx, size_t& cy, size_t& cz) const;
81+
double cellMinDist(double x, double y, double z, size_t ix, size_t iy, size_t iz) const;
82+
83+
double sdfGridDouble_(double x, double y, double z, size_t& bestAtom );
84+
core::T_mv sdfGridDouble(double x, double y, double z );
85+
core::T_mv sdfGridSingle(float x, float y, float z );
86+
87+
void buildAtomGrid(NVector_sp positions,
88+
core::SimpleVector_double_sp effRadii,
89+
double cellSize,
90+
double padding = 0.0);
91+
};
92+
93+
94+
95+
96+
};
97+
#endif //]

‎include/cando/chem/energyKernels/anchor.c‎

Lines changed: 278 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,278 @@
1+
template <typename HESSIAN>
2+
struct Anchor {
3+
static constexpr size_t PositionSize = 3;
4+
static std::string description() { return "mathkernel-anchor"; };
5+
void energy(double ka, double xa, double ya, double za, size_t i3x1, double* position, double* energy_accumulate, double* force, HESSIAN hessian, double* dvec, double* hdvec) {
6+
/* !BASE */
7+
double dx;
8+
double dy;
9+
double dz;
10+
double energy;
11+
double r2;
12+
DOUBLE x1 = position[i3x1 + 0];
13+
DOUBLE y1 = position[i3x1 + 1];
14+
DOUBLE z1 = position[i3x1 + 2];
15+
{
16+
/* !BASE */
17+
dx = (x1 + (-(xa)));
18+
dy = (y1 + (-(ya)));
19+
dz = (z1 + (-(za)));
20+
r2 = ((dx * dx) + (dy * dy) + (dz * dz));
21+
energy = (ka * r2);
22+
*energy_accumulate += energy;
23+
}
24+
}
25+
void energy_fd(double ka, double xa, double ya, double za, size_t i3x1, double* position, double* energy_accumulate, double* force, HESSIAN hessian, double* dvec, double* hdvec)
26+
{
27+
energy(ka, xa, ya, za, i3x1, position, energy_accumulate, force, hessian, dvec, hdvec);
28+
}
29+
30+
void gradient(double ka, double xa, double ya, double za, size_t i3x1, double* position, double* energy_accumulate, double* force, HESSIAN hessian, double* dvec, double* hdvec) {
31+
/* !BASE */
32+
double cse_p1_t1_g16996;
33+
double dx;
34+
double dy;
35+
double dz;
36+
double energy;
37+
double g_x1;
38+
double g_y1;
39+
double g_z1;
40+
double r2;
41+
DOUBLE x1 = position[i3x1 + 0];
42+
DOUBLE y1 = position[i3x1 + 1];
43+
DOUBLE z1 = position[i3x1 + 2];
44+
{
45+
/* !BASE */
46+
dx = (x1 + (-(xa)));
47+
dy = (y1 + (-(ya)));
48+
dz = (z1 + (-(za)));
49+
r2 = ((dx * dx) + (dy * dy) + (dz * dz));
50+
energy = (ka * r2);
51+
*energy_accumulate += energy;
52+
cse_p1_t1_g16996 = (2.0000000000000000 * ka);
53+
g_x1 = (cse_p1_t1_g16996 * dx);
54+
KernelGradientAcc(i3x1, 0, g_x1);
55+
g_y1 = (cse_p1_t1_g16996 * dy);
56+
KernelGradientAcc(i3x1, 1, g_y1);
57+
g_z1 = (cse_p1_t1_g16996 * dz);
58+
KernelGradientAcc(i3x1, 2, g_z1);
59+
}
60+
}
61+
void gradient_fd(double ka, double xa, double ya, double za, size_t i3x1, double* position, double* energy_accumulate, double* force, HESSIAN hessian, double* dvec, double* hdvec)
62+
{
63+
constexpr size_t PositionSize = 3;
64+
const double h = 1.0e-5;
65+
const double inv2h = 1.0/(2.0*h);
66+
double e0 = 0.0;
67+
energy(ka, xa, ya, za, i3x1, position, &e0, 0, 0, 0, 0);
68+
if (energy_accumulate) { *energy_accumulate += e0; }
69+
{
70+
double saved = position[i3x1 + 0];
71+
double e_plus = 0.0;
72+
double e_minus = 0.0;
73+
position[i3x1 + 0] = saved + h;
74+
energy(ka, xa, ya, za, i3x1, position, &e_plus, 0, 0, 0, 0);
75+
position[i3x1 + 0] = saved - h;
76+
energy(ka, xa, ya, za, i3x1, position, &e_minus, 0, 0, 0, 0);
77+
position[i3x1 + 0] = saved;
78+
double d = (e_plus - e_minus) * inv2h;
79+
KernelGradientAcc(i3x1, 0, d);
80+
}
81+
{
82+
double saved = position[i3x1 + 1];
83+
double e_plus = 0.0;
84+
double e_minus = 0.0;
85+
position[i3x1 + 1] = saved + h;
86+
energy(ka, xa, ya, za, i3x1, position, &e_plus, 0, 0, 0, 0);
87+
position[i3x1 + 1] = saved - h;
88+
energy(ka, xa, ya, za, i3x1, position, &e_minus, 0, 0, 0, 0);
89+
position[i3x1 + 1] = saved;
90+
double d = (e_plus - e_minus) * inv2h;
91+
KernelGradientAcc(i3x1, 1, d);
92+
}
93+
{
94+
double saved = position[i3x1 + 2];
95+
double e_plus = 0.0;
96+
double e_minus = 0.0;
97+
position[i3x1 + 2] = saved + h;
98+
energy(ka, xa, ya, za, i3x1, position, &e_plus, 0, 0, 0, 0);
99+
position[i3x1 + 2] = saved - h;
100+
energy(ka, xa, ya, za, i3x1, position, &e_minus, 0, 0, 0, 0);
101+
position[i3x1 + 2] = saved;
102+
double d = (e_plus - e_minus) * inv2h;
103+
KernelGradientAcc(i3x1, 2, d);
104+
}
105+
}
106+
107+
void hessian(double ka, double xa, double ya, double za, size_t i3x1, double* position, double* energy_accumulate, double* force, HESSIAN hessian, double* dvec, double* hdvec) {
108+
/* !BASE */
109+
double cse_p1_t1_g16997;
110+
double dx;
111+
double dy;
112+
double dz;
113+
double energy;
114+
double g_x1;
115+
double g_y1;
116+
double g_z1;
117+
double r2;
118+
DOUBLE x1 = position[i3x1 + 0];
119+
DOUBLE y1 = position[i3x1 + 1];
120+
DOUBLE z1 = position[i3x1 + 2];
121+
{
122+
/* !BASE */
123+
dx = (x1 + (-(xa)));
124+
dy = (y1 + (-(ya)));
125+
dz = (z1 + (-(za)));
126+
r2 = ((dx * dx) + (dy * dy) + (dz * dz));
127+
energy = (ka * r2);
128+
*energy_accumulate += energy;
129+
cse_p1_t1_g16997 = (2.0000000000000000 * ka);
130+
g_x1 = (cse_p1_t1_g16997 * dx);
131+
KernelGradientAcc(i3x1, 0, g_x1);
132+
g_y1 = (cse_p1_t1_g16997 * dy);
133+
KernelGradientAcc(i3x1, 1, g_y1);
134+
g_z1 = (cse_p1_t1_g16997 * dz);
135+
KernelGradientAcc(i3x1, 2, g_z1);
136+
}
137+
}
138+
void hessian_fd(double ka, double xa, double ya, double za, size_t i3x1, double* position, double* energy_accumulate, double* force, HESSIAN hessian, double* dvec, double* hdvec)
139+
{
140+
constexpr size_t PositionSize = 3;
141+
const double h = 1.0e-5;
142+
const double inv2h = 1.0/(2.0*h);
143+
const double invh2 = 1.0/((h*h));
144+
double e0 = 0.0;
145+
energy(ka, xa, ya, za, i3x1, position, &e0, 0, 0, 0, 0);
146+
if (energy_accumulate) { *energy_accumulate += e0; }
147+
{
148+
double saved = position[i3x1 + 0];
149+
double e_plus = 0.0;
150+
double e_minus = 0.0;
151+
position[i3x1 + 0] = saved + h;
152+
energy(ka, xa, ya, za, i3x1, position, &e_plus, 0, 0, 0, 0);
153+
position[i3x1 + 0] = saved - h;
154+
energy(ka, xa, ya, za, i3x1, position, &e_minus, 0, 0, 0, 0);
155+
position[i3x1 + 0] = saved;
156+
double d = (e_plus - e_minus) * inv2h;
157+
KernelGradientAcc(i3x1, 0, d);
158+
}
159+
{
160+
double saved = position[i3x1 + 1];
161+
double e_plus = 0.0;
162+
double e_minus = 0.0;
163+
position[i3x1 + 1] = saved + h;
164+
energy(ka, xa, ya, za, i3x1, position, &e_plus, 0, 0, 0, 0);
165+
position[i3x1 + 1] = saved - h;
166+
energy(ka, xa, ya, za, i3x1, position, &e_minus, 0, 0, 0, 0);
167+
position[i3x1 + 1] = saved;
168+
double d = (e_plus - e_minus) * inv2h;
169+
KernelGradientAcc(i3x1, 1, d);
170+
}
171+
{
172+
double saved = position[i3x1 + 2];
173+
double e_plus = 0.0;
174+
double e_minus = 0.0;
175+
position[i3x1 + 2] = saved + h;
176+
energy(ka, xa, ya, za, i3x1, position, &e_plus, 0, 0, 0, 0);
177+
position[i3x1 + 2] = saved - h;
178+
energy(ka, xa, ya, za, i3x1, position, &e_minus, 0, 0, 0, 0);
179+
position[i3x1 + 2] = saved;
180+
double d = (e_plus - e_minus) * inv2h;
181+
KernelGradientAcc(i3x1, 2, d);
182+
}
183+
{
184+
double saved = position[i3x1 + 0];
185+
double e_plus = 0.0;
186+
double e_minus = 0.0;
187+
position[i3x1 + 0] = saved + h;
188+
energy(ka, xa, ya, za, i3x1, position, &e_plus, 0, 0, 0, 0);
189+
position[i3x1 + 0] = saved - h;
190+
energy(ka, xa, ya, za, i3x1, position, &e_minus, 0, 0, 0, 0);
191+
position[i3x1 + 0] = saved;
192+
double hval = (e_plus + e_minus - (2.0*e0)) * invh2;
193+
KernelHessDiagAcc( PositionSize, hessian, dvec, hdvec, i3x1, 0, i3x1, 0, hval);
194+
}
195+
{
196+
double saved = position[i3x1 + 1];
197+
double e_plus = 0.0;
198+
double e_minus = 0.0;
199+
position[i3x1 + 1] = saved + h;
200+
energy(ka, xa, ya, za, i3x1, position, &e_plus, 0, 0, 0, 0);
201+
position[i3x1 + 1] = saved - h;
202+
energy(ka, xa, ya, za, i3x1, position, &e_minus, 0, 0, 0, 0);
203+
position[i3x1 + 1] = saved;
204+
double hval = (e_plus + e_minus - (2.0*e0)) * invh2;
205+
KernelHessDiagAcc( PositionSize, hessian, dvec, hdvec, i3x1, 1, i3x1, 1, hval);
206+
}
207+
{
208+
double saved = position[i3x1 + 2];
209+
double e_plus = 0.0;
210+
double e_minus = 0.0;
211+
position[i3x1 + 2] = saved + h;
212+
energy(ka, xa, ya, za, i3x1, position, &e_plus, 0, 0, 0, 0);
213+
position[i3x1 + 2] = saved - h;
214+
energy(ka, xa, ya, za, i3x1, position, &e_minus, 0, 0, 0, 0);
215+
position[i3x1 + 2] = saved;
216+
double hval = (e_plus + e_minus - (2.0*e0)) * invh2;
217+
KernelHessDiagAcc( PositionSize, hessian, dvec, hdvec, i3x1, 2, i3x1, 2, hval);
218+
}
219+
{
220+
double saved_i = position[i3x1 + 1];
221+
double saved_j = position[i3x1 + 0];
222+
double e_pp = 0.0;
223+
double e_pm = 0.0;
224+
double e_mp = 0.0;
225+
double e_mm = 0.0;
226+
position[i3x1 + 1] = saved_i + h; position[i3x1 + 0] = saved_j + h;
227+
energy(ka, xa, ya, za, i3x1, position, &e_pp, 0, 0, 0, 0);
228+
position[i3x1 + 0] = saved_j - h;
229+
energy(ka, xa, ya, za, i3x1, position, &e_pm, 0, 0, 0, 0);
230+
position[i3x1 + 1] = saved_i - h; position[i3x1 + 0] = saved_j + h;
231+
energy(ka, xa, ya, za, i3x1, position, &e_mp, 0, 0, 0, 0);
232+
position[i3x1 + 0] = saved_j - h;
233+
energy(ka, xa, ya, za, i3x1, position, &e_mm, 0, 0, 0, 0);
234+
position[i3x1 + 1] = saved_i; position[i3x1 + 0] = saved_j;
235+
double hval = (e_pp - e_pm - e_mp + e_mm) * (0.25*invh2);
236+
KernelHessOffDiagAcc( PositionSize, hessian, dvec, hdvec, i3x1, 1, i3x1, 0, hval);
237+
}
238+
{
239+
double saved_i = position[i3x1 + 2];
240+
double saved_j = position[i3x1 + 0];
241+
double e_pp = 0.0;
242+
double e_pm = 0.0;
243+
double e_mp = 0.0;
244+
double e_mm = 0.0;
245+
position[i3x1 + 2] = saved_i + h; position[i3x1 + 0] = saved_j + h;
246+
energy(ka, xa, ya, za, i3x1, position, &e_pp, 0, 0, 0, 0);
247+
position[i3x1 + 0] = saved_j - h;
248+
energy(ka, xa, ya, za, i3x1, position, &e_pm, 0, 0, 0, 0);
249+
position[i3x1 + 2] = saved_i - h; position[i3x1 + 0] = saved_j + h;
250+
energy(ka, xa, ya, za, i3x1, position, &e_mp, 0, 0, 0, 0);
251+
position[i3x1 + 0] = saved_j - h;
252+
energy(ka, xa, ya, za, i3x1, position, &e_mm, 0, 0, 0, 0);
253+
position[i3x1 + 2] = saved_i; position[i3x1 + 0] = saved_j;
254+
double hval = (e_pp - e_pm - e_mp + e_mm) * (0.25*invh2);
255+
KernelHessOffDiagAcc( PositionSize, hessian, dvec, hdvec, i3x1, 2, i3x1, 0, hval);
256+
}
257+
{
258+
double saved_i = position[i3x1 + 2];
259+
double saved_j = position[i3x1 + 1];
260+
double e_pp = 0.0;
261+
double e_pm = 0.0;
262+
double e_mp = 0.0;
263+
double e_mm = 0.0;
264+
position[i3x1 + 2] = saved_i + h; position[i3x1 + 1] = saved_j + h;
265+
energy(ka, xa, ya, za, i3x1, position, &e_pp, 0, 0, 0, 0);
266+
position[i3x1 + 1] = saved_j - h;
267+
energy(ka, xa, ya, za, i3x1, position, &e_pm, 0, 0, 0, 0);
268+
position[i3x1 + 2] = saved_i - h; position[i3x1 + 1] = saved_j + h;
269+
energy(ka, xa, ya, za, i3x1, position, &e_mp, 0, 0, 0, 0);
270+
position[i3x1 + 1] = saved_j - h;
271+
energy(ka, xa, ya, za, i3x1, position, &e_mm, 0, 0, 0, 0);
272+
position[i3x1 + 2] = saved_i; position[i3x1 + 1] = saved_j;
273+
double hval = (e_pp - e_pm - e_mp + e_mm) * (0.25*invh2);
274+
KernelHessOffDiagAcc( PositionSize, hessian, dvec, hdvec, i3x1, 2, i3x1, 1, hval);
275+
}
276+
}
277+
278+
};

0 commit comments

Comments
 (0)