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"
#include "types/LennardJonesAdapter.hpp"
#include "types/GayBerneAdapter.hpp"

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) {}
    
  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)) {
      
      LennardJonesAdapter lja = LennardJonesAdapter(at);
      GayBerneAdapter gba = GayBerneAdapter(at);

      if (gba.isGayBerne() || lja.isLennardJones())
        addType(at);
    }
   
    initialized_ = true;
  }
      
  void GB::addType(AtomType* atomType){
    // add it to the map:

    pair<map<int,AtomType*>::iterator,bool> ret;    
    ret = GBMap.insert( pair<int, AtomType*>(atomType->getIdent(), atomType) );
    if (ret.second == false) {
      sprintf( painCave.errMsg,
               "GB already had a previous entry with ident %d\n",
               atomType->getIdent() );
      painCave.severity = OPENMD_INFO;
      painCave.isFatal = 0;
      simError();         
    }
    
    RealType d1, l1, eX1, eS1, eE1, dw1;
        
    LennardJonesAdapter lja1 = LennardJonesAdapter(atomType);
    GayBerneAdapter gba1 = GayBerneAdapter(atomType);
    if (gba1.isGayBerne()) {
      d1 = gba1.getD();
      l1 = gba1.getL();
      eX1 = gba1.getEpsX();
      eS1 = gba1.getEpsS();
      eE1 = gba1.getEpsE();
      dw1 = gba1.getDw();
    } else if (lja1.isLennardJones()) {
      d1 = lja1.getSigma() / sqrt(2.0);
      l1 = d1;
      eX1 = lja1.getEpsilon();
      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;
      LennardJonesAdapter lja2 = LennardJonesAdapter(atype2);
      GayBerneAdapter gba2 = GayBerneAdapter(atype2);
      RealType d2, l2, eX2, eS2, eE2, dw2;
      
      if (gba2.isGayBerne()) {
        d2 = gba2.getD();
        l2 = gba2.getL();
        eX2 = gba2.getEpsX();
        eS2 = gba2.getEpsS();
        eE2 = gba2.getEpsE();
        dw2 = gba2.getDw();
      } else if (lja2.isLennardJones()) {
        d2 = lja2.getSigma() / sqrt(2.0);
        l2 = d2;
        eX2 = lja2.getEpsilon();
        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 (gba1.isGayBerne() || gba2.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;

    // This is not good.  We should store this in the mixing map, and not
    // query atom types in calc force.
    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;

    LennardJonesAdapter lja1 = LennardJonesAdapter(atypes.first);
    GayBerneAdapter gba1 = GayBerneAdapter(atypes.first);
    LennardJonesAdapter lja2 = LennardJonesAdapter(atypes.second);
    GayBerneAdapter gba2 = GayBerneAdapter(atypes.second);

    if (gba1.isGayBerne()) {
      RealType d1 = gba1.getD();
      RealType l1 = gba1.getL();
      // sigma is actually sqrt(2)*l  for prolate ellipsoids 
      cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d1, l1));
    } else if (lja1.isLennardJones()) {
      cut = max(cut, RealType(2.5) * lja1.getSigma());
    }

    if (gba2.isGayBerne()) {
      RealType d2 = gba2.getD();
      RealType l2 = gba2.getL();
      cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d2, l2));
    } else if (lja2.isLennardJones()) {
      cut = max(cut, RealType(2.5) * lja2.getSigma());
    }
   
    return cut;
  }
}

