Skip to content

Commit dc40a37

Browse files
authored
Fix ghost atom handling (#25)
The ghost atom handling did not work if the ghost atoms did not appear strictly after the regular coordinates.
1 parent 1d00e3b commit dc40a37

17 files changed

+816
-199
lines changed

.gitignore

+5
Original file line numberDiff line numberDiff line change
@@ -35,3 +35,8 @@
3535
/build*/
3636
/install*/
3737
/_*/
38+
39+
# project-specific
40+
/scripts/old/*
41+
/scripts/*.h
42+
/scripts/*.cpp

include/dftd_ncoord.h

+1-1
Original file line numberDiff line numberDiff line change
@@ -200,4 +200,4 @@ extern inline double log_cn_cut(const double cn_max, const double cn);
200200

201201
extern inline double dlog_cn_cut(const double cn_max, const double cn);
202202

203-
}; // namespace dftd4
203+
} // namespace dftd4

include/dftd_parameters.h

+3-3
Original file line numberDiff line numberDiff line change
@@ -63,10 +63,10 @@ static const double zeff[119]{
6363
* These values are actually never used in the code. Only r4r2 is used.
6464
*/
6565
static const double r2r4[119]{
66-
0.0, // dummy
67-
8.0589, 3.4698, // H,He
66+
0.0, // dummy
67+
8.0589, 3.4698, // H,He
6868
29.0974, 14.8517, 11.8799, 7.8715, 5.5588, 4.7566, 3.8025,
69-
3.1036, // Li-Ne
69+
3.1036, // Li-Ne
7070
26.1552, 17.2304, 17.7210, 12.7442, 9.5361, 8.1652, 6.7463,
7171
5.6004, // Na-Ar
7272
29.2012, 22.3934, // K,Ca

scripts/README.md

+3
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
# Scripts
2+
3+
Scripts for creating the ORCA source files. This mainly concatenates the corresponding files and adds headers and footers.

scripts/make-qcdftd4.sh

+317
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,317 @@
1+
#!/bin/bash
2+
3+
INCLUDE="../include"
4+
SRC="../src"
5+
6+
cat > qcdftd4param.h << 'EOT'
7+
/*
8+
* This is the DFT-D4 equivalent of copyc6, since this is an *empirical*
9+
* dispersion correction we can hardly avoid having this big, scary blob
10+
* of data lying around somewhere.
11+
*
12+
* Files like this are usually computer generated so avoid modifying it,
13+
* but have the appropriate tool generate it from dftd4.
14+
*
15+
* Responsible for this mess:
16+
* Sebastian Ehlert <[email protected]> (SAW190521)
17+
*
18+
* Update/Fix for parameters from periodic extension by:
19+
* Marvin Friede <[email protected]> (MF161222)
20+
*
21+
* Extension for Fr, Ra and Actinides by:
22+
* Marvin Friede <[email protected]> (MF121223)
23+
*/
24+
25+
#ifndef QCDFTD4PARAM_H
26+
#define QCDFTD4PARAM_H
27+
28+
EOT
29+
30+
sed -n '/namespace/,$p' "${INCLUDE}/dftd_parameters.h" >> qcdftd4param.h
31+
32+
# Append the Footer:
33+
printf "\n#endif // QCDFTD4PARAM_H\n" >> qcdftd4param.h
34+
35+
36+
##################################################################
37+
##################################################################
38+
39+
dosed(){
40+
sed -n '/namespace dftd4 {/,/} \/\/ namespace dftd4/p' $1 | \
41+
sed '/namespace dftd4 {/d; /} \/\/ namespace dftd4/d' | \
42+
awk '/./,EOF' > $2
43+
}
44+
45+
dosed "${INCLUDE}/dftd_cutoff.h" cutoff.txt
46+
dosed "${INCLUDE}/dftd_model.h" model.txt
47+
dosed "${INCLUDE}/dftd_dispersion.h" dispersion.txt
48+
dosed "${INCLUDE}/damping/dftd_rational.h" rational.txt
49+
dosed "${INCLUDE}/damping/dftd_atm.h" atm.txt
50+
51+
cat > qcdftd4.h << 'EOT'
52+
/*
53+
* D4(EEQ)-ATM implementation
54+
*
55+
* Sebastian Ehlert <[email protected]> (SAW190521)
56+
* Marvin Friede <[email protected]> (MF161222)
57+
*/
58+
59+
#ifndef __QCDFTD4_H
60+
#define __QCDFTD4_H
61+
62+
#include "qcinpdat.h" // TGeomInput class
63+
#include "qcmat1.h" // TRVector and TRMatrix class
64+
65+
namespace dftd4 {
66+
67+
EOT
68+
69+
######
70+
71+
cat >> qcdftd4.h << 'EOT'
72+
/* --------------------------------------------------------------------------
73+
74+
Cutoff (dftd_cutoff.h)
75+
https://github.com/dftd4/cpp-d4/blob/main/include/dftd_cutoff.h
76+
77+
-------------------------------------------------------------------------- */
78+
79+
EOT
80+
81+
cat cutoff.txt >> qcdftd4.h
82+
rm cutoff.txt
83+
84+
######
85+
86+
cat >> qcdftd4.h << 'EOT'
87+
/* --------------------------------------------------------------------------
88+
89+
D4 model (dftd_model.h)
90+
https://github.com/dftd4/cpp-d4/blob/main/include/dftd_model.h
91+
92+
-------------------------------------------------------------------------- */
93+
94+
EOT
95+
96+
cat model.txt >> qcdftd4.h
97+
rm model.txt
98+
99+
######
100+
101+
cat >> qcdftd4.h << 'EOT'
102+
/* --------------------------------------------------------------------------
103+
104+
Dispersion (dftd_dispersion.h)
105+
https://github.com/dftd4/cpp-d4/blob/main/include/dftd_dispersion.h
106+
107+
-------------------------------------------------------------------------- */
108+
109+
EOT
110+
111+
cat dispersion.txt >> qcdftd4.h
112+
rm dispersion.txt
113+
114+
######
115+
116+
cat >> qcdftd4.h << 'EOT'
117+
/* --------------------------------------------------------------------------
118+
119+
Dispersion (damping/dftd_rational.h)
120+
https://github.com/dftd4/cpp-d4/blob/main/include/damping/dftd_rational.h
121+
122+
-------------------------------------------------------------------------- */
123+
124+
EOT
125+
126+
cat rational.txt >> qcdftd4.h
127+
rm rational.txt
128+
129+
######
130+
131+
cat >> qcdftd4.h << 'EOT'
132+
/* --------------------------------------------------------------------------
133+
134+
Dispersion (damping/dftd_atm.h)
135+
https://github.com/dftd4/cpp-d4/blob/main/include/damping/dftd_atm.h
136+
137+
-------------------------------------------------------------------------- */
138+
139+
EOT
140+
141+
cat atm.txt >> qcdftd4.h
142+
rm atm.txt
143+
144+
######
145+
146+
cat >> qcdftd4.h << 'EOT'
147+
} // namespace dftd4
148+
149+
/* --------------------------------------------------------------------------------------
150+
// Calculates the EEQ charges according to the D4 paper
151+
// https://doi.org/10.1063/1.5090222
152+
// Bernardo de Souza, 14/09/2023
153+
-------------------------------------------------------------------------------------- */
154+
int CalcEEQCharges(TRMatrix &XYZ, TIVector &ATNO, int NAtoms, int totalcharge, TRVector &q, bool printerror=true);
155+
156+
#endif // __QCDFTD4_H
157+
158+
EOT
159+
160+
161+
##################################################################
162+
##################################################################
163+
164+
165+
dosed "${SRC}/dftd_cutoff.cpp" cutoff.txt
166+
dosed "${SRC}/dftd_model.cpp" model.txt
167+
dosed "${SRC}/dftd_dispersion.cpp" dispersion.txt
168+
dosed "${SRC}/damping/dftd_rational.cpp" rational.txt
169+
dosed "${SRC}/damping/dftd_atm.cpp" atm.txt
170+
171+
cat > qcdftd4.cpp << 'EOT'
172+
/*
173+
* D4(EEQ)-ATM implementation
174+
*
175+
* Sebastian Ehlert <[email protected]> (SAW190521)
176+
* Marvin Friede <[email protected]> (MF161222)
177+
*/
178+
#include <cmath>
179+
#include <limits>
180+
181+
#include "qcinpdat.h" // TGeomInput class
182+
#include "qcmat2.h" // BLAS routines
183+
#include "qcmath.h" // TVector and TMatrix class
184+
185+
// we cannot avoid one big, scary parameter file...
186+
#include "qcdftd4param.h"
187+
188+
#include "qceeq.h"
189+
#include "qcncoord.h"
190+
191+
// always include self
192+
#include "qcdftd4.h"
193+
194+
/* --------------------------------------------------------------------------------------
195+
// Calculates the EEQ charges according to the D4 paper
196+
// https://doi.org/10.1063/1.5090222
197+
// Bernardo de Souza, 14/09/2023
198+
-------------------------------------------------------------------------------------- */
199+
int CalcEEQCharges(TRMatrix &XYZ, TIVector &ATNO, int NAtoms, int totalcharge, TRVector &q, bool printerror){
200+
201+
// initialize the charges to zero
202+
q.Init();
203+
204+
// check atomic numbers to guarantee we have all parameters
205+
// ghost atoms (ATNO=0) will have no charge
206+
for (int i=0;i<NAtoms;i++)
207+
if (ATNO(i)>86){
208+
if (printerror) printMessage("Atomic number %d detected. EEQ charges can not be calculated.\n",ATNO(i));
209+
return 1;
210+
}
211+
212+
TGeomInput molecule;
213+
molecule.NAtoms=NAtoms;
214+
molecule.CC.CopyMat(XYZ);
215+
molecule.ATNO.CopyVec(ATNO);
216+
TIVector Index(NAtoms);
217+
for (int i=0;i<NAtoms;i++) Index(i)=i;
218+
TRMatrix dist;
219+
dist.NewMatrix(NAtoms, NAtoms);
220+
dftd4::calc_distances(molecule, Index, dist);
221+
222+
TRMatrix dqdr;
223+
224+
return dftd4::get_charges(molecule, Index, dist, totalcharge, 25, q, dqdr, false);
225+
}
226+
227+
namespace dftd4 {
228+
229+
EOT
230+
231+
######
232+
233+
cat >> qcdftd4.cpp << 'EOT'
234+
/* --------------------------------------------------------------------------
235+
236+
Cutoff (dftd_cutoff.cpp)
237+
https://github.com/dftd4/cpp-d4/blob/main/src/dftd_cutoff.cpp
238+
239+
-------------------------------------------------------------------------- */
240+
241+
EOT
242+
243+
cat cutoff.txt >> qcdftd4.cpp
244+
rm cutoff.txt
245+
246+
######
247+
248+
cat >> qcdftd4.cpp << 'EOT'
249+
/* --------------------------------------------------------------------------
250+
251+
D4 model (dftd_model.cpp)
252+
https://github.com/dftd4/cpp-d4/blob/main/src/dftd_model.cpp
253+
254+
-------------------------------------------------------------------------- */
255+
256+
EOT
257+
258+
cat model.txt >> qcdftd4.cpp
259+
rm model.txt
260+
261+
######
262+
263+
cat >> qcdftd4.cpp << 'EOT'
264+
/* --------------------------------------------------------------------------
265+
266+
Dispersion (dftd_dispersion.cpp)
267+
https://github.com/dftd4/cpp-d4/blob/main/src/dftd_dispersion.cpp
268+
269+
-------------------------------------------------------------------------- */
270+
271+
EOT
272+
273+
cat dispersion.txt >> qcdftd4.cpp
274+
rm dispersion.txt
275+
276+
######
277+
278+
cat >> qcdftd4.cpp << 'EOT'
279+
/* --------------------------------------------------------------------------
280+
281+
Dispersion (damping/dftd_rational.cpp)
282+
https://github.com/dftd4/cpp-d4/blob/main/src/damping/dftd_rational.cpp
283+
284+
-------------------------------------------------------------------------- */
285+
286+
EOT
287+
288+
cat rational.txt >> qcdftd4.cpp
289+
rm rational.txt
290+
291+
######
292+
293+
cat >> qcdftd4.cpp << 'EOT'
294+
/* --------------------------------------------------------------------------
295+
296+
Dispersion (damping/dftd_atm.cpp)
297+
https://github.com/dftd4/cpp-d4/blob/main/src/damping/dftd_atm.cpp
298+
299+
-------------------------------------------------------------------------- */
300+
301+
EOT
302+
303+
cat atm.txt >> qcdftd4.cpp
304+
rm atm.txt
305+
306+
######
307+
308+
cat >> qcdftd4.cpp << 'EOT'
309+
} // namespace dftd4
310+
EOT
311+
312+
313+
##########################################################################
314+
315+
316+
sed -i "s/TMolecule/TGeomInput/" *.h
317+
sed -i "s/TMolecule/TGeomInput/" *.cpp

scripts/make-qcdftd4param.sh

+34
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
#!/bin/bash
2+
3+
INCLUDE="../include"
4+
5+
cat > qcdftd4param.h << 'EOT'
6+
/*
7+
* This is the DFT-D4 equivalent of copyc6, since this is an *empirical*
8+
* dispersion correction we can hardly avoid having this big, scary blob
9+
* of data lying around somewhere.
10+
*
11+
* Files like this are usually computer generated so avoid modifying it,
12+
* but have the appropriate tool generate it from dftd4.
13+
*
14+
* Responsible for this mess:
15+
* Sebastian Ehlert <[email protected]> (SAW190521)
16+
*
17+
* Update/Fix for parameters from periodic extension by:
18+
* Marvin Friede <[email protected]> (MF161222)
19+
*
20+
* Extension for Fr, Ra and Actinides by:
21+
* Marvin Friede <[email protected]> (MF121223)
22+
*/
23+
24+
#ifndef QCDFTD4PARAM_H
25+
#define QCDFTD4PARAM_H
26+
27+
EOT
28+
29+
sed -n '/namespace/,$p' "${INCLUDE}/dftd_parameters.h" >> qcdftd4param.h
30+
31+
printf "\n#endif // QCDFTD4PARAM_H\n" >> qcdftd4param.h
32+
33+
34+
sed -i "s/TMolecule/TGeomInput/" qcdftd4param.h

0 commit comments

Comments
 (0)