Commit aa12f4f7 authored by Anne Gunn's avatar Anne Gunn
Browse files

Per request from Axel, Revert revert one commit too far: 

"Per suggestion from Axel, reformat code to what I believe to be"

This reverts commit b72c1c0e.
parent b72c1c0e

src/USER-BOCS/compute_pressure_bocs.cpp

+146 −157
Original line number Diff line number Diff line
@@ -16,10 +16,10 @@ 


#include "compute_pressure_bocs.h"



#include <mpi.h>

#include <cmath>

#include <cstdlib>

#include <cstring>

#include <mpi.h>

#include <cstdlib>



#include "angle.h"

#include "atom.h"
@@ -37,12 +37,14 @@ 
#include "pair.h"

#include "update.h"





using namespace LAMMPS_NS;



/* ---------------------------------------------------------------------- */



ComputePressureBocs::ComputePressureBocs(LAMMPS *lmp, int narg, char **arg) :

    Compute(lmp, narg, arg), vptr(NULL), id_temp(NULL)

  Compute(lmp, narg, arg),

  vptr(NULL), id_temp(NULL)

{

  if (narg < 4) error->all(FLERR,"Illegal compute pressure/bocs command");

  if (igroup) error->all(FLERR,"Compute pressure/bocs must use group all");
@@ -60,15 +62,15 @@ ComputePressureBocs::ComputePressureBocs(LAMMPS *lmp, int narg, char **arg) : 
  // store temperature ID used by pressure computation

  // insure it is valid for temperature computation



  if (strcmp(arg[3], "NULL") == 0)

    id_temp = NULL;

  if (strcmp(arg[3],"NULL") == 0) id_temp = NULL;

  else {

    int n = strlen(arg[3]) + 1;

    id_temp = new char[n];

    strcpy(id_temp,arg[3]);



    int icompute = modify->find_compute(id_temp);

    if (icompute < 0) error->all(FLERR, "Could not find compute pressure/bocs temperature ID");

    if (icompute < 0)

      error->all(FLERR,"Could not find compute pressure/bocs temperature ID");

    if (modify->compute[icompute]->tempflag == 0)

      error->all(FLERR,"Compute pressure/bocs temperature ID does not "

                 "compute temperature");
@@ -81,37 +83,26 @@ ComputePressureBocs::ComputePressureBocs(LAMMPS *lmp, int narg, char **arg) : 
    pairflag = 1;

    bondflag = angleflag = dihedralflag = improperflag = 1;

    kspaceflag = fixflag = 1;

  }

  else {

  } else {

    keflag = 0;

    pairflag = 0;

    bondflag = angleflag = dihedralflag = improperflag = 0;

    kspaceflag = fixflag = 0;

    int iarg = 4;

    while (iarg < narg) {

      if (strcmp(arg[iarg], "ke") == 0)

        keflag = 1;

      else if (strcmp(arg[iarg], "pair") == 0)

        pairflag = 1;

      else if (strcmp(arg[iarg], "bond") == 0)

        bondflag = 1;

      else if (strcmp(arg[iarg], "angle") == 0)

        angleflag = 1;

      else if (strcmp(arg[iarg], "dihedral") == 0)

        dihedralflag = 1;

      else if (strcmp(arg[iarg], "improper") == 0)

        improperflag = 1;

      else if (strcmp(arg[iarg], "kspace") == 0)

        kspaceflag = 1;

      else if (strcmp(arg[iarg], "fix") == 0)

        fixflag = 1;

      if (strcmp(arg[iarg],"ke") == 0) keflag = 1;

      else if (strcmp(arg[iarg],"pair") == 0) pairflag = 1;

      else if (strcmp(arg[iarg],"bond") == 0) bondflag = 1;

      else if (strcmp(arg[iarg],"angle") == 0) angleflag = 1;

      else if (strcmp(arg[iarg],"dihedral") == 0) dihedralflag = 1;

      else if (strcmp(arg[iarg],"improper") == 0) improperflag = 1;

      else if (strcmp(arg[iarg],"kspace") == 0) kspaceflag = 1;

      else if (strcmp(arg[iarg],"fix") == 0) fixflag = 1;

      else if (strcmp(arg[iarg],"virial") == 0) {

        pairflag = 1;

        bondflag = angleflag = dihedralflag = improperflag = 1;

        kspaceflag = fixflag = 1;

      }

      else

        error->all(FLERR, "Illegal compute pressure/bocs command");

      } else error->all(FLERR,"Illegal compute pressure/bocs command");

      iarg++;

    }

  }
