IGMPlot 3.16
Optimized IGMplot version able to use wfn/wfx/xyz files
Loading...
Searching...
No Matches
NCISolver.h
Go to the documentation of this file.
1/*
2 * Copyright University of Reims Champagne-Ardenne
3 * Authors and Contributors: Akilan RAJAMANI, Corentin LEFEBVRE, Johanna KLEIN,
4 * Emmanuel PLUOT, Hugo ROUSSEL, Gaetan RUBEZ, Hassan KHARTABIL,
5 * Jean-Charles BOISSON and Eric HENON
6 * (24/07/2017)
7 * jean-charles.boisson@univ-reims.fr, eric.henon@univ-reims.fr
8 *
9 * This software is a computer program whose purpose is to
10 * detect and quantify interactions from electron density
11 * obtained either internally from promolecular density or
12 * calculated from an input wave function input file. It also
13 * prepares for the visualization of isosurfaces representing
14 * several descriptors (dg) coming from the IGM methodology.
15 *
16 * This software is governed by the CeCILL-C license under French law and
17 * abiding by the rules of distribution of free software. You can use,
18 * modify and/ or redistribute the software under the terms of the CeCILL-C
19 * license as circulated by CEA, CNRS and INRIA at the following URL
20 * "http://www.cecill.info".
21 *
22 * As a counterpart to the access to the source code and rights to copy,
23 * modify and redistribute granted by the license, users are provided only
24 * with a limited warranty and the software's author, the holder of the
25 * economic rights, and the successive licensors have only limited
26 * liability.
27 *
28 * In this respect, the user's attention is drawn to the risks associated
29 * with loading, using, modifying and/or developing or reproducing the
30 * software by the user in light of its specific status of free software,
31 * that may mean that it is complicated to manipulate, and that also
32 * therefore means that it is reserved for developers and experienced
33 * professionals having in-depth computer knowledge. Users are therefore
34 * encouraged to load and test the software's suitability as regards their
35 * requirements in conditions enabling the security of their systems and/or
36 * data to be ensured and, more generally, to use and operate it in the
37 * same conditions as regards security.
38 *
39 * The fact that you are presently reading this means that you have had
40 * knowledge of the CeCILL-C license and that you accept its terms.
41 *
42 * */
43
48
49#ifndef _NCI_SOLVER_H_
50#define _NCI_SOLVER_H_
51
52// Local
53#include <output.h>
54//#include <reader.h>
55#include <Node.h>
56
58#define THREE std::setw(3) << std::setfill(' ')
59
61#define FOUR std::setw(4) << std::setfill(' ')
62
64#define FIVE std::setw(5) << std::setfill(' ')
65
67#define SIX std::setw(6) << std::setfill(' ')
68
70#define SEVEN std::setw(7) << std::setfill(' ')
71
73#define HEIGHT std::setw(8) << std::setfill(' ')
74
76#define NINE std::setw(9) << std::setfill(' ')
77
78
79
85{
86
87 private :
88
89 // ENTRY DATA
90
93
96
97 // MAIN DATA
98
101
104
107
108 // MISCEALENOUS DATA
109
111 struct timeval start, interm, stop;
112
114 // printCurrentState
116
119
121 unsigned long int fullSize;
122
125
127 bool log;
128
130 double total;
131
133 double max;
134
136 unsigned int LI[35][3];
137
139 unsigned int L3[35];
140
143
145 unsigned int nbThreads;
146
148 unsigned int fragNbPrim;
149
150
152 unsigned int* fragPrim; // to store the primitives indexes;
153 // primitives indexes are in the range [0:n-1]
154 // They are stored at the fragPrim index in the range [0:nA+nB-1]
155 // !! nA+nB <= n
156
158 unsigned int* fragAPrim; // to store the primitives indexes of fragment A
159 // primitives indexes are in the range [0:n-1]
160 // They are stored at the fragAPrim index in the range [0:nA-1]
161
163 unsigned int* fragBPrim; // to store the primitives indexes of fragment B
164 // primitives indexes are in the range [0:n-1]
165 // They are stored at the fragBPrim index in the range [0:nB-1]
166
168 unsigned int* allPrim ; // to store the primitives indexes of the whole system
169 // primitives indexes are in the range [0:n-1]
170 // i.e. allPrim[i] = i, simply (to be used for dgscaled option)
171
173 unsigned int* fragAIndexToFragPrimIndex; // input = integer in the range [0:nA-1]
174 // returns an index in the range[0:nA+nB-1]
175
177 unsigned int* fragBIndexToFragPrimIndex; // input = integer in the range [0:nB-1]
178 // returns an index in the range[0:nA+nB-1]
179
181 static double maxElectron;
182
185
188
190 unsigned int nbPrimInA;
191
193 unsigned int nbPrimInB;
194
197 unsigned int* nprimbond;
198
200 double* primRmin;
201
202 // NOT USED
203 // The number of core ATOMIC orbitals (BDA calculations only account for valence orbitals)
204 //unsigned int coreAONumber;
205
206 // NOT USED
207 // The number of core orbitals read from WFX/WFN (BDA calculations only account for valence orbitals)
208 //unsigned int coreMOWFN;
209
212
214 double elecNumb;
215
218
221
224
225 // NOT USED
226 //A boolean to say if the Mos have been reorganized or not by IGMPlot
227 //double MOSorted;
228
229
231 std::vector<int> atomTypeList; // a list of atom types (integer from IGM atom listR in the range 0:SIZE_ATOM-1)
232
234 std::vector<int> atomTypeNb;
235
237 unsigned int** primbond;
238
240 double *maxc;
241
243 unsigned int npri;
244
246 std::vector<moleculeOrbital> molecularOrbitals;
247
248
249 // a vector containing properties of all found critical points
250 std::vector<criticalpoint> cpList;
251
252 // a boolean to tell if the IGM promolecular determination of seeds is activated
253 bool IGMSeeds;
254
258 void calcprops_wfn(); // IGM dg dgInter dgIntra CUBE calculations
259
263 void calcprops_wfn_cyl();
264
269 void setDensityMatrixForPauli(double** D);
270
274 void setPrimList();
275
279 void setIBSIPrimList();
280
284 void setCOREMOandELEC(); // count the number of core MOs and electrons
285
289 void setChemFormula();
290
306 void IGM(double** chi, double** phi,
307 double* dx, double* dy, double* dz, double** gradPrim,
308 double &deltagIntraCURRENT, double &deltagInterCURRENT, double &qgIntraCURRENT, double &qgInterCURRENT);
309
310
329void gradAtomHirsh(bool fullRHOGRAD, unsigned int threadID, double rho, double gradrho[3],
330 double*** dx, double*** dy, double*** dz, int ip,
331 double* rhoFree, double *pregradFree, double** gradAtom, double* rhoAtom);
332
340void gradAtomGBP(unsigned int nbPrim, unsigned int* prim, double** gradPrim, double** gradAtom);
341
342
343
352void IGMH(double** gradAtom, double &deltagIntraCURRENT, double &deltagInterCURRENT,double &qgIntraCURRENT, double &qgInterCURRENT);
353
360void IGMBDA(double *gradAtom1QM, double *gradAtom2QM, double &bda);
361
368void rhoFree__pregradFree(double* dist, double* rhoFree, double* pregradFree);
369
377void rhoFree__pregradFree_sumSquareBRho(double* dist, double* rhoFree, double* pregradFree, double* sumSquareBRho);
378
379
384double rhopromolFRAG1FRAG2(double* rhoFree);
385
386// rhoQM : input --> the true rho obtained by means of QM for the WHOLE system
387// rhoFree : input --> the free promol ED for every atom (calculated previously for every atom in the WHOLE system)
388// RETURN : estimated QM ED for the FRAG1+FRAG2 system from HIRSHFELD partition
389
390
391
397double rhoHIRSHFRAG1FRAG2(double rhoQM, double* rhoFree);
398
399
400 public :
401
402// Define a GTO function type f(r, alpha) (to sort non-nul atomic orbital(with exponent alpha) at a given
403// distance r of the current grid node)
404using FunctionType = double (*)(double, double);
405static FunctionType functions[];
406static FunctionType derivatives[];
407
408
410typedef struct point
411 {
413 double x,y,z;
415
421 NCISolver(const char* parameterFileName, bool logParam = true);
422
423
427 ~NCISolver();
428
432 void solve();
433
440 void lineProcessing(int posY, int posZ, int cubePosYZ);
441
445 void output();
446
451 std::string score();
452
458 unsigned int* index0(const unsigned int i);
459
465 unsigned int index1(const unsigned int i);
466
470 void initializeLI();
471
475 void initializeL3();
476
477
485 void getLambdaOfHessian(double** hessian,double* eigenValues, double** eigenVect);
486
491 double getPredictedSize();
492
497 bool getRuninfostate();
498
507void calcQMAO(unsigned int nbPrim, unsigned int* prim, double* d2, double* dx, double* dy, double* dz, double** chi);
508
509
520void calcQMPROP(unsigned int nbPrim, unsigned int* prim, double** chi, double** phi, double &rho, double* grad, double** hess);
521
525void findmaxc();
526
538bool NewtonRaphson(unsigned int nbPrim, unsigned int* prim, double rcurr[3], double L[3], double &G, double &rhocp, double &gradrhocp);
539
545int cptype(double L123[3]);
546
557void basischange(unsigned int nbPrim, unsigned int* prim, double matrixHessianEigVec[3][3],
558 double **chi, double* dx, double* dy, double* dz, double **phi, double gg[3]);
559
568void IGMPRO(double *dx, double *dy, double *dz, double &deltag, double &qg);
569
577void gradPRO(double *dx, double *dy, double *dz, double **gradAtom);
578
579
594void rhogradlap(unsigned int atIndex1, unsigned int atIndex2, double* rhoFree, double* pregradFree,
595 double* sumSquareBRho, double* dx, double* dy, double* dz, double* dist,
596 double &rhopro, double &normgradpro, double &lap);
597// ========= PROMOLECULAR MODE ONLY ================ //
598
599// functions used to limit the number of GTO to be calculated
600// (speedup test, Hugo ROUSSEL work)
601// Function 0 : for GTO(0,0,0) <-> s
606static double func0(double r, double alpha);
607
608
609
610// Function 1 : for GTO(1,0,0) ot (0,1,0) or (0,0,1) <-> p
615static double func1(double r, double alpha);
616
617
618
619
620// Function 2 : for GTO(2,0,0) or (1,1,0) or ... <-> d
625static double func2(double r, double alpha);
626
627
628// Function 3 : for GTO(3,0,0) or (1,1,1) or ... <-> f
633static double func3(double r, double alpha);
634
635// Function 4 : for GTO(4,0,0) or (2,2,0) or ... <-> g
640static double func4(double r, double alpha);
641
642
643
648double findRmin(
649 int GTOType, // integer [0:4] coding for s,p,d,f or g INPUT
650 double alpha, // alpha exponent of the GTO INPUT
651 double tolerance = 0.010, // in bohr, convergence threshold INPUT
652 int maxIterations = 100, // Maximum number of iteration INPUT
653 double initialGuess = 0.01 // initial guess of rmin (bohr) for which f(r,alpha) = 0.999 INPUT
654);
655
660
661#endif
662
663
664};
Worker for one Node of the NCI grid.
double getPredictedSize()
Return the predicted size of the ouput in MB.
Definition NCISolver.cpp:7814
void IGMPRO(double *dx, double *dy, double *dz, double &deltag, double &qg)
compute the deltag descriptor at current point
Definition NCISolver.cpp:10247
double * primRmin
Array to store the minimum radius below which a primitive has to be computed (scalProd speedup test)
Definition NCISolver.h:200
void solve()
Solves the problem.
Definition NCISolver.cpp:6641
double MOoccupied
The number of MO occupied with occupancy >=0.1.
Definition NCISolver.h:223
void rhoFree__pregradFree(double *dist, double *rhoFree, double *pregradFree)
Compute quantities needed to compute next rho and gradrho at the promol level.
Definition NCISolver.cpp:8930
unsigned int * nprimbond
Definition NCISolver.h:197
void IGMBDA(double *gradAtom1QM, double *gradAtom2QM, double &bda)
Compute the deltagInter quantity between two given atoms of a bond and the associated BDA based on HI...
Definition NCISolver.cpp:9309
bool runinfostate
runinfo file state (timing info for users printed by
Definition NCISolver.h:115
int cptype(double L123[3])
determine the signature of the critical point (algebrix sum of the signs of the three eigenvalues)
Definition NCISolver.cpp:10133
unsigned int * fragAIndexToFragPrimIndex
An array converting indexes of fragAPrim to indexes of fragPrim.
Definition NCISolver.h:173
void IGMH(double **gradAtom, double &deltagIntraCURRENT, double &deltagInterCURRENT, double &qgIntraCURRENT, double &qgInterCURRENT)
Compute the deltagIntra and deltagInter quantities.
Definition NCISolver.cpp:9187
~NCISolver()
Destructor.
Definition NCISolver.cpp:1867
param_t params
Problem's parameters.
Definition NCISolver.h:92
unsigned int LI[35][3]
For index0.
Definition NCISolver.h:136
void findmaxc()
find the maximum value for primitive coefficients accross the MOs
Definition NCISolver.cpp:9852
bool * inMoleculeA
Array to store the membership of each primitives of Fragment A.
Definition NCISolver.h:184
std::string score()
Return a string with the interaction's score.
Definition NCISolver.cpp:7617
Node ** nodes
A program node used by threads.
Definition NCISolver.h:103
static double maxElectron
The accuracy required for atomic orbitals calculation (speed-up procedure by J. Pilme)
Definition NCISolver.h:181
std::vector< int > atomTypeNb
List of the number of atoms for each atom type.
Definition NCISolver.h:234
double rhopromolFRAG1FRAG2(double *rhoFree)
Compute the free promol ED for FRAG1+FRAG2 system.
Definition NCISolver.cpp:8963
unsigned int * fragBPrim
the list (array) of primitives of fragment B
Definition NCISolver.h:163
unsigned int * allPrim
the list (array) of ALL primitives
Definition NCISolver.h:168
NCISolver(const char *parameterFileName, bool logParam=true)
Main constructor.
Definition NCISolver.cpp:58
unsigned int nbPrimInA
the number of primitives in fragment A (Cube mode)
Definition NCISolver.h:190
unsigned int * index0(const unsigned int i)
get the triplet (i,j,k) defining a GTO (or STO) from a code in the range [0:34]
Definition NCISolver.cpp:7624
int currentPercentage
Current completion in percent.
Definition NCISolver.h:124
void calcQMAO(unsigned int nbPrim, unsigned int *prim, double *d2, double *dx, double *dy, double *dz, double **chi)
Compute primitives from WFX/WFN data at current grid point.
Definition NCISolver.cpp:9326
double rhoHIRSHFRAG1FRAG2(double rhoQM, double *rhoFree)
estimated QM ED for the FRAG1+FRAG2 system from HIRSHFELD partition
Definition NCISolver.cpp:9111
Results * results
Results' values.
Definition NCISolver.h:100
void gradPRO(double *dx, double *dy, double *dz, double **gradAtom)
compute the promolecular ED gradient at current point
Definition NCISolver.cpp:10314
void calcQMPROP(unsigned int nbPrim, unsigned int *prim, double **chi, double **phi, double &rho, double *grad, double **hess)
Compute primitives and MOs from WFX/WFN data at current grid point.
Definition NCISolver.cpp:9739
void initializeL3()
Tool procedure that initializes L3 array (s,p,d,f,g GTOs)
Definition NCISolver.cpp:7738
struct NCISolver::point point
structure to save 3D point (triplet of doubles)
double coreElecWFN
The number of core electrons (BDA calculations only account for valence orbitals)
Definition NCISolver.h:220
axis_t * posGrid
Array with current grid position.
Definition NCISolver.h:106
unsigned int COREHeavyAtomWarning
Boolean flag to indicate if one atom exceeds the upper atomic number limit for BDA treatment.
Definition NCISolver.h:211
void rhoFree__pregradFree_sumSquareBRho(double *dist, double *rhoFree, double *pregradFree, double *sumSquareBRho)
Compute quantities needed to compute next rho and gradrho at the promol level.
Definition NCISolver.cpp:8891
void basischange(unsigned int nbPrim, unsigned int *prim, double matrixHessianEigVec[3][3], double **chi, double *dx, double *dy, double *dz, double **phi, double gg[3])
changes the basis for several vectors of the QM calculation
Definition NCISolver.cpp:10157
unsigned int nbThreads
The number of threads for the current run.
Definition NCISolver.h:145
void setChemFormula()
establish the chemical formula of the studied system
Definition NCISolver.cpp:8351
double total
Total dg inter.
Definition NCISolver.h:130
unsigned int * fragPrim
the list (array) of primitives to be considered in the IGM QM calculations
Definition NCISolver.h:152
double elecNumb
The total theoretical number of electrons.
Definition NCISolver.h:214
bool NewtonRaphson(unsigned int nbPrim, unsigned int *prim, double rcurr[3], double L[3], double &G, double &rhocp, double &gradrhocp)
performs a Newton-Raphson procedure to get the Critical point closest to the guess passed in paramete...
Definition NCISolver.cpp:9872
void setCOREMOandELEC()
count the number of core MOs and electrons
Definition NCISolver.cpp:8144
ProgData * data
Problem's data.
Definition NCISolver.h:95
void calcprops_wfn_cyl()
IBSI + BDA calculations, no cube generated.
Definition NCISolver.cpp:5015
int zyAxisSize
Number of row in the grid.
Definition NCISolver.h:118
unsigned int * fragAPrim
the list (array) of primitives of fragment A
Definition NCISolver.h:158
unsigned int * fragBIndexToFragPrimIndex
An array converting indexes of fragAPrim to indexes of fragPrim.
Definition NCISolver.h:177
void getLambdaOfHessian(double **hessian, double *eigenValues, double **eigenVect)
Using the eig3 implementation of diagonalisation.
Definition NCISolver.cpp:7784
std::vector< moleculeOrbital > molecularOrbitals
the set of expressions of MOs read from the WFN/WFX
Definition NCISolver.h:246
void findRminAllPrim()
function that compute the minimum radius for all primitives
unsigned int L3[35]
For index1.
Definition NCISolver.h:139
double * maxc
the maximum value of primitive coefficient accross the MOs
Definition NCISolver.h:240
void setIBSIPrimList()
filtering primitives according to the cutoff radius around atom pairs
Definition NCISolver.cpp:8007
unsigned int npri
the size of the basis set read from WFN
Definition NCISolver.h:243
void lineProcessing(int posY, int posZ, int cubePosYZ)
Process all the information for an entire line of the grid.
Definition NCISolver.cpp:7438
void gradAtomGBP(unsigned int nbPrim, unsigned int *prim, double **gradPrim, double **gradAtom)
Compute the atomic contribution to the ED gradient according to the GBP partition.
Definition NCISolver.cpp:9143
unsigned int nbPrimInB
the number of primitives in fragment B (Cube mode)
Definition NCISolver.h:193
void setDensityMatrixForPauli(double **D)
compute matrix density elements from wave function
Definition NCISolver.cpp:7825
struct timeval start interm stop
start and stop time to provide timing info for the user
Definition NCISolver.h:111
unsigned long int fullSize
Number of nodes in the grid.
Definition NCISolver.h:121
unsigned int index1(const unsigned int i)
identify which kind of AO it is
Definition NCISolver.cpp:7642
double WFNelecNumb
The real number of electrons reac from WFN.
Definition NCISolver.h:217
unsigned int ** primbond
the list (array) of primitives to be considered in the IBSI calculation
Definition NCISolver.h:237
void initializeLI()
Tool procedure that initializes LI array (36 GTO types)
Definition NCISolver.cpp:7664
std::vector< int > atomTypeList
List of atom types present in the system.
Definition NCISolver.h:231
void setPrimList()
filtering primitives according to fragment definitions
Definition NCISolver.cpp:7869
void IGM(double **chi, double **phi, double *dx, double *dy, double *dz, double **gradPrim, double &deltagIntraCURRENT, double &deltagInterCURRENT, double &qgIntraCURRENT, double &qgInterCURRENT)
compute the deltagIntra and deltagIntrer in cartesian coordinates (not cylindrical)
Definition NCISolver.cpp:8386
bool * inMoleculeB
Array to store the membership of each primitives of Fragment B.
Definition NCISolver.h:187
unsigned int fragNbPrim
The number of primitives to be considered in the IGM WFN( WFX) calculations (cube mode)
Definition NCISolver.h:148
void calcprops_wfn()
IGM dg dgInter dgIntra CUBE calculations.
Definition NCISolver.cpp:1897
bool getRuninfostate()
Return the runinfo state (true = built by printCurrentState)
Definition NCISolver.cpp:7819
double max
Maximum percent value.
Definition NCISolver.h:133
void gradAtomHirsh(bool fullRHOGRAD, unsigned int threadID, double rho, double gradrho[3], double ***dx, double ***dy, double ***dz, int ip, double *rhoFree, double *pregradFree, double **gradAtom, double *rhoAtom)
Compute the atomic contribution to the ED gradient according to the Hirschfeld partition.
Definition NCISolver.cpp:8716
void output()
Outputs according the output type given in the parameters.
Definition NCISolver.cpp:7467
bool log
Activates logs.
Definition NCISolver.h:127
double predictedSize
Predicted size of the data on HDD.
Definition NCISolver.h:142
double findRmin(int GTOType, double alpha, double tolerance=0.010, int maxIterations=100, double initialGuess=0.01)
function that takes as input the type (s,p,d,f,g) of GTO and its alpha exponent and returns a radius ...
Definition NCISolver.cpp:10364
Class used to handle every node data storage and processing ; implemented for multiprocessing impleme...
Definition Node.h:59
Class designed to store the program's main dataprovide some utilities concerning those.
Definition ProgData.h:65
Class designed to store all results that can be used as output by the program.
Definition Results.h:67
Files writing functions.
structure to save 3D point (triplet of doubles)
Definition NCISolver.h:411
double x
x,y,z coordinates
Definition NCISolver.h:413
Structure describing one axis.
Definition general.h:694
Structure used for the storage of the parameters read from the parameter file param....
Definition general.h:560