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