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 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, (RealType(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;

      cerr << "mixer" << er1 << " " << er2 << " " << mu_ << " " << ermu << " " << xp <<" " << ap2 << "\n";

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

    cerr << "in GB, d = " << *(idat.d) << "\n";
    cerr << "abg = " << a << " " << b << " " << g <<"\n";

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

    RealType sigma = sigma0 / sqrt(1.0 - H);
    RealType e1 = 1.0 / sqrt(1.0 - x2*g2);
    RealType e2 = 1.0 - Hp;
    RealType eps = eps0 * pow(e1,nu_) * pow(e2,mu_);
    cerr << "eps = " << eps0 << " " << e1 << " " << nu_ << " " << e2 << " " << mu_ << "\n";
    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);

    cerr << "R12, R6, eps = " << R12 << " " << R6 << " " << eps << " " <<  *(idat.vdwMult) << "\n";

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

    cerr << pref1 << " " << pref2 << " " << dUdr <<" " << dUda << " " << dUdb << dUdg << "\n";

    Vector3d rhat = *(idat.d) / *(idat.rij);   
    Vector3d rxu1 = cross(*(idat.d), ul1);
    Vector3d rxu2 = cross(*(idat.d), ul2);
    Vector3d uxu = cross(ul1, ul2);

    cerr << "U = " << U << "\n";
    cerr << "f1 = " << dUdr * rhat + dUda * ul1 + dUdb * ul2 << "\n";
    cerr << "t1 = " << dUda * rxu1 - dUdg * uxu << "\n";
    cerr << "t2 = " << dUdb * rxu2 - dUdg * uxu << "\n";

    
    (*(idat.pot))[VANDERWAALS_FAMILY] += U *  *(idat.sw);
    *(idat.f1) += dUdr * rhat + dUda * ul1 + dUdb * ul2;    
    *(idat.t1) += dUda * rxu1 - dUdg * uxu;
    *(idat.t2) += dUdb * rxu2 - dUdg * uxu;
    *(idat.vpair) += U * *(idat.sw);

    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:	1668
Committed:	Fri Jan 6 19:03:05 2012 UTC (13 years, 6 months ago) by gezelter
File size:	14139 byte(s)
Log Message:	Some fixes for CMake and single precision builds
#	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			GBInteractionData mixer;
227
228			// Cleaver paper uses sqrt of squares to get sigma0 for
229			// mixed interactions.
230
231			mixer.sigma0 = sqrt(d1d1 + d2d2);
232			mixer.xa2 = (l1l1 - d1d1)/(l1l1 + d2d2);
233			mixer.xai2 = (l2l2 - d2d2)/(l2l2 + d1d1);
234			mixer.x2 = (l1l1 - d1d1) * (l2l2 - d2d2) /
235			((l2l2 + d1d1) * (l1l1 + d2d2));
236
237			// assumed LB mixing rules for now:
238
239			mixer.dw = 0.5 * (dw1 + dw2);
240			mixer.eps0 = sqrt(e1 * e2);
241
242			RealType er = sqrt(er1 * er2);
243	gezelter	1668	RealType ermu = pow(er, (RealType(1.0) / mu_));
244	gezelter	1483	RealType xp = (1.0 - ermu) / (1.0 + ermu);
245			RealType ap2 = 1.0 / (1.0 + ermu);
246
247			mixer.xp2 = xp * xp;
248			mixer.xpap2 = xp * ap2;
249			mixer.xpapi2 = xp / ap2;
250
251	gezelter	1587	cerr << "mixer" << er1 << " " << er2 << " " << mu_ << " " << ermu << " " << xp <<" " << ap2 << "\n";
252
253	gezelter	1483	// only add this pairing if at least one of the atoms is a Gay-Berne atom
254
255			if (atomType->isGayBerne() \|\| atype2->isGayBerne()) {
256
257			pair<AtomType, AtomType> key1, key2;
258			key1 = make_pair(atomType, atype2);
259			key2 = make_pair(atype2, atomType);
260
261			MixingMap[key1] = mixer;
262			if (key2 != key1) {
263			MixingMap[key2] = mixer;
264			}
265			}
266			}
267			}
268	gezelter	1502
269	gezelter	1536	void GB::calcForce(InteractionData &idat) {
270	gezelter	1483
271			if (!initialized_) initialize();
272
273	gezelter	1571	GBInteractionData mixer = MixingMap[idat.atypes];
274	gezelter	1483
275			RealType sigma0 = mixer.sigma0;
276			RealType dw = mixer.dw;
277			RealType eps0 = mixer.eps0;
278			RealType x2 = mixer.x2;
279			RealType xa2 = mixer.xa2;
280			RealType xai2 = mixer.xai2;
281			RealType xp2 = mixer.xp2;
282			RealType xpap2 = mixer.xpap2;
283			RealType xpapi2 = mixer.xpapi2;
284
285	gezelter	1554	Vector3d ul1 = idat.A1->getRow(2);
286			Vector3d ul2 = idat.A2->getRow(2);
287	gezelter	1483
288			RealType a, b, g;
289
290	gezelter	1571	bool i_is_LJ = idat.atypes.first->isLennardJones();
291			bool j_is_LJ = idat.atypes.second->isLennardJones();
292	gezelter	1554
293	gezelter	1483	if (i_is_LJ) {
294			a = 0.0;
295			ul1 = V3Zero;
296			} else {
297	gezelter	1554	a = dot(*(idat.d), ul1);
298	gezelter	1483	}
299
300			if (j_is_LJ) {
301			b = 0.0;
302			ul2 = V3Zero;
303			} else {
304	gezelter	1554	b = dot(*(idat.d), ul2);
305	gezelter	1483	}
306
307			if (i_is_LJ \|\| j_is_LJ)
308			g = 0.0;
309			else
310			g = dot(ul1, ul2);
311
312	gezelter	1587	cerr << "in GB, d = " << *(idat.d) << "\n";
313			cerr << "abg = " << a << " " << b << " " << g <<"\n";
314
315	gezelter	1554	RealType au = a / *(idat.rij);
316			RealType bu = b / *(idat.rij);
317	gezelter	1483
318			RealType au2 = au * au;
319			RealType bu2 = bu * bu;
320			RealType g2 = g * g;
321
322			RealType H = (xa2 * au2 + xai2 * bu2 - 2.0x2aubug) / (1.0 - x2*g2);
323			RealType Hp = (xpap2au2 + xpapi2bu2 - 2.0xp2aubug) / (1.0 - xp2*g2);
324	gezelter	1587	cerr << "xa2, xai2 " << xa2 << " " << xai2 << "\n";
325			cerr << "xpap2, xpapi2 " << xpap2 << " " << xpapi2 << "\n";
326			cerr << "H Hp = " << H << " " << Hp << "\n";
327	gezelter	1483
328			RealType sigma = sigma0 / sqrt(1.0 - H);
329			RealType e1 = 1.0 / sqrt(1.0 - x2*g2);
330			RealType e2 = 1.0 - Hp;
331			RealType eps = eps0 * pow(e1,nu_) * pow(e2,mu_);
332	gezelter	1587	cerr << "eps = " << eps0 << " " << e1 << " " << nu_ << " " << e2 << " " << mu_ << "\n";
333	gezelter	1554	RealType BigR = dwsigma0 / ((idat.rij) - sigma + dw*sigma0);
334	gezelter	1483
335			RealType R3 = BigRBigRBigR;
336			RealType R6 = R3*R3;
337			RealType R7 = R6 * BigR;
338			RealType R12 = R6*R6;
339			RealType R13 = R6*R7;
340
341	gezelter	1554	RealType U = (idat.vdwMult) 4.0 * eps * (R12 - R6);
342	gezelter	1483
343	gezelter	1587	cerr << "R12, R6, eps = " << R12 << " " << R6 << " " << eps << " " << *(idat.vdwMult) << "\n";
344
345	gezelter	1483	RealType s3 = sigmasigmasigma;
346			RealType s03 = sigma0sigma0sigma0;
347
348	gezelter	1554	RealType pref1 = - (idat.vdwMult) 8.0 * eps * mu_ * (R12 - R6) /
349			(e2 * *(idat.rij));
350	gezelter	1483
351	gezelter	1554	RealType pref2 = (idat.vdwMult) 8.0 * eps * s3 * (6.0R13 - 3.0R7) /
352			(dw* (idat.rij) s03);
353	gezelter	1483
354	gezelter	1554	RealType dUdr = - (pref1 * Hp + pref2 * (sigma0 * sigma0 *
355			*(idat.rij) / s3 + H));
356	gezelter	1483
357			RealType dUda = pref1 * (xpap2au - xp2bug) / (1.0 - xp2 g2)
358			+ pref2 * (xa2 * au - x2 bug) / (1.0 - x2 * g2);
359
360			RealType dUdb = pref1 * (xpapi2bu - xp2aug) / (1.0 - xp2 g2)
361			+ pref2 * (xai2 * bu - x2 aug) / (1.0 - x2 * g2);
362
363			RealType dUdg = 4.0 * eps * nu_ * (R12 - R6) * x2 * g / (1.0 - x2*g2)
364			+ 8.0 * eps * mu_ * (R12 - R6) * (xp2aubu - Hpxp2g) /
365			(1.0 - xp2 * g2) / e2 + 8.0 * eps * s3 * (3.0 * R7 - 6.0 * R13) *
366			(x2 * au * bu - H * x2 * g) / (1.0 - x2 * g2) / (dw * s03);
367
368	gezelter	1587	cerr << pref1 << " " << pref2 << " " << dUdr <<" " << dUda << " " << dUdb << dUdg << "\n";
369
370	gezelter	1554	Vector3d rhat = (idat.d) / (idat.rij);
371			Vector3d rxu1 = cross(*(idat.d), ul1);
372			Vector3d rxu2 = cross(*(idat.d), ul2);
373	gezelter	1483	Vector3d uxu = cross(ul1, ul2);
374	gezelter	1587
375			cerr << "U = " << U << "\n";
376			cerr << "f1 = " << dUdr * rhat + dUda * ul1 + dUdb * ul2 << "\n";
377			cerr << "t1 = " << dUda * rxu1 - dUdg * uxu << "\n";
378			cerr << "t2 = " << dUdb * rxu2 - dUdg * uxu << "\n";
379
380	gezelter	1483
381	gezelter	1582	((idat.pot))[VANDERWAALS_FAMILY] += U *(idat.sw);
382	gezelter	1554	(idat.f1) += dUdr rhat + dUda * ul1 + dUdb * ul2;
383			(idat.t1) += dUda rxu1 - dUdg * uxu;
384			(idat.t2) += dUdb rxu2 - dUdg * uxu;
385			(idat.vpair) += U *(idat.sw);
386	gezelter	1483
387			return;
388
389			}
390	gezelter	1505
391	gezelter	1545	RealType GB::getSuggestedCutoffRadius(pair<AtomType, AtomType> atypes) {
392	gezelter	1505	if (!initialized_) initialize();
393
394			RealType cut = 0.0;
395
396	gezelter	1545	if (atypes.first->isGayBerne()) {
397			GayBerneParam gb1 = getGayBerneParam(atypes.first);
398	gezelter	1505	RealType d1 = gb1.GB_d;
399			RealType l1 = gb1.GB_l;
400			// sigma is actually sqrt(2)*l for prolate ellipsoids
401	gezelter	1668	cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d1, l1));
402	gezelter	1545	} else if (atypes.first->isLennardJones()) {
403	gezelter	1668	cut = max(cut, RealType(2.5) * getLJSigma(atypes.first));
404	gezelter	1505	}
405
406	gezelter	1545	if (atypes.second->isGayBerne()) {
407			GayBerneParam gb2 = getGayBerneParam(atypes.second);
408	gezelter	1505	RealType d2 = gb2.GB_d;
409			RealType l2 = gb2.GB_l;
410	gezelter	1668	cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d2, l2));
411	gezelter	1545	} else if (atypes.second->isLennardJones()) {
412	gezelter	1668	cut = max(cut, RealType(2.5) * getLJSigma(atypes.second));
413	gezelter	1505	}
414
415			return cut;
416			}
417	gezelter	1483	}
418