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 "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,(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;
        }
      }
    }      
  }
   
  void GB::calcForce(InteractionData idat) {

    if (!initialized_) initialize();
    
    pair<AtomType*, AtomType*> key = make_pair(idat.atype1, idat.atype2);
    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 = idat.A1.getRow(2);
    Vector3d ul2 = idat.A2.getRow(2);

    RealType a, b, g;

    bool i_is_LJ = idat.atype1->isLennardJones();
    bool j_is_LJ = idat.atype2->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 += 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;

  }
}

Revision:	1502
Committed:	Sat Oct 2 19:53:32 2010 UTC (14 years, 7 months ago) by gezelter
File size:	12125 byte(s)
Log Message:	Many changes, and we're not quite done with this phase yet.
#	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 "utils/simError.h"
48
49			using namespace std;
50			namespace OpenMD {
51
52	gezelter	1502	GB::GB() : name_("GB"), initialized_(false), mu_(2.0), nu_(1.0), forceField_(NULL) {}
53	gezelter	1483
54			GayBerneParam GB::getGayBerneParam(AtomType* atomType) {
55
56			// Do sanity checking on the AtomType we were passed before
57			// building any data structures:
58			if (!atomType->isGayBerne()) {
59			sprintf( painCave.errMsg,
60			"GB::getGayBerneParam was passed an atomType (%s) that does\n"
61			"\tnot appear to be a Gay-Berne atom.\n",
62			atomType->getName().c_str());
63			painCave.severity = OPENMD_ERROR;
64			painCave.isFatal = 1;
65			simError();
66			}
67
68			DirectionalAtomType* daType = dynamic_cast<DirectionalAtomType*>(atomType);
69			GenericData* data = daType->getPropertyByName("GayBerne");
70			if (data == NULL) {
71			sprintf( painCave.errMsg, "GB::getGayBerneParam could not find\n"
72			"\tGay-Berne parameters for atomType %s.\n",
73			daType->getName().c_str());
74			painCave.severity = OPENMD_ERROR;
75			painCave.isFatal = 1;
76			simError();
77			}
78
79			GayBerneParamGenericData* gbData = dynamic_cast<GayBerneParamGenericData*>(data);
80			if (gbData == NULL) {
81			sprintf( painCave.errMsg,
82			"GB::getGayBerneParam could not convert GenericData to\n"
83			"\tGayBerneParamGenericData for atom type %s\n",
84			daType->getName().c_str());
85			painCave.severity = OPENMD_ERROR;
86			painCave.isFatal = 1;
87			simError();
88			}
89
90			return gbData->getData();
91			}
92
93	gezelter	1502	LJParam GB::getLJParam(AtomType* atomType) {
94
95			// Do sanity checking on the AtomType we were passed before
96			// building any data structures:
97			if (!atomType->isLennardJones()) {
98			sprintf( painCave.errMsg,
99			"GB::getLJParam was passed an atomType (%s) that does not\n"
100			"\tappear to be a Lennard-Jones atom.\n",
101			atomType->getName().c_str());
102			painCave.severity = OPENMD_ERROR;
103			painCave.isFatal = 1;
104			simError();
105			}
106
107			GenericData* data = atomType->getPropertyByName("LennardJones");
108			if (data == NULL) {
109			sprintf( painCave.errMsg, "GB::getLJParam could not find Lennard-Jones\n"
110			"\tparameters for atomType %s.\n", atomType->getName().c_str());
111			painCave.severity = OPENMD_ERROR;
112			painCave.isFatal = 1;
113			simError();
114			}
115
116			LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
117			if (ljData == NULL) {
118			sprintf( painCave.errMsg,
119			"GB::getLJParam could not convert GenericData to LJParam for\n"
120			"\tatom type %s\n", atomType->getName().c_str());
121			painCave.severity = OPENMD_ERROR;
122			painCave.isFatal = 1;
123			simError();
124			}
125
126			return ljData->getData();
127			}
128
129			RealType GB::getLJEpsilon(AtomType* atomType) {
130			LJParam ljParam = getLJParam(atomType);
131			return ljParam.epsilon;
132			}
133			RealType GB::getLJSigma(AtomType* atomType) {
134			LJParam ljParam = getLJParam(atomType);
135			return ljParam.sigma;
136			}
137
138	gezelter	1483	void GB::initialize() {
139	gezelter	1502
140	gezelter	1485	ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
141			mu_ = fopts.getGayBerneMu();
142			nu_ = fopts.getGayBerneNu();
143	gezelter	1483	ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
144			ForceField::AtomTypeContainer::MapTypeIterator i;
145			AtomType* at;
146
147			// GB handles all of the GB-GB interactions as well as GB-LJ cross
148			// interactions:
149
150			for (at = atomTypes->beginType(i); at != NULL;
151			at = atomTypes->nextType(i)) {
152
153			if (at->isGayBerne() \|\| at->isLennardJones())
154			addType(at);
155			}
156
157			initialized_ = true;
158			}
159
160			void GB::addType(AtomType* atomType){
161			// add it to the map:
162			AtomTypeProperties atp = atomType->getATP();
163
164			pair<map<int,AtomType*>::iterator,bool> ret;
165			ret = GBMap.insert( pair<int, AtomType*>(atp.ident, atomType) );
166			if (ret.second == false) {
167			sprintf( painCave.errMsg,
168			"GB already had a previous entry with ident %d\n",
169			atp.ident);
170			painCave.severity = OPENMD_INFO;
171			painCave.isFatal = 0;
172			simError();
173			}
174
175			RealType d1, l1, e1, er1, dw1;
176
177			if (atomType->isGayBerne()) {
178			GayBerneParam gb1 = getGayBerneParam(atomType);
179			d1 = gb1.GB_d;
180			l1 = gb1.GB_l;
181			e1 = gb1.GB_eps;
182			er1 = gb1.GB_eps_ratio;
183			dw1 = gb1.GB_dw;
184			} else if (atomType->isLennardJones()) {
185	gezelter	1502	d1 = getLJSigma(atomType) / sqrt(2.0);
186			e1 = getLJEpsilon(atomType);
187	gezelter	1483	l1 = d1;
188			er1 = 1.0;
189			dw1 = 1.0;
190			} else {
191			sprintf( painCave.errMsg,
192			"GB::addType was passed an atomType (%s) that does not\n"
193			"\tappear to be a Gay-Berne or Lennard-Jones atom.\n",
194			atomType->getName().c_str());
195			painCave.severity = OPENMD_ERROR;
196			painCave.isFatal = 1;
197			simError();
198			}
199
200
201			// Now, iterate over all known types and add to the mixing map:
202
203			map<int, AtomType*>::iterator it;
204			for( it = GBMap.begin(); it != GBMap.end(); ++it) {
205
206			AtomType* atype2 = (*it).second;
207
208			RealType d2, l2, e2, er2, dw2;
209
210			if (atype2->isGayBerne()) {
211			GayBerneParam gb2 = getGayBerneParam(atype2);
212			d2 = gb2.GB_d;
213			l2 = gb2.GB_l;
214			e2 = gb2.GB_eps;
215			er2 = gb2.GB_eps_ratio;
216			dw2 = gb2.GB_dw;
217			} else if (atype2->isLennardJones()) {
218	gezelter	1502	d2 = getLJSigma(atype2) / sqrt(2.0);
219			e2 = getLJEpsilon(atype2);
220	gezelter	1483	l2 = d2;
221			er2 = 1.0;
222			dw2 = 1.0;
223			}
224
225			GBInteractionData mixer;
226
227			// Cleaver paper uses sqrt of squares to get sigma0 for
228			// mixed interactions.
229
230			mixer.sigma0 = sqrt(d1d1 + d2d2);
231			mixer.xa2 = (l1l1 - d1d1)/(l1l1 + d2d2);
232			mixer.xai2 = (l2l2 - d2d2)/(l2l2 + d1d1);
233			mixer.x2 = (l1l1 - d1d1) * (l2l2 - d2d2) /
234			((l2l2 + d1d1) * (l1l1 + d2d2));
235
236			// assumed LB mixing rules for now:
237
238			mixer.dw = 0.5 * (dw1 + dw2);
239			mixer.eps0 = sqrt(e1 * e2);
240
241			RealType er = sqrt(er1 * er2);
242			RealType ermu = pow(er,(1.0 / mu_));
243			RealType xp = (1.0 - ermu) / (1.0 + ermu);
244			RealType ap2 = 1.0 / (1.0 + ermu);
245
246			mixer.xp2 = xp * xp;
247			mixer.xpap2 = xp * ap2;
248			mixer.xpapi2 = xp / ap2;
249
250			// only add this pairing if at least one of the atoms is a Gay-Berne atom
251
252			if (atomType->isGayBerne() \|\| atype2->isGayBerne()) {
253
254			pair<AtomType, AtomType> key1, key2;
255			key1 = make_pair(atomType, atype2);
256			key2 = make_pair(atype2, atomType);
257
258			MixingMap[key1] = mixer;
259			if (key2 != key1) {
260			MixingMap[key2] = mixer;
261			}
262			}
263			}
264			}
265	gezelter	1502
266			void GB::calcForce(InteractionData idat) {
267	gezelter	1483
268			if (!initialized_) initialize();
269
270	gezelter	1502	pair<AtomType, AtomType> key = make_pair(idat.atype1, idat.atype2);
271	gezelter	1483	GBInteractionData mixer = MixingMap[key];
272
273			RealType sigma0 = mixer.sigma0;
274			RealType dw = mixer.dw;
275			RealType eps0 = mixer.eps0;
276			RealType x2 = mixer.x2;
277			RealType xa2 = mixer.xa2;
278			RealType xai2 = mixer.xai2;
279			RealType xp2 = mixer.xp2;
280			RealType xpap2 = mixer.xpap2;
281			RealType xpapi2 = mixer.xpapi2;
282
283	gezelter	1502	Vector3d ul1 = idat.A1.getRow(2);
284			Vector3d ul2 = idat.A2.getRow(2);
285	gezelter	1483
286			RealType a, b, g;
287
288	gezelter	1502	bool i_is_LJ = idat.atype1->isLennardJones();
289			bool j_is_LJ = idat.atype2->isLennardJones();
290	gezelter	1483
291			if (i_is_LJ) {
292			a = 0.0;
293			ul1 = V3Zero;
294			} else {
295	gezelter	1502	a = dot(idat.d, ul1);
296	gezelter	1483	}
297
298			if (j_is_LJ) {
299			b = 0.0;
300			ul2 = V3Zero;
301			} else {
302	gezelter	1502	b = dot(idat.d, ul2);
303	gezelter	1483	}
304
305			if (i_is_LJ \|\| j_is_LJ)
306			g = 0.0;
307			else
308			g = dot(ul1, ul2);
309
310	gezelter	1502	RealType au = a / idat.rij;
311			RealType bu = b / idat.rij;
312	gezelter	1483
313			RealType au2 = au * au;
314			RealType bu2 = bu * bu;
315			RealType g2 = g * g;
316
317			RealType H = (xa2 * au2 + xai2 * bu2 - 2.0x2aubug) / (1.0 - x2*g2);
318			RealType Hp = (xpap2au2 + xpapi2bu2 - 2.0xp2aubug) / (1.0 - xp2*g2);
319
320			RealType sigma = sigma0 / sqrt(1.0 - H);
321			RealType e1 = 1.0 / sqrt(1.0 - x2*g2);
322			RealType e2 = 1.0 - Hp;
323			RealType eps = eps0 * pow(e1,nu_) * pow(e2,mu_);
324	gezelter	1502	RealType BigR = dwsigma0 / (idat.rij - sigma + dwsigma0);
325	gezelter	1483
326			RealType R3 = BigRBigRBigR;
327			RealType R6 = R3*R3;
328			RealType R7 = R6 * BigR;
329			RealType R12 = R6*R6;
330			RealType R13 = R6*R7;
331
332	gezelter	1502	RealType U = idat.vdwMult * 4.0 * eps * (R12 - R6);
333	gezelter	1483
334			RealType s3 = sigmasigmasigma;
335			RealType s03 = sigma0sigma0sigma0;
336
337	gezelter	1502	RealType pref1 = - idat.vdwMult * 8.0 * eps * mu_ * (R12 - R6) / (e2 * idat.rij);
338	gezelter	1483
339	gezelter	1502	RealType pref2 = idat.vdwMult * 8.0 * eps * s3 * (6.0R13 - 3.0R7) /(dwidat.rijs03);
340	gezelter	1483
341	gezelter	1502	RealType dUdr = - (pref1 * Hp + pref2 * (sigma0sigma0idat.rij/s3 + H));
342	gezelter	1483
343			RealType dUda = pref1 * (xpap2au - xp2bug) / (1.0 - xp2 g2)
344			+ pref2 * (xa2 * au - x2 bug) / (1.0 - x2 * g2);
345
346			RealType dUdb = pref1 * (xpapi2bu - xp2aug) / (1.0 - xp2 g2)
347			+ pref2 * (xai2 * bu - x2 aug) / (1.0 - x2 * g2);
348
349			RealType dUdg = 4.0 * eps * nu_ * (R12 - R6) * x2 * g / (1.0 - x2*g2)
350			+ 8.0 * eps * mu_ * (R12 - R6) * (xp2aubu - Hpxp2g) /
351			(1.0 - xp2 * g2) / e2 + 8.0 * eps * s3 * (3.0 * R7 - 6.0 * R13) *
352			(x2 * au * bu - H * x2 * g) / (1.0 - x2 * g2) / (dw * s03);
353
354
355	gezelter	1502	Vector3d rhat = idat.d / idat.rij;
356			Vector3d rxu1 = cross(idat.d, ul1);
357			Vector3d rxu2 = cross(idat.d, ul2);
358	gezelter	1483	Vector3d uxu = cross(ul1, ul2);
359
360	gezelter	1502	idat.pot += U*idat.sw;
361			idat.f1 += dUdr * rhat + dUda * ul1 + dUdb * ul2;
362			idat.t1 += dUda * rxu1 - dUdg * uxu;
363			idat.t2 += dUdb * rxu2 - dUdg * uxu;
364			idat.vpair += U*idat.sw;
365	gezelter	1483
366			return;
367
368			}
369			}
370