00001
00007 #include "InfoStream.h"
00008 #include "ComputeNonbondedCUDAExcl.h"
00009 #include "Molecule.h"
00010 #include "Parameters.h"
00011 #include "LJTable.h"
00012 #include "Node.h"
00013 #include "ReductionMgr.h"
00014 #include "Lattice.h"
00015 #include "PressureProfile.h"
00016 #include "Debug.h"
00017
00018
00019
00020 int ExclElem::pressureProfileSlabs = 0;
00021 int ExclElem::pressureProfileAtomTypes = 1;
00022 BigReal ExclElem::pressureProfileThickness = 0;
00023 BigReal ExclElem::pressureProfileMin = 0;
00024
00025 void ExclElem::getMoleculePointers
00026 (Molecule* mol, int* count, int32*** byatom, Exclusion** structarray)
00027 {
00028 #ifdef MEM_OPT_VERSION
00029 NAMD_die("Should not be called in ExclElem::getMoleculePointers in memory optimized version!");
00030 #else
00031 *count = mol->numExclusions;
00032 *byatom = mol->exclusionsByAtom;
00033 *structarray = mol->exclusions;
00034 #endif
00035 }
00036
00037 void ExclElem::getParameterPointers(Parameters *p, const int **v) {
00038 *v = 0;
00039 }
00040
00041 void ExclElem::computeForce(BigReal *reduction,
00042 BigReal *pressureProfileData)
00043 {
00044 const CompAtom &p_i = p[0]->x[localIndex[0]];
00045 const CompAtom &p_j = p[1]->x[localIndex[1]];
00046
00047
00048 const Lattice & lattice = p[0]->p->lattice;
00049 const Vector r12 = lattice.delta(p_i.position, p_j.position);
00050 BigReal r2 = r12.length2();
00051
00052 if ( r2 > cutoff2 ) return;
00053
00054 r2 += r2_delta;
00055
00056 union { double f; int64 i; } r2i;
00057 r2i.f = r2;
00058 const int r2_delta_expc = 64 * (r2_delta_exp - 1023);
00059 int table_i = (r2i.i >> (32+14)) + r2_delta_expc;
00060
00061 const BigReal* const table_four_i = table_noshort + 16*table_i;
00062
00063 BigReal diffa = r2 - r2_table[table_i];
00064
00065 const int doFull = p[0]->p->flags.doFullElectrostatics;
00066
00067 BigReal fast_a, fast_b, fast_c, fast_d;
00068 BigReal slow_a, slow_b, slow_c, slow_d;
00069
00070 if ( modified ) {
00071
00072 const LJTable::TableEntry * lj_pars =
00073 ljTable->table_row(p_i.vdwType) + 2 * p_j.vdwType;
00074
00075
00076 const BigReal A = scaling * ( (lj_pars+1)->A - lj_pars->A );
00077 const BigReal B = scaling * ( (lj_pars+1)->B - lj_pars->B );
00078
00079 BigReal vdw_d = A * table_four_i[0] - B * table_four_i[2];
00080 BigReal vdw_c = A * table_four_i[1] - B * table_four_i[3];
00081 BigReal vdw_b = A * table_four_i[4] - B * table_four_i[6];
00082 BigReal vdw_a = A * table_four_i[5] - B * table_four_i[7];
00083
00084 const BigReal kqq = (1.0 - scale14) *
00085 COLOUMB * p_i.charge * p_j.charge * scaling * dielectric_1;
00086
00087 fast_a = vdw_a + kqq * fast_table[4*table_i+0];
00088 fast_b = vdw_b + 2. * kqq * fast_table[4*table_i+1];
00089 fast_c = vdw_c + 4. * kqq * fast_table[4*table_i+2];
00090 fast_d = vdw_d + 6. * kqq * fast_table[4*table_i+3];
00091
00092 if ( doFull ) {
00093 slow_a = kqq * slow_table[4*table_i+0];
00094 slow_b = 2. * kqq * slow_table[4*table_i+1];
00095 slow_c = 4. * kqq * slow_table[4*table_i+2];
00096 slow_d = 6. * kqq * slow_table[4*table_i+3];
00097 }
00098
00099 } else if ( doFull ) {
00100
00101 const BigReal kqq =
00102 COLOUMB * p_i.charge * p_j.charge * scaling * dielectric_1;
00103
00104 slow_d = kqq * ( table_four_i[8] - table_four_i[12] );
00105 slow_c = kqq * ( table_four_i[9] - table_four_i[13] );
00106 slow_b = kqq * ( table_four_i[10] - table_four_i[14] );
00107 slow_a = kqq * ( table_four_i[11] - table_four_i[15] );
00108
00109 }
00110
00111 register BigReal fast_dir =
00112 (diffa * fast_d + fast_c) * diffa + fast_b;
00113
00114 const Force f12 = fast_dir * r12;
00115
00116
00117 p[0]->r->f[Results::nbond][localIndex[0]] += f12;
00118 p[1]->r->f[Results::nbond][localIndex[1]] -= f12;
00119
00120
00121 reduction[virialIndex_XX] += f12.x * r12.x;
00122
00123
00124
00125 reduction[virialIndex_YY] += f12.y * r12.y;
00126
00127
00128
00129 reduction[virialIndex_ZZ] += f12.z * r12.z;
00130
00131 if ( doFull ) {
00132 register BigReal slow_dir =
00133 (diffa * slow_d + slow_c) * diffa + slow_b;
00134
00135 const Force slow_f12 = slow_dir * r12;
00136
00137 p[0]->r->f[Results::slow][localIndex[0]] += slow_f12;
00138 p[1]->r->f[Results::slow][localIndex[1]] -= slow_f12;
00139
00140
00141 reduction[slowVirialIndex_XX] += slow_f12.x * r12.x;
00142
00143
00144
00145 reduction[slowVirialIndex_YY] += slow_f12.y * r12.y;
00146
00147
00148
00149 reduction[slowVirialIndex_ZZ] += slow_f12.z * r12.z;
00150 }
00151
00152 }
00153
00154 void ExclElem::submitReductionData(BigReal *data, SubmitReduction *reduction)
00155 {
00156
00157 ADD_TENSOR(reduction,REDUCTION_VIRIAL_NBOND,data,virialIndex);
00158 ADD_TENSOR(reduction,REDUCTION_VIRIAL_SLOW,data,slowVirialIndex);
00159 }
00160