Commit d3439bc8 authored by eric pellegrini's avatar eric pellegrini
Browse files

Solved the problem with erfc

Removed unused preprocessor generated files
parent d727f294
......@@ -17,7 +17,7 @@
* Include erfc code, unless user wants to use an erfc from libm.
*/
#ifndef LIBM_HAS_ERFC
#include "ndtr.h"
#include "MMTK/ndtr.h"
#endif
#ifndef M_PI
#define M_PI 3.14159265358979323846264338327950288
......
for (k = 0; k < nblist->npairs; k++, pair++) {
if (cutoff_sq == 0. || pair->r_sq <= cutoff_sq) {
int i = pair->i1;
int j = pair->i2;
double r_sq = pair->r_sq;
double r = sqrt(r_sq);
double f1 = erfc(beta*r);
double f2 = erfc_cutoff;
vector3 rij;
if (pair->one_four)
f1 += one_four_factor-1.;
#ifdef NBLIST_WITH_VECTORS
vector_copy(rij, pair->d);
#else
distance_vector_function(rij, x[j], x[i], universe_data);
#endif
# include "ewald2.i"
}
}
for (k = 0; k < n_ex; k += 2) {
int i = excluded[k];
int j = excluded[k+1];
double r_sq, r, f1, f2;
vector3 rij;
distance_vector_function(rij, x[j], x[i], universe_data);
r_sq = vector_length_sq(rij);
r = sqrt(r_sq);
f1 = erfc(beta*r)-1.;
f2 = erfc_cutoff;
# include "ewald2.i"
}
{
double qiqj = charge[i]*charge[j]*electrostatic_energy_factor;
double ef = 2./sqrt(M_PI);
self->energy += qiqj*(f1/r-f2*inv_cutoff);
if (gradients != NULL || force_constants != NULL) {
double deriv = -qiqj*(f1/r_sq-f2*inv_cutoff
+ef*beta*(exp(-beta*beta*r_sq)/r
-inv_cutoff*exp(-beta*beta*cutoff_sq)))/r;
if (gradients != NULL) {
vector3 grad;
grad[0] = deriv*rij[0];
grad[1] = deriv*rij[1];
grad[2] = deriv*rij[2];
if (gradient_fn != NULL) {
(*gradient_fn)(self, gradients, i, grad);
vector_changesign(grad);
(*gradient_fn)(self, gradients, j, grad);
}
else {
vector3 *f = (vector3 *)((PyArrayObject *)gradients)->data;
f[i][0] += grad[0];
f[i][1] += grad[1];
f[i][2] += grad[2];
f[j][0] -= grad[0];
f[j][1] -= grad[1];
f[j][2] -= grad[2];
}
}
if (force_constants != NULL) {
double deriv2 = 2.*qiqj*(f1/(r*r_sq) - f2*cube(inv_cutoff)
+ ef*beta*exp(-beta*beta*r_sq)
* (beta*beta+1./r_sq));
if (inv_cutoff > 0.)
deriv2 -= 2.*qiqj*ef*beta*exp(-beta*beta*cutoff_sq)
* (beta*beta+sqr(inv_cutoff));
add_pair_fc(self, force_constants, fc_fn, i, j,
rij, r_sq, deriv, deriv2);
}
}
}
......@@ -145,8 +145,9 @@ Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
#include "polevl.h"
#include "ndtr.h"
#include "MMTK/mconf.h"
#include "MMTK/ndtr.h"
#include "MMTK/polevl.h"
#ifdef UNK
static double P[] = {
......
......@@ -8,7 +8,7 @@
#define PY_ARRAY_UNIQUE_SYMBOL PyArray_MMTKFF_API
#ifndef LIBM_HAS_ERFC
#include "ndtr.h"
#include "MMTK/ndtr.h"
#endif
#include "MMTK/forcefield.h"
......
for (k = 0; k < n_ex; k += 2) {
int i = excluded[k];
int j = excluded[k+1];
double qiqj = charge[i]*charge[j]*electrostatic_energy_factor;
pair_es_energy(-qiqj, distance_vector_function);
}
for (k = 0; k < n_14; k += 2) {
int i = one_four[k];
int j = one_four[k+1];
double qiqj = charge[i]*charge[j]*electrostatic_energy_factor;
pair_es_energy(qiqj*(es_one_four_factor-1.), distance_vector_function);
}
......@@ -50,7 +50,7 @@ Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
#include "polevl.h"
#include "MMTK/polevl.h"
double polevl( double x, double coef[], int N )
{
......
......@@ -337,7 +337,7 @@ simulations.
Extension('MMTK_forcefield',
['Src/MMTK_forcefield.c',
'Src/bonded.c', 'Src/nonbonded.c',
'Src/ewald.c', 'Src/sparsefc.c'],
'Src/ewald.c', 'Src/sparsefc.c','Src/ndtr.c','Src/polevl.c'],
extra_compile_args = compile_args + high_opt,
include_dirs=include_dirs + ['Src'],
define_macros = [('SERIAL', None),
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment