NCISolver.h

Go to the documentation of this file.

/*
 * Copyright University of Reims Champagne-Ardenne
 * Authors and Contributors: Akilan RAJAMANI, Corentin LEFEBVRE, Johanna KLEIN,
 *                           Emmanuel PLUOT, Gaetan RUBEZ, Hassan KHARTABIL,
 *                           Jean-Charles BOISSON and Eric HENON
 * (24/07/2017)
 * jean-charles.boisson@univ-reims.fr, eric.henon@univ-reims.fr
 *
 * This software is a computer program whose purpose is to
 * detect and quantify interactions from electron density
 * obtained either internally from promolecular density or
 * calculated from an input wave function input file. It also
 * prepares for the visualization of isosurfaces representing
 * several descriptors (dg) coming from the IGM methodology.
 *
 * This software is governed by the CeCILL-C license under French law and
 * abiding by the rules of distribution of free software.  You can  use,
 * modify and/ or redistribute the software under the terms of the CeCILL-C
 * license as circulated by CEA, CNRS and INRIA at the following URL
 * "http://www.cecill.info".
 *
 * As a counterpart to the access to the source code and  rights to copy,
 * modify and redistribute granted by the license, users are provided only
 * with a limited warranty  and the software's author,  the holder of the
 * economic rights,  and the successive licensors  have only  limited
 * liability.
 *
 * In this respect, the user's attention is drawn to the risks associated
 * with loading,  using,  modifying and/or developing or reproducing the
 * software by the user in light of its specific status of free software,
 * that may mean  that it is complicated to manipulate,  and  that  also
 * therefore means  that it is reserved for developers  and  experienced
 * professionals having in-depth computer knowledge. Users are therefore
 * encouraged to load and test the software's suitability as regards their
 * requirements in conditions enabling the security of their systems and/or
 * data to be ensured and,  more generally, to use and operate it in the
 * same conditions as regards security.
 *
 * The fact that you are presently reading this means that you have had
 * knowledge of the CeCILL-C license and that you accept its terms.
 *
 * */

 
#ifndef _NCI_SOLVER_H_
#define _NCI_SOLVER_H_
 
// Local
#include <output.h>
//#include <reader.h>
#include <Node.h>
#include <map>

#define THREE std::setw(3) << std::setfill(' ')

#define FOUR std::setw(4) << std::setfill(' ')

#define FIVE std::setw(5) << std::setfill(' ')

#define SIX  std::setw(6) << std::setfill(' ')

#define SEVEN  std::setw(7) << std::setfill(' ')

#define HEIGHT std::setw(8) << std::setfill(' ')

#define NINE std::setw(9) << std::setfill(' ')
 
 

class NCISolver
{
    
 private : 
 
  // ENTRY DATA

  param_t params;

  ProgData * data;
        
  // MAIN DATA

  Results * results;

  Node **nodes;

  axis_t *posGrid;
        
  // MISCEALENOUS DATA
 
  struct timeval start, interm, stop;

  //  printCurrentState
  bool runinfostate; 

  int zyAxisSize;

  unsigned long int fullSize;

  int currentPercentage;

  bool log;

  double total;

  double max;

  unsigned int LI[35][3];

  unsigned int L3[35];

  double predictedSize;
  
  unsigned int nbThreads;

  unsigned int  fragNbPrim;
 
 
  unsigned int* fragPrim; // to store the primitives indexes; 
                         // primitives indexes are in the range [0:n-1]
                         // They are stored at the fragPrim index in the range [0:nA+nB-1]
                         // !! nA+nB <= n

  unsigned int* fragAPrim; // to store the primitives indexes of fragment A
                         // primitives indexes are in the range [0:n-1]
                         // They are stored at the fragAPrim index in the range [0:nA-1]

  unsigned int* fragBPrim; // to store the primitives indexes of fragment B 
                         // primitives indexes are in the range [0:n-1]
                         // They are stored at the fragBPrim index in the range [0:nB-1]

