59 for (
int ituple=0; ituple<ntuple; ++ituple ) {
72 #ifdef CALCULATE_THOLE_CORRECTION
107 BigReal polyauaa = 1 + 0.5*auaa;
108 BigReal polyauad = 1 + 0.5*auad;
109 BigReal polyauda = 1 + 0.5*auda;
110 BigReal polyaudd = 1 + 0.5*audd;
114 ethole += qq * raa_invlen * (1 - polyauaa * expauaa);
115 ethole += -qq * rad_invlen * (1 - polyauad * expauad);
116 ethole += -qq * rda_invlen * (1 - polyauda * expauda);
117 ethole += qq * rdd_invlen * (1 - polyaudd * expaudd);
119 polyauaa = 1 + auaa*polyauaa;
120 polyauad = 1 + auad*polyauad;
121 polyauda = 1 + auda*polyauda;
122 polyaudd = 1 + audd*polyaudd;
124 BigReal raa_invlen3 = raa_invlen * raa_invlen * raa_invlen;
125 BigReal rad_invlen3 = rad_invlen * rad_invlen * rad_invlen;
126 BigReal rda_invlen3 = rda_invlen * rda_invlen * rda_invlen;
127 BigReal rdd_invlen3 = rdd_invlen * rdd_invlen * rdd_invlen;
130 BigReal dfaa = qq * raa_invlen3 * (polyauaa*expauaa - 1);
131 BigReal dfad = -qq * rad_invlen3 * (polyauad*expauad - 1);
132 BigReal dfda = -qq * rda_invlen3 * (polyauda*expauda - 1);
133 BigReal dfdd = qq * rdd_invlen3 * (polyaudd*expaudd - 1);
142 if ( simParams->
alchOn ) {
147 mol->get_fep_type(
atomID[0]));
149 mol->get_fep_type(
atomID[2]));
152 if ( 0 < typeSum && typeSum < order ) {
155 ethole *= elec_lambda_1;
156 faa *= elec_lambda_1;
157 fad *= elec_lambda_1;
158 fda *= elec_lambda_1;
159 fdd *= elec_lambda_1;
160 }
else if ( 0 > typeSum && typeSum > -order ) {
163 ethole *= elec_lambda_2;
164 faa *= elec_lambda_2;
165 fad *= elec_lambda_2;
166 fda *= elec_lambda_2;
167 fdd *= elec_lambda_2;
177 DebugM(3,
"::computeForce() -- ending with delta energy " << ethole
181 reduction[virialIndex_XX] += faa.x * raa.x + fad.x * rad.x
182 + fda.x * rda.x + fdd.x * rdd.x;
183 reduction[virialIndex_XY] += faa.x * raa.y + fad.x * rad.y
184 + fda.x * rda.y + fdd.x * rdd.y;
185 reduction[virialIndex_XZ] += faa.x * raa.z + fad.x * rad.z
186 + fda.x * rda.z + fdd.x * rdd.z;
187 reduction[virialIndex_YX] += faa.y * raa.x + fad.y * rad.x
188 + fda.y * rda.x + fdd.y * rdd.x;
189 reduction[virialIndex_YY] += faa.y * raa.y + fad.y * rad.y
190 + fda.y * rda.y + fdd.y * rdd.y;
191 reduction[virialIndex_YZ] += faa.y * raa.z + fad.y * rad.z
192 + fda.y * rda.z + fdd.y * rdd.z;
193 reduction[virialIndex_ZX] += faa.z * raa.x + fad.z * rad.x
194 + fda.z * rda.x + fdd.z * rdd.x;
195 reduction[virialIndex_ZY] += faa.z * raa.y + fad.z * rad.y
196 + fda.z * rda.y + fdd.z * rdd.y;
197 reduction[virialIndex_ZZ] += faa.z * raa.z + fad.z * rad.z
198 + fda.z * rda.z + fdd.z * rdd.z;
200 if (pressureProfileData) {
219 pai, paj, pn, faa.x * raa.x, faa.y * raa.y, faa.z * raa.z,
220 pressureProfileData);
222 pai, pdj, pn, fad.x * rad.x, fad.y * rad.y, fad.z * rad.z,
223 pressureProfileData);
225 pdi, paj, pn, fda.x * rda.x, fda.y * rda.y, fda.z * rda.z,
226 pressureProfileData);
228 pdi, pdj, pn, fdd.x * rdd.x, fdd.y * rdd.y, fdd.z * rdd.z,
229 pressureProfileData);
void pp_reduction(int nslabs, int n1, int n2, int atype1, int atype2, int numtypes, BigReal vxx, BigReal vyy, BigReal vzz, BigReal *reduction)
SimParameters * simParameters
static int pressureProfileSlabs
void pp_clamp(int &n, int nslabs)
BigReal getCurrentLambda2(const int)
Vector delta(const Position &pos1, const Position &pos2) const
static BigReal pressureProfileMin
unsigned char get_fep_type(int anum) const
BigReal getCurrentLambda(const int)
static int pressureProfileAtomTypes
static BigReal pressureProfileThickness
BigReal getElecLambda(const BigReal)