src/nonbonded/GB.cpp

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

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

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

using namespace std;
namespace OpenMD {

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

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

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

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

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

    // Now, iterate over all known types and add to the mixing map:
    
    map<int, AtomType*>::iterator it;
    for( it = GBMap.begin(); it != GBMap.end(); ++it) {
      
      AtomType* atype2 = (*it).second;
      
      RealType d2, l2, e2, er2, dw2;
      
      if (atype2->isGayBerne()) {
        GayBerneParam gb2 = getGayBerneParam(atype2);
        d2 = gb2.GB_d;
        l2 = gb2.GB_l;
        e2 = gb2.GB_eps;
        er2 = gb2.GB_eps_ratio;
        dw2 = gb2.GB_dw;
      } else if (atype2->isLennardJones()) {
        d2 = getLJSigma(atype2) / sqrt(2.0);
        e2 = getLJEpsilon(atype2);
        l2 = d2;
        er2 = 1.0;
        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(e1 * e2);

      mixer2.dw = mixer1.dw;
      mixer2.eps0 = mixer1.eps0;
      
      RealType er = sqrt(er1 * er2);
      RealType ermu = pow(er, (RealType(1.0) / mu_));
      RealType xp = (1.0 - ermu) / (1.0 + ermu);
      RealType ap2 = 1.0 / (1.0 + ermu);
      
      mixer1.xp2 = xp * xp;
      mixer1.xpap2 = xp * ap2;
      mixer1.xpapi2 = xp / ap2;

      mixer2.xp2 = mixer1.xp2;
      mixer2.xpap2 = mixer1.xpap2;
      mixer2.xpapi2 = mixer1.xpapi2;

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

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

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

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

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

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

    RealType a, b, g;

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

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

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

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

    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;

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

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

Revision:	1686
Committed:	Sat Mar 10 04:21:44 2012 UTC (13 years, 4 months ago) by gezelter
File size:	13729 byte(s)
Log Message:	Fixes for compilation under LLVM compiler Bug-fix for switching function in Gay-Berne
#	User	Rev	Content
1	gezelter	1483	/*
2			* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9			* 1. Redistributions of source code must retain the above copyright
10			* notice, this list of conditions and the following disclaimer.
11			*
12			* 2. Redistributions in binary form must reproduce the above copyright
13			* notice, this list of conditions and the following disclaimer in the
14			* documentation and/or other materials provided with the
15			* distribution.
16			*
17			* This software is provided "AS IS," without a warranty of any
18			* kind. All express or implied conditions, representations and
19			* warranties, including any implied warranty of merchantability,
20			* fitness for a particular purpose or non-infringement, are hereby
21			* excluded. The University of Notre Dame and its licensors shall not
22			* be liable for any damages suffered by licensee as a result of
23			* using, modifying or distributing the software or its
24			* derivatives. In no event will the University of Notre Dame or its
25			* licensors be liable for any lost revenue, profit or data, or for
26			* direct, indirect, special, consequential, incidental or punitive
27			* damages, however caused and regardless of the theory of liability,
28			* arising out of the use of or inability to use software, even if the
29			* University of Notre Dame has been advised of the possibility of
30			* such damages.
31			*
32			* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33			* research, please cite the appropriate papers when you publish your
34			* work. Good starting points are:
35			*
36			* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37			* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38			* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	gezelter	1665	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	gezelter	1483	*/
42
43			#include <stdio.h>
44			#include <string.h>
45
46			#include <cmath>
47			#include "nonbonded/GB.hpp"
48			#include "utils/simError.h"
49
50			using namespace std;
51			namespace OpenMD {
52
53	gezelter	1502	GB::GB() : name_("GB"), initialized_(false), mu_(2.0), nu_(1.0), forceField_(NULL) {}
54	gezelter	1483
55			GayBerneParam GB::getGayBerneParam(AtomType* atomType) {
56
57			// Do sanity checking on the AtomType we were passed before
58			// building any data structures:
59			if (!atomType->isGayBerne()) {
60			sprintf( painCave.errMsg,
61			"GB::getGayBerneParam was passed an atomType (%s) that does\n"
62			"\tnot appear to be a Gay-Berne atom.\n",
63			atomType->getName().c_str());
64			painCave.severity = OPENMD_ERROR;
65			painCave.isFatal = 1;
66			simError();
67			}
68
69			DirectionalAtomType* daType = dynamic_cast<DirectionalAtomType*>(atomType);
70			GenericData* data = daType->getPropertyByName("GayBerne");
71			if (data == NULL) {
72			sprintf( painCave.errMsg, "GB::getGayBerneParam could not find\n"
73			"\tGay-Berne parameters for atomType %s.\n",
74			daType->getName().c_str());
75			painCave.severity = OPENMD_ERROR;
76			painCave.isFatal = 1;
77			simError();
78			}
79
80			GayBerneParamGenericData* gbData = dynamic_cast<GayBerneParamGenericData*>(data);
81			if (gbData == NULL) {
82			sprintf( painCave.errMsg,
83			"GB::getGayBerneParam could not convert GenericData to\n"
84			"\tGayBerneParamGenericData for atom type %s\n",
85			daType->getName().c_str());
86			painCave.severity = OPENMD_ERROR;
87			painCave.isFatal = 1;
88			simError();
89			}
90
91			return gbData->getData();
92			}
93
94	gezelter	1502	LJParam GB::getLJParam(AtomType* atomType) {
95
96			// Do sanity checking on the AtomType we were passed before
97			// building any data structures:
98			if (!atomType->isLennardJones()) {
99			sprintf( painCave.errMsg,
100			"GB::getLJParam was passed an atomType (%s) that does not\n"
101			"\tappear to be a Lennard-Jones atom.\n",
102			atomType->getName().c_str());
103			painCave.severity = OPENMD_ERROR;
104			painCave.isFatal = 1;
105			simError();
106			}
107
108			GenericData* data = atomType->getPropertyByName("LennardJones");
109			if (data == NULL) {
110			sprintf( painCave.errMsg, "GB::getLJParam could not find Lennard-Jones\n"
111			"\tparameters for atomType %s.\n", atomType->getName().c_str());
112			painCave.severity = OPENMD_ERROR;
113			painCave.isFatal = 1;
114			simError();
115			}
116
117			LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
118			if (ljData == NULL) {
119			sprintf( painCave.errMsg,
120			"GB::getLJParam could not convert GenericData to LJParam for\n"
121			"\tatom type %s\n", atomType->getName().c_str());
122			painCave.severity = OPENMD_ERROR;
123			painCave.isFatal = 1;
124			simError();
125			}
126
127			return ljData->getData();
128			}
129
130			RealType GB::getLJEpsilon(AtomType* atomType) {
131			LJParam ljParam = getLJParam(atomType);
132			return ljParam.epsilon;
133			}
134			RealType GB::getLJSigma(AtomType* atomType) {
135			LJParam ljParam = getLJParam(atomType);
136			return ljParam.sigma;
137			}
138
139	gezelter	1483	void GB::initialize() {
140	gezelter	1502
141	gezelter	1485	ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
142			mu_ = fopts.getGayBerneMu();
143			nu_ = fopts.getGayBerneNu();
144	gezelter	1483	ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
145			ForceField::AtomTypeContainer::MapTypeIterator i;
146			AtomType* at;
147
148			// GB handles all of the GB-GB interactions as well as GB-LJ cross
149			// interactions:
150
151			for (at = atomTypes->beginType(i); at != NULL;
152			at = atomTypes->nextType(i)) {
153
154			if (at->isGayBerne() \|\| at->isLennardJones())
155			addType(at);
156			}
157
158			initialized_ = true;
159			}
160
161			void GB::addType(AtomType* atomType){
162			// add it to the map:
163			AtomTypeProperties atp = atomType->getATP();
164
165			pair<map<int,AtomType*>::iterator,bool> ret;
166			ret = GBMap.insert( pair<int, AtomType*>(atp.ident, atomType) );
167			if (ret.second == false) {
168			sprintf( painCave.errMsg,
169			"GB already had a previous entry with ident %d\n",
170			atp.ident);
171			painCave.severity = OPENMD_INFO;
172			painCave.isFatal = 0;
173			simError();
174			}
175
176			RealType d1, l1, e1, er1, dw1;
177
178			if (atomType->isGayBerne()) {
179			GayBerneParam gb1 = getGayBerneParam(atomType);
180			d1 = gb1.GB_d;
181			l1 = gb1.GB_l;
182			e1 = gb1.GB_eps;
183			er1 = gb1.GB_eps_ratio;
184			dw1 = gb1.GB_dw;
185			} else if (atomType->isLennardJones()) {
186	gezelter	1502	d1 = getLJSigma(atomType) / sqrt(2.0);
187			e1 = getLJEpsilon(atomType);
188	gezelter	1483	l1 = d1;
189			er1 = 1.0;
190			dw1 = 1.0;
191			} else {
192			sprintf( painCave.errMsg,
193			"GB::addType was passed an atomType (%s) that does not\n"
194			"\tappear to be a Gay-Berne or Lennard-Jones atom.\n",
195			atomType->getName().c_str());
196			painCave.severity = OPENMD_ERROR;
197			painCave.isFatal = 1;
198			simError();
199			}
200
201
202			// Now, iterate over all known types and add to the mixing map:
203
204			map<int, AtomType*>::iterator it;
205			for( it = GBMap.begin(); it != GBMap.end(); ++it) {
206
207			AtomType* atype2 = (*it).second;
208
209			RealType d2, l2, e2, er2, dw2;
210
211			if (atype2->isGayBerne()) {
212			GayBerneParam gb2 = getGayBerneParam(atype2);
213			d2 = gb2.GB_d;
214			l2 = gb2.GB_l;
215			e2 = gb2.GB_eps;
216			er2 = gb2.GB_eps_ratio;
217			dw2 = gb2.GB_dw;
218			} else if (atype2->isLennardJones()) {
219	gezelter	1502	d2 = getLJSigma(atype2) / sqrt(2.0);
220			e2 = getLJEpsilon(atype2);
221	gezelter	1483	l2 = d2;
222			er2 = 1.0;
223			dw2 = 1.0;
224			}
225
226	gezelter	1674	GBInteractionData mixer1, mixer2;
227	gezelter	1483
228			// Cleaver paper uses sqrt of squares to get sigma0 for
229			// mixed interactions.
230
231	gezelter	1674	mixer1.sigma0 = sqrt(d1d1 + d2d2);
232			mixer1.xa2 = (l1l1 - d1d1)/(l1l1 + d2d2);
233			mixer1.xai2 = (l2l2 - d2d2)/(l2l2 + d1d1);
234			mixer1.x2 = (l1l1 - d1d1) * (l2l2 - d2d2) /
235	gezelter	1483	((l2l2 + d1d1) * (l1l1 + d2d2));
236	gezelter	1674
237			mixer2.sigma0 = mixer1.sigma0;
238			// xa2 and xai2 for j-i pairs are reversed from the same i-j pairing.
239			// Swapping the particles reverses the anisotropy parameters:
240			mixer2.xa2 = mixer1.xai2;
241			mixer2.xai2 = mixer1.xa2;
242			mixer2.x2 = mixer1.x2;
243	gezelter	1483
244			// assumed LB mixing rules for now:
245
246	gezelter	1674	mixer1.dw = 0.5 * (dw1 + dw2);
247			mixer1.eps0 = sqrt(e1 * e2);
248
249			mixer2.dw = mixer1.dw;
250			mixer2.eps0 = mixer1.eps0;
251	gezelter	1483
252			RealType er = sqrt(er1 * er2);
253	gezelter	1668	RealType ermu = pow(er, (RealType(1.0) / mu_));
254	gezelter	1483	RealType xp = (1.0 - ermu) / (1.0 + ermu);
255			RealType ap2 = 1.0 / (1.0 + ermu);
256
257	gezelter	1674	mixer1.xp2 = xp * xp;
258			mixer1.xpap2 = xp * ap2;
259			mixer1.xpapi2 = xp / ap2;
260	gezelter	1483
261	gezelter	1674	mixer2.xp2 = mixer1.xp2;
262			mixer2.xpap2 = mixer1.xpap2;
263			mixer2.xpapi2 = mixer1.xpapi2;
264	gezelter	1587
265	gezelter	1483	// only add this pairing if at least one of the atoms is a Gay-Berne atom
266
267			if (atomType->isGayBerne() \|\| atype2->isGayBerne()) {
268
269			pair<AtomType, AtomType> key1, key2;
270			key1 = make_pair(atomType, atype2);
271			key2 = make_pair(atype2, atomType);
272
273	gezelter	1674	MixingMap[key1] = mixer1;
274	gezelter	1483	if (key2 != key1) {
275	gezelter	1674	MixingMap[key2] = mixer2;
276	gezelter	1483	}
277			}
278			}
279			}
280	gezelter	1502
281	gezelter	1536	void GB::calcForce(InteractionData &idat) {
282	gezelter	1483
283			if (!initialized_) initialize();
284
285	gezelter	1571	GBInteractionData mixer = MixingMap[idat.atypes];
286	gezelter	1483
287			RealType sigma0 = mixer.sigma0;
288			RealType dw = mixer.dw;
289			RealType eps0 = mixer.eps0;
290			RealType x2 = mixer.x2;
291			RealType xa2 = mixer.xa2;
292			RealType xai2 = mixer.xai2;
293			RealType xp2 = mixer.xp2;
294			RealType xpap2 = mixer.xpap2;
295			RealType xpapi2 = mixer.xpapi2;
296
297	gezelter	1554	Vector3d ul1 = idat.A1->getRow(2);
298			Vector3d ul2 = idat.A2->getRow(2);
299	gezelter	1483
300			RealType a, b, g;
301
302	gezelter	1571	bool i_is_LJ = idat.atypes.first->isLennardJones();
303			bool j_is_LJ = idat.atypes.second->isLennardJones();
304	gezelter	1554
305	gezelter	1483	if (i_is_LJ) {
306			a = 0.0;
307			ul1 = V3Zero;
308			} else {
309	gezelter	1554	a = dot(*(idat.d), ul1);
310	gezelter	1483	}
311
312			if (j_is_LJ) {
313			b = 0.0;
314			ul2 = V3Zero;
315			} else {
316	gezelter	1554	b = dot(*(idat.d), ul2);
317	gezelter	1483	}
318
319			if (i_is_LJ \|\| j_is_LJ)
320			g = 0.0;
321			else
322			g = dot(ul1, ul2);
323
324	gezelter	1554	RealType au = a / *(idat.rij);
325			RealType bu = b / *(idat.rij);
326	gezelter	1483
327			RealType au2 = au * au;
328			RealType bu2 = bu * bu;
329			RealType g2 = g * g;
330
331			RealType H = (xa2 * au2 + xai2 * bu2 - 2.0x2aubug) / (1.0 - x2*g2);
332			RealType Hp = (xpap2au2 + xpapi2bu2 - 2.0xp2aubug) / (1.0 - xp2*g2);
333
334			RealType sigma = sigma0 / sqrt(1.0 - H);
335			RealType e1 = 1.0 / sqrt(1.0 - x2*g2);
336			RealType e2 = 1.0 - Hp;
337			RealType eps = eps0 * pow(e1,nu_) * pow(e2,mu_);
338	gezelter	1554	RealType BigR = dwsigma0 / ((idat.rij) - sigma + dw*sigma0);
339	gezelter	1483
340			RealType R3 = BigRBigRBigR;
341			RealType R6 = R3*R3;
342			RealType R7 = R6 * BigR;
343			RealType R12 = R6*R6;
344			RealType R13 = R6*R7;
345
346	gezelter	1554	RealType U = (idat.vdwMult) 4.0 * eps * (R12 - R6);
347	gezelter	1483
348			RealType s3 = sigmasigmasigma;
349			RealType s03 = sigma0sigma0sigma0;
350
351	gezelter	1554	RealType pref1 = - (idat.vdwMult) 8.0 * eps * mu_ * (R12 - R6) /
352			(e2 * *(idat.rij));
353	gezelter	1483
354	gezelter	1554	RealType pref2 = (idat.vdwMult) 8.0 * eps * s3 * (6.0R13 - 3.0R7) /
355			(dw* (idat.rij) s03);
356	gezelter	1483
357	gezelter	1554	RealType dUdr = - (pref1 * Hp + pref2 * (sigma0 * sigma0 *
358			*(idat.rij) / s3 + H));
359	gezelter	1483
360			RealType dUda = pref1 * (xpap2au - xp2bug) / (1.0 - xp2 g2)
361			+ pref2 * (xa2 * au - x2 bug) / (1.0 - x2 * g2);
362
363			RealType dUdb = pref1 * (xpapi2bu - xp2aug) / (1.0 - xp2 g2)
364			+ pref2 * (xai2 * bu - x2 aug) / (1.0 - x2 * g2);
365
366			RealType dUdg = 4.0 * eps * nu_ * (R12 - R6) * x2 * g / (1.0 - x2*g2)
367			+ 8.0 * eps * mu_ * (R12 - R6) * (xp2aubu - Hpxp2g) /
368			(1.0 - xp2 * g2) / e2 + 8.0 * eps * s3 * (3.0 * R7 - 6.0 * R13) *
369			(x2 * au * bu - H * x2 * g) / (1.0 - x2 * g2) / (dw * s03);
370
371	gezelter	1554	Vector3d rhat = (idat.d) / (idat.rij);
372			Vector3d rxu1 = cross(*(idat.d), ul1);
373			Vector3d rxu2 = cross(*(idat.d), ul2);
374	gezelter	1483	Vector3d uxu = cross(ul1, ul2);
375	gezelter	1587
376	gezelter	1582	((idat.pot))[VANDERWAALS_FAMILY] += U *(idat.sw);
377	gezelter	1686	(idat.f1) += (dUdr rhat + dUda * ul1 + dUdb * ul2) * *(idat.sw);
378			(idat.t1) += (dUda rxu1 - dUdg * uxu) * *(idat.sw);
379			(idat.t2) += (dUdb rxu2 + dUdg * uxu) * *(idat.sw);
380			*(idat.vpair) += U;
381	gezelter	1483
382			return;
383
384			}
385	gezelter	1505
386	gezelter	1545	RealType GB::getSuggestedCutoffRadius(pair<AtomType, AtomType> atypes) {
387	gezelter	1505	if (!initialized_) initialize();
388
389			RealType cut = 0.0;
390
391	gezelter	1545	if (atypes.first->isGayBerne()) {
392			GayBerneParam gb1 = getGayBerneParam(atypes.first);
393	gezelter	1505	RealType d1 = gb1.GB_d;
394			RealType l1 = gb1.GB_l;
395			// sigma is actually sqrt(2)*l for prolate ellipsoids
396	gezelter	1668	cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d1, l1));
397	gezelter	1545	} else if (atypes.first->isLennardJones()) {
398	gezelter	1668	cut = max(cut, RealType(2.5) * getLJSigma(atypes.first));
399	gezelter	1505	}
400
401	gezelter	1545	if (atypes.second->isGayBerne()) {
402			GayBerneParam gb2 = getGayBerneParam(atypes.second);
403	gezelter	1505	RealType d2 = gb2.GB_d;
404			RealType l2 = gb2.GB_l;
405	gezelter	1668	cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d2, l2));
406	gezelter	1545	} else if (atypes.second->isLennardJones()) {
407	gezelter	1668	cut = max(cut, RealType(2.5) * getLJSigma(atypes.second));
408	gezelter	1505	}
409
410			return cut;
411			}
412	gezelter	1483	}
413