Line data    Source code

       1             : /*----------------------------------------------------------------------------*/
       2             : /*  CP2K: A general program to perform molecular dynamics simulations         */
       3             : /*  Copyright 2000-2024 CP2K developers group <https://cp2k.org>              */
       4             : /*                                                                            */
       5             : /*  SPDX-License-Identifier: BSD-3-Clause                                     */
       6             : /*----------------------------------------------------------------------------*/
       7             : 
       8             : #include <assert.h>
       9             : #include <float.h>
      10             : #include <limits.h>
      11             : #include <math.h>
      12             : #include <stdio.h>
      13             : #include <stdlib.h>
      14             : #include <string.h>
      15             : 
      16             : #define GRID_DO_COLLOCATE 1
      17             : #include "../common/grid_common.h"
      18             : #include "grid_cpu_collint.h"
      19             : #include "grid_cpu_collocate.h"
      20             : #include "grid_cpu_integrate.h"
      21             : #include "grid_cpu_prepare_pab.h"
      22             : 
      23             : /*******************************************************************************
      24             :  * \brief Writes the given arguments into a .task file.
      25             :  *        See grid_replay.h for details.
      26             :  * \author Ole Schuett
      27             :  ******************************************************************************/
      28           0 : static void write_task_file(
      29             :     const bool orthorhombic, const int border_mask, const enum grid_func func,
      30             :     const int la_max, const int la_min, const int lb_max, const int lb_min,
      31             :     const double zeta, const double zetb, const double rscale,
      32             :     const double dh[3][3], const double dh_inv[3][3], const double ra[3],
      33             :     const double rab[3], const int npts_global[3], const int npts_local[3],
      34             :     const int shift_local[3], const int border_width[3], const double radius,
      35             :     const int o1, const int o2, const int n1, const int n2,
      36           0 :     const double pab[n2][n1], const double *grid) {
      37             : 
      38           0 :   static int counter = 1;
      39           0 :   int my_number;
      40             : 
      41           0 : #pragma omp critical
      42           0 :   my_number = counter++;
      43             : 
      44           0 :   char filename[100];
      45           0 :   snprintf(filename, sizeof(filename), "grid_collocate_%05i.task", my_number);
      46             : 
      47           0 :   const int D = DBL_DECIMAL_DIG;
      48           0 :   FILE *fp = fopen(filename, "w+");
      49           0 :   fprintf(fp, "#Grid task v10\n");
      50           0 :   fprintf(fp, "orthorhombic %i\n", orthorhombic);
      51           0 :   fprintf(fp, "border_mask %i\n", border_mask);
      52           0 :   fprintf(fp, "func %i\n", func);
      53           0 :   fprintf(fp, "la_max %i\n", la_max);
      54           0 :   fprintf(fp, "la_min %i\n", la_min);
      55           0 :   fprintf(fp, "lb_max %i\n", lb_max);
      56           0 :   fprintf(fp, "lb_min %i\n", lb_min);
      57           0 :   fprintf(fp, "zeta %.*e\n", D, zeta);
      58           0 :   fprintf(fp, "zetb %.*e\n", D, zetb);
      59           0 :   fprintf(fp, "rscale %.*e\n", D, rscale);
      60           0 :   for (int i = 0; i < 3; i++)
      61           0 :     fprintf(fp, "dh %i %.*e %.*e %.*e\n", i, D, dh[i][0], D, dh[i][1], D,
      62           0 :             dh[i][2]);
      63           0 :   for (int i = 0; i < 3; i++)
      64           0 :     fprintf(fp, "dh_inv %i %.*e %.*e %.*e\n", i, D, dh_inv[i][0], D,
      65           0 :             dh_inv[i][1], D, dh_inv[i][2]);
      66           0 :   fprintf(fp, "ra %.*e %.*e %.*e\n", D, ra[0], D, ra[1], D, ra[2]);
      67           0 :   fprintf(fp, "rab %.*e %.*e %.*e\n", D, rab[0], D, rab[1], D, rab[2]);
      68           0 :   fprintf(fp, "npts_global %i %i %i\n", npts_global[0], npts_global[1],
      69             :           npts_global[2]);
      70           0 :   fprintf(fp, "npts_local %i %i %i\n", npts_local[0], npts_local[1],
      71             :           npts_local[2]);
      72           0 :   fprintf(fp, "shift_local %i %i %i\n", shift_local[0], shift_local[1],
      73             :           shift_local[2]);
      74           0 :   fprintf(fp, "border_width %i %i %i\n", border_width[0], border_width[1],
      75             :           border_width[2]);
      76           0 :   fprintf(fp, "radius %.*e\n", D, radius);
      77           0 :   fprintf(fp, "o1 %i\n", o1);
      78           0 :   fprintf(fp, "o2 %i\n", o2);
      79           0 :   fprintf(fp, "n1 %i\n", n1);
      80           0 :   fprintf(fp, "n2 %i\n", n2);
      81             : 
      82           0 :   for (int i = 0; i < n2; i++) {
      83           0 :     for (int j = 0; j < n1; j++) {
      84           0 :       fprintf(fp, "pab %i %i %.*e\n", i, j, D, pab[i][j]);
      85             :     }
      86             :   }
      87             : 
      88           0 :   const int npts_local_total = npts_local[0] * npts_local[1] * npts_local[2];
      89             : 
      90           0 :   int ngrid_nonzero = 0;
      91           0 :   for (int i = 0; i < npts_local_total; i++) {
      92           0 :     if (grid[i] != 0.0) {
      93           0 :       ngrid_nonzero++;
      94             :     }
      95             :   }
      96           0 :   fprintf(fp, "ngrid_nonzero %i\n", ngrid_nonzero);
      97             : 
      98           0 :   for (int k = 0; k < npts_local[2]; k++) {
      99           0 :     for (int j = 0; j < npts_local[1]; j++) {
     100           0 :       for (int i = 0; i < npts_local[0]; i++) {
     101           0 :         const double val =
     102           0 :             grid[k * npts_local[1] * npts_local[0] + j * npts_local[0] + i];
     103           0 :         if (val != 0.0) {
     104           0 :           fprintf(fp, "grid %i %i %i %.*e\n", i, j, k, D, val);
     105             :         }
     106             :       }
     107             :     }
     108             :   }
     109             : 
     110           0 :   double hab[n2][n1];
     111           0 :   double forces[2][3] = {0};
     112           0 :   double virials[2][3][3] = {0};
     113           0 :   memset(hab, 0, n2 * n1 * sizeof(double));
     114             : 
     115           0 :   const bool compute_tau = (func == GRID_FUNC_DADB);
     116             : 
     117           0 :   grid_cpu_integrate_pgf_product(orthorhombic, compute_tau, border_mask, la_max,
     118             :                                  la_min, lb_max, lb_min, zeta, zetb, dh, dh_inv,
     119             :                                  ra, rab, npts_global, npts_local, shift_local,
     120             :                                  border_width, radius, o1, o2, n1, n2, grid,
     121             :                                  hab, pab, forces, virials, NULL, NULL, NULL);
     122             : 
     123           0 :   for (int i = o2; i < ncoset(lb_max) + o2; i++) {
     124           0 :     for (int j = o1; j < ncoset(la_max) + o1; j++) {
     125           0 :       fprintf(fp, "hab %i %i %.*e\n", i, j, D, hab[i][j]);
     126             :     }
     127             :   }
     128           0 :   fprintf(fp, "force_a %.*e %.*e %.*e\n", D, forces[0][0], D, forces[0][1], D,
     129             :           forces[0][2]);
     130           0 :   fprintf(fp, "force_b %.*e %.*e %.*e\n", D, forces[1][0], D, forces[1][1], D,
     131             :           forces[1][2]);
     132           0 :   for (int i = 0; i < 3; i++)
     133           0 :     fprintf(fp, "virial %i %.*e %.*e %.*e\n", i, D,
     134           0 :             virials[0][i][0] + virials[1][i][0], D,
     135           0 :             virials[0][i][1] + virials[1][i][1], D,
     136           0 :             virials[0][i][2] + virials[1][i][2]);
     137             : 
     138           0 :   fprintf(fp, "#THE_END\n");
     139           0 :   fclose(fp);
     140           0 :   printf("Wrote %s\n", filename);
     141           0 : }
     142             : 
     143             : /*******************************************************************************
     144             :  * \brief Collocates a single product of primitiv Gaussians.
     145             :  *        See grid_cpu_collocate.h for details.
     146             :  * \author Ole Schuett
     147             :  ******************************************************************************/
     148    86290985 : static void collocate_internal(
     149             :     const bool orthorhombic, const int border_mask, const enum grid_func func,
     150             :     const int la_max, const int la_min, const int lb_max, const int lb_min,
     151             :     const double zeta, const double zetb, const double rscale,
     152             :     const double dh[3][3], const double dh_inv[3][3], const double ra[3],
     153             :     const double rab[3], const int npts_global[3], const int npts_local[3],
     154             :     const int shift_local[3], const int border_width[3], const double radius,
     155             :     const int o1, const int o2, const int n1, const int n2,
     156    86290985 :     const double pab[n2][n1], double *grid) {
     157             : 
     158    86290985 :   int la_min_diff, la_max_diff, lb_min_diff, lb_max_diff;
     159    86290985 :   grid_cpu_prepare_get_ldiffs(func, &la_min_diff, &la_max_diff, &lb_min_diff,
     160             :                               &lb_max_diff);
     161             : 
     162    86290985 :   const int la_min_cab = imax(la_min + la_min_diff, 0);
     163    86290985 :   const int lb_min_cab = imax(lb_min + lb_min_diff, 0);
     164    86290985 :   const int la_max_cab = la_max + la_max_diff;
     165    86290985 :   const int lb_max_cab = lb_max + lb_max_diff;
     166    86290985 :   const int n1_cab = ncoset(la_max_cab);
     167    86290985 :   const int n2_cab = ncoset(lb_max_cab);
     168             : 
     169    86290985 :   const size_t cab_size = n2_cab * n1_cab;
     170    86290985 :   double cab[cab_size];
     171    86290985 :   memset(cab, 0, cab_size * sizeof(double));
     172             : 
     173    86290985 :   grid_cpu_prepare_pab(func, o1, o2, la_max, la_min, lb_max, lb_min, zeta, zetb,
     174    86290985 :                        n1, n2, pab, n1_cab, n2_cab, (double(*)[n1_cab])cab);
     175    86290985 :   cab_to_grid(orthorhombic, border_mask, la_max_cab, la_min_cab, lb_max_cab,
     176             :               lb_min_cab, zeta, zetb, rscale, dh, dh_inv, ra, rab, npts_global,
     177             :               npts_local, shift_local, border_width, radius, cab, grid);
     178    86290985 : }
     179             : 
     180             : /*******************************************************************************
     181             :  * \brief Public entry point. A thin wrapper with the only purpose of calling
     182             :  *        write_task_file when DUMP_TASKS = true.
     183             :  * \author Ole Schuett
     184             :  ******************************************************************************/
     185    86290985 : void grid_cpu_collocate_pgf_product(
     186             :     const bool orthorhombic, const int border_mask, const enum grid_func func,
     187             :     const int la_max, const int la_min, const int lb_max, const int lb_min,
     188             :     const double zeta, const double zetb, const double rscale,
     189             :     const double dh[3][3], const double dh_inv[3][3], const double ra[3],
     190             :     const double rab[3], const int npts_global[3], const int npts_local[3],
     191             :     const int shift_local[3], const int border_width[3], const double radius,
     192             :     const int o1, const int o2, const int n1, const int n2,
     193    86290985 :     const double pab[n2][n1], double *grid) {
     194             : 
     195             :   // Set this to true to write each task to a file.
     196    86290985 :   const bool DUMP_TASKS = false;
     197             : 
     198    86290985 :   double *grid_before = NULL;
     199    86290985 :   const size_t npts_local_total = npts_local[0] * npts_local[1] * npts_local[2];
     200             : 
     201    86290985 :   if (DUMP_TASKS) {
     202             :     const size_t sizeof_grid = sizeof(double) * npts_local_total;
     203             :     grid_before = malloc(sizeof_grid);
     204             :     assert(grid_before != NULL);
     205             :     memcpy(grid_before, grid, sizeof_grid);
     206             :     memset(grid, 0, sizeof_grid);
     207             :   }
     208             : 
     209    86290985 :   collocate_internal(orthorhombic, border_mask, func, la_max, la_min, lb_max,
     210             :                      lb_min, zeta, zetb, rscale, dh, dh_inv, ra, rab,
     211             :                      npts_global, npts_local, shift_local, border_width, radius,
     212             :                      o1, o2, n1, n2, pab, grid);
     213             : 
     214    86290985 :   if (DUMP_TASKS) {
     215             :     write_task_file(orthorhombic, border_mask, func, la_max, la_min, lb_max,
     216             :                     lb_min, zeta, zetb, rscale, dh, dh_inv, ra, rab,
     217             :                     npts_global, npts_local, shift_local, border_width, radius,
     218             :                     o1, o2, n1, n2, pab, grid);
     219             : 
     220             :     for (size_t i = 0; i < npts_local_total; i++) {
     221             :       grid[i] += grid_before[i];
     222             :     }
     223             :     free(grid_before);
     224             :   }
     225    86290985 : }
     226             : 
     227             : // EOF