  unsigned int* allPrim  ; // to store the primitives indexes of the whole system
                         // primitives indexes are in the range [0:n-1]
                         // i.e. allPrim[i] = i, simply (to be used for dgscaled option)

  unsigned int* fragAIndexToFragPrimIndex;  // input = integer in the range [0:nA-1]
                                            // returns an index in the range[0:nA+nB-1]

  unsigned int* fragBIndexToFragPrimIndex;  // input = integer in the range [0:nB-1]
                                            // returns an index in the range[0:nA+nB-1]

  static double maxElectron;

  bool* inMoleculeA;

  bool* inMoleculeB;

  unsigned int nbPrimInA; 

  unsigned int nbPrimInB;

  unsigned int*   nprimbond;

  double* primRmin;     
 
  // NOT USED
  // The number of core ATOMIC orbitals (BDA calculations only account for valence orbitals)
  //unsigned int coreAONumber;
 
  // NOT USED
  // The number of core orbitals read from WFX/WFN (BDA calculations only account for valence orbitals)
  //unsigned int coreMOWFN;

  unsigned int COREHeavyAtomWarning;
  
  double elecNumb;

  double WFNelecNumb;

  double coreElecWFN;

  double MOoccupied;
 
  // NOT USED
  //A boolean to say if the Mos have been reorganized or not by IGMPlot                
  //double MOSorted;
 

  std::vector<int> atomTypeList; // a list of atom types (integer from IGM atom listR in the range 0:SIZE_ATOM-1)

  std::vector<int> atomTypeNb;

  unsigned int**  primbond;

  double *maxc;

  unsigned int  npri; 

  std::vector<moleculeOrbital> molecularOrbitals;
 
 
  // a vector containing properties of all found critical points
  std::vector<criticalpoint>  cpList;
 
  // a boolean to tell if the IGM promolecular determination of seeds is activated
  bool IGMSeeds;

  void calcprops_wfn(); // IGM dg dgInter dgIntra CUBE calculations

  void calcprops_wfn_cyl(); 

  void setDensityMatrixForPauli(double** D);

  void setPrimList(); 

  void setIBSIPrimList();

  void setCOREMOandELEC(); // count the number of core MOs and electrons 

  void setChemFormula();

  void IGM(double** chi, double** phi,
           double* dx, double* dy, double* dz, double** gradPrim,
           double &deltagIntraCURRENT, double &deltagInterCURRENT, double &qgIntraCURRENT, double &qgInterCURRENT);
 

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);

void gradAtomGBP(unsigned int nbPrim, unsigned int* prim, double** gradPrim, double** gradAtom);  
 
 

void IGMH(double** gradAtom, double &deltagIntraCURRENT, double &deltagInterCURRENT,double &qgIntraCURRENT, double &qgInterCURRENT);

void IGMBDA(double *gradAtom1QM, double *gradAtom2QM, double &bda);

void rhoFree__pregradFree(double* dist, double* rhoFree, double* pregradFree);

// P R O M O L E C U L A R
void rhoFree_gradFree_lapFree(double* dist, double* dx, double* dy, double* dz,
                                                   double* rhoFree, double** gradFree, double* lapFree);
 

void rhoFree__pregradFree_sumSquareBRho(double* dist, double* rhoFree, double* pregradFree, double* sumSquareBRho);
 

double rhopromolFRAG1FRAG2(double* rhoFree);
 
// rhoQM    : input --> the true rho obtained by means of QM for the WHOLE system 
// rhoFree  : input --> the free promol ED for every atom (calculated previously for every atom in the WHOLE system) 
// RETURN   : estimated QM ED for the FRAG1+FRAG2 system from HIRSHFELD partition
 
 

double rhoHIRSHFRAG1FRAG2(double rhoQM, double* rhoFree);
 
  
 public : 
 
// Define a GTO function type f(r, alpha) (to sort non-nul atomic orbital(with exponent alpha) at a given 
// distance r of the current grid node)
using FunctionType = double (*)(double, double);
static FunctionType functions[];
static FunctionType derivatives[];
 

