src/nonbonded/SHAPES.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/SHAPES.hpp"
#include "nonbonded/LJ.hpp"
#include "utils/simError.h"

using namespace std;
namespace OpenMD {
 
  SHAPES::SHAPES() {
    initialized_ = false;
    lMax_ = 64;
    mMax_ = 64;
    forceField_ = NULL;
  }
  
  void SHAPES::initialize() {    
    
    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
    ForceField::AtomTypeContainer::MapTypeIterator i;
    AtomType* at;

    // SHAPES handles all of the SHAPES-SHAPES interactions as well as
    // SHAPES-LJ cross interactions:
    
    for (at = atomTypes->beginType(i); at != NULL;
         at = atomTypes->nextType(i)) {
      
      if (at->isShape()) 
        addShape(dynamic_cast<ShapeAtomType*>(at));

      if (at->isLennardJones())
        addLJ(at);
      
    }
    
    initialized_ = true;
  }
  
  void SHAPES::addShape(ShapeAtomType* atomType){
    // add it to the map:
    AtomTypeProperties atp = atomType->getATP();    

    pair<map<int,ShapeAtomType*>::iterator, bool> ret;
    ret = shapesMap.insert( pair<int, ShapeAtomType*>(atp.ident, atomType));
    if (ret.second == false) {
      sprintf( painCave.errmsg,
               "SHAPES already had a previous entry with ident %d\n",
               atp.ident);
      painCave.severity = OPENMD_INFO;
      painCave.isFatal = 0;
      simError();         
    }   
    
    ShapesMap.insert( pair<int, ShapeAtomType*>(atp.ident, sAtomType) );
    
  } else if (atomType->isLennardJones()) {
      d1 = getLJSigma(atomType) / sqrt(2.0);
      e1 = getLJEpsilon(atomType);
    } else {
      sprintf( painCave.errMsg,
               "SHAPES::addType was passed an atomType (%s) that does not\n"
               "\tappear to be a SHAPES or Lennard-Jones atom.\n",
               atomType->getName().c_str());
      painCave.severity = OPENMD_ERROR;
      painCave.isFatal = 1;
      simError();
    }   
  }
  

  LJParam SHAPES::getLJParam(AtomType* atomType) {
    
    // Do sanity checking on the AtomType we were passed before
    // building any data structures:
    if (!atomType->isLennardJones()) {
      sprintf( painCave.errMsg,
               "SHAPES::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, "SHAPES::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,
               "SHAPES::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 SHAPES::getLJEpsilon(AtomType* atomType) {    
    LJParam ljParam = getLJParam(atomType);
    return ljParam.epsilon;
  }
  RealType SHAPES::getLJSigma(AtomType* atomType) {    
    LJParam ljParam = getLJParam(atomType);
    return ljParam.sigma;
  }

  RealType SHAPES::getGayBerneCut(int atid) { 
    if (!initialized_) initialize();
    std::map<int, AtomType*> :: const_iterator it;
    it = SHAPESMap.find(atid);
    if (it == SHAPESMap.end()) {
      sprintf( painCave.errMsg,
               "SHAPES::getGayBerneCut could not find atid %d in SHAPESMap\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.SHAPES_l, gb.SHAPES_d);

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

 
  void SHAPES::calcForce(AtomType* at1, AtomType* at2, Vector3d d, 
                         RealType r, RealType r2, RealType sw, 
                         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);
    SHAPESInteractionData mixer = MixingMap[key];

    RealType r3 = r2 * r;
    RealType r5 = r3 * r2;
    
    Vector3d drdi  = -d / r;
    Vector3d drdui = V3Zero;
    Vector3d drdj = d / r;
    Vector3d drduj = V3Zero;
    
    bool i_is_LJ = at1->isLennardJones();
    bool j_is_LJ = at2->isLennardJones();

    RealType sigma_i;
    RealType s_i;
    RealType eps_i;
    Vector3d dsigmaidr;
    Vector3d disgmaidu;
    Vector3d dsidr;
    Vector3d dsidu;
    Vector3d depsidr;
    Vector3d depsidu;

    if (i_is_LJ) {
      sigma_i = LJ::Instance()->getSigma(at1);
      s_i = sigma_i;
      epsilon_i = LJ::Instance()->getEpsilon(at1);
      dsigmaidr = V3Zero;
      dsigmaidu = V3Zero;
      dsidr = V3Zero;
      dsidu = V3Zero;
      depsidr = V3Zero;
      depsidu = V3Zero;
    } else {
      
      // rotate the inter-particle separation into the two different
      // body-fixed coordinate systems:

      Vector3d ri = A1 * d;
      
      RealType xi = ri.x() / r;
      RealType yi = ri.y() / r;
      RealType zi = ri.z() / r;
      RealType xi2 = xi * xi;
      RealType yi2 = yi * yi;
      RealType zi2 = zi * zi;
      RealType cti = zi / r;

      if (cti > 1.0) cti = 1.0;
      if (cti < -1.0_dp) cti = -1.0;

      Vector3d dctidr(-zi * xi / r3, 
                      -zi * yi / r3, 
                      1.0 / r - zi2 / r3);

      Vector3d dctidu(yi / r, 
                      -zi / r,
                      0.0);

      // this is an attempt to try to truncate the singularity when
      // sin(theta) is near 0.0:

      RealType sti2 = 1.0 - cti*cti;
      if (fabs(sti2) < 1.0e-12) {
        RealType proji = sqrt(r * 1.0e-12);
        Vector3d dcpidx(1.0 / proji,
                        0.0,
                        
        dcpidx = 1.0_dp / proji
          dcpidy = 0.0_dp
          dcpidux = xi / proji
          dcpiduy = 0.0_dp
          dspidx = 0.0_dp
          dspidy = 1.0_dp / proji
          dspidux = 0.0_dp
          dspiduy = yi / proji
       else
          proji = sqrt(xi2 + yi2)
          proji3 = proji*proji*proji
          dcpidx = 1.0_dp / proji - xi2 / proji3
          dcpidy = - xi * yi / proji3
          dcpidux = xi / proji - (xi2 * xi) / proji3
          dcpiduy = - (xi * yi2) / proji3
          dspidx = - xi * yi / proji3
          dspidy = 1.0_dp / proji - yi2 / proji3
          dspidux = - (yi * xi2) / proji3
          dspiduy = yi / proji - (yi2 * yi) / proji3
       endif

       cpi = xi / proji
       dcpidz = 0.0_dp
       dcpiduz = 0.0_dp

       spi = yi / proji
       dspidz = 0.0_dp
       dspiduz = 0.0_dp


    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;


    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 SHAPES::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 = SHAPESMap[*atid1];
    AtomType* atype2 = SHAPESMap[*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 fortranDoSHAPESPair FC_FUNC(do_gb_pair, DO_SHAPES_PAIR)
  
  RealType fortranGetGayBerneCut(int* atid) {
    return OpenMD::SHAPES::Instance()->getGayBerneCut(*atid);
  }

  void fortranDoSHAPESPair(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::SHAPES::Instance()->do_gb_pair(atid1, atid2, d, r, r2, sw,
                                              vdwMult, vpair, pot, A1, A2, f1,
                                              t1, t2);
  }
}
Revision:	1505
Committed:	Sun Oct 3 22:18:59 2010 UTC (14 years, 7 months ago) by gezelter
File size:	14019 byte(s)
Log Message:	Less busted than it was on last check-in, but still won't completely build.
#	User	Rev	Content
1	gezelter	1501	/*
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/SHAPES.hpp"
47			#include "nonbonded/LJ.hpp"
48			#include "utils/simError.h"
49
50			using namespace std;
51			namespace OpenMD {
52	gezelter	1505
53	gezelter	1502	SHAPES::SHAPES() {
54			initialized_ = false;
55			lMax_ = 64;
56			mMax_ = 64;
57			forceField_ = NULL;
58	gezelter	1501	}
59	gezelter	1502
60	gezelter	1501	void SHAPES::initialize() {
61
62			ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
63			ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
64			ForceField::AtomTypeContainer::MapTypeIterator i;
65			AtomType* at;
66
67			// SHAPES handles all of the SHAPES-SHAPES interactions as well as
68			// SHAPES-LJ cross interactions:
69	gezelter	1502
70	gezelter	1501	for (at = atomTypes->beginType(i); at != NULL;
71			at = atomTypes->nextType(i)) {
72
73	gezelter	1505	if (at->isShape())
74			addShape(dynamic_cast<ShapeAtomType*>(at));
75
76			if (at->isLennardJones())
77			addLJ(at);
78
79	gezelter	1501	}
80	gezelter	1502
81	gezelter	1501	initialized_ = true;
82			}
83	gezelter	1502
84	gezelter	1505	void SHAPES::addShape(ShapeAtomType* atomType){
85	gezelter	1501	// add it to the map:
86			AtomTypeProperties atp = atomType->getATP();
87	gezelter	1505
88			pair<map<int,ShapeAtomType*>::iterator, bool> ret;
89			ret = shapesMap.insert( pair<int, ShapeAtomType*>(atp.ident, atomType));
90	gezelter	1501	if (ret.second == false) {
91	gezelter	1505	sprintf( painCave.errmsg,
92	gezelter	1501	"SHAPES already had a previous entry with ident %d\n",
93			atp.ident);
94			painCave.severity = OPENMD_INFO;
95			painCave.isFatal = 0;
96			simError();
97			}
98
99	gezelter	1505	ShapesMap.insert( pair<int, ShapeAtomType*>(atp.ident, sAtomType) );
100
101			} else if (atomType->isLennardJones()) {
102	gezelter	1502	d1 = getLJSigma(atomType) / sqrt(2.0);
103			e1 = getLJEpsilon(atomType);
104	gezelter	1501	} else {
105			sprintf( painCave.errMsg,
106			"SHAPES::addType was passed an atomType (%s) that does not\n"
107	gezelter	1502	"\tappear to be a SHAPES or Lennard-Jones atom.\n",
108	gezelter	1501	atomType->getName().c_str());
109			painCave.severity = OPENMD_ERROR;
110			painCave.isFatal = 1;
111			simError();
112	gezelter	1505	}
113	gezelter	1501	}
114
115
116	gezelter	1502	LJParam SHAPES::getLJParam(AtomType* atomType) {
117
118			// Do sanity checking on the AtomType we were passed before
119			// building any data structures:
120			if (!atomType->isLennardJones()) {
121			sprintf( painCave.errMsg,
122			"SHAPES::getLJParam was passed an atomType (%s) that does not\n"
123			"\tappear to be a Lennard-Jones atom.\n",
124			atomType->getName().c_str());
125			painCave.severity = OPENMD_ERROR;
126			painCave.isFatal = 1;
127			simError();
128			}
129
130			GenericData* data = atomType->getPropertyByName("LennardJones");
131			if (data == NULL) {
132			sprintf( painCave.errMsg, "SHAPES::getLJParam could not find Lennard-Jones\n"
133			"\tparameters for atomType %s.\n", atomType->getName().c_str());
134			painCave.severity = OPENMD_ERROR;
135			painCave.isFatal = 1;
136			simError();
137			}
138
139			LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
140			if (ljData == NULL) {
141			sprintf( painCave.errMsg,
142			"SHAPES::getLJParam could not convert GenericData to LJParam for\n"
143			"\tatom type %s\n", atomType->getName().c_str());
144			painCave.severity = OPENMD_ERROR;
145			painCave.isFatal = 1;
146			simError();
147			}
148
149			return ljData->getData();
150			}
151
152			RealType SHAPES::getLJEpsilon(AtomType* atomType) {
153			LJParam ljParam = getLJParam(atomType);
154			return ljParam.epsilon;
155			}
156			RealType SHAPES::getLJSigma(AtomType* atomType) {
157			LJParam ljParam = getLJParam(atomType);
158			return ljParam.sigma;
159			}
160
161	gezelter	1501	RealType SHAPES::getGayBerneCut(int atid) {
162			if (!initialized_) initialize();
163			std::map<int, AtomType*> :: const_iterator it;
164			it = SHAPESMap.find(atid);
165			if (it == SHAPESMap.end()) {
166			sprintf( painCave.errMsg,
167			"SHAPES::getGayBerneCut could not find atid %d in SHAPESMap\n",
168			(atid));
169			painCave.severity = OPENMD_ERROR;
170			painCave.isFatal = 1;
171			simError();
172			}
173
174			AtomType* atype = it->second;
175
176			RealType gbCut;
177
178			if (atype->isGayBerne()) {
179			GayBerneParam gb = getGayBerneParam(atype);
180
181			// sigma is actually sqrt(2) * l for prolate ellipsoids
182			gbCut = 2.5 * sqrt(2.0) * max(gb.SHAPES_l, gb.SHAPES_d);
183
184			} else if (atype->isLennardJones()) {
185			gbCut = 2.5 * LJ::Instance()->getSigma(atype);
186			}
187
188			return gbCut;
189			}
190
191
192			void SHAPES::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
193			RealType r, RealType r2, RealType sw,
194			RealType &vpair, RealType &pot,
195			RotMat3x3d A1, RotMat3x3d A2, Vector3d &f1,
196			Vector3d &t1, Vector3d &t2) {
197
198			if (!initialized_) initialize();
199
200			pair<AtomType, AtomType> key = make_pair(at1, at2);
201			SHAPESInteractionData mixer = MixingMap[key];
202
203			RealType r3 = r2 * r;
204			RealType r5 = r3 * r2;
205
206			Vector3d drdi = -d / r;
207			Vector3d drdui = V3Zero;
208			Vector3d drdj = d / r;
209			Vector3d drduj = V3Zero;
210
211			bool i_is_LJ = at1->isLennardJones();
212			bool j_is_LJ = at2->isLennardJones();
213
214			RealType sigma_i;
215			RealType s_i;
216			RealType eps_i;
217			Vector3d dsigmaidr;
218			Vector3d disgmaidu;
219			Vector3d dsidr;
220			Vector3d dsidu;
221			Vector3d depsidr;
222			Vector3d depsidu;
223
224			if (i_is_LJ) {
225			sigma_i = LJ::Instance()->getSigma(at1);
226			s_i = sigma_i;
227			epsilon_i = LJ::Instance()->getEpsilon(at1);
228			dsigmaidr = V3Zero;
229			dsigmaidu = V3Zero;
230			dsidr = V3Zero;
231			dsidu = V3Zero;
232			depsidr = V3Zero;
233			depsidu = V3Zero;
234			} else {
235
236			// rotate the inter-particle separation into the two different
237			// body-fixed coordinate systems:
238
239			Vector3d ri = A1 * d;
240
241			RealType xi = ri.x() / r;
242			RealType yi = ri.y() / r;
243			RealType zi = ri.z() / r;
244			RealType xi2 = xi * xi;
245			RealType yi2 = yi * yi;
246			RealType zi2 = zi * zi;
247			RealType cti = zi / r;
248
249			if (cti > 1.0) cti = 1.0;
250			if (cti < -1.0_dp) cti = -1.0;
251
252			Vector3d dctidr(-zi * xi / r3,
253			-zi * yi / r3,
254			1.0 / r - zi2 / r3);
255
256			Vector3d dctidu(yi / r,
257			-zi / r,
258			0.0);
259
260			// this is an attempt to try to truncate the singularity when
261			// sin(theta) is near 0.0:
262
263			RealType sti2 = 1.0 - cti*cti;
264			if (fabs(sti2) < 1.0e-12) {
265			RealType proji = sqrt(r * 1.0e-12);
266			Vector3d dcpidx(1.0 / proji,
267			0.0,
268
269			dcpidx = 1.0_dp / proji
270			dcpidy = 0.0_dp
271			dcpidux = xi / proji
272			dcpiduy = 0.0_dp
273			dspidx = 0.0_dp
274			dspidy = 1.0_dp / proji
275			dspidux = 0.0_dp
276			dspiduy = yi / proji
277			else
278			proji = sqrt(xi2 + yi2)
279			proji3 = projiprojiproji
280			dcpidx = 1.0_dp / proji - xi2 / proji3
281			dcpidy = - xi * yi / proji3
282			dcpidux = xi / proji - (xi2 * xi) / proji3
283			dcpiduy = - (xi * yi2) / proji3
284			dspidx = - xi * yi / proji3
285			dspidy = 1.0_dp / proji - yi2 / proji3
286			dspidux = - (yi * xi2) / proji3
287			dspiduy = yi / proji - (yi2 * yi) / proji3
288			endif
289
290			cpi = xi / proji
291			dcpidz = 0.0_dp
292			dcpiduz = 0.0_dp
293
294			spi = yi / proji
295			dspidz = 0.0_dp
296			dspiduz = 0.0_dp
297
298
299
300
301			RealType sigma0 = mixer.sigma0;
302			RealType dw = mixer.dw;
303			RealType eps0 = mixer.eps0;
304			RealType x2 = mixer.x2;
305			RealType xa2 = mixer.xa2;
306			RealType xai2 = mixer.xai2;
307			RealType xp2 = mixer.xp2;
308			RealType xpap2 = mixer.xpap2;
309			RealType xpapi2 = mixer.xpapi2;
310
311			Vector3d ul1 = A1.getRow(2);
312			Vector3d ul2 = A2.getRow(2);
313
314			RealType a, b, g;
315
316
317			if (i_is_LJ) {
318			a = 0.0;
319			ul1 = V3Zero;
320			} else {
321			a = dot(d, ul1);
322			}
323
324			if (j_is_LJ) {
325			b = 0.0;
326			ul2 = V3Zero;
327			} else {
328			b = dot(d, ul2);
329			}
330
331			if (i_is_LJ \|\| j_is_LJ)
332			g = 0.0;
333			else
334			g = dot(ul1, ul2);
335
336			RealType au = a / r;
337			RealType bu = b / r;
338
339			RealType au2 = au * au;
340			RealType bu2 = bu * bu;
341			RealType g2 = g * g;
342
343			RealType H = (xa2 * au2 + xai2 * bu2 - 2.0x2aubug) / (1.0 - x2*g2);
344			RealType Hp = (xpap2au2 + xpapi2bu2 - 2.0xp2aubug) / (1.0 - xp2*g2);
345
346			RealType sigma = sigma0 / sqrt(1.0 - H);
347			RealType e1 = 1.0 / sqrt(1.0 - x2*g2);
348			RealType e2 = 1.0 - Hp;
349			RealType eps = eps0 * pow(e1,nu_) * pow(e2,mu_);
350			RealType BigR = dwsigma0 / (r - sigma + dwsigma0);
351
352			RealType R3 = BigRBigRBigR;
353			RealType R6 = R3*R3;
354			RealType R7 = R6 * BigR;
355			RealType R12 = R6*R6;
356			RealType R13 = R6*R7;
357
358			RealType U = vdwMult * 4.0 * eps * (R12 - R6);
359
360			RealType s3 = sigmasigmasigma;
361			RealType s03 = sigma0sigma0sigma0;
362
363			RealType pref1 = - vdwMult * 8.0 * eps * mu_ * (R12 - R6) / (e2 * r);
364
365			RealType pref2 = vdwMult * 8.0 * eps * s3 * (6.0R13 - 3.0R7) /(dwrs03);
366
367			RealType dUdr = - (pref1 * Hp + pref2 * (sigma0sigma0r/s3 + H));
368
369			RealType dUda = pref1 * (xpap2au - xp2bug) / (1.0 - xp2 g2)
370			+ pref2 * (xa2 * au - x2 bug) / (1.0 - x2 * g2);
371
372			RealType dUdb = pref1 * (xpapi2bu - xp2aug) / (1.0 - xp2 g2)
373			+ pref2 * (xai2 * bu - x2 aug) / (1.0 - x2 * g2);
374
375			RealType dUdg = 4.0 * eps * nu_ * (R12 - R6) * x2 * g / (1.0 - x2*g2)
376			+ 8.0 * eps * mu_ * (R12 - R6) * (xp2aubu - Hpxp2g) /
377			(1.0 - xp2 * g2) / e2 + 8.0 * eps * s3 * (3.0 * R7 - 6.0 * R13) *
378			(x2 * au * bu - H * x2 * g) / (1.0 - x2 * g2) / (dw * s03);
379
380
381			Vector3d rhat = d / r;
382			Vector3d rxu1 = cross(d, ul1);
383			Vector3d rxu2 = cross(d, ul2);
384			Vector3d uxu = cross(ul1, ul2);
385
386			pot += U*sw;
387			f1 += dUdr * rhat + dUda * ul1 + dUdb * ul2;
388			t1 += dUda * rxu1 - dUdg * uxu;
389			t2 += dUdb * rxu2 - dUdg * uxu;
390			vpair += U*sw;
391
392			return;
393
394			}
395
396			void SHAPES::do_gb_pair(int atid1, int atid2, RealType d, RealType r,
397			RealType r2, RealType sw, RealType *vdwMult,
398			RealType vpair, RealType pot, RealType *A1,
399			RealType A2, RealType f1, RealType t1, RealType t2) {
400
401			if (!initialized_) initialize();
402
403			AtomType* atype1 = SHAPESMap[*atid1];
404			AtomType* atype2 = SHAPESMap[*atid2];
405
406			Vector3d disp(d);
407			Vector3d frc(f1);
408			Vector3d trq1(t1);
409			Vector3d trq2(t2);
410			RotMat3x3d Ai(A1);
411			RotMat3x3d Aj(A2);
412
413			// Fortran has the opposite matrix ordering from c++, so we'll use
414			// transpose here. When we finish the conversion to C++, this wrapper
415			// will disappear, as will the transpose below:
416
417			calcForce(atype1, atype2, disp, r, r2, sw, vdwMult, vpair, pot,
418			Ai, Aj, frc, trq1, trq1);
419
420			f1[0] = frc.x();
421			f1[1] = frc.y();
422			f1[2] = frc.z();
423
424			t1[0] = trq1.x();
425			t1[1] = trq1.y();
426			t1[2] = trq1.z();
427
428			t2[0] = trq2.x();
429			t2[1] = trq2.y();
430			t2[2] = trq2.z();
431
432			return;
433			}
434			}
435
436			extern "C" {
437
438			#define fortranGetGayBerneCut FC_FUNC(getgaybernecut, GETGAYBERNECUT)
439			#define fortranDoSHAPESPair FC_FUNC(do_gb_pair, DO_SHAPES_PAIR)
440
441			RealType fortranGetGayBerneCut(int* atid) {
442			return OpenMD::SHAPES::Instance()->getGayBerneCut(*atid);
443			}
444
445			void fortranDoSHAPESPair(int atid1, int atid2, RealType d, RealType r,
446			RealType r2, RealType sw, RealType *vdwMult,
447			RealType vpair, RealType pot, RealType *A1,
448			RealType A2, RealType f1, RealType t1, RealType t2){
449
450			return OpenMD::SHAPES::Instance()->do_gb_pair(atid1, atid2, d, r, r2, sw,
451			vdwMult, vpair, pot, A1, A2, f1,
452			t1, t2);
453			}
454			}