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, 234107 (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() {    
    
    GBtypes.clear();
    GBtids.clear();
    MixingMap.clear();
    nGB_ = 0;

    GBtids.resize( forceField_->getNAtomType(), -1);

    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
    mu_ = fopts.getGayBerneMu();
    nu_ = fopts.getGayBerneNu();

    // GB handles all of the GB-GB interactions as well as GB-LJ cross
    // interactions:
    set<AtomType*>::iterator at;
    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
      if ((*at)->isGayBerne()) nGB_++;
      if ((*at)->isLennardJones()) nGB_++;
    }

    MixingMap.resize(nGB_);
    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
      if ((*at)->isGayBerne() || (*at)->isLennardJones()) addType( *at );
    }
    
    initialized_ = true;
  }
      
  void GB::addType(AtomType* atomType){

    // add it to the map:
    int atid = atomType->getIdent();
    int gbtid = GBtypes.size();
  
    pair<set<int>::iterator,bool> ret;    
    ret = GBtypes.insert( atid );
    if (ret.second == false) {
      sprintf( painCave.errMsg,
               "GB already had a previous entry with ident %d\n",
               atid) ;
      painCave.severity = OPENMD_INFO;
      painCave.isFatal = 0;
      simError();         
    }

    GBtids[atid] = gbtid;
    MixingMap[gbtid].resize( nGB_ );
    
    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:
    
    std::set<int>::iterator it;
    for( it = GBtypes.begin(); it != GBtypes.end(); ++it) {
      
      int gbtid2 = GBtids[ (*it) ];
      AtomType* atype2 = forceField_->getAtomType( (*it) );

      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;
      } else {
        sprintf( painCave.errMsg,
                 "GB::addType found an atomType (%s) that does not\n"
                 "\tappear to be a Gay-Berne or Lennard-Jones atom.\n",
                 atype2->getName().c_str());
        painCave.severity = OPENMD_ERROR;
        painCave.isFatal = 1;
        simError();
      }
      
      
      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;
      // keep track of who is the LJ atom:
      mixer1.i_is_LJ = atomType->isLennardJones();
      mixer1.j_is_LJ = atype2->isLennardJones();
      mixer2.i_is_LJ = mixer1.j_is_LJ;
      mixer2.j_is_LJ = mixer1.i_is_LJ;


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

      if (gba1.isGayBerne() || gba2.isGayBerne()) {
        MixingMap[gbtid2].resize( nGB_ );        
        MixingMap[gbtid][gbtid2] = mixer1;
        if (gbtid2 != gbtid)  {
          MixingMap[gbtid2][gbtid] = mixer2;
        }          
      }
    }
  }
   
  void GB::calcForce(InteractionData &idat) {

    if (!initialized_) initialize();
    
    GBInteractionData &mixer = MixingMap[GBtids[idat.atid1]][GBtids[idat.atid2]];

    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;

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

    RealType a, b, g;
    
    if (mixer.i_is_LJ) {
      a = 0.0;
      ul1 = V3Zero;
    } else {
      a = dot(*(idat.d), ul1);
    }

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

    if (mixer.i_is_LJ || mixer.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);

    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;

    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);
    
    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;

    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:	1930
