src/nonbonded/GB.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */

#include <stdio.h>
#include <string.h>

#include <cmath>
#include "nonbonded/GB.hpp"
#include "utils/simError.h"

using namespace std;
namespace OpenMD {

  /* GB is the Gay-Berne interaction for ellipsoidal particles.  The original 
   * paper (for identical uniaxial particles) is:
   *    J. G. Gay and B. J. Berne, J. Chem. Phys., 74, 3316-3319 (1981).
   * A more-general GB potential for dissimilar uniaxial particles:
   *    D. J. Cleaver, C. M. Care, M. P. Allen and M. P. Neal, Phys. Rev. E, 
   *    54, 559-567 (1996).
   * Further parameterizations can be found in:
   *    A. P. J. Emerson, G. R. Luckhurst and S. G. Whatling, Mol. Phys., 
   *    82, 113-124 (1994).
   * And a nice force expression:
   *    G. R. Luckhurst and R. A. Stephens, Liq. Cryst. 8, 451-464 (1990).
   * Even clearer force and torque expressions:
   *    P. A. Golubkov and P. Y. Ren, J. Chem. Phys., 125, 64103 (2006).
   * New expressions for cross interactions of strength parameters:
   *    J. Wu, X. Zhen, H. Shen, G. Li, and P. Ren, J. Chem. Phys., 
   *    135, 155104 (2011).
   *
   * In this version of the GB interaction, each uniaxial ellipsoidal type 
   * is described using a set of 6 parameters:
   *  d:  range parameter for side-by-side (S) and cross (X) configurations
   *  l:  range parameter for end-to-end (E) configuration
   *  epsilon_X:  well-depth parameter for cross (X) configuration
   *  epsilon_S:  well-depth parameter for side-by-side (S) configuration
   *  epsilon_E:  well depth parameter for end-to-end (E) configuration
   *  dw: "softness" of the potential
   * 
   * Additionally, there are two "universal" paramters to govern the overall 
   * importance of the purely orientational (nu) and the mixed
   * orientational / translational (mu) parts of strength of the interactions. 
   * These parameters have default or "canonical" values, but may be changed
   * as a force field option:
   * nu_: purely orientational part : defaults to 1
   * mu_: mixed orientational / translational part : defaults to 2
   */


  GB::GB() : name_("GB"), initialized_(false), mu_(2.0), nu_(1.0), forceField_(NULL) {}
  
  GayBerneParam GB::getGayBerneParam(AtomType* atomType) {
    
    // Do sanity checking on the AtomType we were passed before
    // building any data structures:
    if (!atomType->isGayBerne()) {
      sprintf( painCave.errMsg,
               "GB::getGayBerneParam was passed an atomType (%s) that does\n"
               "\tnot appear to be a Gay-Berne atom.\n",
               atomType->getName().c_str());
      painCave.severity = OPENMD_ERROR;
      painCave.isFatal = 1;
      simError();
    }
    
    DirectionalAtomType* daType = dynamic_cast<DirectionalAtomType*>(atomType);
    GenericData* data = daType->getPropertyByName("GayBerne");
    if (data == NULL) {
      sprintf( painCave.errMsg, "GB::getGayBerneParam could not find\n"
               "\tGay-Berne parameters for atomType %s.\n", 
               daType->getName().c_str());
      painCave.severity = OPENMD_ERROR;
      painCave.isFatal = 1;
      simError(); 
    }
    
    GayBerneParamGenericData* gbData = dynamic_cast<GayBerneParamGenericData*>(data);
    if (gbData == NULL) {
      sprintf( painCave.errMsg,
               "GB::getGayBerneParam could not convert GenericData to\n"
               "\tGayBerneParamGenericData for atom type %s\n", 
               daType->getName().c_str());
      painCave.severity = OPENMD_ERROR;
      painCave.isFatal = 1;
      simError();          
    }
    
    return gbData->getData();
  }

  LJParam GB::getLJParam(AtomType* atomType) {
    
    // Do sanity checking on the AtomType we were passed before
    // building any data structures:
    if (!atomType->isLennardJones()) {
      sprintf( painCave.errMsg,
               "GB::getLJParam was passed an atomType (%s) that does not\n"
               "\tappear to be a Lennard-Jones atom.\n",
               atomType->getName().c_str());
      painCave.severity = OPENMD_ERROR;
      painCave.isFatal = 1;
      simError();
    }
    
    GenericData* data = atomType->getPropertyByName("LennardJones");
    if (data == NULL) {
      sprintf( painCave.errMsg, "GB::getLJParam could not find Lennard-Jones\n"
               "\tparameters for atomType %s.\n", atomType->getName().c_str());
      painCave.severity = OPENMD_ERROR;
      painCave.isFatal = 1;
      simError(); 
    }
    
    LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
    if (ljData == NULL) {
      sprintf( painCave.errMsg,
               "GB::getLJParam could not convert GenericData to LJParam for\n"
               "\tatom type %s\n", atomType->getName().c_str());
      painCave.severity = OPENMD_ERROR;
      painCave.isFatal = 1;
      simError();          
    }
    
    return ljData->getData();
  }
  
  RealType GB::getLJEpsilon(AtomType* atomType) {    
    LJParam ljParam = getLJParam(atomType);
    return ljParam.epsilon;
  }
  RealType GB::getLJSigma(AtomType* atomType) {    
    LJParam ljParam = getLJParam(atomType);
    return ljParam.sigma;
  }
  
