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.
#	Content
1	/*
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	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40	* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	*/
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	/* 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	GB::GB() : name_("GB"), initialized_(false), mu_(2.0), nu_(1.0), forceField_(NULL) {}
90
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	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	void GB::initialize() {
176
177	ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
178	mu_ = fopts.getGayBerneMu();
179	nu_ = fopts.getGayBerneNu();
180	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	RealType d1, l1, eX1, eS1, eE1, dw1;
213
214	if (atomType->isGayBerne()) {
215	GayBerneParam gb1 = getGayBerneParam(atomType);
216	d1 = gb1.GB_d;
217	l1 = gb1.GB_l;
218	eX1 = gb1.GB_eps_X;
219	eS1 = gb1.GB_eps_S;
220	eE1 = gb1.GB_eps_E;
221	dw1 = gb1.GB_dw;
222	} else if (atomType->isLennardJones()) {
223	d1 = getLJSigma(atomType) / sqrt(2.0);
224	l1 = d1;
225	eX1 = getLJEpsilon(atomType);
226	eS1 = eX1;
227	eE1 = eX1;
228	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	RealType d2, l2, eX2, eS2, eE2, dw2;
248
249	if (atype2->isGayBerne()) {
250	GayBerneParam gb2 = getGayBerneParam(atype2);
251	d2 = gb2.GB_d;
252	l2 = gb2.GB_l;
253	eX2 = gb2.GB_eps_X;
254	eS2 = gb2.GB_eps_S;
255	eE2 = gb2.GB_eps_E;
256	dw2 = gb2.GB_dw;
257	} else if (atype2->isLennardJones()) {
258	d2 = getLJSigma(atype2) / sqrt(2.0);
259	l2 = d2;
260	eX2 = getLJEpsilon(atype2);
261	eS2 = eX2;
262	eE2 = eX2;
263	dw2 = 1.0;
264	}
265
266	GBInteractionData mixer1, mixer2;
267
268	// Cleaver paper uses sqrt of squares to get sigma0 for
269	// mixed interactions.
270
271	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	((l2l2 + d1d1) * (l1l1 + d2d2));
276
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
284	// assumed LB mixing rules for now:
285
286	mixer1.dw = 0.5 * (dw1 + dw2);
287	mixer1.eps0 = sqrt(eX1 * eX2);
288
289	mixer2.dw = mixer1.dw;
290	mixer2.eps0 = mixer1.eps0;
291
292	RealType mi = RealType(1.0)/mu_;
293
294	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
299	// 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	mixer2.xp2 = mixer1.xp2;
304
305	// 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	MixingMap[key1] = mixer1;
314	if (key2 != key1) {
315	MixingMap[key2] = mixer2;
316	}
317	}
318	}
319	}
320
321	void GB::calcForce(InteractionData &idat) {
322
323	if (!initialized_) initialize();
324
325	GBInteractionData mixer = MixingMap[idat.atypes];
326
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	// 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	Vector3d ul1 = idat.A1->getRow(2);
349	Vector3d ul2 = idat.A2->getRow(2);
350
351	// cerr << "ul1 = " <<ul1<<"\n";
352	// cerr << "ul2 = " <<ul2<<"\n";
353
354	RealType a, b, g;
355
356	bool i_is_LJ = idat.atypes.first->isLennardJones();
357	bool j_is_LJ = idat.atypes.second->isLennardJones();
358
359	if (i_is_LJ) {
360	a = 0.0;
361	ul1 = V3Zero;
362	} else {
363	a = dot(*(idat.d), ul1);
364	}
365
366	if (j_is_LJ) {
367	b = 0.0;
368	ul2 = V3Zero;
369	} else {
370	b = dot(*(idat.d), ul2);
371	}
372
373	if (i_is_LJ \|\| j_is_LJ)
374	g = 0.0;
375	else
376	g = dot(ul1, ul2);
377
378	RealType au = a / *(idat.rij);
379	RealType bu = b / *(idat.rij);
380
381	RealType au2 = au * au;
382	RealType bu2 = bu * bu;
383	RealType g2 = g * g;
384
385	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	// 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	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	RealType BigR = dwsigma0 / ((idat.rij) - sigma + dw*sigma0);
399
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	RealType U = (idat.vdwMult) 4.0 * eps * (R12 - R6);
407
408	RealType s3 = sigmasigmasigma;
409	RealType s03 = sigma0sigma0sigma0;
410
411	// 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	RealType pref1 = - (idat.vdwMult) 8.0 * eps * mu_ * (R12 - R6) /
425	(e2 * *(idat.rij));
426
427	RealType pref2 = (idat.vdwMult) 8.0 * eps * s3 * (6.0R13 - 3.0R7) /
428	(dw* (idat.rij) s03);
429
430	RealType dUdr = - (pref1 * Hp + pref2 * (sigma0 * sigma0 *
431	*(idat.rij) / s3 + H));
432
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	// cerr << "pref = " << pref1 << " " << pref2 << "\n";
445	// cerr << "dU = " << dUdr << " " << dUda <<" " << dUdb << " " << dUdg << "\n";
446
447	Vector3d rhat = (idat.d) / (idat.rij);
448	Vector3d rxu1 = cross(*(idat.d), ul1);
449	Vector3d rxu2 = cross(*(idat.d), ul2);
450	Vector3d uxu = cross(ul1, ul2);
451
452	((idat.pot))[VANDERWAALS_FAMILY] += U *(idat.sw);
453	(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
458	// 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	return;
464
465	}
466
467	RealType GB::getSuggestedCutoffRadius(pair<AtomType, AtomType> atypes) {
468	if (!initialized_) initialize();
469
470	RealType cut = 0.0;
471
472	if (atypes.first->isGayBerne()) {
473	GayBerneParam gb1 = getGayBerneParam(atypes.first);
474	RealType d1 = gb1.GB_d;
475	RealType l1 = gb1.GB_l;
476	// sigma is actually sqrt(2)*l for prolate ellipsoids
477	cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d1, l1));
478	} else if (atypes.first->isLennardJones()) {
479	cut = max(cut, RealType(2.5) * getLJSigma(atypes.first));
480	}
481
482	if (atypes.second->isGayBerne()) {
483	GayBerneParam gb2 = getGayBerneParam(atypes.second);
484	RealType d2 = gb2.GB_d;
485	RealType l2 = gb2.GB_l;
486	cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d2, l2));
487	} else if (atypes.second->isLennardJones()) {
488	cut = max(cut, RealType(2.5) * getLJSigma(atypes.second));
489	}
490
491	return cut;
492	}
493	}
494