typedef struct point          
  {
    double x,y,z;
  } point;
 

typedef struct basinProperties 
  {
    double electrons;
    double volume;
    double KinExcess;
    
    // Constructor to facilitate initialization     
    basinProperties(double e = 0.0, double v = 0.0, double k = 0.0) : electrons(e), volume(v), KinExcess(k) {}
  } basinProperties;

  NCISolver(const char* parameterFileName, bool logParam = true);
 
        
  ~NCISolver();

  void solve();
              
  void lineProcessing(int posY, int posZ, int cubePosYZ);
        
  void output();
        
  std::string score();
        
  unsigned int* index0(const unsigned int i);

  unsigned int index1(const unsigned int i);

  void initializeLI();

  void initializeL3();
 

  void getLambdaOfHessian(double** hessian,double* eigenValues, double** eigenVect);
  
  double getPredictedSize();

  bool getRuninfostate();

void calcQMAO(unsigned int nbPrim, unsigned int* prim, double* d2, double* dx, double* dy, double* dz, double** chi);
 

void calcQMPROP(unsigned int nbPrim, unsigned int* prim, double** chi,  double** phi, double &rho, double* grad, double** hess);

void findmaxc();

bool NewtonRaphson(unsigned int nbPrim, unsigned int* prim, double rcurr[3], double L[3], double &G, double &rhocp, double &gradrhocp);

int cptype(double L123[3]); 

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]);

void IGMPRO(double *dx, double *dy, double *dz, double &deltag, double &qg);

void dgdgAtPRO(double *dx, double *dy, double *dz, double &deltag, std::vector<double>& dgAt);
 

void gradPRO(double *dx, double *dy, double *dz, double **gradAtom);
 
 

// ========= PROMOLECULAR G2 MODE ONLY ================ //
void rhogradlap(unsigned int at1, unsigned int at2, double* rhoFree, double** gradFree, double* lapFree,
                           double &rhopro, double* gradpro, double &normgradpro, double &lap);
 
// ========= PROMOLECULAR MODE ONLY ================ //
 
// functions used to limit the number of GTO to be calculated
// (speedup test, Hugo ROUSSEL work)
// Function 0 : for GTO(0,0,0)                        <-> s
static double func0(double r, double alpha); 
 
 
 
// Function 1 : for GTO(1,0,0) ot (0,1,0) or (0,0,1)  <-> p
static double func1(double r, double alpha);
 
 
 
 
// Function 2 : for GTO(2,0,0) or (1,1,0) or ...      <-> d
static double func2(double r, double alpha);
 
 
// Function 3 : for GTO(3,0,0) or (1,1,1) or ...      <-> f
static double func3(double r, double alpha);
 
// Function 4 : for GTO(4,0,0) or (2,2,0) or ...      <-> g
static double func4(double r, double alpha);
 
 

double findRmin(
    int GTOType,               // integer [0:4] coding for s,p,d,f or g     INPUT
    double alpha,              // alpha exponent of the GTO                 INPUT
    double tolerance = 0.010,  // in bohr, convergence threshold             INPUT
    int maxIterations = 100,   // Maximum number of iteration               INPUT
    double initialGuess = 0.01 // initial guess of rmin (bohr) for which f(r,alpha) = 0.999 INPUT
);

void findRminAllPrim();

double encodeBasin(int a, int b);

std::vector<int> decodeBasin(double n);

std::map<std::pair<int, int>, basinProperties> computeBasinProperties(double*** ELFbasin, double*** KE, double*** rho, const unsigned int dimX, const unsigned int dimY, const unsigned int dimZ);

std::pair<size_t, size_t> find2max(const std::vector<double>& dgAt);

// ========= PROMOLECULAR MODE ONLY ================ //
double gradnormIGM(double* pregradFree,
                 double* dx, double* dy, double* dz);
 
#endif
 
 
};