@@ -153,7 +144,8 @@ void ComputePressureBocs::init() 


  if (keflag) {

    int icompute = modify->find_compute(id_temp);

    if (icompute < 0) error->all(FLERR, "Could not find compute pressure/bocs temperature ID");

    if (icompute < 0)

      error->all(FLERR,"Could not find compute pressure/bocs temperature ID");

    temperature = modify->compute[icompute];

  }
@@ -179,19 +171,20 @@ void ComputePressureBocs::init() 
    if (pairflag && force->pair) vptr[nvirial++] = force->pair->virial;

    if (bondflag && force->bond) vptr[nvirial++] = force->bond->virial;

    if (angleflag && force->angle) vptr[nvirial++] = force->angle->virial;

    if (dihedralflag && force->dihedral) vptr[nvirial++] = force->dihedral->virial;

    if (improperflag && force->improper) vptr[nvirial++] = force->improper->virial;

    if (dihedralflag && force->dihedral)

      vptr[nvirial++] = force->dihedral->virial;

    if (improperflag && force->improper)

      vptr[nvirial++] = force->improper->virial;

    if (fixflag)

      for (int i = 0; i < modify->nfix; i++)

        if (modify->fix[i]->virial_flag) vptr[nvirial++] = modify->fix[i]->virial;

        if (modify->fix[i]->virial_flag)

          vptr[nvirial++] = modify->fix[i]->virial;

  }



  // flag Kspace contribution separately, since not summed across procs



  if (kspaceflag && force->kspace)

    kspace_virial = force->kspace->virial;

  else

    kspace_virial = NULL;

  if (kspaceflag && force->kspace) kspace_virial = force->kspace->virial;

  else kspace_virial = NULL;

}



/* Extra functions added for BOCS */
@@ -199,12 +192,14 @@ void ComputePressureBocs::init() 
/* ----------------------------------------------------------------------

   Compute the pressure correction for the analytical basis set

------------------------------------------------------------------------- */

double ComputePressureBocs::get_cg_p_corr(int N_basis, double *phi_coeff, int N_mol, double vavg,

                                          double vCG)

double ComputePressureBocs::get_cg_p_corr(int N_basis, double *phi_coeff,

                                      int N_mol, double vavg, double vCG)