  void GB::initialize() {    
    
    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
    mu_ = fopts.getGayBerneMu();
    nu_ = fopts.getGayBerneNu();
    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
    ForceField::AtomTypeContainer::MapTypeIterator i;
    AtomType* at;

    // GB handles all of the GB-GB interactions as well as GB-LJ cross
    // interactions:

    for (at = atomTypes->beginType(i); at != NULL;
         at = atomTypes->nextType(i)) {
      
      if (at->isGayBerne() || at->isLennardJones())
        addType(at);
    }
   
    initialized_ = true;
  }
      
  void GB::addType(AtomType* atomType){
    // add it to the map:
    AtomTypeProperties atp = atomType->getATP();    

    pair<map<int,AtomType*>::iterator,bool> ret;    
    ret = GBMap.insert( pair<int, AtomType*>(atp.ident, atomType) );
    if (ret.second == false) {
      sprintf( painCave.errMsg,
               "GB already had a previous entry with ident %d\n",
               atp.ident);
      painCave.severity = OPENMD_INFO;
      painCave.isFatal = 0;
      simError();         
    }
    
    RealType d1, l1, eX1, eS1, eE1, dw1;
    
    if (atomType->isGayBerne()) {
      GayBerneParam gb1 = getGayBerneParam(atomType);
      d1 = gb1.GB_d;
      l1 = gb1.GB_l;
      eX1 = gb1.GB_eps_X;
      eS1 = gb1.GB_eps_S;
      eE1 = gb1.GB_eps_E;
      dw1 = gb1.GB_dw;
    } else if (atomType->isLennardJones()) {
      d1 = getLJSigma(atomType) / sqrt(2.0);
      l1 = d1;
      eX1 = getLJEpsilon(atomType);
      eS1 = eX1;
      eE1 = eX1;
      dw1 = 1.0;      
    } else {
      sprintf( painCave.errMsg,
               "GB::addType was passed an atomType (%s) that does not\n"
               "\tappear to be a Gay-Berne or Lennard-Jones atom.\n",
               atomType->getName().c_str());
      painCave.severity = OPENMD_ERROR;
      painCave.isFatal = 1;
      simError();
    }
      

    // Now, iterate over all known types and add to the mixing map:
    
    map<int, AtomType*>::iterator it;
    for( it = GBMap.begin(); it != GBMap.end(); ++it) {
      
      AtomType* atype2 = (*it).second;
      
      RealType d2, l2, eX2, eS2, eE2, dw2;
      
      if (atype2->isGayBerne()) {
        GayBerneParam gb2 = getGayBerneParam(atype2);
        d2 = gb2.GB_d;
        l2 = gb2.GB_l;
        eX2 = gb2.GB_eps_X;
        eS2 = gb2.GB_eps_S;
        eE2 = gb2.GB_eps_E;
        dw2 = gb2.GB_dw;
      } else if (atype2->isLennardJones()) {
        d2 = getLJSigma(atype2) / sqrt(2.0);
        l2 = d2;
        eX2 = getLJEpsilon(atype2);
        eS2 = eX2;
        eE2 = eX2;
        dw2 = 1.0;
      } 
                       
      GBInteractionData mixer1, mixer2;     
      
      //  Cleaver paper uses sqrt of squares to get sigma0 for
      //  mixed interactions.
            
      mixer1.sigma0 = sqrt(d1*d1 + d2*d2);
      mixer1.xa2 = (l1*l1 - d1*d1)/(l1*l1 + d2*d2);
      mixer1.xai2 = (l2*l2 - d2*d2)/(l2*l2 + d1*d1);
      mixer1.x2 = (l1*l1 - d1*d1) * (l2*l2 - d2*d2) /
        ((l2*l2 + d1*d1) * (l1*l1 + d2*d2));

      mixer2.sigma0 = mixer1.sigma0;
      // xa2 and xai2 for j-i pairs are reversed from the same i-j pairing.
      // Swapping the particles reverses the anisotropy parameters:
      mixer2.xa2 = mixer1.xai2;
      mixer2.xai2 = mixer1.xa2;
      mixer2.x2 = mixer1.x2;
 
      // assumed LB mixing rules for now:

      mixer1.dw = 0.5 * (dw1 + dw2);
      mixer1.eps0 = sqrt(eX1 * eX2);

      mixer2.dw = mixer1.dw;
      mixer2.eps0 = mixer1.eps0;

      RealType mi = RealType(1.0)/mu_;
      
      mixer1.xpap2  = (pow(eS1, mi) - pow(eE1, mi)) / (pow(eS1, mi) + pow(eE2, mi));
      mixer1.xpapi2 = (pow(eS2, mi) - pow(eE2, mi)) / (pow(eS2, mi) + pow(eE1, mi));
      mixer1.xp2    = (pow(eS1, mi) - pow(eE1, mi)) * (pow(eS2, mi) - pow(eE2, mi))  / 
        (pow(eS2, mi) + pow(eE1, mi)) / (pow(eS1, mi) + pow(eE2, mi)) ;

      // xpap2 and xpapi2 for j-i pairs are reversed from the same i-j pairing.
      // Swapping the particles reverses the anisotropy parameters:
      mixer2.xpap2 = mixer1.xpapi2;
      mixer2.xpapi2 = mixer1.xpap2;
      mixer2.xp2 = mixer1.xp2;

      // only add this pairing if at least one of the atoms is a Gay-Berne atom

      if (atomType->isGayBerne() || atype2->isGayBerne()) {

        pair<AtomType*, AtomType*> key1, key2;
        key1 = make_pair(atomType, atype2);
        key2 = make_pair(atype2, atomType);
        
        MixingMap[key1] = mixer1;
        if (key2 != key1) {
          MixingMap[key2] = mixer2;
        }
      }
    }      
  }
   
