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]  Vardeman & Gezelter, in progress (2009).                        
 */

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

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

using namespace std;
namespace OpenMD {

  bool GB::initialized_ = false;
  RealType GB::mu_ = 2.0;
  RealType GB::nu_ = 1.0;
  ForceField* GB::forceField_ = NULL;
  map<int, AtomType*> GB::GBMap;
  map<pair<AtomType*, AtomType*>, GBInteractionData> GB::MixingMap;
  
  GB* GB::_instance = NULL;

  GB* GB::Instance() {
    if (!_instance) {
      _instance = new GB();
    }
    return _instance;
  }

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

  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 = LJ::Instance()->getSigma(atomType) / sqrt(2.0);
      e1 = LJ::Instance()->getEpsilon(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 = LJ::Instance()->getSigma(atype2) / sqrt(2.0);
        e2 = LJ::Instance()->getEpsilon(atype2);
        l2 = d2;
        er2 = 1.0;
        dw2 = 1.0;
      } 
                       
      GBInteractionData mixer;        
      
      //  Cleaver paper uses sqrt of squares to get sigma0 for
      //  mixed interactions.
            
      mixer.sigma0 = sqrt(d1*d1 + d2*d2);
      mixer.xa2 = (l1*l1 - d1*d1)/(l1*l1 + d2*d2);
      mixer.xai2 = (l2*l2 - d2*d2)/(l2*l2 + d1*d1);
      mixer.x2 = (l1*l1 - d1*d1) * (l2*l2 - d2*d2) /
        ((l2*l2 + d1*d1) * (l1*l1 + d2*d2));
 
      // assumed LB mixing rules for now:

      mixer.dw = 0.5 * (dw1 + dw2);
      mixer.eps0 = sqrt(e1 * e2);
      
      RealType er = sqrt(er1 * er2);
      RealType ermu = pow(er,(1.0 / mu_));
      RealType xp = (1.0 - ermu) / (1.0 + ermu);
      RealType ap2 = 1.0 / (1.0 + ermu);
      
      mixer.xp2 = xp * xp;
      mixer.xpap2 = xp * ap2;
      mixer.xpapi2 = xp / ap2;

