ComputeNonbondedTI.C File Reference

#include "ComputeNonbondedInl.h"
#include "ComputeNonbondedAlch.h"
#include "ComputeNonbondedBase.h"

Go to the source code of this file.

Defines

#define TIFLAG
#define CALCENERGY
#define NBTYPE   NBPAIR
#define FULLELECT
#define MERGEELECT
#define SLOWONLY
#define NBTYPE   NBSELF
#define FULLELECT
#define MERGEELECT
#define SLOWONLY
#define NBTYPE   NBPAIR
#define FULLELECT
#define MERGEELECT
#define SLOWONLY
#define NBTYPE   NBSELF
#define FULLELECT
#define MERGEELECT
#define SLOWONLY

Functions

void ti_vdw_force_energy_dUdl (BigReal A, BigReal B, BigReal r2, BigReal myVdwShift, BigReal switchdist2, BigReal cutoff2, BigReal switchfactor, Bool vdwForceSwitching, BigReal myVdwLambda, BigReal alchVdwShiftCoeff, Bool alchWCAOn, BigReal myRepLambda, BigReal *alch_vdw_energy, BigReal *alch_vdw_force, BigReal *alch_vdw_dUdl)

Define Documentation

#define CALCENERGY

Definition at line 144 of file ComputeNonbondedTI.C.

#define FULLELECT

Definition at line 189 of file ComputeNonbondedTI.C.

#define FULLELECT

Definition at line 189 of file ComputeNonbondedTI.C.

#define FULLELECT

Definition at line 189 of file ComputeNonbondedTI.C.

#define FULLELECT

Definition at line 189 of file ComputeNonbondedTI.C.

#define MERGEELECT

Definition at line 191 of file ComputeNonbondedTI.C.

#define MERGEELECT

Definition at line 191 of file ComputeNonbondedTI.C.

#define MERGEELECT

Definition at line 191 of file ComputeNonbondedTI.C.

#define MERGEELECT

Definition at line 191 of file ComputeNonbondedTI.C.

#define NBTYPE   NBSELF

Definition at line 187 of file ComputeNonbondedTI.C.

#define NBTYPE   NBPAIR

Definition at line 187 of file ComputeNonbondedTI.C.

#define NBTYPE   NBSELF

Definition at line 187 of file ComputeNonbondedTI.C.

#define NBTYPE   NBPAIR

Definition at line 187 of file ComputeNonbondedTI.C.

#define SLOWONLY

Definition at line 194 of file ComputeNonbondedTI.C.

#define SLOWONLY

Definition at line 194 of file ComputeNonbondedTI.C.

#define SLOWONLY

Definition at line 194 of file ComputeNonbondedTI.C.

#define SLOWONLY

Definition at line 194 of file ComputeNonbondedTI.C.

#define TIFLAG

Definition at line 143 of file ComputeNonbondedTI.C.


Function Documentation

void ti_vdw_force_energy_dUdl ( BigReal  A,
BigReal  B,
BigReal  r2,
BigReal  myVdwShift,
BigReal  switchdist2,
BigReal  cutoff2,
BigReal  switchfactor,
Bool  vdwForceSwitching,
BigReal  myVdwLambda,
BigReal  alchVdwShiftCoeff,
Bool  alchWCAOn,
BigReal  myRepLambda,
BigReal alch_vdw_energy,
BigReal alch_vdw_force,
BigReal alch_vdw_dUdl 
) [inline]

Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved.

Definition at line 17 of file ComputeNonbondedTI.C.

References vdw_forceandenergy(), vdw_fswitch_forceandenergy(), vdw_fswitch_shift(), and vdw_switch().