  void GB::calcForce(InteractionData &idat) {

    if (!initialized_) initialize();
    
    GBInteractionData mixer = MixingMap[idat.atypes];

    RealType sigma0 = mixer.sigma0;
    RealType dw     = mixer.dw;
    RealType eps0   = mixer.eps0;  
    RealType x2     = mixer.x2;    
    RealType xa2    = mixer.xa2;   
    RealType xai2   = mixer.xai2;  
    RealType xp2    = mixer.xp2;   
    RealType xpap2  = mixer.xpap2; 
    RealType xpapi2 = mixer.xpapi2;

    // cerr << "atypes = " << idat.atypes.first->getName() << " " << idat.atypes.second->getName() << "\n";
    // cerr << "sigma0 = " <<mixer.sigma0 <<"\n";
    // cerr << "dw     = " <<mixer.dw <<"\n";
    // cerr << "eps0   = " <<mixer.eps0 <<"\n";  
    // cerr << "x2     = " <<mixer.x2 <<"\n";    
    // cerr << "xa2    = " <<mixer.xa2 <<"\n";   
    // cerr << "xai2   = " <<mixer.xai2 <<"\n";  
    // cerr << "xp2    = " <<mixer.xp2 <<"\n";   
    // cerr << "xpap2  = " <<mixer.xpap2 <<"\n"; 
    // cerr << "xpapi2 = " <<mixer.xpapi2 <<"\n";

    Vector3d ul1 = idat.A1->getRow(2);
    Vector3d ul2 = idat.A2->getRow(2);

    // cerr << "ul1 = " <<ul1<<"\n";
    // cerr << "ul2 = " <<ul2<<"\n";

    RealType a, b, g;

    bool i_is_LJ = idat.atypes.first->isLennardJones();
    bool j_is_LJ = idat.atypes.second->isLennardJones();
    
    if (i_is_LJ) {
      a = 0.0;
      ul1 = V3Zero;
    } else {
      a = dot(*(idat.d), ul1);
    }

    if (j_is_LJ) {
      b = 0.0;
      ul2 = V3Zero;
    } else {
      b = dot(*(idat.d), ul2);
    }

    if (i_is_LJ || j_is_LJ) 
      g = 0.0;
    else
      g = dot(ul1, ul2);

    RealType au = a / *(idat.rij);
    RealType bu = b / *(idat.rij);
    
    RealType au2 = au * au;
    RealType bu2 = bu * bu;
    RealType g2 = g * g;

    RealType H  = (xa2 * au2 + xai2 * bu2 - 2.0*x2*au*bu*g)  / (1.0 - x2*g2);
    RealType Hp = (xpap2*au2 + xpapi2*bu2 - 2.0*xp2*au*bu*g) / (1.0 - xp2*g2);

    // cerr << "au2 = " << au2 << "\n";
    // cerr << "bu2 = " << bu2 << "\n";
    // cerr << "g2 = " << g2 << "\n";
    // cerr << "H = " << H << "\n";
    // cerr << "Hp = " << Hp << "\n";

    RealType sigma = sigma0 / sqrt(1.0 - H);
    RealType e1 = 1.0 / sqrt(1.0 - x2*g2);
    RealType e2 = 1.0 - Hp;
    RealType eps = eps0 * pow(e1,nu_) * pow(e2,mu_);
    RealType BigR = dw*sigma0 / (*(idat.rij) - sigma + dw*sigma0);
    
    RealType R3 = BigR*BigR*BigR;
    RealType R6 = R3*R3;
    RealType R7 = R6 * BigR;
    RealType R12 = R6*R6;
    RealType R13 = R6*R7;

    RealType U = *(idat.vdwMult) * 4.0 * eps * (R12 - R6);

    RealType s3 = sigma*sigma*sigma;
    RealType s03 = sigma0*sigma0*sigma0;

    // cerr << "vdwMult = " << *(idat.vdwMult) << "\n";
    // cerr << "eps = " << eps <<"\n";
    // cerr << "mu = " << mu_ << "\n";
    // cerr << "R12 = " << R12 << "\n";
    // cerr << "R6 = " << R6 << "\n";
    // cerr << "R13 = " << R13 << "\n";
    // cerr << "R7 = " << R7 << "\n";
    // cerr << "e2 = " << e2 << "\n";
    // cerr << "rij = " << *(idat.rij) << "\n";
    // cerr << "s3 = " << s3 << "\n";
    // cerr << "s03 = " << s03 << "\n";
    // cerr << "dw = " << dw << "\n";

    RealType pref1 = - *(idat.vdwMult) * 8.0 * eps * mu_ * (R12 - R6) / 
      (e2 * *(idat.rij));

    RealType pref2 = *(idat.vdwMult) * 8.0 * eps * s3 * (6.0*R13 - 3.0*R7) /
      (dw*  *(idat.rij) * s03);

    RealType dUdr = - (pref1 * Hp + pref2 * (sigma0 * sigma0 *  
                                             *(idat.rij) / s3 + H));
    
    RealType dUda = pref1 * (xpap2*au - xp2*bu*g) / (1.0 - xp2 * g2) 
      + pref2 * (xa2 * au - x2 *bu*g) / (1.0 - x2 * g2);
    
    RealType dUdb = pref1 * (xpapi2*bu - xp2*au*g) / (1.0 - xp2 * g2) 
      + pref2 * (xai2 * bu - x2 *au*g) / (1.0 - x2 * g2);

    RealType dUdg = 4.0 * eps * nu_ * (R12 - R6) * x2 * g / (1.0 - x2*g2)
      + 8.0 * eps * mu_ * (R12 - R6) * (xp2*au*bu - Hp*xp2*g) / 
      (1.0 - xp2 * g2) / e2 + 8.0 * eps * s3 * (3.0 * R7 - 6.0 * R13) * 
      (x2 * au * bu - H * x2 * g) / (1.0 - x2 * g2) / (dw * s03);

    // cerr << "pref = " << pref1 << " " << pref2 << "\n";
    // cerr << "dU = " << dUdr << " " << dUda <<" " << dUdb << " " << dUdg << "\n";

    Vector3d rhat = *(idat.d) / *(idat.rij);   
    Vector3d rxu1 = cross(*(idat.d), ul1);
    Vector3d rxu2 = cross(*(idat.d), ul2);
    Vector3d uxu = cross(ul1, ul2);

    (*(idat.pot))[VANDERWAALS_FAMILY] += U *  *(idat.sw);
    *(idat.f1) += (dUdr * rhat + dUda * ul1 + dUdb * ul2) * *(idat.sw);
    *(idat.t1) += (dUda * rxu1 - dUdg * uxu) * *(idat.sw);
    *(idat.t2) += (dUdb * rxu2 + dUdg * uxu) * *(idat.sw);
    *(idat.vpair) += U;

    // cerr << "f1 term = " <<  (dUdr * rhat + dUda * ul1 + dUdb * ul2) * *(idat.sw) << "\n";
    // cerr << "t1 term = " << (dUda * rxu1 - dUdg * uxu) * *(idat.sw) << "\n";
    // cerr << "t2 term = " << (dUdb * rxu2 + dUdg * uxu) * *(idat.sw) << "\n";
    // cerr << "vp term = " << U << "\n";

    return;

  }