Committed:	Mon Aug 19 13:51:04 2013 UTC (11 years, 11 months ago) by gezelter
File size:	13999 byte(s)
Log Message:	region fixes, performance boosts
#	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	gezelter	1879	* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (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	1895	GBtypes.clear();
96			GBtids.clear();
97			MixingMap.clear();
98			nGB_ = 0;
99
100			GBtids.resize( forceField_->getNAtomType(), -1);
101
102	gezelter	1485	ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
103			mu_ = fopts.getGayBerneMu();
104			nu_ = fopts.getGayBerneNu();
105	gezelter	1483
106			// GB handles all of the GB-GB interactions as well as GB-LJ cross
107			// interactions:
108	gezelter	1895	set<AtomType*>::iterator at;
109			for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
110			if ((*at)->isGayBerne()) nGB_++;
111			if ((*at)->isLennardJones()) nGB_++;
112			}
113	gezelter	1483
114	gezelter	1895	MixingMap.resize(nGB_);
115			for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
116			if ((at)->isGayBerne() \|\| (at)->isLennardJones()) addType( *at );
117	gezelter	1483	}
118	gezelter	1895
119	gezelter	1483	initialized_ = true;
120			}
121
122			void GB::addType(AtomType* atomType){
123	gezelter	1895
124	gezelter	1483	// add it to the map:
125	gezelter	1895	int atid = atomType->getIdent();
126			int gbtid = GBtypes.size();
127
128			pair<set<int>::iterator,bool> ret;
129			ret = GBtypes.insert( atid );
130	gezelter	1483	if (ret.second == false) {
131			sprintf( painCave.errMsg,
132			"GB already had a previous entry with ident %d\n",
133	gezelter	1895	atid) ;
134	gezelter	1483	painCave.severity = OPENMD_INFO;
135			painCave.isFatal = 0;
136			simError();
137			}
138	gezelter	1895
139			GBtids[atid] = gbtid;
140			MixingMap[gbtid].resize( nGB_ );
141	gezelter	1483
142	gezelter	1688	RealType d1, l1, eX1, eS1, eE1, dw1;
143	gezelter	1895
144	gezelter	1710	LennardJonesAdapter lja1 = LennardJonesAdapter(atomType);
145			GayBerneAdapter gba1 = GayBerneAdapter(atomType);
146			if (gba1.isGayBerne()) {
147			d1 = gba1.getD();
148			l1 = gba1.getL();
149			eX1 = gba1.getEpsX();
150			eS1 = gba1.getEpsS();
151			eE1 = gba1.getEpsE();
152			dw1 = gba1.getDw();
153			} else if (lja1.isLennardJones()) {
154			d1 = lja1.getSigma() / sqrt(2.0);
155	gezelter	1483	l1 = d1;
156	gezelter	1710	eX1 = lja1.getEpsilon();
157	gezelter	1688	eS1 = eX1;
158			eE1 = eX1;
159	gezelter	1483	dw1 = 1.0;
160			} else {
161			sprintf( painCave.errMsg,
162			"GB::addType was passed an atomType (%s) that does not\n"
163			"\tappear to be a Gay-Berne or Lennard-Jones atom.\n",
164			atomType->getName().c_str());
165			painCave.severity = OPENMD_ERROR;
166			painCave.isFatal = 1;
167			simError();
168			}
169
170	gezelter	1895
171	gezelter	1483	// Now, iterate over all known types and add to the mixing map:
172
173	gezelter	1895	std::set<int>::iterator it;
174			for( it = GBtypes.begin(); it != GBtypes.end(); ++it) {
175	gezelter	1483
176	gezelter	1895	int gbtid2 = GBtids[ (*it) ];
177			AtomType* atype2 = forceField_->getAtomType( (*it) );
178
179	gezelter	1710	LennardJonesAdapter lja2 = LennardJonesAdapter(atype2);
180			GayBerneAdapter gba2 = GayBerneAdapter(atype2);
181	gezelter	1688	RealType d2, l2, eX2, eS2, eE2, dw2;
182	gezelter	1483
183	gezelter	1710	if (gba2.isGayBerne()) {
184			d2 = gba2.getD();
185			l2 = gba2.getL();
186			eX2 = gba2.getEpsX();
187			eS2 = gba2.getEpsS();
188			eE2 = gba2.getEpsE();
189			dw2 = gba2.getDw();
190			} else if (lja2.isLennardJones()) {
191			d2 = lja2.getSigma() / sqrt(2.0);
192	gezelter	1483	l2 = d2;
193	gezelter	1710	eX2 = lja2.getEpsilon();
194	gezelter	1688	eS2 = eX2;
195			eE2 = eX2;
196	gezelter	1483	dw2 = 1.0;
197	gezelter	1879	} else {
198			sprintf( painCave.errMsg,
199			"GB::addType found an atomType (%s) that does not\n"
200			"\tappear to be a Gay-Berne or Lennard-Jones atom.\n",
201			atype2->getName().c_str());
202			painCave.severity = OPENMD_ERROR;
203			painCave.isFatal = 1;
204			simError();
205			}
206	gezelter	1895
207
208	gezelter	1674	GBInteractionData mixer1, mixer2;
209	gezelter	1483
210			// Cleaver paper uses sqrt of squares to get sigma0 for
211			// mixed interactions.
212	gezelter	1895
213	gezelter	1674	mixer1.sigma0 = sqrt(d1d1 + d2d2);
214			mixer1.xa2 = (l1l1 - d1d1)/(l1l1 + d2d2);
215			mixer1.xai2 = (l2l2 - d2d2)/(l2l2 + d1d1);
216			mixer1.x2 = (l1l1 - d1d1) * (l2l2 - d2d2) /
217	gezelter	1483	((l2l2 + d1d1) * (l1l1 + d2d2));
218	gezelter	1895
219	gezelter	1674	mixer2.sigma0 = mixer1.sigma0;
220			// xa2 and xai2 for j-i pairs are reversed from the same i-j pairing.
221			// Swapping the particles reverses the anisotropy parameters:
222			mixer2.xa2 = mixer1.xai2;
223			mixer2.xai2 = mixer1.xa2;
224			mixer2.x2 = mixer1.x2;
225	gezelter	1895
226	gezelter	1483	// assumed LB mixing rules for now:
227	gezelter	1895
228	gezelter	1674	mixer1.dw = 0.5 * (dw1 + dw2);
229	gezelter	1688	mixer1.eps0 = sqrt(eX1 * eX2);
230	gezelter	1674
231			mixer2.dw = mixer1.dw;
232			mixer2.eps0 = mixer1.eps0;
233	gezelter	1895
234	gezelter	1688	RealType mi = RealType(1.0)/mu_;
235	gezelter	1483
236	gezelter	1688	mixer1.xpap2 = (pow(eS1, mi) - pow(eE1, mi)) / (pow(eS1, mi) + pow(eE2, mi));
237			mixer1.xpapi2 = (pow(eS2, mi) - pow(eE2, mi)) / (pow(eS2, mi) + pow(eE1, mi));
238			mixer1.xp2 = (pow(eS1, mi) - pow(eE1, mi)) * (pow(eS2, mi) - pow(eE2, mi)) /
239			(pow(eS2, mi) + pow(eE1, mi)) / (pow(eS1, mi) + pow(eE2, mi)) ;
240	gezelter	1895
241	gezelter	1688	// xpap2 and xpapi2 for j-i pairs are reversed from the same i-j pairing.
242			// Swapping the particles reverses the anisotropy parameters:
243			mixer2.xpap2 = mixer1.xpapi2;
244			mixer2.xpapi2 = mixer1.xpap2;
245	gezelter	1674	mixer2.xp2 = mixer1.xp2;
246	gezelter	1930	// keep track of who is the LJ atom:
247			mixer1.i_is_LJ = atomType->isLennardJones();
248			mixer1.j_is_LJ = atype2->isLennardJones();
249			mixer2.i_is_LJ = mixer1.j_is_LJ;
250			mixer2.j_is_LJ = mixer1.i_is_LJ;
251	gezelter	1587
252	gezelter	1930
253	gezelter	1483	// only add this pairing if at least one of the atoms is a Gay-Berne atom
254
255	gezelter	1710	if (gba1.isGayBerne() \|\| gba2.isGayBerne()) {
256	gezelter	1895	MixingMap[gbtid2].resize( nGB_ );
257			MixingMap[gbtid][gbtid2] = mixer1;
258	gezelter	1930	if (gbtid2 != gbtid) {
259	gezelter	1895	MixingMap[gbtid2][gbtid] = mixer2;
260	gezelter	1930	}
261	gezelter	1483	}
262	gezelter	1895	}
263	gezelter	1483	}
264	gezelter	1502
265	gezelter	1536	void GB::calcForce(InteractionData &idat) {
266	gezelter	1483
267			if (!initialized_) initialize();
268
269	gezelter	1895	GBInteractionData &mixer = MixingMap[GBtids[idat.atid1]][GBtids[idat.atid2]];
270	gezelter	1483
271			RealType sigma0 = mixer.sigma0;
272			RealType dw = mixer.dw;
273			RealType eps0 = mixer.eps0;
274			RealType x2 = mixer.x2;
275			RealType xa2 = mixer.xa2;
276			RealType xai2 = mixer.xai2;
277			RealType xp2 = mixer.xp2;
278			RealType xpap2 = mixer.xpap2;
279			RealType xpapi2 = mixer.xpapi2;
280
281	gezelter	1554	Vector3d ul1 = idat.A1->getRow(2);
282			Vector3d ul2 = idat.A2->getRow(2);
283	gezelter	1483
284			RealType a, b, g;
285	gezelter	1554
286	gezelter	1930	if (mixer.i_is_LJ) {
287	gezelter	1483	a = 0.0;
288			ul1 = V3Zero;
289			} else {
290	gezelter	1554	a = dot(*(idat.d), ul1);
291	gezelter	1483	}
292
293	gezelter	1930	if (mixer.j_is_LJ) {
294	gezelter	1483	b = 0.0;
295			ul2 = V3Zero;
296			} else {
297	gezelter	1554	b = dot(*(idat.d), ul2);
298	gezelter	1483	}
299
300	gezelter	1930	if (mixer.i_is_LJ \|\| mixer.j_is_LJ)
301	gezelter	1483	g = 0.0;
302			else
303			g = dot(ul1, ul2);
304
305	gezelter	1554	RealType au = a / *(idat.rij);
306			RealType bu = b / *(idat.rij);
307	gezelter	1483
308			RealType au2 = au * au;
309			RealType bu2 = bu * bu;
310			RealType g2 = g * g;
311	gezelter	1688
312	gezelter	1483	RealType H = (xa2 * au2 + xai2 * bu2 - 2.0x2aubug) / (1.0 - x2*g2);
313			RealType Hp = (xpap2au2 + xpapi2bu2 - 2.0xp2aubug) / (1.0 - xp2*g2);
314
315			RealType sigma = sigma0 / sqrt(1.0 - H);
316			RealType e1 = 1.0 / sqrt(1.0 - x2*g2);
317			RealType e2 = 1.0 - Hp;
318			RealType eps = eps0 * pow(e1,nu_) * pow(e2,mu_);
319	gezelter	1554	RealType BigR = dwsigma0 / ((idat.rij) - sigma + dw*sigma0);
320	gezelter	1483
321			RealType R3 = BigRBigRBigR;
322			RealType R6 = R3*R3;
323			RealType R7 = R6 * BigR;
324			RealType R12 = R6*R6;
325			RealType R13 = R6*R7;
326
327	gezelter	1554	RealType U = (idat.vdwMult) 4.0 * eps * (R12 - R6);
328	gezelter	1483
329			RealType s3 = sigmasigmasigma;
330			RealType s03 = sigma0sigma0sigma0;
331
332	gezelter	1554	RealType pref1 = - (idat.vdwMult) 8.0 * eps * mu_ * (R12 - R6) /
333			(e2 * *(idat.rij));
334	gezelter	1483
335	gezelter	1554	RealType pref2 = (idat.vdwMult) 8.0 * eps * s3 * (6.0R13 - 3.0R7) /
336			(dw* (idat.rij) s03);
337	gezelter	1483
338	gezelter	1554	RealType dUdr = - (pref1 * Hp + pref2 * (sigma0 * sigma0 *
339			*(idat.rij) / s3 + H));
340	gezelter	1483
341			RealType dUda = pref1 * (xpap2au - xp2bug) / (1.0 - xp2 g2)
342			+ pref2 * (xa2 * au - x2 bug) / (1.0 - x2 * g2);
343
344			RealType dUdb = pref1 * (xpapi2bu - xp2aug) / (1.0 - xp2 g2)
345			+ pref2 * (xai2 * bu - x2 aug) / (1.0 - x2 * g2);
346	gezelter	1895
347	gezelter	1483	RealType dUdg = 4.0 * eps * nu_ * (R12 - R6) * x2 * g / (1.0 - x2*g2)
348			+ 8.0 * eps * mu_ * (R12 - R6) * (xp2aubu - Hpxp2g) /
349			(1.0 - xp2 * g2) / e2 + 8.0 * eps * s3 * (3.0 * R7 - 6.0 * R13) *
350			(x2 * au * bu - H * x2 * g) / (1.0 - x2 * g2) / (dw * s03);
351	gezelter	1895
352	gezelter	1554	Vector3d rhat = (idat.d) / (idat.rij);
353			Vector3d rxu1 = cross(*(idat.d), ul1);
354			Vector3d rxu2 = cross(*(idat.d), ul2);
355	gezelter	1483	Vector3d uxu = cross(ul1, ul2);
356	gezelter	1895
357	gezelter	1582	((idat.pot))[VANDERWAALS_FAMILY] += U *(idat.sw);
358	gezelter	1686	(idat.f1) += (dUdr rhat + dUda * ul1 + dUdb * ul2) * *(idat.sw);
359			(idat.t1) += (dUda rxu1 - dUdg * uxu) * *(idat.sw);
360			(idat.t2) += (dUdb rxu2 + dUdg * uxu) * *(idat.sw);
361			*(idat.vpair) += U;
362	gezelter	1483
363			return;
364
365			}
366	gezelter	1505
367	gezelter	1545	RealType GB::getSuggestedCutoffRadius(pair<AtomType, AtomType> atypes) {
368	gezelter	1505	if (!initialized_) initialize();
369
370			RealType cut = 0.0;
371
372	gezelter	1710	LennardJonesAdapter lja1 = LennardJonesAdapter(atypes.first);
373			GayBerneAdapter gba1 = GayBerneAdapter(atypes.first);
374			LennardJonesAdapter lja2 = LennardJonesAdapter(atypes.second);
375			GayBerneAdapter gba2 = GayBerneAdapter(atypes.second);
376
377			if (gba1.isGayBerne()) {
378			RealType d1 = gba1.getD();
379			RealType l1 = gba1.getL();
380	gezelter	1505	// sigma is actually sqrt(2)*l for prolate ellipsoids
381	gezelter	1668	cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d1, l1));
382	gezelter	1710	} else if (lja1.isLennardJones()) {
383			cut = max(cut, RealType(2.5) * lja1.getSigma());
384	gezelter	1505	}
385
386	gezelter	1710	if (gba2.isGayBerne()) {
387			RealType d2 = gba2.getD();
388			RealType l2 = gba2.getL();
389	gezelter	1668	cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d2, l2));
390	gezelter	1710	} else if (lja2.isLennardJones()) {
391			cut = max(cut, RealType(2.5) * lja2.getSigma());
392	gezelter	1505	}
393
394			return cut;
395			}
396	gezelter	1483	}
397