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() || at->isLennardJones())
        addType(at);
    }
    
    initialized_ = true;
  }
  
  void SHAPES::addType(AtomType* atomType){
    // add it to the map:
    AtomTypeProperties atp = atomType->getATP();    
    
    pair<map<int,AtomType*>::iterator,bool> ret;    
    ret = ShapesMap.insert( pair<int, AtomType*>(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();         
    }   
    
    if (atomType->isShape()) {
      ShapeAtomType* sAtomType = dynamic_cast<ShapeAtomType*>(atomType);
      if (sAtomType == NULL) {
        sprintf(painCave.errMsg,
                "SHAPES:: Can't cast to ShapeAtomType");
        painCave.severity = OPENMD_ERROR;
        painCave.isFatal = 1;
        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();
    }
    

    // Now, iterate over all known types and add to the mixing map:
    
    map<int, AtomType*>::iterator it;
    for( it = ShapesMap.begin(); it != SHAPESMap.end(); ++it) {
      
      AtomType* atype2 = (*it).second;
      
      RealType d2, l2, e2, er2, dw2;
      
      if (atype2->isGayBerne()) {
        GayBerneParam gb2 = getGayBerneParam(atype2);
        d2 = gb2.SHAPES_d;
        l2 = gb2.SHAPES_l;
        e2 = gb2.SHAPES_eps;
        er2 = gb2.SHAPES_eps_ratio;
        dw2 = gb2.SHAPES_dw;
      } else if (atype2->isLennardJones()) {
        d2 = getLJSigma(atype2) / sqrt(2.0);
        e2 = getLJEpsilon(atype2);
        l2 = d2;
        er2 = 1.0;
        dw2 = 1.0;
      } 
                       
      SHAPESInteractionData 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;
        }
      }
    }      
  }
  

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