Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

ComputeNonbondedBase2.h

Go to the documentation of this file.
00001 
00007 EXCLUDED( FAST( foo bar ) )
00008 EXCLUDED( MODIFIED( foo bar ) )
00009 EXCLUDED( NORMAL( foo bar ) )
00010 NORMAL( MODIFIED( foo bar ) )
00011 ALCHPAIR( NOT_ALCHPAIR( foo bar ) )
00012 
00013 ALCHPAIR(
00014   // get alchemical nonbonded scaling parameters (once per pairlist)
00015   myLambda = ALCH1(lambdaUp) ALCH2(lambdaDown);
00016   FEP(myLambda2 = ALCH1(lambda2Up) ALCH2(lambda2Down);)
00017   myElecLambda =  ALCH1(elecLambdaUp) ALCH2(elecLambdaDown); 
00018   FEP(myElecLambda2 = ALCH1(elecLambda2Up) ALCH2(elecLambda2Down);)
00019   myVdwLambda =  ALCH1(vdwLambdaUp) ALCH2(vdwLambdaDown); 
00020   FEP(myVdwLambda2 = ALCH1(vdwLambda2Up) ALCH2(vdwLambda2Down);) 
00021   myVdwShift =  ALCH1(vdwShiftUp) ALCH2(vdwShiftDown); 
00022   FEP(myVdwShift2 =  ALCH1(vdwShift2Up) ALCH2(vdwShift2Down);) 
00023 )
00024 
00025 #ifdef ARCH_POWERPC
00026      __alignx(16, table_four);
00027      __alignx(16, p_1);
00028 #pragma unroll(1)
00029 #endif
00030 
00031 #ifndef ARCH_POWERPC
00032 #pragma ivdep
00033 #endif
00034 
00035 
00036   
00037     for (k=0; k<npairi; ++k) {      
00038       int table_i = (r2iilist[2*k] >> 14) + r2_delta_expc;  // table_i >= 0 
00039       const int j = pairlisti[k];
00040       register const CompAtom *p_j = p_1 + j;
00041       
00042       BigReal diffa = r2list[k] - r2_table[table_i];
00043       const BigReal* const table_four_i = table_four + 16*table_i;
00044 
00045       FAST(
00046       const LJTable::TableEntry * lj_pars = 
00047               lj_row + 2 * vdwtype_array[j]  MODIFIED(+ 1);
00048       ) 
00049 
00050 #if ( SHORT( FAST( 1+ ) ) 0 ) 
00051       Force *f_j = f_1 + j;
00052 #endif
00053         
00054 #if ( FULL( 1+ ) 0 )
00055       Force *fullf_j = fullf_1 + j;
00056 #endif
00057 
00058       //Power PC aliasing and alignment constraints
00059 #ifdef ARCH_POWERPC
00060       __alignx(16, table_four_i);
00061       FAST (
00062       __alignx(16, lj_pars);
00063       )
00064       __alignx(16, p_j);
00065       
00066 #if ( FULL( 1+ ) 0 )
00067 #pragma disjoint (*table_four_i, *fullf_j)
00068 #pragma disjoint (*p_j,          *fullf_j)
00069 #if ( SHORT( FAST( 1+ ) ) 0 ) 
00070 #pragma disjoint (*f_j    , *fullf_j)
00071 #pragma disjoint (*fullf_j, *f_j)
00072 #endif   //Short + fast
00073 #endif   //Full
00074 
00075 #if ( SHORT( FAST( 1+ ) ) 0 ) 
00076 #pragma disjoint (*table_four_i, *f_j)
00077 #pragma disjoint (*p_j,          *f_j)
00078 #pragma disjoint (*lj_pars,      *f_j)
00079       __prefetch_by_load ((void *)&f_j->x);
00080 #endif //Short + Fast
00081 
00082 #endif   //ARCH_POWERPC
00083 
00084       /*
00085       BigReal modf = 0.0;
00086       int atom2 = p_j->id;
00087       register char excl_flag = ( (atom2 >= excl_min && atom2 <= excl_max) ?
00088                                         excl_flags[atom2-excl_min] : 0 );
00089       if ( excl_flag ) { ++exclChecksum; }
00090       SELF( if ( j < j_hgroup ) { excl_flag = EXCHCK_FULL; } )
00091       if ( excl_flag ) {
00092         if ( excl_flag == EXCHCK_FULL ) {
00093           lj_pars = lj_null_pars;
00094           modf = 1.0;
00095         } else {
00096           ++lj_pars;
00097           modf = modf_mod;
00098         }
00099       }
00100       */
00101 
00102       BigReal kqq = kq_i * p_j->charge;
00103 
00104       LES( BigReal lambda_pair = lambda_table_i[p_j->partition]; )
00105 
00106       register const BigReal p_ij_x = p_i_x - p_j->position.x;
00107       register const BigReal p_ij_y = p_i_y - p_j->position.y;
00108       register const BigReal p_ij_z = p_i_z - p_j->position.z;
00109       
00110 #if ( FAST(1+) 0 )
00111       const BigReal A = scaling * lj_pars->A;
00112       const BigReal B = scaling * lj_pars->B;
00113 
00114       BigReal vdw_d = A * table_four_i[0] - B * table_four_i[2];
00115       BigReal vdw_c = A * table_four_i[1] - B * table_four_i[3];
00116       BigReal vdw_b = A * table_four_i[4] - B * table_four_i[6];
00117       BigReal vdw_a = A * table_four_i[5] - B * table_four_i[7];
00118 
00119       ALCHPAIR (
00120         // Alchemical free energy calculation
00121         // Pairlists are separated so that lambda-coupled pairs are handled
00122         // independently from normal nonbonded (inside ALCHPAIR macro).
00123         // The separation-shifted van der Waals potential and a shifted 
00124         // electrostatics potential for decoupling are calculated explicitly.
00125         // Would be faster with lookup tables but because only a small minority
00126         // of nonbonded pairs are lambda-coupled the impact is minimal. 
00127         // Explicit calculation also makes things easier to modify.
00128         
00129         const BigReal r2 = r2list[k] - r2_delta;
00130 
00131         // These are now inline functions (in ComputeNonbondedFep.C) to 
00132         // tidy the code
00133 
00134         FEP(fep_vdw_forceandenergies(A,B,r2,myVdwShift,myVdwShift2,switchdist2,
00135           cutoff2, myVdwLambda, myVdwLambda2, switchfactor, &alch_vdw_energy, 
00136           &alch_vdw_force,&alch_vdw_energy_2);)
00137         TI(ti_vdw_force_energy_dUdl(A,B,r2,myVdwShift,switchdist2,
00138           cutoff2, myVdwLambda, fepVdwShiftCoeff, switchfactor, 
00139           &alch_vdw_energy, &alch_vdw_force, &alch_vdw_dUdl);)
00140       )
00141       
00142       NOT_ALCHPAIR(ENERGY(
00143         register BigReal vdw_val =
00144           ( ( diffa * vdw_d * (1/6.)+ vdw_c * (1/4.)) * diffa + vdw_b *(1/2.)) * diffa + vdw_a;
00145         vdwEnergy -= LAM(lambda_pair *) vdw_val;
00146         FEP(vdwEnergy_s -= vdw_val;)
00147       ))
00148 
00149       ALCHPAIR(
00150         ENERGY(vdwEnergy   += alch_vdw_energy;)
00151         FEP(vdwEnergy_s += alch_vdw_energy_2;)
00152         TI(ALCH1(vdwEnergy_ti_1) ALCH2(vdwEnergy_ti_2) += alch_vdw_dUdl;)
00153       ) // ALCHPAIR
00154       
00155 #endif // FAST
00156 
00157 #if ( FAST(1+) 0 )
00158       INT( 
00159       register BigReal vdw_dir =
00160       ( diffa * vdw_d + vdw_c ) * diffa + vdw_b;
00161       //BigReal force_r =  LAM(lambda_pair *) vdw_dir;
00162       reduction[pairVDWForceIndex_X] += force_sign * vdw_dir * p_ij_x;
00163       reduction[pairVDWForceIndex_Y] += force_sign * vdw_dir * p_ij_y;
00164       reduction[pairVDWForceIndex_Z] += force_sign * vdw_dir * p_ij_z;
00165       )
00166 
00167 #if ( SHORT(1+) 0 ) // Short-range electrostatics
00168 
00169       NORMAL(
00170       BigReal fast_d = kqq * table_four_i[8];
00171       BigReal fast_c = kqq * table_four_i[9];
00172       BigReal fast_b = kqq * table_four_i[10];
00173       BigReal fast_a = kqq * table_four_i[11];
00174       )
00175       MODIFIED(
00176       BigReal modfckqq = (1.0-modf_mod) * kqq;
00177       BigReal fast_d = modfckqq * table_four_i[8];
00178       BigReal fast_c = modfckqq * table_four_i[9];
00179       BigReal fast_b = modfckqq * table_four_i[10];
00180       BigReal fast_a = modfckqq * table_four_i[11];
00181       )
00182     
00183       {
00184       ENERGY(
00185       register BigReal fast_val =
00186         ( ( diffa * fast_d * (1/6.)+ fast_c * (1/4.)) * diffa + fast_b *(1/2.)) * diffa + fast_a;
00187       NOT_ALCHPAIR (
00188         electEnergy -=  LAM(lambda_pair *) fast_val;
00189         FEP(electEnergy_s -= fast_val;)
00190       )
00191       ) //ENERGY
00192       ALCHPAIR(
00193         ENERGY(electEnergy   -= myElecLambda * fast_val;)
00194         FEP(electEnergy_s -= myElecLambda2 * fast_val;)
00195         TI(
00196           NOENERGY(register BigReal fast_val = 
00197             ( ( diffa * fast_d * (1/6.)+ fast_c * (1/4.)) * diffa + fast_b *(1/2.)) * diffa + fast_a;)
00198           ALCH1(electEnergy_ti_1) ALCH2(electEnergy_ti_2)  -= fast_val;
00199         )
00200       )
00201 
00202       INT(
00203       register BigReal fast_dir =
00204       ( diffa * fast_d + fast_c ) * diffa + fast_b;
00205       // force_r -= -1.0 * LAM(lambda_pair *) fast_dir;
00206       reduction[pairElectForceIndex_X] +=  force_sign * fast_dir * p_ij_x;
00207       reduction[pairElectForceIndex_Y] +=  force_sign * fast_dir * p_ij_y;
00208       reduction[pairElectForceIndex_Z] +=  force_sign * fast_dir * p_ij_z;
00209       )
00210       }
00211 
00212       // Combined short-range electrostatics and VdW force:
00213       NOT_ALCHPAIR(
00214         fast_d += vdw_d;
00215         fast_c += vdw_c;
00216         fast_b += vdw_b;
00217         fast_a += vdw_a;  // not used!
00218       )
00219       register BigReal fast_dir =
00220                   (diffa * fast_d + fast_c) * diffa + fast_b;
00221       BigReal force_r =  LAM(lambda_pair *) fast_dir;
00222       ALCHPAIR(
00223         force_r *= myElecLambda; 
00224         force_r += alch_vdw_force;
00225         // special ALCH forces already multiplied by relevant lambda
00226       )
00227           
00228       register BigReal tmp_x = force_r * p_ij_x;
00229       PAIR( virial_xx += tmp_x * p_ij_x; )
00230       PAIR( virial_xy += tmp_x * p_ij_y; )
00231       PAIR( virial_xz += tmp_x * p_ij_z; )
00232 
00233       f_i_x += tmp_x;
00234       f_j->x -= tmp_x;
00235 
00236       register BigReal tmp_y = force_r * p_ij_y;
00237       PAIR( virial_yy += tmp_y * p_ij_y; )
00238       PAIR( virial_yz += tmp_y * p_ij_z; )
00239       f_i_y += tmp_y;
00240       f_j->y -= tmp_y;
00241       
00242       register BigReal tmp_z = force_r * p_ij_z;
00243       PAIR( virial_zz += tmp_z * p_ij_z; )
00244       f_i_z += tmp_z;
00245       f_j->z -= tmp_z;
00246 
00247       PPROF(
00248         const BigReal p_j_z = p_j->position.z;
00249         int n2 = (int)floor((p_j_z-pressureProfileMin)*invThickness);
00250         pp_clamp(n2, pressureProfileSlabs);
00251         int p_j_partition = p_j->partition;
00252 
00253         pp_reduction(pressureProfileSlabs, n1, n2, 
00254                      p_i_partition, p_j_partition, pressureProfileAtomTypes,
00255                      tmp_x*p_ij_x, tmp_y * p_ij_y, tmp_z*p_ij_z,
00256                      pressureProfileReduction);
00257 
00258       )
00259 
00260 #endif // SHORT
00261 #endif // FAST
00262 
00263 #if ( FULL (EXCLUDED( SHORT ( 1+ ) ) ) 0 ) 
00264       const BigReal* const slow_i = slow_table + 4*table_i;
00265 
00266 #ifdef ARCH_POWERPC  //Alignment and aliasing constraints
00267       __alignx (16, slow_i);
00268 #if ( SHORT( FAST( 1+ ) ) 0 ) 
00269 #pragma disjoint (*slow_i, *f_j)
00270 #endif
00271 #pragma disjoint (*slow_i, *fullf_j)
00272 #endif  //ARCH_POWERPC
00273 
00274 #endif //FULL 
00275 
00276 
00277 #if ( FULL (MODIFIED( SHORT ( 1+ ) ) ) 0 ) 
00278       const BigReal* const slow_i = slow_table + 4*table_i;
00279 
00280 #ifdef ARCH_POWERPC //Alignment and aliasing constraints
00281       __alignx (16, slow_i);
00282 #if ( SHORT( FAST( 1+ ) ) 0 ) 
00283 #pragma disjoint (*slow_i, *f_j)
00284 #endif
00285 #pragma disjoint (*slow_i, *fullf_j)
00286 #endif //ARCH_POWERPC
00287 
00288 #endif //FULL
00289       
00290       FULL(
00291       BigReal slow_d = table_four_i[8 SHORT(+ 4)];
00292       BigReal slow_c = table_four_i[9 SHORT(+ 4)];
00293       BigReal slow_b = table_four_i[10 SHORT(+ 4)];
00294       BigReal slow_a = table_four_i[11 SHORT(+ 4)];
00295       EXCLUDED(
00296       SHORT(
00297       slow_a +=    slow_i[0];
00298       slow_b += 2.*slow_i[1];
00299       slow_c += 4.*slow_i[2];
00300       slow_d += 6.*slow_i[3];
00301       )
00302       NOSHORT(
00303       slow_d -= table_four_i[12];
00304       slow_c -= table_four_i[13];
00305       slow_b -= table_four_i[14];
00306       slow_a -= table_four_i[15];
00307       )
00308       )
00309       MODIFIED(
00310       SHORT(
00311       slow_a +=    modf_mod * slow_i[0];
00312       slow_b += 2.*modf_mod * slow_i[1];
00313       slow_c += 4.*modf_mod * slow_i[2];
00314       slow_d += 6.*modf_mod * slow_i[3];
00315       )
00316       NOSHORT(
00317       slow_d -= modf_mod * table_four_i[12];
00318       slow_c -= modf_mod * table_four_i[13];
00319       slow_b -= modf_mod * table_four_i[14];
00320       slow_a -= modf_mod * table_four_i[15];
00321       )
00322       )
00323       slow_d *= kqq;
00324       slow_c *= kqq;
00325       slow_b *= kqq;
00326       slow_a *= kqq;
00327 
00328       ENERGY(
00329       register BigReal slow_val =
00330                 ( ( diffa * slow_d *(1/6.)+ slow_c * (1/4.)) * diffa + slow_b *(1/2.)) * diffa + slow_a;
00331       
00332       NOT_ALCHPAIR (
00333         fullElectEnergy -= LAM(lambda_pair *) slow_val;
00334         FEP(fullElectEnergy_s -= slow_val;) 
00335       )
00336       ) // ENERGY
00337           
00338       ALCHPAIR(
00339         ENERGY(fullElectEnergy   -= myElecLambda * slow_val;)
00340         FEP(fullElectEnergy_s -= myElecLambda2 * slow_val;)
00341         TI(
00342           NOENERGY(register BigReal slow_val =
00343                 ( ( diffa * slow_d *(1/6.)+ slow_c * (1/4.)) * diffa + slow_b *(1/2.)) * diffa + slow_a;)
00344           ALCH1(fullElectEnergy_ti_1) ALCH2(fullElectEnergy_ti_2) -= slow_val;
00345         )
00346       )
00347 
00348       INT( {
00349       register BigReal slow_dir =
00350         ( diffa * slow_d + slow_c ) * diffa + slow_b;
00351       reduction[pairElectForceIndex_X] += force_sign * slow_dir * p_ij_x;
00352       reduction[pairElectForceIndex_Y] += force_sign * slow_dir * p_ij_y;
00353       reduction[pairElectForceIndex_Z] += force_sign * slow_dir * p_ij_z;
00354       } )
00355 
00356 
00357       NOT_ALCHPAIR ( 
00358         FAST(
00359         NOSHORT(
00360         slow_d += vdw_d;
00361         slow_c += vdw_c;
00362         slow_b += vdw_b;
00363         slow_a += vdw_a; // unused!
00364         )
00365         )
00366       )
00367       register BigReal slow_dir = (diffa * slow_d + slow_c) * diffa + slow_b;
00368       BigReal fullforce_r = slow_dir LAM(* lambda_pair);
00369       ALCHPAIR (
00370         fullforce_r *= myElecLambda;
00371         FAST( NOSHORT(
00372           fullforce_r += alch_vdw_force;
00373         ))
00374       )
00375           
00376       {
00377       register BigReal tmp_x = fullforce_r * p_ij_x;
00378       PAIR( fullElectVirial_xx += tmp_x * p_ij_x; )
00379       PAIR( fullElectVirial_xy += tmp_x * p_ij_y; )
00380       PAIR( fullElectVirial_xz += tmp_x * p_ij_z; )
00381       fullf_i_x += tmp_x;
00382       fullf_j->x -= tmp_x;
00383       register BigReal tmp_y = fullforce_r * p_ij_y;
00384       PAIR( fullElectVirial_yy += tmp_y * p_ij_y; )
00385       PAIR( fullElectVirial_yz += tmp_y * p_ij_z; )
00386       fullf_i_y += tmp_y;
00387       fullf_j->y -= tmp_y;
00388       register BigReal tmp_z = fullforce_r * p_ij_z;
00389       PAIR( fullElectVirial_zz += tmp_z * p_ij_z; )
00390       fullf_i_z += tmp_z;
00391       fullf_j->z -= tmp_z;
00392 
00393       PPROF(
00394         const BigReal p_j_z = p_j->position.z;
00395         int n2 = (int)floor((p_j_z-pressureProfileMin)*invThickness);
00396         pp_clamp(n2, pressureProfileSlabs);
00397         int p_j_partition = p_j->partition;
00398 
00399         pp_reduction(pressureProfileSlabs, n1, n2, 
00400                      p_i_partition, p_j_partition, pressureProfileAtomTypes,
00401                      tmp_x*p_ij_x, tmp_y * p_ij_y, tmp_z*p_ij_z,
00402                      pressureProfileReduction);
00403 
00404       )
00405 
00406       }
00407       )
00408 
00409    } // for pairlist
00410 

Generated on Sun Sep 7 04:07:40 2008 for NAMD by  doxygen 1.3.9.1