Revision:	1710
Committed:	Fri May 18 21:44:02 2012 UTC (13 years, 2 months ago) by gezelter
File size:	15246 byte(s)
Log Message:	Added an adapter layer between the AtomType and the rest of the code to handle the bolt-on capabilities of new types. Fixed a long-standing bug in how storageLayout was being set to the maximum possible value. Started to add infrastructure for Polarizable and fluc-Q calculations.
#	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	gezelter	1710	#include "types/LennardJonesAdapter.hpp"
50			#include "types/GayBerneAdapter.hpp"
51	gezelter	1483
52			using namespace std;
53			namespace OpenMD {
54
55	gezelter	1688	/* GB is the Gay-Berne interaction for ellipsoidal particles. The original
56			* paper (for identical uniaxial particles) is:
57			* J. G. Gay and B. J. Berne, J. Chem. Phys., 74, 3316-3319 (1981).
58			* A more-general GB potential for dissimilar uniaxial particles:
59			* D. J. Cleaver, C. M. Care, M. P. Allen and M. P. Neal, Phys. Rev. E,
60			* 54, 559-567 (1996).
61			* Further parameterizations can be found in:
62			* A. P. J. Emerson, G. R. Luckhurst and S. G. Whatling, Mol. Phys.,
63			* 82, 113-124 (1994).
64			* And a nice force expression:
65			* G. R. Luckhurst and R. A. Stephens, Liq. Cryst. 8, 451-464 (1990).
66			* Even clearer force and torque expressions:
67			* P. A. Golubkov and P. Y. Ren, J. Chem. Phys., 125, 64103 (2006).
68			* New expressions for cross interactions of strength parameters:
69			* J. Wu, X. Zhen, H. Shen, G. Li, and P. Ren, J. Chem. Phys.,
70			* 135, 155104 (2011).
71			*
72			* In this version of the GB interaction, each uniaxial ellipsoidal type
73			* is described using a set of 6 parameters:
74			* d: range parameter for side-by-side (S) and cross (X) configurations
75			* l: range parameter for end-to-end (E) configuration
76			* epsilon_X: well-depth parameter for cross (X) configuration
77			* epsilon_S: well-depth parameter for side-by-side (S) configuration
78			* epsilon_E: well depth parameter for end-to-end (E) configuration
79			* dw: "softness" of the potential
80			*
81			* Additionally, there are two "universal" paramters to govern the overall
82			* importance of the purely orientational (nu) and the mixed
83			* orientational / translational (mu) parts of strength of the interactions.
84			* These parameters have default or "canonical" values, but may be changed
85			* as a force field option:
86			* nu_: purely orientational part : defaults to 1
87			* mu_: mixed orientational / translational part : defaults to 2
88			*/
89
90
91	gezelter	1502	GB::GB() : name_("GB"), initialized_(false), mu_(2.0), nu_(1.0), forceField_(NULL) {}
92	gezelter	1483
93			void GB::initialize() {
94	gezelter	1502
95	gezelter	1485	ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
96			mu_ = fopts.getGayBerneMu();
97			nu_ = fopts.getGayBerneNu();
98	gezelter	1483	ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
99			ForceField::AtomTypeContainer::MapTypeIterator i;
100			AtomType* at;
101
102			// GB handles all of the GB-GB interactions as well as GB-LJ cross
103			// interactions:
104
105			for (at = atomTypes->beginType(i); at != NULL;
106			at = atomTypes->nextType(i)) {
107
108	gezelter	1710	LennardJonesAdapter lja = LennardJonesAdapter(at);
109			GayBerneAdapter gba = GayBerneAdapter(at);
110
111			if (gba.isGayBerne() \|\| lja.isLennardJones())
112	gezelter	1483	addType(at);
113			}
114
115			initialized_ = true;
116			}
117
118			void GB::addType(AtomType* atomType){
119			// add it to the map:
120
121			pair<map<int,AtomType*>::iterator,bool> ret;
122	gezelter	1710	ret = GBMap.insert( pair<int, AtomType*>(atomType->getIdent(), atomType) );
123	gezelter	1483	if (ret.second == false) {
124			sprintf( painCave.errMsg,
125			"GB already had a previous entry with ident %d\n",
126	gezelter	1710	atomType->getIdent() );
127	gezelter	1483	painCave.severity = OPENMD_INFO;
128			painCave.isFatal = 0;
129			simError();
130			}
131
132	gezelter	1688	RealType d1, l1, eX1, eS1, eE1, dw1;
133	gezelter	1710
134			LennardJonesAdapter lja1 = LennardJonesAdapter(atomType);
135			GayBerneAdapter gba1 = GayBerneAdapter(atomType);
136			if (gba1.isGayBerne()) {
137			d1 = gba1.getD();
138			l1 = gba1.getL();
139			eX1 = gba1.getEpsX();
140			eS1 = gba1.getEpsS();
141			eE1 = gba1.getEpsE();
142			dw1 = gba1.getDw();
143			} else if (lja1.isLennardJones()) {
144			d1 = lja1.getSigma() / sqrt(2.0);
145	gezelter	1483	l1 = d1;
146	gezelter	1710	eX1 = lja1.getEpsilon();
147	gezelter	1688	eS1 = eX1;
148			eE1 = eX1;
149	gezelter	1483	dw1 = 1.0;
150			} else {
151			sprintf( painCave.errMsg,
152			"GB::addType was passed an atomType (%s) that does not\n"
153			"\tappear to be a Gay-Berne or Lennard-Jones atom.\n",
154			atomType->getName().c_str());
155			painCave.severity = OPENMD_ERROR;
156			painCave.isFatal = 1;
157			simError();
158			}
159
160
161			// Now, iterate over all known types and add to the mixing map:
162
163			map<int, AtomType*>::iterator it;
164			for( it = GBMap.begin(); it != GBMap.end(); ++it) {
165
166			AtomType* atype2 = (*it).second;
167	gezelter	1710	LennardJonesAdapter lja2 = LennardJonesAdapter(atype2);
168			GayBerneAdapter gba2 = GayBerneAdapter(atype2);
169	gezelter	1688	RealType d2, l2, eX2, eS2, eE2, dw2;
170	gezelter	1483
171	gezelter	1710	if (gba2.isGayBerne()) {
172			d2 = gba2.getD();
173			l2 = gba2.getL();
174			eX2 = gba2.getEpsX();
175			eS2 = gba2.getEpsS();
176			eE2 = gba2.getEpsE();
177			dw2 = gba2.getDw();
178			} else if (lja2.isLennardJones()) {
179			d2 = lja2.getSigma() / sqrt(2.0);
180	gezelter	1483	l2 = d2;
181	gezelter	1710	eX2 = lja2.getEpsilon();
182	gezelter	1688	eS2 = eX2;
183			eE2 = eX2;
184	gezelter	1483	dw2 = 1.0;
185			}
186
187	gezelter	1674	GBInteractionData mixer1, mixer2;
188	gezelter	1483
189			// Cleaver paper uses sqrt of squares to get sigma0 for
190			// mixed interactions.
191
192	gezelter	1674	mixer1.sigma0 = sqrt(d1d1 + d2d2);
193			mixer1.xa2 = (l1l1 - d1d1)/(l1l1 + d2d2);
194			mixer1.xai2 = (l2l2 - d2d2)/(l2l2 + d1d1);
195			mixer1.x2 = (l1l1 - d1d1) * (l2l2 - d2d2) /
196	gezelter	1483	((l2l2 + d1d1) * (l1l1 + d2d2));
197	gezelter	1674
198			mixer2.sigma0 = mixer1.sigma0;
199			// xa2 and xai2 for j-i pairs are reversed from the same i-j pairing.
200			// Swapping the particles reverses the anisotropy parameters:
201			mixer2.xa2 = mixer1.xai2;
202			mixer2.xai2 = mixer1.xa2;
203			mixer2.x2 = mixer1.x2;
204	gezelter	1483
205			// assumed LB mixing rules for now:
206
207	gezelter	1674	mixer1.dw = 0.5 * (dw1 + dw2);
208	gezelter	1688	mixer1.eps0 = sqrt(eX1 * eX2);
209	gezelter	1674
210			mixer2.dw = mixer1.dw;
211			mixer2.eps0 = mixer1.eps0;
212	gezelter	1688
213			RealType mi = RealType(1.0)/mu_;
214	gezelter	1483
215	gezelter	1688	mixer1.xpap2 = (pow(eS1, mi) - pow(eE1, mi)) / (pow(eS1, mi) + pow(eE2, mi));
216			mixer1.xpapi2 = (pow(eS2, mi) - pow(eE2, mi)) / (pow(eS2, mi) + pow(eE1, mi));
217			mixer1.xp2 = (pow(eS1, mi) - pow(eE1, mi)) * (pow(eS2, mi) - pow(eE2, mi)) /
218			(pow(eS2, mi) + pow(eE1, mi)) / (pow(eS1, mi) + pow(eE2, mi)) ;
219	gezelter	1483
220	gezelter	1688	// xpap2 and xpapi2 for j-i pairs are reversed from the same i-j pairing.
221			// Swapping the particles reverses the anisotropy parameters:
222			mixer2.xpap2 = mixer1.xpapi2;
223			mixer2.xpapi2 = mixer1.xpap2;
224	gezelter	1674	mixer2.xp2 = mixer1.xp2;
225	gezelter	1587
226	gezelter	1483	// only add this pairing if at least one of the atoms is a Gay-Berne atom
227
228	gezelter	1710	if (gba1.isGayBerne() \|\| gba2.isGayBerne()) {
229	gezelter	1483
230			pair<AtomType, AtomType> key1, key2;
231			key1 = make_pair(atomType, atype2);
232			key2 = make_pair(atype2, atomType);
233
234	gezelter	1674	MixingMap[key1] = mixer1;
235	gezelter	1483	if (key2 != key1) {
236	gezelter	1674	MixingMap[key2] = mixer2;
237	gezelter	1483	}
238			}
239			}
240			}
241	gezelter	1502
242	gezelter	1536	void GB::calcForce(InteractionData &idat) {
243	gezelter	1483
244			if (!initialized_) initialize();
245
246	gezelter	1571	GBInteractionData mixer = MixingMap[idat.atypes];
247	gezelter	1483
248			RealType sigma0 = mixer.sigma0;
249			RealType dw = mixer.dw;
250			RealType eps0 = mixer.eps0;
251			RealType x2 = mixer.x2;
252			RealType xa2 = mixer.xa2;
253			RealType xai2 = mixer.xai2;
254			RealType xp2 = mixer.xp2;
255			RealType xpap2 = mixer.xpap2;
256			RealType xpapi2 = mixer.xpapi2;
257
258	gezelter	1688	// cerr << "atypes = " << idat.atypes.first->getName() << " " << idat.atypes.second->getName() << "\n";
259			// cerr << "sigma0 = " <<mixer.sigma0 <<"\n";
260			// cerr << "dw = " <<mixer.dw <<"\n";
261			// cerr << "eps0 = " <<mixer.eps0 <<"\n";
262			// cerr << "x2 = " <<mixer.x2 <<"\n";
263			// cerr << "xa2 = " <<mixer.xa2 <<"\n";
264			// cerr << "xai2 = " <<mixer.xai2 <<"\n";
265			// cerr << "xp2 = " <<mixer.xp2 <<"\n";
266			// cerr << "xpap2 = " <<mixer.xpap2 <<"\n";
267			// cerr << "xpapi2 = " <<mixer.xpapi2 <<"\n";
268
269	gezelter	1554	Vector3d ul1 = idat.A1->getRow(2);
270			Vector3d ul2 = idat.A2->getRow(2);
271	gezelter	1483
272	gezelter	1688	// cerr << "ul1 = " <<ul1<<"\n";
273			// cerr << "ul2 = " <<ul2<<"\n";
274
275	gezelter	1483	RealType a, b, g;
276
277	gezelter	1710	// This is not good. We should store this in the mixing map, and not
278			// query atom types in calc force.
279	gezelter	1571	bool i_is_LJ = idat.atypes.first->isLennardJones();
280			bool j_is_LJ = idat.atypes.second->isLennardJones();
281	gezelter	1554
282	gezelter	1483	if (i_is_LJ) {
283			a = 0.0;
284			ul1 = V3Zero;
285			} else {
286	gezelter	1554	a = dot(*(idat.d), ul1);
287	gezelter	1483	}
288
289			if (j_is_LJ) {
290			b = 0.0;
291			ul2 = V3Zero;
292			} else {
293	gezelter	1554	b = dot(*(idat.d), ul2);
294	gezelter	1483	}
295
296			if (i_is_LJ \|\| j_is_LJ)
297			g = 0.0;
298			else
299			g = dot(ul1, ul2);
300
301	gezelter	1554	RealType au = a / *(idat.rij);
302			RealType bu = b / *(idat.rij);
303	gezelter	1483
304			RealType au2 = au * au;
305			RealType bu2 = bu * bu;
306			RealType g2 = g * g;
307	gezelter	1688
308	gezelter	1483	RealType H = (xa2 * au2 + xai2 * bu2 - 2.0x2aubug) / (1.0 - x2*g2);
309			RealType Hp = (xpap2au2 + xpapi2bu2 - 2.0xp2aubug) / (1.0 - xp2*g2);
310
311	gezelter	1688	// cerr << "au2 = " << au2 << "\n";
312			// cerr << "bu2 = " << bu2 << "\n";
313			// cerr << "g2 = " << g2 << "\n";
314			// cerr << "H = " << H << "\n";
315			// cerr << "Hp = " << Hp << "\n";
316
317	gezelter	1483	RealType sigma = sigma0 / sqrt(1.0 - H);
318			RealType e1 = 1.0 / sqrt(1.0 - x2*g2);
319			RealType e2 = 1.0 - Hp;
320			RealType eps = eps0 * pow(e1,nu_) * pow(e2,mu_);
321	gezelter	1554	RealType BigR = dwsigma0 / ((idat.rij) - sigma + dw*sigma0);
322	gezelter	1483
323			RealType R3 = BigRBigRBigR;
324			RealType R6 = R3*R3;
325			RealType R7 = R6 * BigR;
326			RealType R12 = R6*R6;
327			RealType R13 = R6*R7;
328
329	gezelter	1554	RealType U = (idat.vdwMult) 4.0 * eps * (R12 - R6);
330	gezelter	1483
331			RealType s3 = sigmasigmasigma;
332			RealType s03 = sigma0sigma0sigma0;
333
334	gezelter	1688	// cerr << "vdwMult = " << *(idat.vdwMult) << "\n";
335			// cerr << "eps = " << eps <<"\n";
336			// cerr << "mu = " << mu_ << "\n";
337			// cerr << "R12 = " << R12 << "\n";
338			// cerr << "R6 = " << R6 << "\n";
339			// cerr << "R13 = " << R13 << "\n";
340			// cerr << "R7 = " << R7 << "\n";
341			// cerr << "e2 = " << e2 << "\n";
342			// cerr << "rij = " << *(idat.rij) << "\n";
343			// cerr << "s3 = " << s3 << "\n";
344			// cerr << "s03 = " << s03 << "\n";
345			// cerr << "dw = " << dw << "\n";
346
347	gezelter	1554	RealType pref1 = - (idat.vdwMult) 8.0 * eps * mu_ * (R12 - R6) /
348			(e2 * *(idat.rij));
349	gezelter	1483
350	gezelter	1554	RealType pref2 = (idat.vdwMult) 8.0 * eps * s3 * (6.0R13 - 3.0R7) /
351			(dw* (idat.rij) s03);
352	gezelter	1483
353	gezelter	1554	RealType dUdr = - (pref1 * Hp + pref2 * (sigma0 * sigma0 *
354			*(idat.rij) / s3 + H));
355	gezelter	1483
356			RealType dUda = pref1 * (xpap2au - xp2bug) / (1.0 - xp2 g2)
357			+ pref2 * (xa2 * au - x2 bug) / (1.0 - x2 * g2);
358
359			RealType dUdb = pref1 * (xpapi2bu - xp2aug) / (1.0 - xp2 g2)
360			+ pref2 * (xai2 * bu - x2 aug) / (1.0 - x2 * g2);
361
362			RealType dUdg = 4.0 * eps * nu_ * (R12 - R6) * x2 * g / (1.0 - x2*g2)
363			+ 8.0 * eps * mu_ * (R12 - R6) * (xp2aubu - Hpxp2g) /
364			(1.0 - xp2 * g2) / e2 + 8.0 * eps * s3 * (3.0 * R7 - 6.0 * R13) *
365			(x2 * au * bu - H * x2 * g) / (1.0 - x2 * g2) / (dw * s03);
366
367	gezelter	1688	// cerr << "pref = " << pref1 << " " << pref2 << "\n";
368			// cerr << "dU = " << dUdr << " " << dUda <<" " << dUdb << " " << dUdg << "\n";
369
370	gezelter	1554	Vector3d rhat = (idat.d) / (idat.rij);
371			Vector3d rxu1 = cross(*(idat.d), ul1);
372			Vector3d rxu2 = cross(*(idat.d), ul2);
373	gezelter	1483	Vector3d uxu = cross(ul1, ul2);
374	gezelter	1587
375	gezelter	1582	((idat.pot))[VANDERWAALS_FAMILY] += U *(idat.sw);
376	gezelter	1686	(idat.f1) += (dUdr rhat + dUda * ul1 + dUdb * ul2) * *(idat.sw);
377			(idat.t1) += (dUda rxu1 - dUdg * uxu) * *(idat.sw);
378			(idat.t2) += (dUdb rxu2 + dUdg * uxu) * *(idat.sw);
379			*(idat.vpair) += U;
380	gezelter	1483
381	gezelter	1688	// cerr << "f1 term = " << (dUdr * rhat + dUda * ul1 + dUdb * ul2) * *(idat.sw) << "\n";
382			// cerr << "t1 term = " << (dUda * rxu1 - dUdg * uxu) * *(idat.sw) << "\n";
383			// cerr << "t2 term = " << (dUdb * rxu2 + dUdg * uxu) * *(idat.sw) << "\n";
384			// cerr << "vp term = " << U << "\n";
385
386	gezelter	1483	return;
387
388			}
389	gezelter	1505
390	gezelter	1545	RealType GB::getSuggestedCutoffRadius(pair<AtomType, AtomType> atypes) {
391	gezelter	1505	if (!initialized_) initialize();
392
393			RealType cut = 0.0;
394
395	gezelter	1710	LennardJonesAdapter lja1 = LennardJonesAdapter(atypes.first);
396			GayBerneAdapter gba1 = GayBerneAdapter(atypes.first);
397			LennardJonesAdapter lja2 = LennardJonesAdapter(atypes.second);
398			GayBerneAdapter gba2 = GayBerneAdapter(atypes.second);
399
400			if (gba1.isGayBerne()) {
401			RealType d1 = gba1.getD();
402			RealType l1 = gba1.getL();
403	gezelter	1505	// sigma is actually sqrt(2)*l for prolate ellipsoids
404	gezelter	1668	cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d1, l1));
405	gezelter	1710	} else if (lja1.isLennardJones()) {
406			cut = max(cut, RealType(2.5) * lja1.getSigma());
407	gezelter	1505	}
408
409	gezelter	1710	if (gba2.isGayBerne()) {
410			RealType d2 = gba2.getD();
411			RealType l2 = gba2.getL();
412	gezelter	1668	cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d2, l2));
413	gezelter	1710	} else if (lja2.isLennardJones()) {
414			cut = max(cut, RealType(2.5) * lja2.getSigma());
415	gezelter	1505	}
416
417			return cut;
418			}
419	gezelter	1483	}
420