fcs_p3m_p.h

Go to the documentation of this file.

/*
  Copyright (C) 2011,2012 Rene Halver, Olaf Lenz

  This file is part of ScaFaCoS.

  ScaFaCoS is free software: you can redistribute it and/or modify
  it under the terms of the GNU Lesser Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  ScaFaCoS is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU Lesser Public License for more details.

  You should have received a copy of the GNU Lesser Public License
  along with this program.  If not, see <http://www.gnu.org/licenses/>. 
*/

#ifndef _FCS_P3M_P_H
#define _FCS_P3M_P_H

#include "fcs_definitions.h"
#include "fcs_result_p.h"
#include "fcs_interface_p.h"
#include <math.h>

#ifdef __cplusplus
extern "C" {
#endif

/* TODO: test whether to use it or not */
#define FCS_P3M_USE_ERFC_APPROXIMATION 1
  /* #define FCS_P3M_USE_ERFC_APPROXIMATION 0 */

/* This struct is defined open so that an MD implementation can use
   the parameters to perform the near field computation in its own
   code. */
typedef struct {fcs_float alpha; fcs_float potentialOffset;} fcs_p3m_near_parameters_t;

FCSResult 
fcs_p3m_get_near_parameters(FCS handle,
                fcs_p3m_near_parameters_t *near_params);

/* This function is defined inline for maximal performance! */
static inline fcs_float 
fcs_p3m_compute_near_potential(fcs_p3m_near_parameters_t params, fcs_float dist) {
  const fcs_float alpha = params.alpha;
  const fcs_float potentialOffset = params.potentialOffset;  
  fcs_float adist = alpha * dist;

#if FCS_P3M_USE_ERFC_APPROXIMATION

  /* approximate \f$ \exp(d^2) \mathrm{erfc}(d)\f$ by applying a formula from:
     Abramowitz/Stegun: Handbook of Mathematical Functions, Dover
     (9. ed.), chapter 7 */
  fcs_float t = 1.0 / (1.0 + 0.3275911 * adist);
  fcs_float erfc_part_ri = exp(-adist*adist) * 
    (t * (0.254829592 + 
      t * (-0.284496736 + 
           t * (1.421413741 + 
            t * (-1.453152027 + 
             t * 1.061405429))))) 
    / dist;
  
#else

  /* use erf instead of erfc to fix ICC performance problems */
  fcs_float erfc_part_ri = (1.0 - erf(adist)) / dist; 

#endif

  return erfc_part_ri-potentialOffset;

}

static inline fcs_float 
fcs_p3m_compute_near_field(fcs_p3m_near_parameters_t params, fcs_float dist) {
  const fcs_float alpha = params.alpha;
  fcs_float adist = alpha * dist;

#if FCS_P3M_USE_ERFC_APPROXIMATION

  /* approximate \f$ \exp(d^2) \mathrm{erfc}(d)\f$ by applying a formula from:
     Abramowitz/Stegun: Handbook of Mathematical Functions, Dover
     (9. ed.), chapter 7 */
  fcs_float t = 1.0 / (1.0 + 0.3275911 * adist);
  fcs_float erfc_part_ri = exp(-adist*adist) * 
    (t * (0.254829592 + 
      t * (-0.284496736 + 
           t * (1.421413741 + 
            t * (-1.453152027 + 
             t * 1.061405429))))) 
    / dist;
  /* 1/sqrt(PI) = 0.56418958354775627928034964498 */
  return -(erfc_part_ri + 2.0*alpha*0.56418958354775627928034964498*exp(-adist*adist)) 
    / dist;
  
#else

  fcs_float erfc_part_ri = (1.0 - erf(adist)) / dist; /* use erf instead of erfc to fix ICC performance problems */
  /* 1/sqrt(PI) = 0.56418958354775627928034964498 */
  return -(erfc_part_ri + 2.0*alpha*0.56418958354775627928034964498*exp(-adist*adist))
    / dist;

#endif

}


static inline void
fcs_p3m_compute_near(fcs_p3m_near_parameters_t params, fcs_float dist, 
             fcs_float *potential, fcs_float *field) {
  const fcs_float alpha = params.alpha;
  const fcs_float potentialOffset = params.potentialOffset;
  fcs_float adist = alpha * dist;

#if FCS_P3M_USE_ERFC_APPROXIMATION

  /* approximate \f$ \exp(d^2) \mathrm{erfc}(d)\f$ by applying a formula from:
     Abramowitz/Stegun: Handbook of Mathematical Functions, Dover
     (9. ed.), chapter 7 */
  fcs_float t = 1.0 / (1.0 + 0.3275911 * adist);
  fcs_float erfc_part_ri = exp(-adist*adist) * 
    (t * (0.254829592 + 
      t * (-0.284496736 + 
           t * (1.421413741 + 
            t * (-1.453152027 + 
             t * 1.061405429))))) 
    / dist;

  *potential = erfc_part_ri-potentialOffset;
  *field = -(erfc_part_ri + 2.0*alpha*0.56418958354775627928034964498*exp(-adist*adist))
    / dist;
  
#else

  fcs_float erfc_part_ri = (1.0 - erf(adist)) / dist; /* use erf instead of erfc to fix ICC performance problems */
  *potential = erfc_part_ri-potentialOffset;
  *field = -(erfc_part_ri + 2.0*alpha*0.56418958354775627928034964498*exp(-adist*adist)) 
    / dist;

#endif

}

FCSResult fcs_p3m_distribute_parameters();

FCSResult fcs_p3m_set_r_cut(FCS handle, fcs_float r_cut);
FCSResult fcs_p3m_set_r_cut_tune(FCS handle);
FCSResult fcs_p3m_get_r_cut(FCS handle, fcs_float *r_cut);

FCSResult fcs_p3m_set_alpha(FCS handle, fcs_float alpha);
FCSResult fcs_p3m_set_alpha_tune(FCS handle);
FCSResult fcs_p3m_get_alpha(FCS handle, fcs_float *alpha);

FCSResult fcs_p3m_set_grid(FCS handle, fcs_int grid);
FCSResult fcs_p3m_set_grid_tune(FCS handle);
FCSResult fcs_p3m_get_grid(FCS handle, fcs_int *grid);

FCSResult fcs_p3m_set_cao(FCS handle, fcs_int cao);
FCSResult fcs_p3m_set_cao_tune(FCS handle);
FCSResult fcs_p3m_get_cao(FCS handle, fcs_int *cao);

FCSResult fcs_p3m_require_total_energy(FCS handle, fcs_int total_energy);
FCSResult fcs_p3m_get_total_energy(FCS handle, fcs_float *total_energy);

FCSResult fcs_p3m_set_potential_shift(FCS handle, fcs_int flag);
fcs_int fcs_p3m_get_potential_shift(FCS handle);

FCSResult fcs_p3m_set_tolerance_field(FCS handle, fcs_float tolerance_field);
/* FORTRAN wrapper */
void fcs_p3m_set_tolerance_field_f(void *handle, fcs_float tolerance_field, fcs_int *return_value);
FCSResult fcs_p3m_set_tolerance_field_tune(FCS handle);
FCSResult fcs_p3m_get_tolerance_field(FCS handle, fcs_float *tolerance_field);



  /*
FCSResult fcs_p3m_set_variant(FCS handle, fcs_int variant);
FCSResult fcs_p3m_get_variant(FCS handle, fcs_int *variant);
  */

#ifdef __cplusplus
}
#endif

#endif