      // 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] = mixer;
        if (key2 != key1) {
          MixingMap[key2] = mixer;
        }
      }
    }      
  }
  

  RealType GB::getGayBerneCut(int atid) { 
    if (!initialized_) initialize();
    std::map<int, AtomType*> :: const_iterator it;
    it = GBMap.find(atid);
    if (it == GBMap.end()) {
      sprintf( painCave.errMsg,
               "GB::getGayBerneCut could not find atid %d in GBMap\n",
               (atid));
      painCave.severity = OPENMD_ERROR;
      painCave.isFatal = 1;
      simError();          
    }

    AtomType* atype = it->second;

    RealType gbCut;
    
    if (atype->isGayBerne()) {
      GayBerneParam gb = getGayBerneParam(atype);

      // sigma is actually sqrt(2) * l for prolate ellipsoids
      gbCut = 2.5 * sqrt(2.0) * max(gb.GB_l, gb.GB_d);

    } else if (atype->isLennardJones()) {
      gbCut = 2.5 * LJ::Instance()->getSigma(atype);
    }
    
    return gbCut;
  }

 
  void GB::calcForce(AtomType* at1, AtomType* at2, Vector3d d, 
                     RealType r, RealType r2, RealType sw, 
                     RealType vdwMult, RealType &vpair, RealType &pot, 
                     RotMat3x3d A1, RotMat3x3d A2, Vector3d &f1,
                     Vector3d &t1, Vector3d &t2) {

    if (!initialized_) initialize();
    
    pair<AtomType*, AtomType*> key = make_pair(at1, at2);
    GBInteractionData mixer = MixingMap[key];

    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 = A1.getRow(2);
    Vector3d ul2 = A2.getRow(2);

    RealType a, b, g;

    bool i_is_LJ = at1->isLennardJones();
    bool j_is_LJ = at2->isLennardJones();

    if (i_is_LJ) {
      a = 0.0;
      ul1 = V3Zero;
    } else {
      a = dot(d, ul1);
    }

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

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

    RealType au = a / r;
    RealType bu = b / r;
    
    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 / (r - 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 = vdwMult * 4.0 * eps * (R12 - R6);

    RealType s3 = sigma*sigma*sigma;
    RealType s03 = sigma0*sigma0*sigma0;

    RealType pref1 = - vdwMult * 8.0 * eps * mu_ * (R12 - R6) / (e2 * r);

    RealType pref2 = vdwMult * 8.0 * eps * s3 * (6.0*R13 - 3.0*R7) /(dw*r*s03);

    RealType dUdr = - (pref1 * Hp + pref2 * (sigma0*sigma0*r/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 = d / r;   
    Vector3d rxu1 = cross(d, ul1);
    Vector3d rxu2 = cross(d, ul2);
    Vector3d uxu = cross(ul1, ul2);
    
    pot += U*sw;
    f1 += dUdr * rhat + dUda * ul1 + dUdb * ul2;    
    t1 += dUda * rxu1 - dUdg * uxu;
    t2 += dUdb * rxu2 - dUdg * uxu;
    vpair += U*sw;

    return;

  }

  void GB::do_gb_pair(int *atid1, int *atid2, RealType *d, RealType *r, 
                      RealType *r2, RealType *sw, RealType *vdwMult,
                      RealType *vpair, RealType *pot, RealType *A1,
                      RealType *A2, RealType *f1, RealType *t1, RealType *t2) {
    
    if (!initialized_) initialize();
    
    AtomType* atype1 = GBMap[*atid1];
    AtomType* atype2 = GBMap[*atid2];
    
    Vector3d disp(d);
    Vector3d frc(f1);
    Vector3d trq1(t1);
    Vector3d trq2(t2);
    RotMat3x3d Ai(A1);
    RotMat3x3d Aj(A2);
   
    // Fortran has the opposite matrix ordering from c++, so we'll use 
    // transpose here.  When we finish the conversion to C++, this wrapper 
    // will disappear, as will the transpose below:

    calcForce(atype1, atype2, disp, *r, *r2, *sw, *vdwMult, *vpair, *pot, 
              Ai, Aj, frc, trq1, trq1);
      
    f1[0] = frc.x();
    f1[1] = frc.y();
    f1[2] = frc.z();

    t1[0] = trq1.x();
    t1[1] = trq1.y();
    t1[2] = trq1.z();

    t2[0] = trq2.x();
    t2[1] = trq2.y();
    t2[2] = trq2.z();

    return;    
  }
}

extern "C" {
  
#define fortranGetGayBerneCut FC_FUNC(getgaybernecut, GETGAYBERNECUT)
#define fortranDoGBPair FC_FUNC(do_gb_pair, DO_GB_PAIR)
  
  RealType fortranGetGayBerneCut(int* atid) {
    return OpenMD::GB::Instance()->getGayBerneCut(*atid);
  }

  void fortranDoGBPair(int *atid1, int *atid2, RealType *d, RealType *r, 
                       RealType *r2, RealType *sw, RealType *vdwMult, 
                       RealType *vpair, RealType *pot, RealType *A1, 
                       RealType *A2, RealType *f1, RealType *t1, RealType *t2){
    
    return OpenMD::GB::Instance()->do_gb_pair(atid1, atid2, d, r, r2, sw,
                                              vdwMult, vpair, pot, A1, A2, f1,
                                              t1, t2);
  }
}
Revision:	1485
Committed:	Wed Jul 28 19:52:00 2010 UTC (14 years, 10 months ago) by gezelter
File size:	13724 byte(s)
Log Message:	Converting Sticky over to C++
#	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			* [4] Vardeman & Gezelter, in progress (2009).
40			*/
41
42			#include <stdio.h>
43			#include <string.h>
44
45			#include <cmath>
46			#include "nonbonded/GB.hpp"
47			#include "nonbonded/LJ.hpp"
48			#include "utils/simError.h"
49
50			using namespace std;
51			namespace OpenMD {
52
53			bool GB::initialized_ = false;
54			RealType GB::mu_ = 2.0;
55			RealType GB::nu_ = 1.0;
56			ForceField* GB::forceField_ = NULL;
57			map<int, AtomType*> GB::GBMap;
58			map<pair<AtomType, AtomType>, GBInteractionData> GB::MixingMap;
59
60			GB* GB::_instance = NULL;
61
62			GB* GB::Instance() {
63			if (!_instance) {
64			_instance = new GB();
65			}
66			return _instance;
67			}
68
69			GayBerneParam GB::getGayBerneParam(AtomType* atomType) {
70
71			// Do sanity checking on the AtomType we were passed before
72			// building any data structures:
73			if (!atomType->isGayBerne()) {
74			sprintf( painCave.errMsg,
75			"GB::getGayBerneParam was passed an atomType (%s) that does\n"
76			"\tnot appear to be a Gay-Berne atom.\n",
77			atomType->getName().c_str());
78			painCave.severity = OPENMD_ERROR;
79			painCave.isFatal = 1;
80			simError();
81			}
82
83			DirectionalAtomType* daType = dynamic_cast<DirectionalAtomType*>(atomType);
84			GenericData* data = daType->getPropertyByName("GayBerne");
85			if (data == NULL) {
86			sprintf( painCave.errMsg, "GB::getGayBerneParam could not find\n"
87			"\tGay-Berne parameters for atomType %s.\n",
88			daType->getName().c_str());
89			painCave.severity = OPENMD_ERROR;
90			painCave.isFatal = 1;
91			simError();
92			}
93
94			GayBerneParamGenericData* gbData = dynamic_cast<GayBerneParamGenericData*>(data);
95			if (gbData == NULL) {
96			sprintf( painCave.errMsg,
97			"GB::getGayBerneParam could not convert GenericData to\n"
98			"\tGayBerneParamGenericData for atom type %s\n",
99			daType->getName().c_str());
100			painCave.severity = OPENMD_ERROR;
101			painCave.isFatal = 1;
102			simError();
103			}
104
105			return gbData->getData();
106			}
107
108			void GB::initialize() {
109	gezelter	1485
110			ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
111			mu_ = fopts.getGayBerneMu();
112			nu_ = fopts.getGayBerneNu();
113	gezelter	1483	ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
114			ForceField::AtomTypeContainer::MapTypeIterator i;
115			AtomType* at;
116
117			// GB handles all of the GB-GB interactions as well as GB-LJ cross
118			// interactions:
119
120			for (at = atomTypes->beginType(i); at != NULL;
121			at = atomTypes->nextType(i)) {
122
123			if (at->isGayBerne() \|\| at->isLennardJones())
124			addType(at);
125			}
126
127			initialized_ = true;
128			}
129
130			void GB::addType(AtomType* atomType){
131			// add it to the map:
132			AtomTypeProperties atp = atomType->getATP();
133
134			pair<map<int,AtomType*>::iterator,bool> ret;
135			ret = GBMap.insert( pair<int, AtomType*>(atp.ident, atomType) );
136			if (ret.second == false) {
137			sprintf( painCave.errMsg,
138			"GB already had a previous entry with ident %d\n",
139			atp.ident);
140			painCave.severity = OPENMD_INFO;
141			painCave.isFatal = 0;
142			simError();
143			}
144
145			RealType d1, l1, e1, er1, dw1;
146
147			if (atomType->isGayBerne()) {
148			GayBerneParam gb1 = getGayBerneParam(atomType);
149			d1 = gb1.GB_d;
150			l1 = gb1.GB_l;
151			e1 = gb1.GB_eps;
152			er1 = gb1.GB_eps_ratio;
153			dw1 = gb1.GB_dw;
154			} else if (atomType->isLennardJones()) {
155			d1 = LJ::Instance()->getSigma(atomType) / sqrt(2.0);
156			e1 = LJ::Instance()->getEpsilon(atomType);
157			l1 = d1;
158			er1 = 1.0;
159			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
171			// Now, iterate over all known types and add to the mixing map:
172
173			map<int, AtomType*>::iterator it;
174			for( it = GBMap.begin(); it != GBMap.end(); ++it) {
175
176			AtomType* atype2 = (*it).second;
177
178			RealType d2, l2, e2, er2, dw2;
179
180			if (atype2->isGayBerne()) {
181			GayBerneParam gb2 = getGayBerneParam(atype2);
182			d2 = gb2.GB_d;
183			l2 = gb2.GB_l;
184			e2 = gb2.GB_eps;
185			er2 = gb2.GB_eps_ratio;
186			dw2 = gb2.GB_dw;
187			} else if (atype2->isLennardJones()) {
188			d2 = LJ::Instance()->getSigma(atype2) / sqrt(2.0);
189			e2 = LJ::Instance()->getEpsilon(atype2);
190			l2 = d2;
191			er2 = 1.0;
192			dw2 = 1.0;
193			}
194
195			GBInteractionData mixer;
196
197			// Cleaver paper uses sqrt of squares to get sigma0 for
198			// mixed interactions.
199
200			mixer.sigma0 = sqrt(d1d1 + d2d2);
201			mixer.xa2 = (l1l1 - d1d1)/(l1l1 + d2d2);
202			mixer.xai2 = (l2l2 - d2d2)/(l2l2 + d1d1);
203			mixer.x2 = (l1l1 - d1d1) * (l2l2 - d2d2) /
204			((l2l2 + d1d1) * (l1l1 + d2d2));
205
206			// assumed LB mixing rules for now:
207
208			mixer.dw = 0.5 * (dw1 + dw2);
209			mixer.eps0 = sqrt(e1 * e2);
210
211			RealType er = sqrt(er1 * er2);
212			RealType ermu = pow(er,(1.0 / mu_));
213			RealType xp = (1.0 - ermu) / (1.0 + ermu);
214			RealType ap2 = 1.0 / (1.0 + ermu);
215
216			mixer.xp2 = xp * xp;
217			mixer.xpap2 = xp * ap2;
218			mixer.xpapi2 = xp / ap2;
219
220			// only add this pairing if at least one of the atoms is a Gay-Berne atom
221
222			if (atomType->isGayBerne() \|\| atype2->isGayBerne()) {
223
224			pair<AtomType, AtomType> key1, key2;
225			key1 = make_pair(atomType, atype2);
226			key2 = make_pair(atype2, atomType);
227
228			MixingMap[key1] = mixer;
229			if (key2 != key1) {
230			MixingMap[key2] = mixer;
231			}
232			}
233			}
234			}
235
236
237			RealType GB::getGayBerneCut(int atid) {
238			if (!initialized_) initialize();
239			std::map<int, AtomType*> :: const_iterator it;
240			it = GBMap.find(atid);
241			if (it == GBMap.end()) {
242			sprintf( painCave.errMsg,
243			"GB::getGayBerneCut could not find atid %d in GBMap\n",
244			(atid));
245			painCave.severity = OPENMD_ERROR;
246			painCave.isFatal = 1;
247			simError();
248			}
249
250			AtomType* atype = it->second;
251
252			RealType gbCut;
253
254			if (atype->isGayBerne()) {
255			GayBerneParam gb = getGayBerneParam(atype);
256
257			// sigma is actually sqrt(2) * l for prolate ellipsoids
258			gbCut = 2.5 * sqrt(2.0) * max(gb.GB_l, gb.GB_d);
259
260			} else if (atype->isLennardJones()) {
261			gbCut = 2.5 * LJ::Instance()->getSigma(atype);
262			}
263
264			return gbCut;
265			}
266
267
268			void GB::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
269			RealType r, RealType r2, RealType sw,
270			RealType vdwMult, RealType &vpair, RealType &pot,
271			RotMat3x3d A1, RotMat3x3d A2, Vector3d &f1,
272			Vector3d &t1, Vector3d &t2) {
273
274			if (!initialized_) initialize();
275
276			pair<AtomType, AtomType> key = make_pair(at1, at2);
277			GBInteractionData mixer = MixingMap[key];
278
279			RealType sigma0 = mixer.sigma0;
280			RealType dw = mixer.dw;
281			RealType eps0 = mixer.eps0;
282			RealType x2 = mixer.x2;
283			RealType xa2 = mixer.xa2;
284			RealType xai2 = mixer.xai2;
285			RealType xp2 = mixer.xp2;
286			RealType xpap2 = mixer.xpap2;
287			RealType xpapi2 = mixer.xpapi2;
288
289	gezelter	1485	Vector3d ul1 = A1.getRow(2);
290			Vector3d ul2 = A2.getRow(2);
291	gezelter	1483
292			RealType a, b, g;
293
294			bool i_is_LJ = at1->isLennardJones();
295			bool j_is_LJ = at2->isLennardJones();
296
297			if (i_is_LJ) {
298			a = 0.0;
299			ul1 = V3Zero;
300			} else {
301			a = dot(d, ul1);
302			}
303
304			if (j_is_LJ) {
305			b = 0.0;
306			ul2 = V3Zero;
307			} else {
308			b = dot(d, ul2);
309			}
310
311			if (i_is_LJ \|\| j_is_LJ)
312			g = 0.0;
313			else
314			g = dot(ul1, ul2);
315
316			RealType au = a / r;
317			RealType bu = b / r;
318
319			RealType au2 = au * au;
320			RealType bu2 = bu * bu;
321			RealType g2 = g * g;
322
323			RealType H = (xa2 * au2 + xai2 * bu2 - 2.0x2aubug) / (1.0 - x2*g2);
324			RealType Hp = (xpap2au2 + xpapi2bu2 - 2.0xp2aubug) / (1.0 - xp2*g2);
325
326			RealType sigma = sigma0 / sqrt(1.0 - H);
327			RealType e1 = 1.0 / sqrt(1.0 - x2*g2);
328			RealType e2 = 1.0 - Hp;
329			RealType eps = eps0 * pow(e1,nu_) * pow(e2,mu_);
330			RealType BigR = dwsigma0 / (r - sigma + dwsigma0);
331
332			RealType R3 = BigRBigRBigR;
333			RealType R6 = R3*R3;
334			RealType R7 = R6 * BigR;
335			RealType R12 = R6*R6;
336			RealType R13 = R6*R7;
337
338			RealType U = vdwMult * 4.0 * eps * (R12 - R6);
339
340			RealType s3 = sigmasigmasigma;
341			RealType s03 = sigma0sigma0sigma0;
342
343			RealType pref1 = - vdwMult * 8.0 * eps * mu_ * (R12 - R6) / (e2 * r);
344
345			RealType pref2 = vdwMult * 8.0 * eps * s3 * (6.0R13 - 3.0R7) /(dwrs03);
346
347			RealType dUdr = - (pref1 * Hp + pref2 * (sigma0sigma0r/s3 + H));
348
349			RealType dUda = pref1 * (xpap2au - xp2bug) / (1.0 - xp2 g2)
350			+ pref2 * (xa2 * au - x2 bug) / (1.0 - x2 * g2);
351
352			RealType dUdb = pref1 * (xpapi2bu - xp2aug) / (1.0 - xp2 g2)
353			+ pref2 * (xai2 * bu - x2 aug) / (1.0 - x2 * g2);
354
355			RealType dUdg = 4.0 * eps * nu_ * (R12 - R6) * x2 * g / (1.0 - x2*g2)
356			+ 8.0 * eps * mu_ * (R12 - R6) * (xp2aubu - Hpxp2g) /
357			(1.0 - xp2 * g2) / e2 + 8.0 * eps * s3 * (3.0 * R7 - 6.0 * R13) *
358			(x2 * au * bu - H * x2 * g) / (1.0 - x2 * g2) / (dw * s03);
359
360
361			Vector3d rhat = d / r;
362			Vector3d rxu1 = cross(d, ul1);
363			Vector3d rxu2 = cross(d, ul2);
364			Vector3d uxu = cross(ul1, ul2);
365
366			pot += U*sw;
367			f1 += dUdr * rhat + dUda * ul1 + dUdb * ul2;
368			t1 += dUda * rxu1 - dUdg * uxu;
369			t2 += dUdb * rxu2 - dUdg * uxu;
370			vpair += U*sw;
371
372			return;
373
374			}
375
376			void GB::do_gb_pair(int atid1, int atid2, RealType d, RealType r,
377			RealType r2, RealType sw, RealType *vdwMult,
378			RealType vpair, RealType pot, RealType *A1,
379			RealType A2, RealType f1, RealType t1, RealType t2) {
380
381			if (!initialized_) initialize();
382
383			AtomType* atype1 = GBMap[*atid1];
384			AtomType* atype2 = GBMap[*atid2];
385
386			Vector3d disp(d);
387			Vector3d frc(f1);
388			Vector3d trq1(t1);
389			Vector3d trq2(t2);
390			RotMat3x3d Ai(A1);
391			RotMat3x3d Aj(A2);
392
393			// Fortran has the opposite matrix ordering from c++, so we'll use
394			// transpose here. When we finish the conversion to C++, this wrapper
395			// will disappear, as will the transpose below:
396
397			calcForce(atype1, atype2, disp, r, r2, sw, vdwMult, vpair, pot,
398	gezelter	1485	Ai, Aj, frc, trq1, trq1);
399	gezelter	1483
400			f1[0] = frc.x();
401			f1[1] = frc.y();
402			f1[2] = frc.z();
403
404			t1[0] = trq1.x();
405			t1[1] = trq1.y();
406			t1[2] = trq1.z();
407
408			t2[0] = trq2.x();
409			t2[1] = trq2.y();
410			t2[2] = trq2.z();
411
412			return;
413			}
414			}
415
416			extern "C" {
417
418			#define fortranGetGayBerneCut FC_FUNC(getgaybernecut, GETGAYBERNECUT)
419			#define fortranDoGBPair FC_FUNC(do_gb_pair, DO_GB_PAIR)
420
421			RealType fortranGetGayBerneCut(int* atid) {
422			return OpenMD::GB::Instance()->getGayBerneCut(*atid);
423			}
424
425			void fortranDoGBPair(int atid1, int atid2, RealType d, RealType r,
426			RealType r2, RealType sw, RealType *vdwMult,
427			RealType vpair, RealType pot, RealType *A1,
428			RealType A2, RealType f1, RealType t1, RealType t2){
429
430			return OpenMD::GB::Instance()->do_gb_pair(atid1, atid2, d, r, r2, sw,
431			vdwMult, vpair, pot, A1, A2, f1,
432			t1, t2);
433			}
434			}