  RealType GB::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
    if (!initialized_) initialize();   

    RealType cut = 0.0;

    if (atypes.first->isGayBerne()) {
      GayBerneParam gb1 = getGayBerneParam(atypes.first);
      RealType d1 = gb1.GB_d;
      RealType l1 = gb1.GB_l;
      // sigma is actually sqrt(2)*l  for prolate ellipsoids 
      cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d1, l1));
    } else if (atypes.first->isLennardJones()) {
      cut = max(cut, RealType(2.5) * getLJSigma(atypes.first));
    }

    if (atypes.second->isGayBerne()) {
      GayBerneParam gb2 = getGayBerneParam(atypes.second);
      RealType d2 = gb2.GB_d;
      RealType l2 = gb2.GB_l;
      cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d2, l2));
    } else if (atypes.second->isLennardJones()) {
      cut = max(cut, RealType(2.5) * getLJSigma(atypes.second));
    }
   
    return cut;
  }
}

Revision:	1688
Committed:	Wed Mar 14 17:56:01 2012 UTC (13 years, 4 months ago) by gezelter
File size:	17722 byte(s)
Log Message:	Bug fixes for GB. Now using strength parameter mixing ideas from Wu et al. [J. Chem. Phys. 135, 155104 (2011)]. This helps get the dissimilar particle mixing behavior to be the same whichever order the two particles come in. This does require that the force field file to specify explicitly the values for epsilon in the cross (X), side-by-side (S), and end-to-end (E) configurations.
#	User	Rev	Content
1	gezelter	1483	/*
2			* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9			* 1. Redistributions of source code must retain the above copyright
10			* notice, this list of conditions and the following disclaimer.
11			*
12			* 2. Redistributions in binary form must reproduce the above copyright
13			* notice, this list of conditions and the following disclaimer in the
14			* documentation and/or other materials provided with the
15			* distribution.
16			*
17			* This software is provided "AS IS," without a warranty of any
18			* kind. All express or implied conditions, representations and
19			* warranties, including any implied warranty of merchantability,
20			* fitness for a particular purpose or non-infringement, are hereby
21			* excluded. The University of Notre Dame and its licensors shall not
22			* be liable for any damages suffered by licensee as a result of
23			* using, modifying or distributing the software or its
24			* derivatives. In no event will the University of Notre Dame or its
25			* licensors be liable for any lost revenue, profit or data, or for
26			* direct, indirect, special, consequential, incidental or punitive
27			* damages, however caused and regardless of the theory of liability,
28			* arising out of the use of or inability to use software, even if the
29			* University of Notre Dame has been advised of the possibility of
30			* such damages.
31			*
32			* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33			* research, please cite the appropriate papers when you publish your
34			* work. Good starting points are:
35			*
36			* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37			* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38			* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	gezelter	1665	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	gezelter	1483	*/
42
43			#include <stdio.h>
44			#include <string.h>
45
46			#include <cmath>
47			#include "nonbonded/GB.hpp"
48			#include "utils/simError.h"
49
50			using namespace std;
51			namespace OpenMD {
52
53	gezelter	1688	/* GB is the Gay-Berne interaction for ellipsoidal particles. The original
54			* paper (for identical uniaxial particles) is:
55			* J. G. Gay and B. J. Berne, J. Chem. Phys., 74, 3316-3319 (1981).
56			* A more-general GB potential for dissimilar uniaxial particles:
57			* D. J. Cleaver, C. M. Care, M. P. Allen and M. P. Neal, Phys. Rev. E,
58			* 54, 559-567 (1996).
59			* Further parameterizations can be found in:
60			* A. P. J. Emerson, G. R. Luckhurst and S. G. Whatling, Mol. Phys.,
61			* 82, 113-124 (1994).
62			* And a nice force expression:
63			* G. R. Luckhurst and R. A. Stephens, Liq. Cryst. 8, 451-464 (1990).
64			* Even clearer force and torque expressions:
65			* P. A. Golubkov and P. Y. Ren, J. Chem. Phys., 125, 64103 (2006).
66			* New expressions for cross interactions of strength parameters:
67			* J. Wu, X. Zhen, H. Shen, G. Li, and P. Ren, J. Chem. Phys.,
68			* 135, 155104 (2011).
69			*
70			* In this version of the GB interaction, each uniaxial ellipsoidal type
71			* is described using a set of 6 parameters:
72			* d: range parameter for side-by-side (S) and cross (X) configurations
73			* l: range parameter for end-to-end (E) configuration
74			* epsilon_X: well-depth parameter for cross (X) configuration
75			* epsilon_S: well-depth parameter for side-by-side (S) configuration
76			* epsilon_E: well depth parameter for end-to-end (E) configuration
77			* dw: "softness" of the potential
78			*
79			* Additionally, there are two "universal" paramters to govern the overall
80			* importance of the purely orientational (nu) and the mixed
81			* orientational / translational (mu) parts of strength of the interactions.
82			* These parameters have default or "canonical" values, but may be changed
83			* as a force field option:
84			* nu_: purely orientational part : defaults to 1
85			* mu_: mixed orientational / translational part : defaults to 2
86			*/
87
88
89	gezelter	1502	GB::GB() : name_("GB"), initialized_(false), mu_(2.0), nu_(1.0), forceField_(NULL) {}
90	gezelter	1483
91			GayBerneParam GB::getGayBerneParam(AtomType* atomType) {
92
93			// Do sanity checking on the AtomType we were passed before
94			// building any data structures:
95			if (!atomType->isGayBerne()) {
96			sprintf( painCave.errMsg,
97			"GB::getGayBerneParam was passed an atomType (%s) that does\n"
98			"\tnot appear to be a Gay-Berne atom.\n",
99			atomType->getName().c_str());
100			painCave.severity = OPENMD_ERROR;
101			painCave.isFatal = 1;
102			simError();
103			}
104
105			DirectionalAtomType* daType = dynamic_cast<DirectionalAtomType*>(atomType);
106			GenericData* data = daType->getPropertyByName("GayBerne");
107			if (data == NULL) {
108			sprintf( painCave.errMsg, "GB::getGayBerneParam could not find\n"
109			"\tGay-Berne parameters for atomType %s.\n",
110			daType->getName().c_str());
111			painCave.severity = OPENMD_ERROR;
112			painCave.isFatal = 1;
113			simError();
114			}
115
116			GayBerneParamGenericData* gbData = dynamic_cast<GayBerneParamGenericData*>(data);
117			if (gbData == NULL) {
118			sprintf( painCave.errMsg,
119			"GB::getGayBerneParam could not convert GenericData to\n"
120			"\tGayBerneParamGenericData for atom type %s\n",
121			daType->getName().c_str());
122			painCave.severity = OPENMD_ERROR;
123			painCave.isFatal = 1;
124			simError();
125			}
126
127			return gbData->getData();
128			}
129
130	gezelter	1502	LJParam GB::getLJParam(AtomType* atomType) {
131
132			// Do sanity checking on the AtomType we were passed before
133			// building any data structures:
134			if (!atomType->isLennardJones()) {
135			sprintf( painCave.errMsg,
136			"GB::getLJParam was passed an atomType (%s) that does not\n"
137			"\tappear to be a Lennard-Jones atom.\n",
138			atomType->getName().c_str());
139			painCave.severity = OPENMD_ERROR;
140			painCave.isFatal = 1;
141			simError();
142			}
143
144			GenericData* data = atomType->getPropertyByName("LennardJones");
145			if (data == NULL) {
146			sprintf( painCave.errMsg, "GB::getLJParam could not find Lennard-Jones\n"
147			"\tparameters for atomType %s.\n", atomType->getName().c_str());
148			painCave.severity = OPENMD_ERROR;
149			painCave.isFatal = 1;
150			simError();
151			}
152
153			LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
154			if (ljData == NULL) {
155			sprintf( painCave.errMsg,
156			"GB::getLJParam could not convert GenericData to LJParam for\n"
157			"\tatom type %s\n", atomType->getName().c_str());
158			painCave.severity = OPENMD_ERROR;
159			painCave.isFatal = 1;
160			simError();
161			}
162
163			return ljData->getData();
164			}
165
166			RealType GB::getLJEpsilon(AtomType* atomType) {
167			LJParam ljParam = getLJParam(atomType);
168			return ljParam.epsilon;
169			}
170			RealType GB::getLJSigma(AtomType* atomType) {
171			LJParam ljParam = getLJParam(atomType);
172			return ljParam.sigma;
173			}
174
175	gezelter	1483	void GB::initialize() {
176	gezelter	1502
177	gezelter	1485	ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
178			mu_ = fopts.getGayBerneMu();
179			nu_ = fopts.getGayBerneNu();
180	gezelter	1483	ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
181			ForceField::AtomTypeContainer::MapTypeIterator i;
182			AtomType* at;
183
184			// GB handles all of the GB-GB interactions as well as GB-LJ cross
185			// interactions:
186
187			for (at = atomTypes->beginType(i); at != NULL;
188			at = atomTypes->nextType(i)) {
189
190			if (at->isGayBerne() \|\| at->isLennardJones())
191			addType(at);
192			}
193
194			initialized_ = true;
195			}
196
197			void GB::addType(AtomType* atomType){
198			// add it to the map:
199			AtomTypeProperties atp = atomType->getATP();
200
201			pair<map<int,AtomType*>::iterator,bool> ret;
202			ret = GBMap.insert( pair<int, AtomType*>(atp.ident, atomType) );
203			if (ret.second == false) {
204			sprintf( painCave.errMsg,
205			"GB already had a previous entry with ident %d\n",
206			atp.ident);
207			painCave.severity = OPENMD_INFO;
208			painCave.isFatal = 0;
209			simError();
210			}
211
212	gezelter	1688	RealType d1, l1, eX1, eS1, eE1, dw1;
213	gezelter	1483
214			if (atomType->isGayBerne()) {
215			GayBerneParam gb1 = getGayBerneParam(atomType);
216			d1 = gb1.GB_d;
217			l1 = gb1.GB_l;
218	gezelter	1688	eX1 = gb1.GB_eps_X;
219			eS1 = gb1.GB_eps_S;
220			eE1 = gb1.GB_eps_E;
221	gezelter	1483	dw1 = gb1.GB_dw;
222			} else if (atomType->isLennardJones()) {
223	gezelter	1502	d1 = getLJSigma(atomType) / sqrt(2.0);
224	gezelter	1483	l1 = d1;
225	gezelter	1688	eX1 = getLJEpsilon(atomType);
226			eS1 = eX1;
227			eE1 = eX1;
228	gezelter	1483	dw1 = 1.0;
229			} else {
230			sprintf( painCave.errMsg,
231			"GB::addType was passed an atomType (%s) that does not\n"
232			"\tappear to be a Gay-Berne or Lennard-Jones atom.\n",
233			atomType->getName().c_str());
234			painCave.severity = OPENMD_ERROR;
235			painCave.isFatal = 1;
236			simError();
237			}
238
239
240			// Now, iterate over all known types and add to the mixing map:
241
242			map<int, AtomType*>::iterator it;
243			for( it = GBMap.begin(); it != GBMap.end(); ++it) {
244
245			AtomType* atype2 = (*it).second;
246
247	gezelter	1688	RealType d2, l2, eX2, eS2, eE2, dw2;
248	gezelter	1483
249			if (atype2->isGayBerne()) {
250			GayBerneParam gb2 = getGayBerneParam(atype2);
251			d2 = gb2.GB_d;
252			l2 = gb2.GB_l;
253	gezelter	1688	eX2 = gb2.GB_eps_X;
254			eS2 = gb2.GB_eps_S;
255			eE2 = gb2.GB_eps_E;
256	gezelter	1483	dw2 = gb2.GB_dw;
257			} else if (atype2->isLennardJones()) {
258	gezelter	1502	d2 = getLJSigma(atype2) / sqrt(2.0);
259	gezelter	1483	l2 = d2;
260	gezelter	1688	eX2 = getLJEpsilon(atype2);
261			eS2 = eX2;
262			eE2 = eX2;
263	gezelter	1483	dw2 = 1.0;
264			}
265
266	gezelter	1674	GBInteractionData mixer1, mixer2;
267	gezelter	1483
268			// Cleaver paper uses sqrt of squares to get sigma0 for
269			// mixed interactions.
270
271	gezelter	1674	mixer1.sigma0 = sqrt(d1d1 + d2d2);
272			mixer1.xa2 = (l1l1 - d1d1)/(l1l1 + d2d2);
273			mixer1.xai2 = (l2l2 - d2d2)/(l2l2 + d1d1);
274			mixer1.x2 = (l1l1 - d1d1) * (l2l2 - d2d2) /
275	gezelter	1483	((l2l2 + d1d1) * (l1l1 + d2d2));
276	gezelter	1674
277			mixer2.sigma0 = mixer1.sigma0;
278			// xa2 and xai2 for j-i pairs are reversed from the same i-j pairing.
279			// Swapping the particles reverses the anisotropy parameters:
280			mixer2.xa2 = mixer1.xai2;
281			mixer2.xai2 = mixer1.xa2;
282			mixer2.x2 = mixer1.x2;
283	gezelter	1483
284			// assumed LB mixing rules for now:
285
286	gezelter	1674	mixer1.dw = 0.5 * (dw1 + dw2);
287	gezelter	1688	mixer1.eps0 = sqrt(eX1 * eX2);
288	gezelter	1674
289			mixer2.dw = mixer1.dw;
290			mixer2.eps0 = mixer1.eps0;
291	gezelter	1688
292			RealType mi = RealType(1.0)/mu_;
293	gezelter	1483
294	gezelter	1688	mixer1.xpap2 = (pow(eS1, mi) - pow(eE1, mi)) / (pow(eS1, mi) + pow(eE2, mi));
295			mixer1.xpapi2 = (pow(eS2, mi) - pow(eE2, mi)) / (pow(eS2, mi) + pow(eE1, mi));
296			mixer1.xp2 = (pow(eS1, mi) - pow(eE1, mi)) * (pow(eS2, mi) - pow(eE2, mi)) /
297			(pow(eS2, mi) + pow(eE1, mi)) / (pow(eS1, mi) + pow(eE2, mi)) ;
298	gezelter	1483
299	gezelter	1688	// xpap2 and xpapi2 for j-i pairs are reversed from the same i-j pairing.
300			// Swapping the particles reverses the anisotropy parameters:
301			mixer2.xpap2 = mixer1.xpapi2;
302			mixer2.xpapi2 = mixer1.xpap2;
303	gezelter	1674	mixer2.xp2 = mixer1.xp2;
304	gezelter	1587
305	gezelter	1483	// only add this pairing if at least one of the atoms is a Gay-Berne atom
306
307			if (atomType->isGayBerne() \|\| atype2->isGayBerne()) {
308
309			pair<AtomType, AtomType> key1, key2;
310			key1 = make_pair(atomType, atype2);
311			key2 = make_pair(atype2, atomType);
312
313	gezelter	1674	MixingMap[key1] = mixer1;
314	gezelter	1483	if (key2 != key1) {
315	gezelter	1674	MixingMap[key2] = mixer2;
316	gezelter	1483	}
317			}
318			}
319			}
320	gezelter	1502
321	gezelter	1536	void GB::calcForce(InteractionData &idat) {
322	gezelter	1483
323			if (!initialized_) initialize();
324
325	gezelter	1571	GBInteractionData mixer = MixingMap[idat.atypes];
326	gezelter	1483
327			RealType sigma0 = mixer.sigma0;
328			RealType dw = mixer.dw;
329			RealType eps0 = mixer.eps0;
330			RealType x2 = mixer.x2;
331			RealType xa2 = mixer.xa2;
332			RealType xai2 = mixer.xai2;
333			RealType xp2 = mixer.xp2;
334			RealType xpap2 = mixer.xpap2;
335			RealType xpapi2 = mixer.xpapi2;
336
337	gezelter	1688	// cerr << "atypes = " << idat.atypes.first->getName() << " " << idat.atypes.second->getName() << "\n";
338			// cerr << "sigma0 = " <<mixer.sigma0 <<"\n";
339			// cerr << "dw = " <<mixer.dw <<"\n";
340			// cerr << "eps0 = " <<mixer.eps0 <<"\n";
341			// cerr << "x2 = " <<mixer.x2 <<"\n";
342			// cerr << "xa2 = " <<mixer.xa2 <<"\n";
343			// cerr << "xai2 = " <<mixer.xai2 <<"\n";
344			// cerr << "xp2 = " <<mixer.xp2 <<"\n";
345			// cerr << "xpap2 = " <<mixer.xpap2 <<"\n";
346			// cerr << "xpapi2 = " <<mixer.xpapi2 <<"\n";
347
348	gezelter	1554	Vector3d ul1 = idat.A1->getRow(2);
349			Vector3d ul2 = idat.A2->getRow(2);
350	gezelter	1483
351	gezelter	1688	// cerr << "ul1 = " <<ul1<<"\n";
352			// cerr << "ul2 = " <<ul2<<"\n";
353
354	gezelter	1483	RealType a, b, g;
355
356	gezelter	1571	bool i_is_LJ = idat.atypes.first->isLennardJones();
357			bool j_is_LJ = idat.atypes.second->isLennardJones();
358	gezelter	1554
359	gezelter	1483	if (i_is_LJ) {
360			a = 0.0;
361			ul1 = V3Zero;
362			} else {
363	gezelter	1554	a = dot(*(idat.d), ul1);
364	gezelter	1483	}
365
366			if (j_is_LJ) {
367			b = 0.0;
368			ul2 = V3Zero;
369			} else {
370	gezelter	1554	b = dot(*(idat.d), ul2);
371	gezelter	1483	}
372
373			if (i_is_LJ \|\| j_is_LJ)
374			g = 0.0;
375			else
376			g = dot(ul1, ul2);
377
378	gezelter	1554	RealType au = a / *(idat.rij);
379			RealType bu = b / *(idat.rij);
380	gezelter	1483
381			RealType au2 = au * au;
382			RealType bu2 = bu * bu;
383			RealType g2 = g * g;
384	gezelter	1688
385	gezelter	1483	RealType H = (xa2 * au2 + xai2 * bu2 - 2.0x2aubug) / (1.0 - x2*g2);
386			RealType Hp = (xpap2au2 + xpapi2bu2 - 2.0xp2aubug) / (1.0 - xp2*g2);
387
388	gezelter	1688	// cerr << "au2 = " << au2 << "\n";
389			// cerr << "bu2 = " << bu2 << "\n";
390			// cerr << "g2 = " << g2 << "\n";
391			// cerr << "H = " << H << "\n";
392			// cerr << "Hp = " << Hp << "\n";
393
394	gezelter	1483	RealType sigma = sigma0 / sqrt(1.0 - H);
395			RealType e1 = 1.0 / sqrt(1.0 - x2*g2);
396			RealType e2 = 1.0 - Hp;
397			RealType eps = eps0 * pow(e1,nu_) * pow(e2,mu_);
398	gezelter	1554	RealType BigR = dwsigma0 / ((idat.rij) - sigma + dw*sigma0);
399	gezelter	1483
400			RealType R3 = BigRBigRBigR;
401			RealType R6 = R3*R3;
402			RealType R7 = R6 * BigR;
403			RealType R12 = R6*R6;
404			RealType R13 = R6*R7;
405
406	gezelter	1554	RealType U = (idat.vdwMult) 4.0 * eps * (R12 - R6);
407	gezelter	1483
408			RealType s3 = sigmasigmasigma;
409			RealType s03 = sigma0sigma0sigma0;
410
411	gezelter	1688	// cerr << "vdwMult = " << *(idat.vdwMult) << "\n";
412			// cerr << "eps = " << eps <<"\n";
413			// cerr << "mu = " << mu_ << "\n";
414			// cerr << "R12 = " << R12 << "\n";
415			// cerr << "R6 = " << R6 << "\n";
416			// cerr << "R13 = " << R13 << "\n";
417			// cerr << "R7 = " << R7 << "\n";
418			// cerr << "e2 = " << e2 << "\n";
419			// cerr << "rij = " << *(idat.rij) << "\n";
420			// cerr << "s3 = " << s3 << "\n";
421			// cerr << "s03 = " << s03 << "\n";
422			// cerr << "dw = " << dw << "\n";
423
424	gezelter	1554	RealType pref1 = - (idat.vdwMult) 8.0 * eps * mu_ * (R12 - R6) /
425			(e2 * *(idat.rij));
426	gezelter	1483
427	gezelter	1554	RealType pref2 = (idat.vdwMult) 8.0 * eps * s3 * (6.0R13 - 3.0R7) /
428			(dw* (idat.rij) s03);
429	gezelter	1483
430	gezelter	1554	RealType dUdr = - (pref1 * Hp + pref2 * (sigma0 * sigma0 *
431			*(idat.rij) / s3 + H));
432	gezelter	1483
433			RealType dUda = pref1 * (xpap2au - xp2bug) / (1.0 - xp2 g2)
434			+ pref2 * (xa2 * au - x2 bug) / (1.0 - x2 * g2);
435
436			RealType dUdb = pref1 * (xpapi2bu - xp2aug) / (1.0 - xp2 g2)
437			+ pref2 * (xai2 * bu - x2 aug) / (1.0 - x2 * g2);
438
439			RealType dUdg = 4.0 * eps * nu_ * (R12 - R6) * x2 * g / (1.0 - x2*g2)
440			+ 8.0 * eps * mu_ * (R12 - R6) * (xp2aubu - Hpxp2g) /
441			(1.0 - xp2 * g2) / e2 + 8.0 * eps * s3 * (3.0 * R7 - 6.0 * R13) *
442			(x2 * au * bu - H * x2 * g) / (1.0 - x2 * g2) / (dw * s03);
443
444	gezelter	1688	// cerr << "pref = " << pref1 << " " << pref2 << "\n";
445			// cerr << "dU = " << dUdr << " " << dUda <<" " << dUdb << " " << dUdg << "\n";
446
447	gezelter	1554	Vector3d rhat = (idat.d) / (idat.rij);
448			Vector3d rxu1 = cross(*(idat.d), ul1);
449			Vector3d rxu2 = cross(*(idat.d), ul2);
450	gezelter	1483	Vector3d uxu = cross(ul1, ul2);
451	gezelter	1587
452	gezelter	1582	((idat.pot))[VANDERWAALS_FAMILY] += U *(idat.sw);
453	gezelter	1686	(idat.f1) += (dUdr rhat + dUda * ul1 + dUdb * ul2) * *(idat.sw);
454			(idat.t1) += (dUda rxu1 - dUdg * uxu) * *(idat.sw);
455			(idat.t2) += (dUdb rxu2 + dUdg * uxu) * *(idat.sw);
456			*(idat.vpair) += U;
457	gezelter	1483
458	gezelter	1688	// cerr << "f1 term = " << (dUdr * rhat + dUda * ul1 + dUdb * ul2) * *(idat.sw) << "\n";
459			// cerr << "t1 term = " << (dUda * rxu1 - dUdg * uxu) * *(idat.sw) << "\n";
460			// cerr << "t2 term = " << (dUdb * rxu2 + dUdg * uxu) * *(idat.sw) << "\n";
461			// cerr << "vp term = " << U << "\n";
462
463	gezelter	1483	return;
464
465			}
466	gezelter	1505
467	gezelter	1545	RealType GB::getSuggestedCutoffRadius(pair<AtomType, AtomType> atypes) {
468	gezelter	1505	if (!initialized_) initialize();
469
470			RealType cut = 0.0;
471
472	gezelter	1545	if (atypes.first->isGayBerne()) {
473			GayBerneParam gb1 = getGayBerneParam(atypes.first);
474	gezelter	1505	RealType d1 = gb1.GB_d;
475			RealType l1 = gb1.GB_l;
476			// sigma is actually sqrt(2)*l for prolate ellipsoids
477	gezelter	1668	cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d1, l1));
478	gezelter	1545	} else if (atypes.first->isLennardJones()) {
479	gezelter	1668	cut = max(cut, RealType(2.5) * getLJSigma(atypes.first));
480	gezelter	1505	}
481
482	gezelter	1545	if (atypes.second->isGayBerne()) {
483			GayBerneParam gb2 = getGayBerneParam(atypes.second);
484	gezelter	1505	RealType d2 = gb2.GB_d;
485			RealType l2 = gb2.GB_l;
486	gezelter	1668	cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d2, l2));
487	gezelter	1545	} else if (atypes.second->isLennardJones()) {
488	gezelter	1668	cut = max(cut, RealType(2.5) * getLJSigma(atypes.second));
489	gezelter	1505	}
490
491			return cut;
492			}
493	gezelter	1483	}
494