00022                             {
00023   BigReal U, F, dU, switchmul, switchmul2;
00024   /*
00025    * The variable naming here is unavoidably awful. A number immediately after
00026    * a variable implies a distance to that power. The suffix "_2" indicates an
00027    * evaluation at alchLambda2. The prefix "inv" means the inverse (previously
00028    * this also used underscores - it was awful).
00029    */
00030   if (alchWCAOn) {
00031     /* THIS BRANCH IS PARTIALLY INCOMPLETE AND BLOCKED INSIDE SimParameters
00032      *
00033      * The problem is that WCA creates TWO energy terms and thus two different
00034      * TI derivatives, but the reduction code is currently only set up for
00035      * one derivative. There is no such problem with FEP because the energy
00036      * is not decomposed. In principle, the contribution to the free energy
00037      * from each derivative is zero while the other is changing (e.g. while
00038      * repulsion is changing dispersion is not), but the derivatives are still
00039      * non-zero and it is NAMD convention to report them even as such.
00040      * However, even if we were willing to break this convention, both
00041      * derivatives still compute to the free energy at the boundary where
00042      * repulsion is completely on and dispersion is exactly zero - two
00043      * derivatives are thus required for a complete working implementation.
00044      *
00045      */
00046 
00047     // WCA-on, auxilliary, lambda-dependent cutoff based on Rmin
00048     //
00049     // Avoid divide by zero - correctly zeroes interaction below.
00050     const BigReal Rmin2 = (B <= 0.0 ? 0.0 : powf(2.0*A/B, 1.f/3));
00051     if (myRepLambda < 1.0) {
00052       // modified repulsive soft-core
00053       // All of the scaling is baked into the shift and cutoff values. There
00054       // are no linear coupling terms out front.
00055       const BigReal WCAshift = Rmin2*(1 - myRepLambda)*(1 - myRepLambda);
00056       if (r2 <= Rmin2 - WCAshift) {
00057         const BigReal epsilon = B*B/(4.0*A);
00058         vdw_forceandenergy(A, B, r2 + WCAshift, &U, &F);
00059         *alch_vdw_energy = U + epsilon;
00060         *alch_vdw_force = F;
00061         *alch_vdw_dUdl = Rmin2*(1 - myRepLambda)*F;
00062       } else {
00063         *alch_vdw_energy = 0.0;
00064         *alch_vdw_force = 0.0;
00065         *alch_vdw_dUdl = 0.0;
00066       }
00067     } else { // myRepLambda == 1.0
00068       if (vdwForceSwitching) {
00069         // force and potential switching
00070         if (r2 <= Rmin2) {
00071           // normal LJ, potential w/two shifts 
00072           const BigReal epsilon = B*B/(4.0*A); 
00073           vdw_fswitch_shift(A, B, switchdist2, cutoff2, &dU);
00074           vdw_forceandenergy(A, B, r2, &U, &F);
00075           *alch_vdw_energy = U + (1 - myVdwLambda)*epsilon + myVdwLambda*dU;
00076           *alch_vdw_force = F;
00077           *alch_vdw_dUdl = dU - epsilon;
00078         } else if (r2 <= switchdist2) {
00079           // normal LJ potential w/shift
00080           vdw_fswitch_shift(A, B, switchdist2, cutoff2, &dU);
00081           vdw_forceandenergy(A, B, r2, &U, &F);
00082           *alch_vdw_energy = myVdwLambda*(U + dU);
00083           *alch_vdw_force = myVdwLambda*F;
00084           *alch_vdw_dUdl = U + dU;
00085         } else { // r2 > switchdist
00086           // normal LJ potential with linear coupling
00087           vdw_fswitch_forceandenergy(A, B, r2, switchdist2, cutoff2, &U, &F);
00088           *alch_vdw_energy = myVdwLambda*U;
00089           *alch_vdw_force = myVdwLambda*F;
00090           *alch_vdw_dUdl = U;
00091         }
00092       } else {
00093         // potential switching - also correct for no switching
00094         if (r2 <= Rmin2) {
00095           // normal LJ potential w/shift
00096           const BigReal epsilon = B*B/(4.0*A);
00097           vdw_forceandenergy(A, B, r2, &U, &F);
00098           *alch_vdw_energy = U + (1 - myVdwLambda)*epsilon;
00099           *alch_vdw_force = F;
00100           *alch_vdw_dUdl = -epsilon;
00101         } else { // r2 > Rmin2
00102           // normal LJ potential with linear coupling
00103           vdw_switch(r2, switchdist2, cutoff2, switchfactor, &switchmul, \
00104               &switchmul2);
00105           vdw_forceandenergy(A, B, r2, &U, &F);
00106           *alch_vdw_energy = myVdwLambda*switchmul*U;
00107           *alch_vdw_force = myVdwLambda*(switchmul*F + switchmul2*U);
00108           *alch_vdw_dUdl = switchmul*U;
00109         }
00110       }
00111     }
00112   } else { // WCA-off
00113     if (vdwForceSwitching) {
00114       // force and potential switching
00115       if (r2 <= switchdist2) {
00116         // normal LJ potential w/shift
00117         vdw_forceandenergy(A, B, r2 + myVdwShift, &U, &F);
00118         vdw_fswitch_shift(A, B, switchdist2 + myVdwShift, cutoff2, &dU);
00119         *alch_vdw_energy = myVdwLambda*(U + dU);
00120         *alch_vdw_force = myVdwLambda*F;
00121         *alch_vdw_dUdl = U + 0.5*myVdwLambda*alchVdwShiftCoeff*F + dU;
00122       } else { // r2 > switchdist2
00123         vdw_fswitch_forceandenergy(A, B, r2 + myVdwShift, \
00124             switchdist2 + myVdwShift, cutoff2, &U, &F);
00125         *alch_vdw_energy = myVdwLambda*U;
00126         *alch_vdw_force = myVdwLambda*F;
00127         *alch_vdw_dUdl = U + 0.5*myVdwLambda*alchVdwShiftCoeff*F;
00128       }
00129     } else {
00130       // potential switching - also correct for no switching
00131       vdw_switch(r2, switchdist2, cutoff2, switchfactor, &switchmul, \
00132           &switchmul2);
00133       vdw_forceandenergy(A, B, r2 + myVdwShift, &U, &F);
00134       *alch_vdw_energy = myVdwLambda*switchmul*U;
00135       *alch_vdw_force = myVdwLambda*(switchmul*F + switchmul2*U);
00136       *alch_vdw_dUdl = switchmul*(U + 0.5*myVdwLambda*alchVdwShiftCoeff*F);
00137     }
00138   }
00139 }


Generated on 20 Sep 2019 for NAMD by  doxygen 1.6.1