{

  double correction = 0.0;

  for (int i = 1; i <= N_basis; ++i) {

    correction -= phi_coeff[i - 1] * (N_mol * i / vavg) * pow((1 / vavg) * (vCG - vavg), i - 1);

  for (int i = 1; i <= N_basis; ++i)

  {

    correction -= phi_coeff[i-1] * ( N_mol * i / vavg ) *

                                   pow( ( 1 / vavg ) * ( vCG - vavg ),i-1);

  }

  return correction;

}
@@ -217,18 +212,16 @@ double ComputePressureBocs::find_index(double *grid, double value) 
  int i;

  double spacing = fabs(grid[1]-grid[0]);

  int gridsize = spline_length;

  for (i = 0; i < (gridsize - 1); ++i) {

    if (value >= grid[i] && value <= grid[i + 1]) {

      return i;

    }

  for (i = 0; i < (gridsize-1); ++i)

  {

    if (value >= grid[i] && value <= grid[i+1]) { return i; }

  }



  if (value >= grid[i] && value <= (grid[i] + spacing)) {

    return i;

  }

  if (value >= grid[i] && value <= (grid[i] + spacing)) { return i; }



  error->all(FLERR, fmt::format("find_index could not find value in grid for value: {}", value));

  for (int i = 0; i < gridsize; ++i) {

  for (int i = 0; i < gridsize; ++i)

  {

    fprintf(stderr, "grid %d: %f\n",i,grid[i]);

  }
@@ -239,20 +232,24 @@ double ComputePressureBocs::find_index(double *grid, double value) 
   Compute the pressure correction for a spline basis set

------------------------------------------------------------------------- */



double ComputePressureBocs::get_cg_p_corr(double **grid, int basis_type, double vCG)

double ComputePressureBocs::get_cg_p_corr(double ** grid, int basis_type,

                                                               double vCG)

{

  int i = find_index(grid[0],vCG);

  double correction, deltax = vCG - grid[0][i];



  if (basis_type == BASIS_LINEAR_SPLINE) {

    correction =

        grid[1][i] + (deltax) * (grid[1][i + 1] - grid[1][i]) / (grid[0][i + 1] - grid[0][i]);

  if (basis_type == BASIS_LINEAR_SPLINE)

  {

    correction = grid[1][i] + (deltax) *

          ( grid[1][i+1] - grid[1][i] ) / ( grid[0][i+1] - grid[0][i] );

  }

  else if (basis_type == BASIS_CUBIC_SPLINE) {

    correction = grid[1][i] + (grid[2][i] * deltax) + (grid[3][i] * pow(deltax, 2)) +

                 (grid[4][i] * pow(deltax, 3));

  else if (basis_type == BASIS_CUBIC_SPLINE)

  {

    correction = grid[1][i] + (grid[2][i] * deltax) +

            (grid[3][i] * pow(deltax,2)) + (grid[4][i] * pow(deltax,3));

  }

  else {

  else

  {

    error->all(FLERR,"bad spline type passed to get_cg_p_corr()\n");

  }

  return correction;
@@ -262,13 +259,12 @@ double ComputePressureBocs::get_cg_p_corr(double **grid, int basis_type, double 
   send cg info from fix_bocs to compute_pressure_bocs for the analytical

   basis set

------------------------------------------------------------------------- */

void ComputePressureBocs::send_cg_info(int basis_type, int sent_N_basis, double *sent_phi_coeff,

                                       int sent_N_mol, double sent_vavg)

void ComputePressureBocs::send_cg_info(int basis_type, int sent_N_basis,

                double *sent_phi_coeff, int sent_N_mol, double sent_vavg)

{

  if (basis_type == BASIS_ANALYTIC) { p_basis_type = BASIS_ANALYTIC; }

  else

  {

  if (basis_type == BASIS_ANALYTIC) {

    p_basis_type = BASIS_ANALYTIC;

  }

  else {

    error->all(FLERR,"Incorrect basis type passed to ComputePressureBocs\n");

  }
@@ -277,9 +273,7 @@ void ComputePressureBocs::send_cg_info(int basis_type, int sent_N_basis, double 
  N_basis = sent_N_basis;

  if (phi_coeff) free(phi_coeff);

  phi_coeff = ((double *) calloc(N_basis, sizeof(double)) );

  for (int i = 0; i < N_basis; i++) {

    phi_coeff[i] = sent_phi_coeff[i];

  }

  for (int i=0; i<N_basis; i++) { phi_coeff[i] = sent_phi_coeff[i]; }



  N_mol = sent_N_mol;

  vavg = sent_vavg;
@@ -289,15 +283,13 @@ void ComputePressureBocs::send_cg_info(int basis_type, int sent_N_basis, double 
   send cg info from fix_bocs to compute_pressure_bocs for a spline basis

   set

------------------------------------------------------------------------- */

void ComputePressureBocs::send_cg_info(int basis_type, double **in_splines, int gridsize)

void ComputePressureBocs::send_cg_info(int basis_type,

                                         double ** in_splines, int gridsize)

{

  if (basis_type == BASIS_LINEAR_SPLINE) { p_basis_type = BASIS_LINEAR_SPLINE; }

  else if (basis_type == BASIS_CUBIC_SPLINE) { p_basis_type = BASIS_CUBIC_SPLINE; }

  else

  {

  if (basis_type == BASIS_LINEAR_SPLINE) {

    p_basis_type = BASIS_LINEAR_SPLINE;

  }

  else if (basis_type == BASIS_CUBIC_SPLINE) {

    p_basis_type = BASIS_CUBIC_SPLINE;

  }

  else {

    error->all(FLERR,"Incorrect basis type passed to ComputePressureBocs\n");

  }

  splines = in_splines;
@@ -323,8 +315,7 @@ double ComputePressureBocs::compute_scalar() 
  if (keflag) {

    if (temperature->invoked_scalar != update->ntimestep)

      t = temperature->compute_scalar();

    else

      t = temperature->scalar;

    else t = temperature->scalar;

  }



  if (dimension == 3) {
@@ -332,30 +323,34 @@ double ComputePressureBocs::compute_scalar() 
    volume = (domain->xprd * domain->yprd * domain->zprd);



    /* MRD NJD if block */

    if (p_basis_type == BASIS_ANALYTIC) {

    if ( p_basis_type == BASIS_ANALYTIC )

    {

      correction = get_cg_p_corr(N_basis,phi_coeff,N_mol,vavg,volume);

    }

    else if (p_basis_type == BASIS_LINEAR_SPLINE || p_basis_type == BASIS_CUBIC_SPLINE) {

    else if ( p_basis_type == BASIS_LINEAR_SPLINE || p_basis_type == BASIS_CUBIC_SPLINE )

    {

      correction = get_cg_p_corr(splines, p_basis_type, volume);

    }



    virial_compute(3,3);

    if (keflag)

      scalar = (temperature->dof * boltz * t + virial[0] + virial[1] + virial[2]) / 3.0 *

                   inv_volume * nktv2p +

               (correction);

      scalar = (temperature->dof * boltz * t +

                virial[0] + virial[1] + virial[2]) / 3.0 *

                inv_volume * nktv2p + (correction);

    else

      scalar = (virial[0] + virial[1] + virial[2]) / 3.0 * inv_volume * nktv2p + (correction);

  }

  else {

    if (p_match_flag) {

      scalar = (virial[0] + virial[1] + virial[2]) / 3.0 *

               inv_volume * nktv2p + (correction);

  } else {

    if (p_match_flag)

    {

      error->all(FLERR,"Pressure matching not implemented in 2-d.\n");

      exit(1);

    } // The rest of this can probably be deleted.

    inv_volume = 1.0 / (domain->xprd * domain->yprd);

    virial_compute(2,2);

    if (keflag)

      scalar = (temperature->dof * boltz * t + virial[0] + virial[1]) / 2.0 * inv_volume * nktv2p;

      scalar = (temperature->dof * boltz * t +

                virial[0] + virial[1]) / 2.0 * inv_volume * nktv2p;

    else

      scalar = (virial[0] + virial[1]) / 2.0 * inv_volume * nktv2p;

  }
@@ -382,7 +377,8 @@ void ComputePressureBocs::compute_vector() 


  double *ke_tensor;

  if (keflag) {

    if (temperature->invoked_vector != update->ntimestep) temperature->compute_vector();

    if (temperature->invoked_vector != update->ntimestep)

      temperature->compute_vector();

    ke_tensor = temperature->vector;

  }
@@ -392,12 +388,10 @@ void ComputePressureBocs::compute_vector() 
    if (keflag) {

      for (int i = 0; i < 6; i++)

        vector[i] = (ke_tensor[i] + virial[i]) * inv_volume * nktv2p;

    }

    else

    } else

      for (int i = 0; i < 6; i++)

        vector[i] = virial[i] * inv_volume * nktv2p;

  }

  else {

  } else {

    inv_volume = 1.0 / (domain->xprd * domain->yprd);

    virial_compute(4,2);

    if (keflag) {
@@ -405,8 +399,7 @@ void ComputePressureBocs::compute_vector() 
      vector[1] = (ke_tensor[1] + virial[1]) * inv_volume * nktv2p;

      vector[3] = (ke_tensor[3] + virial[3]) * inv_volume * nktv2p;

      vector[2] = vector[4] = vector[5] = 0.0;

    }

    else {

    } else {

      vector[0] = virial[0] * inv_volume * nktv2p;

      vector[1] = virial[1] * inv_volume * nktv2p;

      vector[3] = virial[3] * inv_volume * nktv2p;
@@ -422,15 +415,13 @@ void ComputePressureBocs::virial_compute(int n, int ndiag) 
  int i,j;

  double v[6],*vcomponent;



  for (i = 0; i < n; i++)

    v[i] = 0.0;

  for (i = 0; i < n; i++) v[i] = 0.0;



  // sum contributions to virial from forces and fixes



  for (j = 0; j < nvirial; j++) {

    vcomponent = vptr[j];

    for (i = 0; i < n; i++)

      v[i] += vcomponent[i];

    for (i = 0; i < n; i++) v[i] += vcomponent[i];

  }



  // sum virial across procs
@@ -440,14 +431,12 @@ void ComputePressureBocs::virial_compute(int n, int ndiag) 
  // KSpace virial contribution is already summed across procs



  if (kspace_virial)

    for (i = 0; i < n; i++)

      virial[i] += kspace_virial[i];

    for (i = 0; i < n; i++) virial[i] += kspace_virial[i];



  // LJ long-range tail correction, only if pair contributions are included



  if (force->pair && pairflag && force->pair->tail_flag)

    for (i = 0; i < ndiag; i++)

      virial[i] += force->pair->ptail * inv_volume;

    for (i = 0; i < ndiag; i++) virial[i] += force->pair->ptail * inv_volume;

}



/* ---------------------------------------------------------------------- */

src/USER-BOCS/compute_pressure_bocs.h

+45 −44
Original line number Diff line number Diff line
@@ -20,6 +20,7 @@ ComputeStyle(PRESSURE/BOCS, ComputePressureBocs) 


#else





#ifndef LMP_COMPUTE_PRESSURE_BOCS_H

#define LMP_COMPUTE_PRESSURE_BOCS_H
@@ -76,7 +77,7 @@ class ComputePressureBocs : public Compute { 
  void virial_compute(int, int);

};



} // namespace LAMMPS_NS

}



#endif

#endif

src/USER-BOCS/fix_bocs.cpp

+743 −746

File changed.

Preview size limit exceeded, changes collapsed.

CRA Git | Maintained and supported by SUSTech CRA and CCSE