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::initForce() {
        if (!initialized_) initialize();
  }
 
  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:	1614
Committed:	Tue Aug 23 20:55:51 2011 UTC (14 years, 2 months ago) by mciznick
File size:	14088 byte(s)
Log Message:	Updated scalability of OpenMP threads.
#	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	mciznick	1614	void SHAPES::initForce() {
192			if (!initialized_) initialize();
193			}
194	gezelter	1501
195			void SHAPES::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
196			RealType r, RealType r2, RealType sw,
197			RealType &vpair, RealType &pot,
198			RotMat3x3d A1, RotMat3x3d A2, Vector3d &f1,
199			Vector3d &t1, Vector3d &t2) {
200
201			if (!initialized_) initialize();
202
203			pair<AtomType, AtomType> key = make_pair(at1, at2);
204			SHAPESInteractionData mixer = MixingMap[key];
205
206			RealType r3 = r2 * r;
207			RealType r5 = r3 * r2;
208
209			Vector3d drdi = -d / r;
210			Vector3d drdui = V3Zero;
211			Vector3d drdj = d / r;
212			Vector3d drduj = V3Zero;
213
214			bool i_is_LJ = at1->isLennardJones();
215			bool j_is_LJ = at2->isLennardJones();
216
217			RealType sigma_i;
218			RealType s_i;
219			RealType eps_i;
220			Vector3d dsigmaidr;
221			Vector3d disgmaidu;
222			Vector3d dsidr;
223			Vector3d dsidu;
224			Vector3d depsidr;
225			Vector3d depsidu;
226
227			if (i_is_LJ) {
228			sigma_i = LJ::Instance()->getSigma(at1);
229			s_i = sigma_i;
230			epsilon_i = LJ::Instance()->getEpsilon(at1);
231			dsigmaidr = V3Zero;
232			dsigmaidu = V3Zero;
233			dsidr = V3Zero;
234			dsidu = V3Zero;
235			depsidr = V3Zero;
236			depsidu = V3Zero;
237			} else {
238
239			// rotate the inter-particle separation into the two different
240			// body-fixed coordinate systems:
241
242			Vector3d ri = A1 * d;
243
244			RealType xi = ri.x() / r;
245			RealType yi = ri.y() / r;
246			RealType zi = ri.z() / r;
247			RealType xi2 = xi * xi;
248			RealType yi2 = yi * yi;
249			RealType zi2 = zi * zi;
250			RealType cti = zi / r;
251
252			if (cti > 1.0) cti = 1.0;
253			if (cti < -1.0_dp) cti = -1.0;
254
255			Vector3d dctidr(-zi * xi / r3,
256			-zi * yi / r3,
257			1.0 / r - zi2 / r3);
258
259			Vector3d dctidu(yi / r,
260			-zi / r,
261			0.0);
262
263			// this is an attempt to try to truncate the singularity when
264			// sin(theta) is near 0.0:
265
266			RealType sti2 = 1.0 - cti*cti;
267			if (fabs(sti2) < 1.0e-12) {
268			RealType proji = sqrt(r * 1.0e-12);
269			Vector3d dcpidx(1.0 / proji,
270			0.0,
271
272			dcpidx = 1.0_dp / proji
273			dcpidy = 0.0_dp
274			dcpidux = xi / proji
275			dcpiduy = 0.0_dp
276			dspidx = 0.0_dp
277			dspidy = 1.0_dp / proji
278			dspidux = 0.0_dp
279			dspiduy = yi / proji
280			else
281			proji = sqrt(xi2 + yi2)
282			proji3 = projiprojiproji
283			dcpidx = 1.0_dp / proji - xi2 / proji3
284			dcpidy = - xi * yi / proji3
285			dcpidux = xi / proji - (xi2 * xi) / proji3
286			dcpiduy = - (xi * yi2) / proji3
287			dspidx = - xi * yi / proji3
288			dspidy = 1.0_dp / proji - yi2 / proji3
289			dspidux = - (yi * xi2) / proji3
290			dspiduy = yi / proji - (yi2 * yi) / proji3
291			endif
292
293			cpi = xi / proji
294			dcpidz = 0.0_dp
295			dcpiduz = 0.0_dp
296
297			spi = yi / proji
298			dspidz = 0.0_dp
299			dspiduz = 0.0_dp
300
301
302
303
304			RealType sigma0 = mixer.sigma0;
305			RealType dw = mixer.dw;
306			RealType eps0 = mixer.eps0;
307			RealType x2 = mixer.x2;
308			RealType xa2 = mixer.xa2;
309			RealType xai2 = mixer.xai2;
310			RealType xp2 = mixer.xp2;
311			RealType xpap2 = mixer.xpap2;
312			RealType xpapi2 = mixer.xpapi2;
313
314			Vector3d ul1 = A1.getRow(2);
315			Vector3d ul2 = A2.getRow(2);
316
317			RealType a, b, g;
318
319
320			if (i_is_LJ) {
321			a = 0.0;
322			ul1 = V3Zero;
323			} else {
324			a = dot(d, ul1);
325			}
326
327			if (j_is_LJ) {
328			b = 0.0;
329			ul2 = V3Zero;
330			} else {
331			b = dot(d, ul2);
332			}
333
334			if (i_is_LJ \|\| j_is_LJ)
335			g = 0.0;
336			else
337			g = dot(ul1, ul2);
338
339			RealType au = a / r;
340			RealType bu = b / r;
341
342			RealType au2 = au * au;
343			RealType bu2 = bu * bu;
344			RealType g2 = g * g;
345
346			RealType H = (xa2 * au2 + xai2 * bu2 - 2.0x2aubug) / (1.0 - x2*g2);
347			RealType Hp = (xpap2au2 + xpapi2bu2 - 2.0xp2aubug) / (1.0 - xp2*g2);
348
349			RealType sigma = sigma0 / sqrt(1.0 - H);
350			RealType e1 = 1.0 / sqrt(1.0 - x2*g2);
351			RealType e2 = 1.0 - Hp;
352			RealType eps = eps0 * pow(e1,nu_) * pow(e2,mu_);
353			RealType BigR = dwsigma0 / (r - sigma + dwsigma0);
354
355			RealType R3 = BigRBigRBigR;
356			RealType R6 = R3*R3;
357			RealType R7 = R6 * BigR;
358			RealType R12 = R6*R6;
359			RealType R13 = R6*R7;
360
361			RealType U = vdwMult * 4.0 * eps * (R12 - R6);
362
363			RealType s3 = sigmasigmasigma;
364			RealType s03 = sigma0sigma0sigma0;
365
366			RealType pref1 = - vdwMult * 8.0 * eps * mu_ * (R12 - R6) / (e2 * r);
367
368			RealType pref2 = vdwMult * 8.0 * eps * s3 * (6.0R13 - 3.0R7) /(dwrs03);
369
370			RealType dUdr = - (pref1 * Hp + pref2 * (sigma0sigma0r/s3 + H));
371
372			RealType dUda = pref1 * (xpap2au - xp2bug) / (1.0 - xp2 g2)
373			+ pref2 * (xa2 * au - x2 bug) / (1.0 - x2 * g2);
374
375			RealType dUdb = pref1 * (xpapi2bu - xp2aug) / (1.0 - xp2 g2)
376			+ pref2 * (xai2 * bu - x2 aug) / (1.0 - x2 * g2);
377
378			RealType dUdg = 4.0 * eps * nu_ * (R12 - R6) * x2 * g / (1.0 - x2*g2)
379			+ 8.0 * eps * mu_ * (R12 - R6) * (xp2aubu - Hpxp2g) /
380			(1.0 - xp2 * g2) / e2 + 8.0 * eps * s3 * (3.0 * R7 - 6.0 * R13) *
381			(x2 * au * bu - H * x2 * g) / (1.0 - x2 * g2) / (dw * s03);
382
383
384			Vector3d rhat = d / r;
385			Vector3d rxu1 = cross(d, ul1);
386			Vector3d rxu2 = cross(d, ul2);
387			Vector3d uxu = cross(ul1, ul2);
388
389			pot += U*sw;
390			f1 += dUdr * rhat + dUda * ul1 + dUdb * ul2;
391			t1 += dUda * rxu1 - dUdg * uxu;
392			t2 += dUdb * rxu2 - dUdg * uxu;
393			vpair += U*sw;
394
395			return;
396
397			}
398
399			void SHAPES::do_gb_pair(int atid1, int atid2, RealType d, RealType r,
400			RealType r2, RealType sw, RealType *vdwMult,
401			RealType vpair, RealType pot, RealType *A1,
402			RealType A2, RealType f1, RealType t1, RealType t2) {
403
404			if (!initialized_) initialize();
405
406			AtomType* atype1 = SHAPESMap[*atid1];
407			AtomType* atype2 = SHAPESMap[*atid2];
408
409			Vector3d disp(d);
410			Vector3d frc(f1);
411			Vector3d trq1(t1);
412			Vector3d trq2(t2);
413			RotMat3x3d Ai(A1);
414			RotMat3x3d Aj(A2);
415
416			// Fortran has the opposite matrix ordering from c++, so we'll use
417			// transpose here. When we finish the conversion to C++, this wrapper
418			// will disappear, as will the transpose below:
419
420			calcForce(atype1, atype2, disp, r, r2, sw, vdwMult, vpair, pot,
421			Ai, Aj, frc, trq1, trq1);
422
423			f1[0] = frc.x();
424			f1[1] = frc.y();
425			f1[2] = frc.z();
426
427			t1[0] = trq1.x();
428			t1[1] = trq1.y();
429			t1[2] = trq1.z();
430
431			t2[0] = trq2.x();
432			t2[1] = trq2.y();
433			t2[2] = trq2.z();
434
435			return;
436			}
437			}
438
439			extern "C" {
440
441			#define fortranGetGayBerneCut FC_FUNC(getgaybernecut, GETGAYBERNECUT)
442			#define fortranDoSHAPESPair FC_FUNC(do_gb_pair, DO_SHAPES_PAIR)
443
444			RealType fortranGetGayBerneCut(int* atid) {
445			return OpenMD::SHAPES::Instance()->getGayBerneCut(*atid);
446			}
447
448			void fortranDoSHAPESPair(int atid1, int atid2, RealType d, RealType r,
449			RealType r2, RealType sw, RealType *vdwMult,
450			RealType vpair, RealType pot, RealType *A1,
451			RealType A2, RealType f1, RealType t1, RealType t2){
452
453			return OpenMD::SHAPES::Instance()->do_gb_pair(atid1, atid2, d, r, r2, sw,
454			vdwMult, vpair, pot, A1, A2, f1,
455			t1, t2);
456			}
457			}