From e9fce479e70e6b4e1e16e34705cf0870030e17a3 Mon Sep 17 00:00:00 2001 From: Owen Hickey Date: Thu, 19 Jan 2012 16:01:42 +0100 Subject: [PATCH] neutral particles out of ES short range --- src/forces.h | 37 +++--- src/iccp3m.c | 36 ++++++ src/iccp3m.h | 74 ++---------- src/initialize.c | 25 +--- src/interaction_data.c | 1 - src/mmm1d.c | 5 +- src/mmm1d.h | 2 +- src/mmm2d.c | 258 ++++++++++++++++++++--------------------- 8 files changed, 201 insertions(+), 237 deletions(-) diff --git a/src/forces.h b/src/forces.h index 93178f2cef1..48a4dc1ce18 100644 --- a/src/forces.h +++ b/src/forces.h @@ -299,43 +299,50 @@ MDINLINE void add_non_bonded_pair_force(Particle *p1, Particle *p2, #ifdef ELECTROSTATICS /* real space coulomb */ + double q1q2 = p1->p.q*p2->p.q; switch (coulomb.method) { #ifdef P3M case COULOMB_ELC_P3M: { - p3m_add_pair_force(p1->p.q*p2->p.q,d,dist2,dist,force); + if (q1q2) { + p3m_add_pair_force(q1q2,d,dist2,dist,force); - // forces from the virtual charges - // they go directly onto the particles, since they are not pairwise forces - if (elc_params.dielectric_contrast_on) - ELC_P3M_dielectric_layers_force_contribution(p1, p2, p1->f.f, p2->f.f); + // forces from the virtual charges + // they go directly onto the particles, since they are not pairwise forces + if (elc_params.dielectric_contrast_on) + ELC_P3M_dielectric_layers_force_contribution(p1, p2, p1->f.f, p2->f.f); + } break; } case COULOMB_P3M: { #ifdef NPT - double eng = p3m_add_pair_force(p1->p.q*p2->p.q,d,dist2,dist,force); - if(integ_switch == INTEG_METHOD_NPT_ISO) - nptiso.p_vir[0] += eng; + if (q1q2) { + double eng = p3m_add_pair_force(q1q2,d,dist2,dist,force); + if(integ_switch == INTEG_METHOD_NPT_ISO) + nptiso.p_vir[0] += eng; + } #else - p3m_add_pair_force(p1->p.q*p2->p.q,d,dist2,dist,force); + if (q1q2) p3m_add_pair_force(q1q2,d,dist2,dist,force); #endif break; } #endif case COULOMB_EWALD: { #ifdef NPT - double eng = add_ewald_coulomb_pair_force(p1,p2,d,dist2,dist,force); - if(integ_switch == INTEG_METHOD_NPT_ISO) - nptiso.p_vir[0] += eng; + if (q1q2) { + double eng = add_ewald_coulomb_pair_force(p1,p2,d,dist2,dist,force); + if(integ_switch == INTEG_METHOD_NPT_ISO) + nptiso.p_vir[0] += eng; + } #else - add_ewald_coulomb_pair_force(p1,p2,d,dist2,dist,force); + if (q1q2) add_ewald_coulomb_pair_force(p1,p2,d,dist2,dist,force); #endif break; } case COULOMB_MMM1D: - add_mmm1d_coulomb_pair_force(p1,p2,d,dist2,dist,force); + if (q1q2) add_mmm1d_coulomb_pair_force(q1q2,d,dist2,dist,force); break; case COULOMB_MMM2D: - add_mmm2d_coulomb_pair_force(p1->p.q*p2->p.q,d,dist2,dist,force); + if (q1q2) add_mmm2d_coulomb_pair_force(q1q2,d,dist2,dist,force); break; case COULOMB_NONE: break; diff --git a/src/iccp3m.c b/src/iccp3m.c index 4cc81ae3f70..979ef12493d 100644 --- a/src/iccp3m.c +++ b/src/iccp3m.c @@ -148,6 +148,8 @@ int iccp3m_iteration() { char* errtxt; double globalmax; + iccp3m_sanity_check(); + l_b = coulomb.bjerrum; if((iccp3m_cfg.eout <= 0)) { errtxt = runtime_error(128); @@ -719,5 +721,39 @@ void iccp3m_store_forces() { } } +int iccp3m_sanity_check() +{ + switch (coulomb.method) { + case COULOMB_ELC_P3M: { + if (elc_params.dielectric_contrast_on) { + char *errtxt = runtime_error(128); + ERROR_SPRINTF(errtxt, "ICCP3M conflicts with ELC dielectric constrast"); + return 1; + } + break; + } + case COULOMB_DH: { + char *errtxt = runtime_error(128); + ERROR_SPRINTF(errtxt, "ICCP3M does not work with Debye-Hueckel iccp3m.h"); + return 1; + } + case COULOMB_RF: { + char *errtxt = runtime_error(128); + ERROR_SPRINTF(errtxt, "ICCP3M does not work with COULOMB_RF iccp3m.h"); + return 1; + } + } + +#ifdef NPT + if(integ_switch == INTEG_METHOD_NPT_ISO) { + char *errtxt = runtime_error(128); + ERROR_SPRINTF(errtxt, "ICCP3M does not work in the NPT ensemble"); + return 1; + } +#endif + + return 0; +} + #endif diff --git a/src/iccp3m.h b/src/iccp3m.h index f545742846c..5b3501c173a 100644 --- a/src/iccp3m.h +++ b/src/iccp3m.h @@ -138,6 +138,9 @@ int iccp3m_iteration(); */ void iccp3m_init(void); +/** check sanity of parameters for use with ICCP3M + */ +int iccp3m_sanity_check(); /** The short range part of the electrostatic interation between two particles. * The appropriate function from the underlying electrostatic method is called. */ @@ -147,86 +150,33 @@ MDINLINE void add_non_bonded_pair_force_iccp3m(Particle *p1, Particle *p2, /* IA_parameters *ia_params = get_ia_param(p1->p.type,p2->p.type);*/ double force[3] = { 0, 0, 0 }; int j; - char* errtxt; FORCE_TRACE(fprintf(stderr, "%d: interaction %d<->%d dist %f\n", this_node, p1->p.identity, p2->p.identity, dist)); - /***********************************************/ - /* short range electrostatics */ - /***********************************************/ - - if (coulomb.method == COULOMB_DH) { - errtxt = runtime_error(128); - ERROR_SPRINTF(errtxt, "ICCP3M does not work with Debye-Hueckel iccp3m.h"); - } - if (coulomb.method == COULOMB_RF) { - errtxt = runtime_error(128); - ERROR_SPRINTF(errtxt, "ICCP3M does not work with COULOMB_RF iccp3m.h"); - } - - /*********************************************************************/ - /* everything before this contributes to the virial pressure in NpT, */ - /* but nothing afterwards */ - /*********************************************************************/ -#ifdef NPT - for (j = 0; j < 3; j++) - if(integ_switch == INTEG_METHOD_NPT_ISO) - if (coulomb.method == COULOMB_RF) { - errtxt = runtime_error(128); - ERROR_SPRINTF(errtxt, "ICCP3M does not work in the NPT ensemble"); - } -#endif - /***********************************************/ /* long range electrostatics */ /***********************************************/ /* real space coulomb */ + double q1q2 = p1->p.q*p2->p.q; switch (coulomb.method) { - #ifdef P3M - case COULOMB_ELC_P3M: { - p3m_add_pair_force(p1->p.q*p2->p.q,d,dist2,dist,force); - if (elc_params.dielectric_contrast_on) { - errtxt = runtime_error(128); - ERROR_SPRINTF(errtxt, "{ICCP3M conflicts with ELC dielectric constrast"); - } + case COULOMB_ELC_P3M: + if (q1q2) p3m_add_pair_force(q1q2,d,dist2,dist,force); break; - - } - - case COULOMB_P3M: { - -#ifdef NPT /* ICCP3M does not work in NPT ensemble */ - if(integ_switch == INTEG_METHOD_NPT_ISO){ - errtxt = runtime_error(128); - ERROR_SPRINTF(errtxt, "{ICCP3M cannot be used with pressure coupling} "); - } -#endif - -#ifdef DIPOLES /* ICCP3M still does not work with dipoles, so abort if compiled in */ - errtxt = runtime_error(128); - ERROR_SPRINTF(errtxt, "{ICCP3M and dipoles not implemented yet} "); - -#else /* If it is just electrostatic P3M */ - p3m_add_pair_force(p2->p.q* p1->p.q,d,dist2,dist,force); -#endif /* DIPOLES */ - + case COULOMB_P3M: + if (q1q2) p3m_add_pair_force(q1q2,d,dist2,dist,force); break; - - } #endif /* P3M */ - case COULOMB_MMM1D: - add_mmm1d_coulomb_pair_force(p1,p2,d,dist2,dist,force); + if (q1q2) add_mmm1d_coulomb_pair_force(q1q2,d,dist2,dist,force); break; case COULOMB_MMM2D: - add_mmm2d_coulomb_pair_force(p1->p.q*p2->p.q,d,dist2,dist,force); + if (q1q2) add_mmm2d_coulomb_pair_force(q1q2,d,dist2,dist,force); break; case COULOMB_NONE: break; } - /***********************************************/ /* add total nonbonded forces to particle */ @@ -234,8 +184,8 @@ MDINLINE void add_non_bonded_pair_force_iccp3m(Particle *p1, Particle *p2, for (j = 0; j < 3; j++) { p1->f.f[j] += force[j]; p2->f.f[j] -= force[j]; - } - /***********************************************/ + } + /***********************************************/ } #endif /* ELECTROSTATICS */ diff --git a/src/initialize.c b/src/initialize.c index 42d654f9909..7afe31607dc 100644 --- a/src/initialize.c +++ b/src/initialize.c @@ -28,27 +28,19 @@ #include "utils.h" #include "initialize.h" #include "global.h" -#include "tcl/global_tcl.h" #include "particle_data.h" -#include "tcl/particle_data_tcl.h" #include "interaction_data.h" -#include "tcl/interaction_data_tcl.h" -#include "tcl/binary_file_tcl.h" #include "integrate.h" #include "statistics.h" #include "energy.h" #include "pressure.h" #include "polymer.h" #include "imd.h" -#include "tcl/imd_tcl.h" #include "random.h" -#include "tcl/random_tcl.h" #include "communication.h" #include "cells.h" #include "grid.h" -#include "tcl/grid_tcl.h" #include "thermostat.h" -#include "tcl/thermostat_tcl.h" #include "rotation.h" #include "p3m.h" #include "p3m-dipolar.h" @@ -62,34 +54,19 @@ #include "debye_hueckel.h" #include "reaction_field.h" #include "forces.h" -#include "tcl/uwerr_tcl.h" #include "nsquare.h" #include "nemd.h" -#include "tcl/nemd_tcl.h" #include "domain_decomposition.h" #include "errorhandling.h" #include "rattle.h" #include "bin.h" -#include "tcl/bin_tcl.h" #include "lattice.h" #include "iccp3m.h" /* -iccp3m- */ -#include "tcl/iccp3m_tcl.h" #include "adresso.h" -#include "tcl/adresso_tcl.h" #include "metadynamics.h" -#include "tcl/integrate_tcl.h" -#include "tcl/cells_tcl.h" -#include "tcl/statistics_tcl.h" -#include "tcl/blockfile_tcl.h" -#include "tcl/iccp3m_tcl.h" -#include "tcl/polymer_tcl.h" -#include "tcl/lb-boundaries_tcl.h" #include "lb-boundaries.h" -#include "tcl/initialize_interpreter.h" -#ifdef CUDA -#include "tcl/cuda_init_tcl.h" -#endif +#include "tcl/initialize_interpreter.h" // import function from scriptsdir.c char *get_default_scriptsdir(); diff --git a/src/interaction_data.c b/src/interaction_data.c index b26d873f1f1..5046b44eebd 100644 --- a/src/interaction_data.c +++ b/src/interaction_data.c @@ -45,7 +45,6 @@ #include "hertzian.h" #include "buckingham.h" #include "soft_sphere.h" -#include "tcl/tab_tcl.h" #include "tab.h" #include "overlap.h" #include "ljcos.h" diff --git a/src/mmm1d.c b/src/mmm1d.c index f631a54d05e..c58deaf81c5 100644 --- a/src/mmm1d.c +++ b/src/mmm1d.c @@ -167,17 +167,14 @@ void MMM1D_init() MMM1D_recalcTables(); } -void add_mmm1d_coulomb_pair_force(Particle *p1, Particle *p2, double d[3], double r2, double r, double force[3]) +void add_mmm1d_coulomb_pair_force(double chpref, double d[3], double r2, double r, double force[3]) { int dim; double F[3]; - double chpref = p1->p.q*p2->p.q; double rxy2, rxy2_d, z_d; double pref; double Fx, Fy, Fz; - if (chpref == 0) - return; rxy2 = d[0]*d[0] + d[1]*d[1]; rxy2_d = rxy2*uz2; diff --git a/src/mmm1d.h b/src/mmm1d.h index 7df987593bb..07d1088cc46 100644 --- a/src/mmm1d.h +++ b/src/mmm1d.h @@ -73,7 +73,7 @@ int MMM1D_sanity_checks(); void MMM1D_init(); /// -void add_mmm1d_coulomb_pair_force(Particle *p1, Particle *p2, double d[3], double dist2, +void add_mmm1d_coulomb_pair_force(double chprf, double d[3], double dist2, double dist, double force[3]); /// diff --git a/src/mmm2d.c b/src/mmm2d.c index d9f2b98cd07..dfc3333fb14 100644 --- a/src/mmm2d.c +++ b/src/mmm2d.c @@ -1533,161 +1533,159 @@ void add_mmm2d_coulomb_pair_force(double charge_factor, } #endif - if (pref != 0.0) { - F[0] = F[1] = F[2] = 0; - - /* Bessel sum */ - { - int p, l; - double k0, k1; - double k0Sum, k1ySum, k1Sum; - double freq; - double rho_l, ypl; - double c, s; - - for (p = 1; p < besselCutoff.n; p++) { - k0Sum = 0; - k1ySum = 0; - k1Sum = 0; - - freq = C_2PI*ux*p; - - for (l = 1; l < besselCutoff.e[p-1]; l++) { - ypl = d[1] + l*box_l[1]; - rho_l = sqrt(ypl*ypl + z2); + F[0] = F[1] = F[2] = 0; + + /* Bessel sum */ + { + int p, l; + double k0, k1; + double k0Sum, k1ySum, k1Sum; + double freq; + double rho_l, ypl; + double c, s; + + for (p = 1; p < besselCutoff.n; p++) { + k0Sum = 0; + k1ySum = 0; + k1Sum = 0; + + freq = C_2PI*ux*p; + + for (l = 1; l < besselCutoff.e[p-1]; l++) { + ypl = d[1] + l*box_l[1]; + rho_l = sqrt(ypl*ypl + z2); #ifdef BESSEL_MACHINE_PREC - k0 = K0(freq*rho_l); - k1 = K1(freq*rho_l); + k0 = K0(freq*rho_l); + k1 = K1(freq*rho_l); #else - LPK01(freq*rho_l, &k0, &k1); + LPK01(freq*rho_l, &k0, &k1); #endif - k1 /= rho_l; - k0Sum += k0; - k1Sum += k1; - k1ySum += k1*ypl; + k1 /= rho_l; + k0Sum += k0; + k1Sum += k1; + k1ySum += k1*ypl; - ypl = d[1] - l*box_l[1]; - rho_l = sqrt(ypl*ypl + z2); + ypl = d[1] - l*box_l[1]; + rho_l = sqrt(ypl*ypl + z2); #ifdef BESSEL_MACHINE_PREC - k0 = K0(freq*rho_l); - k1 = K1(freq*rho_l); + k0 = K0(freq*rho_l); + k1 = K1(freq*rho_l); #else - LPK01(freq*rho_l, &k0, &k1); + LPK01(freq*rho_l, &k0, &k1); #endif - k1 /= rho_l; - k0Sum += k0; - k1Sum += k1; - k1ySum += k1*ypl; - } - - /* the ux is multiplied in to bessel, complex and psi at once, not here */ - c = 4*freq*cos(freq*d[0]); - s = 4*freq*sin(freq*d[0]); - F[0] += s*k0Sum; - F[1] += c*k1ySum; - F[2] += d[2]*c*k1Sum; + k1 /= rho_l; + k0Sum += k0; + k1Sum += k1; + k1ySum += k1*ypl; } - // fprintf(stderr, " bessel force %f %f %f\n", F[0], F[1], F[2]); + + /* the ux is multiplied in to bessel, complex and psi at once, not here */ + c = 4*freq*cos(freq*d[0]); + s = 4*freq*sin(freq*d[0]); + F[0] += s*k0Sum; + F[1] += c*k1ySum; + F[2] += d[2]*c*k1Sum; } + // fprintf(stderr, " bessel force %f %f %f\n", F[0], F[1], F[2]); + } - /* complex sum */ - { - double zeta_r, zeta_i; - double zet2_r, zet2_i; - double ztn_r, ztn_i; - double tmp_r; - int end, n; - - ztn_r = zeta_r = uy*d[2]; - ztn_i = zeta_i = uy*d[1]; - zet2_r = zeta_r*zeta_r - zeta_i*zeta_i; - zet2_i = 2*zeta_r*zeta_i; - - end = (int)ceil(COMPLEX_FAC*uy2*rho2); - if (end > COMPLEX_STEP) { - end = COMPLEX_STEP; - fprintf(stderr, "MMM2D: some particles left the assumed slab, precision might be lost\n"); - } - if (end < 0) { - char *errtxt = runtime_error(100); - ERROR_SPRINTF(errtxt, "{MMM2D: distance was negative, coordinates probably out of range } "); - end = 0; - } - end = complexCutoff[end]; + /* complex sum */ + { + double zeta_r, zeta_i; + double zet2_r, zet2_i; + double ztn_r, ztn_i; + double tmp_r; + int end, n; - for (n = 0; n < end; n++) { - F[1] -= bon.e[n]*ztn_i; - F[2] += bon.e[n]*ztn_r; + ztn_r = zeta_r = uy*d[2]; + ztn_i = zeta_i = uy*d[1]; + zet2_r = zeta_r*zeta_r - zeta_i*zeta_i; + zet2_i = 2*zeta_r*zeta_i; - tmp_r = ztn_r*zet2_r - ztn_i*zet2_i; - ztn_i = ztn_r*zet2_i + ztn_i*zet2_r; - ztn_r = tmp_r; - } - // fprintf(stderr, "complex force %f %f %f %d\n", F[0], F[1], F[2], end); + end = (int)ceil(COMPLEX_FAC*uy2*rho2); + if (end > COMPLEX_STEP) { + end = COMPLEX_STEP; + fprintf(stderr, "MMM2D: some particles left the assumed slab, precision might be lost\n"); } + if (end < 0) { + char *errtxt = runtime_error(100); + ERROR_SPRINTF(errtxt, "{MMM2D: distance was negative, coordinates probably out of range } "); + end = 0; + } + end = complexCutoff[end]; - /* psi sum */ - { - int n; - double uxx = ux*d[0]; - double uxrho2 = ux2*rho2; - double uxrho_2n, uxrho_2nm2; /* rho^{2n-2} */ - double mpe, mpo; - - /* n = 0 inflicts only Fx and pot */ - /* one ux is multiplied in to bessel, complex and psi at once, not here */ - F[0] += ux*mod_psi_odd(0, uxx); - - uxrho_2nm2 = 1.0; - for (n = 1;n < n_modPsi; n++) { - mpe = mod_psi_even(n, uxx); - mpo = mod_psi_odd(n, uxx); - uxrho_2n = uxrho_2nm2*uxrho2; - - F[0] += ux *uxrho_2n *mpo; - F[1] += 2*n*ux2*uxrho_2nm2*mpe*d[1]; - F[2] += 2*n*ux2*uxrho_2nm2*mpe*d[2]; - - /* y < rho => ux2*uxrho_2nm2*d[1] < ux*uxrho_2n */ - if (fabs(2*n*ux*uxrho_2n*mpe) < part_error) - break; + for (n = 0; n < end; n++) { + F[1] -= bon.e[n]*ztn_i; + F[2] += bon.e[n]*ztn_r; - uxrho_2nm2 = uxrho_2n; - } - // fprintf(stderr, " psi force %f %f %f %d\n", F[0], F[1], F[2], n); + tmp_r = ztn_r*zet2_r - ztn_i*zet2_i; + ztn_i = ztn_r*zet2_i + ztn_i*zet2_r; + ztn_r = tmp_r; } + // fprintf(stderr, "complex force %f %f %f %d\n", F[0], F[1], F[2], end); + } + /* psi sum */ + { + int n; + double uxx = ux*d[0]; + double uxrho2 = ux2*rho2; + double uxrho_2n, uxrho_2nm2; /* rho^{2n-2} */ + double mpe, mpo; - for (i = 0; i < 3; i++) - F[i] *= ux; + /* n = 0 inflicts only Fx and pot */ + /* one ux is multiplied in to bessel, complex and psi at once, not here */ + F[0] += ux*mod_psi_odd(0, uxx); - /* explicitly added potentials r_{-1,0} and r_{1,0} */ - { - double cx = d[0] + box_l[0]; - double rinv2 = 1.0/(cx*cx + rho2), rinv = sqrt(rinv2); - double rinv3 = rinv*rinv2; - F[0] += cx*rinv3; - F[1] += d[1]*rinv3; - F[2] += d[2]*rinv3; + uxrho_2nm2 = 1.0; + for (n = 1;n < n_modPsi; n++) { + mpe = mod_psi_even(n, uxx); + mpo = mod_psi_odd(n, uxx); + uxrho_2n = uxrho_2nm2*uxrho2; - cx = d[0] - box_l[0]; - rinv2 = 1.0/(cx*cx + rho2); rinv = sqrt(rinv2); - rinv3 = rinv*rinv2; - F[0] += cx*rinv3; - F[1] += d[1]*rinv3; - F[2] += d[2]*rinv3; + F[0] += ux *uxrho_2n *mpo; + F[1] += 2*n*ux2*uxrho_2nm2*mpe*d[1]; + F[2] += 2*n*ux2*uxrho_2nm2*mpe*d[2]; - rinv3 = 1/(dl2*dl); - F[0] += d[0]*rinv3; - F[1] += d[1]*rinv3; - F[2] += d[2]*rinv3; + /* y < rho => ux2*uxrho_2nm2*d[1] < ux*uxrho_2n */ + if (fabs(2*n*ux*uxrho_2n*mpe) < part_error) + break; - // fprintf(stderr, "explcit force %f %f %f\n", F[0], F[1], F[2]); + uxrho_2nm2 = uxrho_2n; } + // fprintf(stderr, " psi force %f %f %f %d\n", F[0], F[1], F[2], n); + } + - for (i = 0; i < 3; i++) - force[i] += pref*F[i]; + for (i = 0; i < 3; i++) + F[i] *= ux; + + /* explicitly added potentials r_{-1,0} and r_{1,0} */ + { + double cx = d[0] + box_l[0]; + double rinv2 = 1.0/(cx*cx + rho2), rinv = sqrt(rinv2); + double rinv3 = rinv*rinv2; + F[0] += cx*rinv3; + F[1] += d[1]*rinv3; + F[2] += d[2]*rinv3; + + cx = d[0] - box_l[0]; + rinv2 = 1.0/(cx*cx + rho2); rinv = sqrt(rinv2); + rinv3 = rinv*rinv2; + F[0] += cx*rinv3; + F[1] += d[1]*rinv3; + F[2] += d[2]*rinv3; + + rinv3 = 1/(dl2*dl); + F[0] += d[0]*rinv3; + F[1] += d[1]*rinv3; + F[2] += d[2]*rinv3; + + // fprintf(stderr, "explcit force %f %f %f\n", F[0], F[1], F[2]); } + + for (i = 0; i < 3; i++) + force[i] += pref*F[i]; } MDINLINE double calc_mmm2d_copy_pair_energy(double d[3])